* Bug #12174 fixed - The function "routh_t" gives incorrect output for all denominato... 45/10345/6
Charlotte HECQUET [Thu, 31 Jan 2013 12:22:42 +0000 (13:22 +0100)]
gain value "k".

Change-Id: I4c0e2a2b4d2869b117c41953d9e9c438c0c50f25

scilab/CHANGES_5.4.X
scilab/modules/cacsd/macros/routh_t.sci
scilab/modules/cacsd/tests/nonreg_tests/bug_12174.dia.ref [new file with mode: 0644]
scilab/modules/cacsd/tests/nonreg_tests/bug_12174.tst [new file with mode: 0644]

index ddc6fba..66e2ff0 100644 (file)
@@ -294,6 +294,9 @@ Bug fixes
 
 * Bug #12168 fixed - matfile_listvar crashed when listing variables of a closed MAT-file.
 
+* Bug #12174 fixed - The function "routh_t" gave incorrect output for all
+                     denominators that include gain value "k".
+
 * Bug #12179 fixed - Fix an incompatibility with MPI version of HDF5.
 
 * Bug #12184 fixed - Performances of the function 'derivat' improved.
index 64fb0a2..7b61f2a 100644 (file)
@@ -9,84 +9,82 @@
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 function r=routh_t(h,k)
-//r=routh_t(h,k) computes routh table of denominator of the
-//system described by transfer matrix SISO continue h with the
-//feedback by the gain k
-//If  k=poly(0,'k') we will have a polynomial matrix with dummy variable
-//k, formal expression of the Routh table.
-//r=routh_t(d) computes Routh table of h :attention ! d=denom of system
+    //r=routh_t(h,k) computes routh table of denominator of the
+    //system described by transfer matrix SISO continue h with the
+    //feedback by the gain k
+    //If  k=poly(0,'k') we will have a polynomial matrix with dummy variable
+    //k, formal expression of the Routh table.
+    //r=routh_t(d) computes Routh table of h :attention ! d=denom of system
 
-//If a zero row appears, it means that there exist a pair of pure imaginary
-//roots (oscillating system) or symmetric real roots. In this case, the pure imaginary roots are the
-//imaginary roots of the bisquare polynomial given by the previous row. The
-//routh table can be continued replacing this row by the derivative of this
-//polynomial.
-//see http://www.jdotec.net/s3i/TD_Info/Routh/Routh.pdf for degenerated
-//cases
-  
-  //see also 
-//Comments on the Routh-Hurwitz criterion, Shamash, Y.,Automatic Control, IEEE T.A.C
-//Volume 25, Issue 1, Feb 1980 Page(s): 132 - 133
-  
-  //http://controls.engin.umich.edu/wiki/index.php/RouthStability
-  [lhs,rhs]=argn(0);
-  h1=h(1);
-  if rhs==2 then
-    if typeof(h)<>'rational' then
-      error(msprintf(gettext("%s: Wrong type for input argument #%d: rational fraction array expected.\n"),"routh_t",1));
-    end
-    [n,d]=h(2:3)
-    if size(n,'*')<>1 then
-      error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"routh_t",1))
-    end
-    nd=max([degree(d) degree(n)])+1;
-    cod=coeff(d,0:nd-1);//coeff du denominateur
-    con=coeff(n,0:nd-1);//coeff du numerateur
-    cobf=cod+k*con //coeff de la boucle fermee
-  else
-    if type(h)>2 then 
-      error(msprintf(gettext("%s: Wrong type for input argument #%d: Polynomial array expected.\n"),"routh_t",1));
-    end
-    if size(h,'*')<>1 then
-      error(msprintf(gettext("%s: Wrong size for input argument #%d: A polynomial expected.\n"),"routh_t",1))
-    end
+    //If a zero row appears, it means that there exist a pair of pure imaginary
+    //roots (oscillating system) or symmetric real roots. In this case, the pure imaginary roots are the
+    //imaginary roots of the bisquare polynomial given by the previous row. The
+    //routh table can be continued replacing this row by the derivative of this
+    //polynomial.
+    //see http://www.jdotec.net/s3i/TD_Info/Routh/Routh.pdf for degenerated
+    //cases
 
-    nd=degree(h)+1;
-    cobf=coeff(h,0:nd-1)
-  end;
-  //
-  r1=cobf(nd:-2:1);
-  r2=cobf(nd-1:-2:1);
-  ncol=size(r1,'*');
-  if size(r2,'*')<>ncol then r2=[r2,0],end
-  r=[r1;r2]
-  if ncol<2 then r=[],return,end;
-  if rhs==2 then
+    //see also
+    //Comments on the Routh-Hurwitz criterion, Shamash, Y.,Automatic Control, IEEE T.A.C
+    //Volume 25, Issue 1, Feb 1980 Page(s): 132 - 133
 
