JP5977207B2 - 状態推定装置、該方法及び該プログラム - Google Patents

Description

本発明は、状態方程式を用いた状態推定方法に関する。

現実のシステムを制御する際には、様々な要因によってその制御量が目標値からばらついてしまうことがある。そのため、そのようなバラツキを抑えるための高精度な制御の実現が必要とされている。

例えば、大規模生産を行う鉄鋼プラントでは、連続鋳造機や転炉、圧延機、鋼板冷却設備、焼鈍炉など、その殆どがバッチプロセスである。連続化されている連続圧延機や連続焼鈍炉でも、板毎に見るとバッチプロセスと見ることができる。このようなバッチ毎の確率変動や製品毎の確率変動の存在下において、高品質な鋼板を安定して生産するが強く要求されている。しかし、種々の要因により製品にばらつきが生じてしまうことがある。これは、現実のシステムでは内部状態を直接観測できない場合があるためである。このため、システムの内部状態を推定することができれば、製品のばらつきを抑制できることが可能となる。

確率システムに対する代表的な状態推定手法として、線形システムにおいて誤差の二乗平均を最小にするように推定量を求めるカルマンフィルタ（非特許文献１）がある。また、非線形確率システムに関する状態推定法として、状態の非線形関数の推定値近傍での線形近似に基づく拡張カルマンフィルタや、期待値計算にモンテカルロ法を用いたパーティクルフィルタが知られている。

更に、確率変数をパラメータに持つ系に対する状態推定問題はＤｅ Ｋｏｎｉｎｇにより考えられている（非特許文献２）。この手法は確率変数パラメータの影響を考慮したうえでカルマンフィルタの問題設定に帰着させるものとなっており、パラメータが確定値であることを前提としている通常のカルマンフィルタでは解くことができない問題に有効である。

Ｒ．Ｅ．Ｋａｌｍａｎ．Ａ ｎｅｗ ａｐｐｒｏａｃｈ ｔｏ ｌｉｎｅａｒ ｆｉｌｔｅｒｉｎｇ ａｎｄ ｐｒｅｄｉｃｔｉｏｎ ｐｒｏｂｌｅｍｓ． Ｊｏｕｒｎａｌ ｏｆ Ｂａｓｉｃ Ｅｎｇｉｎｅｅｒｉｎｇ，Ｖｏｌ．８２，ｐｐ．３５−４５，１９６０． Ｗ．Ｌ．Ｄｅ Ｋｏｎｉｎｇ． Ｏｐｔｉｍａｌ ｅｓｔｉｍａｔｉｏｎ ｏｆ ｌｉｎｅａｒ ｄｉｓｃｒｅｔｅ−ｔｉｍｅ ｓｙｓｔｅｍｓ ｗｉｔｈ ｓｔｏｃｈａｓｔｉｃ ｐａｒａｍｅｔｅｒｓ． Ａｕ−ｔｏｍａｔｉｃａ，Ｖｏｌ．２０，Ｎｏ．１，ｐｐ．１１３−１１５，１９８４．

しかし、このＤｅ Ｋｏｎｉｎｇによる手法では、システムパラメータは毎時刻確率的に変化する値とし扱われ、現実のシステムとは異なる場合がある。つまり、現実のシステムにおいてパラメータは緩やかな時間変化もしくは製造製品ごとに変わる場合が多いと考えられる。

そこで、本発明は、システムパラメータを時間に依存しない確率変数として扱う状態推定法を用いて状態を推定する状態推定装置、状態推定方法及び状態推定プログラムを提供することを目的とする。

上述の目的を達成するために、本発明に係る状態推定装置は、バッチプロセスから観測情報を取得する取得部と、前記観測情報から、状態方程式を用いて、前記バッチプロセスにおける所望の時刻の状態値を推定する状態演算部とを備え、前記状態方程式の状態行列Ａ、入力行列Ｂ及び出力行列Ｃのうちの少なくとも１つ以上が、時間変数に無関係な或る確率分布に従ってばらつくものであり、状態行列Ａ^ｎ（ｎ≧１：ｎは整数）の期待値を、過去の時刻の状態値の係数行列として用い、あるいは、Ａ^ｎを含む項の期待値として用いることを特徴とする。

そして、本発明に係る状態推定方法は、観測情報から、所望の時刻の状態値を推定する状態推定装置で用いられる状態推定方法であって、バッチプロセスから前記観測情報を取得する取得ステップと、前記観測情報から、状態方程式を用いて、前記バッチプロセスにおける所望の時刻の状態値を推定する状態演算ステップとを備え、前記状態方程式の状態行列Ａ、入力行列Ｂ及び出力行列Ｃのうちの少なくとも１つ以上が、時間変数に無関係な或る確率分布に従ってばらつくものであり、状態行列Ａ^ｎ（ｎ≧１：ｎは整数）の期待値を、過去の時刻の状態値の係数行列として用い、あるいは、Ａ^ｎを含む項の期待値として用いることを特徴とする。

