TODO: Redraw diagrams using Python code. TODO: Explanation for setting up animation function.
By the end of this notebook, you should be able to:
# import libraries we need later
import control
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
# https://stackoverflow.com/questions/49722298/drawing-a-nyquist-diagram-animation-on-python
from matplotlib import animation
from IPython.display import HTML
from matplotlib.patches import Arrow
Nyquist plots and stability is an important topic in classical control, which is usually taught as a "bonus" section in undergraduate CHBE control courses. However, it is part of the regular curriculum in MECH/ECE undergrad control, so I consider it necessary to cover it here also in CHBE, at least for those who wish to have a strong background in control.
To motivate the use of Nyquist, here are some of its strengths over other methods:
With that, let's look at the well-known, dreaded -1 point and see how it relates to closed-loop instability:
# TODO replace diagram
Above is an empty Nyquist plot with the blue line representing the imaginary axis, and the red arrow representing a vector from the origin to the point $-1$. Suppose that this vector represents a certain open-loop TF, $L(s)$, at some specific frequency $\omega$. Using simple geometry, we can determine its magnitude and phase as $\lvert L(j \omega) \rvert = 1 = 0dB$ and $\angle L(j \omega) = -180^o$, respectively.
Recall from Bode that this combination implies $L$ is at its ultimate gain, or at the point of marginal stability. In other words, if any slight extra gain or delay is introduced to $L$, the closed-loop would become unstable. Equivalently, at this $-1$ point, the closed-loop TF $\big( \frac{L}{1+L} \big)$ has a pair of poles located at $\pm j \omega_{180}$, which again implies marginal stability. In light of these realizations, as control engineers we strive to stay as far away from the $-1$ point as possible, so as to allow ourselves comfortably large gain and phase margins when designing controllers for our systems.
However, at this point we still haven't justified mathematically why we're concerned with the point $-1$ (and not, for example, $-2$ or $-5$) and where the equation $Z = N + P$ comes from.
To do this, let's discuss the following 2 concepts:
Tracings or mappings of a rational expression on the complex plane, by cycling its argument (in this case frequency) between $- \infty$ and $+\infty$. For example, the expression $F(s) = \frac{1}{s+1}$ has the spectrum $F(j \omega) = \frac{1-j \omega}{1+\omega^2}$ which can be evaluated at the following frequencies:
$\omega$ | $F(j\omega)$ |
---|---|
$-\infty$ | $0+0^+j$ |
$-1$ | $0.5+0.5j$ |
$0$ | $1+0j$ |
$1$ | $0.5-0.5j$ |
$\infty$ | $0+0^-j$ |
Note: Take a look at the Markdown code for the table above. The hidden <img width=200/>
in the last row is a cool trick to size the table.
Knowing these points, the contour of F can be traced, and software such as MATLAB can do this relatively quickly:
We can also do this easily in Python:
# Set up the plot parameters
plt.figure(figsize=(10,5))
plt.grid(True)
plt.title('Nyquist Diagram')
plt.xlabel('Re(s)')
plt.ylabel('Im(s)')
# Nyquist plot
sys = control.tf(1,[1,1])
omega = np.logspace(-3,3,60)
real_arr, imag_arr, freq_arr = control.nyquist_plot(sys,omega)
plt.show()
And even create animations for it!
def animate_Nyquist(sys,omega,xlims,ylims):
# Set up figure params
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(1,1,1)
ax.grid(True)
ax.set_title('Nyquist Diagram')
ax.set_xlabel('Re(s)')
ax.set_ylabel('Im(s)')
# Get the magnitude and phase of the system
mag_tmp, phase_tmp, omega = sys.freqresp(omega)
mag = np.squeeze(mag_tmp)
phase = np.squeeze(phase_tmp)
# Compute the primary curve
# https://www.youtube.com/watch?v=HbuGlYEoYEA
x = np.hstack([mag*sp.cos(phase),np.flip(mag*sp.cos(phase),axis=0)])
y = np.hstack([mag*sp.sin(phase),-mag*sp.sin(phase)]) # negative to get the mirror image
xData, yData = [], []
primary_line, = ax.plot([], [], '-b', animated=True)
primary_arrow = ax.annotate('', xy = (0,0),
xytext = (0,0),
size=15,
arrowprops = {'arrowstyle': "->"})
# Mark the -1 point
ax.plot([-1], [0], 'r+')
def init():
ax.set_xlim(xlims)
ax.set_ylim(ylims)
return (primary_line,primary_arrow)
def animate(i):
xData.append(x[i])
yData.append(y[i])
# plot the actual line
primary_line.set_data(xData, yData)
primary_arrow.set_position((x[i],y[i]))
# This is a crude way to handle the last data point
if i < len(x)-1:
primary_arrow.xy = (x[i+1],y[i+1])
else:
primary_arrow.xy = (x[i],y[i])
# Print out the current frame using "\r", the 'carriage return' character, as our end character.
# This makes Python print the frame on the same line.
print("Frame: {:0d}".format(i), end="\r")
return (primary_line,primary_arrow)
return fig, animate, init, x
sys = control.tf(1,[1,1])
omega = np.logspace(-3,3,60)
xlims = [-1.1, 1.1]
ylims = [-0.55, 0.55]
fig, animate, init, x = animate_Nyquist(sys,omega,xlims,ylims)
# blit=True will only draw the parts that have changed.
# The `interval` parameter is the delay between frames in milliseconds and it controls the speed of the animation.
anim = animation.FuncAnimation(fig, animate, init_func=init, frames=range(len(x)), interval=40, blit=True)
HTML(anim.to_jshtml())
Frame: 119
We can take any rational transfer function and apply the previous procedure to sketch its contour. Knowing this, let's move on to the next important concept, which explains the reason behind the $-1$ point as well as the $Z = N+P$ equation.
The contour of any rational expression will encircle the origin once in a clockwise fashion, for each additional RHP-zero it has compared to RHP-poles. Mathematically, this is expressed as:
$$ N = Z - P $$where $N$ is the number of clockwise encirclements the contour makes about the origin, $Z$ the number of RHP-zeroes in $F$, and $P$ the number of RHP-poles in $F$. This may look confusing, so let's go through several examples as exercise:
$F(s)=\frac{s-1}{s+1}$: $N=Z-P=1-0=1$, So we expect its contour to encircle the origin once clockwise.
sys = control.tf([1,-1],[1,1])
omega = np.logspace(-3,3,60)
xlims = [-1.1, 1.1]
ylims = [-1, 1]
fig, animate, init, x = animate_Nyquist(sys,omega,xlims,ylims)
anim = animation.FuncAnimation(fig, animate, init_func=init, frames=range(len(x)), interval=40, blit=True)
HTML(anim.to_jshtml())
Frame: 119
IOPub data rate exceeded. The notebook server will temporarily stop sending output to the client in order to avoid crashing it. To change this limit, set the config variable `--NotebookApp.iopub_data_rate_limit`.
$F(s)=\frac{(s-2)(s-3)}{(s-1)(s-4)}$: $N=Z-P=2-2=0$, so we expect no clockwise encirclements of the origin.
s = control.tf([1,0],[1])
sys = (s-2)*(s-3)/((s-1)*(s-4))
omega = np.logspace(-3,3,60)
xlims = [-1.1, 1.5]
ylims = [-0.3, 0.3]
fig, animate, init, x = animate_Nyquist(sys,omega,xlims,ylims)
anim = animation.FuncAnimation(fig, animate, init_func=init, frames=range(len(x)), interval=40, blit=True)
HTML(anim.to_jshtml())
Frame: 119
$F(s)=\frac{(s-15)(s-0.8)}{(s+1.1)(s+4.7)}$: $N=Z-P=2-0=2$, so we expect two clockwise encirclements of the origin.
TODO - fix code for mirror, should be:
s = control.tf([1,0],[1])
sys = (s-15)*(s-0.8)/((s+1.1)*(s+4.7))
omega = np.logspace(-3,3,60)
xlims = [-3, 3]
ylims = [-3, 3]
fig, animate, init, x = animate_Nyquist(sys,omega,xlims,ylims)
anim = animation.FuncAnimation(fig, animate, init_func=init, frames=range(len(x)), interval=40, blit=True)
HTML(anim.to_jshtml())
Frame: 119
$F(s)=\frac{s+1}{(s-2)(s-3)}$: $N=Z-P=0-2=-2$; since $N$ is negative, we expect two counterclockwise encirclements of the origin, which is again confirmed.
TODO fix code
s = control.tf([1,0],[1])
sys = (s+1)/((s-2)*(s-3))
omega = np.logspace(-2,2,60)
xlims = [-2.5, 2]
ylims = [-2.5, 2.5]
fig, animate, init, x = animate_Nyquist(sys,omega,xlims,ylims)
anim = animation.FuncAnimation(fig, animate, init_func=init, frames=range(len(x)), interval=40, blit=True)
HTML(anim.to_jshtml())
Frame: 119
So Cauchy's argument principle turns out to be pretty straightforward: it tells us how many more RHP-zeroes than RHP-poles a TF has.
But so far we've discussed encirclements about the origin; how does this translate to the well-known $-1$ point in Nyquist stability? Usually profound realizations (or proofs) are made by connecting a series of trivial statements, and this is no exception.
Consider the following equivalent statements:
The final statement above completes the proof for Nyquist stability.
As a last exercise, let's visually deduce the Gain and Phase Margins on a Nyquist plot, using our understanding from Bode and simple geometry:
In the previous figure, the black solid curve represents the contour of $L(j\omega)$. At the point 1, $\angle L(j\omega)=-180^o$, therefore the corresponding frequency is $\omega = \omega_{180}$. The length of the vector spanning from the origin to the point 1 is $\lvert L(j\omega_{180})\rvert$, which is the inverse of the Gain Margin.
On the other hand, the point 2 is the intersection between the contour of $L(j\omega)$ and the unit circle, so $\lvert L(j\omega)\rvert = 1$ here, therefore the corresponding frequency is $\omega=\omega_{0dB}$. The Phase Margin is therefore the difference between the phase at this point and $-180^o$, i.e. $\angle L(j\omega_{0dB})-(-180^o)$.
From the figure above, we can see that both the Gain and Phase Margins can be geometrically interpreted as ``distances" from the $-1$ point, which we try to stay as far away from as possible.