離散事象シミュレーション最適化｜オペレーションズマネジメント

🎓 レベル：応用　|　重要度：A（必須）

📎 前提：M/M/1 待ち行列モデル（Lindley 再帰の離散事象シミュ）・安全在庫と発注点（σ√L の安全在庫＝L 一定が前提）・確率的在庫モデル（新聞売り子問題）（臨界比 p/(p+h)）　|　確率の土台：正規分布（標準正規・標準化）（CLT・分位点）　|　本章の締め

要点（BLUF）

閉じた式が無いとき——複雑なネットワーク、非標準分布、確率的リードタイム、有限バッファ、相関や追い越し——は、シミュレーションで意思決定変数を最適化します。本テキストはプロセス分析とリトルの法則（リトルの法則を擬似系列で検証）・M/M/1 待ち行列モデル（Lindley 再帰）・ブルウィップ効果（多段直列の増幅）・リーン生産・JIT・かんばん（プッシュ vs プル）で擬似シミュを多用してきました。本稿はそれを**「シミュレーション最適化」という手法**へ昇格させる、第9章＝全テキストの集大成です。
シミュレーションの出力は**確率変数（ノイズ込みの推定値）です。だから最適化は「探索＋統計的配慮」になる。鉄則は3つ——(1) 複製して平均をとり、信頼区間 $\bar X\pm z_{0.975}\dfrac{s}{\sqrt m}$ で精度を示す、(2) ノイズの生の最大/最小をそのまま選ばない、(3) 2案の比較は共通乱数（CRN）**で分散を減らす。
CRN（common random numbers）は、比較したい複数案に同じ乱数列を使う技。共通のショックが相殺して差の分散が激減し、同じ計算量でどちらが良いかをはるかに鋭く判定できます。
例として確率的リードタイム＋確率需要の周期発注（基在庫）を扱います。リードタイムがランダムだと安全在庫と発注点の $\sigma\sqrt{L}$ は前提（ $L$ 一定）ごと崩れ、しかも注文が追い越し合うのでコスト最小の発注水準を閉形式で出せません。シミュレーションで最適 $S^*$ を探索します。

1. 閉形式が効かないとき：擬似シミュを「手法」に昇格する

これまでの章では、まず数式で解き、シミュレーションで裏取りしてきました。M/M/1 は $L=\rho/(1-\rho)$ という閉形式があり（M/M/1 待ち行列モデル）、安全在庫は $z\sigma_d\sqrt{L}$ という閉形式があり（安全在庫と発注点）、シミュは「式が正しい証拠」でした。

ところが現実の運用は、閉形式の前提を平気で破ります。

サービス時間が指数でない・到着が相関する一般の G/G/c の網
確率的リードタイム（ $L$ が確率変数）や、注文が追い越し合う補充
有限バッファでブロッキング／スターベーションが絡む直列ライン（リーン生産・JIT・かんばん）
需要が歪む・断続する非標準分布、複数 SKU の相関

こうなると閉形式は無い（あっても粗い近似）。残る道は、システムを計算機の中で動かして測る——離散事象シミュレーションです。そして「測る」だけでなく、意思決定変数（発注水準・バッファ量・要員数）を動かして最適点を探すのがシミュレーション最適化。本テキストで使ってきた擬似シミュは、実はこの手法の素材だったわけです。

flowchart LR
  PICK["意思決定変数の候補を選ぶ<br/>（例：基在庫水準 S）"] --> SIM["複製してシミュレート<br/>（同じ S を何本も回す）"]
  SIM --> EST["平均と信頼区間を推定<br/>Xbar ± z·s/sqrt(m)"]
  EST --> CMP{"案どうしを比較<br/>（CRNで分散を減らす）"}
  CMP -->|"差がCIで判定できない"| MORE["複製を増やす/CRN"]
  CMP -->|"差が有意"| NEXT["次の候補へ / 採択"]
  MORE --> SIM

2. シミュレーション最適化の作法：推定・複製・信頼区間

