In [22]:
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
import hyperopt
In [23]:
## no variance reduction
def arrival(env,server,Output):
lmbda_a=7
for i in range(10):
inter_arrival=np.random.exponential(1/lmbda_a)
yield env.timeout(inter_arrival)
Output[“arrival”].append(env.now)
env.process(service(env,server,Output))
def service(env,server,Output):
rqt=server.request()
yield rqt
lmbda_s=10
st=np.random.exponential(1/lmbda_s)
service_t=env.timeout(st)
yield service_t
server.release(rqt)
Output[“depart”].append(env.now)
def system():
env=simpy.Environment()
server=simpy.Resource(env,capacity=1)
Output={“arrival”:[],”depart”:[]}
env.process(arrival(env,server,Output))
env.run()
return np.mean(np.array(Output[“depart”])-np.array(Output[“arrival”]))
outputs=[system() for i in range(200)]
np.mean(outputs),np.std(outputs,ddof=1)/200**0.5
Out[23]:
(0.20053266250533128, 0.008949985647650191)
In [24]:
plt.hist(outputs)
plt.show()

In [25]:
##importance sampling
def arrival(env,server,Output,change_a,change_s):
lmbda_a=7*(1+change_a)
for i in range(10):
inter_arrival=np.random.exponential(1/lmbda_a)
yield env.timeout(inter_arrival)
Output[“arrival”].append(env.now)
Output[“inter_arrival”].append(inter_arrival)
env.process(service(env,server,Output,change_s))
def service(env,server,Output,change_s):
rqt=server.request()
yield rqt
lmbda_s=10*(1-change_s)
st=np.random.exponential(1/lmbda_s)
service_t=env.timeout(st)
yield service_t
server.release(rqt)
Output[“depart”].append(env.now)
Output[“service”].append(st)
def system(change_a,change_s):
env=simpy.Environment()
server=simpy.Resource(env,capacity=1)
Output={“arrival”:[],”depart”:[],”inter_arrival”:[],”service”:[]}
env.process(arrival(env,server,Output,change_a,change_s))
env.run()
lmbda_s=10
lmbda_a=7
pdf_ap=spst.expon.pdf(Output[“inter_arrival”],scale=1/lmbda_a)
pdf_aq=spst.expon.pdf(Output[“inter_arrival”],scale=1/(lmbda_a*(1+change_a)))
pdf_sp=spst.expon.pdf(Output[“service”],scale=1/lmbda_s)
pdf_sq=spst.expon.pdf(Output[“service”],scale=1/(lmbda_s*(1-change_s)))
adjust_a=np.prod(pdf_ap/pdf_aq)
adjust_s=np.prod(pdf_sp/pdf_sq)
return (np.mean(np.array(Output[“depart”])-np.array(Output[“arrival”])))*adjust_s*adjust_a
outputs=[system(0.01,0.01) for i in range(200)]
np.mean(outputs),np.std(outputs,ddof=1)/200**0.5
Out[25]:
(0.20408683761470628, 0.008018201034784738)
In [26]:
def obj(params):
a=params[0]
b=params[1]
diffs=[system(a,b) for i in range(200)]
return np.std(diffs,ddof=1)/200**0.5
trials=hyperopt.Trials()
best = hyperopt.fmin(fn=obj,
space=[hyperopt.hp.uniform(‘a’, 0, 0.2),hyperopt.hp.uniform(‘b’, 0, 0.2)],
## random search
#algo=hyperopt.rand.suggest,
##tpe search
algo=hyperopt.tpe.suggest,
max_evals=20,
#timeout=1,
trials=trials)
print(best)
100%|██████████| 20/20 [00:02<00:00, 9.75trial/s, best loss: 0.00408902927619375]
{'a': 0.09340743829322483, 'b': 0.1319408886323093}
In [27]:
outputs=[system(best["a"],best["b"]) for i in range(200)]
np.mean(outputs),np.std(outputs,ddof=1)/200**0.5
Out[27]:
(0.1954104058771727, 0.004037647133702903)
In [28]:
plt.hist(outputs)
Out[28]:
(array([13., 38., 49., 53., 21., 15., 5., 4., 1., 1.]),
array([0.08454616, 0.12000816, 0.15547016, 0.19093216, 0.22639415,
0.26185615, 0.29731815, 0.33278015, 0.36824215, 0.40370414,
0.43916614]),

In [29]:
outputs=[system(best[“a”],best[“b”]) for i in range(10000)]
np.mean(outputs),np.std(outputs,ddof=1)/10000**0.5
Out[29]:
(0.1987896125737201, 0.0006094003377414406)
In [30]:
In [31]:
In [32]:
In [33]:
In [34]:
In [35]:
In [36]:
In [37]:
In [38]:
In [39]:
In [40]: