|
3430
|
1 ## Copyright (C) 1996 Auburn University. All rights reserved. |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
7016
|
6 ## under the terms of the GNU General Public License as published by |
|
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
8 ## your option) any later version. |
|
3430
|
9 ## |
|
7016
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
3430
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
|
17 ## <http://www.gnu.org/licenses/>. |
|
3430
|
18 |
|
|
19 ## -*- texinfo -*- |
|
5016
|
20 ## @deftypefn {Function File} {@var{K} =} hinf_ctr (@var{dgs}, @var{f}, @var{h}, @var{z}, @var{g}) |
|
|
21 ## Called by @code{hinfsyn} to compute the |
|
|
22 ## @iftex |
|
|
23 ## @tex |
|
|
24 ## $ { \cal H }_\infty $ |
|
|
25 ## @end tex |
|
|
26 ## @end iftex |
|
|
27 ## @ifinfo |
|
|
28 ## H-infinity |
|
|
29 ## @end ifinfo |
|
|
30 ## optimal controller. |
|
3430
|
31 ## |
|
|
32 ## @strong{Inputs} |
|
|
33 ## @table @var |
|
|
34 ## @item dgs |
|
|
35 ## data structure returned by @code{is_dgkf} |
|
3502
|
36 ## @item f |
|
|
37 ## @itemx h |
|
3430
|
38 ## feedback and filter gain (not partitioned) |
|
|
39 ## @item g |
|
|
40 ## final gamma value |
|
|
41 ## @end table |
|
|
42 ## @strong{Outputs} |
|
5016
|
43 ## @table @var |
|
|
44 ## @item K |
|
3502
|
45 ## controller (system data structure) |
|
5016
|
46 ## @end table |
|
3430
|
47 ## |
|
|
48 ## Do not attempt to use this at home; no argument checking performed. |
|
|
49 ## @end deftypefn |
|
|
50 |
|
|
51 ## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> |
|
|
52 ## Created: August 1995 |
|
|
53 ## Revised by Kai P. Mueller April 1998 to solve the general H_infinity |
|
|
54 ## problem using unitary transformations Q (on w and z) |
|
|
55 ## and non-singular transformations R (on u and y). |
|
|
56 |
|
|
57 function K = hinf_ctr (dgs, F, H, Z, g) |
|
|
58 |
|
|
59 nw = dgs.nw; |
|
|
60 nu = dgs.nu; |
|
|
61 nz = dgs.nz; |
|
|
62 ny = dgs.ny; |
|
|
63 d22nz = dgs.Dyu_nz; |
|
|
64 |
|
|
65 B1 = dgs.Bw; |
|
|
66 B2 = dgs.Bu; |
|
|
67 C1 = dgs.Cz; |
|
|
68 C2 = dgs.Cy; |
|
|
69 C = [C1; C2]; |
|
|
70 D11 = dgs.Dzw; |
|
|
71 D12 = dgs.Dzu; |
|
|
72 D21 = dgs.Dyw; |
|
|
73 D22 = dgs.Dyu; |
|
|
74 A = dgs.A; |
|
|
75 Ru = dgs.Ru; |
|
|
76 Ry = dgs.Ry; |
|
|
77 |
|
|
78 nout = nz + ny; |
|
|
79 nin = nw + nu; |
|
|
80 nstates = size(A, 1); |
|
|
81 |
|
|
82 F11 = F(1:(nw-ny),:); |
|
|
83 F12 = F((nw-ny+1):nw,:); |
|
|
84 F2 = F((nw+1):nin,:); |
|
|
85 H11 = H(:,1:(nz-nu)); |
|
|
86 H12 = H(:,(nz-nu+1):nz); |
|
|
87 H2 = H(:,(nz+1):nout); |
|
|
88 |
|
|
89 ## D11 partitions |
|
|
90 D1111 = D11(1:(nz-nu),1:(nw-ny)); |
|
|
91 D1112 = D11(1:(nz-nu),(nw-ny+1):nw); |
|
|
92 D1121 = D11((nz-nu+1):nz,1:(nw-ny)); |
|
|
93 D1122 = D11((nz-nu+1):nz,(nw-ny+1):nw); |
|
|
94 |
|
|
95 ## D11ik may be the empty matrix, don't calculate with empty matrices |
|
|
96 [nd1111,md1111] = size(D1111); |
|
|
97 md1112 = length(D1112); |
|
|
98 md1121 = length(D1121); |
|
|
99 |
|
|
100 if ((nd1111 == 0) || (md1112 == 0)) |
|
|
101 d11hat = -D1122; |
|
|
102 else |
|
|
103 xx = inv(g*g*eye(nz-nu) - D1111*D1111'); |
|
|
104 d11hat = -D1121*D1111'*xx*D1112 - D1122; |
|
|
105 endif |
|
|
106 if (md1112 == 0) |
|
|
107 d21hat = eye(ny); |
|
|
108 elseif (nd1111 == 0) |
|
|
109 d21hat = chol(eye(ny) - D1112'*D1112/g/g); |
|
|
110 else |
|
|
111 xx = inv(g*g*eye(nz-nu) - D1111*D1111'); |
|
|
112 xx = eye(ny) - D1112'*xx*D1112; |
|
|
113 d21hat = chol(xx); |
|
|
114 endif |
|
|
115 if (md1121 == 0) |
|
|
116 d12hat = eye(nu); |
|
|
117 elseif (md1111 == 0) |
|
|
118 d12hat = chol(eye(nu) - D1121*D1121'/g/g)'; |
|
|
119 else |
|
|
120 xx = inv(g*g*eye(nw-ny) - D1111'*D1111); |
|
|
121 xx = eye(nu)-D1121*xx*D1121'; |
|
|
122 d12hat = chol(xx)'; |
|
|
123 endif |
|
|
124 |
|
|
125 b2hat = (B2+H12)*d12hat; |
|
|
126 c2hat = -d21hat*(C2+F12)*Z; |
|
|
127 b1hat = -H2 + (b2hat/d12hat)*d11hat; |
|
|
128 c1hat = F2*Z + (d11hat/d21hat)*c2hat; |
|
|
129 ahat = A + H*C + (b2hat/d12hat)*c1hat; |
|
|
130 |
|
|
131 ## rescale controller by Ru and Ry |
|
|
132 b1hat = b1hat/Ry; |
|
|
133 c1hat = Ru\c1hat; |
|
|
134 bhat = [b1hat, b2hat]; |
|
|
135 chat = [c1hat; c2hat]; |
|
|
136 dhat = [Ru\d11hat/Ry, Ru\d12hat; d21hat/Ry, 0*d11hat']; |
|
|
137 |
|
|
138 ## non-zero D22 is a special case |
|
|
139 if (d22nz) |
|
|
140 if (rank(eye(nu) + d11hat*D22) < nu) |
|
|
141 error(" *** cannot compute controller for D22 non-zero."); |
|
|
142 endif |
|
|
143 |
|
|
144 d22new = [D22, zeros(ny,ny); zeros(nu,nu), 0*D22']; |
|
|
145 xx = inv(eye(nu+ny) + d22new*dhat); |
|
|
146 mhat = inv(eye(nu+ny) + dhat*d22new); |
|
|
147 ahat = ahat - bhat*((eye(nu+ny)-xx)/dhat)*chat; |
|
|
148 bhat = bhat*xx; |
|
|
149 chat = mhat*chat; |
|
|
150 dhat = dhat*xx; |
|
|
151 |
|
|
152 endif |
|
|
153 |
|
4771
|
154 K = ss(ahat,bhat(:,1:ny),chat(1:nu,:),dhat(1:nu,1:ny)); |
|
3430
|
155 |
|
|
156 endfunction |
