JP7428288B1 - プラント応答推定装置、プラント応答推定方法、及びプログラム - Google Patents

Abstract

【課題】操業中のプラントの運転データから応答モデルを精度良く推定する技術を提供すること。
【解決手段】本開示の一態様によるプラント応答推定装置は、制御対象プラントの運転データに基づいて、前記制御対象プラントの状態を表す状態ベクトルを作成する状態ベクトル作成部と、前記状態ベクトルを用いて、前記状態ベクトルを学習データとするか否かを判定する学習データ選別部と、前記状態ベクトルを学習データとすると判定された場合、前記状態ベクトルを学習データとして用いて、前記制御対象プラントの応答モデルを表す関数のパラメータを更新するパラメータ更新部と、を有する。
【選択図】図２

Description

本開示は、プラント応答推定装置、プラント応答推定方法、及びプログラムに関する。

制御対象の制御量を目標値に追従させることを目的とする制御方式として、モデル予測制御が知られている。また、モデル予測制御で用いられる線形予測モデルのパラメータを推定する手法として、逐次最小２乗法（ＲＬＳ法：Recursive Least Squares法）等が知られている。例えば、特許文献１には、操業中のプラントの運転データを用いて、ＲＬＳ法によりプラント応答モデルのパラメータを推定する技術が記載されている。

特開２０２３－５３５７３号公報

しかしながら、従来技術では、プラント応答モデルのパラメータ更新に用いられる学習データとして、例えば、定常状態における一定値しか含まれない運転データが用いられ、不適切なモデルパラメータが推定されてしまうことがある。一定値や１つの周波数成分しか含まれない運転データ等は、プラントのダイナミクスの情報が欠落しており、学習データとしては不適切だが、操業中に十分生じ得る。このような不適切な運転データでの学習では、プラント応答モデルのパラメータの推定値が真値と乖離してしまう場合がある。

本開示は、上記の点に鑑みてなされたもので、操業中のプラントの運転データから応答モデルを精度良く推定する技術を提供する。

本開示の一態様によるプラント応答推定装置は、制御対象プラントの運転データに基づいて、前記制御対象プラントの状態を表す状態ベクトルを作成する状態ベクトル作成部と、前記状態ベクトルを用いて、前記状態ベクトルを学習データとするか否かを判定する学習データ選別部と、前記状態ベクトルを学習データとすると判定された場合、前記状態ベクトルを学習データとして用いて、前記制御対象プラントの応答モデルを表す関数のパラメータを更新するパラメータ更新部と、を有する。

操業中のプラントの運転データから応答モデルを精度良く推定する技術が提供される。

第一の実施形態に係るプラント応答推定装置のハードウェア構成の一例を示す図である。 第一の実施形態に係るプラント応答推定装置の機能構成の一例を示す図である。 バッファ部の動作の一例を説明するための図である。 第一の実施形態に係る学習データ選別処理の一例を説明するためのフローチャートである。 第一の実施形態に係る拡大モデルパラメータ推定処理の一例を説明するためのフローチャートである。 パラメータ変換部の動作の一例を説明するための図である。 ステップ応答計算部の動作の一例を説明するための図である。 第一の実施形態に係るステップ応答計算処理の一例を説明するためのフローチャートである。 状態ベクトル更新の一例を説明するための図である。 第二の実施形態に係る学習データ選別処理の一例を説明するためのフローチャートである。 第三の実施形態に係る学習データ選別処理の一例を説明するためのフローチャートである。 第四の実施形態に係る学習データ選別処理の一例を説明するためのフローチャートである。 実施例における対象プラントのステップ応答を示す図である。 学習データ選別部を用いない場合のモデルパラメータ推定結果及び制御結果を示す図（その１）である。 学習データ選別部を用いない場合のモデルパラメータ推定結果及び制御結果を示す図（その２）である。 学習データ選別部を用いた場合のモデルパラメータ推定結果及び制御結果を示す図（その１）である。 学習データ選別部を用いた場合のモデルパラメータ推定結果及び制御結果を示す図（その２）である。

以下、本発明の実施形態について説明する。以下の各実施形態では、操業中のプラントの運転データからプラント応答モデルを精度良く推定することができるプラント応答推定装置１０について説明する。なお、操業中とはプラントが正常に稼働し、通常の運転を行っている状態のことであり、例えば、オンライン中、運転中、運用中等と呼ばれてもよい。また、プラントとは１以上の機械、機器、装置等で構成される産業設備のことであり、プラント応答モデルを用いたモデル予測制御によって制御される制御対象である。プラントの具体例としては、例えば、石油化学プラント、食品プラント、鉄鋼プラント、発電プラント等といったものが挙げられるが、これらは一例であって、これらに限られるものではない。また、一般に、これらのプラントの制御システムでは、プラントの制御量、操作量、外乱量といった入出力データが運転データとして扱われる。

ここで、プラントの運転条件等によっては、例えば、定常状態が継続したり、運転データの値が上限値や下限値に偏ってしまったりすることがある。このため、定常状態における運転データや上限値・下限値に偏った運転データを学習データとして用いてプラント応答モデルのパラメータを更新した場合、それらの運転データによってパラメータが収束し、その結果、パラメータの推定値が真値と乖離してしまうことがある。このように、定常状態における運転データや上限値・下限値に偏った運転データ等は、プラント応答モデルの学習データとしては適切でない運転データであるといえる。なお、学習データとは、プラント応答モデルのパラメータ更新に用いられるデータのことである。

そこで、以下の各実施形態では、定常状態における運転データや上限値・下限値に偏った運転データ等といった適切でない運転データは学習データとしては用いずに、それ以外の運転データを学習データとして用いてプラント応答モデルのパラメータを更新する場合について説明する。これにより、以下の各実施形態に係るプラント応答推定装置１０では、精度の良いプラント応答モデルを推定することが可能となる。

［第一の実施形態］
まず、第一の実施形態について説明する。

＜プラント応答推定装置１０のハードウェア構成例＞
第一の実施形態に係るプラント応答推定装置１０のハードウェア構成例を図１に示す。図１に示すように、第一の実施形態に係るプラント応答推定装置１０は、入力装置１１と、表示装置１２と、外部Ｉ／Ｆ１３と、通信Ｉ／Ｆ１４と、プロセッサ１５と、メモリ装置１６とを有する。これらの各ハードウェアは、それぞれがバス１７を介して通信可能に接続される。

入力装置１１は、例えば、キーボード、マウス、タッチパネル、各種物理ボタン等である。表示装置１２は、例えば、ディスプレイ、表示パネル等である。なお、プラント応答推定装置１０は、例えば、入力装置１１及び表示装置１２のうちの少なくとも一方を有していなくてもよい。

外部Ｉ／Ｆ１３は、記録媒体１３ａ等の外部装置とのインタフェースである。記録媒体１３ａとしては、例えば、ＣＤ（Compact Disc）、ＤＶＤ（Digital Versatile Disk）、ＳＤメモリカード（Secure Digital memory card）、ＵＳＢ（Universal Serial Bus）メモリカード等が挙げられる。

