Introduction to Simpy: discrete event simulation in Python

The python package Simpy provides a convenient way to program discrete event simulation that incorporate resources with limited capacity. This tutorial provides somes of the basics to use the package.

Check the simpy website for instructions on how to install simpy in your computer (be careful not to confuse it with the sympy package, which is different and focus on symbolic algebra).

The tutorial is organized as follows:

Generator functions

You are probably familiar with using functions in Python. Generators are defined using similar syntax as functions but are used in a different way.

When calling a function, it is executed "from scratch" every time it is called. Consider the following function.

This is a simple function that adds all integers from 0..n and returns the result. Now, consider a variation of this function that incrementally adds the integers every time the function is called. This is an example of a generator.

The yield command is similar to a return in that is returns a value when the function is called; however, the function "remembers" the point at which it returned the value and proceed from that point on the next time it is called. The next command tells the generator to resume execution from where it was until the next yield statement is returned.

In simpy, generators are used to create entities in the simulation, which is useful to keep track of the state of the entity during the simulation as the clock moves on.

Basics of simpy

We will use two basic elements in the simpy simulation API:

For our example, consider a bank with two ATM machines. The object env is defined as the enviroment and the atm as a resources with 2 units of capacity (the two machines).

We now define an entity and how it flows through the system. These are defined through a generator that receives as input the simpy environment that will manage the flow of the entity and the ATM resources that he will use. We also provide the arrival time and the service time of the customer (how long he uses the ATM machine), and his name.

Notice the different yield statements inside the generator. The first yield tells the environment to call back the generator later in time, after a timeout of arrival_t time units. When the env calls the generator back again, it resumes from that point onwards. The is a variable of the environment object that keeps track of the clock.

Then, the generator defines req as a request of the ATM resource (atm.request()). A new yield command is called to make the request, and the enviroment will call back the generator as soon as the resource is idle (which could be inmediately). The generator resumes and then a new yield is move the clock while the customer is using the resource. After the generator resumes, the resources is released, and can therefore be used to fill the request of another customer process.

To execute the simulation, we generate customer arrivals and service times. We use the numpy package to generate these random variables according to an exponential distribution. The customer process are "scheduled" using the env.process() function of the environment object.

Finally, we tell the environment to run the simulation with all the processes that were scheduled. The output shows the arrival times, the time each customer starts service and finishes.

Example of a Queuing System

We extend the previous example to study a system with multiple queues, similar to what could be a supermarket with multiple cashiers. We also use the example to illustrate how to store relevant information of the simulation that is used as output for later analysis.

The system is composed by three cashiers with the same service time. Arrival times vary during the day, in four time blocks.

We create a list to store the changes in the occupancy of the system. By convention, -1 indicates empty cashier, 0 is server occupied with no queue and >0 is the length of queue.

Customer entities join the system and choose which cashier to join. In this simulation, customers choose the cashier with the shortest queue. The code has a similar structure to what we had before, with some additions:

In the main program, we initialize the environment and create the cashier resources.

To create arrivals, we consider the non-homogenous arrival rate during different time periods of the day.

Running the simulation... Output is quite long, so its surpressed in the display.

The next chunk uses the pandas packages to export the data collected in the simulation. A csv file is generated showing, for each cashier, the time epochs at which the occupancy of the cashier changed. by convention.

The output table, order by time, looks like this:

The following code transforms the data into a format that makes it easier to visualize in a graph. It uses the matplotliblibrary to do the graph

Analyzing the output of a simulation

Consider the first example of the ATM machine with a single server. We are interested in measuring the expected time at which the last customer finishes processing. What we showed previously is the result of one simulation run. To measure the outcome presented here, we need to run multiple simulations and calculate averages across simulations to estimate these expected values. We recycle some of the previous code and extend it to have multiple simulations.

Now we analize the output of the simulation and construct the 95% confidence interval. Here, the relative precision is defined as the half-width of the confidence interval divided by the midpoint of the interval.

Confidence Interval = $$\bar X \pm z_\alpha S/\sqrt{R}$$ where $R$ is the number of runs.

The relative precision is defined as: $$( z_\alpha S/\sqrt{R} - \bar X)/\bar X$$

Suppose we want to achieve a relative precision of 2%. How many runs should we execute to achieve this precision? The half-width of the interval has to be $d = 0.02\bar X$