論文

パケットネットワークにおける平均系内時間の近似解析
坂本 将清†

中川 健治†† a)

恩田 和幸†††

斎藤

茂††††

Approximation for the Average Soujourn Time in a Packet Network
Masakiyo SAKAMOTO† , Kenji NAKAGAWA††a) , Kazuyuki ONDA††† ,and Shigeru SAITO††††

あらまし 待ち行列からの退去過程は一般に複雑なものとなるため，ネットワークのあるノードからの退去過程をそのまま次ノードへの到着過程とすると解析は困難となる．従来は解析を容易にするため，ネットワークを
Jackson 網と仮定し，すなわち到着をポアソン過程，サービスを指数サービスとして解析が行われていた．しかしポアソン過程では，実際のパケットネットワークにおけるバーストトラヒックを表現するには不十分である．
米倉ら [1] は，到着過程をバースト性が表現できる 2 状態 MMPP（Markov Modulated Poisson Process）へ拡
張し，サービスを指数サービスとするモデルの解析を行った．本研究ではサービス時間分布を相型分布へ拡張し，到着過程を 2 状態 MMPP とするネットワークにおいて，各ノードの平均系内時間を解析する．モデルの拡張に伴い，本研究では米倉らの解析手法を修正，拡張した．本研究の提案法による平均系内時間の近似値は，米倉らによる従来法よりも精度が良く，計算量も少ないことを示す．
キーワード

待ち行列，平均系内時間，MMPP，相型分布

1. ま え が き（Markov Modulated Poisson Process）へ到着過程を本研究ではパケットネットワークにおける end-to-

米倉ら [1] は，複雑な退去過程を 2 状態 MMPP モデ

拡張したモデルの解析が米倉ら [1] によって行われた．

end 遅延時間を評価することを目的として，ネット

ルで近似し解析する方法を提案した．

ワーク上の各ノードを待ち行列でモデル化し平均系

本研究ではサービスモデルの拡張を行う．各ノード

内時間を解析する．従来ではネットワークを Jackson

におけるサービス時間はパケット長に比例するが，例網と仮定し，つまりパケットの到着過程をポアソン過

えば IP トラヒックでは 40 Byte，及び 1500 Byte の程，サービスを指数サービスとして解析が行われてきパケットの割合が他に比べて多いことが知られていた．Jackson 網ではノードからの退去過程もポアソン

る [2]．前者は TCP の最小パケットサイズであり，後過程となるから，各ノードは M/M/1 待ち行列として者は Ethernet における MTU（Maximum Transfer

容易に解析できる．しかし現実のパケットトラヒック

Unit：最大パケット長）のパケットである．このようはバースト的であるためポアソン過程では不十分でにパケットトラヒックが，ある特定のパケット長の割ある．そこでバースト性が表現できる 2 状態 MMPP

合だけ多いという特徴をもつ場合，指数サービスでは不十分である．そこで本研究では，サービス時間分布†

パナソニックコミュニケーションズ株式会社，東京都を任意の確率分布を近似可能な相型分布へ拡張し，到Panasonic Communications Co., Ltd., 2–3–8 Shimomeguro,
Meguro-ku, Tokyo, 153–8687 Japan††

Department of Electrical Engineering, Nagaoka Universityof Technology, Nagaoka-shi, 940–2188 Japan†††近似する．近似において 2 状態 MMPP のパラメータ

NTT アドバンステクノロジ株式会社，武蔵野市を決定するための近似条件が必要である．米倉ら [1] にNTT Advanced Technology Co. Ltd., 2–4–15 Nakamachi,

よって提案された従来の近似条件では，サービス時間

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

着過程を 2 状態 MMPP とするモデルを解析する．
また本研究でも退去過程を 2 状態 MMPP モデルで長岡技術科学大学，長岡市

立命館大学，草津市

分布を相型分布へ拡張した場合には問題がある．そこ

Faculty of Science and Engineering, Ritsumeikan Univer-で本研究では新たな近似条件を提案する．そして従来

sity, 1–1–1 Noji-higashi, Kusatsu-shi, 525–8577 Japana) E-mail: nakagawa@nagaokaut.ac.jp784

電子情報通信学会論文誌

条件と提案条件によって退去過程を 2 状態 MMPP で
B Vol. J87–B No. 6 pp. 784–793 2004 年 6 月

論文／パケットネットワークにおける平均系内時間の近似解析

V [N (t)] = 2νD 1 [exp{Dt} − I − Dt]
×(D − 1 ν)−2 D 1 1 + νD 1 1 t
図 1 直列型ネットワーク
Fig. 1 Tandem network.

(3)となる．また分散と時間の比 V [N (t)]/t の，t → 0,

近似し各ノードにおける平均系内時間を解析したとき，提案条件による近似値が従来よりもシミュレーション
値に近く，また従来よりも高速な解析が可能であることを示す．

2. ネットワークモデル

