多品種組立ラインの投資・運用条件をどう決めるか

シミュレーション最適化で「平均が良い案」から「現場で使える案」へ（No.081〜No.090）

需要変動、設備停止、段取り替えが重なる多品種組立ラインを題材に、人員・中間バッファ・予防保全間隔をどう決めるかを扱います。単一の予測値だけでなく、利益、納期遵守率、仕掛品（WIP）、停止リスクを同時に見ながら、探索手法と意思決定プロセスを段階的に確認します。

[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。

1. はじめに：この記事で扱う製造業の実務課題

生産能力を上げる方法は一つではありません。増員、バッファ拡張、保全頻度の見直しは、それぞれ能力を改善する一方で、人件費、在庫費、保全費を増やします。しかも需要と故障は確率的に発生するため、平均需要だけで決めた案が繁忙日や故障日に崩れることがあります。

本稿の意思決定は、次の制約付き目的として整理します。

$$\max_x\; \mathbb{E}[\Pi(x,\xi)] \quad \text{s.t.}\quad P\{S(x,\xi) \geq 0.95\} \text{を重視},\; x \in \mathcal{X}$$

$x$ は人員数・バッファ容量・保全間隔、$\xi$ は需要・故障などの不確実性、$\Pi$ は日次利益、$S$ は納期遵守率です。期待利益だけを最大化せず、サービス水準と下振れも併記します。

2. 現場でよくある状況

• 各部門が別々のKPIで案を評価し、全体最適の比較表がない
• 過去平均で作った能力計画が、需要上振れや設備停止に耐えない
• 改善案の候補が多く、全組合せの実機試験はできない
• シミュレーションは作ったが、最適化・実績更新・稟議へつながっていない

3. なぜこの問題は判断が難しいのか

シミュレーションの出力には乱数誤差があり、目的関数は滑らかとは限りません。さらに、平均利益が高い案と納期遵守率が安定する案は一致しないことがあります。探索アルゴリズムの結果をそのまま採用するのではなく、制約、再現性、複数シナリオでの頑健性を確認する必要があります。

4. 今回扱うノックの全体像

No.	テーマ	意思決定上の役割
081	シミュレーション最適化	目的・変数・制約・乱数を定義する
082	グリッドサーチ	説明しやすい候補比較を作る
083	遺伝的アルゴリズム	広い離散空間を集団探索する
084	焼きなまし法	局所解を抜ける近傍探索を行う
085	ベイズ最適化	少ない試行で有望条件を探す
086	Optuna入門	探索履歴と再現可能な実装を整える
087	シミュレーション×数理最適化	シナリオ評価と資源配分を分担する
088	強化学習との統合	状態に応じた運用ルールを学ぶ
089	デジタルツイン	実績でモデルを更新し続ける
090	What-if分析	経営シナリオ別の頑健性を比べる

5. Python環境の準備

乱数生成器を引数で渡し、同じシードなら同じ結果を再現できるようにします。グラフはすべて matplotlib で作成します。

import logging
import warnings

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
from scipy.optimize import milp, LinearConstraint, Bounds
from scipy.stats import norm
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, WhiteKernel

warnings.filterwarnings("ignore")
pd.set_option("display.max_columns", 20)
BASE_SEED = 1609
print("NumPy:", np.__version__)
print("pandas:", pd.__version__)

NumPy: 2.5.1
pandas: 3.0.3

6. 架空データの作成

1日を1回の試行とし、需要は正規分布、突発停止時間は指数分布を基礎に生成します。人員 workers は処理能力、buffer は流動の詰まり、pm_interval は故障リスクと保全費に影響する簡略モデルです。これは設備の物理挙動を忠実に再現するモデルではなく、意思決定プロセスを学ぶための代理モデルです。

def simulate_policy(workers, buffer, pm_interval, n_days=240, seed=BASE_SEED,
                    demand_mean=112, failure_scale=1.0):
    rng = np.random.default_rng(seed)
    demand = np.maximum(60, rng.normal(demand_mean, 18, n_days)).round()
    failure_prob = np.clip(0.025 + 0.00045 * (pm_interval - 20) ** 1.35, 0.02, 0.24)
    failure = rng.random(n_days) < failure_prob
    downtime = failure * rng.exponential(2.7 * failure_scale, n_days)
    variation = rng.normal(0, 4, n_days)
    nominal_capacity = 79 + 8.3 * workers + 1.7 * np.sqrt(buffer)
    capacity = np.maximum(0, nominal_capacity + variation - 5.8 * downtime)
    shipped = np.minimum(demand, capacity)
    shortage = demand - shipped
    wip = np.minimum(buffer, np.maximum(0, demand - 0.84 * capacity))
    revenue = 1.65 * shipped
    cost = (0.40 * shipped + 11.5 * workers + 0.075 * buffer
            + 84 / pm_interval + 1.8 * downtime + 2.9 * shortage + 0.055 * wip)
    profit = revenue - cost
    return pd.DataFrame({
        "demand": demand, "shipped": shipped, "shortage": shortage,
        "wip": wip, "downtime_h": downtime, "profit_10k_yen": profit,
        "on_time": shortage <= 2,
    })

def evaluate_policy(workers, buffer, pm_interval, n_days=240, seed=BASE_SEED,
                    demand_mean=112, failure_scale=1.0):
    d = simulate_policy(workers, buffer, pm_interval, n_days, seed,
                        demand_mean, failure_scale)
    return {
        "workers": workers, "buffer": buffer, "pm_interval": pm_interval,
        "mean_profit": d["profit_10k_yen"].mean(),
        "p05_profit": d["profit_10k_yen"].quantile(0.05),
        "on_time_rate": d["on_time"].mean(),
        "mean_wip": d["wip"].mean(),
        "downtime_h": d["downtime_h"].mean(),
    }

baseline_daily = simulate_policy(4, 30, 45)
baseline = pd.DataFrame([evaluate_policy(4, 30, 45)]).round(2)
baseline

	workers	buffer	pm_interval	mean_profit	p05_profit	on_time_rate	mean_wip	downtime_h
0	4	30	45	69.48	24.19	0.67	13.61	0.17

fig, axes = plt.subplots(1, 2, figsize=(11, 3.8))
axes[0].plot(baseline_daily.index[:60], baseline_daily["demand"][:60], label="需要")
axes[0].plot(baseline_daily.index[:60], baseline_daily["shipped"][:60], label="出荷")
axes[0].set_title("基準案：最初の60日の日次推移")
axes[0].set_xlabel("日")
axes[0].set_ylabel("数量（個/日）")
axes[0].grid(True, alpha=.3)
axes[0].legend()
axes[1].hist(baseline_daily["profit_10k_yen"], bins=20, edgecolor="white")
axes[1].axvline(baseline_daily["profit_10k_yen"].quantile(.05), color="red", linestyle="--", label="5%点")
axes[1].set_title("基準案：日次利益の分布")
axes[1].set_xlabel("利益（万円/日）")
axes[1].set_ylabel("日数")
axes[1].grid(True, alpha=.3)
axes[1].legend()
plt.tight_layout()
plt.show()

No.081：シミュレーション最適化とは

実務での意味

シミュレーション最適化は、現実には試せない多数の設備・要員条件を仮想工場で評価し、望ましい案を探す枠組みです。重要なのは「最適化ソフトを使うこと」ではなく、意思決定変数、評価KPI、守るべき制約を合意することです。

分析・モデル化の考え方

ここでは納期遵守率95%未満に大きなペナルティを置くスコアを使います。

$J(x)=E[\Pi]-300\max(0,0.95-E[S])-0.15E[WIP]$

係数300は経営判断を反映する重みであり、統計的に自動決定される値ではありません。

Pythonで確認する

def business_score(row):
    return row["mean_profit"] - 300 * max(0, .95 - row["on_time_rate"]) - .15 * row["mean_wip"]

candidates_081 = [(4, 30, 45), (5, 30, 45), (4, 50, 30), (5, 50, 30)]
result_081 = pd.DataFrame([evaluate_policy(*x, seed=8100+i) for i, x in enumerate(candidates_081)])
result_081["score"] = result_081.apply(business_score, axis=1)
result_081.sort_values("score", ascending=False).round(2)

	workers	buffer	pm_interval	mean_profit	p05_profit	on_time_rate	mean_wip	downtime_h	score
3	5	50	30	69.57	39.39	0.87	8.33	0.02	43.32
1	5	30	45	68.61	31.82	0.80	10.00	0.21	22.11
2	4	50	30	75.09	37.45	0.74	12.95	0.04	9.40
0	4	30	45	74.24	28.65	0.71	13.29	0.10	-0.26

結果の読み取り

候補ごとに平均利益だけでなく、納期遵守率と利益5%点が異なります。スコア首位は「何を重視したか」に依存するため、採用前にペナルティ係数とサービス水準を生産・営業・経理で合意します。

No.082：グリッドサーチ

実務での意味

選択肢が少ない設備仕様や勤務編成では、候補を網羅するグリッドサーチが監査・稟議で説明しやすい方法です。

分析・モデル化の考え方

人員4〜6名、バッファ20〜60個、保全間隔20〜60日を離散候補として総当たりします。同一候補を複数乱数系列で評価するのが理想ですが、ここでは比較を簡潔にするため各候補240日で評価します。

Pythonで確認する

grid_rows = []
for w in range(4, 7):
    for b in range(20, 61, 10):
        for pm in range(20, 61, 10):
            row = evaluate_policy(w, b, pm, seed=8200 + w*100 + b + pm)
            row["score"] = business_score(row)
            grid_rows.append(row)
grid = pd.DataFrame(grid_rows).sort_values("score", ascending=False).reset_index(drop=True)
grid.head(8).round(2)

	workers	buffer	pm_interval	mean_profit	p05_profit	on_time_rate	mean_wip	downtime_h	score
0	6	50	30	63.87	29.45	0.95	4.83	0.05	63.14
1	6	40	40	64.77	30.90	0.95	5.39	0.11	62.71
2	6	60	60	65.68	33.85	0.94	5.29	0.11	62.38
3	6	60	20	61.98	27.30	0.97	4.43	0.04	61.31
4	6	60	40	64.74	30.46	0.94	5.51	0.13	60.17
5	6	30	50	65.14	32.07	0.93	5.98	0.17	59.24
6	6	60	30	59.40	23.64	0.95	4.87	0.17	58.67
7	6	50	20	59.75	24.73	0.95	4.78	0.02	57.78

pivot = grid.groupby(["workers", "buffer"], as_index=False)["score"].max().pivot(index="workers", columns="buffer", values="score")
fig, ax = plt.subplots(figsize=(7, 3.8))
im = ax.imshow(pivot.values, aspect="auto", cmap="viridis")
ax.set_xticks(range(len(pivot.columns)), pivot.columns)
ax.set_yticks(range(len(pivot.index)), pivot.index)
ax.set_title("グリッドサーチ：人員・バッファ別の最高スコア")
ax.set_xlabel("バッファ容量（個）")
ax.set_ylabel("人員（名）")
ax.grid(False)
fig.colorbar(im, ax=ax, label="スコア")
plt.tight_layout()
plt.show()

結果の読み取り

上位表から、人員増が常に得ではなく、バッファと保全間隔との組合せで採算が変わることが分かります。候補数が数千〜数百万になる場合は、次の確率的探索へ切り替えます。

No.083：遺伝的アルゴリズム

実務での意味

レイアウト、設備種、人員など離散変数の組合せが広い場合、遺伝的アルゴリズム（GA）は複数候補を並行して改善できます。

分析・モデル化の考え方

候補を遺伝子 [workers, buffer, pm_interval] とみなし、選択・交叉・突然変異を繰り返します。最適性の証明はできないため、複数シードと最終候補の再評価が必要です。

Pythonで確認する

def clip_gene(g):
    return np.array([np.clip(round(g[0]), 3, 7), np.clip(round(g[1]/5)*5, 10, 80),
                     np.clip(round(g[2]/5)*5, 15, 75)], dtype=int)

rng = np.random.default_rng(8300)
pop = np.column_stack([rng.integers(3, 8, 18), rng.integers(2, 17, 18)*5, rng.integers(3, 16, 18)*5])
ga_history = []
for gen in range(16):
    scored = []
    for i, g in enumerate(pop):
        r = evaluate_policy(*g, n_days=160, seed=83000 + gen*100 + i)
        scored.append((business_score(r), g.copy()))
    scored.sort(key=lambda z: z[0], reverse=True)
    ga_history.append(scored[0][0])
    elite = [g for _, g in scored[:6]]
    children = elite.copy()
    while len(children) < len(pop):
        a, b = rng.choice(len(elite), 2, replace=False)
        child = np.where(rng.random(3) < .5, elite[a], elite[b]).astype(float)
        if rng.random() < .45:
            child += rng.choice([[1,0,0],[-1,0,0],[0,5,0],[0,-5,0],[0,0,5],[0,0,-5]])
        children.append(clip_gene(child))
    pop = np.array(children)
ga_best = scored[0][1]
pd.DataFrame([evaluate_policy(*ga_best, n_days=1000, seed=83999)]).round(2)

	workers	buffer	pm_interval	mean_profit	p05_profit	on_time_rate	mean_wip	downtime_h
0	6	70	40	60.61	22.4	0.96	3.57	0.11

plt.figure(figsize=(7, 3.5))
plt.plot(range(1, len(ga_history)+1), ga_history, marker="o")
plt.title("遺伝的アルゴリズムの探索推移")
plt.xlabel("世代")
plt.ylabel("世代内の最高スコア")
plt.grid(True, alpha=.3)
plt.tight_layout()
plt.show()

結果の読み取り

世代を追って有望領域へ移動しますが、途中の上下はシミュレーションノイズの影響も含みます。最終案は1,000日で再評価し、探索中の偶然の高得点を採用しないようにします。

No.084：焼きなまし法

実務での意味

現行案から一度に少しだけ条件を変える改善活動では、焼きなまし法（SA）の「近傍案を試す」考え方が自然です。

分析・モデル化の考え方

改善案だけでなく、温度 $T$ に応じて悪化案も確率 $\exp(\{\Delta J/T\})$ で受け入れ、局所解からの脱出を狙います。温度の下げ方と近傍定義が結果を左右します。

Pythonで確認する

rng = np.random.default_rng(8400)
current = np.array([4, 30, 45])
current_score = business_score(evaluate_policy(*current, n_days=180, seed=84000))
best, best_score = current.copy(), current_score
sa_history = []
moves = np.array([[1,0,0],[-1,0,0],[0,5,0],[0,-5,0],[0,0,5],[0,0,-5]])
for step in range(120):
    temp = 18 * (0.96 ** step) + .2
    proposal = clip_gene(current + moves[rng.integers(len(moves))])
    score = business_score(evaluate_policy(*proposal, n_days=180, seed=84100 + step))
    if score >= current_score or rng.random() < np.exp((score-current_score)/temp):
        current, current_score = proposal, score
    if current_score > best_score:
        best, best_score = current.copy(), current_score
    sa_history.append(best_score)
pd.DataFrame([evaluate_policy(*best, n_days=1000, seed=84999)]).round(2)

	workers	buffer	pm_interval	mean_profit	p05_profit	on_time_rate	mean_wip	downtime_h
0	6	60	40	61.24	23.15	0.94	4.84	0.17

plt.figure(figsize=(7, 3.5))
plt.plot(sa_history)
plt.title("焼きなまし法：探索中の最高スコア")
plt.xlabel("反復回数")
plt.ylabel("最高スコア")
plt.grid(True, alpha=.3)
plt.tight_layout()
plt.show()

結果の読み取り

現行条件から始めても、悪化案を一時的に許容することで別の有望領域へ移れます。ただし単一経路なので、開始点を変えた複数回実行が実務上の基本です。

No.085：ベイズ最適化

実務での意味

1回の高精度シミュレーションに数時間かかる場合、少ない評価回数で有望条件を見つけるベイズ最適化が有効です。

分析・モデル化の考え方

評価済み点からガウス過程でスコアの平均 $\mu(x)$ と不確実性 $\sigma(x)$ を推定し、Expected Improvement（EI）が大きい点を次に評価します。既知の良好領域の活用と未知領域の探索を両立します。

Pythonで確認する

rng = np.random.default_rng(8500)
pool = np.array([(w,b,p) for w in range(3,8) for b in range(10,81,5) for p in range(15,76,5)])
chosen = list(rng.choice(len(pool), 8, replace=False))
y = [business_score(evaluate_policy(*pool[i], n_days=180, seed=85000+k)) for k,i in enumerate(chosen)]
bo_best = [max(y)]
for t in range(12):
    X = pool[chosen] / np.array([7, 80, 75])
    gp = GaussianProcessRegressor(kernel=Matern(nu=2.5)+WhiteKernel(.5), normalize_y=True,
                                  random_state=BASE_SEED).fit(X, y)
    available = np.array([i for i in range(len(pool)) if i not in chosen])
    mu, sd = gp.predict(pool[available] / np.array([7,80,75]), return_std=True)
    z = (mu - max(y)) / np.maximum(sd, 1e-9)
    ei = (mu-max(y))*norm.cdf(z) + sd*norm.pdf(z)
    idx = available[np.argmax(ei)]
    chosen.append(idx)
    y.append(business_score(evaluate_policy(*pool[idx], n_days=180, seed=85100+t)))
    bo_best.append(max(y))
bo_gene = pool[chosen[int(np.argmax(y))]]
pd.DataFrame([evaluate_policy(*bo_gene, n_days=1000, seed=85999)]).round(2)

	workers	buffer	pm_interval	mean_profit	p05_profit	on_time_rate	mean_wip	downtime_h
0	6	30	45	62.94	26.88	0.91	5.52	0.1

plt.figure(figsize=(7, 3.5))
plt.plot(range(len(bo_best)), bo_best, marker="o")
plt.title("ベイズ最適化：追加評価による最高スコア")
plt.xlabel("追加評価回数")
plt.ylabel("最高スコア")
plt.grid(True, alpha=.3)
plt.tight_layout()
plt.show()

結果の読み取り

全候補を評価せず、代理モデルが選んだ点だけを追加評価しています。次点候補と予測不確実性も残せるため、追加検証の優先順位付けに使えます。

No.086：Optuna入門

実務での意味

Optunaを使うと、試行番号、パラメータ、目的値を一貫して保存でき、担当者の手作業による探索漏れを減らせます。

分析・モデル化の考え方

目的関数内で候補を提案し、シミュレーション結果を返します。本番ではRDBへの永続化、複数ワーカー、試行の停止条件を設計します。

Pythonで確認する

import optuna
optuna.logging.set_verbosity(optuna.logging.WARNING)

def objective(trial):
    w = trial.suggest_int("workers", 3, 7)
    b = trial.suggest_int("buffer", 10, 80, step=5)
    pm = trial.suggest_int("pm_interval", 15, 75, step=5)
    r = evaluate_policy(w, b, pm, n_days=180, seed=86000 + trial.number)
    return business_score(r)

study = optuna.create_study(direction="maximize", sampler=optuna.samplers.TPESampler(seed=8600))
study.optimize(objective, n_trials=35)
optuna_best = [study.best_params[k] for k in ["workers", "buffer", "pm_interval"]]
pd.DataFrame([{**study.best_params, "best_score": study.best_value}]).round(2)

	workers	buffer	pm_interval	best_score
0	5	75	55	66.18

trials = study.trials_dataframe()
plt.figure(figsize=(7, 3.5))
plt.scatter(trials["number"], trials["value"], alpha=.65, label="各試行")
plt.plot(trials["number"], trials["value"].cummax(), color="red", label="累積最高")
plt.title("Optunaの探索履歴")
plt.xlabel("試行番号")
plt.ylabel("スコア")
plt.grid(True, alpha=.3)
plt.legend()
plt.tight_layout()
plt.show()

結果の読み取り

試行履歴から改善の頭打ちや外れ値を確認できます。探索結果だけでなく、コード、乱数シード、入力版数をセットで記録することが再現性につながります。

No.087：シミュレーション×数理最適化

実務での意味

シミュレーションは不確実な能力評価に強く、数理最適化は明示的な資源制約の下で配分を決めることに強みがあります。両者を分担させると、複数ラインへの改善予算配分まで扱えます。

分析・モデル化の考え方

各ライン・改善案の効果をシミュレーションで見積もり、その期待増益を係数として0-1整数計画へ渡します。

$\max \sum_i v_i z_i\quad \text{s.t.}\quad \sum_i c_i z_i\leq B,\; z_i\in\{0,1\}$

Pythonで確認する

projects = pd.DataFrame({
    "project": ["Aライン増員", "Aラインバッファ", "Bライン予防保全", "Bライン増員", "検査自動化"],
    "cost": [55, 28, 36, 50, 68],
    "annual_value": [74, 31, 49, 61, 82],
})
budget = 120
res = milp(c=-projects["annual_value"].to_numpy(), integrality=np.ones(len(projects)),
           bounds=Bounds(np.zeros(len(projects)), np.ones(len(projects))),
           constraints=LinearConstraint(projects["cost"].to_numpy()[None, :], -np.inf, budget))
projects["selected"] = np.rint(res.x).astype(int)
projects

	project	cost	annual_value	selected
0	Aライン増員	55	74	1
1	Aラインバッファ	28	31	1
2	Bライン予防保全	36	49	1
3	Bライン増員	50	61	0
4	検査自動化	68	82	0

selected = projects.query("selected == 1")
print("使用予算:", selected["cost"].sum(), "/", budget)
print("期待年間価値:", selected["annual_value"].sum())

使用予算: 119 / 120
期待年間価値: 154

結果の読み取り

単純な費用対効果順ではなく、予算内で価値合計が最大になる組合せが選ばれます。シミュレーション係数自体に誤差があるため、採択境界の案件は係数を上下させた感度分析が必要です。

No.088：強化学習との統合

実務での意味

固定した最適条件ではなく、「滞留が多い日は応援を呼ぶ」といった状態依存ルールが必要な場合、強化学習（RL）が候補になります。

分析・モデル化の考え方

簡略化した日次環境で、状態を低・中・高の滞留、行動を応援なし・1名・2名、報酬を処理価値から人件費と滞留費を引いた値とします。Q学習で長期報酬を更新します。

Pythonで確認する

rng = np.random.default_rng(8800)
Q = np.zeros((3, 3))
alpha, gamma, epsilon = .12, .92, .18
for episode in range(600):
    backlog = int(rng.integers(0, 25))
    for day in range(50):
        state = min(backlog // 18, 2)
        action = rng.integers(3) if rng.random() < epsilon else int(np.argmax(Q[state]))
        demand = max(70, rng.normal(112, 18))
        capacity = 104 + 10 * action + rng.normal(0, 5)
        processed = min(backlog + demand, capacity)
        backlog = max(0, int(backlog + demand - processed))
        reward = 1.5 * processed - 13 * action - .35 * backlog
        next_state = min(backlog // 18, 2)
        Q[state, action] += alpha * (reward + gamma * Q[next_state].max() - Q[state, action])
policy = pd.DataFrame(Q.round(1), index=["滞留:低", "滞留:中", "滞留:高"],
                      columns=["応援0名", "応援1名", "応援2名"])
policy["選択行動"] = [f"応援{a}名" for a in np.argmax(Q, axis=1)]
policy

	応援0名	応援1名	応援2名	選択行動
滞留:低	1842.6	1825.0	1828.1	応援0名
滞留:中	1828.7	1827.5	1840.8	応援2名
滞留:高	1789.2	1812.2	1839.6	応援2名

fig, ax = plt.subplots(figsize=(6.5, 3.5))
im = ax.imshow(Q, cmap="Blues", aspect="auto")
ax.set_xticks(range(3), ["0名", "1名", "2名"])
ax.set_yticks(range(3), ["低", "中", "高"])
ax.set_title("学習後の状態・行動価値（Q値）")
ax.set_xlabel("応援人数")
ax.set_ylabel("滞留状態")
ax.grid(False)
fig.colorbar(im, ax=ax, label="Q値")
plt.tight_layout()
plt.show()

結果の読み取り

滞留状態ごとに選ばれる応援人数が変わり、固定配置とは異なる運用ルールが得られます。ただし実設備へ直接接続せず、オフライン評価、行動上限、手動介入、異常時の既定動作を先に定めます。

No.089：デジタルツイン

実務での意味

デジタルツインは3D表示そのものではなく、現場実績で状態・パラメータを更新し、将来を試算して意思決定へ返す循環です。

分析・モデル化の考え方

故障確率を固定せず、日々の停止有無でベータ分布のパラメータを更新します。事前分布 $Beta(a,b)$ に故障 $f$、非故障 $n$ を観測すると、事後分布は $Beta(a+f,b+n)$ です。

Pythonで確認する

rng = np.random.default_rng(8900)
true_prob = np.r_[np.repeat(.045, 60), np.repeat(.11, 60)]
observed = rng.random(120) < true_prob
a, b = 2., 38.
posterior_mean = []
for failed in observed:
    a += failed
    b += 1 - failed
    posterior_mean.append(a / (a + b))
twin = pd.DataFrame({"day": np.arange(1,121), "actual_failure_prob": true_prob,
                     "posterior_mean": posterior_mean, "failed": observed})
twin.tail().round(3)

	day	actual_failure_prob	posterior_mean	failed
115	116	0.11	0.064	False
116	117	0.11	0.064	False
117	118	0.11	0.070	True
118	119	0.11	0.069	False
119	120	0.11	0.069	False

plt.figure(figsize=(8, 3.6))
plt.plot(twin["day"], twin["actual_failure_prob"], linestyle="--", label="真の故障確率（架空）")
plt.plot(twin["day"], twin["posterior_mean"], label="ツインの推定値")
plt.scatter(twin.loc[twin.failed, "day"], np.repeat(.145, twin.failed.sum()), marker="x", color="red", label="故障観測")
plt.title("実績データによる故障確率の逐次更新")
plt.xlabel("日")
plt.ylabel("故障確率")
plt.grid(True, alpha=.3)
plt.legend()
plt.tight_layout()
plt.show()

結果の読み取り

60日目以降の劣化を観測に応じて追随します。追随が遅すぎる場合は忘却係数や状態空間モデルを検討します。データ遅延、センサー欠測、設備IDの不一致を監視しなければ、更新される「もっともらしい誤モデル」になります。

No.090：What-if分析

実務での意味

経営会議では単一予測より、「需要が15%増えたら」「故障復旧が長期化したら」という条件別の比較が求められます。What-if分析は、前提を明示して代替案の耐性を確認する方法です。

分析・モデル化の考え方

基準案、グリッド上位案、Optuna案を、通常・需要増・故障長期化・複合ストレスの同じシナリオで評価します。共通乱数により案同士の差を比較しやすくします。

Pythonで確認する

grid_best = grid.loc[0, ["workers", "buffer", "pm_interval"]].astype(int).tolist()
plans = {"基準案": [4,30,45], "グリッド上位案": grid_best, "Optuna案": optuna_best}
scenarios = {"通常": (112,1.0), "需要15%増": (129,1.0), "復旧長期化": (112,1.6), "複合ストレス": (129,1.6)}
rows = []
for s, (dm, fs) in scenarios.items():
    for name, x in plans.items():
        r = evaluate_policy(*x, n_days=1200, seed=9000, demand_mean=dm, failure_scale=fs)
        rows.append({"scenario": s, "plan": name, **r})
whatif = pd.DataFrame(rows)
whatif[["scenario", "plan", "mean_profit", "p05_profit", "on_time_rate", "mean_wip"]].round(2)

	scenario	plan	mean_profit	p05_profit	on_time_rate	mean_wip
0	通常	基準案	71.60	22.69	0.72	12.69
1	通常	グリッド上位案	61.70	25.70	0.94	4.93
2	通常	Optuna案	68.98	32.85	0.88	7.39
3	需要15%増	基準案	56.93	-26.61	0.35	22.05
4	需要15%増	グリッド上位案	71.30	29.88	0.77	14.43
5	需要15%増	Optuna案	72.47	9.93	0.66	18.71
6	復旧長期化	基準案	69.30	17.24	0.72	12.79
7	復旧長期化	グリッド上位案	60.56	24.45	0.94	5.10
8	復旧長期化	Optuna案	66.44	30.29	0.88	7.91
9	複合ストレス	基準案	54.44	-32.06	0.35	22.10
10	複合ストレス	グリッド上位案	70.02	29.45	0.76	14.55
11	複合ストレス	Optuna案	69.48	1.11	0.66	19.23

fig, axes = plt.subplots(1, 2, figsize=(11, 3.8))
for name in plans:
    d = whatif[whatif.plan == name]
    axes[0].plot(d.scenario, d.mean_profit, marker="o", label=name)
    axes[1].plot(d.scenario, d.on_time_rate, marker="o", label=name)
axes[0].set_title("シナリオ別の平均利益")
axes[0].set_xlabel("シナリオ")
axes[0].set_ylabel("利益（万円/日）")
axes[1].set_title("シナリオ別の納期遵守率")
axes[1].set_xlabel("シナリオ")
axes[1].set_ylabel("納期遵守率")
for ax in axes:
    ax.grid(True, alpha=.3)
    ax.tick_params(axis="x", rotation=20)
axes[1].axhline(.95, color="red", linestyle="--", label="目標95%")
axes[1].legend()
plt.tight_layout()
plt.show()

結果の読み取り

通常時の首位案が複合ストレスでも首位とは限らず、納期目標を割るシナリオもあります。経営判断では「最頻シナリオの期待利益」「最悪時の損失」「回復策」をセットで示します。What-ifは未来の予言ではなく、前提が変わったときの意思決定境界を共有する道具です。

7. 対象ノックを通して見える実務上の示唆

1. 目的関数の合意が最初：アルゴリズムより先に、利益、納期、WIP、安全の優先順位を定めます。
2. 探索と検証を分ける：探索時の高得点を、長期間・別シード・ストレス条件で再評価します。
3. 手法は計算費用と変数構造で選ぶ：候補が少なければグリッド、離散組合せならGA/SA、高価な評価ならベイズ最適化が候補です。
4. 固定条件と動的方策を区別する：設備投資は最適化、日々の応援判断はRLなど、意思決定の時間軸を揃えます。
5. モデルを運用資産にする：実績更新、版管理、監視、再計算の責任者がいて初めてデジタルツインになります。

8. 実務導入する場合に必要なこと

• 目的KPI、制約、承認者、利用頻度を記した意思決定定義書
• 製品・設備・工程・停止理由のマスタ整備と時刻同期
• 実績期間を分けた妥当性確認、専門家レビュー、感度分析
• 入力データ、コード、依存環境、乱数シード、実行結果の版管理
• 現場で守れる行動制約、異常時のフォールバック、手動上書き
• 導入後のKPI差分、モデル劣化、データ欠測を監視する運用設計

最適化結果は命令ではなく、前提とリスクを伴う意思決定材料です。安全・品質に関する制約は利益ペナルティへ置き換えず、原則として実行可能性を規定するハード制約にします。

9. まとめ

No.081〜No.090では、シミュレーションの目的設計から探索、資源配分、動的制御、実績更新、シナリオ比較までを一つの架空ラインでつなぎました。実務価値は「最適値を1つ出すこと」ではなく、複数部門が同じ前提で案を比較し、条件変化時に再計算できる意思決定プロセスを作ることにあります。

10. 法人向けのご相談

数理工房では、製造業向けのシミュレーション設計、数理最適化、需要・故障の不確実性評価、PoCから運用基盤化までをご支援しています。課題がまだ数式になっていない段階でも、意思決定とデータの整理からご相談いただけます。

📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。