Lecture 25 (April 23, 2019)¶
In today’s lecture, we discuss how to apply antithetic method to discrete event system simulation.
In [0]:
import numpy as np
import matplotlib.pyplot as plt
import math
import scipy.special as spsp
import scipy.stats as spst
Review antithetic method¶
We will split the outputs into pairs. We will still generate the outcomes from the same distribution. However, we will find a strategy to have two outcomes negatively correlated.
Applying the antithetic method to discrete event system simulation.¶
Discrete event system is nothing but feeding samples from different distributions into a “system” to get the output.
In the basic single queueing system, we have two sets of samples to feed into the system
(1) the inter-arrival time for each customer
(2) the service time for each customer
If the samples are generated using the inverse transform method or from a normal distribution, it will be very easy to generate the original parts and the anththetic parts using the technique we discussed in class.
We the need to think about whether feeding the original part and the anthetic part will lead to negative correlation in the actual system output.
Let’s assume that we are interested in estimating the average time consumers stay in the system.
(1) larger inter-arrival time v.s. smaller inter-arrival time
Having larger inter-arrival time means consumers in this round of simulation are arrival at a lower pacer. In this case, we might be expecting fewer people to line up, which will lead to lower average time in the system.
If the inter-arrival times fed into the system are generally large, we should be expecting the antithetic part to be generally small. Thus, we will be able to produce negatively correlated system outputs.
(2) larger service time v.s. smaller service time
Larger service time has two impacts
a. More people might line up. Thus, consumers might wait for longer time before getting the service. b. The actual service time is part of the time in the system for a consumer. Longer service time will directly contribute to longer time in the system.
Thus, having when service time for two rounds of simulating being negatively correlated can also help to generate negatively correlated system output.
A Queueing system without variance reduction [first strategy]¶
Simulate the single queueing system. Assume that the system will only take the first 20 customers.
Assume
• the inter-arrival time follows $t_a \sim Exp(\lambda=10)$
• the service time follows $t_s \sim Exp(\lambda=11)$
In [0]:
## Arrival: poisson process with λ=100 for this hour.
## Service: Customers will be served in a counter based on first-come first serve with the service time following exponential distribution with λ=110
## Simulation stops when t>T
## Return the average waiting time for customers who finished receiving the service
def single_server():
#initialize clock
t=0
#initialize the states
#current number of customers in the system
n_c=0
#total arrivals so far
N_arrived=0
#total departures so far
N_departed=0
#generate the pending event times
#arrival rate
lmbda_a=10
#number of arrivals
size=20
Arrival_times=np.cumsum(-1/lmbda_a*np.log(np.random.rand(size)))
t_arrival=Arrival_times[N_arrived]
t_depart=np.inf
#service rate
lmbda_s=11
Service_times=-1/lmbda_s*np.log(np.random.rand(size))
#other initialization
##Collecting system output (history)
##array to store Arrival times
#arrival_ts=np.array([])
##array to store Departure times
Depart_times=np.array([])
#T=1
while (t0 and N_departed<20:
t_depart=t+Service_times[N_departed]
else:
t_depart=np.inf
#arrival_ts_cleaned=arrival_ts[arrival_ts<=1]
#epart_ts_cleaned=depart_ts[depart_ts<=1]
#avg_time=np.mean(depart_ts_cleaned-arrival_ts_cleaned[0:len(depart_ts_cleaned)])
##when simulation stops, we compute the average waiting time
return np.mean(Depart_times-Arrival_times)
ATs=[single_server() for i in range(2000)]
print(np.mean(ATs), np.std(ATs, ddof=1)/2000**0.5)
0.27949563921701065 0.003817227104079747
Applying anththetic method¶
• Feed in the inter-arrival time and the service time into the system to produce the output $x_{original}$
• Do the same thing except the inter-arrival time and (or) the service time will be produced using anththetic method. This will allow us to get another system output $x_{anti}$
• compute the average of the two $y=\frac{x_{original}+x_{anti}}{2}$
• repeat the above many times
• construct the confidence interval using all the $y$'s collected
In [0]:
## Arrival: poisson process with λ=100 for this hour.
## Service: Customers will be served in a counter based on first-come first serve with the service time following exponential distribution with λ=110
## Simulation stops when t>T
## Return the average waiting time for customers who finished receiving the service
def single_server(Arrival_times, Service_times):
#initialize clock
t=0
#initialize the states
#current number of customers in the system
n_c=0
#total arrivals so far
N_arrived=0
#total departures so far
N_departed=0
#generate the pending event times
#arrival rate
lmbda_a=10
#number of arrivals
size=20
#Arrival_times=np.cumsum(-1/lmbda_a*np.log(np.random.rand(size)))
t_arrival=Arrival_times[N_arrived]
t_depart=np.inf
#service rate
lmbda_s=11
#Service_times=-1/lmbda_s*np.log(np.random.rand(size))
#other initialization
##Collecting system output (history)
##array to store Arrival times
#arrival_ts=np.array([])
##array to store Departure times
Depart_times=np.array([])
#T=1
while (t0 and N_departed<20:
t_depart=t+Service_times[N_departed]
else:
t_depart=np.inf
#arrival_ts_cleaned=arrival_ts[arrival_ts<=1]
#epart_ts_cleaned=depart_ts[depart_ts<=1]
#avg_time=np.mean(depart_ts_cleaned-arrival_ts_cleaned[0:len(depart_ts_cleaned)])
##when simulation stops, we compute the average waiting time
return np.mean(Depart_times-Arrival_times)
In [0]:
def antithetic():
lmbda_a=10
lmbda_s=11
size=20
U_a=np.random.rand(size)
U_s=np.random.rand(size)
Arrival_times1=np.cumsum(-1/lmbda_a*np.log(U_a))
Service_times1=-1/lmbda_s*np.log(U_s)
Output1=single_server(Arrival_times1, Service_times1)
Arrival_times2=np.cumsum(-1/lmbda_a*np.log(1-U_a))
Service_times2=-1/lmbda_s*np.log(1-U_s)
Output2=single_server(Arrival_times2, Service_times2)
return (Output1+Output2)/2
NewATs=[antithetic() for i in range(1000)]
print(np.mean(NewATs), np.std(NewATs, ddof=1)/1000**0.5)
0.28360426651687404 0.0029601714466365517
In [0]:
## apply antithetic to only service time
def antithetic():
lmbda_a=10
lmbda_s=11
size=20
U_a=np.random.rand(size)
U_s=np.random.rand(size)
Arrival_times1=np.cumsum(-1/lmbda_a*np.log(U_a))
Service_times1=-1/lmbda_s*np.log(U_s)
Output1=single_server(Arrival_times1, Service_times1)
U_a=np.random.rand(size)
Arrival_times2=np.cumsum(-1/lmbda_a*np.log(U_a))
Service_times2=-1/lmbda_s*np.log(1-U_s)
Output2=single_server(Arrival_times2, Service_times2)
return (Output1+Output2)/2
NewATs=[antithetic() for i in range(1000)]
print(np.mean(NewATs), np.std(NewATs, ddof=1)/1000**0.5)
0.2830403173152644 0.003232664948030636
In [0]:
## apply antithetic to only service time
def antithetic():
lmbda_a=10
lmbda_s=11
size=20
U_a=np.random.rand(size)
U_s=np.random.rand(size)
Arrival_times1=np.cumsum(-1/lmbda_a*np.log(U_a))
Service_times1=-1/lmbda_s*np.log(U_s)
Output1=single_server(Arrival_times1, Service_times1)
Arrival_times2=np.cumsum(-1/lmbda_a*np.log(1-U_a))
U_s=np.random.rand(size)
Service_times2=-1/lmbda_s*np.log(U_s)
Output2=single_server(Arrival_times2, Service_times2)
return (Output1+Output2)/2
NewATs=[antithetic() for i in range(1000)]
print(np.mean(NewATs), np.std(NewATs, ddof=1)/1000**0.5)
0.2766721941502086 0.003442687230543482
Extension:¶
What if the number of arrivals is not fixed? For examle, what if we want to simulate the system until T?
In this case, we can first generate the original part like what we did before.
However, this might mean we generated too many or too few arrivals for the antithetic. If we generated too many arrivals, we can simply only keep the arrivals that happened before T. However, if we did not generate enough arrivals for the anththetic part, we will simply generate more arrivals until using the regular method.
For example, let's assume that we would like to simulate the system until t=2 instead of only allowing 20 customers.
In [0]:
def arrival_times():
arrival_times=np.array([])
u_a=np.array([])
t=0
while (t<2):
u=np.random.rand()
lmbda_a=10
t=t-1/lmbda_a*np.log(u)
arrival_times=np.append(arrival_times,t)
u_a=np.append(u_a,u)
return arrival_times[:-1],u_a[:-1]
arrival_timesA, u_a=arrival_times()
arrival_timesA
Out[0]:
array([0.0983999 , 0.18609273, 0.19528048, 0.24121975, 0.37047561,
0.38594442, 0.46595927, 0.4880456 , 0.76304874, 0.81744844,
0.95029794, 0.96242457, 1.18734201, 1.22940473, 1.33911053,
1.44986326, 1.48134867, 1.8771359 , 1.91756292])
If we try to generate the anti-thetic method, you might find we are either generating too few or too many arrivals for the antithetic part. Thus, we will need to modify.
In [0]:
lmbda_a=10
Arrival_times_anti=np.cumsum(-1/lmbda_a*np.log(1-u_a))
Arrival_times_anti
Out[0]:
array([0.04681066, 0.10060619, 0.34389478, 0.44377161, 0.47587031,
0.67013952, 0.72978917, 0.8916503 , 0.89825636, 0.98510733,
1.01587889, 1.23285757, 1.24400507, 1.35090117, 1.39152548,
1.43162921, 1.56252377, 1.56445263, 1.67455325])
We will construct a function to generate the antithetic part. We will feed in the random numbers we generated for the original arrivals. Once we get the anthtietic arrivals, we will possibly generate more arrivals if the last arrival happened before T. However, if the last arrival happened after T, we will only keep the arrivals happened before T.
In [0]:
def arrival_anti(u_a):
lmbda_a=10
Arrival_times_anti=np.cumsum(-1/lmbda_a*np.log(1-u_a))
last_time=np.max(Arrival_times_anti)
if last_time>2:
Arrival_times_anti=Arrival_times_anti[Arrival_times_anti<=2]
else:
additional_arrivals=np.random.poisson(lmbda_a*(2-last_time))
if additional_arrivals>0:
additional_times=np.sort(last_time+(2-last_time)*np.random.rand(additional_arrivals))
Arrival_times_anti=np.append(Arrival_times_anti,additional_times)
return Arrival_times_anti
arrival_anti(u_a)
Out[0]:
array([0.04681066, 0.10060619, 0.34389478, 0.44377161, 0.47587031,
0.67013952, 0.72978917, 0.8916503 , 0.89825636, 0.98510733,
1.01587889, 1.23285757, 1.24400507, 1.35090117, 1.39152548,
1.43162921, 1.56252377, 1.56445263, 1.67455325, 1.85029327,
1.90817644, 1.93891556])
After we determine the number of arrivals, we will need to do the same thing for the service time.
In [0]: