待ち行列理論

待ち行列理論（Queueing theory）は、何らかのサービスを受けるために順番待ちをする人の列（待ち行列）などシステムの混雑現象を、確率の考え方を用いて数学的に理解する理論です。この理論は、デンマークの数学者 Agner Krarup Erlangが、1917年頃にコペンハーゲン電話会社で着信システムの記述モデルに関する研究が起源です。現在では、電気通信や交通、プロジェクトマネジメント、産業工学などの広い分野に応用されています。

　　　　　　　図：待ち行列理論を用いたスーパーマーケットのレジのモデル化

よく挙げられる例は、スーパーマーケットのレジや券売機です。レジにはバラバラに客が到着し、順番を待って会計をしていきます。これをモデル化して客が待たされる確率やその待ち時間、待ち行列の長さなどを評価します。このとき、客が互いに無関係にレジに来るものをポアソン分布で仮定するなど、統計の分野と非常に関係が深いです。

応用分野で特に有名なものが、オペレーションズ・リサーチ（operations research, OR）と呼ばれる、システムや組織の計画・管理・運用上の問題を解決する分野があります。 ORでは待ち行列理論を用いていろいろな問題を数学モデルで表現・解析し、意思決定の判断材料を提供しています。

待ち行列モデル

待ち行列モデルとは、待ち行列という現象を表す数学モデルのことを指します。モデル全体のイメージは、待ち行列システムに客が到着し、レジで会計などのサービスをうけた人から退出する、というものです。この中で、上記の図の待ち行列とレジとの部分を待ち行列システムと呼びます。

待ち行列モデルは、客の到着過程とサービス過程、窓口数、待ち行列の最大長さが変数です。特に最初の2つがこのモデルのポイントです。

到着過程 $X$ ：客の到着の挙動を確率でモデル化したもの。到着時間の間隔を確率分布として取り扱う
サービス過程 $Y$ ：サービスの完了の挙動を確率でモデル化したもの。完了時間の間隔を分布として取り扱う
窓口数 $S$ ：サービスを行う窓口の個数。レジの台数など
待ち行列の最大長さ $N$ ：システムが取り扱える待ち行列の長さ。レジの待機列など

前者の2つ（客の到着過程とサービス過程）には何らかの確率分布を適用し、後者の2つ（窓口数、待ち行列の最大長さ）には数値が入ります。

待ち行列モデルは、ケンドール記号と呼ばれる記法でモデルの構造を表します。

待ち行列モデルの構造（ケンドール記号）： $X/Y/S(N)$

説明	記号	意味
客はランダムに到着。到着過程ならポアソン分布、サービス過程なら指数分布	$M$	マルコフ過程¹
到着/サービスの時間間隔が一定	$D$	連続一様分布
特定の分布に従わず、独立性も仮定しない	$G$	一般の分布
特定の分布に従わず、独立性は仮定する	$G_I$	一般の分布
特殊な場合。 $k$ 個の指数分布の和。 $k=1$ で指数分布、 $k=\infty$ で連続一様分布	$E_k$	位相 $k$ のアーラン分布

また、到着過程の種類 $X$ とサービス過程 $Y$ に使用される記号とその意味は以下の表となります。（説明のため、先にまとめの表を示します）

図：ポアソン分布³　　　　　　　図：指数分布⁴　　　　　　　　　　図：連続一様分布²

たとえば、 $M/M/1(\infty)$ で表される待ち行列モデルは、到着過程がポアソン分布、サービス時間が指数分布、窓口の数が1個、待ち行列の最大長さが $\infty$ なものを示します。

以下に、それぞれの変数を詳しく見ていきます。

到着過程

到着過程は、単位時間あたりに到着する客数がどれだけかをモデル化したものです。客の到着時間の割合を分布として、以下の様に取り扱います。

ランダムに到着：ポアソン分布 $M$
到着の時間間隔が一定時間：連続一様分布² $D$
特定の分布に従わず、独立性も仮定しない分布：一般分布 $G$
特定の分布に従わず、独立性を仮定する分布： $G_I$

