JP2022097222A - 制御装置、制御方法及びプログラム - Google Patents

Abstract

【課題】エッジデバイス上でモデル予測制御を実現しつつ、モデル更新のためのパラメータを推定する方法を提供する。【解決手段】制御対象の制御量を目標値に追従させる制御装置であって、目標値時系列と、目標値を先読みする時間幅を示す先読み長とが入力されると、複数の目標値のうち、先読み目標値を取得する目標値先読み手段と、前記先読み目標値と、制御対象の現在の制御量との差である先読み目標偏差を算出する先読み目標偏差算出手段と、制御対象のプラント応答を表すプラント応答関数と、現在に至るまでの過去の操作量の変化量とに基づいて、先読み長後の制御量の予測値と、先読み目標値との差を示す補正目標偏差を算出する補正目標偏差算出手段と、前記補正目標偏差に基づいて、新たな操作量を算出する操作量算出手段と、制御量と、操作量とに基づいて、プラント応答関数に含まれるモデルパラメータの推定値を算出するモデルパラメータ推定手段を有す。【選択図】図１

Description

本発明は、制御装置、制御方法及びプログラムに関する。

温調制御装置やＰＬＣ（Programmable Logic Controller）、ＤＣＳ（Distributed Control System）等の制御装置、パーソナルコンピュータや組み込み制御機器上で実装される制御装置等が産業上広く利用されている。

また、制御対象の制御量を目標値に追従させることを目的とする制御方式として、ＰＩＤ（Proportional-Integral-Differential）制御、モデル予測制御、内部モデル制御、ＬＱＧ（Linear-Quadratic-Gaussian）制御、Ｈ２制御、Ｈ∞制御等の各種の制御方式が知られている。

モデル予測制御は、制御対象の状態空間モデルや将来の時間応答モデルを用いた最適化計算を逐次的に行うことで望ましい応答を得る方式であり、産業界で広く用いられている（例えば、非特許文献１）。例えば、オンラインで数値最適化アルゴリズムを実行する標準的なモデル予測制御の産業応用として、空調システムの制御等への応用が知られている（例えば、特許文献１）。

また、現在に至るまでの過去の操作量の変化に応じた制御量の収束値の予測値と、目標値との差である補正目標偏差に基づいて新たな操作量を決定する制御装置が提案されている（例えば、特許文献２）。目標偏差の現在値と、操作量の変化量との関係を示す制御用論理式が成立する領域をグラフ上に表示する設計支援装置も提案されている（例えば、特許文献３）。

また、線形予測モデルに対してパラメータ同定を行う手法として、逐次最小２乗法、すなわち、ＲＬＳ（Recursive Least Squares）法やカルマンフィルタ法等が知られている（例えば、非特許文献２及び非特許文献３）。

更に、例えば、発電プラントの最適運転において、発電要素の物理モデルと入力状態量から特性パラメータを推定し、得られた特性パラメータから最適負荷配分を決定する制御方法も知られている。

特開２０１５－１０８４９９号公報 国際公開第２０１６／０９２８７２号 国際公開第２０１５／０６０１４９号

ヤン Ｍ． マチエヨフスキー，「モデル予測制御 制約のもとでの最適制御」，東京電機大学出版局，２００５ 足立 修一，「ＭＡＴＬＡＢによる制御のためのシステム同定」，東京電機大学出版局，１９９６ 足立 修一，「ＭＡＴＬＡＢによる制御のための上級システム同定」，東京電機大学出版局，２００４

しかしながら、従来のモデル予測制御は、逐次的な最適化計算を繰り返し実行する必要があるため、ＣＰＵ（Central Processing Unit）やメモリ等の計算資源が豊富なＰＣ（パーソナルコンピュータ）上で実行するものが多く、比較的計算資源が乏しいエッジデバイス（例えば、ＰＬＣ（Programmable Logic Controller）等の制御装置）上で実行することが困難であった。

また、モデル予測制御はモデルに基づく制御方式であるため、制御対象のモデルを予め同定する必要があるが、一度同定したモデルであっても、制御対象の経年変化や季節性変動等の要因によって変化していく場合がある。

本発明の一実施形態は、上記の点に鑑みてなされたもので、エッジデバイス上でモデル予測制御を実現しつつ、モデル更新のためのパラメータを推定することを目的とする。

上記目的を達成するため、本発明の一実施形態に係る制御装置は、制御対象に対する操作量を出力し、前記制御対象の制御量を目標値に追従させる制御装置であって、前記目標値の時系列である目標値時系列と、前記目標値を先読みする時間幅を示す先読み長とが入力されると、前記目標値時系列に含まれる複数の目標値のうち、前記先読み長後の目標値を示す先読み目標値を取得する目標値先読み手段と、前記先読み目標値と、前記制御対象の現在の制御量との差である先読み目標偏差を算出する先読み目標偏差算出手段と、前記制御対象のプラント応答を表すプラント応答関数と、現在に至るまでの過去の前記操作量の変化量とに基づいて、前記先読み長後の制御量の予測値と、前記先読み目標値との差を示す補正目標偏差を算出する補正目標偏差算出手段と、前記補正目標偏差に基づいて、新たな前記操作量を算出する操作量算出手段と、前記制御量と、前記操作量とに基づいて、前記プラント応答関数に含まれるモデルパラメータの推定値を算出するモデルパラメータ推定手段と、を有することを特徴とする。

エッジデバイス上でモデル予測制御を実現しつつ、モデル更新のためのパラメータを推定することができる。

第一の実施形態に係る制御装置の構成の一例を示す図である。 プラント応答関数の動作の一例を説明するための図である。 プラント応答関数の計算処理の一例を説明するためのフローチャートである。 モデルパラメータの推定処理の一例を説明するためのフローチャートである。 目標値先読み部の動作の一例を説明するための図である。 先読み応答補正部の動作の一例を説明するための図（その１）である。 先読み応答補正部の動作の一例を説明するための図（その２）である。 予測時系列更新部が実行する更新処理の一例を説明するためのフローチャートである。 予測時系列の更新の一例を説明するための図である。 操作変化量算出部の動作の一例を説明するための図である。 一実施形態に係る制御装置のハードウェア構成の一例を示す図である。 第二の実施形態に係る制御装置の構成の一例を示す図である。 制御量推定処理の一例を説明するためのフローチャートである。 第三の実施形態に係る制御装置の構成の一例を示す図である。 第四の実施形態に係る制御装置の構成の一例を示す図である。 実施例１におけるステップ応答を説明するための図である。 実施例１におけるモデルパラメータ推定時の設定を説明するための図である。 実施例１における操作量及び制御量を説明するための図である。 実施例１におけるモデルパラメータ推定値及び制御量推定値を説明するための図である。 実施例２における操作量及び制御量を説明するための図である。 実施例２におけるモデルパラメータ推定値及び制御量推定値を説明するための図である。

以下、本発明の実施の形態について、図面を参照しながら詳細に説明する。以降では、現在から未来への目標値時系列が入力された場合において、制御量を目標値に追従させるための操作量を計算すると共に、モデル更新のためのモデルパラメータを推定する制御装置１０について説明する。なお、以降で説明する各実施形態に係る制御装置１０は、例えば、ＰＬＣ等のエッジデバイスである。

以降で説明する各実施形態に係る制御装置１０は、任意の目標値の時系列データである目標値時系列｛ｒ（ｔ）｝や制御対象プラント２０の状態等を示す制御量ｙ、制御対象プラント２０のプラント応答モデル等に基づいて制御対象プラント２０に対する操作量ｕを算出する共に、モデル更新のためのモデルパラメータを推定する。そして、制御装置１０は、この操作量ｕに応じた制御対象プラント２０の制御量ｙを計測し、目標値時系列｛ｒ（ｔ）｝や制御量ｙ、プラント応答モデル等に基づいて次の操作量ｕを算出すると共に、新たなモデルパラメータを推定する。このように、以降で説明する各実施形態に係る制御装置１０は、制御量ｙを目標値時系列｛ｒ（ｔ）｝に追従させるための操作量ｕの算出と、プラント応答モデルのモデルパラメータの推定とをオンライン実行中（つまり、制御対象プラント２０の制御中）に繰り返し実行する。

以降で説明する各実施形態に係る制御装置１０は、モデル予測制御の高い制御性能を実現しつつも、毎回のオンライン最適化計算を必要としないため、少ない計算量と少ないメモリとで実現することが可能になると共に、プラント応答モデルのモデルパラメータを推定することが可能となる。なお、各実施形態に係る制御装置１０は、例えば、既設の他の制御装置の操作量と制御量とを用いてプラント応答モデルのパラメータを推定することが可能である。このため、例えば、制御装置１０の設置直後は他の制御装置の操作量と制御量とを用いてプラント応答モデルの初期モデルパラメータを推定し、運用開始後は季節変動や制御対象プラント２０の経年劣化等の継時的変化に追従したモデルパラメータを推定する、といったことが可能である。

