71 views (last 30 days)

Show older comments

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!

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.

Paul
on 12 Feb 2021

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.

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

Start Hunting!