05. Pedagogical & Symbolic Control Tools (ctrlpy.pedagogy)¶

Step-by-Step Routh-Hurwitz Analysis, Analytical Root Locus Rules, and Steady-State Error Derivations

This notebook explores the educational and symbolic submodule of ctrlpy (ctrlpy.pedagogy). While the core ctrlpy engine is built for 100% numerical speed with NumPy and SciPy, the ctrlpy.pedagogy module is designed for classroom instruction, homework verification, and pedagogical derivations using SymPy with rich LaTeX output.

In [ ]:

import matplotlib.pyplot as plt
import numpy as np
import sympy as sp

import ctrlpy as cp
from ctrlpy.pedagogy import (
    root_locus_rules,
    routh_table,
    steady_state_analysis,
)

%matplotlib inline

import matplotlib.pyplot as plt import numpy as np import sympy as sp import ctrlpy as cp from ctrlpy.pedagogy import ( root_locus_rules, routh_table, steady_state_analysis, ) %matplotlib inline

---

1. Routh-Hurwitz Stability Criterion¶

The Routh-Hurwitz stability criterion allows us to determine the number of closed-loop poles in the right-half of the complex s-plane (RHP) without explicitly factoring high-order polynomials.

Example 1: Standard 3rd-Order Stable System¶

Consider the characteristic polynomial: $$P(s) = s^3 + 2s^2 + 3s + 4 = 0$$

In [ ]:

# Construct Routh table from polynomial coefficients [s^3, s^2, s^1, s^0]
res1 = routh_table([1, 2, 3, 4])

# Display formatted LaTeX table in Jupyter
res1

# Construct Routh table from polynomial coefficients [s^3, s^2, s^1, s^0] res1 = routh_table([1, 2, 3, 4]) # Display formatted LaTeX table in Jupyter res1

In [ ]:

print("Detailed ASCII Output:")
print(res1)
print(f"\nStrictly Stable: {res1.is_stable}")
print(f"RHP Poles: {res1.num_rhp_poles}")

print("Detailed ASCII Output:") print(res1) print(f"\nStrictly Stable: {res1.is_stable}") print(f"RHP Poles: {res1.num_rhp_poles}")

Example 2: 4th-Order Unstable System¶

Consider the polynomial: $$P(s) = s^4 + 2s^3 + 3s^2 + 4s + 5 = 0$$

In [ ]:

s = sp.Symbol("s")
poly_expr = s**4 + 2 * s**3 + 3 * s**2 + 4 * s + 5

res2 = routh_table(poly_expr)
res2

s = sp.Symbol("s") poly_expr = s**4 + 2 * s**3 + 3 * s**2 + 4 * s + 5 res2 = routh_table(poly_expr) res2

---

2. Handling Routh-Hurwitz Special Cases¶

Textbook Routh-Hurwitz analysis features two classical edge cases:

1. Special Case 1 (Zero in First Column): A leading zero causes division by zero. It is replaced with symbolic $\epsilon > 0$, and column 1 signs are evaluated in the limit $\epsilon \to 0^+$.
2. Special Case 2 (Row of All Zeros): Occurs when roots are symmetrically distributed about the origin (e.g., pairs on the $j\omega$ axis). An auxiliary polynomial $A(s)$ is formed from the previous row and differentiated $\frac{dA}{ds}$ to complete the table.

Special Case 1: Epsilon Substitution ($\epsilon > 0$)¶

Consider the polynomial: $$P(s) = s^5 + 2s^4 + 2s^3 + 4s^2 + 11s + 10 = 0$$

Notice that row $s^3$ has a zero in the first position $((2\cdot 2 - 1\cdot 4)/2 = 0)$ while the next entry is non-zero $((2\cdot 11 - 1\cdot 10)/2 = 6)$.

In [ ]:

res_eps = routh_table([1, 2, 2, 4, 11, 10])
res_eps

res_eps = routh_table([1, 2, 2, 4, 11, 10]) res_eps

In [ ]:

print("Step-by-step notes:")
for note in res_eps.steps:
    print(f" - {note}")

print("Step-by-step notes:") for note in res_eps.steps: print(f" - {note}")

Special Case 2: Row of All Zeros & Auxiliary Polynomial¶

Consider the polynomial with marginal stability (imaginary axis poles at $s = \pm 2j$): $$P(s) = s^3 + 2s^2 + 4s + 8 = (s + 2)(s^2 + 4) = 0$$

Row $s^1$ is computed as $(2\cdot 4 - 1\cdot 8)/2 = 0$, resulting in an entire row of zeros.

In [ ]:

res_zeros = routh_table([1, 2, 4, 8])
res_zeros

res_zeros = routh_table([1, 2, 4, 8]) res_zeros

In [ ]:

print("Auxiliary Polynomials:", res_zeros.auxiliary_polynomials)
print("Construction Notes:")
for s_note in res_zeros.steps:
    print(" *", s_note)

print("Auxiliary Polynomials:", res_zeros.auxiliary_polynomials) print("Construction Notes:") for s_note in res_zeros.steps: print(" *", s_note)

---

3. Parametric Gain Stability Bounds ($K$-Range Solver)¶

In control design, we often need to find the range of feedback gains $K > 0$ that ensure closed-loop stability.

Let an open-loop plant be: $$G(s) = \frac{1}{s(s+1)(s+2)} = \frac{1}{s^3 + 3s^2 + 2s}$$

The closed-loop characteristic equation under proportional gain $K$ is: $$1 + K G(s) = 0 \implies s^3 + 3s^2 + 2s + K = 0$$

In [ ]:

G_plant = cp.tf([1], [1, 3, 2, 0])

# Solve parametric stability range for gain K directly from the TransferFunction
res_k = routh_table(G_plant, k_symbol="K")
res_k

G_plant = cp.tf([1], [1, 3, 2, 0]) # Solve parametric stability range for gain K directly from the TransferFunction res_k = routh_table(G_plant, k_symbol="K") res_k

In [ ]:

print(f"Calculated Stability Range for K: {res_k.k_range}")

# Verify with numerical simulations in ctrlpy:
# 1. Subcritical gain K = 3 (Stable)
T_stable = cp.feedback(3.0 * G_plant, 1.0)
print("Poles for K=3:", T_stable.poles())

# 2. Critical gain K = 6 (Marginally stable, poles on jw axis)
T_crit = cp.feedback(6.0 * G_plant, 1.0)
print("Poles for K=6:", T_crit.poles())

# 3. Supercritical gain K = 10 (Unstable)
T_unstable = cp.feedback(10.0 * G_plant, 1.0)
print("Poles for K=10:", T_unstable.poles())

print(f"Calculated Stability Range for K: {res_k.k_range}") # Verify with numerical simulations in ctrlpy: # 1. Subcritical gain K = 3 (Stable) T_stable = cp.feedback(3.0 * G_plant, 1.0) print("Poles for K=3:", T_stable.poles()) # 2. Critical gain K = 6 (Marginally stable, poles on jw axis) T_crit = cp.feedback(6.0 * G_plant, 1.0) print("Poles for K=6:", T_crit.poles()) # 3. Supercritical gain K = 10 (Unstable) T_unstable = cp.feedback(10.0 * G_plant, 1.0) print("Poles for K=10:", T_unstable.poles())

---

4. Analytical Root Locus Derivations¶

The root_locus_rules function automates the analytical Evans rules taught in university control courses:

1. Number of poles ($n$), zeros ($m$), and branches to infinity ($n - m$).
2. Real-axis root locus segments.
3. Asymptote centroid $\sigma_a = \frac{\sum p_i - \sum z_i}{n - m}$ and angles $\theta_k = \frac{(2k+1)180^\circ}{n-m}$.
4. Breakaway and break-in points via $\frac{dK}{ds} = 0$.
5. Departure angles from complex poles and arrival angles to complex zeros.
6. Imaginary axis crossings ($s = \pm j\omega$) and critical gain $K_{\text{crit}}$.

In [ ]:

# Derive analytical rules for G(s) = 1 / (s(s+1)(s+2))
rl_res = root_locus_rules(G_plant)
rl_res

