home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / TZREDUCE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.6 KB  |  147 lines
  1. function [abcd,mu,nu] = tzreduce(abcd,m,n,p,Zeps,ro,sigma)
  2. %TZREDUCE Extract from the system (A,B,C,D) the reduced system 
  3. %    (a,b,c,d) but with d full row rank.  Used in TZERO.
  4. %
  5. %    [ABCD.MU,NU] = TZREDUCE(ABCD,M,N,P,Zeps,RO,SIGMA)
  6.  
  7. %    Clay M. Thompson  7-23-90
  8. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  9.  
  10. %  Extracts from the (N+P)x(M+N) system [ B A ],  a (NU+MU)x(M+NU) 'reduced' 
  11. %         [ B' A' ]                     [ D C ]
  12. %  system [ D' C' ] having the same transmission zeros but with D' Full Row 
  13. %  Rank.  The system [ A' B' C' D' ] overwrites the old system.
  14. %
  15. % Reference: Adapted from "Computation of Zeros of Linear Multivariable
  16. %            Systems", A. Emami-Naeini, and P. Van Dooren; Automatica
  17. %            Vol. 18, No. 4, pp. 415-430, 1982.
  18.  
  19. % *** initialize ***
  20. Sum2 = zeros(1,max(p,m));
  21. mu=p;
  22. nu=n;
  23. while mu~=0,
  24.  
  25.   ro1=ro;
  26.   mnu=m+nu;
  27.   numu=nu+mu;
  28.  
  29.   if m~=0,
  30.     ro1=ro1+1;
  31.     irow=nu;
  32.  
  33.     % *** Compress rows of D(*).  First exploit triangular shape ***
  34.     for icol=1:sigma-1
  35.       rows = irow + [1:ro1];
  36.       dummy = abcd(rows,icol);
  37.       [dummy,s,zero] = housh(dummy,1,Zeps);
  38.       abcd(rows,icol:mnu) = (eye(ro1)-s*dummy*dummy')*abcd(rows,icol:mnu);
  39.       irow=irow+1;
  40.     end;
  41.  
  42.     % *** Continue householder with pivoting ***
  43.  
  44.     if sigma==0,
  45.       sigma=1;
  46.       ro1=ro1-1;
  47.     end
  48.  
  49.     if sigma~=m,
  50.       if ro1==1,
  51.         Sum2(sigma:m) = abcd(irow+1,sigma:m).*abcd(irow+1,sigma:m);
  52.       else
  53.         Sum2(sigma:m) = sum(abcd(irow+1:irow+ro1,sigma:m).*abcd(irow+1:irow+ro1,sigma:m));
  54.       end
  55.     end
  56.  
  57.     for icol=sigma:m;
  58.  
  59.       % *** Pivot if necessary ***
  60.  
  61.       if icol~=m,
  62.         Rows = 1:numu;
  63.  
  64.         [dum,ibar] = max(Sum2(icol:m));
  65.         ibar=ibar+icol-1;
  66.         if ibar~=icol,
  67.           Sum2(ibar)=Sum2(icol); Sum2(icol)=dum;
  68.           dum=abcd(Rows,icol);
  69.           abcd(Rows,icol)=abcd(Rows,ibar);
  70.          abcd(Rows,ibar)=dum;
  71.         end
  72.       end
  73.  
  74.       % *** Perform Householder transformation ***
  75.  
  76.       dummy=abcd(irow+1:irow+ro1,icol);
  77.       [dummy,s,zero] = housh(dummy,1,Zeps);
  78.       if zero, break, end
  79.       if ro1==1, return, end
  80.       abcd(irow+1:irow+ro1,icol:mnu) = (eye(ro1)-s*dummy*dummy')*abcd(irow+1:irow+ro1,icol:mnu);
  81.       irow=irow+1;
  82.       ro1=ro1-1;
  83.       Sum2(icol:m) = Sum2(icol:m) - abcd(irow,icol:m) .* abcd(irow,icol:m);
  84.     end % if
  85.  
  86.   end % for
  87.   tau=ro1;
  88.   sigma=mu-tau;
  89.  
  90.   % *** Compress the columns of C(*) ***
  91.   if (nu<=0), mu=sigma; nu=0; return, end
  92.  
  93.   i1=nu+sigma;
  94.   mm1=m+1;
  95.   n1=nu;
  96.   if tau~=1,
  97.     if mm1==mnu,
  98.       Sum2(1:tau) = (abcd(i1+1:i1+tau,mm1).*abcd(i1+1:i1+tau,mm1))';
  99.     else
  100.       Sum2(1:tau) = sum((abcd(i1+1:i1+tau,mm1:mnu).*abcd(i1+1:i1+tau,mm1:mnu))');
  101.     end
  102.   end
  103.  
  104.   for ro1=1:tau;
  105.     ro=ro1-1;
  106.     i=tau-ro;
  107.     i2=i+i1;
  108.  
  109.     % *** Pivot if necessary ***
  110.  
  111.     if i~=1,
  112.       [dum,ibar] = max(Sum2(1:i));
  113.       if ibar~=i,
  114.         Sum2(ibar) = Sum2(i); Sum2(i) = dum;
  115.         dum = abcd(i2,mm1:mnu);
  116.         abcd(i2,mm1:mnu) = abcd(ibar+i1,mm1:mnu);
  117.         abcd(ibar+i1,mm1:mnu) = dum;
  118.       end
  119.     end
  120.  
  121.     % *** Perform Householder Transformation ***
  122.     
  123.     cols = m + [1:n1];
  124.     dummy = abcd(i2,cols);
  125.     [dummy,s,zero] = housh(dummy,n1,Zeps);
  126.     if zero, break, end
  127.     if n1~=1
  128.       abcd(1:i2,cols) = abcd(1:i2,cols)*(eye(n1)-s*dummy*dummy');
  129.       mn1=m+n1;
  130.       abcd(1:n1,1:mn1) = (eye(n1)-s*dummy*dummy')*abcd(1:n1,1:mn1);
  131.       Sum2(1:i) = Sum2(1:i)-abcd(i1+1:i1+i,mn1)' .* abcd(i1+1:i1+i,mn1)';
  132.       mnu=mnu-1;
  133.     end
  134.     n1=n1-1;
  135.  
  136.   end % for
  137.  
  138.   if ~zero, ro=tau; end
  139.  
  140.   nu=nu-ro;
  141.   mu=sigma+ro;
  142.   
  143.   if ro==0, return, end
  144. end % while
  145.  
  146.  
  147.