include("tally.jl")
include("tallystore.jl")
using SimJulia
using Distributions
using RandomStreams
const SEED = 12345
rand_dist(rng::MRG32k3a, Dist::Distribution) = quantile(Dist, rand(rng))
seeds = [SEED, SEED, SEED, SEED, SEED, SEED]
gen = MRG32k3aGen(seeds)
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
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
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
println("M/M/1 with processes")
# rates
λ = 150.0
μ = 300.0
ρ = λ/μ
# Working periods
periods = 13
# Period length
period_length = 1.0
# occupancy factor
B = [0.8, 1.0, 1.2, 1.4]
pB = [0.25, 0.55, 0.15, 0.05]
# We consider independent days, so there is no warmup to apply.
warmup = 0;
type System
W
B
arrival
service
counter
arrgen
servgen
day
calls
limit
function System(text::String)
s = new()
s.W = Tally(text)
s.B = 1.0
s.arrival = Exponential(1.0/λ)
s.service = Exponential(1.0/μ)
s.arrgen = next_stream(gen)
s.servgen = next_stream(gen)
s.day = next_stream(gen)
s.calls = 0
s.limit = 1.0/180.0
return s;
end
end
function restart(s::System)
init(s.W)
s.calls = 0
next_substream!(s.arrgen)
next_substream!(s.servgen)
next_substream!(s.day)
end
function set_busynessFactor(s::System, B::Float64)
s.B = B
s.arrival = Exponential(1.0/(B*λ))
end
# Compute the expectation of B
EB = pB'*B
# Expected arrivals
EA = periods*period_length*λ
s = System("Waiting Times")
# 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)
if (wait <= s.limit)
s.calls += 1
end
end
yield(Timeout(env, rand_dist(s.servgen, s.service)))
yield(Release(counter))
end
function gen_day(s::System)
u = rand(s.day)
i = 1
p = pB[1]
while (p < u)
i += 1
p += pB[i]
end
s.B = B[i]
s.arrival = Exponential(1.0/(s.B*λ))
# println("B: $(s.B), u: $u, i: $i")
end
function onesim(s::System)
env = Environment()
s.counter = Resource(env, 1)
# We consider a system over 24 hours
Process(env, source, s, periods*period_length)
run(env)
end
# We genereate the day-dependent quantities
gen_day(s)
onesim(s)
average(s.W)
nobs(s.W)
s.calls
# Service levels
SL = s.calls/nobs(s.W)
meanWaits = TallyStore("Average waiting times")
totWaits = TallyStore("Cumulative waiting times")
nObs = TallyStore("Number of observations")
nCalls = TallyStore("Number of calls answered in less than $(s.limit*3600) seconds")
nOver = TallyStore("Number of calls answered in more than $(s.limit*3600) seconds")
nSim = 100
for n=1:nSim
restart(s)
gen_day(s)
onesim(s)
w = average(s.W)
N = nobs(s.W)
add(meanWaits, w)
add(totWaits, sum(s.W))
add(nObs, Float64(N))
add(nCalls, Float64(s.calls))
add(nOver, Float64(N-s.calls))
println("Sim $n. Average waiting time: $w. Number of observations: $N")
end
average(nObs.t)
average(nCalls.t)
variance(nCalls.t)
variance(nCalls.t)/nSim
# Indirect estimator
EA-average(nOver.t)
variance(nOver.t)
covariance = cov(nCalls.obs, nObs.obs)
# Estimator with control variate
average(nCalls.t)-covariance*(average(nObs.t)-EA)/variance(nObs.t)
# Construct the new estimator
x_vc = Tally("Control Variate estimator")
for i = 1:nSim
add(x_vc, nCalls.obs[i]-covariance*(nObs.obs[i]-EA)/variance(nObs.t))
end
average(x_vc)
variance(x_vc)
# We first define the number of observations per strate
nB = length(B)
nDays = Array(Int64, nB)
for i = 1:nB
nDays[i] = round(nSim*pB[i])
end
nDays
Note: this allocation could be improved by ensuring that each strate has enough observations and that the sum is equal to the number of simulations previously set.
servedCalls = Array(Tally, nB)
for i=1:nB
servedCalls[i] = Tally("Served calls in time for strate %i")
end
for i = 1:nB
set_busynessFactor(s, B[i])
for n=1:nDays[i]
restart(s)
onesim(s)
add(servedCalls[i], Float64(s.calls))
end
end
# Strate details
for i = 1:nB
println("Strate $i. Mean: ", average(servedCalls[i]),
" Variance: ", variance(servedCalls[i]))
end
mean = 0.0
for i = 1:nB
mean += pB[i]*average(servedCalls[i])
end
mean
σ2 = 0.0
for i = 1:nB
σ2 += pB[i]*pB[i]*variance(servedCalls[i])/nobs(servedCalls[i])
end
σ2 *= nSim
total = 0
nSim = 200
for i=1:nB
nDays[i] = round(pB[i]*stdev(servedCalls[i]))
total += nDays[i]
end
for i = 1:nB
nDays[i] = round(nDays[i]*nSim/total)
end
nDays
for i=1:nB
servedCalls[i] = Tally("Served calls in time for strate %i")
end
for i = 1:nB
set_busynessFactor(s, B[i])
for n=1:nDays[i]
restart(s)
onesim(s)
add(servedCalls[i], Float64(s.calls))
end
end
# Strate details
for i = 1:nB
println("Strate $i. Mean: ", average(servedCalls[i]),
" Variance: ", variance(servedCalls[i]))
end
σ2 = 0.0
for i = 1:nB
σ2 += pB[i]*pB[i]*variance(servedCalls[i])/nobs(servedCalls[i])
end
σ2 *= nSim