In [30]:
import numpy as np
import matplotlib.pyplot as plt
import math
import scipy.special as spsp
import scipy.stats as spst
import copy
import simpy
%matplotlib inline
In [48]:
def arrival(env,desk,Outputs):
lmbda_a=4
i=0
while True:
inter_arrival=-1/lmbda_a*np.log(np.random.rand())
yield env.timeout(inter_arrival)
#one new arrival
i+=1
#this id will be shared with the serivce process
env.process(service(i,env,desk,Outputs ))
def service(i,env,desk,Outputs ):
lmbda_s=5
Outputs[“arrivals”].append(env.now)
rqt=desk.request()
yield rqt
service_time=-1/lmbda_s*np.log(np.random.rand())
yield env.timeout(service_time)
desk.release(rqt)
Outputs[“departures”].append(env.now)
Outputs[“Depart_ID”].append(i)
#np.random.seed(50)
def system():
Outputs={“arrivals”:[], “departures”:[], “VIP”:[], “Depart_ID”:[]}
env=simpy.Environment()
desk=simpy.Resource(env)
env.process(arrival(env,desk,Outputs))
T=5
env.run(until=T)
arrivals=np.array(Outputs[“arrivals”])
departs=np.array(Outputs[“departures”])
Depart_ID=np.array(Outputs[“Depart_ID”])
return np.mean(departs-arrivals[Depart_ID-1])
N=1000
ATs=[system() for i in range(N)]
ATs=np.array(ATs)
mean=np.mean(ATs)
std=np.std(ATs,ddof=1)
lower=mean-spst.norm.ppf(0.975)*std/N**0.5
upper=mean+spst.norm.ppf(0.975)*std/N**0.5
print(“95% confidence interval”, lower,upper)
print(“mean”, mean)
print(“width of confidence interval”,upper-lower)
95% confidence interval 0.4579061911689261 0.4950377619489201
mean 0.4764719765589231
width of confidence interval 0.03713157077999396
In [0]:
##apply antithetic to interarrival times
def arrival(env,desk,Outputs,RS_s, anti):
lmbda_a=4
i=0
while True:
if anti==1:
inter_arrival=-1/lmbda_a*np.log(1-RS_s.rand())
else:
inter_arrival=-1/lmbda_a*np.log(RS_s.rand())
inter_arrival=-1/lmbda_a*np.log(np.random.rand())
yield env.timeout(inter_arrival)
#one new arrival
i+=1
#this id will be shared with the serivce process
env.process(service(i,env,desk,Outputs,RS_s, anti ))
def service(i,env,desk,Outputs,RS_s, anti ):
lmbda_s=5
Outputs[“arrivals”].append(env.now)
rqt=desk.request()
yield rqt
u=RS_s.rand()
if anti==1:
u=1-u
service_time=-1/lmbda_s*np.log(u)
yield env.timeout(service_time)
desk.release(rqt)
Outputs[“departures”].append(env.now)
Outputs[“Depart_ID”].append(i)
#np.random.seed(50)
def system(RS_s, anti):
Outputs={“arrivals”:[], “departures”:[], “VIP”:[], “Depart_ID”:[]}
env=simpy.Environment()
desk=simpy.Resource(env)
env.process(arrival(env,desk,Outputs,RS_s, anti))
T=5
env.run(until=T)
arrivals=np.array(Outputs[“arrivals”])
departs=np.array(Outputs[“departures”])
Depart_ID=np.array(Outputs[“Depart_ID”])
return np.mean(departs-arrivals[Depart_ID-1])
def system_avg():
RS_s_original=np.random.RandomState()
Rs_s_anti=copy.copy(RS_s_original)
output1=system(RS_s_original, 0)
output2=system(Rs_s_anti,1)
return (output1+output2)/2
N=1000
N=int(1/2*(N))
ATs=[system_avg() for i in range(N)]
ATs=np.array(ATs)
mean=np.mean(ATs)
std=np.std(ATs,ddof=1)
lower=mean-spst.norm.ppf(0.975)*std/N**0.5
upper=mean+spst.norm.ppf(0.975)*std/N**0.5
print(“95% confidence interval”, lower,upper)
print(“mean”, mean)
print(“width of confidence interval”,upper-lower)
In [0]:
##apply antithetic to interarrival times and service times
def arrival(env,desk,Outputs,RS_a,RS_s, anti):
lmbda_a=4
i=0
while True:
u=RS_a.rand()
if anti==1:
u=1-u
inter_arrival=-1/lmbda_a*np.log(u)
yield env.timeout(inter_arrival)
#one new arrival
i+=1
#this id will be shared with the serivce process
env.process(service(i,env,desk,Outputs,RS_s, anti ))
def service(i,env,desk,Outputs,RS_s, anti ):
lmbda_s=5
Outputs[“arrivals”].append(env.now)
rqt=desk.request()
yield rqt
u=RS_s.rand()
if anti==1:
u=1-u
service_time=-1/lmbda_s*np.log(u)
yield env.timeout(service_time)
desk.release(rqt)
Outputs[“departures”].append(env.now)
Outputs[“Depart_ID”].append(i)
#np.random.seed(50)
def system(RS_s, RS_a,anti):
Outputs={“arrivals”:[], “departures”:[], “VIP”:[], “Depart_ID”:[]}
env=simpy.Environment()
desk=simpy.Resource(env)
env.process(arrival(env,desk,Outputs,RS_s,RS_a, anti))
T=5
env.run(until=T)
arrivals=np.array(Outputs[“arrivals”])
departs=np.array(Outputs[“departures”])
Depart_ID=np.array(Outputs[“Depart_ID”])
return np.mean(departs-arrivals[Depart_ID-1])
def system_avg():
RS_s_original=np.random.RandomState()
Rs_s_anti=copy.copy(RS_s_original)
RS_a_original=np.random.RandomState()
Rs_a_anti=copy.copy(RS_a_original)
output1=system(RS_s_original,RS_a_original, 0)
output2=system(Rs_s_anti,Rs_a_anti,1)
return (output1+output2)/2
N=1000
N=int(1/2*(N))
ATs=[system_avg() for i in range(N)]
ATs=np.array(ATs)
mean=np.mean(ATs)
std=np.std(ATs,ddof=1)
lower=mean-spst.norm.ppf(0.975)*std/N**0.5
upper=mean+spst.norm.ppf(0.975)*std/N**0.5
print(“95% confidence interval”, lower,upper)
print(“mean”, mean)
print(“width of confidence interval”,upper-lower)
In [0]: