JP5457737B2 - プラント制御情報生成装置及び方法、並びにそのためのコンピュータプログラム - Google Patents

Description

この発明は、化学プラント、半導体製造設備、製鉄プラント等の制御技術に関し、特に、過去の複数の稼動時の稼動条件とその結果とに基づいて、プラントの制御の精度を高めるための技術に関する。

化学プラントでは、複数の装置（反応炉等）を使用して同一の製品を生産することが多い。この場合、生産される製品が均一となるように、各装置の稼動条件を適切に設定することが必要である。しかし、装置内におけるプロセスの進行とともにプロセスの特性が変化すること、及び複数の装置を並列で運転している場合に、それぞれの装置のカタログスペックが同一でも実際には個体差があること等により、生産される製品が均一となるような条件を装置ごとに精度よく決定したり、生産される製品の組成を稼動条件から精度よく予測したりすることはむずかしい。同様の問題は、化学プラントに限らず、半導体製造設備、金属精錬プラント、水処理施設、及びガスプラント等でも生じる。また、何らかのものを製造する施設だけではなく、発電・変電施設のように非有体物の特性を変化させるための施設でもこうした問題は生じ得る。以下、本明細書では、何らかのものを製造する施設だけでなく、非有体物の特性を変化させる施設をも含めて「プラント」と呼ぶ。また、複合的な施設だけでなく、単独の設備であって、稼動条件によって異なる挙動を示すものもここでいう「プラント」に含めるものとする。

これまで、本願発明者は、プロセスの特性変化が測定パラメータ間の相関関係の違いとして表現されることに着目し、パラメータ間の相関関係を考慮したモデリングがソフトセンサ設計に有効であることを示してきた。本願発明者はさらに、装置間の個体差に対応できる手法として、パラメータ間の相関関係の類似度を考慮して、過去の装置の稼動パラメータをクラスタリングする手法（向井法）と、相関識別法（ＮｅａｒｅＳｔＣｏｒｒｅｌａｔｉｏｎＭｅｔｈｏｄ；ＮＣ法）とを開発した（非特許文献１、２）。

向井、藤原、加納、長谷部：変数間の相関に着目したクラスタリング手法の開発およびその多変量統計モデリングへの利用、計測自動制御学会第8回制御部門大会、京都大学（２００８） 藤原、加納、長谷部：相関関係を考慮したパターン認識手法の開発とソフトセンサ設計への適用、計測自動制御学会第9回制御部門大会、広島大学（２００９）

向井法は、パラメータ間の相関関係を指標としてデータをクラスタリングできるが、繰返計算を行なう必要がある。一方、ＮＣ法は、線形空間におけるクエリとサンプルの幾何学的な関係に着目した手法であり、教師信号を用いることなくクエリと類似の相関関係を有するサンプルを検出できるという特徴を持つが、クラスタリングに用いることはできなかった。

ここでいうパラメータ間の相関関係とは、以下のようなものをいう。以下、説明を理解しやすくするために、２つのパラメータのみからなるサンプルを扱う。図１を参照して、互いに相関関係を持つサンプルの集合が３つあるものとする。図１においては、第１の集合は黒丸で、第２の集合は白丸で、第３の集合は「＋」記号で、それぞれ表されている。ここでは、各集合を「クラス」と呼ぶ。各クラスに属するサンプルの間の相関関係はクラスごとに異なっている。あるクラスについて新たなサンプルが得られれば、新たなサンプルを使用して、そのクラスに適した制御を行なうことができる。したがって、ある新たなサンプルが得られたときに、そのサンプルを的確なクラスに分類することが必要である。

通常、このような場合には、予めサンプルの集合をクラスに分類しておき、各クラスとサンプルとの間の「距離」を定義しておく手法が用いられる。新たなサンプルが得られると、そのサンプルと各クラスとの間の「距離」を算出して、最も距離が小さなクラスに新たなサンプルを分類する。

得られたサンプルの集合がどのクラスに属するかが明確であれば、上記した手法を用いることができる。しかし実際のプラントの稼動時に得られる実測サンプルがどのクラスに属するかは明らかでない。つまり、実際のプラントから得られるサンプルのクラスを図１のように３つのクラスに分類（クラスタリング）する手法が必要である。

クラスタリング手法として最も頻繁に使用されるものにｋ−平均法がある。図１に示されるサンプルの集合をｋ−平均法を用いてクラスタリングした結果を図２に示す。図２では、クラスタリングによりサンプルを３つのクラスに分類している。各クラスに分類されたサンプルをそれぞれ黒丸、白丸及び「＋」記号で表している。

図２により明らかなように、ｋ−平均法を用いた場合、サンプルを正確にクラスタリングすることはできない。これはｋ−平均法が基本的にサンプルの間の距離に基づいてクラスタリングを行なっているためである。すなわち、各クラスのサンプルが集中する中央部付近ではサンプルが全て同じクラスタに分類され、その両側に散在するサンプルはそれぞれ別々のクラスタに分類される。

ｋ−平均法を用いたクラスタリングとは別のクラスタリング手法も存在する。図３は、スペクトラル・クラスタリングと呼ばれるクラスタリング手法を用いて図１に示されるサンプルをクラスタリングした結果を示す。このスペクトラル・クラスタリングの手法についてはその詳細を後述する。図３において、各クラスに属するサンプルの表示方法は図１及び図２と同様である。このスペクトラル・クラスタリングは、後述するようにグラフ理論に基づくものであって、ｋ−平均法と比較すると単純に距離に基づいてクラスタリングするものではない。したがって、クラスタリングをより正確に行なえる可能性がある。しかし、図３に示される例ではクラスタリングが正しく行なわれているとはいえない。

過去の実測データを統計処理することにより、プラントの制御を行なう場合、実測データを正しくクラスタリングできなければ、制御の前提が失われることになり、正しい制御を実現できない。したがって、実測データを正しくクラスタリングすることができる技術が必要である。

したがって本発明の目的は、サンプル間の相関関係を用いて、スペクトラル・クラスタリング手法によりクラスタリングしたサンプルの集合を用いて、精度よくプラントの制御を行なうことができるプラント制御情報生成装置及び方法、並びにそのためにコンピュータプログラムを提供することである。

本発明に係るプラント制御情報生成装置は、プラントの稼動条件と、当該稼動条件でのプラントの過去の稼動結果から得られるデータ項目の実測値との組合せからなる実測サンプルを複数個記憶するためのデータベースと、データベースに記憶された複数個の実測サンプルを、複数個のクラスにクラスタリングするためのクラスタリング手段と、クラスタリング手段によりクラスタリングされた各クラスに対して、プラントにおける稼動条件とデータ項目の値との関係を示す統計的モデルを構築するためのモデル構築手段と、新たな稼動条件が与えられたことに応答して、統計的モデルと、新たな稼動条件とに基づいて、プラントに関する制御情報を算出するための制御情報算出手段とを含む。クラスタリング手段は、複数個の実測サンプルの任意の２個の間の相関の程度を成分とする行列を、データベースに記憶された実測サンプルから算出するための行列算出手段と、行列算出手段により算出された行列を類似度行列として、複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、複数個の実測サンプルを複数個のクラスにクラスタリングするための手段とを含む。

このプラント制御情報生成装置においては、クラスタリング手段が、データベースに記憶された複数個の実測サンプルを複数個のクラスにクラスタリングし、クラスタリングされた後のクラスごとに、統計的モデルが算出される。クラスタリングの際、複数個の実測サンプルの任意の２個の間の相関の程度を成分とする行列を類似度行列として、スペクトラル・クラスタリングが行なわれる。サンプルの間の相関の相違に基づいて複数個のサンプルがクラスタリングされるため、単純にサンプル間の距離に基づいてクラスタリングされる場合と比較して、異なる相関関係にあるサンプルが正しく別々のクラスにクラスタリングされる確率が高くなる。その結果、正しいクラスタリングが可能になり、それらクラスから得られる統計的モデルは正確なものとなる。こうした処理は、新たな稼動条件が与えられ、新たな制御情報の生成が要求されるまでに行なっておけばよく、複雑な挙動を示すプラントについても、稼動時に、従来よりも正確な制御情報を生成することが可能になる。その結果、サンプル間の相関関係を用いて、スペクトラル・クラスタリング手法によりクラスタリングしたサンプルの集合を用いて、精度よくプラントの制御を行なうことができるプラント制御情報生成装置を提供することができる。

好ましくは、行列算出手段は、プラントの稼動条件とデータ項目とにより定義された実測ベクトルのベクトル空間において、複数個の実測サンプルにより定められる複数個の実測ベクトルのうち、基準となる実測サンプルに対応する実測ベクトルを他の実測ベクトルから減算するための減算手段と、減算手段により減算された後の他の実測ベクトルの任意の２個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関を求めるための相関算出手段と、減算手段による減算と、相関算出手段による相関の算出とを、基準となる実測サンプルを変えながら複数個の実測サンプルの全てに対して行なうための繰返算出手段と、繰返算出手段により算出された、複数個の実測サンプルの任意の２個の間で算出される複数個の相関に基づいて、行列を決定するための手段とを含む。

このように、実測サンプルの各々を基準ベクトルとして他の実測ベクトルから減算することにより、そのベクトルは実測ベクトル空間の原点に移動し、それに伴って他の実測ベクトルも平行移動する。その結果、基準実測ベクトルを通る平面は実測ベクトル空間の線形部分空間となり、その上に存在する実測ベクトルが互いにある相関関係を持つと判定できる。こうした判定を全ての実測ベクトルに対して行ない、さらに基準ベクトルとして全ての実測ベクトルを用いてこの処理を繰返すことで、実測ベクトルの間の相関関係がより明確となる。その結果、そうした相関に基づいて算出される類似度行列が、実測ベクトルの間の相関関係の相違を明確に表すものとなって、スペクトル・クラスタリングが正確に行なえる。

