Documentation of gnns

Global Index (short | long) | Local contents | Local Index (short | long)

Function Synopsis

[g,R,grad]=gnns(z,th,sspm,T,modarg,sqrlam,lim,maxsize,index)

Help text

GNNS    Computes the Gauss-Newton direction for PE-criteria.

    [GN,R,G] = gnns(Z,TH,SSPM,MODARG,SQRLAM,LIM,MAXSIZE,INDEX)

    GN : The Gauss-Newton direction
    R : The Hessian of the criterion (using the GN-approximation)
    G: The gradient direction
    The routine is intended primarily as a subroutine to 
    PEMSS
    See PEMSS for an explanation of arguments.

Cross-Reference Information

Listing of function gnns

function [g,R,grad]=gnns(z,th,sspm,T,modarg,sqrlam,lim,maxsize,index)

%    L. Ljung 10-1-89, 9-25-93
%   Copyright (c) 1986-98 by The MathWorks, Inc.
%   $Revision: 2.5 $  $Date: 1998/05/22 01:18:23 $

% *** Set up the model orders ***


nd=length(th);
n=length(index);
[Ncap,nz]=size(z);

% *** Prepare for gradients calculations ***
[A,B,C,D,K,X0]=feval(sspm,th,T,modarg);
[ny,nx]=size(C);

% *** Compute the gradient PSI. If N>M do it in portions ***
rowmax=max(n*ny,nx+nz);
M=floor(maxsize/rowmax);
R=zeros(n);Fcap=zeros(n,1);dt=nuderst(th);
for kc=1:M:Ncap
      jj=(kc:min(Ncap,kc-1+M));
      if jj(length(jj))<Ncap,jjz=[jj,jj(length(jj))+1];else jjz=jj;end
      psitemp=zeros(length(jj),ny);
      psi=zeros(ny*length(jj),n);
      x=ltitr(A-K*C,[K B-K*D],z(jjz,:),X0);
      yh=(C*x(1:length(jj),:).'+[zeros(ny,ny) D]*z(jj,:).').';
      e=z(jj,1:ny)-yh;
      [nxr,nxc]=size(x);X0=x(nxr,:).';
      if lim==0
         el=e*sqrlam;
      else 
         ll=ones(length(jj),1)*lim; 
         la=abs(e)+eps*ll;regul=sqrt(min(la,ll)./la);el=e.*regul*sqrlam;
      end

      evec=el(:);
      kkl=1;
        for kl=index
drawnow
             th1=th;th1(kl)=th1(kl)+dt(kl)/2;
             th2=th;th2(kl)=th2(kl)-dt(kl)/2;
             [A1,B1,C1,D1,K1,X1]=feval(sspm,th1,T,modarg);
             [A2,B2,C2,D2,K2,X2]=feval(sspm,th2,T,modarg);
             dA=(A1-A2)/dt(kl);dB=(B1-B2)/dt(kl);dC=(C1-C2)/dt(kl);
             dD=(D1-D2)/dt(kl); dK=(K1-K2)/dt(kl);
             if kc==1,dX=(X1-X2)/dt(kl);else dX=dXk(:,kl);end
        psix=ltitr(A-K*C,[dA-dK*C-K*dC,dK,dB-K*dD-dK*D],[x,z(jjz,:)],dX);
        [rpsix,cpsix]=size(psix);
        dXk(:,kl)=psix(rpsix,:)';
        psitemp=(C*psix(1:length(jj),:).' + ...
          [dC,zeros(ny,ny),dD]*[x(1:length(jj),:),z(jj,:)].').'*sqrlam;
        if ~(lim==0),psitemp=psitemp.*regul;end
     psi(:,kkl)=psitemp(:);kkl=kkl+1;   
   end
  
    R=R+psi'*psi;  Fcap=Fcap+ psi'*evec;
end

% *** Compute the Gauss-Newton direction g ***

if Ncap>M, g=R\Fcap; else g=psi\evec;end
g1=zeros(nd,1);g1(index)=g;g=g1;
grad=zeros(nd,1);grad(index)=Fcap;