JP6904473B1 - モデル作成支援装置、モデル作成支援方法及びプログラム - Google Patents

Abstract

【課題】モデル同定試験に要する時間を削減すること。【解決手段】一実施形態に係るモデル作成支援装置１０は、制御対象となるプラントのモデルを同定するモデル作成支援装置であって、前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を計算するパラメータ推定部と、前記パラメータの推定値を加速法により加速した加速推定値を計算する加速計算部と、を有することを特徴とする。【選択図】図１

Description

本発明は、モデル作成支援装置、モデル作成支援方法及びプログラムに関する。

温調制御装置やＰＬＣ（Programmable Logic Controller）、ＤＣＳ（Distributed Control System）等の制御装置、パーソナルコンピュータや組み込み制御機器上で実装される制御装置等が産業上広く利用されている。また、それに関連する技術として、制御対象の制御量を目標値に追従させることを目的とする様々な制御方式（例えば、ＰＩＤ（Proportional-Integral-Differential）制御、モデル予測制御、内部モデル制御、ＬＱＧ（Linear-Quadratic-Gaussian）制御、Ｈ２制御、Ｈ∞制御等）が知られている。特に、モデル予測制御は、制御対象の状態空間モデルや将来の時間応答モデルを用いた最適化計算を逐次的に行うことで望ましい応答を得る方式であり、産業界で広く用いられている（例えば、非特許文献１）。

なお、モデル予測制御等で用いる線形予測モデルに関してパラメータ同定を行う手法として、逐次最小２乗法、すなわち、ＲＬＳ（Recursive Least Squares）法やカルマンフィルタ法等が知られている（例えば、非特許文献２及び非特許文献３）。また、収束する数列の収束速度を加速する手法として、エイトケンのΔ２乗加速法が知られている（非特許文献４）。

ヤン Ｍ． マチエヨフスキー，「モデル予測制御 制約のもとでの最適制御」，東京電機大学出版局，２００５ 足立 修一，「ＭＡＴＬＡＢによる制御のためのシステム同定」，東京電機大学出版局，１９９６ 足立 修一，「ＭＡＴＬＡＢによる制御のための上級システム同定」，東京電機大学出版局，２００４ 長田直樹，「収束の加速法」，数理解析研究所講究録，京都大学 (1994), 880: 28-43

ところで、プラントの制御に用いる応答モデルのパラメータを計測データから推定するために、操作量に対して或る形状の信号（例えば、ステップ信号やランプ信号、Ｍ系列信号等）を印加するモデル同定試験が従来から一般的に行われている。

しかしながら、例えば、設備の規模が大きかったり、長時間の熱的・化学的反応を伴ったりするような時定数の長いプラントを対象とする場合、ステップ応答が収束するのに数時間以上を要する場合がある。

また、操作量が複数存在する場合に、操作量毎にステップ応答試験を行うと、（操作変数の数）×（１操作変数あたりのステップ応答試験時間）のモデル同定試験時間が必要となる。

更に、プラントのプロセスに投入する材料や外気温等の条件毎にモデルを分けて作成する必要がある場合、（操作変数の数）×（１操作変数あたりのステップ応答試験時間）×（条件数）のモデル同定試験時間が必要となる。

このように、プラントの応答モデルを得るための同定試験には非常に多くの時間が掛かり、運用上のコストや負荷が大きい。

本発明の一実施形態は、上記の点に鑑みてなされたもので、モデル同定試験に要する時間を削減することを目的とする。

上記目的を達成するため、一実施形態に係るモデル作成支援装置１０は、制御対象となるプラントのモデルを同定するモデル作成支援装置であって、前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を計算するパラメータ推定部と、前記パラメータの推定値を加速法により加速した加速推定値を計算する加速計算部と、を有することを特徴とする。

モデル同定試験に要する時間を削減することができる。

本実施形態に係るモデル作成支援システムの全体構成の一例を示す図である。 本実施形態に係るモデル作成支援装置のハードウェア構成の一例を示す図である。 本実施形態に係るモデルパラメータ推定処理の一例を示すフローチャートである。 本実施形態に係る加速値推定処理の一例を示すフローチャートである。 実施例１における計測データを示す図である。 実施例１におけるモデルパラメータ推定値及び加速推定値を示す図である。 実施例１におけるモデルパラメータ推定値及び加速推定値を同時表示した場合を示す図である。 実施例２における計測データを示す図である。 実施例２におけるモデルパラメータ推定値及び加速推定値を示す図である。 実施例２におけるモデルパラメータ推定値及び加速推定値を同時表示した場合を示す図である。

以下、本発明の一実施形態について説明する。本実施形態では、モデル同定試験（つまり、各種プラント等の制御対象の制御に用いる応答モデルのパラメータを推定するための試験）に要する時間を削減することができるモデル作成支援システム１について説明する。

＜全体構成＞
まず、本実施形態に係るモデル作成支援システム１の全体構成について、図１を参照しながら説明する。図１は、本実施形態に係るモデル作成支援システム１の全体構成の一例を示す図である。

