100本ノック / 行列 / 行列100本ノック

工場シミュレーションをPythonで実践｜在庫・待ち行列・生産計画の意思決定

2020-03-11

変動に強い加工ラインを設計する

状態遷移・在庫・待ち行列・工場シミュレーション：100本ノック No.071〜No.080

製造現場では、設備劣化、需要変動、工程待ち、在庫切れが互いに影響します。本記事では、架空の精密部品ラインを題材に、状態を行列で表し、個別モデルを組み合わせ、改善案を工場全体のKPIで比較するまでを一つの意思決定ストーリーとして扱います。

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

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

架空の工場では、旋削・研削・検査の3工程で精密部品を生産しています。需要の波に合わせて投入量を増やすと仕掛品と待ち時間が増え、在庫を絞ると欠品が増えます。設備は稼働するほど劣化し、突発停止は後工程にも波及します。

そこで、状態遷移行列、在庫モデル、生産計画、待ち行列、マルコフモデル、離散事象シミュレーション、モンテカルロ法、デジタルツイン、最適化を段階的に使います。目的は精密な未来を言い当てることではなく、増員・在庫・保全・設備投資の選択肢を同じKPIで比較することです。

現場でよくある状況

設備ごとの稼働率は見えるが、劣化が翌週以降へどう蓄積するか分からない
安全在庫を経験で決め、欠品費用と保管費用の比較がない
月次生産計画と日々の待ち時間・故障停止が別々に管理されている
平均サイクルタイムだけで能力を判断し、ばらつきによる滞留を見落とす
改善案を単一の想定値で評価し、需要・故障・加工時間の不確実性を扱えていない
シミュレーション結果が現場実績と照合されず、意思決定に使われない

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

製造ラインは、ある時点の状態が次の状態へ影響する動的システムです。平均値が同じでもばらつきが大きい工程は待ち時間が長くなり、局所工程の高速化が工場全体の出荷増につながらないこともあります。さらに、需要、設備故障、修理時間は確率的です。

モデルを細かくすれば現実に近づくとは限りません。入力値の根拠、モデルの境界、KPIの定義、比較条件を明示し、単純な計算から始めて必要な複雑さだけを加えることが重要です。本記事では、同じ架空ラインを複数の解像度で表し、モデルごとの役割と限界を区別します。

今回扱うノックの全体像

No.	テーマ	製造業での判断
071	状態遷移行列	設備の劣化・故障状態が将来どう分布するか
072	在庫モデル	発注点と安全在庫をどう設定するか
073	生産計画	能力制約の中で何を何個作るか
074	待ち行列	検査工程の増員・能力増強が必要か
075	マルコフモデル	長期的な稼働・停止割合と保全効果を比較する
076	離散事象シミュレーション	到着・加工・故障を時系列イベントとして再現する
077	モンテカルロ法	不確実性を含む利益・納期達成確率を評価する
078	デジタルツイン	実績でモデルを校正し、差異を監視する
079	最適化	限られた改善予算をどこへ配分するか
080	工場シミュレーション	個別施策を統合し、全体KPIで比較する

Python 環境の準備

NumPyで行列・乱数計算、pandasで表、SciPyで最適化、Matplotlibで可視化を行います。外部データや専用シミュレーション製品には依存しません。乱数生成器は np.random.default_rng(71) で固定し、施策比較では同じ乱数系列を使えるよう関数ごとにseedを渡します。

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

import heapq
import math
import platform
import sys
from itertools import combinations

import japanize_matplotlib
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy
from IPython.display import display
from scipy.optimize import linprog

rng = np.random.default_rng(71)
pd.set_option("display.precision", 3)

print("Python:", sys.version.split()[0])
print("OS:", platform.platform())
print("NumPy:", np.__version__)
print("pandas:", pd.__version__)
print("SciPy:", scipy.__version__)
print("Matplotlib:", matplotlib.__version__)

Python: 3.13.1
OS: macOS-26.3-arm64-arm-64bit-Mach-O
NumPy: 2.5.1
pandas: 3.0.3
SciPy: 1.18.0
Matplotlib: 3.11.0

架空データの作成

対象は2製品（標準品A、高精度品B）を旋削・研削・検査するラインです。20営業日分の需要と実績サイクルタイムを生成します。需要には曜日変動、加工時間には右裾の長いばらつきを持たせます。

実務では、停止時間に段取り・休憩を含めるか、加工開始と完了のどちらを記録するか、不良による再加工をどう扱うかを先に定義します。ここでは各モデルの関係を見通せるよう、単位を「個」「分」「営業日」に統一します。

products = ["標準品A", "高精度品B"]
processes = ["旋削", "研削", "検査"]
days = np.arange(1, 21)

base_demand = np.array([42, 27])
weekday_factor = np.array([1.08, 0.96, 1.02, 1.12, 0.82])
demand = np.vstack([
    rng.poisson(base_demand * weekday_factor[(d - 1) % 5]) for d in days
])
demand_df = pd.DataFrame(demand, index=[f"Day {d:02d}" for d in days], columns=products)

# 行：工程、列：製品。標準サイクルタイム（分/個）
cycle_minutes = pd.DataFrame(
    [[6.0, 8.0], [7.5, 10.0], [4.0, 6.5]], index=processes, columns=products
)
observed_cycle = {
    p: rng.lognormal(np.log(cycle_minutes.loc[p].mean()) - 0.5 * 0.18**2, 0.18, 160)
    for p in processes
}

