Control Systems Cheatsheet Cheatsheet

🎯

Control System Basics

CORE

Open-Loop vs Closed-Loop

Aspect	Open-Loop	Closed-Loop
Feedback	None	Output measured and fed back
Accuracy	Low — affected by disturbances	High — compensates disturbances
Stability	Always stable (no feedback)	Must be designed for stability
Complexity	Simple, cheap	Complex, more expensive
Example	Toaster, washing machine timer	Cruise control, thermostat
Sensitivity	High to parameter changes	Low — reduced by loop gain

Feedback Advantages

S_G = 1/(1 + GH) — sensitivity to forward path G

Y_d = G₂D / (1 + G₁G₂H) — feedback reduces disturbance effect

Steady-state error ∝ 1/(1 + K_p) for type 0 system

BW_closed > BW_open (faster response)

Higher loop gain → better accuracy but reduced stability margin

open_vs_closed_loop.py

import numpy as np

# Open-loop vs Closed-loop response comparison
# Plant: G(s) = K / (s + a)
K = 10;  a = 2
dt = 0.01;  t = np.arange(0, 5, dt)

# Open-loop: u = r = 1 (no feedback)
def simulate_open(K, a, r=1.0, dt=0.01, n=500):
    y = np.zeros(n)
    for i in range(1, n):
        y[i] = y[i-1] + dt * (-a * y[i-1] + K * r)
    return y

# Closed-loop: u = Kc(r - y), plant same
def simulate_closed(K, a, Kc, r=1.0, dt=0.01, n=500):
    y = np.zeros(n)
    for i in range(1, n):
        u = Kc * (r - y[i-1])
        y[i] = y[i-1] + dt * (-a * y[i-1] + K * u)
    return y

y_open = simulate_open(K, a)
y_closed = simulate_closed(K, a, Kc=2)

# Compare steady-state
ss_open = y_open[-1]
ss_closed = y_closed[-1]
print(f"Open-loop steady state:  {ss_open:.4f}  (error = {1-ss_open:.4f})")
print(f"Closed-loop SS (Kc=2):   {ss_closed:.4f}  (error = {1-ss_closed:.4f})")
print(f"Closed-loop SS (Kc=5):   {simulate_closed(K,a,5)[-1]:.4f}  (error = {1-simulate_closed(K,a,5)[-1]:.4f})")

💡

The characteristic equation 1 + G(s)H(s) = 0 determines closed-loop stability. The roots of this equation are the closed-loop poles — they must all lie in the left-half s-plane for stability.

📐

Transfer Functions

CORE

Block Diagram Algebra

Configuration	Equivalent TF	Note
Series: G₁ → G₂	G₁G₂	Multiply transfer functions
Parallel: G₁ + G₂	G₁ + G₂	Sum transfer functions
Negative feedback	G/(1+GH)	Most common configuration
Positive feedback	G/(1−GH)	Unstable if GH > 1
Minor feedback loop	G₁G₂/(1+G₂H)	Reduce inner loop first

Standard Transfer Functions

G(s) = K / (τs + 1); τ = time constant

G(s) = ω_n² / (s² + 2ζω_n s + ω_n²)

Steady-state output / steady-state input

Determines overshoot and oscillation

Frequency of oscillation when ζ = 0

s = −ζω_n ± ω_n√(ζ²−1)

block_diagram_reduction.py

# Block diagram algebra — Mason's gain formula
# T(s) = (Σ Δ_k P_k) / Δ
# Δ = 1 - (Σ L_i) + (Σ L_i L_j) - (Σ L_i L_j L_k) + ...
# Δ_k = Δ with loop k removed

# Example: Two overlapping feedback loops
#         ┌──H₁──┐
#   R ─→G₁─→G₂─→G₃─→ Y
#         └──H₂──┘

G1 = "G1"; G2 = "G2"; G3 = "G3"
H1 = "H1"; H2 = "H2"

# Forward paths
P1 = "G1 * G2 * G3"  # only one forward path

# Individual loops
L1 = "-G2 * H1"      # inner loop
L2 = "-G1 * G2 * G3 * H2"  # outer loop

# Δ = 1 - L1 - L2 + L1*L2 (no non-touching loops)
# Δ1 = 1 (no non-touching loops to path 1)

print("Mason's Gain Formula:")
print(f"Forward path: P1 = {P1}")
print(f"Loop 1: L1 = {L1}")
print(f"Loop 2: L2 = {L2}")
print(f"Δ = 1 - L1 - L2 = 1 + G2*H1 + G1*G2*G3*H2")
print(f"T(s) = G1*G2*G3 / (1 + G2*H1 + G1*G2*G3*H2)")
print()
print("With numerical values: G1=5, G2=10, G3=2, H1=0.1, H2=0.05")
G1n, G2n, G3n = 5, 10, 2
H1n, H2n = 0.1, 0.05
T = (G1n * G2n * G3n) / (1 + G2n*H1n + G1n*G2n*G3n*H2n)
print(f"T = {T:.4f}")

Poles & Zeros

Location	System Type	Response Characteristic
LHP real	Stable	Exponential decay
RHP real	Unstable	Exponential growth
LHP complex conjugate	Stable	Damped oscillation
RHP complex conjugate	Unstable	Growing oscillation
jω axis (simple)	Marginally stable	Sustained oscillation
jω axis (repeated)	Unstable	Growing oscillation
Origin (s = 0)	Type N system	N integrators; N pole at origin

⚠️

Zeros affect the transient response shape (overshoot direction, settling time) but do NOT affect stability. Poles alone determine absolute stability.

📈

Time Domain Analysis

CORE

First Order System

G(s) = K / (τs + 1)

y(t) = K(1 − e^(−t/τ))

Time to reach 63.2% of final value

t_s = 4τ

t_s = 3τ

t_r = 2.2τ

First-order systems never overshoot

Second Order System

G(s) = ω_n² / (s² + 2ζω_n s + ω_n²)

ζ > 1: overdamped | ζ = 1: critical | 0 < ζ < 1: underdamped

t_p = π / (ω_n √(1−ζ²))

Mp = e^(−πζ/√(1−ζ²)) × 100%

t_s = 4 / (ζω_n)

t_r ≈ (π − φ) / (ω_d); ω_d = ω_n√(1−ζ²)

second_order_response.py

import numpy as np
import math

def second_order_specs(zeta, wn):
    """Calculate second-order system specifications."""
    if zeta < 1:
        wd = wn * math.sqrt(1 - zeta**2)
        tp = math.pi / wd
        Mp = math.exp(-math.pi * zeta / math.sqrt(1 - zeta**2)) * 100
        tr = (math.pi - math.atan2(math.sqrt(1-zeta**2), zeta)) / wd
    else:
        tp = float('inf')
        Mp = 0.0
        tr = float('inf')
    ts = 4.0 / (zeta * wn)  # 2% criterion
    return tp, Mp, tr, ts

print(f"{'ζ':>5s}  {'ω_n':>5s}  {'t_p (s)':>8s}  {'Mp (%)':>7s}  {'t_r (s)':>8s}  {'t_s (2%)(s)':>11s}")
print("-" * 60)

wn = 10  # rad/s
for zeta in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.707, 0.8, 1.0, 1.5]:
    tp, Mp, tr, ts = second_order_specs(zeta, wn)
    if tp == float('inf'):
        tp_str = "  N/A   "
    else:
        tp_str = f"{tp:>8.4f}"
    if tr == float('inf'):
        tr_str = "  N/A   "
    else:
        tr_str = f"{tr:>8.4f}"
    print(f"{zeta:>5.3f}  {wn:>5.1f}  {tp_str}  {Mp:>7.2f}  {tr_str}  {ts:>11.4f}")

print(f"\nNote: ζ = 0.707 → Mp = {second_order_specs(0.707, 10)[1]:.2f}% (critically damped, optimal)")

Steady-State Error (Step Input)

System Type	Step (r=1/s)	Ramp (r=1/s²)	Parabola (r=1/s³)
Type 0 (no pole at origin)	1/(1+Kp)	∞	∞
Type 1 (1 pole at origin)	0	1/Kv	∞
Type 2 (2 poles at origin)	0	0	1/Ka

Kp = lim(s→0) G(s)H(s)

Kv = lim(s→0) sG(s)H(s)

Ka = lim(s→0) s²G(s)H(s)

💡

For a well-designed system, aim for ζ = 0.4–0.8. ζ = 0.707 gives Mp ≈ 4.3% (often called "critically damped" in practice) and offers the best trade-off between speed and overshoot.

⚖️

Stability Analysis

CRITICAL

Routh-Hurwitz Criterion

All coefficients of characteristic equation must be positive

All elements in first column of Routh array must be positive

Number of sign changes = number of RHP poles

Replace zero with ε → 0, continue table

Form auxiliary polynomial, differentiate, replace row

routh_hurwitz.py

# Routh-Hurwitz stability criterion
def routh_hurwitz(coeffs):
    """
    coeffs: list of polynomial coefficients [a_n, a_{n-1}, ..., a_0]
    Returns: (stable: bool, routh_table, rhp_poles: int)
    """
    n = len(coeffs) - 1
    rows = []

    # First two rows from coefficients
    rows.append(coeffs[0::2])     # even index
    rows.append(coeffs[1::2])     # odd index
    # Pad to same length
    max_len = max(len(rows[0]), len(rows[1]))
    rows[0] += [0] * (max_len - len(rows[0]))
    rows[1] += [0] * (max_len - len(rows[1]))

    while len(rows) < n + 1:
        prev2 = rows[-2]
        prev1 = rows[-1]
        row = []
        for j in range(len(prev1) - 1):
            elem = (prev2[0] * prev1[j+1] - prev1[0] * prev2[j+1]) / prev1[0]
            row.append(round(elem, 6))
        if all(abs(x) < 1e-12 for x in row):
            row = [0] * len(prev1)
        rows.append(row)
        if not row:
            break

    # Count sign changes in first column
    first_col = [r[0] if r else 0 for r in rows]
    sign_changes = 0
    for i in range(1, len(first_col)):
        if first_col[i] != 0 and first_col[i-1] != 0:
            if first_col[i] * first_col[i-1] < 0:
                sign_changes += 1

    stable = sign_changes == 0
    return stable, rows, sign_changes

# Example: s⁴ + 2s³ + 3s² + 4s + 5 = 0
coeffs = [1, 2, 3, 4, 5]
stable, table, rhp = routh_hurwitz(coeffs)
print(f"Characteristic equation: s⁴ + 2s³ + 3s² + 4s + 5 = 0")
print(f"Stable: {stable}")
print(f"RHP poles: {rhp}")
print("\nRouth array:")
for row in table:
    print([f"{x:>8.4f}" if abs(x) > 1e-12 else "       0" for x in row])

# Stable example: s³ + 6s² + 11s + 6 = 0 (poles: -1, -2, -3)
stable2, _, rhp2 = routh_hurwitz([1, 6, 11, 6])
print(f"\ns³ + 6s² + 11s + 6: Stable={stable2}, RHP={rhp2}")

Stability Margins

GM = 1/|G(jω_gc)| in dB at phase crossover freq ω_gc

PM = 180° + ∠G(jω_pc) at gain crossover freq ω_pc

≥ 6 dB (factor of 2)

≥ 30° (typically 45°–60° for good design)

Frequency where |G(jω)| = 1 (0 dB)

Frequency where ∠G(jω) = −180°

🚫

A system can be stable but have poor relative stability (margins too small). Always check BOTH gain margin and phase margin. In practice, GM ≥ 6 dB and PM ≥ 45° provides robust performance with ~15–25% overshoot.

🌳

Root Locus

ANALYSIS

Root Locus Construction Rules

Rule	Statement
#1 Branches	Number of branches = max(n, m); n = poles, m = zeros
#2 Symmetry	Root locus is symmetric about real axis
#3 Real axis segments	Locus exists to the LEFT of odd number of poles+zeros
#4 Asymptotes	Angles: (2k+1)π/(n−m), k=0,1,... Centroid: (Σp−Σz)/(n−m)
#5 Breakaway/Break-in	On real axis; solve dK/ds = 0
#6 Departure angle	From complex pole: θ_dep = 180° + Σ∠(to other poles) − Σ∠(to zeros)
#7 Arrival angle	To complex zero: θ_arr = 180° − Σ∠(to poles) + Σ∠(from other zeros)
#8 Imaginary crossing	Routh array: auxiliary polynomial → jω crossing frequency
#9 K at any point	K = \|product(s−pᵢ)\| / \|product(s−zⱼ)\|

root_locus_sketch.py

import numpy as np
import matplotlib.pyplot as plt

# Root locus for G(s) = K / (s(s+2)(s+4))
# Characteristic equation: s³ + 6s² + 8s + K = 0
# 1 + K/s(s+2)(s+4) = 0

poles = np.array([0, -2, -4], dtype=float)
zeros = np.array([], dtype=float)

n = len(poles);  m = len(zeros)

# Asymptotes
n_m = n - m
angles = [(2*k + 1) * 180 / n_m for k in range(n_m)]
centroid = (poles.sum() - zeros.sum()) / n_m

print(f"Poles: {poles}")
print(f"Zeros: {zeros}")
print(f"Asymptote angles: {angles}°")
print(f"Asymptote centroid: {centroid:.2f}")

# Breakaway point on real axis (between s=0 and s=-2)
# 1/(s-0) + 1/(s+2) + 1/(s+4) = 0  →  3s² + 12s + 8 = 0
a, b, c = 3, 12, 8
disc = b**2 - 4*a*c
s1 = (-b + np.sqrt(disc)) / (2*a)  # not between 0 and -2
s2 = (-b - np.sqrt(disc)) / (2*a)  # between 0 and -2
print(f"Breakaway points: {s1:.3f}, {s2:.3f}")
print(f"Valid breakaway: s = {s2:.3f} (between 0 and -2)")

# K at breakaway
K_break = abs(s2) * abs(s2+2) * abs(s2+4)
print(f"K at breakaway = {K_break:.3f}")

# Imaginary axis crossing: s = jω in s³+6s²+8s+K=0
# Routh: 6ω²−K = 0 and −ω²(6) + 8 = 0 ... or auxiliary poly from row s¹
# From Routh: row s¹ = [(48-K)/6, 0], crossing when K = 48
# Row s² = 6s² + K = 0 → s² = -K/6 → ω = √(48/6) = 2√2
omega_cross = np.sqrt(48/6)
print(f"Imaginary crossing: jω = ±j{omega_cross:.3f}, K = 48")
print(f"System unstable for K > 48")

Root Locus Interpretation

Closest to jω axis; dictate response shape

Add zeros (lead compensator) — improves stability

Add poles (lag compensator) — reduces stability

Draw line at angle θ = cos⁻¹(ζ) from origin; intersection gives K

Place compensator to reshape locus through desired poles

⚠️

For a desired closed-loop response with ζ = 0.5 and ω_n = 5 rad/s, the desired poles are at −2.5 ± j4.33. On the root locus, find K that places poles at (or near) this location. If not possible, add compensator poles/zeros.

📊

Bode Plot

ANALYSIS

Standard Bode Plot Elements

Element	Magnitude (dB)	Phase (degrees)
Gain K	20log₁₀(K)	0°
Integrator 1/s	−20 dB/decade slope	−90°
Differentiator s	+20 dB/decade slope	+90°
Real pole 1/(τs+1)	−20 dB/dec after ω=1/τ	0° → −90°
Real zero (τs+1)	+20 dB/dec after ω=1/τ	0° → +90°
Complex poles	−40 dB/dec after ω_n	0° → −180° (depends on ζ)
Time delay e^(−sT)	0 dB	−ωT × (180/π)°

Bode Plot Construction Steps

Express G(s) in standard form: K(τ₁s+1)... / (s)(τₐs+1)...

Identify all break frequencies: ω = 1/τ for each pole/zero

Start with K (and slope from integrators/differentiators)

At each corner freq, add ±20 or ±40 dB/decade

Add phase contributions: ±45° at corner, ±90° one decade away

bode_plot.py

import numpy as np

def bode_magnitude(omega, K, poles_tau, zeros_tau, n_integrators=0):
    """Compute Bode magnitude in dB."""
    mag = K * omega**(-n_integrators)
    for tau in poles_tau:
        mag *= 1.0 / np.sqrt(1 + (omega * tau)**2)
    for tau in zeros_tau:
        mag *= np.sqrt(1 + (omega * tau)**2)
    return 20 * np.log10(mag)

def bode_phase(omega, poles_tau, zeros_tau, n_integrators=0):
    """Compute Bode phase in degrees."""
    phase = -n_integrators * 90
    for tau in poles_tau:
        phase -= np.arctan(omega * tau) * 180 / np.pi
    for tau in zeros_tau:
        phase += np.arctan(omega * tau) * 180 / np.pi
    return phase

# G(s) = 10 / [s(0.5s+1)(0.1s+1)]
K = 10
pole_taus = [0.5, 0.1]
zero_taus = []
n_int = 1

omega = np.logspace(-2, 3, 1000)
mag_db = bode_magnitude(omega, K, pole_taus, zero_taus, n_int)
phase_deg = bode_phase(omega, pole_taus, zero_taus, n_int)

# Find gain crossover (|G| = 0 dB)
idx_gc = np.argmin(np.abs(mag_db))
omega_gc = omega[idx_gc]
pm = phase_deg[idx_gc] + 180

# Find phase crossover (phase = -180°)
idx_pc = np.argmin(np.abs(phase_deg + 180))
omega_pc = omega[idx_pc]
gm = -mag_db[idx_pc]

print(f"Gain crossover: ω_gc = {omega_gc:.3f} rad/s")
print(f"Phase margin:   PM = {pm:.1f}°")
print(f"Phase crossover: ω_pc = {omega_pc:.3f} rad/s")
print(f"Gain margin:    GM = {gm:.1f} dB")
print(f"\nCorner frequencies: {', '.join(f'{1/tau:.1f} rad/s' for tau in pole_taus)}")

Nyquist Criterion

N = P − Z; N = encirclements of (−1, 0) by GH(jω) plot

N + Z = P → Z = P − N = 0 for stability

Open-loop RHP poles

Closed-loop RHP poles

Read directly from Nyquist distance to (−1, 0)

💡

Bode plots are faster to sketch than Nyquist for systems with no open-loop RHP poles. However, Nyquist criterion handles non-minimum phase systems (RHP zeros) where Bode plots can be misleading.

🔢

State Space Representation

ADVANCED

State Space Model

ẋ = Ax + Bu

y = Cx + Du

x ∈ ℝⁿ; u ∈ ℝᵐ (inputs); y ∈ ℝᵖ (outputs)

n = number of state variables = number of energy storage elements

Controllability & Observability

U = [B AB A²B ... Aⁿ⁻¹B]; rank(U) = n → controllable

V = [C CA CA² ... CAⁿ⁻¹]; rank(V) = n → observable

All states can be reached from input in finite time

All states can be determined from output in finite time

Uncontrollable modes are stable

Unobservable modes are stable

state_space_analysis.py

import numpy as np

# State space: Mass-spring-damper system
# m ẍ + c ẋ + k x = F
# States: x₁ = x, x₂ = ẋ

m = 2.0;  c = 1.5;  k = 10.0

A = np.array([
    [0,  1],
    [-k/m, -c/m]
])
B = np.array([[0], [1/m]])
C = np.array([[1, 0]])
D = np.array([[0]])

# Eigenvalues → poles
eigenvalues = np.linalg.eigvals(A)
print(f"Eigenvalues (poles): {eigenvalues}")
print(f"Natural freq: {np.sqrt(k/m):.3f} rad/s")
print(f"Damping ratio: {c/(2*np.sqrt(k*m)):.3f}")

# Controllability check
n = A.shape[0]
U = np.hstack([B] + [A @ B for _ in range(n - 1)])
rank_U = np.linalg.matrix_rank(U)
print(f"\nControllability rank: {rank_U} (need {n}) → {'Controllable' if rank_U == n else 'NOT controllable'}")

# Observability check
V = np.vstack([C] + [C @ np.linalg.matrix_power(A, i) for i in range(1, n)])
rank_V = np.linalg.matrix_rank(V)
print(f"Observability rank:  {rank_V} (need {n}) → {'Observable' if rank_V == n else 'NOT observable'}")

# Transfer function from state space
# G(s) = C(sI - A)^(-1)B + D
from scipy import signal
sys = signal.StateSpace(A, B, C, D)
num, den = signal.tf2ss(*signal.ss2tf(A, B, C, D))
tf_sys = signal.TransferFunction(A, B, C, D)
print(f"\nTransfer function:")
print(f"  Numerator:   {tf_sys.num}")
print(f"  Denominator: {tf_sys.den}")

State Space → Transfer Function

G(s) = C(sI − A)⁻¹B + D

det(sI − A) = 0 → same as denominator of G(s)

A = companion matrix from denominator

A = transpose of controllable canonical form

ẑ = T⁻¹x → Ã = T⁻¹AT; B̄ = T⁻¹B; C̄ = CT

State Feedback & Observers

u = −Kx → closed-loop A_cl = A − BK; place poles via K

K = [0 ... 0 1] U⁻¹ φ(A); U = controllability matrix

ẋ̂ = Aẋ̂ + Bu + L(y − Cẋ̂); L places observer poles

Controller and observer can be designed independently

Place 2–5× faster than controller poles for good estimation

⚠️

A system is both controllable AND observable if and only if there are no pole-zero cancellations in its transfer function. Any cancellation hides a mode that is either uncontrollable or unobservable.

🎛️

PID Controllers

PRACTICAL

PID Controller Forms

