【作業分析2】M/M/1モデルとM/G/1モデルの待ち時間をPythonでシミュレーション
前回はM/M/1とM/G/1の待ち行列モデルをExcelでシミュレーションしました。
【作業分析1】M/M/1モデルとM/G/1モデルの待ち時間をExcelでシミュレーション
今回はこれと全く同じアルゴリズムを、Numpyを使ってPythonで実装してみました。
シミュレーションのアルゴリズム
ここで前回のExcelによるアルゴリズムを復習しておきましょう。
M/M/1モデルの場合には以下のようになります。
貨物がランダムに到着します。
但し、長い期間取った集計結果から１５分間に平均して１個の貨物が到着するということがわかっています。
その貨物を作業する窓口は１つです。
作業にかかる時間も貨物によりバラバラですが、長い期間取った集計結果から１５分間に平均して１.２５個の貨物をさばけるということがわかっています。
ランダムに到着する貨物の到着間隔は指数分布になります。
またランダムな作業時間の長さも指数分布になります。
これをシミュレーションするために指数分布の乱数をそれぞれ１,０００個ずつ生成します。
但し、到着時間は到着間隔の累積になります。
各イベントの終了時間は
到着時間＋作業時間＋待ち時間
になります。
この終了時間が次の貨物の到着時間（次のイベントの開始時間）より遅ければ、すぐに次の貨物の作業を始められるので待ち時間はゼロです。
逆の場合は
次の貨物の到着時間―終了時間
が待ち時間になります。
以上をまとめると次のようになります。
M/G/1モデルの場合は作業時間を正規分布乱数に変えるだけです。
M/M/1モデルのシミュレーション
まずはNumpyをインポートしておきます。
import numpy as np次にExcelでやったのと同じように、到着時間と作業時間の指数分布乱数を生成するためのパラメータλとμに値を代入します。
lamda = 1
myu = 1.25そしてこのλとμを使って到着時間と作業時間の指数分布乱数をnp.randomの中のexponentialを使って生成します。
それぞれ１,０００個の乱数を生成します。
到着時間は指数分布の累積になるため、最後にcumsum()を付けることに注意です。
arrival_time = np.random.exponential(1 / lamda, 1000).cumsum()
service_time = np.random.exponential(1 / myu, 1000)またもう一つ注意すべき点として、パラメータにはλやμの逆数を指定しないといけません。
λやμは単位時間（この場合は１５分）に平均何回イベントが起こるかを表していますが、これは裏返せばイベントの平均発生間隔が1/λや1/μということです。
np.random.exponentialのパラメータは、この平均発生間隔を入れます。
これに対してPythonの標準ライブラリーとして入っているrandomモジュールの中のexpovariateで指数分布乱数を発生させる方法もあります。
この場合にはλやμをそのままパラメータとして指定します。
またexponentialのように一度に複数の乱数を生成することができないので、for文を使って１,０００回繰り返します。
arrival_time = []
for i in range(1000):
arr = random.expovariate(lamda)
arrival_time.append(arr)
service_time = []
for i in range(1000):
ser = random.expovariate(myu)
service_time.append(ser)これで到着時間として１,０００個の値が入ったリストと、作業時間として１,０００個の値が入った２つのリストができました。
次に、これら到着時間と作業時間のリストから平均待ち時間を計算する関数を定義します。
ここでExcelで作ったシミュレーションを思い出しましょう。
１つのイベントの終了時間は
到着時間＋作業時間＋待ち時間
でした。
またこの終了時間が次の貨物の到着時間（次のイベントの開始時間）より遅ければ、すぐに次の貨物の作業を始められるので待ち時間はゼロです。
逆の場合は
次の貨物の到着時間―終了時間
が待ち時間になりました。
これをPythonで表現するために、まずは１,０００回のイベントそれぞれの待ち時間と終了時間を入れる空のリストを作っておきます。
def waiting_queue(arrival_time, service_time):
waiting_time = []
leaving_time = []初回のイベントは待ち時間ゼロで、終了時間は
最初の到着時間＋最初の作業時間＋待ち時間（ゼロ）
になるので、リストの最初の成分は次のようになります。
waiting_time.append(0)
leaving_time.append(arrival_time[0] + service_time[0] + waiting_time[0])そして、それ以降の成分は同じ考え方で次のように計算できます。
for i in range(1, len(arrival_time)):
waiting_time.append(max(0, leaving_time[i-1] - arrival_time[i]))
leaving_time.append(arrival_time[i] + service_time[i] + waiting_time[i])あとは、waiting_timeにできた１,０００個の待ち時間の平均を計算すれば、それが１,０００回のイベントでシミュレーションした時の平均待ち時間です。
まとめると次のようになります。
import numpy as np
lamda = 1
myu = 1.25
arrival_time = np.random.exponential(1 / lamda, 1000).cumsum()
service_time = np.random.exponential(1 / myu, 1000)
def waiting_queue(arrival_time, service_time):
waiting_time = []
leaving_time = []
waiting_time.append(0)
leaving_time.append(arrival_time[0] + service_time[0] + waiting_time[0])
for i in range(1, len(arrival_time)):
waiting_time.append(max(0, leaving_time[i-1] - arrival_time[i]))
leaving_time.append(arrival_time[i] + service_time[i] + waiting_time[i])
waiting_mean = np.mean(waiting_time)
return waiting_mean
waiting_queue(arrival_time, service_time)計算結果は理論値である３.２に近い値にはなりますが、１.５から５.０くらいの間でばらつきます。
１,０００回の繰り返し回数ではばらつきが大きいということです。
ヒストグラムにしてみると次のようになりました。
import matplotlib.pyplot as plt
plt.hist(waiting_time, bins=10, density=True)
plt.show()
M/G/1モデルのシミュレーション
１,０００個の作業時間を正規分布乱数で生成する点だけが異なります。
次のようになります。
import numpy as np
lamda = 1
myu = 1.25
siguma = 0.2
arrival_time = np.random.exponential(1 / lamda, 1000).cumsum()
service_time = np.random.normal(1 / myu, siguma, 1000)
def waiting_queue(arrival_time, service_time):
waiting_time = []
leaving_time = []
waiting_time.append(0)
leaving_time.append(arrival_time[0] + service_time[0] + waiting_time[0])
for i in range(1, len(arrival_time)):
waiting_time.append(max(0, leaving_time[i-1] - arrival_time[i]))
leaving_time.append(arrival_time[i] + service_time[i] + waiting_time[i])
waiting_mean = np.mean(waiting_time)
return waiting_mean
waiting_queue(arrival_time, service_time)結果は理論値である１.６４近くになりますが、やはりばらつきます。
ヒストグラムにすると次のようになりました。
M/M/1よりバラツキが小さくなっていることがわかります。