LQR & Riccati Equation

When urban dynamics can be linearised and costs are quadratic, the optimal control problem admits a beautiful closed-form solution: the Linear Quadratic Regulator (LQR). The optimal feedback gains come from solving the algebraic Riccati equation, which we derive and apply to traffic signal control of queue lengths.

1. Linear State-Space Model

Consider a network of $n$ signalised intersections with queue lengths assembled into the state vector $\mathbf{x}(t) \in \mathbb{R}^n$ and green time allocations forming the control vector $\mathbf{u}(t) \in \mathbb{R}^m$. Linearising the store-and-forward model about a nominal operating point gives:

$$\frac{d\mathbf{x}}{dt} = A\mathbf{x} + B\mathbf{u}$$

For a corridor with three intersections, the system matrix $A$ captures queue interactions and the input matrix $B$ maps green time to queue discharge:

$$A = \begin{pmatrix} -\mu_1 & 0 & 0 \\ \mu_1 & -\mu_2 & 0 \\ 0 & \mu_2 & -\mu_3 \end{pmatrix}, \quad B = \begin{pmatrix} -\mu_1 & 0 & 0 \\ 0 & -\mu_2 & 0 \\ 0 & 0 & -\mu_3 \end{pmatrix}$$

The diagonal of $A$ represents queue discharge at saturation flow, while the sub-diagonal entries capture platoon transfer from upstream to downstream intersections.

2. Quadratic Cost Functional

The LQR objective balances two competing goals: minimise queue lengths (delay) and minimise control effort (signal changes). The infinite-horizon cost functional is:

$$J = \int_0^\infty \left( \mathbf{x}^T Q \mathbf{x} + \mathbf{u}^T R \mathbf{u} \right) dt$$

The weight matrices encode engineering priorities:

State weight $Q \succeq 0$: Penalises queue lengths. Diagonal entries weight each intersection; off-diagonal entries can penalise differences between neighbouring queues (equity).
Control weight $R \succ 0$: Penalises aggressive signal changes. Larger $R$ produces smoother, less responsive control. Must be positive definite.

The ratio $Q/R$ determines how aggressively the controller acts: large ratio means fast queue reduction at the cost of frequent signal changes; small ratio means gentle control with larger residual queues.

3. Deriving the Optimal Feedback Law

We apply the PMP from the previous chapter. The Hamiltonian is:

$$H = \mathbf{x}^T Q \mathbf{x} + \mathbf{u}^T R \mathbf{u} + \boldsymbol{\lambda}^T (A\mathbf{x} + B\mathbf{u})$$

Setting $\partial H / \partial \mathbf{u} = 0$:

$$2R\mathbf{u} + B^T \boldsymbol{\lambda} = 0 \implies \mathbf{u}^* = -\tfrac{1}{2} R^{-1} B^T \boldsymbol{\lambda}$$

The key ansatz: assume the costate is linear in the state, $\boldsymbol{\lambda} = 2P\mathbf{x}$, where $P$ is a symmetric positive definite matrix. Then:

$$\mathbf{u}^* = -R^{-1} B^T P \mathbf{x}$$

This is a state feedback law: the optimal control is a linear function of the current queue lengths, with gain matrix $K = R^{-1} B^T P$.

4. The Algebraic Riccati Equation

Substituting the feedback law into the costate equation and requiring consistency gives the Continuous Algebraic Riccati Equation (CARE):

$$PA + A^T P - PBR^{-1}B^T P + Q = 0$$

This is a matrix equation for the unknown $P \in \mathbb{R}^{n \times n}$. Key properties:

Existence and uniqueness: A unique stabilising solution $P \succeq 0$ exists if and only if $(A, B)$ is stabilisable and $(A, Q^{1/2})$ is detectable.
Closed-loop stability: The eigenvalues of $A - BR^{-1}B^T P$ all have negative real parts.
Optimality certificate: The optimal cost from state $\mathbf{x}_0$ is $J^* = \mathbf{x}_0^T P \mathbf{x}_0$.

To derive CARE, substitute $\boldsymbol{\lambda} = 2P\mathbf{x}$ into the costate ODE:

$$\dot{\boldsymbol{\lambda}} = 2P\dot{\mathbf{x}} = 2P(A - BR^{-1}B^TP)\mathbf{x}$$

From PMP, $\dot{\boldsymbol{\lambda}} = -2Q\mathbf{x} - 2A^TP\mathbf{x}$. Equating and requiring the relation to hold for all $\mathbf{x}$ gives CARE.

5. Solving the Riccati Equation with SciPy

We solve CARE for a 4-intersection network, compute the LQR gain, simulate the closed-loop response, and compare it to an uncontrolled system.

LQR Traffic Signal Control via Riccati Equation

Python

script.py136 lines