図１に示すように、本実施形態に係るモデル作成支援システム１には、モデル作成支援装置１０と、操作装置２０と、プラント３０と、計測装置４０と、表示装置５０とが含まれる。

モデル作成支援装置１０は、所定の制御周期毎に計測装置４０から計測データ（つまり、制御量ｙ（ｔ）、操作量ｕ（ｔ）及び外乱ｖ（ｔ））を受信して、これらの計測データを用いてプラント３０の応答モデルのパラメータ（以下、「モデルパラメータ」ともいう。）を推定するコンピュータ又はコンピュータシステムである。なお、プラント３０の応答モデルを表す関数を「プラント応答関数」ともいう（つまり、モデルパラメータとはプラント応答関数のパラメータのことを意味する。）。ここで、ｔは時刻ステップを表す。例えば、制御周期をＴ_ｃ、現在の時刻ステップ（つまり、最新の時刻ステップ）をｔ_０とすれば、時刻ステップは、ｔ_０，ｔ_０−Ｔ_ｃ，ｔ_０−２Ｔ_ｃ，・・・と表される。

このとき、モデル作成支援装置１０は、モデルパラメータの推定値（以下、「モデルパラメータ推定値」ともいう。）を加速法により加速させることで、通常よりも早くモデルパラメータを推定する。なお、以下、モデルパラメータ推定値を加速法により加速させた値を「加速推定値」ともいう。この加速推定値をモデルパラメータ値と扱うことで、通常のモデル同定試験よりも早期にモデル同定が可能となり、その試験時間を削減することが可能となる。

操作装置２０は、操作信号に応じて、所定の制御周期毎にプラント３０に対して操作量ｕを出力する。プラント３０は、制御対象の各種設備や工場、発電所、機器等である。計測装置４０は、所定の制御周期毎にプラント３０の制御量ｙ（ｔ）、当該プラント３０に対する操作量ｕ（ｔ）及び当該プラント３０に対する外乱ｖ（ｔ）を計測し、これら制御量ｙ（ｔ）、操作量ｕ（ｔ）及び外乱ｖ（ｔ）を計測データとしてモデル作成支援装置１０に送信するセンサ等である。表示装置５０は、モデル作成支援装置１０で推定されたモデルパラメータ推定や加速推定値を表示する。なお、操作信号は、例えば、プラント３０を制御する制御装置等から送信される。

ここで、図１に示すように、本実施形態に係るモデル作成支援装置１０は、フィルタ部１０１と、モデルパラメータ推定部１０２と、加速計算部１０３とを有する。

フィルタ部１０１は、所定の制御周期毎に計測装置４０から計測データ（制御量ｙ（ｔ）、操作量ｕ（ｔ）及び外乱ｖ（ｔ））を受信し、フィルタパラメータを用いて、これらの制御量ｙ（ｔ）、操作量ｕ（ｔ）及び外乱ｖ（ｔ）の時系列データ（つまり、｛ｙ（ｔ）｝、｛ｕ（ｔ）｝及び｛ｖ（ｔ）｝）に対して所定のフィルタ処理を行う。このようなフィルタ処理を行うことで、計測データがデータ処理に適した形式に変換される。なお、フィルタ処理としては、例えば、移動平均やダウンサンプリング等が挙げられる。フィルタパラメータは、例えば、移動平均の区間やダウンサンプリングの間隔等の種々のパラメータである。以下、制御量ｙ（ｔ）、操作量ｕ（ｔ）及び外乱ｖ（ｔ）はフィルタ処理後の計測データであるものとする。また、時刻ステップｔ＝ｔ_０（つまり、現在の時刻ステップ）における制御量ｙ（ｔ_０）、操作量ｕ（ｔ_０）及び外乱ｖ（ｔ_０）をそれぞれｙ_０、ｕ_０及びｖ_０とも表す。

モデルパラメータ推定部１０２は、制御量ｙ（ｔ）と操作量ｕ（ｔ）と外乱ｖ（ｔ）とを入力として、モデルパラメータ推定値θ_ｅｓｔ（ｋ）を逐次的に算出する。ここで、モデルパラメータ推定部１０２には、モデルパラメータ推定値の初期値θ_０が与えられる。以下、モデルパラメータ推定値の初期値を「モデルパラメータ初期値」ともいう。なお、ｋは初期値として１を取り、モデルパラメータ推定値θ_ｅｓｔ（ｋ）が算出される毎に１ずつ加算されるインデックスである。

加速計算部１０３は、モデルパラメータ推定値θ_ｅｓｔ（ｋ）を入力として、加速法により、モデルパラメータの加速推定値θ_ａｃｃ（ｋ）を逐次的に算出する。

