Pedagogical Case Study: Routh-Hurwitz Stability Criterion¶

The Routh-Hurwitz Stability Criterion is a classical algebraic tool used to determine how many roots of a characteristic polynomial lie in the open Right-Half Plane (RHP) without having to factor the polynomial.

In this notebook, we demonstrate how ctrlpy.symbolic.routh.RouthArray can be used to:

1. Construct complete, analytical Routh arrays for n-th order polynomials.
2. Automatically handle textbook edge cases:
   • Special Case 1: Zero in the first column (epsilon > 0 substitution).
   • Special Case 2: Row of all zeros (Auxiliary polynomial A(s) and dA/ds).
3. Solve for closed-loop parametric gain boundaries (K_min < K < K_max).
4. Generate step-by-step pedagogical explanations for classroom lecture notes and homework grading.

In [ ]:

import matplotlib.pyplot as plt
import numpy as np

import ctrlpy as cp
from ctrlpy.symbolic import RouthArray

print(f"ctrlpy version: {cp.__version__}")

import matplotlib.pyplot as plt import numpy as np import ctrlpy as cp from ctrlpy.symbolic import RouthArray print(f"ctrlpy version: {cp.__version__}")

1. Standard Stable Polynomial¶

Consider a 3rd-order system with characteristic polynomial:

$$P(s) = s^3 + 2s^2 + 3s + 4 = 0$$

In [ ]:

# Construct Routh array
ra_stable = RouthArray([1, 2, 3, 4])

print("ASCII Representation:")
print(ra_stable)
print(f"\nStrictly Stable: {ra_stable.is_stable}")
print(f"Number of RHP Poles: {ra_stable.num_rhp_poles}")

# Jupyter rich LaTeX representation
ra_stable

# Construct Routh array ra_stable = RouthArray([1, 2, 3, 4]) print("ASCII Representation:") print(ra_stable) print(f"\nStrictly Stable: {ra_stable.is_stable}") print(f"Number of RHP Poles: {ra_stable.num_rhp_poles}") # Jupyter rich LaTeX representation ra_stable

2. Unstable Polynomial with Sign Changes¶

Consider the characteristic polynomial:

$$P(s) = s^3 + s^2 + 2s + 24 = 0$$

The first column signs are [+1, +1, -1, +1], showing 2 sign changes (2 RHP poles).

In [ ]:

ra_unstable = RouthArray([1, 1, 2, 24])

print(f"Strictly Stable: {ra_unstable.is_stable}")
print(f"Number of RHP Poles: {ra_unstable.num_rhp_poles}")

ra_unstable

ra_unstable = RouthArray([1, 1, 2, 24]) print(f"Strictly Stable: {ra_unstable.is_stable}") print(f"Number of RHP Poles: {ra_unstable.num_rhp_poles}") ra_unstable

3. Special Case 1: First Element in a Row is Zero¶

Consider the 5th-order polynomial:

$$P(s) = s^5 + 2s^4 + 2s^3 + 4s^2 + 11s + 10 = 0$$

In row s^3, the first element is zero, but subsequent entries are non-zero. RouthArray automatically substitutes epsilon > 0.

In [ ]:

ra_eps = RouthArray([1, 2, 2, 4, 11, 10])

print(f"Strictly Stable: {ra_eps.is_stable}")
print(f"Number of RHP Poles: {ra_eps.num_rhp_poles}")
print("\nConstruction notes:")
for note in ra_eps.steps:
    print(f"- {note}")

ra_eps

ra_eps = RouthArray([1, 2, 2, 4, 11, 10]) print(f"Strictly Stable: {ra_eps.is_stable}") print(f"Number of RHP Poles: {ra_eps.num_rhp_poles}") print("\nConstruction notes:") for note in ra_eps.steps: print(f"- {note}") ra_eps

4. Special Case 2: Row of All Zeros & Auxiliary Polynomials¶

Consider the polynomial:

$$P(s) = s^3 + 2s^2 + 4s + 8 = (s + 2)(s^2 + 4) = 0$$

