Discrete-event simulation engine

Build a resource-constrained model with DESModel, for a scarce shared resource that makes patients queue.

Keywords

discrete-event simulation, DES, queueing model, resource-constrained model, patient waiting time, cost-effectiveness modeling, health economics, Python

This tutorial shows how to build a resource-constrained model with DESModel. Use this model type when a shared, scarce resource makes patients queue, so that one patient’s outcome depends on the others in the system. Neither a cohort model nor a microsimulation of independent individuals can represent that dependence. The model then runs through the same cost-effectiveness and value-of-information analysis as any other. The microsimulation tutorial covers two ideas this one builds on: building the model once and running it on draws, and accruing discounted costs and effects over time. This page walks through examples/des.py step by step. A discrete-event simulation without resource constraints is a continuous-time state-transition model, which the continuous clock of MicrosimModel runs directly. The discrete-event simulation replication reproduces a published example of that kind.

Modeling the queueing process

The model is a specialist clinic operating over a four-year horizon. Patients arrive as a Poisson process at rate ARRIVAL_RATE, so successive inter-arrival gaps are independent and exponentially distributed with mean $1/\lambda$,

\[\text{gap} \sim \text{Exponential}(\lambda),\]

which the patients function draws and accumulates into arrival times. Each patient waits for a specialist, a SimPy Resource. While waiting, untreated disease costs c_wait_year per year at utility u_wait. Once seen, a one-time treatment cost is incurred and the patient spends the rest of the horizon treated at the higher u_treated. DESModel builds on SimPy. The environment, the process, and the resource remain your own SimPy code, and the toolkit argument adds discounted accrual, seeding, and an event log.

import numpy as np
import pandas as pd
import simpy

HORIZON = 4.0
ARRIVAL_RATE = 15.0  # patients per year

def patients(rng, n):
    return pd.DataFrame({"arrival": np.cumsum(rng.exponential(1.0 / ARRIVAL_RATE, n))})

def resources(env, params, intervention):
    return {"specialist": simpy.Resource(env, capacity=int(params["n_servers"]))}

def clinic(env, patient, params, intervention, toolkit):
    arrival = float(patient["arrival"])
    if arrival >= toolkit.horizon:
        return
    yield env.timeout(arrival)
    toolkit.state("waiting")
    with toolkit.request("specialist") as slot:
        result = yield slot | env.timeout(toolkit.horizon - env.now)  # race the horizon
        served = slot in result
        toolkit.accrue_over(arrival, env.now, params["c_wait_year"], params["u_wait"])
        if served:
            toolkit.accrue_cost(params["c_treat"] + params["c_capacity"])
            toolkit.accrue_over(env.now, toolkit.horizon, params["c_followup_year"],
                                params["u_treated"])
            yield env.timeout(toolkit.rng.gamma(2.0, params["service_time"] / 2.0))

toolkit.accrue_over(start, end, cost_rate, utility_rate) integrates and discounts a flow over an absolute interval. Once the specialist is available, one call bills the time the patient spent waiting. The process reads the run’s time horizon from toolkit.horizon rather than repeating the value the engine already holds. A patient still waiting at the horizon is never seen, and only the waiting segment is billed.

Configuring once, evaluating on draws

The model goes to the engine when you build it, and the draw matrix goes to run_psa. The two interventions differ only in capacity. Each is an Intervention carrying the numeric values its arm sets: n_servers, plus a per-patient overhead c_capacity for the expanded arm. Capacity is a parameter the model already reads, not a flag standing in for which arm is which. The two interventions use the same patients and service draws through common random numbers, so the incremental result reflects the capacity change, not sampling noise.

from heormodel.models import DESModel, Intervention
from heormodel.params import Beta, Gamma, ParameterSet
from heormodel.run import SeedManager, run_psa

seeds = SeedManager(20260704)
parameters = ParameterSet({
    "service_time": Gamma.from_mean_se(0.08, 0.015),
    "u_wait": Beta.from_mean_se(0.60, 0.05),
    "u_treated": Beta.from_mean_se(0.85, 0.03),
    "c_wait_year": Gamma.from_mean_se(9_000.0, 1_500.0),
    "c_treat": Gamma.from_mean_se(5_000.0, 800.0),
    "c_followup_year": Gamma.from_mean_se(1_500.0, 300.0),
    "c_capacity": Gamma.from_mean_se(3_500.0, 600.0),
})
draws = parameters.sample(128, seed=seeds.generator())

engine = DESModel(
    process=clinic, population=patients, n_individuals=30, resources=resources,
    interventions=[Intervention("Standard capacity", {"n_servers": 1, "c_capacity": 0.0}),
                Intervention("Expanded capacity", {"n_servers": 2})],
    horizon=HORIZON,
)
outcomes = run_psa(engine, draws, seed=seeds.entropy, sequential=True).outcomes
outcomes.summary().round(2)

	cost	qaly
