1 ## Copyright (C) 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{K}, @var{N}, @var{gamma}, @var{info}] =} ncfsyn (@var{G}, @var{W1}, @var{W2}, @var{factor})
20 ## Loop shaping H-infinity synthesis. Compute positive feedback controller using
21 ## the McFarlane/Glover normalized coprime factor (NCF) loop shaping design procedure.
26 ## LTI model of plant.
28 ## LTI model of precompensator. Model must be SISO or of appropriate size.
29 ## An identity matrix is taken if @var{W1} is not specified or if an empty model
30 ## @code{[]} is passed.
32 ## LTI model of postcompensator. Model must be SISO or of appropriate size.
33 ## An identity matrix is taken if @var{W2} is not specified or if an empty model
34 ## @code{[]} is passed.
36 ## @code{factor = 1} implies that an optimal controller is required.
37 ## @code{factor > 1} implies that a suboptimal controller is required,
38 ## achieving a performance that is @var{factor} times less than optimal.
39 ## Default value is 1.
45 ## State-space model of the H-infinity loop-shaping controller.
47 ## State-space model of the closed loop depicted below.
49 ## L-infinity norm of @var{N}. @code{gamma = norm (N, inf)}.
51 ## Structure containing additional information.
53 ## Nugap robustness. @code{emax = inv (gamma)}.
55 ## Shaped plant. @code{Gs = W2 * G * W1}.
57 ## Controller for shaped plant. @code{Ks = ncfsyn (Gs)}.
59 ## Estimates of the reciprocal condition numbers of the Riccati equations
60 ## and a few other things. For details, see the description of the
61 ## corresponding SLICOT algorithm.
64 ## @strong{Block Diagram of N}
70 ## w1 + | +--------+ | +--------+
71 ## ----->(+)---+-->| Ks |----+--->(+)---->| Gs |----+
72 ## ^ + +--------+ ^ +--------+ |
75 ## +-------------------------------------------------+
79 ## @strong{Algorithm}@*
80 ## Uses SLICOT SB10ID, SB10KD and SB10ZD by courtesy of
81 ## @uref{http://www.slicot.org, NICONET e.V.}
84 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
88 function [K, varargout] = ncfsyn (G, W1 = [], W2 = [], factor = 1.0)
90 if (nargin == 0 || nargin > 4)
95 error ("ncfsyn: first argument must be an LTI system");
98 if (! is_real_scalar (factor) || factor < 1.0)
99 error ("ncfsyn: fourth argument invalid");
104 W1 = __adjust_weighting__ (W1, m);
105 W2 = __adjust_weighting__ (W2, p);
107 Gs = W2 * G * W1; # shaped plant
109 [a, b, c, d, tsam] = ssdata (Gs);
112 if (isct (Gs)) # continuous-time
113 [ak, bk, ck, dk, rcond] = slsb10id (a, b, c, d, factor);
114 elseif (any (d(:))) # discrete-time, d != 0
115 [ak, bk, ck, dk, rcond] = slsb10zd (a, b, c, d, factor, 0.0);
116 else # discrete-time, d == 0
117 [ak, bk, ck, dk, rcond] = slsb10kd (a, b, c, factor);
121 Ks = ss (ak, bk, ck, dk, tsam);
126 ## FIXME: is this really the same thing as the dark side does?
127 N = append (eye (p), Ks, Gs);
128 M = [zeros(p,p), zeros(p,m), eye(p);
129 eye(p), zeros(p,m), zeros(p,p);
130 zeros(m,p), eye(m), zeros(m,p)];
131 in_idx = [1:p, 2*p+(1:m)];
133 N = mconnect (N, M, in_idx, out_idx);
136 gamma = norm (N, inf);
137 varargout{2} = gamma;
139 varargout{3} = struct ("emax", inv (gamma), "Gs", Gs, "Ks", Ks, "rcond", rcond);
147 function W = __adjust_weighting__ (W, s)
153 ## if (! isstable (W))
154 ## error ("ncfsyn: %s must be stable", inputname (1));
156 ## if (! isminimumphase (W))
157 ## error ("ncfsyn: %s must be minimum-phase", inputname (1));
160 if (m == s && p == s) # model is of correct size
162 elseif (m == 1 && p == 1) # model is SISO
165 W = append (W, tmp); # stack SISO model s times
167 else # model is invalid
168 error ("ncfsyn: %s must have 1 or %d inputs and outputs", inputname (1), s);
175 ## continuous-time case, direct access to sb10id
176 %!shared AK, BK, CK, DK, RCOND, AKe, BKe, CKe, DKe, RCONDe
177 %! A = [ -1.0 0.0 4.0 5.0 -3.0 -2.0
178 %! -2.0 4.0 -7.0 -2.0 0.0 3.0
179 %! -6.0 9.0 -5.0 0.0 2.0 -1.0
180 %! -8.0 4.0 7.0 -1.0 -3.0 0.0
181 %! 2.0 5.0 8.0 -9.0 1.0 -4.0
182 %! 3.0 -5.0 8.0 0.0 2.0 -6.0 ];
191 %! C = [ 1.0 -1.0 2.0 -4.0 0.0 -3.0
192 %! -3.0 0.0 5.0 -1.0 1.0 1.0
193 %! -7.0 5.0 0.0 -8.0 2.0 -2.0 ];
201 %! [AK, BK, CK, DK, RCOND] = slsb10id (A, B, C, D, FACTOR);
203 %! AKe = [ -39.0671 9.9293 22.2322 -27.4113 43.8655
204 %! -6.6117 3.0006 11.0878 -11.4130 15.4269
205 %! 33.6805 -6.6934 -23.9953 14.1438 -33.4358
206 %! -32.3191 9.7316 25.4033 -24.0473 42.0517
207 %! -44.1655 18.7767 34.8873 -42.4369 50.8437 ];
209 %! BKe = [ -10.2905 -16.5382 -10.9782
210 %! -4.3598 -8.7525 -5.1447
211 %! 6.5962 1.8975 6.2316
212 %! -9.8770 -14.7041 -11.8778
213 %! -9.6726 -22.7309 -18.2692 ];
215 %! CKe = [ -0.6647 -0.0599 -1.0376 0.5619 1.7297
216 %! -8.4202 3.9573 7.3094 -7.6283 10.6768 ];
218 %! DKe = [ 0.8466 0.4979 -0.6993
219 %! -1.2226 -4.8689 -4.5056 ];
221 %! RCONDe = [ 0.13861D-01 0.90541D-02 ].';
223 %!assert (AK, AKe, 1e-4);
224 %!assert (BK, BKe, 1e-4);
225 %!assert (CK, CKe, 1e-4);
226 %!assert (DK, DKe, 1e-4);
227 %!assert (RCOND, RCONDe, 1e-4);
230 ## continuous-time case
231 %!shared AK, BK, CK, DK, RCOND, AKe, BKe, CKe, DKe, RCONDe
232 %! A = [ -1.0 0.0 4.0 5.0 -3.0 -2.0
233 %! -2.0 4.0 -7.0 -2.0 0.0 3.0
234 %! -6.0 9.0 -5.0 0.0 2.0 -1.0
235 %! -8.0 4.0 7.0 -1.0 -3.0 0.0
236 %! 2.0 5.0 8.0 -9.0 1.0 -4.0
237 %! 3.0 -5.0 8.0 0.0 2.0 -6.0 ];
246 %! C = [ 1.0 -1.0 2.0 -4.0 0.0 -3.0
247 %! -3.0 0.0 5.0 -1.0 1.0 1.0
248 %! -7.0 5.0 0.0 -8.0 2.0 -2.0 ];
256 %! G = ss (A, B, C, D);
257 %! K = ncfsyn (G, [], [], FACTOR);
258 %! [AK, BK, CK, DK] = ssdata (K);
260 %! AKe = [ -39.0671 9.9293 22.2322 -27.4113 43.8655
261 %! -6.6117 3.0006 11.0878 -11.4130 15.4269
262 %! 33.6805 -6.6934 -23.9953 14.1438 -33.4358
263 %! -32.3191 9.7316 25.4033 -24.0473 42.0517
264 %! -44.1655 18.7767 34.8873 -42.4369 50.8437 ];
266 %! BKe = [ -10.2905 -16.5382 -10.9782
267 %! -4.3598 -8.7525 -5.1447
268 %! 6.5962 1.8975 6.2316
269 %! -9.8770 -14.7041 -11.8778
270 %! -9.6726 -22.7309 -18.2692 ];
272 %! CKe = [ -0.6647 -0.0599 -1.0376 0.5619 1.7297
273 %! -8.4202 3.9573 7.3094 -7.6283 10.6768 ];
275 %! DKe = [ 0.8466 0.4979 -0.6993
276 %! -1.2226 -4.8689 -4.5056 ];
278 %! RCONDe = [ 0.13861D-01 0.90541D-02 ];
280 %!assert (AK, AKe, 1e-4);
281 %!assert (BK, BKe, 1e-4);
282 %!assert (CK, CKe, 1e-4);
283 %!assert (DK, DKe, 1e-4);
286 ## discrete-time case D==0, direct access to sb10kd
287 %!shared AK, BK, CK, DK, RCOND, AKe, BKe, CKe, DKe, RCONDe
288 %! A = [ 0.2 0.0 0.3 0.0 -0.3 -0.1
289 %! -0.3 0.2 -0.4 -0.3 0.0 0.0
290 %! -0.1 0.1 -0.1 0.0 0.0 -0.3
291 %! 0.1 0.0 0.0 -0.1 -0.1 0.0
292 %! 0.0 0.3 0.6 0.2 0.1 -0.4
293 %! 0.2 -0.4 0.0 0.0 0.2 -0.2 ];
302 %! C = [ 1.0 -1.0 2.0 -2.0 0.0 -3.0
303 %! -3.0 0.0 1.0 -1.0 1.0 -1.0 ];
307 %! [AK, BK, CK, DK, RCOND] = slsb10kd (A, B, C, FACTOR);
309 %! AKe = [ 0.0337 0.0222 0.0858 0.1264 -0.1872 0.1547
310 %! 0.4457 0.0668 -0.2255 -0.3204 -0.4548 -0.0691
311 %! -0.2419 -0.2506 -0.0982 -0.1321 -0.0130 -0.0838
312 %! -0.4402 0.3654 -0.0335 -0.2444 0.6366 -0.6469
313 %! -0.3623 0.3854 0.4162 0.4502 0.0065 0.1261
314 %! -0.0121 -0.4377 0.0604 0.2265 -0.3389 0.4542 ];
316 %! BKe = [ 0.0931 -0.0269
323 %! CKe = [ -0.3677 0.2188 0.0403 -0.0854 0.3564 -0.3535
324 %! 0.1624 -0.0708 0.0058 0.0606 -0.2163 0.1802 ];
326 %! DKe = [ -0.0857 -0.0246
329 %! RCONDe = [ 0.11269D-01 0.17596D-01 0.18225D+00 0.75968D-03 ].';
331 %!assert (AK, AKe, 1e-4);
332 %!assert (BK, BKe, 1e-4);
333 %!assert (CK, CKe, 1e-4);
334 %!assert (DK, DKe, 1e-4);
335 %!assert (RCOND, RCONDe, 1e-4);
338 ## discrete-time case D==0
339 %!shared AK, BK, CK, DK, RCOND, AKe, BKe, CKe, DKe, RCONDe
340 %! A = [ 0.2 0.0 0.3 0.0 -0.3 -0.1
341 %! -0.3 0.2 -0.4 -0.3 0.0 0.0
342 %! -0.1 0.1 -0.1 0.0 0.0 -0.3
343 %! 0.1 0.0 0.0 -0.1 -0.1 0.0
344 %! 0.0 0.3 0.6 0.2 0.1 -0.4
345 %! 0.2 -0.4 0.0 0.0 0.2 -0.2 ];
354 %! C = [ 1.0 -1.0 2.0 -2.0 0.0 -3.0
355 %! -3.0 0.0 1.0 -1.0 1.0 -1.0 ];
359 %! G = ss (A, B, C, [], 1); # value of sampling time doesn't matter
360 %! K = ncfsyn (G, [], [], FACTOR);
361 %! [AK, BK, CK, DK] = ssdata (K);
363 %! AKe = [ 0.0337 0.0222 0.0858 0.1264 -0.1872 0.1547
364 %! 0.4457 0.0668 -0.2255 -0.3204 -0.4548 -0.0691
365 %! -0.2419 -0.2506 -0.0982 -0.1321 -0.0130 -0.0838
366 %! -0.4402 0.3654 -0.0335 -0.2444 0.6366 -0.6469
367 %! -0.3623 0.3854 0.4162 0.4502 0.0065 0.1261
368 %! -0.0121 -0.4377 0.0604 0.2265 -0.3389 0.4542 ];
370 %! BKe = [ 0.0931 -0.0269
377 %! CKe = [ -0.3677 0.2188 0.0403 -0.0854 0.3564 -0.3535
378 %! 0.1624 -0.0708 0.0058 0.0606 -0.2163 0.1802 ];
380 %! DKe = [ -0.0857 -0.0246
383 %! RCONDe = [ 0.11269D-01 0.17596D-01 0.18225D+00 0.75968D-03 ].';
385 %!assert (AK, AKe, 1e-4);
386 %!assert (BK, BKe, 1e-4);
387 %!assert (CK, CKe, 1e-4);
388 %!assert (DK, DKe, 1e-4);
391 ## discrete-time case D!=0, direct access to sb10zd
392 %!shared AK, BK, CK, DK, RCOND, AKe, BKe, CKe, DKe, RCONDe
393 %! A = [ 0.2 0.0 3.0 0.0 -0.3 -0.1
394 %! -3.0 0.2 -0.4 -0.3 0.0 0.0
395 %! -0.1 0.1 -1.0 0.0 0.0 -3.0
396 %! 1.0 0.0 0.0 -1.0 -1.0 0.0
397 %! 0.0 0.3 0.6 2.0 0.1 -0.4
398 %! 0.2 -4.0 0.0 0.0 0.2 -2.0 ];
407 %! C = [ 1.0 -1.0 2.0 -2.0 0.0 -3.0
408 %! -3.0 0.0 1.0 -1.0 1.0 -1.0
409 %! 2.0 4.0 -3.0 0.0 5.0 1.0 ];
417 %! [AK, BK, CK, DK, RCOND] = slsb10zd (A, B, C, D, FACTOR, 0.0);
419 %! AKe = [ 1.0128 0.5101 -0.1546 1.1300 3.3759 0.4911
420 %! -2.1257 -1.4517 -0.4486 0.3493 -1.5506 -1.4296
421 %! -1.0930 -0.6026 -0.1344 0.2253 -1.5625 -0.6762
422 %! 0.3207 0.1698 0.2376 -1.1781 -0.8705 0.2896
423 %! 0.5017 0.9006 0.0668 2.3613 0.2049 0.3703
424 %! 1.0787 0.6703 0.2783 -0.7213 0.4918 0.7435 ];
426 %! BKe = [ 0.4132 0.3112 -0.8077
427 %! 0.2140 0.4253 0.1811
428 %! -0.0710 0.0807 0.3558
429 %! -0.0121 -0.2019 0.0249
430 %! 0.1047 0.1399 -0.0457
431 %! -0.2542 -0.3472 0.0523 ];
433 %! CKe = [ -0.0372 -0.0456 -0.0040 0.0962 -0.2059 -0.0571
434 %! 0.1999 0.2994 0.1335 -0.0251 -0.3108 0.2048 ];
436 %! DKe = [ 0.0629 -0.0022 0.0363
437 %! -0.0228 0.0195 0.0600 ];
439 %! RCONDe = [ 0.27949D-03 0.66679D-03 0.45677D-01 0.23433D-07 0.68495D-01 0.76854D-01 ].';
441 %!assert (AK, AKe, 1e-4);
442 %!assert (BK, BKe, 1e-4);
443 %!assert (CK, CKe, 1e-4);
444 %!assert (DK, DKe, 1e-4);
445 %!assert (RCOND, RCONDe, 1e-4);
448 ## discrete-time case D!=0
449 %!shared AK, BK, CK, DK, RCOND, AKe, BKe, CKe, DKe, RCONDe
450 %! A = [ 0.2 0.0 3.0 0.0 -0.3 -0.1
451 %! -3.0 0.2 -0.4 -0.3 0.0 0.0
452 %! -0.1 0.1 -1.0 0.0 0.0 -3.0
453 %! 1.0 0.0 0.0 -1.0 -1.0 0.0
454 %! 0.0 0.3 0.6 2.0 0.1 -0.4
455 %! 0.2 -4.0 0.0 0.0 0.2 -2.0 ];
464 %! C = [ 1.0 -1.0 2.0 -2.0 0.0 -3.0
465 %! -3.0 0.0 1.0 -1.0 1.0 -1.0
466 %! 2.0 4.0 -3.0 0.0 5.0 1.0 ];
474 %! G = ss (A, B, C, D, 1); # value of sampling time doesn't matter
475 %! K = ncfsyn (G, [], [], FACTOR);
476 %! [AK, BK, CK, DK] = ssdata (K);
478 %! AKe = [ 1.0128 0.5101 -0.1546 1.1300 3.3759 0.4911
479 %! -2.1257 -1.4517 -0.4486 0.3493 -1.5506 -1.4296
480 %! -1.0930 -0.6026 -0.1344 0.2253 -1.5625 -0.6762
481 %! 0.3207 0.1698 0.2376 -1.1781 -0.8705 0.2896
482 %! 0.5017 0.9006 0.0668 2.3613 0.2049 0.3703
483 %! 1.0787 0.6703 0.2783 -0.7213 0.4918 0.7435 ];
485 %! BKe = [ 0.4132 0.3112 -0.8077
486 %! 0.2140 0.4253 0.1811
487 %! -0.0710 0.0807 0.3558
488 %! -0.0121 -0.2019 0.0249
489 %! 0.1047 0.1399 -0.0457
490 %! -0.2542 -0.3472 0.0523 ];
492 %! CKe = [ -0.0372 -0.0456 -0.0040 0.0962 -0.2059 -0.0571
493 %! 0.1999 0.2994 0.1335 -0.0251 -0.3108 0.2048 ];
495 %! DKe = [ 0.0629 -0.0022 0.0363
496 %! -0.0228 0.0195 0.0600 ];
498 %! RCONDe = [ 0.27949D-03 0.66679D-03 0.45677D-01 0.23433D-07 0.68495D-01 0.76854D-01 ].';
500 %!assert (AK, AKe, 1e-4);
501 %!assert (BK, BKe, 1e-4);
502 %!assert (CK, CKe, 1e-4);
503 %!assert (DK, DKe, 1e-4);