Performance 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.

import numpy as np
import pandas as pd
import simpy
from sim_tools.distributions import Exponential
library(dplyr)
library(simmer)
library(tidyr)


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 = verbose

Patient 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)
            self.patients.append(patient)

            # Print arrival time
            if self.param.verbose:
                print(f"{patient.period} Patient {patient.patient_id} " +
                      f"arrives at time: {self.env.now:.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)

    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.

runner = Runner(param=Parameters())
result = runner.run_single(run=0)
print(result["patient"])
   patient_id period  run
0           1   πŸ”Ή DC    0
1           2   πŸ”Ή DC    0
2           3   πŸ”Ή DC    0
3           4   πŸ”Ή DC    0
4           5   πŸ”Ή DC    0
5           6   πŸ”Ή DC    0
6           7   πŸ”Ή DC    0
print(result["run"])
{'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 simmer add_generator get_mon_arrivals get_mon_resources run
#' @importFrom simmer simmer timeout trajectory

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

param <- create_params()
result <- model(param = param, run_number = 1L)
knitr::kable(result[["arrivals"]] |> arrange(start_time))
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 
knitr::kable(result[["resources"]])
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 
knitr::kable(result[["run_results"]])
replication
1

At the moment, result[["run_results"]] is pretty empty (just has the replication number). Let’s add some performance measures!

πŸšΆβ€β™€οΈ 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

runner = Runner(param=Parameters())
result = runner.run_single(run=0)
print(result["patient"])
   patient_id period  run
0           1   πŸ”Ή DC    0
1           2   πŸ”Ή DC    0
2           3   πŸ”Ή DC    0
3           4   πŸ”Ή DC    0
4           5   πŸ”Ή DC    0
5           6   πŸ”Ή DC    0
6           7   πŸ”Ή DC    0
print(result["run"])
{'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.
#'
#' @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.

param <- create_params()
result <- model(param = param, run_number = 1L)
knitr::kable(result[["arrivals"]] |> arrange(start_time))
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 
knitr::kable(result[["run_results"]])
replication arrivals
1 6

⏱️ Mean wait time

Corresponding term in queueing theory: Average wait time in queue, \(\\W_q\)

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.
    wait_time : float
        Time spent waiting for the doctor (minutes).
    """
    def __init__(self, patient_id, period):
        """
        Initialises a new patient.

        Parameters
        ----------
        patient_id : int
            Unique patient identifier.
        period : str
            Arrival period (warm up or data collection) with emoji.
        """
        self.patient_id = patient_id
        self.period = period
        self.wait_time = np.nan

Model 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)
            self.patients.append(patient)

            # Print arrival time
            if self.param.verbose:
                print(f"{patient.period} Patient {patient.patient_id} " +
                      f"arrives at time: {self.env.now:.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)

    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"].

runner = Runner(param=Parameters())
result = runner.run_single(run=0)
print(result["patient"])
   patient_id period  wait_time  run
0           1   πŸ”Ή DC   0.000000    0
1           2   πŸ”Ή DC   0.000000    0
2           3   πŸ”Ή DC   4.987255    0
3           4   πŸ”Ή DC  12.574406    0
4           5   πŸ”Ή DC   1.800804    0
5           6   πŸ”Ή DC   0.000000    0
6           7   πŸ”Ή DC        NaN    0
print(result["run"])
{'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 tidyr drop_na
#' 
#' @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(resource) |>
    summarise(mean_wait_time = mean(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 mutate
#' @importFrom simmer add_generator get_mon_arrivals get_mon_resources run
#' @importFrom simmer set_attribute simmer timeout trajectory

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) |>
    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"]] <- result[["arrivals"]] |>
    mutate(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

param <- create_params()
result <- model(param = param, run_number = 1L)
knitr::kable(result[["arrivals"]] |> arrange(start_time))
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 
knitr::kable(result[["run_results"]])
replication mean_wait_time_doctor
1 0.3869803

🩺 Mean time in consultation

Corresponding term in queueing theory: Service time, 1/\(\mu\)

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.
    time_with_doctor : float
        Time spent in consultation with a doctor (minutes).
    """
    def __init__(self, patient_id, period):
        """
        Initialises a new patient.

        Parameters
        ----------
        patient_id : int
            Unique patient identifier.
        period : str
            Arrival period (warm up or data collection) with emoji.
        """
        self.patient_id = patient_id
        self.period = period
        self.time_with_doctor = np.nan

Model 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)
            self.patients.append(patient)

            # Print arrival time
            if self.param.verbose:
                print(f"{patient.period} Patient {patient.patient_id} " +
                      f"arrives at time: {self.env.now:.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)

    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"].

runner = Runner(param=Parameters())
result = runner.run_single(run=0)
print(result["patient"])
   patient_id period  time_with_doctor  run
0           1   πŸ”Ή DC         19.610317    0
1           2   πŸ”Ή DC          9.489737    0
2           3   πŸ”Ή DC         41.665475    0
3           4   πŸ”Ή DC          5.874479    0
4           5   πŸ”Ή DC         27.882444    0
5           6   πŸ”Ή DC         24.914615    0
6           7   πŸ”Ή DC               NaN    0
print(result["run"])
{'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.
#'
#' @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(resource) |>
    summarise(mean_serve_time = mean(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 mutate
#' @importFrom simmer add_generator get_mon_arrivals get_mon_resources run
#' @importFrom simmer set_attribute simmer timeout trajectory

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) |>
    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.

param <- create_params()
result <- model(param = param, run_number = 1L)
knitr::kable(result[["arrivals"]] |> arrange(start_time))
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 
knitr::kable(result[["run_results"]])
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.

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)
            self.patients.append(patient)

            # Print arrival time
            if self.param.verbose:
                print(f"{patient.period} Patient {patient.patient_id} " +
                      f"arrives at time: {self.env.now:.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)

    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

runner = Runner(param=Parameters())
result = runner.run_single(run=0)
print(result["patient"])
   patient_id period  run
0           1   πŸ”Ή DC    0
1           2   πŸ”Ή DC    0
2           3   πŸ”Ή DC    0
3           4   πŸ”Ή DC    0
4           5   πŸ”Ή DC    0
5           6   πŸ”Ή DC    0
6           7   πŸ”Ή DC    0
print(result["run"])
{'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.

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.


def init_results(self):
    """
    Reset monitoring attributes.
    """
    self.time_last_event = [self._env.now]
    self.area_resource_busy = [0.0]

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:

  1. The doctor attribute is now set up using MonitoredResource (instead of simpy.Resource).

  2. In reset_results(), we call the MonitoredResource method init_results() to reset their monitoring at the end of a warm-up period.

  3. In run(), once the simulation ends, the MonitoredResource method update_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)
            self.patients.append(patient)

            # Print arrival time
            if self.param.verbose:
                print(f"{patient.period} Patient {patient.patient_id} " +
                      f"arrives at time: {self.env.now:.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)

    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

runner = Runner(param=Parameters())
result = runner.run_single(run=0)
print(result["patient"])
   patient_id period  run
0           1   πŸ”Ή DC    0
1           2   πŸ”Ή DC    0
2           3   πŸ”Ή DC    0
3           4   πŸ”Ή DC    0
4           5   πŸ”Ή DC    0
5           6   πŸ”Ή DC    0
6           7   πŸ”Ή DC    0
print(result["run"])
{'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.

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 summarise If TRUE, return overall utilisation. If FALSE, just return
#' the resource dataframe with the additional columns interval_duration,
#' effective_capacity and utilisation.
#'
#' @return Tibble with columns containing result for each resource.
#' @export

calc_utilisation <- function(resources, summarise = TRUE) {

  # Calculate utilisation
  util_df <- resources |>
    group_by(resource) |>
    mutate(
      # Time between this row and the next
      interval_duration = 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 summarise If TRUE, return overall utilisation. If FALSE, just return
#' the resource dataframe with the additional columns interval_duration,
#' effective_capacity and utilisation.
#'
#' @return Tibble with columns containing result for each resource.
#' @export

calc_utilisation <- function(resources, summarise = TRUE) {

The function takes the resources dataframe from get_mon_resources().

The summarise parameter controls whether to return the overall time-weighted utilisation (TRUE), or a detailed dataframe with per-interval metrics (FALSE).


# Calculate utilisation
util_df <- resources |>
  group_by(resource) |>
  mutate(
    ...

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…


    ...
    # Time between this row and the next
    interval_duration = lead(.data[["time"]]) - .data[["time"]],
    ...

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.


    ...
    # 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"]]),
    ...

The effective_capcaity ensures the denominator in utilisation calculations is always at least as large as the number of servers in use.


    ...
    # 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_)
  )

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.


} else {
  # If summarise = FALSE, just return the util_df with no further processing
  ungroup(util_df)
}

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().

#' 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_utilisation(results[["resources"]])
  )
  dplyr::bind_cols(tibble(replication = run_number), metrics)
}

Model function

No changes required.

Run the model

param <- create_params()
result <- model(param = param, run_number = 1L)
knitr::kable(result[["run_results"]])
replication utilisation_doctor
1 0.5628561


We can also view the result with summarise = FALSE:

util_df <- calc_utilisation(resources = result[["resources"]],
                            summarise = FALSE)
knitr::kable(util_df)
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 NA 3 0.6666667

πŸ“ Mean queue length

Corresponding term in queueing theory: Average number in queue, \(L_q\)

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

runner = Runner(param=Parameters())
result = runner.run_single(run=0)
print(result["patient"])
   patient_id period  run
0           1   πŸ”Ή DC    0
1           2   πŸ”Ή DC    0
2           3   πŸ”Ή DC    0
3           4   πŸ”Ή DC    0
4           5   πŸ”Ή DC    0
5           6   πŸ”Ή DC    0
6           7   πŸ”Ή DC    0
print(result["run"])
{'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.
#'
#' @return Tibble with column containing mean queue length.
#' @export

calc_mean_queue_length <- function(arrivals) {
  arrivals |>
    group_by(resource) |>
    # Sort by arrival time
    arrange(.data[["start_time"]]) |>
    # Calculate time between this row and the next
    mutate(
      interval_duration = (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().

#' 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_queue_length(results[["arrivals"]])
  )
  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).

#' Run simulation.
#'
#' @param param List. Model parameters.
#' @param run_number Numeric. Run number for random seed generation.
#'
#' @importFrom simmer add_generator get_mon_arrivals get_mon_resources
#' @importFrom simmer get_queue_count run set_attribute simmer
#' @importFrom simmer timeout trajectory

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") |>
    # 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 = (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

param <- create_params()
result <- model(param = param, run_number = 1L)
knitr::kable(result[["arrivals"]] |> arrange(start_time))
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 
knitr::kable(result[["run_results"]])
replication mean_queue_length_doctor
1 0

⏳ Mean time in system

Corresponding term in queueing theory: Average time in system, \(\\W\)

Patient class

We create two new attributes for recording the arrival_time and 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 that patient arrives (minutes).
    end_time : float
        Time that patient leaves (minutes), or NaN if remain in system.
    """
    def __init__(self, patient_id, period):
        """
        Initialises a new patient.

        Parameters
        ----------
        patient_id : int
            Unique patient identifier.
        period : str
            Arrival period (warm up or data collection) with emoji.
        """
        self.patient_id = patient_id
        self.period = period
        self.arrival_time = np.nan
        self.end_time = np.nan

Model class

For each patient, the time when they arrive is now stored as their arrival_time, and 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)
            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} " +
                      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

    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

runner = Runner(param=Parameters())
result = runner.run_single(run=0)
print(result["patient"])
   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
print(result["run"])
{'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.
#'
#' @return Tibble with column containing the mean time in the system.

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 simmer add_generator get_mon_arrivals get_mon_resources run
#' @importFrom simmer simmer timeout trajectory

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

param <- create_params()
result <- model(param = param, run_number = 1L)
knitr::kable(result[["arrivals"]] |> arrange(start_time))
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 
knitr::kable(result[["run_results"]])
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\)

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:

  1. Calculates the time since the last update, then appends n_in_system x duration to area_n_in_system.
  2. Then updates time_last_n_insystem and n_in_system accordingly.

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)
            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: {self.env.now:.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)

    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

runner = Runner(param=Parameters())
result = runner.run_single(run=0)
print(result["patient"])
   patient_id period  run
0           1   πŸ”Ή DC    0
1           2   πŸ”Ή DC    0
2           3   πŸ”Ή DC    0
3           4   πŸ”Ή DC    0
4           5   πŸ”Ή DC    0
5           6   πŸ”Ή DC    0
6           7   πŸ”Ή DC    0
print(result["run"])
{'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.
#'
#' @return Tibble with column containing mean number of patients in the system.
#' @export

calc_mean_patients_in_system <- function(patient_count) {
  patient_count |>
    # Sort by time
    arrange(.data[["time"]]) |>
    # Calculate time between this row and the next
    mutate(interval_duration = (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().

#' 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_patients_in_system(results[["patients_in_system"]])
  )
  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).

#' Run simulation.
#'
#' @param param List. Model parameters.
#' @param run_number Numeric. Run number for random seed generation.
#'
#' @importFrom simmer add_generator get_mon_arrivals get_mon_resources run
#' @importFrom simmer simmer timeout trajectory

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"]])

  # 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)

  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"]].

param <- create_params()
result <- model(param = param, run_number = 1L)
knitr::kable(result[["patients_in_system"]])
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 
knitr::kable(result[["run_results"]])
replication mean_patients_in_system
1 1.141111

πŸ™ˆ 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

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.
    time_with_doctor : float
        Time spent in consultation with a doctor (minutes).
    """
    def __init__(self, patient_id, period):
        """
        Initialises a new patient.

        Parameters
        ----------
        patient_id : int
            Unique patient identifier.
        period : str
            Arrival period (warm up or data collection) with emoji.
        """
        self.patient_id = patient_id
        self.period = period
        self.time_with_doctor = np.nan

Model 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)
            self.patients.append(patient)

            # Print arrival time
            if self.param.verbose:
                print(f"{patient.period} Patient {patient.patient_id} " +
                      f"arrives at time: {self.env.now:.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)

    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.

runner = Runner(param=Parameters())
result = runner.run_single(run=0)
print(result["patient"])
   patient_id period  time_with_doctor  run
0           1   πŸ”Ή DC         19.610317    0
1           2   πŸ”Ή DC          9.489737    0
2           3   πŸ”Ή DC         41.665475    0
3           4   πŸ”Ή DC          5.874479    0
4           5   πŸ”Ή DC         27.882444    0
5           6   πŸ”Ή DC         24.914615    0
6           7   πŸ”Ή DC               NaN    0
print(result["run"])
{'run_number': 0, 'unseen_count': np.int64(1)}

Backlogged patient mean wait time

Patient class

Here, we add new attributes just as we have done for other measures:

  • arrival_time (like for mean time in system).
  • time_with_doctor (like for mean time in consultation and unseen patient count).
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 that patient arrives (minutes).
    time_with_doctor : float
        Time spent in consultation with a doctor (minutes).
    """
    def __init__(self, patient_id, period):
        """
        Initialises a new patient.

        Parameters
        ----------
        patient_id : int
            Unique patient identifier.
        period : str
            Arrival period (warm up or data collection) with emoji.
        """
        self.patient_id = patient_id
        self.period = period
        self.arrival_time = np.nan
        self.time_with_doctor = np.nan

Model class

We also repeat previous changes to Model:

  • Storing the arrival time in Model.generate_arrivals() as patient.arrival_time (like for mean time in system).
  • Storing the sampled consultation length in Model.consultation() as patient.time_with_doctor (like for mean time in consultation and unseen patient count).
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)
            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)

    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

runner = Runner(param=Parameters())
result = runner.run_single(run=0)
print(result["patient"])
   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
print(result["run"])
{'run_number': 0, 'unseen_wait_time': np.float64(5.045706616721617)}

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.
#'
#' @return Tibble with columns containing result for each resource.
#' @export

calc_unseen_n <- function(arrivals, groups = NULL) {
  arrivals |>
    group_by(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.
#'
#' @return Tibble with columns containing result for each resource.
#' @export

calc_unseen_mean <- function(arrivals) {
  arrivals |>
    group_by(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 mutate
#' @importFrom simmer add_generator get_mon_arrivals get_mon_resources run
#' @importFrom simmer set_attribute simmer timeout trajectory

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) |>
    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

param <- create_params()
result <- model(param = param, run_number = 1L)
knitr::kable(result[["arrivals"]])
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 
knitr::kable(result[["run_results"]])
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

GitHub Click to visit pydesrap_mms repository

Key files simulation/patient.py
simulation/monitoredresource.py
simulation/model.py
simulation/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.

GitHub Click to visit rdesrap_mms repository

Key files R/get_run_results.R
R/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

GitHub Click to visit pydesrap_stroke repository

Key files simulation/patient.py
simulation/model.py
simulation/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.

GitHub Click to visit rdesrap_stroke repository

Key files R/get_occupancy_stats.R
R/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.