1 ## Copyright (C) 2009, 2010, 2011 Lukas F. Reichlin
3 ## This file is part of LTI Syncope.
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.
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.
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/>.
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.
30 ## LTI model of forward transmission. @code{[p1, m1] = size (sys1)}.
32 ## LTI model of backward transmission.
33 ## If not specified, an identity matrix of appropriate size is taken.
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.
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.
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.
51 ## Resulting LTI model.
54 ## @strong{Block Diagram}
58 ## ------>(+)----->| sys1 |-------+------->
62 ## +-------| sys2 |<------+
68 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
69 ## Created: October 2009
72 function sys = feedback (sys1, sys2, feedin, feedout, fbsign = -1)
74 [p1, m1] = size (sys1);
77 case 1 # sys = feedback (sys)
79 error ("feedback: argument must be a square system");
86 if (ischar (sys2)) # sys = feedback (sys, "+")
88 error ("feedback: first argument must be a square system");
90 fbsign = __check_fbsign__ (sys2);
92 endif # sys = feedback (sys1, sys2)
96 case 3 # sys = feedback (sys1, sys2, "+")
97 fbsign = __check_fbsign__ (feedin);
101 case 4 # sys = feedback (sys1, sys2, feedin, feedout)
102 ## nothing needs to be done here
103 ## case 4 required to prevent "otherwise"
105 case 5 # sys = feedback (sys1, sys2, feedin, feedout, "+")
106 fbsign = __check_fbsign__ (fbsign);
113 [p2, m2] = size (sys2);
115 l_feedin = length (feedin);
116 l_feedout = length (feedout);
119 error ("feedback: feedin indices: %d, outputs sys2: %d", l_feedin, p2);
123 error ("feedback: feedout indices: %d, inputs sys2: %d", l_feedout, m2);
126 if (any (feedin > m1 | feedin < 1))
127 error ("feedback: range of feedin indices exceeds dimensions of sys1");
130 if (any (feedin > p1 | feedin < 1))
131 error ("feedback: range of feedout indices exceeds dimensions of sys1");
134 M11 = zeros (m1, p1);
135 M22 = zeros (m2, p2);
137 M12 = full (sparse (feedin, 1:l_feedin, fbsign, m1, p2));
138 M21 = full (sparse (1:l_feedout, feedout, 1, m2, p1));
140 ## NOTE: for-loops do NOT the same as
141 ## M12(feedin, 1:l_feedin) = fbsign;
142 ## M21(1:l_feedout, feedout) = 1;
144 ## M12 = zeros (m1, p2);
145 ## M21 = zeros (m2, p1);
147 ## for k = 1 : l_feedin
148 ## M12(feedin(k), k) = fbsign;
151 ## for k = 1 : l_feedout
152 ## M21(k, feedout(k)) = 1;
161 sys = __sys_group__ (sys1, sys2);
162 sys = __sys_connect__ (sys, M);
163 sys = __sys_prune__ (sys, out_idx, in_idx);
168 function fbsign = __check_fbsign__ (fbsign)
170 if (is_real_scalar (fbsign))
171 fbsign = sign (fbsign);
172 elseif (ischar (fbsign))
173 if (strcmp (fbsign, "+"))
175 elseif (strcmp (fbsign, "-"))
178 error ("feedback: invalid feedback sign string");
181 error ("feedback: invalid feedback sign type");
187 ## Feedback inter-connection of two systems in state-space form
188 ## Test from SLICOT AB05ND
190 %! A1 = [ 1.0 0.0 -1.0
194 %! B1 = [ 1.0 1.0 0.0
197 %! C1 = [ 3.0 -2.0 1.0
203 %! A2 = [-3.0 0.0 0.0
207 %! B2 = [ 0.0 -1.0 0.0
210 %! C2 = [ 1.0 1.0 0.0
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);
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 ];
229 %! Be = [ 0.5000 0.7500
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 ];
239 %! De = [ 0.5000 -0.2500
242 %! Me = [Ae, Be; Ce, De];
244 %!assert (M, Me, 1e-4);
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!
254 %! P = ss (-2, 3, 4, 5); # meaningless numbers
255 %! C = ss (-1, 1, 1, 0); # ditto
257 %! I = eye (size (L));
258 %! S1 = feedback (I, L*-I, "+"); # draw a block diagram for better understanding
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);