JP6834295B2

JP6834295B2 - 状態予測装置

Info

Publication number: JP6834295B2
Application number: JP2016187019A
Authority: JP
Inventors: 謙一 ▲濱▼口; 周平江本
Original assignee: IHI Corp
Current assignee: IHI Corp
Priority date: 2016-09-26
Filing date: 2016-09-26
Publication date: 2021-02-24
Anticipated expiration: 2036-09-26
Also published as: JP2018055169A

Description

本発明は、状態予測装置に関する。

熱交換等を行うプラントでは、プラントの動きを再現するモデルを用いて、プラントにおけるセンサが設置されていない場所の状態量を推定し、さらに将来の状態量を予測することが行われている。このようなプラントの将来の状態量を予測する状態予測装置が、例えば、特許文献１に記載されている。

特許文献１に記載された状態予測装置では、実際のプラントと同期して動作するシミュレーションモデルを用いて、将来の状態量を予測している。また、この予測装置では、実際のプラントの測定結果に基づいて、シミュレーションモデルの係数パラメータを補正することで、シミュレーションモデルの出力と実際のプラントの測定値とを合致させている。

特許４７８９２７７号

このような状態予測装置の分野においては、実際のプラントの動きを再現するモデルを実際のプラントに精度良く一致させ、実際のプラントの将来の状態量をより精度良く予測することが求められている。しかし、一般的なプラントは機械系などと比べてモデルが複雑であり、全ての反応を再現できる完全なモデルを得るのは難しい。また、装置の個体差や、触媒の劣化などの経年変化によっても、機器の特性は変化する。このように、物理又は化学の理論に基づく動特性モデルでは表現できない機器特性の違いが、予測精度を低下させる。

本発明は、対象装置の将来の状態量をより精度良く予測することができる状態予測装置を提供する。

本発明の一側面に係る状態予測装置は、状態量の予測対象とする対象装置の動特性モデル、及び対象装置の機器特性を表す係数パラメータに関する係数パラメータ変動モデルを用いてカルマンフィルタを設計し、設計されたアルゴリズムに対して対象装置の観測値を入力してカルマンフィルタ処理を行うことで、対象装置の状態量及び係数パラメータを推定する推定部と、推定された状態量、推定された係数パラメータ、動特性モデル、及び係数パラメータ変動モデルに基づいて、推定部において状態量の推定の対象とした時点に対して時間的に後の時点における対象装置の状態量を予測する予測部と、を備える。

この状態予測装置において推定部は、対象装置の動特性モデル等を用いてカルマンフィルタを設計し、カルマンフィルタ処理を行うことで、対象装置の状態量及び係数パラメータを推定する。このように、推定部は、カルマンフィルタを用いて推定を行うことにより、対象装置の状態量に加え、係数パラメータについても精度良く推定することができる。予測部は、カルマンフィルタを用いて推定された対象装置の状態量及び係数パラメータ等を用いて、対象装置の状態量を予測する。このように、予測部は、カルマンフィルタを用いて精度良く推定された対象装置の状態量及び係数パラメータを用いて予測を行うことで、対象装置の将来の状態量をより精度良く予測することができる。

推定部は、状態量及び係数パラメータを推定するとともに、推定された状態量と実際の状態量との推定誤差の大きさを推定し、予測部は、推定誤差の大きさに基づいて、予測された状態量の予測幅を算出してもよい。この場合、対象装置のオペレータ等は、予測された対象装置の将来の状態量に加えて、予測された状態量の予測幅を把握することができる。例えば、対象装置の状態量が予め定められた範囲を超えないように対象装置の制御を行う際に、予測幅が大きい場合と小さい場合とでは、例えばオペレータによる対象装置の制御の実施タイミング、及び制御アルゴリズムにおける制御量（操作量、バルブ開閉の指令値など）の大きさ等が異なることがある。このような場合であっても、予測された対象装置の将来の状態量に加えて、状態量の予測幅も算出することで、対象装置のオペレータ等は予測幅に応じて対象装置の制御を行うことができる。

本発明によれば、対象装置の将来の状態量をより精度良く予測することができる。

実施形態に係る状態予測装置及びプラントシステムの概略構成を示す図である。図２は、プラントとして用いる熱交換型の流通反応器を模式的に示す図である。予測された状態量の予測値と予測幅とを示す図である。予測幅の算出方法の一例を示す図である。予測幅の算出方法の一例を示す図である。図１の推定部が行う状態量の推定処理の流れを示すフローチャートである。図１の予測部が行う将来の状態量の予測処理の流れを示すフローチャートである。

