(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）モデルを採用した場合について説明する。なお、Ｎ及びＭは、例えば、制御装置１０のユーザや管理者（以降、「ユーザ等」とも表す。）によって予め設定される。ただし、これは一例であって、例えば、ＡＲＭＡＸモデル等の他のモデルを採用してもよい。
【００３９】
このとき、モデルパラメータ設定値θを
【００４０】
【数１】
0006927446.tif 000002
とする。モデルパラメータ設定値θの各要素は以下の数２で示す、ＡＲＭＡモデルの係数である。
【００４１】
【数２】
0006927446.tif 000003
ここで、ｋはインデックスである。ただし、これは一例であって、例えば、過去Ｌ点（Ｌは予め設定された１以上の整数）の外乱ｖを考慮したモデルが用いられてもよい。
【００４２】
ステップＳ１１：操作量更新部１０４は、時刻を表す変数をτ、インデックスｋにおける状態ベクトルをφ（ｋ）として、時刻τと状態ベクトルφ（０）とを初期化する。また、操作量更新部１０４は、インデックスｋをｋ＝０に初期化する。
【００４３】
ここで、状態ベクトルφ（ｋ）は、
【００４４】
【数３】
0006927446.tif 000004
と表される。
【００４５】
操作量更新部１０４は、例えば、τ＝０と初期化すると共に、
【００４６】
【数４】
0006927446.tif 000005
と初期化する。すなわち、φ（０）はｕ（ｋ−１）に相当する要素のみ１、他の要素は０に初期化される。これは、単位ステップ応答を模擬する際の初期値として設定していることを意味する。
【００４７】
ステップＳ１２：次に、操作量更新部１０４は、モデルパラメータ設定値θと状態ベクトルφ（ｋ）とに基づいて制御量予測値ｙ（ｋ）を計算する。操作量更新部１０４は、例えば、ｙ（ｋ）＝φ（ｋ）Τθにより制御量予測値ｙ（ｋ）を計算する。ここで、Τは転置を表す。
【００４８】
ステップＳ１３：次に、操作量更新部１０４は、制御量予測値ｙ（ｋ）と状態ベクトルφ（ｋ）とに基づいて、次のインデックス（つまり、ｋ＋１）における状態ベクトルφ（ｋ＋１）を更新する。操作量更新部１０４は、例えば、以下により状態ベクトルφ（ｋ＋１）を更新する。
【００４９】
【数５】
0006927446.tif 000006
ここで、ｙ（ｋ）には上記のステップＳ１２で計算した制御量予測値を設定し、ｕ（ｋ）には１を設定する。また、ｙ（ｋ−１），・・・，ｙ（ｋ−Ｎ＋１）及びｕ（ｋ−１），・・・，ｕ（ｋ−Ｍ＋１）には、状態ベクトルφ（ｋ）と同じ値を設定する（すなわち、状態ベクトルφ（ｋ）に含まれるｙ（ｋ−１），・・・，ｙ（ｋ−Ｎ＋１）及びｕ（ｋ−１），・・・，ｕ（ｋ−Ｍ＋１）をそれぞれ設定する。）。
【００５０】
ステップＳ１４：次に、操作量更新部１０４は、時刻τをτ＋ΔＴに更新すると共に、インデックスｋをｋ＋１に更新する。ここで、ΔＴはＡＲＭＡモデルの１ステップの時間幅を表す。ΔＴは制御対象プラント２０の時定数に応じて任意に設定することができるが、例えば、制御周期Ｔｃと同一にすることが好適である。
【００５１】
ステップＳ１５：次に、操作量更新部１０４は、τ≧ｔであるか否かを判定する。そして、τ≧ｔであると判定されなかった場合（ステップＳ１５でＮＯ）、操作量更新部１０４は、ステップＳ１２に戻る。これにより、τ≧ｔとなるまで、ステップＳ１２〜ステップＳ１４が繰り返し実行される。
【００５２】
一方で、τ≧ｔであると判定された場合（ステップＳ１５でＹＥＳ）、操作量更新部１０４は、処理を終了する。これにより、最終的に計算されたｙ（ｋ）が単位ステップ応答Ｓθ（ｔ）として得られる（つまり、Ｓθ（ｔ）＝ｙ（ｋ）が、ステップ応答関数｛Ｓθ（ｔ）｝の時刻ｔにおける単位ステップ応答として計算される。）。
【００５３】
＜モデルパラメータθの推定＞
次に、プラント応答関数｛Ｓθ（ｔ）｝のモデルパラメータθを推定する処理（つまり、モデルパラメータ推定値θｅｓｔを算出する処理）について、図４を参照しながら説明する。図４は、モデルパラメータθの推定処理の一例を説明するためのフローチャートである。なお、図４では、制御量現在値ｙ０と操作量現在値ｕ０とが入力されたものとして、或るインデックスｋにおけるモデルパラメータ推定値θｅｓｔを算出する場合について説明する。なお、このインデックスｋは、図３でプラント応答関数の計算に用いたインデックスｋとは独立した値であり、モデルパラメータ推定処理の実行時インデックスを表す。
【００５４】
ステップＳ２１：モデルパラメータ推定部１０１は、モデルパラメータ推定値θｅｓｔ（ｋ）と共分散行列Ｐ（ｋ）とを初期化するか否かを判定する。ここで、初期化すると判定される場合としては、例えば、モデルパラメータ推定値θｅｓｔを初めて算出する場合（つまり、モデルパラメータ推定値θｅｓｔの初回算出時）、ユーザ等により初期化指示が行われた場合等が挙げられる。
【００５５】
モデルパラメータ推定値θｅｓｔ（ｋ）と共分散行列Ｐ（ｋ）とを初期化すると判定した場合（ステップＳ２１でＹＥＳ）、モデルパラメータ推定部１０１は、ステップＳ２２に進む。一方で、モデルパラメータ推定値θｅｓｔ（ｋ）と共分散行列Ｐ（ｋ）とを初期化すると判定しなかった場合（ステップＳ２１でＮＯ）、モデルパラメータ推定部１０１は、ステップＳ２３に進む。
【００５６】
ステップＳ２２：上記のステップＳ２１でモデルパラメータ推定値θｅｓｔ（ｋ）と共分散行列Ｐ（ｋ）とを初期化すると判定された場合、モデルパラメータ推定部１０１は、θｅｓｔ（０）＝θ０及びＰ（０）＝Ｉと初期化し、更にモデル推定処理の実行時インデックスｋをｋ＝１と初期化する。ここで、θ０としては、既にプラント応答関数｛Ｓθ（ｔ）｝に設定されているモデルパラメータ設定値θを用いてもよいし、ユーザ等が事前に想定した初期値を用いてもよいし、全ての要素が０のベクトル等の固定の初期値を用いてもよい。また、Ｉとしては、単位行列としてもよいし、予め決められた任意の行列としてもよい。
【００５７】
ステップＳ２３：モデルパラメータ推定部１０１は、モデルパラメータ推定値θｅｓｔ（ｋ−１）と状態ベクトルφ（ｋ）と制御量現在値ｙ０とに基づいて、予測誤差ε（ｋ）を計算する。なお、モデルパラメータ推定値θｅｓｔ（ｋ−１）は前回（つまり、ｋ−１のとき）推定したモデルパラメータθの推定値である。また、状態ベクトルφ（ｋ）は前回更新した状態ベクトル（つまり、ｋ−１のときに図３のステップＳ１３で得られた状態ベクトル又はｋ−１のときに後述するステップＳ２７で得られた状態ベクトル）である。
【００５８】
モデルパラメータ推定部１０１は、例えば、以下により予測誤差ε（ｋ）を計算する。
【００５９】
ｙ（ｋ）＝ｙ０
ε（ｋ）＝ｙ（ｋ）−φ（ｋ）Τθｅｓｔ（ｋ−１）
ステップＳ２４：次に、モデルパラメータ推定部１０１は、共分散行列Ｐ（ｋ）を更新する。モデルパラメータ推定部１０１は、例えば、以下の数６により共分散行列Ｐ（ｋ）を更新する。
【００６０】
【数６】
0006927446.tif 000007
ここで、Ｐ（ｋ−１）は前回（つまり、ｋ−１のとき）得られた共分散行列である。
【００６１】
ステップＳ２５：次に、モデルパラメータ推定部１０１は、モデルパラメータ推定値θｅｓｔを更新するか否かを判定する。ここで、モデルパラメータ推定値θｅｓｔを更新すると判定される場合としては、例えば、上記のステップＳ２２での初期化時から所定の期間が経過するまでモデルパラメータ推定値θｅｓｔが更新されなかった場合、ユーザ等により更新指示が行われた場合等が挙げられる。
【００６２】
モデルパラメータ推定値θｅｓｔを更新すると判定した場合（ステップＳ２５でＹＥＳ）、モデルパラメータ推定部１０１は、ステップＳ２６に進む。一方で、モデルパラメータ推定値θｅｓｔを更新すると判定しなかった場合（ステップＳ２５でＮＯ）、モデルパラメータ推定部１０１は、ステップＳ２６を実行せずに、ステップＳ２７に進む。
【００６３】
ステップＳ２６：モデルパラメータ推定部１０１は、モデルパラメータ推定値θｅｓｔ（ｋ）を更新する。モデルパラメータ推定部１０１は、例えば、以下の数７によりモデルパラメータ推定値θｅｓｔ（ｋ）を更新する。
【００６４】
【数７】
0006927446.tif 000008
ステップＳ２７：モデルパラメータ推定部１０１は、制御量現在値ｙ０と操作量現在値ｕ０と状態ベクトルφ（ｋ）とに基づいて、次のインデックス（つまり、ｋ＋１）における状態ベクトルφ（ｋ＋１）を更新する。モデルパラメータ推定部１０１は、例えば、以下の数８により状態ベクトルφ（ｋ＋１）を更新する。
【００６５】
【数８】
0006927446.tif 000009
ここで、ｙ（ｋ）には制御量現在値ｙ０を設定し、ｕ（ｋ）には操作量現在値ｕ０を設定する。また、ｙ（ｋ−１），・・・，ｙ（ｋ−Ｎ＋１）及びｕ（ｋ−１），・・・，ｕ（ｋ−Ｍ＋１）には、状態ベクトルφ（ｋ）と同じ値を設定する（すなわち、状態ベクトルφ（ｋ）に含まれるｙ（ｋ−１），・・・，ｙ（ｋ−Ｎ＋１）及びｕ（ｋ−１），・・・，ｕ（ｋ−Ｍ＋１）をそれぞれ設定する。）。このように、共分散行列の更新とモデルパラメータ推定値θｅｓｔの算出とは逐次最小２乗法に基づく処理で実現される。ただし、逐次最小２乗法は一例であって、例えば、カルマンフィルタ法が用いられてもよい。
【００６６】
ステップＳ２８：モデルパラメータ推定部１０１は、モデルパラメータ推定処理の実行時インデックスｋをｋ＋１に更新する。
【００６７】
なお、例えば、過去Ｌ点の外乱ｖを考慮する場合には、状態ベクトルφ（ｋ）を、ｙ（ｋ）及びｕ（ｋ）に加えてＬ個の外乱ｖ（ｋ）を要素として持つ状態ベクトルに拡大し、外乱に対する共分散行列とモデルパラメータの推定とを行うようにすればよい。
【００６８】
＜先読み応答補正部１１２の動作＞
次に、先読み応答補正部１１２の動作について、図５を参照しながら説明する。図５は、先読み応答補正部１１２の動作の一例を説明するための図である。
【００６９】
図５に示すように、先読み応答補正部１１２は、プラント応答関数｛Ｓθ（ｔ）｝と、目標偏差ｅ０（ｔ）と、操作変化量時系列｛ｄｕ（ｔ）｝とが入力されると、制御量が過去の操作変化量によって現在時刻ｔからＴｐ経過後に変化すると予測される値を先読み応答補正値ｙｎ（ｔ）として算出する。なお、先読み長Ｔｐは制御パラメータ計算部１１１により計算される。
【００７０】
そして、先読み応答補正部１１２は、先読み応答補正値ｙｎ（ｔ）により目標偏差ｅ０（ｔ）を補正した補正目標偏差ｅ＊（ｔ）を出力する。ここで、補正目標偏差ｅ＊（ｔ）は、ｅ＊（ｔ）＝ｒ（ｔ）−（ｙ０（ｔ）＋ｙｎ（ｔ））＝ｅ０（ｔ）−ｙｎ（ｔ）で計算される。なお、先読み応答補正値ｙｎ（ｔ）は、例えば、上記の特許文献２や特許文献３に記載されている手法により算出することができる。
【００７１】
＜操作変化量計算部１１３の動作＞
次に、操作変化量計算部１１３の動作について、図６を参照しながら説明する。図６は、操作変化量計算部１１３の動作の一例を説明するための図である。
【００７２】
図６に示すように、操作変化量計算部１１３は、補正目標偏差ｅ＊（ｔ）が入力されると、この補正目標偏差ｅ＊（ｔ）に対して、制御パラメータ計算部１１１で計算された制御ゲインｋＩを乗じて操作変化量ｄｕ（ｔ）を出力する。すなわち、操作変化量ｄｕ（ｔ）は、ｄｕ（ｔ）＝ｋＩ×ｅ＊（ｔ）で算出される。
【００７３】
ただし、補正目標偏差ｅ＊（ｔ）に対して制御ゲインｋＩを乗じた結果が上限値ｄｕｍａｘを超える場合、操作変化量計算部１１３は、ｄｕｍａｘを操作変化量ｄｕ（ｔ）する。同様に、補正目標偏差ｅ＊（ｔ）に対して制御ゲインｋＩを乗じた結果が下限値ｄｕｍｉｎを下回る場合、操作変化量計算部１１３は、ｄｕｍｉｎを操作変化量ｄｕ（ｔ）する。これにより、上下限範囲に対するリミッターを設けることができる。
【００７４】
＜制御パラメータの計算＞
次に、制御パラメータとして制御ゲインｋＩと先読み長Ｔｐを計算する処理について、図７を参照しながら説明する。図７は、制御パラメータの計算処理の一例を説明するためのフローチャートである。なお、図７では、プラント応答関数｛Ｓθ（ｔ）｝が入力されたものとする。
【００７５】
ステップＳ３１：制御パラメータ計算部１１１は、プラント応答ゲインｇＰを計算する。プラント応答ゲインｇＰは、プラント応答関数Ｓθ（ｔ）に対して所定の時間の最大値Ｔｍａｘを代入した値として計算する。すなわち、プラント応答ゲインｇＰは、ｇＰ＝Ｓθ（Ｔｍａｘ）により計算される。なお、Ｔｍａｘは、例えば、先読み応答補正部１１２で予測可能な期間の最大値（言い換えると、先読み可能な期間の最終時刻）とすること等が考えられる。
【００７６】
ステップＳ３２：次に、制御パラメータ計算部１１１は、予め設定した所定の係数α（ただし、０＜α≦１）とプラント応答ゲインｇＰとに基づいて、制御ゲインｋＩを計算する。制御ゲインｋＩは、係数αに対して２の平方根を乗じると共にプラント応答ゲインｇＰで割った値により計算される。すなわち、制御ゲインｋＩは、以下の数９により計算される。
【００７７】
【数９】
0006927446.tif 000010
ここで、上記の数９に示す制御ゲインｋＩは、制御対象プラント２０の伝達関数Ｐ（ｓ）が以下の数１０に示すような１次遅れで近似できるときの考察から得られた。
【００７８】
【数１０】
0006927446.tif 000011
なお、Ｔは時定数である。
【００７９】
このとき、制御器の伝達関数Ｃ（ｓ）を以下の数１１に示す式で近似する。
【００８０】
【数１１】
0006927446.tif 000012
すると、一巡伝達関数Ｌ（ｓ）は、以下の数１２に示す式となる。
【００８１】
【数１２】
0006927446.tif 000013
したがって、周波数ω＝１／Ｔにおいて、
【００８２】
【数１３】
0006927446.tif 000014
となるため、αＴが１以下であれば、位相余裕が４５度以上となる。
【００８３】
一方で、αＴが１以上であっても、必ず
【００８４】
【数１４】
0006927446.tif 000015
であるため、ゲイン余裕は無限大であり、閉ループは不安定にはならない。
【００８５】
なお、係数αは０＜α≦１に限られないが、０＜α≦１とすることが好適である。また、制御対象プラント２０が１次遅れよりも高次数で、位相余裕が小さくなることが想定される場合は、係数αで予め調整しておくことができる。
【００８６】
ステップＳ３３：次に、制御パラメータ計算部１１１は、制御対象プラント２０のプラント時定数Ｔｅｓｔを計算する。プラント時定数Ｔｅｓｔは、プラント応答関数Ｓθ（ｔ）が予め設定した所定の値となる時刻ｔとして計算する。すなわち、プラント時定数Ｔｅｓｔは、以下の数１５により計算される。
【００８７】
【数１５】
0006927446.tif 000016
つまり、Ｓθ（ｔ）＝０．６３２ｇＰとなるｔをＴｅｓｔとする。
【００８８】
なお、Ｔｅｓｔを計算する方法は様々に考えられるが、例えば、Ｓθ（ｔ）に対してｔ＝１，２，３，・・・と候補の値を順に代入し、初めて０．６３２ｇＰを超えたときのｔをＴｅｓｔとすること等が考えらえる。
【００８９】
ステップＳ３４：最後に、制御パラメータ計算部１１１は、予め設定した所定の係数βとプラント時定数Ｔｅｓｔとに基づいて、先読み長Ｔｐを計算する。先読み長Ｔｐは、所定の係数βに対してプラント時定数Ｔｅｓｔを乗じた値として計算される。すなわち、先読み長Ｔｐは、Ｔｐ＝βＴｅｓｔにより計算される。なお、係数βは任意の値とすることが可能であるが、０＜β≦１とすることが好適である。
【００９０】
［第二の実施形態］
以下、第二の実施形態について説明する。なお、第二の実施形態では、第一の実施形態との相違点について説明し、第一の実施形態と同一の構成要素についてはその説明を省略する。
【００９１】
＜制御装置１０の全体構成＞
まず、本実施形態に係る制御装置１０の全体構成について、図８を参照しながら説明する。図８は、第二の実施形態に係る制御装置１０の全体構成の一例を示す図である。
【００９２】
図８に示すように、本実施形態に係る制御装置１０は、第一の実施形態で説明した各部に加えて、別置制御器１０６と、切替器１０７とを有する。これらは、例えば、制御装置１０にインストールされた１以上のプログラムが、ＣＰＵ等のプロセッサに実行させる処理により実現される。なお、第二の実施形態では、操作量更新部１０４から出力される操作量をｕａとする。
【００９３】
別置制御器１０６は、操作量更新部１０４とは別の制御器であり、例えば、ＰＩＤ制御やＭＰＣ等の既知の技術に基づく制御器である。別置制御器１０６は、目標偏差ｅ０を入力して、操作量ｕｂを出力する。
【００９４】
切替器１０７は、操作量更新部１０４から出力された操作量ｕａと、別置制御器１０６から出力された操作量ｕｂとのいずれを操作量ｕとして制御対象プラント２０に出力するかを切り替える。この切り替えは、例えば、ユーザの操作等によって行うことができる。
【００９５】
なお、図８に示す例では制御装置１０が別置制御器１０６を有しているが、これに限られず、例えば、別置制御器１０６は、制御装置１０の外部に存在してもよい。
【００９６】
これにより、本実施形態に係る制御装置１０は、例えば、初期段階では既設のＰＩＤ制御等を別置制御器１０６として従来の制御方式で運用しつつ、プラント応答モデルが十分に更新されたタイミングでモデル予測制御に切り替える、といったことが可能になる。このため、本実施形態に係る制御装置１０では、従来の運用方法に干渉することなく、自然にモデル予測制御に切り替えることが可能となり、プラント運用者の負担を低減することができる。
【００９７】
＜制御装置１０のハードウェア構成＞
次に、第一及び第二の実施形態に係る制御装置１０のハードウェア構成について、図９を参照しながら説明する。図９は、一実施形態に係る制御装置１０のハードウェア構成の一例を示す図である。
【００９８】
図９に示すように、第一及び第二の実施形態に係る制御装置１０は、入力装置２０１と、表示装置２０２と、外部Ｉ／Ｆ２０３と、通信Ｉ／Ｆ２０４と、プロセッサ２０５と、メモリ装置２０６とを有する。これらの各ハードウェアは、それぞれがバス２０７により通信可能に接続される。
【００９９】
入力装置２０１は、例えば、タッチパネルや各種ボタン等である。表示装置２０２は、例えば、表示パネル等である。なお、制御装置１０は、入力装置２０１及び表示装置２０２のうちの少なくとも一方を有していなくてもよい。
【０１００】
外部Ｉ／Ｆ２０３は、記録媒体２０３ａ等の外部装置とのインタフェースである。記録媒体２０３ａとしては、例えば、ＳＤメモリカード（Secure Digital memory card）やＵＳＢ（Universal Serial Bus）メモリカード等が挙げられる。
【０１０１】
通信Ｉ／Ｆ２０４は、制御装置１０を通信ネットワークに接続するためのインタフェースである。プロセッサ２０５は、例えば、ＣＰＵ等の各種演算装置である。メモリ装置２０６は、例えば、ＳＳＤ（Solid State Drive）、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）、フラッシュメモリ等の各種記憶装置である。
【０１０２】
なお、図９に示すハードウェア構成は一例であって、制御装置１０は、他のハードウェア構成を有していてもよい。
【０１０３】
［実施例］
次に、実施例について説明する。本実施例では、第二の実施形態に係る制御装置１０によって制御対象プラント２０を制御する場合について説明する。
【０１０４】
まず、真のプラント応答は、ΔＴを時間メッシュとして、以下で表されるものとする。
【０１０５】
ｙ（ｔ）＝ａ１ｙ（ｔ−ΔＴ）＋ａ２ｙ（ｎ−２ΔＴ）＋ｂ１ｕ（ｎ−ΔＴ）＋ｂ２ｕ（ｎ−２ΔＴ）＋δ（ｔ）
ここで、δ（ｔ）はノイズを表す。
【０１０６】
ΔＴ＝Ｔｃとすれば、プラント応答関数｛Ｓθ（ｔ）｝は２次のＡＲＭＡモデルとなり、ｙ（ｋ）＝ａ１ｙ（ｋ−１）＋ａ２ｙ（ｋ−２）＋ｂ１ｕ（ｋ−１）＋ｂ２ｕ（ｋ−２）と表せる。
【０１０７】
このとき、状態ベクトルは、
【０１０８】
【数１６】
0006927446.tif 000017
と表せる。
【０１０９】
モデルパラメータ推定値は、
【０１１０】
【数１７】
0006927446.tif 000018
と表せる。
【０１１１】
モデルパラメータの真値は、
【０１１２】
【数１８】
0006927446.tif 000019
と表せる。
【０１１３】
また、制御量予測値は、ｙ（ｋ）＝φ（ｋ）Τθ（ｋ−１）により計算される。
【０１１４】
本実施例における制御対象プラント２０の応答を図１０に示す。図１０に示すように、本実施例における制御対象プラント２０は、操作量ｕとしてステップ入力を入力すると、制御量ｙは滑らかに増加し、およそ時刻２００付近で定常値になる。なお、外乱ｖは０で一定である。
【０１１５】
次に、本実施例における制御装置１０の動作チャートを図１１に示す。本実施例では、制御対象プラント２０をバッチ運転し、３つのバッチ（バッチ１、バッチ２、バッチ３）が連続するものとする。本実施例では、バッチ１を別置制御器１０６から出力された操作量で制御し、バッチ２及びバッチ３を操作量更新部１０４から出力された操作量で制御するものとする。また、別置制御器１０６にはＰＩ制御を用いる。なお、図中の「提案手法」が操作量更新部１０４により制御する手法のことを表す。
【０１１６】
このとき、運転開始前の時刻ｔ０で共分散行列とモデルパラメータの初期化を行う。
【０１１７】
モデルパラメータの初期値は、
【０１１８】
【数１９】
0006927446.tif 000020
とした。また、共分散行列の初期値は、
【０１１９】
【数２０】
0006927446.tif 000021
とした。
【０１２０】
共分散行列とモデルパラメータの初期化後は、モデルパラメータ推定部１０１を動作させて共分散行列とモデルパラメータの逐次更新を行う。
【０１２１】
時刻ｔ０でバッチ１を開始し、別置制御器１０６により更新される操作量で制御対象プラント２０を制御する。
【０１２２】
時刻ｔ１でバッチ１が終了し、バッチ１で推定されたモデルパラメータを用いて、制御パラメータ計算部１１１で制御ゲインｋＩと先読み長Ｔｐを計算する。そして、そのモデルパラメータと制御ゲインｋＩと先読み長Ｔｐとを用いて、操作量更新部１０４による制御に切り替える。
【０１２３】
時刻ｔ１でバッチ２を開始し、操作量更新部１０４により更新される操作量で制御対象プラント２０を制御する。
【０１２４】
時刻ｔ２でバッチ２が終了し、バッチ２で推定されたモデルパラメータを用いて、制御パラメータ計算部１１１で再度制御ゲインｋＩと先読み長Ｔｐを計算する。そして、そのモデルパラメータと制御ゲインｋＩと先読み長Ｔｐとを用いて、操作量更新部１０４による制御を継続する。
【０１２５】
時刻ｔ２でバッチ３を開始し、操作量更新部１０４により更新される操作量で制御対象プラント２０を制御する。そして、時刻ｔ３でバッチ３が終了する。
【０１２６】
ここで、上述したように、本実施例における別置制御器１０６はＰＩ制御であり、以下により操作量ｕを計算する。
【０１２７】
【数２１】
0006927446.tif 000022
なお、ｋｐ＝０．０５，ＴＩ＝１５０とした。
【０１２８】
一方で、操作量更新部１０４は第一の実施形態で説明した通りであるが、係数はα＝１，β＝１とした。
【０１２９】
このとき、制御応答とモデルパラメータ推定値を図１２に示す。なお、上段の左図はバッチ１での制御応答、上段の真中の図はバッチ２での制御応答、上段の右図はバッチ３での制御応答である。また、下段の左図はバッチ１で推定されたモデルパラメータ推定値、下段の真中の図はバッチ２で推定されたモデルパラメータ推定値である。
【０１３０】
図１２に示すように、バッチ１では、目標値ｒに対して制御量ｙがオーバーシュート、アンダーシュートしており、応答性が良くない。一方で、パラメータ推定値は、制御応答に従って時刻０からパラメータ推定値が変化し、およそ時刻２００で一定の値に収束している。
【０１３１】
バッチ２では、バッチ１で推定したモデルパラメータを用いて、操作量更新部１０４による制御に切り替える。なお、このモデルパラメータは、バッチ２が開始する前に、プラント応答関数（ステップ応答関数）に設定される。
【０１３２】
図１２に示すように、バッチ２では、バッチ１よりも明らかに制御応答が改善しており、制御量がより素早く目標値に追従し、振動も抑制されている。また、パラメータ推定値は、バッチ１で推定したモデルパラメータからスタートし、途中でやや変化するものの、およそ時刻２００で一定の値に収束している。
【０１３３】
バッチ３では、バッチ２で推定したモデルパラメータを用いて、操作量更新部１０４による制御を継続する。図１２に示すように、バッチ３の制御応答は、バッチ２とほぼ同等の制御性能となっている。
【０１３４】
また、制御パラメータ計算部１１１は、バッチ１終了時とバッチ２終了時のそれぞれで制御パラメータを計算した。バッチ１終了時では、このときのモデルパラメータ推定値θを設定したプラント応答関数Ｓθ（ｔ）によりプラント応答ゲインｇＰとプラント時定数Ｔｅｓｔを計算し、制御ゲインｋＩと先読み長Ｔｐを計算した。同様に、バッチ２終了時でも、このときのモデルパラメータ推定値を設定したプラント応答関数によりプラント応答ゲインｇＰとプラント時定数Ｔｅｓｔを計算し、制御ゲインｋＩと先読み長Ｔｐを計算した。なお、プラント時定数Ｔｅｓｔを計算する際には、ｔ＝１０，２０，３０，・・・と１０刻みの値をプラント応答関数Ｓθ（ｔ）に代入し、０．６３２ｇＰを初めて超えたときのｔの値をＴｅｓｔとした。
【０１３５】
バッチ１〜バッチ３における制御性能の変化を図１３に示す。図１３に示すように、バッチ２とバッチ３では、バッチ１に比べて操作量の立ち上がりが早く、かつ、操作量の収束も早い。また、バッチ３はバッチ２と比較して、オーバーシュート量が改善している。したがって、提案手法によってより高い制御性能が実現できることがわかる。
【０１３６】
なお、本実施例ではバッチ毎にモデルパラメータ推定値の設定と制御ゲイン及び先読み長の計算とを行ったが、これに限られず、例えば、モデルパラメータ推定値の設定と制御ゲイン及び先読み長の計算とを都度行ってもよい。
【０１３７】
また、運用開始前の事前検討等によりモデルパラメータが推定できる場合には、そのモデルパラメータを初期値として設定することが好適である。
【０１３８】
また、第一及び第二の実施形態に係る制御装置１０で推定又は計算したモデルパラメータや制御パラメータ（制御ゲイン及び先読み長）は、外部に取り出し可能であってもよい。
【０１３９】
また、第二の実施形態に係る制御装置１０の切替器１０７は、例えば、モデルパラメータ推定値の収束傾向に基づいて自動的に切り替わるようにしてもよい。例えば、初期段階では別置制御器１０６の操作量を制御対象プラント２０に出力するようにしておき、モデルパラメータの更新による変化量が一定値以下になった場合に、操作量更新部１０４の操作量を制御対象プラント２０に出力するように切り替えてもよい。
【０１４０】
本発明は、具体的に開示された上記の実施形態に限定されるものではなく、特許請求の範囲から逸脱することなく、種々の変形や変更、既知の技術との組み合わせ等が可能である。
【符号の説明】
【０１４１】
１０        制御装置
２０        制御対象プラント
１０１      モデルパラメータ推定部
１０２      計測部
１０３      差分器
１０４      操作量更新部
１０５      タイマ
１０６      別置制御器
１０７      切替器
１１１      制御パラメータ計算部
１１２      先読み応答補正部
１１３      操作変化量計算部
１１４      加算器
２０１      入力装置
２０２      表示装置
２０３      外部Ｉ／Ｆ
２０３ａ    記録媒体
２０４      通信Ｉ／Ｆ
２０５      プロセッサ
２０６      メモリ装置
２０７      バス
【要約】      （修正有）【課題】エッジデバイス上でモデル予測制御を実現しつつ、モデル予測制御に必要な制御パラメータも自動調整する装置を提供すること。【解決手段】制御装置は、制御対象に対する操作量を出力し、制御対象の制御量を目標値に追従させる制御装置であって、制御量と、操作量とに基づいて、制御対象のプラント応答を表すプラント応答関数が持つモデルパラメータの推定値を算出するモデルパラメータ推定部と、プラント応答関数に基づいて、制御ゲインと先読み長とを含む制御パラメータを算出する制御パラメータ算出部と、目標偏差を算出する目標偏差算出部と、先読み長経過後における制御量の予測値により目標偏差を補正した補正目標偏差を算出する補正目標偏差算出部と、補正目標偏差に対して制御ゲインを乗じて操作量の変化量を算出する操作変化量算出部と、操作量の変化量と、現在の操作量とを加算して新たな操作量を算出する操作量算出部と、を有する。【選択図】図１【図１】
0006927446.tif 000023
【図２】
0006927446.tif 000024
【図３】
0006927446.tif 000025
【図４】
0006927446.tif 000026
【図５】
0006927446.tif 000027
【図６】
0006927446.tif 000028
【図７】
0006927446.tif 000029
【図８】
0006927446.tif 000030
【図９】
0006927446.tif 000031
【図１０】
0006927446.tif 000032
【図１１】
0006927446.tif 000033
【図１２】
0006927446.tif 000034
【図１３】
0006927446.tif 000035