|
3381
|
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. |
|
3346
|
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 |
|
3213
|
65 |
|
3385
|
66 function [g, gmin, gmax] = hinfnorm (sys, tol, gmin, gmax, ptol) |
|
3213
|
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 |
|
3381
|
74 ## set defaults where applicable |
|
3213
|
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 |
|
3381
|
87 ## eigenvalues of A must all be stable |
|
3213
|
88 if(!is_stable(sys)) |
|
|
89 warning(["hinfnorm: unstable system (is_stable, ptol=",num2str(ptol), ... |
|
|
90 "), returning Inf"]); |
|
3228
|
91 g = Inf; |
|
3213
|
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 |
|
3381
|
105 ## check for scalar input arguments 2...5 |
|
3213
|
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); |
|
3381
|
114 ## find the Hinf norm via binary search |
|
3213
|
115 while((gmax/gmin - 1) > tol) |
|
|
116 g = (gmax+gmin)/2; |
|
|
117 |
|
|
118 if(dflg) |
|
3381
|
119 ## multiply g's through in formulas to avoid extreme magnitudes... |
|
3213
|
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 |
|
3381
|
124 ## set up symplectic generalized eigenvalue problem per Iglesias & Glover |
|
3213
|
125 s1 = [Ak , zeros(nz) ; -Ck, In ]; |
|
|
126 s2 = [In, -(B/Rg)*B' ; zeros(nz) , Ak' ]; |
|
|
127 |
|
3381
|
128 ## guard against roundoff again: zero out extremely small values |
|
|
129 ## prior to balancing |
|
3213
|
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]); |
|
3238
|
136 Hb = [s1, s2]; |
|
3213
|
137 |
|
3381
|
138 ## check R - B' X B condition (Iglesias and Glover's paper) |
|
3213
|
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)']; |
|
3381
|
145 ## guard against roundoff: zero out extremely small values prior |
|
|
146 ## to balancing |
|
3213
|
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 |
