(57)【特許請求の範囲】
【請求項１】
制御対象に対する操作量を出力し、前記制御対象の制御量を目標値に追従させる制御装置であって、
前記目標値の時系列である目標値時系列と、前記目標値を先読みする時間幅を示す先読み長とが入力されると、前記目標値時系列に含まれる複数の目標値のうち、前記先読み長後の目標値を示す先読み目標値を取得する目標値先読み手段と、
前記先読み目標値と、前記制御対象の現在の制御量との差である先読み目標偏差を算出する先読み目標偏差算出手段と、
前記制御対象のプラント応答モデルを表す関数であって、かつ、モデルパラメータが含まれる関数であるプラント応答関数と、現在に至るまでの過去の前記操作量の変化量とに基づいて、前記先読み長後の制御量の予測値と、前記先読み目標値との差を示す補正目標偏差を算出する補正目標偏差算出手段と、
前記補正目標偏差に基づいて、新たな前記操作量を算出する操作量算出手段と、
前記制御量と、前記操作量とに基づいて、前記プラント応答関数に含まれるモデルパラメータの推定値を算出するモデルパラメータ推定手段と、
を有し、
前記目標値先読み手段による前記先読み目標値の取得と、前記先読み目標偏差算出手段による前記先読み目標偏差の算出と、前記補正目標偏差算出手段による前記補正目標偏差の算出と、前記操作量算出手段による新たな前記操作量の算出と、前記モデルパラメータ推定手段による前記モデルパラメータの推定値の算出とを制御周期毎に逐次的に繰り返す、ことを特徴とする制御装置。
【請求項２】
ＯＮ又はＯＦＦのいずれかの値を取るフラグがＯＮになった時点における前記制御量を初期値として、前記モデルパラメータ推定手段により算出された推定値と、前記操作量とに基づいて、前記時点以降の制御量の推定値を算出する制御量推定手段、
を有することを特徴とする請求項１に記載の制御装置。
【請求項３】
前記制御量推定手段は、
前記モデルパラメータ推定手段により算出された推定値を用いた前記制御量の推定値の算出と、前記プラント応答関数に含まれるモデルパラメータを用いた前記制御量の推定値の算出とを行う、
ことを特徴とする請求項２に記載の制御装置。
【請求項４】
ＯＮ又はＯＦＦのいずれかの値を取る更新トリガに基づいて、前記更新トリガの値がＯＮである場合は、前記モデルパラメータ推定手段により算出された推定値を、前記プラント応答関数に含まれるパラメータに設定するモデルパラメータ更新手段、
を有することを特徴とする請求項１乃至３の何れか一項に記載の制御装置。
【請求項５】
前記モデルパラメータ推定手段は、
前記制御量と、前記操作量と、前記制御対象に対する外乱の計測値とに基づいて、前記プラント応答関数に含まれるモデルパラメータの推定値を算出する、
ことを特徴とする請求項１乃至４の何れか一項に記載の制御装置。
【請求項６】
前記プラント応答関数は、ＡＲＭＡモデル又はＡＲＭＡＸモデルから計算するステップ応答モデルで表され、
前記モデルパラメータは、前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの係数であり、
前記先読み長後の制御量の予測値は、前記先読み長後の時点におけるステップ応答を前記ＡＲＭＡモデル又はＡＲＭＡＸモデルにより計算することで算出される、ことを特徴とする請求項１乃至５の何れか一項に記載の制御装置。
【請求項７】
前記モデルパラメータ推定手段は、
逐次最小２乗法又はカルマンフィルタの手法により、共分散行列の更新と前記モデルパラメータの推定値の算出とを行う、
ことを特徴とする請求項１乃至６の何れか一項に記載の制御装置。
【請求項８】
前記モデルパラメータ推定手段は、
前記共分散行列の更新において、前記共分散行列の非対角項の絶対値を減少させる補正処理を行う、
ことを特徴とする請求項７に記載の制御装置。
【請求項９】
制御対象に対する操作量を出力し、前記制御対象の制御量を目標値に追従させるコンピュータが、
前記目標値の時系列である目標値時系列と、前記目標値を先読みする時間幅を示す先読み長とが入力されると、前記目標値時系列に含まれる複数の目標値のうち、前記先読み長後の目標値を示す先読み目標値を取得する目標値先読み手順と、
前記先読み目標値と、前記制御対象の現在の制御量との差である先読み目標偏差を算出する先読み目標偏差算出手順と、
前記制御対象のプラント応答モデルを表す関数であって、かつ、モデルパラメータが含まれる関数であるプラント応答関数と、現在に至るまでの過去の前記操作量の変化量とに基づいて、前記先読み長後の制御量の予測値と、前記先読み目標値との差を示す補正目標偏差を算出する補正目標偏差算出手順と、
前記補正目標偏差に基づいて、新たな前記操作量を算出する操作量算出手順と、
前記制御量と、前記操作量とに基づいて、前記プラント応答関数に含まれるモデルパラメータの推定値を算出するモデルパラメータ推定手順と、
を実行し、
前記目標値先読み手順による前記先読み目標値の取得と、前記先読み目標偏差算出手順による前記先読み目標偏差の算出と、前記補正目標偏差算出手順による前記補正目標偏差の算出と、前記操作量算出手順による新たな前記操作量の算出と、前記モデルパラメータ推定手順による前記モデルパラメータの推定値の算出とを制御周期毎に逐次的に繰り返す、ことを特徴とする制御方法。
【請求項１０】
制御対象に対する操作量を出力し、前記制御対象の制御量を目標値に追従させるコンピュータを、
前記目標値の時系列である目標値時系列と、前記目標値を先読みする時間幅を示す先読み長とが入力されると、前記目標値時系列に含まれる複数の目標値のうち、前記先読み長後の目標値を示す先読み目標値を取得する目標値先読み手段、
前記先読み目標値と、前記制御対象の現在の制御量との差である先読み目標偏差を算出する先読み目標偏差算出手段、
前記制御対象のプラント応答モデルを表す関数であって、かつ、モデルパラメータが含まれる関数であるプラント応答関数と、現在に至るまでの過去の前記操作量の変化量とに基づいて、前記先読み長後の制御量の予測値と、前記先読み目標値との差を示す補正目標偏差を算出する補正目標偏差算出手段、
前記補正目標偏差に基づいて、新たな前記操作量を算出する操作量算出手段、
前記制御量と、前記操作量とに基づいて、前記プラント応答関数に含まれるモデルパラメータの推定値を算出するモデルパラメータ推定手段、
として機能させ、
前記目標値先読み手段による前記先読み目標値の取得と、前記先読み目標偏差算出手段による前記先読み目標偏差の算出と、前記補正目標偏差算出手段による前記補正目標偏差の算出と、前記操作量算出手段による新たな前記操作量の算出と、前記モデルパラメータ推定手段による前記モデルパラメータの推定値の算出とを制御周期毎に逐次的に繰り返す、プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、制御装置、制御方法及びプログラムに関する。
【背景技術】
【０００２】
温調制御装置やＰＬＣ（Programmable Logic Controller）、ＤＣＳ（Distributed Control System）等の制御装置、パーソナルコンピュータや組み込み制御機器上で実装される制御装置等が産業上広く利用されている。
【０００３】
また、制御対象の制御量を目標値に追従させることを目的とする制御方式として、ＰＩＤ（Proportional-Integral-Differential）制御、モデル予測制御、内部モデル制御、ＬＱＧ（Linear-Quadratic-Gaussian）制御、Ｈ２制御、Ｈ∞制御等の各種の制御方式が知られている。
【０００４】
モデル予測制御は、制御対象の状態空間モデルや将来の時間応答モデルを用いた最適化計算を逐次的に行うことで望ましい応答を得る方式であり、産業界で広く用いられている（例えば、非特許文献１）。例えば、オンラインで数値最適化アルゴリズムを実行する標準的なモデル予測制御の産業応用として、空調システムの制御等への応用が知られている（例えば、特許文献１）。
【０００５】
また、現在に至るまでの過去の操作量の変化に応じた制御量の収束値の予測値と、目標値との差である補正目標偏差に基づいて新たな操作量を決定する制御装置が提案されている（例えば、特許文献２）。目標偏差の現在値と、操作量の変化量との関係を示す制御用論理式が成立する領域をグラフ上に表示する設計支援装置も提案されている（例えば、特許文献３）。
【０００６】
また、線形予測モデルに対してパラメータ同定を行う手法として、逐次最小２乗法、すなわち、ＲＬＳ（Recursive Least Squares）法やカルマンフィルタ法等が知られている（例えば、非特許文献２及び非特許文献３）。
【０００７】
更に、例えば、発電プラントの最適運転において、発電要素の物理モデルと入力状態量から特性パラメータを推定し、得られた特性パラメータから最適負荷配分を決定する制御方法も知られている。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開２０１５−１０８４９９号公報
【特許文献２】国際公開第２０１６／０９２８７２号
【特許文献３】国際公開第２０１５／０６０１４９号
【非特許文献】
【０００９】
【非特許文献１】ヤン Ｍ． マチエヨフスキー，「モデル予測制御  制約のもとでの最適制御」，東京電機大学出版局，２００５
【非特許文献２】足立 修一，「ＭＡＴＬＡＢによる制御のためのシステム同定」，東京電機大学出版局，１９９６
【非特許文献３】足立 修一，「ＭＡＴＬＡＢによる制御のための上級システム同定」，東京電機大学出版局，２００４
【発明の概要】
【発明が解決しようとする課題】
【００１０】
しかしながら、従来のモデル予測制御は、逐次的な最適化計算を繰り返し実行する必要があるため、ＣＰＵ（Central Processing Unit）やメモリ等の計算資源が豊富なＰＣ（パーソナルコンピュータ）上で実行するものが多く、比較的計算資源が乏しいエッジデバイス（例えば、ＰＬＣ（Programmable Logic Controller）等の制御装置）上で実行することが困難であった。
【００１１】
また、モデル予測制御はモデルに基づく制御方式であるため、制御対象のモデルを予め同定する必要があるが、一度同定したモデルであっても、制御対象の経年変化や季節性変動等の要因によって変化していく場合がある。
【００１２】
本発明の一実施形態は、上記の点に鑑みてなされたもので、エッジデバイス上でモデル予測制御を実現しつつ、モデル更新のためのパラメータを推定することを目的とする。
【課題を解決するための手段】
【００１３】
上記目的を達成するため、本発明の一実施形態に係る制御装置は、制御対象に対する操作量を出力し、前記制御対象の制御量を目標値に追従させる制御装置であって、前記目標値の時系列である目標値時系列と、前記目標値を先読みする時間幅を示す先読み長とが入力されると、前記目標値時系列に含まれる複数の目標値のうち、前記先読み長後の目標値を示す先読み目標値を取得する目標値先読み手段と、前記先読み目標値と、前記制御対象の現在の制御量との差である先読み目標偏差を算出する先読み目標偏差算出手段と、前記制御対象のプラント応答を表すプラント応答関数と、現在に至るまでの過去の前記操作量の変化量とに基づいて、前記先読み長後の制御量の予測値と、前記先読み目標値との差を示す補正目標偏差を算出する補正目標偏差算出手段と、前記補正目標偏差に基づいて、新たな前記操作量を算出する操作量算出手段と、前記制御量と、前記操作量とに基づいて、前記プラント応答関数に含まれるモデルパラメータの推定値を算出するモデルパラメータ推定手段と、を有することを特徴とする。
【発明の効果】
【００１４】
エッジデバイス上でモデル予測制御を実現しつつ、モデル更新のためのパラメータを推定することができる。
【図面の簡単な説明】
【００１５】
【図１】第一の実施形態に係る制御装置の構成の一例を示す図である。
【図２】プラント応答関数の動作の一例を説明するための図である。
【図３】プラント応答関数の計算処理の一例を説明するためのフローチャートである。
【図４】モデルパラメータの推定処理の一例を説明するためのフローチャートである。
【図５】目標値先読み部の動作の一例を説明するための図である。
【図６】先読み応答補正部の動作の一例を説明するための図（その１）である。
【図７】先読み応答補正部の動作の一例を説明するための図（その２）である。
【図８】予測時系列更新部が実行する更新処理の一例を説明するためのフローチャートである。
【図９】予測時系列の更新の一例を説明するための図である。
【図１０】操作変化量算出部の動作の一例を説明するための図である。
【図１１】一実施形態に係る制御装置のハードウェア構成の一例を示す図である。
【図１２】第二の実施形態に係る制御装置の構成の一例を示す図である。
【図１３】制御量推定処理の一例を説明するためのフローチャートである。
【図１４】第三の実施形態に係る制御装置の構成の一例を示す図である。
【図１５】第四の実施形態に係る制御装置の構成の一例を示す図である。
【図１６】実施例１におけるステップ応答を説明するための図である。
【図１７】実施例１におけるモデルパラメータ推定時の設定を説明するための図である。
【図１８】実施例１における操作量及び制御量を説明するための図である。
【図１９】実施例１におけるモデルパラメータ推定値及び制御量推定値を説明するための図である。
【図２０】実施例２における操作量及び制御量を説明するための図である。
【図２１】実施例２におけるモデルパラメータ推定値及び制御量推定値を説明するための図である。
【発明を実施するための形態】
【００１６】
以下、本発明の実施の形態について、図面を参照しながら詳細に説明する。以降では、現在から未来への目標値時系列が入力された場合において、制御量を目標値に追従させるための操作量を計算すると共に、モデル更新のためのモデルパラメータを推定する制御装置１０について説明する。なお、以降で説明する各実施形態に係る制御装置１０は、例えば、ＰＬＣ等のエッジデバイスである。
【００１７】
以降で説明する各実施形態に係る制御装置１０は、任意の目標値の時系列データである目標値時系列｛ｒ（ｔ）｝や制御対象プラント２０の状態等を示す制御量ｙ、制御対象プラント２０のプラント応答モデル等に基づいて制御対象プラント２０に対する操作量ｕを算出する共に、モデル更新のためのモデルパラメータを推定する。そして、制御装置１０は、この操作量ｕに応じた制御対象プラント２０の制御量ｙを計測し、目標値時系列｛ｒ（ｔ）｝や制御量ｙ、プラント応答モデル等に基づいて次の操作量ｕを算出すると共に、新たなモデルパラメータを推定する。このように、以降で説明する各実施形態に係る制御装置１０は、制御量ｙを目標値時系列｛ｒ（ｔ）｝に追従させるための操作量ｕの算出と、プラント応答モデルのモデルパラメータの推定とをオンライン実行中（つまり、制御対象プラント２０の制御中）に繰り返し実行する。
【００１８】
以降で説明する各実施形態に係る制御装置１０は、モデル予測制御の高い制御性能を実現しつつも、毎回のオンライン最適化計算を必要としないため、少ない計算量と少ないメモリとで実現することが可能になると共に、プラント応答モデルのモデルパラメータを推定することが可能となる。なお、各実施形態に係る制御装置１０は、例えば、既設の他の制御装置の操作量と制御量とを用いてプラント応答モデルのパラメータを推定することが可能である。このため、例えば、制御装置１０の設置直後は他の制御装置の操作量と制御量とを用いてプラント応答モデルの初期モデルパラメータを推定し、運用開始後は季節変動や制御対象プラント２０の経年劣化等の継時的変化に追従したモデルパラメータを推定する、といったことが可能である。
【００１９】
なお、制御量ｙとしては、例えば、制御対象プラント２０の温度、目標値ｒとしては、例えば、設定温度等が挙げられる。ただし、制御量ｙ及び目標値ｒは、温度及び設定温度に限られず、制御対象プラント２０における任意の制御量及び当該制御量の目標となる目標値を用いることができる。
【００２０】
［第一の実施形態］
＜制御装置１０の構成＞
まず、本実施形態に係る制御装置１０の構成について、図１を参照しながら説明する。図１は、第一の実施形態に係る制御装置１０の構成の一例を示す図である。
【００２１】
図１に示すように、本実施形態に係る制御装置１０は、モデルパラメータ推定部１０１と、目標値先読み部１０２と、計測部１０３と、差分器１０４と、操作量更新部１０５と、タイマ１０６とを有する。これら各機能部は、例えば、制御装置１０にインストールされた１以上のプログラムがＣＰＵ等に実行させる処理により実現される。
【００２２】
目標値先読み部１０２は、所定の制御周期Ｔｃ毎に、目標値時系列｛ｒ（ｔ）｝と、先読み長Ｔｐとが入力されると、現在時刻ｔから先読み長後の時刻ｔ＋Ｔｐにおける目標値ｒ（ｔ＋Ｔｐ）を出力する。先読み長Ｔｐとは、目標値時系列｛ｒ（ｔ）｝のうち、先読みする目標値ｒ（ｔ＋Ｔｐ）を決定するための時間長である。なお、以降では、目標値ｒ（ｔ＋Ｔｐ）を「先読み目標値ｒ（ｔ＋Ｔｐ）」とも表す。
【００２３】
計測部１０３は、制御周期Ｔｃ毎に、制御対象プラント２０の制御量ｙを計測する。そして、計測部１０３は、計測した制御量ｙの最新の値を、制御量現在値ｙ０として出力する。なお、制御対象プラント２０の制御量ｙは、操作量ｕと外乱ｖとに応じて決定される。外乱ｖとしては、例えば、制御量ｙが温度である場合には外気温の低下又は上昇等が挙げられる。
【００２４】
また、計測部１０３は、操作量更新部１０５から出力された操作量ｕを取得し、取得した操作量ｕの最新の値を、操作量現在値ｕ０として出力する。
【００２５】
差分器１０４は、目標値先読み部１０２から出力された先読み目標値ｒ（ｔ＋Ｔｐ）と、制御量現在値ｙ０との差（偏差）を目標偏差ｅ０として出力する。目標偏差ｅ０は、ｅ０＝ｒ（ｔ＋Ｔｐ）−ｙ０（ｔ）で算出される。なお、以降では、目標偏差ｅ０を「先読み目標偏差ｅ０」とも表す。
【００２６】
操作量更新部１０５は、制御周期Ｔｃ毎に、制御対象プラント２０に対する操作量ｕを出力する。操作量更新部１０５には、先読み応答補正部１１１と、操作変化量算出部１１２と、加算器１１３とが含まれる。
【００２７】
先読み応答補正部１１１は、プラント応答関数｛Ｓθ（ｔ）｝と、先読み目標偏差ｅ０（ｔ）と、先読み長Ｔｐと、過去の操作量ｕの変化量ｄｕの時系列データである操作変化量時系列｛ｄｕ（ｔ）｝とに基づいて、先読み目標偏差ｅ０（ｔ）を補正した補正目標偏差ｅ＊（ｔ）を算出する。補正目標偏差ｅ＊（ｔ）の算出方法の詳細については後述する。ここで、プラント応答関数｛Ｓθ（ｔ）｝とはモデルパラメータθを含む関数であり、制御対象プラント２０のプラント応答モデルである。なお、以降では、操作量ｕの変化量ｄｕを「操作変化量ｄｕ」とも表す。
【００２８】
操作変化量算出部１１２は、制御周期Ｔｃ毎に、先読み応答補正部１１１により算出された補正目標偏差ｅ＊（ｔ）に基づいて、操作変化量ｄｕを算出する。操作変化量算出部１１２は、例えば、ｄｕ（ｔ−３Ｔｃ）、ｄｕ（ｔ−２Ｔｃ）、ｄｕ（ｔ−Ｔｃ）という順序で操作変化量ｄｕ（ｔ）を算出し、出力する。なお、操作変化量ｄｕは、制御周期Ｔｃ毎に操作量ｕが変化した量である。
【００２９】
加算器１１３は、計測部１０３から出力された操作量現在値ｕ０と、操作変化量算出部１１２から出力された操作変化量ｄｕとを加算して、新たな操作量ｕを算出する。そして、加算器１１３は、この操作量ｕを制御対象プラント２０に出力する。この操作量ｕは、ｕ（ｔ）＝ｕ０＋ｄｕ（ｔ）＝ｕ（ｔ−Ｔｃ）＋ｄｕ（ｔ）で算出される。
【００３０】
タイマ１０６は、制御周期Ｔｃ毎に、目標値先読み部１０２と計測部１０３とを動作させる。なお、目標値先読み部１０２と計測部１０３とが制御周期Ｔｃ毎に動作することで、モデルパラメータ推定部１０１及び操作量更新部１０５も制御周期Ｔｃ毎に動作する。
【００３１】
モデルパラメータ推定部１０１は、制御量現在値ｙ０と、操作量現在値ｕ０と、現在のプラント応答関数｛Ｓθ（ｔ）｝に設定されているモデルパラメータθとが入力されると、プラント応答関数｛Ｓθ（ｔ）｝のモデルパラメータの推定値θｅｓｔを算出し、出力する。なお、以降では、現在のプラント応答関数｛Ｓθ（ｔ）｝に設定されているモデルパラメータθを「モデルパラメータ現在値θ」とも表し、プラント応答関数｛Ｓθ（ｔ）｝のモデルパラメータの推定値θｅｓｔを「モデルパラメータ推定値θｅｓｔ」とも表す。
【００３２】
＜プラント応答関数｛Ｓθ（ｔ）｝の動作＞
次に、プラント応答関数｛Ｓθ（ｔ）｝の動作について、図２を参照しながら説明する。図２は、プラント応答関数｛Ｓθ（ｔ）｝の動作の一例を説明するための図である。
【００３３】
図２に示すように、プラント応答関数｛Ｓθ（ｔ）｝は、モデルパラメータ設定値θと時間ｔとが入力されると、初期時刻０から時間ｔ経過後の時刻ｔにおける単位ステップ応答Ｓθ（ｔ）を出力する。なお、単位ステップ応答とは、操作量ｕを単位ステップ入力とした場合における応答（つまり、制御対象プラント２０のプラント応答モデルの出力）のことである。
【００３４】
＜プラント応答関数｛Ｓθ（ｔ）｝の計算＞
次に、プラント応答関数｛Ｓθ（ｔ）｝の単位ステップ応答Ｓθ（ｔ）を計算する処理について、図３を参照しながら説明する。図３は、プラント応答関数｛Ｓθ（ｔ）｝の計算処理の一例を説明するためのフローチャートである。なお、図３では、モデルパラメータ設定値θと時間ｔとが入力されたものとする。
【００３５】
ここで、本実施形態では、プラント応答関数Ｓθ（ｔ）の計算モデルとして、制御量ｙについては過去Ｎ点の自己回帰、操作量ｕについては過去Ｍ点の移動平均を用いたＡＲＭＡ（autoregressive moving average）モデルを採用した場合について説明する。なお、Ｎ及びＭは、例えば、制御装置１０のユーザや管理者（以降、「ユーザ等」とも表す。）によって予め設定される。
【００３６】
このとき、モデルパラメータ設定値θを
【００３７】
【数１】
0006901037.tif 000002
とする。モデルパラメータ設定値θの各要素は数２で示す、ＡＲＭＡモデルの係数である。
【００３８】
【数２】
0006901037.tif 000003
ここで、ｋはインデックスである。
【００３９】
ステップＳ１１：操作量更新部１０５は、時刻を表す変数をτ、インデックスｋにおける状態ベクトルをφ（ｋ）として、時刻τと状態ベクトルφ（０）とを初期化する。また、操作量更新部１０５は、インデックスｋをｋ＝０に初期化する。
【００４０】
ここで、状態ベクトルφ（ｋ）は、
【００４１】
【数３】
0006901037.tif 000004
と表される。
【００４２】
操作量更新部１０５は、例えば、τ＝０と初期化すると共に、
【００４３】
【数４】
0006901037.tif 000005
と初期化する。すなわち、φ（０）はｕ（ｋ−１）に相当する要素のみ１、他の要素は０に初期化される。これは、単位ステップ応答を模擬する際の初期値として設定していることを意味する。
【００４４】
ステップＳ１２：次に、操作量更新部１０５は、モデルパラメータ設定値θと状態ベクトルφ（ｋ）とに基づいて制御量予測値ｙ（ｋ）を計算する。操作量更新部１０５は、例えば、ｙ（ｋ）＝φ（ｋ）Τθにより制御量予測値ｙ（ｋ）を計算する。ここで、Τは転置を表す。
【００４５】
ステップＳ１３：次に、操作量更新部１０５は、制御量予測値ｙ（ｋ）と状態ベクトルφ（ｋ）とに基づいて、次のインデックス（つまり、ｋ＋１）における状態ベクトルφ（ｋ＋１）を更新する。操作量更新部１０５は、例えば、以下により状態ベクトルφ（ｋ＋１）を更新する。
【００４６】
【数５】
0006901037.tif 000006
ここで、ｙ（ｋ）には上記のステップＳ１２で計算した制御量予測値を設定し、ｕ（ｋ）には１を設定する。また、ｙ（ｋ−１），・・・，ｙ（ｋ−Ｎ＋１）及びｕ（ｋ−１），・・・，ｕ（ｋ−Ｍ＋１）には、状態ベクトルφ（ｋ）と同じ値を設定する（すなわち、状態ベクトルφ（ｋ）に含まれるｙ（ｋ−１），・・・，ｙ（ｋ−Ｎ＋１）及びｕ（ｋ−１），・・・，ｕ（ｋ−Ｍ＋１）をそれぞれ設定する。）。
【００４７】
ステップＳ１４：次に、操作量更新部１０５は、時刻τをτ＋ΔＴに更新すると共に、インデックスｋをｋ＋１に更新する。ここで、ΔＴはＡＲＭＡモデルの１ステップの時間幅を表す。ΔＴは制御対象プラント２０の時定数に応じて任意に設定することができるが、例えば、制御周期Ｔｃと同一にすることが好適である。
【００４８】
ステップＳ１５：次に、操作量更新部１０５は、τ≧ｔであるか否かを判定する。そして、τ≧ｔであると判定されなかった場合（ステップＳ１５でＮＯ）、操作量更新部１０５は、ステップＳ１２に戻る。これにより、τ≧ｔとなるまで、ステップＳ１２〜ステップＳ１４が繰り返し実行される。
【００４９】
一方で、τ≧ｔであると判定された場合（ステップＳ１５でＹＥＳ）、操作量更新部１０５は、処理を終了する。これにより、最終的に計算されたｙ（ｋ）が単位ステップ応答Ｓθ（ｔ）として得られる（つまり、Ｓθ（ｔ）＝ｙ（ｋ）が、ステップ応答関数｛Ｓθ（ｔ）｝の時刻ｔにおける単位ステップ応答として計算される。）。
【００５０】
＜モデルパラメータθの推定＞
次に、プラント応答関数｛Ｓθ（ｔ）｝のモデルパラメータθを推定する処理（つまり、モデルパラメータ推定値θｅｓｔを算出する処理）について、図４を参照しながら説明する。図４は、モデルパラメータθの推定処理の一例を説明するためのフローチャートである。なお、図４では、制御量現在値ｙ０と操作量現在値ｕ０とモデルパラメータ設定値θとが入力されたものとして、或るインデックスｋにおけるモデルパラメータ推定値θｅｓｔを算出する場合について説明する。なお、該インデックスｋは前記プラント応答関数の計算におけるインデックスｋとは独立した値であり、モデルパラメータ推定処理の実行時インデックスを表す。
【００５１】
ステップＳ２１：モデルパラメータ推定部１０１は、モデルパラメータ推定値θｅｓｔ（ｋ）と共分散行列Ｐ（ｋ）とを初期化するか否かを判定する。ここで、初期化すると判定される場合としては、例えば、モデルパラメータ推定値θｅｓｔを初めて算出する場合（つまり、モデルパラメータ推定値θｅｓｔの初回算出時）、ユーザ等により初期化指示が行われた場合等が挙げられる。
【００５２】
モデルパラメータ推定値θｅｓｔ（ｋ）と共分散行列Ｐ（ｋ）とを初期化すると判定した場合（ステップＳ２１でＹＥＳ）、モデルパラメータ推定部１０１は、ステップＳ２２に進む。一方で、モデルパラメータ推定値θｅｓｔ（ｋ）と共分散行列Ｐ（ｋ）とを初期化すると判定しなかった場合（ステップＳ２１でＮＯ）、モデルパラメータ推定部１０１は、ステップＳ２３に進む。
【００５３】
ステップＳ２２：上記のステップＳ２１でモデルパラメータ推定値θｅｓｔ（ｋ）と共分散行列Ｐ（ｋ）とを初期化すると判定された場合、モデルパラメータ推定部１０１は、θｅｓｔ（０）＝θ０及びＰ（０）＝Ｉと初期化し、さらにモデル推定処理の実行時インデックスｋをｋ＝１と初期化する。ここで、θ０としては、既にプラント応答関数｛Ｓθ（ｔ）｝に設定されているモデルパラメータ設定値θを用いてもよいし、ユーザ等が事前に想定した初期値を用いてもよいし、全ての要素が０のベクトル等の固定の初期値を用いてもよい。また、Ｉとしては、単位行列としてもよいし、予め決められた任意の行列としてもよい。
【００５４】
ステップＳ２３：モデルパラメータ推定部１０１は、モデルパラメータ推定値θｅｓｔ（ｋ−１）と状態ベクトルφ（ｋ）と制御量現在値ｙ０とに基づいて、予測誤差ε（ｋ）を計算する。なお、モデルパラメータ推定値θｅｓｔ（ｋ−１）は前回（つまり、ｋ−１のとき）推定したモデルパラメータθの推定値である。また、状態ベクトルφ（ｋ）は前回更新した状態ベクトル（つまり、ｋ−１のときに図３のステップＳ１３で得られた状態ベクトル又はｋ−１のときに後述するステップＳ２７で得られた状態ベクトル）である。
【００５５】
モデルパラメータ推定部１０１は、例えば、以下により予測誤差ε（ｋ）を計算する。
【００５６】
ｙ（ｋ）＝ｙ０
ε（ｋ）＝ｙ（ｋ）−φ（ｋ）Τθｅｓｔ（ｋ−１）
ステップＳ２４：次に、モデルパラメータ推定部１０１は、共分散行列Ｐ（ｋ）を更新する。モデルパラメータ推定部１０１は、例えば、以下の式（２）により共分散行列Ｐ（ｋ）を更新する。
【００５７】
【数６】
0006901037.tif 000007
ここで、Ｐ（ｋ−１）は前回（つまり、ｋ−１のとき）得られた共分散行列である。
【００５８】
ステップＳ２５：次に、モデルパラメータ推定部１０１は、モデルパラメータ推定値θｅｓｔを更新するか否かを判定する。ここで、モデルパラメータ推定値θｅｓｔを更新すると判定される場合としては、例えば、上記のステップＳ２２での初期化時から所定の期間が経過するまでモデルパラメータ推定値θｅｓｔが更新されなかった場合、ユーザ等により更新指示が行われた場合等が挙げられる。
【００５９】
モデルパラメータ推定値θｅｓｔを更新すると判定した場合（ステップＳ２５でＹＥＳ）、モデルパラメータ推定部１０１は、ステップＳ２６に進む。一方で、モデルパラメータ推定値θｅｓｔを更新すると判定しなかった場合（ステップＳ２５でＮＯ）、モデルパラメータ推定部１０１は、ステップＳ２６を実行せずに、ステップＳ２７に進む。
【００６０】
ステップＳ２６：モデルパラメータ推定部１０１は、モデルパラメータ推定値θｅｓｔ（ｋ）を更新する。モデルパラメータ推定部１０１は、例えば、以下の式（３）によりモデルパラメータ推定値θｅｓｔ（ｋ）を更新する。
【００６１】
【数７】
0006901037.tif 000008
ステップＳ２７：モデルパラメータ推定部１０１は、制御量現在値ｙ０と操作量現在値ｕ０と状態ベクトルφ（ｋ）とに基づいて、次のインデックス（つまり、ｋ＋１）における状態ベクトルφ（ｋ＋１）を更新する。モデルパラメータ推定部１０１は、例えば、以下により状態ベクトルφ（ｋ＋１）を更新する。
【００６２】
【数８】
0006901037.tif 000009
ここで、ｙ（ｋ）には制御量現在値ｙ０を設定し、ｕ（ｋ）には操作量現在値ｕ０を設定する。また、ｙ（ｋ−１），・・・，ｙ（ｋ−Ｎ＋１）及びｕ（ｋ−１），・・・，ｕ（ｋ−Ｍ＋１）には、状態ベクトルφ（ｋ）と同じ値を設定する（すなわち、状態ベクトルφ（ｋ）に含まれるｙ（ｋ−１），・・・，ｙ（ｋ−Ｎ＋１）及びｕ（ｋ−１），・・・，ｕ（ｋ−Ｍ＋１）をそれぞれ設定する。）。このように、共分散行列の更新とモデルパラメータ推定値θｅｓｔの算出とは逐次最小２乗法に基づく処理で実現される。ただし、逐次最小２乗法は一例であって、例えば、カルマンフィルタ法が用いられてもよい。
【００６３】
なお、上記のステップＳ２４における共分散行列Ｐ（ｋ）の更新の他の例として、モデルパラメータ推定部１０１は、上記の式（２）を計算した後に、以下の式（４）により共分散行列Ｐ（ｋ）を更新してもよい。
【００６４】
【数９】
0006901037.tif 000010
ここで、０≦α≦１である。また、ｄｉａｇ（｛Ｐｉｉ（ｋ）｝ｉ）は行列Ｐ（ｋ）の対角成分のみで構成される行列（非対角成分以外は０）を表す。
【００６５】
上記のαは対角化調整パラメータであり、１に近いほど共分散行列Ｐ（ｋ）の非対角成分の絶対値を減少させる効果が大きい。つまり、上記の式（４）は、共分散行列Ｐ（ｋ）の非対角成分の絶対値を減少させる補正を行っていることを意味する。このような補正により、共分散行列Ｐ（ｋ）の計算を安定化させることが可能（つまり、共分散行列Ｐ（ｋ）の計算をより早く収束させることができると共に、その計算精度を向上させることが可能）となる。
【００６６】
ステップＳ２８：モデルパラメータ推定部１０１は、モデルパラメータ推定処理の実行時インデックスｋをｋ＋１に更新する。
【００６７】
＜目標値先読み部１０２の動作＞
次に、目標値先読み部１０２の動作について、図５を参照しながら説明する。図５は、目標値先読み部１０２の動作の一例を説明するための図である。
【００６８】
図５に示すように、目標値先読み部１０２は、目標値時系列｛ｒ（ｔ）｝と、先読み長Ｔｐとが入力されると、現在時刻ｔから先読み長後の時刻ｔ＋Ｔｐにおける先読み目標値ｒ（ｔ＋Ｔｐ）を出力する。このように、目標値先読み部１０２は、目標値時系列｛ｒ（ｔ）｝のうち、先読み長後の時刻ｔ＋Ｔｐにおける目標値ｒ（ｔ＋Ｔｐ）を出力する。
【００６９】
なお、図５に示す例では、目標値時系列｛ｒ（ｔ）｝が直線によって表される場合を示しているが、これに限られない。目標値時系列｛ｒ（ｔ）｝は任意の曲線や矩形等によって表されてもよい。特に、目標値時系列｛ｒ（ｔ）｝は、時刻ｔに応じて周期的に変化する曲線によって表されてもよい。また、目標値時系列｛ｒ（ｔ）｝は、予め設定されていてもよいし、未来の目標値が随時更新されてもよい。更に、未来の目標値は、例えば、現在の目標値ｒ（ｔ）が一定で継続するものとして現在の目標値を用いてもよいし、目標値が一定の速度で変化するものとして現在の目標値ｒ（ｔ）が先読み長後の時刻ｔ＋Ｔｐに到達する目標値としてもよい。
【００７０】
＜先読み応答補正部１１１の動作＞
次に、先読み応答補正部１１１の動作について、図６を参照しながら説明する。図６は、先読み応答補正部１１１の動作の一例を説明するための図（その１）である。
【００７１】
図６に示すように、先読み応答補正部１１１は、プラント応答関数｛Ｓθ（ｔ）｝と、先読み目標偏差ｅ０（ｔ）と、先読み長Ｔｐと、操作変化量時系列｛ｄｕ（ｔ）｝とが入力されると、まず、過去の操作変化量ｄｕによって現在の制御量ｙ０がＴｐ後に変化すると予測される値を先読み応答補正値ｙｎ（ｔ）として算出する。なお、例えば、現在時刻をｔとした場合、過去の操作変化量ｄｕは、ｄｕ（ｔ−Ｔｃ）、ｄｕ（ｔ−２Ｔｃ）等と表される。
【００７２】
次に、先読み応答補正部１１１は、先読み目標偏差ｅ０（ｔ）を先読み応答補正量ｙｎ（ｔ）によって補正した補正目標偏差ｅ＊（ｔ）を算出し、算出した補正目標偏差ｅ＊（ｔ）を出力する。図６に示すように、補正目標偏差ｅ＊（ｔ）は、ｅ＊（ｔ）＝ｒ（ｔ＋Ｔｐ）−（ｙ０（ｔ）＋ｙｎ（ｔ））＝ｅ０（ｔ）−ｙｎ（ｔ）で表される。
【００７３】
以降では、先読み応答補正値ｙｎ（ｔ）の計算方法について説明する。操作変化量時系列｛ｄｕ（ｔ）｝のみによる先読み時刻ｔ＋Ｔｐにおける制御量ｙの予測値（先読み応答予測値）ｙｎ，Ａ（ｔ）を、
【００７４】
【数１０】
0006901037.tif 000011
とする。ここで、Ｍ´は先読み応答予測値の計算に使用するモデルの長さ（モデル区間）である。
【００７５】
また、操作変化量時系列｛ｄｕ（ｔ）｝のみによる現在時刻ｔにおける制御量ｙの予測値（自由応答予測値）ｙｎ，Ｂ（ｔ）を、
【００７６】
【数１１】
0006901037.tif 000012
とする。ここで、上記と同様に、Ｍ´は自由応答予測値の計算に使用するモデル区間である。
【００７７】
そして、先読み応答補正部１１１は、先読み応答予測値ｙｎ，Ａ（ｔ）と、自由応答予測値ｙｎ，Ｂ（ｔ）との差を先読み応答補正値ｙｎ（ｔ）とする。すなわち、先読み応答補正部１１１は、ｙｎ（ｔ）＝ｙｎ，Ａ（ｔ）−ｙｎ，Ｂ（ｔ）とする。
【００７８】
また、一例として、予測時系列を記憶する予測時系列記憶部１２１を先読み応答補正部１１１が利用して先読み応答補正値ｙｎ（ｔ）を算出する場合について、図７を参照しながら説明する。図７は、先読み応答補正部１１１の動作の一例を説明するための図（その２）である。なお、予測時系列記憶部１２１は、例えば、補助記憶装置やＲＡＭ（Random Access Memory）等の記憶装置を用いて実現可能である。
【００７９】
図７に示すように、予測時系列記憶部１２１には、現在時刻をｔとして、時刻ｔ−Δｔから未来の時刻ｔ＋Ｔｂまでの制御量ｙの予測値ｙｎ，Ｃ（ｔ−Δｔ｜ｔ），ｙｎ，Ｃ（ｔ｜ｔ），ｙｎ，Ｃ（ｔ＋Δｔ｜ｔ），・・・，ｙｎ，Ｃ（ｔ＋Ｔｂ｜ｔ）を、時系列として記憶している。なお、Ｔｂは、予測時系列記憶部１２１に記憶される予測値ｙｎ，Ｃの長さ（時間長）を決める定数であり、Ｔｂ＝Ｎ´Δｔ（Ｎ´は、任意の正の整数）であるものとする。また、Δｔは、予測間隔であり、Δｔ＝Ｔｃであるものとする。以降では、予測値ｙｎ，Ｃを「一般化予測値ｙｎ，Ｃ」とも表す。
【００８０】
時刻ｔで予測した時刻ｓにおける一般化予測値ｙｎ，Ｃ（ｓ｜ｔ）は、
【００８１】
【数１２】
0006901037.tif 000013
で算出される。
【００８２】
予測時系列記憶部１２１に記憶された一般化予測値ｙｎ，Ｃを用いると、上述した自由応答予測値ｙｎ，Ｂ（ｔ）は一般化予測値ｙｎ，Ｃ（ｔ｜ｔ）と一致し、上述した先読み応答予測値ｙｎ，Ａ（ｔ）は一般化予測値ｙｎ，Ｃ（ｔ＋Ｔｐ｜ｔ）と一致する。
【００８３】
したがって、先読み応答補正部１１１は、予測時系列記憶部１２１に記憶された一般化予測値ｙｎ，Ｃを用いて、ｙｎ（ｔ）＝ｙｎ，Ａ（ｔ）−ｙｎ，Ｂ（ｔ）＝ｙｎ，Ｃ（ｔ＋Ｔｐ｜ｔ）−ｙｎ，Ｃ（ｔ｜ｔ）により先読み応答補正値ｙｎ（ｔ）を算出することができる。これにより、補正目標偏差ｅ＊（ｔ）＝ｅ０（ｔ）−ｙｎ（ｔ）を算出及び出力することができる。このように、予測時系列記憶部１２１を用いることで、先読み応答補正部１１１は、少ない計算量と少ないメモリとで補正目標偏差ｅ＊を計算することができる。
【００８４】
また、予測時系列記憶部１２１は、操作変化量ｄｕ（ｔ）が先読み応答補正部１１１に入力される度に、予測時系列更新部１２２によって更新される。ここで、ｍ´を自然数として、時刻ｔで予測したｍ´Δｔ先における一般化予測値ｙｎ，Ｃ（ｔ＋ｍ´Δｔ｜ｔ）は、
【００８５】
【数１３】
0006901037.tif 000014
で表される。このため、時刻ｔ＋Δｔで予測したｍ´Δｔ先における一般化予測値ｙｎ，Ｃ（ｔ＋Δｔ＋ｍ´Δｔ｜ｔ＋Δｔ）は、
【００８６】
【数１４】
0006901037.tif 000015
となる。すなわち、時刻ｔ＋Δｔで予測したｍ´Δｔ先における一般化予測値ｙｎ，Ｃ（ｔ＋Δｔ＋ｍ´Δｔ｜ｔ＋Δｔ）は、ｙｎ，Ｃ（ｔ＋Δｔ＋ｍ´Δｔ｜ｔ＋Δｔ）＝Ｓθ（ｍ´Δｔ）ｄｕ（ｔ＋Δｔ）＋ｙｎ，Ｃ（ｔ＋（ｍ´＋１）Δｔ｜ｔ）となる。
【００８７】
したがって、時刻ｔ＋Δｔで予測したｍ´Δｔ先における一般化予測値ｙｎ，Ｃ（ｔ＋Δｔ＋ｍ´Δｔ｜ｔ＋Δｔ）は、時刻ｔで予測した（ｍ´＋１）Δｔ先の一般化予測値ｙｎ，Ｃに対して、時刻ｔ＋Δｔにおける操作変化量ｄｕ（ｔ＋Δｔ）による影響を加えた値として更新される。
【００８８】
＜予測時系列更新部１２２が実行する更新処理＞
ここで、予測時系列更新部１２２が予測時系列記憶部１２１に保持されたデータを更新する処理について、図８を参照しながら説明する。図８は、予測時系列更新部１２２が実行する更新処理の一例を説明するためのフローチャートである。以降では、予測時系列記憶部１２１に記憶されている時刻ｔでの一般化予測値ｙｎ，Ｃを、時刻ｔ＋Δｔでの一般化予測値ｙｎ，Ｃに更新する場合について説明する。
【００８９】
ステップＳ４１：予測時系列更新部１２２は、予測時系列記憶部１２１に記憶されている時刻ｔでの一般化予測値ｙｎ，Ｃを時刻Δｔだけシフトする。
【００９０】
例えば、図９（ａ）に示すように、ｙｎ，Ｃ（ｔ−Δ｜ｔ），ｙｎ，Ｃ（ｔ｜ｔ），ｙｎ，Ｃ（ｔ＋Δ｜ｔ），・・・，ｙｎ，Ｃ（ｔ＋Ｔｐ｜ｔ），・・・，ｙｎ，Ｃ（ｔ＋Ｔｂ｜ｔ）が予測時系列記憶部１２１に記憶されており、それぞれの相対位置ｍ´が−１，０，１，・・・，Ｎ´であるものとする。この場合、予測時系列更新部１２２は、相対位置ｍ´≧０の一般化予測値ｙｎ，Ｃの相対位置ｍ´を−１する。すなわち、予測時系列更新部１２２は、相対位置ｍ´＝０であるｙｎ，Ｃ（ｔ｜ｔ）を相対位置ｍ´＝−１に、相対位置ｍ´＝１であるｙｎ，Ｃ（ｔ＋Δｔ｜ｔ）を相対位置ｍ´＝０に順にシフトする。以降も同様に、予測時系列更新部１２２は、相対位置ｍ´＝Ｎ´であるｙｎ，Ｃ（ｔ＋Ｔｂ｜ｔ）を相対位置ｍ´＝Ｎ´−１にシフトするまで、各ｙｎ，Ｃの相対位置ｍ´を順にシフトする。
【００９１】
ステップＳ４２：予測時系列更新部１２２は、最終時刻（すなわち、相対位置ｍ´＝Ｎ´）における一般化予測値ｙｎ，Ｃを更新する。時刻ｔで予測した最終時刻ｔ＋Ｎ´Δｔ＋Δｔにおける一般化予測値ｙｎ，Ｃ（ｔ＋Ｎ´Δｔ＋Δｔ｜ｔ）は、例えば、時刻ｔ＋Ｎ´Δｔ以前における一般化予測値ｙｎ，Ｃを用いて、
【００９２】
【数１５】
0006901037.tif 000016
により推定する。これは、時刻ｔでの時刻ｔ＋Ｎ´Δｔ−（ｊ−１）Δｔにおける一般化予測値ｙｎ，Ｃに対して重みｗｊを掛けて和を取った式である。
【００９３】
一般化予測値ｙｎ，Ｃ（ｔ＋Ｎ´Δｔ＋Δｔ｜ｔ）は、これ以外にも、例えば、時刻ｔで予測した時刻ｔ＋Ｎ´Δｔの一般化予測値ｙｎ，Ｃをそのまま使う場合、すなわち、ｙｎ，Ｃ（ｔ＋Ｎ´Δｔ＋Δｔ｜ｔ）＝ｙｎ，Ｃ（ｔ＋Ｎ´Δｔ｜ｔ）とする場合も考えられる。また、例えば、速度変化を考慮して、時刻ｔ＋Ｎ´Δｔにおける一般化予測値ｙｎ，Ｃと、時刻ｔ＋Ｎ´Δｔ−Δｔの一般化予測値ｙｎ，Ｃとを外挿して、
【００９４】
【数１６】
0006901037.tif 000017
を用いる場合も考えられる。例えば、この数１６で、ｒ＝１とすれば、時刻ｔ＋Ｎ´Δｔにおける一般化予測値ｙｎ，Ｃと、時刻ｔ＋Ｎ´Δｔ−Δｔにおける一般化予測値ｙｎ，Ｃとの差がΔｔだけ継続すると推定することに相当する。
【００９５】
これにより、図９（ｂ）に示すように、予測時系列記憶部１２１の相対位置Ｎ´がｙｎ，Ｃ（ｔ＋Ｔｂ＋Δｔ｜ｔ）＝ｙｎ，Ｃ（ｔ＋Ｎ´Δｔ＋Δｔ｜ｔ）に更新される。
【００９６】
ステップＳ４３：予測時系列更新部１２２は、最新の操作変化量ｄｕ（ｔ＋Δｔ）の影響を反映する。すなわち、図９（ｃ）に示すように、予測時系列更新部１２２は、相対位置ｍ´＝０〜Ｎ´までの各一般化予測値ｙｎ，Ｃ（ｔ＋（ｍ´＋１）Δｔ｜ｔ）に対して、Ｓθ（ｍ´Δｔ）ｄｕ（ｔ＋Δｔ）を加える。
【００９７】
具体的には、相対位置ｍ´＝０の一般化予測値ｙｎ，Ｃ（ｔ＋Δｔ｜ｔ）に対しては、Ｓθ（０）ｄｕ（ｔ＋Δｔ）を加える。また、相対位置ｍ´＝１の一般化予測値ｙｎ，Ｃ（ｔ＋２Δｔ｜ｔ）に対しては、Ｓθ（Δｔ）ｄｕ（ｔ＋Δｔ）を加える。ｍ´≧３の場合も同様である。
【００９８】
ステップＳ４４：予測時系列更新部１２２は、予測時系列記憶部１２１に記憶されている各一般化予測値ｙｎ，Ｃを時刻ｔ＋Δｔの予測時系列とする。すなわち、上記のＳ４３により、ｍ´＝０〜Ｎ´までの各一般化予測値ｙｎ，Ｃは、ｙｎ，Ｃ（ｔ＋Δｔ＋ｍ´Δｔ｜ｔ＋Δｔ）＝Ｓθ（ｍ´Δｔ）ｄｕ（ｔ＋Δｔ）＋ｙｎ，Ｃ（ｔ＋（ｍ´＋１）Δｔ｜ｔ）と表せる。
【００９９】
そこで、ｔ´＝ｔ＋Δｔとすれば、図９（ｄ）に示すように、相対位置ｍ´＝−１の一般化予測値はｙｎ，Ｃ（ｔ´−Δｔ｜ｔ´）、相対位置ｍ´＝０の一般化予測値はｙｎ，Ｃ（ｔ´｜ｔ´）、相対位置ｍ´＝１の一般化予測値はｙｎ，Ｃ（ｔ´＋Δｔ｜ｔ´）と表せる。ｍ´≧２の場合も同様に、ｙｎ，Ｃ（ｔ´＋ｍ´Δｔ｜ｔ´）と表せる。これにより、予測時系列記憶部１２１に記憶されている時刻ｔでの一般化予測値ｙｎ，Ｃが、時刻ｔ＋Δｔでの一般化予測値ｙｎ，Ｃに更新される。
【０１００】
＜操作変化量算出部１１２の動作＞
次に、操作変化量算出部１１２の動作について、図１０を参照しながら説明する。図１０は、操作変化量算出部１１２の動作の一例を説明するための図である。
【０１０１】
図１０に示すように、操作変化量算出部１１２は、補正目標偏差ｅ＊（ｔ）が入力されると、この補正目標偏差ｅ＊（ｔ）に対して所定の制御ゲインを乗じることで操作変化量ｄｕ（ｔ）を算出し、算出したｄｕ（ｔ）を出力する。例えば、所定の制御ゲインとして積分ゲインｋＩを用いる場合、操作変化量ｄｕ（ｔ）は、ｄｕ（ｔ）＝ｋＩ×ｅ＊（ｔ）で算出される。
【０１０２】
ただし、補正目標偏差ｅ＊（ｔ）に対して所定の制御ゲインを乗じた結果、操作変化量ｄｕの上限値ｄｕｍａｘを超える場合、操作変化量算出部１１２は、ｄｕｍａｘを操作変化量ｄｕ（ｔ）とする。同様に、補正目標偏差ｅ＊（ｔ）に対して所定の制御ゲインを乗じた結果、操作変化量ｄｕの下限値ｄｕｍｉｎを下回る場合、操作変化量算出部１１２は、ｄｕｍｉｎを操作変化量ｄｕ（ｔ）とする。
【０１０３】
＜制御装置１０のハードウェア構成＞
次に、一実施形態に係る制御装置１０のハードウェア構成について、図１１を参照しながら説明する。図１１は、一実施形態に係る制御装置１０のハードウェア構成の一例を示す図である。
【０１０４】
図１１に示すように、一実施形態に係る制御装置１０は、入力装置２０１と、表示装置２０２と、外部Ｉ／Ｆ２０３と、通信Ｉ／Ｆ２０４と、ＲＯＭ（Read Only Memory）２０５と、ＲＡＭ２０６と、ＣＰＵ２０７と、補助記憶装置２０８とを有する。これら各ハードウェアは、バス２０９により相互に通信可能に接続されている。
【０１０５】
入力装置２０１は、例えば各種ボタンやタッチパネル、キーボード、マウス等であり、制御装置１０に各種の操作を入力するのに用いられる。表示装置２０２は、例えばディスプレイ等であり、制御装置１０による各種の処理結果を表示する。なお、制御装置１０は、入力装置２０１及び表示装置２０２の少なくとも一方を有していなくてもよい。
【０１０６】
外部Ｉ／Ｆ２０３は、外部装置とのインタフェースである。外部装置には、記録媒体２０３ａ等がある。制御装置１０は、外部Ｉ／Ｆ２０３を介して、記録媒体２０３ａの読み取りや書き込みを行うことができる。記録媒体２０３ａには、例えば、ＳＤメモリカード(SD memory card）やＵＳＢメモリ、ＣＤ（Compact Disk）、ＤＶＤ（Digital Versatile Disk）等がある。なお、制御装置１０が有する各機能部を実現する１以上のプログラムは、記録媒体２０３ａに格納されていてもよい。
【０１０７】
通信Ｉ／Ｆ２０４は、制御装置１０が他の装置とデータ通信を行うためのインタフェースである。なお、制御装置１０が有する各機能部を実現する１以上のプログラムは、通信Ｉ／Ｆ２０４を介して、所定のサーバ等から取得（ダウンロード）されてもよい。
【０１０８】
ＲＯＭ２０５は、電源を切ってもデータを保持することができる不揮発性の半導体メモリである。ＲＡＭ２０６は、プログラムやデータを一時保持する揮発性の半導体メモリである。ＣＰＵ２０７は、例えば補助記憶装置２０８やＲＯＭ２０５からプログラムやデータをＲＡＭ２０６に読み出して、各種処理を実行する演算装置である。
【０１０９】
補助記憶装置２０８は、例えばＨＤＤ（Hard Disk Drive）やＳＳＤ（Solid State Drive）等であり、プログラムやデータを格納している不揮発性のメモリである。補助記憶装置２０８に格納されているプログラムやデータには、例えば、制御装置１０が有する各機能部を実現する１以上のプログラムや基本ソフトウェアであるＯＳ（Operating System）、ＯＳ上で動作する各種アプリケーションプログラム等がある。
【０１１０】
一実施形態に係る制御装置１０は、図１１に示すハードウェア構成を有することにより、上述した各種処理を実現することができる。なお、図１１では、制御装置１０が１台のコンピュータで実現されている場合のハードウェア構成例を示したが、これに限られず、制御装置１０は複数台のコンピュータで実現されていてもよい。
【０１１１】
［第二の実施形態］
次に、第二の実施形態について説明する。第二の実施形態では、モデルパラメータ推定値θｅｓｔを用いて制御量を推定する場合について説明する。これにより、モデルパラメータ推定値θｅｓｔをプラント応答関数｛Ｓθ（ｔ）｝に設定した場合における制御量の推定値を得ることができるため、例えば、ユーザ等は、実際の制御量と比較することで、このモデルパラメータ推定値θｅｓｔをプラント応答関数｛Ｓθ（ｔ）｝に設定すべきか否かを判断することができるようになる。言い換えれば、ユーザ等は、制御量の推定値と実際の制御量とを比較することで、現在のプラント応答関数｛Ｓθ（ｔ）｝の良否を視覚的に容易に判断することができるようになる。
【０１１２】
なお、第二の実施形態では、主に、第一の実施形態との相違点について説明し、第一の実施形態と実質的に同一の構成要素については、その説明を省略又は簡略化する。
【０１１３】
＜制御装置１０の構成＞
まず、本実施形態に係る制御装置１０の構成について、図１２を参照しながら説明する。図１２は、第二の実施形態に係る制御装置１０の構成の一例を示す図である。
【０１１４】
図１２に示すように、本実施形態に係る制御装置１０は、更に、制御量推定部１０７を有する。制御量推定部１０７は、例えば、制御装置１０にインストールされた１以上のプログラムがＣＰＵ等に実行させる処理により実現される。
【０１１５】
制御量推定部１０７は、モデルパラメータ推定部１０１により出力されたモデルパラメータ推定値θｅｓｔと、制御量現在値ｙ０と、操作量現在値ｕ０とに基づいて、制御量の推定値ｙｅｓｔを算出し、出力する。また、制御量推定部１０７は、推定初期化フラグの値に応じて計算の初期化（具体的には状態ベクトルφ（ｋ）の初期化）を行う。なお、以降では、制御量の推定値ｙｅｓｔを「制御量推定値ｙｅｓｔ」とも表す。
【０１１６】
＜制御量の推定＞
次に、制御量を推定する処理（つまり、制御量推定値ｙｅｓｔを算出する処理）について、図１３を参照しながら説明する。図１３は、制御量推定処理の一例を説明するためのフローチャートである。なお、図１４では、モデルパラメータ推定値θｅｓｔと制御量現在値ｙ０と操作量現在値ｕ０とが入力されたものとして、或るインデックスｋにおける制御量推定値ｙｅｓｔを算出する場合について説明する。なお、該インデックスｋは前記プラント応答関数の計算におけるインデックスｋとは独立した値であり、制御量推定処理の実行時インデックスを表す。
【０１１７】
ステップＳ５１：制御量推定部１０７は、推定初期化フラグがＯＮであるか否かを判定する。なお、推定初期化フラグがＯＮである場合には、制御量の推定は行わずに、状態ベクトルφ（ｋ＋１）が初期化される。推定初期化フラグのＯＮ又はＯＦＦは、例えば、ユーザ等によって適宜設定されればよい。
【０１１８】
推定初期化フラグがＯＮであると判定した場合（ステップＳ５１でＹＥＳ）、制御量推定部１０７は、ステップＳ５２に進む。一方で、推定初期化フラグがＯＮであると判定しなかった場合（ステップＳ５１でＮＯ）、制御量推定部１０７は、ステップＳ５３に進む。
【０１１９】
ステップＳ５２：上記のステップＳ５１で推定初期化フラグがＯＮであると判定された場合、制御量推定部１０７は、状態ベクトルφ（ｋ＋１）を初期化して、処理を終了する。制御量推定部１０７は、例えば、以下により状態ベクトルφ（ｋ＋１）を初期化する。
【０１２０】
【数１７】
0006901037.tif 000018
ここで、ｙ（ｋ）には制御量現在値ｙ０を設定し、ｕ（ｋ）には操作量現在値ｕ０を設定する。また、ｙ（ｋ−１），・・・，ｙ（ｋ−Ｎ＋１）及びｕ（ｋ−１），・・・，ｕ（ｋ−Ｍ＋１）には、状態ベクトルφ（ｋ）と同じ値を設定する（すなわち、状態ベクトルφ（ｋ）に含まれるｙ（ｋ−１），・・・，ｙ（ｋ−Ｎ＋１）及びｕ（ｋ−１），・・・，ｕ（ｋ−Ｍ＋１）をそれぞれ設定する。）。
【０１２１】
推定初期化フラグがＯＮの状態が一定時間継続することで、上記のステップＳ５２により、状態ベクトルφ（ｋ＋１）の各要素が、計測部１０３によって計測された制御量及び操作量（つまり、制御量の観測値及び操作量の観測値）に初期化される。
【０１２２】
ステップＳ５３：上記のステップＳ５１で推定初期化フラグがＯＮであると判定されなかった場合、制御量推定部１０７は、モデルパラメータ推定値θｅｓｔと状態ベクトルφ（ｋ）とに基づいて制御量推定値ｙｅｓｔ（ｋ）を計算する。制御量推定部１０７は、例えば、ｙｅｓｔ（ｋ）＝φ（ｋ）Τθｅｓｔにより制御量推定値ｙｅｓｔ（ｋ）を計算する。
【０１２３】
ステップＳ５４：次に、制御量推定部１０７は、上記のステップＳ５３で計算された制御量推定値ｙｅｓｔ（ｋ）と、操作量現在値ｕ０と、状態ベクトルφ（ｋ）とに基づいて、次のインデックス（つまり、ｋ＋１）における状態ベクトルφ（ｋ＋１）を更新する。制御量推定部１０７は、例えば、以下により状態ベクトルφ（ｋ＋１）を更新する。
【０１２４】
【数１８】
0006901037.tif 000019
ここで、ｙ（ｋ）には制御量推定値ｙｅｓｔ（ｋ）を設定し、ｕ（ｋ）には操作量現在値ｕ０を設定する。また、ｙ（ｋ−１），・・・，ｙ（ｋ−Ｎ＋１）及びｕ（ｋ−１），・・・，ｕ（ｋ−Ｍ＋１）には、状態ベクトルφ（ｋ）と同じ値を設定する。このように、制御量推定値ｙｅｓｔ（ｋ）を用いて状態ベクトルφ（ｋ＋１）を更新することで、上記のステップＳ５３での制御量推定値の計算を、制御量の観測値及び操作量の観測値のうちの操作量の観測値のみから行うことが可能となる（つまり、制御量の観測値を用いずに制御量推定値を計算することが可能となる。）。
【０１２５】
ステップＳ５５：制御量推定部１０７は、制御量推定処理の実行時インデックスｋをｋ＋１に更新する。
【０１２６】
［第三の実施形態］
次に、第三の実施形態について説明する。第三の実施形態では、モデルパラメータ推定値θｅｓｔにより、プラント応答関数｛Ｓθ（ｔ）｝に設定されているモデルパラメータ設定値θを更新する場合について説明する。
【０１２７】
なお、第三の実施形態では、主に、第一の実施形態との相違点について説明し、第一の実施形態と実質的に同一の構成要素については、その説明を省略又は簡略化する。
【０１２８】
＜制御装置１０の構成＞
本実施形態に係る制御装置１０の構成について、図１４を参照しながら説明する。図１４は、第三の実施形態に係る制御装置１０の構成の一例を示す図である。
【０１２９】
図１４に示すように、本実施形態に係る制御装置１０は、更に、モデルパラメータ更新部１０８を有する。モデルパラメータ更新部１０８は、例えば、制御装置１０にインストールされた１以上のプログラムがＣＰＵ等に実行させる処理により実現される。
【０１３０】
モデルパラメータ更新部１０８は、モデル更新トリガの値に応じて、モデルパラメータ推定部１０１により出力されたモデルパラメータ推定値θｅｓｔをプラント応答関数｛Ｓθ（ｔ）｝に設定する。例えば、モデル更新トリガがＯＮである場合、モデルパラメータ更新部１０８は、モデルパラメータ推定値θｅｓｔをプラント応答関数｛Ｓθ（ｔ）｝に設定する（つまり、モデルパラメータ推定値θｅｓｔをモデルパラメータ設定値θとする。）。なお、例えば、モデル更新トリガがＯＦＦである場合は、モデルパラメータ更新部１０８は、何もしない。
【０１３１】
これにより、ユーザ等は、例えば、モデルパラメータ推定値θｅｓｔをプラント応答関数｛Ｓθ（ｔ）｝に設定することでより高い精度で制御が可能であると判断した場合には、モデル更新トリガをＯＮにすることで、より高い精度の制御を実現することが可能となる。
【０１３２】
［第四の実施形態］
次に、第四の実施形態について説明する。第四の実施形態では、外乱ｖも考慮してモデルパラメータ推定値θｅｓｔを算出する場合について説明する。
【０１３３】
なお、第四の実施形態では、主に、第一の実施形態との相違点について説明し、第一の実施形態と実質的に同一の構成要素については、その説明を省略又は簡略化する。
【０１３４】
＜制御装置１０の構成＞
本実施形態に係る制御装置１０の構成について、図１５を参照しながら説明する。図１５は、第四の実施形態に係る制御装置１０の構成の一例を示す図である。
【０１３５】
図１５に示すように、本実施形態に係る制御装置１０の構成は第一の実施形態と同様であるが、主に、モデルパラメータ推定部１０１及び計測部１０３の機能が異なる。
【０１３６】
本実施形態に係る計測部１０３は、外乱ｖも計測し、計測した外乱ｖをモデルパラメータ推定部１０１に出力する。なお、計測部１０３により計測された外乱ｖを「観測外乱ｖ」とも表す。
【０１３７】
また、本実施形態に係るモデルパラメータ推定部１０１は、制御量現在値ｙ０と、操作量現在値ｕ０と、観測外乱ｖと、モデルパラメータ設定値θとが入力されると、モデルパラメータの推定値θｅｓｔを算出し、出力する。ここで、観測外乱ｖを用いる場合、プラント応答関数Ｓθ（ｔ）の計算モデルとしては、例えば、ＡＲＭＡＸモデル（autoregressive moving average model with exogenous variables）を採用すればよい。一例として、制御量ｙについては過去Ｎ点の自己回帰、操作量ｕについては過去Ｍ点の移動平均、観測外乱ｖについては過去Ｌ点の攪乱項を用いた場合、プラント応答関数Ｓθ（ｔ）の計算モデルは、以下の式（５）で表される。なお、Ｌは、例えば、制御装置１０のユーザ等によって予め設定される。
【０１３８】
【数１９】
0006901037.tif 000020
なお、モデルパラメータ設定値θは、
【０１３９】
【数２０】
0006901037.tif 000021
と表される。
【０１４０】
したがって、状態ベクトルφ（ｋ）を
【０１４１】
【数２１】
0006901037.tif 000022
と読み替えることで、モデルパラメータ推定部１０１は、図４のステップＳ２１〜ステップＳ２７によりモデルパラメータ推定値θｅｓｔを算出することができる。ただし、上記のステップＳ２７で状態ベクトルφ（ｋ＋１）を更新する際には、ｖ（ｋ）には観測外乱ｖを設定する。これにより、外乱ｖも考慮したモデルパラメータθを推定することが可能となる。
【０１４２】
［実施例１］
次に、実施例１について説明する。実施例１では、第二の実施形態に係る制御装置１０によって制御対象プラント２０を制御する場合について説明する。
【０１４３】
まず、真のプラント応答は、ΔＴを時間メッシュとして、ｙ（ｔ）＝−ａ１ｙ（ｔ−ΔＴ）−ａ２ｙ（ｎ−２ΔＴ）＋ｂ１ｕ（ｎ−ΔＴ）＋ｂ２ｕ（ｎ−２ΔＴ）＋δ（ｔ）で表されるものとする。ここで、δ（ｔ）はノイズを表す。
【０１４４】
ΔＴ＝Ｔｃとすれば、プラント応答関数｛Ｓθ（ｔ）｝は２次のＡＲＭＡモデルとなり、ｙ（ｋ）＝−ａ１ｙ（ｋ−１）−ａ２ｙ（ｋ−２）＋ｂ１ｙ（ｋ−１）＋ｂ２ｙ（ｋ−２）と表せる。
【０１４５】
このとき、状態ベクトルは、
【０１４６】
【数２２】
0006901037.tif 000023
と表せる。
【０１４７】
モデルパラメータ推定値は、
【０１４８】
【数２３】
0006901037.tif 000024
と表せる。
【０１４９】
モデルパラメータの真値は、
【０１５０】
【数２４】
0006901037.tif 000025
と表せる。また、制御量推定値は、ｙｅｓｔ（ｋ）＝φ（ｋ）Τθ（ｋ−１）と表すことができる。
【０１５１】
実施例１では、モデルパラメータの真値は、
【０１５２】
【数２５】
0006901037.tif 000026
であるものとする。このときの制御対象プラント２０のステップ応答を図１６に示す。図１６に示すように、実施例１における制御対象プラント２０のステップ応答は、滑らかに増加する形状となっている。
【０１５３】
実施例１では、制御対象プラント２０をバッチ運転し、図１７に示すように、３つのバッチ（バッチ１、バッチ２、バッチ３）が連続するものとする。このとき、実施例１では、図１７に示すように、運転開始前の時刻ｔ０で共分散行列の初期化とモデルパラメータの初期化とを行う。
【０１５４】
モデルパラメータの初期値は、
【０１５５】
【数２６】
0006901037.tif 000027
とした。
【０１５６】
また、共分散行列の初期値は、
【０１５７】
【数２７】
0006901037.tif 000028
とした。
【０１５８】
更に、対角化調整パラメータはα＝０．７とした。
【０１５９】
実施例１では、図１７に示すように、制御装置１０を運用するものとする。すなわち、バッチ開始後は共分散行列の更新を継続して行うと共に、バッチ１が開始され、或る程度共分散行列が更新された時刻ｔ１でモデルパラメータ推定値更新をＯＮにする（つまり、図４のステップＳ２５でＹＥＳとなるようにする）ことで、以降は継続してモデルパラメータ推定値を更新する。また、或る程度モデルパラメータ推定値が更新され、バッチ１が終了した時点ｔ２で推定初期化フラグをＯＮにし、状態ベクトルを初期化する。その後、推定初期化フラグをＯＦＦにして、バッチ２の間、バッチ２の制御量推定を行う。同様に、バッチ２が終了した時点ｔ３で推定初期化フラグをＯＮにし、状態ベクトルを初期化する。その後、推定初期化フラグをＯＦＦにして、バッチ３の間、バッチ３の制御量推定を行う。
【０１６０】
このとき、操作量ｕ及び制御量ｙが、それぞれ図１８（ａ）及び図１８（ｂ）に示すものであったとする。すなわち、実施例１では、操作量ｕに応じて制御量ｙが変動しており、制御量ｙに１つの山があるものとする。この場合、モデルパラメータ推定値θｅｓｔ及び制御量推定値ｙｅｓｔは、図１９に示すようになる。図１９の上段がモデルパラメータ推定値θｅｓｔ、下段が制御量推定値ｙｅｓｔである。なお、θ１＊、θ２＊、θ３＊及びθ４＊は、それぞれθ１の真値、θ２の真値、θ３の真値及びθ４の真値である。
【０１６１】
≪バッチ１≫
まず、時刻０〜５０までは共分散行列の更新を行い、モデルパラメータ推定値は更新しない。時刻５０からモデルパラメータ推定値の更新を開始しているが、モデルパラメータ推定値は、バッチ１の間、その値が変化し、バッチ１終了近くの時刻２３０付近ではほぼ真値に落ち着いている。バッチ１終了後、推定初期化フラグをＯＮにして、状態ベクトルを初期化する。
【０１６２】
≪バッチ２≫
モデルパラメータ推定値の更新は継続しているが、ほぼ一定に近い値で落ち着いている。また、制御量推定部１０７によって算出された制御量推定値ｙｅｓｔと、計測部１０３によって計測された制御量ｙ（つまり、制御量の観測値）とを比較すると、ピークの位置はほぼ一致しているものの、制御量推定値ｙｅｓｔが観測値ｙよりもやや高めの値となっている。
【０１６３】
≪バッチ３≫
モデルパラメータ推定値の更新は継続しているが、ほぼ一定に近い値で落ち着いている。バッチ２と比較すると変動量は少ない。また、制御量推定部１０７によって算出された制御量推定値ｙｅｓｔと、計測部１０３によって計測された制御量ｙ（観測値ｙ）とを比較すると、ピークの位置はほぼ一致している。バッチ２と比較すると、ピークの位置でｙｅｓｔが観測値ｙよりもやや低めの値となっているが、後半の時刻２００〜２３０では制御量推定値ｙｅｓｔと観測値ｙとがよく一致しているといえる。
【０１６４】
［実施例２］
次に、実施例２について説明する。実施例２では、操作量ｕ及び制御量ｙが実施例１と異なり、それ以外は実施例１と同様であるものとする。実施例２では、操作量ｕ及び制御量ｙが、それぞれ図２０（ａ）及び図２０（ｂ）に示すものであったとする。すなわち、実施例２では、操作量ｕに応じて制御量ｙが変動しており、制御量ｙに２つの山があるものとする。この場合、モデルパラメータ推定値θｅｓｔ及び制御量推定値ｙｅｓｔは、図２１に示すようになる。図２１の上段がモデルパラメータ推定値θｅｓｔ、下段が制御量推定値ｙｅｓｔである。なお、θ１＊、θ２＊、θ３＊及びθ４＊は、それぞれθ１の真値、θ２の真値、θ３の真値及びθ４の真値である。
【０１６５】
≪バッチ１≫
まず、時刻０〜５０までは共分散行列の更新を行い、モデルパラメータ推定値は更新しない。時刻５０からモデルパラメータ推定値の更新を開始しているが、モデルパラメータ推定値は、バッチ１の間、その値が変化し、バッチ１終了近くの時刻２３０付近ではほぼ真値に落ち着いている。バッチ１終了後、推定初期化フラグをＯＮにして、状態ベクトルを初期化する。
【０１６６】
≪バッチ２≫
モデルパラメータ推定値の更新は継続しているが、ほぼ一定に近い値で落ち着いている。また、制御量推定部１０７によって算出された制御量推定値ｙｅｓｔと、計測部１０３によって計測された制御量ｙとを比較すると、ピークの位置で若干のずれはあるものの、ほぼ近い値となっている。
【０１６７】
≪バッチ３≫
モデルパラメータ推定値の更新は継続しているが、ほぼ一定に近い値で落ち着いている。バッチ２と比較すると変動量は少ない。また、制御量推定部１０７によって算出された制御量推定値ｙｅｓｔと、計測部１０３によって計測された制御量ｙ（観測値ｙ）とを比較すると、ピークの位置で若干のずれはあるものの、ほぼ近い値となっている。
【０１６８】
本発明は、具体的に開示された上記の実施形態に限定されるものではなく、特許請求の範囲から逸脱することなく、種々の変形や変更、組み合わせ等が可能である。
【符号の説明】
【０１６９】
１０        制御装置
２０        制御対象プラント
１０１      モデルパラメータ推定部
１０２      目標値先読み部
１０３      計測部
１０４      差分器
１０５      操作量更新部
１０６      タイマ
１１１      先読み応答補正部
１１２      操作変化量算出部
１１３      加算器
【要約】      （修正有）【課題】エッジデバイス上でモデル予測制御を実現しつつ、モデル更新のためのパラメータを推定する方法を提供する。【解決手段】制御対象の制御量を目標値に追従させる制御装置であって、目標値時系列と、目標値を先読みする時間幅を示す先読み長とが入力されると、複数の目標値のうち、先読み目標値を取得する目標値先読み手段と、前記先読み目標値と、制御対象の現在の制御量との差である先読み目標偏差を算出する先読み目標偏差算出手段と、制御対象のプラント応答を表すプラント応答関数と、現在に至るまでの過去の操作量の変化量とに基づいて、先読み長後の制御量の予測値と、先読み目標値との差を示す補正目標偏差を算出する補正目標偏差算出手段と、前記補正目標偏差に基づいて、新たな操作量を算出する操作量算出手段と、制御量と、操作量とに基づいて、プラント応答関数に含まれるモデルパラメータの推定値を算出するモデルパラメータ推定手段を有す。【選択図】図１【図１】
0006901037.tif 000029
【図２】
0006901037.tif 000030
【図３】
0006901037.tif 000031
【図４】
0006901037.tif 000032
【図５】
0006901037.tif 000033
【図６】
0006901037.tif 000034
【図７】
0006901037.tif 000035
【図８】
0006901037.tif 000036
【図９】
0006901037.tif 000037
【図１０】
0006901037.tif 000038
【図１１】
0006901037.tif 000039
【図１２】
0006901037.tif 000040
【図１３】
0006901037.tif 000041
【図１４】
0006901037.tif 000042
【図１５】
0006901037.tif 000043
【図１６】
0006901037.tif 000044
【図１７】
0006901037.tif 000045
【図１８】
0006901037.tif 000046
【図１９】
0006901037.tif 000047
【図２０】
0006901037.tif 000048
【図２１】
0006901037.tif 000049