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}] =} bstmodred (@var{G}, @dots{})
20 ## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} bstmodred (@var{G}, @var{nr}, @dots{})
21 ## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} bstmodred (@var{G}, @var{opt}, @dots{})
22 ## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} bstmodred (@var{G}, @var{nr}, @var{opt}, @dots{})
24 ## Model order reduction by Balanced Stochastic Truncation (BST) 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 ## BST is a relative error method which tries to minimize
32 ## $$ || G^{-1} (G-G_r) ||_{\\infty} = min $$
48 ## LTI model to be reduced.
50 ## The desired order of the resulting reduced order system @var{Gr}.
51 ## If not specified, @var{nr} is chosen automatically according
52 ## to the description of key @var{'order'}.
54 ## Optional pairs of keys and values. @code{"key1", value1, "key2", value2}.
56 ## Optional struct with keys as field names.
57 ## Struct @var{opt} can be created directly or
58 ## by command @command{options}. @code{opt.key1 = value1, opt.key2 = value2}.
64 ## Reduced order state-space model.
66 ## Struct containing additional information.
69 ## The order of the original system @var{G}.
71 ## The order of the @var{alpha}-stable subsystem of the original system @var{G}.
73 ## The Hankel singular values of the phase system corresponding
74 ## to the @var{alpha}-stable part of the original system @var{G}.
75 ## The @var{ns} Hankel singular values are ordered decreasingly.
77 ## The order of the @var{alpha}-unstable subsystem of both the original
78 ## system @var{G} and the reduced-order system @var{Gr}.
80 ## The order of the obtained reduced order system @var{Gr}.
84 ## @strong{Option Keys and Values}
86 ## @item 'order', 'nr'
87 ## The desired order of the resulting reduced order system @var{Gr}.
88 ## If not specified, @var{nr} is the sum of NU and the number of
89 ## Hankel singular values greater than @code{MAX(TOL1,NS*EPS)};
90 ## @var{nr} can be further reduced to ensure that
91 ## @code{HSV(NR-NU) > HSV(NR+1-NU)}.
94 ## Approximation method for the H-infinity norm.
95 ## Valid values corresponding to this key are:
97 ## @item 'sr-bta', 'b'
98 ## Use the square-root Balance & Truncate method.
99 ## @item 'bfsr-bta', 'f'
100 ## Use the balancing-free square-root Balance & Truncate method. Default method.
101 ## @item 'sr-spa', 's'
102 ## Use the square-root Singular Perturbation Approximation method.
103 ## @item 'bfsr-spa', 'p'
104 ## Use the balancing-free square-root Singular Perturbation Approximation method.
108 ## Specifies the ALPHA-stability boundary for the eigenvalues
109 ## of the state dynamics matrix @var{G.A}. For a continuous-time
110 ## system, ALPHA <= 0 is the boundary value for
111 ## the real parts of eigenvalues, while for a discrete-time
112 ## system, 0 <= ALPHA <= 1 represents the
113 ## boundary value for the moduli of eigenvalues.
114 ## The ALPHA-stability domain does not include the boundary.
115 ## Default value is 0 for continuous-time systems and
116 ## 1 for discrete-time systems.
119 ## Use @code{[G, beta*I]} as new system @var{G} to combine
120 ## absolute and relative error methods.
121 ## BETA > 0 specifies the absolute/relative error weighting
122 ## parameter. A large positive value of BETA favours the
123 ## minimization of the absolute approximation error, while a
124 ## small value of BETA is appropriate for the minimization
125 ## of the relative error.
126 ## BETA = 0 means a pure relative error method and can be
127 ## used only if rank(G.D) = rows(G.D) which means that
128 ## the feedthrough matrice must not be rank-deficient.
129 ## Default value is 0.
132 ## If @var{'order'} is not specified, @var{tol1} contains the tolerance for
133 ## determining the order of reduced system.
134 ## For model reduction, the recommended value of @var{tol1} lies
135 ## in the interval [0.00001, 0.001]. @var{tol1} < 1.
136 ## If @var{tol1} <= 0 on entry, the used default value is
137 ## @var{tol1} = NS*EPS, where NS is the number of
138 ## ALPHA-stable eigenvalues of A and EPS is the machine
140 ## If @var{'order'} is specified, the value of @var{tol1} is ignored.
143 ## The tolerance for determining the order of a minimal
144 ## realization of the phase system (see METHOD) corresponding
145 ## to the ALPHA-stable part of the given system.
146 ## The recommended value is TOL2 = NS*EPS. TOL2 <= TOL1 < 1.
147 ## This value is used by default if @var{'tol2'} is not specified
148 ## or if TOL2 <= 0 on entry.
150 ## @item 'equil', 'scale'
151 ## Boolean indicating whether equilibration (scaling) should be
152 ## performed on system @var{G} prior to order reduction.
153 ## Default value is true if @code{G.scaled == false} and
154 ## false if @code{G.scaled == true}.
155 ## Note that for @acronym{MIMO} models, proper scaling of both inputs and outputs
156 ## is of utmost importance. The input and output scaling can @strong{not}
157 ## be done by the equilibration option or the @command{prescale} command
158 ## because these functions perform state transformations only.
159 ## Furthermore, signals should not be scaled simply to a certain range.
160 ## For all inputs (or outputs), a certain change should be of the same
161 ## importance for the model.
165 ## BST is often suitable to perform model reduction in order to obtain
166 ## low order design models for controller synthesis.
168 ## Approximation Properties:
171 ## Guaranteed stability of reduced models
173 ## Approximates simultaneously gain and phase
175 ## Preserves non-minimum phase zeros
177 ## Guaranteed a priori error bound
180 ## $$ || G^{-1} (G-G_r) ||_{\\infty} \\leq 2 \\sum_{j=r+1}^{n} \\frac{1+\\sigma_j}{1-\\sigma_j} - 1 $$
185 ## @strong{Algorithm}@*
186 ## Uses SLICOT AB09HD by courtesy of
187 ## @uref{http://www.slicot.org, NICONET e.V.}
190 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
191 ## Created: October 2011
194 function [Gr, info] = bstmodred (G, varargin)
200 if (! isa (G, "lti"))
201 error ("bstmodred: first argument must be an LTI system");
204 if (nargin > 1) # bstmodred (G, ...)
205 if (is_real_scalar (varargin{1})) # bstmodred (G, nr)
206 varargin = horzcat (varargin(2:end), {"order"}, varargin(1));
208 if (isstruct (varargin{1})) # bstmodred (G, opt, ...), bstmodred (G, nr, opt, ...)
209 varargin = horzcat (__opt2cell__ (varargin{1}), varargin(2:end));
211 ## order placed at the end such that nr from bstmodred (G, nr, ...)
212 ## and bstmodred (G, nr, opt, ...) overrides possible nr's from
213 ## key/value-pairs and inside opt struct (later keys override former keys,
214 ## nr > key/value > opt)
217 nkv = numel (varargin); # number of keys and values
220 error ("bstmodred: keys and values must come in pairs");
223 [a, b, c, d, tsam, scaled] = ssdata (G);
227 alpha = __modred_default_alpha__ (dt);
235 ## handle keys and values
237 key = lower (varargin{k});
241 [nr, ordsel] = __modred_check_order__ (val, rows (a));
244 tol1 = __modred_check_tol__ (val, "tol1");
247 tol2 = __modred_check_tol__ (val, "tol2");
250 alpha = __modred_check_alpha__ (val, dt);
253 if (! issample (val, 0))
254 error ("bstmodred: argument %s must be BETA >= 0", varargin{k});
258 case "method" # approximation method
259 switch (tolower (val))
260 case {"sr-bta", "b"} # 'B': use the square-root Balance & Truncate method
262 case {"bfsr-bta", "f"} # 'F': use the balancing-free square-root Balance & Truncate method
264 case {"sr-spa", "s"} # 'S': use the square-root Singular Perturbation Approximation method
266 case {"bfsr-spa", "p"} # 'P': use the balancing-free square-root Singular Perturbation Approximation method
269 error ("bstmodred: '%s' is an invalid approximation method", val);
272 case {"equil", "equilibrate", "equilibration", "scale", "scaling"}
273 scaled = __modred_check_equil__ (val);
276 warning ("bstmodred: invalid property name '%s' ignored", key);
280 ## perform model order reduction
281 [ar, br, cr, dr, nr, hsv, ns] = slab09hd (a, b, c, d, dt, scaled, job, nr, ordsel, alpha, beta, \
284 ## assemble reduced order model
285 Gr = ss (ar, br, cr, dr, tsam);
287 ## assemble info struct
290 info = struct ("n", n, "ns", ns, "hsv", hsv, "nu", nu, "nr", nr);
295 %!shared Mo, Me, Info, HSVe
296 %! A = [ -0.04165 0.0000 4.9200 -4.9200 0.0000 0.0000 0.0000
297 %! -5.2100 -12.500 0.0000 0.0000 0.0000 0.0000 0.0000
298 %! 0.0000 3.3300 -3.3300 0.0000 0.0000 0.0000 0.0000
299 %! 0.5450 0.0000 0.0000 0.0000 -0.5450 0.0000 0.0000
300 %! 0.0000 0.0000 0.0000 4.9200 -0.04165 0.0000 4.9200
301 %! 0.0000 0.0000 0.0000 0.0000 -5.2100 -12.500 0.0000
302 %! 0.0000 0.0000 0.0000 0.0000 0.0000 3.3300 -3.3300 ];
304 %! B = [ 0.0000 0.0000
312 %! C = [ 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
313 %! 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000
314 %! 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 ];
316 %! D = [ 0.0000 0.0000
320 %! G = ss (A, B, C, D, "scaled", true);
322 %! [Gr, Info] = bstmodred (G, "beta", 1.0, "tol1", 0.1, "tol2", 0.0);
323 %! [Ao, Bo, Co, Do] = ssdata (Gr);
325 %! Ae = [ 1.2729 0.0000 6.5947 0.0000 -3.4229
326 %! 0.0000 0.8169 0.0000 2.4821 0.0000
327 %! -2.9889 0.0000 -2.9028 0.0000 -0.3692
328 %! 0.0000 -3.3921 0.0000 -3.1126 0.0000
329 %! -1.4767 0.0000 -2.0339 0.0000 -0.6107 ];
331 %! Be = [ 0.1331 -0.1331
337 %! Ce = [ -0.6907 -0.6882 0.0779 0.0958 -0.0038
338 %! 0.0676 0.0000 0.6532 0.0000 -0.7522
339 %! 0.6907 -0.6882 -0.0779 0.0958 0.0038 ];
341 %! De = [ 0.0000 0.0000
345 %! HSVe = [ 0.8803 0.8506 0.8038 0.4494 0.3973 0.0214 0.0209 ].';
347 %! Mo = [Ao, Bo; Co, Do];
348 %! Me = [Ae, Be; Ce, De];
350 %!assert (Mo, Me, 1e-4);
351 %!assert (Info.hsv, HSVe, 1e-4);