-    for i=3:nd,
-      r(i,1:ncol-1)=[r(i-1,1),-r(i-2,1)]*[r(i-2,2:ncol);r(i-1,2:ncol)]
+    //http://controls.engin.umich.edu/wiki/index.php/RouthStability
+    [lhs,rhs]=argn(0);
+    h1=h(1);
+    if rhs==2 then
+        if typeof(h)<>'rational' then
+            error(msprintf(gettext("%s: Wrong type for input argument #%d: rational fraction array expected.\n"),"routh_t",1));
+        end
+        [n,d]=h(2:3)
+        if size(n,'*')<>1 then
+            error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"routh_t",1))
+        end
+        nd=max([degree(d) degree(n)])+1;
+        cod=coeff(d,0:nd-1);//coeff du denominateur
+        con=coeff(n,0:nd-1);//coeff du numerateur
+        cobf=cod+k*con //coeff de la boucle fermee
+    else
+        if type(h)>2 then
+            error(msprintf(gettext("%s: Wrong type for input argument #%d: Polynomial array expected.\n"),"routh_t",1));
+        end
+        if size(h,'*')<>1 then
+            error(msprintf(gettext("%s: Wrong size for input argument #%d: A polynomial expected.\n"),"routh_t",1))
+        end
+
+        nd=degree(h)+1;
+        cobf=coeff(h,0:nd-1)
     end;
-  else
-    for i=3:nd,
-     if and(r(i-1,:)==0) then
-       naux=nd-i+2 //order of previous polynomial
-       exponents=naux:-2:0
-       ncoeff=size(exponents,'*')
-       r(i-1,1:ncoeff)=r(i-2,1:ncoeff).*exponents //derivative of previous polynomial
-      end
-      if r(i-1,1)==0 then 
-       if rhs==1 then
-         if typeof(r)=='rational' then 
-           //scilab is not able to handle multivariable polynomials
-           r=horner(r,%eps^2); 
-         end
-         r(i-1,1)=poly(0,'eps')
-       else
-         r(i-1,1)=%eps^2,
-       end
-      end
-      r(i,1:ncol-1)=[1.,-r(i-2,1)/r(i-1,1)]*[r(i-2,2:ncol);r(i-1,2:ncol)]
-      
+    //
+    r1=cobf(nd:-2:1);
+    r2=cobf(nd-1:-2:1);
+    ncol=size(r1,'*');
+    if size(r2,'*')<>ncol then r2=[r2,0],end
+    r=[r1;r2]
+    if ncol<2 then r=[],return,end;
+    if rhs==2 then
+        for i=3:nd,
+            r(i,1:ncol-1)=[1.,-r(i-2,1)/r(i-1,1)]*[r(i-2,2:ncol);r(i-1,2:ncol)]
+        end;
+    else
+        for i=3:nd,
+            if and(r(i-1,:)==0) then
+                naux=nd-i+2 //order of previous polynomial
+                exponents=naux:-2:0
+                ncoeff=size(exponents,'*')
+                r(i-1,1:ncoeff)=r(i-2,1:ncoeff).*exponents //derivative of previous polynomial
+            end
+            if r(i-1,1)==0 then
+                if rhs==1 then
+                    if typeof(r)=='rational' then
+                        //scilab is not able to handle multivariable polynomials
+                        r=horner(r,%eps^2);
+                    end
+                    r(i-1,1)=poly(0,'eps')
+                else
+                    r(i-1,1)=%eps^2,
+                end
+            end
+            r(i,1:ncol-1)=[1.,-r(i-2,1)/r(i-1,1)]*[r(i-2,2:ncol);r(i-1,2:ncol)]
+        end;
     end;
-  end;
 endfunction
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_12174.dia.ref b/scilab/modules/cacsd/tests/nonreg_tests/bug_12174.dia.ref
new file mode 100644 (file)
index 0000000..dd5e591
--- /dev/null
@@ -0,0 +1,28 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- Non-regression test for bug 12174 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12174
+//
+// <-- Short Description -->
+// The function "routh_t" gives incorrect output for all denominators that include
+//gain value "k".
+s = poly(0, 's');
+k=poly(0,'k');
+P=s*(s+7)*(s+11);
+h=1/P;
+r=routh_t(h,k);
+ref=[1 77; 18 k;(1386-k)/18 0; k 0];
+numref=numer(ref);
+numref(3)=1386 - k;
+denr=denom(r);
+denr=horner(denr,1);
+denref=denom(ref);
+denref(3)=18;
+assert_checkequal(numer(r)/denr,numref/denref);
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_12174.tst b/scilab/modules/cacsd/tests/nonreg_tests/bug_12174.tst
new file mode 100644 (file)
index 0000000..118ec0b
--- /dev/null
@@ -0,0 +1,31 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- Non-regression test for bug 12174 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12174
+//
+// <-- Short Description -->
+// The function "routh_t" gives incorrect output for all denominators that include
+//gain value "k".
+
+s = poly(0, 's');
+k=poly(0,'k');
+P=s*(s+7)*(s+11);
+h=1/P;
+r=routh_t(h,k);
+
+ref=[1 77; 18 k;(1386-k)/18 0; k 0];
+numref=numer(ref);
+numref(3)=1386 - k;
+denr=denom(r);
+denr=horner(denr,1);
+denref=denom(ref);
+denref(3)=18;
+
+assert_checkequal(numer(r)/denr,numref/denref);