1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) INRIA - Serge STEER
3 // Copyright (C) 1999 - Lucien.Povy@eudil.fr (to get a good table)
4 //
5 // This file must be used under the terms of the CeCILL.
6 // This source file is licensed as described in the file COPYING, which
7 // you should have received as part of this distribution.  The terms
8 // are also available at
9 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11 function r=routh_t(h,k)
12     //r=routh_t(h,k) computes routh table of denominator of the
13     //system described by transfer matrix SISO continue h with the
14     //feedback by the gain k
15     //If  k=poly(0,'k') we will have a polynomial matrix with dummy variable
16     //k, formal expression of the Routh table.
17     //r=routh_t(d) computes Routh table of h :attention ! d=denom of system
19     //If a zero row appears, it means that there exist a pair of pure imaginary
20     //roots (oscillating system) or symmetric real roots. In this case, the pure imaginary roots are the
21     //imaginary roots of the bisquare polynomial given by the previous row. The
22     //routh table can be continued replacing this row by the derivative of this
23     //polynomial.
24     //see http://www.jdotec.net/s3i/TD_Info/Routh/Routh.pdf for degenerated
25     //cases
27     //see also
28     //Comments on the Routh-Hurwitz criterion, Shamash, Y.,Automatic Control, IEEE T.A.C
29     //Volume 25, Issue 1, Feb 1980 Page(s): 132 - 133
31     //http://controls.engin.umich.edu/wiki/index.php/RouthStability
32     [lhs,rhs]=argn(0);
33     h1=h(1);
34     if rhs==2 then
35         if typeof(h)<>'rational' then
36             error(msprintf(gettext("%s: Wrong type for input argument #%d: rational fraction array expected.\n"),"routh_t",1));
37         end
38         [n,d]=h(2:3)
39         if size(n,'*')<>1 then
40             error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"routh_t",1))
41         end
42         nd=max([degree(d) degree(n)])+1;
43         cod=coeff(d,0:nd-1);//coeff du denominateur
44         con=coeff(n,0:nd-1);//coeff du numerateur
45         cobf=cod+k*con //coeff de la boucle fermee
46     else
47         if type(h)>2 then
48             error(msprintf(gettext("%s: Wrong type for input argument #%d: Polynomial array expected.\n"),"routh_t",1));
49         end
50         if size(h,'*')<>1 then
51             error(msprintf(gettext("%s: Wrong size for input argument #%d: A polynomial expected.\n"),"routh_t",1))
52         end
54         nd=degree(h)+1;
55         cobf=coeff(h,0:nd-1)
56     end;
57     //
58     r1=cobf(nd:-2:1);
59     r2=cobf(nd-1:-2:1);
60     ncol=size(r1,'*');
61     if size(r2,'*')<>ncol then r2=[r2,0],end
62     r=[r1;r2]
63     if ncol<2 then r=[],return,end;
64     if rhs==2 then
65         for i=3:nd,
66             r(i,1:ncol-1)=[1.,-r(i-2,1)/r(i-1,1)]*[r(i-2,2:ncol);r(i-1,2:ncol)]
67         end;
68     else
69         for i=3:nd,
70             if and(r(i-1,:)==0) then
71                 naux=nd-i+2 //order of previous polynomial
72                 exponents=naux:-2:0
73                 ncoeff=size(exponents,'*')
74                 r(i-1,1:ncoeff)=r(i-2,1:ncoeff).*exponents //derivative of previous polynomial
75             end
76             if r(i-1,1)==0 then
77                 if rhs==1 then
78                     if typeof(r)=='rational' then
79                         //scilab is not able to handle multivariable polynomials
80                         r=horner(r,%eps^2);
81                     end
82                     r(i-1,1)=poly(0,'eps')
83                 else
84                     r(i-1,1)=%eps^2,
85                 end
86             end
87             r(i,1:ncol-1)=[1.,-r(i-2,1)/r(i-1,1)]*[r(i-2,2:ncol);r(i-1,2:ncol)]
88         end;
89     end;
90 endfunction