以下、本発明の実施形態について、図面を参照しながら説明する。なお、図面の説明において同一要素には同一符号を付し、重複する説明を省略する。

図１は、実施形態に係る状態予測装置の概略構成を示す図である。図１に示す状態予測装置２は、プラントシステム１のプラント（対象装置）１１の現在の状態量をカルマンフィルタを用いて推定し、推定した状態量に基づいて将来の状態量を予測する装置である。なお、明細書において、プラント１１の状態量を「推定」することとは、プラント１１の「現在の状態量、又は過去のある時点の状態量を推定」することをいう。すなわち、既に起こっている状態量を推定することである。また、プラント１１の状態量を「予測」することとは、プラント１１の「現在から将来の状態量、又は過去のある時点から将来の状態量を予測」することをいう。すなわち、現在又は過去のある時点から将来の未知の状態量を予測することである。

プラントシステム１と状態予測装置２とは、ネットワークＮＷを介して互いにデータの送受信を行うことができる。状態予測装置２は、プラントシステム１とは異なる場所に設置されている。

プラントシステム１は、プラント１１、観測部１２、操作量取得部１３、表示部１４、及び通信部１５を備えている。

プラント１１は、状態予測装置２による将来の状態量の予測等の予測対象となる対象装置である。本実施形態では、プラント１１として、図２に示す熱交換型の流通反応器を用いる。図２に示すように、プラント１１は、原料ガスが流通する触媒層１０１、及び加熱ガス（加熱ガスに代えて加熱流体でもよい）が流通する加熱路１０２を備えている。加熱路１０２を流れる加熱ガスの熱は、触媒層１０１と加熱路１０２との間の壁１０３を介して触媒層１０１に伝達される。原料ガスは、触媒層１０１を通過することによって化学的に反応し、原料ガスから反応ガスに変化する。なお、本実施形態では、化学反応が吸熱反応の場合の反応器を例示したが、化学反応が発熱反応の場合には加熱ガスが流通する加熱路１０２に代えて、冷却水等が流れる冷却路を用いる。

このプラント１１において、加熱ガスから反応ガスへの単位面積当たりの伝熱量ｑ［Ｗ／ｍ^２］は、ｑ＝Ｋ（Ｔ_ｈ−Ｔ_ｃ）で与えられる。ただし、Ｔ_ｈは加熱ガスの代表温度［Ｋ］、Ｔ_ｃは反応ガスの代表温度［Ｋ］である。また、Ｋは、加熱路１０２と触媒層１０１との間の総括伝熱係数［Ｗ／（ｍ^２Ｋ）］である。

図１に戻り、観測部１２は、プラント１１の種々の状態量を観測する。観測部１２は、プラント１１の各部に取り付けられた複数のセンサを備えている。観測部１２は、プラント１１の状態量として、例えば、プラント１１内部の温度、ガスの流量、及びガスの密度等を観測する。なお、以下では、観測部１２で観測されたプラント１１の状態量を「観測値」という。

操作量取得部１３は、プラント１１の運転時における各種の操作量を取得する。操作量取得部１３は、プラント１１の操作量として、例えば、流通させるガスの温度、濃度、及び圧力等を取得する。操作量取得部１３は、プラント１１の操作量として、例えば、プラントシステム１を操作するオペレータの操作量（オペレータが設定した設定値）、プラント１１の制御を行う制御装置が設定した設定値等を取得する。

表示部１４は、プラントシステム１のオペレータ等に情報を提示する機器である。表示部１４は、状態予測装置２によって予測されたプラント１１の将来の状態量等を表示する。また、表示部１４は、観測部１２によって観測されたプラント１１の観測値を表示してもよい。

通信部１５は、ネットワークＮＷを介して状態予測装置２の通信部２１と通信を行う。通信部１５は、観測部１２で観測されたプラント１１の観測値、及び操作量取得部１３で取得された操作量を状態予測装置２へ送信する。また、通信部１５は、状態予測装置２で予測されたプラント１１の将来の状態量等を状態予測装置２から受信する。

状態予測装置２は、通信部２１、演算部２２、及び表示部２３を備えている。通信部２１は、ネットワークＮＷを介してプラントシステム１の通信部１５と通信を行う。通信部２１は、演算部２２で予測されたプラント１１の将来の状態量等をプラントシステム１へ送信する。また、通信部２１は、プラント１１の観測値及びプラント１１の操作量をプラントシステム１から受信する。