import numpy as np
from scipy.linalg import solve_continuous_are
from scipy.integrate import solve_ivp
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# --- 4-Intersection Corridor ---
n = 4
mu = np.array([0.5, 0.45, 0.5, 0.48])  # saturation flows

# System matrix: queue dynamics with upstream platoon transfer
A = np.zeros((n, n))
for i in range(n):
    A[i, i] = -mu[i] * 0.3  # partial self-discharge
    if i > 0:
        A[i, i-1] = mu[i-1] * 0.2  # upstream transfer

# Input matrix: green time -> queue reduction
B = np.diag(-mu)

# Cost matrices
Q = np.diag([4.0, 3.0, 3.0, 2.0])  # prioritize upstream intersections
R = np.eye(n) * 0.5                  # moderate control cost

print("=== System Matrices ===")
print(f"A =\n{A}")
print(f"\nB = diag({np.diag(B)})")
print(f"\nQ = diag({np.diag(Q)})")
print(f"R = diag({np.diag(R)})")

# --- Solve CARE ---
P = solve_continuous_are(A, B, Q, R)
print(f"\n=== Riccati Solution P ===")
print(f"P =\n{np.round(P, 4)}")

# LQR gain
K = np.linalg.solve(R, B.T @ P)
print(f"\n=== LQR Gain K ===")
print(f"K =\n{np.round(K, 4)}")

# Closed-loop eigenvalues
A_cl = A - B @ K
eigs = np.linalg.eigvals(A_cl)
print(f"\nClosed-loop eigenvalues: {np.round(eigs, 4)}")
print(f"All stable: {np.all(np.real(eigs) < 0)}")

# --- Simulate ---
x0 = np.array([15.0, 12.0, 10.0, 8.0])  # initial queues
T_sim = 60.0

# Controlled (LQR)
def lqr_dynamics(t, x):
    u = -K @ x
    u = np.clip(u, -1, 1)  # saturate control
    return A @ x + B @ u

sol_ctrl = solve_ivp(lqr_dynamics, [0, T_sim], x0, max_step=0.5, dense_output=True)

# Uncontrolled
def open_dynamics(t, x):
    return A @ x

sol_open = solve_ivp(open_dynamics, [0, T_sim], x0, max_step=0.5, dense_output=True)

# Evaluate on fine grid
t_plot = np.linspace(0, T_sim, 500)
x_ctrl = sol_ctrl.sol(t_plot)
x_open = sol_open.sol(t_plot)

# Optimal cost
J_star = x0 @ P @ x0
print(f"\nOptimal cost J* = x0^T P x0 = {J_star:.2f}")

# --- Plot ---
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
colors = ['#10b981', '#14b8a6', '#06b6d4', '#0ea5e9']
labels = [f'Int {i+1}' for i in range(n)]

# Queue lengths
ax = axes[0, 0]
for i in range(n):
    ax.plot(t_plot, x_ctrl[i], color=colors[i], linewidth=2, label=labels[i])
    ax.plot(t_plot, x_open[i], color=colors[i], linewidth=1, linestyle='--', alpha=0.5)
ax.set_xlabel('Time (s)', color='white')
ax.set_ylabel('Queue Length', color='white')
ax.set_title('Queue Evolution: LQR (solid) vs Open-loop (dashed)', color='white', fontweight='bold')
ax.legend(fontsize=9)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Control effort
ax = axes[0, 1]
u_ctrl = np.zeros((n, len(t_plot)))
for j, t_val in enumerate(t_plot):
    x_val = sol_ctrl.sol(t_val)
    u_ctrl[:, j] = np.clip(-K @ x_val, -1, 1)
for i in range(n):
    ax.plot(t_plot, u_ctrl[i], color=colors[i], linewidth=2, label=labels[i])
ax.set_xlabel('Time (s)', color='white')
ax.set_ylabel('Control u(t)', color='white')
ax.set_title('Optimal Control Signals', color='white', fontweight='bold')
ax.legend(fontsize=9)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Riccati matrix heatmap
ax = axes[1, 0]
im = ax.imshow(P, cmap='BuGn', aspect='auto')
ax.set_title('Riccati Solution P', color='white', fontweight='bold')
ax.set_xlabel('State index', color='white')
ax.set_ylabel('State index', color='white')
ax.tick_params(colors='white')
plt.colorbar(im, ax=ax)

# Phase portrait (first two states)
ax = axes[1, 1]
ax.plot(x_ctrl[0], x_ctrl[1], color='#10b981', linewidth=2, label='LQR')
ax.plot(x_open[0], x_open[1], color='#ef4444', linewidth=2, linestyle='--', label='Open-loop')
ax.plot(x0[0], x0[1], 'wo', markersize=10, label='Start')
ax.set_xlabel('Queue 1', color='white')
ax.set_ylabel('Queue 2', color='white')
ax.set_title('Phase Portrait (Int 1 vs Int 2)', color='white', fontweight='bold')
ax.legend(fontsize=9)
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: LQR traffic control simulation")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

