JP2019003389A - 異常診断装置、異常診断方法及び異常診断プログラム - Google Patents

Abstract

【課題】異常が発生する確率を高い精度で予想する。【解決手段】本発明の異常診断装置１は、診断対象機器から過去における特徴量を時系列で取得し、取得した特徴量のばらつきを示す異常度を時系列で算出し、診断対象機器に異常が実際に発生した確率である異常発生確率と、異常度との相互関係を決定し、決定した相互関係に基づいて、異常発生確率を時系列で算出する事前情報処理部２１と、診断対象機器から診断対象期間における特徴量を時系列で取得し、取得した診断対象期間における特徴量のばらつきを示す異常度を時系列で算出し、決定した相互関係に基づいて、診断対象期間において診断対象機器に異常が発生する確率である異常発生予想確率を時系列で算出する運用情報処理部２２と、を備えることを特徴とする。【選択図】図７

Description

本発明は、異常診断装置、異常診断方法及び異常診断プログラムに関する。

稼働中の機器が故障すると、当該機器が組み込まれている生産設備等は機能を停止する。予めその機器の故障が予想されている場合、機器の保守・更新等を計画的に行うことができる。しかしながら、機器の故障が予想されていない場合、故障が発生して初めて復旧作業が行われることになる。特に生産設備等の規模が大きいと、復旧に要する時間及びコストは、膨大になる。したがって、機器の異常をいち早く予想することが重要である。近時、機器から取得したセンサ値に現れる異常の予兆を検出する技術が普及している。

特許文献１の減速機異常診断装置は、モータのトルクの伝達効率を診断する。当該減速機異常診断装置は、予め正常時の伝達効率を記憶しておき、診断時の伝達効率と、記憶している伝達効率とのマハラノビス距離を算出する。そして、当該減速機異常診断装置は、マハラノビス距離が所定の閾値を超えた場合、減速機に異常が発生していると判断する。

特許文献２のプラント監視装置は、プラントの運転データを診断する。当該プラント監視装置は、予め正常時の運転データを記憶しておき、診断時の運転データと、記憶している運転データとのマハラノビス距離を算出する。運転データは、複数の項目を含み、当該プラント監視装置は、マハラノビス距離の低下に最も寄与する項目（引き下げ因子）を削除したうえで、診断時の運転データと、記憶している運転データとのマハラノビス距離（修正マハラノビス距離）を算出する。そして、当該プラント監視装置は、修正マハラノビス距離が所定の閾値を超えた場合、プラントに異常が発生していると判断する。

特許文献３の管理システムは、半導体製造装置を診断する。当該管理システムは、正常な動作状態においてサンプリングされたデータ群に基づきマハラノビス空間を作成する。そして、当該管理システムは、動作状態で得られるデータのマハラノビス距離を算出し、マハラノビス距離が所定の閾値を超えた場合、半導体製造装置に異常が発生していると判断する。

特開２０１２−１９４０３５号公報 特許第５６１０６９５号公報 特許第３３４９４５５号公報

前記した３つの先行技術に共通して以下のことがあてはまる。すなわち、サンプルとする正常状態のデータに異常値が含まれてしまうことがある。これを回避しようとして閾値を大きくすると、本来検出するべき異常が検出されなくなってしまう（図２として詳細後記）。さらに、これらを回避しようとして、多くのデータの分布を代表値（平均値、分散、最大値、最小値等）に換算し、代表値同士のマハラノビス距離を算出すると、有意な細部情報が消失してしまう。また、サンプルとなるデータ母数に比して異常値の数が極端に少ない場合、又は、正常期間と異常期間との間で、データのばらつきに有意な差異がない場合、ばらつきを判断基準とする判定は困難になる。
そこで、本発明は、異常が発生する確率を高い精度で予想することを目的とする。

本発明の異常診断装置は、診断対象機器から過去における特徴量を時系列で取得し、取得した特徴量のばらつきを示す異常度を時系列で算出し、診断対象機器に異常が実際に発生した確率である異常発生確率と、異常度との相互関係を決定し、決定した相互関係に基づいて、異常発生確率を時系列で算出する事前情報処理部と、診断対象機器から診断対象期間における特徴量を時系列で取得し、取得した診断対象期間における特徴量のばらつきを示す異常度を時系列で算出し、決定した相互関係に基づいて、診断対象期間において診断対象機器に異常が発生する確率である異常発生予想確率を時系列で算出する運用情報処理部と、を備えることを特徴とする。
その他の手段については、発明を実施するための形態のなかで説明する。

本発明によれば、異常が発生する確率を高い精度で予想することができる。

マハラノビス距離と母集団の標準偏差との関係を説明する図である。 既存技術におけるマハラノビス距離による異常判定の例を説明する図である。 異常度を説明する図である。 （ａ）、（ｂ）及び（ｃ）は、時間幅と診断周期との関係を説明するである。 （ａ）、（ｂ）及び（ｃ）は、時間幅に応じた異常度の現れ方を説明する図である。 （ａ）、（ｂ）及び（ｃ）は、異常度及び異常発生確率の相互関係を説明する図である。 異常診断装置の構成を説明する図である。 計測値情報の一例を示す図である。 時間幅情報の一例を示す図である。 重み情報の一例を示す図である。 異常原因情報の一例を示す図である。 事前情報処理手順のフローチャートである。 運用情報処理手順のフローチャートである。 事前情報入力画面の一例である。 （ａ）、（ｂ）、（ｃ）及び（ｄ）は、事前情報処理手順において画面表示されるデータ例である。 事前情報表示画面の一例である。 （ａ）、（ｂ）及び（ｃ）は、運用情報処理手順において画面表示されるデータ例である。 運用情報表示画面の一例である。

以降、本発明を実施するための形態（“本実施形態”と呼ぶ）を、図等を参照しながら詳細に説明する。本実施形態の異常診断装置が診断対象とする機器（“診断対象機器”と呼ぶ）は、その機器から制御データ又は計測データを取得することができるような任意の機器である。そして、これらの制御データ又は計測データを診断対象機器の“特徴量”と呼ぶ。診断対象機器の種類に応じて、様々な特徴量が存在する。例えば、診断対象機器が冷凍機である場合、その特徴量は、圧縮機の電圧及び電流、特定の箇所における冷媒の温度及び圧力、圧縮機の回転速度及び振動、冷凍機から所定の距離を隔てた位置における騒音等である。そして特徴量は、これらの制御データ又は観測データを加工して得られるデータ（例えば、騒音の時系列波形をフーリエ変換して得られる、特定周波数のスペクトル強度）も含む。

（マハラノビス距離）
図１に沿って、マハラノビス距離と母集団の標準偏差との関係を説明する。周知のように、複数のデータから構成される母集団を代表する基準値（平均値、重心等）と個々のデータとの間のマハラノビス距離を算出することが一般的に行われる。ユークリッド距離とは異なり、マハラノビス距離は、母集団の標準偏差（分散の平方根）を反映している。

いま、異常診断装置が診断対象機器の特徴量（例えば冷媒の温度）を、１分の計測周期で１０回計測する例を想定する。このような計測を行った結果が、図１（ａ）の１０個の点に示されている。１０回の計測時点のうち、ｉ回目に計測された特徴量をｘ_ｉとする。異常診断装置は、ｘ_ｉのマハラノビス距離Ｄｍ（ｘ_ｉ）を、次の式１を使用して算出する。ここで、μは、１０個の特徴量の平均値であり、σは、１０個の特徴量の標準偏差である。
Ｄｍ（ｘ_ｉ）＝｜ｘ_ｉ−μ｜／σ （式１）

このような計測を行った他の例が、図１（ｂ）に示されている。ここでも、異常診断装置は、診断対象機器の特徴量を、１分の計測周期で１０回計測している。図１（ａ）及び図１（ｂ）を比較すると、特徴量の平均値μ及びｉ回目に計測された特徴量ｘ_ｉの値は同じである。しかしながら、図１（ａ）よりも図１（ｂ）の方が、標準偏差σは大きい。すると、式１から明らかなように、図１（ｂ）の方がσが大きい分、図１（ａ）よりも図１（ｂ）の方が、式１で算出されるＤｍ（ｘ_ｉ）は小さくなる。このことは、あるｘ_ｉの値が平均値から相当程度乖離していても、母集団のσが充分に大きい場合は、当該ｘ_ｉは目立たなくなることを意味している。

（特徴量の次元）
前記した特徴量は、１次元である。今度は、異常診断装置が診断対象機器のｎ（ｎ＝１、２、３、・・・）種類の特徴量（例えば、冷媒の温度、冷媒の圧力、圧縮機の回転速度、・・・）を、１分の計測周期で１０回計測する例を想定する。この例では、ｉ回目に計測された特徴量ｘ_ｉは、ｎ次元のベクトルとなる。すると、異常診断装置は、ｎ次元空間内に１０個の特徴量をドットした場合のそれらの重心と、ｘ_ｉとの間のマハラノビス距離Ｄｍ（ｘ_ｉ）を、次の式２を使用して算出することができる。ここで、μは、重心を示すｎ次元のベクトルであり、Ｒは、分散共分散行列である。“√”は平方根を示し、添え字“Ｔ”は転置を示し、添え字“−１”は逆行列を示す。
Ｄｍ（ｘ_ｉ）＝√（（ｘ_ｉ−μ）^Ｔ×Ｒ^−１×（ｘ_ｉ−μ）） （式２）

なお、Ｒは、ｎ次元空間内における母集団のいわば“形状”を示しており、２つの特徴量同士の共分散が大きい方向のＤｍ（ｘ_ｉ）は小さくなる。いま、説明を単純化するために、ｎ＝２が成立し、２種類の特徴量は、冷媒の圧力及び圧縮機の回転速度であるとする。これら２種類の特徴量は、強い正の相関を有する。両者の関係を平面にドットすると、その外縁は楕円となり、その楕円の長軸は、右上りの直線となる。このとき、長軸上のマハラノビス距離のスケール（１目盛の刻み幅）は、短軸上のマハラノビス距離のスケールよりも大きくなる。つまり、楕円と長軸との交点の（重心からの）マハラノビス距離と、楕円と短軸との交点の（重心からの）マハラノビス距離と比較する場合、見かけの長さの相違にかかわらず、両者は、ほぼ一致する。

以降では単純化のため、１次元の特徴量のマハラノビス距離を例として説明する。しかしながら、本実施形態を、ｎ次元の特徴量のマハラノビス距離に一般化することも可能である。

（マハラノビス距離による異常判定）
図２に沿って、既存技術におけるマハラノビス距離による異常判定の例を説明する。図２の座標平面の横軸は時間であり、縦軸はある特徴量（例えば温度）である。いま、時点ｔ_ａ〜時点ｔ_ｂの期間においては、診断対象機器が正常であることが既知であり、時点ｔ_ｂ〜時点ｔ_ｃの期間４１においては、診断対象機器が異常であることが既知であるとする。そして、特徴量の平均値μを基準にして、マハラノビス距離の上限の閾値がＴｈ_Ｕであり、下限の閾値がＴｈ_Ｌであるとする。

すると、正常である期間において、特徴量４２a及び４２ｂが“異常”と判定されてしまう。そこで、上限Ｔｈ_Ｕをより上方に移動し、下限Ｔｈ_Ｌをより下方に移動する。すると、異常である期間４１において、特徴量４３a、４３ｂ、４３ｃ及び４３ｄが“異常”とは判定されなくなってしまう。このように、特徴量のマハラノビス距離に対し単純に閾値を適用するだけでは、異常診断の精度は向上しない。

（異常度）
図３に沿って、異常度を説明する。異常診断装置は、例えば１０分の時間幅を有する時間帯４４において、１分の計測周期で特徴量を計測している。異常診断装置は、１０個の特徴量のそれぞれについて、マハラノビス距離を算出する。異常診断装置は、算出した個々のマハラノビス距離の大きさに応じた重みを、そのマハラノビス距離に乗算する。図３において、マハラノビス距離が大きいほど、換言すれば特徴量が平均値μから乖離するほど、乗算される重みは大きい。

図３では、紙面の制約からすべての特徴量に重みが付されている訳ではないが、平均値μから大きく乖離した特徴点４５ａ及び４５ｂには、重み“１００”が付されている。平均値μに近い特徴点４７には、重み“１”が付されている。その中間の特徴点４６ａ及び４６ｂには、重み“１０”が付されている。

異常診断装置は、重みを乗算されたマハラノビス距離の平均値を、その時間帯についての“異常度”とする。つまり、図３において、式１が成立し、次の式３も成立する。
異常度＝Σ（ｗ_ｉ×Ｄｍ（ｘ_ｉ））／ｎ （式３）
ここで、ｗ_ｉは、Ｄｍ（ｘ_ｉ）の大きさに応じた重みであり、Ｄｍ（ｘ_ｉ）が大きいほど、ｗ_ｉは大きい。ｎは、時間帯に含まれるマハラノビス距離の数（つまり特徴量の数）であり、前記の例では、ｎ＝１０である。“Σ”は、和を示す。
このような処理を例えば１日（１４４０分）繰り返すと、異常診断装置は、連続する１４４個の時間帯のそれぞれに対応する１４４個の異常度を取得することになる。

異常度は、式３が定義するものに限定されない。異常度は、特徴量のばらつきを示す数値であればどのようなものでもよい。標準偏差、分散、最大値と最小値との差異、特徴量が正規分布に従うと仮定した場合の所定の基準を満たす（上位５％又は下位５％に入る）はずれ値の個数、そのようなはずれ値の比率等、様々なものが異常度となり得る。

（用語の定義）
計測周期とは、異常診断装置が特徴量を計測する周期であり、例えば１分である。計測時点とは、所定の計測周期で実際に特徴量が計測される時点である。診断周期とは、異常診断装置が診断対象機器を診断する周期であり、例えば１０分である。診断時点とは、所定の診断周期で実際に診断対象装置が診断される時点である。１つの診断対象時点に複数（例えば１０個）の計測時点が対応する。時間帯とは、複数の特徴量が含まれる連続した時間区分であり、異常診断装置は、１つの時間帯ごとに１つの異常度を算出する。１つの診断時点に１つの時間帯が対応しているので、異常診断装置は、１つの診断時点について１つの異常度を算出することになる。時間幅とは、時間帯の長さであり、例えば１０分である。

（時間幅と診断周期）
図４に沿って、時間幅と診断周期との関係を説明する。図４（ａ）は、時間幅が診断周期に一致する例である。この場合、複数の時間帯（破線の長方形）が時系列で隙間がなく、かつ、重複もない状態で並ぶ。図４（ｂ）は、時間幅が診断周期よりも長い例である。この場合、複数の時間帯が、その一部同士を重複させながら時系列で並ぶ。図４（ａ）及び図４（ｂ）においては、時間幅は固定されている。図４（ｃ）においては、時間幅は可変長であり、時間的に後の時間幅ほど長い。そして、各時間幅の始点がある時点（ｔ_０）に固定されている。異常診断装置は、診断対象機器の特性、特徴量の数等に応じて、任意に時間幅及び診断周期を設定することができる。以降では、図４（ａ）を想定して説明を続ける。

図５に沿って、時間幅に応じた異常度の現れ方を説明する。前記したように、異常診断装置は、１つの診断時点について１つの異常度を算出する。すると、異常診断装置は、診断周期を充分小さくすることによって、図５（ａ）〜（ｃ）のような、時系列の異常度の波形を描画することができる。図５（ａ）は、時間幅が１時間である場合の時系列の異常度である。図５（ｂ）は、時間幅が２４時間である場合の時系列の異常度である。図５（ｃ）は、時間幅が４８時間である場合の時系列の異常度である。一般に、時間幅が短いほど、異常度は高周波的に現れる。

（異常発生確率）
図６に沿って、異常発生確率を説明する。図６（ａ）は、時系列の異常度の波形の一例である。いま、過去の期間４８において、診断対象機器が異常であったことが既知であるとする。異常診断装置は、異常度のある水準（例えば“１０”）において、横軸に水平な直線を引く。すると、直線は、１０箇所（点１〜点１０）で異常度の波形と交差又は接する。このうち、期間４８に含まれる点は２つ存在し（点７及び点８）、期間４８に含まれない点は８つ存在する（点１〜点６、点９及び点１０）。このとき、異常診断装置は、２／（２＋８）×１００＝２０％を“異常発生確率”とする。

異常診断装置は、異常度の他のある水準（例えば“４０”）において、横軸に水平な直線を引く。すると、直線は、４箇所（点１１〜点１４）で異常度の波形と交差又は接する。このうち、期間４８に含まれる点は２つ存在し（点１３及び点１４）、期間４８に含まれない点は２つ存在する（点１１及び点１２）。このとき、異常診断装置は、２／（２＋２））×１００＝５０％を“異常発生確率”とする。

図６（ｂ）の波形は、図６（ａ）と同じ波形である。また、図６（ａ）と同様に、過去の期間４８において、診断対象機器が異常であったことが既知であるとする。異常診断装置は、異常度のある水準（例えば“４０”）において、横軸に水平な直線を引く。すると、直線は、４箇所（点１１〜点１４）で異常度の波形と交差又は接する。この直線のうち、波形が直線を上側に超える又は直線に接する部分の長さは、ｄ_１＋ｄ_２である。さらにそのうち、期間４８に含まれる部分の長さは、ｄ_２であり、期間４８に含まれない部分の長さは、ｄ_１である。このとき、異常診断装置は、ｄ_２／（ｄ_１＋ｄ_２）×１００の値（％）を“異常発生確率”としてもよい。請求項における“頻度”には、ここでの点１１〜点１４の数、又は、部分の長さｄ_１、ｄ_２が相当する。

異常診断装置は、異常度の水準を少しずつ変化させたうえでこのような処理を繰り返す。すると、異常診断装置は、“（異常度，異常発生確率）＝（１０，２０％）”のような異常度及び異常発生確率の組合せを多く取得することになる。図６（ｃ）は、このような異常度及び異常発生確率の相互関係を、横軸を異常度とし縦軸を異常発生確率とする座標平面に表したものである。異常診断装置は、公知の方法で、当該相互関係を滑らかな曲線に近似してもよい。近似した相互関係の曲線は、例えば、上下（１００％及び０％）に漸近線を有する右上りのシグモイド曲線のような曲線になる。

（異常診断装置の構成）
図７に沿って、異常診断装置１の構成を説明する。異常診断装置１は、一般的なコンピュータである。異常診断装置１は、中央制御装置１１、入力装置１２、出力装置１３、主記憶装置１４及び補助記憶装置１５を備える。これらはバスで接続されている。補助記憶装置１５は、計測値情報３１、時間幅情報３２、重み情報３３及び異常原因情報３４（詳細後記）を格納している。主記憶装置１４における事前情報処理部２１及び運用情報処理部２２は、プログラムである。以降の説明において、“○○部は”と動作主体を記した場合、それは、中央制御装置１１が補助記憶装置１５から○○部を読み出し、主記憶装置１４にロードしたうえで○○部の機能（詳細後記）を実現することを意味する。

