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

Abstract

【課題】エッジデバイス上でモデル予測制御を実現しつつ、モデル予測制御に必要な制御パラメータも自動調整すること。【解決手段】制御対象の制御量を目標値に追従させる制御装置であって、プラント応答関数のモデルパラメータの推定値を算出するモデルパラメータ推定部と、前記プラント応答関数と、現在に至るまでの過去の操作量の変化量とに基づいて、先読み長経過後における制御量の予測値により目標偏差を補正した補正目標偏差を算出する補正目標偏差算出部と、操作量の変化量と、現在の操作量とを加算して新たな操作量を算出する操作量算出部と、を有し、操作量算出部による前記新たな操作量の算出とを所定の制御周期毎に逐次的に繰り返すと共に、前記モデルパラメータ推定部による前記モデルパラメータの推定値の算出と、前記制御パラメータ算出部による前記制御パラメータの算出とを所定の調整周期毎に逐次的に繰り返す。【選択図】図１

Description

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

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

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

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

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

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

更に、セルフチューニング方式によりＰＩＤ制御のパラメータを決定する技術が知られている（例えば、非特許文献４、非特許文献５）。

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

ヤン Ｍ． マチエヨフスキー，「モデル予測制御 制約のもとでの最適制御」，東京電機大学出版局，２００５ 足立 修一，「ＭＡＴＬＡＢによる制御のためのシステム同定」，東京電機大学出版局，１９９６ 足立 修一，「ＭＡＴＬＡＢによる制御のための上級システム同定」，東京電機大学出版局，２００４ 山本 透，藤井 憲三，「実用化に向けたセルフチューニング制御系の設計」，電気学会論文誌Ｄ，１２３巻９号（２００３） 山本 透，「セルフチューニング方式の新展開 一般化最小分散制御とＰＩＤ制御」，計測と制御，４０巻１０号（２００１）

非特許文献４、非特許文献５に記載されている技術は、比較的計算資源が乏しいエッジデバイス（例えば、ＰＬＣ等の制御装置）上でパラメータのチューニングを実現し得る技術として有望であるが、以下のような課題がある。

（１）ＰＩＤ制御では３つのパラメータしか自由度がないため理想的な応答にマッチングできない場合がある。

（２）一般化最小分散制御を用いる場合、予めプラント時定数を十分考慮した参照モデル関数を所与とする必要がある。

（３）ＰＩＤ制御は微分信号を用いるためノイズに弱い。

（４）ＰＩＤ制御ではむだ時間が長いプラントや逆応答を有するプラントの制御を十分に実現することができない。

一方で、モデル予測制御は、プラント応答モデルに比較的自由度が高く、ＰＩＤ制御が苦手とするむだ時間の長いプラントや逆応答を有するプラントにも対応し得る。

しかしながら、モデル予測制御の制御パラメータ（例えば、制御ゲイン等）を決定するためには、モデルに基づくシミュレーションやデータ解析といった事前のオフライン調整作業（この作業はエンジニアリングとも呼ばれる。）が必要である。しかも、このエンジニアリングは、一般に、それに関する技能を有する者が専用のソフトウェアを用いて行う必要がある。

加えて、モデル予測制御はモデルに基づく制御方式であるため制御対象のモデルを予め同定する必要があり、このモデル同定と上記のエンジニアリングとをエッジデバイス上で実行することは困難である。

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

上記目的を達成するため、一実施形態に係る制御装置は、制御対象に対する操作量を出力し、前記制御対象の制御量を目標値に追従させる制御装置であって、前記制御量と、前記操作量とに基づいて、前記制御対象のプラント応答モデルを表す関数であって、かつ、モデルパラメータが含まれる関数であるプラント応答関数の前記モデルパラメータの推定値を算出するモデルパラメータ推定部と、前記プラント応答関数に基づいて、制御ゲインと先読み長とを含む制御パラメータを算出する制御パラメータ算出部と、前記目標値と、現在の制御量との差である目標偏差を算出する目標偏差算出部と、前記プラント応答関数と、現在に至るまでの過去の前記操作量の変化量とに基づいて、前記先読み長経過後における前記制御量の予測値により前記目標偏差を補正した補正目標偏差を算出する補正目標偏差算出部と、前記補正目標偏差に対して前記制御ゲインを乗じて前記操作量の変化量を算出する操作変化量算出部と、前記操作量の変化量と、前記現在の操作量とを加算して新たな操作量を算出する操作量算出部と、を有し、前記目標偏差算出部による前記目標偏差の算出と、前記補正目標偏差算出部による前記補正目標偏差の算出と、前記操作量算出部による前記新たな操作量の算出とを所定の制御周期毎に逐次的に繰り返すと共に、前記モデルパラメータ推定部による前記モデルパラメータの推定値の算出と、前記制御パラメータ算出部による前記制御パラメータの算出とを所定の調整周期毎に逐次的に繰り返す。

エッジデバイス上でモデル予測制御を実現しつつ、モデル予測制御に必要な制御パラメータも自動調整することができる。

本実施形態に係る制御装置の全体構成の一例を示す図である。 プラント応答関数の動作の一例を説明するための図である。 プラント応答関数の計算処理の一例を説明するためのフローチャートである。 モデルパラメータの推定処理の一例を説明するためのフローチャートである。 先読み応答補正部の動作の一例を説明するための図である。 操作変化量計算の動作の一例を説明するための図である。 制御パラメータの計算処理の一例を説明するためのフローチャートである。 先読み長の計算概要の一例を説明するための図である。 先読み長の計算詳細の一例を説明するためのフローチャートである。 リミットサイクル制御器の動作の一例を説明するための図（その１）である。 リミットサイクル制御器の動作の一例を説明するための図（その２）である。 本実施形態に係る制御装置のハードウェア構成の一例を しめす図である。 実施例を説明するための図である。 実施例におけるモデルパラメータ推定結果を説明するための図（その１）である。 実施例における制御パラメータ計算結果を説明するための図（その１）である。 実施例における制御応答及び操作量を説明するための図（その１）である。 実施例におけるモデルパラメータ推定結果を説明するための図（その２）である。 実施例における制御パラメータ計算結果を説明するための図（その２）である。 実施例における制御応答及び操作量を説明するための図（その２）である。

以下、本発明の一実施形態について、図面を参照しながら詳細に説明する。以降では、ＰＩＤ制御が苦手とするむだ時間の長いプラントや逆応答を有するプラントを制御対象として、その目標値が与えられた場合に、モデル予測制御により制御量を目標値に追従させるための操作量を計算すると共に、モデル更新のためのモデルパラメータの推定とモデル予測制御に必要な制御パラメータの計算とを行う制御装置１０について説明する。なお、本実施形態に係る制御装置１０は、例えば、ＰＣ（パーソナルコンピュータ）等と比較して計算資源が乏しいエッジデバイス（ＰＬＣやＤＣＳ等）であるものとする。

本実施形態に係る制御装置１０は、任意の目標値ｒや制御対象プラント２０の状態等を示す制御量ｙ、制御対象プラント２０のプラント応答モデル等に基づいて、当該制御対象プラント２０に対する操作量ｕを計算すると共に、モデルパラメータの推定と制御パラメータの計算とを行う。そして、制御装置１０は、この操作量ｕに応じた制御対象プラント２０の制御量ｙを計測し、目標値ｒや制御量ｙ、プラント応答モデル等に基づいて次の操作量ｕを計算すると共に、新たなモデルパラメータの推定と新たな制御パラメータの計算とを行う。このように、本実施形態に係る制御装置１０は、制御量ｙを目標値ｒに追従させるための操作量ｕの計算と、プラント応答モデルのモデルパラメータの推定と、制御パラメータの計算とをオンライン実行中（つまり、制御対象プラント２０の制御中）に繰り返し実行する。すなわち、本実施形態に係る制御装置１０は、毎回のオンライン最適化計算を必要とせずにモデル予測制御を実現すると共に、それに必要な制御パラメータも自動で調整する。

これにより、本実施形態に係る制御装置１０は、比較的乏しい計算資源でもモデル予測制御による高い制御性能を実現しつつ、制御パラメータの調整に必要な作業を不要とすることができる。また、制御パラメータの調整作業（例えば、オフラインにおけるデータ解析やシミュレーション等のモデル予測制御に関するエンジニアリング作業）が不要となることで、それに要する時間・コストを削減することができる。

また、制御対象プラント２０の運用開始後は、プラント特性の季節変動や経年劣化等の継時的変化に追従したモデルパラメータを推定しつつ制御パラメータも自動で調整されるため、制御性能を安定して維持することができる。

更に、事前に調整された参照モデル関数が不要で、ＰＩＤ制御が苦手とするプラント（むだ時間の長いプラントや逆応答を有するプラント）にも対応し得るモデル予測制御を実現することができる。

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

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

図１に示すように、本実施形態に係る制御装置１０は、モデルパラメータ推定部１０１と、計測部１０２と、差分器１０３と、操作量更新部１０４と、制御パラメータ計算部１０５と、第１のタイマ１０６と、第２のタイマ１０７と、リミットサイクル制御器１０８と、切替器１０９とを有する。これらは、例えば、制御装置１０にインストールされた１以上のプログラムが、ＣＰＵ（Central Processing Unit）やＭＰＵ（Micro-Processing Unit）等のプロセッサに実行させる処理により実現される。

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

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

差分器１０３は、目標値ｒと、制御量現在値ｙ_０との差（偏差）を目標偏差ｅ_０として出力する。時刻ｔにおける目標偏差ｅ_０（ｔ）は、ｅ_０（ｔ）＝ｒ（ｔ）－ｙ_０（ｔ）で計算される。なお、本実施形態では、目標値は一定、すなわち、ｒ（ｔ）＝定数であるとする。

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

先読み応答補正部１１１は、プラント応答関数｛Ｓ_θ（ｔ）｝と、目標偏差ｅ_０（ｔ）と、過去の操作量ｕの変化量ｄｕの時系列データである操作変化量時系列｛ｄｕ（ｔ）｝と、先読み長Ｔ_ｐとに基づいて、目標偏差ｅ_０（ｔ）を補正した補正目標偏差ｅ^＊（ｔ）を計算する。

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

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

モデルパラメータ推定部１０１は、調整周期Ｔ_ｔ毎に、制御量現在値ｙ_０（ｔ）と操作量現在値ｕ_０（ｔ）とを入力し、プラント応答関数｛Ｓ_θ（ｔ）｝のモデルパラメータの推定値θ_ｅｓｔを算出し、出力する。プラント応答関数｛Ｓ_θ（ｔ）｝とはモデルパラメータθを含む関数であり、制御対象プラント２０のプラント応答モデルである。以降では、プラント応答関数｛Ｓ_θ（ｔ）｝に設定されているモデルパラメータθを「モデルパラメータ設定値θ」ともいう。なお、モデルパラメータの推定値θ_ｅｓｔを算出する際に、モデルパラメータの初期値θ_０が入力されてもよい。