表示部２３は、状態予測装置２のオペレータ等に情報を提示する機器である。表示部２３は、演算部２２で予測されたプラント１１の将来の状態量等を表示する。

演算部２２は、プラント１１の将来の状態量等の演算処理を行う。演算部２２は、物理的には、１又は複数のプロセッサ、主記憶装置であるＲＡＭ（Random Access Memory）及びＲＯＭ（Read Only Memory)等の記憶装置、並びに、ハードディスク装置等の補助記憶装置等のハードウェアを備えるコンピュータとして構成され得る。図１に示される演算部２２の各機能は、記憶装置等のハードウェアに１又は複数の所定のコンピュータプログラムを読み込ませることにより、１又は複数のプロセッサの制御のもとで各ハードウェアを動作させるとともに、記憶装置及び補助記憶装置におけるデータの読み出し及び書き込みを行うことで実現される。

演算部２２は、機能的には、推定部２２Ａ、及び予測部２２Ｂを備えている。推定部２２Ａは、プラント１１の動特性モデル、プラント１１の機器特性を表す係数パラメータに関する係数パラメータ変動モデル、及びノイズ情報を用いてカルマンフィルタを設計する。推定部２２Ａは、設計したアルゴリズムに対してプラント１１の観測値及び操作量を入力してカルマンフィルタ処理を行うことで、プラント１１の状態量及び係数パラメータを推定する。

まず、動特性モデル、係数パラメータ変動モデル、及びノイズ情報の詳細についてそれぞれ説明する。

動特性モデルとは、プラント１１の状態量の変動を数式で表した数式モデルである。動特性モデルは、次の式（１）及び式（２）で表現される。

ここで、式（１）は、状態方程式である。式（２）は、観測方程式である。
式（１）及び式（２）における記号の意味は、次のとおりである。ｆ、ｈは、線形又は非線形の関数である。ｆ、ｈは、それぞれ状態遷移関数及び観測関数である。ｘ_（ｎ）は、離散時刻ｎにおけるプラント１１の状態量（温度、流量、密度など）である。ｙ_（ｎ）は、離散時刻ｎにおけるプラント１１の観測値（温度、流量、ガス組成など）である。ｕ_（ｎ）は、離散時刻ｎにおけるプラント１１の操作量（投入物質の温度や濃度、圧力など）である。θ_（ｎ）は、離散時刻ｎにおける機器特性を表す係数パラメータである。なお、係数パラメータは、時間に依存していなくてもよい。

ｖ_ｘ（ｎ）は、離散時刻ｎにおけるプロセスノイズ（操作量に対するノイズなど）である。ｗ_（ｎ）は、離散時刻ｎにおける観測ノイズ（観測部１２のセンサに対するノイズ）である。Ｂ_ｘは、プロセスノイズｖ_ｘ（ｎ）の係数行列である。

なお、上記の式（１）及び式（２）で表わされるプラント１１の動特性モデルは、既知の種々の方法によって生成することができる。

係数パラメータ変動モデルとは、プラント１１の機器特性を表す係数パラメータの変動を表す数式モデルである。なお、ここで機器特性とは、プラント１１内部の種々の機器の特性をいう。本実施形態において、機器特性とは、例えば、プラント１１の触媒層１０１の劣化度合い、加熱路１０２と触媒層１０１との間の総括伝熱係数等が挙げられる。本実施形態では、カルマンフィルタによって、プラント１１の状態量と係数パラメータとを同時に推定する。係数パラメータは、多次元とすることができる。

係数パラメータθの係数パラメータ変動モデルは次の式（３）で表現される。

式（３）における記号の意味は、次のとおりである。ｇは、係数パラメータの動的ダイナミクスを表す線形又は非線形の関数である。ｖ_θ（ｎ）は、離散時刻ｎにおける、θ_（ｎ）の更新に対するプロセスノイズである。Ｂ_θは、プロセスノイズｖ_θ（ｎ）の係数行列である。

なお、係数パラメータの変動に関して、物理的知見又は化学的知見、経験的な知見がない場合、次の式（４）で表されるランダムウォークモデルを用いてもよい。

ノイズ情報とは、カルマンフィルタを設計するために事前情報として必要な、プロセスノイズ、及び観測ノイズに関する情報である。プロセスノイズ及び観測ノイズは、次の式（５）によって仮定される。

