100本ノック / 行列 / 行列100本ノック
製造業の品質・需要・生産計画をPythonで統合｜実務で学ぶ行列分析
品質・需要・供給制約を統合し、週次の生産判断を設計する
品質データ、設備センサー、需要、景況感、部品表を行列として捉え、異常の発見から生産配分までを一つの判断プロセスにまとめます。架空の精密部品メーカーを題材に、分析手法と Python ライブラリを「何を決めるために使うのか」という視点で確認します。
[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。
はじめに：この記事で扱う製造業の実務課題
生産会議では、「不良率が上がった」「設備の挙動がいつもと違う」「来月の受注が読みにくい」「材料が足りない」といった情報が同時に持ち込まれます。個別の分析が正しくても、単位・時間軸・判断基準が揃わなければ、増産、点検、外注、在庫積み増しの選択にはつながりません。
本記事では、品質保証、生産管理、営業、調達が持つデータを、共通の行列計算と KPI に落とし込みます。目的は予測精度だけを競うことではなく、制約と損失を明示し、再現可能な週次判断を作ることです。
現場でよくある状況
- 品質管理図は更新されているが、設備点検や生産配分のルールと接続していない
- センサーごとの閾値監視では、複数項目が同時に少しずつ崩れる異常を見逃す
- 需要予測が一点予測のままで、欠品費用と余剰在庫費用の非対称性を扱えない
- 営業の景況感、実受注、生産能力、部材制約が別々の表で管理されている
- 高度なライブラリを導入しても、入力・出力・責任者が定義されず定着しない
そこで本 notebook では、各分析の出力を「次に取る行動」と「判断を見直す条件」まで結び付けます。
なぜこの問題は判断が難しいのか
第一に、品質・需要・設備は確率的に変動します。管理限界を超えた一点が必ず故障を意味するわけではなく、予測値にも誤差があります。第二に、KPI 間にトレードオフがあります。検査を増やせば流出は抑えられますが、納期と検査費用に影響します。第三に、相関と因果は異なります。気温と振動が同時に上がっても、設備劣化とは限りません。
したがって、分析結果は自動的な結論ではなく、観測 → 確率評価 → 制約付き選択 → 現場確認 → 学習という閉ループの判断材料として設計する必要があります。
今回扱うノックの全体像
| No. | テーマ | 今回の問い | 主な出力 |
|---|---|---|---|
| 081 | 品質管理 | 不良率の変化は通常変動か | p 管理図・管理限界 |
| 082 | 異常検知 | 複数センサーの組合せは異常か | マハラノビス距離 |
| 083 | 需要予測 | 13週先の需要と誤差はどの程度か | 回帰予測・MAE |
| 084 | 意思決定支援 | 需要シナリオ別にどの行動を選ぶか | 期待損失行列 |
| 085 | 製造業DI | 現場の方向感は改善か悪化か | DI と移動平均 |
| 086 | NumPy | SKU×ライン計画を高速に試算できるか | 良品数行列 |
| 087 | SciPy | 能力・需要制約下の配分は何か | 線形計画 |
| 088 | PyTorch | 非線形な異常パターンを学習できるか | 異常確率モデル |
| 089 | Sparse Matrix | 大規模 BOM を省メモリで計算できるか | CSR 行列・所要量 |
| 090 | NetworkX | 供給網の重要ノードはどこか | 中心性・依存グラフ |
No.081〜085 は意思決定の論点、No.086〜090 はそれを実装する道具です。後半もライブラリ紹介に留めず、前半の業務判断に接続します。
Python 環境の準備
外部データは使いません。乱数生成器の seed を固定し、同じ環境なら同じ表とグラフを再現できるようにします。グラフ内の英語表記は、Markdown 変換先での日本語フォント差による文字化けを避けるためです。
%matplotlib inline
import platform
import warnings
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd
import scipy
from scipy import optimize, sparse, stats
import torch
SEED = 42
rng = np.random.default_rng(SEED)
np.set_printoptions(precision=3, suppress=True)
pd.options.display.float_format = "{:,.3f}".format
plt.rcParams["figure.figsize"] = (9, 4.5)
print(f"Python {platform.python_version()}")
print(f"NumPy {np.__version__} / pandas {pd.__version__} / SciPy {scipy.__version__}")
print(f"PyTorch {torch.__version__} / NetworkX {nx.__version__}")Python 3.13.1
NumPy 2.5.1 / pandas 3.0.3 / SciPy 1.18.0
PyTorch 2.13.0 / NetworkX 3.6.1架空データの作成
架空の精密部品メーカーが、製品 A・B・C を2ラインで生産しているとします。品質日報180日分、設備センサー240点、週次需要104週分を生成します。後半30日には工程条件の変化を、センサーの一部には複合的な異常を意図的に含めます。
現実には MES、QMS、設備 PLC、販売管理から取得する項目ですが、ここでは列の意味と粒度を明示した一枚の分析用データに整えます。
# 品質日報：後半にラインL2の不良率が上昇する架空設定
dates = pd.date_range("2025-01-01", periods=180, freq="D")
lines = np.where(np.arange(len(dates)) % 2 == 0, "L1", "L2")
inspected = rng.integers(180, 260, len(dates))
temperature = 23 + 3 * np.sin(np.arange(len(dates)) / 18) + rng.normal(0, 0.8, len(dates))
base_rate = 0.018 + 0.002 * (lines == "L2")
shift = ((np.arange(len(dates)) >= 150) & (lines == "L2")) * 0.025
defect_rate_true = np.clip(base_rate + shift + 0.0015 * np.maximum(temperature - 25, 0), 0, 0.12)
defects = rng.binomial(inspected, defect_rate_true)
quality_df = pd.DataFrame({
"date": dates, "line": lines, "inspected": inspected,
"defects": defects, "temperature_c": temperature,
})
quality_df["defect_rate"] = quality_df["defects"] / quality_df["inspected"]
# 設備センサー：正常データを作り、一部に多変量のずれを加える
n_sensor = 240
sensor_mean = np.array([2.2, 58.0, 7.5]) # vibration, temperature, current
sensor_cov = np.array([[0.10, 0.35, 0.08], [0.35, 5.0, 0.45], [0.08, 0.45, 0.30]])
sensor_values = rng.multivariate_normal(sensor_mean, sensor_cov, n_sensor)
true_anomaly = np.zeros(n_sensor, dtype=int)
anomaly_idx = np.array([211, 218, 225, 232, 237])
sensor_values[anomaly_idx] += np.array([1.0, 3.5, 1.4])
true_anomaly[anomaly_idx] = 1
sensor_df = pd.DataFrame(sensor_values, columns=["vibration", "temp_c", "current_a"])
sensor_df["time"] = pd.date_range("2025-06-01", periods=n_sensor, freq="h")
sensor_df["true_anomaly"] = true_anomaly
# 週次需要：傾向、年周期、販促効果、ノイズを持つ3製品
weeks = pd.date_range("2024-01-01", periods=104, freq="W-MON")
t = np.arange(len(weeks))
product_params = {"A": (430, 1.1, 55), "B": (310, 0.5, 38), "C": (220, 1.5, 30)}
rows = []
for product, (level, trend, amp) in product_params.items():
promo = ((t % 26) == 20).astype(int)
demand = level + trend * t + amp * np.sin(2 * np.pi * t / 52) + 45 * promo + rng.normal(0, 22, len(t))
rows.extend(zip(weeks, [product] * len(t), np.maximum(np.rint(demand), 0).astype(int), promo))
demand_df = pd.DataFrame(rows, columns=["week", "product", "demand", "promo"])
display(quality_df.tail(6))
display(sensor_df.loc[anomaly_idx, ["time", "vibration", "temp_c", "current_a"]])
display(demand_df.groupby("product")["demand"].agg(["mean", "std", "min", "max"]))| date | line | inspected | defects | temperature_c | defect_rate | |
|---|---|---|---|---|---|---|
| 174 | 2025-06-24 | L1 | 235 | 7 | 22.467 | 0.030 |
| 175 | 2025-06-25 | L2 | 216 | 11 | 21.677 | 0.051 |
| 176 | 2025-06-26 | L1 | 237 | 2 | 22.340 | 0.008 |
| 177 | 2025-06-27 | L2 | 192 | 5 | 22.618 | 0.026 |
| 178 | 2025-06-28 | L1 | 252 | 3 | 21.781 | 0.012 |
| 179 | 2025-06-29 | L2 | 220 | 17 | 21.792 | 0.077 |
| time | vibration | temp_c | current_a | |
|---|---|---|---|---|
| 211 | 2025-06-09 19:00:00 | 3.741 | 63.986 | 8.928 |
| 218 | 2025-06-10 02:00:00 | 3.254 | 59.705 | 8.781 |
| 225 | 2025-06-10 09:00:00 | 3.232 | 62.033 | 9.007 |
| 232 | 2025-06-10 16:00:00 | 3.367 | 61.797 | 9.584 |
| 237 | 2025-06-10 21:00:00 | 2.821 | 59.198 | 8.830 |
| mean | std | min | max | |
|---|---|---|---|---|
| product | ||||
| A | 485.712 | 47.753 | 376 | 581 |
| B | 336.740 | 37.064 | 243 | 425 |
| C | 300.548 | 46.951 | 195 | 414 |
No.081：品質管理 — p 管理図で工程変化を見つける
実務での意味
日々の不良率の上下をすべて問題視すると、現場は調査疲れを起こします。一方、平均だけを見ると小さな工程変化を見逃します。p 管理図は検査数が日ごとに異なる計数データについて、偶然のばらつきと調査すべき変化を分ける一次スクリーニングです。
分析・モデル化の考え方
基準期間の総不良数を総検査数で割った中心線を とし、日 の検査数を とすると、3シグマ管理限界は次式です。
最初の150日を安定していた基準期間と仮定し、直近30日を監視します。管理限界超過は原因確定ではなく、ロット、材料、設備、測定系を確認する合図です。
Pythonで確認する
baseline = quality_df.iloc[:150]
monitor = quality_df.iloc[150:].copy()
p_bar = baseline["defects"].sum() / baseline["inspected"].sum()
sigma = np.sqrt(p_bar * (1 - p_bar) / monitor["inspected"])
monitor["ucl"] = p_bar + 3 * sigma
monitor["lcl"] = np.maximum(0, p_bar - 3 * sigma)
monitor["signal"] = (monitor["defect_rate"] > monitor["ucl"]) | (monitor["defect_rate"] < monitor["lcl"])
summary_081 = monitor.groupby("line").agg(
days=("date", "size"), mean_defect_rate=("defect_rate", "mean"), signals=("signal", "sum")
)
display(summary_081)
fig, ax = plt.subplots()
ax.plot(monitor["date"], monitor["defect_rate"], marker="o", label="Daily defect rate")
ax.plot(monitor["date"], monitor["ucl"], "r--", label="UCL (3 sigma)")
ax.plot(monitor["date"], monitor["lcl"], "r--", label="LCL (3 sigma)")
ax.axhline(p_bar, color="black", linestyle=":", label="Baseline mean")
signals = monitor[monitor["signal"]]
ax.scatter(signals["date"], signals["defect_rate"], color="red", s=70, zorder=3, label="Signal")
ax.set_title("p-Chart for Recent Quality Performance")
ax.set_xlabel("Date")
ax.set_ylabel("Defect rate")
ax.grid(True, alpha=0.3)
ax.legend()
fig.tight_layout()
plt.show()| days | mean_defect_rate | signals | |
|---|---|---|---|
| line | |||
| L1 | 15 | 0.019 | 0 |
| L2 | 15 | 0.045 | 8 |
結果の読み取り
直近期間では L2 を中心に管理限界超過が現れます。これは「L2 が故障した」という結論ではなく、基準期間と同じ工程状態とは考えにくいという信号です。超過日の材料ロット、段取り、温度、測定器を層別し、原因確認まで不適合品の隔離や検査強化を検討します。基準期間自体に異常が混じると限界が広がるため、中心線の承認プロセスも必要です。
No.082：異常検知 — センサーの組合せで設備を監視する
実務での意味
振動、温度、電流は互いに関係します。各項目が単独閾値内でも、通常は起きにくい組合せなら点検価値があります。多変量異常検知は、その「組合せの違和感」を定量化します。
分析・モデル化の考え方
正常期間の平均 と共分散行列 を使い、観測 のマハラノビス距離を
で測ります。3変量の正常分布を仮定した参考閾値として 分布の99%点を使います。実務では設備状態別に基準を分け、誤警報率と見逃し費用から閾値を校正します。
Pythonで確認する
features = ["vibration", "temp_c", "current_a"]
reference = sensor_df.loc[:179, features].to_numpy()
mu = reference.mean(axis=0)
cov_inv = np.linalg.pinv(np.cov(reference, rowvar=False))
delta = sensor_df[features].to_numpy() - mu
sensor_df["mahalanobis_sq"] = np.einsum("ij,jk,ik->i", delta, cov_inv, delta)
threshold = stats.chi2.ppf(0.99, df=len(features))
sensor_df["detected"] = sensor_df["mahalanobis_sq"] > threshold
confusion_082 = pd.crosstab(sensor_df["true_anomaly"], sensor_df["detected"],
rownames=["Actual"], colnames=["Detected"])
display(confusion_082)
display(sensor_df.nlargest(8, "mahalanobis_sq")[["time", *features, "mahalanobis_sq", "detected"]])
fig, ax = plt.subplots()
ax.plot(sensor_df["time"], sensor_df["mahalanobis_sq"], label="Squared distance")
ax.axhline(threshold, color="red", linestyle="--", label="99% threshold")
ax.scatter(sensor_df.loc[sensor_df["detected"], "time"],
sensor_df.loc[sensor_df["detected"], "mahalanobis_sq"], color="red", s=45)
ax.set_title("Multivariate Sensor Anomaly Score")
ax.set_xlabel("Time")
ax.set_ylabel("Squared Mahalanobis distance")
ax.grid(True, alpha=0.3)
ax.legend()
fig.tight_layout()
plt.show()| Detected | False | True |
|---|---|---|
| Actual | ||
| 0 | 232 | 3 |
| 1 | 1 | 4 |
| time | vibration | temp_c | current_a | mahalanobis_sq | detected | |
|---|---|---|---|---|---|---|
| 211 | 2025-06-09 19:00:00 | 3.741 | 63.986 | 8.928 | 24.886 | True |
| 232 | 2025-06-10 16:00:00 | 3.367 | 61.797 | 9.584 | 20.604 | True |
| 170 | 2025-06-08 02:00:00 | 2.120 | 57.330 | 5.488 | 15.740 | True |
| 218 | 2025-06-10 02:00:00 | 3.254 | 59.705 | 8.781 | 13.740 | True |
| 22 | 2025-06-01 22:00:00 | 2.021 | 64.686 | 8.007 | 13.585 | True |
| 225 | 2025-06-10 09:00:00 | 3.232 | 62.033 | 9.007 | 13.268 | True |
| 234 | 2025-06-10 18:00:00 | 2.133 | 54.265 | 5.851 | 11.882 | True |
| 205 | 2025-06-09 13:00:00 | 1.281 | 52.630 | 6.376 | 10.775 | False |
結果の読み取り
注入した複合異常は距離上位に現れ、単一センサーの固定閾値よりも相関構造を使えることが分かります。ただし、統計的異常と故障は同義ではありません。上位時刻の運転モード、工具交換、暖機運転を照合し、点検結果をラベルとして戻す運用が必要です。共分散が不安定な場合に逆行列ではなく擬似逆行列を使っても、データ不足そのものが解消されるわけではありません。
No.083：需要予測 — 13週先と予測誤差を見積もる
実務での意味
需要予測は売上を当てるためだけでなく、材料発注、要員計画、安全在庫の前提です。重要なのは予測値と実績の差を継続評価し、どの時間軸なら意思決定に使えるかを把握することです。
分析・モデル化の考え方
製品ごとに、切片、時間トレンド、年周期の正弦・余弦、販促フラグを説明変数とする線形モデルを作ります。
最後の13週をテスト期間として学習に使わず、MAE と WAPE を計算します。ここでは説明可能な基準モデルを採用し、複雑なモデルへ進む前の比較軸にします。
Pythonで確認する
forecast_rows, metric_rows = [], []
for product, group in demand_df.groupby("product"):
group = group.sort_values("week").reset_index(drop=True)
tt = np.arange(len(group))
X = np.column_stack([
np.ones(len(group)), tt,
np.sin(2 * np.pi * tt / 52), np.cos(2 * np.pi * tt / 52), group["promo"]
])
split = len(group) - 13
beta, *_ = np.linalg.lstsq(X[:split], group.loc[:split-1, "demand"], rcond=None)
pred = X[split:] @ beta
actual = group.loc[split:, "demand"].to_numpy()
metric_rows.append({
"product": product,
"MAE": np.mean(np.abs(actual - pred)),
"WAPE_pct": 100 * np.sum(np.abs(actual - pred)) / np.sum(actual),
})
forecast_rows.extend(zip(group.loc[split:, "week"], [product] * 13, actual, pred))
forecast_df = pd.DataFrame(forecast_rows, columns=["week", "product", "actual", "forecast"])
metrics_083 = pd.DataFrame(metric_rows).set_index("product")
display(metrics_083)
fig, ax = plt.subplots()
for product, group in forecast_df.groupby("product"):
ax.plot(group["week"], group["actual"], marker="o", label=f"{product} actual")
ax.plot(group["week"], group["forecast"], linestyle="--", label=f"{product} forecast")
ax.set_title("13-Week Holdout Demand Forecast")
ax.set_xlabel("Week")
ax.set_ylabel("Units")
ax.grid(True, alpha=0.3)
ax.legend(ncol=2)
fig.tight_layout()
plt.show()| MAE | WAPE_pct | |
|---|---|---|
| product | ||
| A | 19.536 | 3.856 |
| B | 16.462 | 4.875 |
| C | 26.440 | 7.686 |
結果の読み取り
製品別 MAE は「週当たり何個外れるか」、WAPE は製品規模をならした誤差です。計画担当者は一点予測を確定値として扱わず、少なくとも MAE 相当の変動を能力・在庫の感度分析へ渡します。販促予定の未確定、欠品による観測需要の打ち切り、新製品の履歴不足は別途扱う必要があります。時系列分割を使ったのは、未来情報を学習へ混ぜないためです。
No.084：意思決定支援 — シナリオと損失行列で行動を選ぶ
実務での意味
予測が同じでも、欠品を重く見る企業と在庫を重く見る企業では最適行動が異なります。意思決定支援では予測値だけでなく、選択肢、起こり得る状態、確率、損失を並べ、判断理由を監査可能にします。
分析・モデル化の考え方
行動 と需要シナリオ の損失を 、シナリオ確率を とすると、期待損失は
です。製品 A の次週需要を低・中・高の3状態とし、通常生産、検査強化、残業増産を比較します。損失には欠品、余剰、追加運用、不良流出の代理費用を含めます。金額は架空の相対値です。
Pythonで確認する
actions = pd.DataFrame({
"action": ["通常生産", "検査強化", "残業増産"],
"planned_units": [525, 500, 590],
"yield_rate": [0.965, 0.985, 0.960],
"fixed_cost": [0, 55, 90],
})
states = np.array([480, 535, 600])
state_prob = np.array([0.25, 0.50, 0.25])
good_units = actions["planned_units"].to_numpy() * actions["yield_rate"].to_numpy()
# 損失単位：欠品3、余剰1、流出代理費用6、固定費
shortage = np.maximum(states[None, :] - good_units[:, None], 0)
surplus = np.maximum(good_units[:, None] - states[None, :], 0)
escape_proxy = actions["planned_units"].to_numpy()[:, None] * (1 - actions["yield_rate"].to_numpy()[:, None]) * 6
loss_matrix = 3 * shortage + surplus + escape_proxy + actions["fixed_cost"].to_numpy()[:, None]
expected_loss = loss_matrix @ state_prob
decision_084 = pd.DataFrame(loss_matrix, index=actions["action"],
columns=["需要_低", "需要_中", "需要_高"])
decision_084["期待損失"] = expected_loss
display(decision_084.round(1).sort_values("期待損失"))
fig, ax = plt.subplots()
im = ax.imshow(loss_matrix, cmap="YlOrRd")
for i in range(loss_matrix.shape[0]):
for j in range(loss_matrix.shape[1]):
ax.text(j, i, f"{loss_matrix[i, j]:.0f}", ha="center", va="center")
ax.set_xticks(range(3), ["Low", "Mid", "High"])
ax.set_yticks(range(3), ["Normal", "Inspect", "Overtime"])
ax.set_title("Loss Matrix by Demand Scenario")
ax.set_xlabel("Demand scenario")
ax.set_ylabel("Production action")
ax.grid(False)
fig.colorbar(im, ax=ax, label="Relative loss")
fig.tight_layout()
plt.show()| 需要_低 | 需要_中 | 需要_高 | 期待損失 | |
|---|---|---|---|---|
| action | ||||
| 通常生産 | 136.900 | 195.400 | 390.400 | 229.500 |
| 検査強化 | 112.500 | 227.500 | 422.500 | 247.500 |
| 残業増産 | 318.000 | 263.000 | 332.400 | 294.100 |
結果の読み取り
期待損失が最小の行動が、設定した確率と費用の下での推奨案です。ただし、これは責任者の判断を置き換える答えではありません。高需要確率、欠品単価、検査による歩留まり改善を変え、推奨がどこで切り替わるかを確認することが重要です。労務上限や顧客別優先度のような絶対制約は、損失へ埋め込まず制約条件として扱います。
No.085：製造業DI — 方向感を共通言語にする
実務での意味
受注額がまだ変化していなくても、引き合い、見積、納期照会などの先行情報は現場にあります。DI（Diffusion Index）は、改善と回答した割合から悪化と回答した割合を引き、方向感を一つの指標にします。営業・調達・生産の対話を始める早期指標として使えます。
分析・モデル化の考え方
回答企業数を 、改善回答数を 、悪化回答数を とすると、
です。「横ばい」は分母には含めますが差には直接入りません。架空の60社アンケートを24か月生成し、3か月移動平均で短期ノイズを抑えます。標本構成が変わると時系列比較が崩れるため、継続回答率も管理対象です。
Pythonで確認する
months = pd.date_range("2024-01-01", periods=24, freq="MS")
latent = 0.12 * np.sin(np.arange(24) / 3.2) + np.linspace(-0.06, 0.12, 24)
survey_rows = []
for month, score in zip(months, latent):
p_improve = np.clip(0.28 + score, 0.08, 0.65)
p_worsen = np.clip(0.27 - score, 0.08, 0.65)
p_same = 1 - p_improve - p_worsen
counts = rng.multinomial(60, [p_improve, p_same, p_worsen])
survey_rows.append((month, *counts))
di_df = pd.DataFrame(survey_rows, columns=["month", "improve", "same", "worsen"])
di_df["DI"] = 100 * (di_df["improve"] - di_df["worsen"]) / di_df[["improve", "same", "worsen"]].sum(axis=1)
di_df["DI_3m"] = di_df["DI"].rolling(3).mean()
display(di_df.tail(8))
fig, ax = plt.subplots()
colors = np.where(di_df["DI"] >= 0, "tab:blue", "tab:red")
ax.bar(di_df["month"], di_df["DI"], width=20, color=colors, alpha=0.55, label="Monthly DI")
ax.plot(di_df["month"], di_df["DI_3m"], color="black", marker="o", label="3-month average")
ax.axhline(0, color="gray", linewidth=1)
ax.set_title("Manufacturing Sentiment Diffusion Index")
ax.set_xlabel("Month")
ax.set_ylabel("DI (points)")
ax.grid(True, axis="y", alpha=0.3)
ax.legend()
fig.tight_layout()
plt.show()| month | improve | same | worsen | DI | DI_3m | |
|---|---|---|---|---|---|---|
| 16 | 2025-05-01 | 13 | 31 | 16 | -5.000 | -2.778 |
| 17 | 2025-06-01 | 15 | 28 | 17 | -3.333 | -2.778 |
| 18 | 2025-07-01 | 15 | 26 | 19 | -6.667 | -5.000 |
| 19 | 2025-08-01 | 18 | 30 | 12 | 10.000 | 0.000 |
| 20 | 2025-09-01 | 18 | 31 | 11 | 11.667 | 5.000 |
| 21 | 2025-10-01 | 23 | 30 | 7 | 26.667 | 16.111 |
| 22 | 2025-11-01 | 25 | 28 | 7 | 30.000 | 22.778 |
| 23 | 2025-12-01 | 35 | 21 | 4 | 51.667 | 36.111 |
結果の読み取り
DI が0を上回れば改善回答が悪化回答より多く、3か月平均の上向きは方向感の持続を示します。ただし DI は変化の広がりであり、受注額の増減幅ではありません。したがって、実受注、キャンセル、見積転換率と並べて先行性を検証します。回答数が少ない部門別 DI は振れやすいため、数値と回答件数を必ず併記します。
No.086：NumPy — SKU×ライン計画を行列で試算する
実務での意味
生産計画は製品、ライン、シフト、歩留まりの組合せです。セルごとの手計算では条件変更のたびに式が増え、転記ミスも起きます。NumPy の配列演算なら、計画数量と歩留まりを同じ形の行列として扱い、複数シナリオを一括比較できます。
分析・モデル化の考え方
計画数量行列を 、歩留まり行列を とすると、要素積
が SKU×ラインの期待良品数です。行方向に合計すればライン別、列方向に合計すれば SKU 別の良品数になります。行列積と要素積を取り違えないことが重要です。
Pythonで確認する
product_names = np.array(["A", "B", "C"])
line_names = np.array(["L1", "L2"])
plan = np.array([[300, 180, 100], [240, 190, 150]]) # rows: lines, columns: products
yield_matrix = np.array([[0.985, 0.975, 0.965], [0.970, 0.960, 0.950]])
good_matrix = plan * yield_matrix
numpy_086 = pd.DataFrame(good_matrix, index=line_names, columns=product_names)
numpy_086["Line total"] = numpy_086.sum(axis=1)
numpy_086.loc["Product total"] = numpy_086.sum(axis=0)
display(numpy_086)
yield_scenarios = np.stack([yield_matrix, yield_matrix - 0.01, yield_matrix + 0.005])
scenario_totals = (plan[None, :, :] * yield_scenarios).sum(axis=(1, 2))
display(pd.DataFrame({"scenario": ["base", "downside", "improved"],
"expected_good_units": scenario_totals}))| A | B | C | Line total | |
|---|---|---|---|---|
| L1 | 295.500 | 175.500 | 96.500 | 567.500 |
| L2 | 232.800 | 182.400 | 142.500 | 557.700 |
| Product total | 528.300 | 357.900 | 239.000 | 1,125.200 |
| scenario | expected_good_units | |
|---|---|---|
| 0 | base | 1,125.200 |
| 1 | downside | 1,113.600 |
| 2 | improved | 1,131.000 |
結果の読み取り
期待良品数をライン別・製品別に同じ計算から得られるため、会議中の条件変更にも一貫して対応できます。下振れシナリオとの差は、追加検査や予備能力を検討する数量的な根拠です。ただし期待値は個々のロットの保証値ではありません。歩留まりの推定期間、品種切替、再加工可能数を明記し、元データの単位と軸順をテストします。
No.087：SciPy — 制約付き生産配分を最適化する
実務での意味
需要を満たしたくても、ライン時間、材料、適合設備には上限があります。経験的な配分は説明しにくく、条件が変わるたびに再調整が必要です。線形計画は、利益と制約を明文化し、実行可能な基準案を短時間で作ります。
分析・モデル化の考え方
各ライン・製品の生産量を 、限界利益を とし、 を解きます。SciPy の linprog は最小化問題なので係数の符号を反転します。ライン時間と需要上限を不等式制約 として与えます。
Pythonで確認する
# 変数順：L1-A, L1-B, L1-C, L2-A, L2-B, L2-C
margin = np.array([8.5, 7.0, 6.2, 8.0, 7.4, 6.6])
hours = np.array([0.08, 0.10, 0.13, 0.09, 0.09, 0.11])
demand_cap = np.array([560, 390, 270])
A_ub = []
b_ub = []
A_ub.extend([
[hours[0], hours[1], hours[2], 0, 0, 0],
[0, 0, 0, hours[3], hours[4], hours[5]],
])
b_ub.extend([52, 50])
for p in range(3):
row = np.zeros(6)
row[p] = 1
row[p + 3] = 1
A_ub.append(row)
b_ub.append(demand_cap[p])
result = optimize.linprog(-margin, A_ub=np.array(A_ub), b_ub=np.array(b_ub),
bounds=[(0, None)] * 6, method="highs")
allocation = result.x.reshape(2, 3)
allocation_df = pd.DataFrame(allocation, index=line_names, columns=product_names)
allocation_df["Used hours"] = (allocation * hours.reshape(2, 3)).sum(axis=1)
display(allocation_df.round(1))
print(f"Optimization success: {result.success}")
print(f"Maximum contribution margin: {-result.fun:,.1f}")
print("Demand utilization:", np.round(allocation.sum(axis=0) / demand_cap, 3))| A | B | C | Used hours | |
|---|---|---|---|---|
| L1 | 560.000 | 72.000 | 0.000 | 52.000 |
| L2 | 0.000 | 318.000 | 194.400 | 50.000 |
Optimization success: True
Maximum contribution margin: 8,900.0
Demand utilization: [1. 1. 0.72]結果の読み取り
最適解は限られたライン時間を限界利益と加工時間の組合せに基づいて配分します。需要上限まで作られない製品があれば、能力制約下での機会損失を示します。実務では整数ロット、段取り回数、最低生産量、要員技能、保全時間を追加します。また、数理解が現場で実行できない場合は現場を責めるのではなく、欠けている制約をモデルへ戻すことが改善です。
No.088：PyTorch — 非線形な異常確率を学習する
実務での意味
設備異常は「振動が高い、かつ電流も高い」のような相互作用を持つことがあります。PyTorch はニューラルネットワークと自動微分を用い、単純な線形境界では表しにくいパターンを学習できます。一方、教師ラベルが乏しい現場では過学習に注意が必要です。
分析・モデル化の考え方
センサーを標準化し、3入力→8中間ユニット→1出力の小さなネットワークで異常確率を推定します。損失は二値交差エントロピーです。架空データの異常を増幅複製して学習例を確保しますが、実務では時系列をまたいだ検証と、点検結果に基づく信頼できるラベルが必要です。
Pythonで確認する
torch.manual_seed(SEED)
normal_x = rng.multivariate_normal(sensor_mean, sensor_cov, 500)
abnormal_x = rng.multivariate_normal(sensor_mean + np.array([1.0, 3.5, 1.4]), sensor_cov, 140)
X_raw = np.vstack([normal_x, abnormal_x])
y_raw = np.r_[np.zeros(len(normal_x)), np.ones(len(abnormal_x))]
order = rng.permutation(len(y_raw))
X_raw, y_raw = X_raw[order], y_raw[order]
split = 500
mean_train, std_train = X_raw[:split].mean(0), X_raw[:split].std(0)
X = torch.tensor((X_raw - mean_train) / std_train, dtype=torch.float32)
y = torch.tensor(y_raw[:, None], dtype=torch.float32)
model = torch.nn.Sequential(
torch.nn.Linear(3, 8), torch.nn.ReLU(), torch.nn.Linear(8, 1)
)
optimizer_t = torch.optim.Adam(model.parameters(), lr=0.03)
loss_fn = torch.nn.BCEWithLogitsLoss()
loss_history = []
for epoch in range(250):
optimizer_t.zero_grad()
loss = loss_fn(model(X[:split]), y[:split])
loss.backward()
optimizer_t.step()
loss_history.append(loss.item())
with torch.no_grad():
test_prob = torch.sigmoid(model(X[split:])).numpy().ravel()
test_pred = (test_prob >= 0.5).astype(int)
accuracy = (test_pred == y_raw[split:]).mean()
print(f"Final training loss: {loss_history[-1]:.4f}")
print(f"Holdout accuracy: {accuracy:.3f}")
display(pd.crosstab(pd.Series(y_raw[split:], name="Actual"),
pd.Series(test_pred, name="Predicted")))
fig, ax = plt.subplots()
ax.plot(loss_history)
ax.set_title("PyTorch Training Loss")
ax.set_xlabel("Epoch")
ax.set_ylabel("Binary cross-entropy loss")
ax.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()Final training loss: 0.0603
Holdout accuracy: 0.943| Predicted | 0 | 1 |
|---|---|---|
| Actual | ||
| 0.000 | 103 | 5 |
| 1.000 | 3 | 29 |
結果の読み取り
損失が低下し、保留データでも異常を識別できています。ただし、この精度は生成規則が明確な架空データでの値であり、本番性能の根拠にはなりません。実務導入では設備・期間を分けた外部検証、再現率と適合率、誤警報対応工数、確率校正を確認します。No.082 の説明可能な距離指標を基準モデルとして残し、複雑なモデルの追加価値を測ることが重要です。
No.089：Sparse Matrix — 大規模 BOM を省メモリで計算する
実務での意味
全製品×全部材の BOM 行列は、実際には大半が0です。0まで格納する密行列は、品目数が増えるとメモリと計算時間を浪費します。疎行列は「使用される部材」だけを保存し、所要量展開や影響範囲検索を現実的にします。
分析・モデル化の考え方
BOM 行列 の要素 を製品 1個当たりの部材 使用量、生産計画を とすれば、部材所要量は
です。CSR 形式は行方向のアクセスと行列ベクトル積に向きます。今回は小さな業務例に加え、2,000製品×800部材の架空 BOM で保存メモリを比較します。
Pythonで確認する
materials = ["鋼材", "樹脂", "ベアリング", "センサー", "梱包材"]
bom_dense = np.array([
[2.0, 0.0, 1.0, 0.0, 1.0],
[1.5, 0.4, 2.0, 1.0, 1.0],
[0.8, 0.7, 0.0, 1.0, 1.0],
])
bom_csr = sparse.csr_matrix(bom_dense)
production_plan = np.array([520, 360, 250])
requirements = bom_csr.T @ production_plan
display(pd.DataFrame({"material": materials, "required_quantity": requirements}))
large_bom = sparse.random(2000, 800, density=0.005, format="csr", random_state=SEED,
data_rvs=lambda n: rng.integers(1, 5, n).astype(float))
dense_bytes = np.prod(large_bom.shape) * np.dtype(float).itemsize
sparse_bytes = large_bom.data.nbytes + large_bom.indices.nbytes + large_bom.indptr.nbytes
memory_089 = pd.DataFrame({
"format": ["Dense estimate", "CSR actual"],
"memory_MB": [dense_bytes / 1024**2, sparse_bytes / 1024**2],
"stored_values": [np.prod(large_bom.shape), large_bom.nnz],
})
display(memory_089)
print(f"CSR memory reduction: {(1 - sparse_bytes / dense_bytes) * 100:.1f}%")| material | required_quantity | |
|---|---|---|
| 0 | 鋼材 | 1,780.000 |
| 1 | 樹脂 | 319.000 |
| 2 | ベアリング | 1,240.000 |
| 3 | センサー | 610.000 |
| 4 | 梱包材 | 1,130.000 |
| format | memory_MB | stored_values | |
|---|---|---|---|
| 0 | Dense estimate | 12.207 | 1600000 |
| 1 | CSR actual | 0.099 | 8000 |
CSR memory reduction: 99.2%結果の読み取り
小さな BOM では転置行列との積だけで部材所要量を展開できます。大規模例では、非ゼロ要素だけを持つ CSR が密行列より大幅にメモリを削減します。一方、密度が高い行列や要素ごとの頻繁な更新では疎形式が不利になることもあります。BOM の版、有効日、代替部材、歩留まり、端数発注を別途管理し、単純な所要量をそのまま発注量にしないことが重要です。
No.090：NetworkX — 供給・生産ネットワークの重要点を探す
実務での意味
同じ部材を複数製品が共有すると、一社の供給停止が複数ラインへ波及します。表形式の BOM だけでは、迂回経路、共有部材、橋渡しとなる工程を把握しにくいことがあります。グラフ分析は、ネットワーク構造から優先監視対象を見つけます。
分析・モデル化の考え方
供給先→部材→ライン→製品を有向グラフとして表します。次数中心性は直接接続の多さ、媒介中心性は他ノード間の最短経路上にどれだけ現れるかを示します。中心性は事業影響額や代替日数を含まない構造指標なので、停止損失と併用します。
Pythonで確認する
G = nx.DiGraph()
edges = [
("Supplier-X", "Bearing"), ("Supplier-Y", "Sensor"), ("Supplier-Z", "Steel"),
("Bearing", "L1"), ("Bearing", "L2"), ("Sensor", "L2"),
("Steel", "L1"), ("Steel", "L2"),
("L1", "Product-A"), ("L1", "Product-B"),
("L2", "Product-A"), ("L2", "Product-B"), ("L2", "Product-C"),
]
G.add_edges_from(edges)
degree = nx.degree_centrality(G)
betweenness = nx.betweenness_centrality(G, normalized=True)
centrality_090 = pd.DataFrame({"degree": degree, "betweenness": betweenness}).sort_values(
["betweenness", "degree"], ascending=False
)
display(centrality_090.head(8))
pos = nx.spring_layout(G, seed=SEED)
node_colors = ["tab:orange" if n.startswith("Supplier") else
"tab:green" if n.startswith("Product") else
"tab:blue" if n.startswith("L") else "tab:gray" for n in G.nodes]
fig, ax = plt.subplots(figsize=(10, 6))
nx.draw_networkx(G, pos=pos, ax=ax, node_color=node_colors, node_size=1500,
font_size=8, arrows=True, edge_color="gray")
ax.set_title("Supply-to-Product Dependency Network")
ax.set_xlabel("Network layout (structural, no physical unit)")
ax.set_ylabel("Network layout (structural, no physical unit)")
ax.grid(True, alpha=0.15)
fig.tight_layout()
plt.show()| degree | betweenness | |
|---|---|---|
| L2 | 0.600 | 0.156 |
| Bearing | 0.300 | 0.056 |
| Steel | 0.300 | 0.056 |
| L1 | 0.400 | 0.044 |
| Sensor | 0.200 | 0.044 |
| Product-A | 0.200 | 0.000 |
| Product-B | 0.200 | 0.000 |
| Supplier-X | 0.100 | 0.000 |
結果の読み取り
媒介中心性の高いラインや共通部材は、複数経路をつなぐ構造上の要所です。そこから、複数購買、代替品認定、安全在庫、予防保全の優先候補を作れます。ただし中心性が高いだけで投資を決めず、停止時の粗利影響、復旧時間、検出可能性を重ねます。グラフはマスタの欠落に敏感なため、購買先・代替関係・工程経路の更新責任を決める必要があります。
対象ノックを通して見える実務上の示唆
10テーマを通じて重要なのは、手法を個別に導入するのではなく出力を次工程へ渡すことです。
- p 管理図と多変量異常検知で、通常変動から調査候補を絞る
- 需要予測と DI で、数量見通しと方向感を分けて把握する
- 損失行列で不確実性を行動選択へ翻訳する
- NumPy と SciPy で、歩留まり・能力・需要を整合させた基準計画を作る
- PyTorch は単純モデルを上回る場合に限定して採用し、疎行列とグラフで全社規模の部材・依存関係へ拡張する
分析精度だけでなく、検知から確認までの時間、推奨の採用率、欠品・流出・在庫の総損失を運用 KPI として追うことで、分析が経営成果へ接続しているかを評価できます。
実務導入する場合に必要なこと
- 判断定義：誰が、どの頻度で、どの出力を見て、何を決めるかを先に定義する
- データ契約：設備・品目・ロット ID、単位、時刻、欠損、改訂履歴、責任部門を揃える
- 基準管理：管理限界、異常閾値、予測誤差、損失単価、制約の承認者と見直し周期を決める
- 検証設計：未来情報の混入を避け、設備・期間・製品を分けた検証と既存ルール比較を行う
- 人との分担：モデルは候補提示とし、現場確認、例外承認、停止判断の責任を明確にする
- 運用監視：データ遅延、分布変化、誤警報、未実行理由、改善効果を継続記録する
- 段階導入：一つのラインと一つの会議体で試し、業務効果と保守負荷を確認してから広げる
PoC の成功条件を「モデルが動く」ではなく、「意思決定時間または総損失が改善し、その効果を再測定できる」と定義することが実装の出発点です。
まとめ
No.081〜No.090では、品質管理、異常検知、需要予測、意思決定支援、製造業DIを、NumPy、SciPy、PyTorch、疎行列、NetworkXによる実装へ接続しました。行列は単なる計算形式ではなく、製品×ライン、観測×特徴量、行動×シナリオ、製品×部材という異なる業務を共通の構造で記述する道具です。
実務で価値を生むのは最も複雑なモデルではなく、変動と制約を可視化し、判断理由を共有し、結果から基準を更新できる仕組みです。本 notebook の数値は架空ですが、データ粒度、評価指標、制約、責任分界を自社の文脈へ置き換えることで、小さな検証の設計図として利用できます。
法人向けのご相談
数理工房では、製造業における品質分析、設備異常検知、需要予測、生産計画最適化、サプライチェーン可視化について、課題整理から PoC、運用設計、人材育成まで支援しています。データが十分に整っていない段階でも、意思決定と利用可能なデータを棚卸しし、着手範囲を一緒に設計できます。
📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。