放電加工機の異常検知
１．はじめに-------------------------------------------------日常の経験で、モノが壊れそうかなという勘がはたらくことがあると思います。
「洗濯機の回転音が変だ」 とか 「冷蔵庫がブーンと鳴り出す」とかこのような変化を定量的に（定性でもいいですが）データ化して、可視化できれば、故障の事前対応も可能かと思いました。
たしかに、極めるには、「線形代数学」とか、「統計学」への造詣が必要で、私は不足しているのですが、とりあえず、データを採取して検知してみて、使い物になるか否かを検討してみようと思います。
故に、2 項目は飛ばしていただいて、3 項目を見ていただく方が、話しが早いかと思います。
また、お願いですが、他の装置機器で、物理量（真空度、温度、電圧等）のデータを採取できないか、探しています。思い当たるものあれば、是非、一報ください。
----------------------------------------------------技術部の仕事として、装置の保守管理が挙げられ、装置が必要な時に、必要なだけ利用できる状態にするのは、技術部の使命ともいえます。
また、装置には、故障という使用不能になる状態が発生し、これを無くすることは、不可能に近いです。
部品の寿命、利用者の誤使用、外乱等 回避できないな事象があります。
そこで、少なくとも、故障を予期して、早めの部品交換、装置の状態確認ができれば、修理というプロセスは避けられないかもしれませんが、非稼働を最小限にすることは可能と考えます。
すなわち、機械学習を使って、装置の出力を、客観的な基準で変化の兆候をとらえて（＝異常検知）、対応することを、考えてみました。
また、対象する装置は、担当している放電加工機で行いました。
参考にした本は、「入門 機械学習による異常検知―R による実践ガイドー」 井出剛 著 コロナ社です。 R のプログラム構文は この本を、参照してます。

1２．時系列データの異常検知
異常検知の方法として、基礎的なものとして、・近傍法による異常部位検出 （検証データと訓練データでその差異を評価する）
・特異スペクトル変異法 （時系列データの変化を検知し、変化の程度を評価する）の 二つで試みてみました。
2.1 考え方2.1.1 近傍法による異常部位検出
観測値として、長さ T の時系列があります。ξ()は観測値で、T は 20 個とします。
D={ξ(1), ξ(2), ξ(3), ξ(4), …………….., ξ(17), ξ(18), ξ(19), ξ(20)} とします。
隣接した観測値データをｗ個として、ベクトルの集まりとします。w=3 とします。
ξ(1)

ξ(2)

ξ(17)

ξ(18)

x(1)= ξ(2)x(2)= ξ(3)…………………………x(17)= ξ(18)x(18)= ξ(19){ξ(3)}{ξ(4)}{ξ(19)}{ξ(20)}

18 個の 3 次元ベクトルに変換します。w 個の集まりを、スライド窓と言います。

(3)
ξ(1) (2)(4) (5)(11)(9)(6) (7) (8)(10)(14) (15)(12) (13)(17) (19)(16) (18)
(20)

x(1)
18 個のスライド窓

x(18)時系列データを
D={x(1),x(2),x(3),-----------------------------,x(17),x(18)} として、扱います。
次に、データを、訓練データと検証データを用意します。訓練データとは、異常ないと考えるデータです。
検証データとは、異常が無いかを確認するデータです。
各々の検証データについて、訓練データからどれくらいのズレがあるか距離で算出します。
距離の算出は、R が持つパッケージ FNN の knnx.dist 関数 で近傍距離を算出します。

2

近傍距離の概略を説明します。
3 次元ベクトルで進めてましたが、図にすると分かりにくいので、２次元で進めます。
黒丸が訓練データのグループです。
青丸の検証データの最も接近した黒丸との距離を近傍距離とします。

近傍距離を異常度として、時系列グラフを作成します。
閾値を超えた異常度が発生した時点を、検証することになります。

2.1.2 特異スペクトル変異法
これは、変化の度合いで、異常を判定する方法です。
観測値として、長さ T の時系列があります。ξ()は観測値で、T は 20 個とします。
D={ξ(-10), ξ(-9), ξ(-8), ξ(-7), …………….., ξ(6), ξ(7), ξ(8), ξ(9)} とします。
隣接した観測値データをｗ個として、ベクトルの集まりとします。w=3 とします。
ξ(0)を現在の値として、ξ(-1)からξ(-10)を過去の履歴行列 X とします。
ξ(0)からξ(9)をテスト行列 Z とします。
ξ(-2)とξ(-1)が両行列で重なりますが、相互位置を定める整数で L(ラグ)とします。
ξ(-10)
(-9)

(-8)

(-6)
(-7)

(0)(-2)(-5)(-1)(-4) (-3)(3) (4)(1) (2)

(6)
(5)

(8)
(7)

(9)

3

現在時点での X と Z の行列は以下のようになります。
𝜉(−10) 𝜉(−9) 𝜉(−8) 𝜉(−7) 𝜉(−6) 𝜉(−5) 𝜉(−4) 𝜉(−3)X(0)=( 𝜉(−9) 𝜉(−8) 𝜉(−7) 𝜉(−6) 𝜉(−5) 𝜉(−4) 𝜉(−3) 𝜉(−2))
𝜉(−8) 𝜉(−7) 𝜉(−6) 𝜉(−5) 𝜉(−4) 𝜉(−3) 𝜉(−3) 𝜉(−1)
𝜉(−2) 𝜉(−1) 𝜉(0) 𝜉(1) 𝜉(2) 𝜉(3) 𝜉(4) 𝜉(5) 𝜉(6) 𝜉(7)Z(0)=(𝜉(−1) 𝜉(0) 𝜉(1) 𝜉(2) 𝜉(3) 𝜉(4) 𝜉(5) 𝜉(6) 𝜉(7) 𝜉(8))
𝜉(0)𝜉(1) 𝜉(2) 𝜉(3) 𝜉(4) 𝜉(5) 𝜉(6) 𝜉(7) 𝜉(8) 𝜉(9)Z の X に対する変化度合いを評価するのですが、各行列のパータンを 特異値分解 という手法でもとめます。
特異値分解の説明です。

