# 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 604.392493
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
7 697.157028
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
7 697.157028
8 221.530923
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
7 697.157028
8 221.530923
9 603.924123
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
7 697.157028
8 221.530923
9 603.924123
10 522.444605
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
7 697.157028
8 221.530923
9 603.924123
10 522.444605
11 410.161875
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
7 697.157028
8 221.530923
9 603.924123
10 522.444605
11 410.161875
12 425.966603
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
7 697.157028
8 221.530923
9 603.924123
10 522.444605
11 410.161875
12 425.966603
13 507.233729
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
7 697.157028
8 221.530923
9 603.924123
10 522.444605
11 410.161875
12 425.966603
13 507.233729
14 503.829046
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
7 697.157028
8 221.530923
9 603.924123
10 522.444605
11 410.161875
12 425.966603
13 507.233729
14 503.829046
15 452.236705
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
7 697.157028
8 221.530923
9 603.924123
10 522.444605
11 410.161875
12 425.966603
13 507.233729
14 503.829046
15 452.236705
16 681.681550
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
7 697.157028
8 221.530923
9 603.924123
10 522.444605
11 410.161875
12 425.966603
13 507.233729
14 503.829046
15 452.236705
16 681.681550
17 281.565863
Mean Q Time Nurse
Run Number
0 604.392493
1 392.966328
2 499.065886
3 847.607318
4 692.436899
5 534.788088
6 404.362230
7 697.157028
8 221.530923
9 603.924123
10 522.444605
11 410.161875
12 425.966603
13 507.233729
14 503.829046
15 452.236705
16 681.681550
17 281.565863
18 532.896310 Mean Q Time Nurse
Run Number
0 0.0
Mean Q Time Nurse
Run Number
0 548.990807
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
7 410.618491
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
7 410.618491
8 522.823278
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
7 410.618491
8 522.823278
9 423.640239
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
7 410.618491
8 522.823278
9 423.640239
10 283.635002
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
7 410.618491
8 522.823278
9 423.640239
10 283.635002
11 275.398061
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
7 410.618491
8 522.823278
9 423.640239
10 283.635002
11 275.398061
12 340.053119
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
7 410.618491
8 522.823278
9 423.640239
10 283.635002
11 275.398061
12 340.053119
13 307.091205
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
7 410.618491
8 522.823278
9 423.640239
10 283.635002
11 275.398061
12 340.053119
13 307.091205
14 523.840251
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
7 410.618491
8 522.823278
9 423.640239
10 283.635002
11 275.398061
12 340.053119
13 307.091205
14 523.840251
15 420.605314
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
7 410.618491
8 522.823278
9 423.640239
10 283.635002
11 275.398061
12 340.053119
13 307.091205
14 523.840251
15 420.605314
16 605.624993
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
7 410.618491
8 522.823278
9 423.640239
10 283.635002
11 275.398061
12 340.053119
13 307.091205
14 523.840251
15 420.605314
16 605.624993
17 447.096471
Mean Q Time Nurse
Run Number
0 548.990807
1 430.365736
2 418.413786
3 342.890572
4 199.991928
5 370.084390
6 463.615546
7 410.618491
8 522.823278
9 423.640239
10 283.635002
11 275.398061
12 340.053119
13 307.091205
14 523.840251
15 420.605314
16 605.624993
17 447.096471
18 451.1269119.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 | 1.42 | 0 |
| 2 | 0.00 | 16.20 | 0 |
| 3 | 15.47 | 4.67 | 0 |
| 4 | 13.25 | 4.30 | 0 |
| 5 | 11.79 | 14.66 | 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 | ||
| 217 | 345.05 | 0 |
| 218 | 348.83 | 0 |
| 219 | 349.12 | 0 |
| 220 | 355.78 | 0 |
| 221 | 376.72 | 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 | 548.99 |
| 1 | 430.37 |
| 2 | 418.41 |
| 3 | 342.89 |
| 4 | 199.99 |
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 | 604.39 |
| 1 | 392.97 |
| 2 | 499.07 |
| 3 | 847.61 |
| 4 | 692.44 |
9.3.3.3 Overall results
Without the warm up, our overall average wait time is
'403.89 minutes'With the warm up, our overall average wait time is
'510.61 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()