 |
|
|
|
|||||||||||||||
|
 |
|
||||||||||||||||
|
||||||||||||||||||
 |
 |
Issue 7 - Revision 6 / February 27, 2005
|
|||
|
|
|
|
|
|
|
|
Simulating with SimPy - - - - - - - - - - - - By Dr. Klaus G. Müller | January 16, 2005 
The behavior and performance of a complex system influenced by internal and environmental factors can often not be predicted analytically. In some cases, experimentation can yield useful results, but in other cases, experimentation with the actual system may be undesirable or even impossible, because of risks or costs, or because the system is still under development. One useful alternative to experimentation is simulation, building an approximation, or model, of a complex system, including its resources, processes, and rules, and executing the model on a computer. While such models are typically simplifications of the real system, addressing only those aspects that are relevant to the issue being studied, simulation is often equivalent to performing a sequence of experiments. Good simulation software packages aid the model builder by presenting a strong conceptual framework for modeling a wide class of phenomena. Desireable features of a simulation package include: a notation to describe dynamic behaviors; programming tools that help prevent and rapidly detect and correct errors; and collection tools that capture, analyze, and display statistical data. Additionally, simulation software should facilitate the generation of stochastic variates from a wide range of statistical distributions ("throwing dice") and should provide insight into system processes by tracing the dynamic behavior of model programs. System simulations can be continuous simulations or discrete event simulations. In the former, simulated time progresses by small, constant increments, and state variables change incrementally. In the latter, state variables change by arbitrary amounts only at specified time points. Events are the instantaneous change of one or more state variables. Which type of simulation to use depends on the nature of the system being studied.
A Tasty Bite of SimPy
SimPy (short for Simulation in Python, available for download from http://sourceforge.net/projects/simpy/) is an open source, discrete event simulation package written completely in Python. (In fact, it's the only such simulation package in Python.) Used in many universities, research organizations, and industry for a wide range of activities, SimPy allows the investigation and optimization of existing or planned complex systems. SimPy satisfies all of the design criteria listed above. It provides the system modeler with key components, such as processes and resources, for building discrete event simulation models. SimPy includes modules for data collection, analysis, plotting, and for developing simulation graphical user interfaces (GUIs). Furthermore, SimPy can easily be extended, as can models implemented in SimPy, allowing the development of application libraries and model families. Key to SimPy (and to any discrete event simulation package) is its sequencing mechanism, the way control is passed between model parts as simulated time advances. SimPy uses the process interaction approach. Model segments representing processes in the real world are semi-coroutines that are reactivated by the SimPy supervisor (SimPy's simulate() function) at specified epochs (events). Following David Mertz's great idea for Python microthreads (http://www-106.ibm.com/developerworks/linux/library/l-pythrd.html), these semi-coroutines are implemented as Python generators, passing control back to the supervisor by yield statements with a cargo. Effectively, the supervisor interprets a simulation language. Figure One shows the concept.
Figure One: Event sequencing in SimPy, using semi-coroutines At the start of a simulation run, the simulation clock is set to 0. Processes are scheduled to be activated at a time in the future, and that time is recorded in a list sorted by activation time, the event list. In Figure One, three activated processes are shown. The supervisor looks up the process with the earliest activation time, sets the clock to that time and passes control to that process. Say that the selected process is the left-most one. The supervisor effectively does a call equivalent to driveprocess.Next(), where driveprocess is the reference to the generator instance. That executes the generator code to the first yield statement encountered, changing the process' state variable speed to 50. The yield statement returns the tuple hold,self,25, a wait command. The supervisor schedules the next activation of this process 25 time units in the future and looks up the next process scheduled to activate. This goes on until either all processes have terminated or simulation time has reached endtime.
The SimPy API
The SimPy simulation library provides three abstract classes and a number of commands and functions to control process interactions and simulation execution. In addition, there are modules for tracing, plotting, real-time synchronization, and GUI-building. The first abstract class is class Process. It allows the modeling of dynamic real world phenomena. Process instances can be activated, passivated, or made to wait for a period or an event. The next abstract class is class Resource. It supports the modeling of an individual resource or a group of resources that a Process requires to complete its task. A Resource instance limits access to a resource to one process at a time. If more than one Process requests access to a resource in use, the processes must queue to wait. (In SimPy, the requesting process alone determines how long it wants to retain a resource; the resource itself is just passive.) Finally, there is class Monitor. Monitors are not essential to simulation, but they're very convenient for collecting data dynamically and providing basic statistics on that data. Monitors can also be plotted easily. All process interactions in SimPy, such as waiting, requesting and releasing resources, and so on, are programmed by yield statements with a function-specific keyword in the cargo, like hold, request, or release. Whenever a yield statement is executed, the simulate() functions activates the next event in time sequence. This should be enough of a look at SimPy's basic design. Let's take a simple scenario to illustrate the use of SimPy.
Simulating a Pizza Shop
Joe Python, a software engineer, has a mid-life crisis (after yet another successful, but tiring Python project). He wants to become a business man in a completely different field, far away from software. Since Joe likes pizza, as does every friend and colleague of his, and since Joe knows that no good software gets written without lots of pizza and Red Bull, he conceives of opening a pizza takeaway aimed at computer folks. He already has a name for it, "PizzaPy." However, before investing any time or money, Joe wants to understand the pizza business better, so he decides to simulate PizzaPy. After spying in and around his favorite pizza places and collecting some data on the market and on pizza making, he comes up with the following scenario:
Some of the other important variables are:
Joe wants to investigate the influence of the various parameters on his planned business. He looks for two figures of merit to come out of his simulation: the number of pizzas sold per day and the daily number of lost sales caused by customers leaving without ordering because of excessive wait times. JJoe maps this simple scenario into two SimPy Process sub-classes: class CustomerArrival(Process) and class Pizza(Process). The former generates the stream of customers arriving at the shop. As there is only one stream of customers in this scenario, there will only be one instance of this class. The latter class models what happens to a pizza order, and there will be as many instances of Pizza as there are pizza orders. As pizza orders result from customer arrivals, the Pizza instances can simply be generated in class CustomerArrival. The SimPy Process class is the component that models any active process. Classes derived from it must have one generator method, the Process Execution Method (PEM), which models the process dynamics. Joe models customers, pizzas, staff, oven, and lost sales as passive entities. (Other design decisions are also possible, such as representing staff members or the oven by Process sub-classes and treating pizza orders as passive entities being processed.) To get a feeling for the maximum number of sales his shop could achieve, Joe puts together a first SimPy program (shown below and found in file PizzaPy0.py), which totally ignores any resource constraints, such as staff or oven size.
1 #PizzaPy0.py
2 from SimPy.SimulationTrace import *
3 import random
4
5 class Params:pass
6 params=Params()
7 #defaults (all times in minutes)
8 params.open_time=480.0
9 params.order_time=2.0
10 params.topping_time=2.0
11 params.loading_time=0.5
12 params.baking_time=10.0
13 params.pack_time=1.0
14 params.inter_arrival=10.0
15 params.nr_repetitions=100
16
17 class Pizza(Process):
18 def __init__(self):
19 Process.__init__(self,name="Pizza order %s"%Pizza.id)
20 Pizza.id+=1
21
22 def pizza_process(self):
23 yield hold,self,params.order_time
24 yield hold,self,params.topping_time
25 yield hold,self,params.loading_time
26 yield hold,self,params.baking_time
27 yield hold,self,params.pack_time
28 Pizza.done+=1
29
30 class CustomerArrivals(Process):
31 def __init__(self):
32 Process.__init__(self,name="CustomerArrivals")
33 def arrival_process(self):
34 while True:
35 if 120<=now()<=180: #lunch rush hour
36 inter_time=params.inter_arrival/4.0
37 else:
38 inter_time=params.inter_arrival
39 arrival_delay=random.expovariate(1.0/inter_time)
40 if now()+arrival_delay<=params.open_time:
41 yield hold,self,arrival_delay
42 p=Pizza()
43 activate(p,p.pizza_process())
44 else:
45 break
46
47 def model():
48 obs_pizzas=0
49 print '\nSimulation parameters:\n',params.__dict__
50 for reps in range(params.nr_repetitions):
51 Pizza.id=0
52 Pizza.done=0
53 initialize()
54 c=CustomerArrivals()
55 activate(c,c.arrival_process())
56 simulate(until=params.open_time+30)
57 obs_pizzas+=Pizza.done
58 print "Results:\nMean number of pizzas per day: %s (%s runs)"\
59 %(obs_pizzas/params.nr_repetitions,params.nr_repetitions)
60
61 params.open_time=10.0
62 params.nr_repetitions=1
63 model()
Here's what's going on in the program: Line 2 imports the simulation capabilities with automatic tracing to help with initial debugging. (To suppress tracing, SimPy.Simulation should be imported. Either one or the other is always needed.) Lines 17-27 define the Pizza order process class, with pizza_process() as the Process Execution Method, or PEM (a generator). In lines 22-26, each yield statement causes a simulated delay of the size given by the last parameter in the statement, representing the consecutive activities required to make a pizza. Lines 30-42 define the CustomerArrivals process class, with arrival_process() as the PEM. Line 38 generates a random number from a negative exponential distribution, representing the stochastic inter-arrival times. Line 41 generates one instance of the pizza order process and activates it at the current simulation time.. Line 51 initializes the SimPy simulation machinery. It sets the simulation time to 0. Lines 52 and 53 generate the customer arrival process and activates it. Line 54 starts the simulation and runs it until either the simulated time has reached 30 minutes past closing time or until no more events are scheduled. To get a look into the dynamics of his model implementation, Joe runs the model in SimulationTrace mode, for an opening time of 10 minutes and for one repetition only. The result looks like this:
>>>python PizzaPy0.py
Simulation parameters:
{'open_time': 10, 'order_time': 2.0, 'topping_time': 2.0, 'pack_time': 1.0,
'inter_arrival': 10.0,'loading_time': 0.5, 'baking_time': 10.0,
'nr_repetitions': 1}
0 activate <CustomerArrivals> at time: 0 prior: 0
0 hold <CustomerArrivals> delay: 0.95747783121
0.95747783121 activate <Pizza order 0> at time: 0.95747783121 prior: 0
0.95747783121 hold <CustomerArrivals> delay: 2.78813992697
0.95747783121 hold <Pizza order 0> delay: 2.0
2.95747783121 hold <Pizza order 0> delay: 2.0
3.74561775818 activate <Pizza order 1> at time: 3.74561775818 prior: 0
3.74561775818 <CustomerArrivals> terminated
3.74561775818 hold <Pizza order 1> delay: 2.0
4.95747783121 hold <Pizza order 0> delay: 0.5
5.45747783121 hold <Pizza order 0> delay: 10.0
5.74561775818 hold <Pizza order 1> delay: 2.0
7.74561775818 hold <Pizza order 1> delay: 0.5
8.24561775818 hold <Pizza order 1> delay: 10.0
15.4574778312 hold <Pizza order 0> delay: 1.0
16.4574778312 <Pizza order 0> terminated
18.2456177582 hold <Pizza order 1> delay: 1.0
19.2456177582 <Pizza order 1> terminated
Results:
Mean number of pizzas per day: 2 (1 runs)
(The result will look different every time, because of the stochastic inter-arrival times.) The trace shows the activation of the CustomerArrivals process and two arrivals, resulting in two pizza orders, the progress of which can be followed step-by-step. At about t =3.75 minutes, the CustomerArrival process terminates, because the next arrival would be after closing time.
Pizza for Python Programmers, Please
By changing line 2 and the last few lines, Joe next sets up the model without tracing for the expected (default), optimistic, and pessimistic inter-arrival rates, for the full opening time of 8 hours, and for 100 repetitions. (The code is shown below and can be found in PizzaPy1.py). As Joe's more interested in the expected number of pizza sales, rather than in what happens on one random day, he needs many repetitions to get a statistically meaningful result.
#PizzaPy0.py
from SimPy.Simulation import *
import random
class Params:pass
params=Params()
#defaults (all times in minutes)
params.open_time=480.0
params.order_time=2.0
params.topping_time=2.0
params.loading_time=0.5
params.baking_time=10.0
params.pack_time=1.0
params.inter_arrival=10.0
params.nr_repetitions=100
class Pizza(Process):
def __init__(self):
Process.__init__(self,name="Pizza order %s"%Pizza.id)
Pizza.id+=1
def pizza_process(self):
yield hold,self,params.order_time
yield hold,self,params.topping_time
yield hold,self,params.loading_time
yield hold,self,params.baking_time
yield hold,self,params.pack_time
Pizza.done+=1
class CustomerArrivals(Process):
def __init__(self):
Process.__init__(self,name="CustomerArrivals")
def arrival_process(self):
while True:
if 120<=now()<=180: #lunch rush hour
inter_time=params.inter_arrival/4.0
else:
inter_time=params.inter_arrival
arrival_delay=random.expovariate(1.0/inter_time)
if now()+arrival_delay<=params.open_time:
yield hold,self,arrival_delay
p=Pizza()
activate(p,p.pizza_process())
else:
break
def model():
obs_pizzas=0
print '\nSimulation parameters:\n',params.__dict__
for reps in range(params.nr_repetitions):
Pizza.id=0
Pizza.done=0
initialize()
c=CustomerArrivals()
activate(c,c.arrival_process())
simulate(until=params.open_time+30)
obs_pizzas+=Pizza.done
print "Results:\nMean number of pizzas per day: %s (%s runs)"\
%(obs_pizzas/params.nr_repetitions,params.nr_repetitions)
model()
params.inter_arrival=5.0
model()
params.inter_arrival=15.0
model()
He runs the revised program and gets:
>>>python PizzaPy1.py
Simulation parameters:
{'open_time': 480.0, 'order_time': 2.0, 'topping_time': 2.0,
'pack_time': 1.0, 'inter_arrival': 10.0, 'loading_time': 0.5,
'baking_time': 10.0, 'nr_repetitions': 100}
Results:
Mean number of pizzas per day: 62 (100 runs)
Simulation parameters:
{'open_time': 480.0, 'order_time': 2.0, 'topping_time': 2.0,
'pack_time': 1.0, 'inter_arrival': 5.0, 'loading_time': 0.5,
'baking_time': 10.0, 'nr_repetitions': 100}
Results:
Mean number of pizzas per day: 127 (100 runs)
Simulation parameters:
{'open_time': 480.0, 'order_time': 2.0, 'topping_time': 2.0,
'pack_time': 1.0, 'inter_arrival': 15.0, 'loading_time': 0.5,
'baking_time': 10.0, 'nr_repetitions': 100}
Results:
Mean number of pizzas per day: 41 (100 runs)
Now Joe knows that for the mean interarrival delay range 5 ... 15 minutes, he can expect at most about 127 ... 41 pizza sales per business day, with some small statistical error. This assumes that he does not lose any customer due to excessive waiting.
Simulating shared resources staff and oven
Joe knows, however, that the resources required by all Pizza process instances, namely the staff and the oven, may become overwhelmed by the customer stream. Insufficient staff or an oven that's too small can lead to customer queues building up. If the waiting times get too long, new customers entering the shop may leave without ordering. On the other hand, if he has too much staff or buys too big an oven, his results are sub-optimal; he can even lose money. Joe suspects that the staff is the bottleneck, rather than the oven: one pizza takes 5.5 minutes of staff time and 11 minutes of oven time (loading plus baking). Thus, one person can produce a maximum of 480/5.5 = 87 pizzas per day. The smallest professional pizza ovens have space for at least 10 pizzas. Such ovens can bake (480/11)*10 = 436 pizzas per day. He therefore decides to model these resources, too, and to experiment with their influence on the number of sales and lost business. Joe expands the model to include shared resources and, while he is at it, data collection and plotting of the lost sales over time for a random day (see the file PizzaPy2.py).
1 #PizzaPy2.py
2 from SimPy.Simulation import *
3 import random
4
5 class Params:pass
6 params=Params()
7 #defaults (all times in minutes)
8 params.open_time=480.0
9 params.order_time=2.0
10 params.topping_time=2.0
11 params.loading_time=0.5
12 params.baking_time=10.0
13 params.pack_time=1.0
14 params.inter_arrival=10.0
15 params.nr_repetitions=100
16 params.oven_size=10
17 params.employee_num=1
18 params.tolerance=4
19
20 class Pizza(Process):
21 def __init__(self):
22 Process.__init__(self,name="Pizza %s"%Pizza.id)
23 Pizza.id+=1
24 def pizza_process(self):
25 yield request,self,staff,1
26 yield hold,self,params.order_time
27 yield hold,self,params.topping_time
28 yield request,self,oven
29 yield hold,self,params.loading_time
30 yield release,self,staff
31 yield hold,self,params.baking_time
32 yield request,self,staff,10
33 yield release,self,oven
34 yield hold,self,params.pack_time
35 yield release,self,staff
36 Pizza.done+=1
37 Pizza.waiting_customers-=1
38
39 class CustomerArrivals(Process):
40 def __init__(self):
41 Process.__init__(self,name="CustomerArrivals")
42 def arrival_process(self):
43 while True:
44 if 120<=now()<=180: #lunch rush hour
45 inter_time=params.inter_arrival/4.0
46 else:
47 inter_time=params.inter_arrival
48 arrival_delay=random.expovariate(1.0/inter_time)
49 if now()+arrival_delay<=params.open_time:
50 yield hold,self,arrival_delay
51 if Pizza.waiting_customers<=params.tolerance:
52 p=Pizza()
53 Pizza.waiting_customers+=1
54 activate(p,p.pizza_process())
55 else: #lost sale
56 Pizza.lost+=1
57 lost_sales.observe(Pizza.lost,t=now())
58 else:
59 break
60 lost_sales.observe(Pizza.lost,t=now())
61
62 def plot_a_day():
63 from SimPy.SimPlot import *
64 plt=SimPlot()
65 plt.plotStep(lost_sales,color='red',title="Lost sales on a "+\
66 "random day "+ "(%s employees, %s min mean interarrival time)"\
67 %(params.employee_num,params.inter_arrival))
68 plt.mainloop()
69
70 def model():
71 global oven,staff,lost_sales
72 obs_pizzas=0
73 obs_lost=0
74 print '\nSimulation parameters:\n',params.__dict__
75 for reps in range(params.nr_repetitions):
76 oven=Resource(name="Oven", capacity=params.oven_size)
77 staff=Resource(name="Staff",unitName="employee",
78 capacity=params.employee_num, qType=PriorityQ)
79 lost_sales=Monitor(name="lost sales",
80 tlab="minutes since opening",ylab="nr lost sales")
81 Pizza.id=0
82 Pizza.lost=0
83 Pizza.done=0
84 Pizza.waiting_customers=0
85 initialize()
86 c=CustomerArrivals()
87 activate(c,c.arrival_process())
88 simulate(until=params.open_time+30)
89 obs_pizzas+=Pizza.done
90 obs_lost+=Pizza.lost
91 print "Results:\nMean number of pizzas sold per day: %s (%s runs)"\
92 %(obs_pizzas/params.nr_repetitions,params.nr_repetitions)
93 print "Observed mean number of sales lost per day: %s"\
94 %(obs_lost/params.nr_repetitions)
95 plot_a_day()
96
97 params.inter_arrival=5.0
98 model()
99 params.employee_num=3
100 model()
The major changes in PizzaPy2.py can be found in the pizza_process() method (lines 25-37) and in the mode in lines 76-95. Line 25 requests a staff member (with normal priority, 1) to take and prepare a customer's order. Line 28 requests a place in the oven. On line 30, the staff member is released for other work after loading the pizza, while line 32 requests a staff member to unload the pizza and complete the transaction. This request is done with high priority (10) to prevent pizzas from burning in the oven. Lines 76-77 introduces oven and staff as shared resources with capacity place and employee, respectively. For the oven, resource requests are handled first in/first out (FIFO). For staff, a priority regime allows certain staff tasks to be given higher priority than others. Line 79 adds lost_sales for tracking the daily time series of lost sales. On line 95, the call to plot_a_day() provides insight into when lost sales occur. plot_a_day( ) calls a function to plot lost_sales data collected for the last simulated day. When the number of customer waiting exceeds a tolerance threshold, the customer leaves without ordering. SimPySimPy has a simple plotting library (based on Tk) that can easily plot Monitor instances. This is used to plot the lost_sales time series as a stepped line, shown in lines 62-68. When Joe runs this version of the model, he gets this output:
>>> python PizzaPy2.py
Simulation parameters:
{'open_time': 480.0, 'order_time': 2.0, 'tolerance': 4,
'topping_time': 2.0, 'pack_time': 1.0, 'employee_num': 1,
'oven_size': 10, 'inter_arrival': 5.0, 'loading_time': 0.5,
'baking_time': 10.0, 'nr_repetitions': 100}
Results:
Mean number of pizzas sold per day: 78 (100 runs)
Mean number of sales lost per day: 52
Simulation parameters:
{'open_time': 480.0, 'order_time': 2.0, 'tolerance': 4,
'topping_time': 2.0, 'pack_time': 1.0, 'employee_num': 3,
'oven_size': 10, 'inter_arrival': 5.0, 'loading_time': 0.5,
'baking_time': 10.0, 'nr_repetitions': 100}
Results:
Mean number of pizzas sold per day: 91 (100 runs)
Mean number of sales lost per day: 38
Joe sees that for the high arrival rate, his pizza shop would lose a high percentage of sales, even with a staff of three. The run also generates two plots, shown in Figure Two and Figure Three, showing when during a day the lost sales occur. They can be saved as Postscript files.
Figure 2: Plot of lost sales on a random day (1 employee)
Figure 3: Plot of lost sales on a random day (3 employees) As expected, Joe sees from the slopes of the curves that he loses most sales during the lunch rush hour. Well, now Joe knows enough: it's going to be difficult to make money with the PizzaPy shop. At least he found out before investing any money! So, Joe is back on the market, looking for his next Python software project. But before that, he quickly puts together a user-friendly GUI version of his PizzaPy simulator, with help, interactive parameter changes, and lots more (PizzaPyGUI.py). Who knows, maybe somebody else wants to look at how to set up a pizza takeaway!
Simply SimPy
SimPy is a free, fully-fledged simulation development package that competes with many (expensive) commercial simulation packages. Because its written in Python, SimPy is a great toolbox for developing new simulation systems or model libraries, via new yield statements (to provide new inter-process scheduling constructs) or via sub-classing and inheritance. The SimPy web site shows how easy that is. SimPy's homepage is http://simpy.sourceforge.net where documentation, papers, tutorials etc. can be found. The current production version is 1.5. There is also a fledgling SimPy wiki at http://www.mcs.vuw.ac.nz/cgi-bin/wiki/SimPy. The public mailing list for SimPy is simpy-users@lists.sourceforge.net. Dr. Klaus G. Müller |
|
|||||||||||||||||||||||||||||||||||||||||||||
|
|
|
|
|
|
|
||||||||||||||||||||||||||||||||||||||||||
|
||||||||||||||||||||||||||||||||||||||||||||||||
|
Py is committed to bringing you great Python Articles. | |
||||||||||||||||||||||||||||||||||||||||||||||
|
||||||||||||||||||||||||||||||||||||||||||||||||
 |
Reproduction of material from any of PyZine's pages without prior written permission is strictly prohibited. Copyright 2003 - 2005 PyZine |
|