工程ネットワークの安定性を数値監査する — 製造業の線形代数100本ノック No.051〜No.060

設備間の影響を行列で表したとき、「理論上は安定」だけでは現場判断に足りません。本稿では、4設備からなる架空の循環工程を題材に、半正定値行列から行列の安定性までを一つの監査手順として扱います。

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

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

切削、洗浄、熱処理、検査の負荷偏差が翌シフトへ持ち越されるラインを考えます。状態を $\mathbf{x}_t$、伝播行列を $A$ とすると、$\mathbf{x}_{t+1}=A\mathbf{x}_t+\boldsymbol{\varepsilon}_t$ です。

現場でよくある状況

• 工程別KPIはあるが、工程間の波及が評価されていない
• シミュレーション結果はあるが、パラメータ誤差への感度が不明
• 大規模設備データで固有値計算が遅い

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

安定性は行列全体の固有構造で決まります。非対称行列では長期安定でも短期増幅があり、測定誤差や丸め誤差も判断へ影響します。

今回扱うノックの全体像

No.051〜054で構造から固有値を評価し、No.055で数値信頼性を点検します。No.056〜058で計算法を比較し、No.059〜060で時間応答と安定判定へ統合します。

Python 環境の準備

NumPy、pandas、Matplotlib、SciPyを使います。乱数seedを固定し、表示の丸めと内部計算を分けます。

import sys, numpy as np, pandas as pd, matplotlib, matplotlib.pyplot as plt
from scipy.linalg import expm
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh
from IPython.display import display
rng=np.random.default_rng(51060)
np.set_printoptions(precision=4,suppress=True)
plt.rcParams.update({"figure.figsize":(7.2,4.2),"axes.grid":True})
print("Python",sys.version.split()[0],"NumPy",np.__version__,"pandas",pd.__version__,"Matplotlib",matplotlib.__version__)

Python 3.13.1 NumPy 2.5.1 pandas 3.0.3 Matplotlib 3.11.0

架空データの作成

4設備の翌シフトへの影響を非負行列 $A$ で表し、120シフトを生成します。対角成分は自己持越し、非対角成分は工程間波及です。実務では推定期間とモデル版を記録します。

equipment=["Cutting","Washing","HeatTreat","Inspection"]
A=np.array([[.56,.10,.03,.02],[.08,.50,.11,.03],[.04,.09,.61,.10],[.03,.04,.12,.48]])
states=np.zeros((120,4)); states[0]=[1.2,-.5,.8,.3]
for t in range(1,120): states[t]=A@states[t-1]+rng.normal(0,.12,4)
display(pd.DataFrame(A,index=equipment,columns=equipment).round(2))
plt.plot(states[:40]); plt.title("Standardized equipment deviations"); plt.xlabel("Shift"); plt.ylabel("Standardized deviation"); plt.legend(equipment,ncol=2); plt.grid(True); plt.tight_layout(); plt.show()

	Cutting	Washing	HeatTreat	Inspection
Cutting	0.56	0.10	0.03	0.02
Washing	0.08	0.50	0.11	0.03
HeatTreat	0.04	0.09	0.61	0.10
Inspection	0.03	0.04	0.12	0.48

No.051：半正定値行列

実務での意味

工程間の不均衡損失や共分散行列が負の評価を返さないための土台です。ただしゼロ固有値があると、評価しない偏差方向が残ります。

分析・モデル化の考え方

対称行列 $Q$ が任意の $\mathbf{x}$ に対して $\mathbf{x}^TQ\mathbf{x}\ge0$ を満たすことと、全固有値が非負であることは同値です。工程差を罰するLaplacian $L=D-W$ を使います。

Pythonで確認する

W=np.array([[0,1,.2,0],[1,0,.7,.1],[.2,.7,0,.9],[0,.1,.9,0]])
L=np.diag(W.sum(1))-W; ev=np.linalg.eigvalsh(L); loss=np.einsum("ij,jk,ik->i",states,L,states)
display(pd.DataFrame({"eigenvalue":ev}).round(6)); print("minimum loss",loss.min())
plt.hist(loss,bins=18); plt.title("Cross-equipment imbalance loss"); plt.xlabel(r"$x^T L x$"); plt.ylabel("Shifts"); plt.grid(True); plt.tight_layout(); plt.show()

	eigenvalue
0	-0.000000
1	0.725724
2	2.206091
3	2.868186

minimum loss 0.005229475856508912

結果の読み取り

固有値と損失は非負です。ゼロ固有値は全設備が同じだけずれる共通モードなので、工程間差のKPIだけでなく全体の中心ずれKPIを併用します。

No.052：Rayleigh商

実務での意味

候補となる偏差パターンがどれほど残りやすいかを一値で比較し、改善案の優先順位に使えます。

分析・モデル化の考え方

対称部 $S=(A+A^T)/2$ に対する $R(\mathbf{x})=\mathbf{x}^TS\mathbf{x}/\mathbf{x}^T\mathbf{x}$ は最小・最大固有値の間にあります。

Pythonで確認する

S=(A+A.T)/2; X=rng.normal(size=(2000,4)); X/=np.linalg.norm(X,axis=1,keepdims=True)
rq=np.einsum("ij,jk,ik->i",X,S,X); sev=np.linalg.eigvalsh(S)
print("sample",rq.min(),rq.max(),"theory",sev[0],sev[-1])
plt.hist(rq,bins=30); plt.axvline(sev[-1],color="red",ls="--"); plt.title("Rayleigh quotients"); plt.xlabel("Rayleigh quotient"); plt.ylabel("Count"); plt.grid(True); plt.tight_layout(); plt.show()

sample 0.4034253637601614 0.7513199375074551 theory 0.3994312035852489 0.7546030144172374

結果の読み取り

標本値は理論範囲内です。最大値へ近い工程組合せを重点監視できますが、非対称な $A$ の長期安定性は別途固有値で確認します。

No.053：Gershgorinの定理

実務での意味

係数表だけから固有値範囲を見積もる高速スクリーニングで、モデル更新時の日常監視に向きます。

分析・モデル化の考え方

各行の中心 $a_{ii}$、半径 $r_i=\sum_{j\ne i}|a_{ij}|$ の円盤内に全固有値が入ります。全円盤が単位円内なら安定の十分条件です。

Pythonで確認する

c=np.diag(A); r=np.sum(abs(A),1)-abs(c); g=pd.DataFrame({"center":c,"radius":r,"right":c+r},index=equipment); display(g.round(3))
th=np.linspace(0,2*np.pi,300)
for name,ci,ri in zip(equipment,c,r): plt.plot(ci+ri*np.cos(th),ri*np.sin(th),label=name)
e=np.linalg.eigvals(A); plt.scatter(e.real,e.imag,c="black",marker="x"); plt.title("Gershgorin disks and eigenvalues"); plt.xlabel("Real part"); plt.ylabel("Imaginary part"); plt.legend(ncol=2); plt.grid(True); plt.tight_layout(); plt.show()

	center	radius	right
Cutting	0.56	0.15	0.71
Washing	0.50	0.22	0.72
HeatTreat	0.61	0.23	0.84
Inspection	0.48	0.19	0.67

結果の読み取り

全円盤の右端が1未満なので安定性を安価に保証できます。円盤がはみ出す場合は不安定と断定せず、精密固有値計算へ回します。

No.054：Perron-Frobeniusの定理

実務での意味

非負の工程行列で、長期に残る共通変動の大きさと設備構成を示し、調査優先度を決めます。

分析・モデル化の考え方

正の行列ではスペクトル半径 $\rho(A)$ が正の単純固有値で、対応固有ベクトルは正成分を持ちます。和を1に正規化して寄与を読みます。

Pythonで確認する

vals,vecs=np.linalg.eig(A); ii=np.argmax(abs(vals)); rho=vals[ii].real; v=abs(vecs[:,ii].real); v/=v.sum()
pf=pd.DataFrame({"equipment":equipment,"share":v}).sort_values("share",ascending=False); display(pf.round(4)); print("spectral radius",rho)
plt.bar(pf.equipment,pf.share); plt.title("Dominant propagation mode"); plt.xlabel("Equipment"); plt.ylabel("Normalized share"); plt.xticks(rotation=15); plt.grid(True); plt.tight_layout(); plt.show()

	equipment	share
2	HeatTreat	0.3509
1	Washing	0.2395
3	Inspection	0.2104
0	Cutting	0.1993

spectral radius 0.754099720331183

結果の読み取り

支配固有値は1未満で、熱処理の寄与が大きい結果です。これは因果の断定ではなく、センサー精査や改善実験の優先順位です。

No.055：行列の条件数

実務での意味

測定や係数の小誤差が、原因寄与の逆算結果へ何倍程度増幅され得るかを示します。

分析・モデル化の考え方

2ノルム条件数は $\kappa_2(M)=\sigma_{max}/\sigma_{min}$ です。良条件と、列がほぼ重複する悪条件を比較します。

Pythonで確認する

Mg=np.array([[1,.2,.1],[.1,1.1,.2],[.2,.1,.9]]); Mb=np.array([[1,1.001,.1],[.5,.5004,.2],[.2,.2003,1]])
y=np.array([1,.7,.4]); d=np.array([1e-4,-1e-4,1e-4]); rows=[]
for name,M in [("well-conditioned",Mg),("ill-conditioned",Mb)]:
 z=np.linalg.solve(M,y); zp=np.linalg.solve(M,y+d); rows.append([name,np.linalg.cond(M),np.linalg.norm(zp-z)/np.linalg.norm(z)])
