Top > Search of Japanese Patents > DEVICE AND METHOD FOR GENERATING PLANT CONTROL INFORMATION, AND COMPUTER PROGRAM THEREFOR > Specification

Specification :(In Japanese)プラント制御情報生成装置及び方法、並びにそのためのコンピュータプログラム

Country (In Japanese)日本国特許庁(JP)
Gazette (In Japanese)特許公報(B2)
Patent Number P5457737
Publication number P2011-008562A
Date of registration Jan 17, 2014
Date of issue Apr 2, 2014
Date of publication of application Jan 13, 2011
Title of the invention, or title of the device (In Japanese)プラント制御情報生成装置及び方法、並びにそのためのコンピュータプログラム
IPC (International Patent Classification) G05B  23/02        (2006.01)
FI (File Index) G05B 23/02 P
G05B 23/02 301X
Number of claims or invention 8
Total pages 42
Application Number P2009-151745
Date of filing Jun 26, 2009
Date of request for substantive examination Jun 20, 2012
Patentee, or owner of utility model right (In Japanese)【識別番号】504132272
【氏名又は名称】国立大学法人京都大学
Inventor, or creator of device (In Japanese)【氏名】加納 学
【氏名】藤原 幸一
Representative (In Japanese)【識別番号】100099933、【弁理士】、【氏名又は名称】清水 敏
Examiner (In Japanese)【審査官】青山 純
Document or reference (In Japanese)特開2007-272897(JP,A)
特開2006-276924(JP,A)
向井洋介,変量間の相関に着目したクラスタリング手法およびその多変量統計モデリングへの利用,化学工学会 研究発表講演要旨集 化学工学会第40回秋季大会,日本,社団法人 化学工学会,2008年 8月24日,K302頁,[平成25年7月4日検索],インターネット<URL:https://www.jstage.jst.go.jp/article/scej/2008f/0/2008f_0_499/_pdf>
Field of search G05B 23/02
Scope of claims (In Japanese)【請求項1】
プラントの稼動条件と、当該稼動条件での前記プラントの過去の稼動結果から得られるデータ項目の実測値との組合せからなる実測サンプルを複数個記憶するためのデータベースと、
前記データベースに記憶された前記複数個の実測サンプルを、複数個のクラスにクラスタリングするためのクラスタリング手段と、
前記クラスタリング手段によりクラスタリングされた各クラスに対して、前記プラントにおける稼動条件と前記データ項目の値との関係を示す統計的モデルを構築するためのモデル構築手段と、
新たな稼動条件が与えられたことに応答して、前記統計的モデルと、前記新たな稼動条件とに基づいて、前記プラントの制御に関連する情報を算出するための制御情報算出手段とを含む、プラント制御情報生成装置であって、
前記クラスタリング手段は、
前記複数個の実測サンプルの任意の2個の間の相関の程度を成分とする行列を、前記データベースに記憶された前記実測サンプルから算出するための行列算出手段と、
前記行列算出手段により算出された前記行列を類似度行列として、前記複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、前記複数個の実測サンプルを前記複数個のクラスにクラスタリングするための手段とを含み、
前記行列算出手段は、
前記プラントの稼動条件と前記データ項目とにより定義された実測ベクトルのベクトル空間において、前記複数個の実測サンプルにより定められる複数個の実測ベクトルのうち、基準となる実測サンプルに対応する実測ベクトルを他の実測ベクトルから減算するための減算手段と、
前記減算手段により減算された後の前記他の実測ベクトルの任意の2個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関を求めるための相関算出手段と、
前記減算手段による減算と、前記相関算出手段による相関の算出とを、前記基準となる実測サンプルを変えながら前記複数個の実測サンプルの全てに対して行なうための繰返算出手段と、
前記繰返算出手段により算出された、前記複数個の実測サンプルの任意の2個の間で算出される複数個の相関に基づいて、前記行列を決定するための手段とを含む、プラント制御情報生成装置。
【請求項2】
前記行列を決定するための手段は、前記繰返算出手段により、前記複数個の実測サンプルの任意の2個の間で算出された前記複数個の相関を合計した値を成分とする行列を前記行列として出力するための手段を含む、請求項1に記載のプラント制御情報生成装置。
【請求項3】
前記相関算出手段は、
前記減算手段により減算された後の前記他の実測ベクトルの任意の2個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関係数を求めるための相関係数算出手段と、
前記相関係数算出手段により算出された相関係数が所定のしきい値より大きいか否かにしたがって、前記相関の値を2値化するための手段とを含む、請求項1又は請求項2に記載のプラント制御情報生成装置。
【請求項4】
前記制御情報算出手段は、
新たな稼動条件が与えられたことに応答して、当該稼動条件が前記複数個のクラスのいずれに属するかを判定するための手段と、
当該クラスに対応する統計的モデルを、前記新たな稼動条件に適用することにより、前記新たな稼動条件の下での前記プラントにより生成される物質の属性を予測し、予測情報を出力するための手段とを含む、請求項1~請求項3のいずれかに記載のプラント制御情報生成装置。
【請求項5】
前記プラント制御情報生成装置は前記プラントの稼動状態を測定するための測定手段をさらに含み、
前記制御情報算出手段は、前記測定手段による測定結果が新たな稼動条件として与えられたことに応答して、当該稼動条件と前記複数個のクラスとを比較することにより、前記プラントの状態が異常か否かを判定するための異常判定手段を含む、請求項1~請求項3のいずれかに記載のプラント制御情報生成装置。
【請求項6】
前記異常判定手段は、
前記測定手段による測定結果が新たな稼動条件として与えられたことに応答して、前記複数個のクラスの統計的モデルのうち、当該稼動条件が適合するものを判定するための適合モデル判定手段と、
前記適合モデル判定手段により、前記新たな稼動条件が前記複数個のクラスの統計的モデルのいずれにも適合しないと判定されたことに応答して、前記プラントの異常を示す信号を生成し出力するための手段とを含む、請求項5に記載のプラント制御情報生成装置。
【請求項7】
コンピュータにより実行されると、当該コンピュータを、
プラントの稼動条件と、当該稼動条件での前記プラントの過去の稼動結果から得られるデータ項目の実測値との組合せからなる実測サンプルを複数個記憶するためのデータベースと、
前記データベースに記憶された前記複数個の実測サンプルを、複数個のクラスにクラスタリングするためのクラスタリング手段と、
前記クラスタリング手段によりクラスタリングされた各クラスに対して、前記プラントにおける稼動条件と前記データ項目の値との関係を示す統計的モデルを構築するためのモデル構築手段と、
新たな稼動条件が与えられたことに応答して、前記統計的モデルと、前記新たな稼動条件とに基づいて、前記プラントの制御に関連する情報を算出するための制御情報算出手段として機能させるコンピュータプログラムであって、
前記クラスタリング手段は、
前記複数個の実測サンプルの任意の2個の間の相関の程度を成分とする行列を、前記データベースに記憶された前記実測サンプルから算出するための行列算出手段と、
前記行列算出手段により算出された前記行列を類似度行列として、前記複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、前記複数個の実測サンプルを前記複数個のクラスにクラスタリングするための手段とを含み、
前記行列算出手段は、
前記プラントの稼動条件と前記データ項目とにより定義された実測ベクトルのベクトル空間において、前記複数個の実測サンプルにより定められる複数個の実測ベクトルのうち、基準となる実測サンプルに対応する実測ベクトルを他の実測ベクトルから減算するための減算手段と、
前記減算手段により減算された後の前記他の実測ベクトルの任意の2個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関を求めるための相関算出手段と、
前記減算手段による減算と、前記相関算出手段による相関の算出とを、前記基準となる実測サンプルを変えながら前記複数個の実測サンプルの全てに対して行なうための繰返算出手段と、
前記繰返算出手段により算出された、前記複数個の実測サンプルの任意の2個の間で算出される複数個の相関に基づいて、前記行列を決定するための手段とを含む、コンピュータプログラム。
【請求項8】
プラントの稼動条件と、当該稼動条件での前記プラントの過去の稼動結果から得られるデータ項目の実測値との組合せからなる実測サンプルを複数個記憶するためのデータベースと、
前記データベースに記憶された前記複数個の実測サンプルを、複数個のクラスにクラスタリングするためのクラスタリング手段と、
前記クラスタリング手段によりクラスタリングされた各クラスに対して、前記プラントにおける稼動条件と前記データ項目の値との関係を示す統計的モデルを構築するためのモデル構築手段と、
新たな稼動条件が与えられたことに応答して、前記統計的モデルと、前記新たな稼動条件とに基づいて、前記プラントの制御に関連する情報を算出するための制御情報算出手段とを含む、プラント制御情報生成装置におけるプラント制御情報生成方法であって、
前記クラスタリング手段は、
前記複数個の実測サンプルの任意の2個の間の相関の程度を成分とする行列を、前記データベースに記憶された前記実測サンプルから算出するための行列算出手段と、
前記行列算出手段により算出された前記行列を類似度行列として、前記複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、前記複数個の実測サンプルを前記複数個のクラスにクラスタリングするための手段とを含み、
前記行列算出手段は、
前記プラントの稼動条件と前記データ項目とにより定義された実測ベクトルのベクトル空間において、前記複数個の実測サンプルにより定められる複数個の実測ベクトルのうち、基準となる実測サンプルに対応する実測ベクトルを他の実測ベクトルから減算するための減算手段と、
前記減算手段により減算された後の前記他の実測ベクトルの任意の2個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関を求めるための相関算出手段と、
前記減算手段による減算と、前記相関算出手段による相関の算出とを、前記基準となる実測サンプルを変えながら前記複数個の実測サンプルの全てに対して行なうための繰返算出手段と、
前記繰返算出手段により算出された、前記複数個の実測サンプルの任意の2個の間で算出される複数個の相関に基づいて、前記行列を決定するための手段とを含み、
前記プラント制御情報生成方法は、
前記クラスタリング手段が、前記データベースに記憶された前記複数個の実測サンプルを、複数個のクラスにクラスタリングするステップと、
前記モデル構築手段が、前記クラスタリングするステップにおいてクラスタリングされた各クラスに対して、前記プラントにおける稼動条件と前記データ項目の値との関係を示す統計的モデルを構築するステップと、
前記制御情報算出手段が新たな稼動条件が与えられたことに応答して、前記統計的モデルと、前記新たな稼動条件とに基づいて、前記プラントの制御に関連する情報を算出するステップとを含み、
前記クラスタリングするステップは、
前記行列算出手段が、前記複数個の実測サンプルの任意の2個の間の相関の程度を成分とする行列を、前記データベースに記憶された前記実測サンプルから算出するステップと、
前記クラスタリングするための手段が、前記算出するステップにおいて算出された前記行列を類似度行列として、前記複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、前記複数個の実測サンプルを前記複数個のクラスにクラスタリングするステップとを含み、
前記算出するステップは、
前記減算手段が、前記プラントの稼動条件と前記データ項目とにより定義された実測ベクトルのベクトル空間において、前記複数個の実測サンプルにより定められる複数個の実測ベクトルのうち、基準となる実測サンプルに対応する実測ベクトルを他の実測ベクトルから減算するステップと、
前記相関算出手段が、前記減算するステップにおいて減算された後の前記他の実測ベクトルの任意の2個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関を求めるステップと、
前記繰返算出手段が、前記減算するステップにおける減算と、前記相関を求めるステップにおける相関の算出とを、前記基準となる実測サンプルを変えながら前記複数個の実測サンプルの全てに対して行なう繰返算出ステップと、
前記繰返算出ステップにおいて算出された、前記複数個の実測サンプルの任意の2個の間で算出される複数個の相関に基づいて、前記行列を決定するステップとを含む、プラント制御情報生成方法。
Detailed description of the invention (In Japanese)【技術分野】
【0001】
この発明は、化学プラント、半導体製造設備、製鉄プラント等の制御技術に関し、特に、過去の複数の稼動時の稼動条件とその結果とに基づいて、プラントの制御の精度を高めるための技術に関する。
【背景技術】
【0002】
化学プラントでは、複数の装置(反応炉等)を使用して同一の製品を生産することが多い。この場合、生産される製品が均一となるように、各装置の稼動条件を適切に設定することが必要である。しかし、装置内におけるプロセスの進行とともにプロセスの特性が変化すること、及び複数の装置を並列で運転している場合に、それぞれの装置のカタログスペックが同一でも実際には個体差があること等により、生産される製品が均一となるような条件を装置ごとに精度よく決定したり、生産される製品の組成を稼動条件から精度よく予測したりすることはむずかしい。同様の問題は、化学プラントに限らず、半導体製造設備、金属精錬プラント、水処理施設、及びガスプラント等でも生じる。また、何らかのものを製造する施設だけではなく、発電・変電施設のように非有体物の特性を変化させるための施設でもこうした問題は生じ得る。以下、本明細書では、何らかのものを製造する施設だけでなく、非有体物の特性を変化させる施設をも含めて「プラント」と呼ぶ。また、複合的な施設だけでなく、単独の設備であって、稼動条件によって異なる挙動を示すものもここでいう「プラント」に含めるものとする。
【0003】
これまで、本願発明者は、プロセスの特性変化が測定パラメータ間の相関関係の違いとして表現されることに着目し、パラメータ間の相関関係を考慮したモデリングがソフトセンサ設計に有効であることを示してきた。本願発明者はさらに、装置間の個体差に対応できる手法として、パラメータ間の相関関係の類似度を考慮して、過去の装置の稼動パラメータをクラスタリングする手法(向井法)と、相関識別法(NeareStCorrelationMethod;NC法)とを開発した(非特許文献1、2)。
【先行技術文献】
【0004】

