Introduction to Discrete-Event Simulation and the SimPy Language
Norm Matloff
February 13, 2008
c2006-2008, N.S. Matloff
Contents
1 What Is Discrete-Event Simulation (DES)?
2 World Views in DES Programming
2.1 The Activity-Oriented Paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 The Event-Oriented Paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 The Process-Oriented Paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
Introduction to the SimPy Simulation Language
3.1 SimPy Overview .
.
.
.
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2
Introduction to SimPy Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
3
3
4
6
7
8
9
3.2.1 MachRep1.py: Our First SimPy Program . . . . . . . . . . . . . . . . . . . . . . . 10
3.2.2 MachRep2.py: Introducing the Resource Class . . . . . . . . . . . . . . . . . . . . 14
3.2.3 MachRep3.py: Introducing Passivate/Reactivate Operations
. . . . . . . . . . . . . 16
3.2.4 MMk.py: “Do It Yourself” Queue Management . . . . . . . . . . . . . . . . . . . . 18
3.2.5
SMP.py: Simultaneous Possession of Resources . . . . . . . . . . . . . . . . . . . . 20
3.2.6 Cell.py: Dynamic Creation of Threads . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.3 Note These Restrictions on PEMs
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4 SimPy Data Collection and Display . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4.1
Introduction to Monitors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4.2 Time Averages .
.
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4.3 The Function Monitor.timeAverage()
. . . . . . . . . . . . . . . . . . . . . . . . . 27
3.4.4 But I Recommend That You Not Use This Function . . . . . . . . . . . . . . . . . . 27
1
3.4.5 Little’s Rule .
.
3.5 Other SimPy Features .
.
.
.
.
.
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
A How to Obtain and Install SimPy
B Debugging and Verifying SimPy Programs
29
30
B.1 Debugging Tools
.
.
.
.
.
.
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
B.2 Know How Control Transfers in SimPy Programs . . . . . . . . . . . . . . . . . . . . . . . 30
B.3 Always Know What (Simulated) Time It Is
. . . . . . . . . . . . . . . . . . . . . . . . . . 31
B.4 Starting Over
.
B.5 Repeatability .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
B.6 Peeking at the SimPy’s Internal Event List . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
B.7 SimPy’s Invaluable Tracing Library . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
C Online Documentation for SimPy
33
2
1 What Is Discrete-Event Simulation (DES)?
Consider simulation of some system which evolves through time. There is a huge variety of such applica-
tions. One can simulate a weather system, for instance. A key point, though, is that in that setting, the events
being simulated would be continuous, meaning for example that if we were to graph temperature against
time, the curve would be continuous, no breaks.
By contrast, suppose we simulate the operation of a warehouse. Purchase orders come in and are filled,
reducing inventory, but inventory is replenished from time to time. Here a typical variable would be the
inventory itself, i.e. the number of items currently in stock for a given product. If we were to graph that
number against time, we would get what mathematicians call a step function, i.e. a set of flat line seg-
ments with breaks between them. The events here—decreases and increases in the inventory—are discrete
variables, not continuous ones.
DES involves simulating such systems.
2 World Views in DES Programming
Simulation programming can often be difficult—difficult to write the code, and difficult to debug. The
reason for this is that it really is a form of parallel programming, with many different activities in progress
simultaneously, and parallel programming can be challenging.
For this reason, many people have tried to develop separate simulation languages, or at least simulation
paradigms (i.e. programming styles) which enable to programmer to achieve clarity in simulation code.
Special simulation languages have been invented in the past, notably SIMULA, which was invented in the
1960s and has significance today in that it was the language which invented the concept of object-oriented
programmg that is so popular today. However, the trend today is to simply develop simulation libraries
which can be called from ordinary languages such as C++, instead of inventing entire new languages.1 So,
the central focus today is on the programming paradigms, not on language. In this section we will present
an overview of the three major discrete-event simulation paradigms.
Several world views have been developed for DES programming, as seen in the next few sections.
2.1 The Activity-Oriented Paradigm
Let us think of simulating a queuing system. Jobs arrive at random times, and the job server takes a ran-
dom time for each service. The time between arrivals of jobs, and the time needed to serve a job, will be
continuous random variables, possibly having exponential or other continuous distributions.
For concreteness, think of an example in which the server is an ATM cash machine and the jobs are cus-
tomers waiting in line.
Under the activity-oriented paradigm, we would break time into tiny increments. If for instance the mean
interarrival time were, say 20 seconds, we might break time into increments of size 0.001. At each time
point, our code would look around at all the activities, e.g. currently-active job servicing, and check for the
possible occurrence of events, e.g. completion of service. Our goal is to find the long-run average job wait
1These libraries are often called “languages” anyway, and I will do so too.
3
time.
Let SimTime represent current simulated time. Our simulation code in the queue example above would look
something like this:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
QueueLength = 0
NJobsServed = 0
SumResidenceTimes = 0
ServerBusy = false
generate NextArrivalTime // random # generation
NIncrements = MaxSimTime / 0.001
for SimTime = 1*0.001 to NIncrements*0.001 do
if SimTime = NextArrivalTime then
add new jobobject to queue
QueueLength++
generate NextArrivalTime // random # generation
if not ServerBusy then
ServerBusy = true
jobobject.ArrivalTime = SimTime
generate ServiceFinishedtime
currentjob = jobobject
delete head of queue and assign to currentjob
QueueLength--
else
if SimTime = ServiceFinishedtime then
NJobsServed++
SumResidenceTimes += SimTime - currentjob.ArrivalTime
if QueueLength > 0 then
// random # generation
generate ServiceFinishedtime
delete currentjob from queue
QueueLength--
else
ServerBusy = false
print out SumResidenceTimes / NJobsServed
2.2 The Event-Oriented Paradigm
Clearly, an activity-oriented simulation program is going to be very slow to execute. Most time increments
will produce no state change to the system at all, i.e. no new arrivals to the queue and no completions of
service by the server. Thus the activity checks will be wasted processor time. This is a big issue, because
in general simulation code often needs a very long time to run. (Electronic chip manufacturers use DES for
chip simulation. A simulation can take days to run.)
Inspection of the above pseudocode, though, shows a way to dramatically increase simulation speed. Instead
of having time “creep along” so slowly, why not take a “shortcut” to the next event? What we could do is
something like the following:
Instead of having the simulated time advance via the code
1
for SimTime = 1*0.001 to NIncrements*0.001 do
we could advance simulated time directly to the time of the next event:
4
1
2
3
4
5
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
if ServerBusy and NextArrivalTime < ServiceFinishedtime or
not ServerBusy then
SimTime = NextArrivalTime
else
SimTime = ServiceFinishedtime
(The reason for checking ServerBusy is that ServiceFinishedtime will be undefined if ServerBusy is false.)
The entire pseudocode would then be
QueueLength = 0
NJobsServed = 0
SumResidenceTimes = 0
ServerBusy = false
generate NextArrivalTime
SimTime = 0.0;
while (1) do
if ServerBusy and NextArrivalTime < ServiceFinishedtime or
not ServerBusy then
SimTime = NextArrivalTime
else
SimTime = ServiceFinishedtime
if SimTime > MaxSimTime then break
if SimTime = NextArrivalTime then
QueueLength++
generate NextArrivalTime
if not ServerBusy then
ServerBusy = true
jobobject.ArrivalTime = SimTime
currentjob = jobobject
generate ServiceFinishedtime
QueueLength--
else // the case SimTime = ServiceFinishedtime
NJobsServed++
SumResidenceTimes += SimTime - currentjob.ArrivalTime
if QueueLength > 0 then
generate ServiceFinishedtime
QueueLength--
else
ServerBusy = false
print out SumResidenceTimes / NJobsServed
The event-oriented paradigm formalizes this idea. We store an event set, which is the set of all pending
events. In our queue example above, for instance, there will always be at least one event pending, namely
the next arrival, and sometimes a second pending event, namely the completion of a service. Our code above
simply inspects the scheduled event times of all pending events (again, there will be either one or two of
them in our example here), and updates SimTime to the minimum among them.
In the general case, there may be many events in the event set, but the principle is still the same—in each
iteration of the while loop, we update SimTime to the minimum among the scheduled event times. Note
also that in each iteration of the while loop, a new event is generated and added to the set; be sure to look at
the pseudocode above and verify this.
5
Thus a major portion of the execution time for the program will consist of a find-minimum operation within
the event set. Accordingly, it is desirable to choose a data structure for the set which will facilitate this
operation, such as a heap-based priority queue. In many event-oriented packages, though, the event set is
implemented simply as a linearly-linked list. This will be sufficiently efficient as long as there usually aren’t
too many events in the event set; again, in the queue example above, the maximum size of the event set is 2.
(We will return to the issue of efficient event lists in a later unit.)
Again, note the contrast between this and continuous simulation models. The shortcut which is the heart
of the event-oriented paradigm was only possible because of the discrete nature of system change. So this
paradigm is not possible in models in which the states are continuous in nature.
The event-oriented paradigm was common in the earlier years of simulation, used in packages in which code
in a general-purpose programming language such as C called functions in a simulation library. It still has
some popularity today. Compared to the main alternative, the process-oriented paradigm, the chief virtues
of the event-oriented approach are:
• Ease of implementation. The process-oriented approach requires something like threads, and in those
early days there were no thread packages available. One needed to write one’s own threads mecha-
nisms, by writing highly platform-dependent assembly-language routines for stack manipulation.
• Execution speed. The threads machinery of process-oriented simulation really slows down execution
speed (even if user-level threads are used).
• Flexibility. If for example one event will trigger two others, it is easy to write this into the application
code.
2.3 The Process-Oriented Paradigm
Here each simulation activity is modeled by a process. The idea of a process is similar to the notion by
the same name in Unix, and indeed one could write process-oriented simulations using Unix processes.
However, these would be inconvenient to write, difficult to debug, and above all they would be slow.
As noted earlier, the old process-oriented software such as SIMULA and later CSIM were highly platform-
dependent, due to the need for stack manipulation. However, these days this problem no longer exists, due
to the fact that modern systems include threads packages (e.g. pthreads in Unix, Java threads, Windows
threads and so on). Threads are sometimes called “lightweight” processes.
If we were to simulate a queuing system as above, but using the process-oriented paradigm, we would have
two threads, one simulating the arrivals and the other simulating the operation of the server. Those would
be the application-specific threads (so NumActiveAppThreads = 2 in the code below), and we would also
have a general thread to manage the event set.
Our arrivals thread would look something like
NumActiveAppThreads++
while SimTime < MaxSimTime do
generate NextArrivalTime
add an arrival event for time NextArrivalTime to the event set
sleep until wakened by the event-set manager
jobobject.ArrivalTime = SimTime
1
2
3
4
5
6
6
7
8
1
2
3
4
5
6
7
8
9
10
11
12
1
2
3
4
5
6
1
2
3
4
5
6
7
add jobobject to the machine queue
thread exit
The server thread would look something like
NumActiveAppThreads++
while SimTime < MaxSimTime do
sleep until QueueLength > 0
while QueueLength > 0 do
remove queue head and assign to jobobject
QueueLength--
generate ServiceFinishedtime
add a service-done event for time ServiceFinishedtime to the event set
sleep until wakened by the event-set manager
SumResidenceTimes += SimTime - jobobject.ArrivalTime
NJobsServed++
thread exit
The event set manager thread would look something like
while SimTime < MaxSimTime do
sleep until event set is nonempty
delete the minimum-time event E from the event set
update SimTime to the time scheduled for E
wake whichever thread had added E to the event set
thread exit
The function main() would look something like this:
QueueLength = 0
NJobsServed = 0
SumResidenceTimes = 0
ServerBusy = false
start the 3 threads
sleep until all 3 threads exit
print out SumResidenceTimes / NJobsServed
Note that the event set manager would be library code, while the other modules shown above would be
application code.
Two widely used oper-source process-oriented packages are C++SIM, available at http://cxxsim.
ncl.ac.uk and SimPy, available at http://simpy.sourceforge.net.
The process-oriented paradigm produces more modular code. This is probably easier to write and easier for
others to read. It is considered more elegant, and is the more popular of the two main world views today.
3 Introduction to the SimPy Simulation Language
SimPy (rhymes with “Blimpie”) is a package for process-oriented discrete-event simulation. It is written in,
and called from, Python. I like the clean manner in which it is designed, and the use of Python generators—
7
and for that matter, Python itself—is a really strong point. If you haven’t used Python before, you can learn
enough about it to use SimPy quite quickly; see my quick introduction to Python, at my Python tutorials
page, http://heather.cs.ucdavis.edu/˜matloff/python.html.
Instructions on how to obtain and install SimPy are given in Appendix A.
Instead of using threads, as is the case for most process-oriented simulation packages, SimPy makes novel
use of Python’s generators capability.2 Generators allow the programmer to specify that a function can be
prematurely exited and then later re-entered at the point of last exit, enabling coroutines, meaning functions
that alternate execution with each other. The exit/re-entry points are marked by Python’s yield keyword.
Each new call to the function causes a resumption of execution of the function at the point immediately
following the last yield executed in that function. As you will see below, that is exactly what we need for
DES.
For convenience, I will refer to each coroutine (or, more accurately, each instance of a coroutine), as a
thread.3
3.1 SimPy Overview
Here are the major SimPy classes which we will cover in this introduction:4
• Process: simulates an entity which evolves in time, e.g. one customer who needs to be served by an
ATM machine; we will refer to it as a thread, even though it is not a formal Python thread
• Resource: simulates something to be queued for, e.g. the machine
Here are the major SimPy operations/function calls we will cover in this introduction:
• activate(): used to mark a thread as runnable when it is first created
• simulate(): starts the simulation
• yield hold: used to indicate the passage of a certain amount of time within a thread; yield is a Python
operator whose first operand is a function to be called, in this case a code for a function that performs
the hold operation in the SimPy library
• yield request: used to cause a thread to join a queue for a given resource (and start using it immedi-
ately if no other jobs are waiting for the resource)
• yield release: used to indicate that the thread is done using the given resource, thus enabling the next
thread in the queue, if any, to use the resource
• yield passivate: used to have a thread wait until “awakened” by some other thread
2Python 2.2 or better is required. See my Python generators tutorial at the above URL if you wish to learn about generators, but
you do not need to know about them to use SimPy.
3This tutorial does not assume the reader has a background in threads programming. In fact, readers who do have that back-
ground will have to unlearn some of what they did before, because our threads here will be non-preemptive, unlike the preemptive
type one sees in most major threads packages.
4Others will be covered in our followup tutorial at AdvancedSimpy.pdf.
8