t → ∞ における極限値は次のようになる．
V [N (t)](4)= νD 1 1tV [N (t)]lim= −2νD 1 D(D − 1 ν)−2 D 1 1
t→∞
t
lim

t→0

+νD 1 1

本研究ではトラヒックの合流がないパケットネット
ワークを考える．このとき end-to-end 遅延時間は，
図 1 に示す直列型ネットワークを解析することで得られる．本研究ではノード 1 へのパケットの到着過程を 2状態マルコフ変調ポアソン過程（Markov Modulated

Poisson Process, MMPP）とし，各ノードでのサー
ビス時間分布を相型分布（phase type distribution,PH），バッファサイズを無限大とする．ノード 1 からの退去過程をマルコフ型到着過程（Markovian Arrival
Process, MAP）によって近似表現する．そして，ノード 2 以降では前段のノードからの退去過程が次ノード

2. 2 マルコフ変調ポアソン過程
マルコフ変調ポアソン過程（Markov Modulated
Poisson Process, MMPP）は，MAP において到着を伴う {I(t)} の状態推移を禁止した到着過程であり，D 1の非対角要素がすべて 0 となる．2 状態 MMPP は，






−λ0 − q0
λ0 0
q0
D0 =
, D1 =(6)0 λ1
q1−λ1 − q1の表現をもつ．2 状態 MMPP の定常分布 ν は，

ν=への到着過程となる．以下において，マルコフ型到着過程（MAP），マルコフ変調ポアソン過程（MMPP），
相型分布（PH）について説明する．

1[q1 , q0 ]q0 + q1

マルコフ型到着過程（Markovian Arrival Process,

MAP）は，連続時間マルコフ連鎖 {I(t)} の各状態
{0, 1, · · · , M − 1} にそれぞれ一つのポアソン過程が対応した到着過程である．{I(t)} の状態推移には，到着を伴わない推移と到着を伴う推移の 2 種類があり，それぞれの推移速度行列を D 0 , D 1 とする．(D 0 , D 1 )は MAP の表現と呼ばれ，以降は MAP(D 0 , D 1 ) とい(7)である．平均到着率，分散は式 (2), (3) より

λavg =

2. 1 マルコフ型到着過程(5)

λ1 q0 + λ0 q1
q0 + q1

V [N (t)] = λavg t + 2At −

(8)

2A 
1 − e−(q0 +q1 )t
q0 + q1(9)となる．ここで，A=(λ0 − λ1 )2 q0 q1(q0 + q1 )3う表記方法を用いる．また D = D 0 + D 1 は {I(t)}とおいた．

の推移速度行列である．νD = 0, ν1 = 1 を満たす

ν を，MAP(D 0 , D 1 ) の定常分布という．
MAP(D 0 , D 1 ) における (0, t] での到着数を N (t)2. 3 連続相型分布
連続相型分布（phase type distribution, PH）は，
状態 1, 2, · · · , H を一時的状態，状態 0 を吸収状態ととする．文献 [3] より平均 E[N (t)], 平均到着率 λavg

する連続時間吸収マルコフ連鎖 {J(t)} の吸収時間分

は，布として定義される．{J(t)} の推移速度行列を（注 1）

E[N (t)] = νD 1 1 t
E[N (t)]= νD 1 1
t→∞t

λavg = limとなり，分散 V [N (t)] は


(1)
(2)

∗

Q =

0
s

0
S


(10)

（注 1）：0 は全要素が 0 の行ベクトル．1 は全要素が 1 の列ベクトル．

785電子情報通信学会論文誌 2004/6 Vol. J87–B No. 6とする．ここで s = −S1 である．時刻 t = 0 にで与えられる．ここで π 0 , π 1 は，おける一時的状態の分布を α = [a1 , a2 , · · · , aH ] とする．(α, S) は相型分布の表現と呼ばれる．例えば平均µ−1 ，ステージ数 r のアーラン分布は
−rµ rµO.
−rµ . .

α = [1, 0, · · · , 0], S = 

..



. rµ 
O−rµπ 0 B 1 + π 1 B 2 = 0, π 0 1 + π 1 (I − R)−1 1 = 1(13)の解である．また行列 R は

R 2 A2 + RA1 + A0 = O(14)によって定まる相型分布である．ここで α の要素数
は r, S は r × r の正方行列であり，r = 1 とすれば指数分布が得られる．またアーラン分布においてステーを満たす非負最小行列であり，率行列と呼ばれる [3]．

ジ数 r を大きくすると一点分布に近づくから，サービ

3. 2 MAP/PH/1 の平均系内時間

ス時間一定の場合も相型分布によって任意の精度で近MAP/PH/1 の平均系内時間 W̄ は，平均系内数 N̄

似が可能である．更に [0, ∞) 上の一般分布は一点分と平均到着率 λavg を用いてリトルの公式から，布の混合によって近似できる．相型分布は混合について閉じているから，[0, ∞) 上の任意の分布を相型分布により近似できる．

3. MAP/PH/1 待ち行列についてW̄ = N̄ /λavg(15)で計算できる．N̄ = π 1 (I − R)2 1 で与えられる [3]．

3. 1 MAP/PH/1 の定常分布
MAP(D 0 , D 1 ) を到着過程とし，サービス時間分布を表現 (α, S) の相型分布とする待ち行列 MAP/PH/1

3. 3 MAP/PH/1 からの退去過程
図 1 においてノード 1 の待ち行列は式 (11) の推移
速度行列 Q をもつ MAP/PH/1 となる．ノード 1 かを考える．時刻 t での系内数を L(t)，MAP の状態

らの退去過程は Q を退去を伴わない推移速度行列 U 0

を I(t), 相型 分布の状 態を J(t) と表 し，3 次元のと伴う推移速度行列 U 1 に分解することで導かれる．

状 態 空 間 を も つ 連 続 時 間 マ ル コ フ 連 鎖 {X(t)} ={L(t), I(t), J(t)} を考える．以下では L(t) の値をレベル，I(t), J(t) の値をフェーズと呼ぶ．状態を辞書式順序に並べたとき，{X(t)} の推移速度行列は次式の構造をもつ．


B 1B 0OB 2 A 1 A 0





A2 A1 A0 
Q=



.. .. .. 
. . .
O(11)

ここで I M , I H はそれぞれ M × M, H × H 単位行列である．πQ = 0, π1 = 1 を満たす定常分布 π は
行列幾何形式，

786


O
B1B0
 A1 A0
A1 A 0 
U0 = 

.. .. 

. .
O

(16)

U1 = Q − U0

(17)

MAP(U 0 , U 1 ) は無限大の状態数をもつため，ノーB 0 = D 1 ⊗ α, B 1 = D 0 , B 2 = I M ⊗ s
A0 = D 1 ⊗ I H , A1 = I M ⊗ S + D 0 ⊗ I H
A 2 = I M ⊗ s απ = π 0 , π 1 , π 1 R, π 1 R 2 , · · ·



(12)

ド 2 以降の解析を行うためには，MAP(U 0 , U 1 ) を有限状態の MAP で近似する必要がある．

3. 4 レベル K で切詰め
MAP/PH/1 においてレベル K, K + 1, · · · , (K >1) をまとめて K̂ とする．MAP/PH/1 において L(t)が K 以上で推移している間の退去間隔はサービスと同じ相型分布に従うから，退去過程は次の MAP(Û 0 , Û 1 )で近似できる [4]．

論文／パケットネットワークにおける平均系内時間の近似解析


B 1B 0O
A 1 A0


.. ..

. . 

A1 E 0 
O
E1


O
O
B 2 O





A2 O

Û 1 = 




.. ..

. . 
O
E2E3



Û 0 = 




(18)

Fig. 2(19)

ノードへの到着過程を 2 状態 MMPP とすることができるため，解析が容易になる．
図 2 に示すように，ノード 1 への到着過程である

E 0 = D 1 1 ⊗ I, E 1 = S, E 2 = s y K−1
E 3 = (1 − y K−1 1 )s α
ここで y K−1 は，レベル K 以上でそれぞれのサービスが終了し，レベル K − 1 へ帰還した直後の各フェーズの分布を表し，次式によって与えられる．

y K−1 =

図 2 2 状態 MMPP による退去過程の近似
Approximation of the departure process by 2-state
MMPP.

xK−1
xK−1 1

(20)

式 (20) の xK−1 は，退去発生という条件付きでのレベル K − 1 のフェーズの分布を表し，次式で与えられる．

2 状態 MMPP を MAP(D 0 , D 1 ) と表し，退去過程
である MAP(U 0 , U 1 ) を式 (18), (19) で与えられる
MAP(Û 0 , Û 1 ) で近似する．更に MAP(Û 0 , Û 1 ) を 2
状態 MMPP である MAP(D ∗0 , D ∗1 ) で近似する．ここで D ∗0 , D ∗1 は−λ∗0 − q0∗
λ∗0 0
q0∗
∗
∗
D0 =
, D1 =0 λ∗1
q1∗−λ∗1 − q1∗
であり，パラメータ λ∗0 , λ∗1 , q0∗ , q1∗ を決定するための近似条件が必要である．

4. 1 必要十分条件
Bean [4] によると，MAP(Û 0 , Û 1 ) と MAP(D ∗0 ,xK−1 = π 1 RK−1 A2 (νD 1 1 )−1 , K > 1

D ∗1 ) が確率過程として一致するための必要十分条件
は，任意の自然数 k と任意の実数 t1 , · · · , tk ∈ [0, ∞)MAP(Û 0 , Û 1 ) の定常分布 ν̂ は，に対して次式が成り立つことである．

ν̂ = [π 0 , π 1 , π 1 R, · · · , π 1 RK−2 ,
π 1R

K−1

(I M H − R)−1

(1 ⊗ I M )](21)

として得られる．ここで R は式 (14) で与えられる率行列であり，1 は要素数 M の列ベクトル，I M H は

M · H × M · H 単位行列である．
図 1 でノード 1 からの退去過程を MAP(Û 0 , Û 1 )で近似したとき，その状態数は m1 = 2+2KH −H とν̂eÛ 0 t1 Û 1 eÛ 0 t2 · · · Û 1 eÛ 0 tk Û 1 1ν̂eÛ 0 t1 Û 1 eÛ 0 t2 · · · Û 1 eÛ 0 tk 1=∗∗∗ν ∗ eD 0 t1 D ∗1 eD 0 t2 · · · D ∗1 eD 0 tk D ∗1 1∗∗∗ν ∗ eD 0 t1 D ∗ eD 0 t2 · · · D ∗ eD 0 tk 1 

1

(22)

1

ここで ν ∗ は MAP(D ∗0 , D ∗1 ) の定常分布を表す．米倉ら [1] は，式 (22) において k と t1 , · · · , tk を特別な値に限定したものを近似条件としている．

なる．正確な近似を得るには，トラヒック強度が高いほ4. 2 従来の近似条件（従来法）

ど K を大きくしなければならない．またノード 2 以降

米倉ら [1] により提案された近似条件を示す．

でも MAP(Û 0 , Û 1 ) によって退去過程を近似したとき，ノード n での状態数は mn = mn−1 + mn−1 KH − Hとなる．このように後のノードほど状態数が増加する
ため，計算量が増加し解析が困難となる．

4. 2 状態 MMPP による退去過程の近似

条件 i ） ν̂ Û 1 1 = ν ∗ D ∗1 1
2

条件 ii ） ν̂ Û 1 1 = ν ∗ D ∗1 2 1
条件 iii ） g(Û 0 ) = g(D ∗0 )条件 iv ） W̄1 = W̄2
ここで g(Û 0 ), g(D ∗0 ) は，それぞれ行列 Û 0 , D ∗0 の

米倉ら [1] は 2 状態 MMPP モデルによって退去過

最大固有値を表す．また W̄1 は MAP(Û 0 , Û 1 ) をノー

程を近似する手法を提案した．この近似によって全

ド 2 への到着過程としたときのノード 2 における平787

電子情報通信学会論文誌 2004/6 Vol. J87–B No. 6

均系内時間を表し，W̄2 は MAP(D ∗0 , D ∗1 ) を到着過程

れの分散と時間の比を次のように定義する．

としたときの同様の値を表す．以下で各条件についてv̂(t) ≡

説明する．

4. 2. 1 条 件 i ）条件 i ）は式 (22) で k = 1, t1 = 0 として得られる．条件 i ）は MAP(Û 0 , Û 1 ) と MAP(D ∗0 , D ∗1 ) の平均到着率が一致することを要請する．

4. 2. 2 条 件 ii ）条件 ii ）は k = 2, t1 = 0, t2 = 0 とし，条件 i ）4. 2. 3 条 件 iii ）条件 iii ）は k = 1, t1 → ∞ として得られる．
limv ∗ (t) ≡V [N ∗ (t)]t(23)

更に v̂(t), v ∗ (t) それぞれの t → ∞ における極限値をv̂(∞), v ∗ (∞) とする．これらは式 (5) で与えられる．
本研究で提案する条件を次に示す．
条件 1 ） ν̂ Û 1 1 = ν ∗ D ∗1 1
条件 2 ） g(Û 0 ) = g(D ∗0 )条件 3 ） v̂(t1 ) = v ∗ (t1 )を適用して得られる．

t1 →∞V [N̂ (t)],tν ∗ eD D ∗1 1ν̂eÛ 0 t1 Û 1 1
= lim
∗
t
→∞

1
ν ∗ eD 0 t1 1ν̂eÛ 0 t1 1∗0 t1

上式の両辺は，それぞれ行列 Û 0 , D ∗0 の最大固有値に収束する [1]．よって条件 iii ）は，行列 Û 0 , D ∗0 の最大固有値が一致することを要請する．

4. 2. 4 条 件 iv ）MAP(Û 0 , Û 1 ) を ノ ー ド 2 へ の 到 着 過 程 と し たと きの ，ノー ド 2 に おけ る 平均 系 内時 間 W̄1 と，
MAP(D ∗0 , D ∗1 ) としたときの平均系内時間 W̄2 が等
しくなることを要請する．

4. 2. 5 パラメータの決定
λ∗0 , λ∗1 , q0∗ , q1∗ のうち一つを x とおき，残りの三つを条件 i ）ii ）iii ）から x の関数として表す．そして，条件 iv ）を満たす x をニュートン法により探索する．

4. 2. 6 従来法の問題点
サービス時間分布を相型分布へ拡張したときの，従来法の問題点について述べる．まず条件 ii ）において右辺の D ∗1 は対角行列であるから D ∗1 2 = O である．
2
対して左辺では Û 1 = O となるケースがあり（注 2），こ条件 4 ） v̂(∞) = v ∗ (∞)

