From 118wiki

Jump to: navigation, search



One of the most important applications of computer science is Computer simulation. It's an enormous area used all over in the sciences, government, economics, finance, numerous areas of commerce, and yes, even in computer networks (for example, to analyze queueing networks or to predict how new protocols will behave).

The two main branches of simulation are continuous models and discrete models.

Warmup: Probability Question

Suppose we paint green two sides of a dice. What is the probability, if you toss the dice sixteen times, that you will roll green two consecutive times? Of course we can work this out with probability theory, but another way to get some idea of the answer (really just an estimate) is to simulate rolling such a dice. Here's a code snippet to do this:

def roll():
 import random  # use Python's random.choice function to simulate a roll
 result = "notGreen"
 for c in range(16):  # try 16 times
   previous = result # save result of previous roll
   result = random.choice( (["Green"]*2) + (["notGreen"]*4) )
   if result == previous == "Green":  return True     
 return False

The code above simulates rolling the dice sixteen times, returning True if two consecutive rolls were green. Clearly, simulating just one execution of roll() cannot be a reliable estimate of the probability -- we need to sample over many experiments. Here, for instance, is a sampling of 10,000 experiments:

def sample(): 
 S = 0  # let S be the sum of "successful" experiments
 for i in range(10000):  
    if roll():  S += 1
 return float(S)/10000.0

The sampling above will return the fraction of the 10,000 samples in which green turned up twice consecutively in sixteen tosses -- this is an estimate of the probability in question.

This type of "probability" simulation can sometimes be useful to get an estimate, which confirms any mathematical derivation of the probability. (Not all probability question are so easy as this one, see this, for example.)

Discrete Event Simulation

Discrete event simulation is appropriate for many cases of simulating computer networks, games, operating systems, business applications where the items are discrete objects, and so on. Many software packages exist for discrete event simulation (some are even in Python, such as SimPy). These packages can be complex, to support many options and to optimize performance, we we'll start with a rather simple one.

The main components of a discrete event simulator are explained below.

Simulation Time

Simulated "time" can run faster than real time or slower than real time. Simulators commonly keep track of simulation time using a "simulation clock" (which often turns out just to be some variable or counter). The rate at which simulation time advances need not be regular. It can increase by very small amounts as thousands of events are simulated, then jump by a large amount, as needed to efficiently run the simulation.


A "simulated event" is an object used by the simulator. Events contain timestamps (which say when, in terms of simulation time, than an event should be simulated) and other data describing application-dependent details about the event being simulated.

Event Queue

During a simulation run, there is typically some queue or list containing events to be simulated. When an event is simulated, or fired, it is removed from the queue.


The heart of the simulator is a scheduling algorithm. It consists of a loop to advance the simulation clock to the time of the first event on the event queue, remove that event from the queue, and perform some calculations appropriate for that event: in other words, the loop continually fires events.

Event Firing

What exactly happens when an event fires depends on what is being simulated. If the simulation models motion of a battery-powered robot, then an event could drain some small amount of power from a battery --- this would likely be simulated by subtracting an amount from a "battery" variable. For simulating applications related to queueing theory, an event could represent the arrival of a customer. In this case, a convenient programming technique is to have each such arrival event schedule the next customer arrival event. That is, one event can create another event and put it on the queue to be simulated in the "future" (with respect to simulated time). More generally, a simulated event can create many other events.

Example: Simulated Air Traffic

First, let's make a formal model of air traffic. The "objects" to be simulated are airplanes, airport locations, and perhaps passengers and freight. We make data structures to represent all these objects; most likely, inventing a class for airplane, a class for airport, etc, and then instantiating that class for each object will result in concise code.

Events for this example would be takeoffs, landing, loading, taxi-ing, and so on -- depending on the level of detail we want to simulate.

Simulation time could start at zero, or at some more meaningfull value (to represent 5am, for instance). Before launching the simulation, we'll want to put some initial events in the queue.

Example: Simulation of an M/M/1 Queue

Notation M/M/1 means that we have one server, a FIFO Queue for that server, the service time is exponentially distributed, and the interarrival rate (the "gap" time between customers) is also exponential. This implies that the customer arrival rate is Poisson distributed.

In other words, this is just what we have studied already in Queueing Theory.

Below is a Python program to simulate the queue. (You can also download or view it: media:Shoqave.txt.) Below, I've broken the program into parts to add some commentary between the parts.

import math, random

# Define global variable for current simulation time

simTime = 0.0
simLimit = 10000.0   # this will time for ending the simulation

# Define global variable for dictionary of objects being simulated
# each item in simObjects is of the form  name:reference  where 'name'
# is the name of a simulated object and 'reference' is the object's
# handle (a reference to the object)

simObjects = { }

The simObjects dictionary may seem strange, here is how it is set up below: simObjects = { 'S0':S0 , 'Q0':Q0 , 'V':V }. This is really just a way to simplify some of the coding later. Notice the comment above says that this is a "global" variable; that means global in the sense of Python's global statement, which allows functions and methods to refer to variables outside of the local scope (just as can be done in C, without needing this kind of declaration).

# Define global variable for simulation event list
#    items in this list are tuples of the form
#            (simTime,object,parameter)  
#    where object is a simulated object to be fired
#    when the simulation time reaches simTime.  The
#    firing of the object will be done by executing
#    the method 
eventList = [ ] 

The items in eventList will always be sorted in increasing order with respect to simTime (the first element of the tuple). This will make scheduling easier.

# Define schedule() function that adds an item to
# the simulation event list.  It maintains the invariant
# that the event list is sorted by simulation time

def schedule(t,s,p): 
   global eventList

   # special case for empty list
   if len(eventList) == 0:  
      eventList = [ (t,s,p) ]

   # partition eventList into two parts:  list 
   # with times smaller than t and list with times
   # greater than or equal to t

   smaller, greaterEqual, i = ( [], [], 0 )
   while i < len(eventList):  
      a,b,c = eventList[i]
      if a < t:  
	 i += 1
         greaterEqual = eventList[i:]
   # new eventList just puts (t,s,p) in the middle

   eventList = smaller + [ (t,s,p) ] + greaterEqual

You may wonder why I didn't just Python's built-in "sort" method; the answer is that tuples are being sorted, and there isn't a built-in way to compare tuples. So the code above does the sorting by putting a new event into an already-sorted list at just the right place.

# Define customer server class

class customerServe:

  def __init__(self,v,qsize,name):

    # establish mean service time (v) and initial
    # queue of customers (qsize), with "name" as 
    # the object name for this virtual server

    self.meanServiceTime = v
    self.queue = [ ]

    # for each customer in the initial queue
    # schedule a "done" time for that customer
    for i in range(qsize):  self.enqueue(name)  
  def enqueue(self,name):

    global simTime 

    # a new customer has arrived; so we need 
    # to add this new customer to the current 
    # queue and and also generate a future 
    # doneTime for the new customer

    # first, because we use first-come-first-serve
    # queueing discipline, find end of current queue
    # and its scheduled doneTime

    if len(self.queue) > 0:
       t = self.queue[-1]
    else: t = simTime

    # generate done time by a random choice
    # from an exponential distribution with
    # mean from meanServiceTime
    r = random.random()
    r = - math.log(r)
    r *= self.meanServiceTime
    doneTime = t + r

    # add done time to list in this queue
    # schedule simulator to fire at doneTime

  def fire(self,parameter):

    # this function is invoked by the simulator
    # when a scheduled event (ie, customer finishes)
    # takes place;  return True, otherwise the 
    # entire simulation will stop

    # remove first customer from current queue
    self.queue = self.queue[1:]
    return True

Above, the "server" class is a generic class for servers. We can create as many servers as desired by creating objects from this class. The constructor has parameters for the service rate (μ), initial queue, and symbolic name for this server. The enqueue method precalculates how long the service will be at the instant a customer is added to the queue. Notice how this service time is generated:

   r = random.random()
   r = - math.log(r)
   r *= self.meanServiceTime

That mysterious code does random service time generation that simulates an exponentially distributed service time (with μ = self.meanServiceTime). Take my word for it, or study this if you really want to understand the theory behind it.

