M3-T7 — Digital Control: Sampling, Delay and Anti-Windup

Learning objectives

After this lesson, you should be able to:

• Discretize a compensator with the bilinear (Tustin) transform and implement it as a difference equation.
• Compute the phase lag of sampling and choose a sample rate for a loop.
• Explain integrator windup and cure it with back-calculation.
• Design a notch filter and state its cost in phase.

The intuition

The controller of M3-T5: Loop Shaping: Lead, Lag and PID as Compensators is a transfer function. Nothing in it says how often it is computed, or what happens when its output asks for more than the actuator can give. Both questions have to be answered before it can fly.

From a transfer function to a program

A flight computer reads the sensors at fixed instants, computes a new output and writes it to the actuator, where it stays until the next instant. To turn a compensator $C(s)$ into a program, replace each $s$ by an approximation of a derivative from samples. The Tustin (bilinear) transform $s \to (2/T)(z - 1)/(z + 1)$ is the standard choice: it treats an integrator as the trapezoid rule, keeps stability, and turns the compensator into a difference equation, a weighted sum of recent inputs and outputs, which is a few lines of code.

Sampling costs half a period

Because the output is held for one period $T$, it is on average half a period late. At the crossover frequency $\omega_{gc}$ that is a phase lag of $\omega_{gc}T/2$. A designer treats it like any other lag: it comes out of the phase margin. A guide is to sample 10 to 20 times faster than the crossover, so that the lag is a few degrees. The example below measures it: at 100 Hz the loss is 1.7° of a 55.5° margin, at 10 Hz 17°, and at 3.3 Hz it is 52° and the loop is on the edge.

Windup

The actuator has limits. When the controller asks for more than the gimbal can give, the output is clipped, the pitch error does not shrink as the controller expects, and the integral term, which adds up the error, keeps growing. The integrator has stored a large command that cannot be delivered. When the vehicle finally nears its target, the stored command is still there and drives it past, and the overshoot is far larger than the linear design promised. This is windup, and a loop with an integral term and a saturating actuator always has it.

The cure is to tell the integrator what the actuator did. Back-calculation feeds the difference between the clipped output and the commanded output back into the integrator, so that while the actuator saturates the integrator is pulled back to a value that just reaches the limit. It costs one extra gain and one line of code.

Filters for what you should not follow

Some frequencies should be kept out of the loop: the bending mode of a slender rocket, a motor vibration seen by the gyroscope. A notch filter is a narrow dip in the gain at one frequency. It works, at the price of a phase lag below the notch, which comes out of the same margin as the sampling lag.

Sampling, saturation and structural modes each take something from the loop: phase, linearity, or bandwidth. The recipes are short, and the price of each is a number you can compute before flying.

Definitions

Sample period $T$ and rate $f_s = 1/T$

The time between updates of the controller, and its inverse. The sample frequency in rad/s is $\omega_s = 2\pi/T$.

Zero-order hold (ZOH)

The conversion of a computed value into a continuous signal by holding it constant until the next sample.

Difference equation

The program form of a discrete compensator: $y_k = \sum_i b_i u_{k-i} - \sum_{i \ge 1} a_i y_{k-i}$.

Windup

The growth of the integral state while the actuator is saturated, so that the stored command exceeds what the actuator can deliver.

Back-calculation time constant $T_t$

The speed with which the anti-windup feedback pulls the integrator back: smaller means faster.

Notch filter

A filter $(s^2 + 2\zeta_z\omega_0 s + \omega_0^2)/(s^2 + 2\zeta_p\omega_0 s + \omega_0^2)$ with a dip of depth $\zeta_z/\zeta_p$ at $\omega_0$.

A zero-order hold followed by a continuous plant is equivalent to the plant preceded by the transfer function:

The Tustin substitution and the frequency warping that it causes are:

A PID controller with back-calculation updates its integral state $I$ by:

Derivation: the lag of sampling, Tustin, and the equilibrium of anti-windup

Each of the three recipes rests on a short calculation.

The delay of a zero-order hold

1. $$u(t) = u_k \ \ \text{for}\ \ kT \le t < (k+1)T$$

   A zero-order hold writes each computed value as a constant for one period.