通信Ｉ／Ｆ１４は、プラント応答推定装置１０を通信ネットワークに接続するためのインタフェースである。プロセッサ１５は、例えば、ＣＰＵ（Central Processing Unit）等の各種演算装置である。メモリ装置１６は、例えば、ＳＳＤ（Solid State Drive）、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）、フラッシュメモリ等の各種記憶装置である。

なお、図１に示すハードウェア構成は一例であって、プラント応答推定装置１０は、他のハードウェア構成であってもよい。例えば、プラント応答推定装置１０は、複数のプロセッサ１５や複数のメモリ装置１６を有していてもよいし、図示したハードウェア以外の種々のハードウェアを有していてもよい。

＜プラント応答推定装置１０の機能構成例＞
第一の実施形態に係るプラント応答推定装置１０の機能構成例を図２に示す。図２に示すように、第一の実施形態に係るプラント応答推定装置１０は、モデルパラメータ推定部１０１と、ステップ応答計算部１０２とを有する。これら各部は、例えば、プラント応答推定装置１０にインストールされた１以上のプログラムが、プロセッサ１５等に実行させる処理により実現される。

モデルパラメータ推定部１０１は、サンプリング周期Δ毎に、制御対象のプラント（又はその運転状態を計測するセンサ等の機器）から運転データの観測値（制御量ｙ、操作量ｕ及び外乱ｖ）を受信し、そのサンプリング周期Δに応じて、プラント応答モデルのモデルパラメータを推定する。なお、運転データは、例えば、計測データや観測データ等と呼ばれてもよい。

ここで、時刻をｔとすれば制御量ｙ、操作量ｕ及び外乱ｖはそれぞれｙ（ｔ）、ｕ（ｔ）及びｖ（ｔ）と表され、各サンプリング時刻ｔ_ｋ（ｋはサンプリング時刻を表すインデックスであり、定周期でインクリメントされる。）に関してｔ_ｋ＋１－ｔ_ｋ＝Δが成り立つ。なお、ｋは０以上の整数を取り得るものとし、またサンプリング時刻ｔ_ｋはｔ_１から開始するものとする。以下では、特に断らない限り、各サンプリング時刻ｔ_ｋをそのインデックスｋと同一視し、ｙ（ｋ）、ｕ（ｋ）及びｖ（ｋ）とも表すことにする。

また、以下では、プラント応答モデルは、パラメータθを持つプラント応答関数Ｓ_θ（ｔ）で表されるものとして、パラメータθをモデルパラメータということにする。プラント応答モデルとしては様々なプラント応答関数Ｓ_θで表現されるモデルを採用することが可能であるが、以下では、主に、ＡＲＭＡＸモデルといった多項式モデルを想定する。プラント応答モデルがＡＲＭＡＸモデルである場合、モデルパラメータθはＡＲＭＡＸモデルの係数となる。なお、プラント応答モデルとしてＡＲＭＡＸモデル以外の様々なモデルが採用可能であることは言うまでもない。

更に、以下では、ＡＲＭＡＸモデルで表されたプラント応答モデルに対して、制御対象プラントの定常状態を表すオフセット項（定数項）を追加したモデルを拡大プラント応答モデルと呼ぶことにする。

ステップ応答計算部１０２は、モデルパラメータ推定部１０１によって推定されたモデルパラメータ（以下、モデルパラメータ推定値ともいう。）θを用いて、与えられた時刻ｔにおけるプラントのステップ応答の推定値（以下、ステップ応答推定値ともいう。）Ｓ_θ（ｔ）を計算する。なお、ステップ応答とは、操作量としてステップ信号がプラントに印加されたときの制御量のことである。

ただし、ステップ応答計算部１０２はモデルパラメータ推定値θの利用方法の１つであって、本実施形態に係るプラント応答推定装置１０は、モデルパラメータ推定部１０１のみを有する構成であってもよい。また、本実施形態に係るプラント応答推定装置１０は、ステップ応答以外のプラント応答を計算する機能部を有する構成であってもよい。例えば、インパルス応答を計算する「インパルス応答計算部」を有する構成であってもよいし、ランプ応答を計算する「ランプ応答計算部」を有する構成であってもよい。これら以外にも、その他、本実施形態の主旨を逸脱しない範囲で様々な利用方法に対して適用可能である。

ここで、モデルパラメータ推定部１０１には、バッファ部１１１と、状態ベクトル変換部１１２と、学習データ選別部１１３と、逐次推定計算部１１４と、パラメータ変換部１１５とが含まれる。

バッファ部１１１は、或る所定の期間における制御量ｙ（ｋ）、操作量ｕ（ｋ）及び外乱ｖ（ｋ）をメモリ装置１６に蓄積（バッファ）する。状態ベクトル変換部１１２は、メモリ装置１６にバッファされている制御量ｙ（ｋ）、操作量ｕ（ｋ）及び外乱ｖ（ｋ）の再サンプリングを行って、それら再サンプリングされた制御量ｙ（ｋ）、操作量ｕ（ｋ）及び外乱ｖ（ｋ）が含まれる状態ベクトルξ（ｋ）を作成する。学習データ選別部１１３は、状態ベクトルξ（ｋ）を用いて、再サンプリングされた制御量ｙ（ｋ）、操作量ｕ（ｋ）及び外乱ｖ（ｋ）を学習データとして用いるか否かを判定し、学習データとして用いる場合は更新フラグをＯＮ、学習データとして用いない場合は更新フラグをＯＦＦとする。逐次推定計算部１１４は、更新フラグがＯＮである場合、状態ベクトルξ（ｋ）を用いて、拡大プラント応答モデルのモデルパラメータである拡大モデルパラメータを推定する。パラメータ変換部１１５は、拡大モデルパラメータの推定値をモデルパラメータ推定値θに変換する。なお、状態ベクトルとは制御の分野で用いられる用語であり、制御対象やそれを含む系（システム）の状態を表すベクトルである。

＜バッファ部１１１の動作＞
サンプリング時刻ｔ_ｋにおける制御量バッファをＹ（ｋ）、操作量バッファをＵ（ｋ）、外乱バッファをＶ（ｋ）として、これらの各バッファは以下のベクトルで表されるものとする。

すなわち、制御量バッファＹ（ｋ）にはｋ－Ｂ_１からｋまでの制御量ｙ、操作量バッファＵ（ｋ）にはｋ－Ｂ_２からｋまでの操作量ｕ、外乱バッファＶ（ｋ）にはｋ－Ｂ_３からｋまでの外乱ｖがそれぞれ格納されているものとする。ここで、Ｂ_１、Ｂ_２及びＢ_３はそれぞれ０以上の整数であり、制御量バッファ、操作量バッファ及び外乱バッファの大きさを決めるパラメータである。これらのＢ_１、Ｂ_２及びＢ_３はメモリ装置１６のサイズ等に応じて、適宜、その値が決定される。

