* Bug #12829 fixed - New optional output argument for routh_t function (+fix for... 32/12432/3
Charlotte HECQUET [Thu, 5 Sep 2013 12:28:58 +0000 (14:28 +0200)]
Change-Id: I125014ea2a23679189e3f022aa8a3972875831f2

SEP/INDEX
SEP/SEP_104_routh_t_optional_output.odt [new file with mode: 0644]
scilab/CHANGES_5.5.X
scilab/modules/cacsd/help/en_US/routh_t.xml
scilab/modules/cacsd/macros/routh_t.sci
scilab/modules/cacsd/tests/nonreg_tests/bug_12828.dia.ref [new file with mode: 0644]
scilab/modules/cacsd/tests/nonreg_tests/bug_12828.tst [new file with mode: 0644]
scilab/modules/cacsd/tests/nonreg_tests/bug_12829.dia.ref [new file with mode: 0644]
scilab/modules/cacsd/tests/nonreg_tests/bug_12829.tst [new file with mode: 0644]
scilab/modules/cacsd/tests/unit_tests/routh_t.dia.ref
scilab/modules/cacsd/tests/unit_tests/routh_t.tst

index 4aa290c..0b62131 100644 (file)
--- a/SEP/INDEX
+++ b/SEP/INDEX
@@ -99,5 +99,5 @@ SEP #100: Extends history size
 SEP #101: File browser
 SEP #102: members function.
 SEP #103: Localization : allow multiple domains and add addlocalizationdomain function.
-SEP #104: /* RESERVED */
+SEP #104: New optional output argument for routh_t (num)
 SEP #105: Three new optional output arguments to optim().
diff --git a/SEP/SEP_104_routh_t_optional_output.odt b/SEP/SEP_104_routh_t_optional_output.odt
new file mode 100644 (file)
index 0000000..cf43ffd
Binary files /dev/null and b/SEP/SEP_104_routh_t_optional_output.odt differ
index 559e537..3c38015 100644 (file)
@@ -83,6 +83,9 @@ Improvements
 * nthroot is now vectorizable.
   See bug #12678.
 
+* New optional output argument for routh_t.
+  See bug #12829.
+
 
 Obsolete
 =========
@@ -654,6 +657,8 @@ Bug Fixes
 
 * Bug #12800 fixed - Typo fixed in Polynomials help page.
 
+* Bug #12804 fixed - Typos fixed in routh_t help page.
+
 * Bug #12807 fixed - Display of showprofile improved.
 
 * Bug #12808 fixed - Add missing </td> in documentation generation (note, warning, ...).
@@ -672,6 +677,10 @@ Bug Fixes
 
 * Bug #12827 fixed - noisgen help page improved.
 
+* Bug #12828 fixed - routh_t gave a wrong result if the first element of a row was zero.
+
+* Bug #12829 fixed - New optional output argument added for routh_t function.
+
 * Bug #12830 fixed - In SciNotes, it was not possible to execute a replace action from the caret position.
 
 * Bug #12831 fixed - In SciNotes toolbar, there was no button to open code navigator.
index 5123456..193806a 100644 (file)
@@ -7,7 +7,7 @@
  * This source file is licensed as described in the file COPYING, which
  * you should have received as part of this distribution.  The terms
  * are also available at
- * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+ * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
  *
  -->
 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="routh_t" xml:lang="en">
     </refnamediv>
     <refsynopsisdiv>
         <title>Calling Sequence</title>
-        <synopsis>r=routh_t(p)
-            r=routh_t(h [,k])
-            r=routh_t(h [,k] [,normalized])
+        <synopsis>
+            [r [,num] ]=routh_t(p)
+            [r [,num] ]=routh_t(h ,kp)
+            r=routh_t(h ,k [,normalized])
         </synopsis>
     </refsynopsisdiv>
     <refsection>
                 </listitem>
             </varlistentry>
             <varlistentry>
+                <term>kp</term>
+                <listitem>
+                    <para>a scalar: proportional controller constant</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
                 <term>normalized</term>
                 <listitem>
                     <para>a boolean (%t (default value) or %f)</para>
             <varlistentry>
                 <term>r</term>
                 <listitem>
-                    <para>a matrix</para>
+                    <para>a matrix or a list: Routh array elements</para>
+                </listitem>
+            </varlistentry>
+            <varlistentry>
+                <term>num</term>
+                <listitem>
+                    <para>a scalar: the number of sign changes</para>
                 </listitem>
             </varlistentry>
         </variablelist>
@@ -61,7 +74,7 @@
         <title>Description</title>
         <para>
             <literal>r=routh_t(p)</literal> computes Routh's table of the
