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