home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
OS/2 Shareware BBS: 10 Tools
/
10-Tools.zip
/
octa21eb.zip
/
octave
/
SCRIPTS.ZIP
/
scripts.fat
/
control
/
is_dgkf.m
< prev
next >
Wrap
Text File
|
1999-04-29
|
8KB
|
229 lines
# Copyright (C) 1996,1998 A. Scottedward Hodel
#
# This file is part of Octave.
#
# Octave is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2, or (at your option) any
# later version.
#
# Octave is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License
# along with Octave; see the file COPYING. If not, write to the Free
# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
function [retval,dgkf_struct] = is_dgkf(Asys,nu,ny,tol)
# function [retval,dgkf_struct] = is_dgkf(Asys,nu,ny{,tol})
# Determine whether a continuous time state space system meets
# assumptions of DGKF algorithm.
# Partitions system into: [z] = [A | Bw Bu ][w]
# [y] [Cz | Dzw Dzu ][u]
# [Cy | Dyw Dyu ]
#If necessary, orthogonal transformations Qw, Qz and nonsingular
# transformations Ru, Ry are applied to respective vectors w, z, u, y
# in order to satisfy DGKF assumptions. Loop shifting is used if Dyu block
# is nonzero.
#
# inputs:
# Asys: structured system
# nu: number of controlled inputs
# ny: number of measured outputs
# tol: threshhold for 0. Default: 200*eps
# outputs:
# retval: true(1) if system passes check, false(0) otherwise
# dgkf_struct: data structure of is_dgkf results. Entries:
# nw, nz: dimensions of w, z
# A: system A matrix
# Bw: (n x nw) Qw-transformed disturbance input matrix
# Bu: (n x nu) Ru-transformed controlled input matrix;
# Note: B = [Bw Bu]
# Cz: (nz x n) Qz-transformed error output matrix
# Cy: (ny x n) Ry-transformed measured output matrix
# Note: C = [Cz; Cy]
# Dzu, Dyw: off-diagonal blocks of transformed D matrix that enter
# z, y from u, w respectively
# Ru: controlled input transformation matrix
# Ry: observed output transformation matrix
# Dyu_nz: nonzero if the Dyu block is nonzero.
# Dyu: untransformed Dyu block
# dflg: nonzero if the system is discrete-time
#
# is_dgkf exits with an error if the system is mixed discrete/continuous
#
# [1] Doyle, Glover, Khargonekar, Francis, "State Space Solutions
# to Standard H2 and Hinf Control Problems," IEEE TAC August 1989
#
# [2] Maciejowksi, J.M.: "Multivariable feedback design,"
# Addison-Wesley, 1989, ISBN 0-201-18243-2
#
# [3] Keith Glover and John C. Doyle: "State-space formulae for all
# stabilizing controllers that satisfy and h-infinity-norm bound
# and relations to risk sensitivity,"
# Systems & Control Letters 11, Oct. 1988, pp 167-172.
# [4] P. A. Iglesias and K. Glover, "State-space approach to discrete-time
# H-infinity control." Int. J. Control, 1991, V. 54, #5, 1031-1073.
#
# Written by A. S. Hodel
# Updated by John Ingram July 1996 to accept structured systems
#
# Revised by Kai P Mueller April 1998 to solve the general H_infinity
# problem using unitary transformations Q (on w and z)
# and non-singular transformations R (on u and y) such
# that the Dzu and Dyw matrices of the transformed plant
#
# ~
# P (the variable Asys here)
#
# become
#
# ~ -1 T
# D = Q D R = [ 0 I ] or [ I ],
# 12 12 12 12
#
# ~ T
# D = R D Q = [ 0 I ] or [ I ].
# 21 21 21 21
#
# This transformation together with the algorithm in [1] solves
# the general problem (see [2] for example).
if (nargin < 3) | (nargin > 4)
usage("[retval,dgkf_struct] = is_dgkf(Asys,nu,ny{,tol})");
elseif (! is_scal(nu) | ! is_scal(ny) )
error("is_dgkf: arguments 2 and 3 must be scalars")
elseif (! is_struct(Asys) )
error("Argument 1 must be a system data structure");
endif
if(nargin < 4)
tol = 200*eps;
elseif( !is_sampl(tol) )
error("is_dgkf: tol must be a positive scalar")
endif
retval = 1; # assume passes test
dflg = is_digit(Asys);
[Anc, Anz, nin, nout ] = sysdimen(Asys);
if( Anz == 0 & Anc == 0 )
error("is_dgkf: no system states");
elseif( nu >= nin )
error("is_dgkf: insufficient number of disturbance inputs");
elseif( ny >= nout )
error("is_dgkf: insufficient number of regulated outputs");
endif
nw = nin - nu; nw1 = nw + 1;
nz = nout - ny; nz1 = nz + 1;
[A,B,C,D] = sys2ss(Asys);
# scale input/output for numerical reasons
if(norm(C,'fro')*norm(B,'fro') == 0)
error("||C||*||B|| = 0; no dynamic connnection from inputs to outputs");
endif
xx = sqrt(norm(B, Inf) / norm(C, Inf));
B = B / xx; C = C * xx;
# partition matrices
Bw = B(:,1:nw); Bu = B(:,nw1:nin);
Cz = C(1:nz,:); Dzw = D(1:nz,1:nw); Dzu = D(1:nz,nw1:nin);
Cy = C(nz1:nout,:); Dyw = D(nz1:nout,1:nw); Dyu = D(nz1:nout,nw1:nin);
# Check for loopo shifting
Dyu_nz = (norm(Dyu,Inf) != 0);
if (Dyu_nz)
warning("is_dgkf: D22 nonzero; performing loop shifting");
endif
# 12 - rank condition at w = 0
xx =[A, Bu; Cz, Dzu];
[nr, nc] = size(xx);
irank = rank(xx);
if (irank != nc)
retval = 0;
warning(sprintf("rank([A Bu; Cz Dzu]) = %d, need %d; n=%d, nz=%d, nu=%d", ...
irank,nc,(Anc+Anz),nz,nu));
warning(" *** 12-rank condition violated at w = 0.");
endif
# 21 - rank condition at w = 0
xx =[A, Bw; Cy, Dyw];
[nr, nc] = size(xx);
irank = rank(xx);
if (irank != nr)
retval = 0;
warning(sprintf("rank([A Bw; Cy Dyw]) = %d, need %d; n=%d, ny=%d, nw=%d", ...
irank,nr,(Anc+Anz),ny,nw));
warning(" *** 21-rank condition violated at w = 0.");
endif
# can Dzu be transformed to become [0 I]' or [I]?
# This ensures a normalized weight
[Qz, Ru] = qr(Dzu);
irank = rank(Ru);
if (irank != nu)
retval = 0;
warning(sprintf("*** rank(Dzu(%d x %d) = %d", nz, nu, irank));
warning(" *** Dzu does not have full column rank.");
endif
if (nu >= nz)
Qz = Qz(:,1:nu)';
else
Qz = [Qz(:,(nu+1):nz), Qz(:,1:nu)]';
endif
Ru = Ru(1:nu,:);
# can Dyw be transformed to become [0 I] or [I]?
# This ensures a normalized weight
[Qw, Ry] = qr(Dyw');
irank = rank(Ry);
if (irank != ny)
retval = 0;
warning(sprintf("*** rank(Dyw(%d x %d) = %d", ny, nw, irank));
warning(" *** Dyw does not have full row rank.");
endif
if (ny >= nw)
Qw = Qw(:,1:ny);
else
Qw = [Qw(:,(ny+1):nw), Qw(:,1:ny)];
endif
Ry = Ry(1:ny,:)';
# transform P by Qz/Ru and Qw/Ry
Bw = Bw*Qw;
Bu = Bu/Ru;
B = [Bw, Bu];
Cz = Qz*Cz;
Cy = Ry\Cy;
C = [Cz; Cy];
Dzw = Qz*Dzw*Qw;
Dzu = Qz*Dzu/Ru;
Dyw = Ry\Dyw*Qw;
# pack the return structure
dgkf_struct.nw = nw;
dgkf_struct.nu = nu;
dgkf_struct.nz = nz;
dgkf_struct.ny = ny;
dgkf_struct.A = A;
dgkf_struct.Bw = Bw;
dgkf_struct.Bu = Bu;
dgkf_struct.Cz = Cz;
dgkf_struct.Cy = Cy;
dgkf_struct.Dzw = Dzw;
dgkf_struct.Dzu = Dzu;
dgkf_struct.Dyw = Dyw;
dgkf_struct.Dyu = Dyu;
dgkf_struct.Ru = Ru;
dgkf_struct.Ry = Ry;
dgkf_struct.Dyu_nz = Dyu_nz;
dgkf_struct.dflg = dflg;
endfunction