PID Controller Simulation
Goal: Simulate a PID controller in Python controlling a thermal system. Tune it with Ziegler-Nichols, observe the effect of each gain, and deal with real-world problems (integral windup, derivative noise).
Prerequisites: PID Controller, Open Loop vs Closed Loop, Digital Control
The Plant: A Thermal System
We’re controlling the temperature of a heated box. The heater adds energy; the box loses heat to the environment proportionally to the temperature difference.
Setpoint (desired temp)
↓
[PID] → heater power (0-100%) → [Box] → measured temp
↑ |
└──────────── thermocouple ──────────────┘
Physics: dT/dt = (heater_power / C) - (T - T_ambient) / (R * C)
Where R = thermal resistance, C = thermal capacitance. This is a first-order system with a time constant τ = R * C.
Step 1: Simulate the Plant
import numpy as np
import matplotlib.pyplot as plt
class ThermalPlant:
"""First-order thermal system with delay."""
def __init__(self, R=1.0, C=10.0, T_ambient=20.0, dt=0.1):
self.R = R # thermal resistance (°C/W)
self.C = C # thermal capacitance (J/°C)
self.T_ambient = T_ambient
self.dt = dt
self.T = T_ambient # initial temperature
self.tau = R * C # time constant = 10s
def step(self, heater_power):
"""Advance one time step. heater_power in watts (0-100)."""
heater_power = np.clip(heater_power, 0, 100)
dT = (heater_power / self.C - (self.T - self.T_ambient) / (self.R * self.C)) * self.dt
self.T += dT
return self.TTime constant τ = 10s means the system reaches 63% of its final value in 10 seconds.
Step 2: PID Controller
class PIDController:
def __init__(self, Kp=0, Ki=0, Kd=0, dt=0.1, out_min=0, out_max=100):
self.Kp, self.Ki, self.Kd = Kp, Ki, Kd
self.dt = dt
self.out_min, self.out_max = out_min, out_max
self.integral = 0.0
self.prev_error = 0.0
self.prev_measurement = None
def update(self, setpoint, measurement):
error = setpoint - measurement
# Proportional
P = self.Kp * error
# Integral
self.integral += error * self.dt
I = self.Ki * self.integral
# Derivative on measurement (avoids kick on setpoint change)
if self.prev_measurement is None:
D = 0.0
else:
D = -self.Kd * (measurement - self.prev_measurement) / self.dt
self.prev_measurement = measurement
output = P + I + D
# Anti-windup: clamp and back-calculate
if output > self.out_max:
self.integral -= (output - self.out_max) / self.Ki if self.Ki else 0
output = self.out_max
elif output < self.out_min:
self.integral -= (output - self.out_min) / self.Ki if self.Ki else 0
output = self.out_min
self.prev_error = error
return outputStep 3: Run the Simulation
def simulate(Kp, Ki, Kd, setpoint=60.0, duration=100.0, dt=0.1):
plant = ThermalPlant(dt=dt)
pid = PIDController(Kp=Kp, Ki=Ki, Kd=Kd, dt=dt)
t = np.arange(0, duration, dt)
temps = np.zeros_like(t)
powers = np.zeros_like(t)
for i in range(len(t)):
power = pid.update(setpoint, plant.T)
plant.step(power)
temps[i] = plant.T
powers[i] = power
return t, temps, powers
# P-only
t, temps, powers = simulate(Kp=5.0, Ki=0, Kd=0)
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(t, temps, label='Temperature')
plt.axhline(y=60, color='r', linestyle='--', label='Setpoint')
plt.xlabel('Time (s)'); plt.ylabel('°C'); plt.legend(); plt.grid()
plt.title('P-only: steady-state error')
plt.subplot(1, 2, 2)
plt.plot(t, powers)
plt.xlabel('Time (s)'); plt.ylabel('Heater Power (W)'); plt.grid()
plt.tight_layout(); plt.show()Observation: P-only reaches ~52°C but never 60°C. The steady-state error persists because as the error shrinks, the proportional output decreases.
Step 4: Add Integral — Eliminate Offset
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, (Kp, Ki, Kd, label) in zip(axes, [
(5, 0, 0, 'P only'),
(5, 0.5, 0, 'PI'),
(5, 0.5, 2.0, 'PID'),
]):
t, temps, _ = simulate(Kp, Ki, Kd)
ax.plot(t, temps)
ax.axhline(y=60, color='r', linestyle='--')
ax.set_title(label)
ax.set_xlabel('Time (s)'); ax.set_ylabel('°C')
ax.set_ylim(15, 75); ax.grid()
plt.tight_layout(); plt.show()- PI: eliminates steady-state error, but overshoots
- PID: reduces overshoot by damping with the derivative term
Step 5: Ziegler-Nichols Tuning
- Set Ki=Kd=0. Increase Kp until the system oscillates with constant amplitude.
# Find ultimate gain Ku
for Ku in [5, 10, 15, 20, 25]:
t, temps, _ = simulate(Kp=Ku, Ki=0, Kd=0, duration=200)
plt.plot(t, temps, label=f'Kp={Ku}')
plt.axhline(y=60, color='r', linestyle='--')
plt.legend(); plt.grid(); plt.xlabel('Time (s)'); plt.show()- At Ku (sustained oscillation), measure the period Tu.
- Apply Ziegler-Nichols formulas:
# Suppose Ku=20 gives sustained oscillation with period Tu=12s
Ku, Tu = 20, 12
Kp = 0.6 * Ku # 12
Ki = 2 * Kp / Tu # 2
Kd = Kp * Tu / 8 # 18
t, temps, _ = simulate(Kp, Ki, Kd, duration=200)
plt.plot(t, temps)
plt.axhline(y=60, color='r', linestyle='--')
plt.title(f'Ziegler-Nichols: Kp={Kp:.1f}, Ki={Ki:.1f}, Kd={Kd:.1f}')
plt.grid(); plt.show()Ziegler-Nichols gives aggressive tuning (~25% overshoot). Reduce gains by 30-50% for a smoother response.
Step 6: Setpoint Change and Disturbance
def simulate_scenario(Kp, Ki, Kd, duration=200, dt=0.1):
plant = ThermalPlant(dt=dt)
pid = PIDController(Kp=Kp, Ki=Ki, Kd=Kd, dt=dt)
t = np.arange(0, duration, dt)
temps, setpoints = np.zeros_like(t), np.zeros_like(t)
for i in range(len(t)):
# Setpoint: 60°C for first half, 40°C for second half
sp = 60.0 if t[i] < 100 else 40.0
setpoints[i] = sp
# Disturbance: door opens at t=50, ambient drops 10°C for 20s
if 50 < t[i] < 70:
plant.T_ambient = 10.0
else:
plant.T_ambient = 20.0
power = pid.update(sp, plant.T)
plant.step(power)
temps[i] = plant.T
plt.plot(t, temps, label='Temperature')
plt.plot(t, setpoints, 'r--', label='Setpoint')
plt.axvspan(50, 70, alpha=0.1, color='blue', label='Disturbance')
plt.xlabel('Time (s)'); plt.ylabel('°C'); plt.legend(); plt.grid()
plt.title('PID: setpoint change + disturbance rejection')
plt.show()
simulate_scenario(Kp=8, Ki=0.8, Kd=10)The PID controller should recover from both the setpoint change and the disturbance. This is why closed loop control matters — see Open Loop vs Closed Loop.
Exercises
-
Derivative noise: Add Gaussian noise to the temperature measurement (
T + np.random.randn() * 0.5). Watch the derivative term go crazy. Implement a low-pass filter on the derivative. -
Anti-windup test: Set a very high setpoint (200°C) that the heater can’t reach (max 100W). Without anti-windup, the integral term accumulates massively. Compare with and without clamping.
-
Cascaded PID: Add a second integrator to the plant (thermal mass + transport delay). Implement inner and outer PID loops.
-
Discrete implementation in C: Port the PID to C as a struct + update function (like Digital Control). Run it in a timer interrupt at 10 Hz.
Next: 11 - Kalman Filter for Noisy Sensor Data — optimal estimation from noisy measurements.