目次
待ち行列理論とは、不特定多数の人が限られた施設を利用する時に生じる不確実な混雑現象を数理的に分析する理論として生まれた。混雑を滞留と言い換えれば、車の渋滞や、インターネットの繋がりにくさなど、人に限らず、モノや情報であっても、このような現象は観測されるので、今では、分析の対象はあらゆるものに広がっている。
共通の用語として、「客」が「窓口」に「到着」して「サービス」を受け「退去」する、窓口が先客に占有されていれば「待ち行列」を作って待つ、という用語を使う。単位時間当たりの到着客数を「到着率」、客が無数にいて窓口がサービスをし続けるとして、単位時間当たりの退去客数を「サービス率」という。到着してまだ退去していない客の数を「滞在客数」、滞在客数からサービス中の客を除いたものを「待ち行列長」という。客が到着してからサービスが開始されるまでの時間を「待ち時間」と呼び、それにサービス時間を加えたものを「滞在時間」と呼ぶ。
時間的にランダムに推移するシステムを数理的に分析する道具として、確率過程モデルがある。窓口を利用する客個々の事情を考えれば到着時刻はそれぞれ決まっているのかもしれないが、窓口からみれば次の客がいつ来るか、どれくらいのサービス処理量があるかは予測がつかないので、それを確率的な事象とみなすことは不自然ではない。そこで、客の到着間隔やサービス時間を確率変数と考え、それらによって決まる滞在客数や待ち時間を確率過程によってモデル化する。
モデル分析では、主として確率変動と混雑の関係を調べることにして、設定する条件が時間的に変化しないという定常性の仮定を置く。客が人の場合、例えば朝夕のラッシュのように生活のリズムに合わせて到着率は時間的に変動するのが普通であるが、これについては別途考えることにする。
混雑を左右するのは窓口への到着率と、窓口の処理能力（サービス率）の兼ね合いである。定常性を仮定したことにより、処理能力以上の客が到着すれば、平均的には混雑は増す一方だが、平均的に処理能力より少ない客が到着した場合でも、確率変動により一時的に滞留が発生することがあり、その数量分析に興味がある。サービス率が到着率を上回るという条件を定常条件という。
混雑を視覚的に表現するために、累積グラフを使う。時刻 t までに到着した客の数（累積到着数）を A(t)、退去した客の数（累積退去数）を D(t) として、A(t) と D(t) を同じグラフの上に重ねて描く。A(0) = D(0) = 0 とすれば、A(t) は D(t) 以上となり、その差は滞在客数を表す。したがって、A(t) と D(t) が離れていればいるほど混雑しているということになる。
以下では、確率過程モデル化したものに対し、乱数を使ってそのサンプルパスを生成し、数値的な分析を行う。
もっとも単純な例として、ランダムに到着する客が一定時間のサービス処理を受けて退去する、という例を考える。
「ランダムに到着する」ということを、いつ客が来てもおかしくない、すなわち、次の瞬間に客が到着する確率が常に一定、ということと考えると、数学的処理により、到着の点列はポアソン過程と呼ばれる確率過程に従うことが導かれる。このとき、到着間隔は指数分布に従い、単位時間の到着数はポアソン分布に従う。このことから、ランダム到着のことをポアソン到着ともいう。
一定時間のサービス処理ということは、到着後一定時間はそこに止まらなければいけないということである。もし、処理する人（窓口と呼ぶ）が無数にいるならば、あるいは客が来た時に窓口が常に空いているほど客がまばらにしか来なければ、客の退去過程は到着過程を一定時間ずらしただけなので、あまりむづかしい問題は生じない。窓口が一つしかない場合は、到着時に直ちにサービス処理が開始されない場合があるので、直前の客の退去時刻情報が必要になる。n 番目の客のサービス処理の開始時刻は、その客の到着時刻（X(n) と記す）と直前の客の退去時刻（Z(n-1) と記す）の大きい方で与えられる、ということから、到着時点列 {X(n)} が与えられれば、退去時点列 {Z(n)} は次のようにして「計算」することができることが分かる。ここで、C をサービス時間とした。
Z(n) = max(X(n), Z(n-1)) + C
Rを使って数値実験しよう。ポアソン過程は指数分布に従う間隔で１ずつ増加する確率過程なので、指数分布に従う乱数を生成し、その累積を計算すればポアソン到着客の点列を生成できる。ある客の退去時刻は、到着時刻と直前の客の退去時刻のどちらか大きいものにサービス時間を足したもの、として計算できる。その差が（もし正ならば）その客の待ち時間になる。また、ある時点で、累積到着数 A(t) から累積退去数 D(t) を引いたものがその時点における滞在客数である。
次のプログラム例は、指数乱数を使って到着時点列を生成し、退去時点列を計算するものである。rexp(n, a) はパラメータ a （平均の逆数）の指数分布に従う乱数を n 個生成する。cumsum 関数は引数の累積和を計算する。
到着時点で１増加する階段関数は累積到着数の時間変化を表し、客の到着の様子を視覚化する。同じグラフに退去時点で１増加する階段関数を重ねて描くと興味深いグラフが得られる。plot 関数の type="s" は、階段関数を描くオプション指定である。numeric(n) は n 個のゼロからなるベクトルを生成する。par(new=TRUE) は直前に描いたグラフの上に重ね描きすることを指定する。
> EA = 1.2 # 平均到着間隔 > ES = 1 # サービス時間 > n = 20 > arv = cumsum(rexp(n, 1/EA)) # 到着時点列 > # 退去時点列（０番目の客から数える） > dep = numeric(n+1) > for(i in 2:(n+1)) dep[i] = max(arv[i-1], dep[i-1]) + ES > # 累積到着数、累積退去数のグラフ > tmax = max(arv) > plot(c(0,arv), 0:n, xlim=c(0,tmax), ylim=c(0,n), type="s") > par(new=TRUE) > plot(dep, 0:n, xlim=c(0,tmax), ylim=c(0,n), type="s", col=2)到着時点列と退去時点列を整列し、到着時点では +1、退去時点では -1 という数を累積すると滞在客数が計算できる。次のプログラムは、滞在客数の時間変化を計算してグラフ化したものである。
> # 滞在客数 > ttab = cbind(c(arv, dep), c(rep(1, n), 0, rep(-1, n))) > ttab = ttab[order(ttab[,1]),] > stay = cumsum(ttab[,2]) > plot(ttab[,1], stay, type="s")
実行例を載せる。左図が累積到着客数 A(t) と累積退去数 (D(t) を重ねて描いたもの、右図は滞在客数の変化を描いたものである。
滞在客数の計算は次のようにすることもできる。A(t) は到着時点で１増加させ、D(t) は退去時点で１増加させるので、A(t) - D(t) は到着時点で +1、退去時点で -1 という数を累積したものと同じ、すなわち、左図の黒線と赤線の垂直方向の距離は右図のグラフの高さに等しいことが分かる。
時々刻々の滞在客数の変動は、たまたま観測された客の到着時点列から得られたもので特別な意味はなく、乱数を変えて別の到着時点列を生成して同じような実験を実施すれば全く違うパタンになる。さらに、ポアソン過程（指数分布）の性質から、過去のデータから将来の動きを予測することは不可能なので、この動きそのものを調べることはあまり意味がないが集団的な性質を調べることには意味がある。
客の滞在時間は退去時点から到着時点を引けば良いので、計算した時点列から平均滞在時間が計算できる。滞在客数は、到着あるいは退去が起きなければ客数は変わらないという事実を利用すれば、到着時点列と退去時点列を整列してその時点間の長さを滞在客数ごとに集計すれば、滞在客数の分布が計算できる。次のプログラムで、tapply(a,t,f) は、長さが同じ２つのベクトル a, t に対し、t の値が同じ a の要素を集めて関数 f を適用するという処理を実行する。この場合 t は客数、f は和なので、同じ客数をもつ区間を集計している。
> # 滞在時間 > (EW = mean(dep[-1] - arv)) [1] 1.758178 > hist(dep[-1] - arv, breaks=25, freq=FALSE) > # 滞在客数 > m = length(stay) > (EL = sum(stay[-m] * (ttab[-1,1] - ttab[-m,1])) / ttab[m,1]) [1] 1.546702 > Slen = tapply(ttab[-1,1] - ttab[-m,1], stay[-m], sum) / ttab[m,1] > plot(0:max(stay), Slen, type="h", ylim=c(0, max(Slen)), lwd=3)
次の実行例は、到着客数を 10000 とした場合の滞在時間（左図）と滞在客数（右図）の分布の例である。滞在客数の方は、客数を増やして実験してみても、安定した分布は得られない。
サービス処理が s 個の窓口で行われる場合も同様の計算が可能である。客は s 個の窓口に順番に割り振られ、ある窓口は s 人ごとの客を処理するとすれば、ある客のサービス開始時刻は、到着時刻と s 人前の客の退去時刻の大きい方となることを考慮すれば良い。つまり、到着頻度が s 分の１になった単一窓口モデルを s 個同時に扱うようなモノである。
s が大きければ到着時点列をずらしただけなのでランダムに近くなる。
一般に、窓口が一つしかないサービス施設に、（ランダムに）客が到着し、窓口が空いていれば直ちにサービスを受け、さもなければ、一列に行列を作ってサービスの順番を待つ、というサービス方式を単一窓口モデルといい、G/G/1 と略記する。最初と2番目の G は、それぞれ到着間隔、サービス時間の確率規則を表す記号、1 は窓口が一つ、ということを表す。 ポアソン到着の場合は M、サービス時間が一定の場合は D と書く約束なので、窓口一つで定型業務を行う施設は M/D/1 と書かれる。
平均到着間隔の逆数、あるいは単位時間当たりの到着数を到着率、平均サービス時間の逆数をサービス率と言い、到着率をサービス率で割ったものをトラフィック密度と言う。到着率がサービス率を超える、言い変えれば、客の要求量がサービス処理能力を超えるようであれば、待ち行列の長さは発散する。G/G/1 分析ではトラフィック密度が１未満であることを仮定する。定型業務処理の場合は、ある客の滞在時間を退去時刻と到着時刻の差として計算したが、到着時点列 {X(n} と退去時点列 {Z(n)} の関係を使って直接「計算」することもできる。n 番目の客の待ち時間を w(n)、サービス時間を s(n) とすると、w(n) = Z(n) - s(n) - X(n) なので、先に求めた退去時点列を計算する式を、待ち時間を表す式に書き換えることができる。
w(n) = max(0, Z(n-1) - X(n)) = max(max(0, w(n-1) + s(n-1) - x(n))
ここで、x(n) = X(n) - X(n-1) と置いた。これは n-1 番目の客と n 番目の客の到着間隔を表す。この待ち時間に関する漸化式はリンドレーの公式と呼ばれる。従って、客の到着間隔 {x(n)} と各客のサービス時間 {s(n)} が与えられれば、適当な初期値 w(1) を与えて、各客の待ち時間が計算できる。待ち時間が計算できれば、滞在時間は w(n) + s(n) で与えられることはいうまでもない。
一般の単一窓口モデルの数量分析を R で実行するプログラムは、定型業務処理の場合とほとんど同じになるが、ここではリンドレーの公式を導入したものに書き換えている。加えて、上の問題では待ち時間は滞在時間から定数のサービス時間を引いただけなので特に言及していないが、この問題では待ち時間と滞在時間は別々の分析対象になる。同様にして、サービス待ちの待ち行列の長さは滞在客数がいればサービス中の客を除けば良いのだが、いなければゼロのままなので、単純な計算にはならない。
関数にまとめた例を示す。
GG1 = function(rho, n, arrival, service) {
ES = 1
EA = ES / rho
# 乱数生成
iarv = arrival(n, EA)
serv = service(n, ES)
# 到着時点列
arv = cumsum(iarv)
# 待ち時間（リンドレーの公式）
wait = numeric(n)
for(i in 2:n) wait[i] = max(0, wait[i-1] + serv[i-1] - iarv[i])
# 退去時点列（０番目の客から数える）
dep = c(0, arv + wait + serv)
# 滞在客数のサンプルパス
ttab = cbind(c(arv, dep), c(rep(1, n), 0, rep(-1, n)))
ttab = ttab[order(ttab[,1]),]
stay = cumsum(ttab[,2])
m = length(stay)
# 滞在客数、待ち客数
Slen = tapply(ttab[-1,1] - ttab[-m,1], stay[-m], sum) / ttab[m,1]
EL = sum(Slen * (0:(length(Slen)-1)))
sdL = sqrt(sum(Slen * (0:(length(Slen)-1))^2) - EL^2)
Qlen = c(sum(Slen[1:2]), Slen[-(1:2)])
ELq = sum(Qlen * (0:(length(Qlen)-1)))
sdLq = sqrt(sum(Qlen * (0:(length(Qlen)-1))^2) - ELq^2)
# 滞在時間、待ち時間
EW = mean(dep[-1] - arv)
sdW = sd(dep[-1] - arv)
EWq = mean(wait)
sdWq = sd(wait)
stats = cbind(c(EW, sdW), c(EWq, sdWq), c(EL, sdL), c(ELq, sdLq),
c(mean(iarv), sd(iarv)), c(mean(serv), sd(serv)))
colnames(stats) = c("滞在時間", "待ち時間", "滞在客数", "待ち客数", "到着間隔", "処理時間")
rownames(stats) = c("平均", "標準偏差")
return(list(Stats=stats, NoWait=mean(wait==0), Qdist=Qlen, Wait=wait, Stay=wait+serv))
}
実行例は、例えば次のようになる。
> atime = function(n, EA) return(rexp(n, 1/EA)) > stime = function(n, ES) return(rexp(n, 1/ES)) > rho = 0.8 > n = 10000 > gg1 = GG1(rho, n, atime, stime)
> gg1$Stats 滞在時間 待ち時間 滞在客数 待ち客数 到着間隔 処理時間 平均 4.645231 3.663183 3.695586 2.914303 1.256073 0.9820475 標準偏差 4.417268 4.307087 4.043901 3.861016 1.283117 0.9621668
> gg1$NoWait [1] 0.2124
> gg1$Qdist[1:10] 2 3 4 5 6 7 0.38592368 0.12874946 0.09776885 0.07638170 0.06337529 0.05206056 0.04283518 8 9 10 0.03358814 0.02561130 0.02019374
> plot(0:(length(gg1$Qdist)-1), gg1$Qdist, type="h", lwd=3) > hist(gg1$Wait, breaks=25, freq=FALSE)
単一窓口モデルの窓口の数を複数にしたものを複数窓口モデルと言い、G/G/s と略記する。G/G/1 と同様、到着間隔、サービス時間の確率規則を G と表記している。s は窓口数である。サービスを待つ客は到着順に一列に並び、サービスは先着順に行われるものとする。単一窓口モデルと同様に、サービス能力が客の要求量を上回ることを仮定する。この場合のサービス能力は、各窓口の能力の s 倍となる。到着率を（一つの窓口の）サービス率を s 倍したもので割ったものをトラフィック密度と言い、トラフィック密度が１未満であることを仮定する。
客の待ち時間は、到着時点での先客の状態に依存する。窓口は s 個あるので、それらが最初に手空きになる時刻の最小値が比較の対象になる。到着時刻がそれより後ならば待ち時間はゼロ、さもなければ、その時刻まで待たされる。リンドレーの公式のように表そうとすると、n 番目の客の待ち時間は n-1 番目までの客の退去時刻を大きさの順に大きいもの並べて s 番目の客の退去時刻を D(n,s) として、次のように表される。
W(n) = max{0, D(n,s) - A(n)}
実際の計算では、n が大きくなるにつれて並べ替えの手間が増大し、効率的とは言えなくなる。
待ち時間の計算に必要なのは窓口の手空きになる時刻なので、s 個の補助変数を用意し、手空きになる時刻を更新していけばもう少し簡単に計算できる。B(n,i) を n 人の到着客を処理する時、窓口 i が初めて手空きになる時刻を表す変数とする。そうすると、n 番目の客の待ち時間は、その客の到着時刻と {B(n-1,i)} の最小値との差で決まる。
W(n) = max{0, B(n,k) - A(n)}, B(n,k) = A(n) + W(n) + S(n) ただし k = arg mini B(n-1,i)
ここで、arg mini g(i) は、g(i) が最小値をとる i の値、という意味である。
R でのプログラムは単一窓口の場合の待ち時間の計算だけを変えれば良い。その差分は次のようなものになる。B(n,i) を記憶するために server[1,] を定義する（server[2:3,] は稼働率を計算するために使っている）。which 関数は条件式を満たす配列の添え字を返す。server[2,] は窓口の手空き時間を累積している。あとは定義通りの計算。
server = matrix(0, 3, s)
departure = numeric(N) # 退去時点列
wait = numeric(N) # 待ち時間
for(i in 2:N) {
k = which(server[1,] == min(server[1,]))[1] # 最初にサービス終了するサーバ番号
if(arrival[i] > server[1,k]) {
server[2,k] = server[2,k] + arrival[i] - server[1,k] # サーバの空き時間統計
wait[i] = 0
} else {
wait[i] = server[1,k] - arrival[i]
}
departure[i] = arrival[i] + wait[i] + service[i] # 退去時刻
server[1,k] = departure[i] # サーバkの稼働終了時刻
}
これらを取り込んだ複数窓口モデルのシミュレーションを実施する関数は、例えば次のように書ける。
GGs = function(s, rho, N, arv, serv) {
arate = 1
srate = arate / s / rho
intarv = arv(N-1, 1/arate) # 到着間隔
service = serv(N, 1/srate) # サービス時間
server = matrix(0, 3, s) # サーバの稼働終了時刻と空き時間
# 待ち時間の計算
arrival = c(0, cumsum(intarv)) # 到着時点列
departure = numeric(N) # 退去時点列
wait = numeric(N) # 待ち時間
for(i in 2:N) {
k = which(server[1,] == min(server[1,]))[1] # 最初にサービス終了するサーバ番号
if(arrival[i] > server[1,k]) {
server[2,k] = server[2,k] + arrival[i] - server[1,k] # サーバの空き時間統計
wait[i] = 0
} else {
wait[i] = server[1,k] - arrival[i]
}
departure[i] = arrival[i] + wait[i] + service[i] # 退去時刻
server[1,k] = departure[i] # サーバkの稼働終了時刻
}
Qt = cbind(c(arrival, departure),c(rep(1,N), rep(-1,N)))
Qt = Qt[order(Qt[,1]),] # イベント時刻でソート
Qt[,2] = cumsum(Qt[,2]) # 系内客数過程
intv = Qt[-1,1] - Qt[-(2*N),1] # イベント間隔
Qlength = tapply(intv, Qt[-(2*N),2], sum) / Qt[2*N,1] # 系内客数の分布
# 統計量の計算
EQ = sum((0:(length(Qlength)-1)) * Qlength) # 平均系内客数
SQ = sqrt(sum((0:(length(Qlength)-1))^2 * Qlength) - EQ^2) # 系内客数の標準偏差
stats = rbind(c(mean(wait), EQ, mean(intarv), mean(service)),
c(sd(wait), SQ, sd(intarv), sd(service)))
colnames(stats) = c("待ち時間", "系内客数", "到着間隔", "サービス時間")
rownames(stats) = c("平均", "標準偏差")
server[3,] = 1 - server[2,] / departure[N] # サーバの稼働率
return(list(stats=stats, Qlength=Qlength, Ocrate=server[3,], wait=wait, Qt=Qt))
}
# 実行例
> arv = function(N, a) rexp(N, 1/a)
> serv = function(N, b) rep(b, N)
> s = 3
> rho = 0.8
> N = 1000
> MDs = GGs(s, rho, N, arv, serv)
> MDs$stats
待ち時間 系内客数 到着間隔 サービス時間
平均 1.324 3.712 1.0012 2.4
標準偏差 1.634 2.504 0.9554 0.0
> MDs$Qlength
0 1 2 3 4 5 6 7 8 9
0.0530396 0.1345377 0.1805632 0.1666669 0.1363197 0.1168377 0.0749325 0.0653193 0.0305763 0.0161519
10 11 12 13 14 15 16
0.0079264 0.0076051 0.0026246 0.0031374 0.0016305 0.0019526 0.0001786
> MDs$Ocrate
[1] 0.7989 0.7974 0.7971
サンプル１（系内客数）
サンプル２（待ち時間）
タクシーのように、サービスする側も移動するという場合を考える。ランダムに到着する客とランダムに到着するタクシーの間で生じる客とタクシーの待ち行列長の時間変化を実感してみる。
ある場所に２種類のモノがポアソン到着に従って到着し、２種類のモノが揃ったら同時に立ち去る、揃わないモノはその場所で、マッチング出来るまで待つ、という状況を考える。このとき、その場所にいる各種類のモノの数がどのように時間変化するか知りたい。２種類のモノを、タクシーとその利用客（一人客）と考えると、冒頭の問題になる。
Rを使って数値実験しよう。ポアソン過程は指数分布に従う間隔で点が生成されるので、指数分布に従う乱数を生成し、その累積を計算すればポアソン到着客の点列を生成できる。２つのポアソン過程を合わせたものもポアソン過程になるので、一時に生成することができる。合成したものを到着率の違いで比例配分すれば、個々のポアソン到着点列が得られる。
客が到着したら +1、タクシーが到着したら -1 という数を累積すると、プラスならば客が待ち行列を作り、マイナスならばタクシーが待ち行列を作ることになる。結局、ランダム時点で、両者の到着率に比例した確率で増減するランダムウォーク（ただし、ジャンプ時点以外は状態を変えない）と同じ動きをする。
次のプログラム例は、２つのポアソン到着を重ねたポアソン過程のサンプルパスと、到着時点で両者の到着率の違いに比例した確率で決まる +1, -1 の数列を生成し、後者の数列を累積したものを計算するものである。
う２種類の客がいて、
# 客の到着率λとタクシーの到着率μを与えて、n人、n台到着したときの客の待ち行列長を計算する
# 出力は到着時点直後の行列長（ベクトル）長さ2n+1
taxi = function(lmd=1,mu=1,n=100) {
customer = c(0,cumsum(rexp(n,lmd)))
taxi = c(0,cumsum(rexp(n,mu)))
tt = cbind(c(customer,taxi[-1]),c(0,rep(1,n),rep(-1,n)))
tt = tt[order(tt[,1]),] # 時間列でソート
tt = cbind(tt,cumsum(tt[,2])) # 行列長を計算
return (tt)
}
##### 一様化（ポアソン到着ならば、時点列の和集合もポアソン到着する）
taxiunif = function(lmd=1,mu=1,n=20) {
event = sample(c(-1,1),n,prob=c(lmd,mu),replace=T) # 事象種類
etime = c(0,cumsum(rexp(n,lmd+mu))) # 事象時刻
tt = cbind(etime, c(0,event))
tt = cbind(tt,cumsum(tt[,2])) # 行列長を計算
return(tt)
}
taxiplot = function(tt) {
plot(tt[,1], tt[,3], type="s", main="タクシー乗り場の行列") # 階段関数でプロット
abline(h=0, col=2)
}
## サンプル
n = 200
taxiplot(taxi(1,1,n))
taxiplot(taxiunif(n=100))
## 重ね描き
plot(0,0,xlim=c(0,n), ylim=c(-3,3)*sqrt(n), type="n")
tt = taxi(1,1,n)
par(new="T")
plot(tt[,1], tt[,3],xlim=c(0,n), ylim=c(-3,3)*sqrt(n), type="s")
# 上の３行を繰り返す
### 最終結果の度数分布（２項分布）
q = sapply(rep(200,1000), function(n) taxi(1,1,n)[n,3])
hist(q,main=paste(mean(q),sd(q)))
サンプル１
サンプル２（集団的性質）
不特定多数の人が巨大な施設（たとえばテーマパーク）を利用する場面を想定し、客の到着率や客の滞在時間が時間帯ごとに変化する場合、施設に滞在している客がどのように変化するか、混雑のピークはいつなのか、を調べたい。
理論的な分析はむづかしいので、シナリオを設定して、滞在客数の時間推移を計算する。細かいランダム要因は無視して、時間帯を区切り、各時間帯ごとの到着数を与える。
到着客数は昼頃にピーク（ガンマ分布で表現）
滞在時間はガンマ分布に従う。
## 巨大施設の滞在客数 ## 時間帯ごとの到着客数はガンマ分布に従う a = 3; b = 1 # 昼頃に来場客数のピークが来る t = 8; dt = 1/6 # 全８時間、集計は１０分ごとに n = t/dt N = 40000 # 一日の入場者数（概算） nc = N * (pgamma((2:n)*dt,a,b) - pgamma((1:(n-1))*dt,a,b)) plot(c(1:(n-1))*dt,nc,type="h") ## 滞在時間はガンマ分布に従う aa = 3; bb = 3 mc = nc[1] * (1-pgamma(c(1:(n-1))*dt,aa,bb)) ## 各時間帯ごとの滞在者数を累積 for(i in 2:(n-1)) mc = mc + c(rep(0,i-1), nc[i] * (1-pgamma(c(1:(n-i))*dt,aa,bb))) t0 = 10 plot(t0+(1:(n-1))*dt, mc, type="h") # 滞在者数 points(t0+(1:(n-1))*dt, nc) # 到着数 text(14.5,7000,"滞在客数") legend(16,5000,"到着客数", pch=1)
実行結果：
サービス率が到着率を上回らなければいけないという定常条件は、窓口が複数になったとしても変わらない。
混雑がほどほどに収まるために、「定常条件」を満たしているにもかかわらず、待ち行列の長さが発散するという場合が起こりうる。
客は４種類のサービス（P,Q,R,S としよう）を順番に受けて退去する。窓口は２つ（A, B としよう）あり、窓口 A は P と S、窓口 B は Q と R を処理する。S, Q はそれぞれ P, R に優先してサービスされる（ただし、非割り込み優先）。客が（P から処理されるので）窓口 A に到着し、A で S の処理が終わると退去する。
例えば、４種類のサービスを４台の機械、窓口をそれらの機械を扱う作業員とすると、U 字ラインを２人の作業員が担当する、という状況を想定すれば良い。ただし、機械と言っても作業員が付きっ切りで面倒を見なければいけないので、作業員がいない機械は停止している。
説明の簡単のために処理時間は一定として、それぞれ p,q,r,s と記す。
１）窓口 A が処理 S を実行中に新たな客が到着したとすると、処理 P を実行できないので、客は窓口 A の前で待たされ、窓口 B へ移動できない。処理 Q, R を待つ客がいないとすれば、窓口 B は空いたままである。
２）逆に、窓口 B が処理 Q を実行中で、窓口 A には処理 S を待つ客がいないとすると、到着した客は窓口 A で P の処理を受けて窓口 B で処理待ちになる。処理 Q が終了した時、もし窓口 B に処理 Q を待つ客がいれば、処理 Q が優先されるため、処理 Q を終えた客は処理 R を待つために、窓口 A へ移動できない。このとき、窓口 A は処理 S を実行することができない。
さらに分かりやすくするために、処理 P, R が瞬時に終わるようなものとすれば、処理 Q が続いている間、（処理 R が実行されないので）窓口 B から窓口 A への流入は止まり、処理 S が続いている間、（処理 P が実行されないので）窓口 B への流入が止まる、という構造になっている。したがって、各窓口の処理能力は両窓口とも他の窓口の稼働状況によって処理能力が阻害され、待ち客がいて窓口が空いているのにサービスされずに待たされる状況が起こりうる。
次のプログラム例は、
A は空いたままになる。
一つの窓口が２種類のサービスを提供する場合、
同じ機械で別工程を加工する場合、優先権の与え方によって、定常条件を満たしているのに、待ち行列長が発散することがあり得る。
4つの工程P,Q,R,Sが必要な生産ラインを考える。機械はA,Bの2台で、部品はAに到着する。AでPの処理、BでQの処理をした後、再びBでRの処理を して、最後にAでSの処理をした後退去する。
機械AではSが優先処理（非割り込み）され、機械BではRが優先処理される。
このとき、しかかり部品はどこにどれくらいたまるか？
interarrival = function(ar) return(-ar*log(runif(1)))
# interarrival = function(ar) return(ar)
# リエントラント待ち行列モデルの計算
reentrant = function(ar=10,sv1=1,sv2=7,n=100) {
sv = c(sv1,sv2,sv1,sv2)
dp = c(Inf,Inf,Inf,Inf,0)
Q = c(0,0,0,0)
rcd = rep(0,5)
for(i in 1:n) {
now = min(dp) # 事象時刻
k = which(now == dp)[1] # 2イベントが同時に起きた場合の対応
if(k == 5) { # 新規客の到着
dp[5] = dp[5]+interarrival(ar) # 次の到着時刻
Q[1] = Q[1]+1
if(Q[1]+Q[4] == 1) dp[1] = now + sv[1] # 待ち行列無し
} else {Q[k] = Q[k]-1; if(Q[k] == 0) dp[k] = Inf} # サービス終了
if(k < 4) Q[k+1] = Q[k+1]+1 # 次のステージに進む
if(k == 1 | k == 4) { # 入り口出口の窓口
if(Q[4] > 0) {dp[4] = now + sv[4] # 優先客あり
} else if(Q[1] > 0) dp[1] = now + sv[1] # 非優先客あり（優先客無し）
}
if(k == 1 && Q[2]+Q[3] == 1) dp[2] = now + sv[2] # 次の窓口でサービス開始
if(k == 2 | k == 3) { # 中間の窓口
if(Q[2] > 0) {dp[2] = now + sv[2] # 優先客あり
} else if(Q[3] > 0) dp[3] = now + sv[3] # 非優先客あり（優先客無し）
}
if(k == 3 && Q[1]+Q[4] == 1) dp[4] = now + sv[4] # 次の窓口でサービス開始
rcd = rbind(rcd,c(now,Q))
}
return(rcd)
}
# テスト計算
ar = 10
sv1 = 1; sv2 = 8
n = 1000
rq = reentrant(ar,sv1,sv2,n)
# 非優先クラスの仕掛かり在庫の時間的推移
plot(rq[,1],rq[,2],type="s")
plot(rq[,1],rq[,4],type="s")
lines(rq[,1],rq[,2],col=2)
もし、処理Sの仕掛品がある場合は、到着した部品は処理されずに待ち行列で待つので、到着過程がそのまま待ち行列に累積されていく。
機械Bで処理すべき新たな仕掛品が到着しないので、処理Q,Rが実行され、機械Bは仕掛品がなくなる。
逆に処理Qの仕掛品がある場合は、処理Rが実行されないので、機械Bに処理R待ちの仕掛品がたまり、機械Aは処理Sがないので、到着した部品を処理して機械B に送るため、機械Aの仕掛品はあっても１。
実行例：
処理R待ちの行列(黒）、処理P待ちの行列（赤）を重ねて描いたもの。どちらかの優先処理が始まると、どちらか一方の非優先待ち行列がゼロになり、他方の非優 先待ち行列が着実に増えていくことが分かる。
定常条件を満たしているのに、待ち行列長が発散する
左右二つの窓口があり、タイプ１は左から右へ、タイプ２の客は右から左へ進む。
どちらのタイプの客も、２番目の処理が優先される。
窓口の能力より下回る到着量でも、待ち行列長が発散する場合がある。
# タイプ１の客は０(=a1)から１０(=aa)おきに窓口１に到着
# タイプ２の客は５(=a2)から１０おきに窓口２に到着
# 最初のサービスは１(=m1)，２番目のサービス時間は７(=m2)で、２番目のサービス時間が優先権を持つ
# タイプ１の客のサービス１（窓口１）の終了時刻はt11、サービス２（窓口２）の終了時刻はt12
# タイプ２の客のサービス１（窓口２）の終了時刻はt21、サービス２（窓口１）の終了時刻はt22
# タイプ１の客のサービス１の待ち人数（窓口１）をq11、サービス２の待ち人数（窓口２）をq12
# タイプ２の客のサービス１の待ち人数（窓口２）をq21、サービス２の待ち人数（窓口１）をq22
# イベントの生起時点(T1(=T2))で、窓口１，２の客数を記録(Q1,Q2)
m1 = 1; m2 = 7; tmax = 100000;
a1 = 0; a2 = 5; aa = 10 # a1 -> t11(q11) -> t12(q12)
t11 = t12 = t21 = t22 = tmax; # t22(q22) <- t21(q21) <- a2
q11 = q12 = q21 = q22 = 0;
T1 = c(); T2 = c(); Q1 = c(); Q2 = c()
for(i in 1:400) {
ne = min(a1,a2,t11,t12,t21,t22) # 次の事象時刻
if(a1 == ne) { # タイプ１到着
q11 = q11 + 1; # 待ち行列に加える
a1 = a1 + aa; # 次の到着時刻をセット
if(q11 + q22 == 1) t11 = ne + m1; # 誰もいなければサービス開始
}
else if(a2 == ne) { # タイプ2到着（タイプ１の到着と同様）
q21 = q21 + 1;
a2 = a2 + aa;
if(q12 + q21 == 1) t21 = ne + m1; # 誰もいなければサービス開始
}
else if(t11 == ne) { # タイプ１が窓口１から退去、窓口２へ
q11 = q11 - 1; # 待ち行列から削除
q12 = q12 + 1; # 窓口２のタイプ１客が増える
if(q12 + q21 == 1) t12 = ne + m2 # 誰もいなければ、直ちにサービス開始
t11 = tmax; # 事象時刻をクリア（最大値にセット）
if(q22 > 0) t22 = ne + m2 # タイプ２客優先でサービス開始
else if(q11 > 0) t11 = ne + m1; # タイプ２客なし、タイプ１客ありの場合
}
else if(t21 == ne) { # タイプ２が窓口２から退去、窓口１へ
q21 = q21 - 1;
q22 = q22 + 1;
if(q22 + q11 == 1) t22 = ne + m2
t21 = tmax;
if(q12 > 0) t12 = ne + m2
else if(q21 > 0) t21 = ne + m1;
}
else if(t12 == ne) { # タイプ１が窓口２から退去、おしまい
q12 = q12 - 1; # 待ち行列から削除
t12 = tmax; # 事象時刻をクリア（最大値をセット）
if(q12 > 0) t12 = ne + m2 # タイプ１客優先でサービス開始
else if(q21 > 0) t21 = ne + m1; # タイプ１客なし、タイプ２客ありの場合
}
else if(t22 == ne) { # タイプ２が窓口１から退去、おしまい
q22 = q22 - 1;
t22 = tmax
if(q22 > 0) t22 = ne + m2
else if(q11 > 0) t11 = ne + m1
}
Q1 = c(Q1,q11+q22); T1 = c(T1,ne); # 窓口１の統計量
Q2 = c(Q2,q12+q21); T2 = c(T2,ne); # 窓口２の統計量
cat(ne,q11,q12,"//",q21,q22,"//",a1,a2,t11,t12,"//",t21,t22,"\n")
}
qm = max(Q1,Q2)
plot(T1,Q1, ylim=c(0,qm),type="l")
lines(T2,Q2,col=2)
実行結果：
一旦優先客がたまると、非優先客は確実に増え続ける。しかし、優先客がそのうち確実にいなくなるので、非優先客が処理されて次の窓口に進み、優先処理され る。その窓口で非優先客は確実に増え続ける。優先客は確実にいなくなるが、それまでの時間は非優先客として待たされた時間より確実に長い。ということを繰り返 すので、長い方の待ち行列長は無限に増え続ける。
非定常ポアソン過程のサンプルパスは、定常ポアソン過程のサンプルパスからランダムに（時間依存で）「間引く」ことによって生成できる（Lewis Schedler1979）
到着率の時間変化を記述する（最大値が１となる）関数を rate() とする。ピークの到着率（最大値が1の時間帯の到着率）を peakrate とする。
（到着）率 peakrate のポアソン過程の点列を作り、時刻 t で到着があったとした場合、1 - rate(t) の確率で削除する（その到着は無かったものとする）。
到着率が peakrate の斉時ポアソン過程の点列を生成するには、到着間隔が指数分布になることを利用して、パラメータ peakrate の指数乱数を累積すればよい。
時間区間 [0, T] の間の到着過程を生成するには、累積の到着時刻が T になるまで到着間隔を順番に一つずつ生成すれば良い。しかし、ｒはインタープリタ方式なので、なるべく配列処理した方が速くなる。そこで、ランダムな到着間隔を少し多めに 生成して、累積したものが T を超えたものを切り捨てることにする。これで数百倍速くなる。
斉時ポアソン過程の到着時点 t 毎に rate(t) 関数を評価し、1 - rate(t) の確率で、その到着をなかったことにする（間引く）。そのために、一様乱数を生成し、それが rate(t) 以上ならばその到着は無かったものにして棄却する。
残されたものが非斉時ポアソン過程の到着時点列となる。
次に示す nonhPois(T, arate) は、時刻の上限 T、ピーク到着率 arate が与えられたとき、区間 [0, T] の間の非斉時ポアソン過程を生成する関数である。ただし、到着率の変動は、別に rate() 関数によって与えられているものとする。
数値例として、ここではt = 3 から 6 までをピークとした台形状に増減する rate() 関数を仮定した。以下は、使われている関数の説明である。
rexp(n, a) は、パラメータ a（平均 a 分の１）の指数分布に従う乱数を n 個生成する関数である。which(A) は、論理型ベクトル A の TRUE となる要素の番号を返す関数である。which(z<=T) とした場合は、ベクトル z の中で、その値が T 以下となる要素の添字の集合を返す。runif() は、一様乱数を戻り値とする関数である。length() は、ベクトルの要素数を戻り値とする関数である。ifelse(A,b,c) は、条件式 A が正しければ b、さもなければ c を返す関数である。c に ifelse() を使えば、多重比較が出来る。したがって、rate() は、t0, t1, t2, t3 で折れ曲がる折れ線関数を定義している。nonhPois <- function(T, arate) {
# 平均＋3σ個の到着間隔生成
z = cumsum(rexp((T+3*sqrt(T/arate))*arate, arate))
# Tで打ち切り
az = z[which(z<=T)]
# 間引き
az[rate(az) > runif(length(az))]
}
rate <- function(t, t0=0,t1=3,t2=6,t3=8) {
ifelse(t<0, 0, ifelse(t<t1, 1/(t1-t0)*(t-t0), ifelse(t<t2, 1,
ifelse(t<t3, (t-t3)/(t2-t3), 0))))
}
実行例：rate() 関数を図のような折れ線にしたとき、黒、赤の細い線は到着率を 20 とした2通りの結果を示し、太い線は到着率を 100 とした場合の一つの例である。
T = 8 curve(rate(x), 0, T) peakrate = 20 arvt <- nonhPois(T, peakrate) lines(c(0, arvt), (0:length(arvt))/length(arvt))
不特定多数の客の集まるイベント会場の混雑を予測したい。
ある時刻 t に会場に滞在する客の数 L(t) は、累積入場者数 I(t) から累積退去者数 O(t) を引けば求められるので、I(t), O(t) が分かれば良い。
累積入場者数は、気まぐれな不特定多数の客の到着過程と考え、非斉時ポアソン過程を仮定すれば良い。
累積退去者数については、それを直接扱う代わりに、各客の滞在時間をモデル化し、各到着客にそのモデルを適用して退去時刻を求め、それらを順番に並べて累積退 去者数を計算したほうがわかりやすい。
到着過程の非斉時ポアソン過程に対しては、強度関数 λ(t) を与える必要がある。平均的には、開場時から徐々に増加して昼頃ピークを迎え、徐々に減少していく、と仮定するのが自然なので、ここでは、強度関数としてガンマ分布の密度 関数（を打ち切ったもの）を当てはめる。例えば次のように。
滞在時間については、最短でも30分程度、最長は閉館まで、というばらつきを仮定し、平均滞在時間は単調減少すると考え、ベータ分布の密度関数を当てはめる。 ベータ分布のパラメータは例えば (3,3) とする。
非斉時ポアソン過程の到着時点列は、このプログラムを使って生成する。
rate() 関数のガンマ分布は、dgamma() 関数を利用する。
時刻 t に到着した客の滞在時間を与える関数を stay(t) とする。これは rbeta() 関数を利用する。ba,bb はベータ分布のパラメータである。rbeta() 関数は 0 以上 1 以下の乱数を生成するので、それを T-(t+0.5) 倍して 0.5 を足す。T-(t+0.5) がマイナスになる場合は、一律に 0.5 とする。その場合分けのために、ifelse() 関数を使う。if() ... else を使うと、sapply() のようなベクトル計算の時に不具合が生じる。
stay = function(t, ba=2.5, bb=3, T=10) {
ifelse(t > T-0.5, T-t, rbeta(1,ba,bb)* (T - t - 0.5) + 0.5)
}
退去時刻は、到着時刻に滞在時間を足したものなので、到着時点列の各値に対して、stay() 関数値を計算すればよい。到着時点列と関数を与えると、退去時点列ベクトルを生成する関数として、sapply() がある。sapply(A,f) はデータの集合 A の各要素に対して関数 f を適用し、その関数値をまとめたベクトルを返す関数である。
> dp <- sapply(av, stay) + av
このようにして計算した退去時点列は、早いもの順に整列しているわけではないので、累積退去時点列を生成するには、それらを昇順に並べ替える必要がある。昇順 に並べ替えるには sort() 関数を使う。
累積到着客数 I(t)、累積退去客数 O(t) のグラフを描くには、plot(), lines() という関数を利用する。plot(x,y) は散布図を描く関数で、(x,y) は点の x,y 座標ベクトルである。type="l" オプションを使うと、点を折れ線で結ぶ。lwd=2 オプションは線の太さを2倍にすることを指定するオプションである。
lines(x,y) は、現在表示されているグラフの上に重ねて折れ線を描くことが出来る。x,y は折れ曲がる点の (x,y) 座標ベクトルである。col=2 は線を赤で描くことを指定するオプションである。
> plot(av, 1:length(av), type="l", lwd=2, main="累積到着数と累積退去数") > lines(sort(dp), 1:length(dp), col=2, lwd=2)
時刻 t での累積退去数を計算するには、そうやって計算した退去時点列の中から、t 以下の度数を計算すれば良い（length(which(dp<t))）。
時刻 t における滞在客数を計算するには、到着時点列の中の t 以下のものの個数と、退去時点列の中の t 以下のものの個数を計算し、その差を求めれば良い（which() + length()）。
> T <- 10 > arate <- 1000 > av <- nonhPois(T, arate) > dp <- sapply(av, stay) + av > plot(av, 1:length(av), type="l", lwd=2, main="累積到着数と累積退去数") > lines(sort(dp), 1:length(dp), col=2, lwd=2) > tm <- seq(0,T,by=0.01) > Q <- sapply(tm, function(t) length(which(av<t))-length(which(dp<t))) > plot(tm, Q, type="l", lwd=3, main="滞在客数")
逆瀬川(2014)「待ち行列現象のシミュレーション分析」オペレーションズリサーチ 59.4, 198-204。