Bode plot with right half plane zero

Hi, while using the bode function in matlab script for a transfer function which has right half plane zero, I am finding that the phase is starting from 360 instead of starting from 0. How can this be corrected. The code is inserted.

clc

clear all

s=tf('s')

fs=1e6;

ws=2*pi*fs;

L=1.8e-6;

C=220e-6;

R=5;

Vin=9;

Kg=1/5;

Vcarr=1;

Vbase=5;

Ibase=5;

Zb=Vbase/Ibase

%% Plant TFs

Vo=5;

D=0.357;

wrz=((1-D)^2)*R/(D*L); %% position of the RHP zero

wc=(1-D)/sqrt(L*C);

Q=(1-D)*R*sqrt(C/L);

Gdo=Vo/(D*(1-D));

Gvd=Gdo*(1-(s/wrz))/(((s/wc)^2)+(s/(Q*wc))+1); % control to output

%% bode plots

wct=1:10:10e7;

bode(Gvd,wct);

Answer 1

[Paul]

Hi @Anindya,

The bodeplot function offers options on controlling the appearance of the plot, including to force the phase between +-180 deg.

A bigger hammer would be to use ctrlpref to set the default phase wrapping preference for bode and bodeplot (actually, I'm not 100% sure it applies to bodeplot, some experimentation may be in order).

[Mathieu NOE]

hello @Paul

+1

see my suggestion above

you can force the phase between +-180 deg with the regular bode function

Answer 2

[Andrew Ouellette]

Hi @Anindya,

You are looking for the phase matching options in bodeplot, which you can interpret as "If PhaseMatchingEnabled, then keep phase close to PhaseMatchingValue at the PhaseMatchingFrequency". Here is an example using your code:

s=tf('s');
fs=1e6;
ws=2*pi*fs;
L=1.8e-6;
C=220e-6;
R=5;
Vin=9;
Kg=1/5; 
Vcarr=1; 
Vbase=5;
Ibase=5;
Zb=Vbase/Ibase;
%% Plant TFs
Vo=5;
D=0.357;
wrz=((1-D)^2)*R/(D*L); %% position of the RHP zero
wc=(1-D)/sqrt(L*C);
Q=(1-D)*R*sqrt(C/L);
Gdo=Vo/(D*(1-D));
Gvd=Gdo*(1-(s/wrz))/(((s/wc)^2)+(s/(Q*wc))+1);  % control to output
%% bode plots
wct={1 10e7};
h = bodeplot(Gvd,wct);
h.PhaseMatchingEnabled = true;

●

h.PhaseMatchingValue = 0;
h.PhaseMatchingFrequency = 0;

Answer 3

[Claire]

When you use the bode function to plot the Bode diagram of a transfer function that has a right-half plane (RHP) zero, you may notice that the phase starts from 360 degrees instead of 0 degrees. This is because MATLAB does not normalize the phase when calculating it; therefore, when the phase exceeds 180 degrees, it will be displayed as a value greater than 180 degrees.

[Anindya]

yes. Is there any way to avoid that

Answer 4

[Mathieu NOE]

hello

you can do this way :

BTW , you pulstaion vector was huge and there is a better approach for Bode plots (as the x axis is log scaled ,use therefore a log spaced frequency / pulstaion vector)

here 500 values are enough to make a correct plot

clc
clear all
s=tf('s')

s =
 
  s
 
Continuous-time transfer function.

fs=1e6;
ws=2*pi*fs;
L=1.8e-6;
C=220e-6;
R=5;
Vin=9;
Kg=1/5; 
Vcarr=1; 
Vbase=5;
Ibase=5;
Zb=Vbase/Ibase

Zb = 1

%% Plant TFs
Vo=5;
D=0.357;
wrz=((1-D)^2)*R/(D*L); %% position of the RHP zero
wc=(1-D)/sqrt(L*C);
Q=(1-D)*R*sqrt(C/L);
Gdo=Vo/(D*(1-D));
Gvd=Gdo*(1-(s/wrz))/(((s/wc)^2)+(s/(Q*wc))+1);  % control to output
%% bode plot
freq=logspace(1,7,500);
[g,p] = bode(Gvd,2*pi*freq);
g = g(:);
p = p(:) - 360 ;

figure
subplot(2,1,1),semilogx(freq,20*log10(g))
ylabel('Modulus (dB)')
title('Bode Plot')
grid on
subplot(2,1,2),semilogx(freq,p)
xlabel('Frequency')
ylabel('Phase(°)')
grid on

[Anindya]

Agree, a log scale pulsation vector would be good. Thank you for pointing out that. Is there a way to do the bode plot with fewer lines of code? Thank you.

[Mathieu NOE]

you can create your own Bode function - and store it in your path (no need to keep it in the main script as below)

clc
clear all
s=tf('s');
fs=1e6;
ws=2*pi*fs;
L=1.8e-6;
C=220e-6;
R=5;
Vin=9;
Kg=1/5; 
Vcarr=1; 
Vbase=5;
Ibase=5;
Zb=Vbase/Ibase;
Zb = 1

Zb = 1

%% Plant TFs
Vo=5;
D=0.357;
wrz=((1-D)^2)*R/(D*L); %% position of the RHP zero
wc=(1-D)/sqrt(L*C);
Q=(1-D)*R*sqrt(C/L);
Gdo=Vo/(D*(1-D));
Gvd=Gdo*(1-(s/wrz))/(((s/wc)^2)+(s/(Q*wc))+1);  % control to output
%% bode plot
freq=logspace(1,7,500);
myBode(Gvd,freq);

%% function myBode %%%%%
function myBode(Gvd,freq)

    [g,p] = bode(Gvd,2*pi*freq);
    g = g(:);
    p = p(:);
    if p(1)> 1
        p = p - 360 ;
    end
    figure
    subplot(2,1,1),semilogx(freq,20*log10(g))
    ylabel('Modulus (dB)')
    title('Bode Plot')
    grid on
    subplot(2,1,2),semilogx(freq,p)
    xlabel('Frequency')
    ylabel('Phase(°)')
    grid on

end

[Paul]

Hi Mathieu,

I assume this function is only intended for SISO systems?

The check "if p(1) > 1" isn't robust. The p(1) could be anything depending on freq(1), and subtracting the 360 might not be warranted.

[Mathieu NOE]

sure

it was a quick and dirty solution for this specific case usage

a more robust solution can be to test the number of data for which p>0

[Mathieu NOE]

we could also use this (better) solution :

% force phase to [-180, 180] range : phase = mod(phase + 180, 360) - 180;

combined with unwrap function

clc
clear all
s=tf('s');
fs=1e6;
ws=2*pi*fs;
L=1.8e-6;
C=220e-6;
R=5;
Vin=9;
Kg=1/5; 
Vcarr=1; 
Vbase=5;
Ibase=5;
Zb=Vbase/Ibase;
Zb = 1

Zb = 1

%% Plant TFs
Vo=5;
D=0.357;
wrz=((1-D)^2)*R/(D*L); %% position of the RHP zero
wc=(1-D)/sqrt(L*C);
Q=(1-D)*R*sqrt(C/L);
Gdo=Vo/(D*(1-D));
Gvd=Gdo*(1-(s/wrz))/(((s/wc)^2)+(s/(Q*wc))+1);  % control to output
%% bode plot
freq=logspace(1,7,500);
myBode(Gvd,freq);

%% function myBode %%%%%
function myBode(Gvd,freq)
    [g,p] = bode(Gvd,2*pi*freq);
    g = g(:);
    p = p(:);
    p = mod(p + 180, 360) - 180;
    p = 180/pi*unwrap(p*pi/180);
    
    figure
    subplot(2,1,1),semilogx(freq,20*log10(g))
    ylabel('Modulus (dB)')
    title('Bode Plot')
    grid on
    subplot(2,1,2),semilogx(freq,p)
    xlabel('Frequency')
    ylabel('Phase(°)')
    grid on
end

[Paul]

In some weird cases, unwrap can have problems, IIRC.

If I want to manually force the output of bode to -180 <= p <= 180, I just do

p = angle(exp(1j*p*pi/180))*180/pi;

Though with the more "modern" functions perhaps

p = atan2d(sind(p),cosd(p));

[Mathieu NOE]

that's good !