Added 18/09/2025
Transportation / evacuation
siouxfalls_network
Dimension
{
"x": 40,
"y": 3040,
"F": 1,
"G": 80,
"H": 4,
"f": 1,
"g": 3040,
"h": 960
}Solution
{
"optimality": "infeasible",
"F": 0,
"H": [-2000,-9000,-7000,-2000],
"f": 0
}\subsection{siouxfalls\_network}
\label{sec:siouxfalls_network}
% Description
\paragraph{Description}
This problem is based on the work by Ng~\cite{Ng2010}.
The case study uses the Sioux Falls transportation network with designated origin and shelter nodes.
% Model Formulation
\paragraph{Model Formulation}
Let $G=(V,A)$ be the transportation network with node set $V$ and arc set $A$.
The Bureau of Public Roads (BPR) function defines arc travel times:
\[
t_a(v_a) = t_a^0\left(1+\alpha_a\left(\frac{v_a}{c_a}\right)^{\beta_a}\right),
\quad a\in A,
\]
where $t_a^0$ is the free-flow travel time, $c_a$ the capacity, and $\alpha_a,\beta_a$ calibration parameters.
% Upper-level
\paragraph{Upper-level (Central Planner)}
The central authority assigns evacuees from origins $r\in R \subseteq V$ to shelters $s\in S \subseteq V$:
\[
\min_{q_{rs}} \;\; \sum_{a\in A} v_a t_a(v_a)
\]
subject to
\[
\sum_{s\in S} q_{rs} = O_r, \quad \forall r\in R,
\]
\[
0 \le q_{rs} \le \overline{q}_{rs}, \quad \forall r\in R, \; s\in S,
\]
where $O_r$ is the population at origin $r$ and $\overline{q}_{rs}$ the shelter capacity bound.
% Lower-level
\paragraph{Lower-level (User Equilibrium Traffic Assignment)}
Given $q_{rs}$, evacuees choose shortest paths selfishly:
\[
\min_{f_{rs}^k} \;\; \sum_{a\in A} \int_0^{v_a} t_a(x) \, dx
\]
subject to
\[
\sum_{k} f_{rs}^k = q_{rs}, \quad \forall r\in R, \; s\in S,
\]
\[
v_a = \sum_{r,s}\sum_{k} \delta_{ars}^k f_{rs}^k, \quad \forall a\in A,
\]
\[
f_{rs}^k \ge 0, \quad \forall r,s,k,
\]
where $f_{rs}^k$ is the flow of evacuees from origin $r$ to shelter $s$ on path $k$, and $\delta_{ars}^k=1$ if arc $a$ is on path $k$ connecting $(r,s)$ and $0$ otherwise.
% Data
\paragraph{Data}
The Sioux Falls network has 24 nodes and 76 arcs. Origins are nodes 14, 15, 22, and 23 (with $O_r=2000, 9000, 7000, 2000$).
Shelters are nodes 4, 5, 6, 8, 9, 10, 11, 16, 17, and 18, with capacities
5000, 4000, 6000, 5000, 4000, 4000, 4000, 4000, 1000, 5000 respectively.
import os.path
import json
import numpy as np
from numpy.typing import NDArray
"""
A bilevel model for optimal shelter assignment in emergency evacuations.
Ng, Park, and Waller (2010) - A Hybrid Bilevel Model for the Optimal Shelter Assignment in Emergency Evacuations.
"""
# Properties
name: str = "siouxfalls_network"
category: str = "transportation"
subcategory: str = "evacuation"
datasets: list = ['siouxfalls.json']
paths: list = [
os.path.join('bolib3', 'data', 'network', 'siouxfalls.json'),
]
# Upper-level objective: minimize total travel time
def F(x: NDArray, y: NDArray, data: dict) -> float:
"""
Upper-level objective function (minimize total evacuation travel time).
"""
num_arcs: int = data['num_arcs']
num_comm: int = data['num_commodities']
# Reshape flows y into (num_comm x num_arcs)
y2d = np.reshape(y, (num_comm, num_arcs))
# Total flow on each arc
v = np.sum(y2d, axis=0)
cost = data['cost']
return np.dot(cost, v)
# Upper-level constraints: q >= 0, q <= capacity
def G(x: NDArray, y: NDArray, data: dict) -> NDArray:
"""
Upper-level inequality constraints (bounds on assignment q).
Ensures 0 <= q_rs <= C_s for each origin-shelter pair.
"""
q_upper: NDArray = data['q_upper']
return np.hstack([x, q_upper - x])
# Upper-level equalities: sum_s q_rs = O_r
def H(x: NDArray, y: NDArray, data: dict) -> NDArray:
"""
Upper-level equality constraints (evacuee assignment).
For each origin r: sum_s q_rs = O_r.
"""
origins = data['origins']
shelters = data['shelters']
num_orig = len(origins)
num_shel = len(shelters)
x2d = np.reshape(x, (num_orig, num_shel))
O = np.array(data['O'])
return np.array([np.sum(x2d[i, :]) - O[i] for i in range(num_orig)])
# Lower-level objective: total travel time (user equilibrium)
def f(x: NDArray, y: NDArray, data: dict) -> float:
"""
Lower-level objective function (total travel time under user equilibrium).
"""
num_arcs: int = data['num_arcs']
num_comm: int = data['num_commodities']
y2d = np.reshape(y, (num_comm, num_arcs))
v = np.sum(y2d, axis=0)
cost = data['cost']
return np.dot(cost, v)
# Lower-level flow nonnegativity
def g(x: NDArray, y: NDArray, data: dict) -> NDArray:
"""
Lower-level inequality constraints (flow >= 0).
"""
return y
# Lower-level flow conservation at each node for each commodity
def h(x: NDArray, y: NDArray, data: dict) -> NDArray:
"""
Lower-level equality constraints (flow conservation).
For each commodity j (origin r to shelter s):
outflow - inflow = +q at origin r, -q at shelter s.
"""
num_nodes = data['num_nodes']
num_arcs = data['num_arcs']
num_comm = data['num_commodities']
y2d = np.reshape(y, (num_comm, num_arcs))
inflow = np.array(data['inflow'])
outflow = np.array(data['outflow'])
# Net flow at each node for each commodity
netflow = outflow.dot(y2d.T) - inflow.dot(y2d.T) # shape: (num_nodes x num_comm)
origins = data['origins']
shelters = data['shelters']
idx = 0
for i, r in enumerate(origins):
for j, s in enumerate(shelters):
netflow[r, idx] -= x[idx] # origin r: outflow - inflow - q = 0
netflow[s, idx] += x[idx] # shelter s: outflow - inflow + q = 0
idx += 1
return netflow.flatten()
# Read Sioux Falls network and evacuation data
def read_data(filepath=paths[0]) -> dict:
"""
Reads the Sioux Falls network JSON data and sets up evacuation parameters.
"""
with open(filepath, 'r') as f:
network = json.load(f)
num_nodes = network['num_nodes']
num_arcs = network['num_arcs']
inflow = np.array(network['inflow'])
outflow = np.array(network['outflow'])
cost = np.array(network['free_flow_time'], dtype=float)
# Origins (0-based nodes 13,14,21,22) with evacuees 2000,9000,7000,2000
origins = [13, 14, 21, 22]
O = [2000.0, 9000.0, 7000.0, 2000.0]
# Shelter nodes (0-based 3,4,5,7,8,9,10,15,16,17) and their capacities
shelters = [3, 4, 5, 7, 8, 9, 10, 15, 16, 17]
shelter_capacities = [5000.0, 4000.0, 6000.0, 5000.0, 4000.0,
4000.0, 4000.0, 4000.0, 1000.0, 5000.0]
# Construct list of commodities (origin-shelter pairs)
origin_nodes = []
shelter_nodes = []
for r in origins:
for s in shelters:
origin_nodes.append(r)
shelter_nodes.append(s)
num_comm = len(origin_nodes)
# Upper bounds on q (capacity of shelter for each commodity)
q_upper = []
for r in origins:
for cap in shelter_capacities:
q_upper.append(cap)
return {
'num_nodes': num_nodes,
'num_arcs': num_arcs,
'num_commodities': num_comm,
'inflow': inflow,
'outflow': outflow,
'cost': cost,
'origins': origins,
'O': O,
'shelters': shelters,
'shelter_capacities': shelter_capacities,
'comm_origin': origin_nodes,
'comm_dest': shelter_nodes,
'q_upper': np.array(q_upper, dtype=float),
}
def dimension(key: str = '', data=None) -> int:
"""
Return dimensions for variables/functions.
"""
num_origins = len(data['origins'])
num_shelters = len(data['shelters'])
num_arcs = data['num_arcs']
num_comm = data['num_commodities']
dims = {
'x': num_origins*num_shelters,
'y': num_arcs*num_comm,
'F': 1,
'G': 2*num_origins*num_shelters,
'H': num_origins,
'f': 1,
'g': num_arcs*num_comm,
'h': data['num_nodes']*num_comm,
}
return dims.get(key, dims)