論文

バーストトラヒックを入力とするネットワークにおける平均系内時間
の近似解析
米倉

誠†

中川 健治†† a)

恩田 和幸†††

斎藤

茂††††

Approximation Analysis for the Average Sojourn Time in a Network with Bursty
Input Traﬃc
Makoto YONEKURA† , Kenji NAKAGAWA†† a) , Kazuyuki ONDA††† ,and Shigeru SAITO††††

あらまし 待ち行列からの退去過程は一般に到着過程よりも複雑になるため，他ノードからの退去過程が到着過程となるネットワーク内の各ノードの解析は一般に困難である．従来は解析を容易にするため，ネットワーク
を Jackson 網と仮定して，すなわち，到着をポアソン過程，サービスを指数サービスと仮定してネットワークを解析・設計していた．しかし，ポアソン過程では実際のパケットネットワークにおけるパケットトラヒックを表現するには不十分である．本研究では，バースト性を表現できる MMPP（Markov Modulated Poisson Process）
を到着過程とし，指数サーバ多段系列の各ノードでの平均系内時間の近似値を得る方法を提案する．平均系内時間の計算のために，各ノードからの退去過程を，MMPP を含みより一般的なトラヒックモデルである MAP（Markovian Arrival Process）によって表現し，更にそれを MMPP で近似する．これによって，全ノードにおいて同じアルゴリズムで平均系内時間を計算できる．提案法による平均系内時間の近似値とシミュレーション値を比較して，提案法による近似値が安全側で精度良い値であることを示す．
キーワード

待ち行列，平均系内時間，MMPP，MAP，退去過程

1. ま え が き

ビスと仮定して解析すれば，ネットワークの特性値を通信ネットワーク内のノードは待ち行列でモデル化

信ネットワーク内を流れるトラヒックを表現するのに得ることは容易であるが，ポアソン過程では実際に通され解析される．ネットワーク内では，あるノードかは不十分である．

らの退去過程が他ノードへの到着過程となるが，一般

本研究では，到着過程を 2 状態 MMPP（Markovに待ち行列からの退去過程は到着過程よりも複雑にな

Modulated Poisson Process）として，各ノードでのるために，ネットワーク内の一つのノードの解析でさサービスを指数サービスとする直列のネットワークをえ困難になる．ネットワークを Jackson 網と仮定して，考え，各ノードでの平均系内時間の近似値を得ること

すなわち，到着をポアソン過程，サービスを指数サーを目的とする．対象とするモデルを図 1 に示す．ノード n (n = 1, . . . , L) はサイズ無限大のバッファとサー†

松下電送システム株式会社，東京都
Matsushita Graphic Communication Systems, Inc., Meguroku, Tokyo, 153–8687 Japan

††

長岡技術科学大学，長岡市
Department of Electrical Engineering, Nagaoka University†††

本研究は，パケット通信網のネットワークリソースの最適配分法の設計を目的として行われている．パof Technology, Nagaoka-shi, 940–2188 Japan

ケットトラヒックのバースト性を考慮して，バースト

NTT アドバンステクノロジ株式会社，武蔵野市

性をもつ基本的なモデルである MMPP をトラヒック

NTT Advance Technology Co. Ltd., 2–4–15 Nakamachi,モデルとした．また，サービスは必ずしも指数サービ

Musashino-shi, 180–0006 Japan††††

ビス率 µn の指数サーバからなる．

ATR 適応コミュニケーション研究所，京都府

スでよく表現できるとは限らないが，パケット長は一ATR Adaptive Communications Research Labs., Seika-cho,

定でもない．そこで，近似的なサービスモデルとして

Souraku-gun, Kyoto-fu, 619–0288 Japana) E-mail: nakagawa@nagaokaut.ac.jp

電子情報通信学会論文誌 B Vol. J85–B

解析的に取り扱いやすい指数サービスを仮定した．直No. 7 pp. 1021–1030 2002 年 7 月

1021

電子情報通信学会論文誌 2002/7 Vol. J85–B No. 7

起こる．本研究で到着過程として用いる 2 状態 MMPP
は 4 パラメータをもち，そのパラメータは，状態 0，1
における到着率 λ0 ，λ1 と状態 0，1 から他方の状態への遷移率 q0 ，q1 である．
図 1 ネットワークモデル
Fig. 1 Network model.

2. 2 MAP
MAP（Markovian Arrival Process）は MMPP と同様に，連続時間マルコフ連鎖の各状態がそれぞれ一列ネットワークを仮定する理由は，1 本 1 本のコネク

つのポアソン過程に対応している．MAP は状態遷移

ションの end-to-end の遅延時間を調べることを目的に到着を伴うことを許すモデルである [2]．到着を伴としているためである．
は MAP（Markovian Arrival Process）の特別な場合

わない遷移率行列 C と，到着を伴う遷移率行列 D が
MAP のパラメータとなる．これら C と D を MAPのパラメータ行列という．MMPP では到着を伴う状とみなすことができ，MMPP を MAP として表現す

態遷移は考えないが，到着時には状態が変化しない

ることによってノード 1 を 2 状態 MAP を入力とすると考える，あるいは，同じ状態に遷移すると考えればMAP/M/1 として解析することができる．MAP/M/1

MMPP を MAP とみなすことができる．こう考えて2 状態 MMPP を MAP として表現すると，パラメー
タ行列 C ，D は次のようになる．

ノード 1 への入力は 2 状態 MMPP である．MMPPの平均系内時間の解析は文献 [1] に示されている．その方法は，MAP のパラメータ行列を直接用いた行列の繰返し計算を主たる計算としており，アルゴリズム
が簡単である．しかし，我々の問題にこの解析法を直C=接適用するのは計算量の点から困難である．なぜならと，ノード n (n = 1, . . . , L) への入力となる MAP のn−1

状態数は 2(K + 1)

この場合には，ノード L = 10 での状態数は 2 × 219となってしまう．
そこで，本研究では，各ノードからの退去過程を
ノード 1 への到着過程と同じモデル，すなわち 2 状態

q1

−q1 − λ1

λ0
0

D=となってしまうからである．後に述べるように，K = 20, L = 10 の場合を扱うので，

q0



ば，各ノードのバッファサイズを K として近似する

−q0 − λ0



(1)

0
λ1

(2)

3. MAP/M/1 の解析
本研究で各ノードにおける平均系内時間の計算のために利用する MAP/M/1 の解析を示す（文献 [1]参照）．

すべてのノードにおいて同じ計算量で平均系内時間

3. 1 MAP/M/1 の定常分布
パラメータ行列 C ，D をもつ MAP を到着過程を計算することができる．そして，実際にシミュレーとし，サービス率 µ をもつ指数サービスの待ち行列

ション値と比較して，提案法による近似値が安全側のMAP/M/1 を考える．MAP の状態数を d とする．C ，
D はともに d × d 行列である．MAP/M/1 のキュー
長を N とし，MAP の状態を S とすると，組 (N, S)

MMPP で近似し次ノードを解析する．それによって，精度の良い値であることを示す．

2. で本論文で扱うトラヒックモデルの説明を行う．
3.，4. では文献 [1] による MAP/M/1 の平均系内時間の計算アルゴリズムを紹介する．5. で提案する近似法を説明し，6. で数値例を示し，7. でそれに対する考察を行う．

2. 到着過程について
2. 1 MMPP
MMPP（Markov Modulated Poisson Process）は連続時間マルコフ連鎖の各状態がそれぞれ一つのポア
ソン過程に対応し，そのポアソン過程に従って到着が1022がなす連続時間マルコフ連鎖の遷移率行列 Q は



C
 A


0
Q=

 0

..
.

D
B
A
0
..
.

0
D
B
A
..
.

0
0
D
B
..
.



···
··· 


··· 

··· 

..
.

(3)

によって与えられる．ただし，A = µI ，I は d × d 単位行列，B = C − A である．

MAP/M/1 の定常分布 x は xQ = 0，x1 = 1 の解

論文／バーストトラヒックを入力とするネットワークにおける平均系内時間の近似解析である．ここで，1 はすべての要素が 1 である列ベクを伴わない遷移率行列となり，式 (10) の退去を伴う

トルを表す．定常分布 x は遷移率行列 V が次ノードへの到着を伴う遷移率行列

x = [x0 , x0 R, x0 R2 , . . . ]

(4)

となる．
そこで，3. 1 の MAP/M/1 の定常分布の解析において，C = U ，D = V として式 (8) によって次ノーという表示をもつ．ここで R は，

ドでの平均系内時間を求めたいのであるが，この場合，µR2 + RB + D = 0(5)

繰返し計算を実行することができない．状態数が無限満たす最小非負の行列である．また，x0 は

x0 = π(I − R)

大となるのは，バッファサイズを無限大としているた(6)と与えられる．ここで，π は入力の MAP の定常分布である．文献 [1] では，R を数値的に得るために，式(5) を変形して得られる漸化式
Rk+1 = (µRk2 + D)(−B)−1 ,C ，D がサイズ無限大の行列であるため漸化式 (7) のめである．
（注意） 状態数が無限の場合，式 (9)，(10) の遷移率行列をもつ Markov 連鎖を一般には MAP といわないが，本論文では簡単のため状態数無限の MAP という．

R0 = 0R = limk→∞ Rk(7)

4. 有限バッファによる近似
ノード 1 のバッファサイズを有限値 K で近似すれ
ば，ノード 2 への到着過程は 2(K + 1) 状態の MAPを用いる．

3. 2 MAP/M/1 の平均系内時間で表現できる．その MAP の到着を伴わない遷移率行

MAP/M/1 の平均系内時間 W̄ は，MAP/M/1 の列 Û と，到着を伴う遷移率行列 V̂ は下記のように与平均キュー長 N̄ と MAP の平均到着率 λ を用いて，次のリトルの公式 (8) で計算できる．

W̄ = N̄ /λ(8)

ここで，N̄ = π(I − R)−1

R1，λ = πD1 と与えられる．

3. 3 MAP/M/1 からの退去過程
MAP/M/1 のキュー長 N と状態 S の組 (N, S)

えられる．



C
 0




Û = 


を考え，組としての状態遷移を考えることによってB̂ = B + D

MAP/M/1 からの退去過程の遷移率行列を求める．客の退去に注目する．キュー長 N が一つ減少するとき

0
 A


V̂ =  0

が退去であり，N が変化しないかまたは増加のときは退去でない．したがって，式 (3) の遷移率行列 Q を
退去を伴わない遷移率行列 U と退去を伴う遷移率行
列 V とに分けると以下のようになる．



C

 0
U =
 0

..
.
V =Q−U

D
B
0
..
.

0
D
B
..
.

0
0
D
..
.

0
D
..
.
0


0
..

.

B
0

..

D
B
0







.







0 

D 

(11)

B̂


0
A
..
.

0
..
.

..

0

A









.

(12)

0

Û ，V̂ は 2(K + 1) × 2(K + 1) 行列である．行列サイズが有限なのでノード 2 を MAP/M/1 として 3. 1 の



...

... 

... 

..
.

D
B
..
.

(9)

解析法が使える．以降のノードも同様に，ノード n ま
でのバッファサイズを有限で近似すれば，ノード n + 1
を MAP/M/1 として解析できる．しかし，このよう(10)な近似では後のノードほど到着過程となる MAP の状
態数が大きくなり，ノード L における MAP の状態

退去過程が次ノードへの到着過程となるので，式 (9)数は 2(K + 1)L−1 となってしまう．したがって，対の退去を伴わない遷移率行列 U が次ノードへの到着

象モデルのノード数が多いときに全ノードでこの解析
1023

電子情報通信学会論文誌 2002/7 Vol. J85–B No. 7

法を実行するのは計算量の点で困難である．
そこで，退去過程を簡単なモデルで近似することを考える．特に各ノードからの退去過程をノード 1 へのCn ，Dn は 2 × 2 行列である．MAP (Cn , Dn ) はノー
ド n + 1 への入力となる．以下，同様のことを繰り返し行う．

到着過程と同じモデルである 2 状態 MMPP で近似す

ここで特に，MAP (Ûn , V̂n ) の遷移率行列 Ûn と V̂nれば，各ノードでの計算アルゴリズムが同じになる．

を示しておく．これらはともに 2(K + 1) × 2(K + 1)

5. 退去過程の 2 状態 MMPP による近似（提案法）
上記で述べたように，各ノードのバッファサイズを有限値 K で近似するだけでは，各ノードへの到着過程となる MAP の状態数がノード番号に対して K + 1のべき乗で増加し，解析が困難になる．そこで本研究
では，各ノードからの退去過程を 2 状態 MMPP で近行列である．

Cn−1 Dn−10 0 Bn−1 Dn−1 0





..
..
..
..


.
.
.
.


Ûn =

00 Bn−1 Dn−1

0 Bn−1 Dn−1 
0 B̂n−1

似する方法を提案する．それによってすべてのノード
で同一のタイプの待ち行列を解析すればよいことに
なる．

(13)
Bn−1 = Cn−1 − µn I, B̂n−1 = Bn−1 + Dn−1



0
 µ I
 n

ノード n − 1 からの退去過程を 2 状態 MMPP で近似し，その 2 状態 MMPP を MAP とみなしたときのパラメータ行列を Cn−1 ，Dn−1 とする（図 2 参照）．

Cn−1 ，Dn−1 はともに 2 × 2 行列である．この MAP
を MAP (Cn−1 , Dn−1 ) と表す．ノード n − 1 からの退去はノード n への入力となる．そこで，ノード n での待ち行列は MAP (Cn−1 , Dn−1 )/M/1 となる．3. 3
で述べたように，MAP/M/1 からの退去過程は式 (9)，(10) の U ，V をパラメータ行列とする状態数無限のMAP となる．その状態数無限の MAP のパラメータ
行列を Un ，Vn と表し，その MAP を MAP (Un , Vn )と表す．そして，MAP (Un , Vn ) を 4. の有限バッファによる近似法を用いて近似して得られる 2(K + 1) 状
態 MAP を MAP (Ûn , V̂n ) と表す．
更に，MAP (Ûn , V̂n ) を以下に示す方法で 2 状態

MMPP で近似し，それを MAP (Cn , Dn ) と表す．

V̂n = 




0










0
µn I
..
.

0
..
.

..

0

µn I

.

(14)

0

5. 1 近似のための条件
さ て ，Bean [3] に よ る と ，一 般 に 二 つ の MAP，
MAP (Un , Vn ) と MAP (Cn , Dn ) が確率過程として一致するための必要十分条件は，任意の自然数 k と任意の実数 t1 , . . . , tk ∈ [0, ∞) に対して次式が成り立つことである．

xeUn t1 Vn eUn t2 · · · Vn eUn tk Vn 1
xeUn t1 Vn eUn t2 · · · Vn eUn tk 1=

ψeCn t1 Dn eCn t2 · · · Dn eCn tk Dn 1
ψeCn t1 Dn eCn t2 · · · Dn eCn tk 1

(15)

ここで，x は MAP (Un , Vn ) の定常分布であり，ψ は

MAP (Cn , Dn ) の定常分布である．式 (15) の左辺はMAP (Un , Vn ) における過去の履歴を条件とする時刻
T = t1 + · · · + tk での到着率を表し，式 (15) の右辺
は MAP (Cn , Dn ) における同様の値を表す．
式 (15) が任意の k と任意の t1 , . . . , tk ∈ [0, ∞) に対して成り立てば MAP (Un , Vn ) を MAP (Cn , Dn )で完全に表現できたことになるのだが，行列 Un ，Vnはサイズ無限大で行列 Cn ，Dn は 2 × 2 なのでそれは期待できない．そこで，式 (15) の k や t1 , . . . , tk
図 2 退去過程の 2 状態 MMPP による近似
Fig. 2 Approximation of the departure process by
2-state MMPP.

1024

を特別な値に限定したものを Cn ，Dn が満たすべき条件とする．

論文／バーストトラヒックを入力とするネットワークにおける平均系内時間の近似解析

以上の考察に基づき，Cn ，Dn を決定するための条件として下記の条件 1〜4 を提案する．

様に式 (17) より次のように計算できる．

xVn2 1 = µ2 πR2 1

条件 1：xVn 1 = ψDn 1
2
条件 2：xVn2 1 = ψDn1(18)

5. 4 条件 3 について
条件 3 を考えるにあたって，まず，式 (13) において

条件 3：g(Ûn ) = g(Cn )条件 4：W̄1 = W̄2
ここで，Ûn は十分大きな K により式 (13) で与えg(Ûn ) = g(Cn−1 )(19)られる行列であり，g(Ûn )，g(Cn ) は，それぞれ，Ûn ，であること，すなわち，Ûn の最大固有値と Cn−1 のCn の最大固有値を表す．また，W̄1 は MAP (Ûn , V̂n )を到着過程とし，次ノードのサービス率 µn+1 をもつ

最大固有値が等しいことを示す．式 (13) より，Ûn の

指数サーバをサーバとする待ち行列の平均系内時間を表し，W̄2 は MAP (Cn , Dn ) を到着とする同様の待提案した条件 1〜4 について詳細に検討する．
5. 2 条件 1 について
条件 1 は式 (15) で k = 1，t1 = 0 として得られる．
このとき条件 1 の左辺は，MAP (Un , Vn ) の平均到着率を表し，右辺は MAP (Cn , Dn ) の平均到着率を表す．条件 1 はこれらが等しいことを要請する．
条件 1 の左辺はサイズ無限大の行列 Vn を含んでいるが，3. 1 における 2 × 2 行列 R を用いて次式 (16)で計算できる．

xVn 1 = µπR1

(16)xVnm 1= xVn Vnm−1 1= [π(I − R), π(I − R)R, · · · ]
0

..
.

0
µI
..
.


 µI
×
 0


0

...



..

= µ [π(I − R)R , π(I − R)R

m+1, · · · ]1

∞

= µm

= µ πR 1

ま ず，Bn−1 = Cn−1 − µn I よ り g(Bn−1 ) =g(Cn−1 ) − µn であるから g(Bn−1 ) < g(Cn−1 ) で
ある．

(21)

となる．式 (21) の左辺を θ の 2 次関数とみて f (θ) とおく．f (−λ) ≤ 0 であることがわかるから，式 (21)の大きいほうの解，すなわち θ = g(Cn−1 ) は −λ ≤ θを満たす．待ち行列の安定条件 λ < µn より，結局，π(I − R)Rm+k 1

k=0
m
m

下のようにしてわかる．

=0

.

m

I は 2 × 2 単位行列である．式 (20) より，Ûn の最大固有値 g(Ûn ) は Cn−1 ，Bn−1 ，B̂n−1 の最大固有値
g(Cn−1 )，g(Bn−1 )，g(B̂n−1 ) のいずれかである．実際，これらのなかで g(Cn−1 ) が最大であることが以θ2 + (λ0 + λ1 + q0 + q1 )θ + λ0 λ1 + λ(q0 + q1 )= µ[π(I − R)R, π(I − R)R2 , · · · ]Vnm−1 1
mとなる．ここで，I  は 2(K + 1) × 2(K + 1) 単位行列，

λ = (λ0 q1 + λ1 q0 )/(q0 + q1 ) であることに注意して，


0 ...  m−1 1
0 . . .  Vn

..
.

(20)

次 に ，g(B̂n−1 ) < g(Cn−1 ) を 示 す．B̂n−1 =Cn−1 + Dn−1 − µn I よ り g(B̂n−1 ) = −µn である．θ = g(Cn−1 ) とおいて −µn < θ を示す．Cn−1
を式 (1) のように表すと，固有値 θ が満たすべき方程式は，MAP (Cn−1 , Dn−1 ) の平均到着率 λ がなぜならば，m ≥ 1 に対して

0

|Ûn − θI  |= |Cn−1 − θI||Bn−1 − θI|K−1 |B̂n−1 − θI|

ち行列の平均系内時間を表す．



固有方程式は，(17)g(B̂n−1 ) = −µn < −λ ≤ g(Cn−1 ) を得る．
以上から，g(Ûn ) = g(Cn−1 ) であることがわかった．いいかえると，到着を伴わない遷移率行列の最大となるからである．

固有値は各ノードの入力と出力で不変となる．

5. 3 条件 2 について
条件 2 は式 (15) で k = 2，t1 = t2 = 0 とし，条件1 を適用して得られる．条件 2 の左辺は，条件 1 と同

Vn を Ûn ，V̂n で置き換え，k = 1 の場合を考える．

そこで，条件 3 を得るために，式 (15) の左辺の Un ，
すなわち，
1025

電子情報通信学会論文誌 2002/7 Vol. J85–B No. 7

ψeCn t Dn 1x̂eÛn t V̂n 1=
ψeCn t 1x̂eÛn t 1

ることを要請する．しかし，実際の計算においては，(22)である．ここで，x̂ は MAP (Ûn , V̂n ) の定常分布である．また，簡単のため t = t1 とおいた．

(Ûn + V̂n )1 = 0，(Cn + Dn )1 = 0 より式 (22) はψex̂e Ûn 1
Cn 1=
ψeCn t 1x̂eÛn t 1
Ûn t

Cn t

(23)

となる．式 (23) において，左辺の Ûn の中に含まれる Cn−1 及び右辺の Cn がともに既約であるとする．
ここで，式 (1) の形の行列 Cn−1 が既約であるとは，ある正数 α > 0 に対して Cn−1 + αI が既約であることである．Cn についても同様．
式 (23) の右辺について，付録の定理より

lim

ψe
Cn 1
= g(Cn )ψeCn t 1(24)を得る．次に，式 (23) の左辺を考える．式 (19) 及びCn−1 が既約であるとの仮定から，g(Ûn ) は Ûn の重
複度 1 の最大固有値である．そして，g(Cn−1 ) に対す
る正の固有列ベクトル x1 > 0 をとって Ûn の固有ベ
クトル



(25)

x̂eÛn t Ûn 1= g(Ûn )t→∞ x̂eÛn t 1

個のパラメータ λ1 , q0 , q1 を簡単な式で表すことができる．したがって，MAP (Cn , Dn ) による平均系内時
間 W̄2 は 1 変数 λ0 の関数となる．そこで，1 変数のニュートン法によって |W̄2 − W̄1 | をなるべく小さくなるようにし，その最小値を与えるパラメータを求めるパラメータとする．

6. 提案する近似値とシミュレーション値の比較のシミュレーション値（sim）と，近似バッファサイズ

K = 20 として退去過程を 2 状態 MMPP である MAP(Cn , Dn ) で近似して計算した近似値（proposed），及び M/M/1 の平均系内時間（M/M/1）を表 1，表 2の六つの場合（case 1〜6）において比較する．それぞれの場合の，ノード 1 への到着過程である MMPP のパラメータ (λ0 , λ1 , q0 , q1 ) を表 1 に示し，各ノードの結果，図 4 は case 2 の結果，· · · ，図 8 は case 6 の結果を表している．提案する近似値（proposed）の計算では，漸化式 (7) の収束及び条件 1〜4 を満たしているとの判断は，ともに誤差率 10−5 以下とした．従を考えれば，定理の証明と同様の議論によって

lim

条件 1〜3 より MAP (Cn , Dn ) の 4 個のパラメー
タのうちの 1 個，例えば，到着率 λ0 を用いて他の 3n のサービス率 µn を表 2 に示している．
比較した結果を図 3〜図 8 に示す．図 3 は case 1



x1


 0 
 . 
 . 
 . 
0

ことはできない．そこで，以下のような近似を行う．

L = 10 としたときのノード n での平均系内時間 Wn

Cn t

t→∞

状態数の違いから必ずしも W̄1 と W̄2 を一致させる

来のように到着過程をポアソン過程とみなすと，平均系内時間は図中の M/M/1 となり，シミュレーション(26)

値に対して危険側であることがわかる．

を得る．式 (23)，(24)，(26) より g(Ûn ) = g(Cn ) となる．これが条件 3 である．

表 1 入力の 2 状態 MMPP のパラメータ値
Table 1 Parameter values of input 2-state MMPP.

条件 3 及び式 (19) より

g(Cn ) = g(Cn−1 )casecase 1, 2, 3
case 4, 5, 6

(27)

λ0
1.0
1.0

λ1
0.0
0.0

q0
0.1
0.1

q1
0.01
0.1

となり，2 状態 MMPP の到着を伴わない遷移率行列の最大固有値が，ノード 1 からノード L まで不変に保たれる．

表 2 各ノードのサービス率
Table 2 Service rate of each node.

5. 5 条件 4 について
本研究では各ノードでの平均系内時間を評価しているので，MAP (Ûn , V̂n ) のよる平均系内時間 W̄1 とMAP (Cn , Dn ) による平均系内時間 W̄2 が等しくな
1026case 1, 4case 2, 5case 3, 6

µ1
1
1
1

µ2
1
3
2

µ3
1
1
3

µ4
1
3
4

µ5
1
1
5

µ6
1
3
5

µ7
1
1
4

µ8
1
3
3

µ9
1
1
2

µ10
1
3
1

論文／バーストトラヒックを入力とするネットワークにおける平均系内時間の近似解析

Fig. 3図 3 各ノードでの平均系内時間（case 1）
Average sojourn time at each node (case 1).

Fig. 6図 6 各ノードでの平均系内時間（case 4）
Average sojourn time at each node (case 4).

Fig. 4図 4 各ノードでの平均系内時間（case 2）
Average sojourn time at each node (case 2).

Fig. 7図 7 各ノードでの平均系内時間（case 5）
Average sojourn time at each node (case 5).

Fig. 5図 5 各ノードでの平均系内時間（case 3）
Average sojourn time at each node (case 3).

Fig. 8図 8 各ノードでの平均系内時間（case 6）
Average sojourn time at each node (case 6).

7. 考

察

提案法ではノードを通過するたびに近似を実行するので，後のノードになるほど近似精度が悪くなるが，こ7. 1 近似精度についてのような場合でもシミュレーション値に対して安全側提案法による近似値は，多くの場合において平均系であることが多くの数値例から確認できている．

内時間のシミュレーション値に近い値を与えている．

筆者らは，退去過程の近似のために四つの条件（条1027電子情報通信学会論文誌 2002/7 Vol. J85–B No. 7

件 1〜4）を一致させることを提案した．特に条件 1の平均到着率と条件 3 の到着を伴わない遷移率行列の最大固有値は，ノードを通過する前と後で常に保存

表 3 近似バッファサイズ K に対する提案法の計算時間
T とノード 10 の平均系内時間の近似値 W̄
Table 3 Computation time T and the averagesojourn time W̄ for the buﬀer size K.

される量である．このことから，条件 1 と条件 3 を

K
10
20
30
40
50

要請することは自然なことである．また，条件 4 では，各ノードからの退去過程とそれを近似する 2 状態

MMPP の次ノードのサービス率よる平均系内時間が一致することを要請している．評価すべき平均系内時

T
0.02 (s)0.17 (s)0.66 (s)1.75 (s)3.73 (s)

W̄
2.21
2.27
2.27
2.27
2.27

間を一致するようにしているので，精度の良い特性が得られたと考えられる．条件 1〜3 で退去過程の統計提案法では 0.17 秒，シミュレーションでは 11.9 秒か的な性質を取り入れ，条件 4 で平均系内時間に対するかった．

性質を取り入れ，これらを組み合わせたことによって良い結果が得られたといえる．
筆者らは，最終的に条件 1〜4 に到達するまでにいくつかの条件の組を試みた．例えば，条件 1〜3 に加えて，条件

xVn3 1 = ψDn3 1(28)

すなわち，式 (15) において k = 3, t1 = t2 = t3 = 0

また，同じパラメータ条件で，近似バッファサイズ

K に対する提案法の計算時間 T 及び 10 番目のノードの平均系内時間 W̄ を調べ，表 3 の結果を得た．シ
ミュレーションの実行時間は K に依存しない．
この表より，T
K 3.2 であることがわかる．また，W̄ の値を見ると K = 20 程度で十分収束しているとみなせるので，提案法の計算時間はシミュレーションと比較して高速であるといえる．

K の決め方についてとおいた条件の四つの条件によって到着過程を近似す

7. 3

ることを考えた．しかし，この場合，得られた平均系

バッファサイズを近似する K の値が小さいと，退内時間の近似値はすべて危険側の値となってしまった．

去過程を近似する MMPP のパラメータ値を決定でき

あるいはまた，退去過程を近似する MMPP のパラない場合がある．K の値が小さいと平均系内時間の近メータを遷移率 q0 ，q1 をノード 1 への到着過程の遷似値が M/M/1 の平均系内時間よりも小さくなる．条

移率と等しくし，到着率 λ0 , λ1 を条件 1，4 を満たす

件 1 で平均到着率を一致させているので，到着率の等ように決定するという近似法でも，ノード 3 以降でしいどんな MMPP を用いてもその平均系内時間 W̄2は危険側となってしまった．このことから，条件 4 をは M/M/1 の平均系内時間よりも大きい．この場合に満たすだけでは精度良い近似とはならないことがわか

上記のパラメータ決定法を実行すると到着率や遷移率

る．筆者らが提案する条件 1〜4 を組み合わせることが負の値になるなど不都合が生じる．K の値を大きくで，条件 4 を満たす MMPP のうちで，部分的に退去することによってこの問題を回避することができる場

過程と同じ性質をもち，同時に部分的に到着過程と同合もあるが，そうすると計算量が多くなるので，K のじ性質をもつ MMPP を選ぶことで精度良く近似でき，値を決める基準が問題になる．

安全側になると考える．

7. 2 実行時間の比較について

8. む すび

シミュレーション値及び提案法による近似値を得る

2 状態 MMPP を到着過程とする指数サーバ多段系ために必要な計算時間の比較を行う．case 1 の 10 ノー

列の各ノードでの平均系内時間の近似値の計算方法を

ドの平均系内時間を得る計算時間を比較する．K = 20

提案した．提案法による平均系内時間の近似値がシュとする．ほぼ同等な条件で比較するために，提案法

ミレーション値に対して安全側で精度が良いことを確における条件 4 のニュートン法の収束条件の誤差の 2

認した．本論文には示してないが，筆者らが行った 60

乗を 10−6 とし，一方，シミュレーションでは 10 番

種類の例においても同様のことを確認した．

目のノードの平均系内時間に対する推定値の分散が各ノードのバッファサイズを単に有限値 K で近似

2.8 × 10−5 となるようにシミュレーション時間を決めすると，ノード L への到着過程となる MAP の状態

た．用いた計算機の性能は，Pentium III，1 GHz，メ

数が 2(K + 1)L−1 となるので，MAP/M/1 の解析法

モリ 256 MByte，HDD 80 GByte である．その結果，をそのまま本論文の問題に適用するのは計算量の点で1028

論文／バーストトラヒックを入力とするネットワークにおける平均系内時間の近似解析困難である．
提 案 法に より各 ノー ドか らの 退去 過程を 2 状 態

MMPP で近似することで，ノード 2 以降に対して同
じ計算アルゴリズムで各ノードを MMPP/M/1/K で近似して解析することができた．提案法では，退去過程を近似する 2 状態 MMPP のパラメータ値を決定するための四つの条件を与えた．バッファサイズ K の値が十分に大きい場合には平均系内時間の近似値を精度良く計算できた．K の値が不十分で退去過程を近似する MMPP のパラメータ値を決定できない場合もあったので，K の値の決め方が今後の課題となる．
また，一般のネットワークではトラヒックの合流があるが，このような場合に本研究を拡張することも今後の課題である．合流がある場合には，注目するコネクション以外を背景トラヒックと考えて，それらをまとめてポアソントラヒックとして扱うことを検討中である．
謝辞

成 11〜13 年度の助成研究費によって行われた．また，シミュレーションと近似値の計算時間の比較に関して御協力頂いた長岡技術科学大学修士 2 年生の坂本将清氏に感謝致します．
文

献M.F. Neuts, Matrix-geometric Solutions in Stochastic Models, Johns Hopkins Univ. Press, Baltimore,
MD, 1981.

[2]

D.M. Lucantoni, K.S. Meier-Hellstern, and M.F.
Neuts, “A Single Server Queue with Server Vacations
and a Class of Non Renewal Arrival,” Adv. Appl.
Prob., vol.22, pp.675–705, 1990.

[3]

yeCt Cz=θt→∞ yeCt z
lim

N. Bean, D. Green, and P. Taylor, “The output（証明） α > 0 を十分大きな正数とすると，C + αI
は既約な非負行列となる．Perron-Frobenius の定理より，C + αI は重複度 1 の正の固有値 θα をもち，θα以外の固有値 ηα,2 , . . . , ηα,d に対して

|ηα,i | ≤ θα(i = 2, . . . , d)

えたものに等しいので，定理の θ は θ = θα − α であり，C の他の固有値 ηi は ηi = ηα,i − α (i = 2, . . . , d)である．式 (A·2) より，(i = 2, . . . , d)ﬁnite-demensional ﬁlters for Markov-modulated trafﬁc,” J. Appl. Prob., vol.31, pp.515–525, 1994.
[5]

さて，Perron-Frobenius の定理より，C の固有値

θ に対する固有列ベクトル x1 として x1 > 0 となる
ものがとれる．そして，列ベクトル x2 , . . . , xd が存在して，P = (x1 , x2 , . . . , xd ) とおくと，θ の重複度が1 であることに注意して，C の Jordan 標準形



θ

0

P −1 CP = 
 ..
 .
0

0

···

C̃

0








(A·4)

が得られる．したがって，式 (A·4) より，



e

θ


 0
e C=P
 ..
 .

0

···

−θt Cte(C̃−θI)t C̃

0



 −1
P



0

米倉 誠，“MMPP/M/1 からの退去過程の MMPP に” 長岡技科大修士論文，March 2001.
よる近似，

[6]

(A·3)

が成り立つ．ここで，Re(ηi ) は ηi の実部を表す．

vol.35, pp.998–1002, 1998.
C. Olivier and J. Walrand, “On the existence of(A·2)

が成り立つ．C + αI の固有値は C の固有値に α を加process of an MMPP/M/1 queue,” J. Appl. Prob.,
[4]

(A·1)

が成り立つ．

Re(ηi ) < θ
本研究の一部は，ATR 環境適応通信研究所（現，ATR 適応コミュニケーション研究所）からの平[1]のものを θ とすると，任意の正の d 次元行ベクトルy と正の d 次元列ベクトル z に対して

(A·5)

佐竹一郎，線形代数学，裳華房，1957.

となる．行列 C̃ − θI の固有値 ηi − θ (i = 2, . . . , d)付

録

対角以外の要素が非負である行列 C が既約であるとは，ある正数 α > 0 が存在して C + αI が既約となることである．

[定理] 行列 C は対角以外の要素が非負である既約
な d × d 行列とする．C の実数の固有値のうちで最大は式 (A·3) より Re(ηi − θ) < 0 を満たすので，下記の補題より，t → ∞ のとき

e(C̃−θI)t → O(A·6)となる．したがって，t → ∞ のとき式 (A·5) より

e−θt eCt C → θx1 xT1

(A·7)
1029

電子情報通信学会論文誌 2002/7 Vol. J85–B No. 7を得る．ここで，添字 T は転置を表す．同様にして，e−θt eCt → x1 xT1(A·8)

を得るので，式 (A·7)，(A·8) より，結局，任意の y, z
に対して

lim

t→∞

ye Cz
ye e Cz= limt→∞ ye−θt eCt z
yeCt z
=θ

を得る．

健治

（正員）

昭 55 東工大・理・数学卒．昭 60 同大大学院博士課程満期退学．昭 60 NTT 研究所入社．平 4 長岡技科大工学部助教授．待ち行列理論，ネットワーク特性評価，大偏差理論等の研究に従事．理博．IEEE，情報理論とその応用学会，日本数学会各会員．

−θt Ct

Ct

中川

(A·9)
（証明終）

恩田

和幸

（正員）

平 4 長岡技科大・工・電気・電子システ
ム卒．平 6 同大大学院電気・電子システム
工学専攻修士課程了．同年 NTT アドバン

[補題] 行列 A のすべての固有値の実部が負である
とき，

ステクノロジ（株）入社．平 10 ATR 環境
適応通信研究所研究員．平 13 NTT アド
バンステクノロジ（株）に復帰，現在に至る．

lim eAt = O(A·10)t→∞斎藤

が成り立つ．
（証明） A = N + S という分解を考える．ここで，

S は準単純，N は巾零であり，SN = N S である（文献 [6], p.146 例 1 参照）（この証明中で使用する記号は本文中の記号とは無関係である．むしろ，文献 [6] のまた，文献 [6], p.148 例 2 より

eAt = e(S+N )t1
1
1
N + N 2 +. . .+
N ν−1
1!
2!
(ν − 1)!

となる．A の固有値を α1 , . . . , αn とすると，それらの実部は負であることから




eα1 t

O
..

lim eSt = lim 

t→∞

t→∞

.





eαn t

O
= O

となり，これより補題の結論を得る．

（証明終）（平成 13 年 9 月 13 日受付，14 年 1 月 28 日再受付）

米倉誠

平 11 長岡技科大・工・電子機器卒．平13 同大大学院修士課程了．同年松下電送

システム（株）入社．

1030

（正員）

昭 49 東北大・工・通信卒．昭 54 同大大学院博士課程了．同年 NTT 研究所入所．
平 11 ATR 環境適応通信研究所室長，現在に至る．コヒーレント光伝送，光増幅中継伝送の研究に従事．工博．平 8 本会論文賞受賞．IEEE 会員．

記号に従った）．

= eSt I +茂