待ち行列のシミュレーション｜シミュレーション

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：離散事象シミュレーションとは　|　理論：M/M/1 待ち行列モデル（オペレーションズ）

要点（BLUF）

SimPy を使うと、待ち行列の離散事象シミュレーションを「客＝プロセス、窓口＝Resource」として直感的に書けます。
到着とサービスを指数分布で生成した M/M/1 シミュは、平均待ち時間が理論 $W_q = \rho/(\mu(1-\rho))$ と一致します。
利用率 $\rho$ を上げると待ち時間が急増し、 $\rho \to 1$ で爆発します（ $\rho=0.9$ で $W_q=9$ ）。待ち行列の理論はオペレーションズへ。

1. SimPy でM/M/1を書く

前トピックでは FEL を手書きしましたが、SimPy なら本質だけ書けます。客をジェネレータ関数にし、with server.request() で窓口を取り合い、yield env.timeout(...) で時間を経過させます。

import numpy as np
import simpy

def run_mm1(lam, mu, sim_time, seed):
    rng = np.random.default_rng(seed)
    env = simpy.Environment()
    server = simpy.Resource(env, capacity=1)   # 窓口1つ
    waits = []

    def customer(env, server):
        arrive = env.now
        with server.request() as req:          # 窓口を要求（空くまで待つ）
            yield req
            waits.append(env.now - arrive)     # 待ち時間を記録
            yield env.timeout(rng.exponential(1/mu))   # サービス時間

    def source(env, server):                   # 客の到着源
        while True:
            yield env.timeout(rng.exponential(1/lam))  # 到着間隔
            env.process(customer(env, server))

    env.process(source(env, server))
    env.run(until=sim_time)
    return np.array(waits)

lam, mu = 0.8, 1.0
waits = run_mm1(lam, mu, 200_000, seed=42)
rho = lam / mu
print(f"処理客数 = {len(waits)}")
print(f"平均待ち時間（シミュ） = {waits.mean():.3f}")
print(f"M/M/1 理論 Wq = {rho/(mu*(1-rho)):.3f}")
print(f"M/M/1 理論 Lq = {rho**2/(1-rho):.3f}（行列内の平均人数）")

出力：

処理客数 = 160257
平均待ち時間（シミュ） = 4.138
M/M/1 理論 Wq = 4.000
M/M/1 理論 Lq = 3.200（行列内の平均人数）

出力の意味： $\lambda=0.8, \mu=1.0$ （ $\rho=0.8$ ）で平均待ち時間 4.138、理論 4.000 とよく一致（差はシミュレーション誤差と過渡の影響、出力解析（過渡・定常・複製））。M/M/1 の主要公式は次の通りで、導出はオペレーションズに譲ります（ここは実装と検証が役割）：

W_q = \frac{\rho}{\mu(1-\rho)},\qquad L_q = \frac{\rho^2}{1-\rho},\qquad \rho = \frac{\lambda}{\mu}

2. 利用率を上げると待ち時間が爆発する

待ち行列の最も重要な性質は、利用率 $\rho$ が1に近づくと待ち時間が発散すること。シミュレーションで確かめます。

import numpy as np
import simpy

def run_mm1(lam, mu, sim_time, seed):
    rng = np.random.default_rng(seed)
    env = simpy.Environment()
    server = simpy.Resource(env, capacity=1)
    waits = []
    def customer(env, server):
        arrive = env.now
        with server.request() as req:
            yield req
            waits.append(env.now - arrive)
            yield env.timeout(rng.exponential(1/mu))
    def source(env, server):
        while True:
            yield env.timeout(rng.exponential(1/lam))
            env.process(customer(env, server))
    env.process(source(env, server))
    env.run(until=sim_time)
    return np.array(waits)

mu = 1.0
for rho in [0.5, 0.7, 0.9]:
    waits = run_mm1(rho, mu, 300_000, seed=7)   # lam = rho（mu=1）
    print(f"rho={rho}: シミュ Wq={waits.mean():.3f}  理論 Wq={rho/(1-rho):.3f}")

出力：

rho=0.5: シミュ Wq=0.986  理論 Wq=1.000
rho=0.7: シミュ Wq=2.379  理論 Wq=2.333
rho=0.9: シミュ Wq=9.709  理論 Wq=9.000

出力の意味（ $\mu=1$ なので $W_q = \rho/(1-\rho)$ ）： $\rho$ が $0.5 \to 0.7 \to 0.9$ と上がると、待ち時間は $1 \to 2.33 \to 9$ と非線形に急増。 $\rho=0.9$ では既に9倍。 $\rho \to 1$ で $1/(1-\rho)$ が発散します。これが「サーバをギリギリまで稼働させると待ちが爆発する」という、容量設計の核心的直観です。 $\rho=0.9$ でシミュ値が理論よりやや高いのは、混雑時ほど過渡（warm-up）の影響と分散が大きいため（出力解析（過渡・定常・複製））。

3. 拡張：複数窓口・有限容量

SimPy なら capacity=c にするだけで M/M/c（窓口 $c$ 個）になり、有限の待合室や優先度も PriorityResource などで表現できます。「窓口を増やすべきか、速くすべきか」といった設計問題は、理論（M/M/c）と DES の両輪で解きます。理論が効かない複雑な系（一般分布のサービス時間 G/G/c、割り込み、故障）こそ DES の出番です。

数式の直観的意味

$W_q = \rho/(\mu(1-\rho))$ の $1/(1-\rho)$ が爆発の源です。 $1-\rho$ は「サーバが空いている割合」。これがゼロに近づくと、新しい客が来ても前の客が片付く隙間がほとんどないため、行列が積み上がる。待ち時間が $\rho$ に比例ではなく $1/(1-\rho)$ で増えるのは、混雑が混雑を呼ぶ正のフィードバック。シミュレーションがこの非線形を素直に再現するのは、指数到着・指数サービスのメモリレス性（代表的分布の生成（正規・指数・ポアソン））が M/M/1 の仮定そのものだから。サービス時間のばらつきが大きいほど待ちは長くなり（G/G/1 では変動係数が効く）、そこが理論の限界＝シミュレーションの価値です。

⚠️ よくある誤解・落とし穴

「 $\rho < 1$ なら待ちは大したことない」ではない： $\rho=0.95$ で $W_q=19$ 、 $\rho=0.99$ で $W_q=99$ 。1に近いと激増します。
「シミュと理論がピタリ一致しないとバグ」ではない：有限時間・空状態スタートの過渡で数%ずれます。長時間化・warm-up 除去で縮みます（出力解析（過渡・定常・複製））。
「待ち時間とシステム滞在時間は同じ」ではない： $W_q$ （待ち）＋サービス時間が滞在時間 $W$ 。混同しないこと。
「 $\rho \ge 1$ でも定常値が出る」ではない： $\rho \ge 1$ は不安定で行列が無限に伸び、定常分布が存在しません。シミュは発散します。
「理論をここで導出すべき」ではない：待ち行列理論の導出はオペレーションズの領分。ここは実装と数値検証に徹します（境界）。

対応シミュレーション参照

本文の SimPy M/M/1（default_rng(42)）と利用率スイープ（default_rng(7)）。理論はM/M/1 待ち行列モデルへ。