bug 9285 fix 17/4117/2
Serge Steer [Wed, 25 May 2011 09:13:45 +0000 (11:13 +0200)]
Change-Id: I6fdc6ecb98816b69a9912554ba508564079fe9aa

scilab/CHANGES_5.3.X
scilab/modules/cacsd/macros/g_margin.sci
scilab/modules/cacsd/tests/nonreg_tests/bug_9285.dia.ref [new file with mode: 0644]
scilab/modules/cacsd/tests/nonreg_tests/bug_9285.tst [new file with mode: 0644]

index c3167ee..c0e858b 100644 (file)
@@ -270,6 +270,8 @@ CACSD Module:
                    the symetric part and a new function added for high
                    definition zoom.
 
+* Bug 9285 fixed - g_margin() returned an erroneous result.
+
 
 Bug Fixes:
 ==========
index 13d9bc0..b68ce88 100644 (file)
@@ -20,13 +20,21 @@ function [gm,fr]=g_margin(h)
     error(msprintf(_("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"g_margin",1))
   end
   //
-  eps=1.e-7;// threshold used for testing if complex numbers are real or pure imaginary
+  epsr=1.e-7;//used for testing if complex numbers are real
+  eps1=1.e-7;//used for testing if complex numbers have a modulus near 1
+  epssing=1e-10; //used for testing if arguments are not singular points of h
   if h.dt=='c' then  //continuous time case
     // get s such as h(s)=h(-s) and s=iw 
      s=%i*poly(0,"w");
-     w=roots(imag(horner(h.num,s)*conj(horner(h.den,s))),"e")
-     ws=real(w(abs(imag(w))<eps&real(w)<=0)) //points where phase is -180°
-     ws(abs(horner(h.den,%i*ws))==0)=[];
+     //compute h(s)-h(-s)=num/den
+     num=imag(horner(h.num,s)*conj(horner(h.den,s)))
+     den=real(horner(h.den,s)*conj(horner(h.den,s)))
+     //necessary condition
+     w=roots(num,"e");
+     ws=real(w(abs(imag(w))<epsr&real(w)<=0)) //points where phase is -180°
+                                             
+     //remove nearly singular points     
+     ws(abs(horner(num,ws))>=epssing*abs(horner(den,ws)))=[]
      if ws==[] then gm=%inf,fr=[],return,end
      mingain=real(freq(h.num,h.den,%i*ws))
   else  //discrete time case
@@ -37,12 +45,17 @@ function [gm,fr]=g_margin(h)
     sm=simp_mode();simp_mode(%f);hh=h-horner(h,1/z);simp_mode(sm)
     //find the numerator roots
     z=roots(hh.num,"e");
-    z(abs(abs(z)-1)>eps)=[]// retain only roots with modulus equal to 1
+    z(abs(abs(z)-1)>eps1)=[]// retain only roots with modulus equal to 1
+    
+    //remove nearly singular points                                       
+    z(abs(horner(hh.num,z))>=epssing*abs(horner(hh.den,z)))=[];
+
     w=log(z)/(%i*dt)
-    ws=real(w(abs(imag(w))<eps)) //points where phase is -180°
+    ws=real(w(abs(imag(w))<epsr)) //points where phase is -180°
     if ws==[] then gm=%inf,fr=[],return,end
     mingain=real(horner(h,exp(%i*ws*dt)))
   end
+  
   k=find(mingain<0)
   if k==[] then gm=%inf,fr=[],return,end
   mingain=abs(mingain(k));
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_9285.dia.ref b/scilab/modules/cacsd/tests/nonreg_tests/bug_9285.dia.ref
new file mode 100644 (file)
index 0000000..06f3979
--- /dev/null
@@ -0,0 +1,63 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- JVM NOT MANDATORY -->
+//
+// <-- Non-regression test for bug 9285 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=9285 
+//
+// <-- Short Description -->
+// g_margin returned bad results in some special cases.
+s=%s;
+Kp=2.054;
+Ki=0.13383297644539613;
+Td=0.605;
+Kd=5.0;
+ts=0.25;
+k=0.001712328003555509;
+tw=1.4731525037117905;
+ta=8.678149822856799;
+au=0.0;
+ty=0.1;
+Process=(k*s^2-tw*s+1)/((k*s^2+(tw/2)*s+1)*(1+ts*s)^2*(ta*s+au))
+ Process  =
+                                                2                    
+                     1 - 1.4731525s + 0.0017123s                     
+    -------------------------------------------------------------    
+                           2            3            4            5  
+    8.6781498s + 10.731194s + 3.7533037s + 0.4069374s + 0.0009287s   
+FeedBack=1/((ty*ty)*s^2+(ty*2^0.5)*s+1) * (1 + (Kd*Td*s)/(Kd+Td*s))
+ FeedBack  =
+                  5 + 3.63s                  
+    -------------------------------------    
+                               2          3  
+    5 + 1.3121068s + 0.1355599s + 0.00605s   
+REG=(Kp * ( 1 + Ki/s));
+FT20_0_ = syslin('c', REG * Process * FeedBack);
+[gm,fr]=g_margin(FT20_0_)
+ fr  =
+    0.1173008  
+ gm  =
+    6.6932617  
+if ~assert_checkalmostequal(fr,0.11730083,0,6e-9) then bugmes();quit;end
+if ~assert_checkalmostequal(gm,6.69326173,0,6e-9) then bugmes();quit;end
+FT20_0 = tf2ss(FT20_0_);
+[gm,fr]=g_margin(FT20_0)
+ fr  =
+    0.1173008  
+ gm  =
+    6.6932617  
+if ~assert_checkalmostequal(fr,0.11730083,0,6e-9) then bugmes();quit;end
+if ~assert_checkalmostequal(gm,6.69326173,0,6e-9) then bugmes();quit;end
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_9285.tst b/scilab/modules/cacsd/tests/nonreg_tests/bug_9285.tst
new file mode 100644 (file)
index 0000000..f164a2b
--- /dev/null
@@ -0,0 +1,41 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - INRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- JVM NOT MANDATORY -->
+//
+// <-- Non-regression test for bug 9285 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=9285 
+//
+// <-- Short Description -->
+// g_margin returned bad results in some special cases.
+s=%s;
+Kp=2.054;
+Ki=0.13383297644539613;
+Td=0.605;
+Kd=5.0;
+ts=0.25;
+k=0.001712328003555509;
+tw=1.4731525037117905;
+ta=8.678149822856799;
+au=0.0;
+ty=0.1;
+
+Process=(k*s^2-tw*s+1)/((k*s^2+(tw/2)*s+1)*(1+ts*s)^2*(ta*s+au))
+FeedBack=1/((ty*ty)*s^2+(ty*2^0.5)*s+1) * (1 + (Kd*Td*s)/(Kd+Td*s))
+REG=(Kp * ( 1 + Ki/s));
+
+FT20_0_ = syslin('c', REG * Process * FeedBack);
+[gm,fr]=g_margin(FT20_0_)
+if ~assert_checkalmostequal(fr,0.11730083,0,6e-9) then pause,end
+if ~assert_checkalmostequal(gm,6.69326173,0,6e-9) then pause,end
+
+FT20_0 = tf2ss(FT20_0_);
+[gm,fr]=g_margin(FT20_0)
+if ~assert_checkalmostequal(fr,0.11730083,0,6e-9) then pause,end
+if ~assert_checkalmostequal(gm,6.69326173,0,6e-9) then pause,end