より好ましくは、行列を決定するための手段は、繰返算出手段により、複数個の実測サンプルの任意の２個の間で算出された複数個の相関を合計した値を成分とする行列を行列として出力するための手段を含む。

行列の各成分について、単純に複数個の相関を加算することで、複数個の実測サンプルの任意の２個の間の相関の程度を行列として表現できる。必要な計算量が少なくて済み、処理を早めてリアルタイムで実行することができる。

相関算出手段は、減算手段により減算された後の他の実測ベクトルの任意の２個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関係数を求めるための相関係数算出手段と、相関係数算出手段により算出された相関係数が所定のしきい値より大きいか否かにしたがって、相関の値を２値化するための手段とを含んでもよい。

２値化することによって、加算処理がより単純なものとなり、計算量をより削減することができる。

本発明のある実施の形態では、制御情報算出手段は、新たな稼動条件が与えられたことに応答して、当該稼動条件が複数個のクラスのいずれに属するかを判定するための手段と、当該クラスに対応する統計的モデルを、新たな稼動条件に適用することにより、新たな稼動条件の下でのプラントにより生成される物質の属性を予測し、予測情報を出力するための手段とを含む。

本発明の別の実施の形態では、プラント制御情報生成装置はプラントの稼動状態を測定するための測定手段をさらに含む。制御情報算出手段は、測定手段による測定結果が新たな稼動条件として与えられたことに応答して、当該稼動条件と複数個のクラスとを比較することにより、プラントの状態が異常か否かを判定するための異常判定手段を含む。

好ましくは、異常判定手段は、測定手段による測定結果が新たな稼動条件として与えられたことに応答して、複数個のクラスの統計的モデルのうち、当該稼動条件が適合するものを判定するための適合モデル判定手段と、適合モデル判定手段により、新たな稼動条件が複数個のクラスの統計的モデルのいずれにも適合しないと判定されたことに応答して、プラントの異常を示す信号を生成し出力するための手段とを含む。

本発明の第２の局面に係るコンピュータプログラムは、コンピュータにより実行されると、当該コンピュータを、プラントの稼動条件と、当該稼動条件でのプラントの過去の稼動結果から得られるデータ項目の実測値との組合せからなる実測サンプルを複数個記憶するためのデータベースと、データベースに記憶された複数個の実測サンプルを、複数個のクラスにクラスタリングするためのクラスタリング手段と、クラスタリング手段によりクラスタリングされた各クラスに対して、プラントにおける稼動条件とデータ項目の値との関係を示す統計的モデルを構築するためのモデル構築手段と、新たな稼動条件が与えられたことに応答して、統計的モデルと、新たな稼動条件とに基づいて、プラントに関する制御情報を算出するための制御情報算出手段として機能させる。このプログラムは、クラスタリング手段として、コンピュータを、複数個の実測サンプルの任意の２個の間の相関の程度を成分とする行列を、データベースに記憶された実測サンプルから算出するための行列算出手段と、行列算出手段により算出された行列を類似度行列として、複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、複数個の実測サンプルを複数個のクラスにクラスタリングするための手段として機能させる。

本発明の第３の局面に係るプラント制御情報生成方法は、プラントの稼動条件と、当該稼動条件でのプラントの過去の稼動結果から得られるデータ項目の実測値との組合せからなる実測サンプルを複数個記憶するためのデータベースと、データベースに記憶された複数個の実測サンプルを、複数個のクラスにクラスタリングするためのクラスタリング手段と、クラスタリング手段によりクラスタリングされた各クラスに対して、プラントにおける稼動条件とデータ項目の値との関係を示す統計的モデルを構築するためのモデル構築手段と、新たな稼動条件が与えられたことに応答して、統計的モデルと、新たな稼動条件とに基づいて、プラントに関する制御情報を算出するための制御情報算出手段とを含む、プラント制御情報生成装置におけるプラント制御情報生成方法である。クラスタリング手段は、複数個の実測サンプルの任意の２個の間の相関の程度を成分とする行列を、データベースに記憶された実測サンプルから算出するための行列算出手段と、行列算出手段により算出された行列を類似度行列として、複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、複数個の実測サンプルを複数個のクラスにクラスタリングするための手段とを含む。このプラント制御情報生成方法は、クラスタリング手段が、データベースに記憶された複数個の実測サンプルを、複数個のクラスにクラスタリングするステップと、モデル構築手段が、クラスタリングするステップにおいてクラスタリングされた各クラスに対して、プラントにおける稼動条件とデータ項目の値との関係を示す統計的モデルを構築するステップと、制御情報算出手段が新たな稼動条件が与えられたことに応答して、統計的モデルと、新たな稼動条件とに基づいて、プラントに関する制御情報を算出するステップとを含む。クラスタリングするステップは、行列算出手段が、複数個の実測サンプルの任意の２個の間の相関の程度を成分とする行列を、データベースに記憶された実測サンプルから算出するステップと、クラスタリングするための手段が、算出するステップにおいて算出された行列を類似度行列として、複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、複数個の実測サンプルを複数個のクラスにクラスタリングするステップとを含む。

図１は、相関関係の違いに基づき実測サンプルを真のクラスにクラスタリングした結果を示すグラフである。 図２は、図１に示す実測サンプルを対象にｋ−平均法によってクラスタリングした結果を示すグラフである。 図３は、図１に示す実測サンプルを対象に従来のスペクトル・クラスタリング法によってクラスタリングした結果を示すグラフである。 図４は、８個のノードを持つ無向グラフの例を示す図である。 図５は、アフィン部分空間Ｐ内におけるサンプル間の相関関係を説明するための模式図である。 図６は、図５に示す変数群を、サンプルｘ１が原点と一致するように平行移動した後のサンプルの分布を説明するための模式図である。 図７は、図５に示す変数群を、サンプルｘ３が原点と一致するように平行移動した後のサンプルの分布を説明するための模式図である。 変数間の相関関係を指標として、スペクトラル・クラスタリングのための類似度行列を求める、本発明の位置実施の形態に係る装置で採用したアルゴリズムを説明するためのフローチャートである。 図９は、図１に示す測定サンプルを対象に本願発明の一実施の形態によってクラスタリングした結果を示すグラフである。 図１０は、本発明の一実施の形態に係るバッチプロセスの概略構成図を示す模式図である。 図１１は、図１０に示すバッチプロセスにおいて、製品の分子量を推定する、本願発明の一実施の形態に係るソフトセンサを実現するためのコンピュータのブロック図である。 図１２は、本願発明の一実施の形態において、実測サンプルをクラスタリングするためのプログラムの制御構造を示すフローチャートである。 図１３は、本願発明の一実施の形態において、新たなサンプルにしたがって、バッチプロセスで産生される製品の分子量を推定するプログラムの制御構造を示すフローチャートである。 図１４は、バッチプロセスにおける反応炉温度とジャケット温度との時系列変化を示すグラフである。 図１５は、ｋ−平均法によるクラスタリング結果を示すグラフである。 図１６は、本発明の一実施の形態に係るシステムにおけるクラスタリング結果を示すグラフである。 図１７は、ｋ−平均法及び本発明の第１の実施の形態に係るソフトセンサの推定性能を示すグラフでる。 図１８は、本願発明の第２の実施の形態に係るプロセス異常の検出装置における準備処理を実現するプログラムの制御構造を示すフローチャートである。 図１９は、本願発明の第２の実施の形態に係るプロセス異常の検出装置における、プロセス異常を検出するためのプログラムの制御構造を示すフローチャートである。

１．最初に
まず、スペクトラル・クラスタリングの基礎であるグラフ理論を簡単に紹介し、次にスペクトラル・クラスタリングの代表的なアルゴリズムを述べる。

１．１ グラフ理論
電車での乗換について考える。乗換案内図を見ると、駅間の距離及び配置、並びに路線の形状は、実際の地理とは異なって表現されている。これは、乗換を考える場合、駅と駅とがどのように路線で結ばれているかが問題であって、線路が具体的にどこを通っているか、又はその形状がどのようであるかは本質的ではないためである。つまり、電車の乗換では、駅と駅とのつながり方が重要であって、その経路がどのような形状かは問題ではない。その結果、乗換案内図では、駅は点で、駅間の路線は簡略化された単なる線として描かれることが多い。このように、点のつながり方を、点とそれらを結ぶ線とからなる構造で抽象化したものをグラフと呼ぶ。グラフが持つ様々な性質を探求するのがグラフ理論である。

グラフ理論の始まりは一筆書きの問題にある。オイラーはあるとき、ケーニヒスベルグという町で「町の中心を流れるプレーゲル川に架かっている７つの橋を２度通らずに、全て渡って元の所に帰ってくることができるか。ただし、どこから出発してもよい。」という問題を考えた。これをケーニヒスベルグの問題という。オイラーはこの問題を解決するために、川と橋との位置関係をグラフとして表現し、このグラフが一筆書き可能であれば、ケーニヒスベルクの橋を全て１度ずつ通って戻ってくるルートが存在することになると考えた。そして彼は、このグラフが一筆書きできないことを証明し、ケーニヒスベルクの問題を否定的に解決した。一筆書き可能なグラフをオイラーグラフと呼ぶ。

数学的なグラフの定義を与える。集合Ｖ，Ｅと、Ｅの元に２つのＶの元を対応させる写像ｆ：Ｅ→Ｖ×Ｖがあるとき、Ｇ＝（Ｖ，Ｅ）を有向グラフという。Ｖの元をＧのノード（頂点）、Ｅの元をエッジ（辺）と呼ぶ。一方、Ｐ（Ｖ）をＶのべき集合とし、Ｅの元にＰ（Ｖ）の部分集合を対応させる写像ｇ：Ｅ→Ｐ（Ｖ）があって、ｅ∈Ｅに対しｇ（ｅ）＝｛ｖ_１，ｖ_２｝（ｖ_１，ｖ_２∈Ｖ）であるとき、Ｇ＝（Ｖ，Ｅ）を無向グラフと呼ぶ。ｇ（ｅ）＝｛ｖ_１，ｖ_２｝を満たすｖ_１とｖ_２のペアを「互いに隣接している」といい、ｖ_１〜ｖ_２で表す。グラフによっては、エッジに重みがついている場合がある。これを重み付きグラフという。８個のノードを持つ無向グラフの例を、図４に示す。

図４に示すグラフでは、各ノードに１〜８のインデックスを付してあり、各ノードをインデックスで呼ぶものとする。すなわち、Ｖ＝｛１，２，３，４，５，６．７，８｝、Ｅ＝｛｛１，２｝、｛１，５｝、｛２，５｝{２，３}、{３，４}、{４，５}、{４，６}、｛６，７｝、｛６．８｝、｛７，８｝｝と表すことができる。

さらに、２つのグラフＧ＝（Ｖ，Ｅ）とＧ′＝（Ｖ′，Ｅ′）とに対して、全単射ｈ：Ｖ→Ｖ′でｖ_１〜ｖ_２←→ｈ（ｖ_１）〜ｈ（ｖ_２）を満たすものが存在するとき、Ｇ、Ｇ′は同型であるという。そのような写像ｈを同型写像と呼ぶ。同型写像の例としては、ノードのインデックスのみを入れ替える写像がある。

グラフの表現には行列表現が便利である。特に、グラフは「隣接行列」と呼ばれる行列で表現することが多い。無向グラフＧ＝（Ｖ，Ｅ）のノードに重複のないインデックスが与えられているとき、次の式によってこのグラフＧの隣接行列Ａを定義する。

この定義から明らかなように、無向グラフの場合には隣接行列は対称行列である。有向グラフの場合、ノードｖ_１からｖ_２に向かうエッジが存在するときに行列の（Ａ）_{ｖ１，ｖ２}＝１、そうでないときは（Ａ）_{ｖ１，ｖ２}＝０とすればよい。また、重み付きグラフの場合は（Ａ）_{ｖ１，ｖ２}にノードｖ_１からｖ_２に向かうエッジに与えられた重みを代入する。なお、隣接行列の対角成分には０を割り当てることが多い。たとえば、図４で表現される無向グラフの隣接行列Ａ∈Ｒ^８×８は、次式のように表現できる。

無向グラフの隣接行列は対称行列であるからその固有値は全て実数であり、固有ベクトルは互いに直交する。無向グラフＧの隣接行列の固有値をλ_１＜λ_２…＜λ_Ｎとし、それらの重複度をｍ_１，…，ｍ_Ｎとして得られる以下の表をグラフＧのスペクトルと呼ぶ。

同型なグラフのスペクトルは一致するという性質があり（ただし、逆は成り立たない）、スペクトルはグラフの性質の分類に使われる。

グラフで表現できる対象は、乗換案内図及び一筆書きだけではなく数多く存在する。たとえば、電気回路及びＷｅｂのハイパーリンクの構造はグラフで表現できる。計算機のデータ構造は木構造（ツリー）で表現されるが、ツリーもグラフの一種である。巡回セールスマン問題、及び四色問題に代表される彩色問題等、グラフ理論では組合せ最適化問題として定式化される問題が多く、アルゴリズム論の分野で活発に研究されている。

１．２ Ｍａｘ−Ｍｉｎｃｕｔ法
スペクトラル・クラスタリングは、重み付きグラフのグラフ分割法であり、いくつかのエッジをカットしてグラフをサブグラフに分割することで、ノードのクラスタリングを行なう。たとえば図４のグラフの場合、ｖ_４−Ｖ_６間のエッジをカットすることで、Ｇ_１＝｛ｖ_１，…，ｖ_５｝、Ｇ_２＝｛ｖ_６，ｖ_７，ｖ_８｝の２つの部分グラフに分割でき、Ｖ＝｛ｖ_１，…，ｖ_８｝がＶ_１＝｛ｖ_１，…，ｖ_５｝とＶ_２＝｛ｖ_６，ｖ_７，ｖ_８｝とにクラスタリングされる。スペクトラル・クラスタリングにはいくつかのアルゴリズムがあるが、ここではグラフを２つのサブグラフに分割する方法であるＭａｘ−ＭｉｎＣｕｔ（Ｍｃｕｔ）法について述べる。

いま、重み付きグラフＧ（Ｖ，Ｅ）とその隣接行列Ｗが与えられたとし、グラフＧを２つのサブグラフＡ，Ｂに分割する問題を考える。サブグラフＡ，Ｂ間の類似度ｃｕｔ（Ａ，Ｂ）を次の式（４）により定義する。

ただし

である。すなわち、サブグラフ間の類似度ｃｕｔ（Ａ，Ｂ）とは、Ａ，Ｂ間のエッジの重みの総和である。Ｍｃｕｔ法はｃｕｔ（Ａ，Ｂ）を最小、かつＷ（Ａ）及びＷ（Ｂ）を最大化するようなサブグラフＡ，Ｂを探索する方法である。Ｍｃｕｔ法の目的関数Ｊは次式で定義される。

グラフＧ（Ｖ，Ｅ）のノードのインデックスの与え方は任意であるから、隣接行列Ｗを以下のように定義しても一般性は失わない。

ここで、Ｗ_Ａ、Ｗ_ＢはサブグラフＡ，Ｂ内の隣接行列であり、Ｗ_Ａ，ＢはＡ，Ｂ間の隣接行列である。隣接行列の定義より、Ｗ^Ｔ＝Ｗ_Ｂ，Ａである。さらに、サブクラスへの分割を表現するベクトルｘ＝［１，…，１，０，…，０］、ｙ＝［０，…，０，１，…，１］が次式を満たすようにこれらを定義する。

ただし、Ｄ＝ｄｉａｇ（Ｗｅ）、ｅ＝［１，…，１］^Ｔである。これらの式を用いて、目的関数（６）式を書き換えると

が得られる。また、ｘ^ＴＤｙ＝０、ｘ^ＴＷｘ＞０、ｙ^ＴＷｙ＞０であることがわかる。Ｊの各項はｘ，ｙについて対称であるから、以降、第１項のみについて考える。

目的関数（１１）式の目的変数ｘ，ｙはバイナリ変数であるから、この最小化問題はＮＰ困難である。そこで、（１１）式のバイナリ変数を緩和する。ここでは、ノードのインデックスに対応する指標ベクトルｑを導入する。ｑの、ノードｕに対応する要素をｑ_ｕ＝｛ａ，−ｂ｝（ａ，ｂ＞０）とする。ｑ_ｕ＝ａのときｕ∈Ａ、ｑ_ｕ＝−ｂのときｕ∈Ｂである。目的関数は次の（１２）式のように書き直される。

さらに（１２）式を次のようにスケーリングする。

ここで、^〜Ｗ＝Ｄ^−１／２ＷＤ^−１／２、^〜ｑ＝Ｄ^１／２ｑ／｜Ｄ｜^１／２である（各変数の前の「^〜」は、式中では変数の上に記載されている。）。また、^〜ｑ^Ｔ〜Ｗ^〜ｑ＞０である。^〜Ｗの要素（^〜Ｗ）_ｉ，ｊは１≧（^〜Ｗ）_ｉ，ｊ≧０、各列の和は１であるから、〜Ｗは確率推移行列とみなすことができる。

さらに目的関数^〜Ｊ_ｑにおいて、Ｐ≡Ｉ−^〜Ｗとおくと^〜Ｊ_ｑはＰのＲａｙｌｅｉｇｈ商で表現できることがわかる。Ｒａｙｌｅｉｇｈの原理によると、行列Ａの商Ｒ（ｘ）は最小固有値λ_１に対応する固有ベクトルｘ_１によって最小となり、その最小値はλ_１である（本明細書では、固有値はλ_１≦…≦λ_Ｎのように並べるものとする。）。よって、緩和した最小化問題はＰの固有値問題に帰着することがわかる。Ｐのある固有ベクトルをＤ^１／２ｅと仮定する。Ｗの第ｉ行をｗ_ｉとすると、（Ｄ）_ｉ，ｉ＝ｄ_ｉ＝ｗ_ｉｅであるから

であり、確かにＰのある固有ベクトルはＤ^１／２ｅで、その固有値は０であることが示される。ところで、^〜Ｗは確率推移行列とみなすことができたから、その最大固有値は１である。よって、ｘ^ＴＰｘ＝ｘ^Ｔ（Ｉ−^〜Ｗ）ｘ≧０であるから、Ｐは半正定値行列である。これより、固有ベクトルＤ^１／２ｅはＰの最小固有値λ＝０に対応する固有ベクトルｚ_１であることがわかる。

したがって、Ｒａｙｌｅｉｇｈの原理より最小固有値に対応する固有ベクトルｚ_１が（１３）式を最小とする解である。しかし、この解は全ての要素が正であり、求めるｑの要素がｑ_ｕ＝｛ａ，−ｂ｝であるという仮定に反する。したがってこの解は採用できない。次に小さな解はＲａｙｌｅｉｇｈ商のｍａｘ−ｍｉｎ定理より、２番目に小さな固有値λ_２である。λ_２に対応する固有ベクトルｚ_２は、Ｐの対称性よりｚ_１ ^Ｔｚ_２＝０であるから、必ず正負の要素を持つことがわかる。すなわち、この解は最初の仮定に反しない。したがって、λ_２に対応する固有ベクトルｚ_２が求める最適な^〜ｑ^＊である。

