100本ノック / 数理最適化 / 数理最適化100本ノック
製造業の確率・ロバスト最適化入門｜需要変動に強い生産計画をPythonで作る
需要変動に負けない生産計画：利益・欠品・頑健性を両立する確率最適化10本ノック
架空の精密部品工場を題材に、需要や歩留まりが確定しない段階で生産量をどう決めるかを検討します。対象は No.071〜No.080（確率最適化・ロバスト最適化） です。期待利益だけでなく、欠品確率、最悪ケース、CVaR、利益とCO2の両立まで、同じデータから意思決定の違いを比較します。
[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。
はじめに：この記事で扱う製造業の実務課題
月次の生産量は、受注がすべて確定する前に決めることがあります。多く作れば欠品を避けられますが、余剰在庫・残業・廃棄が増えます。少なく作れば効率はよく見えても、需要上振れ時に機会損失が発生します。本稿の目的は「予測値を一点で当てる」ことではなく、予測誤差を含めて説明可能な計画を作ることです。
現場でよくある状況
- 営業予測は月ごとに変わり、製造リードタイムより確定が遅い
- 歩留まり、設備停止、特急材の単価も変動する
- 利益最大化だけでなく、納期遵守率、残業、CO2にも目標がある
- 平均需要だけで作った計画が、繁忙月には欠品し、閑散月には滞留する
なぜこの問題は判断が難しいのか
不確実性に対する姿勢で「最良」の計画が変わるためです。平均的な利益を重視するのか、一定確率で需要を満たすのか、観測範囲の最悪ケースを守るのかを先に合意する必要があります。また、過度に保守的な計画は平時のコストを押し上げます。したがって、最適値だけでなく、利益分布、欠品率、前提への感度を並べて判断します。
今回扱うノックの全体像
| No. | テーマ | 現場で答える問い |
|---|---|---|
| 071 | 不確実性を含む最適化 | 平均値だけの計画は何を見落とすか |
| 072 | 確率計画法 | シナリオ別の損益をどう統合するか |
| 073 | シナリオ最適化 | 需要・歩留まりの組合せをどう評価するか |
| 074 | チャンス制約 | 欠品確率を許容範囲に抑えるには |
| 075 | ロバスト最適化 | 想定範囲の最悪需要を守るには |
| 076 | 分布ロバスト最適化 | 分布推定のずれにどう備えるか |
| 077 | 多目的最適化 | 利益とCO2をどう同時に扱うか |
| 078 | パレート最適 | 経営が選べる代替案をどう示すか |
| 079 | リスク尺度とCVaR | 悪いケースの損失をどう測るか |
| 080 | 需要変動を考慮した生産計画 | 方針を比較し推奨計画をどう決めるか |
Python 環境の準備
外部データは使用しません。乱数シードを固定し、同じ結果を再現できるようにします。金額は万円、数量はロット、CO2は kg-CO2 とします。グラフはすべてMatplotlibで作成します。
import platform
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar
from scipy.stats import norm
from IPython.display import display
rng = np.random.default_rng(20260712)
plt.rcParams["figure.figsize"] = (8, 4.5)
plt.rcParams["axes.unicode_minus"] = False
print("Python :", platform.python_version())
print("NumPy :", np.__version__)
print("pandas :", pd.__version__)
print("Matplotlib :", matplotlib.__version__)Python : 3.13.1
NumPy : 2.5.1
pandas : 3.0.3
Matplotlib : 3.11.0架空データの作成
精密部品Aの月次計画を想定します。通常生産量を 、需要を 、歩留まりを とすると良品数は です。売価は 1.8 万円、通常生産費は 0.9 万円、余剰品の処分価値は 0.35 万円、欠品の機会損失は 0.8 万円とします。
シナリオ における利益を
と定義します。ここで です。需要と歩留まりには弱い負の相関を持たせ、繁忙時に歩留まりがやや低下する架空状況を表現します。
n_scenarios = 3000
z_d = rng.normal(size=n_scenarios)
z_y = rng.normal(size=n_scenarios)
demand = np.clip(100 + 18*z_d, 45, 160)
yield_rate = np.clip(0.94 - 0.018*z_d + 0.018*z_y, 0.84, 0.99)
scenarios = pd.DataFrame({"demand": demand, "yield_rate": yield_rate})
price, unit_cost, salvage, shortage_penalty = 1.80, 0.90, 0.35, 0.80
def profit(q, d=demand, y=yield_rate):
good = q * y
sold = np.minimum(good, d)
surplus = np.maximum(good-d, 0)
shortage = np.maximum(d-good, 0)
return price*sold + salvage*surplus - unit_cost*q - shortage_penalty*shortage
display(scenarios.describe(percentiles=[.05, .5, .95]).round(2))
fig, ax = plt.subplots()
ax.scatter(demand, yield_rate, s=10, alpha=.25)
ax.set_title("Demand and yield scenarios")
ax.set_xlabel("Demand [lots/month]")
ax.set_ylabel("Yield rate")
ax.grid(True, alpha=.3)
plt.tight_layout()
plt.show()| demand | yield_rate | |
|---|---|---|
| count | 3000.00 | 3000.00 |
| mean | 99.77 | 0.94 |
| std | 18.01 | 0.03 |
| min | 45.00 | 0.85 |
| 5% | 69.90 | 0.90 |
| 50% | 99.73 | 0.94 |
| 95% | 129.11 | 0.98 |
| max | 160.00 | 0.99 |
No.071：不確実性を含む最適化
実務での意味
平均需要だけで計画すると、需要分布の裾と歩留まり低下を無視します。不確実性を含む最適化では、意思決定前に未知な量を確率変数または区間として扱い、計画の成績を分布で評価します。
分析・モデル化の考え方
まず、平均値だけを代入する確定計画と、複数シナリオで期待利益を最大化する計画を比べます。期待値計画は
です。ただし期待利益が高くても、下方リスクが小さいとは限りません。
Pythonで確認する
q_grid = np.arange(70, 151)
deterministic_profit = [profit(q, np.array([demand.mean()]), np.array([yield_rate.mean()]))[0] for q in q_grid]
expected_profit = np.array([profit(q).mean() for q in q_grid])
q_det = q_grid[np.argmax(deterministic_profit)]
q_sto = q_grid[np.argmax(expected_profit)]
comparison_071 = pd.DataFrame({
"plan": ["mean-value", "stochastic"],
"production_q": [q_det, q_sto],
"expected_profit": [profit(q_det).mean(), profit(q_sto).mean()],
"shortage_probability": [(q_det*yield_rate < demand).mean(), (q_sto*yield_rate < demand).mean()]
})
display(comparison_071.round(3))
fig, ax = plt.subplots()
ax.plot(q_grid, expected_profit, label="Expected profit")
ax.axvline(q_det, color="tab:orange", linestyle="--", label="Mean-value plan")
ax.axvline(q_sto, color="tab:green", linestyle="--", label="Stochastic plan")
ax.set_title("Production quantity and expected profit")
ax.set_xlabel("Production quantity q [lots]")
ax.set_ylabel("Expected profit [10,000 JPY]")
ax.grid(True, alpha=.3)
ax.legend()
plt.tight_layout()
plt.show()| plan | production_q | expected_profit | shortage_probability | |
|---|---|---|---|---|
| 0 | mean-value | 106 | 65.951 | 0.501 |
| 1 | stochastic | 119 | 69.028 | 0.283 |
結果の読み取り
平均値計画と確率計画では推奨量や欠品確率が変わります。差は「予測精度の差」ではなく、非対称な余剰費と欠品費を分布全体で評価した効果です。現場では平均利益と欠品率を必ず併記します。
No.072：確率計画法
実務での意味
需要が判明した後に追加生産や外注ができる場合、事前の通常生産と事後対応を分けると実態に近づきます。これは「今決める量」と「後で調整する量」を区別する二段階確率計画です。
分析・モデル化の考え方
第一段階で通常生産 を決め、需要判明後に不足分を特急調達します。特急単価を とすると、サンプル平均近似により
を計算します。ここでは既存の利益式の欠品ペナルティが特急対応・失注影響をまとめて表すと解釈します。
Pythonで確認する
rows = []
for q in q_grid:
p = profit(q)
rows.append([q, p.mean(), p.std(ddof=1), np.quantile(p, .10)])
table_072 = pd.DataFrame(rows, columns=["q", "mean_profit", "profit_sd", "p10_profit"])
best_072 = table_072.loc[table_072["mean_profit"].idxmax()]
display(table_072.sort_values("mean_profit", ascending=False).head(8).round(2))
print("期待利益最大の生産量:", int(best_072["q"]))| q | mean_profit | profit_sd | p10_profit | |
|---|---|---|---|---|
| 49 | 119 | 69.03 | 17.85 | 43.22 |
| 50 | 120 | 69.03 | 18.20 | 42.75 |
| 48 | 118 | 68.99 | 17.51 | 43.51 |
| 51 | 121 | 68.99 | 18.54 | 42.27 |
| 52 | 122 | 68.93 | 18.88 | 41.75 |
| 47 | 117 | 68.93 | 17.16 | 44.00 |
| 53 | 123 | 68.83 | 19.21 | 41.43 |
| 46 | 116 | 68.82 | 16.81 | 44.42 |
期待利益最大の生産量: 119結果の読み取り
候補量ごとに平均だけでなく標準偏差と10%点を保存すると、期待利益の差がわずかな候補をリスク面から比較できます。実装時はシナリオ発生確率、特急能力、調達リードタイムも制約として加えます。
No.073：シナリオ最適化
実務での意味
「好況」「標準」「不況」のような少数ケースだけでは、需要と歩留まりの組合せを見落とす場合があります。多数の整合したシナリオを作り、計画がどの条件で崩れるかを確認します。
分析・モデル化の考え方
シナリオ最適化では、過去データや予測モデルから を生成し、同じ を全シナリオへ適用します。訓練シナリオで選んだ計画は、別に生成した検証シナリオで評価し、過適合を点検します。
Pythonで確認する
train_sizes = [50, 200, 1000, 3000]
rows = []
for n in train_sizes:
idx = rng.choice(n_scenarios, size=n, replace=False)
train_scores = [profit(q, demand[idx], yield_rate[idx]).mean() for q in q_grid]
q_hat = q_grid[np.argmax(train_scores)]
rows.append([n, q_hat, profit(q_hat).mean(), (q_hat*yield_rate < demand).mean()])
stability_073 = pd.DataFrame(rows, columns=["scenario_count", "selected_q", "validation_profit", "validation_shortage_rate"])
display(stability_073.round(3))
fig, ax = plt.subplots()
ax.plot(stability_073["scenario_count"], stability_073["selected_q"], marker="o")
ax.set_title("Scenario count and selected production quantity")
ax.set_xlabel("Number of training scenarios")
ax.set_ylabel("Selected q [lots]")
ax.grid(True, alpha=.3)
plt.tight_layout()
plt.show()| scenario_count | selected_q | validation_profit | validation_shortage_rate | |
|---|---|---|---|---|
| 0 | 50 | 124 | 68.695 | 0.206 |
| 1 | 200 | 117 | 68.926 | 0.318 |
| 2 | 1000 | 119 | 69.028 | 0.283 |
| 3 | 3000 | 119 | 69.028 | 0.283 |
結果の読み取り
シナリオ数が少ないと、偶然含まれた需要構成で推奨量が動きます。推奨量と検証KPIが十分安定する数を選び、生成ロジック、乱数シード、対象期間を記録します。
No.074：チャンス制約
実務での意味
重要顧客向けには「平均的に足りる」ではなく「少なくとも95%の確率で需要を満たす」といったサービス水準が必要です。
分析・モデル化の考え方
欠品を許す確率を とすると、チャンス制約は
です。ここではシナリオ中の充足率で近似し、制約を満たす最小生産量を求めます。有限サンプルでの推定誤差があるため、実運用では信頼区間や安全余裕も検討します。
Pythonで確認する
service = np.array([(q*yield_rate >= demand).mean() for q in q_grid])
chance_rows = []
for target in [.90, .95, .99]:
feasible = q_grid[service >= target]
q_req = int(feasible[0]) if len(feasible) else np.nan
chance_rows.append([target, q_req, profit(q_req).mean() if np.isfinite(q_req) else np.nan])
chance_074 = pd.DataFrame(chance_rows, columns=["target_service", "minimum_q", "expected_profit"])
display(chance_074.round(2))
fig, ax = plt.subplots()
ax.plot(q_grid, service, label="Estimated service probability")
for target in [.90, .95, .99]: ax.axhline(target, linestyle="--", alpha=.6)
ax.set_title("Production quantity and service probability")
ax.set_xlabel("Production quantity q [lots]")
ax.set_ylabel("P(good quantity >= demand)")
ax.grid(True, alpha=.3)
plt.tight_layout()
plt.show()| target_service | minimum_q | expected_profit | |
|---|---|---|---|
| 0 | 0.90 | 135.0 | 65.61 |
| 1 | 0.95 | 143.0 | 62.15 |
| 2 | 0.99 | NaN | NaN |
結果の読み取り
サービス目標を上げるほど必要量が非線形に増えます。99%を一律に課すのではなく、顧客・品番別の欠品影響と追加費用を比較して目標を定めるべきです。
No.075：ロバスト最適化
実務での意味
予測分布を信頼しにくい新製品や災害対応品では、発生確率より「想定範囲のどこでも守れるか」を優先する場合があります。
分析・モデル化の考え方
需要上限 と歩留まり下限 から、最悪条件でも充足する制約
を置きます。確率を必要としない反面、不確実集合が広すぎると過剰生産になります。
Pythonで確認する
robust_rows = []
for d_upper, y_lower, label in [(125,.91,"moderate"),(140,.88,"strict"),(160,.84,"extreme")]:
q_req = int(np.ceil(d_upper/y_lower))
robust_rows.append([label,d_upper,y_lower,q_req,profit(q_req).mean(),(q_req*yield_rate>=demand).mean()])
robust_075 = pd.DataFrame(robust_rows, columns=["uncertainty_set","demand_upper","yield_lower","q","expected_profit","service_rate"])
display(robust_075.round(2))| uncertainty_set | demand_upper | yield_lower | q | expected_profit | service_rate | |
|---|---|---|---|---|---|---|
| 0 | moderate | 125 | 0.91 | 138 | 64.40 | 0.93 |
| 1 | strict | 140 | 0.88 | 160 | 53.15 | 0.99 |
| 2 | extreme | 160 | 0.84 | 191 | 35.62 | 1.00 |
結果の読み取り
不確実集合を広げるほどサービス率は上がりますが、平均利益は低下し得ます。上限・下限は恣意的に置かず、過去の予測誤差、設備故障履歴、BCP方針から合意します。
No.076：分布ロバスト最適化
実務での意味
確率計画は推定した分布が正しいことに依存します。サンプルが少ない、新製品で構成が変わる、といった場合は分布そのもののずれにも備える必要があります。
分析・モデル化の考え方
分布ロバスト最適化（DRO）は候補となる分布集合 の中で最悪の期待利益を評価します。
ここでは説明用に、基準・需要上振れ・歩留まり悪化という3つの分布モデルを作り、その最小期待利益を最大化します。
Pythonで確認する
models = {
"baseline": (demand, yield_rate),
"demand_shift": (demand+8, yield_rate),
"yield_shift": (demand, np.clip(yield_rate-.025,.80,1.0)),
}
model_profit = pd.DataFrame({name:[profit(q,d,y).mean() for q in q_grid] for name,(d,y) in models.items()}, index=q_grid)
worst_expected = model_profit.min(axis=1)
q_dro = int(worst_expected.idxmax())
display(pd.DataFrame({"model":model_profit.columns,"expected_profit_at_DRO":model_profit.loc[q_dro].values}).round(2))
print("DRO風の推奨生産量:", q_dro, " 最悪モデル期待利益:", round(worst_expected.loc[q_dro],2))
fig, ax = plt.subplots()
for col in model_profit: ax.plot(q_grid, model_profit[col], label=col)
ax.plot(q_grid, worst_expected, color="black", linewidth=2, label="Worst model")
ax.axvline(q_dro, color="tab:red", linestyle="--")
ax.set_title("Expected profit under distribution shifts")
ax.set_xlabel("Production quantity q [lots]")
ax.set_ylabel("Expected profit [10,000 JPY]")
ax.grid(True, alpha=.3)
ax.legend()
plt.tight_layout()
plt.show()| model | expected_profit_at_DRO | |
|---|---|---|
| 0 | baseline | 68.93 |
| 1 | demand_shift | 75.04 |
| 2 | yield_shift | 66.05 |
DRO風の推奨生産量: 122 最悪モデル期待利益: 66.05
結果の読み取り
DRO風計画は一つの推定分布への依存を弱めます。実務では分布集合の大きさが保守性を決めるため、バックテストで「どの程度の分布ずれを覆うか」を検証します。
No.077：多目的最適化
実務での意味
製造計画には利益、納期、CO2、残業など複数の目的があります。すべてを金額換算できない場合、単一の答えを先に決めず、重みを変えた候補を比較します。
分析・モデル化の考え方
通常生産のCO2をロット当たり 5.2 kg、欠品時の特急対応を不足ロット当たり 11 kg とし、期待CO2を計算します。重み付き和
で を変えます。目的値の単位とスケールを明示することが重要です。
Pythonで確認する
def expected_co2(q):
shortage = np.maximum(demand-q*yield_rate,0)
return 5.2*q + 11.0*shortage.mean()
co2 = np.array([expected_co2(q) for q in q_grid])
multi_rows=[]
for lam in [0, .02, .05, .10, .20]:
score=expected_profit-lam*co2
i=np.argmax(score)
multi_rows.append([lam,q_grid[i],expected_profit[i],co2[i],service[i]])
multi_077=pd.DataFrame(multi_rows,columns=["carbon_weight","q","expected_profit","expected_co2","service_rate"])
display(multi_077.round(2))| carbon_weight | q | expected_profit | expected_co2 | service_rate | |
|---|---|---|---|---|---|
| 0 | 0.00 | 119 | 69.03 | 656.43 | 0.72 |
| 1 | 0.02 | 118 | 68.99 | 654.19 | 0.70 |
| 2 | 0.05 | 117 | 68.93 | 652.11 | 0.68 |
| 3 | 0.10 | 114 | 68.54 | 646.79 | 0.64 |
| 4 | 0.20 | 112 | 68.11 | 644.05 | 0.60 |
結果の読み取り
炭素の重みを変えると推奨量とサービス率が変わります。重みは分析者が密室で決めず、炭素価格、顧客要件、SBT等の社内方針と結び付けます。
No.078：パレート最適
実務での意味
利益を改善しつつCO2も減らせる候補は、明らかに劣る候補です。どちらか一方を改善すると他方が悪化する候補だけを残せば、経営判断に必要な選択肢を絞れます。
分析・モデル化の考え方
候補AがB以上の利益、B以下のCO2で、少なくとも一方が厳密に優れるとき、AはBを支配します。どの候補にも支配されない点がパレート最適です。パレート集合は「唯一の正解」ではなく、価値判断を行うための境界です。
Pythonで確認する
pareto=[]
for i,q in enumerate(q_grid):
dominated=np.any((expected_profit>=expected_profit[i]) & (co2<=co2[i]) & ((expected_profit>expected_profit[i]) | (co2<co2[i])))
if not dominated: pareto.append(i)
pareto_078=pd.DataFrame({"q":q_grid[pareto],"expected_profit":expected_profit[pareto],"expected_co2":co2[pareto],"service_rate":service[pareto]})
display(pareto_078.iloc[::max(1,len(pareto_078)//8)].round(2))
fig, ax = plt.subplots()
ax.scatter(co2, expected_profit, alpha=.45, label="Candidates")
ax.plot(co2[pareto], expected_profit[pareto], color="tab:red", marker="o", label="Pareto frontier")
ax.set_title("Profit-CO2 Pareto frontier")
ax.set_xlabel("Expected CO2 [kg-CO2]")
ax.set_ylabel("Expected profit [10,000 JPY]")
ax.grid(True, alpha=.3)
ax.legend()
plt.tight_layout()
plt.show()| q | expected_profit | expected_co2 | service_rate | |
|---|---|---|---|---|
| 0 | 106 | 65.95 | 640.16 | 0.50 |
| 1 | 107 | 66.40 | 640.35 | 0.52 |
| 2 | 108 | 66.82 | 640.75 | 0.54 |
| 3 | 109 | 67.20 | 641.29 | 0.55 |
| 4 | 110 | 67.54 | 642.01 | 0.57 |
| 5 | 111 | 67.85 | 642.94 | 0.59 |
| 6 | 112 | 68.11 | 644.05 | 0.60 |
| 7 | 113 | 68.34 | 645.34 | 0.62 |
| 8 | 114 | 68.54 | 646.79 | 0.64 |
| 9 | 115 | 68.70 | 648.42 | 0.66 |
| 10 | 116 | 68.82 | 650.20 | 0.67 |
| 11 | 117 | 68.93 | 652.11 | 0.68 |
| 12 | 118 | 68.99 | 654.19 | 0.70 |
| 13 | 119 | 69.03 | 656.43 | 0.72 |
結果の読み取り
パレート曲線の傾きから、追加利益のために増えるCO2を読めます。意思決定会議では代表的な3〜5案に絞り、利益・CO2・サービス率を同じ表で提示すると合意しやすくなります。
No.079：リスク尺度とCVaR
実務での意味
平均利益が同じでも、まれに大赤字となる計画と、損失が安定した計画は同じではありません。CVaRは悪いケース群の平均を測り、事業継続上の下方リスクを表します。
分析・モデル化の考え方
利益の下位 分位を とし、本稿では下位 の平均利益を lower-tail CVaR と呼びます。
値が大きいほど悪いケースでも利益を維持します。損失で定義するCVaRとは符号が逆になる点に注意します。
Pythonで確認する
alpha=.10
lower_cvar=[]
for q in q_grid:
p=profit(q); threshold=np.quantile(p,alpha)
lower_cvar.append(p[p<=threshold].mean())
lower_cvar=np.array(lower_cvar)
risk_rows=[]
for gamma in [0,.25,.5,.75,1.0]:
score=(1-gamma)*expected_profit+gamma*lower_cvar
i=np.argmax(score)
risk_rows.append([gamma,q_grid[i],expected_profit[i],lower_cvar[i],service[i]])
risk_079=pd.DataFrame(risk_rows,columns=["risk_weight","q","expected_profit","lower_10pct_CVaR","service_rate"])
display(risk_079.round(2))
fig, ax = plt.subplots()
ax.plot(q_grid, expected_profit, label="Expected profit")
ax.plot(q_grid, lower_cvar, label="Lower 10% CVaR")
ax.set_title("Expected profit and downside profit")
ax.set_xlabel("Production quantity q [lots]")
ax.set_ylabel("Profit [10,000 JPY]")
ax.grid(True, alpha=.3)
ax.legend()
plt.tight_layout()
plt.show()| risk_weight | q | expected_profit | lower_10pct_CVaR | service_rate | |
|---|---|---|---|---|---|
| 0 | 0.00 | 119 | 69.03 | 31.92 | 0.72 |
| 1 | 0.25 | 115 | 68.70 | 33.81 | 0.66 |
| 2 | 0.50 | 110 | 67.54 | 35.68 | 0.57 |
| 3 | 0.75 | 108 | 66.82 | 36.13 | 0.54 |
| 4 | 1.00 | 106 | 65.95 | 36.30 | 0.50 |
結果の読み取り
リスク重みを大きくすると、平均利益を一部譲って下位ケースを改善する計画が選ばれます。CVaRの水準と重みは、資金余力、違約金、重要顧客への影響など経営上の損失許容度から決めます。
No.080：需要変動を考慮した生産計画
実務での意味
最後に、期待利益最大、95%サービス、ロバスト、DRO風、CVaR重視を同じKPIで比較します。方式名から選ぶのではなく、自社の責任範囲に合うKPIから選ぶことが重要です。
分析・モデル化の考え方
各方針で決めた を共通の検証シナリオへ適用し、期待利益、下位10% CVaR、サービス率、余剰量、CO2を評価します。本稿では重要顧客向けで95%サービスを必須とし、その中で期待利益が最大の計画を推奨します。
Pythonで確認する
q_chance=int(chance_074.loc[chance_074.target_service==.95,"minimum_q"].iloc[0])
q_robust=int(robust_075.loc[robust_075.uncertainty_set=="moderate","q"].iloc[0])
q_cvar=int(risk_079.loc[risk_079.risk_weight==1.0,"q"].iloc[0])
plans={"Expected profit":int(q_sto),"95% service":q_chance,"Robust":q_robust,"DRO-like":q_dro,"CVaR-focused":q_cvar}
rows=[]
for name,q in plans.items():
p=profit(q); var=np.quantile(p,.10); good=q*yield_rate
rows.append([name,q,p.mean(),p[p<=var].mean(),(good>=demand).mean(),np.maximum(good-demand,0).mean(),expected_co2(q)])
final_080=pd.DataFrame(rows,columns=["plan","q","expected_profit","lower_10pct_CVaR","service_rate","expected_surplus","expected_co2"])
display(final_080.round(2).sort_values("expected_profit",ascending=False))
eligible=final_080[final_080.service_rate>=.95]
recommended=eligible.loc[eligible.expected_profit.idxmax()]
print(f"推奨: {recommended['plan']} / q={int(recommended['q'])} lots")
fig, ax = plt.subplots()
x=np.arange(len(final_080))
ax.bar(x, final_080["expected_profit"], color="tab:blue", alpha=.75)
ax.set_xticks(x, final_080["plan"], rotation=20, ha="right")
ax.set_title("Expected profit by planning policy")
ax.set_xlabel("Planning policy")
ax.set_ylabel("Expected profit [10,000 JPY]")
ax.grid(True, axis="y", alpha=.3)
plt.tight_layout()
plt.show()| plan | q | expected_profit | lower_10pct_CVaR | service_rate | expected_surplus | expected_co2 | |
|---|---|---|---|---|---|---|---|
| 0 | Expected profit | 119 | 69.03 | 31.92 | 0.72 | 15.54 | 656.43 |
| 3 | DRO-like | 122 | 68.93 | 30.35 | 0.76 | 17.65 | 664.14 |
| 4 | CVaR-focused | 106 | 65.95 | 36.30 | 0.50 | 7.99 | 640.16 |
| 2 | Robust | 138 | 64.40 | 21.40 | 0.93 | 30.64 | 724.82 |
| 1 | 95% service | 143 | 62.15 | 18.60 | 0.95 | 35.08 | 747.89 |
推奨: 95% service / q=143 lots
結果の読み取り
同じデータでも、方針により生産量とKPIが変わります。この例では95%サービスをゲート条件とし、その制約内で期待利益が高い案を推奨しました。実務では候補量を設備能力・最小ロット・段取り単位へ丸め、丸め後に全KPIを再評価します。
対象ノックを通して見える実務上の示唆
- 平均需要は入力の要約であり、意思決定に必要なリスク情報を残していません。
- 確率計画、チャンス制約、ロバスト最適化、DROは、守りたいものと前提が異なります。
- シナリオ数を増やすだけでなく、需要と歩留まりの相関、季節性、構造変化を反映する必要があります。
- サービス率を上げる追加費用と、欠品が生む損失を同じ会議資料で比較します。
- 多目的問題ではパレート候補を示し、目的の重みを経営判断として扱います。
- CVaRは平均では見えない悪いケースを可視化しますが、対象確率と符号の定義を明示します。
実務導入する場合に必要なこと
- 意思決定の境界を決める：対象品番、期間、通常生産・外注・在庫の範囲を定義します。
- 損益を棚卸しする：製造費だけでなく、欠品、特急輸送、廃棄、残業、顧客影響を整理します。
- シナリオを検証する：予測誤差、相関、繁忙期、設備停止を含め、未使用期間でバックテストします。
- リスク方針を合意する：サービス水準、不確実集合、CVaR水準を営業・製造・財務で決めます。
- 現場制約を実装する：能力、段取り、最小ロット、保管、原料、要員をモデルへ追加します。
- 運用を監視する：推奨値、採用値、上書き理由、実績KPI、分布ずれを保存し、再学習条件を定めます。
まとめ
不確実性下の生産計画では、予測を一点へ潰さず、複数シナリオで利益とサービスを評価することが出発点です。確率計画は平均的な成果、チャンス制約は達成確率、ロバスト最適化は想定範囲の保証、DROは分布ずれ、CVaRは下方リスクに焦点を当てます。最適化手法を選ぶ前に、何をどの程度守るのかを定義することが、説明可能で運用できる計画につながります。
法人向けのご相談
数理工房では、需要予測の不確実性を織り込んだ生産・在庫計画、サービス水準の設計、シナリオ生成、既存システムと連携する意思決定ロジックの構築をご支援します。PoCでは、現行ルールとの比較、KPI設計、バックテスト、現場制約の整理から始められます。
📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。