|
7017
|
1 ## Copyright (C) 1996, 1998, 2000, 2002, 2004, 2005, 2006, 2007 |
|
|
2 ## Auburn University. All rights reserved. |
|
3440
|
3 ## |
|
|
4 ## This file is part of Octave. |
|
|
5 ## |
|
|
6 ## Octave is free software; you can redistribute it and/or modify it |
|
7016
|
7 ## under the terms of the GNU General Public License as published by |
|
|
8 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
9 ## your option) any later version. |
|
3440
|
10 ## |
|
7016
|
11 ## Octave is distributed in the hope that it will be useful, but |
|
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
14 ## General Public License for more details. |
|
3440
|
15 ## |
|
|
16 ## You should have received a copy of the GNU General Public License |
|
7016
|
17 ## along with Octave; see the file COPYING. If not, see |
|
|
18 ## <http://www.gnu.org/licenses/>. |
|
3440
|
19 |
|
|
20 ## -*- texinfo -*- |
|
3502
|
21 ## @deftypefn {Function File} {[@var{retval}, @var{dgkf_struct} ] =} is_dgkf (@var{asys}, @var{nu}, @var{ny}, @var{tol} ) |
|
3440
|
22 ## Determine whether a continuous time state space system meets |
|
5016
|
23 ## assumptions of @acronym{DGKF} algorithm. |
|
3440
|
24 ## Partitions system into: |
|
|
25 ## @example |
|
5016
|
26 ## [dx/dt] [A | Bw Bu ][w] |
|
|
27 ## [ z ] = [Cz | Dzw Dzu ][u] |
|
3440
|
28 ## [ y ] [Cy | Dyw Dyu ] |
|
|
29 ## @end example |
|
|
30 ## or similar discrete-time system. |
|
3502
|
31 ## If necessary, orthogonal transformations @var{qw}, @var{qz} and nonsingular |
|
|
32 ## transformations @var{ru}, @var{ry} are applied to respective vectors |
|
5016
|
33 ## @var{w}, @var{z}, @var{u}, @var{y} in order to satisfy @acronym{DGKF} assumptions. |
|
3502
|
34 ## Loop shifting is used if @var{dyu} block is nonzero. |
|
3440
|
35 ## |
|
|
36 ## @strong{Inputs} |
|
|
37 ## @table @var |
|
3502
|
38 ## @item asys |
|
3440
|
39 ## system data structure |
|
|
40 ## @item nu |
|
|
41 ## number of controlled inputs |
|
|
42 ## @item ny |
|
|
43 ## number of measured outputs |
|
|
44 ## @item tol |
|
5016
|
45 ## threshold for 0; default: 200*@code{eps}. |
|
3440
|
46 ## @end table |
|
|
47 ## @strong{Outputs} |
|
|
48 ## @table @var |
|
|
49 ## @item retval |
|
|
50 ## true(1) if system passes check, false(0) otherwise |
|
|
51 ## @item dgkf_struct |
|
5016
|
52 ## data structure of @command{is_dgkf} results. Entries: |
|
3440
|
53 ## @table @var |
|
|
54 ## @item nw |
|
|
55 ## @itemx nz |
|
|
56 ## dimensions of @var{w}, @var{z} |
|
3502
|
57 ## @item a |
|
|
58 ## system @math{A} matrix |
|
|
59 ## @item bw |
|
|
60 ## (@var{n} x @var{nw}) @var{qw}-transformed disturbance input matrix |
|
|
61 ## @item bu |
|
|
62 ## (@var{n} x @var{nu}) @var{ru}-transformed controlled input matrix; |
|
3440
|
63 ## |
|
4946
|
64 ## @math{B = [Bw Bu]} |
|
3502
|
65 ## @item cz |
|
3440
|
66 ## (@var{nz} x @var{n}) Qz-transformed error output matrix |
|
3502
|
67 ## @item cy |
|
|
68 ## (@var{ny} x @var{n}) @var{ry}-transformed measured output matrix |
|
3440
|
69 ## |
|
4946
|
70 ## @math{C = [Cz; Cy]} |
|
3502
|
71 ## @item dzu |
|
|
72 ## @item dyw |
|
|
73 ## off-diagonal blocks of transformed system @math{D} matrix that enter |
|
3440
|
74 ## @var{z}, @var{y} from @var{u}, @var{w} respectively |
|
3502
|
75 ## @item ru |
|
3440
|
76 ## controlled input transformation matrix |
|
3502
|
77 ## @item ry |
|
3440
|
78 ## observed output transformation matrix |
|
3502
|
79 ## @item dyu_nz |
|
|
80 ## nonzero if the @var{dyu} block is nonzero. |
|
|
81 ## @item dyu |
|
|
82 ## untransformed @var{dyu} block |
|
3440
|
83 ## @item dflg |
|
|
84 ## nonzero if the system is discrete-time |
|
|
85 ## @end table |
|
|
86 ## @end table |
|
|
87 ## @code{is_dgkf} exits with an error if the system is mixed |
|
5016
|
88 ## discrete/continuous. |
|
3440
|
89 ## |
|
|
90 ## @strong{References} |
|
|
91 ## @table @strong |
|
|
92 ## @item [1] |
|
5016
|
93 ## Doyle, Glover, Khargonekar, Francis, @cite{State Space Solutions to Standard} |
|
|
94 ## @iftex |
|
|
95 ## @tex |
|
|
96 ## $ { \cal H }_2 $ @cite{and} $ { \cal H }_\infty $ |
|
|
97 ## @end tex |
|
|
98 ## @end iftex |
|
|
99 ## @ifinfo |
|
|
100 ## @cite{H-2 and H-infinity} |
|
|
101 ## @end ifinfo |
|
|
102 ## @cite{Control Problems}, @acronym{IEEE} @acronym{TAC} August 1989. |
|
3440
|
103 ## @item [2] |
|
5016
|
104 ## Maciejowksi, J.M., @cite{Multivariable Feedback Design}, Addison-Wesley, 1989. |
|
3440
|
105 ## @end table |
|
|
106 ## @end deftypefn |
|
|
107 |
|
|
108 ## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> |
|
|
109 ## Updated by John Ingram July 1996 to accept structured systems |
|
|
110 |
|
|
111 ## Revised by Kai P. Mueller April 1998 to solve the general H_infinity |
|
|
112 ## problem using unitary transformations Q (on w and z) |
|
|
113 ## and non-singular transformations R (on u and y) such |
|
|
114 ## that the Dzu and Dyw matrices of the transformed plant |
|
|
115 ## |
|
|
116 ## ~ |
|
|
117 ## P (the variable Asys here) |
|
|
118 ## |
|
|
119 ## become |
|
|
120 ## |
|
|
121 ## ~ -1 T |
|
|
122 ## D = Q D R = [ 0 I ] or [ I ], |
|
|
123 ## 12 12 12 12 |
|
|
124 ## |
|
|
125 ## ~ T |
|
|
126 ## D = R D Q = [ 0 I ] or [ I ]. |
|
|
127 ## 21 21 21 21 |
|
|
128 ## |
|
|
129 ## This transformation together with the algorithm in [1] solves |
|
|
130 ## the general problem (see [2] for example). |
|
|
131 |
|
|
132 function [retval, dgkf_struct] = is_dgkf (Asys, nu, ny, tol) |
|
|
133 |
|
|
134 if (nargin < 3) | (nargin > 4) |
|
6046
|
135 print_usage (); |
|
4030
|
136 elseif (! isscalar(nu) | ! isscalar(ny) ) |
|
3440
|
137 error("is_dgkf: arguments 2 and 3 must be scalars") |
|
4030
|
138 elseif (! isstruct(Asys) ) |
|
3440
|
139 error("Argument 1 must be a system data structure"); |
|
|
140 endif |
|
|
141 if(nargin < 4) |
|
|
142 tol = 200*eps; |
|
|
143 elseif( !is_sample(tol) ) |
|
|
144 error("is_dgkf: tol must be a positive scalar") |
|
|
145 endif |
|
|
146 |
|
|
147 retval = 1; # assume passes test |
|
|
148 |
|
|
149 dflg = is_digital(Asys); |
|
|
150 [Anc, Anz, nin, nout ] = sysdimensions(Asys); |
|
|
151 |
|
|
152 if( Anz == 0 & Anc == 0 ) |
|
|
153 error("is_dgkf: no system states"); |
|
|
154 elseif( nu >= nin ) |
|
|
155 error("is_dgkf: insufficient number of disturbance inputs"); |
|
|
156 elseif( ny >= nout ) |
|
|
157 error("is_dgkf: insufficient number of regulated outputs"); |
|
|
158 endif |
|
|
159 |
|
|
160 nw = nin - nu; nw1 = nw + 1; |
|
|
161 nz = nout - ny; nz1 = nz + 1; |
|
|
162 |
|
|
163 [A,B,C,D] = sys2ss(Asys); |
|
|
164 ## scale input/output for numerical reasons |
|
|
165 if(norm (C, "fro") * norm (B, "fro") == 0) |
|
|
166 error("||C||*||B|| = 0; no dynamic connnection from inputs to outputs"); |
|
|
167 endif |
|
|
168 xx = sqrt(norm(B, Inf) / norm(C, Inf)); |
|
|
169 B = B / xx; C = C * xx; |
|
|
170 |
|
|
171 ## partition matrices |
|
|
172 Bw = B(:,1:nw); Bu = B(:,nw1:nin); |
|
|
173 Cz = C(1:nz,:); Dzw = D(1:nz,1:nw); Dzu = D(1:nz,nw1:nin); |
|
|
174 Cy = C(nz1:nout,:); Dyw = D(nz1:nout,1:nw); Dyu = D(nz1:nout,nw1:nin); |
|
|
175 |
|
|
176 ## Check for loopo shifting |
|
|
177 Dyu_nz = (norm(Dyu,Inf) != 0); |
|
|
178 if (Dyu_nz) |
|
|
179 warning("is_dgkf: D22 nonzero; performing loop shifting"); |
|
|
180 endif |
|
|
181 |
|
|
182 ## 12 - rank condition at w = 0 |
|
|
183 xx =[A, Bu; Cz, Dzu]; |
|
|
184 [nr, nc] = size(xx); |
|
|
185 irank = rank(xx); |
|
|
186 if (irank != nc) |
|
|
187 retval = 0; |
|
|
188 warning(sprintf("rank([A Bu; Cz Dzu]) = %d, need %d; n=%d, nz=%d, nu=%d", ... |
|
|
189 irank,nc,(Anc+Anz),nz,nu)); |
|
|
190 warning(" *** 12-rank condition violated at w = 0."); |
|
|
191 endif |
|
|
192 |
|
|
193 ## 21 - rank condition at w = 0 |
|
|
194 xx =[A, Bw; Cy, Dyw]; |
|
|
195 [nr, nc] = size(xx); |
|
|
196 irank = rank(xx); |
|
|
197 if (irank != nr) |
|
|
198 retval = 0; |
|
|
199 warning(sprintf("rank([A Bw; Cy Dyw]) = %d, need %d; n=%d, ny=%d, nw=%d", ... |
|
|
200 irank,nr,(Anc+Anz),ny,nw)); |
|
|
201 warning(" *** 21-rank condition violated at w = 0."); |
|
|
202 endif |
|
|
203 |
|
|
204 ## can Dzu be transformed to become [0 I]' or [I]? |
|
|
205 ## This ensures a normalized weight |
|
|
206 [Qz, Ru] = qr(Dzu); |
|
|
207 irank = rank(Ru); |
|
|
208 if (irank != nu) |
|
|
209 retval = 0; |
|
|
210 warning(sprintf("*** rank(Dzu(%d x %d) = %d", nz, nu, irank)); |
|
|
211 warning(" *** Dzu does not have full column rank."); |
|
|
212 endif |
|
|
213 if (nu >= nz) |
|
|
214 Qz = Qz(:,1:nu)'; |
|
|
215 else |
|
|
216 Qz = [Qz(:,(nu+1):nz), Qz(:,1:nu)]'; |
|
|
217 endif |
|
|
218 Ru = Ru(1:nu,:); |
|
|
219 |
|
|
220 ## can Dyw be transformed to become [0 I] or [I]? |
|
|
221 ## This ensures a normalized weight |
|
|
222 [Qw, Ry] = qr(Dyw'); |
|
|
223 irank = rank(Ry); |
|
|
224 if (irank != ny) |
|
|
225 retval = 0; |
|
|
226 warning(sprintf("*** rank(Dyw(%d x %d) = %d", ny, nw, irank)); |
|
|
227 warning(" *** Dyw does not have full row rank."); |
|
|
228 endif |
|
|
229 |
|
|
230 if (ny >= nw) |
|
|
231 Qw = Qw(:,1:ny); |
|
|
232 else |
|
|
233 Qw = [Qw(:,(ny+1):nw), Qw(:,1:ny)]; |
|
|
234 endif |
|
|
235 Ry = Ry(1:ny,:)'; |
|
|
236 |
|
|
237 ## transform P by Qz/Ru and Qw/Ry |
|
|
238 Bw = Bw*Qw; |
|
|
239 Bu = Bu/Ru; |
|
|
240 B = [Bw, Bu]; |
|
|
241 Cz = Qz*Cz; |
|
|
242 Cy = Ry\Cy; |
|
|
243 C = [Cz; Cy]; |
|
|
244 Dzw = Qz*Dzw*Qw; |
|
|
245 Dzu = Qz*Dzu/Ru; |
|
|
246 Dyw = Ry\Dyw*Qw; |
|
|
247 |
|
|
248 ## pack the return structure |
|
|
249 dgkf_struct.nw = nw; |
|
|
250 dgkf_struct.nu = nu; |
|
|
251 dgkf_struct.nz = nz; |
|
|
252 dgkf_struct.ny = ny; |
|
|
253 dgkf_struct.A = A; |
|
|
254 dgkf_struct.Bw = Bw; |
|
|
255 dgkf_struct.Bu = Bu; |
|
|
256 dgkf_struct.Cz = Cz; |
|
|
257 dgkf_struct.Cy = Cy; |
|
|
258 dgkf_struct.Dzw = Dzw; |
|
|
259 dgkf_struct.Dzu = Dzu; |
|
|
260 dgkf_struct.Dyw = Dyw; |
|
|
261 dgkf_struct.Dyu = Dyu; |
|
|
262 dgkf_struct.Ru = Ru; |
|
|
263 dgkf_struct.Ry = Ry; |
|
|
264 dgkf_struct.Dyu_nz = Dyu_nz; |
|
|
265 dgkf_struct.dflg = dflg; |
|
|
266 |
|
|
267 endfunction |
