(57)【要約】［課題］個体差のある未知の特性を持つ制御機器等の異常検出を高い精度で行う。［解決手段］本発明は、第1の期間において第1〜第K制御機器に与えられた指示値と、第1〜第Kセンサによる計測値とを含む第1時系列データを記憶する第1記憶部と、前記第1〜第K制御機器における2つの制御機器のペア毎に所定の目的モデルの診断モデルインスタンスを生成し、前記診断モデルインスタンスのパラメータを同定した最適化された診断モデルインスタンスを得るモデル最適化部と、第2の期間に取得された第2時系列データを記憶する第2記憶部と、前記最適化された診断モデルインスタンス毎に前記第1および第2時系列データを用いてそれぞれ関連する制御機器およびセンサに対する判定スコアを計算する判定スコア計算部と、前記第1〜第K制御機器および前記第1〜Kセンサ毎に前記判定スコアに基づき異常有無の判定を行う総合異常判定部とを備える。
2011121726.tif 000002
【特許請求の範囲】
【請求項１】
媒体を流通させる媒体経路における異なる箇所に配置され、それぞれ外部から与えられた指示値に従って前記媒体経路内の前記媒体を共通に出力制御する第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閾値より大きい総合スコアをもつ制御機器またはセンサに異常が発生したことを判定する総合異常判定部と、
を備えた異常検出装置。
【請求項２】
前記モデル最適化部は、前記第1〜第K制御機器毎に、前記セクションに前記指示値が含まれるデータ数をカウントし、データ数が多いほど大きな信頼値を前記セクションに設定し、
前記判定スコア計算部は、
前記第1時系列データにおける前記XiおよびXjの値を含むセクションの信頼値を乗算することにより信頼度を計算し、第2閾値以上の前記信頼度が計算されたデータから求められた第1残差のみを用いて前記第1確率分布を生成し、
前記第2時系列データにおける前記XiおよびXjの値を含むセクションの信頼値を乗算することにより信頼度を計算し、前記第2閾値以上の前記信頼度が計算されたデータから求められた第2残差のみを用いて前記第2確率分布を生成し、
前記変更後の時系列データにおける前記XiおよびXjの値を含むセクションの信頼値を乗算することにより信頼度を計算し、前記第2閾値以上の前記信頼度が計算されたデータから求められた第3残差のみを用いて前記第3確率分布を生成する、
ことを特徴とする請求項1に記載の異常検出装置。
【請求項３】
前記モデル最適化部は、さらに、前記Xiに他方の制御機器の指示値、前記Ciに前記他方の制御機器に関する未知の定数、前記Qiに前記他方の制御機器に対応するセンサの計測値、前記Xjに前記一方の制御機器の指示値、前記Cjに前記一方の制御機器に関する未知の定数、前記Qjに前記一方の制御機器に対応するセンサの計測値を割り当てることによって診断モデルインスタンスを生成する
ことを特徴とする請求項1に記載の異常検出装置。
【請求項４】
前記判定スコア計算部は、
前記第1残差の二乗を求め、前記第1残差の二乗の確率分布を前記第1確率分布として生成し、
前記第2残差の二乗を求め、前記第2残差の二乗の確率分布を第2確率分布として生成し、
前記第3残差の二乗を求め、前記第3残差の二乗の確率分布である第3確率分布を生成する
ことを特徴とする請求項1に記載の異常検出装置。
【請求項５】
前記判定スコア計算部は、前記第1および第2拡張KL情報量の合計に対する前記第2拡張KL情報量の比率に応じて前記判定スコアを計算し、
前記総合異常判定部は、前記判定スコアの最大値を前記総合スコアとして求め、前記第1閾値より小さい総合スコアをもつ制御機器またはセンサに異常が発生したことを判定する
ことを特徴とする請求項1に記載の異常検出装置。
【請求項６】
前記目的モデル格納部は、前記目的モデルとして、前記(Ci + βi(Xi) ) × (Qi/Qj)^2 − (Cj + βj(Xj)に代えて、(Ci + βi(Xi)) × f(Qi,Qj) − (Cj + βj(Xj) ) (f(Qi,Qj)は、前記Qi,Qjを入力変数とする所定の関数)を記憶する
ことを特徴とする請求項1に記載の異常検出装置。
【請求項７】
前記第1〜第K制御機器は、可変風量装置または固定風量装置であり、前記第1〜第K制御機器に与えられる指示値は、前記可変風量装置または前記固定風量装置の弁の開度指示値であり、前記第1〜Kセンサは風量センサであることを特徴とする請求項1に記載の異常検出装置。
【請求項８】
前記第1〜第K制御機器は、水量調整装置であり、前記第1〜第K制御機器に与えられる指示値は、前記水量調整装置の弁の開度指示値であり、前記第1〜Kセンサは水量センサであることを特徴とする請求項1に記載の異常検出装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、制御機器の異常検出を行う異常検出装置に関する。
【背景技術】
【０００２】
制御機器の異常検出の方法の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が配置されたダクト構成を例に、目的モデルの導出過程を簡単に説明する。
【数１】
2011121726.tif 000003
【００４０】
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を式変形することで、上述の一般化された目的モデルを導出できる。
【数２】
2011121726.tif 000004
【００４８】
図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参照)。
【表１】
2011121726.tif 000005
【００６５】
ステップ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，ずれが大きいほど値が大きくなり、負の値をとらないことが知られている。
【数３】
2011121726.tif 000006
【０１６７】
このKL情報量に基づき、l（エル）番目の最適化された診断モデルインスタンスにおいて、 k 番目のセンサ・アクチュエータに関する異常スコアScore(l,k) を以下のように求める。
【０１６８】
まず、l（エル）番目の最適化された診断モデルインスタンスにおいて、確率分布Pl（エル）と、確率分布Pn,lのずれ(拡張KL情報量)を、下記の数式10により求める。
【数４】
2011121726.tif 000007
【０１６９】
またl（エル）番目の最適化された診断インデックスモデルにおいて、確率分布P l（エル）と、確率分布Pa,l（エル）,kのずれ(拡張KL情報量)を、下記の数式11により求める。
【数５】
2011121726.tif 000008
【０１７０】
以上に基づき、k 番目のセンサに関する異常スコアScore(l,k) は，以下のように定義することができる。
【数６】
2011121726.tif 000009
【０１７１】
図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（エル）で得られた異常スコアのうち最小値を総合スコアとする。
【数７】
2011121726.tif 000010
【０１７６】
図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により求める。そして、センサ・アクチュエータ毎に最大の正常スコアを選択し、総合スコアが判定閾値より小さいセンサ・アクチュエータ（ポイントオブジェクト名）に、異常ラベルを設定し、総合スコアが判定閾値より大きいセンサ・アクチュエータに正常ラベルを設定してもよい。
【数８】
2011121726.tif 000011
【０１８０】
異常出部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
【０２０１】
以上は一例であり、本発明は他の種々の事例にも適用可能である。
【図１】
2011121726.tif 000012
【図２】
2011121726.tif 000013
【図３】
2011121726.tif 000014
【図４】
2011121726.tif 000015
【図５】
2011121726.tif 000016
【図６】
2011121726.tif 000017
【図７】
2011121726.tif 000018
【図８】
2011121726.tif 000019
【図９】
2011121726.tif 000020
【図１０】
2011121726.tif 000021
【図１１】
2011121726.tif 000022
【図１２】
2011121726.tif 000023
【図１３】
2011121726.tif 000024
【図１４】
2011121726.tif 000025
【図１５】
2011121726.tif 000026
【図１６】
2011121726.tif 000027
【図１７】
2011121726.tif 000028
【図１８】
2011121726.tif 000029
【図１９】
2011121726.tif 000030
【図２０】
2011121726.tif 000031
【図２１】
2011121726.tif 000032
【図２２】
2011121726.tif 000033
【図２３】
2011121726.tif 000034
【図２４】
2011121726.tif 000035
【図２５】
2011121726.tif 000036
【図２６】
2011121726.tif 000037
【図２７】
2011121726.tif 000038
【図２８】
2011121726.tif 000039
【図２９】
2011121726.tif 000040
【国際調査報告】
2011121726.tif 000041
2011121726.tif 000042