論 文
DRR スケジューラのパケット廃棄確率算出における NS-2 へのイン

ポータンスサンプリングシミュレーションの適用について
和泉 光紀†

中川 健治†† a)

横谷 哲也†††

Application of Importance Sampling Simulation to NS-2 for Packet Loss
Probability of DRR Scheduler
Kouki IZUMI† , Kenji NAKAGAWA††a) , and Tetsuya YOKOTANI†††

あらまし ネットワークシミュレータ NS-2 を用いて DRR (Deficit Round Robin) のパケット廃棄確率を推定する．IS (Importance Sampling) シミュレーション法を適用し，シミュレーションを高速化し，かつ推定精度を向上させる．その際，推定値の分散を最小にする分布，いわゆる最適シミュレーション分布を決定して使用する．MC (Monte Carlo) 法と IS 法によるシミュレーション結果を比較し，考察する．
キーワード 不足ラウンドロビン，パケット廃棄確率，NS-2，インポータンスサンプリングシミュレーション

1. ま え が き

通常の FIFO (First-In-Firct-Out) によるパケット処

現在のインターネットでは，パケットの転送品質が特大きくなり著しく不公平な結果となる．このような理では UDP のスループットが TCP のそれよりもに保証されない，いわゆるベストエフォート型のサー

不公平性をなくすため，RR (Round Robin)，WRR

ビスが提供されている．これに対して，次世代ネット(Weighted Round Robin) [5]，DRR (Deﬁclit Round

ワークでは，パケットの遅延や遅延変動，パケット廃

Robin) [9] 等のスケジューラが提案されている．これ

棄確率等の QoS (Quality of Service) を保証するサーらのスケジューラでは，各フローを個別のキューに収ビスが検討されている [4]．各種のアプリケーションを容し，キューから順番にパケットを読み出す．RR で

複数の QoS クラスに分類し，各クラスごとに QoS をは各キューは同じ割合で読み出される．WRR ではフ保証する．ITU-T Rec. Y.1541 [4] では，最もパケッ

ローの重みに従って異なる割合で読み出される．また，

ト廃棄確率の厳しいクラスでは 10−5 の値が想定されDRR では各フローに対してクォンタムというパケッている．したがって，実際のネットワークにおいてこ

ト処理のための Byte 数が割り当てられ，クォンタムのような小さな値を評価する必要がある．

サイズに比例したスループットが各フローに与えられまた，一つのルータにおいて複数のフローが処理る．DRR のアルゴリズムについての詳細は後述する．

されるとき，それぞれのフローに提供されるサービ

DRR はフローごとのスループットの公平性を達成

ス量の公平性も重要である．例えば，ルータで TCPするアルゴリズムとして優れた特性をもち，したがっフローと UDP フローが同時に処理されている場合，て，そのパケット廃棄確率を評価することは実用上，重要なことである．RR や WRR 等の他のスケジュー†（株）コナミディジタルエンタテインメント，東京都
Konami Digital Entertainment Co., Ltd., Tokyo, 176–0022

††

Japan

要な課題である．しかし，例えば，RR や WRR では

長岡技術科学大学，長岡市

DRR におけるクォンタムを適当な値に設定すること

Nagaoka University of Technology, Nagaoka-shi, 940–2188

によって模擬でき，本研究を適用して評価することが

Japan
†††

ラに対してそのパケット廃棄確率を評価することも重

三菱電機株式会社情報技術総合研究所，鎌倉市
Mitsubishi Electric Corporation, Kamakura-shi, 247–8501
Japana) E-mail: nakagawa@nagaokaut.ac.jp698

電子情報通信学会論文誌

できると考える．

DRR では，各フローを個別キューに収容したり，
round robin で個別キューを順番に回ってクォンタム

c （社）電子情報通信学会 2011
B Vol. J94–B No. 5 pp. 698–707 

論文／DRR スケジューラのパケット廃棄確率における NS-2 への IS シミュレーションの適用についてに基づく Deﬁcit Counter 値とパケット長を比較して

2. DRR についてパケットをサーバに送る処理をしたり，と処理プロト
コルが複雑なため，パケット廃棄確率を理論的に解析

DRR (Deﬁcit Round Robin) スケジューラは，各することは困難である．そこで，実際のネットワーク

フローに対して事前に与えられたクォンタムという値

設計は計算機シミュレーションによって行われている．

に比例したパケット送信量を保証するスケジューリン

しかし，C 言語による待ち行列シミュレーションでは，

グアルゴリズムである [9]．DRR のアルゴリズムにつ

パケットに対する各層の処理プロトコルを実装する

いて述べる．
入力フロー数を K とし，各フローのフロー番号をことは難しい．そこで，ネットワークの設計のために
多くの研究者，開発者によって利用されているネットk, k = 1, · · · , K とする．フロー k のパケットはバッ

ワークシミュレータ NS-2 による特性評価を考える．

ファ内の仮想キュー k に収容され，仮想キューの先頭

NS-2 を用いて DRR のパケット廃棄確率シミュレー

パケットが round robin 的に k = 1, 2, · · · の順にスションを行う場合，例えば，DRR へのパケット到着

ケジューラによってアクセスされる．各フロー k にはをポアソン過程などの確率過程とすると，パケット廃

クォンタム δk が割り当てられ，不足カウンタ DCk

棄確率の推定は乱数を使用するモンテカルロ (Monte(Deﬁcit Counter) の値によって先頭パケットをサーCarlo, MC) シミュレーション法によるのが一般的でバに送るかどうかを決める．仮想キュー k が空でなある．NS-2 に MC 法を適用すると，シミュレーション

いときにアクセスのたびに DCk にクォンタム δk がモデルの確率的構造を直接利用することができ，簡単加算され，DCk が先頭パケットのパケット長 lk よにプログラムを作成することができる．しかし，10−5

りも大きいときパケットはサーバに送られて処理され以下の小さい廃棄確率の場合，注目するパケット廃棄

る．DRR の動作アルゴリズムを下記に示す．ただし，事象がなかなか発生せず，シミュレーションに時間がlk ≥ δk とする．

掛かり，また推定値の信頼性も低い．

［DRR の動作アルゴリズム］

そこで，NS-2 によるシミュレーションの高速化と1

Initialize DCk = 0, k = 1, · · · , K

推定値の推定精度の向上のためにインポータンスサン

2

k=1

プリング (Importance Sampling, IS) シミュレーショ

3

if queue k is empty

ン法を適用する．IS 法は実際の確率変数と異なる確率

then DCk = 0

変数を使用して小さい確率の事象をより多く発生させ:goto 6て，得られた値を補正して目標とする確率を得る方法

4

else DCk = DCk + δk

である．

5

if DCk ≥ lk

パケット廃棄確率に関する従来の IS 法では，抽象

then send packet to server

的な待ち行列システムが用いられ，C 言語等でプログ:process packet at server:DCk = DCk − lk

ラムが作成された．しかし，実際のネットワークでの特性評価のためにはパケットのプロトコル処理を含む

6

k = k + 1 (mod K)

様々なアルゴリズムの環境のもとで評価を行う必要が7

go to 3

ある．そこで，本研究の目的を以下のとおりとする．

(I) ネットワークシミュレータ NS-2 を用いて DRRのパケット廃棄確率を推定する．

(II) IS 法を適用し，シミュレーションを高速化し各フローの送信量の公平性に関して次の評価が得られる [9]．

推定精度を向上させる．その際，推定値の分散を最小［定理］ DRR によって t 時間に送信されるフロー kにする分布，いわゆる最適シミュレーション分布を決の Byte 量を sk (t) [Byte] とすると次式が成り立つ．

定して使用する．

(III) MC 法と IS 法によるシミュレーション結果を比較し，考察する．

 sk (t) sj (t)  < 1 + 2lmax
max 
−
k,j
δk
δj 
δmin

ただし，上式において δmin = min{δ1 , · · · , δK } 及びlmax = max{l1 , · · · , lK } とおいた．
699

電子情報通信学会論文誌 2011/5 Vol. J94–B No. 5

3. インポータンスサンプリングシミュレーション法

による推定値の分散よりも小さくすることができる．

minp V [P̂IS ] を達成する p = p∗ が最も精度の良い推定値を与えることになる．標本平均の発生回数を

確率的シミュレーションにおいて，通常のモンテ

M とすると，分散 V [P̂IS ] は M −1 に比例するから，

カルロ (Monte Carlo, MC) 法では得られない小さなIS 推定値の分散を MC 推定値の分散と同程度にすれ確率に対する推定値を得る方法としてインポータン

ば M は小さくて済む．すなわち，シミュレーション

スサンプリング (Importance Sampling, IS) 法があが高速化される．IS 法では minp V [P̂IS ] を達成する

る [1], [8], [10]．ここで，IS 法の概要を説明するために

p = p∗ を得ることが重要である．この分布 p∗ を最正規分布の標本平均の裾確率に対するシミュレーショ

適シミュレーション分布と呼ぶ．正規分布の裾確率にンを考える．

対する IS 法の場合，最適シミュレーション分布 p∗ は

3. 1 正規分布の裾確率に対する IS 法

大偏差理論を利用して次のように得られる [10]．

標準正規分布 N (0, 1) に従う互いに独立な確率変数を X(i), i = 1, · · · , n とし，それらの標本平均を

X̄ ≡

n

i=1

p0 (x) のモーメント母関数を ϕX (θ) とし，p0 (x) に対する指数的測度変換 pθ (x) をX(i)/n とおく．a > 0 を定数として，X̄

の裾確率 P (X̄ > a) をシミュレーションによって推pθ (x) =定したい．a が大きいと通常の MC 法では {X̄ > a}eθx p0 (x), θ ∈ R, −∞ < x < ∞ϕX (θ)(2)となる事象はなかなか発生しないのでシミュレーショ

と定義する [10]．ただし，R は実数全体を表す．J(p ) =ンに膨大な時間が掛かる．また，乱数の精度や周期のlimn→∞ (1/n) log V [P̂IS ] とおくと，minp J(p ) を達問題もある [6]．そこで，IS 法を適用する．IS 法では，成する p = p∗ はX(i) の代わりに他の正規分布に従う確率変数 X (i)を使って，事象 {X̄  > a} が発生しやすくなるようにp∗ = pθ∗ , θ∗ は

する．ここで，X̄ は X (i), i = 1, · · · , n の標本平均を表す．元の X(i) 及びシミュレーションに使う X  (i)の確率密度関数を，それぞれ，p0 (x) 及び p (x) とおϕX (θ)= a の解ϕX (θ)(3)となる [10]．ϕX (θ) は θ による微分を表す．実際に方程式 ϕX (θ)/ϕX (θ) = a を解くと θ∗ = a となり，最適

く．p (x) をシミュレーション分布と呼ぶ．P (X̄ > a)

シミュレーション分布 p∗ = pa は平均 a，分散 1 の正規に対する IS 推定値 P̂IS は次のように与えられる．

分布 p∗ = N (a, 1) であることが分かる．p∗ = N (a, 1)

p0 (x (1)) · · · p0 (xm (n))
1 
1a (x̄m )  mM
p (xm (1)) · · · p (xm (n))

をシミュレーション分布とするとき重み関数 (1) は，

M

P̂IS =

m=1

ここで，xm (i), i = 1, · · · , n, m = 1, · · · , M は m
番目の標本平均を発生させるときの Xi の実現値を表し，その標本平均を x̄m としている．p0 ( ), p ( ) は，それぞれ，確率分布 p0 , p による評価を表す．また，
1a (x̄m ) は条件 x̄m > a が成り立つとき 1，そうでな
いとき 0 となる関数を表す．X̄ によって多く発生している事象 {X̄  > a} の確率を補正するために「重み
関数」

Wm ≡

p0 (xm (1)) · · · p0 (xm (n))
p (xm (1)) · · · p (xm (n))



簡単な計算から Wm = exp −ax̄n (m) + a2 /2 n となる．

3. 2 待ち行列長の裾確率に対する IS 法
IS 法を待ち行列に適用する．定常状態におけるキュー
長を Q とし，裾確率 P (Q > q) に対する IS 法を考える．q が大きいと事象 {Q > q} はなかなか発生しないので，到着率がより大きな値となるシミュレーショ
ン分布 P  を用いてパケットを発生させ，重み関数によって補正して P (Q > q) の推定値を得る．
待ち行列における IS 法では確率 P (Q > q) の推定をサイクル単位で考える [3]．サイクルとはキュー長が(1)

正である一つの時間区間のことである．

3. 3 キュー長の裾確率に対する IS 推定値

が掛け算されている．この補正によって IS 推定値 P̂IS

m 番 目 の サ イ ク ル の 中 で パ ケット の 処 理 終 了は P (X̄ > a) に対する不偏推定量になっている [1]．

時点を tm (1), tm (2), · · · , tm (i), · · · , tm (Im ) とする．

IS 推定値 P̂IS の分散を V [P̂IS ] と表す．適当なシ

P (Q > q) に対する IS 推定値 P̂IS は以下のように与

ミュレーション分布 p を用いると V [P̂IS ] を MC 法
700

えられる [3]．

論文／DRR スケジューラのパケット廃棄確率における NS-2 への IS シミュレーションの適用について

1 1 
1q (Q (tm (i)))Wm (i) (4)·D M
MP̂IS =

Im

4. DRR の パ ケット 廃 棄 確 率 に 対 す る
IS 法

m=1 i=1M

1 
Wm (i)M

D=

Im

4. 1 DRR のモデル化(5)m=1 i=1


P Q (tm (1)), · · · , Q (tm (i))

Wm (i) =P  Q (tm (1)), · · · , Q (tm (i))

DRR の Byte 単位のパケット処理をモデル化する．
フロー数を K とする．DRR に入力されるフロー(6)k, k = 1, · · · , K のパケットは，バッファ内で個別の仮想キュー k に収容されクォンタムサイズ δk [Byte]

ここで，Wm (i) は重み関数である．D は平均サイクが割り当てられる．フロー k のパケットはすべて同一

ル長の IS 推定値であり，(4) の P̂IS の 1/D を除く部の長さ lk [Byte] と仮定する．また，lk ≥ δk とする．

分は 1 サイクル中でキュー長が q を超えている時間のフロー k の到着レートを Λk [Byte/s] とし，パケット

長さの IS 推定値である [3]．

の到着はポアソン過程に従うとする．サーバの処理率

IS 推定値の分散 V [P̂IS ] を最小にする分布 P  = P ∗∗の到着率 λ は一般に処理率よりも大きな値となるのを μ [Byte/s] とし，バッファサイズは無限大とする．
したがって，DRR はパラメータで，P  = P ∗ によるキュー長 Q (t) = Q∗ (t) は時間

S = ((Λk , δk , lk )k=1,...,K , μ)t とともに発散してしまう．そこで，サイクルを強制的に有限時間で打ち切る．打ち切る方法は以下のよう
にする．

i

1 (Q∗ (tm (i)))Wm (i) は十分大きな i
i=1 q

(9)

によって規定される．
本論文では，有限のバッファサイズ q をもつシステに対して収束するので，キュー長 Q∗ (t) が q を超えムにおけるパケット廃棄確率を，バッファサイズ無限る事象を何回か観測したら収束とみなしてそこでサイ

大のシステムにおけるキュー長の裾確率 P (Q > q) でクルを終了させる．経験的には 10〜20 回程度の観測

近似している．そこで，本論文では便宜上，P (Q > q)

で打ち切ると良い結果が得られている．文献 [3] でもをパケット廃棄確率と呼ぶことにする．また，本論文同様のことが示されている．この方法をダイナミック

では検討の最初の段階として，パケット到着を最も基IS という [3]．

本的な到着過程であるポアソン過程としている．今後，

3. 4 最適シミュレーション分布

