home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / HINF.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  16.6 KB  |  576 lines
  1. function [acp,bcp,ccp,dcp,acl,bcl,ccl,dcl,hinfo,ak,bk1,bk2,ck1,ck2,dk11,dk12,dk21,dk22]=hinf(Z1,Z2,Z3,Z4,Z5,Z6,Z7,Z8,Z9,Z10,Z11,Z12,Z13,Z14)
  2.  
  3. %     << 4-block Optimal H-inf Controller (all solution formulae) >>
  4. %               (Safonov, Limebeer and Chiang, 1989 IJC)
  5. %    
  6. %   [SS_CP,SS_CL,HINFO,TSS_K]=HINF(TSS_P,SS_U,VERBOSE) computes the
  7. %   H-infinity controller F(s) and the controller parameterization K(s)
  8. %   using the "loop-shifting formulae". Given any stable U(s) with 
  9. %   || U ||_inf <=1, F(s) is formed by wrapping feedback U(s) around K(s).
  10. %
  11. %   Required inputs --
  12. %      Augmented plant P(s): TSS_P=MKSYS(A,B1,B2,C1,C2,D11,D12,D21,D22)
  13. %
  14. %   Optional inputs --
  15. %          Stable contraction U(s): SS_U=MKSYS(AU,BU,CU,DU) (default: U=0)
  16. %          VERBOSE:  if  1,  a verbose display results (default),
  17. %                    if  0, HINF operates silently
  18. %
  19. %   Output data: Controller F(s):        SS_CP=MKSYS(ACP,BCP,CCP,DCP)
  20. %                Closed-Loop Ty1u1(s):   SS_CL=MKSYS(ACL,BCL,CCL,DCL) 
  21. %                hinfo = (hinflag,RHP_cl,lamps_max) ("hinflag":existence flag)
  22. %                All-solutions controller parameterization K(s):
  23. %                     TSS_K=MKSYS(A,BK1,BK2,CK1,CK2,DK11,DK12,DK21,DK22)
  24.  
  25. % R. Y. Chiang & M. G. Safonov 9/18/1988
  26. % Copyright (c) 1988 by the MathWorks, Inc.
  27. % All Rights Reserved.
  28. % ------------------------------------------------------------------------
  29.  
  30. % COMMENT:  This function also supports the following format for
  31. % input and output data which is not described above.  Possible
  32. % alternative usages include the following:
  33. %     [ACP,BCP,CCP,DCP,ACP,BCL,CCL,DCL,HINFO,AK,BK1,BK2,CK1,CK2,DK11,DK12,...
  34. %          DK21,DK22]=HINF(A,B1,B2,C1,C2,D11,D12,D21,D22,AU,BU,CU,DU,VERBOSE)
  35. %     [ACP,BCP,CCP,DCP,ACP,BCL,CCL,DCL,HINFO,AK,BK1,BK2,CK1,CK2,DK11,DK12,...
  36. %          DK21,DK22]=HINF(A,B1,B2,C1,C2,D11,D12,D21,D22,VERBOSE)
  37. % In all cases, the input data for U(s) and for VERBOSE is optional with
  38. % defaults U(s)=0 and VERBOSE=1. If no output arguments are supplied the
  39. % variables (ss_cp,ss_cl,acp,bcp,ccp, dcp,acl,bcl,ccl,dcl) are automatically
  40. % returned in the main workspace.
  41.  
  42. %
  43. % ----- Initialize all the flags:
  44. %
  45. hexist  = ' OK ';               
  46. D11flag = ' OK ';
  47. Pexist  = ' OK ';
  48. Pflag   = ' OK ';
  49. Sexist  = ' OK ';
  50. Sflag   = ' OK ';
  51. PSflag  = ' OK ';
  52. %
  53.  
  54. %%%%
  55.  
  56. %  Expand input arguments or, if none present, load them from main workspace
  57. nag1 = nargin;
  58. nag2 = nargout;
  59. sysuflag=0;       % This flag will be set to one if U(s) data is input
  60.  
  61.   inargs='(A,B1,B2,C1,C2,D11,D12,D21,D22,AU,BU,CU,DU,disp_hinf)';
  62.   eval(mkargs(inargs,nag1,'tss,ss'))  % Get input args, NARGIN, and XSFLAG 
  63.   nag1=nargin; % NARGIN may have changed if U(s) is given as system SS_U 
  64.                % or if plant is given is system TSS_P
  65.   if nag1==9
  66.      % Case A,B1,B2,C1,C2,D11,D12,D21,D22
  67.      disp_hinf=1;
  68.   elseif nag1==10   % If only 10 input args after expansion
  69.      % Case A,B1,B2,C1,C2,D11,D12,D21,D22,disp_hinf
  70.      disp_hinf=AU;
  71.   elseif nag1==13
  72.      % Case A,B1,B2,C1,C2,D11,D12,D21,D22,AU,BU,CU,DU
  73.      sysuflag=1;
  74.      disp_hinf=1;
  75.   elseif nag1==14
  76.      % Case A,B1,B2,C1,C2,D11,D12,D21,D22,AU,BU,CU,DU,disp_hinf
  77.      sysuflag=1;
  78.   else
  79.      error('Incompatible number or type of input arguments')
  80.   end
  81.  
  82. %  Create SYSP and DIMP, just in case not already present
  83. sysp=[A B1 B2;C1 D11 D12;C2 D21 D22];
  84. [dimx,dimu1]=size(B1);
  85. [dimy1,dimu2]=size(D12);
  86. [dimy2,dimu2]=size(D22);
  87. dimp=[dimx dimu1 dimu2 dimy1 dimy2];
  88.  
  89. %%%
  90.  
  91. %
  92. if disp_hinf == 1;
  93.    disp('   ')
  94.    disp('                << H-inf Optimal Control Synthesis >>')
  95.    disp(' ')
  96.    disp('            Computing the 4-block H-inf optimal controller ')
  97.    disp('          using the S-L-C loop-shifting/descriptor formulae ')
  98. end
  99. %
  100.  
  101. %  If VERBOSE mode set, display message about U(s)
  102. if sysuflag==0,
  103.    if disp_hinf == 1   
  104.       disp('  ')
  105.       disp('     Solving for the H-inf controller F(s) using U(s) = 0 (default)');
  106.    end
  107. else
  108.    if disp_hinf == 1
  109.       disp('  ')
  110.       disp('     Solving for the H-inf controller F(s) using U(s) = (AU,BU,CU,DU)');
  111.    end
  112. end
  113. %
  114. % --------------------------------------------------------
  115. %
  116. % ---------- Checking stabilizability and detectability:
  117. %
  118. [junk,junk,junk,junk,aa,ba,ca,da,m]=...
  119.    stabproj(A,B2,C2,D22);  % get antistable part
  120. if m<min(size(A));
  121.    na=min(size(aa));
  122.    hsv = hksv(-aa-eye(na),ba,ca); 
  123.    % Note: -EYE(NA) term in above prevents Inf grammians
  124.    if min(hsv)<(1.e-15)*max(hsv)
  125.       error('Your system is not stabilizable and/or detectable !!');
  126.    end
  127. end
  128.  
  129. % ----------------------------------------------------------------------------
  130. % ------------------------------------------- Pre-process the augmented plant:
  131. %
  132. % ------------------------------------------- 1) Zero out the D22 term:
  133. %
  134. [mD22,nD22] = size(D22);
  135. d22 = zeros(mD22,nD22);
  136. %
  137. % ------------------------------------------- 2) Do the scaling:
  138. %
  139. [p1,m2] = size(D12);
  140. [uu,sig12,vv] = svd(D12);
  141. if m2>p1, error('      ERROR:  Must have dim(y1) >= dim(u2)'),end
  142. if sig12(m2,m2)<eps, error('ERROR:  D12 must be full rank'),end
  143. u1 = uu(:,1:m2);
  144. u2 = uu(:,m2+1:p1);
  145. SU2 = vv*inv(sig12(1:m2,1:m2));
  146. SY1 = [u2';u1'];
  147. %
  148. [p2,m1] = size(D21);
  149. [uuu,sig21,vvv] = svd(D21);
  150. if p2>m1, error('      ERROR:  Must have dim(u1) >= dim(y2)'),end
  151. if sig21(p2,p2)<eps, error('      ERROR:  D21 must be full rank'),end
  152. v1 = vvv(:,1:p2);
  153. v2 = vvv(:,p2+1:m1);
  154. SY2 = inv(sig21(1:p2,1:p2))*uuu';
  155. SU1 = [v2 v1];
  156. %
  157. b1 = B1*SU1;
  158. b2 = B2*SU2;
  159. c1 = SY1*C1;
  160. c2 = SY2*C2;  
  161. d11 = SY1*D11*SU1;
  162. d12 = SY1*D12*SU2;
  163. d21 = SY2*D21*SU1;
  164. %
  165. if disp_hinf == 1
  166.    disp('')
  167.    disp('     Solving Riccati equations and performing H-infinity')
  168.    disp('     existence tests:')
  169. end
  170. %
  171. % ----------------------------------- 3) Checking the frequency @ inf:
  172. %
  173. d1111 = d11(1:p1-m2,1:m1-p2);
  174. d1121 = d11(p1-m2+1:p1,1:m1-p2);
  175. d1112 = d11(1:p1-m2,m1-p2+1:m1);
  176. d1122 = d11(p1-m2+1:p1,m1-p2+1:m1);
  177. mxsvd11r = max([svd([d1111 d1112]);0]);
  178. mxsvd11c = max([svd([d1111;d1121]);0]);
  179. gamainf = max([mxsvd11r,mxsvd11c]);
  180. %
  181. % ----------------------------------- Fix D11 term:
  182. %
  183. if m1 == p2 | p1 == m2
  184.    fopt = -d1122;
  185. else
  186.    fopt  = -(d1122+d1121*inv(eye(m1-p2)-d1111'*d1111)*d1111'*d1112);
  187. end
  188. %
  189. % ---- Fix the "Singular FOPT (if det(I+fopt*SY2*D22*SU2) = 0):
  190. %
  191. imax = 23;
  192. foptscale = 2^(-imax);
  193. fopt1 = fopt;
  194. for i = 1:imax+1
  195.     temp = eye(size(fopt)*[1;0])+fopt1*SY2*D22*SU2;
  196.     svtemp = svd(temp);
  197.     nulltemp = size(find(svtemp<1.e-8))*[0;1]; 
  198.     if nulltemp > 0
  199.        svfopt = svd(fopt);
  200.        fopt1 = (1+foptscale*(1-gamainf)/svfopt(1,1))*fopt;
  201.     else
  202.        break
  203.     end
  204.     foptscale = 2*foptscale;
  205. end
  206. fopt = fopt1;
  207. %
  208. d11ok = 0;
  209. if gamainf < 1 & foptscale < 1.99
  210.     if disp_hinf == 1
  211.        disp('        1.  Is D11 small enough?                       OK')
  212.     end
  213. else
  214.     if disp_hinf == 1
  215.        disp('        1.  Is D11 small enough?                       FAIL')
  216.     end
  217.     D11flag = 'FAIL'; hexist = 'FAIL';
  218. end   
  219. %
  220. a   = A  + b2*fopt*c2;
  221. b1  = b1 + b2*fopt*d21;
  222. c1  = c1 + d12*fopt*c2;
  223. d11 = d11+ d12*fopt*d21;
  224. %
  225. % ----------------------------------- 4) Zero out D11 term via loop shifting:
  226. %
  227. X = d11;
  228. [rX,cX] = size(X);
  229. IXX1 = eye(rX)-X*X';
  230. IXX2 = eye(cX)-X'*X;
  231. IIXX1 = inv(IXX1);
  232. a   = a+b1*X'*IIXX1*c1;
  233. b2  = b2 + b1*X'*IIXX1*d12;
  234. b1  = b1*IXX2^(-0.5);
  235. c2  = c2+d21*X'*IIXX1*c1;
  236. c1  = IXX1^(-0.5)*c1;
  237. d22 = d21*X'*IIXX1*d12;
  238. [rd11,cd11] = size(d11);
  239. d11 = zeros(rd11,cd11);
  240. d12 = IXX1^(-0.5)*d12;
  241. d21 = d21*IXX2^(-0.5);
  242. %
  243. % ----------------------------------- 5) Zero out the new D22 term:
  244. %
  245. d22_4 = d22;
  246. [rd22,cd22] = size(d22);
  247. d22 = zeros(rd22,cd22);
  248. %
  249. % ----------------------------------- 6) Do the scaling again:
  250. %
  251. [p1,m2] = size(d12);
  252. [uu,sig12,vv] = svd(d12);
  253. u1 = uu(:,1:m2);
  254. u2 = uu(:,m2+1:p1);
  255. su2 = vv*inv(sig12(1:m2,1:m2));
  256. sy1 = [u2';u1'];
  257. %
  258. [p2,m1] = size(d21);
  259. [uuu,sig21,vvv] = svd(d21);
  260. v1 = vvv(:,1:p2);
  261. v2 = vvv(:,p2+1:m1);
  262. sy2 = inv(sig21(1:p2,1:p2))*uuu';
  263. su1 = [v2 v1];
  264. %
  265. b1 = b1*su1;
  266. b2 = b2*su2;
  267. c1 = sy1*c1;
  268. c2 = sy2*c2;  
  269. d11 = sy1*d11*su1;
  270. d12 = sy1*d12*su2;
  271. d21 = sy2*d21*su1;
  272. % 
  273. % ----------------------------------------------------------------------------
  274. % 
  275. % ---------------------------------------- Start to compute the controller:
  276. siz = [1 0]';
  277. n = size(A)*siz;
  278. [rd12,cd12] = size(d12);
  279. c1til = (eye(rd12)-d12*d12')*c1;
  280. [rd21,cd21] = size(d21);
  281. b1til = b1*(eye(cd21)-d21'*d21);
  282. %                                                      
  283. % ---------------------------------------- Solving the P Riccati:
  284. %
  285. if disp_hinf == 1
  286.    disp('        2.  Solving state-feedback (P) Riccati ...')
  287. end
  288. %          
  289. % ------ P =  0, if the Hamiltonian A matrix is stable and D12 is square
  290. %
  291. [rd12,cd12] = size(D12);
  292. ar1 = a-b2*d12'*c1;
  293. if rd12 == cd12
  294.    lamar1 = eig(ar1);
  295.    if max(real(lamar1)) < 0
  296.        ar1flag=1;
  297.    else
  298.        ar1flag=0;
  299.    end
  300. else 
  301.    ar1flag=0;
  302. end
  303. if ar1flag == 1
  304.    P2 = zeros(n,n);
  305.    P1 = eye(n);
  306.    P  = zeros(n,n);
  307.    lamp_inf = lamar1;
  308.    wellposed = 'TRUE ';
  309. else
  310.    q1  = c1til'*c1til;
  311.    r1  = -(b1*b1'-b2*b2');
  312.    [P1,P2,lamp,Perr,wellposed] = aresolv(ar1,q1,r1);
  313.    [vpinf,dpinf] = eig(ar1*P1-b2*b2'*P2,P1);
  314.    lamp_inf = diag(dpinf);
  315.    for i = 1:n
  316.      if lamp_inf(i) == inf | lamp_inf(i) == NaN
  317.         lamp_inf(i) = 1.e15;
  318.      end
  319.    end
  320. end
  321. if wellposed=='FALSE',
  322.     if disp_hinf == 1
  323.        disp('             a.  No Hamiltonian jw-axis roots?         FAIL')
  324.     end
  325.     Pexist='FAIL'; hexist = 'FAIL';
  326. else
  327.     if disp_hinf == 1
  328.        disp('             a.  No Hamiltonian jw-axis roots?         OK')
  329.     end
  330. end
  331. if max(real(lamp_inf)) > 0
  332.     if disp_hinf == 1
  333.        disp('             b.  A-B2*F stable (P >= 0)?               FAIL')
  334.     end
  335.     Pflag ='FAIL'; hexist = 'FAIL';
  336. else
  337.     if disp_hinf == 1
  338.        disp('             b.  A-B2*F stable (P >= 0)?               OK')
  339.     end
  340. end
  341. %                                                      
  342. % ---------------------------------------- Solving the S Riccati:
  343. %
  344. if disp_hinf == 1
  345.    disp('        3.  Solving output-injection (S) Riccati ...')
  346. end
  347. %
  348. % ------ S =  0, if the Hamiltonian A matrix is stable and D21 is square
  349. %
  350. [rd21,cd21] = size(D21);
  351. ar2 = a-b1*d21'*c2;
  352. lamar2 = eig(ar2);
  353. if rd21 == cd21 & max(real(lamar2)) < 0
  354.       S2 = zeros(n,n); Serr = S2;
  355.       S1 = eye(n); S  = zeros(n,n);
  356.       lams_inf = lamar2;
  357.       wellposed = 'TRUE ';
  358. else
  359.    q2  = b1til*b1til';
  360.    r2  = -(c1'*c1-c2'*c2);
  361.    [S1,S2,lams,Serr,wellposed] = aresolv(ar2',q2,r2);
  362.    S1 = S1';
  363.    S2 = S2';
  364.    [vsinf,dsinf] = eig(S1*ar2-S2*c2'*c2,S1);
  365.    lams_inf = diag(dsinf);
  366.    for i = 1:n
  367.       if lams_inf(i) == inf | lams_inf(i) == NaN
  368.          lams_inf(i) = 1.e15;
  369.       end
  370.    end
  371. end
  372. if wellposed=='FALSE',
  373.     if disp_hinf == 1
  374.        disp('             a.  No Hamiltonian jw-axis roots?         FAIL')
  375.     end
  376.     Sexist='FAIL'; hexist = 'FAIL';
  377. else
  378.     if disp_hinf == 1
  379.        disp('             a.  No Hamiltonian jw-axis roots?         OK')
  380.     end
  381. end
  382. if max(real(lams_inf)) > 0
  383.     if disp_hinf == 1
  384.        disp('             b.  A-G*C2 stable (S >= 0)?               FAIL')
  385.     end
  386.     Sflag ='FAIL'; hexist = 'FAIL';
  387. else
  388.     if disp_hinf == 1
  389.        disp('             b.  A-G*C2 stable (S >= 0)?               OK')
  390.     end
  391. end
  392. %
  393. [vps,lamps] = eig(P2'*S2',P1'*S1');
  394. lamps = diag(lamps);
  395. for i = 1:n
  396.    if lamps(i) == inf | lamps(i) == NaN
  397.       lamps(i) = 1.e15;
  398.    end
  399. end
  400. %
  401. lamps_max = max(real(lamps));
  402. if lamps_max > 1
  403.     if disp_hinf == 1
  404.        disp('        4.  max eig(P*S) < 1 ?                         FAIL')
  405.     end
  406.     PSflag ='FAIL'; hexist = 'FAIL';
  407. else
  408.     if disp_hinf == 1
  409.        disp('        4.  max eig(P*S) < 1 ?                         OK')
  410.     end
  411. end
  412. %
  413. % --------------------------------------- Checking the H-inf solution:
  414. %
  415. if disp_hinf == 1
  416.    disp('    -------------------------------------------------------');
  417. end
  418. if hexist == 'FAIL',
  419.    if disp_hinf == 1
  420.          disp('       NO STABILIZING CONTROLLER MEETS THE SPEC. !!')
  421.    end
  422. else
  423.    if disp_hinf == 1
  424.          disp('       all tests passed -- computing H-inf controller ...')
  425.    end
  426. end
  427. %
  428. % --------------- Controller parameterization K(s) in descriptor form:
  429. %
  430. E = S1*P1-S2*P2;
  431. %
  432. ak   = a-b2*d12'*c1-b1*d21'*c2;
  433. ak   = S1*ak*P1+S2*ak'*P2 + S1*(b1til*b1til'-b2*b2')*P2 + ...
  434.        S2*(c1til'*c1til-c2'*c2)*P1;
  435. bk1  = S2*c2' + S1*b1*d21'; 
  436. bk2  = S1*b2 + S2*c1'*d12;
  437. ck1  = -(b2'*P2+d12'*c1*P1);
  438. ck2  = -(c2*P1 + d21*b1'*P2);
  439. [rbk1,cbk1] = size(bk1);
  440. [rbk2,cbk2] = size(bk2);
  441. [rck1,cck1] = size(ck1);
  442. [rck2,cck2] = size(ck2);
  443. dk11 = zeros(rck1,cbk1);
  444. dk12 = eye(rck1);
  445. dk21 = eye(rck2);
  446. dk22 = zeros(rck2,cbk2);
  447. % 
  448. % ------------------------------------------------------------------------
  449. %
  450. % ------------------------------------ Post-process the controller:
  451. %
  452. % ------------------------------------ 1) Reverse the controller scaling:
  453. %
  454. bk1  = bk1*sy2;
  455. ck1  = su2*ck1;
  456. dk11 = su2*dk11*sy2;
  457. dk12 = su2*dk12;
  458. dk21 = dk21*sy2;
  459. %
  460. % ------------------------------------ 2) Shifting D22_4 term:
  461. %
  462. temp = inv(eye(size(dk11)*[1;0]) + dk11*d22_4);
  463. ak   = ak-bk1*d22_4*temp*ck1;
  464. bk2  = bk2-bk1*d22_4*temp*dk12;
  465. ck2  = ck2-dk21*d22_4*temp*ck1;
  466. ck1  = temp*ck1;
  467. dk22 = dk22-dk21*d22_4*temp*dk12;
  468. dk12 = temp*dk12;
  469. dk11 = temp*dk11;
  470. temp = eye(size(d22_4)*[1;0])-d22_4*dk11;
  471. bk1  = bk1*temp;
  472. dk21 = dk21*temp;
  473. %
  474. % ------------------------------------ 3) Reverse the "fopt" term:
  475. %
  476. dk11 = dk11 + fopt;
  477. %
  478. % ------------------------------------ 4) Reverse the controller scaling again:
  479. %
  480. bk1  = bk1*SY2;
  481. ck1  = SU2*ck1;
  482. dk11 = SU2*dk11*SY2;
  483. dk12 = SU2*dk12;
  484. dk21 = dk21*SY2;
  485. %
  486. % ------------------------------------ 5) Shifting D22 term:
  487. %
  488. temp = inv(eye(size(dk11)*[1;0]) + dk11*D22);
  489. ak   = ak-bk1*D22*temp*ck1;
  490. bk2  = bk2-bk1*D22*temp*dk12;
  491. ck2  = ck2-dk21*D22*temp*ck1;
  492. ck1  = temp*ck1;
  493. dk22 = dk22-dk21*D22*temp*dk12;
  494. dk12 = temp*dk12;
  495. dk11 = temp*dk11;
  496. temp = eye(size(D22)*[1;0])-D22*dk11;
  497. bk1  = bk1*temp;
  498. dk21 = dk21*temp;
  499. %
  500. % ----------------  Putting descriptor controller K(s) in state space form:
  501. %
  502. [ak,bk,ck,dk] = des2ss(ak,[bk1,bk2],[ck1;ck2],...
  503.                          [dk11,dk12;dk21,dk22],E);
  504. bk1=bk(:,1:cbk1);
  505. bk2=bk(:,cbk1+1:cbk1+cbk2);
  506. ck1=ck(1:rck1,:);
  507. ck2=ck(rck1+1:rck1+rck2,:);
  508. dk11=dk(1:rck1,1:cbk1);
  509. dk12=dk(1:rck1,cbk1+1:cbk1+cbk2);
  510. dk21=dk(rck1+1:rck1+rck2,1:cbk1);
  511. dk22=dk(rck1+1:rck1+rck2,cbk1+1:cbk1+cbk2);
  512. %
  513. % --------------------- computing controller F(s) = (acp,bcp,ccp,dcp):
  514. %
  515. sysk = [ak,bk1,bk2;ck1,dk11,dk12;ck2,dk21,dk22];
  516. dimk = [size(ak)*[1 0]',size(bk1)*[0 1]',size(bk2)*[0 1]',...
  517.             size(ck1)*[1 0]',size(ck2)*[1 0]'];
  518. if sysuflag,
  519.    [acp,bcp,ccp,dcp] = lftf(sysk,dimk,AU,BU,CU,DU);
  520. else
  521.    acp = ak; bcp = bk1; ccp = ck1; dcp = dk11;
  522. end
  523. %
  524. % ------------------------------------ Closed-loop TF (Ty1u1):
  525. %
  526. [acl,bcl,ccl,dcl] = lftf(sysp,dimp,acp,bcp,ccp,dcp);
  527. lamcl = eig(acl);
  528. RHP_cl = size(find(real(lamcl)>0))*[1;0];
  529. if RHP_cl > 0
  530.    if disp_hinf == 1
  531.       disp(''),          disp('               -- CLOSED-LOOP UNSTABLE -- ')
  532.    end
  533.    if hexist == ' OK '
  534.       if disp_hinf == 1
  535.                       disp('               CANNOT COMPUTE CONTROLLER')
  536.                       disp('            DUE TO NUMERICAL ROUNDOFF ERRORS')
  537.       end
  538.    end
  539. else 
  540.    if disp_hinf == 1
  541.       disp('                         DONE!!!')
  542.       disp('    -------------------------------------------------------');
  543.    end
  544. end
  545. syscp = [acp bcp;ccp dcp]; xcp = size(acp)*[1;0];
  546. syscl = [acl bcl;ccl dcl]; xcl = size(acl)*[1;0];
  547. hinflag = [hexist D11flag Pexist Pflag Sexist Sflag PSflag];
  548. hinfo = [abs(hinflag),RHP_cl,lamps_max];
  549. if nag2 == 0
  550.    % Case when used as script file with no output arguments
  551.    ss_cl=mksys(acl,bcl,ccl,dcl);
  552.    ss_cp=mksys(acp,bcp,ccp,dcp);
  553.    save hinf2 acp bcp ccp dcp acl bcl ccl dcl hinfo ss_cl ss_cp
  554.    hinfload
  555.    acp = 'Controller: (acp,bcp,ccp,dcp); Cost: (acl,bcl,ccl,dcl); Test Info: hinfo';   
  556.    delete hinf1.mat
  557.    delete hinf2.mat
  558. else
  559.    if xsflag
  560.      % Case     [SS_CP,SS_CL,HINFO,TSS_K]=HINF(TSS_P,SS_U,..)
  561.      acp=mksys(acp,bcp,ccp,dcp);
  562.      if nag2>1,
  563.            bcp=mksys(acl,bcl,ccl,dcl);
  564.            ccp=hinfo;
  565.            if nag2>3,
  566.              dcp=mksys(ak,bk1,bk2,ck1,ck2,dk11,dk12,dk21,dk22,'tss');
  567.            end
  568.      end
  569.    end
  570. end
  571. % If none of the above, then output args remain:
  572. %     [acp,bcp,ccp,dcp,acp,bcl,ccl,dcl,hinfo,ak,bk1,bk2,ck1,ck2,dk11,dk12,...
  573. %          dk21,dk22]=hinf(A,B1,B2,C1,C2,D11,D12,D21,D22,AU,BU,CU,DU,VERBOSE)
  574. %
  575. % ---------- End of HINF.M --- RYC/MGS -- %
  576.