octave_packages/control-2.3.52/h2syn.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}] =} h2syn (@var{P}, @var{nmeas}, @var{ncon})
  20 ## H-2 control synthesis for LTI plant.
  21 ##
  22 ## @strong{Inputs}
  23 ## @table @var
  24 ## @item P
  25 ## Generalized plant.  Must be a proper/realizable LTI model.
  26 ## @item nmeas
  27 ## Number of measured outputs v.  The last @var{nmeas} outputs of @var{P} are connected to the
  28 ## inputs of controller @var{K}.  The remaining outputs z (indices 1 to p-nmeas) are used
  29 ## to calculate the H-2 norm.
  30 ## @item ncon
  31 ## Number of controlled inputs u.  The last @var{ncon} inputs of @var{P} are connected to the
  32 ## outputs of controller @var{K}.  The remaining inputs w (indices 1 to m-ncon) are excited
  33 ## by a harmonic test signal.
  34 ## @end table
  35 ##
  36 ## @strong{Outputs}
  37 ## @table @var
  38 ## @item K
  39 ## State-space model of the H-2 optimal controller.
  40 ## @item N
  41 ## State-space model of the lower LFT of @var{P} and @var{K}.
  42 ## @item gamma
  43 ## H-2 norm of @var{N}.
  44 ## @item rcond
  45 ## Vector @var{rcond} contains estimates of the reciprocal condition
  46 ## numbers of the matrices which are to be inverted and
  47 ## estimates of the reciprocal condition numbers of the
  48 ## Riccati equations which have to be solved during the
  49 ## computation of the controller @var{K}.  For details,
  50 ## see the description of the corresponding SLICOT algorithm.
  51 ## @end table
  52 ##
  53 ## @strong{Block Diagram}
  54 ## @example
  55 ## @group
  56 ##
  57 ## gamma = min||N(K)||             N = lft (P, K)
  58 ##          K         2
  59 ##
  60 ##                +--------+
  61 ##        w ----->|        |-----> z
  62 ##                |  P(s)  |
  63 ##        u +---->|        |-----+ v
  64 ##          |     +--------+     |
  65 ##          |                    |
  66 ##          |     +--------+     |
  67 ##          +-----|  K(s)  |<----+
  68 ##                +--------+
  69 ##
  70 ##                +--------+
  71 ##        w ----->|  N(s)  |-----> z
  72 ##                +--------+
  73 ## @end group
  74 ## @end example
  75 ##
  76 ## @strong{Algorithm}@*
  77 ## Uses SLICOT SB10HD and SB10ED by courtesy of
  78 ## @uref{http://www.slicot.org, NICONET e.V.}
  79 ##
  80 ## @seealso{augw, lqr, dlqr, kalman}
  81 ## @end deftypefn
  82
  83 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  84 ## Created: December 2009
  85 ## Version: 0.1
  86
  87 function [K, varargout] = h2syn (P, nmeas, ncon)
  88
  89   ## check input arguments
  90   if (nargin != 3)
  91     print_usage ();
  92   endif
  93
  94   if (! isa (P, "lti"))
  95     error ("h2syn: first argument must be an LTI system");
  96   endif
  97
  98   if (! is_real_scalar (nmeas))
  99     error ("h2syn: second argument invalid");
 100   endif
 101
 102   if (! is_real_scalar (ncon))
 103     error ("h2syn: third argument invalid");
 104   endif
 105
 106   [a, b, c, d, tsam] = ssdata (P);
 107
 108   ## check assumptions A1 - A3
 109   m = columns (b);
 110   p = rows (c);
 111
 112   m1 = m - ncon;
 113   p1 = p - nmeas;
 114
 115   d11 = d(1:p1, 1:m1);
 116
 117   if (isct (P) && any (d11(:)))
 118     warning ("h2syn: setting matrice D11 to zero");
 119     d(1:p1, 1:m1) = 0;
 120   endif
 121
 122   if (! isstabilizable (P(:, m1+1:m)))
 123     error ("h2syn: (A, B2) must be stabilizable");
 124   endif
 125
 126   if (! isdetectable (P(p1+1:p, :)))
 127     error ("h2syn: (C2, A) must be detectable");
 128   endif
 129
 130   ## H-2 synthesis
 131   if (isct (P))             # continuous plant
 132     [ak, bk, ck, dk, rcond] = slsb10hd (a, b, c, d, ncon, nmeas);
 133   else                      # discrete plant
 134     [ak, bk, ck, dk, rcond] = slsb10ed (a, b, c, d, ncon, nmeas);
 135   endif
 136
 137   ## controller
 138   K = ss (ak, bk, ck, dk, tsam);
 139
 140   if (nargout > 1)
 141     N = lft (P, K);
 142     varargout{1} = N;
 143     if (nargout > 2)
 144       varargout{2} = norm (N, 2);
 145       if (nargout > 3)
 146         varargout{3} = rcond;
 147       endif
 148     endif
 149   endif
 150
 151 endfunction
 152
 153
 154 ## continuous-time case
 155 %!shared M, M_exp
 156 %! A = [-1.0  0.0  4.0  5.0 -3.0 -2.0
 157 %!      -2.0  4.0 -7.0 -2.0  0.0  3.0
 158 %!      -6.0  9.0 -5.0  0.0  2.0 -1.0
 159 %!      -8.0  4.0  7.0 -1.0 -3.0  0.0
 160 %!       2.0  5.0  8.0 -9.0  1.0 -4.0
 161 %!       3.0 -5.0  8.0  0.0  2.0 -6.0];
 162 %!
 163 %! B = [-3.0 -4.0 -2.0  1.0  0.0
 164 %!       2.0  0.0  1.0 -5.0  2.0
 165 %!      -5.0 -7.0  0.0  7.0 -2.0
 166 %!       4.0 -6.0  1.0  1.0 -2.0
 167 %!      -3.0  9.0 -8.0  0.0  5.0
 168 %!       1.0 -2.0  3.0 -6.0 -2.0];
 169 %!
 170 %! C = [ 1.0 -1.0  2.0 -4.0  0.0 -3.0
 171 %!      -3.0  0.0  5.0 -1.0  1.0  1.0
 172 %!      -7.0  5.0  0.0 -8.0  2.0 -2.0
 173 %!       9.0 -3.0  4.0  0.0  3.0  7.0
 174 %!       0.0  1.0 -2.0  1.0 -6.0 -2.0];
 175 %!
 176 %! D = [ 0.0  0.0  0.0 -4.0 -1.0
 177 %!       0.0  0.0  0.0  1.0  0.0
 178 %!       0.0  0.0  0.0  0.0  1.0
 179 %!       3.0  1.0  0.0  1.0 -3.0
 180 %!      -2.0  0.0  1.0  7.0  1.0];
 181 %!
 182 %! P = ss (A, B, C, D);
 183 %! K = h2syn (P, 2, 2);
 184 %! M = [K.A, K.B; K.C, K.D];
 185 %!
 186 %! KA = [  88.0015  -145.7298   -46.2424    82.2168   -45.2996   -31.1407
 187 %!         25.7489   -31.4642   -12.4198     9.4625    -3.5182     2.7056
 188 %!         54.3008  -102.4013   -41.4968    50.8412   -20.1286   -26.7191
 189 %!        108.1006  -198.0785   -45.4333    70.3962   -25.8591   -37.2741
 190 %!       -115.8900   226.1843    47.2549   -47.8435   -12.5004    34.7474
 191 %!         59.0362  -101.8471   -20.1052    36.7834   -16.1063   -26.4309];
 192 %!
 193 %! KB = [   3.7345     3.4758
 194 %!         -0.3020     0.6530
 195 %!          3.4735     4.0499
 196 %!          4.3198     7.2755
 197 %!         -3.9424   -10.5942
 198 %!          2.1784     2.5048];
 199 %!
 200 %! KC = [  -2.3346     3.2556     0.7150    -0.9724     0.6962     0.4074
 201 %!          7.6899    -8.4558    -2.9642     7.0365    -4.2844     0.1390];
 202 %!
 203 %! KD = [   0.0000     0.0000
 204 %!          0.0000     0.0000];
 205 %!
 206 %! M_exp = [KA, KB; KC, KD];
 207 %!
 208 %!assert (M, M_exp, 1e-4);
 209
 210
 211 ## discrete-time case
 212 %!shared M, M_exp
 213 %! A = [-0.7  0.0  0.3  0.0 -0.5 -0.1
 214 %!      -0.6  0.2 -0.4 -0.3  0.0  0.0
 215 %!      -0.5  0.7 -0.1  0.0  0.0 -0.8
 216 %!      -0.7  0.0  0.0 -0.5 -1.0  0.0
 217 %!       0.0  0.3  0.6 -0.9  0.1 -0.4
 218 %!       0.5 -0.8  0.0  0.0  0.2 -0.9];
 219 %!
 220 %! B = [-1.0 -2.0 -2.0  1.0  0.0
 221 %!       1.0  0.0  1.0 -2.0  1.0
 222 %!      -3.0 -4.0  0.0  2.0 -2.0
 223 %!       1.0 -2.0  1.0  0.0 -1.0
 224 %!       0.0  1.0 -2.0  0.0  3.0
 225 %!       1.0  0.0  3.0 -1.0 -2.0];
 226 %!
 227 %! C = [ 1.0 -1.0  2.0 -2.0  0.0 -3.0
 228 %!      -3.0  0.0  1.0 -1.0  1.0  0.0
 229 %!       0.0  2.0  0.0 -4.0  0.0 -2.0
 230 %!       1.0 -3.0  0.0  0.0  3.0  1.0
 231 %!       0.0  1.0 -2.0  1.0  0.0 -2.0];
 232 %!
 233 %! D = [ 1.0 -1.0 -2.0  0.0  0.0
 234 %!       0.0  1.0  0.0  1.0  0.0
 235 %!       2.0 -1.0 -3.0  0.0  1.0
 236 %!       0.0  1.0  0.0  1.0 -1.0
 237 %!       0.0  0.0  1.0  2.0  1.0];
 238 %!
 239 %! P = ss (A, B, C, D, 1);  # value of sampling time doesn't matter
 240 %! K = h2syn (P, 2, 2);
 241 %! M = [K.A, K.B; K.C, K.D];
 242 %!
 243 %! KA = [-0.0551  -2.1891  -0.6607  -0.2532   0.6674  -1.0044
 244 %!       -1.0379   2.3804   0.5031   0.3960  -0.6605   1.2673
 245 %!       -0.0876  -2.1320  -0.4701  -1.1461   1.2927  -1.5116
 246 %!       -0.1358  -2.1237  -0.9560  -0.7144   0.6673  -0.7957
 247 %!        0.4900   0.0895   0.2634  -0.2354   0.1623  -0.2663
 248 %!        0.1672  -0.4163   0.2871  -0.1983   0.4944  -0.6967];
 249 %!
 250 %! KB = [-0.5985  -0.5464
 251 %!        0.5285   0.6087
 252 %!       -0.7600  -0.4472
 253 %!       -0.7288  -0.6090
 254 %!        0.0532   0.0658
 255 %!       -0.0663   0.0059];
 256 %!
 257 %! KC = [ 0.2500  -1.0200  -0.3371  -0.2733   0.2747  -0.4444
 258 %!        0.0654   0.2095   0.0632   0.2089  -0.1895   0.1834];
 259 %!
 260 %! KD = [-0.2181  -0.2070
 261 %!        0.1094   0.1159];
 262 %!
 263 %! M_exp = [KA, KB; KC, KD];
 264 %!
 265 %!assert (M, M_exp, 1e-4);