home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / THD2THC.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.7 KB  |  113 lines
  1. function thc=thd2thc(th,delays,nop)
  2. %THD2THC  converts a model to continuous-time form.
  3. %
  4. %    THC = thd2thc(TH)
  5. %
  6. %    TH: A discrete time model defined in the format described by HELP THETA
  7. %    THC: The continuous-time counterpart, stored in THETA format
  8. %
  9. %    If the model TH has extra delays from some inputs (nk>1), these are
  10. %    first removed. They should be appended to the continuous-time model
  11. %    as a pure time-delay (deadtime) (this is done by TH2FF).
  12. %    To have the delays approximated by THC, enter THC = thd2thc(TH,'del');
  13. %
  14. %    The covariance matrix P of TH is translated by the use of numerical
  15. %    derivatives. The step sizes used for the differences are given by the
  16. %     m-file nuderst. To inhibit the translation of P, which saves time,
  17. %    enter THC = thd2thc(TH,c,'nop'), where c is 'del' or 'nodel'. 
  18. %
  19. %    CAUTION: The transformation could be ill-conditioned. Negative
  20. %    real poles, poles in the origin and in 1 may cause problems. See the
  21. %    manual.
  22. %    See also THC2THD
  23.  
  24. %    L. Ljung 10-1-86,4-20-87,1-13-90
  25. %    Copyright (c) 1986-90 by the MathWorks, Inc.
  26. %    All Rights Reserved.
  27.  
  28.  
  29. T=th(1,2);
  30. if T<=0, error('This model is already denoted as continuous time!'),end
  31. if isthss(th), thc=th;thc(1,2)=-T;
  32. if thc(2,8)==3, if ~exist('minreal')
  33. disp('WARNING: Transformation of multioutput arx-models will be supported only if you have the CONTROL SYSTEMS TOOLBOX'),
  34. end,end
  35. return,end
  36. nu=th(1,3);nk=th(1,7+2*nu:6+3*nu);
  37. nn=th(1,3:6+2*nu);na=nn(2);nb=nn(3:2+nu);nc=nn(3+nu);nddd=nn(4+nu);nf=nn(5+nu:4+2*nu);
  38. % *** Set up default values ***
  39. if nargin<3,nop=[];end
  40. if nargin<2,delays=[];end
  41.  
  42. ndd=na+sum(nb)+nc+nddd+sum(nf);
  43.  
  44. if isempty(delays),delays='nodel';end
  45.  
  46. if delays(1)=='d',delays=0;else delays=1;end
  47. [par,P]=th2par(th);
  48. if norm(P)==0 | isempty(P) | ~isempty(nop), nkder=0;else nkder=[0:ndd];end
  49. dt=nuderst(par);
  50. for kder=nkder
  51. th1=th;if kder>0,th1(3,kder)=th(3,kder)+dt(kder);end    
  52.  
  53. [a,b,c,d,f]=th2poly(th1);
  54. b(nu+1,1:length(c))=c;
  55. f(nu+1,1:length(d))=d;
  56. nk(nu+1)=0;nb(nu+1)=nc+1;nf(nu+1)=nddd;
  57. ss=1;
  58. thb=[];thcc=[];thd=[];thf=[];ncn=0;ndn=0;
  59. for k=[nu+1,1:nu]
  60. den=conv(a,f(k,1:nf(k)+1)); 
  61. if delays~=0,nkmod=max(nk(k),1);else nkmod=1;end
  62. num=b(k,nkmod:nb(k)+nk(k));
  63. nd=length(den);, nn=length(num);
  64. n=max(nd,nn);
  65. de=zeros(1,n);, de(1:nd)=den;
  66. if kder==0
  67. if abs(de(n))<eps, disp('WARNING: Pole in the origin. Result may be unreliable!'),end
  68. rp=roots(de);
  69. if any(abs(rp-1))<eps disp('WARNING: Pure integration. Result may be unreliable!'),end
  70. if any((rp==real(rp)).*(real(rp)<0)),disp('WARNING: Pole on the negative real axis. Result may be unreliable!'),end
  71. end
  72. nume=zeros(1,n);, nume(1:nn)=num;
  73. A=[-de(2:n);eye(n-2,n-1)]; % Transform to state-space
  74. B=eye(n-1,1);
  75. C=nume(2:n)-nume(1)*de(2:n);
  76. D=nume(1); 
  77. s=real(logm([[A B];zeros(1,n-1) eye(1)]))/T;  % Sample
  78. Ac=s(1:n-1,1:n-1);
  79. Bcc=s(1:n-1,n);
  80. ff=poly(Ac);                % Transform to i/o form
  81. bb=poly(Ac-Bcc*C)+(D-1)*ff;
  82. if nk(k)>0,bb=bb(2:length(bb));nkn(ss)=1;else nkn(ss)=0;end
  83.  
  84. if k<nu+1,
  85.     nbn(ss)=length(bb);nfn(ss)=length(ff)-1;ss=ss+1; 
  86.     thb=[thb bb]; thf=[thf ff(2:length(ff))];
  87. end
  88. if k==nu+1, 
  89.     ncn=length(bb)-1;ndn=length(ff)-1;
  90.     thcc=bb(2:length(bb)); thd=ff(2:length(ff));
  91. end
  92.  
  93. end
  94. thh=[thb thcc thd thf];
  95. if kder==0,
  96. [ntr,ntc]=size(th);ntcc=max([ntc length(thh)]);
  97. thc=zeros(3+length(thh),ntcc);
  98. thc(1:ntr,1:ntc)=th;
  99. thc(1,2)=-T;
  100. if nu>0,thc(1,3:6+3*nu)=[nu 0 nbn ncn ndn nfn nkn];
  101. else thc(1,3:6)=[nu 0 ncn ndn];end
  102. thc(3,1:length(thh))=thh;thh0=thh;
  103. end
  104. if kder>0,
  105. dif(:,kder)=(thh-thh0)'/dt(kder);
  106. end
  107. end
  108.  
  109. if ~isempty(dif),Pc=dif*P*dif';else Pc=zeros(length(thh),length(thh));end
  110. thc(4:length(thh)+3,1:length(thh))=Pc;
  111. if delays, thc(length(thh)+4,1:length(nk))=(max(nk,1)-1)*T;end
  112.  
  113.