octave_packages/control-2.3.52/hinfsyn.m

   1 ## Copyright (C) 2009   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{K}, @var{N}, @var{gamma}, @var{rcond}] =} hinfsyn (@var{P}, @var{nmeas}, @var{ncon})
  20 ## @deftypefnx{Function File} {[@var{K}, @var{N}, @var{gamma}, @var{rcond}] =} hinfsyn (@var{P}, @var{nmeas}, @var{ncon}, @var{gmax})
  21 ## H-infinity control synthesis for LTI plant.
  22 ##
  23 ## @strong{Inputs}
  24 ## @table @var
  25 ## @item P
  26 ## Generalized plant.  Must be a proper/realizable LTI model.
  27 ## @item nmeas
  28 ## Number of measured outputs v.  The last @var{nmeas} outputs of @var{P} are connected to the
  29 ## inputs of controller @var{K}.  The remaining outputs z (indices 1 to p-nmeas) are used
  30 ## to calculate the H-infinity norm.
  31 ## @item ncon
  32 ## Number of controlled inputs u.  The last @var{ncon} inputs of @var{P} are connected to the
  33 ## outputs of controller @var{K}.  The remaining inputs w (indices 1 to m-ncon) are excited
  34 ## by a harmonic test signal.
  35 ## @item gmax
  36 ## The maximum value of the H-infinity norm of @var{N}.  It is assumed that @var{gmax} is
  37 ## sufficiently large so that the controller is admissible.
  38 ## @end table
  39 ##
  40 ## @strong{Outputs}
  41 ## @table @var
  42 ## @item K
  43 ## State-space model of the H-infinity (sub-)optimal controller.
  44 ## @item N
  45 ## State-space model of the lower LFT of @var{P} and @var{K}.
  46 ## @item gamma
  47 ## L-infinity norm of @var{N}.
  48 ## @item rcond
  49 ## Vector @var{rcond} contains estimates of the reciprocal condition
  50 ## numbers of the matrices which are to be inverted and
  51 ## estimates of the reciprocal condition numbers of the
  52 ## Riccati equations which have to be solved during the
  53 ## computation of the controller @var{K}.  For details,
  54 ## see the description of the corresponding SLICOT algorithm.
  55 ## @end table
  56 ##
  57 ## @strong{Block Diagram}
  58 ## @example
  59 ## @group
  60 ##
  61 ## gamma = min||N(K)||             N = lft (P, K)
  62 ##          K         inf
  63 ##
  64 ##                +--------+
  65 ##        w ----->|        |-----> z
  66 ##                |  P(s)  |
  67 ##        u +---->|        |-----+ v
  68 ##          |     +--------+     |
  69 ##          |                    |
  70 ##          |     +--------+     |
  71 ##          +-----|  K(s)  |<----+
  72 ##                +--------+
  73 ##
  74 ##                +--------+
  75 ##        w ----->|  N(s)  |-----> z
  76 ##                +--------+
  77 ## @end group
  78 ## @end example
  79 ##
  80 ## @strong{Algorithm}@*
  81 ## Uses SLICOT SB10FD and SB10DD by courtesy of
  82 ## @uref{http://www.slicot.org, NICONET e.V.}
  83 ##
  84 ## @seealso{augw, mixsyn}
  85 ## @end deftypefn
  86
  87 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  88 ## Created: December 2009
  89 ## Version: 0.1
  90
  91 ## TODO: improve compatibility for nargin >= 4
  92
  93 function [K, varargout] = hinfsyn (P, nmeas, ncon, gmax = 1e15)
  94
  95   ## check input arguments
  96   if (nargin < 3 || nargin > 4)
  97     print_usage ();
  98   endif
  99
 100   if (! isa (P, "lti"))
 101     error ("hinfsyn: first argument must be an LTI system");
 102   endif
 103
 104   if (! is_real_scalar (nmeas))
 105     error ("hinfsyn: second argument invalid");
 106   endif
 107
 108   if (! is_real_scalar (ncon))
 109     error ("hinfsyn: third argument invalid");
 110   endif
 111
 112   if (! is_real_scalar (gmax) || gmax < 0)
 113     error ("hinfsyn: fourth argument invalid");
 114   endif
 115
 116   [a, b, c, d, tsam] = ssdata (P);
 117
 118   ## check assumption A1
 119   m = columns (b);
 120   p = rows (c);
 121
 122   m1 = m - ncon;
 123   p1 = p - nmeas;
 124
 125   if (! isstabilizable (P(:, m1+1:m)))
 126     error ("hinfsyn: (A, B2) must be stabilizable");
 127   endif
 128
 129   if (! isdetectable (P(p1+1:p, :)))
 130     error ("hinfsyn: (C2, A) must be detectable");
 131   endif
 132
 133   ## H-infinity synthesis
 134   if (isct (P))             # continuous plant
 135     [ak, bk, ck, dk, rcond] = slsb10fd (a, b, c, d, ncon, nmeas, gmax);
 136   else                      # discrete plant
 137     [ak, bk, ck, dk, rcond] = slsb10dd (a, b, c, d, ncon, nmeas, gmax);
 138   endif
 139
 140   ## controller
 141   K = ss (ak, bk, ck, dk, tsam);
 142
 143   if (nargout > 1)
 144     N = lft (P, K);
 145     varargout{1} = N;
 146     if (nargout > 2)
 147       varargout{2} = norm (N, inf);
 148       if (nargout > 3)
 149         varargout{3} = rcond;
 150       endif
 151     endif
 152   endif
 153
 154 endfunction
 155
 156
 157 ## continuous-time case
 158 %!shared M, M_exp
 159 %! A = [-1.0  0.0  4.0  5.0 -3.0 -2.0
 160 %!      -2.0  4.0 -7.0 -2.0  0.0  3.0
 161 %!      -6.0  9.0 -5.0  0.0  2.0 -1.0
 162 %!      -8.0  4.0  7.0 -1.0 -3.0  0.0
 163 %!       2.0  5.0  8.0 -9.0  1.0 -4.0
 164 %!       3.0 -5.0  8.0  0.0  2.0 -6.0];
 165 %!
 166 %! B = [-3.0 -4.0 -2.0  1.0  0.0
 167 %!       2.0  0.0  1.0 -5.0  2.0
 168 %!      -5.0 -7.0  0.0  7.0 -2.0
 169 %!       4.0 -6.0  1.0  1.0 -2.0
 170 %!      -3.0  9.0 -8.0  0.0  5.0
 171 %!       1.0 -2.0  3.0 -6.0 -2.0];
 172 %!
 173 %! C = [ 1.0 -1.0  2.0 -4.0  0.0 -3.0
 174 %!      -3.0  0.0  5.0 -1.0  1.0  1.0
 175 %!      -7.0  5.0  0.0 -8.0  2.0 -2.0
 176 %!       9.0 -3.0  4.0  0.0  3.0  7.0
 177 %!       0.0  1.0 -2.0  1.0 -6.0 -2.0];
 178 %!
 179 %! D = [ 1.0 -2.0 -3.0  0.0  0.0
 180 %!       0.0  4.0  0.0  1.0  0.0
 181 %!       5.0 -3.0 -4.0  0.0  1.0
 182 %!       0.0  1.0  0.0  1.0 -3.0
 183 %!       0.0  0.0  1.0  7.0  1.0];
 184 %!
 185 %! P = ss (A, B, C, D);
 186 %! K = hinfsyn (P, 2, 2, 15);
 187 %! M = [K.A, K.B; K.C, K.D];
 188 %!
 189 %! KA = [ -2.8043  14.7367   4.6658   8.1596   0.0848   2.5290
 190 %!         4.6609   3.2756  -3.5754  -2.8941   0.2393   8.2920
 191 %!       -15.3127  23.5592  -7.1229   2.7599   5.9775  -2.0285
 192 %!       -22.0691  16.4758  12.5523 -16.3602   4.4300  -3.3168
 193 %!        30.6789  -3.9026  -1.3868  26.2357  -8.8267  10.4860
 194 %!        -5.7429   0.0577  10.8216 -11.2275   1.5074 -10.7244];
 195 %!
 196 %! KB = [ -0.1581  -0.0793
 197 %!        -0.9237  -0.5718
 198 %!         0.7984   0.6627
 199 %!         0.1145   0.1496
 200 %!        -0.6743  -0.2376
 201 %!         0.0196  -0.7598];
 202 %!
 203 %! KC = [ -0.2480  -0.1713  -0.0880   0.1534   0.5016  -0.0730
 204 %!         2.8810  -0.3658   1.3007   0.3945   1.2244   2.5690];
 205 %!
 206 %! KD = [  0.0554   0.1334
 207 %!        -0.3195   0.0333];
 208 %!
 209 %! M_exp = [KA, KB; KC, KD];
 210 %!
 211 %!assert (M, M_exp, 1e-4);
 212
 213
 214 ## discrete-time case
 215 %!shared M, M_exp
 216 %! A = [-0.7  0.0  0.3  0.0 -0.5 -0.1
 217 %!      -0.6  0.2 -0.4 -0.3  0.0  0.0
 218 %!      -0.5  0.7 -0.1  0.0  0.0 -0.8
 219 %!      -0.7  0.0  0.0 -0.5 -1.0  0.0
 220 %!       0.0  0.3  0.6 -0.9  0.1 -0.4
 221 %!       0.5 -0.8  0.0  0.0  0.2 -0.9];
 222 %!
 223 %! B = [-1.0 -2.0 -2.0  1.0  0.0
 224 %!       1.0  0.0  1.0 -2.0  1.0
 225 %!      -3.0 -4.0  0.0  2.0 -2.0
 226 %!       1.0 -2.0  1.0  0.0 -1.0
 227 %!       0.0  1.0 -2.0  0.0  3.0
 228 %!       1.0  0.0  3.0 -1.0 -2.0];
 229 %!
 230 %! C = [ 1.0 -1.0  2.0 -2.0  0.0 -3.0
 231 %!      -3.0  0.0  1.0 -1.0  1.0  0.0
 232 %!       0.0  2.0  0.0 -4.0  0.0 -2.0
 233 %!       1.0 -3.0  0.0  0.0  3.0  1.0
 234 %!       0.0  1.0 -2.0  1.0  0.0 -2.0];
 235 %!
 236 %! D = [ 1.0 -1.0 -2.0  0.0  0.0
 237 %!       0.0  1.0  0.0  1.0  0.0
 238 %!       2.0 -1.0 -3.0  0.0  1.0
 239 %!       0.0  1.0  0.0  1.0 -1.0
 240 %!       0.0  0.0  1.0  2.0  1.0];
 241 %!
 242 %! P = ss (A, B, C, D, 1);  # value of sampling time doesn't matter
 243 %! K = hinfsyn (P, 2, 2, 111.294);
 244 %! M = [K.A, K.B; K.C, K.D];
 245 %!
 246 %! KA = [-18.0030  52.0376  26.0831  -0.4271 -40.9022  18.0857
 247 %!        18.8203 -57.6244 -29.0938   0.5870  45.3309 -19.8644
 248 %!       -26.5994  77.9693  39.0368  -1.4020 -60.1129  26.6910
 249 %!       -21.4163  62.1719  30.7507  -0.9201 -48.6221  21.8351
 250 %!        -0.8911   4.2787   2.3286  -0.2424  -3.0376   1.2169
 251 %!        -5.3286  16.1955   8.4824  -0.2489 -12.2348   5.1590];
 252 %!
 253 %! KB = [ 16.9788  14.1648
 254 %!       -18.9215 -15.6726
 255 %!        25.2046  21.2848
 256 %!        20.1122  16.8322
 257 %!         1.4104   1.2040
 258 %!         5.3181   4.5149];
 259 %!
 260 %! KC = [ -9.1941  27.5165  13.7364  -0.3639 -21.5983   9.6025
 261 %!         3.6490 -10.6194  -5.2772   0.2432   8.1108  -3.6293];
 262 %!
 263 %! KD = [  9.0317   7.5348
 264 %!        -3.4006  -2.8219];
 265 %!
 266 %! M_exp = [KA, KB; KC, KD];
 267 %!
 268 %!assert (M, M_exp, 1e-4);