待ち行列のM/M/1モデルとM/G/1モデルをExcelでシミュレーション
M/M/1モデルとは？
物流の作業分析に待ち行列は大変有用で、中でもM/M/1モデルは最も基礎となる待ち行列です。
簡単にいうと、作業場が一つしかないモデルのことです。
作業したい貨物はランダムに到着します。
ランダムとはいっても、１分間に平均２個とかいうように平均はわかっているものとします。
最初の１分間は２個、次の１分間は１個、次の１分間は０個、次の１分間は４個というような感じです。
作業時間もランダムに変わります。
簡単に作業できる貨物は３０秒で終わるけれども、難しい貨物は５分かかる。
でも平均すると１分間に２.５個というような感じです。
このような状況下では、難しい貨物を処理している間に貨物がたくさん到着すれば、すぐに作業場の前に行列ができてしまうことが容易に想像できます。
この行列の待ち時間が平均何分であるとか、最大何分といったことをM/M/1モデルを使えば推定することができます。
詳しくはこちらで解説しています。
このM/M/1モデルは、例えば宅配便の集荷センターの作業をシミュレーションするのに向いています。
宅配便の貨物サイズはある程度限定されていますが、小さなものから大きなものまで種々雑多です。
また伝票の手書き文字が読めないとか、間違っていたりすると作業に思わぬ時間がかかってしまうこともあります。
そのため作業時間がランダムという前提のM/M/1モデルがしっくりきます。
M/G/1モデルとは？
一方でメーカーのB to B物流ではどうでしょうか？
メーカー１社で取り扱う貨物の種類は限られています。
従って作業時間もある程度読めます。
完全には一定にならないかもしれませんが、１分±２０秒という程度は読めるかもしれません。
このようなケースに有用なのがM/G/1モデルです。
貨物の到着はランダムだけれども、作業時間は平均±標準偏差で表すモデルです。
M/G/1のMはMarkov（マルコフ）過程を表しており、到着はマルコフ過程（平たくいうとランダム）、GはGeneral distribution（一般分布）を表しており、作業時間は一般分布という意味です。
一般分布は平均と標準偏差（または分散）で表します。
最後の１は作業場（サービス窓口）が１つという意味です。
待ち時間を求める公式をExcelシミュレーションで検証する
M/M/1モデル、M/G/1モデルともに待ち時間を求める公式が知られています。
当然、待ち時間は状況によって変化するため、この待ち時間Wqは待ち時間の期待値（平たくいうと平均値）です。
M/M/1モデル
Wq＝ρ／（１－ρ）
M/G/1モデル
Lq＝（λ2σ2＋ρ2）／２（１－ρ）
Wq＝Lq／λ
ここで
λ：単位時間当たりに到着する貨物の平均
μ：単位時間当たりに処理できる貨物の平均
ρ：λ／μ
σ：単位時間当たりに処理できる貨物の標準偏差
Lq：行列に並んでいる台数の期待値
です。
M/G/1 queue – TKK より抜粋
このようにM/G/1モデルを解析的に解くのは少々難しいです。
ところがExcelの数値シミュレーションなら高度な数学の知識は不要です。
そこでExcelシミュレーションでどこまでできるか検証してみましょう。
シミュレーションの想定
M/M/1モデル
物流センターにトラックが１５分間に平均１台、ランダムに到着します。
貨物を荷卸しするための荷卸し場は１箇所しかなく、１５分間当たりの荷卸し台数はランダム、但し平均１.２５台ということが過去実績から判明しています。
トラックの待ち時間は平均で何分になるでしょうか？
M/G/1モデル
物流センターにトラックが１５分間に平均１台、ランダムに到着します。
貨物を荷卸しするための荷卸し場は１箇所しかなく、１５分間当たりの荷卸し台数は平均１.２５台、標準偏差０.２の正規分布で近似できることが過去実績から判明しています。
トラックの待ち時間は平均で何分になるでしょうか？
M/M/1モデルをシミュレーション
この場合は１５分間を単位時間とします。
すると、単位時間当たりに到着するトラックの平均台数は１台、つまりλ＝１です。
また単位時間当たりに処理できるトラックの平均台数は１.２５台、つまりμ＝１.２５です。
ρ＝λ／μ＝１／１.２５＝０.８になります。
その上で次のように入力していきます。
イベント
１,０００回のシミュレーションを行うために、１から１,０００までの通し番号を振ります。
到着間隔
トラックは単位時間（１５分）当たり平均１台、ランダムに到着します。
ランダムに発生する事象の時間間隔は指数分布に従います。
指数分布の式は次の通りです。
f(x) = λe–λx
x：単位時間数
λ：単位時間で事象が起こる平均回数
これをｘについて積分すると累積分布関数が求まり、次のようになります。
F(x) = \(\int λe^{-λx} dx\)
これを計算するとF(x)=Xとして、
\(X=1-e^{-λx}\)
となります。
そしてこれをｘについて解くと
x=-ln(1-X)/λ
となります。
Xは累積確率なので０～１の値を取ります。
従ってXに０～１の乱数を入れれば、指数分布ｘの乱数ができます。
Xが０～１の乱数であればXでも１－Xでも同じなので、
x=-ln(X)/λ
で指数分布の乱数を生成できます。
到着時間
到着間隔の累積が到着時間になります。
作業時間
これも到着間隔と同様の式で指数分布に従う乱数を生成できます。
但し、λの代わりにμを代入します。
終了時間
到着時間に作業時間と待ち時間を加えたものが、作業の終了時間になります。
待ち時間
待ち時間が発生するのはどんな時でしょうか？
前のトラックの作業が終了したのに、次のトラックが到着していない時ですね。
つまり、ある行の終了時間が次の行の到着時間より遅い時です。
そのような場合には、
終了時間―次の行の到着時間
で待ち時間が計算できます。
逆の場合、つまり上式が負になる時は待ち時間はゼロです。
従って、
Max（０, 終了時間―次の行の到着時間）
で待ち時間が計算できます。
これで１,０００台のトラックが到着した場合のシミュレーションができました。
待ち時間F列の平均を計算すると、３.１５５になりました。
単位は(１５分)です。
約４７分の待ち時間ということです。
これをM/M/1待ち行列モデルの公式で計算すると、
Wq
＝ρ／（μ－λ）
＝３.２
なので、近い数字になりました。
M/G/1モデルをシミュレーション
これもM/M/1モデルと同様にシミュレーションできますが、異なる点は下記の２点です。
平均作業速度
これはC3セルにあるMeanの計算式のことです。
１５分間に平均トラック１.２５台分の貨物を処理できるので、１台の貨物を処理するスピードは１／１.２５＝０.８（単位時間／台）です。
作業時間
作業時間は正規分布で近似できることが判明しているため、平均と標準偏差を使って正規分布の乱数を生成します。
NORM.INV関数を使って正規分布の乱数を生成する方法については、こちらで解説しています。
残りの式はM/M/1と同じです。
Ｆ列に計算される１,０００個の待ち時間の平均を計算すると１.５９１になりました。
単位は(１５分)です。
約２４分の待ち時間ということです。
これを先ほど出てきたM/G/1モデルの公式で計算してみると、
Lq
＝（λ2σ2＋ρ2）／２（１－ρ）
＝１.６４
Wq＝Lq／λ
＝１.６４
このように近しい値になりました。
同時に、M/M/1モデルと比べると、待ち時間がほぼ半減していることがわかります。
M/M/1では作業時間がランダムですが、M/G/1ではある程度の振れ幅の中で作業時間がわかり、不確実性が小さくなるためです。
まとめ
M/M/1の待ち行列を解析的に解く公式は広く知られていて比較的簡単です。
M/G/1の公式も知られていますが、公式を導くのは難解です。
ところが、到着時間や作業時間を乱数を使ってExcelシミュレーションすることは比較的容易です。
今回は１,０００回のシミュレーションを行いましたが、待ち時間の期待値はシミュレーションをするたびに結構変わりました。
これはシミュレーション回数を増やすことにより解決できます。
１,０００回を１,０００,０００回に増やしてもいいですし、１,０００回のシミュレーションを１,０００回繰り返して平均を取るのでも良いでしょう。
そうすることにより振れ幅は確実に小さくなりますが、Excelでは計算時間がかかりすぎます。
そこで次回はPythonで同様のシミュレーションを行ってみます。