home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / MKETAARX.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.5 KB  |  91 lines
  1. function eta=mketaarx(nn,par,lam,Tsamp)
  2. %MKETAARX constructs an ETA-model structure for (multivariate) ARX-models
  3. %    Auxiliary routine to ARX2TH
  4. %    TH = mketaarx([NA NB NK])
  5. %
  6. %    TH: The resulting structure in the THETA-format (see HELP THETA)
  7. %    NA: an ny|ny matrix whose i-j entry is the order of the
  8. %        polynomial (in the delay operator) that relates the
  9. %        j:th output to the i:th output
  10. %    NB: an ny|nu matrix whose i-j entry is the order of the 
  11. %        polynomial that related the j:th input to the i:th
  12. %        output
  13. %    NK: an ny|nu matrix whose i-j entry is the delay from
  14. %        the j:th input to the i:th output
  15. %    (ny: the number of outputs, nu: the number of inputs)
  16. %    
  17. %    Each row of [NA NB NK] is thus consistent with the way
  18. %    single output ARX structures is defined (See Help ARX)
  19. %    
  20. %    With TH = mketaarx([NA NB NK], PAR, LAM, T) the nominal values
  21. %    of the free parameters are set to PAR (default zeros). LAM will
  22. %       be the innovations covariance matrix, and T the sampling interval.
  23.  
  24. %    L. Ljung 10-2-90,25-1-92
  25. %    Copyright (c) 1990-92 by the MathWorks, Inc.
  26. %    All Rights Reserved.
  27.  
  28. if nargin<4, Tsamp=[];end
  29. if nargin<3, lam=[];end
  30. if nargin<2, par=[];end
  31.  
  32. [ny,nc]=size(nn);
  33. if isempty(Tsamp),Tsamp=1;end
  34. if isempty(lam),lam=eye(ny);end
  35. nu=(nc-ny)/2;
  36. na=nn(:,1:ny);
  37. if nu>0,nb=nn(:,ny+1:ny+nu);else nb=zeros(ny,1);end
  38. if nu>0,nk=nn(:,ny+nu+1:nc);else nk=ones(ny,1);end
  39. mna=max(max(na)');
  40. nb1=nb+(nk-1);
  41. mnb=max(max(nb1)');
  42.  
  43. mnk=min(min(nk)');
  44. nx=mna*ny+mnb*nu;
  45. C=zeros(ny,nx);
  46. for ky=1:ny
  47. for kr=1:ny
  48. for kc=1:na(ky,kr)
  49. C(ky,(kc-1)*ny+kr)=NaN;
  50. end
  51. end
  52. for kr=1:nu
  53. for kc=max(1,nk(ky,kr)):nb1(ky,kr)
  54. C(ky,(kc-1)*nu+mna*ny+kr)=NaN;
  55. end
  56. end
  57. end
  58.  
  59. if mna>0
  60. A0=C;A1=[eye(ny*(mna-1)),zeros(ny*(mna-1),ny+nu*mnb)];
  61. else A0=[];A1=[];
  62. end
  63. if mnb>0
  64. A3=[zeros(nu,nx);[zeros(nu*(mnb-1),ny*mna),eye(nu*(mnb-1)),zeros(nu*(mnb-1),nu)]];
  65. else A3=[];
  66. end
  67. A=[A0;A1;A3];
  68. if nu>0
  69. B=zeros(nx,nu);D=zeros(ny,nu);
  70. bc=(nk==0).*(nb>0);
  71. for kbc=1:ny
  72. for kbr=1:nu
  73. if bc(kbc,kbr), B(kbc,kbr)=NaN;end
  74. end
  75. end
  76. D=B(1:ny,1:nu);
  77. if mnb>0,B(ny*mna+1:ny*mna+nu,:)=eye(nu);end
  78. else B=[];end
  79. if nx==0, B=zeros(ny,nu);A=zeros(ny,ny);C=A;end
  80. if nu>0,if mna>0,D=zeros(ny,nu);end,else D=[];end
  81. K=zeros(nx,ny);
  82.  
  83. if mna>0;C=zeros(ny,nx);K(1:ny,1:ny)=eye(ny);end
  84. x0=zeros(nx,1);if nx==0,x0=zeros(ny,1);end
  85. ms=modstruc(A,B,C,D,K,x0);
  86. if isempty(par), par=zeros(1,sum(sum(na)')+sum(sum(nb)'));end
  87. [rms,cms]=size(ms);ms(1,cms)=-ms(1,cms);
  88. eta=ms2th(ms,'d',par,lam,Tsamp);
  89. eta(2,7)=39;eta(2,8)=3;
  90.  
  91.