Bode plots, Nyquist plots, ans stability margins¶
$$ \newcommand{\eg}{{\it e.g.}} \newcommand{\ie}{{\it i.e.}} \newcommand{\argmin}{\operatornamewithlimits{argmin}} \newcommand{\mc}{\mathcal} \newcommand{\mb}{\mathbb} \newcommand{\mf}{\mathbf} \newcommand{\minimize}{{\text{minimize}}} \newcommand{\diag}{{\text{diag}}} \newcommand{\cond}{{\text{cond}}} \newcommand{\rank}{{\text{rank }}} \newcommand{\range}{{\mathcal{R}}} \newcommand{\null}{{\mathcal{N}}} \newcommand{\tr}{{\text{trace}}} \newcommand{\dom}{{\text{dom}}} \newcommand{\dist}{{\text{dist}}} \newcommand{\R}{\mathbf{R}} \newcommand{\SM}{\mathbf{S}} \newcommand{\ball}{\mathcal{B}} \newcommand{\bmat}[1]{\begin{bmatrix}#1\end{bmatrix}} $$ ASE3093: Automatic control, Inha University.
Jong-Han Kim (jonghank@inha.ac.kr)
In [1]:
!pip install control
Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/ Requirement already satisfied: control in /usr/local/lib/python3.10/dist-packages (0.9.3.post2) Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from control) (1.22.4) Requirement already satisfied: scipy>=1.3 in /usr/local/lib/python3.10/dist-packages (from control) (1.10.1) Requirement already satisfied: matplotlib in /usr/local/lib/python3.10/dist-packages (from control) (3.7.1) Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib->control) (1.0.7) Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib->control) (0.11.0) Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib->control) (4.39.3) Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib->control) (1.4.4) Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib->control) (23.1) Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib->control) (8.4.0) Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib->control) (3.0.9) Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib->control) (2.8.2) Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.7->matplotlib->control) (1.16.0)
In [2]:
import numpy as np
import matplotlib.pyplot as plt
from control.matlab import *
In [3]:
omega = 1
zeta = 0.5
num = np.array([omega**2])
den = np.array([1, 2*zeta*omega, omega**2])
G = tf(num, den)
r = 15
bode(G, Hz=False)
fig = plt.gcf().set_size_inches(8, 8)
plt.show()
In [4]:
num = np.array([1, 2])
den = np.array([1, 0, 0])
G = tf(num, den)
r = 10
rlocus(G)
fig = plt.gcf().set_size_inches(8, 8)
plt.xlim([-r, r]), plt.ylim([-r, r])
plt.grid()
plt.show()
In [5]:
bode(G, Hz=False)
fig = plt.gcf().set_size_inches(8, 8)
plt.show()
In [6]:
r = 4
nyquist(G)
fig = plt.gcf().set_size_inches(8, 8)
plt.xlim([-r, r]), plt.ylim([-r, r])
plt.show()
r = 40
nyquist(G)
fig = plt.gcf().set_size_inches(8, 8)
plt.xlim([-r, r]), plt.ylim([-r, r])
plt.show()
In [7]:
pole(feedback(G,4,-1))
Out[7]:
array([-2.+2.j, -2.-2.j])
In [8]:
bode(4*G, Hz=False)
fig = plt.gcf().set_size_inches(8, 8)
plt.show()
In [9]:
num = np.array([15])
den = np.array([1, 6, 5, 0])
G = tf(num, den)
r = 10
rlocus(G)
fig = plt.gcf().set_size_inches(8, 8)
plt.xlim([-r, r]), plt.ylim([-r, r])
plt.grid()
plt.show()
In [10]:
bode(G, Hz=False)
fig = plt.gcf().set_size_inches(8, 8)
plt.show()
In [11]:
nyquist(G)
r = 4
fig = plt.gcf().set_size_inches(8, 8)
plt.xlim([-r, r]), plt.ylim([-r, r])
plt.show()
nyquist(G)
r = 160
fig = plt.gcf().set_size_inches(8, 8)
plt.xlim([-r, r]), plt.ylim([-r, r])
plt.show()
In [12]:
margin(G)
Out[12]:
(2.0000000000000004, 15.552687434223344, 2.23606797749979, 1.5519216327158523)
In [13]:
def plot_margins(sys):
# mag,phase,omega = bode(sys,dB=True,Plot=False)
mag,phase,omega = bode(sys,dB=True,plot=False)
magdB = 20*np.log10(mag)
phase_deg = phase*180.0/np.pi
Gm,Pm,Wcg,Wcp = margin(sys)
GmdB = 20*np.log10(Gm)
##Plot Gain and Phase
f,(ax1,ax2) = plt.subplots(2,1)
ax1.semilogx(omega,magdB)
ax1.grid(which="both")
ax1.set_xlabel('Frequency (rad/s)')
ax1.set_ylabel('Magnitude (dB)')
ax2.semilogx(omega,phase_deg)
ax2.grid(which="both")
ax2.set_xlabel('Frequency (rad/s)')
ax2.set_ylabel('Phase (deg)')
ax1.set_title('Gm = '+str(np.round(GmdB,2))+' dB (at '+str(np.round(Wcg,2))+' rad/s), Pm = '+str(np.round(Pm,2))+' deg (at '+str(np.round(Wcp,2))+' rad/s)')
###Plot the zero dB line
ax1.plot(omega,0*omega,'k--',linewidth=2)
###Plot the -180 deg lin
ax2.plot(omega,-180+0*omega,'k--',linewidth=2)
##Plot the vertical line from -180 to 0 at Wcg
ax2.plot([Wcg,Wcg],[-180,0],'r--',linewidth=2)
##Plot the vertical line from -180+Pm to 0 at Wcp
ax2.plot([Wcp,Wcp],[-180+Pm,0],'g--',linewidth=2)
##Plot the vertical line from min(magdB) to 0-GmdB at Wcg
ax1.plot([Wcg,Wcg],[np.min(magdB),0-GmdB],'r--',linewidth=2)
##Plot the vertical line from min(magdB) to 0db at Wcp
ax1.plot([Wcp,Wcp],[np.min(magdB),0],'g--',linewidth=2)
return Gm,Pm,Wcg,Wcp
In [14]:
K = zpk([-2],[-10],5)
bode(G,G*K, Hz=False)
fig = plt.gcf().set_size_inches(8, 8)
plt.show()
In [16]:
plot_margins(G)
plot_margins(G*K)
Out[16]:
(6.969340028069279, 40.90168542948257, 5.982790434368795, 1.8102803080172076)
In [17]:
G2 = zpk([],[-1,-2,-10],20)
K2 = zpk([-1],[-0.01],1)
bode(G2,G2*K2, Hz=False)
fig = plt.gcf().set_size_inches(8, 8)
plt.show()
