Simpyを使ってM/M/Cモデルの待ち行列をシミュレーションしてみた
待ち時間はサービス窓口を増やせば短くなりますが、その分コストは上がります。
一般的に窓口を増やせば待ち時間は短くなりますが、段々とその効果は小さくなりいずれ頭打ちになります。
また窓口を増やす効果を測る指標（KPI）には待機台数や窓口の稼働率などもあります。
待機台数は平均だけではなく、最大待機台数も重要です。
なぜなら駐車場のキャパをオーバーしてしまうと、近隣に駐車して近所迷惑になってしまうからです。
窓口の稼働率は１００％がいいわけではありません。
トラックはランダムに到着するので、稼働率が１００％ということはかなりの待ち行列が発生していることになります。
でもだからといって低い方がいいわけでもありません。
７０―８０％くらいが妥当な線です。
実際にはいろいろな制約条件を満たす最適な窓口数を求めることになるわけですが、これをPythonのSimpyを使って試してみました。
Simpyを使ってM/M/1モデルの待ち行列をシミュレーションしてみた
ではトラックドックが１つ、つまり窓口が１つしかないM/M/1モデルでのシミュレーションを行いましたが、これでは最適な窓口数を求めることができないので、今回は複数窓口であるM/M/Cモデルのシミュレーションを行いました。
M/M/Cモデルについては【Excelで解く待ち行列M/M/Cモデル】物流現場で応用する方法を具体的に解説で公式を具体例で確認したので、同じ例を使いました。
「トラックが平均６台／時で物流センターに到着し、トラックドック１か所当たり５人のスタッフが担当していて、平均１.７台／時で処理できます。
最適なトラックドックの数は何ヶ所でしょうか？」
M/M/1モデルの延長でシミュレーションしてみる
Simpyを使ってM/M/1モデルの待ち行列をシミュレーションしてみたで作ったPythonのコードでは、トラックドックの数、トラックの平均到着率、平均作業スピードをパラメータで指定するようにしているので、そのままM/M/Cモデルのシミュレーションができます。
今回は次のようにパラメータを設定します。
# パラメータ設定
NUMBER_OF_BAYS = 4 # トラックドックの数
ARRIVAL_MEAN = 6 # トラックの平均到着率（台／時）
SERVICE_MEAN = 1.7 # 平均作業スピード（台／時）これでシミュレーションを実行すると、次のような結果が得られます。
1 arrived at 0.05
1 entered queue at 0.05
1 left queue at 0.05
2 arrived at 0.17
2 entered queue at 0.17
2 left queue at 0.17
3 arrived at 0.32
3 entered queue at 0.32
3 left queue at 0.32
3 finished unloading at 0.41
1 finished unloading at 0.64
4 arrived at 0.68
…行列での待機時間はリストduration_in_queueに、待機台数はリストlen_in_queueに入るようにコーディングしているのでこれらの要素の平均を計算すると次のようになります。
np.mean(duration_in_queue), np.mean(len_in_queue)
>>
(0.917848062217733, 6.308962618225491)この値はシミュレーションするたびに変わりますが、大体こんなもんです。
これらは平均行列待ち時間Wqと平均待機台数Lqとして公式で計算することができ、
Wq = 0.94
Lq = 5.64
になります。
詳しくはこちらの記事を参照下さい。
【Excelで解く待ち行列M/M/Cモデル】物流現場で応用する方法を具体的に解説
シミュレーション結果と比較してみると、平均待機台数が１割くらい大きくなっています。
これはなぜでしょうか？
いろいろ悩んだ結果、次のようにすれば理論通りの結果になることがわかりました。
M/M/Cモデルでは時間も考慮した待機台数を考えるべき
ドラックドックの数、つまりリソースは次のようにオブジェクト化して定義しています。
# リソースを定義する
operation_capa = simpy.Resource(env, capacity = NUMBER_OF_BAYS)このオブジェクトを使って
length = len(operation_capa.queue)により、その時の行列の長さ（待機台数）が取得できます。
前回のコードではこれを行列に並んだ時に取得していました。
しかしこれでは「どのくらいの時間、その待機台数だったか？」を考慮していません。
ある時は３台が待っていたけれども、すぐに次が来て４台になり、すべてのドックで荷卸しに時間がかかって３台になったのは１５分後だったという場合には、平均待機台数は単純に３と４の平均を取って３.５台ではなく、極めて４台に近いというべきでしょう。
このように時間も考慮した平均待機台数を求めるには、待機台数に時間を掛けた加重平均を取るのが適切です。
またそうするためには行列に並んだ時に待機台数を取得するだけでは不足で、行列から出た時の待機台数も取得すべきです。
このように待機台数が変化する瞬間の待機台数を時間と一緒に取得しておけば、あとは時間の差分を取ることによって「その待機台数だったのは何分間だったか」がわかります。
従って、まずは次のようにして行列に並んだ時と行列から出た時の待機台数をリストlen_in_queueに、またその時の時間をリストtm_in_queueに取得しておきます。
def operation(env, operation_capa, truck_number, time_of_arrival):
with operation_capa.request() as req:
print("{:.0f} entered queue at {:.2f}".format(truck_number, env.now))
length = len(operation_capa.queue) # 行列に入った時の待機台数
len_in_queue.append(length)
queue_in = env.now # 行列に入った時間
tm_in_queue.append(queue_in)
yield req
print("{:.0f} left queue at {:.2f}".format(truck_number, env.now))
length = len(operation_capa.queue) # 行列から出た時の待機台数
len_in_queue.append(length)
queue_out = env.now # 行列から出た時間
tm_in_queue.append(queue_out)次に、これら２つのリストは同じ長さなので、関連付けた計算がやりやすいように１つのデータフレームに合体します。
# 待機台数が変化した時の時間と台数を一つのデータフレームにまとめる
df1 = pd.DataFrame(tm_in_queue, columns = ['time'])
df2 = pd.DataFrame(len_in_queue, columns = ['len'])
df_length = pd.concat([df1, df2], axis = 1)できたデータフレームdf_lengthの中身を見てみると次のようになっています。
次に時間の差分をとって、delta_timeとして新たな列を追加します。
df_length['delta_time'] = df_length['time'].shift(-1) - df_length['time'] # 時間の差分を求める
df_length = df_length[0:-1] # 最終行を削除次の行の時間を引くことによって差分を計算していて、最終行は意味がなくなるため削除しています。
新たなデータフレームの中身は次のようになっています。
これで各待機台数の持続時間がわかったので、持続時間を重みとする待機台数の加重平均を計算します。
avg_len = np.average(df_length['len'], weights = df_length['delta_time']) # 時間を重みとする待機台数の加重平均を求めるこのようにして計算したavg_lenを見てみると、公式で計算した理論値である５.６４に近くなりました。
待機台数が７台以上になる確率を求めてみる
この物流センターの駐車場にはトラック６台分のスペースしかありません。
待機台数の平均が５.６４なので、駐車場のキャパオーバーになる可能性があります。
そこで待機台数が７台以上になる確率を求めてみます。
これは先程作ったデータフレームで、待機台数が７台以上になる時間の合計を求めて、全体のシミュレーション時間に対する割合を求めればよいので次のようにコーディングできます。
sum_over_queue = df_length[df_length['len']>= 7]['delta_time'].sum() # 待機台数オーバーの時間合計を求める
over_ratio = round((sum_over_queue / SIM_TIME) * 100, 2) # 上記時間の総時間に対する割合を求める
over_ratio
>>
29.99このように３０％近い時間帯は駐車場のキャパオーバーになることがわかります。
設備稼働率を求める
次に４つあるトラックドックの稼働状況を調べてみましょう。
何ヶ所のトラックドックが使われているかは、リソースオブジェクトoperation_capaを使って次のように取得することができます。
rsc = operation_capa.count # 行列から出た時に埋まっているドック数
resource_in_use.append(rsc)リソースの稼働状況が変化するのは行列から出た時なので、
yield reqの後に
queue_out = env.now # 行列から出た時間
tm_out_queue.append(queue_out)と一緒に記述することによって、行列から出た時間がリストtm_out_queueに、その時のリソースの専有数（トラックドックの稼働数）がリストresource_in_useに記録されます。
あとは先程と同じように２つのリストを１つのデータフレームに統合してから、時間差分の列を追加した後に、リソースの専有数と時間差分の加重平均を計算します。
# 稼働ドック数が変化した時の時間と稼働ドック数を一つのデータフレームにまとめる
df3 = pd.DataFrame(tm_out_queue, columns = ['time'])
df4 = pd.DataFrame(resource_in_use, columns = ['use'])
df_use = pd.concat([df3, df4], axis = 1)
df_use['delta_time'] = df_use['time'].shift(-1) - df_use['time'] # 時間の差分を求める
df_use = df_use[0:-1] # 最終行を削除
avg_use = np.average(df_use['use'], weights = df_use['delta_time']) # 時間を重みとする稼働ドック数の荷重平均を求める
avg_use
>>
3.784960300825919トラックドックは４ヶ所なので稼働率９４.５％、ほぼフル稼働であることがわかります。
最適なトラックドック数を求めてみる
以上のように、トラックドックが４ヶ所では明らかにキャパオーバーであることがわかります。
それでは最適なトラックドックは何ヶ所でしょうか？
これは何を基準にするかで変わってきますが、ここでは
- 平均待機時間が１０分以下
- 待機台数が７台以上になる時間が１％以下
- トラックドックの稼働率が７０％以上
を満たす最小トラックドック数が最適であるとします。
今までのシミュレーションからトラックドックが４ヶ所では少なすぎることが明らかなので、トラックドック数を４から７まで変化させて
- 平均待機時間
- 平均待機台数
- 待機台数が７台以上になる確率
- 設備稼働率
がどうなるかを調べてみます。
そこでfor文でトラックドックの数を４から７まで変化させながら、上記のKPIをデータフレームに記録していきます。
KPIを記録するデータフレームをresultsとすると、次のようにコーディングできます。
import simpy
import numpy as np
import random
import pandas as pd
# パラメータ設定
ARRIVAL_MEAN = 6 # トラックの平均到着率（台／時）
SERVICE_MEAN = 1.7 # 平均作業スピード（台／時）
SIM_TIME = 10000 # シミュレーション時間（時）
OVER_QUEUE = 7 # 許容できる待機台数
def truck_arrival(env, operation_capa):
truck_number = 0
while True:
yield env.timeout(random.expovariate(ARRIVAL_MEAN))
time_of_arrival = env.now # トラックが到着した時間
arrivals.append(time_of_arrival)
truck_number += 1
print("{:.0f} arrived at {:.2f}".format(truck_number, env.now))
env.process(operation(env, operation_capa, truck_number, time_of_arrival)) # 荷卸しプロセスを定義
def operation(env, operation_capa, truck_number, time_of_arrival):
with operation_capa.request() as req:
print("{:.0f} entered queue at {:.2f}".format(truck_number, env.now))
length = len(operation_capa.queue) # 行列に入った時の待機台数
len_in_queue.append(length)
queue_in = env.now # 行列に入った時間
tm_in_queue.append(queue_in)
yield req
print("{:.0f} left queue at {:.2f}".format(truck_number, env.now))
length = len(operation_capa.queue) # 行列から出た時の待機台数
len_in_queue.append(length)
queue_out = env.now # 行列から出た時間
tm_in_queue.append(queue_out)
tm_out_queue.append(queue_out)
rsc = operation_capa.count # 行列から出た時に埋まっているドック数
resource_in_use.append(rsc)
yield env.timeout(random.expovariate(SERVICE_MEAN))
time_of_departure = env.now # ドックから出た時間
departures.append(time_of_departure)
print("{:.0f} finished unloading at {:.2f}".format(truck_number, env.now))
time_in_system = time_of_departure - time_of_arrival # システム内の滞在時間（到着してからドックを離れるまでの時間）
duration_in_system.append(time_in_system)
time_in_service = time_of_departure - queue_out # サービス時間（ドックに入ってからドックを離れるまでの時間）
duration_in_service.append(time_in_service)
time_in_queue = queue_out - queue_in # 待機時間（行列に入ってから出るまでの時間）
duration_in_queue.append(time_in_queue)
results = pd.DataFrame(columns=['NUM_BAY', 'AVR_TIME', 'AVR_LEN', 'OVER_RATIO', 'OCC_RATE']) # トラックドック数、平均待機時間、平均待機台数、７台以上待機台数になる確率、稼働率を格納するデータフレーム
c = 0
for bay in range(4,8): # トラックドック数を４から７まで変化させる
c += 1
# リスト初期化
arrivals, departures = [], []
duration_in_queue, duration_in_service, duration_in_system = [], [], []
tm_in_queue, tm_out_queue = [], []
len_in_queue = []
resource_in_use = []
# シミュレーション環境を作る
env = simpy.Environment()
# リソースを定義する
operation_capa = simpy.Resource(env, capacity = bay)
# 到着プロセスを定義する
env.process(truck_arrival(env, operation_capa))
# シミュレーションを実行する
env.run(until = SIM_TIME)
results.loc[c, 'NUM_BAY'] = bay
results.loc[c, 'AVR_TIME'] = np.mean(duration_in_queue)
df1 = pd.DataFrame(tm_in_queue, columns = ['time'])
df2 = pd.DataFrame(len_in_queue, columns = ['len'])
df_length = pd.concat([df1, df2], axis = 1)
df_length['delta_time'] = df_length['time'].shift(-1) - df_length['time'] # 時間の差分を求める
df_length = df_length[0:-1] # 最終行を削除
avg_len = np.average(df_length['len'], weights = df_length['delta_time']) # 時間を重みとする待機台数の加平均
results.loc[c, 'AVR_LEN'] = avg_len
sum_over_queue = df_length[df_length['len']>= OVER_QUEUE]['delta_time'].sum() # 待機台数オーバーの時間合計を求める
over_ratio = round((sum_over_queue / SIM_TIME) * 100, 2) # 上記時間の総時間に対する割合を求める
results.loc[c, 'OVER_RATIO'] = over_ratio
df3 = pd.DataFrame(tm_out_queue, columns = ['time'])
df4 = pd.DataFrame(resource_in_use, columns = ['use'])
df_use = pd.concat([df3, df4], axis = 1)
df_use['delta_time'] = df_use['time'].shift(-1) - df_use['time'] # 時間の差分を求める
df_use = df_use[0:-1] # 最終行を削除
avg_use = np.average(df_use['use'], weights = df_use['delta_time']) # 時間を重みとする稼働ドック数の荷重平均を求める
occ_rate = avg_use / bay
results.loc[c, 'OCC_RATE'] = occ_rate結果を表示すると次のようになりました。
すべての基準を満たす最小トラックドック数は６ヶ所であることがわかりました。