【非特許文献1】向井、藤原、加納、長谷部:変数間の相関に着目したクラスタリング手法の開発およびその多変量統計モデリングへの利用、計測自動制御学会第8回制御部門大会、京都大学(2008)
【非特許文献2】藤原、加納、長谷部:相関関係を考慮したパターン認識手法の開発とソフトセンサ設計への適用、計測自動制御学会第9回制御部門大会、広島大学(2009)
【発明の概要】
【発明が解決しようとする課題】
【0005】
向井法は、パラメータ間の相関関係を指標としてデータをクラスタリングできるが、繰返計算を行なう必要がある。一方、NC法は、線形空間におけるクエリとサンプルの幾何学的な関係に着目した手法であり、教師信号を用いることなくクエリと類似の相関関係を有するサンプルを検出できるという特徴を持つが、クラスタリングに用いることはできなかった。
【0006】
ここでいうパラメータ間の相関関係とは、以下のようなものをいう。以下、説明を理解しやすくするために、2つのパラメータのみからなるサンプルを扱う。図1を参照して、互いに相関関係を持つサンプルの集合が3つあるものとする。図1においては、第1の集合は黒丸で、第2の集合は白丸で、第3の集合は「+」記号で、それぞれ表されている。ここでは、各集合を「クラス」と呼ぶ。各クラスに属するサンプルの間の相関関係はクラスごとに異なっている。あるクラスについて新たなサンプルが得られれば、新たなサンプルを使用して、そのクラスに適した制御を行なうことができる。したがって、ある新たなサンプルが得られたときに、そのサンプルを的確なクラスに分類することが必要である。
【0007】
通常、このような場合には、予めサンプルの集合をクラスに分類しておき、各クラスとサンプルとの間の「距離」を定義しておく手法が用いられる。新たなサンプルが得られると、そのサンプルと各クラスとの間の「距離」を算出して、最も距離が小さなクラスに新たなサンプルを分類する。
【0008】
得られたサンプルの集合がどのクラスに属するかが明確であれば、上記した手法を用いることができる。しかし実際のプラントの稼動時に得られる実測サンプルがどのクラスに属するかは明らかでない。つまり、実際のプラントから得られるサンプルのクラスを図1のように3つのクラスに分類(クラスタリング)する手法が必要である。
【0009】
クラスタリング手法として最も頻繁に使用されるものにk-平均法がある。図1に示されるサンプルの集合をk-平均法を用いてクラスタリングした結果を図2に示す。図2では、クラスタリングによりサンプルを3つのクラスに分類している。各クラスに分類されたサンプルをそれぞれ黒丸、白丸及び「+」記号で表している。
【0010】
図2により明らかなように、k-平均法を用いた場合、サンプルを正確にクラスタリングすることはできない。これはk-平均法が基本的にサンプルの間の距離に基づいてクラスタリングを行なっているためである。すなわち、各クラスのサンプルが集中する中央部付近ではサンプルが全て同じクラスタに分類され、その両側に散在するサンプルはそれぞれ別々のクラスタに分類される。
【0011】
k-平均法を用いたクラスタリングとは別のクラスタリング手法も存在する。図3は、スペクトラル・クラスタリングと呼ばれるクラスタリング手法を用いて図1に示されるサンプルをクラスタリングした結果を示す。このスペクトラル・クラスタリングの手法についてはその詳細を後述する。図3において、各クラスに属するサンプルの表示方法は図1及び図2と同様である。このスペクトラル・クラスタリングは、後述するようにグラフ理論に基づくものであって、k-平均法と比較すると単純に距離に基づいてクラスタリングするものではない。したがって、クラスタリングをより正確に行なえる可能性がある。しかし、図3に示される例ではクラスタリングが正しく行なわれているとはいえない。
【0012】
過去の実測データを統計処理することにより、プラントの制御を行なう場合、実測データを正しくクラスタリングできなければ、制御の前提が失われることになり、正しい制御を実現できない。したがって、実測データを正しくクラスタリングすることができる技術が必要である。
【0013】
したがって本発明の目的は、サンプル間の相関関係を用いて、スペクトラル・クラスタリング手法によりクラスタリングしたサンプルの集合を用いて、精度よくプラントの制御を行なうことができるプラント制御情報生成装置及び方法、並びにそのためにコンピュータプログラムを提供することである。
【0014】
本発明に係るプラント制御情報生成装置は、プラントの稼動条件と、当該稼動条件でのプラントの過去の稼動結果から得られるデータ項目の実測値との組合せからなる実測サンプルを複数個記憶するためのデータベースと、データベースに記憶された複数個の実測サンプルを、複数個のクラスにクラスタリングするためのクラスタリング手段と、クラスタリング手段によりクラスタリングされた各クラスに対して、プラントにおける稼動条件とデータ項目の値との関係を示す統計的モデルを構築するためのモデル構築手段と、新たな稼動条件が与えられたことに応答して、統計的モデルと、新たな稼動条件とに基づいて、プラントに関する制御情報を算出するための制御情報算出手段とを含む。クラスタリング手段は、複数個の実測サンプルの任意の2個の間の相関の程度を成分とする行列を、データベースに記憶された実測サンプルから算出するための行列算出手段と、行列算出手段により算出された行列を類似度行列として、複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、複数個の実測サンプルを複数個のクラスにクラスタリングするための手段とを含む。
【0015】
このプラント制御情報生成装置においては、クラスタリング手段が、データベースに記憶された複数個の実測サンプルを複数個のクラスにクラスタリングし、クラスタリングされた後のクラスごとに、統計的モデルが算出される。クラスタリングの際、複数個の実測サンプルの任意の2個の間の相関の程度を成分とする行列を類似度行列として、スペクトラル・クラスタリングが行なわれる。サンプルの間の相関の相違に基づいて複数個のサンプルがクラスタリングされるため、単純にサンプル間の距離に基づいてクラスタリングされる場合と比較して、異なる相関関係にあるサンプルが正しく別々のクラスにクラスタリングされる確率が高くなる。その結果、正しいクラスタリングが可能になり、それらクラスから得られる統計的モデルは正確なものとなる。こうした処理は、新たな稼動条件が与えられ、新たな制御情報の生成が要求されるまでに行なっておけばよく、複雑な挙動を示すプラントについても、稼動時に、従来よりも正確な制御情報を生成することが可能になる。その結果、サンプル間の相関関係を用いて、スペクトラル・クラスタリング手法によりクラスタリングしたサンプルの集合を用いて、精度よくプラントの制御を行なうことができるプラント制御情報生成装置を提供することができる。
【0016】
好ましくは、行列算出手段は、プラントの稼動条件とデータ項目とにより定義された実測ベクトルのベクトル空間において、複数個の実測サンプルにより定められる複数個の実測ベクトルのうち、基準となる実測サンプルに対応する実測ベクトルを他の実測ベクトルから減算するための減算手段と、減算手段により減算された後の他の実測ベクトルの任意の2個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関を求めるための相関算出手段と、減算手段による減算と、相関算出手段による相関の算出とを、基準となる実測サンプルを変えながら複数個の実測サンプルの全てに対して行なうための繰返算出手段と、繰返算出手段により算出された、複数個の実測サンプルの任意の2個の間で算出される複数個の相関に基づいて、行列を決定するための手段とを含む。
【0017】
このように、実測サンプルの各々を基準ベクトルとして他の実測ベクトルから減算することにより、そのベクトルは実測ベクトル空間の原点に移動し、それに伴って他の実測ベクトルも平行移動する。その結果、基準実測ベクトルを通る平面は実測ベクトル空間の線形部分空間となり、その上に存在する実測ベクトルが互いにある相関関係を持つと判定できる。こうした判定を全ての実測ベクトルに対して行ない、さらに基準ベクトルとして全ての実測ベクトルを用いてこの処理を繰返すことで、実測ベクトルの間の相関関係がより明確となる。その結果、そうした相関に基づいて算出される類似度行列が、実測ベクトルの間の相関関係の相違を明確に表すものとなって、スペクトル・クラスタリングが正確に行なえる。
【0018】
より好ましくは、行列を決定するための手段は、繰返算出手段により、複数個の実測サンプルの任意の2個の間で算出された複数個の相関を合計した値を成分とする行列を行列として出力するための手段を含む。
【0019】
行列の各成分について、単純に複数個の相関を加算することで、複数個の実測サンプルの任意の2個の間の相関の程度を行列として表現できる。必要な計算量が少なくて済み、処理を早めてリアルタイムで実行することができる。
【0020】
相関算出手段は、減算手段により減算された後の他の実測ベクトルの任意の2個について、共通の線形部分空間上に存在するか否かを基準として、互いの相関係数を求めるための相関係数算出手段と、相関係数算出手段により算出された相関係数が所定のしきい値より大きいか否かにしたがって、相関の値を2値化するための手段とを含んでもよい。
【0021】
2値化することによって、加算処理がより単純なものとなり、計算量をより削減することができる。
【0022】
本発明のある実施の形態では、制御情報算出手段は、新たな稼動条件が与えられたことに応答して、当該稼動条件が複数個のクラスのいずれに属するかを判定するための手段と、当該クラスに対応する統計的モデルを、新たな稼動条件に適用することにより、新たな稼動条件の下でのプラントにより生成される物質の属性を予測し、予測情報を出力するための手段とを含む。
【0023】
本発明の別の実施の形態では、プラント制御情報生成装置はプラントの稼動状態を測定するための測定手段をさらに含む。制御情報算出手段は、測定手段による測定結果が新たな稼動条件として与えられたことに応答して、当該稼動条件と複数個のクラスとを比較することにより、プラントの状態が異常か否かを判定するための異常判定手段を含む。
【0024】
好ましくは、異常判定手段は、測定手段による測定結果が新たな稼動条件として与えられたことに応答して、複数個のクラスの統計的モデルのうち、当該稼動条件が適合するものを判定するための適合モデル判定手段と、適合モデル判定手段により、新たな稼動条件が複数個のクラスの統計的モデルのいずれにも適合しないと判定されたことに応答して、プラントの異常を示す信号を生成し出力するための手段とを含む。
【0025】
本発明の第2の局面に係るコンピュータプログラムは、コンピュータにより実行されると、当該コンピュータを、プラントの稼動条件と、当該稼動条件でのプラントの過去の稼動結果から得られるデータ項目の実測値との組合せからなる実測サンプルを複数個記憶するためのデータベースと、データベースに記憶された複数個の実測サンプルを、複数個のクラスにクラスタリングするためのクラスタリング手段と、クラスタリング手段によりクラスタリングされた各クラスに対して、プラントにおける稼動条件とデータ項目の値との関係を示す統計的モデルを構築するためのモデル構築手段と、新たな稼動条件が与えられたことに応答して、統計的モデルと、新たな稼動条件とに基づいて、プラントに関する制御情報を算出するための制御情報算出手段として機能させる。このプログラムは、クラスタリング手段として、コンピュータを、複数個の実測サンプルの任意の2個の間の相関の程度を成分とする行列を、データベースに記憶された実測サンプルから算出するための行列算出手段と、行列算出手段により算出された行列を類似度行列として、複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、複数個の実測サンプルを複数個のクラスにクラスタリングするための手段として機能させる。
【0026】
本発明の第3の局面に係るプラント制御情報生成方法は、プラントの稼動条件と、当該稼動条件でのプラントの過去の稼動結果から得られるデータ項目の実測値との組合せからなる実測サンプルを複数個記憶するためのデータベースと、データベースに記憶された複数個の実測サンプルを、複数個のクラスにクラスタリングするためのクラスタリング手段と、クラスタリング手段によりクラスタリングされた各クラスに対して、プラントにおける稼動条件とデータ項目の値との関係を示す統計的モデルを構築するためのモデル構築手段と、新たな稼動条件が与えられたことに応答して、統計的モデルと、新たな稼動条件とに基づいて、プラントに関する制御情報を算出するための制御情報算出手段とを含む、プラント制御情報生成装置におけるプラント制御情報生成方法である。クラスタリング手段は、複数個の実測サンプルの任意の2個の間の相関の程度を成分とする行列を、データベースに記憶された実測サンプルから算出するための行列算出手段と、行列算出手段により算出された行列を類似度行列として、複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、複数個の実測サンプルを複数個のクラスにクラスタリングするための手段とを含む。このプラント制御情報生成方法は、クラスタリング手段が、データベースに記憶された複数個の実測サンプルを、複数個のクラスにクラスタリングするステップと、モデル構築手段が、クラスタリングするステップにおいてクラスタリングされた各クラスに対して、プラントにおける稼動条件とデータ項目の値との関係を示す統計的モデルを構築するステップと、制御情報算出手段が新たな稼動条件が与えられたことに応答して、統計的モデルと、新たな稼動条件とに基づいて、プラントに関する制御情報を算出するステップとを含む。クラスタリングするステップは、行列算出手段が、複数個の実測サンプルの任意の2個の間の相関の程度を成分とする行列を、データベースに記憶された実測サンプルから算出するステップと、クラスタリングするための手段が、算出するステップにおいて算出された行列を類似度行列として、複数個の実測サンプルに対するスペクトル・クラスタリングを行なって、複数個の実測サンプルを複数個のクラスにクラスタリングするステップとを含む。
【図面の簡単な説明】
【0027】
【図1】図1は、相関関係の違いに基づき実測サンプルを真のクラスにクラスタリングした結果を示すグラフである。
【図2】図2は、図1に示す実測サンプルを対象にk-平均法によってクラスタリングした結果を示すグラフである。
【図3】図3は、図1に示す実測サンプルを対象に従来のスペクトル・クラスタリング法によってクラスタリングした結果を示すグラフである。
【図4】図4は、8個のノードを持つ無向グラフの例を示す図である。
【図5】図5は、アフィン部分空間P内におけるサンプル間の相関関係を説明するための模式図である。
【図6】図6は、図5に示す変数群を、サンプルx1が原点と一致するように平行移動した後のサンプルの分布を説明するための模式図である。
【図7】図7は、図5に示す変数群を、サンプルx3が原点と一致するように平行移動した後のサンプルの分布を説明するための模式図である。
【図8】変数間の相関関係を指標として、スペクトラル・クラスタリングのための類似度行列を求める、本発明の位置実施の形態に係る装置で採用したアルゴリズムを説明するためのフローチャートである。
【図9】図9は、図1に示す測定サンプルを対象に本願発明の一実施の形態によってクラスタリングした結果を示すグラフである。
【図10】図10は、本発明の一実施の形態に係るバッチプロセスの概略構成図を示す模式図である。
【図11】図11は、図10に示すバッチプロセスにおいて、製品の分子量を推定する、本願発明の一実施の形態に係るソフトセンサを実現するためのコンピュータのブロック図である。
【図12】図12は、本願発明の一実施の形態において、実測サンプルをクラスタリングするためのプログラムの制御構造を示すフローチャートである。
【図13】図13は、本願発明の一実施の形態において、新たなサンプルにしたがって、バッチプロセスで産生される製品の分子量を推定するプログラムの制御構造を示すフローチャートである。
【図14】図14は、バッチプロセスにおける反応炉温度とジャケット温度との時系列変化を示すグラフである。
【図15】図15は、k-平均法によるクラスタリング結果を示すグラフである。
【図16】図16は、本発明の一実施の形態に係るシステムにおけるクラスタリング結果を示すグラフである。
【図17】図17は、k-平均法及び本発明の第1の実施の形態に係るソフトセンサの推定性能を示すグラフでる。
【図18】図18は、本願発明の第2の実施の形態に係るプロセス異常の検出装置における準備処理を実現するプログラムの制御構造を示すフローチャートである。
【図19】図19は、本願発明の第2の実施の形態に係るプロセス異常の検出装置における、プロセス異常を検出するためのプログラムの制御構造を示すフローチャートである。
【発明を実施するための形態】
【0028】
1.最初に
まず、スペクトラル・クラスタリングの基礎であるグラフ理論を簡単に紹介し、次にスペクトラル・クラスタリングの代表的なアルゴリズムを述べる。

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

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

【0031】
数学的なグラフの定義を与える。集合V,Eと、Eの元に2つのVの元を対応させる写像f:E→V×Vがあるとき、G=(V,E)を有向グラフという。Vの元をGのノード(頂点)、Eの元をエッジ(辺)と呼ぶ。一方、P(V)をVのべき集合とし、Eの元にP(V)の部分集合を対応させる写像g:E→P(V)があって、e∈Eに対しg(e)={v1,v2}(v1,v2∈V)であるとき、G=(V,E)を無向グラフと呼ぶ。g(e)={v1,v2}を満たすv1とv2のペアを「互いに隣接している」といい、v1~v2で表す。グラフによっては、エッジに重みがついている場合がある。これを重み付きグラフという。8個のノードを持つ無向グラフの例を、図4に示す。

【0032】
図4に示すグラフでは、各ノードに1~8のインデックスを付してあり、各ノードをインデックスで呼ぶものとする。すなわち、V={1,2,3,4,5,6.7,8}、E={{1,2}、{1,5}、{2,5}{2,3}、{3,4}、{4,5}、{4,6}、{6,7}、{6.8}、{7,8}}と表すことができる。

【0033】
さらに、2つのグラフG=(V,E)とG′=(V′,E′)とに対して、全単射h:V→V′でv1~v2←→h(v1)~h(v2)を満たすものが存在するとき、G、G′は同型であるという。そのような写像hを同型写像と呼ぶ。同型写像の例としては、ノードのインデックスのみを入れ替える写像がある。

