home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / IVSTRUC.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.4 KB  |  103 lines
  1. function V=ivstruc(z,zv,nn,p,maxsize)
  2. %IVSTRUC    Computes the output error fit for families of single-output
  3. %    ARX-models.
  4. %
  5. %    V = ivstruc(ZE,ZV,NN)
  6. %
  7. %    V: The first row of V is returned as the output error fit for 
  8. %    models of the structures, defined by the rows of NN.
  9. %    These models are computed by the IV-method. The next rows
  10. %    of V are returned as NN'. The last row of V contains the logarithms 
  11. %    of the conditioning numbers of the IV-matrix for the structures in 
  12. %    question. The last column of V contains the number of data points in Z
  13. %
  14. %    ZE : the output-input data Z=[y u], with y and u as column vectors.This
  15. %    data is used for estimation. For multi-input systems u=[u1 u2 ... un]
  16. %    ZV : output-input data used to evaluate the models (could equal ZE).
  17. %    NN: Each row of NN is of the format [na nb nk], the orders and 
  18. %    delays associated with the corresponding model. (See HELP IV4 and 
  19. %    HELP STRUC.)
  20. %
  21. %    See HELP SELSTRUC for analysis of V. With
  22. %    V = ivstruc(Z,NN,p,maxsize) the calculation of conditioning numbers
  23. %    can be suppressed (p=0), and the speed/memory trade-off affected.
  24. %    See HELP AUXVAR.
  25.  
  26. %    L. Ljung 4-12-87,12-27-91
  27. %    Copyright (c) 1987-92 by the MathWorks, Inc.
  28. %    All Rights Reserved.
  29.  
  30. [Ncap,nz]=size(z); nu=nz-1;[nm,nl]=size(nn);
  31. na=nn(:,1);nb=nn(:,2:1+nu);nk=nn(:,2+nu:1+2*nu);
  32. [snbr,snbc]=size(nb);if any(any(nb==zeros(snbr,snbc))'),
  33. error('nb=0 is not feasible for ivstruc!'),end % Mod 911227
  34. nma=max(na);nbkm=max(nb+nk)-ones(1,nu);nkm=min(nk);
  35. n=nma+sum((nbkm-nkm))+nu;nbm=max(nb);
  36. % 
  37. % *** Some initial tests on the input arguments ***
  38. %
  39. maxsdef=idmsize(Ncap,n);
  40.  
  41. if nargin<5,maxsize=maxsdef;,end
  42. if nargin<4, p=1;,end
  43.  
  44. % *** construct instruments (see (7.111)-(7.112)) ***
  45. tha=arx(z,[nma nbm nkm],maxsize,1,0);
  46. NF=fstab([1 tha(1:nma)']);
  47.      
  48.      x=zeros(Ncap,1);rs=nma;
  49.      for k=1:nu
  50.           MF=[zeros(1,nkm(k)) tha(rs+1:rs+nbm(k))'];rs=rs+nbm(k);
  51.          x=x+filter(MF,NF,z(:,k+1));
  52.      end
  53. %
  54. % construct regression matrix
  55. %
  56. nmax=max(max([na+ones(nm,1) nb+nk]'))-1;
  57. M=floor(maxsize/n);nnuu=sum(nbkm)-sum(nkm)+nu;
  58. Rxx=zeros(nma);Ruu=zeros(nnuu);Rxu=zeros(nma,nnuu);Rxy=zeros(nma);
  59. Ruy=zeros(nnuu,nma); F=zeros(n,1);
  60. for k=nmax:M:Ncap
  61.      jj=(k+1:min(Ncap,k+M));
  62.      phix=zeros(length(jj),nma); phiy=phix; phiu=zeros(length(jj),nnuu);
  63.      for kl=1:nma, phiy(:,kl)=-z(jj-kl,1); phix(:,kl)=-x(jj-kl); end
  64.      ss=0;
  65.      for ku=1:nu
  66.             for kl=nkm(ku):nbkm(ku), phiu(:,ss+kl+1-nkm(ku))=z(jj-kl,ku+1);end
  67.              ss=ss+nbkm(ku)-nkm(ku)+1;
  68.      end
  69.      Rxy=Rxy+phix'*phiy;
  70.      Ruy=Ruy+phiu'*phiy;
  71.      Rxu=Rxu+phix'*phiu;
  72.      Ruu=Ruu+phiu'*phiu;
  73.      Rxx=Rxx+phix'*phix;
  74.      if nma>0, F(1:nma)=F(1:nma)+phix'*z(jj,1);end
  75.       F(nma+1:n)=F(nma+1:n)+phiu'*z(jj,1);
  76. end
  77. clear phiu, clear phix, clear phiy, 
  78. %
  79. % compute estimate
  80.  
  81. for j=1:nm
  82.     s=[];rs=0;
  83.     for ku=1:nu
  84.         s=[s,rs+nk(j,ku)-nkm(ku)+1:rs+nb(j,ku)+nk(j,ku)-nkm(ku)];
  85.         rs=rs+nbkm(ku)-nkm(ku)+1;
  86.     end
  87.     RR=[Rxy(1:na(j),1:na(j)) Rxu(1:na(j),s);Ruy(s,1:na(j)) Ruu(s,s)];
  88.     TH=RR\F([1:na(j) nma+s]);
  89.     a=TH(1:na(j))';b=TH(na(j)+1:length(TH))';f=fstab([1 a]);
  90.     [mzv,nzv] = size(zv(:,1));
  91.     yhat=zeros(mzv,nzv);rs=0;
  92.     for ku=1:nu
  93.     bdum=[zeros(1,nk(j,ku)) b(rs+1:rs+nb(j,ku))]; rs=rs+nb(j,ku);
  94.     yhat=yhat+filter(bdum,f,zv(:,ku));
  95.     end
  96.     V(1,j)=(zv(:,1)-yhat)'*(zv(:,1)-yhat)/(Ncap-nmax);
  97.     V(2:nl+1,j)=nn(j,:)';
  98.     if p~=0, V(nl+2,j)=log(cond(RR));,end
  99. end
  100. V(1,nm+1)=Ncap-nmax;
  101.  
  102.  
  103.