octave_packages/control-2.3.52/cfconred.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}] =} cfconred (@var{G}, @var{F}, @var{L}, @dots{})
  20 ## @deftypefnx{Function File} {[@var{Kr}, @var{info}] =} cfconred (@var{G}, @var{F}, @var{L}, @var{ncr}, @dots{})
  21 ## @deftypefnx{Function File} {[@var{Kr}, @var{info}] =} cfconred (@var{G}, @var{F}, @var{L}, @var{opt}, @dots{})
  22 ## @deftypefnx{Function File} {[@var{Kr}, @var{info}] =} cfconred (@var{G}, @var{F}, @var{L}, @var{ncr}, @var{opt}, @dots{})
  23 ##
  24 ## Reduction of state-feedback-observer based controller by coprime factorization (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}.
  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-bta', 'b'
  75 ## Use the square-root Balance & Truncate method.
  76 ## @item 'bfsr-bta', 'f'
  77 ## Use the balancing-free square-root Balance & Truncate method.  Default method.
  78 ## @item 'sr-spa', 's'
  79 ## Use the square-root Singular Perturbation Approximation method.
  80 ## @item 'bfsr-spa', 'p'
  81 ## Use the balancing-free square-root Singular Perturbation Approximation method.
  82 ## @end table
  83 ##
  84 ## @item 'cf'
  85 ## Specifies whether left or right coprime factorization is
  86 ## to be used as follows:
  87 ## @table @var
  88 ## @item 'left', 'l'
  89 ## Use left coprime factorization.  Default method.
  90 ## @item 'right', 'r'
  91 ## Use right coprime factorization.
  92 ## @end table
  93 ##
  94 ## @item 'feedback'
  95 ## Specifies whether @var{F} and @var{L} are fed back positively or negatively:
  96 ## @table @var
  97 ## @item '+'
  98 ## A+BK and A+LC are both Hurwitz matrices.
  99 ## @item '-'
 100 ## A-BK and A-LC are both Hurwitz matrices.  Default value.
 101 ## @end table
 102 ##
 103 ## @item 'tol1'
 104 ## If @var{'order'} is not specified, @var{tol1} contains the tolerance for
 105 ## determining the order of the reduced system.
 106 ## For model reduction, the recommended value of @var{tol1} is
 107 ## c*info.hsv(1), where c lies in the interval [0.00001, 0.001].
 108 ## Default value is n*eps*info.hsv(1).
 109 ## If @var{'order'} is specified, the value of @var{tol1} is ignored.
 110 ##
 111 ## @item 'tol2'
 112 ## The tolerance for determining the order of a minimal
 113 ## realization of the coprime factorization controller.
 114 ## TOL2 <= TOL1.
 115 ## If not specified, n*eps*info.hsv(1) is chosen.
 116 ##
 117 ## @item 'equil', 'scale'
 118 ## Boolean indicating whether equilibration (scaling) should be
 119 ## performed on system @var{G} prior to order reduction.
 120 ## Default value is true if @code{G.scaled == false} and
 121 ## false if @code{G.scaled == true}.
 122 ## Note that for @acronym{MIMO} models, proper scaling of both inputs and outputs
 123 ## is of utmost importance.  The input and output scaling can @strong{not}
 124 ## be done by the equilibration option or the @command{prescale} command
 125 ## because these functions perform state transformations only.
 126 ## Furthermore, signals should not be scaled simply to a certain range.
 127 ## For all inputs (or outputs), a certain change should be of the same
 128 ## importance for the model.
 129 ## @end table
 130 ##
 131 ## @strong{Algorithm}@*
 132 ## Uses SLICOT SB16BD by courtesy of
 133 ## @uref{http://www.slicot.org, NICONET e.V.}
 134 ## @end deftypefn
 135
 136 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
 137 ## Created: December 2011
 138 ## Version: 0.1
 139
 140 function [Kr, info] = cfconred (G, F, L, varargin)
 141
 142   if (nargin < 3)
 143     print_usage ();
 144   endif
 145
 146   if (! isa (G, "lti"))
 147     error ("cfconred: first argument must be an LTI system");
 148   endif
 149
 150   if (! is_real_matrix (F))
 151     error ("cfconred: second argument must be a real matrix");
 152   endif
 153
 154   if (! is_real_matrix (L))
 155     error ("cfconred: third argument must be a real matrix");
 156   endif
 157
 158   if (nargin > 3)                                  # cfconred (G, F, L, ...)
 159     if (is_real_scalar (varargin{1}))              # cfconred (G, F, L, nr)
 160       varargin = horzcat (varargin(2:end), {"order"}, varargin(1));
 161     endif
 162     if (isstruct (varargin{1}))                    # cfconred (G, F, L, opt, ...), cfconred (G, F, L, nr, opt, ...)
 163       varargin = horzcat (__opt2cell__ (varargin{1}), varargin(2:end));
 164     endif
 165     ## order placed at the end such that nr from cfconred (G, F, L, nr, ...)
 166     ## and cfconred (G, F, L, nr, opt, ...) overrides possible nr's from
 167     ## key/value-pairs and inside opt struct (later keys override former keys,
 168     ## nr > key/value > opt)
 169   endif
 170
 171   nkv = numel (varargin);                          # number of keys and values
 172
 173   if (rem (nkv, 2))
 174     error ("cfconred: keys and values must come in pairs");
 175   endif
 176
 177   [a, b, c, d, tsam, scaled] = ssdata (G);
 178   [p, m] = size (G);
 179   n = rows (a);
 180   [mf, nf] = size (F);
 181   [nl, pl] = size (L);
 182   dt = isdt (G);
 183   jobd = any (d(:));
 184
 185   if (mf != m || nf != n)
 186     error ("cfconred: dimensions of state-feedback matrix (%dx%d) and plant (%dx%d, %d states) don't match", \
 187            mf, nf, p, m, n);
 188   endif
 189
 190   if (nl != n || pl != p)
 191     error ("cfconred: dimensions of observer matrix (%dx%d) and plant (%dx%d, %d states) don't match", \
 192            nl, pl, p, m, n);
 193   endif
 194
 195   ## default arguments
 196   tol1 = 0.0;
 197   tol2 = 0.0;
 198   jobcf = 0;
 199   jobmr = 2;                                       # balancing-free BTA
 200   equil = scaled;                                  # equil: 0 means "S", 1 means "N"
 201   ordsel = 1;
 202   ncr = 0;
 203   negfb = true;                                    # A-BK, A-LC Hurwitz
 204
 205
 206   ## handle keys and values
 207   for k = 1 : 2 : nkv
 208     key = lower (varargin{k});
 209     val = varargin{k+1};
 210     switch (key)
 211       case {"order", "ncr", "nr"}
 212         [ncr, ordsel] = __modred_check_order__ (val, n);
 213
 214       case "tol1"
 215         tol1 = __modred_check_tol__ (val, "tol1");
 216
 217       case "tol2"
 218         tol2 = __modred_check_tol__ (val, "tol2");
 219
 220       case "cf"
 221         switch (lower (val(1)))
 222           case "l"
 223             jobcf = 0;
 224           case "r"
 225             jobcf = 1;
 226           otherwise
 227             error ("cfconred: '%s' is an invalid coprime factorization", val);
 228         endswitch
 229
 230       case "method"                                # approximation method
 231         switch (tolower (val))
 232           case {"sr-bta", "b"}                     # 'B':  use the square-root Balance & Truncate method
 233             jobmr = 0;
 234           case {"bfsr-bta", "f"}                   # 'F':  use the balancing-free square-root Balance & Truncate method
 235             jobmr = 1;
 236           case {"sr-spa", "s"}                     # 'S':  use the square-root Singular Perturbation Approximation method
 237             jobmr = 2;
 238           case {"bfsr-spa", "p"}                   # 'P':  use the balancing-free square-root Singular Perturbation Approximation method
 239             jobmr = 3;
 240           otherwise
 241             error ("cfconred: '%s' is an invalid approach", val);
 242         endswitch
 243
 244       case {"equil", "equilibrate", "equilibration", "scale", "scaling"}
 245         equil = __modred_check_equil__ (val);
 246
 247       case "feedback"
 248         negfb = __conred_check_feedback_sign__ (val);
 249
 250       otherwise
 251         warning ("cfconred: invalid property name '%s' ignored", key);
 252     endswitch
 253   endfor
 254
 255
 256   ## A - B*F --> A + B*F  ;    A - L*C --> A + L*C
 257   if (negfb)
 258     F = -F;
 259     L = -L;
 260   endif
 261
 262   ## perform model order reduction
 263   [acr, bcr, ccr, dcr, ncr, hsv] = slsb16bd (a, b, c, d, dt, equil, ncr, ordsel, jobd, jobmr, \
 264                                              F, L, jobcf, tol1, tol2);
 265
 266   ## assemble reduced order controller
 267   Kr = ss (acr, bcr, ccr, dcr, tsam);
 268
 269   ## assemble info struct
 270   info = struct ("ncr", ncr, "hsv", hsv);
 271
 272 endfunction
 273
 274
 275 %!shared Mo, Me, Info, HSVe
 276 %! A =  [       0    1.0000         0         0         0         0         0        0
 277 %!              0         0         0         0         0         0         0        0
 278 %!              0         0   -0.0150    0.7650         0         0         0        0
 279 %!              0         0   -0.7650   -0.0150         0         0         0        0
 280 %!              0         0         0         0   -0.0280    1.4100         0        0
 281 %!              0         0         0         0   -1.4100   -0.0280         0        0
 282 %!              0         0         0         0         0         0   -0.0400    1.850
 283 %!              0         0         0         0         0         0   -1.8500   -0.040 ];
 284 %!
 285 %! B =  [  0.0260
 286 %!        -0.2510
 287 %!         0.0330
 288 %!        -0.8860
 289 %!        -4.0170
 290 %!         0.1450
 291 %!         3.6040
 292 %!         0.2800 ];
 293 %!
 294 %! C =  [  -.996 -.105 0.261 .009 -.001 -.043 0.002 -0.026 ];
 295 %!
 296 %! D =  [  0.0 ];
 297 %!
 298 %! G = ss (A, B, C, D);  % "scaled", false
 299 %!
 300 %! F =  [  4.4721e-002  6.6105e-001  4.6986e-003  3.6014e-001  1.0325e-001 -3.7541e-002 -4.2685e-002  3.2873e-002 ];
 301 %!
 302 %! L =  [  4.1089e-001
 303 %!         8.6846e-002
 304 %!         3.8523e-004
 305 %!        -3.6194e-003
 306 %!        -8.8037e-003
 307 %!         8.4205e-003
 308 %!         1.2349e-003
 309 %!         4.2632e-003 ];
 310 %!
 311 %! [Kr, Info] = cfconred (G, F, L, 4, "method", "bfsr-bta", "cf", "left", "feedback", "+");
 312 %! [Ao, Bo, Co, Do] = ssdata (Kr);
 313 %!
 314 %! Ae = [  0.5946  -0.7336   0.1914  -0.3368
 315 %!         0.5960  -0.0184  -0.1088   0.0207
 316 %!         1.2253   0.2043   0.1009  -1.4948
 317 %!        -0.0330  -0.0243   1.3440   0.0035 ];
 318 %!
 319 %! Be = [  0.0015
 320 %!        -0.0202
 321 %!         0.0159
 322 %!        -0.0544 ];
 323 %!
 324 %! Ce = [  0.3534   0.0274   0.0337  -0.0320 ];
 325 %!
 326 %! De = [  0.0000 ];
 327 %!
 328 %! HSVe = [  4.9078   4.8745   3.8455   3.7811   1.2289   1.1785   0.5176   0.1148 ].';
 329 %!
 330 %! Mo = [Ao, Bo; Co, Do];
 331 %! Me = [Ae, Be; Ce, De];
 332 %!
 333 %!assert (Mo, Me, 1e-4);
 334 %!assert (Info.hsv, HSVe, 1e-4);