2. $$H(s) = \int_0^T e^{-st}\,dt = \frac{1 - e^{-sT}}{s}$$

   The response of the hold to a unit impulse of the sampled sequence is a rectangle of length $T$, whose Laplace transform is this.

3. $$1 - e^{-j\omega T} = e^{-j\omega T/2}\bigl(e^{j\omega T/2} - e^{-j\omega T/2}\bigr) = 2j\,e^{-j\omega T/2}\sin\tfrac{\omega T}{2}$$

   Factor out half the exponent. The bracket is twice $j$ times a sine.

4. $$H(j\omega) = T\,\frac{\sin(\omega T/2)}{\omega T/2}\,e^{-j\omega T/2}, \qquad \angle H = -\frac{\omega T}{2}$$

   Divide by $j\omega$. The sine ratio is real and close to 1 for small $\omega T$, so the phase is exactly $-\omega T/2$: a delay of half a period at every frequency.

Why Tustin works, and where it warps

1. $$z = e^{sT} \approx \frac{1 + sT/2}{1 - sT/2} \;\Rightarrow\; s \approx \frac{2}{T}\,\frac{z - 1}{z + 1}$$

   Approximate the exponential by the ratio of the first-order terms of its numerator and denominator (the Padé approximant). Solving for $s$ gives the bilinear substitution, which is the trapezoid rule for integration.

2. $$z = e^{j\omega T} \;\Rightarrow\; \frac{2}{T}\,\frac{e^{j\omega T} - 1}{e^{j\omega T} + 1} = \frac{2}{T}\,j\tan\frac{\omega T}{2}$$

   On the unit circle the substitution gives a purely imaginary $s$, so stable stays stable. The frequency seen by the continuous design is $\Omega = (2/T)\tan(\omega T/2)$.

3. $$\Omega \approx \omega\left(1 + \frac{(\omega T)^2}{12}\right)$$

   For small $\omega T$ the warping is second order: at $\omega T = 0.4$ it is 1.3%, which is why the Tustin lead of the example changes its gain by under 1% at 100 Hz.

Where anti-windup settles

1. $$\dot{I} = 0 \;\Rightarrow\; K_i\,e = \frac{u - u_{\text{sat}}}{T_t}$$

   In a steady saturation the integrator stops moving when its two inputs balance.

2. $$u - u_{\text{sat}} = T_t\,K_i\,e$$

   The command can exceed the limit by only $T_t K_i e$, which is small for a small $T_t$. Without the feedback $T_t \to \infty$ and the excess is limited by nothing.

3. $$\Delta I \approx K_i\,e\,t_{\text{sat}}$$

   Without anti-windup the integrator grows linearly during the time $t_{\text{sat}}$ of saturation. This is the stored command that must later be unwound.

Sampling delays the loop by half a period, saturation breaks the integrator, and a filter takes phase for the band it removes. Each has a number: ωT/2, T_t K_i e, and the notch phase lag.

A worked example

Run the pitch loop on a computer with the real servo and see what each of the three costs.

