1 ## Copyright (C) 2011 Lukas F. Reichlin
3 ## This file is part of LTI Syncope.
5 ## LTI Syncope is free software: you can redistribute it and/or modify
6 ## it under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation, either version 3 of the License, or
8 ## (at your option) any later version.
10 ## LTI Syncope is distributed in the hope that it will be useful,
11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 ## GNU General Public License for more details.
15 ## You should have received a copy of the GNU General Public License
16 ## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
19 ## @deftypefn{Function File} {[@var{Gr}, @var{info}] =} hnamodred (@var{G}, @dots{})
20 ## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} hnamodred (@var{G}, @var{nr}, @dots{})
21 ## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} hnamodred (@var{G}, @var{opt}, @dots{})
22 ## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} hnamodred (@var{G}, @var{nr}, @var{opt}, @dots{})
24 ## Model order reduction by frequency weighted optimal Hankel-norm (HNA) method.
25 ## The aim of model reduction is to find an LTI system @var{Gr} of order
26 ## @var{nr} (nr < n) such that the input-output behaviour of @var{Gr}
27 ## approximates the one from original system @var{G}.
29 ## HNA is an absolute error method which tries to minimize
32 ## $$ || G - G_r ||_H = min $$
33 ## $$ || V \\ (G - G_r) \\ W ||_H = min $$
41 ## ||V (G-Gr) W|| = min
45 ## where @var{V} and @var{W} denote output and input weightings.
51 ## LTI model to be reduced.
53 ## The desired order of the resulting reduced order system @var{Gr}.
54 ## If not specified, @var{nr} is chosen automatically according
55 ## to the description of key @var{"order"}.
57 ## Optional pairs of keys and values. @code{"key1", value1, "key2", value2}.
59 ## Optional struct with keys as field names.
60 ## Struct @var{opt} can be created directly or
61 ## by command @command{options}. @code{opt.key1 = value1, opt.key2 = value2}.
67 ## Reduced order state-space model.
69 ## Struct containing additional information.
72 ## The order of the original system @var{G}.
74 ## The order of the @var{alpha}-stable subsystem of the original system @var{G}.
76 ## The Hankel singular values corresponding to the projection @code{op(V)*G1*op(W)},
77 ## where G1 denotes the @var{alpha}-stable part of the original system @var{G}.
78 ## The @var{ns} Hankel singular values are ordered decreasingly.
80 ## The order of the @var{alpha}-unstable subsystem of both the original
81 ## system @var{G} and the reduced-order system @var{Gr}.
83 ## The order of the obtained reduced order system @var{Gr}.
88 ## @strong{Option Keys and Values}
90 ## @item 'order', 'nr'
91 ## The desired order of the resulting reduced order system @var{Gr}.
92 ## If not specified, @var{nr} is the sum of @var{info.nu} and the number of
93 ## Hankel singular values greater than @code{max(tol1, ns*eps*info.hsv(1)};
96 ## Specifies the computational approach to be used.
97 ## Valid values corresponding to this key are:
100 ## Use the inverse free descriptor system approach.
102 ## Use the inversion based standard approach.
104 ## Switch automatically to the inverse free
105 ## descriptor approach in case of badly conditioned
106 ## feedthrough matrices in V or W. Default method.
111 ## LTI model of the left/output frequency weighting.
112 ## The weighting must be antistable.
114 ## @math{|| V \\ (G-G_r) \\dots ||_H = min}
118 ## || V (G-Gr) . || = min
123 ## @item 'right', 'w'
124 ## LTI model of the right/input frequency weighting.
125 ## The weighting must be antistable.
127 ## @math{|| \\dots (G-G_r) \\ W ||_H = min}
131 ## || . (G-Gr) W || = min
137 ## @item 'left-inv', 'inv-v'
138 ## LTI model of the left/output frequency weighting.
139 ## The weighting must have only antistable zeros.
141 ## @math{|| inv(V) \\ (G-G_r) \\dots ||_H = min}
145 ## || inv(V) (G-Gr) . || = min
150 ## @item 'right-inv', 'inv-w'
151 ## LTI model of the right/input frequency weighting.
152 ## The weighting must have only antistable zeros.
154 ## @math{|| \\dots (G-G_r) \\ inv(W) ||_H = min}
158 ## || . (G-Gr) inv(W) || = min
164 ## @item 'left-conj', 'conj-v'
165 ## LTI model of the left/output frequency weighting.
166 ## The weighting must be stable.
168 ## @math{|| conj(V) \\ (G-G_r) \\dots ||_H = min}
172 ## || V (G-Gr) . || = min
177 ## @item 'right-conj', 'conj-w'
178 ## LTI model of the right/input frequency weighting.
179 ## The weighting must be stable.
181 ## @math{|| \\dots (G-G_r) \\ conj(W) ||_H = min}
185 ## || . (G-Gr) W || = min
191 ## @item 'left-conj-inv', 'conj-inv-v'
192 ## LTI model of the left/output frequency weighting.
193 ## The weighting must be minimum-phase.
195 ## @math{|| conj(inv(V)) \\ (G-G_r) \\dots ||_H = min}
199 ## || V (G-Gr) . || = min
204 ## @item 'right-conj-inv', 'conj-inv-w'
205 ## LTI model of the right/input frequency weighting.
206 ## The weighting must be minimum-phase.
208 ## @math{|| \\dots (G-G_r) \\ conj(inv(W)) ||_H = min}
212 ## || . (G-Gr) W || = min
219 ## Specifies the ALPHA-stability boundary for the eigenvalues
220 ## of the state dynamics matrix @var{G.A}. For a continuous-time
221 ## system, ALPHA <= 0 is the boundary value for
222 ## the real parts of eigenvalues, while for a discrete-time
223 ## system, 0 <= ALPHA <= 1 represents the
224 ## boundary value for the moduli of eigenvalues.
225 ## The ALPHA-stability domain does not include the boundary.
226 ## Default value is 0 for continuous-time systems and
227 ## 1 for discrete-time systems.
230 ## If @var{'order'} is not specified, @var{tol1} contains the tolerance for
231 ## determining the order of the reduced model.
232 ## For model reduction, the recommended value of @var{tol1} is
233 ## c*info.hsv(1), where c lies in the interval [0.00001, 0.001].
235 ## If @var{'order'} is specified, the value of @var{tol1} is ignored.
238 ## The tolerance for determining the order of a minimal
239 ## realization of the ALPHA-stable part of the given
240 ## model. @var{tol2} <= @var{tol1} < 1.
241 ## If not specified, ns*eps*info.hsv(1) is chosen.
243 ## @item 'equil', 'scale'
244 ## Boolean indicating whether equilibration (scaling) should be
245 ## performed on system @var{G} prior to order reduction.
246 ## Default value is true if @code{G.scaled == false} and
247 ## false if @code{G.scaled == true}.
248 ## Note that for @acronym{MIMO} models, proper scaling of both inputs and outputs
249 ## is of utmost importance. The input and output scaling can @strong{not}
250 ## be done by the equilibration option or the @command{prescale} command
251 ## because these functions perform state transformations only.
252 ## Furthermore, signals should not be scaled simply to a certain range.
253 ## For all inputs (or outputs), a certain change should be of the same
254 ## importance for the model.
258 ## Approximation Properties:
261 ## Guaranteed stability of reduced models
263 ## Lower guaranteed error bound
265 ## Guaranteed a priori error bound
268 ## $$ \\sigma_{r+1} \\leq || (G-G_r) ||_{\\infty} \\leq 2 \\sum_{j=r+1}^{n} \\sigma_j $$
273 ## @strong{Algorithm}@*
274 ## Uses SLICOT AB09JD by courtesy of
275 ## @uref{http://www.slicot.org, NICONET e.V.}
278 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
279 ## Created: October 2011
282 function [Gr, info] = hnamodred (G, varargin)
288 if (! isa (G, "lti"))
289 error ("hnamodred: first argument must be an LTI system");
292 if (nargin > 1) # hnamodred (G, ...)
293 if (is_real_scalar (varargin{1})) # hnamodred (G, nr)
294 varargin = horzcat (varargin(2:end), {"order"}, varargin(1));
296 if (isstruct (varargin{1})) # hnamodred (G, opt, ...), hnamodred (G, nr, opt, ...)
297 varargin = horzcat (__opt2cell__ (varargin{1}), varargin(2:end));
299 ## order placed at the end such that nr from hnamodred (G, nr, ...)
300 ## and hnamodred (G, nr, opt, ...) overrides possible nr's from
301 ## key/value-pairs and inside opt struct (later keys override former keys,
302 ## nr > key/value > opt)
305 nkv = numel (varargin); # number of keys and values
308 error ("hnamodred: keys and values must come in pairs");
311 [a, b, c, d, tsam, scaled] = ssdata (G);
316 alpha = __modred_default_alpha__ (dt);
317 av = bv = cv = dv = [];
319 aw = bw = cw = dw = [];
327 ## handle keys and values
329 key = lower (varargin{k});
332 case {"left", "v", "wo"}
333 [av, bv, cv, dv, jobv] = __modred_check_weight__ (val, dt, p, p);
334 ## TODO: correct error messages for non-square weights
336 case {"right", "w", "wi"}
337 [aw, bw, cw, dw, jobw] = __modred_check_weight__ (val, dt, m, m);
339 case {"left-inv", "inv-v"}
340 [av, bv, cv, dv] = __modred_check_weight__ (val, dt, p, p);
343 case {"right-inv", "inv-w"}
344 [aw, bw, cw, dw] = __modred_check_weight__ (val, dt, m, m);
347 case {"left-conj", "conj-v"}
348 [av, bv, cv, dv] = __modred_check_weight__ (val, dt, p, p);
351 case {"right-conj", "conj-w"}
352 [aw, bw, cw, dw] = __modred_check_weight__ (val, dt, m, m);
355 case {"left-conj-inv", "conj-inv-v"}
356 [av, bv, cv, dv] = __modred_check_weight__ (val, dt, p, p);
359 case {"right-conj-inv", "conj-inv-w"}
360 [aw, bw, cw, dw] = __modred_check_weight__ (val, dt, m, m);
364 [nr, ordsel] = __modred_check_order__ (val, rows (a));
367 tol1 = __modred_check_tol__ (val, "tol1");
370 tol2 = __modred_check_tol__ (val, "tol2");
373 alpha = __modred_check_alpha__ (val, dt);
376 switch (tolower (val(1)))
377 case {"d", "n"} # "descriptor"
379 case {"s", "i"} # "standard"
381 case "a" # {"auto", "automatic"}
384 error ("hnamodred: invalid computational approach");
387 case {"equil", "equilibrate", "equilibration", "scale", "scaling"}
388 scaled = __modred_check_equil__ (val);
391 warning ("hnamodred: invalid property name '%s' ignored", key);
396 ## perform model order reduction
397 [ar, br, cr, dr, nr, hsv, ns] = slab09jd (a, b, c, d, dt, scaled, nr, ordsel, alpha, \
398 jobv, av, bv, cv, dv, \
399 jobw, aw, bw, cw, dw, \
402 ## assemble reduced order model
403 Gr = ss (ar, br, cr, dr, tsam);
405 ## assemble info struct
408 info = struct ("n", n, "ns", ns, "hsv", hsv, "nu", nu, "nr", nr);
413 %!shared Mo, Me, Info, HSVe
414 %! A = [ -3.8637 -7.4641 -9.1416 -7.4641 -3.8637 -1.0000
419 %! 0 0 0 0 1.0000 0 ];
428 %! C = [ 0 0 0 0 0 1 ];
432 %! G = ss (A, B, C, D); # "scaled", false
434 %! AV = [ 0.2000 -1.0000
440 %! CV = [ -1.8000 0 ];
444 %! V = ss (AV, BV, CV, DV);
446 %! [Gr, Info] = hnamodred (G, "left", V, "tol1", 1e-1, "tol2", 1e-14);
447 %! [Ao, Bo, Co, Do] = ssdata (Gr);
449 %! Ae = [ -0.2391 0.3072 1.1630 1.1967
450 %! -2.9709 -0.2391 2.6270 3.1027
451 %! 0.0000 0.0000 -0.5137 -1.2842
452 %! 0.0000 0.0000 0.1519 -0.5137 ];
459 %! Ce = [ -0.4466 0.0143 -0.4780 -0.2013 ];
463 %! HSVe = [ 2.6790 2.1589 0.8424 0.1929 0.0219 0.0011 ].';
465 %! Mo = [Ao, Bo; Co, Do];
466 %! Me = [Ae, Be; Ce, De];
468 %!assert (Mo, Me, 1e-4);
469 %!assert (Info.hsv, HSVe, 1e-4);