home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / hinf_ctr.m < prev    next >
Text File  |  1999-04-29  |  4KB  |  136 lines
  1. # Copyright (C) 1996 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 = hinf_ctr(dgs,F,H,Z,g)
  20.   # K = hinf_ctr(dgs,F,H,Z,g)
  21.   #
  22.   # Called by hinfsyn to compute the H_inf optimal controller.
  23.   # 
  24.   # inputs:
  25.   #   dgs: data structure returned by is_dgkf
  26.   #           F, H:         feedback and filter gain (not partitioned)
  27.   #           g:            final gamma value
  28.   # outputs: 
  29.   #           controller K (system data structure)
  30.   #
  31.   # Do not attempt to use this at home; no argument checking performed.
  32.  
  33.   # A. S. Hodel August 1995
  34.   # Revised by Kai P Mueller April 1998 to solve the general H_infinity
  35.   # problem using unitary transformations Q (on w and z)
  36.   # and non-singular transformations R (on u and y).
  37.  
  38.   nw = dgs.nw;
  39.   nu = dgs.nu;
  40.   nz = dgs.nz;
  41.   ny = dgs.ny;
  42.   d22nz = dgs.Dyu_nz;
  43.   
  44.   B1  = dgs.Bw;
  45.   B2  = dgs.Bu;
  46.   C1  = dgs.Cz;
  47.   C2  = dgs.Cy;
  48.   C = [C1; C2];
  49.   D11 = dgs.Dzw;
  50.   D12 = dgs.Dzu;
  51.   D21 = dgs.Dyw;
  52.   D22 = dgs.Dyu;
  53.   A = dgs.A;
  54.   Ru = dgs.Ru;
  55.   Ry = dgs.Ry;
  56.  
  57.   nout = nz + ny;
  58.   nin = nw + nu;
  59.   nstates = size(A, 1);
  60.  
  61.   F11 = F(1:(nw-ny),:);
  62.   F12 = F((nw-ny+1):nw,:);
  63.   F2  = F((nw+1):nin,:);
  64.   H11 = H(:,1:(nz-nu));
  65.   H12 = H(:,(nz-nu+1):nz);
  66.   H2  = H(:,(nz+1):nout);
  67.  
  68.   # D11 partitions
  69.   D1111 = D11(1:(nz-nu),1:(nw-ny));
  70.   D1112 = D11(1:(nz-nu),(nw-ny+1):nw);
  71.   D1121 = D11((nz-nu+1):nz,1:(nw-ny));
  72.   D1122 = D11((nz-nu+1):nz,(nw-ny+1):nw);
  73.  
  74.   # D11ik may be the empty matrix, don't calculate with empty matrices
  75.   [nd1111,md1111] = size(D1111);
  76.   md1112 = length(D1112);
  77.   md1121 = length(D1121);
  78.  
  79.   if ((nd1111 == 0) || (md1112 == 0))
  80.     d11hat = -D1122;
  81.   else
  82.     xx = inv(g*g*eye(nz-nu) - D1111*D1111');
  83.     d11hat = -D1121*D1111'*xx*D1112 - D1122;
  84.   endif
  85.   if (md1112 == 0)
  86.     d21hat = eye(ny);
  87.   elseif (nd1111 == 0)
  88.     d21hat = chol(eye(ny) - D1112'*D1112/g/g);
  89.   else
  90.     xx = inv(g*g*eye(nz-nu) - D1111*D1111');
  91.     xx = eye(ny) - D1112'*xx*D1112;
  92.     d21hat = chol(xx);
  93.   endif
  94.   if (md1121 == 0)
  95.     d12hat = eye(nu);
  96.   elseif (md1111 == 0)
  97.     d12hat = chol(eye(nu) - D1121*D1121'/g/g)';
  98.   else
  99.     xx = inv(g*g*eye(nw-ny) - D1111'*D1111);
  100.     xx = eye(nu)-D1121*xx*D1121';
  101.     d12hat = chol(xx)';
  102.   endif
  103.  
  104.   b2hat = (B2+H12)*d12hat;
  105.   c2hat = -d21hat*(C2+F12)*Z;
  106.   b1hat = -H2 + (b2hat/d12hat)*d11hat;
  107.   c1hat =  F2*Z + (d11hat/d21hat)*c2hat;
  108.   ahat  =  A + H*C + (b2hat/d12hat)*c1hat;
  109.  
  110.   # rescale controller by Ru and Ry
  111.   b1hat = b1hat/Ry;
  112.   c1hat = Ru\c1hat;
  113.   bhat  = [b1hat, b2hat];
  114.   chat  = [c1hat; c2hat];
  115.   dhat  = [Ru\d11hat/Ry, Ru\d12hat; d21hat/Ry, 0*d11hat'];
  116.  
  117.   # non-zero D22 is a special case
  118.   if (d22nz)
  119.     if (rank(eye(nu) + d11hat*D22) < nu)
  120.       error(" *** cannot compute controller for D22 non-zero.");
  121.     endif
  122.  
  123.     d22new = [D22, zeros(ny,ny); zeros(nu,nu), 0*D22'];
  124.     xx = inv(eye(nu+ny) + d22new*dhat);
  125.     mhat = inv(eye(nu+ny) + dhat*d22new);
  126.     ahat = ahat - bhat*((eye(nu+ny)-xx)/dhat)*chat;
  127.     bhat = bhat*xx;
  128.     chat = mhat*chat;
  129.     dhat = dhat*xx;
  130.  
  131.   endif
  132.  
  133.   K = ss2sys(ahat,bhat(:,1:ny),chat(1:nu,:),dhat(1:nu,1:ny));
  134.  
  135. endfunction
  136.