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

Heathrow Terminal 2 cost $3.2 billion to build — and before a single steel beam went up, engineers ran discrete event simulation models of passengers walking, queueing, and scanning for years. The simulations saved an estimated $200 million by flagging checkpoint layouts that would have melted down during morning peaks. Amazon does the same thing at a different scale: every new fulfillment center is simulated with ten billion synthetic package routes before a single conveyor belt is installed. And if you have ever sat in an emergency room where the wait felt suspiciously predictable — it probably was. Mayo Clinic, Cleveland Clinic, and most large hospital systems use DES to design triage flow so carefully that moving a single bed can shave thirty minutes off average patient waits.

Discrete event simulation is one of those quietly powerful techniques that shapes billions of dollars of infrastructure, millions of patient-hours, and the back-end of nearly every large logistics operation in the world — yet most software engineers have never written a line of DES code. That ends today. In this guide we will build real, working simulations in Python using the SimPy library, cover the statistical machinery that turns simulation noise into confident decisions, and connect DES to the adjacent worlds of optimization and agent-based modeling so you know when to reach for which tool.

The Big Idea Behind Discrete Event Simulation

At its heart, DES answers a question that analytical math often cannot: how does a complex system with randomness, queues, and shared resources actually behave over time? Instead of writing a closed-form equation, you build a computer model of the system and let simulated time march forward — but only by jumping from one interesting moment (an “event”) to the next.

Imagine 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, gets served starting at 5.1, finishes at 9.4. Between events, nothing changes — so the simulation clock simply leaps forward to the next scheduled event. That leap is the secret to DES’s efficiency: you can simulate a week of activity in milliseconds because you never waste cycles on “idle” time between events.

DES vs Monte Carlo, System Dynamics, and Agent-Based Modeling

Newcomers often confuse DES with Monte Carlo simulation. The easiest way to separate them: Monte Carlo samples random outcomes from a distribution and aggregates statistics, but there is no evolving system state. If you estimate the value of π by dropping random points into a square, that is Monte Carlo — beautiful, but time-less. DES, by contrast, tracks how entities (customers, packets, patients) move through shared resources as simulated time advances.

System dynamics (SD) is another cousin. SD models continuous flows using differential equations — think of water levels in tanks representing “population” or “inventory.” SD is great for strategic, aggregate questions like “how does advertising spend translate into market share over five years?” But SD cannot see individuals, so it cannot answer “how long did patient #417 wait for the CT scanner?” DES can.

Agent-based modeling (ABM) goes further than DES: each agent has autonomous behavior, memory, and often geography. ABM is ideal for modeling crowd evacuation, epidemics, or economic actors who learn. DES agents, by contrast, are usually passive — they arrive, request a resource, get served, and leave. You can think of DES as “ABM-lite with a global event queue.”

When DES Shines and When It Doesn’t

DES dominates wherever you have queues, shared resources, and randomness. Hospitals, call centers, manufacturing lines, supply chains, airports, data center networks, and traffic corridors are all DES’s natural habitat. If your question involves “how long will people or things wait?” or “what utilization will this resource hit?” or “what happens during peak demand?” — DES is your tool.

DES is not the right tool when the underlying physics is continuous (fluid dynamics, electromagnetics — use PDE solvers), when the system is deterministic and small enough for a spreadsheet, or when a closed-form queueing result already exists. Classic M/M/1 queues, for example, have elegant analytical solutions: mean wait W = ρ/(μ(1−ρ)) where ρ = λ/μ. Simulating M/M/1 is mostly useful as a pedagogical exercise or a sanity check on your simulation engine.

Core DES Concepts You Must Know

Every DES model, whether written in SimPy or a $30,000 commercial tool, shares the same vocabulary. Master these six concepts and you can read any simulation paper in the literature.

Entities are the “things” flowing through the system. Customers in a bank, packets in a router, patients in an ER, 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 specialized resources — triage nurses, ER doctors, beds, CT scanners. When an entity requests a busy resource, it joins a queue.

Events are moments when the system state changes: an arrival, a service completion, a machine breakdown, a shift change. Everything between events is nothing — the clock skips straight through.

The future event list (FEL) is the priority queue (ordered by simulation time) that drives the whole engine. At each step the simulator pops the earliest event, executes its logic, which 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 just a float — it has nothing to do with wall-clock time. A twenty-four-hour call-center simulation may take 200 ms to run; a single second of a network-packet simulation may take an hour.

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

