home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / ETA2SS.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.6 KB  |  89 lines
  1. function [a,b,c,d,k,x0,da,db,dc,dd,dk,dx0]=eta2ss(eta)
  2. %ETA2SS  Auxiliary routine to TH2SS
  3. %
  4. %    [A,B,C,D,K,X0]=eta2ss(TH)
  5. %
  6. %    A,B etc: The state space matrices corresponding to TH. 
  7. %
  8. %    with
  9. %    [A,B,C,D,K,X0,dA,dB,dC,dD,dK,dX0]=eta2ss(TH)
  10. %    
  11. %    also the standard deviations dA, dB etc are obtained
  12.  
  13. %    L. Ljung 7-10-88
  14. %    Copyright (c) 1988-90 The MathWorks, Inc.
  15. %    All Rights Reserved
  16.  
  17. [nr,nc]=size(eta);
  18. T=eta(1,2);
  19.  
  20. sspm=getmfth(eta);
  21.  
  22. nd=eta(1,5);
  23. carg=eta(1,7);
  24. etapar=eta(3,1:nd);if nd==0,etapar=0;end
  25. arg=getargth(eta);
  26.  
  27. if any(eta(2,8)==[2 3 4]) & T<0,flag=1;else flag=0;end
  28. if eta(2,8)==3 & T<0,flaga=1;else flaga=0;end
  29. if any(eta(2,8)==[2 3]),Tmod=-1;else Tmod=T;end
  30. [a,b,c,d,k,x0]=feval(sspm,etapar,Tmod,arg);
  31. [nx,nu]=size(b);[ny,nx]=size(c);
  32. if flag,
  33.     [a,bk,c,d]=contth(a,[b k],c,d,abs(T),flaga);
  34.     b=bk(:,1:eta(1,3));k=bk(:,eta(1,3)+1:eta(1,3)+eta(1,4));
  35. end
  36. input=1;if isempty(b),input=0;end
  37. if nargout>6
  38.   if nd==0,da=zeros(nx,nx);db=zeros(nx,nu);dc=zeros(ny,nx);dd=zeros(ny,nu);
  39.            dk=zeros(nx,ny);dx0=zeros(nx,1);
  40.   else   
  41.   if (eta(2,8)==1 & T<0) | (any(eta(2,8)==[2 3]) & T>0)
  42.       p=diag(eta(4:3+nd,1:nd));p=sqrt(p')+eps;
  43.       etap1=etapar+p;
  44.       [a1,b1,c1,d1,k1,x01]=feval(sspm,etap1,Tmod,arg);
  45.       if flag, 
  46.          [a1,bk1,c1,d1]=contth(a1,[b1 k1],c1,d1,abs(T),flaga);
  47.      b1=bk1(:,1:eta(1,3));k1=bk1(:,eta(1,3)+1:eta(1,3)+eta(1,4));
  48.       end
  49.       da=abs(a1-a);
  50.       db=abs(b1-b);
  51.       dc=abs(c1-c);
  52.       dd=abs(d1-d);
  53.       dk=abs(k1-k);
  54.       dx0=abs(x01-x0);
  55.   else
  56.       [etapar,p]=th2par(eta);
  57.       if isempty(p) | norm(p)==0,
  58.             [nx,nu]=size(b);[ny,nx]=size(c);    
  59.           da=zeros(nx,nx);db=zeros(nx,nu);dc=zeros(ny,nx);
  60.       dd=zeros(ny,nu);dk=zeros(nx,ny);dx0=zeros(nx,1);
  61.       return,end
  62.       dt=nuderst(etapar); 
  63.       for kk=1:length(dt)
  64.         etap1=etapar;etap1(kk)=etap1(kk)+dt(kk);
  65.     [a1,b1,c1,d1,k1,x01]=feval(sspm,etap1,Tmod,arg);
  66.     if flag, 
  67.       [a1,bk1,c1,d1]=contth(a1,[b1 k1],c1,d1,abs(T),flaga);
  68.       b1=bk1(:,1:eta(1,3));k1=bk1(:,eta(1,3)+1:eta(1,3)+eta(1,4));
  69.     end
  70.     da=a1-a;cola(:,kk)=da(:)/dt(kk);
  71.     dc=c1-c;colc(:,kk)=dc(:)/dt(kk);
  72.     dk=k1-k;colk(:,kk)=dk(:)/dt(kk);dx0=x01-x0;colx0(:,kk)=dx0(:)/dt(kk);
  73.     if input
  74.        db=b1-b;colb(:,kk)=db(:)/dt(kk);dd=d1-d;cold(:,kk)=dd(:)/dt(kk);
  75.         end
  76.       end
  77.       pa=cola*p*cola';da(:)=sqrt(diag(pa));
  78.       pc=colc*p*colc';dc(:)=sqrt(diag(pc));
  79.       pk=colk*p*colk';dk(:)=sqrt(diag(pk));
  80.       px0=colx0*p*colx0';dx0(:)=sqrt(diag(px0));
  81.       if input
  82.         pb=colb*p*colb';db(:)=sqrt(diag(pb));
  83.         pd=cold*p*cold';dd(:)=sqrt(diag(pd));
  84.       else db=[];dd=[];
  85.       end
  86.    end
  87.    end
  88. end
  89.