課題2: sample2.cpp — M/M/1 待ち行列

構成

負荷は \( \rho = \lambda / \mu \)。sample2.cpp 既定では \(\mu = 1\) 固定、\(\rho\) を 0.1〜0.9 まで 0.05 刻みで変化させる。

ソースコード（配布版）

下のテキストエリアには 配布/sample2.cpp の中身がそのまま貼り付けてある。// *** の箇所を埋めて、 「sample2.cpp をダウンロード」で自分の環境に保存できる。

実験手順

1. 準備 上のテキストエリア（または 配布/sample2.cpp）の // *** を埋める:
   • server::departure / server::arrive 冒頭: r_total の積分更新
     \(r_\text{total} \mathrel{+}= (t_c - t_\text{current}) \cdot \mathrm{queue}\)
   • server::service の dist==3（指数）: \(t_\text{finish} = t_c - t_\text{ave}\ln(1-u),\; u = \texttt{ran2(idum)}\)
   • server::service の dist==4（対数正規）: 平均 t_ave・分散 2*t_ave*t_ave の設定で 正規パラメタを計算して exp する
   • server::cal_average_wait: 各パケットの滞在時間を総和 \(\texttt{total\_W} \mathrel{+}= \texttt{wait\_time[i][1]} - \texttt{wait\_time[i][0]}\)
2. 実行 rho を 0.1〜0.9 まで 0.05 刻み（17 点） で変えて毎回ビルド・実行、出力を 1 つのファイルに追記:

   ./a.out >> results.txt

   17 点ぶんの rho = … / L = … / W = … が連結される。
3. 確認 専用ツール (queueing_result.html) で results.txt を読み込み、 サービス分布（指数 / 対数正規）を選んで「集計」。L–ρ, W–ρ のグラフで理論カーブとシミュ点が重なるか確認する。
   学習ポイント: Little の式 \(L = \lambda W\) が成立するかも併せて確認（別々に求めた値の整合性チェック）。 高負荷（ρ ≥ 0.85）では観測時間不足で W が理論より大きく出る傾向がある — これもレポートで触れる価値あり。
4. オプション ソース先頭の #define TRACE_QUEUE 1 にして 1 回実行すると queue_trace.csv が出力される。 同 HTML の「4. バッファ状態」で Q(t) の階段グラフ + パケットがバッファに溜まるアニメーション を確認できる。 ρ=0.5 と ρ=0.9 で比較すると高負荷時の詰まり方が視覚的に分かる。

出力形式と HTML との連携

rho = 0.5, lambda = 0.5, mu = 1
Time-average Packet Number L = 1.00394
Average Sojourn Time  W = 2.07682

1 回の実行で 1 セット出力される。rho を変えて複数回実行し、1 ファイルに追記:

./a.out >> results.txt

queueing_result.html が自動で全 ρ を抽出して L–ρ, W–ρ を理論カーブと重ね描きする。

主要な定数

IDUM	高品質乱数生成器 ran2 の初期値（負のシード）。
END_TIME	シミュレーション終了時刻（例 100,000）。大きいほど統計が安定。
DIST_TYPE	サービス時間の分布。3:指数, 4:対数正規。
rho, mu	負荷・サービス率。lambda = rho * mu が到着率。
MAX_NUM	滞在時間を記録できるパケット数の上限（配列サイズ）。

イベントドリブン方式

時間を 1 単位ずつ進めるタイムドリブンに対し、イベントドリブンは「次に何かが起きる時刻までジャンプ」する方式。待ち行列のように離散イベントの系では、無駄な処理がほぼなくなる。

// メインループ (擬似コード)
while (top != NULL) {
  top->check();              // 最先頭イベント（離脱）を処理
  top = top->get_next();     // 次のイベントへ
  if (top->get_time() > END_TIME) break;
}

scheduler クラス

終了時刻 time でソート済みのリンクトリスト。add() でイベントを時刻順に挿入、check() で先頭の処理を呼ぶ、get_next() で次のスケジュールを返す。

server クラス

• arrive(t): 到着イベント処理。バッファに入れ、必要なら service() を起動。
• service(t): サービス開始。サービス時間を発生させ、t_finish 時刻に離脱イベントをスケジュールに追加。
• departure(t): 離脱イベント処理。order==1（ルータ）なら総数カウント、queue==0 で終了 or 次のパケットをサービス。

