BurstGPT Queueing Simulation Study

A simulation study of LLM inference request scheduling using real traces from the BurstGPT dataset. Covers analytical benchmarks, DES simulation, factorial experiments, statistical analysis, variance reduction, and sensitivity analysis.

---

The Data & How We Model It

• Source: BurstGPT_1.csv, ~1.4 M real LLM requests with timestamps and response-token counts. Loaded via pandas.read_csv, sampled down to 50 000 rows for fitting.
• Inter-arrival time (IAT): np.diff(sorted_timestamps), then filtered to strictly positive values.
• Service time: response_tokens / TOKEN_THROUGHPUT, with TOKEN_THROUGHPUT = 14 tokens/s calibrated so baseline utilisation ρ ≈ 0.41 (a stable, moderately loaded server).
• Fitted parameters: λ = 1 / mean(IAT) ≈ 0.0574 req/s, μ = 1 / mean(svc) ≈ 0.1392 req/s.

Distribution fitting (src/data/loader.py). Five candidate distributions from scipy.stats are fit by MLE and ranked by Kolmogorov–Smirnov statistic:

_CANDIDATE_DISTS = {
    "expon":       stats.expon,
    "gamma":       stats.gamma,
    "weibull_min": stats.weibull_min,
    "lognorm":     stats.lognorm,
    "pareto":      stats.pareto,
}
# For each: params = dist.fit(data); ks, p = stats.kstest(data, dist(*params).cdf)
# AIC computed as 2k - 2·ΣlogPDF for reporting.

Both IAT and service time are best fit by lognormal (IAT KS ≈ 0.248, service KS ≈ 0.205). The tails are heavy — a handful of extreme requests dominate the system behavior.

---

Task 1 — How Bad Does Heavy-Tailed Service Make Things?

Goal: Get a back-of-envelope prediction of queueing delay before we simulate anything.

We use two closed-form queueing formulas (src/queueing/analytical.py):

• M/M/1 — assumes exponential service (CV² = 1). Simple, classical.
• M/G/1 (Pollaczek–Khinchine) — allows any service distribution. The key input is E[S²], the second moment of service time.

Implementation. Both are plain algebra on NumPy scalars, wrapped in a QueueMetrics dataclass:

# M/M/1
rho = lam / mu
L   = rho / (1 - rho)
Wq  = rho / (mu - lam)
Lq  = rho**2 / (1 - rho)

# M/G/1 Pollaczek–Khinchine
rho = lam * ES
Lq  = (lam**2 * ES2) / (2 * (1 - rho))      # ES2 = empirical E[S²]
Wq  = Lq / lam
W   = Wq + ES;  L = lam * W

We deliberately give both models the same E[S] = 1 / μ (empirical mean), so ρ is identical — only E[S²] differs. For M/M/1 we use E[S²] = 2 / μ² (exponential identity); for M/G/1 we use ES2 = np.mean(svc_data ** 2) — the raw empirical second moment. This isolates the effect of service-time variability as the sole difference.

Finding. The empirical CV² ≈ 1.85, much higher than exponential's CV² = 1. Plugging this into P-K gives a mean queue length ≈ 1.42× higher than M/M/1 would predict. In other words: the heavy tail alone inflates delay by 40%, before we even pick a scheduler.

📊 results/task1_distribution_fits.png, task1_analytical_benchmarks.png

---

Task 2 — Build a Simulator (and Trust It)

Goal: Build a discrete-event simulator flexible enough to test any scheduler — but first, prove it's correct.

Core design (src/simulation/des_engine.py). A minimal event-driven simulator using Python's heapq as the priority queue:

event_heap = []                              # (time, event_type, payload)
heapq.heappush(event_heap, (t, ARRIVAL, req_id))
t_now, etype, payload = heapq.heappop(event_heap)

Events are encoded as tuples with two types, ARRIVAL = 0 and DEPARTURE = 1. Per-request state (arrival / service-start / departure timestamps, dropped flag) is stored in a _Record dataclass, and aggregate metrics are computed in SimResult.analyse() after warm-up trimming via Little's Law:

self.mean_W  = np.mean(W_all)
self.mean_Wq = np.mean(Wq_all)
self.mean_L  = lam_eff * self.mean_W     # Little's Law
self.mean_Lq = lam_eff * self.mean_Wq
self.rho     = busy_time / total_time    # time-average utilisation

Schedulers (src/simulation/schedulers.py). All three implement the same enqueue / dequeue(n) / __len__ interface:

Policy	Data structure	Dispatch
FCFS	collections.deque	popleft() — O(1)
SJF	heapq min-heap keyed on service_time	heappop() — O(log n)
Batching	collections.deque	pop up to MAX_BATCH_SIZE = 8 in one dispatch

The Batching policy models GPU-style continuous batching: a batch departs at t_now + max(service_times in batch) (the slowest job holds up the batch — realistic for token-parallel LLM decoding).

RNG injection for CRN. run_simulation(...) accepts pre-generated iat_rvs and svc_rvs arrays rather than sampling internally. This is critical for Task 5: by passing in the same streams to different policies we get paired observations that dramatically tighten comparisons.

Validation. We feed the simulator exponential IAT and exponential service times generated via rng.exponential(1/λ, size=N) — i.e. the M/M/1 setting — and check that it reproduces the closed-form values for W, Wq, L, Lq, and ρ. All five metrics match theory within < 2 % relative error.

📊 results/task2_simulator_validation.png

---

Task 3 — A Controlled Factorial Experiment

Goal: Compare the three scheduling policies fairly across a range of loads and queue capacities.

Design: 3 × 3 × 2 = 18 configurations × 10 replications = 180 runs.

Factor	Levels
A — Scheduling policy	FCFS, SJF, Batching
B — Load multiplier on λ	0.7 ×, 1.0 ×, 1.3 ×
C — Queue capacity	∞, 100

Implementation (src/experiment/runner.py). The factorial runner builds a list of ExperimentConfig dataclasses, groups them by (load_mult, queue_cap), and iterates replications as the outer loop so that all three policies in a cell share the exact same random streams.

Inverse-CDF sampling (enables CRN). Rather than using the built-in rng.exponential / dist.rvs, we draw uniform variates and invert:

u_iat = rng.random(n)                     # U(0,1)^n  (PCG64 via default_rng)
u_svc = rng.random(n)
iat   = iat_dist.ppf(u_iat) / load_mult   # scale λ via the time axis
svc   = svc_dist.ppf(u_svc)
iat   = np.clip(iat, 1e-9, None)          # numerical floor
svc   = np.clip(svc, 1e-9, None)

where iat_dist and svc_dist are frozen scipy.stats lognormal distributions from Task 1's fit. Using .ppf() (the inverse CDF) guarantees that two different policies, given the same uniform stream, will encounter the same sequence of arrival/service events — the textbook requirement for Common Random Numbers to work as a variance reducer.

Seeding. Every np.random.default_rng(seed) is deterministic, with seeds derived as base_seed + rep * 1000 + int(load_mult * 100) + cap_tag so that different cells are independent but each (cell, rep) is reproducible.

Crucial design choice — Common Random Numbers (CRN). Within a given (load, capacity, replication) cell, all three policies see the same arrival and service streams. Any observed difference is caused purely by scheduling logic — not by one policy getting luckier random draws.

Finding. At baseline load (λ × 1.0, infinite queue):

Policy	Mean Wq (s)
SJF	≈ 147
Batching	≈ 252
FCFS	≈ 839

FCFS is ≈ 5.7× slower than SJF. Why? Under heavy-tailed service a single extreme request blocks everyone behind it — classical head-of-line blocking. SJF short-circuits this by serving cheap requests first.

📊 results/task3_factorial_results.png, task3_heatmap.png

---

Task 4 — Are the Differences Real? (Statistics)

Goal: Put proper confidence intervals around every policy mean and test whether observed differences are significant.

All statistical routines live in src/analysis/output_analysis.py.

4a. Confidence intervals

For each (policy, load, capacity) cell we compute a 95 % t-distribution CI over the 10 replications:

n  = len(data)
m  = np.mean(data)
se = np.std(data, ddof=1) / np.sqrt(n)        # sample std with Bessel correction
t  = scipy.stats.t.ppf(1 - alpha/2, df=n-1)    # two-sided critical value
half_width = t * se

