octave_packages/control-2.3.52/fwcfconred.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{Kr}, @var{info}] =} fwcfconred (@var{G}, @var{F}, @var{L}, @dots{})
  20 ## @deftypefnx{Function File} {[@var{Kr}, @var{info}] =} fwcfconred (@var{G}, @var{F}, @var{L}, @var{ncr}, @dots{})
  21 ## @deftypefnx{Function File} {[@var{Kr}, @var{info}] =} fwcfconred (@var{G}, @var{F}, @var{L}, @var{opt}, @dots{})
  22 ## @deftypefnx{Function File} {[@var{Kr}, @var{info}] =} fwcfconred (@var{G}, @var{F}, @var{L}, @var{ncr}, @var{opt}, @dots{})
  23 ##
  24 ## Reduction of state-feedback-observer based controller by frequency-weighted coprime factorization (FW CF).
  25 ## Given a plant @var{G}, state feedback gain @var{F} and full observer gain @var{L},
  26 ## determine a reduced order controller @var{Kr} by using stability enforcing frequency weights.
  27 ##
  28 ## @strong{Inputs}
  29 ## @table @var
  30 ## @item G
  31 ## LTI model of the open-loop plant (A,B,C,D).
  32 ## It has m inputs, p outputs and n states.
  33 ## @item F
  34 ## Stabilizing state feedback matrix (m-by-n).
  35 ## @item L
  36 ## Stabilizing observer gain matrix (n-by-p).
  37 ## @item ncr
  38 ## The desired order of the resulting reduced order controller @var{Kr}.
  39 ## If not specified, @var{ncr} is chosen automatically according
  40 ## to the description of key @var{'order'}.
  41 ## @item @dots{}
  42 ## Optional pairs of keys and values.  @code{"key1", value1, "key2", value2}.
  43 ## @item opt
  44 ## Optional struct with keys as field names.
  45 ## Struct @var{opt} can be created directly or
  46 ## by command @command{options}.  @code{opt.key1 = value1, opt.key2 = value2}.
  47 ## @end table
  48 ##
  49 ## @strong{Outputs}
  50 ## @table @var
  51 ## @item Kr
  52 ## State-space model of reduced order controller.
  53 ## @item info
  54 ## Struct containing additional information.
  55 ## @table @var
  56 ## @item info.hsv
  57 ## The Hankel singular values of the extended system?!?.
  58 ## The @var{n} Hankel singular values are ordered decreasingly.
  59 ## @item info.ncr
  60 ## The order of the obtained reduced order controller @var{Kr}.
  61 ## @end table
  62 ## @end table
  63 ##
  64 ## @strong{Option Keys and Values}
  65 ## @table @var
  66 ## @item 'order', 'ncr'
  67 ## The desired order of the resulting reduced order controller @var{Kr}.
  68 ## If not specified, @var{ncr} is chosen automatically such that states with
  69 ## Hankel singular values @var{info.hsv} > @var{tol1} are retained.
  70 ##
  71 ## @item 'method'
  72 ## Order reduction approach to be used as follows:
  73 ## @table @var
  74 ## @item 'sr', 'b'
  75 ## Use the square-root Balance & Truncate method.
  76 ## @item 'bfsr', 'f'
  77 ## Use the balancing-free square-root Balance & Truncate method.  Default method.
  78 ## @end table
  79 ##
  80 ## @item 'cf'
  81 ## Specifies whether left or right coprime factorization is
  82 ## to be used as follows:
  83 ## @table @var
  84 ## @item 'left', 'l'
  85 ## Use left coprime factorization.
  86 ## @item 'right', 'r'
  87 ## Use right coprime factorization.  Default method.
  88 ## @end table
  89 ##
  90 ## @item 'feedback'
  91 ## Specifies whether @var{F} and @var{L} are fed back positively or negatively:
  92 ## @table @var
  93 ## @item '+'
  94 ## A+BK and A+LC are both Hurwitz matrices.
  95 ## @item '-'
  96 ## A-BK and A-LC are both Hurwitz matrices.  Default value.
  97 ## @end table
  98 ##
  99 ## @item 'tol1'
 100 ## If @var{'order'} is not specified, @var{tol1} contains the tolerance for
 101 ## determining the order of the reduced system.
 102 ## For model reduction, the recommended value of @var{tol1} is
 103 ## c*info.hsv(1), where c lies in the interval [0.00001, 0.001].
 104 ## Default value is n*eps*info.hsv(1).
 105 ## If @var{'order'} is specified, the value of @var{tol1} is ignored.
 106 ## @end table
 107 ##
 108 ## @strong{Algorithm}@*
 109 ## Uses SLICOT SB16CD by courtesy of
 110 ## @uref{http://www.slicot.org, NICONET e.V.}
 111 ## @end deftypefn
 112
 113 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
 114 ## Created: December 2011
 115 ## Version: 0.1
 116
 117 function [Kr, info] = fwcfconred (G, F, L, varargin)
 118
 119   if (nargin < 3)
 120     print_usage ();
 121   endif
 122
 123   if (! isa (G, "lti"))
 124     error ("fwcfconred: first argument must be an LTI system");
 125   endif
 126
 127   if (! is_real_matrix (F))
 128     error ("fwcfconred: second argument must be a real matrix");
 129   endif
 130
 131   if (! is_real_matrix (L))
 132     error ("fwcfconred: third argument must be a real matrix");
 133   endif
 134
 135   if (nargin > 3)                                  # fwcfconred (G, F, L, ...)
 136     if (is_real_scalar (varargin{1}))              # fwcfconred (G, F, L, nr)
 137       varargin = horzcat (varargin(2:end), {"order"}, varargin(1));
 138     endif
 139     if (isstruct (varargin{1}))                    # fwcfconred (G, F, L, opt, ...), fwcfconred (G, F, L, nr, opt, ...)
 140       varargin = horzcat (__opt2cell__ (varargin{1}), varargin(2:end));
 141     endif
 142     ## order placed at the end such that nr from fwcfconred (G, F, L, nr, ...)
 143     ## and fwcfconred (G, F, L, nr, opt, ...) overrides possible nr's from
 144     ## key/value-pairs and inside opt struct (later keys override former keys,
 145     ## nr > key/value > opt)
 146   endif
 147
 148   nkv = numel (varargin);                          # number of keys and values
 149
 150   if (rem (nkv, 2))
 151     error ("fwcfconred: keys and values must come in pairs");
 152   endif
 153
 154   [a, b, c, d, tsam] = ssdata (G);
 155   [p, m] = size (G);
 156   n = rows (a);
 157   [mf, nf] = size (F);
 158   [nl, pl] = size (L);
 159   dt = isdt (G);
 160   jobd = any (d(:));
 161
 162   if (mf != m || nf != n)
 163     error ("fwcfconred: dimensions of state-feedback matrix (%dx%d) and plant (%dx%d, %d states) don't match", \
 164            mf, nf, p, m, n);
 165   endif
 166
 167   if (nl != n || pl != p)
 168     error ("fwcfconred: dimensions of observer matrix (%dx%d) and plant (%dx%d, %d states) don't match", \
 169            nl, pl, p, m, n);
 170   endif
 171
 172   ## default arguments
 173   tol1 = 0.0;
 174   jobcf = 1;
 175   jobmr = 1;                                       # balancing-free BTA
 176   ordsel = 1;
 177   ncr = 0;
 178   negfb = true;                                    # A-BK, A-LC Hurwitz
 179
 180
 181   ## handle keys and values
 182   for k = 1 : 2 : nkv
 183     key = lower (varargin{k});
 184     val = varargin{k+1};
 185     switch (key)
 186       case {"order", "ncr", "nr"}
 187         [ncr, ordsel] = __modred_check_order__ (val, n);
 188
 189       case {"tol1", "tol"}
 190         tol1 = __modred_check_tol__ (val, "tol1");
 191
 192       case "cf"
 193         switch (lower (val(1)))
 194           case "l"
 195             jobcf = 0;
 196           case "r"
 197             jobcf = 1;
 198           otherwise
 199             error ("cfconred: '%s' is an invalid coprime factorization", val);
 200         endswitch
 201
 202       case "method"                                # approximation method
 203         switch (tolower (val))
 204           case {"sr-bta", "b", "sr"}               # 'B':  use the square-root Balance & Truncate method
 205             jobmr = 0;
 206           case {"bfsr-bta", "f", "bfsr"}           # 'F':  use the balancing-free square-root Balance & Truncate method
 207             jobmr = 1;
 208           otherwise
 209             error ("cfconred: '%s' is an invalid approach", val);
 210         endswitch
 211
 212       case "feedback"
 213         negfb = __conred_check_feedback_sign__ (val);
 214
 215       otherwise
 216         warning ("fwcfconred: invalid property name '%s' ignored", key);
 217     endswitch
 218   endfor
 219
 220
 221   ## A - B*F --> A + B*F  ;    A - L*C --> A + L*C
 222   if (negfb)
 223     F = -F;
 224     L = -L;
 225   endif
 226
 227   ## perform model order reduction
 228   [acr, bcr, ccr, ncr, hsv] = slsb16cd (a, b, c, d, dt, ncr, ordsel, jobd, jobmr, \
 229                                         F, L, jobcf, tol1);
 230
 231   ## assemble reduced order controller
 232   Kr = ss (acr, bcr, ccr, [], tsam);
 233
 234   ## assemble info struct
 235   info = struct ("ncr", ncr, "hsv", hsv);
 236
 237 endfunction
 238
 239
 240 %!shared Mo, Me, Info, HSVe
 241 %! A =  [       0    1.0000         0         0         0         0         0        0
 242 %!              0         0         0         0         0         0         0        0
 243 %!              0         0   -0.0150    0.7650         0         0         0        0
 244 %!              0         0   -0.7650   -0.0150         0         0         0        0
 245 %!              0         0         0         0   -0.0280    1.4100         0        0
 246 %!              0         0         0         0   -1.4100   -0.0280         0        0
 247 %!              0         0         0         0         0         0   -0.0400    1.850
 248 %!              0         0         0         0         0         0   -1.8500   -0.040 ];
 249 %!
 250 %! B =  [  0.0260
 251 %!        -0.2510
 252 %!         0.0330
 253 %!        -0.8860
 254 %!        -4.0170
 255 %!         0.1450
 256 %!         3.6040
 257 %!         0.2800 ];
 258 %!
 259 %! C =  [  -.996 -.105 0.261 .009 -.001 -.043 0.002 -0.026 ];
 260 %!
 261 %! D =  [  0.0 ];
 262 %!
 263 %! G = ss (A, B, C, D);  % "scaled", false
 264 %!
 265 %! F =  [  4.472135954999638e-002    6.610515358414598e-001    4.698598960657579e-003  3.601363251422058e-001    1.032530880771415e-001   -3.754055214487997e-002  -4.268536964759344e-002    3.287284547842979e-002 ];
 266 %!
 267 %! L =  [  4.108939884667451e-001
 268 %!         8.684600000000012e-002
 269 %!         3.852317308197148e-004
 270 %!        -3.619366874815911e-003
 271 %!        -8.803722876359955e-003
 272 %!         8.420521094001852e-003
 273 %!         1.234944428038507e-003
 274 %!         4.263205617645322e-003 ];
 275 %!
 276 %! [Kr, Info] = fwcfconred (G, F, L, 2, "method", "bfsr", "cf", "right", "feedback", "+");
 277 %! [Ao, Bo, Co, Do] = ssdata (Kr);
 278 %!
 279 %! Ae = [  -0.4334   0.4884
 280 %!         -0.1950  -0.1093 ];
 281 %!
 282 %! Be = [  -0.4231
 283 %!         -0.1785 ];
 284 %!
 285 %! Ce = [  -0.0326  -0.2307 ];
 286 %!
 287 %! De = [  0.0000 ];
 288 %!
 289 %! HSVe = [  3.3073   0.7274   0.1124   0.0784   0.0242   0.0182   0.0101   0.0094 ].';
 290 %!
 291 %! Mo = [Ao, Bo; Co, Do];
 292 %! Me = [Ae, Be; Ce, De];
 293 %!
 294 %!assert (Mo, Me, 1e-4);
 295 %!assert (Info.hsv, HSVe, 1e-4);