home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts / control / hinfsyn.m < prev    next >
Encoding:
Text File  |  1999-12-15  |  9.7 KB  |  314 lines
  1. ## Copyright (C) 1996,1998 Auburn University.  All Rights Reserved
  2. ##
  3. ## This file is part of Octave. 
  4. ##
  5. ## Octave is free software; you can redistribute it and/or modify it 
  6. ## under the terms of the GNU General Public License as published by the 
  7. ## Free Software Foundation; either version 2, or (at your option) any 
  8. ## later version. 
  9. ## 
  10. ## Octave is distributed in the hope that it will be useful, but WITHOUT 
  11. ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
  12. ## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
  13. ## for more details.
  14. ## 
  15. ## You should have received a copy of the GNU General Public License 
  16. ## along with Octave; see the file COPYING.  If not, write to the Free 
  17. ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. 
  18.  
  19. ## -*- texinfo -*-
  20. ## @deftypefn {Function File } {[@var{K}, @var{g}, @var{GW}, @var{Xinf}, @var{Yinf}] =} hinfsyn(@var{Asys}, @var{nu}, @var{ny}, @var{gmin}, @var{gmax}, @var{gtol}@{, @var{ptol}, @var{tol}@})
  21. ## 
  22. ## @strong{Inputs} input system is passed as either
  23. ## @table @var
  24. ## @item Asys
  25. ## system data structure (see ss2sys, sys2ss)
  26. ## @itemize @bullet
  27. ## @item controller is implemented for continuous time systems 
  28. ## @item controller is NOT implemented for discrete time systems  (see
  29. ## bilinear transforms in @code{c2d}, @code{d2c})
  30. ## @end itemize
  31. ## @item nu
  32. ## number of controlled inputs
  33. ## @item ny
  34. ## number of measured outputs
  35. ## @item gmin
  36. ## initial lower bound on H-infinity optimal gain
  37. ## @item gmax
  38. ## initial upper bound on H-infinity optimal gain
  39. ## @item gtol
  40. ## gain threshhold.  Routine quits when gmax/gmin < 1+tol
  41. ## @item ptol
  42. ## poles with abs(real(pole)) < ptol*||H|| (H is appropriate
  43. ## Hamiltonian) are considered to be on the imaginary axis.  
  44. ## Default: 1e-9
  45. ## @item tol
  46. ## threshhold for 0.  Default: 200*eps
  47. ## 
  48. ## @var{gmax}, @var{min}, @var{tol}, and @var{tol} must all be postive scalars.
  49. ## @end table 
  50. ## @strong{Outputs}
  51. ## @table @var
  52. ## @item K
  53. ## system controller
  54. ## @item g
  55. ## designed gain value
  56. ## @item GW
  57. ## closed loop system
  58. ## @item Xinf
  59. ## ARE solution matrix for regulator subproblem
  60. ## @item Yinf
  61. ## ARE solution matrix for filter subproblem
  62. ## @end table
  63. ## 
  64. ## @enumerate
  65. ## @item Doyle, Glover, Khargonekar, Francis, "State Space Solutions
  66. ##      to Standard H2 and Hinf Control Problems," IEEE TAC August 1989
  67. ## 
  68. ## @item Maciejowksi, J.M., "Multivariable feedback design,"
  69. ##      Addison-Wesley, 1989, ISBN 0-201-18243-2
  70. ## 
  71. ## @item Keith Glover and John C. Doyle, "State-space formulae for all
  72. ##      stabilizing controllers that satisfy and h-infinity-norm bound
  73. ##      and relations to risk sensitivity,"
  74. ##      Systems & Control Letters 11, Oct. 1988, pp 167-172.
  75. ## @end enumerate
  76. ## @end deftypefn
  77.  
  78. function [K, g, GW, Xinf, Yinf] = hinfsyn (Asys, nu, ny, gmin, gmax, gtol, ptol, tol)
  79.  
  80.   ## A. S. Hodel August 1995
  81.   ## Updated for Packed system structures December 1996 by John Ingram
  82.   ## 
  83.   ## Revised by Kai P Mueller April 1998 to solve the general H_infinity
  84.   ## problem using unitary transformations Q (on w and z)
  85.   ## and non-singular transformations R (on u and y).
  86.  
  87.   if( (nargin < 1) | (nargin > 8) )
  88.     usage("[K,g,GW,Xinf,Yinf] = hinfsyn(Asys,nu,ny,gmin,gmax,gtol,ptol,tol)");
  89.   endif
  90.   ## set default arguments
  91.   if(nargin < 8)
  92.     tol = 200*eps;
  93.   elseif(!is_sample(tol))
  94.     error("tol must be a positive scalar.")
  95.   endif
  96.   if(nargin < 7)
  97.     ptol = 1e-9;
  98.   elseif(!is_sample(ptol))
  99.     error("hinfsyn: ptol must be a positive scalar");
  100.   endif
  101.     
  102.   if(!is_sample(gmax) | !is_sample(gmin) | !is_sample(gtol) )
  103.     error(["hinfsyn: gmax=",num2str(gmax),", gmin=",num2str(gmin), ...
  104.       "gtol=",num2str(gtol), " must be positive scalars."])
  105.   endif
  106.  
  107.   [chkdgkf,dgs] = is_dgkf(Asys,nu,ny,tol);
  108.  
  109.   if (! chkdgkf )
  110.     disp("hinfsyn: system does not meet required assumptions")
  111.     help is_dgkf
  112.     error("hinfsyn: exit");
  113.   endif
  114.  
  115.   ## extract dgs information
  116.               nw = dgs.nw;    nu = dgs.nu;
  117.   A = dgs.A;        B1 = dgs.Bw;    B2 = dgs.Bu;
  118.   C1 = dgs.Cz;        D11 = dgs.Dzw;    D12 = dgs.Dzu;        nz = dgs.nz;
  119.   C2 = dgs.Cy;        D21 = dgs.Dyw;    D22 = dgs.Dyu;        ny = dgs.ny;
  120.   d22nz = dgs.Dyu_nz;
  121.   dflg = dgs.dflg;
  122.  
  123.   ## recover i/o transformations
  124.   R12 = dgs.Ru;        R21 = dgs.Ry;
  125.   [ncstates, ndstates, nin, nout] = sysdimensions(Asys);
  126.   Atsam = sysgettsam(Asys);
  127.   [Ast, Ain, Aout] = sysgetsignals(Asys);
  128.  
  129.   BB = [B1, B2];
  130.   CC = [C1 ; C2];
  131.   DD = [D11, D12 ; D21,  D22];
  132.  
  133.   if (dflg == 0)
  134.     n = ncstates;
  135.     ## perform binary search to find gamma min
  136.     ghi = gmax;
  137.     ## derive a lower lower bound for gamma from D matrices
  138.     xx1 = norm((eye(nz) - (D12/(D12'*D12))*D12')*D11);
  139.     xx2 = norm(D11*(eye(nw)-(D21'/(D21*D21'))*D21));
  140.     glo = max(xx1, xx2);
  141.     if (glo > gmin)
  142.       disp(" *** D matrices indicate a greater value of gamma min.");
  143.       fprintf("     gamma min (%f) superseeded by %f.", gmin, glo);
  144.       glo = xx1;
  145.     else
  146.       glo = gmin;
  147.     endif
  148.     if (glo > ghi)
  149.       fprintf(" *** lower bound of gamma greater than upper bound(%f)", ...
  150.           glo, ghi);
  151.       disp(" *** unable to continue, Goodbye.");
  152.       return;
  153.     endif
  154.  
  155.     de = ghi - glo;
  156.     g = glo;
  157.     search_state = 0;
  158.     iteration_finished = 0;
  159.     disp(" o structural tests passed, start of iteration...");
  160.     disp("        o <-> test passed   # <-> test failed   - <-> cannot test");
  161.     printf("----------------------------------------");
  162.     printf("--------------------------------------\n");
  163.  
  164.     ## ------123456789012345678901234567890123456789012345678901234567890
  165.     printf("           .........X......... .........Y......... ");
  166.     printf(".Z. PASS REMARKS\n");
  167.     printf("        ga iax nev ene sym pos iax nev ene sym pos ");
  168.     printf("rho  y/n ======>\n");
  169.     printf("----------------------------------------");
  170.     printf("--------------------------------------\n");
  171.  
  172.     ## set up error messages
  173.     errmesg = list(" o   o   o   o   o  ", ...
  174.     " #   -   -   -   -  ", ...
  175.     " o   #   -   -   -  ", ...
  176.     " o   o   #   -   -  ", ...
  177.     " o   o   o   #   -  ", ...
  178.     " o   o   o   o   #  ", ...
  179.     " -   -   -   -   -  ");
  180.     errdesx = list("", ...
  181.     "X im eig.", ...
  182.     "Hx not Ham.", ...
  183.     "X inf.eig", ...
  184.     "X not symm.", ...
  185.     "X not pos", ...
  186.     "R singular");
  187.  
  188.     errdesy = list(" ", ...
  189.     "Y im eig.", ...
  190.     "Hy not Ham.", ...
  191.     "Y inf.eig", ...
  192.     "Y not symm.", ...
  193.     "Y not pos", ...
  194.     "Rtilde singular");
  195.  
  196.  
  197.     ## now do the search
  198.     while (!iteration_finished)
  199.       switch (search_state)
  200.         case (0)     g = ghi;
  201.         case (1)     g = glo;
  202.         case (2)     g = 0.5 * (ghi + glo);
  203.         otherwise     error(" *** This should never happen!");
  204.       endswitch
  205.       printf("%10.4f ", g);
  206.  
  207.       ## computing R and R~
  208.       d1dot = [D11, D12];
  209.       R = zeros(nin, nin);
  210.       R(1:nw,1:nw) = -g*g*eye(nw);
  211.       R = R + d1dot' * d1dot;
  212.       ddot1 = [D11; D21];
  213.       Rtilde = zeros(nout, nout);
  214.       Rtilde(1:nz,1:nz) = -g*g*eye(nz);
  215.       Rtilde = Rtilde + ddot1 * ddot1';
  216.  
  217.       [Xinf,x_ha_err] = hinfsyn_ric(A,BB,C1,d1dot,R,ptol);
  218.       [Yinf,y_ha_err] = hinfsyn_ric(A',CC',B1',ddot1',Rtilde,ptol);
  219.  
  220.       ## assume failure for this gamma
  221.       passed = 0;
  222.       rerr="";
  223.       if (!x_ha_err && !y_ha_err)
  224.     ## test spectral radius condition
  225.     rho = max(abs(eig(Xinf * Yinf)));
  226.     if (rho < g*g)
  227.       ## spectral radius condition passed
  228.       passed = 1;
  229.         else
  230.           rerr = sprintf("rho=%f",rho);
  231.     endif
  232.       endif
  233.  
  234.       if(x_ha_err >= 0 & x_ha_err <= 6)
  235.         printf("%s",nth(errmesg,x_ha_err+1));
  236.         xerr = nth(errdesx,x_ha_err+1);
  237.       else
  238.         error(" *** Xinf fail: this should never happen!");
  239.       endif
  240.  
  241.       if(y_ha_err >= 0 & y_ha_err <= 6)
  242.         printf("%s",nth(errmesg,y_ha_err+1));
  243.         yerr = nth(errdesy,y_ha_err+1);
  244.       else
  245.         error(" *** Yinf fail: this should never happen!");
  246.       endif
  247.  
  248.       if(passed)  printf("  y all tests passed.\n");
  249.       else        printf("  n %s/%s%s\n",xerr,yerr,rerr);          endif
  250.  
  251.       if (passed && (de/g < gtol))
  252.     search_state = 3;
  253.       endif
  254.  
  255.       switch (search_state)
  256.         case (0)
  257.       if (!passed)
  258.         ## upper bound must pass but did not
  259.         fprintf(" *** the upper bound of gamma (%f) is too small.\n", g);
  260.         iteration_finished = 2;
  261.       else
  262.             search_state = 1;
  263.       endif
  264.         case (1)
  265.       if (!passed)      search_state = 2;
  266.       else
  267.         ## lower bound must not pass but passed
  268.         fprintf(" *** the lower bound of gamma (%f) passed.\n", g);
  269.         iteration_finished = 3;
  270.       endif
  271.         case (2)
  272.           ## Normal case; must check that singular R, Rtilde wasn't the problem.
  273.       if ((!passed) & (x_ha_err != 6) & (y_ha_err != 6) ) glo = g;
  274.       else         ghi = g;        endif
  275.       de = ghi - glo;
  276.         case (3)       iteration_finished = 1;        # done
  277.         otherwise      error(" *** This should never happen!");
  278.       endswitch
  279.     endwhile
  280.  
  281.     printf("----------------------------------------");
  282.     printf("--------------------------------------\n");
  283.     if (iteration_finished != 1)
  284.       K = [];
  285.     else
  286.       ## success: compute controller
  287.       fprintf("   hinfsyn final: glo=%f ghi=%f, test gain g=%f\n", \
  288.               glo, ghi, g);
  289.       printf("----------------------------------------");
  290.       printf("--------------------------------------\n");
  291.       Z = inv(eye(ncstates) - Yinf*Xinf/g/g);
  292.       F = -R \ (d1dot'*C1 + BB'*Xinf);
  293.       H = -(B1*ddot1' + Yinf*CC') / Rtilde;
  294.       K = hinf_ctr(dgs,F,H,Z,g);
  295.  
  296.       Kst = strappend(Ast,"_K");
  297.       Kin = strappend(Aout((nout-ny+1):(nout)),"_K");
  298.       Kout = strappend(Ain((nin-nu+1):(nin)),"_K");
  299.       [Ac, Bc, Cc, Dc] = sys2ss(K);
  300.       K = ss2sys(Ac,Bc,Cc,Dc,Atsam,ncstates,ndstates,Kst,Kin,Kout);
  301.       if (nargout >= 3)
  302.     GW = starp(Asys, K);
  303.       endif
  304.     endif
  305.     
  306.   elseif(ndstates)
  307.  
  308.     ## discrete time solution
  309.     error("hinfsyn: discrete-time case not yet implemented")
  310.  
  311.   endif
  312.  
  313. endfunction
  314.