home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts / control / hinfnorm.m < prev    next >
Encoding:
Text File  |  1999-12-15  |  5.1 KB  |  161 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{g}, @var{gmin}, @var{gmax}] =} hinfnorm(@var{sys}@{, @var{tol}, @var{gmin}, @var{gmax}, @var{ptol}@})
  21. ##  Computes the H infinity norm of a system data structure.
  22. ## 
  23. ## @strong{Inputs}
  24. ## @table @var
  25. ## @item sys 
  26. ## system data structure
  27. ## @item tol 
  28. ## H infinity norm search tolerance (default: 0.001)
  29. ## @item gmin 
  30. ## minimum value for norm search (default: 1e-9)
  31. ## @item gmax 
  32. ## maximum value for norm search (default: 1e+9)
  33. ## @item ptol
  34. ##  pole tolerance:
  35. ## @itemize @bullet
  36. ## @item if sys is continuous, poles with 
  37. ## |real(pole)| < ptol*||H|| (H is appropriate Hamiltonian)
  38. ## are considered to be on the imaginary axis.  
  39. ## 
  40. ## @item if sys is discrete, poles with
  41. ## |abs(pole)-1| < ptol*||[s1,s2]|| (appropriate symplectic pencil)
  42. ## are considered to be on the unit circle
  43. ## 
  44. ## @item Default: 1e-9
  45. ## @end itemize
  46. ## @end table
  47. ## 
  48. ## @strong{Outputs}
  49. ## @table @var
  50. ## @item g
  51. ## Computed gain, within @var{tol} of actual gain.  @var{g} is returned as Inf 
  52. ## if the system is unstable.
  53. ## @item gmin, gmax
  54. ## Actual system gain lies in the interval [@var{gmin}, @var{gmax}]
  55. ## @end table
  56. ## 
  57. ##  References:
  58. ##  Doyle, Glover, Khargonekar, Francis, "State space solutions to standard
  59. ##     H2 and Hinf control problems", IEEE TAC August 1989
  60. ##  Iglesias and Glover, "State-Space approach to discrete-time Hinf control,"
  61. ##     Int. J. Control, vol 54, #5, 1991
  62. ##  Zhou, Doyle, Glover, "Robust and Optimal Control," Prentice-Hall, 1996
  63. ##  $Revision: 1.9 $
  64. ## @end deftypefn
  65.  
  66. function [g, gmin, gmax] = hinfnorm (sys, tol, gmin, gmax, ptol)
  67.  
  68.   if((nargin == 0) || (nargin > 4))
  69.     usage("[g gmin gmax] = hinfnorm(sys[,tol,gmin,gmax,ptol])");
  70.   elseif(!is_struct(sys))
  71.     error("Sys must be a system data structure");
  72.   endif
  73.  
  74.   ## set defaults where applicable
  75.   if(nargin < 5)
  76.     ptol = 1e-9;    # pole tolerance
  77.   endif
  78.   if(nargin < 4)
  79.     gmax = 1e9;        # max gain value
  80.   endif
  81.  
  82.   dflg = is_digital(sys);
  83.   sys = sysupdate(sys,"ss");
  84.   [A,B,C,D] = sys2ss(sys);
  85.   [n,nz,m,p] = sysdimensions(sys);
  86.  
  87.   ## eigenvalues of A must all be stable
  88.   if(!is_stable(sys))
  89.     warning(["hinfnorm: unstable system (is_stable, ptol=",num2str(ptol), ...
  90.       "), returning Inf"]);
  91.     g = Inf;
  92.   endif
  93.  
  94.   Dnrm = norm(D);
  95.   if(nargin < 3)
  96.     gmin = max(1e-9,Dnrm);     # min gain value
  97.   elseif(gmin < Dnrm)
  98.     warning(["hinfnorm: setting Gmin=||D||=",num2str(Dnrm)]);
  99.   endif
  100.  
  101.   if(nargin < 2)
  102.     tol = 0.001;    # convergence measure for gmin, gmax
  103.   endif
  104.  
  105.   ## check for scalar input arguments 2...5
  106.   if( ! (is_scalar(tol) && is_scalar(gmin) 
  107.     && is_scalar(gmax) && is_scalar(ptol)) )
  108.     error("hinfnorm: tol, gmin, gmax, ptol must be scalars");
  109.   endif
  110.  
  111.   In = eye(n+nz);
  112.   Im = eye(m);
  113.   Ip = eye(p);
  114.   ## find the Hinf norm via binary search
  115.   while((gmax/gmin - 1) > tol)
  116.     g = (gmax+gmin)/2;
  117.  
  118.     if(dflg)
  119.       ## multiply g's through in formulas to avoid extreme magnitudes...
  120.       Rg = g^2*Im - D'*D;
  121.       Ak = A + (B/Rg)*D'*C;
  122.       Ck = g^2*C'*((g^2*Ip-D*D')\C);
  123.  
  124.       ## set up symplectic generalized eigenvalue problem per Iglesias & Glover
  125.       s1 = [Ak , zeros(nz) ; -Ck, In ];
  126.       s2 = [In, -(B/Rg)*B' ; zeros(nz) , Ak' ];
  127.  
  128.       ## guard against roundoff again: zero out extremely small values
  129.       ## prior to balancing
  130.       s1 = s1 .* (abs(s1) > ptol*norm(s1,"inf"));
  131.       s2 = s2 .* (abs(s2) > ptol*norm(s2,"inf"));
  132.       [cc,dd,s1,s2] = balance(s1,s2);
  133.       [qza,qzb,zz,pls] = qz(s1,s2,"S");    # ordered qz decomposition
  134.       eigerr = abs(abs(pls)-1);
  135.       normH = norm([s1,s2]);
  136.       Hb = [s1, s2];
  137.  
  138.       ## check R - B' X B condition (Iglesias and Glover's paper)
  139.       X = zz((nz+1):(2*nz),1:nz)/zz(1:nz,1:nz);
  140.       dcondfailed = min(real( eig(Rg - B'*X*B)) < ptol);
  141.     else
  142.       Rinv = inv(g*g*Im - (D' * D));
  143.       H = [A + B*Rinv*D'*C,        B*Rinv*B'; ...
  144.            -C'*(Ip + D*Rinv*D')*C, -(A + B*Rinv*D'*C)'];
  145.       ## guard against roundoff: zero out extremely small values prior 
  146.       ## to balancing
  147.       H = H .* (abs(H) > ptol*norm(H,"inf"));
  148.       [DD,Hb] = balance(H);
  149.       pls = eig(Hb);
  150.       eigerr = abs(real(pls));
  151.       normH = norm(H);
  152.       dcondfailed = 0;        # digital condition; doesn't apply here
  153.     endif
  154.     if( (min(eigerr) <= ptol * normH) | dcondfailed)
  155.       gmin = g;
  156.     else
  157.       gmax = g;
  158.     endif
  159.   endwhile
  160. endfunction
  161.