ここで、Ｅは期待値を算出することを意味し、ノイズ情報Ｑ_ｘ、Ｑ_θ、Ｒはすべて半正定行列である。カルマンフィルタを設計する際、ノイズ情報Ｑ_ｘ、Ｑ_θ、Ｒは設計者が定める。ノイズ情報Ｑ_ｘ、Ｑ_θ、Ｒは、それぞれ次のような意味を持つ。Ｑ_ｘは、状態量のプロセスノイズの分散・共分散を表し、フィルタ設計の段階では、チューニングパラメータとすることができる。Ｑ_θは、推定する係数パラメータのプロセスノイズの分散・共分散を表し、係数パラメータの推定の収束速度に関係するチューニングパラメータである。Ｒは観測ノイズの分散・共分散を表し、一般に観測部１２のセンサの仕様によって定まる。

なお、式（１）に示す動特性モデルと式（３）で示す係数パラメータ変動モデルとを、次の式（６）で示すようにまとめて記述することができる。

また、式（１）で示す動特性モデルの状態量と係数パラメータとを、まとめて拡張モデルの状態量Ｘ_（ｎ）＝［ｘ_（ｎ），θ_（ｎ）］^Ｔとすると、式（２）及び式（６）から、拡張モデルを次の式（７）で表すことができる。この拡張モデルとは、係数パラメータを状態量に含めたモデルである。

推定部２２Ａは、式（７）で表される方程式に対し、ノイズ情報Ｑ_ｘ、Ｑ_θ、Ｒを用いてカルマンフィルタを設計する。推定部２２Ａは、設計したアルゴリズムに対して、通信部２１で取得されたプラント１１の観測値及び操作量を入力してカルマンフィルタ処理を行う。これにより、推定部２２Ａは、現在のプラント１１における拡張モデルの状態量、及び誤差共分散行列（推定誤差の大きさ）を推定する。カルマンフィルタを用いた状態量等の推定は、既知の方法を用いて行うことができる。

なお、式（７）を用いて推定される拡張モデルの状態量は、上述のように、Ｘ_（ｎ）＝［ｘ_（ｎ），θ_（ｎ）］^Ｔとなっている。すなわち、推定部２２Ａが拡張モデルの状態量を推定することとは、プラント１１の動特性モデルにおける状態量、及び係数パラメータを推定することを意味する。このように、カルマンフィルタを用いて拡張モデルの状態量を推定することで、プラント１１の動特性モデルにおける状態量、及び係数パラメータを同時に推定できる。

予測部２２Ｂは、推定部２２Ａにおいて推定されたプラント１１の状態量及び係数パラメータと、プラント１１の動特性モデル及び係数パラメータ変動モデルとに基づいて、プラント１１の将来の状態量を予測する。また、予測部２２Ｂは、将来の状態量に加え、推定部２２Ａにおいて推定された誤差共分散行列に基づいて、予測した将来の状態量の予測幅を算出する。なお、ここでの「将来」とは、推定部２２Ａにおいて状態量の推定の対象とした時点に対して、時間的に後の時点である。

ここで、推定部２２Ａは、現在の時点におけるプラント１１の状態量等を推定することに加え、過去のある時点におけるプラント１１の状態量等を推定することもできる。すなわち、推定部２２Ａが現在の時点におけるプラント１１の状態量等を推定している場合には、現在の時点に対して時間的に後の時点におけるプラント１１の状態量を推定する。また、推定部２２Ａが過去のある時点におけるプラント１１の状態量等を推定している場合には、過去のある時点に対して時間的に後の時点におけるプラント１１の状態量を推定する。すなわち、過去のある時点に対して時間的に後の時点とは、過去のある時点以降の時点であり且つ現在から見て過去の時点、現在、及び現在から見て将来の時点を含む。

まず、予測部２２Ｂが将来の状態量を予測する構成について具体的に説明する。予測部２２Ｂは、次の式（８）を用いて将来の状態量を算出する。但し、予測時における現在の離散時刻をｎとする。また、ｋ＝０，１，２，３，・・・とする。推定部２２Ａは、推定部２２Ａで推定された拡張モデルの状態量を初期値として、予測したい所望の時刻（予測時間）となるまで、式（８）を繰り返し計算する。

