Social Force Model

Dirk Helbing's Social Force Model treats pedestrians as particles subject to Newtonian-like forces: a driving force toward their destination, repulsive forces from other pedestrians and walls, and optional attractive forces. The model famously reproduces emergent phenomena like lane formation, oscillation at bottlenecks, and faster-is-slower effects in evacuations.

1. Equation of Motion

Each pedestrian $i$ of mass $m_i$ has position $\mathbf{r}_i$ and velocity $\mathbf{v}_i$. The equation of motion is:

$$m_i \frac{d\mathbf{v}_i}{dt} = \mathbf{F}_i^{\text{desired}} + \sum_{j \neq i} \mathbf{F}_{ij}^{\text{rep}} + \sum_w \mathbf{F}_{iw}^{\text{wall}} + \boldsymbol{\xi}_i$$

where $\boldsymbol{\xi}_i$ is a stochastic fluctuation term representing individual variations in behaviour. Each force component has a distinct physical interpretation.

2. Desired (Driving) Force

The desired force models the pedestrian's intention to walk toward their destination at a preferred speed $v_i^0$ in direction $\mathbf{e}_i$:

$$\mathbf{F}_i^{\text{desired}} = m_i \frac{v_i^0 \, \mathbf{e}_i - \mathbf{v}_i}{\tau_i}$$

Here $\tau_i$ is the relaxation time, typically 0.5 seconds, representing how quickly the pedestrian adjusts to their desired velocity. The direction $\mathbf{e}_i$ points from the current position to the next waypoint:

$$\mathbf{e}_i = \frac{\mathbf{r}_i^{\text{goal}} - \mathbf{r}_i}{|\mathbf{r}_i^{\text{goal}} - \mathbf{r}_i|}$$

This force continuously accelerates the pedestrian toward $v_i^0 \mathbf{e}_i$. Without other forces, the velocity exponentially relaxes:$\; \mathbf{v}_i(t) = v_i^0 \mathbf{e}_i + (\mathbf{v}_i(0) - v_i^0 \mathbf{e}_i) e^{-t/\tau_i}$.

3. Pedestrian Repulsion

Pedestrians maintain a personal space. The repulsive force between pedestrians $i$ and $j$ decays exponentially with distance:

$$\mathbf{F}_{ij}^{\text{rep}} = A \exp\!\left(\frac{r_{ij} - d_{ij}}{B}\right) \hat{\mathbf{n}}_{ij}$$

where:

$r_{ij} = r_i + r_j$ is the sum of body radii
$d_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ is the centre-to-centre distance
$\hat{\mathbf{n}}_{ij} = (\mathbf{r}_i - \mathbf{r}_j) / d_{ij}$ is the unit normal from $j$ to $i$
$A \approx 2000 \text{ N}$ is the interaction strength
$B \approx 0.08 \text{ m}$ is the range of the repulsion

An anisotropic version accounts for the fact that pedestrians react more strongly to people in their field of view:

$$\mathbf{F}_{ij}^{\text{rep}} = A \exp\!\left(\frac{r_{ij} - d_{ij}}{B}\right) \hat{\mathbf{n}}_{ij} \cdot w(\varphi_{ij})$$

with the anisotropy weight $w(\varphi) = \lambda + (1-\lambda)\frac{1+\cos\varphi}{2}$, where $\varphi_{ij}$ is the angle between $\mathbf{e}_i$ and $-\hat{\mathbf{n}}_{ij}$, and $\lambda \approx 0.5$.

4. Wall Repulsion

Walls exert a similar repulsive force:

$$\mathbf{F}_{iw}^{\text{wall}} = A_w \exp\!\left(\frac{r_i - d_{iw}}{B_w}\right) \hat{\mathbf{n}}_{iw}$$

where $d_{iw}$ is the shortest distance from pedestrian $i$ to wall $w$, and $\hat{\mathbf{n}}_{iw}$ is the normal pointing away from the wall. In dense crowds where pedestrians are pushed against walls, additional body compression and sliding friction forces are included:

$$\mathbf{F}_{iw}^{\text{contact}} = k \, g(r_i - d_{iw}) \, \hat{\mathbf{n}}_{iw} + \kappa \, g(r_i - d_{iw}) \, (\mathbf{v}_i \cdot \hat{\mathbf{t}}_{iw}) \, \hat{\mathbf{t}}_{iw}$$

