ipde.m

国外经典书籍MULTIVARIABLE FEEDBACK CONTROL-多变量反馈控制 的源码

%   [Pi, Up, Xpi, Spi] = ipde(G,EPP,RF)
%
%   The ``Input-Pole-Direction-Eigenvalue'' ipde-function computes the 
%        output pole directions through the use of  egenvalue decomposition.
%   NOTE:  Avoids the change in ordering of eigenvalues that may occur if 
%        the input and output directions are computed separately.
%
%   Outputs: Pi   - Vector containing the poles.
%            Up   - Matrix containing the iinput pole vectors/directions (normalized = length 1), 
%                  column i corresponds to pole Pi(i).
%            Xpi   - Input state directions  = left eigenvectors of A (non-normalized).
%                    NOTE: These are scaled such that pole vector = pole directions
%            Spi   - Scalars used to scale the state directions (eigenvectors).
%
%   Inputs:  G   - System matrix on mutools format: G=pck(A,B,C,D)
%            EPP - (optional) Tolerance, default value is 1e-12.
%                  It is used to check the norm of the input directions.
%                  If the norm of input direction is less than EPP,
%                  zero is stored in the corresponding input direction.
%            RF  - (optional) If FLG_RF==1 then the first element of 
%                  the input direction is made real.
%
%   Written by: Kjetil Havre, 12/9-1996, e-mail: [e-mail]
%   Requires subroutines: zerm.m and sorte.m
%
%   See also:  OPDE, PDFST, IZDE, OZDE, ZDSVD and PDSVD.
%
%   Reference: Havre K. and S. Skogestad, 1996, ``Effect of RHP Zeros and
%              Poles on Performance in Multivariable Systems''.
%
   function [P,U,X,S] = ipde(G,epp,rf);

   P = []; U = []; X = []; 
   [mt,ny,nu, nx] = minfo(G);
   S = ones(nx,1);
   
   if strcmp(mt, 'syst')==0
      disp( 'System matrix G is required, usage: [] = ipde(G,epp)' );
      return
   end

   if nargin < 3
      rf = 0;
   end
   if nargin < 2
      epp = 1e-12;
   end
   [A, B, C, D] = unpck(G);
   
   if norm(D) > 1/epp
      disp('Warning: Directions may be inaccurate due to large effect from D.');
   end
   
   [V, D] = eig(A.','nobalance');   
   P = diag(D); X = V;
   [P,Is] = sorte(P); X = X(:,Is);
   Hm = A.'*X-X*diag(P);
   if norm(Hm) > 100*epp
      disp('Warning: inaccurate eigenvalue computation.')
      Hm = Hm
      Test = 5
   end
   X = conj(X);
   Hm = X'*A-diag(P)*X';
   if norm(Hm) > 100*epp
      disp('Warning: inaccurate left-eigenvalue computation.')
      Hm = Hm
   end
   
   Npd = max(size(P));
   for i=1:Npd	 
      if abs(P(i)-real(P(i))) < epp 
         P(i) = real(P(i));
	 if isreal(X(:,i)) == 0 % Rotate eigenvector.
	    inz = find(abs(X(:,i))>epp); 
	    if isempty(inz) == 0     
	       vh = angle(X(inz(1),i));
               X(:,i)=X(:,i)*exp(-sqrt(-1)*vh);
	    end
            X(:,i) = zerm(X(:,i),epp/10);
         end
      end
      
      uph = B'*X(:,i); nuph = norm(uph);
      
      if nuph > epp
         U(:,i) = uph/nuph;
	 X(:,i) = X(:,i)/nuph;	 
	 S(i,1) = 1/nuph;
	 if (isreal(U(:,i))==0) & (rf==1)   % Rotate output direction. 
            inz = find(abs(U(:,i))>epp); 
	    if isempty(inz) == 0
	       vh = angle(U(inz(1),i));
               U(:,i) = U(:,i)*exp(-sqrt(-1)*vh); 
               X(:,i)=X(:,i)*exp(-sqrt(-1)*vh);
               S(i,1) = S(i,1)*exp(-sqrt(-1)*vh);
            end
         end
      else
         U(:,i) = zeros(nu,1);    % Norm is small uncontrollable.
	 S(i,1) = 0;
      end
   end
   P = zerm(P,epp/10); X = zerm(X,epp/10); U = zerm(U,epp/10);
end