Jacksonネットワークのシミュレーションをします。シミュレーションの流れは次のようになります。
(1) 初期設定(ネットワーク構成:推移確率行列、到着率、サービス率)
int N = 3;
double p[][] = {{0,0,1},{0,0,0.6},{0.5,0,0}};
double lambda[] = {2,1,0};
double mu[] = {5,4,6};
(2) 最初の到着、サービスの設定
サービスは到着の後に起こるので、サービス時間に到着までの時間を足してあります。
for(int i = 0; i < N; i++) {
arrival[i] = this.getExponential(lambda[i]);
service[i] = arrival[i] + this.getExponential(mu[i]);
}
(3) シミュレーション実施(イベント算出)
シミュレーションはイベントが発生した時に処理を実行するイベントドリブン型で実施します。次のイベントまで発生する時間が何かを到着と退去にわけて、イベント発生までの最小時間を算出します。
//最小値の算出(次のイベント)
min_arrival = 100;
min_service = 100;
arrival_index = N;
service_index = N;
for(int i = 0; i < N; i++) {
if(min_arrival > arrival[i]) {
min_arrival = arrival[i];
arrival_index = i;
}
if(min_service > service[i]) {
min_service = service[i];
service_index = i;
}
}
(4) シミュレーション(到着の場合)
到着と退去に分けて算出したイベントまでの最小時間を比較して、到着と退去のどちらが起きたかを確認します。到着が起きた場合は下記のようになります。
total_queueは延べ系内人数、queueは現在の系内人数です。queue*min_arrivalでそれまでの平均系内人数を延べ系内人数に加えています。service、arrivalで時間を進め、到着があったノードには次の到着の時間を設定します。
if( min_arrival < min_service ) { //外部or内部到着が発生
System.out.println("Arrival");
for(int i = 0; i < N; i++) {
total_queue[i] += queue[i] * min_arrival; //延べ系内人数
service[i] -= min_arrival;
arrival[i] -= min_arrival;
}
queue[arrival_index]++; //現在の系内人数
arrival[arrival_index] = this.getExponential(lambda[arrival_index]);
elapse += min_arrival;
System.out.println("Index = "+ arrival_index);
for(int i = 0; i < N; i++) queuelength[i].add(queue[i]);
}
(5) シミュレーション(退去の場合)
退去の場合も到着とほとんど同じですが、退去後、まだ待ち人数が1以上ならサービス時間のみ設定しますが、待ち人数が0ならば、サービス時間に到着時間を加えておきます。
else if(min_arrival >= min_service ){ //退去が発生
System.out.println("Departure");
for(int i = 0; i < N; i++) {
total_queue[i] += queue[i] * min_service; //延べ系内人数
arrival[i] -= min_service;
service[i] -= min_service;
}
queue[service_index]--; //現在の系内人数
if(queue[service_index] > 0) service[service_index] = this.getExponential(mu[service_index]);
else service[service_index] = arrival[service_index] + this.getExponential(mu[service_index]);
elapse += min_service;
System.out.println("Index = "+ service_index);
for(int i = 0; i < N; i++) queuelength[i].add(queue[i]);
さらに退去する客は、外部に退去する場合もありますが他のノードに移動する場合があります。推移確率行列からどこに移動するかを決定します。
//退去する客の行き先決定
sum_p = 0;
dep_p = rnd.nextDouble();
dep_index = N;
for(int i = 0; i < N; i++) {
sum_p += p[service_index][i];
if(dep_p < sum_p) {
dep_index = i;
break;
}
}
if(dep_index != N) {
arrival[dep_index] = 0; //行き先は到着が発生, Nの時は外部へ退去
service[dep_index] = this.getExponential(mu[dep_index]); //サービス時間を再設定
}
注意ですが、他のノードに移動する場合、移動時間は0で移動し、次のイベントとしてすぐに到着が発生します。行き先のノードでの到着時間を0として、その客のサービス時間を再設定します。
(6) 平均系内人数の算出
延べ系内人数をシミュレーション時間で割り、平均系内人数を算出します。
for(int i = 0; i < N; i++) total_queue[i] = total_queue[i] / time;
return total_queue;
シミュレーション時間を10000とすると、下記の結果が得られます。
Simulation : 系内人数 = [10.020321642112451, 0.32465843337530376, 5.849312685287813]
シミュレーション時間を増やせば精度も上がっていきます。
(7) グラフ描画
時系列で変化する情報をグラフに描いていきます。例えば系内人数の動きはグラフのようになります。
(7) ソースコード
https://github.com/smzn/Jackson_Simulation
2019年8月15日木曜日
2019年8月10日土曜日
待ち行列:Jacksonネットワークの計算
待ち行列理論の中でJacksonネットワークの計算をしてみます。
(1)推移確率行列
ノードがN個あるネットワークを考える。推移確率行列は次の式になります。通常、行和が1になりますが、1にならない場合がそのノードから外部への退去と考えます。例えばi行の行和が0.7の場合、ノードiから外部への退去は0.3となります。
推移確率行列が上の場合、ノード2は行和が0.6のため、ノード2から外部への退去は0.4となります。
(2) トラフィック方程式
次にノードiへの到着率αを求めます。トラフィック方程式は次のようになります。
ベクトル表記
要素で表記
行列表記
これを式変形して計算しやすい形にします。
ベクトル表記
この線形方程式を解くことでαが得られます。ρは次の式で求められます。μはサービス率です。
(3) 定常分布の算出
ノード i の人数を ni とすると、全ノードのそれぞれの人数のベクトル表現を n = (n1, n2, · · · , nN ) とする。 この時、定常分布は下記になる。ただし全てのノードで ρi < 1 となることが必要である。周辺分布は、それぞれのノードの定常分布となる
(4) 平均系内人数
各ノードの平均系内人数は M/M/1 待ち行列モデルに従うので、それぞれの平均系内人数 Li は下記になる
(5) 数値計算
簡単な例で数値計算をしてみます。加藤、小沢、ORの基礎、実教出版、P153の例を使ってみます。
この場合の推移確率行列は、下記です。
外部からの到着率は
サービス率は
これからトラフィック方程式を解くと、αは
これからρは
平均系内人数は
ρが1に近くなると系内人数は発散します。
定常分布を求めます。今回はn = [6, 1, 2]と人数分布になる場合のそれぞれのノードの定常分布は、下記になります。
n = [6, 1, 2]となる場合の結合分布は
となります。各ノードの定常分布のグラフは下記になります。
[ソースコード:Matlab]
(1)推移確率行列
ノードがN個あるネットワークを考える。推移確率行列は次の式になります。通常、行和が1になりますが、1にならない場合がそのノードから外部への退去と考えます。例えばi行の行和が0.7の場合、ノードiから外部への退去は0.3となります。
推移確率行列が上の場合、ノード2は行和が0.6のため、ノード2から外部への退去は0.4となります。
(2) トラフィック方程式
次にノードiへの到着率αを求めます。トラフィック方程式は次のようになります。
ベクトル表記
要素で表記
行列表記
これを式変形して計算しやすい形にします。
ベクトル表記
この線形方程式を解くことでαが得られます。ρは次の式で求められます。μはサービス率です。
(3) 定常分布の算出
ノード i の人数を ni とすると、全ノードのそれぞれの人数のベクトル表現を n = (n1, n2, · · · , nN ) とする。 この時、定常分布は下記になる。ただし全てのノードで ρi < 1 となることが必要である。周辺分布は、それぞれのノードの定常分布となる
(4) 平均系内人数
各ノードの平均系内人数は M/M/1 待ち行列モデルに従うので、それぞれの平均系内人数 Li は下記になる
(5) 数値計算
簡単な例で数値計算をしてみます。加藤、小沢、ORの基礎、実教出版、P153の例を使ってみます。
この場合の推移確率行列は、下記です。
外部からの到着率は
サービス率は
これからトラフィック方程式を解くと、αは
これからρは
平均系内人数は
ρが1に近くなると系内人数は発散します。
定常分布を求めます。今回はn = [6, 1, 2]と人数分布になる場合のそれぞれのノードの定常分布は、下記になります。
n = [6, 1, 2]となる場合の結合分布は
となります。各ノードの定常分布のグラフは下記になります。
[ソースコード:Matlab]
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
%パラメタ設定 | |
N = 3; | |
P = [0 0 1; 0 0 0.6; 0.5 0 0] | |
lambda = [2; 1; 0] | |
mu = [5; 4; 6] | |
%[p^t-E]α = -λ | |
A = P.' - eye(N) | |
alpha = linsolve(A, -lambda) | |
rho = alpha ./ mu | |
%平均系内人数 | |
L = rho ./ (1 - rho) | |
fplot(@(x) x/(1-x), [0 0.97]) | |
%人数分布 | |
n = [6, 1, 2] | |
%nとした場合の定常分布(周辺)を求める | |
for i = 1:N | |
i | |
rho(i) | |
n(i) | |
pi(i) = (1-rho(i))*rho(i)^n(i); | |
end | |
pi | |
%rho(i)を使い、人数を変化させた時の定常分布のグラフ | |
figure | |
hold on | |
for i = 1:N | |
fplot(@(x) (1-rho(i))*rho(i)^x, [0 30]) | |
end | |
legend | |
hold off | |
%Jackson定常分布 | |
jpi = 1; | |
for i = 1:N | |
jpi = jpi * pi(i); | |
end | |
jpi |
登録:
投稿 (Atom)