Congestion Pricing & Model Predictive Control

Congestion pricing internalises the externality that each driver imposes on others. The Pigouvian toll emerges from the BPR travel-time function as the marginal external cost. Model Predictive Control (MPC) then provides the real-time computational framework for adaptive signal timing, receding the optimisation horizon at each step.

1. Bureau of Public Roads (BPR) Travel-Time Function

The BPR function is the standard link performance model in transportation engineering. It relates travel time to traffic flow on a link with capacity $C$:

$$t(f) = t_0 \left( 1 + \alpha \left( \frac{f}{C} \right)^\beta \right)$$

The standard calibration parameters are $\alpha = 0.15$ and $\beta = 4$, where:

$t_0$ is the free-flow travel time.
$f$ is the traffic flow (vehicles per hour).
$C$ is the practical capacity.
The exponent $\beta = 4$ produces sharp congestion onset near capacity.

The total cost experienced by all users on the link is:

$$\text{TC}(f) = f \cdot t(f) = f \, t_0 \left( 1 + \alpha \left( \frac{f}{C} \right)^\beta \right)$$

2. The Pigouvian Toll

Each driver perceives only their private cost $t(f)$ but imposes an additional delay on all other $f$ drivers by marginally increasing congestion. The marginal social cost is:

$$\text{MSC}(f) = \frac{d}{df}\big[f \cdot t(f)\big] = t(f) + f \cdot \frac{dt}{df}$$

The second term is the marginal external cost—the damage each additional vehicle inflicts on everyone else. The Pigouvian toll equals this externality:

$$\tau^* = f \cdot \frac{dt}{df} = f \cdot t_0 \frac{\alpha \beta}{C} \left( \frac{f}{C} \right)^{\beta - 1} = t_0 \alpha \beta \left( \frac{f}{C} \right)^\beta \cdot \frac{f}{C}$$

Simplifying with the BPR parameters:

$$\tau^* = t_0 \alpha \beta \frac{f^{\beta+1}}{C^{\beta+1}}$$

At the social optimum, each driver pays the toll $\tau^*$ on top of their travel time, so their perceived cost equals the marginal social cost. This aligns private incentives with social welfare:

$$\underbrace{t(f) + \tau^*}_{\text{perceived cost}} = \underbrace{t(f) + f \cdot t'(f)}_{\text{marginal social cost}} = \text{MSC}(f)$$

User Equilibrium vs Social Optimum

Without tolls, drivers reach a user equilibrium (Wardrop's first principle) where no driver can unilaterally reduce their travel time. With the Pigouvian toll, the user equilibrium coincides with the system optimum (Wardrop's second principle), minimising total network delay.

3. Model Predictive Control Framework

MPC is a receding-horizon strategy that solves a finite-horizon optimal control problem at each time step, applies only the first control action, then re-solves with updated measurements. The MPC algorithm at time $t_k$:

$$\min_{\mathbf{u}_{k}, \dots, \mathbf{u}_{k+N-1}} \sum_{j=0}^{N-1} \ell(\mathbf{x}_{k+j}, \mathbf{u}_{k+j}) + V_f(\mathbf{x}_{k+N})$$

subject to system dynamics, state constraints (non-negative queues, maximum capacities), and control constraints (minimum/maximum green times). The components are:

Stage cost $\ell(\mathbf{x}, \mathbf{u})$: Weighted sum of delay and control effort, as in the LQR formulation.
Terminal cost $V_f(\mathbf{x})$: Approximation of cost-to-go beyond the horizon. Often set to $\mathbf{x}^T P \mathbf{x}$ where $P$ is the LQR Riccati solution.
Prediction horizon $N$: Longer horizons improve optimality but increase computation. Typical: 5-20 signal cycles.

The receding horizon provides feedback: by re-solving at each step with fresh measurements, MPC naturally compensates for model errors and demand fluctuations. It bridges the gap between the theoretical optimality of PMP/LQR and practical real-time implementation.

4. Control Strategy Comparison

We compare three signal control strategies for a corridor under time-varying demand:

Fixed-time: Pre-computed green splits based on average demand. No feedback. The standard baseline.
Actuated: Local feedback using vehicle detectors. Extends green for detected vehicles up to a maximum. No coordination.
MPC: Centralised optimisation using predicted demand over a rolling horizon. Full coordination with feedback.

The simulation below implements all three and plots the resulting queue evolution, total delay, and congestion pricing toll under each strategy.

MPC vs Fixed-Time vs Actuated Control with Congestion Pricing

Python

script.py193 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

np.random.seed(42)

# --- Network Setup ---
n_int = 3
T_sim = 300.0  # seconds
dt = 1.0
N_steps = int(T_sim / dt)
t = np.arange(0, T_sim, dt)

mu = np.array([0.5, 0.45, 0.5])  # saturation flows
C = mu  # capacity = saturation flow
t0_free = np.array([20.0, 25.0, 22.0])  # free-flow travel times

# Time-varying demand (morning peak)
def demand(t_val):
    peak = 0.35 + 0.25 * np.exp(-((t_val - 150)**2) / (2 * 50**2))
    return np.array([peak, peak * 0.8, peak * 0.7])

# BPR function
alpha_bpr, beta_bpr = 0.15, 4

def bpr_time(f, t0, cap):
    return t0 * (1 + alpha_bpr * (f / cap)**beta_bpr)

def pigouvian_toll(f, t0, cap):
    return t0 * alpha_bpr * beta_bpr * (f / cap)**(beta_bpr) * f / cap

# --- Strategy 1: Fixed-Time Control ---
def simulate_fixed(cycle=60, green_ratio=0.5):
    x = np.zeros((n_int, N_steps))
    x[:, 0] = 5.0
    total_delay = 0
    for k in range(N_steps - 1):
        phase = (t[k] % cycle) / cycle
        g = np.ones(n_int) * (1.0 if phase < green_ratio else 0.0)
        d = demand(t[k])
        for i in range(n_int):
            s_i = d[i]
            d_i = g[i] * mu[i] * (1 if x[i, k] > 0.01 else 0)
            x[i, k+1] = max(0, x[i, k] + dt * (s_i - d_i))
        total_delay += np.sum(x[:, k]) * dt
    return x, total_delay

# --- Strategy 2: Actuated Control ---
def simulate_actuated(min_green=10, max_green=40, gap_threshold=3.0):
    x = np.zeros((n_int, N_steps))
    x[:, 0] = 5.0
    green_timer = np.zeros(n_int)
    is_green = np.ones(n_int, dtype=bool)
    total_delay = 0
    for k in range(N_steps - 1):
        d = demand(t[k])
        g = np.zeros(n_int)
        for i in range(n_int):
            green_timer[i] += dt
            if is_green[i]:
                if x[i, k] > 0.5 and green_timer[i] < max_green:
                    g[i] = 1.0
                elif green_timer[i] >= min_green:
                    is_green[i] = False
                    green_timer[i] = 0
                else:
                    g[i] = 1.0
            else:
                if green_timer[i] >= min_green:
                    is_green[i] = True
                    green_timer[i] = 0
                    g[i] = 1.0

s_i = d[i]
            d_i = g[i] * mu[i] * (1 if x[i, k] > 0.01 else 0)
            x[i, k+1] = max(0, x[i, k] + dt * (s_i - d_i))
        total_delay += np.sum(x[:, k]) * dt
    return x, total_delay

# --- Strategy 3: MPC Control ---
def simulate_mpc(horizon=20):
    x = np.zeros((n_int, N_steps))
    x[:, 0] = 5.0
    total_delay = 0
    for k in range(N_steps - 1):
        d = demand(t[k])
        # Simple MPC: look ahead and choose green for most congested
        # Greedy horizon optimisation
        best_g = np.zeros(n_int)
        # Allocate green proportional to queue + predicted arrivals
        queue_priority = x[:, k] + d * horizon * dt
        if np.sum(queue_priority) > 0:
            best_g = queue_priority / np.sum(queue_priority)
            best_g = np.clip(best_g, 0, 1)

for i in range(n_int):
            s_i = d[i]
            d_i = best_g[i] * mu[i] * (1 if x[i, k] > 0.01 else 0)
            x[i, k+1] = max(0, x[i, k] + dt * (s_i - d_i))
        total_delay += np.sum(x[:, k]) * dt
    return x, total_delay

# Run all strategies
x_fixed, delay_fixed = simulate_fixed()
x_act, delay_act = simulate_actuated()
x_mpc, delay_mpc = simulate_mpc()

print("=== Total Delay Comparison ===")
print(f"Fixed-time:  {delay_fixed:.1f} veh-seconds")
print(f"Actuated:    {delay_act:.1f} veh-seconds")
print(f"MPC:         {delay_mpc:.1f} veh-seconds")
print(f"\nMPC improvement over fixed: {(1 - delay_mpc/delay_fixed)*100:.1f}%")
print(f"MPC improvement over actuated: {(1 - delay_mpc/delay_act)*100:.1f}%")

# Compute tolls
flows = np.array([demand(tk) for tk in t])  # shape (N_steps, n_int)
tolls = np.zeros_like(flows)
for k in range(N_steps):
    for i in range(n_int):
        tolls[k, i] = pigouvian_toll(flows[k, i], t0_free[i], C[i])

print(f"\n=== Pigouvian Toll at Peak ===")
peak_idx = np.argmax(flows[:, 0])
for i in range(n_int):
    print(f"Intersection {i+1}: toll = {tolls[peak_idx, i]:.2f} time-units at f/C = {flows[peak_idx, i]/C[i]:.2f}")

# --- Plot ---
fig, axes = plt.subplots(2, 2, figsize=(12, 9))
colors = ['#10b981', '#14b8a6', '#06b6d4']

# Total queue comparison
ax = axes[0, 0]
ax.plot(t, np.sum(x_fixed, axis=0), color='#ef4444', linewidth=2, label='Fixed-time')
ax.plot(t, np.sum(x_act, axis=0), color='#f59e0b', linewidth=2, label='Actuated')
ax.plot(t, np.sum(x_mpc, axis=0), color='#10b981', linewidth=2, label='MPC')
ax.fill_between(t, 0, np.sum(x_mpc, axis=0), color='#10b981', alpha=0.1)
ax.set_xlabel('Time (s)', color='white')
ax.set_ylabel('Total Queue (veh)', color='white')
ax.set_title('Total Network Queue', color='white', fontweight='bold')
ax.legend(fontsize=10)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Per-intersection MPC queues
ax = axes[0, 1]
for i in range(n_int):
    ax.plot(t, x_mpc[i], color=colors[i], linewidth=2, label=f'Int {i+1}')
ax.set_xlabel('Time (s)', color='white')
ax.set_ylabel('Queue Length (veh)', color='white')
ax.set_title('MPC: Per-Intersection Queues', color='white', fontweight='bold')
ax.legend(fontsize=10)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# BPR curve and toll
ax = axes[1, 0]
f_range = np.linspace(0, 1.2, 200)
t_bpr = 1 + alpha_bpr * f_range**beta_bpr
toll_bpr = alpha_bpr * beta_bpr * f_range**(beta_bpr + 1)
msc = t_bpr + toll_bpr
ax.plot(f_range, t_bpr, color='#10b981', linewidth=2, label='Private cost t(f)')
ax.plot(f_range, msc, color='#ef4444', linewidth=2, label='Marginal social cost')
ax.fill_between(f_range, t_bpr, msc, color='#f59e0b', alpha=0.3, label='Pigouvian toll')
ax.axvline(x=1.0, color='white', linestyle='--', alpha=0.3, label='Capacity')
ax.set_xlabel('Flow / Capacity (f/C)', color='white')
ax.set_ylabel('Cost (time units)', color='white')
ax.set_title('BPR Curve & Pigouvian Toll', color='white', fontweight='bold')
ax.legend(fontsize=9)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Toll over time
ax = axes[1, 1]
for i in range(n_int):
    ax.plot(t, tolls[:, i], color=colors[i], linewidth=2, label=f'Int {i+1}')
ax.set_xlabel('Time (s)', color='white')
ax.set_ylabel('Pigouvian Toll', color='white')
ax.set_title('Dynamic Congestion Pricing', color='white', fontweight='bold')
ax.legend(fontsize=10)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

fig.patch.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=120, bbox_inches='tight', facecolor='#0f172a')
plt.close()
print("\nPlot saved: control strategy comparison with congestion pricing")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

