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

6 views (last 30 days)
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)

Accepted Answer

Paul
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 Comments
Paul
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).

Sign in to comment.

More Answers (0)

Categories

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

Community Treasure Hunt

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

Start Hunting!

Translated by