traffic.js
133 lines · in natura-src-1.0.0.zip as natura/traffic.js · a module · its entry in the module reference · its test: traffic-test.js
// traffic.js -- Natura TRAFFIC (society / flow): the microscopic Intelligent Driver Model on a ring. Each driver
// accelerates toward a desired speed and brakes on the gap + closing-rate to the car ahead -- purely local, no
// coordination. Below a critical density the road runs free; above it a single brake-tap NUCLEATES a PHANTOM JAM,
// a stop-and-go wave that travels BACKWARD through forward-moving cars, though no one hit anything and no one
// stopped for a reason. The fundamental diagram (flow vs density) rises then falls -- the emergence. city.js has
// only a macroscopic volume-delay curve; this is the missing car-following layer. Deterministic. UMD.
(function (root, factory) {
if (typeof module !== 'undefined' && module.exports) module.exports = factory();
else root.Traffic = factory();
})(typeof self !== 'undefined' ? self : this, function () {
'use strict';
const clamp = (v, a, b) => v < a ? a : v > b ? b : v;
function mulberry32(a) { return function () { a |= 0; a = a + 0x6D2B79F5 | 0; let t = Math.imul(a ^ a >>> 15, 1 | a); t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t; return ((t ^ t >>> 14) >>> 0) / 4294967296; }; }
// config: { N(cars), roadLen, v0, aMax, bComf, T, s0, carLen, delta, jitter, dtInner, maxSub, secPerDay, seed, name }
function makeTraffic(config) {
config = config || {};
const N = config.N != null ? config.N : 30;
const L = config.roadLen != null ? config.roadLen : 400; // ring length (m)
const v0 = config.v0 != null ? config.v0 : 25; // desired speed (m/s)
const aMax = config.aMax != null ? config.aMax : 1.3; // max accel
const bComf = config.bComf != null ? config.bComf : 1.8; // comfortable decel
const T = config.T != null ? config.T : 1.3; // safe time headway (s)
const s0 = config.s0 != null ? config.s0 : 2.0; // min bumper gap (m)
const carLen = config.carLen != null ? config.carLen : 5; // car length (m)
const delta = config.delta != null ? config.delta : 4;
const jitter = config.jitter != null ? config.jitter : 0.6; // initial position jitter (seeds the instability)
const dtInner = config.dtInner != null ? config.dtInner : 0.3; // IDM needs small steps
const maxSub = config.maxSub != null ? config.maxSub : 400;
const secPerDay = config.secPerDay != null ? config.secPerDay : 86400;
// --- the ANALYTIC fundamental diagram: the equilibrium speed-density relation these car-following PARAMETERS
// imply, independent of the ring's momentary occupancy. At equilibrium (dv=0, leader matched) the IDM accel is
// zero when s0 + v*T = gap*sqrt(1 - (v/v0)^delta); the LHS rises and the RHS falls in v, so a bisection finds
// the unique equilibrium speed at each spacing. This is the car-following layer city.js asks for by name --
// a volume-delay curve GROUNDED in the microsim rather than an assumed BPR power law.
function equilibriumSpeed(gap) {
if (gap <= s0) return 0; // packed to the min gap -> jammed, v=0
let lo = 0, hi = v0;
for (let it = 0; it < 44; it++) {
const v = 0.5 * (lo + hi);
const rhs = gap * Math.sqrt(Math.max(0, 1 - Math.pow(v / v0, delta)));
const lhs = s0 + v * T;
if (lhs > rhs) hi = v; else lo = v; // f = lhs-rhs is increasing in v; seek f=0
}
return 0.5 * (lo + hi);
}
// build a VOLUME-DELAY table: cost multiplier (travel-time inflation v0/v_eq) as a function of utilization
// x = flow/capacity. Free branch (density up to the flow-peak) maps x in [0,1]; the congested branch beyond
// the peak is folded to x in (1,2] (x = 2 - flow/capacity) so the curve stays single-valued and monotonic
// past capacity -- the standard way to make a volume-delay function usable when demand exceeds capacity.
function buildFD() {
const rhoMax = 1 / (carLen + s0); // jam density (veh/m): bumper-to-bumper at the min gap
const M = 120, pts = []; let cap = 0, critRho = rhoMax * 0.5;
for (let k = 1; k <= M; k++) {
const rho = (k / M) * rhoMax * 0.999;
const v = equilibriumSpeed(1 / rho - carLen), q = rho * v; // veh/s
pts.push({ rho, v, q }); if (q > cap) { cap = q; critRho = rho; }
}
const xs = [], cs = [];
for (const p of pts) {
const cost = p.v > 1e-6 ? v0 / p.v : 40; // travel-time inflation, raw-capped so it stays finite
(p.rho <= critRho ? (xs.push(p.q / cap), cs.push(cost))
: (xs.push(2 - p.q / cap), cs.push(cost)));
}
const S = 33, vdCost = new Float64Array(S); // resample onto a uniform x-grid [0,2]
for (let s = 0; s < S; s++) {
const x = (s / (S - 1)) * 2;
if (x <= xs[0]) { vdCost[s] = cs[0]; continue; }
if (x >= xs[xs.length - 1]) { vdCost[s] = cs[cs.length - 1]; continue; }
let j = 1; while (j < xs.length && xs[j] < x) j++;
const t = (x - xs[j - 1]) / ((xs[j] - xs[j - 1]) || 1);
vdCost[s] = cs[j - 1] + t * (cs[j] - cs[j - 1]);
}
return { vdCost, vdXmax: 2, capacity: +(cap * 3600).toFixed(1), critDensity: +(critRho * 1000).toFixed(2),
jamDensity: +(rhoMax * 1000).toFixed(2), freeSpeed: v0 };
}
return {
name: config.name || 'traffic',
domain: 'society',
scale: 'local',
params: config,
init() {
const rng = mulberry32(config.seed || 1), x = new Float64Array(N), v = new Float64Array(N);
for (let i = 0; i < N; i++) { x[i] = (i / N) * L + (rng() - 0.5) * jitter; v[i] = v0 * 0.8; }
return { N, L, x, v, t: 0, _rng: rng, fd: buildFD() }; // fd = the analytic volume-delay curve (computed once; params fixed)
},
step(st, ctx, dt) {
const nsub = clamp(Math.round(dt / dtInner), 1, maxSub); // integrate dt at dtInner resolution
const { x, v } = st, a = st._a || (st._a = new Float64Array(N));
for (let s = 0; s < nsub; s++) {
// order cars around the ring to find each one's leader
const idx = st._idx || (st._idx = new Int32Array(N)); for (let i = 0; i < N; i++) idx[i] = i;
idx.sort((p, q) => x[p] - x[q]);
for (let o = 0; o < N; o++) {
const i = idx[o], lead = idx[(o + 1) % N];
let gap = x[lead] - x[i] - carLen; if (o === N - 1) gap += L; // wrap to the leader ahead
gap = Math.max(0.1, gap);
const dv = v[i] - v[lead]; // closing rate
const sStar = s0 + Math.max(0, v[i] * T + v[i] * dv / (2 * Math.sqrt(aMax * bComf)));
a[i] = aMax * (1 - Math.pow(v[i] / v0, delta) - (sStar / gap) * (sStar / gap));
}
for (let i = 0; i < N; i++) { v[i] = Math.max(0, v[i] + a[i] * dtInner); x[i] = (x[i] + v[i] * dtInner) % L; if (x[i] < 0) x[i] += L; }
}
st.t += dt / secPerDay;
return st;
},
// OUT: positions x, speeds v; meanSpeed, flow (veh/hr through a point), density (veh/km), jamCount, minGap
fields(st) {
const { x, v, N, L } = st; let sv = 0, jam = 0;
for (let i = 0; i < N; i++) { sv += v[i]; if (v[i] < v0Frac()) jam++; }
const mean = sv / N, density = N / (L / 1000), flow = density * mean * 3.6; // veh/km * m/s -> veh/hr
// min gap (collision-free check)
const idx = Array.from({ length: N }, (_, i) => i).sort((p, q) => x[p] - x[q]); let mg = Infinity;
for (let o = 0; o < N; o++) { const i = idx[o], lead = idx[(o + 1) % N]; let g = x[lead] - x[i] - carLen; if (o === N - 1) g += L; if (g < mg) mg = g; }
function v0Frac() { return (config.v0 != null ? config.v0 : 25) * 0.35; }
const fd = st.fd || (st.fd = buildFD());
return { N, positions: x, speeds: v, meanSpeed: +mean.toFixed(3), flow: +flow.toFixed(1),
density: +density.toFixed(2), jamCount: jam, minGap: +mg.toFixed(3),
// the analytic fundamental diagram this road implies -- a coupled system (city.js) reads vdCost as the
// volume-delay curve to ground its congestion in this microsim instead of an assumed BPR exponent.
vdCost: fd.vdCost, vdXmax: fd.vdXmax, capacity: fd.capacity, critDensity: fd.critDensity,
jamDensity: fd.jamDensity, freeSpeed: fd.freeSpeed };
},
events(st) { const f = this.fields(st); return f.jamCount > st.N * 0.2 ? [{ type: 'jam', jammed: f.jamCount }] : []; }
};
}
return { makeTraffic };
});