Python PID Control Systems — Deep Dive

Discrete PID Implementation

Real controllers run at fixed intervals, not in continuous time. The discrete form uses summation instead of integration and differences instead of derivatives:

class PIDController:
    def __init__(self, kp, ki, kd, dt, output_limits=(-1, 1)):
        self.kp = kp
        self.ki = ki
        self.kd = kd
        self.dt = dt
        self.output_min, self.output_max = output_limits

        self.integral = 0.0
        self.prev_error = 0.0
        self.prev_measurement = 0.0

    def update(self, setpoint, measurement):
        error = setpoint - measurement

        # Proportional
        p_term = self.kp * error

        # Integral with anti-windup (clamping)
        self.integral += error * self.dt
        i_term = self.ki * self.integral

        # Derivative on measurement (avoids derivative kick on setpoint change)
        d_measurement = (measurement - self.prev_measurement) / self.dt
        d_term = -self.kd * d_measurement

        # Total output
        output = p_term + i_term + d_term

        # Output clamping
        if output > self.output_max:
            output = self.output_max
            # Back-calculate integral to prevent windup
            self.integral = (output - p_term - d_term) / self.ki if self.ki != 0 else self.integral
        elif output < self.output_min:
            output = self.output_min
            self.integral = (output - p_term - d_term) / self.ki if self.ki != 0 else self.integral

        self.prev_error = error
        self.prev_measurement = measurement

        return output

Two critical details distinguish this from a textbook implementation:

Derivative on Measurement

The standard formula computes Kd × d(error)/dt. When the setpoint changes suddenly (step change), the error jumps, producing a massive derivative spike called “derivative kick.” Computing the derivative on the measurement signal instead avoids this because the measurement changes smoothly.

Anti-Windup

When the output saturates (hits its limit), the integral keeps accumulating. When conditions change, the inflated integral causes severe overshoot. The back-calculation method above recomputes the integral to match the clamped output, preventing accumulation during saturation.

Derivative Filtering

Raw derivatives amplify high-frequency noise. A first-order low-pass filter smooths the derivative:

class FilteredPID(PIDController):
    def __init__(self, kp, ki, kd, dt, output_limits=(-1, 1), tau_d=0.1):
        super().__init__(kp, ki, kd, dt, output_limits)
        self.tau_d = tau_d  # filter time constant
        self.filtered_derivative = 0.0

    def update(self, setpoint, measurement):
        error = setpoint - measurement

        p_term = self.kp * error
        self.integral += error * self.dt
        i_term = self.ki * self.integral

        # Filtered derivative
        raw_derivative = -(measurement - self.prev_measurement) / self.dt
        alpha = self.dt / (self.tau_d + self.dt)
        self.filtered_derivative = alpha * raw_derivative + (1 - alpha) * self.filtered_derivative
        d_term = self.kd * self.filtered_derivative

        output = p_term + i_term + d_term
        output = np.clip(output, self.output_min, self.output_max)

        self.prev_measurement = measurement
        return output

The filter time constant tau_d controls the tradeoff: larger values remove more noise but add delay to the derivative response. A common rule: set tau_d to 1/8 to 1/10 of the derivative time constant.

Simulating Plant Dynamics

To test PID controllers, you need a simulated “plant” (the system being controlled). Common models:

First-Order System (Temperature Control)

import numpy as np

class ThermalPlant:
    """First-order thermal system: tau * dT/dt = -T + K*u + T_ambient"""
    def __init__(self, tau=60, gain=50, ambient=20, noise_std=0.1):
        self.tau = tau         # time constant (seconds)
        self.gain = gain       # steady-state gain (°C per unit input)
        self.ambient = ambient
        self.noise_std = noise_std
        self.temperature = ambient

    def step(self, control_input, dt):
        dT = (-self.temperature + self.ambient + self.gain * control_input) / self.tau
        self.temperature += dT * dt
        return self.temperature + np.random.normal(0, self.noise_std)

Second-Order System (Motor Position)

class MotorPlant:
    """DC motor position control: J*theta_ddot + b*theta_dot = K*u"""
    def __init__(self, J=0.01, b=0.1, K=0.01):
        self.J = J     # moment of inertia
        self.b = b     # damping
        self.K = K     # motor constant
        self.theta = 0
        self.omega = 0

    def step(self, voltage, dt):
        alpha = (self.K * voltage - self.b * self.omega) / self.J
        self.omega += alpha * dt
        self.theta += self.omega * dt
        return self.theta

Complete Simulation Loop

import matplotlib.pyplot as plt

def simulate_pid(plant, controller, setpoint, duration, dt):
    times = np.arange(0, duration, dt)
    measurements = []
    outputs = []
    setpoints = []

    for t in times:
        sp = setpoint(t) if callable(setpoint) else setpoint
        measurement = plant.step(controller.update(sp, plant.temperature), dt)
        measurements.append(measurement)
        outputs.append(controller.update(sp, measurement))
        setpoints.append(sp)

    return times, measurements, outputs, setpoints

# Example: tune a thermal controller
plant = ThermalPlant(tau=30, gain=40, ambient=20)
pid = PIDController(kp=2.0, ki=0.1, kd=5.0, dt=0.1, output_limits=(0, 1))

times, temps, outputs, sps = simulate_pid(
    plant, pid, setpoint=60.0, duration=300, dt=0.1
)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6))
ax1.plot(times, temps, label='Temperature')
ax1.axhline(y=60, color='r', linestyle='--', label='Setpoint')
ax1.set_ylabel('°C')
ax1.legend()

ax2.plot(times, outputs, label='Control Output')
ax2.set_xlabel('Time (s)')
ax2.set_ylabel('Output')
ax2.legend()
plt.tight_layout()

Auto-Tuning: Relay Feedback Method

The relay feedback method automates Ziegler-Nichols tuning without manual intervention:

def relay_auto_tune(plant, relay_amplitude=1.0, duration=200, dt=0.1):
    """Relay feedback auto-tuning to find ultimate gain and period."""
    setpoint = plant.temperature + 10  # target 10° above current
    output = relay_amplitude
    measurements = []
    outputs = []
    times = np.arange(0, duration, dt)

    for t in times:
        measurement = plant.step(output, dt)
        measurements.append(measurement)

        # Relay: full positive when below setpoint, full negative when above
        error = setpoint - measurement
        output = relay_amplitude if error > 0 else -relay_amplitude
        outputs.append(output)

    # Find oscillation parameters from last half of data
    measurements = np.array(measurements)
    half = len(measurements) // 2
    recent = measurements[half:]

    # Peak-to-peak amplitude of oscillation
    amplitude = (recent.max() - recent.min()) / 2

    # Find period by counting zero crossings
    mean_val = recent.mean()
    crossings = np.where(np.diff(np.sign(recent - mean_val)))[0]
    if len(crossings) >= 4:
        periods = np.diff(crossings[::2]) * dt * 2
        Tu = np.mean(periods)
    else:
        return None  # could not determine period

    # Ultimate gain from describing function analysis
    Ku = 4 * relay_amplitude / (np.pi * amplitude)

    # Ziegler-Nichols PID tuning
    return {
        'Ku': Ku, 'Tu': Tu,
        'Kp': 0.6 * Ku,
        'Ki': 1.2 * Ku / Tu,
        'Kd': 0.075 * Ku * Tu
    }

Using the python-control Library

For transfer function analysis and advanced control design:

import control

# Define plant as transfer function: G(s) = K / (tau*s + 1)
K, tau = 40, 30
plant_tf = control.tf([K], [tau, 1])

# PID controller: C(s) = Kp + Ki/s + Kd*s
kp, ki, kd = 2.0, 0.1, 5.0
pid_tf = control.tf([kd, kp, ki], [1, 0])

# Closed-loop transfer function
closed_loop = control.feedback(pid_tf * plant_tf)

# Step response analysis
t, y = control.step_response(closed_loop, T=np.linspace(0, 300, 3000))

# Stability margins
gm, pm, wgc, wpc = control.margin(pid_tf * plant_tf)
print(f"Gain margin: {20*np.log10(gm):.1f} dB, Phase margin: {pm:.1f}°")

Bode Plot Analysis

# Bode plot of open-loop system
control.bode_plot(pid_tf * plant_tf, dB=True)

# Root locus
control.root_locus(pid_tf * plant_tf)

Cascaded PID

Complex systems use multiple PID loops. A motor position controller typically has:

• Outer loop: Position PID (setpoint = desired angle, measurement = encoder)
• Inner loop: Velocity PID (setpoint = outer loop output, measurement = tachometer)

The inner loop runs faster (1 kHz) than the outer loop (100 Hz). This cascade responds faster and handles disturbances better than a single position loop.

class CascadePID:
    def __init__(self, outer_gains, inner_gains, dt_outer, dt_inner):
        self.outer = PIDController(*outer_gains, dt=dt_outer)
        self.inner = PIDController(*inner_gains, dt=dt_inner)
        self.inner_ratio = int(dt_outer / dt_inner)

    def update(self, position_setpoint, position, velocity, inner_plant_step):
        velocity_setpoint = self.outer.update(position_setpoint, position)
        for _ in range(self.inner_ratio):
            output = self.inner.update(velocity_setpoint, velocity)
            velocity = inner_plant_step(output)
        return output

Performance Metrics

When comparing tuning configurations, measure:

Metric	Formula	Good For
Rise time	Time from 10% to 90% of setpoint	Speed
Overshoot	(peak - setpoint) / setpoint × 100%	Safety
Settling time	Time to stay within ±2% of setpoint	Stability
ISE	∫e² dt	Penalizes large errors
ITAE	∫t·|e| dt	Penalizes persistent errors

Automated tuning can minimize these metrics using scipy.optimize:

from scipy.optimize import minimize

def cost_function(params):
    kp, ki, kd = params
    plant = ThermalPlant()
    pid = PIDController(kp, ki, kd, dt=0.1)
    _, temps, _, _ = simulate_pid(plant, pid, 60, 300, 0.1)
    errors = np.array(temps) - 60
    times = np.arange(len(errors)) * 0.1
    return np.sum(times * np.abs(errors) * 0.1)  # ITAE

result = minimize(cost_function, x0=[1, 0.1, 1], bounds=[(0, 10), (0, 1), (0, 20)])
print(f"Optimal gains: Kp={result.x[0]:.3f}, Ki={result.x[1]:.3f}, Kd={result.x[2]:.3f}")

One thing to remember: A production PID controller needs derivative-on-measurement to avoid kick, anti-windup to prevent integral saturation, and derivative filtering to handle sensor noise — three details that separate textbook PID from working PID.