このとき、バッファ部１１１は、制御量ｙ（ｋ）、操作量ｕ（ｋ）及び外乱ｖ（ｋ）を受信すると、図３に示すように、制御量ｙ（ｋ）と制御量バッファＹ（ｋ－１）から制御量バッファＹ（ｋ）、操作量ｕ（ｋ）と操作量バッファＵ（ｋ－１）から操作量バッファＵ（ｋ）、外乱ｖ（ｋ）と外乱バッファＶ（ｋ－１）から外乱バッファＶ（ｋ）にそれぞれ更新する。

具体的には、バッファ部１１１は、制御量バッファＹ（ｋ－１）に格納されているｙ（ｋ－Ｂ_１－１）を削除した上で、新たに観測された制御量ｙ（ｋ）を格納することで、ｙ（ｋ－Ｂ_１）からｙ（ｋ）までの制御量が格納された制御量バッファＹ（ｋ）に更新する。同様に、バッファ部１１１は、操作量バッファＵ（ｋ－１）に格納されている操作量ｕ（ｋ－Ｂ_２－１）を削除した上で、新たに観測された操作量ｕ（ｋ）を格納することで、ｕ（ｋ－Ｂ_２）からｕ（ｋ）までの操作量が格納された操作量バッファＵ（ｋ）に更新する。同様に、バッファ部１１１は、外乱バッファＶ（ｋ－１）に格納されている外乱ｖ（ｋ－Ｂ_３－１）を削除した上で、新たに観測された外乱ｖ（ｋ）を格納することで、ｖ（ｋ－Ｂ_３）からｖ（ｋ）までの外乱が格納された外乱バッファＶ（ｋ）に更新する。

なお、制御量バッファＹ（０）、操作量バッファＵ（０）、外乱バッファＶ（０）は適切に初期化（例えば、すべて０等に初期化）されているものとする。また、上記では外乱ｖがある場合を想定して説明したが、外乱ｖがない場合は外乱バッファＶ（ｋ）がなくてもよい。

＜状態ベクトル変換部１１２の動作＞
再サンプリング周期をＤとする。このとき、状態ベクトル変換部１１２は、制御量バッファＹ（ｋ）、操作量バッファＵ（ｋ）及び外乱バッファＶ（ｋ）が更新されると、これらの各バッファＹ（ｋ）、Ｕ（ｋ）及びＶ（ｋ）から再サンプリング周期Ｄで再サンプリングを行って、以下の再サンプリング制御量ベクトルＹ_Ｄ（ｋ）、再サンプリング操作量ベクトルＵ_Ｄ（ｋ）及び再サンプリング外乱ベクトルＶ_Ｄ（ｋ）をそれぞれ作成する。

ここで、Ｎ、Ｍ及びＬはプラント応答モデルに応じて決定されるパラメータ（具体的には、ＮはＡＲＭＡＸモデルの制御量ｙに関する項の係数の数、Ｍは操作量ｕに関する項の係数の数、Ｌは外乱ｖに関する項の係数の数）である。また、再サンプリング制御量ベクトルＹ_Ｄ（ｋ）ではＡＲＭＡＸモデルにならい、現在のサンプリング時刻を表す要素ｙ（ｋ）がスキップされる（つまり、ｙ（ｋ）は再サンプリングされない。）。

なお、一般に、Ｎ、Ｍ及びＬの値は大きい方が多様な表現が可能で、プラント応答の高精度な予測が期待できるが、モデルパラメータθの推定のために多くの計算資源やメモリ量が必要となる。

そして、状態ベクトル変換部１１２は、再サンプリング制御量ベクトルＹ_Ｄ（ｋ）、再サンプリング操作量ベクトルＵ_Ｄ（ｋ）及び再サンプリング外乱ベクトルＶ_Ｄ（ｋ）から以下により状態ベクトルξ（ｋ）を作成する。

すなわち、状態ベクトルξ（ｋ）とは、Ｙ_Ｄ（ｋ）、Ｕ_Ｄ（ｋ）及びＶ_Ｄ（ｋ）の要素に対して、オフセット項を表す要素として「１」を追加したベクトルである。ここで、オフセット項とは、上述したように、ＡＲＭＡＸモデルで定常状態を表す項として出現する定数項のことである。これにより、再サンプリング制御量ベクトルＹ_Ｄ（ｋ）、再サンプリング操作量ベクトルＵ_Ｄ（ｋ）及び再サンプリング外乱ベクトルＶ_Ｄ（ｋ）が状態ベクトルξ（ｋ）に変換されたことになる。この状態ベクトルξ（ｋ）は「拡大状態ベクトル」等と呼ばれてもよい。

なお、状態ベクトルξ（０）は適切な値に初期化（例えば、零ベクトル等に初期化）されているものとする。また、上記では外乱ｖがある場合を想定して説明したが、外乱ｖがない場合は再サンプリング外乱ベクトルＶ_Ｄ（ｋ）がなくてもよい。

＜学習データ選別部１１３の動作＞
学習データ選別部１１３は、状態ベクトルξ（ｋ）が作成されると、この状態ベクトルξ（ｋ）を用いて、再サンプリング制御量ベクトルＹ_Ｄ（ｋ）、再サンプリング操作量ベクトルＵ_Ｄ（ｋ）及び再サンプリング外乱ベクトルＶ_Ｄ（ｋ）を学習データとして用いるか否かを判定した上で、その判定結果に応じて更新フラグをＯＮ又はＯＦＦにする。これにより、学習データとして用いる運転データ（の観測値）が選別される。

以下、或るｋに関して学習データとして用いる運転データを選別する処理（第一の実施形態に係る学習データ選別処理）について、図４を参照しながら説明する。ただし、ｋ＝０である場合は更新フラグ＝ＯＮとして処理を終了する。

ステップＳ１０１：まず、学習データ選別部１１３は、状態ベクトルξ（ｋ）と１時刻前の状態ベクトルξ（ｋ－１）との差のノルムｄ（ｋ）を計算する。すなわち、学習データ選別部１１３は、ｄ（ｋ）＝||ξ（ｋ）－ξ（ｋ－１）||_ｐを計算する。ここで、ｐは、ｐ＝１、ｐ＝２又はｐ＝∞のいずれかであり、||・||_１はＬ１ノルム、||・||_２はＬ２ノルム、||・||_∞はＬ∞ノルムをそれぞれ表す。以下、ｄ（ｋ）をベクトル差ノルムともいう。

ステップＳ１０２：次に、学習データ選別部１１３は、ベクトル差ノルムｄ（ｋ）が予め設定された閾値γ以上であるか否かを判定する。

ｄ（ｋ）≧γであると判定した場合（ステップＳ１０２でＹＥＳ）、学習データ選別部１１３は、ステップＳ１０３に進む。一方で、ｄ（ｋ）≧γであると判定しなかった場合（ステップＳ１０２でＮＯ）、学習データ選別部１１３は、ステップＳ１０４に進む。

