CompRGA.m000064400136200001130000000030150575143002600131260ustar00kjetilentek00000000000000% RGA = CompRGA( Gjw, outp, inp )
%           or
% CompRGA( Gjw, outp, inp )
% 
% Computation of frequency dependent relative gain array.
% 
% Outputs: RGA  - relative gain array.
% If no return arguments are supplied (equivalent to second calling
% syntax) plots of RGA elements in each row are given.
%
% Inputs:  Gjw  - varying matrix,
%          outp - vector which describes which rows to select,
%                 if omitted all rows are selected.
%          inp  - vector which describes which columns to select,
%                 if omitted all columns are selected.
%
function RGA = CompRGA( m, outp, inp )
   [mt, rd, cd, n] = minfo( m );

   if (strcmp( mt, 'vary' ) == 1) & ( max(size(outp)) > 1 )
   
      if nargin < 3
         inp = ':';
      end
      if nargin < 2
         outp = ':';
      end
      
      if nargin < 1
         disp( 'Usage: RGA = CompRGA( m, outp, inp )' )
      end

      ms = sel(m, outp, inp);
      ms_inv_trsp = vinv( ms );
      ms_inv_trsp = transp( ms_inv_trsp );
      RGA = veval( '.*', ms, ms_inv_trsp );

      % Plot results.
      if nargout < 1
         for i=1:max(size(outp))
            % plot some of the results.
            figure( i );
   	    iv = getiv( ms );
	    vplot( 'liv,m', sel(RGA, i, ':') );	grid; 
	    title( [ 'Frequency dependent RGA for output ' int2str(i) '.' ] );
            xlabel ('Frequency [rad/s].');
            ylabel ('Magnitude');
            v = axis;
	    v(1) = iv(1);
	    v(2) = iv(max(size(iv)));
	    axis( v );
	 end
      end
      
   end
   
end
, ms, ms_inv_trsp );

      % Plot results.
      if nargout < 1
         for i=1:max(size(outp))
            % plot some of the results.
            figure( i );
   	    iv = getiv( ms );
	    vplot( 'liv,m', sel(RGA, i, ':') );	grid; 
	    title( [ 'Frequency dependent RGA for output ' int2str(i) '.' ] );
            xlabel ('Frequency [rad/s].');
            ylabel ('Magnitude');
            v = axis;
	    v(1) = iv(1);
	    v(2) = iv(max(size(iv)));
	    axis( v );
	 end
      end
      
  CompSingCond.m000064400136200001130000000033070575142012500142240ustar00kjetilentek00000000000000% function [s, c, u, v] = CompSingCond( Gjw, outp, inp )
% 
% Computation of singular value decomposition and condition 
% number for a given VARYING or Constant matrix.
% 
% Outputs: s    - singular values,
%          c    - condition number,
%          u    - output singular vectors,
%          v    - input singular vectors.
% G = u*s*v';
%
% Inputs:  Gjw  - varying matrix,
%          outp - vector which describes which rows to select,
%                 if omitted all rows are selected.
%          inp  - vector which describes which columns to select,
%                 if omitted all columns are selected.


function [sin, con, u, v] = CompSingCond( G_jw, outp, inp )
   [mt, rd, cd, n] = minfo( G_jw );
   if strcmp( mt, 'vary' ) == 1

      if nargin < 3
         inp = ':';
      end
      if nargin < 2
         outp = ':';
      end
      
      if nargin < 1
         disp( 'Usage: [s, c, u, v] = CompSingCond( Gjw, outp, inp )' )
      end
      
      [u, sin, v] = vsvd( sel(G_jw, outp, inp) );
      con = vcond( sel(G_jw, outp, inp) );

      % Plot the results if no return variables.
      if nargout < 1
         figure( 1 );
         vplot( 'liv,lm', sin ); grid;
	 v = axis;
	 iv = getiv( G_jw );
	 v(1) = iv(1);
	 v(2) = iv( max(size(iv)) );
	 axis( v );
	 title ( 'Singular values as function of frequency.');
         xlabel ('Frequency [rad/s].');
         ylabel ('Magnitude');

	 if max(size(outp)) > 1
  	    figure( 2 );
	    vplot( 'liv,m', con ); grid;
	    v = axis;
	    v(1) = iv(1);
	    v(2) = iv( max(size(iv)) );
	    axis( v );
	    title ( 'Condition number as function of frequency.');
	    xlabel ('Frequency [rad/s].');
	    ylabel ('Magnitude');
         end
      end
      
   end
end
 v(1) = iv(1);
	 v(2) = iv( max(size(iv)) );
	 axis( v );
	 title ( 'Singular values as function of frequency.');
         xlabel ('Frequency [rad/s].');
         ylabel ('Magnitude');

	 if max(size(outp)) > 1
  	    figure( 2 );
	    vplot( 'liv,m', con ); grid;
	    v = axis;
	    v(1) = iv(1);
	    v(2) = ivSubSysZero.m000064400136200001130000000010530577723714600140060ustar00kjetilentek00000000000000function Z = SubSysZero(G,outp,inp,tol)
   if nargin < 3
      inp = ':';
   end
   if nargin < 2
      outp = ':';
   end
   if nargin < 1
     disp('Usage: SubSysZero(sys,outp, inp)')
     return
   end
    
   [mt, rd, cd, n] = minfo( G );
   if strcmp( mt, 'syst' ) == 0
     disp('Usage: SubSysZero(sys,outp, inp)')
     return      
   end

   Gs = sel(G,outp,inp);
   [A,B,C,D] = unpck(Gs);
   if nargin < 4
      tol = 10*max(size(A))*norm(A,1)*eps;
   end
   [Am,Bm,Cm,Dm] = minreal(A,B,C,D,tol);
   Gm = pck(Am,Bm,Cm,Dm);
   Z = szeros(Gm);
end
if nargin < 3
      inp = ':';
   end
   if nargin < 2
      outp = ':';
   end
   if nargin < 1
     disp('Usage: SubSysZero(sys,outp, inp)')
     return
   end
    
   [mt, rd, cd, n] = minfo( G );
   if strcmp( mt, 'syst' ) == 0
     disp('Usage: SubSysZero(sys,outp, inp)')
     return      
   end

   Gs = sel(G,outp,inp);
   [A,B,C,D] = unpck(Gs);
   if nargin < 4
      tol = 10*max(size(A))*norm(A,1)*eps;
   end
   [Am,Bm,Cm,Dm] = minreal(A,B,C,D,tol);
   Gm allbpz.m000064400136200001130000000037550605734014400131760ustar00kjetilentek00000000000000%   function [abpz, seq] = allbpz(yp,p,z,lv)
%
%   The ``All-Bp(z)'' allbpz-function calculates all possible filters
%   Bp(z) as described in Havre and Skogestad (1995) by manipulating 
%   the sequence of factorization of the indvidual filters.
%   Note! that the function is recursive.
%   
%   Inputs:  yp   - output directions of RHP-poles stored as a column 
%                   vectors, yp(:,i) corresponds to pole p(i).
%            p    - poles stores in a singel vector p, p(i) is expected 
%                   to corrrespond to yp(:,i).
%            z    - single zero where Bp(z).
%            lv   - Input argument only used by the function itself.
% 
%   Outputs: abpz - Bp(z) given as one large matrix of size (ny x np!*ny)
%                   where ny - number outputs, np - numper of poles in p.
%            seq  - vector with variables describing the sequence of 
%                   factorizations of the poles p(i).
%
%   Written by: Kjetil Havre, 24/11-1995, e-mail: kjetil@ife.no
%
%   Reference: Havre K. and S. Skogestad, 1995, ``Effect of RHP Zeros and
%              Poles on Performance in Multivariable Systems''.
%

function [abpz, seq] = allbpz(yp,p,z,lv)

   [n, np] = size(p);
   if n > 1
      p = vtp(p);
      [n, np] = size(p);
   end

   if isempty(z)
      disp('Error: z is empty.');
      return
   end
   
   [nr,nc] = size(yp);
   if nc ~= np
      disp(['Error: Number of poles:', int2num(np),...
            ' is not equal to number of output zero directions: ',...
	    int2num(nc)] );
      return
   end

   if nargin < 4
      lv = 1:np;
   end

   if np==1;
      abpz = eye(nr) + 2*real(p)/(z-p)*yp*yp';
      seq = lv(1);
      return;
   end
   
   abpz = []; seq = [];
   sif = 10^(np-1);
   for i=1:np
      bi      = eye(nr) + 2*real(p(i))/(z-p(i))*yp(:,i)*yp(:,i)';
      [Bl sl] = allbpz( [yp(:,1:i-1) yp(:,i+1:np)], [p(1:i-1) p(i+1:np)], z, [lv(1:i-1) lv(i+1:np)] );
      B    = bi*Bl;
      abpz = [abpz B];
      seq  = [seq lv(i)*sif*ones(size(sl))+sl];
   end
end

nt2num(nc)] );
    allbzp.m000064400136200001130000000040400605763751600131760ustar00kjetilentek00000000000000%   function [abzp, seq] = allbzp(yz,z,p,lv)
%
%   The ``All-Bz(p)'' allbzp-function calculates all possible filters
%   Bz(p) as described in Havre and Skogestad (1995) by manipulating 
%   the sequence of factorization of the indvidual filters.
%   Note! that the function is recursive.
%   
%   Inputs:  yz   - output directions of RHP-zeros stored as a column 
%                   vectors, yz(:,i) corresponds to zero z(i).
%            z    - zeros stores in a singel vector z, z(i) is expected 
%                   to corrrespond to yz(:,i).
%            p    - single pole where Bz(p).
%            lv   - Input argument only used by the function itself.
% 
%   Outputs: abzp - Bz(p) given as one large matrix of size (nz!*ny x ny)
%                   where ny - number outputs, nz - numper of zeros in z.
%            seq  - vector with variables describing the sequence of 
%                   factorizations of the zeros z(i).
%
%   Written by: Kjetil Havre, 30/11-1995, e-mail: kjetil@ife.no
%
%   Reference: Havre K. and S. Skogestad, 1995, ``Effect of RHP Zeros and
%              Poles on Performance in Multivariable Systems''.
%

function [abzp, seq] = allbzp(yz,z,p,lv)

   [n, nz] = size(z);
   if n > 1
      z = vtp(z);
      [n, nz] = size(z);
   end

   if isempty(p)
      disp('Error: p is empty.');
      return
   end
   
   [nr,nc] = size(yz);
   if nc ~= nz
      disp(['Error: Number of zeros:', int2num(nz),...
            ' is not equal to number of output zero directions: ',...
	    int2num(nc)] );
      return
   end

   if nargin < 4
      lv = 1:nz;
   end

   if nz==1;
      abzp = eye(nr) + 2*real(z)/(p-z)*yz*yz';
      seq = lv(1);
      return;
   end
   
   abzp = []; seq = [];
   sif = 10^(nz-1);
   for i=1:nz
      bi      = eye(nr) + 2*real(z(i))/(p-z(i))*yz(:,i)*yz(:,i)';
      [Bl sl] = allbzp( [yz(:,1:i-1) yz(:,i+1:nz)], [z(1:i-1) z(i+1:nz)],...
                                                p, [lv(1:i-1) lv(i+1:nz)] );
      B    = bi*Bl;
      abzp = [abzp B];
      seq  = [seq lv(i)*sif*ones(size(sl))+sl];
   end
end

4
      lv = 1:nz;
   end

   if nz==1;
      abzp = eye(nr) + 2*real(z)/(p-z)*yz*yz';
      seq = lv(1);
      return;
   end
   
   abzp = []; seq = [];
   sif = 10^(nz-1);
   for i=1:nz
      bi      = eye(nr) + 2*real(z(i))/(p-z(i))*yz(:,i)*yz(:,i)';
      [Bl sl] = allbzp( [yz(:,1:i-1) yz(:,i+1:nz)], [z(1:i-1) z(i+1:nz)],...
                                                p, [lv(1:i-1) lv(i+1:nz)] );
      B    = bi*Bl;
      abzp = [abzp B];
      seq  = [seq lv(i)*sif*condmin.m000064400136200001130000000010750603624463200133340ustar00kjetilentek00000000000000
% Computes the minimized condition number using the upper mu-bound
% S.Skogestad, 14May95

function gmin = condmin(Gf)

[dum1,n,m,dum2]=minfo(Gf);
zero = 0*eye(n);
H = abv( sbs(zero, minv(Gf)), sbs(Gf, zero));

blk0=[1 1]; blk1=blk0;
for k=2:n
  blk1 = [blk1; blk0];
end
blk2=blk1; blk=[blk1; blk2];

mu1=mu(H,blk,'Uc'); gmin1=sel(mu1,':',1); 
gmin=mmult(gmin1,gmin1); 

% Remark: The upper mu-bound with the present software is not always
% reliable. Andy Packard has a better bound for the case when all blocks
% are 1x1 complex:     [d,sclmat,gmin1] = sbscl(H,1.e-3)


-bound
% S.Skogestad, 14May95

function gmin = condmin(Gf)

[dum1,n,m,dum2]=minfo(Gf);
zero = 0*eye(n);
H = abv( sbs(zero, minv(Gf)), sbs(Gf, zero));

blk0=[1 1]; blk1=blk0;
for k=2:n
  blk1 = [blk1; blk0];
end
blk2=blk1; blk=[blk1; blk2];

mu1=mu(H,blk,'Uc'); gmin1=sel(mu1,':',1); 
gmin=mmult(gmin1,gmin1); 

% Remark: The upper mu-bound with the present software is not always
% reliable. Andy Packard has a better bound for the case when all blockcondmini.m000064400136200001130000000004220603624464100135000ustar00kjetilentek00000000000000
function gmin = condmini(Gf)

[dum1,n,m,dum2]=minfo(Gf);
zero = 0*eye(n);
H = abv( sbs(zero, minv(Gf)), sbs(Gf, zero));

blk0=[1 1]; blk1=blk0;
for k=2:n
  blk1 = [blk1; blk0];
end
blk2=[n n]; blk=[blk1; blk2];

gmin1=sel(mu(H,blk,'Uc'),':',1); gmin=mmult(gmin1,gmin1); 


k2=blk1; blk=[blk1; blk2];

mu1=mu(H,blk,'Uc'); gmin1=sel(mu1,':',1); 
gmin=mmult(gmin1,gmin1); 

% Remark: The upper mu-bound with the present software is not always
% reliable. Andy Packard has a better bound for the case when all blockcondmino.m000064400136200001130000000004220603624465100135070ustar00kjetilentek00000000000000
function gmin = condmino(Gf)

[dum1,n,m,dum2]=minfo(Gf);
zero = 0*eye(n);
H = abv( sbs(zero, minv(Gf)), sbs(Gf, zero));

blk0=[1 1]; blk1=blk0;
for k=2:n
  blk1 = [blk1; blk0];
end
blk2=[n n]; blk=[blk2; blk1];

gmin1=sel(mu(H,blk,'Uc'),':',1); gmin=mmult(gmin1,gmin1); 


k2=blk1; blk=[blk1; blk2];

mu1=mu(H,blk,'Uc'); gmin1=sel(mu1,':',1); 
gmin=mmult(gmin1,gmin1); 

% Remark: The upper mu-bound with the present software is not always
% reliable. Andy Packard has a better bound for the case when all blockcoprimeunc.m000064400136200001130000000036450605213471600140550ustar00kjetilentek00000000000000function [K, X, Z, gamopt, A1, R1, Q1] = coprimeunc(Gs, gammarel)
%
% Finds the controller which optimally robustifies a given
% shaped plant in terms of tolerating maximum coprime uncertainty.
% Inputs:
%   Gs       - System containig the shaped plant.
%   gammarel - Final gamma used is gamma is: gamme = gammarel*gamopt
% Outputs:
%   K        - "Robustifying" controller assuming positive feedback.
%   X        - Solution to the ricatti eqn:
%                             A1'*X + X*A1 
%                               - X*R1*X          + Q1          == 0 
%              (A-B*inv(S)*D'*C)'*X + X*(A-B*inv(S)*D'*C) 
%                               - X*B*inv(S)*B'*X + C'*inv(R)*C == 0
%   Z        - Solution to the ricatti eqn:
%                             A1*X + X*A1' 
%                               - X*Q1*X          + R1          == 0 
%              (A-B*inv(S)*D'*C)*Z + Z*(A-B*inv(S)*D'*C)' 
%                               - Z*C'*inv(R)*C*Z + B*inv(S)*B' == 0
%
%
%

% find normalized coprime factors of the shaped plant

[A,B,C,D] = unpck(Gs);
S = eye(size(D'*D)) + D'*D; R = eye(size(D*D')) + D*D';
Rinv = inv(R); Sinv = inv(S);

A1 = A-B*Sinv*D'*C; Q1 = C'*Rinv*C; R1 = B*Sinv*B';
ham1 = [A1 R1;Q1 -A1'];
[X1,X2,fail,reig_min] = ric_schr(ham1);
if (fail>=1)
   disp('*** WARNING! X-Riccati Equation may be ill-posed ***')
else
   X = X2/X1;
end

A2 = A1'; Q2 = R1; R2 = Q1;
ham2 = [A2 R2; Q2 -A2']; 
[Z1,Z2,fail,reig_min] = ric_schr(ham2);
if (fail>=1)
   disp('*** WARNING! Z-Riccati Equation may be ill-posed ***')
else
   Z = Z2/Z1;
end

XZ = X*Z; gamopt=sqrt(1+max(eig(XZ)));
disp(['Optimal/minimal gamma value: ', num2str(gamopt), '.']);

% Use higher gamma, gamma = gammarel*gamopt;
gam = gammarel*gamopt; gam2 = gam*gam;  gamconst = (1-gam2)*eye(size(XZ));
Lc  = gamconst + XZ;   Li   = inv(Lc'); Fc       = -Sinv*(D'*C+B'*X);
Ac  = A+B*Fc + gam2*Li*Z*C'*(C+D*Fc);
Bc = gam2*Li*Z*C'; Cc = B'*X; Dc = -D';

K = pck(Ac,Bc,Cc,Dc);




_schr(ham2);
if (fail>=1)
   disp('*** WARNING! Z-Riccati Equation may be ill-posed ***')
ecwcc1.m000064400136200001130000000073050605754224200127100ustar00kjetilentek00000000000000%   function [c1, z, seq, fc1, fseq] = cwcc1(G, ztol, pst, zst)
%
%   The ``Compute-Worst-Case-c1'' cwcc1-function maximize the
%   c1 value by considering all possible factorizations of 
%   RHP-poles in Bp(s) and all RHP-zeros.
%
%   Inputs:  G    - System containing RHP-zeros and poles.
%            ztol - Upper value on zeros which are treated as 
%                   RHP-zeros. Typical 1e12, default value 
%                   1e12.
%            pst  - SST argument passed to PDFST.
%            zst  - SST argument passed to ZDSVD.
%
%   Outputs: c1   - Worst case c1 value.
%            z    - RHP-zero for which c1 occure.
%            seq  - Describes the sequence of factorization for c1.
%            fc1  - All possible values of c1. Matrix containing
%                   dim(Zrhp)-rows and fact{dim(Prhp)} columns.
%                   Element fc1(i,j) is given for RHP-zero Zrhp(i)
%                   with factorization described by fseq(j).
%            fseq - Describes the sequence of factorization for the
%                   columns in fc1.
%            Zrhp - RHP-zeros.
%            Prhp - RHP-poles.
%
%   Written by: Kjetil Havre, 30/11-1995, e-mail: kjetil@ife.no
%
%   See also:  ZDSVD, PDFST ALLBPZ and CWCC2.
%
%   Reference: Havre K. and S. Skogestad, 1995, ``Effect of RHP Zeros and
%              Poles on Performance in Multivariable Systems''.
%
function [c1, z, seq, fc1, fseq, Zrhp, Prhp] = cwcc1(G, ztol, pst, zst)

   c1 = []; z = []; seq = []; fc1 = []; fseq = []; Zrhp = []; Prhp = [];
   if nargin < 4
      zst = 1e-14;
   end
   if nargin < 3
      pst = 1e-12;
   end
   if nargin < 2
      ztol = 1e12;
   end
   if nargin < 1
      disp('Usage: [c1, z, seq, fc1, fseq] = cwcc1(G, ztol, pst, zst)');
      return
   end

   [mt,ny,nu,nx] = minfo(G);

   if strcmp(mt, 'syst')==0
      disp( 'System matrix G is required, usage: [I,O] = cwcc1(G, ztol)' );
      return
   end

   Z = szeros(G);                        % Zeros.
   P = spoles(G);                        % Poles.
   IPrhp = find(P>0);                    % Finds indices for poles in RHP 
   IZrhp = find(real(Z)>0&abs(Z)<ztol);  % Finds indices for zeros in RHP
   Prhp  = P(IPrhp);                     % Store RHP-poles.
   Zrhp  = Z(IZrhp);                     % Store RHP-zeros.
   np = max(size(IPrhp));                % Number of RHP-poles.
   nz = max(size(IZrhp));                % Number of RHP-zeros.

   if np == 0
      disp('Be happy: No RHP-poles.');
      return;
   end
   if nz == 0
      disp('Be happy: No RHP-zeros.');
      return;
   end
   
   [UP, YP, erp] = pdfst(G, Prhp, pst);
   if erp 
      if erp > 0
         disp('Error: In pole-direction calculation PDFST.')
         disp('       No direction with gain larger than PST')
         disp('       No direction with gain larger than PST')	 
         disp('       Decrease PST argument.')
	 return
      else
         disp('Warning: Two directions are possible pole directions.')
         disp('         May try to increase the PST argument.')
      end
   end
   
   [UZ, YZ, erz] = zdsvd(G, Zrhp, zst);
   if erz
      if erz > 0
         disp('Error: In zero-direction calculation zdsvd.')
         disp('       No direction with gain smaller than ZST.')
         disp('       Increace ZST argument.')
         return;
      else
         disp('Warning: Two directions are possible pole directions.')      
         disp('         May try to decrease the ZST argument.')
      end
   end

   for j=1:nz
      [Bpzs fseq]  = allbpz(YP,Prhp,Zrhp(j));
      yzBpz = YZ(:,j)'*Bpzs;
      for i=1:fact(np)
         fc1(j,i) = norm(yzBpz(1,ny*(i-1)+1:ny*i));
      end
   end

   c1         = max(max(fc1));
   [Ic1, Jc1] = find(fc1==c1);
   seq        = fseq(Jc1);
   z          = Zrhp(Ic1);

end
reace ZST argument.')
         return;
      else
         disp('Warning: Two directions are possible pole directions.')      
         disp('         May try to decrease the ZST argument.')
      end
   end

   for j=1:nz
      [Bpzs fseq]  = allbpz(YP,Prhp,Zrhp(j));
      yzBpz = YZ(:,j)'*Bpzs;
      for i=1:faccwcc2.m000064400136200001130000000073050605763777200127260ustar00kjetilentek00000000000000%   [c2, p, seq, fc2, fseq, Zrhp, Prhp] = cwcc2(G, ztol, pst, zst)
%
%   The ``Compute-Worst-Case-c2'' cwcc2-function maximize the
%   c2 value by considering all possible factorizations of 
%   RHP-poles in Bz(s) and all RHP-poles.
%
%   Inputs:  G    - System containing RHP-zeros and poles.
%            ztol - Upper value on zeros which are treated as 
%                   RHP-zeros. Typical 1e12, default value 
%                   1e12.
%            pst  - SST argument passed to PDFST.
%            zst  - SST argument passed to ZDSVD.
%
%   Outputs: c2   - Worst case c2 value.
%            p    - RHP-pole for which c2 occure.
%            seq  - Describes the sequence of factorization for c2.
%            fc2  - All possible values of c2. Matrix containing
%                   dim(Prhp)-rows and fact{dim(Zrhp)} columns.
%                   Element fc2(i,j) is given for RHP-pole Prhp(i)
%                   with factorization described by fseq(j).
%            fseq - Describes the sequence of factorization for the
%                   columns in fc2.
%            Zrhp - RHP-zeros.
%            Prhp - RHP-poles.
%
%   Written by: Kjetil Havre, 30/11-1995, e-mail: kjetil@ife.no
%
%   See also:  ZDSVD, PDFST ALLBZP and CWCC1.
%
%   Reference: Havre K. and S. Skogestad, 1995, ``Effect of RHP Zeros and
%              Poles on Performance in Multivariable Systems''.
%
function [c2, p, seq, fc2, fseq, Zrhp, Prhp] = cwcc2(G, ztol, pst, zst)

   c2 = []; p = []; seq = []; fc2 = []; fseq = []; Zrhp = []; Prhp = [];
   if nargin < 4
      zst = 1e-14;
   end
   if nargin < 3
      pst = 1e-12;
   end
   if nargin < 2
      ztol = 1e12;
   end
   if nargin < 1
      disp('Usage: [c2, p, seq, fc2, fseq] = cwcc2(G, ztol, pst, zst)');
      return
   end

   [mt,ny,nu,nx] = minfo(G);

   if strcmp(mt, 'syst')==0
      disp( 'System matrix G is required, usage: [c2, p, seq, fc2, fseq] = cwcc2(G, ztol, pst, zst)' );
      return
   end

   Z = szeros(G);                   % Zeros.
   P = spoles(G);                   % Poles.
   IPrhp = find(P>0);               % Finds indices for poles in RHP 
   IZrhp = find(real(Z)>0&abs(Z)<ztol);  % Finds indices for zeros in RHP
   Prhp  = P(IPrhp);                % Store RHP-poles.
   Zrhp  = Z(IZrhp);                % Store RHP-zeros.
   np = max(size(IPrhp));           % Number of RHP-poles.
   nz = max(size(IZrhp));           % Number of RHP-zeros.

   if np == 0
      disp('Be happy: No RHP-poles.');
      return;
   end
   if nz == 0
      disp('Be happy: No RHP-zeros.');
      return;
   end
   
   [UP, YP, erp] = pdfst(G, Prhp, pst);
   if erp 
      if erp > 0
         disp('Error: In pole-direction calculation PDFST.')
         disp('       No direction with gain larger than PST')
         disp('       No direction with gain larger than PST')	 
         disp('       Decrease PST argument.')
	 return
      else
         disp('Warning: Two directions are possible pole directions.')
         disp('         May try to increase the PST argument.')
      end
   end
   
   [UZ, YZ, erz] = zdsvd(G, Zrhp, zst);
   if erz
      if erz > 0
         disp('Error: In zero-direction calculation zdsvd.')
         disp('       No direction with gain smaller than ZST.')
         disp('       Increace ZST argument.')
         return;
      else
         disp('Warning: Two directions are possible pole directions.')      
         disp('         May try to decrease the ZST argument.')
      end
   end

   for j=1:np
      [Bzps fseq]  = allbzp(YZ,Zrhp,Prhp(j));
      for i=1:fact(nz)
         Bzpyp = Bzps(:,ny*(i-1)+1:ny*i)*YP(:,j);
         fc2(j,i) = norm( Bzpyp );
      end
   end

   c2         = max(max(fc2));
   [Ic2, Jc2] = find(fc2==c2);
   seq        = fseq(Jc2);
   p          = Prhp(Ic2);

end
e ZST argument.')
         return;
      else
         disp('Warning: Two directions are possible pole directions.')      
         disp('         May try to decrease the ZST argument.')
      end
   end

   for j=1:np
      [Bzps fseq]  = allbzp(YZ,Zrhp,Prhp(j));
      for i=1:fact(nz)
         Bzpyp = Bzps(:,ny*delay.m000064400136200001130000000012500605462207500127770ustar00kjetilentek00000000000000
% function delay_g = delay(theta,omega,nu,epp)
%
% delay_g = delay(theta,omega,nu,epp)
% delay_g : frequency respones of  a delay
%           delay_g = exp(-s*theta)
% theta : time delay
% omega : frequency vector
% nu    : dimension of varying matrix     (default = 1)
% epp   : small negative A-matrix element (default = 0)

function delay_g = delay(theta,omega,nu,epp)

if (nargin == 0) | (nargin > 4),
   disp('usage: delay_g = delay(theta,omega,nu,epp)')
   return
end
if nargin < 4,
   epp = 0;
end
if nargin < 3,
   nu = 1;
end

if epp > 0,
   epp = -epp;
end

i = pck(epp,1,1,0);
tmp = frsp(i,omega,theta);
delay_g = tmp;
for j=2:nu,
   delay_g = daug(delay_g,tmp);
end

 theta : time delay
% omega : frequency vector
% nu    : dimension of varying matrix     (default = 1)
% epp   : small negative A-matrix element (default = 0)

function delay_g = delay(theta,omega,nu,epp)

if (nargin == 0) | (nargin > 4),
   disp('usage: delay_g = delay(theta,omega,nu,epp)')
   return
end
if nargin < 4,
   epp = 0;
end
if narfact.m000064400136200001130000000002170611361337700126210ustar00kjetilentek00000000000000function y = fact(x)
   if x == 1
      y = 1;
      return
   end
   
   if x <= 0
      y = 1;
      return
   end
   
   y=x*fact(x-1);
end
 daug(delay_g,tmp);
end

 theta : time delay
% omega : frequency vector
% nu    : dimension of varying matrix     (default = 1)
% epp   : small negative A-matrix element (default = 0)

function delay_g = delay(theta,omega,nu,epp)

if (nargin == 0) | (nargin > 4),
   disp('usage: delay_g = delay(theta,omega,nu,epp)')
   return
end
if nargin < 4,
   epp = 0;
end
if narinpcomp.m000064400136200001130000000017760602752112400133540ustar00kjetilentek00000000000000%
%  c = inpcomp( G, ip );
%
%  Returns the complement set of inputs described in ip. All input indices
%  not in ip. Example: Number of inputs is four and ip=[2 3]; then c=[1 4];
% 
%  Inputs:  G  - Matrix describing the system, can bd system, constant or 
%                varying.
%           ip - Set of inputs to be complemented.
%
%  Outputs: c  - complement set.
%
%  Written by: Kjetil Havre, e-mail: kjetil@ife.no
%  
%  Note: This is for sure not the most effective way to find the complemet 
%        set in matlab.
%

function c = inpcomp( G, ip )
   c = [];
   [mt, rd, cd, nd] = minfo(G);
   no = max(size(ip));
% Added 9/8-95 KH
   if strcmp(ip,':') == 1 
      return
   end
% to this spot 9/8-95.
   if no <= cd
      ipc = [];
      for i=1:cd,                     
	 found = 0;
	 for j=1:no,
	    if ip(j) == i;
	       found = 1;
	    end
	 end
	 if found == 0
	    ipc = [ipc i];
	 end
      end
      c = ipc;
   else
      disp('Number of outputs is larger than columns' );
      return;
   end
end
  inpcompfact.m000064400136200001130000000146240606376517500142260ustar00kjetilentek00000000000000%
%  Factorize a pair of complex RHP-zeros in input factorization.
%
%  [GI, BI, Xf, Uf, Zf] = inpcompfact( G, Z, epp )
%  
%  The INPCOMPFACT function is used to obtain an input factorization of
%  a pair of complex RHP-zeros so that G(s) = Gm(s)*BI(s).
%  The input argument Z contains the RHP-zeros. The optional EPP argument 
%  is used to decide the difference between the zero given in Z and the 
%  zeros calculated by the function. If the function does not calculate
%  any zero Ze so that abs(Z(i) - Ze) < EPP the function returns with 
%  BI = eye(nu) and GI = G. The default EPP value is 1e-12.
%  The IDZE function is used to obtain the zero state direction Xz and the
%  input zero direction Uz for RHP-zeros Z.
%  
%  BI(s) is an all pass transfer function for s=jw. All singular values
%  are equal to one, and the absolute value of all eigenvalues are 
%  equal to one. 
%

%  GI is constructed as follows:
%           G(s) = |  A  |  B  |      GI(s) = |  A  |   Bm  |
%                  |-----|-----|              |-----|-------|
%                  |  C  |  D  |              |  C  |   D   |
%  where Bm is obtained by:
%      Bm = B - 2*Z*Xz*Uz'.
%
%  Xz and Uz are the state and the input direction corresponding to the 
%  RHP-zero Z for the system G(s) = |  A  |  B  |  (see IZDE)
%                                   |-----|-----|
%                                   |  C  |  D  |
%
%  Possible malfunctions: Multiple zeros located at the same place in the 
%                         complex plane.
%
%  Written by: Kjetil Havre,  12/12 - 1995, e-mail: kjetil@ife.no                         
%
%  See also: IZDE, INPCOMPFACT, OUTREALFACT, OUTCOMPFACT, INPFACT, OUTFACT.
%
%
%
%
%

function [GI, BI, Xz, Uz, Zf] = inpcompfact( G, Z, epp )

   if nargin < 2
      disp('Usage: inpcompfact(G, Z)');
      return
   end
   if nargin < 3
      epp = 1e-12;
   end

   [mtype,ny,nu,nx] = minfo(G);
   if mtype ~= 'syst'
      disp('Error: System matrix G is required.');
      disp('usage: inpcompfact(G, Z)');
      GI = G; BI = eye( nu ); Xz = []; Uz = []; Zf = [];
      return
   end

   [nz, nc] = size(Z);
   if nz < 2
      disp('Error: Two zeros are required.');
      GI = G; BI = eye(nu); Xz = []; Uz = []; Zf = [];
      return     
   end
   if nz > 2
      disp('Warning: Only the first zero and its complexed conjugate is factorized.')
   end
   z1 = Z(1);
   if isreal(z1)
      disp(['Error: Zero is supposed to be complex Z: ', num2str(z), '.']);
      GI = G; BI = eye( nu ); Xz = []; Uz = []; Zf = [];
      return;
   end

   Izc=find(abs(Z-conj(z1))<epp);
   if isempty(Izc)
      disp( ['Error: Unable to find complex conjugate zero ', num2str(conj(z1)),' in Z.'] );
      GI = G; BI = eye(nu); Xz = []; Uz = []; Zf = [];
      return;      
   end
   
   [Ze, uiz, xiz] = izde( G );              % Compute directions.
   
   Iz1 = find( abs(Ze-z1) < epp );  
   if isempty(Iz1)
      disp( ['Error: Requested zero Z(1): ',...
	           num2str(z1), ' is not in system.' ] );      
      GI = G; BI = eye(nu); Xz = []; Uz = []; Zf = [];
      return;
   end

   Xz1 = xiz(:,Iz1);    Uz1 = uiz(:,Iz1);
   [A,B,C,D] = unpck( G );
   evk1 = [Xz1; Uz1];
   M = [A-eye(size(A))*z1   B;  C    D ];
   epvk1=M*evk1;
   nepvk1 = norm(epvk1);
   
   if nepvk1 > epp
      disp(['Error: Unable to solve generalized eigenvalue problem: ', num2str(nepvk1),'.'] )
      GI = G; BI = eye(nu); Xz = []; Uz = []; Zf = [];
      return;
   end
   
   
   Bn = B - 2*real(z1)*Xz1*Uz1';
   Gn = pck(A,Bn,C,D);

   [Ze, uiz, xiz] = izde( Gn );             % Compute directions.
   Iz2 = find( abs(Ze-conj(z1)) < epp );

   if isempty(Iz2)
      disp( ['Error: Requested zero conj(Z(1)): ',...
	           num2str(conj(z1)), ' is not in system.' ] );      
      GI = G; BI = eye(nu); Xz = []; Uz = []; Zf = [];
      return;
   end
   Xz2 = xiz(:,Iz2);
   Uz2 = uiz(:,Iz2);
   evk2 = [Xz2; Uz2];
   M = [A-eye(size(A))*conj(z1)   Bn;  C    D ];
   epvk2=M*evk2;
   nepvk2 = norm(epvk2);
   
   if nepvk2 > epp
      disp(['Error: Unable to solve generalized eigenvalue problem: ', num2str(nepvk2),'.'] )
      GI = G; BI = eye(nu); Xz = []; Uz = []; Zf = [];
      return;
   end
   
   Xz2t = conj(Xz1);
   Uz2t = conj(Uz1);
   
   EUz = Uz2 - conj(Uz1)
   

   Xz     = [Xz1 Xz2];
   Uz     = [Uz1 Uz2];
   Bn     = B - 2*real(z1)*Xz*Uz'
   Bnr    = real(Bn);
   eppBnr = norm(Bn - Bnr);
   
   if eppBnr > epp
      disp( 'Warning: A real realization of GI is not found.' );
      Bnr = Bn;       
   end
   
   GI = pck(A,Bn,C,D);
   Ab  = -conj(z1); Bb = 2*real(z1)*Uz1'; Cb  = -Uz1; Db = eye(nu);
   B1  = pck(Ab,Bb,Cb,Db);
   Ab  = -z1; Bb = 2*real(z1)*Uz2'; Cb  = -Uz2; Db = eye(nu);   
   B2  = pck(Ab,Bb,Cb,Db);
   
   Ac = [-real(z1) imag(z1); -imag(z1) -real(z1)];
   Cc  = [eye(min(ny,2)); zeros(ny-min(ny,2),2)];
   Xo = -pinv(Cc)*Uz;
   
   BI  = mmult(B2, B1);
   [Ab,Bb,Cb,Db] = unpck(BI);
   
   Ab
   
   Af = [-real(z1) imag(z1); -imag(z1) -real(z1)]
   
   T = are(zeros(2,2), Ab, Af)
   
   inv(T)*Af*T-Ab
   
   a1 = Ab(1,1); a2 = Ab(1,2); a3 = Ab(2,1); a4 = Ab(2,2);
   im = imag(z1);  re = real(z1);
   
   M = [ a1+re   0     a2   im
         -im   a1+re    0   a2
	  a3    0    a4+re  im
	  0     a3    -im a4+re ];
   rkm = rank(M)
   
   if rkm == 4
      disp('Warning: Can not find similarity transformation to convet to real realization.')
      T = eye(2);
   else
      [U,S,V] = svd(M);
      X = V(:,4)
      T(1,1) = X(1); T(1,2) = X(2); T(2,1) = X(3); T(2,2) = X(4);
      T
   end
   
   
   
   Af
   An = real(inv(T)*Ab*T)
   eig(An)
   Bn = inv(T)*Bb
   % angle(Bn)*180/pi
   Cn =  Cb*T
   
   % angle(Cn)*180/pi
   
   Dn = Db
   
   BI = pck(An,Bn,Cn,Dn);

%   BI0 = Cn*inv(-An)*Bn + Dn   
%   inv(BI0)  
%   Cn*Bn
%
%   Since BI0 is real and inv(B10) is real and Cn*Bn is real
%   It shuld be possible to find a unitary matrix UB such that
%   B in BI is real.
%
%
   
%   f1 = nd2sys( [1 -conj(z1)], [1 2*real(z1) z1*conj(z1)], 4*real(z1) )
%   f2 = nd2sys( [1 -z1], [1 2*real(z1) z1*conj(z1)], 4*real(z1) )
%   b1 = mmult(Uz1,f1,Uz1');
%   b2 = mmult(Uz2,f2,Uz2');
%   BI = msub(eye(nu), b1, b2);
      
%    f1 = nd2sys( [1 -2*real(z1) z1*conj(z1)], [1 2*real(z1) z1*conj(z1)], 1 );
%    Df = daug(f1, f1);
%    UzUzH  = Uz*Uz';
%    UzUz = real(UzUzH);
%    if norm(UzUzH-UzUz) > epp
%       disp('Warning: Real realization of Uz*Uz^H is not acheived.')
%       UzUz = UzUzH;
%    end
%    BI = mmult(Df, UzUz);
%
%    This BI does not work. 
%    Kjetil Havre 14/12-1995.
%

end
*real(z1) z1*conj(z1)], 4*real(z1) )
%   b1 = mmult(Uz1,f1,Uz1');
%   b2 = mmult(Uz2,f2,Uz2');
%   BI = msubinpcprfact.m000064400136200001130000000014230606403357500140360ustar00kjetilentek00000000000000%
%   Function INPCPRFACT finds factorization such that G(s) = W(s)^{-1}*V(s),
%   where V(jw)^H*V(jw) = I, and both W and V are stable. That means
%   all RHP-zeros in V and all RHP-poles as RHP-zeros in W. 
%   System G must be square and have no jw-axis poles. 
%   This is performed by factorizing G^T(s) = N(s)*M^{-1}(s) with
%   INPCRPFACT. The filters are then given by N(s) = V^T(s) and
%   Q = M^T(s).
%
%   The outputs: Q, R, E, X and Xerr are outputs returned from OUTCRPFACT.
%
%   Reference: Morari M. ``Robust Process Control'', Prentice Hall 1989.
%

function [V, W, Wi, Q, R, E, X, Xerr] = inpcprfact( G, deps, epp )

   GT = transp(G);
   
   [N, M, Mi, Q, R, E, X, Xerr] = outcprfact( GT, deps, epp );
   
   V  = transp(N);
   W  = transp(M);
   Wi = transp(Mi);

end
s performed by factorizing G^T(s) = N(s)*M^{-1}(s) with
%   INPCRPFACT. The filters are then given by N(s) = V^T(s) and
%   Q = M^T(s).
%
%   The outputs: Q, R, E, X and Xerr are outputs returned from OUTCRPFACT.
%
%   Reference: Morari inprealfact.m000064400136200001130000000065510606330065700142010ustar00kjetilentek00000000000000%
%  Factorize one single real RHP-zero in input factorization.
%
%  [GI, BI, Xf, Uf, Zf] = inprealfact( G, Z, epp )
%  
%  The INPREALFACT function is used to obtain an input factorization of
%  a single RHP-zero so that G(s) = Gm(s)*BI(s).
%  The input argument Z contains the RHP-zero. The optional EPP argument 
%  is used to decide the difference between the zero given in Z and the 
%  zeros calculated by the function. If the function does not calculate
%  any zero Ze so that abs(Z(i) - Ze) < EPP the function returns with 
%  BI = eye(nu) and GI = G. The default EPP value is 1e-12.
%  The IDZE function is used to obtain the zero state direction Xz and the
%  input zero direction Uz for RHP-zero Z.
%  
%  BI(s) is an all pass transfer function for s=jw. All singular values
%  are equal to one, and the absolute value of all eigenvalues are 
%  equal to one. 
%
%  GI is constructed as follows:
%           G(s) = |  A  |  B  |      GI(s) = |  A  |   Bm  |
%                  |-----|-----|              |-----|-------|
%                  |  C  |  D  |              |  C  |   D   |
%  where Bm is obtained by:
%      Bm = B - 2*Z*Xz*Uz'.
%
%  Xz and Uz are the state and the input direction corresponding to the 
%  RHP-zero Z for the system G(s) = |  A  |  B  |  (see IZDE)
%                                   |-----|-----|
%                                   |  C  |  D  |
%
%  Possible malfunctions: Multiple zeros located at the same place in the 
%                         complex plane.
%
%  Written by: Kjetil Havre,  12/12 - 1995, e-mail: kjetil@ife.no                         
%
%  See also: IZDE, INPCOMPFACT, OUTREALFACT, OUTCOMPFACT, INPFACT, OUTFACT.
%

function [GI, BI, Xz, Uz, Zf] = inprealfact( G, Z, epp )

   if nargin < 2
      disp('Usage: inprealfact(G, Z)');
      return
   end
   if nargin < 3
      epp = 1e-12;
   end
   
   [mtype,ny,nu,nx] = minfo(G);
   if mtype ~= 'syst'
      disp('Error: System matrix G is required.');
      disp('usage: inprealfact(G, Z)');
      GI = G; BI = eye( nu ); Xz = []; Uz = []; Zf = [];
      return
   end

   [nz, nc] = size(Z);
   z = Z(1);
   if nz < 1
      disp('Message: No zeros?');
      GI = G; BI = eye(nu); Xz = []; Uz = []; Zf = [];
      return     
   end
   if nz > 1
      disp('Warning: Only the first zero is factorized.')
   end
   
   if ~isreal(z)
      disp(['Error: Zero is supposed to be real Z: ', num2str(z), '.']);
      GI = G; BI = eye( nu ); Xz = []; Uz = []; Zf = [];
      return;
   end

   [Ze, uiz, xiz] = izde( G );   % Compute directions.
   Iz = find( abs(Ze-z) < epp );  
   if isempty(Iz)
      disp( ['Warning: Requested zero Z(', int2str(lv),'): ',...
	           num2str(Z(1)), ' is not in system.' ] );      
      GI = G; BI = eye(nu); Xz = []; Uz = []; Zf = [];
      return;
   end
   z = real(Ze(Iz));
   Xz = xiz(:,Iz); Uz = uiz(:,Iz);
   inz = find( abs(Uz)>1000*epp );
   vh = angle(Uz(inz(1)))
   Uz
   Xz
   Uz = exp(-vh*sqrt(-1))*Uz
   Xz = exp(-vh*sqrt(-1))*Xz
   err = norm([Xz; Uz] - real([Xz; Uz]));
   if err > epp
      disp('Error: State and input direction is not real.');
      Xz
      Uz
      GI = G; BI = eye( nu ); Xz = []; Uz = []; Zf = [];
      return;      
   end
   Xz = real(Xz); Uz = real(Uz);
   
   [A,B,C,D] = unpck( G );
   Bn  = B - 2*z*Xz*vtp(Uz);
   GI  = pck(A,Bn,C,D);
   Ab  = -z; Bb = 2*z*vtp(Uz); Cb  = -Uz; Db = eye(nu);
   BI  = pck(Ab,Bb,Cb,Db);
   Zf  = z; 
end
  

   Uz
   Xz
   Uz = exp(-vh*sqrt(-1))*Uz
   Xz = exp(-vh*sqrt(-1))*Xz
   err = norm([Xz; Uz] - real([Xz; Uz]));
   if err > epp
      disp('Error: Staipf.m000064400136200001130000000066450610345370200124660ustar00kjetilentek00000000000000%  function [GI, BI, Up, Xp, Pf] = opf( G, P, epp )
%  
%  The IPF function is used to obtain an output factorization of RHP-poles 
%  so that G(s) = GI(s)*BI(s) where BI contains the RHP-poles.
%  The input argument P is a vector containing the RHP-poles. 
%  The PDFST function is used to obtain the pole state direction Xpi and 
%  the output zero direction Upi for the RHP-pole P(i). The optional EPP 
%  argument is sendt to the PDFST function as the SST argument. The default 
%  EPP value is 1e-12.
%  
%  BI(s) is an all pass transfer function for s=jw. All singular values
%  are equal to one, and the absolute value of all eigenvalues are 
%  equal to one. 
%
%  GI is constructed as follows:
%           G(s) = |  A  |  B  |      GI(s) = |  Am  |  B  |
%                  |  C  |  D  |              |  Cm  |  D  |
%  where Am and Cm are obtained by running the following formulas repeatedly:
%      Ai = Ai_1 - 2*Re(Pi)*B*Upi*Xpi',  Ao = A for i=1 to Np, where Np
%                                        is the number of RHP-poles.
%                                        Am = Ai, with i = Np.
%      Ci = Ci_1 - 2*Re(Pi)*D*Upi*Xpi'   Co = C for i=1 to Np, where Np
%                                        is the number of RHP-poles.
%                                        Cm = Ci, with i = Np.
%
%  BI(s) is given by:
%      BI(s) = PROD Bi(s), where Bi(s) = I + 2*Re(Pi)/(s-Pi') *Upi*Upi'
%                          the index runs from i=1 to Np.
%  Xpi and Upi are the state and the output direction corresponding to the 
%  RHP-pole Pi for the system Gi_1(s) = |  Ai_1  |  B  |
%                                       |--------|-----|
%                                       |  Ci_1  |  D  |
%
%
%  Written  by: Kjetil Havre,  30/1 - 1996, e-mail: kjetil@ife.no
%
function [GI, BI, Up, Xp, Pf] = opf( G, P, epp, lv )
   if nargin < 2
      disp('usage: opf(G, Z)')
      GO = G; BO = eye( nu );
      Yp = []; Xp = []; Pf = [];
      return
   end
   if nargin < 3
      epp = 1e-12;
   end
   if nargin < 4
      lv = 1;
   end
  
   [mtype,ny,nu,nx] = minfo(G);
   if mtype ~= 'syst'
      disp('Error: System matrix G is required.');
      disp('usage: opf(G, Z)');
      GI = G; BI = eye( nu );
      Up = []; Xp = []; Pf = [];
      return
   end

   [np, nc] = size(P);
   if nc > 1                % Vector  is horizontal.
      P = vtp(P);
      [np, nc] = size(P);
   end

   if np < 1
      disp('Message: No poles?');
      GI = G; BI = eye(nu);
      Up = []; Xp = []; Pf = [];
      return     
   end
   
   p = P(1);
   [up, yp, xip, xop, err] = pdfst( G, p, epp );   % Compute directions.
   if err != 0 
      disp( ['Warning: Requested pole P(', int2str(lv),'): ',...
	           num2str(P(1)), ' is not in system.' ] );      
      disp(['Error when computing pole direction, ',...
            'error number: ', int2str(err), '.']);
      GI = G; BI = eye(nu); Up = []; Xp = []; Pf = [];
      return;
   else
      [A,B,C,D] = unpck( G );
      An  = A - 2*real( p )*B*up*xip';
      Cn  = C - 2*real( p )*D*up*xip';
      GI  = pck(An,B,Cn,D);
      Ab  = p;
      Bb  = 2*real(p)*up';
      Cb  = up;
      Db  = eye(nu);
      BI  = pck(Ab,Bb,Cb,Db);
      Up  = up;
      Xp  = xip;
      Pf  = p;
   end
   
   if np > 1   % Call recursive for factorizing additional RHP-zeros.
      [Gn, Bn, Un, Xn, Pn] = ipf(GI, P(2:np), epp, lv+1);
      GI = Gn;
      BI = mmult(Bn,BI);
      Up = [Up Un];
      Xp = [Xp Xn];
      Pf = [Pf; Pn];
   end
end

   An  = A - 2*real( p )*B*up*xip';
      Cn  = C - 2*real( p )*D*up*xip';
      GI  = pck(isstable.m000064400136200001130000000006060606005344100135030ustar00kjetilentek00000000000000%   stable = isstable(G)
%   
%   Function istable returns one if all poles of G have real part 
%   less than zero.
%
%   See also: spoles and eig.
%
%   Written by: Kjetil Havre, 1/12-1995, e-mail: kjetil@ife.no.
%

function stable = isstable(G)

   P = spoles(G);
    
   IPrhp = find(real(P)>=0);
    
   if isempty(IPrhp)
      stable = 1;
   else
      stable = 0;
   end
end
       
    Pf = [Pf; Pn];
   end
end

   An  = A - 2*real( p )*B*up*xip';
      Cn  = C - 2*real( p )*D*up*xip';
      GI  = pck(izde.m000064400136200001130000000135140611041716500126350ustar00kjetilentek00000000000000%   [Z, U, X] = izde(sys,epp)
%
%   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 u 
%   are normalized so that U'*U = I. 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.
%   
%   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. 
%   Comments, corrections and malfunctions, can be e-mailed to:
%                               kjetil@ife.no
%
%   See also: EIG, SZEROS, OZDE and SPOLES. 
%   Algorithm based on Laub & Moore 1978 paper, Automatica
%

function [Z, U, X] = izde(sys,epp)
  if nargin < 1
    disp('usage: izde(sys)')
    return
  end
  if nargin == 1
    epp = eps;
  end
  [mtype,ny,nu,nx] = minfo(sys);
  if mtype == 'syst'
    [a,b,c,d] = unpck(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));
	  Inz = find( abs(vu(:,i)) > 1000*epp );
          U(:,i) = vu(:,i)/nrmu * exp( -angle(vu(Inz(1),i))*sqrt(-1) );
	  X(:,i) = vx(:,i)/nrmu * exp( -angle(vu(Inz(1),i))*sqrt(-1) );
      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)); 
	  Inz = find( abs(vu(:,i)) > 1000*epp );
          U(:,i) = vu(:,i)/nrmu * exp( -angle(vu(Inz(1),i))*sqrt(-1) );
	  X(:,i) = vx(:,i)/nrmu * exp( -angle(vu(Inz(1),i))*sqrt(-1) );
      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
%
 );
      vu = vech3(nx+1:nx+nu,:);
% Normalize columns.
      [nvr, nvc] = size( vu ); 
      for i=1:nvc,
          nrmu = norm(vu(:,i)); 
	  Inz = find( abs(vu(:,i)) > 1000*epp izf.m000064400136200001130000000077600610467035500125050ustar00kjetilentek00000000000000%  function [GI, BI, Uz, Xz, Zf] = izf( G, Z, epp )
%  
%  The IZF function is used to obtain an input factorization of
%  RHP - zeros so that G(s) = GI(s)*BI(s).
%  The input argument Z is a vector containing the RHP-zeros. The optional
%  EPP argument is used to decide the difference between the zeros given
%  in Z and the zeros calculated by the IZF function. If the function does
%  not calculate any zeros Ze so that abs(Z(i) - Ze) < EPP the function 
%  skipps Z(i) in the factorization. The default EPP value is 1e-12.
%  The IDZE function is used to obtain the zero state direction Xzi and the
%  input zero direction Uzi for RHP-zero Z(i).
%  
%  BI(s) is an all pass transfer function for s=jw. All singular values
%  are equal to one, and the absolute value of all eigenvalues are 
%  equal to one. 
%
%  GI is constructed as follows:
%           G(s) = |  A  |  B  |      GI(s) = |  A  |   Bm  |
%                  |-----|-----|              |-----|-------|
%                  |  C  |  D  |              |  C  |   D   |
%  where Bm is obtained by running the following formula repeatedly:
%      Bi = Bi_1 - 2*Re(Zi)*Xzi*Uzi',   Bo = B for i=1 to Nz, where Nz
%                                       is the number of RHP - zeros.
%                                       Bm = Bi, with i = Nz.
%
%  BI(s) is given by:
%      BI(s) = PROD Bi(s), where Bi(s) = I - 2*re(Zi)/(s+Zi') *Uzi*Uzi'
%                          the index runs from i=Nz to one.
%  Xzi and Uzi are the state and the input  direction corresponding to the 
%  RHP-zero Zi for the system Gi_1(s) = |  A  |  Bi_1  |
%                                       |-----|--------|
%                                       |  C  |    D   |
%
%  Note In order to obtain an real input factorization for complex RHP-zeros
%  the two etc.
%
%  Possible malfunctions: Multiple zeros located at the same place in the 
%                         complex plane.
%                         RHP-zeros on the jw axis, pure imaginary zeros.
%                         
%
%
%  Written  by: Kjetil Havre,  15/5 - 1995, e-mail: kjetil@ife.no 
%  Modified by: Kjetil Havre,  30/1 - 1998, e-mail: kjetil@ife.no
%     Modification was to also return the state direction, and change of name.
%

function [GI, BI, Uz, Xz, Zf] = izf( G, Z, epp, lv )
   if nargin < 2
      disp('usage: izf(G, Z)');
      return
   end
   if nargin < 3
      epp = 1e-12;
   end
   if nargin < 4
      lv = 1;
   end
   
   [mtype,ny,nu,nx] = minfo(G);
   if mtype ~= 'syst'
      disp('Error: System matrix G is required.');
      disp('usage: izf(G, Z)');
      GI = G;
      BI = eye( nu );
      Uz = []; Xz = []; Zf = [];
      return
   end

   [nz, nc] = size(Z);
   if nc > 1                % Vector  is horizontal.
      Z = vtp(Z);
      [nz, nc] = size(Z);
   end
  
   if nz < 1
      disp('Message: No zeros?');
      GI = G;
      BI = eye(nu);
      Uz = []; Xz = []; Zf = [];
      return     
   end
  
   z = Z(1);
   [Ze, uiz, xiz] = izde( G );   % Compute directions.

   Iz = find( abs(Ze-z) < epp );
   if isempty(Iz)
      disp( ['Warning: Requested zero Z(', int2str(lv),'): ',...
	           num2str(Z(1)), ' is not in system.' ] );      
      GI = G;
      BI = eye(nu);
      Uz = []; Xz = []; Zf = [];
   else
                     % Iz contains the index to the zero in Ze. 
                     % Factorize this zero.
      if isreal(z)
         xiz(:,Iz) = real(xiz(:,Iz));
         uiz(:,Iz) = real(uiz(:,Iz));
	 Ze(Iz) = real(Ze(Iz));
      end
      [A,B,C,D] = unpck( G );
      Bn  = B - 2*real( Ze(Iz) )*xiz(:,Iz)*uiz(:,Iz)';
      GI  = pck(A,Bn,C,D);
      Ab  = -conj(Ze(Iz));
      Bb  = 2*real(Ze(Iz))*uiz(:,Iz)';
      Cb  = -uiz(:,Iz);
      Db  = eye(nu);
      BI  = pck(Ab,Bb,Cb,Db);
      Uz  = uiz(:,Iz);
      Xz  = xiz(:,Iz);
      Zf  = Ze(Iz); 
   end
  

   if nz > 1   % Call recursive for factorizing additional RHP-zeros.
      [Gn, Bn, Un, Xn, Zn] = izf(GI, Z(2:nz), epp, lv+1);
      GI = Gn;
      BI = mmult(Bn, BI);
      Uz = [Uz Un];
      Xz = [Xz Xn];
      Zf = [Zf; Zn];
   end

end

*xiz(:,Iz)*uiz(:mincond.m000064400136200001130000000007740602506227000133330ustar00kjetilentek00000000000000function gmin = condmin(Gf)
   [mt,ny,nu,nx] = minfo(Gf); gmin = []; 

%   [dum1,n,m,dum2]=minfo(Gf);
   zeroy = zeros(ny,ny);
   zerou = zeros(nu,nu);
   H = abv( sbs(zerou, vpinv(Gf)), sbs(Gf, zeroy));
%
%  H = |  [0]_{nuxnu}    G^{-1}    |
%      |      G       [O]_{nyxny}  |
%
%  
   blk0=[1 1]; blk1=blk0;
   for k=2:nu+ny
     blk1 = [blk1; blk0];   % All of size 1x1.
   end
%   blk2=blk1; blk=[blk1; blk2];
   blk = blk1;
   mu1=mu(H,blk,'Uc'); gmin1=sel(mu1,':',1);
   gmin=mmult(gmin1,gmin1);
end
iz(:minicomd.m000064400136200001130000000007360603624256400135110ustar00kjetilentek00000000000000function gmin = minicond(Gf)
   [mt,ny,nu,nx] = minfo(Gf); gmin = []; 

   zeroy = zeros(ny,ny);
   zerou = zeros(nu,nu);
   H = abv( sbs(zerou, vpinv(Gf)), sbs(Gf, zeroy));
%
%  H = |  [0]_{nuxnu}    G^{-1}    |
%      |      G       [O]_{nyxny}  |
%
%  
   blk0=[1 1]; blk1=blk0;
   for k=2:nu+ny
     blk1 = [blk1; blk0];   % All of size 1x1.
   end
%   blk2=blk1; blk=[blk1; blk2];
   blk = blk1;
   mu1=mu(H,blk,'Uc'); gmin1=sel(mu1,':',1);
   gmin=mmult(gmin1,gmin1);
end
 gmin=mmult(gmin1,gmin1);
end
iz(:mkspr.m000064400136200001130000000016630603671200500130360ustar00kjetilentek00000000000000function Gm = mkspr(G,z);
%
% Virker darlig
%
  [mt,ny,nu,nx] = minfo(G);
  if strcmp(mt, 'syst') == 0
     disp('Usage: Gm = mkspr(G,z); where G is system');
     return
  end  
  
  n = nu;
  
  g1  = sel(G,1,1:nu);
  
  Gs  = sel(G,2:ny,1:nu);
  for i=1:n,
     g1i = sel(g1,1,i);
     [ai,bi,ci,di] = unpck(g1i);
%     [ai,bi,ci,ci] = minreal(ai,bi,ci,di);
     [zi,pi,ki] = ss2zp(ai,bi,ci,di);
     nzi = max(size(zi)); npi = max(size(pi));
     if nzi > npi
        for l=1:npi-nzi,
	   zi(nzi+l) = z;
	end
     end
     gn = zp2sys(zi,pi,ki/z);
     g1 = sbs(g1,gn);
%     [a,b,c,d] = unpck(g1); [a,b,c,d] = minreal(a,b,c,d); g1 = pck(a,b,c,d);
     g1 = sysbal(g1);
  end
  [mt,nrg1,ncg1,nxg1] = minfo(g1)
  g1n = sel(g1, 1,nu+1:nu+nu);
%  [a,b,c,d] = unpck(g1n); [a,b,c,d] = minreal(a,b,c,d); g1n = pck(a,b,c,d);
  Gm  = abv(g1n, Gs);
  Gm = sysbal(Gm);
%  [a,b,c,d] = unpck(Gm); [a,b,c,d] = minreal(a,b,c,d); Gm = pck(a,b,c,d);
end
    
size(pi));
     if nzi > npi
        for l=1:npi-nzi,
	   zi(nzi+l) = z;
	endmusynflpS.m000064400136200001130000000114070605462722500137100ustar00kjetilentek00000000000000% function [dsysl,dsysr] = musynfit(pre_dsysl,Ddata,sens,blk,nmeas,...
%                                    ncntrl,clpg,upbd,wt)
%
%   Fits the magnitude curve obtained by multiplying the old D 
%   frequency response (from PRE_DSYSL) with the DDATA data. 
%   Returns stable, minimum phase system matrices DSYSL and 
%   DSYSR, which can be absorbed into the original interconnection 
%   structure with multiplications (using MMULT), and one inverse 
%   (using MINV). Once absorbed, a H_infinity design can be 
%   performed, completing another iteration of MU-SYNTHESIS.
%
%   For the first MU-SYNTHESIS iteration, the variable
%   PRE_DSYSL should be set to the string  'first'
%   In subsequent iterations, PRE_DSYSL should be the previous
%   (left) rational D-scaling system matrix, DSYSL.
%
%   The (optional) variable CLPG is the matrix that produced the DDATA, SENS
%   and UPBD data. Upon fitting the DDATA, the stable, minimum phase
%   system matrices DSYSL and DSYSR, are absorbed into the original
%   matrix CLPG and plotted along with the UPBD data. A comparison of
%   these two plots provides insight into how well the rational system
%   matrices DSYSL and DSYSR fit the DDATA. The newly scaled matrix DMDI, 
%   which is CLPG with DSYSL and DSYSR wrapped is then return. If CLPG 
%   and UPBD are not provided, the default is to plot the SENS 
%   variable instead. 
%
%   The optional WT variable allows you to weight the DDATA, the default 
%   is no additional weighting on the DDATA. 
%
%   Identical to MUSYNFIT except that FITMAGLP is called instead of
%   FITMAG, and this calls MAGFIT. 
%
%   See also: DKIT, FITMAGLP, FITSYS, MAGFIT, and MUFTBTCH and MUSYNFIT.

%[dsysl,dsysr]=musynflp(dsysl,rowd,sens,blk,nmeas,nu);
%osysl=dsysl;magdata=rowd;ncntrl=nu;nargin=6;

function [dsysl,dsysr,newdata] = musynflpS(osysl,magdata,sens,blk,nmeas,ncntrl,...
                                    clpg,upbd,wt)
% CLPG has the Previous rational D's
% MAGDATA is the new incremental D's
% The returned DMDI is  DSYSL*CLPG*DSYSR^{-1}

if nargin < 6 | nargin == 7
 disp('usage: [dsysl,dsysr] = musynflp(pre_dsysl,Ddata,sens,blk,nmeas,ncntrl)')
 return
end

[nblk,x] = size(blk);
leftdim = nmeas;
 for i=1:nblk
  if blk(i,2) == 0
    leftdim = leftdim + blk(i,1);
    if blk(i,1) > 1
      disp('full Ds not implemented yet')
      return
    else
      blk(i,2) = 1;
    end
  else
    leftdim = leftdim + blk(i,2);
  end
 end

 [mtype,mrows,mcols,mnum] = minfo(magdata);
 omega = magdata(1:mnum,mcols+1);
 [stype,srows,scols,snum] = minfo(sens);
 if mnum ~= snum
   error('SENS & MAG data have diff # of points - rerun MU')
   return
 end
 if nargin == 9
   [wtype,wrows,wcols,wnum] = minfo(wt);
   if wnum ~= mnum
    error('NOMG & MAG data have diff # of points - rerun MU ')
    return
   end
   supwt = wt(1:mnum,1);
 else
   supwt = ones(mnum,1);
 end
 if nargin == 8
   if ~indvcmp(clpg,magdata)
    error('Closed-loop freq response (CLPG) & MAGDATA have different IVs')
    return
   end
 end

 if nargin > 6
   [dleft,dright] = unwrapd(magdata,blk);
   dmdi = mmult(dleft,clpg,minv(dright)); % rational D's, VABS(data)
 end
%  NORMALIZE THE DATA BY THE LAST DSCALE: This rapidly scales
%  all of the D-scalings by the inverse of the last D scaling.
  
   [ddd,ddrp,ddindv,dder] = vunpck(magdata);
   [ddrows,ddcols] = size(ddd);
   dddn = ddd(:,ddcols*ones(1,ddcols)).\ddd;
   magdata(1:ddrows,1:ddcols) = dddn;

 if isstr(osysl)
   if strcmp(osysl,'first');
     osysl_g = eye(leftdim);
   else
     error('First argument should be the string   ''first''   or a SYSTEM matrix')
     return  
   end       
 else        
   [dum1,dum2,dum3,dum4] = minfo(osysl);
   if strcmp(dum1,'syst')
     osysl_g = frsp(osysl,omega);
   elseif strcmp(dum1,'cons')
     osysl_g = osysl;
   else      
     error('First argument should be the string   ''first''   or a SYSTEM matrix')
     return  
   end       
 end         

 loc = 1;
 dsysl = [];
 dsysr = [];
 for i=1:nblk-1
   wt = vpck((abs(sens(1:mnum,i) .* supwt)),omega);
   osysl_g_loc = sel(osysl_g,loc,loc);
   newdata = vabs(mmult(osysl_g_loc,sel(magdata,1,i)));
   heading = ['FITTING D SCALING #' int2str(i) ' of ' int2str(nblk-1)];
   if nargin == 6
     [scalard] = fitmaglp(newdata,wt,heading,osysl_g_loc);
   else
     [scalard] = fitmaglp(newdata,wt,heading,osysl_g_loc,dmdi,upbd,blk,i);
   end
   tmpsys = [];
   for j=1:min([blk(i,1) blk(i,2)])
     tmpsys = daug(tmpsys,scalard);
   end
   dsysl = daug(dsysl,tmpsys);
   dsysr = daug(dsysr,tmpsys);
   if blk(i,1) < blk(i,2)
     for j=1:(blk(i,2)-blk(i,1))
       dsysl = daug(dsysl,scalard);
     end
   elseif blk(i,2) < blk(i,1)
     for j=1:(blk(i,1)-blk(i,2))
       dsysr = daug(dsysr,scalard);
     end
   end
   loc = loc + blk(i,2);
 end
 dsysl = daug(dsysl,eye(blk(nblk,2)+nmeas));
 dsysr = daug(dsysr,eye(blk(nblk,1)+ncntrl));
blk,i);
   end
   tmpsys = [];
   for j=1:min([blk(i,1) blk(i,2)])
     tmpsys = daug(tmpsys,scalard);
   end
   dsysl = daug(dsysl,tmpsys);
   dsysr = daug(dsysr,tmpsys);
   if blk(i,1) < blk(i,2)
     for j=1:(blk(i,2)-blk(i,1))
       dsysl = dauncomb.m000064400136200001130000000000710611361340400127670ustar00kjetilentek00000000000000function y=ncomb(a,b)

   y = fact(a)/fact(b)/fact(a-b);
     for j=1:(blk(i,1)-blk(i,2))
       dsysr = daug(dsysr,scalard);
     end
   end
   loc = loc + blk(i,2);
 end
 dsysl = daug(dsysl,eye(blk(nblk,2)+nmeas));
 dsysr = daug(dsysr,eye(blk(nblk,1)+ncntrl));
blk,i);
   end
   tmpsys = [];
   for j=1:min([blk(i,1) blk(i,2)])
     tmpsys = daug(tmpsys,scalard);
   end
   dsysl = daug(dsysl,tmpsys);
   dsysr = daug(dsysr,tmpsys);
   if blk(i,1) < blk(i,2)
     for j=1:(blk(i,2)-blk(i,1))
       dsysl = daunewaxis.m000064400136200001130000000005320604272616600133640ustar00kjetilentek00000000000000function newaxis(lim, ori)

axis('off');
hold on
xh=plot([lim(1) lim(2)],[ori(2) ori(2)]); set(xh, 'Color', [1 1 1]);
yh=plot([ori(1) ori(1)],[lim(3) lim(4)]); set(yh, 'Color', [1 1 1]);
arrow(lim(2),ori(2),0,0.05*(lim(2)-lim(1)),0.03*(lim(4)-lim(3)),[1 1 1]);
arrow(ori(1),lim(4),90,0.05*(lim(4)-lim(3)),0.03*(lim(2)-lim(1)),[1 1 1]);
hold off

ug(tmpsys,scalard);
   end
   dsysl = daug(dsysl,tmpsys);
   dsysr = daug(dsysr,tmpsys);
   if blk(i,1) < blk(i,2)
     for j=1:(blk(i,2)-blk(i,1))
       dsysl = dauonealpha.m000064400136200001130000000027560605411712100134730ustar00kjetilentek00000000000000function [out,delta] = onealpha(sysg,blk,alphav,pertsize,perctol)

 [nblk,dum] = size(blk);
 [mtype,mrows,mcols,mnum] = minfo(sysg);
 cperf = blk(nblk,1);
 rperf = blk(nblk,2);
 m11row = mrows - rperf;
 m11col = mcols - cperf;
 if blk(nblk,2) == 0
   rperf = blk(nblk,1);
 end
 alphav = abs(alphav);

 xl = [];
 xu = [];
 yl = [];
 yu = [];
 delta = [];
 bestdelta = [];
 dist = inf;
 for i=1:length(alphav)
    alpha = alphav(i);
    lfac = diag([ones(1,mrows-rperf) sqrt(alpha)*ones(1,rperf)]);
    rfac = diag([ones(1,mcols-cperf) sqrt(alpha)*ones(1,cperf)]);
    malpha = mmult(lfac,sysg,rfac);
    [bnds,dvec,sens,pvec] = mu(malpha,blk,'csw');
    uppbnd = pkvnorm(sel(bnds,1,1));
    lowbnd = pkvnorm(sel(bnds,1,2));
    if abs(1/lowbnd-pertsize)/pertsize < perctol
        if strcmp(mtype,'vary')
            delta = dypert(pvec,blk,bnds,[1:nblk-1]);
            ndelta = hinfnorm(delta);
            if abs(ndelta(1,1)-pertsize) < dist
                bestdelta = delta; 
                dist = abs(ndelta(1,1)-pertsize);
            end
        elseif strcmp(mtype,'cons')
            delta = unwrapp(pvec,blk);
            delta = delta(1:m11col,1:m11row);
            ndelta = norm(delta);
            if abs(ndelta-pertsize) < dist
                bestdelta = delta; 
                dist = abs(ndelta-pertsize);
            end        
        end
    end
    xl = [xl;1/lowbnd];
    yl = [yl;lowbnd/alpha];
    xu = [xu;1/uppbnd];
    yu = [yu;uppbnd/alpha];
 end
delta = bestdelta;
out = [xl yl xu yu];
(ndelta(1,1)-pertsopf.m000064400136200001130000000066610610342443300124710ustar00kjetilentek00000000000000%  function [GO, BO, Yp, Xp, Pf] = opf( G, P, epp )
%  
%  The OPF function is used to obtain an output factorization of RHP-poles 
%  so that G(s) = BO(s)*GO(s) where BO contains the RHP-poles.
%  The input argument P is a vector containing the RHP-poles. 
%  The PDFST function is used to obtain the pole state direction Xpi and 
%  the output zero direction Ypi for the RHP-pole P(i). The optional EPP 
%  argument is sendt to the PDFST function as the SST argument. The default 
%  EPP value is 1e-12.
%  
%  BO(s) is an all pass transfer function for s=jw. All singular values
%  are equal to one, and the absolute value of all eigenvalues are 
%  equal to one. 
%
%  GO is constructed as follows:
%           G(s) = |  A  |  B  |      GO(s) = |  Am  |  Bm  |
%                  |  C  |  D  |              |  C   |  D   |
%  where Am and Bm are obtained by running the following formulas repeatedly:
%      Ai = Ai_1 - 2*Re(Pi)*Xpi*Ypi'*C,  Ao = A for i=1 to Np, where Np
%                                        is the number of RHP-poles.
%                                        Am = Ai, with i = Np.
%      Bi = Bi_1 - 2*Re(Pi)*Xpi*Ypi'*D   Bo = B for i=1 to Np, where Np
%                                        is the number of RHP-poles.
%                                        Bm = Bi, with i = Np.
%
%  BO(s) is given by:
%      BO(s) = PROD Bi(s), where Bi(s) = I + 2*Re(Pi)/(s-Pi') *Ypi*Ypi'
%                          the index runs from i=1 to Np.
%  Xpi and Ypi are the state and the output direction corresponding to the 
%  RHP-pole Pi for the system Gi_1(s) = |  Ai_1  |  Bi_1  |
%                                       |--------|--------|
%                                       |    C   |    D   |
%
%
%  Written  by: Kjetil Havre,  30/1 - 1996, e-mail: kjetil@ife.no
%
function [GO, BO, Yp, Xp, Pf] = opf( G, P, epp, lv )
   if nargin < 2
      disp('usage: opf(G, Z)')
      GO = G; BO = eye( ny );
      Yp = []; Xp = []; Pf = [];
      return
   end
   if nargin < 3
      epp = 1e-12;
   end
   if nargin < 4
      lv = 1;
   end
  
   [mtype,ny,nu,nx] = minfo(G);
   if mtype ~= 'syst'
      disp('Error: System matrix G is required.');
      disp('usage: opf(G, Z)');
      GO = G; BO = eye( ny );
      Yp = []; Xp = []; Pf = [];
      return
   end

   [np, nc] = size(P);
   if nc > 1                % Vector  is horizontal.
      P = vtp(P);
      [np, nc] = size(P);
   end

   if np < 1
      disp('Message: No poles?');
      GO = G; BO = eye(ny);
      Yp = []; Xp = []; Pf = [];
      return     
   end
   
   p = P(1);
   [up, yp, xip, xop, err] = pdfst( G, p, epp );   % Compute directions.
   if err != 0 
      disp( ['Warning: Requested pole P(', int2str(lv),'): ',...
	           num2str(P(1)), ' is not in system.' ] );      
      disp(['Error when computing pole direction, ',...
            'error number: ', int2str(err), '.']);
      GO = G; BO = eye(ny); Yp = []; Xp = []; Pf = [];
      return;
   else
      [A,B,C,D] = unpck( G );
      An  = A - 2*real( p )*xop*yp'*C;
      Bn  = B - 2*real( p )*xop*yp'*D;
      GO  = pck(An,Bn,C,D);
      Ab  = p;
      Bb  = 2*real(p)*yp';
      Cb  = yp;
      Db  = eye(ny);
      BO  = pck(Ab,Bb,Cb,Db);
      Yp  = yp;
      Xp  = xop;
      Pf  = p;
   end
   
   if np > 1   % Call recursive for factorizing additional RHP-zeros.
      [Gn, Bn, Yn, Xn, Pn] = opf(GO, P(2:np), epp, lv+1);
      GO = Gn;
      BO = mmult(BO, Bn);
      Yp = [Yp Yn];
      Xp = [Xp Xn];
      Pf = [Pf; Pn];
   end
end

  An  = A - 2*real( p )*xop*yp'*C;
      Bn  = B - 2*real( p )*xop*yp'*D;
     orpf.m000064400136200001130000000032140606457650600126610ustar00kjetilentek00000000000000%
%  Factorize one single real RHP-pole p in output factorization.
%
%  [GI, BI, Xf, Uf, Zf] = orpf( G, p, epp )
%  
%
%  Written by: Kjetil Havre,  16/12 - 1995, e-mail: kjetil@ife.no                         
%
%  See also: IZDE, INPCOMPFACT, OUTREALFACT, OUTCOMPFACT, INPFACT, OUTFACT.
%

function [GO, BO, Yp, Xp, Pf] = orpf( G, p, epp )

   if nargin < 2
      disp('Usage: orpf(G, p)');
      return
   end
   if nargin < 3
      epp = 1e-12;
   end
   
   [mtype,ny,nu,nx] = minfo(G);
   if mtype ~= 'syst'
      disp('Error: System matrix G is required.');
      disp('usage: inprealfact(G, Z)');
      GO = G; BO = eye( ny ); Xp = []; Up = []; Pf = [];
      return
   end

   [np, nc] = size(p);
   p = p(1);
   if np < 1
      disp('Message: No pole?');
      GO = G; BO = eye( ny ); Xp = []; Up = []; Pf = [];
      return     
   end
   if np > 1
      disp('Warning: Only the first pole is factorized.')
   end
   
   if ~isreal(p)
      disp(['Error: Zero is supposed to be real Z: ', num2str(z), '.']);
      GI = G; BI = eye( nu ); Xz = []; Uz = []; Zf = [];
      return;
   end

   [Up, Yp, XI, XO, err] = pdfst(G, p);
   
   pr = real(p); Ypr = real(Yp); Xpr = real(XO);

   err = norm([XO; Yp] - [Xpr; Yp]);
   if err > epp
      disp('Error: State and input direction is not real.');
      GO = G; BO = eye( ny ); Xp = []; Up = []; Pf = [];
      return;      
   end
   Yp = Ypr; Xp = Xpr;

   [A,B,C,D] = unpck( G );
   An = A - 2*real(p)*Xp*Yp'*C;
   Bn = B - 2*real(p)*Xp*Yp'*D;
   GO  = pck(An,Bn,C,D);

   Ab  = -p; Bb = 2*p*vtp(Yp); Cb  = -Yp; Db = eye(ny);
   BO  = minv(pck(Ab,Bb,Cb,Db));
   Pf  = p; 
end


% This works OK as long as D is zero.
real(Yp); Xpr = real(XO);

   err = norm([XO; Yp] - [Xpr; Yp]);
   if err > epp
      disp('Error: State and input direction is not real.');
      GO = G; BO = eye( ny ); Xp = []; Up = []; Pf = [];
      return;      
   end
   Yp = Ypr; Xp = Xpr;

   [A,B,C,D] = unpck( G );
   An = A - 2*real(p)*Xp*Yp'*C;
   Bn = B - 2*real(p)*Xp*Yp'*D;
   GO  = pck(An,Bn,C,D);

   Ab outcomp.m000064400136200001130000000017760602752131200133740ustar00kjetilentek00000000000000%
%  c = outcomp( G, op );
%
%  Returns the complement set of outputs described in op. All input indices
%  not in op. Example: Number of outputs is four and op=[2 3]; then c=[1 4];
% 
%  Inputs:  G  - Matrix describing the system, can bd system, constant or 
%                varying.
%           op - Set of outputs to be complemented.
%
%  Outputs: c  - complement set.
%
%  Written by: Kjetil Havre, e-mail: kjetil@ife.no
%  
%  Note: This is for sure not the most effective way to find the complemet 
%        set in matlab.
%

function c = outcomp( G, op )
   c = [];
   [mt, rd, cd, nd] = minfo(G);
   no = max(size(op));
% Added 9/8-95 KH
   if strcmp(op,':') == 1 
      return
   end
% to this spot 9/8-95.
   if no <= rd
      opc = [];
      for i=1:rd,                     
	 found = 0;
	 for j=1:no,
	    if op(j) == i;
	       found = 1;
	    end
	 end
	 if found == 0
	    opc = [opc i];
	 end
      end
      c = opc;
   else
      disp('Number of outputs is larger than rows' );
      return;
   end
end
  outcprfact.m000064400136200001130000000042740606404524000140560ustar00kjetilentek00000000000000%
%   Function OUTCPRFACT finds factorization such that G(s) = N(s)*M(s)^{-1},
%   where N(jw)^H*N(jw) = I, and both N and M are stable. That means
%   all RHP-zeros in N and all RHP poles as zeros in M. 
%   System G must be square and have no jw-axis poles. 
%
%   Reference: Morari M. ``Robust Process Control'', Prentice Hall 1989.
%

function [N, M, Mi, Q, R, E, X, Xerr] = outcprfact( G, deps, epp )

   if nargin < 2
      deps = 1e-6;
   end
   if nargin < 3
      epp = 1e-12;
   end

   [mt, ny, nu, nx] = minfo(G);
   N = []; M =[];

   if mt ~= 'syst'
      disp('Error: System matrix is requierd');
      return
   end
   
   if ny ~= nu
      disp('Error: Square system matrix is required.');
      return
   end
   
   
   [A,B,C,D] = unpck(G);
   Dn = D;
   SD = svd(D);
   if SD(1,1) < deps
      disp(['Warning: adding diagonal part to D of size: ', num2str(deps), 'I.']);
      Dn = D + deps*eye(size(D));
   end
   
   [Q,R,E] = qr(Dn)
   
   Q'*Q
   
   SQ = svd(vtp(Q)*Q-eye(nu,nu));
   if SQ(1,1) > epp
      disp('Warning: QR decomposition of D is not acheived.');
      return
   end
   
   Aric = A - B*inv(R)*vtp(Q)*C;
   Bric = B*inv(R)*vtp(B*inv(R));
   Cric = zeros(nx,nx);
   

   X = are(Aric, Bric, Cric);
   Xerr = vtp(Aric)*X + X*Aric - X*Bric*X + Cric;

%   SXerr = svd(Xerr)
%
%   Control system Toolbox: tolerance not an input argument.
%   Resulting error is in order od 1e-3.
%

%   H = [ Aric -Bric; -Cric -A' ];
%   eig(H)
%   [x1, x2, fail, reig_min] = ric_eig(H)
%   if fail==1
%      disp('Error: jw - axis eigenvalues.')
%      return
%   end
%   X = x2/x1;
%   Xerr = vtp(Aric)*X + X*Aric - X*Bric*X + Cric;
%   SXerr = svd(Xerr)
%
%
%   Mu - tools are not usable, it returns that H has jw axis 
%        eigenvalues althoug reig_min = 7.5326e-02.
%
%   Conclusion: We need robust control toolbox. 
%

   F = vtp(Q)*C + vtp(B*inv(R))*X;
%   F1 = Q'*C+inv(R)'*B'*X;
%   EF = F-F1
   An = A-B*inv(R)*F;
   Bn = B*inv(R);
   Cn = C-Q*F;
   Dn = Q;
   N = pck(An,Bn,Cn,Dn);
   
   Ami = A;
   Bmi = B;
   Cmi = F;
   Dmi = R;
   Mi = pck(Ami, Bmi, Cmi, Dmi);
   
   Am = Ami-Bmi*inv(Dmi)*Cmi;
   Bm = -Bmi*inv(Dmi);
   Cm = inv(Dmi)*Cmi;
   Dm = inv(Dmi);
   M  = pck(Am, Bm, Cm, Dm);
   
end
hat H has jw axis 
%        eigenvalues althoug reig_min = 7.5326e-02.
%
%   Conclusion: We need robust control toolbox. 
%

   F = vtp(Q)*C + vtp(B*inv(R))*X;
%   F1 = Q'*C+inv(R)'*B'*X;
%   EF = F-F1
   An = A-B*inv(R)*F;
   Bn = B*inv(R);
   Cn = C-Q*F;
   Dn = Q;
   N = pck(An,Bn,Cn,Dn);
   
   Ami = A;
   Bmi = B;
   ozde.m000064400136200001130000000150400611041764400126410ustar00kjetilentek00000000000000% function [Z, Y, X] = ozde(sys,epp)
%
%   This Output-Zero-Direction-"through genaralized Eigenvalue" decomposition
%   OZDE function is a modification of the szeros function contained in mu 
%   - toolbox. The modification consists of returning the output zero 
%   directions and the state zero directions in addition to the zeros.
%   The output zero directions are defined as:
%              y'*G(z) = 0,
%   where s = z is a zero of G(s) and ' denotes conjugate transposed.
%   This is done by solving the generalized eigenvalue problem of the 
%   transposed system:
%
%           | x' | y' | * | A-Iz |  B  |  =  |  0 | 0  |
%                         |------------|
%                         |  C   |  D  |
%
%   are solved by genaralized eigenvalues:
%
%           | A"-Iz |  C"  | * | xi |  =  | 0 |
%           |--------------|   |----|     |---|
%           |  B"   |  D"  |   | yi |     | 0 |
%
%   x = conj(xi); y = conj(yi);
%
%   The prime ' denotes the conjugate transposed and " denotes the ordinary 
%   transposed.
%
%   OZDE 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 eigenvalue 
%   problems. EPP (optional) defines if the difference between two generalized
%   eigenvalues is small. OZDE also finds the output y and the state x 
%   directions of the zeros.
%
%   The output zero directions are stored as column vectors in Y, and the y's 
%   are normalized so that Y'*Y = I. The state zero directions are stored as
%   columns in X. The degree of freedom to normalize the generalized eigen-
%   vectors is used to normalize the y part of the vectors. So the length of 
%   x is not equal to one. Each column in Y and X corresponds to the element
%   in Z with same place.

%   Note that when the number of outputs is larger than the number of inputs,
%   the output zero direction is not complete. A bit clearer: If z is a zero 
%   of a non-square plant with number of outputs greater than number of inputs,
%   the zero direction is not a line but a surface. As an example, consider 
%   G(s) with dimensions 3x2 (3 rows and 2 columns). Let s=z be a zero of 
%   G(s) such that Yz'*G(z) = [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.
%   y 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 zero space on the output 
%   completely. Two basis vectors are requiered. 
%
%   Modification for square systems was made by:      Kjetil Havre 24/4-1995.
%   Modification for non square systems was made by:  Kjetil Havre 2/5 -1995.
%   Including the state zero dircetions was made by:  Kjetil Havre 15/5-1995.
%   Modified so that first element of U(:,i) is real: Kjetil Havre  3/2-1996. 
%   Comments, corrections and malfunctions, can be e-mailed to: 
%                             kjetil@ife.no
%
%   See also: EIG, SZEROS, IZDE and SPOLES. 
%   Algorithm based on Laub & Moore 1978 paper, Automatica
%

function [Z, Y, X] = ozde(sys,epp)
  if nargin < 1
    disp('usage: szeros(sys)')
    return
  end
  if nargin == 1
    epp = eps;
  end
  [mtype,ny,nu,nx] = minfo(sys);
  if mtype == 'syst'
    [a,b,c,d] = unpck(sys);
    if nx == 0
      disp('SYSTEM has no states')
    end
    sysu = [a b; c d];
    sysu = vtp( sysu );   % Transpose system.
    
% 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);

      % Add something here to check for not a number or infinite.
      % remove corresponding vectors.
      % Also remove top part corresponding to the states
      % and normalize bottom part. KH 21/4 - 1995.

      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 in x and y.
      vx = vech2(1:nx, : );
      vy = vech2(nx+1:nx+ny,:);
% Normalize columns.
      [nvr, nvc] = size( vy ); 
      for i=1:nvc,
          nrmy   = norm(vy(:,i));
	  Inz = find( abs(vy(:,i)) > 1000*epp );
          X(:,i) = conj(vx(:,i)/nrmy * exp( -angle(vy(Inz(1),i))*sqrt(-1) ) );
          Y(:,i) = conj(vy(:,i)/nrmy * exp( -angle(vy(Inz(1),i))*sqrt(-1) ) );
      end
      
    else                                  % Non-square systems 
      nrm = norm(sysu,1);
      if nu > ny
        x1 = [ sysu  nrm*(rand(nx+nu,nu-ny)-.5)];
        x2 = [ sysu  nrm*(rand(nx+nu,nu-ny)-.5)];
      else
        x1 = [ sysu; nrm*(rand(ny-nu,nx+ny)-.5)];
        x2 = [ sysu; nrm*(rand(ny-nu,nx+ny)-.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      
      if isempty(vech3)
         Z = []; Y = []; X = [];
	 return;
      end
% Split columns in x and y.
      vx = vech3(1:nx,:);
      vy = vech3(nx+1:nx+ny,:);
% Normalize columns.
      [nvr, nvc] = size( vy ); 
      for i=1:nvc,
          nrmy   = norm(vy(:,i));
	  Inz = find( abs(vy(:,i)) > 1000*epp );
          X(:,i) = conj(vx(:,i)/nrmy * exp( -angle(vy(Inz(1),i))*sqrt(-1) ) );
          Y(:,i) = conj(vy(:,i)/nrmy * exp( -angle(vy(Inz(1),i))*sqrt(-1) ) );
      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
%
+ny,:);
% Normalize columns.
      [nvr, nvc] = size( vy ); 
      for i=1:nvc,
          nrmy   = norm(vy(:,i));
	  Inz = find( abs(vy(:,i)) > 1000*epp );
          X(:,i) = conj(vx(:,i)/nrmy * exp( -angle(vy(Inz(1),i))*sqrt(-1) ) );
          Y(:,i) = conj(vy(:,i)/nrmy * exp( -angle(vy(Inz(1),i))*sqrt(-1) ) );
      end
      
    end

  else
     error('matrix is not a SYSTEM matrix')
     return
  end 
%
% Copyright MUSYN INC 1991,  All Rights Reserved
% Copyright MUSYN Iozf.m000064400136200001130000000110020610467100500124640ustar00kjetilentek00000000000000%  function [GO, BO, Yz, Xz, Zf] = ozf( G, Z, epp )
%  
%  The OZF function is used to obtain an output factorization of RHP-zeros 
%  so that G(s) = BO(s)*GO(s).
%  The input argument Z is a vector containing the RHP-zeros. The optional
%  EPP argument is used to decide the difference between the zeros given
%  in Z and the zeros calculated by the OZF function. If the function does
%  not calculate any zeros Zf so that abs(Z(i) - Zf) < EPP the function 
%  skipps Z(i) in the factorization. The default EPP value is 1e-12.
%  The ODZE function is used to obtain the zero state direction Xzi and 
%  the output zero direction Yzi for the RHP-zero Z(i).
%  
%  BO(s) is an all pass transfer function for s=jw. All singular values
%  are equal to one, and the absolute value of all eigenvalues are 
%  equal to one. 
%
%  GO is constructed as follows:
%           G(s) = |  A  |  B  |      GO(s) = |  A   |   B  |
%                  |  C  |  D  |              |  Cm  |   D  |
%  where Cm is obtained by running the following formula repeatedly:
%      Ci = Ci_1 - 2*Re(Zi)*Yzi*Xzi',   Co = C for i=1 to Nz, where Nz
%                                       is the number of RHP - zeros.
%                                       Cm = Ci, with i = Nz.
%
%  BO(s) is given by:
%      BO(s) = PROD Bi(s), where Bi(s) = I - 2*re(Zi)/(s+Zi') *Yzi*Yzi'
%                          the index runs from i=Nz to one.
%  Xzi and Yzi are the state and the output direction corresponding to the 
%  RHP-zero Zi for the system Gi_1(s) = |    A   |  B  |
%                                       |--------|-----|
%                                       |  Ci_1  |  D  |
%
%  Possible malfunctions: Multiple zeros located at the same place in the 
%                         complex plane.
%                         RHP-zeros on the jw axis, pure imaginary zeros.
%                         The function has not been tested for number of 
%                         zeros Nz greater than one. However, it should
%                         work OK with Nz > 1. Can also be used for complex
%                         zeros.
%
%  Written  by: Kjetil Havre,  15/5 - 1995, e-mail: kjetil@ife.no
%  Modified by: Kjetil Havre,  30/1 - 1996, e-mail: kjetil@ife.no
%      Modified to also return the state direction.
%
function [GO, BO, Yz, Xz, Zf] = ozf( G, Z, epp, lv )
   if nargin < 2
      disp('usage: ozf(G, Z)')
      GO = G; BO = eye( ny );
      Yz = []; Xz = []; Zf = [];
      return
   end
   if nargin < 3
      epp = 1e-12;
   end
   if nargin < 4
      lv = 1;
   end
  
   [mtype,ny,nu,nx] = minfo(G);
   if mtype ~= 'syst'
      disp('Error: System matrix G is required.');
      disp('usage: ozf(G, Z)');
      GO = G; BO = eye( ny );
      Yz = []; Xz = []; Zf = [];
      return
   end

   [nz, nc] = size(Z);
   if nc > 1                % Vector  is horizontal.
      Z = vtp(Z);
      [nz, nc] = size(Z);
   end

   if nz < 1
      disp('Message: No zeros?');
      GO = G; BO = eye(ny);
      Yz = []; Xz = []; Zf = [];
      return     
   end
   
   z = Z(1);
   [Ze, yiz, xiz] = ozde( G );   % Compute directions.
   
   Iz = find( abs(Ze-z) < epp );
   if isempty(Iz)
      disp( ['Warning: Requested zero Z(', int2str(lv),'): ',...
	           num2str(Z(1)), ' is not in system.' ] );      
      GO = G; BO = eye(ny); Yz = []; Xz = []; Zf = [];
   else
      % Iz contains the index to the zero in Ze. 
      % Factorize this zero.
      
      [A,B,C,D] = unpck( G );
      
      H = [A-z*eye(size(A))  B;   C   D];
      errvk  = [xiz(:,Iz)' yiz(:,Iz)']*H;
      nerrvk = norm(errvk);
      if norm(errvk) > epp
         disp('Error: Could not find generalized eigenvalues.')
         GO = G; BO = eye(ny); Yz = []; Xz = []; Zf = [];
      end

      if isreal(z)
         xiz(:,Iz) = real(xiz(:,Iz));
         yiz(:,Iz) = real(yiz(:,Iz));
	 Ze(Iz) = real(Ze(Iz));
      end

      Cn  = C - 2*real( Ze(Iz) )*yiz(:,Iz)*xiz(:,Iz)';
      GO  = pck(A,B,Cn,D);
      Ab  = -conj(Ze(Iz));
      Bb  = 2*real(Ze(Iz))*yiz(:,Iz)';
      Cb  = -yiz(:,Iz);
      Db  = eye(ny);
      BO  = pck(Ab,Bb,Cb,Db);
%      br  = nd2sys( [1], [1 conj(Ze(Iz))], 2*real(Ze(Iz)) );
%      brm = mmult( yiz(:,Iz), br, yiz(:,Iz)' );
%      BO = msub( eye(ny), brm );
      Yz  = yiz(:,Iz);
      Xz  = xiz(:,Iz);
      Zf  = Ze(Iz);
   end
   
   if nz > 1   % Call recursive for factorizing additional RHP-zeros.
      [Gn, Bn, Yn, Xn, Zn] = ozf(GO, Z(2:nz), epp, lv+1);
      GO = Gn;
      BO = mmult(BO, Bn);
      Yz = [Yz Yn];
      Xz = [Xz Xn];
      Zf = [Zf; Zn];
   end
end


   
   
   
   
   
   
   
   
   
   
   
);
%      br  = nd2sys( [1], [1 conj(Ze(Iz))], 2*real(Ze(Iz)) );
%      brm = mmult( yiz(:,Iz), br, yiz(:,Iz)' );
%      BO = msub( eye(ny), brm );
      Yz  = yiz(:,Iz);
      Xz  = xiz(:,Iz);
      Zf  = Ze(Iz);
   end
   
   if nz > 1   % Call recursive for factorizing additional RHP-zeros.
      [Gn, Bn, Yn, Xn, Zn] = ozf(GO, Z(2:nz), epp, lv+1);
      GO = Gn;
      BO = mmult(BO, Bn);
      Yz = [Yz Yn];
      Xz = [Xz Xn];
      Zf = [Zf; Zn];
   end
end


   
   
   
   
   
   
   
   
   
   
  partcont.m000064400136200001130000000065430602751757300135530ustar00kjetilentek00000000000000% [Pd, Pr, Pu, G11, G12, G21, G22, Gd1, Gd2, R1, R2] =...
%                              partcont( G, Gd, R, op, ip, dip )
% 
% Computation of frequency dependent partial gains.
% 
% Outputs: Pd   - Partial disturbance gain.
%          Pr   - Partial reference gain.
%          Pu   - Partial control gain of unused inputs.
%          G11  - Subsystem of uncontrolled outputs and unused inputs.
%          G12  - Subsystem of uncontrolled outputs and used contol inputs.
%          G21  - Subsystem of controlled outputs and unused inputs.
%          G22  - Subsystem of controlled outputs and used inputs.
%          Gd1  - Disturbance model for uncontrolled outputs.
%          Gd2  - Disturbance model for ontrolled outputs.
%          R1   - Part of R matrix corrosponding to uncontrolled outputs.
%          R2   - Part of R matrix corrosponding to controlled outputs.
%
% Inputs:  G    - System as a varying matrix,
%          Gd   - Disturbance system as a varying matrix or constant matrix,
%          R    - Reference matrix as a varying matrix or constant matrix,
%          op   - Vector specifying the outputs,
%          ip   - Vector specifying the inputs,
%          dinp - Vector specifying the disturbance inputs.
%
% The function uses Mu-tools and the following three functions:
%      outcomp, inpcomp and selc.
%
% Written by: Kjetil Havre, e-mail: kjetil@ife.no
%
%

function [Pd, Pr, Pu, G11, G12, G21, G22, Gd1, Gd2, R1, R2] = partcont( G, Gd, R, op, ip, dip )
   [mt, rd, cd, n]     = minfo( G );
   [mtd, rdd, cdd, nd] = minfo( Gd );
   [mtr, rdr, cdr, nr] = minfo( R );
   if rd ~= rdd
      disp( 'Number of rows in G and Gd does not agree.' );
      return
   end

   
   if strcmp( mt, 'vary' ) ~= 1
       disp( 'Only implemented for varying matrices.' );
       return;
   end
   
   if nargin < 6
      dip = ':';
   end
   if nargin < 5
     ip = ':';
   end
   if nargin < 4           %  Perfect control 
      op = ':';
   end
      
   if nargin < 2
      disp( 'Usage: [Pd, Pr, Pu] = partcont( G, Gd, op, inp, dip )' )
      return
   end
      
   opc  = outcomp( G, op );   % Find the complement of the outputset.
   ipc  = inpcomp( G, ip );   % Find the complement of the input set.
   nipc = max(size(ipc));
   dipc = inpcomp( Gd, dip ); % Find the complement of the dist. set.
   if strcmp(mtd, 'cons') == 1
     Gd   = selc( Gd, ':', dip );
     Gd1  = selc( Gd, opc, dip );
     Gd2  = selc( Gd,  op, dip );
   else
     Gd   = sel( Gd, ':', dip );
     Gd1  = sel( Gd, opc, dip );
     Gd2  = sel( Gd,  op, dip );
   end
   if strcmp(mtr, 'cons') == 1
      R1 = selc(R, opc, opc);
      R2 = selc(R, op, op );
   else
      R1 = sel(R, opc, opc);
      R2 = sel(R, op, op );
   end
   G22  = sel(  G,  op, ip   );
   G12  = sel(  G, opc,  ip  );
   if nipc > 0
      G11  = sel(  G, opc, ipc  );
      G21  = sel(  G,  op, ipc  );
   else
      G11 = [];
      G21 = [];
   end
	 
   [mt, nr22, nc22, nx22] = minfo( G22 );
   if nr22 < rd
      if nr22 == nc22
         G22i = minv(G22);
      else
	 disp('WARNING: Not implemented for G22 non-square.' )
	 return
      end
      Pd = msub( Gd1, mmult(G12,G22i,Gd2) );
      if nipc > 0
         Pu = msub( G11, mmult(G12,G22i,G21) );
      else
         Pu = [];
      end
      Pr = mmult(G12,G22i,R2);
   else                                    %  Perfect control
      disp(' Not implemented.' );
      return
   end      
end
;
   end
	 
   [mt, nr22, nc22, nx22] = minfo( G22 );
   if nr22 < rd
      if nr22 == nc22
         G22i = minv(G22);
      else
	 disp('WARNING: Not implempartcontsys.m000064400136200001130000000063070601506143600142750ustar00kjetilentek00000000000000% [Pd, Pr, Pu, G11, G12, G21, G22, Gd1, Gd2, R1, R2] = 
%                              partcont( G, Gd, R, op, ip, dip )
% 
% Computation of partial gains.
% 
% Outputs: Pd   - Partial disturbance gain, Pd = Gd1 - G12*minv(G22)*Gd2.
%          Pr   - Partial reference gain,   Pr =       G12*minv(G22)*R2.
%          Pu   - Partial control gain,     Pu = G11 - G12*minv(G22)*G21.
%          G11  - Subsystem of uncontrolled outputs and unused inputs.
%          G12  - Subsystem of uncontrolled outputs and used contol inputs.
%          G21  - Subsystem of controlled outputs and unused inputs.
%          G22  - Subsystem of controlled outputs and used inputs.
%
% Inputs:  G    - System as a varying matrix.
%          Gd   - Disturbance system as a varying matrix or constant matrix.
%          R    - Reference matrix as a varying matrix or constant matrix.
%          op   - Vector specifying the outputs.
%          ip   - Vector specifying the inputs.
%          dinp - Vector specifying the disturbance inputs.
%
% Rearranges G = | G11  G12 |,  R = | R1   0 |,  Gd = | Gd1 |
%                | G21  G22 |       |  0  R2 |        | Gd2 |
% according to op, ip and dip.
% 
%
% Written by: Kjetil Havre, e-mail: kjetil@ife.no
%
%
function [Pd, Pr, Pu, G11, G12, G21, G22, Gd1, Gd2, R1, R2] = partcont( G, Gd, R, op, ip, dip )
   [mt, rd, cd, n] = minfo( G );
   [mtd, rdd, cdd, nd] = minfo( Gd );
   [mtr, rdr, cdr, nr] = minfo( R );
   if rd ~= rdd
      disp( 'Number of rows in G and Gd does not agree.' );
      return
   end

   
   if strcmp( mt, 'vary' ) ~= 1
       disp( 'Only implemented for varying matrices.' );
       return;
   end
   
   if nargin < 6
      dip = ':';
   end
   if nargin < 5
     ip = ':';
   end
   if nargin < 4           %  Perfect control 
      op = ':';
   end
      
   if nargin < 2
      disp( 'Usage: [Pd, Pr, Pu] = partcont( G, Gd, op, inp, dip )' )
      return
   end
      
   opc  = outcomp( G, op );   % Find the complement of the outputset.
   ipc  = inpcomp( G, ip );   % Find the complement of the input set.
   nipc = max(size(ipc));
   dipc = inpcomp( Gd, dip ); % Find the complement of the dist. set.
   if strcmp(mtd, 'cons') == 1
     Gd   = selc( Gd, ':', dip );
     Gd1  = selc( Gd, opc, dip );
     Gd2  = selc( Gd,  op, dip );
   else
     Gd   = sel( Gd, ':', dip );
     Gd1  = sel( Gd, opc, dip );
     Gd2  = sel( Gd,  op, dip );
   end
   if strcmp(mtr, 'cons') == 1
      R1 = selc(R, opc, opc);
      R2 = selc(R, op, op );
   else
      R1 = sel(R, opc, opc);
      R2 = sel(R, op, op );
   end
   G22  = sel(  G,  op, ip   );
   G12  = sel(  G, opc,  ip  );
   if nipc > 0
      G11  = sel(  G, opc, ipc  );
      G21  = sel(  G,  op, ipc  );
   else
      G11 = [];
      G21 = [];
   end
	 
   [mt, nr22, nc22, nx22] = minfo( G22 );
   if nr22 < rd
      if nr22 == nc22
         G22i = minv(G22);
      else
	 disp('WARNING: Not implemented for G22 non-square.' )
	 return
      end
      Pd = msub( Gd1, mmult(G12,G22i,Gd2) );
      if nipc > 0
         Pu = msub( G11, mmult(G12,G22i,G21) );
      else
         Pu = [];
      end
      Pr = mmult(G12,G22i,R2);
   else                                    %  Perfect control
      disp(' Not implemented.' );
      return
   end      
end
;
   end
	 
   [mt, nr22, nc22, nx22] = minfo( G22 );
   if nr22 < rd
      if nr22 == nc22
         G22i = minv(G22);
      else
	 disp('WARNING: Not implemented for G22 non-square.' )
	 return
      end
      Pd = msub( Gd1, mmult(G12,G22i,Gd2) );
      if nipc > 0
         Pu = msub( G11, mmult(G12,G22i,G21) pdfst.m000064400136200001130000000154160611061675000130260ustar00kjetilentek00000000000000%   [I, O, XI, XO, err] = pdfst(G, P, SST, NDL)
%
%   The ``Pole-Direction-From-State'' pdfst-function computes the 
%   input and output pole directions through the use of Singular 
%   Value Decomposition of (pI-A) where p is a pole.
%   If smallest singular value of (pI-A) is smaller than SST then 
%   the output pole direction yp can be computed as yp = C*xr,
%   where xr is the eigenvektor of A corresponding to p, A*xr = p*xr. 
%   We note that xr is similar to the input singular vektor 
%   corresponding to the zero singular value.
%   The input pole direction up is computed from the transposed system 
%   i.e.  up = B.'*xl, where xl is the solution to the eigenvalue 
%   problem A^T*xl = p*xl, which corresponds to the output singular 
%   vector with zero singular value of (pI-A).
%   To calculate the zero direction in the state space, the function 
%   uses singular decomposition of (pI-A) and not the the eigen value 
%   formulation. The pole directions are stored as column vectors in 
%   I and O when the smallest singular value is less than SST.
%   If the smallest singular value is larger than SST, the zero vectors
%   are stored in the columns of I and O.
%   A warning message is printed if also the second singular value is 
%   greater than SST. 
%
%   Note! PDFST is numerically more robust than PDSVD.
%   
%   
%   Inputs:  G   - system matrix.
%            P   - vector of poles.
%            SST - tolerance (optional) for singular values, 
%                  default value is 1e-12.
%            NDL - if the norm of the D matrix is large, the directions
%                  computed in PDFST may be wrong. A warning message is
%                  displayed if || D ||_2 >= NDL. Default value is 
%                  1000/SST.
%
%   Outputs: I   - Matrix containing the input directions, 
%                  column i corresponds to pole P(i).
%            O   - Matrix containing the output directions, 
%                  column i corresponds to pole P(i).
%            XI  - Input state directions.
%            XO  - Output state directions.
%            err - Error indicator:
%                     err =  0, all is OK.
%                     err =  1, No state direction with gain smaller than 
%                               SST.
%                     err = -1, One basis vector is not enough to describe 
%                               the subspace of states with gains smaller 
%                               than SST.
%                     err = -2, One basis vector is not enough
%                               to describe the subspace with 
%                               gains larger than EPP.
%
%   Written by: Kjetil Havre, 30/11-1996, e-mail: kjetil@ife.no
%
%   See also:  IZDE, OZDE, ZDSVD and PDSVD.
%
%   Reference: Havre K. and S. Skogestad, 1995, ``Effect of RHP Zeros and
%              Poles on Performance in Multivariable Systems''.
%
   function [I, O, XI, XO, err] = pdfst(G, P, sst, ndl)

   I = []; O = []; XI = []; XO = []; err = 0;
   [mt,ny,nu, nx] = minfo(G);
   
   if strcmp(mt, 'syst')==0
      disp( 'System matrix G is required, usage: [I, O, XI, XO, err] = pdfst(G, P, SST, NDL)' );
      return
   end
   if nargin < 2
       P = spoles(G);
%      disp( 'Usage: [I, O, XI, XO, err] = pdsvd(G, P, SST, NDL)' );
%      return
   end   
   if nargin < 3
      sst = 1e-12;
   end
   if nargin < 4
      ndl = 1000/sst;
   end
   [A, B, C, D] = unpck(G);
   if norm(D) > ndl
      disp('Warning: Directions may be inaccurate due to large effect from D.');
      err = -2;
   end
   
   for i=1:max(size(P))
      [U, S, V] = svd( A-eye(nx)*P(i) );
      if S(nx,nx) < sst
         if nx > 1
	    if S(nx-1,nx-1) < sst
               disp(['Warning: Second smallest singular value is smaller than SST: ', num2str(sst),','])
               disp(['         S(',int2str(nx-1),',',int2str(nx-1),'): ', num2str(s(nr-1,nr-1)),'.'])
               err = -1;
	    end
	 end
	 
%         [Vr, Dr] = eig(A);
%	 [Vl, Dl] = eig(A.');
%	 Ir = find(diag(Dr)==P(i))
%	 Il = find(diag(Dl)==P(i))
%	 XO(:,i) = Vr(:,Ir);
%	 XI(:,i) = Vl(:,Il);
         	 
	 XO(:,i) = V(:,nx);
	 yph = C*XO(:,i);  O(:,i) = yph/norm(yph); XO(:,i) = XO(:,i)/norm(yph);	 
         inz = find(abs(O(:,i))>sst); vh = angle(O(inz(1),i));
	 O(:,i) = O(:,i)*exp(-sqrt(-1)*vh); XO(:,i)=XO(:,i)*exp(-sqrt(-1)*vh);

%         XI(:,i) = U(:,nx); % I need this to get the error small 
         XI(:,i) = U(:,nx);
	 uph=B'*XI(:,i); 
%
%        Comment: Almost always is the B matrix real, however when
%                 facorizing first on RHP pole then the B matrix 
%                 becomes complex if the pole is complex.
%        in order that up^H*G^{-1}(p) = 0 we need to define
%        u_p = B^H x_p
%        This make on a difference when B is complex.
%        Kjetil Havre 15/2 - 1996
%
         I(:,i) = uph/norm(uph); XI(:,i)=XI(:,i)/norm(uph);
         inz=find(abs(I(:,i))>sst); vh = angle(I(inz(1),i));
	 I(:,i)=I(:,i)*exp(-sqrt(-1)*vh); XI(:,i)=XI(:,i)*exp(-sqrt(-1)*vh);
	 
	 wp = P(i)/sqrt(-1)+max(10*eps,0.01*sst);
	 Gp = frsp(G, wp);
	 ypn = vnorm( mmult(vpinv(Gp), O(:,i)) );
	 upn = vnorm( mmult(I(:,i)',vpinv(Gp)) );
	 
	 if ypn(1,1) > 100*sst
	    ypn= ypn
	    Ph = P(i)
	    disp(['Warning:   Norm || G^{-1}(p)*yp ||_2 is larger than 10*sst: ',...
	          num2str(ypn(1,1)), '.']);
%	    disp('         Storing zero vector for output direction.');
%	    O(:,i) = zeros(ny,1);
%	    err = 2;
	 end
	 
	 if upn(1,1) > 100*sst
	    upn = upn
	    P(i) = P(i)
	    isrB = isreal(B)
	    if isrB
	       disp('B is real')
	    else 
	       disp('B is not real')
	       B = B
	    end
	    disp(['Warning:   Norm || up^H*G^{-1}(p) ||_2 is larger than 10*sst: ',...
	          num2str(upn(1,1)), '.']);
%	    disp('         Storing zero vector for input direction.');
%	    I(:,i) = zeros(nu,1);
%	    err = 3;
	 end
%
%   This check is removed due to numerical problems.
%
	 
	 Amp  = A-eye(nx)*P(i);
	 xon  = norm(Amp*XO(:,i));
	 xin  = norm(XI(:,i)'*Amp);
	 
	 if xon > 10*sst
	    disp(['Error:   Norm || (A-p*I)*xo ||_2 is larger than 10*sst: ',...
	          num2str(xon(1,1)), '.']);
	    disp('         Storing zero vector for output state direction.');
	    XO(:,i) = zeros(nx,1);
	    err = 4;
	 end
	 
	 if xin > 10*sst
	    disp(['Error:   Norm || xi^H*(A-p*I) ||_2 is larger than 10*sst: ',...
	          num2str(xin(1,1)), '.']);
	    disp('         Storing zero vector for input state direction.');
	    XI(:,i) = zeros(nx,1);
	    err = 4;
	 end

%
%   Kjetil Havre 16/12 - 1995.
%
	 
      else
         disp(['Error:   Matrix (I*p-A) is not singular for i=',...
	        int2str(i),', P(i)=', num2str(P(i)), '.'] );
         disp(['         Smallest singular value is: ', num2str(S(nx,nx)),'. Storing zero vectors.'] );
	 O(:,i)  = zeros(ny,1);  I(:,i)  = zeros(nu,i);
	 XO(:,i) = zeros(nx,1);  XI(:,i) = zeros(nx,1);
	 err = 1;
      end
   end
end

   
   Storing zero vector for input state direction.');
	    XI(:,i) = zeros(nx,1);
	    err = 4;
	 end

%
%   Kjetil Havre 16/12 - 1995.
%
	 
      else
         disp(['Error:   Matrix (I*p-A) is not singular for i=',...
	        int2str(i),', pdisvd.m000064400136200001130000000070550605707760200132060ustar00kjetilentek00000000000000%   [I, O, err] = pdisvd(G, P, )
%
%   The ``Pole-Direction-Inverse-SVD'' function computes the input 
%   and output pole directions through the use of Singular Value 
%   Decomposition of G^{-1}(p) where p is a pole.
%   The system matrix G(s) is inverted before evaluation of G(p).
%   If largest singular value of G^{-1}(p) is smaller than EPP then 
%   the corrosponding input and output singular vectors are stored 
%   in columns of I and O. Otherwise, zero vectors are stored in the 
%   corrosponding columns.
%   A warning message is printed if also the second singular 
%   value is grater than EPP.
%   
%   Inputs:  G   - system matrix.
%            P   - vector of poles.
%            epp - tolerance (optional) for singular values, 
%                  default value is 1/eps which most often is
%                  to large typical value around 1e14.
%
%   Outputs: I   - Matrix containing the input directions, 
%                  column i corresponds to pole P(i).
%            O   - Matrix containing the output directions, 
%                  column i corresponds to pole P(i).
%            err - Error indicator:
%                     err =  0, all is OK.
%                     err =  1, No direction with gain larger 
%                               than EPP.
%                     err = -1, One basis vector is not enough
%                               to describe the subspace with 
%                               gains larger than EPP.
%
%   Written by: Kjetil Havre, e-mail: kjetil@ife.no
%
%   See also: izde, ozde, zdsvd.
%
   function [I, O, err] = pdisvd(G, P, vtol, mst, maxitr, delts)

   I = []; O = []; err = 0;

   [mt,ny,nu, nx] = minfo(G);
   
   if strcmp(mt, 'syst')==0
      disp( 'System matrix G is required, usage: [I,O] = pdisvd(G, P, epp)' );
      return
   end
   
   if nargin < 2
      disp( 'Usage: [I,O] = pdisvd(G, P, epp)' );
      return
   end
   
   if nargin < 6
      delts = 100;
   end
   if nargin < 5
      maxitr = 14;
   end
   if nargin < 4
      mst = 1e-10;
   end
   if nargin < 3
      vtol = 1e-6;
   end

   [A,B,C,D] = unpck(G); Bm = B; Cm = C; Db = D;
   if ny ~= nu
      if ny > nu
	 Bm = [B zeros(nx, ny-nu)];   % Add inputs.
         Db = [D zeros(ny, ny-nu)];   % Add inputs.
      else
         Cm = [C; zeros(nu-ny, nx)];  % Add outputs.
         Db = [D; zeros(nu-ny, nu)];  % Add outputs.
      end
   end
      
   nm = max(ny,nu);

   for i=1:max(size(P))
      erc = maxitr;
      delt = delts;
      for j=1:maxitr
         if rank(Db) < nm
            Dm = Db + delt*10*vtol*eye(nm);
         else
            Dm = Db;
         end
	 
         Gm  = pck(A, Bm, Cm, Dm);
         Gi  = minv(Gm);
         wp  = P(i)/sqrt(-1);
         Gip = frsp(Gi, wp);
         [u,s,v] = svd(Gip(1:nm,1:nm));
	 if j > 1;
 	    ypo = yp; upo = up;
	 end
         yp = v(:,1);	 up = u(:,1);
	 if j > 1
            if (norm(yp-ypo) + norm(up-upo)) < 2*sqrt(nm)*vtol
	       erc = 0
	       break;
	    end
	 end
	 delt = delt/2
      end            
      
      % Check and save directions
      nr = min(nu,ny);
      if (s(nr,nr) > mst) | (erc~=0)
         disp('Error: Required tolerance not acheived.');
         I(:,i) = zeros(nu,1);
         O(:,i) = zeros(ny,1);
	 err = 1;
      else
	 if nr > 2 
	    if s(nr-1,nr-1) < mst
               disp(['Warning: Second smallest singular value is less than eep: ', num2str(epp),','])
               disp(['         S(',int2str(nr-1),',',int2str(nr-1),'): ', num2str(s(nr-1,nr-1)),'.'])
               err = -1;
	    end
	 end
         I(:,i) = up(1:nu);
         O(:,i) = yp;(1:ny);
      end
   end
end

   
0)
         disp('Error: Required tolerance not acheived.');
         I(:,i) = zeros(nu,1);
         O(:,i) = zeros(ny,1);
	 err = 1;
      else
	 if nr > 2 
	    if s(nr-1,nr-1) < mst
               disp(['Warning: Second smallest singular value is less than eep: ', num2str(epp),','])
               disp(['         S(',int2str(nr-1),',',int2str(nr-1),'): ', num2str(s(nr-1,nr-1)),'.'])
               err = -1;
	    end
	 end
         I(:,i) = up(1:nu);
         Opdsvd.m000064400136200001130000000061700605760564100130320ustar00kjetilentek00000000000000%   [I, O, err] = pdsvd(G, P, epp)
%
%   The ``Pole-Direction-SVD'' function computes the input and
%   output pole directions through the use of Singular Value 
%   Decomposition of G(p) where p is a pole. Note that this is an
%   ill-posed problem. Better way of calculating the pole directions 
%   are given in PDFST-function.
%   If largest singular value of G(p) is larger than EPP then the 
%   corrosponding input and output singular vectors are stored in
%   columns of I and O. Otherwise, zero vectors are stored in the 
%   corrosponding columns.
%   A warning message is printed if also the second singular 
%   value is grater than EPP.
%   
%   Inputs:  G   - system matrix.
%            P   - vector of poles.
%            epp - tolerance (optional) for singular values, 
%                  default value is 1/eps which most often is
%                  to large typical value around 1e14.
%
%   Outputs: I   - Matrix containing the input directions, 
%                  column i corresponds to pole P(i).
%            O   - Matrix containing the output directions, 
%                  column i corresponds to pole P(i).
%            err - Error indicator:
%                     err =  0, all is OK.
%                     err =  1, No direction with gain larger 
%                               than EPP.
%                     err = -1, One basis vector is not enough
%                               to describe the subspace with 
%                               gains larger than EPP.
%
%   Written by: Kjetil Havre, 24/11-1995, e-mail: kjetil@ife.no
%
%   See also: IZDE, OZDE, ZDSVD and PDFST.
%
%   NOTE! This function should not be used for anything else than 
%         checking the results returned by PDFST, since it depends 
%         poor numerical properties. 
%         In the function PDFST the pole directions are calculated 
%         avoiding these poor numerical properties.
%
   function [I, O, err] = pdsvd(G, P, epp)

   I = []; O = []; err = 0;

   [mt,ny,nu, nx] = minfo(G);
   
   if strcmp(mt, 'syst')==0
      disp( 'System matrix G is required, usage: [I,O] = pdsvd(G, P, epp)' );
      return
   end
   
   if nargin < 2
      disp( 'Usage: [I,O] = pdsvd(G, P, epp)' );
      return
   end
   
   if nargin < 3
      epp = 1/eps;
   end
   
   for i=1:max(size(P))
      wp  = P(i)/sqrt(-1);
      Gp = frsp(G, wp);              % This may cause problems.
      [u,s,v] = svd(Gp(1:ny,1:nu));
      nr = min(size(s));
      if s(1, 1) > epp
         if nr > 1
            if s(2, 2) > epp
               disp(['Warning: Second singular value is larger than eep: ', num2str(epp),','])
               disp(['         S(2,2): ', num2str(s(nr-1,nr-1)),'.'])
	       err = -1;
	    end
	 end
	 idx = find(abs(v(:,nr)) > epp);
	 vn = angle(v(idx(1),nr));
         I(:,i) = v(:,nr)*exp(-vn*sqrt(-1));
	 
	 idx = find(abs(u(:,nr)) > epp);
	 vn = angle(u(idx(1),nr));
	 O(:,i) = u(:,nr)*exp(-vn*sqrt(-1));
      else
         disp(['Warning: Gain of G(p) is not larger than EPP: ', num2str(epp)]);
	 disp(['         for i=',int2str(i),', p=',num2str(P(i)),'.'])
         I(:,i) = zeros(nu,1);
	 O(:,i) = zeros(ny,1);
         err = 1;
      end
   end
end

   
,nr-1)),'.'])
	       err = -1;
	    end
	 end
	 idx = find(abs(v(:,nr)) > epp);
	 vn = angle(v(idx(1),nr));
         I(:,i) = v(:,nr)*exp(-vn*sqrt(-1));
	 
	 idx = find(abs(u(:,nr)) > epp);
	 vn = angle(u(idx(1),nr));
	 O(:,i) = u(:,nr)*exp(-vn*sqrt(-1));
      else
         disp(['Warning: Gain of G(p) is not larger than EPP: ', num2str(epp)]);
	 disp(['         for i=',int2str(i),', p='poledir.m000064400136200001130000000006420575144150100133370ustar00kjetilentek00000000000000function [up,sp,vp,Gip] = CompPolDir( G, p );
   [A,B,C,D] = unpck( G );
   sd = svd( D );
   ms = max( max( sd ) );
   if ms == 0
      disp( 'Making system proper by adding a small term in D.' );
      Dm = [0.0001 0; 0 0.0001];
      G = pck( A, B, C, Dm );
   end
   Gi= minv( G );
%   Gi = sysbal( Gi );
   [Ai,Bi,Ci,Di] = unpck( Gi );
   Gip = Ci*inv(eye(size(Ai))*p-Ai)*Bi+Di;
   [up, sp, vp] = svd( Gip );
end
of G(p) is not larger than EPP: ', num2str(epp)]);
	 disp(['         for i=',int2str(i),', p='selc.m000064400136200001130000000007150602752167700126420ustar00kjetilentek00000000000000%
%  Mo = selc(M,r,c);
%
%  Selects rows and columns of M as described in r and c. This is similar 
%  to sel. However it works on constant matrices. 
%

function Mo = selc(M,r,c)
   if strcmp(r, ':') == 1
      Mr = M;
   else
      Mr = [];
      for i=1:max(size(r))
         Mr = [Mr; M(r(i),:)];
      end
   end
   if strcmp(c, ':') == 1
      Mo = Mr;
   else
      Mo = [];
      for i=1:max(size(c))
         Mo = [Mo Mr(:,c(i))];
      end
   end
end
epp)]);
	 disp(['         for i=',int2str(i),', p='sewdmu.m000064400136200001130000000047370605312111700132100ustar00kjetilentek00000000000000% function [delta,lowerbnd,upperbnd] = wcperf(sysg,blk,pertsize)
function [delta,lowerbnd,upperbnd] = wcperf(sysg,blk,pertsize,startalpha)

    if nargin == 0
        disp('usage: [delta,lowerbound,upperbound] = wcperf(sysg,blk,pertsize)')
        return
    elseif nargin == 3
        startalpha = [.3 1 1/.3];
    end

 perctol = 0.001;
 itermax = 10;
    

 [nblk,dum] = size(blk);
 [mtype,mrows,mcols,mnum] = minfo(sysg);
 repblks = find(blk(:,2)==0);
 ablk = abs(blk)
 if length(repblks)>0
    ablk(repblks,2) = ablk(repblks,1);
 end
 rdimblk = sum(ablk(:,1));
 cdimblk = sum(ablk(:,2));
 if rdimblk<mcols & cdimblk<mrows
    blk = [blk;[mcols-rdimblk mrows-cdimblk]];
    nblk = nblk + 1;
 elseif rdimblk~=mcols | cdimblk~=mrows
    error('dimension error in Uncertainty Structure')
    return
 end
 
 
 cperf = blk(nblk,1);
 rperf = blk(nblk,2);
 if blk(nblk,2) == 0
   rperf = blk(nblk,1);
 end

% FIRST POINT
 nominal = pkvnorm(sel(sysg,mrows-rperf+1:mrows,mcols-cperf+1:mcols));
% LAST POINT
 if any(blk(1:nblk-1,1)>0)
    [bnds,dvec,sens,pvec] = mu(sel(sysg,1:mrows-rperf,1:mcols-cperf),blk(1:nblk-1,:),'csw');
    xdestab = 1/pkvnorm(sel(bnds,1,2)); 
 else
    xdestab = inf;
    disp('Warning: Uncertainty Structure is All Real - There may be problems')
 end

    if pertsize>=xdestab
        disp('Perturbation Size is destabilizing')
        delta = dypert(pvec,blk(1:nblk-1,:),bnds,[1:nblk-1]);
        return
    end

    xl = 0;
    yl = nominal;
    xu = 0;
    yu = nominal;

    xmax = xdestab;;

    [out,delta] = onealpha(sysg,blk,startalpha,pertsize,perctol);
    xl = [xl;out(:,1)]
    yl = [yl;out(:,2)] 
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
 
 if isempty(delta)
    dist = inf;
    bestdelta = []; 
 else
    ndelta = hinfnorm(delta);
    ndelta = ndelta(1);
    bestdelta = delta;
    dist = abs(ndelta-pertsize);
 end
 cnt = 0;

 while cnt<itermax & isempty(delta)
    p = polyfit(xl,yl,min([4 length(xl)-2]));
    newalpha = 1/pertsize/polyval(p,pertsize);
    [out,delta] = onealpha(sysg,blk,newalpha,pertsize,perctol);
    if ~isempty(delta)
        ndelta = hinfnorm(delta);
        ndelta = ndelta(1);
        if abs(ndelta-pertsize)<dist
            bestdelta = delta;
            dist = abs(ndelta-pertsize);
        end
    end
    xl = [xl;out(:,1)];
    yl = [yl;out(:,2)]; 
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
    cnt = cnt+1;
end
delta = bestdelta;
[xls,index] = sort(xl);
lowerbnd = vpck(yl(index),xls);
[xus,index] = sort(xu);
upperbnd = vpck(yu(index),xus);
= onealpha(sysg,blk,newalpha,pertskewmuandy.m000064400136200001130000000031370605436116500140760ustar00kjetilentek00000000000000function [bnds,pert] = skewmu(matin,blk,fixnorm)

    [nblk,dum] = size(blk);
    blkp = ptrs(blk);
    if dum~=2
        error('Invalid Block structure');
        return
    end

    fixnorm = fixnorm(:);
    [nr,nc] = size(fixnorm);
    if nr~=nblk & nc~=1
        error('Invalid Fix-Norm Specification');
        return
    end
    
    [mtype,mrows,mcols,mnum] = minfo(matin);
    leftone = [];
    lefttwo = ones(mrows,1);
    rightone = [];
    righttwo = ones(mcols,1);
    for i=1:nblk
        ldim = abs(blk(i,1));
        if blk(i,2)==0
            rdim = ldim;
        else
            rdim = blk(i,2);
        end
        if isinf(fixnorm(i))
            leftone = [leftone blkp(i,2):blkp(i+1,2)-1];
            rightone = [rightone blkp(i,1):blkp(i+1,1)-1];
            lefttwo(blkp(i,2):blkp(i+1,2)-1) = ones(1,rdim);
            righttwo(blkp(i,1):blkp(i+1,1)-1) = ones(1,ldim);
        else
            lefttwo(blkp(i,2):blkp(i+1,2)-1) = sqrt(fixnorm(i))*ones(1,rdim);
            righttwo(blkp(i,1):blkp(i+1,1)-1) = sqrt(fixnorm(i))*ones(1,ldim);
        end    
    end        
    am = diag(lefttwo)*m*diag(righttwo);    
    while 
        sm = am;
        sm(leftone,:) = sqrt(alpha)*sm(leftone,:);
        sm(:,rightone) = sqrt(alpha)*sm(:,rightone);
        bb = mu(sm,blk);
        xdata = [xdata;alpha];
        ylow = [ylow;bb(2)];
        yupp = [yupp;bb(1)];        
        if bb(1)<1
            if alpha>alphau
                alphau = alpha;
            end
        end
        if bb(2)>=1
            if alpha<alphal
                alphal = alpha;
            end
        end
    end
            
m*diag(righttwo);    
    while 
        sm = am;
        sm(leftone,:) = sqrt(alpha)*sm(leftone,:);
        sm(:,rightone) = sqrt(alpha)*sm(:,rightone);
        bb = mu(sm,blk);
        xdata = [xdata;alpha];
        ylow = [ylow;bb(2)];
        yupp = [yupp;bb(1)];        
        if bb(1)<1
            if alpha>alphau
                alphau = alpha;
            end
        end
        if bb(2)>=1
            ifskewmuold.m000064400136200001130000000034450605437245500137270ustar00kjetilentek00000000000000function [bnds,pert] = skewmu(matin,blk,fixnorm)

    [nblk,dum] = size(blk);
    blkp = ptrs(blk);
    if dum~=2
        error('Invalid Block structure');
        return
    end

    fixnorm = fixnorm(:);
    [nr,nc] = size(fixnorm);
    if nr~=nblk & nc~=1
        error('Invalid Fix-Norm Specification');
        return
    end
    
    [mtype,mrows,mcols,mnum] = minfo(matin);
    leftone = [];
    lefttwo = ones(mrows,1);
    rightone = [];
    righttwo = ones(mcols,1);
    for i=1:nblk
        ldim = abs(blk(i,1));
        if blk(i,2)==0
            rdim = ldim;
        else
            rdim = blk(i,2);
        end
        if isinf(fixnorm(i))
            leftone = [leftone blkp(i,2):blkp(i+1,2)-1];
            rightone = [rightone blkp(i,1):blkp(i+1,1)-1];
            lefttwo(blkp(i,2):blkp(i+1,2)-1) = ones(1,rdim);
            righttwo(blkp(i,1):blkp(i+1,1)-1) = ones(1,ldim);
        else
            lefttwo(blkp(i,2):blkp(i+1,2)-1) = sqrt(fixnorm(i))*ones(1,rdim);
            righttwo(blkp(i,1):blkp(i+1,1)-1) = sqrt(fixnorm(i))*ones(1,ldim);
        end    
    end  

    leftone
    rightone
    
    lefttwo
    righttwo

    
    [iv] = getiv(matin);
    
    mr = 1;
    alpha = 1;
    m = xtract(matin, iv(1))
    diag(lefttwo)
    diag(righttwo)
    am = mmult(diag(lefttwo),m,diag(righttwo));    
    while mr
        sm = am;
        sm(leftone,:) = sqrt(alpha)*sm(leftone,:)
        sm(:,rightone) = sqrt(alpha)*sm(:,rightone)
        [bb, rd, sn, rp, rg] = mu(sm,blk)
        xdata = [xdata;alpha]
        ylow = [ylow;bb(2)]
        yupp = [yupp;bb(1)]      
        if bb(1)<1
            if alpha>alphau
                alphau = alpha;
            end
        end
        if bb(2)>=1
            if alpha<alphal
                alphal = alpha;
            end
        end
    end
            
wo));    
    while mr
        sm = am;
        sm(leftone,:) = sqrt(alpha)*sm(leftone,:)
        sm(:,rightone) = sqrt(alpha)*sm(:,rightone)
        [bb, rd, sn, rp, rg] = mu(sm,blk)
        xdata = [xdata;alpha]
     skrpmu.m000064400136200001130000000113040605461055600132240ustar00kjetilentek00000000000000function [skm, bnds,rowd,sens,rowp,rowg] = skrpmu(matin,blk,tol,opt,maxit)
%
%  [skm, bnds,rowd,sens,rowp,rowg] = skrpmu(matin,blk,tol,opt)
%
%  Function skrpmu computes skewed mu by scaling the part of the plant
%  corrosponding to the last perturbation block so that mu(Km*M) = 1,
%  where Km = diag(I(1), I(2), ... , I(nblk-1), km*I(np) )
%  I(i) denotes the identity of delta block i. 
%  It is assumed that diag(I(1), I(2), ..., I(nblk-1)) is square.
%  The last last delt block can be nonsquare.
%  The skrpmu function iterates on km until th LOWER bound from the 
%  mu(Km*M) calculation is equal to one within the specified tolerance tol.
%  Nessecary requirement is that mu(M11) < 1 for all frequencies with 
%  perturbation equal to the n-1 first blocks.
%
%  Inputs:
%     matin       varying matrix for the frequenct respons.
%     blk         description of perturbation block, same as for mu.
%                 The size of matin must be consistent with the specified
%                 block perturbation blk.
%     tol         tolerance for exiting iteration, default 1e-6.
%     opt         options to be passe on to mu-Tools function mu.
%     maxit       maximum number of iteration, default value 20.
%
%  Outputs:
%     skm = 1/km, skewd mu with last perturbation held constant.
%     bnds,       bnds from the mu calculation, both lower and upper 
%                 should be equal to one.
%     rowd,       rowd from last mu calculation.
%     sens,       sens from last mu calculation.
%     rowp,       rowp from last mu calculation.
%     rowg,       rowg from last mu calculation.
%
%  Note! skewed mu calculation by iterating on mu calculations in 
%  mu tools is very slow. However, it is easy to implement.
%  More efficient algorithms for calculation of skewed mu
%  are likely to be provided by mu tools in next version.
%  
%  Written by Kjetil Havre, e-mail: kjetil@ife.no after proposal  
%  by Sigurd Skogestad, e-mail: skoge@kjemi.unit.no.
%

    if nargin < 5
       maxit = 20;
    end    
    if nargin < 4
       opt = 'lu';
    end
    if nargin < 3
       tol = 1e-6;
    end
    if nargin < 2
       disp( '   Error in skrpmu: To few input arguments.')
       disp( '      [skm, bnds,rowd,sens,rowp,rowg] = skewedRPmu(matin,blk)' )
       return
    end
    

    [nblk,dum] = size(blk);
    blkp = ptrs(blk);
    if dum~=2
        error('Invalid Block structure');
        return
    end

    blkRP = blk(nblk, :);
    nI2 = min(blkRP);
    
    I2 = eye(nI2); 
    R2=abv(I2,zeros(blkRP(1)-nI2,nI2));
    R2=sbs(R2, zeros(blkRP(1), blkRP(2)-nI2));
    
    nr1 = sum(blk(1:nblk-1,1));
    R1 = eye(nr1);
    Z1 = zeros(nr1, blkRP(2));
    Z2 = zeros(blkRP(1), nr1);
    
    iv = getiv(matin);
    np = max(size(iv));
    ma  = xtract(matin, iv(1));
    [b, rd, s, rp, rg] = mu(ma,blk, opt); 
    fp = b(1) - 1;xp = 1; 
    x  = 1/b(1);
    
    skm = zeros(np+1, 2);
    [nr, nc] = size( b); bnds=zeros(np+1,nc);
    [nr, nc] = size(rd); rowd=zeros(np+1,nc);
    [nr, nc] = size( s); sens=zeros(np+1,nc);
    [nr, nc] = size(rp); rowp=zeros(np+1,nc);
    if isempty(rg)
       rowg = [];
    else
       [nr, nc] = size(rg); rowg=zeros(np+1,nc);
    end
    
    for i=1:np
       disp(['Frequency point ', int2str(i)] );
       if i > 1
          ma = xtract(matin, iv(i));
          Km = [R1 Z1; Z2 xp*R2];
          [b, rd, s, rp, rg] = mu(mmult(Km,ma),blk,opt); fp = b(1)-1;
       end
       
       Km = [R1 Z1; Z2 x*R2];
       [b, rd, s, rp, rg] = mu(mmult(Km,ma),blk,opt); 
       f = b(1) - 1;
       mr = 1; antit = 1;
       while (abs(f) > tol) & (mr==1)
          xn = x-f*(x-xp)/(f-fp); xn = abs(xn);
	  Km = [R1 Z1; Z2 xn*R2];
	  m = mmult(Km,ma);
          [b, rd, s, rp, rg] = mu(m,blk,opt);
	  fp = f; f = b(1)-1;
	  xp = x; x = xn;
	  if antit >= maxit
	     disp('Warning: Maximum number of iterations, exiting without convergence.')
             mr = 0;
	  end
	  antit = antit + 1;
       end
       disp(['Frequency point ', int2str(i),'. Solution for  km = ', num2str(x), ', residual = ', num2str(f)])
       skm(i,:)   = [1/x iv(i)];
       bnds(i,:) =  b(1,:);
       rowd(i,:) = rd(1,:);
       sens(i,:) =  s(1,:);
       rowp(i,:) = rp(1,:);
       if isempty(rowg) == 0
           rowg(i,:) = rg(1,:);
       end
    end
    skm(np+1,:) =  [np Inf];  
    [nr, nc] = size( b); bnds(np+1,:) = zeros(1, nc); bnds(np+1, nc-1:nc) = [np Inf];
    [nr, nc] = size(rd); rowd(np+1,:) = zeros(1, nc); rowd(np+1, nc-1:nc) = [np Inf];
    [nr, nc] = size( s); sens(np+1,:) = zeros(1, nc); sens(np+1, nc-1:nc) = [np Inf];
    [nr, nc] = size(rp); rowp(np+1,:) = zeros(1, nc); rowp(np+1, nc-1:nc) = [np Inf];
    if isempty(rowg) == 0
       [nr, nc] = size(rg); rowg(np+1,:) = zeros(1, nc); rowg(np+1, nc-1:nc) = [np Inf];
    end
        
 end
    skm(np+1,:) =  [np Inf];  
    [nr, nc] = size( b); bnds(np+1,:) = zeros(1, nc); bnds(np+1, nc-1:nc) = [np Inf];
    [nr, nc] = size(rd); rowd(np+1,:) = zeros(1, nc); rowd(np+1, nc-1:nc) = [np Inf];
    [nr, nc] = size( s); sens(np+1,:) = zeros(1, nc); sens(np+1, nc-1:nc) = [np Inf];
    [nr, nc] = size(rpsminreal.m000064400136200001130000000004740606404671400135230ustar00kjetilentek00000000000000%
%
%   Function SMINREAL is uses minreal to obtain a minimal realization.
%
%   See: MINREAL
%

function Gm = minreal(G, tol)

   if nargin < 2
      tol = 1e-12;
   end

   [A,B,C,D] = unpck(Gm);
   
   [Am,Bm,Cm,Dm] = minreal(A, B, C, D, tol);
   Gm = pck(Am,Bm,Cm,Dm);
   if isempty(Gm)
      Gm = G;
   end
end
;
    [nr, nc] = size(rd); rowd(np+1,:) = zeros(1, nc); rowd(np+1, nc-1:nc) = [np Inf];
    [nr, nc] = size( s); sens(np+1,:) = zeros(1, nc); sens(np+1, nc-1:nc) = [np Inf];
    [nr, nc] = size(rpvrga.m000064400136200001130000000003420604670010300126270ustar00kjetilentek00000000000000function rga = vrga(G)
   [mt,ny,nu,nx] = minfo(G); rga = []; 

   if strcmp( mt, 'syst' ) == 1 
       disp( 'Not implemented for system matrices.' );
       return;      
   end
   
   rga = veval('.*',G,vpinv(vtp(G)));
end
al(A, B, C, D, tol);
   Gm = pck(Am,Bm,Cm,Dm);
   if isempty(Gm)
      Gm = G;
   end
end
;
    [nr, nc] = size(rd); rowd(np+1,:) = zeros(1, nc); rowd(np+1, nc-1:nc) = [np Inf];
    [nr, nc] = size( s); sens(np+1,:) = zeros(1, nc); sens(np+1, nc-1:nc) = [np Inf];
    [nr, nc] = size(rpvsum.m000064400136200001130000000003050603054456500126740ustar00kjetilentek00000000000000function o = vsum(A)
   [mt,nr,nc,n] = minfo(A);
   if strcmp( mt, 'vary' ) ~= 1
       disp( 'Only implemented for varying matrices.' );
       return;      
   end
   
   o = veval('sum',A);
end
l('.*',G,vpinv(vtp(G)));
end
al(A, B, C, D, tol);
   Gm = pck(Am,Bm,Cm,Dm);
   if isempty(Gm)
      Gm = G;
   end
end
;
    [nr, nc] = size(rd); rowd(np+1,:) = zeros(1, nc); rowd(np+1, nc-1:nc) = [np Inf];
    [nr, nc] = size( s); sens(np+1,:) = zeros(1, nc); sens(np+1, nc-1:nc) = [np Inf];
    [nr, nc] = size(rpwcperf.m000064400136200001130000000150470605461071200131730ustar00kjetilentek00000000000000%function [delta,lowerbnd,upperbnd] = wcperf(Mg,blk,pertsizew,npts)
%
%   Calculates worst-case performance for an upper-loop
%       LinearFractionalTransformation.
%       
%   Inputs: Mg - input matrix, CONSTANT/VARYING
%       BLK   - block structure
%       OPT   - optional argument, character string
%
%   Outputs:    BNDS  - upper and lower bounds
%       SENS  - sensitivity of || D M D^{-1} || to D scaling
%       ROWP  - perturbation from lower bound
%       ROWD  - D scaling from upper bound
%       ROWG  - G scaling from upper bound
%           ROWD, SENS, ROWP, ROWG are in compressed form
%
%   See also:  MU, DYPERT

function [delta,lowerbnd,upperbnd] = wcperf(sysg,blk,pertsizew,npts,opt)

    diagnose = 0;

    if nargin == 0 | nargin == 1 | nargin == 2
        disp('usage: [delta,lowerbound,upperbound] = wcperf(sysg,blk,pertsize)')
        return
    elseif nargin == 3
        npts = 1;
        opt = 'd';
    elseif nargin == 4 | nargin == 5
        if isempty(npts)
            npts = 5;
        else
            if floor(npts(1))==ceil(npts(1)) & npts(1)>0
                npts = npts(1);
            else
                error('Npts should be a positive integer');
                return
            end
        end
    end
    if nargin == 5
        if ~isstr(opt)
            error('Options should be a string');
            return
        end
    elseif nargin == 4
        opt = 'd';  % default
    end

    startalpha = [.3 1 1/.3];
    perctol = 0.01;
    itermax = npts + 10;

    [nblk,dum] = size(blk);
    if dum~=2
        error('Invalid Block structure');
        return
    end
    [mtype,mrows,mcols,mnum] = minfo(sysg);
    repblks = find(blk(:,2)==0);
    ablk = abs(blk);
    if length(repblks)>0
        ablk(repblks,2) = ablk(repblks,1);
    end
    rdimblk = sum(ablk(:,1));
    cdimblk = sum(ablk(:,2));
    if rdimblk<mcols & cdimblk<mrows
        blk = [blk;[mcols-rdimblk mrows-cdimblk]];
        nblk = nblk + 1;
        ablk = abs(blk);
    elseif rdimblk~=mcols | cdimblk~=mrows
        error('Dimension error in Uncertainty Structure')
        return
    end
    cperf = ablk(nblk,1);
    rperf = ablk(nblk,2);
    if rperf == 0
        disp('Performance block is Repeated Scalar...');
        rperf = cperf;      % get dimension correct
    end

    nominal = pkvnorm(sel(sysg,mrows-rperf+1:mrows,mcols-cperf+1:mcols));
    if any(blk(1:nblk-1,1)>0)
        [bnds,dvec,sens,pvec] = mu(sel(sysg,1:mrows-rperf,1:mcols-cperf),...
            blk(1:nblk-1,:),'csw');
        mubnd = pkvnorm(sel(bnds,1,2));
        if mubnd==0
            disp('        Guaranteed: No Perturbation Destabilizes');
            xdestab = inf;
            topbnd = 2*pertsizew;
        else
            xdestab = 1/pkvnorm(sel(bnds,1,2)); 
            topbnd = xdestab;
        end
    else
        [bnds,dvec,sens,pvec] = mu(sel(sysg,1:mrows-rperf,1:mcols-cperf),...
            blk(1:nblk-1,:),'csw');
        xdestabi = pkvnorm(sel(bnds,1,2)); 
        if xdestabi==0
            disp('Warning: Uncertainty Structure is All Real');
            [bnds,dvec,sens,pvec] = mu(sel(sysg,1:mrows-rperf,1:mcols-cperf),...
                abs(blk(1:nblk-1,:)),'csw');
            cmubnd = pkvnorm(sel(bnds,1,2)); 
            if cmubnd==0
                disp('        Guaranteed: No Perturbation Destabilizes');
                xdestab = inf;
                topbnd = 2*pertsizew;
            else
                disp('        Unable to find destabilizing Delta');
                xdestab = inf;
                topbnd = 1/cmubnd;
            end
        else
            xdestab = 1/xdestabi;
            topbnd = xdestab;
        end
    end

    deltafound = 0;
    if pertsizew>=xdestab
        disp('Perturbation Size is destabilizing')
        deltadestab = dypert(pvec,blk(1:nblk-1,:),bnds,[1:nblk-1]);
        pertsizew = xdestab;
        deltafound = 1;
    end

    xl = 0;
    yl = nominal;
    xu = 0;
    yu = nominal;


    [out,delta] = onealpha(sysg,blk,startalpha,pertsizew,perctol);
    xl = [xl;out(:,1)];
    yl = [yl;out(:,2)];
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
    [xls,index] = sort(xl);
    yls = yl(index); 
    if isempty(delta)
        dist = inf;
    else
        ndelta = hinfnorm(delta);
        ndelta = ndelta(1);
        bestdelta = delta;
        dist = abs(ndelta-pertsizew);
    end
    cnt = 0;
    gap = max(diff(xls));
    gap = max([gap topbnd-max(xls)]);

 while cnt<itermax & (gap>topbnd/npts | deltafound==0) & nargin >= 3
    newalpha = [];
    if length(xl)==2
        p = polyfit(xl,yl,1);
    elseif length(xl)==3
        p = polyfit(xl,yl,1);
    elseif length(xl)==4
        p = polyfit(xl,yl,2);
    else
        p = polyfit(xl,yl,min([5 length(xl)-3]));
    end
    if diagnose==1
        xx = linspace(0,1.2*topbnd,100);
        plot(xl,yl,'+',xx,polyval(p,xx));
        drawnow
    end
    f1=0;f2=0;
    if gap>topbnd/npts
        f1 = 1;
        loc = find(gap==diff(xls));
        if isempty(loc)
            pertsize = 1.2*max(xls);
            pdt = polyval(p,pertsize);
            alpha = 1/pertsize/pdt;
            newalpha = [newalpha alpha];
        else
            pertsize = (xls(loc)+xls(loc+1))/2;
            pdt = polyval(p,pertsize);
            fb = yls(loc);
            ft=fb+(yls(loc+1)-yls(loc))/(xls(loc+1)-xls(loc))*(pertsize-xls(loc));
            if pdt<fb | pdt>ft
                pdt = 0.5*(ft+fb);
            end
            alpha = 1/pertsize/pdt;
            newalpha = [newalpha alpha];
        end
    end
    if deltafound==0
        xli = max(find(xls<pertsizew));
        xleft = xls(xli);
        xri = min(find(xls>pertsizew));
        xright = xls(xri);
        pdt = polyval(p,pertsizew);
        fb = yls(xli);
        ft=fb+(yls(xri)-yls(xli))/(xls(xri)-xls(xli))*(pertsizew-xls(xli));
        if pdt<fb | pdt>ft
            pdt = 0.5*(ft+fb);
        end
        f2 = 1;
        newalpha = [newalpha 1/pertsizew/pdt];
    end
    [newalpha f1 f2];
    [out,delta] = onealpha(sysg,blk,newalpha,pertsizew,perctol);
    if ~isempty(delta)
        deltafound = 1;
        ndelta = hinfnorm(delta);
        ndelta = ndelta(1);
        if abs(ndelta-pertsizew)<dist
            bestdelta = delta;
            dist = abs(ndelta-pertsizew);
        end
    end
    xl = [xl;out(:,1)];
    yl = [yl;out(:,2)]; 
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
    cnt = cnt+1;
    [xls,index] = sort(xl);
    gap = max(diff(xls));
    gap = max([gap topbnd-max(xls)]);
    yls = yl(index);
end
delta = bestdelta;
if nargout > 1
    [xls,index] = sort(xl);
    lowerbnd = vpck(yl(index),xls);
    [xus,index] = sort(xu);
    upperbnd = vpck(yu(index),xus);
end
elta-pertsizew)<dist
            bestdelta = delta;
            dist = abs(ndelta-pertsizew);
        end
    end
    xl = [xl;out(:,1)];
    yl = [yl;out(:,2)]; 
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
    cnt = cnt+1;
    [xls,index] = sort(xl);
    gap = max(diff(xls));
    gap = max([gap topbnd-max(xls)]);
    yls = yl(index);
end
delta = bestdelta;
if nargout > 1
    [xls,index] = sort(xl);
    lowerbnd = vpck(yl(index),xls);
    [xus,index] = sort(xu);
 wcperfV2.m000064400136200001130000000047370605015426100134050ustar00kjetilentek00000000000000% function [delta,lowerbnd,upperbnd] = wcperf(sysg,blk,pertsize)
function [delta,lowerbnd,upperbnd] = wcperf(sysg,blk,pertsize,startalpha)

    if nargin == 0
        disp('usage: [delta,lowerbound,upperbound] = wcperf(sysg,blk,pertsize)')
        return
    elseif nargin == 3
        startalpha = [.3 1 1/.3];
    end

 perctol = 0.001;
 itermax = 10;
    

 [nblk,dum] = size(blk);
 [mtype,mrows,mcols,mnum] = minfo(sysg);
 repblks = find(blk(:,2)==0);
 ablk = abs(blk)
 if length(repblks)>0
    ablk(repblks,2) = ablk(repblks,1);
 end
 rdimblk = sum(ablk(:,1));
 cdimblk = sum(ablk(:,2));
 if rdimblk<mcols & cdimblk<mrows
    blk = [blk;[mcols-rdimblk mrows-cdimblk]];
    nblk = nblk + 1;
 elseif rdimblk~=mcols | cdimblk~=mrows
    error('dimension error in Uncertainty Structure')
    return
 end
 
 
 cperf = blk(nblk,1);
 rperf = blk(nblk,2);
 if blk(nblk,2) == 0
   rperf = blk(nblk,1);
 end

% FIRST POINT
 nominal = pkvnorm(sel(sysg,mrows-rperf+1:mrows,mcols-cperf+1:mcols));
% LAST POINT
 if any(blk(1:nblk-1,1)>0)
    [bnds,dvec,sens,pvec] = mu(sel(sysg,1:mrows-rperf,1:mcols-cperf),blk(1:nblk-1,:),'csw');
    xdestab = 1/pkvnorm(sel(bnds,1,2)); 
 else
    xdestab = inf;
    disp('Warning: Uncertainty Structure is All Real - There may be problems')
 end

    if pertsize>=xdestab
        disp('Perturbation Size is destabilizing')
        delta = dypert(pvec,blk(1:nblk-1,:),bnds,[1:nblk-1]);
        return
    end

    xl = 0;
    yl = nominal;
    xu = 0;
    yu = nominal;

    xmax = xdestab;;

    [out,delta] = onealpha(sysg,blk,startalpha,pertsize,perctol);
    xl = [xl;out(:,1)]
    yl = [yl;out(:,2)] 
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
 
 if isempty(delta)
    dist = inf;
    bestdelta = []; 
 else
    ndelta = hinfnorm(delta);
    ndelta = ndelta(1);
    bestdelta = delta;
    dist = abs(ndelta-pertsize);
 end
 cnt = 0;

 while cnt<itermax & isempty(delta)
    p = polyfit(xl,yl,min([4 length(xl)-2]));
    newalpha = 1/pertsize/polyval(p,pertsize);
    [out,delta] = onealpha(sysg,blk,newalpha,pertsize,perctol);
    if ~isempty(delta)
        ndelta = hinfnorm(delta);
        ndelta = ndelta(1);
        if abs(ndelta-pertsize)<dist
            bestdelta = delta;
            dist = abs(ndelta-pertsize);
        end
    end
    xl = [xl;out(:,1)];
    yl = [yl;out(:,2)]; 
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
    cnt = cnt+1;
end
delta = bestdelta;
[xls,index] = sort(xl);
lowerbnd = vpck(yl(index),xls);
[xus,index] = sort(xu);
upperbnd = vpck(yu(index),xus);
= onealpha(sysg,blk,newalpha,pertwcperfV3b.m000064400136200001130000000150470605411661500135510ustar00kjetilentek00000000000000%function [delta,lowerbnd,upperbnd] = wcperf(Mg,blk,pertsizew,npts)
%
%   Calculates worst-case performance for an upper-loop
%       LinearFractionalTransformation.
%       
%   Inputs: Mg - input matrix, CONSTANT/VARYING
%       BLK   - block structure
%       OPT   - optional argument, character string
%
%   Outputs:    BNDS  - upper and lower bounds
%       SENS  - sensitivity of || D M D^{-1} || to D scaling
%       ROWP  - perturbation from lower bound
%       ROWD  - D scaling from upper bound
%       ROWG  - G scaling from upper bound
%           ROWD, SENS, ROWP, ROWG are in compressed form
%
%   See also:  MU, DYPERT

function [delta,lowerbnd,upperbnd] = wcperf(sysg,blk,pertsizew,npts,opt)

    diagnose = 0;

    if nargin == 0 | nargin == 1 | nargin == 2
        disp('usage: [delta,lowerbound,upperbound] = wcperf(sysg,blk,pertsize)')
        return
    elseif nargin == 3
        npts = 1;
        opt = 'd';
    elseif nargin == 4 | nargin == 5
        if isempty(npts)
            npts = 5;
        else
            if floor(npts(1))==ceil(npts(1)) & npts(1)>0
                npts = npts(1);
            else
                error('Npts should be a positive integer');
                return
            end
        end
    end
    if nargin == 5
        if ~isstr(opt)
            error('Options should be a string');
            return
        end
    elseif nargin == 4
        opt = 'd';  % default
    end

    startalpha = [.3 1 1/.3];
    perctol = 0.01;
    itermax = npts + 10;

    [nblk,dum] = size(blk);
    if dum~=2
        error('Invalid Block structure');
        return
    end
    [mtype,mrows,mcols,mnum] = minfo(sysg);
    repblks = find(blk(:,2)==0);
    ablk = abs(blk);
    if length(repblks)>0
        ablk(repblks,2) = ablk(repblks,1);
    end
    rdimblk = sum(ablk(:,1));
    cdimblk = sum(ablk(:,2));
    if rdimblk<mcols & cdimblk<mrows
        blk = [blk;[mcols-rdimblk mrows-cdimblk]];
        nblk = nblk + 1;
        ablk = abs(blk);
    elseif rdimblk~=mcols | cdimblk~=mrows
        error('Dimension error in Uncertainty Structure')
        return
    end
    cperf = ablk(nblk,1);
    rperf = ablk(nblk,2);
    if rperf == 0
        disp('Performance block is Repeated Scalar...');
        rperf = cperf;      % get dimension correct
    end

    nominal = pkvnorm(sel(sysg,mrows-rperf+1:mrows,mcols-cperf+1:mcols));
    if any(blk(1:nblk-1,1)>0)
        [bnds,dvec,sens,pvec] = mu(sel(sysg,1:mrows-rperf,1:mcols-cperf),...
            blk(1:nblk-1,:),'csw');
        mubnd = pkvnorm(sel(bnds,1,2));
        if mubnd==0
            disp('        Guaranteed: No Perturbation Destabilizes');
            xdestab = inf;
            topbnd = 2*pertsizew;
        else
            xdestab = 1/pkvnorm(sel(bnds,1,2)); 
            topbnd = xdestab;
        end
    else
        [bnds,dvec,sens,pvec] = mu(sel(sysg,1:mrows-rperf,1:mcols-cperf),...
            blk(1:nblk-1,:),'csw');
        xdestabi = pkvnorm(sel(bnds,1,2)); 
        if xdestabi==0
            disp('Warning: Uncertainty Structure is All Real');
            [bnds,dvec,sens,pvec] = mu(sel(sysg,1:mrows-rperf,1:mcols-cperf),...
                abs(blk(1:nblk-1,:)),'csw');
            cmubnd = pkvnorm(sel(bnds,1,2)); 
            if cmubnd==0
                disp('        Guaranteed: No Perturbation Destabilizes');
                xdestab = inf;
                topbnd = 2*pertsizew;
            else
                disp('        Unable to find destabilizing Delta');
                xdestab = inf;
                topbnd = 1/cmubnd;
            end
        else
            xdestab = 1/xdestabi;
            topbnd = xdestab;
        end
    end

    deltafound = 0;
    if pertsizew>=xdestab
        disp('Perturbation Size is destabilizing')
        deltadestab = dypert(pvec,blk(1:nblk-1,:),bnds,[1:nblk-1]);
        pertsizew = xdestab;
        deltafound = 1;
    end

    xl = 0;
    yl = nominal;
    xu = 0;
    yu = nominal;


    [out,delta] = onealpha(sysg,blk,startalpha,pertsizew,perctol);
    xl = [xl;out(:,1)];
    yl = [yl;out(:,2)];
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
    [xls,index] = sort(xl);
    yls = yl(index); 
    if isempty(delta)
        dist = inf;
    else
        ndelta = hinfnorm(delta);
        ndelta = ndelta(1);
        bestdelta = delta;
        dist = abs(ndelta-pertsizew);
    end
    cnt = 0;
    gap = max(diff(xls));
    gap = max([gap topbnd-max(xls)]);

 while cnt<itermax & (gap>topbnd/npts | deltafound==0) & nargin >= 3
    newalpha = [];
    if length(xl)==2
        p = polyfit(xl,yl,1);
    elseif length(xl)==3
        p = polyfit(xl,yl,1);
    elseif length(xl)==4
        p = polyfit(xl,yl,2);
    else
        p = polyfit(xl,yl,min([5 length(xl)-3]));
    end
    if diagnose==1
        xx = linspace(0,1.2*topbnd,100);
        plot(xl,yl,'+',xx,polyval(p,xx));
        drawnow
    end
    f1=0;f2=0;
    if gap>topbnd/npts
        f1 = 1;
        loc = find(gap==diff(xls));
        if isempty(loc)
            pertsize = 1.2*max(xls);
            pdt = polyval(p,pertsize);
            alpha = 1/pertsize/pdt;
            newalpha = [newalpha alpha];
        else
            pertsize = (xls(loc)+xls(loc+1))/2;
            pdt = polyval(p,pertsize);
            fb = yls(loc);
            ft=fb+(yls(loc+1)-yls(loc))/(xls(loc+1)-xls(loc))*(pertsize-xls(loc));
            if pdt<fb | pdt>ft
                pdt = 0.5*(ft+fb);
            end
            alpha = 1/pertsize/pdt;
            newalpha = [newalpha alpha];
        end
    end
    if deltafound==0
        xli = max(find(xls<pertsizew));
        xleft = xls(xli);
        xri = min(find(xls>pertsizew));
        xright = xls(xri);
        pdt = polyval(p,pertsizew);
        fb = yls(xli);
        ft=fb+(yls(xri)-yls(xli))/(xls(xri)-xls(xli))*(pertsizew-xls(xli));
        if pdt<fb | pdt>ft
            pdt = 0.5*(ft+fb);
        end
        f2 = 1;
        newalpha = [newalpha 1/pertsizew/pdt];
    end
    [newalpha f1 f2];
    [out,delta] = onealpha(sysg,blk,newalpha,pertsizew,perctol);
    if ~isempty(delta)
        deltafound = 1;
        ndelta = hinfnorm(delta);
        ndelta = ndelta(1);
        if abs(ndelta-pertsizew)<dist
            bestdelta = delta;
            dist = abs(ndelta-pertsizew);
        end
    end
    xl = [xl;out(:,1)];
    yl = [yl;out(:,2)]; 
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
    cnt = cnt+1;
    [xls,index] = sort(xl);
    gap = max(diff(xls));
    gap = max([gap topbnd-max(xls)]);
    yls = yl(index);
end
delta = bestdelta;
if nargout > 1
    [xls,index] = sort(xl);
    lowerbnd = vpck(yl(index),xls);
    [xus,index] = sort(xu);
    upperbnd = vpck(yu(index),xus);
end
elta-pertsizew)<dist
            bestdelta = delta;
            dist = abs(ndelta-pertsizew);
        end
    end
    xl = [xl;out(:,1)];
    yl = [yl;out(:,2)]; 
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
    cnt = cnt+1;
    [xls,index] = sort(xl);
    gap = max(diff(xls));
    gap = max([gap topbnd-max(xls)]);
    yls = yl(index);
end
delta = bestdelta;
if nargout > 1
    [xls,index] = sort(xl);
    lowerbnd = vpck(yl(index),xls);
    [xus,index] = sort(xu);
 zdsvd.m000064400136200001130000000053030605760540000130320ustar00kjetilentek00000000000000%   [I, O, err] = zdsvd(G, Z, epp)
%
%   The ``Zero-Direction-SVD'' function computes the input and
%   output zero directions through the use of Singular Value 
%   Decomposition of G(z) where z is a zero.
%   If singular value number nr (nr=min{ny,nu}, ny-number of outputs, 
%   nu-number of inputs) of G(z) is less than EPP then the 
%   corrosponding input and output singular vectors are stored in
%   columns of I and O. Otherwise, zero vectors are stored in the 
%   corrosponding columns.
%   A warning message is printed if also the second smallest singular 
%   value is less than EPP.
%   
%   Inputs:  G   - system matrix.
%            Z   - vector of zeros.
%            epp - tolerance (optional) for singular values, 
%                  default value is 1e-14.
%
%   Outputs: I   - Matrix containing the input directions, 
%                  column i corresponds to zeros Z(i).
%            O   - Matrix containing the output directions, 
%                  column i corresponds to zeros Z(i).
%            err - Error indicator:
%                     err =  0, all is OK.
%                     err =  1, No direction with gain smaller 
%                               than EPP.
%                     err = -1, One basis vector is not enough
%                               to describe the subspace with 
%                               gains smaler than EPP.
%
%   Written by: Kjetil Havre, 24/11-1995, e-mail: kjetil@ife.no
%
%   See also: izde, ozde, pdsvd.
%
   function [I, O, err] = zdsvd(G, Z, epp)

   I = []; O = []; err = 0;
   [mt,ny,nu, nx] = minfo(G);
   
   if strcmp(mt, 'syst')==0
      disp( 'System matrix G is required, usage: [I,O] = zdsvd(G, Z, epp)' );
      return
   end
   
   if nargin < 2
      disp( 'Usage: [I,O] = zdsvd(G, Z, epp)' );
      return
   end
   
   if nargin < 3
      epp = 1e-14;
   end

   for i=1:max(size(Z))
      wz  = Z(i)/sqrt(-1);
      Gz = frsp(G, wz);
      [u,s,v] = svd(Gz(1:ny,1:nu));
      nr = min(size(s));
      if s(nr, nr) < epp
         if nr > 1
            if s(nr-1,nr-1) < epp
               disp(['Warning: Second smallest singular value is less than eep: ', num2str(epp),','])
               disp(['         S(',int2str(nr-1),',',int2str(nr-1),'): ', num2str(s(nr-1,nr-1)),'.'])
               err = -1;
	    end
	 end
	 idx = find(abs(v(:,nr)) > epp);
	 vn = angle(v(idx(1),nr));
         I(:,i) = v(:,nr)*exp(-vn*sqrt(-1));
	 
	 idx = find(abs(u(:,nr)) > epp);       % Rotate elements of vector.
	 vn = angle(u(idx(1),nr));
	 O(:,i) = u(:,nr)*exp(-vn*sqrt(-1));
      else
         disp(['Warning: G(z) is not singular for Z(i): ',num2str(Z(i)), ' for i=',int2str(i),'.'])
         I(:,i) = zeros(nu,1);
	 O(:,i) = zeros(ny,1);
	 err = 1;
      end
   end
end

   
str(s(nr-1,nr-1)),'.'])
               err = -1;
	    end
	 end
	 idx = find(abs(v(:,nr)) > epp);
	 vn = angle(v(idx(1),nr));
         I(:,i) = v(:,nr)*exp(-vn*sqrt(-1));
	 
	 idx = find(abs(u(:,nr)) > epp);       % Rotate elements of vector.
	 vn = angle(u(idx(1),nr));
	 O(:,i) = u(:,nr)*exp(-vn*sqrt(-1));
      elzerodir.m000064400136200001130000000057660575607370600134100ustar00kjetilentek00000000000000%  [u, s, v, Gz] = zerodir( G, z, epp );
%           or
%  [u, s, v, Gz] = zerodir( G, z );
%
%  This function computes the singular value decomposition of a given
%  system at s = z. The outputs are the output directions U, the singular
%  values S along the diagonal, the input directions V and the transfer 
%  function G(z).
%
%  The last ny-r columns of u, corresponds to the singular output directions.
%  ny - number of outputs, r - is the rank of G(z), if z is an zero then r<ny.
%  This should is a safe way to compute the directions of a zero. However, 
%  the zeros must be known in advance and if more than one singular value is 
%  zero or close to zero it can be difficult to pick out the corrosponding
%  direction.
%
%  If G(s) has a zero on the same location as a pole, then the matrix (zI-A)
%  becomes singular (This is also true when the directions of the pole and 
%  the zero are not the same). The function then adds small singular values 
%  to the zero singular values of (zI-A). The singular values of (zI-A)
%  which are zero are set to the EPP times the smallest non-zero singular 
%  value. If the parameter EPP is not supplied in the call to the function
%  the default value is 1e-8.
%
%  The method with generalized eigenvalues does not suffer from this weakness.
%
%  See also:
%              OZDE, IZDE ans SZEROS.
%
%  Written by: Kjetil Havre 20/4 - 1995.
%  e-mail:     kjetil@ife.no
%

%
%  Malfunctions: When the system has a zero at the same location of a pole,
%                the SVD function promts warning message about singular
%                matrix for G(z). At the moment I can not see why.
%                Maybe it is the inv function. It can be better to use the 
%                pseudo inverse. I think so.
%  The reason is that the matrix (zI-A) becomes singular. This is corrected 
%  by adding small singular values for the zero singular values. See 
%  description above.
%
%
%
%   4/5-1995     I tested the pseudo inverse but then the invers if G(z)
%                does not seems to agree with the zero direction computed 
%                eigenvalue decomposition. This malfunction only occures 
%                when z = p = 2. KH.
%
%

function [u, s, v, Gz] = zerodir( G, z, epp )

   if nargin < 3
      epp = 1e-8;
   end


   [A,B,C,D] = unpck( G );
   
   % May expect problems if Iz - A is singular. Have to look into this.
   H = z*eye(size(A)) - A;
   [U, S, V]  = svd( H );
   Si = zeros(size(S));
   nr = rank(S);
   [nxr, nx] = size( A );
   if nr < max(size(A))
      disp( ['Warning: Matrix (zI-A) is singular adds small singular values' ] );
      disp( ['         of size: ', num2str(S(nr,nr)*epp), ' to singular part.' ] );
      disp( ['         Result may be inaccurate, the result should be compared to'] );
      disp( ['         the result obtained with generalized eigenvalue methode.'] );
      for i=nr+1:nx,
         S(i,i) = S(nr,nr)*epp;
      end
   end
   for i=1:nx
      Si(i,i) = 1/S(i,i);
   end
   Hi = V*Si*U';
   Gz = C*Hi*B+D;
   [u, s, v] = svd( Gz );
end



 ['Warning  newalpha = [newalpha alpha];
        else
            pertsize = (xls(loc)+xls(loc+1))/2;
            pdt = polyval(p,pertsize);
            fb = yls(loc);
            ft=fb+(yls(loc+1)-yls(loc))/(xls(loc+1)-xls(loc))*(pertsize-xls(loc));
            if pdt<fb | pdt>ft
                pdt = 0.5*(ft+fb);
            end
            alpha = 1/pertsize/pdt;
            newalpha = [newalpha alpha];
        end
    end
    if deltafound==0
        xli = max(find(xls<pertsizew));
        xleft = xls(xli);
        xri = min(find(xls>pertsizew));
        xright = xls(xri);
        pdt = polyval(p,pertsizew);
        fb = yls(xli);
        ft=fb+(yls(xri)-yls(xli))/(xls(xri)-xls(xli))*(pertsizew-xls(xli));
        if pdt<fb | pdt>ft
            pdt = 0.5*(ft+fb);
        end
        f2 = 1;
        newalpha = [newalpha 1/pertsizew/pdt];
    end
    [newalpha f1 f2];
    [out,delta] = onealpha(sysg,blk,newalpha,pertsizew,perctol);
    if ~isempty(delta)
        deltafound = 1;
        ndelta = hinfnorm(delta);
        ndelta = ndelta(1);
        if abs(ndelta-pertsizew)<dist
            bestdelta = delta;
            dist = abs(ndelta-pertsizew);
        end
    end
    xl = [xl;out(:,1)];
    yl = [yl;out(:,2)]; 
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
    cnt = cnt+1;
    [xls,index] = sort(xl);
    gap = max(diff(xls));
    gap = max([gap topbnd-max(xls)]);
    yls = yl(index);
end
delta = bestdelta;
if nargout > 1
    [xls,index] = sort(xl);
    lowerbnd = vpck(yl(index),xls);
    [xus,index] = sort(xu);
    upperbnd = vpck(yu(index),xus);
end
elta-pertsizew)<dist
            bestdelta = delta;
            dist = abs(ndelta-pertsizew);
        end
    end
    xl = [xl;out(:,1)];
    yl = [yl;out(:,2)]; 
    xu = [xu;out(:,3)];
    yu = [yu;out(:,4)]; 
    cnt = cnt+1;
    [xls,index] = sort(xl);
    gap = max(diff(xls));
    gap = max([gap topbnd-max(xls)]);
    yls = yl(index);
end
delta = bestdelta;
if nargout > 1
    [xls,index] = sort(xl);
    lowerbnd = vpck(yl(index),xls);
    [xus,index] = sort(xu);
 zdsvd.m000064400136200001130000000053030605760540000130320ustar00kjetilentek00000000000000%   [I, O, err] = zdsvd(G, Z, epp)
%
%   The ``Zero-Direction-SVD'' function computes the input and
%   output zero directions through the use of Singular Value 
%   Decomposition of G(z) where z is a zero.
%   If singular value number nr (nr=min{ny,nu}, ny-number of outputs, 
%   nu-number of inputs) of G(z) is less than EPP then the 
%   corrosponding input and output singular vectors are stored in
%   columns of I and O. Otherwise, zero vectors are stored in the 
%   corrosponding columns.
%   A warning message is printed if also the second smallest singular 
%   value is less than EPP.
%   
%   Inputs:  G   - system matrix.
%            Z   - vector of zeros.
%            epp - tolerance (optional) for singular values, 
%                  default value is 1e-14.
%
%   Outputs: I   - Matrix containing the input directions, 
%                  column i corresponds to zeros Z(i).
%            O   - Matrix containing the output directions, 
%                  column i corresponds to zeros Z(i).
%            err - Error indicator:
%                     err =  0, all is OK.
%                     err =  1, No direction with gain smaller 
%                               than EPP.
%                     err = -1, One basis vector is not enough
%                               to describe the subspace with 
%                               gains smaler than EPP.
%
%   Written by: Kjetil Havre, 24/11-1995, e-mail: kjetil@ife.no
%
%   See also: izde, ozde, pdsvd.
%
   function [I, O, err] = zdsvd(G, Z, epp)

   I = []; O = []; err = 0;
   [mt,ny,nu, nx] = minfo(G);
   
   if strcmp(mt, 'syst')==0
      disp( 'System matrix G is required, usage: [I,O] = zdsvd(G, Z, epp)' );
      return
   end
   
   if nargin < 2
      disp( 'Usage: [I,O] = zdsvd(G, Z, epp)' );
      return
   end
   
   if nargin < 3
      epp = 1e-14;
   end

   for i=1:max(size(Z))
      wz  = Z(i)/sqrt(-1);
      Gz = frsp(G, wz);
      [u,s,v] = svd(Gz(1:ny,1:nu));
      nr = min(size(s));
      if s(nr, nr) < epp
         if nr > 1
            if s(nr-1,nr-1) < epp
               disp(['Warning: Second smallest singular value is less than eep: ', num2str(epp),','])
               disp(['         S(',int2str(nr-1),',',int2str(nr-1),'): ', num2str(s(nr-1,nr-1)),'.'])
               err = -1;
	    end
	 end
	 idx = find(abs(v(:,nr)) > epp);
	 vn = angle(v(idx(1),nr));
         I(:,i) = v(:,nr)*exp(-vn*sqrt(-1));
	 
	 idx = find(abs(u(:,nr)) > epp);       % Rotate elements of vector.
	 vn = angle(u(idx(1),nr));
	 O(:,i) = u(:,nr)*exp(-vn*sqrt(-1));
      el