JP5337909B2 - 異常検出装置 - Google Patents

異常検出装置 Download PDF

Info

Publication number
JP5337909B2
JP5337909B2 JP2012507951A JP2012507951A JP5337909B2 JP 5337909 B2 JP5337909 B2 JP 5337909B2 JP 2012507951 A JP2012507951 A JP 2012507951A JP 2012507951 A JP2012507951 A JP 2012507951A JP 5337909 B2 JP5337909 B2 JP 5337909B2
Authority
JP
Japan
Prior art keywords
value
series data
probability distribution
model
kth
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2012507951A
Other languages
English (en)
Other versions
JPWO2011121726A1 (ja
Inventor
山 明 弘 酢
田 雄 一 花
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Toshiba Corp
Original Assignee
Toshiba Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Toshiba Corp filed Critical Toshiba Corp
Publication of JPWO2011121726A1 publication Critical patent/JPWO2011121726A1/ja
Application granted granted Critical
Publication of JP5337909B2 publication Critical patent/JP5337909B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B23/00—Testing or monitoring of control systems or parts thereof
    • G05B23/02—Electric testing or monitoring
    • G05B23/0205—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults
    • G05B23/0218—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults
    • G05B23/0224—Process history based detection method, e.g. whereby history implies the availability of large amounts of data
    • G05B23/0227—Qualitative history assessment, whereby the type of data acted upon, e.g. waveforms, images or patterns, is not relevant, e.g. rule based assessment; if-then decisions
    • G05B23/0232—Qualitative history assessment, whereby the type of data acted upon, e.g. waveforms, images or patterns, is not relevant, e.g. rule based assessment; if-then decisions based on qualitative trend analysis, e.g. system evolution
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B23/00—Testing or monitoring of control systems or parts thereof
    • G05B23/02—Electric testing or monitoring
    • G05B23/0205—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults
    • G05B23/0218—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults
    • G05B23/0243—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults model based detection method, e.g. first-principles knowledge model
    • G05B23/0254—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults model based detection method, e.g. first-principles knowledge model based on a quantitative model, e.g. mathematical relationships between inputs and outputs; functions: observer, Kalman filter, residual calculation, Neural Networks
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B2219/00—Program-control systems
    • G05B2219/20—Pc systems
    • G05B2219/24—Pc safety
    • G05B2219/24042—Signature analysis, compare recorded with current data, if error then alarm

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Physics (AREA)
  • Testing And Monitoring For Control Systems (AREA)

Description

本発明は、制御機器の異常検出を行う異常検出装置に関する。
制御機器の異常検出の方法の1つとして、モデルベース異常検出がある。モデルベース異常検出とは、予め監視対象プロセスをモデル化しておき、ある入力に対して、プロセスの出力として観測されたデータと、モデルによって予測された出力との残差をとり、その値に基づいて制御機器の異常度合を出力する方法である。モデルベース異常検出では、未知パラメータを持つ近似モデルを構築し、異常を含まない実績データから未知パラメータを推定するモデル学習機能を含んでいる。
特許第3254624号
モデルベースで異常検出を行う場合、バルブの開度に対するCv特性(流量特性)など、監視対象とするプロセスに用いられる機器特性を、パラメトリックな関数式として定義し、近似モデルを構築する必要が多々ある。
しかしながら、機器特性は、数式として与えられることは稀であるため、新しい機器が導入される毎に、精度が高い近似モデルを与えることはコスト高である。また、各機器には固体差が存在するため、特性グラフからバイアスがかかった形ではずれる問題も存在する。
たとえば特許文献1(特許第3254624号)では、調節弁のスティックスリップの異常検出を行うための近似モデルとして、摺動特性を x’’ = A * x’ + B として定義し、パラメータA,Bを実績データから推定している。
本発明は、個体差のある未知の特性を持つ制御機器の異常検出を高い精度で行うことを可能とした異常検出装置を提供する。
本発明の一態様としての異常検出装置は、
媒体を流通させる媒体経路における異なる箇所に配置され、それぞれ外部から与えられた指示値に従って前記媒体経路内の前記媒体を共通に出力制御する第1〜第K制御機器と、前記第1〜第K制御機器によって出力制御される媒体の出力量を計測する第1〜第Kセンサに対する異常検出装置であって、
(Ci + βi(Xi) ) × (Qi/Qj)^2 − (Cj + βj(Xj) (Xiは第1観測変数、Bi(Xi)は Xiを入力変数とする未知の関数である第1未知特性項、Ciは第1定数項、Qiは第2観測変数、Xjは第3観測変数、Bj(Xj)は Xjを入力変数とする未知の関数である第2未知特性項、Cjは第2定数項、Qjは第4観測変数)によって定義される目的モデルを記憶する目的モデル格納部と、
第1の期間において前記第1〜第K制御機器に与えられた前記指示値と、前記第1〜第Kセンサによる計測値とを含む第1時系列データを所定時間間隔で記憶する第1記憶部と、
前記第1〜第K制御機器における2つの制御機器のペア毎に、前記目的モデルの前記Xiに一方の制御機器の指示値、前記Ciに前記一方の制御機器の配置箇所に応じた未知の定数、前記Qiに前記一方の制御機器に対応するセンサの計測値、前記Xjに他方の制御機器の指示値、前記Cjに前記他方の制御機器の配置箇所に応じた未知の定数、前記Qjに前記他方の制御機器に対応するセンサの計測値を割り当てることによって、それぞれ診断モデルインスタンスを生成し、
前記ペア毎に生成したすべての診断モデルインスタンスのそれぞれにおける前記Bi(Xi)およびBj(Xj)を、前記XiおよびXjの値の範囲をそれぞれ分割したセクション毎に一定の値を取る関数と定義し、前記第1記憶部内の前記第1時系列データを用いて、前記すべての診断モデルインスタンスの値の二乗を最小化するように前記Bi(Xi)、前記Ci、前記Bj(Xj)、前記Cjを同定することにより、最適化された診断モデルインスタンスをそれぞれ得るモデル最適化部と、
前記第1の期間と異なる第2の期間において、前記第1〜第K制御機器に与えられた前記指示値と、前記第1〜第Kセンサによる計測値とを含む第2時系列データを所定時間間隔で記憶する第2記憶部と、
前記第1〜第Kセンサによる最大観測誤差と、前記指示値に対する前記第1〜第K制御機器の最大応答誤差とを表す精度情報を記憶する第3記憶部と、
前記最適化された診断モデルインスタンス毎に、
前記第1記憶部内の前記第1時系列データを用いて前記最適化された診断モデルインスタンスの値である第1残差を求め、前記第1残差の確率分布である第1確率分布を生成し、
前記第2記憶部内の前記第2時系列データを用いて前記最適化された診断モデルインスタンスの値である第2残差を求め、前記第2残差の確率分布である第2確率分布を生成し、
前記最適化された診断モデルインスタンスを生成する元となった2つの制御機器および2つのセンサのそれぞれ毎に、前記第1時系列データにおけるそれぞれに対応する値をそれぞれの前記精度情報に応じて変更し、変更後の時系列データを用いて前記最適化された診断モデルインスタンスの値である第3残差を求め、前記第3残差の確率分布である第3確率分布をそれぞれ生成し、
前記第2確率分布と前記第1確率分布間の距離を表す第1の拡張KL(Kulback-Leibler divergence)情報量と、前記第2確率分布と各前記第3確率分布のそれぞれとの距離を表す第2の拡張KL情報量とを求め、前記第1および第2拡張KL情報量の合計に対する前記第1拡張KL情報量の比率に応じて判定スコアをそれぞれ計算する判定スコア計算部と、
前記第1〜第K制御機器および前記第1〜Kセンサ毎に、前記最適化された診断モデルインスタンスから得られた判定スコアの最小値である総合スコアを求め、第1閾値より大きい総合スコアをもつ制御機器またはセンサに異常が発生したことを判定する総合異常判定部と、
を備える。
本発明により、個体差のある未知の特性を持つ制御機器の異常検出を高い精度で行うことが可能となる。
本発明の実施の形態としての異常検出装置の構成を概略的に示すブロック図である。 本実施の形態を説明するための事例を説明するための図である。 図2と異なる事例を示す図である。 入力画面の一例を示す図である。 インデックステーブルの一例を示す図である。 データテーブルの一例を示す図である。 目的モデル格納部に格納された複数の目的モデル例を示す図である。 モデル生成部の処理の流れを示すフローチャートである。 モデル生成部の出力であるモデルリストの例を示す図である。 モデル情報の例を示す図である。 モデル最適化部による処理の流れを示すフローチャートである。 図11におけるパラメータ同定演算過程におけるメモリ内容を示す図である。 パラメータ同定演算の様子を示す図である。 パラメータ同定演算の様子を示す図である。 パラメータ同定演算の様子を示す図である。 モデル情報の一例を示す図である。 時系列データ記憶部に記憶された時系列データの例を示す図である。 異常検出部よる処理の流れを示すフローチャートである。 異常スコア計算部による異常スコアの計算処理の流れを示すフローチャートである。 残差および信頼度の計算例を説明する図である。 残差および信頼度の計算例を説明する図である。 異常スコアの計算方法を説明する図である。 精度情報に基づき時系列データを補正する例を示す図である。 総合異常判定部による総合異常判定処理の流れを示すフローチャートである。 計算された異常スコア、および総合スコアの計算例を示す図である。 結果表示の画面例を示す図である。 結果表示の画面例を示す図である。 結果表示の画面例を示す図である。 解析結果画面の例を示す図である。
図1は本発明の実施の形態としての異常検出装置の構成を概略的に示すブロック図である。
まず、理解の簡単のため、図1の装置を説明する前に、異常検出を行う対象となる事例について説明する。
図2は、本実施の形態を説明するための事例を説明するための図である。
風または空気(媒体)を送り込むファンをダクト内に持たない空調システムを想定する。ダクト本流が1つ存在し、本流ダクトから複数の支流ダクトに分かれる。本例では3つの支流ダクトに分かれる。これら本流および支流は、本発明の媒体経路の一例に相当する。
各支流ダクトの途中で風量を調節する可変風量装置(VAV：Variable Air Volume)あるいは固定風量装置(CAV: Constant Air Volume)が設置される。
本例ではVAV1〜VAV3が設置される。VAV1〜VAV3はそれぞれ異なる箇所に設置される。VAV1〜VAV3は、各支流ダクトに流れる風量を調整する弁X1〜X3を備える。VAV1〜VAV3は、k (k=1..n) 番目の支流の開度指示値x(k)に応じて、弁X1〜X3の開きを調整するアクチュエータを含む。弁の開きを調整して支流ダクトに流れる風量を外部に出力する。開度指示値は、図示しないコントローラによって生成され、VAV1〜VAV3に与えられる。
またVAV1〜VAV3の近傍に風量センサAF1〜AF3が配置される。各支流ダクトの出口が1つの室内に存在し、k (k=1..n) 番目の支流の風量センサは、風量測定値 Q(k)を観測する。風量センサはVAV（制御機器）によって出力制御された風量（出力量）を観測する。
本装置は、k (k=1..n) 番目の支流のVAVの開度指示値 x(k) と風量測定値 Q(k) を逐次、観測時刻とともに内部の記憶部に記憶し、記憶したデータを元に、アクチュエータ・センサ異常の検出を行う。すなわち、VAV開度指示値x(k)と実際の開度の間に乖離があるか否か(開度異常があるか否か)、風量センサ読み取り値Q(k)と実際の風量の間に乖離があるか否か(風量センサに異常があるか否か)を検出する。
なおダクト構成は図2に限定されず、様々な態様が考えられる。たとえば図3(A)のように支流ダクトの個数は2個以上の任意の個数nでよい。また図3(B)のように各支流ダクトの出口が複数の室内に存在してもよい。図3(B)の例では、各室内は、解放可能なドアまたは窓によって隣接している。
以下、図1は異常検出装置について詳細に説明する。
異常検出装置1は、モデル生成部11、モデル最適化部12、異常検出部13、目的モデル格納部16、基準時系列データ記憶部（第1記憶部）17、モデル格納部18、時系列データ記憶部（第2記憶部）19、結果格納部20、センサ情報記憶部（第3記憶部）21と、条件入力部15、表示部16と、を備える。
要素11〜13は、たとえばプログラムモジュールをCPUに実行させることによって実現される。CPUは、当該プログラムモジュールを、図示しないプログラム記憶部から読み出し、メモリに展開して、実行する。プログラム記憶部は、ハードディスク装置、メモリ装置、CD-R、DVD-R、DVD-RW、フレキシブルディスク等の任意の記憶媒体によって実現できる。
要素16〜21は、ハードディスク装置、メモリ装置、CD-R、DVD-ROM、等の任意の記憶媒体によって構成される。記憶媒体は、異常検出装置1の内部に設けられてもよいし、脱着可能な外部装置でもよい。
条件入力部15は、ユーザによりデータ入力を行うための入力インターフェースであり、キーボード、マウスなどの任意の入力手段によって構成されることができる。
表示部16は、ユーザとの出力インターフェースであり、液晶表示装置、プラズマ表示装置、有機EL表示装置、CRT装置などの任意の表示装置によって構成されることができる。表示部16は、たとえばモデル生成部11または異常検出部13から指定されたデータの表示を行う。表示部がタッチパネル機能を有する場合は、条件入力部15は、当該タッチパネル機能を含んでも良い。
以下、本装置の説明を、モデル生成、モデル最適化、異常検出の順に行う。
[モデル生成]
図4は、条件入力部15による入力を行うための入力画面の一例を示す。
入力画面の左側において、建物としてXビル、フロアとして2階〜14階のそれぞれの部屋A〜Dが選択されている。なお選択するフロア、選択する部屋はそれぞれ1つでもよい。入力画面の右側において、部屋内における異常検出対象となるセンサおよびアクチュエータに対するパラメータ設定を行う。図示の例では、2階〜14階のそれぞれの部屋A〜Dに対して、共通の設定を行っている。ここでは図2の事例を想定し、3つの支流ダクトが存在し、各支流ダクトに配置されたVAV1〜VAV3に対し、VAV開度指示値と、風量をそれぞれ観測パラメータとして選択している。
条件入力部15による以上の入力により、以下の形式を有するポイントオブジェクト名が設定される。
<ビルID.フロアID.ポイント大種別.ポイント小種別.ポイントID.状態量>
ポイントオブジェクト名は、区切り文字(ここでは “.” )によって区分された6つのフィールドを含む。
ビルIDは、選択した建物のIDである。
フロアIDは、選択した階および部屋を特定するIDである。
ポイント大種別は、部屋内での場所種別である。たとえば「I」(Interior)は部屋の内部、「P」(Perimeter)は窓際を表す。
ポイント小種別は、機器種別である。たとえばVAV、CAVがある。
ポイントIDは、ポイント小種別で指定された機器のIDである。
状態量は、当該指定された機器に対して観測される物理量または指示値を示す。
図4の例では、たとえばXビルの12階の部屋Cに対して、1番目の支流ダクト(あるいは1番目のVAV)に対して、
<Xビル.12C.I.VAV.1.開度指示値>と、<Xビル.12C.I.VAV.1.風量>が設定される。
また、2番目の支流ダクト(あるいは2番目のVAV)に対して、
<Xビル.12C.I.VAV.2.開度指示値>と、<Xビル.12C.I.VAV.2.風量>が設定され、
3目の支流ダクト(あるいは3番目のVAV)に対して、
<Xビル.12C.I.VAV.3.開度指示値>と、<Xビル.12C.I.VAV.3.風量>が設定される。
図1の基準時系列データ記憶部（第1記憶部）17は、インデックステーブルと、データテーブルをあらかじめ記憶している。
図5はインデックステーブル171の一例を示す。
インデックステーブル171は、ハンドルIDと、ポイントオブジェクト名から構成される。インデックステーブル171は、たとえば図4の入力画面で設定可能なすべてのポイントオブジェクト名とそのハンドルIDとを記憶している。
図6はデータテーブル172の一例を示す。
データテーブル172は、インデックステーブル171のハンドルIDごとに、該当するポイントオブジェクト名で特定されるセンサ・アクチュエータの観測値の履歴を、基準時系列データ（第1時系列データ）として記憶している。データテーブル171に記憶された基準時系列データ（第1時系列データ）は、いずれもセンサ・アクチュエータが正常であるときの第1の期間に取得されたものである。各ハンドルIDに対応する観測値は、それぞれ一定時間間隔(本例では10分ごと)に、同時に記録されている。
目的モデル格納部16には、ビル空調システムにおいて下記の数式1が適用可能な目的モデルを、複数格納している。目的モデルは診断対象システムにおける部分構造(本例ではビル空調システムにおける任意の2つの制御機器間（あるいは支流ダクト間)）に適用可能なものである。
数式1に含まれる (Ci + βi(Xi) ) × (Qi/Qj)^2− (Cj + βj(Xj) ) が、本発明の目的モデルに対応する。なお「^」はべき乗を意味する。数式1は、目的モデルのCi, Cj, βi, βjを、当該目的モデルの値(e)の二乗(R) を最小にするように求めることを示している。
(数式 1)
min_{Ci, Cj, βi, βj}
R = { (Ci + βi(Xi) ) × (Qi/Qj) ^2 − (Cj + βj(Xj) ) }^2, ただしCi,Cj>0
Ci, Cj ：定数項、βi, βj ：未知特性項、Xi：第1の観測変数、Qi：第2の観測変数,Xj：第3の観測変数,Qj：第4の観測変数
なお数式1をより一般化した以下の式1aを用いてもよい。式1aのf(Y)は第2および第4観測変数の演算項であり、既知の形の関数である。f(Y)はf(Qi,Qj)と記述されてもよい。この場合、(Ci + βi(Xi) ] ) × f(Y) − (Cj + βj(Xj) )が目的モデルに対応する。
(数式 1a)
min_{Ci, Cj, βi, βj}
R = { (Ci + βi(Xi) ] ) × f(Y) − (Cj + βj(Xj) ) }^2
上記の数式1に含まれる目的モデルは任意の2つの支流ダクト内に流れる媒体(流体または気体)のエネルギー保存則をモデル化したものである。このエネルギー保存則は以下の数式1bに示される。図2のような3つのVAVが配置されたダクト構成を例に、目的モデルの導出過程を簡単に説明する。
Figure 0005337909
Ki1,Ki2,Kj1,Kj2は、対象ゾーンのインテリイア/ペリメータエリアに設置されたダクトの圧力損失係数を表す定数である。
i,jは異なる2つの支流を示し、ダクト(またはVAV)のIDである。
Ki1,Kj1は、各支流への分岐点から各支流端に至る管路で定常的に損失する圧力損失係数(支流ダクトの固有の圧力損失係数)である。
Ki2,Kj2は各支流ダクトのVAVの開度に応じて可変的に損失する各支流ダクトの圧力損失係数を示す。
Qi,Qjは風量センサが示す風量計測値(m3/s)である。
vi,vjはダクト風量を調整する弁の開度指示値(0〜1の実数で単位無し)である。
数式1bは、支流iのVAV指示値(VAV開度)とVAV風量、ならびに支流jのVAV指示値(VAV開度)から、支流jのVAV風量を推定できことを意味している。この関係式はVAVの数が変更しても成立するとともに、外乱の影響を受けにくいため汎用性が高い。
これらのパラメータを用いて、ベルヌーイの定理から、任意の2つの支流ダクトに対して、数式1bが導出される。ベルヌーイの定理に基づく数式1bの導出過程は、説明の簡略化のため、割愛する。この数式1bを式変形することで、上述の一般化された目的モデルを導出できる。
Figure 0005337909
図7は、目的モデル格納部16に格納された複数の目的モデル例を示す。
F-01,F02,F03,F04,F05,F06….の目的モデル(エントリ)が格納されている、各目的モデルに含まれる< >内に記述される文字列は、ポイントオブジェクト名に対する抽象表現を意味している。
[s1]、 [s2] 、 [s3] 、 [s4] および [s5]の変数 はそれぞれ任意の文字列に対応づけられる。ただし、1つのエントリに同じ変数が複数存在する場合、各変数は同じ文字列に対応づけられなければならない。たとえばF-01の目的モデルには、［s1］が3つ含まれるが、これらの［s1］はそれぞれ同一の文字列に対応づけられなければならない。
モデル生成部11は、入力情報解析部111とインスタンス生成部112とを備える。モデル生成部11は、これらの要素111，112を用いて、条件入力部15で入力された情報(ポイントオブジェクト名)と、目的モデル格納部16に基づき、数式1のインスタンス生成を行う。
ここにおけるインスタンスとは、数式1の任意の観測変数Xi, Xj, Qi,Qj に対して、具体的な観測変数を割り当てることを意味する。たとえば数式1のXiに “Xビル.2A.I.VAV.1.開度指示値” 、Qiに “Xビル.2A.I.VAV.1風量” を割り当てるなどである。
以下、図8を参照して、モデル生成部11の処理の詳細を説明する。
図8は、モデル生成部11の処理の流れを示すフローチャートである。
モデル生成部11は、条件入力部15から利用者が情報を入力したことをトリガーとして処理を開始する。
入力情報読込部111が、条件入力部15により登録された情報(複数のポイントオブジェクト名)と、基準時系列データ記憶部17のインデックステーブル171をメモリに読み込む(S11)。
各ポイントオブジェクト名における上位3つの情報(ビルID、フロアID、ポイント大分類ID)を抽出し、抽出した上位3つの情報を含む要素のリスト(部分構成リスト)を作成する(S12)。
図4の入力事例の場合、{“Xビル.2A.I”, “Xビル.2B.I”, …, “Xビル.14C.I”, “Xビル.14D.I”} のような52個の要素を含む部分構成リストが作成される。
次に、部分構成リストにおける各要素に対して、以下のステップS13〜S16を実行する。
ステップS13で、部分構成リストから、まだ選択されていない要素を選択する。
ステップS14では、ステップS13で選択された要素に対して、診断モデルインスタンスを生成する(S14)。
たとえば選択された要素 が“Xビル.2A.I”であるとすると以下のようになる。
条件入力部15から入力されたポイントオブジェクト名のうち、この要素“Xビル.2A.I”を含むポイントオブジェクト名を列挙すると、以下の6個のポイントオブジェクト名が得られる。
Xビル.2A.I.VAV.1.開度指示値
Xビル.2A.I.VAV.1.風量
Xビル.2A.I.VAV.2開度指示値
Xビル.2A.I.VAV.2.風量
Xビル.2A.I.VAV.3.開度指示値
Xビル.2A.I.VAV.3.風量
これらのポイントオブジェクト名を、目的モデル格納部16に登録されている各モデルと照合すると、下記のF-01のモデルにマッチすることが分かる。照合は文字列マッチングで行えばよい。つまり、[s1],[s2]等の各変数に矛盾が発生しない(同じ変数には同じ文字列が入る)ようなモデルを選択する。
F-01: (C1 + β1(<[s1].[s2].[s3].VAV.[s4].開度指示値>))
*(<[s1].[s2].[s3].VAV.[s4].風量>/<[s1].[s2].[s3].VAV.[s5].風量>) 2 -(C2 + β2(<[s1].[s2].[s3].VAV.[s5].開度指示値>))
このF-01のモデルにおいて、[s1]〜[s3]は一意に決まるが、[s4], [s5]に関しては6通りの組合せが発生する。
具体的に、[s1]は“Xビル”、[s2]は“2A”、[s3]は“I” で一意に決まるが、[s4], [s5]の組み合わせは、以下の6通り存在する。なお、［s4］,[s5]はいずれもポイントIDに対応する。
([s4], [s5]) =(1, 2), (2, 1), (1, 3), (3, 1), (2, 3), (3, 2)
よって、F-01のモデルに対して、下記の6通りのインスタンス(診断モデルインスタンス)<1>〜<6>が得られる(図9参照)。
Figure 0005337909
ステップS15では、ステップS14で生成した各診断モデルインスタンスに対応して、目的モデルに照合するための部分診断リスト(パラメータテーブルとアクセステーブル)を登録する。
ここでパラメータテーブルは、各診断モデルインスタンスの係数(Ci, Cj, βi, βj)を格納したテーブルである。パラメータテーブルの例を図9の右に示す(「PARMS」属性にリンクされている)。
たとえば上記表の最下のインスタンス<6>に対しては、F-01の係数C1として係数C［3］、F-01の係数C2として係数C［2］、F-01の係数β1として係数β［3］、F-01の係数β2としてβ[2]が格納される。すなわち、当該最下のインスタンスに対して、F-01のC1,C2,β1,β2(あるいは数式1内のCi,Cj,β2, β2)に対応して、“C[3] ,C[2],β[3],β[2]”がパラメータテーブルに登録される。なお現時点ではこれらのパラメータの値は同定されていない(未知である)。
同様にして他の残りの5つのインスタンス<1>〜<5>に対しても同様にパラメータが登録される。
アクセステーブルは、各診断モデルインスタンスのそれぞれに含まれる観測変数(開度指示値、風量)に対応するIDを格納する。アクセステーブルの例が図9の右に示される(「ACCESS」属性にリンクされている)。
アクセステーブルは、基準時系列データ記憶部17に記憶された基準時系列データへのマッピング、および後述する異常診断において用いる時系列データ記憶部19に記憶された時系列データへのマッピングのために用いられる。
たとえばインスタンス<6>の場合、
<Xビル.12C.I.VAV.3.開度指示値>
<Xビル.12C.I.VAV.2.開度指示値>
<Xビル.12C.I.VAV.3.風量>
<Xビル.12C.I.VAV.3.風量>
の4つの観測変数が存在し、
それぞれ対応するハンドルIDは、135,134,145,144である(図5のインデックステーブル171を参照)。したがって、上記インスタンス<6>に対して、“135,134,145,144”がアクセステーブルに登録される。
以上に説明したステップS14,S15の処理を、部分構成リストの残りのすべての要素についても同様に行う。
部分構成リストに含まれるすべての要素についての処理が終了すると、ビルの階およびフロア毎のモデル情報を含むモデルリストが、モデル格納部18に格納される(ステップS17)。
モデルリストはモデル生成部の出力であり、図9に示される。
図9の右には、「12C」(12階の部屋C)に対して生成されたモデル情報の例が取り出して示されている。他の残りの階および部屋についても同様の形式のモデル情報が存在する。これらのモデル情報の集合がモデルリストに相当する。
図9に示すように、モデルリストは、条件入力部で選択された各階および各部屋のモデル情報を含み、1〜Nmの番号がモデル情報の番号である。ビルIDとフロアIDの組ごとに番号が設定され、モデル情報が生成されるごとに順次1〜Nmの番号が増える。図10に「2A」(2階の部屋A)に対して生成されたモデル情報の例を示す。「2A」のモデル情報は上記の処理フローで最初に作成されるため番号1が付与され、この時点ではまだ他のモデル情報は存在しない。「12C」に対するモデル情報の生成は39番目に行われ、番号39が付与されている(図9参照)。なお、本例では最後の番号Nmは52となる。
ここでモデル情報はそれぞれ以下の属性を含む。
“ID”、“Fig”、“BIND”、“PARMS”、“ACCESS”、“C”、“segment”,“B”,“w”
たとえば、Xビルの12階の部屋Cのモデル情報の番号は39番であり、“ID”は、モデル情報の識別子“Xビル.12C”へのリンクを示す。
“Fig”は、条件入力部15による入力画面の画面データが保存されたファイルへのリンクを示す。
“Bind”は、前述した診断モデルインスタンスへのリンクを示す。
“PARMS”は、先に少し説明したが、パラメータテーブルにリンクされ、“ACCESS”がアクセステーブルにリンクされる。
属性“C”(定数項),“segment”(セグメント)、“B”(未知特性項)、“w”(信頼値)に対しては現時点ではまだ何にも設定されていない。詳細は後述する。
なお図4の入力画面ではポイント大種別が“I”の場合を示したが、1つの部屋に“I”と“P”が混在する場合もある。この場合は、“I”と“P”それぞれごとに別個に処理が行われる(上述したように目的モデルへのマッチングは、ポイントオブジェクト名における上位3つの情報(ビルID、フロアID、ポイント大分類ID)により行われる)。したがって、この場合、1つの部屋(たとえば12階の部屋C)に対して、モデル情報が、“I”と“P”のそれぞれごとに生成される。
[モデル最適化]
図1のモデル最適化部12は、データ読込部121、未知特性項初期演算部122、定数項収束演算部123、未知特性項収束演算部124を備える。
モデル最適化部12は、モデル情報毎に、基準時系列データ記憶部17に記憶された基準時系列データに対して、診断モデルインスタンスの値の二乗(数式1のR)を最小化するように当該診断モデルインスタンスのパラメータを学習する。
データ読込部121は、基準時系列データ記憶部17に記憶された基準時系列データをメモリに読み込む。
未知特性項初期演算部122は、各診断モデルインスタンスに含まれる定数項を固定するとともに未知特性項を共通とする条件下で、メモリに読み込んだ基準時系列データに対して誤差(診断モデルインスタンスの値)の二乗が最小になるように、未知特性項を同定し、メモリに記憶する。
定数項収束演算部123では、現在メモリに記憶されている状態で未知特性項を固定した上、メモリに読み込んだ基準時系列データに対して誤差の二乗が最小になるように定数項を同定しメモリに記憶する。
未知特性項収束演算部124は、現在メモリに記憶されている状態で定数項を固定した上、メモリに読み込んだ基準時系列データに対して誤差の二乗が最小になるように、未知特性項を同定し、メモリに記憶する。
定数項収束演算部123と未知特性項収束演算部124は、終了条件を満たすまで交互に処理を繰り返し実行し、終了条件を満たした段階で、そのときのパラメータをモデル格納部18に格納する。このように求められたパラメータを設定した診断モデルインスタンスは、最適化された診断モデルインスタンスと称される。
図11は、モデル最適化部12による処理の流れを示すフローチャートである。モデル最適化部12は、モデル生成部11の後に実行される。したがって、モデル最適化部12の実行前には、Xビルのユーザ選択された階・部屋に対するモデル情報がモデル格納部18に記憶されている。
まず、データ読込部21により以下のステップS21〜S23の処理を実行する。
モデル格納部18に記憶されたモデルリストにおいて任意の番号のモデル情報をメモリに読み込む (S21)。たとえば39番の「Xビル、12階、部屋C」のモデル情報を読み込む。
読み込んだモデル情報に含まれる診断モデルインスタンス毎に、ACCESS属性を利用して、基準時系列データ記憶部17の基準時系列データから該当するデータ(ACCESS属性のハンドルIDに対応する観測値と、タイムスタンプ(時刻情報))をメモリに読み込む(S22)。
ステップS22で読み込んだデータに対して欠損値処理を行う(S23)。欠損値処理は、既存の方法で行う。たとえば、入力されたデータに対して大きな値変更をすることがなく、欠損値や外れ値の除去能力が比較的高く、かつ計算量が少ない修正メジアンフィルタを適用することが望ましい。
次に、未知特性項初期演算部122、定数項収束演算部123および未知特性項収束演算部124によりステップS24〜S28を行う。
ステップS24は未知特性項初期演算部122により行う。まず、未知特性項であるβk(X) は既定の関数で与えられていないため、βk(X)を、α[k, s(X)] と置き換える。α[k, s(X)]は、Xの取り得る値の範囲を一定幅で分割したM個のセグメントのそれぞれ毎に定数を与える関数である。
すなわち、s(X)は、Xが与えられた場合に、Xが属するセグメントのインデックスを返す関数であり、与えられたXに対して1〜Mのいずれかの自然数を出力する(M：セグメント数)。kは、数式1の例ではiまたはjに相当する未知特性項識別子である。したがって、α[k, s(X)]は、Xが与えられたとき、Xが属するセグメントの定数を返す。
上記の置き換えを用いることにより、最適化を示す式である数式1は、次の数式2で表すことができる。
(数式 2)
min_{C[i], C[j], α[i,1], …,α[i,M], α[j,1], …, α[j,M]}
R = { (C[i] + α[i, s(Xi) ] ) × (Qi/Qj)^2 − (C[j] + α[j, s(Xj) ] ) }^2, i≠j
この数式2は Rを最小化するようにC[i], C[j], α[i,1], …,α[i,M], α[j,1], …, α[j,M] の計(2×M + 2)個のパラメータを求める(最適化する)ことを意味する。C[i]>0, C[j]>0を制約条件とする。
この問題を一度に求解するのは困難である。このような問題を解く場合の1つとして、固定する目的変数と、最適化する目的変数とを交互に入れ替えながら反復演算で求解する方法が知られている。ここでは、以下の数式3、数式4、数式5による3つのステップの反復演算によって、数式2を求解する。
(数式 3)
min_{α[0,1], …, α[0,M] }
R = { (C0 + α[0, s(Xi) ] ) * (Qi/Qj)^2 − (C0 + α[0, s(Xj) ] ) }^2