6. Pseudo-Time Integration for the Riccati Equation

An alternative to direct algebraic solvers is pseudo-time integration: we introduce a fictitious time $\tau$ and solve the matrix differential equation:

$$\frac{dP}{d\tau} = -\left( PA + A^T P - PBR^{-1}B^T P + Q \right)$$

Starting from any $P(0) \succeq 0$, this flow converges to the steady state $dP/d\tau = 0$, which is precisely the CARE solution. The right-hand side is the negative of the Riccati residual, so convergence corresponds to the residual vanishing.

This approach is implemented in Fortran below using explicit Euler integration with symmetry enforcement at each step.

Pseudo-Time Integration for Riccati Equation (Fortran-style approach in Python)

Python

script.py94 lines

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

# --- Pseudo-time integration for Riccati equation ---
n = 3
mu = np.array([0.5, 0.45, 0.5])

A = np.array([
    [-0.15, 0.0,  0.0],
    [ 0.10, -0.14, 0.0],
    [ 0.0,   0.09, -0.15]
])
B = np.diag(-mu)
Q = np.diag([3.0, 3.0, 2.0])
R = np.eye(n) * 0.5
R_inv = np.linalg.inv(R)
BRBt = B @ R_inv @ B.T

# Pseudo-time parameters
d_tau = 0.01
n_steps = 5000
P = np.eye(n) * 0.1  # initial guess

residual_history = []

for step in range(n_steps):
    # Riccati residual
    residual = P @ A + A.T @ P - P @ BRBt @ P + Q

# Update
    P = P - d_tau * residual

# Enforce symmetry
    P = 0.5 * (P + P.T)

res_norm = np.linalg.norm(residual, 'fro')
    residual_history.append(res_norm)

if res_norm < 1e-10:
        print(f"Converged at step {step+1}, residual = {res_norm:.2e}")
        break

print(f"\n=== Pseudo-Time Riccati Solution ===")
print(f"P =\n{np.round(P, 6)}")

# Verify against scipy
from scipy.linalg import solve_continuous_are
P_exact = solve_continuous_are(A, B, Q, R)
print(f"\nP (scipy) =\n{np.round(P_exact, 6)}")
print(f"\nMax difference: {np.max(np.abs(P - P_exact)):.2e}")

# LQR gain from pseudo-time solution
K = np.linalg.solve(R, B.T @ P)
print(f"\nLQR Gain K =\n{np.round(K, 4)}")
eigs = np.linalg.eigvals(A - B @ K)
print(f"Closed-loop eigenvalues: {np.round(eigs, 4)}")

# --- Plot convergence ---
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
ax.semilogy(residual_history, color='#10b981', linewidth=2)
ax.set_xlabel('Pseudo-time step', fontsize=12, color='white')
ax.set_ylabel('Riccati residual (Frobenius norm)', fontsize=12, color='white')
ax.set_title('Convergence of Pseudo-Time Integration', fontsize=14, color='white', fontweight='bold')
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.2)

# Plot P matrix evolution would require storing history; instead show gain comparison
ax = axes[1]
K_labels = []
K_pseudo = K.flatten()
K_scipy = (np.linalg.solve(R, B.T @ P_exact)).flatten()
x_pos = np.arange(len(K_pseudo))
width = 0.35
ax.bar(x_pos - width/2, K_pseudo, width, color='#10b981', label='Pseudo-time', alpha=0.8)
ax.bar(x_pos + width/2, K_scipy, width, color='#06b6d4', label='SciPy CARE', alpha=0.8)
ax.set_xlabel('Gain matrix entry (flattened)', fontsize=12, color='white')
ax.set_ylabel('Value', fontsize=12, color='white')
ax.set_title('LQR Gains: Pseudo-time vs SciPy', fontsize=14, 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: pseudo-time Riccati convergence")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7. Robustness Properties of LQR

The LQR controller enjoys remarkable robustness guarantees, often called the “Kalman robustness margins”:

Gain margin: The closed-loop system remains stable for any gain perturbation in the range $[\frac{1}{2}, \infty)$.
Phase margin: At least $60°$ per channel.
Disk margin: The Nyquist plot of each loop avoids a disk of radius $\frac{1}{2}$ centred at $-1$.

These margins mean the LQR signal controller tolerates substantial modelling errors: if the actual saturation flow rates differ by up to a factor of 2 from the design values, the system remains stable. This is critical for urban applications where traffic demand is inherently uncertain.

Key Takeaway

LQR provides the optimal linear feedback law for quadratic cost problems. The Riccati equation encodes the full future cost into a single matrix $P$, enabling real-time computation of optimal green times from measured queue lengths. The built-in robustness margins make it practical for the uncertain environment of real traffic networks.

← Pontryagin Maximum Principle Congestion Pricing & MPC →