また、本発明に係る、コンピュータに実行させるための状態推定プログラムは、観測情報から、所望の時刻の状態値を推定する状態推定装置で用いられる状態推定プログラムであって、バッチプロセスから前記観測情報を取得する取得手段と、前記観測情報から、状態方程式を用いて、前記バッチプロセスにおける所望の時刻の状態値を推定する状態演算手段としてコンピュータを機能させ、前記状態方程式の状態行列Ａ、入力行列Ｂ及び出力行列Ｃのうちの少なくとも１つ以上が、時間変数に無関係な或る確率分布に従ってばらつくものであり、状態行列Ａ^ｎ（ｎ≧１：ｎは整数）の期待値を、過去の時刻の状態値の係数行列として用い、あるいは、Ａ^ｎを含む項の期待値として用いることを特徴とする。

また、上述の状態推定装置において、前記状態演算部は、ｗハットをシステムノイズの推定値、ｕを観測情報、Ｅ[ ]を期待値演算子、Ｎを過去の状態値の数とするとき、以下の状態方程式により、時刻ｋの時の状態値ｘハットを算出することが好ましい。

このような構成の状態推定装置、状態推定方法及び状態推定プログラムによれば、状態方程式の状態行列Ａ、入力行列Ｂ及び出力行列Ｃのうちの少なくとも１つ以上を、或る確率分布に従ってばらつくものとして扱い、その期待値を用いて状態値を推定するので、より確からしい状態値を推定することが可能となる。尚、Ａ^ｎを含む項とは、ｎステップ過去の状態値が、現在の状態値に与える影響を表す係数行列である。

また、上述の状態推定装置において、前記確率分布は、状態値を推定する製品に応じて設定されることが好ましい。

この構成によれば、製品ごとに適した状態推定方程式を用いて状態値を推定することが可能となる。

このような構成の状態推定装置、状態推定方法及び状態推定プログラムでは、システムパラメータを時間に依存しない確率変数として扱う状態推定法を用いて状態を推定することができる。

状態推定装置の構成を示すブロック図である。 図１の状態推定装置の状態推定処理のフローチャートである。 板の内部温度の変化を示す図である。 遷移領域での熱伝導率の確率分布を示す図である。 図５（ａ）は、確率変数ａが従う一様分布の確率密度関数を示す図であり、図５（ｂ）は、確率変数ｂ、ｃが従う一様分布の確率密度関数を示す図である。 １００サンプルの確率変数ａの値を示す図である。 １００サンプルの確率変数ａに対して行列Ａの固有値の大きさを示す図である。 カルマンフィルタを用いた場合に誤差ノルムの値を示す図である。 実施形態の推定方法を用いた場合の誤差ノルムの値を示す図であり、図９（ａ）は、設計パラメータＮ＝０のときの誤差ノルムを示し、図９（ｂ）は、設計パラメータＮ＝１のときの誤差ノルムを示し、図９（ｃ）は、設計パラメータＮ＝２のときの誤差ノルムを示し、図９（ｄ）は、設計パラメータＮ＝３のときの誤差ノルムを示す図である。 サンプルの平均値及び標準偏差の表を示す図である。 実施形態の方法を用いた場合の誤差ノルムの平均値の遷移を示す図である。 取り鍋の保温能力の確率分布を示す図である。 圧延材の塑性係数の確率分布を示す図である。

＜実施形態１＞
以下、本発明に係る実施形態を図面に基づいて説明する。
＜構成＞
図１は、状態推定装置の構成を示すブロック図である。図１において、状態推定装置１は、演算処理部１１、入力部１２、出力部１３、内部記憶部１４、補助記憶部１６及びバス１８を備えて構成される。

演算処理部１１は、例えば、マイクロプロセッサ等を備えて構成され、機能的に、制御部１１１、初期設定部１１２、前時刻推定値算出部１１３、及び、最適推定値算出部１１４を備える。

制御部１１１は、初期設定部１１２、前時刻推定値算出部１１３、及び、最適推定値算出部１１４を制御して状態推定の演算を行い、また、制御プログラムに従い入力部１２、出力部１３、内部記憶部１４及び補助記憶部１６を制御する。

初期設定部１１２は、状態推定処理を行う際の初期値を設定する機能を有する。

前時刻推定値算出部１１３は、推定しようとする時刻よりも前の時刻の状態の推定値及びシステムノイズの推定値を算出する機能を有する。

最適推定値算出部１１４は、前時刻推定値算出部１１３が算出した推定値を用いて最適推定値を算出する。これらの詳細については、＜動作＞の項で説明する。

入力部１２は、本状態推定装置１の演算開始指示等の各種コマンドやシステムパラメータ等の各種データを状態推定装置１に入力する機器であり、例えば、キーボードやマウス等である。出力部１３は、入力部１２から入力されたコマンド及びデータや本状態推定装置１の演算結果等を出力する機器であり、例えばＣＲＴ（Ｃａｔｈｏｄｅ Ｒａｙ Ｔｕｂｅ）ディスプレイ、ＬＣＤ（Ｌｉｑｕｉｄ Ｃｒｙｓｔａｌ Ｄｉｓｐｌａｙ）、有機ＥＬ（Ｅｌｅｃｔｒｏ Ｌｕｍｉｎｅｓｅｎｃｅ）ディスプレイ又はプラズマディスプレイ等の表示装置やプリンタ等の印字装置等である。