u(t) = Kp·e(t) + Ki·∫e dt + Kd·de/dt

G_c(s) = Kp + Ki/s + Kd·s = Kp(1 + 1/(Tᵢs) + Td s)

PB = 100/Kp (%) — older nomenclature

Tᵢ = Kp/Ki

Td = Kd/Kp

G_c(s) = Kp(1 + 1/Tᵢs + Td·s/(1 + Td·s/N)); N = 8–20

Effect of Each Term

Term	Response Speed	Overshoot	Steady-State Error	Stability
P only	Faster	Increases with Kp	Nonzero (Type 0)	Degrades
PI	Moderate	Increases	Zero (Type 1)	Degrades further
PD	Much faster	Reduces	Nonzero	Improves
PID	Fast	Controlled	Zero	Good (tuned)

pid_tuning.py

import numpy as np

def ziegler_nichols_pid(ku, tu, mode='pid'):
    """Ziegler-Nichols tuning rules.
    ku = ultimate gain (at stability limit)
    tu = ultimate period (oscillation period at ku)
    """
    rules = {
        'p':   {'Kp': 0.5 * ku,    'Ki': 0,           'Kd': 0},
        'pi':  {'Kp': 0.45 * ku,   'Ki': 0.54*ku/tu,  'Kd': 0},
        'pid': {'Kp': 0.6 * ku,    'Ki': 1.2*ku/tu,   'Kd': 0.075*ku*tu},
    }
    return rules[mode]

# Example: plant with Ku = 4, Tu = 2 seconds
Ku = 4.0;  Tu = 2.0

for mode in ['p', 'pi', 'pid']:
    params = ziegler_nichols_pid(Ku, Tu, mode)
    Ti = params['Kp'] / params['Ki'] if params['Ki'] != 0 else float('inf')
    Td = params['Kd'] / params['Kp'] if params['Kp'] != 0 else 0
    print(f"{mode.upper():>4s}: Kp={params['Kp']:.3f}  Ki={params['Ki']:.3f}  Kd={params['Kd']:.3f}"
          f"  (Ti={Ti:.3f}  Td={Td:.3f})")

# Cohen-Coon tuning (for first-order + dead-time: K·e^(-sL)/(τs+1))
# Parameters: K=1, tau=2, L=0.5
K_plant = 1.0;  tau = 2.0;  L = 0.5
ratio = L / tau

# Cohen-Coon PID
Kp_cc = (tau / (K_plant * L)) * (4/3 + ratio/4)
Ti_cc = L * (32 + 6*ratio) / (13 + 8*ratio)
Td_cc = L * (4 / (11 + 2*ratio))
Ki_cc = Kp_cc / Ti_cc
Kd_cc = Kp_cc * Td_cc

print(f"\nCohen-Coon PID:")
print(f"  Kp={Kp_cc:.3f}  Ki={Ki_cc:.3f}  Kd={Kd_cc:.3f}")
print(f"  Ti={Ti_cc:.3f}  Td={Td_cc:.3f}")
print(f"  (Good for processes with significant dead time)")

Compensator Design

Type	Transfer Function	Effect	Application
Lead	Gc = K(1+sτ₁)/(1+sτ₂); τ₁>τ₂	Adds PM, speeds up	Improve transient response
Lag	Gc = K(1+sτ₁)/(1+sτ₂); τ₁<τ₂	Increases low-freq gain	Reduce steady-state error
Lead-Lag	Combination	Both improvements	Comprehensive design
PD	Gc = Kp + Kd·s ≈ lead	Adds derivative action	Reduce overshoot
PI	Gc = Kp + Ki/s ≈ lag	Adds integral action	Eliminate steady-state error

PID Implementation Notes

Limit integrator when actuator saturates; clamping or back-calc

Filter on error; use derivative on measurement, not error

Use different gains on setpoint vs measurement: b·Kp·(r−y)

Ts ≤ 0.1·τ_min; typically 10–100 ms for process control

Initialize integrator to match current output when switching modes

🚫

The Ziegler-Nichols method often produces aggressive tuning (high overshoot ~25%). For less overshoot, use modified ZN: Kp = 0.33Ku, Ti = Tu/2, Td = Tu/3 for PID. Modern auto-tune methods (relay method, IMC) generally outperform classical ZN.

⏳

Loading cheatsheet...

Aspect

Open-Loop

Closed-Loop

Feedback

None

Output measured and fed back

Accuracy

Low — affected by disturbances

High — compensates disturbances

Stability

Always stable (no feedback)

Must be designed for stability

Complexity

Simple, cheap

Complex, more expensive

Example

Toaster, washing machine timer

Cruise control, thermostat

Sensitivity

High to parameter changes

Low — reduced by loop gain

import numpy as np # Open-loop vs Closed-loop response comparison # Plant: G(s) = K / (s + a) K = 10; a = 2 dt = 0.01; t = np.arange(0, 5, dt) # Open-loop: u = r = 1 (no feedback) def simulate_open(K, a, r=1.0, dt=0.01, n=500): y = np.zeros(n) for i in range(1, n): y[i] = y[i-1] + dt * (-a * y[i-1] + K * r) return y # Closed-loop: u = Kc(r - y), plant same def simulate_closed(K, a, Kc, r=1.0, dt=0.01, n=500): y = np.zeros(n) for i in range(1, n): u = Kc * (r - y[i-1]) y[i] = y[i-1] + dt * (-a * y[i-1] + K * u) return y y_open = simulate_open(K, a) y_closed = simulate_closed(K, a, Kc=2) # Compare steady-state ss_open = y_open[-1] ss_closed = y_closed[-1] print(f"Open-loop steady state: {ss_open:.4f} (error = {1-ss_open:.4f})") print(f"Closed-loop SS (Kc=2): {ss_closed:.4f} (error = {1-ss_closed:.4f})") print(f"Closed-loop SS (Kc=5): {simulate_closed(K,a,5)[-1]:.4f} (error = {1-simulate_closed(K,a,5)[-1]:.4f})")

Configuration

Equivalent TF

Note

Series: G₁ → G₂

G₁G₂

Multiply transfer functions

Parallel: G₁ + G₂

G₁ + G₂

Sum transfer functions

Negative feedback

G/(1+GH)

Most common configuration

Positive feedback

G/(1−GH)

Unstable if GH > 1

Minor feedback loop

G₁G₂/(1+G₂H)

Reduce inner loop first

# Block diagram algebra — Mason's gain formula # T(s) = (Σ Δ_k P_k) / Δ # Δ = 1 - (Σ L_i) + (Σ L_i L_j) - (Σ L_i L_j L_k) + ... # Δ_k = Δ with loop k removed # Example: Two overlapping feedback loops # ┌──H₁──┐ # R ─→G₁─→G₂─→G₃─→ Y # └──H₂──┘ G1 = "G1"; G2 = "G2"; G3 = "G3" H1 = "H1"; H2 = "H2" # Forward paths P1 = "G1 * G2 * G3" # only one forward path # Individual loops L1 = "-G2 * H1" # inner loop L2 = "-G1 * G2 * G3 * H2" # outer loop # Δ = 1 - L1 - L2 + L1*L2 (no non-touching loops) # Δ1 = 1 (no non-touching loops to path 1) print("Mason's Gain Formula:") print(f"Forward path: P1 = {P1}") print(f"Loop 1: L1 = {L1}") print(f"Loop 2: L2 = {L2}") print(f"Δ = 1 - L1 - L2 = 1 + G2*H1 + G1*G2*G3*H2") print(f"T(s) = G1*G2*G3 / (1 + G2*H1 + G1*G2*G3*H2)") print() print("With numerical values: G1=5, G2=10, G3=2, H1=0.1, H2=0.05") G1n, G2n, G3n = 5, 10, 2 H1n, H2n = 0.1, 0.05 T = (G1n * G2n * G3n) / (1 + G2n*H1n + G1n*G2n*G3n*H2n) print(f"T = {T:.4f}")

Location

System Type

Response Characteristic

LHP real

Stable

Exponential decay

RHP real

Unstable

Exponential growth

LHP complex conjugate

Stable

Damped oscillation

RHP complex conjugate

Unstable

Growing oscillation

jω axis (simple)

Marginally stable

Sustained oscillation

jω axis (repeated)

Unstable

Growing oscillation

Origin (s = 0)

Type N system

N integrators; N pole at origin

import numpy as np import math def second_order_specs(zeta, wn): """Calculate second-order system specifications.""" if zeta < 1: wd = wn * math.sqrt(1 - zeta**2) tp = math.pi / wd Mp = math.exp(-math.pi * zeta / math.sqrt(1 - zeta**2)) * 100 tr = (math.pi - math.atan2(math.sqrt(1-zeta**2), zeta)) / wd else: tp = float('inf') Mp = 0.0 tr = float('inf') ts = 4.0 / (zeta * wn) # 2% criterion return tp, Mp, tr, ts print(f"{'ζ':>5s} {'ω_n':>5s} {'t_p (s)':>8s} {'Mp (%)':>7s} {'t_r (s)':>8s} {'t_s (2%)(s)':>11s}") print("-" * 60) wn = 10 # rad/s for zeta in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.707, 0.8, 1.0, 1.5]: tp, Mp, tr, ts = second_order_specs(zeta, wn) if tp == float('inf'): tp_str = " N/A " else: tp_str = f"{tp:>8.4f}" if tr == float('inf'): tr_str = " N/A " else: tr_str = f"{tr:>8.4f}" print(f"{zeta:>5.3f} {wn:>5.1f} {tp_str} {Mp:>7.2f} {tr_str} {ts:>11.4f}") print(f"\nNote: ζ = 0.707 → Mp = {second_order_specs(0.707, 10)[1]:.2f}% (critically damped, optimal)")

System Type

Step (r=1/s)

Ramp (r=1/s²)

Parabola (r=1/s³)

Type 0 (no pole at origin)

1/(1+Kp)

∞

Type 1 (1 pole at origin)

1/Kv

∞

Type 2 (2 poles at origin)

1/Ka

# Routh-Hurwitz stability criterion def routh_hurwitz(coeffs): """ coeffs: list of polynomial coefficients [a_n, a_{n-1}, ..., a_0] Returns: (stable: bool, routh_table, rhp_poles: int) """ n = len(coeffs) - 1 rows = [] # First two rows from coefficients rows.append(coeffs[0::2]) # even index rows.append(coeffs[1::2]) # odd index # Pad to same length max_len = max(len(rows[0]), len(rows[1])) rows[0] += [0] * (max_len - len(rows[0])) rows[1] += [0] * (max_len - len(rows[1])) while len(rows) < n + 1: prev2 = rows[-2] prev1 = rows[-1] row = [] for j in range(len(prev1) - 1): elem = (prev2[0] * prev1[j+1] - prev1[0] * prev2[j+1]) / prev1[0] row.append(round(elem, 6)) if all(abs(x) < 1e-12 for x in row): row = [0] * len(prev1) rows.append(row) if not row: break # Count sign changes in first column first_col = [r[0] if r else 0 for r in rows] sign_changes = 0 for i in range(1, len(first_col)): if first_col[i] != 0 and first_col[i-1] != 0: if first_col[i] * first_col[i-1] < 0: sign_changes += 1 stable = sign_changes == 0 return stable, rows, sign_changes # Example: s⁴ + 2s³ + 3s² + 4s + 5 = 0 coeffs = [1, 2, 3, 4, 5] stable, table, rhp = routh_hurwitz(coeffs) print(f"Characteristic equation: s⁴ + 2s³ + 3s² + 4s + 5 = 0") print(f"Stable: {stable}") print(f"RHP poles: {rhp}") print("\nRouth array:") for row in table: print([f"{x:>8.4f}" if abs(x) > 1e-12 else " 0" for x in row]) # Stable example: s³ + 6s² + 11s + 6 = 0 (poles: -1, -2, -3) stable2, _, rhp2 = routh_hurwitz([1, 6, 11, 6]) print(f"\ns³ + 6s² + 11s + 6: Stable={stable2}, RHP={rhp2}")

Rule

Statement

#1 Branches

Number of branches = max(n, m); n = poles, m = zeros

#2 Symmetry

Root locus is symmetric about real axis

#3 Real axis segments

Locus exists to the LEFT of odd number of poles+zeros

#4 Asymptotes

Angles: (2k+1)π/(n−m), k=0,1,... Centroid: (Σp−Σz)/(n−m)

#5 Breakaway/Break-in

On real axis; solve dK/ds = 0

#6 Departure angle

From complex pole: θ_dep = 180° + Σ∠(to other poles) − Σ∠(to zeros)

#7 Arrival angle

To complex zero: θ_arr = 180° − Σ∠(to poles) + Σ∠(from other zeros)

#8 Imaginary crossing

Routh array: auxiliary polynomial → jω crossing frequency

#9 K at any point

K = |product(s−pᵢ)| / |product(s−zⱼ)|

import numpy as np import matplotlib.pyplot as plt # Root locus for G(s) = K / (s(s+2)(s+4)) # Characteristic equation: s³ + 6s² + 8s + K = 0 # 1 + K/s(s+2)(s+4) = 0 poles = np.array([0, -2, -4], dtype=float) zeros = np.array([], dtype=float) n = len(poles); m = len(zeros) # Asymptotes n_m = n - m angles = [(2*k + 1) * 180 / n_m for k in range(n_m)] centroid = (poles.sum() - zeros.sum()) / n_m print(f"Poles: {poles}") print(f"Zeros: {zeros}") print(f"Asymptote angles: {angles}°") print(f"Asymptote centroid: {centroid:.2f}") # Breakaway point on real axis (between s=0 and s=-2) # 1/(s-0) + 1/(s+2) + 1/(s+4) = 0 → 3s² + 12s + 8 = 0 a, b, c = 3, 12, 8 disc = b**2 - 4*a*c s1 = (-b + np.sqrt(disc)) / (2*a) # not between 0 and -2 s2 = (-b - np.sqrt(disc)) / (2*a) # between 0 and -2 print(f"Breakaway points: {s1:.3f}, {s2:.3f}") print(f"Valid breakaway: s = {s2:.3f} (between 0 and -2)") # K at breakaway K_break = abs(s2) * abs(s2+2) * abs(s2+4) print(f"K at breakaway = {K_break:.3f}") # Imaginary axis crossing: s = jω in s³+6s²+8s+K=0 # Routh: 6ω²−K = 0 and −ω²(6) + 8 = 0 ... or auxiliary poly from row s¹ # From Routh: row s¹ = [(48-K)/6, 0], crossing when K = 48 # Row s² = 6s² + K = 0 → s² = -K/6 → ω = √(48/6) = 2√2 omega_cross = np.sqrt(48/6) print(f"Imaginary crossing: jω = ±j{omega_cross:.3f}, K = 48") print(f"System unstable for K > 48")

Element

Magnitude (dB)

Phase (degrees)

Gain K

20log₁₀(K)

0°

Integrator 1/s

−20 dB/decade slope

−90°

Differentiator s

+20 dB/decade slope

+90°

Real pole 1/(τs+1)

−20 dB/dec after ω=1/τ

0° → −90°

Real zero (τs+1)

+20 dB/dec after ω=1/τ

0° → +90°

Complex poles

−40 dB/dec after ω_n

0° → −180° (depends on ζ)

Time delay e^(−sT)

0 dB

−ωT × (180/π)°

import numpy as np def bode_magnitude(omega, K, poles_tau, zeros_tau, n_integrators=0): """Compute Bode magnitude in dB.""" mag = K * omega**(-n_integrators) for tau in poles_tau: mag *= 1.0 / np.sqrt(1 + (omega * tau)**2) for tau in zeros_tau: mag *= np.sqrt(1 + (omega * tau)**2) return 20 * np.log10(mag) def bode_phase(omega, poles_tau, zeros_tau, n_integrators=0): """Compute Bode phase in degrees.""" phase = -n_integrators * 90 for tau in poles_tau: phase -= np.arctan(omega * tau) * 180 / np.pi for tau in zeros_tau: phase += np.arctan(omega * tau) * 180 / np.pi return phase # G(s) = 10 / [s(0.5s+1)(0.1s+1)] K = 10 pole_taus = [0.5, 0.1] zero_taus = [] n_int = 1 omega = np.logspace(-2, 3, 1000) mag_db = bode_magnitude(omega, K, pole_taus, zero_taus, n_int) phase_deg = bode_phase(omega, pole_taus, zero_taus, n_int) # Find gain crossover (|G| = 0 dB) idx_gc = np.argmin(np.abs(mag_db)) omega_gc = omega[idx_gc] pm = phase_deg[idx_gc] + 180 # Find phase crossover (phase = -180°) idx_pc = np.argmin(np.abs(phase_deg + 180)) omega_pc = omega[idx_pc] gm = -mag_db[idx_pc] print(f"Gain crossover: ω_gc = {omega_gc:.3f} rad/s") print(f"Phase margin: PM = {pm:.1f}°") print(f"Phase crossover: ω_pc = {omega_pc:.3f} rad/s") print(f"Gain margin: GM = {gm:.1f} dB") print(f"\nCorner frequencies: {', '.join(f'{1/tau:.1f} rad/s' for tau in pole_taus)}")

