pm-td — Peeta & Mahmassani’s (1995) time-dependent UE and SO

What. pm-td is the methodology paper of simulation-based dynamic traffic assignment: time-dependent, EXPERIENCED-time system-optimal AND user- equilibrium route choice, solved by method-of-successive-averages (MSA) iterations around the repo’s own CTM/LTM loading (playing the role of the paper’s mesoscopic simulator). One machinery, two objectives — TD-UE assigns to least EXPERIENCED time, TD-SO to least time-dependent MARGINAL time — this notebook covers both (pm-td-ue and pm-td-so).

Why it is in the benchmark. It is the first ITERATIVE, simulation-based time-dependent route equilibrium: multi-origin/multi-destination, path-based, FIFO-enforced by the loader, with no closed form. Its headline result — the system optimum strictly beats user equilibrium’s total travel time — is made executable here. See the model compendium (Peeta & Mahmassani 1995) and docs/design/adr-031-peeta-mahmassani.md (P1).

Scope. This notebook solves TD-UE on a symmetric two-route diamond (an exact 50/50 anchor) and both TD-UE/TD-SO on an SO != UE wedge, certifying every result. The MSA solvers are explicitly NON-certified research code — the certifier is the arbiter, never the solver’s own claim.

Canon. [peeta1995system], docs/REFERENCES.md / docs/references.bib.

How this notebook is graded

A notebook never claims a number it does not compute in that cell. Every scored quantity below is recomputed live by the P1 TDTAEvaluator from the emitted per-path departure plan alone — the harness reruns its OWN per-path loading and recomputes the experienced-time route-swap residual (tdue_gap) and the LP-bound gap (so_bound_gap) against a freshly-solved lp-so-dta- style optimum, never trusting the MSA solver’s self-reported iterate (README, Certified, not self-reported).

# Setup. `pm-td` is a core model: a plain `pip install -e .` suffices — no
# optional extra, so no guard cell. The inline backend is Agg-based (headless
# CI renders into the notebook); NEVER matplotlib.use("Agg") in-kernel — it
# silently suppresses inline figure capture.
%matplotlib inline
import numpy as np

from tabench import (
    PathLoader,
    TDPathFlows,
    TDTAEvaluator,
    pm_diamond_scenario,
    pm_wedge_scenario,
    solve_td_so,
    solve_td_ue,
    viz,
)

pm-td-ue: an exact 50/50 anchor

The symmetric two-route diamond: origin 1 to destination 2 by two byte-identical routes (fast feeder + Q=1 bottleneck each). By symmetry the exact TD-UE is the 50/50 split every interval; an all-on-one control queues one bottleneck while its twin sits idle — a hand-computable positive gap. TDTAScenario is frozen and content-hashed (P2).

diamond = pm_diamond_scenario()
print(f"scenario      : {diamond.name}")
print(f"content hash  : {diamond.content_hash()[:16]}…")
print(f"paths         : {[p.links for p in diamond.paths]}")
print(f"grid          : dt={diamond.grid.dt}, n_steps={diamond.grid.n_steps}")
print("task          : time-dependent user equilibrium, experienced travel time")

scenario      : pm-diamond-ctm
content hash  : 1ed2c58061543b1b…
paths         : [(0, 1), (2, 3)]
grid          : dt=1.0, n_steps=16
task          : time-dependent user equilibrium, experienced travel time

ue_flows = solve_td_ue(diamond, iters=40)
evaluator = TDTAEvaluator(diamond)
metrics = evaluator.certify(ue_flows)
print(f"feasible      : {metrics['feasible']:.0f}")
print(f"tdue_gap      : {metrics['tdue_gap']:.3e}")
print(f"tstt          : {metrics['tstt']:.3f}")
assert metrics["feasible"] == 1.0
assert metrics["tdue_gap"] < 1e-6
assert np.isclose(metrics["tstt"], 14.0, atol=1e-2)

# The 50/50 split, recomputed from the emitted departures directly.
per_path = ue_flows.departures.sum(axis=1)
print(f"departures per path : {per_path}")
np.testing.assert_allclose(per_path, [2.0, 2.0], atol=1e-6)