The fire method of the server is very simple, it just removes the first item of the queue. If this were a more detailed simulation, say of aircraft landing, then the fire method would probably do more (it might enqueue a new event for another server, simulating the handing over from traffic control to ground control of an aircraft).

# Define customer generator class

class customerGenerate:

  def __init__(self,v,servers,name):

    # establish mean interarrival time (v) and a  
    # list of candidates (servers) to handle the generated
    # customers -- if customers have a choice ---
    # that will handle the generated customers,
    # with "name" as the object name for this generator 

    self.meanInterArrivalTime = v
    self.servers = servers = name

    # generate one initial customer 
  def generate(self,name):

    global simTime 

    # create a new customer sometime in future 

    # generate arrival time by a random choice
    # from an exponential distribution with
    # mean from meanInterArrivalTime
    r = random.random()
    r = - math.log(r)
    r *= self.meanInterArrivalTime
    arriveTime = simTime + r

    # schedule simulator to fire at arriveTime

  def fire(self,parameter):

    global simObjects 

    # this function is invoked by the simulator
    # when a customer "arrives" (as we previously scheduled)

    # currently, it just puts the customer on one server queue
    server = self.servers[0]
    servRef = simObjects[server]

    # finally, schedule another customer to arrive
    return True

Above is the "customer" class, from which we can instantiate any number of objects to represent streams of customers. That is, each object will represent a stream of customers that feed into the same queue. The constructor enables tailoring the stream according to the mean interarrival rate (which is 1/λ). The constructor also has a parameter for a list of potential servers for this customer. The generate method simulates the creation of one new customer by adding a random "gap" time to the current simulation clock; it invokes the schedule function to add a new event representing the customer's arrival to the queue.

The fire method is invoked at the simulated time when the customer arrives, and this customer should queue up on a server. The fire method chooses one server and enqueues this customer on that server. Finally, the fire method also schedules a future time at which the next customer will arrive.

# Define visualize class
#   This class only exists to periodically examine
#   the state of the simulation.  In this case, our
#   only interest is the length of the queue

class visualize:

  def __init__(self,name,limit):

    global simTime = name
    self.limit = limit
    self.history = [ ]

    # generate first event, which will reschedule
    # itself thereafter 

  def fire(self,parameter):

    global simTime, simObjects 

    # record queue size at current simTime
    qsize = len( simObjects['S0'].queue )


    if simTime >= self.limit: return False
    else: return True

Many simulators have some built-in way to save the result of events as they are being simulated, either for graphing (visual display) or for saving some statistics. Here, instead, we just define another kind of event called "the recording event", which regularly schedules itself, also building a list of the queue sizes observed at each firing. In a way, this is cheating: the visualize event is not actually an event in what we are simulating, but the pattern of simulated event scheduling nicely fits with what we want to do.

# Define simCycle() function that runs next event
# of the simulation.  It returns False iff the 
# simulation is over, that is, it returns True if
# the simulation should be continued

def simCycle(): 
   global eventList, simTime, simObjects

   if len(eventList) == 0:  return False 

   # get next event, which defines the next simulation
   # time, and fire the corresponding action

   simTime, object, parameter = eventList[0]
   eventList = eventList[1:] 

   # print "> ", object, " fired @ ", simTime 

   result = simObjects[object].fire(parameter) 
   return result 

# Main simulator logic:  do simCycle forever
def simulator():
   go = True
   while go:  go = simCycle()

The scheduler, shown above, turns out to be really simple. The code that follows creates initial server and customer objects, then invokes the scheduler; when the scheduler finishes, it prints the observed mean queue length.

#----------- Example: simulate server ---------------------

S0 = customerServe(2,0,"S0")
Q0 = customerGenerate(2,["S0"],"Q0")
V  = visualize("V",simLimit)
simObjects = { 'S0':S0 , 'Q0':Q0 , 'V':V }
random.seed(1)  # force simulation to be repeatable
print "Running simulation up to time ", simLimit 

# after simulation is done,
# compute average queue length
qs = 0.0 
for i in range(len(V.history)): qs += V.history[i]
qave = qs / float(len(V.history))
print "Average queue length was ", qave