到着過程で最も一般的なケースである、客がランダムに到着する場合は、客の到着時間の割合をポアソン分布でモデル化できます。これは、ポアソン分布がある離散的な現象の単位時間あたりの生起回数の確率を示すことを、そのまま利用しています。

ポアソン分布の到着過程を「ポアソン到着」などと呼び、以下の3つの特徴を持ちます。

定常性：事象が起こる確率はどの時間帯でも同一
無記憶性：ある時間にある事象が起こる確率は、その時刻以前にその事象が起こった回数に依存しない
希少性：微小時間の間に事象が2回以上起こる確率は無視できるほど小さい

ポアソン分布の定式化は以下のとおりで、単位時間あたりに到着する客数 $x$ 、単位時間あたりに到着する平均客数を平均到着率 $\lambda$ [人/単位時間]、とします。この式はポアソン分布の定義そのものですので、ポアソン分布の平均と分散、最頻値は $\lambda$ となります。

$f(x) = \frac{{\lambda}^{x}}{x!} {e^{-\lambda}}$

たとえば、1時間に平均5人の客が来る場合、1時間あたりの客数は平均5のポアソン分布に従います。また、到着間隔は平均0.2時間の指数分布に従います。

サービス過程

サービス過程とは、サービスの完了の挙動をモデル化したものです。到着過程と同様に、サービスの完了を時間の確率分布として扱います。ただ、サービス過程では、単位時間あたりに何人分のサービスが完了できるかを重視するため、利用する確率分布が異なります。これを平均サービス率[人/単位時間]: $\mu$ と呼び、平均サービス時間はこの値の逆数となります。

サービスを受ける順番を先着順(FIFO, First In First Out)とすると、平均サービス率は指数分布に従うと考えられています。

指数分布

$g(x) = \lambda \times {e^{-\lambda x}}$

窓口数

待ち行列モデル自体は無限個の窓口数を取り扱うことができるが、現実的にはそうではない。窓口数を増やすことで待ち行列の長さを低減させたりサービスの質を向上させられるが、機材や人員コストが上昇する、というトレードオフがこれに当たる。なお、サービスの内容が違う場合にはサービス過程の種類を増やし、窓口数を増やすべきである。

待ち行列の最大長さ

待合室で待つことができる人数とサービスを受けている人数の総和であり、システムの容量のこと。待ち行列モデルでは、ことわりがない場合はこの長さを無限大 $N=\infty$ し、ケンドール記号の $N$ を省略することがある。

待ち行列システム

待ち行列システムは、前述のように、待ち行列モデルの待ち行列とサービスの部分を指すものです。ケンドール記号： $X/Y/S(N)$ の $S(N)$ の部分を指しており、実際のトレードオフを考える部分になります。たとえば、1時間に平均5人の客が来る場合、どのような待ち行列システムを設計すると最も平均待ち時間が少ないか？などを考えます。

待ち行列システムの良し悪しを判断するにはいくつかの指標があります。本講義資料では代表的な評価指標である、窓口利用率、待たされる確率、平均滞在客数、平均待ち行列長さ、平均滞在時間、平均待ち時間について見ていきます。なお、待ち行列モデルは最もシンプルなM/M/1( $\infty$ )を前提としています。 M/M/1(N)やM/M/S( $\infty$ )では値が異なります。

※ 講義資料ではすべての値の導出は行っていないので、詳細は参考書 9.1章待ち行列理論 P.137 を参照してください。また、参考書も微分方程式を用いた導出を省いているので、詳細は参考文献⁵などを参考にしてください。

代表的な評価指標

窓口利用率 $\rho$

窓口利用率は、以下に示すように、平均到着率 ${\lambda}$ を平均サービス率 ${\mu}$ で割ったものです。この値が1以上になると待ち行列は発散してしまいます。つまり、平均到着率が平均サービス率を上回る＝単位時間あたりの客の到着数がサービスの完了数を上回ると、待ち行列は解消されないと言えます。

$\displaystyle \rho = \frac{\lambda}{\mu}$

この指標は、待ち行列システムの良し悪しを左右する重要なもので、後述するいくつかの値がこの指標に依存しています。

待たされる確率 $P_q$

この指標は、とある客が待ち行列システムに到着したときに、先客がおりサービスを待たされる確率を示しています。先客がいるということは、システムに1人以上の客がいるという状態です。この状態をとる確率 $P_q$ は、以下の様に全確率(=1)から客が1人もいない確率 $P_0$ を引くことで求められます。また、 $P_0$ は $P_0 = 1 - \rho$ です。

$P_q = 1 - P_0 = 1 - (1 - \rho) = \rho$

平均滞在客数 $L$

この指標は、システム内に存在する客数の平均を示し、その客数が $n$ 人になる確率 $P_n$ は以下の様な期待値で表されます。 $...$ の部分の導出は、参考書 P.137 9.1 待ち行列理論を参照してください。

$L = \sum_{n=0}^{\infty}n P_n = P_0 \sum_{n=1}^{\infty} n \rho^{n} = ... = \frac{\rho}{1 - \rho}$

平均待ち行列長さ $L_q$

この指標は、待ち行列の長さの平均を示し、サービス中の1人を除いて客数が $(n - 1)$ 人になるときの確率 $P_n$ の期待値を用いて、以下の様に表されます。また、平均滞在客数 $L$ とは $L_q = \rho \times L$ もしくは $L_q = L - \rho$ の関係が成り立ちます。これは待ち行列の平均長さ $L_q$ が、システム内の平均滞在客数 $L$ から $\rho$ を引いたものであると言えます。

$L_q = \sum_{n=0}^{\infty} (n-1) P_n = P_0\sum_{n=0}^{\infty}n (n-1) \rho = ... = \frac{\rho^2}{1-\rho}$

平均滞在時間 $W$

この指標は、とある客がシステムに滞在する全体の時間を示したものです。具体的には、システムに到着してから、サービスを受けてシステムから出ていくまでの時間の平均です。これは、平均滞在客数 $L$ を平均到着率 $\lambda$ で割ったものになります。

$W = \frac{L}{\lambda} = \frac{\rho}{\lambda(1 - \rho)} = \frac{1}{\mu (1 - \rho)}$

平均待ち時間 $W_q$

この指標は、とある客がシステムに到着してからサービスを受けるまでに待ち行列に滞在する時間の平均を示します。これは、平均待ち時間 $W_q$ から平均サービス時間 $\frac{1}{\mu}$ を引いたもので表せます。

$W_q = W - \frac{1}{\mu} = \frac{1}{{\mu (1 - \rho)}} - \frac{1}{\mu } = \frac{\rho}{\mu (1 - \rho)}$

リトルの公式

平均待ち時間 $W_q$ と平均待ち行列長さ $L_q$ の間と待ち時間 $W$ と待ち行列長さ $L$ には、以下の様な関係式が成立します。これをリトルの公式（Little's formula）と呼び、システムに到着した客がいつかのサービスを受ける限り、どのような待ち行列モデルにおいても成立します。

$\begin{aligned} L_q =& \lambda \times W_q \\ L =& \lambda \times W \end{aligned}$

サンプルプログラム

参考文献⁶⁷で解説されているPythonプログラムを、Juliaに移植したものを以下に示します。到着時間とサービス時間、待ち時間、終了時間を、ヒストグラムと時系列データでプロットします。

MM1な待ち行列モデルのシミュレーション

# Simulate M/M/1(∞) waiting queue model
using Random
using Plots
using Distributions

# global variables
# 平均到着率：単位時間あたりに到着する平均客数
lambda = 5
# 平均サービス率：単位時間あたりに何人分のサービスが完了するか
mu = 10
# 窓口利用率
rho = lambda / mu
# the number of events
num_of_events = 1000;

# 到着過程は指数分布
arrival_time_non_cumsum = rand(Exponential(1 / lambda), num_of_events)
arrival_time = cumsum(arrival_time_non_cumsum)
# サービス過程は指数分布
service_time = rand(Exponential(1 / mu), num_of_events)
waiting_time = Float32[]
leaving_time = Float32[]
leaving_time_non_cumsum = Float32[]