上記の X(0)は 3 列 8 行の行列ですが、X(0) = UΣV* に分解します。 U は 3 列 3 行の行列、V*は 8 列 8 行の行列です。
Σは対角行列で、対角成分のσ1,・・・・・σq と並び、部分行列は、０行列です。
また、σ1 ≧…≧ σq ≧0 の関係にします。σ1 が最大値で、行列 X(0)に最も大きく影響を与えていることになります。σの値を大きい方から、m 個採取して、Um(0)の行列を求めます。
同様に、Z(0) においても Qm(0)の行列 を求めて、変化度を算出します。
（文字がダブルので Q としました。 Z 側の U 行列とみなしてください）
次に、変化度を求めます。
「行列２ノルム」という量で求めます。 U と Q との空間距離を、||Um(0)*Qm(0)||2^2 で算出します。（詳細は、別途に譲ります）
変化の大きさ = 1- ||Um(0)*Qm(0)||2^2 = 1-( Um(0)*Qm(0)の最大特異点)^2とします。
特異値分解、変化度は、R の svd 関数で算出します。

4

3. 放電加工機での取り組み
3.1 メッセージデータの内容①加工機は、自己診断機能を有しており、3 段階のレベルで、それに応じたメッセージを表示します。
・エラーメッセージ --- 続行不可能な障害が発生した時に表示され、動作は中断します。 170 種類・ハルトメッセージ --- 再開可能な障害が発生した時に表示され、動作は一時停止します。78 種類
・コメントメッセージ --- 続行可能な障害が発生した時に、注意として表示されます。 58 種類（注意というより、プロセスを表示する事が多いです）②発生時間[2022-06-29 20:17:43] のような 形式で表示されます。
③ その他
XYZ 座標値、プログラム名（プログラム稼働時）等 表示されます。
実際の 1 回分の実データは以下です。エクセルに変換できます。

53.2 前処理されたデータ（＝元データ）の内容
今回は 4 月 1 日から 6 月 20 日の分について、データを作成しました。6317 個あります。
実データから下記の元データを作成するのに、手動では面倒ですので、マクロを活用してます。
Point とは メッセージに対応する点数で、障害の程度で付与してます。
エラーメッセージ 3～４ ハルトメッセージ ２～2.5 コメントメッセージ 0.5～1.5 としました。
同じメッセージ内でも、数値を加減してます。
以下の元データを作成します。（エクセルから CSV(コンマ区切り)に変換します）

whenpoint

2022/04/01 08:28:01

0.5

2022/06/20 08:35:57

12022/06/20 08:36:10

1

元データから、異常度を評価します。
計算は R 言語 で行いました。 R 言語とは、「オープンソース・フリーソフトウエアの統計解析向けのプログラミング言語」 とのことです。
Python の機械学習もありますが、参考書が R による実践ガイド だったこともあり、R で検証してます。
また、Python も試行しましたが、R の方が異常度が顕著でした。

元データを Rstudio（R の開発環境）に取込んで、・近傍法による異常部位検出 と ・特異スペクトル変異法 との 2 種について、計算して、プロットしました。

6

3.3 結果のプロット
3.3.1 近傍法による異常部位検出（スライド窓 w=10 で 計算しました。）
元データは黒のグラフです。6317 個のデータに対して、0.5 から 4.0 のポイントが付与されてます。
異常度のグラフは赤のグラフです。1000 から 2000 の間で０になってるのは、この間は安定していたので訓練データとして採用しました。訓練データ同士を比較したので、異常度は０です。

２

１354

異常度 5 以上に注目して、メッセージを確認しました。
1⃣、2⃣、4⃣ は [H00207: () ﾃﾝｼｮﾝ制御装置が異常です。ENT ｽｲｯﾁを押して下さい。]

3⃣ は [E00245: () 機械座標が確定していないため実行できません。]

5⃣ は [E00350: (V) ﾌｨｰﾄﾞﾊﾞｯｸ速度異常です。]

が 多発していました。

7

3.3.2 特異スペクトル変異法（w=50, svd 次元数 m=1, svd 算出用部分時系列 k=w/2, ラグ L=k/2 で計算しました、）
黒の元データは 3.3.1 と同じです。異常度のグラフは赤のグラフです１

異常度 0.10 以上に注目して、メッセージを確認しました。
1⃣ は [H00207: () ﾃﾝｼｮﾝ制御装置が異常です。ENT ｽｲｯﾁを押して下さい。]

が多発してました。

8

3.3.3 プロットへの評価
元データの時点で、異常度のランク付けをしてるので、当然の結果かもしれませんが、異常度の高い時点で、不具合が発生しそうな事象が起きています。
近傍法のグラフで1⃣、2⃣、4⃣ で テンション装置が異常とあるので、予防修理が必要です。
3⃣ の 機械座標の件は、加工機が使用不可になりました。
（コールセンターからの対応で復帰しましたが）
5⃣ のフィードバック異常は、今までにないメッセージです。
元データを眺めてるだけでは、ここまでは、分からないです。
異常度グラフの方が判別しやすいのは、明白です。

3.3.4 今後の取組継続して、データ採取して、どれくらい役に立つのか検証していきます。
時系列データの異常検知は、他の手法のあるよいですので、検討してみます。

9