octave_packages/control-2.3.52/btamodred.m

   1 ## Copyright (C) 2011   Lukas F. Reichlin
   2 ##
   3 ## This file is part of LTI Syncope.
   4 ##
   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.
   9 ##
  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.
  14 ##
  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/>.
  17
  18 ## -*- texinfo -*-
  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{})
  23 ##
  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}.
  28 ##
  29 ## BTA is an absolute error method which tries to minimize
  30 ## @iftex
  31 ## @tex
  32 ## $$ || G - G_r ||_{\\infty} = min $$
  33 ## $$ || V \\ (G - G_r) \\ W ||_{\\infty} = min $$
  34 ## @end tex
  35 ## @end iftex
  36 ## @ifnottex
  37 ## @example
  38 ## ||G-Gr||    = min
  39 ##         inf
  40 ##
  41 ## ||V (G-Gr) W||    = min
  42 ##               inf
  43 ## @end example
  44 ## @end ifnottex
  45 ## where @var{V} and @var{W} denote output and input weightings.
  46 ##
  47 ##
  48 ## @strong{Inputs}
  49 ## @table @var
  50 ## @item G
  51 ## LTI model to be reduced.
  52 ## @item nr
  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'}.
  56 ## @item @dots{}
  57 ## Optional pairs of keys and values.  @code{"key1", value1, "key2", value2}.
  58 ## @item opt
  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}.
  62 ## @end table
  63 ##
  64 ## @strong{Outputs}
  65 ## @table @var
  66 ## @item Gr
  67 ## Reduced order state-space model.
  68 ## @item info
  69 ## Struct containing additional information.
  70 ## @table @var
  71 ## @item info.n
  72 ## The order of the original system @var{G}.
  73 ## @item info.ns
  74 ## The order of the @var{alpha}-stable subsystem of the original system @var{G}.
  75 ## @item info.hsv
  76 ## The Hankel singular values of the @var{alpha}-stable part of
  77 ## the original system @var{G}, ordered decreasingly.
  78 ## @item info.nu
  79 ## The order of the @var{alpha}-unstable subsystem of both the original
  80 ## system @var{G} and the reduced-order system @var{Gr}.
  81 ## @item info.nr
  82 ## The order of the obtained reduced order system @var{Gr}.
  83 ## @end table
  84 ## @end table
  85 ##
  86 ##
  87 ## @strong{Option Keys and Values}
  88 ## @table @var
  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.
  93 ##
  94 ## @item 'left', 'output'
  95 ## LTI model of the left/output frequency weighting @var{V}.
  96 ## Default value is an identity matrix.
  97 ##
  98 ## @item 'right', 'input'
  99 ## LTI model of the right/input frequency weighting @var{W}.
 100 ## Default value is an identity matrix.
 101 ##
 102 ## @item 'method'
 103 ## Approximation method for the L-infinity norm to be used as follows:
 104 ## @table @var
 105 ## @item 'sr', 'b'
 106 ## Use the square-root Balance & Truncate method.
 107 ## @item 'bfsr', 'f'
 108 ## Use the balancing-free square-root Balance & Truncate method.  Default method.
 109 ## @end table
 110 ##
 111 ## @item 'alpha'
 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.
 121 ##
 122 ## @item 'tol1'
 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.
 129 ##
 130 ## @item 'tol2'
 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.
 135 ##
 136 ## @item 'gram-ctrb'
 137 ## Specifies the choice of frequency-weighted controllability
 138 ## Grammian as follows:
 139 ## @table @var
 140 ## @item 'standard'
 141 ## Choice corresponding to a combination method [4]
 142 ## of the approaches of Enns [1] and Lin-Chiu [2,3].  Default method.
 143 ## @item 'enhanced'
 144 ## Choice corresponding to the stability enhanced
 145 ## modified combination method of [4].
 146 ## @end table
 147 ##
 148 ## @item 'gram-obsv'
 149 ## Specifies the choice of frequency-weighted observability
 150 ## Grammian as follows:
 151 ## @table @var
 152 ## @item 'standard'
 153 ## Choice corresponding to a combination method [4]
 154 ## of the approaches of Enns [1] and Lin-Chiu [2,3].  Default method.
 155 ## @item 'enhanced'
 156 ## Choice corresponding to the stability enhanced
 157 ## modified combination method of [4].
 158 ## @end table
 159 ##
 160 ## @item 'alpha-ctrb'
 161 ## Combination method parameter for defining the
 162 ## frequency-weighted controllability Grammian.
 163 ## abs(alphac) <= 1.
 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.
 169 ##
 170 ## @item 'alpha-obsv'
 171 ## Combination method parameter for defining the
 172 ## frequency-weighted observability Grammian.
 173 ## abs(alphao) <= 1.
 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.
 179 ##
 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.
 193 ## @end table
 194 ##
 195 ##
 196 ## Approximation Properties:
 197 ## @itemize @bullet
 198 ## @item
 199 ## Guaranteed stability of reduced models
 200 ## @item
 201 ## Lower guaranteed error bound
 202 ## @item
 203 ## Guaranteed a priori error bound
 204 ## @iftex
 205 ## @tex
 206 ## $$ \\sigma_{r+1} \\leq || (G-G_r) ||_{\\infty} \\leq 2 \\sum_{j=r+1}^{n} \\sigma_j $$
 207 ## @end tex
 208 ## @end iftex
 209 ## @end itemize
 210 ##
 211 ##
 212 ## @strong{References}@*
 213 ## [1] Enns, D.
 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.
 217 ##
 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.
 222 ##
 223 ## [3] Sreeram, V., Anderson, B.D.O and Madievski, A.G.
 224 ## New results on frequency weighted balanced reduction
 225 ## technique.
 226 ## Proc. ACC, Seattle, Washington, pp. 4004-4009, 1995.
 227 ##
 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)
 232 ##
 233 ##
 234 ## @strong{Algorithm}@*
 235 ## Uses SLICOT AB09ID by courtesy of
 236 ## @uref{http://www.slicot.org, NICONET e.V.}
 237 ## @end deftypefn
 238
 239 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
 240 ## Created: November 2011
 241 ## Version: 0.1
 242
 243 function [Gr, info] = btamodred (varargin)
 244
 245   [Gr, info] = __modred_ab09id__ ("bta", varargin{:});
 246
 247 endfunction
 248
 249
 250 %!shared Mo, Me, Info, HSVe
 251 %! A =  [ -26.4000,    6.4023,    4.3868;
 252 %!         32.0000,         0,         0;
 253 %!               0,    8.0000,         0 ];
 254 %!
 255 %! B =  [       16
 256 %!               0
 257 %!               0 ];
 258 %!
 259 %! C =  [   9.2994     1.1624     0.1090 ];
 260 %!
 261 %! D =  [        0 ];
 262 %!
 263 %! G = ss (A, B, C, D);  % "scaled", false
 264 %!
 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 ];
 271 %!
 272 %! BV = [        1
 273 %!               1
 274 %!               1
 275 %!               0
 276 %!               0
 277 %!               0 ];
 278 %!
 279 %! CV = [        1          1          1          0          0          0 ];
 280 %!
 281 %! DV = [        0 ];
 282 %!
 283 %! V = ss (AV, BV, CV, DV);
 284 %!
 285 %! [Gr, Info] = btamodred (G, 2, "left", V);
 286 %! [Ao, Bo, Co, Do] = ssdata (Gr);
 287 %!
 288 %! Ae = [  9.1900   0.0000
 289 %!         0.0000 -34.5297 ];
 290 %!
 291 %! Be = [ 11.9593
 292 %!        16.9329 ];
 293 %!
 294 %! Ce = [  2.8955   6.9152 ];
 295 %!
 296 %! De = [  0.0000 ];
 297 %!
 298 %! HSVe = [  3.8253   0.2005 ].';
 299 %!
 300 %! Mo = [Ao, Bo; Co, Do];
 301 %! Me = [Ae, Be; Ce, De];
 302 %!
 303 %!assert (Mo, Me, 1e-4);
 304 %!assert (Info.hsv, HSVe, 1e-4);