なお、式（８）において、操作量ｕ_{（ｎ＋ｋ）}：（ｋ＝０，１，２，・・・）は任意に設定することができる。例えば、現在の操作量を維持したときに状態量がどのように推移するかを予測したい場合には、ｕ_{（ｎ＋ｋ）}＝ｕ_（ｎ）とする。また、操作量の切り替えに対する応答を予測したい場合には、その切り替え後の操作量を用いる。

また、式（８）は、式（７）で示すカルマンフィルタの拡張モデルから、ノイズ成分を取り除いた方程式に相当する。すなわち、式（８）は、プラント１１の動特性モデル及び係数パラメータ変動モデルに基づいて生成された数式モデルである。予測部２２Ｂは、推定部２２Ａで推定された拡張モデルの状態量を初期値として用いている。上述したように、拡張モデルの状態量には、プラント１１の動特性モデル（式（１））の状態量、及び係数パラメータが含まれている。従って、予測部２２Ｂが拡張モデルの状態量及び式（８）に基づいてプラント１１の将来の状態量を予測することとは、推定部２２Ａにおいて推定された現在のプラント１１の状態量及び係数パラメータと、プラント１１の動特性モデル及び係数パラメータ変動モデルとに基づいて、プラント１１の将来の状態量を予測することを意味する。

式（８）は、式（７）で示すカルマンフィルタの拡張モデルから、ノイズ成分を取り除いた方程式としたが、ノイズ成分を取り除かなくてもよい。ノイズ成分を残す場合、例えば、乱数を用いてノイズ成分を設定してもよい。また、式（８）を用い、ノイズ成分を変えて将来の状態量を複数予測してもよい。この場合、ノイズ成分の変化によってどの程度、状態量が変動するかを把握することができる。

図３に、予測部２２Ｂによって予測された状態量の一例を示す。図３では、現在の状態量から将来の状態量を推定した場合を示している。現在よりも過去の時点の状態量（実線部分）は、実測値（観測値）である。また、現在よりも将来の時点の状態量（破線部分）は、予測部２２Ｂによって予測された予測値である。

次に、予測部２２Ｂが将来の状態量の予測幅を算出する構成について具体的に説明する。予測部２２Ｂは、将来の状態量の予測幅を算出する際に、推定部２２Ａにおいて推定された誤差共分散行列を用いる。この誤差共分散行列は、推定部２２Ａで推定された状態量の不確かさを表している。このため、予測部２２Ｂは、推定部２２Ａで推定された誤差共分散行列が表す不確かさに基づいて、予測される将来の状態量の予測幅を算出する。この予測幅とは、推定された状態量の確からしさが、所定値以上である範囲を表している。図３に、予測部２２Ｂで算出された予測幅Ｗの一例を点線で示す。

予測部２２Ｂは、誤差共分散行列が表す状態量の不確かさに基づいて、種々の方法によって状態量の予測幅を算出することができる。以下、予測幅の具体的な算出方法について説明する。なお、ここでは、現在の状態量から将来の状態量が予測されたときに、将来の状態量の予測幅を算出する場合について説明する。

まず、予測幅の第１の算出方法について説明する。予測部２２Ｂは、推定部２２Ａで推定された誤差共分散行列に基づいて、推定された状態量の不確かさを取得する。この取得された不確かさとは、推定された状態量の現在の時点における不確かさである。予測部２２Ｂは、現在の時点における不確かさが、状来の時点においても続くと仮定して予測幅を算出する。この場合、図４に示すように、現在の時点以降、予測幅Ｗの範囲は一定となり、且つ、予測された状態量の値に対する予測幅Ｗの位置が一定となる。

次に、予測幅の第２の算出方法について説明する。予測部２２Ｂは、推定部２２Ａで推定された誤差共分散行列に基づいて、推定された状態量の不確かさを取得する。図５に示すように、予測部２２Ｂは、取得した不確かさに基づいて、推定された現在の時点における状態量の確からしさが所定値以上である範囲Ｗ１を設定する。予測部２２Ｂは、範囲Ｗ１の中で、状態量の初期値を変えて複数回、式（８）を用いて状態量の予測を行う。予測部２２Ｂは、複数回予測した状態量の広がりに基づいて、予測幅を算出する。図５に示す例では、範囲Ｗ１の上限値を状態量の初期値として予測を行った場合を予測値Ｓ１とし、範囲Ｗ１の下限値を状態量の初期値として予測を行った場合を予測値Ｓ２としている。この場合、例えば、予測部２２Ｂは、予測値Ｓ１と予測値Ｓ２との間を、予測幅Ｗとして算出することができる。

