home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / TZERO2.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.8 KB  |  67 lines
  1. function z = tzero(a,b,c,d)
  2. %TZERO2    Transmission zeros.
  3. %    TZERO2(A,B,C,D) returns the transmission zeros of the state-space
  4. %    system:
  5. %        .
  6. %        x = Ax + Bu
  7. %        y = Cx + Du
  8. %
  9. %    Care must be exercised to disregard any large zeros that may
  10. %    actually be at infinity.  It is best to use FORMAT SHORT E so
  11. %    large zeros don't swamp out the display of smaller finite zeros.
  12. %
  13. %    This M-file finds the transmission zeros using an algorithm based
  14. %    on QZ techniques.  It is not as numerically reliable as the 
  15. %    algorithm in TZERO.
  16. %
  17. %    See also: TZERO.
  18.  
  19. % For more information on this algorithm, see:
  20. %   [1] A.J Laub and B.C. Moore, Calculation of Transmission Zeros Using
  21. %   QZ Techniques, Automatica, 14, 1978, p557
  22. %
  23. % For a better algorithm, not implemented here, see:
  24. %   [2] A. Emami-Naeini and P. Van Dooren, Computation of Zeros of Linear 
  25. %   Multivariable Systems, Automatica, Vol. 14, No. 4, pp. 415-430, 1982.
  26. %
  27.  
  28. %    J.N. Little 4-21-85
  29. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  30.  
  31. error(nargchk(4,4,nargin));
  32. error(abcdchk(a,b,c,d));
  33.  
  34. [n,m] = size(b);
  35. [r,n] = size(c);
  36. aa = [a b;c d];
  37.  
  38. if m == r
  39.     [mcols, nrows] = size(aa);
  40.     bb = zeros(mcols, nrows);
  41.     bb(1:n,1:n) = eye(n);
  42.     z = eig(aa,bb);
  43.     z = z(~isnan(z) & finite(z)); % Punch out NaN's and Inf's
  44. else
  45.     nrm = norm(aa,1);
  46.     if m > r
  47.         aa1 = [aa;nrm*(rand(m-r,n+m)-.5)];
  48.         aa2 = [aa;nrm*(rand(m-r,n+m)-.5)];
  49.     else
  50.         aa1 = [aa nrm*(rand(n+r,r-m)-.5)];
  51.         aa2 = [aa nrm*(rand(n+r,r-m)-.5)];
  52.     end
  53.     [mcols, nrows] = size(aa1);
  54.     bb = zeros(mcols, nrows);
  55.     bb(1:n,1:n) = eye(n);
  56.     z1 = eig(aa1,bb);
  57.     z2 = eig(aa2,bb);
  58.     z1 = z1(~isnan(z1) & finite(z1));   % Punch out NaN's and Inf's
  59.     z2 = z2(~isnan(z2) & finite(z2));
  60.     nz = length(z1);
  61.     for i=1:nz
  62.         if any(abs(z1(i)-z2) < nrm*sqrt(eps))
  63.             z = [z;z1(i)];
  64.         end
  65.     end
  66. end
  67.