PID Control, Tuning & Anti-Windup

PID is the default answer to feedback control and has been for eighty years: Åström puts it at more than 95% of industrial control loops, most of them actually PI. The law itself is three lines of arithmetic. Everything that makes a PID controller work — where the derivative is taken, how it is filtered, what the integrator does when the actuator is pinned at its limit, how fast you sample — is not in that equation, and a faithful implementation of the textbook formula is a bad controller. This entry is mostly about the parts that are not in the formula.

The law, and what each term buys #

The non-interacting ("standard", "ISA") form:

$$u(t) = K\left( e(t) + \frac{1}{T_i}\int_0^t e(\tau)\,d\tau + T_d\,\frac{de(t)}{dt} \right), \qquad e = r - y$$

with $r$ the setpoint, $y$ the measurement, $K$ the proportional gain, $T_i$ the integral time and $T_d$ the derivative time. The equivalent parallel form, which is what you actually implement because the parameters enter linearly and each term can be switched off with a finite value:

$$u(t) = k_p\,e(t) + k_i\!\int_0^t e(\tau)\,d\tau + k_d\,\frac{de(t)}{dt}, \qquad k_p = K,\ \ k_i = \frac{K}{T_i},\ \ k_d = K T_d$$

A third parameterization, the interacting (series) form $K_P(1 + 1/sT_i^P)(1 + sT_d^P)$, is common in commercial controllers. It is not the same controller: an interacting equivalent of a non-interacting one exists only when $T_i \ge 4T_d$. If you copy tuning constants between two vendors' boxes without checking which form they implement, you have changed the controller.

Derivative on the measurement, and always filtered #

Two independent problems, both fatal, both fixed at the same place in the code.

Derivative kick. A step in the setpoint produces a step in $e$, hence an impulse in $de/dt$, hence a spike in $u$ that is limited only by the actuator. Nobody wants that. The fix is setpoint weighting: let the P and D terms act on weighted fractions of the reference,

$$u = K\left(b\,r - y + \frac{1}{T_i}\int_0^t e\,d\tau + T_d\left(c\,\frac{dr}{dt} - \frac{dy}{dt}\right)\right), \qquad 0 \le b,c \le 1$$

The integral term must stay on the true error, or steady-state accuracy is gone. In practice $c = 0$ always — that is what "derivative on measurement" means, and it removes the kick entirely. $b$ is the useful knob: it is a second degree of freedom that changes setpoint response without touching the load-disturbance and noise response. Overshoot to a setpoint step is smallest at $b = 0$ (the I-PD controller) and grows with $b$. Tune $K$, $T_i$, $T_d$ for disturbances first; then set $b$ for the reference response.

Noise amplification. Differentiating $\sin t + a_n\sin\omega_n t$ turns a signal-to-noise ratio of $1/a_n$ into $\omega_n / a_n$ — unbounded as $\omega_n$ grows. So the ideal derivative $sT_d$ is never implemented. Replace it with

$$D(s) = \frac{s K T_d}{1 + sT_d/N}\,Y(s)$$

an ideal derivative cascaded with a first-order filter of time constant $T_f = T_d/N$. High-frequency gain of the D-term is now capped at $KN$ instead of infinite. Typical $N$ is 5–20 (Åström & Murray say 5–20, Åström's Ch. 6 says 8–20); $N = 10$ is the usual default.

Discretizing it: positional and incremental #

Two forms, both used in production, with genuinely different operational properties.

Positional (direct) form — the controller outputs the absolute command. Integral by forward rectangular rule, derivative by backward difference:

$$P(t_k) = k_p\big(b\,r(t_k) - y(t_k)\big)$$

$$I(t_{k+1}) = I(t_k) + k_i h\, e(t_k)$$

$$D(t_k) = \frac{T_f}{T_f + h}\,D(t_{k-1}) - \frac{k_d}{T_f + h}\big(y(t_k) - y(t_{k-1})\big)$$

The backward difference is not an aesthetic choice: it makes the pole coefficient $T_f/(T_f+h)$ land in $[0,1)$ for every $h > 0$, so the difference equation is unconditionally stable. A forward difference there can be unstable under discretization for large $h$.

Incremental (velocity) form — the controller outputs $\Delta u$ and the actuator (or a trailing integrator) accumulates it: $u(t) = u(t-h) + \Delta u(t)$, with

$$\Delta u_p = k_p\,\Delta e, \qquad \Delta u_i = k_i h\, e, \qquad \Delta u_d = k_d\left(\frac{\Delta e(t)}{h} - \frac{\Delta e(t-h)}{h}\right)$$

The incremental form's inability to express a P or PD controller is the real constraint. A velocity-form proportional controller reaches stationarity at any error, because a differentiator's output is zero for any constant input. If your loop needs no integral action, use the positional form with an explicit bias term $u_0$ (the control signal corresponding to zero error) — otherwise the steady-state error is zero only at $u = 0$.

Integral windup #

When the actuator saturates, the feedback loop is open: the plant input no longer responds to the controller output. The error stays non-zero, and the integrator keeps integrating it. The integral state grows without bound — it "winds up" — and the command stays past the saturation limit long after the error has changed sign. Recovery requires an equally long period of opposite-sign error just to unwind the state. The visible symptom is a large overshoot followed by a slow, bang-bang-looking oscillation as the command flips between the limits.

Any controller with integral action can wind up. Setpoint changes are the classic trigger, but so are load disturbances, equipment faults, feedforward signals added downstream of the controller, and a controller left running while a selector gave the actuator to somebody else.

flowchart LR
  R["Setpoint r"] --> ES((sum))
  Y["Measurement y"] --> ES
  ES --> I["Integrator ki/s"]
  R --> P["P term, weight b"]
  Y --> P
  Y --> D["Filtered derivative on y"]
  P --> S((sum))
  I --> S
  D --> S
  S --> V["v: nominal command"]
  V --> SAT["Actuator limits"]
  SAT --> U["u: applied command"]
  U --> AW["es = u - v, gain 1/Tt"]
  V --> AW
  AW --> I

Back-calculation (tracking) #

Feed the saturation error $e_s = u - v$ (applied minus nominal command) back into the integrator through a gain $1/T_t$. When nothing saturates, $e_s = 0$ and the extra path does nothing at all. When the actuator saturates, the integrator input becomes

$$\frac{1}{T_t}e_s + \frac{K}{T_i}e = 0 \implies e_s = -\frac{K T_t}{T_i}e \implies v = u_{\text{lim}} + \frac{K T_t}{T_i}e$$

so the nominal command settles just outside the limit rather than far past it, and $u$ responds the instant the error changes sign. $T_t$ is the tracking time constant: it sets how fast the integrator is reset.

The other two schemes #

• Clamping (control-signal clamping). Saturate the output and store the saturated value as the state. In the incremental form this is the whole implementation — it is what you get for free. It is also exactly back-calculation with $T_t = h$.
• Integrator clamping / conditional integration. Freeze the integrator ($\Delta u_i = 0$) whenever the output is saturated. The refined version is better and worth the few extra lines: (1) integrate to saturation — if only the increment pushes the output over the limit, reduce the increment so it lands exactly on the limit rather than dropping it; (2) allow integration in the right direction — while saturated high, block $\Delta u_i > 0$ but permit $\Delta u_i < 0$. Applying this clamp before back-calculation makes the saturation error decay exactly at $T_t$ instead of at a rate shifted by the increment.

A correct discrete PID #

import math


class PID:
    """Positional PID: setpoint weighting, filtered derivative on the
    measurement, back-calculation anti-windup. Fixed sample period h."""

    def __init__(self, kp, ki, kd, h, u_min, u_max, N=10.0, b=1.0):
        if ki <= 0.0:
            raise ValueError("no integral action; use a P/PD form with a bias term")
        if not u_min < u_max:
            raise ValueError("actuator limits must bracket a non-empty range")

        Ti = kp / ki
        Td = kd / kp
        Tf = Td / N
        Tt = math.sqrt(Ti * Td) if Td > 0.0 else Ti

        if Td > 0.0 and Tf < 2.0 * h:
            raise ValueError(f"sample period {h} too coarse for filter Tf={Tf}")

        self.kp, self.kd, self.b = kp, kd, b
        self.u_min, self.u_max = u_min, u_max
        self.bi = ki * h            # integral gain per step
        self.ad = Tf / (Tf + h)     # backward difference: 0 <= ad < 1 for all h
        self.bd = kd / (Tf + h)
        self.br = h / Tt            # anti-windup reset gain
        self.integral = 0.0
        self.deriv = 0.0
        self.y_prev = None

    def step(self, r, y):
        if self.y_prev is None:
            self.y_prev = y

        p = self.kp * (self.b * r - y)
        self.deriv = self.ad * self.deriv - self.bd * (y - self.y_prev)

        v = p + self.integral + self.deriv                  # nominal command
        u = min(max(v, self.u_min), self.u_max)             # actuator model

        self.integral += self.bi * (r - y) + self.br * (u - v)
        self.y_prev = y
        return u

Note the ordering: the output is computed and applied before the states are updated, which is the sequence that minimizes latency (read input → compute → write output → update state). Both the derivative and the anti-windup path read $y$ and $v$ from this cycle only. The Tf < 2h check is the one non-obvious guard — see below.

Choosing the sample rate #

Sampling is pure loss: it buys nothing and adds dead time. With the ideal loop ordering above, the average measurement delay is $h/2$ and the computation is usually negligible, so the loop inherits roughly $h/2$ of extra delay. Most controllers and instrument systems do not order their calculations that way and end up with several sample periods of delay.

A delay of $h/2$ costs phase margin $\approx \omega_c h / 2$ radians at crossover. Working that arithmetic against the range of $\omega_c h$ commonly quoted in the sampled-data literature:

$$\omega_c h \in [0.15,\ 0.5] \implies \text{phase-margin loss of } 4.3^\circ \text{ to } 14.3^\circ$$

The $[0.15,\ 0.5]$ range itself is convention rather than a result — treat the phase-margin cost as the real constraint and the range as a starting point. Expressed against the ultimate period $T_u$ (with $\omega_u = 2\pi/T_u$, and crossover typically near $\omega_u$), that is roughly 13 to 42 samples per ultimate period — a useful number because $T_u$ is exactly what the tuning experiment below measures. Faster is not free either: the derivative's noise amplification is $\Delta n / h$, inversely proportional to $h$, so shrinking $h$ makes the measurement filter more necessary, not less.

Three further constraints:

• Anti-aliasing is mandatory, in analog, before the ADC. Aliasing does not merely lose high-frequency content; it folds it down into the control bandwidth where the loop will faithfully act on it. A second-order Butterworth is the common choice. Åström's sizing example: to attenuate by 16 at $\omega_s/2$ with a second-order filter, solve $(\omega_s / 2\omega_b)^2 = 16$, giving $\omega_b = \omega_s/8$. The prefilter's own dynamics are significant and belong in the plant model for design purposes.
• The derivative filter must be resolvable at $h$. In $D_k = \frac{T_f}{T_f+h}D_{k-1} - \frac{k_d}{T_f+h}\Delta y$, if $h \gg T_f$ then $T_f/(T_f+h) \to 0$ and the recursion degenerates into a raw first difference with gain $k_d/h$ — the filter you specified does not exist. Keep $T_f$ at least a small multiple of $h$; that is what the guard in the code above enforces.
• Jitter. Fixed $h$ is the normal assumption and small jitter is normally ignored. Where jitter is large or deliberate (event-based control), compute $\Delta t$ as the actual elapsed time — and remember that a changed $\Delta t$ also rescales any rate limits expressed per sample.

Tuning #

Ziegler-Nichols, with the real numbers #

Two experiments, published in 1942, and still the reference point everyone argues against.

Frequency-response (ultimate-gain) method. Set $T_i = \infty$, $T_d = 0$, raise $K$ until the loop oscillates with constant amplitude. That gain is the ultimate gain $K_u$ and the oscillation period is the ultimate period $T_u$; by Nyquist, the experiment has located the point where the process has exactly 180° of phase lag, at $\omega_{180} = 2\pi/T_u$.

Step-response (bump test) method. Drive the process open-loop with a step, draw the tangent at the point of steepest slope, and read off its intercepts: $L$ (apparent dead time, the time-axis intercept) and $a$ (the value-axis intercept for a unit step).

Softer variants of the frequency-response rules circulate widely (a "no overshoot" variant at $K = 0.20K_u$, $T_i = 0.50T_u$, $T_d = 0.33T_u$; Pessen's integral rule at $0.7K_u$, $0.4T_u$, $0.15T_u$). They are the same experiment with different constants and the same fundamental limitation: too little information about the process.

Relay auto-tuning #

Running a proportional loop up to sustained oscillation is unpleasant on real equipment — you are deliberately driving the plant to its stability boundary and hoping the amplitude stays bounded. Åström and Hägglund's 1984 alternative replaces the proportional controller with a relay of amplitude $d$. The relay-plus-plant loop settles into a limit cycle automatically, at the frequency where the plant has 180° of lag. Describing-function analysis of the ideal relay (keeping only the fundamental of the square wave) gives

$$K_u = \frac{4d}{\pi a}, \qquad T_u = \frac{2\pi}{\omega_u}$$

where $a$ is the measured amplitude of the output oscillation. The experiment is short, self-limiting (amplitude is controlled by adjusting $d$), and automatable — push a button, the controller runs the relay test, computes $K_u$ and $T_u$, and installs a PID. This is what "autotune" means on most commercial controllers.

The catch is the describing-function approximation: it discards all harmonics above the fundamental, and the resulting error in $K_u$ is typically 5–20% for common process transfer functions. The frequency estimate is considerably better than the amplitude estimate. A relay with hysteresis $\epsilon$ has describing function $\frac{4d}{\pi a}\left(\sqrt{1 - \epsilon^2/a^2} - j\,\epsilon/a\right)$, which lets a sweep of $\epsilon$ identify points at known phase shifts rather than only the 180° point.

Manual loop shaping #

The order that actually works, and matches how the terms interact:

1. P alone. $T_i = \infty$, $T_d = 0$. Raise $K$ until the response is as fast as you can accept, then back off. Accept the steady-state offset for now.
2. Add I. Reduce $T_i$ from a large value until steady-state error is removed at an acceptable settling time. Oscillation tendency rises as $T_i$ falls; when it starts oscillating you have gone too far. Most loops stop here — most industrial loops are PI.
3. Add D only if it helps. Derivative action is worth adding for lag-dominated processes (normalized delay $\tau_n = L/(L+T)$ near 0) and is close to worthless for delay-dominated ones ($\tau_n$ near 1). Increase $T_d$ to damp the oscillation; damping improves, then gets worse again once $T_d$ exceeds roughly $T_u/6$.
4. Then set $b$. Only after the disturbance response is right, reduce $b$ from 1 toward 0 to trim setpoint overshoot. It cannot destabilize the loop — it is outside the feedback path.

PID in software systems #

Feedback controllers show up in infrastructure code far more often than they are labelled as such, and the ones that are labelled are worth reading.

PIE (RFC 8033) is an active queue management scheme that is explicitly a PI controller on queueing latency. Every T_UPDATE it computes an increment to the packet drop probability:

p = alpha * (current_qdelay - QDELAY_REF) +
        beta * (current_qdelay - PIE->qdelay_old_);
PIE->drop_prob_ += p;

with QDELAY_REF = 15 ms (the setpoint), T_UPDATE = 15 ms (the sample period), alpha = 1/8 Hz and beta = 1 + 1/4 Hz. This is exactly the incremental form from above: the accumulated state is the output (drop_prob_) rather than a separate integrator, the alpha term is the integral action on the error, and the beta term — the trend, $q_k - q_{k-1} = e_k - e_{k-1}$ — is the proportional action. The RFC's own scaling rule confirms the reading: "if we cut T_UPDATE in half, we should also cut alpha by half", which is what $k_i h$ does. The bound $0 \le p \le 1$ is clamping anti-windup, and because clamping alone leaves residual state after a congestion episode ends, PIE adds an explicit exponential decay (×0.98 per update, a time constant of about $50 \times$ T_UPDATE) to unwind it. RFC 8033 also documents the deliberate detuning: alpha and beta are divided by up to 2048 at low drop probabilities, a gain schedule on the operating point.

Kubernetes HPA is the instructive counter-example — a loop that deliberately refuses integral action:

$$\text{desiredReplicas} = \left\lceil \text{currentReplicas} \times \frac{\text{currentMetricValue}}{\text{desiredMetricValue}} \right\rceil$$

That is proportional control with a ratio instead of a difference, and nothing else. There is no integrator because the plant — pods starting, warming up, joining a load balancer — has long, variable, poorly-known dead time, and integral action against unmodelled dead time is how you get sustained oscillation. Instead the loop is stabilized the way delay-dominated loops usually are: a tolerance deadband (default 0.1 — no action while the ratio is within 10% of 1.0) and a stabilization window on scale-down (default 300 s, taking the maximum recommendation over the window, a rolling max that refuses to act on transients). Deadband plus rate limiting instead of derivative damping, and no integrator at all.