【0034】
グラフの表現には行列表現が便利である。特に、グラフは「隣接行列」と呼ばれる行列で表現することが多い。無向グラフG=(V,E)のノードに重複のないインデックスが与えられているとき、次の式によってこのグラフGの隣接行列Aを定義する。

【0035】
【数1】
JP0005457737B2_000002t.gif
この定義から明らかなように、無向グラフの場合には隣接行列は対称行列である。有向グラフの場合、ノードv1からv2に向かうエッジが存在するときに行列の(A)v1,v2=1、そうでないときは(A)v1,v2=0とすればよい。また、重み付きグラフの場合は(A)v1,v2にノードv1からv2に向かうエッジに与えられた重みを代入する。なお、隣接行列の対角成分には0を割り当てることが多い。たとえば、図4で表現される無向グラフの隣接行列A∈R8×8は、次式のように表現できる。

【0036】
【数2】
JP0005457737B2_000003t.gif

無向グラフの隣接行列は対称行列であるからその固有値は全て実数であり、固有ベクトルは互いに直交する。無向グラフGの隣接行列の固有値をλ1<λ2…<λNとし、それらの重複度をm1,…,mNとして得られる以下の表をグラフGのスペクトルと呼ぶ。

【0037】
【数3】
JP0005457737B2_000004t.gif
同型なグラフのスペクトルは一致するという性質があり(ただし、逆は成り立たない)、スペクトルはグラフの性質の分類に使われる。

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

【0039】
1.2 Max-Mincut法
スペクトラル・クラスタリングは、重み付きグラフのグラフ分割法であり、いくつかのエッジをカットしてグラフをサブグラフに分割することで、ノードのクラスタリングを行なう。たとえば図4のグラフの場合、v4-V6間のエッジをカットすることで、G1={v1,…,v5}、G2={v6,v7,v8}の2つの部分グラフに分割でき、V={v1,…,v8}がV1={v1,…,v5}とV2={v6,v7,v8}とにクラスタリングされる。スペクトラル・クラスタリングにはいくつかのアルゴリズムがあるが、ここではグラフを2つのサブグラフに分割する方法であるMax-MinCut(Mcut)法について述べる。

【0040】
いま、重み付きグラフG(V,E)とその隣接行列Wが与えられたとし、グラフGを2つのサブグラフA,Bに分割する問題を考える。サブグラフA,B間の類似度cut(A,B)を次の式(4)により定義する。

【0041】
【数4】
JP0005457737B2_000005t.gif
ただし

【0042】
【数5】
JP0005457737B2_000006t.gif
である。すなわち、サブグラフ間の類似度cut(A,B)とは、A,B間のエッジの重みの総和である。Mcut法はcut(A,B)を最小、かつW(A)及びW(B)を最大化するようなサブグラフA,Bを探索する方法である。Mcut法の目的関数Jは次式で定義される。

【0043】
【数6】
JP0005457737B2_000007t.gif
グラフG(V,E)のノードのインデックスの与え方は任意であるから、隣接行列Wを以下のように定義しても一般性は失わない。

【0044】
【数7】
JP0005457737B2_000008t.gif
ここで、WA、WBはサブグラフA,B内の隣接行列であり、WA,BはA,B間の隣接行列である。隣接行列の定義より、WT=WB,Aである。さらに、サブクラスへの分割を表現するベクトルx=[1,…,1,0,…,0]、y=[0,…,0,1,…,1]が次式を満たすようにこれらを定義する。

【0045】
【数8】
JP0005457737B2_000009t.gif
ただし、D=diag(We)、e=[1,…,1]Tである。これらの式を用いて、目的関数(6)式を書き換えると

【0046】
【数9】
JP0005457737B2_000010t.gif
が得られる。また、xTDy=0、xTWx>0、yTWy>0であることがわかる。Jの各項はx,yについて対称であるから、以降、第1項のみについて考える。

【0047】
目的関数(11)式の目的変数x,yはバイナリ変数であるから、この最小化問題はNP困難である。そこで、(11)式のバイナリ変数を緩和する。ここでは、ノードのインデックスに対応する指標ベクトルqを導入する。qの、ノードuに対応する要素をqu={a,-b}(a,b>0)とする。qu=aのときu∈A、qu=-bのときu∈Bである。目的関数は次の(12)式のように書き直される。

【0048】
【数10】
JP0005457737B2_000011t.gif
さらに(12)式を次のようにスケーリングする。

【0049】
【数11】
JP0005457737B2_000012t.gif
ここで、W=D-1/2WD-1/2q=D1/2q/|D|1/2である(各変数の前の「」は、式中では変数の上に記載されている。)。また、qT~Wq>0である。Wの要素(W)i,jは1≧(W)i,j≧0、各列の和は1であるから、~Wは確率推移行列とみなすことができる。

【0050】
さらに目的関数Jqにおいて、P≡I-WとおくとJqはPのRayleigh商で表現できることがわかる。Rayleighの原理によると、行列Aの商R(x)は最小固有値λ1に対応する固有ベクトルx1によって最小となり、その最小値はλ1である(本明細書では、固有値はλ1≦…≦λNのように並べるものとする。)。よって、緩和した最小化問題はPの固有値問題に帰着することがわかる。Pのある固有ベクトルをD1/2eと仮定する。Wの第i行をwiとすると、(D)i,i=di=wieであるから

【0051】
【数12】
JP0005457737B2_000013t.gif
であり、確かにPのある固有ベクトルはD1/2eで、その固有値は0であることが示される。ところで、Wは確率推移行列とみなすことができたから、その最大固有値は1である。よって、xTPx=xT(I-W)x≧0であるから、Pは半正定値行列である。これより、固有ベクトルD1/2eはPの最小固有値λ=0に対応する固有ベクトルz1であることがわかる。

【0052】
したがって、Rayleighの原理より最小固有値に対応する固有ベクトルz1が(13)式を最小とする解である。しかし、この解は全ての要素が正であり、求めるqの要素がqu={a,-b}であるという仮定に反する。したがってこの解は採用できない。次に小さな解はRayleigh商のmax-min定理より、2番目に小さな固有値λ2である。λ2に対応する固有ベクトルz2は、Pの対称性よりz1Tz2=0であるから、必ず正負の要素を持つことがわかる。すなわち、この解は最初の仮定に反しない。したがって、λ2に対応する固有ベクトルz2が求める最適なqである。

【0053】
以上の考察より、最終的な定式化は
【数13】
JP0005457737B2_000014t.gif
において、固有値問題
【数14】
JP0005457737B2_000015t.gif
の2番目の大きさの固有値λ2に対応する固有ベクトルz2を選択すればよいことになり、求めるqはq=D-1/2z2で表される。なお、このλ2、z2をそれぞれFielder値、Fielderベクトルと呼ぶ。

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

【0055】
理論的な考察により、緩和問題の目的関数が最小となるとき、(6)式においてほぼW(A)=W(B)が成立することが示されている。すなわち、Mcut法ではサブグラフ内の類似度がバランスすると期待できる。他のスペクトラル・クラスタリングのアルゴリズムでは必ずしもサブグラフ内の類似度がバランスしないため、このことがMcut法の利点とされる。