なお、推定部２２Ａで推定された状態量には、複数の項目（温度、流量、密度など）が含まれている。予測部２２Ｂは、推定部２２Ａで推定された状態量の中で、不確かさが所定値以上の項目について、各項目の不確かさに基づいて設定された範囲の中で各項目の初期値を変えて予測幅を算出してもよい。このように、予測部２２Ｂは、推定部２２Ａで推定された状態量に含まれる全ての項目について初期値を変えて予測幅を行うのではなく、予測幅の変化に与える影響が大きい項目についてのみ初期値を変えて予測幅の算出を行う。これにより、予測幅の算出の計算量を効果的に削減することができる。

次に、予測幅の第３の算出方法について説明する。予測部２２Ｂは、現在から過去の所定の時点までの間に推定部２２Ａで推定された誤差共分散行列を用い、これらの誤差共分散行列が表す状態量の不確かさの変動に基づいて、将来の状態量の予測幅を算出する。例えば、予測部２２Ｂは、現在までの不確かさの変動に基づいて数式を設定し、設定した数式を用いて将来の状態量の予測幅を算出してもよい。

なお、予測幅の算出方法として、第１の算出方法から第３の算出方法を例に説明したが、これらの方法以外の方法を用いて、予測幅を算出してもよい。

次に、推定部２２Ａが行う拡張モデルの状態量等の推定処理の流れについて説明する。ここでは、プラント１１と同期しながら、カルマンフィルタにより状態量のオンライン推定を行う。式（５）で示すノイズ情報は、予め設定されているとする。現在の時刻を（ｎ）で表す。

図６に示すように、推定部２２Ａは、拡張モデルの現在の状態量の推定値Ｘ_ｈ（ｎ）及び誤差共分散行列Ｐ_（ｎ）と、式（７）の拡張モデルと用いて、拡張モデルの次の時刻における状態量の事前推定値Ｘ_ｈ ⁻ _{（ｎ＋１）}及び事前誤差共分散行列Ｐ⁻ _{（ｎ＋１）}を算出する（Ｓ１０１）。なお、Ｓ１０１の処理を初めて行う場合、拡張モデルの現在の状態量の推定値Ｘ_ｈ（ｎ）及び誤差共分散行列Ｐ_（ｎ）として、予め設定された拡張モデルの状態量の初期推定値Ｘ_ｈ（０）及び誤差共分散行列の初期値Ｐ_（０）を用いる。また、Ｓ１０１の処理が２回目以降である場合、拡張モデルの現在の状態量の推定値Ｘ_ｈ（ｎ）及び誤差共分散行列Ｐ_（ｎ）として、後述のＳ１０４で算出された拡張モデルの状態量の事後推定値Ｘ_{ｈ（ｎ＋１）}及び事後誤差共分散行列Ｐ_{（ｎ＋１）}を用いる。

推定部２２Ａは、カルマンゲインＧを算出する（Ｓ１０２）。事前にカルマンゲインＧが決定（設計）している場合、このＳ１０２の処理を省略してもよい。推定部２２Ａは、通信部２１を介してプラント１１の観測値ｙ_（ｎ）及び操作量ｕ_（ｎ）を取得する（Ｓ１０３）。

推定部２２Ａは、観測値ｙ_（ｎ）、操作量ｕ_（ｎ）、及びカルマンゲインＧを用いて、状態量の事前推定値Ｘ_ｈ ⁻ _{（ｎ＋１）}及び事前誤差共分散行列Ｐ⁻ _{（ｎ＋１）}から、状態量の事後推定値Ｘ_{ｈ（ｎ＋１）}及び事後誤差共分散行列Ｐ_{（ｎ＋１）}を算出する（Ｓ１０４）。状態量の事後推定値Ｘ_{ｈ（ｎ＋１）}及び事後誤差共分散行列Ｐ_{（ｎ＋１）}の算出後、推定部２２Ａは、ｎ＝ｎ＋１として、予め定められたタイミングで上述したＳ１０１から再び処理を行う。