なお、制御量ｙとしては、例えば、制御対象プラント２０の温度、目標値ｒとしては、例えば、設定温度等が挙げられる。ただし、制御量ｙ及び目標値ｒは、温度及び設定温度に限られず、制御対象プラント２０における任意の制御量及び当該制御量の目標となる目標値を用いることができる。

［第一の実施形態］
＜制御装置１０の構成＞
まず、本実施形態に係る制御装置１０の構成について、図１を参照しながら説明する。図１は、第一の実施形態に係る制御装置１０の構成の一例を示す図である。

図１に示すように、本実施形態に係る制御装置１０は、モデルパラメータ推定部１０１と、目標値先読み部１０２と、計測部１０３と、差分器１０４と、操作量更新部１０５と、タイマ１０６とを有する。これら各機能部は、例えば、制御装置１０にインストールされた１以上のプログラムがＣＰＵ等に実行させる処理により実現される。

目標値先読み部１０２は、所定の制御周期Ｔ_ｃ毎に、目標値時系列｛ｒ（ｔ）｝と、先読み長Ｔ_ｐとが入力されると、現在時刻ｔから先読み長後の時刻ｔ＋Ｔ_ｐにおける目標値ｒ（ｔ＋Ｔ_ｐ）を出力する。先読み長Ｔ_ｐとは、目標値時系列｛ｒ（ｔ）｝のうち、先読みする目標値ｒ（ｔ＋Ｔ_ｐ）を決定するための時間長である。なお、以降では、目標値ｒ（ｔ＋Ｔ_ｐ）を「先読み目標値ｒ（ｔ＋Ｔ_ｐ）」とも表す。

計測部１０３は、制御周期Ｔ_ｃ毎に、制御対象プラント２０の制御量ｙを計測する。そして、計測部１０３は、計測した制御量ｙの最新の値を、制御量現在値ｙ_０として出力する。なお、制御対象プラント２０の制御量ｙは、操作量ｕと外乱ｖとに応じて決定される。外乱ｖとしては、例えば、制御量ｙが温度である場合には外気温の低下又は上昇等が挙げられる。

また、計測部１０３は、操作量更新部１０５から出力された操作量ｕを取得し、取得した操作量ｕの最新の値を、操作量現在値ｕ_０として出力する。

差分器１０４は、目標値先読み部１０２から出力された先読み目標値ｒ（ｔ＋Ｔ_ｐ）と、制御量現在値ｙ_０との差（偏差）を目標偏差ｅ_０として出力する。目標偏差ｅ_０は、ｅ_０＝ｒ（ｔ＋Ｔ_ｐ）－ｙ_０（ｔ）で算出される。なお、以降では、目標偏差ｅ_０を「先読み目標偏差ｅ_０」とも表す。

操作量更新部１０５は、制御周期Ｔ_ｃ毎に、制御対象プラント２０に対する操作量ｕを出力する。操作量更新部１０５には、先読み応答補正部１１１と、操作変化量算出部１１２と、加算器１１３とが含まれる。

先読み応答補正部１１１は、プラント応答関数｛Ｓ_θ（ｔ）｝と、先読み目標偏差ｅ_０（ｔ）と、先読み長Ｔ_ｐと、過去の操作量ｕの変化量ｄｕの時系列データである操作変化量時系列｛ｄｕ（ｔ）｝とに基づいて、先読み目標偏差ｅ_０（ｔ）を補正した補正目標偏差ｅ^＊（ｔ）を算出する。補正目標偏差ｅ^＊（ｔ）の算出方法の詳細については後述する。ここで、プラント応答関数｛Ｓ_θ（ｔ）｝とはモデルパラメータθを含む関数であり、制御対象プラント２０のプラント応答モデルである。なお、以降では、操作量ｕの変化量ｄｕを「操作変化量ｄｕ」とも表す。

操作変化量算出部１１２は、制御周期Ｔ_ｃ毎に、先読み応答補正部１１１により算出された補正目標偏差ｅ^＊（ｔ）に基づいて、操作変化量ｄｕを算出する。操作変化量算出部１１２は、例えば、ｄｕ（ｔ－３Ｔ_ｃ）、ｄｕ（ｔ－２Ｔ_ｃ）、ｄｕ（ｔ－Ｔ_ｃ）という順序で操作変化量ｄｕ（ｔ）を算出し、出力する。なお、操作変化量ｄｕは、制御周期Ｔ_ｃ毎に操作量ｕが変化した量である。

加算器１１３は、計測部１０３から出力された操作量現在値ｕ_０と、操作変化量算出部１１２から出力された操作変化量ｄｕとを加算して、新たな操作量ｕを算出する。そして、加算器１１３は、この操作量ｕを制御対象プラント２０に出力する。この操作量ｕは、ｕ（ｔ）＝ｕ_０＋ｄｕ（ｔ）＝ｕ（ｔ－Ｔ_ｃ）＋ｄｕ（ｔ）で算出される。

タイマ１０６は、制御周期Ｔ_ｃ毎に、目標値先読み部１０２と計測部１０３とを動作させる。なお、目標値先読み部１０２と計測部１０３とが制御周期Ｔ_ｃ毎に動作することで、モデルパラメータ推定部１０１及び操作量更新部１０５も制御周期Ｔ_ｃ毎に動作する。

モデルパラメータ推定部１０１は、制御量現在値ｙ_０と、操作量現在値ｕ_０と、現在のプラント応答関数｛Ｓ_θ（ｔ）｝に設定されているモデルパラメータθとが入力されると、プラント応答関数｛Ｓ_θ（ｔ）｝のモデルパラメータの推定値θ_ｅｓｔを算出し、出力する。なお、以降では、現在のプラント応答関数｛Ｓ_θ（ｔ）｝に設定されているモデルパラメータθを「モデルパラメータ現在値θ」とも表し、プラント応答関数｛Ｓ_θ（ｔ）｝のモデルパラメータの推定値θ_ｅｓｔを「モデルパラメータ推定値θ_ｅｓｔ」とも表す。

＜プラント応答関数｛Ｓ_θ（ｔ）｝の動作＞
次に、プラント応答関数｛Ｓ_θ（ｔ）｝の動作について、図２を参照しながら説明する。図２は、プラント応答関数｛Ｓ_θ（ｔ）｝の動作の一例を説明するための図である。

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

＜プラント応答関数｛Ｓ_θ（ｔ）｝の計算＞
次に、プラント応答関数｛Ｓ_θ（ｔ）｝の単位ステップ応答Ｓ_θ（ｔ）を計算する処理について、図３を参照しながら説明する。図３は、プラント応答関数｛Ｓ_θ（ｔ）｝の計算処理の一例を説明するためのフローチャートである。なお、図３では、モデルパラメータ設定値θと時間ｔとが入力されたものとする。

ここで、本実施形態では、プラント応答関数Ｓ_θ（ｔ）の計算モデルとして、制御量ｙについては過去Ｎ点の自己回帰、操作量ｕについては過去Ｍ点の移動平均を用いたＡＲＭＡ（autoregressive moving average）モデルを採用した場合について説明する。なお、Ｎ及びＭは、例えば、制御装置１０のユーザや管理者（以降、「ユーザ等」とも表す。）によって予め設定される。

このとき、モデルパラメータ設定値θを

とする。モデルパラメータ設定値θの各要素は数２で示す、ＡＲＭＡモデルの係数である。

ここで、ｋはインデックスである。

ステップＳ１１：操作量更新部１０５は、時刻を表す変数をτ、インデックスｋにおける状態ベクトルをφ（ｋ）として、時刻τと状態ベクトルφ（０）とを初期化する。また、操作量更新部１０５は、インデックスｋをｋ＝０に初期化する。

ここで、状態ベクトルφ（ｋ）は、

と表される。

操作量更新部１０５は、例えば、τ＝０と初期化すると共に、

と初期化する。すなわち、φ（０）はｕ（ｋ－１）に相当する要素のみ１、他の要素は０に初期化される。これは、単位ステップ応答を模擬する際の初期値として設定していることを意味する。

ステップＳ１２：次に、操作量更新部１０５は、モデルパラメータ設定値θと状態ベクトルφ（ｋ）とに基づいて制御量予測値ｙ（ｋ）を計算する。操作量更新部１０５は、例えば、ｙ（ｋ）＝φ（ｋ）^Τθにより制御量予測値ｙ（ｋ）を計算する。ここで、Τは転置を表す。

ステップＳ１３：次に、操作量更新部１０５は、制御量予測値ｙ（ｋ）と状態ベクトルφ（ｋ）とに基づいて、次のインデックス（つまり、ｋ＋１）における状態ベクトルφ（ｋ＋１）を更新する。操作量更新部１０５は、例えば、以下により状態ベクトルφ（ｋ＋１）を更新する。

ここで、ｙ（ｋ）には上記のステップＳ１２で計算した制御量予測値を設定し、ｕ（ｋ）には１を設定する。また、ｙ（ｋ－１），・・・，ｙ（ｋ－Ｎ＋１）及びｕ（ｋ－１），・・・，ｕ（ｋ－Ｍ＋１）には、状態ベクトルφ（ｋ）と同じ値を設定する（すなわち、状態ベクトルφ（ｋ）に含まれるｙ（ｋ－１），・・・，ｙ（ｋ－Ｎ＋１）及びｕ（ｋ－１），・・・，ｕ（ｋ－Ｍ＋１）をそれぞれ設定する。）。

ステップＳ１４：次に、操作量更新部１０５は、時刻τをτ＋ΔＴに更新すると共に、インデックスｋをｋ＋１に更新する。ここで、ΔＴはＡＲＭＡモデルの１ステップの時間幅を表す。ΔＴは制御対象プラント２０の時定数に応じて任意に設定することができるが、例えば、制御周期Ｔ_ｃと同一にすることが好適である。

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

一方で、τ≧ｔであると判定された場合（ステップＳ１５でＹＥＳ）、操作量更新部１０５は、処理を終了する。これにより、最終的に計算されたｙ（ｋ）が単位ステップ応答Ｓ_θ（ｔ）として得られる（つまり、Ｓ_θ（ｔ）＝ｙ（ｋ）が、ステップ応答関数｛Ｓ_θ（ｔ）｝の時刻ｔにおける単位ステップ応答として計算される。）。

＜モデルパラメータθの推定＞
次に、プラント応答関数｛Ｓ_θ（ｔ）｝のモデルパラメータθを推定する処理（つまり、モデルパラメータ推定値θ_ｅｓｔを算出する処理）について、図４を参照しながら説明する。図４は、モデルパラメータθの推定処理の一例を説明するためのフローチャートである。なお、図４では、制御量現在値ｙ_０と操作量現在値ｕ_０とモデルパラメータ設定値θとが入力されたものとして、或るインデックスｋにおけるモデルパラメータ推定値θ_ｅｓｔを算出する場合について説明する。なお、該インデックスｋは前記プラント応答関数の計算におけるインデックスｋとは独立した値であり、モデルパラメータ推定処理の実行時インデックスを表す。

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

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

ステップＳ２２：上記のステップＳ２１でモデルパラメータ推定値θ_ｅｓｔ（ｋ）と共分散行列Ｐ（ｋ）とを初期化すると判定された場合、モデルパラメータ推定部１０１は、θ_ｅｓｔ（０）＝θ_０及びＰ（０）＝Ｉと初期化し、さらにモデル推定処理の実行時インデックスｋをｋ＝１と初期化する。ここで、θ_０としては、既にプラント応答関数｛Ｓ_θ（ｔ）｝に設定されているモデルパラメータ設定値θを用いてもよいし、ユーザ等が事前に想定した初期値を用いてもよいし、全ての要素が０のベクトル等の固定の初期値を用いてもよい。また、Ｉとしては、単位行列としてもよいし、予め決められた任意の行列としてもよい。

ステップＳ２３：モデルパラメータ推定部１０１は、モデルパラメータ推定値θ_ｅｓｔ（ｋ－１）と状態ベクトルφ（ｋ）と制御量現在値ｙ_０とに基づいて、予測誤差ε（ｋ）を計算する。なお、モデルパラメータ推定値θ_ｅｓｔ（ｋ－１）は前回（つまり、ｋ－１のとき）推定したモデルパラメータθの推定値である。また、状態ベクトルφ（ｋ）は前回更新した状態ベクトル（つまり、ｋ－１のときに図３のステップＳ１３で得られた状態ベクトル又はｋ－１のときに後述するステップＳ２７で得られた状態ベクトル）である。

モデルパラメータ推定部１０１は、例えば、以下により予測誤差ε（ｋ）を計算する。

ｙ（ｋ）＝ｙ_０
ε（ｋ）＝ｙ（ｋ）－φ（ｋ）^Τθ_ｅｓｔ（ｋ－１）
ステップＳ２４：次に、モデルパラメータ推定部１０１は、共分散行列Ｐ（ｋ）を更新する。モデルパラメータ推定部１０１は、例えば、以下の式（２）により共分散行列Ｐ（ｋ）を更新する。

ここで、Ｐ（ｋ－１）は前回（つまり、ｋ－１のとき）得られた共分散行列である。

ステップＳ２５：次に、モデルパラメータ推定部１０１は、モデルパラメータ推定値θ_ｅｓｔを更新するか否かを判定する。ここで、モデルパラメータ推定値θ_ｅｓｔを更新すると判定される場合としては、例えば、上記のステップＳ２２での初期化時から所定の期間が経過するまでモデルパラメータ推定値θ_ｅｓｔが更新されなかった場合、ユーザ等により更新指示が行われた場合等が挙げられる。

モデルパラメータ推定値θ_ｅｓｔを更新すると判定した場合（ステップＳ２５でＹＥＳ）、モデルパラメータ推定部１０１は、ステップＳ２６に進む。一方で、モデルパラメータ推定値θ_ｅｓｔを更新すると判定しなかった場合（ステップＳ２５でＮＯ）、モデルパラメータ推定部１０１は、ステップＳ２６を実行せずに、ステップＳ２７に進む。

ステップＳ２６：モデルパラメータ推定部１０１は、モデルパラメータ推定値θ_ｅｓｔ（ｋ）を更新する。モデルパラメータ推定部１０１は、例えば、以下の式（３）によりモデルパラメータ推定値θ_ｅｓｔ（ｋ）を更新する。

ステップＳ２７：モデルパラメータ推定部１０１は、制御量現在値ｙ_０と操作量現在値ｕ_０と状態ベクトルφ（ｋ）とに基づいて、次のインデックス（つまり、ｋ＋１）における状態ベクトルφ（ｋ＋１）を更新する。モデルパラメータ推定部１０１は、例えば、以下により状態ベクトルφ（ｋ＋１）を更新する。

ここで、ｙ（ｋ）には制御量現在値ｙ_０を設定し、ｕ（ｋ）には操作量現在値ｕ_０を設定する。また、ｙ（ｋ－１），・・・，ｙ（ｋ－Ｎ＋１）及びｕ（ｋ－１），・・・，ｕ（ｋ－Ｍ＋１）には、状態ベクトルφ（ｋ）と同じ値を設定する（すなわち、状態ベクトルφ（ｋ）に含まれるｙ（ｋ－１），・・・，ｙ（ｋ－Ｎ＋１）及びｕ（ｋ－１），・・・，ｕ（ｋ－Ｍ＋１）をそれぞれ設定する。）。このように、共分散行列の更新とモデルパラメータ推定値θ_ｅｓｔの算出とは逐次最小２乗法に基づく処理で実現される。ただし、逐次最小２乗法は一例であって、例えば、カルマンフィルタ法が用いられてもよい。

なお、上記のステップＳ２４における共分散行列Ｐ（ｋ）の更新の他の例として、モデルパラメータ推定部１０１は、上記の式（２）を計算した後に、以下の式（４）により共分散行列Ｐ（ｋ）を更新してもよい。

ここで、０≦α≦１である。また、ｄｉａｇ（｛Ｐ_ｉｉ（ｋ）｝_ｉ）は行列Ｐ（ｋ）の対角成分のみで構成される行列（非対角成分以外は０）を表す。

上記のαは対角化調整パラメータであり、１に近いほど共分散行列Ｐ（ｋ）の非対角成分の絶対値を減少させる効果が大きい。つまり、上記の式（４）は、共分散行列Ｐ（ｋ）の非対角成分の絶対値を減少させる補正を行っていることを意味する。このような補正により、共分散行列Ｐ（ｋ）の計算を安定化させることが可能（つまり、共分散行列Ｐ（ｋ）の計算をより早く収束させることができると共に、その計算精度を向上させることが可能）となる。

ステップＳ２８：モデルパラメータ推定部１０１は、モデルパラメータ推定処理の実行時インデックスｋをｋ＋１に更新する。

＜目標値先読み部１０２の動作＞
次に、目標値先読み部１０２の動作について、図５を参照しながら説明する。図５は、目標値先読み部１０２の動作の一例を説明するための図である。

図５に示すように、目標値先読み部１０２は、目標値時系列｛ｒ（ｔ）｝と、先読み長Ｔ_ｐとが入力されると、現在時刻ｔから先読み長後の時刻ｔ＋Ｔ_ｐにおける先読み目標値ｒ（ｔ＋Ｔ_ｐ）を出力する。このように、目標値先読み部１０２は、目標値時系列｛ｒ（ｔ）｝のうち、先読み長後の時刻ｔ＋Ｔ_ｐにおける目標値ｒ（ｔ＋Ｔ_ｐ）を出力する。