制御パラメータ計算部１０５は、調整周期Ｔ_ｔ毎に、プラント応答関数｛Ｓ_θ（ｔ）｝と、調整係数α，βとを入力し、制御パラメータとして制御ゲインｋ_Ｉと先読み長Ｔ_ｐとを計算し、操作変化量計算部１１２と先読み応答補正部１１１とにそれぞれ出力する。

第１のタイマ１０６は、制御周期Ｔ_ｃ毎に、計測部１０２と操作量更新部１０４とを動作させる。制御周期Ｔ_ｃは、制御対象プラント２０を制御する周期のことであり、その値は予め設定される。

第２のタイマ１０７は、調整周期Ｔ_ｔ毎に、モデルパラメータ推定部１０１と制御パラメータ計算部１０５とを動作させる。調整周期Ｔ_ｔは、モデルパラメータの推定と制御パラメータの計算とを行う周期のことであり、その値は予め設定される。ただし、Ｔ_ｔ≧Ｔ_ｃであるものとする。なお、ｎを予め決められた任意の自然数として、Ｔ_ｔ＝ｎＴ_ｃと設定されることが好ましい。調整周期Ｔ_ｔを制御周期Ｔ_ｃよりも長くすることで、モデルパラメータと制御パラメータの計算量を削減することができると共に、制御対象プラント２０の時定数が長い場合であってもプラント応答モデルの同定精度（つまり、モデルパラメータの同定精度）を向上させることが可能となる。

リミットサイクル制御器１０８は、操作量更新部１０４とは別の制御器であり、リミットサイクル法により制御対象プラント２０を制御する制御器である。リミットサイクル制御器１０８は、目標偏差ｅ_０を入力して、操作量ｕを出力する。

切替器１０９は、操作量更新部１０４から出力される操作量と、リミットサイクル制御器１０８から出力される操作量とのいずれの操作量を制御対象プラント２０に出力するかを切り替える。以下では、操作量更新部１０４から出力される操作量をｕ_ａ、リミットサイクル制御器１０８から出力される操作量をｕ_ｂとも表す。このとき、切替器１０９は、操作量ｕ_ａと操作量ｕ_ｂとのいずれかを操作量ｕとして制御対象プラント２０に出力する。操作量ｕ_ａと操作量ｕ_ｂとのいずれを操作量ｕとして制御対象プラント２０に出力するかの切り替えは、例えば、ユーザの操作等によって行われてもよいし、所定の期間経過後に自動的に行われてもよいし、モデルパラメータが十分に更新されたタイミングで自動的に行われてもよい。なお、モデルパラメータが十分に更新されたタイミングとは、例えば、モデルパラメータが収束したタイミング、つまりモデルパラメータ設定値を更新する際における更新前後の変化量が所定の値以下となったタイミング等のことである。これにより、本実施形態に係る制御装置１０は、例えば、初期段階ではリミットサイクル制御器１０８により制御対象プラント２０に対する操作量を出力し、モデルパラメータの更新と制御パラメータの計算とを行い、モデルパラメータが十分に更新されたタイミングでモデル予測制御に切り替える、といったことが可能になる。

なお、本実施形態に係る制御装置１０は、例えば、調整周期Ｔ_ｔ毎に計算されたモデルパラメータや制御パラメータ（制御ゲイン、先読み長）をディスプレイ等の表示装置に逐次表示する表示部を有していてもよいし、これらのモデルパラメータや制御パラメータを通信ネットワークにより接続されるＰＣ等といった端末の表示装置上に逐次表示させる表示制御部を有していてもよい。これにより、制御装置１０のユーザや管理者（以降、「ユーザ等」とも表す。）は、モデルパラメータや制御パラメータのトレンドを知ることが可能となり、プラント応答モデルの同定具合や制御パラメータの調整具合を確認することができるようになる。また、プラント応答モデルに急激な変動が起きる場合にその兆候を発見することも可能となる。

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

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

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

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

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

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

ここで、ｋはインデックスである。ただし、これは一例であって、例えば、過去Ｌ点（Ｌは予め設定された１以上の整数）の外乱ｖを考慮したモデルが用いられてもよい。

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

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

と表される。

なお、例えば、過去Ｌ点の外乱ｖを考慮する場合には、状態ベクトルφ（ｋ）を、ｙ（ｋ）及びｕ（ｋ）に加えてＬ個の外乱ｖ（ｋ）を要素として持つ状態ベクトルに拡大すればよい。

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

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

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

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

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

ステップＳ１４：次に、操作量更新部１０４は、時刻τをτ＋ΔＴに更新すると共に、インデックスｋをｋ＋１に更新する。ここで、ΔＴはＡＲＭＡモデルの１ステップの時間幅を表す。プラント応答関数の時間刻みΔＴは、調整周期Ｔ_ｔと制御周期Ｔ_ｃに応じて変更することができる。例えば、Ｔ_ｔ＝４［ｓｅｃ］、Ｔ_ｃ＝１［ｓｅｃ］であるとする。このとき、モデルパラメータは４［ｓｅｃ］周期のサンプリングに基づくＡＲＭＡモデル（又はＡＲＭＡＸモデル）で更新されるため、制御周期１［ｓｅｔ］で見たときに４倍のサンプリングで計算されることとなる。このため、例えば制御周期単位でｔ＝４０［ｓｅｃ］先を予測したい場合、調整周期では１０サンプル分の計算に相当するため、時間刻みΔＴを４に設定すればよい。このように、Ｔ_ｔ＝ｎＴ_ｃと表現される場合には、ΔＴ＝ｎと設定すればよい。これにより、ｎ倍分の時間補正を行って制御量予測値が計算されることになる。

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

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

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

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

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

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

