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