内部記憶部１４は、演算処理部１１が実行する状態推定プログラムや制御プログラムを補助記憶部１６から読み込むと共に、状態推定プログラム実行中の各データを一時的に記憶する所謂ワーキングメモリであり、例えば揮発性の記憶素子であるＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）を備えて構成される。

補助記憶部１６は、例えばＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）及びＥＥＰＲＯＭ（Ｅｌｅｃｔｒｉｃａｌｌｙ Ｅｒａｓａｂｌｅ Ｐｒｏｇｒａｍｍａｂｌｅ Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）等の不揮発性の記憶素子やハードディスク等のデータやプログラムを記憶する装置であり、測定値記憶部１６１、状態空間モデル記憶部１６２、及び、推定値記憶部１６３を備える他、状態推定プログラムや状態推定装置１を動作させるための制御プログラム等の各プログラム（不図示）、及び、各プログラムの実行に必要なデータ（不図示）等を記憶する。

測定値記憶部１６１は、観測データ等の推定値の演算に必要なデータ、例えば、厚板の冷却工程における厚板の表面温度等を記憶する。

状態空間モデル記憶部１６２は、推定対象の状態空間モデル、例えば、厚板の冷却工程における内部温度を推定する厚板冷却モデルを記憶する機能を有する。

推定値記憶部１６３は、最適推定値算出部１１４が算出した時刻毎の推定値や、前時刻推定値算出部１１３が算出した前時刻の状態の推定値、システムノイズの推定値等を必要に応じて記憶しておく機能を有する。

これら演算処理部１１、入力部１２、出力部１３、内部記憶部１４及び補助記憶部１６は、データを相互に交換することができるようにバス１８にそれぞれ接続される。

なお、必要に応じて分子軌道演算装置１は、図１に破線で示す外部記憶部１５や通信インターフェース部（不図示）をさらに備えてもよい。

外部記憶部１５は、例えば、フレキシブルディスク、ＣＤ−ＲＯＭ（Ｃｏｍｐａｃｔ Ｄｉｓｃ Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、ＣＤ−Ｒ（Ｃｏｍｐａｃｔ Ｄｉｓｃ Ｒｅｃｏｒｄａｂｌｅ）及びＤＶＤ−Ｒ（Ｄｉｇｉｔａｌ Ｖｅｒｓａｔｉｌｅ Ｄｉｓｃ Ｒｅｃｏｒｄａｂｌｅ）等の記録媒体との間でデータを読み込み及び／又は書き込みを行う装置であり、例えば、フレキシブルディスクドライブ、ＣＤ−ＲＯＭドライブ、ＣＤ−Ｒドライブ及びＤＶＤ−Ｒドライブ等である。通信インターフェース部は、通信網に接続され、通信網を介して他のクライアント装置との間で通信信号を送受信するための機器である。

各プログラムが格納されていない場合には、これらを記録した記録媒体から外部記憶部１５を介して補助記憶部１６にインストールされるように構成してもよく、また、これらプログラムを管理するサーバコンピュータ（不図示）から通信網及び通信インターフェース部を介して各プログラムがダウンロードされるように構成してもよい。また、状態推定の演算に当たって状態推定装置１に入力すべきデータは、このデータを記憶した記録媒体によって外部記憶部１５を介して状態推定装置１に入力されるように構成してもよく、また、クライアント装置から通信網及び通信インターフェース部を介して状態推定装置１に入力されるように構成してもよい。

次に、本実施形態の状態推定方法について説明する。尚、以下、式中で英文字記号（変数）に＾（ハット）を付すものは、「（英文字記号）ハット」と記載し、英文文字に￣（オーバーライン）を付すものは、「（英文字記号）オーバーライン」と記載するものとする。
＜状態推定方法＞
現実のシステムにおいては、システムパラメータ（Ａ，Ｂ，Ｃの要素）が或るばらつきをもって変化（確率変動）することを考慮し、実施形態の状態推定方法は、その影響が最も顕著に現れる期待値Ｅ[Ａ^ｎ]を用いる。

まず、既存の状態推定方法であるカルマンフィルタを考える。通常のカルマンフィルタは、以下の離散時間線形システムを対象とする。
ｘ_ｋ＋１＝Ａｘ_ｋ＋Ｂｕ_ｋ＋ｗ_ｋ ・・・（Ｋ１）
ｙ_ｋ＝Ｃｘ_ｋ＋ｖ_ｋ ・・・（Ｋ２）
ｐ（ｗ_ｋ）＝Ｎ（ｗ_ｋ｜０，Ｑ） ・・・（Ｋ３）
ｐ（ｖ_ｋ）＝Ｎ（ｖ_ｋ｜０，Ｒ） ・・・（Ｋ４）
ここで、ｘ∈Ｒ^ｎが状態、ｕ∈Ｒ^ｍが制御入力、ｙ∈Ｒ^ｐが出力、ｗ∈Ｒ^ｎがシステムノイズ、ｖ∈Ｒ^ｐが観測ノイズである。Ａ、Ｂ、Ｃは、それぞれｎ×ｎ、ｎ×ｍ、ｐ×ｎの大きさをもつ行列であり、Ａを状態行列、Ｂを入力行列、Ｃを出力行列という。Ｒ^ｎは、ｎ次元ユークリッド空間を表す。

ｐ（ａ）は、ａの確率密度関数（確率密度）を表し、ｐ（ａ）＝Ｎ（ａ｜μ，Σ）は、確率変数ａがガウス分布に従うことを示す。μは平均を示し、Σは共分散行列を示す。従って、式（Ｋ３）は、ｗ_ｋの確率分布ｐ（ｗ_ｋ）が、平均０、分散Ｑの正規分布に従うことを示し、式（Ｋ４）は、ｖ_ｋの確率分布ｐ（ｖ_ｋ）が、平均０、分散Ｒの正規分布に従うことを示す。

ｘハットは状態ｘの推定値であり、ｋは自然数であって時刻を示すものとする。時刻ｋでの状態の最適推定値をｘ_ｋハットとし、時刻ｋよりも１時点前での出力観測後の推定値の平均値をｘ_{ｋ｜ｋ−１}ハットとし、共分散をΠ_{ｋ｜ｋ−１}とすると、内部状態ベクトルｘ_ｋの推定値ｘ_ｋ｜ｋハットは、以下の式（Ｋ５）、（Ｋ６）、（Ｋ７）、（Ｋ８）を繰り返し計算することで求められる。

尚、（Ａ）^Ｔは、行列Ａ又はベクトルＡの転置を表す作用素とし、（Ａ）^−１は、行列Ａの逆行列を表す作用素とする。

このように、カルマンフィルタなどの従来の状態推定手法では、観測ノイズｖ_ｋやシステム外乱ｗ_ｋなどの加法的な外乱が確率的な変動をするものとして取り扱うことができる。しかし、実際のシステムでは、外乱のみが確率的な変動をするのではなく、パラメータそのものもあるばらつきを持っている。

最も代表的なカルマンフィルタではｖ_ｋやｗ_ｋが正規分布に従い、その分散が分かっていれば、もっともらしい推定値を与えることが示されているが、Ａ、Ｂ、Ｃ行列の要素（システムパラメータ）が確率的に変動する場合は扱うことができない（推定することができるが、その推定値が確からしいかどうかは、保証されない）。

そこで、実施形態では、以下のように、Ａ、Ｂ、Ｃ行列の要素（システムパラメータ）を、確率的な変動を有する確率変数として取り扱う。

実施形態の状態推定は、以下の離散時間線形システムを対象とする。

通常のカルマンフィルタが対象とするシステムと異なる点は、Ａ（ω）∈Ｒ^ｎ×ｎ、Ｂ（ω）∈Ｒ^ｎ×ｍ、Ｃ（ω）∈Ｒ^ｐ×ｎが確率変数である点である。尚、ωは偶発性を表す変数である。Ａ（ω）、Ｂ（ω）、Ｃ（ω）、とｗ_ｋ、ｖ_ｋは独立であると仮定する。

そして、以下の式（３）ように条件付き確率密度関数が最大となるときを最適推定値と考えるものとする。時刻ｋでの状態の最適推定値をｘ_ｋハットとし、考慮する過去の推定値の数をＮとする。ａｒｇ ｍａｘ ｆ（ａ，ｂ，ｃ，…）は、関数ｆ（ａ，ｂ，ｃ，…）が最大値をとるときのパラメータａ，ｂ，ｃ，…を示している。例えばａｒｇ ｍａｘ ｐ（ｘ｜ｙ）は、出力ｙの出力制約付き確率密度関数ｐ（ｘ｜ｙ）が最大値をとるときの状態ｘを示している。

式（１）、及び、式（２）より、システムのマルコフ性から状態ｘ_ｋの同時確率密度関数は、式（４）となり、出力ｙ_ｋの条件付き確率密度関数は、式（５）となる。

式（４）、及び、式（２）より、ベイズの定理を用いると式（３）式の確率密度関数は、式（６）とできる。

この関数は、対数を考えても最大値は変わらないため、最適推定値は、式（７）のように表せる。

ここで、状態ｘ_ｋ−Ｎは、平均ｘ_ｋ−１オーバーラインで分散Π_ｋ−Ｎの正規分布に従い、ノイズｗ_ｋ、ｖ_ｋは、平均０で分散がそれぞれＱ、Ｒの正規分布に従うと仮定すると、最適推定値ｘハットを求める際の評価関数は、次式（８）で表される。

従って、式（３）の最適推定値ｘハットは、式（９）式で求められることになる。ａｒｇ ｍｉｎ ｆ（ａ，ｂ，ｃ，…）は、関数ｆ（ａ，ｂ，ｃ，…）が最小値をとるときのパラメータａ，ｂ，ｃ，…を示している。例えばａｒｇ ｍｉｎ Ｅ［Ｊ（ｘ；ｕ）］は、入力ｕが入力時の評価関数Ｊ（ｘ；ｕ）の期待値が最小値をとるときの状態ｘを示している。

期待値Ｅ［ｘ］は、ｘの期待値演算を表し、具体的には、Ｅ[ｘ]＝∫ｘｐ（ｘ）ｄｘである。式（９）中の期待値演算Ｅ［ａ］は、ωに関する期待値を表しており、ここでは確率変数Ａ、Ｂ、Ｃの統計量を用いて計算を行う。

式（１）及び式（２）で表すシステムに対して、式（１０）で示す二次形式評価関数を考える。

そして、この評価関数を最小にする最適推定値ｘ_ｋハット（ｋ＝Ｎ、Ｎ＋１・・・）は、式（１１）で得ることができる。

ここで、Ｎ個前の推定値ｘ_ｋ−Ｎハットは、以下の式（１２）で求める。この式（１２）中の実測値の平均値ｘ_ｋ−Ｎオーバーラインは式（１３）で求め、分散Π_ｋ−Ｎを式（１４）で求める。

また、式（１１）中のシステムノイズの推定値ｗ_ｋ−ｑハットは、以下の式（１５）で求める。

但し、Ｙ_ｋ：ｑ、Ｏ_ｋ、Ｒ_ｋは、それぞれ以下の式（１６）、式（１７）、式（１８）で表される。尚、式（１８）は、対角成分のみを表し、他の成分はすべて０である。

ここで、Ｂ行列が確率変数ではない場合、すなわち、確定的である場合には、上述の式（１１）に代えて、以下の式（１９）としてもよい。すなわち、ｎ≧１のＡ^ｎの期待値Ｅ[Ａ^ｎ]を係数行列として算出し、各ベクトルｘ_ｋ−Ｎハット、Ｂ・ｕ_ｊ、ｗ_ｊハットに乗じることで現在の状態推定値ｘ_ｋハットを推定することができる。

＜状態推定モデルの例＞
ここでは、鋼鉄プラントにおいて熱間圧延により厚鋼板を製造する際の冷却工程（バッチプロセスの１例）で、板の表面温度から板内部の温度を推定する場合を例に説明する。

厚板の冷却工程においては、板内部の温度にて管理・制御することが重要であり、内部温度は表面温度とは大きく異なっている。従って、内部温度をより正確に推定する必要がある。

そこで、熱伝達率α_ｊなどを確率変数として扱い、板の厚み方向の熱伝導を考慮した状態方程式にて定式化し、Ａ^ｎの期待値Ｅ[Ａ^ｎ]を考慮したうえで状態推定することで、板の内部温度の変化を精度良く推定する。

図３に、板の内部温度の変化を示すグラフを示す。横軸は時間を表し、縦軸は内部温度を表す。このグラフは、板を厚み方向に１５に等分割し、表面温度と、等分割されたそれぞれ厚み位置における板の内部温度の変化を示している。板の内部の温度分布は、板の中央を中心に対称となっているため、板中央以外の温度変化のグラフは２本ずつ重なっている。従って、全部で９本のグラフとして見えている。また、図４は、遷移領域での熱伝達率の確立分布を示すグラフである。図４のグラフは、横軸が熱伝達率を表し、縦軸が確率を表す。遷移領域とは、膜沸騰から核沸騰に遷移する温度範囲をいい、この例では５００℃前後で遷移している。熱伝達率α_ｊは、膜沸騰と核沸騰で値が大きく異なる。高温では膜沸騰、低温では核沸騰となり、膜沸騰から核沸騰に遷移する温度域５００℃前後では、核沸騰への遷移が始まると一気に核沸騰に移行するため、熱伝達率は、図４のグラフのように、ふたこぶの確率分布を示すことになる。

熱伝達率及び熱伝導率を確率変数として扱った場合、厚板冷却モデルは、以下のように導出される。

物体内の熱伝導と境界で熱伝達は、以下の式（２０）、式（２１）で表される。

ここで、Ｔ_δは表面温度、Ｔ_ｗは冷媒温度であり、αは熱伝達率を表す。また、ｋは熱拡散率を表す。

熱拡散率ｋは、ｋ＝λ／ρｃであり、λは熱伝導率を表し、ρは密度を表し、ｃは比熱を表す。これは、一般的に温度の関数となるが、温度の変動範囲が小さい場合には定数として扱うことができる。この時、式（２１）は、以下の式（２２）とできる。

ここで、式（２１）、式（２２）を差分近似して代数方程式を作る。板厚をＮ等分して、表面温度をＴ_０、Ｔ_Ｎ−１、板内部の温度をＴ_ｉ（１≦ｉ≦Ｎ−２）とすると、ｊ＋１時刻の各点の温度は、ｊ時刻の温度を用いて、式（２３）、式（２４）、式（２５）のように書ける。

ここでは、γ＝ｋΔｔ／（Δｘ）２として、時不変の定数として扱う。そして、λ及びαは、時変とする。

上述の式（２３）、式（２４）、式（２５）をまとめると、以下の式（２６）となる。

ここで、表面温度は観測できる値であるので、以下の式（２７）のように表せる。

以上より、厚板冷却状態空間モデルは、以下の式（２８）、式（２９）で表せる。

ｘ_ｊが、内部温度の推定値であり、ｕ_ｊが、冷媒温度であり、ｙ_ｊが、表面温度であり、Ａ_ｊ、Ｂ_ｊ、Ｃ_ｊが、システムパラメータである。以下に、ｘ_ｊ、ｕ_ｊ、ｙ_ｊ、Ａ_ｊ、Ｂ_ｊ、Ｃ_ｊを示す。

＜動作＞
次に、本実施形態の動作について説明する。ここでは、上述の厚板冷却状態空間モデルを用いて、厚板の内部温度を推定する場合を説明する。

図２は、状態推定装置１の状態推定処理のフローチャートである。

状態推定処理を開始する前に、厚板冷却状態空間モデルは、状態空間モデル記憶部１６２に記憶され、厚板の表面温度等の推定処理に必要なデータは測定値記憶部１６１に記憶されているものとする。

また、推定値記憶部１６３には、今までに推定した推定値が記憶されているものとする。少なくとも、時刻ｋの最適推定値を求める場合には、時刻ｋ−Ｎ−１までの推定値及び分散値が記憶されているものとする。例えば、Ｎ（設計パラメータ）が５の場合は、時刻ｋから６個前の時刻までの推定値が記憶されているものとする。状態推定装置１の状態推定処理では、或る時刻の推定値を算出する際に、１つ前の推定値を用いる。つまり、時刻ｋ−６の推定値を用いて、時刻ｋ−５の推定値を算出し、時刻ｋ−５の推定値を用いて、時刻ｋ−４の推定値を算出することを繰り返して、時刻ｋの推定値を算出する。従って、Ｎは、対象となる状態空間モデルにおいて最適な推定値が求まるまでさかのぼる回数を設定する。尚、最初の推定値を算出する場合は、１つ前の推定値が存在しないため、経験則等から決定した値を推定値として用いる。

フローチャートでは、ｋに１を加算することで、必要な回数の繰り返しが実現される。

まず、ユーザは、入力部１２を用いて、状態推定処理の開始を指示するコマンド等を入力する。

入力部１２を介して、推定処理の開始を指示するコマンドが入力されたことを検知した１１の制御部１１１は、初期設定部１１２に初期値の取得を依頼する。依頼を受けた初期設定部１１２は、出力部１３に初期値の入力を要求するメッセージを表示し、入力部１２を介して、ユーザから設計パラメータＮを取得すると、推定値記憶部１６３から状態初期値ｘ_{ｋ−Ｎ−１}ハット、分散初期値Π_{ｋ−Ｎ−１}を読み出す。初期設定部１１２は、読み出した値及びシステムパラメータＮを、作業メモリに記憶する（ステップＳ１０）。

次に、制御部１１１は、前時刻推定値算出部１１３に前時間の推定値の算出を依頼する。依頼を受けた前時刻推定値算出部１１３は、補助記憶部１６の状態空間モデル記憶部１６２に記憶されている空間モデルを参照し、また、必要に応じて測定値記憶部１６１及び作業メモリからデータを読み出し、上述の式（１３）を用いて、状態の平均値ｘ_ｋ−Ｎオーバーラインを算出する。また、前時刻推定値算出部１１３は、上述の式（１４）を用いて、分散Π_ｋ−Ｎを算出する（ステップＳ１１）。

次に、前時刻推定値算出部１１３は、算出した状態の平均値ｘ_ｋ−Ｎオーバーライン、及び、分散Π_ｋ−Ｎを用いて、上述の式（１２）から、Ｎ個前の状態推定値ｘ_ｋ−Ｎハットを算出し、更に、式（１５）を用いて、ｋ−Ｎからｋ−１までそれぞれの外乱推定値ｗ_ｋ−Ｎハット、・・・、ｗ_ｋ−１ハットを算出する（ステップＳ１２）。尚、厚板冷却状態空間モデルでは、上述の式（２８）に示すように、システムノイズ（外乱）は用いないことから、外乱推定値は算出しない。

前時刻推定値算出部１１３が推定値を算出すると、次に、制御部１１１は、最適推定値算出部１１４に最適推定値の算出を依頼する。依頼を受けた最適推定値算出部１１４は、前時刻推定値算出部１１３が算出した状態推定値ｘ_ｋ−Ｎハット、及び、外乱推定値ｗ_ｋ−Ｎハット、・・・、ｗ_ｋ−１ハットと、測定値記憶部１６１から読出したｋ−Ｎからｋ−１までそれぞれの制御入力ｕ_ｋ−Ｎ、・・・、ｕ_ｋ−１を用いて、上述の式（１１）から、最適推定値ｘ_ｋハットを算出する（ステップＳ１３）。

制御部１１１は、ｋに１を加算し、ステップＳ１１からの処理を繰り返す（ステップＳ１５：Ｎｏ）。すなわち、ｋに１を加算する毎に、終端時刻（推定値を求めたい時刻）からＮ−ｋ個前の状態推定値を求めることになる。言い換えれば、Ｎ個前から、Ｎ−１個前、Ｎ−２個というように、順番に終端時刻に近い時刻の推定値を求めていく。

ｋが終端時刻（ｋの初期値が０の場合は、ｋ＝Ｎ）となったら（ステップＳ１５：Ｙｅｓ）、制御部１１１は、直近に最適推定値算出部１１４が算出した最適推定値ｘ_ｋハットを、出力部１３に表示し（ステップＳ１６）、処理を終了する。
＜実施形態２＞
実施形態１では、確率変数を正規分布にて与えていたが、実施形態２では、一様分布にて与える場合を説明する。カルマンフィルタは正規分布しか扱えなかったが、実施形態の手法では、正規分布に限らず、種々の確率分布を扱うことが可能である。

ここでは、以下の式（３１）、式（３２）、式（３３）で示すパラメータを持つ、上述の式（１）、式（２）で表される線形システムを考える。

確率変数ａ、ｂ、ｃは、以下のような一様分布に従うとする。Ｕ（ｚ｜ａ，ｂ）は、ｚがａ〜ｂの範囲で一定値であり、ｚがａ未満又はｙを超えると０であることを表す。
ｐ（ａ）＝Ｕ（ａ｜−０．５，０．５）
ｐ（ｂ）＝Ｕ（ｂ｜−３，３）
ｐ（ｃ）＝Ｕ（ｃ｜−３，３）
確率変数ａが従う一様分布の確率密度関数を、図５（ａ）に示し、確率変数ｂ、ｃが従う一様分布の確率密度関数を、図５（ｂ）に示す。

図５（ａ）、（ｂ）の確率分布に従って、確率変数ａ、ｂ、ｃのサンプルを１００通り作成する。そのときの１００サンプルでの確率変数ａの値を、図６に示す。横軸はサンプルの番号を表し、縦軸は確率変数ａの値を示す。図中の破線は、全サンプルの平均値を表す。また、この１００サンプルのａに対して、上述の式（３１）で示す行列Ａの固有値の大きさを、図７に示す。式（３１）の行列Ａは２行×２列の行列となっており、２行×２列の行列の固有値は２つ存在する。この２つの固有値双方を最大固有値（実線のグラフ）と最小固有値（破線のグラフ）として図７で図示してる。この図７から、全サンプルにおいて、固有値が１より小さくなっているので、上述の線形システムは安定システムであることが分かる。

また、各サンプルに対して、５０ステップの状態推定（ｘ_１，ｘ_２，ｘ_３，…，ｘ_５０，の５０ステップの状態推定）を行った。以下の式（３４）を用いて、１００サンプルについて算出した誤差ノルムの値を計算した結果を、図８、図９に示す。図８、図９の各グラフにおいて、横軸が、サンプルの番号を表し、縦軸が、誤差ノルムを表す。

図８は、カルマンフィルタを用いて推定した場合の誤差ノルムの値であり、図９は、実施形態の方法を用いて推定した場合の誤差ノルムの値である。図９（ａ）は、設計パラメータＮ＝０のときの誤差ノルムを示し、図９（ｂ）は、設計パラメータＮ＝１のときの誤差ノルムを示し、図９（ｃ）は、設計パラメータＮ＝２のときの誤差ノルムを示し、図９（ｄ）は、設計パラメータＮ＝３のときの誤差ノルムを示す。尚、推定の初期値はｘ_０ハット＝（０，０）^Ｔとしている。また、図８、図９において、実線が各サンプルでの誤差ノルムの値を示す。太線の破線が前サンプルの平均値を示し、細線の破線が、平均値に標準偏差を加減算した値を示す。

図８、図９に示すサンプルの平均値及び標準偏差の具体的な数値を図１０に示す。図１０の表において、「Ｋａｌｍａｎ Ｆｉｌｔｅｒ」は、カルマンフィルタを用いた場合の平均値と標準偏差を示し、「Ｐｒｏｐｏｓｅｄ Ｍｅｔｈｏｄ」は、実施形態の方法を用いた場合の、Ｎの値ごとの平均値と標準偏差を示す。また、「ｄｅｉｆｆｅｒｅｎｃｅ」は、実施形態の方法を用いた場合の、カルマンフィルタと比較した平均値の減少率を示す。

図１１は、実施形態の方法を用いた場合の誤差ノルムの平均値の遷移を示す。横軸が、設計パラメータＮの値を表し、縦軸が、誤差ノルムの平均値を表す。設計パラメータＮの値が大きくなる程、誤差ノルムの平均値が小さくなることが分かる。