ここで、モデルパラメータ推定値θ_ｅｓｔ（ｋ）及び加速推定値θ_ａｃｃ（ｋ）は表示装置５０に出力され、当該表示装置５０上に表示される。また、必要に応じて、当該表示装置５０上でアラームが発報される。アラームが発報される場合としては、例えば、モデルパラメータ推定値θ_ｅｓｔ（ｋ）と加速推定値θ_ａｃｃ（ｋ）とを比較結果、所定の閾値以上の乖離が存在する場合等である。これは、モデルパラメータ推定値θ_ｅｓｔ（ｋ）と加速推定値θ_ａｃｃ（ｋ）とで乖離が大きい場合は、モデル同定試験のデータが不十分であることを示唆するため、当該試験を不適切に早期に打ち切ってしまうことを防止するためである。これにより、モデル同定試験の担当者等は、その試験時間を削減（短縮）することができると共に、短縮し過ぎない適切な試験時間の長さを知ることができる。なお、モデルパラメータ推定値θ_ｅｓｔ（ｋ）と加速推定値θ_ａｃｃ（ｋ）との間に所定の閾値以上の乖離が存在するか否かの比較判定は、モデル作成支援装置１０で行われてもよいし、表示装置５０が演算装置を備える場合には当該表示装置５０で行われてもよい。

また、上記以外にも、例えば、加速推定値θ_ａｃｃ（ｋ）の変化が一定期間の間、微小である場合（例えば、当該加速推定値θ_ａｃｃ（ｋ）の全ての要素の変化が予め設定された範囲内に収まっている場合等）に当該表示装置５０上にアラームが発報されてもよい。これは、この場合、加速推定値θ_ａｃｃ（ｋ）が収束したと考えられるためである。

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

図２に示すように、本実施形態に係るモデル作成支援装置１０は、外部Ｉ／Ｆ２０１と、通信Ｉ／Ｆ２０２と、ＲＡＭ（Random Access Memory）２０３と、ＲＯＭ（Read Only Memory）２０４と、ＣＰＵ（Central Processing Unit）２０５と、補助記憶装置２０６とを有する。これら各ハードウェアは、バス２０７を介して相互に通信可能に接続されている。

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

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

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

本実施形態に係るモデル作成支援装置１０は、図２に示すハードウェア構成を有することにより、後述する各種処理を実現することができる。なお、図２に示すハードウェア構成は一例であって、モデル作成支援装置１０は、他のハードウェア構成を有していてもよい。例えば、モデル作成支援装置１０は、ＧＰＵ（Graphics Processing Unit）等の演算装置を有していてもよいし、複数のＣＰＵ２０５や複数の補助記憶装置２０６を有していてもよい。

＜モデルパラメータ推定処理＞
次に、モデルパラメータ推定部１０２がモデルパラメータ推定値θ_ｅｓｔ（ｋ）を逐次的に算出する場合について説明する。モデルパラメータ推定部１０２は、ｋ＝１，２，・・・に対してモデルパラメータ推定処理を繰り返し実行することで、モデルパラメータ推定値θ_ｅｓｔ（ｋ）を逐次的に算出する。以下、プラント応答関数としてはＡＲＭＡ（autoregressive and moving average）モデルを想定した上で、或るｋにおけるモデルパラメータ推定処理について、図３を参照しながら説明する。図３は、本実施形態に係るモデルパラメータ推定処理の一例を示すフローチャートである。なお、モデルパラメータ推定部１０２は、例えば、所定の制御周期（又は、予め決められた時間幅）毎にモデルパラメータ推定処理を繰り返し実行すればよい。

まず、モデルパラメータ推定部１０２は、モデルパラメータ推定値θ_ｅｓｔと共分散行列Ｐとを初期化するか否かを判定する（ステップＳ１０１）。なお、初期化が必要な場合としては、例えば、モデルパラメータ推定処理を初回に実行する場合（つまり、ｋ＝１のとき）、外部（例えばオペレータ等のユーザ）からの初期化指示があった場合等が挙げられる。

上記のステップＳ１０１で初期化すると判定された場合、モデルパラメータ推定部１０２は、モデルパラメータ推定値θ_ｅｓｔ（ｋ−１）＝θ_０、共分散行列Ｐ（ｋ−１）＝Ｉと初期化する（ステップＳ１０２）。特に、ｋ＝１のとき、θ_ｅｓｔ（０）＝θ_０、Ｐ（０）＝Ｉである。

なお、モデルパラメータ初期値θ_０としては、例えば、プラント応答関数に既に設定されているパラメータを用いてもよいし、固定の初期値（例えば、零ベクトル等）を用いてもよいし、その他の任意の方法で決定された初期値を用いてもよい。また、Ｉとしては、例えば、単位行列を用いてもよいし、予め決められた任意の行列を用いてもよいし、その他の任意の方法で決定された行列を用いてもよい。

上記のステップＳ１０１で初期化すると判定されなかった場合又は上記のステップＳ１０２に続いて、モデルパラメータ推定部１０２は、以下により予測誤差ε（ｋ）を計算する（ステップＳ１０３）。

ｙ（ｋ）＝ｙ_０
ε（ｋ）＝ｙ（ｋ）−φ（ｋ）^Ｔθ_ｅｓｔ（ｋ−１）
ここで、Ｔは行列又はベクトルの転置を表す記号である。また、φ（ｋ）は状態ベクトルであり、例えば、以下で表される。