All entries in row s^1 vanish. RouthArray extracts the Auxiliary Polynomial A(s) = 2s^2 + 8, differentiates dA/ds = 4s, and completes the array.

In [ ]:

ra_aux = RouthArray([1, 2, 4, 8])

print(f"Strictly Stable: {ra_aux.is_stable}")
print(f"Number of RHP Poles: {ra_aux.num_rhp_poles}")
print(f"Auxiliary Polynomial: {ra_aux.auxiliary_polynomials}")

ra_aux

ra_aux = RouthArray([1, 2, 4, 8]) print(f"Strictly Stable: {ra_aux.is_stable}") print(f"Number of RHP Poles: {ra_aux.num_rhp_poles}") print(f"Auxiliary Polynomial: {ra_aux.auxiliary_polynomials}") ra_aux

5. Closed-Loop Parametric Gain Stability Bounds (K)¶

Open-loop plant: $$G(s) = \frac{1}{s(s+1)(s+2)} = \frac{1}{s^3 + 3s^2 + 2s}$$

Closed-loop characteristic equation with gain K: $$s^3 + 3s^2 + 2s + K = 0$$

In [ ]:

G = cp.tf([1], [1, 3, 2, 0])
ra_param = RouthArray(G, k_symbol="K")
print(f"Solved Stability Range for K: {ra_param.k_range}")
ra_param

G = cp.tf([1], [1, 3, 2, 0]) ra_param = RouthArray(G, k_symbol="K") print(f"Solved Stability Range for K: {ra_param.k_range}") ra_param

Time-Domain Verification of Stability Bounds¶

We verify the predicted critical gain K_cr = 6:

• K = 3 (Stable: damped step response)
• K = 6 (Marginally Stable: sustained oscillations)
• K = 10 (Unstable: diverging response)

In [ ]:

fig, ax = plt.subplots(figsize=(9, 5))
t_span = np.linspace(0, 15, 1000)

for K_val, color, label in [
    (3.0, "blue", "K = 3.0 (Stable)"),
    (6.0, "orange", "K = 6.0 (Critical / Marginally Stable)"),
    (10.0, "red", "K = 10.0 (Unstable)"),
]:
    L = K_val * G
    T_cl = cp.feedback(L, 1.0)
    resp = cp.step_response(T_cl, T=t_span)
    ax.plot(resp.t, resp.y, label=label, color=color, linewidth=1.8)

ax.set_title("Closed-Loop Step Response vs Controller Gain K", fontsize=13, fontweight="bold")
ax.set_xlabel("Time (s)", fontsize=11)
ax.set_ylabel("Output y(t)", fontsize=11)
ax.set_ylim(-3, 4)
ax.grid(True, linestyle="--", alpha=0.6)
ax.legend(fontsize=10, loc="upper left")
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(9, 5)) t_span = np.linspace(0, 15, 1000) for K_val, color, label in [ (3.0, "blue", "K = 3.0 (Stable)"), (6.0, "orange", "K = 6.0 (Critical / Marginally Stable)"), (10.0, "red", "K = 10.0 (Unstable)"), ]: L = K_val * G T_cl = cp.feedback(L, 1.0) resp = cp.step_response(T_cl, T=t_span) ax.plot(resp.t, resp.y, label=label, color=color, linewidth=1.8) ax.set_title("Closed-Loop Step Response vs Controller Gain K", fontsize=13, fontweight="bold") ax.set_xlabel("Time (s)", fontsize=11) ax.set_ylabel("Output y(t)", fontsize=11) ax.set_ylim(-3, 4) ax.grid(True, linestyle="--", alpha=0.6) ax.legend(fontsize=10, loc="upper left") plt.tight_layout() plt.show()

6. Step-by-Step Pedagogical Explanation (.explain_steps())¶

Instructors and students can call .explain_steps() to obtain a detailed mathematical derivation breakdown:

In [ ]:

for i, step in enumerate(ra_stable.explain_steps(), 1):
    print(f"{i}. {step}")

for i, step in enumerate(ra_stable.explain_steps(), 1): print(f"{i}. {step}")