Kalman Filtering & Model Predictive Control

Estimation and optimisation at the top of the control pyramid. Learn the Kalman filter (predict/update, gain, covariance, Q and R tuning), sensor fusion and the extended Kalman filter, and model predictive control (prediction horizon, cost function, constraints, receding horizon, offset-free MPC) — with a running Kalman filter and a constrained MPC in plain Python.

Measurements are noisy and models are imperfect. The Kalman filter combines both optimally, weighting each by its uncertainty. It navigates aircraft, phones and rockets, and estimates unmeasured process variables.

Model predictive control (MPC) uses a model to predict the future and optimises a sequence of moves that respect limits (valve ranges, rate limits, product specs). It is the standard advanced process control technology in refineries and chemical plants.

An analogy

A Kalman filter is a navigator mixing dead reckoning (the model) with occasional landmark sightings (measurements), trusting whichever is more reliable. MPC is a chess player: it plans several moves ahead within the rules, plays only the first move, then re-plans after seeing the opponent's reply.

Why estimation?

Sensors are noisy, slow or missing (you cannot measure every temperature inside a reactor). Averaging reduces noise but adds lag. The Kalman filter uses a model to predict and then corrects with measurements, giving the minimum-variance estimate for linear systems with Gaussian noise.

Tuning Q and R, and sensor fusion

R reflects sensor noise (from datasheets or data); Q reflects how much the true state can change unpredictably. Large Q/R → trust measurements (fast, noisy); small Q/R → trust the model (smooth, slower). Fusing sensors (GPS + accelerometer, flow + level) is natural: each measurement updates the state with its own R. For nonlinear systems use the extended (EKF) or unscented (UKF) Kalman filter.

MPC basics

An MPC controller at each sample:

Industrial APC

Refinery MPCs (AspenTech DMC3, Honeywell Profit Controller, Yokogawa, Emerson) handle dozens of manipulated variables and constraints, pushing units to their economic optimum: typically 2–6% throughput or energy benefit. Key practical elements: step-test identification of dynamic models, steady-state LP/QP optimiser layer, constraint handling and priority, disturbance models for offset-free control, and sustained monitoring (APC benefits decay without maintenance).

Formulas

Worked example: Kalman gain

Given: P⁻ = 0.04, R = 0.25

Find: Gain.

Answer: 0.138 — mostly trust the prediction

Worked example: Fusing two sensors

Given: σ₁ = 1 °C, σ₂ = 2 °C (independent)

Find: Fused σ.

Answer: σ = 0.894 °C

Worked example: Horizon length

Given: Process settles in 60 min, sample 1 min

Find: Prediction horizon.

Answer: N_p ≈ 60–90 samples

1-D Kalman filter estimating a temperature

import random, math
random.seed(1)
# 1-D Kalman filter estimating a constant temperature from noisy readings
true_T, R = 25.0, 0.5**2          # measurement noise variance (sigma 0.5 C)
Q = 1e-5                          # process noise variance
x, P = 20.0, 10.0                 # initial guess and its variance
errs = []
for k in range(1, 51):
    z = true_T + random.gauss(0, 0.5)
    P = P + Q                     # predict
    Kk = P / (P + R)              # Kalman gain
    x = x + Kk * (z - x)          # update
    P = (1 - Kk) * P
    if k in (1, 5, 20, 50):
        print(f"k={k:2d}: z={z:6.2f}  estimate={x:6.3f}  gain={Kk:.3f}  sigma={math.sqrt(P):.3f}")
# k= 1: z= 25.64  estimate=25.506  gain=0.976  sigma=0.494
# k= 5: z= 24.45  estimate=25.069  gain=0.199  sigma=0.223
# k=20: z= 24.93  estimate=25.000  gain=0.050  sigma=0.112
# k=50: z= 24.35  estimate=25.056  gain=0.021  sigma=0.072

At first the filter trusts measurements (gain ≈ 1) because its own estimate is poor; as confidence grows the gain falls and single noisy readings barely move the estimate. The estimate's σ (0.07 °C) ends up far smaller than the sensor's 0.5 °C.

Constrained MPC with receding horizon

# Receding-horizon MPC for a first-order process x+ = a x + b u,
# horizon 10, input limits |u| <= 1, move-rate limit |du| <= 0.3, by projected gradient.
a, b, N = 0.9, 0.1, 10
q, r = 1.0, 0.01
def cost(x0, us, sp):
    x, J = x0, 0.0
    for u in us:
        x = a * x + b * u
        J += q * (x - sp) ** 2 + r * u * u
    return J
def solve(x0, u_prev, sp):
    us = [u_prev] * N
    for it in range(400):
        base = cost(x0, us, sp)
        grad = []
        for k in range(N):
            d = us[:]; d[k] += 1e-5
            grad.append((cost(x0, d, sp) - base) / 1e-5)
        us = [u - 0.5 * g for u, g in zip(us, grad)]
        prev = u_prev
        for k in range(N):                            # project onto constraints
            lo, hi = max(-1.0, prev - 0.3), min(1.0, prev + 0.3)
            us[k] = min(hi, max(lo, us[k])); prev = us[k]
    return us
x, u, sp = 0.0, 0.0, 0.8
for t in range(40):
    u = solve(x, u, sp)[0]                            # apply only the first move
    x = a * x + b * u
    if t in (0, 1, 2, 3, 5, 10, 20, 39):
        print(f"k={t:2d}: u={u:5.2f}  x={x:5.3f}")
print(f"steady-state input needed: (1-a)*sp/b = {(1-a)*sp/b:.2f} (inside the |u| <= 1 limit)")
# k= 0: u= 0.30  x=0.030
# k= 1: u= 0.60  x=0.087
# k= 2: u= 0.90  x=0.168
# k= 3: u= 1.00  x=0.251
# k= 5: u= 1.00  x=0.394
# k=10: u= 1.00  x=0.642
# k=20: u= 0.79  x=0.792
# k=39: u= 0.79  x=0.792
# steady-state input needed: (1-a)*sp/b = 0.80 (inside the |u| <= 1 limit)

The controller ramps the input at the allowed 0.3 per step, rides the upper limit, then settles — constraints are handled by the optimiser, not by clamping afterwards. Note the small offset (0.792 vs 0.8): penalising u itself trades a little error for effort. Industrial MPC penalises Δu and adds a disturbance model (offset-free MPC). Real MPC uses QP solvers (OSQP, qpOASES) instead of this simple gradient method.

Common mistakes

Applications

Key takeaways

Going further