Randomness: The Heart of Stochastic Simulation

Real systems are noisy. Inter-arrival times between customers are not exactly every 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 are the usual sources.

Two concepts that trip up every beginner: the warm-up period and replications. When a simulation starts, it’s in an unrealistic empty state — no customers in queue, all servers idle. Statistics gathered during this warm-up are biased toward low values. Professionals discard the first X events (or X time units) before computing KPIs. And because every run uses different random numbers, a single simulation run is just one realization of a random process. You need replications (typically 20–100 independent runs with different seeds) and confidence intervals to say anything meaningful.

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. Install with pip install simpy. The core idea: every entity is a generator that yields timeouts or resource requests. SimPy’s environment orchestrates the event queue under the hood. If you love clean, readable code, you will love SimPy — and for more on writing code your future self will thank you for, see our guide on clean code principles for maintainable software.

Example 1: The M/M/1 Queue

Let us start with the textbook M/M/1 queue: one server, Poisson arrivals (mean inter-arrival 6 minutes), exponential service (mean 5 minutes). Utilization ρ = 5/6 ≈ 0.83, which analytical queueing theory says should give a mean wait of about 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"

Notice the elegance: twenty lines and you have 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 for you.

Example 2: Hospital Emergency Room

A real ER has multiple resource pools and priority-based routing. Patients go through triage first, 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: cutting → assembly → packing, with a buffer between stations. Machines break down randomly and are repaired. This is a classic supply-chain question, and the outputs feed directly into financial models — many teams couple DES with time-series demand forecasting 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 this reveals a classic lesson: the bottleneck (assembly, 5-minute mean) dictates throughput. Adding a second cutter does nothing. Adding a second assembly station or reducing assembly’s mean time by 20% is where the money is. This is the kind of insight you never get from a spreadsheet.

Example 4: Call Center with Abandonment

Call centers have time-varying arrival rates (morning peaks, lunch lulls), multi-skill routing, and — crucially — callers who hang up if they wait too long. 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 beautiful trick here is req | env.timeout(patience) — SimPy’s | operator waits for either event, whichever fires first. That one line of code captures the entire logic of impatient callers.

Statistical Analysis of DES Output

This is where most beginner simulations fall apart. You run the M/M/1 model once, see “avg wait = 22.1 min,” and report it. But run it again with a different seed and you might see 28.4. Which is right? Neither. They are samples from a random process, and a single sample is nearly useless.

Replications and Confidence Intervals

The standard remedy: 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 your CI width is too wide to distinguish scenarios, increase the number of replications or the simulation length. A handy rule: to halve the CI width, quadruple the number of replications.

Warm-Up Bias and Terminating vs Steady-State

Two flavors of simulation require different analysis. Terminating simulations have a natural end (a bank open 9 to 5, a single baseball game) — just replicate and average. Steady-state simulations are meant to describe long-run behavior (a 24/7 data center). For steady-state, always discard the warm-up. Welch’s method (plot the moving average and eyeball when it stabilizes) is the standard technique.

Comparing Scenarios

“Should I hire two more agents or upgrade the phone system?” To compare Scenario A vs B, use common random numbers: run A and B with the same random seeds so the only difference between them is the scenario itself. Then a paired t-test is far more powerful than comparing two independent samples. This variance reduction trick alone can cut required replications by 5–10×.

Real-World Applications That Shape Your Life

Every time you have waited somewhere, a DES model probably influenced the layout. Here are the domains where DES is industry standard, each with the KPIs practitioners obsess over.

A few stories worth telling. Mayo Clinic’s ER simulation reduced door-to-doctor time by 27% by reallocating triage nurses across shifts — zero new hires, just better scheduling informed by DES. Toyota pioneered simulation-driven production line design in the 1980s, which is part of why their lines still out-throughput competitors. TSMC simulates every new fab layout at the individual wafer level before construction; a single 3-nanometer fab costs $20 billion, and a layout error could cost billions in lost throughput. Amazon’s operations research team uses DES to decide how many robots to deploy per zone, balancing capex against peak-season throughput. FedEx’s Memphis superhub — the beating heart of overnight shipping — was simulated down to the conveyor level before a single package moved through it.

