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