ステップＳ１０３：学習データ選別部１１３は、更新フラグをＯＮにする。ここで、更新フラグとは、再サンプリング操作量ベクトルＵ_Ｄ（ｋ）及び再サンプリング外乱ベクトルＶ_Ｄ（ｋ）を学習データとして用いるか否かを表すフラグのことである。更新フラグがＯＮである場合、再サンプリング操作量ベクトルＵ_Ｄ（ｋ）及び再サンプリング外乱ベクトルＶ_Ｄ（ｋ）が学習データとして用いられて拡大モデルパラメータの更新が行われる。一方で、更新フラグがＯＦＦである場合、再サンプリング操作量ベクトルＵ_Ｄ（ｋ）及び再サンプリング外乱ベクトルＶ_Ｄ（ｋ）は学習データとして用いられず、拡大モデルパラメータの更新は行われない。

ステップＳ１０４：学習データ選別部１１３は、更新フラグをＯＦＦにする。

＜逐次推定計算部１１４の動作＞
以下では、拡大モデルパラメータをθ_ｅと表す。逐次推定計算部１１４は、学習データ選別部１１３によって更新フラグがＯＮ又はＯＦＦに設定されると、状態ベクトルξ（ｋ）を用いて、拡大モデルパラメータθ_ｅの値を推定する。すなわち、逐次推定計算部１１４は、サンプリング時刻ｔ_ｋ毎に、逐次的に拡大モデルパラメータθ_ｅの値を推定する。

以下、或るｋに関して拡大モデルパラメータθ_ｅの値を推定する処理（第一の実施形態に係る拡大モデルパラメータ推定処理）について、図５を参照しながら説明する。なお、以下では、サンプリング時刻ｔ_ｋにおける拡大モデルパラメータθ_ｅの推定値をθ_ｅ（ｋ）と表す。

ステップＳ２０１：まず、逐次推定計算部１１４は、更新フラグがＯＮであるか否かを判定する。

更新フラグがＯＮであると判定した場合（ステップＳ２０１でＹＥＳ）、逐次推定計算部１１４は、ステップＳ２０２に進む。一方で、更新フラグがＯＮであると判定しなかった場合（ステップＳ２０１でＮＯ）、逐次推定計算部１１４は、拡大モデルパラメータ推定処理を終了する。すなわち、更新フラグがＯＦＦである場合、逐次推定計算部１１４は、拡大モデルパラメータθ_ｅの値を推定及び更新せずに処理を終了する。

ステップＳ２０２：逐次推定計算部１１４は、拡大モデルパラメータθ_ｅと共分散行列Ｐとを初期化するか否かを判定する。ここで、初期化すると判定される場合としては、例えば、初回計算時（つまり、ｋ＝１のとき）、ユーザ等により初期化指示が行われたとき等が挙げられる。

拡大モデルパラメータθ_ｅと共分散行列Ｐとを初期化すると判定した場合（ステップＳ２０２でＹＥＳ）、逐次推定計算部１１４は、ステップＳ２０３に進む。一方で、拡大モデルパラメータθ_ｅと共分散行列Ｐとを初期化すると判定しなかった場合（ステップＳ２０２でＮＯ）、逐次推定計算部１１４は、ステップＳ２０４に進む。

ステップＳ２０３：逐次推定計算部１１４は、サンプリング時刻ｔ_ｋのインデックスｋをｋ＝１に初期化すると共に、θ_ｅ（０）＝θ_０及びＰ（０）＝Ｉと初期化する。ここで、θ_０は予め設定された拡大モデルパラメータの初期値、Ｉは予め設定された任意の行列（例えば、単位行列等）である。

ステップＳ２０４：逐次推定計算部１１４は、拡大モデルパラメータ推定値θ_ｅ（ｋ－１）と状態ベクトルξ（ｋ）と制御量ｙ（ｋ）とを用いて、予測誤差ε（ｋ）を計算する。予測誤差ε（ｋ）は、例えば、ε（ｋ）＝ｙ（ｋ）－ξ（ｋ）^Τθ_ｅ（ｋ－１）により計算される。なお、Τは転置を表す。

ステップＳ２０５：逐次推定計算部１１４は、共分散行列Ｐ（ｋ－１）を以下により更新して共分散行列Ｐ（ｋ）を得る。

ここで、λ＞０は調整係数であり、予め設定された値である。調整係数λはその値に応じて以下の２つの使い方をすることができる。

１つ目は、０＜λ≦１とすることにより、過去の運転データを忘却するための忘却係数として用いることができる。λが０に近いほど急激に過去の運転データの影響が減少し、１に近いほど過去の運転データの影響が保持される。

２つ目は、λ＞１とすることにより、パラメータ更新を加速させるための係数として用いることができる。λの値が大きいほどパラメータの収束を加速させることができる。

ステップＳ２０６：そして、逐次推定計算部１１４は、拡大モデルパラメータ推定値θ_ｅ（ｋ－１）を以下により更新して拡大モデルパラメータ推定値θ_ｅ（ｋ）を得る。

これにより、サンプリング時刻ｔ_ｋにおける拡大モデルパラメータ推定値θ_ｅ（ｋ）が得られる。

＜パラメータ変換部１１５の動作＞
拡大モデルパラメータ推定値θ_ｅ（ｋ）を以下で表す。

ここで、θ_Ｙ（ｋ）はサンプリング時刻ｔ_ｋにおけるＡＲＭＡＸモデルの制御量ｙに関する項の係数を要素とするＮ次元ベクトル、θ_Ｕ（ｋ）は操作量ｕに関する項の係数を要素とするＭ次元ベクトル、θ_Ｖ（ｋ）は外乱ｖに関する項の係数を要素とするＬ次元ベクトル、θ_Ｃ（ｋ）はオフセット項を表すスカラー値である。

このとき、パラメータ変換部１１５は、拡大モデルパラメータ推定値θ_ｅ（ｋ）が得られると、図６に示すように、θ_Ｃ（ｋ）を除いたベクトルに変換することで、モデルパラメータ推定値θ（ｋ）を得る。これにより、以下のモデルパラメータ推定値が得られる。

このように、制御対象プラントのステップ応答を推定する際には、拡大モデルパラメータ推定値θ_ｅ（ｋ）からθ_Ｃ（ｋ）を除くことで、モデルパラメータ推定値θ（ｋ）が得られる。この理由については、例えば、特許文献１等を参照されたい。

＜ステップ応答計算部１０２の動作＞
以下、ステップ応答の推定に用いられるモデルパラメータ推定値（つまり、プラント応答関数Ｓ_θに設定されるモデルパラメータ推定値）θ＝θ（ｋ）のことを「モデルパラメータ設定値θ」ともいう。また、以下では、簡単のため、ステップ信号として単位ステップ信号を想定する。

ステップ応答計算部１０２は、図７に示すように、モデルパラメータ設定値θと時刻ｔとが与えられると、初期時刻０から時間ｔ経過後の時刻ｔにおける単位ステップ応答Ｓ_θ（ｔ）を計算する。なお、単位ステップ応答とは、操作量ｕとして単位ステップ信号を印加した場合における応答（つまり、制御対象プラントのプラント応答モデルの出力）のことである。

或る時刻ｔと或るモデルパラメータ設定値θとが与えられたときに、ステップ応答Ｓ_θ（ｔ）を計算する処理（第一の実施形態に係るステップ応答計算処理）について、図８を参照しながら説明する。なお、以下では、インデックスｋ'における状態ベクトルφ（ｋ'）を以下で表すものとする。

なお、インデックスｋ'は、サンプリング時刻ｔ_ｋのインデックスｋと同じ値を取り得る変数であるが、本処理の中でのみ利用され、インデックスｋとは独立に値が更新されることに留意されたい。

ステップＳ３０１：ステップ応答計算部１０２は、本処理の中でのみ利用する時刻を表すインデックスをτとして、τ＝０、ｋ'＝０と初期化すると共に、状態ベクトルφ（０）を以下のように初期化する。

すなわち、ｕ（０）のみ１、それ以外の要素は０と状態ベクトルφ（０）を初期化する。

ステップＳ３０２：ステップ応答計算部１０２は、ｙ（ｋ'）＝φ（ｋ'）^Τθにより制御量予測値ｙ（ｋ'）を計算する。

ステップＳ３０３：ステップ応答計算部１０２は、制御量予測値ｙ（ｋ'）を用いて、状態ベクトルφ（ｋ'）を、次のインデックスｋ'＋１における状態ベクトルφ（ｋ'＋１）に更新する。このとき、ステップ応答計算部１０２は、図９に示すように、状態ベクトルφ（ｋ'＋１）のｙ（ｋ'）には上記のステップＳ３０２で計算した制御量予測値ｙ（ｋ'）を設定し、ｙ（ｋ'－Ｎ＋１）～ｙ（ｋ'－１）には状態ベクトルφ（ｋ'）と同じ値を設定する。また、状態ベクトルφ（ｋ'＋１）のｕ（ｋ'＋１）には１を設定し、ｕ（ｋ'－Ｍ＋１）～ｕ（ｋ'）には状態ベクトルφ（ｋ'）と同じ値を設定する。更に、状態ベクトルφ（ｋ'＋１）のｖ（ｋ'＋１）には０を設定し、ｖ（ｋ'－Ｌ＋１）～ｖ（ｋ'）には状態ベクトルφ（ｋ'）と同じ値を設定する。

ステップＳ３０４：ステップ応答計算部１０２は、時刻τをτ＋Δに更新すると共に、インデックスｋ'をｋ'＋１に更新する。

ステップＳ３０５：ステップ応答計算部１０２は、τ≧ｔであるか否かを判定する。そして、τ≧ｔであると判定されなかった場合（ステップＳ３０５でＮＯ）、ステップ応答計算部１０２は、ステップＳ３０２に戻る。これにより、τ≧ｔとなるまで、ステップＳ３０２～ステップＳ３０４が繰り返し実行される。

一方で、τ≧ｔであると判定された場合（ステップＳ３０５でＹＥＳ）、ステップ応答計算部１０２は、ステップ応答計算処理を終了する。これにより、最終的に計算された制御量予測値ｙ（ｋ'）が単位ステップ応答Ｓ_θ（ｔ）として得られる（つまり、Ｓ_θ（ｔ）＝ｙ（ｋ）が、プラント応答モデルの単位ステップ応答として得られる。）。

［第二の実施形態］
次に、第二の実施形態について説明する。第二の実施形態では、学習データ選別処理が第一の実施形態と異なる。なお、第二の実施形態で特に言及しない事項については、第一の実施形態と同様としてよい。

そこで、以下、或るｋに関して学習データとして用いる運転データを選別する処理（第二の実施形態に係る学習データ選別処理）について、図１０を参照しながら説明する。ただし、ｋ＝０である場合は更新フラグ＝ＯＮとして処理を終了する。

ステップＳ４０１：まず、学習データ選別部１１３は、図４のステップＳ１０１と同様に、状態ベクトルξ（ｋ）と１時刻前の状態ベクトルξ（ｋ－１）との差のノルムｄ（ｋ）＝||ξ（ｋ）－ξ（ｋ－１）||_ｐをベクトル差ノルムとして計算する。

ステップＳ４０２：次に、学習データ選別部１１３は、状態ベクトルξ（ｋ）を用いて、閾値γ（ｋ）を計算する。すなわち、学習データ選別部１１３は、以下により閾値γ（ｋ）を計算する。

ここで、Ｅ_ξ（ｋ）は状態ベクトルξ（１），・・・，ξ（ｋ）の平均値、Ｖ_ξ（ｋ）は状態ベクトルξ（１），・・・，ξ（ｋ）の分散、γ（ｋ）は状態ベクトルξ（１），・・・，ξ（ｋ）の標準偏差のα倍を表す。また、αは調整係数であり、予め設定された値である。

ステップＳ４０３：次に、学習データ選別部１１３は、ベクトル差ノルムｄ（ｋ）が閾値γ（ｋ）以上であるか否かを判定する。

ｄ（ｋ）≧γ（ｋ）であると判定した場合（ステップＳ４０３でＹＥＳ）、学習データ選別部１１３は、ステップＳ４０４に進む。一方で、ｄ（ｋ）≧γ（ｋ）であると判定しなかった場合（ステップＳ４０３でＮＯ）、学習データ選別部１１３は、ステップＳ４０５に進む。

ステップＳ４０４：学習データ選別部１１３は、図４のステップＳ１０３と同様に、更新フラグをＯＮにする。

ステップＳ４０５：学習データ選別部１１３は、図４のステップＳ１０４と同様に、更新フラグをＯＦＦにする。

［第三の実施形態］
次に、第三の実施形態について説明する。第三の実施形態では、学習データ選別処理が第一の実施形態と異なる。また、閾値γ（ｋ）の計算方法が第二の実施形態と異なる。なお、第三の実施形態で特に言及しない事項については、第一の実施形態と同様としてよい。

そこで、以下、或るｋに関して学習データとして用いる運転データを選別する処理（第三の実施形態に係る学習データ選別処理）について、図１１を参照しながら説明する。ただし、ｋ＝０である場合は更新フラグ＝ＯＮとして処理を終了する。

ステップＳ５０１：まず、学習データ選別部１１３は、ｋ＝１であるか否かを判定する。

ｋ＝１であると判定した場合（ステップＳ５０１でＹＥＳ）、学習データ選別部１１３は、ステップＳ５０２に進む。一方で、ｋ＝１であると判定しなかった場合（ステップＳ５０１でＮＯ）、学習データ選別部１１３は、ステップＳ５０３に進む。

ステップＳ５０２：学習データ選別部１１３は、Ｅ_ξ（０）＝０、Ｓ_ξ（０）＝０と初期化する。ここで、Ｅ_ξ（０）は状態ベクトルのノルムの平均値の初期値、Ｓ_ξ（０）は状態ベクトルのノルムの２乗の平均値の初期値である。

ステップＳ５０３：学習データ選別部１１３は、図４のステップＳ１０１と同様に、状態ベクトルξ（ｋ）と状態ベクトルξ（ｋ－１）との差のノルムｄ（ｋ）＝||ξ（ｋ）－ξ（ｋ－１）||_ｐをベクトル差ノルムとして計算する。

ステップＳ５０４：次に、学習データ選別部１１３は、状態ベクトルξ（ｋ）を用いて、閾値γ（ｋ）を計算する。すなわち、学習データ選別部１１３は、以下により閾値γ（ｋ）を計算する。

ここで、Ｅ_ξ（ｋ）は状態ベクトルξ（１），・・・，ξ（ｋ）のノルムの平均値、Ｓ_ξ（ｋ）は状態ベクトルξ（１），・・・，ξ（ｋ）のノルムの２乗の平均値、Ｖ_ξ（ｋ）は状態ベクトルξ（１），・・・，ξ（ｋ）のノルムの分散、γ（ｋ）は状態ベクトルξ（１），・・・，ξ（ｋ）のノルムの標準偏差のα倍を表す。また、αは調整係数であり、予め設定された値である。

ステップＳ５０５：次に、学習データ選別部１１３は、ベクトル差ノルムｄ（ｋ）が閾値γ（ｋ）以上であるか否かを判定する。

ｄ（ｋ）≧γ（ｋ）であると判定した場合（ステップＳ５０５でＹＥＳ）、学習データ選別部１１３は、ステップＳ５０６に進む。一方で、ｄ（ｋ）≧γ（ｋ）であると判定しなかった場合（ステップＳ５０５でＮＯ）、学習データ選別部１１３は、ステップＳ５０７に進む。

ステップＳ５０６：学習データ選別部１１３は、図４のステップＳ１０３と同様に、更新フラグをＯＮにする。

ステップＳ５０７：学習データ選別部１１３は、図４のステップＳ１０４と同様に、更新フラグをＯＦＦにする。

［第四の実施形態］
次に、第四の実施形態について説明する。第四の実施形態では、学習データ選別処理が第一の実施形態～第三の実施形態と異なる。なお、第四の実施形態で特に言及しない事項については、第一の実施形態と同様としてよい。

そこで、以下、或るｋに関して学習データとして用いる運転データを選別する処理（第四の実施形態に係る学習データ選別処理）について、図１２を参照しながら説明する。ただし、ｋ＝０である場合は更新フラグ＝ＯＮとして処理を終了する。

ステップＳ６０１：まず、学習データ選別部１１３は、ｋ＝１であるか否かを判定する。

ｋ＝１であると判定した場合（ステップＳ６０１でＹＥＳ）、学習データ選別部１１３は、ステップＳ６０２に進む。一方で、ｋ＝１であると判定しなかった場合（ステップＳ６０１でＮＯ）、学習データ選別部１１３は、ステップＳ６０３に進む。

ステップＳ６０２：学習データ選別部１１３は、Ｅ_ｙ（０）＝０、Ｓ_ｙ（０）＝０と初期化する。ここで、Ｅ_ｙ（０）は制御量の平均値の初期値、Ｓ_ｙ（０）は制御量の２乗の平均値の初期値である。

ステップＳ６０３：学習データ選別部１１３は、状態ベクトルξ（ｋ）を用いて、予測誤差（の絶対値）ｄ（ｋ）を計算する。予測誤差ｄ（ｋ）は、図５のステップＳ２０４と同様に、ｄ（ｋ）＝｜ｙ（ｋ）－ξ（ｋ）^Τθ_ｅ（ｋ－１）｜により計算される。

ステップＳ６０４：次に、学習データ選別部１１３は、制御量ｙ（ｋ）を用いて、閾値γ（ｋ）を計算する。すなわち、学習データ選別部１１３は、以下により閾値γ（ｋ）を計算する。

ここで、Ｅ_ｙ（ｋ）は制御量ｙ（１），・・・，ｙ（ｋ）の平均値、Ｓ_ｙ（ｋ）は制御量ｙ（１），・・・，ｙ（ｋ）の２乗の平均値、Ｖ_ｙ（ｋ）は制御量ｙ（１），・・・，ｙ（ｋ）の分散、γ（ｋ）は制御量ｙ（１），・・・，ｙ（ｋ）の標準偏差のα倍を表す。また、αは調整係数であり、予め設定された値である。

ステップＳ６０５：次に、学習データ選別部１１３は、予測誤差ｄ（ｋ）が閾値γ（ｋ）以上であるか否かを判定する。

ｄ（ｋ）≧γ（ｋ）であると判定した場合（ステップＳ６０５でＹＥＳ）、学習データ選別部１１３は、ステップＳ６０６に進む。一方で、ｄ（ｋ）≧γ（ｋ）であると判定しなかった場合（ステップＳ６０５でＮＯ）、学習データ選別部１１３は、ステップＳ６０７に進む。

ステップＳ６０６：学習データ選別部１１３は、図４のステップＳ１０３と同様に、更新フラグをＯＮにする。

ステップＳ６０７：学習データ選別部１１３は、図４のステップＳ１０４と同様に、更新フラグをＯＦＦにする。

［実施例］
以下、第一の実施形態に係るプラント応答推定装置１０の実施例について説明する。

本実施例では、制御対象のプラントの拡大プラント応答モデルは次式で表されるものとする。

ｙ（ｋ）＝ａ_１ｙ（ｋ－１）＋ａ_２ｙ（ｋ－２）＋ｂ_０ｕ（ｋ）＋ｂ_１ｕ（ｋ－１）＋ｂ_２ｕ（ｋ－２）＋ｃ_０
すなわち、制御量ｙに関しては２次元、操作量ｕに関しては３次元の要素を有するものとする。なお、本実施例では外乱ｖがないものとする。

このとき、状態ベクトルξ（ｋ）及び拡大モデルパラメータθ_ｅ（ｋ）は以下で表される。

また、制御量ｙの推定値はｙ_ｅｓｔ（ｋ）＝ξ（ｋ）^Τθ_ｅ（ｋ－１）で表される。

更に、本実施例では、学習データ選別部１１３がベクトル差ノルムｄ（ｋ）を計算するときのノルムとしてはＬ１ノルムを用いると共に閾値はγ＝０．１とした。また、本実施例における制御対象プラントのステップ応答としては、図１３に示すステップ応答を想定する。

以上の設定の下で、学習データ選別部１１３を用いない場合と学習データ選別部１１３を用いた場合とを比較した。

図１４は、学習データ選別部１１３を用いない場合で調整係数をλ＝０．９９９としたときのモデルパラメータ推定結果及び制御結果である。また、図１５は、学習データ選別部１１３を用いない場合で調整係数をλ＝１．１としたときのモデルパラメータ推定結果及び制御結果である。

図１４及び図１５はいずれも時刻０から１０００までは制御量ｙと操作量ｕが一定であるため、モデルパラメータも時刻０付近で初期値から変動した後、時刻１０００まで一定となっている。

時刻１０００以降に関しては、図１４では、モデルパラメータは変化し、時刻４０００付近まで変化がみられる。一方で、図１５では、モデルパラメータは変化していない。これは、調整係数λが大きい場合、共分散行列の更新が早まり、ゼロ付近になるためである。

制御応答に関しては、図１４では、モデルパラメータが変化したことにより、制御量が目標値に追従できている。一方で、図１５では、モデルパラメータが初期に推定した値から変化していないため、モデルが正しく推定されておらず、制御量が目標値に追従できずに応答が不安定になっている。

