home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / OE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.7 KB  |  146 lines
  1. function th=oe(z,nn,maxiter,tol,lim,maxsize,Tsamp)
  2. %OE    Computes the prediction error estimate of an output-error model.
  3. %
  4. %    TH=oe(Z,NN)
  5. %
  6. %    TH: Returned as the estimated parameters of the output-error model
  7. %    y(t) = [B(q)/F(q)] u(t-nk) + e(t)
  8. %    along with estimated covariances and structure information.
  9. %    For the exact format of TH, see HELP THETA.
  10. %
  11. %    Z : The output-input data Z=[y u], with y and u being column vectors.
  12. %    This routine only works for single-input systems. Use PEM otherwise.
  13. %
  14. %    NN: Initial value and structure information, given either as
  15. %    NN=[nb nf nk], the orders and delay of the above model, or as
  16. %    NN=THI, with THI being a theta matrix of the standrad format. In
  17. %    the latter case the criterion minimization is initialized at THI.
  18. %
  19. %    Some parameters associated with the algorithm are accessed by
  20. %    TH = oe(Z,NN,maxiter,tol,lim,maxsize,T)
  21. %    See HELP AUXVAR for an explanation of these, and their default values.
  22.  
  23. %    L. Ljung 10-1-86,4-7-92
  24. %    Copyright (c) 1986-92 by the MathWorks, Inc.
  25. %    All Rights Reserved.
  26.  
  27. % *** Set up default values ***
  28. [nr,cc]=size(nn); [Ncap,nz]=size(z); nu=nz-1;
  29. maxsdef=idmsize(Ncap);
  30.  
  31. if nargin<7, Tsamp=[];end
  32. if nargin<6, maxsize=[];end
  33. if nargin<5, lim=[];end
  34. if nargin<4, tol=[];, end
  35. if nargin<3, maxiter=[]; end
  36. if isempty(Tsamp),if nr>1,Tsamp=gett(nn);else Tsamp=1;end,end
  37. if isempty(maxsize),maxsize=maxsdef;end
  38. if isempty(lim),lim=1.6;end,if isempty(tol),tol=0.01;end
  39. if isempty(maxiter),maxiter=10;end
  40. if Tsamp<0,Tsamp=1;end, if maxsize<0,maxsize=maxsdef;end, if lim<0,lim=1.6;end
  41. if tol<0,tol=0.01;end,if maxiter<0, maxiter=10;end
  42.  
  43.  
  44. %  *** Do some consistency tests ***
  45.  
  46. if nz>Ncap, error('The data should be organized in column vectors')
  47. return,end
  48. if nu>1, error('This routine only works for single-input systems. Use PEM instead!')
  49. return,end
  50. if nu==0, error('OE makes sense only if an input is present!')
  51. return,end
  52. if nr==1 & cc~=3,
  53.          error('Incorrect number of orders specified: Should be  nn=[nb nf nk]')
  54. return,end
  55. %
  56. %*** if nn has more than 1 row, then it is a theta matrix
  57. %     and we jump directly to the minimization ***
  58. %
  59. if nr>1 nu=nn(1,3);
  60.    if nu>1,
  61.       error('This routine only works for single-input systems')
  62.    return,end
  63.    nf=nn(1,8); nb=nn(1,5); nk=nn(1,9); n=nf+nb;ni=max(nf,nb+nk-1);
  64.    t=nn(3,1:n)';
  65.    if nf>0; f=[1 t(nb+1:n)'];else f=1;end
  66. end
  67. %
  68. %
  69. if nr==1
  70.    nb=nn(1); nf=nn(2); nk=nn(3); n=nb+nf;ni=max(nf,nb+nk-1);
  71. %
  72. %
  73. % *** compute initial estimate via LS and IV ***
  74. %
  75.   t=arx(z,[nf nb nk],maxsize,Tsamp,0);
  76.   f=fstab([1 t(1:nf)']);
  77.         
  78.   t=iv(z,[nf nb nk],f,[zeros(1,nk) t(nf+1:n)'],maxsize,Tsamp,0);
  79.   f=fstab([1 t(1:nf)']);
  80.         
  81. %
  82.   if nf>0, t=[t(nf+1:n); t(1:nf)];end
  83. end
  84. % *** Display initial estimate ***
  85. %
  86. w=pefilt([zeros(1,nk) t(1:nb)'],f,z(:,2),z(1:ni,1));
  87. v=z(:,1)-w; Vcap=v'*v/(Ncap-ni);
  88. if isnan(Vcap)|~finite(Vcap)
  89.     error('Your data or model structure has some deficiencies. Try a different input or a lower order model!'),end
  90. clc, disp(['  INITIAL ESTIMATE'])
  91. disp(['Current loss:' num2str(Vcap)])
  92. disp(['th-vector'])
  93. disp(t)
  94. %
  95. % *** start minimizing ***
  96. %
  97. % ** determine limit for robustification **
  98.  if lim~=0, lim=median(abs(v-median(v)))*lim/0.7;end
  99.  if lim==0, lim=1000*Vcap;end
  100.  [mmv,nnv]=size(v);
  101.  ll=ones(mmv,nnv)*lim;vl=max(min(v,ll),-ll); clear ll
  102.  g=ones(n,1); l=0; st=0; nmax=max(nf,nb+nk-1);
  103. %
  104. % ** the minimization loop **
  105. %
  106. while [norm(g)>tol l<maxiter st~=1]
  107.       l=l+1;
  108. %     * compute gradient *
  109.       wf=filter(-1,f,w); uf=filter(1,f,z(:,2));
  110.    M=floor(maxsize/n);
  111.    R=zeros(n);Fcap=zeros(n,1);
  112.    for k=nmax:M:Ncap-1
  113.       jj=(k+1:min(Ncap,k+M));
  114.       psi=zeros(length(jj),n);
  115.       for k=1:nb, psi(:,k)=uf(jj-k-nk+1);end
  116.       for k=1:nf, psi(:,k+nb)=wf(jj-k);end
  117.       if Ncap>M,R=R+psi'*psi; Fcap=Fcap+psi'*vl(jj);end
  118.    end
  119.    if Ncap>M, g=R\Fcap; else g=psi\vl(jj);end
  120. %
  121. %     * search along the g-direction *
  122. %
  123.       [t1,w,vl,v,V1,f,st]=searchoe(z,t,g,lim,Vcap,nb,nf,nk,ni);
  124.       home
  125.       disp(['      ITERATION # ' int2str(l)])
  126.       disp(['Current loss: ' num2str(V1) '  Previous loss: ' num2str(Vcap)])
  127.       disp(['Current th  prev. th   gn-dir'])
  128.       disp([t1 t g])
  129.       disp(['Norm of gn-vector: ' num2str(norm(g))])
  130.       if st==1
  131.      disp('No improvement of the criterion possible along the search direction')
  132.      disp('Iterations therefore terminated'),end
  133.       t=t1; Vcap=V1;
  134. end
  135. Vcap=v'*v/(length(v)-ni);
  136. th(1,:)=[Vcap Tsamp 1 0 nb 0 0 nf nk];
  137. th(2,1)=Vcap*(1+n/Ncap)/(1-n/Ncap);
  138. ti=fix(clock); ti(1)=ti(1)/100;
  139. th(2,2:6)=ti(1:5);
  140. th(2,7)=8;
  141. th(3,1:n)=t';
  142. if maxiter==0, return,end
  143. if Ncap>M, PP=inv(R); else PP=inv(psi'*psi);end
  144. th(4:3+n,1:n)=Vcap*PP;
  145.  
  146.