4. 3. 1 条件 3 ）と条件 4 ）について
条件 1 ），2 ）は米倉らの条件 i ），iii ）と同じである．トラヒックのもつバースト性が平均系内時間に影響を与えるので，MAP(Û 0 , Û 1 ) と MAP(D ∗0 , D ∗1 )のバースト性を近似することを考える．バースト性は
分散 V [N̂ (t)] と V [N ∗ (t)] によって規定されると考えられるので，理想的にはすべての時間区間 0 ≤ t < ∞にわたって V [N̂ (t)] と V [N ∗ (t)] を一致させることが望ましい．しかしそれは困難であるため，条件 3 ），4 ）では MAP(Û 0 , Û 1 ) と MAP(D ∗0 , D ∗1 ) の分散と時間の比を t = t1 と t = ∞ の 2 点で一致させる．ここでt1 は，曲線 v̂(t) に v ∗ (t) が最も近似できる位置に設定する．
また本論文では，トラヒック強度に関係なく K = 10として MAP(Û 0 , Û 1 ) を得る．

4. 3. 2 t1 の決定方法
式 (4) より，式 (23) で与えられる v̂(t) は t → 0 のとき MAP(Û 0 , Û 1 ) の平均到着率 λ̂avg に収束する．
本研究では t1 を，v̂(t1 ) =λ̂avg + v̂(∞)2(24)の場合は従来法により退去過程を 2 状態 MMPP モデから求める．式 (24) は経験的に得た決定方法であり，ルで近似することは不可能である．

指数サービスの場合に良好な近似結果が得られる．ま更に条件 iv）では平均系内時間を解析する必要があた本論文では指数サービス以外のサービスモデルにつり，このための計算量は相型分布の状態数が大きいほいても式 (24) から t1 を決定する．また式 (24) を t1

ど増大するため，近似が困難となる．

について解くことは困難であるが，指数サービスの場

4. 3 提案する近似条件（提案法）

合は以下のような方法で t1 の近似値が得られる．

本 研 究 で 提 案 す る 近 似 条 件 に つ い て 述 べ る ．本研究では問題がある従来の条件 ii）条件 iv）を変更
する．以下では (0, t] における MAP(Û 0 , Û 1 ) 及びMAP(D ∗0 , D ∗1 ) の到着数をそれぞれ N̂ (t), N ∗ (t) とする．MAP(Û 0 , Û 1 ) 及び MAP(D ∗0 , D ∗1 ) のそれぞ
788

指数サービスにおける v̂(t) の t → 0 における接線

f1 (t) を求める．この接線の傾き ρ は，2

ρ = ν̂ Û 1 1 − (ν̂ Û 1 1 )2（注 2）：例えばサービス時間分布がアーラン分布のとき．

(25)

論文／パケットネットワークにおける平均系内時間の近似解析∗∗

1 − e−(q0 +q1 )t1v̂(∞) − v̂(t1 )=(q0∗ + q1∗ )t1v̂(∞) − λ̂avg(28)2(λ∗0 − λ∗1 )2 q0∗ q1∗= v̂(∞) − λ̂avg(q0∗ + q1∗ )3(29)

ここで，式 (27) では g(Û 0 ) = θ̂ とした．またq ∗ = q0∗ + q1∗(30)とおくと式 (28) は∗

図 3 指数サービスにおける t1 の決定方法
Fig. 3 The determination method of t1 in case of
exponential service.

v̂(∞) − v̂(t1 )1 − e−q t1=
q ∗ t1v̂(∞) − λ̂avg(28 )となり，ニュートン法によって q ∗ が決定できる．式(26), (27), (29), (30) の連立方程式を λ∗0 , λ∗1 , q0∗ , q1∗ にと計算され，切片は λ̂avg である．そして図 3 のよ
うに，ついて解き，λ∗0 =

f1 (t) = ρt + λ̂avg

λ∗1 =

λ̂avg + v̂(∞)
f2 (t) =2a4(λ̂avg + θ̂)
2λ̂avg (λ∗0 − λ̂avg ) − q ∗ c2(λ∗0 − λ̂avg )2(λ∗0 − λ̂avg )2 q ∗2(λ∗0 − λ̂avg )2 + q ∗ cが交差する点 t1 を t1 の代わりに用いる．アーランq0∗ =サービスの場合は U 21 = O となるから，式 (25) より

q1∗ = q ∗ − q0∗

ρ < 0 である．よってアーランサービスの場合の v̂(t)は，指数サービスの場合とは性質が異なる．アーラン
サービスの場合も式 (24) から t1 を決定するが，指数サービスのように t → 0 における接線を利用することができないため，数値的に t1 を決定する．
アーランサービスに対しても式 (24) を使う理由は，
指数サービスに対して精度の良い近似が得られている
ことと，アーランサービスはステージ数 r = 1 とすれば指数サービスとなるのでアルゴリズムの一貫性のためである．しかし，t1 の決め方については指数サービ
ス，アーランサービスともにもっと良い方法があれば精度が向上する．これについては 6. 1 で議論する．