なお、Ｓ１０１、Ｓ１０２、及びＳ１０４の具体的な計算方法は、用いるカルマンフィルタの種類によって異なる。一般に、拡張モデル（係数パラメータを状態量に含めたモデル）は非線形となるので、非線形カルマンフィルタを用いることが好ましい。特に、非線形カルマンフィルタとしては、ＵＫＦ（Unscented Kalman Filter：アンセンテッドカルマンフィルタ）を用いることが有効である。これは、従来のようなシステムの線形近似（勾配情報）を用いた、ＥＫＦ（Extended Kalman Filter：拡張カルマンフィルタ）と異なり、勾配情報を必要としないため、より実装が容易となる。また、一般に、ＵＫＦはＥＫＦよりも、安定性や推定性能が高い。

次に、予測部２２Ｂが行う将来の状態量の予測及び予測幅の算出処理の流れについて説明する。なお、図７に示す処理は、将来の状態量の予測を実施したい任意のタイミングで開始される。図７に示すように、予測部２２Ｂは、推定部２２Ａで推定された拡張モデルの現在の状態量及び誤差共分散行列を取得する（Ｓ２０１）。推定部２２Ａは、操作量ｕ及び予測時間を設定する（Ｓ２０２）。この操作量ｕ及び予測時間は、状態予測装置２のオペレータ等が設定してもよく、予め設定された値であってもよい。

予測部２２Ｂは、取得した拡張モデルの現在の状態量及び誤差共分散行列と、設定された操作量ｕ及び予測時間と、式（８）とに基づいて、拡張モデルの将来の状態量を予測する（Ｓ２０４）。また、予測部２２Ｂは、取得した誤差共分散行列に基づいて、状態量の予測幅を算出する（Ｓ２０４）。予測部２２Ｂは、予測した将来の状態量及び算出した予測幅を出力する（Ｓ２０５）。

予測部２２Ｂから出力された将来の状態量及び予測幅は、状態予測装置２の通信部２１からネットワークＮＷを介してプラントシステム１に送信される。表示部１４は、通信部１５で受信された将来の状態量及び予測幅を表示する。プラント１１のオペレータ等は、表示部１４に表示された将来の状態量及び予測幅を見て、プラント１１の制御を行うことができる。また、状態予測装置２の表示部２３は、予測部２２Ｂから出力された将来の状態量及び予測幅を表示する。

本実施形態は以上のように構成され、状態予測装置２の推定部２２Ａは、プラント１１の動特性モデル等を用いてカルマンフィルタを設計し、カルマンフィルタ処理を行うことで、プラント１１の状態量及び係数パラメータを推定する。このように、推定部２２Ａは、カルマンフィルタを用いて推定を行うことにより、プラント１１の状態量及び係数パラメータを精度良く推定することができる。予測部２２Ｂは、カルマンフィルタを用いて推定されたプラント１１の状態量及び係数パラメータ等を用いて、プラント１１の将来の状態量を予測する。このように、予測部２２Ｂは、カルマンフィルタを用いて精度良く推定されたプラント１１の状態量及び係数パラメータを用いて予測を行うことで、プラント１１の将来の状態量をより精度良く予測することができる。

また、予測部２２Ｂは、将来の状態量に加えて状態量の予測幅を算出する。これにより、プラント１１のオペレータ等は、予測されたプラント１１の将来の状態量に加えて、予測された状態量の予測幅を把握することができる。例えば、プラント１１の状態量が予め定められた範囲を超えないようにプラント１１の制御を行う際に、予測幅が大きい場合と小さい場合とでは、例えばオペレータによるプラント１１の制御の実施タイミング、及び制御アルゴリズムにおける制御量（操作量、バルブ開閉の指令値など）の大きさ等が異なることがある。例えば、予測幅が小さい場合には、プラント１１に対して大きな制御量を与えてプラント１１を早く動かすことができる。また、予測幅が大きい場合には、プラント１１に対して小さな制御量を与えて慎重にゆっくり動かすことができる。このように、予測されたプラント１１の将来の状態量に加えて、状態量の予測幅も算出することで、プラント１１のオペレータ等は、予測幅に応じてプラント１１の制御を行うことができる。

プラントシステム１と状態予測装置２とは、ネットワークＮＷを介して互いに接続されている。このため、状態予測装置２は、プラントシステム１の設置場所とは異なる場所において、プラント１１の状態量等を予測することができる。また、状態量等の予測処理を行う場合、予測処理における計算量が多いため、処理能力の高いコンピュータを用いる必要がある。このような場合であっても、プラントシステム１側において処理能力の高いコンピュータを設置する必要が無い。また、１台の状態予測装置２に対して、複数のプラントシステム１が接続されていてもよい。この場合、１台の状態予測装置２において複数のプラントシステム１の将来の状態量等をまとめて予測することができる。