なお、図５に示す例では、目標値時系列｛ｒ（ｔ）｝が直線によって表される場合を示しているが、これに限られない。目標値時系列｛ｒ（ｔ）｝は任意の曲線や矩形等によって表されてもよい。特に、目標値時系列｛ｒ（ｔ）｝は、時刻ｔに応じて周期的に変化する曲線によって表されてもよい。また、目標値時系列｛ｒ（ｔ）｝は、予め設定されていてもよいし、未来の目標値が随時更新されてもよい。更に、未来の目標値は、例えば、現在の目標値ｒ（ｔ）が一定で継続するものとして現在の目標値を用いてもよいし、目標値が一定の速度で変化するものとして現在の目標値ｒ（ｔ）が先読み長後の時刻ｔ＋Ｔ_ｐに到達する目標値としてもよい。

＜先読み応答補正部１１１の動作＞
次に、先読み応答補正部１１１の動作について、図６を参照しながら説明する。図６は、先読み応答補正部１１１の動作の一例を説明するための図（その１）である。

図６に示すように、先読み応答補正部１１１は、プラント応答関数｛Ｓ_θ（ｔ）｝と、先読み目標偏差ｅ_０（ｔ）と、先読み長Ｔ_ｐと、操作変化量時系列｛ｄｕ（ｔ）｝とが入力されると、まず、過去の操作変化量ｄｕによって現在の制御量ｙ_０がＴ_ｐ後に変化すると予測される値を先読み応答補正値ｙ_ｎ（ｔ）として算出する。なお、例えば、現在時刻をｔとした場合、過去の操作変化量ｄｕは、ｄｕ（ｔ－Ｔ_ｃ）、ｄｕ（ｔ－２Ｔ_ｃ）等と表される。

次に、先読み応答補正部１１１は、先読み目標偏差ｅ_０（ｔ）を先読み応答補正量ｙ_ｎ（ｔ）によって補正した補正目標偏差ｅ^＊（ｔ）を算出し、算出した補正目標偏差ｅ^＊（ｔ）を出力する。図６に示すように、補正目標偏差ｅ^＊（ｔ）は、ｅ^＊（ｔ）＝ｒ（ｔ＋Ｔ_ｐ）－（ｙ_０（ｔ）＋ｙ_ｎ（ｔ））＝ｅ_０（ｔ）－ｙ_ｎ（ｔ）で表される。

以降では、先読み応答補正値ｙ_ｎ（ｔ）の計算方法について説明する。操作変化量時系列｛ｄｕ（ｔ）｝のみによる先読み時刻ｔ＋Ｔ_ｐにおける制御量ｙの予測値（先読み応答予測値）ｙ_ｎ，Ａ（ｔ）を、

とする。ここで、Ｍ´は先読み応答予測値の計算に使用するモデルの長さ（モデル区間）である。

また、操作変化量時系列｛ｄｕ（ｔ）｝のみによる現在時刻ｔにおける制御量ｙの予測値（自由応答予測値）ｙ_ｎ，Ｂ（ｔ）を、

とする。ここで、上記と同様に、Ｍ´は自由応答予測値の計算に使用するモデル区間である。

そして、先読み応答補正部１１１は、先読み応答予測値ｙ_ｎ，Ａ（ｔ）と、自由応答予測値ｙ_ｎ，Ｂ（ｔ）との差を先読み応答補正値ｙ_ｎ（ｔ）とする。すなわち、先読み応答補正部１１１は、ｙ_ｎ（ｔ）＝ｙ_ｎ，Ａ（ｔ）－ｙ_ｎ，Ｂ（ｔ）とする。

また、一例として、予測時系列を記憶する予測時系列記憶部１２１を先読み応答補正部１１１が利用して先読み応答補正値ｙ_ｎ（ｔ）を算出する場合について、図７を参照しながら説明する。図７は、先読み応答補正部１１１の動作の一例を説明するための図（その２）である。なお、予測時系列記憶部１２１は、例えば、補助記憶装置やＲＡＭ（Random Access Memory）等の記憶装置を用いて実現可能である。

図７に示すように、予測時系列記憶部１２１には、現在時刻をｔとして、時刻ｔ－Δｔから未来の時刻ｔ＋Ｔ_ｂまでの制御量ｙの予測値ｙ_ｎ，Ｃ（ｔ－Δｔ｜ｔ），ｙ_ｎ，Ｃ（ｔ｜ｔ），ｙ_ｎ，Ｃ（ｔ＋Δｔ｜ｔ），・・・，ｙ_ｎ，Ｃ（ｔ＋Ｔ_ｂ｜ｔ）を、時系列として記憶している。なお、Ｔ_ｂは、予測時系列記憶部１２１に記憶される予測値ｙ_ｎ，Ｃの長さ（時間長）を決める定数であり、Ｔ_ｂ＝Ｎ´Δｔ（Ｎ´は、任意の正の整数）であるものとする。また、Δｔは、予測間隔であり、Δｔ＝Ｔ_ｃであるものとする。以降では、予測値ｙ_ｎ，Ｃを「一般化予測値ｙ_ｎ，Ｃ」とも表す。

時刻ｔで予測した時刻ｓにおける一般化予測値ｙ_ｎ，Ｃ（ｓ｜ｔ）は、

で算出される。

予測時系列記憶部１２１に記憶された一般化予測値ｙ_ｎ，Ｃを用いると、上述した自由応答予測値ｙ_ｎ，Ｂ（ｔ）は一般化予測値ｙ_ｎ，Ｃ（ｔ｜ｔ）と一致し、上述した先読み応答予測値ｙ_ｎ，Ａ（ｔ）は一般化予測値ｙ_ｎ，Ｃ（ｔ＋Ｔ_ｐ｜ｔ）と一致する。

したがって、先読み応答補正部１１１は、予測時系列記憶部１２１に記憶された一般化予測値ｙ_ｎ，Ｃを用いて、ｙ_ｎ（ｔ）＝ｙ_ｎ，Ａ（ｔ）－ｙ_ｎ，Ｂ（ｔ）＝ｙ_ｎ，Ｃ（ｔ＋Ｔ_ｐ｜ｔ）－ｙ_ｎ，Ｃ（ｔ｜ｔ）により先読み応答補正値ｙ_ｎ（ｔ）を算出することができる。これにより、補正目標偏差ｅ^＊（ｔ）＝ｅ_０（ｔ）－ｙ_ｎ（ｔ）を算出及び出力することができる。このように、予測時系列記憶部１２１を用いることで、先読み応答補正部１１１は、少ない計算量と少ないメモリとで補正目標偏差ｅ^＊を計算することができる。

また、予測時系列記憶部１２１は、操作変化量ｄｕ（ｔ）が先読み応答補正部１１１に入力される度に、予測時系列更新部１２２によって更新される。ここで、ｍ´を自然数として、時刻ｔで予測したｍ´Δｔ先における一般化予測値ｙ_ｎ，Ｃ（ｔ＋ｍ´Δｔ｜ｔ）は、

で表される。このため、時刻ｔ＋Δｔで予測したｍ´Δｔ先における一般化予測値ｙ_ｎ，Ｃ（ｔ＋Δｔ＋ｍ´Δｔ｜ｔ＋Δｔ）は、

となる。すなわち、時刻ｔ＋Δｔで予測したｍ´Δｔ先における一般化予測値ｙ_ｎ，Ｃ（ｔ＋Δｔ＋ｍ´Δｔ｜ｔ＋Δｔ）は、ｙ_ｎ，Ｃ（ｔ＋Δｔ＋ｍ´Δｔ｜ｔ＋Δｔ）＝Ｓ_θ（ｍ´Δｔ）ｄｕ（ｔ＋Δｔ）＋ｙ_ｎ，Ｃ（ｔ＋（ｍ´＋１）Δｔ｜ｔ）となる。

したがって、時刻ｔ＋Δｔで予測したｍ´Δｔ先における一般化予測値ｙ_ｎ，Ｃ（ｔ＋Δｔ＋ｍ´Δｔ｜ｔ＋Δｔ）は、時刻ｔで予測した（ｍ´＋１）Δｔ先の一般化予測値ｙ_ｎ，Ｃに対して、時刻ｔ＋Δｔにおける操作変化量ｄｕ（ｔ＋Δｔ）による影響を加えた値として更新される。

＜予測時系列更新部１２２が実行する更新処理＞
ここで、予測時系列更新部１２２が予測時系列記憶部１２１に保持されたデータを更新する処理について、図８を参照しながら説明する。図８は、予測時系列更新部１２２が実行する更新処理の一例を説明するためのフローチャートである。以降では、予測時系列記憶部１２１に記憶されている時刻ｔでの一般化予測値ｙ_ｎ，Ｃを、時刻ｔ＋Δｔでの一般化予測値ｙ_ｎ，Ｃに更新する場合について説明する。

ステップＳ４１：予測時系列更新部１２２は、予測時系列記憶部１２１に記憶されている時刻ｔでの一般化予測値ｙ_ｎ，Ｃを時刻Δｔだけシフトする。

例えば、図９（ａ）に示すように、ｙ_ｎ，Ｃ（ｔ－Δ｜ｔ），ｙ_ｎ，Ｃ（ｔ｜ｔ），ｙ_ｎ，Ｃ（ｔ＋Δ｜ｔ），・・・，ｙ_ｎ，Ｃ（ｔ＋Ｔ_ｐ｜ｔ），・・・，ｙ_ｎ，Ｃ（ｔ＋Ｔ_ｂ｜ｔ）が予測時系列記憶部１２１に記憶されており、それぞれの相対位置ｍ´が－１，０，１，・・・，Ｎ´であるものとする。この場合、予測時系列更新部１２２は、相対位置ｍ´≧０の一般化予測値ｙ_ｎ，Ｃの相対位置ｍ´を－１する。すなわち、予測時系列更新部１２２は、相対位置ｍ´＝０であるｙ_ｎ，Ｃ（ｔ｜ｔ）を相対位置ｍ´＝－１に、相対位置ｍ´＝１であるｙ_ｎ，Ｃ（ｔ＋Δｔ｜ｔ）を相対位置ｍ´＝０に順にシフトする。以降も同様に、予測時系列更新部１２２は、相対位置ｍ´＝Ｎ´であるｙ_ｎ，Ｃ（ｔ＋Ｔ_ｂ｜ｔ）を相対位置ｍ´＝Ｎ´－１にシフトするまで、各ｙ_ｎ，Ｃの相対位置ｍ´を順にシフトする。

ステップＳ４２：予測時系列更新部１２２は、最終時刻（すなわち、相対位置ｍ´＝Ｎ´）における一般化予測値ｙ_ｎ，Ｃを更新する。時刻ｔで予測した最終時刻ｔ＋Ｎ´Δｔ＋Δｔにおける一般化予測値ｙ_ｎ，Ｃ（ｔ＋Ｎ´Δｔ＋Δｔ｜ｔ）は、例えば、時刻ｔ＋Ｎ´Δｔ以前における一般化予測値ｙ_ｎ，Ｃを用いて、

により推定する。これは、時刻ｔでの時刻ｔ＋Ｎ´Δｔ－（ｊ－１）Δｔにおける一般化予測値ｙ_ｎ，Ｃに対して重みｗ_ｊを掛けて和を取った式である。

一般化予測値ｙ_ｎ，Ｃ（ｔ＋Ｎ´Δｔ＋Δｔ｜ｔ）は、これ以外にも、例えば、時刻ｔで予測した時刻ｔ＋Ｎ´Δｔの一般化予測値ｙ_ｎ，Ｃをそのまま使う場合、すなわち、ｙ_ｎ，Ｃ（ｔ＋Ｎ´Δｔ＋Δｔ｜ｔ）＝ｙ_ｎ，Ｃ（ｔ＋Ｎ´Δｔ｜ｔ）とする場合も考えられる。また、例えば、速度変化を考慮して、時刻ｔ＋Ｎ´Δｔにおける一般化予測値ｙ_ｎ，Ｃと、時刻ｔ＋Ｎ´Δｔ－Δｔの一般化予測値ｙ_ｎ，Ｃとを外挿して、

を用いる場合も考えられる。例えば、この数１６で、ｒ＝１とすれば、時刻ｔ＋Ｎ´Δｔにおける一般化予測値ｙ_ｎ，Ｃと、時刻ｔ＋Ｎ´Δｔ－Δｔにおける一般化予測値ｙ_ｎ，Ｃとの差がΔｔだけ継続すると推定することに相当する。

これにより、図９（ｂ）に示すように、予測時系列記憶部１２１の相対位置Ｎ´がｙ_ｎ，Ｃ（ｔ＋Ｔ_ｂ＋Δｔ｜ｔ）＝ｙ_ｎ，Ｃ（ｔ＋Ｎ´Δｔ＋Δｔ｜ｔ）に更新される。

ステップＳ４３：予測時系列更新部１２２は、最新の操作変化量ｄｕ（ｔ＋Δｔ）の影響を反映する。すなわち、図９（ｃ）に示すように、予測時系列更新部１２２は、相対位置ｍ´＝０～Ｎ´までの各一般化予測値ｙ_ｎ，Ｃ（ｔ＋（ｍ´＋１）Δｔ｜ｔ）に対して、Ｓ_θ（ｍ´Δｔ）ｄｕ（ｔ＋Δｔ）を加える。

具体的には、相対位置ｍ´＝０の一般化予測値ｙ_ｎ，Ｃ（ｔ＋Δｔ｜ｔ）に対しては、Ｓ_θ（０）ｄｕ（ｔ＋Δｔ）を加える。また、相対位置ｍ´＝１の一般化予測値ｙ_ｎ，Ｃ（ｔ＋２Δｔ｜ｔ）に対しては、Ｓ_θ（Δｔ）ｄｕ（ｔ＋Δｔ）を加える。ｍ´≧３の場合も同様である。

ステップＳ４４：予測時系列更新部１２２は、予測時系列記憶部１２１に記憶されている各一般化予測値ｙ_ｎ，Ｃを時刻ｔ＋Δｔの予測時系列とする。すなわち、上記のＳ４３により、ｍ´＝０～Ｎ´までの各一般化予測値ｙ_ｎ，Ｃは、ｙ_ｎ，Ｃ（ｔ＋Δｔ＋ｍ´Δｔ｜ｔ＋Δｔ）＝Ｓ_θ（ｍ´Δｔ）ｄｕ（ｔ＋Δｔ）＋ｙ_ｎ，Ｃ（ｔ＋（ｍ´＋１）Δｔ｜ｔ）と表せる。

そこで、ｔ´＝ｔ＋Δｔとすれば、図９（ｄ）に示すように、相対位置ｍ´＝－１の一般化予測値はｙ_ｎ，Ｃ（ｔ´－Δｔ｜ｔ´）、相対位置ｍ´＝０の一般化予測値はｙ_ｎ，Ｃ（ｔ´｜ｔ´）、相対位置ｍ´＝１の一般化予測値はｙ_ｎ，Ｃ（ｔ´＋Δｔ｜ｔ´）と表せる。ｍ´≧２の場合も同様に、ｙ_ｎ，Ｃ（ｔ´＋ｍ´Δｔ｜ｔ´）と表せる。これにより、予測時系列記憶部１２１に記憶されている時刻ｔでの一般化予測値ｙ_ｎ，Ｃが、時刻ｔ＋Δｔでの一般化予測値ｙ_ｎ，Ｃに更新される。

＜操作変化量算出部１１２の動作＞
次に、操作変化量算出部１１２の動作について、図１０を参照しながら説明する。図１０は、操作変化量算出部１１２の動作の一例を説明するための図である。

図１０に示すように、操作変化量算出部１１２は、補正目標偏差ｅ^＊（ｔ）が入力されると、この補正目標偏差ｅ^＊（ｔ）に対して所定の制御ゲインを乗じることで操作変化量ｄｕ（ｔ）を算出し、算出したｄｕ（ｔ）を出力する。例えば、所定の制御ゲインとして積分ゲインｋ_Ｉを用いる場合、操作変化量ｄｕ（ｔ）は、ｄｕ（ｔ）＝ｋ_Ｉ×ｅ^＊（ｔ）で算出される。

ただし、補正目標偏差ｅ^＊（ｔ）に対して所定の制御ゲインを乗じた結果、操作変化量ｄｕの上限値ｄｕ_ｍａｘを超える場合、操作変化量算出部１１２は、ｄｕ_ｍａｘを操作変化量ｄｕ（ｔ）とする。同様に、補正目標偏差ｅ^＊（ｔ）に対して所定の制御ゲインを乗じた結果、操作変化量ｄｕの下限値ｄｕ_ｍｉｎを下回る場合、操作変化量算出部１１２は、ｄｕ_ｍｉｎを操作変化量ｄｕ（ｔ）とする。

＜制御装置１０のハードウェア構成＞
次に、一実施形態に係る制御装置１０のハードウェア構成について、図１１を参照しながら説明する。図１１は、一実施形態に係る制御装置１０のハードウェア構成の一例を示す図である。

図１１に示すように、一実施形態に係る制御装置１０は、入力装置２０１と、表示装置２０２と、外部Ｉ／Ｆ２０３と、通信Ｉ／Ｆ２０４と、ＲＯＭ（Read Only Memory）２０５と、ＲＡＭ２０６と、ＣＰＵ２０７と、補助記憶装置２０８とを有する。これら各ハードウェアは、バス２０９により相互に通信可能に接続されている。

