Level 4 · Signals, Systems & Control: Transforms and Stability

Level 2 wrote systems as difference equations; Level 3 measured how a filter treats each frequency. This level joins them. Write z^-1 for "delay by one sample" and a difference equation becomes a fraction of two polynomials, the transfer function H(z). Connecting systems becomes list arithmetic, and putting a number in for z answers physical questions: z = 1 gives the gain for a constant, the unit circle the frequency response, and the roots of the denominator, the poles, say whether a system is stable and how it rings.

The second half completes the controller. Level 3's integral closed P control's gap; the derivative adds a brake, with its own problems of kick and noise. A relay test gives first gains in minutes, and the gain and phase margins say how close a loop is to oscillating and how much delay it can take.

The other paths meet these tools: polynomial products are the Algorithms path's convolution; an autoregressive model in the Statistics path is stationary exactly when its poles lie inside the unit circle; and in the Machine Learning path, too large a learning rate for gradient descent with momentum pushes a pole outside the circle, and training diverges.

Polynomials in the delay z^-1

An FIR filter with weights b is the polynomial

B(z) = b_0 + b_1 z^{-1} + b_2 z^{-2} + \dots = \sum_i b_i \, z^{-i},

and b is its coefficient list, constant term first. Filters side by side (one input, outputs added) add the coefficients of equal powers, a missing one counting as 0: [0.5, 0.5] plus [1, -1] is [1.5, -0.5]. Filters in a row multiply: every term times every term, with z^{-i} z^{-j} = z^{-(i+j)}, so coefficient k of the product is

c_k = \sum_{i + j = k} p_i \, q_j,

Level 2's convolution. The two indices under the sum become the two nested for loops below, each pair adding into out[i + j].

def mul(p, q):                               # a polynomial product: Level 2's convolution
    out = [0.0] * (len(p) + len(q) - 1)
    for i, u in enumerate(p):
        for j, v in enumerate(q):
            out[i + j] += u * v
    return out
smooth, change = [0.5, 0.5], [1, -1]         # the mean of the last two readings; their difference
print(mul(smooth, change), mul(change, change))

[0.5, 0.0, -0.5] [1.0, -2.0, 1.0]

Smoothing then differencing is one filter, half the change over two samples; differencing twice turns positions into an acceleration. (Beware zip(p, q) for a sum: it drops the longer list's tail.)

Practice: Two Echo Units, Side by Side and in a Row

Transfer functions and difference equations

A bathroom scale smooths its display with y[k] = 0.7 * y[k-1] + 0.3 * x[k]. With delays as z^{-1}, Y = 0.7 z^{-1} Y + 0.3 X, so Y (1 - 0.7 z^{-1}) = 0.3 X and

H(z) = \frac{Y}{X} = \frac{B(z)}{A(z)} = \frac{b_0 + b_1 z^{-1} + \dots}{a_0 + a_1 z^{-1} + \dots}

with b = [0.3], a = [1, -0.7]. The convention is \sum_i a_i \, y_{k-i} = \sum_i b_i \, x_{k-i}: output terms move left, so a feedback weight changes sign in a. A radiator heating a room, one sample a minute:

y[k] = 1.3 * y[k-1] - 0.4 * y[k-2] + 0.05 * x[k-1] + 0.05 * x[k-2]
a = [1, -1.3, 0.4]        b = [0, 0.05, 0.05]

The leading zero in b is the minute the hot water needs: a dead time of d samples is d zeros at the front of b. Back again, solve for the newest output,

y_k = \frac{1}{a_0}\left(\sum_{i \ge 0} b_i \, x_{k-i} - \sum_{i \ge 1} a_i \, y_{k-i}\right)

For example, b = [3, 0, 1], a = [2, -1] is y[k] = 1.5 * x[k] + 0.5 * x[k-2] + 0.5 * y[k-1]. Now let z be a number, any complex one but 0, and H(z) is a complex number. For the scale, H(1) = 0.3 / 0.3 = 1: a steady weight is shown exactly. At z = 0.71, H = 0.3 / (1 - 0.7 / 0.71) = 21.3, and at z = 0.7 the denominator is 0: a pole, where H does not exist. Code that evaluates H must check A before dividing.

Common mistakes: copying the feedback signs into a. Delay zeros at the end of b, where they do nothing.

Practice: The Greenhouse Model on Paper, Probing a Hearing-Aid Filter at Any z, Running the Seismometer Filter From Its Datasheet

DC and Nyquist gain

A delay leaves a constant unchanged, so z^{-1} acts as 1 and the DC gain is

H(1) = \frac{\sum_i b_i}{\sum_i a_i},

sum(b) / sum(a), Level 2's steady gain. The fastest signal, +1, -1, \dots at f_s / 2, flips sign at every delay, so the Nyquist gain H(-1) = \sum_i (-1)^i b_i \big/ \sum_i (-1)^i a_i divides the alternating sums b_0 - b_1 + b_2 - \dots. For b = [0.6, 0.2], a = [1, -0.5]: H(1) = 0.8 / 0.5 = 1.6, H(-1) = 0.4 / 1.5 = 0.27, a treble cut. A negative gain inverts; a zero denominator is a pole there, like the integrator a = [1, -1].

Practice: Bass and Treble at a Glance

The frequency response on the unit circle

A sine at f hertz turns by \omega = 2\pi f / f_s radians per sample, and a delay turns it back, multiplying by e^{-j\omega}. So the frequency response is H at z = e^{j\omega},

H(e^{j\omega}) = \frac{\sum_i b_i \, e^{-j\omega i}}{\sum_i a_i \, e^{-j\omega i}},

on the unit circle from z = 1 (0 Hz) to z = -1 (f_s / 2). The gain is |H|, 20\log_{10}|H| in decibels; the phase cmath.phase(H) is negative for a lag. Level 2's wind-gauge display, 10 readings a second:

import cmath, math
b, a, fs = [0.2], [1, -0.8], 10
def response(f):
    back = cmath.exp(-2j * math.pi * f / fs)                     # z^-1 on the unit circle
    return sum(c * back ** k for k, c in enumerate(b)) / sum(c * back ** k for k, c in enumerate(a))
for f in (0, 0.2, 1, 5):
    H = response(f)
    print(f"{f:3} Hz: {20 * math.log10(abs(H)):6.2f} dB, phase {round(math.degrees(cmath.phase(H)), 1) + 0.0:5.1f} degrees")
freqs = [0.02 * 1.1 ** i for i in range(58)]                     # 0.02 to 5 Hz
plot(freqs, [20 * math.log10(abs(response(f))) for f in freqs], logx=True, xlabel="Hz", ylabel="gain (dB)")

  0 Hz:   0.00 dB, phase   0.0 degrees
0.2 Hz:  -1.19 dB, phase -25.9 degrees
  1 Hz:  -9.36 dB, phase -53.1 degrees
  5 Hz: -19.08 dB, phase   0.0 degrees

Press Run for the Bode-style chart: slow changes pass, flicker is cut, -3 dB near 0.36 Hz.

Common mistakes: exp(+1j * w) for z^-1, which flips every phase. 10 * log10 for a gain. math.atan(imag / real), which loses the quadrant. A null, where neither decibels nor phase exist.

Practice: The Subwoofer Filter's Bode Table

Poles and zeros

Multiply A(z) = a_0 + a_1 z^{-1} + \dots + a_N z^{-N} by z^N and it is a_0 z^N + \dots + a_N: the same list, read as powers of z, highest first. Its roots are the poles; the roots of b are the zeros, where H is 0. Factored, with zeros q_i and poles p_i,

H(z) = g \, \frac{(1 - q_1 z^{-1})(1 - q_2 z^{-1}) \cdots}{(1 - p_1 z^{-1})(1 - p_2 z^{-1}) \cdots}

and each factor is the list [1, -root], so chosen poles and zeros multiply out into a and b. Non-real roots of real lists come in conjugate pairs. A zero on the unit circle at angle w removes that frequency; a pole just inside lifts the response near its angle. Degree 2 has the quadratic formula:

import cmath, math
def quadratic_roots(p):                      # p0 z^2 + p1 z + p2 = 0
    root = cmath.sqrt(p[1] ** 2 - 4 * p[0] * p[2])               # complex when negative
    return [(-p[1] + root) / (2 * p[0]), (-p[1] - root) / (2 * p[0])]
zeros, poles = quadratic_roots([1, 0, -1]), quadratic_roots([1, -1.6, 0.89])   # a bat detector
print("poles", [complex(round(p.real, 3), round(p.imag, 3)) for p in poles], "radius", round(abs(poles[0]), 4),
      "angle", round(cmath.phase(poles[0]), 4))
circle = [cmath.exp(2j * math.pi * k / 100) for k in range(101)]
plot([z.real for z in circle], [z.imag for z in circle], label="unit circle")
plot([z.real for z in zeros], [z.imag for z in zeros], kind="scatter", label="zeros")
plot([p.real for p in poles], [p.imag for p in poles], kind="scatter", label="poles")

poles [(0.8+0.5j), (0.8-0.5j)] radius 0.9434 angle 0.5586

Press Run for the pole-zero map. The zeros at +1 and -1 block DC and Nyquist; the poles pass a band around 0.5586 / (2 * pi) of the sample rate, a bat's call near 22 kHz at 250,000 samples a second.

Practice: Placing the Poles and Zeros of a Hum Filter

Finding roots numerically

Beyond degree 4 there is no formula. Horner's rule evaluates a highest-first list with one multiplication per coefficient (start at 0, repeat s = s * z + c), and Newton's method slides a guess down the tangent, z \leftarrow z - P(z) / P'(z): for z^3 - 1.6 z^2 + 1.1 z - 0.3, thirty steps from 1 reach the root 0.6, and from 1j and -1j the pair 0.5 +- 0.5j.

Each start finds one root. The Durand-Kerner method improves guesses for all roots together, without a derivative: each guess z_i moves by P(z_i) (the list first divided by its leading entry) over the product of its distances to the others,

z_i \leftarrow z_i - \frac{P(z_i)}{\prod_{j \ne i} (z_i - z_j)},

a product that math.prod or a loop computes. Were the others already the other roots, that quotient would be exactly z_i - root_i, and one step would land. Start with different, non-real guesses (a real guess stays real, as Newton's start at 1 did), and sweep until the largest step is tiny, with a cap on the sweeps.

Numerical roots are never exact: check them by evaluating P or multiplying the factors (z - root) back together. A double root comes out only to about 1e-8, so tests on roots need tolerances.

Practice: Finding the Poles: a Root Finder in Plain Python

Stability from the poles

The impulse response splits like H: each pole p contributes a multiple of p^k, a repeated pole also k\,p^k.

import cmath, math
w = cmath.exp(1j * math.pi / 5)              # a pair on the circle: a tenth of a turn per sample
terms = {"0.9": lambda k: 0.9 ** k, "1.02": lambda k: 1.02 ** k, "1": lambda k: 1.0,
         "exp(+-j pi/5)": lambda k: 2 * (w ** k).real, "0.9 twice": lambda k: (k + 1) * 0.9 ** k,
         "1 twice": lambda k: k + 1.0}
for name, h in terms.items():
    print(f"pole {name:13}  k = 5: {h(5):6.3f}   k = 300: {h(300):.3g}")

pole 0.9            k = 5:  0.590   k = 300: 1.87e-14
pole 1.02           k = 5:  1.104   k = 300: 380
pole 1              k = 5:  1.000   k = 300: 1
pole exp(+-j pi/5)  k = 5: -2.000   k = 300: 2
pole 0.9 twice      k = 5:  3.543   k = 300: 5.64e-12
pole 1 twice        k = 5:  6.000   k = 300: 301

Every pole inside the unit circle: stable. Even a repeated pole dies away, after rising first.
Any pole outside: unstable.
Poles on the circle, none outside or repeated: marginally stable. A pole at 1 holds a value for ever (an integrator), a pair rings for ever (an oscillator). A repeated pole on the circle grows like k: unstable.

A marginal system is not BIBO stable: feed the integrator a constant, or the oscillator its own frequency, and the output grows without limit. The largest abs(p), the pole radius, says how fast the slowest part fades. A pole on the circle may be computed as 1.0000000001: compare with a tolerance.

Practice: Stable, Unstable or Ringing For Ever?

Resonance and decay from a pole pair

A pair r e^{\pm j\theta} contributes r^k \cos(\theta k + \varphi), an arrow turning \theta per sample and shrinking by r:

f = \frac{\theta f_s}{2\pi} \text{ hertz}, \qquad T_{\text{decay}} = \frac{-1}{f_s \ln r} \text{ seconds},

the time to fall to 1/e (from r^n = e^{-1}), and the pair multiplied out is 1 - 2r\cos\theta \, z^{-1} + r^2 z^{-2}, the list [1, -2 * r * cos(angle), r * r]. A pole near the circle rings long and makes a sharp peak in the frequency response: a resonance. Backwards, the formulas design a sound: a 440 Hz chime at 8000 samples a second, fading to 1/e in 0.25 s, needs angle = 0.3456 and r = exp(-1 / 2000) = 0.9995, so a = [1, -1.8808, 0.999]; run y[k] = 1.8808 * y[k-1] - 0.999 * y[k-2] from one tap and plot it. A real pole only decays (a negative one flips sign every sample); several pairs are several resonances, like a guitar body's.

Practice: The Guitar Body's Resonances, Read From Its Poles

Series and feedback connections

Transfer functions combine like fractions: series multiplies numerators and denominators, parallel gives (b1 * a2 + b2 * a1) / (a1 * a2). For feedback, controller and plant in series form the loop L = C * G = bL / aL; the error R - Y goes in, so Y = L\,(R - Y) and

T = \frac{L}{1 + L} = \frac{b_L}{a_L + b_L}

one polynomial sum (pad the shorter list at the end, then divide both by aT[0]). Its roots are the closed-loop poles; if all are inside the circle, a unit step settles at T(1). An incubator gains 0.1 °C a minute per watt and keeps 95 % of its warmth, under P control with K watts per °C:

bG, aG = [0, 0.1], [1, -0.95]
for K in (1, 5, 10, 21):
    bL = [K * v for v in bG]                 # in series with the gain K (aL = aG)
    aT = [aG[0] + bL[0], aG[1] + bL[1]]      # aT = aL + bL
    y = 0.0
    for k in range(300):                     # the loop itself, asked for one degree more
        y = 0.95 * y + 0.1 * K * (1 - y)
    print(f"K = {K:2}: pole {-aT[1]:5.2f}   T(1) = {sum(bL) / sum(aT):.4f}   simulated {y:.4g}")

K =  1: pole  0.85   T(1) = 0.6667   simulated 0.6667
K =  5: pole  0.45   T(1) = 0.9091   simulated 0.9091
K = 10: pole -0.05   T(1) = 0.9524   simulated 0.9524
K = 21: pole -1.15   T(1) = 0.9767   simulated -1.582e+18

Gain moves the pole from 0.95 (sluggish) past 0 (done in a sample) to beyond -1, where the loop flips sign every minute and explodes, though T(1) still gives a number. A PI controller puts 1 - z^-1 into aL, so sum(aL) = 0 and T(1) = 1: Level 3's zero offset, seen in the algebra.

Practice: Closing the Loop Around a Conveyor Motor

The derivative term: a brake

A plant with inertia arrives at speed: P control pushes until the error is zero, so it overshoots. The derivative term watches the error's slope, a difference over one step:

D_k = K_d \, \frac{e_k - e_{k-1}}{\Delta t}

While the error shrinks, D opposes the P term and brakes before arrival. A sliding door, sent 1 m:

import math
dt = 0.01
for Kd in (0, 1):
    x = v = 0.0
    e_prev, xs = 1.0, []
    for step in range(800):                  # 8 seconds
        e = 1.0 - x
        u = 3 * e + Kd * (e - e_prev) / dt   # Kp = 3, plus the D term
        e_prev = e
        v += dt * (4 * u - 0.5 * v)          # motor push minus friction
        x += dt * v
        xs.append(x)
    zeta = (0.5 + 4 * Kd) / (2 * math.sqrt(4 * 3))
    print(f"Kd = {Kd}: zeta = {zeta:.2f}, overshoot {100 * (max(xs) - 1):4.1f} cm, after 8 s {x:.3f} m")
    plot([k * dt for k in range(800)], xs, label=f"Kd = {Kd}", xlabel="seconds", ylabel="door (m)")

Kd = 0: zeta = 0.07, overshoot 79.7 cm, after 8 s 1.107 m
Kd = 1: zeta = 0.65, overshoot  6.5 cm, after 8 s 1.000 m

With a fixed target the error's slope is minus the door's speed, so D = -Kd * v adds 4 * Kd to the rail's friction: a damper made of software, with Level 2's wn = sqrt(4 * Kp) and zeta = (0.5 + 4 * Kd) / (2 * wn). Once settled the slope is 0, so the final position is unchanged.

Practice: The Gimbal's Braking Term, Stopping a Robot Joint Without a Wobble

Derivative on the measurement

When the setpoint jumps, the error's slope over that sample is the whole jump over dt: the derivative kick. Told to open, the resting door's D term fires 1 * 1.0 / 0.01 = 100 for one sample, 33 times the P term: a hammer blow to the gears. Between setpoint changes the error moves exactly as fast as the measurement, with the opposite sign, so differentiate the measurement:

\text{on the error: } D_k = K_d \, \frac{e_k - e_{k-1}}{\Delta t}, \qquad \text{on the measurement: } D_k = -K_d \, \frac{y_k - y_{k-1}}{\Delta t}

The two agree except at the jump, so braking and the answer to a push from outside are unchanged.

Common mistakes: dropping the minus sign. A first sample with no previous measurement, whose slope must be 0.

Practice: The Drone That Lurched When Told to Climb

A filtered derivative

A derivative divides every change by dt, noise included. The cure is Level 2's first-order smoother on the slope, with time constant T_f:

\text{raw}_k = \frac{y_k - y_{k-1}}{\Delta t}, \qquad D_k = D_{k-1} + \alpha\,(\text{raw}_k - D_{k-1}), \qquad \alpha = \frac{\Delta t}{T_f + \Delta t}

import random, statistics
random.seed(3)
dt = 0.01
pos = [0.3 * k * dt + random.gauss(0, 0.002) for k in range(300)]    # 0.3 m/s, sensor good to 2 mm
for Tf in (0, 0.03, 0.1):
    D, rates = 0.0, []
    for k in range(1, 300):
        D += dt / (Tf + dt) * ((pos[k] - pos[k - 1]) / dt - D)
        rates.append(D)
    print(f"Tf = {Tf:4}: speed {statistics.mean(rates[100:]):.3f} ± {statistics.pstdev(rates[100:]):.3f} m/s")
    plot(rates, label=f"Tf = {Tf}", xlabel="sample", ylabel="m/s")

Tf =    0: speed 0.300 ± 0.276 m/s
Tf = 0.03: speed 0.300 ± 0.050 m/s
Tf =  0.1: speed 0.300 ± 0.017 m/s

Two millimetres of noise make the raw speed wrong by 0.28 m/s; Tf = 0.1 cuts that sixteen times. Press Run for the price: the filtered rate takes about Tf seconds to catch up, and a late brake is weak. Tf is often a fifth to a tenth of the derivative time Td = Kd / Kp.

Practice: Calming the Barometer Before Differentiating It

A full PID update, step by step

Without limits or filter, a PID controller computes

u_k = K_p e_k + K_i \sum_{i=0}^{k} e_i \, \Delta t - K_d \, \frac{y_k - y_{k-1}}{\Delta t}:

a proportional term, a running sum and a difference. It keeps the integral, the filtered derivative and the last measurement between samples, and runs a fixed sequence on each call. This one heats a pottery kiln: three lags of a minute (element, walls, air), and 10 °C in the end for each per cent of power.

class PID:
    def __init__(self, Kp, Ki, Kd, Tf, dt, lo, hi):
        self.k, self.Tf, self.dt, self.lo, self.hi = (Kp, Ki, Kd), Tf, dt, lo, hi
        self.I, self.D, self.y_prev = 0.0, 0.0, None

    def update(self, r, y):
        Kp, Ki, Kd = self.k
        e = r - y                                                         # 1. error
        raw = 0.0 if self.y_prev is None else -(y - self.y_prev) / self.dt   # 2. D on the measurement
        self.y_prev = y
        self.D += self.dt / (self.Tf + self.dt) * (raw - self.D)          # 3. filtered
        I_try = self.I + Ki * e * self.dt                                 # 4. integrate, kept only if
        if self.lo <= Kp * e + I_try + Kd * self.D <= self.hi:            #    the heater could obey
            self.I = I_try
        return max(self.lo, min(self.hi, Kp * e + self.I + Kd * self.D))  # 5. sum and limit

dt, pid = 0.01, PID(0.3, 0.08, 0.25, Tf=0.05, dt=0.01, lo=0, hi=100)
x1 = x2 = y = top = 0.0
for step in range(4000):                     # 40 minutes, from cold to 600 °C above the room
    u = pid.update(600, y)
    x1 += dt * (10 * u - x1)
    x2 += dt * (x1 - x2)
    y += dt * (x2 - y)
    top = max(top, y)
print(f"highest {top:.1f} °C; after 40 minutes {y:.1f} °C on {u:.1f} %, of which the integral gives {pid.I:.1f}")

highest 600.0 °C; after 40 minutes 600.0 °C on 60.0 %, of which the integral gives 60.0

On the first call the P term is 180 %, so the integral step is refused (Level 3's clamping) and the command limited to 100 %; the integral grows only near the target, and ends supplying all 60 %. Every line is a choice: integrating before or after the limit test, clamping the integral itself, updating y_prev before or after use. Each changes the numbers, so a specification states its order, and replaying a log means following it exactly.

Practice: Replaying a 3D Printer's PID Log, Design the Drone's Altitude Hold

Relay tuning and Ziegler-Nichols

The relay test replaces the controller by a switch: d above the working power while the measurement is low, d below while high. The loop settles into an oscillation where the plant lags half a cycle, with the ultimate period Tu. The switch's strongest sine component, 4 * d / pi, becomes a swing of ±a, so the ultimate gain, the P gain at which the loop would oscillate for ever, is about K_u = \frac{4d}{\pi a}.

import math
dt, d = 0.01, 10                             # minutes; relay step, % of power
x1 = x2 = y = 600.0                          # settled at 600 °C on 60 %
temps = []
for step in range(3000):
    u = 60 + d if y < 600 else 60 - d        # the relay
    x1 += dt * (10 * u - x1)
    x2 += dt * (x1 - x2)
    y += dt * (x2 - y)
    temps.append(y)
late = temps[1500:]                          # the last 15 minutes, once steady
a = (max(late) - min(late)) / 2
peaks = [k for k in range(1, len(late) - 1) if late[k - 1] < late[k] >= late[k + 1]]
Tu = (peaks[-1] - peaks[0]) * dt / (len(peaks) - 1)
print(f"swing ±{a:.2f} °C, period {Tu:.3f} min, Ku = {4 * d / (math.pi * a):.3f} % per °C")
plot([k * dt for k in range(3000)], temps, xlabel="minutes", ylabel="°C above the room")

swing ±16.10 °C, period 3.640 min, Ku = 0.791 % per °C

The exact values are 0.8 and 3.628. The Ziegler-Nichols table turns them into gains, with Ki = Kp / Ti and Kd = Kp * Td:

rule    Kp          Ti          Td
P       0.5 * Ku    -           -
PI      0.45 * Ku   Tu / 1.2    -
PID     0.6 * Ku    Tu / 2      Tu / 8

Here P gets Kp = 0.40, PI 0.36 with Ki = 0.12, PID 0.48, 0.26 and 0.22. The rules aim for each swing to be a quarter of the last: lively, a starting point to be calmed.

Practice: Tuning the Mash Tun From a Relay Test

Gain and phase margins

A continuous-time block's response to a sine of w radians per unit time is its formula at s = 1j * w, as z = exp(1j * w * dt) is for sampled ones: a lag 1 / (1 + s * T), an integrator 1 / s (-90 degrees), a derivative s (+90), a delay exp(-s * tau). A wave sent round the loop L returns scaled by abs(L) and shifted by its angle. If it returns upside down at full size, the feedback's minus sign turns it upright and it feeds itself: an oscillation nobody drives. The margins measure the distance from there:

at the gain crossover \omega_c, where |L| falls through 1, the phase margin \text{PM} = 180^\circ + \angle L(\omega_c) is the extra lag the loop could take;
at the phase crossover \omega_p, where the angle falls through -180^\circ, the gain margin \text{GM} = -20\log_{10}|L(\omega_p)| dB is the extra gain it could take.

Scan a log-spaced grid, adding the angles block by block: cmath.phase of the product jumps to +180 once the lag passes 180 degrees. The kiln with Kp = 0.4:

import cmath, math
ws = [0.01 * 1000 ** (k / 1999) for k in range(2000)]            # 0.01 to 10 rad/min
size, angle = [], []
for w in ws:
    blocks = [0.4 * 10] + [1 / (1 + 1j * w * T) for T in (1, 1, 1)]
    size.append(math.prod(abs(c) for c in blocks))
    angle.append(sum(math.degrees(cmath.phase(c)) for c in blocks))
gc = next(k for k in range(1, 2000) if size[k - 1] > 1 >= size[k])
pc = next(k for k in range(1, 2000) if angle[k - 1] > -180 >= angle[k])
print(f"gain crossover {ws[gc]:.3f} rad/min, phase margin {180 + angle[gc]:.1f} degrees")
print(f"phase crossover {ws[pc]:.3f} rad/min, gain margin {-20 * math.log10(size[pc]):.2f} dB")
plot(ws, [20 * math.log10(v) for v in size], logx=True, title="Loop gain (dB)", xlabel="rad/min")
plot(ws, angle, logx=True, title="Loop phase (degrees)", xlabel="rad/min")

gain crossover 1.236 rad/min, phase margin 26.9 degrees
phase crossover 1.734 rad/min, gain margin 6.05 dB

Exactly: 1.233 rad/min with 27.1 degrees, and sqrt(3) with 6.02 dB, where each lag gives 60 degrees and the loop gain is 4 / 8, so the relay's Ku, twice this Kp, would oscillate. Designers like 45 to 60 degrees and 6 dB or more; this loop rings.

Practice: How Close Is the Shower to Oscillating?

How a delay eats phase margin

A delay tau has size 1 and angle -w * tau radians. It leaves the gain crossover alone and takes \omega_c \tau radians off the phase margin, so the loop oscillates once \tau reaches the delay margin \tau_{\max} = \text{PM} / \omega_c, with PM in radians (math.radians(PM) / wc). The kiln's thermocouple, read late:

import math
print(f"delay margin: {math.radians(27.14) / 1.2328:.3f} min")
dt = 0.005
for delay in (0, 0.3, 0.45):                 # minutes
    x1 = x2 = y = 0.0
    seen, ys = [0.0] * (round(delay / dt) + 1), []                # readings on their way
    for step in range(8000):                 # 40 minutes, asked for 100 °C more
        u = 0.4 * (100 - seen[0])
        x1 += dt * (10 * u - x1)
        x2 += dt * (x1 - x2)
        y += dt * (x2 - y)
        seen = seen[1:] + [y]
        ys.append(y)
    print(f"delay {delay:4.2f} min: swing {max(ys[6000:]) - min(ys[6000:]):6.1f} °C in the last 10 minutes")
    plot([k * dt for k in range(0, 8000, 10)], ys[::10], label=f"delay {delay}", xlabel="minutes")

delay margin: 0.384 min
delay 0.00 min: swing    0.2 °C in the last 10 minutes
delay 0.30 min: swing   48.3 °C in the last 10 minutes
delay 0.45 min: swing  373.8 °C in the last 10 minutes

Twenty-three seconds is all this loop can take: press Run, and at 18 the swings shrink, slowly, at 27 they grow. Delays hide in networks, sensor sheaths, derivative filters, even sampling (about half a sample). Lower gain lowers wc, so each second of delay costs fewer degrees, at the price of speed.

Practice: How Much Delay Can the Telescope Take?

P, PI or PID?

The kiln is asked for 100 °C more under the Ziegler-Nichols gains and a calmer PID. IAE, the integrated absolute error, \text{IAE} = \sum_k |e_k| \, \Delta t, adds up how wrong the temperature was, and for how long:

dt = 0.01
gains = {"P": (0.40, 0, 0), "PI": (0.36, 0.12, 0), "PID": (0.48, 0.26, 0.22), "calm": (0.3, 0.08, 0.25)}
for name, (Kp, Ki, Kd) in gains.items():
    x1 = x2 = y = y_prev = I = iae = 0.0
    ys = []
    for step in range(4000):                 # 40 minutes
        e = 100 - y
        I += e * dt
        u = Kp * e + Ki * I - Kd * (y - y_prev) / dt
        y_prev = y
        x1 += dt * (10 * u - x1)
        x2 += dt * (x1 - x2)
        y += dt * (x2 - y)
        ys.append(y)
        iae += abs(e) * dt
    print(f"{name:4}: overshoot {max(ys) - 100:4.1f} °C, after 10 min {ys[999]:5.1f}, after 40 {y:6.2f}, IAE {iae:5.1f}")
    plot([k * dt for k in range(4000)], ys, label=name, xlabel="minutes", ylabel="°C above the start")

P   : overshoot 23.0 °C, after 10 min  72.4, after 40  80.01, IAE 894.5
PI  : overshoot 55.7 °C, after 10 min  77.2, after 40  99.35, IAE 469.5
PID : overshoot 50.3 °C, after 10 min  98.8, after 40 100.00, IAE 213.9
calm: overshoot  4.2 °C, after 10 min 101.5, after 40 100.00, IAE 177.6

P stops 1 / (1 + 4) of the step short. PI removes the offset, but the integral adds lag, so it swings for half an hour. The derivative gives phase back near the crossover, so PID is faster and better damped, and calmer gains turn 50 °C of overshoot into 4. Use P where an offset does not matter; PI for most temperatures, flows and pressures, where noise makes a derivative more trouble than it is worth; PID where inertia or a chain of lags makes the brake pay: motors, kilns, drones.

Practice: P, PI or PID for the Lifting Arm?

Practice

These problems are required on this path at this level; the path adds one problem from each of the level's topics, picked for you and kept.

Checklist

You can now:

Add and multiply coefficient lists, and turn a difference equation into b and a and back.
Evaluate H(z) anywhere: DC and Nyquist gains, decibels and phase on the unit circle.
Find and check poles and zeros, explain Durand-Kerner, and draw a pole-zero map.
Judge stability from the poles, and read a pair's frequency and decay.
Build a closed loop bL / (aL + bL) and find its poles and final value.
Use the derivative as damping, on the measurement and filtered.
Write a PID update in a stated order, tune it from a relay test, and compare P, PI and PID.
Find gain, phase and delay margins on a frequency grid.