On this page:
16.1 Cellphone
16.1.1 Simpy Example – Cellphone
16.1.2 Simulation Collection Example – Cellphone
16.2 Jackson
16.2.1 Simpy Example – Jackson
16.2.2 Simulation Collection Example – Jackson
Version: 4.1.3

16 SimPy and Simulation Collection Example Comparisons

    16.1 Cellphone

      16.1.1 Simpy Example – Cellphone

      16.1.2 Simulation Collection Example – Cellphone

    16.2 Jackson

      16.2.1 Simpy Example – Jackson

      16.2.2 Simulation Collection Example – Jackson

This chapter had some example from the SimPy web site and equivalent examples in the simulation collection.

16.1 Cellphone

16.1.1 Simpy Example – Cellphone

#!/usr/bin/env python
""" Simulate the operation of a BCC cellphone system using SimPy
"""
from __future__ import generators   # not needed for Python 2.3+
from SimPy.Simulation import *
from random import Random,expovariate
class Generator(Process):
 """ generates a sequence of calls """
 def __init__(self, maxN, lam):
     global busyEndTime
     Process.__init__(self)
     self.name = "generator"
     self.maxN = maxN
     self.lam = lam
     self.iatime = 1.0/lam
     self.rv = Random(gSeed)
     busyEndTime = now() # simulation start time
 def execute(self):
     for i in range(self.maxN):
         j = Job(i)
         activate(j,j.execute())
         yield hold,self,self.rv.expovariate(lam)
     self.trace("WARNING generator finished -- ran out of jobs")
 def trace(self,message):
     if GTRACING: print "%7.4f \t%s"%(self.time(), message)
class Job(Process):
    """ instances of the Job class represent calls arriving at
    random at the chnnels.
    """
    def __init__(self,i):
        Process.__init__(self)
        self.i = i
        self.name = "Job"+`i`
    def execute(self):
        global busyStartTime,totalBusyVisits,totalBusyTime
        global Nfree,busyEndTime,Jrv
        self.trace("arrived ")
        if Nfree == 0: self.trace("blocked and left")
        else:
             self.trace("got a channel")
             Nfree -=  1
             if Nfree == 0:
                 self.trace("start busy period======")
                 busyStartTime = now()
                 totalBusyVisits += 1
                 interIdleTime = now() - busyEndTime
             yield hold,self,Jrv.expovariate(mu)
             self.trace("finished")
             if Nfree == 0:
                 self.trace("end   busy period++++++")
                 busyEndTime = now()
                 busy = now() - busyStartTime
                 self.trace("         busy  = %9.4f"%(busy,))
                 totalBusyTime +=busy
             Nfree += 1
    def trace(self,message):
         if TRACING: print "%7.4f \t%s %s "%(now(), message , self.name)
class Statistician(Process):
     """ observes the system at intervals """
     def __init__(self,Nhours,interv,gap):
         Process.__init__(self)
         self.Nhours = Nhours
         self.interv = interv
         self.gap = gap
         self.name="Statistician"
     def execute(self):
         global busyStartTime,totalBusyTime,totalBusyVisits
         global hourBusyTime,hourBusyVisits
         for i in range(self.Nhours):
             yield hold,self,self.gap
             totalBusyTime = 0.0
             totalBusyVisits = 0
             if Nfree == 0: busyStartTime = now()
             yield hold,self,self.interv
             if Nfree == 0: totalBusyTime += now()-busyStartTime
             if STRACING:  print "%7.3f %5d"%(totalBusyTime,totalBusyVisits)
             m.tally(totalBusyTime)
             bn.tally(totalBusyVisits)
         print("Busy Time:   mean = %10.5f var= %10.5f"%(m.mean(),m.var()))
         print("Busy Number: mean = %10.5f var= %10.5f"%(bn.mean(),bn.var()))
     def trace(self,message):
         if STRACING: print "%7.4f \t%s"%(self.time(), message)
