100本ノック / シミュレーション / シミュレーション100本ノック
製造業の連続シミュレーション入門｜在庫・需要・生産をPythonで可視化
生産・在庫の変化を先読みする：製造業の連続シミュレーション実践10題
本記事では、架空の精密部品工場を題材に、受注・生産・在庫・設備状態が時間とともにどう変化するかを連続シミュレーションで捉えます。No.061〜No.070を通じて、微分方程式の設計、数値計算、需要と在庫のフィードバック、計算の安定性までを一つの意思決定ストーリーとして確認します。
狙いは、数式を解くこと自体ではありません。「増産すると欠品はいつ解消するか」「発注ルールが在庫振動を起こさないか」「計算結果を会議の判断材料として信頼できるか」を、再現可能な形で検討できるようにすることです。
[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。
はじめに：この記事で扱う製造業の実務課題
対象工場は、1時間単位で連続運転する組立ラインを持ちます。営業は需要増を見込み、製造は能力増強を検討しています。一方で、仕掛品を増やしすぎればリードタイムと資金負担が膨らみ、完成品在庫を絞りすぎれば欠品が増えます。さらに、設備劣化や制御ルールの遅れが、計画にはない変動を生みます。
そこで「現時点の集計値」だけでなく、状態が次の状態を生む仕組みをモデル化します。連続シミュレーションは、実機で試しにくい増産、保全周期、在庫基準の変更を、時間軸上の仮想実験として比較する手段です。
現場でよくある状況
- 月平均では能力に余裕があるのに、需要の立ち上がり時に欠品する
- 在庫を見て生産量を調整した結果、増産と減産が交互に起きる
- 表計算の刻み幅や計算式が違うだけで、数日後の予測が変わる
- 部門ごとに「需要」「能力」「安全在庫」の前提が異なる
- シミュレーション結果はあるが、現場KPIや実行条件につながっていない
なぜこの問題は判断が難しいのか
製造システムには、蓄積、遅れ、非線形性、フィードバックがあります。在庫は一瞬で目標値にならず、流入と流出の差が時間をかけて積み上がります。設備劣化が進むと同じ指示量でも実生産量は下がり、需要変化に対する判断の遅れは過剰反応を招きます。このため、平均値の比較だけでは過渡的な欠品や振動を見落とします。
今回扱うノックの全体像
| No. | テーマ | 工場での問い | 主な確認指標 |
|---|---|---|---|
| 061 | 微分方程式モデル | 設備状態はどの速度で変わるか | 温度、平衡値 |
| 062 | オイラー法 | 時間刻みで結果はどれだけ変わるか | 近似誤差 |
| 063 | Runge-Kutta法 | 精度と計算量をどう両立するか | 終端誤差 |
| 064 | ロトカ・ヴォルテラモデル | 競合する資源が周期変動を生むか | 仕掛量、作業余力 |
| 065 | SIRモデル | 不良要因は工程間でどう波及するか | 影響中工程数、ピーク時点 |
| 066 | システムダイナミクス | 蓄積と流量をどう経営指標へ結ぶか | 仕掛品、スループット |
| 067 | 在庫ダイナミクス | 増産ルールで欠品を抑えられるか | 在庫、欠品時間 |
| 068 | 需要ダイナミクス | 需要ショックへの追随遅れは何を生むか | 需要、生産指示、累積不足 |
| 069 | フィードバックループ | 制御の強さは在庫を安定させるか | オーバーシュート、整定 |
| 070 | シミュレーションの安定性 | 数値上の振動を現象と誤認していないか | 安定条件、刻み幅感度 |
前半で数値計算の信頼性を押さえ、後半で在庫・需要・制御の意思決定へ展開します。
Python 環境の準備
外部データには依存せず、numpy、pandas、matplotlib、scipyを使います。乱数生成器の seed は固定します。グラフの日本語表示には japanize_matplotlib を利用します。
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import platform
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import japanize_matplotlib
from scipy.integrate import solve_ivp
SEED = 20260712
rng = np.random.default_rng(SEED)
pd.set_option("display.float_format", "{:.3f}".format)
print("Python :", platform.python_version())
print("numpy :", np.__version__)
print("pandas :", pd.__version__)
print("matplotlib :", matplotlib.__version__)
print("random seed:", SEED)Python : 3.11.9
numpy : 1.26.4
pandas : 2.2.2
matplotlib : 3.9.2
random seed: 20260712架空データの作成
14日間（336時間）の需要を想定します。基準需要は毎時100個、日中の周期性と週後半の販促増を加え、予測誤差に相当する小さな乱れを持たせます。値はすべて架空です。以降では、この需要系列と工場パラメータを共通の前提として使います。
hours = np.arange(14 * 24)
daily_cycle = 8 * np.sin(2 * np.pi * (hours - 6) / 24)
promotion = np.where(hours >= 7 * 24, 15.0, 0.0)
noise = rng.normal(0, 2.5, len(hours))
demand = np.clip(100 + daily_cycle + promotion + noise, 80, None)
factory_data = pd.DataFrame({
"時刻[h]": hours,
"需要[個/h]": demand,
"販促期間": np.where(hours >= 7 * 24, "販促後", "通常"),
})
display(factory_data.head())
display(factory_data.groupby("販促期間")["需要[個/h]"].agg(["mean", "std", "min", "max"]).round(2))
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(hours, demand, color="tab:blue", linewidth=1.3)
ax.axvline(7 * 24, color="tab:red", linestyle="--", label="販促開始")
ax.set_title("架空工場の時間別需要")
ax.set_xlabel("経過時間 [h]")
ax.set_ylabel("需要 [個/h]")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 時刻[h] | 需要[個/h] | 販促期間 | |
|---|---|---|---|
| 0 | 0 | 93.847 | 通常 |
| 1 | 1 | 93.977 | 通常 |
| 2 | 2 | 94.116 | 通常 |
| 3 | 3 | 95.370 | 通常 |
| 4 | 4 | 97.036 | 通常 |
| mean | std | min | max | |
|---|---|---|---|---|
| 販促期間 | ||||
| 販促後 | 114.780 | 6.200 | 101.920 | 129.560 |
| 通常 | 99.960 | 6.270 | 84.510 | 111.230 |
No.061：微分方程式モデル
実務での意味
設備温度、槽内濃度、摩耗量、仕掛品量などは、現在値だけでなく「単位時間あたりの変化」で管理すると、将来の上限到達時刻や定常状態を予測できます。ここでは設備温度を例に、加熱と放熱の釣り合いを確認します。
分析・モデル化の考え方
状態変数を設備温度 とし、一定の発熱 と外気への放熱を表す係数 を置きます。右辺が正なら温度は上昇し、ゼロになる温度が平衡点です。モデル境界、単位、初期値を明示することが実務利用の第一歩です。
Pythonで確認する
# 設備温度 T の一次遅れモデル: dT/dt = q - k(T - T_env)
T_env, heat_input, cooling = 25.0, 6.0, 0.18
t_eval = np.linspace(0, 24, 241)
def temperature_ode(t, y):
return [heat_input - cooling * (y[0] - T_env)]
sol_061 = solve_ivp(temperature_ode, [0, 24], [25.0], t_eval=t_eval)
T_eq = T_env + heat_input / cooling
summary_061 = pd.DataFrame({
"指標": ["初期温度", "24時間後温度", "理論平衡温度"],
"温度[℃]": [sol_061.y[0, 0], sol_061.y[0, -1], T_eq],
})
display(summary_061.round(2))
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(sol_061.t, sol_061.y[0], label="設備温度")
ax.axhline(T_eq, color="tab:red", linestyle="--", label="平衡温度")
ax.set_title("設備温度の連続時間モデル")
ax.set_xlabel("経過時間 [h]")
ax.set_ylabel("温度 [℃]")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 指標 | 温度[℃] | |
|---|---|---|
| 0 | 初期温度 | 25.000 |
| 1 | 24時間後温度 | 57.890 |
| 2 | 理論平衡温度 | 58.330 |
結果の読み取り
温度は初期の25℃から上昇し、理論平衡温度へ漸近します。24時間後の値だけでなく、上昇速度と上限が見えるため、警報温度、暖機時間、冷却能力の検討に使えます。ただし、実機適用前に発熱量と放熱係数を運転ログから同定する必要があります。
No.062：オイラー法
実務での意味
微分方程式を表計算や制御周期に落とすとき、最も直感的なのがオイラー法です。一方、刻み幅が粗いと安全限界への到達時刻を早く、または遅く見積もる危険があります。
分析・モデル化の考え方
時刻 の傾きを使って次の状態を直線的に予測します。刻み幅 を半分にして結果がほぼ変わらないかを確かめる「刻み幅感度分析」を、モデル検証の最低条件とします。
Pythonで確認する
def euler(f, y0, t):
y = np.empty(len(t), dtype=float)
y[0] = y0
for i in range(len(t) - 1):
dt = t[i + 1] - t[i]
y[i + 1] = y[i] + dt * f(t[i], y[i])
return y
f_temp = lambda t, T: heat_input - cooling * (T - T_env)
rows = []
fig, ax = plt.subplots(figsize=(8, 4))
for dt in [2.0, 1.0, 0.25]:
t = np.arange(0, 24 + dt, dt)
y = euler(f_temp, 25.0, t)
exact = T_eq + (25.0 - T_eq) * np.exp(-cooling * t)
rows.append({"刻み幅[h]": dt, "計算点数": len(t), "最大絶対誤差[℃]": np.max(np.abs(y - exact))})
ax.plot(t, y, marker="o", markersize=2, label=f"Euler dt={dt}h")
ax.plot(t_eval, T_eq + (25.0 - T_eq) * np.exp(-cooling * t_eval), color="black", linestyle="--", label="解析解")
display(pd.DataFrame(rows).round(4))
ax.set_title("オイラー法の刻み幅と近似精度")
ax.set_xlabel("経過時間 [h]")
ax.set_ylabel("温度 [℃]")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 刻み幅[h] | 計算点数 | 最大絶対誤差[℃] | |
|---|---|---|---|
| 0 | 2.000 | 13 | 2.582 |
| 1 | 1.000 | 25 | 1.194 |
| 2 | 0.250 | 97 | 0.281 |
結果の読み取り
刻み幅を小さくするほど解析解に近づきます。毎時データがあるからといって計算刻みも1時間で十分とは限りません。温度変化が速い起動直後や制御変更直後では細かくし、計算負荷との釣り合いを取ります。
No.063：Runge-Kutta法
実務での意味
長期予測や非線形モデルでは、オイラー法の誤差が累積します。設備投資や安全余裕の判断では、モデル構造の不確実性と数値計算の誤差を分ける必要があります。
分析・モデル化の考え方
4次Runge-Kutta法（RK4）は、1ステップ内の複数地点で傾きを評価し、重み付き平均で更新します。同じ刻み幅でも高精度になりやすい一方、傾き評価は4回必要です。精度だけでなく計算時間、説明可能性、再現性で手法を選びます。
Pythonで確認する
def rk4(f, y0, t):
y = np.empty(len(t), dtype=float)
y[0] = y0
for i in range(len(t) - 1):
dt = t[i + 1] - t[i]
k1 = f(t[i], y[i])
k2 = f(t[i] + dt/2, y[i] + dt*k1/2)
k3 = f(t[i] + dt/2, y[i] + dt*k2/2)
k4 = f(t[i] + dt, y[i] + dt*k3)
y[i + 1] = y[i] + dt * (k1 + 2*k2 + 2*k3 + k4) / 6
return y
t_1h = np.arange(0, 25, 1.0)
exact_1h = T_eq + (25.0 - T_eq) * np.exp(-cooling * t_1h)
euler_1h = euler(f_temp, 25.0, t_1h)
rk4_1h = rk4(f_temp, 25.0, t_1h)
comparison_063 = pd.DataFrame({
"手法": ["Euler", "RK4"],
"24時間後誤差[℃]": [abs(euler_1h[-1]-exact_1h[-1]), abs(rk4_1h[-1]-exact_1h[-1])],
"最大絶対誤差[℃]": [np.max(abs(euler_1h-exact_1h)), np.max(abs(rk4_1h-exact_1h))],
})
display(comparison_063.round(6))
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(t_1h, exact_1h, color="black", linewidth=2, label="解析解")
ax.plot(t_1h, euler_1h, "o--", label="Euler (1h)")
ax.plot(t_1h, rk4_1h, "s:", label="RK4 (1h)")
ax.set_title("オイラー法と4次Runge-Kutta法の比較")
ax.set_xlabel("経過時間 [h]")
ax.set_ylabel("温度 [℃]")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 手法 | 24時間後誤差[℃] | 最大絶対誤差[℃] | |
|---|---|---|---|
| 0 | Euler | 0.159 | 1.194 |
| 1 | RK4 | 0.000 | 0.000 |
結果の読み取り
1時間刻みでもRK4の誤差はオイラー法より大幅に小さくなります。高精度な手法は粗い入力データを正確にするものではありません。入力パラメータの推定誤差が支配的なら、数値精度を上げる前に計測設計を見直すべきです。
No.064：ロトカ・ヴォルテラモデル
実務での意味
仕掛品が増えると応援要員を投入し、仕掛品が減ると応援を解除する、といった相互作用は周期的な増減を生みます。古典的な捕食者・被食者モデルを、そのまま予測器ではなく「循環が起こる構造」の教材として使います。
分析・モデル化の考え方
仕掛品指数 と作業余力指数 が互いの増減率に影響すると仮定します。係数は因果の向きと反応速度を表します。実務では非負制約、能力上限、シフト切替を追加し、実測データで妥当性を検証します。
Pythonで確認する
# x: 処理待ち仕掛品、y: 応援可能な作業余力（概念モデル）
alpha, beta, delta, gamma = 0.55, 0.025, 0.012, 0.35
def lv_factory(t, z):
x, y = z
return [alpha*x - beta*x*y, delta*x*y - gamma*y]
t_lv = np.linspace(0, 80, 1601)
sol_064 = solve_ivp(lv_factory, [0, 80], [24, 15], t_eval=t_lv, rtol=1e-8, atol=1e-10)
lv_summary = pd.DataFrame({
"状態": ["仕掛品", "作業余力"],
"最小": sol_064.y.min(axis=1), "最大": sol_064.y.max(axis=1),
"平均": sol_064.y.mean(axis=1),
})
display(lv_summary.round(2))
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(t_lv, sol_064.y[0], label="仕掛品")
axes[0].plot(t_lv, sol_064.y[1], label="作業余力")
axes[0].set_title("仕掛品と作業余力の周期変動")
axes[0].set_xlabel("経過時間 [h]"); axes[0].set_ylabel("指数")
axes[0].grid(True, alpha=0.3); axes[0].legend()
axes[1].plot(sol_064.y[0], sol_064.y[1], color="tab:purple")
axes[1].set_title("状態空間上の循環")
axes[1].set_xlabel("仕掛品指数"); axes[1].set_ylabel("作業余力指数")
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()| 状態 | 最小 | 最大 | 平均 | |
|---|---|---|---|---|
| 0 | 仕掛品 | 17.120 | 45.850 | 29.940 |
| 1 | 作業余力 | 14.490 | 31.740 | 21.860 |
結果の読み取り
時系列と状態空間の両方で循環が確認できます。仕掛品ピークの後に作業余力が増える位相差は、遅れて応援する運用の典型です。周期が見えたからモデルが正しいとは限らず、応援発令・解除の実績時刻と照合して構造仮説を検証します。
No.065：SIRモデル
実務での意味
共通治具のずれ、誤った作業標準、材料ロットの異常などは複数工程へ波及します。SIRモデルを使うと、未影響・影響中・対策済みの工程数を分け、封じ込め速度の必要水準を議論できます。
分析・モデル化の考え方
全40工程を、未影響 、影響中 、対策済み に分類します。波及係数 は接触・共通要因の強さ、対策係数 は検知から封じ込めまでの速度に対応します。母数保存 も確認します。
Pythonで確認する
# 工程群における不良要因の波及をSIR型で表現
N = 40
beta_sir, gamma_sir = 0.42, 0.18
def sir_quality(t, z):
S, I, R = z
new_affected = beta_sir * S * I / N
contained = gamma_sir * I
return [-new_affected, new_affected-contained, contained]
t_sir = np.linspace(0, 40, 401)
sol_065 = solve_ivp(sir_quality, [0, 40], [39, 1, 0], t_eval=t_sir)
peak_i = np.argmax(sol_065.y[1])
display(pd.DataFrame({
"KPI": ["影響工程のピーク", "ピーク時点[日]", "40日後の未対策工程"],
"値": [sol_065.y[1, peak_i], sol_065.t[peak_i], sol_065.y[0, -1]],
}).round(2))
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(t_sir, sol_065.y[0], label="未影響 S")
ax.plot(t_sir, sol_065.y[1], label="影響中 I")
ax.plot(t_sir, sol_065.y[2], label="対策済み R")
ax.set_title("不良要因の工程間波及（SIR型モデル）")
ax.set_xlabel("経過日数 [日]")
ax.set_ylabel("工程数")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| KPI | 値 | |
|---|---|---|
| 0 | 影響工程のピーク | 8.800 |
| 1 | ピーク時点[日] | 16.000 |
| 2 | 40日後の未対策工程 | 5.500 |
結果の読み取り
影響中工程にはピークがあり、その時点までに検査能力や代替工程を確保する必要があります。 を高める施策は、初動連絡、トレーサビリティ、標準作業の即時更新です。独立でない工程を一様に扱う単純化があるため、実務では工程ネットワーク別の波及率へ拡張します。
No.066：システムダイナミクス
実務での意味
システムダイナミクスは、仕掛品や在庫のようなストックと、投入・完成のようなフローを明示します。KPIを個別最適で追うのではなく、因果関係と時間遅れを含む全体像で捉えるために有効です。
分析・モデル化の考え方
仕掛品 は投入率とスループットの差で変化します。ここでは、仕掛品が少ないと部材待ちで能力を使い切れず、多いほど最大能力へ近づく飽和関数を用います。パラメータは投入実績、完成実績、仕掛棚卸から推定できます。
Pythonで確認する
# 仕掛品 stock に、投入 inflow と完成 outflow が作用する
arrival_rate, max_capacity, half_saturation = 108.0, 125.0, 180.0
def wip_ode(t, z):
wip = z[0]
throughput = max_capacity * wip / (half_saturation + wip)
return [arrival_rate - throughput]
t_sd = np.linspace(0, 120, 481)
sol_066 = solve_ivp(wip_ode, [0, 120], [200], t_eval=t_sd)
throughput_066 = max_capacity * sol_066.y[0] / (half_saturation + sol_066.y[0])
display(pd.DataFrame({
"KPI": ["初期仕掛品", "120時間後仕掛品", "120時間後スループット"],
"値": [sol_066.y[0,0], sol_066.y[0,-1], throughput_066[-1]],
"単位": ["個", "個", "個/h"],
}).round(2))
fig, ax1 = plt.subplots(figsize=(8, 4))
ax1.plot(t_sd, sol_066.y[0], color="tab:blue", label="仕掛品")
ax1.set_xlabel("経過時間 [h]"); ax1.set_ylabel("仕掛品 [個]", color="tab:blue")
ax2 = ax1.twinx()
ax2.plot(t_sd, throughput_066, color="tab:orange", label="スループット")
ax2.set_ylabel("スループット [個/h]", color="tab:orange")
ax1.set_title("ストック・フローとして見た仕掛品")
ax1.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()| KPI | 値 | 単位 | |
|---|---|---|---|
| 0 | 初期仕掛品 | 200.000 | 個 |
| 1 | 120時間後仕掛品 | 1035.760 | 個 |
| 2 | 120時間後スループット | 106.490 | 個/h |
結果の読み取り
投入率が当初のスループットを上回るため仕掛品は増えますが、仕掛増に伴って完成率も上がり、やがて均衡へ向かいます。仕掛品を増やせば常に良いわけではなく、最大能力付近では追加仕掛の効果が小さく、滞留・品質・運転資金の負担が残ります。
No.067：在庫ダイナミクス
実務での意味
完成品在庫は需要変動を吸収しますが、増産には上限があります。在庫基準だけでなく、最低在庫、欠品継続時間、生産平準化を同時に見ることで、サービス水準と操業負荷のトレードオフを判断できます。
分析・モデル化の考え方
在庫 は生産 と需要 の差で更新します。生産は基準量に在庫偏差の補正を加え、設備能力の上下限で切ります。需要を満たせない場合、負在庫は受注残として解釈します。
Pythonで確認する
def simulate_inventory(gain, initial_inventory=900.0):
inv = np.empty(len(hours) + 1)
prod = np.empty(len(hours))
inv[0] = initial_inventory
target = 1000.0
for i, d in enumerate(demand):
prod[i] = np.clip(108 + gain * (target - inv[i]), 85, 130)
inv[i+1] = inv[i] + prod[i] - d
return inv, prod
inv_067, prod_067 = simulate_inventory(gain=0.04)
shortage_hours = int(np.sum(inv_067[1:] < 0))
display(pd.DataFrame({
"KPI": ["期末在庫", "最低在庫", "欠品時間", "平均生産量"],
"値": [inv_067[-1], inv_067.min(), shortage_hours, prod_067.mean()],
"単位": ["個", "個", "h", "個/h"],
}).round(2))
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(hours, inv_067[1:], label="完成品在庫")
ax.axhline(1000, color="tab:green", linestyle="--", label="目標在庫")
ax.axhline(0, color="tab:red", linewidth=1, label="欠品境界")
ax.set_title("需要変動下の完成品在庫ダイナミクス")
ax.set_xlabel("経過時間 [h]")
ax.set_ylabel("在庫 [個]")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| KPI | 値 | 単位 | |
|---|---|---|---|
| 0 | 期末在庫 | 843.920 | 個 |
| 1 | 最低在庫 | 796.080 | 個 |
| 2 | 欠品時間 | 0.000 | h |
| 3 | 平均生産量 | 107.200 | 個/h |
結果の読み取り
販促後の需要上昇に対して在庫が低下し、制御ルールが増産を促します。欠品時間がゼロでも最低在庫が小さければ予測誤差への余裕は不足しています。安全在庫、能力上限、残業コストを変えた複数シナリオで評価すべきです。
No.068：需要ダイナミクス
実務での意味
需要変化は即座に生産へ伝わりません。集計、承認、計画凍結、調達リードタイムの遅れが重なると、需要が増えた後もしばらく旧水準で生産し、累積不足が拡大します。
分析・モデル化の考え方
実需要 から認識需要 、さらに生産指示 へ、一次遅れを二段接続します。時定数 は反応の速さを表し、小さいほど追随が速くなります。累積需給差は必要な在庫バッファや受注残の規模に対応します。
Pythonで確認する
# 実需要に対し、指数平滑型の認識需要と生産指示が遅れて追随する
perceived = np.empty(len(hours))
order = np.empty(len(hours))
perceived[0], order[0] = demand[0], 100.0
tau_perception, tau_order = 18.0, 12.0
for i in range(1, len(hours)):
perceived[i] = perceived[i-1] + (demand[i-1] - perceived[i-1]) / tau_perception
order[i] = order[i-1] + (perceived[i-1] - order[i-1]) / tau_order
cumulative_gap = np.cumsum(demand - order)
display(pd.DataFrame({
"KPI": ["販促後平均需要", "販促後平均生産指示", "最大累積需給差"],
"値": [demand[168:].mean(), order[168:].mean(), cumulative_gap.max()],
"単位": ["個/h", "個/h", "個"],
}).round(2))
fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=True)
axes[0].plot(hours, demand, alpha=0.55, label="実需要")
axes[0].plot(hours, perceived, label="認識需要")
axes[0].plot(hours, order, label="生産指示")
axes[0].set_title("需要ショックと意思決定の遅れ")
axes[0].set_ylabel("数量 [個/h]"); axes[0].grid(True, alpha=0.3); axes[0].legend()
axes[1].plot(hours, cumulative_gap, color="tab:red")
axes[1].set_title("累積需給差（正は供給不足）")
axes[1].set_xlabel("経過時間 [h]"); axes[1].set_ylabel("累積差 [個]"); axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()| KPI | 値 | 単位 | |
|---|---|---|---|
| 0 | 販促後平均需要 | 114.780 | 個/h |
| 1 | 販促後平均生産指示 | 112.180 | 個/h |
| 2 | 最大累積需給差 | 598.950 | 個 |
結果の読み取り
販促開始後、認識需要と生産指示が段階的に遅れ、累積不足が増えます。対策は単純な増産だけでなく、販促情報の先行共有、凍結期間の短縮、部材手配の前倒しです。時定数は会議頻度ではなく、情報発生から実行までの実測リードタイムで設定します。
No.069：フィードバックループ
実務での意味
在庫偏差を見て生産を調整する負のフィードバックは、欠品を抑える基本です。しかし反応を強くしすぎると、生産量の変動、残業、段取り替えが増え、遅れがある場合は在庫振動を拡大します。
分析・モデル化の考え方
在庫補正係数 を政策変数として、弱・中・強の3ケースを比較します。評価は在庫の上下限だけでなく、生産量の標準偏差も含めます。実務では能力制約と意思決定遅れを入れた上で、頑健な係数範囲を選びます。
Pythonで確認する
feedback_rows = []
fig, ax = plt.subplots(figsize=(10, 4))
for gain in [0.01, 0.04, 0.12]:
inv, prod = simulate_inventory(gain=gain)
feedback_rows.append({
"フィードバック係数": gain,
"最低在庫[個]": inv.min(),
"最大在庫[個]": inv.max(),
"生産量標準偏差[個/h]": prod.std(),
"欠品時間[h]": np.sum(inv[1:] < 0),
})
ax.plot(hours, inv[1:], label=f"gain={gain}")
display(pd.DataFrame(feedback_rows).round(2))
ax.axhline(1000, color="black", linestyle="--", linewidth=1, label="目標在庫")
ax.set_title("フィードバック強度による在庫応答の違い")
ax.set_xlabel("経過時間 [h]")
ax.set_ylabel("在庫 [個]")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| フィードバック係数 | 最低在庫[個] | 最大在庫[個] | 生産量標準偏差[個/h] | 欠品時間[h] | |
|---|---|---|---|---|---|
| 0 | 0.010 | 542.260 | 1653.650 | 3.360 | 0 |
| 1 | 0.040 | 796.080 | 1248.940 | 6.410 | 0 |
| 2 | 0.120 | 900.000 | 1108.140 | 7.590 | 0 |
結果の読み取り
係数を大きくすると目標在庫への復帰は速くなりやすい一方、生産量の変動が増えます。最適な係数は単一KPIでは決まりません。欠品損失、残業費、段取り損失、在庫金利を共通の費用単位へ換算して比較する必要があります。
No.070：シミュレーションの安定性
実務での意味
グラフの振動や発散が、現場の不安定性ではなく計算方法のせいで生じることがあります。誤った設備投資や制御変更を避けるため、モデル妥当性とは別に数値安定性を検証します。
分析・モデル化の考え方
減衰系 を陽的オイラー法で解くと、更新ごとの増幅率は です。その絶対値が1未満なら誤差は減衰します。刻み幅を変えてKPIが収束することも確認します。
Pythonで確認する
# dy/dt = -lambda*y を陽的オイラー法で計算。安定条件は |1-lambda*dt| < 1
lam, y0, horizon = 1.2, 1.0, 10.0
stability_rows = []
fig, ax = plt.subplots(figsize=(8, 4))
for dt in [0.25, 1.0, 1.8]:
t = np.arange(0, horizon + dt, dt)
y = euler(lambda t, y: -lam*y, y0, t)
factor = 1 - lam*dt
stability_rows.append({
"刻み幅dt": dt, "増幅率(1-λdt)": factor,
"理論上安定": abs(factor) < 1, "最大絶対値": np.max(np.abs(y)),
})
ax.plot(t, y, marker="o", label=f"dt={dt}")
t_exact = np.linspace(0, horizon, 301)
ax.plot(t_exact, np.exp(-lam*t_exact), color="black", linestyle="--", label="解析解")
display(pd.DataFrame(stability_rows).round(3))
ax.set_title("時間刻みと数値安定性")
ax.set_xlabel("経過時間")
ax.set_ylabel("正規化した状態 y")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()| 刻み幅dt | 増幅率(1-λdt) | 理論上安定 | 最大絶対値 | |
|---|---|---|---|---|
| 0 | 0.250 | 0.700 | True | 1.000 |
| 1 | 1.000 | -0.200 | True | 1.000 |
| 2 | 1.800 | -1.160 | False | 2.436 |
結果の読み取り
細かい刻みでは滑らかに減衰し、中程度では符号を変えながら減衰します。安定条件を外れた粗い刻みでは、実際は減衰する系が発散して見えます。公開・稟議用モデルでは、採用ソルバー、許容誤差、刻み幅感度、再現環境を記録します。
対象ノックを通して見える実務上の示唆
- 状態と流量を分ける：在庫・仕掛品・設備状態はストックであり、受注・生産・劣化・回復はフローです。会議資料でも両者を混同しないことが重要です。
- 遅れをKPI化する：需要認識、承認、調達、増産の各遅れを時間で測ると、欠品の原因を「予測精度だけ」に帰さず改善できます。
- 平均ではなく軌跡を見る：平均在庫が適正でも、一時的な欠品や能力上限への張り付きがあれば運用リスクは残ります。
- 制御の強さには副作用がある：在庫回復を速める施策は、生産変動や現場負荷を増やす可能性があります。複数KPIで評価します。
- 数値誤差を分離する：刻み幅と解法を変えて結果が収束するか確認し、現象の不安定性と計算上の不安定性を区別します。
実務導入する場合に必要なこと
| 項目 | 最低限そろえるもの | 実務上の注意 |
|---|---|---|
| 目的 | 対象意思決定、KPI、比較シナリオ | 「現状再現」だけで終わらせない |
| データ | 時刻、在庫、投入、完成、需要、停止 | タイムスタンプと単位を統一する |
| モデル | 境界、状態、フロー、制約、遅れ | 現場レビューで因果の向きを確認する |
| 推定 | パラメータ根拠、期間、誤差 | 正常時と異常時を分ける |
| 検証 | ホールドアウト期間、感度分析、極端条件 | 数値安定性も別途確認する |
| 運用 | 更新頻度、責任者、判断基準、版管理 | 結果から実行までの導線を決める |
PoCでは、1つの製品群・1つの意思決定に範囲を絞り、過去期間の再現、反実仮想シナリオ、現場レビューを短い周期で回すと進めやすくなります。
まとめ
No.061〜No.070では、微分方程式を立てるところから、オイラー法とRK4、相互作用・波及・ストック＆フロー、在庫と需要のダイナミクス、フィードバック、数値安定性までを確認しました。連続シミュレーションの価値は、未来を一点で当てることではなく、前提と因果を共有し、変更案の結果とリスクを実機投入前に比較できることにあります。
分析結果を意思決定へつなげるには、モデル精度だけでなく、現場の制約、KPI間のトレードオフ、更新・承認プロセスまで設計することが欠かせません。
法人向けのご相談
数理工房では、製造業向けのシミュレーション設計、在庫・生産計画、設備投資のWhat-if分析、デジタルツイン構築、社内研修をご支援します。課題整理やデータ確認から、意思決定に使えるモデルのPoC、運用定着まで段階的に進めることができます。
📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。