《待ち行列に対するアルゴリズム的解法》

{まちぎょうれつにたいするあるごりずむてきかいほう}{algorithmic solution of queues}

7.1 背景 従来の待ち行列研究では，多くのモデルに共通する普遍的性質の解明や各種モデルにおける陽な解の導出に主眼が置かれていたが，計算機が身近になると共に，数値計算による待ち行列モデルの分析を意図する研究が活発になってきた[6]．ここで共通する考え方は，(i) 待ち行列の挙動を，系内客数を表す変数と系の挙動を表現するために必要な補助的な変数の組で構成される２変数マルコフ連鎖で定式化し，(ii) その２変数マルコフ連鎖がもつ構造に注目して，マルコフ連鎖の定常分布を計算するための数値計算アルゴリズムを構築する，と言うものである[3, 7, 8]．このようなアプローチをアルゴリズム的解法（algorithmic solution）あるいは行列解析法（matrix-analytic method）と呼ぶ．

7.2 相型分布とMAP 状態 0 が吸収状態，状態 $i \in \M =\{1,2,\ldots,M\}$ が一時的状態である連続時間有限状態吸収マルコフ連鎖 $X(t)$ を考える．このとき $X(t)$ の推移率行列は % \[ \left(\begin{array}{cc} 0 & 0 \\ -\bm{T} \bm{e} & \bm{T} \end{array} \right) \] % の形に書ける．ただし $\bm{e}$ はすべての要素が 1 の列ベクトルである．このマルコフ連鎖の初期状態分布を $(0, \bm{\alpha})$ としたとき，$X(t)$ が吸収されるまでの時間は $(0,\infty)$ 上の確率分布を定め，これを表現$(\bm{\alpha}, \bm{T})$ をもつ相型分布（phase-type distribution）という[1, 4]．以下ではマルコフ連鎖 $X(t)$ の一時的状態を相と呼ぶ． 定義より，表現 $(\bm{\alpha}, \bm{T})$ をもつ相型分布の分布関数は $F(x) = 1 - \bm{\alpha} \exp(\bm{T} x) \bm{e}$ となり，その平均は$\bm{\alpha} (-\bm{T})^{-1}\bm{e}$ である．相の数 $M$ を十分に大きく取ることにより，相型分布は $(0,\infty)$ で定義されたあらゆる確率分布を任意の精度で近似できることが知られており，指数分布，アーラン分布，超指数分布など，待ち行列で頻繁に用いられる確率分布を特別な場合として含む． 表現 $(\bm{\alpha}, \bm{T})$ をもつ相型分布を用いて客の到着間隔の列，すなわち相型再生過程（phase-type renewal process）を生成するには次のようにすればよい．時刻 0 において初期状態分布 $(0,\bm{\alpha})$ に従い相を選ぶ．その後初めて吸収が起こったとき，1 番目の客が到着する．一般に，$n$ 番目（$n=1,2,\ldots$）の客の到着が起こると，直ちに過去の履歴とは独立に初期状態分布 $(0, \bm{\alpha})$ に従い相を選び，その後吸収が起こった時点で $n+1$ 番目の客が到着する．このようにして生成された客の到着間隔は独立で同一な確率変数列となる． 相型分布に従う確率変数の列に対して相関を導入可能にしたものにマルコフ型到着過程（Markovian arrival process：MAP）がある[2, 3]．相型分布を用いた確率変数列の生成では，到着が起こる度に過去の履歴とは独立に初期分布を選んでいたが，MAP では一時的状態 $i \in \M$ から吸収されたとき，確率 $p_{i,j}$ ($j \in \M$) で次の初期状態 $j$ を選ぶ．ただし，全ての $i \in \M$ に対して$\sum_{j \in \M} p_{i,j} =1$ である． 通常，MAP は二つの $M \times M$ 行列 $\bm{C},\bm{D}$ を用いて表現される．ここで行列 $\bm{C}$ は一時的状態間の推移を表し，相型分布の行列 $\bm{T}$ に対応する．一方，行列 $\bm{D}$ の $(i,j)$ 要素は状態 $i$ で吸収が起こり，次の初期状態が $j$ である率 $[-\bm{C}\bm{e}]_i p_{i,j}$ で与えられる．なお，$\bm{C}+\bm{D}$ は規約であり，$\bm{D} \neq \bm{O}$ である．$\bm{\pi}$ を $\bm{\pi} (\bm{C} + \bm{D}) = {\bf 0}$ を満たす確率ベクトルとしたとき，定常な MAP における客の到着間隔分布は $F(x)=1- \bm{\pi} \bm{D} \exp(\bm{C}x) \bm{e} / \bm{\pi} \bm{D} \bm{e}$ となり，その平均は $1/\bm{\pi}\bm{D}\bm{e}$ で与えられる．MAP はあらゆる定常点過程を任意の精度で近似できることが知られており，マルコフ変調ポワソン過程（Markov modulated Poisson process）や独立な相型再生過程の重畳などを特別な場合として含む．特に到着間隔が独立同一な表現 $(\bm{\alpha},\bm{T})$ をもつ相型分布に従う場合は $\bm{C}=\bm{T}$, $\bm{D} =-\bm{T} \bm{e} \bm{\alpha}$ である．

