bugGNU Octave - Bugs: bug #63393, [octave forge] (control) Bode plot...

 
 

bug #63393: [octave forge] (control) Bode plot problem. Different results between MATLAB and Octave.

Submitter:  Luiz Antonio Maccari Junior <luiz>
Submitted:  Mon 21 Nov 2022 01:57:52 PM UTC
   
 
Category:  Octave Forge Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Confirmed Assigned to:  None
Originator Name:  Luiz Antonio Maccari Jr. Open/Closed:  * Open
Release:  * 7.2.0 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 06 Feb 2023 05:08:14 PM UTC, comment #23: 

Yes these two issues are happening here and another one is related with the function unwrap which uses a summation and when we have an ideal resonance the argument of the frequency response is NaN. This behaviour causes differences in the stability margins between MATLAB and OCTAVE. This last one is described in bug #61355.




Luiz Antonio Maccari Junior <luiz>
Fri 03 Feb 2023 08:21:20 PM UTC, comment #22: 

It seems that the bode plot of proportional but unstable systems starts with a phase not equal to zero in Matlab, like, e.g.


G = tf (1,[1 -2 1]);
bode (G)


Therefore, the phase for your example in Octave looks correct, except for the phase at 400 < w < 1000 when not using a sufficiently high number of frequency samples.

Actually, when applying Dimitri's patch in comment #17 and the frequency range from comment #13, the bode plot looks good on my system.

To summarize, we have the following two issues:

  • Evaluation of polynomials not numerically robust
  • Numerical issues (especially in the phase) depend on the number frequency samples of the bode plot

Or am I missing something?

Torsten Lilge <ttl>
Group Member
Fri 03 Feb 2023 06:34:35 PM UTC, comment #21: 

Good point.

This system came from a discretization of 4 resonant controllers in series with a system which also has a resonance.

In a first sight the initial phase should be zero. In this case it seems that MATLAB shows a wrong phase starting point.


Luiz Antonio Maccari Junior <luiz>
Fri 03 Feb 2023 02:26:26 PM UTC, comment #20: 

Which phase is correct for f -> 0? Isn't it 0°? Or are there really 8 integrators in the system?

Torsten Lilge <ttl>
Group Member
Fri 03 Feb 2023 12:30:40 PM UTC, comment #19: 

The comment #13 seem to be related with the bug #61355.

https://savannah.gnu.org/bugs/?61355

Using the suggested frequency range the phase plot assumes the same pattern of the MATLAB one (b6.png). But it still begins on zero in octave and -720 in matlab.

The noise is not affected but the phase plot is improved, although not yet equal as the MATLAB.





Luiz Antonio Maccari Junior <luiz>
Fri 03 Feb 2023 12:11:40 PM UTC, comment #18: 

Hi,

I've made the modifications cited in comment #17.

I've got the results in b4.png and b5.png . However the noise is still there and shape also is not correct. Look to error on the dc gain.




Luiz Antonio Maccari Junior <luiz>
Fri 03 Feb 2023 04:23:54 AM UTC, comment #17: 

I made the following change in @tf/__freqresp__.m :



36d35
<
38c37,44
<     H = polyval (num, s) ./ polyval (den, s);
---
>     Hnum0 = polyval (num, s);
>     [num1,~,mu] = polyfit(s, Hnum0, length(num));
>     Hnum1 = polyval (num1, s, [], mu);
>     Hden0 = polyval (den, s);
>     [den1,~,mu] = polyfit(s, Hden0, length(den));
>     Hden1 = polyval (den1, s, [], mu);
>
>     H = Hnum1 ./ Hden1 ;


(and also compiled octave and control package with AVX instructions instead of default sse -- perhaps this is important too) and the result is attached. Still not perfect, but looks better.

Dmitri.
--



Dmitri A. Sergatskov <dasergatskov>
Fri 03 Feb 2023 03:18:34 AM UTC, comment #16: 

I would like to withdraw  the previous comment. It is not that easy. Octave should improve accuracy of polyval() by replacing
classical Horner algorithm with more modern ones (like compensated Horner).

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Fri 03 Feb 2023 02:40:27 AM UTC, comment #15: 

Replacing polyval (P, X) calls with polyval (P, X, [], MU)
may improve numerical stability.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Fri 03 Feb 2023 02:13:23 AM UTC, comment #14: 

It seems that the errors result from the polyval() function applied to the denominator of the system.

Torsten Lilge <ttl>
Group Member
Thu 02 Feb 2023 10:37:19 PM UTC, comment #13: 

I can confirm the issue from comment #12. Besides the numerical errors leading to the "noise", the phase plot is different between 400 and 1000 rad/s compared to the result from Matlab. However, this effect disappears if you select another range of frequencies, like, e.g.


