home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / GNNS.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.2 KB  |  67 lines
  1. function [g,R,grad]=gnns(z,th,sspm,T,modarg,sqrlam,lim,maxsize,index)
  2. %GNNS        Computes the Gauss-Newton direction for PE-criteria
  3. %
  4. %        [GN,R,G] = gnns(Z,TH,SSPM,MODARG,SQRLAM,LIM,MAXSIZE,INDEX)
  5. %
  6. %        GN : The Gauss-Newton direction
  7. %        R : The Hessian of the criterion (using the GN-approximation)
  8. %     G: The gradient direction
  9. %        The routine is intended primarily as a subroutine to 
  10. %     PEMSS
  11. %        See PEMSS for an explanation of arguments.
  12.  
  13. %        L. Ljung 10-1-89
  14. %        Copyright (c) 1989-90 by the MathWorks, Inc.
  15. %        All Rights Reserved.
  16.  
  17. % *** Set up the model orders ***
  18.  
  19.  
  20. nd=length(th);
  21. n=length(index);
  22. [Ncap,nz]=size(z);
  23.  
  24. % *** Prepare for gradients calculations ***
  25. [A,B,C,D,K,X0]=feval(sspm,th,T,modarg);
  26. [ny,nx]=size(C);
  27.  
  28. % *** Compute the gradient PSI. If N>M do it in portions ***
  29. rowmax=max(n*ny,nx+nz);
  30. M=floor(maxsize/rowmax);
  31. R=zeros(n);Fcap=zeros(n,1);dt=nuderst(th);
  32. for kc=1:M:Ncap-1
  33.       jj=(kc:min(Ncap,kc-1+M));
  34.       if jj(length(jj))<Ncap,jjz=[jj,jj(length(jj))+1];else jjz=jj;end
  35.       psitemp=zeros(length(jj),ny);
  36.       psi=zeros(ny*length(jj),n);
  37.       x=ltitr(A-K*C,[K B-K*D],z(jjz,:),X0);
  38.       yh=(C*x(1:length(jj),:)'+[zeros(ny,ny) D]*z(jj,:)')';
  39.       e=z(jj,1:ny)-yh;
  40.       X0=x(length(x),:)';
  41.       ll=ones(length(jj),1)*lim; el=min(e,ll); el=max(el,-ll)*sqrlam; 
  42.       evec=el(:);
  43.       kkl=1;
  44.     for kl=index
  45.              th1=th;th1(kl)=th1(kl)+dt(kl)/2;
  46.              th2=th;th2(kl)=th2(kl)-dt(kl)/2;
  47.              [A1,B1,C1,D1,K1,X1]=feval(sspm,th1,T,modarg);
  48.              [A2,B2,C2,D2,K2,X2]=feval(sspm,th2,T,modarg);
  49.              dA=(A1-A2)/dt(kl);dB=(B1-B2)/dt(kl);dC=(C1-C2)/dt(kl);
  50.          dD=(D1-D2)/dt(kl); dK=(K1-K2)/dt(kl);
  51.          if kc==1,dX=(X1-X2)/dt(kl);else dX=dXk(:,kl);end
  52.         psix=ltitr(A-K*C,[dA-dK*C-K*dC,dK,dB-K*dD-dK*D],[x,z(jjz,:)],dX);
  53.     dXk(:,kl)=psix(length(psix),:)';
  54.     psitemp=(C*psix(1:length(jj),:)' + ...
  55.     [dC,zeros(ny,ny),dD]*[x(1:length(jj),:),z(jj,:)]')'*sqrlam;
  56.      psi(:,kkl)=psitemp(:);kkl=kkl+1;    
  57.      end
  58.  
  59.     R=R+psi'*psi;  Fcap=Fcap+ psi'*evec;
  60. end
  61.  
  62. % *** Compute the Gauss-Newton direction g ***
  63.  
  64. if Ncap>M, g=R\Fcap; else g=psi\evec;end
  65. g1=zeros(nd,1);g1(index)=g;g=g1;
  66. grad=zeros(nd,1);grad(index)=Fcap;
  67.