【0056】
1.3 その他のアルゴリズム
前節で紹介したMcut法に類似の手法として、Normalized Cut(Ncut)法がある。Ncut法の評価関数は次の(17)式で表される。

【0057】
【数15】
JP0005457737B2_000016t.gif
ここで、V=A∪B、deg(A,V)=Σu∈A,v∈VWu,vであり、deg(B,V)も同様に定義される。解析の詳細はMcut法とほぼ同等である。

【0058】
Mcut/Ncut法では類似度の与え方は任意であり、類似度行列が無向グラフの隣接行列としての性質を満たしていればよい。

【0059】
一方、類似度行列の構成にGaussianカーネルを用いることが従来の一手法として提案されている。この手法では、サンプルsi,sjの類似度(W)i,jとして

【0060】
【数16】
JP0005457737B2_000017t.gif
を用いる。ここでd(・,・)は距離関数、σは調整パラメータである。この手法に関連してさらに、グラフの2分割を繰返すMcut法と異なり、スケーリング後の類似度行列の固有空間にてサンプルの集合を直接任意の個数のクラスタに分割する手法が提案されている。この場合、カーネルを用いて非線形写像を行なった後に固有値問題を解いていることになるから、カーネル主成分分析(主成分分析を以下「PCA」と呼ぶ。)と類似の手法であるとされる。この手法をNg-Jordan-Weiss(NJW)アルゴリズムと呼ぶ。図3に示したのはNJWアルゴリズムを用いたスペクトル・クラスタリングの結果である。

【0061】
NJWアルゴリズムの(18)式のσの決め方において、ローカルスケーリングという方法が提案されている。(W)i,jは次のように書き換えられる。

【0062】
【数17】
JP0005457737B2_000018t.gif
ここで、σi=d(Si,SK)であり、SKはSiのK番目に近傍のサンプルを示す。Kは調整パラメータであるが、数値実験ではK=7で概ね良好な結果が得られたと報告されている。

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

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

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

【0066】
いま、サンプルx1に着目する。まず、サンプルx1が原点となるように空間全体を平行移動させる。つまり、全サンプルxi(i=1,…,7)からx1を減算する。この操作によってアフィン部分空間Pは原点を含むことになるため線形部分空間の定義を満たす。この線形部分空間をVとする。次に図6に示すように、任意のサンプルと原点を結ぶ直線を引き、この直線上で別のサンプルが発見できたとする。図6に示す例では、サンプルx2及びx5並びにサンプルx3及びx4がこのような関係を満たしている。このとき、これらのサンプルのペアの相関係数の絶対値は1である。このような関係が成り立つとき、「x2~x5」のように書く。一方、部分空間Vの要素ではないサンプルx6及びx7の相関係数の絶対値は1未満である。それゆえ、相関係数が±1であるペアのサンプルは、同一の相関関係を有していると判定できる。

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

【0068】
ここで、各サンプルをノードとし、サンプル間が同一の相関関係を有しているか否かをそれらのサンプルに対応するノードが隣接しているかとして表現した無向グラフ及びその隣接行列S1∈R7×7を考える。この例では、x2~x5及びx3~x4であるから、隣接行列S1は、(S1)2,5=(S1)5,2=1及び(S1)3,4=(S1)4,3=1とし、それ以外の要素を0とする行列として表現できる。

【0069】
次に、サンプルx2が空間の原点と一致するように平面を平行移動させて同様の手順を行ない、隣接行列S2を求める。これを全てのサンプルについて繰返し、最終的にS1,…,S7を求める。

【0070】
この手続きによって隣接行列が7つ得られる。これら7つのグラフは、同一のノードの集合を持ち、エッジの集合の部分集合が異なる7つの無向グラフと考えることができる。すなわち、G1(V,E(G1)),…,G7(V,E(G7))である。ここで、これらのグラフの和のグラフG′=G1∪…∪G7を考えると、グラフG′のエッジの集合はグラフG1からG7のエッジの和集合E(G1)∪…∪E(G7)であることがわかる。これは隣接行列Si(i=1,…,7)のブール和(各行列の対応する要素の論理和を要素とする行列)を求めることに等しいから、グラフG′の隣接行列はS′=S1+…+S7となる。ただし、ここで+は行列のブール和を表す。以上の手続きより、サンプル間の相関関係を表現したグラフを得ることができる。

【0071】
しかし、サンプルの分布と原点の選び方とによっては、異なった相関関係を持つサンプルのペアであっても、偶々同一の相関関係を有すると判断されてしまうことがある。たとえば図5及び図6に示した例について考える。図7のようにサンプルx3が原点と一致するように平面を平行移動させると、x1~x3以外にx6~x7と判定されてしまう(サンプルx5と原点とを結んだ直線上にサンプルx7が存在する。)。つまり、グラフG′では真に同一の相関関係を有するサンプル間のエッジと、偶々同一の相関関係を有すると判断されてしまったサンプル間のエッジとが区別できない。

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

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

【0074】
図8を参照して、本実施の形態に係るNC-スペクトル・クラスタリングを実現するコンピュータプログラムの制御構造について説明する。いま、サンプルxn∈RM(n=1,…,N)がデータベースに保存されており、これより変数間の相関関係を指標として、スペクトラル・クラスタリングのための類似度行列を求めたいとする。

【0075】
図8を参照して、このプログラムは、以下の繰返し処理のため繰返し制御変数Lに1を代入するステップ30と、類似度行列Sのための2次元配列S[N,N]の領域を確保するステップ32と、配列S[N,N]の全ての要素をゼロでクリアするステップ34とを含む。

【0076】
このプログラムはさらに、L番目(1≦L≦N)のサンプルについて、全サンプルとの間の隣接行列SLを保存するための2次元配列SL[N,N]の領域を確保するステップ36と、配列SL[N,N]の各要素をゼロでクリアするステップ38と、L≦n≦Nを満たす全てのnについてサンプルxnからサンプルxLを減算して新たなサンプルxnを算出するステップ42と、ステップ42で算出されたサンプルxnに関し、全てのxkとxlとの組合せ(1≦k、l≦N、ただしk≠l)について、相関係数Ck,lを算出するステップ46と、ステップ48で算出された相関係数Ck,lのうちで|Ck,l|≧γ(γは予め定めるしきい値)を満たすk及びlを求めるステップ48と、SL[k、l]とSL[l、k]とに1を代入するステップ50と、配列Sの各要素に、配列SLの対応する要素の値を加算するステップ52とを含む。

【0077】
このプログラムはさらに、変数Lの値がN-1より小さいか否かを判定し、判定結果に応じて制御の流れを分岐させるステップ54と、ステップ54における判定結果がYESのときに実行され、変数Lに1を加算して制御をステップ36に戻すステップ56と、ステップ54における判定結果がNOのときに実行され、類似度行列を表す配列として配列Sを出力し、処理を終了するステップ58とを含む。

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

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

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