ステップＳ２３：モデルパラメータ推定部１０１は、状態ベクトルφ（ｋ－１）と操作量現在値ｕ_０と制御量前回値ｙ_－１とに基づいて、状態ベクトルφ（ｋ）を更新する。すなわち、モデルパラメータ推定部１０１は、以下により状態ベクトルφ（ｋ）を更新する。

ここで、ｙ（ｋ－１）にはｙ_－１を設定し、ｕ（ｋ）にはｕ_０を設定する。また、ｙ（ｋ－２），・・・，ｙ（ｋ－Ｎ）及びｕ（ｋ－１），・・・，ｕ（ｋ－Ｍ）には、状態ベクトルφ（ｋ－１）と同じ値を設定する（すなわち、状態ベクトルφ（ｋ－１）に含まれるｙ（ｋ－２），・・・，ｙ（ｋ－Ｎ）及びｕ（ｋ－１），・・・，ｕ（ｋ－Ｍ）をそれぞれ設定する。）。更に、モデルパラメータ推定部１０１は、制御量前回値ｙ_－１←ｙ_０と更新する。

なお、例えば、過去Ｌ点の外乱ｖを考慮する場合には、上述したように、状態ベクトルφ（ｋ）を、ｙ（ｋ）及びｕ（ｋ）に加えてＬ個の外乱ｖ（ｋ）を要素として持つ状態ベクトルに拡大し、外乱に対する共分散行列とモデルパラメータの推定とを行うようにすればよい。

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

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

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

ここで、Ｐ（ｋ－１）は前回（つまり、ｋ－１のとき）得られた共分散行列である。また、λは０＜λ≦１を取る忘却係数であり、予め設定された値である。忘却係数λは過去データを忘却するための係数であり、０に近いほど急激に過去データの影響が減少し、１に近いほど過去データの影響が保持される。

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

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

ステップＳ２７：上記のステップＳ２６でモデルパラメータ推定値θ_ｅｓｔを更新すると判定された場合、モデルパラメータ推定部１０１は、モデルパラメータ推定値θ_ｅｓｔ（ｋ）を更新する。モデルパラメータ推定部１０１は、例えば、以下によりモデルパラメータ推定値θ_ｅｓｔ（ｋ）を更新する。

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

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

図５に示すように、先読み応答補正部１１１は、プラント応答関数｛Ｓ_θ（ｔ）｝と、目標偏差ｅ_０（ｔ）と、操作変化量時系列｛ｄｕ（ｔ）｝とが入力されると、制御量が過去の操作変化量によって現在時刻ｔからＴ_ｐ経過後に変化すると予測される値を先読み応答補正値ｙ_ｎ（ｔ）として算出する。なお、先読み長Ｔ_ｐは制御パラメータ計算部１０５により計算される。

そして、先読み応答補正部１１１は、先読み応答補正値ｙ_ｎ（ｔ）により目標偏差ｅ_０（ｔ）を補正した補正目標偏差ｅ^＊（ｔ）を出力する。ここで、補正目標偏差ｅ^＊（ｔ）は、ｅ^＊（ｔ）＝ｒ（ｔ）－（ｙ_０（ｔ）＋ｙ_ｎ（ｔ））＝ｅ_０（ｔ）－ｙ_ｎ（ｔ）で計算される。なお、先読み応答補正値ｙ_ｎ（ｔ）は、例えば、上記の特許文献２や特許文献３に記載されている手法により算出することができる。

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

図６に示すように、操作変化量計算部１１２は、補正目標偏差ｅ^＊（ｔ）が入力されると、この補正目標偏差ｅ^＊（ｔ）に対して、制御パラメータ計算部１０５で計算された制御ゲインｋ_Ｉを乗じて操作変化量ｄｕ（ｔ）を出力する。すなわち、操作変化量ｄｕ（ｔ）は、ｄｕ（ｔ）＝ｋ_Ｉ×ｅ^＊（ｔ）で算出される。

ただし、補正目標偏差ｅ^＊（ｔ）に対して制御ゲインｋ_Ｉを乗じた結果が上限値ｄｕ_ｍａｘを超える場合、操作変化量計算部１１２は、ｄｕ_ｍａｘを操作変化量ｄｕ（ｔ）とする。同様に、補正目標偏差ｅ^＊（ｔ）に対して制御ゲインｋ_Ｉを乗じた結果が下限値ｄｕ_ｍｉｎを下回る場合、操作変化量計算部１１２は、ｄｕ_ｍｉｎを操作変化量ｄｕ（ｔ）とする。これにより、上下限範囲に対するリミッターを設けることができる。なお、操作量現在値ｕ_０と操作変化量ｄｕとが加算器１１３によって加算された後の操作量ｕ_ａが所定の上下限範囲を逸脱しないように、ｄｕ_ｍａｘ及びｄｕ_ｍｉｎを都度設定してもよい。

＜制御パラメータの計算＞
次に、制御パラメータとして制御ゲインｋ_Ｉと先読み長Ｔ_ｐを計算する処理について、図７を参照しながら説明する。図７は、制御パラメータの計算処理の一例を説明するためのフローチャートである。なお、図７では、プラント応答関数｛Ｓ_θ（ｔ）｝と調整係数α，βとが入力されたものとする。

ステップＳ３１：制御パラメータ計算部１０５は、プラント応答関数｛Ｓ_θ（ｔ）｝と調整係数βとに基づいて、先読み長Ｔ_ｐを計算する。制御パラメータ計算部１０５は、プラント応答関数｛Ｓ_θ（ｔ）｝を用いて、十分長い値として予め設定された最終時刻Ｔ_ｍａｘのプラント応答関数値Ｓ_θ（Ｔ_ｍａｘ）に対して調整係数βを乗じた値にプラント応答関数値が等しくなる時点を先読み長Ｔ_ｐ時点とする。すなわち、制御パラメータ計算部１０５は、以下により先読み長Ｔ_ｐを計算する。

