Mercurial > octave-nkf

comparison scripts/control/hinf/is_dgkf.m @ 7131:a184bc985c37

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-11-08 15:55:02 by jwe]
author jwe
date Thu, 08 Nov 2007 15:55:02 +0000
parents a1dbe9d80eee
children df9519e9990c
comparison
equal deleted inserted replaced
7130:5eeb46c784d7 7131:a184bc985c37
129 ## This transformation together with the algorithm in [1] solves 129 ## This transformation together with the algorithm in [1] solves
130 ## the general problem (see [2] for example). 130 ## the general problem (see [2] for example).
131 131
132 function [retval, dgkf_struct] = is_dgkf (Asys, nu, ny, tol) 132 function [retval, dgkf_struct] = is_dgkf (Asys, nu, ny, tol)
133 133
134 if (nargin < 3) | (nargin > 4) 134 if (nargin < 3) || (nargin > 4)
135 print_usage (); 135 print_usage ();
136 elseif (! isscalar(nu) | ! isscalar(ny) ) 136 elseif (! isscalar (nu) || ! isscalar (ny))
137 error("is_dgkf: arguments 2 and 3 must be scalars") 137 error ("is_dgkf: arguments 2 and 3 must be scalars")
138 elseif (! isstruct(Asys) ) 138 elseif (! isstruct (Asys))
139 error("Argument 1 must be a system data structure"); 139 error ("Argument 1 must be a system data structure");
140 endif 140 endif
141 if(nargin < 4) 141 if (nargin < 4)
142 tol = 200*eps; 142 tol = 200*eps;
143 elseif( !is_sample(tol) ) 143 elseif (! is_sample (tol))
144 error("is_dgkf: tol must be a positive scalar") 144 error ("is_dgkf: tol must be a positive scalar")
145 endif 145 endif
146 146
147 retval = 1; # assume passes test 147 retval = 1; # assume passes test
148 148
149 dflg = is_digital(Asys); 149 dflg = is_digital (Asys);
150 [Anc, Anz, nin, nout ] = sysdimensions(Asys); 150 [Anc, Anz, nin, nout ] = sysdimensions (Asys);
151 151
152 if( Anz == 0 & Anc == 0 ) 152 if (Anz == 0 && Anc == 0)
153 error("is_dgkf: no system states"); 153 error ("is_dgkf: no system states");
154 elseif( nu >= nin ) 154 elseif (nu >= nin)
155 error("is_dgkf: insufficient number of disturbance inputs"); 155 error ("is_dgkf: insufficient number of disturbance inputs");
156 elseif( ny >= nout ) 156 elseif (ny >= nout)
157 error("is_dgkf: insufficient number of regulated outputs"); 157 error ("is_dgkf: insufficient number of regulated outputs");
158 endif 158 endif
159 159
160 nw = nin - nu; nw1 = nw + 1; 160 nw = nin - nu;
161 nz = nout - ny; nz1 = nz + 1; 161 nz = nout - ny;
162 162
163 [A,B,C,D] = sys2ss(Asys); 163 nw1 = nw + 1;
164 nz1 = nz + 1;
165
166 [A, B, C, D] = sys2ss (Asys);
164 ## scale input/output for numerical reasons 167 ## scale input/output for numerical reasons
165 if(norm (C, "fro") * norm (B, "fro") == 0) 168 if (norm (C, "fro") * norm (B, "fro") == 0)
166 error("||C||*||B|| = 0; no dynamic connnection from inputs to outputs"); 169 error ("||C||*||B|| = 0; no dynamic connnection from inputs to outputs");
167 endif 170 endif
168 xx = sqrt(norm(B, Inf) / norm(C, Inf)); 171 xx = sqrt (norm (B, Inf) / norm (C, Inf));
169 B = B / xx; C = C * xx; 172 B = B / xx; C = C * xx;
170 173
171 ## partition matrices 174 ## partition matrices
172 Bw = B(:,1:nw); Bu = B(:,nw1:nin); 175 Bw = B(:,1:nw);
173 Cz = C(1:nz,:); Dzw = D(1:nz,1:nw); Dzu = D(1:nz,nw1:nin); 176 Bu = B(:,nw1:nin);
174 Cy = C(nz1:nout,:); Dyw = D(nz1:nout,1:nw); Dyu = D(nz1:nout,nw1:nin); 177
178 Cz = C(1:nz,:);
179 Cy = C(nz1:nout,:);
180
181 Dzw = D(1:nz,1:nw);
182 Dzu = D(1:nz,nw1:nin);
183
184 Dyw = D(nz1:nout,1:nw);
185 Dyu = D(nz1:nout,nw1:nin);
175 186
176 ## Check for loopo shifting 187 ## Check for loopo shifting
177 Dyu_nz = (norm(Dyu,Inf) != 0); 188 Dyu_nz = (norm (Dyu, Inf) != 0);
178 if (Dyu_nz) 189 if (Dyu_nz)
179 warning("is_dgkf: D22 nonzero; performing loop shifting"); 190 warning ("is_dgkf: D22 nonzero; performing loop shifting");
180 endif 191 endif
181 192
182 ## 12 - rank condition at w = 0 193 ## 12 - rank condition at w = 0
183 xx =[A, Bu; Cz, Dzu]; 194 xx =[A, Bu; Cz, Dzu];
184 [nr, nc] = size(xx); 195 [nr, nc] = size (xx);
185 irank = rank(xx); 196 irank = rank (xx);
186 if (irank != nc) 197 if (irank != nc)
187 retval = 0; 198 retval = 0;
188 warning(sprintf("rank([A Bu; Cz Dzu]) = %d, need %d; n=%d, nz=%d, nu=%d", ... 199 warning ("rank([A Bu; Cz Dzu]) = %d, need %d; n=%d, nz=%d, nu=%d",
189 irank,nc,(Anc+Anz),nz,nu)); 200 irank, nc, Anc+Anz, nz, nu);
190 warning(" *** 12-rank condition violated at w = 0."); 201 warning (" *** 12-rank condition violated at w = 0");
191 endif 202 endif
192 203
193 ## 21 - rank condition at w = 0 204 ## 21 - rank condition at w = 0
194 xx =[A, Bw; Cy, Dyw]; 205 xx =[A, Bw; Cy, Dyw];
195 [nr, nc] = size(xx); 206 [nr, nc] = size (xx);
196 irank = rank(xx); 207 irank = rank (xx);
197 if (irank != nr) 208 if (irank != nr)
198 retval = 0; 209 retval = 0;
199 warning(sprintf("rank([A Bw; Cy Dyw]) = %d, need %d; n=%d, ny=%d, nw=%d", ... 210 warning ("rank([A Bw; Cy Dyw]) = %d, need %d; n=%d, ny=%d, nw=%d",
200 irank,nr,(Anc+Anz),ny,nw)); 211 irank, nr, Anc+Anz, ny, nw);
201 warning(" *** 21-rank condition violated at w = 0."); 212 warning (" *** 21-rank condition violated at w = 0");
202 endif 213 endif
203 214
204 ## can Dzu be transformed to become [0 I]' or [I]? 215 ## can Dzu be transformed to become [0 I]' or [I]?
205 ## This ensures a normalized weight 216 ## This ensures a normalized weight
206 [Qz, Ru] = qr(Dzu); 217 [Qz, Ru] = qr (Dzu);
207 irank = rank(Ru); 218 irank = rank (Ru);
208 if (irank != nu) 219 if (irank != nu)
209 retval = 0; 220 retval = 0;
210 warning(sprintf("*** rank(Dzu(%d x %d) = %d", nz, nu, irank)); 221 warning ("*** rank(Dzu(%d x %d) = %d", nz, nu, irank);
211 warning(" *** Dzu does not have full column rank."); 222 warning ("*** Dzu does not have full column rank");
212 endif 223 endif
213 if (nu >= nz) 224 if (nu >= nz)
214 Qz = Qz(:,1:nu)'; 225 Qz = Qz(:,1:nu)';
215 else 226 else
216 Qz = [Qz(:,(nu+1):nz), Qz(:,1:nu)]'; 227 Qz = [Qz(:,(nu+1):nz), Qz(:,1:nu)]';
217 endif 228 endif
218 Ru = Ru(1:nu,:); 229 Ru = Ru(1:nu,:);
219 230
220 ## can Dyw be transformed to become [0 I] or [I]? 231 ## can Dyw be transformed to become [0 I] or [I]?
221 ## This ensures a normalized weight 232 ## This ensures a normalized weight
222 [Qw, Ry] = qr(Dyw'); 233 [Qw, Ry] = qr (Dyw');
223 irank = rank(Ry); 234 irank = rank (Ry);
224 if (irank != ny) 235 if (irank != ny)
225 retval = 0; 236 retval = 0;
226 warning(sprintf("*** rank(Dyw(%d x %d) = %d", ny, nw, irank)); 237 warning ("*** rank(Dyw(%d x %d) = %d", ny, nw, irank);
227 warning(" *** Dyw does not have full row rank."); 238 warning (" *** Dyw does not have full row rank");
228 endif 239 endif
229 240
230 if (ny >= nw) 241 if (ny >= nw)
231 Qw = Qw(:,1:ny); 242 Qw = Qw(:,1:ny);
232 else 243 else
244 Dzw = Qz*Dzw*Qw; 255 Dzw = Qz*Dzw*Qw;
245 Dzu = Qz*Dzu/Ru; 256 Dzu = Qz*Dzu/Ru;
246 Dyw = Ry\Dyw*Qw; 257 Dyw = Ry\Dyw*Qw;
247 258
248 ## pack the return structure 259 ## pack the return structure
249 dgkf_struct.nw = nw; 260 dgkf_struct.nw = nw;
250 dgkf_struct.nu = nu; 261 dgkf_struct.nu = nu;
251 dgkf_struct.nz = nz; 262 dgkf_struct.nz = nz;
252 dgkf_struct.ny = ny; 263 dgkf_struct.ny = ny;
253 dgkf_struct.A = A; 264 dgkf_struct.A = A;
254 dgkf_struct.Bw = Bw; 265 dgkf_struct.Bw = Bw;
255 dgkf_struct.Bu = Bu; 266 dgkf_struct.Bu = Bu;
256 dgkf_struct.Cz = Cz; 267 dgkf_struct.Cz = Cz;
257 dgkf_struct.Cy = Cy; 268 dgkf_struct.Cy = Cy;
258 dgkf_struct.Dzw = Dzw; 269 dgkf_struct.Dzw = Dzw;
259 dgkf_struct.Dzu = Dzu; 270 dgkf_struct.Dzu = Dzu;
260 dgkf_struct.Dyw = Dyw; 271 dgkf_struct.Dyw = Dyw;
261 dgkf_struct.Dyu = Dyu; 272 dgkf_struct.Dyu = Dyu;
262 dgkf_struct.Ru = Ru; 273 dgkf_struct.Ru = Ru;
263 dgkf_struct.Ry = Ry; 274 dgkf_struct.Ry = Ry;
264 dgkf_struct.Dyu_nz = Dyu_nz; 275 dgkf_struct.Dyu_nz = Dyu_nz;
265 dgkf_struct.dflg = dflg; 276 dgkf_struct.dflg = dflg;
266 277
267 endfunction 278 endfunction