In everyday life we often encounter and engage with queue based systems. This is the case where a system or process necessitates that entities, such as customers or patients, flow through the system while queuing for limited resources. We also often find ourselves wondering if the system we are engaged with is optimal, or speculate about the amount of time it will take to be served.
Within the domain of queuing theory, Discrete Event Simulation (DES) is a method by which we can address these issue in a programatic and mathematical way.
I previously wrote an introduction to DES, which outlines some of the key ideas behind the methodology, see here.
Here, we will go through an implementation of DES using the Python package SimPy. Further, we will demonstrate on object-orientated approach, highly useful as our models grow in complexity.
Let’s Propose a Problem
Put simply, DES can help us optimise any system characterised by entities queuing to request a particular service. This is applicable in numerous fields, such as manufacturing, telecommunications and healthcare. A typically highly pressurised environment, emergency departments offer a good use-case to test this method. Below, we have a process flow diagram depicting a simplified emergency department system. The entities, in this case patients, flow into the system, are triaged, and subsequently go through one of two processes before either going home, or being admitted into hospital.
Note that this proposed system can indeed be altered and added to depending on the needs and wants of the project. The level of detail required will need to be be carefully considered. An overly simplified depiction of a system will produce unrealistic results and won’t reflect reality, while we need to be careful not to excessively model aspects of a system that go beyond our scope.
There are up to four points in our diagram where a queue can develop, while certain elements are also subject to a degree of randomness. These are highlighted by red dashed boarders and are mainly associated with processes in our system. For instance, the interval time between patient arrivals, and the time spent with a nurse or doctor, can all be stochastically modelled. In other words, be characterised be a degree of randomness in our model.
The system depicted above greatly lends itself to being modelled via DES. We will utilise the Python package SimPy to do this.
Code Template
Here, we will go through a template one can follow to code a DES model in SimPy. We will utilise an object orientated approach to solve this coding problem, as this facilitates the development of user friendly, readable code, particularly as the model gets bigger.
The code will constitute three main classes, as follows:
- A class simply holding key variables used to manipulate the model, or test what-if scenarios.
- An entity class used to hold key information about entities passing through the system. In this case, a patient class.
- And a class representing our system, and the processes that characterise it.
In addition to the above, we may wish to include further classes for utility like purposes, such as tracking metrics.
Here, we will include only key elements of the code, displaying the three main classes above. The full code can be viewed here.
Parameters and Variables
We will include a class that simply acts as a placeholder for key variables and values that define our model. This allows for ease of access and manipulation. Note we have included information around resource levels and mean times, which will feed into the stochastic elements of the model.
Patients
The patient class represents patients coming into our system. It is instantiated with a patient ID and a placeholder for the time spent in our system. It will also hold two methods, one to determine the priority the patients holds, indicative of severity, and a triage decision method, utilised during the triage process to determine which pathway the patient flows down. These methods will be used at appropriate places in our main model code.
System Class
The final defining class of our code can be described as the system, or model class. This is where key processes are defined and represents the pathway(s) through which our entities flow.
I will break this class up into a few main chunks, given it is somewhat longer than the other classes. Firstly, upon instantiation…
The first important step here is to set up our SimPy environment. The environment manages the simulation, stepping through time on an event to event basis. In addition to time, the environment manages the scheduling and processing of events. These are all critical components of a DES, which you would otherwise need to code manually. SimPy thus greatly simplifies implementation.
In addition to the environment, we assign a patient counter, which also acts as the patient ID.
The resources (the yellow boxes in our flow chart) are also defined at this point. These are critical elements of our model, and largely determine how efficient our system is at processing entities. We have two types of resources defined, which are derived from functionality offered by SimPy. The first, most basic type, are the Resource class. We set a certain capacity for these resources, and they work on a first in first out (FIFO) basis. That is, entities are dealt with on a first come, first serve basis. In our system, the triage nurse is one such resource, hence self.nurse = simpy.Resource. Another type of resource SimPy allows you to define is a priority resource. These resources deal with entities based on a level of priority the entity is associated with, as opposed to a FIFO queue system. In our model, we ranked entities in terms of priority from 1 to 5, with 1 being the most critical and 5 the least. Examples of priority resources in our system are doctors and cubicles.
Following the instantiation of our key model variables, we move onto the first process function. The first process we want to define is that of arrivals into our system, and we can do that as follows.
We generate patients with the use of an infinite while loop, which iterates on the patient counter to give us a new patient ID for each patient. We then initiate the patient with our AEPatient class, defined above, and subsequently send the patient through our main emergency department processes, which we will define later. The other main component of our function is the samples_interarrival variable, which randomly generates the time to the next entity arrival, utilising an exponential distribution to sample from. This distribution is commonly used to model time between arrivals. We feed into the exponential function the mean time between arrivals, which in our case we have predetermined, but in reality you would want to determine from real data.
The arrival process is then frozen for the sampled amount of time, following which the next patient is generated. This is achieved using yield self.env.timeout.
The while loop runs until the SimPy environment ends, which we will define later.
The next process in our system class is the attend_ed process, in which we define the possible journeys a patient might take through a system.
The comments in the code below describe the pathways through which the patient might travel, and follow the flow chart defined above. I will therefore below explain some of the key components of the code, which you will notice are utilised multiple times through the processes.
In the context of SimPy, we can use a with block to request a specific resource, as shown in the code. One of the first things we do within this block is yield the request, which tells the process to request the resource, and if it is not available, the process will freeze in place until it is available (handily, SimPy takes care of this for us). In our system, once the resource is available, we typically then sample a consultation time, whether it be with the nurse or the doctor. For this purpose, we utilise a lognormal distribution. A random variable that is lognormally distributed will take only positive real values, and as such is useful for modelling process times in fields such as engineering, where negative values are not possible, as well as in finance to model income distributions. It is characterised by a right skew, in other words, with a long tail towards the right¹. This article gives a great overview. Its characteristics also make it well suited to model our consultation times. As such, you will note that our process durations utilise a custom function² that generates a random number from a lognormal distribution. We then utilise the timeout function of the SimPy environment class to pause the process for the sampled amount of time.
Notice that if a patient goes through to the emergency department or to the Minor Injury Unit (MIU), a request is first made for a cubicle. Within this with block, is another with block where a request is made for a doctor. This means that only once the patient is in a cubicle will a doctor be requested, and only once the doctor is finished will the cubicle be released, in effect, where the cubicle with clause finishes. This is a really neat way to build a relationship between two related resources. Another thing to note is as the cubicle and doctor are priority resources, we indicate what priority we are using in the resource function parameter. The end of the attend_ae function represents the end of our system, any point at which a patient can leave is depicted in our flow chart above in green.
The last method in our class is run and is shown below.
This method is what we use later to begin our environment. We start the first process, which is the generator function that produces patients and use the environment’s run function to run our simulation for a specified time.
The complete code includes additional aspects predominantly related to tracking some key metrics and plotting results. I have omitted these here for the sake of brevity.
Outside of the code blocks above, we would utilise a for loop to run the simulation a specified number of times, and include any additional elements to track progress as required. A simple example is as follows.
Outputs
In our complete code, we track some metrics that help us determine the overall performance of our system. These include average wait in minutes (across all runs) for the resources in our system, including the nurse, cubicle and doctors. In addition, we also track queue lengths throughout simulation time for some key services, in an attempt to identify bottlenecks.
Given we set a priority to each patient as they enter our system, we can also visualise metrics based on priority.
We can begin with our overall results, showing average wait times in minutes, for model parameters set as above (as a reminder, we have 3 doctors, 2 nurses and 7 accident & emergency (AE) cubicles, see above for other inputs).
# average results across all runs of our model
The average number of patients served by the system was 55
The overall average wait across all runs for a triage nurse was 31.4 minutes
The overall average wait across all runs for a cubicle was 62.3 minutes
The overall average wait across all runs for a doctor was 17.2 minutes
The overall average wait across all runs for a MIU doctor was 0.2 minutes
The mean patient time in the system across all runs was 132.5 minutesWe can see the largest wait is for a cubicle, while on average each patients spends about 130 minutes getting through the system.
There is virtually no pressure on our MIU, this is because with our default settings, our simulation isn’t utilising this unit a great deal.
A nice way to visualise bottlenecks is to chart the numbers of patients waiting for a resource at any given time in our simulator. The chart below shows this, focusing only on AE components.
We see quite a queue build up for a cubicle, which is clearly one of the bottlenecks in our system and thus contributing to the long time it takes for a patient to pass through.
Again, with our default parameters, we can also visualise the mean time in our system by patient priority.
Given the way we set up our model, most patients will be priority 3. The priority based queue, as well as our system dictating that priority 5 patients can go to the MIU or home, mean high priorities and very low priorities get processed quicker.
What-If?
Now we have a simulated system, we can test what-if scenarios, prior to real-world application. In addition, we can test operational changes in an attempt to optimise our system. This is one of the great uses of DES, and here, we will look at potential scenarios.
Our default system is producing quite an extensive total process time of on average over two hours!
Imagine we are approached by an operational manager, hoping to reduce the amount of time it takes for patients to pass through the AE department. It is known that patients are waiting a long time to access a cubicle. Will simply adding more cubicles resolve the problem?
Let’s find out by adding two more cubicles and comparing the results.
The average number of patients served by the system was 56
The overall average wait across all runs for a triage nurse was 30.1 minutes
The overall average wait across all runs for a cubicle was 59.2 minutes
The overall average wait across all runs for a doctor was 30.5 minutes
The mean patient time in the system across all runs was 137.2 minutesThis hasn’t helped much! The average wait for a cubicle has very slightly gone done, while the overall time has actually slightly increased.
However, experience leads us to think that perhaps the number of cubicles isn’t the problem, but the problems lies in how efficiently we process patients once they are in a cubicle. Perhaps you noticed that only adding more cubicles actually increased the time it takes to be seen by a doctor, as the same number of doctors now have more cubicles to get through.
We propose that the problem lies upstream of the cubicles. In other words, having more doctors, will allow us to clear cubicles at a quicker rate, and thus get patients home quicker. So lets revert the number of cubicles to their default value and instead increase the number of doctors.
Below we see the results of increasing the number of doctors by one, so we now have four doctors.
The average number of patients served by the system was 63
The overall average wait across all runs for a triage nurse was 29.6 minutes
The overall average wait across all runs for a cubicle was 44.5 minutes
The overall average wait across all runs for a doctor was 5.2 minutes
The overall average wait across all runs for a MIU doctor was 0.2 minutes
The mean patient time in the system across all runs was 117.8 minutesBy adding one doctor we have reduced average total process times by about 20 minutes. It has also relieved the pressure on the cubicles.
You may have noticed that another upstream pressure in our default system is quite significant. This being the average amount of time it takes for a patient who requires admission to be found an inpatient bed, and thus (in our system), leave an AE cubicle. Although only 30% of our priority 4 and below patients require this, they are having to wait on average 90 minutes to be allocated an inpatient bed. To demonstrate the impact of upstream pressures on our AE system, we simulate a reduction in this average time to 45 minutes. The results are show below.
The average number of patients served by the system was 72
The overall average wait across all runs for a triage nurse was 33.8 minutes
The overall average wait across all runs for a cubicle was 33.0 minutes
The overall average wait across all runs for a doctor was 10.3 minutes
The overall average wait across all runs for a MIU doctor was 0.3 minutes
The mean patient time in the system across all runs was 110.1 minutesThe average wait times for a cubicle have reduced down to about 30 minutes! remember this was at about 60 minutes with our default system. We can see below also that the number of patients queuing for cubicles throughout our simulator has also reduced, from almost 40 patients waiting to now about 20.
Quite remarkable given we have only added one AE resource to our model! Clearly, inefficiencies in the inpatient systems are having a bearing on AE performance, in our simulator. Indeed this is a well known phenomenon³.
We could go on, but we can conclude here having demonstrated the immense usefulness of SimPy in developing DES models with relative ease. To take this further, we could certainly add further complexity to our system to achieve a better reflection of a real-world emergency care system.
Try out this Streamlit app I developed for an interactive demonstration of the simulation!
Acknowledgments
All images unless otherwise noted are by the author.
Thanks to colleagues from the Peninsula Collaboration for Health Operational Research and Data Science (PenCHORD) for their support during our work in the Health Service Modelling Programme (HSMA) we are currently enrolled on. Some HSMA resources can be viewed at the following links.
GitHub repos:
GitHub – hsma5/3a_introduction_to_discrete_event_simulation
There are also some YouTube videos available: https://www.youtube.com/channel/UCCY9_Gxg6kM-xjk9vV0mzIQ/videos
References
- https://towardsdatascience.com/log-normal-distribution-a-simple-explanation-7605864fb67c
- Custom function that produces a random number from a lognormal distribution kindly shared by Thomas Monks, Associate Professor of Health Data Science at The University of Exeter.
- Spill Over Effects of Inpatient Bed Capacity on Accident and Emergency Performance in England. R Friebel and R M.Juarez.https://www.sciencedirect.com/science/article/abs/pii/S0168851020301901#!
