""" bcc.py
Queue with blocked customers cleared
Jobs (e.g messages) arrive randomly at rate 1.0 per minute at a
2-server system. Mean service time is 0.75 minutes and the
service-time distribution is (1) exponential, (2) Erlang-5, or (3)
hyperexponential with p=1/8,m1=2.0, and m2=4/7. However no
queue is allowed; a job arriving when all the servers are busy is
rejected.
Develop and run a simulation program to estimate the probability of
rejection (which, in steady-state, is the same as p(c)) Measure
and compare the probability for each service time distribution.
Though you should test the program with a trace, running just a few
jobs, the final runs should be of 10000 jobs without a trace. Stop
the simulation when 10000 jobs have been generated.
"""
from SimPy.Simulation import *
from random import seed,Random,expovariate,uniform
## Model components ------------------------
dist = ""
def bcc(lam, mu,s):
""" bcc - blocked customers cleared model
- returns p[i], i = 0,1,..s.
- ps = p[s] = prob of blocking
- lameff = effective arrival rate = lam*(1-ps)
See Winston 22.11 for Blocked Customers Cleared Model (Erlang B formula)
"""
rho = lam/mu
n = range(s+1)
p = [0]*(s+1)
p[0] = 1
sump = 1.0
for i in n[1:]:
p[i] = (rho/i)*p[i-1]
sump = sump + p[i]
p0 = 1.0/sump
for i in n:
p[i] = p[i]*p0
p0 = p[0]
ps = p[s]
lameff = lam*(1-ps)
L = rho*(1-ps)
return {'lambda':lam,'mu':mu,'s':s,
'p0':p0,'p[i]':p,'ps':ps, 'L':L}
def ErlangVariate(mean,K):
""" Erlang random variate
mean = mean
K = shape parameter
g = rv to be used
"""
sum = 0.0 ; mu = K/mean
for i in range(K):
sum += expovariate(mu)
return (sum)
def HyperVariate(p,m1,m2):
""" Hyperexponential random variate
p = prob of branch 1
m1 = mean of exponential, branch 1
m2 = mean of exponential, branch 2
g = rv to be used
"""
if random() < p:
return expovariate(1.0/m1)
else: return expovariate(1.0/m2)
def testHyperVariate():
""" tests the HyerVariate rv generator"""
ERR=0
x = (1.0981,1.45546,5.7470156)
p = 0.0, 1.0 ,0.5
g = Random(1113355)
for i in range(3):
x1 = HyperVariate(p[i],1.0,10.0,g)
#print p[i], x1
assert abs(x1 - x[i]) < 0.001,'HyperVariate error'
def erlangB(rho,c):
""" Erlang's B formula for probabilities in no-queue
Returns p[n] list
see also SPlus and R version in que.q mmcK
que.py has bcc.
"""
n = range(c+1) ; pn = range(c+1)
term = 1
pn[0] = 1
sum = 1 ; term = 1.0
i=1
while i < (c+1):
term *= rho/i
pn[i] = term
sum += pn[i]
i += 1
for i in n: pn[i] = pn[i]/sum
return(pn)
class JobGen(Process):
""" generates a sequence of Jobs
"""
def execute(self,JobRate,MaxJob,mu):
global NoInService, Busy
for i in range(MaxJob):
j = Job()
activate(j,j.execute(i,mu),delay=0.0)
t = expovariate(JobRate)
MT.tally(t)
yield hold,self,t
self.trace("Job generator finished")
def trace(self,message):
if JobGenTRACING: print "%8.4f \t%s"%(now(), message)
class Job(Process):
""" Jobs that are either accepted or rejected
"""
def execute(self,i,mu):
""" Job execution, only if accepted"""
global NoInService,Busy,dist,NoRejected
if NoInService < c:
self.trace("Job %2d accepted b=%1d"%(i,Busy))
NoInService +=1
if NoInService == c:
Busy =1
try: BM.accum(Busy,now())
except: "accum error BM=",BM
#yield hold,self,Job.g.expovariate(self.mu); dist= "Exponential"
yield hold,self,ErlangVariate(1.0/mu,5); dist= "Erlang "
#yield hold,self,HyperVariate(1.0/8,m1=2.0,m2=4.0/7,g=Job.g); dist= "HyperExpon "
NoInService -=1
Busy =0
BM.accum(Busy,now())
self.trace("Job %2d leaving b=%1d"%(i,Busy))
else:
self.trace("Job %2d REJECT b=%1d"%(i,Busy))
NoRejected +=1
def trace(self,message):
if JobTRACING: print "%8.4f \t%s"%(now(), message)
## Experiment data -------------------------
c = 2
lam = 1.0 ## per minute
mu = 1.0/0.75 ## per minute
p = 1.0/8 ; m1= 2.0; m2 = 4.0/7.0
K = 5
rho = lam/mu
NoRejected = 0
NoInService = 0
Busy = 0
JobRate = lam
JobMax = 10000
JobTRACING = 0
JobGenTRACING = 0
## Model/Experiment ------------------------------
seed(111333)
BM=Monitor()
MT = Monitor()
initialize()
jbg = JobGen()
activate(jbg,jbg.execute(1.0,JobMax,mu),0.0)
simulate(until=20000.0)
## Analysis/output -------------------------
print 'bcc'
print "time at the end =",now()
print "now=",now(), " startTime ",BM.startTime
print "No Rejected = %d, ratio= %s"%(NoRejected,(1.0*NoRejected)/JobMax)
print "Busy proportion (%6s) = %8.6f"%(dist,BM.timeAverage(now()),)
print "Erlang pc (th) = %8.6f"%(erlangB(rho,c)[c],)
|