図５〜図１１に示すように、確率変数がガウス分布に従う場合と同様に、実施形態の方法では、平均するとカルマンフィルタよりも誤差が小さくなっており、設計パラメータＮの値を大きくすることで、実施形態の方法の性能が向上することが分かる。また、図１０の表より、確率変数がガウス分布に従う場合に比べて、実施形態の方法とカルマンフィルタとの違いが、設計パラメータＮの値を大きくなるにつれて大きくなっていることが分かる。
＜実施形態３＞
実施形態１では、鋼鉄プラントにおいて熱間圧延により厚鋼板を製造する際の冷却工程で、板の表面温度から板内部の温度を推定する場合を例に説明した。その際、熱伝達率α_ｊを確率変数として扱うことして、図４に、遷移領域での熱伝達率の確率分布を示すグラフを示した。

熱伝達率α_ｊを確率変数として扱う以外に、厚鋼板の原料となる溶鋼を保温する取り鍋の保温能力を確率変数として扱うこととしてもよい。図１２に、取り鍋の保温能力の確率分布を示す。図１２のグラフは、横軸が熱伝達率を表し、縦軸が確率を表す。取り鍋の保温能力は、図１２に示すように正規分布ではなく、偏りを持った分布を示す。

また、圧延材の塑性係数を確率変数として扱うこととしてもよい。図１３に、圧延材の塑性係数の確率分布を示す。図１３のグラフは、横軸が塑性係数を表し、縦軸が確率を表す。圧延材は、温度低下とともに硬くなる。そして、加熱炉から出されて圧延機まで搬送される最短時間は設備能力などで決定されるが、待ち時間等の様々な流動的な要件により、加熱炉から圧延機で処理されるまでの時間は、延長されることになる。従って、圧延材の塑性係数は、図１３に示すように、温度低下（硬い）方向に延びたような確率分布で示されることが多くなる。

１ 状態推定装置
１１ 演算処理部
１２ 入力部
１３ 出力部
１６ 補助記憶部
１１１ 制御部
１１２ 初期設定部
１１３ 前時刻推定値算出部
１１４ 最適推定値算出部
１６１ 測定値記憶部
１６２ 状態空間モデル記憶部
１６３ 推定値記憶部

Claims (5)

1. バッチプロセスから観測情報を取得する取得部と、
   前記観測情報から、状態方程式を用いて、前記バッチプロセスにおける所望の時刻の状態値を推定する状態演算部とを備え、
   前記状態方程式の状態行列Ａ、入力行列Ｂ及び出力行列Ｃのうちの少なくとも１つ以上が、時間変数に無関係な或る確率分布に従ってばらつくものであり、状態行列Ａ^ｎ（ｎ≧１：ｎは整数）の期待値を、過去の時刻の状態値の係数行列として用い、あるいは、Ａ^ｎを含む項の期待値として用いる
   ことを特徴とする状態推定装置。
2. 前記状態演算部は、ｗハットをシステムノイズの推定値、ｕを観測情報、Ｅ[ ]を期待値演算子、Ｎを過去の状態値の数とするとき、以下の状態方程式により、時刻ｋの時の状態値ｘハットを算出すること
   

   

   を特徴とする請求項１に記載の状態推定装置。
3. 前記確率分布は、状態値を推定する製品に応じて設定されること
   を特徴とする請求項１に記載の状態推定装置。
4. 観測情報から、所望の時刻の状態値を推定する状態推定装置で用いられる状態推定方法であって、
   バッチプロセスから前記観測情報を取得する取得ステップと、
   前記観測情報から、状態方程式を用いて、前記バッチプロセスにおける所望の時刻の状態値を推定する状態演算ステップとを備え、
   前記状態方程式の状態行列Ａ、入力行列Ｂ及び出力行列Ｃのうちの少なくとも１つ以上が、時間変数に無関係な或る確率分布に従ってばらつくものであり、状態行列Ａ^ｎ（ｎ≧１：ｎは整数）の期待値を、過去の時刻の状態値の係数行列として用い、あるいは、Ａ^ｎを含む項の期待値として用いる
   ことを特徴とする状態推定方法。
5. 観測情報から、所望の時刻の状態値を推定する状態推定装置で用いられる状態推定プログラムであって、
   バッチプロセスから前記観測情報を取得する取得手段と、
   前記観測情報から、状態方程式を用いて、前記バッチプロセスにおける所望の時刻の状態値を推定する状態演算手段としてコンピュータを機能させ、
   前記状態方程式の状態行列Ａ、入力行列Ｂ及び出力行列Ｃのうちの少なくとも１つ以上が、時間変数に無関係な或る確率分布に従ってばらつくものであり、状態行列Ａ^ｎ（ｎ≧１：ｎは整数）の期待値を、過去の時刻の状態値の係数行列として用い、あるいは、Ａ^ｎを含む項の期待値として用いる
   ことを特徴とする状態推定プログラム。

Images