待ち行列：平均値解析法での重複組合せ

BCMPネットワークで重複組合せの総組み合わせに対して計算が必要です。そのために、重複組み合わせの計算をやっておきます。
参考資料として、
https://www.elsevier.com/books/probability-statistics-and-queueing-theory/allen/978-0-08-057105-8
を参考にしていきます。

(1) 重複組合せの計算(全ての組合せの算出)
まず、計算の難しさの部分で出てくる重複組合せの計算をしていきます。

重複組合せを再帰計算でやってみます。
重複組合せは
c種類のものから重複を許してn個選ぶ時の組合せ
cHn = c+n-1Cn
2H5 = 6C5 = 6通り(5,0),(4,1),(3,2),(2,3),(1,4),(0,5)
となります。この時、2種類でやる時は２重ループ、3種類でやる場合は３重ループになります。
public int[][] getCombination() {
int value[][] = new int[6][c];
int index = 0;
for(int i = 0; i <= n; i++) {//第１層
value[index][0] = i;
for(int j = 0; j <= n-i;j++) {//第２層
value[index][1] = j;
//最下層
value[index][2] = n - i - j;
index++;
}
}
return value;
}
これだと、種類が増えるごとにループが増えてしまうので、通常こういう場合は再帰計算でやります。
public void GetRecursive(int n, int c_index, int pre_i) {
if(c_index == c - 1) {
value[index][c_index] = n;
index ++;
return;
}
else {
for(int i = 0; i <= n; i++) {
if(index == this.combi)break;
if(c_index > 0) value[index][c_index - 1] = pre_i;
value[index][c_index] = i;
this.GetRecursive(n - i, c_index + 1, i);
if(c_index == c-1 ) c_index = 0;
}
}
}

例えば、n=5,c=3でやると

重複組合せ:個数21

重複組合せ:[[0, 0, 5], [0, 1, 4], [0, 2, 3], [0, 3, 2], [0, 4, 1], [0, 5, 0], [1, 0, 4], [1, 1, 3], [1, 2, 2], [1, 3, 1], [1, 4, 0], [2, 0, 3], [2, 1, 2], [2, 2, 1], [2, 3, 0], [3, 0, 2], [3, 1, 1], [3, 2, 0], [4, 0, 1], [4, 1, 0], [5, 0, 0]]

となります。

[追記(2019/10/12)]
今回の場合は、クラスごとの最大数を制約として設定する必要があったので、追加しました。一旦出た解に対し、制約を満たしているものだけ取り出しました。ソースコードに関数として追加しました。

ソースコード