(数式4)
min_{C[i], C[j]}
R = { (C[i] + α[i, s(Xi) ] ) * (Qi/Qj)^2 − (C[j] + α[j, s(Xj)] ) }^2

(数式5)
min_{α[i,1], …,α[i, M], α[j, 1], …, α[j, M]}
R = { (C[i] + α[i, s(Xi) ] ) * (Qi/Qj)^2 − (C[j] + α[j, s(Xj) ] ) }^2

(ステップ1)1つめのステップは、未知特性項初期演算部122において、数式3のインスタンスである下記の数式3.1を求解する(S25)。
(数式3.1)
min_{α[0,1], …, α[0,M]}
R = { (C0 + α[0, s(xi) ] ) × (Qi/Qj)^2 − (C0 + α[0, s(xj) ] ) }^2, i≠j
以下、事例に従って (数式3.1)を説明する。
3つの各支流の固有のダクト圧力損失係数に相当するC1, C2, C3 が同じ値C0( = 1)を持つと仮定する。すなわち各支流のダクトの長さ・形状・材質が同じと仮定する。なおC1, C2, C3は、VAV１〜３が配置された箇所（ダクト）に応じた未知の定数項であるといえる。
さらに、VAV1〜VAV3のそれぞれの開度指示値から推定されるダクト圧力損失係数に相当するα[1,1:M], α[2,1:M], α[3,1:M] も、機器ごとの個体差がなく共通の支流ダクト圧力損失係数α[0,1:M] (すなわちα[0,1], …, α[0,M])を持つと仮定する。つまり、α[1,1], α[2,1], α[3,1]はいずれも同じ値α[0,1]をもち、α[1,2], α[2,2], α[3,2]はいずれも同じ値α[0,2]をもち、・・・α[1,M], α[2,M], α[3,M]はいずれも同じ値α[0,M]をもつと仮定する。
このとき、基準時系列データ記憶部17から読み込んだデータに対する二乗誤差が最小となるように、共通の支流ダクト圧力損失係数に相当するα[0,1:M]を求める。

