home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / CANSTART.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.4 KB  |  69 lines
  1. function th = canstart(z,pobs,nu,oe)
  2. %CANSTART Produces a multivariable theta-model with initial parameter estimates
  3. %    THETA = canstart(Z,PSOBS,NU,OE)
  4. %    
  5. %    THETA: The resulting model, parametrized according to the pseudoobser-
  6. %              vability indices and with initial parameter estimates calculated
  7. %           from data. 
  8. %    Z :    The output-input data [y u], with each input/output component
  9. %           as a column vector.
  10. %    PSOBS: The pseudo-observability indices that define the multivariable
  11. %           state-space structure. This is a row vector with as many entries
  12. %           as there are outputs. Basically, PSOBS(k) gives the number of
  13. %           delays of output # k that are used. The order of the system
  14. %           is the sum of the psobs-indices. See also HELP CANFORM.
  15. %    NU:    The number of inputs.
  16. %    OE:    If OE='oe', then the K-matrix is THETA is fixed to zero, so
  17. %           that an output-error model structure is obtained. Otherwise
  18. %           (default) K is initialized at a value that gives a stable
  19. %           predictor.
  20. %
  21. %    THETA is just an initial estimate of the system, and should be further
  22. %          improved by TH = pem(z,THETA);
  23.  
  24. %    L. Ljung 21-9-90,11-2-92
  25. %    Copyright (c) 1990-92 by the MathWorks, Inc.
  26. %    All Rights Reserved.
  27.  
  28. if nargin<4,oe='no';end
  29. if nu==0, error('Sorry, this routine does not work for time series!'),end
  30. ny=length(pobs);
  31. [Ncap,nz]=size(z);
  32. if nz~=nu+ny, error(['The number of inputs and outputs in the data vector is not'...
  33. '\nconsistent with the number of observability indices and the number of inputs!']),end
  34. mp=max(pobs);
  35. r=cumsum(pobs);
  36. stateind=[]; newstate=[];
  37. for kl=1:ny
  38.     for ni=0:pobs(kl)-1
  39.         stateind=[stateind, kl+ni*ny];
  40.     end    
  41.     newstate=[newstate,kl+(ni+1)*ny];
  42. end
  43. nx=length(stateind);
  44. th = iv4mv(z,[mp*ones(ny,ny),mp*ones(ny,nu),ones(ny,nu)]);
  45. [a,b,c,d]=th2ss(th);
  46. [a,b,c]=minreal(a,b,c,d);
  47. [nxa,ndum]=size(a);
  48. co=ctrb(a,b);ob=obsv(a,c); ob=[ob;ob((nxa-1)*ny+1:nxa*ny,:)*a];
  49. hm=ob*co;
  50. apar=(hm(newstate,:)/hm(stateind,:))';
  51. bpar=(hm(stateind,1:nu))';
  52. ms=canform(pobs,nu);
  53. th=ms2th(ms,'d',[apar(:)' bpar(:)' zeros(1,ny*nx)]);
  54. if oe=='oe',
  55.     kind=[];
  56.     for kx=1:nx,for ky=1:ny
  57.         kind=[kind;kx,ky];
  58.     end,end
  59.     th=fixpar(th,'K',kind);
  60. else    
  61.     [a,b,c]=th2ss(th);
  62.     k=real(dlqe(a,eye(nx),c,eye(nx),eye(ny)));
  63.     k=(a*k)';
  64.     th=ms2th(ms,'d',[apar(:)' bpar(:)' k(:)']);
  65. end
  66. th(2,7)=24;
  67.  
  68.  
  69.