|
3430
|
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 |
|
7016
|
6 ## under the terms of the GNU General Public License as published by |
|
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
8 ## your option) any later version. |
|
3430
|
9 ## |
|
7016
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
3430
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
|
17 ## <http://www.gnu.org/licenses/>. |
|
3430
|
18 |
|
|
19 ## -*- texinfo -*- |
|
3501
|
20 ## @deftypefn {Function File} {[@var{g}, @var{gmin}, @var{gmax}] =} hinfnorm (@var{sys}, @var{tol}, @var{gmin}, @var{gmax}, @var{ptol}) |
|
5016
|
21 ## Computes the |
|
|
22 ## @iftex |
|
|
23 ## @tex |
|
|
24 ## $ { \cal H }_\infty $ |
|
|
25 ## @end tex |
|
|
26 ## @end iftex |
|
|
27 ## @ifinfo |
|
|
28 ## H-infinity |
|
|
29 ## @end ifinfo |
|
|
30 ## norm of a system data structure. |
|
3430
|
31 ## |
|
|
32 ## @strong{Inputs} |
|
|
33 ## @table @var |
|
|
34 ## @item sys |
|
|
35 ## system data structure |
|
|
36 ## @item tol |
|
5016
|
37 ## @iftex |
|
|
38 ## @tex |
|
|
39 ## $ { \cal H }_\infty $ |
|
|
40 ## @end tex |
|
|
41 ## @end iftex |
|
|
42 ## @ifinfo |
|
|
43 ## H-infinity |
|
|
44 ## @end ifinfo |
|
|
45 ## norm search tolerance (default: 0.001) |
|
3430
|
46 ## @item gmin |
|
|
47 ## minimum value for norm search (default: 1e-9) |
|
|
48 ## @item gmax |
|
|
49 ## maximum value for norm search (default: 1e+9) |
|
|
50 ## @item ptol |
|
|
51 ## pole tolerance: |
|
|
52 ## @itemize @bullet |
|
|
53 ## @item if sys is continuous, poles with |
|
5016
|
54 ## @iftex |
|
|
55 ## @tex |
|
|
56 ## $ \vert {\rm real}(pole) \vert < ptol \Vert H \Vert $ |
|
|
57 ## @end tex |
|
|
58 ## @end iftex |
|
|
59 ## @ifinfo |
|
|
60 ## @math{ |real(pole))| < ptol*||H|| } |
|
|
61 ## @end ifinfo |
|
|
62 ## (@var{H} is appropriate Hamiltonian) |
|
3430
|
63 ## are considered to be on the imaginary axis. |
|
|
64 ## |
|
|
65 ## @item if sys is discrete, poles with |
|
5016
|
66 ## @iftex |
|
|
67 ## @tex |
|
|
68 ## $ \vert { \rm pole } - 1 \vert < ptol \Vert [ s_1 s_2 ] \Vert $ |
|
|
69 ## @end tex |
|
|
70 ## @end iftex |
|
|
71 ## @ifinfo |
|
|
72 ## @math{|abs(pole)-1| < ptol*||[s1,s2]||} |
|
|
73 ## @end ifinfo |
|
|
74 ## (appropriate symplectic pencil) |
|
|
75 ## are considered to be on the unit circle. |
|
3430
|
76 ## |
|
5016
|
77 ## @item Default value: 1e-9 |
|
3430
|
78 ## @end itemize |
|
|
79 ## @end table |
|
|
80 ## |
|
|
81 ## @strong{Outputs} |
|
|
82 ## @table @var |
|
|
83 ## @item g |
|
|
84 ## Computed gain, within @var{tol} of actual gain. @var{g} is returned as Inf |
|
|
85 ## if the system is unstable. |
|
|
86 ## @item gmin |
|
|
87 ## @itemx gmax |
|
5016
|
88 ## Actual system gain lies in the interval [@var{gmin}, @var{gmax}]. |
|
3430
|
89 ## @end table |
|
|
90 ## |
|
|
91 ## References: |
|
5016
|
92 ## Doyle, Glover, Khargonekar, Francis, @cite{State-space solutions to standard} |
|
|
93 ## @iftex |
|
|
94 ## @tex |
|
|
95 ## $ { \cal H }_2 $ @cite{and} $ { \cal H }_\infty $ |
|
|
96 ## @end tex |
|
|
97 ## @end iftex |
|
|
98 ## @ifinfo |
|
|
99 ## @cite{H-2 and H-infinity} |
|
|
100 ## @end ifinfo |
|
|
101 ## @cite{control problems}, @acronym{IEEE} @acronym{TAC} August 1989; |
|
|
102 ## Iglesias and Glover, @cite{State-Space approach to discrete-time} |
|
|
103 ## @iftex |
|
|
104 ## @tex |
|
|
105 ## $ { \cal H }_\infty $ |
|
|
106 ## @end tex |
|
|
107 ## @end iftex |
|
|
108 ## @ifinfo |
|
|
109 ## @cite{H-infinity} |
|
|
110 ## @end ifinfo |
|
|
111 ## @cite{control}, Int. J. Control, vol 54, no. 5, 1991; |
|
|
112 ## Zhou, Doyle, Glover, @cite{Robust and Optimal Control}, Prentice-Hall, 1996. |
|
3430
|
113 ## @end deftypefn |
|
|
114 |
|
|
115 function [g, gmin, gmax] = hinfnorm (sys, tol, gmin, gmax, ptol) |
|
|
116 |
|
|
117 if((nargin == 0) || (nargin > 4)) |
|
6046
|
118 print_usage (); |
|
4030
|
119 elseif(!isstruct(sys)) |
|
3430
|
120 error("Sys must be a system data structure"); |
|
|
121 endif |
|
|
122 |
|
|
123 ## set defaults where applicable |
|
|
124 if(nargin < 5) |
|
|
125 ptol = 1e-9; # pole tolerance |
|
|
126 endif |
|
|
127 if(nargin < 4) |
|
|
128 gmax = 1e9; # max gain value |
|
|
129 endif |
|
|
130 |
|
|
131 dflg = is_digital(sys); |
|
|
132 sys = sysupdate(sys,"ss"); |
|
|
133 [A,B,C,D] = sys2ss(sys); |
|
|
134 [n,nz,m,p] = sysdimensions(sys); |
|
|
135 |
|
|
136 ## eigenvalues of A must all be stable |
|
|
137 if(!is_stable(sys)) |
|
|
138 warning(["hinfnorm: unstable system (is_stable, ptol=",num2str(ptol), ... |
|
|
139 "), returning Inf"]); |
|
|
140 g = Inf; |
|
|
141 endif |
|
|
142 |
|
|
143 Dnrm = norm(D); |
|
|
144 if(nargin < 3) |
|
|
145 gmin = max(1e-9,Dnrm); # min gain value |
|
|
146 elseif(gmin < Dnrm) |
|
|
147 warning(["hinfnorm: setting Gmin=||D||=",num2str(Dnrm)]); |
|
|
148 endif |
|
|
149 |
|
|
150 if(nargin < 2) |
|
|
151 tol = 0.001; # convergence measure for gmin, gmax |
|
|
152 endif |
|
|
153 |
|
|
154 ## check for scalar input arguments 2...5 |
|
4030
|
155 if( ! (isscalar(tol) && isscalar(gmin) |
|
|
156 && isscalar(gmax) && isscalar(ptol)) ) |
|
3430
|
157 error("hinfnorm: tol, gmin, gmax, ptol must be scalars"); |
|
|
158 endif |
|
|
159 |
|
|
160 In = eye(n+nz); |
|
|
161 Im = eye(m); |
|
|
162 Ip = eye(p); |
|
|
163 ## find the Hinf norm via binary search |
|
|
164 while((gmax/gmin - 1) > tol) |
|
|
165 g = (gmax+gmin)/2; |
|
|
166 |
|
|
167 if(dflg) |
|
|
168 ## multiply g's through in formulas to avoid extreme magnitudes... |
|
|
169 Rg = g^2*Im - D'*D; |
|
|
170 Ak = A + (B/Rg)*D'*C; |
|
|
171 Ck = g^2*C'*((g^2*Ip-D*D')\C); |
|
|
172 |
|
|
173 ## set up symplectic generalized eigenvalue problem per Iglesias & Glover |
|
|
174 s1 = [Ak , zeros(nz) ; -Ck, In ]; |
|
|
175 s2 = [In, -(B/Rg)*B' ; zeros(nz) , Ak' ]; |
|
|
176 |
|
|
177 ## guard against roundoff again: zero out extremely small values |
|
|
178 ## prior to balancing |
|
|
179 s1 = s1 .* (abs(s1) > ptol*norm(s1,"inf")); |
|
|
180 s2 = s2 .* (abs(s2) > ptol*norm(s2,"inf")); |
|
|
181 [cc,dd,s1,s2] = balance(s1,s2); |
|
|
182 [qza,qzb,zz,pls] = qz(s1,s2,"S"); # ordered qz decomposition |
|
|
183 eigerr = abs(abs(pls)-1); |
|
|
184 normH = norm([s1,s2]); |
|
|
185 Hb = [s1, s2]; |
|
|
186 |
|
|
187 ## check R - B' X B condition (Iglesias and Glover's paper) |
|
|
188 X = zz((nz+1):(2*nz),1:nz)/zz(1:nz,1:nz); |
|
|
189 dcondfailed = min(real( eig(Rg - B'*X*B)) < ptol); |
|
|
190 else |
|
|
191 Rinv = inv(g*g*Im - (D' * D)); |
|
|
192 H = [A + B*Rinv*D'*C, B*Rinv*B'; ... |
|
|
193 -C'*(Ip + D*Rinv*D')*C, -(A + B*Rinv*D'*C)']; |
|
|
194 ## guard against roundoff: zero out extremely small values prior |
|
|
195 ## to balancing |
|
|
196 H = H .* (abs(H) > ptol*norm(H,"inf")); |
|
|
197 [DD,Hb] = balance(H); |
|
|
198 pls = eig(Hb); |
|
|
199 eigerr = abs(real(pls)); |
|
|
200 normH = norm(H); |
|
|
201 dcondfailed = 0; # digital condition; doesn't apply here |
|
|
202 endif |
|
|
203 if( (min(eigerr) <= ptol * normH) | dcondfailed) |
|
|
204 gmin = g; |
|
|
205 else |
|
|
206 gmax = g; |
|
|
207 endif |
|
|
208 endwhile |
|
|
209 endfunction |
