Queues, servers, conveyors, vehicles and resources simulated event by event: utilization, waiting times, throughput and their confidence ranges.
How the precision is kept
Event times and statistics are 64-bit. The results are estimates from random draws: how close they are is set by the run length and the number of replications, which the confidence half-width reports. Uniform draws carry 32 random bits, percentiles come from a stored sample of at most 50,000 values, and the Student t factor is a 3-decimal table.
Open item. Scale the time offsets on shift, schedule and conveyor edges to the model's own time unit.
Calculations under this policy
15 calculations in 2 source files. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
engine.jsShared shell, every app64-bit11 calculationsReviewed
The average of a value that changes at events: queue length, units busy, work in process.
average = (Σ vᵢ × Δtᵢ) ÷ (t − t₀)
The area is a running sum of value × interval in 64-bit, added in event order. Each interval is a difference of two clock times, so its rounding grows with the size of the clock value.
Waiting time, service time and time in system: mean, standard deviation, minimum, maximum, P50 and P95.
m ← m + (x − m) ÷ n; s ← s + (x − m_old)(x − m_new); sd = √(s ÷ (n − 1)); P = sorted[⌊p × k⌋]
Welford's update avoids subtracting two large sums. The 50th and 95th percentiles are read from a stored copy of the first 50,000 observations, with no interpolation between neighbors; later observations still count in the mean and deviation but not in the percentiles.
Each event time is the current clock plus a delay, so times are built by repeated addition and carry one rounding per event in the chain. Ties are ordered by priority, then by a sequence number.
The 1e-9 added before rounding down keeps a length that is a whole number of pitches from losing a slot to rounding. It is absolute, not scaled to the ratio.
x = t − ⌊t ÷ period⌋ × period; on shift when a ≤ x < b
The position in the period is found by subtracting a multiple of the period from the clock, which loses digits as the clock grows. The next edge is the first one more than 1e-12 after now, and the state is sampled 1e-9 after the edge; both offsets are absolute in the model's time unit.
engine.js
lines 151 to 157
const sh = s.shifts, inWin = (t) => { const x = t - Math.floor(t / sh.period) * sh.period; return sh.windows.some(([a, b]) => x >= a && x < b); };
const edges = [...new Set(sh.windows.flatMap(([a, b]) => [a, b]))].filter((x) => x < sh.period).sort((a, b) => a - b);
const next = (t) => { const base = Math.floor(t / sh.period) * sh.period; for (const x of edges) if (base + x > t + 1e-12) return base + x; return base + sh.period + edges[0]; };
const set = (on) => { if (on === this.onShift) return; this.onShift = on; this.refresh(); };
this.onShift = inWin(0); this.mark();
const tick = () => { set(inWin(this.sim.t + 1e-9)); this.sim.at(next(this.sim.t) - this.sim.t, tick, -1); };
this.sim.at(next(0), tick, -1); this.refresh();
Server utilization
64-bit
The share of time a server's units were busy, blocked, down and idle.
utilization = average busy units ÷ units; idle = max(0, 1 − busy − blocked − down)
Each share is a time average divided by the unit count. Idle is the remainder, clamped at zero because the three shares are rounded separately.
The table has 3 decimals. Above 30 degrees of freedom the factor is a step function, for example 2.0 at 31 where the exact value is about 2.040, so the half-width is good to about 2 digits there.
The random numbers every arrival, service time and route choice is drawn from.
u = (k + 0.5) ÷ 2³², k a 32-bit integer; z = √(−2 ln u) × cos(2π v)
A uniform draw takes one of 2³² evenly spaced values and is never 0 or 1. The normal draw is the Box–Muller formula on two such uniforms, so its largest possible size is √(−2 ln(0.5 ÷ 2³²)), about 6.8.
Inverse-transform formulas in 64-bit on a 32-bit uniform. A normal draw below zero is redrawn up to 50 times and then replaced by the mean. The Erlang draw multiplies k uniforms before one logarithm, which can underflow to zero for k of several hundred.
rng.js
lines 42 to 57
export function sample(d, rng) {
if (typeof d === 'number') return d;
const name = Object.keys(d)[0], v = d[name];
switch (name) {
case 'const': return v;
case 'exp': return -v * Math.log(rng.next());
case 'uniform': return v[0] + (v[1] - v[0]) * rng.next();
case 'tri': { const [a, c, b] = v, u = rng.next(), f = (c - a) / (b - a); return u < f ? a + Math.sqrt(u * (b - a) * (c - a)) : b - Math.sqrt((1 - u) * (b - a) * (b - c)); }
case 'normal': { for (let i = 0; i < 50; i++) { const x = v[0] + v[1] * rng.normal(); if (x >= 0) return x; } return v[0]; }
case 'lognormal': { const s2 = Math.log(1 + (v[1] / v[0]) ** 2), mu = Math.log(v[0]) - s2 / 2; return Math.exp(mu + Math.sqrt(s2) * rng.normal()); }
case 'erlang': { let p = 1; for (let i = 0; i < v[0]; i++) p *= rng.next(); return -(v[1] / v[0]) * Math.log(p); }
case 'weibull': return v[1] * Math.pow(-Math.log(rng.next()), 1 / v[0]);
case 'empirical': { const tot = v.reduce((s, p) => s + p[1], 0); let u = rng.next() * tot; for (const p of v) { u -= p[1]; if (u <= 0) return p[0]; } return v[v.length - 1][0]; }
}
return NaN;
}
Mean of a distribution
64-bit
The mean used for load estimates and for routing to the least loaded target.
tri: (a + m + b) ÷ 3; weibull: scale × Γ(1 + 1 ÷ shape)
Closed forms. For the normal distribution the given mean is returned although draws are cut at zero, so the mean of the draws is higher than this value when the deviation is large against the mean.
rng.js
lines 58 to 73
/** Mean of a distribution (used for load estimates and for checks). */
export function mean(d) {
if (typeof d === 'number') return d;
const name = Object.keys(d)[0], v = d[name];
switch (name) {
case 'const': case 'exp': return v;
case 'uniform': return (v[0] + v[1]) / 2;
case 'tri': return (v[0] + v[1] + v[2]) / 3;
case 'normal': return v[0];
case 'lognormal': return v[0];
case 'erlang': return v[1];
case 'weibull': { return v[1] * gammaFn(1 + 1 / v[0]); }
case 'empirical': { const tot = v.reduce((s, p) => s + p[1], 0); return v.reduce((s, p) => s + p[0] * p[1], 0) / tot; }
}
return NaN;
}
Gamma function
64-bit
Γ(x), used for the mean of a Weibull distribution.
Lanczos approximation, g = 7, 9 coefficients
A fixed 9-term series, so its error is that of the approximation, not of 64-bit arithmetic; the same routine in the risk engine carries a comment giving about 1e-14 relative. Arguments below 0.5 use the reflection formula.
rng.js
lines 74 to 80
/** Lanczos approximation of Γ(x) for x > 0. */
export function gammaFn(x) {
const g = 7, c = [0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7];
if (x < 0.5) return Math.PI / (Math.sin(Math.PI * x) * gammaFn(1 - x));
x -= 1; let a = c[0]; const t = x + g + 0.5; for (let i = 1; i < g + 2; i++) a += c[i] / (x + i);
return Math.sqrt(2 * Math.PI) * Math.pow(t, x + 0.5) * Math.exp(-t) * a;
}
In detail
Method
Event scheduling on a binary-heap calendar (the file cites Law, Simulation Modeling and Analysis). Time averages are integrals of a step function; means and variances use Welford's running update; replications use the seeds seed, seed + 1 and so on.
Number types
64-bit numbers for time and statistics. Random integers are 32-bit (xoshiro128**), turned into uniforms with 2³² possible values.
Tolerances
Shift and schedule edges use absolute offsets of 1e-12 and 1e-9 time units; conveyor timing 1e-12; the conveyor slot count adds 1e-9 before rounding down. Time has no fixed unit, so these are not scaled to the model.
Checked by
The engines test runs a single-server queue example through replicate and compares utilization, wait, time in system and queue length with the closed forms, within 3 to 20 percent.
Limits
Percentiles use the first 50,000 observations only. The t factor is 2.0, 1.98 or 1.96 above 30 degrees of freedom.
Approach: 64-bit, stated limit
The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.