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