シミュレーションの出力（平均コスト・平均待ち時間）は1回ごとにばらつく確率変数です。乱数の種を変えれば違う値が出る。だから「1回回して一番良かった案」を選ぶのは危険——たまたまノイズで良く出ただけかもしれません。

正しくは、各案を** $m$ 回複製して平均 $\bar X$ をとり、その精度を信頼区間**で表します。複製は独立なので、中心極限定理（正規分布（標準正規・標準化））から平均の95%信頼区間は

\bar X\pm z_{0.975}\,\frac{s}{\sqrt{m}}\qquad(z_{0.975}=1.96,\ s=\text{複製の標本標準偏差})

複製を4倍にすると区間は半分（ $\sqrt m$ の効き）。「どれくらい精度があるか」を区間で示さないシミュ結果は、点推定だけの裸の数字で、比較の根拠になりません。

そしてもう一つの鉄則——ノイズの生の最大/最小をそのまま選ばない。多数の案を1回ずつ評価して「最小コストの案」を選ぶと、運悪く下振れした案を拾いがちで、選んだ案の真のコストは見かけより悪い（最適化バイアス／winner’s curse）。複製と信頼区間で「本当に差があるか」を確かめてから選びます。

3. 共通乱数（CRN）：比較の分散を減らす

2つの案 $S_a,S_b$ のどちらが良いかを知りたいとき、見たいのは差 $\text{cost}(S_b)-\text{cost}(S_a)$ です。差の分散は

\mathrm{Var}\bigl(\text{cost}_a-\text{cost}_b\bigr)=\mathrm{Var}_a+\mathrm{Var}_b-2\,\mathrm{Cov}(\text{cost}_a,\text{cost}_b)

独立な乱数で別々に回すと $\mathrm{Cov}=0$ で、差にはノイズが2案分のります。ところが同じ乱数列（同じ需要列・同じリードタイム列）で両案を回すと、両者は同じショックに同じように反応するので $\mathrm{Cov}>0$ 。第3項が効いて差の分散が縮みます。これが共通乱数（CRN, common random numbers）——「条件を揃えて比較する」対照実験の発想です。

直感は明快です。同じ需要の山谷を両案にぶつければ、「需要が荒れたから両方コスト高」という共通変動は差し引きで消え、方策の違いだけが残る。同じ計算量でも、CRN は比較を桁違いに鋭くします（第5節で差の標準偏差が約1/4になるのを見ます）。

4. 確率的リードタイム＋確率需要の基在庫を探索する（コード）

周期発注の基在庫（base-stock）方策：毎期、在庫ポジションを $S$ まで引き上げる。需要は Poisson、リードタイムも確率的（ $\{1..5\}$ 一様）。 $L$ がランダムだと $\sigma_d\sqrt{L}$ の安全在庫式（安全在庫と発注点）は前提が崩れ、しかも注文が追い越し合うのでコスト最小の $S$ は閉形式で出せません。 $S$ のグリッドを複製して信頼区間つきでシミュし、最適 $S^*$ を探します。冒頭で、2つの“もっともらしい閉形式”が食い違って真値を挟むことも見せます。

import numpy as np
import pandas as pd
from scipy.stats import norm

rng = np.random.default_rng(20260627)

# 周期発注の基在庫(base-stock)方策：毎期 在庫ポジションを S まで引き上げる。
# 需要は確率的（Poisson）・リードタイムも確率的（{1..5}一様）。
# L が確率変数だと σ_d*sqrt(L) の安全在庫式は前提（L 一定）が崩れ、
# しかも周期発注で注文が追い越し合うので、コスト最小の S は閉形式で出せない → 探索する。
mu_D = 10.0          # 1期あたり平均需要
L_choices = np.array([1, 2, 3, 4, 5])   # リードタイム（期）：一様, 平均3
h = 1.0              # 保管コスト（円/個/期, 期末の手持ち在庫に課金）
pen = 9.0            # 欠品ペナルティ（円/個/期, 期末のバックオーダーに課金）
# 望ましいサービス水準の目安は臨界比 pen/(pen+h) = 0.9 だが、S* は探索で決める。

Lmax = int(L_choices.max())

# 参考：2つの“もっともらしい閉形式”は食い違い、真値を挟む（だから探索が要る）
z = norm.ppf(pen / (pen + h))            # 臨界比 0.9 の正規分位点
EL = L_choices.mean()                    # 平均リードタイム = 3
varL = L_choices.var()                   # リードタイムの分散 = 2
# A: L 一定とみなす素朴式（保護区間 EL+1 期, 需要分散だけ）→ L のばらつきを無視
sigmaA = np.sqrt(mu_D * (EL + 1))
S_naiveA = mu_D * (EL + 1) + z * sigmaA
# B: L のばらつきも入れるが注文の追い越しを無視（保護区間需要の正規近似）
sigmaB = np.sqrt(mu_D * (EL + 1) + mu_D**2 * varL)
S_naiveB = mu_D * (EL + 1) + z * sigmaB
print(f"素朴式A（L一定とみなす）        S ~ {S_naiveA:.1f}（リードタイムのばらつきを無視→過少）")
print(f"素朴式B（L変動・追い越し無視）  S ~ {S_naiveB:.1f}（追い越しを無視→過大）")
print("どちらが正しいか閉形式では決まらない → シミュレーションで探索する。")
print()

def simulate(S, demand, leadtime, warmup):
    """基在庫水準 S を1本のサンプルパスで評価し、ウォームアップ後の平均コストを返す"""
    n = len(demand)
    receipts = np.zeros(n + Lmax + 2)    # 到着期ごとの入荷量（配列パイプライン）
    net_inv = float(S)                   # 正味在庫（手持ち−バックオーダー, 負なら欠品）
    on_order = 0.0                       # 発注残（輸送中）
    cost = 0.0
    for t in range(n):
        net_inv += receipts[t]; on_order -= receipts[t]   # 入荷
        IP = net_inv + on_order                           # 在庫ポジション
        order = S - IP                                    # S まで引き上げる
        if order > 1e-9:
            arr = t + int(leadtime[t])
            receipts[arr] += order; on_order += order
        net_inv -= demand[t]                              # 当期需要
        if t >= warmup:
            cost += h * max(net_inv, 0.0) + pen * max(-net_inv, 0.0)
    return cost / (n - warmup)

# --- S のグリッドを、独立な複製でシミュし、平均コストと95%CIを出す ---
S_grid = np.arange(34, 61, 2)
R = 30               # 複製数
warmup, horizon = 200, 2000
n_periods = warmup + horizon

means, halfwidths = [], []
for S in S_grid:
    rep_costs = np.empty(R)
    for r in range(R):
        demand = rng.poisson(mu_D, n_periods)
        leadtime = rng.choice(L_choices, n_periods)
        rep_costs[r] = simulate(S, demand, leadtime, warmup)
    m = rep_costs.mean()
    hw = 1.96 * rep_costs.std(ddof=1) / np.sqrt(R)   # CLT の95%信頼区間 半幅
    means.append(m); halfwidths.append(hw)
means = np.array(means); halfwidths = np.array(halfwidths)

S_star = S_grid[np.argmin(means)]
df = pd.DataFrame({"S": S_grid, "平均コスト": means, "CI半幅": halfwidths})
print(df.to_string(index=False, float_format=lambda x: f"{x:.3f}"))
print()
print(f"推定された最適基在庫 S* = {S_star}（平均コスト最小, 95%CI 半幅 {halfwidths[np.argmin(means)]:.3f}）")
print(f"  S* の平均コスト = {means.min():.3f} 円/期")
print("閉形式が無いので、複製の平均とCIで“ノイズ込みの推定値”として最小化している。")

出力：

素朴式A（L一定とみなす）        S ~ 48.1（リードタイムのばらつきを無視→過少）
素朴式B（L変動・追い越し無視）  S ~ 59.9（追い越しを無視→過大）
どちらが正しいか閉形式では決まらない → シミュレーションで探索する。

 S  平均コスト  CI半幅