バースト性を表現可能な更に複雑なモデルについても離散時刻 t(i) におけるもとの分布 P によるキュー

検討してゆきたい．

長を Q(t(i)), i = 1, · · · , I とする．Q(t(i)) がマルコ

4. 2 DRR のキュー長変化の確率的構造

フ連鎖をなすとし，Q(t(i − 1)) > 0 のとき i.i.d.（独立同分布）確率変数 X によってすべての仮想キューにパケットが 1 個以上あるとす
ると，仮想キュー k の先頭パケットが，先頭に位置

Q(t(i)) = Q(t(i − 1)) + X, i = 1, · · · , I(7)

してからサーバに送られるまでに掛かる時間はおよそ

lk /δk に比例する．したがって，長い時間で考えると，と表されているとする．つまり，キュー長の確率は X

サーバで処理されるフロー k のパケット数の割合 πkによって決まる．このとき，P̂IS の分散 V [P̂IS ] を最は小にする X ∗ は次のように与えられる [1]．

X のモーメント母関数を ϕX (θ) =xeθx P (X =x) とし，P (X = x) に対する指数的測度変換を
eθx P (X = x)
P (X = x) =, θ∈RϕX (θ)
δk /lkπk = K, k = 1, · · · , Kδ /l
j=1 j j

(10)

となる．

(8)

と定義する．次の定理が成り立つ [1]．
［定理 A］ V [P̂IS ] を最小にする確率変数 X  = X ∗∗時刻 t における仮想キュー k 内の全パケットの合計

Byte 数を Qk (t) [Byte] とし，Q(t) =

K

k=1

Qk (t)[Byte] とおく．ある k に対して Qk (t) > 0 となる時
間区間，すなわち，サイクルを考える．サイクルにおは指数的測度変換 P (X ∗ = x) = eθ x P (X = x) にいてサーバは連続的に稼働している．一つのサイクル

よって与えられる．ここで，θ∗ は方程式 ϕX (θ∗ ) =C の中で，パケットのサーバでの処理終了時刻を，順1, θ∗ > 0 の解である．

に，t(1), · · · , t(i), · · · , t(I) ∈ C とし，t(1) に処理が701

電子情報通信学会論文誌 2011/5 Vol. J94–B No. 5

終了するパケットの処理開始時刻を t(0) とする．時以上の記号の準備のもとで，X の確率 P (X = x)

間区間 t(i − 1) < t ≤ t(i) に到着するフロー k の及び X のモーメント母関数 ϕX (θ) を計算する．具体

パケット数を nk (i) 個とすると，全フローの Byte 数的な計算過程は付録 1. で述べる．計算結果は以下のA(i) は A(i) =ようになる．

Kn (i)lk [Byte] となる．時間区
k=1 k

間 t(i − 1) < t ≤ t(i) で処理されるパケットのパケッ
ト長を L(i) [Byte] とすると，次の漸化式が成り立つ．

P (X = x) =に つ い て 考 え る ．変 数 i を 省 略 し て A = A(i),ϕX (θ) =KK

R(nj ; λkj ), (16)n·l=x+lk j=1

πk exp f (θ)lk



(17)

k=1

ここで，

より

1  Λj lj θ(e − 1)μlj
K

P (L = lk ) = πk , k = 1, · · · , K,
P (X = x) =

K


(11)

πk P (X = x|L = lk ),(12)

f (θ) = −θ +

(18)

j=1

である．f (θ) は以下の性質をもつ．(i) f (0) = 0,

k=1(ii) f  (0) < 0, (iii) f  (θ) > 0, (iv) limθ→∞ f (θ) =x ≥ −lmax∞ 実際，(i)，(iii)，(iv) はやさしい．(ii) はトラヒッとなる．X の確率 (12) 及び X のモーメント母関数

ϕX (θ) を具体的に計算したい．そのために，いくつかの記号を準備する．

R(n; λ) =ク密度 ρ ≡

j

Λj /μ < 1 より分かる．性質 (i)〜(iv)より f (θ∗ ) = 0, θ∗ > 0 を満たす θ = θ∗ がただ一つ存在する．したがって，方程式 ϕX (θ) = 1 は唯一の解

ϕX (θ∗ ) = 1, θ∗ > 0

まず，λ > 0 と n = 0, 1, · · · に対して，

(19)

n

λ −λ
e
n!

(13)

と お く．こ れ は ポ ア ソ ン 確 率 で あ る ．時 間 区 間t(i − 1) < t ≤ t(i) に処理されるパケットがフロー kのパケット，すなわち，長さ lk [Byte] のパケットであるとする．その処理時間を τk [s] とすると，サーバの処理速度が μ [Byte/s] なので，

lk
τk =
[s]
μ

をもつ．

4. 3 DRR の最適シミュレーション分布
DRR のパケット廃棄確率に対する IS 法の最適なシミュレーション分布及びその分布をもつ確率変数 X ∗を決定する．それは 3. 4 の定理 A より，(19) の θ∗ を用いて指数的測度変換によって与えられる．すなわち，∗P (X ∗ = x) = eθ x P (X = x), x ≥ −lmax (20)(14)

となる．長さ τk の時間区間におけるフロー j のパ
ケットの到着率を λkj [packet/τk ] とおく．フロー jである．(19) より，(20) の X ∗ が確率変数であるこP (X ∗ = x) = ϕX (θ∗ ) = 1x≥−lmax∗lj θ ∗が確認される．λkj = e
λkj とおくと，付録 2. よりと，つまり，の 1 [s] 当りの平均到着レートが Λj [Byte/s] であり，

1 パケット当りのパケット長が lj [Byte/packet] な

P (X ∗ = x) =ので，

λkj = τkx ≥ −lmax ,L = L(i) として，X = A − L とおく．また，最大パケット長を lmax = max(l1 , · · · , lK ) とする．(10)

πk

k=1

Q(t(i)) = Q(t(i − 1))+A(i)−L(i), i = 1, · · · , I
ここで，キュー長差 Q(t(i))−Q(t(i−1)) = A(i)−L(i)

K


K

k=1

Λj
lk Λj
=
[packet/τk ]
lj
μ lj

(15)

πk



K


R(nj ; λ∗kj ) (21)n·l=x+lk j=1となる．

(16) と (21) を比較すると，到着率 λkj が λ∗kj にを得る．したがって，長さ τk の時間区間においてフ置き換えられ，πk , lk は変わらないことが分かる．しロー j のパケットが nj 個到着する確率は R(nj ; λkj )たがって，(15) より，Λ∗k = elk θ Λk , k = 1, · · · , K∗と表される．記号の簡単化のために，n1 , · · · , nK ,とおくと，X ∗ における到着率は X における到着

l1 , · · · , lK に対して，n · l =率 Λk を Λ∗k に置き換えたものである．もとの DRR

702

K

k=1

nk lk と書く．

論文／DRR スケジューラのパケット廃棄確率における NS-2 への IS シミュレーションの適用についてのキュー長変化規則 X を表すシステムパラメータS = ((Λk , lk , δk )k=1,··· ,K , μ) に対して，IS 法の最適5. シミュレーションによる比較評価

なシミュレーションの X ∗ を表すシステムパラメータ

様々な条件のもとで MC 法と IS 法によるパケット

を S ∗ = ((Λ∗k , lk∗ , δk∗ )k=1,··· ,K , μ∗ ) とおくと，上で述

廃棄確率のシミュレーションを行い，推定値の精度，

べたことから S と S ∗ との関係は得られる推定値の限界，シミュレーション時間を比較∗Λ∗k = elk θ Λk , lk∗ = lk , δk∗ = δk , k = 1, · · · , K,μ∗ = μまた，トラヒックのプロトコルとしては UDP を想定している．

4. 4 IS 推定値と重み関数システム S∗による時刻 t での仮想キュー k 内

Q∗k (t)の全パケットの合計 Byte 数を

Q (t) =

率，処理率，パケットサイズ，クォンタムの値を設定
し，NS-2 によってパケット廃棄確率を推定している．

となることが分かる．

∗

する．本研究ではパラメータとして，フロー数，到着

K[Byte] とし，Q∗k (t) [Byte] とおく．3. 3 と同様に，
k=1

m 番目のサイクル Cm においてパケットの処理終了
時刻を tm (1), · · · , tm (i), · · · , tm (Im ) とする．tm (1)5. 1 シミュレーション 1
表 1，表 2 にシミュレーションのパラメータを示す．
表 2 に四つの異なる case における，もとの到着率 Λkと IS で用いる到着率 Λ∗k を示す．シミュレーションの終了条件は，MC 法ではサイクル数が 107 , IS 法ではに処理が終了するパケットの処理開始時刻を tm (0) と104 とした．これらのサイクル数で MC 推定値の分散する．DRR のパケット廃棄確率に対して，システムと IS 推定値の分散がほぼ同程度となる．

∗

S をシミュレーション分布とする IS 推定値 P̂IS は
3. 3 と同様に次のように与えられる．
1 1 
P̂IS = ·1q (Q∗ (tm (i)))Wm (i) (22)
D M
Im

M

m=1 i=1

M

D=

1 
Wm (i)
M
Im

h=1
i



よってはパケット廃棄事象が発生しないためである．
サイクル数が 107 なので，これをもっと大きくすれば，すなわちシミュレーション時間をもっと長くすれはない．MC 法による推定の限界は，生成するポアソ
ン乱数の実数としての桁の問題，乱数の周期の問題等

∗

∗



P Q (tm (h))|Q (tm (h − 1))

Wm (i) =る．これは，10−5 程度以下の廃棄確率では MC 法にば更に小さい廃棄確率まで推定できる，というわけで(23)

m=1 i=1i

図 1 より MC 法では 10−5 程度が推定の限界であによって決まる [6]．一方，IS 法では，非常に小さい
廃棄確率の値も得られていて推定の限界はない．これP ∗ Q∗ (tm (h))|Q∗ (tm (h − 1))h=1

は，IS 法によって廃棄事象が大きな確率で起こるように到着率を大きくしているからであるが，到着率を単に大きくすればいいというものではない．最適な到着(24)率 Λ∗k を用いることによって分散が小さい安定した推Q∗ (tm (0)) = 0 とすると，(24) の重み関数 Wm (i) は計算によって次式 (25) のように簡単な形になる．計
算過程は付録 3. に示す．





Wm (i) = exp − θ∗ Q∗ (tm (i))(25)

その結果，(22) の IS 推定値 P̂IS は簡略化されて∗ ∗1 1 
1q (Q∗ (tm (i)))e−θ Q (tm (i))
·
D M

M

P̂IS =

Im

m=1 i=1MD=

1   −θ∗ Q∗ (tm (i))
eMIm

m=1 i=1となる．

表 1 シミュレーション 1 の各パラメータ
Table 1 Parameters in simulation 1.
フロー数
処理率 [kByte/s]
パケットサイズ [Byte]
クォンタム [Byte]

K=2
μ = 125
l1 = 100, l2 = 200
δ1 = 10, δ2 = 10

表 2 元の到着率 Λk と IS の到着率 Λ∗k , k = 1, 2
Table 2 Original arrival rates and IS arrival rates.
case 1case 2case 3case 4

Λ1 [kByte/s]

12.5

25.0

37.5

50.0

Λ∗
1 [kByte/s]

62.8

69.1

69.0

66.3

Λ2 [kByte/s]

12.5

25.0

37.5

50.0

Λ∗
2 [kByte/s]

315.4

191.0

126.9

88.0

703

電子情報通信学会論文誌 2011/5 Vol. J94–B No. 5

図 1 MC 推定値と IS 推定値（シミュレーション 1）
Fig. 1 MC estimates and IS estimates. (simulation 1)

図 2 MC 推定値と IS 推定値（シミュレーション 2）
Fig. 2 MC estimates and IS estimates. (simulation 2)

表 3 シミュレーション 2 の各パラメータ
Table 3 Parameters in simulation 2.

表 5 シミュレーション 3 の各パラメータ
Table 5 Parameters in simulation 3.

フロー数
K = 100
処理率 [kByte/s]
μ = 125
パケットサイズ [Byte] l1 = · · · = l50 = 100
l51 = · · · = l100 = 200
クォンタム [Byte]
δ1 = · · · = δ100 = 10元の到着率 Λk と IS の到着率 Λ∗k

表4
Table 4 Original arrival rates and IS arrival rates.
case 1case 2case 3case 4フロー数
処理率 [kByte/s]
パケットサイズ [Byte]
クォンタム [Byte]
もとの到着率 [kByte/s]
IS の到着率 [kByte/s]
キュー長しきい値 q [Byte]
標本数

K = 100
μ = 125
l1 = · · · = l100 = 200
δ1 = · · · = δ100 = 10
Λ1 = · · · = Λ100 = 0.375∗
Λ∗
1 = · · · = Λ100 = 2.965

800
20

Λ1 , · · · , Λ50 [kByte/s]

0.125

0.250

0.375

0.500

∗
Λ∗
1 , · · · , Λ50 [kByte/s]

1.057

1.256

1.344

1.382

い値まで推定できている．一方，シミュレーション 1

Λ51 , · · · , Λ100 [kByte/s]

0.125

0.250

0.375

0.500

∗
Λ∗
51 , · · · , Λ100 [kByte/s]と同様に，MC 法では約 10−5 が推定限界である．

8.935

6.309

4.820

3.819

5. 3 シミュレーション 3 高速性能評価
IS 法の高速性能を評価するため，MC 法と IS 法のシ定値が得られている．必ずしも最適でない到着率を用ミュレーション時間を比較する．表 5 のパラメータをいたときの真の確率からのずれについて [3], [7] で検討

考える．評価する値は，q = 800 として P (Q > 800)

されている．

である．標本数は 20 である．すなわち，同一の条件

なお，図 1 の IS において点が二つずつ横に並んでで 20 回のシミュレーションを行って，それらの標本プロットされているが，これは x 軸の q が 50 刻み

平均と標本分散を計算する．P (Q > 800)  10−4 なであり，一方，パケット長が 100 または 200 なのでので MC 法でも精度の良い推定ができ，MC 法と IS

キュー長 Q が 100 の整数倍であるためである．この法のシミュレーション時間を直接比較することができことから，NS-2 によってパケット長が考慮された処

る．次の手順 (a)〜(c) によって比較する．

理が行われていることが分かる．

(a) MC 法によって，1 回のシミュレーションでの5. 2 シミュレーション 2

サイクル数を 107 として 20 回シミュレーションを行

次に，表 3，表 4 に示すパラメータでシミュレー
ションを行う．このシミュレーションではフロー数

い，パケット廃棄確率の 20 個の標本値を得る．それらの標本分散を計算する．

K = 100 と大きな値に設定している．このような条(b) IS 法で，サイクル数を何種類か変えて (a) と件でも図 2 に示すように，IS 法によって非常に小さ同様にそれぞれ 20 回のシミュレーションを行い，標704

論文／DRR スケジューラのパケット廃棄確率における NS-2 への IS シミュレーションの適用について

V [P̂MC ] = (1 − )/n  /n となる．一方，文献 [1]より，V [P̂IS ] 2/n となる．したがって，V [P̂MC ] とV [P̂IS ] を同程度にすると，サイクル数の比では MC :IS = 1 : となる．表 6 によると，1 サイクル当りの時
間は，MC 法では 106×10−7 [秒] であり，IS 法（サイクル数 103 の場合）では 2.5×10−3 [秒] である．したがって，上記の見積りによると，シミュレーション時間の比は MC : IS = 1×106×10−7 : ×2.5×10−3  4 : 103となる．例えば，パケット廃棄確率= 10−5 に対しては，シミュレーション時間の比は MC : IS = 400 : 1となり，大幅にシミュレーション時間が短縮される．
図 3 MC 推定値と IS 推定値のばらつきの比較
Fig. 3 Comparison of sample variance between MC
and IS estimates.

6. む す び
本論文では，DRR のパケット廃棄確率を NS-2 を用いた IS シミュレーションによって推定した．推定

表6サイクル数を変化させたときの標本分散，計算時間の比較
Table 6 Comparison of sample variance and computation
time by varying the number of cycles.

MCIS

値の分散を最小にする最適パラメータを決定し，実際にシミュレーションを行って得られた推定値の精度とシミュレーションの高速性について考察した．

サイクル数

標本平均

標本分散

時間 [秒][6] でも述べられているように，IS 法はシミュレー

107

1.3 × 10−4

2.1 × 10−11

106.0

ションの確率的構造と数学的構造のトレードオフに基

102

1.2 × 10−4

2.5 × 10−10

2.4

3

−4

2.7 × 10−11

2.5

4.1 × 10−12

3.0

10

104

1.2 × 10

1.2 × 10−4づく方法といえる．その意味で IS 法は MC 法と数学的解法の中間的な方法である．確率的構造のみに基づいているのが MC 法であり，数学的構造のみに基づいているのが待ち行列理論による数学的解法である．

本分散を計算する．

MC 法では数学的構造を利用せずに乱数を発生させて(c) 上記 (a)，(b) で得られた分散が等しいとき推いるので無駄が多く，シミュレーションに時間が掛か定値の精度が同じと考えて，シミュレーション時間のる [6]．一方，待ち行列解析では，すべてを数学的に解比較を行う．

こうとするために複雑な問題では解析の困難さに遭遇

図 3 と表 6 より，標本分散が MC 法と同程度となする．IS 法では，大まかな数学的構造に関する知識を3

る IS 法はサイクル数が 10 であり，サイクル数では利用して無駄な乱数の発生を少なくし，その一方で数104 倍，シミュレーション時間では約 40 倍の高速化が学的に複雑な部分はシミュレーションによって推定値実現されている．

を得ている．本論文においては，問題の大まかな数学

5. 4 小さい廃棄確率でのシミュレーション時間の比較

的構造はキュー長の裾確率の大偏差理論的な構造である．特に，重み関数が (25) で表されていて裾確率の大

上記では P (Q > q) = 10−4 付近の値で比較してい偏差理論的な構造がとらえられている．その構造ではるが，もっと小さい値のパケット廃棄確率での比較をとらえられない部分，例えば，NS-2 による DRR の考える．その場合，MC 法では値が得られないので，アルゴリズムやルータにおける各層のパケット処理プ

シミュレーション時間の比は ∞ ともいえるが，もうロトコルの部分，等を乱数による確率シミュレーショ

少し考えてみよう．

ンによって推定値を得ている．

= P (Q > q) とおいて，パケット廃棄事象を確
率のランダム事象（ベルヌーイ事象）と考える

文[1]と，観測回数 n 回中の廃棄事象の発生回数 N は二項分布 Bi(n, ) に従う．この場合，廃棄確率に対する MC 推定値 P̂MC は P̂MC = N/n となるので，

献

J.A. Bucklew, Large Deviation Techniques in Decision, Simulation, and Estimation, Wiley, New York,
1990.

[2]

M. Cottrell, J-C. Fort, and G. Malgouyres, “Large705電子情報通信学会論文誌 2011/5 Vol. J94–B No. 5
deviations and rare events in the study of stochastic
algorithms,” IEEE Trans. Autom. Control, vol.AC[3]=K28, no.9, pp.907–920, Sept. 1983.

k=1

M. Devetsikiotis and J.K. Townsend, “StatisticalK

optimization of dynamic importance sampling parameters for eﬃcient simulation of communication=

pp.293–305, June 1993.

=

K


[4]

ITU-T Recommendation Y.1541,
http://www.itu.int/rec/T-REC-Y.1541-200602-I/en

k=1

[5]

M. Katevenis, S. Sidiropoulos, and C. Courcoubetis,K=purpose ATM switch chip,” IEEE J. Sel. Areas Com-中川健治，“モンテカルロシミュレーション基礎 —推定精度評価の問題点とその克服,” 信学通誌，no.6, pp.11–20,
Sept. 2008.

[7]小川耕司，中川健治，“MMPP/D/1 キューイングにおけ” 信学論（B-I），vol.J80る最適 IS シミュレーション分布，
B-I, no.2, pp.64–73, Feb. 1997.

[8]

=

=

queue,” IEEE Trans. Autom. Control, vol.36, no.12,
pp.1383–1394, 1991.
[9]

M. Shreedhar and G. Varghese, “Eﬃcient fair queuing using deﬁcit round-robin,” IEEE/ACM Trans.
Netw., vol.4, no.3, pp.375–385, June 1996.

[10]

K


K


πk

e

n1 ,··· ,nK ∈Nπk e−lk θ=R(nj ; λkj )

j=1
K

n l θ

e j j R(nj ; λkj )n1 ,··· ,nK ∈N j=1

πk e−lk θ

K




j=1
−lk θ

πk e

exp

K








λkj elj θ − 1

j=1

πk exp

exp λkj elj θ − 1
1  Λj  lj θ
e −1
μ
lj
K

lk −θ +

k=1
K
R(nj ; λkj )Kθ(n·l−lk )k=1

J.S. Sadowsky, “Large deviations theory and eﬃcient simulation of excessive backlogs in a GI/GI/m

