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}] =} btamodred (@var{G}, @dots{})
20 ## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} btamodred (@var{G}, @var{nr}, @dots{})
21 ## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} btamodred (@var{G}, @var{opt}, @dots{})
22 ## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} btamodred (@var{G}, @var{nr}, @var{opt}, @dots{})
24 ## Model order reduction by frequency weighted Balanced Truncation Approximation (BTA) 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 ## BTA is an absolute error method which tries to minimize
32 ## $$ || G - G_r ||_{\\infty} = min $$
33 ## $$ || V \\ (G - G_r) \\ W ||_{\\infty} = 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 of the @var{alpha}-stable part of
77 ## the original system @var{G}, ordered decreasingly.
79 ## The order of the @var{alpha}-unstable subsystem of both the original
80 ## system @var{G} and the reduced-order system @var{Gr}.
82 ## The order of the obtained reduced order system @var{Gr}.
87 ## @strong{Option Keys and Values}
89 ## @item 'order', 'nr'
90 ## The desired order of the resulting reduced order system @var{Gr}.
91 ## If not specified, @var{nr} is chosen automatically such that states with
92 ## Hankel singular values @var{info.hsv} > @var{tol1} are retained.
94 ## @item 'left', 'output'
95 ## LTI model of the left/output frequency weighting @var{V}.
96 ## Default value is an identity matrix.
98 ## @item 'right', 'input'
99 ## LTI model of the right/input frequency weighting @var{W}.
100 ## Default value is an identity matrix.
103 ## Approximation method for the L-infinity norm to be used as follows:
106 ## Use the square-root Balance & Truncate method.
108 ## Use the balancing-free square-root Balance & Truncate method. Default method.
112 ## Specifies the ALPHA-stability boundary for the eigenvalues
113 ## of the state dynamics matrix @var{G.A}. For a continuous-time
114 ## system, ALPHA <= 0 is the boundary value for
115 ## the real parts of eigenvalues, while for a discrete-time
116 ## system, 0 <= ALPHA <= 1 represents the
117 ## boundary value for the moduli of eigenvalues.
118 ## The ALPHA-stability domain does not include the boundary.
119 ## Default value is 0 for continuous-time systems and
120 ## 1 for discrete-time systems.
123 ## If @var{'order'} is not specified, @var{tol1} contains the tolerance for
124 ## determining the order of the reduced model.
125 ## For model reduction, the recommended value of @var{tol1} is
126 ## c*info.hsv(1), where c lies in the interval [0.00001, 0.001].
127 ## Default value is info.ns*eps*info.hsv(1).
128 ## If @var{'order'} is specified, the value of @var{tol1} is ignored.
131 ## The tolerance for determining the order of a minimal
132 ## realization of the ALPHA-stable part of the given
133 ## model. TOL2 <= TOL1.
134 ## If not specified, ns*eps*info.hsv(1) is chosen.
137 ## Specifies the choice of frequency-weighted controllability
138 ## Grammian as follows:
141 ## Choice corresponding to a combination method [4]
142 ## of the approaches of Enns [1] and Lin-Chiu [2,3]. Default method.
144 ## Choice corresponding to the stability enhanced
145 ## modified combination method of [4].
149 ## Specifies the choice of frequency-weighted observability
150 ## Grammian as follows:
153 ## Choice corresponding to a combination method [4]
154 ## of the approaches of Enns [1] and Lin-Chiu [2,3]. Default method.
156 ## Choice corresponding to the stability enhanced
157 ## modified combination method of [4].
160 ## @item 'alpha-ctrb'
161 ## Combination method parameter for defining the
162 ## frequency-weighted controllability Grammian.
164 ## If alphac = 0, the choice of
165 ## Grammian corresponds to the method of Enns [1], while if
166 ## alphac = 1, the choice of Grammian corresponds
167 ## to the method of Lin and Chiu [2,3].
168 ## Default value is 0.
170 ## @item 'alpha-obsv'
171 ## Combination method parameter for defining the
172 ## frequency-weighted observability Grammian.
174 ## If alphao = 0, the choice of
175 ## Grammian corresponds to the method of Enns [1], while if
176 ## alphao = 1, the choice of Grammian corresponds
177 ## to the method of Lin and Chiu [2,3].
178 ## Default value is 0.
180 ## @item 'equil', 'scale'
181 ## Boolean indicating whether equilibration (scaling) should be
182 ## performed on system @var{G} prior to order reduction.
183 ## This is done by state transformations.
184 ## Default value is true if @code{G.scaled == false} and
185 ## false if @code{G.scaled == true}.
186 ## Note that for @acronym{MIMO} models, proper scaling of both inputs and outputs
187 ## is of utmost importance. The input and output scaling can @strong{not}
188 ## be done by the equilibration option or the @command{prescale} command
189 ## because these functions perform state transformations only.
190 ## Furthermore, signals should not be scaled simply to a certain range.
191 ## For all inputs (or outputs), a certain change should be of the same
192 ## importance for the model.
196 ## Approximation Properties:
199 ## Guaranteed stability of reduced models
201 ## Lower guaranteed error bound
203 ## Guaranteed a priori error bound
206 ## $$ \\sigma_{r+1} \\leq || (G-G_r) ||_{\\infty} \\leq 2 \\sum_{j=r+1}^{n} \\sigma_j $$
212 ## @strong{References}@*
214 ## Model reduction with balanced realizations: An error bound
215 ## and a frequency weighted generalization.
216 ## Proc. 23-th CDC, Las Vegas, pp. 127-132, 1984.
218 ## [2] Lin, C.-A. and Chiu, T.-Y.
219 ## Model reduction via frequency-weighted balanced realization.
220 ## Control Theory and Advanced Technology, vol. 8,
221 ## pp. 341-351, 1992.
223 ## [3] Sreeram, V., Anderson, B.D.O and Madievski, A.G.
224 ## New results on frequency weighted balanced reduction
226 ## Proc. ACC, Seattle, Washington, pp. 4004-4009, 1995.
228 ## [4] Varga, A. and Anderson, B.D.O.
229 ## Square-root balancing-free methods for the frequency-weighted
230 ## balancing related model reduction.
231 ## (report in preparation)
234 ## @strong{Algorithm}@*
235 ## Uses SLICOT AB09ID by courtesy of
236 ## @uref{http://www.slicot.org, NICONET e.V.}
239 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
240 ## Created: November 2011
243 function [Gr, info] = btamodred (varargin)
245 [Gr, info] = __modred_ab09id__ ("bta", varargin{:});
250 %!shared Mo, Me, Info, HSVe
251 %! A = [ -26.4000, 6.4023, 4.3868;
259 %! C = [ 9.2994 1.1624 0.1090 ];
263 %! G = ss (A, B, C, D); % "scaled", false
265 %! AV = [ -1.0000, 0, 4.0000, -9.2994, -1.1624, -0.1090;
266 %! 0, 2.0000, 0, -9.2994, -1.1624, -0.1090;
267 %! 0, 0, -3.0000, -9.2994, -1.1624, -0.1090;
268 %! 16.0000, 16.0000, 16.0000, -26.4000, 6.4023, 4.3868;
269 %! 0, 0, 0, 32.0000, 0, 0;
270 %! 0, 0, 0, 0, 8.0000, 0 ];
279 %! CV = [ 1 1 1 0 0 0 ];
283 %! V = ss (AV, BV, CV, DV);
285 %! [Gr, Info] = btamodred (G, 2, "left", V);
286 %! [Ao, Bo, Co, Do] = ssdata (Gr);
288 %! Ae = [ 9.1900 0.0000
289 %! 0.0000 -34.5297 ];
294 %! Ce = [ 2.8955 6.9152 ];
298 %! HSVe = [ 3.8253 0.2005 ].';
300 %! Mo = [Ao, Bo; Co, Do];
301 %! Me = [Ae, Be; Ce, De];
303 %!assert (Mo, Me, 1e-4);
304 %!assert (Info.hsv, HSVe, 1e-4);