Simulation of a M/M/1 queue

Confidence interval on average waiting time using batch means

In [1]:
using SimJulia
using Distributions
using RandomStreams
using Distributions

const SEED = 12345

rand_dist(rng::MRG32k3a, Dist::Distribution) = quantile(Dist, rand(rng))

seeds = [SEED, SEED, SEED, SEED, SEED, SEED]
gen = MRG32k3aGen(seeds)
Out[1]:
Seed for next MRG32k3a generator:
[12345,12345,12345,12345,12345,12345]
In [2]:
include("tally.jl")
Out[2]:
stdev (generic function with 1 method)
In [3]:
include("tallystore.jl")
Out[3]:
add (generic function with 2 methods)
In [4]:
println("M/M/1 with processes")

#rates
λ = 1.9
μ = 2.0
ρ = λ/μ

warmup = 500 # this should be set with a formal procedure

type System
    W
    arrival
    service
    counter
    arrgen
    servgen
    
    function System(text::String)
        s = new()
        s.W = Tally(text)
        s.arrival = Exponential(1.0/λ)
        s.service = Exponential(1.0/μ)
        s.arrgen = next_stream(gen)
        s.servgen = next_stream(gen)
        return s;
    end
end

function restart(s::System)
    init(s.W)
    next_substream!(s.arrgen)
    next_substream!(s.servgen)
end

M/M/1 with processes
Out[4]:
restart (generic function with 1 method)
In [5]:
s = System("Waiting Times")
Out[5]:
System(Tally("Waiting Times",0,0.0,Inf,-Inf,0.0,0.0,0.0),Distributions.Exponential{Float64}(θ=0.5263157894736842),Distributions.Exponential{Float64}(θ=0.5),#undef,Full state of MRG32k3a generator:
Cg = [12345,12345,12345,12345,12345,12345]
Bg = [12345,12345,12345,12345,12345,12345]
Ig = [12345,12345,12345,12345,12345,12345],Full state of MRG32k3a generator:
Cg = [3692455944,1366884236,2968912127,335948734,4161675175,475798818]
Bg = [3692455944,1366884236,2968912127,335948734,4161675175,475798818]
Ig = [3692455944,1366884236,2968912127,335948734,4161675175,475798818])
In [9]:
meanWaits = Tally("Mean waiting times")

k = 30            # number of batches
batchSize = 10000 # number of clients per batch
Out[9]:
10000
In [6]:
# Allow a simulation with a fixed horizon

function source(env::Environment, s::System, limit::Float64)
    i = 0
    while (true)
        yield(Timeout(env, rand_dist(s.arrgen, s.arrival)))
        if (now(env) > limit) break end
        i += 1
        Process(env, customer, s, i, s.counter)
    end
end

# Allow a simulation with a fixed number of clients

function source(env::Environment, s::System, nCusts::Int64)
    for i = 1:nCusts
        yield(Timeout(env, rand_dist(s.arrgen, s.arrival)))
        Process(env, customer, s, i, s.counter)
    end
end

function customer(env::Environment, s::System, idx::Int, counter::Resource)
    # Record the arrival time in the system
    arrive = now(env)
    yield(Request(counter))
    # The simulation clock now contains the time when the client goes to the server.
    # Record the waiting time
    if (idx > warmup)
        wait = now(env) - arrive
        add(s.W, wait)
        # Check if we have a new batch
        if ((idx - warmup) % batchSize == 0)
            # Collect the information of the elapsed batch
            println("Batch. Average waiting time: ", average(s.W),
                    "Number of observations: ", nobs(s.W))
            add(meanWaits, average(s.W))
            init(s.W)
        end
    end
    yield(Timeout(env, rand_dist(s.servgen, s.service)))
    yield(Release(counter))
end
Out[6]:
customer (generic function with 1 method)
In [7]:
function onesim(s::System)    
    env = Environment()
    s.counter = Resource(env, 1)
    
    Process(env, source, s, warmup+k*batchSize)
    run(env)