4. 4 パラメータの決定
条件 1 ）〜4 ）よりそれぞれ次の 4 式が得られる．

λ∗1 q0∗ + λ∗0 q1∗= λ̂avg(26)
q0∗ + q1∗1
−λ∗0 − λ∗1 − q0∗ − q1∗ +2
(λ∗0 +λ∗1 +q0∗ +q1∗ )2 − 4(λ∗0 λ∗1 +λ∗1 q0∗ +λ∗0 q1∗ )
= θ̂

(27)

(31)
(32)
(33)
(34)

が得られる．ここで，



a=b+

b2 + 8(λ̂avg + θ̂){2dλ̂avg + θ̂q ∗ c}b = 2(λ̂avg + θ̂)(λ̂avg − θ̂ − q ∗ ) + q ∗ c
c = v̂(∞) − λ̂avg
d = (λ̂avg + θ̂)(q ∗ + θ̂)
とおいた．

5. 提案法と従来法の比較
図 1 において N = 10 としたときの，各ノードの平均系内時間のシミュレーション値（sim）と，従来法（yonekura），提案法（proposed）を比較する．ノード 1 への到着である 2 状態 MMPP のパラメータを(λ0 , λ1 , q0 , q1 ) = (1.0, 0.0, 0.01, 0.01) とする．サービスとして指数サービスとステージ数 r = 2, 10 のアーランサービスを考える．ノード n におけるサービ
スとして表 1 に示す M1 と M2 の 2 種類について調べる．M1 は基本的な検討のためにすべてのノードのサービス率が等しい場合であり，M2 はネットワークの end 側のサービス率が低くコア側のサービス率が高い場合をモデル化している．結果を図 4〜図 9 に示す．
789

電子情報通信学会論文誌 2004/6 Vol. J87–B No. 6
表 1 各ノードのサービス率
Table 1 Service rate of each node.
M1
M2

µ1
1
1

µ2
1
2

µ3
1
3

µ4
1
4

µ5
1
5

µ6
1
5

µ7
1
4

µ8
1
3

µ9
1
2

µ10
1
1

図 7 各ノードの平均系内時間（指数サービス，M2 ）Fig. 7 Average sojourn time at each node in case of
exponential service and service rate is M2 .

図 4 各ノードの平均系内時間（指数サービス，M1 ）Fig. 4 Average sojourn time at each node in case of
exponential service and service rate is M1 .

図 8 各ノードの平均系内時間（アーランサービス，M2 ）Fig. 8 Average sojourn time at each node in case ofErlang service (r = 2) and service rate is M2 .

図 5 各ノードの平均系内時間（アーランサービス，M1 ）Fig. 5 Average sojourn time at each node in case ofErlang service (r = 2) and service rate is M1 .

図 9 各ノードの平均系内時間（アーランサービス，M2 ）Fig. 9 Average sojourn time at each node in case ofErlang service (r = 10) and service rate isM2 .

図 6 各ノードの平均系内時間（アーランサービス，M1 ）Fig. 6 Average sojourn time at each node in case ofErlang service (r = 10) and service rate is
M1 .

6. 考

察

近似値はシミュレーション値とほぼ一致しており従来法よりも精度が良い．また近似値はシミュレーション
値に対して安全側である．提案法の精度が良い理由
は，分散–時間曲線をフィッティングする近似条件を採用したためである．一般に待ち行列では，平均到着率

6. 1 近似精度についてが一致していても到着過程のバースト性が強いほど平指数サービスの場合（図 4，図 7），提案法による

均キュー長が伸びるため平均系内時間が大きくなる．

790

論文／パケットネットワークにおける平均系内時間の近似解析

よって退去過程の近似においてバースト性を無視すると近似精度が悪い．バースト性は分散–時間曲線によって定義されるから，本研究ではこれをフィッティング
するための近似条件を設けた．従来法における近似条
件はバースト性を考慮していなかったため誤差が大きくなった．なお従来法では次ノードにおける平均系内時間を一致させる条件を設けていたため，若いノードでの近似値は誤差が小さい．しかしノードごとに近似を行うため誤差が累積し，後のノードほど近似値の誤差が拡大する．
アーランサービスの場合（図 5，図 6，図 8，図 9），
ステージ数が大きいほどノードごとの平均系内時間の減少が急となる．これはステージ数が大きいほど一定

図 10 アーランサービスにおける v̂(t), v ∗ (t) のフィッティ
ングの一例 (r = 10)
Fig. 10 An example of the fitting of v̂(t) and v∗ (t)in case of Erlang service (r = 10).

サービスに近づき，また完全な一定サービスの場合にはノード 2 以降での待ち時間は 0 となるためである．
従来法では退去過程を近似できないため，図 5，図 6，
図 8，図 9 における従来法は，アーランサービスと同じサービス率をもつ指数サービスとして解析を行っ

