Control systems analysis and PID design with python-control and scipy.signal: transfer functions, Bode/Nyquist plots, root locus, stability margins, and step-response simulation for SISO and discrete-time systems.
Control systems analysis and PID design with python-control and scipy.signal: transfer functions, Bode/Nyquist plots, root locus, stability margins, and step-response simulation for SISO and discrete-time systems.
One-line summary: This Skill helps engineers model, analyse, and design feedback control systems using python-control and scipy.signal, covering everything from transfer-function definition through PID synthesis, frequency-domain plots, and discrete-time simulation.
When to Use This Skill
When you need to model a plant as a transfer function or state-space representation
When you need to assess stability via gain/phase margin, Nyquist criterion, or pole locations
When you need to design a PID or lead-lag compensator using Ziegler-Nichols or analytical methods
When you need to visualise frequency response (Bode plot, Nyquist diagram, root locus)
When you need to simulate closed-loop step response and compute time-domain specifications
When you need to discretise a continuous controller for embedded implementation (ZOH, Tustin)
Trigger keywords: transfer function, state space, Bode plot, Nyquist plot, root locus, PID tuning, Ziegler-Nichols, gain margin, phase margin, step response, lead compensator, lag compensator, ZOH discretization, stability analysis, python-control
Background & Key Concepts
Transfer Function Representation
A linear time-invariant (LTI) system is described in the Laplace domain by its transfer function:
Ziegler-Nichols ultimate-gain method: bring the system to sustained oscillation with a P-only controller to find the ultimate gain $K_u$ and ultimate period $T_u$, then use the ZN table to compute $K_p$, $K_i$, $K_d$.
Stability Margins
The gain margin (GM) and phase margin (PM) quantify robustness against model uncertainty:
where $\omega_{pc}$ is the phase crossover frequency and $\omega_{gc}$ is the gain crossover frequency. A rule of thumb: GM > 6 dB and PM > 30° for adequate robustness.
import control
import scipy
import numpy as np
import matplotlib
print(f"python-control : {control.__version__}")
print(f"scipy : {scipy.__version__}")
print(f"numpy : {np.__version__}")
print(f"matplotlib : {matplotlib.__version__}")
# Expected: python-control : 0.9.x or higher
Core Workflow
Step 1: Define the Plant Model
import control
import numpy as np
import matplotlib.pyplot as plt
# ── 1a. Transfer function: second-order plant with time-delay approximation ───# G(s) = K / (s(τs + 1)) — a DC motor / velocity loop prototype
K_plant = 2.0
tau = 0.5# time constant (s)# Numerator and denominator polynomial coefficients (highest power first)
num = [K_plant]
den = [tau, 1, 0] # τs² + s = s(τs + 1)
G = control.tf(num, den)
print("Plant transfer function:")
print(G)
# ── 1b. State-space representation of the same plant ──────────────────────────# Canonical controllable form for G(s) = 2 / (0.5s² + s)
A = np.array([[0, 1],
[0, -1/tau]])
B = np.array([[0],
[K_plant/tau]])
C = np.array([[1, 0]])
D = np.array([[0]])
sys_ss = control.ss(A, B, C, D)
print("\nState-space representation:")
print(sys_ss)
# Verify equivalence
G_from_ss = control.ss2tf(sys_ss)
print("\nTransfer function from SS (should match):")
print(G_from_ss)
Nyquist: if the $(-1, 0)$ point is not encircled, the closed loop is stable (for open-loop stable $G$)
Step 3: PID Controller Design — Ziegler-Nichols
defziegler_nichols_pid(Ku: float, Tu: float, variant: str = "classic") -> dict:
"""
Compute PID gains using the Ziegler-Nichols ultimate-gain method.
Parameters
----------
Ku : float
Ultimate proportional gain (sustains oscillation).
Tu : float
Ultimate period (seconds) at Ku.
variant : str
'classic' — original ZN table (aggressive, slight overshoot)
'some_overshoot' — modified for ~20 % overshoot
'no_overshoot' — modified for near-zero overshoot
Returns
-------
dict with keys Kp, Ki, Kd
"""if variant == "classic":
Kp = 0.6 * Ku
Ti = 0.5 * Tu
Td = 0.125 * Tu
elif variant == "some_overshoot":
Kp = 0.33 * Ku
Ti = 0.5 * Tu
Td = 0.33 * Tu
elif variant == "no_overshoot":
Kp = 0.2 * Ku
Ti = 0.5 * Tu
Td = 0.33 * Tu
else:
raise ValueError(f"Unknown variant: {variant!r}")
Ki = Kp / Ti
Kd = Kp * Td
return {"Kp": Kp, "Ki": Ki, "Kd": Kd}
# Simulate proportional-only closed loop to find Ku, Tu ─────────────────────# For this plant we can compute analytically: at phase = -180°# G(jω) phase = -90° - arctan(τω) = -180° → arctan(0.5ω) = 90° → no finite ω# The plant G = 2/[s(0.5s+1)] has inherent -90° from the integrator, so we add# a lag correction factor and use gain-based approach from margin analysis# Demonstration: design PID using the margins computed earlier# For a real process: run a relay experiment or apply proportional-only control# Here we approximate: Ku from gain margin, Tu from phase crossover period
Ku_approx = 2.5# would be found experimentally
Tu_approx = 2 * np.pi / wpc if wpc > 0else2.0
gains = ziegler_nichols_pid(Ku_approx, Tu_approx, variant="classic")
print(f"\nZiegler-Nichols PID gains (classic):")
print(f" Kp = {gains['Kp']:.4f}")
print(f" Ki = {gains['Ki']:.4f}")
print(f" Kd = {gains['Kd']:.4f}")
# Build the PID controller transfer function: C(s) = Kp + Ki/s + Kd*s# With derivative filter to prevent differentiator wind-up: Kd*s/(s/N + 1), N=10
N = 10
Kp, Ki, Kd = gains["Kp"], gains["Ki"], gains["Kd"]
C_num = [Kd + Kp/N, Kp + Ki/N, Ki]
C_den = [1, N, 0]
C = control.tf(C_num, C_den)
print("\nPID controller transfer function C(s):")
print(C)
Cause: Attempting to create a state-space system with incompatible D matrix dimensions when the plant has more outputs than inputs.
Fix:
# Ensure D has shape (n_outputs, n_inputs)
n_outputs, n_inputs = 1, 1
D = np.zeros((n_outputs, n_inputs))
sys = control.ss(A, B, C, D)
Error: control.matlab.sisotool not found
Cause: sisotool is only available in the MATLAB Control System Toolbox; python-control provides control.root_locus and interactive design via control.sisotool in newer releases.
Fix:
pip install control>=0.9.4 # sisotool added in 0.9# Then in Python:# control.sisotool(G) # opens interactive root-locus GUI
Issue: Nyquist plot looks incomplete or truncated
Cause: The default omega range may not span the full frequency range of interest.
Fix:
import numpy as np
omega = np.logspace(-3, 3, 2000) # 0.001 to 1000 rad/s
control.nyquist_plot(G, omega=omega)
Version Compatibility
Package
Tested versions
Known issues
control
0.9.4, 0.10.x
API change in 0.9: bode_plot returns (mag, phase, omega) not (mag, phase)
Interpreting these results: A phase margin above 45° with the PID controller confirms robust stability. The rise time target of < 50 ms at 50 rad/s bandwidth is met. Adjust Kd_pid to reduce overshoot and Ki_pid to improve steady-state tracking of ramp references.
Example 2: Lead-Lag Compensator for Unstable Plant
Scenario: Design a lead-lag compensator for a marginally-stable inverted pendulum linearisation and verify performance.
# =============================================# End-to-end example 2: lead-lag compensator# =============================================import control
import numpy as np
import matplotlib.pyplot as plt
# Linearised inverted pendulum on cart: G(s) = 1 / (s² - ω_n²)
omega_n = 3.0# rad/s (unstable pole at ±3 rad/s)
G_inv = control.tf([1], [1, 0, -omega_n**2])
print("Inverted pendulum plant:")
print(G_inv)
print(f"Open-loop poles: {control.poles(G_inv)}")
# Phase lead to stabilise and add PM ≈ 45°# Required PM increase: 45° - (current PM) ≈ large value for unstable plant
alpha_l = 0.1
T_l = 0.5
C_lead_ip = control.tf([T_l, 1], [alpha_l * T_l, 1])
# Gain to place gain crossover at 6 rad/s# |C_lead * G| at ω=6 should be 1 → compute and scale
omega_target = 6.0
freq_resp = control.freqresp(C_lead_ip * G_inv, [omega_target])
mag_at_target = abs(freq_resp.fresp[0, 0, 0])
K_scale = 1.0 / mag_at_target
C_total = K_scale * C_lead_ip
L_ip = C_total * G_inv
gm_ip, pm_ip, _, wgc_ip = control.stability_margins(L_ip)
print(f"\nWith lead: PM = {pm_ip:.1f}°, bandwidth ≈ {wgc_ip:.2f} rad/s")
T_cl_ip = control.feedback(L_ip, 1)
print(f"Closed-loop poles: {control.poles(T_cl_ip)}")
t_ip = np.linspace(0, 5, 1000)
t_ip_out, y_ip_out = control.step_response(T_cl_ip, T=t_ip)
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
control.bode_plot(L_ip, dB=True, margins=True, omega_limits=(0.1, 100),
ax=axes)
axes[0].set_title("Lead Compensated Loop")
fig2, ax2 = plt.subplots(figsize=(8, 4))
ax2.plot(t_ip_out, y_ip_out, linewidth=2)
ax2.axhline(1, color="r", linestyle="--", alpha=0.6, label="Reference")
ax2.set_xlabel("Time (s)")
ax2.set_ylabel("Output")
ax2.set_title("Inverted Pendulum — Lead Compensated Step Response")
ax2.legend()
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("inverted_pendulum_step.png", dpi=150, bbox_inches="tight")
plt.show()
print("Analysis complete.")
print(f" Phase margin : {pm_ip:.1f}°")
print(f" All closed-loop poles in LHP: "f"{all(p.real < 0for p in control.poles(T_cl_ip))}")
Interpreting these results: The key check is that all closed-loop poles have negative real parts (LHP), confirming stabilisation. A phase margin above 30° ensures robustness to model uncertainty.
Last updated: 2026-03-17 | Maintainer: @xjtulycIssues: GitHub Issues
