Chapter 12.2: The Mean Field Game System

Coupling Optimality with Consistency

A Mean Field Game (MFG) couples two PDEs: the HJB equation (backward in time), describing each agent’s optimal strategy given the population density, and the Fokker-Planck (FP) equation (forward in time), describing how the population density evolves under the optimal strategy. The fixed-point condition—the density that agents optimise against is the same density that their optimal actions produce—constitutes the Nash equilibrium of the continuum game.

This framework, introduced independently by Lasry-Lions (2006) and Huang-Caines-Malhamé (2006), provides the mathematical foundation for understanding how individual rational decisions aggregate into collective urban patterns.

12.2.1 The HJB-FP Coupled System

The MFG system consists of two coupled PDEs. The first is the Hamilton-Jacobi-Bellman equation, solved backward from the terminal time$T$:

$$-\frac{\partial u}{\partial t} - \nu \Delta u + \frac{1}{2}|\nabla u|^2 = F(x, \rho)$$

$$u(x, T) = g(x) \quad \text{(terminal condition)}$$

The second is the Fokker-Planck equation, solved forward from the initial time:

$$\frac{\partial \rho}{\partial t} - \nu \Delta \rho - \operatorname{div}(\rho \nabla u) = 0$$

$$\rho(x, 0) = \rho_0(x) \quad \text{(initial condition)}$$

The coupling is through $F(x, \rho)$ in the HJB (the cost depends on the density) and through$\nabla u$ in the FP (the velocity field is the optimal control $\alpha^* = -\nabla u$).

The Fixed-Point Structure

Given a density $\rho$, solve HJB backward to get$u$. Then use$\alpha^* = -\nabla u$ in FP to get a new$\rho'$. At equilibrium,$\rho' = \rho$: the density is self-consistent.

12.2.2 Three Coupling Functions

The coupling function $F(x, \rho)$determines how congestion affects costs. Three canonical choices:

Local Power-Law Coupling

$$F_{\text{local}}(x, \rho) = \rho(x)^\gamma, \qquad \gamma > 0$$

Higher local density means higher cost. The exponent$\gamma$ controls the sensitivity to congestion. For $\gamma = 1$, cost is linear in density (BPR-like). For$\gamma > 1$, cost grows superlinearly, modelling severe congestion effects.

Nonlocal Coupling

$$F_{\text{nonlocal}}(x, \rho) = (K * \rho)(x) = \int K(x - y) \, \rho(y) \, dy$$

Agents care about density in a neighbourhood, not just their exact location. The kernel$K$ represents the range of perception.

Logarithmic Coupling

$$F_{\text{log}}(x, \rho) = \ln \rho(x)$$

This arises when the coupling represents an entropy term. It has the special property of making the MFG system equivalent to a Schrödinger system via Cole-Hopf.

12.2.3 Variational Structure and Well-Posedness

When $F$ is monotone increasing in$\rho$ (congestion aversion), the MFG system has a variational structure. The solution$(u, \rho)$ is a saddle point of the functional:

$$\mathcal{L}[u, \rho] = \int_0^T \!\!\int \left[\rho\!\left(\frac{\partial u}{\partial t} + \nu\Delta u - \frac{1}{2}|\nabla u|^2\right) + \mathcal{F}(x, \rho)\right] dx \, dt$$

where $\mathcal{F}$ satisfies$\partial_\rho \mathcal{F} = F$. For$F = \rho^\gamma$, we have$\mathcal{F} = \rho^{\gamma+1}/(\gamma+1)$.

The monotonicity of $F$ ensuresuniqueness of the MFG equilibrium. This is the Lasry-Lions monotonicity condition. When $F$ is decreasing (agents prefer crowded areas—agglomeration), multiple equilibria can exist.

12.2.4 Stationary MFG: Bid-Rent as Value Function

In the stationary (time-independent) case, setting$\partial_t u = -\bar{H}$ (the ergodic constant), the MFG system reduces to:

$$-\nu \Delta u + \frac{1}{2}|\nabla u|^2 = F(x, \rho) + \bar{H}$$

$$-\nu \Delta \rho - \operatorname{div}(\rho \nabla u) = 0$$

The stationary value function $u(x)$represents the long-run cost of being at location$x$. This is precisely the bid-rent function from urban economics (Module 4): the maximum rent an agent is willing to pay at location $x$ in equilibrium.

The density $\rho(x)$ is the equilibrium population distribution. The MFG system simultaneously determines both the rent gradient and the population distribution—extending Alonso-Muth-Mills to a fully dynamic framework with strategic agents.

12.2.5 Price of Anarchy and Pigouvian Tax

The MFG equilibrium is a Nash equilibrium, but it is generally not socially optimal. Each agent ignores the externality they impose on others through congestion. The Price of Anarchy (PoA) quantifies this inefficiency:

$$\text{PoA} = \frac{\text{Social cost at Nash equilibrium}}{\text{Minimum social cost}} \leq e^{1/\gamma}$$

The bound depends on the congestion exponent$\gamma$. For linear congestion ($\gamma = 1$), PoA$\leq e \approx 2.718$—the Nash equilibrium can be nearly three times worse than the social optimum.