where $g(x) = \max(0, x)$ activates only during physical contact,$k \approx 1.2 \times 10^5$ N/m is the body stiffness, and $\kappa \approx 2.4 \times 10^5$ kg/(m·s) is the friction coefficient.

5. Emergent Lane Formation

One of the most striking emergent phenomena in the Social Force Model is lane formation. When two groups of pedestrians walk in opposite directions in a corridor, they spontaneously organise into lanes -- without any explicit coordination.

The mechanism is intuitive: a pedestrian walking to the right is repelled less by other rightward-walking pedestrians (who are moving alongside) than by leftward-walking ones (who approach head-on with higher closing speed). Over time, same-direction pedestrians cluster. The lane order parameter is:

$$\Phi = \left\langle \frac{1}{N_i} \sum_{j \in \mathcal{N}_i} \text{sign}(v_j^x) \cdot \text{sign}(v_i^x) \right\rangle_i$$

where the average is over all pedestrians and their neighbours $\mathcal{N}_i$. When $\Phi \to 1$, perfect lanes have formed; when $\Phi \approx 0$, pedestrians are mixed.

6. Python: Corridor Flow with Lane Formation

We simulate two groups of pedestrians walking in opposite directions in a corridor. The simulation uses the Social Force Model with periodic boundary conditions along the corridor axis.

Social Force Model: Corridor Lane Formation

Python

script.py153 lines

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# --- Parameters ---
N_ped = 40          # total pedestrians (20 each direction)
L = 20.0            # corridor length
W = 4.0             # corridor width
dt = 0.05           # time step
T_total = 15.0      # total simulation time
tau = 0.5           # relaxation time
v0 = 1.3            # desired speed (m/s)
A = 2000.0          # repulsion strength (N)
B = 0.08            # repulsion range (m)
A_w = 2000.0        # wall repulsion strength
B_w = 0.08          # wall repulsion range
r_ped = 0.25        # pedestrian radius

# --- Initialise ---
x = np.random.uniform(0, L, N_ped)
y = np.random.uniform(r_ped + 0.1, W - r_ped - 0.1, N_ped)
vx = np.zeros(N_ped)
vy = np.zeros(N_ped)