import numpy as np # State space: Mass-spring-damper system # m ẍ + c ẋ + k x = F # States: x₁ = x, x₂ = ẋ m = 2.0; c = 1.5; k = 10.0 A = np.array([ [0, 1], [-k/m, -c/m] ]) B = np.array([[0], [1/m]]) C = np.array([[1, 0]]) D = np.array([[0]]) # Eigenvalues → poles eigenvalues = np.linalg.eigvals(A) print(f"Eigenvalues (poles): {eigenvalues}") print(f"Natural freq: {np.sqrt(k/m):.3f} rad/s") print(f"Damping ratio: {c/(2*np.sqrt(k*m)):.3f}") # Controllability check n = A.shape[0] U = np.hstack([B] + [A @ B for _ in range(n - 1)]) rank_U = np.linalg.matrix_rank(U) print(f"\nControllability rank: {rank_U} (need {n}) → {'Controllable' if rank_U == n else 'NOT controllable'}") # Observability check V = np.vstack([C] + [C @ np.linalg.matrix_power(A, i) for i in range(1, n)]) rank_V = np.linalg.matrix_rank(V) print(f"Observability rank: {rank_V} (need {n}) → {'Observable' if rank_V == n else 'NOT observable'}") # Transfer function from state space # G(s) = C(sI - A)^(-1)B + D from scipy import signal sys = signal.StateSpace(A, B, C, D) num, den = signal.tf2ss(*signal.ss2tf(A, B, C, D)) tf_sys = signal.TransferFunction(A, B, C, D) print(f"\nTransfer function:") print(f" Numerator: {tf_sys.num}") print(f" Denominator: {tf_sys.den}")

Term

Response Speed

Overshoot

Steady-State Error

Stability

P only

Faster

Increases with Kp

Nonzero (Type 0)

Degrades

Moderate

Increases

Zero (Type 1)

Degrades further

Much faster

Reduces

Nonzero

Improves

PID

Fast

Controlled

Zero

Good (tuned)

import numpy as np def ziegler_nichols_pid(ku, tu, mode='pid'): """Ziegler-Nichols tuning rules. ku = ultimate gain (at stability limit) tu = ultimate period (oscillation period at ku) """ rules = { 'p': {'Kp': 0.5 * ku, 'Ki': 0, 'Kd': 0}, 'pi': {'Kp': 0.45 * ku, 'Ki': 0.54*ku/tu, 'Kd': 0}, 'pid': {'Kp': 0.6 * ku, 'Ki': 1.2*ku/tu, 'Kd': 0.075*ku*tu}, } return rules[mode] # Example: plant with Ku = 4, Tu = 2 seconds Ku = 4.0; Tu = 2.0 for mode in ['p', 'pi', 'pid']: params = ziegler_nichols_pid(Ku, Tu, mode) Ti = params['Kp'] / params['Ki'] if params['Ki'] != 0 else float('inf') Td = params['Kd'] / params['Kp'] if params['Kp'] != 0 else 0 print(f"{mode.upper():>4s}: Kp={params['Kp']:.3f} Ki={params['Ki']:.3f} Kd={params['Kd']:.3f}" f" (Ti={Ti:.3f} Td={Td:.3f})") # Cohen-Coon tuning (for first-order + dead-time: K·e^(-sL)/(τs+1)) # Parameters: K=1, tau=2, L=0.5 K_plant = 1.0; tau = 2.0; L = 0.5 ratio = L / tau # Cohen-Coon PID Kp_cc = (tau / (K_plant * L)) * (4/3 + ratio/4) Ti_cc = L * (32 + 6*ratio) / (13 + 8*ratio) Td_cc = L * (4 / (11 + 2*ratio)) Ki_cc = Kp_cc / Ti_cc Kd_cc = Kp_cc * Td_cc print(f"\nCohen-Coon PID:") print(f" Kp={Kp_cc:.3f} Ki={Ki_cc:.3f} Kd={Kd_cc:.3f}") print(f" Ti={Ti_cc:.3f} Td={Td_cc:.3f}") print(f" (Good for processes with significant dead time)")

Type

Transfer Function

Effect

Application

Lead

Gc = K(1+sτ₁)/(1+sτ₂); τ₁>τ₂

Adds PM, speeds up

Improve transient response

Lag

Gc = K(1+sτ₁)/(1+sτ₂); τ₁<τ₂

Increases low-freq gain

Reduce steady-state error

Lead-Lag

Combination

Both improvements

Comprehensive design

Gc = Kp + Kd·s ≈ lead

Adds derivative action

Reduce overshoot

Gc = Kp + Ki/s ≈ lag

Adds integral action

Eliminate steady-state error