For CRN-paired differences we additionally provide a paired-t CI on a − b (using the standard deviation of the paired differences) and a Welch CI for the independent-samples case (unequal-variance Satterthwaite df).

4b. One-way ANOVA + Bonferroni pairwise

Null hypothesis: H₀: µ_FCFS = µ_SJF = µ_Batching

F, p = scipy.stats.f_oneway(*groups)          # one-way ANOVA

Result: F = 12.36, p = 0.00015 — reject H₀ decisively.

Pairwise comparisons use a Bonferroni correction: with k = 3 pairs and overall α = 0.10, each pairwise test uses α / k ≈ 0.033. Implemented by calling scipy.stats.ttest_ind (Welch, equal_var=False) on every pair and flagging p < alpha/k:

Pair	Δ (s)	Significant?
FCFS vs SJF	+691	✅ yes
FCFS vs Batching	+587	✅ yes
Batching vs SJF	+104	❌ no (n.s.)

4c. Warm-up detection (MSER-m)

Every simulation has a transient period where the queue is empty and the metrics are biased low. We use MSER-m (Minimum Squared Error with batch means), which finds the truncation point d that minimises the estimated variance of the tail mean:

1. Chop the per-request Wq time series into non-overlapping batches of size batch_sz = 50.
2. Compute batch means to reduce noise.
3. Search over d ≤ n/2 for the d* minimising Σ(Yᵢ − Ȳ_{d:n})² / (n−d)².

We then discard the first N_WARMUP = 1000 arrivals from every replication when computing aggregate metrics (via SimResult.analyse(warmup_n=…)).

📊 results/task4_ci_comparison.png, task4_bonferroni.png, task4_warmup_detection.png

---

Task 5 — Squeeze More Information Out of the Same Runs (CRN)

Goal: Demonstrate a classical variance-reduction technique — Common Random Numbers — and quantify its benefit.

The idea. When comparing policy A vs policy B, feeding them the same random arrivals and service times means they face identical workloads. Any difference in their output is attributable entirely to policy — the noise from the random draws cancels out in the subtraction.

Mathematically, for paired observations (A_i, B_i):

Var[A − B] = Var[A] + Var[B] − 2·Cov[A, B]

With independent streams Cov[A, B] = 0. With CRN, Cov[A, B] > 0 (A and B tend to be large / small together), so the variance of the difference shrinks.

Implementation (run_crn_vs_independent in src/experiment/runner.py). Two nested loops over rep:

# CRN branch — one shared (iat, svc) stream per replication
rng_shared = np.random.default_rng(base_seed + rep*100)
iat_s, svc_s = _generate_rvs(rng_shared, n, iat_dist, svc_dist, load_mult)
for pol in policies:
    run_simulation(pol, iat_s.copy(), svc_s.copy(), ...)

# Independent branch — distinct seed per (rep, policy)
for i, pol in enumerate(policies):
    rng_ind = np.random.default_rng(base_seed + rep*100 + (i+1)*7919)
    iat_i, svc_i = _generate_rvs(rng_ind, n, ...)
    run_simulation(pol, iat_i, svc_i, ...)

The .copy() is necessary because run_simulation may mutate the service-time array internally. The 7919 prime offset keeps the independent seeds well separated.

We compare CIs on paired policy differences under:

• CRN: shared random streams within a replication.
• Independent: each policy gets its own random seed.

Finding. CRN shrinks the 95% CI width on paired differences by 20% – 71%, with the largest gain (–70.6%) for Batching vs SJF.

📊 results/task5_crn_comparison.png

---

Task 6 — Sensitivity Analysis: How Far Can We Trust Theory?

Goal: Show how the gap between theoretical and simulated waiting time depends on a single critical parameter — system load ρ.

Setup (src/experiment/sensitivity.py).

• Policy: FCFS; sweep ρ ∈ {0.50, 0.55, …, 0.95}; 40 replications per ρ.
• Arrivals: exponential IAT via rng.exponential(1/lam) — Poisson process, satisfying M/G/1 P-K theory's arrival assumption exactly.
• Service: empirical bootstrap via rng.choice(svc_data, replace=True). Bootstrap (not inverse-CDF of fitted lognormal) is used deliberately so the simulator's realized E[S] and E[S²] match the global moments used in the P-K formula — making sim vs theory a fair comparison.
• λ is set to rho * mu per ρ; each (ρ, rep) pair has a unique seed.