intervention
Standard capacity	11314.54	2.27
Expanded capacity	12700.84	2.35

Analyzing cost-effectiveness, value of information, and the queue

From here nothing is engine-specific. The same icer_table, evpi, and evppi_ranking calls used for a cohort or microsimulation model apply to this queueing model’s Outcomes. Expanded capacity costs more and increases quality-adjusted life-years, because patients spend less time waiting at the low utility.

from heormodel.cea import icer_table
icer_table(outcomes).round(2)

	cost	cost_lo	cost_hi	effect	effect_lo	effect_hi	inc_cost	inc_cost_lo	inc_cost_hi	inc_effect	inc_effect_lo	inc_effect_hi	icer	icer_lo	icer_hi	status
intervention
Standard capacity	11314.54	7805.64	17120.89	2.27	1.93	2.58	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	ND
Expanded capacity	12700.84	10073.07	15489.75	2.35	1.98	2.62	1386.3	-2627.67	3628.68	0.07	0.01	0.2	18690.46	-17954.49	410553.66	ND

from heormodel.voi import evpi, evppi_ranking

WTP = 30_000.0
print(f"EVPI at WTP {WTP:,.0f}: {evpi(outcomes, WTP):,.1f}")
evppi_ranking(outcomes, draws, WTP).round(1)

EVPI at WTP 30,000: 783.6

service_time       563.9
c_wait_year         59.4
c_capacity          43.1
u_treated           24.6
c_treat             24.3
u_wait              18.3
c_followup_year      7.0
Name: evppi, dtype: float64

Service time leads the value-of-information ranking. It sets how fast the queue clears, and so determines whether the extra capacity is worth its cost. Passing collect="events" to run_psa returns the event log alongside the outcomes, and queue_waits converts it into per-request waits.

from heormodel.models import queue_waits

trace = run_psa(engine, draws.iloc[[0]], seed=seeds.entropy, collect="events").events
(queue_waits(trace).groupby("intervention", sort=False)["wait"].mean() * 365).round(1)

intervention
Standard capacity    27.3
Expanded capacity     1.1
Name: wait, dtype: float64

With a single specialist, patients wait weeks. A second specialist cuts the typical wait to about a day. The other model engine tutorials cover the cohort state-transition, microsimulation, and compartmental engines.

Next: the full pipeline tutorial walks through the whole workflow, from ParameterSet sampling through run_psa to the analysis functions.

Keywords: discrete-event simulation, DES, queueing model, resource-constrained model, patient waiting time, cost-effectiveness modeling, health economics, Python