Ｆｉｎｄ Ｔ_ｐ，ｗｈｅｒｅ Ｓ_θ（Ｔ_ｐ）＝β×Ｓ_θ（Ｔ_ｍａｘ）
ここで、βは、０＜β≦１であることが好ましい。なお、閉ループの応答が遅い場合には、調整係数βをより小さい値とすることで応答の早さを調整することができる。

また、制御パラメータ計算部１０５は、先読み長Ｔ_ｐ時点でのゲインｇ_ｐを以下により計算する。

ｇ_ｐ＝Ｓ_θ（Ｔ_ｐ）
このように、先読み長Ｔ_ｐに基づいて制御ゲインｇ_ｐを計算する。これにより、例えば、むだ時間と時定数が長い場合であっても、より早く応答させることが可能となる。

なお、Ｔ_ｍａｘは、プラント応答関数が収束するのに十分な値としてもよいし、メモリ上の制約等から定めた値であってもよい。又は、先読み応答補正部１１１で予測可能な期間の最大値をＴ_ｍａｘとしてもよい。

ステップＳ３２：そして、制御パラメータ計算部１０５は、先読み長Ｔ_ｐ時点でのゲインｇ_ｐと調整係数αとに基づいて、制御ゲインｋ_Ｉを計算する。制御パラメータ計算部１０５は、以下により制御ゲインｋ_Ｉを計算する。

すなわち、先読み長Ｔ_ｐ時点でのゲインｇ_ｐに非負の微小値εを加えた値の逆数に対して調整係数αを乗じた値を制御ゲインｋ_Ｉとする。ここで、αは、０＜α≦１であることが好ましい。また、εはゲインが大きくなりすぎるのを回避するための微小な非負の値である。なお、閉ループの応答でのオーバーシュートが多い場合には、調整係数αをより小さい値とすることでオーバーシュートの多さを調整することができる。

＜先読み長Ｔ_ｐの計算概要＞
次に、上記のステップＳ３１における先読み長Ｔ_ｐの計算概要について、図８を参照しながら説明する。図８は、先読み長Ｔ_ｐの計算概要の一例を説明するための図である。

図８に示すように、プラント応答関数Ｓ_θ（ｔ）の値が時刻０からＴ_ｍａｘまで徐々に変化していくが、或る時刻Ｔ_ｐにおいてＳ_θ（Ｔ_ｐ）＝β×Ｓ_θ（Ｔ_ｍａｘ）が成り立つ。そこで、このような時刻Ｔ_ｐを先読み長とする。これにより、先読み長Ｔ_ｐ時点におけるプラントの挙動が最終的なプラントの挙動と連動するため、無駄な操作量の変化を抑制することが可能となり、オーバーシュートやアンダーシュートの発生を抑制することができる。また、それに加えて、むだ時間や逆応答が生じる応答初期の期間を回避した先読み長を決定することが可能となる。

なお、従来技術として公知のＰＩＤ制御においては、むだ時間の大きさに応じて微分ゲインを大きくする必要があり、ノイズに対して弱いという欠点がある。例えば、Ziegler-Nicholsのステップ応答法による制御則では、ｋ_ｐ＝１．２／（ＲＬ）、Ｔ_Ｉ＝２Ｌ、Ｔ_Ｄ＝０．５Ｌ、Ｒ＝ｇ_ｐ／Ｔ_ｐが用いられており、むだ時間Ｌに比例して微分時間Ｔ_Ｄが長く、微分強度が大きくなるため、結果的にノイズやインパルス状の外乱に脆弱な制御となってしまう。これに対して、本実施形態では、観測値の微分信号を用いないため、ノイズやむだ時間に対してロバストな制御を実現することができる。

＜先読み長Ｔ_ｐの計算詳細＞
次に、上記のステップＳ３１における先読み長Ｔ_ｐの計算詳細について、図９を参照しながら説明する。図９は、先読み長Ｔ_ｐの計算詳細の一例を説明するためのフローチャートである。

ステップＳ４１：まず、制御パラメータ計算部１０５は、先読み長Ｔ_ｐを探索するための変数をＴ_１として、Ｔ_１←０とする。

ステップＳ４２：次に、制御パラメータ計算部１０５は、Ｔ_１≧Ｔ_ｍａｘであるか否かを判定する。

Ｔ_１≧Ｔ_ｍａｘであると判定した場合（ステップＳ４２でＹＥＳ）、制御パラメータ計算部１０５は、ステップＳ４６に進む。一方で、Ｔ_１≧Ｔ_ｍａｘであると判定しなかった場合（ステップＳ４２でＮＯ）、制御パラメータ計算部１０５は、ステップＳ４３に進む。

ステップＳ４３：制御パラメータ計算部１０５は、Ｓ_θ（Ｔ_１）≧β×Ｓ_θ（Ｔ_ｍａｘ）を満たすか否かを判定する。

Ｓ_θ（Ｔ_１）≧β×Ｓ_θ（Ｔ_ｍａｘ）を満たすと判定した場合（ステップＳ４３でＹＥＳ）、制御パラメータ計算部１０５は、ステップＳ４５に進む。一方で、Ｓ_θ（Ｔ_１）≧β×Ｓ_θ（Ｔ_ｍａｘ）を満たすと判定しなかった場合（ステップＳ４３でＮＯ）、制御パラメータ計算部１０５は、ステップＳ４４に進む。

