Written and maintained by CASRAI Editorial Board
Published
On this page
SimPy is an open-source, process-based discrete-event simulation (DES) framework for Python that models systems where operational states change at discrete points in time. Developed by Klaus Müller and Tony Vignaux and maintained by Team SimPy, the framework uses Python generator coroutines to model concurrent activities—such as customer queues, manufacturing cells, hospital clinics, and inventory buffers—without multithreading overhead. By organizing interactions around an event-driven clock and shared capacity-constrained resources, SimPy enables operations researchers and management analysts to evaluate bottlenecks, queue dynamics, and stochastic service systems.
Definition: Discrete-Event Simulation with SimPy
Discrete-event simulation models systems where state variables change instantaneously at distinct event timestamps, remaining static between events. Unlike continuous engines such as MATLAB and Simulink or domain solvers like EnergyPlus building energy simulation, SimPy does not step through fixed time slices (Δt) or solve differential equations. Instead, the simulation clock jumps directly to the scheduled time of the next earliest event on an internal priority queue.
Under SimPy’s process-oriented approach, active entities are defined as generator functions that run until encountering a yield expression. Yielding an event—such as a delay (env.timeout()), a resource request (resource.request()), or a signal—suspends the coroutine and returns control to simpy.Environment. When that event triggers, the environment reactivates the generator at the exact point of suspension, preserving all local variables.
Formula and Discrete-Event System Dynamics
The state of a service system at simulation time t ≥ 0 is defined by the state vector S(t):
S(t) = (N(t), Q(t), B(t))
Where N(t) is entities in system, Q(t) = max(0, N(t) – c) is queue length, B(t) = min(N(t), c) is busy servers, and c is server capacity (c = 1 for M/M/1).
SimPy advances the simulation clock t using the Next-Event Time Advance algorithm. For event calendar E = {e1, e2, …, em} with event timestamps t(ej), the next clock time is:
tk+1 = min { t(e) : e ∈ E }
For queueing systems with arrival rate λ and service rate μ, outputs are validated against steady-state queueing formulas:
- Traffic Intensity (Server Utilization, ρ): ρ = λ / (c × μ) (stability requires ρ < 1).
- Little’s Law: L = λ × W and Lq = λ × Wq, relating entity counts to waiting times.
- M/M/1 Mean Queue Wait Time (Wq): Wq = λ / [μ(μ – λ)] = [ρ / (1 – ρ)] × (1 / μ)
- M/M/1 Mean Queue Length (Lq): Lq = λ2 / [μ(μ – λ)] = ρ2 / (1 – ρ)
- M/M/1 Mean System Time (W): W = Wq + (1 / μ) = 1 / (μ – λ)
Worked Example with Real Numbers and Arithmetic Shown
Consider a customer service workstation modeled as an M/M/1 queue with c = 1 server. Customers arrive at rate λ = 2.0 per hour (mean interarrival time = 0.50 hours, or 30.0 minutes), and the server processes customers at rate μ = 3.0 per hour (mean service duration = 0.3333 hours, or 20.0 minutes).
1. Analytical Steady-State Baseline
- Server Utilization (ρ): ρ = 2.0 / (1 × 3.0) = 2/3 ≈ 0.6667 (66.67%).
- Mean Waiting Time in Queue (Wq): Wq = 2.0 / [3.0 × (3.0 – 2.0)] = 0.6667 hours (40.00 minutes).
- Mean Time in System (W): W = 1.0 / (3.0 – 2.0) = 1.0000 hour (60.00 minutes).
- Mean Queue Length (Lq): Lq = (0.6667)2 / (1 – 0.6667) = 1.3333 customers.
- Mean Number in System (L): L = 2.0 / (3.0 – 2.0) = 2.0000 customers.
2. Chronological Event Trace for 5 Sample Customers
The table below traces 5 customer arrivals through the single server, recording clock advances, waiting times, and departures (times in minutes):
| Customer | Interarrival (ΔA) | Arrival (A) | Service Start (S) | Wait Time (Wq = S – A) | Service Duration (D) | Departure (C = S + D) | System Time (T = C – A) |
|---|---|---|---|---|---|---|---|
| 1 | 10.0 | 10.0 | 10.0 | 0.0 | 15.0 | 25.0 | 15.0 |
| 2 | 8.0 | 18.0 | 25.0 | 7.0 | 20.0 | 45.0 | 27.0 |
| 3 | 12.0 | 30.0 | 45.0 | 15.0 | 10.0 | 55.0 | 25.0 |
| 4 | 30.0 | 60.0 | 60.0 | 0.0 | 18.0 | 78.0 | 18.0 |
| 5 | 10.0 | 70.0 | 78.0 | 8.0 | 12.0 | 90.0 | 20.0 |
3. Sample Summary Arithmetic
- Total Queue Wait Time: 0.0 + 7.0 + 15.0 + 0.0 + 8.0 = 30.0 minutes.
- Mean Sample Queue Wait (W̄q): 30.0 / 5 = 6.00 minutes (0.1000 hours).
- Total Time in System: 15.0 + 27.0 + 25.0 + 18.0 + 20.0 = 105.0 minutes.
- Mean Sample System Time (W̄): 105.0 / 5 = 21.00 minutes (0.3500 hours).
- Total Server Busy Time: 15.0 + 20.0 + 10.0 + 18.0 + 12.0 = 75.0 minutes.
- Total Server Idle Time: (10.0 – 0.0) + (60.0 – 55.0) = 15.0 minutes.
- Sample Server Utilization (ρ̂): 75.0 / 90.0 = 0.8333 (83.33%).
Interpretation and Reading Guidance
Interpreting simulation outputs requires distinguishing between small-sample transient fluctuations and long-run steady-state convergence. In the 5-customer sample, sample utilization (83.33%) exceeds steady-state utilization (66.67%) due to transient arrival clustering.
| Performance Metric | Simulated Value Range | Operational Interpretation & Management Action |
|---|---|---|
| Server Utilization (ρ) | < 0.50 | Underutilized capacity. Server is idle majority of operating time; consider resource consolidation. |
| Server Utilization (ρ) | 0.70 – 0.85 | Target operational zone. Balances high equipment utilization with stable, acceptable queue lengths. |
| Server Utilization (ρ) | > 0.90 | Congestion risk threshold. Small arrival bursts generate severe nonlinear increases in queue wait times. |
| Mean Queue Wait (Wq) | Within 95% CI of Wq,theory | Model verification established. Simulated queueing dynamics match theoretical assumptions. |
Assumptions and Common Mistakes
- Cold-Start Initialization Bias: Starting with empty queues biases early metrics downward; discard a warm-up period (t0) before logging data.
- Missing Generator Yields: Omitting
yieldcauses coroutines to run synchronously without waiting, freezing the event loop. - Resource Deadlocks: Failing to release resources blocks waiting processes; always use context managers (
with resource.request() as req: yield req). - Correlated Random Streams: Unseeded generators compromise replication independence; apply rigorous random sampling in Python with explicit seeds.
- Conflating Simulation Clock with Wall Time: SimPy’s
env.timeout()advances virtual time, whereastime.sleep()pauses the OS thread; never usetime.sleep(). - Underpowered Replication Count: Single runs lack statistical validity; run multiple replications and compute confidence intervals following sample size determination formulas.
How to Report Simulation Results (APA Example)
When reporting discrete-event simulation experiments in academic management journals adhering to APA 7th edition standards, document the software version, pseudorandom generator, warm-up period, replication count, and parameter confidence intervals:
“A discrete-event simulation of the single-server service facility was conducted using SimPy (Version 4.1.1; Müller et al., 2018) in Python 3.12. The system was parameterized with a Poisson arrival process (λ = 2.00 arrivals/hr) and exponential service times (μ = 3.00 services/hr), yielding theoretical traffic intensity ρ = .67. To eliminate initialization transient bias, an initial warm-up period of 200 simulation hours was discarded prior to data collection. Across n = 30 independent replications of 1,000 operational hours each (Mersenne Twister PRNG), simulated mean queue waiting time was Wq = 0.669 hr (95% CI [0.648, 0.690]), converging with the theoretical steady-state value of 0.667 hr, t(29) = 0.22, p = .828. Average server utilization was ρ̂ = .668 (95% CI [.659, .677]), verifying model validity.”
How to Calculate and Simulate in Excel, R or Python
Below are verified implementations of an M/M/1 discrete-event queue across Python, R, and Excel.
Python (SimPy 4 Implementation)
import random
import simpy
def customer(env, name, server, wait_times, mu):
arrive_time = env.now
with server.request() as req:
yield req
wait_times.append(env.now - arrive_time)
yield env.timeout(random.expovariate(mu))
def arrival_generator(env, server, wait_times, lambd, mu):
i = 0
while True:
yield env.timeout(random.expovariate(lambd))
i += 1
env.process(customer(env, f"Cust_{i}", server, wait_times, mu))
def run_simulation(sim_time=1000.0, lambd=2.0, mu=3.0, seed=42):
random.seed(seed)
env = simpy.Environment()
server = simpy.Resource(env, capacity=1)
wait_times = []
env.process(arrival_generator(env, server, wait_times, lambd, mu))
env.run(until=sim_time)
return sum(wait_times) / len(wait_times) if wait_times else 0.0
avg_wait = run_simulation()
print(f"Mean Wait Time: {avg_wait:.4f} hours")
R (Discrete-Event Simulation with simmer)
In R, discrete-event queueing models can be implemented using the simmer framework, complementing broader statistical modeling in R and econometric workflows such as the GARCH time-series model:
library(simmer)
set.seed(42)
env <- simmer("MM1_System")
customer_traj <- trajectory() %>%
seize("server", 1) %>%
timeout(function() rexp(1, rate = 3.0)) %>%
release("server", 1)
env %>%
add_resource("server", 1) %>%
add_generator("arrival", customer_traj, function() rexp(1, rate = 2.0)) %>%
run(until = 1000)
arrivals <- get_mon_arrivals(env)
avg_wait <- mean(arrivals$end_time - arrivals$start_time - arrivals$activity_time)
cat(sprintf("Simulated Mean Wait: %.4f hoursn", avg_wait))
Excel (Discrete-Event Recurrence Model)
In Excel, a single-server queue is evaluated with row-wise formulas:
- Row 2 (Customer 1): Col B (Interarrival):
=-LN(RAND())/2.0; Col C (Arrival):=B2; Col D (Duration):=-LN(RAND())/3.0; Col E (Service Start):=C2; Col F (Wait):=E2-C2; Col G (Departure):=E2+D2. - Row 3 (Customer i ≥ 2, fill to Row 1001): Col C (Arrival):
=C2+B3; Col E (Service Start):=MAX(C3, G2); Col F (Wait):=E3-C3; Col G (Departure):=E3+D3. - Summary: Mean Queue Wait:
=AVERAGE(F2:F1001); Server Utilization:=SUM(D2:D1001)/MAX(G2:G1001).
Frequently Asked Questions
What is SimPy used for in quantitative research?
SimPy models complex operations research systems—such as supply chains, manufacturing cells, and hospital clinics—where stochastic arrivals and shared resource constraints create queueing bottlenecks that cannot be solved analytically.
How does SimPy differ from continuous simulation tools?
Continuous tools advance in fixed steps (Δt) to solve differential equations. SimPy advances via next-event scheduling, jumping instantaneously to the next event timestamp without calculating inactive periods.
Why does SimPy use Python generators instead of threading?
Generator coroutines (yield) eliminate OS thread context-switching overhead, memory usage, and race conditions, allowing thousands of concurrent entities to execute safely within a single thread.
What types of shared resources does SimPy support?
SimPy provides Resource for discrete servers with fixed capacity, Container for continuous homogeneous bulk materials, and Store for discrete inventory item pools.
How do you validate a SimPy discrete-event simulation model?
Validation requires verifying outputs against analytical queueing laws (such as Little’s Law and M/M/1 formulas), discarding initial warm-up transient data, and computing confidence intervals across independent replications.
Sources
- SimPy Development Team. (2023). SimPy: Discrete-Event Simulation in Python (Version 4.1.1). Read the Docs. Retrieved November 9, 2026, from https://simpy.readthedocs.io/en/latest/
- PyPI. (2023). simpy 4.1.1: Event Discrete, Process Based Simulation for Python. Python Package Index. Retrieved November 9, 2026, from https://pypi.org/project/simpy/
- Matloff, N. (2008). Introduction to Discrete-Event Simulation and the SimPy Language. Department of Computer Science, University of California, Davis. Retrieved November 9, 2026, from https://web.cs.ucdavis.edu/~matloff/matloff/public_html/156/PLN/DESimIntro.pdf
- Banks, J., Carson, J. S., Nelson, B. L., & Nicol, D. M. (2010). Discrete-Event System Simulation (5th ed.). Pearson.
- Law, A. M. (2015). Simulation Modeling and Analysis (5th ed.). McGraw-Hill.







