知財求人 - 知財ポータルサイト「IP Force」
特開2024-44583プラント応答推定装置、プラント応答推定方法、及びプログラム
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2024044583
(43)【公開日】2024-04-02
(54)【発明の名称】プラント応答推定装置、プラント応答推定方法、及びプログラム
(51)【国際特許分類】
G05B 13/04 20060101AFI20240326BHJP
【ＦＩ】
G05B13/04
【審査請求】有
【請求項の数】6
【出願形態】ＯＬ
(21)【出願番号】P 2022150202
(22)【出願日】2022-09-21
(11)【特許番号】
(45)【特許公報発行日】2023-03-01
(71)【出願人】
【識別番号】000005234
【氏名又は名称】富士電機株式会社
(74)【代理人】
【識別番号】100107766
【弁理士】
【氏名又は名称】伊東 忠重
(74)【代理人】
【識別番号】100070150
【弁理士】
【氏名又は名称】伊東 忠彦
(72)【発明者】
【氏名】丹下 吉雄
【テーマコード（参考）】
5H004
【Ｆターム（参考）】
5H004GA08
5H004GB02
5H004GB03
5H004GB04
5H004KB01
5H004KC28
5H004LA12
(57)【要約】
【課題】操業中のプラントの運転データから応答モデルのパラメータを安定的に推定することができる技術を提供すること。
【解決手段】本開示の一態様によるプラント応答推定装置は、制御対象プラントの運転データに基づいて、前記制御対象プラントのプラント応答モデルを表す関数のパラメータであるモデルパラメータの推定値を逐次的に計算するように構成されている推定計算部と、前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている安定化計算部と、前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算するように構成されている応答計算部と、を有し、前記推定計算部は、前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算部によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算するように構成されている。
【選択図】図２
【解決手段】本開示の一態様によるプラント応答推定装置は、制御対象プラントの運転データに基づいて、前記制御対象プラントのプラント応答モデルを表す関数のパラメータであるモデルパラメータの推定値を逐次的に計算するように構成されている推定計算部と、前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている安定化計算部と、前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算するように構成されている応答計算部と、を有し、前記推定計算部は、前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算部によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算するように構成されている。
【選択図】図２
【特許請求の範囲】
【請求項1】
制御対象プラントの運転データに基づいて、前記制御対象プラントのプラント応答モデルを表す関数のパラメータであるモデルパラメータの推定値を逐次的に計算するように構成されている推定計算部と、
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている安定化計算部と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算するように構成されている応答計算部と、
を有し、
前記推定計算部は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算部によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算するように構成されている、プラント応答推定装置。
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている安定化計算部と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算するように構成されている応答計算部と、
を有し、
前記推定計算部は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算部によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算するように構成されている、プラント応答推定装置。
【請求項2】
前記安定化計算部は、
前記モデルパラメータの推定値を表すベクトルの各要素から所定の統計量を計算し、前記統計量により前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている、請求項１に記載のプラント応答推定装置。
前記モデルパラメータの推定値を表すベクトルの各要素から所定の統計量を計算し、前記統計量により前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている、請求項１に記載のプラント応答推定装置。
【請求項3】
前記プラント応答モデルは多項式モデルであり、
前記安定化計算部は、
前記多項式モデルの次元毎に、前記モデルパラメータの推定値を表すベクトルの各要素のうち、前記次元方向の要素から前記統計量を計算し、
前記次元毎に、前記次元方向の要素を前記統計量で置換した第１のパラメータを計算し、
前記第１のパラメータと前記モデルパラメータとを所定の比率で混合した第２のパラメータを計算し、
前記第２のパラメータを用いて前記制御対象プラントの応答を予測したときの第１の予測誤差が、前記モデルパラメータを用いて前記制御対象プラントの応答を予測したときの第２の予測誤差よりも所定以上改善している場合は前記第２のパラメータを前記安定化モデルパラメータとし、前記第１の予測誤差が前記第２の予測誤差よりも所定以上改善していない場合は前記モデルパラメータを前記安定化モデルパラメータとする、ように構成されている、請求項２に記載のプラント応答推定装置。
前記安定化計算部は、
前記多項式モデルの次元毎に、前記モデルパラメータの推定値を表すベクトルの各要素のうち、前記次元方向の要素から前記統計量を計算し、
前記次元毎に、前記次元方向の要素を前記統計量で置換した第１のパラメータを計算し、
前記第１のパラメータと前記モデルパラメータとを所定の比率で混合した第２のパラメータを計算し、
前記第２のパラメータを用いて前記制御対象プラントの応答を予測したときの第１の予測誤差が、前記モデルパラメータを用いて前記制御対象プラントの応答を予測したときの第２の予測誤差よりも所定以上改善している場合は前記第２のパラメータを前記安定化モデルパラメータとし、前記第１の予測誤差が前記第２の予測誤差よりも所定以上改善していない場合は前記モデルパラメータを前記安定化モデルパラメータとする、ように構成されている、請求項２に記載のプラント応答推定装置。
【請求項4】
前記多項式モデルは、ＡＲＭＡモデル又はＡＲＭＡＸモデルであり、
前記モデルパラメータは、前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの係数を要素とするベクトルで表される、請求項３に記載のプラント応答推定装置。
前記モデルパラメータは、前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの係数を要素とするベクトルで表される、請求項３に記載のプラント応答推定装置。
【請求項5】
前記安定化計算部は、
前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの自己回帰係数の平均値を表す第１の重心と、前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの移動平均係数の平均値を表す第２の重心とを前記統計量として計算し、
前記自己回帰係数を前記第１の重心、前記移動平均係数を前記第２の重心でそれぞれ置換することで、前記第２のパラメータを計算する、ように構成されている、請求項４に記載のプラント応答推定装置。
前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの自己回帰係数の平均値を表す第１の重心と、前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの移動平均係数の平均値を表す第２の重心とを前記統計量として計算し、
前記自己回帰係数を前記第１の重心、前記移動平均係数を前記第２の重心でそれぞれ置換することで、前記第２のパラメータを計算する、ように構成されている、請求項４に記載のプラント応答推定装置。
【請求項6】
前記ＡＲＭＡモデル又はＡＲＭＡＸモデルには、前記制御対象プラントの定常状態を表すオフセット項が含まれる、請求項４又は５に記載のプラント応答推定装置。
【請求項7】
制御対象プラントの運転データに基づいて、前記制御対象プラントのプラント応答モデルを表す関数のパラメータであるモデルパラメータの推定値を逐次的に計算する推定計算手順と、
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算する安定化計算手順と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算する応答計算手順と、
をコンピュータが実行し、
前記推定計算手順は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算手順によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算する、プラント応答推定方法。
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算する安定化計算手順と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算する応答計算手順と、
をコンピュータが実行し、
前記推定計算手順は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算手順によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算する、プラント応答推定方法。
【請求項8】
制御対象プラントの運転データに基づいて、前記制御対象プラントのプラント応答モデルを表す関数のパラメータであるモデルパラメータの推定値を逐次的に計算する推定計算手順と、
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算する安定化計算手順と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算する応答計算手順と、
をコンピュータに実行させ、
前記推定計算手順は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算手順によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算する、プログラム。
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算する安定化計算手順と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算する応答計算手順と、
をコンピュータに実行させ、
前記推定計算手順は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算手順によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算する、プログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本開示は、プラント応答推定装置、プラント応答推定方法、及びプログラムに関する。
【背景技術】
【0002】
温調制御装置やＰＬＣ（Programmable Logic Controller）、ＤＣＳ（Distributed Control System）等の制御装置、パーソナルコンピュータや組み込み制御機器上で実装される制御装置等が産業上広く利用されている。
【0003】
また、制御対象の制御量を目標値に追従させることを目的とする制御方式として、ＰＩＤ（Proportional-Integral-Differential）制御、モデル予測制御、内部モデル制御、ＬＱＧ（Linear-Quadratic-Gaussian）制御、Ｈ２制御、Ｈ∞制御等の各種の制御方式が知られている。
【0004】
モデル予測制御は、制御対象の状態空間モデルや将来の時間応答モデルを用いた最適化計算を逐次的に行うことで望ましい応答を得る方式であり、産業界で広く用いられている（例えば、特許文献１、非特許文献１）。このようなモデル予測制御等で用いられる線形予測モデルに関してそのモデルパラメータを推定（同定）する手法として、逐次最小２乗法（ＲＬＳ法：Recursive Least Squares法）やカルマンフィルタ法等が知られている（例えば、非特許文献２、非特許文献３）。
【0005】
なお、制御対象プラントの定常状態を表すオフセット項が含まれるモデルパラメータを推定する技術も知られている（例えば、特許文献３）。
【先行技術文献】
【特許文献】
【0006】
【特許文献1】特開２０２０－２１４１１号公報
【特許文献2】特許第７０４７９６６号公報
【非特許文献】
【0007】
【非特許文献1】ヤン Ｍ． マチエヨフスキー，「モデル予測制御 制約のもとでの最適制御」，東京電機大学出版局，２００５
【非特許文献2】足立 修一，「ＭＡＴＬＡＢによる制御のためのシステム同定」，東京電機大学出版局，１９９６
【非特許文献3】足立 修一，「ＭＡＴＬＡＢによる制御のための上級システム同定」，東京電機大学出版局，２００４
【発明の概要】
【発明が解決しようとする課題】
【0008】
逐次最小２乗法やカルマンフィルタ法は、操業中（運転中）のプラントのデータから逐次的にモデルパラメータを推定できる点で大変有益であるが、モデルパラメータの推定値が逐次的に更新されるため以下の２つの課題がある。
【0009】
１つ目の課題は、真値に収束するまでの途中の段階ではモデルパラメータの推定値の信頼性が担保されないという点である。２つ目の課題は、真のモデルパラメータの次数が事前に未知であるため、推定対象のモデルパラメータが過剰な次数を持つ場合には多重共線性が生じ得るという点である。
【0010】
上記の２つの課題により、例えば、本来安定であるはずのプラントに対して不安定極を持つようなモデルパラメータが途中で計算されたり、観測ノイズの影響によってモデルパラメータが大きく変動したり、多重共線性によって複数の局所最適解が生じたりする等といった不安定性が生じる恐れがある。
【0011】
オンラインでモデルパラメータの更新と制御対象の制御とを実現する場合、モデルパラメータの不安定性が制御に重要な影響を与えかねないため、収束途中におけるモデルパラメータの安定性を確保することは重要である。
【0012】
本開示は、上記の点に鑑みてなされたもので、操業中のプラントの運転データから応答モデルのパラメータを安定的に推定することができる技術を提供する。
【課題を解決するための手段】
【0013】
本開示の一態様によるプラント応答推定装置は、制御対象プラントの運転データに基づいて、前記制御対象プラントのプラント応答モデルを表す関数のパラメータであるモデルパラメータの推定値を逐次的に計算するように構成されている推定計算部と、前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている安定化計算部と、前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算するように構成されている応答計算部と、を有し、前記推定計算部は、前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算部によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算するように構成されている。
【発明の効果】
【0014】
操業中のプラントの運転データから応答モデルのパラメータを安定的に推定することができる技術が提供される。
【図面の簡単な説明】
【0015】
【図1】第一の実施形態に係るプラント応答推定装置のハードウェア構成の一例を示す図である。
【図2】第一の実施形態に係るプラント応答推定装置の機能構成の一例を示す図である。
【図3】バッファ部の動作の一例を説明するための図である。
【図4】モデルパラメータ推定処理の一例を説明するためのフローチャートである。
【図5】安定化モデルパラメータ計算処理の一例を説明するためのフローチャートである。
【図6】ステップ応答計算部の動作の一例を説明するための図である。
【図7】ステップ応答計算処理の一例を説明するためのフローチャートである。
【図8】第二の実施形態に係るプラント応答推定装置の機能構成の一例を示す図である。
【図9】パラメータ変換部の動作の一例を説明するための図である。
【図10】一実施例におけるプラントのステップ応答を説明するための図である。
【図11】一実施例で対象とするデータ１の一例を説明するための図である。
【図12】一実施例におけるモデルパラメータの推定結果（安定化モデルパラメータ計算部なし）の一例を説明するための図（その１）である。
【図13】一実施例におけるモデル予測制御（安定化モデルパラメータ計算部なし）の制御結果の一例を説明するための図（その１）である。
【図14】一実施例におけるモデルパラメータの推定結果（安定化モデルパラメータ計算部あり）の一例を説明するための図（その１）である。
【図15】一実施例におけるモデル予測制御（安定化モデルパラメータ計算部あり）の制御結果の一例を説明するための図（その１）である。
【図16】一実施例で対象とするデータ２の一例を説明するための図である。
【図17】一実施例におけるモデルパラメータの推定結果（安定化モデルパラメータ計算部なし）の一例を説明するための図（その２）である。
【図18】一実施例におけるモデル予測制御（安定化モデルパラメータ計算部なし）の制御結果の一例を説明するための図（その２）である。
【図19】一実施例におけるモデルパラメータの推定結果（安定化モデルパラメータ計算部あり）の一例を説明するための図（その２）である。
【図20】一実施例におけるモデル予測制御（安定化モデルパラメータ計算部あり）の制御結果の一例を説明するための図（その２）である。
【発明を実施するための形態】
【0016】
以下、本発明の一実施形態について説明する。以下の各実施形態では、操業中のプラントの運転データから応答モデルのパラメータ（モデルパラメータ）を安定的に推定すると共に、そのパラメータによりプラント応答を推定することができるプラント応答推定装置１０について説明する。なお、操業中とはプラントが正常に稼働し、通常の運転を行っている状態のことであり、例えば、オンライン中、運転中、運用中等と呼ばれてもよい。また、プラントとは１以上の機械、機器、装置等で構成される産業設備のことであり、プラント応答モデルを用いたモデル予測制御によって制御される制御対象である。プラントの具体例としては、例えば、石油化学プラント、食品プラント、鉄鋼プラント、発電プラント等といったものが挙げられるが、これらは一例であって、これらに限られるものではない。
【0017】
［第一の実施形態］
以下、第一の実施形態について説明する。
以下、第一の実施形態について説明する。
【0018】
＜プラント応答推定装置１０のハードウェア構成例＞
本実施形態に係るプラント応答推定装置１０のハードウェア構成例を図１に示す。図１に示すように、本実施形態に係るプラント応答推定装置１０は、入力装置１１と、表示装置１２と、外部Ｉ／Ｆ１３と、通信Ｉ／Ｆ１４と、プロセッサ１５と、メモリ装置１６とを有する。これらの各ハードウェアは、それぞれがバス１７を介して通信可能に接続される。
本実施形態に係るプラント応答推定装置１０のハードウェア構成例を図１に示す。図１に示すように、本実施形態に係るプラント応答推定装置１０は、入力装置１１と、表示装置１２と、外部Ｉ／Ｆ１３と、通信Ｉ／Ｆ１４と、プロセッサ１５と、メモリ装置１６とを有する。これらの各ハードウェアは、それぞれがバス１７を介して通信可能に接続される。
【0019】
入力装置１１は、例えば、キーボード、マウス、タッチパネル、各種物理ボタン等である。表示装置１２は、例えば、ディスプレイ、表示パネル等である。なお、プラント応答推定装置１０は、例えば、入力装置１１及び表示装置１２のうちの少なくとも一方を有していなくてもよい。
【0020】
外部Ｉ／Ｆ１３は、記録媒体１３ａ等の外部装置とのインタフェースである。記録媒体１３ａとしては、例えば、ＣＤ（Compact Disc）、ＤＶＤ（Digital Versatile Disk）、ＳＤメモリカード（Secure Digital memory card）、ＵＳＢ（Universal Serial Bus）メモリカード等が挙げられる。
【0021】
通信Ｉ／Ｆ１４は、プラント応答推定装置１０を通信ネットワークに接続するためのインタフェースである。プロセッサ１５は、例えば、ＣＰＵ（Central Processing Unit）等の各種演算装置である。メモリ装置１６は、例えば、ＨＤＤ（Hard Disk Drive）、ＳＳＤ（Solid State Drive）、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）、フラッシュメモリ等の各種記憶装置である。
【0022】
なお、図１に示すハードウェア構成は一例であって、プラント応答推定装置１０はこれに限られるものではなく、他のハードウェア構成であってもよい。例えば、プラント応答推定装置１０は、複数のプロセッサ１５や複数のメモリ装置１６を有していてもよいし、図示したハードウェアのうちの一部を有していなくてもよいし、図示したハードウェア以外の種々のハードウェアを有していてもよい。
【0023】
＜プラント応答推定装置１０の機能構成例＞
本実施形態に係るプラント応答推定装置１０の機能構成例を図２に示す。図２に示すように、本実施形態に係るプラント応答推定装置１０は、モデルパラメータ推定部１０１と、ステップ応答計算部１０２とを有する。これら各部は、例えば、プラント応答推定装置１０にインストールされた１以上のプログラムが、プロセッサ１５等に実行させる処理により実現される。
本実施形態に係るプラント応答推定装置１０の機能構成例を図２に示す。図２に示すように、本実施形態に係るプラント応答推定装置１０は、モデルパラメータ推定部１０１と、ステップ応答計算部１０２とを有する。これら各部は、例えば、プラント応答推定装置１０にインストールされた１以上のプログラムが、プロセッサ１５等に実行させる処理により実現される。
【0024】
モデルパラメータ推定部１０１は、サンプリング周期Δ毎に、制御対象のプラント（又はその運転状態を計測するセンサ等の機器）から運転データの観測値（制御量ｙ、操作量ｕ及び外乱ｖ）を受信し、そのサンプリング周期Δに応じて、プラント応答モデルのモデルパラメータを推定する。なお、運転データは、例えば、計測データや観測データ等と呼ばれてもよい。
【0025】
ここで、時刻をｔとすれば制御量ｙ、操作量ｕ及び外乱ｖはそれぞれｙ（ｔ）、ｕ（ｔ）及びｖ（ｔ）と表され、各サンプリング時刻ｔｋ（ｋはサンプリング時刻を表すインデックス）に関してｔｋ＋１－ｔｋ＝Δが成り立つ。なお、ｋは０以上の整数である。以下では、特に断らない限り、各サンプリング時刻ｔｋをそのインデックスｋと同一視し、ｙ（ｋ）、ｕ（ｋ）及びｖ（ｋ）とも表すことにする。
【0026】
また、以下では、プラント応答モデルは、パラメータθを持つプラント応答関数Ｓθ（ｔ）で表されるものとして、パラメータθをモデルパラメータということにする。プラント応答モデルとしては様々なプラント応答関数Ｓθで表現されるモデルを採用することが可能であるが、以下では、主に、ＡＲＭＡモデル（Autoregressive Moving Average Model）といった多項式モデルを想定する。プラント応答モデルがＡＲＭＡモデル等といった多項式モデルである場合、モデルパラメータθはその多項式の係数を要素として持つベクトルとなる。なお、プラント応答モデルとしてＡＲＭＡモデル以外の様々なモデルが採用可能であることは言うまでもない。
【0027】
ステップ応答計算部１０２は、モデルパラメータ推定部１０１によって推定されたモデルパラメータ（以下、モデルパラメータ推定値ともいう。）θを用いて、与えられた時刻ｔにおけるプラントのステップ応答の推定値（以下、ステップ応答推定値ともいう。）Ｓθ（ｔ）を計算する。なお、ステップ応答とは、操作量としてステップ信号がプラントに印加されたときの制御量のことである。
【0028】
ここで、モデルパラメータ推定部１０１には、バッファ部１１１と、状態ベクトル変換部１１２と、逐次推定計算部１１３と、安定化モデルパラメータ計算部１１４とが含まれる。
【0029】
バッファ部１１１は、或る所定の期間における制御量ｙ（ｋ）、操作量ｕ（ｋ）及び外乱ｖ（ｋ）をメモリ装置１６に蓄積（バッファ）する。状態ベクトル変換部１１２は、メモリ装置１６にバッファされている制御量ｙ（ｋ）、操作量ｕ（ｋ）及び外乱ｖ（ｋ）の再サンプリングを行って、状態ベクトルと呼ぶベクトルを作成する。逐次推定計算部１１３は、状態ベクトルを用いて、ＲＬＳ法により、プラント応答モデルのモデルパラメータを推定する。安定化モデルパラメータ計算部１１４は、モデルパラメータから安定化モデルパラメータを計算する。安定化モデルパラメータとは、モデルパラメータを或る統計量（例えば、自己回帰係数の平均値、移動平均係数の平均値等）により安定化したベクトルのことである。以下では、安定化モデルパラメータをθＳと表し、またサンプリング時刻ｔｋにおける安定化モデルパラメータをθＳ（ｋ）と表す。
【0030】
＜バッファ部１１１の動作例＞
サンプリング時刻ｔｋにおける制御量バッファをＹ（ｋ）、操作量バッファをＵ（ｋ）、外乱バッファをＶ（ｋ）として、これらの各バッファは以下のベクトルで表されるものとする。
サンプリング時刻ｔｋにおける制御量バッファをＹ（ｋ）、操作量バッファをＵ（ｋ）、外乱バッファをＶ（ｋ）として、これらの各バッファは以下のベクトルで表されるものとする。
【0031】
【数1】
すなわち、制御量バッファＹ（ｋ）にはｋ－Ｂ１からｋまでの制御量ｙ、操作量バッファＵ（ｋ）にはｋ－Ｂ２からｋまでの操作量ｕ、外乱バッファＶ（ｋ）にはｋ－Ｂ３からｋまでの外乱ｖがそれぞれ格納されているものとする。ここで、Ｂ１、Ｂ２及びＢ３はそれぞれ０以上の整数であり、制御量バッファ、操作量バッファ及び外乱バッファの大きさを決めるパラメータである。これらのＢ１、Ｂ２及びＢ３はメモリ装置１６のサイズ等に応じて、適宜、その値が決定される。
【0032】
このとき、バッファ部１１１は、制御量ｙ（ｋ）、操作量ｕ（ｋ）及び外乱ｖ（ｋ）を受信すると、図３に示すように、制御量ｙ（ｋ）と制御量バッファＹ（ｋ－１）から制御量バッファＹ（ｋ）、操作量ｕ（ｋ）と操作量バッファＵ（ｋ－１）から操作量バッファＵ（ｋ）、外乱ｖ（ｋ）と外乱バッファＶ（ｋ－１）から外乱バッファＶ（ｋ）にそれぞれ更新する。
【0033】
具体的には、バッファ部１１１は、制御量バッファＹ（ｋ－１）に格納されているｙ（ｋ－Ｂ１－１）を削除した上で、新たに観測された制御量ｙ（ｋ）を格納することで、ｙ（ｋ－Ｂ１）からｙ（ｋ）までの制御量が格納された制御量バッファＹ（ｋ）に更新する。同様に、バッファ部１１１は、操作量バッファＵ（ｋ－１）に格納されている操作量ｕ（ｋ－Ｂ２－１）を削除した上で、新たに観測された操作量ｕ（ｋ）を格納することで、ｕ（ｋ－Ｂ２）からｕ（ｋ）までの操作量が格納された操作量バッファＵ（ｋ）に更新する。同様に、バッファ部１１１は、外乱バッファＶ（ｋ－１）に格納されている外乱ｖ（ｋ－Ｂ３－１）を削除した上で、新たに観測された外乱ｖ（ｋ）を格納することで、ｖ（ｋ－Ｂ３）からｖ（ｋ）までの外乱が格納された外乱バッファＶ（ｋ）に更新する。
【0034】
＜状態ベクトル変換部１１２の動作例＞
再サンプリング周期をＤとする。このとき、状態ベクトル変換部１１２は、制御量バッファＹ（ｋ）、操作量バッファＵ（ｋ）及び外乱バッファＶ（ｋ）が更新されると、これらの各バッファＹ（ｋ）、Ｕ（ｋ）及びＶ（ｋ）から再サンプリング周期Ｄで再サンプリングを行って、以下の再サンプリング制御量ベクトルＹＤ（ｋ）、再サンプリング操作量ベクトルＵＤ（ｋ）及び再サンプリング外乱ベクトルＶＤ（ｋ）をそれぞれ作成する。
再サンプリング周期をＤとする。このとき、状態ベクトル変換部１１２は、制御量バッファＹ（ｋ）、操作量バッファＵ（ｋ）及び外乱バッファＶ（ｋ）が更新されると、これらの各バッファＹ（ｋ）、Ｕ（ｋ）及びＶ（ｋ）から再サンプリング周期Ｄで再サンプリングを行って、以下の再サンプリング制御量ベクトルＹＤ（ｋ）、再サンプリング操作量ベクトルＵＤ（ｋ）及び再サンプリング外乱ベクトルＶＤ（ｋ）をそれぞれ作成する。
【0035】
【数2】
【0036】
なお、一般に、Ｎ、Ｍ及びＬの値は大きい方が多様な表現が可能で、プラント応答の高精度な予測が期待できるが、モデルパラメータθの推定のために多くの計算資源やメモリ量が必要となる。
【0037】
そして、状態ベクトル変換部１１２は、再サンプリング制御量ベクトルＹＤ（ｋ）、再サンプリング操作量ベクトルＵＤ（ｋ）及び再サンプリング外乱ベクトルＶＤ（ｋ）から以下により状態ベクトルξ（ｋ）を作成する。
【0038】
【数3】
【0039】
＜逐次推定計算部１１３の動作例＞
逐次推定計算部１１３は、状態ベクトルξ（ｋ）が作成されると、この状態ベクトルξ（ｋ）を用いて、モデルパラメータθの値を推定する。すなわち、逐次推定計算部１１３は、サンプリング時刻ｔｋ毎に、逐次的にモデルパラメータθの値を推定する。
逐次推定計算部１１３は、状態ベクトルξ（ｋ）が作成されると、この状態ベクトルξ（ｋ）を用いて、モデルパラメータθの値を推定する。すなわち、逐次推定計算部１１３は、サンプリング時刻ｔｋ毎に、逐次的にモデルパラメータθの値を推定する。
【0040】
或るｋに関してモデルパラメータθの値を推定する処理について、図４を参照しながら説明する。なお、以下では、サンプリング時刻ｔｋにおけるモデルパラメータθ又はその推定値をθ（ｋ）とも表す。
【0041】
ステップＳ１０１：まず、逐次推定計算部１１３は、モデルパラメータθと共分散行列Ｐとを初期化するか否かを判定する。ここで、初期化すると判定される場合としては、例えば、初回計算時（つまり、ｋ＝０のとき）、ユーザ等により初期化指示が行われたとき等が挙げられる。
【0042】
モデルパラメータθと共分散行列Ｐとを初期化すると判定した場合（ステップＳ１０１でＹＥＳ）、逐次推定計算部１１３は、ステップＳ１０２に進む。一方で、モデルパラメータθと共分散行列Ｐとを初期化すると判定しなかった場合（ステップＳ１０１でＮＯ）、逐次推定計算部１１３は、ステップＳ１０３に進む。
【0043】
ステップＳ１０２：逐次推定計算部１１３は、サンプリング時刻ｔｋのインデックスｋをｋ＝０に初期化すると共に、θ（０）＝θ０及びＰ（０）＝Ｉと初期化する。ここで、θ０は予め設定されたモデルパラメータの初期値、Ｉは予め設定された任意の行列（例えば、単位行列等）である。
【0044】
ステップＳ１０３：逐次推定計算部１１３は、安定化モデルパラメータ計算部１１４を呼び出して安定化モデルパラメータθＳ（ｋ－１）を計算する。なお、安定化モデルパラメータθＳ（ｋ－１）の計算方法については後述する。
【0045】
ステップＳ１０４：逐次推定計算部１１３は、安定化モデルパラメータθＳ（ｋ－１）をモデルパラメータ推定値θ（ｋ－１）に設定する。すなわち、逐次推定計算部１１３は、θ（ｋ－１）←θＳ（ｋ－１）とする。これは、モデルパラメータ推定値θ（ｋ－１）を安定化モデルパラメータθＳ（ｋ－１）に変換している、ということもできる。
【0046】
ステップＳ１０５：逐次推定計算部１１３は、モデルパラメータ推定値θ（ｋ－１）と状態ベクトルξ（ｋ）と制御量ｙ（ｋ）とを用いて、予測誤差ε（ｋ）を計算する。予測誤差ε（ｋ）は、例えば、ε（ｋ）＝ｙ（ｋ）－ξ（ｋ）Τθ（ｋ－１）により計算される。なお、Τは転置を表す。
【0047】
ステップＳ１０６：逐次推定計算部１１３は、共分散行列Ｐ（ｋ－１）を以下により更新して共分散行列Ｐ（ｋ）を得る。
【0048】
【数4】
【0049】
ステップＳ１０７：そして、逐次推定計算部１１３は、モデルパラメータ推定値θ（ｋ－１）を以下により更新してモデルパラメータ推定値θ（ｋ）を得る。
【0050】
【数5】
【0051】
＜安定化モデルパラメータ計算部１１４の動作例＞
安定化モデルパラメータ計算部１１４は、逐次推定計算部１１３から呼び出されると、安定化モデルパラメータθＳ（ｋ－１）を計算する。すなわち、安定化モデルパラメータ計算部１１４は、サンプリング時刻ｔｋで逐次推定計算部１１３から呼び出される毎に、モデルパラメータ推定値θ（ｋ－１）から安定化モデルパラメータθＳ（ｋ－１）を計算する。
安定化モデルパラメータ計算部１１４は、逐次推定計算部１１３から呼び出されると、安定化モデルパラメータθＳ（ｋ－１）を計算する。すなわち、安定化モデルパラメータ計算部１１４は、サンプリング時刻ｔｋで逐次推定計算部１１３から呼び出される毎に、モデルパラメータ推定値θ（ｋ－１）から安定化モデルパラメータθＳ（ｋ－１）を計算する。
【0052】
或るサンプリング時刻ｔｋで逐次推定計算部１１３から呼び出されたときに、安定化モデルパラメータ計算部１１４が安定化モデルパラメータを計算する処理について、図５を参照しながら説明する。ただし、以下では、一例として、プラント応答モデルとしては以下のＡＲＭＡモデルを想定する。
【0053】
ｙ（ｋ）＝ａ１ｙ（ｋ－１）＋・・・＋ａＮｙ（ｋ－Ｎ）＋ｂ０ｕ（ｋ）＋・・・＋ｂＭｕ（ｋ－Ｍ）
また、このとき、モデルパラメータ推定値θ（ｋ－１）は以下で表される。
また、このとき、モデルパラメータ推定値θ（ｋ－１）は以下で表される。
【0054】
θ（ｋ－１）＝［ａ１，・・・，ａＮ，ｂ０，・・・，ｂＭ］Τ
ここで、Ｎは制御量ｙに関する次数、Ｍは操作量ｕに関する次数である。すなわち、上記のプラント応答モデルは（Ｎ，Ｍ）次のＡＲＭＡモデルである。なお、ａ１，・・・，ａＮは自己回帰係数、Ｎは自己回帰係数の次数とも呼ばれる。また、ｂ０，・・・，ｂＭは移動平均係数、Ｍは移動平均係数の次数とも呼ばれる。
ここで、Ｎは制御量ｙに関する次数、Ｍは操作量ｕに関する次数である。すなわち、上記のプラント応答モデルは（Ｎ，Ｍ）次のＡＲＭＡモデルである。なお、ａ１，・・・，ａＮは自己回帰係数、Ｎは自己回帰係数の次数とも呼ばれる。また、ｂ０，・・・，ｂＭは移動平均係数、Ｍは移動平均係数の次数とも呼ばれる。
【0055】
ステップＳ２０１：まず、安定化モデルパラメータ計算部１１４は、重心化パラメータθｃを計算する。安定化モデルパラメータ計算部１１４は、自己回帰係数ａ１，・・・，ａＮの重心ａｃと移動平均係数ｂ０，・・・，ｂＭの重心ｂｃとをそれぞれ計算した上で、モデルパラメータ推定値θ（ｋ－１）の各自己回帰係数を重心ａｃ、各移動平均係数を重心ｂｃにそれぞれ置き換えたものを重心化パラメータθｃとする。すなわち、重心化パラメータθｃは以下で表される。
【0056】
θｃ＝［ａｃ，・・・，ａｃ，ｂｃ，・・・，ｂｃ］Τ
ここで、ａｃ＝（ａ１＋・・・＋ａＮ）／Ｎ、ｂｃ＝（ｂ０，・・・，ｂＭ）／（Ｍ＋１）である。
ここで、ａｃ＝（ａ１＋・・・＋ａＮ）／Ｎ、ｂｃ＝（ｂ０，・・・，ｂＭ）／（Ｍ＋１）である。
【0057】
ステップＳ２０２：次に、安定化モデルパラメータ計算部１１４は、安定化モデルパラメータの候補となるパラメータ（以下、候補パラメータともいう。）計算する。安定化モデルパラメータ計算部１１４は、以下により候補パラメータθｒを計算する。
【0058】
θｒ＝（１－β）×θ（ｋ－１）＋β×θｃ
ここで、βは０＜β≦１を取るブレンド係数であり、予め設定された値である。ブレンド係数βが１に近いほど候補パラメータθｒは重心化パラメータθｃに近づき、０に近いほど候補パラメータθｒはモデルパラメータ推定値θ（ｋ－１）に近づく。
ここで、βは０＜β≦１を取るブレンド係数であり、予め設定された値である。ブレンド係数βが１に近いほど候補パラメータθｒは重心化パラメータθｃに近づき、０に近いほど候補パラメータθｒはモデルパラメータ推定値θ（ｋ－１）に近づく。
【0059】
ステップＳ２０３：次に、安定化モデルパラメータ計算部１１４は、２種類の予測誤差を計算する。すなわち、安定化モデルパラメータ計算部１１４は、モデルパラメータ推定値θ（ｋ－１）を用いたときの予測誤差ｅ０と、候補パラメータθｒを用いたときの予測誤差ｅｒとをそれぞれ以下により計算する。
【0060】
ｅ０＝ｙ０－ξ（ｋ）Τθ（ｋ－１）
ｅｒ＝ｙ０－ξ（ｋ）Τθｒ
ここで、ｙ０は制御量の現在値、すなわちｙ０＝ｙ（ｋ）である。また、ξ（ｋ）は状態ベクトルである。すなわち、ｅ０は、モデルパラメータ推定値θ（ｋ－１）を用いたときの制御量予測値と、制御量現在値との誤差を表している。同様に、ｅｒは、候補パラメータθｒを用いたときの制御量予測値と、制御量現在値との誤差を表している。
ｅｒ＝ｙ０－ξ（ｋ）Τθｒ
ここで、ｙ０は制御量の現在値、すなわちｙ０＝ｙ（ｋ）である。また、ξ（ｋ）は状態ベクトルである。すなわち、ｅ０は、モデルパラメータ推定値θ（ｋ－１）を用いたときの制御量予測値と、制御量現在値との誤差を表している。同様に、ｅｒは、候補パラメータθｒを用いたときの制御量予測値と、制御量現在値との誤差を表している。
【0061】
ステップＳ２０４：次に、安定化モデルパラメータ計算部１１４は、予測誤差ｅ０の絶対値と予測誤差ｅｒの絶対値とを比較する。安定化モデルパラメータ計算部１１４は、例えば、予測誤差ｅ０の絶対値と予測誤差ｅｒの絶対値とを比較し、以下を満たすか否かを判定する。
【0062】
｜ｅｒ｜＜α｜ｅ０｜
ここで、αは０＜α≦１を取る比較係数であり、予め設定された値である。比較係数αは、候補パラメータθｒを用いたときの予測誤差ｅｒをどの程度厳しく評価するかを制御しており、０に近いほど、予測誤差ｅ０に対して予測誤差ｅｒの改善が顕著であることが要求される。
ここで、αは０＜α≦１を取る比較係数であり、予め設定された値である。比較係数αは、候補パラメータθｒを用いたときの予測誤差ｅｒをどの程度厳しく評価するかを制御しており、０に近いほど、予測誤差ｅ０に対して予測誤差ｅｒの改善が顕著であることが要求される。
【0063】
上記の比較の結果、｜ｅｒ｜＜α｜ｅ０｜を満たす場合、安定化モデルパラメータ計算部１１４は、ステップＳ２０５に進む。一方で、｜ｅｒ｜＜α｜ｅ０｜を満たさない場合（つまり、｜ｅｒ｜≧α｜ｅ０｜を満たす場合）、安定化モデルパラメータ計算部１１４は、ステップＳ２０６に進む。
【0064】
ステップＳ２０５：｜ｅｒ｜＜α｜ｅ０｜を満たす場合、安定化モデルパラメータ計算部１１４は、候補パラメータθｒを安定化モデルパラメータθＳ（ｋ－１）とする。すなわち、安定化モデルパラメータ計算部１１４は、θＳ（ｋ－１）←θｒとする。
【0065】
ステップＳ２０６：｜ｅｒ｜＜α｜ｅ０｜を満たさない場合、安定化モデルパラメータ計算部１１４は、モデルパラメータ推定値θ（ｋ－１）を安定化モデルパラメータθＳ（ｋ－１）とする。すなわち、安定化モデルパラメータ計算部１１４は、θＳ（ｋ－１）←θ（ｋ－１）とする。
【0066】
なお、上記では、一例として、外乱を考慮していないが、外乱が存在する場合であっても同様に適用できることはいうまでもない。外乱を考慮する場合、例えば、ＡＲＭＡモデルには外乱に関する移動平均モデルが追加され、その結果、モデルパラメータには外乱に関する自己回帰係数が追加される。
【0067】
また、上記では、一例として、自己回帰係数ａ１，・・・，ａＮの平均値をａｃ、移動平均係数ｂ０，・・・，ｂＭの平均値をｂｃとしたが、これに限られるものではない。ａｃとして、自己回帰係数ａ１，・・・，ａＮから計算可能な任意の統計量を用いることが可能である。同様に、ｂｃとして、移動平均係数ｂ０，・・・，ｂＭから計算可能な任意の統計量を用いることが可能である。より一般には、次元毎に、その次元方向の係数から計算可能な任意の統計量を用いることが可能である。このような統計量としては、平均の他、例えば、加重平均、中央値、最頻値等といったものが挙げられる。
【0068】
＜ステップ応答計算部１０２の動作例＞
以下、ステップ応答の推定に用いられるモデルパラメータ推定値（つまり、プラント応答関数Ｓθに設定されるモデルパラメータ推定値）θ＝θ（ｋ）のことを「モデルパラメータ設定値θ」ともいう。また、以下では、簡単のため、ステップ信号として単位ステップ信号を想定する。
以下、ステップ応答の推定に用いられるモデルパラメータ推定値（つまり、プラント応答関数Ｓθに設定されるモデルパラメータ推定値）θ＝θ（ｋ）のことを「モデルパラメータ設定値θ」ともいう。また、以下では、簡単のため、ステップ信号として単位ステップ信号を想定する。
【0069】
ステップ応答計算部１０２は、図６に示すように、モデルパラメータ設定値θと時刻ｔとが与えられると、初期時刻０から時間ｔ経過後の時刻ｔにおける単位ステップ応答Ｓθ（ｔ）を計算する。なお、単位ステップ応答とは、操作量ｕとして単位ステップ信号を印加した場合における応答（つまり、制御対象プラントのプラント応答モデルの出力）のことである。
【0070】
或る時刻ｔと或るモデルパラメータ設定値θとが与えられたときに、ステップ応答Ｓθ（ｔ）を計算する処理について、図７を参照しながら説明する。ただし、以下では、一例として、プラント応答モデルとしては以下のＡＲＭＡモデルを想定する。
【0071】
ｙ（ｋ'）＝ａ１ｙ（ｋ'－１）＋・・・＋ａＮｙ（ｋ'－Ｎ）＋ｂ０ｕ（ｋ'）＋・・・＋ｂＭｕ（ｋ'－Ｍ）
また、このとき、インデックスｋ'における状態ベクトルφ（ｋ'）を以下で表すものとする。
また、このとき、インデックスｋ'における状態ベクトルφ（ｋ'）を以下で表すものとする。
【0072】
φ（ｋ'）＝［ｙ（ｋ'－１），ｙ（ｋ'－２），・・・，ｙ（ｋ'－Ｎ），ｕ（ｋ'），ｕ（ｋ'－１），・・・，ｕ（ｋ'－Ｍ）］Τ
なお、インデックスｋ'は、サンプリング時刻ｔｋのインデックスｋと同じ値を取り得る変数であるが、本処理の中でのみ利用され、インデックスｋとは独立に値が更新されることに留意されたい。
なお、インデックスｋ'は、サンプリング時刻ｔｋのインデックスｋと同じ値を取り得る変数であるが、本処理の中でのみ利用され、インデックスｋとは独立に値が更新されることに留意されたい。
【0073】
ステップＳ３０１：ステップ応答計算部１０２は、本処理の中でのみ利用する時刻を表すインデックスをτとして、τ＝０、ｋ'＝０と初期化すると共に、状態ベクトルφ（０）を以下のように初期化する。
【0074】
φ（０）＝［０，０，・・・，０，１，０，・・・，０］Τ
すなわち、ｕ（０）のみ１、それ以外の要素は０と状態ベクトルφ（０）を初期化する。これは、単位ステップ応答を模擬する際の初期値を表している。
すなわち、ｕ（０）のみ１、それ以外の要素は０と状態ベクトルφ（０）を初期化する。これは、単位ステップ応答を模擬する際の初期値を表している。
【0075】
ステップＳ３０２：次に、ステップ応答計算部１０２は、ｙ（ｋ'）＝φ（ｋ'）Τθにより制御量予測値ｙ（ｋ'）を計算する。
【0076】
ステップＳ３０３：次に、ステップ応答計算部１０２は、制御量予測値ｙ（ｋ'）を用いて、状態ベクトルφ（ｋ'）を、次のインデックスｋ'＋１における状態ベクトルφ（ｋ'＋１）に更新する。このとき、ステップ応答計算部１０２は、状態ベクトルφ（ｋ'＋１）のｙ（ｋ'）には上記のステップＳ３０２で計算した制御量予測値ｙ（ｋ'）を設定し、ｙ（ｋ'－Ｎ＋１）～ｙ（ｋ'－１）には状態ベクトルφ（ｋ'）と同じ値を設定する。また、状態ベクトルφ（ｋ'＋１）のｕ（ｋ'＋１）には１を設定し、ｕ（ｋ'－Ｍ＋１）～ｕ（ｋ'）には状態ベクトルφ（ｋ'）と同じ値を設定する。
【0077】
ステップＳ３０４：次に、ステップ応答計算部１０２は、時刻τをτ＋Δに更新すると共に、インデックスｋ'をｋ'＋１に更新する。
【0078】
ステップＳ３０５：次に、ステップ応答計算部１０２は、τ≧ｔであるか否かを判定する。そして、τ≧ｔであると判定されなかった場合（ステップＳ３０５でＮＯ）、ステップ応答計算部１０２は、ステップＳ３０２に戻る。これにより、τ≧ｔとなるまで、ステップＳ３０２～ステップＳ３０４が繰り返し実行される。
【0079】
一方で、τ≧ｔであると判定された場合（ステップＳ３０５でＹＥＳ）、ステップ応答計算部１０２は、処理を終了する。これにより、最終的に計算された制御量予測値ｙ（ｋ'）が単位ステップ応答Ｓθ（ｔ）として得られる（つまり、Ｓθ（ｔ）＝ｙ（ｋ）が、プラント応答モデルの単位ステップ応答として得られる。）。
【0080】
［第二の実施形態］
以下、第二の実施形態について説明する。第二の実施形態では、制御対象のプラントの定常状態を表すオフセット項が含まれるモデルパラメータを対象とする。なお、プラントの定常状態を表すオフセット項が含まれるモデルパラメータとそれに関連する技術については、例えば、上記の特許文献３等を参照されたい。
以下、第二の実施形態について説明する。第二の実施形態では、制御対象のプラントの定常状態を表すオフセット項が含まれるモデルパラメータを対象とする。なお、プラントの定常状態を表すオフセット項が含まれるモデルパラメータとそれに関連する技術については、例えば、上記の特許文献３等を参照されたい。
【0081】
具体的には、第二の実施形態では、プラント応答モデルとして以下のオフセット項が含まれるＡＲＭＡＸモデル（Autoregressive Moving Average Model with Exogenous Variables）を想定する。
【0082】
ｙ（ｋ）＝ａ１ｙ（ｋ－１）＋・・・＋ａＮｙ（ｋ－Ｎ）＋ｂ０ｕ（ｋ）＋・・・＋ｂＭｕ（ｋ－Ｍ）＋ｃ０
ここで、ｃ０がオフセット項である。
ここで、ｃ０がオフセット項である。
【0083】
このとき、モデルパラメータ（及びその推定値）θは、以下で表されるモデルパラメータに拡大することができる。
【0084】
θｅ＝［ａ１，・・・，ａＮ，ｂ０，・・・，ｂＭ，ｃ０］Τ
このモデルパラメータ（及びその推定値）θｅを「拡大モデルパラメータ」とも呼ぶ。
このモデルパラメータ（及びその推定値）θｅを「拡大モデルパラメータ」とも呼ぶ。
【0085】
なお、第二の実施形態では、主に、第一の実施形態と相違点について説明し、第一の実施形態と同様としてよい箇所についてはその説明を省略する。
【0086】
＜プラント応答推定装置１０の機能構成例＞
本実施形態に係るプラント応答推定装置１０の機能構成例を図８に示す。図８に示すように、本実施形態に係るプラント応答推定装置１０は、パラメータ変換部１１５がモデルパラメータ推定部１０１に含まれている点が第一の実施形態と異なる。
本実施形態に係るプラント応答推定装置１０の機能構成例を図８に示す。図８に示すように、本実施形態に係るプラント応答推定装置１０は、パラメータ変換部１１５がモデルパラメータ推定部１０１に含まれている点が第一の実施形態と異なる。
【0087】
パラメータ変換部１１５は、拡大モデルパラメータθｅをモデルパラメータθに変換する。
【0088】
＜状態ベクトル変換部１１２の動作例＞
状態ベクトル変換部１１２は、状態ベクトルを拡大した拡大状態ベクトルと呼ぶベクトルを作成する点が第一の実施形態と異なる。その他の点は第一の実施形態と同様としてよい。
状態ベクトル変換部１１２は、状態ベクトルを拡大した拡大状態ベクトルと呼ぶベクトルを作成する点が第一の実施形態と異なる。その他の点は第一の実施形態と同様としてよい。
【0089】
拡大状態ベクトルとは、状態ベクトルξ（ｋ）をオフセット項に対応する要素を加えて拡大したものであり、以下で表される。
【0090】
【数6】
【0091】
＜逐次推定計算部１１３の動作例＞
逐次推定計算部１１３は、モデルパラメータθの代わりに拡大モデルパラメータθｅを用いる点、状態ベクトルξ（ｋ）の代わりに拡大状態ベクトルξｅ（ｋ）を用いる点が第一の実施形態と異なる。その他の点は第一の実施形態と同様としてよい。
逐次推定計算部１１３は、モデルパラメータθの代わりに拡大モデルパラメータθｅを用いる点、状態ベクトルξ（ｋ）の代わりに拡大状態ベクトルξｅ（ｋ）を用いる点が第一の実施形態と異なる。その他の点は第一の実施形態と同様としてよい。
【0092】
＜安定化モデルパラメータ計算部１１４の動作例＞
安定化モデルパラメータ計算部１１４は、モデルパラメータθの代わりに拡大モデルパラメータθｅを用いる点、状態ベクトルξ（ｋ）の代わりに拡大状態ベクトルξｅ（ｋ）を用いる点が第一の実施形態と異なる。その他の点は第一の実施形態と同様としてよい。
安定化モデルパラメータ計算部１１４は、モデルパラメータθの代わりに拡大モデルパラメータθｅを用いる点、状態ベクトルξ（ｋ）の代わりに拡大状態ベクトルξｅ（ｋ）を用いる点が第一の実施形態と異なる。その他の点は第一の実施形態と同様としてよい。
【0093】
なお、上記のステップＳ２０１で重心化パラメータθｃを計算する際に、オフセット項は置き換えの対象外であることに留意されたい。すなわち、モデルパラメータ推定値θ（ｋ－１）がθｅ（ｋ－１）＝［ａ１，・・・，ａＮ，ｂ０，・・・，ｂＭ，ｃ０］Τである場合、重心化パラメータθｃはθｃ＝［ａｃ，・・・，ａｃ，ｂｃ，・・・，ｂｃ，ｃ０］Τとなる。
【0094】
＜パラメータ変換部１１５の動作例＞
以下、一例として、拡大モデルパラメータ推定値θｅ（ｋ）は以下で表されるものとする。
以下、一例として、拡大モデルパラメータ推定値θｅ（ｋ）は以下で表されるものとする。
【0095】
θｅ（ｋ）＝［θＹ（ｋ），θＵ（ｋ），θＶ（ｋ），θＣ（ｋ）］Τ
ここで、θＹ（ｋ）はサンプリング時刻ｔｋにおけるＡＲＭＡＸモデルの制御量ｙに関する項の係数を要素とするＮ次元ベクトル、θＵ（ｋ）は操作量ｕに関する項の係数を要素とするＭ次元ベクトル、θＶ（ｋ）は外乱ｖに関する項の係数を要素とするＬ次元ベクトル、θＣ（ｋ）はオフセット項を表すスカラー値である。
ここで、θＹ（ｋ）はサンプリング時刻ｔｋにおけるＡＲＭＡＸモデルの制御量ｙに関する項の係数を要素とするＮ次元ベクトル、θＵ（ｋ）は操作量ｕに関する項の係数を要素とするＭ次元ベクトル、θＶ（ｋ）は外乱ｖに関する項の係数を要素とするＬ次元ベクトル、θＣ（ｋ）はオフセット項を表すスカラー値である。
【0096】
このとき、パラメータ変換部１１５は、拡大モデルパラメータ推定値θｅ（ｋ）が得られると、図９に示すように、θＣ（ｋ）を除いたベクトルに変換することで、モデルパラメータ推定値θ（ｋ）を得る。これにより、以下のモデルパラメータ推定値が得られる。
【0097】
θ（ｋ）＝［θＹ（ｋ），θＵ（ｋ），θＶ（ｋ）］Τ
このように、制御対象プラントのステップ応答を推定する際には、拡大モデルパラメータ推定値θｅ（ｋ）からθＣ（ｋ）を除くことで、モデルパラメータ推定値θ（ｋ）が得られる。この理由は、上記の特許文献１に記載されている通り、制御量ｙの変化量に対してオフセット項は影響しないためである。
このように、制御対象プラントのステップ応答を推定する際には、拡大モデルパラメータ推定値θｅ（ｋ）からθＣ（ｋ）を除くことで、モデルパラメータ推定値θ（ｋ）が得られる。この理由は、上記の特許文献１に記載されている通り、制御量ｙの変化量に対してオフセット項は影響しないためである。
【0098】
［実施例］
以下、一実施例について説明する。本実施例では、第二の実施形態に係るプラント応答推定装置１０によりモデルパラメータとステップ応答を推定した。
以下、一実施例について説明する。本実施例では、第二の実施形態に係るプラント応答推定装置１０によりモデルパラメータとステップ応答を推定した。
【0099】
制御対象プラントのプラント応答モデルとしては、以下のＡＲＭＡＸモデルを採用した。
【0100】
ｙ（ｋ）＝ａ１ｙ（ｋ－１）＋ａ２ｙ（ｋ－２）＋ｂ０ｕ（ｋ）＋ｂ１ｕ（ｋ－１）＋ｂ２ｕ（ｋ－２）＋ｃ０
すなわち、制御量ｙに関する２次元の要素と、操作量ｕに関する３次元の要素と、オフセット項とを有するＡＲＭＡＸモデルを採用した。
すなわち、制御量ｙに関する２次元の要素と、操作量ｕに関する３次元の要素と、オフセット項とを有するＡＲＭＡＸモデルを採用した。
【0101】
このとき、拡大状態ベクトルξｅ（ｋ）は以下で表される。
【0102】
ξｅ（ｋ）＝［ｙ（ｋ－１），ｙ（ｋ－２），ｕ（ｋ），ｕ（ｋ－１），ｕ（ｋ－２），１］Τ
また、拡大モデルパラメータθｅは以下で表される。
また、拡大モデルパラメータθｅは以下で表される。
【0103】
θｅ（ｋ）＝［ａ１（ｋ），ａ２（ｋ），ｂ０（ｋ），ｂ１（ｋ），ｂ２（ｋ），ｃ０（ｋ）］Τ
制御量の予測値（以下、これをｙｅｓｔとも表す。）は以下で表される。
制御量の予測値（以下、これをｙｅｓｔとも表す。）は以下で表される。
【0104】
ｙｅｓｔ（ｋ）＝ξｅ（ｋ）Τθｅ（ｋ－１）
ここで、本実施例における制御対象プラントのステップ応答は、図１０に示すようなものであるとする。なお、ｒは目標値、ｙは制御量、ｕは操作量、ｖは外乱を表す。また、本実施例では、再サンプリング周期Ｄ＝２、忘却係数λ＝０．９９９、比較係数α＝０．９とした。
ここで、本実施例における制御対象プラントのステップ応答は、図１０に示すようなものであるとする。なお、ｒは目標値、ｙは制御量、ｕは操作量、ｖは外乱を表す。また、本実施例では、再サンプリング周期Ｄ＝２、忘却係数λ＝０．９９９、比較係数α＝０．９とした。
【0105】
このとき、本実施例では、２つのデータ（後述するデータ１、データ２）をそれぞれ用いて、モデルパラメータとステップ応答を推定した。
【0106】
・データ１を用いた場合
制御対象プラントをＰＩＤ制御で制御したデータ１を図１１に示す。このとき、安定化モデルパラメータ計算部１１４を利用しないでモデルパラメータを推定した結果とモデル予測制御の制御結果とをそれぞれ図１２及び図１３に示す。また、第二の実施形態に係るプラント応答推定装置１０によりモデルパラメータを推定した結果とモデル予測制御の制御結果とをそれぞれ図１４及び図１５に示す。図１２及び図１４では、モデルパラメータのａ１、ａ２、ｂ０、ｂ１、ｂ２を時系列で示している。図１２（安定化モデルパラメータ計算部１１４なし）と図１４（安定化モデルパラメータ計算部１１４あり）とを比較すると、自己回帰係数ａ１、ａ２は安定化モデルパラメータ計算部１１４ありの方が重心に近く、移動平均係数ｂ０、ｂ１、ｂ２は安定化モデルパラメータ計算部１１４ありの方が重心に近いことがわかる。一方で、図１３（安定化モデルパラメータ計算部１１４なし）と図１５（安定化モデルパラメータ計算部１１４あり）とを比較すると、大きな違いはなく、安定化モデルパラメータ計算部１１４によって制御性能の劣化は生じていないことがわかる。
制御対象プラントをＰＩＤ制御で制御したデータ１を図１１に示す。このとき、安定化モデルパラメータ計算部１１４を利用しないでモデルパラメータを推定した結果とモデル予測制御の制御結果とをそれぞれ図１２及び図１３に示す。また、第二の実施形態に係るプラント応答推定装置１０によりモデルパラメータを推定した結果とモデル予測制御の制御結果とをそれぞれ図１４及び図１５に示す。図１２及び図１４では、モデルパラメータのａ１、ａ２、ｂ０、ｂ１、ｂ２を時系列で示している。図１２（安定化モデルパラメータ計算部１１４なし）と図１４（安定化モデルパラメータ計算部１１４あり）とを比較すると、自己回帰係数ａ１、ａ２は安定化モデルパラメータ計算部１１４ありの方が重心に近く、移動平均係数ｂ０、ｂ１、ｂ２は安定化モデルパラメータ計算部１１４ありの方が重心に近いことがわかる。一方で、図１３（安定化モデルパラメータ計算部１１４なし）と図１５（安定化モデルパラメータ計算部１１４あり）とを比較すると、大きな違いはなく、安定化モデルパラメータ計算部１１４によって制御性能の劣化は生じていないことがわかる。
【0107】
・データ２を用いた場合
制御対象プラントをＰＩＤ制御で制御したデータ２を図１６に示す。データ２では、制御量ｙに観測ノイズが加わっており、その結果、操作量ｕにもノイズが加わっている。なお、ノイズの影響以外はデータ１とデータ２とで条件は同一である。このとき、安定化モデルパラメータ計算部１１４を利用しないでモデルパラメータを推定した結果とモデル予測制御の制御結果とをそれぞれ図１７及び図１８に示す。また、第二の実施形態に係るプラント応答推定装置１０によりモデルパラメータを推定した結果とモデル予測制御の制御結果とをそれぞれ図１９及び図２０に示す。図１７及び図１９では、モデルパラメータのａ１、ａ２、ｂ０、ｂ１、ｂ２を時系列で示している。図１７（安定化モデルパラメータ計算部１１４なし）と図１９（安定化モデルパラメータ計算部１１４あり）とを比較すると、自己回帰係数ａ１、ａ２は安定化モデルパラメータ計算部１１４ありの方が重心に近く、移動平均係数ｂ０、ｂ１、ｂ２は安定化モデルパラメータ計算部１１４ありの方が重心に近いことがわかる。また、安定化モデルパラメータ計算部１１４なしの場合（図１７）では、自己回帰係数ａ１、ａ２が徐々に互いに離れた値となっていくが、安定化モデルパラメータ計算部１１４ありの場合（図１９）では、その動きが抑制できていることがわかる。すなわち、ノイズのあるデータ２ではノイズによって自己回帰係数及び移動平均係数が重心から離れた値で推定されるが、安定化モデルパラメータ計算部１１４によって重心から離れる動きが抑制できている。
制御対象プラントをＰＩＤ制御で制御したデータ２を図１６に示す。データ２では、制御量ｙに観測ノイズが加わっており、その結果、操作量ｕにもノイズが加わっている。なお、ノイズの影響以外はデータ１とデータ２とで条件は同一である。このとき、安定化モデルパラメータ計算部１１４を利用しないでモデルパラメータを推定した結果とモデル予測制御の制御結果とをそれぞれ図１７及び図１８に示す。また、第二の実施形態に係るプラント応答推定装置１０によりモデルパラメータを推定した結果とモデル予測制御の制御結果とをそれぞれ図１９及び図２０に示す。図１７及び図１９では、モデルパラメータのａ１、ａ２、ｂ０、ｂ１、ｂ２を時系列で示している。図１７（安定化モデルパラメータ計算部１１４なし）と図１９（安定化モデルパラメータ計算部１１４あり）とを比較すると、自己回帰係数ａ１、ａ２は安定化モデルパラメータ計算部１１４ありの方が重心に近く、移動平均係数ｂ０、ｂ１、ｂ２は安定化モデルパラメータ計算部１１４ありの方が重心に近いことがわかる。また、安定化モデルパラメータ計算部１１４なしの場合（図１７）では、自己回帰係数ａ１、ａ２が徐々に互いに離れた値となっていくが、安定化モデルパラメータ計算部１１４ありの場合（図１９）では、その動きが抑制できていることがわかる。すなわち、ノイズのあるデータ２ではノイズによって自己回帰係数及び移動平均係数が重心から離れた値で推定されるが、安定化モデルパラメータ計算部１１４によって重心から離れる動きが抑制できている。
【0108】
一方で、図１８（安定化モデルパラメータ計算部１１４なし）と図２０（安定化モデルパラメータ計算部１１４あり）とを比較すると、図１８では操作量ｕが比較的鋭く変化している。これに対して、図２０では、上記のような鋭い変化は少なく、ノイズのないデータ１での制御結果（図１５）の応答波形と比較的良く似た結果となっている。したがって、安定化モデルパラメータ計算部１１４により、ノイズの影響を抑制しており、モデル予測制御の安定化に寄与しているといえる。
【0109】
［まとめ］
以上のように、第一及び第二の実施形態に係るプラント応答推定装置１０は、操業中のプラントの運転データからプラント応答モデルのモデルパラメータを安定的に推定することができる。すなわち、重心化パラメータによってモデルパラメータの推定には重心方向に推定されるような力が働くため、例えば、図１７のように、ａ１、ａ２が徐々に互いに離れた値となっていく、というような動き（つまり、モデルパラメータが発散するような動き）を抑制することが可能となり、その結果、モデルパラメータを安定的に推定することが可能となる。これにより、真値に収束するまでの途中の段階でもモデルパラメータ推定値の信頼性が担保され、上記の１つ目の課題が解決される。このため、第一及び第二の実施形態に係るプラント応答推定装置１０によれば、推定したモデルパラメータを制御（例えば、モデル予測制御）に容易に用いることが可能となる。
以上のように、第一及び第二の実施形態に係るプラント応答推定装置１０は、操業中のプラントの運転データからプラント応答モデルのモデルパラメータを安定的に推定することができる。すなわち、重心化パラメータによってモデルパラメータの推定には重心方向に推定されるような力が働くため、例えば、図１７のように、ａ１、ａ２が徐々に互いに離れた値となっていく、というような動き（つまり、モデルパラメータが発散するような動き）を抑制することが可能となり、その結果、モデルパラメータを安定的に推定することが可能となる。これにより、真値に収束するまでの途中の段階でもモデルパラメータ推定値の信頼性が担保され、上記の１つ目の課題が解決される。このため、第一及び第二の実施形態に係るプラント応答推定装置１０によれば、推定したモデルパラメータを制御（例えば、モデル予測制御）に容易に用いることが可能となる。
【0110】
また、第一及び第二の実施形態に係るプラント応答推定装置１０では、モデルパラメータの次数が未知であり、必要以上に高い次数が設定された場合であっても、重心方向に推定されるような力が働くため、重心に近いモデルパラメータが推定されやすくなり、多重共線性等の高次数に由来する問題を回避することができるようになる。これにより、上記の２つ目の課題が解決される。
【0111】
更に、第一及び第二の実施形態に係るプラント応答推定装置１０では、安定化モデルパラメータ計算部１１４にて現状のモデルパラメータを用いた場合の予測誤差と候補パラメータを用いた場合の予測誤差とを比較し、予測精度が改善する場合にのみ候補パラメータを安定化モデルパラメータとして採用する。このため、予測精度の改善も期待できる。
【0112】
加えて、第一及び第二の実施形態に係るプラント応答推定装置１０では、観測値に対するノイズに対してもモデルパラメータ推定値が安定するため、ノイズの影響を抑制した推定が可能となる。
【0113】
本発明は、具体的に開示された上記の実施形態に限定されるものではなく、特許請求の範囲の記載から逸脱することなく、種々の変形や変更、既知の技術との組み合わせ等が可能である。
【符号の説明】
【0114】
１０ プラント応答推定装置
１１ 入力装置
１２ 表示装置
１３ 外部Ｉ／Ｆ
１３ａ 記録媒体
１４ 通信Ｉ／Ｆ
１５ プロセッサ
１６ メモリ装置
１７ バス
１０１ モデルパラメータ推定部
１０２ ステップ応答計算部
１１１ バッファ部
１１２ 状態ベクトル変換部
１１３ 逐次推定計算部
１１４ 安定化モデルパラメータ計算部
１１５ パラメータ変換部
１１ 入力装置
１２ 表示装置
１３ 外部Ｉ／Ｆ
１３ａ 記録媒体
１４ 通信Ｉ／Ｆ
１５ プロセッサ
１６ メモリ装置
１７ バス
１０１ モデルパラメータ推定部
１０２ ステップ応答計算部
１１１ バッファ部
１１２ 状態ベクトル変換部
１１３ 逐次推定計算部
１１４ 安定化モデルパラメータ計算部
１１５ パラメータ変換部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【手続補正書】
【提出日】2022-12-20
【手続補正1】
【補正対象書類名】特許請求の範囲
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【特許請求の範囲】
【請求項1】
制御対象プラントの運転データに基づいて、前記制御対象プラントのプラント応答モデルを表す関数のパラメータであるモデルパラメータの推定値を逐次的に計算するように構成されている推定計算部と、
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている安定化計算部と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算するように構成されている応答計算部と、
を有し、
前記推定計算部は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算部によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算するように構成されており、
前記プラント応答モデルは多項式モデルであり、
前記安定化計算部は、
前記多項式モデルの次元毎に、前記モデルパラメータの推定値を表すベクトルの各要素のうち、前記次元方向の要素から所定の統計量を計算し、
前記次元毎に、前記次元方向の要素を前記統計量で置換した第１のパラメータを計算し、
前記第１のパラメータと前記モデルパラメータとを所定の比率で混合した第２のパラメータを計算し、
前記第２のパラメータを用いて前記制御対象プラントの応答を予測したときの第１の予測誤差が、前記モデルパラメータを用いて前記制御対象プラントの応答を予測したときの第２の予測誤差よりも所定以上改善している場合は前記第２のパラメータを前記安定化モデルパラメータとし、前記第１の予測誤差が前記第２の予測誤差よりも所定以上改善していない場合は前記モデルパラメータを前記安定化モデルパラメータとする、ように構成されている、プラント応答推定装置。
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算するように構成されている安定化計算部と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算するように構成されている応答計算部と、
を有し、
前記推定計算部は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算部によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算するように構成されており、
前記プラント応答モデルは多項式モデルであり、
前記安定化計算部は、
前記多項式モデルの次元毎に、前記モデルパラメータの推定値を表すベクトルの各要素のうち、前記次元方向の要素から所定の統計量を計算し、
前記次元毎に、前記次元方向の要素を前記統計量で置換した第１のパラメータを計算し、
前記第１のパラメータと前記モデルパラメータとを所定の比率で混合した第２のパラメータを計算し、
前記第２のパラメータを用いて前記制御対象プラントの応答を予測したときの第１の予測誤差が、前記モデルパラメータを用いて前記制御対象プラントの応答を予測したときの第２の予測誤差よりも所定以上改善している場合は前記第２のパラメータを前記安定化モデルパラメータとし、前記第１の予測誤差が前記第２の予測誤差よりも所定以上改善していない場合は前記モデルパラメータを前記安定化モデルパラメータとする、ように構成されている、プラント応答推定装置。
【請求項2】
前記多項式モデルは、ＡＲＭＡモデル又はＡＲＭＡＸモデルであり、
前記モデルパラメータは、前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの係数を要素とするベクトルで表される、請求項１に記載のプラント応答推定装置。
前記モデルパラメータは、前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの係数を要素とするベクトルで表される、請求項１に記載のプラント応答推定装置。
【請求項3】
前記安定化計算部は、
前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの自己回帰係数の平均値を表す第１の重心と、前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの移動平均係数の平均値を表す第２の重心とを前記統計量として計算し、
前記自己回帰係数を前記第１の重心、前記移動平均係数を前記第２の重心でそれぞれ置換することで、前記第２のパラメータを計算する、ように構成されている、請求項２に記載のプラント応答推定装置。
前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの自己回帰係数の平均値を表す第１の重心と、前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの移動平均係数の平均値を表す第２の重心とを前記統計量として計算し、
前記自己回帰係数を前記第１の重心、前記移動平均係数を前記第２の重心でそれぞれ置換することで、前記第２のパラメータを計算する、ように構成されている、請求項２に記載のプラント応答推定装置。
【請求項4】
前記ＡＲＭＡモデル又はＡＲＭＡＸモデルには、前記制御対象プラントの定常状態を表すオフセット項が含まれる、請求項２又は３に記載のプラント応答推定装置。
【請求項5】
制御対象プラントの運転データに基づいて、前記制御対象プラントのプラント応答モデルを表す関数のパラメータであるモデルパラメータの推定値を逐次的に計算する推定計算手順と、
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算する安定化計算手順と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算する応答計算手順と、
をコンピュータが実行し、
前記推定計算手順は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算手順によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算し、
前記プラント応答モデルは多項式モデルであり、
前記安定化計算手順は、
前記多項式モデルの次元毎に、前記モデルパラメータの推定値を表すベクトルの各要素のうち、前記次元方向の要素から所定の統計量を計算し、
前記次元毎に、前記次元方向の要素を前記統計量で置換した第１のパラメータを計算し、
前記第１のパラメータと前記モデルパラメータとを所定の比率で混合した第２のパラメータを計算し、
前記第２のパラメータを用いて前記制御対象プラントの応答を予測したときの第１の予測誤差が、前記モデルパラメータを用いて前記制御対象プラントの応答を予測したときの第２の予測誤差よりも所定以上改善している場合は前記第２のパラメータを前記安定化モデルパラメータとし、前記第１の予測誤差が前記第２の予測誤差よりも所定以上改善していない場合は前記モデルパラメータを前記安定化モデルパラメータとする、プラント応答推定方法。
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算する安定化計算手順と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算する応答計算手順と、
をコンピュータが実行し、
前記推定計算手順は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算手順によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算し、
前記プラント応答モデルは多項式モデルであり、
前記安定化計算手順は、
前記多項式モデルの次元毎に、前記モデルパラメータの推定値を表すベクトルの各要素のうち、前記次元方向の要素から所定の統計量を計算し、
前記次元毎に、前記次元方向の要素を前記統計量で置換した第１のパラメータを計算し、
前記第１のパラメータと前記モデルパラメータとを所定の比率で混合した第２のパラメータを計算し、
前記第２のパラメータを用いて前記制御対象プラントの応答を予測したときの第１の予測誤差が、前記モデルパラメータを用いて前記制御対象プラントの応答を予測したときの第２の予測誤差よりも所定以上改善している場合は前記第２のパラメータを前記安定化モデルパラメータとし、前記第１の予測誤差が前記第２の予測誤差よりも所定以上改善していない場合は前記モデルパラメータを前記安定化モデルパラメータとする、プラント応答推定方法。
【請求項6】
制御対象プラントの運転データに基づいて、前記制御対象プラントのプラント応答モデルを表す関数のパラメータであるモデルパラメータの推定値を逐次的に計算する推定計算手順と、
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算する安定化計算手順と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算する応答計算手順と、
をコンピュータに実行させ、
前記推定計算手順は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算手順によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算し、
前記プラント応答モデルは多項式モデルであり、
前記安定化計算手順は、
前記多項式モデルの次元毎に、前記モデルパラメータの推定値を表すベクトルの各要素のうち、前記次元方向の要素から所定の統計量を計算し、
前記次元毎に、前記次元方向の要素を前記統計量で置換した第１のパラメータを計算し、
前記第１のパラメータと前記モデルパラメータとを所定の比率で混合した第２のパラメータを計算し、
前記第２のパラメータを用いて前記制御対象プラントの応答を予測したときの第１の予測誤差が、前記モデルパラメータを用いて前記制御対象プラントの応答を予測したときの第２の予測誤差よりも所定以上改善している場合は前記第２のパラメータを前記安定化モデルパラメータとし、前記第１の予測誤差が前記第２の予測誤差よりも所定以上改善していない場合は前記モデルパラメータを前記安定化モデルパラメータとする、プログラム。
前記モデルパラメータの推定値を安定化した安定化モデルパラメータを計算する安定化計算手順と、
前記モデルパラメータの推定値に基づいて、前記プラント応答モデルを表す関数により前記制御対象プラントの応答を計算する応答計算手順と、
をコンピュータに実行させ、
前記推定計算手順は、
前記運転データと、前回計算したモデルパラメータの推定値を前記安定化計算手順によって安定化した安定化モデルパラメータとに基づいて、前記モデルパラメータの推定値を計算し、
前記プラント応答モデルは多項式モデルであり、
前記安定化計算手順は、
前記多項式モデルの次元毎に、前記モデルパラメータの推定値を表すベクトルの各要素のうち、前記次元方向の要素から所定の統計量を計算し、
前記次元毎に、前記次元方向の要素を前記統計量で置換した第１のパラメータを計算し、
前記第１のパラメータと前記モデルパラメータとを所定の比率で混合した第２のパラメータを計算し、
前記第２のパラメータを用いて前記制御対象プラントの応答を予測したときの第１の予測誤差が、前記モデルパラメータを用いて前記制御対象プラントの応答を予測したときの第２の予測誤差よりも所定以上改善している場合は前記第２のパラメータを前記安定化モデルパラメータとし、前記第１の予測誤差が前記第２の予測誤差よりも所定以上改善していない場合は前記モデルパラメータを前記安定化モデルパラメータとする、プログラム。