freq = logspace(0,5,10000);


The poles (and most of the zeros) of the systems are all very close to 1, which might be the cause of the numerical problems in Octave. What results do you get when you discretize the original continuous-time system with a larger sampling time?

Torsten Lilge <ttl>
Group Member
Thu 02 Feb 2023 02:13:13 PM UTC, comment #12: 

Now I was working again with the bode function an I have realized that the different behaviour is not only because Matlab reducing the order of the model or something like that.

Please see the following example:


num =[                   0;
                         0;
     3.297563545330486e-05;
    -6.662437436850499e-05;
    -2.565172417884272e-04;
     1.098159870753534e-03;
    -1.616164103992111e-03;
     1.098654441170365e-03;
    -2.570132570435638e-04;
    -6.641072266842919e-05;
     3.293973876807465e-05;
     1.834542499710346e-11;];

den =[      1.000000000000000e+00;
    -1.079593440596339e+01;
     5.316681632945271e+01;
    -1.576732947973972e+02;
     3.128986175251184e+02;
    -4.363042796274709e+02;
     4.362158752637891e+02;
    -3.127087912934111e+02;
     1.575144575020965e+02;
    -5.309223626372895e+01;
     1.077659943863232e+01;
    -9.978296711173387e-01 ;];

TF = tf(num',den',5e-05)

freq =[1:1:2*pi*0.5/5e-05];
bode(TF,freq)


Results are in image b1.png and a zoom can be seen in b2.png. Matlab result is in the left and octave result in the right.

As one can notice the "noise" appears only in Octave image. The order of both results are the same however the octave result shows a lot of numerical errors.

This type of behaviour for the case high-order systems makes impossible to correct evaluate the results from octave. Takes from example other system which I working now the figure b3.png. The results from octave seems to be completely wrong.

Best Regards.






Luiz Antonio Maccari Junior <luiz>
Mon 21 Nov 2022 08:47:32 PM UTC, comment #11: 

I took a look at what minreal was doing, and It was reducing the order of your filter. I think you should reduce the order of the filter  at the start rather than using minreal.

If you want to stay with that order of filter for your work then just use it and don't worry about what bode shows you. But also realize that the is a limit on using floating point resolution.

Doug Stewart <dastew>
Mon 21 Nov 2022 04:38:08 PM UTC, comment #10: 

I am sorry -- I misunderstood the report.
I can reproduce the problem.

The figure I attached earlier is
from the sample code  in comment #2 (with "minreal").

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Mon 21 Nov 2022 04:34:40 PM UTC, comment #9: 

I have tested here using the following command:

LD_PRELOAD=/lib/x86_64-linux-gnu/blas/libblas.so /home/luiz/my_octave/bin/octave --gui

The results are the same. I've got the result with noise.




Luiz Antonio Maccari Junior <luiz>
Mon 21 Nov 2022 04:08:54 PM UTC, comment #8: 

I am using now Linux Mint 20.3. I have built it following the instruction on octave wiki.

https://wiki.octave.org/Building

It is important to notice that the problem that I've reported was found also on windows 10 with the version 7.3.




Luiz Antonio Maccari Junior <luiz>
Mon 21 Nov 2022 04:05:18 PM UTC, comment #7: 

On Windows -- there is a blas switcher utility program is included with the installation. On Linux you can start octave as



LD_PRELOAD=/usr/lib64/libblas.so octave


(adjust the path to libblas.so as appropriate)

I think nn MacOS an equivalent would be DYLD_INSERT_LIBRARIES,
but I am not 100% sure.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Mon 21 Nov 2022 03:55:24 PM UTC, comment #6: 

for switching BLAS, which operating system are you using and how did you install octave?

Nicholas Jankowski <nrjank>
Group Member
Mon 21 Nov 2022 03:53:18 PM UTC, comment #5: 

regrading minreal not having help text or being viewable in edit, I assume you're referring to the Octave implementations?

Because of the way Octave has implemented these as class functions, you can't just check the functions with 'help minreal', even after loading the package.

'help @lti/minreal'  will diplay the (somewhat minimal) package help.  'edit @lti/minreal' will open the function for editing.  Note that it's handing off most of the functionality to @ss/___minreal___.m, which in turn will call other helper functions (located under @ss/) depending on the form of the input.

Also, note that if you wanted to debug/step through the functions, setting breakpoints is problematic as you need to be in the same folder as the function with the breakpoint, and that doesn't work well for class functions in a @/ folder. What you can do instead is have all of the relevant functions open in your editor, then add a 'keyboard' command in the first one called. from there you can step through the functions and into subfunctions using the toolbar or F10/F11 keys as needed.