ステップＳ４４：制御パラメータ計算部１０５は、Ｔ_１←Ｔ_１＋ΔＴと更新し、ステップＳ４２に戻る。これにより、Ｔ_１≧Ｔ_ｍａｘ又はＳ_θ（Ｔ_１）≧β×Ｓ_θ（Ｔ_ｍａｘ）のいずれかを満たすまで、ステップＳ４４が繰り返し実行される。

ステップＳ４５：Ｓ_θ（Ｔ_１）≧β×Ｓ_θ（Ｔ_ｍａｘ）を満たすと判定した場合、制御パラメータ計算部１０５は、Ｔ_ｐ←Ｔ_１とする。すなわち、制御パラメータ計算部１０５は、Ｓ_θ（Ｔ_１）≧β×Ｓ_θ（Ｔ_ｍａｘ）を満たすＴ_１を先読み長Ｔ_ｐとする。

ステップＳ４６：Ｔ_１≧Ｔ_ｍａｘであると判定しなかった場合、制御パラメータ計算部１０５は、Ｔ_ｐ←γ×Ｔ_ｍａｘとする。ここで、０＜γ＜１である。すなわち、Ｓ_θ（Ｔ_１）≧β×Ｓ_θ（Ｔ_ｍａｘ）を満たすＴ_１が探索できなかった場合は、Ｔ_ｍａｘに対してγを乗じた値を先読み長Ｔ_ｐとする。

＜リミットサイクル制御器１０８の動作＞
次に、リミットサイクル制御器１０８の動作について、図１０及び図１１を参照しながら説明する。図１０及び図１１は、リミットサイクル制御器１０８の動作の一例を説明するための図（その１）である。

図１０に示すように、リミットサイクル制御器１０８は、目標偏差ｅ_０（ｔ）が入力されると、ｅ_０（ｔ）≧０の場合は所定の値ＭＶＨを操作量ｕ_ｂ（ｔ）として出力し、ｅ_０（ｔ）＜０の場合は所定の値ＭＶＬを操作量ｕ_ｂ（ｔ）として出力する。なお、図１０に示す例では、制御対象プラント２０が正動作する（つまり、操作量の上昇が制御量の上昇の方向に動作する）と仮定し、ＭＶＨ＞ＭＶＬとなっているが、制御対象プラント２０が逆動作する場合には条件判定の符号を反転（つまり、ｅ_０（ｔ）≦０の場合はＭＶＨを操作量ｕ_ｂ（ｔ）として出力し、ｅ_０（ｔ）＞０の場合はＭＶＬを操作量ｕ_ｂ（ｔ）として出力）する必要がある。

図１１に示すように、リミットサイクル制御器１０８を閉ループ制御に用いることで、一定周期の振動（リミットサイクル）を引き起こすことができる。

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

図１２に示すように、本実施形態に係る制御装置１０は、入力装置２０１と、表示装置２０２と、外部Ｉ／Ｆ２０３と、通信Ｉ／Ｆ２０４と、プロセッサ２０５と、メモリ装置２０６とを有する。これらの各ハードウェアは、それぞれがバス２０７により通信可能に接続される。

入力装置２０１は、例えば、タッチパネルや各種ボタン等である。表示装置２０２は、例えば、表示パネル等である。なお、制御装置１０は、入力装置２０１及び表示装置２０２のうちの少なくとも一方を有していなくてもよい。

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

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

なお、図１２に示すハードウェア構成は一例であって、制御装置１０は、他のハードウェア構成を有していてもよい。例えば、制御装置１０は、図示したハードウェア以外にも種々のハードウェアを有していてもよい。

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

本実施形態における制御装置１０と制御対象プラント２０とで構成される制御系全体を図１３に示す。図１３に示すように、金属ブロックにヒーターと熱電対とが取り付けられており、ヒーターはＡＣ電源と制御装置１０（ＰＬＣ）のＤＯモジュールに接続され、熱電対は制御装置１０のＡＩモジュールに接続されている。このとき、制御量は熱電対から検知される金属ブロックの温度、操作量は制御装置１０のＤＯモジュールによるＰＷＭ変調のデューティー比（ＯＮ時間／ＰＷＭ周期の割合［％］）であるものとする。

また、制御対象プラント２０の真のプラント応答は、ΔＴを時間メッシュとして、

と表されるものとする。ここで、δ（ｔ）はノイズを表す。

一方で、プラント応答関数は、モデルパラメータθを持つ以下の２次のＡＲＭＡモデルとする。

よって、状態ベクトルは、

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

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

となる。また、制御量推定値は、ｙ_ｅｓｔ（ｋ）＝φ（ｋ）^Τθ（ｋ－１）で計算される。

第１のタイマ１０６と第２のタイマ１０７には異なる値が設定されており、制御周期Ｔ_ｃは５［ｓｅｃ］、調整周期Ｔ_ｔは２０［ｓｅｃ］とした。すなわち、調整周期Ｔ_ｔは制御周期Ｔ_ｃの４倍であり、制御周期の４回に１回の頻度でモデルパラメータ推定部１０１と制御パラメータ計算部１０５が動作する。よって、プラント応答関数の時間刻みΔＴは４である。このように、調整の頻度を抑制することで、制御装置１０（ＰＬＣ）の計算負荷を低減している。

また、操作量の上限を９０、下限を０とした。調整係数はα＝０．９，β＝０．３とした。忘却係数はλ＝０．９５とした。リミットサイクルのパラメータはＭＶＨ＝６０、ＭＶＬ＝４０とした。

このとき、リミットサイクル制御器１０８を用いない場合の制御動作例について説明する。時刻０から操作量更新部１０４による制御を開始、時刻１０００まで制御したときのモデルパラメータ推定結果を図１４、制御パラメータ計算結果を図１５、制御応答及び操作量を図１６にそれぞれ示す。