# function definition
function waiting_queue(arrival_time, service_time)
    push!(waiting_time, 0)
    push!(leaving_time, arrival_time[1] + service_time[1] + waiting_time[1])

    for i in 2:length(arrival_time)
        push!(waiting_time, max(0, leaving_time[i-1] - arrival_time[i]))
        push!(leaving_time, arrival_time[i] + service_time[i] + waiting_time[i])
        push!(leaving_time_non_cumsum, arrival_time_non_cumsum[i] + service_time[i] + waiting_time[i])
    end
end

# Run the simulation
println("平均到着率（単位時間あたりに到着する平均客数）：", lambda)
println("平均サービス率（単位時間あたりに何人分のサービスが完了するか）：", mu)
println("The number of events：", num_of_events)

waiting_queue(arrival_time, service_time)
waiting_mean = mean(waiting_time)
println("Values from theory")
println("窓口利用率    : ",      )
println("待たされる確率 P_q theory    : ",       rho )
println("平均滞在客数 L    : ",  rho / (1 - rho))
println("平均待ち行列長さ L_q    : ",    rho*rho / (1 - rho))
println("平均待ち時間 W_q (theory)    : ",       rho / (mu * (1 - rho)))
println("Values from simulation")
println("平均待ち時間 W_q (simulation): ", waiting_mean)

# plots the simulation results
plt = []

push!(plt, histogram(arrival_time_non_cumsum, bins=100, title="arrival time (hist)"))
push!(plt, plot(arrival_time_non_cumsum, title="arrival time"))
push!(plt, histogram(leaving_time_non_cumsum, bins=100, title="leaving time (hist)"))
push!(plt, plot(leaving_time_non_cumsum, title="leaving time"))
push!(plt, histogram(service_time, bins=100, title="service time (hist)"))
push!(plt, plot(service_time, title="service time"))
push!(plt, histogram(waiting_time, bins=100, title="waiting time (hist)"))
push!(plt, plot(waiting_time, title="waiting time"))

plot(plt..., layout = (4, 2), size=(1600,1200))
savefig("queueing_theory.png")

実行結果

julia> include("queueing_theory.jl")
平均到着率（単位時間あたりに到着する平均客数）：5
平均サービス率（単位時間あたりに何人分のサービスが完了するか）：10
The number of events：1000
Values from theory
窓口利用率    :
待たされる確率 P_q theory    : 0.5
平均滞在客数 L    : 1.0
平均待ち行列長さ L_q    : 0.5
平均待ち時間 W_q (theory)    : 0.1
Values from simulation
平均待ち時間 W_q (simulation): 0.09337941
"/root/queueing_theory.png"

ヒストグラムから、それぞれが指数分布やポアソン分布になっていることがわかります。また、時系列データから、客それぞれの挙動やシステム内部の挙動がわかります。

演習

MM1な待ち行列モデルを以下のパラメータでシミュレーションを行い、その結果を示せ。
- 平均到着率=番号
- 平均サービス率=組
窓口数Sを2としたとき（M/M/2(∞)）のシミュレーションを行い、その結果を示せ。変更した部分とその理由を説明せよ。なお、窓口数を増やすと、それに反比例して平均サービス時間も減少することとする。

以上

参考文献

マルコフ過程とは、未来の挙動が現在のみに関係し、過去の挙動とは無関係である性質を持つ確率過程 ↩
連続一様分布 - Wikipedia ↩↩
ポアソン分布 - Wikipedia ↩
指数分布 - Wikipedia ↩
待ち行列理論（MM1モデル）の定理とその証明 - 高校数学の美しい物語 ↩
【作業分析1】M_M_1モデルとM_G_1モデルの待ち時間をExcelでシミュレーション - 物流業界の歩き方 ↩
【作業分析2】M_M_1モデルとM_G_1モデルの待ち時間をPythonでシミュレーション - 物流業界の歩き方 ↩