* Bug 15368 fixed: freson() issues on some cases 84/21184/5
Samuel GOUGEON [Sat, 21 Dec 2019 22:19:42 +0000 (23:19 +0100)]
  http://bugzilla.scilab.org/15368

Change-Id: Ic6385968a3c71aa627fdc68ae237c4424bb95abd

scilab/CHANGES.md
scilab/modules/cacsd/macros/freson.sci
scilab/modules/cacsd/tests/nonreg_tests/bug_15368.dia.ref [new file with mode: 0644]
scilab/modules/cacsd/tests/nonreg_tests/bug_15368.tst

index 19edb2c..887fe3d 100644 (file)
@@ -187,6 +187,7 @@ Bug Fixes
 * [#15269](http://bugzilla.scilab.org/show_bug.cgi?id=15269): `xgetech` was poor and stiff compared to any combination of `gca()` properties `.axes_bounds`, `.data_bounds`, `.log_flags`, and `.margins`. It is removed.
 * [#15271](http://bugzilla.scilab.org/show_bug.cgi?id=15271): `bitget` needed to be upgraded.
 * [#15321](http://bugzilla.scilab.org/show_bug.cgi?id=15321): `lu()` was leaking memory.
+* [#15368](http://bugzilla.scilab.org/show_bug.cgi?id=15368): `freson()` silently returned frequencies not corresponding to a maximum, or returned [] instead of some still computable maxima frequencies.
 * [#15425](http://bugzilla.scilab.org/show_bug.cgi?id=15425): The Kronecker product `a.*.b` failed when `a` or `b` or both are hypermatrices, with one or both being polynomials or rationals.
 * [#15523](http://bugzilla.scilab.org/show_bug.cgi?id=15523): `%ODEOPTIONS(1)=2` didn't work with solvers 'rk' and 'rkf'
 * [#15248](http://bugzilla.scilab.org/show_bug.cgi?id=15248): `lsq()`was leaking memory.
index adff61a..9140372 100644 (file)
@@ -29,13 +29,13 @@ function fr = freson(h)
         msg = gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n")
         error(msprintf(msg, "freson", 1))
     end
-
+    h0 = h
+    dt = h.dt
     if typeof(h)=="state-space" then
         h=ss2tf(h)
     elseif typeof(h)=="zpk" then
         h=zpk2tf(h)
     end
-    dt = h.dt
     [n,d] = h(["num","den"]);
     if type(n)==1 then
         n=poly(n,varn(d),"c");
@@ -48,23 +48,35 @@ function fr = freson(h)
     //look for  omega such that derivative of magn. is zero
     if dt=="c" then
         //find frequencies which zeros the magnitude derivative
-        hh = h*horner(h,-%s)
+        hh = h*horner(h, -%s);
+        hh = clean(hh, 0, %eps);    // http://bugzilla.scilab.org/15368
         r = roots(derivat(hh).num)
-        // k=find(imag(r)>0&abs(real(r))<%eps*abs(r)); // http://bugzilla.scilab.org/15368
-        k = find(imag(r)>0 & abs(real(r))<1.e-14*abs(r));
+        k = find(imag(r)>0 & abs(real(r)) < %eps*abs(r));
         fr = imag(r(k))/(2*%pi)
     else
         if dt=="d" | dt==[] then dt = 1; end
         //find frequencies which zeros the magnitude derivative
-        hh = h*horner(h,1/%z)
-        r = roots(derivat(hh).num);
-        k = find(abs(abs(r)-1)<=sqrt(%eps)*abs(r));
-        r = imag(log(r(k)));
-        fr = r(r>0&r<0.999*%pi)/(2*%pi*dt)
+        hh = h*horner(h, 1/%z)
+        r = roots(derivat(hh).num)
+        k = find(abs(abs(r)-1) <= %eps^(1/size(r,1))*abs(r))
+        r = imag(log(r(k)))
+        fr = r(r>0 & r < 0.999*%pi) / (2*%pi*dt)
     end
     if fr==[] then return;end
 
     //find frequencies that correspond to a magnitude peak and sort them:
-    k = find(abs(repfreq(h,fr))-abs(repfreq(h,fr*0.999))>0)
+    k = find(abs(repfreq(h,fr)) - abs(repfreq(h,fr*0.999))>0)
     fr = gsort(fr(k),"g","d");
+
+    // Checking                 // http://bugzilla.scilab.org/15368
+    f = fr(:)'
+    f = [f*0.99 ; f ; f*1.01];
+    [f,repf] = repfreq(h0, f(:));
+    [?,mag] = phasemag(repf);
+    mag = matrix(mag, 3, -1);
+    i = find(mag(2,:)>mag(1,:) & mag(3,:)<mag(2,:))
+    if length(i) < length(fr) then
+        fr = fr(i)
+        warning(_("freson: Some maxima frequencies may be undetected."))
+    end
 endfunction
diff --git a/scilab/modules/cacsd/tests/nonreg_tests/bug_15368.dia.ref b/scilab/modules/cacsd/tests/nonreg_tests/bug_15368.dia.ref
new file mode 100644 (file)
index 0000000..109be11
--- /dev/null
@@ -0,0 +1,46 @@
+// =============================================================================\r
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab\r
+// Copyright (C) 2018, 2019 - Lucien POVY\r
+// Copyright (C) 2018, 2019 - Samuel GOUGEON\r
+//\r
+//  This file is distributed under the same license as the Scilab package.\r
+// =============================================================================\r
+//\r
+// <-- CLI SHELL MODE -->\r
+// <-- ENGLISH IMPOSED -->\r
+//\r
+// <-- Non-regression test for bug 15368 -->\r
+//\r
+// <-- Bugzilla URL -->\r
+// http://bugzilla.scilab.org/15368\r
+//\r
+// <-- Short Description -->\r
+// For some continuous time systems, freson() wrongly returned []\r
+s = %s;\r
+rtol = 1e-14;\r
+sl1 = syslin("c",1,s*(s+1)*(1+s/3));\r
+[K,P] = kpure(sl1);\r
+sl = K*sl1;\r
+taud = 1/imag(P);\r
+sltaud = sl*(1+taud*s);         //first open loop\r
+sltaud2 = sl*(1+2*taud*s);      //second open loop\r
+ff = freson(sltaud/(1+sltaud)); //freson of first close loop\r
+assert_checkalmostequal(ff, 0.32447218051387, rtol);\r
+SLBF = sltaud2/(1+sltaud2); //THE CLOSE LOOP OF sltau2\r
+fr6 = freson(SLBF);\r
+assert_checkalmostequal(fr6, 0.48257117694410, rtol);\r
+// ----------------------\r
+num = s + 1.4*s^2 + 1.4*s^3 + 0.4*s^4 ;\r
+den = 0.5 + 0.8*s + 1.4*s^2 + 1.4*s^3 + 0.4*s^4 ;\r
+h = syslin("c", num, den);\r
+assert_checkalmostequal(freson(h), 0.1111642261457733, rtol);\r
+// ----------------------\r
+h = syslin('c',-1+%s,(3+2*%s+%s^2)*(50+0.1*%s+%s^2));\r
+hd = dscr(h, 0.05);\r
+assert_checkalmostequal(freson(hd), [1.125279232661004 ; 0.2615570266570646]);\r
+hd = dscr(h, 0.02);\r
+assert_checkalmostequal(freson(hd), 1.123778, 1e-4);\r
+WARNING: freson: Some maxima frequencies may be undetected.\r
+hd = dscr(h, 0.01);\r
+assert_checkequal(freson(hd), []);\r
+WARNING: freson: Some maxima frequencies may be undetected.\r
index 3cd4b65..8d04694 100644 (file)
@@ -1,13 +1,13 @@
 // =============================================================================
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2018 - Lucien POVY
-// Copyright (C) 2018 - Samuel GOUGEON
+// Copyright (C) 2018, 2019 - Lucien POVY
+// Copyright (C) 2018, 2019 - Samuel GOUGEON
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
 //
 // <-- CLI SHELL MODE -->
-// <-- NO CHECK REF -->
+// <-- ENGLISH IMPOSED -->
 //
 // <-- Non-regression test for bug 15368 -->
 //
 // For some continuous time systems, freson() wrongly returned []
 
 s = %s;
-sl1 = syslin("c",1,s*(s+1)*(1+s/3))
-[K,P] = kpure(sl1)
-sl = K*sl1
-taud = 1/imag(P)
-sltaud=sl*(1+taud*s)//first open loop
-sltaud2=sl*(1+2*taud*s)//second open loop
-ff = freson(sltaud/(1+sltaud))//freson of first close loop
-assert_checkalmostequal(ff, 0.32447218051387, 1e-14);
+rtol = 1e-14;
+sl1 = syslin("c",1,s*(s+1)*(1+s/3));
+[K,P] = kpure(sl1);
+sl = K*sl1;
+taud = 1/imag(P);
+sltaud = sl*(1+taud*s);         //first open loop
+sltaud2 = sl*(1+2*taud*s);      //second open loop
+ff = freson(sltaud/(1+sltaud)); //freson of first close loop
+assert_checkalmostequal(ff, 0.32447218051387, rtol);
 
 SLBF = sltaud2/(1+sltaud2); //THE CLOSE LOOP OF sltau2
 fr6 = freson(SLBF);
-assert_checkalmostequal(fr6, 0.48257117694410, 1e-14);
+assert_checkalmostequal(fr6, 0.48257117694410, rtol);
 
+// ----------------------
+
+num = s + 1.4*s^2 + 1.4*s^3 + 0.4*s^4 ;
+den = 0.5 + 0.8*s + 1.4*s^2 + 1.4*s^3 + 0.4*s^4 ;
+h = syslin("c", num, den);
+assert_checkalmostequal(freson(h), 0.1111642261457733, rtol);
+
+// ----------------------
+
+h = syslin('c',-1+%s,(3+2*%s+%s^2)*(50+0.1*%s+%s^2));
+
+hd = dscr(h, 0.05);
+assert_checkalmostequal(freson(hd), [1.125279232661004 ; 0.2615570266570646]);
+
+hd = dscr(h, 0.02);
+assert_checkalmostequal(freson(hd), 1.123778, 1e-4);
+
+hd = dscr(h, 0.01);
+assert_checkequal(freson(hd), []);