-            polynomial <literal>h</literal>.
+            polynomial <literal>p</literal>.
         </para>
         <para>
             <literal>r=routh_t(h,k)</literal> computes Routh's table of
@@ -71,7 +84,7 @@
         </para>
         <para>
             If <literal>k=poly(0,'k')</literal> we will have a polynomial or
-            a rational matrix with dummy variable <literal>k</literal>, 
+            a rational matrix with dummy variable <literal>k</literal>,
             formal expression of the Routh table.
         </para>
         <para>
             with non normalized elements. In the other case, we will have a rational
             and normalized matrix.
         </para>
+        <para>
+            The second output argument <literal>num</literal> returns the number of sign changes
+            in the first column of Routh's table. Note that, this argument value will only have sense
+            when the table is normalized.
+        </para>
+        <para>
+            <note>
+                Hint: If <literal>h=1/p</literal>, then <literal>routh_t(h, kp)</literal> is equivalent to
+                <literal>routh_t(p+kp)</literal> .
+            </note>
+        </para>
     </refsection>
     <refsection>
         <title>Examples</title>
-        <programlisting role="example"><![CDATA[ 
+        <programlisting role="example"><![CDATA[
 s=%s;
 P=5*s^3-10*s^2+7*s+20;
 routh_t(P)
 
 // Transfer function with formal feedback, normalized case
 routh_t((1+s)/P,poly(0,'k'))
-    
+
 // Transfer function with formal feedback, non normalized case
 routh_t((1+s)/P,poly(0,'k'),%f)
 
 // One of the coefficients in the polynomial equals zero
-P1=2*s^3-24*s+32; 
+P1=2*s^3-24*s+32;
 routh_t(P1)
 
 // A row full of zeros
 P2=s^4-6*s^3+10*s^2-6*s+9;
 routh_t(P2)
 
+// The number of roots in the rhp could be retrieved as the second output argument
+P3=5*s^3-10*s^2+7*s;
+[r,num]=routh_t(1/P3,20)
+if num==0
+   disp("System is stable")
+else
+   mprintf("There is %g sign changes in entries of first column.\nTherefore, system is unstable.", num)
+end
+
 //
  ]]></programlisting>
     </refsection>
@@ -111,10 +144,16 @@ routh_t(P2)
                 <link linkend="roots">roots</link>
             </member>
             <member>
+                <link linkend="spec">spec</link>
+            </member>
+            <member>
                 <link linkend="evans">evans</link>
             </member>
             <member>
-                <link linkend="spec">spec</link>
+                <link linkend="kpure">kpure</link>
+            </member>
+            <member>
+                <link linkend="krac2">krac2</link>
             </member>
         </simplelist>
     </refsection>
@@ -134,8 +173,12 @@ routh_t(P2)
         <title>History</title>
         <revhistory>
             <revision>
+                <revnumber>5.5.0</revnumber>
+                <revremark>New output argument added: num (SEP #104).</revremark>
+            </revision>
+            <revision>
                 <revnumber>5.4.0</revnumber>
-                <revremark> A new parameter added: normalized (SEP 89).</revremark>
+                <revremark>New input argument added: normalized (SEP #89).</revremark>
             </revision>
         </revhistory>
     </refsection>
index 68d3ed6..bc2fa69 100644 (file)
@@ -2,14 +2,15 @@
 // Copyright (C) INRIA - Serge STEER
 // Copyright (C) 1999 - Lucien.Povy@eudil.fr (to get a good table)
 // Copyright (C) 2013 - Charlotte HECQUET (new option)
-//
+// Copyright (C) 2013 - A. Khorshidi (to define a new optional output argument)
+// 
 // This file must be used under the terms of the CeCILL.
 // This source file is licensed as described in the file COPYING, which
-// you should have received as part of this distribution.  The terms
-// are also available at
-// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+// you should have received as par of this distribution.  The terms
+// are also available at    
+// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
-function r=routh_t(h,k,normalized)
+function [r,num]=routh_t(h,k,normalized)
     //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
@@ -25,7 +26,7 @@ function r=routh_t(h,k,normalized)
     //see http://www.jdotec.net/s3i/TD_Info/Routh/Routh.pdf for degenerated
     //cases
 
-    //see also
+    //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
 
@@ -42,11 +43,11 @@ function r=routh_t(h,k,normalized)
         rhs=2;
     end
     if rhs==2 then
-        if typeof(h)<>"rational" 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
+        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;
@@ -54,10 +55,10 @@ function r=routh_t(h,k,normalized)
         con=coeff(n,0:nd-1);//coeff du numerateur
         cobf=cod+k*con //coeff de la boucle fermee
     else
-        if type(h)>2 then
+        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
+        if size(h,'*')<>1 then
             error(msprintf(gettext("%s: Wrong size for input argument #%d: A polynomial expected.\n"),"routh_t",1))
         end
 
@@ -67,34 +68,42 @@ function r=routh_t(h,k,normalized)
     //
     r1=cobf(nd:-2:1);
     r2=cobf(nd-1:-2:1);
-    ncol=size(r1,"*");
-    if size(r2,"*")<>ncol then r2=[r2,0],end
+    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,
-            if flag==0 then
-                r(i,1:ncol-1)=[r(i-1,1),-r(i-2,1)]*[r(i-2,2:ncol);r(i-1,2:ncol)]
-            else
-                r(i,1:ncol-1)=[1.,-r(i-2,1)/r(i-1,1)]*[r(i-2,2:ncol);r(i-1,2:ncol)]
+            if flag==0 then // for a non-normalized table
+                r(i,1:ncol-1)=[r(i-1,1),-r(i-2,1)]*[r(i-2,2:ncol);r(i-1,2:ncol)] 
+            else // for a normalized table
+                if r(i-1,1)==0 then
+                    if type(k)<>1 then
+                        error(msprintf(gettext("%s: Wrong type for input argument #%d: An scalar expected.\n"),"routh_t",2));
+                    end
+                    r(i-1,1)=poly(0,'eps') 
+                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;
     else
         for i=3:nd,
+            // Special Case: Row of zeros detected: 
             if and(r(i-1,:)==0) then
                 naux=nd-i+2 //order of previous polynomial
                 exponents=naux:-2:0
-                ncoeff=size(exponents,"*")
+                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
+            // Special Case: First element of the 2nd row or upper is zero and is replaced with %eps:
+            if r(i-1,1)==0 then 
                 if rhs==1 then
-                    if typeof(r)=="rational" then
+                    if typeof(r)=='rational' then 
                         //scilab is not able to handle multivariable polynomials
-                        r=horner(r,%eps^2);
+                        r=horner(r,%eps^2); 
                     end
-                    r(i-1,1)=poly(0,"eps")
+                    r(i-1,1)=poly(0,'eps')
                 else
                     r(i-1,1)=%eps^2,
                 end
@@ -103,4 +112,25 @@ function r=routh_t(h,k,normalized)
 
         end;
     end;
+
+    if lhs==2  then
+        if rhs==2 & type(k)<>1 then
+            error(msprintf(gettext("%s: Wrong type for input argument #%d: A scalar expected.\n"),"routh_t",2));
+        end
+        nrow=size(r,'r')
+        num = 0;
+        c = 0;
+        for i = 1:nrow
+            if horner(r(i,1),%eps) >= 0 then
+                if c == 1 then
+                    num= num+1;
+                    c = 0;
+                end
+            elseif c == 0 then
+                num= num+1;
+                c = 1;
+            end;
+        end;
+    end
+
 endfunction
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_12828.dia.ref b/scilab/modules/cacsd/tests/nonreg_tests/bug_12828.dia.ref
new file mode 100644 (file)
index 0000000..e2629da
--- /dev/null
@@ -0,0 +1,19 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- Non-regression test for bug 12828 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12828
+//
+// <-- Short Description -->
+// routh_t() function with two input argument (routh_t(h,kp)) gives a wrong table for the special case
+// in which the first element of a row is zero.
+s=%s;
+P=s^3+3*s;
+kp=2;
+routh_t(1/P, kp);
+assert_checkequal(routh_t(1/P, kp), routh_t(P+kp));
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_12828.tst b/scilab/modules/cacsd/tests/nonreg_tests/bug_12828.tst
new file mode 100644 (file)
index 0000000..8da7c45
--- /dev/null
@@ -0,0 +1,20 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- Non-regression test for bug 12828 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12828
+//
+// <-- Short Description -->
+// routh_t() function with two input argument (routh_t(h,kp)) gives a wrong table for the special case
+// in which the first element of a row is zero.
+
+s=%s;
+P=s^3+3*s;
+kp=2;
+routh_t(1/P, kp);
+assert_checkequal(routh_t(1/P, kp), routh_t(P+kp));
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_12829.dia.ref b/scilab/modules/cacsd/tests/nonreg_tests/bug_12829.dia.ref
new file mode 100644 (file)
index 0000000..75085a4
--- /dev/null
@@ -0,0 +1,19 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- Non-regression test for bug 12829 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12829
+//
+// <-- Short Description -->
+// routh_t() function with two input argument (routh_t(h,kp)) gives a wrong table for the special case
+// in which the first element of a row is zero.
+s=%s;
+P=5*s^3-10*s^2+7*s+20;
+routh_t((1+s)/P,10);
+[r, num]=routh_t(P);
+assert_checkequal(num, 2);
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_12829.tst b/scilab/modules/cacsd/tests/nonreg_tests/bug_12829.tst
new file mode 100644 (file)
index 0000000..8f543e5
--- /dev/null
@@ -0,0 +1,20 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- Non-regression test for bug 12829 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=12829
+//
+// <-- Short Description -->
+// routh_t() function with two input argument (routh_t(h,kp)) gives a wrong table for the special case
+// in which the first element of a row is zero.
+
+s=%s;
+P=5*s^3-10*s^2+7*s+20;
+routh_t((1+s)/P,10);
+[r, num]=routh_t(P);
+assert_checkequal(num, 2);
index a8e259d..c629767 100644 (file)
@@ -11,7 +11,10 @@ s=poly(0,'s');
 k=poly(0,'k');
 //Test with a polynomial
 P1=s*(s+7)*(s+11);
-assert_checkequal(routh_t(P1),[1 77; 18 0; 77 0; 0 0]);
+assert_checkequal(routh_t(P1), [1 77; 18 0; 77 0; 0 0]);
+[r, num] = routh_t(P1);
+assert_checkequal(r, [1 77; 18 0; 77 0; 0 0]);
+assert_checkequal(num, 0);
 //Test with a transfert function
 h=1/P1;
 ref1=[1 77; 18 k;(1386-k)/18 0; k 0];
@@ -24,10 +27,16 @@ assert_checkequal(numer(routh_t(h,k,%t))/horner(denom(routh_t(h,k,%t)),1),numref
 assert_checkequal(routh_t(h,k,%f),[1 77; 18 k; 1386-k 0; 1386*k-k*k 0]);
 // One of the coefficients in the polynomial equals zero
 P2=2*s^2-24;
-assert_checkequal(routh_t(P2),[2 -24; 4 0; -24 0]);
+assert_checkequal(routh_t(P2), [2 -24; 4 0; -24 0]);
+[r, num] = routh_t(P2);
+assert_checkequal(r, [2 -24; 4 0; -24 0]);
+assert_checkequal(num, 1);
 // A row full of zeros
 P3=s^4-6*s^3+10*s^2-6*s+9;
 assert_checkequal(routh_t(P3),[1 10 9; -6 -6 0; 9 9 0; 18 0 0; 9 0 0]);
+[r, num] = routh_t(P3);
+assert_checkequal(r,[1 10 9; -6 -6 0; 9 9 0; 18 0 0; 9 0 0]);
+assert_checkequal(num, 2);
 //Error messages
 assert_checkfalse(execstr("routh_t(P,%t)","errcatch")==0);
 assert_checkfalse(execstr("routh_t(P,k)","errcatch")==0);
index f9b96b8..427d08b 100644 (file)
@@ -15,7 +15,10 @@ k=poly(0,'k');
 
 //Test with a polynomial
 P1=s*(s+7)*(s+11);
-assert_checkequal(routh_t(P1),[1 77; 18 0; 77 0; 0 0]);
+assert_checkequal(routh_t(P1), [1 77; 18 0; 77 0; 0 0]);
+[r, num] = routh_t(P1);
+assert_checkequal(r, [1 77; 18 0; 77 0; 0 0]);
+assert_checkequal(num, 0);
 
 //Test with a transfert function
 h=1/P1;
@@ -30,11 +33,17 @@ assert_checkequal(routh_t(h,k,%f),[1 77; 18 k; 1386-k 0; 1386*k-k*k 0]);
 
 // One of the coefficients in the polynomial equals zero
 P2=2*s^2-24;
-assert_checkequal(routh_t(P2),[2 -24; 4 0; -24 0]);
+assert_checkequal(routh_t(P2), [2 -24; 4 0; -24 0]);
+[r, num] = routh_t(P2);
+assert_checkequal(r, [2 -24; 4 0; -24 0]);
+assert_checkequal(num, 1);
 
 // A row full of zeros
 P3=s^4-6*s^3+10*s^2-6*s+9;
 assert_checkequal(routh_t(P3),[1 10 9; -6 -6 0; 9 9 0; 18 0 0; 9 0 0]);
+[r, num] = routh_t(P3);
+assert_checkequal(r,[1 10 9; -6 -6 0; 9 9 0; 18 0 0; 9 0 0]);
+assert_checkequal(num, 2);
 
 //Error messages
 assert_checkfalse(execstr("routh_t(P,%t)","errcatch")==0);