bug 6829 fixed
Serge Steer [Mon, 29 Mar 2010 13:31:20 +0000 (15:31 +0200)]
scilab/CHANGES_5.2.X
scilab/modules/cacsd/help/en_US/kpure.xml
scilab/modules/cacsd/macros/kpure.sci
scilab/modules/cacsd/tests/nonreg_tests/bug_6829.dia.ref [new file with mode: 0644]
scilab/modules/cacsd/tests/nonreg_tests/bug_6829.tst [new file with mode: 0644]

index cb712b8..98687b8 100644 (file)
@@ -182,6 +182,8 @@ Bug Fixes:
 
 * bug 6817 fixed - The "xstring" function shifted the start position of the string relative to the string length.
 
+* bug 6829 fixed - kpure fails to compute  when applied to an high degree system
+
 
                 Changes between Versions 5.2.0 and 5.2.1 of Scilab
                 ==================================================
index fa74ae2..b7cae81 100644 (file)
       <varlistentry>
         <term>tol</term>
         <listitem>
-          <para>vector with 2 elements <literal>[epsK epsI]</literal>. <literal>epsK</literal> is a tolerance used to determine if two values of  <literal>K</literal> can be considered as  equal <literal>epsI</literal> is a tolerance used to determine if a root is imaginary or not. The default value is <literal>[1e-6 1e-6 ]</literal> </para>
+          <para>a positive scalar.  tolerance used to determine if a
+          root is imaginary or not. The default value is
+          <literal>1e-6</literal> </para>
         </listitem>
       </varlistentry>
       <varlistentry>
         <term>K</term>
         <listitem>
-          <para>Real vector, the vector of gains for which at least one closed loop pole is imaginary. </para>
+          <para>Real vector, the vector of gains for which at least
+          one closed loop pole is imaginary. </para>
         </listitem>
       </varlistentry>
       <varlistentry>
         <term>R</term>
         <listitem>
-          <para>Complex vector, the imaginary closed loop poles associated with the  values of  <literal>K</literal>.</para>
+          <para>Complex vector, the imaginary closed loop poles
+          associated with the values of <literal>K</literal>.</para>
         </listitem>
       </varlistentry>
     </variablelist>
   <refsection>
     <title>Examples</title>
     <programlisting role="example"><![CDATA[ 
-s=poly(0,'s');
-h=syslin('c',(s-1)/(1+5*s+s^2+s^3))
+
+num=real(poly([-1+%i, -1-%i, -1+8*%i  -1-8*%i],'s'));
+den=real(poly([0.5 0.5  -6+7*%i  -6-7*%i  -3 -7 -11],'s'));
+h=num/den;
+
+[K,Y]=kpure(h)
 clf();evans(h)
-K=kpure(h)
-hf=h/.K(1)
-roots(denom(hf))
+plot(real(Y),imag(Y),'+r')
+
+num=real(poly([-1+%i*1, -1-%i*1, 2+%i*8  2-%i*8 -2.5+%i*13 -2.5-%i*13],'s'));
+den=real(poly([1 1 3+%i*3 3-%i*3 -15+%i*7  -15-%i*7  -3 -7 -11],'s'));
+h=num/den;
+
+[K,Y]=kpure(h)
+clf();evans(h,100000)
+plot(real(Y),imag(Y),'+r')
+
+
+
  ]]></programlisting>
   </refsection>
   <refsection>
+    <title>Author</title>
+    <para>Serge Steer, INRIA</para>
+  </refsection>
+  <refsection>
     <title>See Also</title>
     <simplelist type="inline">
       <member>
index 2f1d4c8..5df9503 100644 (file)
@@ -7,59 +7,43 @@
 // are also available at    
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
-function [y,R]=kpure(sl,eps)
+function [K,R]=kpure(sl,eps)
+//if sl is a transfert function N(S)/D(s) kpure looks for K producing
+//pure  imaginary roots for
+//  P(s)=D(s)+K*N(s)
+//There is a pair of pure imaginary poles if and only if
+//  P(%i*q)=0  (1)
+// and
+//  P(-%i*q)=0 (2)
+// because N and D are polynomials with real coefficients.
+  
+//Author: Serge Steer, INRIA
   y=[];R=[];
-  if argn(2)==1 then eps=[1e-6,1e-6],end
+  msg=_("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n")
+  if argn(2)==1 then eps=1e-6,end
+  if size(eps,'*')==2 then  eps=eps(2),end //compatibility
   select typeof(sl)
   case 'rational' then
-    if size(sl.num,'*') > 1 then error(95,1),end
+    if size(sl.num,'*') <> 1 then 
+      error(msprintf(msg,"kpure",1))
+    end
   case 'state-space' then