Pigouvian Congestion Tax

To align individual incentives with social welfare, we add a Pigouvian taxthat internalises the congestion externality. The optimal tax at the socially optimal density$\rho^*$ is:

$$\tau(x) = \gamma \, \rho^*(x)^\gamma$$

This equals the marginal externality: the additional cost one agent imposes on all others. With this tax, the Nash equilibrium coincides with the social optimum—the Price of Anarchy becomes 1.

12.2.6 MFG on Networks and Wardrop Limit

The MFG framework extends naturally to networks (graphs). On a graph with nodes$i = 1, \ldots, N$ and edges$e_{ij}$, the discrete MFG system is:

$$u_i = \min_{j \sim i} \left\{c_{ij}(\rho_{ij}) + u_j\right\} \qquad \text{(discrete HJB)}$$

$$\rho_i = \sum_{j: i \to j \text{ optimal}} \rho_j \qquad \text{(discrete FP / flow conservation)}$$

In the limit $\nu \to 0$ (no randomness), agents follow deterministic shortest paths. The MFG equilibrium reduces to the Wardrop equilibrium of traffic assignment: all used routes have equal cost, and no agent can improve by switching. The MFG framework thus provides a stochastic generalisation of the classical Wardrop equilibrium, with$\nu$ playing the role of a “temperature” controlling the randomness of route choice.

12.2.7 Simulation: Solving the MFG System

We solve the stationary MFG system via fixed-point iteration: solve HJB for$u$ given$\rho$, then solve FP for$\rho$ given$u$, and iterate until convergence.

MFG Solver: Fixed-Point Iteration (HJB ↔ FP)

Python

script.py169 lines

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

# ── Stationary MFG solver: 1D periodic domain ──
# HJB: -nu * u_xx + (1/2)|u_x|^2 = rho^gamma + H_bar
# FP: -nu * rho_xx - (rho * u_x)_x = 0

Nx = 200
x = np.linspace(0, 2*np.pi, Nx, endpoint=False)
dx = x[1] - x[0]
nu = 0.15
gamma = 1.0

# External potential (city center attraction)
V_ext = -0.8 * np.cos(x)  # attractive at x = pi

def solve_hjb(rho, nu, gamma, V_ext, n_iter=500, dt_relax=0.001):
    """Solve stationary HJB via relaxation."""
    u = np.zeros(Nx)
    for _ in range(n_iter):
        # Periodic derivatives
        u_x = np.zeros(Nx)
        u_x[1:-1] = (u[2:] - u[:-2]) / (2 * dx)
        u_x[0] = (u[1] - u[-1]) / (2 * dx)
        u_x[-1] = (u[0] - u[-2]) / (2 * dx)

u_xx = np.zeros(Nx)
        u_xx[1:-1] = (u[2:] - 2*u[1:-1] + u[:-2]) / dx**2
        u_xx[0] = (u[1] - 2*u[0] + u[-1]) / dx**2
        u_xx[-1] = (u[0] - 2*u[-1] + u[-2]) / dx**2

# HJB residual: -nu*u_xx + 0.5*|u_x|^2 - rho^gamma - V_ext = H_bar
        rhs = nu * u_xx - 0.5 * u_x**2 + np.maximum(rho, 1e-8)**gamma + V_ext
        H_bar = np.mean(rhs)
        u = u + dt_relax * (rhs - H_bar)
        u = u - np.mean(u)  # fix gauge

return u, H_bar

def solve_fp(u, nu, n_iter=500, dt_relax=0.005):
    """Solve stationary FP: -nu*rho_xx - div(rho*u_x) = 0, periodic."""
    rho = np.ones(Nx) / (2 * np.pi)  # uniform initial

u_x = np.zeros(Nx)
    u_x[1:-1] = (u[2:] - u[:-2]) / (2 * dx)
    u_x[0] = (u[1] - u[-1]) / (2 * dx)
    u_x[-1] = (u[0] - u[-2]) / (2 * dx)

for _ in range(n_iter):
        rho_xx = np.zeros(Nx)
        rho_xx[1:-1] = (rho[2:] - 2*rho[1:-1] + rho[:-2]) / dx**2
        rho_xx[0] = (rho[1] - 2*rho[0] + rho[-1]) / dx**2
        rho_xx[-1] = (rho[0] - 2*rho[-1] + rho[-2]) / dx**2

# div(rho * u_x)
        flux = rho * u_x
        div_flux = np.zeros(Nx)
        div_flux[1:-1] = (flux[2:] - flux[:-2]) / (2 * dx)
        div_flux[0] = (flux[1] - flux[-1]) / (2 * dx)
        div_flux[-1] = (flux[0] - flux[-2]) / (2 * dx)

rho = rho + dt_relax * (nu * rho_xx + div_flux)
        rho = np.maximum(rho, 1e-10)
        rho = rho / (np.sum(rho) * dx)  # normalise

return rho

# ── Fixed-point iteration ──
rho = np.ones(Nx) / (2 * np.pi)
convergence = []