end
Out[7]:
onesim (generic function with 1 method)
In [10]:
onesim(s)

Batch. Average waiting time: 8.621755467420913Number of observations: 10000
Batch. Average waiting time: 7.398029305875097Number of observations: 10000
Batch. Average waiting time: 9.046167606371142Number of observations: 10000
Batch. Average waiting time: 8.161309524502816Number of observations: 10000
Batch. Average waiting time: 10.867060883206099Number of observations: 10000
Batch. Average waiting time: 8.251634028351802Number of observations: 10000
Batch. Average waiting time: 10.74186672496718Number of observations: 10000
Batch. Average waiting time: 15.372885210881343Number of observations: 10000
Batch. Average waiting time: 6.5273934272794945Number of observations: 10000
Batch. Average waiting time: 10.858159453177263Number of observations: 10000
Batch. Average waiting time: 6.775052429493656Number of observations: 10000
Batch. Average waiting time: 8.583728105781306Number of observations: 10000
Batch. Average waiting time: 13.761091699137285Number of observations: 10000
Batch. Average waiting time: 11.761739529831425Number of observations: 10000
Batch. Average waiting time: 7.839947241178255Number of observations: 10000
Batch. Average waiting time: 13.175341674302198Number of observations: 10000
Batch. Average waiting time: 7.561304927962882Number of observations: 10000
Batch. Average waiting time: 12.375469083609321Number of observations: 10000
Batch. Average waiting time: 11.657675031947543Number of observations: 10000
Batch. Average waiting time: 9.88746221743163Number of observations: 10000
Batch. Average waiting time: 11.829701392170858Number of observations: 10000
Batch. Average waiting time: 6.581724195480052Number of observations: 10000
Batch. Average waiting time: 16.681128014725978Number of observations: 10000
Batch. Average waiting time: 6.591455453355905Number of observations: 10000
Batch. Average waiting time: 6.43485254499561Number of observations: 10000
Batch. Average waiting time: 5.867099371725206Number of observations: 10000
Batch. Average waiting time: 6.720846523451426Number of observations: 10000
Batch. Average waiting time: 9.056576222625335Number of observations: 10000
Batch. Average waiting time: 10.68069589932742Number of observations: 10000
Batch. Average waiting time: 9.090939486755488Number of observations: 10000
In [11]:
average(meanWaits)
Out[11]:
9.625336422577398
In [12]:
nobs(meanWaits)
Out[12]:
30
In [13]:
function ci_normal(n::Int64, mean::Float64, stdev::Float64, α::Float64)
    z = quantile(Normal(), 1-α/2)
    w = z*stdev/sqrt(n)
    # Lower bound
    l = mean - w
    # Upper bound
    u = mean + w
    
    return l, u
end
Out[13]:
ci_normal (generic function with 1 method)
In [14]:
ci_normal(nobs(meanWaits), average(meanWaits), stdev(meanWaits), 0.05)
Out[14]:
(8.62299505704682,10.627677788107976)
In [15]:
function ci_tdist(n::Int64, mean::Float64, stdev::Float64, α::Float64)
    z = quantile(TDist(n-1), 1-α/2)
    w = z*stdev/sqrt(n)
    # Lower bound
    l = mean - w
    # Upper bound
    u = mean + w
    
    return l, u
end
Out[15]:
ci_tdist (generic function with 1 method)
In [16]:
ci_tdist(nobs(meanWaits)-1, average(meanWaits), stdev(meanWaits), 0.05)
Out[16]:
(8.55985600482606,10.690816840328736)

Average waiting time according to M/M/1 formulas. In stationnary regime, the mean waiting time is $$ \overline{W} = \frac{\frac{\rho}{\mu}}{1-\rho} = \frac{\rho}{\mu-\lambda}. $$

In [17]:
w = (ρ/μ)/(1-ρ)
Out[17]:
9.499999999999991
In [ ]: