(57)【特許請求の範囲】
【請求項１】
制御対象となるプラントのモデルを同定するモデル作成支援装置であって、
前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を逐次計算するパラメータ推定部と、
前記パラメータの推定値をエイトケンのΔ２乗加速法により加速した加速値を逐次計算する加速計算部と、を有し
前記パラメータの推定値と、前記加速値から決定された加速推定値とを表示装置に逐次出力する、ことを特徴とするモデル作成支援装置。
【請求項２】
前記加速計算部は、
前記パラメータの推定値と前記加速値との間に所定の比率以上の乖離が存在する場合は、前記パラメータの推定値を前記加速推定値に決定し、
前記パラメータの推定値と前記加速値との間に前記比率以上の乖離が存在しない場合は、前記加速値を前記加速推定値に決定する、ことを特徴とする請求項１に記載のモデル作成支援装置。
【請求項３】
前記加速計算部は、
現在時点の前記パラメータの推定値である第３の推定値ｘ３と、１つ前の時点の前記パラメータの推定値である第２の推定値ｘ２と、２つ前の時点の前記パラメータの推定値である第１の推定値ｘ１とを用いて、ｄ：＝ｘ３−２ｘ２＋ｘ１とｎ：＝（ｘ３−ｘ２）２とを計算し、
前記ｄの絶対値が予め設定された調整パラメータεより小さい場合は、前記加速値を前記第３の推定値ｘ３とし、
前記ｄの絶対値が前記調整パラメータε以上である場合は、前記加速値を、ｘ３−ｎ／ｄとする、ことを特徴とする請求項２に記載のモデル作成支援装置。
【請求項４】
前記モデル作成支援装置は、
前記パラメータの推定値と前記加速推定値との間に所定の閾値以上の乖離が生じた場合、前記表示装置上にアラームを発報させる、ことを特徴とする請求項１乃至３の何れか一項に記載のモデル作成支援装置。
【請求項５】
前記モデル作成支援装置は、
前記加速推定値の変化が一定期間の間、予め設定された範囲内に収まっている場合、前記表示装置上に前記アラームを発報させる、ことを特徴とする請求項４に記載のモデル作成支援装置。
【請求項６】
制御対象となるプラントのモデルを同定するモデル作成支援装置が、
前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を逐次計算するパラメータ推定手順と、
前記パラメータの推定値をエイトケンのΔ２乗加速法により加速した加速値を逐次計算する加速計算手順と、を実行し、
前記パラメータの推定値と、前記加速値から決定された加速推定値とを表示装置に逐次出力する、ことを特徴とするモデル作成方法。
【請求項７】
制御対象となるプラントのモデルを同定するモデル作成支援装置に、
前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を逐次計算するパラメータ推定手順と、
前記パラメータの推定値をエイトケンのΔ２乗加速法により加速した加速値を逐次計算する加速計算手順と、を実行させ、
前記パラメータの推定値と、前記加速値から決定された加速推定値とを表示装置に逐次出力する、ことを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、モデル作成支援装置、モデル作成支援方法及びプログラムに関する。
【背景技術】
【０００２】
温調制御装置やＰＬＣ（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）モデルを想定した上で、或るｋにおけるモデルパラメータ推定処理について、図３を参照しながら説明する。図３は、本実施形態に係るモデルパラメータ推定処理の一例を示すフローチャートである。なお、モデルパラメータ推定部１０２は、例えば、所定の制御周期（又は、予め決められた時間幅）毎にモデルパラメータ推定処理を繰り返し実行すればよい。
【００３３】
まず、モデルパラメータ推定部１０２は、モデルパラメータ推定値θｅｓｔと共分散行列Ｐとを初期化するか否かを判定する（ステップＳ１０１）。なお、初期化が必要な場合としては、例えば、モデルパラメータ推定処理を初回に実行する場合（つまり、ｋ＝１のとき）、外部（例えばオペレータ等のユーザ）からの初期化指示があった場合等が挙げられる。
【００３４】
上記のステップＳ１０１で初期化すると判定された場合、モデルパラメータ推定部１０２は、モデルパラメータ推定値θｅｓｔ（ｋ−１）＝θ０、共分散行列Ｐ（ｋ−１）＝Ｉと初期化する（ステップＳ１０２）。特に、ｋ＝１のとき、θｅｓｔ（０）＝θ０、Ｐ（０）＝Ｉである。
【００３５】
なお、モデルパラメータ初期値θ０としては、例えば、プラント応答関数に既に設定されているパラメータを用いてもよいし、固定の初期値（例えば、零ベクトル等）を用いてもよいし、その他の任意の方法で決定された初期値を用いてもよい。また、Ｉとしては、例えば、単位行列を用いてもよいし、予め決められた任意の行列を用いてもよいし、その他の任意の方法で決定された行列を用いてもよい。
【００３６】
上記のステップＳ１０１で初期化すると判定されなかった場合又は上記のステップＳ１０２に続いて、モデルパラメータ推定部１０２は、以下により予測誤差ε（ｋ）を計算する（ステップＳ１０３）。
【００３７】
ｙ（ｋ）＝ｙ０
ε（ｋ）＝ｙ（ｋ）−φ（ｋ）Ｔθｅｓｔ（ｋ−１）
ここで、Ｔは行列又はベクトルの転置を表す記号である。また、φ（ｋ）は状態ベクトルであり、例えば、以下で表される。
【００３８】
【数１】
0006904473.tif 000002
すなわち、状態ベクトルφ（ｋ）は、過去のＮ個の制御量ｙ（ｋ−１）〜ｙ（ｋ−Ｎ）と過去のＭ個の操作量ｕ（ｋ−１）〜ｕ（ｋ−Ｍ）とを順に並べたＮ＋Ｍ次元ベクトルである。なお、Ｎ及びＭは予め設定される１以上の整数である。
【００３９】
ただし、上記の数１に示す状態ベクトルφ（ｋ）は外乱ｖを考慮しない場合の状態ベクトルφ（ｋ）であって、上記の数１に示す状態ベクトルφ（ｋ）の代わりに、外乱ｖを考慮した状態ベクトルφ（ｋ）が用いられてもよい。外乱ｖを考慮する場合、状態ベクトルφ（ｋ）は、過去のＮ個の制御量ｙ（ｋ−１）〜ｙ（ｋ−Ｎ）と過去のＭ個の操作量ｕ（ｋ−１）〜ｕ（ｋ−Ｍ）と過去のＫ個の外乱ｖ（ｋ−１）〜ｖ（ｋ−Ｋ）とを順に並べたＮ＋Ｍ＋Ｋ次元ベクトルとなる。なお、Ｋは予め設定される１以上の整数である。
【００４０】
なお、状態ベクトルφ（ｋ）の次元数をＬ（Ｌ＝Ｎ＋Ｍ又はＬ＝Ｎ＋Ｍ＋Ｋ）とすれば、モデルパラメータ推定値θｅｓｔ（ｋ）はＬ次元ベクトル、共分散行列Ｐ（ｋ）はＬ×Ｌ行列である。
【００４１】
次に、モデルパラメータ推定部１０２は、状態ベクトルφ（ｋ）を用いて、以下により共分散行列Ｐ（ｋ）を更新する（ステップＳ１０４）。
【００４２】
【数２】
0006904473.tif 000003
次に、モデルパラメータ推定部１０２は、モデルパラメータ推定値θｅｓｔを更新するか否かを判定する（ステップＳ１０５）。なお、モデルパラメータ推定値θｅｓｔを更新する場合としては、例えば、上記のステップＳ１０２の初期化が行われた後に所定の時間が経過した場合（つまり、当該初期化後に予め決められた時刻ステップが経過した場合）、外部（例えばオペレータ等のユーザ）からの更新指示があった場合等が挙げられる。
【００４３】
上記のステップＳ１０５でモデルパラメータ推定値θｅｓｔを更新すると判定されなかった場合、モデルパラメータ推定部１０２は、モデルパラメータ推定値θｅｓｔ（ｋ）を前回の値のままとする（ステップＳ１０６）。すなわち、モデルパラメータ推定部１０２は、θｅｓｔ（ｋ）＝θｅｓｔ（ｋ−１）とする。
【００４４】
一方で、上記のステップＳ１０５でモデルパラメータ推定値θｅｓｔを更新すると判定された場合、モデルパラメータ推定部１０２は、状態ベクトルφ（ｋ）と共分散行列Ｐ（ｋ−１）と予測誤差ε（ｋ）と前回のモデルパラメータ推定値θｅｓｔ（ｋ−１）とを用いて、以下によりモデルパラメータ推定値θｅｓｔ（ｋ）を更新（算出）する（ステップＳ１０７）。
【００４５】
【数３】
0006904473.tif 000004
そして、モデルパラメータ推定部１０２は、状態ベクトルφ（ｋ）を状態ベクトルφ（ｋ＋１）に更新する（ステップＳ１０８）。すなわち、モデルパラメータ推定部１０２は、以下の状態ベクトルφ（ｋ＋１）を作成する。
【００４６】
【数４】
0006904473.tif 000005
ここで、状態ベクトルφ（ｋ＋１）に含まれるｙ（ｋ）にはｙ０が設定される。また、状態ベクトルφ（ｋ＋１）に含まれるｙ（ｋ−１）〜ｙ（ｋ−Ｎ＋１）には、状態ベクトルφ（ｋ）に含まれるｙ（ｋ−１）〜ｙ（ｋ−Ｎ＋１）がそれぞれ設定される。
【００４７】
同様に、状態ベクトルφ（ｋ＋１）に含まれるｕ（ｋ）にはｕ０が設定される。また、状態ベクトルφ（ｋ＋１）に含まれるｕ（ｋ−１）〜ｕ（ｋ−Ｍ＋１）には、状態ベクトルφ（ｋ）に含まれるｕ（ｋ−１）〜ｕ（ｋ−Ｍ＋１）がそれぞれ設定される。
【００４８】
なお、外乱ｖを考慮する場合、状態ベクトルφ（ｋ＋１）に含まれるｖ（ｋ）にはｖ０が設定される。また、状態ベクトルφ（ｋ＋１）に含まれるｖ（ｋ−１）〜ｖ（ｋ−Ｋ＋１）には、状態ベクトルφ（ｋ）に含まれるｖ（ｋ−１）〜ｖ（ｋ−Ｋ＋１）がそれぞれ設定される。
【００４９】
以上により、本実施形態に係るモデル作成支援装置１０は、ｋ＝１，２，・・・に対してモデルパラメータ推定値θｅｓｔ（ｋ）を逐次的に算出することができる。これらのモデルパラメータ推定値θｅｓｔ（ｋ）は逐次的に算出される都度、加速計算部１０３と表示装置５０にそれぞれ出力される。
【００５０】
＜加速値推定処理＞
次に、加速計算部１０３が加速推定値θａｃｃ（ｋ）を逐次的に算出する場合について説明する。加速計算部１０３は、ｋ＝３，４，・・・に対して加速値推定処理を繰り返し実行することで、加速推定値θａｃｃ（ｋ）を逐次的に算出する。以下、或るｋ（ただし、ｋは３以上）における加速値推定処理について、図４を参照しながら説明する。図４は、本実施形態に係る加速値推定処理の一例を示すフローチャートである。なお、加速計算部１０３は、例えば、モデルパラメータ推定部１０２からモデルパラメータ推定値θｅｓｔ（ｋ）（ただし、ｋは３以上）を受信する毎に加速値推定処理を繰り返し実行すればよい。
【００５１】
ここで、以下では、少なくとも過去３点のモデルパラメータ推定値θｅｓｔ（ｋ−１）、θｅｓｔ（ｋ−２）及びθｅｓｔ（ｋ−３）がモデルパラメータ履歴として補助記憶装置２０６等に格納されているものとする。ただし、加速値推定処理の初回実行時（つまり、ｋ＝３のとき）は、θｅｓｔ（２）及びθｅｓｔ（１）の過去２点のみであってもよい。
【００５２】
まず、加速計算部１０３は、以下によりモデルパラメータ履歴［θｅｓｔ（ｋ−１），θｅｓｔ（ｋ−２），θｅｓｔ（ｋ−３）］を［θｅｓｔ（ｋ），θｅｓｔ（ｋ−１），θｅｓｔ（ｋ−２）］に更新する（ステップＳ２０１）。
【００５３】
次に、加速計算部１０３は、モデルパラメータ推定値θｅｓｔ（ｋ）の各要素の番号を表すインデックスｉを１に初期化する（ステップＳ２０２）。
【００５４】
次に、加速計算部１０３は、モデルパラメータ推定値θｅｓｔ（ｋ）のｉ番目の要素の加速値を計算する（ステップＳ２０３）ここで、加速計算部１０３は、以下のＳｔｅｐ２−１〜Ｓｔｅｐ２−３により、モデルパラメータ推定値θｅｓｔ（ｋ）のｉ番目の要素の加速値を計算する。なお、以下、モデルパラメータ推定値θｅｓｔ（ｋ）のｉ番目の要素を「θｅｓｔ，ｉ（ｋ）」、その加速値を「θａｃｃ，ｉ（ｋ）」とも表す。
【００５５】
Ｓｔｅｐ２−１：加速計算部１０３は、モデルパラメータ履歴［θｅｓｔ（ｋ），θｅｓｔ（ｋ−１），θｅｓｔ（ｋ−２）］を用いて、以下のｘ１，ｘ２，ｘ３，ｄ及びｎを定義する。
【００５６】
ｘ３：＝θｅｓｔ，ｉ（ｋ）
ｘ２：＝θｅｓｔ，ｉ（ｋ−１）
ｘ１：＝θｅｓｔ，ｉ（ｋ−２）
ｄ：＝ｘ３−２ｘ２＋ｘ１
ｎ：＝（ｘ３−ｘ２）２
Ｓｔｅｐ２−２：加速計算部１０３は、エイトケンのΔ２乗加速法により加速値ｘ４を計算する。このとき、加速計算部１０３は、｜ｄ｜＜εを満たす場合はｘ４：＝ｘ３とし、そうでない場合はｘ４：＝ｘ３−ｎ／ｄとする。ここで、εは予め設定された調整パラメータである。なお、本Ｓｔｅｐで、絶対値｜ｄ｜が調整パラメータεよりも小さい場合にｘ４：＝ｘ３としたのは、ゼロ割の回避と加速値ｘ４が過剰に大きくなるのを防ぐためである。
【００５７】
Ｓｔｅｐ２−３：加速計算部１０３は、加速値θａｃｃ，ｉ（ｋ）を計算する。このとき、加速計算部１０３は、加速値ｘ４がｘ３（つまり、θｅｓｔ，ｉ（ｋ））と所定の比率α以上の乖離がある場合はｘ４をθａｃｃ，ｉ（ｋ）とし、そうでない場合はｘ３をθａｃｃ，ｉ（ｋ）とする。すなわち、加速計算部１０３は、｜（ｘ４−ｘ３）／ｘ３｜＜αを満たす場合はθａｃｃ，ｉ（ｋ）：＝ｘ４とし、そうでない場合はθａｃｃ，ｉ（ｋ）：＝ｘ３とする。なお、αは予め設定されたパラメータである。本Ｓｔｅｐは、加速法により計算された加速値ｘ４とモデルパラメータ推定処理で推定された値ｘ３＝θｅｓｔ，ｉ（ｋ）との乖離が大きい場合は、当該加速値ｘ４をθａｃｃ，ｉ（ｋ）として採用しないことを意味している。
【００５８】
なお、本実施形態では過去３点のモデルパラメータ推定値θｅｓｔ（ｋ−１）、θｅｓｔ（ｋ−２）及びθｅｓｔ（ｋ−３）を用いてエイトケンのΔ２乗加速法により加速値を計算したが、これは一例であって、加速法はエイトケンのΔ２乗加速法に限られない。エイトケンのΔ２乗加速法以外にも任意の数列の加速法を用いることが可能である。
【００５９】
上記のステップＳ２０３に続いて、加速計算部１０３は、モデルパラメータ推定値θｅｓｔ（ｋ）の全要素に対して加速値θａｃｃ，ｉ（ｋ）を計算したか否かを判定する（ステップＳ２０４）。
【００６０】
上記のステップＳ２０４でモデルパラメータ推定値θｅｓｔ（ｋ）の全要素に対して加速値θａｃｃ，ｉ（ｋ）を計算したと判定されなかった場合、加速計算部１０３は、ｉに１を加算して（ステップＳ２０５）、上記のステップＳ２０３に戻る。これにより、モデルパラメータ推定値θｅｓｔ（ｋ）のｉ番目の要素に対して加速値θａｃｃ，ｉ（ｋ）が繰り返し計算される。
【００６１】
一方で、上記のステップＳ２０４でモデルパラメータ推定値θｅｓｔ（ｋ）の全要素に対して加速値θａｃｃ，ｉ（ｋ）を計算したと判定された場合、加速計算部１０３は、各加速値θａｃｃ，ｉ（ｋ）をまとめて以下の加速推定値θａｃｃ（ｋ）を得る（ステップＳ２０６）。
【００６２】
【数５】
0006904473.tif 000006
なお、Ｌはモデルパラメータ推定値θｅｓｔ（ｋ）の要素数である。
【００６３】
以上により、本実施形態に係るモデル作成支援装置１０は、ｋ＝３，４，・・・に対してモデルパラメータ推定値θｅｓｔ（ｋ）を加速法により加速させた加速推定値θａｃｃ（ｋ）を逐次的に算出することができる。これらの加速推定値θａｃｃ（ｋ）は逐次的に算出される都度、表示装置５０に出力される。
【００６４】
＜実施例１＞
次に、実施例１について説明する。実施例１では、外乱ｖを考慮せずに、モデルパラメータ推定値θｅｓｔ（ｋ）と加速推定値θａｃｃ（ｋ）を得る場合について説明する。
【００６５】
制御対象であるプラント３０の真のプラント応答は、ΔＴを時間メッシュとして、以下で表されるものとする。
【００６６】
ｙ（ｔ）＝ａ１ｙ（ｔ−ΔＴ）＋ａ２ｙ（ｎ−２ΔＴ）＋ｂ１ｕ（ｎ−ΔＴ）＋ｂ２ｕ（ｎ−２ΔＴ）＋δ（ｔ）
ここで、δ（ｔ）はノイズを表す。
【００６７】
このとき、ΔＴ＝Ｔｃ（ただし、Ｔｃは制御周期）とすれば、プラント応答関数は２次のＡＲＭＡモデルとなり、
ｙ（ｋ）＝ａ１ｙ（ｋ−１）＋ａ２ｙ（ｋ−２）＋ｂ１ｙ（ｋ−１）＋ｂ２ｙ（ｋ−２）と表せる。
【００６８】
このとき、状態ベクトルは、
【００６９】
【数６】
0006904473.tif 000007
となる。
【００７０】
また、モデルパラメータ推定値は、
【００７１】
【数７】
0006904473.tif 000008
となる。
【００７２】
また、モデルパラメータの真値は、
【００７３】
【数８】
0006904473.tif 000009
となる。
【００７４】
更に、制御量の推定値は、ｙｅｓｔ（ｋ）＝φ（ｋ）Ｔθｅｓｔ（ｋ−１）で計算される。
【００７５】
実施例１における計測データ（つまり、制御量ｙ、操作量ｕ及び外乱ｖ）を図５に示す。図５（ａ）が制御量ｙ、図５（ｂ）が操作量ｕ及び外乱ｖである。図５（ｂ）に示すように、実施例１では時刻０で操作量ｕをステップ状に変化させるステップ応答試験を行っている。このとき、図５（ａ）に示すように、制御量ｙは操作量ｕがステップ入力された後、徐々に上昇しているが、時刻１７５００時点では未だ収束していない。なお、実施例１では外乱ｖを考慮しないため、ｖ＝０で固定としている。また、時刻の単位は分とし、時刻０を基準に過去の時刻は負の時刻とする。
【００７６】
このとき、本実施形態に係るモデル作成支援装置１０で推定したモデルパラメータ推定値θｅｓｔと加速推定値θａｃｃを図６に示す。図６（ａ）がθｅｓｔ，１〜θｅｓｔ，４、図６（ｂ）がθａｃｃ，１〜θａｃｃ，４である。
【００７７】
従来より行われてきたステップ応答モデルを作成するステップ応答試験では、制御量ｙが収束するまでの波形が必要であるため、時刻１７５００時点では不十分で、更にモデル同定試験を継続する必要がある。
【００７８】
一方で、本実施形態に係るモデル作成支援装置１０で推定したモデルパラメータ推定値θｅｓｔは時刻５０００付近でほぼ収束し、更に、加速推定値θａｃｃは時刻３０００付近で収束している。
【００７９】
ここで、本実施形態に係るモデル作成支援装置１０で推定したモデルパラメータ推定値θｅｓｔと加速推定値θａｃｃとを表示装置５０に同時表示した結果を図７に示す。実施例１では、４つのパラメータ（つまり、その真値がそれぞれａ１，ａ２，ｂ１，ｂ２である４つのパラメータ）を推定しているが、いずれも一定時間後に収束している。なお、モデルパラメータ推定値θｅｓｔの１番目の要素θｅｓｔ，１〜４番目の要素θｅｓｔ，４がそれぞれａ１，ａ２，ｂ１及びｂ２の推定値であり、加速推定値θａｃｃの１番目の要素θａｃｃ，１〜４番目の要素θａｃｃ，４がそれぞれθｅｓｔ，１〜θｅｓｔ，４の加速値である。
【００８０】
特に、パラメータθｅｓｔ，２及びθａｃｃ，２は収束までに要する時間が長いが、パラメータθｅｓｔ，２がおよそ時刻６５００で収束するのに対して、パラメータθａｃｃ，２はおよそ時刻３５００で収束している。すなわち、加速値θａｃｃ，２は、θｅｓｔ，２よりもおよそ３０００分程度早く収束しており、しかもその収束値は一致している。
【００８１】
したがって、本実施形態に係るモデル作成支援装置１０により推定された加速推定値θａｃｃを用いることで、モデルパラメータ推定値θｅｓｔを用いる場合と比較しておよそ３０００分、比率にして４６％もの時間を削減することができる。この時間削減効果は、通常のステップ応答試験と比較すれば、更に顕著であることを言うまでもない。
【００８２】
＜実施例２＞
次に、実施例２について説明する。実施例２では、外乱ｖを考慮して、モデルパラメータ推定値θｅｓｔ（ｋ）と加速推定値θａｃｃ（ｋ）を得る場合について説明する。なお、特に言及しない点については実施例１と同様であるものとする。
【００８３】
実施例２における計測データ（つまり、制御量ｙ、操作量ｕ及び外乱ｖ）を図８に示す。図８（ａ）が制御量ｙ、図８（ｂ）が操作量ｕ及び外乱ｖである。図８（ｂ）に示すように、実施例２では操作量ｕが波状に変動し、更に外乱ｖも波状に変動して同時に加わっている。このとき、図８（ａ）に示すように、制御量ｙは波状に変動しつつ徐々に上昇し、時刻１７５００では収束傾向すら見えていない。したがって、波形のみではモデル同定試験を終了して良いのか否かを判断することができない状況と言える。
【００８４】
このとき、本実施形態に係るモデル作成支援装置１０で推定したモデルパラメータ推定値θｅｓｔと加速推定値θａｃｃを図９に示す。図９（ａ）がθｅｓｔ，１〜θｅｓｔ，４、図９（ｂ）がθａｃｃ，１〜θａｃｃ，４である。図９に示すように、モデルパラメータ推定値θｅｓｔは時刻７０００付近でほぼ収束し、更に加速推定値θａｃｃは時刻５０００で収束している。
【００８５】
ここで、本実施形態に係るモデル作成支援装置１０で推定したモデルパラメータ推定値θｅｓｔと加速推定値θａｃｃとを表示装置５０に同時表示した結果を図１０に示す。実施例２でも４つのパラメータを推定しているが、いずれも一定時間後に収束している。
【００８６】
特に、パラメータθｅｓｔ，２及びθａｃｃ，２は収束までに要する時間が長いが、パラメータθｅｓｔ，２がおよそ時刻７０００で収束するのに対して、パラメータθａｃｃ，２はおよそ時刻５０００で収束している。すなわち、加速値θａｃｃ，２は、θｅｓｔ，２よりもおよそ２０００分程度早く収束しており、しかもその収束値は一致している。
【００８７】
したがって、本実施形態に係るモデル作成支援装置１０により推定された加速推定値θａｃｃを用いることで、モデルパラメータ推定値θｅｓｔを用いる場合と比較しておよそ２０００分、比率にして４０％もの時間を削減することができる。
【００８８】
ここで、実施例２では、時刻３５００付近でモデルパラメータ推定値θｅｓｔと加速推定値θａｃｃとが一致するものの、時刻３７００〜時刻５０００にかけてその値の乖離が生じている様子が確認できる。これは、パラメータθｅｓｔ，２及びθａｃｃ，２が一旦０．２５付近として推定されつつあったものの、実はデータが不十分であり、真値とは遠い値に落ち着きつつあったことを示している。
【００８９】
したがって、モデルパラメータ推定値θｅｓｔと加速推定値θａｃｃとの間で所定の閾値以上の乖離が存在する場合、表示装置５０上等にアラームを発報することで、モデル同定試験を継続する必要があることをユーザ等に知らせることができる。このため、例えば、アラーム発報によってモデル同定試験のデータが不十分であること等が報知されるため、当該試験を早期に打ち切ってしまうような誤りを防止することが可能となる。実施例２の場合は閾値の値にもよるが、例えば、時刻４０００付近程度でアラームが発報されることで、モデル同定試験を継続する必要があることをユーザ等に知らせることができる。
【００９０】
なお、上述したように、モデルパラメータ推定値θｅｓｔと加速推定値θａｃｃとの間に所定の閾値以上の乖離が存在するか否かの比較判定は、モデル作成支援装置１０で行われてもよいし、表示装置５０が演算装置を備える場合には当該表示装置５０で行われてもよい。
【００９１】
本発明は、具体的に開示された上記の実施形態に限定されるものではなく、特許請求の範囲の記載から逸脱することなく、種々の変形や変更、既知の技術との組み合わせ等が可能である。
【符号の説明】
【００９２】
１        モデル作成支援システム
１０      モデル作成支援装置
２０      操作装置
３０      プラント
４０      計測装置
５０      表示装置
１０１    フィルタ部
１０２    モデルパラメータ推定部
１０３    加速計算部
２０１    外部Ｉ／Ｆ
２０１ａ  記録媒体
２０２    通信Ｉ／Ｆ
２０３    ＲＡＭ
２０４    ＲＯＭ
２０５    ＣＰＵ
２０６    補助記憶装置
２０７    バス
【要約】
【課題】モデル同定試験に要する時間を削減すること。
【解決手段】一実施形態に係るモデル作成支援装置１０は、制御対象となるプラントのモデルを同定するモデル作成支援装置であって、前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を計算するパラメータ推定部と、前記パラメータの推定値を加速法により加速した加速推定値を計算する加速計算部と、を有することを特徴とする。
【選択図】図１【図１】
0006904473.tif 000010
【図２】
0006904473.tif 000011
【図３】
0006904473.tif 000012
【図４】
0006904473.tif 000013
【図５】
0006904473.tif 000014
【図６】
0006904473.tif 000015
【図７】
0006904473.tif 000016
【図８】
0006904473.tif 000017
【図９】
0006904473.tif 000018
【図１０】
0006904473.tif 000019