bugGNU Octave - Bugs: bug #61355, [octave forge] (control) Error on...

 
 

bug #61355: [octave forge] (control) Error on phase value of bode and margin plots

Submitter:  Luiz Antonio Maccari Junior <luiz>
Submitted:  Tue 19 Oct 2021 04:23:03 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  None Assigned to:  None
Originator Name:  Luiz Antonio Maccari Jr. Open/Closed:  * Open
Release:  * 6.3.0 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 19 Dec 2023 05:48:16 PM UTC, comment #15: 

Okay, I should have said that my changes were targeting the issue described in comment #10.

Torsten Lilge <ttl>
Group Member
Tue 19 Dec 2023 10:06:04 AM UTC, comment #14: 

Hi,

Making a small test with a pure oscillator the phase value still present problems. The phase goes to +180º but should be -180º.

Code:

FT1 = tf(1,[1 0 100])

bode(FT1)


Luiz Antonio Maccari Junior <luiz>
Mon 18 Dec 2023 09:44:32 PM UTC, comment #13: 

Sorry, I forgot to attach the file _frequency_response_.m, which is also updated for the fix.

(file #55460, file #55461)

Torsten Lilge <ttl>
Group Member
Mon 18 Dec 2023 09:31:08 PM UTC, comment #12: 

While testing the patch I came across the issue that the phase correction might not work when the user specifies an omega range that starts at high frequencies where the phase might be already far away from the initial phase. This effect also affects bode plots without poles or zeros at the origin: the phase for an omega range at high frequencies might be wring since the "history" of the phase starting at low frequencies is missing.

I have attached a version of bode.m which should fix this (including Juan's fix for zeros or poles at the origin). Some tests and feedback are very much appreciated.

(file #55459)

Torsten Lilge <ttl>
Group Member
Sun 19 Nov 2023 09:11:09 AM UTC, comment #11: 

Juan, thanks for the proposed solutions. The one for the initial phase looks very promising and I will have a closer look to this. Regarding the cases with purely imaginary pole/zero pairs, I am not sure if it is really necessary to handle all subsystems separately. But I do not have any other ideas right now.

Torsten Lilge <ttl>
Group Member
Sun 13 Aug 2023 11:57:14 AM UTC, comment #10: 

I've been looking at the implementation of the bode function.
As I understand it:
- Takes the system and calls _frequency_response_, which
- Uses _frequency_vector_ to compute a set of frequencies in a representative range given system's singularities
- Evaluates system for said frequencies and returns a vector of complex numbers
- Magnitude is simply calculated with the abs() function
- Phase is calculated using the arg() function (returns phase from -pi/2 to pi/2)
- Phase is "corrected" with the unwrap() function, which adds +- 2*pi until the difference between consecutive phase values is less than pi.
- As unwrap() can't correct the initial phase, additional +-2pi increments are added as a function of the number of poles and zeros at the origin
(This last step is skipped by margin(). The bode-drawing code is otherwise duplicated in bode.m and margin.m. This is the reason why bode() and margin() phases differ)

This approach has several problems:

Initial phase
--------------

Example: Two transfer functions with two poles at the origin. Low frequency phase should be: 2 poles at the origin * -90º = -180º for both.

G = (0.1*s+1)/(s^2*(0.01*s+1))
H = (0.1*s+1)/(s^2*(s+1))

1. Bode starts evaluating G and H at 0.1 rad/s. unwrap() can't do anyhting about the first value it is fed.

G(0.1*i) = -1.0000e+02 - 9.0000e-01i --> phase -179.48º OK
H(0.1*i) = -9.9109e+01 + 8.9109e+00i --> phase 174.86º Wrong, should be -185.14º.

2. Two poles at the origin result in a correction of -360º to the phase of both functions (as implemented).

G(0.1*i) --> phase -179.48º - 360º = -539.48º Wrong
H(0.1*i) --> phase 174.86º - 360º = -185.14º OK

Proposed solution:

Low frequency (asymptotic) phase is given by (zeros_at_origin - poles_at_origin) * 90º.
Computed phase should be corrected by adding +-360º increments so its first value is closest to the asymptotic one.
For example: phase += round((asymptotic_phase - first_phase)/360)*360

I have implemented and lightly tested this solution on a copy of the bode.m file found in https://github.com/gnu-octave/pkg-control/tree/main/inst

---------------------------------------------
## check for poles and zeroes at the origin for each of the numsys systems

% Initial phase correction so it is closest to asymptotic value
  initial_phase = cellfun(@(x) x(1), pha)
  for h = 1:numsys
sys1 = varargin{sys_idx(h)};
[num,den] = tfdata (sys1,'v');

n_poles_at_origin = sum (roots (den) == 0);
n_zeros_at_origin = sum (roots (num) == 0);
asymptotic_low_freq_phase = (n_zeros_at_origin - n_poles_at_origin)*90;
pha(h)={cell2mat(pha(h))+round((asymptotic_low_freq_phase - initial_phase(h))/360)*360};
  endfor
---------------------------------------------

Resonances
-----------

Pure imaginary poles or zeros pose an additional difficulty as phase changes abruptly from 0 to -180º for pole pairs and from 0 to +180º for zero pairs.
As unwrap() has a threshold of pi radians (180º) for incrementing +- 2pi, numerical errors can may it flip the phase the wrong way, as seen on the examples presented by Luiz. Another example: bode(1 / ( (s^2+1)*(s + 1) )), high frequency phase should go towards 3 poles*-90º = -270º but octave gives +90º.

Even Matlab's bode has trouble with this, in R2023a:
G = 1 / ( (s^2+1)*(0.1*s^2+1) ) is plotted correctly, with phase 0º -> -180º -> -360º
G = 1 / ( (s^2+1)*(0.01*s^2+1) ) is plotted incorrectly, with phase -360º -> -180º -> -360º

A further problem, as Luiz pointed out, is the phase at resonance, which comes out as NaN, and putting this through unwrap() results in all further phases being NaN.

To resolve this, my first instinct would be to compute poles and zeros of the transfer function, get a frequency response for each real or complex singularity (handling resonance appropriately), and then adding up all the responses together. That way, phase is guaranteed to be correct. Something tells me this may have some numerical issues that may result in inaccuracies though.

Please let me know your thoughts on these matters.
    Juan

Juan <juang>
Sat 12 Aug 2023 01:45:16 PM UTC, comment #9: 

Here Is what I think is the problem.

The original question is about having multiple roots on the imaginary axis. (2 oscillators at the same frequency.) What is happening is:
1) the code does all the math and comes up with a combined math statement.
2) Bode now has to decompose this complex statement back into it individual roots.
3) Using floating point arithmetic there is round off errors.
4) These round off errors tend to answers that are close to the . but not the same.
5) since the original roots were right on the imaginary axis then some of the round off error roots end up in the right hand plane.
6) this can be demonstrated by a=roots(poly ([ 1j 1j 1j 1j 1j ]))
7) all the roots should be 1j but some have a positive real part!
a =

   0.0008 + 0.9995i
  -0.0002 + 0.9990i
   0.0008 + 1.0006i
  -0.0010 + 0.9999i
  -0.0004 + 1.0009i

 One way that I can see to get this fixed is to use more digits etc. and then round out the errors.
 The second way would be to somehow keep the information used in the original part of the code and bypass the the complexity of have to fined the roots. But I don't have any ideas on how this would work.

