% [Z, U, X] = izde(G,EPP); % % Inputs: G - system matrix in mu-tools format. % EPP - tolerance, see below, default value EPS. % Outputs: Z - zeros. % U - input zero directions, stored as column vectors. % X - state directions, stored as column vectors. % Each column X(:,i) and U(:,i) corresponds to the zeros Z(i). % % This is a modifaction of szeros.m written by Kjetil Havre % % This Input-Zero-Direction-"through genaralized Eigenvalue" decomposition % IZDE function is a modification of the szeros function contained in mu % - toolbox. The modification consists of returning the input zero % directions and the state zero directions in addition to the zeros. The % input zero directions u are defined as: % G(z)*u = 0, % where s = z is a zero of G(s). This is done by solving the generalized % eigenvalue problem: % % | A-Iz | B | * | x | = | 0 | % |------------| |---| |---| % | C | D | | u | | 0 | % % IZDE finds the transmission zeros z of a SYSTEM matrix. Occasionally, % large zeros are included which may actually be at infinity. Solving % for the transmission zeros of a system involves two generalized eigen- % value problems. EPP (optional) defines if the difference between two % generalize eigenvalues is small. IZDE also finds the input u and the % state x directions of the zeros. % % The input zero directions are stored as column vectors in U, and each % of the columns are normalized. The state zero directions are stored as % columns in X. The degree of freedom to normalize the generalized eigen- % vector is used to normalize the u part. So the length of x is not equal % to one. Each column in U and X corresponds to the element in Z with % same place. % % For systems with more inputs than outputs the input zero direction % is not a complete basis for the nullspace of G(z). % Zeros with multiplicity greater than one (rare cases which may % occure in non-minimal realizations), may (not sure) cause wrong % directions. % % Comments, corrections and malfunctions, can be e-mailed to: % havre@kjemi.unit.no or skoge@kjemi.unit.no % % See also: EIG, SZEROS, OZDE and SPOLES. % Algorithm based on Laub & Moore 1978 paper, Automatica % % Note that when the number of inputs is larger than the number of outputs, % the input zero direction is not complete. A bit clearer: If z is a zero of % a non-square plant with number of inputs greater than number of outputs, % the zero direction is not a line but a surface. As an example, consider % G(s) with dimensions 2x3 (2 rows and 3 columns). Let s=z be a zero of % G(s) such that G(z)*u = [0;0]; Since s=z is a zero then the rank of % G(z) has to be less than the normal rank of G(s), which at maximum can % be 2. This implies that rank of G(z) must be less than 2. % u is element in the three dimensional field of real numbers. Since the % rank of G(z) is maximum one the zero direction is a actually a subspace % in this three dimensional field of real numbers given by two basis % vectors. Since this function only gives one zero direction for a given % zero this direction does not describe the input zero space completely. % Two basis vectors are requiered. % % Modification for square systems was made by: Kjetil Havre 14/5-1995. % Modification for non square systems was made by: Kjetil Havre 14/5-1995. % Inclusion of state directions was made by: Kjetil Havre 14/5-1995. % Modified so that first element of U(:,i) is real: Kjetil Havre 3/2-1996. % Copyright 1996-2003 Sigurd Skogestad & Ian Postlethwaite % $Id: izde.m,v 1.2 2004/01/19 14:52:11 aske Exp $ function [Z, U, X] = izde(sys,epp) if nargin < 1 disp('usage: [Z,Y,X] = ozde(G) ') return end if nargin == 1 epp = eps; end [ny,nu,nx]=size(sys); if class(sys) == 'tf'|'ss'|'zpk'|'frd' [a,b,c,d] = ssdata(sys); if nx == 0 disp('SYSTEM has no states') end sysu = [a b; c d]; % find generalized eigenvalues of a square system matrix if ny == nu x = zeros(nx+nu,nx+nu); x(1:nx,1:nx) = eye(nx); [vech, ev] = eig(sysu, x); z = diag(ev); % Extract the eigenvalues. kc=0; % Counter for eigenvalues. for k=1:max(size(z)), logic = ~isnan(z(k)) & finite(z(k)); if logic kc= kc+1; Z(kc,1) = z(k); vech2(:,kc) = vech(:,k); end end % Split x and u. vx = vech2(1:nx, : ); vu = vech2(nx+1:nx+nu,:); % Normalize columns. [nvr, nvc] = size( vu ); for i=1:nvc, nrmu = norm(vu(:,i)); if nrmu > 1000*epp vx(:,i) = vx(:,i)/nrmu; vu(:,i) = vu(:,i)/nrmu; else vu(:,i) = zeros(nu,1); end Inz = find( abs(vu(:,i)) > 1000*epp ); if isempty(Inz) == 0 U(:,i) = vu(:,i) * exp( -angle(vu(Inz(1),i))*sqrt(-1) ); X(:,i) = vx(:,i) * exp( -angle(vu(Inz(1),i))*sqrt(-1) ); else X(:,i) = vx(:,i); Y(:,i) = vu(:,i); end end else % Non-square systems nrm = norm(sysu,1); if nu < ny x1 = [ sysu nrm*(rand(nx+ny,ny-nu)-.5)]; x2 = [ sysu nrm*(rand(nx+ny,ny-nu)-.5)]; else x1 = [ sysu; nrm*(rand(nu-ny,nx+nu)-.5)]; x2 = [ sysu; nrm*(rand(nu-ny,nx+nu)-.5)]; end [x]= zeros(size(x1)); x(1:nx,1:nx) = eye(nx); [v1h z1h] = eig(x1,x); % Compute the genaralized eigenvalues [v2h z2h] = eig(x2,x); % for the two augumented systems. z1h2 = diag( z1h ); z2h2 = diag( z2h ); z2 = z2h2(~isnan(z2h2) & finite(z2h2)); kc=0; % Counter for eigenvalues. for k=1:max(size(z1h2)), logic = ~isnan(z1h2(k)) & finite(z1h2(k)); if logic kc= kc+1; z1(kc,1) = z1h2(k); vech2(:,kc) = v1h(:,k); end end nz = length(z1); vech3 = []; Z = []; for i=1:nz, if any(abs(z1(i)-z2) < nrm*sqrt(epp)) Z = [Z; z1(i)]; vech3 = [vech3 vech2(:,i)]; end end % Split in ux and xz if isempty( vech3 ) Z = []; U = []; X = []; return; end vx = vech3(1:nx, : ); vu = vech3(nx+1:nx+nu,:); % Normalize columns. [nvr, nvc] = size( vu ); for i=1:nvc, nrmu = norm(vu(:,i)); if nrmu > 1000*epp vx(:,i) = vx(:,i)/nrmu; vu(:,i) = vu(:,i)/nrmu; else vu(:,i) = zeros(nu,1); end Inz = find( abs(vu(:,i)) > 1000*epp ); if isempty(Inz) == 0 U(:,i) = vu(:,i) * exp( -angle(vu(Inz(1),i))*sqrt(-1) ); X(:,i) = vx(:,i) * exp( -angle(vu(Inz(1),i))*sqrt(-1) ); else X(:,i) = vx(:,i); Y(:,i) = vu(:,i); end end end else error('matrix is not a SYSTEM matrix') return end % % Copyright MUSYN INC 1991, All Rights Reserved % Copyright MUSYN INC 1995, All Rights Reserved %