【0081】
【数19】
JP0005457737B2_000020t.gif
ここでS~N(0,10)であり、ai∈R2は係数行列である。また、n=[n1 n2T、ni~N(0,0.1)である。ただし、N(m,σ)は平均m、標準偏差σの正規分布にしたがう乱数を表す。係数ベクトルはそれぞれ、a1=[1 2]T、a2=[2 2]T、a3=[2 1]Tとし、各クラスにおいて100サンプル生成した。

【0082】
本シミュレーション例において、本実施の形態のステップ48で用いるパラメータ(しきい値)γは、0.999とした。生成したサンプル及びk-平均法、従来のスペクトラル・クラスタリングでのクラスタリング結果が、それぞれ図1~図3に示したものである。本願実施の形態にしたがったクラスタリング結果を図9に示す。

【0083】
図1~図3と、図9とを参照して、k-平均法でクラスタリングした結果(図2)及び従来のスペクトル・クラスタリングでクラスタリングした結果(図3)は、いずれも図1と比較すると正しいクラスタリングとはいえないことがわかる。一方、本実施の形態によるクラスタリング結果(図9)では、原点付近のサンプルの属するクラスは正しく識別できていないが、それ以外の領域では正しく識別できていることがわかる。

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

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

【0086】
【数21】
JP0005457737B2_000022t.gif

ここでAi∈R2×5は係数行列であり、S=[s1 s2]、Si~N(0,10)である。また、n=[n1 n2 n3 n4 n5T、ni~N(0,0.1)である。係数行列Ai

【0087】
【数22】
JP0005457737B2_000023t.gif


とし、各クラスについて100サンプル生成した。すなわち、K=100である。本シミュレーション例において、本実施の形態のステップ48で用いるパラメータ(しきい値)γは、0.999とした。

【0088】
生成したサンプルに対し、k-平均法(KM)、スペクトラル・クラスタリング(SC)、及び本実施の形態に係るNC-スペクトラル・クラスタリング(NC)のクラスタリングにより各クラスについてどの程度の精度で識別できたかをテーブル1に示す。テーブル1を参照して、NC-スペクトル・クラスタリングによれば、k-平均法及び従来のスペクトル・クラスタリングのいずれと比較しても、どのクラスについてもより高い精度でクラスを識別できていることが分かる。以上の結果より、本実施の形態に係るNC-スペクトラル・クラスタリングの有効性が示されたといえる。

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

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

【0091】
【数23】
JP0005457737B2_000025t.gif
(ただし、A,Bは原料であり、Cは製品、Dは副製品)を行なう5台の反応炉82,84,…,86と、これら反応炉、及びこれら反応路のジャケットの温度をそれぞれ測定するための図示しない複数個のセンサと、これら複数個のセンサの出力にしたがって、所定のプロセス制御を反応炉ごとに行なう複数のコントローラ92,94,…,96と、これら反応炉で過去に行なわれたバッチプロセスの結果から得られた生成物に関する情報、及び各プロセスを開始するときの各反応炉に関する所定のプロセスパラメータ(ここではジャケット温度Tjと反応炉内温度Trとする。)に基づき、新たなバッチプロセスで得られる製品Cの分子量MCを推定するための、コンピュータからなる製品分子量推定装置98とを含む。本プラントでは、目的生成物は製品Cである。プロセスの目的は副製品Dの生成を抑制し、できるだけ多くの製品Cを得ることである。バッチプロセス運転開始前に原料A及びBが室温で反応炉に仕込まれる。

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

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

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

【0095】
図11を参照して、製品分子量推定装置98は、コンピュータ120と、全てコンピュータ120に接続された、モニタ122と、キーボード126と、マウス128と、プリンタ124と、を含む。さらに、コンピュータ120はDVD-ROM(Digital Versatile Disk Read-Only-Memory:ディジタル多用途ディスク読出専用メモリ)ドライブ150と、半導体メモリポート152とを含む。

【0096】
図11を参照して、コンピュータ120はさらに、DVD-ROMドライブ150と半導体メモリポート152とに接続されたバス142と、全てバス142に接続された、CPU140と、コンピュータ120のブートアッププログラムを記憶するROM144と、CPU140によって使用される作業領域を提供するとともにCPU140によって実行されるプログラムのための記憶領域となるRAM146と、音声データ、音響モデル、言語モデル、レキシコン、及びマッピングテーブルを記憶するためのハードディスクドライブ148と、ネットワーク164への接続を提供するネットワークインターフェイス154とを含む。

【0097】
上述の実施の形態に係る製品分子量推定装置98を実現するソフトウェアは、DVD-ROM162又は半導体メモリ160等の媒体に記録されたオブジェクトコード又はスクリプトの形で流通し、DVD-ROMドライブ150又は半導体メモリポート152等の読出装置を介してコンピュータ120に提供され、ハードディスクドライブ148に記憶される。CPU140がプログラムを実行する際には、プログラムはハードディスクドライブ148から読出されてRAM146にロードされる。図示しないプログラムカウンタによって指定されたアドレスから命令がフェッチされ、その命令が実行される。CPU140はハードディスクドライブ148から処理すべきデータを読出し、処理の結果をこれもまたハードディスクドライブ148に記憶する。推定された製品分子量のハードコピーが、プリンタ124により出力される。

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

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

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

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

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

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

【0104】
このプログラムはさらに、ステップ182に続き、ステップ182により得られた実測データに対して、図8を参照して説明したNC-スペクトラル・クラスタリングのためのプログラムを適用することにより、実測データを複数のクラスに分類するステップ184と、ステップ184により得られたクラスごとに、予め選択された構築方法により統計的モデルを構築するステップ186と、ステップ186で得られた各クラスの統計的モデルをHDD148に保存して処理を終了するステップ188とを含む。

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

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

【0107】
4.3 シミュレーション
本実施の形態について、以下のような条件でその性能についての数値シミュレーションを行ない、従来の方法による結果と比較した。以下の数値シミュレーションでは、バッチプロセス運転開始前に原料A及びBが室温20℃で反応炉に仕込まれるものとする。このシミュレーションでは、原料A及びBの仕込み量はそれぞれN(20,0.1)にしたがう乱数として変化するものとする。

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

【0109】
本プロセスにおいて運転中に測定されるのは,ジャケット温度Tjと反応炉内温度Trであり、それぞれサンプリング周期は1分である。バッチ終了時間は120分で固定する。あるバッチの運転データを図14に示す。

【0110】
本シミュレーションでは5台の反応炉82,84,…,86が並列に運転されている。バッチ終了時にはこれら反応炉のメンテナンスが行なわれるため、バッチ毎に反応炉の伝熱係数が変化する。反応炉それぞれの伝熱係数Ui(i=1,…,5)は、U1=U2=U(40.60,40.62),U3=U4=U(40.57,40.59),U5=U(40.54,40.56)として。ただし、U(a,b)は閉区間[a,b]上の一様乱数であることを示す。すなわち,本プロセスでは5台の反応炉が存在するが,それらの個体差としては3種である。本シミュレーションでは、各反応炉についてプロセス20バッチ、計100バッチのデータがDBに保存されているとし、これらのデータより製品の分子量推定を行なった。

【0111】
4.4 クラスタリング
データベースに保存されている100バッチ分の運転データのクラスタリングを,従来法であるk-平均法と、本実施の形態に係るNC-スペクトラル・クラスタリングを用いて実施した。シミュレーションでは、図12のステップ182から処理を開始した。すなわち、100データに対してPCAを行なって各データの次元を圧縮した。ここで採用する主成分数は2、主成分スコアt1,t2を入力としてクラスタリングを実施した。クラス数は3とし,NC-スペクトラル・クラスタリングにおけるしきい値γ=0.99とした。

【0112】
クラスタリング結果を図15及び図16に示す。ここではt1-t2平面上でのサンプルの分布を示している。k-平均法ではさらに、各クラスの中心を○印で示している。この結果より,k-平均法ではクラス中心からの距離に基づいてサンプルを分類しているが、本実施の形態では各クラスは相関関係の違いを表現していることがわかる。

【0113】
4.5 ソフトセンサ設計
製品Cの分子量MC[kmol]を推定するソフトセンサを設計する。設計手順は図12及び図13に示したとおりであり、図13のステップ204で新たなプロセスでの製品の分子量を推定した。なお、対比のために、k-平均法を用いたクラスタリングも行ない、その結果を用いて図13に示す方法で製品の分子量を推定した。クラスタリング時のクラス数は3である。

【0114】
以下のように、3クラスそれぞれで統計的モデルを構築した。ここで、統計的モデルの入力は反応炉温度Tr(k)及びジャケット温度Tj(k)(k=0,…,120)、出力は製品Cの分子量MCである。モデル構築にはPLS(Partial Least Square)回帰を用いた。採用した潜在変数の数はともに3である。新たに測定されたサンプルの識別には、k-平均法ではクラス中心からの距離、本実施の形態では前述の通り各クラスとの間で求めたQ統計量を指標として用いた。

【0115】
本シミュレーションでは、新たな反応炉R6のMCを推定した。R6の伝熱係数はU6=U(40.60,40.62)とした。検証には20サンプルを用い、k-平均法及び本実施の形態それぞれを用いたソフトセンサの推定性能を比較した。

【0116】
結果を図17(A)(B)にそれぞれ示す。図中、グラフの横軸が真値、縦軸が予測値を示している。「RMSE」は根平均2乗誤差、「R」は真値と予測値との相関係数である。この結果より、k-平均法を用いた場合は精度の良い推定ができていないが、本実施の形態を用いた場合は高い推定性能を達成していることが分かる。たとえば、k-平均法を用いた場合のRMSEが0.007であるのと比較して、本実施の形態ではRMSEは0.004となり、RMSEにして43%改善した。このことより、本実施の形態によれば、k-平均法を用いた場合よりも装置の個体差を正しくクラスタリングできていることがわかる。

【0117】
5. 第2の実施の形態
5.1 構成
第1の実施の形態と同様のシステムハードウェア構成で、プロセス異常を検出する実施の形態について、説明する。

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

【0119】
図18を参照して、本実施の形態では、上記したNC-スペクトラル・クラスタリングと多変量統計的プロセス管理手法(Multivariate Statistical Process Control;MSPC)を組合せた異常検出を行なう。なお、比較のためにNC-スペクトラル・クラスタリングに代えてk-平均法を用いた場合の結果についても後に示す。以下、k-平均法又はNC-スペクトル・クラスタリングと、MSPCとを組合せた異常検出法をそれぞれKM-MSPC、NC-MSPCと呼ぶことにする。

【0120】
図18を参照して、本実施の形態に係るプロセス異常の検出装置を実現するためのプログラムの制御構造は以下のとおりである。すなわち、このプログラムは、これまでにDBに記憶されていた全サンプルに対してNC-スペクトル・クラスタリングを用いたクラスタリングを実施するステップ230と、ステップ230に続き、全クラスについてT2統計量とQ統計量との管理限界T2SQS(S=1,2,…,S)(各記号の前の「」は、通常は各記号の直上に記載されるべきものである。)とを設定するステップ232とを含む。ただし、Sはクラス数である。

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

【0122】
図19を参照して、このプログラムは、所定のインターバル(たとえば10分)が経過するまで待機するステップ250と、10分が経過するごとに、対象となる反応炉の反応炉温度、ジャケット温度等の新たな測定値xqをサンプリングするステップ252と、ステップ252でサンプリングされた測定値xqと、新たなサンプルxqと、全てのクラスとの間でT2統計量とQ統計量Tq2、Qqとを求めるステップ254と、全てのクラスに対して第1の判定式Tq2T2Sが成立しているか否かを判定し、判定結果に応じて制御の流れを分岐させるステップ256と、ステップ256で第1の判定式が成立していないと判定されたときに、全てのクラスに対して第2の判定式QqQSが成立しているか否かを判定し、判定結果に応じて制御の流れを分岐させるステップ258とを含む。

【0123】
ステップ256又は258における判定結果がYESのとき、すなわち新たにサンプリングされた測定値xqがいずれのクラスにも属していないと判定されたときには、ステップ262でプロセス異常が発生したと判定され、必要な処理が実行される。ステップ256及び258の双方において判定結果がNOのときには、ステップ260でプロセスは正常であると判定される。この後、制御の流れはステップ250に戻り、さらに10分が経過するまで待機した後、上記した処理を再び実行する。

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

【0125】
本ケーススタディでは、対象とする異常として、ジャケット温度測定センサのドリフトを考慮した。すなわち時間と共にセンサのゼロ点が移動し、センサにて測定される測定値はTj(k)=Tj0(k)+δT×(k-k0)のように変化する。ここでTj0はゼロ点が固定されているときの測定値、kは測定時刻、k0はドリフト開始時間、δTは単位時間当たりのドリフト量(K/分)である。
いま、新たな反応炉R6を対象にドリフト開始時間k0=20分とし、δTを変化させ異常を発生させた。対象とするデータは後掲のテーブルに示す6ケースであるが、このうちケース1、ケース2は正常データ(δT=0)である。本実施の形態に係る装置の目的は、異常発生後できるだけ早く、正確に異常を検出することである。なお、T2統計量、Q統計量の管理限界はそれぞれ信頼区間95%として設定し、採用した主成分数は5である。

【0126】
KM-MSPC及び提案するNC-MSPCを用いた異常検出結果を表2-表17に示す。表では、11分~20分(第2インターバル)、21分~30分(第3インターバル)、31分~40分(第4インターバル)、41分~50分(第5インターバル)におけるそれぞれの管理限界(CL)と各ケースでのT2統計量、Q統計量を示す。

【0127】
この結果より、KM-MSPC、NC-MSPC双方ともに、未だ異常の発生していない第2インターバルの全てのケースと、正常データであるケース1、ケース2とについては、その後も正しく正常であると判定している。しかし、KM-MSPCでは第3インターバル、第4インターバルでケース4の異常を検出できていない。一方、NC-MSPCでは第3インターバルでケース4の異常を正しく検出できている。以上より、NC-MSPCを採用した本実施の形態がKM-MSPCより素早く異常を検出することができることがわかる。

【0128】
【表2】
JP0005457737B2_000026t.gif

【0129】
【表3】
JP0005457737B2_000027t.gif

【0130】
【表4】
JP0005457737B2_000028t.gif

【0131】
【表5】
JP0005457737B2_000029t.gif

【0132】
【表6】
JP0005457737B2_000030t.gif

【0133】
【表7】
JP0005457737B2_000031t.gif

【0134】
【表8】
JP0005457737B2_000032t.gif

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

【0136】
Appendix A
バッチプロセスモデル
シミュレーションにて用いたバッチプロセスのモデルを示す。プロセスにおける収支式等は以下の通りである。

【0137】
【数24】
JP0005457737B2_000034t.gif

【0138】
【数25】
JP0005457737B2_000035t.gif

【0139】
【数26】
JP0005457737B2_000036t.gif
ここで、a(k)~N(0,0.1)であり、kは測定ステップ数を表す。ただし、N(m,σ)は平均0標準偏差σの正規分布にしたがう乱数である。温度コントローラからジャケット温度への動特性は、(45)式で与えられる。

【0140】
【数27】
JP0005457737B2_000037t.gif
K=1/16、τ=2分である。モデルのパラメータをテーブル18に、初期状態をテーブル19に、それぞれ示す。

【0141】
【表10】
JP0005457737B2_000038t.gif
【符号の説明】
【0142】
80 プラント
82,84,86 反応炉
92,94,96 コントローラ
98 製品分子量推定装置
140 CPU
Drawing
(In Japanese)【図8】
0
(In Japanese)【図10】
1
(In Japanese)【図19】
2
(In Japanese)【図1】
3
(In Japanese)【図2】
4
(In Japanese)【図3】
5
(In Japanese)【図4】
6
(In Japanese)【図5】
7
(In Japanese)【図6】
8
(In Japanese)【図7】
9
(In Japanese)【図9】
10
(In Japanese)【図11】
11
(In Japanese)【図12】
12
(In Japanese)【図13】
13
(In Japanese)【図14】
14
(In Japanese)【図15】
15
(In Japanese)【図16】
16
(In Japanese)【図17】
17
(In Japanese)【図18】
18