図１６は、学習データ選別部１１３を用いた場合で調整係数をλ＝０．９９９としたときのモデルパラメータ推定結果及び制御結果である。また、図１７は、学習データ選別部１１３を用いた場合で調整係数をλ＝１．１としたときのモデルパラメータ推定結果及び制御結果である。

図１６及び図１７はいずれも時刻０から１０００までは制御量ｙと操作量ｕが一定であるため、モデルパラメータも時刻０付近で初期値から変動した後、時刻１０００まで一定となっている。

時刻１０００以降に関しては、図１６では、モデルパラメータは変化し、時刻１５００付近で収束している。一方で、図１７では、図１５と異なり、モデルパラメータは変化し、時刻１２００付近で収束している。これは、学習データ選別部１１３の働きによって時刻０～１０００までの定常状態の運転データが更新に用いられないため共分散行列が定常状態でもゼロに収束せず、また時刻１０００付近で制御量ｙと操作量ｕが変化した運転データは更新に用いられるためパラメータ更新が加速され真値に近くなり、収束が早まったためである。

制御応答に関しては、図１６及び図１７のいずれもモデルパラメータが比較的正しく推定されているため、制御量が目標値に追従できている。

よって、調整係数λを１より大きくしてモデルパラメータの収束を早めた際に、学習データ選別部１１３の働きにより、定常状態における運転データでパラメータが収束してしまう事態を回避できているといえる。

以上のように、学習データ選別部１１３により、例えば、定常状態が継続するだけ等といった運転データに十分な情報が含まれない場合にはパラメータ更新が停止され、制御量ｙや操作量ｕが変化する場合にはパラメータ更新が行われる。これにより、例えば、定常状態で真値とは遠い局所的な値にパラメータが不用意に収束してしまう事態を回避しつつ、学習データを選別しない場合よりも真値に近い値に素早く収束する、という利点が生じている。

なお、上記の各実施形態では、主に、プラント応答モデルのパラメータを推定する場合について説明したが、これに限られず、例えば、上記の各実施形態に係るプラント応答推定装置１０は、当該パラメータを設定したプラント応答モデルにより制御対象プラントを実際に制御する制御装置として機能してもよい。このとき、上記の各実施形態に係るプラント応答推定装置１０は、モデル予測制御等といった既知の制御手法を用いて、当該パラメータを設定したプラント応答モデルにより制御対象プラントを制御すればよい。このように、上記の各実施形態に係るプラント応答推定装置１０は、例えば、モデル予測制御等といった既知の制御手法と組み合わせることにより、制御の追従性と応答の安定化にも寄与することができる。

本発明は、具体的に開示された上記の実施形態に限定されるものではなく、特許請求の範囲の記載から逸脱することなく、種々の変形や変更、既知の技術との組み合わせ等が可能である。

１０ プラント応答推定装置
１１ 入力装置
１２ 表示装置
１３ 外部Ｉ／Ｆ
１３ａ 記録媒体
１４ 通信Ｉ／Ｆ
１５ プロセッサ
１６ メモリ装置
１７ バス
１０１ モデルパラメータ推定部
１０２ ステップ応答計算部
１１１ バッファ部
１１２ 状態ベクトル変換部
１１３ 学習データ選別部
１１４ 逐次推定計算部
１１５ パラメータ変換部

Claims (6)

1. 制御対象プラントの運転データに基づいて、前記制御対象プラントの状態を表す状態ベクトルを作成する状態ベクトル作成部と、
   前記状態ベクトルを用いて、前記状態ベクトルを学習データとするか否かを判定する学習データ選別部と、
   前記状態ベクトルを学習データとすると判定された場合、前記状態ベクトルを学習データとして用いて、前記制御対象プラントの応答モデルを表す関数のパラメータを更新するパラメータ更新部と、
   を有し、
   前記学習データ選別部は、
   前記状態ベクトルの現在値と前回値との差のノルムを計算し、
   前記状態ベクトル作成部によって現在までに作成されたすべての前記状態ベクトルの標準偏差と、０より大きい値を取る所定の係数との積を閾値として計算し、
   前記ノルムが前記閾値以上である場合、前記運転データが定常状態、上限値及び下限値のいずれでもないことにより、前記状態ベクトルを学習データとすると判定し、
   前記ノルムが前記閾値未満である場合、前記運転データが定常状態、上限値又は下限値のいずれかであることにより、前記状態ベクトルを学習データとしないと判定し、
   前記パラメータ更新部は、
   前記学習データとすると判定された前記状態ベクトルを用いて、逐次最小２乗法により共分散行列と前記パラメータの推定値とを、前記学習データとすると判定される毎に逐次的に更新し、
   前記状態ベクトル作成部による前記状態ベクトルの作成と、前記学習データ選別部による前記閾値の計算及び前記学習データとするか否かの判定と、前記パラメータ更新部による前記パラメータの更新とは、前記運転データのサンプリング周期毎に逐次的に実行される、プラント応答推定装置。
2. 制御対象プラントの運転データに基づいて、前記制御対象プラントの状態を表す状態ベクトルを作成する状態ベクトル作成部と、
   前記状態ベクトルを用いて、前記状態ベクトルを学習データとするか否かを判定する学習データ選別部と、
   前記状態ベクトルを学習データとすると判定された場合、前記状態ベクトルを学習データとして用いて、前記制御対象プラントの応答モデルを表す関数のパラメータを更新するパラメータ更新部と、
   を有し、
   前記学習データ選別部は、
   前記状態ベクトルの現在値と前記パラメータとを用いて、前記制御対象プラントの制御量の推定値を計算し、
   前記制御量の推定値と前記制御量の現在値との誤差を表す予測誤差を計算し、
   前記状態ベクトル作成部によって現在までに作成されたすべての前記状態ベクトルの標準偏差と０より大きい値を取る所定の係数との積、又は、現在までのすべての前記制御量の標準偏差と０より大きい値を取る所定の係数との積を、閾値として計算し、
   前記予測誤差が前記閾値以上である場合、前記運転データが定常状態、上限値及び下限値のいずれでもないことにより、前記状態ベクトルを学習データとすると判定し、
   前記予測誤差が前記閾値未満である場合、前記運転データが定常状態、上限値又は下限値のいずれかであることにより、前記状態ベクトルを学習データとしないと判定し、
   前記パラメータ更新部は、
   前記学習データとすると判定された前記状態ベクトルを用いて、逐次最小２乗法により共分散行列と前記パラメータの推定値とを、前記学習データとすると判定される毎に逐次的に更新し、
   前記状態ベクトル作成部による前記状態ベクトルの作成と、前記学習データ選別部による前記閾値の計算及び前記学習データとするか否かの判定と、前記パラメータ更新部による前記パラメータの更新とは、前記運転データのサンプリング周期毎に逐次的に実行される、プラント応答推定装置。