1. Tustin. Discretize the lead of M3-T5: Loop Shaping: Lead, Lag and PID as Compensators (zero 1.234, pole 12.966, gain 0.2606) at three periods and compare its response at 4 rad/s. Continuous: gain 0.8447, phase 55.71°. At 10 ms the Tustin version has gain 0.8448, at 50 ms 0.8471 and at 100 ms 0.8542, all with the same phase to two decimals: the warping is small because $\omega T$ is at most 0.4.
2. Sampling lag. Simulate the PD loop ($K_p = 0.84, K_d = 0.285$) with the plant and servo continuous and the controller updated every $T$. At 100 Hz the overshoot of a 0.1 rad step is 8.0%, at 20 Hz 8.3%, at 10 Hz 9.4%, at 5 Hz 19.7% and at 3.3 Hz 49.3%. The predicted loss of phase $\omega_{gc}T/2$ is 1.7°, 8.7°, 17.4°, 34.8° and 52.2° of the 55.5° available: the loop is almost out of margin at 3.3 Hz, and would be unstable at about 3.1 Hz.
3. Windup. Give the pitch loop an integral term (poles at −1.5 and −2.8 ± 2.86j, $K_p = 1.27, K_i = 1.22, K_d = 0.362$) and the saturating servo, update it at 100 Hz, and command steps of 0.1, 0.4 and 0.8 rad. The gimbal saturates for the larger two. With a plain integrator the overshoot is 31.5%, 41.6% and 101.5%, settling in 2.5, 2.9 and 5.0 s. With back-calculation ($T_t = 0.2$ s) it is 31.5%, 7.1% and 0%, settling in 2.5, 1.7 and 2.0 s. The first, unsaturated, case is identical.
4. Notch. A notch at 30 rad/s with zero damping 0.05 and pole damping 0.7 has a gain of 0.071 (−22.9 dB) at 30 rad/s, exactly $\zeta_z/\zeta_p$. Its price: a phase of −15.1° at 6 rad/s, the crossover of the PD design, and −31° at 12 rad/s. Nearly a third of the 55.5° margin is spent on the notch, which is why a structural mode is first moved, if possible, and filtered only if not.

import cmath, math
import gnckit as g

b, lam, wn, zeta = 19.62, 0.7, 4.0, 0.7
Kd, Kp = 2 * zeta * wn / b, (wn ** 2 + lam ** 2) / b

# 1. a lead compensator, continuous against its Tustin discretization
z, p, K = 1.234, 12.966, 0.2606
num, den = [K / z, K], [1 / p, 1]
print('lead compensator at 4 rad/s, continuous and discretized (Tustin):')
c = K * (1 + 1j * 4 / z) / (1 + 1j * 4 / p)
print(f'  continuous           : gain {abs(c):.4f}, phase {math.degrees(cmath.phase(c)):6.2f} deg')
for T in (0.01, 0.05, 0.1):
    bz, az = g.tustin(num, den, T)
    zz = cmath.exp(1j * 4 * T)
    d = sum(bk * zz ** (-k) for k, bk in enumerate(bz)) / sum(ak * zz ** (-k) for k, ak in enumerate(az))
    print(f'  Tustin, T = {T:4.2f} s : gain {abs(d):.4f}, phase {math.degrees(cmath.phase(d)):6.2f} deg')

# 2. the sampled PD loop: simulate the plant continuously and update the controller every T seconds
def run(T_s, step=0.1, t_end=10.0, dt=0.001):
    act, th, q, cmd, peak = g.Actuator(40.0, 0.15, 3.0), 0.0, 0.0, 0.0, 0.0
    sub = int(round(T_s / dt))
    for i in range(int(t_end / dt)):
        if i % sub == 0:
            cmd = Kp * (step - th) - Kd * q
        y = act.step(cmd, dt)
        th, q = g.rk4_step(lambda t, s: [s[1], lam ** 2 * s[0] + b * y], 0.0, [th, q], dt)
        peak = max(peak, th)
    return peak / step - 1
print('sampled PD loop, step of 0.1 rad; the lag of sampling at the 6.07 rad/s crossover is wgc T / 2:')
for T_s in (0.01, 0.05, 0.1, 0.2, 0.3):
    print(f'  T = {T_s:4.2f} s ({1 / T_s:5.1f} Hz): overshoot {100 * run(T_s):5.1f}%, phase lost {math.degrees(6.07 * T_s / 2):5.1f} deg of 55.5')

