Help with creating my own Nyquist plotting function

1 visualización (últimos 30 días)
Daniel Beeson
Daniel Beeson el 24 de Mayo de 2020
Respondida: Daniel Beeson el 24 de Mayo de 2020
Hi everyone,
I am trying to make a function that takes a transfer function L(s) in the form of an array of numerator and denominator coefficients and spits out the nyquist plot, avoiding imaginary poles by drawing a semi-circle on the right half plane with radius epsilon around those imaginary poles and ending with an s-plane contour that is a semi-circle on the right half plane with radius R (see attached image)
So far, this is my code:
function n = nyquill(N,D,R,epsilon)
%N is numerator coefficients, D is denominator coefffs
L = tf(N,D);
t = pole(L)
r = NaN(3,1);
d = zeros(3,1);
for n = 1:size(t)
if real(t(n)) == 0
t(n) = r(n);
d(n) = NaN;
else
t(n) = d(n);
r(n) = NaN;
end
end
if r(1)==NaN & r(2)==NaN & r(3)==NaN
g1 = [-R:0.1:R]*i;
g2 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2];
elseif r(1)==NaN & r(2)==NaN
g1 = [-R;0.1:imag(r(3))-epsilon]*i;
g2 = imag(r(3)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(3)) + epsilon:.01:R]*i;
g4 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4];
elseif r(1)==NaN
g1 = [-R;0.1:imag(r(2))-epsilon]*i;
g2 = imag(r(2)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(2)) + epsilon:.01:imag(r(3))-epsilon]*i;
g4 = imag(r(3)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g5 = [imag(r(3)) + epsilon:.01:R]*i;
g6 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4 g5 g6];
elseif r(2)==NaN
g1 = [-R;0.1:imag(r(1))-epsilon]*i;
g2 = imag(r(1)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(1)) + epsilon:.01:imag(r(3))-epsilon]*i;
g4 = imag(r(3)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g5 = [imag(r(3)) + epsilon:.01:R]*i;
g6 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4 g5 g6];
elseif r(3)==NaN
g1 = [-R;0.1:imag(r(1))-epsilon]*i;
g2 = imag(r(1)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(1)) + epsilon:.01:imag(r(2))-epsilon]*i;
g4 = imag(r(2)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g5 = [imag(r(2)) + epsilon:.01:R]*i;
g6 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4 g5 g6];
elseif r(2)==NaN & r(3)==NaN
g1 = [-R;0.1:imag(r(1))-epsilon]*i;
g2 = imag(r(1)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(1)) + epsilon:.01:R]*i;
g4 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4];
elseif r(1)==NaN & r(3)==NaN
g1 = [-R;0.1:imag(r(2))-epsilon]*i;
g2 = imag(r(2)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(2)) + epsilon:.01:R]*i;
g4 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4];
elseif r(1)~=NaN & r(2)~=NaN & r(3)~=NaN
g1 = [-R;0.1:imag(r(1))-epsilon]*i;
g2 = imag(r(1)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g3 = [imag(r(1)) + epsilon:.01:imag(r(2))-epsilon]*i;
g4 = imag(r(2)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g5 = [imag(r(2)) + epsilon:.01:imag(r(3))-epsilon]*i;
g6 = imag(r(3)) + epsilon*exp(i*[-pi/2:pi/90:pi/2]);
g7 = [imag(r(3)) + epsilon:.01:R]*i;
g8 = R*exp(i*[pi/2:-pi/90:-pi/2]);
g = [g1 g2 g3 g4 g5 g6 g7 g8];
end
c = poly2sym(N);
d = poly2sym(D);
ln = c/d;
x = g;
subs(ln);
plot(real(ln),imag(ln));grid on;axis('equal')
hold on
plot(-1,0, '-o')
end
And I end up receiving the error:
Error using plot
Data must be numeric, datetime, duration or an array convertible to double.
Error in nyquill (line 83)
plot(real(ln),imag(ln));grid on;axis('equal')
I don't know what this error means or exactly how to fix it, any help would be appreciated!
  2 comentarios
Daniel Beeson
Daniel Beeson el 24 de Mayo de 2020
The inputs N and D are the numerator and denominator coefficients
Daniel Beeson
Daniel Beeson el 24 de Mayo de 2020
The idea is to avoid the poles that lie on the imaginary axis by going around them with a small semi circle of radius epsilon that only goes into the right half plane

Iniciar sesión para comentar.

Respuesta aceptada

Daniel Beeson
Daniel Beeson el 24 de Mayo de 2020
Hey everyone, I figured it out. My conditions for NaN were incorrect. instead of using ==NaN, I changed it to isnan(r(...))==1 for if the value of r was NaN or ==0 for if it is not. Thank you everyone for the help!

Más respuestas (1)

rubindan
rubindan el 24 de Mayo de 2020

Categorías

Más información sobre Scatter Plots en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by