サービス時間の生成 埋める

service() 内で分布に応じてサービス時間を発生させる。dist==3 の指数分布は逆関数法で:

\[ \text{service time} = -t_{\text{ave}} \ln(1 - u), \quad u = \texttt{ran2(idum)} \]

dist==4 の対数正規は、sample1.cpp の lognormal_dist と同じロジック。このコード中では平均が t_ave、分散は 2*t_ave*t_ave に設定する規約。

void server::service(double t_c){
    t_current = t_c;
    if(dist == 1) t_finish = t_c + t_ave;              // 一定
    else if(dist == 2) t_finish = t_c + t_ave*2*ran2(idum); // 一様
    else if(dist == 3){
      //*** 次のイベント時間 t_finish を設定．時間間隔は指数分布
    }
    else if(dist == 4){
      //*** 次のイベント時間 t_finish を設定．時間間隔は対数正規分布
    }
    ...
}

統計量の計算 埋める

平均パケット数 L

時間平均 \( L = \tfrac{1}{T}\int_0^T Q(t)\,dt \) を r_total に積分で蓄積する。 キュー長が変わるのは到着時と離脱時のみなので、その直前に「前回時刻から現時刻までの区間の queue * dt」を足しこめばよい。

void server::departure(double t_c){
  if(t_current > begin_time){
    // *** r_total を更新（t_current から t_c の間は queue は一定）
    // **** r_total += (t_c - t_current)*queue; が典型例
  }
  ...
}

滞在時間 W

各パケットの到着時刻と離脱時刻を wait_time[i][0], wait_time[i][1] に記録しておき、差の総和をパケット数で割る。

double server::cal_average_wait(){
  double total_W = 0.0;
  for(int i = 0; i < dep_num; i++){
    //*** total_W に i 番目のパケットの滞在時間を足しこむ
    //**** wait_time[i][0]: 到着時刻, wait_time[i][1]: 離脱時刻
  }
  return total_W / dep_num;
}

理論式

M/M/1 (指数サービス)

\[ L = \frac{\rho}{1-\rho}, \qquad W = \frac{1}{\mu(1-\rho)} \]

M/G/1 (一般サービス・Pollaczek–Khinchine 式)

サービス時間の変動係数の 2 乗を \(C^2 = \mathrm{Var}(S)/\mathbb{E}[S]^2\) とおくと:

\[ L_q = \frac{\rho^2 (1+C^2)}{2(1-\rho)}, \quad L = \rho + L_q, \quad W = L/\lambda \]

sample2 の対数正規サービスは「平均 \(t_{\rm ave}\)・分散 \(2\, t_{\rm ave}^2\)」という規約なので:

\[ C^2 = \frac{2\, t_{\rm ave}^2}{t_{\rm ave}^2} = 2 \quad\Rightarrow\quad L = \rho + \frac{3\rho^2}{2(1-\rho)} \]

指数分布は \(C^2 = 1\) なので M/M/1 の式と一致する（整合性チェックになる）。

ran2() — 高品質な一様乱数

sample2, sample3 で使用される ran2(long *idum) は Numerical Recipes の手法で、周期が約 \(2\times 10^{18}\) と非常に長く、統計的品質の高い一様乱数を返す。長時間シミュレーションでも偏りが出にくい。

初期値 IDUM は負の整数で与える（内部状態のリセットのトリガー）。このシード固定により実験の再現性を確保する。

sample1 の学習用 LCG (周期 \(\le 65536\)) とは役割分担していることに注意。

埋める箇所のまとめ ***

• server::departure 冒頭: r_total += (t_c - t_current)*queue（Q(t) 積分）
• server::arrive 冒頭: 同上（到着側でも積分を更新）
• server::service の dist==3: 指数分布のサービス時間 t_finish = t_c + (-t_ave)*log(1 - ran2(idum))
• server::service の dist==4: 対数正規の正規パラメタを計算して normal_dist 相当を生成し、exp してから t_c に加算
• server::cal_average_wait: total_W += wait_time[i][1] - wait_time[i][0]