octave_packages/control-2.3.52/hnamodred.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}] =} 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{})
  23 ##
  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}.
  28 ##
  29 ## HNA is an absolute error method which tries to minimize
  30 ## @iftex
  31 ## @tex
  32 ## $$ || G - G_r ||_H = min $$
  33 ## $$ || V \\ (G - G_r) \\ W ||_H = min $$
  34 ## @end tex
  35 ## @end iftex
  36 ## @ifnottex
  37 ## @example
  38 ## ||G-Gr||  = min
  39 ##         H
  40 ##
  41 ## ||V (G-Gr) W||  = min
  42 ##               H
  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 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.
  79 ## @item info.nu
  80 ## The order of the @var{alpha}-unstable subsystem of both the original
  81 ## system @var{G} and the reduced-order system @var{Gr}.
  82 ## @item info.nr
  83 ## The order of the obtained reduced order system @var{Gr}.
  84 ## @end table
  85 ## @end table
  86 ##
  87 ##
  88 ## @strong{Option Keys and Values}
  89 ## @table @var
  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)};
  94 ##
  95 ## @item 'method'
  96 ## Specifies the computational approach to be used.
  97 ## Valid values corresponding to this key are:
  98 ## @table @var
  99 ## @item 'descriptor'
 100 ## Use the inverse free descriptor system approach.
 101 ## @item 'standard'
 102 ## Use the inversion based standard approach.
 103 ## @item 'auto'
 104 ## Switch automatically to the inverse free
 105 ## descriptor approach in case of badly conditioned
 106 ## feedthrough matrices in V or W.  Default method.
 107 ## @end table
 108 ##
 109 ##
 110 ## @item 'left', 'v'
 111 ## LTI model of the left/output frequency weighting.
 112 ## The weighting must be antistable.
 113 ## @iftex
 114 ## @math{|| V \\ (G-G_r) \\dots ||_H = min}
 115 ## @end iftex
 116 ## @ifnottex
 117 ## @example
 118 ## || V (G-Gr) . ||  = min
 119 ##                 H
 120 ## @end example
 121 ## @end ifnottex
 122 ##
 123 ## @item 'right', 'w'
 124 ## LTI model of the right/input frequency weighting.
 125 ## The weighting must be antistable.
 126 ## @iftex
 127 ## @math{|| \\dots (G-G_r) \\ W ||_H = min}
 128 ## @end iftex
 129 ## @ifnottex
 130 ## @example
 131 ## || . (G-Gr) W ||  = min
 132 ##                 H
 133 ## @end example
 134 ## @end ifnottex
 135 ##
 136 ##
 137 ## @item 'left-inv', 'inv-v'
 138 ## LTI model of the left/output frequency weighting.
 139 ## The weighting must have only antistable zeros.
 140 ## @iftex
 141 ## @math{|| inv(V) \\ (G-G_r) \\dots ||_H = min}
 142 ## @end iftex
 143 ## @ifnottex
 144 ## @example
 145 ## || inv(V) (G-Gr) . ||  = min
 146 ##                      H
 147 ## @end example
 148 ## @end ifnottex
 149 ##
 150 ## @item 'right-inv', 'inv-w'
 151 ## LTI model of the right/input frequency weighting.
 152 ## The weighting must have only antistable zeros.
 153 ## @iftex
 154 ## @math{|| \\dots (G-G_r) \\ inv(W) ||_H = min}
 155 ## @end iftex
 156 ## @ifnottex
 157 ## @example
 158 ## || . (G-Gr) inv(W) ||  = min
 159 ##                      H
 160 ## @end example
 161 ## @end ifnottex
 162 ##
 163 ##
 164 ## @item 'left-conj', 'conj-v'
 165 ## LTI model of the left/output frequency weighting.
 166 ## The weighting must be stable.
 167 ## @iftex
 168 ## @math{|| conj(V) \\ (G-G_r) \\dots ||_H = min}
 169 ## @end iftex
 170 ## @ifnottex
 171 ## @example
 172 ## || V (G-Gr) . ||  = min
 173 ##                 H
 174 ## @end example
 175 ## @end ifnottex
 176 ##
 177 ## @item 'right-conj', 'conj-w'
 178 ## LTI model of the right/input frequency weighting.
 179 ## The weighting must be stable.
 180 ## @iftex
 181 ## @math{|| \\dots (G-G_r) \\ conj(W) ||_H = min}
 182 ## @end iftex
 183 ## @ifnottex
 184 ## @example
 185 ## || . (G-Gr) W ||  = min
 186 ##                 H
 187 ## @end example
 188 ## @end ifnottex
 189 ##
 190 ##
 191 ## @item 'left-conj-inv', 'conj-inv-v'
 192 ## LTI model of the left/output frequency weighting.
 193 ## The weighting must be minimum-phase.
 194 ## @iftex
 195 ## @math{|| conj(inv(V)) \\ (G-G_r) \\dots ||_H = min}
 196 ## @end iftex
 197 ## @ifnottex
 198 ## @example
 199 ## || V (G-Gr) . ||  = min
 200 ##                 H
 201 ## @end example
 202 ## @end ifnottex
 203 ##
 204 ## @item 'right-conj-inv', 'conj-inv-w'
 205 ## LTI model of the right/input frequency weighting.
 206 ## The weighting must be minimum-phase.
 207 ## @iftex
 208 ## @math{|| \\dots (G-G_r) \\ conj(inv(W)) ||_H = min}
 209 ## @end iftex
 210 ## @ifnottex
 211 ## @example
 212 ## || . (G-Gr) W ||  = min
 213 ##                 H
 214 ## @end example
 215 ## @end ifnottex
 216 ##
 217 ##
 218 ## @item 'alpha'
 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.
 228 ##
 229 ## @item 'tol1'
 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].
 234 ## @var{tol1} < 1.
 235 ## If @var{'order'} is specified, the value of @var{tol1} is ignored.
 236 ##
 237 ## @item 'tol2'
 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.
 242 ##
 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.
 255 ## @end table
 256 ##
 257 ##
 258 ## Approximation Properties:
 259 ## @itemize @bullet
 260 ## @item
 261 ## Guaranteed stability of reduced models
 262 ## @item
 263 ## Lower guaranteed error bound
 264 ## @item
 265 ## Guaranteed a priori error bound
 266 ## @iftex
 267 ## @tex
 268 ## $$ \\sigma_{r+1} \\leq || (G-G_r) ||_{\\infty} \\leq 2 \\sum_{j=r+1}^{n} \\sigma_j $$
 269 ## @end tex
 270 ## @end iftex
 271 ## @end itemize
 272 ##
 273 ## @strong{Algorithm}@*
 274 ## Uses SLICOT AB09JD by courtesy of
 275 ## @uref{http://www.slicot.org, NICONET e.V.}
 276 ## @end deftypefn
 277
 278 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
 279 ## Created: October 2011
 280 ## Version: 0.1
 281
 282 function [Gr, info] = hnamodred (G, varargin)
 283
 284   if (nargin == 0)
 285     print_usage ();
 286   endif
 287
 288   if (! isa (G, "lti"))
 289     error ("hnamodred: first argument must be an LTI system");
 290   endif
 291
 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));
 295     endif
 296     if (isstruct (varargin{1}))                    # hnamodred (G, opt, ...), hnamodred (G, nr, opt, ...)
 297       varargin = horzcat (__opt2cell__ (varargin{1}), varargin(2:end));
 298     endif
 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)
 303   endif
 304
 305   nkv = numel (varargin);                          # number of keys and values
 306
 307   if (rem (nkv, 2))
 308     error ("hnamodred: keys and values must come in pairs");
 309   endif
 310
 311   [a, b, c, d, tsam, scaled] = ssdata (G);
 312   [p, m] = size (G);
 313   dt = isdt (G);
 314
 315   ## default arguments
 316   alpha = __modred_default_alpha__ (dt);
 317   av = bv = cv = dv = [];
 318   jobv = 0;
 319   aw = bw = cw = dw = [];
 320   jobw = 0;
 321   jobinv = 2;
 322   tol1 = 0;
 323   tol2 = 0;
 324   ordsel = 1;
 325   nr = 0;
 326
 327   ## handle keys and values
 328   for k = 1 : 2 : nkv
 329     key = lower (varargin{k});
 330     val = varargin{k+1};
 331     switch (key)
 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
 335
 336       case {"right", "w", "wi"}
 337         [aw, bw, cw, dw, jobw] = __modred_check_weight__ (val, dt, m, m);
 338
 339       case {"left-inv", "inv-v"}
 340         [av, bv, cv, dv] = __modred_check_weight__ (val, dt, p, p);
 341         jobv = 2;
 342
 343       case {"right-inv", "inv-w"}
 344         [aw, bw, cw, dw] = __modred_check_weight__ (val, dt, m, m);
 345         jobv = 2
 346
 347       case {"left-conj", "conj-v"}
 348         [av, bv, cv, dv] = __modred_check_weight__ (val, dt, p, p);
 349         jobv = 3;
 350
 351       case {"right-conj", "conj-w"}
 352         [aw, bw, cw, dw] = __modred_check_weight__ (val, dt, m, m);
 353         jobv = 3
 354
 355       case {"left-conj-inv", "conj-inv-v"}
 356         [av, bv, cv, dv] = __modred_check_weight__ (val, dt, p, p);
 357         jobv = 4;
 358
 359       case {"right-conj-inv", "conj-inv-w"}
 360         [aw, bw, cw, dw] = __modred_check_weight__ (val, dt, m, m);
 361         jobv = 4
 362
 363       case {"order", "nr"}
 364         [nr, ordsel] = __modred_check_order__ (val, rows (a));
 365
 366       case "tol1"
 367         tol1 = __modred_check_tol__ (val, "tol1");
 368
 369       case "tol2"
 370         tol2 = __modred_check_tol__ (val, "tol2");
 371
 372       case "alpha"
 373         alpha = __modred_check_alpha__ (val, dt);
 374
 375       case "method"
 376         switch (tolower (val(1)))
 377           case {"d", "n"}      # "descriptor"
 378             jobinv = 0;
 379           case {"s", "i"}      # "standard"
 380             jobinv = 1;
 381           case "a"             # {"auto", "automatic"}
 382             jobinv = 2;
 383           otherwise
 384             error ("hnamodred: invalid computational approach");
 385         endswitch
 386
 387       case {"equil", "equilibrate", "equilibration", "scale", "scaling"}
 388         scaled = __modred_check_equil__ (val);
 389
 390       otherwise
 391         warning ("hnamodred: invalid property name '%s' ignored", key);
 392     endswitch
 393   endfor
 394
 395
 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, \
 400                                             jobinv, tol1, tol2);
 401
 402   ## assemble reduced order model
 403   Gr = ss (ar, br, cr, dr, tsam);
 404
 405   ## assemble info struct
 406   n = rows (a);
 407   nu = n - ns;
 408   info = struct ("n", n, "ns", ns, "hsv", hsv, "nu", nu, "nr", nr);
 409
 410 endfunction
 411
 412
 413 %!shared Mo, Me, Info, HSVe
 414 %! A =  [ -3.8637   -7.4641   -9.1416   -7.4641   -3.8637   -1.0000
 415 %!         1.0000,         0         0         0         0         0
 416 %!              0    1.0000         0         0         0         0
 417 %!              0         0    1.0000         0         0         0
 418 %!              0         0         0    1.0000         0         0
 419 %!              0         0         0         0    1.0000         0 ];
 420 %!
 421 %! B =  [       1
 422 %!              0
 423 %!              0
 424 %!              0
 425 %!              0
 426 %!              0 ];
 427 %!
 428 %! C =  [       0         0         0         0         0         1 ];
 429 %!
 430 %! D =  [       0 ];
 431 %!
 432 %! G = ss (A, B, C, D);  # "scaled", false
 433 %!
 434 %! AV = [  0.2000   -1.0000
 435 %!         1.0000         0 ];
 436 %!
 437 %! BV = [       1
 438 %!              0 ];
 439 %!
 440 %! CV = [ -1.8000         0 ];
 441 %!
 442 %! DV = [       1 ];
 443 %!
 444 %! V = ss (AV, BV, CV, DV);
 445 %!
 446 %! [Gr, Info] = hnamodred (G, "left", V, "tol1", 1e-1, "tol2", 1e-14);
 447 %! [Ao, Bo, Co, Do] = ssdata (Gr);
 448 %!
 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 ];
 453 %!
 454 %! Be = [ -1.0497
 455 %!        -3.7052
 456 %!         0.8223
 457 %!         0.7435 ];
 458 %!
 459 %! Ce = [ -0.4466   0.0143  -0.4780  -0.2013 ];
 460 %!
 461 %! De = [  0.0219 ];
 462 %!
 463 %! HSVe = [  2.6790   2.1589   0.8424   0.1929   0.0219   0.0011 ].';
 464 %!
 465 %! Mo = [Ao, Bo; Co, Do];
 466 %! Me = [Ae, Be; Ce, De];
 467 %!
 468 %!assert (Mo, Me, 1e-4);
 469 %!assert (Info.hsv, HSVe, 1e-4);