100本ノック / 確率統計 / 確率・統計マーケティング応用100本ノック
製造業シミュレーションをPythonで学ぶ｜モンテカルロからデジタルツインまで10本ノック
変動を先読みする工場経営：産業用ポンプ工場のシミュレーション10本ノック
本記事では、架空の産業用ポンプ工場を題材に、モンテカルロ、離散イベント、エージェントベース、待ち行列、システムダイナミクス、在庫、生産ライン、サプライチェーン、需要、デジタルツインの10テーマを一続きの意思決定として扱います。
目的は手法の実装自体ではなく、「利益計画はどこまで下振れするか」「設備・要員・在庫をどこへ配分するか」「供給途絶や需要変動にどう備えるか」を、再現可能な仮想実験で説明できるようにすることです。掲載データは Python で生成し、外部データには依存しません。
[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。
はじめに：この記事で扱う製造業の実務課題
対象工場は、標準型と高耐圧型の産業用ポンプを、加工・組立・検査して法人顧客へ出荷しています。需要、加工時間、設備故障、部材の納期は日々変動します。生産管理者は、月次計画だけでなく、日々の仕掛、部材在庫、応援要員、特急輸送、顧客への納期回答を同時に判断しなければなりません。
シミュレーションは未来を一点で言い当てる道具ではありません。現実を意思決定に必要な範囲でモデル化し、まだ実行していない施策を安全に比較する「仮想実験」です。本記事では平均値に加え、分布、時間順序、相互作用、フィードバックを段階的にモデルへ入れます。
現場でよくある状況
- 年間利益は黒字計画だが、需要減・材料高・停止が重なる下振れを説明できない
- 各設備の平均稼働率は分かるが、仕掛がどこで何分待つか分からない
- 保全員、作業者、搬送車、仕入先が局所判断し、全体KPIへの影響が見えない
- 安全在庫を一律に増やし、欠品は減ったが在庫金額と廃棄が増えた
- 需要予測、生産計画、調達計画の前提が部門ごとに異なる
- デジタルツインを導入したが、現場実績でモデルを更新する運用がない
共通する課題は、変動と時間的な依存関係を平均値へ押し込めていることです。
なぜこの問題は判断が難しいのか
製造システムには非線形性があります。稼働率が100%へ近づくと待ち時間は急増し、一つの部材欠品が後工程の停止へ波及します。発注量を増やせば欠品は減りますが、在庫費用と陳腐化リスクは増えます。また、現場の主体は同じ情報を持たず、それぞれ異なるルールで行動します。
したがって、平均利益や平均リードタイムだけでは不十分です。本記事では、利益の下方分位点、納期遵守率、待ち時間、仕掛在庫、充足率、バックログ、回復時間などを併記し、期待値・リスク・応答速度・費用を同時に評価します。
今回扱うノックの全体像
| No. | テーマ | 主な問い | 主な出力 |
|---|---|---|---|
| 081 | モンテカルロ | 利益はどこまで下振れするか | 利益分布、赤字確率 |
| 082 | 離散イベント | ジョブがいつ流れ、どこで待つか | イベントログ、リードタイム |
| 083 | エージェントベース | 個々の保全判断が全体へどう波及するか | 停止率、予防保全率 |
| 084 | 待ち行列 | 検査員を何名配置すべきか | 待ち時間、サービス水準 |
| 085 | システムダイナミクス | 受注・能力・バックログはどう循環するか | 在庫・受注残の推移 |
| 086 | 在庫シミュレーション | 発注政策をどう選ぶか | 充足率、在庫費用 |
| 087 | 生産ラインシミュレーション | 改善投資をどの工程へ向けるか | スループット、工程利用率 |
| 088 | サプライチェーンシミュレーション | 供給途絶へどう備えるか | 欠品、総費用、回復時間 |
| 089 | 需要シミュレーション | S&OPで能力案をどう比較するか | 月間需要分布、超過確率 |
| 090 | デジタルツイン | 実績でモデルをどう更新するか | 状態更新、シナリオ比較 |
Python 環境の準備
NumPy は乱数と数値計算、pandas はイベント・KPIの表現、matplotlib は可視化に利用します。乱数生成器は用途ごとに default_rng(SEED + 番号) として固定し、コードを再実行しても同じ結果になるようにします。グラフはすべて matplotlib で作成します。
import heapq
import platform
import japanize_matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import display
SEED = 42
pd.set_option("display.max_columns", 20)
pd.set_option("display.float_format", lambda x: f"{x:,.2f}")
print(f"Python: {platform.python_version()}")
print(f"NumPy: {np.__version__} / pandas: {pd.__version__}")
print(f"乱数seed: {SEED}")Python: 3.13.1
NumPy: 2.5.1 / pandas: 3.0.3
乱数seed: 42架空データの作成
標準型 P-100 と高耐圧型 P-200 の2製品を想定します。原価・価格・需要、主要3工程の標準時間、日次能力を共通マスタとして用意します。個別ノックでは、このマスタを基準に故障、需要、納期の変動を加えます。
実務では、標準時間が正味時間か余裕込みか、需要が受注日か希望納期日か、欠品が受注残になるか失注になるかを最初に定義します。シミュレーションの精度以前に、データの意味の不一致が結論を変えるためです。
products = pd.DataFrame({
"製品": ["P-100 標準型", "P-200 高耐圧型"],
"平均月間需要_台": [620, 310],
"販売価格_万円": [18.0, 27.0],
"材料費_万円": [8.2, 12.5],
"加工_分": [18, 26],
"組立_分": [22, 34],
"検査_分": [14, 24],
}).set_index("製品")
resources = pd.DataFrame({
"工程": ["加工", "組立", "検査"],
"設備・要員数": [2, 3, 2],
"日次正味時間_分": [450, 450, 450],
"変動係数": [0.18, 0.12, 0.22],
}).set_index("工程")
print("製品マスタ")
display(products)
print("工程能力マスタ")
display(resources)製品マスタ| 平均月間需要_台 | 販売価格_万円 | 材料費_万円 | 加工_分 | 組立_分 | 検査_分 | |
|---|---|---|---|---|---|---|
| 製品 | ||||||
| P-100 標準型 | 620 | 18.00 | 8.20 | 18 | 22 | 14 |
| P-200 高耐圧型 | 310 | 27.00 | 12.50 | 26 | 34 | 24 |
工程能力マスタ| 設備・要員数 | 日次正味時間_分 | 変動係数 | |
|---|---|---|---|
| 工程 | |||
| 加工 | 2 | 450 | 0.18 |
| 組立 | 3 | 450 | 0.12 |
| 検査 | 2 | 450 | 0.22 |
No.081：モンテカルロ — 年間利益の下振れを確率で測る
実務での意味
予算の売上・材料単価・停止時間を一点で置くと、計画利益は一つに決まります。しかし経営判断では、平均利益だけでなく、赤字確率や厳しいケースの損失額が必要です。設備投資枠、運転資金、長期契約価格の検討に直結します。
分析・モデル化の考え方
シナリオ の年間利益を
とします。 は販売数量、 は価格、 は材料費、 は変動加工費、 は固定費、 は停止損失です。入力変数を確率分布から繰り返し生成し、 と5%分位点（95%のケースが上回る利益）を評価します。
Pythonで確認する
mc_rng = np.random.default_rng(SEED + 81)
n_mc = 20_000
annual_demand = np.maximum(0, mc_rng.normal(11_200, 1_250, n_mc))
material_cost = mc_rng.triangular(9.2, 10.0, 12.4, n_mc) # 万円/台
downtime_days = mc_rng.poisson(8, n_mc)
selling_price = 21.2
variable_cost = 3.6
fixed_cost = 73_000
capacity = 11_800
sales = np.minimum(annual_demand, capacity - downtime_days * 28)
profit = sales * (selling_price - material_cost - variable_cost) - fixed_cost - downtime_days * 55
profit_kpi = pd.Series({
"平均利益_万円": profit.mean(),
"利益5%分位_万円": np.quantile(profit, 0.05),
"利益95%分位_万円": np.quantile(profit, 0.95),
"赤字確率": (profit < 0).mean(),
}, name="値")
display(profit_kpi.to_frame())
fig, ax = plt.subplots(figsize=(8, 3.8))
ax.hist(profit, bins=45, color="#4472C4", edgecolor="white")
ax.axvline(0, color="#C00000", linestyle="--", label="損益分岐")
ax.axvline(np.quantile(profit, 0.05), color="#ED7D31", linestyle=":", label="5%分位")
ax.set_title("需要・材料費・停止を反映した年間利益分布")
ax.set_xlabel("年間利益（万円）")
ax.set_ylabel("シナリオ数")
ax.grid(axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 値 | |
|---|---|
| 平均利益_万円 | 3,384.79 |
| 利益5%分位_万円 | -13,233.85 |
| 利益95%分位_万円 | 17,770.93 |
| 赤字確率 | 0.36 |
結果の読み取り
平均利益がプラスでも、0より左側に分布があれば赤字リスクは残ります。5%分位点は資金余力や投資上限を議論する保守的な基準です。ただし、分布の形は入力仮定に依存します。実務では需要と材料価格の相関、価格転嫁の時差、能力上限を実績で推定し、赤字シナリオの入力条件も併せて確認します。
No.082：離散イベント — 工程の競合とジョブの待ちを再現する
実務での意味
製品は加工完了後に組立へ進み、設備が空くまで待ちます。このように状態が「到着」「加工完了」などのイベント時点で変わるシステムでは、時間を細かく刻むより離散イベントシミュレーションが効率的です。納期回答、仕掛削減、段取り順の検討に使えます。
分析・モデル化の考え方
ジョブ の工程 の開始・完了時刻を
とします。 は工程への到着、 は資源が空く時刻、 は処理時間です。イベントを時刻順に優先度付きキューで処理し、イベントログから待ち時間とリードタイムを集計します。
Pythonで確認する
des_rng = np.random.default_rng(SEED + 82)
routes = ["加工", "組立", "検査"]
resource_ready = {m: 0.0 for m in routes}
event_queue = []
log = []
n_jobs = 18
arrival_times = np.cumsum(des_rng.exponential(16, n_jobs))
for j, arrival in enumerate(arrival_times):
heapq.heappush(event_queue, (arrival, j, 0))
while event_queue:
arrival, job, step = heapq.heappop(event_queue)
machine = routes[step]
means = [20, 27, 18]
duration = des_rng.lognormal(np.log(means[step]), 0.18)
start = max(arrival, resource_ready[machine])
finish = start + duration
resource_ready[machine] = finish
log.append([job, machine, step + 1, arrival, start, finish, start - arrival])
if step + 1 < len(routes):
heapq.heappush(event_queue, (finish, job, step + 1))
event_log = pd.DataFrame(log, columns=["job", "工程", "工程順", "到着", "開始", "完了", "待ち時間"])
job_kpi = event_log.groupby("job").agg(投入=("到着", "min"), 完了=("完了", "max"), 総待ち=("待ち時間", "sum"))
job_kpi["リードタイム"] = job_kpi["完了"] - job_kpi["投入"]
display(event_log.head(10))
display(job_kpi.describe().loc[["mean", "50%", "max"]])
fig, ax = plt.subplots(figsize=(8, 3.8))
for i, machine in enumerate(routes):
part = event_log[event_log["工程"] == machine]
ax.scatter(part["開始"], [i] * len(part), s=part["待ち時間"] * 5 + 20, label=machine)
ax.set_yticks(range(len(routes)), routes)
ax.set_title("工程別の処理開始時刻（点の大きさ＝直前待ち時間）")
ax.set_xlabel("シミュレーション時刻（分）")
ax.set_ylabel("工程")
ax.grid(axis="x", alpha=0.3)
plt.tight_layout()
plt.show()| job | 工程 | 工程順 | 到着 | 開始 | 完了 | 待ち時間 | |
|---|---|---|---|---|---|---|---|
| 0 | 0 | 加工 | 1 | 17.26 | 17.26 | 41.15 | 0.00 |
| 1 | 1 | 加工 | 1 | 36.88 | 41.15 | 62.93 | 4.27 |
| 2 | 0 | 組立 | 2 | 41.15 | 41.15 | 75.09 | 0.00 |
| 3 | 1 | 組立 | 2 | 62.93 | 75.09 | 102.12 | 12.16 |
| 4 | 2 | 加工 | 1 | 73.18 | 73.18 | 91.40 | 0.00 |
| 5 | 0 | 検査 | 3 | 75.09 | 75.09 | 93.83 | 0.00 |
| 6 | 3 | 加工 | 1 | 81.55 | 91.40 | 112.45 | 9.85 |
| 7 | 4 | 加工 | 1 | 87.67 | 112.45 | 128.27 | 24.78 |
| 8 | 2 | 組立 | 2 | 91.40 | 102.12 | 127.52 | 10.72 |
| 9 | 1 | 検査 | 3 | 102.12 | 102.12 | 117.42 | 0.00 |
| 投入 | 完了 | 総待ち | リードタイム | |
|---|---|---|---|---|
| mean | 177.21 | 322.24 | 79.07 | 145.04 |
| 50% | 182.37 | 325.33 | 77.44 | 144.69 |
| max | 379.36 | 555.73 | 176.19 | 233.52 |
結果の読み取り
同じ標準時間でも、到着の重なりによりジョブ別リードタイムは異なります。点が大きい工程・時間帯は待ちが蓄積した箇所です。平均待ちだけでなく最大値とジョブ別内訳を見て、増員、優先順、ロット分割の候補を作ります。本番モデルでは複数台設備、故障、段取り、休憩、手直しをイベントとして追加します。
No.083：エージェントベース — 保全員の局所判断と設備群の挙動を捉える
実務での意味
設備ごとの劣化状態は異なり、保全員は限られた時間で巡回します。各設備が警告を出し、保全員が優先順位を判断する相互作用は、集計式だけでは表しにくい問題です。予防保全ルール、巡回能力、アラーム閾値の設計に使えます。
分析・モデル化の考え方
各設備を状態 を持つエージェントとし、毎期の劣化 、保全による回復 を
で更新します。保全員エージェントは、閾値未満の設備を健康度の低い順に1日2台まで処置します。閾値を変えて故障停止と保全回数のトレードオフを比較します。
Pythonで確認する
def run_agents(threshold, seed, days=120, n_machines=18):
agent_rng = np.random.default_rng(seed)
health = agent_rng.uniform(0.65, 1.0, n_machines)
failures = maintenance = lost_hours = 0
history = []
for day in range(days):
health -= agent_rng.gamma(shape=2.0, scale=0.015, size=n_machines)
failed = health <= 0.12
failures += failed.sum()
lost_hours += failed.sum() * 7.5
health[failed] = 0.58 # 事後修理
targets = np.where((health < threshold) & (~failed))[0]
targets = targets[np.argsort(health[targets])][:1]
maintenance += len(targets)
lost_hours += len(targets) * 1.2
health[targets] = np.minimum(1.0, health[targets] + 0.38)
history.append([day, health.mean(), failures, maintenance])
return failures, maintenance, lost_hours, pd.DataFrame(history, columns=["日", "平均健康度", "累積故障", "累積予防保全"])
agent_rows = []
histories = {}
for threshold in [0.30, 0.45, 0.60]:
# 同じ劣化乱数で保全ルールだけを比較する（共通乱数法）
f, m, lost, hist = run_agents(threshold, SEED + 83)
agent_rows.append([threshold, f, m, lost])
histories[threshold] = hist
agent_result = pd.DataFrame(agent_rows, columns=["保全閾値", "故障件数", "予防保全件数", "停止換算時間"])
display(agent_result)
fig, ax = plt.subplots(figsize=(8, 3.8))
for threshold, hist in histories.items():
ax.plot(hist["日"], hist["平均健康度"], label=f"閾値 {threshold:.2f}")
ax.set_title("保全ルール別の設備エージェント平均健康度")
ax.set_xlabel("経過日")
ax.set_ylabel("平均健康度")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 保全閾値 | 故障件数 | 予防保全件数 | 停止換算時間 | |
|---|---|---|---|---|
| 0 | 0.30 | 32 | 109 | 370.80 |
| 1 | 0.45 | 28 | 114 | 346.80 |
| 2 | 0.60 | 26 | 119 | 337.80 |
結果の読み取り
高い閾値は早めの予防保全を増やし、故障を抑える一方で計画停止を増やします。評価は故障件数だけでなく停止換算時間や保全費用で行います。結果には確率変動があるため、本来は各ルールを複数seedで反復します。実務では設備重要度、部品在庫、保全員スキル、同時故障時の優先ルールもエージェント属性にします。
No.084：待ち行列 — 最終検査の要員数をサービス水準で決める
実務での意味
検査工程へ製品が不規則に到着すると、平均処理能力が平均到着量を上回っていても待ちが発生します。検査員を増やせば待ちは減りますが、労務費は増えます。出荷締切までの検査完了率を使うと、要員案を納期サービス水準で比較できます。
分析・モデル化の考え方
到着率を 、1人当たり処理率を 、窓口数を とすると利用率は
です。 では待ち行列が長期的に発散します。ここでは指数分布の到着・処理を仮定した M/M/c をシミュレーションし、平均待ちと「10分以内に検査開始」の割合を比較します。
Pythonで確認する
def simulate_queue(servers, seed, n=3000, arrival_rate=5.2, service_rate=3.0):
queue_rng = np.random.default_rng(seed)
arrivals = np.cumsum(queue_rng.exponential(60 / arrival_rate, n))
service = queue_rng.exponential(60 / service_rate, n)
ready = np.zeros(servers)
waits = np.zeros(n)
for i, (arrival, duration) in enumerate(zip(arrivals, service)):
k = np.argmin(ready)
start = max(arrival, ready[k])
waits[i] = start - arrival
ready[k] = start + duration
return waits
queue_rows = []
queue_samples = {}
for c in [2, 3, 4]:
waits = simulate_queue(c, SEED + 840 + c)
queue_samples[c] = waits
queue_rows.append([c, 5.2 / (c * 3.0), waits.mean(), np.quantile(waits, 0.95), (waits <= 10).mean()])
queue_result = pd.DataFrame(queue_rows, columns=["検査員数", "利用率", "平均待ち_分", "待ち95%分位_分", "10分以内開始率"])
display(queue_result)
fig, ax = plt.subplots(figsize=(8, 3.8))
ax.boxplot([queue_samples[c] for c in [2, 3, 4]], tick_labels=["2名", "3名", "4名"], showfliers=False)
ax.axhline(10, color="#C00000", linestyle="--", label="10分基準")
ax.set_title("検査員数別の待ち時間分布")
ax.set_xlabel("検査員数")
ax.set_ylabel("待ち時間（分）")
ax.grid(axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 検査員数 | 利用率 | 平均待ち_分 | 待ち95%分位_分 | 10分以内開始率 | |
|---|---|---|---|---|---|
| 0 | 2 | 0.87 | 48.88 | 156.23 | 0.31 |
| 1 | 3 | 0.58 | 5.24 | 29.70 | 0.82 |
| 2 | 4 | 0.43 | 1.08 | 8.09 | 0.97 |
結果の読み取り
利用率が高い案では待ち時間の上側が大きく伸びます。平均待ちだけを見ず、95%分位や10分以内開始率を出荷締切と結び付けます。4名案の改善幅が小さければ、常時増員ではなく繁忙時間帯だけの応援が候補です。到着がロット状、処理時間が製品別なら、実績分布を直接再標本化します。
No.085：システムダイナミクス — 受注残と能力調整のフィードバックを見る
実務での意味
受注が増えると受注残が積み上がり、残業・応援で出荷能力を増やします。しかし能力調整には遅れがあり、増やし過ぎると在庫が膨らみ、その後に削減する振動が起こり得ます。S&OP、人員計画、増産判断の時間軸を揃えるために使えます。
分析・モデル化の考え方
在庫 、受注残 をストック、生産・出荷をフローとし、1日刻みの差分方程式で
と更新します。能力 は目標受注残との差に反応しますが、調整係数により遅れて動くとします。急反応と平滑反応を比較します。
Pythonで確認する
def system_dynamics(adjustment, days=120):
demand = np.full(days, 42.0)
demand[25:65] = 58.0
inventory = np.zeros(days)
backlog = np.zeros(days)
capacity = np.zeros(days)
shipment = np.zeros(days)
inventory[0], backlog[0], capacity[0] = 90, 10, 44
for t in range(days - 1):
production = capacity[t]
shipment[t] = min(inventory[t] + production, demand[t] + backlog[t])
inventory[t + 1] = max(0, inventory[t] + production - shipment[t])
backlog[t + 1] = max(0, backlog[t] + demand[t] - shipment[t])
target_capacity = 44 + 0.22 * (backlog[t] - 20) - 0.08 * (inventory[t] - 80)
capacity[t + 1] = np.clip(capacity[t] + adjustment * (target_capacity - capacity[t]), 30, 70)
return pd.DataFrame({"日": np.arange(days), "需要": demand, "在庫": inventory, "受注残": backlog, "能力": capacity})
fast_sd = system_dynamics(0.55)
smooth_sd = system_dynamics(0.12)
sd_kpi = pd.DataFrame({
"調整方針": ["急反応", "平滑反応"],
"最大受注残": [fast_sd["受注残"].max(), smooth_sd["受注残"].max()],
"最大在庫": [fast_sd["在庫"].max(), smooth_sd["在庫"].max()],
"最大日次能力変更": [fast_sd["能力"].diff().abs().max(), smooth_sd["能力"].diff().abs().max()],
})
display(sd_kpi)
fig, axes = plt.subplots(1, 2, figsize=(10, 3.8), sharey=True)
for ax, df, title in zip(axes, [fast_sd, smooth_sd], ["急反応", "平滑反応"]):
ax.plot(df["日"], df["在庫"], label="完成品在庫")
ax.plot(df["日"], df["受注残"], label="受注残")
ax.set_title(f"能力調整：{title}")
ax.set_xlabel("日")
ax.set_ylabel("台")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 調整方針 | 最大受注残 | 最大在庫 | 最大日次能力変更 | |
|---|---|---|---|---|
| 0 | 急反応 | 58.05 | 90.00 | 2.97 |
| 1 | 平滑反応 | 110.36 | 91.40 | 1.86 |
結果の読み取り
急反応は受注残を早く抑えられる一方、能力変更と在庫の振れが大きくなりやすい方針です。平滑反応は運用を安定させますが、需要急増時の顧客待ちを増やします。どちらを選ぶかは、納期遅延費用、残業・応援の変更費用、在庫費用で決まります。実務では採用・教育・設備調達の遅れを明示的な遅延として組み込みます。
No.086：在庫シミュレーション — 発注点と発注量を費用・充足率で比較する
実務での意味
主要シール部品の欠品はラインを止めますが、過剰在庫は資金を固定し、設計変更時の廃棄を増やします。需要と調達リードタイムが変動する環境では、平均値の計算だけでなく、発注政策を日次運用として再現する必要があります。
分析・モデル化の考え方
在庫ポジションを 手持在庫＋発注残－受注残とし、 なら 個を発注する 政策を評価します。日次費用は
で、 は保管費、 は欠品費、 は発注費です。共通の需要・納期シナリオで3政策を比較します。
Pythonで確認する
def inventory_sim(R, Q, seed, days=365):
inv_rng = np.random.default_rng(seed)
demand = inv_rng.poisson(18, days)
lead_times = inv_rng.integers(3, 9, days)
on_hand, backlog = 130, 0
pipeline = []
total_demand = filled = holding = shortage = order_count = 0
history = []
for day in range(days):
arrivals = sum(q for due, q in pipeline if due == day)
pipeline = [(due, q) for due, q in pipeline if due > day]
on_hand += arrivals
requested = backlog + demand[day]
shipped = min(on_hand, requested)
on_hand -= shipped
backlog = requested - shipped
total_demand += demand[day]
filled += min(demand[day], max(0, shipped - max(0, requested - demand[day])))
inventory_position = on_hand + sum(q for _, q in pipeline) - backlog
if inventory_position <= R:
pipeline.append((day + int(lead_times[day]), Q))
order_count += 1
holding += on_hand
shortage += backlog
history.append([day, on_hand, backlog, inventory_position])
cost = holding * 35 + shortage * 1200 + order_count * 2500
return filled / total_demand, cost, np.mean([x[1] for x in history]), pd.DataFrame(history, columns=["日", "手持在庫", "受注残", "在庫ポジション"])
policies = [(70, 120), (100, 140), (140, 180)]
inv_rows, inv_hist = [], {}
for R, Q in policies:
fill, cost, avg_inv, hist = inventory_sim(R, Q, SEED + 86)
name = f"R={R}, Q={Q}"
inv_rows.append([name, fill, avg_inv, cost])
inv_hist[name] = hist
inventory_result = pd.DataFrame(inv_rows, columns=["政策", "即納充足率", "平均手持在庫", "年間関連費用_円"])
display(inventory_result)
fig, ax = plt.subplots(figsize=(8, 3.8))
for name, hist in inv_hist.items():
ax.plot(hist["日"], hist["手持在庫"], label=name, alpha=0.8)
ax.set_title("発注政策別の手持在庫推移")
ax.set_xlabel("日")
ax.set_ylabel("手持在庫（個）")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 政策 | 即納充足率 | 平均手持在庫 | 年間関連費用_円 | |
|---|---|---|---|---|
| 0 | R=70, Q=120 | 0.78 | 35.92 | 3951485 |
| 1 | R=100, Q=140 | 0.94 | 75.34 | 1807165 |
| 2 | R=140, Q=180 | 1.00 | 135.67 | 1856865 |
結果の読み取り
発注点と発注量を増やすと欠品は抑えやすい一方、平均在庫と保管費が増えます。費用最小案が会社の目標充足率を満たすかを先に確認し、未達ならサービス水準を制約として政策を選びます。実務ではロット制約、休日カレンダー、発注済み残、最小発注量、欠品時の代替部品も入れます。
No.087：生産ラインシミュレーション — ボトルネック改善案を比較する
実務での意味
加工・組立・検査のどこを改善すれば出荷量が増えるかは、単純な標準時間の比較だけでは決まりません。上流停止で下流が空き、下流渋滞で仕掛が増えるためです。設備増強、サイクル短縮、予防保全の投資効果検証に使えます。
分析・モデル化の考え方
直列3工程の完了時刻をジョブ順に
で更新します。処理時間 は対数正規分布、故障停止は一定確率で付加します。基準、組立10%短縮、検査増強の3案を同じ乱数系列で比較します。
Pythonで確認する
def line_sim(scenario, seed, jobs=240):
line_rng = np.random.default_rng(seed)
means = np.array([18.0, 25.0, 21.0])
if scenario == "組立10%短縮":
means[1] *= 0.90
if scenario == "検査15%短縮":
means[2] *= 0.85
durations = line_rng.lognormal(np.log(means), [0.18, 0.14, 0.22], size=(jobs, 3))
failures = line_rng.random((jobs, 3)) < [0.025, 0.015, 0.035]
durations += failures * line_rng.uniform(20, 55, size=(jobs, 3))
completion = np.zeros((jobs, 3))
busy = durations.sum(axis=0)
for j in range(jobs):
for k in range(3):
prev_job = completion[j - 1, k] if j else 0
prev_stage = completion[j, k - 1] if k else 0
completion[j, k] = max(prev_job, prev_stage) + durations[j, k]
makespan = completion[-1, -1]
return makespan, jobs / (makespan / 450), busy / makespan, completion
line_rows = []
line_completion = {}
for scenario in ["基準", "組立10%短縮", "検査15%短縮"]:
makespan, daily, utilization, completion = line_sim(scenario, SEED + 87)
line_rows.append([scenario, makespan, daily, *utilization])
line_completion[scenario] = completion
line_result = pd.DataFrame(line_rows, columns=["シナリオ", "完了時間_分", "日産換算_台", "加工利用率", "組立利用率", "検査利用率"])
display(line_result)
fig, ax = plt.subplots(figsize=(8, 3.8))
plot_data = line_result.set_index("シナリオ")[["加工利用率", "組立利用率", "検査利用率"]] * 100
plot_data.plot.bar(ax=ax, color=["#4472C4", "#ED7D31", "#70AD47"], rot=0)
ax.set_title("改善シナリオ別の工程利用率")
ax.set_xlabel("シナリオ")
ax.set_ylabel("利用率（%）")
ax.grid(axis="y", alpha=0.3)
plt.tight_layout()
plt.show()| シナリオ | 完了時間_分 | 日産換算_台 | 加工利用率 | 組立利用率 | 検査利用率 | |
|---|---|---|---|---|---|---|
| 0 | 基準 | 6,127.66 | 17.62 | 0.80 | 0.99 | 0.91 |
| 1 | 組立10%短縮 | 5,698.83 | 18.95 | 0.86 | 0.96 | 0.98 |
| 2 | 検査15%短縮 | 6,124.45 | 17.63 | 0.80 | 0.99 | 0.79 |
結果の読み取り
サイクル短縮率と日産向上率は一致しません。制約工程以外を改善しても、全体完了時刻への効果は限定的です。利用率が高い工程を改善候補としつつ、故障による長い停止や工程間仕掛も確認します。投資判断では、複数反復の平均・下方分位点と、増産1台当たりの費用を比較します。
No.088：サプライチェーンシミュレーション — 供給途絶への緩衝策を選ぶ
実務での意味
主要部材の仕入先停止は、数日遅れて工場欠品となり、さらに顧客納期へ波及します。安全在庫、代替調達、特急輸送は回復力を高めますが、平時費用がかかります。BCPを「念のため」ではなく、欠品削減と費用で比較するために使えます。
分析・モデル化の考え方
仕入先から工場への入荷をリードタイム付きイベントとして管理し、工場在庫から日需要を充当します。途絶期間中は通常発注が届かず、対策案では発注点上昇と、欠品時の代替調達を許します。評価指標は充足率、最大受注残、年間総費用、途絶後の回復日数です。
Pythonで確認する
def supply_chain_sim(policy, seed, days=180):
sc_rng = np.random.default_rng(seed)
on_hand, backlog = 280, 0
pipeline = []
total_demand = immediate_total = holding = expedite = 0
max_backlog = 0
history = []
R, Q = ((180, 260) if policy == "標準" else (280, 320))
for day in range(days):
arrivals = sum(q for due, q in pipeline if due == day)
pipeline = [(due, q) for due, q in pipeline if due > day]
on_hand += arrivals
demand = int(sc_rng.poisson(42))
total_demand += demand
old_backlog = backlog
immediate_total += min(demand, max(0, on_hand - old_backlog))
requested = old_backlog + demand
shipped = min(on_hand, requested)
on_hand -= shipped
backlog = requested - shipped
inventory_position = on_hand + sum(q for _, q in pipeline) - backlog
disruption = 55 <= day < 72
if inventory_position <= R and not disruption:
pipeline.append((day + int(sc_rng.integers(5, 10)), Q))
if policy == "BCP対策" and backlog > 80:
emergency = min(120, backlog)
pipeline.append((day + 2, emergency))
expedite += emergency
holding += on_hand
max_backlog = max(max_backlog, backlog)
history.append([day, on_hand, backlog])
cost = holding * 28 + backlog * 1800 + expedite * 650
hist = pd.DataFrame(history, columns=["日", "工場在庫", "受注残"])
after = hist.loc[hist["日"] >= 72]
recovered = after.loc[after["受注残"] <= 5, "日"]
recovery_days = (recovered.iloc[0] - 72) if len(recovered) else np.nan
return immediate_total / total_demand, max_backlog, cost, recovery_days, hist
sc_rows, sc_hist = [], {}
for policy in ["標準", "BCP対策"]:
fill, max_b, cost, recovery, hist = supply_chain_sim(policy, SEED + 88)
sc_rows.append([policy, fill, max_b, cost, recovery])
sc_hist[policy] = hist
supply_result = pd.DataFrame(sc_rows, columns=["方針", "即納充足率", "最大受注残", "総費用_円", "途絶後回復日数"])
display(supply_result)
fig, ax = plt.subplots(figsize=(8, 3.8))
for policy, hist in sc_hist.items():
ax.plot(hist["日"], hist["受注残"], label=policy)
ax.axvspan(55, 72, color="#C00000", alpha=0.12, label="仕入先途絶")
ax.set_title("供給途絶時の受注残推移")
ax.set_xlabel("日")
ax.set_ylabel("受注残（台）")
ax.grid(alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 方針 | 即納充足率 | 最大受注残 | 総費用_円 | 途絶後回復日数 | |
|---|---|---|---|---|---|
| 0 | 標準 | 0.64 | 621 | 277452 | 9 |
| 1 | BCP対策 | 0.88 | 152 | 1448626 | 4 |
結果の読み取り
BCP対策は平時在庫と特急費を増やしますが、途絶時の最大受注残と回復時間を抑えます。平均年間費用だけでなく、重要顧客の停止損失や契約違約金を含めて選びます。一つの固定途絶だけでは結論が偏るため、発生時期・長さ・複数仕入先の同時被災を変えたストレステストが必要です。
No.089：需要シミュレーション — S&OPの能力案を確率で比較する
実務での意味
営業予測が月間1,000台でも、販促効果、季節性、大口案件の受注で実績は上下します。能力を予測平均に合わせると、繁忙月の残業・欠品が増えます。S&OPでは需要シナリオを用い、通常能力・残業・外注の組合せを合意します。
分析・モデル化の考え方
月 の需要を
とし、 は基準需要、 は傾向、 は季節係数、 は販促・大口案件、 は通常変動です。1,000本の年間パスを生成し、能力超過確率と年間不足量を能力案別に評価します。
Pythonで確認する
demand_rng = np.random.default_rng(SEED + 89)
n_paths, months = 3000, 12
season = np.array([0.90, 0.92, 0.98, 1.02, 1.05, 1.08, 0.95, 0.93, 1.00, 1.08, 1.16, 1.22])
trend = np.arange(months) * 8
base = (900 + trend) * season
noise = demand_rng.normal(0, 75, size=(n_paths, months))
campaign = (demand_rng.random((n_paths, months)) < 0.16) * demand_rng.normal(180, 45, (n_paths, months))
large_order = (demand_rng.random((n_paths, months)) < 0.08) * demand_rng.integers(180, 380, (n_paths, months))
demand_paths = np.maximum(0, base + noise + campaign + large_order)
capacity_options = {"通常能力": 1050, "残業枠": 1180, "外注枠込み": 1320}
demand_rows = []
for name, cap in capacity_options.items():
shortage = np.maximum(0, demand_paths - cap)
demand_rows.append([name, cap, (demand_paths > cap).mean(), shortage.sum(axis=1).mean(), np.quantile(shortage.sum(axis=1), 0.95)])
demand_result = pd.DataFrame(demand_rows, columns=["能力案", "月能力_台", "月次能力超過確率", "平均年間不足_台", "年間不足95%分位_台"])
display(demand_result)
q10, q50, q90 = np.quantile(demand_paths, [0.1, 0.5, 0.9], axis=0)
fig, ax = plt.subplots(figsize=(8, 3.8))
ax.fill_between(np.arange(1, 13), q10, q90, color="#4472C4", alpha=0.2, label="需要10–90%範囲")
ax.plot(np.arange(1, 13), q50, color="#4472C4", marker="o", label="需要中央値")
for name, cap in capacity_options.items():
ax.axhline(cap, linestyle="--", linewidth=1, label=name)
ax.set_title("月別需要シナリオと能力案")
ax.set_xlabel("月")
ax.set_ylabel("需要（台/月）")
ax.set_xticks(range(1, 13))
ax.grid(alpha=0.3)
ax.legend(ncol=2)
plt.tight_layout()
plt.show()| 能力案 | 月能力_台 | 月次能力超過確率 | 平均年間不足_台 | 年間不足95%分位_台 | |
|---|---|---|---|---|---|
| 0 | 通常能力 | 1050 | 0.38 | 670.75 | 1,188.59 |
| 1 | 残業枠 | 1180 | 0.18 | 242.43 | 613.50 |
| 2 | 外注枠込み | 1320 | 0.05 | 65.70 | 289.95 |
結果の読み取り
中央値が能力内でも、需要分布の上側では不足します。通常能力を固定費、残業・外注を変動費として、能力超過確率と不足費用の許容範囲を合意します。販促や大口案件は純粋な乱数ではなく営業活動と連動するため、案件確度別シナリオを営業部門と共同で設定し、毎月更新します。
No.090：デジタルツイン — 現場イベントで状態を更新し施策を先回り評価する
実務での意味
デジタルツインは3D表示そのものではなく、現場の現在状態をデータで同期し、その状態から将来シナリオを比較する仕組みです。標準時間が古いままでは、精巧なモデルでも判断を誤ります。日々のサイクル実績でモデルを更新し、当日残数の完了見込みを出します。
分析・モデル化の考え方
工程 の推定サイクル時間を、観測 により指数平滑で
と更新します。最小構成は、(1) 対象KPI、(2) 状態データ、(3) 更新規則、(4)将来シミュレーション、(5)意思決定、(6)実績での検証です。現在の推定値を初期状態として、基準・応援・停止短縮案を比較します。
Pythonで確認する
twin_rng = np.random.default_rng(SEED + 90)
standards = {"加工": 18.0, "組立": 25.0, "検査": 21.0}
events = []
for machine, standard in standards.items():
drift = {"加工": 1.03, "組立": 1.12, "検査": 0.98}[machine]
observed = twin_rng.lognormal(np.log(standard * drift), 0.12, 40)
events.extend([[machine, i + 1, value] for i, value in enumerate(observed)])
event_data = pd.DataFrame(events, columns=["工程", "完了順", "実績CT_分"])
alpha = 0.20
state_rows = []
for machine, group in event_data.groupby("工程", sort=False):
estimate = standards[machine]
for value in group["実績CT_分"]:
estimate = alpha * value + (1 - alpha) * estimate
state_rows.append([machine, standards[machine], estimate, estimate / standards[machine] - 1])
twin_state = pd.DataFrame(state_rows, columns=["工程", "標準CT_分", "更新CT_分", "標準乖離率"]).set_index("工程")
display(twin_state)
def twin_forecast(ct, scenario, seed, remaining=80, reps=3000):
r = np.random.default_rng(seed)
adjusted = ct.copy()
downtime_mean = 38
if scenario == "組立応援":
adjusted["組立"] *= 0.88
if scenario == "停止短縮":
downtime_mean = 22
bottleneck_ct = max(adjusted.values())
finish = remaining * r.lognormal(np.log(bottleneck_ct), 0.06, reps) + r.gamma(2, downtime_mean / 2, reps)
return finish
ct_now = twin_state["更新CT_分"].to_dict()
twin_rows, twin_samples = [], {}
for i, scenario in enumerate(["基準", "組立応援", "停止短縮"]):
samples = twin_forecast(ct_now, scenario, SEED + 900 + i)
twin_samples[scenario] = samples
twin_rows.append([scenario, samples.mean(), np.quantile(samples, 0.90), (samples <= 2100).mean()])
twin_result = pd.DataFrame(twin_rows, columns=["シナリオ", "平均完了見込_分", "完了90%分位_分", "2100分以内完了確率"])
display(twin_result)
fig, ax = plt.subplots(figsize=(8, 3.8))
for scenario, samples in twin_samples.items():
ax.hist(samples, bins=35, alpha=0.35, label=scenario)
ax.axvline(2100, color="#C00000", linestyle="--", label="納期基準")
ax.set_title("更新状態から予測した残り80台の完了時間")
ax.set_xlabel("完了までの時間（分）")
ax.set_ylabel("シミュレーション回数")
ax.grid(axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 標準CT_分 | 更新CT_分 | 標準乖離率 | |
|---|---|---|---|
| 工程 | |||
| 加工 | 18.00 | 18.70 | 0.04 |
| 組立 | 25.00 | 27.68 | 0.11 |
| 検査 | 21.00 | 20.09 | -0.04 |
| シナリオ | 平均完了見込_分 | 完了90%分位_分 | 2100分以内完了確率 | |
|---|---|---|---|---|
| 0 | 基準 | 2,257.37 | 2,433.68 | 0.12 |
| 1 | 組立応援 | 1,987.86 | 2,138.01 | 0.83 |
| 2 | 停止短縮 | 2,239.69 | 2,412.84 | 0.15 |
結果の読み取り
標準乖離率が大きい工程は、劣化、品種構成、作業方法の変化を調べるシグナルです。更新状態から算出した納期内完了確率により、応援と停止短縮を同じ尺度で比較できます。本番では予測誤差を継続監視し、センサー欠損や時刻ずれを検知します。自動提案を導入しても、製造指示の承認権限、モデル停止時の手動運用、変更履歴を明確にします。
対象ノックを通して見える実務上の示唆
- 問いに合う時間表現を選ぶ：総利益の分布ならモンテカルロ、工程順序なら離散イベント、主体間相互作用ならエージェントベースが適します。
- 平均値から分布へ進む：平均利益・平均需要だけでなく、赤字確率、下方分位点、能力超過確率を意思決定に使います。
- 高稼働を目的にしない：待ち行列では利用率上昇に対して待ちが非線形に増えます。納期サービスと費用を併記します。
- 局所改善を全体KPIで評価する：工程短縮、予防保全、安全在庫は、スループット・停止・総費用への効果で比較します。
- 時間遅れを明示する：能力調整、発注、供給途絶の影響には遅れがあります。月次集計だけでは振動と波及を見落とします。
- モデルを現場実績へ閉じる：デジタルツインでは、状態同期、将来比較、実績検証を一つの更新サイクルとして運用します。
実務導入する場合に必要なこと
1. 意思決定・対象範囲・KPIを先に定義する
「工場全体を再現する」から始めず、検査員数、発注点、残業枠、代替調達など、誰がいつ何を決めるかを定義します。KPIには平均だけでなく、納期遵守率、95%分位、最大受注残、総費用を含めます。
2. 入力データと業務ルールを版管理する
品目、工程、設備、シフト、段取り、発注残、停止理由、受注・失注の定義を揃えます。データ抽出日時、単位、欠損補完、標準時間の適用開始日を記録し、再現可能にします。
3. 妥当性確認を段階的に行う
イベント順序や在庫収支の単体確認、過去期間の再現、現場担当者による極端ケースの確認を行います。モデルの数値が実績と合うことに加え、原因と結果の方向が現場知識と整合するかを確認します。
4. 小さな意思決定ループから運用する
1部材・1工程で「データ更新→シナリオ比較→承認→実行→実績評価」を回し、判断時間や欠品削減などの効果を測ります。その後、ERP・MES・設備・調達データとの連携範囲を広げます。
5. 不確実性と責任分界を伝える
予測区間、前提、適用外条件を画面と会議資料に表示します。自動提案の採否、緊急時の上書き権限、モデル停止時の代替手順、監査ログを業務設計へ含めます。
まとめ
No.081〜No.090では、利益リスクのモンテカルロ評価から、工程イベント、設備・保全員の相互作用、待ち行列、需給フィードバック、在庫政策、生産ライン、供給途絶、需要シナリオ、実績連動型デジタルツインまでを扱いました。
シミュレーションの価値は、複雑なモデルを作ることではなく、実行前に選択肢の結果とリスクを比較できることです。小さな対象で収支とイベントを検証し、意思決定に必要な要素だけを段階的に追加します。そして予測と実績の差を残し続けることで、モデルを一度きりの分析から現場の意思決定基盤へ育てられます。
法人向けのご相談
数理工房では、製造業におけるシミュレーション、在庫・生産・サプライチェーン設計、デジタルツインの構想策定から実装・内製化支援までご相談を承ります。
- 工程・物流・保全を対象とした離散イベントシミュレーション
- 需要・供給・利益リスクのモンテカルロ評価とシナリオ設計
- 在庫政策、要員配置、設備投資、BCP施策の比較検証
- ERP・MES・設備データをつなぐデジタルツイン／意思決定基盤
- Python・統計・シミュレーションを扱う法人研修
📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。