すなわち、状態ベクトルφ（ｋ）は、過去のＮ個の制御量ｙ（ｋ−１）〜ｙ（ｋ−Ｎ）と過去のＭ個の操作量ｕ（ｋ−１）〜ｕ（ｋ−Ｍ）とを順に並べたＮ＋Ｍ次元ベクトルである。なお、Ｎ及びＭは予め設定される１以上の整数である。

ただし、上記の数１に示す状態ベクトルφ（ｋ）は外乱ｖを考慮しない場合の状態ベクトルφ（ｋ）であって、上記の数１に示す状態ベクトルφ（ｋ）の代わりに、外乱ｖを考慮した状態ベクトルφ（ｋ）が用いられてもよい。外乱ｖを考慮する場合、状態ベクトルφ（ｋ）は、過去のＮ個の制御量ｙ（ｋ−１）〜ｙ（ｋ−Ｎ）と過去のＭ個の操作量ｕ（ｋ−１）〜ｕ（ｋ−Ｍ）と過去のＫ個の外乱ｖ（ｋ−１）〜ｖ（ｋ−Ｋ）とを順に並べたＮ＋Ｍ＋Ｋ次元ベクトルとなる。なお、Ｋは予め設定される１以上の整数である。

なお、状態ベクトルφ（ｋ）の次元数をＬ（Ｌ＝Ｎ＋Ｍ又はＬ＝Ｎ＋Ｍ＋Ｋ）とすれば、モデルパラメータ推定値θ_ｅｓｔ（ｋ）はＬ次元ベクトル、共分散行列Ｐ（ｋ）はＬ×Ｌ行列である。

次に、モデルパラメータ推定部１０２は、状態ベクトルφ（ｋ）を用いて、以下により共分散行列Ｐ（ｋ）を更新する（ステップＳ１０４）。

次に、モデルパラメータ推定部１０２は、モデルパラメータ推定値θ_ｅｓｔを更新するか否かを判定する（ステップＳ１０５）。なお、モデルパラメータ推定値θ_ｅｓｔを更新する場合としては、例えば、上記のステップＳ１０２の初期化が行われた後に所定の時間が経過した場合（つまり、当該初期化後に予め決められた時刻ステップが経過した場合）、外部（例えばオペレータ等のユーザ）からの更新指示があった場合等が挙げられる。

上記のステップＳ１０５でモデルパラメータ推定値θ_ｅｓｔを更新すると判定されなかった場合、モデルパラメータ推定部１０２は、モデルパラメータ推定値θ_ｅｓｔ（ｋ）を前回の値のままとする（ステップＳ１０６）。すなわち、モデルパラメータ推定部１０２は、θ_ｅｓｔ（ｋ）＝θ_ｅｓｔ（ｋ−１）とする。

一方で、上記のステップＳ１０５でモデルパラメータ推定値θ_ｅｓｔを更新すると判定された場合、モデルパラメータ推定部１０２は、状態ベクトルφ（ｋ）と共分散行列Ｐ（ｋ−１）と予測誤差ε（ｋ）と前回のモデルパラメータ推定値θ_ｅｓｔ（ｋ−１）とを用いて、以下によりモデルパラメータ推定値θ_ｅｓｔ（ｋ）を更新（算出）する（ステップＳ１０７）。

そして、モデルパラメータ推定部１０２は、状態ベクトルφ（ｋ）を状態ベクトルφ（ｋ＋１）に更新する（ステップＳ１０８）。すなわち、モデルパラメータ推定部１０２は、以下の状態ベクトルφ（ｋ＋１）を作成する。

ここで、状態ベクトルφ（ｋ＋１）に含まれるｙ（ｋ）にはｙ_０が設定される。また、状態ベクトルφ（ｋ＋１）に含まれるｙ（ｋ−１）〜ｙ（ｋ−Ｎ＋１）には、状態ベクトルφ（ｋ）に含まれるｙ（ｋ−１）〜ｙ（ｋ−Ｎ＋１）がそれぞれ設定される。

同様に、状態ベクトルφ（ｋ＋１）に含まれるｕ（ｋ）にはｕ_０が設定される。また、状態ベクトルφ（ｋ＋１）に含まれるｕ（ｋ−１）〜ｕ（ｋ−Ｍ＋１）には、状態ベクトルφ（ｋ）に含まれるｕ（ｋ−１）〜ｕ（ｋ−Ｍ＋１）がそれぞれ設定される。

なお、外乱ｖを考慮する場合、状態ベクトルφ（ｋ＋１）に含まれるｖ（ｋ）にはｖ_０が設定される。また、状態ベクトルφ（ｋ＋１）に含まれるｖ（ｋ−１）〜ｖ（ｋ−Ｋ＋１）には、状態ベクトルφ（ｋ）に含まれるｖ（ｋ−１）〜ｖ（ｋ−Ｋ＋１）がそれぞれ設定される。