K




eθx



k=1

mun., vol.9, no.8, pp.1265–1279, Oct. 1991.
[6]x≥−lmax n·l=x+lk j=1

k=1

networks,” IEEE/ACM Trans. Netw., vol.1, no.3,“Weighted round-robin cell multiplexing in a general-

πk

j=1



πk exp f (θ)lk

k=1を得る．ここで，N は自然数全体の集合を表す．また，(15) より λkj = lk Λj /μlj であることを使った．

R. Srinivasan, Importance Sampling, Applications in

2. 指数的測度変換 P (X ∗ = x) の計算

Communications and Detection, Springer, 2002.

(19) の解 θ∗ によって定義される指数的測度変換

付(20) は計算によって (21) となる．

録

1. X の確率とモーメント母関数
引き続く二つのパケット処理終了時刻におけるキュー
長差 X の確率 P (X = x) 及び X のモーメント母関
数 ϕX (θ) は (16)，(17)，(18) で与えられる．
（証明） 処理中のパケットのパケット長を lk とすると，（証明）

P (X ∗ = x)∗= eθ x P (X = x)∗= eθ xA = X + L より P (X = x|L = lk ) = P (A = x + lk )となる．A =

K

k=1



k=1

nk lk = n · l なので，=

πk P (A = x + lk )
=

k=1

=

K

k=1

K

P (X = x) =

K

k=1

πk



K


R(nj ; λkj )n·l=x+lk j=1

K


K


πk



πk




n·l=x+lkπkR(nj ; λkj )n·l=x+lk j=1

k=1

k=1

K


∗

eθ (n·l−lk )n·l=x+lk

K

=

πk

K


R(nj ; λkj )j=1∗

ef (θ )lkK
∗R(nj ; λkj elj θ )

j=1


K

R(nj ; λ∗kj )n·l=x+lk j=1を得る．また，モーメント母関数については，
最後から二つ目の等号を示すために，(18)，(15) より

ϕX (θ)

=eθx P (X = x)x≥−lmax706

f (θ∗ )lk = −lk θ∗ +

K

j=1

∗

λkj elj θ −

K

j=1

λkj

論文／DRR スケジューラのパケット廃棄確率における NS-2 への IS シミュレーションの適用について= −lk θ∗ +

K

∗

K


j=1

j=1

λkj −

和泉

λkj

光紀

平 20 長岡技科大・工・電気電子情報卒．
平 22 同大大学院修士課程了．同年（株）
コナミディジタルエンタテインメント入社．

であることを使った．

3. 最適な重み関数の計算(24) で定義される重み関数 Wm (i) は計算によって(25) となる．
（証明） 時間区間 tm (h − 1) < t ≤ tm (h) において処理されるパケットがフロー k のパケットであるとし，
この時間区間に到着するフロー j のパケット数を njとする．このとき，(24) の各因子は次のようになる．
∗

P Q (tm (h))|Q (tm (h − 1))

j=1

R(nj ; λkj )j=1

R(nj ; λ∗kj )

=

K




exp −nj lj θ∗

（正員：シニア会員）

昭 55 東工大・理・数学卒．昭 60 同大
大学院博士課程満期退学．昭 60 NTT 研
究所入社．情報理論，ATM 容量設計等の従事．理博．IEEE，情報理論とその応用学会，日本 OR 学会，
日本数学会，日本工業教育協会各会員．

P ∗ Q∗ (tm (h))|Q∗ (tm (h − 1))
= K

健治

研究に従事．平 4 長岡技科大工学部助教
授．平 19 同准教授．ネットワーク特性評
価，待ち行列理論，大偏差理論等の研究に

∗

K

中川

横谷

+ λ∗kj − λkj（正員）昭 60 東京理科大・理工・情報科学卒，



昭 62 同大大学院修士課程了．同年三菱電機（株）入社．以来，高速 LAN，ATM 装置，イーサネットスイッチ，PON システ

j=1= exp (−θ∗ (n · l − lk ) + f (θ∗ )lk )ム，ホームゲートウェイ等の開発並びにト
ラヒック制御，シグナリングプロトコル，
性能評価手法に関する研究に従事．現在，同社情報技術総合= exp (−θ∗ (n · l − lk ))= exp (−θ∗ xm (h))

哲也

(A·1)

ここで，xm (h) = n · l − lk とおいた．xm (h) は時間
区間 tm (h − 1) < t ≤ tm (h) における全到着 Byte
数 n · l からこの区間で処理される Byte 数 lk を研究所通信技術部門主管技師長．また，光アクセスシステム
を中心とした標準化活動に参加．ITU-T SG15，FSAN (FullService Access Network) においてエディタを多数経験．現在，
ITU-T SG15 WP1/Q1 副ラポータ．日本 ITU 協会賞国際活動奨励賞，電子情報通信学会活動功労賞等受賞．博士（理学）．

引いたものであるから，すなわち，xm (h) はこの区間の両端でのキュー長差を表す．つまり，xm (h) =Q∗ (tm (h)) − Q∗ (tm (h − 1)) である．Q∗ (tm (0)) = 0と仮定すると，(24)，(A·1) より

Wm (i) =

i

h=1





exp − θ∗ xm (h)= exp − θ∗ Q∗ (tm (i))
となる．

（平成 22 年 9 月 14 日受付，12 月 20 日再受付）

707