34 73.519 1.097
36 61.813 0.943
38 52.421 0.771
40 44.969 0.765
42 36.711 0.625
44 31.425 0.515
46 26.868 0.658
48 23.864 0.427
50 21.774 0.318
52 20.344 0.284
54 20.103 0.266
56 20.289 0.195
58 20.806 0.182
60 21.907 0.162

推定された最適基在庫 S* = 54（平均コスト最小, 95%CI 半幅 0.266）
  S* の平均コスト = 20.103 円/期
閉形式が無いので、複製の平均とCIで“ノイズ込みの推定値”として最小化している。

出力の意味：まず2つの“もっともらしい閉形式”が食い違います。リードタイムを平均で固定して $\sigma\sqrt{L}$ を使う素朴式Aは $S\approx48$ （ $L$ のばらつきを無視するので過少）。 $L$ の分散まで入れるが注文の追い越しを無視する素朴式Bは $S\approx60$ （追い越しで実際は分散が下がるのに過大）。どちらが正しいか式では決まらない——だから探索します。シミュレーションのコストは $S$ について見事なU字で、低い $S$ では欠品ペナルティ（係数9）が、高い $S$ では保管費が効く。最小は $S^*=54$ （平均コスト20.103円/期）で、AとBのちょうど中間に落ち着きました。閉形式の片方を信じていたら6個ぶん過少（A）か過大（B）に在庫を持つところでした。なお $S=52,54,56$ のコスト（20.344／20.103／20.289）は接近しており、信頼区間（半幅0.2〜0.3）が重なります——「54が本当に最小か」は次節の CRN で詰めます。

5. CRN で比較を鋭く・ノイズの最小値を選ばない（コード）

前節で僅差だった $S=54$ と $S=56$ 。独立乱数と**共通乱数（CRN）**で比較し、CRN が差の分散をどれだけ縮めるかを見ます。あわせて、単一複製の生の最小がいかに当てにならないかも確かめます。

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

rng = np.random.default_rng(20260627)

mu_D = 10.0
L_choices = np.array([1, 2, 3, 4, 5])
h, pen = 1.0, 9.0
Lmax = int(L_choices.max())

def simulate(S, demand, leadtime, warmup):
    n = len(demand)
    receipts = np.zeros(n + Lmax + 2)
    net_inv = float(S); on_order = 0.0; cost = 0.0
    for t in range(n):
        net_inv += receipts[t]; on_order -= receipts[t]
        order = S - (net_inv + on_order)
        if order > 1e-9:
            arr = t + int(leadtime[t]); receipts[arr] += order; on_order += order
        net_inv -= demand[t]
        if t >= warmup:
            cost += h * max(net_inv, 0.0) + pen * max(-net_inv, 0.0)
    return cost / (n - warmup)

warmup, horizon = 200, 2000
n = warmup + horizon

# --- 隣接する2案 S_a, S_b を「独立乱数」と「共通乱数(CRN)」で比較する ---
S_a, S_b = 54, 56          # コード1で僅差だった2案（どちらが本当に良い？）
R = 300

# 独立：各案に別々の乱数列を割り当てる → 差はノイズが2つ分のる
diff_ind = np.empty(R)
for r in range(R):
    ca = simulate(S_a, rng.poisson(mu_D, n), rng.choice(L_choices, n), warmup)
    cb = simulate(S_b, rng.poisson(mu_D, n), rng.choice(L_choices, n), warmup)
    diff_ind[r] = cb - ca

# CRN：両案に同じ乱数列を使う → 共通ノイズが相殺し、差の分散が激減
diff_crn = np.empty(R)
ca_crn = np.empty(R); cb_crn = np.empty(R)
for r in range(R):
    demand = rng.poisson(mu_D, n)               # 同じ需要列
    leadtime = rng.choice(L_choices, n)         # 同じリードタイム列
    ca_crn[r] = simulate(S_a, demand, leadtime, warmup)
    cb_crn[r] = simulate(S_b, demand, leadtime, warmup)
    diff_crn[r] = cb_crn[r] - ca_crn[r]

