Discrete Event Simulation (DES) in Python: A Practical Guide with SimPy

Heathrow Terminal 2 cost $3.2 billion to build. Before a single steel beam was raised, engineers ran discrete event simulation models of passengers walking, queueing, and scanning over a period of years. The simulations saved an estimated $200 million by identifying checkpoint layouts that would have failed during morning peaks. Amazon applies the same approach at a different scale: every new fulfilment centre is simulated with ten billion synthetic package routes before a single conveyor belt is installed. An emergency room in which the waiting time feels suspiciously predictable is often the product of similar work. Mayo Clinic, Cleveland Clinic, and most large hospital systems use DES to design triage flow so carefully that moving a single bed can reduce average patient wait times by thirty minutes.

Discrete event simulation is a quietly powerful technique that shapes billions of dollars of infrastructure, millions of patient-hours, and the back end of nearly every large logistics operation in the world. Most software engineers have nevertheless written no DES code. This guide aims to close that gap. It presents real, working simulations in Python using the SimPy library, covers the statistical machinery required to convert simulation noise into confident decisions, and connects DES to the adjacent worlds of optimisation and agent-based modelling so that the appropriate tool can be selected for each problem.

The Big Idea Behind Discrete Event Simulation

At its core, DES answers a question that analytical mathematics often cannot: how does a complex system with randomness, queues, and shared resources behave over time? Rather than writing a closed-form equation, an engineer builds a computer model of the system and lets simulated time advance, but only by jumping from one interesting moment, or “event,” to the next.

Consider a coffee shop. A customer arrives at minute 2.3. The barista starts service immediately. Service finishes at 4.7. Another customer arrives at 5.1, waits, begins service at 5.1, and finishes at 9.4. Between events, nothing changes; the simulation clock leaps forward to the next scheduled event. That leap is the basis of DES’s efficiency: a week of activity can be simulated in milliseconds because no cycles are spent on idle intervals between events.

DES Compared with Monte Carlo, System Dynamics, and Agent-Based Modelling

Newcomers often confuse DES with Monte Carlo simulation. The distinction is straightforward: Monte Carlo samples random outcomes from a distribution and aggregates statistics, but there is no evolving system state. Estimating the value of π by dropping random points into a square is Monte Carlo. It is elegant, but it lacks a time dimension. DES, by contrast, tracks how entities (customers, packets, patients) move through shared resources as simulated time advances.

System dynamics (SD) is a related approach. SD models continuous flows using differential equations: water levels in tanks may represent population or inventory, for example. SD is well suited to strategic, aggregate questions such as how advertising spend translates into market share over five years. SD cannot resolve individuals, however, and cannot answer questions such as how long patient #417 waited for the CT scanner. DES can.

Agent-based modelling (ABM) goes further than DES: each agent has autonomous behaviour, memory, and often geography. ABM is well suited to modelling crowd evacuation, epidemics, or economic actors that learn. DES agents, by contrast, are typically passive: they arrive, request a resource, are served, and leave. DES may be regarded as “ABM-lite with a global event queue.”

When DES Is Appropriate and When It Is Not

DES dominates wherever queues, shared resources, and randomness are present. Hospitals, call centres, manufacturing lines, supply chains, airports, data centre networks, and traffic corridors are all DES’s natural habitats. Questions of the form “how long will people or things wait?” or “what utilisation will this resource achieve?” or “what happens during peak demand?” are well suited to DES.

DES is not the appropriate tool when the underlying physics is continuous (fluid dynamics, electromagnetics, in which PDE solvers should be used), when the system is deterministic and small enough for a spreadsheet, or when a closed-form queueing result already exists. The classic M/M/1 queue, for example, has elegant analytical solutions: mean wait W = ρ/(μ(1−ρ)), where ρ = λ/μ. Simulating M/M/1 is useful primarily as a pedagogical exercise or as a sanity check on the simulation engine.

Core DES Concepts

Every DES model, whether written in SimPy or in a $30,000 commercial tool, shares the same vocabulary. Mastery of the following six concepts is sufficient to read any simulation paper in the literature.

Entities are the “things” that flow through the system: customers in a bank, packets in a router, patients in an ER, or pallets in a warehouse. Entities can have attributes (priority, size, type) that influence their routing.

Resources have limited capacity and hold entities while serving them. A single-teller bank has one resource of capacity 1; a hospital has dozens of specialised resources, including triage nurses, ER doctors, beds, and CT scanners. When an entity requests a busy resource, it joins a queue.

Events are moments at which the system state changes: an arrival, a service completion, a machine breakdown, or a shift change. Nothing happens between events; the clock skips through.

The future event list (FEL) is the priority queue, ordered by simulation time, that drives the entire engine. At each step the simulator pops the earliest event, executes its logic, and may schedule new events onto the FEL. When the FEL is empty or the clock passes the stop time, the simulation ends.

The simulation clock is simply a float. It has no relation to wall-clock time. A 24-hour call centre simulation may complete in 200 ms; a single second of a network-packet simulation may require an hour.

Statistics collection occurs continuously or at events: average wait time, maximum queue length, resource utilisation, throughput per hour, abandonment rate. These are the KPIs that stakeholders care about.

Randomness: The Heart of Stochastic Simulation

Real systems are noisy. Inter-arrival times between customers are not exactly six minutes; they follow a distribution. Service times vary. Machines break down at unpredictable moments. DES uses pseudo-random number generators (PRNGs) to sample from these distributions. Python’s random module or numpy.random is the typical source.

Two concepts confound nearly every beginner: the warm-up period and replications. When a simulation starts, it is in an unrealistic empty state, with no customers in the queue and all servers idle. Statistics gathered during this warm-up are biased toward low values. Practitioners discard the first X events, or X time units, before computing KPIs. Because every run uses different random numbers, a single simulation run is only one realisation of a random process. Replications (typically 20–100 independent runs with different seeds) and confidence intervals are required to support meaningful conclusions.

SimPy in Action: Four Complete Working Examples

SimPy is the Python DES library. It is free, open source, pure Python, and uses generator functions (yield-based) to express what would otherwise be callback spaghetti. Installation is via pip install simpy. The core idea is that every entity is a generator that yields timeouts or resource requests. SimPy’s environment orchestrates the event queue internally. Readers who value clean, readable code will appreciate SimPy. For more on writing code that the author’s future self will appreciate, see the guide on clean code principles for maintainable software.

Example 1: The M/M/1 Queue

The discussion begins with the textbook M/M/1 queue: one server, Poisson arrivals (mean inter-arrival 6 minutes), and exponential service (mean 5 minutes). The utilisation is ρ = 5/6 ≈ 0.83, which analytical queueing theory predicts should produce a mean wait of approximately 25 minutes.

import simpy
import random
import statistics

WAIT_TIMES = []

def customer(env, name, server, mean_service):
    arrival_time = env.now
    with server.request() as req:
        yield req                                   # wait for server
        wait = env.now - arrival_time
        WAIT_TIMES.append(wait)
        yield env.timeout(random.expovariate(1.0 / mean_service))

def arrival_process(env, server, mean_interarrival, mean_service):
    i = 0
    while True:
        yield env.timeout(random.expovariate(1.0 / mean_interarrival))
        i += 1
        env.process(customer(env, f'C{i}', server, mean_service))

def run_mm1(sim_time=10_000, seed=42):
    random.seed(seed)
    WAIT_TIMES.clear()
    env = simpy.Environment()
    server = simpy.Resource(env, capacity=1)
    env.process(arrival_process(env, server, 6, 5))
    env.run(until=sim_time)
    # discard warm-up (first 10%)
    warm = int(0.1 * len(WAIT_TIMES))
    stable = WAIT_TIMES[warm:]
    return statistics.mean(stable), len(stable)

mean_wait, n = run_mm1()
print(f"Avg wait: {mean_wait:.2f} min over {n} customers")
# Typical output: "Avg wait: 24.87 min over ~1500 customers"

The elegance is notable: twenty lines suffice for a full stochastic simulation with event-driven resource contention. The with server.request() as req: yield req pattern is idiomatic SimPy. It acquires the resource, automatically releases it when the with block exits, and handles queueing internally.

Example 2: Hospital Emergency Room

A real ER has multiple resource pools and priority-based routing. Patients undergo triage first and then compete for a doctor and a bed. Severity 1 (critical) patients preempt severity 3 (mild).

import simpy
import random
from collections import defaultdict

class ER:
    def __init__(self, env, n_triage=2, n_doctors=4, n_beds=10):
        self.env = env
        self.triage = simpy.Resource(env, n_triage)
        self.doctors = simpy.PriorityResource(env, n_doctors)
        self.beds = simpy.Resource(env, n_beds)
        self.wait_by_severity = defaultdict(list)
        self.treated = 0

def patient(env, pid, er):
    arrival = env.now
    severity = random.choices([1, 2, 3], weights=[0.1, 0.3, 0.6])[0]

    # Triage (every patient)
    with er.triage.request() as req:
        yield req
        yield env.timeout(random.triangular(2, 4, 8))

    # Bed + doctor — priority by severity (lower int = higher priority)
    with er.beds.request() as bed_req:
        yield bed_req
        with er.doctors.request(priority=severity) as doc_req:
            yield doc_req
            wait = env.now - arrival
            er.wait_by_severity[severity].append(wait)
            # severity-dependent treatment
            mean_treat = {1: 60, 2: 30, 3: 15}[severity]
            yield env.timeout(random.lognormvariate(
                mu=__import__('math').log(mean_treat), sigma=0.4))
            er.treated += 1

def arrivals(env, er, mean_iat=4.0):
    i = 0
    while True:
        yield env.timeout(random.expovariate(1.0 / mean_iat))
        i += 1
        env.process(patient(env, i, er))

random.seed(7)
env = simpy.Environment()
er = ER(env)
env.process(arrivals(env, er))
env.run(until=24 * 60)   # one day in minutes

for sev in sorted(er.wait_by_severity):
    waits = er.wait_by_severity[sev]
    print(f"Severity {sev}: n={len(waits):3d}  avg wait = "
          f"{sum(waits)/len(waits):.1f} min")
print(f"Total treated: {er.treated}")

Example 3: Manufacturing Line with Breakdowns

A three-workstation line is configured as cutting → assembly → packing, with a buffer between stations. Machines break down at random and are repaired. The question is a classic supply-chain problem, and the outputs feed directly into financial models. Many teams couple DES with time-series demand forecasting in order to close the planning loop.

import simpy, random

PROCESS_TIME = {'cut': 3, 'assm': 5, 'pack': 2}
MTBF = 120   # mean time between failures (min)
MTTR = 15    # mean time to repair

class Machine:
    def __init__(self, env, name, proc_time, buffer_in, buffer_out):
        self.env = env
        self.name = name
        self.proc_time = proc_time
        self.in_buf = buffer_in
        self.out_buf = buffer_out
        self.broken = False
        self.processed = 0
        env.process(self.run())
        env.process(self.breakdowns())

    def run(self):
        while True:
            part = yield self.in_buf.get()
            while self.broken:
                yield self.env.timeout(1)
            yield self.env.timeout(random.expovariate(1.0 / self.proc_time))
            yield self.out_buf.put(part)
            self.processed += 1

    def breakdowns(self):
        while True:
            yield self.env.timeout(random.expovariate(1.0 / MTBF))
            self.broken = True
            yield self.env.timeout(random.expovariate(1.0 / MTTR))
            self.broken = False

def raw_material_arrivals(env, buf):
    i = 0
    while True:
        yield env.timeout(random.expovariate(1.0 / 2.5))
        i += 1
        yield buf.put(f'Part-{i}')

random.seed(1)
env = simpy.Environment()
b0 = simpy.Store(env, capacity=20)   # raw
b1 = simpy.Store(env, capacity=10)   # between cut and assembly
b2 = simpy.Store(env, capacity=10)   # between assembly and pack
b3 = simpy.Store(env, capacity=1000) # finished goods

m1 = Machine(env, 'cut',  PROCESS_TIME['cut'],  b0, b1)
m2 = Machine(env, 'assm', PROCESS_TIME['assm'], b1, b2)
m3 = Machine(env, 'pack', PROCESS_TIME['pack'], b2, b3)

env.process(raw_material_arrivals(env, b0))
env.run(until=8 * 60)   # 8-hour shift

print(f"Cut: {m1.processed}   Assembly: {m2.processed}   Pack: {m3.processed}")
print(f"Finished goods: {len(b3.items)}")

Running the simulation reveals a classic lesson: the bottleneck (assembly, with a five-minute mean) dictates throughput. Adding a second cutter has no effect. The economic benefit lies in adding a second assembly station or in reducing assembly’s mean time by 20%. The insight is the kind that a spreadsheet cannot reliably surface.

Example 4: Call Centre with Abandonment

Call centres have time-varying arrival rates (morning peaks and lunch lulls), multi-skill routing, and, crucially, callers who hang up if they wait too long. The abandonment rate is a first-class KPI.

import simpy, random

# Hourly arrival rate (calls/min) for a 12-hour day
LAMBDA = [0.5, 0.8, 1.2, 1.8, 2.0, 1.8, 1.5, 1.3, 1.4, 1.2, 0.9, 0.6]
PATIENCE_MEAN = 3.0   # minutes before abandonment
SERVICE_MEAN  = 4.5

answered, abandoned, waits = 0, 0, []

def caller(env, agents):
    global answered, abandoned
    arrival = env.now
    patience = random.expovariate(1.0 / PATIENCE_MEAN)
    req = agents.request()
    result = yield req | env.timeout(patience)
    if req in result:
        wait = env.now - arrival
        waits.append(wait)
        answered += 1
        yield env.timeout(random.expovariate(1.0 / SERVICE_MEAN))
        agents.release(req)
    else:
        abandoned += 1
        req.cancel()

def arrivals(env, agents):
    while True:
        hour = int(env.now // 60) % 12
        rate = LAMBDA[hour]
        yield env.timeout(random.expovariate(rate))
        env.process(caller(env, agents))

random.seed(2026)
env = simpy.Environment()
agents = simpy.Resource(env, capacity=10)  # 10 agents all day
env.process(arrivals(env, agents))
env.run(until=12 * 60)

total = answered + abandoned
print(f"Answered: {answered}  Abandoned: {abandoned}  "
      f"Abandonment rate: {abandoned/total:.1%}")
print(f"Avg wait (answered): {sum(waits)/len(waits):.2f} min")

The elegant device is req | env.timeout(patience). SimPy’s | operator waits for either event, whichever fires first. A single line of code captures the entire logic of impatient callers.

Statistical Analysis of DES Output

This is the area in which most beginner simulations fail. The M/M/1 model is run once, “avg wait = 22.1 min” is observed, and the figure is reported. A second run with a different seed may yield 28.4. Which is correct? Neither. Both are samples from a random process, and a single sample is essentially useless.

Replications and Confidence Intervals

The standard remedy is to run N independent replications with different seeds, treat each replication’s mean as one observation, and compute the sample mean and 95% confidence interval.

import statistics, math

def replicate(n_reps=30, sim_time=10_000):
    means = []
    for seed in range(n_reps):
        m, _ = run_mm1(sim_time=sim_time, seed=seed)
        means.append(m)
    xbar = statistics.mean(means)
    s = statistics.stdev(means)
    half_width = 1.96 * s / math.sqrt(n_reps)   # 95% CI
    return xbar, (xbar - half_width, xbar + half_width)

mean, ci = replicate()
print(f"Mean wait = {mean:.2f}  95% CI: [{ci[0]:.2f}, {ci[1]:.2f}]")

If the CI width is too wide to distinguish scenarios, the number of replications or the simulation length should be increased. A useful rule of thumb is that halving the CI width requires quadrupling the number of replications.

Warm-Up Bias, Terminating and Steady-State Simulations

Two variants of simulation require different analysis. Terminating simulations have a natural end (a bank open from 9 to 5, or a single baseball game). For these, replication and averaging are sufficient. Steady-state simulations describe long-run behaviour (a 24/7 data centre). For steady-state simulations, the warm-up period should always be discarded. Welch’s method, in which the moving average is plotted and the point of stabilisation is identified visually, is the standard technique.

Comparing Scenarios

Consider the question “should two more agents be hired, or should the phone system be upgraded?” To compare Scenario A and Scenario B, common random numbers should be used: A and B are run with the same random seeds so that the only difference between them is the scenario itself. A paired t-test is then substantially more powerful than a comparison of two independent samples. The variance reduction technique alone can reduce the number of required replications by a factor of 5–10.

Real-World Applications

Many of the queues that one encounters in daily life were shaped by a DES model. The domains in which DES is industry standard are summarised below, together with the KPIs that practitioners focus on.

Several examples illustrate the impact. Mayo Clinic’s ER simulation reduced door-to-doctor time by 27% by reallocating triage nurses across shifts; no new hires were required, only better scheduling informed by DES. Toyota pioneered simulation-driven production line design in the 1980s, which partly explains why its lines continue to outperform competitors. TSMC simulates every new fab layout at the individual wafer level before construction; a single 3-nanometre fab costs $20 billion, and a layout error could cost billions in lost throughput. Amazon’s operations research team uses DES to determine how many robots to deploy per zone, balancing capital expenditure against peak-season throughput. FedEx’s Memphis superhub, the central facility of overnight shipping, was simulated down to the conveyor level before a single package moved through it.

In computer networking, simulators such as ns-3 and OMNeT++ are discrete event simulators at their core. Every paper that proposes a new TCP congestion control algorithm is backed by a DES model. For teams orchestrating large batches of such runs, Apache Airflow is well suited to managing the simulation pipeline.

DES with Optimisation: MIP, GA, and Sim-Opt Loops

DES answers the question “how does the system perform given these parameters?” The relevant business question, however, is usually “what parameters should be chosen?” That is optimisation. The two are complementary, and their combination yields the strongest economic results.

If the system is deterministic and linear, mixed-integer programming (MIP) can often find the global optimum directly. Real systems, however, have stochastic queues and nonlinear wait-time curves that MIP cannot capture. The standard pattern is therefore a simulation-optimisation loop: an outer optimiser proposes candidate parameter sets, and the DES model evaluates each by running replications and reporting KPIs.

For combinatorial search spaces, such as “which 10 of these 50 shift patterns should be used?”, genetic algorithms are a natural fit because they tolerate noisy fitness evaluations and handle discrete decision variables. Bayesian optimisation is well suited to continuous, expensive-to-evaluate parameters (such as the one-hour, three-replication DES evaluations common in industry). Commercial tools such as OptQuest bundle simulated annealing, tabu search, and scatter search into AnyLogic and Simio.

In recent years, reinforcement learning has been added to the mix: the DES model becomes an environment, and an RL agent learns policies (dispatch rules, dynamic pricing, inventory reorder points) that outperform hand-coded heuristics. DES combined with RL is currently among the most active research areas in operations research.

Tools Compared: SimPy, AnyLogic, Arena, and Others

SimPy is well suited to learners, researchers, and data teams that already work in Python. Production environments often use commercial tools for visualisation and GUI model builders. The landscape is summarised below.

For raw speed on very large simulations, Python is not the fastest option. For billions of packets or entities, a C++ framework (OMNeT++ or ns-3) or rewriting the hot path in a faster language should be considered. The Python vs Rust performance comparison discusses when that trade-off is justified. SimPy models nevertheless routinely process more than 100,000 entities per second on a laptop, which covers 95% of business cases.

Practical Tips and Common Pitfalls

Building one DES model is straightforward. Building one that stakeholders trust is more demanding. The following list identifies the practices that distinguish hobbyists from professionals.

Verification compared with validation. Verification asks “does the code do what was intended?”: unit tests, code review, and animation playback. Validation asks “does the model match reality?”: simulated KPIs are compared against historical data. A model can be verified (free of defects) but invalid (built on incorrect assumptions). Both procedures are required.

Use realistic distributions. Beginners default to exponential distributions everywhere because they are memoryless and mathematically convenient. Real service times are often lognormal or gamma, right-skewed with a long tail. Distributions should be fitted from data using scipy.stats or maximum likelihood. For storing and preprocessing historical data at scale, see the guide on databases for preprocessed time series.

Common defects. Forgetting to release a resource (early-return paths require attention). Confusing arrival rate λ with mean inter-arrival time 1/λ, a potential threefold error. Using random.random() without seeding, which produces irreproducible runs. Allowing warm-up bias to enter production reports.

Keep the model legible. DES models are read many more times than they are written, by auditors, new team members, and the original author at a later date. Entities and events should be named descriptively, the source of every distribution parameter should be commented (for example, “service time fitted from Q3 2025 log, n=28,441”), and everything should be version-controlled in accordance with solid Git practices.

Sensitivity analysis. A DES model has dozens of parameters, and stakeholders invariably ask “what if demand increases by 20%?” One parameter at a time should be varied, the response curve plotted, and the few parameters that materially affect KPIs identified. A related concern is anomaly detection on the input data feeding the model, since garbage in produces garbage out; the guide on time-series anomaly detection is a useful companion.

Frequently Asked Questions

DES vs Monte Carlo simulation, what’s the difference?

Monte Carlo samples random outcomes from distributions and aggregates statistics; there is no concept of time-evolving state. DES tracks entities moving through a system over simulated time, with events firing at specific moments and state changing discretely. If your problem has queues, resource contention, or time-dependent behavior, use DES. If it is pure probabilistic risk (e.g., estimating the VaR of a portfolio), Monte Carlo suffices.

How many replications do I need for valid DES results?

A practical rule is to start with 30 replications, compute the 95% confidence interval half-width, and decide whether it is narrow enough to distinguish the scenarios you care about. If not, quadruple the reps to halve the half-width. For high-stakes decisions (hospital layout, $100M facility), 100+ replications with common random numbers across scenarios is standard.

Can SimPy handle large industrial simulations?

Yes, for most business-scale problems—tens of thousands of concurrent entities and millions of events per hour of wall time are routine. For simulations requiring billions of entities or real-time constraints (5G network simulators, substantial wargames), commercial tools or C++ frameworks like ns-3 and OMNeT++ are better choices. Many teams prototype in SimPy and port the core engine to C++ only if profiling proves it necessary.

DES vs Agent-Based Modeling—when to use which?

DES is best when entities are passive, they flow through pre-defined paths, request resources, and depart. ABM is best when individuals make autonomous decisions, interact with neighbors, or have memory and learning. Hospital patient flow is DES. Pandemic spread with individual behavioral choice is ABM. Many modern tools (AnyLogic especially) let you combine both paradigms in one model.

How does DES integrate with optimization (MIP/GA)?

The standard pattern is a simulation-optimization loop: an outer optimizer—MIP for deterministic linear structure, genetic algorithms for combinatorial search, Bayesian optimization for expensive continuous parameters—proposes parameter sets, and the DES model evaluates each by running replications. The optimizer uses the KPI feedback to guide its next proposal. This hybrid approach captures stochastic queueing behavior that pure MIP cannot, while still finding near-optimal designs.

Closing Thoughts

Discrete event simulation is the often-overlooked workhorse behind emergency rooms that feel surprisingly well run, factories that meet throughput targets, and airports that frequently manage to clear security on time. It is the tool that engineers reach for when a system has queues, randomness, and shared resources, and when closed-form mathematics fails. SimPy provides Python with a DES library that is free, readable, and sufficiently capable for most real-world problems.

The recommended approach is to begin modestly. The M/M/1 example should be coded, verified against analytical results, and then extended one concept at a time: priority queues, multi-server resources, breakdowns, and time-varying arrivals. Within a week, models that answer real business questions can be built. Pairing DES with optimisation (MIP for structure and GA for combinatorial search) allows the transition from “how does this system behave?” to “what design should be built?”—and that transition is where DES proves its economic value.

This article is for informational and educational purposes only and should not be treated as financial or engineering advice. Always validate simulation models against real data before making capital-intensive decisions.

References and Further Reading

• SimPy Official Documentation—API reference, tutorials, and community examples.
• Banks, J., Carson, J. S., Nelson, B. L., Nicol, D. M. Discrete-Event System Simulation (5th ed.),the classic textbook for academic DES courses.
• Law, A. M. Simulation Modeling and Analysis (5th ed.)—the practitioner’s bible on input modeling, output analysis, and variance reduction.
• AnyLogic Learning Resources—free tutorials on DES, ABM, and SD modeling.
• INFORMS Simulation Society,the leading professional community for simulation research, with the annual Winter Simulation Conference.