Sampling contrast vs Tasks 3–4. Tasks 3/4 use inverse-CDF sampling (dist.ppf(u)) of the fitted lognormal for both IAT and service times, enabling CRN across policies. Task 6 uses exponential IAT (to satisfy P-K assumptions) and empirical bootstrap for service (to anchor moments), with no CRN needed since only one policy is studied.

Finding. With 40 replications every ρ agrees with M/G/1 P-K to within ±10 %. A notable pattern emerges at high load: the relative error (sim_Wq − theory_Wq) / theory_Wq is systematically positive for ρ ≥ 0.85. This is not noise — it is a real, statistically significant upward bias driven by two compounding effects:

1. Bootstrap variance in E[S²]: each replication's realized second moment fluctuates around the global E[S²]; under heavy-tailed lognormal service the distribution of E[S²] is right-skewed, so the mean bootstrap realization slightly exceeds the global value.
2. Amplification near saturation: P-K Wq ∝ 1/(1−ρ), so even a small positive deviation in E[S²] is magnified greatly as ρ → 1.

With only 5 replications the rel-error curve was non-monotonic and uninformative (per-rep CV reaches 0.76 at ρ=0.9, giving SE ≈ 34 %). Increasing to 40 replications (SE ≈ 12 %) reveals the clean systematic pattern.

📊 results/task6_sensitivity.png

---

How to Run

python run_all.py                # all 6 tasks
python run_all.py --task 1       # Task 1 only
python run_all.py --task 3 4 5   # Tasks 3, 4, 5
python run_all.py --task 6       # Sensitivity analysis only

All figures and CSVs are written to results/.

Key Parameters (configs/default.py)

Parameter	Value	Meaning
SAMPLE_SIZE	50 000	rows used for distribution fitting
TOKEN_THROUGHPUT	14 tok/s	calibrates service time in seconds
N_ARRIVALS	5 000	arrivals per replication
N_WARMUP	1 000	arrivals discarded as transient
N_REPLICATIONS	10	replications per config
MAX_BATCH_SIZE	8	max parallel jobs in Batching policy

Implementation Stack

Concern	Choice	Why
Language	Python 3.13	ubiquity, mature scientific stack
Numerics	numpy	vectorised arrays, random generation
Distribution fitting	scipy.stats (.fit(), .kstest(), frozen dists)	MLE out-of-the-box + inverse CDF via .ppf()
Statistics	scipy.stats.t, .f_oneway, .ttest_ind	textbook tests without reinventing
RNG	numpy.random.default_rng (PCG64)	modern, stream-splittable, reproducible
Sampling scheme	Inverse-CDF: u = rng.random(n); x = dist.ppf(u)	enables Common Random Numbers across policies
Event loop	heapq min-heap of (time, event_type, payload) tuples	O(log n) per event, no external DES framework
FCFS queue	collections.deque	O(1) append / popleft
SJF queue	heapq keyed on service time	O(log n) insert / extract-min
Batching queue	collections.deque, dispatch up to MAX_BATCH_SIZE	mimics GPU continuous batching
Dataframes	pandas	tidy factorial results, groupby summaries, CSV round-trip
Plotting	matplotlib (Agg backend) + custom rcParams	publication-quality figures saved as PNG

File Structure

configs/default.py              ← tunable parameters
src/data/loader.py              ← CSV loader + distribution fitter
src/queueing/analytical.py      ← M/M/1 and M/G/1 closed forms
src/simulation/des_engine.py    ← heapq-based DES
src/simulation/schedulers.py    ← FCFS, SJF, Batching schedulers
src/experiment/runner.py        ← factorial + CRN experiment runner
src/experiment/sensitivity.py   ← Task 6 ρ-sweep
src/analysis/output_analysis.py ← CI, ANOVA, Bonferroni, warm-up
src/visualization/plots.py      ← all figures
run_all.py                      ← master script
results/                        ← generated figures and CSVs