2019年9月22日日曜日

待ち行列:BCMP Network(1)

BCMPネットワークについて考えていきます。
BCMP ネットワークでは待ち行列ネットワークが拡張されています。開放型、閉鎖型の両方が表現でき、複数の客クラスを持てます。またサービスについても、FCFS(First Come First Served)、プロセッサーシェアリング、LCFS-PR(Last Come First Served with Preemption)、無限サーバがありますが、ここではFCFSを扱います。

(1) 今回のモデルの特徴をまとめると下記になります。

  • 閉鎖型待ち行列ネットワーク。ノード数:m、ネットワーク内客数:N
  • 複数の客クラスを持つ。クラス数:J
  • 各ノードでのサービスはFCFS。サービス時間はμi(ni)の指数分布に従う(全ての客クラスで共通)
(2) 推移確率の定義
ノード:i、クラス:μから、ノード:j、クラス:νに移動する確率を次のようにする。
クラスの集合Tを次のように定義しておく。
推移確率は次を満たす。

推移確率行列Rは上記になります。ノードへの移動の推移確率行列のブロックが、クラス推移に対してそれぞれ存在します。

(2) ノードiへのクラスμの客の総到着率:トラフィック方程式を解く(閉鎖型)
ノードiへのクラスμの客の総到着率を求めるために、トラフィック方程式 α = α R (閉鎖型)を解いていきます。ここでαはクラス別にm個のノードを並べたm * J 個のベクトルです。
トラフィック方程式の右辺は下の式になります。

これを計算すると、m * J個のベクトルができますが、それぞれの成分はαとRの列(To Class)の内積になります。簡単に書くと、次のようになります。

最初の第1成分だけ取り出して見ると(上の図の一番左側)、m個のベクトルができます。この時、(1) クラスを固定、(2) ノードkを動かして計算していきます。
To Class 1だけで考えると、上記のm個のベクトルができるので、これをTo Class 全てで考えると、m * J 個のベクトルができます。
つまり、ノードi、クラスμにおけるトラフィック方程式を成分で書くと、
この m * J 個の連立方程式を解くことになります。この場合、解は相対値で出てきますが、解を固定するために、α11 = 1とおきます。
を式変形すると、Rの転置行列を使って、次のように変形できます。
さらに、α11 = 1とすると、式変形ができますが、簡単に4行4列でやった場合が次になります。
つまり、転置した行列の第1行と第1列要素を取り除いた行列を作成し、右辺に転置行列の第1列の第2成分から全てを書きます。転置前の行列では赤の部分になります。


(3) トラフィック方程式を実際に求める
簡単に考えるために、客は2クラスで同じ推移確率で動くとします。ノード数は12です。クラス間の移動は無いものとします。
例えば推移確率行列を図のように作ってみました。

緑の部分はクラス1からクラス1へのノード間推移確率行列です。右上と左下のグレーの部分はそれぞれ、クラス1 -> クラス2、クラス2 -> クラス1へ移動する場合のノード間推移確率です。クラス間の移動が無いため、全て0となります。右下の黄色の部分はクラス2でのノード間推移確率になります。ノード数が12、クラス数が2のため12*2=24、24*24の推移確率行列になります。

次にトラフィック方程式用の行列を用意します。上でも書いた通り、Rを転置して、1行、1列の要素を全て取り除き、対角要素は-1したものです。(m * J -1) × (m * J -1)の大きさになります。
実は、このようにクラス間移動が無い場合、閉鎖型待ち行列ネットワークのモデルとして、うまくないです。このようにクラス間移動がない場合、閉鎖型待ち行列ネットワークが2個あるのと同じになります。ノードの平均系内人数は足し合わせたものになります。(これは今の僕の考え。変わるかも)
実際これをα11=1として解いてみてもクラス1は数値が求められますが、クラス2は全て0となってしまいます。クラス2のα12=1としないと求められないです。
エクセルファイル
[確認事項]
・クラス間の移動が無くても、大丈夫なはずなので、トラフィック方程式を確認(クラス間の移動がない場合、クラス2でもα12=1としないといけないか?式変形の問題なので計算してみる)

BCMPモデルとしてクラス間の移動があった場合の方がモデルとして適していると考えます。そこで次のような推移確率行列を定義します。
クラス1でノード5,10に来た客は退去するとクラス2に移動します(右上の赤い部分)。またクラス2でノード1,6,11に来た客はクラス1に移動します(左下の青い部分)。
推移確率行列csvファイル
これを転置して、対角要素から-1をしておきます。

試しにエクセルのソルバーで解いてみます。今回はα11-α122までの24個のベクトルがあり、次のようになります。
・初期値は全て1
・αベクトルと転置して-1した行列の行との内積をとります。
目的関数は|α11-1|としておきます。変数はαのベクトルです。
・α11=1
・各内積=0
プログラム(Java)で解いた解とほぼ同じになりました。Javaで解いた解は下記
ソルバー付きエクセルファイル
これでトラフィック方程式の解αが求められ、各ノードの到着率がわかりました。

2019年8月15日木曜日

待ち行列: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

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]

2019年7月25日木曜日

サーバ証明書の設定(Let's Encrypt)

Let's Encryptをサーバに設定します。
(1) 前提条件
・OSはCentOS7
・DNSサーバに登録済み
・Webサーバが稼働している(ユーザディレクトリも解放している)
・Virtualhostとしてmzncert.mizunolab.infoがつけられているとします。
http://mzncert.mizunolab.info/test.php  でアクセスすると、アクセス可能
 https://mzncert.mizunolab.info/test.php  ではアクセス出来ない
(2) Let's Encryptの導入
# yum install certbot certbot-apache
# certbot run --apache -d mzncert.mizunolab.info
いくつかの質問に答えます。
https://mzncert.mizunolab.info/test.php  にアクセスすると、httpsで見られるようになっています。