Juan:
   Your  code does not seem to be this problem so I am not sure if you have exposed a real bug or not.
 


Doug Stewart <dastew>
Sat 12 Aug 2023 10:55:29 AM UTC, comment #8: 

Any news on this? I also experienced similar bugs.

For example, for this script, Bode for G and Lead_TF are correct, but not for the product, in which phase is off by 360º
-------------------------------------
  pkg load control

  s = tf('s');

  G = 10e6/(s^2);

  alpha = 10;
  tau = 562.3e-6;
  Lead_TF = (1/alpha)*((alpha*tau*s+1)/(tau*s+1));

  figure
  bode(G, Lead_TF, Lead_TF * G)
-------------------------------------

Juan <juang>
Tue 14 Jun 2022 03:53:57 PM UTC, comment #7: 

Is there any news about this bug report?

Is there some way of using the margin function with correct results in these cases?

Luiz Antonio Maccari Junior <luiz>
Mon 25 Oct 2021 09:08:55 PM UTC, comment #6: 

An additional bug was detected. If you insert the values of the frequencies in the bode  function for frequencies after the resonance the phase values return NaN. Example:

s = tf('s');

TF2 = 1/(s^5*(s+2)*(s+10)*(s^2+100))

w = [9:0.1:11];
[mag, fase] = bode(TF2,w)


However if we use the bode function without specifying the frequencies the bode plot can give a phase plot with finite values.

 Studying the bode function I see the problem related with the function unwrap which uses a summation and when we have an ideal resonance the argument of the frequency response is NaN. To bypass this error it is needed to put a damping values for the resonant. As for example:

s = tf('s');

TF2 = 1/(s^5*(s+2)*(s+10)*(s^2+0.001+100))

w = [9:0.1:11];
[mag, fase] = bode(TF2,w)


Luiz Antonio Maccari Junior <luiz>
Mon 25 Oct 2021 08:57:13 PM UTC, comment #5: 

An additional example of the error with the phase plots. The code below results in the same phase plot. Bode function does not make difference between the phase of doubled poles or doubled zeros.

s = tf('s');

TF1 = (s^2+100)/(s^5*(s+2)*(s+10))

TF2 = 1/(s^5*(s+2)*(s+10)*(s^2+100))

bode(TF1,TF2)

Luiz Antonio Maccari Junior <luiz>
Wed 20 Oct 2021 12:35:50 AM UTC, comment #4: 

I am posting two additional examples of error in margin and bode functions.

Example 4 in file "bug_report_bode_ex4.m".

The bode plot of matlab and octave are the same both are correct. However, the margin function is incorrect in octave. The phase plot in octave margin phase plot starts in +180 and it should be -180.
See files: matlab_output_margin_ex4.jpg and octave_output_margin_ex4.jpg