--- title: Discrete-event simulation engine description: "Build a resource-constrained model with `DESModel`, for a scarce shared resource that makes patients queue." keywords: [discrete-event simulation, DES, queueing model, resource-constrained model, patient waiting time, cost-effectiveness modeling, health economics, Python] ---  <a href="https://colab.research.google.com/github/pedroliman/heormodel/blob/main/docs/_notebooks/des.ipynb"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"></a>  This tutorial shows how to build a resource-constrained model with `DESModel`. Use this model type when a shared, scarce resource makes patients queue, so that one patient's outcome depends on the others in the system. Neither a cohort model nor a microsimulation of independent individuals can represent that dependence. The model then runs through the same cost-effectiveness and value-of-information analysis as any other. The [microsimulation tutorial](microsim.qmd) covers two ideas this one builds on: building the model once and running it on draws, and accruing discounted costs and effects over time. This page walks through [`examples/des.py`](https://github.com/pedroliman/heormodel/blob/main/examples/des.py) step by step. A discrete-event simulation without resource constraints is a continuous-time state-transition model, which the continuous clock of `MicrosimModel` runs directly. The [discrete-event simulation replication](mdm-des.qmd) reproduces a published example of that kind. ## Modeling the queueing process The model is a specialist clinic operating over a four-year horizon. Patients arrive as a Poisson process at rate `ARRIVAL_RATE`, so successive inter-arrival gaps are independent and exponentially distributed with mean $1/\lambda$, $$\text{gap} \sim \text{Exponential}(\lambda),$$ which the `patients` function draws and accumulates into arrival times. Each patient waits for a specialist, a SimPy `Resource`. While waiting, untreated disease costs `c_wait_year` per year at utility `u_wait`. Once seen, a one-time treatment cost is incurred and the patient spends the rest of the horizon treated at the higher `u_treated`. `DESModel` builds on SimPy. The environment, the process, and the resource remain your own SimPy code, and the `toolkit` argument adds discounted accrual, seeding, and an event log. ```{python} import numpy as np import pandas as pd import simpy HORIZON = 4.0 ARRIVAL_RATE = 15.0 # patients per year def patients(rng, n): return pd.DataFrame({"arrival": np.cumsum(rng.exponential(1.0 / ARRIVAL_RATE, n))}) def resources(env, params, intervention): return {"specialist": simpy.Resource(env, capacity=int(params["n_servers"]))} def clinic(env, patient, params, intervention, toolkit): arrival = float(patient["arrival"]) if arrival >= toolkit.horizon: return yield env.timeout(arrival) toolkit.state("waiting") with toolkit.request("specialist") as slot: result = yield slot | env.timeout(toolkit.horizon - env.now) # race the horizon served = slot in result toolkit.accrue_over(arrival, env.now, params["c_wait_year"], params["u_wait"]) if served: toolkit.accrue_cost(params["c_treat"] + params["c_capacity"]) toolkit.accrue_over(env.now, toolkit.horizon, params["c_followup_year"], params["u_treated"]) yield env.timeout(toolkit.rng.gamma(2.0, params["service_time"] / 2.0)) ``` `toolkit.accrue_over(start, end, cost_rate, utility_rate)` integrates and discounts a flow over an absolute interval. Once the specialist is available, one call bills the time the patient spent waiting. The process reads the run's time horizon from `toolkit.horizon` rather than repeating the value the engine already holds. A patient still waiting at the horizon is never seen, and only the waiting segment is billed. ## Configuring once, evaluating on draws The model goes to the engine when you build it, and the draw matrix goes to `run_psa`. The two interventions differ only in capacity. Each is an `Intervention` carrying the numeric values its arm sets: `n_servers`, plus a per-patient overhead `c_capacity` for the expanded arm. Capacity is a parameter the model already reads, not a flag standing in for which arm is which. The two interventions use the same patients and service draws through common random numbers, so the incremental result reflects the capacity change, not sampling noise. ```{python} from heormodel.models import DESModel, Intervention from heormodel.params import Beta, Gamma, ParameterSet from heormodel.run import SeedManager, run_psa seeds = SeedManager(20260704) parameters = ParameterSet({ "service_time": Gamma.from_mean_se(0.08, 0.015), "u_wait": Beta.from_mean_se(0.60, 0.05), "u_treated": Beta.from_mean_se(0.85, 0.03), "c_wait_year": Gamma.from_mean_se(9_000.0, 1_500.0), "c_treat": Gamma.from_mean_se(5_000.0, 800.0), "c_followup_year": Gamma.from_mean_se(1_500.0, 300.0), "c_capacity": Gamma.from_mean_se(3_500.0, 600.0), }) draws = parameters.sample(128, seed=seeds.generator()) engine = DESModel( process=clinic, population=patients, n_individuals=30, resources=resources, interventions=[Intervention("Standard capacity", {"n_servers": 1, "c_capacity": 0.0}), Intervention("Expanded capacity", {"n_servers": 2})], horizon=HORIZON, ) outcomes = run_psa(engine, draws, seed=seeds.entropy, sequential=True).outcomes outcomes.summary().round(2) ``` ## Analyzing cost-effectiveness, value of information, and the queue From here nothing is engine-specific. The same `icer_table`, `evpi`, and `evppi_ranking` calls used for a cohort or microsimulation model apply to this queueing model's `Outcomes`. Expanded capacity costs more and increases quality-adjusted life-years, because patients spend less time waiting at the low utility. ```{python} from heormodel.cea import icer_table icer_table(outcomes).round(2) ``` ```{python} from heormodel.voi import evpi, evppi_ranking WTP = 30_000.0 print(f"EVPI at WTP {WTP:,.0f}: {evpi(outcomes, WTP):,.1f}") evppi_ranking(outcomes, draws, WTP).round(1) ``` Service time leads the value-of-information ranking. It sets how fast the queue clears, and so determines whether the extra capacity is worth its cost. Passing `collect="events"` to `run_psa` returns the event log alongside the outcomes, and `queue_waits` converts it into per-request waits. ```{python} from heormodel.models import queue_waits trace = run_psa(engine, draws.iloc[[0]], seed=seeds.entropy, collect="events").events (queue_waits(trace).groupby("intervention", sort=False)["wait"].mean() * 365).round(1) ``` With a single specialist, patients wait weeks. A second specialist cuts the typical wait to about a day. The other [model engine tutorials](index.qmd) cover the cohort state-transition, microsimulation, and compartmental engines. Next: the [full pipeline tutorial](full-pipeline.qmd) walks through the whole workflow, from `ParameterSet` sampling through `run_psa` to the analysis functions.