octave_packages/control-2.3.52/__modred_ab09id__.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}] =} __modred_ab09id__ (@var{method}, @dots{})
  20 ## Backend for btamodred and spamodred.
  21 ## @end deftypefn
  22
  23 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  24 ## Created: November 2011
  25 ## Version: 0.1
  26
  27 function [Gr, info] = __modred_ab09id__ (method, varargin)
  28
  29   if (nargin < 2)
  30     print_usage ();
  31   endif
  32
  33   if (method != "bta" && method != "spa")
  34     error ("modred: invalid method");
  35   endif
  36
  37   G = varargin{1};
  38   varargin = varargin(2:end);
  39
  40   if (! isa (G, "lti"))
  41     error ("%smodred: first argument must be an LTI system", method);
  42   endif
  43
  44   if (nargin > 2)                                  # *modred (G, ...)
  45     if (is_real_scalar (varargin{1}))              # *modred (G, nr)
  46       varargin = horzcat (varargin(2:end), {"order"}, varargin(1));
  47     endif
  48     if (isstruct (varargin{1}))                    # *modred (G, opt, ...), *modred (G, nr, opt, ...)
  49       varargin = horzcat (__opt2cell__ (varargin{1}), varargin(2:end));
  50     endif
  51     ## order placed at the end such that nr from *modred (G, nr, ...)
  52     ## and *modred (G, nr, opt, ...) overrides possible nr's from
  53     ## key/value-pairs and inside opt struct (later keys override former keys,
  54     ## nr > key/value > opt)
  55   endif
  56
  57   nkv = numel (varargin);                          # number of keys and values
  58
  59   if (rem (nkv, 2))
  60     error ("%smodred: keys and values must come in pairs", method);
  61   endif
  62
  63   [a, b, c, d, tsam, scaled] = ssdata (G);
  64   [p, m] = size (G);
  65   dt = isdt (G);
  66
  67   ## default arguments
  68   alpha = __modred_default_alpha__ (dt);
  69   av = bv = cv = dv = [];
  70   jobv = 0;
  71   aw = bw = cw = dw = [];
  72   jobw = 0;
  73   alphac = alphao = 0.0;
  74   tol1 = 0.0;
  75   tol2 = 0.0;
  76   jobc = jobo = 0;
  77   bf = true;                                # balancing-free
  78   weight = 1;
  79   equil = 0;
  80   ordsel = 1;
  81   nr = 0;
  82
  83   ## handle keys and values
  84   for k = 1 : 2 : nkv
  85     key = lower (varargin{k});
  86     val = varargin{k+1};
  87     switch (key)
  88       case {"left", "output", "v"}
  89         [av, bv, cv, dv, jobv] = __modred_check_weight__ (val, dt, p, []);
  90
  91       case {"right", "input", "w"}
  92         [aw, bw, cw, dw, jobw] = __modred_check_weight__ (val, dt, [], m);
  93
  94       case {"order", "n", "nr"}
  95         [nr, ordsel] = __modred_check_order__ (val, rows (a));
  96
  97       case "tol1"
  98         tol1 = __modred_check_tol__ (val, "tol1");
  99
 100       case "tol2"
 101         tol2 = __modred_check_tol__ (val, "tol2");
 102
 103       case "alpha"
 104         alpha = __modred_check_alpha__ (val, dt);
 105
 106       case "method"
 107         switch (tolower (val))
 108           case "sr"
 109             bf = false;
 110           case "bfsr"
 111             bf = true;
 112           otherwise
 113             error ("modred: '%s' is an invalid approach", val);
 114         endswitch
 115
 116       case {"jobc", "gram-ctrb"}
 117         jobc = __modred_check_gram__ (val, "gram-ctrb");
 118
 119       case {"jobo", "gram-obsv"}
 120         jobo = __modred_check_gram__ (val, "gram-obsv");
 121
 122       case {"alphac", "alpha-ctrb"}
 123         alphac = __modred_check_alpha_gram__ (val, "alpha-ctrb");
 124
 125       case {"alphao", "alpha-obsv"}
 126         alphao = __modred_check_alpha_gram__ (val, "alpha-obsv");
 127
 128       case {"equil", "equilibrate", "equilibration", "scale", "scaling"}
 129         scaled = __modred_check_equil__ (val);
 130
 131       otherwise
 132         warning ("%smodred: invalid property name '%s' ignored", method, key);
 133     endswitch
 134   endfor
 135
 136   ## handle type of frequency weighting
 137   if (jobv && jobw)
 138     weight = 3;                             # 'B':  both left and right weightings V and W are used
 139   elseif (jobv)
 140     weight = 1;                             # 'L':  only left weighting V is used (W = I)
 141   elseif (jobw)
 142     weight = 2;                             # 'R':  only right weighting W is used (V = I)
 143   else
 144     weight = 0;                             # 'N':  no weightings are used (V = I, W = I)
 145   endif
 146
 147   ## handle model reduction approach
 148   if (method == "bta" && ! bf)              # 'B':  use the square-root Balance & Truncate method
 149     job = 0;
 150   elseif (method == "bta" && bf)            # 'F':  use the balancing-free square-root Balance & Truncate method
 151     job = 1;
 152   elseif (method == "spa" && ! bf)          # 'S':  use the square-root Singular Perturbation Approximation method
 153     job = 2;
 154   elseif (method == "spa" && bf)            # 'P':  use the balancing-free square-root Singular Perturbation Approximation method
 155     job = 3;
 156   else
 157     error ("modred: invalid job option");   # this should never happen
 158   endif
 159
 160
 161   ## perform model order reduction
 162   [ar, br, cr, dr, nr, hsv, ns] = slab09id (a, b, c, d, dt, equil, nr, ordsel, alpha, job, \
 163                                             av, bv, cv, dv, \
 164                                             aw, bw, cw, dw, \
 165                                             weight, jobc, jobo, alphac, alphao, \
 166                                             tol1, tol2);
 167
 168   ## assemble reduced order model
 169   Gr = ss (ar, br, cr, dr, tsam);
 170
 171   ## assemble info struct
 172   info = struct ("nr", nr, "ns", ns, "hsv", hsv);
 173
 174 endfunction
 175
 176
 177
 178