home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / IMP2SS.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.4 KB  |  129 lines
  1. function [a,b,c,d,totbnd,hsv] = imp2ss(Z1,Z2,Z3,Z4,Z5)
  2. % [A,B,C,D,TOTBND,HSV] = IMP2SS(Y) or
  3. % [A,B,C,D,TOTBND,HSV] = IMP2SS(Y,TS,NU,NY,TOL) or
  4. % [SS_,TOTBND,HSV] = IMP2SS(IMP_) 
  5. % [SS_,TOTBND,HSV] = IMP2SS(IMP_,TOL) produces an approximate 
  6. % state-space  realization of a given impulse response 
  7. %                 IMP_=MKSYS(Y,TS,NU,NY,'imp')
  8. % using the Hankel SVD method proposed by S. Kung (Proc. Asilomar
  9. % Conf. on Circ. Syst. & Computers, 1978).  A continuous-time 
  10. % realization is computed via the the inverse Tustin transform if
  11. % TS is positive; otherwise, a discrete-time realization is returned.
  12. %
  13. %  INPUTS:  Y --- impulse response H1,...,HN stored rowwise
  14. %                 Y=[H1(:)'; H2(:)'; H3(:)'; ...; HN(:)']
  15. %  OPTIONAL INPUTS:
  16. %           TS  ---  sampling interval (default  TS=0 --- discrete-time)
  17. %           NU  ---  number of inputs  (default NU=1)
  18. %           NY  ---  number of outputs (default NY= size(Y)*[1;0]/nu)
  19. %           TOL ---  Hinfnorm of error bound (default TOL=0.01*S(1))
  20. %  OUTPUTS: (A,B,C,D) state-space realization
  21. %           TOTBND  Infinity norm error bound  2*Sum([S(NX+1),S(NX+2),...])
  22. %           HSV  Hankel singular values [S(1);S(2);... ]
  23.  
  24. % R. Y. Chiang & M. G. Safonov 8/91
  25. % Copyright (c) 1991 by the MathWorks, Inc.
  26. % All Rights Reserved.
  27.  
  28. nag1 = nargin;
  29. nag2 = nargout;
  30.  
  31. inargs='(y,ts,nu,ny,tol)';
  32. eval(mkargs(inargs,nargin,'imp'))
  33.  
  34. y=y';
  35. [mm,nn] = size(y);
  36. if nn==1,
  37.   y=y';
  38.   [mm,nn] = size(y);
  39. end
  40.  
  41. % Do defaults for TS,NU,NY and TOL
  42. if nargin<2,ts=0;end
  43. if nargin<3,nu=1;end
  44. if nargin<4,
  45.   if rem(mm,nu)~=0, 
  46.     error('Incompatible Dimensions---The row dimension of Y must be divisible by NU'),
  47.   end
  48.   ny=mm/nu;
  49. end
  50.  
  51.  
  52. % Check dimensional compatibility:
  53. if ny*nu~=mm,
  54.   error('Incompatable Dimensions---Must have NY*NU = (no. of columns of Y)')
  55. end
  56.  
  57. z=zeros(ny,nu);
  58.  
  59. % Get D matrix D=H(0)
  60. d=z;
  61. d(:) = y(:,1);
  62.  
  63. % Overwrite Y with Y= [H1 H2 ...Hnn]
  64. tmp=y;
  65. m=ny;
  66. n=nn*nu;
  67. y=zeros(m,n);
  68. y(:)=tmp(:);
  69.  
  70. % Build Hankel matrix H = [H(1)   H(2) ... H(nn)
  71. %                          H(2)   H(3) ...   0
  72. %                                  ...
  73. %                          H(nn)    0        0 ]
  74. if nu==1 & ny==1,
  75.   h=hankel(y(2:nn));
  76. else
  77.   for i=1:nn
  78.     y=[y(:,nu+1:n) z];   % y= [H(i), H(i+1),...,H(nn),0,...,0]
  79.     h=[h;y];
  80.   end    
  81. end
  82. [rh,ch]=size(h);
  83.  
  84. % Compute the SVD of the Hankel
  85. [u,hsv,v] = svd(h);
  86. hsv=diag(hsv);
  87. ns=size(hsv);
  88. if nargin<5,tol=0.01*hsv(1);end
  89.  
  90. % Determine NX and TOTBND
  91. tail=[conv(ones(ns,1)',hsv') 0];
  92. tail=tail(ns:ns+ns);
  93. nx=find(2*tail<=tol*ones(tail));
  94. nx=nx(1)-1;
  95. totbnd=2*tail(nx+1);  % TOTBND = 2*(HSV(NX+1)+HSV(NX+2)+...+HSV(NS))
  96.  
  97.  
  98. % 
  99. % Truncate and partition singular vector matrices U,V
  100. u1 = u(1:rh-ny,1:nx);
  101. v1 = v(1:ch-nu,1:nx);
  102. u2 = u(ny+1:rh,1:nx);
  103. u3 = u(rh-ny+1:rh,1:nx);
  104. %
  105. ss = sqrt(hsv(1:nx));
  106. invss = ones(nx,1)./ss;
  107. %
  108.  
  109. % Kung's formula for the reduced model:
  110. %
  111. % The following is equivalent to UBAR = PINV(U)*[U2; 0]:
  112. ubar = u1'*u2;
  113. a = ubar.*(invss*ss');   % A = INV(DIAG(SS))*UBAR*DIAG(SS)
  114. b = (ss*ones(1,nu)).*v1(1:nu,:)';  % B = DIAG(SS)*V1(1:NU,:)'
  115. c = u1(1:ny,:).*(ones(ny,1)*ss');  % C = U1(1:NY,:)*DIAG(SS)
  116. %
  117. if ts>0,
  118.   [a,b,c,d]=bilin(a,b,c,d,-1,'Tustin',ts); % convert to continuous
  119. end  
  120.  
  121. if nag2<4;
  122.    ss_ = mksys(a,b,c,d);
  123.    a = ss_;
  124.    b=totbnd;
  125.    c=hsv;
  126. end
  127. %
  128. % -------- End of IMP2SS.M % RYC/MGS
  129.