待ち行列：Jacksonネットワークのシミュレーション

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