Level 05 · Signals, Systems & Control

Feedback Control

The FFT, filter design, root locus and loop stability.

Level 5 · Signals, Systems & Control: Feedback Control

Level 4 described a system by its transfer function and judged it by its poles. This level designs with them. The first half builds filters: the fast Fourier transform, which makes long spectra cheap; low-pass and band-pass filters from the sinc recipe; recursive biquads, which do similar work with five multiplications a sample; a notch for mains hum; and an honest measurement of what a filter gained.

The second half closes loops on purpose. A gain moves the closed loop's poles; tabulated over many gains they form the root locus, and bisection finds the gain at which the loop breaks. A pole's radius and angle predict settling and overshoot. Then come cascade loops, feedforward of a measured disturbance, and a computer that samples, holds and answers late, its state written as a vector so that one matrix holds the poles of the whole machine.

The other paths meet these ideas. The FFT is the Algorithms path's divide and conquer, with merge sort's recurrence and O(n log n). Finding the largest stable gain is binary search on a real number. And the Machine Learning path's models, built there from plain lists, run on numpy arrays almost everywhere else.

numpy: arrays instead of lists

From this level, and on this path only, programs may import numpy as np. It brings speed, since one operation runs a compiled loop over a whole array, often tens of times faster than Python, and linear algebra: matrix products, eigenvalues and polynomial roots, one call each. None of it is new: np.convolve is Level 2's convolution, np.fft Level 3's DFT, np.roots Level 4's root finder. An array holds numbers of one type, and arithmetic acts on every element:

import numpy as np
depths = [1.20, 1.35, 1.10, 1.42]             # a river, metres, one reading a day
a = np.array(depths)
print(depths + depths, (a + a).tolist())      # a list joins, an array adds
print((100 * a - 50).round(1).tolist(), float(a.mean()), int(np.argmax(a)), a[a > 1.3].tolist())
M = np.array([[1.0, 0.5], [0.0, 1.0]])
v = np.array([2.0, 4.0])
print((M @ v).tolist(), (M * v).tolist())     # matrix product; element by element

[1.2, 1.35, 1.1, 1.42, 1.2, 1.35, 1.1, 1.42] [2.4, 2.7, 2.2, 2.84]
[70.0, 85.0, 60.0, 92.0] 1.2675 3 [1.35, 1.42]
[4.0, 4.0] [[2.0, 2.0], [0.0, 4.0]]

a > 1.3 is an array of True and False that selects elements. @ is the matrix product; * works element by element. Print through .tolist() or float(...): numpy's own numbers print as np.float64(1.2675).

Common mistakes: * where @ was meant. Looping over an array element by element, which throws the speed away.

The FFT: divide and conquer

Level 3's DFT, X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi j k n / N}, costs N^2 products. The fast Fourier transform halves the job: transform the even-numbered samples (E) and the odd-numbered ones (O) separately, each of length N / 2, then combine

X_k = E_k + e^{-2\pi j k / N} O_k, \qquad X_{k + N/2} = E_k - e^{-2\pi j k / N} O_k, \qquad k = 0, \dots, N/2 - 1

The odd samples sit one step later, which turns their arrows by the twiddle factor; half a turn on, it flips sign, so one product serves two bins. Checked on bins 1 and 5 of eight:

import cmath, math
dft = lambda x: [sum(v * cmath.exp(-2j * math.pi * k * n / len(x)) for n, v in enumerate(x)) for k in range(len(x))]
x = [3, 1, 4, 1, 5, 9, 2, 6]
X, E, O = dft(x), dft(x[0::2]), dft(x[1::2])
t = cmath.exp(-2j * math.pi * 1 / 8) * O[1]   # the twiddle for k = 1
print(abs(X[1] - (E[1] + t)) < 1e-12, abs(X[5] - (E[1] - t)) < 1e-12)

True True

Each half splits again, down to single samples, whose transform is the sample itself. For N a power of two there are log2(N) levels of N / 2 products each: the cost obeys merge sort's recurrence T(N) = 2\,T(N/2) + N, so it is O(N \log N). For N = 65536, half a million products instead of 4.3 billion.

Common mistakes: splitting into first and second halves instead of evens and odds. Using the full length in the twiddle at a smaller level.

Practice: The Gearbox Analyser's Fast Transform

numpy.fft and zero padding

np.fft.rfft(x) returns bins 0 to N // 2 of a real signal of any length, and np.fft.rfftfreq(N, 1 / fs) their frequencies; ten minutes of a 100 Hz recording take milliseconds. Amplitudes are 2|X_k| / N, 2 * np.abs(X) / N as in Level 3, or 2|X_k| / \sum_n w_n with a window w, since each sample then counts only as much as its weight.

Zero padding appends zeros before transforming: np.fft.rfft(x, n=L) pads to L samples. The zeros add nothing to any sum, so the underlying spectrum is unchanged, but it is shown every fs / L hertz instead of fs / N. A tuning fork, heard for 50 ms:

import numpy as np
fs, N = 8000, 400                             # bins 20 Hz apart
fork = np.cos(2 * np.pi * 443 * np.arange(N) / fs)
for L in (N, 1600, 16384):
    mags = np.abs(np.fft.rfft(fork, n=L))
    print(f"padded to {L:5}: grid {fs / L:6.3f} Hz, peak at {np.argmax(mags) * fs / L:6.2f} Hz")
    plot(np.fft.rfftfreq(L, 1 / fs), 2 * mags / N, label=f"{L}", kind="stem" if L == N else "line", xlabel="Hz")

padded to   400: grid 20.000 Hz, peak at 440.00 Hz
padded to  1600: grid  5.000 Hz, peak at 445.00 Hz
padded to 16384: grid  0.488 Hz, peak at 442.87 Hz

Press Run: the padded curves pass through the stems and fill the gaps between them. Padding locates one peak precisely, but adds no resolution: two forks closer than fs / N still make one lump, and only a longer recording separates them.

Common mistakes: dividing by N when the padded bins are fs / L apart. Forgetting that argmax of a slice counts from the slice's start.

Practice: Tuning a Guitar String with Zero Padding, A Minute of Vibration from the Engine Room

Windowed-sinc low-pass design

Design asks Level 3's question backwards. The ideal low-pass (gain 1 below the cutoff f_c, 0 above) has a sinc for its impulse response, at t samples from the centre

h_t = \frac{\sin(2\pi f_c t / f_s)}{\pi t} = \frac{2 f_c}{f_s} \operatorname{sinc}\left(\frac{2 f_c t}{f_s}\right), \qquad \operatorname{sinc}(u) = \frac{\sin(\pi u)}{\pi u},

with h_0 = 2 f_c / f_s at the centre, rippling away for ever. A real filter keeps taps weights around the centre (t = k - (taps - 1) / 2), tapers them with a window, since a sinc cut off abruptly leaves ripples in the response, and divides them by their sum so that a constant passes unchanged. numpy has the sinc as np.sinc(u) and the usual windows; the zero-padded FFT of the weights is their response. A load cell under an e-bike's pedals:

import numpy as np
fs, fc, taps = 1000, 50, 61
t = np.arange(taps) - (taps - 1) / 2
ideal = 2 * fc / fs * np.sinc(2 * fc * t / fs)
f = np.fft.rfftfreq(10000, 1 / fs)            # a 0.1 Hz grid
for name, w in (("rectangular", np.ones(taps)), ("hamming", np.hamming(taps)), ("blackman", np.blackman(taps))):
    h = ideal * w / np.sum(ideal * w)
    db = 20 * np.log10(np.maximum(np.abs(np.fft.rfft(h, n=10000)), 1e-12))
    print(f"{name:11}  at 30 Hz {db[300]:5.2f} dB   above 100 Hz at most {db[f >= 100].max():5.1f} dB"
          f"   -40 dB from {f[np.argmax(db < -40)]:4.1f} Hz")
    plot(f[1:], db[1:], label=name, logx=True, xlabel="Hz", ylabel="gain (dB)")

rectangular  at 30 Hz  0.18 dB   above 100 Hz at most -33.5 dB   -40 dB from 60.0 Hz
hamming      at 30 Hz -0.36 dB   above 100 Hz at most -64.7 dB   -40 dB from 74.9 Hz
blackman     at 30 Hz -0.89 dB   above 100 Hz at most -76.2 dB   -40 dB from 83.7 Hz

The window is a trade. Without one (rectangular) the gain dives fastest, but its ripples climb back to -33.5 dB; Hamming stays below -64 dB; Blackman goes deeper, with the widest edge. More taps narrow the edge in proportion.

Common mistakes: taps where taps - 1 belongs, which makes the window lopsided. Testing k == taps // 2 for the centre, wrong for even taps.

Practice: A Low-Pass Filter from the Sinc Recipe

Band-pass from two low-passes

Filters of the same length can be subtracted weight by weight, and since a response is a sum over the weights, their responses subtract too. A low-pass at f_high minus one at f_low keeps what the first keeps and the second does not: a band-pass. An earthquake-warning seismometer at 100 Hz that wants 1 to 10 Hz:

import numpy as np
fs, taps = 100, 101
def lowpass(fc):
    t = np.arange(taps) - (taps - 1) / 2
    h = 2 * fc / fs * np.sinc(2 * fc * t / fs) * np.blackman(taps)
    return h / h.sum()
band = lowpass(10) - lowpass(1)
H = np.abs(np.fft.rfft(band, n=1000))         # a 0.1 Hz grid
print(", ".join(f"{f} Hz {20 * np.log10(H[round(f * 10)]):.1f} dB" for f in (0.2, 3, 6, 25)))

0.2 Hz -36.6 dB, 3 Hz -0.1 dB, 6 Hz 0.0 dB, 25 Hz -105.5 dB

Below 1 Hz both pass and cancel; above 10 Hz neither passes. The weights sum to 0, so the gain at 0 Hz is minus infinity decibels: printing decibels needs a floor. A single 1 at the centre (a pure delay) minus a low-pass is a high-pass.

Practice: A Band-Pass for the Birdsong Recorder

IIR biquads and their cost

A sharp FIR edge at a low frequency needs many weights; feedback is cheaper. A biquad is Level 4's difference equation with two of each kind of memory,

y_n = b_0 x_n + b_1 x_{n-1} + b_2 x_{n-2} - a_1 y_{n-1} - a_2 y_{n-2},

with transfer function H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}, and standard recipes turn a corner frequency and a sharpness Q (0.7071 for the flattest pass band) into its coefficients, each divided by a0. A drone's microphone must lose the rotors' rumble below 100 Hz:

import math, cmath
fs, fc, q = 16000, 100, 0.7071
w0 = 2 * math.pi * fc / fs
c, alpha = math.cos(w0), math.sin(w0) / (2 * q)
b = [(1 + c) / 2, -(1 + c), (1 + c) / 2]      # the high-pass recipe
a = [1 + alpha, -2 * c, 1 - alpha]
b, a = [v / a[0] for v in b], [v / a[0] for v in a]
for f in (25, 50, 100, 400):
    z1 = cmath.exp(-2j * math.pi * f / fs)    # z^-1 on the unit circle
    H = sum(v * z1 ** k for k, v in enumerate(b)) / sum(v * z1 ** k for k, v in enumerate(a))
    print(f"{f} Hz: {round(20 * math.log10(abs(H)), 1) + 0.0} dB")

25 Hz: -24.1 dB
50 Hz: -12.3 dB
100 Hz: -3.0 dB
400 Hz: 0.0 dB

The gain is -3 dB at the corner and falls 12 dB per octave below it. A windowed-sinc FIR needs roughly attenuation_dB * fs / (22 * transition_Hz) weights: 233 for the same -24 dB at 25 Hz, or 117 multiplications a sample if symmetric pairs share one. The biquad needs five; a steeper corner chains a few. The prices: the poles must stay inside the unit circle (rounding in a small chip can push them out), and different frequencies are delayed by different amounts.

Common mistakes when running one: adding the feedback terms instead of subtracting them; forgetting to divide by a0; shifting the memory before computing y.

Practice: A Biquad Low-Pass for the Crane's Tilt Sensor, Multiplications per Sample in a Hearing Aid

Group delay and aligning outputs

Every causal smoother is late. The group delay at a frequency is how many samples late a wave near it, and the bump it carries, comes out: minus the rate at which the phase turns as the frequency rises. A wave buoy's height sensor:

import cmath, math
fs = 10
def response(b, a, f):
    z1 = cmath.exp(-2j * math.pi * f / fs)
    return sum(v * z1 ** k for k, v in enumerate(b)) / sum(v * z1 ** k for k, v in enumerate(a))
def delay(b, a, f, df=1e-6):                  # samples
    return -cmath.phase(response(b, a, f + df) / response(b, a, f)) / (2 * math.pi * df / fs)
tri = [v / 25 for v in (1, 2, 3, 4, 5, 4, 3, 2, 1)]       # symmetric FIR
for f in (0, 0.25, 0.5, 1):
    print(f"{f:4} Hz: triangle {delay(tri, [1], f):4.2f}, EMA {delay([0.2], [1, -0.8], f):4.2f} samples")

   0 Hz: triangle 4.00, EMA 4.00 samples
0.25 Hz: triangle 4.00, EMA 2.52 samples
 0.5 Hz: triangle 4.00, EMA 1.02 samples
   1 Hz: triangle 4.00, EMA 0.02 samples

Symmetric weights always delay by exactly (N - 1) / 2 samples: linear phase, which keeps a waveform's shape. Level 2's exponential average is as late for a slow swell but hardly late at 1 Hz, so it bends a wave group out of shape; ECG monitors often use FIR filters for this reason.

A constant delay can be undone when the whole recording is at hand: drop D = (N - 1) / 2 samples from the front of the full convolution and keep the original length, which is what np.convolve(x, h, mode="same") does for odd N. Compare a filtered signal with the original only after aligning, or the error measured is mostly the delay; and quote delays in the units asked, D / fs * 1000 for milliseconds. A live system lives with D.

Practice: How Late Is the Smoothed Pulse?

A notch for hum

Level 4 placed zeros and poles by hand to remove a frequency. The standard notch biquad puts two zeros on the unit circle at the hum's angle and two poles at the same angle, radius r, just inside. With c = \cos(2\pi f_0 / f_s):

H(z) = g \, \frac{1 - 2c\,z^{-1} + z^{-2}}{1 - 2rc\,z^{-1} + r^2 z^{-2}}, \qquad g = \frac{1 - 2rc + r^2}{2 - 2c}

The numerator's zeros sit on the circle, so the hum's gain is 0; the poles at radius r fill the notch in nearby; and g makes the gain at 0 Hz, H(1), exactly 1. As lists, B = g * [1, -2 * c, 1] and A = [1, -2 * r * c, r * r].

Away from the hum the poles nearly cancel the zeros, so the closer r is to 1, the narrower the notch. A strain gauge on a bridge picks up 50 Hz mains:

import numpy as np
fs, f0 = 1000, 50
c = np.cos(2 * np.pi * f0 / fs)
f = np.linspace(0, fs / 2, 50001)             # a 0.01 Hz grid
z1 = np.exp(-2j * np.pi * f / fs)
for r in (0.95, 0.99):
    g = (1 - 2 * r * c + r * r) / (2 - 2 * c)
    H = np.abs(g * (1 - 2 * c * z1 + z1 ** 2) / (1 - 2 * r * c * z1 + r * r * z1 ** 2))
    notch = f[H < 1 / np.sqrt(2)]             # inside the -3 dB notch
    print(f"r = {r}: gain at 50 Hz {H[5000]:.3f}, at 45 Hz {H[4500]:.3f};"
          f" width {notch.max() - notch.min():.2f} Hz, rule {(1 - r) * fs / np.pi:.2f} Hz")
    plot(f[2000:8001], 20 * np.log10(H[2000:8001] + 1e-9), label=f"r = {r}", xlabel="Hz", ylabel="dB")

r = 0.95: gain at 50 Hz 0.000, at 45 Hz 0.534; width 15.60 Hz, rule 15.92 Hz
r = 0.99: gain at 50 Hz 0.000, at 45 Hz 0.953; width 3.18 Hz, rule 3.18 Hz

The rule of thumb, a width of about (1 - r) f_s / \pi hertz, suits narrow notches; an exact width means solving for r, for instance by bisection. Mains hum has harmonics at 100, 150, ... Hz, so a hum filter chains one notch per harmonic below fs / 2, each fed the last one's output. A narrow notch rings at switch-on for several times 1 / (1 - r) samples.

Common mistakes: forgetting g, which scales the baseline. Carrying one notch's memory into the next.

Practice: Notching the Mains Hum out of an ECG

Measuring noise reduction in decibels

On a test bench, where the true signal clean is known, a filter's benefit is a ratio of error powers, the mean squares of noisy - clean and filtered - clean:

\text{gain} = 10\log_{10}\frac{\frac{1}{N}\sum_n (\text{noisy}_n - \text{clean}_n)^2}{\frac{1}{N}\sum_n (\text{filtered}_n - \text{clean}_n)^2} \text{ dB}

(powers, so 10, not 20). The filtered error includes damage to the signal, so too much smoothing scores badly. Align first, and leave out the edges, where the filter was fed zeros. A rowing machine's handle force:

import numpy as np, random
rng = random.Random(5)
N = 2400                                      # 2 minutes at 20 Hz
clean = 300 + 250 * np.sin(2 * np.pi * 0.45 * np.arange(N) / 20)    # 27 strokes a minute
noisy = clean + np.array([rng.gauss(0, 40) for _ in range(N)])
for w in (3, 7, 15, 31):                      # moving averages of w samples
    out = np.convolve(noisy, np.ones(w) / w, mode="same")
    keep = slice(w // 2, N - w // 2)
    ratio = np.mean((noisy - clean)[keep] ** 2) / np.mean((out - clean)[keep] ** 2)
    print(f"w = {w:2}: {10 * np.log10(ratio):5.2f} dB")

w =  3:  4.74 dB
w =  7:  7.20 dB
w = 15:  1.38 dB
w = 31: -9.08 dB

Averaging 7 samples divides the noise power by 7 (8.5 dB) but starts to flatten the strokes; at 31 the "filtered" force is further from the truth than the raw one. Only a measurement finds the best compromise.

Practice: Tuning the Step Counter's Filter in Decibels

The closed loop's poles

Level 4 closed a loop around a conveyor and found its denominator, aL + bL. With a proportional gain K and a sensor that reports d samples late, u = K * (r - z^-d y), and the characteristic polynomial is

A(z) + K z^{-d} B(z),

where z^{-d} B is the list b with d zeros in front.

Its roots, the closed-loop poles, decide whether the loop settles, rings or runs away. A list in powers of z^-1, constant first and padded to a common length, is also the polynomial in z read highest power first, as np.roots expects. A lift's speed drive:

import numpy as np
A = np.array([1, -1.5, 0.54, 0])              # the drive alone, padded to four entries
B = np.array([0, 0.05, 0.04, 0])
B1 = np.array([0, 0, 0.05, 0.04])             # B one sample later: z^-1 B
for name, p in (("plant", A[:3]), ("K = 4", (A + 4 * B)[:3]), ("K = 4, d = 1", A + 4 * B1)):
    radii = np.sort(np.abs(np.roots(p)))[::-1]
    print(f"{name:12} {np.round(p, 2).tolist()}  radii {np.round(radii, 3).tolist()}")

plant        [1.0, -1.5, 0.54]  radii [0.9, 0.6]
K = 4        [1.0, -1.3, 0.7]  radii [0.837, 0.837]
K = 4, d = 1 [1.0, -1.5, 0.74, 0.16]  radii [1.002, 1.002, 0.159]

The gain pulls the plant's poles, 0.9 and 0.6, into a ringing pair. One sample of delay adds a third pole and pushes the pair just outside the unit circle. (A trailing 0 is a pole at the origin, which is harmless.)

Common mistakes: multiplying A by K too. Putting the delay's zeros at the end, an advance. Comparing real parts with 1 instead of radii.

Practice: Is the Drone's Height Loop Stable?, A Gain Around the Loop

Root locus as a table

Each gain gives new poles; their paths over all gains are the root locus. A loop over gains draws it, and a table summarises it by the slowest pole:

import numpy as np
a, b = np.array([1, -1.5, 0.54]), np.array([0, 0.05, 0.04])
re, im = [], []
for K in np.arange(0, 20, 0.1):
    z = np.roots(a + K * b)
    re, im = re + z.real.tolist(), im + z.imag.tolist()
plot(re, im, kind="scatter", xlabel="real part", ylabel="imaginary part")
for K in (0, 0.5, 4, 11, 12):
    z = max(np.roots(a + K * b), key=abs)
    print(f"K = {K:4}: radius {abs(z):.3f}, angle {abs(np.degrees(np.angle(z))):4.1f} degrees")

K =    0: radius 0.900, angle  0.0 degrees
K =  0.5: radius 0.748, angle  9.8 degrees
K =    4: radius 0.837, angle 39.0 degrees
K =   11: radius 0.990, angle 61.3 degrees
K =   12: radius 1.010, angle 63.5 degrees

Press Run. The real poles slide together, meet near K = 0.29 and split into a pair whose angle (faster swings) and radius (slower decay) grow until it leaves the unit circle between 11 and 12. The quickest settling comes where the poles meet, far below the largest stable gain.

Practice: The Root Locus of a Lawnmower's Steering

The largest stable gain

Stability is a yes-or-no question whose answer here changes once, so bisection finds the switch: keep a stable lo and an unstable hi, test the middle, move the end with the same verdict.

import numpy as np
A = np.array([1, -1.5, 0.54, 0])
for name, Bd in (("no delay", np.array([0, 0.05, 0.04, 0])), ("d = 1", np.array([0, 0, 0.05, 0.04]))):
    lo, hi = 0.0, 100.0
    for _ in range(50):                       # each test halves the gap
        mid = (lo + hi) / 2
        if max(abs(np.roots(A + mid * Bd))) < 1:
            lo = mid
        else:
            hi = mid
    print(f"{name:8}: critical gain {lo:.6f}, gap {hi - lo:.0e}")

no delay: critical gain 11.500000, gap 9e-14
d = 1   : critical gain 3.954370, gap 9e-14

Without delay the pair's radius squared is 0.54 + 0.04 * K, which reaches 1 at exactly 11.5; one sample of delay leaves a third of that. Loops run well below the critical gain: half of it is a gain margin of 2, or 6 dB.

Practice: How Much Gain Can the Robot Arm Take?

Reading a pole

A pole z = r e^{j\theta} shrinks its part of the response by r and turns it by \theta every sample. So, in samples (times \Delta t for seconds):

k_{\text{settle}} = \frac{\ln 0.02}{\ln r} \ \text{(from } r^k = 0.02\text{)}, \qquad k_{\text{swing}} = \frac{2\pi}{\theta}, \qquad \text{overshoot} \approx 100\, r^{\pi/\theta} \ \%,

the overshoot being what is left half a swing later.

They are exact for one pair of poles, and good estimates when one pair is much slower than the rest. The lift at K = 4, against a simulated step:

import cmath, math
dt = 0.05                                     # 20 samples a second
z = complex(0.65, math.sqrt(0.7 - 0.65 ** 2)) # a root of z^2 - 1.3 z + 0.7
r, th = abs(z), abs(cmath.phase(z))
print(f"predicted: settles in {math.log(0.02) / math.log(r) * dt:.2f} s, swings every {2 * math.pi / th * dt:.2f} s,"
      f" overshoot {100 * r ** (math.pi / th):.0f} %")
y = [0.0, 0.2]                                # (A + 4 B) y = 4 B r, with r = 1 from sample 0
for k in range(2, 80):
    y.append(1.3 * y[-1] - 0.7 * y[-2] + 4 * (0.05 + 0.04))
late = [k for k, v in enumerate(y) if abs(v - y[-1]) > 0.02 * y[-1]][-1] + 1
print(f"simulated: settles in {late * dt:.2f} s, overshoot {100 * (max(y) / y[-1] - 1):.0f} %")
plot([k * dt for k in range(80)], y, kind="step", xlabel="seconds", ylabel="speed (m/s)")

predicted: settles in 1.10 s, swings every 0.46 s, overshoot 44 %
simulated: settles in 1.00 s, overshoot 42 %

The zero in B, which the rules ignore, explains the small differences. A positive real pole (theta = 0) never swings; a negative one (theta = pi) flips sign every sample.

Practice: What a Gimbal's Pole Says About Its Swing

Sampling too slowly

A computer that reads every Ts seconds and holds its command until the next reading is a zero-order hold. A held pump flow raises a tank's level in a straight line, so one step is exact, and with P control the loop has one pole, 1 - K T_s / \text{area}:

area, K = 2.0, 0.5                            # a 2 m² tank; the pump gives K * (1.2 - h) m³/s
for Ts in (0.5, 4, 6, 9):                     # seconds between readings
    h, levels = 1.0, []
    for k in range(4):
        h = h + Ts * K * (1.2 - h) / area     # the flow is held for Ts seconds
        levels.append(round(h, 3))
    print(f"Ts = {Ts:3}: pole {1 - K * Ts / area:5.2f}  levels {levels}")

Ts = 0.5: pole  0.88  levels [1.025, 1.047, 1.066, 1.083]
Ts =   4: pole  0.00  levels [1.2, 1.2, 1.2, 1.2]
Ts =   6: pole -0.50  levels [1.3, 1.15, 1.225, 1.188]
Ts =   9: pole -1.25  levels [1.45, 0.887, 1.591, 0.712]

At 4 s one correction lands exactly; at 6 s each overshoots and the level zigzags in; beyond 2 * area / K = 8 s each zigzag outgrows the last, although measured continuously this loop is stable at any gain. For a first-order plant (time constant \tau, gain G) the exact held step is

y_{k+1} = a\,y_k + b\,u_k, \qquad a = e^{-T_s/\tau}, \qquad b = G\,(1 - a),

and the loop fails once a - bK passes -1. Sample 10 to 20 times per closed-loop time constant.

Practice: A Hot Plate Sampled Too Slowly

State as a vector

A mass on a spring stores two quantities, position and speed. As a state vector s, its rate of change is A s + B u for a matrix A. Holding u for T seconds gives the exact step

s_{k+1} = A_d \, s_k + B_d \, u_k,

where A_d is the matrix exponential of AT, the series I + AT + \frac{(AT)^2}{2!} + \frac{(AT)^3}{3!} + \dots with matrix products for powers, and B_d follows from the same series. Many tiny Euler steps, (I + AT/n)^n, reach the same matrix, and u as an extra, unchanging state gives Bd too. Feedback u_k = -K s_k, with a row of gains K, makes the loop one matrix, s_{k+1} = (A_d - B_d K)\, s_k, whose eigenvalues, np.linalg.eigvals(M), are its poles; a command applied a sample late is one more state, waiting. Level 2's mass-spring-damper as an isolation table:

import numpy as np
m, k, c, T = 50, 5000, 100, 0.02              # kg, N/m, N·s/m; a computer every 0.02 s
big = np.array([[0, 1, 0], [-k / m, -c / m, 1 / m], [0, 0, 0]])   # [[A, B], [0, 0]]
step = np.linalg.matrix_power(np.eye(3) + big * T / 2 ** 20, 2 ** 20)
Ad, Bd = step[:2, :2], step[:2, 2]
print("Ad =", np.round(Ad, 4).tolist(), " Bd =", [f"{v:.3e}" for v in Bd])
gains = np.array([15000, 1900])               # PD: u = -(15000 * height + 1900 * speed)
late = np.zeros((3, 3))                       # state [height, speed, waiting u]
late[:2, :2], late[:2, 2], late[2, :2] = Ad, Bd, -gains
for name, M in (("plant", Ad), ("PD at once", Ad - np.outer(Bd, gains)), ("PD one late", late)):
    z = np.linalg.eigvals(M)
    print(f"{name:11} radii {np.round(np.sort(np.abs(z))[::-1], 3).tolist()}, largest angle {np.degrees(np.abs(np.angle(z))).max():.1f}")

Ad = [[0.9803, 0.0195], [-1.9475, 0.9414]]  Bd = ['3.934e-06', '3.895e-04']
plant       radii [0.98, 0.98], largest angle 11.4
PD at once  radii [0.751, 0.371], largest angle 0.0
PD one late radii [0.931, 0.931, 0.786], largest angle 52.4

PD control applied at once gives two real poles and no ringing; the same gains a sample late give a pair at radius 0.931 and 52 degrees, an overshoot near 78 %. A model without the waiting command hides this. A ball rolling on a tilting beam is the same idea, with the tilt as the state that waits.

Practice: Sampling a Motor with a Matrix Exponential, Tune a Ball on a Beam by Its Poles

Cascade loops

When a slow quantity depends on a faster one that can also be measured, two loops beat one. In cascade control the outer loop watches the slow quantity and asks for a value of the fast one; the inner loop delivers it. A village's header tank: the level controller asks for a flow, and a flow controller moves the valve until the meter agrees. After two minutes the supply weakens by 40 %:

def run(cascade, dt=0.1, area=2.0):           # level wanted 1.5 m; 0.05 m³/s drawn off
    h, v, I, Iq, low = 1.5, 0.25, 0.0, 0.0125, 1.5
    for step in range(6000):                  # 10 minutes
        q = (0.2 if step < 1200 else 0.12) * v          # the flow meter
        I += (1.5 - h) * dt
        if cascade:
            q_ref = 0.05 + 0.4 * (1.5 - h) + 0.004 * I  # outer: the level asks for a flow
            Iq += (q_ref - q) * dt
            v_cmd = 5 * (q_ref - q) + 20 * Iq           # inner: PI on the flow
        else:
            v_cmd = 0.25 + 2 * (1.5 - h) + 0.02 * I     # one loop: the level moves the valve
        v += dt * (min(max(v_cmd, 0), 1) - v)           # the valve lags by 1 s
        h += dt * (q - 0.05) / area
        low = min(low, h)
    return low
print(f"lowest level: single loop {run(False):.3f} m, cascade {run(True):.3f} m")

lowest level: single loop 1.427 m, cascade 1.494 m

The outer gains match (2 per metre on a valve worth 0.2 m³/s is 0.4 m³/s per metre). The single loop notices the weak supply only once the level falls; the inner loop sees the flow drop and opens the valve before the tank feels it. The rule: the inner loop must be much faster than the outer, five to ten times, so that the outer loop can treat a requested flow as delivered. Here they settle in about one and five seconds; a sluggish inner loop makes the pair ring.

Practice: Cascade Control of a Camera Slider

Feedforward of a measured disturbance

Feedback acts only after an error appears. If the disturbance is measured and a model says what it costs, the controller can pay in advance: feedforward. A cruise control reads the slope from its map; holding speed on a grade s takes an extra m g s newtons (g = 9.81 m/s²), so the force is F = F_{\text{feedback}} + m g s:

import math
dt, m, drag = 0.1, 1500, 40                   # kg, N per m/s; cruise at 25 m/s
grade = [0.0] * 300 + [0.05] * 600 + [-0.03] * 600 + [0.0] * 300
def drive(m_model):                           # m_model = 0: no feedforward
    v, sq = 25.0, 0.0
    for s in grade:
        F = 800 * (25 - v) + drag * 25 + m_model * 9.81 * s           # feedback + feedforward
        v += dt * (F - drag * v - m * 9.81 * s) / m
        sq += (25 - v) ** 2
    return math.sqrt(sq / len(grade))
print(f"RMS speed error: feedback only {drive(0):.3f}, feedforward {drive(1500):.3f}, mass 20 % low {drive(1200):.3f} m/s")

RMS speed error: feedback only 0.577, feedforward 0.000, mass 20 % low 0.115 m/s

With a perfect model the hills vanish; with the mass 20 % wrong, a fifth of the trouble remains, and feedback cleans it up. Feedforward cannot exceed the actuator: an engine at full throttle cannot pay for a steeper hill.

Common mistakes: feedforward without feedback, which leaves every model error uncorrected. Limiting the feedback part before adding the feedforward: the actuator limits the sum.

Practice: A Greenhouse That Reads the Weather

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:

  • Use numpy arrays, masks and @, and print their values plainly.
  • Explain why the FFT costs O(N log N), use np.fft.rfft, and zero-pad to locate a peak.
  • Design windowed-sinc low-pass and band-pass filters, choose a window, and weigh an FIR filter against a biquad.
  • Measure group delay, align a filtered signal, notch out hum, and measure noise reduction in decibels.
  • Find a loop's poles with gain and delay, tabulate a root locus, and bisect for the critical gain.
  • Predict settling, swing and overshoot from a pole's radius and angle.
  • Discretise with a zero-order hold, spot a sample period that is too long, and get a loop's poles from its state matrix.
  • Design a cascade with a fast inner loop, and add feedforward of a measured disturbance.