# DISTINCTIVE: an all-on-one control certifies a hand-computable POSITIVE gap
# (0.75) -- built here from the SAME total demand, forced onto path 0 alone.
all_on_one = np.zeros_like(ue_flows.departures)
all_on_one[0] = ue_flows.departures.sum(axis=0)
gap_all_on_one = evaluator.certify(TDPathFlows(diamond.content_hash(), all_on_one))["tdue_gap"]
print(f"all-on-one control tdue_gap : {gap_all_on_one:.4f}")
assert np.isclose(gap_all_on_one, 0.75, atol=1e-3)

feasible      : 1
tdue_gap      : 0.000e+00
tstt          : 14.000
departures per path : [2. 2.]
all-on-one control tdue_gap : 0.7500

pm-td-ue vs pm-td-so: SO strictly beats UE (the paper’s headline)

The wedge: a FAST but capacitated route (Q=1 bottleneck) versus a SLOWER uncongested route. Selfish UE overloads the fast route until its queue delay matches the slow route’s extra free-flow time; the system optimum diverts more traffic to the slow route earlier, cutting total queueing.

wedge = pm_wedge_scenario()
print(f"scenario      : {wedge.name}")
print(f"content hash  : {wedge.content_hash()[:16]}…")

ue_wedge = solve_td_ue(wedge, iters=40)
so_wedge = solve_td_so(wedge, iters=40)
ev_wedge = TDTAEvaluator(wedge)
m_ue = ev_wedge.certify(ue_wedge)
m_so = ev_wedge.certify(so_wedge)
print(f"UE: tdue_gap={m_ue['tdue_gap']:.4f}  tstt={m_ue['tstt']:.3f}")
print(f"SO: so_bound_gap={m_so['so_bound_gap']:.3e}  tstt={m_so['tstt']:.3f}")
assert m_ue["feasible"] == 1.0 and m_so["feasible"] == 1.0
# MSA is non-certified research code -- the UE's own residual is disclosed,
# not asserted to zero (the certifier, not the solver, is the arbiter).
assert m_ue["tdue_gap"] < 0.15
assert np.isclose(m_ue["tstt"], 24.0, atol=1e-2)
# SO attains its LP-relaxation bound essentially exactly.
assert m_so["so_bound_gap"] < 1e-6
assert np.isclose(m_so["tstt"], 23.0, atol=1e-2)
# The headline: SO is strictly below UE, recomputed as a plain comparison of
# the two independently certified totals.
assert m_so["tstt"] < m_ue["tstt"]
print(f"SO/UE TSTT ratio : {m_so['tstt'] / m_ue['tstt']:.4f}  (SO strictly lower)")

scenario      : pm-wedge
content hash  : 2fa30ab1fb705308…
UE: tdue_gap=0.0909  tstt=24.000
SO: so_bound_gap=0.000e+00  tstt=23.000
SO/UE TSTT ratio : 0.9583  (SO strictly lower)

Visualize

Both figures come from tabench.viz — pm-td scenarios use the road Network (unlike transit-strategy’s separate TransitNetwork), so the house visualizer applies directly. The plotted flow is each link’s TIME- AVERAGED throughput, recomputed here by re-running PathLoader on the emitted departure plan (the same loading the certifier itself performs).

T = float(wedge.grid.edges[-1])
ue_out = PathLoader(wedge, ue_wedge.departures).run()
so_out = PathLoader(wedge, so_wedge.departures).run()
ue_avg = ue_out.n_out[:, -1] / T
so_avg = so_out.n_out[:, -1] / T

display(viz.plot_network_flows(wedge.network, ue_avg))
display(viz.plot_flow_scatter(("UE (selfish)", ue_avg), {"SO (system optimum)": so_avg}))
../../_images/60e90ecbe04cb5d5072af8ab1bbe703a70cedd4b2e85755bcbbae422c992ddf2.png ../../_images/9d30dd0700bcfa34f1a3b8837eb299b542f8e6b37e0ca84118e07f32e4dd6e14.png

Takeaways & pointers

  • Certified, not self-reported. Every gap and TSTT above came from TDTAEvaluator’s own re-loading and re-scoring — the MSA solver’s best- iterate claim was only diffed against it.

  • The MSA solvers are honestly non-certified. UE’s wedge residual (tdue_gap ~ 0.09) is disclosed, not hidden — a simulation-based MSA has no convergence guarantee (the paper says so); the certifier reports the best iterate, never a convergence claim.

  • Where next. the exit-function / LP twins merchant-nemhauser · lp-so-dta (the SO bound so_bound_gap scores against); the CTM/LTM loading underneath ctm; the lineage in the model compendium.