# 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.
9.3 Coding the model
Tip
Throughout the code, anything new that’s been added will be followed by the comment ##NEW - so look out for that in the following code chunks.
9.3.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.3.2 The patient class
Our patient class is unchanged.
9.3.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.3.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.3.3.2 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.3.3.3 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.3.4 The trial class
Our trial class is unchanged.
9.4 The impact of the warm-up period
Let’s compare the results we get with and without the warm-up period.
9.4.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.4.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 733.547177
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
7 620.019615
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
7 620.019615
8 599.098511
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
7 620.019615
8 599.098511
9 589.739673
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
7 620.019615
8 599.098511
9 589.739673
10 492.675507
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
7 620.019615
8 599.098511
9 589.739673
10 492.675507
11 521.048975
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
7 620.019615
8 599.098511
9 589.739673
10 492.675507
11 521.048975
12 495.923618
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
7 620.019615
8 599.098511
9 589.739673
10 492.675507
11 521.048975
12 495.923618
13 348.453699
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
7 620.019615
8 599.098511
9 589.739673
10 492.675507
11 521.048975
12 495.923618
13 348.453699
14 474.269147
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
7 620.019615
8 599.098511
9 589.739673
10 492.675507
11 521.048975
12 495.923618
13 348.453699
14 474.269147
15 436.753107
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
7 620.019615
8 599.098511
9 589.739673
10 492.675507
11 521.048975
12 495.923618
13 348.453699
14 474.269147
15 436.753107
16 330.760125
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
7 620.019615
8 599.098511
9 589.739673
10 492.675507
11 521.048975
12 495.923618
13 348.453699
14 474.269147
15 436.753107
16 330.760125
17 561.358343
Mean Q Time Nurse
Run Number
0 733.547177
1 378.684014
2 755.920508
3 472.523404
4 578.379365
5 444.701755
6 434.718171
7 620.019615
8 599.098511
9 589.739673
10 492.675507
11 521.048975
12 495.923618
13 348.453699
14 474.269147
15 436.753107
16 330.760125
17 561.358343
18 622.622793 Mean Q Time Nurse
Run Number
0 0.0
Mean Q Time Nurse
Run Number
0 265.826637
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
7 568.610312
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
7 568.610312
8 507.310566
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
7 568.610312
8 507.310566
9 352.069080
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
7 568.610312
8 507.310566
9 352.069080
10 401.147111
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
7 568.610312
8 507.310566
9 352.069080
10 401.147111
11 342.079000
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
7 568.610312
8 507.310566
9 352.069080
10 401.147111
11 342.079000
12 305.569948
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
7 568.610312
8 507.310566
9 352.069080
10 401.147111
11 342.079000
12 305.569948
13 493.448464
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
7 568.610312
8 507.310566
9 352.069080
10 401.147111
11 342.079000
12 305.569948
13 493.448464
14 556.372446
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
7 568.610312
8 507.310566
9 352.069080
10 401.147111
11 342.079000
12 305.569948
13 493.448464
14 556.372446
15 318.082947
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
7 568.610312
8 507.310566
9 352.069080
10 401.147111
11 342.079000
12 305.569948
13 493.448464
14 556.372446
15 318.082947
16 195.727161
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
7 568.610312
8 507.310566
9 352.069080
10 401.147111
11 342.079000
12 305.569948
13 493.448464
14 556.372446
15 318.082947
16 195.727161
17 402.771738
Mean Q Time Nurse
Run Number
0 265.826637
1 433.884590
2 402.346533
3 270.089347
4 437.540951
5 459.819779
6 315.810715
7 568.610312
8 507.310566
9 352.069080
10 401.147111
11 342.079000
12 305.569948
13 493.448464
14 556.372446
15 318.082947
16 195.727161
17 402.771738
18 464.2269599.4.3 Comparing the results
9.4.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.4.3.1.1 Without warm-up
all_results_patient_level.head()| Q Time Nurse | Time with Nurse | run | |
|---|---|---|---|
| Patient ID | |||
| 1 | 0.0 | 1.15 | 0 |
| 2 | 0.0 | 2.15 | 0 |
| 3 | 0.0 | 4.93 | 0 |
| 4 | 0.0 | 0.00 | 0 |
| 5 | 0.0 | 2.89 | 0 |
9.4.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 | ||
| 213 | 413.59 | 0 |
| 214 | 425.32 | 0 |
| 215 | 426.62 | 0 |
| 216 | 421.55 | 0 |
| 217 | 428.35 | 0 |
9.4.3.2 Per-run results
9.4.3.2.1 Without warm-up
df_trial_results.round(2).head()| Mean Q Time Nurse | |
|---|---|
| Run Number | |
| 0 | 265.83 |
| 1 | 433.88 |
| 2 | 402.35 |
| 3 | 270.09 |
| 4 | 437.54 |
9.4.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 | 733.55 |
| 1 | 378.68 |
| 2 | 755.92 |
| 3 | 472.52 |
| 4 | 578.38 |
9.4.3.3 Overall results
Without the warm up, our overall average wait time is
'392.04 minutes'With the warm up, our overall average wait time is
'515.56 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.4.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()