se_ind = diff_ind.std(ddof=1) / np.sqrt(R)       # 平均差の標準誤差
se_crn = diff_crn.std(ddof=1) / np.sqrt(R)
print(f"2案の比較： コスト(S={S_b}) - コスト(S={S_a})  （正なら S={S_a} が安い）")
print(f"  独立乱数 : 平均差 {diff_ind.mean():+.4f}  95%CI [{diff_ind.mean()-1.96*se_ind:+.4f}, {diff_ind.mean()+1.96*se_ind:+.4f}]  差のSD {diff_ind.std(ddof=1):.4f}")
print(f"  共通乱数 : 平均差 {diff_crn.mean():+.4f}  95%CI [{diff_crn.mean()-1.96*se_crn:+.4f}, {diff_crn.mean()+1.96*se_crn:+.4f}]  差のSD {diff_crn.std(ddof=1):.4f}")
print(f"  => 独立のCIは0をまたぎ「差なし」と区別できない。CRNは0を含まず S={S_a} の勝ちと結論できる")
print(f"  差のSD が {diff_ind.std(ddof=1) / diff_crn.std(ddof=1):.1f} 分の1に縮小・corr(cost_a,cost_b)={np.corrcoef(ca_crn, cb_crn)[0,1]:.4f}（高いほどCRN有効）")
print()

# --- ノイズの“生の最小値”を選ぶ危険：単一複製の最小は当てにならない ---
S_grid = np.arange(46, 63, 2)
R_plot = 80
M = np.empty((R_plot, len(S_grid)))
for r in range(R_plot):
    demand = rng.poisson(mu_D, n)               # CRN：1複製の中では全 S で共通列
    leadtime = rng.choice(L_choices, n)
    for j, S in enumerate(S_grid):
        M[r, j] = simulate(S, demand, leadtime, warmup)
mean_curve = M.mean(axis=0)
hw = 1.96 * M.std(axis=0, ddof=1) / np.sqrt(R_plot)
S_star = S_grid[np.argmin(mean_curve)]
single_argmins = S_grid[np.argmin(M, axis=1)]    # 各複製を1本だけで選んだ最小
frac_correct = np.mean(single_argmins == S_star)
example = int(np.argmax(single_argmins != S_star))   # 最小がS*とズレた一例
print(f"複製平均({R_plot}本)の最小 S* = {S_star}（信頼できる）")
print(f"単一複製1本の生の最小は S={single_argmins.min()}..{single_argmins.max()} に散らばり、S*={S_star} を当てるのは {frac_correct:.0%} だけ")
print(f"  例：複製#{example} の生の最小 = {single_argmins[example]}（ノイズで誤誘導）")

# --- 図：コスト vs S（CRN平均＋95%CI帯）と、単一複製のギザギザ ---
plt.figure(figsize=(10, 5.5))
plt.fill_between(S_grid, mean_curve - hw, mean_curve + hw,
                 color="#1f77b4", alpha=0.25, label="95%CI 帯（複製平均）")
plt.plot(S_grid, mean_curve, "o-", color="#1f77b4", label=f"平均コスト（CRN, {R_plot}複製）")
plt.plot(S_grid, M[example], "x--", color="#ff7f0e", alpha=0.8, label=f"単一複製#{example}（生のノイズ）")
plt.scatter([S_star], [mean_curve[np.argmin(mean_curve)]], color="#2ca02c",
            zorder=6, s=140, marker="*", label=f"信頼できる最小 S*={S_star}")
plt.scatter([single_argmins[example]], [M[example].min()], color="#d62728",
            zorder=6, s=70, label=f"単一複製の生の最小={single_argmins[example]}（罠）")
plt.xlabel("基在庫水準 S"); plt.ylabel("平均コスト（円/期）")
plt.title("シミュレーション最適化：複製とCIで“ノイズの最小値”を選ばない（CRNで比較を鋭く）")
plt.legend(); plt.tight_layout(); plt.show()

出力：

