MATLAB Answers

Zero-crossings of Nyquist Diagram

71 views (last 30 days)
Leo
Leo on 2 Feb 2021
Edited: Paul on 12 Feb 2021
Hi,
I have been trying to extract all the zero-crossings of the Nyquist plot of my transfer function zp. So I want to obtain the real part value where the Nyquist plot crosses the Real Axis.
Here is an example of how the Nyquist diagram could look like (I keep changing specific parameters so the Nyquist diagram also keeps changing):
I already tried extracting the zero crossings with the folloing code,
[re,im] = nyquist(zp);
Re = squeeze(re);
Im = squeeze(im);
With the given values I searched where the imaginary part was close to zero. But unfortunately it isn't precise enough and what I need is very detailed data. Increasing the amount of values also didn't get me much further. This is the code I used to increase my amount of values:
upperbound = 2*pi*(1600);
lowerbound = 2*pi*(-1600);
increase_number_of_Values = linspace(lowerbound,upperbound,5000);
[re,im] = nyquist(zp, increase_number_of_Values);
Re = squeeze(re);
Im = squeeze(im);
This code doesn't give me every zero-crossing of the Nyquist plot. I already changed the frequency interval various times, yet it doesn't show me every zero-crossing. Then I tried to extract the zero crossings directly from the transfer function itself. With the following code:
fp = [-1600:0.1:1600];
zp_values = [];
for k = 1 : length(fp)
w= 2*pi*fp(k);
zp_values(k) = (my transfer function); %it calculates the value of zp for every frequency in my array fp
end
However this also did not work. I chose the frequency interval -1600Hz <= fp <= 1600Hz because the values which I want to observe in the Nyquist plot have a frequency between -10000 rad/s to 10000 rad/s.
Has anybody an idea how I can extract the zero-crossings of my Nyquist plot, as accurately as possible?
Thanks very much for your help!

Accepted Answer

Kiran Felix Robert
Kiran Felix Robert on 5 Feb 2021
Hi Leopold,
I understand that you need to find the zero-crossings of the imaginary part of the Nyquist diagram. I assume that you have already extracted the real-part and the imaginary part separately.
In that case you can refer the following answer thread for a solution,
The following is a specific example
H = tf([2 5 1],[1 2 3]);
[re,im] = nyquist(H);
re = squeeze(re);
im = squeeze(im);
plot(re,im)
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);
zind = zci(im);
figure(1)
plot(re,im,'r');
hold on
plot(re(zind),im(zind),'bp')
hold off
grid
You can also try interpolating your data for a more accurate values.
  2 Comments
Paul
Paul on 12 Feb 2021
Wouldn't it be true by symmetry that if the Nyquist plot crosses the real axis at frequency w it must also cross the real axis (at the same point) at -w, as shown in your plot?

Sign in to comment.

More Answers (1)

Paul
Paul on 12 Feb 2021
Edited: Paul on 12 Feb 2021
By definition, the gain margin is defined at the point(s) where the Nyquist plot crosses the real axis. Use the allmargin() function to find the gain margin(s). The actual crossing point is at 1/GM.
>> h=tf([2 5 1],[1 2 3])*tf(1,[1 2 2 1]) % example transfer function
h =
2 s^2 + 5 s + 1
--------------------------------------
s^5 + 4 s^4 + 9 s^3 + 11 s^2 + 8 s + 3
Continuous-time transfer function.
>> nyquist(h)
>> s=allmargin(h)
s =
struct with fields:
GainMargin: 2.250221178446187e+00
GMFrequency: 1.870895057257997e+00
PhaseMargin: [1.675188645227028e+02 5.814149251422266e+01]
PMFrequency: [5.894485378547998e-01 1.287131228874161e+00]
DelayMargin: [4.960154377113239e+00 7.883892905655167e-01]
DMFrequency: [5.894485378547998e-01 1.287131228874161e+00]
Stable: 1
>> realcrossing1 = -1./s.GainMargin % crossings on the negative real axis
realcrossing1 =
-4.444007591691567e-01
>> s=allmargin(-h)
s =
struct with fields:
GainMargin: [3.000000000000005e+00 1.186263434433761e+00]
GMFrequency: [0 4.858495157107252e-01]
PhaseMargin: [-1.248113547729714e+01 -1.218585074857773e+02]
PMFrequency: [5.894485378547998e-01 1.287131228874161e+00]
DelayMargin: [1.028986927474128e+01 3.229160350356917e+00]
DMFrequency: [5.894485378547998e-01 1.287131228874161e+00]
Stable: 1
>> realcrossing2 = 1./s.GainMargin % crossings on the postive real axis
realcrossing2 =
3.333333333333328e-01 8.429830769228168e-01
As shown, these crossings are for positive (and zero) frequencies. By symmetry, the same crossings should occur for negative frequencies (don't count w = 0 or w=inf twice). Presumably you don't care about crossings on the real axis at infinity due to poles on the imaginary axis.

Community Treasure Hunt

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

Start Hunting!