（計測値情報）
図８に沿って、計測値情報３１を説明する。計測値情報３１は、計測時点欄１０１に記憶された計測時点に関連付けて、特徴量１欄１０２には特徴量１（ここでは電圧）、特徴量２欄１０３には特徴量２（ここでは電流）、・・・、特徴量ｎ欄１０９には特徴量ｎが記憶されている。図８において“＃”は、異なる複数の測定値（又は制御値）を省略的に示している。なお、請求項における“項目”には、ここでの電圧、電流、温度、圧力、・・・が相当する。異常診断装置１は、図示しない診断対象機器と図示しないネットワークを介して接続されている。そして異常診断装置１は、過去及び現在における特徴量を診断対象機器から取得し、計測値情報３１として記憶している。

（時間幅情報）
図９に沿って、時間幅情報３２を説明する。時間幅情報３２においては、パタンＩＤ欄１１１に記憶されたパタンＩＤに関連付けて、始点欄１１２には時間幅の始点（時間軸上の１点）が、終点欄１１３には時間幅の終点（時間軸上の１点）が、診断周期欄１１４には診断周期が記憶されている。
パタンＩＤ欄１１１のパタンＩＤは、時間幅の始点、時間幅の終点及び診断周期の組合せ（パタン）を一意に特定する識別子である。

始点欄１１２の時間幅の始点は、時間幅の始点を定義する情報であり、ここでは、時間幅の始点が診断時点に対してどの程度時間的に遡るかを示す情報である。なお、時間幅の始点は、特定の“固定時点”であってもよい。
終点欄１１３の時間幅の終点は、時間幅の終点を定義する情報であり、ここでは一律に前記した診断時点である。診断時点は、例えば図４（ａ）のｔ_１、ｔ_２、ｔ_３、ｔ_４、・・・である。
診断周期欄１１４の診断周期は、前記した診断周期である。
図９のパタン“Ｐ０１”は、図４（ａ）に対応する。パタン“Ｐ０４”は、図４（ｂ）に対応する。パタン“Ｐ０５”は、図４（ｃ）に対応する。

（重み情報）
図１０に沿って、重み情報３３を説明する。重み情報３３においては、マハラノビス距離欄１２１に記憶されたマハラノビス距離に関連付けて、重み欄１２２には重みが記憶されている。
マハラノビス距離欄１２１のマハラノビス距離は、特徴量のマハラノビス距離の区分である。ここで区分の境目となる２、３、４、・・・は、あくまでも一例である。

重み欄１２２の重みは、特徴量のマハラノビス距離に対して乗算される重みである。ここでの重みは、１０のべき乗であり、マハラノビス距離が大きくなるほど、重みの値も大きくなる。ここでの重みもまた、あくまでも一例に過ぎない。他の例として、重みは、マハラノビス距離が正規分布に従って発生すると仮定した場合の発生確率の逆数であってもよい。
図１０の例では、マハラノビス距離がある整数値に達すると重みが階段的に増加し、その後次の整数値に達するまで、重みは同じ水準を維持する。しかしながら、重みは、マハラノビス距離の増加に応じて、連続的に滑らかに増加してもよい。

（異常原因情報）
図１１に沿って、異常原因情報３４を説明する。異常原因情報３４においては、特徴量欄１３１に記憶された特徴量に関連付けて、異常名称欄１３２には異常名称が、異常原因欄１３３には異常原因が記憶されている。
特徴量欄１３１の特徴量は、前記した特徴量（電圧、電流、温度、・・・）である。
異常名称欄１３２の異常名称は、診断対象機器の異常の態様である。
異常原因欄１３３の異常原因は、診断対象機器の異常の原因である。異常原因は、診断対象機器を保守員が回復させるために手掛かりとなる情報であってもよい。

以降、本実施形態の処理手順を説明する。処理手順は２つ存在し、それらは、事前情報処理手順（図１２）及び運用情報処理手順（図１３）である。事前情報処理手順が実行される前提として、計測値情報３１（図８）、時間幅情報３２（図９）及び重み情報３３（図１０）が完成状態で補助記憶装置１５に格納されているものとする。そして、運用情報処理手順が実行される前提として、事前情報処理手順が既に実行されているものとする。

（事前情報処理手順）
図１２に沿って、事前情報処理手順を説明する。
ステップＳ２０１において、異常診断装置１の事前情報処理部２１は、データ取込期間及びパタンを受け付ける。具体的には、第１に、事前情報処理部２１は、出力装置１３に事前情報入力画面５１（図１４）を表示する。

第２に、事前情報処理部２１は、ユーザが事前情報入力画面５１のデータ取込期間欄５１ａに対して開始日時及び終了日時を入力するのを受け付ける。例えば、現時点が２０１７年５月３０日であり、ユーザは、２０１７年１月１日から２０１７年３月３１日までの特徴量を、過去の時点におけるサンプルデータとして使用したいとする。この場合、ユーザは、“開始日時”として“２０１７年１月１日００時００分００秒”を入力し、“終了日時”として“２０１７年４月１日００時００分００秒”を入力する。
第３に、事前情報処理部２１は、時間幅情報３２（図９）を出力装置１３に表示し、ユーザが任意のパタンＩＤを選択するのを受け付ける。

ステップＳ２０２において、事前情報処理部２１は、特徴量、異常発生期間、異常名称及び異常原因を受け付ける。具体的には、第１に、事前情報処理部２１は、ユーザが事前情報入力画面５１の特徴量選択欄５１ｂの候補のうちから１つの特徴量を選択するのを受け付ける。
第２に、事前情報処理部２１は、ユーザが特徴量選択欄５１ｂの特徴量を選択した状態で、異常発生期間欄５１ｃに対して、異常名称、異常原因、開始日時及び終了日時を入力するのを受け付ける。

ユーザは、例えば“２０１７年３月１０日１０時００分００秒から２０１７年３月１５日１２時００分００秒までの期間、圧縮機の電圧が過大となった。その原因は、工程ＡにおけるＢ処理の不良であった”ということを知っているとする。この場合、ユーザは、“特徴量（電圧）”選択した状態で、“開始日時”として“２０１７年３月１０日１０時００分００秒”を入力し、“終了日時”として“２０１７年３月１５日１２時００分００秒”を入力する。そして、ユーザは、“異常名称”として“圧縮機過電圧”を入力し、“異常原因”として、“工程ＡのＢ処理”を入力する。

事前情報処理部２１は、ステップＳ２０２の“第１”及び“第２”の処理を任意の回数繰り返す。すると、事前情報処理部２１は、特徴量ごとに、異常名称、異常原因及び異常発生期間（開始日時及び終了日時）を受け付けることになる。
第３に、事前情報処理部２１は、ステップＳ２０２の“第１”及び“第２”において受け付けた情報に基づいて、異常原因情報３４（図１１）を作成し、補助記憶装置１５に格納する。

ステップＳ２０３において、事前情報処理部２１は、データを取得し、時間帯ごとに分割する。具体的には、第１に、事前情報処理部２１は、計測値情報３１（図８）のうちから、ステップＳ２０１の“第２”において受け付けた開始日時から同じく受け付けた終了日時までの期間のデータを切り出す。
第２に、事前情報処理部２１は、ある未処理の特徴量を、ステップＳ２０１の“第３”において受け付けたパタンＩＤに基づき、時間帯ごとに分割する。例えば、事前情報処理部２１がステップＳ２０１の“第３”においてパタンＩＤ“Ｐ０１”を受け付けたとする。この場合、事前情報処理部２１は、ステップＳ２０３の“第１”において切り出したデータを、隙間がなく重複することもない連続する１時間の時間幅を有する時間帯ごとに分割することになる。

ステップＳ２０４において、事前情報処理部２１は、特徴量の基準値を算出する。具体的には、第１に、事前情報処理部２１は、未処理の（時間的に最も古い）時間帯に含まれる複数の特徴量について、平均値μ及び標準偏差σを算出する。
第２に、事前情報処理部２１は、“＋３σ”以上の特徴量及び“−３σ”以下の特徴量（ノイズ又ははずれ値）を削除する。ここでの“＋３”及び“−３”はあくまでも一例である。
第３に、事前情報処理部２１は、当該時間帯に含まれる、ノイズ等を削除した後の複数の特徴量について、平均値μ_２及び標準偏差σ_２を算出し、算出したσ_２を基準値とする。

ステップＳ２０５において、事前情報処理部２１は、マハラノビス距離を算出する。具体的には、事前情報処理部２１は、式１の“μ”に“μ_２”を代入し、“σ”に基準値“σ_２”を代入し、“ｘ_ｉ”に特徴量を代入することによって、個々の特徴量ごとのマハラノビス距離Ｄｍ（ｘ_ｉ）を算出する。

ステップＳ２０６において、事前情報処理部２１は、マハラノビス距離を加重平均し、異常度を算出する。具体的には、第１に、事前情報処理部２１は、重み情報３３（図１０）を参照し、処理対象のマハラノビス距離に乗算される重みを決定する。
第２に、事前情報処理部２１は、式３を使用して異常度を算出する。このとき、事前情報処理部２１は、マハラノビス距離に対し重みを乗算することに代えて、マハラノビス距離を重みで置換してもよい。この場合、異常度は、重み自体の平均となる。

ステップＳ２０７において、事前情報処理部２１は、未処理の時間帯があるか否かを判断する。具体的には、事前情報処理部２１は、未処理の時間帯がない場合（ステップＳ２０７“Ｎｏ”）、ステップＳ２０８に進み、それ以外の場合（ステップＳ２０７“Ｙｅｓ”）、ステップＳ２０４に戻る。戻った先のステップＳ２０４〜Ｓ２０６において、事前情報処理部２１は、未処理の次に古い時間帯について、同様の処理を繰り返す。

ステップＳ２０８において、事前情報処理部２１は、時系列の異常度を生成する。具体的には、事前情報処理部２１は、ステップＳ２０４〜Ｓ２０６の繰り返し処理において算出した異常度に基づいて、データ取込期間における時系列の異常度（時間帯ごとの異常度）を生成する。ここで生成された時系列の異常度の例が、図１５（ａ）である。

ステップＳ２０９において、事前情報処理部２１は、異常発生確率を算出する。具体的には、事前情報処理部２１は、図１５（ａ）のような時系列の異常度に対して図６で説明した方法を適用することによって、異常度の水準ごとに異常発生確率を算出する。

ステップＳ２１０において、事前情報処理部２１は、異常発生確率と異常度との相互関係を決定する。具体的には、第１に、事前情報処理部２１は、ステップＳ２０９において算出した異常度の水準ごとの異常発生確率に基づいて、図６（ｃ）のような、異常発生確率と異常度との相互関係を決定する。
第２に、事前情報処理部２１は、ステップＳ２１０の“第１”において決定した相互関係を滑らかな曲線に近似する。すると、相互関係は、図１５（ｂ）のような曲線になる。図１５（ｂ）から明らかなように、相互関係は、異常度を入力することによって異常発生確率を出力する一種の変換関数である。

ステップＳ２１１において、事前情報処理部２１は、時系列の異常発生確率を算出する。具体的には、事前情報処理部２１は、相互関係を使用することによって、時系列の異常度を時系列の異常発生確率に変換する。ここで変換された時系列の異常発生確率の例が、図１５（ｃ）である。

ステップＳ２１２において、事前情報処理部２１は、未処理の特徴量があるか否かを判断する。具体的には、事前情報処理部２１は、未処理の特徴量がない場合（ステップＳ２１２“Ｎｏ”）、ステップＳ２１３に進み、それ以外の場合（ステップＳ２１２“Ｙｅｓ”）、ステップＳ２０３に戻る。戻った先のステップＳ２０３〜Ｓ２１１において、事前情報処理部２１は、未処理の次の特徴量について、同様の処理を繰り返す。事前情報処理部２１は、ステップＳ２０２の“第１”において受け付けたすべての特徴量についてＳ２０３〜Ｓ２１１の処理を終了すると、当該繰り返し処理を終了する。

ステップＳ２１３において、事前情報処理部２１は、異常発生確率の代表値を算出する。具体的には、事前情報処理部２１は、各診断時点の複数の異常発生確率（電圧の異常発生確率、電流の異常発生確率、温度の異常発生確率、・・・）のうちの最大値を当該診断時点の異常発生確率の代表値とする。事前情報処理部２１がこのようにして時系列の代表値を算出した例が、図１５（ｄ）である。なお、事前情報処理部２１が最大値を代表値とするのはあくまでも一例である。事前情報処理部２１は、例えば、大きい順に所定の数の異常発生確率を平均した値を代表値としてもよいし、単純にすべての特徴量についての異常発生確率の平均値を代表値としてもよい。

ステップＳ２１４において、事前情報処理部２１は、時系列の異常度等を表示する。具体的には、第１に、事前情報処理部２１は、出力装置１３に事前情報表示画面５２（図１６）を表示する。
第２に、事前情報処理部２１は、ユーザが特徴量を選択するのを受け付ける。ユーザが事前情報表示画面５２の“特徴量選択”ボタン５３を押下すると、図１４の特徴量選択欄５１ｂと同じウインドウが表示される。このとき、ステップＳ２０２の“第１”において受け付けた特徴量が強調表示されている。ユーザは、強調表示されている特徴量のうちの１つを選択する。

第３に、事前情報処理部２１は、ユーザがステップＳ２１４の“第２”において特徴量を選択したままの状態で、ボタン５４ａ〜５４ｃのうちのいずれかを押下するのを受け付ける。“時系列の異常度”ボタン５４ａが押下されると、事前情報処理部２１は、事前情報表示画面５２のグラフ表示欄５５に、時系列の異常度（図１５（ａ））を表示する。同様に、“異常度と異常発生確率との相互関係”ボタン５４ｂが押下されると、事前情報処理部２１は、異常度と異常発生確率との相互関係（図１５（ｂ））を表示する。“時系列の異常発生確率”ボタン５４ｃが押下されると、事前情報処理部２１は、時系列の異常発生確率（図１５（ｃ））を表示する。

第４に、事前情報処理部２１は、ステップＳ２１４の“第２”において選択された特徴量を検索キーとして異常原因情報３４（図１１）を検索し、該当したレコード（行）の異常名称又は異常原因を、事前情報表示画面５２の異常原因候補欄５６に表示する。なお、この例の他にも、事前情報処理部２１は、強調表示されている特徴量のそれぞれについてデータ取込期間の異常発生確率の最大値を取得し、その最大値が最も大きい特徴量を検索キーとしてもよい。事前情報処理部２１は、データ取込期間の異常発生確率の最大値が大きい順に所定の数の特徴量を検索キーとしてもよいし、データ取込期間の異常発生確率の最大値が所定の閾値を超える特徴量を検索キーとしてもよい。すると、１又は複数の特徴量（電圧、電流、・・・）について、異常原因等が表示されることになる。

第５に、事前情報処理部２１は、“異常発生確率の代表値”ボタン５４ｄが押下されると、事前情報表示画面５２のグラフ表示欄５５に、時系列の異常発生確率の代表値（図１５（ｄ））を表示する。
事前情報処理部２１は、各情報を画面表示することに代えて又は併せて、各情報の内容を音声出力してもよいし、他の装置（現場責任者に対するメール作成装置等）に対してデータ出力してもよい。
その後、事前情報処理手順を終了する。

（運用情報処理手順）
図１３に沿って、運用情報処理手順を説明する。運用情報処理手順は、以下の相違点を除けば、全体的に事前情報処理手順（図１２）に類似している。
（第１の相違点）事前情報処理手順において異常診断装置は、異常が発生していることが既知である期間を含む過去の期間（例えば数カ月前）の異常発生確率を算出する。一方、運用情報処理手順において異常診断装置は、異常が発生しているか否かが未知である直近の診断対象期間（例えば、現在から遡った２日間）の異常発生予想確率を算出する。
（第２の相違点）事前情報処理手順において異常診断装置は、異常度と異常発生確率との相互関係を新たに生成する。一方、運用情報処理手順において異常診断装置は、相互関係を新たに生成せずに、事前情報処理手順において生成された相互関係を流用する。

図１３の各ステップの処理は、そのステップの下２桁が同じである図１２のステップの処理と概ね同じである。図１３においてステップＳ３０２、Ｓ３０９及びＳ３１０は欠番である。このことは、対応する図１２の処理が図１３では省略されることを意味する。

ステップＳ３０１において、異常診断装置１の運用情報処理部２２は、データ取込期間及びパタンを受け付ける。具体的には、第１に、運用情報処理部２２は、ユーザが入力装置１２に対して開始日時及び終了日時を入力するのを受け付ける。例えば、現時点が２０１７年５月３０日であり、ユーザは、２０１７年５月２８日から２０１７年５月２９日までの特徴量を、診断対象期間における診断対象データとして使用したいとする。この場合、ユーザは、“開始日時”として“２０１７年５月２８日００時００分００秒”を入力し、“終了日時”として“２０１７年５月３０日００時００分００秒”を入力する。
第２に、運用情報処理部２２は、時間幅情報３２（図９）を出力装置１３に表示し、ユーザが任意のパタンＩＤを選択するのを受け付ける。

なお、運用情報処理部２２は、事前情報処理手順のステップＳ２０２の“第１”において受け付けた特徴量が、運用情報処理手順においても引き続き選択されていると看做す。

ステップＳ３０３において、運用情報処理部２２は、データを取得し、時間帯ごとに分割する。具体的には、第１に、運用情報処理部２２は、計測値情報３１（図８）のうちから、ステップＳ３０１の“第１”において受け付けた開始日時から同じく受け付けた終了日時までの期間のデータを切り出す。
第２に、運用情報処理部２２は、ある未処理の特徴量を、ステップＳ３０１の“第２”において受け付けたパタンＩＤに基づき、時間帯ごとに分割する。

ステップＳ３０４において、運用情報処理部２２は、特徴量の基準値を算出する。具体的には、第１に、運用情報処理部２２は、未処理の（時間的に最も古い）時間帯に含まれる複数の特徴量について、平均値μ及び標準偏差σを算出する。
第２に、運用情報処理部２２は、“＋３σ”以上の特徴量及び“−３σ”以下の特徴量（ノイズ又ははずれ値）を削除する。ここでの“＋３”及び“−３”はあくまでも一例である。
第３に、運用情報処理部２２は、当該時間帯に含まれる、ノイズ等を削除した後の複数の特徴量について、平均値μ_３及び標準偏差σ_３を算出し、算出したσ_３を基準値とする。

ステップＳ３０５において、運用情報処理部２２は、マハラノビス距離を算出する。具体的には、運用情報処理部２２は、式１の“μ”に“μ_３”を代入し、“σ”に基準値“σ_３”を代入し、“ｘ_ｉ”に特徴量を代入することによって、個々の特徴量ごとのマハラノビス距離Ｄｍ（ｘ_ｉ）を算出する。

ステップＳ３０６において、運用情報処理部２２は、マハラノビス距離を加重平均し、異常度を算出する。具体的には、第１に、運用情報処理部２２は、重み情報３３（図１０）を参照し、処理対象のマハラノビス距離に乗算される重みを決定する。
第２に、運用情報処理部２２は、式３を使用して異常度を算出する。このとき、運用情報処理部２２は、マハラノビス距離に対し重みを乗算することに代えて、マハラノビス距離を重みで置換してもよい。この場合、異常度は、重み自体の平均となる。

ステップＳ３０７において、運用情報処理部２２は、未処理の時間帯があるか否かを判断する。具体的には、運用情報処理部２２は、未処理の時間帯がない場合（ステップＳ３０７“Ｎｏ”）、ステップＳ３０８に進み、それ以外の場合（ステップＳ３０７“Ｙｅｓ”）、ステップＳ３０４に戻る。戻った先のステップＳ３０４〜Ｓ３０６において、運用情報処理部２２は、未処理の次に古い時間帯について、同様の処理を繰り返す。

ステップＳ３０８において、運用情報処理部２２は、時系列の異常度を生成する。具体的には、運用情報処理部２２は、ステップＳ３０４〜Ｓ３０６の繰り返し処理において算出した異常度に基づいて、データ取込期間における時系列の異常度（時間帯ごとの異常度）を生成する。ここで生成された時系列の異常度の例が、図１７（ａ）である。

ステップＳ３１１において、運用情報処理部２２は、時系列の異常発生予想確率を算出する。具体的には、運用情報処理部２２は、ステップＳ２１０において決定した同じ特徴量についての相互関係を使用することによって時系列の異常度を時系列の異常発生予想確率に変換する。このとき変換された時系列の異常発生予想確率の例が、図１７（ｂ）である。

ステップＳ３１２において、運用情報処理部２２は、未処理の特徴量があるか否かを判断する。具体的には、運用情報処理部２２は、未処理の特徴量がない場合（ステップＳ３１２“Ｎｏ”）、ステップＳ３１３に進み、それ以外の場合（ステップＳ３１２“Ｙｅｓ”）、ステップＳ３０３に戻る。戻った先のステップＳ３０３〜Ｓ３０８及びＳ３１１において、運用情報処理部２２は、未処理の次の特徴量について、同様の処理を繰り返す。運用情報処理部２２は、事前情報処理手順のステップＳ２０２の“第１”において受け付けたすべての特徴量についてＳ３０３〜Ｓ３０８及びＳ３１１の処理を終了すると、当該繰り返し処理を終了する。

ステップＳ３１３において、運用情報処理部２２は、異常発生予想確率の代表値を算出する。具体的には、運用情報処理部２２は、各診断時点の複数の異常発生予想確率（電圧の異常発生予想確率、電流の異常発生予想確率、温度の異常発生予想確率、・・・）のうちの最大値を当該診断時点の異常発生予想確率の代表値とする。運用情報処理部２２がこのようにして時系列の代表値を算出した例が、図１７（ｃ）である。なお、運用情報処理部２２が最大値を代表値とするのはあくまでも一例である。運用情報処理部２２は、例えば、大きい順に所定の数の異常発生予想確率を平均した値を代表値としてもよいし、単純にすべての特徴量についての異常発生予想確率の平均値を代表値としてもよい。

ステップＳ３１４において、運用情報処理部２２は、時系列の異常度等を表示する。具体的には、第１に、運用情報処理部２２は、出力装置１３に運用情報表示画面６２（図１８）を表示する。
第２に、運用情報処理部２２は、ユーザが特徴量を選択するのを受け付ける。ユーザが運用情報表示画面６２の“特徴量選択”ボタン６３を押下すると、図１４の特徴量選択欄５１ｂと同じウインドウが表示される。このとき、事前情報処理手順のステップＳ２０２の“第１”において受け付けた特徴量が強調表示されている。ユーザは、強調表示されている特徴量のうちの１つを選択する。

第３に、運用情報処理部２２は、ユーザがステップＳ３１４の“第２”において特徴量を選択したままの状態で、ボタン６４ａ及び６４ｂのうちのいずれかを押下するのを受け付ける。“時系列の異常度”ボタン６４ａが押下されると、運用情報処理部２２は、運用情報表示画面６２のグラフ表示欄６５に、時系列の異常度（図１７（ａ））を表示する。同様に、“時系列の異常発生予想確率”ボタン６４ｂが押下されると、運用情報処理部２２は、時系列の異常発生予想確率（図１７（ｂ））を表示する。そして同じタイミングで、運用情報処理部２２は、選択されている特徴量についての最も直近の診断時点における異常発生予想確率を異常発生予想確率欄６７に表示する。

第４に、運用情報処理部２２は、“異常発生予想確率の代表値”ボタン６４ｃが押下されると、時系列の異常発生予想確率の代表値（図１７（ｃ））を表示する。そして同じタイミングで、運用情報処理部２２は、最も直近の診断時点における異常発生予想確率の代表値を異常発生予想確率欄６７に表示する。

第５に、運用情報処理部２２は、“異常原因”ボタン６４ｄが押下されると、ステップＳ３１４の“第２”において受け付けた特徴量を検索キーとして異常原因情報３４（図１１）を検索する。そして、運用情報処理部２２は、該当したレコード（行）の異常名称又は異常原因を、運用情報表示画面６２の異常原因候補欄６６に表示する。

なお、この例の他にも、運用情報処理部２２は、強調表示されている特徴量のそれぞれについてデータ取込期間の異常発生予想確率の最大値を取得し、その最大値が最も大きい特徴量を検索キーとしてもよい。運用情報処理部２２は、データ取込期間の異常発生予想確率の最大値が大きい順に所定の数の特徴量を検索キーとしてもよいし、データ取込期間の異常発生予想確率の最大値が所定の閾値を超える特徴量を検索キーとしてもよい。すると、１又は複数の特徴量（電圧、電流、・・・）について、異常原因等が表示されることになる。

運用情報処理部２２は、各情報を画面表示することに代えて又は併せて、各情報の内容を音声出力してもよいし、他の装置（現場責任者に対するメール作成装置等）に対してデータ出力してもよい。
その後、運用情報処理手順を終了する。

本実施形態の異常診断装置は、以下の効果を奏する。
（１）異常診断装置は、過去の異常発生例に基づき、現在の異常発生予想確率を算出することができる。
（２）異常診断装置は、特徴量のマハラノビス距離のうち特に大きなものを確実に評価することができる。
（３）異常診断装置は、過去に異常が実際に発生した頻度に基づき、高い精度で現在の異常発生予想確率を算出することができる。
（４）異常診断装置は、過去の異常発生例に基づき作成した異常度と異常発生確率との相互関係を、現在の診断に活用し、活用結果を出力することができる。
（５）異常診断装置は、特徴量の種類が複数存在する場合であっても、いずれかが示す異常の予兆を知ることができる。
（６）異常診断装置は、拾ったノイズに対して頑強である。
（７）異常診断装置は、異常の原因を出力することができる。

１ 異常診断装置
１１ 中央制御装置
１２ 入力装置
１３ 出力装置
１４ 主記憶装置
１５ 補助記憶装置
２１ 事前情報処理部
２２ 運用情報処理部
３１ 計測値情報
３２ 時間幅情報
３３ 重み情報
３４ 異常原因情報

Claims (9)

1. 診断対象機器から過去における特徴量を時系列で取得し、
   前記取得した特徴量のばらつきを示す異常度を時系列で算出し、
   前記診断対象機器に異常が実際に発生した確率である異常発生確率と、前記異常度との相互関係を決定し、
   前記決定した相互関係に基づいて、前記異常発生確率を時系列で算出する事前情報処理部と、
   前記診断対象機器から診断対象期間における特徴量を時系列で取得し、
   前記取得した診断対象期間における特徴量のばらつきを示す異常度を時系列で算出し、
   前記決定した相互関係に基づいて、前記診断対象期間において前記診断対象機器に異常が発生する確率である異常発生予想確率を時系列で算出する運用情報処理部と、
   を備えることを特徴とする異常診断装置。
2. 前記事前情報処理部及び前記運用情報処理部は、
   所定の時間帯に属する特徴量の値のそれぞれについてマハラノビス距離を算出し、
   前記算出したマハラノビス距離を、その値の大きさに応じて加重平均することによって前記異常度を算出すること、
   を特徴とする請求項１に記載の異常診断装置。
3. 前記事前情報処理部は、
   前記異常度がある水準に達した際前記診断対象機器に異常が実際に発生した頻度、及び、前記異常度が当該水準に達した際前記診断対象機器に異常が実際に発生しなかった頻度に基づいて、前記異常発生確率を算出すること、
   を特徴とする請求項２に記載の異常診断装置。
4. 前記事前情報処理部は、
   前記算出した異常度、前記決定した相互関係、及び、前記算出した異常発生確率を出力し、
   前記運用情報処理部は、
   前記算出した異常度、及び、前記算出した異常発生予想確率を出力すること、
   を特徴とする請求項３に記載の異常診断装置。
5. 前記特徴量は、
   複数の項目を有し、
   前記事前情報処理部は、
   前記異常発生確率を前記項目ごとに算出し、
   前記算出した項目ごとの異常発生確率を使用して所定の方法で異常発生確率の代表値を時系列で決定するとともに出力し、
   前記運用情報処理部は、
   前記異常発生予想確率を前記項目ごとに時系列で算出し、
   前記算出した項目ごとの異常発生予想確率を使用して所定の方法で異常発生予想確率の代表値を時系列で決定するとともに出力すること、
   を特徴とする請求項４に記載の異常診断装置。
6. 前記事前情報処理部及び前記運用情報処理部は、
   所定の時間帯に属する特徴量の値の集合の標準偏差を算出し、
   前記算出した標準偏差に基づいて、特徴量のはずれ値を削除し、
   前記はずれ値を削除した後の特徴量の値の集合の平均値を前記マハラノビス距離を算出するための基準値とすること、
   を特徴とする請求項５に記載の異常診断装置。
7. 前記異常診断装置は、
   前記複数の項目に関連付けて、前記診断対象機器に実際に発生した異常の内容及び原因を記憶している記憶部を備え、
   前記事前情報処理部は、
   所定の基準を満たす程度に前記異常発生確率が高い特徴量の項目に対応する前記異常の内容又は原因を出力し、
   前記運用情報処理部は、
   所定の基準を満たす程度に前記異常発生予想確率が高い特徴量の項目に対応する前記異常の内容又は原因を出力すること、
   を特徴とする請求項６に記載の異常診断装置。
8. 異常診断装置の事前情報処理部は、
   診断対象機器から過去における特徴量を時系列で取得し、
   前記取得した特徴量のばらつきを示す異常度を時系列で算出し、
   前記診断対象機器に異常が実際に発生した確率である異常発生確率と、前記異常度との相互関係を決定し、
   前記決定した相互関係に基づいて、前記異常発生確率を時系列で算出し、
   前記異常診断装置の運用情報処理部は、
   前記診断対象機器から診断対象期間における特徴量を時系列で取得し、
   前記取得した診断対象期間における特徴量のばらつきを示す異常度を時系列で算出し、
   前記決定した相互関係に基づいて、前記診断対象期間において前記診断対象機器に異常が発生する確率である異常発生予想確率を時系列で算出すること、
   を特徴とする異常診断装置の異常診断方法。
9. 異常診断装置の事前情報処理部に対し、
   診断対象機器から過去における特徴量を時系列で取得し、
   前記取得した特徴量のばらつきを示す異常度を時系列で算出し、
   前記診断対象機器に異常が実際に発生した確率である異常発生確率と、前記異常度との相互関係を決定し、
   前記決定した相互関係に基づいて、前記異常発生確率を時系列で算出する処理を実行させ、
   前記異常診断装置の運用情報処理部に対し、
   前記診断対象機器から診断対象期間における特徴量を時系列で取得し、
   前記取得した診断対象期間における特徴量のばらつきを示す異常度を時系列で算出し、
   前記決定した相互関係に基づいて、前記診断対象期間において前記診断対象機器に異常が発生する確率である異常発生予想確率を時系列で算出するする処理を実行させること、
   を特徴とする異常診断装置を機能させるための異常診断プログラム。

Images