以上により、本実施形態に係るモデル作成支援装置１０は、ｋ＝１，２，・・・に対してモデルパラメータ推定値θ_ｅｓｔ（ｋ）を逐次的に算出することができる。これらのモデルパラメータ推定値θ_ｅｓｔ（ｋ）は逐次的に算出される都度、加速計算部１０３と表示装置５０にそれぞれ出力される。

＜加速値推定処理＞
次に、加速計算部１０３が加速推定値θ_ａｃｃ（ｋ）を逐次的に算出する場合について説明する。加速計算部１０３は、ｋ＝３，４，・・・に対して加速値推定処理を繰り返し実行することで、加速推定値θ_ａｃｃ（ｋ）を逐次的に算出する。以下、或るｋ（ただし、ｋは３以上）における加速値推定処理について、図４を参照しながら説明する。図４は、本実施形態に係る加速値推定処理の一例を示すフローチャートである。なお、加速計算部１０３は、例えば、モデルパラメータ推定部１０２からモデルパラメータ推定値θ_ｅｓｔ（ｋ）（ただし、ｋは３以上）を受信する毎に加速値推定処理を繰り返し実行すればよい。

ここで、以下では、少なくとも過去３点のモデルパラメータ推定値θ_ｅｓｔ（ｋ−１）、θ_ｅｓｔ（ｋ−２）及びθ_ｅｓｔ（ｋ−３）がモデルパラメータ履歴として補助記憶装置２０６等に格納されているものとする。ただし、加速値推定処理の初回実行時（つまり、ｋ＝３のとき）は、θ_ｅｓｔ（２）及びθ_ｅｓｔ（１）の過去２点のみであってもよい。

まず、加速計算部１０３は、以下によりモデルパラメータ履歴［θ_ｅｓｔ（ｋ−１），θ_ｅｓｔ（ｋ−２），θ_ｅｓｔ（ｋ−３）］を［θ_ｅｓｔ（ｋ），θ_ｅｓｔ（ｋ−１），θ_ｅｓｔ（ｋ−２）］に更新する（ステップＳ２０１）。

次に、加速計算部１０３は、モデルパラメータ推定値θ_ｅｓｔ（ｋ）の各要素の番号を表すインデックスｉを１に初期化する（ステップＳ２０２）。

次に、加速計算部１０３は、モデルパラメータ推定値θ_ｅｓｔ（ｋ）のｉ番目の要素の加速値を計算する（ステップＳ２０３）ここで、加速計算部１０３は、以下のＳｔｅｐ２−１〜Ｓｔｅｐ２−３により、モデルパラメータ推定値θ_ｅｓｔ（ｋ）のｉ番目の要素の加速値を計算する。なお、以下、モデルパラメータ推定値θ_ｅｓｔ（ｋ）のｉ番目の要素を「θ_{ｅｓｔ，ｉ}（ｋ）」、その加速値を「θ_{ａｃｃ，ｉ}（ｋ）」とも表す。

Ｓｔｅｐ２−１：加速計算部１０３は、モデルパラメータ履歴［θ_ｅｓｔ（ｋ），θ_ｅｓｔ（ｋ−１），θ_ｅｓｔ（ｋ−２）］を用いて、以下のｘ_１，ｘ_２，ｘ_３，ｄ及びｎを定義する。

ｘ_３：＝θ_{ｅｓｔ，ｉ}（ｋ）
ｘ_２：＝θ_{ｅｓｔ，ｉ}（ｋ−１）
ｘ_１：＝θ_{ｅｓｔ，ｉ}（ｋ−２）
ｄ：＝ｘ_３−２ｘ_２＋ｘ_１
ｎ：＝（ｘ_３−ｘ_２）^２
Ｓｔｅｐ２−２：加速計算部１０３は、エイトケンのΔ２乗加速法により加速値ｘ_４を計算する。このとき、加速計算部１０３は、｜ｄ｜＜εを満たす場合はｘ_４：＝ｘ_３とし、そうでない場合はｘ_４：＝ｘ_３−ｎ／ｄとする。ここで、εは予め設定された調整パラメータである。なお、本Ｓｔｅｐで、絶対値｜ｄ｜が調整パラメータεよりも小さい場合にｘ_４：＝ｘ_３としたのは、ゼロ割の回避と加速値ｘ_４が過剰に大きくなるのを防ぐためである。

Ｓｔｅｐ２−３：加速計算部１０３は、加速値θ_{ａｃｃ，ｉ}（ｋ）を計算する。このとき、加速計算部１０３は、加速値ｘ_４がｘ_３（つまり、θ_{ｅｓｔ，ｉ}（ｋ））と所定の比率α以上の乖離がある場合はｘ_４をθ_{ａｃｃ，ｉ}（ｋ）とし、そうでない場合はｘ_３をθ_{ａｃｃ，ｉ}（ｋ）とする。すなわち、加速計算部１０３は、｜（ｘ_４−ｘ_３）／ｘ_３｜＜αを満たす場合はθ_{ａｃｃ，ｉ}（ｋ）：＝ｘ_４とし、そうでない場合はθ_{ａｃｃ，ｉ}（ｋ）：＝ｘ_３とする。なお、αは予め設定されたパラメータである。本Ｓｔｅｐは、加速法により計算された加速値ｘ_４とモデルパラメータ推定処理で推定された値ｘ_３＝θ_{ｅｓｔ，ｉ}（ｋ）との乖離が大きい場合は、当該加速値ｘ_４をθ_{ａｃｃ，ｉ}（ｋ）として採用しないことを意味している。