以上の考察より、最終的な定式化は
において、固有値問題
の２番目の大きさの固有値λ_２に対応する固有ベクトルｚ_２を選択すればよいことになり、求めるｑ^＊はｑ^＊＝Ｄ^−１／２ｚ_２で表される。なお、このλ₂、ｚ_２をそれぞれＦｉｅｌｄｅｒ値、Ｆｉｅｌｄｅｒベクトルと呼ぶ。

Ｍｃｕｔ法は任意のグラフを２つのサブグラフに分割する方法である。したがって、グラフを任意の個数のサブグラフに分割したいときは、分割後のサブグラフに対して上記手順を繰返せばよい。

理論的な考察により、緩和問題の目的関数が最小となるとき、（６）式においてほぼＷ（Ａ）＝Ｗ（Ｂ）が成立することが示されている。すなわち、Ｍｃｕｔ法ではサブグラフ内の類似度がバランスすると期待できる。他のスペクトラル・クラスタリングのアルゴリズムでは必ずしもサブグラフ内の類似度がバランスしないため、このことがＭｃｕｔ法の利点とされる。

１．３ その他のアルゴリズム
前節で紹介したＭｃｕｔ法に類似の手法として、Ｎｏｒｍａｌｉｚｅｄ Ｃｕｔ（Ｎｃｕｔ）法がある。Ｎｃｕｔ法の評価関数は次の（１７）式で表される。

ここで、Ｖ＝Ａ∪Ｂ、ｄｅｇ（Ａ，Ｖ）＝Σ_{ｕ∈Ａ，ｖ∈Ｖ}Ｗ_ｕ，ｖであり、ｄｅｇ（Ｂ，Ｖ）も同様に定義される。解析の詳細はＭｃｕｔ法とほぼ同等である。

Ｍｃｕｔ/Ｎｃｕｔ法では類似度の与え方は任意であり、類似度行列が無向グラフの隣接行列としての性質を満たしていればよい。

一方、類似度行列の構成にＧａｕｓｓｉａｎカーネルを用いることが従来の一手法として提案されている。この手法では、サンプルｓ_ｉ，ｓ_ｊの類似度（Ｗ）_ｉ，ｊとして

を用いる。ここでｄ（・，・）は距離関数、σは調整パラメータである。この手法に関連してさらに、グラフの２分割を繰返すＭｃｕｔ法と異なり、スケーリング後の類似度行列の固有空間にてサンプルの集合を直接任意の個数のクラスタに分割する手法が提案されている。この場合、カーネルを用いて非線形写像を行なった後に固有値問題を解いていることになるから、カーネル主成分分析（主成分分析を以下「ＰＣＡ」と呼ぶ。）と類似の手法であるとされる。この手法をＮｇ−Ｊｏｒｄａｎ−Ｗｅｉｓｓ（ＮＪＷ）アルゴリズムと呼ぶ。図３に示したのはＮＪＷアルゴリズムを用いたスペクトル・クラスタリングの結果である。

ＮＪＷアルゴリズムの（１８）式のσの決め方において、ローカルスケーリングという方法が提案されている。（Ｗ）_ｉ，ｊは次のように書き換えられる。

ここで、σ_ｉ＝ｄ（Ｓ_ｉ，Ｓ_Ｋ）であり、Ｓ_ＫはＳ_ｉのＫ番目に近傍のサンプルを示す。Ｋは調整パラメータであるが、数値実験ではＫ＝７で概ね良好な結果が得られたと報告されている。

１．４ アプリケーション
既に述べたように、グラフとして表現できる構造は多いため、スペクトラル・クラスタリングも応用は広い。たとえば、Ｗｅｂのハイパーリンク及びソーシャルネットワークの構造解析、ＣＡＤ等のメッシュ分割、タンパク質ネットワークの解析、並びにコンピュータビジョン等にスペクトル・クラスタリングが用いられている。各問題に適した類似度の与え方に工夫がある。

２ 変数間の相関関係を指標とした類似度行列
前節で紹介したＭｃｕｔ法では、類似度の定義は任意であった。一方、ＮＪＷアルゴリズムでは、類似度として、サンプル間の距離をカーネルによって写像した値を用いている。しかし、いずれにせよ従来法では変数間の相関関係を指標としてクラスタリングを行なうことは考慮されていない。そこで本実施の形態では、スペクトラル・クラスタリングにおいて変数間の相関関係を考慮した類似度行列の構成法を実現する。本実施の形態に係る方法は、変数間の相関関係を考慮した類似度行列を構築し、スペクトラル・クラスタリングによってクラスタリングを実施する方法である。この方法を本明細書ではＮＣ−スペクトラル・クラスタリングと呼ぶことにする。

本実施の形態で採用したクラスタリング方法のコンセプトを述べる。なお、以下では、説明を理解しやすくするために、サンプルが３次元空間内に存在するものとする。図５を参照して、この３次元空間内のアフィン部分空間Ｐが変数間の相関関係を表しており、Ｐ上のサンプルは全て同一の相関関係にしたがっているものとする。すなわち、Ｐ上のサンプルｘ_１,…,ｘ_５は同一の相関関係にしたがうサンプルである。一方、Ｐ上には存在しないサンプルｘ_６及びｘ_７は異なった相関関係を有している。

いま、サンプルｘ_１に着目する。まず、サンプルｘ_１が原点となるように空間全体を平行移動させる。つまり、全サンプルｘ_ｉ（ｉ＝１，…，７）からｘ_１を減算する。この操作によってアフィン部分空間Ｐは原点を含むことになるため線形部分空間の定義を満たす。この線形部分空間をＶとする。次に図６に示すように、任意のサンプルと原点を結ぶ直線を引き、この直線上で別のサンプルが発見できたとする。図６に示す例では、サンプルｘ_２及びｘ_５並びにサンプルｘ_３及びｘ_４がこのような関係を満たしている。このとき、これらのサンプルのペアの相関係数の絶対値は１である。このような関係が成り立つとき、「ｘ_２〜ｘ_５」のように書く。一方、部分空間Ｖの要素ではないサンプルｘ_６及びｘ_７の相関係数の絶対値は１未満である。それゆえ、相関係数が±１であるペアのサンプルは、同一の相関関係を有していると判定できる。

実際は相関係数が厳密に±１になるペアは存在しないため、しきい値γ（０＜γ≦１）を用いて相関関係を２値化して、同一の相関関係を有しているかを判定する。すなわち、相関係数の絶対値がしきい値γ以上となるペアを、同一の相関関係を有していると判定する。本実施の形態では、サンプルの幾何学的な配置のみを考慮したここまでの手続きをＮＣ法と呼ぶ。ＮＣ法自体は前述したとおり公知である。

ここで、各サンプルをノードとし、サンプル間が同一の相関関係を有しているか否かをそれらのサンプルに対応するノードが隣接しているかとして表現した無向グラフ及びその隣接行列Ｓ_１∈Ｒ^７×７を考える。この例では、ｘ_２〜ｘ_５及びｘ_３〜ｘ_４であるから、隣接行列Ｓ_１は、（Ｓ１）_２，５＝（Ｓ１）_５，２＝１及び（Ｓ１）_３，４＝（Ｓ１）_４，３＝１とし、それ以外の要素を０とする行列として表現できる。

次に、サンプルｘ_２が空間の原点と一致するように平面を平行移動させて同様の手順を行ない、隣接行列Ｓ_２を求める。これを全てのサンプルについて繰返し、最終的にＳ_１，…，Ｓ_７を求める。

この手続きによって隣接行列が７つ得られる。これら７つのグラフは、同一のノードの集合を持ち、エッジの集合の部分集合が異なる７つの無向グラフと考えることができる。すなわち、Ｇ_１（Ｖ，Ｅ（Ｇ_１））,…,Ｇ_７（Ｖ，Ｅ（Ｇ_７））である。ここで、これらのグラフの和のグラフＧ′＝Ｇ_１∪…∪Ｇ_７を考えると、グラフＧ′のエッジの集合はグラフＧ_１からＧ_７のエッジの和集合Ｅ（Ｇ_１）∪…∪Ｅ（Ｇ_７）であることがわかる。これは隣接行列Ｓ_ｉ（ｉ＝１，…，７）のブール和（各行列の対応する要素の論理和を要素とする行列）を求めることに等しいから、グラフＧ′の隣接行列はＳ′＝Ｓ_１＋…＋Ｓ_７となる。ただし、ここで＋は行列のブール和を表す。以上の手続きより、サンプル間の相関関係を表現したグラフを得ることができる。

しかし、サンプルの分布と原点の選び方とによっては、異なった相関関係を持つサンプルのペアであっても、偶々同一の相関関係を有すると判断されてしまうことがある。たとえば図５及び図６に示した例について考える。図７のようにサンプルｘ_３が原点と一致するように平面を平行移動させると、ｘ_１〜ｘ_３以外にｘ_６〜ｘ_７と判定されてしまう（サンプルｘ_５と原点とを結んだ直線上にサンプルｘ_７が存在する。）。つまり、グラフＧ′では真に同一の相関関係を有するサンプル間のエッジと、偶々同一の相関関係を有すると判断されてしまったサンプル間のエッジとが区別できない。

そこで、上記手続きにおいて、それぞれのサンプルのペアで相関関係を有すると判定された回数をカウントし、それに応じてエッジに重み付けを行なう。このグラフをＧとする。つまり、隣接行列のブール和をとるのではなく、通常の和を用いてサンプル間の相関関係の大きさを決定し、Ｇの隣接行列Ｓを求める。

隣接行列Ｓ_ｉは対称行列であったからその和も対称行列である。したがって、Ｓを類似度行列としてスペクトラル・クラスタリングを行なうことで、変数間の相関関係を指標としたクラスタリングを実施できることになる。

図８を参照して、本実施の形態に係るＮＣ−スペクトル・クラスタリングを実現するコンピュータプログラムの制御構造について説明する。いま、サンプルｘ_ｎ∈Ｒ^Ｍ（ｎ＝１，…，Ｎ）がデータベースに保存されており、これより変数間の相関関係を指標として、スペクトラル・クラスタリングのための類似度行列を求めたいとする。

図８を参照して、このプログラムは、以下の繰返し処理のため繰返し制御変数Ｌに１を代入するステップ３０と、類似度行列Ｓのための２次元配列Ｓ［Ｎ，Ｎ］の領域を確保するステップ３２と、配列Ｓ［Ｎ，Ｎ］の全ての要素をゼロでクリアするステップ３４とを含む。

このプログラムはさらに、Ｌ番目（１≦Ｌ≦Ｎ）のサンプルについて、全サンプルとの間の隣接行列Ｓ_Ｌを保存するための２次元配列Ｓ_Ｌ［Ｎ，Ｎ］の領域を確保するステップ３６と、配列Ｓ_Ｌ［Ｎ，Ｎ］の各要素をゼロでクリアするステップ３８と、Ｌ≦ｎ≦Ｎを満たす全てのｎについてサンプルｘ_ｎからサンプルｘ_Ｌを減算して新たなサンプルｘ^′ _ｎを算出するステップ４２と、ステップ４２で算出されたサンプルｘ^′ _ｎに関し、全てのｘ^′ _ｋとｘ^’ _ｌとの組合せ（１≦ｋ、ｌ≦Ｎ、ただしｋ≠ｌ）について、相関係数Ｃ_ｋ，ｌを算出するステップ４６と、ステップ４８で算出された相関係数Ｃ_ｋ，ｌのうちで｜Ｃ_ｋ，ｌ｜≧γ（γは予め定めるしきい値）を満たすｋ及びｌを求めるステップ４８と、Ｓ_Ｌ［ｋ、ｌ］とＳ_Ｌ［ｌ、ｋ］とに１を代入するステップ５０と、配列Ｓの各要素に、配列Ｓ_Ｌの対応する要素の値を加算するステップ５２とを含む。

このプログラムはさらに、変数Ｌの値がＮ−１より小さいか否かを判定し、判定結果に応じて制御の流れを分岐させるステップ５４と、ステップ５４における判定結果がＹＥＳのときに実行され、変数Ｌに１を加算して制御をステップ３６に戻すステップ５６と、ステップ５４における判定結果がＮＯのときに実行され、類似度行列を表す配列として配列Ｓを出力し、処理を終了するステップ５８とを含む。

このプログラムは類似度行列を算出するためのものであり、スペクトラル・クラスタリング自体のアルゴリズムは任意である。本実施の形態ではＮＪＷアルゴリズムを用いる。

３ クラスタリングのシミュレーション例
シミュレーションにより、従来法と本実施の形態に係るスペクトル・クラスタリングによるクラスタリングの識別性能の比較を行なった。本シミュレーション例では、２次元及び５次元のデータをそれぞれ用いて２通りの比較を行なった。

３．１ ２次元の場合
それぞれ互いに異なる相関関係を有する３クラスのデータを用い、クラスタリングを実施する。各クラスに属するサンプルは次式を用いて生成した。

ここでＳ〜Ｎ（０，１０）であり、ａ_ｉ∈Ｒ^２は係数行列である。また、ｎ＝［ｎ_１ ｎ_２］^Ｔ、ｎ_ｉ〜Ｎ（０，０．１）である。ただし、Ｎ（ｍ，σ）は平均ｍ、標準偏差σの正規分布にしたがう乱数を表す。係数ベクトルはそれぞれ、ａ_１＝［１ ２］^Ｔ、ａ_２＝［２ ２］^Ｔ、ａ_３＝［２ １］^Ｔとし、各クラスにおいて１００サンプル生成した。

本シミュレーション例において、本実施の形態のステップ４８で用いるパラメータ（しきい値）γは、０．９９９とした。生成したサンプル及びｋ−平均法、従来のスペクトラル・クラスタリングでのクラスタリング結果が、それぞれ図１〜図３に示したものである。本願実施の形態にしたがったクラスタリング結果を図９に示す。

図１〜図３と、図９とを参照して、ｋ−平均法でクラスタリングした結果（図２）及び従来のスペクトル・クラスタリングでクラスタリングした結果（図３）は、いずれも図１と比較すると正しいクラスタリングとはいえないことがわかる。一方、本実施の形態によるクラスタリング結果（図９）では、原点付近のサンプルの属するクラスは正しく識別できていないが、それ以外の領域では正しく識別できていることがわかる。

３．２ ５次元の場合
前節と同様にそれぞれ互いに異なる相関関係を有する３クラスのデータを用い、クラスタリングを実施する。各クラスに属するサンプルは次式を用いて生成した。識別率は次式で定義される。

ここでＫは識別すべきクラスのサンプルの数であり、ＬはＫ個のサンプルのうち実際に対象のクラスと同一であったサンプルの数である。各クラスに属するサンプルは次式を用いて生成した。

ここでＡ_ｉ∈Ｒ^２×５は係数行列であり、Ｓ＝［ｓ_１ ｓ_２］、Ｓ_ｉ〜Ｎ（０，１０）である。また、ｎ＝［ｎ_１ ｎ_２ ｎ₃ ｎ₄ ｎ₅］^Ｔ、ｎ_ｉ〜Ｎ（０，０．１）である。係数行列Ａ_ｉは

とし、各クラスについて１００サンプル生成した。すなわち、Ｋ＝１００である。本シミュレーション例において、本実施の形態のステップ４８で用いるパラメータ（しきい値）γは、０．９９９とした。

生成したサンプルに対し、ｋ−平均法（ＫＭ）、スペクトラル・クラスタリング（ＳＣ）、及び本実施の形態に係るＮＣ−スペクトラル・クラスタリング（ＮＣ）のクラスタリングにより各クラスについてどの程度の精度で識別できたかをテーブル１に示す。テーブル１を参照して、ＮＣ−スペクトル・クラスタリングによれば、ｋ−平均法及び従来のスペクトル・クラスタリングのいずれと比較しても、どのクラスについてもより高い精度でクラスを識別できていることが分かる。以上の結果より、本実施の形態に係るＮＣ−スペクトラル・クラスタリングの有効性が示されたといえる。

４．第１の実施の形態
４．１ 概略構成
本実施の形態は、並列化されたバッチプロセスの運転データのクラスタリングを対象とし、クラスタリング結果に基づいてソフトセンサ設計を行なうものである。この実施の形態のシミュレーションで用いたバッチプロセスモデルの詳細を最後に付録として示す。

図１０を参照して、本実施の形態に係るプラント８０は、各々、内部で逐次並列反応

（ただし、Ａ，Ｂは原料であり、Ｃは製品、Ｄは副製品）を行なう５台の反応炉８２，８４，…，８６と、これら反応炉、及びこれら反応路のジャケットの温度をそれぞれ測定するための図示しない複数個のセンサと、これら複数個のセンサの出力にしたがって、所定のプロセス制御を反応炉ごとに行なう複数のコントローラ９２，９４，…，９６と、これら反応炉で過去に行なわれたバッチプロセスの結果から得られた生成物に関する情報、及び各プロセスを開始するときの各反応炉に関する所定のプロセスパラメータ（ここではジャケット温度Ｔｊと反応炉内温度Ｔｒとする。）に基づき、新たなバッチプロセスで得られる製品Ｃの分子量Ｍ_Ｃを推定するための、コンピュータからなる製品分子量推定装置９８とを含む。本プラントでは、目的生成物は製品Ｃである。プロセスの目的は副製品Ｄの生成を抑制し、できるだけ多くの製品Ｃを得ることである。バッチプロセス運転開始前に原料Ａ及びＢが室温で反応炉に仕込まれる。

反応１と反応２との進行は、いずれも温度によって変化するものとする。条件によっては、反応１により生成した製品Ｃが、反応２によって消費されて副製品Ｄが精製してしまう。したがって、反応炉内の温度を調整する必要がある。プラント８０はその目的のため、冷媒と熱媒とを用いたジャケット温度調節部１００を有している。

製品分子量推定装置９８は、実質的にコンピュータからなり、反応炉８２，８４，…，８６による過去のバッチプロセスから得られたデータ（反応炉温度とジャケット温度、及び最終的に得られた製品Ｃの分子量Ｍ_Ｃ）を記憶し、これらをそれぞれ適切なクラスに分類し、各クラス別に、各バッチの開始時のデータと、バッチプロセスにより得られた製品Ｃの分子量とに関する統計的モデルをパラメータの値から構築する機能を持つ。製品分子量推定装置９８はさらに、反応炉温度とジャケット温度とからなる新たなデータが与えられると、それを適切なクラスに分類し、そのクラスの統計的モデルを用いて製品Ｃの分子量を推定する。

４．２ コンピュータによる実現
上述の実施の形態は、コンピュータシステムハードウェア及びコンピュータシステム上で実行されるプログラムによって実現され得る。図１１はこの実施の形態で用いられる製品分子量推定装置９８のブロック図である。なお、ここで示す製品分子量推定装置９８は単なる例であって、他の構成も利用可能である。

図１１を参照して、製品分子量推定装置９８は、コンピュータ１２０と、全てコンピュータ１２０に接続された、モニタ１２２と、キーボード１２６と、マウス１２８と、プリンタ１２４と、を含む。さらに、コンピュータ１２０はＤＶＤ−ＲＯＭ（Ｄｉｇｉｔａｌ Ｖｅｒｓａｔｉｌｅ Ｄｉｓｋ Ｒｅａｄ−Ｏｎｌｙ−Ｍｅｍｏｒｙ：ディジタル多用途ディスク読出専用メモリ）ドライブ１５０と、半導体メモリポート１５２とを含む。

図１１を参照して、コンピュータ１２０はさらに、ＤＶＤ−ＲＯＭドライブ１５０と半導体メモリポート１５２とに接続されたバス１４２と、全てバス１４２に接続された、ＣＰＵ１４０と、コンピュータ１２０のブートアッププログラムを記憶するＲＯＭ１４４と、ＣＰＵ１４０によって使用される作業領域を提供するとともにＣＰＵ１４０によって実行されるプログラムのための記憶領域となるＲＡＭ１４６と、音声データ、音響モデル、言語モデル、レキシコン、及びマッピングテーブルを記憶するためのハードディスクドライブ１４８と、ネットワーク１６４への接続を提供するネットワークインターフェイス１５４とを含む。

上述の実施の形態に係る製品分子量推定装置９８を実現するソフトウェアは、ＤＶＤ−ＲＯＭ１６２又は半導体メモリ１６０等の媒体に記録されたオブジェクトコード又はスクリプトの形で流通し、ＤＶＤ−ＲＯＭドライブ１５０又は半導体メモリポート１５２等の読出装置を介してコンピュータ１２０に提供され、ハードディスクドライブ１４８に記憶される。ＣＰＵ１４０がプログラムを実行する際には、プログラムはハードディスクドライブ１４８から読出されてＲＡＭ１４６にロードされる。図示しないプログラムカウンタによって指定されたアドレスから命令がフェッチされ、その命令が実行される。ＣＰＵ１４０はハードディスクドライブ１４８から処理すべきデータを読出し、処理の結果をこれもまたハードディスクドライブ１４８に記憶する。推定された製品分子量のハードコピーが、プリンタ１２４により出力される。

コンピュータ１２０の一般的動作は周知であるので、詳細な説明は省略する。

ソフトウェアの流通の方法に関して、ソフトウェアは必ずしも記憶媒体上に固定されたものでなくてもよい。たとえば、ソフトウェアはネットワークに接続された別のコンピュータから分配されてもよい。ソフトウェアの一部がハードディスクドライブ１４８に記憶され、ソフトウェアの残りの部分についてはネットワークを介してハードディスクドライブ１４８に取込み、実行の際に統合する様にしてもよい。

また、ソフトウェアの流通形態はオブジェクトコードには限らない。前述したようにスクリプト形式でもよいし、ソースプログラムの形で供給され、コンピュータ１２０にインストールされた適切なコンパイラでオブジェクトコードに変換されるという流通形態もあり得る。

典型的には、現代のコンピュータはコンピュータのオペレーティングシステム（ＯＳ）によって提供される一般的な機能を利用し、所望の目的にしたがって制御された態様で機能を達成する。したがって、ＯＳ又はサードパーティから提供されうる一般的な機能を含まず、一般的な機能の実行順序の組合せのみを指定したプログラムであっても、そのプログラムが全体として所望の目的を達成する制御構造を有する限り、そのプログラムがこの発明の範囲に包含されることは明らかである。

図１２及び図１３を参照して、コンピュータ１２０により実行され、製品Ｃの分子量Ｍ_Ｃを推定するためのプログラムの制御構造について説明する。このプログラムは、バッチプロセスが終了すると、そのバッチプロセスによって得られた実測データを過去の実測データとともにデータベースに保存し、将来の製品分子量推定に備えてクラス別の統計的モデルを構築するための統計的モデル構築プログラム（図１２）と、新たなデータが与えられると、構築済のクラス別の統計的モデルのうちで最適な統計的モデルを選択し、その統計的モデルを新たなデータに適用することにより、最終的な製品分子量Ｍ_Ｃを推定するためのプログラム（図１３）とを含む。

図１２を参照して、統計的モデル構築プログラムは、バッチプロセスが終了し、そのバッチプロセスによって得られた実測データ（プロセス開始時の反応炉温度、ジャケット温度その他の関連パラメータ、並びに最終的に得られた製品Ｃの平均分子量）が与えられると、その実測データをコンピュータ１２０（図１１）内に構築されたデータベース（以下単に「ＤＢ」と呼ぶ。）に追加するステップ１８０と、ステップ１８０に続き、ＤＢ内の実測データに対して主成分分析を行なってクラスタリングのための次元数を削減するステップ１８２とを含む。本実施の形態では、図面を理解しやすくするため、主成分分析により、実測データの各々を第１主成分及び第２主成分からなる２次元のデータに変換するものとする。

このプログラムはさらに、ステップ１８２に続き、ステップ１８２により得られた実測データに対して、図８を参照して説明したＮＣ−スペクトラル・クラスタリングのためのプログラムを適用することにより、実測データを複数のクラスに分類するステップ１８４と、ステップ１８４により得られたクラスごとに、予め選択された構築方法により統計的モデルを構築するステップ１８６と、ステップ１８６で得られた各クラスの統計的モデルをＨＤＤ１４８に保存して処理を終了するステップ１８８とを含む。

一方、図１３を参照して、新たなサンプル（新たなバッチロセス開始時の反応炉温度、ジャケット温度その他の関連パラメータ等）が与えられたことに応答して、図１２のステップ１８２での主成分分析の結果にしたがって主成分スコアを算出するステップ２００と、この主成分スコアに基づき、新たなサンプルがどのクラスに属するかを推定するステップ２０２と、ステップ２０２で推定されたクラスに対する統計的モデルをＨＤＤ１４８から読出し、新たなサンプルを当該統計的モデルに適用することで新たなバッチプロセスにより生成する製品Ｃの分子量を推定して処理を終了するステップ２０４とを含む。

ステップ２０２では、たとえば各クラスとの間で求めたＱ統計量を指標としてクラスを推定することが考えられる。

４．３ シミュレーション
本実施の形態について、以下のような条件でその性能についての数値シミュレーションを行ない、従来の方法による結果と比較した。以下の数値シミュレーションでは、バッチプロセス運転開始前に原料Ａ及びＢが室温２０℃で反応炉に仕込まれるものとする。このシミュレーションでは、原料Ａ及びＢの仕込み量はそれぞれＮ（２０，０．１）にしたがう乱数として変化するものとする。

反応１は９０℃以上で進行し、反応炉温度が１００℃を超えると反応２が急速に進行するものとする。反応１を促進するため、バッチ開始後は速やかに反応炉を室温から９０℃まで加熱する必要がある。反応が進むと反応熱によって反応炉内温度がさらに上昇するが，反応炉内温度が１００℃を超えると反応２が急速に進行し，製品Ｃが消費され，副製品Ｄが生成してしまう。したがって、反応炉内温度が９０℃付近まで上昇した後は、ジャケット温度調節部１００により、反応炉を冷却し反応炉内温度の上昇を抑制する。シミュレーションでは、反応炉内温度制御には、デュアルモード制御が用いられるものとした。

本プロセスにおいて運転中に測定されるのは，ジャケット温度Ｔｊと反応炉内温度Ｔｒであり、それぞれサンプリング周期は１分である。バッチ終了時間は１２０分で固定する。あるバッチの運転データを図１４に示す。

本シミュレーションでは５台の反応炉８２，８４，…，８６が並列に運転されている。バッチ終了時にはこれら反応炉のメンテナンスが行なわれるため、バッチ毎に反応炉の伝熱係数が変化する。反応炉それぞれの伝熱係数Ｕ_ｉ（ｉ＝１，…，５）は、Ｕ_１＝Ｕ_２＝Ｕ（４０．６０，４０．６２），Ｕ_３＝Ｕ_４＝Ｕ（４０．５７，４０．５９），Ｕ_５＝Ｕ（４０．５４，４０．５６）として。ただし、Ｕ（ａ，ｂ）は閉区間［ａ，ｂ］上の一様乱数であることを示す。すなわち，本プロセスでは５台の反応炉が存在するが，それらの個体差としては３種である。本シミュレーションでは、各反応炉についてプロセス２０バッチ、計１００バッチのデータがＤＢに保存されているとし、これらのデータより製品の分子量推定を行なった。

４．４ クラスタリング
データベースに保存されている１００バッチ分の運転データのクラスタリングを，従来法であるｋ−平均法と、本実施の形態に係るＮＣ−スペクトラル・クラスタリングを用いて実施した。シミュレーションでは、図１２のステップ１８２から処理を開始した。すなわち、１００データに対してＰＣＡを行なって各データの次元を圧縮した。ここで採用する主成分数は２、主成分スコアｔ_１，ｔ_２を入力としてクラスタリングを実施した。クラス数は３とし，ＮＣ−スペクトラル・クラスタリングにおけるしきい値γ＝０．９９とした。

クラスタリング結果を図１５及び図１６に示す。ここではｔ_１−ｔ_２平面上でのサンプルの分布を示している。ｋ−平均法ではさらに、各クラスの中心を○印で示している。この結果より，ｋ−平均法ではクラス中心からの距離に基づいてサンプルを分類しているが、本実施の形態では各クラスは相関関係の違いを表現していることがわかる。

４．５ ソフトセンサ設計
製品Ｃの分子量Ｍ_Ｃ［ｋｍｏｌ］を推定するソフトセンサを設計する。設計手順は図１２及び図１３に示したとおりであり、図１３のステップ２０４で新たなプロセスでの製品の分子量を推定した。なお、対比のために、ｋ−平均法を用いたクラスタリングも行ない、その結果を用いて図１３に示す方法で製品の分子量を推定した。クラスタリング時のクラス数は３である。

以下のように、３クラスそれぞれで統計的モデルを構築した。ここで、統計的モデルの入力は反応炉温度Ｔ_ｒ（ｋ）及びジャケット温度Ｔ_ｊ（ｋ）（ｋ＝０，…，１２０）、出力は製品Ｃの分子量Ｍ_Ｃである。モデル構築にはＰＬＳ（Ｐａｒｔｉａｌ Ｌｅａｓｔ Ｓｑｕａｒｅ）回帰を用いた。採用した潜在変数の数はともに３である。新たに測定されたサンプルの識別には、ｋ−平均法ではクラス中心からの距離、本実施の形態では前述の通り各クラスとの間で求めたＱ統計量を指標として用いた。

本シミュレーションでは、新たな反応炉Ｒ_６のＭ_Ｃを推定した。Ｒ_６の伝熱係数はＵ_６＝Ｕ（４０．６０，４０．６２）とした。検証には２０サンプルを用い、ｋ−平均法及び本実施の形態それぞれを用いたソフトセンサの推定性能を比較した。

結果を図１７（Ａ）（Ｂ）にそれぞれ示す。図中、グラフの横軸が真値、縦軸が予測値を示している。「ＲＭＳＥ」は根平均２乗誤差、「Ｒ」は真値と予測値との相関係数である。この結果より、ｋ−平均法を用いた場合は精度の良い推定ができていないが、本実施の形態を用いた場合は高い推定性能を達成していることが分かる。たとえば、ｋ−平均法を用いた場合のＲＭＳＥが０．００７であるのと比較して、本実施の形態ではＲＭＳＥは０．００４となり、ＲＭＳＥにして４３％改善した。このことより、本実施の形態によれば、ｋ−平均法を用いた場合よりも装置の個体差を正しくクラスタリングできていることがわかる。

５． 第２の実施の形態
５．１ 構成
第１の実施の形態と同様のシステムハードウェア構成で、プロセス異常を検出する実施の形態について、説明する。

図１８に、この実施の形態に係るプロセス異常の検出装置を実現するためのプログラムの制御構造をフローチャート形式で示す。このプログラムは、図１１に示すＣＰＵ１４０によって実行され、それによって以下に説明するようなプロセス異常の検出装置を実現することができる。

図１８を参照して、本実施の形態では、上記したＮＣ−スペクトラル・クラスタリングと多変量統計的プロセス管理手法（Ｍｕｌｔｉｖａｒｉａｔｅ Ｓｔａｔｉｓｔｉｃａｌ Ｐｒｏｃｅｓｓ Ｃｏｎｔｒｏｌ；ＭＳＰＣ）を組合せた異常検出を行なう。なお、比較のためにＮＣ−スペクトラル・クラスタリングに代えてｋ−平均法を用いた場合の結果についても後に示す。以下、ｋ−平均法又はＮＣ−スペクトル・クラスタリングと、ＭＳＰＣとを組合せた異常検出法をそれぞれＫＭ−ＭＳＰＣ、ＮＣ−ＭＳＰＣと呼ぶことにする。

図１８を参照して、本実施の形態に係るプロセス異常の検出装置を実現するためのプログラムの制御構造は以下のとおりである。すなわち、このプログラムは、これまでにＤＢに記憶されていた全サンプルに対してＮＣ−スペクトル・クラスタリングを用いたクラスタリングを実施するステップ２３０と、ステップ２３０に続き、全クラスについてＴ^２統計量とＱ統計量との管理限界^￣Ｔ^２ _Ｓ、^￣Ｑ_Ｓ（Ｓ＝１，２，…，Ｓ）（各記号の前の「^￣」は、通常は各記号の直上に記載されるべきものである。）とを設定するステップ２３２とを含む。ただし、Ｓはクラス数である。

さらに、本実施の形態はバッチプロセスを対象とし、サンプルはジャケット温度Ｔ_ｊと反応炉内温度Ｔ_ｒである。したがって、得られたデータは時系列データとなる。得られた時系列データを一定のインターバルで分割し、インターバル終了ごとにプロセスに異常が発生しているかを判定する。以下、そうした判定を行なうプログラムの制御構造について説明する。

図１９を参照して、このプログラムは、所定のインターバル（たとえば１０分）が経過するまで待機するステップ２５０と、１０分が経過するごとに、対象となる反応炉の反応炉温度、ジャケット温度等の新たな測定値ｘ_ｑをサンプリングするステップ２５２と、ステップ２５２でサンプリングされた測定値ｘ_ｑと、新たなサンプルｘ_ｑと、全てのクラスとの間でＴ^２統計量とＱ統計量Ｔ_ｑ ^２、Ｑ_ｑとを求めるステップ２５４と、全てのクラスに対して第１の判定式Ｔ_ｑ ^２≧^￣Ｔ^２ _Ｓが成立しているか否かを判定し、判定結果に応じて制御の流れを分岐させるステップ２５６と、ステップ２５６で第１の判定式が成立していないと判定されたときに、全てのクラスに対して第２の判定式Ｑ_ｑ≧^￣Ｑ_Ｓが成立しているか否かを判定し、判定結果に応じて制御の流れを分岐させるステップ２５８とを含む。

ステップ２５６又は２５８における判定結果がＹＥＳのとき、すなわち新たにサンプリングされた測定値ｘ_ｑがいずれのクラスにも属していないと判定されたときには、ステップ２６２でプロセス異常が発生したと判定され、必要な処理が実行される。ステップ２５６及び２５８の双方において判定結果がＮＯのときには、ステップ２６０でプロセスは正常であると判定される。この後、制御の流れはステップ２５０に戻り、さらに１０分が経過するまで待機した後、上記した処理を再び実行する。

５．２ シミュレーション
本実施の形態について、プロセスの運転時間を１０分間隔にて１２分割し、１０分ごとに異常の判定を実施するシミュレーションを行なった。

本ケーススタディでは、対象とする異常として、ジャケット温度測定センサのドリフトを考慮した。すなわち時間と共にセンサのゼロ点が移動し、センサにて測定される測定値はＴ_ｊ（ｋ）＝Ｔ_ｊ０（ｋ）＋δＴ×（ｋ−ｋ_０）のように変化する。ここでＴ_ｊ０はゼロ点が固定されているときの測定値、ｋは測定時刻、ｋ_０はドリフト開始時間、δＴは単位時間当たりのドリフト量（Ｋ／分）である。
いま、新たな反応炉Ｒ_６を対象にドリフト開始時間ｋ_０＝２０分とし、δＴを変化させ異常を発生させた。対象とするデータは後掲のテーブルに示す６ケースであるが、このうちケース１、ケース２は正常データ（δＴ＝０）である。本実施の形態に係る装置の目的は、異常発生後できるだけ早く、正確に異常を検出することである。なお、Ｔ^２統計量、Ｑ統計量の管理限界はそれぞれ信頼区間９５％として設定し、採用した主成分数は５である。

ＫＭ−ＭＳＰＣ及び提案するＮＣ−ＭＳＰＣを用いた異常検出結果を表２−表１７に示す。表では、１１分〜２０分（第２インターバル）、２１分〜３０分（第３インターバル）、３１分〜４０分（第４インターバル）、４１分〜５０分（第５インターバル）におけるそれぞれの管理限界（ＣＬ）と各ケースでのＴ^２統計量、Ｑ統計量を示す。

この結果より、ＫＭ−ＭＳＰＣ、ＮＣ−ＭＳＰＣ双方ともに、未だ異常の発生していない第２インターバルの全てのケースと、正常データであるケース１、ケース２とについては、その後も正しく正常であると判定している。しかし、ＫＭ−ＭＳＰＣでは第３インターバル、第４インターバルでケース４の異常を検出できていない。一方、ＮＣ−ＭＳＰＣでは第３インターバルでケース４の異常を正しく検出できている。以上より、ＮＣ-ＭＳＰＣを採用した本実施の形態がＫＭ−ＭＳＰＣより素早く異常を検出することができることがわかる。

今回開示された実施の形態は単に例示であって、本発明が上記した実施の形態のみに制限されるわけではない。本発明の範囲は、発明の詳細な説明の記載を参酌した上で、特許請求の範囲の各請求項によって示され、そこに記載された文言と均等の意味及び範囲内での全ての変更を含む。

Ａｐｐｅｎｄｉｘ Ａ
バッチプロセスモデル
シミュレーションにて用いたバッチプロセスのモデルを示す。プロセスにおける収支式等は以下の通りである。

ここで、ａ（ｋ）〜Ｎ（０，０．１）であり、ｋは測定ステップ数を表す。ただし、Ｎ（ｍ，σ）は平均0標準偏差σの正規分布にしたがう乱数である。温度コントローラからジャケット温度への動特性は、（４５）式で与えられる。

Ｋ＝１／１６、τ＝２分である。モデルのパラメータをテーブル１８に、初期状態をテーブル１９に、それぞれ示す。

８０ プラント
８２，８４，８６ 反応炉
９２，９４，９６ コントローラ
９８ 製品分子量推定装置
１４０ ＣＰＵ

Claims (8)