具体的には、基準時系列データ記憶部17から読み込んだデータを、 診断モデルインスタンス<1>〜<6>毎に、数式3.1に代入することで、下記のように、データ数(=n)×6の非線形連立方程式(数式3.2)が立てられる。これを既存のソルバで解くことで、α[0,1:M] (=α[0,1], …,α[0,M])が求まる(S24〜S25)。すなわち各式の左辺が同時にできるだけゼロに近づくように最適化することで、α[0,1:M] (=α[0,1], …,α[0,M])が求まる。数式3.2において、R[a,b]の“a”は、診断モデルインスタンスを識別する番号であり、“b”はn個の各データを識別する番号である。
(数式3.2)
{ (C0 + α[0, s(x1[1]) ] ) × (Q1[1] / Q2[1] )^2 − (C0 + α[0, s(x2[1]) ] ) }^2 −R[1,1] = 0
…
{ (C0 + α[0, s(x1[n]) ] ) × (Q1[n] / Q2[n] )^2 − (C0 + α[0, s(x2[n]) ] ) }^2 − R[1,n] = 0
{ (C0 + α[0, s(x1[1]) ] ) × (Q1[1] / Q3[1] )^2 − (C0 + α[0, s(x3[1]) ] ) }^2 − R[2,1] = 0
…
{ (C0 + α[0, s(x1[n]) ] ) × (Q1[n] / Q3[n] )^2 − (C0 + α[0, s(x3[n]) ] ) }^2 − R[2,n] = 0
{ (C0 + α[0, s(x2[1]) ] ) × (Q2[1] / Q1[1] )^2 − (C0 + α[0, s(x1[1]) ] ) }^2 − R[3,1] = 0
…
{ (C0 + α[0, s(x2[n]) ] ) × (Q2[n] / Q1[n] )^2 − (C0 + α[0, s(x1[n]) ] ) }^2 − R[3,n] = 0
・ ・・・
{ (C0 + α[0, s(x3[1]) ] ) × (Q3[1] / Q2[1] )^2 − (C0 + α[0, s(x2[1]) ] ) }^2 − R[6,1] = 0
…
{ (C0 + α[0, s(x3[n]) ] ) × (Q3[n] / Q2[n] )^2 − (C0 + α[0, s(x2[n]) ] ) }^2 − R[6,n] = 0
図13に、本ステップ1(S25)の演算の様子を示す。
図13の左に示すように、C1, C2, C3 が同じ値C0( = 1)として固定され、またα[0,1:M] (=α[0,1], …,α[0,M])の初期値を1としている。そして、上記数式3.2を解くことで、図13の右のように、α[0,1:M] (=α[0,1], …,α[0,M])が求められている。X軸がVAV開度指示値で、Y軸がα[0,s(xk)]+C0である。
図12は、ステップ1〜3(S25〜S27)の演算過程におけるメモリ内容を示すもので、最適化する変数はメモリ上で配列テーブルとして保持されている。本ステップ1のメモリ内容は図12の上に示される。白抜きにしているセルの変数を固定しながら、網掛けにしているセルの変数を最適化している。本ステップ1では上述のようにC1, C2, C3 が同じ値C0( = 1)を持つとして固定される(白抜きセル)。α[1,1:M], α[2,1:M], α[3,1:M]は共通のα[0,1:M]を持つと仮定した上で、α[0,1:M]が最適化される(網掛けのセル)。
(ステップ2)2つめのステップは、定数項収束演算部123において、数式4のインスタンスである下記の数式4.1を求解する(S26)。
(数式4.1)
min_{C[i], C[j]}
R = { (Ci + α[i, s(xi) ] ) ×(Qi/Qj)^2 − (Cj + α[j, s(xj) ] ) }^2, i≠j
以下、事例に従って(数式4.1)を説明する。
VAV1〜VAV3の開度に対する圧力損失係数に相当するα[1,1:M], α[2,1:M], α[3,1:M]を現在の状態で固定(前述のように最初は、α[1,1:M], α[2,1:M], α[3,1:M]はすべて共通のα[0,1:M])し、基準時系列データ記憶部17から読み込んだデータに対する二乗誤差が最小となるように、各支流のダクト圧力損失係数に相当する C1, C2, C3を求める。
ステップ1と同様に、診断モデルインスタンス<1>〜<6>毎に、読み込んだデータを入力することで、データ数(=nとする)×6の連立方程式が立てられ、これを既存のソルバで解くことでC1,C2,C3が求まる (S26)。
図14に、本ステップ2(S26)の演算の様子を示す。
図左のC1old,C2 old,C3 oldは、本ステップによる最適化前のC1,C2,C3であり、図右のC1new,C2 new,C3 newは、本ステップによる最適化後のC1,C2,C3である。
各支流ダクトに対応するα[k,s(xk)](すなわちα[1,s(x1)]、α[2,s(x2)]、α[3,s(x3)])は、本ステップでは固定される。
なお、X軸は、VAV開度指示値であり、具体的にx1はVAV1の開度指示値、x2はVAV2の開度指示値、x3はVAV3の開度指示値である。Y軸は、α[k,s(x)]+Ckである。
本ステップ2(S26)でのメモリ内容は図12の上から2番目に示される。白抜きにしているセルの変数を固定しながら、網掛けにしているセルの変数を最適化している。本ステップ2では上述のようにα[1,1:M], α[2,1:M], α[3,1:M]を固定しつつ、C1,C2,C3を最適化している。
(ステップ3)3つめのステップは、未知特性項収束演算部124において、数式5のインスタンスである下記の数式5.1を求解する。
(数式5.1)
min_{α[i,1], …,α[i,M], α[j,1], …, α[j,M]}
R = { (C[i] + α[i, s(xi) ] ) ×(Qi/Qj)^2 − (C[j] + α[j, s(xj) ] ) }^2, i, j=1,..,3, i≠j
以下、事例に従って(数式5.1)を説明する。
各支流のダクト圧力損失係数に相当する C1, C2, C3を現在の状態(S26で求めた値)で固定し、基準時系列データ記憶部17から読み込んだデータに対する二乗誤差が最小となるように、VAV1〜VAV3の開度指示値に対するダクト圧力損失係数に相当するα[1,1:M] (= α[1,1], …,α[1,M])、α[2,1:M](=α[2,1], …, α[2,M])、α[3,1:M](=α[3,1], …, α[3,M])を求める。具体的には、これまで同様に、診断モデルインスタンス<1>〜<6>毎に、読み込んデータを入力することで、データ数(=nとする)×6の連立方程式を立て、これを既存のソルバで解くことで、α[1,1:M]、α[2,1:M]、α[3,1:M]が求まる (S27)。
図15に、本ステップ3(S27)の演算の様子を示す。
図左のαold [1,s(x1)]、αold [2,s(x2)]、αold [3,s(x3)]は、本ステップによる最適化前のα[1,s(x1)]、α[2,s(x2)]、α[3,s(x3)]であり、αnew [1,s(x1)]、αnew [2,s(x2)]、αnew [3,s(x3)]は本ステップによる最適化後のα[1,s(x1)]、α[2,s(x2)]、α[3,s(x3)]である。
各支流ダクトに対応するC1,C2,C3は、本ステップでは固定される。
なお、図右の網掛けのブロックは最適化の前後でα[k,s(x)]の値が変化したことを示し、網掛けされていないブロックは最適化の前後でα[k,s(x)]の値が変化しなかったことを示す。
本ステップ3(S27)でのメモリ内容は図12の下から2番目に示される。白抜きにしているセルの変数を固定しながら、網掛けにしているセルの変数を最適化している。本ステップ3では上述のようにC1,C2,C3を固定しつつ、α[1,1:M], α[2,1:M], α[3,1:M]を最適化している。
ここでステップS27では後述の異常判定の信頼性向上のため以下の処理を行っても良い。
ステップS22でメモリに読み込んだデータに含まれるVAV開度指示値の同時分布に疎な部分がある場合、疎な部分のデータにより求められる部分解 α[k,s] は、真に求めるべき値に対して誤差が大きくなる可能性がある。誤差が大きくなると、後述の異常検出部13における異常判定に誤りが生じやすくなる。そこで、後の異常判定においてαの信頼性に基づいた処理を行うことができるように、α[k,s]に対する信頼値 w[k,s]を設定する。
具体的に、α[1,1:M]( =α[1,1]、α[1,2] , …,α[1,M])に対し、信頼値w[1,1:M](= w [1,1]、w [1,2] , …, w [1,M]))を計算し、α[2,1:M]( =α[2,1]、α[2,2] , …,α[2,M])に対し、信頼値w[2,1:M](= w [2,1]、w [2,2] , …, w [2,M]))を計算し、α[3,1:M]( =α[3,1]、α[3,2] , …,α[3,M])に対し、信頼値w[3,1:M](= w [3,1]、w [3,2] , …, w [3,M]))を計算する。
この際、区間s(=1,2,…,M)に含まれるデータ数が多いほど、信頼値が高くなるようにする。たとえばデータ数が多いほど高い信頼値を与える関数またはテーブルを用意しておき、当該関数またはテーブルに、該当区間に含まれるデータ数を与えることにより、当該区間の信頼値を計算する。
本ステップS27で信頼値の計算を行う場合のメモリ内容が図12の下に示される。なお信頼値の計算は、一度行えばよいため、ステップS27の実行のたびに、繰り返し行う必要はない。
以上に説明したステップ2，3(図11のS26〜S27)を、終了条件を満たすまで繰り返し実行することで(S28)、環境に対応した、固有の定数項Ck (k=1,2,3)(例：支流ダクトの固有の圧力損失係数)と、機器特有の未知特性項α[k, s(xk)](k=1,2,3, s(xk)=1,2,3,…M) (例：VAV開度指示値から推定される圧力損失係数)を求めることができる。
ここで終了条件は、利用したソルバに応じて決定されるものであるが、たとえば、「解に収束」「目的変数の変化量が許容範囲より小さい」「繰り返し計算回数が最大数に到達」が挙げられる。
終了条件が満足されたら、求められたモデルパラメータ情報として、Ck (k=1,2,3)、α[k, s(xk)](k=1,2,3, s(xk)=1,2,3,…M))を、モデル格納部18に記憶させる(S29)。
具体的に、モデル格納部18のモデル情報の属性“C”(定数項),“segment”(セグメント)、“B”(未知特性項)、“w”(信頼値)にリンクさせて、Ck、s(xk), α[k, s(xk)],w[k, s(xk)]をそれぞれ記憶させる。
ステップS29でパラメータ情報が追加された後のモデル情報の一例を図16に示す。
[異常検出]
時系列データ記憶部(第2時系列データ記憶部)19は、異常検出の対象とする機器の状態量(指示値、風量など)を、一定時間間隔で順次、時系列データ(第2時系列データ)としてリルタイムに記憶する。具体的に、図5に示したような各ハンドルIDに対応する状態量を逐次記憶する。時系列データ記憶部19に記憶された時系列データの例を図17に示す。時系列データ記憶部19がデータを記憶する期間は、基準時系列データ記憶部17が記憶する第1の期間と異なる第２の期間である。第２の期間は第1の期間より後でもよいし、前でもよい。
センサ情報記憶部21は、センサ・アクチュエータの精度情報を記憶する。精度情報は、たとえば工場出荷でメーカが保証している精度を用いる。本例では風量センサ毎の最大検出誤差が、たとえば±5％、または±0.5(m/s)などとして、センサ情報記憶部21に記憶させている。また開度指示値に対するアクチュエータの最大応答誤差がセンサ情報記憶部21に記憶されている。アクチュエータの最大応答誤差は、指示値に対して最大でどの程度、弁の実動作開度に誤差が生じるかを示す。なおアクチュエータの精度情報は、開度指示値での表現に換算して、センサ情報記憶部21に記憶されてもよい。たとえば弁度の開度誤差範囲に対応する開度指示値範囲を求め、求めた範囲を精度情報として記憶してもよい。
異常検出部13は、モデル読込部131、データ選択部132、異常スコア計算部(判定スコア計算部)133、総合異常判定部134、異常検出提示情報生成部135を備える。
異常検出部13は、時系列データ記憶部19から診断対象とする任意の期間分の時系列データを読み込み、読み込んだ時系列データから制御機器(アクチュエータ)・センサに異常があるか否かを判定する。そして、判定結果を結果格納部20へ格納するとともに、表示部16へ表示する。
これまで説明したモデル生成部11、モデル最適化部12が、異常検出装置1の導入直後や、異常検出の対象となる機器を取り巻く環境に変化があった場合に実行されるブロックであるのに対し、異常検出部13は、常時、もしくは、定期間隔(例えば毎日1回など)で実行される。異常検出部13の処理の実行前には、モデル生成部11およびモデル最適化部12によって生成および更新されたモデル情報がモデル格納部18に格納されていることが前提となる。
以下では、1日1回、VAV開度とVAV風量の異常診断を実行する場合を例にして説明を行う。
図18は、異常検出部13における要素131〜134による処理の流れを示すフローチャートである。
まず、モデル読込部131がモデル格納部18に記憶されたモデル情報を読み込む(S301)。より詳細に、ビルIDをキーにして、モデル格納部18に格納されているモデル情報(各階の各部屋)を読み込む。モデル情報毎に、以降のステップを行って、異常診断を行う。説明の簡単のため、12階の部屋Cのモデル情報を想定して説明を行う。
データ選択部132が、モデル情報の診断モデルインスタンス毎に、ACCESS属性を利用して、時系列データ記憶部19から、対応する時系列データ(ACCESS属性に対応する値と時刻)をメモリに読み込む(S303)。
データ選択部132は、読み込んだ時系列データに対して欠損値処理を行う(S304)。欠損値処理は既存の方法を用いればよい。読み込んだ時系列データに対して大きな値変更をすることがなく、欠損値や外れ値の除去能力が比較的高く、かつ計算量が少ない修正メジアンフィルタを適用することが望ましい。本ステップは省略することも可能である。
次に、異常スコア計算部133による異常スコア(判定スコア)の計算を行う(S305)。この処理を、当該モデル情報における診断モデルインスタンス毎に繰り返し実行する。
図19は異常スコア計算部133による異常スコアの計算処理の流れを示すフローチャートである。
まず、最初の診断モデルインスタンスに対応するPRMS属性、ACCESS属性に基づき、該当するパラメータ情報と、時系列データの時刻毎のデータとを参照し、数式6(残差の計算式)、および数式7(信頼度の計算式)を計算することで、残差ベクトルと信頼度ベクトルを求める(S305a,S305b)。以下図20および図21を用いて、詳細を説明する。
(数式6)
e = ( C[i] + α[i, s(xi) ] ) ×( Qi / Qj )^2 − ( C[j] + α[j, s(xj) ] )
(数式7)
r= w[i,s(xi)] × w[j, s(xj)]
たとえば1番目の診断モデルインスタンス<1>の場合、まず時系列データの先頭データに対して、図20に示す参照を行う。2010/8/22/00:10の時系列データ(診断期間の一番始めの時系列データ)のうち、ACCESS属性に示されるID133,134,143,144の値が参照される。ID133はVAV1の開度指示値、ID 134はVAV2の開度指示値、ID 143はVAV1の風量、ID 144はVAV2の風量である(図5参照)。
またVAV1、VAV2の定数項は、C1=2.31、C2=1.99である。
また、VAV1の未知特性項B1(α[1, s(x1) ])、VAV2の未知特性項B2(α[2, s(x2) ])の値は、それぞれ、0.43、0.44である。すなわち、VAV1の観測変数値(開度指示値)は、0.71(ID133の時系列データ値)であり、当該0.71はsegment属性のセグメント「0．8」に属するから、α[1, s(0.71) ]=α[1, 0.8 ]=0.43となる。同様にVAV2の観測変数値(開度指示値)のは、0.71(ID134の時系列データ値)であり、当該0.71はsegment属性のセグメント「0．8」に属するから、α[2, s(0.71) ]=α[2, 0.8 ]=0.44となる。
また各未知特性項B1,B2の信頼値は、segment属性のセグメント「0．8」に対応する0.98である。
また風量センサAF1,AF2の風量値Q1,Q2は、0.20,0.21(それぞれID143,144の時系列データ値)である。
したがって数式6により、診断モデルインスタンス<1>の残差は以下のように計算される。
e = ( C[i] + α[i, s(xi) ] ) ×( Qi / Qj )^2 − ( C[j] + α[j, s(xj) ] )
= ( C[1] + α[1, s(x1) ] ) ×( Q1 / Q2 )^2 − ( C[2] + α[2, s(x2) ] )
= (2.31 + α[1, s(0.71) ] ) ×( 0.20 / 0.21 )^2 − ( 1.99 + α[2, s(0.71) ] )
= (2.31 + α[1, 0.8 ] ) ×( 0.20 / 0.21 )^2 − ( 1.99 + α[2, 0.8 ] )
= (2.31 + 0.43 ) ×( 0.20 / 0.21 )^2 − ( 1.99 + 0.44 )
=0.06
また数式7により、診断モデルインスタンス<1>の信頼度は以下のように計算される。
r= w[i,s(xi)] × w[j, s(xj)]
= w[1,s(x1)] × w[2, s(x2)]
= w[1,s(0.71)] × w[2, s(0.71)]
= w[1,0.8] × w[2, 08]
= 0.98 × 0.98
=0.96
以上に示した数式6および数式7の計算を、時系列データの2番目以降の時刻のデータについても同様にして行う。2番目の時刻のデータの場合の参照を示すと、図21のようになる。この場合、数式6によりe=-0.03、数式7によりr=0.88が計算される。
以上の処理により、診断モデルインスタンス<1>に対し、時系列データに含まれるデータ数(時刻数)分の残差eを並べた残差ベクトルe<1>と、同個数の信頼度rを並べた信頼度ベクトルr<1>、以下のように得られる。
残差ベクトルe<1>=(e1,e2,e3,…..en) n:はデータの個数
信頼度ベクトルr<1>=(r1,r2,r3,…,rn)
なお<>内の数値は、1番目の診断モデルインスタンス<1>に対応することを示す。