図１４に示すように、モデルパラメータθの推定値の初期値は０であるが、制御開始の時刻０から値が変化し、およそ時刻５０付近で最終値に近い値となっている。図１５に示すように、制御パラメータは時刻０での変化が大きく、制御ゲインｋ_Ｉと先読み長Ｔ_ｐが大きな値となるが、すぐに低い値へと変化している。図１６に示すように、制御応答を見ると、時刻０付近で比較的大きな振動が生じているものの、すぐに制御量ｙが目標値ｒへと収束している。また、目標値ｒが５１→４９→５１→４９→５１と変化するたびに、徐々にオーバーシュートが減少し、より追従性が向上している。これは、モデルパラメータ推定部１０１によるモデル更新と制御パラメータ計算部１０５による制御パラメータの計算とによって、制御中に自動的に調整が行われ、制御性能が向上していくためである。

このように、モデルパラメータの初期値や制御パラメータの初期値に事前の情報を用いたり、一般化分散制御のような規範モデル（参照モデル関数）を指定したりしなくとも、本実施形態に係る制御装置１０では、完全な初期状態から制御対象プラント２０を制御しつつ、制御運転中に自動的に制御パラメータをチューニングすると共に制御性能を自動的に向上させることができる。

次に、リミットサイクル制御器１０８を用いた場合の制御動作例について説明する。時刻０から時刻５０までリミットサイクル制御器１０８による制御を行い、その後、時刻５０から操作量更新部１０４による制御を開始し、時刻１０００まで制御したときのモデルパラメータ推定結果を図１７、制御パラメータ計算結果を図１８、制御応答及び操作量を図１９にそれぞれ示す。

図１７に示すように、モデルパラメータθの推定値の初期値は０であるが、制御開始の時刻０から値が変化し、およそ時刻５０付近で最終値に近い値となっている。図１８に示すように、制御パラメータ計算部１０５は、リミットサイクル制御器１０８による制御が終了する時刻５０から制御パラメータを計算している。図１５と比較すると、制御パラメータ計算部１０５で計算された制御パラメータは比較的安定しており、リミットサイクル制御時（つまり、時刻０から時刻５０までの制御時）に推定されたモデルパラメータを用いることで安定化を実現できていることがわかる。図１９に示すように、時刻０付近での制御量ｙと操作量ｕの振れ幅が軽減されており、プラントに与える影響を比較的抑制できていることがわかる。また、リミットサイクル制御器１０８による制御から操作量更新部１０４による制御に切り替えた後の応答では、図１８と同様に、目標値ｒが５１→４９→５１→４９→５１と変化するたびに、徐々にオーバーシュートが減少し、より追従性が向上している。

このように、リミットサイクル制御器１０８を用いることで、モデルパラメータや制御パラメータを比較的早期に安定化させることが可能となり、初期の変動を抑制することができる。なお、例えば、事前の検討によりモデルパラメータが予め推定できる場合には、その値を初期値として設定することもできる。

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

１０ 制御装置
２０ 制御対象プラント
１０１ モデルパラメータ推定部
１０２ 計測部
１０３ 差分器
１０４ 操作量更新部
１０５ 制御パラメータ計算部
１０６ 第１のタイマ
１０７ 第２のタイマ
１０８ リミットサイクル制御器
１０９ 切替器
１１１ 先読み応答補正部
１１２ 操作変化量計算部
１１３ 加算器
２０１ 入力装置
２０２ 表示装置
２０３ 外部Ｉ／Ｆ
２０３ａ 記録媒体
２０４ 通信Ｉ／Ｆ
２０５ プロセッサ
２０６ メモリ装置
２０７ バス

Claims (12)

1. 制御対象に対する操作量を出力し、前記制御対象の制御量を目標値に追従させる制御装置であって、
   前記制御量と、前記操作量とに基づいて、前記制御対象のプラント応答モデルを表す関数であって、かつ、モデルパラメータが含まれる関数であるプラント応答関数の前記モデルパラメータの推定値を算出するモデルパラメータ推定部と、
   前記プラント応答関数に基づいて、制御ゲインと先読み長とを含む制御パラメータを算出する制御パラメータ算出部と、
   前記目標値と、現在の制御量との差である目標偏差を算出する目標偏差算出部と、
   前記プラント応答関数と、現在に至るまでの過去の前記操作量の変化量とに基づいて、前記先読み長経過後における前記制御量の予測値により前記目標偏差を補正した補正目標偏差を算出する補正目標偏差算出部と、
   前記補正目標偏差に対して前記制御ゲインを乗じて前記操作量の変化量を算出する操作変化量算出部と、
   前記操作量の変化量と、前記現在の操作量とを加算して新たな操作量を算出する操作量算出部と、
   を有し、
   前記目標偏差算出部による前記目標偏差の算出と、前記補正目標偏差算出部による前記補正目標偏差の算出と、前記操作量算出部による前記新たな操作量の算出とを所定の制御周期毎に逐次的に繰り返すと共に、前記モデルパラメータ推定部による前記モデルパラメータの推定値の算出と、前記制御パラメータ算出部による前記制御パラメータの算出とを所定の調整周期毎に逐次的に繰り返す、制御装置。
2. 前記モデルパラメータ推定部は、前記調整周期毎に、前記モデルパラメータの推定値により前記プラント応答関数に含まれるモデルパラメータを更新する、請求項１に記載の制御装置。
3. 前記調整周期は、前記制御周期に対して所定の自然数ｎを乗じた値であり、
   前記補正目標偏差算出部は、
   前記自然数ｎ倍分の時間補正を行って、前記制御量の予測値を算出する、請求項１又は２に記載の制御装置。