1. プラントの稼動条件と、当該稼動条件での前記プラントの過去の稼動結果から得られるデータ項目の実測値との組合せからなる実測サンプルを複数個記憶するためのデータベースと、
   前記データベースに記憶された前記複数個の実測サンプルを、複数個のクラスにクラスタリングするためのクラスタリング手段と、
   前記クラスタリング手段によりクラスタリングされた各クラスに対して、前記プラントにおける稼動条件と前記データ項目の値との関係を示す統計的モデルを構築するためのモデル構築手段と、
   新たな稼動条件が与えられたことに応答して、前記統計的モデルと、前記新たな稼動条件とに基づいて、前記プラントの制御に関連する情報を算出するための制御情報算出手段とを含む、プラント制御情報生成装置であって、
   前記クラスタリング手段は、
   前記複数個の実測サンプルの任意の２個の間の相関の程度を成分とする行列を、前記データベースに記憶された前記実測サンプルから算出するための行列算出手段と、
   前記行列算出手段により算出された前記行列を類似度行列として、前記複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、前記複数個の実測サンプルを前記複数個のクラスにクラスタリングするための手段とを含み、
   前記行列算出手段は、
   前記プラントの稼動条件と前記データ項目とにより定義された実測ベクトルのベクトル空間において、前記複数個の実測サンプルにより定められる複数個の実測ベクトルのうち、基準となる実測サンプルに対応する実測ベクトルを他の実測ベクトルから減算するための減算手段と、
   前記減算手段により減算された後の前記他の実測ベクトルの任意の２個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関を求めるための相関算出手段と、
   前記減算手段による減算と、前記相関算出手段による相関の算出とを、前記基準となる実測サンプルを変えながら前記複数個の実測サンプルの全てに対して行なうための繰返算出手段と、
   前記繰返算出手段により算出された、前記複数個の実測サンプルの任意の２個の間で算出される複数個の相関に基づいて、前記行列を決定するための手段とを含む、プラント制御情報生成装置。
2. 前記行列を決定するための手段は、前記繰返算出手段により、前記複数個の実測サンプルの任意の２個の間で算出された前記複数個の相関を合計した値を成分とする行列を前記行列として出力するための手段を含む、請求項１に記載のプラント制御情報生成装置。
3. 前記相関算出手段は、
   前記減算手段により減算された後の前記他の実測ベクトルの任意の２個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関係数を求めるための相関係数算出手段と、
   前記相関係数算出手段により算出された相関係数が所定のしきい値より大きいか否かにしたがって、前記相関の値を２値化するための手段とを含む、請求項１又は請求項２に記載のプラント制御情報生成装置。
4. 前記制御情報算出手段は、
   新たな稼動条件が与えられたことに応答して、当該稼動条件が前記複数個のクラスのいずれに属するかを判定するための手段と、
   当該クラスに対応する統計的モデルを、前記新たな稼動条件に適用することにより、前記新たな稼動条件の下での前記プラントにより生成される物質の属性を予測し、予測情報を出力するための手段とを含む、請求項１〜請求項３のいずれかに記載のプラント制御情報生成装置。
5. 前記プラント制御情報生成装置は前記プラントの稼動状態を測定するための測定手段をさらに含み、
   前記制御情報算出手段は、前記測定手段による測定結果が新たな稼動条件として与えられたことに応答して、当該稼動条件と前記複数個のクラスとを比較することにより、前記プラントの状態が異常か否かを判定するための異常判定手段を含む、請求項１〜請求項３のいずれかに記載のプラント制御情報生成装置。
6. 前記異常判定手段は、
   前記測定手段による測定結果が新たな稼動条件として与えられたことに応答して、前記複数個のクラスの統計的モデルのうち、当該稼動条件が適合するものを判定するための適合モデル判定手段と、
   前記適合モデル判定手段により、前記新たな稼動条件が前記複数個のクラスの統計的モデルのいずれにも適合しないと判定されたことに応答して、前記プラントの異常を示す信号を生成し出力するための手段とを含む、請求項５に記載のプラント制御情報生成装置。
7. コンピュータにより実行されると、当該コンピュータを、
   プラントの稼動条件と、当該稼動条件での前記プラントの過去の稼動結果から得られるデータ項目の実測値との組合せからなる実測サンプルを複数個記憶するためのデータベースと、
   前記データベースに記憶された前記複数個の実測サンプルを、複数個のクラスにクラスタリングするためのクラスタリング手段と、
   前記クラスタリング手段によりクラスタリングされた各クラスに対して、前記プラントにおける稼動条件と前記データ項目の値との関係を示す統計的モデルを構築するためのモデル構築手段と、
   新たな稼動条件が与えられたことに応答して、前記統計的モデルと、前記新たな稼動条件とに基づいて、前記プラントの制御に関連する情報を算出するための制御情報算出手段として機能させるコンピュータプログラムであって、
   前記クラスタリング手段は、
   前記複数個の実測サンプルの任意の２個の間の相関の程度を成分とする行列を、前記データベースに記憶された前記実測サンプルから算出するための行列算出手段と、
   前記行列算出手段により算出された前記行列を類似度行列として、前記複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、前記複数個の実測サンプルを前記複数個のクラスにクラスタリングするための手段とを含み、
   前記行列算出手段は、
   前記プラントの稼動条件と前記データ項目とにより定義された実測ベクトルのベクトル空間において、前記複数個の実測サンプルにより定められる複数個の実測ベクトルのうち、基準となる実測サンプルに対応する実測ベクトルを他の実測ベクトルから減算するための減算手段と、
   前記減算手段により減算された後の前記他の実測ベクトルの任意の２個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関を求めるための相関算出手段と、
   前記減算手段による減算と、前記相関算出手段による相関の算出とを、前記基準となる実測サンプルを変えながら前記複数個の実測サンプルの全てに対して行なうための繰返算出手段と、
   前記繰返算出手段により算出された、前記複数個の実測サンプルの任意の２個の間で算出される複数個の相関に基づいて、前記行列を決定するための手段とを含む、コンピュータプログラム。
8. プラントの稼動条件と、当該稼動条件での前記プラントの過去の稼動結果から得られるデータ項目の実測値との組合せからなる実測サンプルを複数個記憶するためのデータベースと、
   前記データベースに記憶された前記複数個の実測サンプルを、複数個のクラスにクラスタリングするためのクラスタリング手段と、
   前記クラスタリング手段によりクラスタリングされた各クラスに対して、前記プラントにおける稼動条件と前記データ項目の値との関係を示す統計的モデルを構築するためのモデル構築手段と、
   新たな稼動条件が与えられたことに応答して、前記統計的モデルと、前記新たな稼動条件とに基づいて、前記プラントの制御に関連する情報を算出するための制御情報算出手段とを含む、プラント制御情報生成装置におけるプラント制御情報生成方法であって、
   前記クラスタリング手段は、
   前記複数個の実測サンプルの任意の２個の間の相関の程度を成分とする行列を、前記データベースに記憶された前記実測サンプルから算出するための行列算出手段と、
   前記行列算出手段により算出された前記行列を類似度行列として、前記複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、前記複数個の実測サンプルを前記複数個のクラスにクラスタリングするための手段とを含み、
   前記行列算出手段は、
   前記プラントの稼動条件と前記データ項目とにより定義された実測ベクトルのベクトル空間において、前記複数個の実測サンプルにより定められる複数個の実測ベクトルのうち、基準となる実測サンプルに対応する実測ベクトルを他の実測ベクトルから減算するための減算手段と、
   前記減算手段により減算された後の前記他の実測ベクトルの任意の２個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関を求めるための相関算出手段と、
   前記減算手段による減算と、前記相関算出手段による相関の算出とを、前記基準となる実測サンプルを変えながら前記複数個の実測サンプルの全てに対して行なうための繰返算出手段と、
   前記繰返算出手段により算出された、前記複数個の実測サンプルの任意の２個の間で算出される複数個の相関に基づいて、前記行列を決定するための手段とを含み、
   前記プラント制御情報生成方法は、
   前記クラスタリング手段が、前記データベースに記憶された前記複数個の実測サンプルを、複数個のクラスにクラスタリングするステップと、
   前記モデル構築手段が、前記クラスタリングするステップにおいてクラスタリングされた各クラスに対して、前記プラントにおける稼動条件と前記データ項目の値との関係を示す統計的モデルを構築するステップと、
   前記制御情報算出手段が新たな稼動条件が与えられたことに応答して、前記統計的モデルと、前記新たな稼動条件とに基づいて、前記プラントの制御に関連する情報を算出するステップとを含み、
   前記クラスタリングするステップは、
   前記行列算出手段が、前記複数個の実測サンプルの任意の２個の間の相関の程度を成分とする行列を、前記データベースに記憶された前記実測サンプルから算出するステップと、
   前記クラスタリングするための手段が、前記算出するステップにおいて算出された前記行列を類似度行列として、前記複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、前記複数個の実測サンプルを前記複数個のクラスにクラスタリングするステップとを含み、
   前記算出するステップは、
   前記減算手段が、前記プラントの稼動条件と前記データ項目とにより定義された実測ベクトルのベクトル空間において、前記複数個の実測サンプルにより定められる複数個の実測ベクトルのうち、基準となる実測サンプルに対応する実測ベクトルを他の実測ベクトルから減算するステップと、
   前記相関算出手段が、前記減算するステップにおいて減算された後の前記他の実測ベクトルの任意の２個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関を求めるステップと、
   前記繰返算出手段が、前記減算するステップにおける減算と、前記相関を求めるステップにおける相関の算出とを、前記基準となる実測サンプルを変えながら前記複数個の実測サンプルの全てに対して行なう繰返算出ステップと、
   前記繰返算出ステップにおいて算出された、前記複数個の実測サンプルの任意の２個の間で算出される複数個の相関に基づいて、前記行列を決定するステップとを含む、プラント制御情報生成方法。