なお、本実施形態では過去３点のモデルパラメータ推定値θ_ｅｓｔ（ｋ−１）、θ_ｅｓｔ（ｋ−２）及びθ_ｅｓｔ（ｋ−３）を用いてエイトケンのΔ２乗加速法により加速値を計算したが、これは一例であって、加速法はエイトケンのΔ２乗加速法に限られない。エイトケンのΔ２乗加速法以外にも任意の数列の加速法を用いることが可能である。

上記のステップＳ２０３に続いて、加速計算部１０３は、モデルパラメータ推定値θ_ｅｓｔ（ｋ）の全要素に対して加速値θ_{ａｃｃ，ｉ}（ｋ）を計算したか否かを判定する（ステップＳ２０４）。

上記のステップＳ２０４でモデルパラメータ推定値θ_ｅｓｔ（ｋ）の全要素に対して加速値θ_{ａｃｃ，ｉ}（ｋ）を計算したと判定されなかった場合、加速計算部１０３は、ｉに１を加算して（ステップＳ２０５）、上記のステップＳ２０３に戻る。これにより、モデルパラメータ推定値θ_ｅｓｔ（ｋ）のｉ番目の要素に対して加速値θ_{ａｃｃ，ｉ}（ｋ）が繰り返し計算される。

一方で、上記のステップＳ２０４でモデルパラメータ推定値θ_ｅｓｔ（ｋ）の全要素に対して加速値θ_{ａｃｃ，ｉ}（ｋ）を計算したと判定された場合、加速計算部１０３は、各加速値θ_{ａｃｃ，ｉ}（ｋ）をまとめて以下の加速推定値θ_ａｃｃ（ｋ）を得る（ステップＳ２０６）。

なお、Ｌはモデルパラメータ推定値θ_ｅｓｔ（ｋ）の要素数である。

以上により、本実施形態に係るモデル作成支援装置１０は、ｋ＝３，４，・・・に対してモデルパラメータ推定値θ_ｅｓｔ（ｋ）を加速法により加速させた加速推定値θ_ａｃｃ（ｋ）を逐次的に算出することができる。これらの加速推定値θ_ａｃｃ（ｋ）は逐次的に算出される都度、表示装置５０に出力される。

＜実施例１＞
次に、実施例１について説明する。実施例１では、外乱ｖを考慮せずに、モデルパラメータ推定値θ_ｅｓｔ（ｋ）と加速推定値θ_ａｃｃ（ｋ）を得る場合について説明する。

制御対象であるプラント３０の真のプラント応答は、ΔＴを時間メッシュとして、以下で表されるものとする。

ｙ（ｔ）＝ａ_１ｙ（ｔ−ΔＴ）＋ａ_２ｙ（ｎ−２ΔＴ）＋ｂ_１ｕ（ｎ−ΔＴ）＋ｂ_２ｕ（ｎ−２ΔＴ）＋δ（ｔ）
ここで、δ（ｔ）はノイズを表す。

このとき、ΔＴ＝Ｔ_ｃ（ただし、Ｔ_ｃは制御周期）とすれば、プラント応答関数は２次のＡＲＭＡモデルとなり、
ｙ（ｋ）＝ａ_１ｙ（ｋ−１）＋ａ_２ｙ（ｋ−２）＋ｂ_１ｙ（ｋ−１）＋ｂ_２ｙ（ｋ−２）と表せる。

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

となる。

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

となる。

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

となる。

更に、制御量の推定値は、ｙ_ｅｓｔ（ｋ）＝φ（ｋ）^Ｔθ_ｅｓｔ（ｋ−１）で計算される。

実施例１における計測データ（つまり、制御量ｙ、操作量ｕ及び外乱ｖ）を図５に示す。図５（ａ）が制御量ｙ、図５（ｂ）が操作量ｕ及び外乱ｖである。図５（ｂ）に示すように、実施例１では時刻０で操作量ｕをステップ状に変化させるステップ応答試験を行っている。このとき、図５（ａ）に示すように、制御量ｙは操作量ｕがステップ入力された後、徐々に上昇しているが、時刻１７５００時点では未だ収束していない。なお、実施例１では外乱ｖを考慮しないため、ｖ＝０で固定としている。また、時刻の単位は分とし、時刻０を基準に過去の時刻は負の時刻とする。

