Discrete-event simulation engine

This tutorial shows how to build a resource-constrained model with DESModel, the model type to use when a shared, scarce resource makes entities 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 carry. 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 over time and wait for a specialist, a SimPy Resource. While waiting, untreated disease costs c_wait_year per year at utility u_wait. Once seen, a one-off 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 on top.

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 request loses the race 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 knobs that arm sets: n_servers, and a per-patient overhead c_capacity for the expanded arm. This is the right use of decision levers, a scenario knob the model already reads as a parameter, rather than a flag standing in for which arm is which. They see 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 cutting 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 turns it into per-request waits without touching engine internals.

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 engines concept page for how this model type compares with the others.

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