# Group 1: walk right (+x), Group 2: walk left (-x)
direction = np.ones(N_ped)
direction[N_ped//2:] = -1.0

steps = int(T_total / dt)
history_x = np.zeros((steps + 1, N_ped))
history_y = np.zeros((steps + 1, N_ped))
history_x[0] = x.copy()
history_y[0] = y.copy()

for step in range(steps):
    fx = np.zeros(N_ped)
    fy = np.zeros(N_ped)

# Desired force
    fx += (v0 * direction - vx) / tau
    fy += (0 - vy) / tau

# Pedestrian repulsion
    for i in range(N_ped):
        for j in range(i+1, N_ped):
            dx = x[i] - x[j]
            # Periodic BC in x
            if dx > L/2: dx -= L
            if dx < -L/2: dx += L
            dy = y[i] - y[j]
            dist = np.sqrt(dx**2 + dy**2)
            if dist < 1e-6:
                dist = 1e-6
            rij = 2 * r_ped
            force_mag = A * np.exp((rij - dist) / B)
            nx_ij = dx / dist
            ny_ij = dy / dist
            fx[i] += force_mag * nx_ij
            fy[i] += force_mag * ny_ij
            fx[j] -= force_mag * nx_ij
            fy[j] -= force_mag * ny_ij

# Wall repulsion (top wall y=W, bottom wall y=0)
    for i in range(N_ped):
        # Bottom wall
        d_bottom = y[i]
        f_bottom = A_w * np.exp((r_ped - d_bottom) / B_w)
        fy[i] += f_bottom
        # Top wall
        d_top = W - y[i]
        f_top = A_w * np.exp((r_ped - d_top) / B_w)
        fy[i] -= f_top

# Euler integration
    vx += fx * dt
    vy += fy * dt
    # Speed limit
    speed = np.sqrt(vx**2 + vy**2)
    mask = speed > 2 * v0
    vx[mask] *= 2 * v0 / speed[mask]
    vy[mask] *= 2 * v0 / speed[mask]

x += vx * dt
    y += vy * dt
    # Periodic BC in x
    x = x % L
    # Clamp y
    y = np.clip(y, r_ped, W - r_ped)

history_x[step+1] = x.copy()
    history_y[step+1] = y.copy()

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

# Snapshot at end
ax = axes[0]
group1 = direction > 0
group2 = direction < 0
ax.scatter(x[group1], y[group1], c='#10b981', s=80, alpha=0.8, edgecolors='white', linewidths=0.5, label='Rightward', zorder=5)
ax.scatter(x[group2], y[group2], c='#f59e0b', s=80, alpha=0.8, edgecolors='white', linewidths=0.5, label='Leftward', zorder=5)
# Draw walls
ax.axhline(y=0, color='#64748b', linewidth=3)
ax.axhline(y=W, color='#64748b', linewidth=3)
ax.set_xlim(0, L)
ax.set_ylim(-0.5, W + 0.5)
ax.set_xlabel('x (m)', fontsize=12, color='white')
ax.set_ylabel('y (m)', fontsize=12, color='white')
ax.set_title(f'Final Positions (t={T_total}s)', fontsize=13, color='white')
ax.legend(fontsize=11, facecolor='#1e293b', edgecolor='#10b981', labelcolor='white')
ax.set_aspect('equal')
ax.tick_params(colors='white')

# Trajectories (last 5 seconds)
ax2 = axes[1]
start = max(0, steps - int(5.0/dt))
for i in range(N_ped):
    c = '#10b981' if direction[i] > 0 else '#f59e0b'
    ax2.plot(history_x[start:, i], history_y[start:, i], color=c, alpha=0.3, linewidth=0.8)
ax2.axhline(y=0, color='#64748b', linewidth=3)
ax2.axhline(y=W, color='#64748b', linewidth=3)
ax2.set_xlim(0, L)
ax2.set_ylim(-0.5, W + 0.5)
ax2.set_xlabel('x (m)', fontsize=12, color='white')
ax2.set_ylabel('y (m)', fontsize=12, color='white')
ax2.set_title('Trajectories (last 5s)', fontsize=13, color='white')
ax2.set_aspect('equal')
ax2.tick_params(colors='white')

fig.patch.set_facecolor('#0f172a')
for ax in axes:
    ax.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight', facecolor='#0f172a')
plt.show()

# Lane order parameter
phi_values = []
for i in range(N_ped):
    same_dir = 0
    count = 0
    for j in range(N_ped):
        if i == j: continue
        dx = abs(x[i] - x[j])
        if dx > L/2: dx = L - dx
        dy = abs(y[i] - y[j])
        if np.sqrt(dx**2 + dy**2) < 2.0:
            same_dir += 1 if direction[i] == direction[j] else -1
            count += 1
    if count > 0:
        phi_values.append(same_dir / count)
print(f"Lane order parameter Phi = {np.mean(phi_values):.3f} (1=perfect lanes, 0=mixed)")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

7. Evacuation Dynamics and the Faster-is-Slower Effect

In emergency evacuations, Helbing et al. discovered a counter-intuitive phenomenon: increasing the desired speed (panic parameter) above a critical value decreases the flow through a bottleneck. This is the faster-is-slower effect.

The mechanism: at high desired speeds, pedestrians push harder toward the exit, forming clogging arches. The flow rate $J$ as a function of desired speed $v^0$ is:

$$J(v^0) \propto \begin{cases} v^0 & v^0 < v^0_c \\ v^0 \exp\!\left[-\gamma (v^0 - v^0_c)^2\right] & v^0 \geq v^0_c \end{cases}$$

The panic parameter can be incorporated by modifying the relaxation time:

$$\tau_{\text{panic}} = \tau_0 (1 - n_{\text{panic}}) + \tau_{\min} \cdot n_{\text{panic}}$$

where $n_{\text{panic}} \in [0, 1]$ interpolates between normal ($\tau_0 \approx 0.5$ s) and panicked ($\tau_{\min} \approx 0.1$ s) behaviour. Higher panic also increases the desired speed: $v^0_{\text{panic}} = (1 + n_{\text{panic}}) v^0_{\text{normal}}$.

Evacuation: Faster-is-Slower Effect

Python

script.py128 lines

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(7)

# Evacuation simulation: room with one exit
N_ped = 30
room_W, room_H = 10.0, 10.0
exit_y = room_H / 2       # exit at right wall, centered
exit_width = 1.5
dt = 0.05
tau = 0.5
A, B = 2000.0, 0.08
A_w, B_w = 2000.0, 0.08
r_ped = 0.25

def simulate_evacuation(v0, T_max=30.0):
    x = np.random.uniform(1, room_W - 2, N_ped)
    y = np.random.uniform(1, room_H - 1, N_ped)
    vx = np.zeros(N_ped)
    vy = np.zeros(N_ped)
    active = np.ones(N_ped, dtype=bool)
    escaped_times = []

steps = int(T_max / dt)
    for step in range(steps):
        t = step * dt
        n_active = np.sum(active)
        if n_active == 0:
            break

# Goal: exit at (room_W, exit_y)
        goal_x = room_W
        goal_y = exit_y

fx = np.zeros(N_ped)
        fy = np.zeros(N_ped)

for i in range(N_ped):
            if not active[i]:
                continue
            # Desired force toward exit
            dx_g = goal_x - x[i]
            dy_g = goal_y - y[i]
            dist_g = np.sqrt(dx_g**2 + dy_g**2)
            if dist_g > 0.01:
                ex, ey = dx_g / dist_g, dy_g / dist_g
            else:
                ex, ey = 1.0, 0.0
            fx[i] += (v0 * ex - vx[i]) / tau
            fy[i] += (v0 * ey - vy[i]) / tau

# Repulsion from others
            for j in range(N_ped):
                if i == j or not active[j]:
                    continue
                ddx = x[i] - x[j]
                ddy = y[i] - y[j]
                dist = np.sqrt(ddx**2 + ddy**2)
                if dist < 0.01: dist = 0.01
                f_mag = A * np.exp((2*r_ped - dist) / B)
                fx[i] += f_mag * ddx / dist
                fy[i] += f_mag * ddy / dist

# Wall forces (4 walls, but exit opening on right)
            # Bottom (y=0)
            fy[i] += A_w * np.exp((r_ped - y[i]) / B_w)
            # Top (y=room_H)
            fy[i] -= A_w * np.exp((r_ped - (room_H - y[i])) / B_w)
            # Left (x=0)
            fx[i] += A_w * np.exp((r_ped - x[i]) / B_w)
            # Right wall (x=room_W) except at exit
            if y[i] < exit_y - exit_width/2 or y[i] > exit_y + exit_width/2:
                fx[i] -= A_w * np.exp((r_ped - (room_W - x[i])) / B_w)

# Update
        for i in range(N_ped):
            if not active[i]:
                continue
            vx[i] += fx[i] * dt
            vy[i] += fy[i] * dt
            speed = np.sqrt(vx[i]**2 + vy[i]**2)
            if speed > 3 * v0:
                vx[i] *= 3 * v0 / speed
                vy[i] *= 3 * v0 / speed
            x[i] += vx[i] * dt
            y[i] += vy[i] * dt
            y[i] = np.clip(y[i], r_ped, room_H - r_ped)
            x[i] = max(r_ped, x[i])
            # Check escape
            if x[i] >= room_W and abs(y[i] - exit_y) <= exit_width / 2:
                active[i] = False
                escaped_times.append(t)

return escaped_times

# Sweep desired speeds
v0_values = np.linspace(0.5, 5.0, 12)
mean_escape = []

for v0 in v0_values:
    times = simulate_evacuation(v0, T_max=40.0)
    if len(times) > 0:
        mean_escape.append(np.mean(times))
    else:
        mean_escape.append(40.0)

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(v0_values, mean_escape, 'o-', color='#10b981', linewidth=2.5, markersize=7, markerfacecolor='#34d399')
ax.set_xlabel('Desired Speed $v^0$ (m/s)', fontsize=13, color='white')
ax.set_ylabel('Mean Escape Time (s)', fontsize=13, color='white')
ax.set_title('Faster-is-Slower Effect in Evacuation', fontsize=14, color='white')
ax.tick_params(colors='white')
ax.grid(True, alpha=0.3)

# Mark critical speed
min_idx = np.argmin(mean_escape)
ax.axvline(x=v0_values[min_idx], color='#f59e0b', linestyle='--', linewidth=1.5,
           label=f'Optimal $v^0 \\approx {v0_values[min_idx]:.1f}$ m/s')
ax.legend(fontsize=11, facecolor='#1e293b', edgecolor='#10b981', labelcolor='white')

fig.patch.set_facecolor('#0f172a')
ax.set_facecolor('#0f172a')
plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight', facecolor='#0f172a')
plt.show()
print(f"Optimal desired speed: {v0_values[min_idx]:.2f} m/s")
print(f"Min mean escape time: {mean_escape[min_idx]:.2f} s")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Bid-Rent Gradient Part V: Stochastic Dynamics →