16 SimPy and Simulation Collection Example Comparisons

This chapter has 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 racket/base

;  Cellphone simulation from PySim

(require (planet williams/simulation/simulation))

; 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)

(for ((i (in-range n)))

(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)

(for ((i (in-range Nhours)))

; 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 racket/base

; Jackson Network from PySim

(require (planet williams/simulation/simulation-with-graphics))

; Parameters

(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)

; Global Variables (set in main)

(define n-arrivals 0)

(define n-in-system #f)

(define time-in-system #f)

(define processor #f)

; Data Collection

(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))   ; start node is always 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)

(for ((i (in-range max-messages)))

(set! n-arrivals (+ n-arrivals 1))

(set! n-arrivals (+ n-arrivals 1))

(schedule #:now (message i))

(wait (random-exponential (/ mu)))))

(define (main)

(with-new-simulation-environment

; Initialize global variables

(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