W11_StabilityOnFrequencyResponse

!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 + N$$

$$G(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=1$,

K = 1
result = ctl.nyquist(K*G)

<Figure size 640x480 with 1 Axes>

$N = -2$ and 2 CCW circulations around $(-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.5$,

K = 0.5
result = ctl.nyquist(K*G)

<Figure size 640x480 with 1 Axes>

$N = -2$ and 0 CCW circulations around $(-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) = 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>