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, the model type to use when a shared, scarce resource makes patients queue, so that one patient’s outcome depends on the others in the system, a coupling neither a cohort nor an independent microsimulation model can represent. It runs that model through the same cost-effectiveness and value-of-information analysis as any other model; the microsimulation tutorial covers 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 stay the user’s own 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, so billing the queueing time a patient just endured is one call once the specialist is granted. The process reads the run’s time horizon back as toolkit.horizon, so it need not repeat the value the engine already holds. A patient still waiting at the horizon is never seen: the horizon arrives before the specialist, and only the waiting segment is billed.

Configuring once, evaluating on draws

The model goes to the engine when you build it; the draw matrix goes to run_psa. The two interventions differ only in capacity, so each is an Intervention carrying the numeric values that an arm sets: n_servers, and 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 increases quality-adjusted life-years by reducing the time patients spend waiting at the low utility, at a higher cost.

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

	cost	effect	inc_cost	inc_effect	icer	status
intervention
Standard capacity	11314.54	2.27	NaN	NaN	NaN	ND
Expanded capacity	12700.84	2.35	1386.3	0.07	18690.46	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, so it 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

The single specialist leaves patients waiting weeks; the second clears the queue to about a day. See the other model engine tutorials for the cohort state-transition, microsimulation, and compartmental engines this one sits alongside.

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`, the model type to use when a shared, scarce resource makes patients queue, so that one patient's outcome depends on the others in the system, a coupling neither a cohort nor an independent microsimulation model can represent. It runs that model through the same cost-effectiveness and value-of-information analysis as any other model; the [microsimulation tutorial](microsim.qmd) covers 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 stay the user's own 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, so billing the queueing time a patient just endured is one call once the specialist is granted. The process reads the run's time horizon back as `toolkit.horizon`, so it need not repeat the value the engine already holds. A patient still waiting at the horizon is never seen: the horizon arrives before the specialist, and only the waiting segment is billed. ## Configuring once, evaluating on draws The model goes to the engine when you build it; the draw matrix goes to `run_psa`. The two interventions differ only in capacity, so each is an `Intervention` carrying the numeric values that an arm sets: `n_servers`, and 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 increases quality-adjusted life-years by reducing the time patients spend waiting at the low utility, at a higher cost. ```{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, so it 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) ``` The single specialist leaves patients waiting weeks; the second clears the queue to about a day. See the other [model engine tutorials](index.qmd) for the cohort state-transition, microsimulation, and compartmental engines this one sits alongside. Next: the [full pipeline tutorial](full-pipeline.qmd) walks through the whole workflow, from `ParameterSet` sampling through `run_psa` to the analysis functions.