4. 前記目標偏差に基づいてリミットサイクルにより操作量を出力するリミットサイクル制御部と、
   前記リミットサイクル制御部から出力された操作量と、前記操作量算出部で算出された操作量とのいずれを前記制御対象に出力するかを切り替える切替部と、を有し、
   前記モデルパラメータ推定部は、
   前記制御量と、前記制御対象に出力された操作量に基づいて、前記モデルパラメータの推定値を算出する、請求項１乃至３の何れか一項に記載の制御装置。
5. 前記切替部は、
   所定の初期期間においては前記リミットサイクル制御部から出力された操作量を前記制御対象に出力し、前記初期期間経過後においては前記操作量算出部で算出された操作量を前記制御対象に出力する、請求項４に記載の制御装置。
6. 前記制御パラメータ算出部は、
   前記先読み長を、前記先読み長時点における前記プラント応答関数値と、前記制御対象の最終時点における前記プラント応答関数値に対して所定の調整係数βを乗じた値とが一致する時点までの時間幅として算出する、請求項１乃至５の何れか一項に記載の制御装置。
7. 前記制御パラメータ算出部は、
   前記制御ゲインを、前記先読み長時点における前記プラント応答関数値に所定の非負値を加算した値の逆数に対して所定の調整係数αを乗じた値として算出する、請求項１乃至６の何れか一項に記載の制御装置。
8. 前記プラント応答関数は、ＡＲＭＡモデル又はＡＲＭＡＸモデルであり、
   前記モデルパラメータ推定部は、
   前記ＡＲＭＡモデル又はＡＲＭＡＸモデルの係数を前記モデルパラメータとして前記推定値を算出する、請求項１乃至７の何れか一項に記載の制御装置。
9. 前記モデルパラメータ推定部は、
   忘却要素を含む逐次最小２乗法に基づいて、共分散行列の推定と前記モデルパラメータの推定とを行うことで、前記モデルパラメータの推定値を算出する、請求項１乃至８の何れか一項に記載の制御装置。
10. 前記モデルパラメータ推定部により推定された前記モデルパラメータの推定値と、前記制御パラメータ算出部により算出された前記制御パラメータとの少なくとも一方を前記調整周期毎に所定の表示装置に逐次的に出力する表示制御部を有する、請求項１乃至９の何れか一項に記載の制御装置。
11. 制御対象に対する操作量を出力し、前記制御対象の制御量を目標値に追従させる制御装置が、
    前記制御量と、前記操作量とに基づいて、前記制御対象のプラント応答モデルを表す関数であって、かつ、モデルパラメータが含まれる関数であるプラント応答関数の前記モデルパラメータの推定値を算出するモデルパラメータ推定手順と、
    前記プラント応答関数に基づいて、制御ゲインと先読み長とを含む制御パラメータを算出する制御パラメータ算出手順と、
    前記目標値と、現在の制御量との差である目標偏差を算出する目標偏差算出手順と、
    前記プラント応答関数と、現在に至るまでの過去の前記操作量の変化量とに基づいて、前記先読み長経過後における前記制御量の予測値により前記目標偏差を補正した補正目標偏差を算出する補正目標偏差算出手順と、
    前記補正目標偏差に対して前記制御ゲインを乗じて前記操作量の変化量を算出する操作変化量算出手順と、
    前記操作量の変化量と、前記現在の操作量とを加算して新たな操作量を算出する操作量算出手順と、
    を実行し、
    前記目標偏差算出手順による前記目標偏差の算出と、前記補正目標偏差算出手順による前記補正目標偏差の算出と、前記操作量算出手順による前記新たな操作量の算出とを所定の制御周期毎に逐次的に繰り返すと共に、前記モデルパラメータ推定手順による前記モデルパラメータの推定値の算出と、前記制御パラメータ算出手順による前記制御パラメータの算出とを所定の調整周期毎に逐次的に繰り返す、制御方法。
12. 制御対象に対する操作量を出力し、前記制御対象の制御量を目標値に追従させる制御装置が、
    前記制御量と、前記操作量とに基づいて、前記制御対象のプラント応答モデルを表す関数であって、かつ、モデルパラメータが含まれる関数であるプラント応答関数の前記モデルパラメータの推定値を算出するモデルパラメータ推定手順と、
    前記プラント応答関数に基づいて、制御ゲインと先読み長とを含む制御パラメータを算出する制御パラメータ算出手順と、
    前記目標値と、現在の制御量との差である目標偏差を算出する目標偏差算出手順と、
    前記プラント応答関数と、現在に至るまでの過去の前記操作量の変化量とに基づいて、前記先読み長経過後における前記制御量の予測値により前記目標偏差を補正した補正目標偏差を算出する補正目標偏差算出手順と、
    前記補正目標偏差に対して前記制御ゲインを乗じて前記操作量の変化量を算出する操作変化量算出手順と、
    前記操作量の変化量と、前記現在の操作量とを加算して新たな操作量を算出する操作量算出手順と、
    を実行させ、
    前記目標偏差算出手順による前記目標偏差の算出と、前記補正目標偏差算出手順による前記補正目標偏差の算出と、前記操作量算出手順による前記新たな操作量の算出とを所定の制御周期毎に逐次的に繰り返すと共に、前記モデルパラメータ推定手順による前記モデルパラメータの推定値の算出と、前記制御パラメータ算出手順による前記制御パラメータの算出とを所定の調整周期毎に逐次的に繰り返す、プログラム。

Images