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

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

2020-03-06

工程間のつながりから品質リスクと改善施策を読み解く

共分散・PCA・グラフ行列で学ぶ製造業の意思決定：100本ノック No.041〜No.050

製造ラインの品質問題は、単一設備だけで完結するとは限りません。温度・振動・電流などが連動し、上流の変動が複数工程を経て検査結果へ現れるためです。本記事では、架空の精密部品ラインを題材に、多変量データの構造を捉え、工程ネットワーク上の影響を評価し、改善施策を比較する一連の判断を扱います。

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

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

架空工場では、混合から梱包まで8工程で精密部品を生産しています。ロットごとにセンサー値と不良率を記録していますが、工程責任者には次の問いが残っています。

多数のセンサーのうち、同じ現象を重複して測っているものはどれか
品質変動を少数の軸で説明し、注意ロットを見つけられるか
上流異常が伝わる工程網で、どこから点検すべきか
改善策を導入したとき、停止日数や高リスク工程はどの程度減るか

目標は分析手法を並べることではなく、点検順位、対策候補、投資判断に使える比較材料へ変えることです。

現場でよくある状況

温度、振動、電流、圧力が同時に上がり、どれが主因か分からない
工程別KPIはあるが、工程間の搬送・再加工・品質伝播が評価に入っていない
異常発生工程だけを点検し、実際の起点である上流工程を見逃す
AIのPoCはできたが、学習データ不足や説明責任のため運用へ進めない
改善前後の平均値だけを比べ、停止リスクの分布やばらつきを見ていない

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

製造データには、変数間の相関と工程間の接続関係という二つの構造があります。表形式の集計だけでは、似たセンサーの重複や、迂回・再加工を含む影響経路を十分に表せません。

そこで前半では共分散行列とPCAでセンサー変動を整理します。後半では工程をノード、影響経路をエッジとするグラフを行列で表し、PageRank、マルコフ連鎖、グラフラプラシアン、ランダムウォーク、GCNなどへ展開します。最後に推薦とシミュレーションを用いて、分析結果を具体的な施策候補へ接続します。

今回扱うノックの全体像

No.	テーマ	製造業での判断
041	共分散行列	どのセンサーが連動して変化するか
042	PCA	品質変動を少数の軸へ要約できるか
043	PageRank	影響経路を踏まえて点検工程をどう順位付けするか
044	マルコフ連鎖	設備状態が将来どの比率へ移るか
045	グラフラプラシアン	工程網のまとまりと局所的な不整合はどこか
046	ランダムウォーク	異常の起点から訪れやすい工程はどこか
047	量子ランダムウォーク	干渉を伴う探索が古典法とどう異なるか
048	GCNへの応用	隣接工程の情報を含めてリスク表現を作る
049	推薦システムへの応用	設備ごとに未実施の改善策を候補化する
050	シミュレーションへの応用	改善案による停止リスク低減を比較する

Python 環境の準備

NumPyで行列計算、pandasで表形式の確認、Matplotlibで可視化を行います。外部データは使用しません。乱数生成器は np.random.default_rng(42) で固定し、同じ架空データとシミュレーション結果を再現できるようにします。

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

import platform
import sys

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

rng = np.random.default_rng(42)
pd.set_option("display.precision", 4)

print("Python:", sys.version.split()[0])
print("OS:", platform.platform())
print("NumPy:", np.__version__)
print("pandas:", pd.__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
Matplotlib: 3.11.0

架空データの作成

8工程を持つ精密部品ラインについて、正常寄りの240ロットと品質注意傾向を加えた80ロット、合計320ロットを生成します。センサー値は背後にある「負荷」「熱」「機械状態」の影響を受ける設定です。後半で使う工程ネットワークは、通常の順流に加え、研削から熱処理への再加工経路を持たせます。

注意傾向は説明用に人工的に付与しています。実務では、不良ラベルの定義時点、測定後の選別、再検査、品種構成の違いを整理し、未来の情報が説明変数へ混入するデータリーケージを防ぐ必要があります。

n_lots = 320
sensor_names = ["主軸電流", "軸受温度", "振動RMS", "油圧", "冷却水温", "加工時間"]
stations = ["混合", "成形", "粗加工", "熱処理", "研削", "洗浄", "検査", "梱包"]

# 潜在状態を混合して、相関を持つセンサーデータを生成する
latent = rng.normal(size=(n_lots, 3))
latent[:, 1] += 0.50 * latent[:, 0]
latent[:, 2] += 0.25 * latent[:, 0]
loading = np.array([
    [1.00, 0.20, 0.15], [0.25, 1.05, 0.25], [0.15, 0.25, 1.10],
    [0.75, 0.05, 0.20], [0.05, 0.90, 0.10], [0.80, 0.30, 0.35],
])
baseline = np.array([38.0, 66.0, 2.2, 5.6, 23.0, 41.0])
scale = np.array([4.5, 3.0, 0.42, 0.32, 1.7, 3.8])
sensor_values = baseline + (latent @ loading.T + rng.normal(0, 0.16, (n_lots, 6))) * scale

# 後半80ロットの一部に、熱・振動方向の注意傾向を付与
attention = np.zeros(n_lots, dtype=bool)
attention_rows = rng.choice(np.arange(240, n_lots), size=22, replace=False)
attention[attention_rows] = True
sensor_values[attention] += np.array([1.2, 5.5, 0.85, 0.05, 1.8, 2.0])

lot_df = pd.DataFrame(sensor_values, columns=sensor_names)
lot_df.insert(0, "ロット", [f"LOT-{i+1:03d}" for i in range(n_lots)])
lot_df["品質注意"] = attention
lot_df["不良率_pct"] = np.clip(
    0.8 + 0.22 * latent[:, 0] + 0.28 * latent[:, 2] + 2.2 * attention + rng.normal(0, 0.18, n_lots),
    0.05, None,
)

# 行iから列jへ影響が伝わる強さを表す有向隣接行列
A = np.zeros((len(stations), len(stations)))
edges = [
    (0, 1, 1.0), (1, 2, 0.9), (2, 3, 0.9), (3, 4, 1.0),
    (4, 5, 0.8), (5, 6, 0.9), (6, 7, 0.7), (4, 3, 0.25),
    (2, 5, 0.20), (3, 6, 0.35),
]
for i, j, weight in edges:
    A[i, j] = weight

display(lot_df.head().style.format({name: "{:.2f}" for name in sensor_names + ["不良率_pct"]}))
print("ロット数:", len(lot_df), " / 品質注意ロット数:", int(lot_df["品質注意"].sum()))

	ロット	主軸電流	軸受温度	振動RMS	油圧	冷却水温	加工時間	品質注意	不良率_pct
0	LOT-001	40.29	64.12	2.44	5.71	21.80	41.18	False	1.00
1	LOT-002	40.00	60.88	1.67	5.74	20.56	40.69	False	0.45
2	LOT-003	38.52	65.61	2.12	5.57	22.93	40.87	False	1.02
3	LOT-004	33.93	67.42	2.49	5.36	23.78	40.12	False	0.67
4	LOT-005	39.94	70.37	2.46	5.68	24.78	42.87	False	0.88

ロット数: 320  / 品質注意ロット数: 22

fig, axes = plt.subplots(1, 2, figsize=(11, 4.2))
axes[0].plot(lot_df.index, lot_df["軸受温度"], color="tab:red", linewidth=1)
axes[0].axvline(239.5, color="black", linestyle="--", label="評価期間開始")
axes[0].set_title("ロット別の軸受温度")
axes[0].set_xlabel("ロット順")
axes[0].set_ylabel("軸受温度（℃）")
axes[0].grid(True, alpha=0.3)
axes[0].legend()

im = axes[1].imshow(A, cmap="Blues", vmin=0, vmax=1)
axes[1].set_title("工程間の影響を表す有向隣接行列")
axes[1].set_xlabel("影響先工程")
axes[1].set_ylabel("影響元工程")
axes[1].set_xticks(range(len(stations)), stations, rotation=45, ha="right")
axes[1].set_yticks(range(len(stations)), stations)
axes[1].grid(True, color="white", linewidth=0.3, alpha=0.5)
fig.colorbar(im, ax=axes[1], label="影響の重み")
plt.tight_layout()
plt.show()

svg

No.041：共分散行列

実務での意味

共分散行列は、センサーが単独でどれだけばらつくかと、どの組み合わせが連動するかを一つの表にまとめます。連動するセンサー群の把握、重複計測の見直し、多変量異常検知の基礎になります。

分析・モデル化の考え方

変数 $j,k$ の標本共分散は、

s_{jk}=\frac{1}{n-1}\sum_{i=1}^{n}(x_{ij}-\bar{x}_j)(x_{ik}-\bar{x}_k)

です。対角要素は分散、非対角要素は二変数の共変動を表します。ただし共分散は単位の影響を受けるため、強さを比較するときは標準化した共分散、すなわち相関行列も確認します。相関は因果を意味せず、共通の運転条件が両方を動かす場合があります。

Pythonで確認する

train_sensors = lot_df.loc[:239, sensor_names]
cov_df = train_sensors.cov()
corr_df = train_sensors.corr()
display(cov_df.style.format("{:.3f}").background_gradient(cmap="Blues"))

fig, ax = plt.subplots(figsize=(7.4, 5.4))
im = ax.imshow(corr_df, cmap="coolwarm", vmin=-1, vmax=1)
ax.set_title("正常基準期間のセンサー相関行列")
ax.set_xlabel("センサー")
ax.set_ylabel("センサー")
ax.set_xticks(range(len(sensor_names)), sensor_names, rotation=35, ha="right")
ax.set_yticks(range(len(sensor_names)), sensor_names)
ax.grid(True, color="white", linewidth=0.4, alpha=0.5)
for i in range(len(sensor_names)):
    for j in range(len(sensor_names)):
        ax.text(j, i, f"{corr_df.iloc[i, j]:.2f}", ha="center", va="center", fontsize=8)
fig.colorbar(im, ax=ax, label="相関係数")
plt.tight_layout()
plt.show()

	主軸電流	軸受温度	振動RMS	油圧	冷却水温	加工時間
主軸電流	28.781	17.287	1.444	1.445	6.526	22.949
軸受温度	17.287	17.562	1.217	0.827	7.543	15.373
振動RMS	1.444	1.217	0.246	0.086	0.464	1.526
油圧	1.445	0.827	0.086	0.079	0.300	1.193
冷却水温	6.526	7.543	0.464	0.300	3.434	5.916
加工時間	22.949	15.373	1.526	1.193	5.916	19.812

svg

結果の読み取り

主軸電流・油圧・加工時間、軸受温度・冷却水温などに正の相関が見られます。したがって、すべてを独立した警報として扱うと同じ現象を重複通知する可能性があります。一方で、相関の高いセンサーを即座に削除するのも危険です。故障時だけ関係が崩れる場合があるため、正常時の相関、異常時の残差、センサーの冗長性という三つの役割を分けて評価します。

No.042：PCA（主成分分析）

実務での意味

PCAは、相関する複数センサーを、互いに直交する少数の合成指標へ要約します。多数の管理図を俯瞰図へまとめ、通常運転から外れるロットを発見する入口になります。

分析・モデル化の考え方

標準化データの共分散行列 $C$ を固有値分解し、固有値の大きい固有ベクトルから主成分軸を作ります。第 $k$ 主成分スコアは

\boldsymbol{t}_k=Z\boldsymbol{v}_k

です。寄与率が高い軸は通常変動をよく説明しますが、異常が小さい固有値方向へ現れることもあります。そのためスコア図だけでなく、再構成残差も運用上は確認します。

Pythonで確認する

train_mean = train_sensors.mean()
train_std = train_sensors.std(ddof=1)
Z = (lot_df[sensor_names] - train_mean) / train_std
C = np.cov(Z.iloc[:240], rowvar=False)
eigvals, eigvecs = np.linalg.eigh(C)
order = np.argsort(eigvals)[::-1]
eigvals, eigvecs = eigvals[order], eigvecs[:, order]
scores = Z.to_numpy() @ eigvecs
pca_summary = pd.DataFrame({
    "固有値": eigvals,
    "寄与率": eigvals / eigvals.sum(),
    "累積寄与率": np.cumsum(eigvals / eigvals.sum()),
}, index=[f"PC{i}" for i in range(1, 7)])
display(pca_summary.style.format({"固有値": "{:.3f}", "寄与率": "{:.1%}", "累積寄与率": "{:.1%}"}))

fig, axes = plt.subplots(1, 2, figsize=(11, 4.3))
axes[0].bar(range(1, 7), pca_summary["寄与率"] * 100, color="steelblue")
axes[0].plot(range(1, 7), pca_summary["累積寄与率"] * 100, marker="o", color="tab:orange")
axes[0].set_title("PCAの寄与率と累積寄与率")
axes[0].set_xlabel("主成分")
axes[0].set_ylabel("寄与率（%）")
axes[0].grid(True, axis="y", alpha=0.3)

colors = np.where(lot_df["品質注意"], "tab:red", "steelblue")
axes[1].scatter(scores[:, 0], scores[:, 1], c=colors, alpha=0.65, s=24)
axes[1].set_title("第1・第2主成分によるロット俯瞰")
axes[1].set_xlabel("第1主成分スコア")
axes[1].set_ylabel("第2主成分スコア")
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

	固有値	寄与率	累積寄与率
PC1	4.715	78.6%	78.6%
PC2	0.684	11.4%	90.0%
PC3	0.541	9.0%	99.0%
PC4	0.028	0.5%	99.5%
PC5	0.018	0.3%	99.8%
PC6	0.013	0.2%	100.0%

svg

結果の読み取り

先頭の少数主成分がセンサー変動の大部分を説明し、品質注意ロットの多くが通常ロット群から離れて現れます。これは「注意ロットを自動確定できた」という意味ではなく、複数センサーをまとめたレビュー候補を作れたという意味です。実務では主成分負荷量から寄与センサーを確認し、品種・設備・季節別にも同じ分離が保たれるか検証します。

No.043：PageRank

実務での意味

工程間に複数の影響経路があると、単純な接続本数だけでは重要工程を決められません。PageRankを使うと、重要な工程から影響を受ける工程を高く評価し、品質影響が集まりやすい点検候補を順位付けできます。

分析・モデル化の考え方

出次数で正規化した遷移行列を $P$ 、ダンピング係数を $\alpha$ とすると、PageRankベクトル $\boldsymbol{r}$ は

\boldsymbol{r}=\alpha P^\mathsf{T}\boldsymbol{r}+(1-\alpha)\frac{\boldsymbol{1}}{n}

を満たします。反復で定常解を求めます。ここでの順位は隣接行列の向きと重みの定義に依存するため、「物の流れ」「不良の伝播」「情報参照」のどれを表すかを先に固定します。

Pythonで確認する

row_sum = A.sum(axis=1, keepdims=True)
P_graph = np.divide(A, row_sum, out=np.full_like(A, 1 / len(stations)), where=row_sum > 0)
alpha = 0.85
rank = np.full(len(stations), 1 / len(stations))
for _ in range(200):
    new_rank = alpha * P_graph.T @ rank + (1 - alpha) / len(stations)
    if np.linalg.norm(new_rank - rank, ord=1) < 1e-12:
        break
    rank = new_rank

pagerank_df = pd.DataFrame({"工程": stations, "PageRank": rank}).sort_values("PageRank", ascending=False)
display(pagerank_df.style.format({"PageRank": "{:.4f}"}))

fig, ax = plt.subplots(figsize=(7.8, 4.1))
ax.bar(pagerank_df["工程"], pagerank_df["PageRank"], color="teal")
ax.set_title("品質影響ネットワークのPageRank")
ax.set_xlabel("工程")
ax.set_ylabel("PageRankスコア")
ax.grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()

	工程	PageRank
7	梱包	0.1983
6	検査	0.1864
5	洗浄	0.1371
3	熱処理	0.1365
4	研削	0.1258
2	粗加工	0.1024
1	成形	0.0737
0	混合	0.0398

svg

結果の読み取り

検査、熱処理、研削など、複数経路の影響が集まる工程が上位になります。点検人員が限られる場合の一次順位として利用できますが、PageRankが高い工程が原因工程とは限りません。影響の集約点である可能性もあるため、逆向きグラフによる起点探索、故障頻度、停止損失、検出可能性を合わせて保全優先度を決めます。

No.044：マルコフ連鎖

実務での意味

設備状態を「正常・注意・停止」のような段階に分け、状態遷移確率を推定すると、将来の停止比率、点検負荷、予備品需要を見積もれます。

分析・モデル化の考え方

現在の状態分布を行ベクトル $\boldsymbol{p}_t$ 、遷移行列を $P$ とすると、

\boldsymbol{p}_{t+k}=\boldsymbol{p}_tP^k

です。各行の和は1で、要素は非負である必要があります。通常の一次マルコフ連鎖は「次状態が現在状態だけに依存する」と仮定します。劣化履歴、累積稼働時間、保全内容が効く場合は状態を増やすか、別モデルを検討します。

Pythonで確認する

state_names = ["正常", "注意", "停止"]
P_state = np.array([
    [0.90, 0.09, 0.01],
    [0.35, 0.55, 0.10],
    [0.70, 0.20, 0.10],
])
p0 = np.array([0.92, 0.07, 0.01])
days = np.arange(0, 31)
state_path = np.array([p0 @ np.linalg.matrix_power(P_state, int(day)) for day in days])

evals, evecs = np.linalg.eig(P_state.T)
stationary = np.real(evecs[:, np.argmin(np.abs(evals - 1))])
stationary /= stationary.sum()
display(pd.DataFrame(P_state, index=state_names, columns=state_names).style.format("{:.1%}"))
display(pd.DataFrame({"30日後": state_path[-1], "定常分布": stationary}, index=state_names).style.format("{:.2%}"))

fig, ax = plt.subplots(figsize=(7.6, 4.2))
for i, state in enumerate(state_names):
    ax.plot(days, state_path[:, i] * 100, label=state)
ax.set_title("設備状態分布の30日予測")
ax.set_xlabel("経過日数")
ax.set_ylabel("設備構成比（%）")
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	正常	注意	停止
正常	90.0%	9.0%	1.0%
注意	35.0%	55.0%	10.0%
停止	70.0%	20.0%	10.0%

	30日後	定常分布
正常	79.96%	79.96%
注意	17.24%	17.24%
停止	2.80%	2.80%

svg

結果の読み取り

初期時点では正常設備が多くても、時間とともに遷移行列が定める定常比率へ近づきます。30日後の注意・停止比率は、点検要員や代替設備の必要量を考える基礎値になります。実務では遷移確率の推定件数と信頼区間を示し、保全方針や季節が変わった期間を一つに混ぜないようにします。

No.045：グラフラプラシアン

実務での意味

グラフラプラシアンは、工程ネットワークのつながり方を表し、工程群の分割、隣接工程間のKPIの滑らかさ、局所的な不整合の検出に使えます。

分析・モデル化の考え方

無向隣接行列を $W$ 、次数行列を $D$ とすると、ラプラシアンは

L=D-W

です。任意の工程KPIベクトル $\boldsymbol{x}$ に対し、 $\boldsymbol{x}^\mathsf{T}L\boldsymbol{x}$ は接続された工程間の差が大きいほど増えます。また、2番目に小さい固有値に対応するFiedlerベクトルは、工程網の分割候補を与えます。

Pythonで確認する

W = np.maximum(A, A.T)
D = np.diag(W.sum(axis=1))
L = D - W
lap_evals, lap_evecs = np.linalg.eigh(L)
fiedler = lap_evecs[:, 1]
cluster = np.where(fiedler >= 0, "工程群A", "工程群B")
lap_df = pd.DataFrame({"工程": stations, "Fiedler値": fiedler, "分割候補": cluster})
display(lap_df.style.format({"Fiedler値": "{:+.3f}"}))
print("小さい方から4つの固有値:", np.round(lap_evals[:4], 4))

fig, axes = plt.subplots(1, 2, figsize=(10.8, 4.1))
axes[0].plot(range(1, len(stations) + 1), lap_evals, marker="o")
axes[0].set_title("グラフラプラシアンの固有値")
axes[0].set_xlabel("固有値の順番")
axes[0].set_ylabel("固有値")
axes[0].grid(True, alpha=0.3)
axes[1].bar(stations, fiedler, color=np.where(fiedler >= 0, "steelblue", "tab:orange"))
axes[1].set_title("Fiedlerベクトルによる工程網の分割候補")
axes[1].set_xlabel("工程")
axes[1].set_ylabel("Fiedlerベクトルの要素")
axes[1].tick_params(axis="x", rotation=30)
axes[1].grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()

	工程	Fiedler値	分割候補
0	混合	+0.590	工程群A
1	成形	+0.452	工程群A
2	粗加工	+0.182	工程群A
3	熱処理	-0.046	工程群B
4	研削	-0.143	工程群B
5	洗浄	-0.222	工程群B
6	検査	-0.325	工程群B
7	梱包	-0.487	工程群B

小さい方から4つの固有値: [0.     0.2336 0.601  1.3984]

svg

結果の読み取り

Fiedlerベクトルの符号から、上流側と下流側を中心とする工程群の分割候補が得られます。これは品質会議や監視モデルを工程群単位で設計する手掛かりになります。ただし、分割は数学的な接続構造だけに基づくため、設備責任、物理距離、検査能力、ロット追跡可能性を踏まえて運用単位を決めます。

No.046：ランダムウォーク

実務での意味

異常の影響が接続先へ確率的に伝わると考えると、ランダムウォークで「起点から数ステップ後にどの工程へ到達しやすいか」を評価できます。追加検査の配置や追跡範囲を考える材料になります。

分析・モデル化の考え方

無向グラフの次数で正規化した遷移行列 $P=D^{-1}W$ を用い、初期分布 $\boldsymbol{p}_0$ から

\boldsymbol{p}_k=\boldsymbol{p}_0P^k

を計算します。ここでは接続重みを遷移しやすさとして扱います。現実の不良伝播率を表すには、トレーサビリティデータや工程実験から重みを推定する必要があります。

Pythonで確認する

P_walk = np.divide(W, W.sum(axis=1, keepdims=True), out=np.zeros_like(W), where=W.sum(axis=1, keepdims=True) > 0)
start_idx = stations.index("熱処理")
p_start = np.eye(len(stations))[start_idx]
walk_steps = np.arange(0, 9)
walk_dist = np.array([p_start @ np.linalg.matrix_power(P_walk, int(k)) for k in walk_steps])
display(pd.DataFrame(walk_dist[[0, 2, 4, 8]], index=["0歩", "2歩", "4歩", "8歩"], columns=stations).style.format("{:.1%}"))

fig, ax = plt.subplots(figsize=(8.5, 4.6))
im = ax.imshow(walk_dist.T, aspect="auto", cmap="YlOrRd", vmin=0)
ax.set_title("熱処理を起点とするランダムウォーク")
ax.set_xlabel("ステップ数")
ax.set_ylabel("工程")
ax.set_xticks(range(len(walk_steps)), walk_steps)
ax.set_yticks(range(len(stations)), stations)
ax.grid(True, color="white", linewidth=0.3, alpha=0.5)
fig.colorbar(im, ax=ax, label="到達確率")
plt.tight_layout()
plt.show()

	混合	成形	粗加工	熱処理	研削	洗浄	検査	梱包
0歩	0.0%	0.0%	0.0%	100.0%	0.0%	0.0%	0.0%	0.0%
2歩	0.0%	18.0%	0.0%	45.5%	0.0%	30.9%	0.0%	5.6%
4歩	0.0%	23.0%	0.0%	36.9%	0.0%	30.4%	0.0%	9.8%
8歩	0.0%	26.1%	0.0%	33.8%	0.0%	29.2%	0.0%	10.8%

svg

結果の読み取り

熱処理から始めた影響は、研削、検査、再加工経路などへ確率的に広がります。ステップ数によって重点工程が変わるため、発生直後の隔離と、時間が経過した後の追跡では見る範囲を変える必要があります。実務ではロットの実際の経路、滞留時間、分岐率を使い、有向グラフとして再設計します。

No.047：量子ランダムウォーク

実務での意味

量子ランダムウォークは、確率ではなく複素振幅を伝播させ、経路間の干渉を表します。量子探索・量子アルゴリズムの基礎概念であり、複雑なグラフ探索を考える研究テーマです。本ノックでは、製造現場へ即導入する手法ではなく、古典的ランダムウォークとの違いを小規模行列で確認します。

分析・モデル化の考え方

連続時間量子ウォークでは、グラフラプラシアン $L$ をHamiltonianとして

|\psi(t)\rangle=\exp(-itL)|\psi(0)\rangle

と発展させ、工程 $j$ の観測確率を $|\psi_j(t)|^2$ とします。古典ウォークと異なり確率分布へ単調収束せず、干渉による振動が生じます。以下は理想的な数学モデルで、量子優位や実機効果を示すものではありません。

Pythonで確認する

times = np.linspace(0, 12, 121)
psi0 = np.zeros(len(stations), dtype=complex)
psi0[start_idx] = 1.0

# L = QΛQ^T を使い、行列指数関数を固有値ごとに計算する
quantum_prob = []
for t in times:
    phase = np.exp(-1j * lap_evals * t)
    psi_t = lap_evecs @ (phase * (lap_evecs.T @ psi0))
    quantum_prob.append(np.abs(psi_t) ** 2)
quantum_prob = np.array(quantum_prob)

print("全時刻における確率和の最大誤差:", f"{np.max(np.abs(quantum_prob.sum(axis=1) - 1)):.3e}")
display(pd.DataFrame(quantum_prob[[0, 25, 60, 120]], index=["t=0", "t=2.5", "t=6", "t=12"], columns=stations).style.format("{:.1%}"))

fig, ax = plt.subplots(figsize=(8.5, 4.6))
im = ax.imshow(quantum_prob.T, aspect="auto", origin="lower", cmap="viridis", extent=[times.min(), times.max(), -0.5, 7.5])
ax.set_title("連続時間量子ランダムウォークの観測確率")
ax.set_xlabel("時間パラメータ t")
ax.set_ylabel("工程")
ax.set_yticks(range(len(stations)), stations)
ax.grid(True, color="white", linewidth=0.3, alpha=0.35)
fig.colorbar(im, ax=ax, label="観測確率")
plt.tight_layout()
plt.show()

全時刻における確率和の最大誤差: 6.661e-16

	混合	成形	粗加工	熱処理	研削	洗浄	検査	梱包
t=0	0.0%	0.0%	0.0%	100.0%	0.0%	0.0%	0.0%	0.0%
t=2.5	17.8%	8.0%	5.6%	0.4%	8.0%	4.5%	44.2%	11.6%
t=6	4.5%	7.6%	18.9%	18.3%	7.8%	14.8%	14.5%	13.5%
t=12	4.7%	4.9%	17.9%	0.2%	4.5%	33.5%	28.6%	5.8%

svg

結果の読み取り

確率の総和は数値誤差の範囲で1を保ちつつ、工程ごとの観測確率は時間とともに振動します。No.046の古典ウォークのような単調な拡散とは異なり、経路の干渉が分布を変える点が特徴です。現時点の実務判断では古典法を基準とし、量子法は問題規模、データ入力コスト、古典計算との公平な比較、実機ノイズまで含む検証課題として扱うのが妥当です。

No.048：GCNへの応用

実務での意味

GCN（Graph Convolutional Network）は、自工程の特徴だけでなく隣接工程の情報を集約して学習します。上流・下流の状態が品質へ影響するラインで、工程リスク予測や異常箇所推定に利用できます。

分析・モデル化の考え方

自己ループを加えた隣接行列を $\tilde{A}=W+I$ 、次数行列を $\tilde{D}$ とすると、1層のGCNは

H=\mathrm{ReLU}(\tilde{D}^{-1/2}\tilde{A}\tilde{D}^{-1/2}X\Theta)

と書けます。 $X$ は工程特徴、 $\Theta$ は学習する重みです。ここでは情報集約の仕組みを示すため固定重みを使います。実運用では時系列分割した教師データで学習・評価し、未知設備への一般化を確認します。

Pythonで確認する

# 工程ごとの架空特徴：温度偏差、振動偏差、直近不良率
X_node = np.array([
    [0.2, 0.1, 0.4], [0.4, 0.2, 0.6], [0.8, 0.7, 1.1], [1.7, 1.0, 1.5],
    [1.1, 1.8, 1.8], [0.5, 0.4, 0.8], [0.3, 0.2, 1.4], [0.2, 0.1, 0.5],
])
A_tilde = W + np.eye(len(stations))
d_inv_sqrt = np.diag(1 / np.sqrt(A_tilde.sum(axis=1)))
A_norm = d_inv_sqrt @ A_tilde @ d_inv_sqrt
theta = np.array([[0.8, -0.2], [0.5, 0.7], [0.6, 0.4]])
H = np.maximum(0, A_norm @ X_node @ theta)
gcn_risk = H @ np.array([0.65, 0.35])

gcn_df = pd.DataFrame({
    "工程": stations,
    "自工程特徴の単純合計": X_node.sum(axis=1),
    "GCN集約リスク": gcn_risk,
}).sort_values("GCN集約リスク", ascending=False)
display(gcn_df.style.format({"自工程特徴の単純合計": "{:.2f}", "GCN集約リスク": "{:.3f}"}))

fig, ax = plt.subplots(figsize=(7.8, 4.1))
ax.bar(gcn_df["工程"], gcn_df["GCN集約リスク"], color="slateblue")
ax.set_title("隣接工程情報を集約したGCNリスク表現")
ax.set_xlabel("工程")
ax.set_ylabel("GCN集約リスク（例示値）")
ax.grid(True, axis="y", alpha=0.3)
plt.tight_layout()
plt.show()

	工程	自工程特徴の単純合計	GCN集約リスク
3	熱処理	4.20	1.975
4	研削	4.70	1.836
5	洗浄	1.70	1.393
2	粗加工	2.60	1.307
6	検査	1.90	0.975
1	成形	1.20	0.770
7	梱包	0.80	0.552
0	混合	0.70	0.434

svg

結果の読み取り

熱処理や研削の高い特徴が隣接工程へ反映され、自工程の単純集計とは異なる順位が得られます。これがGCNの情報集約の基本です。ただし、今回の値は未学習の例示スコアであり、予測性能を主張するものではありません。導入時は比較対象としてロジスティック回帰や木モデルを置き、GCNによる改善が工程グラフの保守コストに見合うか確認します。

No.049：推薦システムへの応用

実務での意味

設備ごとに試した改善策と効果評価を行列へ整理すると、未実施の組み合わせを補完し、「次に小規模検証する施策」の候補を推薦できます。熟練者の知見を置き換えるのではなく、見落としを減らす用途です。

分析・モデル化の考え方

設備×施策の評価行列 $R$ には未評価要素があります。欠損を列平均で仮補完し、中心化後の切断SVD

R\approx U_r\Sigma_rV_r^\mathsf{T}

で低次元の類似構造を作ります。実務では欠損をゼロ評価と混同せず、観測済み要素だけを損失関数へ入れる行列分解や、設備仕様・費用・安全制約を加えたランキングを使います。

Pythonで確認する

machine_names = ["旋盤A", "旋盤B", "研削A", "研削B", "熱処理炉A", "熱処理炉B"]
action_names = ["給油周期短縮", "冷却条件変更", "工具交換前倒し", "芯出し調整", "温度監視強化", "振動監視強化", "搬送速度調整"]
ratings = np.array([
    [4.5, 3.0, 4.8, np.nan, 2.8, 4.0, np.nan],
    [4.2, 3.2, 4.6, 4.0, np.nan, 4.3, np.nan],
    [3.0, 3.8, 4.2, 4.7, np.nan, 4.9, 2.5],
    [np.nan, 3.6, 4.0, 4.5, 3.2, 4.8, 2.8],
    [2.5, 4.8, np.nan, 2.8, 4.9, 3.4, 3.6],
    [2.7, 4.6, 3.0, np.nan, 4.7, 3.5, 3.8],
])
observed = ~np.isnan(ratings)
item_mean = np.nanmean(ratings, axis=0)
filled = np.where(observed, ratings, item_mean)
centered = filled - item_mean
U, s, Vt = np.linalg.svd(centered, full_matrices=False)
rank = 2
estimated = (U[:, :rank] * s[:rank]) @ Vt[:rank] + item_mean

recommendations = []
for i, machine in enumerate(machine_names):
    candidates = np.where(~observed[i])[0]
    if len(candidates):
        best = candidates[np.argmax(estimated[i, candidates])]
        recommendations.append((machine, action_names[best], estimated[i, best]))
recommend_df = pd.DataFrame(recommendations, columns=["設備", "次の検証候補", "推定評価"])
display(pd.DataFrame(ratings, index=machine_names, columns=action_names).style.format("{:.1f}", na_rep="未実施"))
display(recommend_df.style.format({"推定評価": "{:.2f}"}))

fig, ax = plt.subplots(figsize=(8.8, 4.8))
im = ax.imshow(estimated, cmap="YlGn", vmin=2, vmax=5)
ax.set_title("低ランク行列分解による設備×施策の推定評価")
ax.set_xlabel("改善施策")
ax.set_ylabel("設備")
ax.set_xticks(range(len(action_names)), action_names, rotation=40, ha="right")
ax.set_yticks(range(len(machine_names)), machine_names)
ax.grid(True, color="white", linewidth=0.3, alpha=0.5)
fig.colorbar(im, ax=ax, label="推定評価")
plt.tight_layout()
plt.show()

	給油周期短縮	冷却条件変更	工具交換前倒し	芯出し調整	温度監視強化	振動監視強化	搬送速度調整
旋盤A	4.5	3.0	4.8	未実施	2.8	4.0	未実施
旋盤B	4.2	3.2	4.6	4.0	未実施	4.3	未実施
研削A	3.0	3.8	4.2	4.7	未実施	4.9	2.5
研削B	未実施	3.6	4.0	4.5	3.2	4.8	2.8
熱処理炉A	2.5	4.8	未実施	2.8	4.9	3.4	3.6
熱処理炉B	2.7	4.6	3.0	未実施	4.7	3.5	3.8

	設備	次の検証候補	推定評価
0	旋盤A	芯出し調整	3.93
1	旋盤B	温度監視強化	3.45
2	研削A	温度監視強化	3.70
3	研削B	給油周期短縮	3.46
4	熱処理炉A	工具交換前倒し	3.74
5	熱処理炉B	芯出し調整	3.67

svg

結果の読み取り

各設備について、未実施の施策から推定評価が高い候補を一つずつ抽出できました。これは実施決定ではなく、次の検証順を決める材料です。安全要件、停止時間、費用、設備メーカーの保証条件を満たさない施策は推薦対象から除外し、推定評価の高いものもA/Bテストや段階導入で実効果を確認します。

No.050：シミュレーションへの応用

実務での意味

改善策の価値は、平均停止率だけでなく「30日で何日止まり得るか」という分布で示すと判断しやすくなります。行列で定義した状態遷移をモンテカルロシミュレーションへ展開し、現行運用と予防保全強化案を比較します。

分析・モデル化の考え方

No.044の遷移行列を現行案とし、注意から正常への回復を増やし、注意から停止への遷移を減らした改善案を作ります。各シナリオで多数の30日経路を発生させ、停止日数の平均、95パーセンタイル、一度以上停止する確率を比較します。

期待値だけでは稀な長期停止を表せません。一方、シミュレーション結果は入力した遷移確率の範囲を超えて正しくならないため、感度分析と実績による更新が必要です。

Pythonで確認する

P_improved = np.array([
    [0.92, 0.075, 0.005],
    [0.48, 0.47, 0.05],
    [0.78, 0.17, 0.05],
])

def simulate_stop_days(P, n_runs=10000, horizon=30):
    states = rng.choice(3, size=n_runs, p=p0)
    stop_days = np.zeros(n_runs, dtype=int)
    for _ in range(horizon):
        u = rng.random(n_runs)
        cumulative = np.cumsum(P[states], axis=1)
        states = (u[:, None] > cumulative).sum(axis=1)
        stop_days += states == 2
    return stop_days

stop_current = simulate_stop_days(P_state)
stop_improved = simulate_stop_days(P_improved)
simulation_df = pd.DataFrame({
    "シナリオ": ["現行運用", "予防保全強化"],
    "平均停止日数": [stop_current.mean(), stop_improved.mean()],
    "95%点": [np.quantile(stop_current, 0.95), np.quantile(stop_improved, 0.95)],
    "1日以上停止する確率": [(stop_current >= 1).mean(), (stop_improved >= 1).mean()],
})
display(simulation_df.style.format({"平均停止日数": "{:.2f}", "95%点": "{:.0f}", "1日以上停止する確率": "{:.1%}"}))

bins = np.arange(-0.5, max(stop_current.max(), stop_improved.max()) + 1.5)
fig, ax = plt.subplots(figsize=(7.8, 4.2))
ax.hist(stop_current, bins=bins, alpha=0.60, density=True, label="現行運用", color="tab:red")
ax.hist(stop_improved, bins=bins, alpha=0.60, density=True, label="予防保全強化", color="steelblue")
ax.set_title("30日間の停止日数シミュレーション")
ax.set_xlabel("30日間の停止日数")
ax.set_ylabel("相対頻度")
ax.grid(True, axis="y", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

	シナリオ	平均停止日数	95%点	1日以上停止する確率
0	現行運用	0.82	3	53.8%
1	予防保全強化	0.34	1	27.9%

svg

結果の読み取り

予防保全強化案では、平均停止日数と1日以上停止する確率の双方が低下します。改善費用と停止1日当たりの損失を組み合わせれば、期待便益や投資回収の試算へ進めます。ただし、改善後の遷移確率は仮定です。まず一部設備で試行し、実測した回復率・停止率で行列を更新してから全体展開を判断します。

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

表と工程網の二つの構造を見る：共分散・PCAはセンサー間、グラフ行列は工程間の関係を整理します。
順位は原因確定ではない：PageRankやランダムウォークは点検候補を絞りますが、物理原因は現場確認で確定します。
モデルの複雑さは比較で正当化する：GCNや量子法は、単純な統計・古典法を基準に追加価値を評価します。
推薦は検証順を作る：未実施施策の推定値は、制約確認と小規模試験を省略する根拠にはなりません。
点予測から分布へ進む：施策効果を平均だけでなく停止日数の分布で示すと、リスク許容度を含む判断ができます。

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

1. 工程・変数・時間の定義をそろえる

工程ID、設備ID、ロット、品種、時刻、再加工経路を一意に追跡できるようにします。センサー単位、集計窓、欠損処理、正常基準期間をデータ辞書へ残します。

2. グラフの意味と更新責任を決める

隣接行列が物の流れ、不良伝播、情報参照のどれを表すか明記します。工程変更、設備増設、搬送経路変更が起きたとき、誰がグラフと重みを更新するか決めます。

3. 単純な基準モデルと時系列で比較する

PCA、ルール、ロジスティック回帰などを基準に、複雑なモデルの精度、先行検知時間、誤警報、計算時間を時系列分割で比較します。品種・設備・期間別の性能差も確認します。

4. 現場KPIと意思決定手順へ接続する

スコア上位を誰が確認し、どの追加測定を行い、どの条件で停止・継続・経過観察とするかを定めます。停止時間、不要点検、仕損費、納期影響、確認工数で効果を評価します。

5. 不確実性と監査可能性を残す

遷移確率や推薦値には推定誤差があります。入力データ、コード、乱数seed、モデル版、判断結果を記録し、人が覆した判断も次回更新へ活かします。

まとめ

No.041〜No.050では、架空の精密部品ラインを題材に、共分散行列とPCAによる多変量センサーの把握、PageRank・マルコフ連鎖・グラフラプラシアン・ランダムウォークによる工程ネットワーク分析、量子ランダムウォークの概念、GCNの情報集約、行列分解による改善策推薦、状態遷移シミュレーションを確認しました。

行列の価値は、複雑な計算そのものではありません。センサー、工程、施策の関係を共通の形式で表し、現場が次に何を確認し、どの案を試すか決められることにあります。

法人向けのご相談

数理工房では、製造業の多変量センサー分析、工程ネットワーク分析、設備状態モデル、異常検知、改善施策シミュレーション、PoCから現場運用までをご支援しています。

「工程別データはあるが全体の品質影響を追えない」「異常検知の点検順位を現場へ説明したい」「改善施策の効果を導入前に比較したい」といった課題について、データと業務フローの棚卸しから小規模検証、運用設計までご相談いただけます。

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

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

工程間のつながりから品質リスクと改善施策を読み解く

共分散・PCA・グラフ行列で学ぶ製造業の意思決定：100本ノック No.041〜No.050

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

現場でよくある状況

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

今回扱うノックの全体像

Python 環境の準備

架空データの作成

No.041：共分散行列

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.042：PCA（主成分分析）

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.043：PageRank

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.044：マルコフ連鎖

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.045：グラフラプラシアン

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.046：ランダムウォーク

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.047：量子ランダムウォーク

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.048：GCNへの応用

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.049：推薦システムへの応用

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.050：シミュレーションへの応用

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

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

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

1. 工程・変数・時間の定義をそろえる

2. グラフの意味と更新責任を決める

3. 単純な基準モデルと時系列で比較する

4. 現場KPIと意思決定手順へ接続する

5. 不確実性と監査可能性を残す

まとめ

法人向けのご相談

Related Posts

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

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

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

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

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