5. MPC Stability and Convergence

A critical question for any receding-horizon scheme is stability: does repeatedly solving a finite-horizon problem and applying only the first step produce a stable closed loop? The answer relies on choosing the terminal cost $V_f$ and terminal constraint set appropriately.

The standard stability result uses the LQR cost-to-go as the terminal ingredient. Define the Lyapunov function candidate:

$$V(\mathbf{x}_k) = \min_{\mathbf{u}} \sum_{j=0}^{N-1} \ell(\mathbf{x}_{k+j}, \mathbf{u}_{k+j}) + \mathbf{x}_{k+N}^T P \mathbf{x}_{k+N}$$

If $P$ satisfies the Riccati equation and the terminal set is the LQR-invariant ellipsoid $\mathcal{X}_f = \{\mathbf{x} : \mathbf{x}^T P \mathbf{x} \leq \gamma\}$, then:

$$V(\mathbf{x}_{k+1}) - V(\mathbf{x}_k) \leq -\ell(\mathbf{x}_k, \mathbf{u}_k^*)$$

This guarantees that $V$ decreases at each step by at least the stage cost, ensuring asymptotic stability of the origin (zero queues).

Practical Implications

MPC bridges theory and practice: it inherits LQR's optimality guarantees for the linearised system while handling nonlinear dynamics, constraints, and time-varying demand through repeated online optimisation. Combined with Pigouvian pricing, it provides a complete framework for real-time, economically efficient traffic management.

6. Welfare Analysis of Congestion Pricing

The welfare gain from congestion pricing can be computed as the area between the private cost and marginal social cost curves, evaluated at the user equilibrium and social optimum flows:

$$\Delta W = \int_{f^{\text{SO}}}^{f^{\text{UE}}} \big[\text{MSC}(f) - t(f)\big] \, df = \int_{f^{\text{SO}}}^{f^{\text{UE}}} f \cdot t'(f) \, df$$

For BPR with $\beta = 4$, this integral has a closed form. The deadweight loss of unpriced congestion grows as $O((f/C)^{\beta+2})$, meaning the welfare case for pricing strengthens dramatically as demand approaches capacity.

Revenue from the toll, $\mathcal{R} = f^{\text{SO}} \cdot \tau^*$, can be recycled to fund transit improvements, reduce other taxes, or provide rebates to low-income travellers—addressing the equity concerns that often accompany congestion pricing proposals.

← LQR & Riccati Equation Part 14: Kuramoto Synchronization →