100本ノック / 確率統計 / 確率・統計Python100本ノック
不確実性を設備・在庫・生産能力の判断へつなぐ：製造業シミュレーション10本ノック
不確実性を設備・在庫・生産能力の判断へつなぐ：製造業シミュレーション10本ノック
製造現場では、平均需要、標準時間、平均故障間隔が分かっていても、欠品、滞留、停止がいつ重なるかまでは分かりません。本記事では、架空の精密部品工場を題材に、モンテカルロ積分、待ち行列、在庫、故障、エージェント、マルコフ連鎖、離散イベント、ブラウン運動、確率微分方程式、デジタルツインをPythonで実装します。
目的はシミュレーション手法を並べることではなく、「検査能力を増やすか」「安全在庫を何個持つか」「予防保全をいつ行うか」「設備投資案をどう比較するか」という意思決定へつなげることです。対象は確率・統計Python実装100本ノックの No.091〜No.100 です。
[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。
はじめに：この記事で扱う製造業の実務課題
架空の精密ポンプ部品工場で、翌四半期の生産計画と設備投資を検討します。受注は日々変動し、部品の補充には時間がかかり、設備は確率的に故障します。さらに、加工後の全数検査がボトルネックになり、仕掛品の滞留が納期を圧迫しています。
こうした問題は、平均値だけでは判断できません。複数の不確実性を時間軸上で再現し、KPIの分布、悪化条件、施策間のトレードオフを比較する必要があります。本記事では小さく監査可能なモデルから始め、最後に現場データで更新する簡易デジタルツインへ統合します。
現場でよくある状況
- 平均処理能力は需要を上回るのに、繁忙時だけ検査待ちが急増する
- 安全在庫を経験則で設定し、欠品損失と保有費の比較ができていない
- 平均故障間隔だけを使い、設備の経年劣化や予防保全時期を考慮していない
- 工程ごとの改善案は評価できても、工場全体のリードタイム効果が見えない
- シミュレーションの数字が現場実績から乖離しても、更新方法と責任者が決まっていない
シミュレーションは未来を一点予言するものではありません。不確実な前提の下で、代替案ごとの結果分布を同じ物差しで比較するための意思決定支援です。
なぜこの問題は判断が難しいのか
需要、処理時間、故障、補充リードタイムは確率変数であり、しかも互いの影響が時間を通じて蓄積します。平均需要が平均能力より小さくても、利用率が1に近づけば待ち時間は非線形に増えます。また、欠品率を下げる施策は在庫費を増やし、予防保全を増やす施策は故障を減らす一方で計画停止を増やします。
シミュレーションで推定する期待値は一般に
と表せます。ただし、有限回の試行にはモンテカルロ誤差があり、入力分布が誤っていれば精密な出力も誤ります。したがって、平均だけでなく分位点、信頼区間、感度、再現性、モデル妥当性を確認します。
今回扱うノックの全体像
| No. | テーマ | 製造業での問い |
|---|---|---|
| 091 | モンテカルロ積分 | 需要変動を含む期待欠品費を推定できるか |
| 092 | 待ち行列シミュレーション | 全数検査の待ち時間とSLA超過率はどの程度か |
| 093 | 在庫シミュレーション | 発注点を変えると欠品と在庫費はどう変わるか |
| 094 | 故障シミュレーション | 予防保全周期を何日に設定すべきか |
| 095 | エージェントシミュレーション | 製品の動的振分けは固定振分けより有効か |
| 096 | マルコフ連鎖 | 設備状態の遷移から長期停止率を見積もれるか |
| 097 | 離散イベントシミュレーション | 多工程ラインのボトルネックはどこか |
| 098 | ブラウン運動 | センサードリフトの閾値到達リスクはどれくらいか |
| 099 | 確率微分方程式 | 温度制御の平均回帰と外乱を同時に表せるか |
| 100 | デジタルツイン入門 | 実績で更新したモデルから改善案を比較できるか |
前半で個別の不確実性を扱い、後半で状態・時間・工程間依存を表現し、最後にデータ更新とシナリオ比較を一つの運用ループへまとめます。
Python 環境の準備
NumPyで乱数と数値計算、pandasで表、SciPyで分布と積分、matplotlibで可視化します。japanize_matplotlib は日本語ラベル表示にのみ使用し、seabornや外部データは使いません。再現性のため、分析全体のseedを固定します。
import platform
import heapq
import numpy as np
import pandas as pd
import scipy
from scipy import integrate, stats
import matplotlib
import matplotlib.pyplot as plt
import japanize_matplotlib
SEED = 20260711
rng = np.random.default_rng(SEED)
plt.rcParams["figure.figsize"] = (8, 4.5)
plt.rcParams["axes.unicode_minus"] = False
print(f"Python : {platform.python_version()}")
print(f"NumPy : {np.__version__}")
print(f"pandas : {pd.__version__}")
print(f"SciPy : {scipy.__version__}")
print(f"Matplotlib : {matplotlib.__version__}")
print(f"random seed: {SEED}")Python : 3.13.1
NumPy : 2.5.1
pandas : 3.0.3
SciPy : 1.18.0
Matplotlib : 3.11.0
random seed: 20260711架空データの作成
対象は、1日8時間・年240日稼働する精密ポンプ部品工場です。過去120日分に相当する需要、検査到着間隔、検査時間、設備状態、炉温をPython内で生成します。ここでの「観測データ」は後のデジタルツイン更新にも使います。
真の生成条件をコード内に置けるのは教材だからです。実務では、欠測、停止中の打切り、品種構成、計画停止、勤務体系を整理し、分布仮定を診断してからモデル入力を作ります。
n_history_days = 120
daily_demand = np.maximum(0, rng.normal(102, 18, n_history_days).round()).astype(int)
inspection_interarrival = rng.exponential(scale=5.2, size=1500) # 分
inspection_service = rng.lognormal(mean=np.log(4.4), sigma=0.25, size=1500) # 分
equipment_states = rng.choice(["健全", "劣化", "停止"], size=n_history_days,
p=[0.78, 0.17, 0.05])
furnace_temperature = 180 + rng.normal(0, 1.8, n_history_days)
factory_data = pd.DataFrame({
"項目": ["日次需要", "検査到着間隔", "検査時間", "設備状態", "炉温"],
"観測数": [len(daily_demand), len(inspection_interarrival), len(inspection_service),
len(equipment_states), len(furnace_temperature)],
"要約": [f"平均 {daily_demand.mean():.1f} 個/日",
f"平均 {inspection_interarrival.mean():.2f} 分",
f"平均 {inspection_service.mean():.2f} 分",
f"停止率 {np.mean(equipment_states == '停止'):.1%}",
f"平均 {furnace_temperature.mean():.2f} ℃"],
})
factory_data| 項目 | 観測数 | 要約 | |
|---|---|---|---|
| 0 | 日次需要 | 120 | 平均 101.7 個/日 |
| 1 | 検査到着間隔 | 1500 | 平均 5.05 分 |
| 2 | 検査時間 | 1500 | 平均 4.48 分 |
| 3 | 設備状態 | 120 | 停止率 5.0% |
| 4 | 炉温 | 120 | 平均 179.99 ℃ |
No.091：モンテカルロ積分 — 需要超過による期待損失を見積もる
実務での意味
日産能力を決めるとき、平均需要だけを満たしても需要の上振れで特急対応や失注が生じます。能力 を超えた数量に1個4,500円の損失が生じるとして、期待損失を評価します。
分析・モデル化の考え方
需要 、損失 とすると、期待損失は
です。モンテカルロ法では需要を繰り返し生成して標本平均を求めます。推定標準誤差は概ね で減少するため、精度を2倍にするには約4倍の試行が必要です。
Pythonで確認する
mu_d, sigma_d, capacity, penalty = 102, 18, 125, 4500
n_mc = 100_000
mc_demand = rng.normal(mu_d, sigma_d, n_mc)
mc_loss = penalty * np.maximum(mc_demand - capacity, 0)
mc_estimate = mc_loss.mean()
mc_se = mc_loss.std(ddof=1) / np.sqrt(n_mc)
integral_value, _ = integrate.quad(
lambda d: penalty * (d - capacity) * stats.norm.pdf(d, mu_d, sigma_d),
capacity, np.inf,
)
mc_table = pd.DataFrame({
"方法": ["Monte Carlo", "数値積分"],
"期待損失 [円/日]": [mc_estimate, integral_value],
"Monte Carlo標準誤差 [円]": [mc_se, np.nan],
})
display(mc_table.round(1))
checkpoints = np.unique(np.logspace(2, 5, 80).astype(int))
cumulative_estimates = np.cumsum(mc_loss)[checkpoints - 1] / checkpoints
plt.plot(checkpoints, cumulative_estimates, label="Monte Carlo推定")
plt.axhline(integral_value, color="tab:red", ls="--", label="数値積分")
plt.xscale("log")
plt.title("試行回数と期待欠品損失の収束")
plt.xlabel("試行回数 [回、対数軸]")
plt.ylabel("期待損失 [円/日]")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| 方法 | 期待損失 [円/日] | Monte Carlo標準誤差 [円] | |
|---|---|---|---|
| 0 | Monte Carlo | 3913.0 | 49.6 |
| 1 | 数値積分 | 3865.5 | NaN |
結果の読み取り
試行回数が増えると推定値が数値積分の値へ近づきます。モンテカルロ法の価値は、この1変数なら積分できることではなく、品種構成、能力低下、残業判断などを含む高次元の損失にも同じ考え方を拡張できる点です。期待損失だけでなく95%分位などの尾部リスクも併記し、能力増強費と比較します。
No.092：待ち行列シミュレーション — 全数検査の滞留を評価する
実務での意味
検査工程では、平均処理時間が到着間隔より短くても、到着と処理のばらつきによって待ちが発生します。平均待ち時間と「15分以内に検査開始」というSLAの超過率を見ます。
分析・モデル化の考え方
1台の検査機に製品が到着するGI/G/1型の待ち行列を、到着時刻 、開始時刻 、終了時刻 として
で逐次計算します。定常性を近似するため、立上がり期間を除外します。稼働率が高いほど入力推定誤差の影響も大きくなります。
Pythonで確認する
def simulate_single_queue(interarrival, service, warmup=200):
arrival = np.cumsum(interarrival)
start = np.empty(len(arrival))
finish = np.empty(len(arrival))
previous_finish = 0.0
for i in range(len(arrival)):
start[i] = max(arrival[i], previous_finish)
finish[i] = start[i] + service[i]
previous_finish = finish[i]
waiting = start - arrival
return arrival[warmup:], waiting[warmup:], finish[warmup:] - arrival[warmup:]
arrival_q, waiting_q, lead_q = simulate_single_queue(
inspection_interarrival, inspection_service
)
queue_result = pd.DataFrame({
"KPI": ["推定利用率", "平均待ち時間", "待ち時間95%点", "15分超過率"],
"値": [inspection_service.mean() / inspection_interarrival.mean(),
waiting_q.mean(), np.quantile(waiting_q, 0.95), np.mean(waiting_q > 15)],
"単位": ["比率", "分", "分", "比率"],
})
display(queue_result.round(3))
plt.hist(waiting_q, bins=35, alpha=0.75, edgecolor="white")
plt.axvline(15, color="tab:red", ls="--", label="SLA 15分")
plt.title("全数検査工程の待ち時間分布")
plt.xlabel("待ち時間 [分]")
plt.ylabel("製品数 [個]")
plt.grid(True, axis="y", alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| KPI | 値 | 単位 | |
|---|---|---|---|
| 0 | 推定利用率 | 0.886 | 比率 |
| 1 | 平均待ち時間 | 16.735 | 分 |
| 2 | 待ち時間95%点 | 58.754 | 分 |
| 3 | 15分超過率 | 0.394 | 比率 |
結果の読み取り
平均値に加えて95%点とSLA超過率を見ると、繁忙時の顧客影響を把握できます。待ち時間の分布は右に長く、平均だけでは長い滞留を隠します。能力追加を決める前に、品種別処理時間、段取り、休憩、再検査、到着の時間帯集中をモデルへ反映し、実績分布との再現性を確認します。
No.093：在庫シミュレーション — 発注点と欠品・保有費を比較する
実務での意味
安全在庫を増やせば欠品は減りますが、保管スペース、資金、陳腐化リスクが増えます。発注点候補を同じ需要系列で比較し、サービス率と総費用の両方で判断します。
分析・モデル化の考え方
在庫ポジションが発注点 以下になれば固定量 を発注し、4日後に入荷する 方策を日次で再現します。費用は保有費と未充足需要のペナルティから構成します。候補間に共通乱数を使うと、需要系列の偶然差を抑えて方策差を比較できます。
Pythonで確認する
def simulate_inventory(demand, reorder_point, order_qty=420, lead_time=4,
initial_stock=420, holding_cost=35, shortage_cost=1800):
on_hand = initial_stock
pipeline = []
rows = []
for day, d in enumerate(demand):
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
shipped = min(on_hand, d)
shortage = d - shipped
on_hand -= shipped
position = on_hand + sum(q for _, q in pipeline)
if position <= reorder_point:
pipeline.append((day + lead_time, order_qty))
rows.append((day, d, arrivals, shipped, shortage, on_hand))
df = pd.DataFrame(rows, columns=["日", "需要", "入荷", "出荷", "欠品", "期末在庫"])
cost = holding_cost * df["期末在庫"].sum() + shortage_cost * df["欠品"].sum()
return df, cost
inventory_demand = np.maximum(0, rng.normal(102, 18, 365).round()).astype(int)
inventory_results = []
inventory_paths = {}
for s in [320, 380, 440, 500]:
path, cost = simulate_inventory(inventory_demand, s)
inventory_paths[s] = path
inventory_results.append({
"発注点": s,
"充足率": path["出荷"].sum() / path["需要"].sum(),
"平均期末在庫": path["期末在庫"].mean(),
"欠品数量": path["欠品"].sum(),
"年間総費用 [円]": cost,
})
inventory_table = pd.DataFrame(inventory_results)
display(inventory_table.round({"充足率": 4, "平均期末在庫": 1}))
for s in [320, 440, 500]:
plt.plot(inventory_paths[s]["日"][:90], inventory_paths[s]["期末在庫"][:90],
label=f"発注点{s}")
plt.title("発注点別の期末在庫推移（最初の90日）")
plt.xlabel("日")
plt.ylabel("期末在庫 [個]")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| 発注点 | 充足率 | 平均期末在庫 | 欠品数量 | 年間総費用 [円] | |
|---|---|---|---|---|---|
| 0 | 320 | 0.8914 | 150.2 | 4123 | 9339995 |
| 1 | 380 | 0.9824 | 179.9 | 668 | 3500500 |
| 2 | 440 | 0.9980 | 233.5 | 75 | 3117595 |
| 3 | 500 | 1.0000 | 289.6 | 0 | 3699325 |
結果の読み取り
発注点を上げるほど一般に充足率は改善しますが、平均在庫と保有費が増えます。表の最小費用案は、設定した欠品単価と保有単価の下での結論であり、絶対解ではありません。実務ではリードタイム自体の変動、最小発注量、複数品目の保管制約、バックオーダー、停止時需要を追加し、費用パラメータの感度も確認します。
No.094：故障シミュレーション — 予防保全周期を設計する
実務での意味
設備が劣化すると故障率が上がる場合、故障してから直す事後保全と、早めに止める予防保全のバランスが重要です。予防交換周期ごとの長期費用率を比較します。
分析・モデル化の考え方
寿命 をWeibull分布とし、形状母数 で経年劣化を表します。周期 より前に故障すれば故障修理費、故障しなければ で予防保全費が発生します。更新過程を多数サイクル生成し、
で評価します。各周期は同じ一様乱数から寿命を生成し、比較のばらつきを抑えます。
Pythonで確認する
shape, scale = 2.2, 170.0
n_cycles = 80_000
u_life = rng.random(n_cycles)
lifetimes = scale * (-np.log(1 - u_life)) ** (1 / shape)
maintenance_results = []
for tau in np.arange(60, 241, 20):
failed = lifetimes < tau
cycle_length = np.minimum(lifetimes, tau)
cycle_cost = np.where(failed, 1_200_000, 280_000)
maintenance_results.append({
"予防保全周期 [日]": tau,
"周期内故障確率": failed.mean(),
"費用率 [円/稼働日]": cycle_cost.sum() / cycle_length.sum(),
})
maintenance_table = pd.DataFrame(maintenance_results)
best_tau = maintenance_table.loc[maintenance_table["費用率 [円/稼働日]"].idxmin()]
display(maintenance_table.round(2))
print(f"最小費用の候補周期: {best_tau['予防保全周期 [日]']:.0f} 日")
plt.plot(maintenance_table["予防保全周期 [日]"],
maintenance_table["費用率 [円/稼働日]"], marker="o")
plt.axvline(best_tau["予防保全周期 [日]"], color="tab:red", ls="--", label="最小費用候補")
plt.title("予防保全周期と長期保全費用率")
plt.xlabel("予防保全周期 [日]")
plt.ylabel("費用率 [円/稼働日]")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| 予防保全周期 [日] | 周期内故障確率 | 費用率 [円/稼働日] | |
|---|---|---|---|
| 0 | 60 | 0.10 | 6365.89 |
| 1 | 80 | 0.17 | 5837.04 |
| 2 | 100 | 0.27 | 5788.45 |
| 3 | 120 | 0.37 | 5957.85 |
| 4 | 140 | 0.48 | 6222.70 |
| 5 | 160 | 0.58 | 6513.43 |
| 6 | 180 | 0.68 | 6807.53 |
| 7 | 200 | 0.76 | 7079.27 |
| 8 | 220 | 0.83 | 7317.14 |
| 9 | 240 | 0.88 | 7509.18 |
最小費用の候補周期: 100 日
結果の読み取り
周期を短くしすぎると予防保全費が増え、長くしすぎると高額な故障が増えるため、費用率は中間で最小になります。この周期はWeibull仮定と費用設定に依存します。故障による品質流出、安全、納期、連鎖停止を費用化できない場合は、費用最小だけでなく故障確率の上限制約を設けます。
No.095：エージェントシミュレーション — 製品ごとの動的振分けを比べる
実務での意味
並列設備へ製品を固定比率で振り分けると、処理時間の偶然差で片側だけが混雑することがあります。各製品を意思決定主体（エージェント）として、到着時に早く終わる設備を選ぶ方策を比較します。
分析・モデル化の考え方
製品エージェントは到着時刻と自身の処理時間を持ちます。固定方策では交互に設備を選び、動的方策では各設備の現在の完了予定時刻を見て、予想完了が早い側を選びます。局所ルールから工場全体のリードタイム分布が生じる点がエージェントモデルの特徴です。
Pythonで確認する
def simulate_routing(arrivals, service_a, service_b, policy):
available = np.zeros(2)
records = []
for i, arrival in enumerate(arrivals):
service = np.array([service_a[i], service_b[i]])
if policy == "固定交互":
machine = i % 2
else:
predicted_finish = np.maximum(arrival, available) + service
machine = int(np.argmin(predicted_finish))
start = max(arrival, available[machine])
finish = start + service[machine]
available[machine] = finish
records.append((i, arrival, machine, start, finish, finish - arrival))
return pd.DataFrame(records, columns=["製品", "到着", "設備", "開始", "完了", "リードタイム"])
n_agents = 1200
agent_arrivals = np.cumsum(rng.exponential(2.15, n_agents))
service_a = rng.lognormal(np.log(3.7), 0.28, n_agents)
service_b = rng.lognormal(np.log(4.1), 0.22, n_agents)
routing_tables = {
policy: simulate_routing(agent_arrivals, service_a, service_b, policy)
for policy in ["固定交互", "最短完了"]
}
routing_result = pd.DataFrame([
{"方策": policy, "平均リードタイム [分]": df["リードタイム"].mean(),
"95%点 [分]": df["リードタイム"].quantile(0.95),
"設備A配分率": np.mean(df["設備"] == 0)}
for policy, df in routing_tables.items()
])
display(routing_result.round(3))
for policy, df in routing_tables.items():
plt.hist(df["リードタイム"], bins=35, density=True, alpha=0.5, label=policy)
plt.title("振分け方策別の製品リードタイム分布")
plt.xlabel("リードタイム [分]")
plt.ylabel("確率密度")
plt.grid(True, axis="y", alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| 方策 | 平均リードタイム [分] | 95%点 [分] | 設備A配分率 | |
|---|---|---|---|---|
| 0 | 固定交互 | 55.83 | 141.602 | 0.500 |
| 1 | 最短完了 | 13.81 | 32.225 | 0.527 |
結果の読み取り
最短完了方策は設備ごとの混雑と速度差を利用し、固定交互よりリードタイムの平均や上側分位点を下げられる場合があります。ただし現場では、品種適合、治工具、段取り、作業者資格、搬送距離を無視できません。複雑なルールを追加するほど検証が難しくなるため、方策の説明可能性と障害時の手動運用も設計します。
No.096：マルコフ連鎖 — 設備状態から長期停止率を見積もる
実務での意味
設備を「健全・劣化・停止」の状態で日次管理すると、単純な故障件数より、劣化から停止へ進むリスクや修理後の復帰を表現できます。長期的な状態構成と、現在状態からの将来分布を求めます。
分析・モデル化の考え方
遷移行列 の要素 は、今日が状態 のとき翌日が状態 となる確率です。状態分布 は
で更新され、定常分布は を満たします。将来が現在状態だけに依存し、滞在時間を明示しないという仮定に注意します。
Pythonで確認する
state_names = np.array(["健全", "劣化", "停止"])
P = np.array([
[0.91, 0.08, 0.01],
[0.24, 0.66, 0.10],
[0.55, 0.15, 0.30],
])
eigenvalues, eigenvectors = np.linalg.eig(P.T)
stationary = np.real(eigenvectors[:, np.argmin(np.abs(eigenvalues - 1))])
stationary = stationary / stationary.sum()
prob = np.array([0.0, 1.0, 0.0]) # 現在は劣化
state_history = [prob.copy()]
for _ in range(30):
prob = prob @ P
state_history.append(prob.copy())
state_history = np.array(state_history)
display(pd.DataFrame(P, index=state_names, columns=state_names).round(2))
display(pd.DataFrame({"状態": state_names, "定常確率": stationary}).round(4))
for i, state in enumerate(state_names):
plt.plot(np.arange(31), state_history[:, i], label=state)
plt.axhline(stationary[2], color="tab:red", ls=":", label="定常停止率")
plt.title("現在『劣化』から始めた設備状態確率の推移")
plt.xlabel("経過日数 [日]")
plt.ylabel("状態確率")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| 健全 | 劣化 | 停止 | |
|---|---|---|---|
| 健全 | 0.91 | 0.08 | 0.01 |
| 劣化 | 0.24 | 0.66 | 0.10 |
| 停止 | 0.55 | 0.15 | 0.30 |
| 状態 | 定常確率 | |
|---|---|---|
| 0 | 健全 | 0.7640 |
| 1 | 劣化 | 0.1970 |
| 2 | 停止 | 0.0391 |
結果の読み取り
現在が劣化でも、日数が進むと初期状態の影響が薄れ、状態確率は定常分布へ近づきます。定常停止率は長期能力計画の入力になりますが、「劣化状態に何日いたか」で翌日の故障率が変わる設備には単純なマルコフ性が不十分です。その場合は状態を細分化するか、セミマルコフモデルや生存時間モデルを検討します。
No.097：離散イベントシミュレーション — 多工程ラインのボトルネックを特定する
実務での意味
プレス、加工、検査が直列につながると、各工程の平均時間だけでは工場全体のリードタイムを評価できません。製品の到着、処理開始、処理終了というイベントだけ時刻を進め、工程別待ち時間を測ります。
分析・モデル化の考え方
離散イベントシミュレーションでは、固定刻みで全状態を更新せず、次のイベント時刻へ時計を進めます。ここでは各工程の設備が空く時刻を管理し、製品ごとに直列3工程を通します。工程間バッファは無制限、故障なし、追越しなしという簡略化です。
Pythonで確認する
def simulate_serial_line(arrivals, service_matrix, machine_counts=(1, 1, 1)):
machine_available = [np.zeros(m) for m in machine_counts]
records = []
for job, arrival in enumerate(arrivals):
ready = arrival
for process in range(service_matrix.shape[1]):
machine = int(np.argmin(machine_available[process]))
start = max(ready, machine_available[process][machine])
finish = start + service_matrix[job, process]
records.append((job, process, machine, ready, start, finish, start - ready))
machine_available[process][machine] = finish
ready = finish
return pd.DataFrame(records, columns=["製品", "工程", "設備", "到着", "開始", "完了", "待ち"])
n_jobs = 900
line_arrivals = np.cumsum(rng.exponential(4.8, n_jobs))
service_matrix = np.column_stack([
rng.lognormal(np.log(3.2), 0.18, n_jobs),
rng.lognormal(np.log(4.3), 0.25, n_jobs),
rng.lognormal(np.log(3.8), 0.20, n_jobs),
])
line_result = simulate_serial_line(line_arrivals, service_matrix)
process_names = {0: "プレス", 1: "加工", 2: "検査"}
process_kpi = line_result.groupby("工程").agg(
平均待ち分=("待ち", "mean"),
待ち95点分=("待ち", lambda x: x.quantile(0.95)),
).reset_index()
process_kpi["工程名"] = process_kpi["工程"].map(process_names)
process_kpi["利用率概算"] = [
service_matrix[:, p].sum() / (line_result["完了"].max() - line_arrivals.min())
for p in range(3)
]
display(process_kpi[["工程名", "平均待ち分", "待ち95点分", "利用率概算"]].round(3))
plot_data = process_kpi.set_index("工程名")
plot_data["平均待ち分"].plot(kind="bar", color="tab:blue", alpha=0.75)
plt.title("直列生産ラインの工程別平均待ち時間")
plt.xlabel("工程")
plt.ylabel("平均待ち時間 [分]")
plt.grid(True, axis="y", alpha=0.3)
plt.xticks(rotation=0)
plt.tight_layout()
plt.show()| 工程名 | 平均待ち分 | 待ち95点分 | 利用率概算 | |
|---|---|---|---|---|
| 0 | プレス | 4.324 | 15.194 | 0.695 |
| 1 | 加工 | 39.328 | 92.109 | 0.962 |
| 2 | 検査 | 0.742 | 3.100 | 0.822 |
結果の読み取り
加工工程の利用率と待ち時間が大きければ、そこがライン全体のボトルネック候補です。ボトルネック以外の工程だけを高速化しても、全体リードタイムの改善は限定的です。実務モデルには有限バッファ、ブロッキング、段取り、故障、作業カレンダー、ロット搬送を加え、工程別実績と製品別リードタイムを別々に検証します。
No.098：ブラウン運動 — センサードリフトの閾値到達リスクを見る
実務での意味
計測器のゼロ点が小さな外乱を累積してずれると、平均が0でも時間とともに不確かさが広がります。校正許容幅へ到達する確率を見積もり、点検周期の候補を作ります。
分析・モデル化の考え方
標準ブラウン運動 は独立増分を持ち、 です。ドリフト量を とすると、分散は で時間に比例します。ここでは日次増分を生成し、複数経路と初回閾値到達を評価します。
Pythonで確認する
n_paths, n_days, sigma_sensor, limit = 4000, 60, 0.018, 0.12
increments = rng.normal(0, sigma_sensor, size=(n_paths, n_days))
brownian_paths = np.column_stack([np.zeros(n_paths), np.cumsum(increments, axis=1)])
crossed = np.any(np.abs(brownian_paths[:, 1:]) >= limit, axis=1)
first_crossing = np.where(
crossed, np.argmax(np.abs(brownian_paths[:, 1:]) >= limit, axis=1) + 1, np.nan
)
brownian_table = pd.DataFrame({
"KPI": ["60日以内閾値到達率", "到達した経路の初回到達日中央値", "60日目標準偏差"],
"値": [crossed.mean(), np.nanmedian(first_crossing), brownian_paths[:, -1].std(ddof=1)],
"単位": ["比率", "日", "mm"],
})
display(brownian_table.round(4))
days = np.arange(n_days + 1)
for path in brownian_paths[:20]:
plt.plot(days, path, alpha=0.35, lw=0.9)
plt.axhline(limit, color="tab:red", ls="--", label="校正許容幅")
plt.axhline(-limit, color="tab:red", ls="--")
plt.title("ブラウン運動で表したセンサードリフト経路")
plt.xlabel("経過日数 [日]")
plt.ylabel("ゼロ点ドリフト [mm]")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| KPI | 値 | 単位 | |
|---|---|---|---|
| 0 | 60日以内閾値到達率 | 0.7058 | 比率 |
| 1 | 到達した経路の初回到達日中央値 | 29.0000 | 日 |
| 2 | 60日目標準偏差 | 0.1406 | mm |
結果の読み取り
経路ごとの方向は予測できませんが、時間とともに分布が広がり、許容幅へ到達する確率を評価できます。この結果から点検周期を決める場合、見逃し損失と校正費を比較します。実際のセンサーには平均回帰、温度依存、段差変化、劣化傾向があり得るため、ブラウン運動は残差の独立性と分散成長が実績に合う場合に限って使います。
No.099：確率微分方程式 — 炉温の平均回帰と外乱を再現する
実務での意味
炉温は外乱で揺れますが、制御によって設定値へ戻ります。この「戻る力」とランダム外乱を同時に表し、温度規格外の確率や制御強化の効果を比較します。
分析・モデル化の考え方
Ornstein–Uhlenbeck過程を
と置きます。 は設定値、 は平均へ戻る速さ、 は外乱強度です。Euler–Maruyama法では
で離散化します。時間刻みを変えた安定性確認が必要です。
Pythonで確認する
def simulate_ou(n_paths, hours, dt, x0, mu, theta, sigma, random_generator):
n_steps = int(hours / dt)
x = np.empty((n_paths, n_steps + 1))
x[:, 0] = x0
z = random_generator.normal(size=(n_paths, n_steps))
for t in range(n_steps):
x[:, t + 1] = (x[:, t] + theta * (mu - x[:, t]) * dt
+ sigma * np.sqrt(dt) * z[:, t])
return x
dt, hours = 0.05, 8
ou_paths = simulate_ou(3000, hours, dt, x0=176, mu=180,
theta=1.1, sigma=2.0, random_generator=rng)
time_grid = np.arange(ou_paths.shape[1]) * dt
out_of_spec = np.mean((ou_paths < 176) | (ou_paths > 184))
ou_table = pd.DataFrame({
"KPI": ["8時間後平均温度", "8時間後標準偏差", "全時点の規格外割合"],
"値": [ou_paths[:, -1].mean(), ou_paths[:, -1].std(ddof=1), out_of_spec],
"単位": ["℃", "℃", "比率"],
})
display(ou_table.round(3))
for path in ou_paths[:15]:
plt.plot(time_grid, path, alpha=0.35, lw=0.9)
plt.axhline(180, color="black", lw=1.5, label="設定温度")
plt.axhline(184, color="tab:red", ls="--", label="管理範囲")
plt.axhline(176, color="tab:red", ls="--")
plt.title("平均回帰型SDEによる炉温シミュレーション")
plt.xlabel("経過時間 [時間]")
plt.ylabel("炉温 [℃]")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()| KPI | 値 | 単位 | |
|---|---|---|---|
| 0 | 8時間後平均温度 | 179.997 | ℃ |
| 1 | 8時間後標準偏差 | 1.372 | ℃ |
| 2 | 全時点の規格外割合 | 0.015 | 比率 |
結果の読み取り
初期温度が設定値より低くても、経路の平均は設定値へ近づきながら外乱で揺れます。規格外割合は制御強度と外乱強度の組合せで決まるため、制御改修案を比較できます。ただし観測ノイズとプロセス外乱を混同せず、推定期間を分けた再現性確認、残差診断、異常時のモデル外挙動を確認します。
No.100：デジタルツイン入門 — 実績で更新し、設備投資案を比較する
実務での意味
デジタルツインは3D表示そのものではなく、現場の状態・データ・モデルを継続的につなぎ、シナリオ比較とフィードバックに使う仕組みです。ここでは実績から到着・工程時間を更新し、現行ライン、加工高速化、加工設備増設を比較します。
分析・モデル化の考え方
簡易ツインを、(1)観測、(2)パラメータ更新、(3)シミュレーション、(4)KPI比較、(5)実績との差による再検証、のループとして構成します。入力は点推定だけでなく分布として保持します。同一の仮想需要と処理時間を施策間で共有し、差分を比較します。
Pythonで確認する
# 観測データからツインの入力を更新（説明用の単純な経験分布）
observed_interarrival = rng.exponential(4.9, 600)
observed_service = np.column_stack([
rng.lognormal(np.log(3.3), 0.18, 600),
rng.lognormal(np.log(4.4), 0.25, 600),
rng.lognormal(np.log(3.7), 0.20, 600),
])
n_twin_jobs = 1600
twin_arrivals = np.cumsum(rng.choice(observed_interarrival, n_twin_jobs, replace=True))
twin_service = np.column_stack([
rng.choice(observed_service[:, p], n_twin_jobs, replace=True) for p in range(3)
])
twin_scenarios = {
"現行": ((1, 1, 1), twin_service.copy()),
"加工10%高速化": ((1, 1, 1), twin_service * np.array([1.0, 0.9, 1.0])),
"加工設備を1台増設": ((1, 2, 1), twin_service.copy()),
}
twin_results = []
for scenario, (counts, services) in twin_scenarios.items():
result = simulate_serial_line(twin_arrivals, services, machine_counts=counts)
completion = result.groupby("製品")["完了"].max()
lead_time = completion.to_numpy() - twin_arrivals
makespan = completion.max() - twin_arrivals.min()
twin_results.append({
"シナリオ": scenario,
"平均リードタイム [分]": lead_time.mean(),
"95%点 [分]": np.quantile(lead_time, 0.95),
"スループット [個/8時間]": n_twin_jobs / makespan * 480,
})
twin_table = pd.DataFrame(twin_results)
twin_table["平均LT改善率"] = 1 - twin_table["平均リードタイム [分]"] / twin_table.loc[0, "平均リードタイム [分]"]
display(twin_table.round(3))
x = np.arange(len(twin_table))
plt.bar(x, twin_table["平均リードタイム [分]"], alpha=0.75)
plt.xticks(x, twin_table["シナリオ"], rotation=10)
plt.title("簡易デジタルツインによる改善シナリオ比較")
plt.xlabel("シナリオ")
plt.ylabel("平均リードタイム [分]")
plt.grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()| シナリオ | 平均リードタイム [分] | 95%点 [分] | スループット [個/8時間] | 平均LT改善率 | |
|---|---|---|---|---|---|
| 0 | 現行 | 39.323 | 78.132 | 98.509 | 0.000 |
| 1 | 加工10%高速化 | 23.699 | 44.927 | 98.517 | 0.397 |
| 2 | 加工設備を1台増設 | 18.906 | 34.792 | 98.509 | 0.519 |
結果の読み取り
改善案ごとに平均、95%点、スループットを比較すると、局所的な加工時間短縮と設備増設がライン全体へ与える効果の違いを確認できます。ただし、投資判断には設備費、保全要員、設置スペース、前後工程の制約を加える必要があります。ツインは一度作って終わりではなく、予測と実績の誤差を監視し、入力分布・ルール・版を承認付きで更新する運用が本体です。
対象ノックを通して見える実務上の示唆
- 平均値ではなく分布で比較する：平均待ち時間や期待費用に加え、95%点、SLA超過率、閾値到達率を示すと、繁忙時や異常側のリスクを判断できます。
- 比較案には同じ乱数条件を使う：共通の需要・寿命・処理時間を使うことで、偶然差ではなく施策差を読み取りやすくします。
- モデルの粒度は意思決定に合わせる：在庫方策なら日次、搬送や加工ならイベント単位、設備健全性なら状態遷移というように、必要な時間・主体・状態だけを表現します。
- ボトルネック以外の改善を過大評価しない：局所KPIの向上がライン全体のリードタイムやスループットへ伝わるかを必ず確認します。
- デジタルツインは更新プロセスまで含む：データ接続、入力推定、妥当性確認、シナリオ承認、結果監視が継続して初めて実務価値を持ちます。
シミュレーション結果は、前提条件付きの比較材料です。出力の桁数を増やすより、前提、誤差、適用範囲、判断ルールを明示する方が重要です。
実務導入する場合に必要なこと
- 意思決定とKPIの定義：誰が、いつ、何を選ぶためのモデルかを定め、費用・納期・品質・安全の優先順位を合意する
- 時刻と状態データの整備：到着、開始、終了、停止理由、設備、品種、ロット、作業カレンダーを一貫した定義で記録する
- 入力分析：分布適合、相関、時間帯、季節性、打切り、外れ値を確認し、点推定の過信を避ける
- VerificationとValidation：コードが仕様通りかを単体テストし、モデルが実績の平均・分位点・滞留箇所を再現するかを現場と確認する
- 感度・不確かさ分析：入力、費用、制約を変えたとき結論が反転しないかを確認する
- 運用ガバナンス：モデル版、データ期間、承認者、更新頻度、予実差、停止条件、手動判断への切替を管理する
まず限定ラインで現行運用と並走し、意思決定前に予測を固定して予実差を蓄積します。精度だけでなく、欠品削減、納期遵守、判断時間、現場負荷を含む業務KPIで評価します。
まとめ
No.091〜No.100では、期待損失のモンテカルロ積分から始め、待ち行列、在庫、故障、エージェント、マルコフ連鎖、離散イベント、ブラウン運動、確率微分方程式を実装し、最後に実績データで入力を更新する簡易デジタルツインへ統合しました。
共通する考え方は、現場の不確実性と時間依存を必要十分な粒度で表し、代替案を同じ条件で比較し、結果分布を業務KPIへ翻訳することです。実務導入では、高度なモデルより先に、イベント時刻、設備状態、品種・ロットのデータ定義と、予実差を継続監視する運用を整えることが重要です。
法人向けのご相談
数理工房では、製造業における生産ラインシミュレーション、在庫・保全方策の設計、設備投資シナリオ評価、デジタルツインPoC、Python研修、現場データ基盤から運用定着までをご支援しています。「平均能力では足りているのに納期が乱れる」「設備増設前に効果を比較したい」「シミュレーションを作ったが現場で更新できない」といった段階からご相談いただけます。
📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。