display(demand_df.head())
display(cycle_minutes)
fig, ax = plt.subplots(figsize=(8.2, 4.1))
ax.plot(days, demand[:, 0], marker="o", label=products[0])
ax.plot(days, demand[:, 1], marker="s", label=products[1])
ax.set_title("日別需要の推移（架空データ）")
ax.set_xlabel("営業日")
ax.set_ylabel("需要（個/日）")
ax.set_xticks(days[::2])
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	標準品A	高精度品B
Day 01	40	34
Day 02	34	26
Day 03	45	33
Day 04	43	39
Day 05	35	19

	標準品A	高精度品B
旋削	6.0	8.0
研削	7.5	10.0
検査	4.0	6.5

svg

No.071：状態遷移行列

実務での意味

設備状態を「正常・注意・故障」に分け、今日の状態から明日の状態へ移る確率を行列で表すと、点検対象台数や予備機の必要数を先読みできます。個別設備の故障日を断定するのではなく、設備群の状態構成を予測する用途に向きます。

分析・モデル化の考え方

状態分布を行ベクトル $\boldsymbol{p}_t$ 、遷移行列を $P$ とすると、 $k$ 日後は

\boldsymbol{p}_{t+k}=\boldsymbol{p}_tP^k, \qquad P_{ij}=\Pr(S_{t+1}=j\mid S_t=i), \qquad \sum_jP_{ij}=1

で求めます。行和が1であること、確率が非負であることを確認します。実績から推定するときは、設備型式、負荷帯、保全方針が変わった期間を無条件に混ぜないことが重要です。

Pythonで確認する

states = ["正常", "注意", "故障"]
P = np.array([
    [0.90, 0.09, 0.01],
    [0.30, 0.58, 0.12],
    [0.65, 0.00, 0.35],
])
p0 = np.array([0.80, 0.15, 0.05])
state_path = np.vstack([p0 @ np.linalg.matrix_power(P, k) for k in range(15)])
state_df = pd.DataFrame(state_path, index=np.arange(15), columns=states)

display(pd.DataFrame(P, index=[f"現在:{s}" for s in states], columns=[f"翌日:{s}" for s in states]))
display(state_df.iloc[[0, 1, 3, 7, 14]].style.format("{:.1%}"))
assert np.allclose(P.sum(axis=1), 1)

fig, ax = plt.subplots(figsize=(7.6, 4.0))
for state in states:
    ax.plot(state_df.index, state_df[state] * 100, marker="o", label=state)
ax.set_title("状態遷移行列による設備状態分布の予測")
ax.set_xlabel("経過日数")
ax.set_ylabel("設備構成比（%）")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	翌日:正常	翌日:注意	翌日:故障
現在:正常	0.90	0.09	0.01
現在:注意	0.30	0.58	0.12
現在:故障	0.65	0.00	0.35

	正常	注意	故障
0	80.0%	15.0%	5.0%
1	79.8%	15.9%	4.3%
3	79.1%	16.7%	4.2%
7	78.8%	16.9%	4.3%
14	78.8%	16.9%	4.3%

svg

結果の読み取り

初日の構成比から、注意・故障状態が時間とともに一定の比率へ近づきます。例えば設備20台なら、予測故障比率に20を掛けて修理対応台数の期待値を見積もれます。ただし期待値は「必ずその台数になる」という意味ではありません。日別の遷移回数、型式別の差、修理後の初期故障を確認し、遷移確率には推定誤差を付けて運用します。

No.072：在庫モデル

実務での意味

在庫は欠品を防ぐ一方、保管費、陳腐化、資金拘束を生みます。発注点方式では、在庫ポジションが閾値を下回ったときに補充し、需要変動と調達リードタイムを吸収します。

分析・モデル化の考え方

日需要の平均を $\mu_d$ 、標準偏差を $\sigma_d$ 、リードタイムを $L$ 日、安全係数を $z$ とすると、独立同分布を仮定した発注点は

R=\mu_dL+z\sigma_d\sqrt{L}

です。前半はリードタイム中の平均需要、後半は安全在庫です。実務では曜日性、需要相関、リードタイム自体の変動、発注ロット、未納残を含めます。

Pythonで確認する

daily_a = demand_df["標準品A"].to_numpy()
mu_d, sigma_d = daily_a.mean(), daily_a.std(ddof=1)
lead_time, z = 3, 1.645  # 片側95%の目安
reorder_point = int(np.ceil(mu_d * lead_time + z * sigma_d * np.sqrt(lead_time)))
order_qty = int(np.ceil(mu_d * 5))

def inventory_run(daily, initial, reorder, quantity, lead=3):
    on_hand, pipeline, rows = initial, [], []
    for day, req in enumerate(daily, 1):
        arrivals = sum(q for arrival, q in pipeline if arrival == day)
        pipeline = [(arrival, q) for arrival, q in pipeline if arrival > day]
        on_hand += arrivals
        shipped = min(on_hand, req)
        shortage = req - shipped
        on_hand -= shipped
        position = on_hand + sum(q for _, q in pipeline)
        ordered = quantity if position <= reorder else 0
        if ordered:
            pipeline.append((day + lead, ordered))
        rows.append([day, req, arrivals, shipped, shortage, on_hand, position, ordered])
    return pd.DataFrame(rows, columns=["日", "需要", "入荷", "出荷", "欠品", "期末在庫", "在庫ポジション", "発注"])

