home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / TZERO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.0 KB  |  108 lines
  1. function [z,gain] = tzero(a,b,c,d)
  2. %TZERO    Transmission zeros.
  3. %    Z = TZERO(A,B,C,D) returns the transmission zeros of the state-
  4. %    space system:   .
  5. %                       x = Ax + Bu   or   x[n+1] = Ax[n] + Bu[n]
  6. %                       y = Cx + Du          y[n] = Cx[n] + Du[n]
  7. %
  8. %    [Z,GAIN] = TZERO(A,B,C,D) also returns the transfer function 
  9. %    gain if the system is SISO.
  10. %
  11. %    See also: PZMAP, SS2ZP, and ZP2SS.
  12.  
  13. %    Clay M. Thompson  7-23-90
  14. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  15.  
  16. % Extracts from the system matrix of a state-space system [ A B C D ] a regular
  17. % pencil [ Lambda*Bf - Af ] which has the NU Invariant Zeros of the system as
  18. % Generalized Eigenvalues.
  19. %  
  20. %  Reference: Adapted from "Computation of Zeros of Linear Multivariable
  21. %             Systems", A. Emami-Naeini, and P. Van Dooren; Automatica 
  22. %             Vol. 18, No. 4, pp. 415-430, 1982.
  23.  
  24. % The transmission zero calculation can be tested by checking that the
  25. % matrix: [A B;C D]-lambda*[I 0; 0 0] is rank deficient, where lambda
  26. % is an element of Z.
  27.  
  28. error(nargchk(4,4,nargin));
  29. error(abcdchk(a,b,c,d));
  30.  
  31. if isempty(b)&isempty(c), z=[]; gain=0; return, end
  32.  
  33. Zeps = 2*eps*norm(a,'fro');    % Epsilon used for transmission zero calculation.
  34.  
  35. [nn,nx]=size(a);
  36. [pp,na]=size(c);
  37. [na,mm]=size(b);
  38. if pp*mm==1, tf=1; else tf=0; end
  39.  
  40. % *** Construct the Compound Matrix [ B A ] of Dimension (N+P)x(M+N) 
  41. %                                   [ D C ]
  42.  
  43. bf = [b,a;d,c];
  44.  
  45. % *** Reduce this system to one with the same invariant zeros and with
  46. %     D(*) full rank MU (The Normal Rank of the original system) ***
  47.  
  48. [bf,mu,nu] = tzreduce(bf,mm,nn,pp,Zeps,pp,0);
  49. Rank=mu;
  50.  
  51. if nu==0, 
  52.   z = [];
  53. else
  54.   %  *** Pretranspose the system *** 
  55.   mnu = mm+nu;
  56.   numu = nu+mu;
  57.   af = zeros(mnu,numu);
  58.   af(mnu:-1:1,numu:-1:1) = bf(1:numu,1:mnu)';
  59.  
  60.   if mu~=mm,
  61.     pp=mm;
  62.     nn=nu;
  63.     mm=mu;
  64.  
  65.     % *** Reduce the system to one with the same invariant zeros and with }
  66.     %     D(*) square invertible *** }
  67.  
  68.     [af,mu,nu] = tzreduce(af,mm,nn,pp,Zeps,pp-mm,mm);
  69.     mnu=mm+nu;
  70.   end
  71.  
  72.   if nu==0,
  73.     z = [];
  74.   else
  75.     % *** Perform a unitary transformation on the columns of [ sI-A B ] in }
  76.     %                           [ sBf-Af X ]                 [   -C D ] }
  77.     %     order to reduce it to [   0    Y ] with Y & Bf square invertible }
  78.  
  79.     bf(1:nu,1:mnu)=[zeros(nu,mm),eye(nu)];
  80.     if (Rank~=0)
  81.       nu1 = nu+1;
  82.       i1=nu+mu;
  83.       i0 = mm;
  84.       for i=1:mm
  85.         i0  = i0-1;
  86.         cols = i0 + [1:nu1];
  87.         dummy = af(i1,cols);
  88.         [dummy,s,zero] = housh(dummy,nu1,Zeps);
  89.         af(1:i1,cols) = af(1:i1,cols)*(eye(nu1)-s*dummy*dummy');
  90.         bf(1:nu,cols) = bf(1:nu,cols)*(eye(nu1)-s*dummy*dummy');
  91.         i1 = i1-1;
  92.       end % for
  93.     end % if
  94.  
  95.     % -- Solve Generalized zeros of sBF - AF
  96.     z = eig(af(1:nu,1:nu)/bf(1:nu,1:nu));
  97.   end
  98.  
  99. end
  100.  
  101. if (nargout==2)&tf, % Compute transfer function gain
  102.   if nu==nx,
  103.     gain=bf(nu+1,1);
  104.   else
  105.     gain=bf(nu+1,1)*prod(diag(bf(nu+2:nx+1,nu+2:nx+1)));
  106.   end
  107. end
  108.