Simpyの使い方｜M/M/1モデルの待ち行列シミュレーションをしてみた
前回はトラックの平均到着間隔や１台当たりの平均荷卸し時間がわかっている場合に、トラックドックでの待機時間が平均何分になるかをExcelやPythonを使ってシミュレーションしました。
【作業分析1】M/M/1モデルとM/G/1モデルの待ち時間をExcelでシミュレーション
【作業分析2】M/M/1モデルとM/G/1モデルの待ち時間をPythonでシミュレーション
ただこれでは待ち台数のシミュレーションができません。
いやできないことはないのですが、トリッキーなロジックをコーディングしないとだめで面倒です。
そこで、離散事象シミュレーションのライブラリーとして有名なSimpyを使って待ち行列のシミュレーションを行ってみました。
これなら各トラックが何分待ったかだけでなく、行列に並んだ時点で前に何台あったかまでシミュレーションすることができます。
Simpyとは？
物流現場では貨物は複数の処理がなされます。
例えば、入荷された貨物は検査された後、保管ロケーションに移動され、何日か後にはピッキングされます。
作業者も同じで、入荷貨物を検査したり、保管ロケーションへ移動したり、ピッキング作業を行ったりします。
つまりモノもヒトも時間的に不連続に、いろいろな処理をされたり処理したりします。
このようなのを離散事象システムといいます。
Simpyは離散事象システムのシミュレーションに使えるPythonのライブラリーです。
決定論的モデルと確率論的モデルの違い
システムには決定論的モデルと確率論的モデルがあります。
何回シミュレーションしても結果が同じモデルを決定論的、シミュレーションするたびに結果が変わるモデルを確率論的と言ったりします。
物流現場では同じ作業者が作業していても、毎日貨物が来るタイミングや物量も異なるので、作業が終わる時間は通常異なります。
つまり確率論的モデルです。
このように現実世界を近似したモデルは通常確率論的モデルになりますが、Simpyの使い方に少しずつ慣れていくために、まずは決定論的モデルのシミュレーションをしてみます。
Simpyで決定論的モデルをシミュレーションしてみる
トラックの荷卸しをする場面を考えます。
トラックバースは１つだけで、トラックバースが空になったらピッタリ１５分後に貨物を積んだトラックが入ってきます。
荷卸しは１５分間でピッタリ１.２５台分できます。
現実的ではありませんが、このようなケースを考えます。
まずSimpyとNumpyをインポートします。
import simpy
import numpy as np次にトラックの到着間隔と荷卸しの時間間隔を定義しますが、ここでは１５分間を単位時間とします。
するとトラックの到着間隔は単位時間当たり１台来るので、
１単位時間／１台＝１
です。
同様に荷卸しの時間間隔は１単位時間／１.２５台＝０.８です。
またシミュレーション時間は１０,０００単位時間として、これもSIM_TIMEとして定義しておきます。
# パラメータ設定
ARRIVAL_MEAN = 1/1
SERVICE_MEAN = 1/1.25
SIM_TIME = 10000次に、トラックが到着してから荷卸しが終わって出発する時間までの時間を記録しておくために、これらの結果を格納するための空のリストを作っておきましょう。
# リスト初期化
arrivals, departures = [], []
duration_in_system = []ここからSimpy 特有のコードを書いていきます。
Simpyには型があって、下記の３行は最低限必要なコードのようです。
# シミュレーション環境を作る
env = simpy.Environment()
# プロセスを定義する
env.process(operation(env))
# シミュレーションを実行する
env.run(until = SIM_TIME)その上でプロセスを関数として定義します。
ここではoperationという関数名で、各工程のフローを書いていきます。
これにも型があるようで、下記のように書きます。
while True:
時間の記録
yield env.timeout(工程１で費やす時間)
時間の記録
yield env.timeout(工程２で費やす時間)
時間の記録これで、工程１が始まった時間、工程１が終わった時間（＝工程２が始まった時間）、工程２が終わった時間が記録されます。
そして先に定義した１０,０００単位時間になるまでシミュレーションが実行されます。
今回はトラックの到着を工程１、荷卸しを工程２として次のようなコードが書けます。
def operation(env):
truck_number = 0
while True:
yield env.timeout(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))
yield env.timeout(SERVICE_MEAN) # 荷卸しにかかる時間
time_of_departure = env.now
departures.append(time_of_departure)
print("{:.0f} finished unloading at {:.2f}".format(truck_number, env.now))以上のコードを実行すると、次のような結果が表示されます。
1 arrived at 1.00
1 finished unloading at 1.80
2 arrived at 2.80
2 finished unloading at 3.60
3 arrived at 4.60
3 finished unloading at 5.40
4 arrived at 6.40
…
5554 finished unloading at 9997.20
5555 arrived at 9998.20
5555 finished unloading at 9999.00
5556 arrived at 10000.00１台目のトラックは１単位時間後に到着し、１.８単位時間後に荷卸しが終了しました。
２台目のトラックは２.８単位時間後に到着し、３.６単位時間後に荷卸しが終了しました。
これが延々と繰り返され、１０,０００単位時間後に５,５５６台目のトラックが到着したところでシミュレーションが終了しています。
尚、この結果は何回シミュレーションしても同じになります。
Simpyで確率論的モデルをシミュレーションしてみる
次に、トラックの到着間隔や荷卸しの時間間隔を確率変数にしたシミュレーションを実行してみましょう。
毎回スピードが確率的に変わる場合のシミュレーションということです。
到着時間と作業時間を確率的にする
そのために乱数を使いますが、全くデタラメの値の乱数ではなく、平均はわかっているものとします。
平均がわかっている場合のランダム乱数を生成する時には、指数分布を使います。
従って、工程フローを定義するoperation関数の中身は次のようになります。
def operation(env):
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))
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))結果は次のようになりました。
1 arrived at 0.78
1 finished unloading at 1.33
2 arrived at 1.51
2 finished unloading at 1.86
3 arrived at 3.51
3 finished unloading at 4.15
4 arrived at 4.24
…
5584 arrived at 9997.37
5584 finished unloading at 9997.49
5585 arrived at 9998.54
5585 finished unloading at 9998.73今度はトラックの到着間隔や荷卸し時間が毎回変わっています。
また、シミュレーションをするたびに、この結果は変わります。
リソースを定義して待ち時間を考慮する
次に待ち時間をシミュレーションできるようにしてみましょう。
待ち時間が発生するのは、リソースが限られているためです。
Simpyではリソースクラスを使って、この限られたリソースを定義できます。
今回の例ではトラックドックがリソースになり、キャパシティはドックが１つしかないので１です。
まずはキャパシティを１として、リソースクラスをインスタンス化しておきます。
変数ドックの数を表すNUMBER_OF_BAYSには、予め１を代入しておきます。
# リソースを定義する
operation_capa = simpy.Resource(env, capacity = NUMBER_OF_BAYS)そしてこのインスタンスを使って、プロセスを記述する関数の中で次のように記述します。
with インスタンス名.request() as req:
行列に並んだ時間の記録
length = len(インスタンス名.queue) # 行列に並んだ時の行列の長さを取得
yield req
行列から出た時間の記録こうすることにより、行列に並び始めた時間、その時の行列の長さ、行列から出た時の時間を取得することができます。
今回の例では次のようにコーディングしました。
def operation(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))
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))
queue_out = env.now
tm_out_queue.append(queue_out)
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))結果は次のようになりました。
1 arrived at 1.02
1 entered queue at 1.02
1 left queue at 1.02
1 finished unloading at 1.11
2 arrived at 2.65
2 entered queue at 2.65
2 left queue at 2.65
2 finished unloading at 3.40
3 arrived at 3.41
…
5570 finished unloading at 9994.87
5571 arrived at 9995.78
5571 entered queue at 9995.78
5571 left queue at 9995.78
5571 finished unloading at 9997.57各トラックについての到着時間、行列INした時間、行列OUTした時間、荷卸しが終わった時間が記録されていますが、何かが違います。
そうです、待ち時間がすべてゼロなのです。
これは前のトラックの荷卸しが終了しないと次のトラックが到着しないからです。
実際には、前のトラックの荷卸しに時間がかかっていれば、次のトラック、いや時には次の次のトラックまでもが容赦なく到着してくるはずです。
M/M/1モデルをシミュレーションする
それをシミュレーションするためには、到着のプロセスを独立させて１０,０００単位時間になるまで延々と繰り返します。
そして、トラックドックに並んで荷卸しが終わるまでのプロセスを別に作り、トラックが到着したらそちらへ受け渡すようにしてやります。
そうすればトラックが容赦なく到着する環境で、ドックが空いていればすぐに荷卸し、空いていなければ空くまで待ってから荷卸しするという状況を再現できます。
具体的にはメインのプログラムでtruck_arrivalというプロセス（トラックが到着するプロセス）を定義して１０,０００単位時間になるまで実行し、truck_arrivalの中で別のoperationというプロセス（ドックが空くまで待ってから荷卸しをするプロセス）を定義して受け渡してやることにします。
# パラメータ設定
NUMBER_OF_BAYS = 1
ARRIVAL_MEAN = 1
SERVICE_MEAN = 1.25
SIM_TIME = 10000
# リスト初期化
arrivals, departures = [], []
duration_in_queue, duration_in_service, duration_in_system = [], [], []
tm_in_queue, tm_out_queue = [], []
len_in_queue = []
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))
queue_out = env.now
tm_out_queue.append(queue_out)
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)
# シミュレーション環境を作る
env = simpy.Environment()
# リソースを定義する
operation_capa = simpy.Resource(env, capacity = NUMBER_OF_BAYS)
# プロセスを定義する
env.process(truck_arrival(env, operation_capa))
# シミュレーションを実行する
env.run(until = SIM_TIME)結果は次のようになりました。
1 arrived at 0.81
1 entered queue at 0.81
1 left queue at 0.81
2 arrived at 1.41
2 entered queue at 1.41
1 finished unloading at 1.56
2 left queue at 1.56
2 finished unloading at 1.87
3 arrived at 6.99
3 entered queue at 6.99
3 left queue at 6.99
4 arrived at 8.34
4 entered queue at 8.34
3 finished unloading at 8.43
…１台目が荷卸しをしている最中の１.４１単位時間に２台目のトラックが到着し、ドックが空く１.５６単位時間まで待っている様子がわかります。
３台目は２台目の荷卸しが終了してから６.９９単位時間に到着していますが、４台目はまたもや３台目の荷卸し中に到着しています。
このようにトラックがドックで待っている様子がシミュレーションできていることがわかります。
先のコードではduration_in_queueのリストに行列での待ち時間、duration_in_serviceに荷卸し時間、duration_in_systemにその合計（各トラック当たりの所要時間）を格納していますので、それらの合計を計算してみましょう。
np.mean(duration_in_queue), np.mean(duration_in_service), np.mean(duration_in_system)
>>
(3.3755474385469806, 0.8009464850021448, 4.176493923549126)これら３つの値の理論値はそれぞれ３.２、０.８、４.０なので、それらに近しい値になっています。
【待ち行列｜M/M/1モデル】５つの公式をわかりやすい言葉で解説します
また各トラックがドックに入るまでに前に何台のトラックが並んでいたかは、len_in_queueのリストに格納されています。
その平均と最大値を求めてみたら、次のようになりました。
np.mean(len_in_queue), np.max(len_in_queue)
>>
(4.192713989534998, 29)平均４台が前に並んでいたけれども、多い時には２９台も並んでいたことがわかりました。
M/M/1モデルの公式では平均４台というのは分かれど、最大２９台というのは分かりません。
シミュレーションの有用性がここにあります。
次回はこの結果を動画で表示させてみます。
FuncAnimationを使って待ち行列シミュレーションをアニメーションしてみた