次にステップS305cにおいて、信頼度ベクトル内の各要素の値を検査し、所定の信頼閾値以上の値の要素をすべて特定する。そして特定した要素と同じ位置(時刻)に対応する残差を残差ベクトルから抽出する。たとえばrx(xは1〜nの任意の値)が信頼閾値以上であれば、残差ベクトルから残差exを抽出する。抽出されなかった残差は信頼性の低いとして以降の処理では用いない。信頼閾値は本発明の第2閾値に対応する。
次に、センサ情報記憶部21から該当する精度情報を読み出し、読み出した精度情報、およびステップS305cで抽出した残差群を利用して、異常スコア(判定スコア)を求める(S305d)。読み出す精度情報は、たとえば診断モデルインスタンス<1>であれば、風量センサAF1,AF2の精度情報と、VAV1、VAV2の精度情報である。
図22は異常スコアの計算方法を説明する図である。
まずステップS305cで抽出した残差群の確率分布Pl（エル）を求める。より詳細に、残差eの平均および分散を計算し、正規分布を仮定して、当該平均および分散に基づいて確率分布Pl（エル）を求める。つまり、Plは、時系列データ記憶部19から読み込んだ時系列データを入力として、l（エル）番目の診断インデックスのパラメータを用いて求められた残差の確率分布である。確率分布Pl（エル）は本発明の第2確率分布に対応する。
また異常がない場合の確率分布Pn,l（エル）を求める。まず基準時系列データ記憶部に記憶された基準時系列データと、上記最適化したl（エル）番目の診断モデルインスタンス<l>とを用いて、ステップS305b,S305cにより残差群を抽出する。そして、抽出した残差群の平均および分散を求め、正規分布を仮定して、当該平均および分散に基づき確率分布Pn,l（エル）を求める。つまり Pn,lは基準時系列データ17から読み込んだ基準時系列データを入力として、l（エル）番目の診断インデックスのパラメータを用いて求められた残差の確率分布である。確率分布Pn,l（エル）は本発明の第1確率分布に対応する。
さらに、センサ・VAV(アクチュエータ)毎に静的に異常が発生したと仮定した場合の確率分布Pa,l（エル）,kを求める。kはアクセステーブルの列番号である。本例では、k=1〜4であるため、確率分布Pa,l（エル）,1、Pa, l（エル）,2、Pa, l（エル）,3、Pa, l（エル）,4を求める。まず基準時系列データ17の基準時系列データにおける、アクセステーブルのl（エル）行k列の値を、該当する精度情報に応じて補正する。補正した基準時系列データと、上記最適化したl（エル）番目の診断モデルインスタンス<l>とを用いて、ステップS305b,S305cにより残差群を抽出する。そして、抽出した残差分の平均および分散を求め、正規分布を仮定して、当該平均および分散に基づき確率分布Pa,l（エル）,kを求める。これをk=1〜4のそれぞれごとに行う。確率分布Pa,l（エル）,kは本発明の第3確率分布に対応する。
k=1〜4のそれぞれに対応するセンサ・アクチュエータの値を補正して残差を求める式を、数式8-1〜数式8-4として以下に示す。
(数式8-1) e(k=1) = (C[i] + α[i, xi +Δxi]) ( Qi / Qj )^2 - (C[j] +α[j, xj])
(数式8-2) e(k=2) = (C[i] + α[i, xi]) ( Qi / Qj )^2 - (C[j] +α[j, xj +Δxj])
(数式8-3) e(k=3) = (C[i] + α[i, xi]) (( Qi + ΔQi )/ Qj )^2 - (C[j] +α[j, xj])
(数式8-4) e(k=4) = (C[i] + α[i, xi]) ( Qi /(Qj + ΔQj))^2 - (C[j] +α[j, xj])
基準時系列データにおいて、アクセステーブルのk列の値を補正する様子を図23に示す。(1)が第1列(k=1)の値を補正、(2)が第2列(k=2)、(3)が第3列(k=3)、(4)が第4列(k=4)の値を補正した結果を示す。たとえば(1)では、数式8-1におけるk=1列の指示値xiをxi +Δxiに補正している。補正の方法としては、ここでは、精度情報のプラス方向(あるいはマイナス方向)の最大誤差を加算する。(1)では最大誤差+0．05を加算している。指示値の補正は、開度の最大誤差に相当する指示値の差分を求め、求めた差分を指示値に補正することで行う。(2)〜(4)でも同様にして補正を行う。精度情報のプラス方向、マイナス方向の両方を考慮する場合は後述する。
以上のように各確率分布Pl、Pn,l（エル）、Pa,l（エル）,1、Pa, l（エル）,2、Pa, l（エル）,3、Pa, l（エル）,4を求めたら、拡張KL情報量(Kulback-Leibler divergence)を用いて、異常スコアの計算を行う。
KL情報量は、確率分布間の差を表す尺度として利用される。KL情報量は対称性がないため距離ではないが、ここで用いる拡張KL情報量は対称性を持つため、確率分布間の距離として定義することができる。
離散確率変数の確率分布P,Qに関する拡張KL情報量は以下の数式9により定義される。拡張KL情報量は，確率分布P,Q が一致する場合0，ずれが大きいほど値が大きくなり、負の値をとらないことが知られている。
Figure 0005337909
このKL情報量に基づき、l（エル）番目の最適化された診断モデルインスタンスにおいて、 k 番目のセンサ・アクチュエータに関する異常スコアScore(l,k) を以下のように求める。
まず、l（エル）番目の最適化された診断モデルインスタンスにおいて、確率分布Pl（エル）と、確率分布Pn,lのずれ(拡張KL情報量)を、下記の数式10により求める。
Figure 0005337909
またl（エル）番目の最適化された診断インデックスモデルにおいて、確率分布P l（エル）と、確率分布Pa,l（エル）,kのずれ(拡張KL情報量)を、下記の数式11により求める。
Figure 0005337909
以上に基づき、k 番目のセンサに関する異常スコアScore(l,k) は，以下のように定義することができる。
Figure 0005337909
図25の上に、診断インデックス<1>〜<6>毎に、k=1〜4のそれぞれに対応するセンサ・アクチュエータに対して計算された異常スコアを示す。
上記説明では、補正の際、精度情報のプラス方向およびマイナス方向の一方のみ考慮したが、プラスとマイナス方向の両方を考慮する場合は、以下のようにすればよい。まず、k=1〜4に対応する数式8-1〜8-4毎に、プラス方向の最大誤差分の補正を行い、計算結果(e)の二乗(e^2)を求める。同様にして、マイナス方向の最大誤差分の補正も行い、計算結果(e)の二乗(e^2)を求める。k=1〜4毎に、プラス方向の補正とマイナス方向の補正で得たe^2群の平均および分散を計算し、カイ二乗分布を仮定して、確率分布Pa,l（エル）,kを求める。一方、基準時系列データに対応する確率分布Pn,l（エル）も、残差eの二乗(e^2)を求め、e^2群の平均および分散を計算し、カイ二乗分布を仮定して、確率分布を求める。また時系列データ記憶部19内の時系列データに対応する確率分布Pl（エル）も、残差eの二乗(e^2)を求め、e^2群の平均および分散を計算し、カイ二乗分布を仮定して、確率分布を求める。以降は先と同様にして、拡張KL情報量を求め、異常スコアを計算する。
次に図18の総合異常判定処理(S306)に進む。
図24は総合異常判定部134による総合異常判定処理の流れを示すフローチャートである。
まずセンサ・アクチュエータ(VAV1〜VAV3、風量センサAF1〜AF3)ごとに、総合スコアを計算する(S306a)。総合スコアは、下記の数式13を用いて計算する。センサ・アクチュエータ毎に、各l（エル）で得られた異常スコアのうち最小値を総合スコアとする。
Figure 0005337909
図25には総合スコアの計算例が示される。VAV1の開度指示値(ID133)の総合スコアは、診断インデックス<1>〜<4>で算出された異常スコア0．96，0.95,0.08,0.09のうちの最小値0．08となる。同様にして他のセンサ・アクチュエータ(ID134,135,143,144,145)についても異常スコアの最小値を選択することにより総合スコアを求める。
次にステップS306bでは、各センサ・アクチュエータの総合スコアを所定の判定閾値と比較し、総合スコアが判定閾値より大きいときは、該当するセンサ・アクチュエータに対して異常ラベルを設定し(S306c)、判定閾値以下のときは正常ラベルを設定する(S306d)。本例では総合スコアが0.92のVAV2に対して異常ラベルが設定され、他のVAV1,3および風量センサAF1〜AF3には正常ラベルが設定される。判定閾値は、本発明の第1閾値に対応する。
ステップS306eでは、異常検出部13の診断結果として、診断対象期間(時系列データを用いた期間)、計算された異常スコア、総合スコア、センサ・アクチュエータに設定したラベルを、結果格納部20へ書き込む (S306e)。
上記説明では診断モデルインスタンス<1>〜<6>毎に、センサ・アクチュエータの異常スコアを求め、センサ・アクチュエータ毎に最小の異常スコアを選択することで総合スコアを求めた。そして、総合スコアが判定閾値より大きいセンサ・アクチュエータに異常ラベルを設定した。別の実施形態として、まず診断モデルインスタンス<1>〜<6>毎に、センサ・アクチュエータの正常スコア(判定スコア)を、下記の数式12aにより求める。そして、センサ・アクチュエータ毎に最大の正常スコアを選択し、総合スコアが判定閾値より小さいセンサ・アクチュエータ（ポイントオブジェクト名）に、異常ラベルを設定し、総合スコアが判定閾値より大きいセンサ・アクチュエータに正常ラベルを設定してもよい。
Figure 0005337909
異常出部13の異常検出提示情報生成部135は、時系列データ記憶部19、基準時系列データ記憶部17、結果格納部20、精度情報記憶部21に記憶されている情報を利用して、予め定められた形式に従って結果表示画面を動的に生成し、生成した画面を表示部16に表示する。
図26、図27，図28は異常検出提示情報生成部136が生成した結果表示の画面例である。
図26は画面内の「リスト」が選択されたときの表示例を示している。総合スコアの大きい順に、各VAVに関するポイントオブジェクト名と、総合スコアが表示される。ただし図示の例では、ポイントオブジェクト名のビル名の表示は省略されている。図左では、異常ラベルが設定されたセンサ・VAV(アクチュエータ)が備え付けられた部屋が網がけ表示されている。
図27は画面内の「配置図」が選択されたとき表示例を示す。
図示の画面は、図26の網がけ表示された部屋のうち12階の部屋Cが選択されたときの表示例を示す。図右では、異常ラベルが設定された、VAV2に関するポイントオブジェクト名「Xビル.12.C.I.VAV.2.開度指示値」が網がけ表示されている。
図28は、画面内の「グラフ」が選択されたときの表示例を示す。
図28では12階の部屋Cが選択されたときの表示例が示される。図右の上では、横軸が時刻、縦軸が開度指示値の座標系が示される。この座標系に、VAV1の開度指示値のグラフG101, VAV2の開度指示値のグラフG201, VAV3の開度指示値のグラフG301が表示されている。
また図右の下では、横軸が時刻、縦軸が風量の座標系が示される。この座標系に、VAV1の風量値(風量センサAF1の値)のグラフG102, VAV2の風量値(風量センサAF2の値)のグラフG202, VAV3の風量値(風量センサAF3の値)のグラフG302が表示されている。
各グラフG101、G201、G301、G102、G202、G302は、診断に用いた日付情報(本例では2009年9月1日0：00〜23：59)を利用して、時系列データ記憶部19からは当該期間の時系列データ(本例では1日分)を読み、読み込んだデータをプロットして得られる。
異常区間Z1はVAV2の開度異常が検出された区間である。詳細を知るため、右下の「解析結果詳細へのリンク」をユーザが入力手段を用いてクリックすると、図29の解析結果画面が生成および表示される。
この画面も、図26〜図28と同様に、異常検出提示情報生成部136が、時系列データ記憶部19、結果格納部20に記憶されている情報を利用して、動的に生成する。
本例は、Xビル.12.C.I.VAV.2.開度指示値に異常ラベルが設定(総合スコア0.92)されたた場合である。この場合、まずXビル.12Cのモデル情報を読み込み、総合スコア0.92を出力した診断インデックス<2>(図25参照)を参照する。診断インデックス<2>のACCESSのインデックス番号133,134,143,144に該当する時系列データ(診断対象期間分)を、時系列データ記憶部19と基準時系列データ記憶部17からそれぞれ読み込む。
画面左上のグラフG101、G201、左下のG102,202は、時系列データ記憶部19から読み込んだデータ(観測値)をプロットしたものであり、図28に示したのと同じである。
画面左上のグラフG211は、VAV2の推定開度指示値である。この推定開度指示値は、数式6を利用して求める。具体的には、数式6のeを0とおき、パラメータを診断モデルインスタンスのPARAMS情報を利用してパラメータを代入する。そして、異常と判断されている x[i] 以外の変数 x[j], Q[i], Q[j] に、時系列データ記憶部19から読み込んだ時系列データを代入してx[i]を求める。x[i] が推定値(期間データ推定値)となる。
画面右側には回帰プロットが示される。この回帰プロットは以下のようにして作成する。上記同様に数式6の変数 x[j], Q[i], Q[j]に、基準時系列データ記憶部17から読み込んだ基準時系列データを代入することにより、x[i] を、比較データの推定値として得る。そして、当該読み込んだ基準時系列データのI.VAV2.開度指示値に関する値と、当該比較データの推定値の交点に「□」をプロットする。一方、時系列データ記憶部19から読み込んだ時系列データのI.VAV2.開度指示値と、上記で求めた期間データ推定値の交点に「×」をプロットする。
このように、正常データ(基準時系列データ)が入力されると、観測値と推定値とが略一致するため横軸値=縦軸値(y=x)に漸近するようにデータがプロットされる。一方、異常データが入力されると、観測値から外れる推定値が出力されるため、y=xから逸脱したところにプロットが行われる。
上記実施形態では、2つのVAVのペア毎に、目的モデルに対し2つの診断モデルインスタンスを生成した。たとえばVAV1とVAV2のペアに対して診断モデルインスタンス＜1＞、＜2＞を生成した。しかしながら、図2のようなシンプルな事例では、ペア毎に、いずれか1つのみを生成しても、異常検出が可能である。したがって、この場合、VAV1とVAV2のペアに対しては、＜1＞および＜2＞のいずれか任意の一方のみを生成し、VAV1とVAV3のペアに対しては、＜3＞および＜4＞のいずれか任意の一方のみを生成し、VAV2とVAV3のペアに対しては、＜5＞および＜6＞のいずれか任意の一方のみを生成してもよい。
また上記実施形態では、目的モデルF01に従って、2つのVAVのすべてのペアを作成し、ペア毎に診断モデルインスタンスを生成したが、本発明はこれに限定されない。たとえば多数のVAVが存在するときに、相互の影響が弱いと事前に分かっているペアについては、診断モデルインスタンスを生成しなくてもよい。診断モデルインスタンスを生成すべきペアは事前に指定しておけばよい。目的モデルF01ではすべてのペアを生成することが暗に示されている。
以上のように、本実施形態によれば、個体差のある未知の特性を持つ機器（センサ・アクチュエータ）の異常検出を高い精度で行うことが可能となる。

