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

[Juan]

Follow-up Comment #10, bug #61355 (project octave):

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


    _______________________________________________________

Reply to this item at:

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

_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/