home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / tzero.m < prev    next >
Text File  |  1999-04-29  |  3KB  |  114 lines
  1. # Copyright (C) 1996 A. Scottedward Hodel 
  2. #
  3. # This file is part of Octave. 
  4. #
  5. # Octave is free software; you can redistribute it and/or modify it 
  6. # under the terms of the GNU General Public License as published by the 
  7. # Free Software Foundation; either version 2, or (at your option) any 
  8. # later version. 
  9. # 
  10. # Octave is distributed in the hope that it will be useful, but WITHOUT 
  11. # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
  12. # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
  13. # for more details.
  14. # 
  15. # You should have received a copy of the GNU General Public License 
  16. # along with Octave; see the file COPYING.  If not, write to the Free 
  17. # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 
  18.  
  19. function [zer, gain] = tzero(A,B,C,D)
  20.   # [zer{,gain}] = tzero(A,B,C,D) -or-
  21.   # [zer{,gain}] = tzero(Asys)
  22.   # Compute transmission zeros of a continuous
  23.   #      .
  24.   #      x = Ax + Bu
  25.   #      y = Cx + Du
  26.   #
  27.   # or discrete
  28.   #      x(k+1) = A x(k) + B u(k)
  29.   #      y(k)   = C x(k) + D u(k)
  30.   #
  31.   # system.
  32.   #
  33.   # outputs: 
  34.   #   zer: transmission zeros of the system
  35.   #   gain: leading coefficient (pole-zero form) of SISO transfer function
  36.   #         returns gain=0 if system is multivariable
  37.   # References:
  38.   # Hodel, "Computation of Zeros with Balancing," 1992 Lin. Alg. Appl.
  39.   
  40.   # R. Bruce Tenison July 4, 1994
  41.   # A. S. Hodel Aug 1995: allow for MIMO and system data structures
  42.  
  43.   # get A,B,C,D and Asys variables, regardless of initial form
  44.   if(nargin == 4)
  45.     Asys = ss2sys(A,B,C,D);
  46.   elseif( (nargin == 1) && (! is_struct(A)))
  47.     usage("[zer,gain] = tzero(A,B,C,D) or zer = tzero(Asys)");
  48.   elseif(nargin != 1)
  49.     usage("[zer,gain] = tzero(A,B,C,D) or zer = tzero(Asys)");
  50.   else
  51.     Asys = A;
  52.     [A,B,C,D] = sys2ss(Asys);
  53.   endif
  54.  
  55.   Ao = Asys;            # save for leading coefficient
  56.   siso = is_siso(Asys);
  57.   digital = is_digit(Asys);    # check if it's mixed or not
  58.  
  59.   # see if it's a gain block
  60.   if(isempty(A))
  61.     zer = [];
  62.     gain = D;
  63.     return;
  64.   endif
  65.  
  66.   # First, balance the system via the zero computation generalized eigenvalue
  67.   # problem balancing method (Hodel and Tiller, Linear Alg. Appl., 1992)
  68.  
  69.   Asys = zgpbal(Asys); [A,B,C,D] = sys2ss(Asys);   # balance coefficients
  70.   meps = 2*eps*norm([A, B; C, D],'fro');
  71.   Asys = zgreduce(Asys,meps);  [A, B, C, D] = sys2ss(Asys); # ENVD algorithm
  72.   if(!isempty(A))
  73.     # repeat with dual system
  74.     Asys = ss2sys(A', C', B', D');   Asys = zgreduce(Asys,meps);
  75.  
  76.     # transform back
  77.     [A,B,C,D] = sys2ss(Asys);    Asys = ss2sys(A', C', B', D');
  78.   endif
  79.  
  80.   zer = [];            # assume none
  81.   [A,B,C,D] = sys2ss(Asys);
  82.   if( !isempty(C) )
  83.     [W,r,Pi] = qr([C, D]');
  84.     [nonz,ztmp] = zgrownor(r,meps);
  85.     if(nonz)
  86.       # We can now solve the generalized eigenvalue problem.
  87.       [pp,mm] = size(D);
  88.       nn = rows(A);
  89.       Afm = [A , B ; C, D] * W';
  90.       Bfm = [eye(nn), zeros(nn,mm); zeros(pp,nn+mm)]*W';
  91.  
  92.       jdx = (mm+1):(mm+nn);
  93.       Af = Afm(1:nn,jdx);
  94.       Bf = Bfm(1:nn,jdx);
  95.       zer = qz(Af,Bf);
  96.     endif
  97.   endif
  98.   
  99.   mz = length(zer);
  100.   [A,B,C,D] = sys2ss(Ao);        # recover original system
  101.   #compute leading coefficient
  102.   if ( (nargout == 2) && siso)
  103.     n = rows(A);
  104.     if ( mz == n)
  105.       gain = D;
  106.     elseif ( mz < n )
  107.       gain = C*(A^(n-1-mz))*B;
  108.     endif
  109.   else
  110.     gain = [];
  111.   endif
  112. endfunction
  113.  
  114.