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

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)

Answer 1

[Paul]

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.

[Alex Seaver]

Sorry, I'm not well-versed in the techniques you used to apply positive and negative feedback to the transfer function via the characteristic polynomial (especially, the addition or subtraction of the numerator/denominator). If you could recommend any resources for learning those techniques, that would be great. Nonetheless, it seems the solution is unfortunately to simply artificially invert the sign of transfer function. Thank you for the thorough response!

[Paul]

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