このとき、本実施形態に係るモデル作成支援装置１０で推定したモデルパラメータ推定値θ_ｅｓｔと加速推定値θ_ａｃｃを図６に示す。図６（ａ）がθ_{ｅｓｔ，１}〜θ_{ｅｓｔ，４}、図６（ｂ）がθ_{ａｃｃ，１}〜θ_{ａｃｃ，４}である。

従来より行われてきたステップ応答モデルを作成するステップ応答試験では、制御量ｙが収束するまでの波形が必要であるため、時刻１７５００時点では不十分で、更にモデル同定試験を継続する必要がある。

一方で、本実施形態に係るモデル作成支援装置１０で推定したモデルパラメータ推定値θ_ｅｓｔは時刻５０００付近でほぼ収束し、更に、加速推定値θ_ａｃｃは時刻３０００付近で収束している。

ここで、本実施形態に係るモデル作成支援装置１０で推定したモデルパラメータ推定値θ_ｅｓｔと加速推定値θ_ａｃｃとを表示装置５０に同時表示した結果を図７に示す。実施例１では、４つのパラメータ（つまり、その真値がそれぞれａ_１，ａ_２，ｂ_１，ｂ_２である４つのパラメータ）を推定しているが、いずれも一定時間後に収束している。なお、モデルパラメータ推定値θ_ｅｓｔの１番目の要素θ_{ｅｓｔ，１}〜４番目の要素θ_{ｅｓｔ，４}がそれぞれａ_１，ａ_２，ｂ_１及びｂ_２の推定値であり、加速推定値θ_ａｃｃの１番目の要素θ_{ａｃｃ，１}〜４番目の要素θ_{ａｃｃ，４}がそれぞれθ_{ｅｓｔ，１}〜θ_{ｅｓｔ，４}の加速値である。

特に、パラメータθ_{ｅｓｔ，２}及びθ_{ａｃｃ，２}は収束までに要する時間が長いが、パラメータθ_{ｅｓｔ，２}がおよそ時刻６５００で収束するのに対して、パラメータθ_{ａｃｃ，２}はおよそ時刻３５００で収束している。すなわち、加速値θ_{ａｃｃ，２}は、θ_{ｅｓｔ，２}よりもおよそ３０００分程度早く収束しており、しかもその収束値は一致している。

したがって、本実施形態に係るモデル作成支援装置１０により推定された加速推定値θ_ａｃｃを用いることで、モデルパラメータ推定値θ_ｅｓｔを用いる場合と比較しておよそ３０００分、比率にして４６％もの時間を削減することができる。この時間削減効果は、通常のステップ応答試験と比較すれば、更に顕著であることを言うまでもない。

＜実施例２＞
次に、実施例２について説明する。実施例２では、外乱ｖを考慮して、モデルパラメータ推定値θ_ｅｓｔ（ｋ）と加速推定値θ_ａｃｃ（ｋ）を得る場合について説明する。なお、特に言及しない点については実施例１と同様であるものとする。

実施例２における計測データ（つまり、制御量ｙ、操作量ｕ及び外乱ｖ）を図８に示す。図８（ａ）が制御量ｙ、図８（ｂ）が操作量ｕ及び外乱ｖである。図８（ｂ）に示すように、実施例２では操作量ｕが波状に変動し、更に外乱ｖも波状に変動して同時に加わっている。このとき、図８（ａ）に示すように、制御量ｙは波状に変動しつつ徐々に上昇し、時刻１７５００では収束傾向すら見えていない。したがって、波形のみではモデル同定試験を終了して良いのか否かを判断することができない状況と言える。

このとき、本実施形態に係るモデル作成支援装置１０で推定したモデルパラメータ推定値θ_ｅｓｔと加速推定値θ_ａｃｃを図９に示す。図９（ａ）がθ_{ｅｓｔ，１}〜θ_{ｅｓｔ，４}、図９（ｂ）がθ_{ａｃｃ，１}〜θ_{ａｃｃ，４}である。図９に示すように、モデルパラメータ推定値θ_ｅｓｔは時刻７０００付近でほぼ収束し、更に加速推定値θ_ａｃｃは時刻５０００で収束している。

ここで、本実施形態に係るモデル作成支援装置１０で推定したモデルパラメータ推定値θ_ｅｓｔと加速推定値θ_ａｃｃとを表示装置５０に同時表示した結果を図１０に示す。実施例２でも４つのパラメータを推定しているが、いずれも一定時間後に収束している。

特に、パラメータθ_{ｅｓｔ，２}及びθ_{ａｃｃ，２}は収束までに要する時間が長いが、パラメータθ_{ｅｓｔ，２}がおよそ時刻７０００で収束するのに対して、パラメータθ_{ａｃｃ，２}はおよそ時刻５０００で収束している。すなわち、加速値θ_{ａｃｃ，２}は、θ_{ｅｓｔ，２}よりもおよそ２０００分程度早く収束しており、しかもその収束値は一致している。

したがって、本実施形態に係るモデル作成支援装置１０により推定された加速推定値θ_ａｃｃを用いることで、モデルパラメータ推定値θ_ｅｓｔを用いる場合と比較しておよそ２０００分、比率にして４０％もの時間を削減することができる。

