# Measure Phase/Gain Margin from a different phase reference using margin()?

6 views (last 30 days)
Alex Seaver on 24 Oct 2022
Commented: Paul on 31 Oct 2022
I'm attempting to plot a bode diagram with phase and gain margins labeled using the margin() function. My transfer function is inverting - thus, its phase starts near 180°. When margin() is called, it plots the system properly. However, it measures the gain and phase margins with respect to 180° instead of 0° or 360° (the inversion points for an initially inverting system). How can I measure the margin distances from a different reference phase?
For instance, in the code below the phase margin is measured as -142° (that is, 142° below 180°). However, the true phase margin is around 38° (that is, the phase distance remaining before the system is inverted at 0°). Is this problem occuring becuase the phase doesn't start at exactly 180°? (That's another problem I'm not sure how to solve).
Here's my code (I'm using symbolic expressions to setup the transfer function):
clc; clear; close all;
% setup symbolic evaluation of transfer function
syms s
T0 = sym(-97.74e3); % T0, the DC loop transmission
f0 = sym(5); % system's first pole in Hz
f1 = sym(5e6); % system's second pole in Hz
fs = sym(286.2e3); % additional stray pole
Ts = T0*(1/(1+s/(2*pi*f0)))*(1/(1+s/(2*pi*f1)))*(1/(1+s/(2*pi*fs))); % loop transmission
pretty(Ts)
97740 - --------------------------------------------------- / s \ / s \ / s \ | ----- + 1 | | --------- + 1 | | ----------- + 1 | \ 10 pi / \ 572400 pi / \ 10000000 pi /
% preprocess the symbolic expression for use in a tf object
Ts = expand(Ts); % expand into polynomial
[symNum,symDen] = numden(Ts); % separate numerator and denominator
num = sym2poly(symNum); % convert to array of coefficients
den = sym2poly(symDen);
% set the options for the Bode plot
options = bodeoptions;
options.FreqUnits = 'Hz'; % set Bode plot units of frequency to Hz
options.Grid = 'On'; % plot with a grid for easier viewing
options.XLim = [1, 1e9]; % set frequency plotting range
% plot with phase and gain margins labeled
sys = tf(num, den);
margin(sys, options)

Paul on 24 Oct 2022
Edited: Paul on 25 Oct 2022
Hi Alex,
I'd look at the problem as follows
syms s
T0 = sym(-97.74e3); % T0, the DC loop transmission
f0 = sym(5); % system's first pole in Hz
f1 = sym(5e6); % system's second pole in Hz
fs = sym(286.2e3); % additional stray pole
Ts = T0*(1/(1+s/(2*pi*f0)))*(1/(1+s/(2*pi*f1)))*(1/(1+s/(2*pi*fs))) % loop transmission
Ts =
% preprocess the symbolic expression for use in a tf object
Ts = expand(Ts); % expand into polynomial
[symNum,symDen] = numden(Ts); % separate numerator and denominator
Compute the characteristic polynomial assuming negative feedback:
Delta = symDen + symNum
Delta =
The constant term is negative, hence the closed-loop system, assuming negative feedback, is unstable.
Compute the characteristic polynomial assuming positive feedback:
Delta = symDen - symNum
Delta =
The closed-loop poles, assuming positive feedback, are:
vpasolve(Delta)
ans =
The poles all have negative real parts, so the closed-loop system is stable. With that information, we can now proceed to compute the gain and phase margin of the system
num = sym2poly(symNum); % convert to array of coefficients
den = sym2poly(symDen);
% set the options for the Bode plot
options = bodeoptions;
options.FreqUnits = 'Hz'; % set Bode plot units of frequency to Hz
options.Grid = 'On'; % plot with a grid for easier viewing
options.XLim = [1, 1e9]; % set frequency plotting range
% plot with phase and gain margins labeled
sys = tf(num, den);
margin assumes that they system uses negative feedback. But we want the system to use positive feedback. Therefore, we negate the open-loop transfer function
margin(-sys, options)
GM is 20.7 dB and P is 37.8 deg (as you expected?).
Compute the actual gain margin to verify it.
gm = margin(-sys);
Compute the roots of the closed-loop characteristic polynomial with the GM applied in the loop with positive feedback
format long e
roots(den - gm*[0 0 0 num])
ans =
-3.321420559583051e+07 + 0.000000000000000e+00i 4.545625532045960e-03 + 7.516292925894370e+06i 4.545625532045960e-03 - 7.516292925894370e+06i
As expected, we have two poles on the imaginary axis at a frequency of
imag(ans(2))/2/pi % Hz
ans =
1.196255172882734e+06
as shown on the plot.
##### 2 CommentsShowHide 1 older comment
Paul on 31 Oct 2022
Assume we have a plant G(s) and controller K(s) in series.
For simplicity, define
L(s) = K(s)*G(s) = N(s)/D(s).
Assuming unity, negative feedback for a closed-loop syste, the closed-loop transfer function is:
H(s) = L(s)/(1 + L(s)) = N(s)/(D(s) + N(s)). The poles of the system are the roots of the denominator: D(s) + N(s) = 0
With positive feedback, the closed-loop TF is
H(s) = L(s) / (1 - L(s)) and the denominator of that is D(s) - N(s).

### Categories

Find more on Get Started with Control System Toolbox in Help Center and File Exchange

R2022b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by