import numpy as np
import pandas as pd
import simpy
from sim_tools.distributions import ExponentialPerformance measures
Learning objectives:
- Learn how to record various performance measures within the simulation model.
- Understand how time-weighted calculations work.
Relevant reproducibility guidelines:
- STARS Reproducibility Recommendations (⭐): Include code to calculate all required model outputs
- STARS Reproducibility Recommendations: Ensure clarity and consistency in the model results tables.
- NHS Levels of RAP (🥈): Data is handled and output in a Tidy data format.
Pre-reading:
This page continues on from: Initialisation bias.
Entity generation → Entity processing → Initialisation bias → Performance measures
Required packages:
These should be available from environment setup in the “Test yourself” section of Environments.
Where relevant, we have included corresponding terms from queueing theory for each measure.
What is queueing theory? Queueing theory is a field of mathematics that analyses queues under strict assumptions, and is suitable for simple well-characterised systems. Discrete event simulation (DES) can be applied in similar ways (e.g., replicating classic queueing models like M/M/s), but more often extends beyond simple theory to handle a wider range of arrival and service patterns and more complex resource constraints.
Preparation for results collection
Before we start calculating any simulation performance measures, some set-up is required. This includes initial changes to Parameters and Model, and the creation of a new class: Runner.
Parameter class
Set verbose to False in the Parameters class. This lets the workflow focus on the calculated results rather than individual patient activities.
class Parameters:
"""
Parameter class.
Attributes
----------
interarrival_time : float
Mean time between arrivals (minutes).
consultation_time : float
Mean length of doctor's consultation (minutes).
number_of_doctors : int
Number of doctors.
warm_up_period : int
Duration of the warm-up period (minutes).
data_collection_period : int
Duration of the data collection period (minutes).
verbose : bool
Whether to print messages as simulation runs.
"""
def __init__(
self, interarrival_time=5, consultation_time=10,
number_of_doctors=3, warm_up_period=30, data_collection_period=40,
verbose=False
):
"""
Initialise Parameters instance.
Parameters
----------
interarrival_time : float
Time between arrivals (minutes).
consultation_time : float
Length of consultation (minutes).
number_of_doctors : int
Number of doctors.
warm_up_period : int
Duration of the warm-up period (minutes).
data_collection_period : int
Duration of the data collection period (minutes).
verbose : bool
Whether to print messages as simulation runs.
"""
self.interarrival_time = interarrival_time
self.consultation_time = consultation_time
self.number_of_doctors = number_of_doctors
self.warm_up_period = warm_up_period
self.data_collection_period = data_collection_period
self.verbose = verbosePatient class
No changes required for Patient at this stage.
Model class
We need to update the Model class so its run method returns a list of patient attributes. This makes patient-level data easy to analyse after simulation runs.
class Model:
"""
Simulation model.
Attributes
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
env : simpy.Environment
The SimPy environment for the simulation.
doctor : simpy.Resource
SimPy resource representing doctors.
patients : list
List of Patient objects.
results_list : list
List of dictionaries with the attributes of each patient.
arrival_dist : Exponential
Distribution used to generate random patient inter-arrival times.
consult_dist : Exponential
Distribution used to generate length of a doctor's consultation.
"""
def __init__(self, param, run_number):
"""
Create a new Model instance.
Parameters
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
"""
self.param = param
self.run_number = run_number
# Create SimPy environment
self.env = simpy.Environment()
# Create resource
self.doctor = simpy.Resource(
self.env, capacity=self.param.number_of_doctors
)
# Create a random seed sequence based on the run number
ss = np.random.SeedSequence(self.run_number)
seeds = ss.spawn(2)
# Set up attributes to store results
self.patients = []
self.results_list = []
# Initialise distributions
self.arrival_dist = Exponential(mean=self.param.interarrival_time,
random_seed=seeds[0])
self.consult_dist = Exponential(mean=self.param.consultation_time,
random_seed=seeds[1])
def generate_arrivals(self):
"""
Process that generates patient arrivals.
"""
while True:
# Sample and pass time to next arrival
sampled_iat = self.arrival_dist.sample()
yield self.env.timeout(sampled_iat)
# Check whether arrived during warm-up or data collection
if self.env.now < self.param.warm_up_period:
period = "\U0001F538 WU"
else:
period = "\U0001F539 DC"
# Create a new patient
patient = Patient(patient_id=len(self.patients)+1,
period=period,
arrival_time=self.env.now)
self.patients.append(patient)
# Print arrival time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"arrives at time: {patient.arrival_time:.3f}")
# Start process of consultation
self.env.process(self.consultation(patient))
def consultation(self, patient):
"""
Process that simulates a consultation.
Parameters
----------
patient :
Instance of the Patient() class representing a single patient.
"""
# Patient requests access to a doctor (resource)
with self.doctor.request() as req:
yield req
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"starts consultation at: {self.env.now:.3f}")
# Sample consultation duration and pass time spent with doctor
time_with_doctor = self.consult_dist.sample()
yield self.env.timeout(time_with_doctor)
# Record end time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"leaves at: {self.env.now:.3f}")
def reset_results(self):
"""
Reset results.
"""
self.patients = []
def warmup(self):
"""
Reset result collection after the warm-up period.
"""
if self.param.warm_up_period > 0:
# Delay process until warm-up period has completed
yield self.env.timeout(self.param.warm_up_period)
# Reset results variables
self.reset_results()
if self.param.verbose:
print(f"Warm up period ended at time: {self.env.now}")
def run(self):
"""
Run the simulation.
"""
# Schedule arrival generator and warm-up
self.env.process(self.generate_arrivals())
self.env.process(self.warmup())
# Run the simulation
self.env.run(until=(self.param.warm_up_period +
self.param.data_collection_period))
# Create list of dictionaries containing each patient's attributes
self.results_list = [x.__dict__ for x in self.patients]Runner class
We’re also introducing a new Runner class, responsible for running simulations and performing results calculations. This keeps calculation logic separate from the core modelling code. Also, this class will be modified to perform multiple model runs on the Replications page.
The basic structure of Runner is as follows. It initialises with simulation parameters, runs the model for a single replication, and organises results into separate patient-level and run-level outputs.
class Runner:
"""
Run the simulation.
Attributes
----------
param : Parameters
Simulation parameters.
"""
def __init__(self, param):
"""
Initialise a new instance of the Runner class.
Parameters
----------
param : Parameters
Simulation parameters.
"""
self.param = param
def run_single(self, run):
"""
Runs the simulation once and performs results calculations.
Parameters
----------
run : int
Run number for the simulation.
Returns
-------
dict
Contains patient-level results and results from each run.
"""
model = Model(param=self.param, run_number=run)
model.run()
# Patient results
patient_results = pd.DataFrame(model.results_list)
patient_results["run"] = run
# Run results
run_results = {
"run_number": run
}
return {
"patient": patient_results,
"run": run_results
}Running the model via Runner produces our new structured results.
patient_id period arrival_time run
0 1 🔹 DC 37.778250 0
1 2 🔹 DC 38.108184 0
2 3 🔹 DC 42.610666 0
3 4 🔹 DC 44.087985 0
4 5 🔹 DC 55.587763 0
5 6 🔹 DC 64.880007 0
6 7 🔹 DC 64.954293 0{'run_number': 0}Now we’re ready to record some performance measures. 😎
Before we start calculating any simulation performance measures, some set-up is required. This includes making a new function get_run_results(), and calling that and filter_warmup() from within model().
Parameter function
No changes required.
Warm-up function
No changes required.
Run results function
Let’s create a new function get_run_results(). This will be used to combine the performance measures we calculate into a single table.
#' Get average results for the provided single run.
#'
#' @param results Named list with `arrivals` from `get_mon_arrivals()` and
#' `resources` from `get_mon_resources()`
#' (`per_resource = TRUE` and `ongoing = TRUE`).
#' @param run_number Integer index of current simulation run.
#'
#' @importFrom dplyr bind_cols
#' @importFrom tibble tibble
#'
#' @return Tibble with processed results from replication.
#' @export
get_run_results <- function(results, run_number) {
metrics <- list()
dplyr::bind_cols(tibble(replication = run_number), metrics)
}Model function
Having introduced the filter_warmup() function on the Initialisation bias page and get_run_results() above, we’ll now call these from within model().
We also need to correct the “replication” column. This is generated by get_mon_arrivals() and get_mon_resources(), and will have been set to 1. This is because these functions assume, when they haven’t been supplied with a list of environments, that there was one replication. We need to set it to the correct run_number.
#' Run simulation.
#'
#' @param param List. Model parameters.
#' @param run_number Numeric. Run number for random seed generation.
#'
#' @importFrom dplyr mutate
#' @importFrom simmer add_generator add_resource get_mon_arrivals
#' @importFrom simmer get_mon_resources release run seize simmer timeout
#' @importFrom simmer trajectory
#'
#' @return Named list with tables `arrivals`, `resources` and `run_results`.
#' @export
model <- function(param, run_number) {
# Set random seed based on run number
set.seed(run_number)
# Create simmer environment
env <- simmer("simulation", verbose = param[["verbose"]])
# Define the patient trajectory
patient <- trajectory("consultation") |>
seize("doctor", 1L) |>
timeout(function() {
rexp(n = 1L, rate = 1L / param[["consultation_time"]])
}) |>
release("doctor", 1L)
env <- env |>
# Add doctor resource
add_resource("doctor", param[["number_of_doctors"]]) |>
# Add patient generator
add_generator("patient", patient, function() {
rexp(n = 1L, rate = 1L / param[["interarrival_time"]])
}) |>
# Run the simulation
simmer::run(until = (param[["warm_up_period"]] +
param[["data_collection_period"]]))
# Extract information on arrivals and resources from simmer environment
result <- list(
arrivals = get_mon_arrivals(env, per_resource = TRUE, ongoing = TRUE),
resources = get_mon_resources(env)
)
# Filter to remove results from the warm-up period
result <- filter_warmup(result, param[["warm_up_period"]])
# Replace "replication" value with appropriate run number
result[["arrivals"]] <- mutate(result[["arrivals"]],
replication = run_number)
result[["resources"]] <- mutate(result[["resources"]],
replication = run_number)
# Calculate the average results for that run
result[["run_results"]] <- get_run_results(result, run_number)
result
}Run the model
| name | start_time | end_time | activity_time | resource | replication |
|---|---|---|---|---|---|
| patient6 | 32.10174 | 44.96905 | 10.545432 | doctor | 1 |
| patient7 | 54.22141 | NA | NA | doctor | 1 |
| patient8 | 59.39763 | 62.76696 | 3.369335 | doctor | 1 |
| patient9 | 62.67136 | NA | NA | doctor | 1 |
| patient10 | 65.61376 | 68.55496 | 2.941204 | doctor | 1 |
| patient11 | 68.82322 | 69.88395 | 1.060726 | doctor | 1 |
| resource | time | server | queue | capacity | queue_size | system | limit | replication |
|---|---|---|---|---|---|---|---|---|
| doctor | 30.00000 | 3 | 0 | 3 | Inf | 3 | Inf | 1 |
| doctor | 32.10174 | 3 | 1 | 3 | Inf | 4 | Inf | 1 |
| doctor | 34.42362 | 3 | 0 | 3 | Inf | 3 | Inf | 1 |
| doctor | 40.66762 | 2 | 0 | 3 | Inf | 2 | Inf | 1 |
| doctor | 41.46371 | 1 | 0 | 3 | Inf | 1 | Inf | 1 |
| doctor | 44.96905 | 0 | 0 | 3 | Inf | 0 | Inf | 1 |
| doctor | 54.22141 | 1 | 0 | 3 | Inf | 1 | Inf | 1 |
| doctor | 59.39763 | 2 | 0 | 3 | Inf | 2 | Inf | 1 |
| doctor | 62.67136 | 3 | 0 | 3 | Inf | 3 | Inf | 1 |
| doctor | 62.76696 | 2 | 0 | 3 | Inf | 2 | Inf | 1 |
| doctor | 65.61376 | 3 | 0 | 3 | Inf | 3 | Inf | 1 |
| doctor | 68.55496 | 2 | 0 | 3 | Inf | 2 | Inf | 1 |
| doctor | 68.82322 | 3 | 0 | 3 | Inf | 3 | Inf | 1 |
| doctor | 69.88395 | 2 | 0 | 3 | Inf | 2 | Inf | 1 |
| replication |
|---|
| 1 |
At the moment, result[["run_results"]] is pretty empty (just has the replication number). Let’s add some performance measures!
Total arrivals
None▶ View how to record the total arrivals
Recording the total number of arrivals is super simple! We just need to find the length of the patient list.
Patient class
No changes required.
Model class
No changes required.
Runner class
class Runner:
"""
Run the simulation.
Attributes
----------
param : Parameters
Simulation parameters.
"""
def __init__(self, param):
"""
Initialise a new instance of the Runner class.
Parameters
----------
param : Parameters
Simulation parameters.
"""
self.param = param
def run_single(self, run):
"""
Runs the simulation once and performs results calculations.
Parameters
----------
run : int
Run number for the simulation.
Returns
-------
dict
Contains patient-level results and results from each run.
"""
model = Model(param=self.param, run_number=run)
model.run()
# Patient results
patient_results = pd.DataFrame(model.results_list)
patient_results["run"] = run
# Run results
run_results = {
"run_number": run,
"arrivals": len(patient_results)
}
return {
"patient": patient_results,
"run": run_results
}Run the model
patient_id period arrival_time run
0 1 🔹 DC 37.778250 0
1 2 🔹 DC 38.108184 0
2 3 🔹 DC 42.610666 0
3 4 🔹 DC 44.087985 0
4 5 🔹 DC 55.587763 0
5 6 🔹 DC 64.880007 0
6 7 🔹 DC 64.954293 0{'run_number': 0, 'arrivals': 7}Arrivals function
We introduce a new function calc_arrivals(). This determines the number of arrivals by counting the number of unique patient names in the arrival records.
#' Calculate the number of arrivals
#'
#' @param arrivals Dataframe with times for each patient with each resource.
#'
#' @importFrom dplyr n_distinct summarise ungroup
#'
#' @return Tibble with column containing total number of arrivals.
#' @export
calc_arrivals <- function(arrivals) {
arrivals |>
summarise(arrivals = n_distinct(.data[["name"]])) |>
ungroup()
}Run results function
In get_run_results(), the table results[["arrivals"]] is passed to calc_arrivals().
#' Get average results for the provided single run.
#'
#' @param results Named list with `arrivals` from `get_mon_arrivals()` and
#' `resources` from `get_mon_resources()`
#' (`per_resource = TRUE` and `ongoing = TRUE`).
#' @param run_number Integer index of current simulation run.
#'
#' @importFrom dplyr bind_cols
#' @importFrom tibble tibble
#'
#' @return Tibble with processed results from replication.
#' @export
get_run_results <- function(results, run_number) {
metrics <- list(
calc_arrivals(results[["arrivals"]])
)
dplyr::bind_cols(tibble(replication = run_number), metrics)
}Model function
No changes required.
Run the model
In arrivals we can see six patients listed, which aligns with our results table listing arrivals = 6.
| name | start_time | end_time | activity_time | resource | replication |
|---|---|---|---|---|---|
| patient6 | 32.10174 | 44.96905 | 10.545432 | doctor | 1 |
| patient7 | 54.22141 | NA | NA | doctor | 1 |
| patient8 | 59.39763 | 62.76696 | 3.369335 | doctor | 1 |
| patient9 | 62.67136 | NA | NA | doctor | 1 |
| patient10 | 65.61376 | 68.55496 | 2.941204 | doctor | 1 |
| patient11 | 68.82322 | 69.88395 | 1.060726 | doctor | 1 |
| replication | arrivals |
|---|---|
| 1 | 6 |
Mean wait time
Corresponding term in queueing theory: Average wait time in queue, \(\\W_q\)
None▶ View how to record the mean wait time
For this measure, we need to record the time each patient spent waiting for the doctor.
Patient class
A new attribute wait_time is add to the Patient class.
class Patient:
"""
Represents a patient.
Attributes
----------
patient_id : int
Unique patient identifier.
period : str
Arrival period (warm up or data collection) with emoji.
arrival_time : float
Time patient entered the system (minutes).
wait_time : float
Time spent waiting for the doctor (minutes).
"""
def __init__(self, patient_id, period, arrival_time):
"""
Initialises a new patient.
Parameters
----------
patient_id : int
Unique patient identifier.
period : str
Arrival period (warm up or data collection) with emoji.
arrival_time : float
Time patient entered the system (minutes).
"""
self.patient_id = patient_id
self.period = period
self.arrival_time = arrival_time
self.wait_time = np.nanModel class
Within Model.consultation(), the start time of the consultation is recorded (start_wait). When the resource request has been granted, the time elapsed is saved under patient.wait_time.
class Model:
"""
Simulation model.
Attributes
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
env : simpy.Environment
The SimPy environment for the simulation.
doctor : simpy.Resource
SimPy resource representing doctors.
patients : list
List of Patient objects.
results_list : list
List of dictionaries with the attributes of each patient.
arrival_dist : Exponential
Distribution used to generate random patient inter-arrival times.
consult_dist : Exponential
Distribution used to generate length of a doctor's consultation.
"""
def __init__(self, param, run_number):
"""
Create a new Model instance.
Parameters
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
"""
self.param = param
self.run_number = run_number
# Create SimPy environment
self.env = simpy.Environment()
# Create resource
self.doctor = simpy.Resource(
self.env, capacity=self.param.number_of_doctors
)
# Create a random seed sequence based on the run number
ss = np.random.SeedSequence(self.run_number)
seeds = ss.spawn(2)
# Set up attributes to store results
self.patients = []
self.results_list = []
# Initialise distributions
self.arrival_dist = Exponential(mean=self.param.interarrival_time,
random_seed=seeds[0])
self.consult_dist = Exponential(mean=self.param.consultation_time,
random_seed=seeds[1])
def generate_arrivals(self):
"""
Process that generates patient arrivals.
"""
while True:
# Sample and pass time to next arrival
sampled_iat = self.arrival_dist.sample()
yield self.env.timeout(sampled_iat)
# Check whether arrived during warm-up or data collection
if self.env.now < self.param.warm_up_period:
period = "\U0001F538 WU"
else:
period = "\U0001F539 DC"
# Create a new patient
patient = Patient(patient_id=len(self.patients)+1,
period=period,
arrival_time=self.env.now)
self.patients.append(patient)
# Print arrival time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"arrives at time: {patient.arrival_time:.3f}")
# Start process of consultation
self.env.process(self.consultation(patient))
def consultation(self, patient):
"""
Process that simulates a consultation.
Parameters
----------
patient :
Instance of the Patient() class representing a single patient.
"""
start_wait = self.env.now
# Patient requests access to a doctor (resource)
with self.doctor.request() as req:
yield req
# Record how long patient waited before consultation
patient.wait_time = self.env.now - start_wait
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} starts" +
f" consultation at: {self.env.now:.3f}")
# Sample consultation duration and pass time spent with doctor
time_with_doctor = self.consult_dist.sample()
yield self.env.timeout(time_with_doctor)
# Record end time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"leaves at: {self.env.now:.3f}")
def reset_results(self):
"""
Reset results.
"""
self.patients = []
def warmup(self):
"""
Reset result collection after the warm-up period.
"""
if self.param.warm_up_period > 0:
# Delay process until warm-up period has completed
yield self.env.timeout(self.param.warm_up_period)
# Reset results variables
self.reset_results()
if self.param.verbose:
print(f"Warm up period ended at time: {self.env.now}")
def run(self):
"""
Run the simulation.
"""
# Schedule arrival generator and warm-up
self.env.process(self.generate_arrivals())
self.env.process(self.warmup())
# Run the simulation
self.env.run(until=(self.param.warm_up_period +
self.param.data_collection_period))
# Create list of dictionaries containing each patient's attributes
self.results_list = [x.__dict__ for x in self.patients]Runner class
The mean wait_time from the patients is calculated.
class Runner:
"""
Run the simulation.
Attributes
----------
param : Parameters
Simulation parameters.
"""
def __init__(self, param):
"""
Initialise a new instance of the Runner class.
Parameters
----------
param : Parameters
Simulation parameters.
"""
self.param = param
def run_single(self, run):
"""
Runs the simulation once and performs results calculations.
Parameters
----------
run : int
Run number for the simulation.
Returns
-------
dict
Contains patient-level results and results from each run.
"""
model = Model(param=self.param, run_number=run)
model.run()
# Patient results
patient_results = pd.DataFrame(model.results_list)
patient_results["run"] = run
# Run results
run_results = {
"run_number": run,
"mean_wait_time": patient_results["wait_time"].mean()
}
return {
"patient": patient_results,
"run": run_results
}Run the model
We can see each patient’s wait time in result["patient"], and the mean wait time in result["run"].
patient_id period arrival_time wait_time run
0 1 🔹 DC 37.778250 0.000000 0
1 2 🔹 DC 38.108184 0.000000 0
2 3 🔹 DC 42.610666 4.987255 0
3 4 🔹 DC 44.087985 12.574406 0
4 5 🔹 DC 55.587763 1.800804 0
5 6 🔹 DC 64.880007 0.000000 0
6 7 🔹 DC 64.954293 NaN 0{'run_number': 0, 'mean_wait_time': np.float64(3.227077467801481)}Wait time function
A new calc_mean_wait() function removes patients still waiting (since their wait times are censored and including them would underestimate the mean), then calculates the mean wait time for each resource.
#' Calculate the mean wait time for each resource
#'
#' @param arrivals Dataframe with times for each patient with each resource.
#'
#' @importFrom dplyr group_by summarise ungroup
#' @importFrom rlang .data
#' @importFrom tidyr drop_na pivot_wider
#' @importFrom tidyselect any_of
#'
#' @return Tibble with columns containing result for each resource.
#' @export
calc_mean_wait <- function(arrivals) {
# Create subset of data that removes patients who were still waiting
complete_arrivals <- drop_na(arrivals, any_of("wait_time"))
# Calculate mean wait time for each resource
complete_arrivals |>
group_by(.data[["resource"]]) |>
summarise(mean_wait_time = mean(.data[["wait_time"]])) |>
pivot_wider(names_from = "resource",
values_from = "mean_wait_time",
names_glue = "mean_wait_time_{resource}") |>
ungroup()
}Run results function
In get_run_results(), the table results[["arrivals"]] is passed to calc_mean_wait().
#' Get average results for the provided single run.
#'
#' @param results Named list with `arrivals` from `get_mon_arrivals()` and
#' `resources` from `get_mon_resources()`
#' (`per_resource = TRUE` and `ongoing = TRUE`).
#' @param run_number Integer index of current simulation run.
#'
#' @importFrom dplyr bind_cols
#' @importFrom tibble tibble
#'
#' @return Tibble with processed results from replication.
#' @export
get_run_results <- function(results, run_number) {
metrics <- list(
calc_mean_wait(results[["arrivals"]])
)
dplyr::bind_cols(tibble(replication = run_number), metrics)
}Model function
In model(), we save the time the patient started their consultation (doctor_serve_start) using set_attribute().
To retrieve this later using simmer’s get_mon_attributes() function, we have to set mon = 2L within the patient generator, as this enables attribute monitoring beyond basic arrival tracking.
When retrieving the attribute, we restructure the returned dataframe using a generic approach that handles any attributes following the pattern resourcename_attribute - for example, transforming doctor_serve_start into serve_start (aligned with the corresponding doctor rows).
We can calculate the wait time for each patient by subtracting their arrival time (start_time) from when service began (serve_start).
#' Run simulation.
#'
#' @param param List. Model parameters.
#' @param run_number Numeric. Run number for random seed generation.
#'
#' @importFrom dplyr left_join mutate
#' @importFrom rlang .data
#' @importFrom simmer add_generator add_resource get_mon_arrivals
#' @importFrom simmer get_mon_resources now release run seize set_attribute
#' @importFrom simmer simmer timeout trajectory
#' @importFrom tidyr pivot_wider
#'
#' @return Named list with tables `arrivals`, `resources` and `run_results`.
#' @export
model <- function(param, run_number) {
# Set random seed based on run number
set.seed(run_number)
# Create simmer environment
env <- simmer("simulation", verbose = param[["verbose"]])
# Define the patient trajectory
patient <- trajectory("consultation") |>
seize("doctor", 1L) |>
# Record time resource is obtained
set_attribute("doctor_serve_start", function() now(env)) |>
timeout(function() {
rexp(n = 1L, rate = 1L / param[["consultation_time"]])
}) |>
release("doctor", 1L)
env <- env |>
# Add doctor resource
add_resource("doctor", param[["number_of_doctors"]]) |>
# Add patient generator (mon=2 so can get manual attributes)
add_generator("patient", patient, function() {
rexp(n = 1L, rate = 1L / param[["interarrival_time"]])
}, mon = 2L) |>
# Run the simulation
simmer::run(until = (param[["warm_up_period"]] +
param[["data_collection_period"]]))
# Extract information on arrivals and resources from simmer environment
result <- list(
arrivals = get_mon_arrivals(env, per_resource = TRUE, ongoing = TRUE),
resources = get_mon_resources(env)
)
# Get the extra arrivals attributes
extra_attributes <- get_mon_attributes(env) |>
dplyr::select("name", "key", "value") |>
# Add column with resource name, and remove that from key
mutate(resource = gsub("_.+", "", .data[["key"]]),
key = gsub("^[^_]+_", "", .data[["key"]])) |>
pivot_wider(names_from = "key", values_from = "value")
# Merge extra attributes with the arrival data
result[["arrivals"]] <- left_join(
result[["arrivals"]], extra_attributes, by = c("name", "resource")
)
# Filter to remove results from the warm-up period
result <- filter_warmup(result, param[["warm_up_period"]])
# Replace "replication" value with appropriate run number
result[["arrivals"]] <- mutate(result[["arrivals"]],
replication = run_number)
result[["resources"]] <- mutate(result[["resources"]],
replication = run_number)
# Calculate wait time for each patient
result[["arrivals"]] <- mutate(
result[["arrivals"]],
wait_time = .data[["serve_start"]] - .data[["start_time"]]
)
# Calculate the average results for that run
result[["run_results"]] <- get_run_results(result, run_number)
result
}Run the model
| name | start_time | end_time | activity_time | resource | replication | serve_start | wait_time |
|---|---|---|---|---|---|---|---|
| patient6 | 32.10174 | 44.96905 | 10.545432 | doctor | 1 | 34.42362 | 2.321881 |
| patient7 | 54.22141 | NA | NA | doctor | 1 | 54.22141 | 0.000000 |
| patient8 | 59.39763 | 62.76696 | 3.369335 | doctor | 1 | 59.39763 | 0.000000 |
| patient9 | 62.67136 | NA | NA | doctor | 1 | 62.67136 | 0.000000 |
| patient10 | 65.61376 | 68.55496 | 2.941204 | doctor | 1 | 65.61376 | 0.000000 |
| patient11 | 68.82322 | 69.88395 | 1.060726 | doctor | 1 | 68.82322 | 0.000000 |
| replication | mean_wait_time_doctor |
|---|---|
| 1 | 0.3869803 |
Mean time in consultation
Corresponding term in queueing theory: Service time, 1/\(\mu\)
None▶ View how to record the mean time spent in the consultation
Patient class
A new attribute time_with_doctor is add to the Patient class.
class Patient:
"""
Represents a patient.
Attributes
----------
patient_id : int
Unique patient identifier.
period : str
Arrival period (warm up or data collection) with emoji.
arrival_time : float
Time patient entered the system (minutes).
time_with_doctor : float
Time spent in consultation with a doctor (minutes).
"""
def __init__(self, patient_id, period, arrival_time):
"""
Initialises a new patient.
Parameters
----------
patient_id : int
Unique patient identifier.
period : str
Arrival period (warm up or data collection) with emoji.
arrival_time : float
Time patient entered the system (minutes).
"""
self.patient_id = patient_id
self.period = period
self.arrival_time = arrival_time
self.time_with_doctor = np.nanModel class
Within Model.consultation(), the sampled time with the doctor is saved under patient.time_with_doctor.
class Model:
"""
Simulation model.
Attributes
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
env : simpy.Environment
The SimPy environment for the simulation.
doctor : simpy.Resource
SimPy resource representing doctors.
patients : list
List of Patient objects.
results_list : list
List of dictionaries with the attributes of each patient.
arrival_dist : Exponential
Distribution used to generate random patient inter-arrival times.
consult_dist : Exponential
Distribution used to generate length of a doctor's consultation.
"""
def __init__(self, param, run_number):
"""
Create a new Model instance.
Parameters
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
"""
self.param = param
self.run_number = run_number
# Create SimPy environment
self.env = simpy.Environment()
# Create resource
self.doctor = simpy.Resource(
self.env, capacity=self.param.number_of_doctors
)
# Create a random seed sequence based on the run number
ss = np.random.SeedSequence(self.run_number)
seeds = ss.spawn(2)
# Set up attributes to store results
self.patients = []
self.results_list = []
# Initialise distributions
self.arrival_dist = Exponential(mean=self.param.interarrival_time,
random_seed=seeds[0])
self.consult_dist = Exponential(mean=self.param.consultation_time,
random_seed=seeds[1])
def generate_arrivals(self):
"""
Process that generates patient arrivals.
"""
while True:
# Sample and pass time to next arrival
sampled_iat = self.arrival_dist.sample()
yield self.env.timeout(sampled_iat)
# Check whether arrived during warm-up or data collection
if self.env.now < self.param.warm_up_period:
period = "\U0001F538 WU"
else:
period = "\U0001F539 DC"
# Create a new patient
patient = Patient(patient_id=len(self.patients)+1,
period=period,
arrival_time=self.env.now)
self.patients.append(patient)
# Print arrival time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"arrives at time: {patient.arrival_time:.3f}")
# Start process of consultation
self.env.process(self.consultation(patient))
def consultation(self, patient):
"""
Process that simulates a consultation.
Parameters
----------
patient :
Instance of the Patient() class representing a single patient.
"""
# Patient requests access to a doctor (resource)
with self.doctor.request() as req:
yield req
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} starts" +
f" consultation at: {self.env.now:.3f}")
# Sample consultation duration and pass time spent with doctor
patient.time_with_doctor = self.consult_dist.sample()
yield self.env.timeout(patient.time_with_doctor)
# Record end time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"leaves at: {self.env.now:.3f}")
def reset_results(self):
"""
Reset results.
"""
self.patients = []
def warmup(self):
"""
Reset result collection after the warm-up period.
"""
if self.param.warm_up_period > 0:
# Delay process until warm-up period has completed
yield self.env.timeout(self.param.warm_up_period)
# Reset results variables
self.reset_results()
if self.param.verbose:
print(f"Warm up period ended at time: {self.env.now}")
def run(self):
"""
Run the simulation.
"""
# Schedule arrival generator and warm-up
self.env.process(self.generate_arrivals())
self.env.process(self.warmup())
# Run the simulation
self.env.run(until=(self.param.warm_up_period +
self.param.data_collection_period))
# Create list of dictionaries containing each patient's attributes
self.results_list = [x.__dict__ for x in self.patients]Runner class
The mean time_with_doctor from the patients is calculated.
class Runner:
"""
Run the simulation.
Attributes
----------
param : Parameters
Simulation parameters.
"""
def __init__(self, param):
"""
Initialise a new instance of the Runner class.
Parameters
----------
param : Parameters
Simulation parameters.
"""
self.param = param
def run_single(self, run):
"""
Runs the simulation once and performs results calculations.
Parameters
----------
run : int
Run number for the simulation.
Returns
-------
dict
Contains patient-level results and results from each run.
"""
model = Model(param=self.param, run_number=run)
model.run()
# Patient results
patient_results = pd.DataFrame(model.results_list)
patient_results["run"] = run
# Run results
run_results = {
"run_number": run,
"mean_time_with_doctor": (
patient_results["time_with_doctor"].mean()
)
}
return {
"patient": patient_results,
"run": run_results
}Run the model
We can see each patient’s consultation length in result["patient"], and the mean time in result["run"].
patient_id period arrival_time time_with_doctor run
0 1 🔹 DC 37.778250 19.610317 0
1 2 🔹 DC 38.108184 9.489737 0
2 3 🔹 DC 42.610666 41.665475 0
3 4 🔹 DC 44.087985 5.874479 0
4 5 🔹 DC 55.587763 27.882444 0
5 6 🔹 DC 64.880007 24.914615 0
6 7 🔹 DC 64.954293 NaN 0{'run_number': 0, 'mean_time_with_doctor': np.float64(21.57284454497007)}Service length function
A new function calc_mean_serve_length() finds the mean length of time patients spent with each resource using a new column we’ll introduce in model() called serve_length.
#' Calculate the mean length of time patients spent with each resource
#'
#' @param arrivals Dataframe with times for each patient with each resource.
#' @param resources Dataframe with times patients use or queue for resources.
#' @param groups Optional list of columns to group by for the calculation.
#'
#' @importFrom dplyr group_by summarise ungroup
#' @importFrom rlang .data
#' @importFrom tidyr drop_na pivot_wider
#' @importFrom tidyselect any_of
#'
#' @return Tibble with columns containing result for each resource.
#' @export
calc_mean_serve_length <- function(arrivals, resources, groups = NULL) {
# Create subset of data that removes patients who were still waiting
complete_arrivals <- drop_na(arrivals, any_of("serve_length"))
# Calculate mean serve time for each resource
complete_arrivals |>
group_by(.data[["resource"]]) |>
summarise(mean_serve_time = mean(.data[["serve_length"]])) |>
pivot_wider(names_from = "resource",
values_from = "mean_serve_time",
names_glue = "mean_time_with_{resource}") |>
ungroup()
}Run results function
In get_run_results(), the table results[["arrivals"]] is passed to calc_mean_serve_length().
#' Get average results for the provided single run.
#'
#' @param results Named list with `arrivals` from `get_mon_arrivals()` and
#' `resources` from `get_mon_resources()`
#' (`per_resource = TRUE` and `ongoing = TRUE`).
#' @param run_number Integer index of current simulation run.
#'
#' @importFrom dplyr bind_cols
#' @importFrom tibble tibble
#'
#' @return Tibble with processed results from replication.
#' @export
get_run_results <- function(results, run_number) {
metrics <- list(
calc_mean_serve_length(results[["arrivals"]])
)
dplyr::bind_cols(tibble(replication = run_number), metrics)
}Model function
In model(), we save the sampled length of the consultation (doctor_serve_length) using set_attribute(). This is then used for timeout(), retrieving it again with get_attribute().
Next, we repeat the steps we did for calculating mean wait time: setting mon = 2L in the patient generator, retrieving results using get_mon_attributes(), and restructuring the returned dataframe to transform doctor_serve_length into serve_length (aligned with the corresponding doctor rows).
#' Run simulation.
#'
#' @param param List. Model parameters.
#' @param run_number Numeric. Run number for random seed generation.
#'
#' @importFrom dplyr left_join mutate
#' @importFrom rlang .data
#' @importFrom simmer add_generator add_resource get_mon_arrivals
#' @importFrom simmer get_attribute get_mon_resources release run seize
#' @importFrom simmer set_attribute simmer timeout trajectory
#' @importFrom tidyr pivot_wider
#'
#' @return Named list with tables `arrivals`, `resources` and `run_results`.
#' @export
model <- function(param, run_number) {
# Set random seed based on run number
set.seed(run_number)
# Create simmer environment
env <- simmer("simulation", verbose = param[["verbose"]])
# Define the patient trajectory
patient <- trajectory("consultation") |>
seize("doctor", 1L) |>
# Record sampled length of consultation
set_attribute("doctor_serve_length", function() {
rexp(n = 1L, rate = 1L / param[["consultation_time"]])
}) |>
timeout(function() get_attribute(env, "doctor_serve_length")) |>
release("doctor", 1L)
env <- env |>
# Add doctor resource
add_resource("doctor", param[["number_of_doctors"]]) |>
# Add patient generator (mon=2 so can get manual attributes)
add_generator("patient", patient, function() {
rexp(n = 1L, rate = 1L / param[["interarrival_time"]])
}, mon = 2L) |>
# Run the simulation
simmer::run(until = (param[["warm_up_period"]] +
param[["data_collection_period"]]))
# Extract information on arrivals and resources from simmer environment
result <- list(
arrivals = get_mon_arrivals(env, per_resource = TRUE, ongoing = TRUE),
resources = get_mon_resources(env)
)
# Get the extra arrivals attributes
extra_attributes <- get_mon_attributes(env) |>
dplyr::select("name", "key", "value") |>
# Add column with resource name, and remove that from key
mutate(resource = gsub("_.+", "", .data[["key"]]),
key = gsub("^[^_]+_", "", .data[["key"]])) |>
pivot_wider(names_from = "key", values_from = "value")
# Merge extra attributes with the arrival data
result[["arrivals"]] <- left_join(
result[["arrivals"]], extra_attributes, by = c("name", "resource")
)
# Filter to remove results from the warm-up period
result <- filter_warmup(result, param[["warm_up_period"]])
# Replace "replication" value with appropriate run number
result[["arrivals"]] <- mutate(result[["arrivals"]],
replication = run_number)
result[["resources"]] <- mutate(result[["resources"]],
replication = run_number)
# Calculate the average results for that run
result[["run_results"]] <- get_run_results(result, run_number)
result
}Run the model
Why do the patient results look different? This measure requires saving sampled times using set_attribute() ahead of resource use. This changes the order in which random numbers are drawn from the simulation’s global random number generator. As a result, even with the same seed and parameters, the results will differ from previous runs.
| name | start_time | end_time | activity_time | resource | replication | serve_length |
|---|---|---|---|---|---|---|
| patient6 | 30.75713 | NA | NA | doctor | 1 | 44.239342 |
| patient7 | 36.02984 | 46.38228 | 10.352439 | doctor | 1 | 10.352439 |
| patient8 | 45.41002 | 51.95749 | 6.547466 | doctor | 1 | 6.547466 |
| patient9 | 47.09469 | 52.97948 | 5.884797 | doctor | 1 | 5.884797 |
| patient10 | 58.91726 | 65.33619 | 6.418926 | doctor | 1 | 6.418926 |
| patient11 | 60.38786 | 66.04652 | 5.658655 | doctor | 1 | 5.658655 |
| patient12 | 60.91823 | NA | NA | doctor | 1 | 11.733121 |
| patient13 | 61.21542 | NA | NA | doctor | 1 | 9.968130 |
| patient14 | 64.10899 | NA | NA | doctor | 1 | NA |
| replication | mean_time_with_doctor |
|---|---|
| 1 | 12.60036 |
Mean doctor utilisation
Corresponding term in queueing theory: Server utilisation, \(\rho\)
Mean utilisation measures the proportion of available doctor time that is actually used during the data collection period. We have suggested two different approaches for measuring utilisation - both return the same result.
The first approach is easier to understand, as it just sums consultation times for each patient. However, it relies on making corrections when consultations overlap the end of the simulation or span the warm-up and data collection periods.
None▶ View how to record utilisation
Patient class
No changes required.
Model class
A new attribute doctor_time_used is used to record the total time patients spend with doctor. This value is reset at the end of any warm-up period, so that only data from the main collection phase is analysed.
Recording the time used requires two corrections:
End of simulation overrun. If a consultation begins during the data collection period and extends past the simulation end, only the time up to simulation end is included. This prevents overestimating usage by excluding extra time beyond the simulation window.
Warm-up overrun. For models with a warm-up, consultations that begin during warm-up but finish during data collection would be missed after the reset, leading to underestimated usage. To correct for this, any time overlapping data collection (for these patients) is recorded in the doctor_time_used_correction attribute and added back in after the warm-up reset.
class Model:
"""
Simulation model.
Attributes
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
env : simpy.Environment
The SimPy environment for the simulation.
doctor : simpy.Resource
SimPy resource representing doctors.
patients : list
List of Patient objects.
doctor_time_used : float
Total time that doctor resources were used for (minutes).
doctor_time_used_correction : float
Adjustment for doctor time. Without this, usage is underestimated
since patients whose consultations began in warm-up but ended
during data collection are excluded.
results_list : list
List of dictionaries with the attributes of each patient.
arrival_dist : Exponential
Distribution used to generate random patient inter-arrival times.
consult_dist : Exponential
Distribution used to generate length of a doctor's consultation.
"""
def __init__(self, param, run_number):
"""
Create a new Model instance.
Parameters
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
"""
self.param = param
self.run_number = run_number
# Create SimPy environment
self.env = simpy.Environment()
# Create resource
self.doctor = simpy.Resource(
self.env, capacity=self.param.number_of_doctors
)
# Create a random seed sequence based on the run number
ss = np.random.SeedSequence(self.run_number)
seeds = ss.spawn(2)
# Set up attributes to store results
self.patients = []
self.results_list = []
self.doctor_time_used = 0
self.doctor_time_used_correction = 0
# Initialise distributions
self.arrival_dist = Exponential(mean=self.param.interarrival_time,
random_seed=seeds[0])
self.consult_dist = Exponential(mean=self.param.consultation_time,
random_seed=seeds[1])
def generate_arrivals(self):
"""
Process that generates patient arrivals.
"""
while True:
# Sample and pass time to next arrival
sampled_iat = self.arrival_dist.sample()
yield self.env.timeout(sampled_iat)
# Check whether arrived during warm-up or data collection
if self.env.now < self.param.warm_up_period:
period = "\U0001F538 WU"
else:
period = "\U0001F539 DC"
# Create a new patient
patient = Patient(patient_id=len(self.patients)+1,
period=period,
arrival_time=self.env.now)
self.patients.append(patient)
# Print arrival time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"arrives at time: {patient.arrival_time:.3f}")
# Start process of consultation
self.env.process(self.consultation(patient))
def consultation(self, patient):
"""
Process that simulates a consultation.
Parameters
----------
patient :
Instance of the Patient() class representing a single patient.
"""
# Patient requests access to a doctor (resource)
with self.doctor.request() as req:
yield req
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"starts consultation at: {self.env.now:.3f}")
# Sample consultation duration
time_with_doctor = self.consult_dist.sample()
# Add to total doctor time used
# If it runs past simulation end, only count the time until end
remaining_time = (
self.param.warm_up_period +
self.param.data_collection_period) - self.env.now
self.doctor_time_used += min(
time_with_doctor, remaining_time)
# During warm-up: check if consultation continues past warm-up.
# If so, record the portion overlapping with data collection in
# doctor_time_used_correction (capped at the simulation end).
remaining_warmup = self.param.warm_up_period - self.env.now
if remaining_warmup > 0:
time_exceeding_warmup = time_with_doctor - remaining_warmup
if time_exceeding_warmup > 0:
self.doctor_time_used_correction += min(
time_exceeding_warmup,
self.param.data_collection_period)
# Pass time spent with the doctor
yield self.env.timeout(time_with_doctor)
# Record end time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"leaves at: {self.env.now:.3f}")
def reset_results(self):
"""
Reset results.
"""
self.patients = []
self.doctor_time_used = 0
def warmup(self):
"""
Reset result collection after the warm-up period.
"""
if self.param.warm_up_period > 0:
# Delay process until warm-up period has completed
yield self.env.timeout(self.param.warm_up_period)
# Reset results variables
self.reset_results()
if self.param.verbose:
print(f"Warm up period ended at time: {self.env.now}")
# Add correction for patients whose consultations began in
# warm-up but continued into data collection.
self.doctor_time_used += self.doctor_time_used_correction
def run(self):
"""
Run the simulation.
"""
# Schedule arrival generator and warm-up
self.env.process(self.generate_arrivals())
self.env.process(self.warmup())
# Run the simulation
self.env.run(until=(self.param.warm_up_period +
self.param.data_collection_period))
# Create list of dictionaries containing each patient's attributes
self.results_list = [x.__dict__ for x in self.patients]Runner class
Mean utilisation can now be calculated as:
\[ \text{mean utilisation} = \frac{\text{doctor\_time\_used}}{\text{number of doctors} \times \text{data collection period}} \]
class Runner:
"""
Run the simulation.
Attributes
----------
param : Parameters
Simulation parameters.
"""
def __init__(self, param):
"""
Initialise a new instance of the Runner class.
Parameters
----------
param : Parameters
Simulation parameters.
"""
self.param = param
def run_single(self, run):
"""
Runs the simulation once and performs results calculations.
Parameters
----------
run : int
Run number for the simulation.
Returns
-------
dict
Contains patient-level results and results from each run.
"""
model = Model(param=self.param, run_number=run)
model.run()
# Patient results
patient_results = pd.DataFrame(model.results_list)
patient_results["run"] = run
# Run results
run_results = {
"run_number": run,
"mean_utilisation": model.doctor_time_used / (
self.param.number_of_doctors *
self.param.data_collection_period
)
}
return {
"patient": patient_results,
"run": run_results
}Run the model
patient_id period arrival_time run
0 1 🔹 DC 37.778250 0
1 2 🔹 DC 38.108184 0
2 3 🔹 DC 42.610666 0
3 4 🔹 DC 44.087985 0
4 5 🔹 DC 55.587763 0
5 6 🔹 DC 64.880007 0
6 7 🔹 DC 64.954293 0{'run_number': 0, 'mean_utilisation': 0.8743911437929205}While this method introduces a new class and may initially seem harder to grasp, it has the advantage of requiring fewer changes to the rest of the simulation code. It is also handy for choosing the length of the warm-up period, as it provides you with utilisation at each timepoint (rather than a final overall utilisation).
None▶ View how to record utilisation (time-weighted approach)
Patient class
No changes required.
Monitored resource class
For this approach, we create a new class called MonitoredResource. This is explained line-by-line below.
class MonitoredResource(simpy.Resource):
"""
SimPy Resource subclass that tracks resource utilisation.
Attributes
----------
time_last_event : list
Time of the most recent resource request or release event.
area_resource_busy : list
Cumulative time resources were busy, for computing utilisation.
Notes
-----
Class adapted from Monks, Harper and Heather 2025.
"""
def __init__(self, *args, **kwargs):
"""
Initialise the MonitoredResource, resetting monitoring statistics.
Parameters
----------
*args :
Positional arguments to be passed to the parent class.
**kwargs :
Keyword arguments to be passed to the parent class.
"""
# Initialise a SimPy Resource
super().__init__(*args, **kwargs)
# Run the init_results() method
self.init_results()
def init_results(self):
"""
Reset monitoring attributes.
"""
self.time_last_event = [self._env.now]
self.area_resource_busy = [0.0]
def request(self, *args, **kwargs):
"""
Update time-weighted statistics before requesting a resource.
Parameters
----------
*args :
Positional arguments to be passed to the parent class.
**kwargs :
Keyword arguments to be passed to the parent class.
Returns
-------
simpy.events.Event
Event for the resource request.
"""
# Update time-weighted statistics
self.update_time_weighted_stats()
# Request a resource
return super().request(*args, **kwargs)
def release(self, *args, **kwargs):
"""
Update time-weighted statistics before releasing a resource.
Parameters
----------
*args :
Positional arguments to be passed to the parent class.
**kwargs :
Keyword arguments to be passed to the parent class.
Returns
-------
simpy.events.Event
Event for the resource release.
"""
# Update time-weighted statistics
self.update_time_weighted_stats()
# Release a resource
return super().release(*args, **kwargs)
def update_time_weighted_stats(self):
"""
Update the time-weighted statistics for resource usage.
Notes
-----
- These sums can be referred to as "the area under the curve".
- They are called "time-weighted" statistics as they account for how
long certain events or states persist over time.
"""
# Calculate time since last event
time_since_last_event = self._env.now - self.time_last_event[-1]
# Add record of current time
self.time_last_event.append(self._env.now)
# Add "area under curve" of resources in use
# self.count is the number of resources in use
self.area_resource_busy.append(self.count * time_since_last_event)Explaining MonitoredResource…
class MonitoredResource(simpy.Resource):
"""
SimPy Resource subclass that tracks resource utilisation.
Attributes
----------
time_last_event : list
Time of the most recent resource request or release event.
area_resource_busy : list
Cumulative time resources were busy, for computing utilisation.
Notes
-----
Class adapted from Monks, Harper and Heather 2025.
"""
def __init__(self, *args, **kwargs):
"""
Initialise the MonitoredResource, resetting monitoring statistics.
Parameters
----------
*args :
Positional arguments to be passed to the parent class.
**kwargs :
Keyword arguments to be passed to the parent class.
"""
# Initialise a SimPy Resource
super().__init__(*args, **kwargs)
# Run the init_results() method
self.init_results()This class inherits SimPy’s Resource class (if new to inheritance, see the Code organisation page for an introduction). This means it has same functionality, but we can modify, etcetc.
In the __init__ method, the class is initialised as it would be for a standard resource class, but then also calls the new init_results method.
This method resets all tracking attributes. This is called when the object is created or when usage statistics need to be cleared, such as after a warm-up period.
def request(self, *args, **kwargs):
"""
Update time-weighted statistics before requesting a resource.
Parameters
----------
*args :
Positional arguments to be passed to the parent class.
**kwargs :
Keyword arguments to be passed to the parent class.
Returns
-------
simpy.events.Event
Event for the resource request.
"""
# Update time-weighted statistics
self.update_time_weighted_stats()
# Request a resource
return super().request(*args, **kwargs)Whenever a resource is requested, this method first updates time-weighted statistics based on the current system state. It then makes the actual resource request via the parent class, ensuring all tracking is up to date.
def release(self, *args, **kwargs):
"""
Update time-weighted statistics before releasing a resource.
Parameters
----------
*args :
Positional arguments to be passed to the parent class.
**kwargs :
Keyword arguments to be passed to the parent class.
Returns
-------
simpy.events.Event
Event for the resource release.
"""
# Update time-weighted statistics
self.update_time_weighted_stats()
# Release a resource
return super().release(*args, **kwargs)Just as we did with the request() method, the release() method updates usage statistics before a resource if released.
def update_time_weighted_stats(self):
"""
Update the time-weighted statistics for resource usage.
Notes
-----
- These sums can be referred to as "the area under the curve".
- They are called "time-weighted" statistics as they account for how
long certain events or states persist over time.
"""
# Calculate time since last event
time_since_last_event = self._env.now - self.time_last_event[-1]
# Add record of current time
self.time_last_event.append(self._env.now)
# Add "area under curve" of resources in use
# self.count is the number of resources in use
self.area_resource_busy.append(self.count * time_since_last_event)Calculates and records usage statistics since the last resource event (request or release). It determines how long the resource usage lasted, then adds its “area under the curve” to the running total.
Model class
In Model, we make three changes:
The
doctorattribute is now set up usingMonitoredResource(instead ofsimpy.Resource).In
reset_results(), we call theMonitoredResourcemethodinit_results()to reset their monitoring at the end of a warm-up period.In
run(), once the simulation ends, theMonitoredResourcemethodupdate_time_weighted_stats()is called one final time, to update with the usage between the last record resource event (request or release) and the simulation end.
class Model:
"""
Simulation model.
Attributes
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
env : simpy.Environment
The SimPy environment for the simulation.
doctor : simpy.Resource
SimPy resource representing doctors.
patients : list
List of Patient objects.
results_list : list
List of dictionaries with the attributes of each patient.
arrival_dist : Exponential
Distribution used to generate random patient inter-arrival times.
consult_dist : Exponential
Distribution used to generate length of a doctor's consultation.
"""
def __init__(self, param, run_number):
"""
Create a new Model instance.
Parameters
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
"""
self.param = param
self.run_number = run_number
# Create SimPy environment
self.env = simpy.Environment()
# Create resource
self.doctor = MonitoredResource(
self.env, capacity=self.param.number_of_doctors
)
# Create a random seed sequence based on the run number
ss = np.random.SeedSequence(self.run_number)
seeds = ss.spawn(2)
# Set up attributes to store results
self.patients = []
self.results_list = []
# Initialise distributions
self.arrival_dist = Exponential(mean=self.param.interarrival_time,
random_seed=seeds[0])
self.consult_dist = Exponential(mean=self.param.consultation_time,
random_seed=seeds[1])
def generate_arrivals(self):
"""
Process that generates patient arrivals.
"""
while True:
# Sample and pass time to next arrival
sampled_iat = self.arrival_dist.sample()
yield self.env.timeout(sampled_iat)
# Check whether arrived during warm-up or data collection
if self.env.now < self.param.warm_up_period:
period = "\U0001F538 WU"
else:
period = "\U0001F539 DC"
# Create a new patient
patient = Patient(patient_id=len(self.patients)+1,
period=period,
arrival_time=self.env.now)
self.patients.append(patient)
# Print arrival time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"arrives at time: {patient.arrival_time:.3f}")
# Start process of consultation
self.env.process(self.consultation(patient))
def consultation(self, patient):
"""
Process that simulates a consultation.
Parameters
----------
patient :
Instance of the Patient() class representing a single patient.
"""
# Patient requests access to a doctor (resource)
with self.doctor.request() as req:
yield req
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"starts consultation at: {self.env.now:.3f}")
# Sample consultation duration and pass time spent with doctor
time_with_doctor = self.consult_dist.sample()
yield self.env.timeout(time_with_doctor)
# Record end time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"leaves at: {self.env.now:.3f}")
def reset_results(self):
"""
Reset results.
"""
self.patients = []
self.doctor.init_results()
def warmup(self):
"""
Reset result collection after the warm-up period.
"""
if self.param.warm_up_period > 0:
# Delay process until warm-up period has completed
yield self.env.timeout(self.param.warm_up_period)
# Reset results variables
self.reset_results()
if self.param.verbose:
print(f"Warm up period ended at time: {self.env.now}")
def run(self):
"""
Run the simulation.
"""
# Schedule arrival generator and warm-up
self.env.process(self.generate_arrivals())
self.env.process(self.warmup())
# Run the simulation
self.env.run(until=(self.param.warm_up_period +
self.param.data_collection_period))
# At simulation end, update time-weighted statistics by accounting
# for the time from the last event up to the simulation finish.
self.doctor.update_time_weighted_stats()
# Create list of dictionaries containing each patient's attributes
self.results_list = [x.__dict__ for x in self.patients]Runner class
Mean utilisation can now be calculated as:
\[ \text{mean utilisation} = \frac{\text{area\_under\_the\_curve}}{\text{number of doctors} \times \text{data collection period}} \]
class Runner:
"""
Run the simulation.
Attributes
----------
param : Parameters
Simulation parameters.
"""
def __init__(self, param):
"""
Initialise a new instance of the Runner class.
Parameters
----------
param : Parameters
Simulation parameters.
"""
self.param = param
def run_single(self, run):
"""
Runs the simulation once and performs results calculations.
Parameters
----------
run : int
Run number for the simulation.
Returns
-------
dict
Contains patient-level results and results from each run.
"""
model = Model(param=self.param, run_number=run)
model.run()
# Patient results
patient_results = pd.DataFrame(model.results_list)
patient_results["run"] = run
# Run results
run_results = {
"run_number": run,
"mean_utilisation_tw": (
sum(model.doctor.area_resource_busy) / (
self.param.number_of_doctors *
self.param.data_collection_period
)
)
}
return {
"patient": patient_results,
"run": run_results
}Run the model
patient_id period arrival_time run
0 1 🔹 DC 37.778250 0
1 2 🔹 DC 38.108184 0
2 3 🔹 DC 42.610666 0
3 4 🔹 DC 44.087985 0
4 5 🔹 DC 55.587763 0
5 6 🔹 DC 64.880007 0
6 7 🔹 DC 64.954293 0{'run_number': 0, 'mean_utilisation_tw': 0.8743911437929205}Mean utilisation measures the proportion of available doctor time that is actually used during the data collection period.
None▶ View how to record utilisation
Utilisation function
Our function for utilisation is a time-weighted calculation: it computes a weighted mean, multiplying each interval’s utilisation by its interval duration.
As described in the docstring, the calculation is adapted from the plot.resources.utilization() function in simmer.plot 0.1.18, which is shared under an MIT Licence (Ucar I, Smeets B (2023). simmer.plot: Plotting Methods for ‘simmer’. https://r-simmer.org https://github.com/r-simmer/simmer.plot).
The function is explained line-by-line below…
#' Calculate the resource utilisation
#'
#' Utilisation is given by the total effective usage time (`in_use`) over
#' the total time intervals considered (`dt`).
#'
#' Credit: The utilisation calculation is adapted from the
#' `plot.resources.utilization()` function in simmer.plot 0.1.18, which is
#' shared under an MIT Licence (Ucar I, Smeets B (2023). simmer.plot: Plotting
#' Methods for 'simmer'. https://r-simmer.org
#' https://github.com/r-simmer/simmer.plot.).
#'
#' @param resources Dataframe with times patients use or queue for resources.
#' @param simulation_end_time Time at end of simulation run.
#' @param summarise If TRUE, return overall utilisation. If FALSE, just return
#' the resource dataframe with the additional columns interval_duration,
#' effective_capacity and utilisation.
#'
#' @importFrom dplyr group_by mutate n row_number summarise ungroup
#' @importFrom rlang .data
#' @importFrom tidyr pivot_wider
#'
#' @return Tibble with columns containing result for each resource.
#' @export
calc_utilisation <- function(
resources, simulation_end_time, summarise = TRUE
) {
# Calculate utilisation
util_df <- resources |>
group_by(.data[["resource"]]) |>
mutate(
# Calculate time between this row and the next. For final row,
# lead(time) returns NA, so use simulation_end_time - time
interval_duration = ifelse(
row_number() == n(),
simulation_end_time - .data[["time"]],
lead(.data[["time"]]) - .data[["time"]]
),
# Ensures effective capacity is never less than number of servers in
# use (in case of situations where servers may exceed "capacity").
effective_capacity = pmax(.data[["capacity"]], .data[["server"]]),
# Divide number of servers in use by effective capacity
# Set to NA if effective capacity is 0
utilisation = ifelse(.data[["effective_capacity"]] > 0L,
.data[["server"]] / .data[["effective_capacity"]],
NA_real_)
)
# If summarise = TRUE, find total utilisation
if (summarise) {
util_df |>
summarise(
# Multiply each utilisation by its own unique duration. The total of
# those is then divided by the total duration of all intervals.
# Hence, we are calculated a time-weighted average utilisation.
utilisation = (sum(.data[["utilisation"]] *
.data[["interval_duration"]], na.rm = TRUE) /
sum(.data[["interval_duration"]], na.rm = TRUE))
) |>
pivot_wider(names_from = "resource",
values_from = "utilisation",
names_glue = "utilisation_{resource}") |>
ungroup()
} else {
# If summarise = FALSE, just return the util_df with no further processing
ungroup(util_df)
}
}Explaining the utilisation function
#' Calculate the resource utilisation
#'
#' Utilisation is given by the total effective usage time (`in_use`) over
#' the total time intervals considered (`dt`).
#'
#' Credit: The utilisation calculation is adapted from the
#' `plot.resources.utilization()` function in simmer.plot 0.1.18, which is
#' shared under an MIT Licence (Ucar I, Smeets B (2023). simmer.plot: Plotting
#' Methods for 'simmer'. https://r-simmer.org
#' https://github.com/r-simmer/simmer.plot.).
#'
#' @param resources Dataframe with times patients use or queue for resources.
#' @param simulation_end_time Time at end of simulation run.
#' @param summarise If TRUE, return overall utilisation. If FALSE, just return
#' the resource dataframe with the additional columns interval_duration,
#' effective_capacity and utilisation.
#'
#' @importFrom dplyr group_by mutate ungroup
#' @importFrom rlang .data
#' @importFrom tidyr pivot_wider
#'
#' @return Tibble with columns containing result for each resource.
#' @export
calc_utilisation <- function(
resources, simulation_end_time, summarise = TRUE
) {The function takes the resources dataframe from get_mon_resources().
The simulation_end_time parameter will be used in the time-weighted calculations.
The summarise parameter controls whether to return the overall time-weighted utilisation (TRUE), or a detailed dataframe with per-interval metrics (FALSE).
The code groups by resource, so we can calculate utilisation of each resource separately. The mutate() call will then add new columns to the dataframe…
The interval_duration column computes how long each state lasts by taking the time difference between consecutive events. This is the basis for time-weighting.
The effective_capacity ensures the denominator in utilisation calculations is always at least as large as the number of servers in use.
The utilisation column calculates, for each interval, the proportion of servers in use relative to available capacity. If capacity is zero, it returns NA.
# If summarise = TRUE, find total utilisation
if (summarise) {
util_df |>
summarise(
# Multiply each utilisation by its own unique duration. The total of
# those is then divided by the total duration of all intervals.
# Hence, we are calculated a time-weighted average utilisation.
utilisation = (sum(.data[["utilisation"]] *
.data[["interval_duration"]], na.rm = TRUE) /
sum(.data[["interval_duration"]], na.rm = TRUE))
) |>
pivot_wider(names_from = "resource",
values_from = "utilisation",
names_glue = "utilisation_{resource}") |>
ungroup()If summarise is true, the final result is a time-weighted average: per-interval utilisations are weighted by their durations and averaged. Output is one value per resource.
When summarise is false, detailed interval-level results are returned for further analysis or plotting.
Run results function
In get_run_results(), the table results[["resources"]] is passed to calc_utilisation().
We also add a new parameter simulation_end_time.
#' Get average results for the provided single run.
#'
#' @param results Named list with `arrivals` from `get_mon_arrivals()` and
#' `resources` from `get_mon_resources()`
#' (`per_resource = TRUE` and `ongoing = TRUE`).
#' @param run_number Integer index of current simulation run.
#' @param simulation_end_time Time at end of simulation run.
#'
#' @importFrom dplyr bind_cols
#' @importFrom tibble tibble
#'
#' @return Tibble with processed results from replication.
#' @export
get_run_results <- function(results, run_number, simulation_end_time) {
metrics <- list(
calc_utilisation(results[["resources"]], simulation_end_time)
)
dplyr::bind_cols(tibble(replication = run_number), metrics)
}Model function
We need to pass the full simulation run time (i.e. warm-up period + data collection period) to get_run_results() in model(). To avoid calculating this twice (for simmer::run() and get_run_results()), we create a new variable full_run_length and use that for both.
#' Run simulation.
#'
#' @param param List. Model parameters.
#' @param run_number Numeric. Run number for random seed generation.
#'
#' @importFrom dplyr mutate
#' @importFrom simmer add_generator add_resource get_mon_arrivals
#' @importFrom simmer get_mon_resources release run seize simmer timeout
#' @importFrom simmer trajectory
#'
#' @return Named list with tables `arrivals`, `resources` and `run_results`.
#' @export
model <- function(param, run_number) {
# Set random seed based on run number
set.seed(run_number)
# Create simmer environment
env <- simmer("simulation", verbose = param[["verbose"]])
# Calculate full run length
full_run_length <- (param[["warm_up_period"]] +
param[["data_collection_period"]])
# Define the patient trajectory
patient <- trajectory("consultation") |>
seize("doctor", 1L) |>
timeout(function() {
rexp(n = 1L, rate = 1L / param[["consultation_time"]])
}) |>
release("doctor", 1L)
env <- env |>
# Add doctor resource
add_resource("doctor", param[["number_of_doctors"]]) |>
# Add patient generator
add_generator("patient", patient, function() {
rexp(n = 1L, rate = 1L / param[["interarrival_time"]])
}) |>
# Run the simulation
simmer::run(until = full_run_length)
# Extract information on arrivals and resources from simmer environment
result <- list(
arrivals = get_mon_arrivals(env, per_resource = TRUE, ongoing = TRUE),
resources = get_mon_resources(env)
)
# Filter to remove results from the warm-up period
result <- filter_warmup(result, param[["warm_up_period"]])
# Replace "replication" value with appropriate run number
result[["arrivals"]] <- mutate(result[["arrivals"]],
replication = run_number)
result[["resources"]] <- mutate(result[["resources"]],
replication = run_number)
# Calculate the average results for that run
result[["run_results"]] <- get_run_results(
result, run_number, simulation_end_time = full_run_length
)
result
}Run the model
| replication | utilisation_doctor |
|---|---|
| 1 | 0.5631573 |
We can also view the result with summarise = FALSE:
| resource | time | server | queue | capacity | queue_size | system | limit | replication | interval_duration | effective_capacity | utilisation |
|---|---|---|---|---|---|---|---|---|---|---|---|
| doctor | 30.00000 | 3 | 0 | 3 | Inf | 3 | Inf | 1 | 2.1017355 | 3 | 1.0000000 |
| doctor | 32.10174 | 3 | 1 | 3 | Inf | 4 | Inf | 1 | 2.3218815 | 3 | 1.0000000 |
| doctor | 34.42362 | 3 | 0 | 3 | Inf | 3 | Inf | 1 | 6.2440047 | 3 | 1.0000000 |
| doctor | 40.66762 | 2 | 0 | 3 | Inf | 2 | Inf | 1 | 0.7960858 | 3 | 0.6666667 |
| doctor | 41.46371 | 1 | 0 | 3 | Inf | 1 | Inf | 1 | 3.5053411 | 3 | 0.3333333 |
| doctor | 44.96905 | 0 | 0 | 3 | Inf | 0 | Inf | 1 | 9.2523579 | 3 | 0.0000000 |
| doctor | 54.22141 | 1 | 0 | 3 | Inf | 1 | Inf | 1 | 5.1762197 | 3 | 0.3333333 |
| doctor | 59.39763 | 2 | 0 | 3 | Inf | 2 | Inf | 1 | 3.2737332 | 3 | 0.6666667 |
| doctor | 62.67136 | 3 | 0 | 3 | Inf | 3 | Inf | 1 | 0.0956016 | 3 | 1.0000000 |
| doctor | 62.76696 | 2 | 0 | 3 | Inf | 2 | Inf | 1 | 2.8467970 | 3 | 0.6666667 |
| doctor | 65.61376 | 3 | 0 | 3 | Inf | 3 | Inf | 1 | 2.9412039 | 3 | 1.0000000 |
| doctor | 68.55496 | 2 | 0 | 3 | Inf | 2 | Inf | 1 | 0.2682591 | 3 | 0.6666667 |
| doctor | 68.82322 | 3 | 0 | 3 | Inf | 3 | Inf | 1 | 1.0607262 | 3 | 1.0000000 |
| doctor | 69.88395 | 2 | 0 | 3 | Inf | 2 | Inf | 1 | 0.1160527 | 3 | 0.6666667 |
Mean queue length
Corresponding term in queueing theory: Average number in queue, \(L_q\)
None▶ View how to record mean queue length
Patient class
No changes required.
Monitored resource class
We need the new MonitoredResource class (introduced and explained under time-weighted utilisation).
For queue length, the class is similar, except that we are monitoring area_n_in_queue (len(self.queue) * time_since_last_event) instead of area_resource_busy. These changes are highlighted.
class MonitoredResource(simpy.Resource):
"""
SimPy Resource subclass that tracks queue length.
Attributes
----------
time_last_event : list
Time of the most recent resource request or release event.
area_n_in_queue : list
Cumulative time patients spent queueing for a resource, for
computing average queue length.
Notes
-----
Class adapted from Monks, Harper and Heather 2025.
"""
def __init__(self, *args, **kwargs):
"""
Initialise the MonitoredResource, resetting monitoring statistics.
Parameters
----------
*args :
Positional arguments to be passed to the parent class.
**kwargs :
Keyword arguments to be passed to the parent class.
"""
# Initialise a SimPy Resource
super().__init__(*args, **kwargs)
# Run the init_results() method
self.init_results()
def init_results(self):
"""
Reset monitoring attributes.
"""
self.time_last_event = [self._env.now]
self.area_n_in_queue = [0.0]
def request(self, *args, **kwargs):
"""
Update time-weighted statistics before requesting a resource.
Parameters
----------
*args :
Positional arguments to be passed to the parent class.
**kwargs :
Keyword arguments to be passed to the parent class.
Returns
-------
simpy.events.Event
Event for the resource request.
"""
# Update time-weighted statistics
self.update_time_weighted_stats()
# Request a resource
return super().request(*args, **kwargs)
def release(self, *args, **kwargs):
"""
Update time-weighted statistics before releasing a resource.
Parameters
----------
*args :
Positional arguments to be passed to the parent class.
**kwargs :
Keyword arguments to be passed to the parent class.
Returns
-------
simpy.events.Event
Event for the resource release.
"""
# Update time-weighted statistics
self.update_time_weighted_stats()
# Release a resource
return super().release(*args, **kwargs)
def update_time_weighted_stats(self):
"""
Update the time-weighted statistics for resource usage.
Notes
-----
- These sums can be referred to as "the area under the curve".
- They are called "time-weighted" statistics as they account for how
long certain events or states persist over time.
"""
# Calculate time since last event
time_since_last_event = self._env.now - self.time_last_event[-1]
# Add record of current time
self.time_last_event.append(self._env.now)
# Add "area under curve" of people in queue
# len(self.queue) is the number of requests queued
self.area_n_in_queue.append(len(self.queue) * time_since_last_event)Model class
Same changes as made for time-weighted utilisation.
Runner class
Mean queue length is calculated as:
\[ \text{mean queue length} = \frac{\text{area\_under\_the\_curve}}{\text{data collection period}} \]
class Runner:
"""
Run the simulation.
Attributes
----------
param : Parameters
Simulation parameters.
"""
def __init__(self, param):
"""
Initialise a new instance of the Runner class.
Parameters
----------
param : Parameters
Simulation parameters.
"""
self.param = param
def run_single(self, run):
"""
Runs the simulation once and performs results calculations.
Parameters
----------
run : int
Run number for the simulation.
Returns
-------
dict
Contains patient-level results and results from each run.
"""
model = Model(param=self.param, run_number=run)
model.run()
# Patient results
patient_results = pd.DataFrame(model.results_list)
patient_results["run"] = run
# Run results
run_results = {
"run_number": run,
"mean_queue_length": (
sum(model.doctor.area_n_in_queue) /
self.param.data_collection_period
)
}
return {
"patient": patient_results,
"run": run_results
}Run the model
patient_id period arrival_time run
0 1 🔹 DC 37.778250 0
1 2 🔹 DC 38.108184 0
2 3 🔹 DC 42.610666 0
3 4 🔹 DC 44.087985 0
4 5 🔹 DC 55.587763 0
5 6 🔹 DC 64.880007 0
6 7 🔹 DC 64.954293 0{'run_number': 0, 'mean_queue_length': 0.6102042855882626}Queue length function
A new function calc_mean_queue_length() finds the time-weighted mean queue length for each resource. It does this by:
- Sorting the arrivals by start time.
- Calculating the duration until the next arrival for each patient.
- Multiplying each observed queue length by its corresponding interval duration (“time weighting”).
The sum of these products is then divided by the total observed time to find the average queue length. This approach ensures that, for example, intervals where the queue was longer and persisted for more time have a greater effect on the mean than short-lived fluctuations.
#' Calculate the time-weighted mean queue length.
#'
#' @param arrivals Dataframe with times for each patient with each resource.
#' @param simulation_end_time Time at end of simulation run.
#'
#' @importFrom dplyr arrange group_by lead mutate n row_number summarise
#' @importFrom dplyr ungroup
#' @importFrom tidyr pivot_wider
#'
#' @return Tibble with column containing mean queue length.
#' @export
calc_mean_queue_length <- function(arrivals, simulation_end_time) {
arrivals |>
group_by(.data[["resource"]]) |>
# Sort by arrival time
arrange(.data[["start_time"]]) |>
# Calculate time between this row and the next. For final row,
# lead(start_time) returns NA, so use simulation_end_time - start_time
mutate(
interval_duration = ifelse(
row_number() == n(),
simulation_end_time - .data[["start_time"]],
lead(.data[["start_time"]]) - .data[["start_time"]]
)
) |>
# Multiply each queue length by its own unique duration. The total of
# those is then divided by the total duration of all intervals.
# Hence, we are calculated a time-weighted average queue length.
summarise(mean_queue_length = (
sum(.data[["queue_on_arrival"]] *
.data[["interval_duration"]], na.rm = TRUE) /
sum(.data[["interval_duration"]], na.rm = TRUE)
)
) |>
# Reshape dataframe
pivot_wider(names_from = "resource",
values_from = "mean_queue_length",
names_glue = "mean_queue_length_{resource}") |>
ungroup()
}Run results function
In get_run_results(), the table results[["arrivals"]] is passed to calc_mean_queue_length(). A new parameter simulation_end_time is add to get_run_results() (as likewise done for utilisation).
#' Get average results for the provided single run.
#'
#' @param results Named list with `arrivals` from `get_mon_arrivals()` and
#' `resources` from `get_mon_resources()`
#' (`per_resource = TRUE` and `ongoing = TRUE`).
#' @param run_number Integer index of current simulation run.
#' @param simulation_end_time Time at end of simulation run.
#'
#' @importFrom dplyr bind_cols
#' @importFrom tibble tibble
#'
#' @return Tibble with processed results from replication.
#' @export
get_run_results <- function(results, run_number, simulation_end_time) {
metrics <- list(
calc_mean_queue_length(results[["arrivals"]], simulation_end_time)
)
dplyr::bind_cols(tibble(replication = run_number), metrics)
}Model function
We need to record the length of the queue when each patient arrives (doctor_queue_on_arrival).
We can achieve this using set_attribute() and get_queue_count().
However, add this attribute to our patient results, we need to make the same changes described above for mean wait time (and mean time in consultation): setting mon = 2L in the patient generator, retrieving results using get_mon_attributes(), and restructuring the returned dataframe to transform doctor_serve_length into serve_length (aligned with the corresponding doctor rows).
We also need to pass the full run time to get_run_results(), so calculate full_run_length (as likewise done for utilisation).
#' Run simulation.
#'
#' @param param List. Model parameters.
#' @param run_number Numeric. Run number for random seed generation.
#'
#' @importFrom dplyr left_join mutate
#' @importFrom rlang .data
#' @importFrom simmer add_generator add_resource get_mon_arrivals
#' @importFrom simmer get_queue_count get_mon_resources release run seize
#' @importFrom simmer set_attribute simmer timeout trajectory
#' @importFrom tidyr pivot_wider
#'
#' @return Named list with tables `arrivals`, `resources` and `run_results`.
#' @export
model <- function(param, run_number) {
# Set random seed based on run number
set.seed(run_number)
# Create simmer environment
env <- simmer("simulation", verbose = param[["verbose"]])
# Calculate full run length
full_run_length <- (param[["warm_up_period"]] +
param[["data_collection_period"]])
# Define the patient trajectory
patient <- trajectory("consultation") |>
# Record queue length on arrival
set_attribute("doctor_queue_on_arrival",
function() get_queue_count(env, "doctor")) |>
seize("doctor", 1L) |>
timeout(function() {
rexp(n = 1L, rate = 1L / param[["consultation_time"]])
}) |>
release("doctor", 1L)
env <- env |>
# Add doctor resource
add_resource("doctor", param[["number_of_doctors"]]) |>
# Add patient generator (mon=2 so can get manual attributes)
add_generator("patient", patient, function() {
rexp(n = 1L, rate = 1L / param[["interarrival_time"]])
}, mon = 2L) |>
# Run the simulation
simmer::run(until = full_run_length)
# Extract information on arrivals and resources from simmer environment
result <- list(
arrivals = get_mon_arrivals(env, per_resource = TRUE, ongoing = TRUE),
resources = get_mon_resources(env)
)
# Get the extra arrivals attributes
extra_attributes <- get_mon_attributes(env) |>
dplyr::select("name", "key", "value") |>
# Add column with resource name, and remove that from key
mutate(resource = gsub("_.+", "", .data[["key"]]),
key = gsub("^[^_]+_", "", .data[["key"]])) |>
pivot_wider(names_from = "key", values_from = "value")
# Merge extra attributes with the arrival data
result[["arrivals"]] <- left_join(
result[["arrivals"]], extra_attributes, by = c("name", "resource")
)
# Filter to remove results from the warm-up period
result <- filter_warmup(result, param[["warm_up_period"]])
# Replace "replication" value with appropriate run number
result[["arrivals"]] <- mutate(result[["arrivals"]],
replication = run_number)
result[["resources"]] <- mutate(result[["resources"]],
replication = run_number)
# Calculate the average results for that run
result[["run_results"]] <- get_run_results(
result, run_number, simulation_end_time = full_run_length
)
result
}Run the model
| name | start_time | end_time | activity_time | resource | replication | queue_on_arrival |
|---|---|---|---|---|---|---|
| patient6 | 32.10174 | 44.96905 | 10.545432 | doctor | 1 | 0 |
| patient7 | 54.22141 | NA | NA | doctor | 1 | 0 |
| patient8 | 59.39763 | 62.76696 | 3.369335 | doctor | 1 | 0 |
| patient9 | 62.67136 | NA | NA | doctor | 1 | 0 |
| patient10 | 65.61376 | 68.55496 | 2.941204 | doctor | 1 | 0 |
| patient11 | 68.82322 | 69.88395 | 1.060726 | doctor | 1 | 0 |
| replication | mean_queue_length_doctor |
|---|---|
| 1 | 0 |
Mean time in system
Corresponding term in queueing theory: Average time in system, \(\\W\)
None▶ View how to record mean time in system
Patient class
We create a new attribute for recording the end_time of each patient.
class Patient:
"""
Represents a patient.
Attributes
----------
patient_id : int
Unique patient identifier.
period : str
Arrival period (warm up or data collection) with emoji.
arrival_time : float
Time patient entered the system (minutes).
end_time : float
Time that patient leaves (minutes), or NaN if remain in system.
"""
def __init__(self, patient_id, period, arrival_time):
"""
Initialises a new patient.
Parameters
----------
patient_id : int
Unique patient identifier.
period : str
Arrival period (warm up or data collection) with emoji.
arrival_time : float
Time patient entered the system (minutes).
"""
self.patient_id = patient_id
self.period = period
self.arrival_time = arrival_time
self.end_time = np.nanModel class
For each patient, the time at the end of their consultation is recorded as their end_time.
class Model:
"""
Simulation model.
Attributes
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
env : simpy.Environment
The SimPy environment for the simulation.
doctor : simpy.Resource
SimPy resource representing doctors.
patients : list
List of Patient objects.
results_list : list
List of dictionaries with the attributes of each patient.
arrival_dist : Exponential
Distribution used to generate random patient inter-arrival times.
consult_dist : Exponential
Distribution used to generate length of a doctor's consultation.
"""
def __init__(self, param, run_number):
"""
Create a new Model instance.
Parameters
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
"""
self.param = param
self.run_number = run_number
# Create SimPy environment
self.env = simpy.Environment()
# Create resource
self.doctor = simpy.Resource(
self.env, capacity=self.param.number_of_doctors
)
# Create a random seed sequence based on the run number
ss = np.random.SeedSequence(self.run_number)
seeds = ss.spawn(2)
# Set up attributes to store results
self.patients = []
self.results_list = []
# Initialise distributions
self.arrival_dist = Exponential(mean=self.param.interarrival_time,
random_seed=seeds[0])
self.consult_dist = Exponential(mean=self.param.consultation_time,
random_seed=seeds[1])
def generate_arrivals(self):
"""
Process that generates patient arrivals.
"""
while True:
# Sample and pass time to next arrival
sampled_iat = self.arrival_dist.sample()
yield self.env.timeout(sampled_iat)
# Check whether arrived during warm-up or data collection
if self.env.now < self.param.warm_up_period:
period = "\U0001F538 WU"
else:
period = "\U0001F539 DC"
# Create a new patient
patient = Patient(patient_id=len(self.patients)+1,
period=period,
arrival_time=self.env.now)
self.patients.append(patient)
# Print the arrival time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"arrives at time: {patient.arrival_time:.3f}")
# Start process of consultation
self.env.process(self.consultation(patient))
def consultation(self, patient):
"""
Process that simulates a consultation.
Parameters
----------
patient :
Instance of the Patient() class representing a single patient.
"""
# Patient requests access to a doctor (resource)
with self.doctor.request() as req:
yield req
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"starts consultation at: {self.env.now:.3f}")
# Sample consultation duration and pass time spent with doctor
time_with_doctor = self.consult_dist.sample()
yield self.env.timeout(time_with_doctor)
# Record end time
patient.end_time = self.env.now
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"leaves at: {patient.end_time:.3f}")
def reset_results(self):
"""
Reset results.
"""
self.patients = []
def warmup(self):
"""
Reset result collection after the warm-up period.
"""
if self.param.warm_up_period > 0:
# Delay process until warm-up period has completed
yield self.env.timeout(self.param.warm_up_period)
# Reset results variables
self.reset_results()
if self.param.verbose:
print(f"Warm up period ended at time: {self.env.now}")
def run(self):
"""
Run the simulation.
"""
# Schedule arrival generator and warm-up
self.env.process(self.generate_arrivals())
self.env.process(self.warmup())
# Run the simulation
self.env.run(until=(self.param.warm_up_period +
self.param.data_collection_period))
# Create list of dictionaries containing each patient's attributes
self.results_list = [x.__dict__ for x in self.patients]Runner class
The time each patient spent in the system is calculated from patient_results["end_time"] - patient_results["arrival_time"]. Backlogged patients will not have an end time so will automatically set to NaN.
class Runner:
"""
Run the simulation.
Attributes
----------
param : Parameters
Simulation parameters.
"""
def __init__(self, param):
"""
Initialise a new instance of the Runner class.
Parameters
----------
param : Parameters
Simulation parameters.
"""
self.param = param
def run_single(self, run):
"""
Runs the simulation once and performs results calculations.
Parameters
----------
run : int
Run number for the simulation.
Returns
-------
dict
Contains patient-level results and results from each run.
"""
model = Model(param=self.param, run_number=run)
model.run()
# Patient results
patient_results = pd.DataFrame(model.results_list)
patient_results["run"] = run
patient_results["time_in_system"] = (
patient_results["end_time"] - patient_results["arrival_time"]
)
# Run results
run_results = {
"run_number": run,
"mean_time_in_system": patient_results["time_in_system"].mean()
}
return {
"patient": patient_results,
"run": run_results
}Run the model
patient_id period arrival_time end_time run time_in_system
0 1 🔹 DC 37.778250 57.388567 0 19.610317
1 2 🔹 DC 38.108184 47.597921 0 9.489737
2 3 🔹 DC 42.610666 NaN 0 NaN
3 4 🔹 DC 44.087985 62.536870 0 18.448885
4 5 🔹 DC 55.587763 NaN 0 NaN
5 6 🔹 DC 64.880007 NaN 0 NaN
6 7 🔹 DC 64.954293 NaN 0 NaN{'run_number': 0, 'mean_time_in_system': np.float64(15.849646506571181)}Time in system function
The calc_mean_time_in_system() function finds the mean of a new column time_in_system.
#' Calculate the mean time in the system for finished patients.
#'
#' @param arrivals Dataframe with times for each patient with each resource.
#'
#' @importFrom dplyr summarise ungroup
#'
#' @return Tibble with column containing the mean time in the system.
#' @export
calc_mean_time_in_system <- function(arrivals) {
arrivals |>
summarise(
mean_time_in_system = mean(.data[["time_in_system"]], na.rm = TRUE)
) |>
ungroup()
}Run results function
In get_run_results(), the table results[["arrivals"]] is passed to calc_mean_time_in_system().
#' Get average results for the provided single run.
#'
#' @param results Named list with `arrivals` from `get_mon_arrivals()` and
#' `resources` from `get_mon_resources()`
#' (`per_resource = TRUE` and `ongoing = TRUE`).
#' @param run_number Integer index of current simulation run.
#'
#' @importFrom dplyr bind_cols
#' @importFrom tibble tibble
#'
#' @return Tibble with processed results from replication.
#' @export
get_run_results <- function(results, run_number) {
metrics <- list(
calc_mean_time_in_system(results[["arrivals"]])
)
dplyr::bind_cols(tibble(replication = run_number), metrics)
}Model function
In model() we create the new time_in_system column by finding the difference between each patient’s start and end time. This metric excludes unfinished patients - their time_in_system will set to NaN by default as they will have no recorded end time.
#' Run simulation.
#'
#' @param param List. Model parameters.
#' @param run_number Numeric. Run number for random seed generation.
#'
#' @importFrom dplyr mutate
#' @importFrom simmer add_generator add_resource get_mon_arrivals
#' @importFrom simmer get_mon_resources release run seize simmer timeout
#' @importFrom simmer trajectory
#'
#' @return Named list with tables `arrivals`, `resources` and `run_results`.
#' @export
model <- function(param, run_number) {
# Set random seed based on run number
set.seed(run_number)
# Create simmer environment
env <- simmer("simulation", verbose = param[["verbose"]])
# Define the patient trajectory
patient <- trajectory("consultation") |>
seize("doctor", 1L) |>
timeout(function() {
rexp(n = 1L, rate = 1L / param[["consultation_time"]])
}) |>
release("doctor", 1L)
env <- env |>
# Add doctor resource
add_resource("doctor", param[["number_of_doctors"]]) |>
# Add patient generator
add_generator("patient", patient, function() {
rexp(n = 1L, rate = 1L / param[["interarrival_time"]])
}) |>
# Run the simulation
simmer::run(until = (param[["warm_up_period"]] +
param[["data_collection_period"]]))
# Extract information on arrivals and resources from simmer environment
result <- list(
arrivals = get_mon_arrivals(env, per_resource = TRUE, ongoing = TRUE),
resources = get_mon_resources(env)
)
# Add time in system (unfinished patients will set to NaN)
result[["arrivals"]][["time_in_system"]] <- (
result[["arrivals"]][["end_time"]] - result[["arrivals"]][["start_time"]]
)
# Filter to remove results from the warm-up period
result <- filter_warmup(result, param[["warm_up_period"]])
# Replace "replication" value with appropriate run number
result[["arrivals"]] <- mutate(result[["arrivals"]],
replication = run_number)
result[["resources"]] <- mutate(result[["resources"]],
replication = run_number)
# Calculate the average results for that run
result[["run_results"]] <- get_run_results(result, run_number)
result
}Run the model
| name | start_time | end_time | activity_time | resource | replication | time_in_system |
|---|---|---|---|---|---|---|
| patient6 | 32.10174 | 44.96905 | 10.545432 | doctor | 1 | 12.867313 |
| patient7 | 54.22141 | NA | NA | doctor | 1 | NA |
| patient8 | 59.39763 | 62.76696 | 3.369335 | doctor | 1 | 3.369335 |
| patient9 | 62.67136 | NA | NA | doctor | 1 | NA |
| patient10 | 65.61376 | 68.55496 | 2.941204 | doctor | 1 | 2.941204 |
| patient11 | 68.82322 | 69.88395 | 1.060726 | doctor | 1 | 1.060726 |
| replication | mean_time_in_system |
|---|---|
| 1 | 5.059645 |
Mean number of patients in the system
Corresponding term in queueing theory: Mean system size or average number in the system, \(L\)
None▶ View how to record mean number of patients in the system
This is a time-weighted calculation, but unlike utilisation and queue length, we record this within Model itself (and not with an additional MonitoredResource class).
Patient class
No changes required.
Model class
We introduce three new attributes:
area_n_in_system- a list storing the area under the curve for the number of patients in the system over time.time_last_n_in_system- records last time the number-in-system statistic was updated.n_in_system- maintains current count of patients (whether waiting or being served)
When a patient arrives or leaves, we call a new update_n_in_system(inc) method. This:
- Calculates the time since the last update, then appends
n_in_system x durationtoarea_n_in_system. - Then updates
time_last_n_insystemandn_in_systemaccordingly.
Warmup logic resets only the necessary area and timing components, preserving the correct patient count for accurate post-warmup statistics.
At the simulation end, a final update with inc=0 ensures statistics reflect the last interval up to the end of the run.
class Model:
"""
Simulation model.
Attributes
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
env : simpy.Environment
The SimPy environment for the simulation.
doctor : simpy.Resource
SimPy resource representing doctors.
patients : list
List of Patient objects.
results_list : list
List of dictionaries with the attributes of each patient.
area_n_in_system : list of float
List containing incremental area contributions used for
time-weighted statistics of the number of patients in the system.
time_last_n_in_system : float
Simulation time at last update of the number-in-system statistic.
n_in_system: int
Current number of patients present in the system, including
waiting and being served.
arrival_dist : Exponential
Distribution used to generate random patient inter-arrival times.
consult_dist : Exponential
Distribution used to generate length of a doctor's consultation.
"""
def __init__(self, param, run_number):
"""
Create a new Model instance.
Parameters
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
"""
self.param = param
self.run_number = run_number
# Create SimPy environment
self.env = simpy.Environment()
# Create resource
self.doctor = simpy.Resource(
self.env, capacity=self.param.number_of_doctors
)
# Create a random seed sequence based on the run number
ss = np.random.SeedSequence(self.run_number)
seeds = ss.spawn(2)
# Set up attributes to store results
self.patients = []
self.results_list = []
self.area_n_in_system = [0]
self.time_last_n_in_system = self.env.now
self.n_in_system = 0
# Initialise distributions
self.arrival_dist = Exponential(mean=self.param.interarrival_time,
random_seed=seeds[0])
self.consult_dist = Exponential(mean=self.param.consultation_time,
random_seed=seeds[1])
def update_n_in_system(self, inc):
"""
Update the time-weighted statistics for number of patients in system.
Parameters
----------
inc : int
Change in the number of patients (+1, 0, -1).
"""
# Compute time since last event and calculate area under curve for that
duration = self.env.now - self.time_last_n_in_system
self.area_n_in_system.append(self.n_in_system * duration)
# Update time and n in system
self.time_last_n_in_system = self.env.now
self.n_in_system += inc
def generate_arrivals(self):
"""
Process that generates patient arrivals.
"""
while True:
# Sample and pass time to next arrival
sampled_iat = self.arrival_dist.sample()
yield self.env.timeout(sampled_iat)
# Check whether arrived during warm-up or data collection
if self.env.now < self.param.warm_up_period:
period = "\U0001F538 WU"
else:
period = "\U0001F539 DC"
# Create a new patient
patient = Patient(patient_id=len(self.patients)+1,
period=period,
arrival_time=self.env.now)
self.patients.append(patient)
# Update the number in the system
self.update_n_in_system(inc=1)
# Print arrival time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"arrives at time: {patient.arrival_time:.3f}")
# Start process of consultation
self.env.process(self.consultation(patient))
def consultation(self, patient):
"""
Process that simulates a consultation.
Parameters
----------
patient :
Instance of the Patient() class representing a single patient.
"""
# Patient requests access to a doctor (resource)
with self.doctor.request() as req:
yield req
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"starts consultation at: {self.env.now:.3f}")
# Sample consultation duration and pass time spent with doctor
time_with_doctor = self.consult_dist.sample()
yield self.env.timeout(time_with_doctor)
# Update number in system
self.update_n_in_system(inc=-1)
# Record end time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"leaves at: {self.env.now:.3f}")
def reset_results(self):
"""
Reset results.
"""
self.patients = []
# For number in system, we reset area and time but not the count, as
# it should include any remaining warm-up patients in the count
self.area_n_in_system = [0]
self.time_last_n_in_system = self.env.now
def warmup(self):
"""
Reset result collection after the warm-up period.
"""
if self.param.warm_up_period > 0:
# Delay process until warm-up period has completed
yield self.env.timeout(self.param.warm_up_period)
# Reset results variables
self.reset_results()
if self.param.verbose:
print(f"Warm up period ended at time: {self.env.now}")
def run(self):
"""
Run the simulation.
"""
# Schedule arrival generator and warm-up
self.env.process(self.generate_arrivals())
self.env.process(self.warmup())
# Run the simulation
self.env.run(until=(self.param.warm_up_period +
self.param.data_collection_period))
# Run final calculation of number in system
self.update_n_in_system(inc=0)
# Create list of dictionaries containing each patient's attributes
self.results_list = [x.__dict__ for x in self.patients]Runner class
Mean number of patients in the system is calculated as:
\[ \text{mean patients in system} = \frac{\text{area\_under\_the\_curve}}{\text{data collection period}} \]
class Runner:
"""
Run the simulation.
Attributes
----------
param : Parameters
Simulation parameters.
"""
def __init__(self, param):
"""
Initialise a new instance of the Runner class.
Parameters
----------
param : Parameters
Simulation parameters.
"""
self.param = param
def run_single(self, run):
"""
Runs the simulation once and performs results calculations.
Parameters
----------
run : int
Run number for the simulation.
Returns
-------
dict
Contains patient-level results and results from each run.
"""
model = Model(param=self.param, run_number=run)
model.run()
# Patient results
patient_results = pd.DataFrame(model.results_list)
patient_results["run"] = run
# Run results
run_results = {
"run_number": run,
"mean_patients_in_system": (
sum(model.area_n_in_system) /
self.param.data_collection_period
)
}
return {
"patient": patient_results,
"run": run_results
}Run the model
patient_id period arrival_time run
0 1 🔹 DC 37.778250 0
1 2 🔹 DC 38.108184 0
2 3 🔹 DC 42.610666 0
3 4 🔹 DC 44.087985 0
4 5 🔹 DC 55.587763 0
5 6 🔹 DC 64.880007 0
6 7 🔹 DC 64.954293 0{'run_number': 0, 'mean_patients_in_system': 3.2333777169670244}Patients in system function
The function calc_mean_patients_in_system() does not use the arrivals or resources dataframes - instead, it uses a new dataframe which contains patient counts over time.
It is a time-weighted calculation which multiplies the patient count in each time interval by the interval duration. The sum of these is then divided by the total simulation time to find the average system occupancy.
#' Calculate the time-weighted mean number of patients in the system.
#'
#' @param patient_count Dataframe with patient counts over time.
#' @param simulation_end_time Time at end of simulation run.
#'
#' @importFrom dplyr arrange lead mutate n row_number summarise ungroup
#'
#' @return Tibble with column containing mean number of patients in the system.
#' @export
calc_mean_patients_in_system <- function(patient_count, simulation_end_time) {
patient_count |>
# Sort by time
arrange(.data[["time"]]) |>
# Calculate time between this row and the next. For final row,
# lead(time) returns NA, so use simulation_end_time - time
mutate(
interval_duration = ifelse(
row_number() == n(),
simulation_end_time - .data[["time"]],
lead(.data[["time"]]) - .data[["time"]]
)
) |>
# Multiply each patient count by its own unique duration. The total of
# those is then divided by the total duration of all intervals.
# Hence, we are calculated a time-weighted average patient count.
summarise(
mean_patients_in_system = (
sum(.data[["count"]] * .data[["interval_duration"]], na.rm = TRUE) /
sum(.data[["interval_duration"]], na.rm = TRUE)
)
) |>
ungroup()
}Run results function
In get_run_results(), the table results[["patients_in_system"]] is passed to calc_mean_patients_in_system(). We also add a simulation_end_time parameter (as likewise done for utilisation and queue length).
#' Get average results for the provided single run.
#'
#' @param results Named list with `arrivals` from `get_mon_arrivals()` and
#' `resources` from `get_mon_resources()`
#' (`per_resource = TRUE` and `ongoing = TRUE`).
#' @param run_number Integer index of current simulation run.
#' @param simulation_end_time Time at end of simulation run.
#'
#' @importFrom dplyr bind_cols
#' @importFrom tibble tibble
#'
#' @return Tibble with processed results from replication.
#' @export
get_run_results <- function(results, run_number, simulation_end_time) {
metrics <- list(
calc_mean_patients_in_system(results[["patients_in_system"]],
simulation_end_time)
)
dplyr::bind_cols(tibble(replication = run_number), metrics)
}Model function
In model(), we create a new dataframe that gives the count of patients in the system at every time there was an arrival or exit.
To do this, we first make two dataframes: one recording arrival times (arrivals_start) and another recording exit times (arrivals_end). Each event is tagged as either a +1 or -1 change.
These event records are then combined into a single dataframe and sorted, allowing calculation of the cumulative sum of the +1 and -1 values to give us the patient count over time.
We also add a column with the run_number (so - if running multiple replications later - the results can be linked to the relevant run).
We also need to pass the full run time to get_run_results(), so calculate full_run_length (as likewise done for utilisation and queue length).
#' Run simulation.
#'
#' @param param List. Model parameters.
#' @param run_number Numeric. Run number for random seed generation.
#'
#' @importFrom dplyr arrange bind_rows desc mutate transmute
#' @importFrom rlang .data
#' @importFrom simmer add_generator add_resource get_mon_arrivals
#' @importFrom simmer get_mon_resources release run seize simmer timeout
#' @importFrom simmer trajectory
#' @importFrom tidyr drop_na
#' @importFrom tidyselect all_of
#'
#' @return Named list with tables `arrivals`, `resources` and `run_results`.
#' @export
model <- function(param, run_number) {
# Set random seed based on run number
set.seed(run_number)
# Create simmer environment
env <- simmer("simulation", verbose = param[["verbose"]])
# Calculate full run length
full_run_length <- (param[["warm_up_period"]] +
param[["data_collection_period"]])
# Define the patient trajectory
patient <- trajectory("consultation") |>
seize("doctor", 1L) |>
timeout(function() {
rexp(n = 1L, rate = 1L / param[["consultation_time"]])
}) |>
release("doctor", 1L)
env <- env |>
# Add doctor resource
add_resource("doctor", param[["number_of_doctors"]]) |>
# Add patient generator
add_generator("patient", patient, function() {
rexp(n = 1L, rate = 1L / param[["interarrival_time"]])
}) |>
# Run the simulation
simmer::run(until = full_run_length)
# Extract information on arrivals and resources from simmer environment
result <- list(
arrivals = get_mon_arrivals(env, per_resource = TRUE, ongoing = TRUE),
resources = get_mon_resources(env)
)
# Filter to remove results from the warm-up period
result <- filter_warmup(result, param[["warm_up_period"]])
# Gather all start and end times, with a row for each, marked +1 or -1
# Drop NA for end time, as those are patients who haven't left system
# at the end of the simulation
arrivals_start <- transmute(
result[["arrivals"]], time = .data[["start_time"]], change = 1L
)
arrivals_end <- result[["arrivals"]] |>
drop_na(all_of("end_time")) |>
transmute(time = .data[["end_time"]], change = -1L)
events <- bind_rows(arrivals_start, arrivals_end)
# Determine the count of patients in the service with each entry/exit
result[["patients_in_system"]] <- events |>
# Sort events by time
arrange(.data[["time"]], desc(.data[["change"]])) |>
# Use cumulative sum to find number of patients in system at each time
mutate(count = cumsum(.data[["change"]])) |>
dplyr::select(c("time", "count"))
# Replace "replication" value with appropriate run number
result[["arrivals"]] <- mutate(result[["arrivals"]],
replication = run_number)
result[["resources"]] <- mutate(result[["resources"]],
replication = run_number)
result[["patients_in_system"]] <- mutate(result[["patients_in_system"]],
replication = run_number)
# Calculate the average results for that run
result[["run_results"]] <- get_run_results(
result, run_number, simulation_end_time = full_run_length
)
result
}Run the model
Here, we see the new result[["patient_in_system"]] dataframe, as well as the result for mean_patients_in_system in result[["run_results"]].
| time | count | replication |
|---|---|---|
| 32.10174 | 1 | 1 |
| 44.96905 | 0 | 1 |
| 54.22141 | 1 | 1 |
| 59.39763 | 2 | 1 |
| 62.67136 | 3 | 1 |
| 62.76696 | 2 | 1 |
| 65.61376 | 3 | 1 |
| 68.55496 | 2 | 1 |
| 68.82322 | 3 | 1 |
| 69.88395 | 2 | 1 |
| replication | mean_patients_in_system |
|---|---|
| 1 | 1.143741 |
Unseen (backlogged) patients
Unseen patients are those still waiting in the queue, having not started their consultation when the simulation ends.
Measuring them is important, otherwise system congestion and long wait times may be missed. Their wait time is not included with completed consultations because these waits are censored - i.e., the full duration is unknown until consultation begins. Tracking them separately avoids distorting the calculated averages, whilst accurately reflecting system backlog at simulation end.
Backlogged patient count
None▶ View how to record the backlogged patient count
Patient class
This requires the same change as for mean time in consultation - adding an attribute time_with_doctor.
class Patient:
"""
Represents a patient.
Attributes
----------
patient_id : int
Unique patient identifier.
period : str
Arrival period (warm up or data collection) with emoji.
arrival_time : float
Time patient entered the system (minutes).
time_with_doctor : float
Time spent in consultation with a doctor (minutes).
"""
def __init__(self, patient_id, period, arrival_time):
"""
Initialises a new patient.
Parameters
----------
patient_id : int
Unique patient identifier.
period : str
Arrival period (warm up or data collection) with emoji.
arrival_time : float
Time patient entered the system (minutes).
"""
self.patient_id = patient_id
self.period = period
self.arrival_time = arrival_time
self.time_with_doctor = np.nanModel class
Again, we make the same change as for mean time in consultation - recording the sampled time with the doctor as patient.time_with_doctor during Model.consultation().
class Model:
"""
Simulation model.
Attributes
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
env : simpy.Environment
The SimPy environment for the simulation.
doctor : simpy.Resource
SimPy resource representing doctors.
patients : list
List of Patient objects.
results_list : list
List of dictionaries with the attributes of each patient.
arrival_dist : Exponential
Distribution used to generate random patient inter-arrival times.
consult_dist : Exponential
Distribution used to generate length of a doctor's consultation.
"""
def __init__(self, param, run_number):
"""
Create a new Model instance.
Parameters
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
"""
self.param = param
self.run_number = run_number
# Create SimPy environment
self.env = simpy.Environment()
# Create resource
self.doctor = simpy.Resource(
self.env, capacity=self.param.number_of_doctors
)
# Create a random seed sequence based on the run number
ss = np.random.SeedSequence(self.run_number)
seeds = ss.spawn(2)
# Set up attributes to store results
self.patients = []
self.results_list = []
# Initialise distributions
self.arrival_dist = Exponential(mean=self.param.interarrival_time,
random_seed=seeds[0])
self.consult_dist = Exponential(mean=self.param.consultation_time,
random_seed=seeds[1])
def generate_arrivals(self):
"""
Process that generates patient arrivals.
"""
while True:
# Sample and pass time to next arrival
sampled_iat = self.arrival_dist.sample()
yield self.env.timeout(sampled_iat)
# Check whether arrived during warm-up or data collection
if self.env.now < self.param.warm_up_period:
period = "\U0001F538 WU"
else:
period = "\U0001F539 DC"
# Create a new patient
patient = Patient(patient_id=len(self.patients)+1,
period=period,
arrival_time=self.env.now)
self.patients.append(patient)
# Print arrival time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"arrives at time: {patient.arrival_time:.3f}")
# Start process of consultation
self.env.process(self.consultation(patient))
def consultation(self, patient):
"""
Process that simulates a consultation.
Parameters
----------
patient :
Instance of the Patient() class representing a single patient.
"""
# Patient requests access to a doctor (resource)
with self.doctor.request() as req:
yield req
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} starts" +
f" consultation at: {self.env.now:.3f}")
# Sample consultation duration and pass time spent with doctor
patient.time_with_doctor = self.consult_dist.sample()
yield self.env.timeout(patient.time_with_doctor)
# Record end time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"leaves at: {self.env.now:.3f}")
def reset_results(self):
"""
Reset results.
"""
self.patients = []
def warmup(self):
"""
Reset result collection after the warm-up period.
"""
if self.param.warm_up_period > 0:
# Delay process until warm-up period has completed
yield self.env.timeout(self.param.warm_up_period)
# Reset results variables
self.reset_results()
if self.param.verbose:
print(f"Warm up period ended at time: {self.env.now}")
def run(self):
"""
Run the simulation.
"""
# Schedule arrival generator and warm-up
self.env.process(self.generate_arrivals())
self.env.process(self.warmup())
# Run the simulation
self.env.run(until=(self.param.warm_up_period +
self.param.data_collection_period))
# Create list of dictionaries containing each patient's attributes
self.results_list = [x.__dict__ for x in self.patients]Runner class
We can count the number of unseen patients by finding how many have np.nan for time_with_doctor.
class Runner:
"""
Run the simulation.
Attributes
----------
param : Parameters
Simulation parameters.
"""
def __init__(self, param):
"""
Initialise a new instance of the Runner class.
Parameters
----------
param : Parameters
Simulation parameters.
"""
self.param = param
def run_single(self, run):
"""
Runs the simulation once and performs results calculations.
Parameters
----------
run : int
Run number for the simulation.
Returns
-------
dict
Contains patient-level results and results from each run.
"""
model = Model(param=self.param, run_number=run)
model.run()
# Patient results
patient_results = pd.DataFrame(model.results_list)
patient_results["run"] = run
# Run results
run_results = {
"run_number": run,
"unseen_count": patient_results["time_with_doctor"].isna().sum()
}
return {
"patient": patient_results,
"run": run_results
}Run the model
In result["patient"] we can see one person had NaN for time_with_doctor. Correspondingly, the unseen_count in result["run"] is 1.
patient_id period arrival_time time_with_doctor run
0 1 🔹 DC 37.778250 19.610317 0
1 2 🔹 DC 38.108184 9.489737 0
2 3 🔹 DC 42.610666 41.665475 0
3 4 🔹 DC 44.087985 5.874479 0
4 5 🔹 DC 55.587763 27.882444 0
5 6 🔹 DC 64.880007 24.914615 0
6 7 🔹 DC 64.954293 NaN 0{'run_number': 0, 'unseen_count': np.int64(1)}Backlogged patient mean wait time
None▶ View how to record the backlogged patient mean wait time
Patient class
As done in mean time in consultation and unseen patient count, we add a new attribute time_with_doctor.
class Patient:
"""
Represents a patient.
Attributes
----------
patient_id : int
Unique patient identifier.
period : str
Arrival period (warm up or data collection) with emoji.
arrival_time : float
Time patient entered the system (minutes).
time_with_doctor : float
Time spent in consultation with a doctor (minutes).
"""
def __init__(self, patient_id, period, arrival_time):
"""
Initialises a new patient.
Parameters
----------
patient_id : int
Unique patient identifier.
period : str
Arrival period (warm up or data collection) with emoji.
arrival_time : float
Time patient entered the system (minutes).
"""
self.patient_id = patient_id
self.period = period
self.arrival_time = arrival_time
self.time_with_doctor = np.nanModel class
We also repeat the change to Model - recording the sampled time with doctor as patient.time_with_doctor during Model.consultation().
class Model:
"""
Simulation model.
Attributes
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
env : simpy.Environment
The SimPy environment for the simulation.
doctor : simpy.Resource
SimPy resource representing doctors.
patients : list
List of Patient objects.
results_list : list
List of dictionaries with the attributes of each patient.
arrival_dist : Exponential
Distribution used to generate random patient inter-arrival times.
consult_dist : Exponential
Distribution used to generate length of a doctor's consultation.
"""
def __init__(self, param, run_number):
"""
Create a new Model instance.
Parameters
----------
param : Parameters
Simulation parameters.
run_number : int
Run number for random seed generation.
"""
self.param = param
self.run_number = run_number
# Create SimPy environment
self.env = simpy.Environment()
# Create resource
self.doctor = simpy.Resource(
self.env, capacity=self.param.number_of_doctors
)
# Create a random seed sequence based on the run number
ss = np.random.SeedSequence(self.run_number)
seeds = ss.spawn(2)
# Set up attributes to store results
self.patients = []
self.results_list = []
# Initialise distributions
self.arrival_dist = Exponential(mean=self.param.interarrival_time,
random_seed=seeds[0])
self.consult_dist = Exponential(mean=self.param.consultation_time,
random_seed=seeds[1])
def generate_arrivals(self):
"""
Process that generates patient arrivals.
"""
while True:
# Sample and pass time to next arrival
sampled_iat = self.arrival_dist.sample()
yield self.env.timeout(sampled_iat)
# Check whether arrived during warm-up or data collection
if self.env.now < self.param.warm_up_period:
period = "\U0001F538 WU"
else:
period = "\U0001F539 DC"
# Create a new patient
patient = Patient(patient_id=len(self.patients)+1,
period=period,
arrival_time=self.env.now)
self.patients.append(patient)
# Record and print the arrival time
patient.arrival_time = self.env.now
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"arrives at time: {patient.arrival_time:.3f}")
# Start process of consultation
self.env.process(self.consultation(patient))
def consultation(self, patient):
"""
Process that simulates a consultation.
Parameters
----------
patient :
Instance of the Patient() class representing a single patient.
"""
# Patient requests access to a doctor (resource)
with self.doctor.request() as req:
yield req
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} starts" +
f" consultation at: {self.env.now:.3f}")
# Sample consultation duration and pass time spent with doctor
patient.time_with_doctor = self.consult_dist.sample()
yield self.env.timeout(patient.time_with_doctor)
# Record end time
if self.param.verbose:
print(f"{patient.period} Patient {patient.patient_id} " +
f"leaves at: {self.env.now:.3f}")
def reset_results(self):
"""
Reset results.
"""
self.patients = []
def warmup(self):
"""
Reset result collection after the warm-up period.
"""
if self.param.warm_up_period > 0:
# Delay process until warm-up period has completed
yield self.env.timeout(self.param.warm_up_period)
# Reset results variables
self.reset_results()
if self.param.verbose:
print(f"Warm up period ended at time: {self.env.now}")
def run(self):
"""
Run the simulation.
"""
# Schedule arrival generator and warm-up
self.env.process(self.generate_arrivals())
self.env.process(self.warmup())
# Run the simulation
self.env.run(until=(self.param.warm_up_period +
self.param.data_collection_period))
# Create list of dictionaries containing each patient's attributes
self.results_list = [x.__dict__ for x in self.patients]Runner class
After the simulation finishes, we check whether patients have been seen yet by a doctor by looking for missing values in the time_with_doctor column.
For those patients, their waiting time (unseen_wait_time) is calculated as the difference between the current simulation time and their recorded arrival time. We can then find the mean of this column.
Patients who have already seen a doctor have their unseen_wait_time set to NaN, indicating that it does not apply to them.
class Runner:
"""
Run the simulation.
Attributes
----------
param : Parameters
Simulation parameters.
"""
def __init__(self, param):
"""
Initialise a new instance of the Runner class.
Parameters
----------
param : Parameters
Simulation parameters.
"""
self.param = param
def run_single(self, run):
"""
Runs the simulation once and performs results calculations.
Parameters
----------
run : int
Run number for the simulation.
Returns
-------
dict
Contains patient-level results and results from each run.
"""
model = Model(param=self.param, run_number=run)
model.run()
# Patient results
patient_results = pd.DataFrame(model.results_list)
patient_results["run"] = run
# For each patient, if they haven't seen a doctor yet, calculate
# their wait as current time minus arrival, else set as missing
patient_results["unseen_wait_time"] = np.where(
patient_results["time_with_doctor"].isna(),
model.env.now - patient_results["arrival_time"], np.nan
)
# Run results
run_results = {
"run_number": run,
"unseen_wait_time": patient_results["unseen_wait_time"].mean()
}
return {
"patient": patient_results,
"run": run_results
}Run the model
patient_id period arrival_time time_with_doctor run unseen_wait_time
0 1 🔹 DC 37.778250 19.610317 0 NaN
1 2 🔹 DC 38.108184 9.489737 0 NaN
2 3 🔹 DC 42.610666 41.665475 0 NaN
3 4 🔹 DC 44.087985 5.874479 0 NaN
4 5 🔹 DC 55.587763 27.882444 0 NaN
5 6 🔹 DC 64.880007 24.914615 0 NaN
6 7 🔹 DC 64.954293 NaN 0 5.045707{'run_number': 0, 'unseen_wait_time': np.float64(5.045706616721617)}
None▶ View how to record the backlogged patient count and mean wait time
We will record two performance measures here: the count of unseen patients and their mean wait time (ie., time between arrival and end of simulation).
Unseen patient count function
The calc_unseen_n function finds the count of patients who have a result in the new wait_time_unseen column.
#' Calculate the number of patients still waiting for resource at end
#'
#' @param arrivals Dataframe with times for each patient with each resource.
#'
#' @importFrom dplyr arrange lead mutate summarise ungroup
#' @importFrom rlang .data
#'
#' @return Tibble with columns containing result for each resource.
#' @export
calc_unseen_n <- function(arrivals, groups = NULL) {
arrivals |>
group_by(.data[["resource"]]) |>
summarise(value = sum(!is.na(.data[["wait_time_unseen"]]))) |>
pivot_wider(names_from = "resource",
values_from = "value",
names_glue = "count_unseen_{resource}") |>
ungroup()
}Unseen patient mean wait time function
The calc_unseen_mean function finds the mean of a new column wait_time_unseen.
#' Calculate the mean wait time of patients who are still waiting for a
#' resource at the end of the simulation
#'
#' @param arrivals Dataframe with times for each patient with each resource.
#'
#' @importFrom dplyr group_by summarise ungroup
#' @importFrom rlang .data
#' @importFrom tidyr pivot_wider
#'
#' @return Tibble with columns containing result for each resource.
#' @export
calc_unseen_mean <- function(arrivals) {
arrivals |>
group_by(.data[["resource"]]) |>
summarise(value = mean(.data[["wait_time_unseen"]], na.rm = TRUE)) |>
pivot_wider(names_from = "resource",
values_from = "value",
names_glue = "mean_waiting_time_unseen_{resource}") |>
ungroup()
}Run results function
In get_run_results(), the table results[["arrivals"]] is passed to calc_unseen_n() and calc_unseen_mean.
#' Get average results for the provided single run.
#'
#' @param results Named list with `arrivals` from `get_mon_arrivals()` and
#' `resources` from `get_mon_resources()`
#' (`per_resource = TRUE` and `ongoing = TRUE`).
#' @param run_number Integer index of current simulation run.
#'
#' @importFrom dplyr bind_cols
#' @importFrom tibble tibble
#'
#' @return Tibble with processed results from replication.
#' @export
get_run_results <- function(results, run_number) {
metrics <- list(
calc_unseen_n(results[["arrivals"]]),
calc_unseen_mean(results[["arrivals"]])
)
dplyr::bind_cols(tibble(replication = run_number), metrics)
}Model function
The modifications in model() are very similar to those described for calculating mean wait time of finished patients.
The only difference here is that we make a column called wait_time_unseen, which records the time elapsed between their arrival and the end of the simulation, for patients who are unseen at the end of the simulation (i.e., their serve_start is NA).
The earlier changes to add serve_start to the arrivals dataframe are the same as before though - saving with set_attribute, setting mon = 2L in the patient generator, retrieving results using get_mon_attributes(), and restructuring the returned dataframe to transform doctor_serve_length into serve_length (aligned with the corresponding doctor rows).
#' Run simulation.
#'
#' @param param List. Model parameters.
#' @param run_number Numeric. Run number for random seed generation.
#'
#' @importFrom dplyr left_join mutate
#' @importFrom rlang .data
#' @importFrom simmer add_generator add_resource get_mon_arrivals
#' @importFrom simmer get_mon_resources now release run seize set_attribute
#' @importFrom simmer simmer timeout trajectory
#' @importFrom tidyr pivot_wider
#'
#' @return Named list with tables `arrivals`, `resources` and `run_results`.
#' @export
model <- function(param, run_number) {
# Set random seed based on run number
set.seed(run_number)
# Create simmer environment
env <- simmer("simulation", verbose = param[["verbose"]])
# Define the patient trajectory
patient <- trajectory("consultation") |>
seize("doctor", 1L) |>
# Record time resource is obtained
set_attribute("doctor_serve_start", function() now(env)) |>
timeout(function() {
rexp(n = 1L, rate = 1L / param[["consultation_time"]])
}) |>
release("doctor", 1L)
env <- env |>
# Add doctor resource
add_resource("doctor", param[["number_of_doctors"]]) |>
# Add patient generator (mon=2 so can get manual attributes)
add_generator("patient", patient, function() {
rexp(n = 1L, rate = 1L / param[["interarrival_time"]])
}, mon = 2L) |>
# Run the simulation
simmer::run(until = (param[["warm_up_period"]] +
param[["data_collection_period"]]))
# Extract information on arrivals and resources from simmer environment
result <- list(
arrivals = get_mon_arrivals(env, per_resource = TRUE, ongoing = TRUE),
resources = get_mon_resources(env)
)
# Get the extra arrivals attributes
extra_attributes <- get_mon_attributes(env) |>
dplyr::select("name", "key", "value") |>
# Add column with resource name, and remove that from key
mutate(resource = gsub("_.+", "", .data[["key"]]),
key = gsub("^[^_]+_", "", .data[["key"]])) |>
pivot_wider(names_from = "key", values_from = "value")
# Merge extra attributes with the arrival data
result[["arrivals"]] <- left_join(
result[["arrivals"]], extra_attributes, by = c("name", "resource")
)
# Filter to remove results from the warm-up period
result <- filter_warmup(result, param[["warm_up_period"]])
# Replace "replication" value with appropriate run number
result[["arrivals"]] <- mutate(result[["arrivals"]],
replication = run_number)
result[["resources"]] <- mutate(result[["resources"]],
replication = run_number)
# Calculate wait time for each patients who remain unseen at the end of
# the simulation
result[["arrivals"]] <- result[["arrivals"]] |>
mutate(
wait_time_unseen = ifelse(
is.na(.data[["serve_start"]]), now(env) - .data[["start_time"]], NA
)
)
# Calculate the average results for that run
result[["run_results"]] <- get_run_results(result, run_number)
result
}Run the model
| name | start_time | end_time | activity_time | resource | replication | serve_start | wait_time_unseen |
|---|---|---|---|---|---|---|---|
| patient6 | 32.10174 | 44.96905 | 10.545432 | doctor | 1 | 34.42362 | NA |
| patient8 | 59.39763 | 62.76696 | 3.369335 | doctor | 1 | 59.39763 | NA |
| patient10 | 65.61376 | 68.55496 | 2.941204 | doctor | 1 | 65.61376 | NA |
| patient11 | 68.82322 | 69.88395 | 1.060726 | doctor | 1 | 68.82322 | NA |
| patient9 | 62.67136 | NA | NA | doctor | 1 | 62.67136 | NA |
| patient7 | 54.22141 | NA | NA | doctor | 1 | 54.22141 | NA |
| replication | count_unseen_doctor | mean_waiting_time_unseen_doctor |
|---|---|---|
| 1 | 0 | NaN |
Time-weighted averages
A time-weighted average calculates the mean value of a performance measure by weighting each observed value by the length of time it persisted, then dividing by the total time observed.
This is used by three of the performance measures above: utilisation, queue length and patients in system. But why do we use it? Let’s use queue length as an example…
Consider a queue with the following time periods and sizes:
| Interval | Queue Size | Duration |
|---|---|---|
| 0 - 14.4 | 0 | 14.4 |
| 14.4 - 15.2 | 1 | 0.8 |
| 15.2 - 16.1 | 2 | 0.9 |
| 16.1 - 17.0 | 3 | 0.9 |
Simple (unweighted) average
This method ignores duration and only averages the observed queue sizes:
\[ \text{Simple average} = \frac{0 + 1 + 2 + 3}{4} = 1.5 \]
It overestimates the typical queue size because most of the time, the queue was empty.
Time-weighted average
This method weights each value by how long it lasted:
\[ \text{Time-weighted average} = \frac{(0 \times 14.4) + (1 \times 0.8) + (2 \times 0.9) + (3 \times 0.9)}{17.0} \approx 0.312 \]
The value is much closer to the true “average” experienced in time, because the queue was empty for most of the observation period.
The same logic applies for the other performance measures. In the code, we track the product of value and duration (“area under the curve”) for each metric, then divide the total by the observation time.
It’s called “area under the curve” because, when you plot the measure (e.g., queue length) over time, the total area of all the rectangles you get under the curve (height=queue x width=time) returns your average.
Explore the example models
Nurse visit simulation
Click to visit pydesrap_mms repository
| Key files | simulation/patient.pysimulation/monitoredresource.pysimulation/model.pysimulation/runner.py |
| What to look for? | All of the metrics above can be seen in this example repository. |
| Why it matters? | Seeing all these metrics together in the example repository can feel overwhelming! So this page breaks them down individually and shows how each is captured in the code. In the nurse visit repository, there is some additional logic in places to allow the model to handle a scenario with zero arrivals. |
Click to visit rdesrap_mms repository
| Key files | R/get_run_results.RR/model.R |
| What to look for? | All of the metrics above can be seen in this example repository. |
| Why it matters? | Seeing all these metrics together in the example repository can feel overwhelming! So this page breaks them down individually and shows how each is captured in the code. In the nurse visit repository, there is some additional logic in places to allow the model to handle a scenario with zero arrivals, and also to allow calculations to be performed by group. |
Stroke pathway simulation
Click to visit pydesrap_stroke repository
| Key files | simulation/patient.pysimulation/model.pysimulation/runner.py |
| What to look for? | This model computes less standard metrics: occupancy frequency, percentages, cumulative percentages, probability of delay, and “1 in n” delays. |
| Why it matters? | While this page focuses on some common DES metrics, the stroke simulation shows examples of some more unusual metrics. Indeed, depending on your model structure, there is a huge array of metrics you could calculate, well beyond those described on this page. |
Click to visit rdesrap_stroke repository
| Key files | R/get_occupancy_stats.RR/model.R |
| What to look for? | This model computes less standard metrics: occupancy frequency, percentages, cumulative percentages, probability of delay, and “1 in n” delays. |
| Why it matters? | While this page focuses on some common DES metrics, the stroke simulation shows examples of some more unusual metrics. Indeed, depending on your model structure, there is a huge array of metrics you could calculate, well beyond those described on this page. |
Test yourself
Try adding one or more of the measures described above to your own DES model.
See how the results change as you adjust the model (e.g. arrival rate, number of resources).
Acknowledgements
The MonitoredResource class is based on Tom Monks, Alison Harper and Amy Heather (2025) An introduction to Discrete-Event Simulation (DES) using Free and Open Source Software (https://github.com/pythonhealthdatascience/intro-open-sim/tree/main) (MIT Licence). They based it on the method described in Law. Simulation Modeling and Analysis 4th Ed. Pages 14 - 17.