n_fp_iter = 40
for it in range(n_fp_iter):
    rho_old = rho.copy()
    u, H_bar = solve_hjb(rho, nu, gamma, V_ext, n_iter=300, dt_relax=0.002)
    rho = solve_fp(u, nu, n_iter=300, dt_relax=0.003)

# Damped update for stability
    rho = 0.5 * rho + 0.5 * rho_old

err = np.max(np.abs(rho - rho_old))
    convergence.append(err)

# ── Also solve for different gammas ──
gammas = [0.5, 1.0, 2.0, 3.0]
rho_results = {}
u_results = {}

for gam in gammas:
    rho_g = np.ones(Nx) / (2 * np.pi)
    for it in range(30):
        rho_old_g = rho_g.copy()
        u_g, _ = solve_hjb(rho_g, nu, gam, V_ext, n_iter=200, dt_relax=0.002)
        rho_g = solve_fp(u_g, nu, n_iter=200, dt_relax=0.003)
        rho_g = 0.5 * rho_g + 0.5 * rho_old_g
    rho_results[gam] = rho_g
    u_results[gam] = u_g

# ── Plotting ──
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.patch.set_facecolor('#0a0a0a')

colors_g = ['#34d399', '#2dd4bf', '#fbbf24', '#f87171']

# Panel 1: Converged density and value function
ax1 = axes[0, 0]
ax1.set_facecolor('#0a0a0a')
ax1.tick_params(colors='white')
for spine in ax1.spines.values():
    spine.set_color('#334155')
ax1_twin = ax1.twinx()
ax1.plot(x, rho, color='#34d399', linewidth=2.5, label='ρ(x)')
ax1_twin.plot(x, u, color='#fbbf24', linewidth=2.5, linestyle='--', label='u(x)')
ax1.plot(x, -V_ext * 0.3, color='#f472b6', linewidth=1, linestyle=':', alpha=0.5, label='V_ext (scaled)')
ax1.set_xlabel('Position x', color='white')
ax1.set_ylabel('Density ρ', color='#34d399')
ax1_twin.set_ylabel('Value u', color='#fbbf24')
ax1.set_title('MFG Equilibrium (γ=1)', color='#6ee7b7', fontsize=13)
ax1.legend(loc='upper left', facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)
ax1_twin.legend(loc='upper right', facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)
ax1_twin.tick_params(colors='#fbbf24')

# Panel 2: Convergence
ax2 = axes[0, 1]
ax2.set_facecolor('#0a0a0a')
ax2.tick_params(colors='white')
for spine in ax2.spines.values():
    spine.set_color('#334155')
ax2.semilogy(range(len(convergence)), convergence, color='#34d399', linewidth=2.5)
ax2.set_xlabel('Fixed-point iteration', color='white')
ax2.set_ylabel('max|ρ_new - ρ_old|', color='white')
ax2.set_title('Convergence of MFG Iteration', color='#6ee7b7', fontsize=13)

# Panel 3: Density for different gammas
ax3 = axes[1, 0]
ax3.set_facecolor('#0a0a0a')
ax3.tick_params(colors='white')
for spine in ax3.spines.values():
    spine.set_color('#334155')
for i, gam in enumerate(gammas):
    ax3.plot(x, rho_results[gam], color=colors_g[i], linewidth=2,
            label=f'γ = {gam}')
ax3.axhline(y=1/(2*np.pi), color='white', linestyle='--', alpha=0.3, label='Uniform')
ax3.set_xlabel('Position x', color='white')
ax3.set_ylabel('Density ρ(x)', color='white')
ax3.set_title('Equilibrium Density vs Congestion γ', color='#6ee7b7', fontsize=13)
ax3.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)

# Panel 4: Value function for different gammas
ax4 = axes[1, 1]
ax4.set_facecolor('#0a0a0a')
ax4.tick_params(colors='white')
for spine in ax4.spines.values():
    spine.set_color('#334155')
for i, gam in enumerate(gammas):
    ax4.plot(x, u_results[gam], color=colors_g[i], linewidth=2,
            label=f'γ = {gam}')
ax4.set_xlabel('Position x', color='white')
ax4.set_ylabel('Value u(x)', color='white')
ax4.set_title('Value Function vs Congestion γ', color='#6ee7b7', fontsize=13)
ax4.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)

plt.tight_layout(pad=2.0)
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("MFG system solved via fixed-point iteration. Equilibrium density and value function plotted.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Key Takeaways

The MFG system couples HJB (backward, optimal control) with FP (forward, density evolution) through the feedback $F(x, \rho)$.
Three canonical coupling functions:$F = \rho^\gamma$ (local),$F = K * \rho$ (nonlocal),$F = \ln\rho$ (entropic).
Monotone coupling (congestion aversion) guarantees uniqueness via the Lasry-Lions condition.
The stationary MFG value function $u(x)$ is the bid-rent function from urban economics.
Price of Anarchy $\leq e^{1/\gamma}$; Pigouvian tax $\tau = \gamma\rho^{*\gamma}$ restores social optimality.
On networks, MFG reduces to Wardrop equilibrium as $\nu \to 0$.

← HJB Equation MFG Numerics →