Skip to content
BranchesTags

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

1 Commit
.gitignore
.gitignore
 
 
README.md
README.md
 
 
oba_assignment.py
oba_assignment.py
 
 

Repository files navigation

User Equilibrium Traffic Assignment — Chicago-Sketch

Wardrop User Equilibrium (UE) traffic assignment for the Chicago-Sketch benchmark network. Solver: Bi-Conjugate Frank-Wolfe (BFW) with exact line search. Convergence criterion: Relative Gap < 1×10⁻⁴.

Table of Contents

  1. Problem Formulation
  2. Algorithm Design
  3. Implementation
  4. Dataset
  5. Dependencies & Installation
  6. Usage
  7. Configuration
  8. Outputs
  9. Performance
  10. References

1. Problem Formulation

The UE assignment problem seeks a feasible flow pattern x on network G(N, A) such that no traveller can reduce their travel time by unilaterally changing routes (Wardrop's first principle).

Beckmann et al. (1956) showed this is equivalent to minimising the convex objective:

$$\min_{\mathbf{x}} ; Z(\mathbf{x}) = \sum_{a \in A} \int_0^{x_a} t_a(u), du$$

subject to:

$$\sum_{k \in K_{rs}} f_k^{rs} = q_{rs}, \quad \forall r,s$$ $$f_k^{rs} \geq 0, \quad \forall k, r, s$$ $$x_a = \sum_{r,s} \sum_k f_k^{rs} \delta_{ak}^{rs}, \quad \forall a$$

where $f_k^{rs}$ is the flow on path $k$ between OD pair $(r,s)$, $q_{rs}$ is the demand, and $\delta_{ak}^{rs} = 1$ if link $a$ is on path $k$.

2. Algorithm Design

2.1 Beckmann Objective

The gradient of the Beckmann function with respect to link flows equals the vector of current link travel times:

$$\nabla Z(\mathbf{x}) = \mathbf{t}(\mathbf{x})$$

This makes every Frank-Wolfe–type direction automatically descent-compatible via the dot-product check $\mathbf{t}^T \mathbf{d} &lt; 0$.

2.2 Generalized Link Cost

Links are costed using a generalized cost combining BPR travel time, distance, and toll:

$$c_a(x_a) = t_a^0 \left[1 + b_a \left(\frac{x_a}{C_a}\right)^{n_a}\right] + w_d \cdot l_a + w_\tau \cdot \tau_a$$

Parameter Symbol Value
Free-flow travel time $t_a^0$ from network file (min)
BPR alpha coefficient $b_a$ from network file
BPR beta exponent $n_a$ from network file
Link capacity $C_a$ from network file (veh/hr)
Link length $l_a$ from network file (mi)
Toll $\tau_a$ from network file (cents)
Distance weight $w_d$ 0.04 min/mi
Toll weight $w_\tau$ 0.02 min/cent

The fixed component $g_a = w_d l_a + w_\tau \tau_a$ is pre-computed once and reused every iteration. This is critical for connector links (fftt = 0) whose entire cost derives from distance.

The BPR derivative is:

$$\frac{\partial c_a}{\partial x_a} = t_a^0 \cdot b_a \cdot n_a \left(\frac{x_a}{C_a}\right)^{n_a - 1} \frac{1}{C_a}$$

2.3 Frank-Wolfe Framework

The standard Frank-Wolfe (FW) iteration is:

  1. Solve the linearised sub-problem: $\mathbf{y}^k = \arg\min_{\mathbf{x} \text{ feasible}} \mathbf{t}(\mathbf{x}^k)^T \mathbf{x}$ → All-or-Nothing (AON) assignment with costs $\mathbf{t}^k$
  2. Set direction $\mathbf{d}^k = \mathbf{y}^k - \mathbf{x}^k$
  3. Line search: $\alpha^k = \arg\min_{\alpha \in [0,1]} Z(\mathbf{x}^k + \alpha \mathbf{d}^k)$
  4. Update: $\mathbf{x}^{k+1} = \mathbf{x}^k + \alpha^k \mathbf{d}^k$

FW converges at O(1/k), which is slow near the optimum. Conjugate and bi-conjugate extensions improve this to O(1/k²).

2.4 Bi-Conjugate Direction

The Bi-Conjugate Frank-Wolfe (BFW) method (Mitradjieva & Lindberg 2013, Liu & Nie 2010) blends the current FW target with two previous iterate vectors to construct a conjugate-like search direction:

$$\mathbf{d}^k = \beta \underbrace{(\mathbf{y}^k - \mathbf{x}^k)}_{\mathbf{f}} + \gamma \underbrace{(\mathbf{x}^{k-2} - \mathbf{x}^k)}_{\mathbf{s}_2} + (1-\beta-\gamma) \underbrace{(\mathbf{x}^{k-1} - \mathbf{x}^k)}_{\mathbf{s}_1}$$

The coefficients $(\beta, \gamma)$ are found by minimising the second-order Taylor expansion of $Z$ along $\mathbf{d}^k$, leading to the 2×2 linear system:

$$\begin{bmatrix} A & B \ B & C \end{bmatrix} \begin{bmatrix} \beta \ \gamma \end{bmatrix} = \begin{bmatrix} -c_1 \ -c_2 \end{bmatrix}$$

where, with $\Delta\mathbf{f} = \mathbf{f} - \mathbf{s}_1$ and $\Delta\mathbf{s} = \mathbf{s}_2 - \mathbf{s}_1$:

Symbol Definition
$A$ $\sum_a t'_a ,(\Delta f_a)^2$
$B$ $\sum_a t'_a ,\Delta f_a ,\Delta s_a$
$C$ $\sum_a t'_a ,(\Delta s_a)^2$
$c_1$ $\sum_a (t_a + t'a s{1,a}) ,\Delta f_a$
$c_2$ $\sum_a (t_a + t'a s{1,a}) ,\Delta s_a$

The solution is:

$$\beta = \frac{B c_2 - C c_1}{\det}, \quad \gamma = \frac{B c_1 - A c_2}{\det}, \quad \det = AC - B^2$$

Projection onto the 2-simplex: $\beta \geq 0$, $\gamma \geq 0$, $\beta + \gamma \leq 1$. If $\beta + \gamma &gt; 1$ normalize: $\beta \leftarrow \beta/(\beta+\gamma)$, $\gamma \leftarrow \gamma/(\beta+\gamma)$.

Iteration 1 (no history): pure FW direction ($\mathbf{d}^1 = \mathbf{f}$). Iteration 2 (one previous iterate): CFW formula (1×1 system, $\beta$ only). Iteration 3+: full BFW.

Descent safeguard: if $\mathbf{t}^T \mathbf{d} \geq 0$, fall back to FW direction.

2.5 Exact Line Search

The optimal step size minimises $Z(\mathbf{x} + \alpha \mathbf{d})$ over $\alpha \in [0,1]$:

$$\frac{dZ}{d\alpha} = \sum_a d_a \cdot c_a(x_a + \alpha d_a) = 0$$

Since $c_a$ are strictly increasing (BPR), $dZ/d\alpha$ is strictly increasing in $\alpha$, guaranteeing a unique root. Bisection on this slope function converges in 50 steps to machine precision:

slope(α) = Σ_a  d_a · [t_a^0·(1 + b_a·((x_a+αd_a)/C_a)^n_a) + g_a]

Special cases: if $\text{slope}(0) \geq 0$ return $\alpha = 0$; if $\text{slope}(1) \leq 0$ return $\alpha = 1$.

2.6 Convergence Criterion

Relative Gap (RelGap):

$$\text{RelGap}^k = \frac{\text{TSTT}^k - \text{SPTT}^k}{\text{SPTT}^k}$$

where:

  • $\text{TSTT}^k = \sum_a x_a^k \cdot c_a(x_a^k)$ — Total System Travel Time at current flows
  • $\text{SPTT}^k = \sum_{r,s} q_{rs} \cdot d_{rs}^k$ — Shortest Path Total Travel Time (sum of minimum OD path costs at current costs)

Both are evaluated with the same post-update cost vector $\mathbf{c}(\mathbf{x}^{k+1})$, ensuring consistency. The AON pass at the end of each iteration provides both $\mathbf{y}^{k+1}$ (for the next direction) and $\text{SPTT}^{k+1}$ simultaneously.

Convergence is declared when $\text{RelGap} &lt; 10^{-4}$.

Theoretical bound: At Wardrop equilibrium, $\text{RelGap} = 0$.

2.7 All-or-Nothing Assignment

For each origin zone $r$, a shortest-path tree is computed with Dijkstra's algorithm under current generalised costs. Demand $q_{rs}$ is loaded entirely onto the minimum-cost path to each destination $s$.

Subtree flow propagation (O(N) per origin, avoids per-OD-pair path tracing):

  1. Initialise sub[v] = Σ_s q_{rs} for destination nodes $v = s$
  2. Process nodes in decreasing distance order (valid since all link costs > 0 → child dist > parent dist)
  3. For each node $v$: y[link(pred(v)→v)] += sub[v]; sub[pred(v)] += sub[v]

This reduces complexity from O(path_length × |destinations|) to O(N + |destinations|) per origin.

3. Implementation

3.1 Module Structure

oba_assignment.py
├── parse_network()          — parse ChicagoSketch_net.tntp
├── parse_trips()            — parse ChicagoSketch_trips.tntp
├── parse_nodes()            — parse ChicagoSketch_node.tntp
├── Network                  — holds arrays: tail, head, cap, fftt, cost, flow, …
│   ├── update_costs()       — BPR + gen_fixed → self.cost
│   ├── cost_deriv()         — ∂c_a/∂x_a (BPR derivative)
│   └── tstt()               — Σ x_a · c_a(x_a)
├── BFW                      — solver
│   ├── _setup_scipy()       — build CSR matrix + edge_to_lid lookup
│   ├── dijkstra()           — single-source Dijkstra (Python, fallback)
│   ├── aon()                — dispatch to _aon_scipy or _aon_python
│   ├── _aon_scipy()         — multi-source scipy Dijkstra + subtree flow
│   ├── _aon_python()        — pure-Python fallback
│   ├── line_search()        — bisection on dZ/dα = 0
│   ├── bfw_direction()      — compute β, γ; construct BFW/CFW/FW direction
│   └── run()                — main iteration loop; returns (gap, history)
├── sp_trees_final()         — build final SP trees for post-processing
├── save_paths()             — write output/paths.csv
├── save_turning_movements() — write output/turning_movements.csv
├── save_network_shp()       — write output/network.*  (ESRI Shapefile)
├── save_nodes_shp()         — write output/nodes.*    (ESRI Shapefile)
└── main()

3.2 Data Structures

All link-indexed arrays are NumPy float64/int32 arrays of length M = 2950:

Array Content
net.tail, net.head tail/head node of each link (int32)
net.cap, net.fftt capacity (veh/hr), free-flow travel time (min)
net.b_coef, net.power BPR α and β parameters
net.length, net.toll link length (mi), toll (cents)
net.gen_fixed pre-computed $w_d l_a + w_\tau \tau_a$
net.cost current generalised cost (updated each iteration)
net.flow current link flow (veh/hr)

Adjacency is stored as Python lists net.out_adj[u] = [(v, lid), …] for Dijkstra traversal.

The CSR sparse matrix BFW.csr_mat (shape N×N, nnz = 2950) is built once during setup. Each iteration, only csr_mat.data[:] is updated (O(M) write), avoiding matrix reconstruction.

The edge-to-link lookup BFW.edge_to_lid_arr is a (N×N) int32 array. edge_to_lid_arr[u, v] gives the link ID for directed edge (u→v), enabling O(1) predecessor-node-to-link conversion during flow propagation.

3.3 AON Acceleration

The AON sub-problem is the computational bottleneck. Three optimisations are applied:

1. scipy multi-source Dijkstra scipy.sparse.csgraph.dijkstra(csr, indices=all_zones) runs all 387 Dijkstra searches in a single C-level call, returning dist_mat (R×N) and pred_mat (R×N). This avoids Python interpreter overhead for each search.

2. Distance-sort subtree propagation Since all generalised costs are positive, the shortest-path tree has the property that dist[child] > dist[parent] strictly. Sorting reachable nodes by decreasing distance gives a valid bottom-up processing order without explicit BFS/DFS tree construction.

3. Pre-allocated reuse buffer A single sub array (length N) is zeroed and reused for each origin, avoiding N×R allocations per iteration.

3.4 Multithreading Strategy

The multithreading infrastructure supports parallel Dijkstra for large networks:

  • Source zones are pre-split into N_THREADS equal chunks
  • A persistent ThreadPoolExecutor submits one scipy.dijkstra task per chunk
  • scipy releases the Python GIL during its C computation → true parallel execution
  • The subtree propagation runs single-threaded on the main thread

For Chicago-Sketch (N = 933): N_THREADS = 1 (single call, no pool). The bottleneck is the Python subtree loop (~0.15 s/iter), which is GIL-bound. Spawning idle threads competes for the GIL and degrades performance by 3–4×.

For larger networks (N ≥ 10K): set N_THREADS = os.cpu_count(). When per-source Dijkstra time dominates (>>0.1 s), parallel chunking yields near-linear speedup on the C-Dijkstra phase.

3.5 Post-processing

Paths (output/paths.csv) Final SP trees are reconstructed from sp_trees_final(), which calls scipy dijkstra once more at converged costs. For each OD pair with positive demand, the predecessor-link array is traced back to reconstruct the node sequence.

Turning movements (output/turning_movements.csv) For each origin, a post-order DFS traversal of the SP tree computes subtree_flow[v] = total demand destined for any descendant of v. At each interior node j, the turn from the incoming link pred_link[j] to each outgoing link pred_link[child] carries subtree_flow[child].

turn_vol[(lid_in, lid_out)] += subtree_flow[child]
  where  lid_in  = pred_link[junction_node]
         lid_out = pred_link[child_node]

Only turns with volume ≥ 0.01 veh are written.

Shapefiles Written with pyshp. Links are polylines connecting tail/head node coordinates. Both network.* and nodes.* include a .prj file with the NAD83 / Illinois State Plane East (FIPS 1201, US survey feet) WKT projection string.

4. Dataset

Chicago-Sketch — a widely-used benchmark for UE assignment algorithms.

Property Value
Nodes 933
Zone centroids 387
Links 2,950
Connector links 774 (387 outgoing + 387 incoming)
OD pairs (positive demand) ~150,000
Coordinate system NAD83 / Illinois State Plane East, US survey feet
Reference UE TSTT ~1.893 × 10⁷ veh·min
Reference max V/C ~2.43

Source files (TNTP format):

File Content
Chicago-Sketch/ChicagoSketch_net.tntp Link geometry and BPR parameters
Chicago-Sketch/ChicagoSketch_trips.tntp OD demand matrix
Chicago-Sketch/ChicagoSketch_node.tntp Node coordinates (State Plane feet)

5. Dependencies & Installation

Python 3.8+ required.

pip install numpy scipy pyshp matplotlib
Package Purpose
numpy All array operations, BPR evaluation, BFW direction
scipy scipy.sparse.csgraph.dijkstra — C-level multi-source Dijkstra
pyshp (shapefile) Writing ESRI Shapefiles
matplotlib Convergence plot (output/convergence.png) — optional

If matplotlib is not installed the plot step is silently skipped; all other outputs are still produced.

6. Usage

Place the three TNTP data files in a subdirectory named Chicago-Sketch/ next to the script, then run:

python oba_assignment.py

Output is written to output/ (created automatically).

Expected console output:

Loading network data...
  Links=2950  Nodes=933  Zones=387

Running BFW assignment...
  AON mode: single-thread scipy Dijkstra
  [init] done.  TSTT=1.1148e+08

   Iter      RelGap            TSTT            SPTT        α   t/iter    total
  ───────────────────────────────────────────────────────────────────────────
      1   3.462e-01   2.5172e+07   1.8699e+07   0.4763   0.30s    0.9s
    ...
    164   9.423e-05   1.8935e+07   1.8933e+07   0.3937   0.30s   53.8s

  Converged at iteration 164: RelGap=9.4227e-05

── Summary ──────────────────────────────────────────────────────────
  Iterations         : 164
  Final Relative Gap : 9.4227e-05
  TSTT               : 1.893513e+07  veh·min
  V/C  max=2.434  mean=0.416
  Wall-clock time    : 54.0s

7. Configuration

Top-level constants in oba_assignment.py:

Constant Default Description
CONV 1e-4 Convergence threshold (Relative Gap)
MAX_ITER 500 Maximum BFW iterations
DIST_WEIGHT 0.04 Distance penalty (min/mi) in generalised cost
TOLL_WEIGHT 0.02 Toll penalty (min/cent) in generalised cost
N_THREADS 1 Parallel Dijkstra workers; set to os.cpu_count() for N ≥ 10K networks
DATA_DIR ./Chicago-Sketch Path to TNTP input files
OUT_DIR ./output Output directory

8. Outputs

output/
├── paths.csv               SP path for every OD pair
├── turning_movements.csv   Intersection turning movement volumes
├── convergence.csv         Per-iteration convergence history
├── convergence.png         Two-panel convergence plot
├── network.shp / .dbf / .shx / .prj    Link Shapefile (polylines)
└── nodes.shp   / .dbf / .shx / .prj    Node Shapefile (points)

paths.csv

Column Content
origin Origin zone ID
destination Destination zone ID
demand OD demand (veh/hr)
path_nodes Space-separated node sequence (origin → destination)

turning_movements.csv

Column Content
junction_node Intersection node ID
from_node Upstream node of the incoming link
junction_in Same as junction_node
in_link_id Incoming link index
to_node Downstream node of the outgoing link
junction_out Same as junction_node
out_link_id Outgoing link index
turn_volume_veh Turning movement volume (veh/hr)

convergence.csv

Column Content
iter Iteration number
rel_gap Relative Gap
tstt Total System Travel Time (veh·min)
sptt Shortest Path Total Travel Time (veh·min)
alpha Line search step size
iter_time_s Wall time for this iteration (s)
cumulative_s Cumulative wall time (s)

network Shapefile attributes

Field Type Content
LINK_ID Integer Link index (0-based)
TAIL / HEAD Integer Tail/head node ID
CAPACITY Float Capacity (veh/hr)
LENGTH_MI Float Length (miles)
FFTT_MIN Float Free-flow travel time (min)
B_COEF Float BPR alpha
POWER Float BPR beta
TOLL Float Toll (cents)
LINK_TYPE Integer Link type code
FLOW Float Assigned flow (veh/hr)
VC_RATIO Float Volume-to-capacity ratio
TRAV_TIME Float BPR travel time at assigned flow (min)
GEN_COST Float Generalised cost at assigned flow (min)

Projection: NAD83 / Illinois State Plane East (FIPS 1201), US survey feet.

9. Performance

Benchmark results on Chicago-Sketch (Intel/AMD desktop, Python 3.11, single thread):

Metric Value
Iterations to RelGap < 1×10⁻⁴ 164
Wall-clock time ~54 s
Time per iteration ~0.30 s
Final RelGap 9.42×10⁻⁵
Final TSTT 1.8935×10⁷ veh·min
Max V/C 2.43

Convergence profile (selected iterations):

Iter RelGap TSTT (veh·min)
1 3.46×10⁻¹ 2.517×10⁷
5 3.38×10⁻² 1.969×10⁷
10 8.28×10⁻³ 1.908×10⁷
20 2.34×10⁻³ 1.896×10⁷
50 5.04×10⁻⁴ 1.894×10⁷
100 3.37×10⁻⁴ 1.894×10⁷
164 9.42×10⁻⁵ 1.894×10⁷

The convergence plot (output/convergence.png) shows both RelGap (log scale) and TSTT versus iteration number.

10. References

  1. Wardrop, J.G. (1952). Some theoretical aspects of road traffic research. Proceedings of the Institution of Civil Engineers, 1(3), 325–362.

  2. Beckmann, M., McGuire, C.B., & Winsten, C.B. (1956). Studies in the Economics of Transportation. Yale University Press.

  3. Bureau of Public Roads (1964). Traffic Assignment Manual. U.S. Department of Commerce. (BPR volume-delay function: $t = t^0[1 + 0.15(v/c)^4]$)

  4. Frank, M. & Wolfe, P. (1956). An algorithm for quadratic programming. Naval Research Logistics Quarterly, 3(1–2), 95–110.

  5. Liu, Z. & Nie, Y. (2010). A new algorithm for traffic assignment: Finding shortest paths and equilibrating with two auxiliary networks. Transportation Science.

  6. Mitradjieva, M. & Lindberg, P.O. (2013). The stiff is moving — conjugate direction Frank-Wolfe methods with applications to traffic assignment. Transportation Science, 47(2), 280–293.

  7. Bar-Gera, H. (2002). Origin-based algorithm for the traffic assignment problem. Transportation Science, 36(4), 398–417.

  8. Sheffi, Y. (1985). Urban Transportation Networks: Equilibrium Analysis with Mathematical Programming Methods. Prentice-Hall.

  9. Transportation Networks for Research (TNTP format & datasets): https://github.com/bstabler/TransportationNetworks

About

oba method for traffic assignment

Resources

Readme
Activity

Stars

0 stars

Watchers

0 watching

Forks

0 forks
Report repository

Releases

Packages

Contributors

Languages