7.3 準出生死滅過程 到着間隔ならびにサービス時間がそれぞれ表現 ($\bm{\alpha}_A, \bm{T}_A$) ならびに表現 ($\bm{\alpha}_B, \bm{T}_B$) をもつ相型分布に従う PH/PH/1 待ち行列を考える．時刻 $t$ における客数を $L(t)$，到着間隔が従う相型分布の相を $S_A(t)$，サービス中の場合はサービス時間が従う相型分布の相を $S_B(t)$ とする．このとき，システムの状態は $(L(t), S_A(t),S_B(t))$ で表される．ただし $L(t)=0$ のときは $S_B(t)=0$ とする． 表現 ($\bm{\alpha}_X, \bm{T}_X$)（$X=A,B$）をもつ相型分布の相の数を $M_X$ としたとき，$(S_A(t),S_B(t))$ の取り得る状態の数は $L(t)=0$ のとき，$M_A$ 通り，$L(t) \geq 1$ のときは $M_A \times M_B$ 通りである．これらに適当な番号を割り当てることで，マルコフ連鎖 $(L(t), S_A(t),S_B(t))$ は２変数マルコフ連鎖 $(L(t), S(t))$ と見なすことができる．さらに，この２変数マルコフ連鎖の状態集合を $L(t)$ の値で分類し，$L(t)=l$ であるような状態からなる部分集合をレベル $l$ と呼ぶ．レベルの変化に注目すると，２変数マルコフ連鎖 $(L(t),S(t))$ は以下のような形の推移率行列 $\bm{Q}$ をもつ連続時間マルコフ連鎖となる． % \[ \bm{Q} = \left(\begin{array}{ccccc} \overline{\bm{Q}}_1 & \overline{\bm{Q}}_0 & \bm{O} & \bm{O} & \cdots \\ \overline{\bm{Q}}_2 & \bm{Q}_1 & \bm{Q}_0 & \bm{O} & \cdots \\ \bm{O} & \bm{Q}_2 & \bm{Q}_1 & \bm{Q}_0 & \cdots \\ \bm{O} & \bm{O} & \bm{Q}_2 & \bm{Q}_1 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array} \right) \] % ここで $\bm{Q}_m$ ならびに $\overline{\bm{Q}}_m$ ($m=1,2,3$)はクロネッカー積 $\otimes$，クロネッカー和 $\oplus$ を用いて次式で与えられる． % \begin{eqnarray*} \bm{Q}_0 = (-\bm{T}_A) \bm{e} \bm{\alpha}_A \otimes \bm{I}_B, & \quad \bm{Q}_1 = \bm{T}_A \oplus \bm{T}_B, \quad & \bm{Q}_2 = \bm{I}_A \otimes (-\bm{T}_B) \bm{e} \bm{\alpha}_B \\ \overline{\bm{Q}}_0 = (-\bm{T}_A) \bm{e} \bm{\alpha}_A \otimes \bm{\alpha}_B, & \quad \overline{\bm{Q}}_1 = \bm{T}_A, \quad & \overline{\bm{Q}}_2 = \bm{I}_A \otimes (-\bm{T}_B)\bm{e} \end{eqnarray*} % ただし $\bm{I}_X$ ($X=A,B$) は $M_X$ 次元の単位行列である． 上記の２変数マルコフ連鎖 $(L(t),S(t))$ は次の性質をもつ．(i) 境界部分を除いて空間的に同質（例えばレベルが1以上（$L(t) \geq 1$）のとき，一回の推移でレベルが一つ増加する率は $\bm{Q}_0$ の要素で与えられる），(ii) 一回の推移でレベルは高々一つしか増減しない．このような性質をもつ２変数連続時間マルコフ連鎖を準出生死滅過程（quasi birth-and-death process:QBD）という[1, 3, 4]． 推移率行列 $\bm{Q}$ をもつ準出生死滅過程において $\bm{Q}_*=\bm{Q}_0+\bm{Q}_1+\bm{Q}_2$ が規約である場合，$\bm{\eta}\bm{Q}_*={\bf0}$, $\bm{\eta}\bm{e}=1$ なる正のベクトル $\bm{\eta}$ が一意に定まり，$\bm{\eta}\bm{Q}_2\bm{e} > \bm{\eta}\bm{Q}_0\bm{e}$ のとき準出生死滅過程 $(L(t),S(t))$ はエルゴード的となる．PH/PH/1 待ち行列の場合，この条件は平均到着間隔 $\bm{\alpha}_A (-\bm{T}_A)^{-1}\bm{e}$ が平均サービス時間 $\bm{\alpha}_B (-\bm{T}_B)^{-1}\bm{e}$ より大きいことと等価である． エルゴード的な準出生死滅過程 $(L(t),S(t))$ の定常分布を $\bm{q} =(\bm{q}_0, \bm{q}_1, \ldots)$ とする．定常分布 $\bm{q}$ は $\bm{q} \bm{Q}={\bf 0}$ を満たすことに注意．このとき， % \begin{equation} \bm{Q}_0 + \bm{R}\bm{Q}_1 + \bm{R}^2 \bm{Q}_2 = \bm{O} \label{rate-matrix-QBD} \end{equation} % を満たす最小の非負行列 $\bm{R}$ を用いて $\bm{q}_k$ ($k=2,3,\ldots$) は % \begin{equation} \bm{q}_k = \bm{q}_1 \bm{R}^{k-1}, \quad k=2,3,\ldots \label{matrix-geometric-QBD} \end{equation} % で与えられる．さらに $\bm{q}_0$, $\bm{q}_1$ は % \[ \bm{q}_0 \overline{\bm{Q}}_1 + \bm{q}_1 \overline{\bm{Q}}_2 = {\bf 0}, \quad \bm{q}_0 \overline{\bm{Q}}_0 + \bm{q}_1 (\bm{Q}_1 + \bm{R}\bm{Q}_2) = {\bf 0} \] % ならびに正規化条件 $\bm{q}_0 + \bm{q}_1(\bm{I}-\bm{R})^{-1}\bm{e}=1$ より求めることができる． 非負行列 $\bm{R}$ は公比行列（rate matrix）と呼ばれており，(2)は行列幾何解（matrix-geometric solution）と呼ばれている．公比行列 $\bm{R}$ は一般には陽に与えられないため，様々な数値計算アルゴリズムが提案されている．最も単純なものは(1)を変形することで得られる % \[ \bm{R}_{n+1} = (\bm{Q}_0 + \bm{R}_n^2 \bm{Q}_2)(-\bm{Q}_1)^{-1}, \quad n=0,1,\ldots \] % に基づいて，$\bm{R}_0=\bm{O}$ を初期値として順次 $\bm{R}_n$ を求めるというものである．このとき $\bm{R}_n$ は要素毎に単調に増加し，公比行列 $\bm{R}$ に収束することが知られている．

7.4 G/M/1型マルコフ連鎖 客の到着間隔 $G$ が独立で同一な分布 $G(x)=\Pr(G \leq x)$ に従い，サービス時間が表現 ($\bm{\alpha}_B, \bm{T}_B$) をもつ相型分布に従う GI/PH/1 待ち行列を考える．$n$ 番目の客の到着直前の系内客数を $X_n$，十 分に多くの客が系内にいるという仮定の下で，$n$ 番目と$n+1$ 番目の客の到着間隔の間にサービス可能な客数を $B_{n+1}$ とすると % \[ X_{n+1} = \max(0, X_n + 1 - B_{n+1}), \quad n=0,1,\ldots \] % が成立する．よって $n$ 番目の客の到着直前におけるサービス時間が従う相型分布の相を $S_n^-$ とすると，$(X_n, S_n^-)$ は下記の遷移確率行列 $\bm{P}_{\rm G/M/1}$ をもつ２変数離散時間マルコフ連鎖となる． % \[ \bm{P}_{\rm G/M/1} = \left(\begin{array}{ccccc} \overline{\bm{B}}_1 & \bm{B}_0 & \bm{O} & \bm{O} & \cdots \\ \overline{\bm{B}}_2 & \bm{B}_1 & \bm{B}_0 & \bm{O} & \cdots \\ \overline{\bm{B}}_3 & \bm{B}_2 & \bm{B}_1 & \bm{B}_0 & \cdots \\ \overline{\bm{B}}_4 & \bm{B}_3 & \bm{B}_2 & \bm{B}_1 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array} \right) \] % ここで $\bm{B}_k$ は % \[ \sum_{k=0}^{\infty} \bm{B}_k z^k = \int_0^{\infty} \exp[(\bm{T}_B + z (-\bm{T}_B)\bm{e}\bm{\alpha}_B) x] dG(x) \] % を満たし，$\overline{\bm{B}}_k$ は $\overline{\bm{B}}_k =\sum_{l=k}^{\infty} \bm{B}_l \bm{e}\bm{\alpha}_B$ で与えられる．準出生死滅過程と同様に，$X_n=l$ であるような状態からなる部分集合をレベル $l$ と呼ぶ． 2変数マルコフ連鎖 $(X_n,S_n^-)$ は(i) 境界部分を除いて空間的に同質（例えばレベルが1以上（$X_n \geq 1$）のとき，一回の遷移でレベルが変化しない確率は $\bm{B}_1$ の要素で与えられる），(ii) 一回の遷移でレベルは高々 一つしか増加しない．このような性質をもつ２変数離散時間マルコフ連鎖をG/M/1型（Markov chain of G/M/1 type）と呼ぶ[3, 4, 7, 8]． 遷移確率行列 $\bm{P}_{\rm G/M/1}$ をもつG/M/1型マルコフ連鎖において，$\bm{B}=\sum_{k=0}^{\infty} \bm{B}_k$ が規約である場合，$\bm{b}\bm{B}=\bm{b}$，$\bm{b}\bm{e}=1$ となるような正のベクトル $\bm{b}$ が一意に定まり，$\sum_{k=2}^{\infty} (k-1)\bm{b}\bm{B}_k\bm{e} > \bm{b}\bm{B}_0\bm{e}$ ならば G/M/1型マルコフ連鎖はエルゴード的である．GI/PH/1 待ち行列の場合，この条件は平均到着間隔 $\E[G]$ が平均サービス時間 $\bm{\alpha}_B(-\bm{T}_B)^{-1}\bm{e}$ より大きいことと等価である． エルゴード的なG/M/1型マルコフ連鎖 $(X_n,S_n^-)$ の定常分布を $\bm{x}=(\bm{x}_0, \bm{x}_1, \ldots)$ とする．定常分布 $\bm{x}$ は $\bm{x} = \bm{x} \bm{P}_{\rm G/M/1}$ を満たすことに注意．このとき， % \begin{equation} \bm{R} = \sum_{k=0}^{\infty} \bm{R}^k \bm{B}_k \label{rate-matrix-GM1} \end{equation} % を満たす最小の非負行列 $\bm{R}$ を用いて $\bm{x}_k$ ($k=1,2,\ldots$) は % \begin{equation} \bm{x}_k = \bm{x}_0 \bm{R}^k, \quad k=1,2,\ldots \label{matrix-geometric-GM1} \end{equation} % で与えられる．さらに $\bm{x}_0$ は % \[ \bm{x}_0 = \bm{x}_0 \sum_{k=0}^{\infty} \bm{R}^k \overline{\bm{B}}_{k+1}, \quad \bm{x}_0 (\bm{I}_B -\bm{R})^{-1} \bm{e} = 1 \] % より求められる． 準出生死滅過程の場合と同様に，非負行列 $\bm{R}$ を公比行列と呼び，(4)を行列幾何解と呼ぶ．公比行列 $\bm{R}$ を求める最も単純な方法は，(3)より得られる % \[ \bm{R}_{n+1} = \sum_{k=0}^{\infty} \bm{R}_n^k \bm{B}_k, \quad n=0,1,\ldots \] % に基づいて，$\bm{R}_0=\bm{O}$ を初期値として順次 $\bm{R}_n$ を求めるというものである．このとき $\bm{R}_n$ は要素毎に単調に増加し，公比行列 $\bm{R}$ に収束することが知られている．

7.5 M/G/1型マルコフ連鎖 客の到着が表現 $(\bm{C},\bm{D})$ をもつ MAP に従い，サービス時間 $H$ が独立で同一な分布 $H(x)=\Pr(H \leq x)$ に従うMAP/G/1 待ち行列を考える．$n$ 番目の客のサービス終了直後の系内客数を $Y_n$，$n$ 番目のサービス中に新たに到着する客数を $A_n$ とすると % \[ Y_{n+1} = \max(0, Y_n -1) + A_{n+1}, \quad n=0,1,\ldots \] % が成立する．よって $n$ 番目のサービス直後における客の到着が従うMAPの相を $S_n^+$ とすると，$(Y_n, S_n^+)$ は下記の遷移確率行列 $\bm{P}_{\rmM/G/1}$ をもつ２変数離散時間マルコフ連鎖となる． % \[ \bm{P}_{\rm M/G/1} = \left(\begin{array}{ccccc} \overline{\bm{A}}_0 & \overline{\bm{A}}_1 & \overline{\bm{A}}_2 & \overline{\bm{A}}_3 & \cdots \\ \bm{A}_0 & \bm{A}_1 & \bm{A}_2 & \bm{A}_3 & \cdots \\ \bm{O} & \bm{A}_0 & \bm{A}_1 & \bm{A}_2 & \cdots \\ \bm{O} & \bm{O} & \bm{A}_0 & \bm{A}_1 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array} \right) \] % ここで $\bm{A}_k$ は % \[ \sum_{k=0}^{\infty} \bm{A}_k z^k = \int_0^{\infty} \exp[(\bm{C} + z \bm{D}) x] dH(x) \] % を満たし，$\overline{\bm{A}}_k$ は $\overline{\bm{A}}_k =(-\bm{C})^{-1} \bm{D} \bm{A}_k$ で与えられる．準出生死滅過程やG/M/1型マルコフ連鎖の場合と同様に，$Y_n=l$ であるような状態からなる部分集合をレベル $l$ と呼ぶ． 2変数マルコフ連鎖 $(Y_n,S_n^+)$ は(i) 境界部分を除いて空間的に同質（例えばレベルが1以上（$Y_n \geq 1$）のとき，一回の遷移でレベルが $k$ 増加する確率は $\bm{A}_{k+1}$ の要素で与えられる），(ii) 一回の遷移でレベ ルは高々一つしか減少しない．このような性質をもつ２変数離散時間マルコフ連鎖をM/G/1型（Markov chain of M/G/1 type）と呼ぶ[3, 5, 7, 8]． 遷移確率行列 $\bm{P}_{\rm M/G/1}$ をもつM/G/1型マルコフ連鎖において，$\bm{A}=\sum_{k=0}^{\infty} \bm{A}_k$ が規約である場合，$\bm{a}\bm{A}=\bm{a}$，$\bm{a}\bm{e}=1$ となるような正のベクトル $\bm{a}$ が一意に定まり，$\sum_{k=1}^{\infty} k \bm{a}\bm{A}_k \bm{e}< 1$ かつ $\sum_{k=1}^{\infty} k \overline{\bm{A}}_k \bm{e}$ の各要素が有限ならば M/G/1型マルコフ連鎖はエルゴード的である．MAP/G/1 待ち行列の場合，この条件は平均到着間隔 $1/\bm{\pi}\bm{D}\bm{e}$ が平均サービス時間 $\E[H]$ より大きいことと等価である． エルゴード的なM/G/1型マルコフ連鎖 $(Y_n,S_n^+)$ の定常分布を $\bm{y}=(\bm{y}_0, \bm{y}_1, \ldots)$ とする．定常分布 $\bm{y}$ は $\bm{y} = \bm{y} \bm{P}_{\rm M/G/1}$ を満たすことに注意．ここで % \begin{equation} \bm{G}=\sum_{k=0}^{\infty} \bm{A}_k \bm{G}^k \label{matrix-G} \end{equation} % を満たす確率行列を $\bm{G}$ とし，$\bm{K}=\sum_{k=0}^{\infty}\overline{\bm{A}}_k \bm{G}^k$ を定義する．行列 $\bm{K}$ は確率行列であることに注意．このとき $\bm{\kappa} = \bm{\kappa} \bm{K}$を満たす確率ベクトル $\bm{\kappa}$ を用いて $\bm{y}_0$ は次式で与えられる． % \[ \bm{y}_0 = \bm{\kappa} / C \] % ただし，定数 $C$ はレベル 0 への平均再帰時間である．さらに $\bm{y}_0$ が与えられたとき，$\bm{y}_k$ ($k=1,2,\ldots$) は次式により順次計算される． % \[ \bm{y}_k = \left( \bm{y}_0 \overline{\bm{C}}_k + \sum_{l=1}^{k-1} \bm{y}_l \bm{C}_{k-l+1} \right) (\bm{I} - \bm{C}_1)^{-1}, \quad k=1,2,\ldots \] % ただし % \[ \bm{C}_k = \sum_{l=k}^{\infty} \bm{A}_l \bm{G}^{l-k}, \quad \overline{\bm{C}}_k = \sum_{l=k}^{\infty} \overline{\bm{A}}_l \bm{G}^{l-k} \] なお，行列 $\bm{G}$ は(5)より得られる再帰式 % \[ \bm{G}_{n+1} = \sum_{k=0}^{\infty} \bm{A}_k \bm{G}_n^k, \quad n=0,1,\ldots \] % に基づいて，$\bm{G}_0=\bm{O}$ を初期値として順次 $\bm{G}_n$ を計算すればよい．このとき $\bm{G}_n$ は要素毎に単調に増加し，確率行列 $\bm{G}$ に収束することが知られている．