# 3. integrator windup: PID on the pitch with the saturating, rate-limited gimbal, updated at 100 Hz
a_, wn_ = 1.5, 4.0
c2, c1, c0 = a_ + 2 * zeta * wn_, wn_ ** 2 + 2 * zeta * wn_ * a_, a_ * wn_ ** 2
Kd2, Kp2, Ki2 = c2 / b, (c1 + lam ** 2) / b, c0 / b
def pid(step, anti_windup, Ts=0.01, Tt=0.2, t_end=12.0, dt=0.001):
    act, th, q, I, cmd, peak, trace = g.Actuator(40.0, 0.15, 3.0), 0.0, 0.0, 0.0, 0.0, 0.0, []
    for i in range(int(t_end / dt)):
        if i % int(round(Ts / dt)) == 0:
            e = step - th
            u = Kp2 * e + I - Kd2 * q
            usat = max(-0.15, min(0.15, u))
            I += Ts * (Ki2 * e + ((usat - u) / Tt if anti_windup else 0.0))
            cmd = u
        y = act.step(cmd, dt)
        th, q = g.rk4_step(lambda t, s: [s[1], lam ** 2 * s[0] + b * y], 0.0, [th, q], dt)
        peak = max(peak, th)
        trace.append(th)
    settle = max([i * dt for i, v in enumerate(trace) if abs(v - step) > 0.02 * step] or [0.0])
    return 100 * (peak / step - 1), settle
print('PID with a saturating gimbal (limit 0.15 rad): overshoot and 2% settling time')
for step in (0.1, 0.4, 0.8):
    o1, s1 = pid(step, False)
    o2, s2 = pid(step, True)
    print(f'  step {step:.1f} rad: plain integrator {o1:6.1f}% / {s1:4.1f} s, with back-calculation {o2:6.1f}% / {s2:4.1f} s')

# 4. a notch filter for a structural mode at 30 rad/s: deep and narrow, at the price of some phase lag below it
w0, zz, zp = 30.0, 0.05, 0.7
notch = g.TF([1, 2 * zz * w0, w0 ** 2], [1, 2 * zp * w0, w0 ** 2])
print('notch at 30 rad/s (zero damping 0.05, pole damping 0.7):')
for w in (6.0, 12.0, 30.0, 100.0):
    v = notch.freq(w)
    print(f'  w = {w:5.1f}: gain {abs(v):.3f} ({20 * math.log10(abs(v)):6.1f} dB), phase {math.degrees(cmath.phase(v)):6.1f} deg')

lead compensator at 4 rad/s, continuous and discretized (Tustin):
  continuous           : gain 0.8447, phase  55.71 deg
  Tustin, T = 0.01 s : gain 0.8448, phase  55.71 deg
  Tustin, T = 0.05 s : gain 0.8471, phase  55.71 deg
  Tustin, T = 0.10 s : gain 0.8542, phase  55.71 deg
sampled PD loop, step of 0.1 rad; the lag of sampling at the 6.07 rad/s crossover is wgc T / 2:
  T = 0.01 s (100.0 Hz): overshoot   8.0%, phase lost   1.7 deg of 55.5
  T = 0.05 s ( 20.0 Hz): overshoot   8.3%, phase lost   8.7 deg of 55.5
  T = 0.10 s ( 10.0 Hz): overshoot   9.4%, phase lost  17.4 deg of 55.5
  T = 0.20 s (  5.0 Hz): overshoot  19.7%, phase lost  34.8 deg of 55.5
  T = 0.30 s (  3.3 Hz): overshoot  49.3%, phase lost  52.2 deg of 55.5
PID with a saturating gimbal (limit 0.15 rad): overshoot and 2% settling time
  step 0.1 rad: plain integrator   31.5% /  2.5 s, with back-calculation   31.5% /  2.5 s
  step 0.4 rad: plain integrator   41.6% /  2.9 s, with back-calculation    7.1% /  1.7 s
  step 0.8 rad: plain integrator  101.5% /  5.0 s, with back-calculation   -0.0% /  2.0 s
notch at 30 rad/s (zero damping 0.05, pole damping 0.7):
  w =   6.0: gain 0.960 (  -0.4 dB), phase  -15.1 deg
  w =  12.0: gain 0.833 (  -1.6 dB), phase  -31.0 deg
  w =  30.0: gain 0.071 ( -22.9 dB), phase    0.0 deg
  w = 100.0: gain 0.908 (  -0.8 dB), phase   22.9 deg

In robotics

These three effects are present in every flight computer, and the recipes are standard.