inventory_df = inventory_run(daily_a, reorder_point + order_qty, reorder_point, order_qty, lead_time)
inventory_kpi = pd.DataFrame({
    "KPI": ["発注点", "発注量", "平均期末在庫", "充足率"],
    "値": [reorder_point, order_qty, inventory_df["期末在庫"].mean(), inventory_df["出荷"].sum() / inventory_df["需要"].sum()],
})
display(inventory_kpi.style.format({"値": "{:.2f}"}))

fig, ax = plt.subplots(figsize=(8.0, 4.0))
ax.step(inventory_df["日"], inventory_df["期末在庫"], where="mid", label="期末在庫")
ax.axhline(reorder_point, color="tab:red", linestyle="--", label=f"発注点 {reorder_point}")
ax.bar(inventory_df["日"], inventory_df["入荷"], alpha=0.25, label="入荷")
ax.set_title("発注点方式による標準品Aの在庫推移")
ax.set_xlabel("営業日")
ax.set_ylabel("数量（個）")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	KPI	値
0	発注点	144.00
1	発注量	201.00
2	平均期末在庫	148.70
3	充足率	1.00

svg

結果の読み取り

発注点は「平均需要×3日」に安全在庫を加えた値です。この架空期間では高い充足率を維持できますが、20日間だけの結果で95%のサービス水準を保証することはできません。欠品費、保管費、最小発注量を含め、複数年相当の需要系列で発注点と発注量を比較します。需要急増時には予測連動、供給途絶時には別の例外ルールが必要です。

No.073：生産計画

実務での意味

受注をすべて作れないとき、売上だけでなく限界利益、納期、重要顧客、後工程の能力を見て製品構成を決めます。線形計画法は「何がボトルネックで、能力を増やす価値がどこにあるか」を明示します。

分析・モデル化の考え方

製品 $j$ の生産量を $x_j$ 、単位利益を $c_j$ 、工程 $i$ の必要時間を $a_{ij}$ 、能力を $b_i$ とすると、

\max_{\boldsymbol{x}}\ \boldsymbol{c}^{\mathsf T}\boldsymbol{x} \quad \text{s.t.}\quad A\boldsymbol{x}\leq\boldsymbol{b}, \quad 0\leq x_j\leq d_j

です。連続解は能力配分の基準になりますが、実際のロット、段取り順、最低生産量は整数・混合整数モデルで扱います。

Pythonで確認する

profit = np.array([3200, 5100])
capacity = np.array([780, 850, 520])  # 分/日
requirements = cycle_minutes.to_numpy()
upper_demand = np.array([70, 50])
plan = linprog(-profit, A_ub=requirements, b_ub=capacity, bounds=list(zip([0, 0], upper_demand)), method="highs")
production_plan = plan.x
used_minutes = requirements @ production_plan

plan_df = pd.DataFrame({
    "製品": products,
    "計画数量": production_plan,
    "需要上限": upper_demand,
    "単位利益（円）": profit,
})
capacity_df = pd.DataFrame({
    "工程": processes, "使用時間（分）": used_minutes,
    "能力（分）": capacity, "負荷率": used_minutes / capacity,
})
display(plan_df.style.format({"計画数量": "{:.1f}", "単位利益（円）": "{:,.0f}"}))
display(capacity_df.style.format({"使用時間（分）": "{:.1f}", "負荷率": "{:.1%}"}))
print("日次限界利益:", f"{profit @ production_plan:,.0f} 円")