(変形例)
これまで説明した実施形態では、ダクトに設置されたVAVの開度指示値と、VAV風量を事例として取り上げたが、本発明は、これに限定されず、たとえば以下のような事例にも適用可能である。
たとえば水量調整装置における冷水弁の開度指示値(vw)と、冷水弁の水量(Qw)の事例も可能である。この場合、数式(1)から、以下のようなモデルインスタンスを生成できる。
min_{Ci,Cj,βi,βj }
R = { (Ci + βi(vwi) ) × (Qwi/Qwj)^2 − (Cj +βj(vwj)) }^2
以上は一例であり、本発明は他の種々の事例にも適用可能である。

Claims (8)

  1. 媒体を流通させる媒体経路における異なる箇所に配置され、それぞれ外部から与えられた指示値に従って前記媒体経路内の前記媒体を共通に出力制御する第1〜第K制御機器と、前記第1〜第K制御機器によって出力制御される媒体の出力量を計測する第1〜第Kセンサに対する異常検出装置であって、
    (Ci + βi(Xi) ) × (Qi/Qj)^2 − (Cj + βj(Xj) (Xiは第1観測変数、Bi(Xi)は Xiを入力変数とする未知の関数である第1未知特性項、Ciは第1定数項、Qiは第2観測変数、Xjは第3観測変数、Bj(Xj)は Xjを入力変数とする未知の関数である第2未知特性項、Cjは第2定数項、Qjは第4観測変数)によって定義される目的モデルを記憶する目的モデル格納部と、
    第1の期間において前記第1〜第K制御機器に与えられた前記指示値と、前記第1〜第Kセンサによる計測値とを含む第1時系列データを所定時間間隔で記憶する第1記憶部と、
    前記第1〜第K制御機器における2つの制御機器のペア毎に、前記目的モデルの前記Xiに一方の制御機器の指示値、前記Ciに前記一方の制御機器の配置箇所に応じた未知の定数、前記Qiに前記一方の制御機器に対応するセンサの計測値、前記Xjに他方の制御機器の指示値、前記Cjに前記他方の制御機器の配置箇所に応じた未知の定数、前記Qjに前記他方の制御機器に対応するセンサの計測値を割り当てることによって、それぞれ診断モデルインスタンスを生成し、
    前記ペア毎に生成したすべての診断モデルインスタンスのそれぞれにおける前記Bi(Xi)およびBj(Xj)を、前記XiおよびXjの値の範囲をそれぞれ分割したセクション毎に一定の値を取る関数と定義し、前記第1記憶部内の前記第1時系列データを用いて、前記すべての診断モデルインスタンスの値の二乗を最小化するように前記Bi(Xi)、前記Ci、前記Bj(Xj)、前記Cjを同定することにより、最適化された診断モデルインスタンスをそれぞれ得るモデル最適化部と、
    前記第1の期間と異なる第2の期間において、前記第1〜第K制御機器に与えられた前記指示値と、前記第1〜第Kセンサによる計測値とを含む第2時系列データを所定時間間隔で記憶する第2記憶部と、
    前記第1〜第Kセンサによる最大観測誤差と、前記指示値に対する前記第1〜第K制御機器の最大応答誤差とを表す精度情報を記憶する第3記憶部と、
    前記最適化された診断モデルインスタンス毎に、
    前記第1記憶部内の前記第1時系列データを用いて前記最適化された診断モデルインスタンスの値である第1残差を求め、前記第1残差の確率分布である第1確率分布を生成し、
    前記第2記憶部内の前記第2時系列データを用いて前記最適化された診断モデルインスタンスの値である第2残差を求め、前記第2残差の確率分布である第2確率分布を生成し、
    前記最適化された診断モデルインスタンスを生成する元となった2つの制御機器および2つのセンサのそれぞれ毎に、前記第1時系列データにおけるそれぞれに対応する値をそれぞれの前記精度情報に応じて変更し、変更後の時系列データを用いて前記最適化された診断モデルインスタンスの値である第3残差を求め、前記第3残差の確率分布である第3確率分布をそれぞれ生成し、
    前記第2確率分布と前記第1確率分布間の距離を表す第1の拡張KL(Kulback-Leibler divergence)情報量と、前記第2確率分布と各前記第3確率分布のそれぞれとの距離を表す第2の拡張KL情報量とを求め、前記第1および第2拡張KL情報量の合計に対する前記第1拡張KL情報量の比率に応じて判定スコアをそれぞれ計算する、
    判定スコア計算部と、
    前記第1〜第K制御機器および前記第1〜Kセンサ毎に、前記最適化された診断モデルインスタンスから得られた判定スコアの最小値である総合スコアを求め、第1閾値より大きい総合スコアをもつ制御機器またはセンサに異常が発生したことを判定する総合異常判定部と、
    を備えた異常検出装置。
  2. 前記モデル最適化部は、前記第1〜第K制御機器毎に、前記セクションに前記指示値が含まれるデータ数をカウントし、データ数が多いほど大きな信頼値を前記セクションに設定し、
    前記判定スコア計算部は、
    前記第1時系列データにおける前記XiおよびXjの値を含むセクションの信頼値を乗算することにより信頼度を計算し、第2閾値以上の前記信頼度が計算されたデータから求められた第1残差のみを用いて前記第1確率分布を生成し、
    前記第2時系列データにおける前記XiおよびXjの値を含むセクションの信頼値を乗算することにより信頼度を計算し、前記第2閾値以上の前記信頼度が計算されたデータから求められた第2残差のみを用いて前記第2確率分布を生成し、
    前記変更後の時系列データにおける前記XiおよびXjの値を含むセクションの信頼値を乗算することにより信頼度を計算し、前記第2閾値以上の前記信頼度が計算されたデータから求められた第3残差のみを用いて前記第3確率分布を生成する、
    ことを特徴とする請求項1に記載の異常検出装置。
  3. 前記モデル最適化部は、さらに、前記Xiに他方の制御機器の指示値、前記Ciに前記他方の制御機器に関する未知の定数、前記Qiに前記他方の制御機器に対応するセンサの計測値、前記Xjに前記一方の制御機器の指示値、前記Cjに前記一方の制御機器に関する未知の定数、前記Qjに前記一方の制御機器に対応するセンサの計測値を割り当てることによって診断モデルインスタンスを生成する
    ことを特徴とする請求項1に記載の異常検出装置。
  4. 前記判定スコア計算部は、
    前記第1残差の二乗を求め、前記第1残差の二乗の確率分布を前記第1確率分布として生成し、
    前記第2残差の二乗を求め、前記第2残差の二乗の確率分布を第2確率分布として生成し、
    前記第3残差の二乗を求め、前記第3残差の二乗の確率分布である第3確率分布を生成する
    ことを特徴とする請求項1に記載の異常検出装置。
  5. 前記判定スコア計算部は、前記第1および第2拡張KL情報量の合計に対する前記第2拡張KL情報量の比率に応じて前記判定スコアを計算し、
    前記総合異常判定部は、前記判定スコアの最大値を前記総合スコアとして求め、前記第1閾値より小さい総合スコアをもつ制御機器またはセンサに異常が発生したことを判定する
    ことを特徴とする請求項1に記載の異常検出装置。
  6. 前記目的モデル格納部は、前記目的モデルとして、前記(Ci + βi(Xi) ) × (Qi/Qj)^2 − (Cj + βj(Xj)に代えて、(Ci + βi(Xi)) × f(Qi,Qj) − (Cj + βj(Xj) ) (f(Qi,Qj)は、前記Qi,Qjを入力変数とする所定の関数)を記憶する
    ことを特徴とする請求項1に記載の異常検出装置。
  7. 前記第1〜第K制御機器は、可変風量装置または固定風量装置であり、前記第1〜第K制御機器に与えられる指示値は、前記可変風量装置または前記固定風量装置の弁の開度指示値であり、前記第1〜Kセンサは風量センサであることを特徴とする請求項1に記載の異常検出装置。
  8. 前記第1〜第K制御機器は、水量調整装置であり、前記第1〜第K制御機器に与えられる指示値は、前記水量調整装置の弁の開度指示値であり、前記第1〜Kセンサは水量センサであることを特徴とする請求項1に記載の異常検出装置。
JP2012507951A 2010-03-30 2010-03-30 異常検出装置 Active JP5337909B2 (ja)

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2010/055695 WO2011121726A1 (ja) 2010-03-30 2010-03-30 異常検出装置

Publications (2)

Publication Number Publication Date
JPWO2011121726A1 JPWO2011121726A1 (ja) 2013-07-04
JP5337909B2 true JP5337909B2 (ja) 2013-11-06

Family

ID=44711517

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2012507951A Active JP5337909B2 (ja) 2010-03-30 2010-03-30 異常検出装置

Country Status (3)

Country Link
US (1) US8577649B2 (ja)
JP (1) JP5337909B2 (ja)
WO (1) WO2011121726A1 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPWO2017134908A1 (ja) * 2016-02-05 2018-06-14 株式会社東芝 センサ故障診断装置、方法、及びプログラム
JP2020035297A (ja) * 2018-08-31 2020-03-05 三菱電機ビルテクノサービス株式会社 機器状態監視装置及びプログラム

Families Citing this family (31)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9197511B2 (en) * 2012-10-12 2015-11-24 Adobe Systems Incorporated Anomaly detection in network-site metrics using predictive modeling
FR3010448B1 (fr) 2013-09-06 2015-08-21 Snecma Procede de surveillance d’une degradation d’un dispositif embarque d’un aeronef avec determination automatique d’un seuil de decision
JP6444885B2 (ja) * 2013-11-26 2018-12-26 株式会社トクヤマ 運動装置の状態監視システム
IL229819A (en) * 2013-12-05 2016-04-21 Deutsche Telekom Ag System and method for anomalously identifying information technology servers through incident consolidation
US10316457B2 (en) * 2014-03-03 2019-06-11 Brent Richard SINGLEY Flood prevention device
JP6356449B2 (ja) * 2014-03-19 2018-07-11 株式会社東芝 センサ診断装置、センサ診断方法、およびコンピュータプログラム
JP2015179443A (ja) * 2014-03-19 2015-10-08 株式会社東芝 診断モデル生成装置、診断用モデル生成方法、及び異常診断装置
US9727533B2 (en) * 2014-05-20 2017-08-08 Facebook, Inc. Detecting anomalies in a time series
US9680646B2 (en) * 2015-02-05 2017-06-13 Apple Inc. Relay service for communication between controllers and accessories
MY192904A (en) 2015-02-17 2022-09-14 Fujitsu Ltd Determination device, determination method, and determination program
US9823160B2 (en) * 2015-04-02 2017-11-21 The Boeing Company Apparatus and methods for testing suction cups mounted to a track
GB201510957D0 (en) * 2015-06-22 2015-08-05 Ge Aviat Systems Group Ltd Systems and Methods For Verification And Anomaly Detection
KR20170007958A (ko) * 2015-07-13 2017-01-23 에스케이하이닉스 주식회사 메모리 시스템
JP6603078B2 (ja) 2015-08-31 2019-11-06 Ｎｔｎ株式会社 車輪用軸受装置
JP6599727B2 (ja) * 2015-10-26 2019-10-30 株式会社Ｓｃｒｅｅｎホールディングス 時系列データ処理方法、時系列データ処理プログラム、および、時系列データ処理装置
US10095757B2 (en) 2015-12-07 2018-10-09 Sap Se Multi-representation storage of time series data
US10685306B2 (en) * 2015-12-07 2020-06-16 Sap Se Advisor generating multi-representations of time series data
US9785495B1 (en) * 2015-12-14 2017-10-10 Amazon Technologies, Inc. Techniques and systems for detecting anomalous operational data
US10642813B1 (en) 2015-12-14 2020-05-05 Amazon Technologies, Inc. Techniques and systems for storage and processing of operational data
US11037066B2 (en) 2016-07-13 2021-06-15 International Business Machines Corporation Estimation of abnormal sensors
JP7017851B2 (ja) * 2016-12-01 2022-02-09 住友重機械工業株式会社 故障診断システムおよび処理ユニット
US10628252B2 (en) * 2017-11-17 2020-04-21 Google Llc Real-time anomaly detection and correlation of time-series data
JP6928669B2 (ja) * 2017-11-28 2021-09-01 株式会社安川電機 制御システム、工場システム、学習システム、推定用モデルの生成方法及びアクチュエータの状態推定方法
US11036715B2 (en) * 2018-01-29 2021-06-15 Microsoft Technology Licensing, Llc Combination of techniques to detect anomalies in multi-dimensional time series
US11792215B1 (en) * 2018-12-04 2023-10-17 Amazon Technologies, Inc. Anomaly detection system for a data monitoring service
EP3796117B1 (de) * 2019-09-18 2021-10-27 Siemens Aktiengesellschaft Diagnoseverfahren und diagnosesystem für eine verfahrenstechnische anlage
CN111459778B (zh) * 2020-03-12 2024-05-07 平安科技（深圳）有限公司 运维***异常指标检测模型优化方法、装置及存储介质
JP7093031B2 (ja) * 2020-09-23 2022-06-29 ダイキン工業株式会社 情報処理装置、情報処理方法、及びプログラム
CN112905671A (zh) * 2021-03-24 2021-06-04 北京必示科技有限公司 时间序列异常处理方法、装置、电子设备及存储介质
US20230259113A1 (en) * 2022-02-11 2023-08-17 Novity, Inc. Subsystem-level model-based diagnosis
US11821324B2 (en) 2022-04-25 2023-11-21 General Electric Company Duct failure detection in a turbine engine

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2003131733A (ja) * 2001-10-26 2003-05-09 Kajima Corp 設備制御監視方法及び設備制御監視装置
JP2003208219A (ja) * 2002-01-10 2003-07-25 Yamatake Building Systems Co Ltd 異常検知装置および方法
JP2008065393A (ja) * 2006-09-04 2008-03-21 Research Organization Of Information & Systems グループ判別装置及びグループ判別方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3254624B2 (ja) 1996-05-31 2002-02-12 株式会社山武 スティックスリップ検出方法および検出装置
US6223544B1 (en) 1999-08-05 2001-05-01 Johnson Controls Technology Co. Integrated control and fault detection of HVAC equipment
US7590513B2 (en) * 2006-01-30 2009-09-15 Nec Laboratories America, Inc. Automated modeling and tracking of transaction flow dynamics for fault detection in complex systems
JP4413915B2 (ja) * 2006-12-13 2010-02-10 株式会社東芝 異常兆候検出装置および方法
JP4782727B2 (ja) * 2007-05-17 2011-09-28 株式会社東芝 機器状態監視装置並びに機器状態監視のための方法およびプログラム
US8515719B2 (en) * 2009-01-14 2013-08-20 Hitachi, Ltd. Apparatus anomaly monitoring method and system

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2003131733A (ja) * 2001-10-26 2003-05-09 Kajima Corp 設備制御監視方法及び設備制御監視装置
JP2003208219A (ja) * 2002-01-10 2003-07-25 Yamatake Building Systems Co Ltd 異常検知装置および方法
JP2008065393A (ja) * 2006-09-04 2008-03-21 Research Organization Of Information & Systems グループ判別装置及びグループ判別方法

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPWO2017134908A1 (ja) * 2016-02-05 2018-06-14 株式会社東芝 センサ故障診断装置、方法、及びプログラム
US10837809B2 (en) 2016-02-05 2020-11-17 Kabushiki Kaisha Toshiba Sensor failure diagnosis device, method and program
JP2020035297A (ja) * 2018-08-31 2020-03-05 三菱電機ビルテクノサービス株式会社 機器状態監視装置及びプログラム
WO2020044898A1 (ja) * 2018-08-31 2020-03-05 三菱電機株式会社 機器状態監視装置及びプログラム
JP7038629B2 (ja) 2018-08-31 2022-03-18 三菱電機ビルテクノサービス株式会社 機器状態監視装置及びプログラム

Also Published As

Publication number Publication date
US20130024172A1 (en) 2013-01-24
US8577649B2 (en) 2013-11-05
JPWO2011121726A1 (ja) 2013-07-04
WO2011121726A1 (ja) 2011-10-06

Similar Documents

Publication Publication Date Title
JP5337909B2 (ja) 異常検出装置
JP6095732B2 (ja) コンピュータベースの方法及びコンピュータベースのモデル
US7917240B2 (en) Univariate method for monitoring and analysis of multivariate data
JP6538057B2 (ja) プロセスユニットの中の全プロセスセグメントの状態を自動的に監視し決定するための、コンピュータで実装される方法およびシステム
EP2791745B1 (en) A method of operating a process or machine
US20210216386A1 (en) Time-sequential data diagnosis device, additional learning method, and recording medium
WO2008085705A1 (en) Method and system for modeling a process in a process plant
US20190188110A1 (en) Industrial control system, and assistance apparatus, control assist method, and program thereof
JP6240050B2 (ja) バイアス推定装置、その方法およびプログラム、並びに故障診断装置、その方法およびプログラム
JP5264796B2 (ja) プラント運転支援装置
CN113287104A (zh) 数据分类装置
CN115129146A (zh) 用于识别关键性能指示符的方法
JP2010146137A (ja) パラメータ調整支援装置
JP5948998B2 (ja) 異常診断装置
JP5811683B2 (ja) 異常診断装置
CN112272804B (zh) 无需动态***模型的工业过程在线故障定位
JP2018044958A (ja) バイアス推定装置、その方法およびプログラム
Pradhan A Dynamic Bayesian Network Framework for Data-Driven Fault Diagnosis and Prognosis of Smart Building Systems
WO2023233858A1 (ja) ニューラルネットワークの計算装置及び計算方法
JP7397623B2 (ja) プラント状態評価装置、プラント状態評価方法およびプログラム
JP2022103931A (ja) モデルを生成する方法、プログラムおよび装置
CN117043700A (zh) 基于模拟的征兆的异常的原因分析
CN115114358A (zh) 信息处理装置、信息处理方法、信息处理***以及计算机程序

Legal Events

Date Code Title Description
TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20130712

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20130805

R151 Written notification of patent or utility model registration

Ref document number: 5337909

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R151