Nicholas Jankowski <nrjank>
Group Member
Mon 21 Nov 2022 03:49:38 PM UTC, comment #4: 

Sorry I don't know how to do it. But if you give me the instructions to change BLAS Library, I could try.

Luiz Antonio Maccari Junior <luiz>
Mon 21 Nov 2022 03:43:47 PM UTC, comment #3: 

It looks OK (attached bode_das1.png) on linux octave 7.3.0 and control 3.4.0

(Control package seems to be sensitive to blas libraries. Could you try to change you blas from, presumably, openblas to reference blas?)

Dmitri.
--



Dmitri A. Sergatskov <dasergatskov>
Mon 21 Nov 2022 02:57:06 PM UTC, comment #2: 

The new code is:
-------------------------
load T.mat

Ts = 1/20000;

TF = tf(num,den,Ts)
TF=minreal(TF,0.1)


freq=[0.01:1:2*pi*(0.5/Ts)];

bode(TF,freq)

------------------------------

With the insertion of minreal the results are apparently the same. It seems that is not a bug on octave just a different type of behavior of the bode function.

Reasoning about the problem, I think that MATLAB is doing something not requested by the user. Maybe an warning or an explanation on the help should be interesting.

Sadly the function minreal does not have a help text and I cant open it using edit command.

Thank by attention.

Luiz Antonio Maccari Junior <luiz>
Mon 21 Nov 2022 02:42:33 PM UTC, comment #1: 

Hi,

I have made some tests using the function minreal with a nonstandard tolerance (tol=0,1). It seems that the bode function from MATLAB does this automatically and Octave bode doesn't.

On my point of view this bug report can be closed.

Thanks.

Luiz Antonio Maccari Junior <luiz>
Mon 21 Nov 2022 01:57:52 PM UTC, original submission:  

Dear developers,

I have found a different behavior from MATLAB and Octave bode function. I am taking the frequency response of a high-order discrete transfer function and getting different results from MATLAB (Bode_bug_matlab.jpg) and from octave (Bode_bug_octave.jpg).

As one can see, the Octave returns a frequency response with high-order (with noise) and MATLAB returns a filtered one (without noise).

On my interpretation, MATLAB is cancelling some pole and zeros inside a unknown precision range. However, I am not sure if the octave result is correct or it is just plotting something due to numeric errors.

All the files to reproduce the error are attached.

Thanks.

        

 





All files are attached.  

Luiz Antonio Maccari Junior <luiz>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #54310:  b6.png added by luiz (67KiB - image/png)
file #54308:  b4.png added by luiz (35KiB - image/png)
file #54309:  b5.png added by luiz (37KiB - image/png)
file #54298:  b1.PNG added by luiz (66KiB - image/png)
file #54299:  b3.PNG added by luiz (73KiB - image/png)
file #54300:  b2.PNG added by luiz (82KiB - image/png)
file #53986:  bode_das1.png added by dasergatskov (40KiB - image/png)
file #53981:  Bode_bug_matlab.jpg added by luiz (20KiB - image/jpeg)
file #53982:  TF.mat added by luiz (648B - application/octet-stream)
file #53983:  Bode_bug_octave.jpg added by luiz (59KiB - image/jpeg)
file #53984:  bode_bug.m added by luiz (102B - text/x-objcsrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by ttl (Posted a comment)
  • -email is unavailable- added by dastew (Posted a comment)
  • -email is unavailable- added by dasergatskov (Updated the item)
  • -email is unavailable- added by nrjank (Updated the item)
  • -email is unavailable- added by luiz (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 15 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2023-02-03 luiz Attached File- Added b6.png, #54310
    2023-02-03 luiz Attached File- Added b4.png, #54308
        Attached File- Added b5.png, #54309
    2023-02-03 dasergatskov Attached File- Added das_control_mods_2023-02-02.png, #54304
    2023-02-02 luiz Attached File- Added b1.PNG, #54298
        Attached File- Added b3.PNG, #54299
        Attached File- Added b2.PNG, #54300
    2022-11-21 nrjank StatusNeed Info Confirmed
    2022-11-21 nrjank StatusNone Need Info
    2022-11-21 dasergatskov Attached File- Added bode_das1.png, #53986
    2022-11-21 nrjank SummaryBode plot problem. Different results between MATLAB and Octave. [octave forge] (control) Bode plot problem. Different results between MATLAB and Octave.
    2022-11-21 luiz Attached File- Added Bode_bug_matlab.jpg, #53981
        Attached File- Added TF.mat, #53982
        Attached File- Added Bode_bug_octave.jpg, #53983
        Attached File- Added bode_bug.m, #53984

    Back to the top

    Powered by Savane 3.12.
    Corresponding source code