Table 2表 2 K の値による計算値の比較
Comparison of computation values withvarying K.

K
Proposed (M1 , r = 2)
Proposed (M1 , r = 10)

10
3.99
3.41

15
4.02
3.42

20
4.02
3.42

simulation
3.21
1.81

ており，従来法による近似値は安全側だが誤差が大きい．提案法では従来法よりもシミュレーション値に近の値を大きくすると後に示すように計算時間が非常にい近似値が得られる．しかし，提案法でも指数サービ

大きくなる．以上のことからレベル K の値としては

スのときと比べると誤差が大きく，ステージ数が大きK = 10 で十分であることが分かる．

いほど誤差が拡大する．アーランサービスのとき誤差

各ノードのサービス率が M2 の場合（図 7，図 8，が拡大する原因の一つとして，分散–時間曲線のフィッ

図 9）における提案法の近似精度が M1 の場合（図 4，

ティングが正確に行えないことが挙げられる．図 10

図 5，図 6) よりも良くなっているように見える理由はは，r = 10 のアーランサービスにおけるノード 1 か

次のとおりである．M2 ではノード 2 からノード 9 にらの退去過程について，式 (23) の v̂(t) と v ∗ (t) をプ

おけるサービス率は M1 のサービス率よりも大きい．

ロットしたものであり，正確なフィッティングが行え

サービス率が大きいとノードにおいてパケットの待ちていないことが分かる．近似モデルの状態数を上げパはほとんど発生せず，系内時間はサービス時間にほぼ

ラメータ数を増やせばより正確なフィッティングが可等しくなる．したがって，この場合，トラヒックのモ能であるが，パラメータは数値計算によってしか決定

デル化が仮に適切でなくても平均到着率を一致させるできなくなり，計算量が増加し実用的でなくなる．ま程度のモデル化であれば系内時間の計算値はほぼ真値た提案条件 3 ）における t1 の選び方を工夫することに等しくなる．つまり，上で議論した 2 状態 MMPP

で，近似精度は向上する．

によるモデル化の限界が目立たなくなっているだけで次に K の値による近似精度への影響について考える．

ある．

図 5 に対応してサービス率 M1 ，ステージ数 r = 2 の

6. 2 計算量について場合，及び図 6 に対応してサービス率 M1 ，ステージ

提案法による退去過程の近似では，式 (23) の v̂(t1 )数 r = 10 の場合において，K の値を K = 10, 15, 20を求めるため，式 (19) で与えられる行列 Û の指数と変化させてノード 10 での平均系内時間を提案法で関数 exp{Û t} を計算する必要がある．本論文では
K = 10 としたから，指数サービスの場合は Û は計算した値，及びシミュレーション値を表 2 に示す．

K の値を大きくしても得られる値はほとんど変わらな
い．これは，提案法が各ノードにおいて前段からの複雑な退去過程をより簡単な MMPP で近似していることによる限界であると考えられる．それに対して，K

22 × 22 の行列となり計算量はほとんど問題とならな
い．式 (24) の t1 は，指数サービスの場合は t → 0 における v̂(t) の接線を利用することで高速にその近似解が得られる．また提案法では式 (28 ) からニュート791

電子情報通信学会論文誌 2004/6 Vol. J87–B No. 6
表 3 解析に要した CPU Time の比較（指数サービス）
Table 3 Comparison of CPU Time required to analysis
in case of exponential serivce.
M1
M2

Proposed
0.21[s]0.13[s]Yonekura
1.21[s]0.51[s]

ratio
17.4%
25.5%

表 4 解析とシミュレーションの CPU Time の比較（指数サービス）
Table 4 Comparison of CPU Time required to analysis and simulation in case of exponential
service.
Proposed
Simulation

10 nodes0.19[s]15.0[s]

20 nodes0.27[s]27.0[s]

30 nodes0.34[s]51.0[s]表 5 K 及び r の値による計算時間の比較
Table 5 Comparison of computation time with varying
K and r.
r ＼K
2
5
10
15

10
6.96[s]71.43[s]1.25 × 103 [s]4.13 × 103 [s]15
17.97[s]2.10 × 102 [s]4.55 × 103 [s]1.50 × 104 [s]20
39.70[s]1.41 × 103 [s]1.08 × 104 [s]3.65 × 104 [s]間分布を相型分布とした直列型ネットワークにおいて，各ノードの平均系内時間を解析した．通常の解析では到着過程の状態数がノードごとに増加し解析が困難となるため，本研究では退去過程を 2 状態 MMPP モデ
ルで近似した．この近似においてサービス時間分布をン法によって q を求める必要があるが収束は速い．更

相型分布としたとき，米倉ら [1] による従来の近似条に q が求まれば式 (31)〜(34) からパラメータが決定