totalBusyVisits = 0
totalBusyTime   = 0.0
NChannels =  4  # number of channels in the cell
Nfree   = NChannels
maxN    = 1000
gSeed   = 11111111
JrvSeed = 3333333
lam = 1.0
mu = 0.6667
meanLifeTime = 1.0/mu
TRACING  = 0
GTRACING = 0
STRACING = 1
Nhours  =  10
interv = 60.0    # monitor
gap    = 15.0    # monitor
print "lambda    mu      s  Nhours interv  gap"
print "%7.4f %6.4f %4d %4d    %6.2f %6.2f"%(lam,mu,NChannels,Nhours,interv,gap)
m = Monitor()
bn=Monitor()
Jrv = Random(JrvSeed)
s = Statistician(Nhours,interv,gap)
initialize()
g = Generator(maxN, lam)
activate(g,g.execute())
activate(s,s.execute())
simulate(until=10000.0)

The output from the program follows.

lambda    mu      s  Nhours interv  gap
 1.0000 0.6667    4   10     60.00  15.00
  3.274     8
  0.444     2
  3.649     8
  0.652     3
  2.277     6
  1.372     9
  0.411     4
  0.067     2
  2.910     8
  5.073     9
Busy Time:   mean =    2.01297 var=    2.55814
Busy Number: mean =    5.90000 var=    7.49000