以上、本発明の実施形態について説明したが、本発明は上記実施形態に限定されない。例えば、上記実施形態では、プラント１１の一例として、熱交換型の流通反応器が用いられているがこれに限られない。例えば、プラント１１として、ガスタービンエンジン、水素製造装置、又は焼却炉等、種々のプラントを用いてもよい。

上記では、有効なカルマンフィルタとしてＵＫＦを挙げたが、カルマンフィルタの拡張概念であるベイズフィルタ（ベイズ理論に基づくフィルタ）に含まれるフィルタであれば、種々のフィルタの利用が可能である。ベイズフィルタとして、例えば、上述したＥＫＦ、アンサンブル・カルマンフィルタ、粒子フィルタ、情報フィルタといったフィルタの利用が考えられる。また、Ｈ∞フィルタは、ベイズフィルタに分類されない場合もあるが、このＨ∞フィルタを用いてもよい。

予測部２２Ｂは、必ずしも式（８）に示すような動的モデルを用いて将来の状態量等を予測するシミュレータである必要はない。例えば、予測部２２Ｂは、過去の運転データから作成された静的な回帰モデルを用いて将来の状態量等を予測してもよい。

プラントシステム１の制御部は、プラント１１に対してモデル予測制御を行ってもよい。このモデル予測制御とは、将来の状態量をモデルによって予測しながら、適切な将来の操作量（式（８）におけるｕ_{（ｎ＋ｋ）}）を決定する技術のことである。状態予測装置２の予測部２２Ｂは、予測幅まで含めて将来の状態量を予測できる。このため、例えば、予測幅を含めた将来の状態量の予測値の領域が、設定されたある領域に入らないよう、適切に将来の操作量を決定することができる（ロバストモデル予測制御）。このようなモデル予測制御を行うことで、最適性を考慮しつつ、より安全なプラントの運転を実現できる。

プラントシステム１と状態予測装置２とは互いに同じ場所に設置されていてもよい。この場合、状態予測装置２は、プラントシステム１の機能の一部として組み込まれていてもよい。

２状態予測装置
１１プラント（対象装置）
２２Ａ推定部
２２Ｂ予測部
Ｗ予測幅

Claims

状態量の予測対象とする対象装置の動特性モデル、及び前記対象装置の機器特性を表す係数パラメータに関する係数パラメータ変動モデルを用いてカルマンフィルタを設計し、設計されたアルゴリズムに対して前記対象装置の観測値を入力してカルマンフィルタ処理を行うことで、前記対象装置の現在の状態量又は過去のある時点の状態量、及び前記係数パラメータを推定する推定部と、
推定された前記状態量、推定された前記係数パラメータ、前記動特性モデル、及び前記係数パラメータ変動モデルに基づいて、前記推定部において前記状態量の推定の対象とした時点に対して時間的に後の時点における前記対象装置の状態量を予測する予測部と、を備え、
前記推定部は、前記状態量及び前記係数パラメータを推定するとともに、推定された前記状態量と実際の状態量との推定誤差の大きさを推定し、
前記予測部は、前記推定誤差の範囲内において前記状態量の初期値を変えて前記状態量の予測を複数回行い、複数回予測した前記状態量の拡がりに基づいて、予測された前記状態量の予測幅を算出する、状態予測装置。
状態量の予測対象とする対象装置の動特性モデル、及び前記対象装置の機器特性を表す係数パラメータに関する係数パラメータ変動モデルを用いてカルマンフィルタを設計し、設計されたアルゴリズムに対して前記対象装置の観測値を入力してカルマンフィルタ処理を行うことで、前記対象装置の現在の状態量又は過去のある時点の状態量、及び前記係数パラメータを推定する推定部と、
推定された前記状態量、推定された前記係数パラメータ、前記動特性モデル、及び前記係数パラメータ変動モデルに基づいて、前記推定部において前記状態量の推定の対象とした時点に対して時間的に後の時点における前記対象装置の状態量を予測する予測部と、を備え、
前記推定部は、前記状態量及び前記係数パラメータを推定するとともに、推定された前記状態量と実際の状態量との推定誤差の大きさを推定し、
前記予測部は、現在又は過去のある時点から過去の所定時間分の前記推定誤差を用い、前記推定誤差が表す前記状態量の不確かさの変動に基づいて、予測された前記状態量の予測幅を算出する、状態予測装置。