# Class to store global parameter values.
class g:
# Inter-arrival times
patient_inter = 5
# Activity times
mean_n_consult_time = 6
# Resource numbers
number_of_nurses = 1
# Simulation meta parameters
sim_duration = 2880
warm_up_period = 1440 ##NEW - this will be in addition to the sim_duration
number_of_runs = 1009 Warm Up Periods
In the models we’ve created so far patients start coming in when the service opens, and then all leave when it closes.
But what if our system isn’t like that? What if we have a system that is never empty - like an Emergency Department?
By default, a DES model assumes that our system is empty at the start of a simulation run. But if we were modelling an ED, that would skew (throw off) our results, as the initial period during which patients were coming in to an empty system wouldn’t represent what’s happening in the real world - known as the Initialisation Bias.
The solution to this in DES modelling is to use a Warm Up Period.
The idea of a warm up period is simple. We run the model as normal - from empty - but for a period of time (the warm up period) we don’t collect results.
The model continues to run as normal, it’s just we don’t count what’s happening.
Warning
If you don’t use a warm-up period, you may find that the average waits you give are a lot lower than the true state - the average will be pulled lower by the earlier period of results before queues build up to their normal levels.
9.1 How long should a warm-up period be?
The length of the warm up period is up to you as the modeller to define.
You could be very precise about analysing it and use statistical testing to identify when the system reaches equilibrium (see https://eudl.eu/pdf/10.4108/ICST.SIMUTOOLS2009.5603 as an example).
Or you could plot what’s happening over time by eye and make an estimate.
Or you could just set your warm up period long enough that it’ll be representative when it starts collecting results.
9.2 Implementing the warm-up period
Implementing a warm up period in SimPy is really easy.
We just simply check the current time whenever we go to calculate / store a result, and see if it’s beyond the warm up period. If it is, we do it. If it’s not, we don’t.
Let’s look at an example. This is a slightly amended version of the model of patients coming in for a nurse consultation with a few tweaks (longer duration, more runs, added trial results calculation)
We’re going to assume this is a system that’s open 24 hours - let’s imagine this is a triage function at an emergency department.
I’ve marked the bits I’ve added to include warm up with ##NEW
9.2.1 The g class
First we add in a new parameter - the length of the warm-up period.
Here, the sim duration has been set to 2880, and the warm-up-period to half of this (1440). You don’t need to stick to this pattern - your warm-up could even be longer than your results collection if you want!
Tip
If you find it easier to keep track of, you could define your warm-up like this instead.
results_collection_period = 2880
warm_up_period = 1440
total_sim_duration = results_collection_period + warm_up_period9.2.2 The patient class
Our patient class is unchanged.
9.2.3 The model class
In the model class, the ‘attend_clinic’ method changes.
We look at the current elapsed simulation time with the attribute self.env.now
Then, whenever a patient attends the clinic and is using a nurse resource, we check whether the current simulation time is later than the number of time units we’ve set as our warm-up.
9.2.3.1 The attend_clinic method
# Generator function representing pathway for patients attending the
# clinic.
def attend_clinic(self, patient):
# Nurse consultation activity
start_q_nurse = self.env.now
with self.nurse.request() as req:
yield req
end_q_nurse = self.env.now
patient.q_time_nurse = end_q_nurse - start_q_nurse
##NEW - this checks whether the warm up period has passed before
# adding any results
if self.env.now > g.warm_up_period:
self.results_df.at[patient.id, "Q Time Nurse"] = (
patient.q_time_nurse
)
sampled_nurse_act_time = random.expovariate(1.0 /
g.mean_n_consult_time)
yield self.env.timeout(sampled_nurse_act_time)For example, if the simulation time is at 840 and our warm_up is 1440, this bit of code - which adds the queuing time for this patient to our records - won’t run:
self.results_df.at[patient.id, "Q Time Nurse"] = (
patient.q_time_nurse
)However, if the simulation time is 1680, for example, it will.
9.2.4 the calculate_run_results method
As we now won’t count the first patient, we need to remove the dummy first patient result entry we created when we set up the dataframe.
# Method to calculate and store results over the run
def calculate_run_results(self):
self.results_df.drop([1], inplace=True) ##NEW
self.mean_q_time_nurse = self.results_df["Q Time Nurse"].mean()9.2.4.1 The run method
Next we need to tweak the duration of our model to reflect the combination of the period we want to collect results for and the warm-up period.
# Method to run a single run of the simulation
def run(self):
# Start up DES generators
self.env.process(self.generator_patient_arrivals())
# Run for the duration specified in g class
##NEW - we need to tell the simulation to run for the specified duration
# + the warm up period if we still want the specified duration in full
self.env.run(until=(g.sim_duration + g.warm_up_period))
# Calculate results over the run
self.calculate_run_results()
# Print patient level results for this run
print (f"Run Number {self.run_number}")
print (self.results_df)9.2.5 The Trial class
Our trial class is unchanged.
9.3 The impact of the warm-up period
Let’s compare the results we get with and without the warm-up period.
9.3.1 Editing our results method
To make it easier to look at the outputs, I’m going to modify two methods slightly.
First, we modify the run method of the Model class slightly to swap from print the patient level dataframes to returning them as an output.
# Method to run a single run of the simulation
def run(self):
# Start up DES generators
self.env.process(self.generator_patient_arrivals())
# Run for the duration specified in g class
# We need to tell the simulation to run for the specified duration
# + the warm up period if we still want the specified duration in full
self.env.run(until=(g.sim_duration + g.warm_up_period))
# Calculate results over the run
self.calculate_run_results()
# Return patient level results for this run
return (self.results_df) ##NEWNext, we modify the run_trial method of the Trial class so that we get multiple outputs: the full patient level dataframes, a summary of results per trial, and an overall average figure for all of the trials.
# Method to run a trial
def run_trial(self):
# Run the simulation for the number of runs specified in g class.
# For each run, we create a new instance of the Model class and call its
# run method, which sets everything else in motion. Once the run has
# completed, we grab out the stored run results and store it against
# the run number in the trial results dataframe. We also return the
# full patient-level dataframes.
# First, create an empty list for storing our patient-level dataframes.
results_dfs = []
for run in range(g.number_of_runs):
my_model = Model(run)
patient_level_results = my_model.run()
print( self.df_trial_results)
# First let's record our mean wait time for this run
self.df_trial_results.loc[run] = [my_model.mean_q_time_nurse]
# Next let's work on our patient-level results dataframes
# We start by rounding everything to 2 decimal places
patient_level_results = patient_level_results.round(2)
# Add a new column recording the run
patient_level_results['run'] = run
# Now we're just going to add this to our empty list (or, after the first
# time we loop through, as an extra dataframe in our list)
results_dfs.append(patient_level_results)
all_results_patient_level = pd.concat(results_dfs)
# This calculates the attribute self.mean_q_time_nurse_trial
self.calculate_means_over_trial()
# Once the trial (ie all runs) has completed, return the results
return self.df_trial_results, all_results_patient_level, self.mean_q_time_nurse_trial9.3.2 The full updated code
import simpy
import random
import pandas as pd
# Class to store global parameter values.
class g:
# Inter-arrival times
patient_inter = 5
# Activity times
mean_n_consult_time = 6
# Resource numbers
number_of_nurses = 1
# Simulation meta parameters
sim_duration = 2880
number_of_runs = 20
warm_up_period = 1440 ##NEW - this will be in addition to the sim_duration
# Class representing patients coming in to the clinic.
class Patient:
def __init__(self, p_id):
self.id = p_id
self.q_time_nurse = 0
# Class representing our model of the clinic.
class Model:
# Constructor
def __init__(self, run_number):
# Set up SimPy environment
self.env = simpy.Environment()
# Set up counters to use as entity IDs
self.patient_counter = 0
# Set up resources
self.nurse = simpy.Resource(self.env, capacity=g.number_of_nurses)
# Set run number from value passed in
self.run_number = run_number
# Set up DataFrame to store patient-level results
self.results_df = pd.DataFrame()
self.results_df["Patient ID"] = [1]
self.results_df["Q Time Nurse"] = [0.0]
self.results_df.set_index("Patient ID", inplace=True)
# Set up attributes that will store mean queuing times across the run
self.mean_q_time_nurse = 0
# Generator function that represents the DES generator for patient arrivals
def generator_patient_arrivals(self):
while True:
self.patient_counter += 1
p = Patient(self.patient_counter)
self.env.process(self.attend_clinic(p))
sampled_inter = random.expovariate(1.0 / g.patient_inter)
yield self.env.timeout(sampled_inter)
# Generator function representing pathway for patients attending the
# clinic.
def attend_clinic(self, patient):
# Nurse consultation activity
start_q_nurse = self.env.now
with self.nurse.request() as req:
yield req
end_q_nurse = self.env.now
patient.q_time_nurse = end_q_nurse - start_q_nurse
##NEW - this checks whether the warm up period has passed before
# adding any results
if self.env.now > g.warm_up_period:
self.results_df.at[patient.id, "Q Time Nurse"] = (
patient.q_time_nurse
)
sampled_nurse_act_time = random.expovariate(1.0 /
g.mean_n_consult_time)
yield self.env.timeout(sampled_nurse_act_time)
# Method to calculate and store results over the run
def calculate_run_results(self):
##NEW - as we now won't count the first patient, we need to remove
# the dummy first patient result entry we created when we set up the
# dataframe
self.results_df.drop([1], inplace=True)
self.mean_q_time_nurse = self.results_df["Q Time Nurse"].mean()
# Method to run a single run of the simulation
def run(self):
# Start up DES generators
self.env.process(self.generator_patient_arrivals())
# Run for the duration specified in g class
##NEW - we need to tell the simulation to run for the specified duration
# + the warm up period if we still want the specified duration in full
self.env.run(until=(g.sim_duration + g.warm_up_period))
# Calculate results over the run
self.calculate_run_results()
# Return patient level results for this run
return (self.results_df)
# Class representing a Trial for our simulation
class Trial:
# Constructor
def __init__(self):
self.df_trial_results = pd.DataFrame()
self.df_trial_results["Run Number"] = [0]
self.df_trial_results["Mean Q Time Nurse"] = [0.0]
self.df_trial_results.set_index("Run Number", inplace=True)
# Method to calculate and store means across runs in the trial
def calculate_means_over_trial(self):
self.mean_q_time_nurse_trial = (
self.df_trial_results["Mean Q Time Nurse"].mean()
)
def run_trial(self):
# Run the simulation for the number of runs specified in g class.
# For each run, we create a new instance of the Model class and call its
# run method, which sets everything else in motion. Once the run has
# completed, we grab out the stored run results and store it against
# the run number in the trial results dataframe. We also return the
# full patient-level dataframes.
# First, create an empty list for storing our patient-level dataframes.
results_dfs = []
for run in range(g.number_of_runs):
my_model = Model(run)
patient_level_results = my_model.run()
print( self.df_trial_results)
# First let's record our mean wait time for this run
self.df_trial_results.loc[run] = [my_model.mean_q_time_nurse]
# Next let's work on our patient-level results dataframes
# We start by rounding everything to 2 decimal places
patient_level_results = patient_level_results.round(2)
# Add a new column recording the run
patient_level_results['run'] = run
# Now we're just going to add this to our empty list (or, after the first
# time we loop through, as an extra dataframe in our list)
results_dfs.append(patient_level_results)
all_results_patient_level = pd.concat(results_dfs)
# This calculates the attribute self.mean_q_time_nurse_trial
self.calculate_means_over_trial()
# Once the trial (ie all runs) has completed, return the results
return self.df_trial_results, all_results_patient_level, self.mean_q_time_nurse_trial
# Method to print trial results, including averages across runs
def print_trial_results(self):
print ("Trial Results")
# EDIT: We are omitting the printouts of the patient level data for now
# print (self.df_trial_results)
print (f"Mean Q Nurse : {self.mean_q_time_nurse_trial:.1f} minutes")
# Create new instance of Trial and run it
my_trial = Trial()
df_trial_results_warmup, all_results_patient_level_warmup, means_over_trial_warmup = my_trial.run_trial() Mean Q Time Nurse
Run Number
0 0.0
Mean Q Time Nurse
Run Number
0 411.550794
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
7 463.991249
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
7 463.991249
8 530.155833
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
7 463.991249
8 530.155833
9 682.078888
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
7 463.991249
8 530.155833
9 682.078888
10 577.423196
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
7 463.991249
8 530.155833
9 682.078888
10 577.423196
11 710.838935
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
7 463.991249
8 530.155833
9 682.078888
10 577.423196
11 710.838935
12 576.129679
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
7 463.991249
8 530.155833
9 682.078888
10 577.423196
11 710.838935
12 576.129679
13 396.039835
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
7 463.991249
8 530.155833
9 682.078888
10 577.423196
11 710.838935
12 576.129679
13 396.039835
14 516.628283
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
7 463.991249
8 530.155833
9 682.078888
10 577.423196
11 710.838935
12 576.129679
13 396.039835
14 516.628283
15 225.503592
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
7 463.991249
8 530.155833
9 682.078888
10 577.423196
11 710.838935
12 576.129679
13 396.039835
14 516.628283
15 225.503592
16 582.607593
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
7 463.991249
8 530.155833
9 682.078888
10 577.423196
11 710.838935
12 576.129679
13 396.039835
14 516.628283
15 225.503592
16 582.607593
17 490.646299
Mean Q Time Nurse
Run Number
0 411.550794
1 467.165767
2 421.307009
3 591.466487
4 587.817937
5 366.734538
6 354.616850
7 463.991249
8 530.155833
9 682.078888
10 577.423196
11 710.838935
12 576.129679
13 396.039835
14 516.628283
15 225.503592
16 582.607593
17 490.646299
18 614.716693 Mean Q Time Nurse
Run Number
0 0.0
Mean Q Time Nurse
Run Number
0 432.968159
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
7 209.933542
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
7 209.933542
8 362.505115
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
7 209.933542
8 362.505115
9 274.078193
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
7 209.933542
8 362.505115
9 274.078193
10 376.897458
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
7 209.933542
8 362.505115
9 274.078193
10 376.897458
11 452.180957
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
7 209.933542
8 362.505115
9 274.078193
10 376.897458
11 452.180957
12 482.076520
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
7 209.933542
8 362.505115
9 274.078193
10 376.897458
11 452.180957
12 482.076520
13 335.217374
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
7 209.933542
8 362.505115
9 274.078193
10 376.897458
11 452.180957
12 482.076520
13 335.217374
14 430.020577
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
7 209.933542
8 362.505115
9 274.078193
10 376.897458
11 452.180957
12 482.076520
13 335.217374
14 430.020577
15 528.598828
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
7 209.933542
8 362.505115
9 274.078193
10 376.897458
11 452.180957
12 482.076520
13 335.217374
14 430.020577
15 528.598828
16 356.191823
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
7 209.933542
8 362.505115
9 274.078193
10 376.897458
11 452.180957
12 482.076520
13 335.217374
14 430.020577
15 528.598828
16 356.191823
17 485.239630
Mean Q Time Nurse
Run Number
0 432.968159
1 234.125223
2 433.900936
3 382.891231
4 291.235410
5 419.142365
6 477.208932
7 209.933542
8 362.505115
9 274.078193
10 376.897458
11 452.180957
12 482.076520
13 335.217374
14 430.020577
15 528.598828
16 356.191823
17 485.239630
18 439.3048619.3.3 Comparing the results
9.3.3.1 Patient-level dataframes
First, let’s look at the first five rows of our patient dataframes.
Without the warm-up, our patient IDs start at 1.
9.3.3.1.1 Without warm-up
all_results_patient_level.head()| Q Time Nurse | Time with Nurse | run | |
|---|---|---|---|
| Patient ID | |||
| 1 | 0.00 | 3.71 | 0 |
| 2 | 0.73 | 0.18 | 0 |
| 3 | 0.00 | 1.87 | 0 |
| 4 | 0.00 | 0.80 | 0 |
| 5 | 0.00 | 1.49 | 0 |
9.3.3.1.2 With warm-up
With the warm-up, our patient IDs start later.
all_results_patient_level_warmup.head()| Q Time Nurse | run | |
|---|---|---|
| Patient ID | ||
| 235 | 194.20 | 0 |
| 236 | 192.72 | 0 |
| 237 | 187.42 | 0 |
| 238 | 177.10 | 0 |
| 239 | 176.15 | 0 |
9.3.3.2 Per-run results
9.3.3.2.1 Without warm-up
df_trial_results.round(2).head()| Mean Q Time Nurse | |
|---|---|
| Run Number | |
| 0 | 432.97 |
| 1 | 234.13 |
| 2 | 433.90 |
| 3 | 382.89 |
| 4 | 291.24 |
9.3.3.2.2 With warm-up
With the warm-up, our patient IDs start later.
df_trial_results_warmup.round(2).head()| Mean Q Time Nurse | |
|---|---|
| Run Number | |
| 0 | 411.55 |
| 1 | 467.17 |
| 2 | 421.31 |
| 3 | 591.47 |
| 4 | 587.82 |
9.3.3.3 Overall results
Without the warm up, our overall average wait time is
'391.52 minutes'With the warm up, our overall average wait time is
'493.98 minutes'You can see overall that the warm-up time can have a very significant impact on our waiting times!
Let’s look at this in a graph.
9.3.3.4 Results over time
import plotly.express as px
df_trial_results = df_trial_results.reset_index()
df_trial_results['Warm Up'] = 'No Warm Up'
df_trial_results_warmup = df_trial_results_warmup.reset_index()
df_trial_results_warmup['Warm Up'] = 'With Warm Up'
fig = px.histogram(
pd.concat([df_trial_results, df_trial_results_warmup]).round(2).reset_index(),
x="Warm Up",
color="Run Number", y="Mean Q Time Nurse",
barmode='group',
title='Average Queue Times per Run - With and Without Warmups')
fig.show()