home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / tzero.m < prev    next >
Text File  |  1999-12-24  |  4KB  |  126 lines
  1. ## Copyright (C) 1996 Auburn University.  All Rights Reserved
  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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. 
  18.  
  19. ##  -*- texinfo -*-
  20. ## @deftypefn {Function File} {} tzero (@var{a}, @var{b}, @var{c}, @var{d}@{, @var{opt}@})
  21. ## @deftypefnx {Function File} {} tzero (@var{sys}@{,@var{opt}@})
  22. ##  Compute transmission zeros of a continuous
  23. ## @example
  24. ## .
  25. ## x = Ax + Bu
  26. ## y = Cx + Du
  27. ## @end example
  28. ## or discrete
  29. ## @example
  30. ## x(k+1) = A x(k) + B u(k)
  31. ## y(k)   = C x(k) + D u(k)
  32. ## @end example
  33. ## system.
  34. ## @strong{Outputs}
  35. ## @table @var
  36. ## @item zer
  37. ##  transmission zeros of the system
  38. ## @item gain
  39. ## leading coefficient (pole-zero form) of SISO transfer function
  40. ## returns gain=0 if system is multivariable
  41. ## @end table
  42. ## @strong{References}
  43. ## @enumerate
  44. ## @item Emami-Naeini and Van Dooren, Automatica, 1982.
  45. ## @item Hodel, "Computation of Zeros with Balancing," 1992 Lin. Alg. Appl.
  46. ## @end enumerate
  47. ## @end deftypefn
  48.  
  49.  
  50. function [zer, gain] = tzero (A, B, C, D)
  51.  
  52.   ## R. Bruce Tenison July 4, 1994
  53.   ## A. S. Hodel Aug 1995: allow for MIMO and system data structures
  54.  
  55.   ## get A,B,C,D and Asys variables, regardless of initial form
  56.   if(nargin == 4)
  57.     Asys = ss2sys(A,B,C,D);
  58.   elseif( (nargin == 1) && (! is_struct(A)))
  59.     usage("[zer,gain] = tzero(A,B,C,D) or zer = tzero(Asys)");
  60.   elseif(nargin != 1)
  61.     usage("[zer,gain] = tzero(A,B,C,D) or zer = tzero(Asys)");
  62.   else
  63.     Asys = A;
  64.     [A,B,C,D] = sys2ss(Asys);
  65.   endif
  66.  
  67.   Ao = Asys;            # save for leading coefficient
  68.   siso = is_siso(Asys);
  69.   digital = is_digit(Asys);    # check if it's mixed or not
  70.  
  71.   ## see if it's a gain block
  72.   if(isempty(A))
  73.     zer = [];
  74.     gain = D;
  75.     return;
  76.   endif
  77.  
  78.   ## First, balance the system via the zero computation generalized eigenvalue
  79.   ## problem balancing method (Hodel and Tiller, Linear Alg. Appl., 1992)
  80.  
  81.   Asys = zgpbal(Asys); [A,B,C,D] = sys2ss(Asys);   # balance coefficients
  82.   meps = 2*eps*norm([A, B; C, D],'fro');
  83.   Asys = zgreduce(Asys,meps);  [A, B, C, D] = sys2ss(Asys); # ENVD algorithm
  84.   if(!isempty(A))
  85.     ## repeat with dual system
  86.     Asys = ss2sys(A', C', B', D');   Asys = zgreduce(Asys,meps);
  87.  
  88.     ## transform back
  89.     [A,B,C,D] = sys2ss(Asys);    Asys = ss2sys(A', C', B', D');
  90.   endif
  91.  
  92.   zer = [];            # assume none
  93.   [A,B,C,D] = sys2ss(Asys);
  94.   if( !isempty(C) )
  95.     [W,r,Pi] = qr([C, D]');
  96.     [nonz,ztmp] = zgrownor(r,meps);
  97.     if(nonz)
  98.       ## We can now solve the generalized eigenvalue problem.
  99.       [pp,mm] = size(D);
  100.       nn = rows(A);
  101.       Afm = [A , B ; C, D] * W';
  102.       Bfm = [eye(nn), zeros(nn,mm); zeros(pp,nn+mm)]*W';
  103.  
  104.       jdx = (mm+1):(mm+nn);
  105.       Af = Afm(1:nn,jdx);
  106.       Bf = Bfm(1:nn,jdx);
  107.       zer = qz(Af,Bf);
  108.     endif
  109.   endif
  110.   
  111.   mz = length(zer);
  112.   [A,B,C,D] = sys2ss(Ao);        # recover original system
  113.   ## compute leading coefficient
  114.   if ( (nargout == 2) && siso)
  115.     n = rows(A);
  116.     if ( mz == n)
  117.       gain = D;
  118.     elseif ( mz < n )
  119.       gain = C*(A^(n-1-mz))*B;
  120.     endif
  121.   else
  122.     gain = [];
  123.   endif
  124. endfunction
  125.  
  126.