# Derive analytical rules for G(s) = 1 / (s(s+1)(s+2)) rl_res = root_locus_rules(G_plant) rl_res

In [ ]:

# Visual comparison with ctrlpy's numerical root locus plotter
fig, ax = cp.plot_root_locus(G_plant, gains=np.linspace(0, 10, 1000))
ax.set_title("Root Locus with Analytical Rules Verification")
plt.show()

# Visual comparison with ctrlpy's numerical root locus plotter fig, ax = cp.plot_root_locus(G_plant, gains=np.linspace(0, 10, 1000)) ax.set_title("Root Locus with Analytical Rules Verification") plt.show()

Complex Poles & Departure Angles Example¶

Consider a plant with complex conjugate open-loop poles: $$G(s) = \frac{s + 2}{s^2 + 2s + 2} = \frac{s + 2}{(s + 1 + j)(s + 1 - j)}$$

In [ ]:

G_complex = cp.tf([1, 2], [1, 2, 2])
rl_complex = root_locus_rules(G_complex)
rl_complex

G_complex = cp.tf([1, 2], [1, 2, 2]) rl_complex = root_locus_rules(G_complex) rl_complex

---

5. System Type & Steady-State Tracking Error Analysis¶

The steady-state error $e_{ss} = \lim_{t \to \infty} e(t)$ of a unity feedback system depends on the number of open-loop integrators (System Type $N$):

System Type	Position Constant $K_p$	Velocity Constant $K_v$	Acceleration Constant $K_a$	Step Input $e_{ss}$	Ramp Input $e_{ss}$	Parabola Input $e_{ss}$
Type 0	Finite	$0$	$0$	$\frac{1}{1 + K_p}$	$\infty$	$\infty$
Type 1	$\infty$	Finite	$0$	$0$	$\frac{1}{K_v}$	$\infty$
Type 2	$\infty$	$\infty$	Finite	$0$	$0$	$\frac{1}{K_a}$

Let us analyze Type 0, Type 1, and Type 2 systems.

In [ ]:

# Type 0 System: G(s) = 10 / (s^2 + 3s + 2)
G_type0 = cp.tf([10], [1, 3, 2])
res_ss0 = steady_state_analysis(G_type0)
res_ss0

# Type 0 System: G(s) = 10 / (s^2 + 3s + 2) G_type0 = cp.tf([10], [1, 3, 2]) res_ss0 = steady_state_analysis(G_type0) res_ss0

In [ ]:

# Type 1 System: G(s) = 10 / (s(s + 2)) = 10 / (s^2 + 2s)
G_type1 = cp.tf([10], [1, 2, 0])
res_ss1 = steady_state_analysis(G_type1)
res_ss1

# Type 1 System: G(s) = 10 / (s(s + 2)) = 10 / (s^2 + 2s) G_type1 = cp.tf([10], [1, 2, 0]) res_ss1 = steady_state_analysis(G_type1) res_ss1

In [ ]:

# Type 2 System: G(s) = 10 / (s^2(s + 2)) = 10 / (s^3 + 2s^2)
G_type2 = cp.tf([10], [1, 2, 0, 0])
res_ss2 = steady_state_analysis(G_type2)
res_ss2

# Type 2 System: G(s) = 10 / (s^2(s + 2)) = 10 / (s^3 + 2s^2) G_type2 = cp.tf([10], [1, 2, 0, 0]) res_ss2 = steady_state_analysis(G_type2) res_ss2

---

6. Summary & Key Takeaways¶

• Architectural Isolation: ctrlpy keeps its core numerical pipeline fast and zero-overhead while offering ctrlpy.pedagogy for deep symbolic, classroom-grade derivations.
• Robust Special Cases: routh_table robustly handles epsilon substitutions for zero entries and auxiliary polynomials for symmetric/imaginary-axis poles.
• Parametric Gain Solvers: Direct computation of closed-loop stability margins for symbolic gain $K$.
• Classroom Analytical Rules: root_locus_rules and steady_state_analysis provide complete step-by-step mathematical reasoning with formatted LaTeX display.