annotate scripts/control/hinf/hinfnorm.m @ 7131:a184bc985c37

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