(57)【特許請求の範囲】
【請求項１】
ネットワークを構成する各機器に対して入出力されるエネルギーの計測値と計算値との誤差を最小化する目的関数を用い、状態変数としてのネットワークに対する入出力エネルギーを、遺伝的アルゴリズム、タブサーチ、パーティクルスウォームオプティマイゼーションの何れかのメタヒューリスティク最適化手法により求めて現在の状態変数として出力する最適化部と、
前記最適化部から入力される現在の状態変数であるネットワークに対する入出力エネルギーのもとで、ネットワークを構成する各機器の機器運用ルールを含む入力エネルギーに対する出力エネルギーの関係が関数によってブランチとして表現された機器モデルを多数のノードにより接続し、各ノードにおいて入力エネルギーと出力エネルギーとがバランスしている非線形連立方程式によって表現されたネットワークモデルを有し、各機器の入出力エネルギーを用いて前記非線形連立方程式を解くことにより、定常状態のある時間断面におけるネットワーク内の各部の入出力エネルギーを入力側の機器から末端に向かい逐次計算により求めて、ネットワーク内の各機器における入出力エネルギーの計算値として出力するプログラムからなる定常ネットワークシミュレータと、を備え、
前記最適化部にネットワークの入出力エネルギーの初期値として各機器における現在の入出力エネルギーの計測値を入力し、
前記定常ネットワークシミュレータに現在の状態変数であるネットワークに対する入出力エネルギーを入力して前記計算値を求め、
前記最適化部に前記定常ネットワークシミュレータにより求められた前記計算値を入力し、当該計算値と初期値として入力された前記計測値とを用いて前記メタヒューリスティク最適化手法により新たな状態変数であるネットワークに対する入出力エネルギーを求め、
当該新たな状態変数を現在の状態変数としてネットワークに対する入出力エネルギーを前記定常ネットワークシミュレータに入力して前記計算値を求める処理を繰り返し、
その繰り返し回数が事前の設定回数に達した時点で前記定常ネットワークシミュレータにより求められた現在の状態変数であるネットワークに対する入出力エネルギーに基づく計算値を、最も確からしい計測補正値として出力することを特徴とするネットワークにおける計測値の計測補正値出力装置。
【請求項２】
同時性が確保された計測値群における各計測値に対し、過去の実績を考慮して事前に設定した計測範囲または各計測器の最大・最小計測範囲を逸脱する計測値、及び、明らかに欠損と判る計測値を計測異常値と判定して計測値群より取り除き出力する初期計測異常出力装置と、
前記初期計測異常出力装置から出力された前記計測値群における各計測値の計測補正値を求めて出力する請求項１に記載の計測補正値出力装置と、
前記計測補正値の計測値に対する補正量が、ネットワーク全体のエネルギーフローや各計測器の公差を用いて事前に設定された最大補正量を超過する場合に、前記計測値群に前記初期計測異常出力装置では取り除けない計測異常値が含まれると判定する計測異常判定装置と、
前記計測異常判定装置にて計測異常値が含まれると判定された場合に、前記初期計測異常出力装置から出力される計測値群を対象として、ネットワークを構成する各機器に対して入出力される電気・ガス・水・蒸気のエネルギーの計測値の状態から正常値及び異常値を指定し、正常値として指定された前記エネルギーの計測値のみを対象として当該計測値と計算値との誤差を最小化する目的関数と、前記計測値が正常である確率を考慮した目的関数とを総合した全体目的関数を用い、この全体目的関数による評価値が最も小さい事象において異常値として指定した計測値を計測異常値として出力する計測異常値出力装置と、
を有することを特徴とするネットワークにおける計測値の精査装置。
【請求項３】
請求項２に記載したネットワークにおける計測値の精査装置において、
前記計測異常値出力装置における正常値及び異常値の指定に列挙法を用いることを特徴とするネットワークにおける計測値の精査装置。
【請求項４】
請求項２に記載したネットワークにおける計測値の精査装置において、
前記計測異常値出力装置における正常値及び異常値の指定、並びに、全体目的関数を用いた評価による計測異常値の出力に遺伝的アルゴリズム、タブサーチ、パーティクルスウォームオプティマイゼーションの何れかを用いることを特徴とするネットワークにおける計測値の精査装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、電気・ガス・水・蒸気などが入出力されるネットワークを対象として、これらのエネルギーの計測値(電流・電圧、ガス・水・蒸気の流量・圧力など)のうち最も確からしい計測値を計測補正値として出力し、また、計測センサーの精度や故障に起因して異常と思われる計測値を計測異常値として出力するようにした、ネットワークにおける計測値の精査装置に関するものである。本発明は、電気・ガス・水・蒸気などが入出力されるネットワークを対象として、これらのエネルギーの計測値(電流・電圧、ガス・水・蒸気の流量・圧力など)のうち最も確からしい計測値を計測補正値として出力し、また、計測センサーの精度や故障に起因して異常と思われる計測値を計測異常値として出力するようにした、ネットワークにおける計測値の計測補正値出力装置及び精査装置に関するものである。
【背景技術】
【０００２】
従来、この種のネットワークにおける計測値に関しては、個別の値に対して、過去の実績から値の範囲を事前に設定し、この範囲外のデータに関しては異常値であるという簡単な判断のみ実行されていた。また、この範囲内のデータに関しては、特に他の計測値との関連は考慮することなく正常値として判断されていた。
また、電力系統分野では、非特許文献１に示すように、系統ネットワークを考慮した計測補正値の決定方法として状態推定技術が開発されてきている。この状態推定技術は、ニュートン法あるいは感度を利用した最小自乗法を利用することを基本的な方法としており、ネットワークの方程式も含めた最適化問題の一種として定式化することができ、かつ微分可能であることを前提としている。
【０００３】
【非特許文献１】深尾毅，豊田淳一著，「電力系統へのコンピュータの応用」，産業図書，1972年2月，ｐ151〜ｐ212
【発明の開示】
【発明が解決しようとする課題】
【０００４】
前述の従来技術では、ネットワーク上の実績値としての計測値相互間の関連を考慮していないため、通常の計測値からかけ離れた値以外は、異常値を判定することができない。また、電力系統分野で開発されてきている状態推定技術は、ネットワークの計算式で微分可能であることを前提にしており、ネットワークシミュレータが別ソフトウェアである場合や、各種非線形な特性により微分不可能な場合には適用できないという問題があった。
【０００５】
そこで本発明は、ネットワーク上の計測値相互間の関連を考慮すると共に、ネットワークに関する方程式を利用し、各種の非線形な機器特性や数式化できない運用ルールなどを全て考慮した別プログラムとして構成される定常ネットワークシミュレータを用いて、メタヒューリスティク最適化手法（遺伝的アルゴリズムまたはその改良手法、タブサーチまたはその改良手法、Particle Swarm Optimizationまたはその改良手法）を利用し、過去の計測値の範囲や計測値相互間のつじつまを考慮した最も確からしい計測補正値及び異常計測値を求め、これらを出力可能にしたネットワークにおける計測値の計測補正値出力装置及び精査装置を提供することを目的とする。
【課題を解決するための手段】
【０００６】
上記課題を解決するため、請求項１に記載したネットワークにおける計測値の計測補正値出力装置は、ネットワークを構成する各機器に対して入出力されるエネルギーの計測値と計算値との誤差を最小化する目的関数を用い、状態変数としてのネットワークに対する入出力エネルギーを、遺伝的アルゴリズム、タブサーチ、パーティクルスウォームオプティマイゼーションの何れかのメタヒューリスティク最適化手法により求めて現在の状態変数として出力する最適化部と、
前記最適化部から入力される現在の状態変数であるネットワークに対する入出力エネルギーのもとで、ネットワークを構成する各機器の機器運用ルールを含む入力エネルギーに対する出力エネルギーの関係が関数によってブランチとして表現された機器モデルを多数のノードにより接続し、各ノードにおいて入力エネルギーと出力エネルギーとがバランスしている非線形連立方程式によって表現されたネットワークモデルを有し、各機器の入出力エネルギーを用いて前記非線形連立方程式を解くことにより、定常状態のある時間断面におけるネットワーク内の各部の入出力エネルギーを入力側の機器から末端に向かい逐次計算により求めて、ネットワーク内の各機器における入出力エネルギーの計算値として出力するプログラムからなる定常ネットワークシミュレータと、を備え、
前記最適化部にネットワークの入出力エネルギーの初期値として各機器における現在の入出力エネルギーの計測値を入力し、
前記定常ネットワークシミュレータに現在の状態変数であるネットワークに対する入出力エネルギーを入力して前記計算値を求め、
前記最適化部に前記定常ネットワークシミュレータにより求められた前記計算値を入力し、当該計算値と初期値として入力された前記計測値とを用いて前記メタヒューリスティク最適化手法により新たな状態変数であるネットワークに対する入出力エネルギーを求め、
当該新たな状態変数を現在の状態変数としてネットワークに対する入出力エネルギーを前記定常ネットワークシミュレータに入力して前記計算値を求める処理を繰り返し、
その繰り返し回数が事前の設定回数に達した時点で前記定常ネットワークシミュレータにより求められた現在の状態変数であるネットワークに対する入出力エネルギーに基づく計算値を、最も確からしい計測補正値として出力するものである。
【０００７】
請求項２に記載した精査装置は、同時性が確保された計測値群における各計測値に対し、過去の実績を考慮して事前に設定した計測範囲または各計測器の最大・最小計測範囲を逸脱する計測値、及び、明らかに欠損と判る計測値を計測異常値と判定して計測値群より取り除き出力する初期計測異常出力装置と、
前記初期計測異常出力装置から出力された前記計測値群における各計測値の計測補正値を求めて出力する請求項１に記載の計測補正値出力装置と、
前記計測補正値の計測値に対する補正量が、ネットワーク全体のエネルギーフローや各計測器の公差を用いて事前に設定された最大補正量を超過する場合に、前記計測値群に前記初期計測異常出力装置では取り除けない計測異常値が含まれると判定する計測異常判定装置と、
前記計測異常判定装置にて計測異常値が含まれると判定された場合に、前記初期計測異常出力装置から出力される計測値群を対象として、ネットワークを構成する各機器に対して入出力される電気・ガス・水・蒸気のエネルギーの計測値の状態から正常値及び異常値を指定し、正常値として指定された前記エネルギーの計測値のみを対象として当該計測値と計算値との誤差を最小化する目的関数と、前記計測値が正常である確率を考慮した目的関数とを総合した全体目的関数を用い、この全体目的関数による評価値が最も小さい事象において異常値として指定した計測値を計測異常値として出力する計測異常値出力装置と、
を有するものである。
【０００８】
請求項３に記載した精査装置は、請求項２に記載したネットワークにおける計測値の精査装置において、
前記計測異常値出力装置における正常値及び異常値の指定に列挙法を用いるものである。
【０００９】
請求項４に記載した精査装置は、請求項２に記載したネットワークにおける計測値の精査装置において、
前記計測異常値出力装置における正常値及び異常値の指定、並びに、全体目的関数を用いた評価による計測異常値の出力に遺伝的アルゴリズム、タブサーチ、パーティクルスウォームオプティマイゼーションの何れかを用いるものである。
【発明の効果】
【００１６】
本発明によれば、ネットワークを構成する機器の非線形な特性や数式化できない運用ルール等を考慮し、各計測値間のつじつまを考慮した最も確からしい計測補正値及び計測異常値を出力することができる。
【発明を実施するための最良の形態】
【００１７】
以下、図に沿って本発明の実施形態を説明する。
まず、本発明において使用される定常ネットワークシミュレータについて述べる。
定常ネットワークシミュレータは、ネットワーク及びその構成機器をモデル化し、各機器及びネットワークのエネルギーフローに関する方程式を解くことにより、ある時間断面におけるネットワーク全体の定常状態を計算するものである。
【００１８】
この種のネットワークシミュレータは、例えば、電力エネルギーに関しては、「電力系統解析理論」（関根泰次著，  電気書院，昭和46年1月）などに潮流計算技術として記述された既存技術であり、ガス流及び水道に関しては、「コンピュータ水理学入門」(C.A. ブレビア・A.J.フェラント著，磯部雅彦訳，サイエンス社，1983年)などに管路網解析技術として記述された既存技術である。更に、空調用の水(温水・冷水)の場合には、熱のバランス方程式となる。
エネルギーフロー方程式は、線形あるいは非線形連立方程式として定式化され、ニュートン法などを用いることにより、解、つまり、ある時間断面におけるネットワークの定常状態を得ることができる。
【００１９】
定常ネットワークシミュレータに用いられる機器モデルは、ブランチとして表現され、その入出力関係を数式１の関数により表現する。
【００２０】
【数１】
0004682553.tif 000002
【００２１】
ここで、機器の運用ルールも、上記関数ｆの中で表現される。なお、入力状態量はブランチに対する入力エネルギーから計算可能であり、出力状態量を用いてブランチから出力される出力エネルギーも計算可能である。また、送電線、水道管、ガス管などは一般的には抵抗分を含むため、この抵抗分によりエネルギー損失が生じるブランチも数式１により表現することとする。
【００２２】
次に、定常ネットワークシミュレータに用いられるネットワークモデルについて述べる。ブランチである機器モデルを多数ノードにより接続したものがネットワークモデルとなる。各ノードにおいて、入力エネルギーと出力エネルギーとがバランスしている場合、以下の関数により表現することができる。
【００２３】
【数２】
0004682553.tif 000003
【００２４】
ネットワークの定常シミュレーション結果は、数式２からなる非線形連立方程式を対象として、ネットワークに対する入力エネルギー及び出力エネルギーを指定して方程式を解くことによって得られる。
すなわち、定常ネットワークシミュレータは、機器モデル及びネットワークのエネルギーフローに関する方程式を内部に有し、各機器の運転・停止状態、または入力エネルギー及び出力エネルギーを指定することにより、ある時間断面におけるネットワーク内の各部の状態量を出力する。
【００２５】
このようにして定常ネットワークシミュレータにより計算された各部の入出力状態量は、最も確からしい計測補正値を出力するための目的関数の計算用に、計測値に対する計算値として最適化部に送られることになる。なお、最適化部及び定常ネットワークシミュレータは、何れも計算機のハードウェア及びソフトウェアによって実現されるものである。
例えば、計測値の精査対象となるネットワークとして、ガスタービン，排ガスボイラ，スチームタービン，コンプレッサ，ボイラ，吸収式冷凍機，ガス焚冷凍機，蓄熱層から構成され、空気負荷，電気負荷，蒸気負荷，熱負荷を得るためのコージェネレーションシステムを想定したとする。
【００２６】
この場合の定常ネットワークシミュレータは、図１に示すように、上記各機器の構成・接続状態及び入出力特性が与えられた、ガスタービンモデル２１、排ガスボイラモデル２２，スチームタービンモデル２３，コンプレッサモデル２４，ボイラモデル２５，吸収式冷凍機モデル２６，ガス焚冷凍機モデル２７，蓄熱層モデル２８，空気負荷モデル１３１，電気負荷モデル１３２，蒸気負荷モデル１３３，熱負荷モデル１３４から構成されている。
そして、ある時間断面(例えば制御時間毎)におけるガスタービン（ガスタービンモデル２１），ボイラ（ボイラモデル２５），ガス焚冷凍機（ガス焚冷凍機モデル２７）に対する燃料注入量、及び各機器（各モデル２１〜２８）の運転・停止状態がネットワーク入力状態量として与えられた時に、例えばガスタービン（ガスタービンモデル２１）から出力される排ガスの量やボイラ（ボイラモデル２５）から出力されるスチームの量を計算する。そして、これらを次段の排ガスボイラ（排ガスボイラモデル２２）や吸収式冷凍機（吸収式冷凍機モデル２６）等の入力として、各機器（各モデル２２，２６，……）の出力を逐次計算していくことにより、この定常ネットワークシミュレータで表されるコージェネレーションシステム全体の定常的な入出力状態を計算する。
【００２７】
この際、コージェネレーションシステムの各構成機器の非線形な入出力特性は、以下のような方法を用いて表現する。
（１）最小二乗法により、その係数が調整された２次までの多項式
（２）応答曲面法により、その係数が調整された一般的な非線形方程式
（３）ニューラルネットワーク
なお、（２）の応答曲面法とは、例えば、Andre I. Khuri, and John A. Cornell, Response Surfaces Designs and Analyses, Second Edition, Marcel Dekker Inc., 1996等に説明されている。
また、ニューラルネットワークは、D. E. Rumelhart and J. L. McCelland(Eds.), Parallel Distributed Processing, MIT Press, 1986.等に説明されているところの、多層パーセプトロンを用いたモデルのように、ニューラルネットワークに対する入力と出力との写像関数を自動的に生成する全てのニューラルネットワークを包含する。
【００２８】
ここで、前記の応答曲面法は、タグチメソッドと同様に、製品のばらつきなどを考慮して製品開発を行う品質工学分野において、特にアメリカで実用化されている手法である。
応答曲面法は、適切な実測値の組合せに対する限られた数の応答を求め、予測変数から予測される応答の近似関数を作成する方法であり、この応答曲面法により作成した近似関数を用いて各機器の入出力特性を表現することができる。ここで、応答曲面を求める方法は、ソフトウェア・シミュレーションを用いても、実際に試作したハードウェアによる実験に基づいても、あるいは実測値を用いても良い。
上記近似関数の関数形に制限はないが、近似関数としては以下のような形態が考えられる。
（１）１次関数
数式３の線形関数を利用する場合であり、統計で利用されている重回帰に相当する。
［数３］
ｙ＝β０＋β１ｘ１＋β２ｘ２＋……
ここで、
βi：係数
ｘi：予測変数（設計変数）
である。
（２）２次多項式
数式４の２次多項式を利用する場合である。
［数４］
ｙ＝β０＋β１ｘ１＋β２ｘ２＋……＋β１１ｘ１２＋β２２ｘ２２＋…
…＋β１２ｘ１ｘ２＋β１３ｘ１ｘ３＋……
同様にして、３次以上の多項式を用いることも可能である。また、２次以上の多項式には、以下に示すように直交実験とするために直交多項式を用いる場合もある。
なお、近似関数は、実験点を選定し、その実験点を用いて作成する。
【００２９】
実験点の選定には、以下のような方法がある。
（１）直交計画
タグチメソッドで用いる直交表を用いる場合であり、通常は１次多項式（線形重回帰）による近似に用いられる。しかし、２次以上の多項式回帰モデルでは、L27などの多水準の直交表を用いる。また、この際に、近似関数として直交関数を用いることにより、近似関数の作成において必要な係数の選定を独立に行えるようになるので、有効である。
（２）中央複合計画(Central Composite Design: CCD)
ＣＣＤは、例えば２変数の場合に図２に示す点を実験点とする方法であり、２次多項式による近似で有効であることが知られている。
（３）Ｄ最適水準
計算機支援の実験計画として、以下の数式５で示されるＤ最適性指標を最大化するｐ個の実験点を用いる方法である。
【００３０】
【数５】
0004682553.tif 000004
【００３１】
ここで、Ｘ：実験点の行列(ｎ×ｋ)
ｎ：実験点の組の総数
ｋ： 設計変数の数
ｐ：近似に必要な実験点数
である。
【００３２】
なお、以下に示すように、実験点の選定と用いる近似関数とには密接な関係がある。
・直交計画による実験点の選定と直交関数による近似は計算が簡単であり、更に高次の多項式にも対応できるという利点を有している。しかし、例えば、Ｌ２７直交表＋２次直交関数のＤｅｆｆは0.462であり、余りよい近似とはいえない。
・ＣＣＤは２次多項式だけに有効である。ＣＣＤ＋２次多項式のＤｅｆｆは0.996であり、かなり良い近似となる。
・上記以外の場合には、Ｄ最適規準が必要である。
２次以上の多項式も変数変換することにより１次多項式にすることができるため、１次多項式、つまり統計の分野で線形重回帰式に利用されている方法を用いて近似関数を作成することができる。
【００３３】
（１）モデルの係数の求め方
数式３，４のモデルの各係数βｉは、重回帰で用いられるように最小２乗法で求めることができる。
（２）モデルの適切性の評価
モデルが適切かどうかは、数式６に示される決定係数（回帰モデルの適合性を示す指標）Ｒ２により評価することができる。
数式６において、ＳＳＲは回帰二乗和、Ｓｙｙは応答ｙの平均値まわりの変動である。回帰式が応答に完全に一致する場合にはＳＳＲ，Ｓｙｙが完全に同一になり、決定係数Ｒ２は１となる。この決定係数Ｒ２は０〜１の値をとり、一般には１に近いほど適切な回帰モデルということができる。
なお、モデルの各係数はｔ検定により有効性を確認することができる。
【００３４】
【数６】
0004682553.tif 000005
【００３５】
上記のように、定常ネットワークシミュレータは、例えばコージェネレーションシステムを対象ネットワークとした場合には、各機器の線形または非線形な入出力特性を最小二乗法や応答曲面法、ニューラルネットワーク等により表現し、これらの入出力特性と、時間帯に応じた各機器の運転パターン（始動後に定常状態に達するまでの時間等の特性も含む）等をデータベースとして予め備えている。そして、ある時間断面におけるネットワークへの入力状態量(制御時間毎の燃料注入量や各機器の運転・停止状態)を入力情報として、ネットワーク入力側の機器から末端方向に向かって、ットワーク内の各機器の入出力状態量（入出力されるガスやスチーム、等の量やこれらの入出力時間）を逐次計算することにより、ネットワーク全体の定常的な入出力状態量を計算することとなる。
【００３６】
次に、ネットワークの計測値に基づいて最も確からしい計測補正値を出力するための定式化について述べる。
まず、状態変数としては、ネットワークに対する入力エネルギー及び出力エネルギーとする。一般的には、入力エネルギーは、発電機、ガス製造所、ボイラへの入力、あるいは燃料などとなる。また、出力エネルギーは通常は負荷となる。例えば、電力ネットワークの場合は電力を消費する各種負荷となり、温水・冷水を用いるネットワークの場合は熱を消費する各種負荷となる。
【００３７】
なお、定式化に当たっては、ネットワークの入力エネルギー、ネットワークの出力エネルギー、ネットワークにおける限定されたノード及びブランチの状態量の３種類の値が計測(測定)できると仮定する。
最も確からしい計測補正値を出力するための定式化として、数式７に示すように、目的関数は計測点における計測値と計算値とのノルムの最小化とする。なお、数式８は制約条件を示す。
【００３８】
【数７】
0004682553.tif 000006
【００３９】
【数８】
0004682553.tif 000007
【００４０】
次に、本実施形態における計測補正値の出力方法の基本的な概念について説明する。
電力系統分野において状態推定として確立されてきた技術は、数式７と同じ定式化において、状態変数ｘが連続変数のみで表現でき、かつ、状態方程式ｈｉが微分可能であることを前提としていた。
【００４１】
しかし、本実施形態では、状態方程式ｈｉは定常ネットワークシミュレータにより独立に計算される。このネットワークシミュレータは、上述のように、機器モデル、機器運用ルール、ネットワークモデルなどを含み、ネットワークの入力・出力エネルギーが入力された時に、ネットワークの各部の状態量が計算される。この状態量を用いて、数式７により計測値及び計算値の評価が可能となる。
本発明では、数式８の制約条件を考慮した数式７の目的関数を最小化する最適化問題を、メタヒューリスティク最適化手法（遺伝的アルゴリズムまたはその改良手法、タブサーチまたはその改良手法、Particle Swarm Optimization或いはその改良手法）を利用して解くこととする。
【００４２】
ここで、数式７における重み係数ｗｍｉの決定方法を以下に述べる。
ｗｍｉは、計測値の確からしさに対する各計測値の相対的な重みを表現している。この重みは、以下の２つの項目により評価できる。
1)計測センサーの精度(SensorPRE)
2)計測センサーの故障率(SensorFAU)
計測センサーの精度は、&plusmn;５％などの％誤差範囲で表現される。また、故障率は、年間の故障回数(故障回数／年)で表現される。これらの値を用いて、各計測値に対する重み係数ｗｍｉを数式９のように決定する。
【００４３】
【数９】
0004682553.tif 000008
【００４４】
上記の重み係数ｗｍｉは、例えば数式１０により表現可能である。なお、数式１０におけるａ，ｂは事前に指定するものとする。
【００４５】
【数１０】
0004682553.tif 000009
【００４６】
次に、計測補正値は、以下のアルゴリズムにより求める。
・ステップ１：状態変数(ネットワークに対する入力・出力エネルギー)の初期値として現在の計測値を最適化部に入力する。
・ステップ２：現在の状態変数値を定常ネットワークシミュレータに入力し、計測値に対する計算値を求める。
・ステップ３：最適化部は、ステップ１で入力された計測値とステップ２で計算された各計測値に対する計算値、及び、数式７〜９を用いて、評価値を求める。
・ステップ４：メタヒューリスティク最適化手法（遺伝的アルゴリズムまたはその改良手法、タブサーチまたはその改良手法、あるいは、ＰＳＯまたはその改良手法）を用いて、数式７，８に基づき現在の状態変数値を変更する。
・ステップ５：事前に設定した探索回数に達したらステップ６へ進み、そうでない場合にはステップ２に戻る。
・ステップ６：定常ネットワークシミュレータは、現在の計測値に対する計算値を計測補正値として出力する。
【００４７】
ステップ２〜ステップ５は最適化アルゴリズムの収束過程であり、ステップ４における具体的な方法を以下に述べる。
ここで、図３は本実施形態における計測補正値の計算概念を示す図である。図３において、１０はメタヒューリスティク最適化手法により計測値と計算値との誤差の最小化を目的として現在の状態変数値を変更する（数式７，８を解く）最適化部、２０は最適化部１０からの状態変数（ネットワークの入出力エネルギー）に対するネットワークの定常状態のある時間断面における状態量を計算する定常ネットワークシミュレータを示す。
【００４８】
次に、メタヒューリスティク最適化手法を利用した計測補正値の出力方法について説明する。
１．遺伝的アルゴリズムまたはその改良手法を用いた計測補正値の出力方法
ここで、遺伝的アルゴリズムまたはその改良手法とは、D. E. Goldberg，Genetic Algorithms in Search，Optimization，and Machine Learning，Addison-Wesley，1989の中で述べられているSimple Genetic Algorithm（以下、ＳＧＡという）とその改良手法を指す。以下では、このＳＧＡを例として計測補正値の出力方法を説明する。
【００４９】
(１)状態変数の表現
対象とするネットワークへの入出力エネルギーを番号付けする。対象とするエネルギー数だけの遺伝子座をもつ遺伝子を用意する。ＳＧＡでは、全て離散量としなければならないことから、各遺伝子座に対応する各状態変数(ネットワークへの入力・出力エネルギー)は離散化する最小エネルギー刻みで下限から上限の間の連続量を離散化した値を利用し、離散化した数値の最小値から最大値までを整数と対応させる。内部で整数で表現し、実際のエネルギー量へは対応表より変換する。つまり、全ての状態を整数で表現し、入力・出力エネルギーの数の合計の長さを持つ遺伝子により表現する。
【００５０】
(２)求解アルゴリズム
・ステップ１：前提条件と初期条件設定
1)各計測値を入力する。
2)状態変数(入力・出力エネルギー)の上下限値を設定し、この上下限値から、状態変数を整数値に変換する表を作成する。
3)ストリング数、交差確率、突然変異確率、最大世代数を設定する。
【００５１】
・ステップ２：初期値の生成
1)各ストリングの遺伝子座について、上限及び下限の間の範囲の可能な整数値の中から、ランダムに整数値を生成する。
2)現在の世代数を１とする。
3)評価値の最良値(Obj*)を無限大とする。
【００５２】
・ステップ３：各ストリングの評価と選択
1)ストリングの遺伝子座の整数値と対応表を用いて入出力エネルギー量を決定する。
2)定常ネットワークモデルに入出力エネルギー量を入力し、各計測値に対する計算値を求め、この値を用いて数式３により評価値を計算する。
3)2)の評価値の逆数をFitness関数値として、ルーレット・ホイール・セレクションによりストリングの選択を行う。
4)選択されたストリングの評価値のうち、最小値が現在のObj*より小さかったら、その評価値をObj*とする。また、その時の状態変数の値をX*とする。また、計測値に対する計算値をMeaCalc*とする。
【００５３】
・ステップ４：ストリング操作
ストリング集合に対して、交差確率、突然変異確率を用いて、交差及び突然変異を実行する。
・ステップ５：終了判定
世代が事前に決定した値に達したら終了し、現在のObj*、X*、MeaCalc*を出力する。このうち、MeaCalc*が計測値の補正値となる。
世代が事前に決定した値に達していない場合には、世代数に１足して、ステップ３へ戻る。
【００５４】
２．タブサーチまたはその改良手法を用いた計測補正値の出力方法
ここで、タブサーチまたはその改良手法とはF. Glover，"Tabu Search Part I"，ORSA Journal of Computing，Vol. 1，NO. 3，Summer 1989で述べられているタブーサーチ（以下、ＴＳという）とその改良手法のことを指す。以下、このＴＳを例として計測補正値の出力方法を説明する。
【００５５】
(１)状態変数の表現
対象とするネットワークへの入出力エネルギーを番号付けする。対象とするエネルギー数だけの要素をもつ配列を用意する。ＴＳでは、全て離散量としなければならないことから、各要素に対応する各状態変数(ネットワークへの入力・出力エネルギー)は離散化する最小エネルギー刻みで下限から上限の間の連続量を離散化した値を利用し、離散化した数値の最小値から最大値までを整数と対応させる。内部で整数で表現し、実際のエネルギー量へは対応表より変換する。つまり、全ての状態を整数で表現し、入力・出力エネルギーの数の合計の長さを持つ配列で表現する。
【００５６】
(２)求解アルゴリズム
・ステップ１：前提条件と初期条件設定
1)各計測値を入力する。
2)状態変数(入力・出力エネルギー)の上下限値を設定し、この上下限値から、状態変数を整数値に変換する表を作成する。
3)タブー長及び最大探索回数を設定する。
【００５７】
・ステップ２：初期値の生成
1)状態表現の配列の各要素について、対応する状態変数の上限及び下限の間の範囲の可能な整数値の中からランダムに整数値を生成し、現在状態とする。
2)現在状態をタブーリストに入れる。
3)現在の探索回数を１とする。
4)現在状態に対して、配列の整数を対応表を用いて入出力エネルギー量を決定する。定常ネットワークシミュレータに、現在状態に対する入出力エネルギー量を入力し、各計測値に対する計算値を求め、この値を用いて数式７により評価値を計算する。
5)現在の状態に対する目的関数値を最良値(Obj*)とする。
【００５８】
・ステップ３：隣接状態の生成と次状態の決定
1)現在状態の各配列インデックスに対して整数値の±１の値（上下限値の場合は、限界値を超える方向には値を生成しない）を隣接状態として生成する。例えば、状態変数の数が３の場合、配列の要素は３つとなるが、この配列の値が（２，３，４）であった場合、第１要素の２の±１、第２要素の３の±１、第３要素の４の±１をそれぞれ個別に作成し、隣接状態を作成するため、隣接状態は、（１，３，４），（３，３，４），（２，２，４），（２，４，４），（２，３，３），（２，３，５）の６つとなる。
2)各隣接状態に対して、配列の整数を対応表を用いて入出力エネルギー量を決定する。
3)定常ネットワークシミュレータに、各隣接状態に対する入出力エネルギー量を入力し、各計測値に対する計算値を求め、この値を用いて数式７により評価値を計算する。
4)各隣接状態の中で、タブーでない最も評価がよいものを次状態とし、次状態の評価値が現在のObj*より小さかったら、その評価値をObj*とする。
また、その時の状態変数の値をX*とする。更に、計測値に対する計算値をMeaCalc*とする。
5)現在状態を次状態とする。
【００５９】
・ステップ４：終了判定
探索回数が事前に設定した値に達したら、終了し、現在のObj*，X*，MeaCalc*を出力する。このうち、MeaCalc*が計測補正値となる。
探索回数が事前に設定した値に達しない場合には、現在の探索回数を１足してステップ３へ戻る。
【００６０】
３．Particle Swarm Optimizationまたはその改良手法を用いた計測補正値の出力方法
ここで、Particle Swarm Optimization(以下、ＰＳＯという)またはその改良手法とは、J. Kennedy and R. Eberhart，Swarm Intelligence，Morgan Kaufmann Publishers，2001で解説されているように、Eberhart氏らにより開発された群れ理論を基にしたParticle Swarm OptimizationのGbestモデル、Lbestモデル、あるいは、Angeline氏により開発されたHybrid Particle Swarm OptimizationなどのParticle Swarm Optimizationの様々なバリエーションを含む手法を指す。以下、Gbestモデル及びHybrid PSOを例として、計測補正値の出力方法を説明する。
【００６１】
Ａ．Gbestモデル
（１）始めに、ＰＳＯ研究の背景について述べる。
例えば、魚の群れや鳥の群れなどの自然界の生物の群れの動きは、簡単なルールで記述できることがわかってきている。つまり、個体(エージェント)の動きを簡単なルールで記述しても、群れの複雑な動きを記述することが可能である。
Reynoldsが開発したboid(Reynolds，C. W. ，"Flocks，Herds，and Schools: A Distributed Behavioral Model in Computer Graphics"，Proceedings of SIGGRAPH '87，pp. 25-34，1987.)は、以下の３つの簡単なルールのみを利用している。
1)隣接するエージェントから離れようとする方向へ進む。
2)群れ全体の向かっている方向に進む。
3)群れの中心に向かって進む。
【００６２】
つまり、群れの中の個々のエージェントの動きは、簡単なベクトルの合成で考えることができる。
BoydとRichersonは人間の意思決定のプロセスを調査し、個人の学習と文化の伝達の概念を開発した(R. Boyd and P. Richerson，Culture and the Evolutionary Process，University of Chicago Press，1985.)。彼らの研究によると、人間は意思決定において、２つの重要な情報を利用している。１つ目は、自分自身の経験である。つまり、自分自身でこれまで選択してきたことから、どのような状況が良かったか、その状況がどの程度良かったかについて知っている。２つ目は、他の人々の経験である。つまり、他の人々がどのように行動してどのような結果になったかということを知っている。このような２つの情報、つまり、自分の経験及び他人の経験上、良かった状況に関する情報を用いて人間は意志決定を行っている。
【００６３】
以上のように、群れの個体が簡単なベクトルを用いて動いていること、及び、人間の意志決定が自分の経験及び他人の経験上、良かった状況に関する情報を用いて行われているという２つの観点が、新しい最適化手法であるＰＳＯの背景となっている。
【００６４】
（２）次に、ＰＳＯによる探索方法について述べる。
ＰＳＯは、上述の観点を用いて以下のように開発された。各エージェントの位置(状態量)をXY座標で表現し、その速度をvx（X方向の速度）、vy（Y方向の速度）で表現する。この位置と速度情報から、次時点の各エージェントの位置を更新することができる。この概念に基づき、鳥の群れ全体が何らかの目的関数を最適化するような行動をとると考えると、以下のような最適化が考えられる。
【００６５】
いま、複数のエージェントが共同して良い解を見つける探索を考える。この時、各エージェントは、過去の探索における目的関数の最良値（パーソナルベスト：pbest）とそのXY位置(状態量)を覚えている。このpbestは、BoydとRichersonらの示した個人の経験に対応している。また、各エージェントはpbestのうち集団の中で最も最良なもの、つまり集団の過去の探索における目的関数の最良値（グループベスト：gbest）情報を共有している。
このgbestは、BoydとRichersonらの示した他人の経験に対応している。各エージェントは、現在のXY値とvx，vy値、及びpbestとgbestとの距離に応じてその方向に位置を変更しようとする。この変更しようとする行動は速度で表現される。現在の速度と、pbest及びgbestを用いて、各エージェントの速度は以下の数式１１により修正される。
【００６６】
【数１１】
0004682553.tif 000010
【００６７】
つまり、最適化を実現する複雑な探索の方向は、Boydと同様にベクトルの合成で表現できる。また、数式１１の右辺第２項は、pbestに収束する方向となっており、第３項は、同様にgbestに収束する方向となっている。つまり、この探索方向は、個人の経験と他人の経験を融合した方向となっている。また、数式１１による探索方向の算出により、ＰＳＯはこれまでの速度を維持しようとする大域探索(数式１１の右辺第１項)と、pbest，gbestを用いてそれに近づこうとする局所探索(数式１１の右辺第２，３項)とをバランスよく行う機構をもった探索手法である。
上記の数式１１を用いることにより、各エージェントのこれまでの最良解及び集団の最良解に確率的に近づくような速度が求められ、これにより各エージェントの現在の位置（探索点）を以下の数式１２により修正する。
【００６８】
【数１２】
0004682553.tif 000011
【００６９】
上述した探索点の修正の概念を図４に示す。この図４において、エージェントの現在の探索点ｓｋ(速度vorig)が、pbestに基づく速度ｖpbest及びgbestに基づく速度ｖgbest により、数式７，８に従って新たな探索点ｓｋ＋１に修正される。
ＰＳＯはＧＡ等と同様に複数の探索点を持った多点探索であり、各探索点のpbestと集団のgbestを用いて各探索点を確率的に変更していくことにより、大域最適解（最良解）を得る方法である。また、ＰＳＯは各ステップで目的関数値を評価する必要があるが、評価の回数は問題の規模によらずエージェント数のみで良いというメリットがある。従って、容易に大規模問題への適用が可能である。
図５に、多次元空間におけるＰＳＯによる解探索の概念を示す。この図は、ｎ次元の問題をｍ個のエージェントにより探索する場合を示している。
【００７０】
（３）ＰＳＯのパラメータ
ＰＳＯにおいて重要なパラメータは、試行回数、エージェント数、反復回数、数式１１の各重み係数である。これらは全て設定値として変更可能とする。
1)試行回数
ＰＳＯは確率的な最適化手法であるため、通常１回の最適化のみでなく乱数の種（seed）を変更して複数回の最適化を行い、その結果の中から最も良い解を最適解とする。その最適化を実行する回数を試行回数とする。
【００７１】
2)エージェント数と反復回数
エージェント数は、何点の同時探索を行うかの数値である。反復回数は、各探索点（エージェント）が何回位置を更新して探索を行うかの数値である。つまり、ＰＳＯのアルゴリズムでは、各エージェントの目的関数値は反復回数分評価される。したがって、１回の最適化試行に対して、目的関数値の評価に必要な計算時間×エージェント数×反復回数の計算時間が必要となる。この計算時間を考慮しながら、これらのパラメータを決定する必要がある。
【００７２】
3)重み係数（数式１１におけるw，c1，c2）
重み係数w，c1，c2は、経験的に以下の数式１３，１４のように設定すると良いことが経験的に分かっている。
【００７３】
【数１３】
0004682553.tif 000012
【００７４】
［数１４］  
ｃ１＝ｃ２＝２．０
【００７５】
数式１３によって決定される重みｗは、図６に示すように、探索が進む（反復回数が増える）に従って小さくなることが分かる。これは、数式１１において重みｗのかかる項は上述のように大域探索に相当する項であり、探索の開始時点では大域探索への比重を比較的大きくしておき、探索が進むに従ってその比重を小さくする、つまり局所探索の比重を大きくするような特性となっているためである。
【００７６】
（４）ＰＳＯによる最適化アルゴリズム
一回の最適化試行におけるアルゴリズムの概略を、図７に示す。
すなわち、エージェント数に応じた複数の初期状態をランダムに生成し（Ｓ１）、各エージェントに対する目的関数値を数式７によって計算する（Ｓ２）。
次に、各エージェントの探索過程で目的関数を評価し、最良解が得られた状態をｐｂｅｓｔとして保存する（Ｓ３）。その後、すべてのエージェントのｐｂｅｓｔを比較し、最良解が得られたものをｇｂｅｓｔとして保存する（Ｓ４）。
次いで、最大反復数に達したか否かを判断し、達した場合には現時点のｇｂｅｓｔ（最良解が得られたエージェント）の状態を最適解として出力する（Ｓ５ｙｅｓ，Ｓ７）。
また、最大反復数に達していない場合には、数式１１，１２により制約条件を満たす範囲で各エージェントの速度を計算し、探索点を移動してステップＳ２に戻る（Ｓ５ｎｏ，Ｓ６）。
【００７７】
Ｂ．Gbestモデルを用いた計測補正値の出力方法
（１）状態変数の表現
対象とするネットワークへの入出力エネルギーを番号付けする。対象とするエネルギー数だけの要素をもつ配列を用意する。ＰＳＯでは、状態変数は連続量として表現可能であり、配列の１つ１つの要素が各入出力エネルギーとなっており、入力・出力エネルギーの数の合計の長さを持つ配列で表現する。
【００７８】
（２）求解アルゴリズム
・ステップ１：データの入力
1)各計測値を入力する。
2)状態変数(入力・出力エネルギー)の上下限値を設定する。
3)エージェント数、各最適化パラメータ値、最大探索回数を設定する。
【００７９】
・ステップ２：初期値の生成
1)各エージェント毎に各状態変数値を、設定した上下限値内でランダムに生成する。
2)pbest及びgbestの初期設定
各エージェント毎のpbest値及びgbest値を、事前に設定した大きな値とする。
【００８０】
・ステップ３：各エージェントの評価
1)各エージェント毎の評価値の計算
現在の各状態変数(ネットワークへの入出力エネルギー)を定常ネットワークシミュレータに入力し、各計測値に対する計算値を求める。この計算値を用いて数式７により評価値を計算する。
2)pbest及びgbestの更新
1)で計算した各エージェント毎の評価値が現在の各エージェント毎のpbest値より良かったら、現在の値をpbestに変更する。
上記で計算したpbestのうちの最良値が現在のgbestより良かったら、その値をgbest値に変更する。
【００８１】
・ステップ４：各エージェントの探索点の修正
設定した上下限値内で、各エージェント毎に状態変数値を数式１１，１２を用いて修正する。
【００８２】
・ステップ５：終了条件のチェック
探索回数が、入力した最大探索回数に達したら終了し、現状のgbestの評価値、状態変数の値、及び、各計測値に対する計算値（計測補正値）を出力する。
探索回数が最大探索回数に達しない場合は、ステップ３へ戻る。
【００８３】
Ｃ．Hybrid  PSO(ＨＰＳＯ)
ＨＰＳＯは、ＰＳＯのメカニズムとＧＡ等で利用される自然淘汰の概念を組み合わせたハイブリッド手法である(P. Angeline，"Using Selection to Improve Particle Swarm Optimization"，Proc. of IEEE International Conference on Evolutionary Computation (ICEC)，Anchorage，May 1998.)。
メタヒューリスティック手法（ＭＨ手法）の多くは、自然淘汰に対応する選択の概念を入れている。選択は、探索空間の中で、他の探索点が最近探索した点と比較して相対的に有効である探索点に探索方向を向け直す機能となる。ＰＳＯもそれまでの探索で有効であった探索点(pbest，gbest)を利用しているが、その点自体にかなり依存してそれ以降の探索点が限定されてしまう。
【００８４】
これに対し、ＨＰＳＯは、通常のＭＨ手法と同様に有効な点は利用するが、その点への依存が徐々に薄らいでいくような機能を実現している。具体的には、現状の探索点の評価値を良い順番に並べ、事前に設定した割合分だけ、現状の探索点に対し評価値の最も悪いエージェントの探索点と速度から順番に、評価値の最も良いエージェントの値にリプレイスする。
この際、各エージェントのこれまで探索してきた最も評価が良い探索点(pbest)の情報は残す。このような方法により、有効な領域への探索の集中とこれまで探索してきた有効な領域への探索の方向の向け直し及びその方向への弱い依存関係という機能を実現することができる。
【００８５】
ＨＰＳＯの一般的なアルゴリズムは、図８のようになる。
すなわち、エージェント数に応じた複数の初期状態をランダムに生成し（Ｓ１）、各エージェントに対する目的関数値を数式７によって計算する（Ｓ２）。
次に、各エージェントの探索過程で目的関数を評価し、最良解が得られた状態をｐｂｅｓｔとして保存する（Ｓ３）。その後、すべてのエージェントのｐｂｅｓｔを比較し、最良解が得られたものをｇｂｅｓｔとして保存する（Ｓ４）。
更に、評価の悪いエージェントの探索点及び速度を、評価の良いエージェントの探索点及び速度に変更することにより、ＨＰＳＯの主要部である選択を実行する（Ｓ４Ａ）。
次いで、最大反復数に達したか否かを判断し、達した場合には現時点のｇｂｅｓｔ（最良解が得られたエージェント）の状態を最適解として出力する（Ｓ５ｙｅｓ，Ｓ７）。
また、最大反復数に達していない場合には、数式１１，１２により制約条件を満たす範囲で各エージェントの速度を計算し、探索点を移動してステップＳ２に戻る（Ｓ５ｎｏ，Ｓ６）。
【００８６】
このアルゴリズムにより、ＨＰＳＯはＰＳＯとＭＨ手法の選択のダイナミクスの混合となる。ＨＰＳＯと従来のＰＳＯのメカニズムとしての差は少ないが、この選択手法を加えることによって、評価の高い領域を集中的に探索するメカニズムが追加され、より良い解を見つけることができる探索メカニズムとなる。
【００８７】
Ｄ．ＨＰＳＯを用いた計測補正値の出力方法
計測データ及びデータの仮定、定式化、状態変数の取り方は、ＰＳＯを用いた方法と同様である。ＨＰＳＯを用いた求解アルゴリズムは以下のようになる。
・ステップ１：データの入力
1)各計測値を入力する。
2)状態変数(入力・出力エネルギー)の上下限値を設定する。
3)エージェント数、各最適化パラメータ値、最大探索回数を設定する。
【００８８】
・ステップ２：初期値の生成
1)各エージェント毎に各状態変数値を、設定した上下限値内でランダムに生成する。
2)pbest及びgbestの初期設定
各エージェント毎のpbest値及びgbest値を、事前に設定した大きな値とする。
【００８９】
・ステップ３：各エージェントの評価
1)各エージェント毎の評価値の計算
現在の各状態変数(ネットワークへの入出力エネルギー)を定常ネットワークシミュレータに入力し、各計測値に対する計算値を求める。この計算値を用いて数式７により評価値を計算する。
2)pbest及びgbestの更新
1)で計算した各エージェント毎の評価値が現在の各エージェント毎のpbest値より良かったら、現在の値をpbestに変更する。
上記で計算したpbestのうちの最良値が現在のgbestより良かったら、その値をgbest値に変更する。
【００９０】
・ステップ４：選択の実行
1)現在の各エージェントの評価値を、良い順番にソーティングする。
2)事前に設定した割合(Sr)を用いて、変更する個体数(Ns)を以下の数式により計算する。
Ns = Round(AgentNum&times;Sr)
3)最も評価の悪いエージェントからNs個だけの探索点及び速度を、最も評価の良いエージェントからNs個だけの探索点及び速度と変更する。この際、変更したエージェントのpbest情報はそのままとする。
【００９１】
・ステップ５：各エージェントの探索点の修正
設定した上下限値内で、各エージェント毎に状態変数値を数式１１，１２を用いて修正する。
【００９２】
・ステップ６：終了条件のチェック
探索回数が、入力した最大探索回数に達したら終了し、現状のgbestの評価値、状態変数の値、及び、各計測値に対する計算値（計測補正値）を出力する。
探索回数が最大探索回数に達しない場合は、ステップ３へ戻る。
【００９３】
４．計測異常値の出力方法の定式化
Ａ．計測異常値を含めた計測値の事象確率の計算方法
対象とする各計測値において、正常値が計測できる確率Ｐｃｉを事前に設定する。この確率を用いて、いま、NumNorMea個の計測値が正常であった場合の事象確率Ｐｋは、以下のように計算することができる。
【００９４】
【数１５】
0004682553.tif 000013
【００９５】
ここで、Ｐｃｉは例えば以下のように設定することができる。
Ｐｃｉは、計測値が正常に計測できる確率であるが、以下の定式化で示すように、各状態量の相対的な比較に用いられるものであり、確率値そのものを絶対値として用いることはない。従って、相対的に大きいか小さいかが判断できればよいことになり、以下のように設定可能である。つまり、この値は以下の２つの項目により評価できる。
1)計測センサーの精度(SensorPRE)
2)計測センサーの平均故障時間(SensorFT)
【００９６】
計測センサーの精度は、±５％などの％誤差範囲で表現される。また、平均故障時間は、１回の故障時の平均故障持続時間(時間／年)で表現される。この値を用いて、各計測値に対するＰｃｉは、以下の数式１６のように計算する。
【００９７】
【数１６】
0004682553.tif 000014
【００９８】
例えば、Ｐｃｉは以下の数式１７により表現可能である。
なお、数式１７における重みｗｃ１，ｗｃ２は事前に指定する。
【００９９】
【数１７】
0004682553.tif 000015
【０１００】
Ｂ．計測異常値の最適化問題としての定式化
(１)状態変数
計測値数だけの状態変数とし、各状態変数は１，０で要素を表現する。１は計測正常値、０は計測異常値とする。このように１，０のバイナリ値で状態変数は表現される。
【０１０１】
(２)目的関数
前記事象確率Ｐｋを用いた場合、計測異常値に対する目的関数は以下のように表される。
【０１０２】
【数１８】
0004682553.tif 000016
【０１０３】
また、計測値と計算値との誤差に関する目的関数（数式７）を利用することにより、全体の評価値は、以下のように表される。
【０１０４】
【数１９】
0004682553.tif 000017
【０１０５】
つまり、目的関数は数式２０のように表現でき、一例を挙げると数式２１のようになる。
【０１０６】
【数２０】
0004682553.tif 000018
【０１０７】
【数２１】
0004682553.tif 000019
【０１０８】
５．計測異常値の出力方法
計測異常値の最適化問題としての定式化を用いた異常値の決定方法を、以下に述べる。
計測値数だけの要素数を持つ配列を用意し、各要素は１または０の数値を取ることとする。各要素は、番号付けされた対象とする各計測値が正常値である場合に１とし、異常値である場合に０とする。
【０１０９】
(１)正常値及び異常値が指定された場合の評価値の計算方法
正常値及び異常値が指定された場合、数式７について、誤差最小化の対象とする計測値を正常と指定された計測値のみとする。すなわち、数式２２により評価する。
【０１１０】
【数２２】
0004682553.tif 000020
【０１１１】
数式２２及び数式１８を用いて、数式２０により評価値を計算することができる。
【０１１２】
(２)評価値を用いた計測異常値の決定方法
状態変数のうち、その状態変数に対する数式２０の評価値が最も小さくなるものの状態変数の異常値と指定したものを、計測異常値として決定する。
【０１１３】
計測異常値の出力方法の全体アルゴリズムは、以下の通りである。
・ステップ１：異常値及び正常値の指定
列挙法、遺伝的アルゴリズムまたはその改良法、タブサーチまたはその改良法、ＰＳＯまたはその改良手法などの方法の各探索過程において、異常値及び正常値の指定(ある事象の指定)を行う。
・ステップ２：指定された事象に対する評価値の計算
数式１６を用いて、ステップ１で指定された事象に対する評価値を計算し、最も評価値が小さい事象を格納する。
・ステップ３：終了条件を満たせば終了する。終了条件を満たしていない場合には、ステップ１へ戻る。
・ステップ４：計測異常値の決定
最も評価値が小さい事象において、異常値と指定した計測値を計測異常値として出力する。
【０１１４】
６．列挙法を用いた計測異常値の出力方法
上述した計測異常値の出力方法を、列挙法を用いて実現する。
(１)状態変数
各計測値を番号付けし、計測値数だけの要素数を持つ配列を用意する。配列の各要素は対応する番号付けされた計測値の状態を表現し、１，０で要素を表現する。１は計測正常値、０は計測異常値とする。このように１，０のバイナリ値で配列値を表現する。
【０１１５】
(２)求解アルゴリズム
・ステップ１：事象の指定とその評価値の計算
配列の各要素に対し、１或いは０の要素を生成し、要素の組み合わせからなる配列の値の２NumMea個の全ての組合せを生成し、各々の組合せに対し、 数式２０を用いて指定された事象に対する評価値を計算する。
・ステップ２：最も評価値が小さい事象の確定
ステップ１で計算した全ての事象の評価値のうち、最も評価値が小さい事象を格納する。
・ステップ３：計測異常値の決定
最も評価値が小さい事象において、異常値と指定した計測値を計測異常値として出力する。
【０１１６】
７．遺伝的アルゴリズムまたはその改良手法を用いた計測異常値の出力方法
前述した計測異常値の出力方法を、遺伝的アルゴリズムまたはその改良手法を用いて実現する。
(１)状態変数
各計測値を番号付けし、計測値数だけの遺伝子座を持つ遺伝子を用意する。遺伝子の各遺伝子座は対応する番号付けされた計測値の状態を表現し、１，０のバイナリ値で表現する。１は計測正常値、０は計測異常値とする。このように１，０のバイナリ値で遺伝子を表現する。
【０１１７】
(２)求解アルゴリズム
・ステップ１：前提条件と初期条件設定
1)各計測値を入力する。
2)ストリング数、交差確率、突然変異確率、最大世代数を設定する。
【０１１８】
・ステップ２：初期値の生成
1)各ストリングの遺伝子座について、ランダムにバイナリ値(１，０)を生成する。
2)現在の世代数を１とする。
3)評価値の最良値(Obj*)を無限大とする。
【０１１９】
・ステップ３：各ストリングの評価と選択
1)各ストリングの遺伝子座の１，０値により正常・異常データを指定し、その指定のもとで数式２２により目的関数を計算することとし、請求項１〜４に示した計測値の補正方法のいずれかを用いて、最も評価値が小さくなる補正値を求める。
2)数式２０により目的関数値を計算する。
3)2)の評価値の逆数をFitness関数値として、ルーレット・ホイール・セレクションによりストリングの選択を行う。
4)選択されたストリングの評価値のうち、最小値が現在のObj*より小さかったら、その評価値をObj*とする。また、その時の状態変数の値をX*とする。
【０１２０】
・ステップ４：ストリング操作
ストリング集合に対して、交差確率、突然変異確率を用いて、交差及び突然変異を実行する。
・ステップ５：ＧＡによる探索の終了判定
世代が事前に決定した値に達したら、ステップ６へ移行し、達しない場合には世代数を１足して、ステップ３へ戻る。
・ステップ６：計測異常値の決定
現在のX*において、異常値と指定された計測値を計測異常値として出力し、終了する。
【０１２１】
８．タブサーチまたはその改良手法を用いた計測異常値の出力方法
前述した計測異常値の出力方法を、タブサーチまたはその改良手法を用いて実現する。
(１)状態変数
各計測値を番号付けし、計測値数だけの要素数を持つ配列を用意する。配列の各要素は対応する番号付けされた計測値の状態を表現し、１，０のバイナリ値で表現する。１は計測正常値、０は計測異常値とする。このように１，０のバイナリ値で配列を表現する。
【０１２２】
(２)求解アルゴリズム
・ステップ１：前提条件と初期条件設定
1)各計測値を入力する。
2)タブー長及び最大探索回数を設定する。
【０１２３】
・ステップ２：初期値の生成
1)状態表現の配列の各要素に対して、ランダムにバイナリ値(１，０)を生成する。
2)現在状態をタブーリストに入れる。
3)現在の探索回数を１とする。
4)現在状態の配列要素の１，０値により正常・異常データを指定し、その指定のもとで、数式２２により目的関数を計算することとし、請求項１〜４に示した計測値の補正方法のいずれかを用いて最も評価値が小さくなる補正値を求める。
5)数式２０により目的関数値を計算する。
6)現在の状態に対する目的関数値を最良値(Obj*)とする。
【０１２４】
・ステップ３：隣接状態の生成と次状態の決定
1)現在状態の各配列要素に対して現在値と異なるバイナリ値を生成し、他の要素数は現在値とする配列値を順番に生成し、これを隣接状態とする。例えば、計測値数が３の場合、配列の要素は３つとなるが、この配列の値が（１，１，０）であった場合、第１要素のみ現在と異なるバイナリ値、第２要素のみ現在値と異なるバイナリ値、第３要素のみ現在値と異なるバイナリ値をそれぞれ個別に作成して隣接状態を作成するため、隣接状態は、（０，１，０），（１，０，０），（１，１，１）の３つとなる。
この方法によると、隣接状態は、必ず対象とする計測値数だけ生成される。
2)隣接状態の各状態の配列要素の１，０値により正常・異常データを指定し、その指定のもとで、数式２２により目的関数を計算することとし、請求項１〜４に示した計測値の補正方法のいずれかを用いて、最も評価値が小さくなる補正値を求める。
3)数式２０により、各隣接状態に対する目的関数値を計算する。
4)各隣接状態の中で、タブーでない最も評価がよいものを次状態とする。
次状態の評価値が現在のObj*より小さかったら、その評価値をObj*とし、その時の状態変数の値をX*とする。Obj*とX*の履歴を保存する。
5)現在状態を次状態とする。
【０１２５】
・ステップ４：タブサーチによる探索の終了判定
探索回数が事前に設定した値に達したら、ステップ５へ移行し、達しない場合には現在の探索回数に１足してステップ３へ戻る。
【０１２６】
・ステップ５：異常データの決定
現在のX*及び評価が悪くなる方向のObj*の履歴のうち、事前に決められた数の組合せ分だけ、異常値と指定された計測値を計測異常値として出力して終了する。
【０１２７】
９．Particle Swarm Optimizationまたはその改良手法を用いた計測異常値の出力方法
前述した計測異常値の出力方法を、Particle Swarm Optimizationまたはその改良手法を用いて実現する。
Ａ．バイナリＰＳＯ
状態変数が１，０で表現された場合のバイナリＰＳＯは、エバーハート等によって開発された(J. Kennedy and R. Eberhart，"A discrete binary version of the particle swarm optimization algorithm"，Proc. of the 1997 conference on Systems，Man and Cybernetics (SMC'97)，pp.4104-4109，1997.)。
【０１２８】
GbestモデルやＨＰＳＯは、連続変数を用いた非線形最適化問題に対して開発された。
しかし、現実の多くの問題は組合せ最適化問題として定式化できる。エバーハート等はエージェントのイエスかノーか、正しいか正しくないか、あるいは他の意思決定を行う確率Ｐが個人と社会要因との関数で表現される数式２３のようなモデルを提案した。
【０１２９】
【数２３】
0004682553.tif 000021
【０１３０】
つまり、オリジナルなＰＳＯで利用されているもの(状態、速度、pbest，gbest)で、次の状態が１になる確率を関数表現できるとした。また、エージェントにおける１つあるいは他の決定の傾向であるパラメータｖは、確率値のしきい値を決定する。もし、ｖが大きかったらエージェントはより１を選択し、小さかったら０を選択する。そのようなしきい値は、[0，1]の範囲で決定される。このような特性を実現するための１つの関数は、数式２４に示すシグモイド関数である。このシグモイド関数は、ニューラルネットワークの出力関数の非線形な特性を表現するためにも利用されている。
【０１３１】
【数２４】
0004682553.tif 000022
【０１３２】
あるエージェントの傾向は、そのエージェントのそのグループの成功のために修正されなければならない。このため、それぞれの速度の式は、エージェントの現在状態と、そのエージェント及びグループのそれまで最も評価が高かった状態との差の関数となる。つまり、オリジナルの連続量に対するＰＳＯのように、バイナリＰＳＯは、以下の数式２５，２６のように記述できる。
【０１３３】
【数２５】
0004682553.tif 000023
【０１３４】
【数２６】
0004682553.tif 000024
【０１３５】
結果的には、前述した数式１１，１２の代わりに数式２５，２６を用いることが異なるだけで、後はオリジナルのＰＳＯと同様のアルゴリズムを利用することができる。速度ｖとしては、sig(vik)が０または１に急激に近くならないように最大値を決定する必要がある。
つまり、速度は常に１から０、あるいは０から１に変更する可能性があるようにする必要がある。この速度の最大値は、事前に決定する必要があるが、例えば[-4.0，+4.0]に設定される。
【０１３６】
Ｂ．バイナリＰＳＯを用いた計測異常値の出力方法
(１)状態変数
各計測値を番号付けし、計測値数だけの要素数を持つ配列でエージェントの状態を表す。配列の各要素は対応する番号付けされた計測値の状態を表現し、１，０のバイナリ値で表現する。１は計測正常値、０は計測異常値とする。このように１，０のバイナリ値で配列を表現する。
【０１３７】
(２)求解アルゴリズム
・ステップ１：前提条件と初期条件設定
1)各計測値を入力する。
2)エージェント数、各最適化パラメータ値、最大探索回数を設定する。
【０１３８】
・ステップ２：初期値の生成
1)各エージェントの状態表現の配列の各要素に対して、ランダムにバイナリ値(１，０)を生成する。
2)pbest及びgbestの初期設定
各エージェント毎のpbest値及びgbest値を事前に設定した大きな値とする。
【０１３９】
・ステップ３：各エージェントの評価
1)各エージェント毎の評価値の計算
現在状態の配列要素の１，０値により正常・異常データを指定し、その指定のもとで、数式２２により目的関数を計算することとし、請求項１〜４に示した計測値の補正方法のいずれかを用いて、最も評価値が小さくなる補正値を求める。
その後、数式２０により目的関数値を計算する。
2)pbest及びgbestの更新
1)で計算した各エージェント毎の評価値が現在の各エージェント毎のpbest値より良かったら、現在の値をpbestに変更する。
上記で計算したpbestのうちの最良値が現在のgbestより良かったら、その値をgbest値に変更する。gbestの履歴は保存しておく。
【０１４０】
・ステップ４：各エージェントの探索点の修正
各エージェント毎に、状態変数値を数式２５，２６を用いて修正する。
・ステップ５：ＰＳＯによる探索の終了条件のチェック
探索回数が入力した最大探索回数に達したら終了し、ステップ６へ移行する。達しない場合には、ステップ３へ戻る。
・ステップ６：異常データの決定
現在のgbestから順番に評価値が低くなる方向のgbestの履歴のうち、事前に決められた数の組合せだけ、異常値と指定された計測値の組合せを計測異常値として出力し、終了する。
【０１４１】
以下に、請求項２〜４に記載したネットワークにおける計測値の精査装置の実施形態について説明する。図９は、この実施形態に係る精査装置の構成図である。
この精査装置１００は、同時性を確保した（同期のとれた）計測値群がデータ入力装置６０を介して入力されると共に、計測補正値出力装置３０により求められた計測補正値と、計測異常値出力装置４０や初期計測異常出力装置７０により求められた計測異常値とを組み合わせて、データ出力装置５０を介して出力するものである。
【０１４２】
図９において、精査装置１００には、ネットワーク中に設置された各計測値から同期のとれた計測値が入力される。この同期のとり方には、例えばGPSなどを利用し、事前に決定した周期（例えば１０分単位）で計測値と同時に時間タグを付けた計測データを装置に入力し、装置側でもGPSに合わせた計測時間を持ち、各計測機器の計測データのうち、この装置側の計測時間に一番近いデータをその時点での計測値とする方法が考えられる。
また、この方法はGPSを用いるのではなく、各計測器をネットワークで接続し、ＮＴＰサーバなどにより、絶対時間の同期をとり、この絶対的な計測時間の時間タグと計測値とを装置に入力し、同じ計測時間の計測値を同期のとれた計測値としてみる方法も考えられる。
【０１４３】
ここで、データ入力装置６０は、キーボード，ＬＡＮ，インターネット等の通信設備で構成されており、計測補正値出力装置３０による計測補正値の計算、計測異常値出力装置４０による計測異常値の計算、および初期計測異常出力装置７０や計測異常判定装置８０による計測異常値判定に必要な各種データやパラメータを入力する装置である。
初期計測異常出力装置７０は、同時性が確保された計測値群における各計測値に対し、過去の実績を考慮して事前に設定した計測範囲または各計測器の最大・最小計測範囲を逸脱する計測値および明らかに欠損と判る計測値を計測異常値と判定して出力し、その計測異常値を計測値群から取り除く装置である。
計測補正値出力装置３０は、前述した最適化部１０および定常ネットワークシミュレータ２０により構成され、初期計測異常出力装置７０により前記初期計測異常値が取り除かれた計測値群を用いて、請求項１によって最適な計測補正値を求めて出力する装置である。
計測異常判定装置８０は、計測補正値出力装置３０により求められた計測補正値の計測値に対する補正量が、ネットワーク全体のエネルギーフローや各計測器の公差等を用いて事前に設定された最大補正量を超過する場合に前記計測値群に初期計測異常出力装置７０では取り除けない計測異常値が含まれると判定し、計測異常値の出力装置４０に実行指令を送る装置である。
計測異常値の出力装置４０は、計測異常判定装置８０により計測異常値が含まれると判断された計測値群に対し、請求項２〜４の何れかにより計測異常値を求めて出力する装置である。
データ保存・蓄積装置９０は、ハードディスク，ＤＶＤ，ＣＤ-Ｒ，ＭＯ，フレキシブルディスク，メモリ等で構成され、データ入力装置６０から入力されたデータやパラメータ、計測補正値出力装置３０から得られた計測補正値、および計測異常値の出力装置４０や初期計測異常出力装置７０から得られた計測異常値等を保存・蓄積する装置である。
データ出力装置５０は、ディスプレイ，伝搬装置，プリンタ等により構成されており、データ入力装置６０から入力されたデータやパラメータ、計測補正値の出力装置３０から得られた計測補正値、および計測異常値出力装置４０や初期計測異常出力装置７０から得られた計測異常値等を表示，伝送，印刷する装置である。
【実施例１】
【０１４４】
次に、本発明の実施例１について述べる。
１．計算条件
図１０に、本実施例が適用されるコジェネレーションシステムを示す。図において、６１はガスタービン、６２はスチームタービン、６３は熱負荷、６４は電気負荷であり、データの計測点としては５つの計測点ｍ１〜ｍ５を仮定している。
【０１４５】
(１)機器モデル
1)ガスタービン６１
ここでは、排熱ボイラも含めてガスタービン６１としてモデル化しており、ガスタービン６１は、燃料を入力することにより、電気と蒸気を出力するモデルとなっている。ガスタービン６１の入出力関係は数式２７，２８に示す通りであり、入出力量は数式２９，３０に示すように各々線形形式で表現されているとする。
［数２７］
ｙ１＝ｆ１（ｘ１）
［数２８］
ｙ２＝ｆ２（ｘ１）
［数２９］
ｙ１＝ａ１ｘ１＋ｂ１
［数３０］
ｙ２＝ａ２ｘ１＋ｂ２
ここで、ｘ１：ガスタービンへの燃料入力量
ｙ１：ガスタービンの電気出力量
ｙ２：ガスタービンの蒸気出力量
ａ１，ｂ１，ａ２，ｂ２：機器特性パラメータ
【０１４６】
2)スチームタービン６２
スチームタービンは、蒸気を入力することにより、電気と熱(蒸気または温水)を出力するモデルとなっている。スチームタービン６２の入出力関係は数式３１，３２に示す通りであり、入出力量は数式３３，３４に示すように各々線形形式で表現されているとする。
［数３１］
ｙ３＝ｆ３（ｘ２）
［数３２］
ｙ４＝ｆ４（ｘ２）
［数３３］
ｙ３＝ａ３ｘ２＋ｂ３
［数３４］
ｙ４＝ａ４ｘ２＋ｂ４
ここで、ｘ２：スチームタービンへの蒸気入力量（＝ｙ２）
ｙ３：スチームタービンの電気出力量
ｙ４：スチームタービンの熱出力量
ａ３，ｂ３，ａ４，ｂ４：機器特性パラメータ
【０１４７】
3)電気合計出力量（電気負荷６４）
電気合計出力量は、数式３５に示すように、ガスタービン６１及びスチームタービン６２で発生した電気出力合計値となる。
［数３５］
ｙ５＝ｆ５（ｙ１，ｙ３）＝ｙ１＋ｙ３
ここで、ｙ５：電気合計出力量
【０１４８】
4)熱合計出力量（熱負荷６３）
熱合計出力量は、スチームタービン６２で出力された熱であり、数式３６で表される。
［数３６］
ｙ６＝ｆ６（ｙ４）＝ｙ４
ここで、ｙ６：熱合計出力量
【０１４９】
(２)定常ネットワークシミュレータ
定常ネットワークシミュレータは、ガスタービン６１への燃料入力量ｘ１によって、ガスタービン６１とスチームタービン６２からの電気出力量ｙ１，ｙ３（すなわち電気合計出力量ｙ５）及びスチームタービン６２からの熱出力量ｙ４（すなわち熱合計出力量ｙ６）を得ることになる。
この例題では、全て線形モデルになっているため、収束計算は必要なく、以下のようなステップで入力から出力を得ることができる。
【０１５０】
・ステップ１：ガスタービンへの燃料入力量ｘ１と数式２９，３０から、ガスタービン出力量ｙ１，ｙ２を計算する。
・ステップ２：ステップ１で計算したガスタービン出力量ｙ２（＝ｘ２）と数式３３，３４より、スチームタービン出力量ｙ３，ｙ４を計算する。
・ステップ３：ステップ１，２で得られたガスタービンの電気出力量ｙ１及びスチームタービンの電気出力量ｙ３と数式３５を用いて、電気合計出力量ｙ５を計算する。
・ステップ４：ステップ２で得られたスチームタービンの熱出力量ｙ４と数式３６から熱合計出力量ｙ６を計算する。
【０１５１】
この実施例における定常ネットワークシミュレータの入力情報はガスタービンへの燃料入力量ｘ１であり、内部演算処理は構成機器の特性方程式である数式２９，３０，３３，３４、及び機器の接続状態を示す数式３５，３６を用いて行われ、出力情報は内部演算処理によって得られるガスタービン出力量ｙ１，ｙ２、スチームタービン出力量ｙ３，ｙ４、電気合計出力量ｙ５、熱合計出力量ｙ６となる。
【０１５２】
ここで、簡単な計算例としてガスタービンへの燃料入力量をｘ１＝５、数式２９，３０，３３，３４における各機器特性パラメータをａ１＝４，ａ２＝３，ａ３＝２，ａ４＝１，ｂ１＝１，ｂ２＝２，ｂ３＝３，ｂ４＝４とした場合の各出力状態量を、前記ステップに従い入力側から末端に向かい逐次計算する。
・ステップ１：ガスタービン出力量ｙ１，ｙ２の計算
ｙ１＝ａ１ｘ１＋ｂ１＝４・５＋１＝２１
ｙ２＝ａ２ｘ１＋ｂ２＝３・５＋２＝１７
・ステップ２：スチームタービン出力量ｙ３，ｙ４の計算
ｘ２＝ｙ２＝１７
ｙ３＝ａ３ｘ２＋ｂ３＝２・１７＋３＝３７
ｙ４＝ａ４ｘ２＋ｂ４＝１・１７＋４＝２１
・ステップ３：電気合計出力量ｙ５の計算
ｙ５＝ｙ１＋ｙ３＝２１＋３７＝５８
・ステップ４：熱合計出力量ｙ６の計算
ｙ６＝ｙ２＋ｙ４＝１７＋２１＝３８
【０１５３】
従って、定常ネットワークシミュレータにより計算される出力情報は、各ステップによって得られたガスタービン出力量ｙ１＝２１，ｙ２＝１７、スチームタービン出力量ｙ３＝３７，ｙ４＝２１、電気合計出力量ｙ５＝５８及び熱合計出力量ｙ６＝３８となり、この出力情報ｙ１〜ｙ６を最適化部に送ることとなる。
なお、ネットワークを構成する機器が増えた場合にも、各機器の特性方程式及び接続状態の関係式を用いて入力側から末端に向かい逐次計算することで、ネットワーク全体の各機器の入出力状態を計算することが可能である。
【０１５４】
(３)その他の計算条件
図１０の計測点ｍ１〜ｍ５における計測値を、それぞれ、ｍ１meas，ｍ２meas，ｍ３meas，ｍ４meas，ｍ５measとする。以下では、ｍｉmeasはｉ番目の計測値を意味するものとする。
また、この計測値ｍｉmeas(ｉ＝１〜５)は同時性が確保されており、過去の実績を考慮して事前に設定した計測範囲または各計測器の最大・最小計測範囲を逸脱する計測値および欠損している計測異常値はないものとする。
ｍｉmin≦ｍｉmeas≦ｍｉmax (ｉ＝１〜５)
ここで，ｍｉmin：ｉ番目の計測値の最小計測値
ｍｉmax：ｉ番目の計測値の最大計測値
【０１５５】
２．計算の実行
Ａ．計測補正値の計算
上述の計算条件を用いることにより、最も確からしい計測補正値を計算することが可能である。
(１)計測点における計算値の計算方法
それぞれの計測値は、数式３７〜４１のように計算することができる。
［数３７］
ｍ１calc＝ｘ１
［数３８］
ｍ２calc＝ｙ１
［数３９］
ｍ３calc＝ｙ３
［数４０］
ｍ４calc＝ｙ５
［数４１］
ｍ５calc＝ｙ６
ここで、ｍｉcalc：計測点ｍｉにおける計算値
【０１５６】
(２)最適化手法を用いた計測補正値の計算方法
機器の上下限制約などの機器制約も考慮した最適化の定式化は以下のようになる。
1)状態変数
状態変数は、ネットワークへの入出力エネルギーとなるため、ｘ１，ｙ５，ｙ６となる。
2)目的関数
目的関数は、数式４２に示すように、計測点における計測値と計算値との誤差の最小化となる。
【０１５７】
【数４２】
0004682553.tif 000025
【０１５８】
3)制約条件
ガスタービン６１及びスチームタービン６２の電気出力、蒸気出力に関する制約条件は、数式４３〜４６に示すとおりである。
【０１５９】
ａ）ガスタービンの電気出力制約
【０１６０】
【数４３】
0004682553.tif 000026
【０１６１】
ｂ）ガスタービンの蒸気出力制約
【０１６２】
【数４４】
0004682553.tif 000027
【０１６３】
ｃ）スチームタービンの電気出力制約
【０１６４】
【数４５】
0004682553.tif 000028
【０１６５】
ｄ）スチームタービンの蒸気出力制約
【０１６６】
【数４６】
0004682553.tif 000029
【０１６７】
つまり、数式４３〜数式４６の制約を満たす範囲で、数式４２の目的関数を最小化する状態変数（ネットワークへの入出力エネルギーであるｘ１，ｙ５，ｙ６）を求めることになり、その際、各計測値に対する計算値が数式２７〜数式４１の方程式から計算されて計測補正値として出力される。
【０１６８】
３．計測補正値を用いた計測異常値の判定
上述の方法により求めた計測補正値(計算値)における計測値に対する補正量を計算し、計測値に計測異常値が含まれているかを判定する。
(１)計測補正値の計測値に対する補正量の計算方法
計測値に対する補正量は、例えば次式のように計算することができる。
Ｃｏi＝｜ｍicalc−ｍimeas｜
ここでｋ、Ｃｏi：ｉ番目の計測値に対する補正量
(２)計測異常判定方法
上述で算出した補正量とネットワーク全体のエネルギーフローや計測器の公差等を用いて事前に設定した閾値とを比較し、閾値を超過する補正を行った計測値が事前に設定された判定数(例えば１つ)以上存在する場合に、計測値群に計測異常値が含まれていると判定する。
また、上記方法により計測値群に計測異常値が含まれていると判定された場合には、請求項２〜４に示した方法を用いて計測値群の中から計測異常値を求める必要がある。
【０１６９】
４．計測異常値の出力
請求項２〜４に示した方法を用いることにより、計測値の中から計測異常値を求めることができる。
（１)状態変数
各計測値が正常か異常かを示す配列を、計測値数に対応させてＴＲＵＥi(ｉ=１，……，５)とし、変数の中身は計測正常値の場合は１、計測異常値の場合は０のバイナリ値とする。
【０１７０】
(２)目的関数
目的関数は、数式４７のようになる。
【０１７１】
【数４７】
0004682553.tif 000030
【０１７２】
数式４７の目的関数を最小化する状態変数を求め、この状態変数において、異常値として指定されている計測値を計測異常値として出力する。
【実施例２】
【０１７３】
次に、本発明の実施例２について述べる。
図１１に、本実施例が適用されるコジェネレーションシステムを示す。図において、１１０はガスタービン、１１１は排ガスボイラ、１１２はスチームタービン、１１３はボイラ、１１４はコンプレッサ、１１５は吸収式冷凍機、１１６はガス焚冷凍機、１１７は蓄熱槽、１１８は空気負荷、１１９は電気負荷、１２０は蒸気負荷、１２１は熱負荷であり、それぞれ１つまたは複数台の機器が連なってプラントを構成している。また、データの計測点は各種プラント構成機器の入出力側であり、各機器の入出力状態量を計測することとなる。
【０１７４】
以下に、本実施例の条件について示す。
(条件１)
本実施例では、任意のネットワークシミュレータ入力状態量に対するネットワークシミュレータ計算結果を正しい計測値(真値)と仮定し、この計測値に任意の計測誤差および計測異常値を設定した計測値群を用いることとする。なお、ネットワークシミュレータ入力状態量としては、対象プラントにおける実計測値を用いることとする。
【０１７５】
(条件２)
本実施例において、各装置における各種計算に用いる計測値および計算値は、数式４８を用いて正規化を行うこととする。
［数４８］
Ｖnorm＝(Ｖ−Ｖmin)／(Ｖmax−Ｖmin)
ここで、Ｖ      ：計測値または計算値
Ｖmin      ：計測器の計測最小値
Ｖmax      ：計測器の計測最大値
Ｖnorm      ：各計測値の正規化後の値
(条件３)
また、各計測値に対する補正量は数式４９を用いて補正率として表現することとする。
［数４９］
Ｃｒ＝Ｃｏ／Ｖ×１００
ここで、Ｃｏ      ：補正量
Ｖ      ：計測値
Ｃｒ      ：補正率［％］
【０１７６】
(条件４)
初期計測異常出力装置７０における計測異常判定値は、各計測値を計測する計測器の計測可能範囲(計測最小値〜計測最大値)とする。つまり、精査装置入力時に計測可能範囲を逸脱する計測値を計測異常値として対象計測値群より除外することとする。
(条件５)
また、計測異常判定装置８０における計測異常判定値は補正率上限３％、判定数は１つとした。つまり、補正率が３％を超過する計測補正値が１つでも存在すれば、対象計測値群に計測異常値が含まれていると判定することとする。
【０１７７】
上記条件の下で、以下の２ケースにおける具体的な実施例について、図１２を参照しつつ説明する。
（ケース１）計測補正値出力例：計測値精査の対象となる計測値群に計測異常値がない場合
（ケース２）計測異常値出力例：計測値精査の対象となる計測値群に計測異常値が含まれる場合
【０１７８】
（ケース１）計測補正値出力の実施例
ここでは、前記条件１に従い、任意のネットワークシミュレータ計算値(真値)の各値に最大±3％のランダムな計測誤差を設定した計測値群を対象とする。
（１）初期計測異常出力ステップＳ１０
本実施例で対象とした計測値群の中には、初期計測異常出力装置７０によって計測異常と判定されて出力される計測異常値はない。なお、初期計測異常判定値(各計測器の計測可能範囲)を超過する計測値がある場合は、その計測値を計測異常値として計測値群から取り除く必要がある。
（２）計測補正値出力ステップＳ１１
前記初期計測異常出力ステップＳ１０により初期計測異常値が取り除かれた計測値群を対象とし、請求項１を用いて目的関数による評価値が最良(計測値とネットワークシミュレータ計算値の誤差が最小)となる計測補正値を求める。
【０１７９】
図１３に、計測補正値出力ステップＳ１１により得られた計測補正値の補正率分布を図に示す。また，補正率の最大・最小値および平均値は表１となる。
【０１８０】
【表１】
0004682553.tif 000031
【０１８１】
（３）計測異常判定ステップＳ１２
図１３および表１に示したように、このケースでは前記条件５に示した計測異常判定値である補正率３％を超過する計測補正値はないため、対象計測値群に計測異常値は含まれていないと判定し、前記計測補正値出力ステップＳ１１で得られたネットワークシミュレータ計算値を計測補正値として出力し、終了する。
【０１８２】
（ケース２）計測異常値出力の実施例
ここでは、前記条件１に従い任意のネットワークシミュレータ計算値(真値)の各値に最大±3％のランダムな計測誤差を設定し、また蒸気出力の計測値(図１１における計測点ｍST)に＋30％の計測異常値を設定した計測値群を対象とする。
（１）初期計測異常出力ステップＳ１０
本実施例で対象とした計測値群の中には初期計測異常出力装置によって計測異常と判定されて出力される計測異常値はない。なお，ケース１の場合と同様に初期計測異常判定値を超過する計測値がある場合はその計測値を計測異常値として計測値群から取り除く必要がある。
（２）計測補正値出力ステップＳ１１
ケース１の場合と同様に、前記初期計測異常出力ステップＳ１０により初期計測異常値が取り除かれた計測値群を対象とし、請求項１を用いて目的関数による評価値が最良(計測値とネットワークシミュレータ計算値の誤差が最小)となる計測補正値を求める。
【０１８３】
図１４に、計測補正値出力ステップＳ１１により得られた各計測補正値の補正率分布を図に示す。また、補正率の最大・最小値および平均値は表２となる。
【０１８４】
【表２】
0004682553.tif 000032
【０１８５】
（３）計測異常判定ステップＳ１２
図１４および表２に示すように、このケースでは前記条件５に示した計測異常判定値である補正率３％を超過する計測補正値が判定数(１つ)以上存在するため、対象計測値群に初期計測異常出力ステップＳ１０では取り除けない計測異常値が含まれていると判定する。
（４）計測異常値出力ステップＳ１３
前記計測異常判定ステップＳ１２により計測異常値が含まれると判定された計測値群を対象とし、請求項２〜４のいずれかの方法を用いることにより対象計測値群の中から計測異常値を求める。
具体的には、計測異常値出力装置４０において計測値状態として正常値及び異常値を指定し、正常値として指定された計測値のみを対象として、計測補正値出力装置３０を利用して計測値とネットワークシミュレータ計算値との誤差で表現される評価値を計算する。この評価値は、数式４７に示した計測異常値出力決定方法における全体目的関数の第１項である。この評価値と、別途計測異常値出力装置４０内で計算した全体目的関数の第２項である計測値が正常値である確率で表現される評価値とを合わせた総合評価値が最良となる事象(正常値及び異常値で表現される状態変数)において異常値として指定した計測値を、計測異常値として出力する。なお、本実施例においては前記条件１に従い設定した蒸気出力量の計測値(図１１における計測点ｍSTの計測値）を計測異常値として出力することとなる。
また、全体目的関数による評価値が最良となる事象における計測補正値出力装置３０において計算したネットワークシミュレータ計算値は、計測異常値を取り除いた計測値群に対する計測補正値となる。
【０１８６】
図１５に、計測異常値出力ステップＳ１３処理中の計測補正値出力装置３０により得られた各計測補正値の補正率分布を示す。また、補正率の最大・最小値および平均値は表３となる。
【０１８７】
【表３】
0004682553.tif 000033
【０１８８】
以上、ネットワークにおける計測値の計測補正値出力装置及び精査装置について説明した。本発明によって、従来技術では困難であったネットワーク上の計測値相互間の関連を考慮すると共に、ネットワークに関する方程式を利用し、各種の非線形な機器特性や数式化できない運用ルールなどを全て考慮した別プログラムとして構成される定常ネットワークシミュレータを用いて、メタヒューリスティク最適化手法（遺伝的アルゴリズムまたはその改良手法、タブサーチまたはその改良手法、Particle Swarm Optimizationまたはその改良手法）を利用し、過去の計測値の範囲や計測値相互間のつじつまを考慮した最も確からしい計測補正値及び計測異常値を求め、これらを出力可能となる。
【図面の簡単な説明】
【０１８９】
【図１】定常ネットワークシミュレータの一例を示す構成図である。
【図２】応答曲面法における中央複合計画の実験点を示す概念図である。
【図３】本発明の実施形態における計測補正値の計算概念を示す図である。
【図４】ＰＳＯにおける探索点の修正の概念を示す図である。
【図５】ＰＳＯにおける解探索の概念を示す図である。
【図６】数式９における重み係数の説明図である。
【図７】ＰＳＯの１回の最適化試行におけるアルゴリズムの概略を示す図である。
【図８】ＨＰＳＯのアルゴリズムの概略を示す図である。
【図９】ネットワーク計測値の精査装置を示す構成図である。
【図１０】本発明の実施例１が適用されるコジェネレーションシステムを示す図である。
【図１１】本発明の実施例２が適用されるコジェネレーションシステムを示す図である。
【図１２】本発明の実施形態における計測値の精査アルゴリズムの概略を示す図である。
【図１３】ケース１の計測補正値出力ステップにおける補正率分布を示す図である。
【図１４】ケース２の計測補正値出力ステップにおける補正率分布を示す図である。
【図１５】ケース２の計測異常値出力ステップ処理中の計測補正値出力装置における補正率分布を示す図である。
【符号の説明】
【０１９０】
１０：最適化部
２０：定常ネットワークシミュレータ
３０：計測補正値の出力装置
４０：計測異常値の出力装置
５０：データ出力装置
６０：データ入力装置
７０：初期計測異常出力装置
８０：計測異常判定装置
９０：データ保存・蓄積装置
１００：精査装置
２１：ガスタービンモデル
２２：排ガスボイラモデル
２３：スチームタービンモデル
２４：コンプレッサモデル
２５：ボイラモデル
２６：吸収式冷凍機モデル
２７：ガス焚冷凍機モデル
２８：蓄熱層モデル
１３１：空気負荷モデル
１３２：電気負荷モデル
１３３：蒸気負荷モデル
１３４：熱負荷モデル
６１：ガスタービン
６２：スチームタービン
６３：熱負荷
６４：電気負荷
１１０：ガスタービン
１１１：排ガスボイラ
１１２：スチームタービン
１１３：ボイラ
１１４：コンプレッサ
１１５：吸収式冷凍機
１１６：ガス焚冷凍機
１１７：蓄熱槽
１１８：空気負荷
１１９：電気負荷
１２０：蒸気負荷
１２１：熱負荷
【図１】
0004682553.tif 000034
【図２】
0004682553.tif 000035
【図３】
0004682553.tif 000036
【図４】
0004682553.tif 000037
【図５】
0004682553.tif 000038
【図６】
0004682553.tif 000039
【図７】
0004682553.tif 000040
【図８】
0004682553.tif 000041
【図９】
0004682553.tif 000042
【図１０】
0004682553.tif 000043
【図１１】
0004682553.tif 000044
【図１２】
0004682553.tif 000045
【図１３】
0004682553.tif 000046
【図１４】
0004682553.tif 000047
【図１５】
0004682553.tif 000048