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