2案の比較： コスト(S=56) - コスト(S=54)  （正なら S=54 が安い）
  独立乱数 : 平均差 +0.0842  95%CI [-0.0047, +0.1732]  差のSD 0.7861
  共通乱数 : 平均差 +0.1404  95%CI [+0.1195, +0.1612]  差のSD 0.1842
  => 独立のCIは0をまたぎ「差なし」と区別できない。CRNは0を含まず S=54 の勝ちと結論できる
  差のSD が 4.3 分の1に縮小・corr(cost_a,cost_b)=0.9644（高いほどCRN有効）

複製平均(80本)の最小 S* = 54（信頼できる）
単一複製1本の生の最小は S=52..56 に散らばり、S*=54 を当てるのは 79% だけ
  例：複製#2 の生の最小 = 56（ノイズで誤誘導）

出力の意味：肝は比較の鋭さです。独立乱数だと「 $S=56$ は $S=54$ より平均 +0.084 円高い」が、その95%信頼区間は $[-0.005,+0.173]$ で 0 をまたぐ——統計的には「差があるとは言えない」。これでは54と56のどちらを採るか決められません。一方CRNでは平均差 +0.140 円、95%信頼区間 $[+0.120,+0.161]$ は 0 を含まず、 $S=54$ の勝ちと自信を持って結論できます。同じ複製数なのに、差の標準偏差が 0.786→0.184 と約4.3分の1に縮んだから。2案のコストの相関が 0.96 と高い（同じ需要・リードタイムに同じように反応する）ので、共通変動が差し引きで消え、方策の違いだけが残ったわけです。後半は「ノイズの生の最小を選ぶ罠」。単一複製1本だけで最小を選ぶと、答えは $S=52\sim56$ に散らばり、本当の $S^*=54$ を当てるのは79%だけ（複製#2 では誤って56を指す）。図でも、CRN平均の青い曲線は細い信頼帯で滑らかに $S^*=54$ （緑星）で底を打つのに、単一複製のオレンジはギザギザで最小が56（赤）にずれています。複製・信頼区間・CRN——この3点セットがシミュレーション最適化の作法です。

⚠️ よくある誤解

「シミュレーションの出力はそのまま答え」ではない：出力はノイズ込みの推定値です。複製して平均・信頼区間を出し、精度を示すのが必須。1回の結果や、点推定だけの数字を「答え」として比較してはいけません。 $\sqrt m$ で効くので、精度がほしければ複製を増やします。
「たくさん試して一番良かった案を選べばよい」ではない：多数案を1回ずつ評価して最小（最大）を選ぶと、運悪く下振れ（上振れ）した案を拾う**最適化バイアス（winner’s curse）**が起きます。選んだ案の真の性能は見かけより悪い。複製と信頼区間で差を確かめてから選ぶ、僅差なら追加複製や CRN で詰める、が鉄則です。
「比較も独立に回せばよい」ではない：2案の優劣を知りたいなら共通乱数（CRN）で条件を揃えて比較します。独立に回すと差にノイズが2案分のり、本当はある差を見逃します（本稿で CI が 0 をまたいだ通り）。CRN は同じ計算量で比較を鋭くする、ほぼ無料の分散減です（ただし両案が乱数に同じ向きで反応するときに効く。逆相関だと逆効果）。
「ウォームアップや定常性は気にしなくていい」ではない：初期状態の影響（過渡期）をウォームアップ期間として捨てる、十分長い地平で定常を測る、入力分布が現実に合うか検証（validation）する——シミュレーションはモデルが正しくて初めて意味を持ちます。「動いた数字」を鵜呑みにせず、既知の閉形式ケース（本テキストでは constant $L$ の新聞売り子）で校正してから使います。
「閉形式があってもシミュで十分」ではない：閉形式があるならそれを使うべきです。式は厳密・即時・洞察（どの変数がどう効くか）を与え、シミュは近似・遅い・ノイズ込み。シミュレーションは閉形式が無い／前提が破れたときの最終手段。本稿も、 $L$ 一定なら $\sigma\sqrt L$ で済む（安全在庫と発注点）ところを、 $L$ がランダムで追い越しが絡むから探索に頼りました。手法の選択も最適化です（要最新確認）。