+    if size(sl.D,'*') <> 1 then 
+      error(msprintf(msg,"kpure",1))
+    end
     sl=ss2tf(sl)
-    if size(sl.num,'*') > 1 then error(95,1),end
   else
-    error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"kpure",1))
+    error(msprintf(msg,"kpure",1))
   end
-  if sl.dt<>'c' then 
-    error(msprintf(gettext("%s: Wrong value for input argument #%d: Continuous time system expected.\n"),"kpure",1))
-  end
-  if size(sl.num,'*')<>1 then
-    error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"kpure",1))
-  end
-
-  //build the Routh table of the given system
-  r=routh_t(sl,poly(0,'k')),
-  [s,t]=size(r);
-  
-  //Check of infinite solution
-  for i=1:s,
-    if and(coeff(r(i,:))==0) then 
-      error(msprintf(gettext("%s: Infinite solution.\n"),"kpure")),
-    end
-  end,
-
-  z=0;u=0;
-
-  for i=1:s,
-    k=roots(gcd(r(i,:)));// k is must be real
-    k=real(k)
-    //k(k<=0)=[]; //remove negative and zero values
-    h=prod(size(k)),
-    if h>0 then
-      for j=1:h,
-       if and(abs(k(j)-y)>eps(1)) then //remove duplicated k's
-         //compute the poles of the closed loop
-         wr=roots(k(j)*sl.num+sl.den) 
-         //retains k if at least one closed loop pole is imaginary            
-         [w,ki]=mini(abs(real(wr)))
-         if w<eps(2) then 
-           y=[y real(k(j))],
-           R=[R wr(ki)]
-         end
-       end
-      end,
-    end;
-  end;
-
-  [y,k]=gsort(y)
-  R=R(k)
+  //(1) give K(s)=-D(s)/N(s)
+  s=poly(0,varn(sl.den))
+  K=-sl.den/sl.num;
+  // replace K by the previous value in (2) and find the roots
+  s=roots(numer(horner(sl.den,-s)+K*horner(sl.num,-s)),'e');
+  //retain pure imaginary roots
+  s=imag(s(abs(real(s))<eps));
+  R=(s(s>0).'*%i);
+  //find the K(s) values K(s)=-D(s)/N(s)
+  K=-real(freq(den,num,R))
 endfunction
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_6829.dia.ref b/scilab/modules/cacsd/tests/nonreg_tests/bug_6829.dia.ref
new file mode 100644 (file)
index 0000000..2616fdd
--- /dev/null
@@ -0,0 +1,26 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2009 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- Non-regression test for bug 6829 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=6829
+//
+// <-- Short Description -->
+//  kpure fails to compute  gains when applied to an high degree system
+num=real(poly([-1+%i*1, -1-%i*1, 2+%i*8  2-%i*8 -2.5+%i*13 -2.5-%i*13],'s'));
+den=real(poly([1 1 3+%i*3 3-%i*3 -15+%i*7  -15-%i*7  -3 -7 -11],'s'));
+h=num/den;
+[K,Y]=kpure(h);
+n=size(K,'*');
+if n<>4 then bugmes();quit;end
+for i=1:n
+  r=roots(denom(h/.K(i)));
+  r=r(abs(real(r))<1e-6);//pure imaginary
+  r=r(imag(r)>0); //retains only positive imaginary part
+  if r==[] then bugmes();quit;end
+  if abs(r-Y(i))>1e-10 then bugmes();quit;end
+end
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_6829.tst b/scilab/modules/cacsd/tests/nonreg_tests/bug_6829.tst
new file mode 100644 (file)
index 0000000..d69cffc
--- /dev/null
@@ -0,0 +1,29 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2009 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- Non-regression test for bug 6829 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=6829
+//
+// <-- Short Description -->
+//  kpure fails to compute  gains when applied to an high degree system
+
+num=real(poly([-1+%i*1, -1-%i*1, 2+%i*8  2-%i*8 -2.5+%i*13 -2.5-%i*13],'s'));
+den=real(poly([1 1 3+%i*3 3-%i*3 -15+%i*7  -15-%i*7  -3 -7 -11],'s'));
+h=num/den;
+
+[K,Y]=kpure(h);
+n=size(K,'*');
+if n<>4 then pause,end
+for i=1:n
+  r=roots(denom(h/.K(i)));
+  r=r(abs(real(r))<1e-6);//pure imaginary
+  r=r(imag(r)>0); //retains only positive imaginary part
+  if r==[] then pause,end
+  if abs(r-Y(i))>1e-10 then pause,end
+end