(57)【要約】  【課題】機器特性式や定常ネットワークシミュレータを用いて、メタヒューリスティク最適化手法により、機器の最適パラメータを自動調整してそのトレンド等と一緒に出力可能とする。  【解決手段】ネットワークを構成する各機器に対して入出力されるエネルギーの計測値と計算値との誤差、パラメータ値の規定値との誤差を最小化する目的関数を用い、状態変数としての機器のパラメータ値をメタヒューリスティク最適化手法により求めて出力する最適化部１０と、機器モデル、機器運用ルール等を内部に有し、前記最適化部から入力されるパラメータ値のもとで、機器の出力量やネットワーク内の各部の状態量を求めて計測値に対する計算値として出力する機器特性式２０または定常ネットワークシミュレータ３０と、を備える。【選択図】  図１
2005050283.tif 000002
【特許請求の範囲】
【請求項１】
ネットワークを構成する各機器に対して入出力される電気・ガス・水・蒸気等のエネルギーの計測値と計算値との誤差、及び、各機器特性パラメータ値の規定値からの誤差を最小化する目的関数を用い、状態変数としての各機器特性パラメータ値をメタヒューリスティク最適化手法により求めて出力する最適化部と、
ネットワークを構成する各機器の機器モデル、機器運用ルール等を内部に有し、前記最適化部から入力される各機器特性パラメータ値のもとで、機器の入力エネルギーに対する出力エネルギーを求めて計測値に対する計算値として出力する機器特性式と、を備え、
最適化部に対象機器の入出力エネルギーの初期値として現在の計測値を入力すると共に、機器特性パラメータ値の規定値を入力し、機器特性式に現在の機器特性パラメータ値を入力して計測値に対する計算値を求め、
最適化部は、入力された計測値と機器特性式により求められた計算値とを用いてメタヒューリスティク最適化手法により機器特性パラメータ値を求める処理を繰り返し、その回数が事前の設定回数に達した時点のパラメータ値を機器特性式の最適パラメータ値として出力することを特徴とする機器特性パラメータ推定装置。
【請求項２】
ネットワークを構成する各機器に対して入出力される電気・ガス・水・蒸気等のエネルギーの計測値と計算値との誤差、及び、各機器特性パラメータ値の規定値からの誤差を最小化する目的関数を用い、状態変数としての各機器特性パラメータ値をメタヒューリスティク最適化手法により求めて出力する最適化部と、
ネットワークを構成する各機器の機器モデル、機器運用ルール、各機器の接続状態、及び、ネットワークのエネルギーフローに関する方程式等を内部に有し、前記最適化部から入力される各機器特性パラメータ値のもとで、定常状態のある時間断面におけるネットワーク内の各部の状態量を求めて計測値に対する計算値として出力する定常ネットワークシミュレータと、を備え、
最適化部にネットワークの入出力エネルギーの初期値として現在の計測値を入力すると共に、各機器特性パラメータ値の規定値を入力し、定常ネットワークシミュレータに現在の各機器特性パラメータ値を入力して計測値に対する計算値を求め、
最適化部は、入力された計測値と定常ネットワークシミュレータにより求められた計算値とを用いてメタヒューリスティク最適化手法により全ての機器の機器特性パラメータ値を求める処理を繰り返し、その回数が事前の設定回数に達した時点のパラメータ値を各機器特性式の最適パラメータ値として出力することを特徴とする機器特性パラメータ推定装置。
【請求項３】
請求項１または２に記載した機器特性パラメータ推定装置において、
最適化部によるメタヒューリスティク最適化手法として、遺伝的アルゴリズムまたはその改良手法を用いることを特徴とする機器特性パラメータ推定装置。
【請求項４】
請求項１または２に記載した機器特性パラメータ推定装置において、
最適化部によるメタヒューリスティク最適化手法として、タブサーチまたはその改良手法を用いることを特徴とする機器特性パラメータ推定装置。
【請求項５】
請求項１または２に記載した機器特性パラメータ推定装置において、
最適化部によるメタヒューリスティク最適化手法として、Particle Swarm Optimizationまたはその改良手法を用いることを特徴とする機器特性パラメータ推定装置。
【請求項６】
同時性が確保された計測値に基づいて請求項１〜５のいずれかの機器特性パラメータ推定装置により得られたパラメータ値、そのトレンド、及び、パラメータ値が急激に変化した機器に関する情報を生成して画面に表示出力させることを特徴とする機器特性パラメータ情報出力装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、電気・ガス・水・蒸気などが入出力されるネットワークを対象として、これらのネットワークを構成する各機器の特性式の最適パラメータ値を推定する機器特性パラメータ推定装置及び機器特性パラメータ情報出力装置に関するものである。
【背景技術】
【０００２】
電気・ガス・空調用の水（温水・冷水）・水道・蒸気などのネットワークの定常状態に対するソフトウェアによるシミュレーションでは、ネットワーク上の機器の特性式のパラメータを求め、ネットワークに対する入力データを与えることにより、ネットワーク内部の各部分に対する定常状態を計算によって求めることが可能となる。
上述したソフトウェアによるネットワークの定常状態のシミュレーションでは、従来、個別の機器毎に、機器に対する入力値及び出力値の計測値から、最小自乗法などによって個々の機器の特性式のパラメータを決定していた。このような方式では、個別の機器特性の誤差及び計測値誤差などの要因から、ネットワークの入力計測値と同時性を持った出力計測値、並びに、これらの計算値が合わない状況が発生する可能性がある。
【０００３】
このような状況に対しては、ネットワークに対する入力計測値・出力計測値を用いて、全てのネットワークの機器の特性式のパラメータを同時にチューニングする必要がある。
一方、プラントを循環する冷却系に関わる物理的または化学的な挙動を微分方程式で記述した場合については、特許文献１に記載された「プラントの自己学習診断，予測方法及び装置」が知られている。
【０００４】
【特許文献１】特開平６−２８９１７９号公報
【発明の開示】
【発明が解決しようとする課題】
【０００５】
ネットワークに対する入力計測値・出力計測値を用いて、全てのネットワークの機器の特性式のパラメータを同時にチューニングする場合、これは非線形最適化問題として定式化することができる。
従来、非線形最適化問題では、勾配法とその改良手法などの非線形計画法によって解が求められる。非線形計画法では、非線形最適化問題の場合、ネットワークのシミュレーションモデルも含めて数式化され、これを最適化プログラムの中の制約条件で考慮しながら解を生成する必要がある。このような状況において、目的関数の微分を計算し勾配方向を決定していく。つまり、対象問題は非線形最適化問題となるが、従来の非線形計画法では、モデルは最適化プログラムの中に記述されており、モデルが独立したプログラムとなっている場合には対応することができなかった。また、対応することができたとしても、ローカルな最適解（局所最適解）しか得ることができず、グローバルな最適解（大域最適解）を得ることができなかった。
【０００６】
前述した特許文献１記載の装置において、計測値とシミュレーションの値とが合致するようにパラメータを変更する部分は、離散的な感度解析を用いた一種の勾配法を用いており、対象となる最適化問題にローカルミニマムが存在する場合には初期パラメータ値から最適な値に収束しないという問題があり、しかも離散値のみで計算しているため、ローカルな最適解に収束することが不可能となる場合もある。
【０００７】
そこで本発明は、ネットワーク上の計測値間の関連を考慮しながらネットワークに関する方程式を利用し、非線形な各種の機器特性や数式化できない運用ルールなどを全て考慮した別プログラムとしてある定常ネットワークシミュレータを用いて、メタヒューリスティク最適化手法（遺伝的アルゴリズムまたはその改良手法、タブサーチまたはその改良手法、Particle Swarm Optimizationまたはその改良手法）を利用し、機器特性式のパラメータ値をチューニングして最適パラメータ値を推定するようにした機器特性パラメータ推定装置を提供することを目的とする。
また、パラメータ値自体やパラメータ値のトレンド、更には、パラメータ値が急速に変化した機器情報を画面に出力することにより、機器の計測器や付随する設備に対するCondition-based Maintenance（CBM）を実現可能とした機器特性パラメータ情報出力装置を提供することも目的とする。
【課題を解決するための手段】
【０００８】
上記課題を解決するため、請求項１に記載した機器特性パラメータ推定装置は、
ネットワークを構成する各機器に対して入出力される電気・ガス・水・蒸気等のエネルギーの計測値と計算値との誤差、及び、各機器特性パラメータ値の規定値からの誤差を最小化する目的関数を用い、状態変数としての各機器特性パラメータ値をメタヒューリスティク最適化手法により求めて出力する最適化部と、
ネットワークを構成する各機器の機器モデル、機器運用ルール等を内部に有し、前記最適化部から入力される各機器特性パラメータ値のもとで、機器の入力エネルギーに対する出力エネルギーを求めて計測値に対する計算値として出力する機器特性式と、を備え、
最適化部に対象機器の入出力エネルギーの初期値として現在の計測値を入力すると共に、機器特性パラメータ値の規定値を入力し、機器特性式に現在の機器特性パラメータ値を入力して計測値に対する計算値を求め、
最適化部は、入力された計測値と機器特性式により求められた計算値とを用いてメタヒューリスティク最適化手法により機器特性パラメータ値を求める処理を繰り返し、その回数が事前の設定回数に達した時点のパラメータ値を機器特性式の最適パラメータ値として出力するものである。
【０００９】
請求項２に記載した機器特性パラメータ推定装置は、
ネットワークを構成する各機器に対して入出力される電気・ガス・水・蒸気等のエネルギーの計測値と計算値との誤差、及び、各機器特性パラメータ値の規定値からの誤差を最小化する目的関数を用い、状態変数としての各機器特性パラメータ値をメタヒューリスティク最適化手法により求めて出力する最適化部と、
ネットワークを構成する各機器の機器モデル、機器運用ルール、各機器の接続状態、及び、ネットワークのエネルギーフローに関する方程式等を内部に有し、前記最適化部から入力される各機器特性パラメータ値のもとで、定常状態のある時間断面におけるネットワーク内の各部の状態量を求めて計測値に対する計算値として出力する定常ネットワークシミュレータと、を備え、
最適化部にネットワークの入出力エネルギーの初期値として現在の計測値を入力すると共に、各機器特性パラメータ値の規定値を入力し、定常ネットワークシミュレータに現在の各機器特性パラメータ値を入力して計測値に対する計算値を求め、
最適化部は、入力された計測値と定常ネットワークシミュレータにより求められた計算値とを用いてメタヒューリスティク最適化手法により全ての機器の機器特性パラメータ値を求める処理を繰り返し、その回数が事前の設定回数に達した時点のパラメータ値を各機器特性式の最適パラメータ値として出力するものである。
【００１０】
請求項３に記載した機器特性パラメータ推定装置は、請求項１または２において、最適化部によるメタヒューリスティク最適化手法として、遺伝的アルゴリズムまたはその改良手法を用いるものである。
【００１１】
請求項４に記載した機器特性パラメータ推定装置は、請求項１または２において、最適化部によるメタヒューリスティク最適化手法として、タブサーチまたはその改良手法を用いるものである。
【００１２】
請求項５に記載した機器特性パラメータ推定装置は、請求項１または２において、最適化部によるメタヒューリスティク最適化手法として、Particle Swarm Optimizationまたはその改良手法を用いるものである。
【００１３】
請求項６に記載した及び機器特性パラメータ情報出力装置は、同時性が確保された計測値に基づいて請求項１〜５のいずれかの機器特性パラメータ推定装置により得られたパラメータ値、そのトレンド、及び、パラメータ値が急激に変化した機器に関する情報を生成して画面に表示出力させるものである。
【発明の効果】
【００１４】
本発明によれば、メタヒューリスティク最適化手法を利用し、機器特性パラメータを自動調整して最適パラメータを推定することができる。また、パラメータ値が異常である機器を特定してそのパラメータ値自体やパラメータ値のトレンド、パラメータ値が急速に変化した機器情報を画面に出力することにより、Condition-based Maintenance（CBM）を実現可能とする。
【発明を実施するための最良の形態】
【００１５】
以下、図に沿って本発明の実施形態を説明する。
１．まず、本発明において使用される定常ネットワークシミュレータについて述べる。
定常ネットワークシミュレータは、ネットワーク及びその構成機器をモデル化し、各機器及びネットワークのエネルギーフローに関する方程式を解くことにより、ある時間断面におけるネットワーク全体の定常状態を計算するものである。
【００１６】
この種のネットワークシミュレータは、例えば、電力エネルギーに関しては、「電力系統解析理論」（関根泰次著，  電気書院，昭和46年1月）などに潮流計算技術として記述された既存技術であり、ガス流及び水道に関しては、「コンピュータ水理学入門」(C.A. ブレビア・A.J.フェラント著，磯部雅彦訳，サイエンス社，1983年)などに管路網解析技術として記述された既存技術である。更に、空調用の水(温水・冷水)の場合には、熱のバランス方程式となる。
エネルギーフロー方程式は、線形あるいは非線形連立方程式として定式化され、ニュートン法などを用いることにより、解、つまり、ある時間断面におけるネットワークの定常状態を得ることができる。
【００１７】
（１）定常ネットワークシミュレータに用いられる機器モデルは、ブランチとして表現され、その入出力関係を数式１の関数（機器特性式）により表現する。
【００１８】
【数１】
2005050283.tif 000003
【００１９】
ここで、機器の運用ルールも、上記関数ｆの中で表現される。なお、入力状態量はブランチに対する入力エネルギーから計算可能であり、出力状態量を用いてブランチから出力される出力エネルギーも計算可能である。また、送電線、水道管、ガス管などは一般的には抵抗分を含むため、この抵抗分によりエネルギー損失が生じるブランチも数式１により表現することとする。
【００２０】
（２）次に、定常ネットワークシミュレータに用いられるネットワークモデルについて述べる。ブランチである機器モデルを多数ノードにより接続したものがネットワークモデルとなる。各ノードにおいて、入力エネルギーと出力エネルギーとがバランスしている場合、以下の関数により表現することができる。
【００２１】
【数２】
2005050283.tif 000004
【００２２】
（３）ネットワークの定常シミュレーション結果は、数式２からなる非線形連立方程式を対象として、ネットワークに対する入力エネルギー及び出力エネルギーを指定して方程式を解くことによって得られる。なお、ネットワーク中の各機器は、数式１に示すように特性関数のパラメータを指定することによって入出力関係を表現可能であり、本発明では最終的にこのパラメータ値ａｊ（機器特性パラメータという）を調整するものである。
定常ネットワークシミュレータは、機器モデル及びネットワークのエネルギーフローに関する方程式を内部に有し、入力エネルギー及び出力エネルギーを指定することにより、ある時間断面におけるネットワーク内の各部の状態量を出力する。
【００２３】
２．次に、請求項１における、機器特性式の最適パラメータを推定するための定式化について述べる。
まず、状態変数としては、対象機器の特性式におけるパラメータ値とする。また、定式化に当たっては、対象機器の入力・出力状態量が計測(測定)できると仮定する。
【００２４】
最適パラメータを推定するための定式化として、数式３に示すように、目的関数は計測点における出力計測値と計算値とのノルム、及び、パラメータの規定値からの偏差のノルムの最小化とする。なお、数式４は制約条件を示す。
【００２５】
【数３】
2005050283.tif 000005
【００２６】
【数４】
2005050283.tif 000006
【００２７】
ここで、Ｗｍｉ１は後述の方法により決定する。ａｒｅｆｌ１は、経験上得られている特性式のパラメータやカタログ値から得られる特性式のパラメータとなる。また、ｗｐｊ１，Ｗ１，Ｗ２，ｐ１，ｑ１，ａｌ１，ｍｉｎ，ａｌ１，ｍａｘ等の値は事前に設定する。
最適化手法としては、計算対象とする全ての計測点の組合せの計測値と計測点における計算値とが合うように、また、各パラメータ値が規定値からできるだけ外れないように、各機器の特性式のパラメータを調整することとなる。
【００２８】
３．次に、本実施形態における最適パラメータの推定方法の基本的な概念について説明する。
本実施形態において、状態方程式ｈｉは対象機器の特性式から計算される。この特性式は、前述のように、機器特性、機器運用ルールなどを含み、機器の入力エネルギーが入力された時に出力エネルギーが計算される。この出力エネルギーを用いて、数式３により計測値及び計算値の評価が可能となる。
本発明では、数式４の制約条件を考慮した数式３の目的関数を最小化する最適化問題を、メタヒューリスティク最適化手法（遺伝的アルゴリズムまたはその改良手法、タブサーチまたはその改良手法、Particle Swarm Optimization或いはその改良手法）を利用して解くこととする。
【００２９】
（１）ここで、数式３における重み係数ｗｍｉ１の決定方法を以下に述べる。
ｗｍｉ１は、計測値の確からしさに対する各計測値の相対的な重みを表現している。この重みは、以下の２つの項目により評価できる。
1)計測センサーの精度(SensorPRE)
2)計測センサーの故障率(SensorFAU)
計測センサーの精度は、±５％などの％誤差範囲で表現される。また、故障率は、年間の故障回数(故障回数／年)で表現される。これらの値を用いて、各計測値に対する重み係数ｗｍｉ１を数式５のように決定する。
【００３０】
【数５】
2005050283.tif 000007
【００３１】
上記の重み係数ｗｍｉ１は、例えば数式６により表現可能である。なお、数式６におけるａ，ｂは事前に指定するものとする。
【００３２】
【数６】
2005050283.tif 000008
【００３３】
（２）次に、機器特性式のパラメータは、以下のアルゴリズムにより求める。
・ステップ１：対象機器の機器特性式情報、機器特性式の各パラメータに対する規定値、対象とする計測値の組合せを最適化部に入力する。
・ステップ２：現在の状態変数値（各機器特性式のパラメータ値）を機器特性、機器運用ルールなどを含んだ機器特性式に入力し、出力計測値に対する計算値を求める。
・ステップ３：最適化部は、ステップ１で入力された出力計測値とステップ２で計算された出力計測値に対する計算値、及び、数式３〜５を用いて、評価値を求める。
・ステップ４：メタヒューリスティク最適化手法（遺伝的アルゴリズムまたはその改良手法、タブサーチまたはその改良手法、あるいは、ＰＳＯまたはその改良手法）を用いて、数式３，４に基づき現在の状態変数値を変更する。
・ステップ５：事前に設定した探索回数に達したらステップ６へ進み、そうでない場合にはステップ２に戻る。
・ステップ６：現在のパラメータ値を最適パラメータ値として出力する。
【００３４】
ここで、図１は本実施形態における最適パラメータ値の計算概念を示す図である。図１において、１０はメタヒューリスティク最適化手法により計測値と計算値との誤差の最小化、状態変数（機器特性式のパラメータ値）とパラメータ規定値との誤差の最小化を目的として、現在の状態変数値を変更する（数式３，４を解く）最適化部、２０は最適化部１０からの状態変数（機器特性式のパラメータ値）を用い、機器特性や運用ルール等を考慮して入力された入力エネルギーに対する出力エネルギーを計算し、これを計算値として出力する機器特性式を示す。
【００３５】
４．次に、メタヒューリスティク最適化手法としての遺伝的アルゴリズムまたはその改良手法を用いた最適パラメータの推定方法を説明する。
ここで、遺伝的アルゴリズムまたはその改良手法とは、D. E. Goldberg，Genetic Algorithms in Search，Optimization，and Machine Learning，Addison-Wesley，1989の中で述べられているSimple Genetic Algorithm（以下、ＳＧＡという）とその改良手法を指す。以下では、このＳＧＡを例として最適パラメータの推定方法を説明する。
【００３６】
(１)状態変数の表現
対象とする機器特性式のパラメータを番号付けする。対象とするパラメータ数だけの遺伝子座をもつ遺伝子を用意する。ＳＧＡでは、全て離散量としなければならないことから、各遺伝子座に対応する各状態変数(パラメータ)は離散化する最小パラメータ値刻みで数式４の下限から上限の間の連続量を離散化した値を利用し、離散化した数値の最小値から最大値までを整数と対応させる。内部で整数で表現し、実際のパラメータ値へは対応表より変換する。つまり、状態変数は全ての状態を整数で表現し、パラメータの数の合計の長さを持つ遺伝子により表現する。
【００３７】
(２)求解アルゴリズム
・ステップ１：前提条件と初期条件設定
1)各計測値を入力する。
2)機器特性式の各パラメータに対する規定値を入力する。
3)状態変数(パラメータ)の上下限値を設定し、この上下限値から、状態変数を整数値に変換する表を作成する。
4)ストリング数、交差確率、突然変異確率、最大世代数を設定する。
【００３８】
・ステップ２：初期値の生成
1)各ストリングの遺伝子座について、上限及び下限の間の範囲の可能な整数値の中から、ランダムに整数値を生成する。
2)現在の世代数を１とする。
3)評価値の最良値(Obj*)を無限大とする。
【００３９】
・ステップ３：各ストリングの評価と選択
1)ストリングの遺伝子座の整数値と対応表を用いて各状態変数値(パラメータ値)を決定する。
2)機器特性式に入力エネルギー計測値を入力し、各出力エネルギー計測値に対する計算値を求め、この値を用いて数式３により評価値を計算する。
3)2)の評価値の逆数をFitness関数値として、ルーレット・ホイール・セレクションによりストリングの選択を行う。
4)選択されたストリングの評価値のうち、最小値が現在のObj*より小さかったら、その評価値をObj*とする。また、その時の状態変数の値をX*とする。更に、計測値に対する計算値をMeaCalc*とする。
【００４０】
・ステップ４：ストリング操作
ストリング集合に対して、交差確率、突然変異確率を用いて、交差及び突然変異を実行する。
・ステップ５：終了判定
世代が事前に決定した値に達したら終了し、現在のObj*，X*，MeaCalc*を出力する。このうち、X*が最適パラメータ値となる。
世代が事前に決定した値に達していない場合には、世代数に１足して、ステップ３へ戻る。
【００４１】
５．次いで、タブサーチまたはその改良手法を用いた最適パラメータの推定方法を説明する。
ここで、タブサーチまたはその改良手法とはF. Glover，"Tabu Search Part I"，ORSA Journal of Computing，Vol. 1，NO. 3，Summer 1989で述べられているタブーサーチ（以下、ＴＳという）とその改良手法のことを指す。以下、このＴＳを例として最適パラメータの推定方法を説明する。
【００４２】
(１)状態変数の表現
対象とする機器特性式のパラメータを番号付けする。対象とするパラメータ数だけの要素をもつ配列を用意する。ＴＳでは、全て離散量としなければならないことから、各要素に対応する各状態変数(パラメータ)は離散化する最小パラメータ値刻みで下限から上限の間の連続量を離散化した値を利用し、離散化した数値の最小値から最大値までを整数と対応させる。内部で整数で表現し、実際のパラメータ値へは対応表より変換する。つまり、全ての状態を整数で表現し、機器特性式のパラメータの数の合計の長さを持つ配列で表現する。
【００４３】
(２)求解アルゴリズム
・ステップ１：前提条件と初期条件設定
1)各計測値を入力する。
2)機器特性式の各パラメータに対する規定値を入力する。
3)状態変数(パラメータ)の上下限値を設定し、この上下限値から、状態変数を整数値に変換する表を作成する。
4)タブー長及び最大探索回数を設定する。
【００４４】
・ステップ２：初期値の生成
1)状態表現の配列の各要素について、対応する状態変数の上限及び下限の間の範囲の可能な整数値の中からランダムに整数値を生成し、現在状態とする。
2)現在状態をタブーリストに入れる。
3)現在の探索回数を１とする。
4)現在状態に対して、配列の整数を対応表を用いて状態変数(パラメータ)を決定する。機器特性式に、現在状態に対する入力エネルギー量を入力し、各出力エネルギー計測値に対する計算値を求め、この値を用いて数式３により評価値を計算する。
5)現在の状態に対する目的関数値を最良値(Obj*)とする。
【００４５】
・ステップ３：隣接状態の生成と次状態の決定
1)現在状態の各配列インデックスに対して整数値の±１の値（上下限値の場合は、限界値を超える方向には値を生成しない）を隣接状態として生成する。例えば、状態変数の数が３の場合、配列の要素は３つとなるが、この配列の値が（２，３，４）であった場合、第１要素の２の±１、第２要素の３の±１、第３要素の４の±１をそれぞれ個別に作成し、隣接状態を作成するため、隣接状態は、（１，３，４），（３，３，４），（２，２，４），（２，４，４），（２，３，３），（２，３，５）の６つとなる。
2)各隣接状態に対して、配列の整数を対応表を用いてパラメータ値を決定する。
3)機器特性式に各隣接状態に対するパラメータ値を入力し、各計測値に対する計算値を求め、この値を用いて数式３により評価値を計算する。
4)各隣接状態の中で、タブーでない最も評価がよいものを次状態とし、次状態の評価値が現在のObj*より小さかったら、その評価値をObj*とする。
また、その時の状態変数の値をX*とする。更に、計測値に対する計算値をMeaCalc*とする。
5)現在状態を次状態とする。
【００４６】
・ステップ４：終了判定
探索回数が事前に設定した値に達したら、終了し、現在のObj*，X*，MeaCalc*を出力する。このうち、X*が最適なパラメータ値となる。
探索回数が事前に設定した値に達しない場合には、現在の探索回数に１足してステップ３へ戻る。
【００４７】
６．次に、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を用いて、各エージェントの速度は以下の数式７により修正される。
【００５３】
【数７】
2005050283.tif 000009
【００５４】
つまり、最適化を実現する複雑な探索の方向は、Boydと同様にベクトルの合成で表現できる。また、数式７の右辺第２項は、pbestに収束する方向となっており、第３項は、同様にgbestに収束する方向となっている。つまり、この探索方向は、個人の経験と他人の経験を融合した方向となっている。また、数式７による探索方向の算出により、ＰＳＯはこれまでの速度を維持しようとする大域探索(数式７の右辺第１項)と、pbest，gbestを用いてそれに近づこうとする局所探索(数式７の右辺第２，３項)とをバランスよく行う機構をもった探索手法である。
上記の数式７を用いることにより、各エージェントのこれまでの最良解及び集団の最良解に確率的に近づくような速度が求められ、これにより各エージェントの現在の位置（探索点）を以下の数式８により修正する。
【００５５】
【数８】
2005050283.tif 000010
【００５６】
上述した探索点の修正の概念を図２に示す。この図２において、エージェントの現在の探索点ｓｋ(速度vorig)が、pbestに基づく速度ｖpbest及びgbestに基づく速度ｖgbest により、数式７，８に従って新たな探索点ｓｋ＋１に修正される。
ＰＳＯはＧＡ等と同様に複数の探索点を持った多点探索であり、各探索点のpbestと集団のgbestを用いて各探索点を確率的に変更していくことにより、大域最適解（最良解）を得る方法である。また、ＰＳＯは各ステップで目的関数値を評価する必要があるが、評価の回数は問題の規模によらずエージェント数のみで良いというメリットがある。従って、容易に大規模問題への適用が可能である。
図３に、多次元空間におけるＰＳＯによる解探索の概念を示す。この図は、ｎ次元の問題をｍ個のエージェントにより探索する場合を示している。
【００５７】
（３）ＰＳＯのパラメータ
ＰＳＯにおいて重要なパラメータは、試行回数、エージェント数、反復回数、数式７の各重み係数である。これらは全て設定値として変更可能とする。
1)試行回数
ＰＳＯは確率的な最適化手法であるため、通常１回の最適化のみでなく乱数の種（seed）を変更して複数回の最適化を行い、その結果の中から最も良い解を最適解とする。その最適化を実行する回数を試行回数とする。
【００５８】
2)エージェント数と反復回数
エージェント数は、何点の同時探索を行うかの数値である。反復回数は、各探索点（エージェント）が何回位置を更新して探索を行うかの数値である。つまり、ＰＳＯのアルゴリズムでは、各エージェントの目的関数値は反復回数分評価される。したがって、１回の最適化試行に対して、目的関数値の評価に必要な計算時間×エージェント数×反復回数の計算時間が必要となる。この計算時間を考慮しながら、これらのパラメータを決定する必要がある。
【００５９】
3)重み係数（数式７におけるw，c1，c2）
重み係数w，c1，c2は、経験的に以下の数式９，１０のように設定すると良いことが経験的に分かっている。
【００６０】
【数９】
2005050283.tif 000011
【００６１】
【数１０】
2005050283.tif 000012
【００６２】
数式９によって決定される重みｗは、図４に示すように、探索が進む（反復回数が増える）に従って小さくなることが分かる。これは、数式７において重みｗのかかる項は上述のように大域探索に相当する項であり、探索の開始時点では大域探索への比重を比較的大きくしておき、探索が進むに従ってその比重を小さくする、つまり局所探索の比重を大きくするような特性となっているためである。
【００６３】
（４）ＰＳＯによる最適化アルゴリズム
一回の最適化試行におけるアルゴリズムの概略を、図５に示す。
すなわち、エージェント数に応じた複数の初期状態をランダムに生成し（Ｓ１）、各エージェントに対する目的関数値を数式３によって計算する（Ｓ２）。
次に、各エージェントの探索過程で目的関数を評価し、最良解が得られた状態をｐｂｅｓｔとして保存する（Ｓ３）。その後、すべてのエージェントのｐｂｅｓｔを比較し、最良解が得られたものをｇｂｅｓｔとして保存する（Ｓ４）。
次いで、最大反復数に達したか否かを判断し、達した場合には現時点のｇｂｅｓｔ（最良解が得られたエージェント）の状態を最適解として出力する（Ｓ５ｙｅｓ，Ｓ７）。
また、最大反復数に達していない場合には、数式７，８により制約条件を満たす範囲で各エージェントの速度を計算し、探索点を移動してステップＳ２に戻る（Ｓ５ｎｏ，Ｓ６）。
【００６４】
Ｂ．Gbestモデルを用いた最適パラメータの推定方法
（１）状態変数の表現
対象とする機器の機器特性式のパラメータ値を番号付けする。対象とするパラメータ数だけの要素をもつ配列を用意する。ＰＳＯでは、状態変数は連続量として表現可能であり、配列の１つ１つの要素が各パラメータ値となっており、パラメータの数の合計の長さを持つ配列で表現する。
【００６５】
（２）求解アルゴリズム
・ステップ１：データの入力
1)各計測値を入力する。
2)機器特性式の各パラメータに対する規定値を入力する。
3)状態変数(機器特性式のパラメータ値)の上下限値を設定する。
4)エージェント数、各最適化パラメータ値、最大探索回数を設定する。
【００６６】
・ステップ２：初期値の生成
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)状態変数(機器特性式のパラメータ値)の上下限値を設定する。
4)エージェント数、各最適化パラメータ値、最大探索回数を設定する。
【００７５】
・ステップ２：初期値の生成
1)各エージェント毎に各状態変数値を、設定した上下限値内でランダムに生成する。
2)pbest及びgbestの初期設定
各エージェント毎のpbest値及びgbest値を、事前に設定した大きな値とする。
【００７６】
・ステップ３：各エージェントの評価
1)各エージェント毎の評価値の計算
現在の各状態変数(機器特性式のパラメータ値）を機器特性式に入力し、各出力エネルギー計測値に対する計算値を求める。この計算値を用いて数式３により評価値を計算する。
2)pbest及びgbestの更新
1)で計算した各エージェント毎の評価値が現在の各エージェント毎のpbest値より良かったら、現在の値をpbestに変更する。
上記で計算したpbestのうちの最良値が現在のgbestより良かったら、その値をgbest値に変更する。
【００７７】
・ステップ４：選択の実行
1)現在の各エージェントの評価値を、良い順番にソーティングする。
2)事前に設定した割合(Sr)を用いて、変更する個体数(Ns)を以下の数式１１により計算する。
［数１１］
Ns = Round(AgentNum×Sr)
3)最も評価の悪いエージェントからNs個だけの探索点及び速度を、最も評価の良いエージェントからNs個だけの探索点及び速度と変更する。この際、変更したエージェントのpbest情報はそのままとする。
【００７８】
・ステップ５：各エージェントの探索点の修正
設定した上下限値内で、各エージェント毎に状態変数値を数式７，８を用いて修正する。
【００７９】
・ステップ６：終了条件のチェック
探索回数が、入力した最大探索回数に達したら終了し、現状のgbestの評価値、状態変数値（最適パラメータ値）、及び、各計測値に対する計算値を出力する。
探索回数が最大探索回数に達しない場合は、ステップ３へ戻る。
【００８０】
７．次に、請求項２における、機器特性式の最適パラメータを推定するための定式化について述べる。つまり、ネットワーク内の全機器の機器特性式のパラメータを同時に最適化する方法に対する定式化を示す。
状態変数としては、対象機器の特性式におけるパラメータ値とする。また、定式化に当たっては、対象機器の入力・出力状態量が計測(測定)できると仮定する。
【００８１】
最適パラメータを推定するための定式化として、数式１２に示すように、目的関数は計測点における計測値と計算値とのノルム、及び、パラメータの規定値からの偏差のノルムの最小化とする。なお、数式１３は制約条件を示す。
【００８２】
【数１２】
2005050283.tif 000013
【００８３】
【数１３】
2005050283.tif 000014
【００８４】
ここで、数式１２におけるＷｍｉは、前述したＷｍｉ１と同様の方法により決定することができる。ａｒｅｆｌ１は、各機器の入力及び出力計測値から請求項１〜４の方法により得られる機器特性式のパラメータやカタログ値から得られる特性式のパラメータとなる。
また、ｗｐｊ，Ｗ１，Ｗ２，ｐ，ｑ，ａｌ，ｍｉｎ，ａｌ，ｍａｘ等の値は事前に設定する。
最適化手法としては、計算対象とする全ての計測点の組合せの計測値と計測点における計算値とが合うように、また、各パラメータ値が規定値からできるだけ外れないように、各機器の特性式のパラメータを調整することとなる。
【００８５】
８．次に、定常ネットワークシミュレータを用いた最適パラメータの推定方法について説明する。
この実施形態において、数式１２における状態方程式ｈｉは、定常ネットワークシミュレータにより独立に計算される。この定常ネットワークシミュレータは、前述のように、機器モデル、機器運用ルール、ネットワークの接続状態などを含み、ネットワークの入力・出力エネルギーが入力された時に定常状態のある時間断面におけるネットワークの各部の状態量が計算される。この状態量を用いて、数式１２により計測値及び計算値の評価が可能となる。
本実施形態では、数式１３の制約条件を考慮した数式１２の目的関数を最小化する最適化問題を、前記同様に各種のメタヒューリスティク最適化手法を利用して解くこととする。
【００８６】
（２）機器特性式のパラメータは、以下のアルゴリズムにより求める。
・ステップ１：ネットワーク内の機器接続情報、機器特性式情報、各機器の機器特性式の各パラメータに対する規定値、対象とする計測値の組合せを最適化部に入力する。
・ステップ２：現在の状態変数値（各機器特性式のパラメータ値）を定常ネットワークシミュレータに入力し、計測値に対する計算値を求める。
・ステップ３：最適化部は、ステップ１で入力された計測値とステップ２で計算された各計測値に対する計算値、及び、数式１２，１３を用いて、評価値を求める。
・ステップ４：メタヒューリスティク最適化手法（遺伝的アルゴリズムまたはその改良手法、タブサーチまたはその改良手法、あるいは、ＰＳＯまたはその改良手法）を用いて、現在の状態変数値を変更する。
・ステップ５：事前に設定した探索回数に達したらステップ６へ進み、そうでない場合にはステップ２に戻る。
・ステップ６：現在のパラメータ値を最適パラメータ値として出力する。
【００８７】
ここで、図７は本実施形態における最適パラメータ値の計算概念を示す図である。図７において、１０はメタヒューリスティク最適化手法により計測値と計算値との誤差の最小化、状態変数（機器特性式のパラメータ値）とパラメータ規定値との誤差の最小化を目的として、現在の状態変数値を変更する（数式１２，１３を解く）最適化部、３０は最適化部１０からの状態変数（各機器特性式のパラメータ値）を用い、機器特性や運用ルール等を考慮して入力された入・出力エネルギーに対するネットワーク各部の状態量を計算し、これを計算値として出力する定常ネットワークシミュレータを示す。
【００８８】
９．次に、メタヒューリスティク最適化手法として遺伝的アルゴリズムまたはその改良手法を用い、かつ、定常ネットワークシミュレータを用いた最適パラメータの推定方法を説明する。
なお、以下の大部分は前述の「４．」項と共通している。
また、遺伝的アルゴリズムまたはその改良手法の定義は、前述したものと同一である。
【００８９】
(１)状態変数の表現
対象とするネットワーク内の機器特性式のパラメータを番号付けする。対象とするパラメータ数だけの遺伝子座をもつ遺伝子を用意する。ＳＧＡでは、全て離散量としなければならないため、各遺伝子座に対応する各状態変数(パラメータ)は離散化する最小パラメータ値刻みで数式１３の下限から上限の間の連続量を離散化した値を利用し、離散化した数値の最小値から最大値までを整数と対応させる。内部で整数で表現し、実際のパラメータ値へは対応表より変換する。つまり、状態変数は全ての状態を整数で表現し、パラメータの数の合計の長さを持つ遺伝子により表現する。
【００９０】
(２)求解アルゴリズム
・ステップ１：前提条件と初期条件設定
1)各計測値を入力する。
2)機器特性式の各パラメータに対する規定値を入力する。
3)状態変数(パラメータ)の上下限値を設定し、この上下限値から、状態変数を整数値に変換する表を作成する。
4)ストリング数、交差確率、突然変異確率、最大世代数を設定する。
【００９１】
・ステップ２：初期値の生成
1)各ストリングの遺伝子座について、上限及び下限の間の範囲の可能な整数値の中から、ランダムに整数値を生成する。
2)現在の世代数を１とする。
3)評価値の最良値(Obj*)を無限大とする。
【００９２】
・ステップ３：各ストリングの評価と選択
1)ストリングの遺伝子座の整数値と対応表を用いて各状態変数値(パラメータ値)を決定する。
2)定常ネットワークモデルに入出力エネルギー量を入力し、各計測値に対する計算値を求め、この値を用いて数式１２により評価値を計算する。
3)2)の評価値の逆数をFitness関数値として、ルーレット・ホイール・セレクションによりストリングの選択を行う。
4)選択されたストリングの評価値のうち、最小値が現在のObj*より小さかったら、その評価値をObj*とする。また、その時の状態変数の値をX*とする。更に、計測値に対する計算値をMeaCalc*とする。
【００９３】
・ステップ４：ストリング操作
ストリング集合に対して、交差確率、突然変異確率を用いて、交差及び突然変異を実行する。
・ステップ５：終了判定
世代が事前に決定した値に達したら終了し、現在のObj*，X*，MeaCalc*を出力する。このうち、X*が最適パラメータ値となる。
世代が事前に決定した値に達していない場合には、世代数に１足して、ステップ３へ戻る。
【００９４】
１０．次いで、メタヒューリスティク最適化手法としてタブサーチまたはその改良手法を用い、かつ、定常ネットワークシミュレータを用いた最適パラメータの推定方法を説明する。
なお、以下の大部分は前述の「５．」項と共通している。
また、タブサーチまたはその改良手法の定義は、前述したものと同一である。
【００９５】
(１)状態変数の表現
対象とするネットワーク内の機器特性式のパラメータを番号付けする。対象とする全パラメータ数だけの要素をもつ配列を用意する。ＴＳでは、全て離散量としなければならないことから、各要素に対応する各状態変数(全パラメータ)は離散化する最小パラメータ値刻みで下限から上限の間の連続量を離散化した値を利用し、離散化した数値の最小値から最大値までを整数と対応させる。内部で整数で表現し、実際のパラメータ値へは対応表より変換する。つまり、状態変数は全ての状態を整数で表現し、機器特性式のパラメータの数の合計の長さを持つ配列で表現する。
【００９６】
(２)求解アルゴリズム
・ステップ１：前提条件と初期条件設定
1)各計測値を入力する。
2)機器特性式の各パラメータに対する規定値を入力する。
3)状態変数(パラメータ)の上下限値を設定し、この上下限値から、状態変数を整数値に変換する表を作成する。
4)タブー長及び最大探索回数を設定する。
【００９７】
・ステップ２：初期値の生成
1)状態表現の配列の各要素について、対応する状態変数の上限及び下限の間の範囲の可能な整数値の中からランダムに整数値を生成し、現在状態とする。
2)現在状態をタブーリストに入れる。
3)現在の探索回数を１とする。
4)現在状態に対して、配列の整数を対応表を用いて状態変数(パラメータ)を決定する。定常ネットワークシミュレータに、現在状態に対する入出力エネルギー量を入力し、各計測値に対する計算値を求め、この値を用いて数式１２により評価値を計算する。
5)現在の状態に対する目的関数値を最良値(Obj*)とする。
【００９８】
・ステップ３：隣接状態の生成と次状態の決定
1)現在状態の各配列インデックスに対して整数値の±１の値（上下限値の場合は、限界値を超える方向には値を生成しない）を隣接状態として生成する。例えば、状態変数の数が３の場合、配列の要素は３つとなるが、この配列の値が（２，３，４）であった場合、第１要素の２の±１、第２要素の３の±１、第３要素の４の±１をそれぞれ個別に作成し、隣接状態を作成するため、隣接状態は、（１，３，４），（３，３，４），（２，２，４），（２，４，４），（２，３，３），（２，３，５）の６つとなる。
2)各隣接状態に対して、配列の整数を対応表を用いてパラメータ値を決定する。
3)定常ネットワークシミュレータに各隣接状態に対するパラメータ値を入力し、各計測値に対する計算値を求め、この値を用いて数式１２により評価値を計算する。
4)各隣接状態の中で、タブーでない最も評価がよいものを次状態とし、次状態の評価値が現在のObj*より小さかったら、その評価値をObj*とする。
また、その時の状態変数の値をX*とする。更に、計測値に対する計算値をMeaCalc*とする。
5)現在状態を次状態とする。
【００９９】
・ステップ４：終了判定
探索回数が事前に設定した値に達したら、終了し、現在のObj*，X*，MeaCalc*を出力する。このうち、X*が最適なパラメータ値となる。
探索回数が事前に設定した値に達しない場合には、現在の探索回数に１足してステップ３へ戻る。
【０１００】
１１．次に、メタヒューリスティク最適化手法としてParticle Swarm Optimizationまたはその改良手法を用い、かつ、定常ネットワークシミュレータを用いた最適パラメータの推定方法を説明する。
なお、以下の大部分は前述の「６．Ｂ」項、「６．Ｄ」項と共通している。
また、Particle Swarm Optimizationまたはその改良手法の定義は、前述したものと同一である。
【０１０１】
Ａ．Gbestモデルを用いた最適パラメータの推定方法
（１）状態変数の表現
対象とする機器の機器特性式のパラメータ値を番号付けする。対象とするパラメータ数だけの要素をもつ配列を用意する。ＰＳＯでは、状態変数は連続量として表現可能であり、配列の１つ１つの要素が各パラメータ値となっており、パラメータの数の合計の長さを持つ配列で表現する。
【０１０２】
（２）求解アルゴリズム
・ステップ１：データの入力
1)各計測値を入力する。
2)機器特性式の各パラメータに対する規定値を入力する。
3)状態変数(機器特性式のパラメータ値)の上下限値を設定する。
4)エージェント数、各最適化パラメータ値、最大探索回数を設定する。
【０１０３】
・ステップ２：初期値の生成
1)各エージェント毎に各状態変数値を、設定した上下限値内でランダムに生成する。
2)pbest及びgbestの初期設定
各エージェント毎のpbest値及びgbest値を、事前に設定した大きな値とする。
【０１０４】
・ステップ３：各エージェントの評価
1)各エージェント毎の評価値の計算
現在の各状態変数(機器特性式のパラメータ値)を定常ネットワークシミュレータに入力し、各計測値に対する計算値を求める。この計算値を用いて数式１２により評価値を計算する。
2)pbest及びgbestの更新
1)で計算した各エージェント毎の評価値が現在の各エージェント毎のpbest値より良かったら、現在の値をpbestに変更する。
上記で計算したpbestのうちの最良値が現在のgbestより良かったら、その値をgbest値に変更する。
【０１０５】
・ステップ４：各エージェントの探索点の修正
設定した上下限値内で、各エージェント毎に状態変数値を数式７，８を用いて修正する。
【０１０６】
・ステップ５：終了条件のチェック
探索回数が、入力した最大探索回数に達したら終了し、現状のgbestの評価値、状態変数値（最適パラメータ値）、及び、各計測値に対する計算値を出力する。
探索回数が最大探索回数に達しない場合は、ステップ３へ戻る。
【０１０７】
Ｂ．ＨＰＳＯを用いた最適パラメータの推定方法
計測データ及びデータの仮定、定式化、状態変数の取り方は、ＰＳＯを用いた方法と同様である。ＨＰＳＯを用いた求解アルゴリズムは以下のようになる。
・ステップ１：データの入力
1)各計測値を入力する。
2)機器特性式の各パラメータに対する規定値を入力する。
3)状態変数(機器特性式のパラメータ値)の上下限値を設定する。
4)エージェント数、各最適化パラメータ値、最大探索回数を設定する。
【０１０８】
・ステップ２：初期値の生成
1)各エージェント毎に各状態変数値を、設定した上下限値内でランダムに生成する。
2)pbest及びgbestの初期設定
各エージェント毎のpbest値及びgbest値を、事前に設定した大きな値とする。
【０１０９】
・ステップ３：各エージェントの評価
1)各エージェント毎の評価値の計算
現在の各状態変数(機器特性式のパラメータ値）を定常ネットワークシミュレータに入力し、各計測値に対する計算値を求める。この計算値を用いて数式１２により評価値を計算する。
2)pbest及びgbestの更新
1)で計算した各エージェント毎の評価値が現在の各エージェント毎のpbest値より良かったら、現在の値をpbestに変更する。
上記で計算したpbestのうちの最良値が現在のgbestより良かったら、その値をgbest値に変更する。
【０１１０】
・ステップ４：選択の実行
1)現在の各エージェントの評価値を、良い順番にソーティングする。
2)事前に設定した割合(Sr)を用いて、変更する個体数(Ns)を前述の数式１１により計算する。
3)最も評価の悪いエージェントからNs個だけの探索点及び速度を、最も評価の良いエージェントからNs個だけの探索点及び速度と変更する。この際、変更したエージェントのpbest情報はそのままとする。
【０１１１】
・ステップ５：各エージェントの探索点の修正
設定した上下限値内で、各エージェント毎に状態変数値を数式７，８を用いて修正する。
【０１１２】
・ステップ６：終了条件のチェック
探索回数が、入力した最大探索回数に達したら終了し、現状のgbestの評価値、状態変数値（最適パラメータ値）、及び、各計測値に対する計算値を出力する。
探索回数が最大探索回数に達しない場合は、ステップ３へ戻る。
【０１１３】
１２．以下に、請求項６に記載した機器特性パラメータ情報出力装置について、図８を参照しつつ説明する。
この機器特性パラメータ情報出力装置８０は、同時性を確保した計測値が入力されると共に、請求項１〜８の何れかにより機器特性式の最適パラメータ値を計算し、このパラメータ値、パラメータ値トレンド、及び、パラメータ値が急激に変化した機器を表示装置９０の画面に表示出力するものである。
【０１１４】
図８において、出力装置８０には、ネットワーク中に設置された各計測値から同期の取れた計測値が入力される。この同期の取り方には、例えばGPSなどを利用し、事前に決定した周期（例えば１０分単位）で計測値と同時に時間タグを付けた計測データを装置に入力し、装置側でもGPSに合わせた計測時間を持ち、各計測機器の計測データのうち、この装置側の計測時間に一番近いデータをその時点での計測値とする方法が考えられる。
また、この方法はGPSを用いるのではなく、各計測器をネットワークで接続し、ＮＴＰサーバなどによって絶対時間の同期を取り、この絶対的な計測時間の時間タグと計測値を装置に入力し、同じ計測時間の計測値を同期のとれた計測値としてみる方法も考えられる。
【０１１５】
このようにして入力された計測値群は加工装置５０に入力され、加工装置５０では、事前に設定した周期毎の計測値の平均値を求め（図８の丸付き数字1）、同じ期間の計測値群を機器特性パラメータの出力装置４０に送る（図８の丸付き数字2）。ここで、機器特性パラメータの推定装置４０は、前述した最適化部１０及び機器特性式２０または定常ネットワークシミュレータ３０により構成されている。
推定装置４０では、入力された計測値群を用いて、事前に設定された周期毎に、請求項１〜５の何れかにより最適なパラメータ値を計算する（図８の丸付き数字3）。この計算されたパラメータ値は、機器パラメータ保存データベース７０及び機器パラメータ情報生成装置６０に送られる（図８の丸付き数字4，5）。
【０１１６】
機器パラメータ情報生成装置６０は、推定装置４０から送られた新たなパラメータ値と、データベース７０から送られた保存済み（過去）のパラメータ値（図８の丸付き数字6）とに基づいて、現在のパラメータ値、パラメータ値トレンド、及び、パラメータ値が急激に変化した機器に関する情報を生成し（図８の丸付き数字7）、表示装置９０に出力する。
以下では、この機器パラメータ情報生成装置６０における処理内容の例を説明する。
【０１１７】
（１）機器パラメータが急激に変化した状況の出力
機器パラメータ情報生成装置６０は、数式１４の関数値ｆ３を計算する。
【０１１８】
【数１４】
2005050283.tif 000015
【０１１９】
数式１４において、パラメータの現在値とは機器パラメータ保存データベース７０から得られる前回の各パラメータ値であり、パラメータの修正値とは推定装置４０から新たに得られた今回の各パラメータ値である。数式１４の関数値ｆ３が事前に設定した値以上になった場合には、パラメータ値が急激に変化したと判断し、機器パラメータ情報生成装置６０は警報を生成して表示装置９０に表示させる。
【０１２０】
（２）機器パラメータが急激に変化した機器情報の出力
上記の（１）により、パラメータ値が急激に変化したと判定された場合、機器パラメータ情報生成装置６０は、個々の機器について数式１５の関数値ｆ４を計算する。
【０１２１】
【数１５】
2005050283.tif 000016
【０１２２】
そして、各機器毎の上記関数値を大きい順番に並べ、各機器の種類（名称）及びその関数値を表示装置９０に表示させる。
【０１２３】
（３）パラメータ値のトレンド出力
機器パラメータ情報生成装置６０は、前記推定装置４０から新たなパラメータ値が送られてきたことをトリガーとして、事前に設定された時点数の過去の複数のパラメータ値をデータベース７０から入力し、全てのパラメータ値の過去から現在までのトレンド情報を生成して表示装置９０に表示させる。
【０１２４】
図９は、上述したパラメータ情報の処理シーケンスを示している。
すなわち、推定装置４０は所定周期で最適パラメータ値を計算し、機器パラメータ情報生成装置６０は、推定装置４０から送られた新たなパラメータ値をトリガーとしてパラメータ値のトレンドを画面出力させる（Ｓ１１）。次いで、パラメータ値に急激な変化があるか否かを前記数式１４により判断し（Ｓ１２）、急激な変化があった場合には前記数式１５による評価値（関数値ｆ４）を計算し、機器の種類（名称）及びその評価値を画面出力させる（Ｓ１２ｙｅｓ，Ｓ１３）。また、急激な変化がない場合には（Ｓ１２ｎｏ）ステップＳ１１に戻る。
【実施例】
【０１２５】
次に、本発明の実施例について述べる。
１．計算条件
図１０に、本実施例が適用されるコジェネレーションシステムを示す。図において、１０１はガスタービン、１０２はボイラ、１０３はスチームタービン、１０４は熱負荷、１０５は電気負荷であり、データの計測点としては１０個の計測点ｍ１〜ｍ１０を仮定している。
【０１２６】
(１)機器モデル
1)ガスタービン１０１
ここでは、排熱ボイラも含めてガスタービン６１としてモデル化しており、ガスタービン１０１は、燃料を入力することにより、電気と蒸気を出力するモデルとなっている。ガスタービン１０１の入出力関係は数式１６，１７に示す通りであり、入出力量は数式１８，１９に示すように各々線形形式で表現されているとする。
［数１６］
ｙ１＝ｆ１（ｘ１）
［数１７］
ｙ２＝ｆ２（ｘ１）
［数１８］
ｙ１＝ａ１ｘ１＋ｂ１
［数１９］
ｙ２＝ａ２ｘ１＋ｂ２
ここで、ｘ１：ガスタービンへの燃料入力量
ｙ１：ガスタービンの電気出力量
ｙ２：ガスタービンの蒸気出力量
ａ１，ｂ１，ａ２，ｂ２：ガスタービンの機器特性パラメータ
【０１２７】
2)スチームタービン１０３
スチームタービン１０３は、蒸気を入力することにより電気と熱(蒸気または温水)を出力するモデルとなっている。スチームタービン１０３の入出力関係は数式２０，２１に示す通りであり、入出力量は数式２２，２３に示すように各々線形形式で表現されているとする。
［数２０］
ｙ３＝ｆ３（ｘ２）
［数２１］
ｙ４＝ｆ４（ｘ２）
［数２２］
ｙ３＝ａ３ｘ２＋ｂ３
［数２３］
ｙ４＝ａ４ｘ２＋ｂ４
ここで、ｘ２：スチームタービンへの蒸気入力量（＝ｙ２）
ｙ３：スチームタービンの電気出力量
ｙ４：スチームタービンの熱出力量
ａ３，ｂ３，ａ４，ｂ４：スチームタービンの機器特性パラメータ
【０１２８】
3)ボイラ１０２
ボイラ１０２は、燃料を入力することにより、熱(蒸気または温水)を出力するモデルとなっている。ボイラ１０２の入出力関係は数式２４に示す通りであり、入出力量は数式２５に示すように線形形式で表現されているとする。
［数２４］
ｙ５＝ｆ５（ｘ３）
［数２５］
ｙ５＝ａ５ｘ３＋ｂ５
ここで、ｘ３：ボイラへの燃料入力量
ｙ５：スチームタービンの熱出力量
ａ５，ｂ５：ボイラの機器特性パラメータ
【０１２９】
4)電気合計出力量（電気負荷１０５）
電気合計出力量は、数式２６に示すように、ガスタービン１０１及びスチームタービン１０３で発生した電気出力合計値となる。
［数２６］
ｙ７＝ｆ７（ｙ１，ｙ３，ｙ６）＝ｙ１＋ｙ３＋ｙ６
ここで、ｙ６：電気購入分
ｙ７：電気合計出力量
【０１３０】
5)熱合計出力量（熱負荷１０４）
熱合計出力量は、スチームタービン１０３及びボイラ１０２で出力された熱であり、数式２７で表される。
［数２７］
ｙ８＝ｆ８（ｙ４，ｙ５）＝ｙ４＋ｙ５
ここで、ｙ８：熱合計出力量
【０１３１】
上記実施例では、ａ１〜ａ５，ｂ１〜ｂ５が最適化するべきパラメータとなる。なお、ここでは簡略化のため、冷却塔、復水器、冷凍機等の機器は無視している。
【０１３２】
(２)定常ネットワークシミュレータ
定常ネットワークシミュレータでは、ガスタービン１０１への燃料入力量ｘ１によって、ガスタービン１０１及びスチームタービン１０３からの電気出力量ｙ１，ｙ３が得られ、これらと電気購入量ｙ６とを合計して電気負荷１０５に対する電気合計出力量ｙ７を得る状態を模擬する。また、スチームタービン６２からの熱出力量ｙ４とボイラ１０２からの熱出力量ｙ５とを合計して熱負荷１０４に対する熱合計出力量ｙ８を得る状態を模擬する。
この例題では、全て線形モデルになっているため、収束計算は必要なく、以下のようなステップで入力から出力を得ることができる。
【０１３３】
・ステップ１：ガスタービンへの燃料入力量ｘ１と数式１８，１９から、ガスタービン出力量ｙ１，ｙ２を計算する。
・ステップ２：ステップ１で計算したガスタービン出力量ｙ２と数式２２，２３より、スチームタービン出力量ｙ３，ｙ４を計算する。
・ステップ３：ボイラ１０２への燃料入力量ｘ３と数式２５から、ボイラ出力量ｙ５を計算する。
・ステップ４：ステップ１，２で得られたガスタービンの電気出力量ｙ１、スチームタービンの電気出力量ｙ３、電気購入量ｙ６、及び数式２６を用いて、電気合計出力量ｙ７を計算する。
・ステップ５：ステップ２で得られたスチームタービンの熱出力量ｙ４とステップ３で得られたボイラの熱出力量ｙ５と数式２７から、熱合計出力量ｙ８を計算する。
【０１３４】
(３)計測点及び計測値、計算値
図１０の計測点ｍ１〜ｍ１０における計測値をｍ１meas〜ｍ１０measとし、ｍｉmeasはｉ番目の計測値を意味するものとする。また、各計測点ｍ１〜ｍ１０における計算値をｍ１calc〜ｍ１０calcとし、ｍｉcalcはｉ番目の計算値を意味するものとする。
各計測値は、数式２８〜３７のように計算することができる。
［数２８］
ｍ１calc＝ｘ１
［数２９］
ｍ２calc＝ｘ３
［数３０］
ｍ３calc＝ｙ６
［数３１］
ｍ４calc＝ｙ１
［数３２］
ｍ５calc＝ｙ２
［数３３］
ｍ６calc＝ｙ５
［数３４］
ｍ７calc＝ｙ３
［数３５］
ｍ８calc＝ｙ７
［数３６］
ｍ９calc＝ｙ４
［数３７］
ｍ１０calc＝ｙ８
【０１３５】
２．計算の実行
Ａ．機器特性式の最適パラメータの機器ごとの推定方法
上述した計算条件を用い、請求項１，３〜５の推定装置によって機器特性式の最適パラメータを推定することができる。ここでは、例としてガスタービン１０１の最適パラメータの推定につき説明する。他の機器についても、同様の方法で最適パラメータを推定可能である。
機器の上下限制約などの機器制約も考慮した最適化の定式化は、以下のようになる。
【０１３６】
（１）状態変数
状態変数は、機器特性式のパラメータとなるため、ガスタービン１０１の機器特性パラメータａ１，ｂ１，ａ２，ｂ２である。
【０１３７】
（２）目的関数
目的関数は、数式３８に示すように、ガスタービン１０１の出力に関する計測値と計算値との誤差の最小化、及び、パラメータの規定値との誤差の最小化となる。
【０１３８】
【数３８】
2005050283.tif 000017
【０１３９】
（３）制約条件
ガスタービン１０１の電気出力、蒸気出力に関する制約条件は、数式３９，４０に示すとおりである。
【０１４０】
1)ガスタービンの電気出力制約
【０１４１】
【数３９】
2005050283.tif 000018
【０１４２】
2)ガスタービンの蒸気出力制約
【０１４３】
【数４０】
2005050283.tif 000019
【０１４４】
従って、計測値ｍ１，ｍ４，ｍ５の組合せをNumComb数だけ入力し、これらの計測値を用いて数式３９，４０の制約を満たす範囲で数式３８の目的関数を最小化するパラメータ値を請求項３〜５の方法により求めることになる。その際、各計測値に対する計算値は、数式１６〜３７のネットワークに関する方程式から計算される。
実際には、例えばボイラが複数台あった場合は、通常台数制御がされており、この台数制御のロジックを入れたボイラモデルを考えると数式によって単純に表現することはできず、if-thenルールなどを利用した記述になってしまう。
このような場合、従来の最小自乗法と数理計画法によって解を得ることは困難であり、その他の運用ルールなども考慮すると、また、シミュレータを独立したプログラムとして考え、この中の機器特性を独立にチューニングすることを考えると、ここで提案しているような方法により最適パラメータを計算する方法が必要であることがわかる。
【０１４５】
Ｂ．機器特性式の最適パラメータの同時推定方法
前述した計算条件を用いることにより、定常ネットワークシミュレータを用いて、請求項２，３〜５の推定装置により全ての機器の機器特性式の最適パラメータ値を同時に決定することが可能である。
機器の上下限制約などの機器制約も考慮した最適化の定式化は、以下のようになる。
【０１４６】
(１)状態変数
状態変数は全ての機器特性式のパラメータ値となるため、ａ１〜ａ５，ｂ１〜ｂ５となる。
【０１４７】
(２)目的関数
目的関数は、数式４１に示すように計測点における計測値と計算値の誤差の最小化となる。ここでは、例として２次ノルム(ユークリッド距離)を利用する。
【０１４８】
【数４１】
2005050283.tif 000020
【０１４９】
(３)制約条件
ガスタービン１０１の制約条件は、数式３９，４０である。その他、スチームタービン１０３及びボイラ１０２の制約条件は、数式４２〜４４となる。
【０１５０】
1)スチームタービンの電気出力制約
【０１５１】
【数４２】
2005050283.tif 000021
【０１５２】
2)スチームタービンの蒸気出力制約
【０１５３】
【数４３】
2005050283.tif 000022
3)ボイラの蒸気出力制約
【０１５４】
【数４４】
2005050283.tif 000023
【０１５５】
従って、計測値ｍ１〜ｍ１０の組合せをNumComb数だけ入力し、これらの計測値を用いて数式３９，４０，４２〜４４の制約を満たす範囲で、数式４１の目的関数を最小化するパラメータ値を請求項３〜５の方法により求めることになる。その際、各計測値に対する計算値は、数式１６〜３７のネットワークに関する方程式から計算される。
【０１５６】
図８に示した機器特性パラメータ情報出力装置８０には、図９のコジェネレーションシステムの１０個の計測値ｍ１〜ｍ１０を入力する。ここでは、調整の周期を３０分毎とする。
図８における計測情報の加工装置５０には、同期がとれた計測値群が入力され、３０分毎の各計測値の平均値が求められ、例えば午前１０時〜午前１０時３０分の間の平均計測値が機器特性パラメータの推定装置４０に送られる。
機器特性パラメータの推定装置４０は、与えられた計測値を用いて、最適化部１０及び機器特性式２０または定常ネットワークシミュレータ３０により最適な機器パラメータ値を計算する。
【０１５７】
計算された最適パラメータ値及びその計測対象時間（上記の例の場合には、日付と午前１０時〜午前１０時３０分）は、機器パラメータ保存データベース７０及び機器パラメータ情報生成装置６０に送られる。
機器パラメータ保存データベース７０は、送られてきた最適パラメータ値とその計測対象時間とを保存する。機器パラメータ情報生成装置６０は、機器特性パラメータの推定装置４０から送られてきた最適パラメータ値及びその計測対象時間と、機器パラメータ保存データベース７０から送られてきた過去の最適パラメータ値及びその計測対象時間に基づき、以下の情報を表示装置９０に出力する。
【０１５８】
1)ガスタービン１０１に関するパラメータａ１，ｂ１，ａ２，ｂ２の最新のトレンド
2)スチームタービン１０３に関するパラメータａ３，ｂ３，ａ４，ｂ４の最新のトレンド
3)ボイラ１０２に関するパラメータａ５，ｂ５の最新のトレンド
【０１５９】
また、ここで対象となる１０個のパラメータａ１，〜ａ５，ｂ１〜ｂ５に対し、数式１４を３０分毎に計算し、パラメータ値の急激な変化があるかどうかを判定する。そして数式１４の関数ｆ３の値が事前に設定した値以上になった場合は、画面上に「機器パラメータが急激な変化をした」というメッセージを出力する。また、このような状況になった場合、個々の機器について数式１５の関数ｆ４を計算する。つまり、以下の計算を行う。
【０１６０】
1)ガスタービン１０１に関するパラメータａ１，ｂ１，ａ２，ｂ２の評価値（関数ｆ４）
2)スチームタービン１０３に関するパラメータａ３，ｂ３，ａ４，ｂ４の評価値（関数ｆ４）
3)ボイラ１０２に関するパラメータａ５，ｂ５の評価値（関数ｆ４）
【０１６１】
そして、機器パラメータ情報生成装置６０は、これらの評価値が大きい順番にソートして表示装置９０の画面に表示出力させる。表示例としては、例えば以下のようになる。
【０１６２】
「機器特性が急激に変化した機器は、変化が大きい順番に以下のとおり
評価値
1)ボイラ                100.0
2)スチームタービン       90.0
3)ガスタービン           75.0」
【０１６３】
この表示を見た運転員は、上記の評価値が大きい機器について、計測値及びその計測値の通信に関する設備などを点検するというCondition-based Maintenance (CBM)を実現することができる。
【図面の簡単な説明】
【０１６４】
【図１】本発明の実施形態における最適パラメータの計算概念を示す図である。
【図２】ＰＳＯにおける探索点の修正の概念を示す図である。
【図３】ＰＳＯにおける解探索の概念を示す図である。
【図４】数式９における重み係数の説明図である。
【図５】ＰＳＯの１回の最適化試行におけるアルゴリズムの概略を示す図である。
【図６】ＨＰＳＯのアルゴリズムの概略を示す図である。
【図７】本発明の実施形態における最適パラメータの計算概念を示す図である。
【図８】機器特性パラメータ情報出力装置の説明図である。
【図９】パラメータ情報の処理シーケンスを示す図である。
【図１０】本発明の実施例が適用されるコジェネレーションシステムを示す図である。
【符号の説明】
【０１６５】
１０：最適化部
２０：機器特性式
３０：定常ネットワークシミュレータ
４０：機器特性パラメータの推定装置
５０：計測情報の加工装置
６０：機器パラメータ情報生成装置
７０：機器パラメータ保存データベース
８０：機器特性パラメータ情報出力装置
９０：表示装置
１０１：ガスタービン
１０２：ボイラ
１０３：スチームタービン
１０４：熱負荷
１０５：電気負荷
【図１】
2005050283.tif 000024
【図２】
2005050283.tif 000025
【図３】
2005050283.tif 000026
【図４】
2005050283.tif 000027
【図５】
2005050283.tif 000028
【図６】
2005050283.tif 000029
【図７】
2005050283.tif 000030
【図８】
2005050283.tif 000031
【図９】
2005050283.tif 000032
【図１０】
2005050283.tif 000033