home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / H2SYN.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  6.1 KB  |  221 lines
  1. % function [k,g,norms,kfi,gfi,hamx,hamy] = h2syn(plant,nmeas,ncon,ricmethod)
  2. %
  3. %   Calculate the H2 optimal controller (K) and the closed loop
  4. %   system (G) for the linear fractional interconnection structure
  5. %   P.  NMEAS and NCON are the dimensions of the measurement outputs
  6. %   from P and the controller inputs to P. RICMETHOD determines the 
  7. %   method  used to solve the Riccati equations.  
  8. %
  9. %   inputs:
  10. %     PLANT   -   system interconnection structure
  11. %     NMEAS   -   measurements output to controller
  12. %     NCON    -   control inputs
  13. %     RICMETHOD   =   1     Eigenvalue decomposition.
  14. %                     2     Schur decomposition. (default)
  15. %
  16. %   outputs:
  17. %     K      -  H2 optimal controller
  18. %     G      -  closed-loop system with H2 optimal controller
  19. %     NORMS  -  norms of 4 different quantities, full information
  20. %               control cost (FI), output estimation cost (OEF), 
  21. %               direct feedback cost (DFL) and full control cost (FC).
  22. %               norms = [FI OEF DFL FC];
  23. %     KFI    -  full information/state feedback control law
  24. %     GFI    -  full information/state feedback closed-loop system
  25. %     HAMX   -  X Hamiltonian matrix
  26. %     HAMY   -  Y Hamiltonian matrix
  27. %
  28. %   See also: H2NORM, HINFSYN, HINFFI, HINFNORM, RIC_EIG and RIC_SCHR.
  29.  
  30. function [k,g,norms,kfi,gfi,hamx,hamy] = h2syn(plant,p2,m2,ricmethod)
  31.  
  32. if nargin ~= 3 & nargin ~=4,
  33.   disp('usage: [k,g,norms,kfi,gfi,hamx,hamy] = h2syn(plant,nmeas,ncon,ricmethod)')
  34.   return
  35. end
  36.  
  37. if nargin == 3,
  38.   ricmethod = 2;
  39. end
  40.  
  41. %    check the dimensions and assumptions of the problem.
  42. %    D11 must be zero, and here we will assume that D22 is
  43. %    also zero.  D12 and D21 must also be of appropriate rank.
  44.  
  45. [a,b,c,d] = unpck(plant);
  46. [anom,bnom,cnom,dnom] = unpck(plant);
  47. nx = max(size(a));
  48. [i,m] = size(b);
  49. [p,i] = size(c);
  50. if m <= m2,
  51.   disp('control input dimension incorrect')
  52.   return
  53. end
  54. if p <= p2,
  55.   disp('measurement output dimension incorrect')
  56.   return
  57. end
  58. p1 = p - p2;
  59. m1 = m - m2;
  60.  
  61. %    now introduce a unitary scaling (qin, qout) on the inputs and outputs,
  62. %    and a change of basis (tin, tout) on the controller measurements
  63. %    and controls.  This means that the d term meets the assumptions
  64. %    of the formulae.
  65.  
  66. %    select out the d12 and d21 terms.
  67.  
  68. d12 = d(1:p1,m1+1:m);
  69. d21 = d(p1+1:p,1:m1);
  70.  
  71. %    get qout and rin such that [qout*d12*rin] = [0;I]
  72. %    This is done by rearranging the result of a qr decomposition.
  73.  
  74. [q,r] = qr(d12);
  75. if rank(r) ~= m2,
  76.  disp('d12 does not have full column rank')
  77.  return
  78. end
  79. qout = [q(:,(m2+1):p1),q(:,1:m2)]';
  80. rin = eye(m2)/(r(1:m2,:));
  81.  
  82. %    get qin and rout such that [rout*d21*qin] = [0,I]
  83. %    Again this is done by rearranging the result of a qr decomposition.
  84.  
  85. [q,r] = qr(d21');
  86. if rank(r) ~= p2,
  87.     disp('d21 does not have full row rank')
  88.     return
  89.     end
  90. qin = [q(:,(p2+1):m1),q(:,1:p2)];
  91. rout = eye(p2)/(r(1:p2,:)');
  92.  
  93. %    now scale the inputs and outputs appropriately
  94. %    Note that the controller scaling (rin, rout) must
  95. %    be included in the calculation of k.  qin and qout
  96. %    do not affect the assumptions of the problem.
  97.  
  98. c = daug(qout,rout)*c;
  99. b = b*daug(qin,rin);
  100. d = daug(qout,rout)*d*daug(qin,rin);
  101.  
  102. %    now decompose the new system into its appropriate parts
  103.  
  104. c1 = c(1:p1,:);
  105. c2 = c(p1+1:p,:);
  106. b1 = b(:,1:m1);
  107. b2 = b(:,m1+1:m);
  108.  
  109. %  d12 and d21 are hardwired to the form obtained by the scaling and
  110. %  unitary transformation above.  This will help eliminate rounding errors.
  111.  
  112. d11 = d(1:p1,1:m1);
  113. d12 = [zeros(p1-m2,m2);eye(m2)];
  114. d21 = [zeros(p2,m1-p2),eye(p2)];
  115. d22 = d(p1+1:p,m1+1:m);
  116.  
  117. %    check for d11 being zero.
  118.  
  119. if any(any(d11 ~= 0)),
  120.   disp('d11 is non zero')
  121.   return
  122. end
  123.  
  124. %    now form the hamiltonian for X2
  125.  
  126. Ah = a - b2*d12'*c1;
  127. Eh = daug(eye(p1-m2),zeros(m2,m2));
  128.  
  129. if ricmethod == 1 | ricmethod == 2,
  130.     hamx = [Ah, -b2*b2'; -c1'*Eh*c1, -Ah'];
  131.     if ricmethod == 1,
  132.       [X1,X2,fail,xeig_min] = ric_eig(hamx,1e-13);
  133.     else
  134.       [X1,X2,fail,xeig_min] = ric_schr(hamx,1e-13);
  135.     end
  136.  
  137.     if fail ~= 0,
  138.       disp('Decomposition of X2 failed')
  139.       return
  140.     end
  141.     X2 = real(X2/X1);
  142. elseif ricmethod == 3,
  143. %    [X2,fail] = dgkfnewt(Ah,b2,c1'*Eh*c1,nx^2*1e-14*norm(Ah,inf),nx);
  144. %   X2 = real((X2+X2')/2);
  145. else
  146.     disp('invalid Riccati method')
  147.     return
  148. end
  149.  
  150. fprintf('minimum eigenvalue of X2: %e\n',xeig_min);
  151.  
  152. %    now form the hamiltonian for Y2
  153.  
  154. Aj = a - b1*d21'*c2;
  155. Ej = daug(eye(m1-p2),zeros(p2,p2));
  156.  
  157. if ricmethod == 1 | ricmethod == 2,
  158.     hamy = [Aj', -c2'*c2; -b1*Ej*b1', -Aj];
  159.     if ricmethod == 1,
  160.         [Y1,Y2,fail,yeig_min] = ric_eig(hamy,1e-13);
  161.     else
  162.         [Y1,Y2,fail,yeig_min] = ric_schr(hamy,1e-13);
  163.         end
  164.     Y2 = real(Y2/Y1);
  165.     if fail ~=0,
  166.       disp('Decomposition of Y2 failed')
  167.       return
  168.     end
  169. elseif ricmethod == 3,
  170. %    [Y2,fail] = dgkfnewt(Aj,c2',b1*Ej*b1',nx^2*1e-14*norm(Aj,inf),nx);
  171. %    Y2 = real((Y2+Y2')/2);
  172. else
  173.     disp('invalid Riccati method')
  174.     return
  175. end
  176.  
  177. fprintf('minimum eigenvalue of Y2: %e\n',yeig_min);
  178.  
  179. % form the general controller (as an LFT)
  180.  
  181. f2 = -d12'*c1 - b2'*X2;
  182. h2 = -b1*d21' - Y2*c2';
  183.  
  184. % include the scaling matrices into d22 feedback problem
  185. ak = a + h2*c2 + b2*f2 + h2*rout*d22*rin*f2;
  186.  
  187. % ak = a + h2*c2 + b2*f2 + h2*d22*f2;
  188. % bk = [-h2, (h2*d22 + b2)];
  189. % ck = [f2; (-d22*f2 - c2)];
  190. % dk = [zeros(m2,p2), eye(m2,m2); eye(p2,p2), -d22];
  191.  
  192. kgen = pck(ak,-h2,f2,zeros(m2,p2));
  193.  
  194. k = mmult(rin,mmult(kgen,rout));
  195. g = starp(plant,k,p2,m2);
  196.  
  197. %
  198. %  construct the full information controller and closed-loop system
  199. %
  200. if nargout >= 4
  201.   afi = anom;
  202.   bfi = bnom;
  203.   cfi = [cnom(1:p1,1:nx); eye(nx)];
  204.   dfi = [zeros(p1,m1) dnom(1:p1,m1+1:m); zeros(nx,m)];
  205.   pfi = pck(afi,bfi,cfi,dfi);
  206.   kfi = rin*f2;
  207.   gfi = starp(pfi,kfi,nx,m2);
  208. end
  209.  
  210. % construct the norm for the full info and output injection cases
  211.  
  212. if nargout >= 3
  213.   FI  = sqrt(trace((b1)'*X2*(b1)));
  214.   OEF = sqrt(trace((f2)*Y2*(f2)'));
  215.   DFL = sqrt(trace((h2)'*X2*(h2)));
  216.   FC  = sqrt(trace((c1)*Y2*(c1)'));
  217. end
  218. norms = [ FI OEF DFL FC];
  219. %
  220. % Copyright MUSYN INC 1991,  All Rights Reserved
  221.