件では問題があり，例えばアーランサービスの場合にできる．指数サービスの場合において，10 ノードをすは，従来法では退去過程を近似不能である．そこで提べて解析し終えるまでの CPU 時間を，提案法，従来案法では，分散–時間曲線をフィッティングする近似条

（注 3）

．本

件を採用した．このフィッティング条件によって，従研究では多くの解析を行い，提案法での解析に要する来法よりもシミュレーション値に近く，かつ安全側のCPU 時間が従来法よりも平均で 70% 以上短縮できる

近似値を得ることができた．また提案法での解析に要ことを確認した．

する CPU 時間は，指数サービスの場合，従来法より法それぞれ測定した．結果を表 3 にまとめる

アーランサービスの場合，行列 Û のサイズはステー

も平均で 70%以上短縮した．アーランサービスの場

ジ数に比例して増加する．また式 (24) から数値的にt1 を決定する必要がある．よってステージ数が大きい

合，提案法ではレベル K 及びステージ数 r の増加にほど計算量が増加してしまう問題がある．

伴い計算量が増加した．
本研究ではトラヒックの合流がないパケットネット

次に，提案法での解析に要する時間とシミュレーショ

ワークを対象として解析を行った．今後は更に一般的

ンに要する時間を比較する．計算時間は，ともに，最終に，トラヒックの合流を考慮したネットワークモデル

ノードにおける平均系内時間の値を相対誤差 10−5 以について解析を進めていく．

下で得られるまでの時間である．ノード 1 への到着を謝辞

本研究の一部は，ATR 環境適応通信研究所

パラメータ (λ0 , λ1 , q0 , q1 ) = (1.0, 0.0, 0.01, 0.01) を（現 ATR 適応コミュニケーション研究所）からの平成

もつ MMPP とし，すべてのノードのサービスはサー11〜13 年度の助成研究費によって行われました．まビス率 1 の指数サービスとする．ノード数 = 10, 20, 30

た，論文改訂の際の数値計算に御協力頂いた長岡技術に対する計算時間の比較を表 4 に示す．提案法による科学大学修士 1 年生の立原和氏に感謝致します．
文

解析はシミュレーションと比較して十分に高速である
といえる．

[1]

ヒックを入力とするネットワークにおける平均系内時間のまた，指数サービス以外のときの提案法の計算時間

近似解析，” 信学論 (B)，vol.J85-B, no.7, pp.1021–1030,
July 2002.

と K, r との関係を考える．表 5 にサービス率が M1の場合に K と r の値を変化させて得られた提案法[2]

三屋光史朗，長健二朗，加藤 朗，村井 純，“パケットサ
イズ分布からみた IP トラフィックの傾向，” インターネッ[3]

牧本直樹，待ち行列アルゴリズム，朝倉書店，2001.
N. Bean, D. Green, and P. Taylor, “Approximationsの計算時間を示す．この結果によると計算時間はほぼ

K 3 r 3 に比例している．したがって現時点では，指数サービス以外の場合には提案法の計算量は大きいといえる．これを改良することは今後の課題である．

7. む す び献米倉 誠，中川健治，恩田和幸，斎藤 茂，“バーストトラ

トコンファレンス 2000 論文集，pp.47–56, Nov. 2000.
[4]

to the output process of MAP/PH/1 queues,” in Advances in Matrix Analytic Methods for Stochastic（注 3）：計 算 機 の 性 能：Intel Pentium III, 871 MHz, Memory:

本研究では到着過程を 2 状態 MMPP，サービス時
792

320 MByte, OS: Windows2000, Application: MatLab.

論文／パケットネットワークにおける平均系内時間の近似解析
Models, ed. A. Alfa and S. Chakravarthy, pp.151–
170, Notable Publications, NJ, 1998.
[5]

N. Bean, D. Green, and P. Taylor, “The output process of an MMPP/M/1 queue,” J.Appl. Prob., vol.35,
pp.998–1002, 1998.

（平成 15 年 6 月 12 日受付，11 月 2 日再受付）

坂本

将清

2001 長岡技科大・工・電気・電子システ

ム卒．2003 同大大学院修士課程了．同年
パナソニックコミュニケーションズ（株）
入社．

中川

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

恩田

和幸 （正員）平 4 長岡技科大・工・電気・電子システ
ム卒．平 6 同大大学院電気・電子システム
工学専攻修士課程了．同年 NTT アドバン
ステクノロジ（株）入社．平 10 ATR 環
境適応通信研究所研究員．平 13 NTT ア
ドバンステクノロジ（株）に復帰，現在に至る．

斎藤

茂 （正員）

昭 49 東北大・工・通信卒．昭 54 同大大学院博士課程了．同年 NTT 研究所入所．
平 11 ATR 環境適応通信研究所室長．平15 立命館大学理工学部教授．コヒーレント

光伝送，光増幅中継伝送，アドホックネッ
トワークの研究に従事．工博．平 8 本学会
論文賞受賞．IEEE 会員．

793