1 ## Copyright (C) 2009 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{rcond}] =} h2syn (@var{P}, @var{nmeas}, @var{ncon})
20 ## H-2 control synthesis for LTI plant.
25 ## Generalized plant. Must be a proper/realizable LTI model.
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.
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.
39 ## State-space model of the H-2 optimal controller.
41 ## State-space model of the lower LFT of @var{P} and @var{K}.
43 ## H-2 norm of @var{N}.
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.
53 ## @strong{Block Diagram}
57 ## gamma = min||N(K)|| N = lft (P, K)
61 ## w ----->| |-----> z
63 ## u +---->| |-----+ v
67 ## +-----| K(s) |<----+
71 ## w ----->| N(s) |-----> z
76 ## @strong{Algorithm}@*
77 ## Uses SLICOT SB10HD and SB10ED by courtesy of
78 ## @uref{http://www.slicot.org, NICONET e.V.}
80 ## @seealso{augw, lqr, dlqr, kalman}
83 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
84 ## Created: December 2009
87 function [K, varargout] = h2syn (P, nmeas, ncon)
89 ## check input arguments
95 error ("h2syn: first argument must be an LTI system");
98 if (! is_real_scalar (nmeas))
99 error ("h2syn: second argument invalid");
102 if (! is_real_scalar (ncon))
103 error ("h2syn: third argument invalid");
106 [a, b, c, d, tsam] = ssdata (P);
108 ## check assumptions A1 - A3
117 if (isct (P) && any (d11(:)))
118 warning ("h2syn: setting matrice D11 to zero");
122 if (! isstabilizable (P(:, m1+1:m)))
123 error ("h2syn: (A, B2) must be stabilizable");
126 if (! isdetectable (P(p1+1:p, :)))
127 error ("h2syn: (C2, A) must be detectable");
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);
138 K = ss (ak, bk, ck, dk, tsam);
144 varargout{2} = norm (N, 2);
146 varargout{3} = rcond;
154 ## continuous-time case
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];
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];
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];
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];
182 %! P = ss (A, B, C, D);
183 %! K = h2syn (P, 2, 2);
184 %! M = [K.A, K.B; K.C, K.D];
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];
193 %! KB = [ 3.7345 3.4758
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];
203 %! KD = [ 0.0000 0.0000
206 %! M_exp = [KA, KB; KC, KD];
208 %!assert (M, M_exp, 1e-4);
211 ## discrete-time case
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];
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];
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];
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];
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];
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];
250 %! KB = [-0.5985 -0.5464
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];
260 %! KD = [-0.2181 -0.2070
263 %! M_exp = [KA, KB; KC, KD];
265 %!assert (M, M_exp, 1e-4);