home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / h2syn.m < prev    next >
Text File  |  1999-04-29  |  5KB  |  151 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,gain, Kc, Kf, Pc,  Pf] = h2syn(Asys,nu,ny,tol)
  20.   #  [K,gain, Kc, Kf, Pc, Pf] = h2syn(Asys,nu,ny,tol)
  21.   #
  22.   # Design H2 optimal controller per procedure in 
  23.   # Doyle, Glover, Khargonekar, Francis, "State Space Solutions to Standard
  24.   # H2 and Hinf Control Problems", IEEE TAC August 1989
  25.   #
  26.   # Discrete time control per Zhou, Doyle, and Glover, ROBUST AND OPTIMAL
  27.   #   CONTROL, Prentice-Hall, 1996
  28.   #
  29.   # inputs: 
  30.   #   Asys: system data structure.
  31.   #   nu: number of controlled inputs
  32.   #   ny: number of measured outputs
  33.   #   tol: threshhold for 0.  Default: 200*eps
  34.   # outputs: 
  35.   #    K: system controller (system data structure)
  36.   #    gain: optimal closed loop gain (see Kc, Kf warning below)
  37.   #    Kc: full information control (system data structure)
  38.   #    Kf: state estimator (system data structure)
  39.   #       WARNING: incorporation of the is_dgkf nonsingular transformations
  40.   #       Ru and Ry into Kc and Kf has not been tested.  
  41.   #    Pc: ARE solution matrix for regulator subproblem
  42.   #    Pf: ARE solution matrix for filter subproblem
  43.  
  44.   # Updated for System structure December 1996 by John Ingram
  45.  
  46.   if ((nargin < 3) | (nargin > 4))
  47.     usage("[K,gain, Kc, Kf, Pc, Pf] = h2syn(Asys,nu,ny[,tol])");
  48.   elseif(nargin == 3 )
  49.     [chkdgkf,dgs] = is_dgkf(Asys,nu,ny);
  50.   elseif(nargin == 4)
  51.     [chkdgkf,dgs] = is_dgkf(Asys,nu,ny,tol);
  52.   endif
  53.  
  54.   if (!chkdgkf )
  55.     disp("h2syn: system does not meet required assumptions")
  56.     help is_dgkf
  57.     error("h2syn: exit");
  58.   endif
  59.  
  60.   # extract dgs information
  61.               nw = dgs.nw;     nu = dgs.nu;
  62.   A = dgs.A;            Bw = dgs.Bw;    Bu = dgs.Bu;
  63.   Cz = dgs.Cz;          Dzw = dgs.Dzw;  Dzu = dgs.Dzu;    nz = dgs.nz;
  64.   Cy = dgs.Cy;          Dyw = dgs.Dyw;  Dyu = dgs.Dyu;    ny = dgs.ny;
  65.   d22nz = dgs.Dyu_nz;
  66.   dflg = dgs.dflg;
  67.  
  68.   if(norm(Dzw,Inf) > norm([Dzw, Dzu ; Dyw, Dyu],Inf)*1e-12)
  69.     warning("h2syn: Dzw nonzero; feedforward not implemented")
  70.     Dzw
  71.     D = [Dzw, Dzu ; Dyw, Dyu]
  72.   endif
  73.  
  74.   # recover i/o transformations
  75.   Ru = dgs.Ru;         Ry = dgs.Ry;
  76.   [ncstates, ndstates, nout, nin] = sysdimen(Asys);
  77.   Atsam = sysgetts(Asys);
  78.   [Ast, Ain, Aout] = sysgetsg(Asys);
  79.  
  80.   if(dgs.dflg == 0)
  81.     Pc = are(A,Bu*Bu',Cz'*Cz);    # solve control, filtering ARE's
  82.     Pf = are(A',Cy'*Cy,Bw*Bw');
  83.     F2 = -Bu'*Pc;          # calculate feedback gains
  84.     L2 = -Pf*Cy';
  85.  
  86.     AF2 = A + Bu*F2;
  87.     AL2 = A + L2*Cy;
  88.     CzF2 = Cz + (Dzu/Ru)*F2;
  89.     BwL2 = Bw+L2*(Ry\Dyw);
  90.  
  91.   else
  92.     # discrete time solution
  93.     error("h2syn: discrete-time case not yet implemented")
  94.     Pc = dare(A,Bu*Bu',Cz'*Cz);
  95.     Pf = dare(A',Cy'*Cy,Bw*Bw');
  96.   endif
  97.  
  98.   nn = ncstates + ndstates;
  99.   In = eye(nn);
  100.   KA = A + Bu*F2 + L2*Cy;
  101.   Kc1 = ss2sys(AF2,Bw,CzF2,zeros(nz,nw));
  102.   Kf1 = ss2sys(AL2,BwL2,F2,zeros(nu,nw));
  103.  
  104.   g1 = h2norm(Kc1);
  105.   g2 = h2norm(Kf1);
  106.   
  107.   # compute optimal closed loop gain
  108.   gain = sqrt ( g1*g1 + g2*g2 );
  109.  
  110.   if(nargout)
  111.     Kst = strappen(Ast,"_K");
  112.     Kin = strappen(Aout((nout-ny+1):(nout)),"_K");
  113.     Kout = strappen(Ain((nin-nu+1):(nin)),"_K");
  114.  
  115.     # compute systems for return
  116.     K = ss2sys(KA,-L2/Ru,Ry\F2,zeros(nu,ny),Atsam,ncstates,ndstates,Kst,Kin,Kout);
  117.   endif
  118.  
  119.   if (nargout > 2)
  120.     #system full information control state names
  121.     stname2 = strappen(Ast,"_FI");
  122.  
  123.    #system full information control input names
  124.    inname2 = strappen(Ast,"_FI_in");
  125.  
  126.     #system full information control output names
  127.     outname2 = strappen(Aout(1:(nout-ny)),"_FI_out");
  128.  
  129.     nz = rows (Cz);
  130.     nw = columns (Bw);
  131.  
  132.     Kc = ss2sys(AF2, In, CzF2, zeros(nz,nn), Atsam, ...
  133.     ncstates, ndstates, stname2, inname2, outname2);
  134.   endif
  135.  
  136.   if (nargout >3)
  137.     #fix system state estimator state names
  138.     stname3 = strappen(Ast,"_Kf");
  139.  
  140.     #fix system state estimator input names
  141.     inname3 = strappen(Ast,"_Kf_noise");
  142.  
  143.     #fix system state estimator output names
  144.     outname3 = strappen(Ast,"_est");
  145.  
  146.     Kf = ss2sys(AL2, BwL2, In, zeros(nn,nw),Atsam,  ...
  147.       ncstates, ndstates, stname3, inname3,outname3);
  148.   endif
  149.  
  150. endfunction
  151.