Example 5 in file "bug_report_bode_ex5.m".

Here the bode function phase plot is wrong. It should start with 0 degree but it is starting with +180.

Sometimes the phase plot is wrong in bode function and sometimes is wrong in margin function. In margin function the error also leads to an error on gain margin. It seems that the error is in the wrapping of the phase plot. In example 5 bode function plotted a wrong phase range and in example 4 the error is in the phase range of margin function.

As I mentioned before, Matlab does not put wrap as default giving an option to choose it. Also in matlab it is important to note that the wrapping does not interfere on the results of gain or phase margin.

Thank for your time.  



(file #52130, file #52132,

Luiz Antonio Maccari Junior <luiz>
Tue 19 Oct 2021 11:21:38 PM UTC, comment #3: 

I can confirm that you have found a bug.
It seems to appear when there are 2 oscillators.
I am doing more tests.

Doug Stewart <dastew>
Tue 19 Oct 2021 08:58:45 PM UTC, comment #2: 

Now I was reviewing the functions in matlab and there they have an option named wrap phase. Using this all the phase plots will be plotted between the range -180 to +180. However, this option does not affect the value of gain margin. And for octave it affects causing errors. As for example in Example 3.

Another thing is that in octave there is not a standard form to plot phase. Sometimes the phase plots using bode and margin are different.

The ideal would be a non wrapped bode and margin functions or an option to choose the phase plot form.

Luiz Antonio Maccari Junior <luiz>
Tue 19 Oct 2021 04:33:15 PM UTC, comment #1: 

I have attached the results of the example 3 for matlab margin function (ex3_matlab.jpg) and octave margin function (ex3_octave.jpg). The result in octave is wrong.


Luiz Antonio Maccari Junior <luiz>
Tue 19 Oct 2021 04:23:03 PM UTC, original submission:  

Hello,

I am using bode function and I got wrong results when the system has two poles in the same position. The phase should be plotted as -180 (as two poles) however in some cases it is plotted as +180 (as two zeros). ? This behavior leads to a wrong plot and phase margins when I call margin function.

Above I give three examples with distinct results.

Example 1:

First one in the file "bug_report_bode.m" is a system with 4 poles, two at the frequency of 10^2 rad/s and two at 10^4 rad/s.

The expected behavior is -180° and after -180°. However the bode plot presents +180 and after -180. Bode and margin functions give the same phase plot. I have tested the same script in matlab and there the results are correct, as you can see on figures attached.

right result = matlab_output_bode.jpg
wrong result = octave_output_bode.jpg

Example 2:

On this example the bode plot gives the correct plot but the margin function does not. The example "bug_report_bode_teste3.m" has two poles at origin and additional one at 10 rad/s.The phase plot of function bode is correct but the margin one is wrong.


Example 3:
These wrong phase plot can lead to a problem to identify the correct gain margin as for example in the file "bug_report_bode_ex3.m" where because the wrong plot there is no crossover at -180 phase angle.

Best Regards

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 #55460:  bode.m added by ttl (6KiB - text/x-objcsrc)
file #55461:  __frequency_response__.m added by ttl (5KiB - text/x-objcsrc)
file #55459:  bode.m added by ttl (6KiB - text/x-objcsrc)
file #52130:  bug_report_bode_ex4.m added by luiz (171B - text/x-objcsrc)
file #52132:  bug_report_bode_ex5.m added by luiz (192B - text/x-objcsrc)
file #52127:  ex3_matlab.jpg added by luiz (25KiB - image/jpeg)
file #52128:  ex3_octave.jpg added by luiz (63KiB - image/jpeg)
file #52124:  octave_output_bode.jpg added by luiz (56KiB - image/jpeg)
file #52125:  bug_report_files.zip added by luiz (844B - application/zip)
file #52126:  matlab_output_bode.jpg added by luiz (19KiB - image/jpeg)

 

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 juang (Posted a comment)
  • -email is unavailable- added by dastew (Posted a comment)
  • -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-12-19 luiz Attached File- Added Screenshot@from@2023-12-19@06-56-52.jpg, #55462
    2023-12-18 ttl Attached File- Added bode.m, #55460
        Attached File- Added _frequency_response_.m, #55461
    2023-12-18 ttl Attached File- Added bode.m, #55459
    2021-10-20 mmuetzel Carbon-CopyRemoved 102357 -
    2021-10-20 mmuetzel SummaryError on phase value of bode and margin plots [octave forge] (control) Error on phase value of bode and margin plots
    2021-10-20 luiz Attached File- Added bug_report_bode_ex4.m, #52130
        Attached File- Added matlab_output_margin_ex4.jpg, #52131
        Attached File- Added bug_report_bode_ex5.m, #52132
        Attached File- Added octave_output_margin_ex4.jpg, #52133
    2021-10-19 luiz Attached File- Added ex3_matlab.jpg, #52127
        Attached File- Added ex3_octave.jpg, #52128
    2021-10-19 luiz Attached File- Added octave_output_bode.jpg, #52124
        Attached File- Added bug_report_files.zip, #52125
        Attached File- Added matlab_output_bode.jpg, #52126

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code