入力装置２０１は、例えば各種ボタンやタッチパネル、キーボード、マウス等であり、制御装置１０に各種の操作を入力するのに用いられる。表示装置２０２は、例えばディスプレイ等であり、制御装置１０による各種の処理結果を表示する。なお、制御装置１０は、入力装置２０１及び表示装置２０２の少なくとも一方を有していなくてもよい。

外部Ｉ／Ｆ２０３は、外部装置とのインタフェースである。外部装置には、記録媒体２０３ａ等がある。制御装置１０は、外部Ｉ／Ｆ２０３を介して、記録媒体２０３ａの読み取りや書き込みを行うことができる。記録媒体２０３ａには、例えば、ＳＤメモリカード(SD memory card）やＵＳＢメモリ、ＣＤ（Compact Disk）、ＤＶＤ（Digital Versatile Disk）等がある。なお、制御装置１０が有する各機能部を実現する１以上のプログラムは、記録媒体２０３ａに格納されていてもよい。

通信Ｉ／Ｆ２０４は、制御装置１０が他の装置とデータ通信を行うためのインタフェースである。なお、制御装置１０が有する各機能部を実現する１以上のプログラムは、通信Ｉ／Ｆ２０４を介して、所定のサーバ等から取得（ダウンロード）されてもよい。

ＲＯＭ２０５は、電源を切ってもデータを保持することができる不揮発性の半導体メモリである。ＲＡＭ２０６は、プログラムやデータを一時保持する揮発性の半導体メモリである。ＣＰＵ２０７は、例えば補助記憶装置２０８やＲＯＭ２０５からプログラムやデータをＲＡＭ２０６に読み出して、各種処理を実行する演算装置である。

補助記憶装置２０８は、例えばＨＤＤ（Hard Disk Drive）やＳＳＤ（Solid State Drive）等であり、プログラムやデータを格納している不揮発性のメモリである。補助記憶装置２０８に格納されているプログラムやデータには、例えば、制御装置１０が有する各機能部を実現する１以上のプログラムや基本ソフトウェアであるＯＳ（Operating System）、ＯＳ上で動作する各種アプリケーションプログラム等がある。

一実施形態に係る制御装置１０は、図１１に示すハードウェア構成を有することにより、上述した各種処理を実現することができる。なお、図１１では、制御装置１０が１台のコンピュータで実現されている場合のハードウェア構成例を示したが、これに限られず、制御装置１０は複数台のコンピュータで実現されていてもよい。

［第二の実施形態］
次に、第二の実施形態について説明する。第二の実施形態では、モデルパラメータ推定値θ_ｅｓｔを用いて制御量を推定する場合について説明する。これにより、モデルパラメータ推定値θ_ｅｓｔをプラント応答関数｛Ｓ_θ（ｔ）｝に設定した場合における制御量の推定値を得ることができるため、例えば、ユーザ等は、実際の制御量と比較することで、このモデルパラメータ推定値θ_ｅｓｔをプラント応答関数｛Ｓ_θ（ｔ）｝に設定すべきか否かを判断することができるようになる。言い換えれば、ユーザ等は、制御量の推定値と実際の制御量とを比較することで、現在のプラント応答関数｛Ｓ_θ（ｔ）｝の良否を視覚的に容易に判断することができるようになる。

なお、第二の実施形態では、主に、第一の実施形態との相違点について説明し、第一の実施形態と実質的に同一の構成要素については、その説明を省略又は簡略化する。

＜制御装置１０の構成＞
まず、本実施形態に係る制御装置１０の構成について、図１２を参照しながら説明する。図１２は、第二の実施形態に係る制御装置１０の構成の一例を示す図である。

図１２に示すように、本実施形態に係る制御装置１０は、更に、制御量推定部１０７を有する。制御量推定部１０７は、例えば、制御装置１０にインストールされた１以上のプログラムがＣＰＵ等に実行させる処理により実現される。

制御量推定部１０７は、モデルパラメータ推定部１０１により出力されたモデルパラメータ推定値θ_ｅｓｔと、制御量現在値ｙ_０と、操作量現在値ｕ_０とに基づいて、制御量の推定値ｙ_ｅｓｔを算出し、出力する。また、制御量推定部１０７は、推定初期化フラグの値に応じて計算の初期化（具体的には状態ベクトルφ（ｋ）の初期化）を行う。なお、以降では、制御量の推定値ｙ_ｅｓｔを「制御量推定値ｙ_ｅｓｔ」とも表す。

＜制御量の推定＞
次に、制御量を推定する処理（つまり、制御量推定値ｙ_ｅｓｔを算出する処理）について、図１３を参照しながら説明する。図１３は、制御量推定処理の一例を説明するためのフローチャートである。なお、図１４では、モデルパラメータ推定値θ_ｅｓｔと制御量現在値ｙ_０と操作量現在値ｕ_０とが入力されたものとして、或るインデックスｋにおける制御量推定値ｙ_ｅｓｔを算出する場合について説明する。なお、該インデックスｋは前記プラント応答関数の計算におけるインデックスｋとは独立した値であり、制御量推定処理の実行時インデックスを表す。

ステップＳ５１：制御量推定部１０７は、推定初期化フラグがＯＮであるか否かを判定する。なお、推定初期化フラグがＯＮである場合には、制御量の推定は行わずに、状態ベクトルφ（ｋ＋１）が初期化される。推定初期化フラグのＯＮ又はＯＦＦは、例えば、ユーザ等によって適宜設定されればよい。

推定初期化フラグがＯＮであると判定した場合（ステップＳ５１でＹＥＳ）、制御量推定部１０７は、ステップＳ５２に進む。一方で、推定初期化フラグがＯＮであると判定しなかった場合（ステップＳ５１でＮＯ）、制御量推定部１０７は、ステップＳ５３に進む。

ステップＳ５２：上記のステップＳ５１で推定初期化フラグがＯＮであると判定された場合、制御量推定部１０７は、状態ベクトルφ（ｋ＋１）を初期化して、処理を終了する。制御量推定部１０７は、例えば、以下により状態ベクトルφ（ｋ＋１）を初期化する。

ここで、ｙ（ｋ）には制御量現在値ｙ_０を設定し、ｕ（ｋ）には操作量現在値ｕ_０を設定する。また、ｙ（ｋ－１），・・・，ｙ（ｋ－Ｎ＋１）及びｕ（ｋ－１），・・・，ｕ（ｋ－Ｍ＋１）には、状態ベクトルφ（ｋ）と同じ値を設定する（すなわち、状態ベクトルφ（ｋ）に含まれるｙ（ｋ－１），・・・，ｙ（ｋ－Ｎ＋１）及びｕ（ｋ－１），・・・，ｕ（ｋ－Ｍ＋１）をそれぞれ設定する。）。

推定初期化フラグがＯＮの状態が一定時間継続することで、上記のステップＳ５２により、状態ベクトルφ（ｋ＋１）の各要素が、計測部１０３によって計測された制御量及び操作量（つまり、制御量の観測値及び操作量の観測値）に初期化される。

ステップＳ５３：上記のステップＳ５１で推定初期化フラグがＯＮであると判定されなかった場合、制御量推定部１０７は、モデルパラメータ推定値θ_ｅｓｔと状態ベクトルφ（ｋ）とに基づいて制御量推定値ｙ_ｅｓｔ（ｋ）を計算する。制御量推定部１０７は、例えば、ｙ_ｅｓｔ（ｋ）＝φ（ｋ）^Τθ_ｅｓｔにより制御量推定値ｙ_ｅｓｔ（ｋ）を計算する。

ステップＳ５４：次に、制御量推定部１０７は、上記のステップＳ５３で計算された制御量推定値ｙ_ｅｓｔ（ｋ）と、操作量現在値ｕ_０と、状態ベクトルφ（ｋ）とに基づいて、次のインデックス（つまり、ｋ＋１）における状態ベクトルφ（ｋ＋１）を更新する。制御量推定部１０７は、例えば、以下により状態ベクトルφ（ｋ＋１）を更新する。

ここで、ｙ（ｋ）には制御量推定値ｙ_ｅｓｔ（ｋ）を設定し、ｕ（ｋ）には操作量現在値ｕ_０を設定する。また、ｙ（ｋ－１），・・・，ｙ（ｋ－Ｎ＋１）及びｕ（ｋ－１），・・・，ｕ（ｋ－Ｍ＋１）には、状態ベクトルφ（ｋ）と同じ値を設定する。このように、制御量推定値ｙ_ｅｓｔ（ｋ）を用いて状態ベクトルφ（ｋ＋１）を更新することで、上記のステップＳ５３での制御量推定値の計算を、制御量の観測値及び操作量の観測値のうちの操作量の観測値のみから行うことが可能となる（つまり、制御量の観測値を用いずに制御量推定値を計算することが可能となる。）。

ステップＳ５５：制御量推定部１０７は、制御量推定処理の実行時インデックスｋをｋ＋１に更新する。

［第三の実施形態］
次に、第三の実施形態について説明する。第三の実施形態では、モデルパラメータ推定値θ_ｅｓｔにより、プラント応答関数｛Ｓ_θ（ｔ）｝に設定されているモデルパラメータ設定値θを更新する場合について説明する。

なお、第三の実施形態では、主に、第一の実施形態との相違点について説明し、第一の実施形態と実質的に同一の構成要素については、その説明を省略又は簡略化する。

＜制御装置１０の構成＞
本実施形態に係る制御装置１０の構成について、図１４を参照しながら説明する。図１４は、第三の実施形態に係る制御装置１０の構成の一例を示す図である。

図１４に示すように、本実施形態に係る制御装置１０は、更に、モデルパラメータ更新部１０８を有する。モデルパラメータ更新部１０８は、例えば、制御装置１０にインストールされた１以上のプログラムがＣＰＵ等に実行させる処理により実現される。

モデルパラメータ更新部１０８は、モデル更新トリガの値に応じて、モデルパラメータ推定部１０１により出力されたモデルパラメータ推定値θ_ｅｓｔをプラント応答関数｛Ｓ_θ（ｔ）｝に設定する。例えば、モデル更新トリガがＯＮである場合、モデルパラメータ更新部１０８は、モデルパラメータ推定値θ_ｅｓｔをプラント応答関数｛Ｓ_θ（ｔ）｝に設定する（つまり、モデルパラメータ推定値θ_ｅｓｔをモデルパラメータ設定値θとする。）。なお、例えば、モデル更新トリガがＯＦＦである場合は、モデルパラメータ更新部１０８は、何もしない。

これにより、ユーザ等は、例えば、モデルパラメータ推定値θ_ｅｓｔをプラント応答関数｛Ｓ_θ（ｔ）｝に設定することでより高い精度で制御が可能であると判断した場合には、モデル更新トリガをＯＮにすることで、より高い精度の制御を実現することが可能となる。

［第四の実施形態］
次に、第四の実施形態について説明する。第四の実施形態では、外乱ｖも考慮してモデルパラメータ推定値θ_ｅｓｔを算出する場合について説明する。

なお、第四の実施形態では、主に、第一の実施形態との相違点について説明し、第一の実施形態と実質的に同一の構成要素については、その説明を省略又は簡略化する。

＜制御装置１０の構成＞
本実施形態に係る制御装置１０の構成について、図１５を参照しながら説明する。図１５は、第四の実施形態に係る制御装置１０の構成の一例を示す図である。

図１５に示すように、本実施形態に係る制御装置１０の構成は第一の実施形態と同様であるが、主に、モデルパラメータ推定部１０１及び計測部１０３の機能が異なる。

本実施形態に係る計測部１０３は、外乱ｖも計測し、計測した外乱ｖをモデルパラメータ推定部１０１に出力する。なお、計測部１０３により計測された外乱ｖを「観測外乱ｖ」とも表す。

また、本実施形態に係るモデルパラメータ推定部１０１は、制御量現在値ｙ_０と、操作量現在値ｕ_０と、観測外乱ｖと、モデルパラメータ設定値θとが入力されると、モデルパラメータの推定値θ_ｅｓｔを算出し、出力する。ここで、観測外乱ｖを用いる場合、プラント応答関数Ｓ_θ（ｔ）の計算モデルとしては、例えば、ＡＲＭＡＸモデル（autoregressive moving average model with exogenous variables）を採用すればよい。一例として、制御量ｙについては過去Ｎ点の自己回帰、操作量ｕについては過去Ｍ点の移動平均、観測外乱ｖについては過去Ｌ点の攪乱項を用いた場合、プラント応答関数Ｓ_θ（ｔ）の計算モデルは、以下の式（５）で表される。なお、Ｌは、例えば、制御装置１０のユーザ等によって予め設定される。

なお、モデルパラメータ設定値θは、

と表される。

したがって、状態ベクトルφ（ｋ）を

と読み替えることで、モデルパラメータ推定部１０１は、図４のステップＳ２１～ステップＳ２７によりモデルパラメータ推定値θ_ｅｓｔを算出することができる。ただし、上記のステップＳ２７で状態ベクトルφ（ｋ＋１）を更新する際には、ｖ（ｋ）には観測外乱ｖを設定する。これにより、外乱ｖも考慮したモデルパラメータθを推定することが可能となる。

［実施例１］
次に、実施例１について説明する。実施例１では、第二の実施形態に係る制御装置１０によって制御対象プラント２０を制御する場合について説明する。

まず、真のプラント応答は、ΔＴを時間メッシュとして、ｙ（ｔ）＝－ａ_１ｙ（ｔ－ΔＴ）－ａ_２ｙ（ｎ－２ΔＴ）＋ｂ_１ｕ（ｎ－ΔＴ）＋ｂ_２ｕ（ｎ－２ΔＴ）＋δ（ｔ）で表されるものとする。ここで、δ（ｔ）はノイズを表す。

ΔＴ＝Ｔ_ｃとすれば、プラント応答関数｛Ｓ_θ（ｔ）｝は２次のＡＲＭＡモデルとなり、ｙ（ｋ）＝－ａ_１ｙ（ｋ－１）－ａ_２ｙ（ｋ－２）＋ｂ_１ｙ（ｋ－１）＋ｂ_２ｙ（ｋ－２）と表せる。

このとき、状態ベクトルは、

と表せる。

モデルパラメータ推定値は、

と表せる。

モデルパラメータの真値は、

と表せる。また、制御量推定値は、ｙ_ｅｓｔ（ｋ）＝φ（ｋ）^Τθ（ｋ－１）と表すことができる。

実施例１では、モデルパラメータの真値は、

であるものとする。このときの制御対象プラント２０のステップ応答を図１６に示す。図１６に示すように、実施例１における制御対象プラント２０のステップ応答は、滑らかに増加する形状となっている。

実施例１では、制御対象プラント２０をバッチ運転し、図１７に示すように、３つのバッチ（バッチ１、バッチ２、バッチ３）が連続するものとする。このとき、実施例１では、図１７に示すように、運転開始前の時刻ｔ_０で共分散行列の初期化とモデルパラメータの初期化とを行う。

モデルパラメータの初期値は、

とした。

また、共分散行列の初期値は、

とした。

更に、対角化調整パラメータはα＝０．７とした。

実施例１では、図１７に示すように、制御装置１０を運用するものとする。すなわち、バッチ開始後は共分散行列の更新を継続して行うと共に、バッチ１が開始され、或る程度共分散行列が更新された時刻ｔ_１でモデルパラメータ推定値更新をＯＮにする（つまり、図４のステップＳ２５でＹＥＳとなるようにする）ことで、以降は継続してモデルパラメータ推定値を更新する。また、或る程度モデルパラメータ推定値が更新され、バッチ１が終了した時点ｔ_２で推定初期化フラグをＯＮにし、状態ベクトルを初期化する。その後、推定初期化フラグをＯＦＦにして、バッチ２の間、バッチ２の制御量推定を行う。同様に、バッチ２が終了した時点ｔ_３で推定初期化フラグをＯＮにし、状態ベクトルを初期化する。その後、推定初期化フラグをＯＦＦにして、バッチ３の間、バッチ３の制御量推定を行う。

このとき、操作量ｕ及び制御量ｙが、それぞれ図１８（ａ）及び図１８（ｂ）に示すものであったとする。すなわち、実施例１では、操作量ｕに応じて制御量ｙが変動しており、制御量ｙに１つの山があるものとする。この場合、モデルパラメータ推定値θ_ｅｓｔ及び制御量推定値ｙ_ｅｓｔは、図１９に示すようになる。図１９の上段がモデルパラメータ推定値θ_ｅｓｔ、下段が制御量推定値ｙ_ｅｓｔである。なお、θ_１ ^＊、θ_２ ^＊、θ_３ ^＊及びθ_４ ^＊は、それぞれθ_１の真値、θ_２の真値、θ_３の真値及びθ_４の真値である。

≪バッチ１≫
まず、時刻０～５０までは共分散行列の更新を行い、モデルパラメータ推定値は更新しない。時刻５０からモデルパラメータ推定値の更新を開始しているが、モデルパラメータ推定値は、バッチ１の間、その値が変化し、バッチ１終了近くの時刻２３０付近ではほぼ真値に落ち着いている。バッチ１終了後、推定初期化フラグをＯＮにして、状態ベクトルを初期化する。

≪バッチ２≫
モデルパラメータ推定値の更新は継続しているが、ほぼ一定に近い値で落ち着いている。また、制御量推定部１０７によって算出された制御量推定値ｙ_ｅｓｔと、計測部１０３によって計測された制御量ｙ（つまり、制御量の観測値）とを比較すると、ピークの位置はほぼ一致しているものの、制御量推定値ｙ_ｅｓｔが観測値ｙよりもやや高めの値となっている。

≪バッチ３≫
モデルパラメータ推定値の更新は継続しているが、ほぼ一定に近い値で落ち着いている。バッチ２と比較すると変動量は少ない。また、制御量推定部１０７によって算出された制御量推定値ｙ_ｅｓｔと、計測部１０３によって計測された制御量ｙ（観測値ｙ）とを比較すると、ピークの位置はほぼ一致している。バッチ２と比較すると、ピークの位置でｙ_ｅｓｔが観測値ｙよりもやや低めの値となっているが、後半の時刻２００～２３０では制御量推定値ｙ_ｅｓｔと観測値ｙとがよく一致しているといえる。

［実施例２］
次に、実施例２について説明する。実施例２では、操作量ｕ及び制御量ｙが実施例１と異なり、それ以外は実施例１と同様であるものとする。実施例２では、操作量ｕ及び制御量ｙが、それぞれ図２０（ａ）及び図２０（ｂ）に示すものであったとする。すなわち、実施例２では、操作量ｕに応じて制御量ｙが変動しており、制御量ｙに２つの山があるものとする。この場合、モデルパラメータ推定値θ_ｅｓｔ及び制御量推定値ｙ_ｅｓｔは、図２１に示すようになる。図２１の上段がモデルパラメータ推定値θ_ｅｓｔ、下段が制御量推定値ｙ_ｅｓｔである。なお、θ_１ ^＊、θ_２ ^＊、θ_３ ^＊及びθ_４ ^＊は、それぞれθ_１の真値、θ_２の真値、θ_３の真値及びθ_４の真値である。

≪バッチ１≫
まず、時刻０～５０までは共分散行列の更新を行い、モデルパラメータ推定値は更新しない。時刻５０からモデルパラメータ推定値の更新を開始しているが、モデルパラメータ推定値は、バッチ１の間、その値が変化し、バッチ１終了近くの時刻２３０付近ではほぼ真値に落ち着いている。バッチ１終了後、推定初期化フラグをＯＮにして、状態ベクトルを初期化する。

≪バッチ２≫
モデルパラメータ推定値の更新は継続しているが、ほぼ一定に近い値で落ち着いている。また、制御量推定部１０７によって算出された制御量推定値ｙ_ｅｓｔと、計測部１０３によって計測された制御量ｙとを比較すると、ピークの位置で若干のずれはあるものの、ほぼ近い値となっている。

≪バッチ３≫
モデルパラメータ推定値の更新は継続しているが、ほぼ一定に近い値で落ち着いている。バッチ２と比較すると変動量は少ない。また、制御量推定部１０７によって算出された制御量推定値ｙ_ｅｓｔと、計測部１０３によって計測された制御量ｙ（観測値ｙ）とを比較すると、ピークの位置で若干のずれはあるものの、ほぼ近い値となっている。

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

１０ 制御装置
２０ 制御対象プラント
１０１ モデルパラメータ推定部
１０２ 目標値先読み部
１０３ 計測部
１０４ 差分器
１０５ 操作量更新部
１０６ タイマ
１１１ 先読み応答補正部
１１２ 操作変化量算出部
１１３ 加算器

Claims (10)

1. 制御対象に対する操作量を出力し、前記制御対象の制御量を目標値に追従させる制御装置であって、
   前記目標値の時系列である目標値時系列と、前記目標値を先読みする時間幅を示す先読み長とが入力されると、前記目標値時系列に含まれる複数の目標値のうち、前記先読み長後の目標値を示す先読み目標値を取得する目標値先読み手段と、
   前記先読み目標値と、前記制御対象の現在の制御量との差である先読み目標偏差を算出する先読み目標偏差算出手段と、
   前記制御対象のプラント応答を表すプラント応答関数と、現在に至るまでの過去の前記操作量の変化量とに基づいて、前記先読み長後の制御量の予測値と、前記先読み目標値との差を示す補正目標偏差を算出する補正目標偏差算出手段と、
   前記補正目標偏差に基づいて、新たな前記操作量を算出する操作量算出手段と、
   前記制御量と、前記操作量とに基づいて、前記プラント応答関数に含まれるモデルパラメータの推定値を算出するモデルパラメータ推定手段と、
   を有することを特徴とする制御装置。
2. ＯＮ又はＯＦＦのいずれかの値を取るフラグがＯＮになった時点における前記制御量を初期値として、前記モデルパラメータ推定手段により算出された推定値と、前記操作量とに基づいて、前記時点以降の制御量の推定値を算出する制御量推定手段、
   を有することを特徴とする請求項１に記載の制御装置。
3. 前記制御量推定手段は、
   前記モデルパラメータ推定手段により算出された推定値を用いた前記制御量の推定値の算出と、前記プラント応答関数に含まれるモデルパラメータを用いた前記制御量の推定値の算出とを行う、
   ことを特徴とする請求項２に記載の制御装置。
4. 予め決められた更新トリガに応じて、前記モデルパラメータ推定手段により算出された推定値を、前記プラント応答関数に含まれるパラメータに設定するモデルパラメータ更新手段、
   を有することを特徴とする請求項１乃至３の何れか一項に記載の制御装置。
5. 前記モデルパラメータ推定手段は、
   前記制御量と、前記操作量と、前記制御対象に対する外乱の計測値とに基づいて、前記プラント応答関数に含まれるモデルパラメータの推定値を算出する、
   ことを特徴とする請求項１乃至４の何れか一項に記載の制御装置。
6. 前記プラント応答関数は、ＡＲＭＡモデル又はＡＲＭＡＸモデルから計算するステップ応答であり、
   前記モデルパラメータは、前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの係数であり、
   前記先読み長後の制御量の予測値は、前記先読み長後の時点におけるステップ応答を前記ＡＲＭＡモデル又はＡＲＭＡＸモデルにより計算することで算出される、ことを特徴とする請求項１乃至５の何れか一項に記載の制御装置。
7. 前記モデルパラメータ推定手段は、
   逐次最小２乗法又はカルマンフィルタの手法により、共分散行列の更新と前記モデルパラメータの推定値の算出とを行う、
   ことを特徴とする請求項１乃至６の何れか一項に記載の制御装置。
8. 前記モデルパラメータ推定手段は、
   前記共分散行列の更新において、前記共分散行列の非対角項の絶対値を減少させる補正処理を行う、
   ことを特徴とする請求項７に記載の制御装置。
9. 制御対象に対する操作量を出力し、前記制御対象の制御量を目標値に追従させるコンピュータが、
   前記目標値の時系列である目標値時系列と、前記目標値を先読みする時間幅を示す先読み長とが入力されると、前記目標値時系列に含まれる複数の目標値のうち、前記先読み長後の目標値を示す先読み目標値を取得する目標値先読み手順と、
   前記先読み目標値と、前記制御対象の現在の制御量との差である先読み目標偏差を算出する先読み目標偏差算出手順と、
   前記制御対象のプラント応答を表すプラント応答関数と、現在に至るまでの過去の前記操作量の変化量とに基づいて、前記先読み長後の制御量の予測値と、前記先読み目標値との差を示す補正目標偏差を算出する補正目標偏差算出手順と、
   前記補正目標偏差に基づいて、新たな前記操作量を算出する操作量算出手順と、
   前記制御量と、前記操作量とに基づいて、前記プラント応答関数に含まれるモデルパラメータの推定値を算出するモデルパラメータ推定手順と、
   を実行することを特徴とする制御方法。
10. 制御対象に対する操作量を出力し、前記制御対象の制御量を目標値に追従させるコンピュータを、
    前記目標値の時系列である目標値時系列と、前記目標値を先読みする時間幅を示す先読み長とが入力されると、前記目標値時系列に含まれる複数の目標値のうち、前記先読み長後の目標値を示す先読み目標値を取得する目標値先読み手段、
    前記先読み目標値と、前記制御対象の現在の制御量との差である先読み目標偏差を算出する先読み目標偏差算出手段、
    前記制御対象のプラント応答を表すプラント応答関数と、現在に至るまでの過去の前記操作量の変化量とに基づいて、前記先読み長後の制御量の予測値と、前記先読み目標値との差を示す補正目標偏差を算出する補正目標偏差算出手段、
    前記補正目標偏差に基づいて、新たな前記操作量を算出する操作量算出手段、
    前記制御量と、前記操作量とに基づいて、前記プラント応答関数に含まれるモデルパラメータの推定値を算出するモデルパラメータ推定手段、
    として機能させるためのプログラム。

Images