vi-due — Friesz et al.’s (1993) variational-inequality dynamic user equilibrium

What. vi-due recasts dynamic user equilibrium as a variational inequality (VI) over path-departure flows: N travelers choose BOTH a route r (a free-flow time f_r followed by its own Vickrey point-queue bottleneck of capacity s_r) AND a departure time, trading queue delay against schedule delay around a common t*. This benchmark instantiates the VI with the generalized Vickrey point-queue loading, closed-form and FIFO by construction.

Why it is in the benchmark. It is the simultaneous route-AND-departure-time DUE (SRDC-DUE) — the multi-route generalization of vickrey, and it reduces EXACTLY to vickrey when there is one route with f = 0, which this notebook verifies numerically rather than asserting. See the model compendium (Friesz et al. 1993) and docs/design/adr-022-vi-due.md (P1).

Scope. This notebook runs the closed-form SRDC-DUE on the built-in two- route worked instance, certifies it, and demonstrates the exact vickrey reduction and the analytic route-inclusion threshold.

Canon. [friesz1993variational], 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 DUEEvaluator from the emitted per-route cumulative departure curves alone — reconstructing the deterministic point queue on EACH route, scoring every used traveler’s cost by level-inversion, AND checking that no traveler could improve by switching ROUTE (not just departure time) via a marginal-insertion sweep. The solver’s own C/split/window provenance is never trusted (README, Certified, not self-reported).

# Setup. `vi-due` 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 (
    DUEEvaluator,
    DUEScenario,
    BottleneckEvaluator,
    due_closed_form,
    friesz_two_route_scenario,
    ue_closed_form,
    vickrey_worked_scenario,
)

The scenario

The built-in worked instance (adr-022): N=6000, two parallel routes — route 0 (f=0.2, s=3000) and route 1 (f=0.7, s=1500) — same schedule-delay slopes as vickrey (alpha=1, beta=0.5, gamma=2, t*=9). DUEScenario is frozen and content-hashed (P2), domain-separated from the single-route BottleneckScenario hash space.

scenario = friesz_two_route_scenario()
print(f"scenario      : {scenario.name}")
print(f"content hash  : {scenario.content_hash()[:16]}…")
print(f"N={scenario.n_travelers}, routes f={list(scenario.route_free_flow)}, "
      f"s={list(scenario.route_capacity)}")
print(f"alpha={scenario.alpha}, beta={scenario.beta}, gamma={scenario.gamma}, "
      f"t*={scenario.t_star}")
c_star, used, n_r = scenario.equilibrium_structure()
print(f"common cost C (method)   : {c_star}")
print(f"used routes / split      : {used}  {n_r}")
print("task          : simultaneous route-and-departure-time user equilibrium")

scenario      : friesz-two-route
content hash  : 74c911fe30d6f52b…
N=6000.0, routes f=[np.float64(0.2), np.float64(0.7)], s=[np.float64(3000.0), np.float64(1500.0)]
alpha=1.0, beta=0.5, gamma=2.0, t*=9.0
common cost C (method)   : 0.9
used routes / split      : [ True  True]  [5250.  750.]
task          : simultaneous route-and-departure-time user equilibrium

Solve

No Budget/RngBundle/Trace — the closed form is a pure function of the scenario, emitting a DUEProfile (per-route cumulative departure curves, the P1-certifiable artifact).

profile = due_closed_form(scenario)
print(f"provenance : {profile.provenance}")
print(f"routes x grid : {profile.cumulative.shape}")

provenance : {'equilibrium_cost': 0.9, 'total_cost': 5400.0, 'n_route_0': 5250.0, 'n_route_1': 750.0000000000001}
routes x grid : (2, 2006)

Certify (P1) — no traveler can improve by switching route OR time

DUEEvaluator reconstructs each route’s point queue from its emitted curve, scores every used traveler’s cost by level-inversion, and additionally sweeps the MARGINAL-INSERTION cost of every route at a dense set of candidate times (an infinitesimal extra traveler does not move the curves) — due_gap = (max cost over used travelers − min marginal cost anywhere) / C, zero iff no traveler can improve by shifting time OR switching route.

evaluator = DUEEvaluator(scenario)
metrics = evaluator.certify(profile)
print(f"feasible          : {metrics['feasible']:.0f}")
print(f"due_gap           : {metrics['due_gap']:.3e}")
print(f"total_cost        : {metrics['total_cost']:.2f}")
print(f"expected_cost (C) : {metrics['expected_cost']:.4f}")
assert metrics["feasible"] == 1.0
assert metrics["due_gap"] < 1e-6
assert np.isclose(metrics["expected_cost"], 0.9, atol=1e-3)
assert np.isclose(metrics["total_cost"], 5400.0, atol=1e-1)

# Recompute the split directly from the emitted curves (route volumes), not
# quoting scenario.equilibrium_structure() a second time.
n_emitted = profile.cumulative[:, -1]
print(f"emitted route volumes : {n_emitted}")
np.testing.assert_allclose(n_emitted, [5250.0, 750.0], atol=1e-2)

