Bub in module CACSD macro freson
Reported by Lucien Povy
BUG DESCRIPTION:Bug in cacsd module, macro "freson.sci", for some continuous system freson returns [].
----------------
ERROR LOG:
----------
HOW TO REPRODUCE THE BUG: Do :
s=%s, h=syslin("c",4+4.6188*s,s+1.333*s*s+0.3333*s^3)
fre=freson(h)
This version of freson returns [] Scilab-6..
The good value is fre=0.4825 Hz in scilab-5... and by black and datatips
-------------------------
OTHER INFORMATION:
The problem is :
At lines 47 and 54 we do hh=h*horner(h,-%s) or hh=h*horner(h,1/%z),
the result must be even rationals, but in the exemple it is not :
One solution is : replace hh=clean(h*horner(h,-%s)) or hh=clean(h*horner(h,1/%z))
------------------Activity
By Samuel GOUGEON (@sgougeon)
Hello Lucien,
This version of freson returns [] Scilab-6..
With Scilab 6.0.1 [2017-12-18], I get an error, not []:
--> fre=freson(h) at line 29 of function freson ( SCI\modules\cacsd\macros\freson.sci line 41 )
freson: Wrong value for input argument #1: infinite gain at zero frequency.
The good value is fre=0.4825 Hz in scilab-5... and by black and datatips
With 5.5.2, i get the same error (but not from the same line): -->fre=freson(h) !--error 10000 freson: Wrong value for input argument #1: infinite gain at zero frequency. at line 8 of function freson called by :
fre=freson(h)So, which version 5.x are you referring to/ on which you got the error?
-
By Samuel GOUGEON (@sgougeon)
The same error occurs with the 5.4.1, 5.3.0, 5.1.0, and even 4.1.2 scilab versions. So, it does not come from Scilab 6. It is an old bug.
OK, i have understood your point about black() and datatips.
-
By Samuel GOUGEON (@sgougeon)
At lines 47 and 54 we do hh=hhorner(h,-%s) or hh=hhorner(h,1/%z), the result must be even rationals, but in the exemple it is not : One solution is : replace hh=clean(hhorner(h,-%s)) or hh=clean(hhorner(h,1/%z))
This can't fix the error, since this one occurs upstream, on line 29.
By the way, the fact that you get [] while i get a proper error is still to be explained. Are you working on Linux? We have already found different behaviors on Windows and Linux, as for the #14162.
-
By Samuel GOUGEON (@sgougeon)
--> fre=freson(h) at line 29 of function freson ( SCI\modules\cacsd\macros\freson.sci line 41 ) freson: Wrong value for input argument #1: infinite gain at zero frequency.
This is true for the given system. So, in my opinion, there is no bug. According to its bode diagram, this system does not show any other magnitude peak. The only thing that we could possibly expect is that freson() returns 0 Hz instead of an error.
So, could you elaborate about
The good value is fre=0.4825 Hz in scilab-5... and by black and datatips
? Thanks Best regards Samuel
-
By Lucien Povy
hello, SAMUEL Excuse me i do a big mistake, "h" in the example, is the open loop model and the problem is : what is the resonance frequency of close loop ? The exercise was: s=%s;sl1=syslin("c",1,s*(s+1)(1+s/3)) [K,P]=kpure(sl1) sl=Ksl1 taud=1/imag(P) sltaud=sl*(1+tauds)//first open loop sltaud2=sl(1+2tauds)//second open loop ff=freson(sltaud/(1+sltaud))//freson of first close loop rep1=repfreq(sltaud/(1+sltaud),ff)//response at this frequency for first close loop [d1,phi1]=dbphi(rep1)//magnitude and phase //NO PROBLEM //BUT SLBF=sltaud2/(1+sltaud2)//THE CLOSE LOOP OF sltau2 fr6=freson(SLBF)//NO resonance frequency //RETURN fr6=[] //HERE IS THE BUG.......... //IF I DO CORRECTION OF freson (line 47 hh=clean(hhorner(h,-s)) freson --> ffreson rre=repfreq(SLBF,ffreson(SLBF))//NO PROBLEM [dbf,phif]=dbphi(rre) //ok fr=ffreson(SLBF)//ok //do black(SLBF,0.1,1)//ok with black and datatips figure(1);bode(SLBF,0.1,1)//ok now do: h=SLBF;hh = hhorner(h,-s)//hh IS NOT EVEN FUNCTION, terms - 2.842D-14s in hh.num and - 3.553D-15s^3 in hh.den. //I thing it's a problem of precision in h*horner(h,-s)
-
By Samuel GOUGEON (@sgougeon)
(In reply to Lucien Povy from comment #5)
hello, SAMUEL Excuse me i do a big mistake, "h" in the example, is the open loop model and the problem is : what is the resonance frequency of close loop ? The exercise was: s=%s;sl1=syslin("c",1,s*(s+1)(1+s/3)) [K,P]=kpure(sl1) sl=Ksl1 taud=1/imag(P) sltaud=sl*(1+tauds)//first open loop sltaud2=sl(1+2tauds)//second open loop ff=freson(sltaud/(1+sltaud))//freson of first close loop rep1=repfreq(sltaud/(1+sltaud),ff)//response at this frequency for first close loop [d1,phi1]=dbphi(rep1)//magnitude and phase //NO PROBLEM
really? I get:
--> sltaud/(1+sltaud) ans =
12 + 6.9282032s
2 312 + 9.9282032s + 4s + s
--> ff=freson(sltaud/(1+sltaud))//freson of first close loop ff = [] <<<<<<<< The peak is not detected
//IF I DO CORRECTION OF freson (line 47 hh=clean(hhorner(h,-s)) freson --> ffreson rre=repfreq(SLBF,ffreson(SLBF))//NO PROBLEM [dbf,phif]=dbphi(rre) //ok fr=ffreson(SLBF)//ok //do black(SLBF,0.1,1)//ok with black and datatips figure(1);bode(SLBF,0.1,1)//ok now do: h=SLBF;hh = hhorner(h,-s)//hh IS NOT EVEN FUNCTION, terms - 2.842D-14s in hh.num and - 3.553D-15s^3 in hh.den. //I thing it's a problem of precision in h*horner(h,-s)
Of numerical residues? It is often quite tricky to use clean(), that may have also some bad side effects. To be tested and confirmed.
-
By Samuel GOUGEON (@sgougeon)
Fix proposed on review for Scilab 6.0.1 @ https://gitlab.com/scilab/legacy_codereview/-/wikis/changes/19730/
-
By Samuel GOUGEON (@sgougeon)
(In reply to Samuel GOUGEON from comment #7)
Fix proposed on review for Scilab 6.0.1 @ https://gitlab.com/scilab/legacy_codereview/-/wikis/changes/19730/
The new relative tolerance 1e-14 instead of %eps is proposed by Lucien after some trials on a set of problematic systems (as discussed in private mails). As for numerical peak detection based on derivative, it is not surprising to have to not be too much severe in tolerance, since derivatives are strongly noise-dependent.
-
By Antoine ELIAS (@antoine.elias)
merged @ be95d4dc
-
By Lucien Povy
Le 25/09/2019 I believe that there are new problems with "freson.sci" program in CACSD module. If you are running two programs, one for continuous system another for sampled system, sometimes you have a result as [ ] : you must improve the accuracy of result ;("ffreson.sci" program is proposed). (Bugs not bugs ???? yes bugs : see in text REMARK.
1 . Example : continuous system. s=%s ; num = s+1.4s^2+1.4s^3+0.4s^4 ; den = 0.5+0.8s+1.4s^2+1.4s^3+0.4*s^4 ; h = syslin("c",num,den) fr = freson(h) fr = [] gainplot(h,0.01,1) with my new function "ffreson.sci" : ffreson(h)//La nouvelle fonction. ans = 0.111164226146
2 . Example : sampled system. The exercise in "help" page of "freson" with sampling period Ts=0.01s : the program gives "fr=[]" . With sampling period Ts=0.04s you have a good result (problem of accuracy) : run the example in help and change sampling period. More explanations in my email "New problems..."on Scilab-users".
-
By Lucien Povy
Le 25/09/2019 I believe that there are new problems with "freson.sci" program in CACSD module. If you are running two programs, one for continuous system another for sampled system, sometimes you have a result as [ ] : you must improve the accuracy of result ;("ffreson.sci" program is proposed). (Bugs not bugs ???? yes bugs : see in text REMARK.
1 . Example : continuous system. s=%s ; num = s+1.4s^2+1.4s^3+0.4s^4 ; den = 0.5+0.8s+1.4s^2+1.4s^3+0.4*s^4 ; h = syslin("c",num,den) fr = freson(h) fr = [] gainplot(h,0.01,1) with my new function "ffreson.sci" : ffreson(h)//La nouvelle fonction. ans = 0.111164226146
2 . Example : sampled system. The exercise in "help" page of "freson" with sampling period Ts=0.01s : the program gives "fr=[]" . With sampling period Ts=0.04s you have a good result (problem of accuracy) : run the example in help and change sampling period. More explanations in my email "New problems..."on Scilab-users".
-
By Samuel GOUGEON (@sgougeon)
(In reply to Lucien Povy from comment #12)
1 . Example : continuous system. s=%s ; num = s+1.4s^2+1.4s^3+0.4s^4 ; den = 0.5+0.8s+1.4s^2+1.4s^3+0.4*s^4 ; h = syslin("c",num,den) fr = freson(h) fr = []
Confirmed
2 . Example : sampled system. The exercise in "help" page of "freson" with sampling period Ts=0.01s : the program gives "fr=[]" .
Confirmed.
With sampling period Ts=0.04s you have a good result (problem of accuracy) : run the example in help and change sampling period. More explanations in my email "New problems..."on Scilab-users".
-
By Samuel GOUGEON (@sgougeon)
(In reply to Lucien Povy from comment #12)
1 . Example : continuous system. s=%s ; num = s+1.4s^2+1.4s^3+0.4s^4 ; den = 0.5+0.8s+1.4s^2+1.4s^3+0.4*s^4 ; h = syslin("c",num,den) fr = freson(h) fr = [] gainplot(h,0.01,1) with my new function "ffreson.sci" : ffreson(h)//La nouvelle fonction. ans = 0.111164226146
Just replacing hh = horner(h,-%s) with hh = clean(h*horner(h,-%s), 0, %eps)
gives the same and passes as well the existing tests:
--> format(18) --> fr = freson(h) fr = 0.111164226145773
2 . Example : sampled system. The exercise in "help" page of "freson" with sampling period Ts=0.01s : the program gives "fr=[]" .
I am no longer able to reproduce it. Since the used h is not part of the discrete example, i think i initially used the current h value without paying attention, and it was not the one defined with the continuous example. Now i get: --> h=syslin('c',-1+%s,(3+2*%s+%s^2)(50+0.1%s+%s^2)); --> hd=dscr(h,0.01); --> fr=freson(hd) fr = 1.188645012628571
-
By Samuel GOUGEON (@sgougeon)
(In reply to Samuel GOUGEON from comment #14)
Just replacing hh = horner(h,-%s) with hh = clean(h*horner(h,-%s), 0, %eps)
gives the same and passes as well the existing tests:
Sorry, bad copy/Paste. The fixing replacement made is
hh = clean(h*horner(h,-%s),0,%eps)replacing hh = h*horner(h,-%s)
-
By Samuel GOUGEON (@sgougeon)
By the way, using hh = clean(hhorner(h,-%s), 0, %eps) is also enough to fix the initial 15368 report, and we can go back to k = find(imag(r)>0 & abs(real(r))<%epsabs(r)); instead of k = find(imag(r)>0 & abs(real(r))<1.e-14*abs(r)); without breaking the test bug_15368.tst
-
By Samuel GOUGEON (@sgougeon)
(In reply to Samuel GOUGEON from comment #14)
2 . Example : sampled system. The exercise in "help" page of "freson" with sampling period Ts=0.01s : the program gives "fr=[]" .
I am no longer able to reproduce it. Since the used h is not part of the discrete example, i think i initially used the current h value without paying attention, and it was not the one defined with the continuous example.
Actually, i reproduce the issue on a computer, but not on another one, both running Scilab 6.0.2 on Win7 64 bits.