Simulation of a M/M/1 queue

Confidence interval on average waiting time using Delta theorem

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)

Seed for next MRG32k3a generator:
[12345,12345,12345,12345,12345,12345]

include("tally.jl")

stdev (generic function with 1 method)

include("tallystore.jl")

add (generic function with 2 methods)

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

restart (generic function with 1 method)

s = System("Waiting Times")

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

# Allow a simulation with a fixed number of clients

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

function source_fixed(env::Environment, s::System, counter::Resource, nCusts::Int64)
    for i = 1:nCusts
        yield(Timeout(env, rand_dist(s.arrgen, s.arrival)))
        Process(env, customer, s, i, 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)
    end
    yield(Timeout(env, rand_dist(s.servgen, s.service)))
    yield(Release(counter))
end

customer (generic function with 1 method)

function onesim(s::System)    
    env = Environment()
    s.counter = Resource(env, 1)
    
    # We consider a system over 24 hours
    Process(env, source, s, 48*60.0)
    run(env)
end

onesim (generic function with 1 method)

onesim(s)

average(s.W)

7.351199208498522

nobs(s.W)

4979

variance(s.W)

40.89354988685545

meanWaits = TallyStore("Average waiting times")

TallyStore(Float64[],Tally("Average waiting times",0,0.0,Inf,-Inf,0.0,0.0,0.0))

totWaits = TallyStore("Cumulative waiting times")

TallyStore(Float64[],Tally("Cumulative waiting times",0,0.0,Inf,-Inf,0.0,0.0,0.0))

nObs = TallyStore("Number of observations")

TallyStore(Float64[],Tally("Number of observations",0,0.0,Inf,-Inf,0.0,0.0,0.0))

nSim = 50

50

for n=1:nSim
    restart(s)
    onesim(s)
    w = average(s.W)
    N = nobs(s.W)
    add(meanWaits, w)
    add(totWaits, sum(s.W))
    add(nObs, Float64(N))
    println("Sim $n. Average waiting time: $w. Number of observations: $N")
end

