04. Practical Case Study: Frequency-Domain Deep Dive & Nyquist Stability Analysis¶
Bode Margins, Cauchy's Argument Principle, and Indented Integrator Contours
In this advanced case study, we dive deeply into frequency-domain stability criteria. We analyze open-loop stability margins ($GM, PM, \omega_{cg}, \omega_{cp}$), investigate Nyquist stability encirclements with indented contours around $s=0$, and correlate frequency-domain metrics with time-domain closed-loop stability.
import matplotlib.pyplot as plt
import numpy as np
import ctrlpy as cp
from ctrlpy.plotting_plotly import (
plot_bode_plotly,
plot_nyquist_plotly,
plot_root_locus_plotly,
)
%matplotlib inline
1. Foundations of Frequency-Domain Stability¶
Consider an open-loop plant and controller in negative feedback: $$T(s) = \frac{L(s)}{1 + L(s)}, \quad \text{where } L(s) = C(s) G(s)$$
The closed-loop poles are the roots of the characteristic equation: $$1 + L(s) = 0 \iff L(s) = -1 + 0j$$
Nyquist Stability Criterion (Cauchy's Argument Principle)¶
$$Z = N + P$$ where:
- $P$: Number of unstable open-loop poles (in right-half plane $\mathrm{Re}(s) > 0$).
- $N$: Number of clockwise encirclements of the critical point $-1 + 0j$ in the complex plane as $s$ traverses the Nyquist contour.
- $Z$: Number of unstable closed-loop poles. For asymptotic stability, we must have $Z = 0 \implies N = -P$ (or $N = 0$ if open-loop is stable).
2. Study System & Analytical Margins¶
We analyze the canonical third-order system with an integrator: $$L(s) = \frac{K}{s (s + 1)(s + 2)} = \frac{K}{s^3 + 3s^2 + 2s}$$
Analytical Margin Derivation:¶
Phase Crossover Frequency $\omega_{cp}$: $$\angle L(j\omega) = -90^\circ - \arctan(\omega) - \arctan(\omega/2)$$ Phase reaches $-180^\circ$ when $\arctan(\omega) + \arctan(\omega/2) = 90^\circ$: $$\tan(90^\circ) = \frac{\omega + \omega/2}{1 - \omega^2/2} = \infty \implies 1 - \frac{\omega^2}{2} = 0 \implies \omega_{cp} = \sqrt{2} \approx 1.4142\text{ rad/s}$$
Critical Gain $K_{crit}$: $$|L(j\sqrt{2})| = \frac{K}{\sqrt{2} \cdot |j\sqrt{2} + 1| \cdot |j\sqrt{2} + 2|} = \frac{K}{\sqrt{2} \cdot \sqrt{3} \cdot \sqrt{6}} = \frac{K}{6}$$ Setting $|L(j\sqrt{2})| = 1 \implies \mathbf{K_{crit} = 6.0}$.
# Define system for three distinct gain regimes
K_sub = 2.0 # Subcritical (Gain < 6): Asymptotically Stable
K_crit = 6.0 # Critical (Gain = 6): Marginally Stable (Undamped Oscillation)
K_super = 12.0 # Supercritical (Gain > 6): Unstable (Diverging Oscillation)
L_sub = cp.tf([K_sub], [1.0, 3.0, 2.0, 0.0])
L_crit = cp.tf([K_crit], [1.0, 3.0, 2.0, 0.0])
L_super = cp.tf([K_super], [1.0, 3.0, 2.0, 0.0])
# Compute exact stability margins
sm_sub = cp.margin(L_sub)
sm_crit = cp.margin(L_crit)
sm_super = cp.margin(L_super)
print("=== Stability Margins Summary ===")
print(
f"K = 2.0 (Subcritical) : GM = {sm_sub.gm_db:+.2f} dB, PM = {sm_sub.pm_deg:+.2f}°, wcg = {sm_sub.wcg:.3f} rad/s, wcp = {sm_sub.wcp:.3f} rad/s"
)
print(
f"K = 6.0 (Critical) : GM = {sm_crit.gm_db:+.2f} dB, PM = {sm_crit.pm_deg:+.2f}°, wcg = {sm_crit.wcg:.3f} rad/s, wcp = {sm_crit.wcp:.3f} rad/s"
)
print(
f"K = 12.0 (Supercritical): GM = {sm_super.gm_db:+.2f} dB, PM = {sm_super.pm_deg:+.2f}°, wcg = {sm_super.wcg:.3f} rad/s, wcp = {sm_super.wcp:.3f} rad/s"
)
3. Bode Diagram Analysis with Margin Annotations¶
# Matplotlib Bode comparison
fig, (ax_mag, ax_phase) = cp.plot_bode(L_sub, margins=True)
ax_mag.set_title("Bode Diagram: K = 2 (Stable, GM = +9.54 dB)")
plt.show()
# Interactive Plotly Bode plot for the critical system K = 6
fig_bode_crit = plot_bode_plotly(L_crit, margins=True)
fig_bode_crit.show()
4. Nyquist Diagram & Integrator Indented Contour Analysis¶
When an open-loop system has poles on the imaginary axis (such as the pole at $s = 0$), the standard Nyquist contour is indented into the right-half plane along a small semicircular arc: $$s = \epsilon\, e^{j\theta}, \quad \theta \in [-\pi/2, +\pi/2]$$
The mapped response $L(\epsilon e^{j\theta}) \approx \frac{K}{2\epsilon} e^{-j\theta}$ sweeps an arc of infinite radius in the right complex plane from $+90^\circ$ through $0^\circ$ to $-90^\circ$.
ctrlpy automatically detects imaginary axis poles and computes this indented arc:
# Compute Nyquist data for subcritical and critical systems
nyq_sub = cp.nyquist_data(L_sub)
nyq_crit = cp.nyquist_data(L_crit)
nyq_super = cp.nyquist_data(L_super)
print("Nyquist Arc Response Detected:")
print(f"Arc points generated: {len(nyq_sub.arc_s) if nyq_sub.arc_s is not None else 0}")
# Interactive Nyquist diagram with critical point (-1, 0)
fig_nyq = plot_nyquist_plotly(L_sub)
fig_nyq.show()
# Interactive Nyquist diagram for supercritical system K = 12
fig_nyq_super = plot_nyquist_plotly(L_super)
fig_nyq_super.show()
5. Corroboration: Root Locus & Closed-Loop Time Response¶
Let's cross-verify our frequency-domain conclusions with the Root Locus trajectory and time-domain simulations:
# Root Locus showing the jw-axis crossing at w = sqrt(2) ~= 1.414 rad/s
fig_rl = plot_root_locus_plotly(L_sub)
fig_rl.show()
# Form closed-loop systems
T_sub = cp.feedback(L_sub, 1.0)
T_crit = cp.feedback(L_crit, 1.0)
T_super = cp.feedback(L_super, 1.0)
print(f"Closed-loop poles (K=2) : {T_sub.poles()}")
print(f"Closed-loop poles (K=6) : {T_crit.poles()}")
print(f"Closed-loop poles (K=12): {T_super.poles()}")
# Simulate closed-loop step response
t_sim = np.linspace(0.0, 15.0, 1000)
resp_sub = cp.step_response(T_sub, T=t_sim)
resp_crit = cp.step_response(T_crit, T=t_sim)
resp_super = cp.step_response(T_super, T=t_sim)
plt.figure(figsize=(10, 5))
plt.plot(t_sim, resp_sub.y, "b-", lw=2, label="K = 2.0 (Stable, Damped)")
plt.plot(
t_sim,
resp_crit.y,
"g--",
lw=2,
label=r"K = 6.0 (Marginally Stable, Undamped $\omega=\sqrt{2}$)",
)
plt.plot(t_sim, resp_super.y, "r-.", lw=2, label="K = 12.0 (Unstable, Diverging)")
plt.axhline(1.0, color="k", linestyle=":", label="Command Input")
plt.title("Closed-Loop Step Response Across Stability Regimes")
plt.xlabel("Time [s]")
plt.ylabel("Output $y(t)$")
plt.ylim(-1.5, 3.5)
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()
Key Takeaways¶
- Exact stability margins ($GM, PM, \omega_{cg}, \omega_{cp}$) quantify the gain and phase delay the system can tolerate before becoming unstable.
- The critical gain $K_{crit} = 6.0$ corresponds to $GM = 0\text{ dB}$ and $PM = 0^\circ$ at the phase crossover frequency $\omega_{cp} = \sqrt{2}\text{ rad/s}$.
- Nyquist diagrams correctly handle open-loop poles on the imaginary axis via indented arc contour mapping.
- Frequency-domain margins, Root Locus imaginary axis crossings, and time-domain transient responses all reinforce the exact same physical stability boundaries.