fig, ax = plt.subplots(figsize=(7.4, 4.0))
ax.bar(capacity_df["工程"], capacity_df["負荷率"] * 100, color=["steelblue", "tab:orange", "tab:green"])
ax.axhline(100, color="tab:red", linestyle="--", label="能力上限")
ax.set_title("最適生産計画における工程別負荷率")
ax.set_xlabel("工程")
ax.set_ylabel("負荷率（%）")
ax.grid(True, axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	製品	計画数量	需要上限	単位利益（円）
0	標準品A	46.7	70	3,200
1	高精度品B	50.0	50	5,100

	工程	使用時間（分）	能力（分）	負荷率
0	旋削	680.0	780	87.2%
1	研削	850.0	850	100.0%
2	検査	511.7	520	98.4%

日次限界利益: 404,333 円

svg

結果の読み取り

負荷率が100%に達する工程が計画上のボトルネックです。その工程の能力追加、サイクル短縮、外注化を検討する根拠になります。ただし、線形モデルは段取り損失、故障、日内の到着順を平均化しています。最適値をそのまま現場指示にせず、整数ロットへ丸めた後に能力制約を再確認し、次の待ち行列・シミュレーションで実現可能性を検証します。

No.074：待ち行列

実務での意味

検査工程は平均能力に余裕があっても、到着と処理時間のばらつきで待ちが発生します。待ち行列モデルを使うと、検査員を1名から2名へ増やした場合の待ち時間と滞留数を概算できます。

分析・モデル化の考え方

到着率 $\lambda$ 、1人当たり処理率 $\mu$ 、窓口数 $c$ の M/M/c を考えます。利用率は

\rho=\frac{\lambda}{c\mu}<1

が安定条件です。Erlang C式から待ち確率と平均待ち時間 $W_q$ を求めます。指数分布仮定が合わない場合でも、能力余裕と待ちの非線形な関係を理解する基準モデルになります。

Pythonで確認する

def mmc_metrics(arrival_rate, service_rate, servers):
    offered = arrival_rate / service_rate
    rho = offered / servers
    if rho >= 1:
        return {"窓口数": servers, "利用率": rho, "待ち確率": 1.0, "平均待ち時間_分": np.inf, "系内平均個数": np.inf}
    terms = sum(offered**n / math.factorial(n) for n in range(servers))
    tail = offered**servers / (math.factorial(servers) * (1 - rho))
    p0 = 1 / (terms + tail)
    p_wait = tail * p0
    wq = p_wait / (servers * service_rate - arrival_rate)
    return {"窓口数": servers, "利用率": rho, "待ち確率": p_wait, "平均待ち時間_分": wq * 60, "系内平均個数": arrival_rate * (wq + 1 / service_rate)}

arrival_rate, service_rate = 8.0, 9.0  # 個/時
queue_df = pd.DataFrame([mmc_metrics(arrival_rate, service_rate, c) for c in [1, 2, 3]])
display(queue_df.style.format({"利用率": "{:.1%}", "待ち確率": "{:.1%}", "平均待ち時間_分": "{:.1f}", "系内平均個数": "{:.2f}"}))

fig, ax = plt.subplots(figsize=(7.2, 4.0))
ax.bar(queue_df["窓口数"].astype(str), queue_df["平均待ち時間_分"], color="tab:orange")
ax.set_title("検査窓口数と理論平均待ち時間（M/M/c）")
ax.set_xlabel("検査窓口数")
ax.set_ylabel("平均待ち時間（分）")
ax.grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()

	窓口数	利用率	待ち確率	平均待ち時間_分	系内平均個数
0	1	88.9%	88.9%	53.3	8.00
1	2	44.4%	27.4%	1.6	1.11
2	3	29.6%	6.8%	0.2	0.92

svg

結果の読み取り

窓口1つでは利用率が高く、わずかな変動でも待ちが急増します。2つに増やすと平均待ち時間は大きく下がりますが、3つ目の改善幅は相対的に小さくなります。実務では平均だけでなく95パーセンタイル、ピーク時間帯、検査員の兼務、製品別検査時間を確認します。増員費用と納期遅延・仕掛在庫の削減額を比べて判断します。

No.075：マルコフモデル

実務での意味

状態遷移行列を長期運用へ広げると、設備が正常・注意・故障にいる時間割合を見積もれます。予防保全によって「注意→正常」を増やし「注意→故障」を減らした場合、長期稼働率がどれだけ変わるかを比較できます。

分析・モデル化の考え方

定常分布 $\boldsymbol{\pi}$ は

\boldsymbol{\pi}=\boldsymbol{\pi}P, \qquad \sum_i\pi_i=1

を満たします。固有値1に対応する左固有ベクトル、または反復計算で求めます。有限期間の投資評価では初期状態からの過渡変化も重要であり、定常分布だけで判断しません。

Pythonで確認する

P_preventive = np.array([
    [0.92, 0.075, 0.005],
    [0.48, 0.47, 0.05],
    [0.72, 0.00, 0.28],
])

def stationary_distribution(matrix):
    values, vectors = np.linalg.eig(matrix.T)
    vector = np.real(vectors[:, np.argmin(np.abs(values - 1))])
    return vector / vector.sum()

pi_base = stationary_distribution(P)
pi_pm = stationary_distribution(P_preventive)
markov_compare = pd.DataFrame({"状態": states, "現状": pi_base, "予防保全後": pi_pm})
markov_compare["差（pt）"] = (markov_compare["予防保全後"] - markov_compare["現状"]) * 100
display(markov_compare.style.format({"現状": "{:.1%}", "予防保全後": "{:.1%}", "差（pt）": "{:+.1f}"}))

fig, ax = plt.subplots(figsize=(7.4, 4.0))
x = np.arange(len(states))
ax.bar(x - 0.18, pi_base * 100, width=0.36, label="現状")
ax.bar(x + 0.18, pi_pm * 100, width=0.36, label="予防保全後")
ax.set_title("マルコフモデルによる長期状態構成の比較")
ax.set_xlabel("設備状態")
ax.set_ylabel("定常構成比（%）")
ax.set_xticks(x, states)
ax.grid(True, axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	状態	現状	予防保全後	差（pt）
0	正常	78.8%	86.3%	+7.5
1	注意	16.9%	12.2%	-4.7
2	故障	4.3%	1.4%	-2.9

svg

結果の読み取り

予防保全シナリオでは正常状態の長期比率が上がり、故障状態の比率が下がります。故障比率の差に設備台数、1日当たり停止損失、年間日数を掛ければ便益の粗い上限を作れます。ただし保全時間中の停止、部品費、保全による初期不良、遷移確率の信頼区間を含めて投資採算を評価します。

No.076：離散事象シミュレーション

実務での意味

離散事象シミュレーションは、製品の到着、加工開始、加工完了を時刻順に処理します。平均値だけでは見えない待ち時間の分布、ピーク時の滞留、先着順ルールの影響を再現できます。

分析・モデル化の考え方

将来イベントを優先度付きキューへ入れ、最も早いイベントへ時計を進めます。単一工程なら、ジョブ $i$ の開始・完了時刻は

B_i=\max(A_i,C_{i-1}),\qquad C_i=B_i+S_i, \qquad W_i=B_i-A_i

です。ここではこの再帰をイベントログとして実装します。イベント駆動にすると、故障・修理・複数工程も同じ考え方で追加できます。

Pythonで確認する

def single_station_des(n_jobs=120, arrival_mean=7.5, service_mean=6.5, seed=76):
    local_rng = np.random.default_rng(seed)
    arrivals = np.cumsum(local_rng.exponential(arrival_mean, n_jobs))
    services = local_rng.lognormal(np.log(service_mean) - 0.5 * 0.25**2, 0.25, n_jobs)
    available = 0.0
    records = []
    event_queue = []
    for job, (arrival, service) in enumerate(zip(arrivals, services), 1):
        heapq.heappush(event_queue, (arrival, "到着", job))
        start = max(arrival, available)
        finish = start + service
        available = finish
        heapq.heappush(event_queue, (finish, "完了", job))
        records.append([job, arrival, start, finish, start - arrival, service])
    return pd.DataFrame(records, columns=["ジョブ", "到着時刻", "開始時刻", "完了時刻", "待ち時間", "加工時間"]), event_queue

des_df, event_queue = single_station_des()
display(des_df.head().style.format({c: "{:.1f}" for c in des_df.columns[1:]}))
display(des_df[["待ち時間", "加工時間"]].describe(percentiles=[0.5, 0.9, 0.95]).round(2))

fig, ax = plt.subplots(figsize=(7.6, 4.0))
ax.hist(des_df["待ち時間"], bins=18, color="steelblue", edgecolor="white")
ax.axvline(des_df["待ち時間"].quantile(0.95), color="tab:red", linestyle="--", label="95%点")
ax.set_title("離散事象シミュレーションによる待ち時間分布")
ax.set_xlabel("待ち時間（分）")
ax.set_ylabel("ジョブ数")
ax.grid(True, axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	ジョブ	到着時刻	開始時刻	完了時刻	待ち時間	加工時間
0	1	13.8	13.8	21.1	0.0	7.3
1	2	28.5	28.5	38.6	0.0	10.1
2	3	34.7	38.6	44.2	3.9	5.6
3	4	43.9	44.2	51.0	0.4	6.8
4	5	44.1	51.0	57.6	6.9	6.6

	待ち時間	加工時間
count	120.00	120.00
mean	17.25	6.46
std	15.44	1.57
min	0.00	3.52
50%	13.54	6.50
90%	42.94	8.33
95%	48.76	9.81
max	58.99	11.55

svg

結果の読み取り

平均待ち時間だけでなく、長く待つジョブの裾が確認できます。納期やバッファ設計には平均より90%点・95%点が有用です。この例は単一設備で、故障や休憩を省略しています。本番モデルでは実績ログから到着間隔・加工時間の分布を検証し、優先品、段取り、再加工、シフト境界を必要な範囲で追加します。

No.077：モンテカルロ法

実務での意味

事業計画では需要、歩留まり、停止時間を1つの想定値に固定しがちです。モンテカルロ法は複数の不確実性を繰り返しサンプリングし、利益の期待値だけでなく赤字確率や下振れ幅を示します。

分析・モデル化の考え方

不確実な入力 $\boldsymbol{X}$ からKPI $Y=g(\boldsymbol{X})$ を計算し、 $N$ 回の標本 $Y^{(1)},\ldots,Y^{(N)}$ で分布を近似します。

\widehat{E[Y]}=\frac{1}{N}\sum_{r=1}^NY^{(r)}, \qquad \widehat{\Pr}(Y<y_0)=\frac{1}{N}\sum_{r=1}^N\mathbf{1}(Y^{(r)}<y_0)

反復回数を増やしても入力分布の誤りは直りません。分布形、変数間相関、極端事象を現場データと知見から定義します。

Pythonで確認する

n_sim = 10_000
mc_rng = np.random.default_rng(77)
monthly_demand = np.maximum(0, mc_rng.normal(1450, 180, n_sim))
yield_rate = np.clip(mc_rng.beta(90, 6, n_sim), 0, 1)
downtime_hours = mc_rng.gamma(shape=2.2, scale=7.0, size=n_sim)
capacity_units = np.maximum(0, 1600 - downtime_hours * 5.5)
good_units = np.minimum(monthly_demand, capacity_units) * yield_rate
sales = good_units * 8_200
variable_cost = np.minimum(monthly_demand, capacity_units) * 4_600
fixed_cost = 4_100_000
profit_mc = sales - variable_cost - fixed_cost

mc_kpi = pd.DataFrame({
    "KPI": ["平均利益", "利益5%点", "利益95%点", "赤字確率"],
    "値": [profit_mc.mean(), np.quantile(profit_mc, 0.05), np.quantile(profit_mc, 0.95), np.mean(profit_mc < 0)],
})
display(mc_kpi.style.format({"値": lambda x: f"{x:,.0f}" if abs(x) > 1 else f"{x:.1%}"}))

fig, ax = plt.subplots(figsize=(7.6, 4.0))
ax.hist(profit_mc / 1e6, bins=35, color="tab:green", edgecolor="white")
ax.axvline(0, color="tab:red", linestyle="--", label="損益分岐")
ax.set_title("需要・歩留まり・停止を含む月次利益分布")
ax.set_xlabel("月次利益（百万円）")
ax.set_ylabel("試行回数")
ax.grid(True, axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	KPI	値
0	平均利益	233,893
1	利益5%点	-657,416
2	利益95%点	948,985
3	赤字確率	29.6%

svg

結果の読み取り

平均利益がプラスでも、下側5%点や赤字確率を見ると意思決定の安全余裕が分かります。設備投資案は平均利益の増加だけでなく、赤字確率や納期未達確率をどこまで下げるかで比較します。実務では需要と価格、停止と歩留まりの相関を検討し、分布外の長期停止はストレスシナリオとして別に評価します。

No.078：デジタルツイン

実務での意味

デジタルツインは、現場の状態を継続的に取り込み、仮想モデルとの差を更新しながら施策を試す仕組みです。3D表示そのものではなく、現実との同期、校正、差異監視、判断への接続が中心です。

分析・モデル化の考え方

工程 $i$ のモデル平均サイクル $m_i$ と実績平均 $\bar{x}_i$ の校正係数を

k_i=\frac{\bar{x}_i}{m_i},\qquad m_i^{\mathrm{updated}}=k_im_i

とします。単純な比率校正でも、どの工程で仮想モデルが楽観的かを特定できます。校正データと妥当性確認データを分け、更新後に別期間を再現できるか確認します。

Pythonで確認する

model_cycle = cycle_minutes.mean(axis=1).to_numpy()
actual_cycle = np.array([observed_cycle[p].mean() for p in processes])
calibration = actual_cycle / model_cycle
updated_cycle = model_cycle * calibration
twin_df = pd.DataFrame({
    "工程": processes,
    "モデル初期値（分）": model_cycle,
    "実績平均（分）": actual_cycle,
    "校正係数": calibration,
    "更新後（分）": updated_cycle,
    "初期差異（%）": (model_cycle / actual_cycle - 1) * 100,
})
display(twin_df.style.format({c: "{:.2f}" for c in twin_df.columns[1:]}))

fig, ax = plt.subplots(figsize=(7.8, 4.1))
x = np.arange(len(processes))
ax.bar(x - 0.26, model_cycle, width=0.26, label="モデル初期値")
ax.bar(x, actual_cycle, width=0.26, label="現場実績")
ax.bar(x + 0.26, updated_cycle, width=0.26, label="校正後")
ax.set_title("工程別サイクルタイムのデジタルツイン校正")
ax.set_xlabel("工程")
ax.set_ylabel("平均サイクルタイム（分/個）")
ax.set_xticks(x, processes)
ax.grid(True, axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	工程	モデル初期値（分）	実績平均（分）	校正係数	更新後（分）	初期差異（%）
0	旋削	7.00	6.92	0.99	6.92	1.14
1	研削	8.75	8.68	0.99	8.68	0.83
2	検査	5.25	5.29	1.01	5.29	-0.78

svg

結果の読み取り

校正係数が1より大きい工程では、初期モデルが現場を楽観的に表しています。更新後に平均は一致しますが、これは妥当性確認の完了を意味しません。分散、時間帯、製品構成、故障頻度について別期間で誤差を確認します。設備ID・時刻・単位の整合、センサー欠測、モデル版、更新承認者を管理して初めて運用可能なツインになります。

No.079：最適化

実務での意味

改善候補が複数あっても、予算と要員は限られます。各施策の費用と期待削減損失を整理し、同時実施できない組合せも含めて、予算内で最大効果となるポートフォリオを選びます。

分析・モデル化の考え方

施策 $j$ を採用するかを二値変数 $x_j\in\{0,1\}$ とし、便益 $v_j$ 、費用 $c_j$ 、予算 $B$ に対して

\max \sum_jv_jx_j \quad\text{s.t.}\quad \sum_jc_jx_j\leq B

を解きます。便益はシミュレーションや実績から推定し、相乗効果・排他条件・実行要員も制約にします。ここでは候補が少ないため全組合せを列挙し、透明に比較します。

Pythonで確認する

actions = pd.DataFrame({
    "施策": ["検査応援", "研削予防保全", "安全在庫追加", "段取り短縮", "センサー増設"],
    "費用_万円": [90, 130, 70, 160, 110],
    "期待年間便益_万円": [170, 240, 105, 260, 150],
})
budget = 300
portfolios = []
for mask in range(1 << len(actions)):
    selected = [i for i in range(len(actions)) if mask & (1 << i)]
    cost = actions.loc[selected, "費用_万円"].sum()
    benefit = actions.loc[selected, "期待年間便益_万円"].sum()
    feasible = cost <= budget
    portfolios.append({"施策": "・".join(actions.loc[selected, "施策"]) or "実施なし", "費用_万円": cost, "便益_万円": benefit, "予算内": feasible})
portfolio_df = pd.DataFrame(portfolios)
best_portfolio = portfolio_df[portfolio_df["予算内"]].sort_values(["便益_万円", "費用_万円"], ascending=[False, True]).iloc[0]
display(actions)
display(best_portfolio.to_frame("最適ポートフォリオ"))

feasible_df = portfolio_df[portfolio_df["予算内"]]
fig, ax = plt.subplots(figsize=(7.6, 4.2))
ax.scatter(feasible_df["費用_万円"], feasible_df["便益_万円"], alpha=0.55, label="実行可能案")
ax.scatter(best_portfolio["費用_万円"], best_portfolio["便益_万円"], s=130, marker="*", color="tab:red", label="最適案")
ax.axvline(budget, color="black", linestyle="--", label="予算上限")
ax.set_title("改善施策ポートフォリオの費用と期待便益")
ax.set_xlabel("費用（万円）")
ax.set_ylabel("期待年間便益（万円）")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	施策	費用_万円	期待年間便益_万円
0	検査応援	90	170
1	研削予防保全	130	240
2	安全在庫追加	70	105
3	段取り短縮	160	260
4	センサー増設	110	150

	最適ポートフォリオ
施策	検査応援・研削予防保全・安全在庫追加
費用_万円	290
便益_万円	515
予算内	True

svg

結果の読み取り

予算内で期待便益が最大の組合せを選べます。ただし表の便益は点推定です。モンテカルロ法で下振れを評価し、施策同士の効果重複、導入リードタイム、現場要員、停止工事期間を制約へ加えます。最適化は経営判断を自動化するものではなく、前提とトレードオフを明示して合意形成を助ける道具です。

No.080：工場シミュレーション

実務での意味

局所改善を工場全体で評価します。研削を速くしても検査が詰まれば出荷は増えません。3工程を流れる製品、工程ごとのばらつき、故障停止を統合し、現状と改善案を同じ需要・乱数条件で比較します。

分析・モデル化の考え方

ジョブ $j$ 、工程 $i$ の完了時刻を

C_{ij}=\max(C_{i-1,j},C_{i,j-1})+S_{ij}+D_{ij}

とします。 $S_{ij}$ は加工時間、 $D_{ij}$ は故障などによる停止です。スループット、リードタイム、工程待ち、納期達成率を複数反復で比較します。共通乱数を使うと、施策差以外のランダムな揺れを抑えられます。

Pythonで確認する

def factory_simulation(n_jobs=90, scenario="現状", seed=80):
    local_rng = np.random.default_rng(seed)
    arrivals = np.cumsum(local_rng.exponential(7.0, n_jobs))
    base = np.array([7.0, 8.8, 5.2])
    sigma = np.array([0.18, 0.24, 0.20])
    failure_prob = np.array([0.015, 0.040, 0.018])
    repair_mean = np.array([22, 38, 18])
    if scenario == "改善案":
        base = base * np.array([1.0, 0.92, 0.88])
        failure_prob = failure_prob * np.array([1.0, 0.45, 1.0])
    available = np.zeros(3)
    records = []
    for job in range(n_jobs):
        previous_finish = arrivals[job]
        waits = []
        for i in range(3):
            start = max(previous_finish, available[i])
            waits.append(start - previous_finish)
            service = local_rng.lognormal(np.log(base[i]) - 0.5 * sigma[i]**2, sigma[i])
            downtime = local_rng.exponential(repair_mean[i]) if local_rng.random() < failure_prob[i] else 0.0
            finish = start + service + downtime
            available[i] = finish
            previous_finish = finish
        lead_time = previous_finish - arrivals[job]
        records.append([job + 1, arrivals[job], previous_finish, lead_time, *waits])
    result = pd.DataFrame(records, columns=["ジョブ", "投入", "完了", "リードタイム", "旋削待ち", "研削待ち", "検査待ち"])
    horizon = result["完了"].max() - result["投入"].min()
    kpi = {
        "スループット（個/時）": n_jobs / horizon * 60,
        "平均リードタイム（分）": result["リードタイム"].mean(),
        "95%リードタイム（分）": result["リードタイム"].quantile(0.95),
        "納期60分以内率": np.mean(result["リードタイム"] <= 60),
    }
    return result, kpi

scenario_rows = []
lead_samples = {"現状": [], "改善案": []}
for scenario in lead_samples:
    for rep in range(120):
        result, kpi = factory_simulation(scenario=scenario, seed=8000 + rep)
        scenario_rows.append({"シナリオ": scenario, "反復": rep, **kpi})
        lead_samples[scenario].append(kpi["平均リードタイム（分）"])
scenario_df = pd.DataFrame(scenario_rows)
factory_summary = scenario_df.groupby("シナリオ").agg({
    "スループット（個/時）": "mean",
    "平均リードタイム（分）": "mean",
    "95%リードタイム（分）": "mean",
    "納期60分以内率": "mean",
})
display(factory_summary.style.format({
    "スループット（個/時）": "{:.2f}", "平均リードタイム（分）": "{:.1f}",
    "95%リードタイム（分）": "{:.1f}", "納期60分以内率": "{:.1%}",
}))

fig, ax = plt.subplots(figsize=(7.6, 4.1))
ax.boxplot([lead_samples["現状"], lead_samples["改善案"]], tick_labels=["現状", "改善案"])
ax.set_title("工場シミュレーションによる平均リードタイム比較")
ax.set_xlabel("シナリオ")
ax.set_ylabel("反復ごとの平均リードタイム（分）")
ax.grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()

	スループット（個/時）	平均リードタイム（分）	95%リードタイム（分）	納期60分以内率
シナリオ
改善案	6.54	129.1	222.6	22.4%
現状	5.72	188.4	336.1	15.7%

svg

結果の読み取り

改善案では研削の故障確率と加工時間、検査時間を同時に改善し、スループット、平均・95%リードタイム、納期内率を比較できます。箱ひげ図の重なりは効果の不確実性も示します。平均差だけで採用せず、改善費用、効果の信頼区間、繁忙期、長期停止、製品ミックスを変えた感度分析を行います。局所KPIが改善しても全体KPIが変わらなければ、別の制約工程を探します。

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

状態を定義すると将来を議論できる：状態遷移行列とマルコフモデルで、設備群の短期変化と長期構成を分けて評価できます。
能力と変動は別に扱う：生産計画は平均能力の配分、待ち行列と離散事象シミュレーションは到着・処理のばらつきを扱います。
在庫はサービス水準との交換条件である：安全在庫を固定値でなく、需要・リードタイム・欠品費から設計します。
平均値だけで投資を決めない：モンテカルロ法で赤字確率、分位点、納期未達など下振れを確認します。
デジタルツインは継続的な校正プロセスである：現場実績との差異を測り、モデル版と入力品質を管理します。
最適化とシミュレーションを往復する：最適化で候補を絞り、シミュレーションで動的な実現可能性を検証します。
全体KPIで局所改善を評価する：工程稼働率だけでなく、出荷量、リードタイム、仕掛品、利益まで追います。

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

1. 判断とKPIを先に定義する

「検査員を増やすか」「安全在庫を何個にするか」「定修でどの設備を整備するか」など意思決定を明確にします。スループット、納期達成率、仕掛在庫、欠品費、停止損失の計算範囲と単位を揃えます。

2. データの時刻・状態・イベントを整える

設備ID、品番、ロット、工程、加工開始・完了、停止理由、修理完了を共通キーで結びます。欠測とゼロ、計画停止と故障停止、再加工と新規加工を区別し、設備・工程変更の履歴を残します。

3. 段階的にモデルを検証する

手計算や待ち行列を基準にし、必要な場合だけ離散事象へ進みます。工程別平均・分散、日別出荷、WIP、停止回数を実績と照合し、校正に使っていない期間・繁忙期・異常時でも再現性を確認します。

4. 不確実性と感度を示す

需要、加工時間、故障、修理時間の分布と相関を検討します。単一結果でなく分位点、信頼区間、ストレスシナリオを提示し、結論を左右する入力を感度分析で特定します。

5. 意思決定と更新運用を設計する

モデルの承認者、入力更新頻度、再校正条件、改善案の実施責任者を決めます。実施後のKPIを記録し、予測便益と実績便益の差を次回モデルへ戻します。シミュレーションは現場と経営の共通言語として管理します。

まとめ

No.071〜No.080では、架空の精密部品ラインを題材に、状態遷移行列、在庫モデル、生産計画、待ち行列、マルコフモデル、離散事象シミュレーション、モンテカルロ法、デジタルツイン、改善施策の最適化、工場全体シミュレーションを確認しました。

製造業でシミュレーションを活かす鍵は、精巧なモデルを作ること自体ではありません。意思決定を定義し、現場実績でモデルを検証し、不確実性を示したうえで、施策後の結果を次の判断へ戻すことです。

法人向けのご相談

数理工房では、製造業の状態遷移分析、在庫・生産計画、待ち行列分析、離散事象シミュレーション、モンテカルロ評価、デジタルツイン、数理最適化について、課題整理からデータ設計、PoC、妥当性確認、現場運用までをご支援しています。

「設備増強前にボトルネックを確認したい」「在庫と欠品のバランスを定量化したい」「工場シミュレーションを投資判断へつなげたい」「既存シミュレーションと現場実績の差を改善したい」といった課題をご相談いただけます。

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

工場シミュレーションをPythonで実践｜在庫・待ち行列・生産計画の意思決定

変動に強い加工ラインを設計する

状態遷移・在庫・待ち行列・工場シミュレーション：100本ノック No.071〜No.080

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

現場でよくある状況

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

今回扱うノックの全体像

Python 環境の準備

架空データの作成

No.071：状態遷移行列

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.072：在庫モデル

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.073：生産計画

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.074：待ち行列

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.075：マルコフモデル

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.076：離散事象シミュレーション

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.077：モンテカルロ法

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.078：デジタルツイン

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.079：最適化

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.080：工場シミュレーション

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

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

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

1. 判断とKPIを先に定義する

2. データの時刻・状態・イベントを整える

3. 段階的にモデルを検証する

4. 不確実性と感度を示す

5. 意思決定と更新運用を設計する

まとめ

法人向けのご相談

Related Posts

製造業のPython高速化入門｜CuPy・JAX・Polarsと行列計算の実務

製造業の品質・需要・生産計画をPythonで統合｜実務で学ぶ行列分析

製造業の予知保全をPythonで学ぶ｜Transformer・GNN・行列分解の実務活用

製造業の品質予測をPythonで学ぶ｜Ridge・Logistic回帰・Attention入門

製造業の工程ネットワーク分析入門｜PCA・PageRank・GCNをPythonで学ぶ