---

参考文献

[1] G.\ Latouche and V.\ Ramaswami,{\it Introduction to Matrix Analytic Methods in Stochastic Modeling},SIAM, 1999.

[2] D.\ M.\ Lucantoni, K.\ S.\ Meier-Hellstern and M.\ F.\ Neuts,``A single-server queue with server vacations and a class of nonrenewal arrival processes,{\it Advances in Applied Probability}, {\bf 22} (1990), 676--705.

[3] 牧本直樹，『待ち行列アルゴリズム ---行列解析アプローチ---』，朝倉書店, 2001.

[4] M.\ F.\ Neuts,{\it Matrix-Geometric Solutions in Stochastic Models, An Algorithmic Approach}, Johns Hopkins University Press, 1981.

[5] M.\ F.\ Neuts,{\it Structured Stochastic Matrices of M/G/1 Type and Their Applications}, Dekker, 1989.

[6] 高橋幸雄，「待ち行列研究の変遷」, 『オペレーションズ・リサーチ』，{\bf 43} (1998), 495--499.

[7] 高橋幸雄，牧本直樹，「相型分布と行列解析法」, 『オペレーションズ・リサーチ』，{\bf 43} (1998), 618--623.

[8] 滝根哲哉，「構造化されたマルコフ連鎖と待ち行列」, 『システム／制御／情報』，{\bf 43} (1999), 135--140.