bcc.py

""" 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], )