|
3228
|
1 # Copyright (C) 1996,1998 A. Scottedward Hodel |
|
3213
|
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"]); |
|
3228
|
66 g = Inf; |
|
3213
|
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 |