feasible          : 1
due_gap           : 7.303e-14
total_cost        : 5400.00
expected_cost (C) : 0.9000
emitted route volumes : [5250.  750.]

Distinctive result 1: the route-inclusion threshold

A route joins the used set only while its free-flow cost alpha*f_r stays below the common cost level. Recomputed here directly from the scenario’s parameters (not quoted): with delta = beta*gamma/(beta+gamma), route 1 joins iff N > alpha*s_0*(f_1 - f_0)/delta.

delta = scenario.beta * scenario.gamma / (scenario.beta + scenario.gamma)
threshold = (
    scenario.alpha * scenario.route_capacity[0]
    * (scenario.route_free_flow[1] - scenario.route_free_flow[0]) / delta
)
print(f"route-inclusion threshold N* : {threshold}")
assert np.isclose(threshold, 3750.0, atol=1e-6)
assert scenario.n_travelers > threshold  # both routes used on the worked instance, verified above

below = DUEScenario(
    name="vi-due-below-threshold", n_travelers=threshold - 1.0,
    alpha=scenario.alpha, beta=scenario.beta, gamma=scenario.gamma, t_star=scenario.t_star,
    route_free_flow=scenario.route_free_flow, route_capacity=scenario.route_capacity,
)
_, used_below, n_r_below = below.equilibrium_structure()
print(f"below threshold (N={below.n_travelers}): used={used_below}, split={n_r_below}")
assert not used_below[1]  # route 1 excluded just below the threshold

route-inclusion threshold N* : 3749.999999999999
below threshold (N=3748.999999999999): used=[ True False], split=[3749.    0.]

Distinctive result 2: the exact reduction to vickrey

vi-due reduces exactly to vickrey when there is one route with f = 0 — verified here by running due_closed_form on a single-route DUEScenario at vickrey’s own parameters and diffing the certified totals against the vickrey anchor, computed independently by ue_closed_form + BottleneckEvaluator (not the same code path).

vk_scenario = vickrey_worked_scenario()
single_route = DUEScenario(
    name="vi-due-single-route", n_travelers=vk_scenario.n_travelers,
    alpha=vk_scenario.alpha, beta=vk_scenario.beta, gamma=vk_scenario.gamma,
    t_star=vk_scenario.t_star,
    route_free_flow=np.array([0.0]), route_capacity=np.array([vk_scenario.capacity]),
)
single_profile = due_closed_form(single_route)
single_metrics = DUEEvaluator(single_route).certify(single_profile)

vk_metrics = BottleneckEvaluator(vk_scenario).certify(ue_closed_form(vk_scenario))

print(f"vi-due (1 route, f=0) : C={single_metrics['expected_cost']:.6f}  "
      f"total={single_metrics['total_cost']:.4f}  max_queue={single_metrics['max_queue']:.4f}")
print(f"vickrey UE             : C={vk_scenario.equilibrium_cost:.6f}  "
      f"total={vk_metrics['total_cost']:.4f}  max_queue={vk_metrics['max_queue']:.4f}")
assert np.isclose(single_metrics["expected_cost"], vk_scenario.equilibrium_cost, atol=1e-6)
assert np.isclose(single_metrics["total_cost"], vk_metrics["total_cost"], atol=1e-2)
assert np.isclose(single_metrics["max_queue"], vk_metrics["max_queue"], atol=1e-2)
print("vi-due (single route, f=0) reduces exactly to the vickrey anchor")

vi-due (1 route, f=0) : C=0.800000  total=4800.0000  max_queue=2400.0000
vickrey UE             : C=0.800000  total=4800.0000  max_queue=2400.0000
vi-due (single route, f=0) reduces exactly to the vickrey anchor

Visualize

As in vickrey, tabench.viz is a road-network visualizer and this is a scalar departure-time model with no network — the house departure-diagram plots the two certified per-route cumulative departure curves directly.

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(5.5, 4.0))
for r in range(scenario.n_routes):
    ax.plot(profile.times, profile.cumulative[r],
            label=f"route {r} (f={scenario.route_free_flow[r]}, N_r={n_emitted[r]:.0f})")
ax.axvline(scenario.t_star, color="0.6", linestyle=":", label="t*")
ax.set_xlabel("time")
ax.set_ylabel("cumulative departures")
ax.set_title("vi-due: per-route SRDC-DUE departure schedules")
ax.legend(fontsize=8)
fig.tight_layout()
display(fig)
plt.close(fig)

Takeaways & pointers

• Certified, not self-reported. The route split and total cost came from DUEEvaluator, which checks BOTH the departure-time margin and the route margin from the emitted curves alone.

• The reduction is exact, not a limiting case with a fudge tolerance — a single f=0 route certifies to the vickrey anchor within closed-form float precision.

• Where next. the single-route special case vickrey; the lineage in the model compendium.