(確認すること)
・クラス間の移動が無いモデルでは、トラフィック方程式が必要無い？(初期のクラス人数だけで解けてしまう？）

待ち行列：BCMP Network (3) 平均値解析法algorithmについて

20191021追記
以下の方法で一つ懸念がどこにトラフィック方程式の解の情報が入るのかがわかりませんでした。これがないと客の移動が表現できません。よくよく読んでいくと、Dcmがaverage total demand of a class c job for service at queue mとなっていてサービスデマンドとなっています。このDckは

で表現されています。ここでVmは、average number of visits a single customer makes to queue mとなっているので、トラフィック方程式のαと考えて良さそうです。これをクラス別にしたものを利用することになります。つまり、よく目にする、ρ = α / μ の式になります。Dckは利用率を表していたことになります。

20191025追記
下記のアルゴリズムの中で for each feasible population n = となっています。これはnを増やした時、それぞれ重複組み合わせが必要となります。その重複組み合わせの計算方法はhttps://mizunolab.sist.ac.jp/2019/10/5.html です。

平均値解析法で再起計算が必要でした。その内容を確認します。
http://www.cs.ucr.edu/~mart/204/MVA.pdf
で紹介されているアルゴリズムを参照します。

ただし、Ac,mは下記で求めます。

この計算例として、

があるので、この計算を追ってみます。

(1) 初期設定
例のところである、D(A,CPU) = 1はクラスAでノードCPUでの平均サービス時間が1となります。同様にD(A,DISK) = 3、D(B,CPU) = 2, D(B,DISK) = 4となります。実行可能母集団として、N(A) = 1、N(B) = 1としています。アルゴリズムの最初のループで(0A,0B)の時のQ(A,CPU),Q(A,Disk),Q(B,CPU),Q(B,Disk)を0としています。

(2) ループ：n = 1の時
この場合の実行可能母集団は(A,B)=(1,0),(0,1)となります。
(i) (1A,0B)の時
[滞在時間の計算]
R(A,CPU,(1,0)) = D(A,CPU)(1+A(A,CPU,(1,0)))
=D(A,CPU)(1+Q(CPU,(1-1,0))
Q(CPU)は、Q(CPU) = Q(A,CPU)+Q(B,CPU)

=D(A,CPU)(1+Q(A,CPU,(0,0)) +Q(B,CPU,(0,0)))
=1*(1+0+0)=1
となる。同様に
R(A,Disk,(1,0)) = D(A,Disk)(1+A(A,Disk,(1,0)))
=D(A,Disk)(1+Q(CPU,(1-1,0)))
=D(A,Disk)(1+Q(A,Disk,(0,0))+Q(B,Disk,(0,0)))
=3(1+0+0)=3
R(B,CPU),R(B,Disk)はnの部分が-1になるので、計算しません。

[スループットの計算]
X(A) = 1 / (R(A,CPU,(1,0)) + R(A,Disk,(1,0))) = 1 / (1 + 3) = 1/4
X(B)はR(B,*)がないので計算しません。

[系内人数の計算]
Q(A,CPU,(1,0)) = X(A) * R(A,CPU,(1,0)) = 1/4 * 1 = 1/4
Q(A,Disk,(1,0)) = X(A) * R(A,Disk,(1,0)) = 1/4 * 3 = 3/4
Q(B,CPU,(1,0)) = X(B) * R(B,CPU,(1,0)) = 0 (前の結果、更新なし?)
Q(B,Disk,(1,0)) = X(B) * R(B,Disk,(1,0)) = 0 (前の結果、更新なし?)

[注意]
アルゴリズムでは、Qm(1,0)つまり、Q(CPU,(1,0))とQ(Disk,(1,0))を求めることになっています。クラスを包含した結果を求めていますが、上のようにそれぞれを求めていきます。アルゴリズムの式を一応上げておくと下のようになります。
Q(CPU,(1,0)) = X(A) * R(A,CPU,(1,0)) + X(B) * R(B,CPU,(1,0))
= 1/4 * 1 + NULL = 1/4
Q(Disk,(1,0)) = X(A) * R(A,Disk,(1,0)) + X(B) * R(B,Disk,(1,0))
= 1/4 * 3 + NULL = 3/4

(ii) (0A,1B)
[滞在時間]
R(B,CPU,(0,1)) = D(B,CPU)(1+A(B,CPU,(0,1)))
=D(B,CPU)(1+Q(CPU,(0,1-1))
=D(B,CPU)(1+Q(A,CPU,(0,0)) +Q(B,CPU,(0,0)))
=2(1+0+0)
=2

R(B,Disk,(0,1)) = D(B,Disk)(1+A(B,Disk,(0,1)))
=D(B,Disk)(1+Q(Disk,(0,1-1))
=D(B,Disk)(1+Q(A,Disk,(0,0)) +Q(B,Disk,(0,0)))
=4(1+0+0)
=4

[スループットの計算]
X(B) = 1 / (R(B,CPU,(0,1)) + R(B,Disk,(0,1))) = 1 / (2 + 4) = 1/6

[系内人数の計算]
Q(A,CPU,(0,1)) = X(A) * R(A,CPU,(0,1)) = 0 (前の結果、更新なし?)
Q(A,Disk,(0,1)) = X(A) * R(A,Disk,(0,1)) = 0 (前の結果、更新なし?)
Q(B,CPU,(0,1)) = X(B) * R(B,CPU,(0,1)) = 1/6*2=1/3
Q(B,Disk,(0,1)) = X(B) * R(B,Disk,(0,1)) = 1/6*4=2/3

(iii) (1A,1B)
[滞在時間]
R(A,CPU,(1,1)) = D(A,CPU)(1+A(A,CPU,(1,1)))
=D(A,CPU)(1+Q(CPU,(1-1,1))
=D(A,CPU)(1+Q(A,CPU,(0,1)) +Q(B,CPU,(0,1)))
=1(1+0+1/3)
=4/3

R(A,Disk,(1,1)) = D(A,Disk)(1+A(A,Disk,(1,1)))
=D(A,Disk)(1+Q(Disk,(1-1,1))
=D(A,Disk)(1+Q(A,Disk,(0,1)) +Q(B,Disk,(0,1)))
=3(1+0+2/3)
=5

R(B,CPU,(1,1)) = D(B,CPU)(1+A(B,CPU,(1,1)))
=D(B,CPU)(1+Q(CPU,(1,1-1))
=D(B,CPU)(1+Q(A,CPU,(1,0)) +Q(B,CPU,(1,0)))
=2(1+1/4+0)
=5/2

R(B,Disk,(1,1)) = D(B,Disk)(1+A(B,Disk,(1,1)))
=D(B,Disk)(1+Q(Disk,(1,1-1))
=D(B,Disk)(1+Q(A,Disk,(1,0)) +Q(B,Disk,(1,0)))
=4(1+3/4+0)
=7

[スループットの計算]
X(A) = 1 / (R(A,CPU,(1,1)) + R(A,Disk,(1,1))) = 1 / (4/3 + 5) = 3/19
X(B) = 1 / (R(B,CPU,(1,1)) + R(B,Disk,(1,1))) = 1 / (5/2 + 7) = 2/19

[系内人数の計算]
Q(A,CPU,(1,1)) = X(A) * R(A,CPU,(1,1)) = 3/19*4/3=4/19
Q(A,Disk,(1,1)) = X(A) * R(A,Disk,(1,1)) = 3/19*5=15/19
Q(B,CPU,(1,1)) = X(B) * R(B,CPU,(1,1)) = 2/19*5/2=5/19
Q(B,Disk,(1,1)) = X(B) * R(B,Disk,(1,1)) = 2/19*7=14/19

このように求めていく。
関数体系として
getQ(int node, int class, int [] number){
(i)系内時間
R(A,CPU,(1,1)) = D(A,CPU)(1+A(A,CPU,(1,1)))
=D(A,CPU)(1+Q(CPU,(1-1,1))
=D(A,CPU)(1+Q(A,CPU,(0,1)) +Q(B,CPU,(0,1))) //ここで再帰呼び出し getQ(CPU,A,[0,1])
(ii)スループット
(iii)系内人数
}
のように作成していく

ソースコード

今回、それ以前の値を利用して、最新の値を求めていかなくてはならなかったので、プログラムとしてどのような動きになるかを確かめた。
例として、C=2(クラス数), N=10(全客数),Nベクトル=(5,5)(クラス人数の分布)、K=3(拠点数)としたときの流れである。

(1) [Lの初期化](Lは上記ではQ)
L[拠点番号][c1クラス数][c2クラス数]
L[0][0][0]=0
L[1][0][0]=0
L[2][0][0]=0

(2) n = 1 組み合わせ(1,0)
(i) w[c][k][n1][n2]
w[0][0][1][0] = D11 (1+L[0][0][1-1][0]+L[1][0][1-1][0]) = D11(1+L[0][1-1][0])=D11(1+L[0][0][0])=D11
w[0][1][1][0] = D12(1+L[0][1][1-1][0]+L[1][1][1-1][0]) = D11(1+L[1][1-1][0])=D12
w[0][2][1][0] = D13(1+L[0][2][1-1][0]+L[1][2][1-1][0]) = D13(1+L[2][1-1][0])=D13
w[1][0][1][0] = D21(1+L[0][0][1][0-1]+L[1][0][1][0-1]) = D21(1+L[0][1][0-1])=D21 負の要素番号の場合は、L=0として扱う
w[1][1][1][0] = D22(1+L[0][1][1][0-1]+L[1][1][1][0-1]) = D22(1+L[1][1][0-1])=D22
w[1][2][1][0] = D23(1+L[0][2][1][0-1]+L[1][2][1][0-1]) = D23(1+L[2][1][0-1])=D23

(ii) lambda[c][n1][n2]
lambda[0][1][0] = N[0] / (w[0][0][1][0] + w[0][1][1][0] + w[0][2][1][0] ) = N[0] / (D11 + D12 + D13 )
lambda[1][1][0] = N[1] / (w[1][0][1][0] + w[1][1][1][0] + w[1][2][1][0] ) = N[1] / (D21 + D22 + D23 )

(iii) L[k][n1][n2]
L[0][1][0] = lambda[0][1][0] * w[0][0][1][0] + lambda[1][1][0] * w[1][0][1][0] = N[0] / ( D11 + D12 + D13 ) * D11 + N[1] / (D21 + D22 + D23 ) * D21
L[1][1][0] = lambda[0][1][0] * w[0][1][1][0] + lambda[1][1][0] * w[1][1][1][0] = N[0] / ( D11 + D12 + D13 ) * D12 + N[1] / (D21 + D22 + D23 ) * D22
L[2][1][0] = lambda[0][1][0] * w[0][2][1][0] + lambda[1][1][0] * w[1][2][1][0] = N[0] / ( D11 + D12 + D13 ) * D13 + N[1] / (D21 + D22 + D23 ) * D23

(3) n = 1 組み合わせ (0,1)
w[0][0][0][1] = D11 (1+L[0][0][0-1][1]+L[1][0][0-1][1]) = D11(1+L[0][0-1][1])=D11(1+L[0][0-1][0])=D11
w[0][1][0][1] = D12(1+L[0][1][0-1][1]+L[1][1][0-1][1]) = D11(1+L[1][0-1][1])=D12
w[0][2][0][1] = D13(1+L[0][2][0-1][1]+L[1][2][0-1][1]) = D13(1+L[2][0-1][1])=D13
w[1][0][0][1] = D21(1+L[0][0][0][1-1]+L[1][0][0][1-1]) = D21(1+L[0][0][1-1])=D21 負の要素番号の場合は、L=0として扱う
w[1][1][0][1] = D22(1+L[0][1][0][1-1]+L[1][1][0][1-1]) = D22(1+L[1][0][1-1])=D22
w[1][2][0][1] = D23(1+L[0][2][0][1-1]+L[1][2][0][1-1]) = D23(1+L[2][0][1-1])=D23

(ii) lambda[c][n1][n2]
lambda[0][0][1] = N[0] / (w[0][0][0][1] + w[0][1][0][1] + w[0][2][0][1] ) = N[0] / (D11 + D12 + D13 )
lambda[1][0][1] = N[1] / (w[1][0][0][1] + w[1][1][0][1] + w[1][2][0][1] ) = N[1] / (D21 + D22 + D23 )

(iii) L[k][n1][n2]
L[0][0][1] = lambda[0][0][1] * w[0][0][0][1] + lambda[1][0][1] * w[1][0][0][1] = N[0] / ( D11 + D12 + D13 ) * D11 + N[1] / (D21 + D22 + D23 ) * D21
L[1][0][1] = lambda[0][0][1] * w[0][1][0][1] + lambda[1][0][1] * w[1][1][0][1] = N[0] / ( D11 + D12 + D13 ) * D12 + N[1] / (D21 + D22 + D23 ) * D22
L[2][0][1] = lambda[0][0][1] * w[0][2][0][1] + lambda[1][0][1] * w[1][2][0][1] = N[0] / ( D11 + D12 + D13 ) * D13 + N[1] / (D21 + D22 + D23 ) * D23

(4) n = 2 組み合わせ(1,1)
w[0][0][1][1] = D11 ( 1 + L[0][1-1][1] ) = D11 ( 1 + L[0][0][1] )
w[0][1][1][1] = D12 ( 1 + L[1][1-1][1] ) = D11 ( 1 + L[1][0][1] )
w[0][2][1][1] = D13 ( 1 + L[2][1-1][1] ) = D11 ( 1 + L[2][0][1] )
w[1][0][1][1] = D21 ( 1 + L[0][1][1-1] ) = D11 ( 1 + L[0][1][0] )
w[1][1][1][1] = D22 ( 1 + L[1][1][1-1] ) = D11 ( 1 + L[1][1][0] )
w[1][2][1][1] = D23 ( 1 + L[2][1][1-1] ) = D11 ( 1 + L[2][1][0] )

lambda[0][1][1] = N[0] / (w[0][0][1][1] + w[0][1][1][1] + w[0][2][1][1] )
lambda[1][1][1] = N[1] / (w[1][0][1][1] + w[1][1][1][1] + w[1][2][1][1] )

L[0][1][1] = lambda[0][1][1] * w[0][0][1][1] + lambda[1][1][1] * w[1][0][1][1]
L[1][1][1] = lambda[0][1][1] * w[0][1][1][1] + lambda[1][1][1] * w[1][1][1][1]
L[2][1][1] = lambda[0][1][1] * w[0][2][1][1] + lambda[1][1][1] * w[1][2][1][1]

以降繰り返し

D(c,k)の取り扱い
D(c,k) = α(c,k) * μ(c,k)
で計算する。D(c,k)は利用率に相当

次やること
クラス間移動なしはできたので、クラス間移動ありの場合でも確認する
応用モデルを考える
プログラムの改良(クラスが増えた場合)