home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / MINREAL.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.2 KB  |  94 lines
  1. function [am,bm,cm,dm] = minreal(a,b,c,d,tol)
  2. %MINREAL Minimal realization and pole-zero cancellation.
  3. %    [Am,Bm,Cm,Dm] = MINREAL(A,B,C,D) returns a minimal realization
  4. %    of the state-space system {A,B,C,D}.  A message is displayed 
  5. %    indicating the number of states removed.
  6. %    [Am,Bm,Cm,Dm] = MINREAL(A,B,C,D,TOL) uses the tolerance TOL
  7. %    in deciding which states to eliminate.
  8. %
  9. %    [Zm,Pm] = MINREAL(Z,P), where Z and P are column vectors
  10. %    containing poles and zeros, cancels the common roots that
  11. %    are within TOL = 10*SQRT(EPS)*ABS(Z(i)) of each other.
  12. %    [Zm,Pm] = MINREAL(Z,P,TOL) uses tolerance TOL.
  13. %
  14. %    For transfer functions, [NUMm,DENm] = MINREAL(NUM,DEN), where
  15. %    NUM and DEN are row vectors of polynomial coefficients, cancels
  16. %    the common roots in the polynomials.
  17. %    [NUMm,DENm] = MINREAL(NUM,DEN,TOL) uses tolerance TOL.
  18.  
  19. %    J.N. Little 7-17-86
  20. %    Revised A.C.W.Grace 12-1-89
  21. %    Copyright (c) 1986-93 by The MathWorks, Inc.
  22.  
  23. disp(' ')
  24.  
  25. if nargin == 2 | nargin == 3
  26.     z = a;
  27.     p = b;
  28.     [mz,nz] = size(z);
  29.     [mp,np] = size(p);
  30.     if (mz == 1) & (mp == 1)
  31.         % If transfer function, convert to zero-pole:
  32.         [z,p,k] = tf2zp(z,p);
  33.     end
  34.  
  35.     % Strip infinities from zeros and throw away.
  36.     z = z(finite(z));
  37.  
  38.     mz = max(size(z));
  39.     mp = max(size(p));
  40.     iz = ones(mz,1);
  41.  
  42.     % Loop through zeros, looking for matching poles:
  43.     for i=1:mz
  44.         zi = z(i);
  45.         if (nargin == 2)
  46.             tol = 10*abs(zi)*sqrt(eps);
  47.         else
  48.             tol=c;
  49.         end
  50.         kk = find(abs(p-zi) <= tol);
  51.         if all(size(kk))
  52.             p(kk(1)) = [];
  53.             iz(i) = 0;
  54.         end
  55.     end
  56.  
  57.     % Eliminate matches in zeros:
  58.     z = z(iz);
  59.  
  60.     % If transfer function, convert back.
  61.     if (nz > 1) | (np > 1)
  62.         [z,p] = zp2tf(z,p,k);
  63.     end
  64.     am = z;
  65.     bm = p;
  66.     disp([int2str(mz-sum(iz)),' pole-zeros cancelled'])
  67.     return
  68. end
  69.  
  70. % Do state-space case
  71. [ns,nu] = size(b);
  72. if nargin == 4
  73.     tol = 10*ns*norm(a,1)*eps;
  74. end
  75. [am,bm,cm,t,k] = ctrbf(a,b,c,tol);
  76. kk = sum(k);
  77. nu = ns - kk;
  78. nn = nu;
  79. am = am(nu+1:ns,nu+1:ns);
  80. bm = bm(nu+1:ns,:);
  81. cm = cm(:,nu+1:ns);
  82. ns = ns - nu;
  83. if ns
  84.     [am,bm,cm,t,k] = obsvf(am,bm,cm,tol);
  85.     kk = sum(k);
  86.     nu = ns - kk;
  87.     nn = nn + nu;
  88.     am = am(nu+1:ns,nu+1:ns);
  89.     bm = bm(nu+1:ns,:);
  90.     cm = cm(:,nu+1:ns);
  91. end
  92. disp([int2str(nn),' states removed'])
  93. dm = d;
  94.