df=pd.DataFrame(rows,columns=["case","condition_number","relative_solution_change"]); display(df.round(6))
plt.bar(df.case,df.condition_number); plt.yscale("log"); plt.title("Condition numbers"); plt.xlabel("Case"); plt.ylabel("Condition number (log)"); plt.grid(True); plt.tight_layout(); plt.show()

	case	condition_number	relative_solution_change
0	well-conditioned	1.664609	0.000185
1	ill-conditioned	22771.393274	0.000957

結果の読み取り

悪条件行列では寄与推定が敏感です。センサー追加、変数統合、正則化を検討し、有効桁を制限します。条件数は実誤差ではなく最悪時の増幅しやすさです。

No.056：冪乗法

実務での意味

大規模モデルで最大固有値と支配モードだけを低コストに更新し、安定余裕を定期監視します。

分析・モデル化の考え方

$\mathbf{v}_{k+1}=A\mathbf{v}_k/\|A\mathbf{v}_k\|$ を反復し、Rayleigh商で固有値を推定します。収束速度は固有値の比に依存します。

Pythonで確認する

v=np.ones(4)/2; hist=[]
for k in range(25):
 w=A@v; v=w/np.linalg.norm(w); hist.append(float(v@(A@v)/(v@v)))
df=pd.DataFrame({"iteration":range(1,26),"estimate":hist,"error":abs(np.array(hist)-rho)}); display(df.iloc[[0,1,2,4,9,24]].round(7))
plt.semilogy(df.iteration,df.error,marker="o"); plt.title("Power-method convergence"); plt.xlabel("Iteration"); plt.ylabel("Absolute error"); plt.grid(True); plt.tight_layout(); plt.show()

	iteration	estimate	error
0	1	0.745741	8.359300e-03
1	2	0.750208	3.892100e-03
2	3	0.752186	1.913300e-03
4	5	0.753591	5.087000e-04
9	10	0.754081	1.890000e-05
24	25	0.754100	1.000000e-07

結果の読み取り

推定値は支配固有値へ収束します。本番では回数だけでなく残差 $\|A\mathbf{v}-\lambda\mathbf{v}\|$ を停止条件にします。

No.057：Lanczos法

実務での意味

多数設備の対称疎行列から端の固有値だけを少ないメモリで求め、更新頻度を上げます。

分析・モデル化の考え方

Krylov部分空間へ射影して三重対角化します。eigsh を使い400設備の上位3固有値を求め、既知の厳密式と照合します。

Pythonで確認する

n=400; LS=diags([.12*np.ones(n-1),.62*np.ones(n),.12*np.ones(n-1)],[-1,0,1],format="csr")
lv=eigsh(LS,k=3,which="LA",return_eigenvectors=False)[::-1]; exact=.62+.24*np.cos(np.arange(1,4)*np.pi/(n+1))
df=pd.DataFrame({"rank":[1,2,3],"Lanczos":lv,"exact":exact,"error":abs(lv-exact)}); display(df.round(10))
plt.plot(df["rank"],df.Lanczos,"o-",label="Lanczos"); plt.plot(df["rank"],df.exact,"x--",label="Exact"); plt.title("Largest sparse-matrix eigenvalues"); plt.xlabel("Rank"); plt.ylabel("Eigenvalue"); plt.legend(); plt.grid(True); plt.tight_layout(); plt.show()

	rank	Lanczos	exact	error
0	1	0.859993	0.859993	0.0
1	2	0.859971	0.859971	0.0
2	3	0.859934	0.859934	0.0

結果の読み取り

密行列へ変換せず厳密値と一致します。非対称行列にはArnoldi法を選び、許容誤差、最大反復、残差を記録します。

No.058：QR法

実務での意味

全固有値を棚卸しし、部分的な反復法を検証する基準にできます。

分析・モデル化の考え方

$A_k=Q_kR_k$、$A_{k+1}=R_kQ_k$ を反復します。相似変換で固有値を保ち、対称行列では対角へ収束します。

Pythonで確認する

Ak=S.copy(); off=[]
for k in range(60):
 Q,R=np.linalg.qr(Ak); Ak=R@Q; off.append(np.linalg.norm(Ak-np.diag(np.diag(Ak)),"fro"))
qv=np.sort(np.diag(Ak))[::-1]; tv=np.sort(np.linalg.eigvalsh(S))[::-1]; display(pd.DataFrame({"QR":qv,"library":tv,"error":abs(qv-tv)}).round(8))
plt.semilogy(range(1,61),off); plt.title("QR off-diagonal convergence"); plt.xlabel("Iteration"); plt.ylabel("Off-diagonal norm"); plt.grid(True); plt.tight_layout(); plt.show()

	QR	library	error
0	0.754603	0.754603	0.000000
1	0.560052	0.560052	0.000000
2	0.435912	0.435914	0.000002
3	0.399433	0.399431	0.000002

結果の読み取り

教育用の無シフト実装でもライブラリ結果へ一致します。本番では高速化・検証済みの eigvalsh などを使い、独自実装を避けます。

No.059：行列関数

実務での意味

連続時間の温度・濃度・振動モデルから任意時間後の状態と復帰時間を計算します。

分析・モデル化の考え方

$\dot{\mathbf{x}}=F\mathbf{x}$ の解は $\mathbf{x}(t)=e^{Ft}\mathbf{x}(0)$ です。行列指数は要素別指数ではありません。

Pythonで確認する

F=np.array([[-.45,.16,0],[.08,-.34,.10],[0,.12,-.28]]); x0=np.array([2,.2,1]); tt=np.linspace(0,16,65); tr=np.array([expm(F*t)@x0 for t in tt])
print("eigenvalues",np.linalg.eigvals(F)); display(pd.DataFrame(tr[[0,16,32,64]],index=tt[[0,16,32,64]],columns=["Zone A","Zone B","Zone C"]).round(4))
plt.plot(tt,tr); plt.title("Thermal deviation response"); plt.xlabel("Time (hours)"); plt.ylabel("Temperature deviation"); plt.legend(["Zone A","Zone B","Zone C"]); plt.grid(True); plt.tight_layout(); plt.show()

eigenvalues [-0.5359+0.j -0.3573+0.j -0.1768+0.j]

	Zone A	Zone B	Zone C
0.0	2.0000	0.2000	1.0000
4.0	0.4321	0.3274	0.4242
8.0	0.1406	0.1839	0.2078
16.0	0.0275	0.0454	0.0518

結果の読み取り

固有値の実部は負で、偏差は減衰します。「安定」だけでなく許容帯への復帰時間を保全計画へ落とし、入力・飽和・遅れも別途検証します。

No.060：行列の安定性

実務での意味

漸近安定性、短期最大増幅、安定余裕を統合し、アラートと再同定の基準を決めます。

分析・モデル化の考え方

離散時間系は $\rho(A)<1$ なら漸近安定です。$\|A^k\|_2$ は最悪方向の過渡増幅を表します。境界1ではなく推定誤差込みの社内閾値を置きます。

Pythonで確認する

hh=np.arange(31); gain=np.array([np.linalg.norm(np.linalg.matrix_power(A,int(k)),2) for k in hh]); pi=int(np.argmax(gain))
display(pd.DataFrame({"metric":["spectral radius","margin","peak gain","peak shift"],"value":[rho,1-rho,gain[pi],pi]}).round(4))
plt.plot(hh,gain,marker="o",ms=3); plt.axhline(1,color="red",ls="--"); plt.title("Worst-case propagation gain"); plt.xlabel("Shifts ahead"); plt.ylabel(r"$\|A^k\|_2$"); plt.grid(True); plt.tight_layout(); plt.show()

	metric	value
0	spectral radius	0.7541
1	margin	0.2459
2	peak gain	1.0000
3	peak shift	0.0000

結果の読み取り

モデルは漸近安定で顕著な過渡増幅もありません。ただし係数の信頼区間でストレステストし、たとえば $ho(A)\ge0.90$ を再同定の社内トリガーにします。

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

今回のモデルは安定ですが、「1未満」だけが結論ではありません。半正定値損失には見えない共通モードがあり、支配モードでは熱処理の寄与が大きく、条件数は原因寄与の逆算が不安定になり得ることを示しました。簡易境界、支配値更新、大規模部分計算、全体監査を用途別に使い分けます。

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

1. 時刻同期、欠測処理、標準化、品種・設備状態を統一する
2. 学習期間外で予測誤差と残差を検証する
3. 条件数、固有対残差、許容誤差、ライブラリ版を記録する
4. 安定余裕、過渡ピーク、復帰時間を品質・保全基準へ対応づける
5. ドリフト検知、再学習条件、承認者、ロールバックを決める

数学的安定性は安全性や品質保証そのものではありません。FMEA、管理図、設備知識、実験計画と組み合わせます。

まとめ

No.051〜060では、行列構造、固有値境界、支配モード、数値感度、計算法、時間応答を一貫した安定性監査へつなげました。重要なのは、どの仮定で何を保証し、どの条件で再評価するかを明文化することです。

法人向けのご相談

数理工房では、設備・品質データ整理から工程伝播モデル、異常兆候監視、シミュレーション、現場説明、運用設計まで支援します。

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