16.1.2 Simulation Collection Example – Cellphone

  #lang scheme/base
  ; Cellphone simulation from PySim
  
  (require (planet williams/simulation/simulation))
  (require (planet williams/science/random-distributions))
  
  ; Parameters
  
  (define Nchannels 4)
  (define maxN 1000)
  (define lam 1.0)
  (define mu 0.6667)
  (define meanLifeTime (/ 1.0 mu))
  (define Nhours 10)
  (define interv 60.0)
  (define gap 15.0)
  
  ; Globals
  
  (define Nfree Nchannels)
  (define totalBusyVisits 0)
  (define totalBusyTime 0.0)
  (define busyStartTime 0.0)
  (define busyEndTime 0.0)
  
  
  ; Statistics
  
  (define m #f)
  (define bn #f)
  
  (define-process (generator n)
    (do ((i 0 (+ i 1)))
        ((= i n) (void))
      (schedule now (job i))
      (wait (random-exponential (/ 1.0 lam)))))
  
  (define-process (job i)
    (when (> Nfree 0)
      (set! Nfree (- Nfree 1))
      (when (= Nfree 0)
        ; Start busy period
        (set! busyStartTime (current-simulation-time))
        (set! totalBusyVisits (+ totalBusyVisits 1)))
      (work (random-exponential (/ 1.0 mu)))
      (when (= Nfree 0)
        ; End busy period
        (set! busyEndTime (current-simulation-time))
        (set! totalBusyTime (+ totalBusyTime
                               (- (current-simulation-time)
                                  busyStartTime))))
      (set! Nfree (+ Nfree 1))))
  
  (define-process (statistician)
    (do ((i 0 (+ i 1)))
        ((= i Nhours) (void))
      ; wait the specified gap time
      (wait gap)
      ; initialize
      (set! totalBusyTime 0.0)
      (set! totalBusyVisits 0)
      (when (= Nfree 0)
        (set! busyStartTime (current-simulation-time)))
      ; wait the specified interval time
      (wait interv)
      ; trace busy time and busy visits
      (when (= Nfree 0)
        (set! totalBusyTime (+ totalBusyTime
                               (- (current-simulation-time)
                                  busyStartTime))))
      (printf "~a: busy time = ~a; busy visits = ~a~n"
              (current-simulation-time) totalBusyTime totalBusyVisits)
      ; Tally statistics
      (set-variable-value! m totalBusyTime)
      (set-variable-value! bn totalBusyVisits))
    ; Print final statistics
    (printf "Busy time:   mean = ~a, variance = ~a~n"
            (variable-mean m) (variable-variance m))
    (printf "Busy number: mean = ~a, variance = ~a~n"
            (variable-mean bn) (variable-variance bn)))
  
  (define (run-simulation)
    (with-new-simulation-environment
     ; Print parameters
     (printf "lambda = ~a~n" lam)
     (printf "mu     = ~a~n" mu)
     (printf "s      = ~a~n" Nchannels)
     (printf "Nhours = ~a~n" Nhours)
     (printf "interv = ~a~n" interv)
     (printf "gap    = ~a~n" gap)
     ; Initialize global variables
     (set! Nfree Nchannels)
     (set! totalBusyVisits 0)
     (set! totalBusyTime 0.0)
     (set! busyEndTime 0.0)
     ; Create statistics
     (set! m (make-variable))
     (tally (variable-statistics m))
     (set! bn (make-variable))
     (tally (variable-statistics bn))
     ; Create initial processes and start the execution
     (schedule (at 0.0) (generator maxN))
     (schedule (at 0.0) (statistician))
     (schedule (at 10000.0) (stop-simulation))
     (start-simulation)))
  
  (run-simulation)

The output from the program follows.

lambda = 1.0
mu     = 0.6667
s      = 4
Nhours = 10
interv = 60.0
gap    = 15.0
75.0: busy time = 3.19230597325965; busy visits = 11
150.0: busy time = 2.4780498393357675; busy visits = 7
225.0: busy time = 3.5391222097609045; busy visits = 12
300.0: busy time = 2.1031429717027947; busy visits = 6
375.0: busy time = 1.3915494625956057; busy visits = 3
450.0: busy time = 2.9838869629721785; busy visits = 6
525.0: busy time = 3.8242483424334637; busy visits = 8
600.0: busy time = 2.7636093265058435; busy visits = 9
675.0: busy time = 1.73398034270042; busy visits = 8
750.0: busy time = 1.957177682358406; busy visits = 9
Busy time:   mean = 2.596707311362503, variance = 0.5790891717346138
Busy number: mean = 7.9, variance = 6.089999999999996

16.2 Jackson

16.2.1 Simpy Example – Jackson

#!/usr/bin/env python
"""
  Simulation of a jackson queue network
"""
from __future__ import generators  # not needed in Python 2.3+
from SimPy.Simulation import *
from random import Random,expovariate,uniform
def sum(P):
    """ Sum of the real vector P """
    sumP = 0.0
    for i in range(len(P)): sumP += P[i]
    return (sumP)
def choose2dA(i,P,g):
    """  return a random choice from a set j = 0..n-1
       with probs held in list of lists P[j] (n by n)
       using row i
       g = random variable
       call:  next = choose2d(i,P,g)
    """
    U = g.random()
    sumP = 0.0
    for j in range(len(P[i])):  # j = 0..n-1
        sumP +=  P[i][j]
        if U < sumP: break
    return(j)
## -----------------------------------------------
class Msg(Process):
    """a message"""
    def __init__(self,i,node):
        Process.__init__(self)
        self.i = i
        self.node = node
    def execute(self):
        """ executing a message """
        global noInSystem
        startTime = now()
        noInSystem += 1
        ##print "DEBUG noInSystm = ",noInSystem
        NoInSystem.accum(noInSystem)
        self.trace("Arrived node  %d"%(self.node,))
        while self.node <> 3:
            yield request,self,processor[self.node]
            self.trace("Got processor %d"%(self.node,))
            st = Mrv.expovariate(1.0/mean[self.node])
            Stime.observe(st)
            yield hold,self,st
            yield release,self,processor[self.node]
            self.trace("Finished with %d"%(self.node,))
            self.node = choose2dA(self.node,P,Mrv)
            self.trace("Transfer to   %d"%(self.node,))
        TimeInSystem.observe(now()-startTime)
        self.trace(    "leaving       %d"%(self.node,),noInSystem)
        noInSystem -= 1
        NoInSystem.accum(noInSystem)
    def trace(self,message,nn=0):
        if MTRACING: print "%7.4f %3d %10s %3d"%(now(),self.i, message,nn)
## -----------------------------------------------
class Generator(Process):
 """ generates a sequence of msgs """
 def __init__(self, rate,maxT,maxN):
     Process.__init__(self)
     self.name = "Generator"
     self.rate = rate
     self.maxN = maxN
     self.maxT = maxT
     self.g = Random(11335577)
     self.i = 0
 def execute(self):
     while (now() < self.maxT)  & (self.i < self.maxN):
         self.i+=1
         p = Msg(self.i,startNode)
         activate(p,p.execute())
         ## self.trace("starting "+p.name)
         yield hold,self,self.g.expovariate(self.rate)
     self.trace("generator finished with "+`self.i`+" ========")
     self.stopSim()
 def trace(self,message):
     if GTRACING: print "%7.4f \t%s"%(now(), message)
## -----------------------------------------------
# generator:
GTRACING = 0
# messages:
rate = 0.5
noInSystem = 0
MTRACING = 0
startNode = 0
maxNumber = 10000
maxTime = 100000.0
Mrv = Random(77777)
TimeInSystem=Monitor()
NoInSystem= Monitor()
Stime = Monitor()
processor = [Resource(1),Resource(1),Resource(1)]
mean=[1.0, 2.0, 1.0]
P = [[0, 0.5, 0.5, 0],[0,0,0.8,0.2], [0.2,0,0,0.8]]
initialize()
g = Generator(rate,maxTime,maxNumber)
activate(g,g.execute())
simulate(until=100.0)
print "Mean numberm= %8.4f"%(NoInSystem.timeAverage(),)
print "Mean delay  = %8.4f"%(TimeInSystem.mean(),)
print "Mean stime  = %8.4f"%(Stime.mean(),)
print "Total jobs  = %8d"%(TimeInSystem.count(),)
print "Total time  = %8.4f"%(now(),)
print "Mean rate   = %8.4f"%(TimeInSystem.count()/now(),)

The following is the output from the program.

Mean numberm=   1.5412
Mean delay  =   3.9369
Mean stime  =   1.0332
Total jobs  =       38
Total time  = 100.0000
Mean rate   =   0.3800

16.2.2 Simulation Collection Example – Jackson

  #lang scheme/base
  
  
  (require (planet williams/simulation/simulation-with-graphics))
  (require (planet williams/science/random-distributions))
  
  
  
  (define n 3)
  (define m #(1.0 2.0 1.0))
  (define p (vector
             (make-discrete #(0.0 0.5 0.5 0.0))
             (make-discrete #(0.0 0.0 0.8 0.2))
             (make-discrete #(0.2 0.0 0.0 0.8))))
  (define np '(1 1 1))
  (define mu (/ 1.0 1.5))
  (define max-messages 10000)
  
  
  
  (define n-arrivals 0)
  (define n-in-system #f)
  (define time-in-system #f)
  (define processor #f)
  
  
  
  (define-process (message i)
    (let ((arrive-time (current-simulation-time))
          (work-duration 0.0))
      (set-variable-value!
       n-in-system (+ (variable-value n-in-system) 1))
      (let loop ((node 0))
        (when (< node n)
          (with-resource ((vector-ref processor node))
            (work (random-exponential (vector-ref m node))))
          (loop (random-discrete (vector-ref p node)))))
      (set-variable-value!
       n-in-system (- (variable-value n-in-system) 1))
      (set-variable-value!
       time-in-system (- (current-simulation-time)
                         arrive-time))))
  
  (define (generator)
    (do ((i 0 (+ i 1)))
        ((= i max-messages) (void))
      (set! n-arrivals (+ n-arrivals 1))
      (schedule now (message i))
      (wait (random-exponential (/ mu)))))
  
  (define (main)
    (with-new-simulation-environment
  
      (set! n-arrivals 0)
      (set! n-in-system (make-variable 0))
      (accumulate (variable-history n-in-system))
      (set! time-in-system (make-variable))
      (tally (variable-statistics time-in-system))
      (set! processor (list->vector
                       (map make-resource np)))
  
      (schedule (at 0.0) (generator))
      (schedule (at 5000.0) (stop-simulation))
      (start-simulation)
      (printf "Mean number in system = ~a~n" (variable-mean n-in-system))
      (printf "Mean delay in system  = ~a~n" (variable-mean time-in-system))
      (printf "Total time run        = ~a~n" (current-simulation-time))
      (printf "Total jobs arrived    = ~a~n" n-arrivals)
      (printf "Total jobs completed  = ~a~n" (variable-n time-in-system))
      (printf "Average arrival rate  = ~a~n" (/ n-arrivals
                                                (current-simulation-time)))
      (write-special (history-plot (variable-history n-in-system)))
      (newline)))
  
  (main)

The following is the output from the program.

Mean number in system = 14.726063392545544
Mean delay in system  = 21.56635336305308
Total time run        = 5000.0
Total jobs arrived    = 3419
Total jobs completed  = 3413
Average arrival rate  = 0.6838