Distributed and Parallel Simulation

The source for notes and lectures on this topic are mainly from a book (ISBN 0471183830) by Richard Fujimoto; also I use Fujimoto's introductory presentation slides (media:Fujimoto.ppt, media:Fujimoto.pdf).

The basic idea of distributed simulation is to partition the work of discrete event simulation, as described above, onto a set of hosts that communicate by sending messages. To be more concrete, each host "owns" some subset of the simulated objects. Each host also has its own clock, its own event list, and its own scheduler.

Consider two hosts A and B participating in the distributed simulation. If none of the objects on host A have anything to do with objects on host B, then it could be that simulations at A and B can proceed without any messages -- we could even see the distributed simulation run twice as fast as putting all the work on one host. But usually, objects and A and B are related. For instance, an object fired on host A could schedule an event for some object on host B.

Therefore, scheduling an event is an operation that sends a message, containing the time, object reference, and parameters, to the host that owns the object.

One of the major differences between traditional single-CPU simulation and distributed simulation is that each host has its own clock. This can potentially allow one host to virtually run ahead (in time) of other hosts. Does this make sense? Not always. We get into trouble in the following situation: let host A's clock be at time 100 and host B's clock at time 200; then suppose host A fires an event that schedules a future event for time 150 (this is in A's future) but for an object of B. Now A sends a message for the event to B, but when B receives the event message for time 150, the clock is already at time 200, so the message is too late: if B fires the event, the simulation would be inconsistent because events would not fire in proper time sequence.

During the lecture, I explain how a conservative strategy of event firing will avoid inconsistent firing of events. This conservative method can result in simulation deadlock; one way to avoid deadlock is to have each host periodically send null messages.

The alternative to conservative event processing is an optimistic strategy for firing events. In the optimistic strategy, each simulation host is allowed to fire events, even without certain knowledge of other host simulation clocks. Because this could result in inconsistency, each simultion host is prepared to "undo" or rollback (data management) its part of the simulation to an earlier state, so that events can be reprocessed in the correct order. This idea, sometimes called time warp, results in some complications, including cascaded rollbacks and anti-messages. During the lecture I show examples of messages, anti-messages, and rollbacks.

Continuous Simulation

The alternative to discrete event simulation could be called "continuous" simulation, because it approximates the treatment of time as a continuous variable. In this type of simulation, clocks do not jump to the next event: in fact, there are no events. Rather, the behavior of simulated processes are described using differential equations; instead of discrete objects, there are quantities. Such simulations are useful to model weather, ocean currents, even vehicle traffic (my modeling the traffic as a flow, like a liquid).

An good way to introduce the topic is to illustrate the famous Lotka-Volterra Equations.

The Lotka-Volterra equations, described here, give an explanation of how predators and prey interact in a biological system. The equations are differential equations governing the rate of population growth (or decline) in two species. While the equilibrium solution for these equations can be found directly, simulation is more versatile because we can modify the equations and easily see the result by simulation. This page has more details about the equations and suggests exploring the results by simulation. Rather than using a simulation package, you can simulate the populations of foxes and rabbits using a Python program that I wrote: media:Volterra.txt. I'll illustrate this in class and show how changing the parameters gives different results.

A somewhat lengthy paper shows how ideas from the Lotka-Volterra equations can be generalized to very complex systems (glance at the diagram on pages 14-15, and you get the idea). In class we'll discuss how one might make the equations more realistic an explore the consequences using simulation.

Sim World, Simulations of Networks?

It's not hard to imagine how a simulation based on equations could be used for guiding the way virtual realities and related games (Sims World) could be implemented. Simulation is also very useful for testing hypotheses about systems before they are built: we would like to know if a system will have unexpected behaviors before spending millions to implement it! Economists, businesses, chip-manufacturers, insurance companies -- are but a few examples of organizations where simulation has become a crucial tool. In the networking arena, simulators such as NS are very widely used to evaluate the performance of network protocols. Generally, most networking analysis uses discrete event simulation, however some large-scale aspects can also benefit from continuous simulations.

Personal tools