home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / PEMSS.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.9 KB  |  120 lines
  1. function eta=pemss(z,etainit,index,maxiter,tol,lim,maxsize,T)
  2. %PEMSS   Minimizes prediction error criteria
  3. %
  4. %        TH = pemss(Z,TH0)
  5. %
  6. %     TH: The result is returned in this matrix. See HELP THETA
  7. %     Z: The datamatrix Z = [Y U], where each column of Y and U
  8. %        corresponds to a particular output and input signal.
  9. %     TH0: This argument provides the model structure,
  10. %         the initial conditions, and the quadratic norm.
  11. %        (See MS2TH, MF2TH and THINIT)
  12. %
  13. %     With TH = pemss(Z,TH0,INDEX) only the adjustable parameters
  14. %     with indices listed in the row vector INDEX will be estimated
  15. %     Remaining parameters will be fixed to their values in ETA0
  16. %
  17. %     Optional arguments can be reached by
  18. %     TH = pemss(Z,TH0,INDEX,MAXITER,TOL,LIM,MAXSIZE,T)
  19. %     See Help AUXVAR for an explanation of the these arguments.
  20. %     In the m-file NUDERSTEP the stepsize to be used for numerical 
  21. %     differentiation is set.
  22.  
  23. %        L. Ljung 10-2-90
  24. %        Copyright (c) 1986-92 by the MathWorks, Inc.
  25. %        All Rights Reserved.
  26.  
  27. [Ncap,nz]=size(z);
  28. maxdef=idmsize(Ncap);
  29. if nargin<8,T=etainit(1,2);end
  30. if nargin<7,maxsize=maxdef;end
  31. if nargin<6,lim=1.7;end
  32. if nargin<5,tol=0.01;end
  33. if nargin<4,maxiter=10;end
  34. if nargin<3,index=[];end
  35. [nr,nc]=size(etainit);
  36. nd=etainit(1,5);
  37. if nd==0,error('This model structure does not have any adjustable parameters!'),end
  38.  
  39. if isempty(index),index=1:nd;end
  40. sspm=getmfth(etainit);
  41. if any(etainit(2,8)==[2 3 4]), Tmod=-abs(T);else Tmod=abs(T);end
  42. [th,P,lambda]=th2par(etainit);
  43. [modarg,rarg,carg]=getargth(etainit);
  44. n=nd;
  45. [A,B,C,D,K,X0]=feval(sspm,th,Tmod,modarg);[ny,nx]=size(C);
  46. [nx,nu]=size(B);[ny,nx]=size(C);
  47. rowmax=max(length(index)*ny,nx+nz);
  48. Mloop=floor(maxsize/rowmax);
  49. if (nu+ny)~=nz, error('Incorrect number of inputs and outputs specified in data vector'),end
  50. ei=eig(A-K*C);
  51. if max(abs(ei))>1.1,error('Initial guess predictor unstable!'),end
  52. for kc=1:Mloop:Ncap-1
  53.       jj=(kc:min(Ncap,kc-1+Mloop));
  54.       if jj(length(jj))<Ncap,jjz=[jj,jj(length(jj))+1];else jjz=jj;end
  55.       xh=ltitr(A-K*C,[K B-K*D],z(jjz,:),X0);
  56.       yh(jj,:)=(C*xh(1:length(jj),:)'+[zeros(ny,ny) D]*z(jj,:)')';
  57.       X0=xh(length(xh),:)';
  58. end
  59. e=z(:,1:ny)-yh;
  60. V=det(e'*e/length(e));
  61.  
  62. g=ones(n,1); 
  63. %** Determine limit for robustification of criterion ***
  64.  
  65. if lim~=0, lim=median(abs(e-ones(Ncap,1)*median(e)))*lim/0.7;end
  66. if lim==0, lim=1000*ones(1,ny)*V;end
  67.         
  68. el=min(e,ones(Ncap,1)*lim); el=max(el,-ones(Ncap,1)*lim);
  69. V=det(el'*e/length(e));
  70. clear e, l=0; st=0;
  71.  
  72. %  ** The minimization loop **
  73.  
  74. while [norm(g)>tol l<maxiter st~=1]
  75.    l=l+1;
  76.  
  77. %   * Compute the Gauss-Newton direction *
  78.    if norm(lambda-eye(ny))==0,lam=eye(ny);else lam=inv(sqrtm(lambda));end
  79.    [g,R,ng]=gnns(z,th,sspm,Tmod,modarg,lam,lim,maxsize,index);
  80.  
  81. %    * Search along the g-direction for a lower value of V *
  82.  
  83.     [Vn,th1,st,lambda]=sess(V,z,th,sspm,Tmod,modarg,lam,g,lim,Mloop);
  84.  
  85.     if st==1,
  86. [Vn,th1,st,lambda]=...
  87. sess(V,z,th,sspm,Tmod,modarg,lam,ng/trace(R)*length(R),lim,Mloop);
  88.     end
  89.  
  90. %  * Display the current values along with the previous ones *
  91.  
  92.      t=zeros(n,3); t(:,1)=th1';t(:,2)=th';t(:,3)=g;
  93.      home,clc
  94.      disp(['    ITERATION # ' int2str(l)])
  95.      disp(['Current loss: ' num2str(Vn) '   Previous loss: ' num2str(V)])
  96.       disp(['Current th  prev. th   gn-dir. '])
  97.       disp(t)
  98.       disp(['Norm of gn-vector: ' num2str(norm(g))])
  99.       if st==1
  100.           disp('No improvement of the criterion possible along the search direction')
  101.           disp('Iterations therefore terminated')
  102.       end
  103.       th=th1;V=Vn;
  104. end
  105.  
  106. % *** Build up the ETA-matrix ***
  107.  
  108. eta = etainit;
  109. P1=zeros(nd,nd);
  110. P1(index,index)=inv(R);
  111. P=P1;
  112. eta(3,1:nd)=th;
  113. eta(1,1)=det(lambda);eta(2,1)=eta(1,1)*(1+nd/Ncap)/(1-nd/Ncap);
  114. eta(4:nd+3,1:nd)=P;
  115. eta(4+nd+etainit(1,6):3+ny+nd+etainit(1,6),1:ny)=lambda;
  116. ti=fix(clock); ti(1)=ti(1)/100;
  117. eta(2,2:6)=ti(1:5); eta(2,7)=23;
  118. if eta(2,8)==3,eta(2,7)=33;end
  119.  
  120.