ここで、実施例２では、時刻３５００付近でモデルパラメータ推定値θ_ｅｓｔと加速推定値θ_ａｃｃとが一致するものの、時刻３７００〜時刻５０００にかけてその値の乖離が生じている様子が確認できる。これは、パラメータθ_{ｅｓｔ，２}及びθ_{ａｃｃ，２}が一旦０．２５付近として推定されつつあったものの、実はデータが不十分であり、真値とは遠い値に落ち着きつつあったことを示している。

したがって、モデルパラメータ推定値θ_ｅｓｔと加速推定値θ_ａｃｃとの間で所定の閾値以上の乖離が存在する場合、表示装置５０上等にアラームを発報することで、モデル同定試験を継続する必要があることをユーザ等に知らせることができる。このため、例えば、アラーム発報によってモデル同定試験のデータが不十分であること等が報知されるため、当該試験を早期に打ち切ってしまうような誤りを防止することが可能となる。実施例２の場合は閾値の値にもよるが、例えば、時刻４０００付近程度でアラームが発報されることで、モデル同定試験を継続する必要があることをユーザ等に知らせることができる。

なお、上述したように、モデルパラメータ推定値θ_ｅｓｔと加速推定値θ_ａｃｃとの間に所定の閾値以上の乖離が存在するか否かの比較判定は、モデル作成支援装置１０で行われてもよいし、表示装置５０が演算装置を備える場合には当該表示装置５０で行われてもよい。

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

１ モデル作成支援システム
１０ モデル作成支援装置
２０ 操作装置
３０ プラント
４０ 計測装置
５０ 表示装置
１０１ フィルタ部
１０２ モデルパラメータ推定部
１０３ 加速計算部
２０１ 外部Ｉ／Ｆ
２０１ａ 記録媒体
２０２ 通信Ｉ／Ｆ
２０３ ＲＡＭ
２０４ ＲＯＭ
２０５ ＣＰＵ
２０６ 補助記憶装置
２０７ バス

Claims (7)

1. 制御対象となるプラントのモデルを同定するモデル作成支援装置であって、
   前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を逐次計算するパラメータ推定部と、
   前記パラメータの推定値をエイトケンのΔ２乗加速法により加速した加速値を逐次計算する加速計算部と、を有し
   前記パラメータの推定値と、前記加速値から決定された加速推定値とを表示装置に逐次出力する、ことを特徴とするモデル作成支援装置。
2. 前記加速計算部は、
   前記パラメータの推定値と前記加速値との間に所定の比率以上の乖離が存在する場合は、前記パラメータの推定値を前記加速推定値に決定し、
   前記パラメータの推定値と前記加速値との間に前記比率以上の乖離が存在しない場合は、前記加速値を前記加速推定値に決定する、ことを特徴とする請求項１に記載のモデル作成支援装置。
3. 前記加速計算部は、
   現在時点の前記パラメータの推定値である第３の推定値ｘ_３と、１つ前の時点の前記パラメータの推定値である第２の推定値ｘ_２と、２つ前の時点の前記パラメータの推定値である第１の推定値ｘ_１とを用いて、ｄ：＝ｘ_３−２ｘ_２＋ｘ_１とｎ：＝（ｘ_３−ｘ_２）^２とを計算し、
   前記ｄの絶対値が予め設定された調整パラメータεより小さい場合は、前記加速値を前記第３の推定値ｘ_３とし、
   前記ｄの絶対値が前記調整パラメータε以上である場合は、前記加速値を、ｘ_３−ｎ／ｄとする、ことを特徴とする請求項２に記載のモデル作成支援装置。
4. 前記モデル作成支援装置は、
   前記パラメータの推定値と前記加速推定値との間に所定の閾値以上の乖離が生じた場合、前記表示装置上にアラームを発報させる、ことを特徴とする請求項１乃至３の何れか一項に記載のモデル作成支援装置。
5. 前記モデル作成支援装置は、
   前記加速推定値の変化が一定期間の間、予め設定された範囲内に収まっている場合、前記表示装置上に前記アラームを発報させる、ことを特徴とする請求項４に記載のモデル作成支援装置。
6. 制御対象となるプラントのモデルを同定するモデル作成支援装置が、
   前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を逐次計算するパラメータ推定手順と、
   前記パラメータの推定値をエイトケンのΔ２乗加速法により加速した加速値を逐次計算する加速計算手順と、を実行し、
   前記パラメータの推定値と、前記加速値から決定された加速推定値とを表示装置に逐次出力する、ことを特徴とするモデル作成方法。
7. 制御対象となるプラントのモデルを同定するモデル作成支援装置に、
   前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を逐次計算するパラメータ推定手順と、
   前記パラメータの推定値をエイトケンのΔ２乗加速法により加速した加速値を逐次計算する加速計算手順と、を実行させ、
   前記パラメータの推定値と、前記加速値から決定された加速推定値とを表示装置に逐次出力する、ことを特徴とするプログラム。

Images