Sim 1. Average waiting time: 7.923026555654201. Number of observations: 4952
Sim 2. Average waiting time: 4.783584620592791. Number of observations: 4848
Sim 3. Average waiting time: 8.98735624690155. Number of observations: 5049
Sim 4. Average waiting time: 4.963218831979749. Number of observations: 4954
Sim 5. Average waiting time: 9.862566899872856. Number of observations: 5061
Sim 6. Average waiting time: 5.035487595806768. Number of observations: 4917
Sim 7. Average waiting time: 7.6128357063099585. Number of observations: 4931
Sim 8. Average waiting time: 8.477087768486514. Number of observations: 4915
Sim 9. Average waiting time: 4.3736718361789055. Number of observations: 4974
Sim 10. Average waiting time: 5.888797320103889. Number of observations: 4869
Sim 11. Average waiting time: 20.348290444282092. Number of observations: 5165
Sim 12. Average waiting time: 12.019813220592987. Number of observations: 4965
Sim 13. Average waiting time: 6.075280686523191. Number of observations: 4978
Sim 14. Average waiting time: 17.306307493201903. Number of observations: 5016
Sim 15. Average waiting time: 6.094500598873442. Number of observations: 5004
Sim 16. Average waiting time: 4.774085492088705. Number of observations: 4918
Sim 17. Average waiting time: 12.423924628891914. Number of observations: 5026
Sim 18. Average waiting time: 12.57002163899818. Number of observations: 4981
Sim 19. Average waiting time: 12.79514951605715. Number of observations: 4898
Sim 20. Average waiting time: 11.866891963317252. Number of observations: 4890
Sim 21. Average waiting time: 8.414237846366767. Number of observations: 4881
Sim 22. Average waiting time: 10.044780856019159. Number of observations: 5017
Sim 23. Average waiting time: 5.7273214107258426. Number of observations: 4891
Sim 24. Average waiting time: 10.527895093453498. Number of observations: 4949
Sim 25. Average waiting time: 11.101721344330779. Number of observations: 5047
Sim 26. Average waiting time: 11.389101793960045. Number of observations: 4915
Sim 27. Average waiting time: 10.874519039599182. Number of observations: 5047
Sim 28. Average waiting time: 10.245226701708907. Number of observations: 5070
Sim 29. Average waiting time: 15.250667547180743. Number of observations: 5025
Sim 30. Average waiting time: 10.626136627909744. Number of observations: 5011
Sim 31. Average waiting time: 13.756591286274134. Number of observations: 5037
Sim 32. Average waiting time: 6.5187064464116355. Number of observations: 4963
Sim 33. Average waiting time: 10.565309628020445. Number of observations: 5034
Sim 34. Average waiting time: 11.703006995239729. Number of observations: 4914
Sim 35. Average waiting time: 7.35196303272571. Number of observations: 4953
Sim 36. Average waiting time: 10.618706380566575. Number of observations: 4984
Sim 37. Average waiting time: 10.423751746173027. Number of observations: 4990
Sim 38. Average waiting time: 16.773892247402678. Number of observations: 5054
Sim 39. Average waiting time: 7.5966351276654684. Number of observations: 5016
Sim 40. Average waiting time: 10.584350945121969. Number of observations: 4999
Sim 41. Average waiting time: 19.380706824729874. Number of observations: 5083
Sim 42. Average waiting time: 7.075997382440564. Number of observations: 5067
Sim 43. Average waiting time: 6.294637995866658. Number of observations: 4952
Sim 44. Average waiting time: 12.008451258896894. Number of observations: 5050
Sim 45. Average waiting time: 10.282145843874222. Number of observations: 4972
Sim 46. Average waiting time: 10.094191702669676. Number of observations: 5002
Sim 47. Average waiting time: 6.463346813055246. Number of observations: 4926
Sim 48. Average waiting time: 6.075763416330862. Number of observations: 4971
Sim 49. Average waiting time: 12.009256604866634. Number of observations: 4953
Sim 50. Average waiting time: 12.841816653591838. Number of observations: 5042

average(meanWaits.t)

9.936054713157848

variance(meanWaits.t)

13.911827919365427

var(meanWaits.obs)

13.911827919365418

totw = average(totWaits.t)

49641.9480761314

no = average(nObs.t)

4982.520000000002

totw/no

9.963221035968019

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

ci_normal (generic function with 1 method)

ci_normal(nobs(meanWaits.t), average(meanWaits.t), stdev(meanWaits.t), 0.05)

(8.902210294758007,10.969899131557689)

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

ci_tdist (generic function with 1 method)

ci_tdist(nobs(meanWaits.t)-1, average(meanWaits.t), stdev(meanWaits.t), 0.05)

(8.86471490016308,11.007394526152616)

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}.$$

w = (ρ/μ)/(1-ρ)

9.499999999999991

covariance = cov(totWaits.obs, nObs.obs)

735748.6102991527

cor(totWaits.obs, nObs.obs)

0.5958331591255

function ci_ratio(x, y::TallyStore, α::Float64)
    n = nobs(x.t)
    z = quantile(TDist(n-1), 1-α/2)

    μ1 = average(x.t)
    μ2 = average(y.t)
    # Ratio estimator
    if (μ2 == 0)
        ν = NaN
    else
        ν = μ1/μ2
    end
    
    # Variance of the ratio
    σ = sqrt((variance(x.t) + variance(y.t)*ν*ν - 2*cov(x.obs, y.obs)*ν)/(μ2*μ2))

    w = z*σ/sqrt(n)

    # Lower bound
    l = ν - w
    # Upper bound
    u = ν + w
    
    return ν, l, u
end

ci_ratio (generic function with 1 method)

ci_ratio(totWaits, nObs, 0.05)

(9.963221035968019,8.895957541375745,11.030484530560292)