• Loop rates. Multirotor flight controllers run their inner loops at hundreds of hertz to several kilohertz, for the reason of this lesson: the phase lost to sampling must be a few degrees at a crossover of tens of rad/s. A slower outer loop can run at tens of hertz. The rate hierarchy of M1-T2: Nested Loops and Separation of Timescales is a sampling hierarchy too.
• Delays beyond the hold. The half-period of the hold is only the beginning: the time to read the sensor, compute and write the output is a further delay, often larger. Jitter, a delay that varies, is worse, because a loop designed for the nominal delay meets longer ones. Timing analysis and deadline monitoring are the subject of UAV M7-T7: Jitter and deadline monitoring.
• Anti-windup in practice. Flight stacks limit the integrator or freeze it when the output is saturated, which is a cruder form of the same remedy. The failure it prevents is familiar: after a large manoeuvre or a motor limit, the vehicle swings past its target and takes seconds to settle.
• Notch filters. Flight controllers filter the gyroscope with notches placed at the motor vibration frequencies, and some move them as the throttle changes. Launch vehicles notch the bending modes of the structure. Each notch spends phase below its frequency, so the filter is the last resort after stiffening the structure or changing the sensor location.
• Verification. A discretized controller must be tested at its real rate with the real actuator model, in simulation first. Software-in-the-loop setups (UAV M6-T14) run the production code at its rate against a model of the vehicle, which is where these effects first show.

This closes the classical part of the course. The tools of the module are the loop gain, the margins, the lead and the limits. The state-space methods of module 4 design several loops at once, and use an observer in place of the rate gyroscope. They inherit every limit of this module: delays, saturation and the unstable pole do not go away because the algebra changes.

Common misconceptions

Key takeaways

1. A flight computer implements a compensator by a difference equation. The Tustin transform $s \to (2/T)(z - 1)/(z + 1)$ does it and keeps stability, with a small frequency warping.
2. A zero-order hold delays by half a period: a phase lag $\omega T/2$, to be paid from the phase margin. Sample 10 to 20 times faster than the crossover, and measure the rest.
3. Hopper's PD loop loses 1.7° at 100 Hz, 17° at 10 Hz and 52° at 3.3 Hz; its overshoot rises from 8% to 49%.
4. Windup is the growth of the integrator while the actuator is saturated; back-calculation, $\dot I = K_i e + (u_{\text{sat}} - u)/T_t$, cures it and leaves the unsaturated loop unchanged.
5. For large steps the plain integrator overshoots by 41.6% and 101.5%; back-calculation by 7.1% and 0%.
6. A notch removes one frequency and lags the phase below it: 15° at the crossover for a notch at 30 rad/s.

References and further reading

Each source is marked free or paid. Free to read does not mean free to reuse: the course resources page lists every source with its licence.

• K. J. Åström and R. M. Murray, Feedback Systems, 2nd edition, Princeton University Press, 2021 , the sections on implementation, saturation and integrator windup. Free to read as the authors' PDF; the book is copyrighted.
• B. Douglas, “Discrete Control” (video series) . Free to watch.
• MATLAB Tech Talks, “Control Systems in Practice” (video series) , with videos on anti-windup. Free to watch.
• G. F. Franklin, J. D. Powell and M. L. Workman, Digital Control of Dynamic Systems, 3rd edition, Ellis-Kagle Press, 1998, the chapters on design by emulation and on sample-rate selection. The standard text. Paid book.

Knowledge check

Three questions. Write an answer, then reveal the explanation. Your answers stay in your browser.

What’s next

Related content

Nothing below is required to finish this lesson.

• Before
  M3-T6: Limits of Performance: Delay, Right-Half-Plane Zeros and the Bode Integral

  Explain why delay and non-minimum-phase zeros bound the achievable bandwidth, and why pushing sensitivity down at one frequency raises it at another.

• Lab
  M3-L2: When the Gimbal Saturates

  A large pitch command saturates the gimbal and the integrator winds up. Implement the controller at 50 Hz with anti-windup and show the overshoot removed.

Other courses on this site treat neighbouring topics from their own angle:

• Lab: Tune a PID Controller
• Lab: Fix the Oscillating Robot
• UAV M5-T6: PID controllers
• UAV M7-T6: End-to-end latency

This course leaves the following to the pages that already teach them:

• What P, I and D do, and tuning them by hand on a step response
• PID cascades for a multirotor