NiiL

W11_StabilityOnFrequencyResponse

English translation

!pip install control
Requirement already satisfied: control in /usr/local/lib/python3.10/dist-packages (0.10.0)
Requirement already satisfied: numpy>=1.23 in /usr/local/lib/python3.10/dist-packages (from control) (1.25.2)
Requirement already satisfied: scipy>=1.8 in /usr/local/lib/python3.10/dist-packages (from control) (1.11.4)
Requirement already satisfied: matplotlib>=3.6 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>=3.6->control) (1.2.1)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.6->control) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.6->control) (4.51.0)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.6->control) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.6->control) (24.0)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.6->control) (9.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.6->control) (3.1.2)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.6->control) (2.8.2)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.7->matplotlib>=3.6->control) (1.16.0)

import numpy as np
import control as ctl
import matplotlib.pyplot as plt

Nyquist Criterion

Z=P+NZ = P + N

ZZ: Stable if Z=0Z = 0 \ PP: # of right half plane poles \ NN: # of circulations around (−1,0)(-1, 0) \ \;\;\;(CW: ++, CCW: −-)

G(s)=K(s2+10s+24)s2−8s+15G(s) = \frac{K(s^2 + 10s + 24)}{s^2 - 8s + 15}

Right half plane poles:

s=ctl.tf('s')
G = (s**2 + 10*s + 24) / (s**2 - 8*s + 15)
poles = ctl.poles(G)
poles
array([5.+0.j, 3.+0.j])
P = len(poles)
P
2
N = 0 - P
N
-2

For K=1K=1,

K = 1
result = ctl.nyquist(K*G)
<Figure size 640x480 with 1 Axes>
실행 결과

N=−2N = -2 and 2 CCW circulations around (−1,0)(-1,0) →\rightarrow Stable

t,y = ctl.step_response(ctl.feedback(K*G,1), np.arange(0,10, 0.1))
f1 = plt.plot(t,y)
<Figure size 640x480 with 1 Axes>
실행 결과

For K=0.5K = 0.5,

K = 0.5
result = ctl.nyquist(K*G)
<Figure size 640x480 with 1 Axes>
실행 결과

N=−2N = -2 and 0 CCW circulations around (−1,0)(-1,0) →\rightarrow Unstable

t,y = ctl.step_response(ctl.feedback(K*G,1), np.arange(0,10, 0.1))
f2 = plt.plot(t,y)
<Figure size 640x480 with 1 Axes>
실행 결과

Gain and phase margin

G(s)=K50s3+9s2+30s+40G(s) = K\frac{50}{s^3 + 9s^2 + 30s + 40}

s = ctl.tf('s')
G = 50/(s**3 + 9*s**2 + 30*s + 40)
G.name='G(s)'
K = 5.0
KG = K*G
KG.name = 'K G(s)'
out = ctl.bode([G, KG])
gm_1, pm_1, Wcg_1, Wcp_1 = ctl.margin(G)
gm_2, pm_2, Wcg_2, Wcp_2 = ctl.margin(KG)
print('For K=1, GM = %.2f dB (at %.2f rad/s) and PM = %.2f dB (at %.2f rad/s)' % (gm_1, Wcg_1, pm_1, Wcp_1))
print('For K=5, GM = %.2f dB (at %.2f rad/s) and PM = %.2f dB (at %.2f rad/s)' % (gm_2, Wcg_2, pm_2, Wcp_2))
For K=1, GM = 4.60 dB (at 5.48 rad/s) and PM = 100.66 dB (at 1.85 rad/s)
For K=5, GM = 0.92 dB (at 5.48 rad/s) and PM = -2.85 dB (at 5.67 rad/s)

<Figure size 640x480 with 2 Axes>
실행 결과
t,y = ctl.step_response(ctl.feedback(G,1), np.arange(0,10, 0.1))
f3 = plt.plot(t,y)
<Figure size 640x480 with 1 Axes>
실행 결과
t,y = ctl.step_response(ctl.feedback(KG,1), np.arange(0,10, 0.1))
f4 = plt.plot(t,y)
<Figure size 640x480 with 1 Axes>
실행 결과