octave_packages/control-2.3.52/@lti/feedback.m

   1 ## Copyright (C) 2009, 2010, 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{sys} =} feedback (@var{sys1})
  20 ## @deftypefnx {Function File} {@var{sys} =} feedback (@var{sys1}, @var{"+"})
  21 ## @deftypefnx {Function File} {@var{sys} =} feedback (@var{sys1}, @var{sys2})
  22 ## @deftypefnx {Function File} {@var{sys} =} feedback (@var{sys1}, @var{sys2}, @var{"+"})
  23 ## @deftypefnx {Function File} {@var{sys} =} feedback (@var{sys1}, @var{sys2}, @var{feedin}, @var{feedout})
  24 ## @deftypefnx {Function File} {@var{sys} =} feedback (@var{sys1}, @var{sys2}, @var{feedin}, @var{feedout}, @var{"+"})
  25 ## Feedback connection of two LTI models.
  26 ##
  27 ## @strong{Inputs}
  28 ## @table @var
  29 ## @item sys1
  30 ## LTI model of forward transmission.  @code{[p1, m1] = size (sys1)}.
  31 ## @item sys2
  32 ## LTI model of backward transmission.
  33 ## If not specified, an identity matrix of appropriate size is taken.
  34 ## @item feedin
  35 ## Vector containing indices of inputs to @var{sys1} which are involved in the feedback loop.
  36 ## The number of @var{feedin} indices and outputs of @var{sys2} must be equal.
  37 ## If not specified, @code{1:m1} is taken.
  38 ## @item feedout
  39 ## Vector containing indices of outputs from @var{sys1} which are to be connected to @var{sys2}.
  40 ## The number of @var{feedout} indices and inputs of @var{sys2} must be equal.
  41 ## If not specified, @code{1:p1} is taken.
  42 ## @item "+"
  43 ## Positive feedback sign.  If not specified, @var{"-"} for a negative feedback interconnection
  44 ## is assumed.  @var{+1} and @var{-1} are possible as well, but only from the third argument
  45 ## onward due to ambiguity.
  46 ## @end table
  47 ##
  48 ## @strong{Outputs}
  49 ## @table @var
  50 ## @item sys
  51 ## Resulting LTI model.
  52 ## @end table
  53 ##
  54 ## @strong{Block Diagram}
  55 ## @example
  56 ## @group
  57 ##  u    +         +--------+             y
  58 ## ------>(+)----->|  sys1  |-------+------->
  59 ##         ^ -     +--------+       |
  60 ##         |                        |
  61 ##         |       +--------+       |
  62 ##         +-------|  sys2  |<------+
  63 ##                 +--------+
  64 ## @end group
  65 ## @end example
  66 ## @end deftypefn
  67
  68 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  69 ## Created: October 2009
  70 ## Version: 0.5
  71
  72 function sys = feedback (sys1, sys2, feedin, feedout, fbsign = -1)
  73
  74   [p1, m1] = size (sys1);
  75
  76   switch (nargin)
  77     case 1                          # sys = feedback (sys)
  78       if (p1 != m1)
  79         error ("feedback: argument must be a square system");
  80       endif
  81       sys2 = eye (p1);
  82       feedin = 1 : m1;
  83       feedout = 1 : p1;
  84
  85     case 2
  86       if (ischar (sys2))            # sys = feedback (sys, "+")
  87         if (p1 != m1)
  88           error ("feedback: first argument must be a square system");
  89         endif
  90         fbsign = __check_fbsign__ (sys2);
  91         sys2 = eye (p1);
  92       endif                         # sys = feedback (sys1, sys2)
  93       feedin = 1 : m1;
  94       feedout = 1 : p1;
  95
  96     case 3                          # sys = feedback (sys1, sys2, "+")
  97       fbsign = __check_fbsign__ (feedin);
  98       feedin = 1 : m1;
  99       feedout = 1 : p1;
 100
 101     case 4                          # sys = feedback (sys1, sys2, feedin, feedout)
 102       ## nothing needs to be done here
 103       ## case 4 required to prevent "otherwise"
 104
 105     case 5                          # sys = feedback (sys1, sys2, feedin, feedout, "+")
 106       fbsign = __check_fbsign__ (fbsign);
 107
 108     otherwise
 109       print_usage ();
 110
 111   endswitch
 112
 113   [p2, m2] = size (sys2);
 114
 115   l_feedin = length (feedin);
 116   l_feedout = length (feedout);
 117
 118   if (l_feedin != p2)
 119     error ("feedback: feedin indices: %d, outputs sys2: %d", l_feedin, p2);
 120   endif
 121
 122   if (l_feedout != m2)
 123     error ("feedback: feedout indices: %d, inputs sys2: %d", l_feedout, m2);
 124   endif
 125
 126   if (any (feedin > m1 | feedin < 1))
 127     error ("feedback: range of feedin indices exceeds dimensions of sys1");
 128   endif
 129
 130   if (any (feedin > p1 | feedin < 1))
 131     error ("feedback: range of feedout indices exceeds dimensions of sys1");
 132   endif
 133
 134   M11 = zeros (m1, p1);
 135   M22 = zeros (m2, p2);
 136
 137   M12 = full (sparse (feedin, 1:l_feedin, fbsign, m1, p2));
 138   M21 = full (sparse (1:l_feedout, feedout, 1, m2, p1));
 139
 140   ## NOTE: for-loops do NOT the same as
 141   ##       M12(feedin, 1:l_feedin) = fbsign;
 142   ##       M21(1:l_feedout, feedout) = 1;
 143   ##
 144   ## M12 = zeros (m1, p2);
 145   ## M21 = zeros (m2, p1);
 146   ##
 147   ## for k = 1 : l_feedin
 148   ##   M12(feedin(k), k) = fbsign;
 149   ## endfor
 150   ##
 151   ## for k = 1 : l_feedout
 152   ##   M21(k, feedout(k)) = 1;
 153   ## endfor
 154
 155   M = [M11, M12;
 156        M21, M22];
 157
 158   in_idx = 1 : m1;
 159   out_idx = 1 : p1;
 160
 161   sys = __sys_group__ (sys1, sys2);
 162   sys = __sys_connect__ (sys, M);
 163   sys = __sys_prune__ (sys, out_idx, in_idx);
 164
 165 endfunction
 166
 167
 168 function fbsign = __check_fbsign__ (fbsign)
 169
 170   if (is_real_scalar (fbsign))
 171     fbsign = sign (fbsign);
 172   elseif (ischar (fbsign))
 173     if (strcmp (fbsign, "+"))
 174       fbsign = +1;
 175     elseif (strcmp (fbsign, "-"))
 176       fbsign = -1;
 177     else
 178       error ("feedback: invalid feedback sign string");
 179     endif
 180   else
 181     error ("feedback: invalid feedback sign type");
 182   endif
 183
 184 endfunction
 185
 186
 187 ## Feedback inter-connection of two systems in state-space form
 188 ## Test from SLICOT AB05ND
 189 %!shared M, Me
 190 %! A1 = [ 1.0   0.0  -1.0
 191 %!        0.0  -1.0   1.0
 192 %!        1.0   1.0   2.0 ];
 193 %!
 194 %! B1 = [ 1.0   1.0   0.0
 195 %!        2.0   0.0   1.0 ].';
 196 %!
 197 %! C1 = [ 3.0  -2.0   1.0
 198 %!        0.0   1.0   0.0 ];
 199 %!
 200 %! D1 = [ 1.0   0.0
 201 %!        0.0   1.0 ];
 202 %!
 203 %! A2 = [-3.0   0.0   0.0
 204 %!        1.0   0.0   1.0
 205 %!        0.0  -1.0   2.0 ];
 206 %!
 207 %! B2 = [ 0.0  -1.0   0.0
 208 %!        1.0   0.0   2.0 ].';
 209 %!
 210 %! C2 = [ 1.0   1.0   0.0
 211 %!        1.0   1.0  -1.0 ];
 212 %!
 213 %! D2 = [ 1.0   1.0
 214 %!        0.0   1.0 ];
 215 %!
 216 %! sys1 = ss (A1, B1, C1, D1);
 217 %! sys2 = ss (A2, B2, C2, D2);
 218 %! sys = feedback (sys1, sys2);
 219 %! [A, B, C, D] = ssdata (sys);
 220 %! M = [A, B; C, D];
 221 %!
 222 %! Ae = [-0.5000  -0.2500  -1.5000  -1.2500  -1.2500   0.7500
 223 %!       -1.5000  -0.2500   0.5000  -0.2500  -0.2500  -0.2500
 224 %!        1.0000   0.5000   2.0000  -0.5000  -0.5000   0.5000
 225 %!        0.0000   0.5000   0.0000  -3.5000  -0.5000   0.5000
 226 %!       -1.5000   1.2500  -0.5000   1.2500   0.2500   1.2500
 227 %!        0.0000   1.0000   0.0000  -1.0000  -2.0000   3.0000 ];
 228 %!
 229 %! Be = [ 0.5000   0.7500
 230 %!        0.5000  -0.2500
 231 %!        0.0000   0.5000
 232 %!        0.0000   0.5000
 233 %!       -0.5000   0.2500
 234 %!        0.0000   1.0000 ];
 235 %!
 236 %! Ce = [ 1.5000  -1.2500   0.5000  -0.2500  -0.2500  -0.2500
 237 %!        0.0000   0.5000   0.0000  -0.5000  -0.5000   0.5000 ];
 238 %!
 239 %! De = [ 0.5000  -0.2500
 240 %!        0.0000   0.5000 ];
 241 %!
 242 %! Me = [Ae, Be; Ce, De];
 243 %!
 244 %!assert (M, Me, 1e-4);
 245
 246
 247 ## sensitivity function
 248 ## Note the correct physical meaning of the states.
 249 ## Test would fail on a commercial octave clone
 250 ## because of wrong signs of matrices B and C.
 251 ## NOTE: Don't use T = I - S for complementary sensitivity,
 252 ##       use T = feedback (L) instead!
 253 %!shared S1, S2
 254 %! P = ss (-2, 3, 4, 5);  # meaningless numbers
 255 %! C = ss (-1, 1, 1, 0);  # ditto
 256 %! L = P * C;
 257 %! I = eye (size (L));
 258 %! S1 = feedback (I, L*-I, "+");  # draw a block diagram for better understanding
 259 %! S2 = inv (I + L);
 260 %!assert (S1.a, S2.a, 1e-4);
 261 %!assert (S1.b, S2.b, 1e-4);
 262 %!assert (S1.c, S2.c, 1e-4);
 263 %!assert (S1.d, S2.d, 1e-4);