In computer networking, simulators like ns-3 and OMNeT++ are actually discrete event simulators under the hood. Every time you read a paper proposing a new TCP congestion control algorithm, there is a DES model backing the numbers. If you are orchestrating large batches of such runs, Apache Airflow can manage the simulation pipeline beautifully.

DES Meets Optimization: MIP, GA, and Sim-Opt Loops

DES answers “how does the system perform given these parameters?” But the real question is usually “what parameters should I choose?” That is optimization. The two are complementary, and combining them is where the serious money gets made.

If your system is deterministic and linear, you can often use mixed-integer programming (MIP) to find the global optimum directly. But real systems have stochastic queues and nonlinear wait-time curves that MIP cannot capture. In that case, the standard pattern is a simulation-optimization loop: an outer optimizer proposes candidate parameter sets, and the DES model evaluates each one by running replications and reporting KPIs.

For combinatorial search spaces — “which 10 of these 50 shift patterns should I use?” — genetic algorithms are a natural fit because they tolerate noisy fitness evaluations and handle discrete decision variables. Bayesian optimization is great for continuous, expensive-to-evaluate parameters (like the one-hour-and-three-rep DES evaluations common in industry). Commercial tools like OptQuest bundle simulated annealing, tabu search, and scatter search into AnyLogic and Simio.

In the last few years, reinforcement learning has entered 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 + RL is currently one of the hottest research areas in operations research.

Tools Compared: SimPy, AnyLogic, Arena, and More

SimPy is perfect for learners, researchers, and data teams already living in Python. But production shops often use commercial tools for the visualization and GUI model-builders. Here is the landscape.

For raw speed on very large simulations, Python is not the fastest option. If you are simulating billions of packets or entities, consider a C++ framework (OMNeT++, ns-3) or even rewriting the hot path in a faster language — see our Python vs Rust performance comparison for when that trade-off is worth it. That said, SimPy models routinely run 100,000+ entities per second on a laptop, which covers 95% of business cases.

Practical Tips and Common Pitfalls

Building one DES model is easy. Building one that stakeholders trust is hard. Here is a curated list of things that separate hobbyists from professionals.

Verification vs Validation. Verification asks “does the code do what I intended?” — unit tests, code review, animation playback. Validation asks “does the model match reality?” — compare simulated KPIs against historical data. A model can be verified (bug-free) but invalid (wrong assumptions). Always do both.

Use real distributions. Beginners default to exponential everything because it is memoryless and mathematically convenient. Real service times are often lognormal or gamma — right-skewed with a long tail. Fit your distributions from data using scipy.stats or maximum likelihood. For storing and preprocessing that historical data at scale, see our guide on databases for preprocessed time series.

Classic bugs. Forgetting to release a resource (watch for early-return paths). Mixing arrival rate λ with mean inter-arrival time 1/λ — a 3× error waiting to happen. Using random.random() without seeding — irreproducible runs. Letting warm-up bias sneak into production reports.

Keep the model legible. DES models are read many more times than they are written — by auditors, new team members, and future you. Name entities and events descriptively, comment the source of every distribution parameter (“service time fit from Q3 2025 log, n=28,441”), and version-control everything with solid Git practices.

Sensitivity analysis. A DES model has dozens of parameters, and stakeholders always ask “what if demand goes up 20%?” Vary one parameter at a time, plot the response curve, and identify the few parameters that move KPIs meaningfully. A related idea is anomaly detection on the input data feeding your model — garbage in, garbage out — and our piece on time-series anomaly detection is a good companion there.

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, massive 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.

Conclusion

Discrete event simulation is the unsung workhorse behind emergency rooms that feel oddly well-run, factories that hit their throughput targets, and airports that almost manage to get you through security on time. It is the tool engineers reach for when a system has queues, randomness, and shared resources — and when closed-form math gives up. With SimPy, Python has a DES library that is free, readable, and powerful enough for most real-world problems.

Start small. Code up the M/M/1 example, verify against analytical results, and then expand one concept at a time: priority queues, multi-server resources, breakdowns, time-varying arrivals. Within a week you can be building models that answer real business questions. Pair DES with optimization (MIP for structure, GA for combinatorial search) and you can move from “how does this system behave?” to “what design should we build?” — and that jump is where DES earns its keep.

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.