octave_packages/control-2.3.52/bstmodred.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}] =} 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{})
  23 ##
  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}.
  28 ##
  29 ## BST is a relative error method which tries to minimize
  30 ## @iftex
  31 ## @tex
  32 ## $$ || G^{-1} (G-G_r) ||_{\\infty} = min $$
  33 ## @end tex
  34 ## @end iftex
  35 ## @ifnottex
  36 ## @example
  37 ##    -1
  38 ## ||G  (G-Gr)||    = min
  39 ##              inf
  40 ## @end example
  41 ## @end ifnottex
  42 ##
  43 ##
  44 ##
  45 ## @strong{Inputs}
  46 ## @table @var
  47 ## @item G
  48 ## LTI model to be reduced.
  49 ## @item nr
  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'}.
  53 ## @item @dots{}
  54 ## Optional pairs of keys and values.  @code{"key1", value1, "key2", value2}.
  55 ## @item opt
  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}.
  59 ## @end table
  60 ##
  61 ## @strong{Outputs}
  62 ## @table @var
  63 ## @item Gr
  64 ## Reduced order state-space model.
  65 ## @item info
  66 ## Struct containing additional information.
  67 ## @table @var
  68 ## @item info.n
  69 ## The order of the original system @var{G}.
  70 ## @item info.ns
  71 ## The order of the @var{alpha}-stable subsystem of the original system @var{G}.
  72 ## @item info.hsv
  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.
  76 ## @item info.nu
  77 ## The order of the @var{alpha}-unstable subsystem of both the original
  78 ## system @var{G} and the reduced-order system @var{Gr}.
  79 ## @item info.nr
  80 ## The order of the obtained reduced order system @var{Gr}.
  81 ## @end table
  82 ## @end table
  83 ##
  84 ## @strong{Option Keys and Values}
  85 ## @table @var
  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)}.
  92 ##
  93 ## @item 'method'
  94 ## Approximation method for the H-infinity norm.
  95 ## Valid values corresponding to this key are:
  96 ## @table @var
  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.
 105 ## @end table
 106 ##
 107 ## @item 'alpha'
 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.
 117 ##
 118 ## @item 'beta'
 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.
 130 ##
 131 ## @item 'tol1'
 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
 139 ## precision.
 140 ## If @var{'order'} is specified, the value of @var{tol1} is ignored.
 141 ##
 142 ## @item 'tol2'
 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.
 149 ##
 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.
 162 ## @end table
 163 ##
 164 ##
 165 ## BST is often suitable to perform model reduction in order to obtain
 166 ## low order design models for controller synthesis.
 167 ##
 168 ## Approximation Properties:
 169 ## @itemize @bullet
 170 ## @item
 171 ## Guaranteed stability of reduced models
 172 ## @item
 173 ## Approximates simultaneously gain and phase
 174 ## @item
 175 ## Preserves non-minimum phase zeros
 176 ## @item
 177 ## Guaranteed a priori error bound
 178 ## @iftex
 179 ## @tex
 180 ## $$ || G^{-1} (G-G_r) ||_{\\infty} \\leq 2 \\sum_{j=r+1}^{n} \\frac{1+\\sigma_j}{1-\\sigma_j} - 1 $$
 181 ## @end tex
 182 ## @end iftex
 183 ## @end itemize
 184 ##
 185 ## @strong{Algorithm}@*
 186 ## Uses SLICOT AB09HD by courtesy of
 187 ## @uref{http://www.slicot.org, NICONET e.V.}
 188 ## @end deftypefn
 189
 190 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
 191 ## Created: October 2011
 192 ## Version: 0.1
 193
 194 function [Gr, info] = bstmodred (G, varargin)
 195
 196   if (nargin == 0)
 197     print_usage ();
 198   endif
 199
 200   if (! isa (G, "lti"))
 201     error ("bstmodred: first argument must be an LTI system");
 202   endif
 203
 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));
 207     endif
 208     if (isstruct (varargin{1}))                    # bstmodred (G, opt, ...), bstmodred (G, nr, opt, ...)
 209       varargin = horzcat (__opt2cell__ (varargin{1}), varargin(2:end));
 210     endif
 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)
 215   endif
 216
 217   nkv = numel (varargin);                          # number of keys and values
 218
 219   if (rem (nkv, 2))
 220     error ("bstmodred: keys and values must come in pairs");
 221   endif
 222
 223   [a, b, c, d, tsam, scaled] = ssdata (G);
 224   dt = isdt (G);
 225
 226   ## default arguments
 227   alpha = __modred_default_alpha__ (dt);
 228   beta = 0;
 229   tol1 = 0;
 230   tol2 = 0;
 231   ordsel = 1;
 232   nr = 0;
 233   job = 1;
 234
 235   ## handle keys and values
 236   for k = 1 : 2 : nkv
 237     key = lower (varargin{k});
 238     val = varargin{k+1};
 239     switch (key)
 240       case {"order", "nr"}
 241         [nr, ordsel] = __modred_check_order__ (val, rows (a));
 242
 243       case "tol1"
 244         tol1 = __modred_check_tol__ (val, "tol1");
 245
 246       case "tol2"
 247         tol2 = __modred_check_tol__ (val, "tol2");
 248
 249       case "alpha"
 250         alpha = __modred_check_alpha__ (val, dt);
 251
 252       case "beta"
 253         if (! issample (val, 0))
 254           error ("bstmodred: argument %s must be BETA >= 0", varargin{k});
 255         endif
 256         beta = val;
 257
 258       case "method"                  # approximation method
 259         switch (tolower (val))
 260           case {"sr-bta", "b"}       # 'B':  use the square-root Balance & Truncate method
 261             job = 0;
 262           case {"bfsr-bta", "f"}     # 'F':  use the balancing-free square-root Balance & Truncate method
 263             job = 1;
 264           case {"sr-spa", "s"}       # 'S':  use the square-root Singular Perturbation Approximation method
 265             job = 2;
 266           case {"bfsr-spa", "p"}     # 'P':  use the balancing-free square-root Singular Perturbation Approximation method
 267             job = 3;
 268           otherwise
 269             error ("bstmodred: '%s' is an invalid approximation method", val);
 270         endswitch
 271
 272       case {"equil", "equilibrate", "equilibration", "scale", "scaling"}
 273         scaled = __modred_check_equil__ (val);
 274
 275       otherwise
 276         warning ("bstmodred: invalid property name '%s' ignored", key);
 277     endswitch
 278   endfor
 279
 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, \
 282                                             tol1, tol2);
 283
 284   ## assemble reduced order model
 285   Gr = ss (ar, br, cr, dr, tsam);
 286
 287   ## assemble info struct
 288   n = rows (a);
 289   nu = n - ns;
 290   info = struct ("n", n, "ns", ns, "hsv", hsv, "nu", nu, "nr", nr);
 291
 292 endfunction
 293
 294
 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 ];
 303 %!
 304 %! B =  [  0.0000   0.0000
 305 %!         12.500   0.0000
 306 %!         0.0000   0.0000
 307 %!         0.0000   0.0000
 308 %!         0.0000   0.0000
 309 %!         0.0000   12.500
 310 %!         0.0000   0.0000 ];
 311 %!
 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 ];
 315 %!
 316 %! D =  [  0.0000   0.0000
 317 %!         0.0000   0.0000
 318 %!         0.0000   0.0000 ];
 319 %!
 320 %! G = ss (A, B, C, D, "scaled", true);
 321 %!
 322 %! [Gr, Info] = bstmodred (G, "beta", 1.0, "tol1", 0.1, "tol2", 0.0);
 323 %! [Ao, Bo, Co, Do] = ssdata (Gr);
 324 %!
 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 ];
 330 %!
 331 %! Be = [  0.1331  -0.1331
 332 %!        -0.0862  -0.0862
 333 %!        -2.6777   2.6777
 334 %!        -3.5767  -3.5767
 335 %!        -2.3033   2.3033 ];
 336 %!
 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 ];
 340 %!
 341 %! De = [  0.0000   0.0000
 342 %!         0.0000   0.0000
 343 %!         0.0000   0.0000 ];
 344 %!
 345 %! HSVe = [  0.8803   0.8506   0.8038   0.4494   0.3973   0.0214   0.0209 ].';
 346 %!
 347 %! Mo = [Ao, Bo; Co, Do];
 348 %! Me = [Ae, Be; Ce, De];
 349 %!
 350 %!assert (Mo, Me, 1e-4);
 351 %!assert (Info.hsv, HSVe, 1e-4);