3. 制御対象プラントの運転データに基づいて、前記制御対象プラントの状態を表す状態ベクトルを作成する状態ベクトル作成手順と、
   前記状態ベクトルを用いて、前記状態ベクトルを学習データとするか否かを判定する学習データ選別手順と、
   前記状態ベクトルを学習データとすると判定された場合、前記状態ベクトルを学習データとして用いて、前記制御対象プラントの応答モデルを表す関数のパラメータを更新するパラメータ更新手順と、
   をコンピュータが実行し、
   前記学習データ選別手順は、
   前記状態ベクトルの現在値と前回値との差のノルムを計算し、
   前記状態ベクトル作成手順によって現在までに作成されたすべての前記状態ベクトルの標準偏差と、０より大きい値を取る所定の係数との積を閾値として計算し、
   前記ノルムが前記閾値以上である場合、前記運転データが定常状態、上限値及び下限値のいずれでもないことにより、前記状態ベクトルを学習データとすると判定し、
   前記ノルムが前記閾値未満である場合、前記運転データが定常状態、上限値又は下限値のいずれかであることにより、前記状態ベクトルを学習データとしないと判定し、
   前記パラメータ更新手順は、
   前記学習データとすると判定された前記状態ベクトルを用いて、逐次最小２乗法により共分散行列と前記パラメータの推定値とを、前記学習データとすると判定される毎に逐次的に更新し、
   前記状態ベクトル作成手順による前記状態ベクトルの作成と、前記学習データ選別手順による前記閾値の計算及び前記学習データとするか否かの判定と、前記パラメータ更新手順による前記パラメータの更新とは、前記運転データのサンプリング周期毎に逐次的に実行される、プラント応答推定方法。
4. 制御対象プラントの運転データに基づいて、前記制御対象プラントの状態を表す状態ベクトルを作成する状態ベクトル作成手順と、
   前記状態ベクトルを用いて、前記状態ベクトルを学習データとするか否かを判定する学習データ選別手順と、
   前記状態ベクトルを学習データとすると判定された場合、前記状態ベクトルを学習データとして用いて、前記制御対象プラントの応答モデルを表す関数のパラメータを更新するパラメータ更新手順と、
   をコンピュータが実行し、
   前記学習データ選別手順は、
   前記状態ベクトルの現在値と前記パラメータとを用いて、前記制御対象プラントの制御量の推定値を計算し、
   前記制御量の推定値と前記制御量の現在値との誤差を表す予測誤差を計算し、
   前記状態ベクトル作成手順によって現在までに作成されたすべての前記状態ベクトルの標準偏差と０より大きい値を取る所定の係数との積、又は、現在までのすべての前記制御量の標準偏差と０より大きい値を取る所定の係数との積を、閾値として計算し、
   前記予測誤差が前記閾値以上である場合、前記運転データが定常状態、上限値及び下限値のいずれでもないことにより、前記状態ベクトルを学習データとすると判定し、
   前記予測誤差が前記閾値未満である場合、前記運転データが定常状態、上限値又は下限値のいずれかであることにより、前記状態ベクトルを学習データとしないと判定し、
   前記パラメータ更新手順は、
   前記学習データとすると判定された前記状態ベクトルを用いて、逐次最小２乗法により共分散行列と前記パラメータの推定値とを、前記学習データとすると判定される毎に逐次的に更新し、
   前記状態ベクトル作成手順による前記状態ベクトルの作成と、前記学習データ選別手順による前記閾値の計算及び前記学習データとするか否かの判定と、前記パラメータ更新手順による前記パラメータの更新とは、前記運転データのサンプリング周期毎に逐次的に実行される、プラント応答推定方法。
5. 制御対象プラントの運転データに基づいて、前記制御対象プラントの状態を表す状態ベクトルを作成する状態ベクトル作成手順と、
   前記状態ベクトルを用いて、前記状態ベクトルを学習データとするか否かを判定する学習データ選別手順と、
   前記状態ベクトルを学習データとすると判定された場合、前記状態ベクトルを学習データとして用いて、前記制御対象プラントの応答モデルを表す関数のパラメータを更新するパラメータ更新手順と、
   をコンピュータに実行させ、
   前記学習データ選別手順は、
   前記状態ベクトルの現在値と前回値との差のノルムを計算し、
   前記状態ベクトル作成手順によって現在までに作成されたすべての前記状態ベクトルの標準偏差と、０より大きい値を取る所定の係数との積を閾値として計算し、
   前記ノルムが前記閾値以上である場合、前記運転データが定常状態、上限値及び下限値のいずれでもないことにより、前記状態ベクトルを学習データとすると判定し、
   前記ノルムが前記閾値未満である場合、前記運転データが定常状態、上限値又は下限値のいずれかであることにより、前記状態ベクトルを学習データとしないと判定し、
   前記パラメータ更新手順は、
   前記学習データとすると判定された前記状態ベクトルを用いて、逐次最小２乗法により共分散行列と前記パラメータの推定値とを、前記学習データとすると判定される毎に逐次的に更新し、
   前記状態ベクトル作成手順による前記状態ベクトルの作成と、前記学習データ選別手順による前記閾値の計算及び前記学習データとするか否かの判定と、前記パラメータ更新手順による前記パラメータの更新とは、前記運転データのサンプリング周期毎に逐次的に実行される、プログラム。
6. 制御対象プラントの運転データに基づいて、前記制御対象プラントの状態を表す状態ベクトルを作成する状態ベクトル作成手順と、
   前記状態ベクトルを用いて、前記状態ベクトルを学習データとするか否かを判定する学習データ選別手順と、
   前記状態ベクトルを学習データとすると判定された場合、前記状態ベクトルを学習データとして用いて、前記制御対象プラントの応答モデルを表す関数のパラメータを更新するパラメータ更新手順と、
   をコンピュータに実行させ、
   前記学習データ選別手順は、
   前記状態ベクトルの現在値と前記パラメータとを用いて、前記制御対象プラントの制御量の推定値を計算し、
   前記制御量の推定値と前記制御量の現在値との誤差を表す予測誤差を計算し、
   前記状態ベクトル作成手順によって現在までに作成されたすべての前記状態ベクトルの標準偏差と０より大きい値を取る所定の係数との積、又は、現在までのすべての前記制御量の標準偏差と０より大きい値を取る所定の係数との積を、閾値として計算し、
   前記予測誤差が前記閾値以上である場合、前記運転データが定常状態、上限値及び下限値のいずれでもないことにより、前記状態ベクトルを学習データとすると判定し、
   前記予測誤差が前記閾値未満である場合、前記運転データが定常状態、上限値又は下限値のいずれかであることにより、前記状態ベクトルを学習データとしないと判定し、
   前記パラメータ更新手順は、
   前記学習データとすると判定された前記状態ベクトルを用いて、逐次最小２乗法により共分散行列と前記パラメータの推定値とを、前記学習データとすると判定される毎に逐次的に更新し、
   前記状態ベクトル作成手順による前記状態ベクトルの作成と、前記学習データ選別手順による前記閾値の計算及び前記学習データとするか否かの判定と、前記パラメータ更新手順による前記パラメータの更新とは、前記運転データのサンプリング周期毎に逐次的に実行される、プログラム。

Images