特許第6904473号(P6904473)IP Force 特許公報掲載プロジェクト 2022.1.31 β版

知財求人 - 知財ポータルサイト「IP Force」

▶ 富士電機株式会社の特許一覧

特許6904473モデル作成支援装置、モデル作成支援方法及びプログラム
<>
  • 特許6904473-モデル作成支援装置、モデル作成支援方法及びプログラム 図000010
  • 特許6904473-モデル作成支援装置、モデル作成支援方法及びプログラム 図000011
  • 特許6904473-モデル作成支援装置、モデル作成支援方法及びプログラム 図000012
  • 特許6904473-モデル作成支援装置、モデル作成支援方法及びプログラム 図000013
  • 特許6904473-モデル作成支援装置、モデル作成支援方法及びプログラム 図000014
  • 特許6904473-モデル作成支援装置、モデル作成支援方法及びプログラム 図000015
  • 特許6904473-モデル作成支援装置、モデル作成支援方法及びプログラム 図000016
  • 特許6904473-モデル作成支援装置、モデル作成支援方法及びプログラム 図000017
  • 特許6904473-モデル作成支援装置、モデル作成支援方法及びプログラム 図000018
  • 特許6904473-モデル作成支援装置、モデル作成支援方法及びプログラム 図000019
< >
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B1)
(11)【特許番号】6904473
(24)【登録日】2021年6月28日
(45)【発行日】2021年7月14日
(54)【発明の名称】モデル作成支援装置、モデル作成支援方法及びプログラム
(51)【国際特許分類】
   G05B 23/02 20060101AFI20210701BHJP
   G05B 13/04 20060101ALI20210701BHJP
【FI】
   G05B23/02 G
   G05B23/02 P
   G05B13/04
【請求項の数】7
【全頁数】18
(21)【出願番号】特願2020-204836(P2020-204836)
(22)【出願日】2020年12月10日
【審査請求日】2020年12月10日
【早期審査対象出願】
(73)【特許権者】
【識別番号】000005234
【氏名又は名称】富士電機株式会社
(74)【代理人】
【識別番号】100107766
【弁理士】
【氏名又は名称】伊東 忠重
(74)【代理人】
【識別番号】100070150
【弁理士】
【氏名又は名称】伊東 忠彦
(72)【発明者】
【氏名】丹下 吉雄
【審査官】 松本 泰典
(56)【参考文献】
【文献】 特開2009−282804(JP,A)
【文献】 特開2020−095352(JP,A)
【文献】 米国特許出願公開第2003/0028263(US,A1)
(58)【調査した分野】(Int.Cl.,DB名)
G05B 23/02
G05B 13/04
(57)【特許請求の範囲】
【請求項1】
制御対象となるプラントのモデルを同定するモデル作成支援装置であって、
前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を逐次計算するパラメータ推定部と、
前記パラメータの推定値をエイトケンのΔ2乗加速法により加速した加速値逐次計算する加速計算部と、を有し
前記パラメータの推定値と、前記加速値から決定された加速推定値とを表示装置に逐次出力する、ことを特徴とするモデル作成支援装置。
【請求項2】
前記加速計算部は、
記パラメータの推定値と前記加速値との間に所定の比率以上の乖離が存在する場合は、前記パラメータの推定値を前記加速推定値に決定し、
前記パラメータの推定値と前記加速値との間に前記比率以上の乖離が存在しない場合は、前記加速値を前記加速推定値に決定する、ことを特徴とする請求項1に記載のモデル作成支援装置。
【請求項3】
前記加速計算部は、
現在時点の前記パラメータの推定値である第3の推定値xと、1つ前の時点の前記パラメータの推定値である第2の推定値xと、2つ前の時点の前記パラメータの推定値である第1の推定値xとを用いて、d:=x−2x+xとn:=(x−xとを計算し、
前記dの絶対値が予め設定された調整パラメータεより小さい場合は、前記加速値を前記第3の推定値xとし、
前記dの絶対値が前記調整パラメータε以上である場合は、前記加速値を、x−n/dとする、ことを特徴とする請求項2に記載のモデル作成支援装置。
【請求項4】
モデル作成支援装置は、
前記パラメータの推定値と前記加速推定値との間に所定の閾値以上の乖離が生じた場合、前記表示装置上にアラームを発報させる、ことを特徴とする請求項1乃至3の何れか一項に記載のモデル作成支援装置。
【請求項5】
前記モデル作成支援装置は、
前記加速推定値の変化が一定期間の間、予め設定された範囲内に収まっている場合、前記表示装置上に前記アラームを発報させる、ことを特徴とする請求項4に記載のモデル作成支援装置。
【請求項6】
制御対象となるプラントのモデルを同定するモデル作成支援装置が、
前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を逐次計算するパラメータ推定手順と、
前記パラメータの推定値をエイトケンのΔ2乗加速法により加速した加速値逐次計算する加速計算手順と、を実行し、
前記パラメータの推定値と、前記加速値から決定された加速推定値とを表示装置に逐次出力する、ことを特徴とするモデル作成方法。
【請求項7】
制御対象となるプラントのモデルを同定するモデル作成支援装置に、
前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を逐次計算するパラメータ推定手順と、
前記パラメータの推定値をエイトケンのΔ2乗加速法により加速した加速値逐次計算する加速計算手順と、を実行させ、
前記パラメータの推定値と、前記加速値から決定された加速推定値とを表示装置に逐次出力する、ことを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、モデル作成支援装置、モデル作成支援方法及びプログラムに関する。
【背景技術】
【0002】
温調制御装置やPLC(Programmable Logic Controller)、DCS(Distributed Control System)等の制御装置、パーソナルコンピュータや組み込み制御機器上で実装される制御装置等が産業上広く利用されている。また、それに関連する技術として、制御対象の制御量を目標値に追従させることを目的とする様々な制御方式(例えば、PID(Proportional-Integral-Differential)制御、モデル予測制御、内部モデル制御、LQG(Linear-Quadratic-Gaussian)制御、H2制御、H∞制御等)が知られている。特に、モデル予測制御は、制御対象の状態空間モデルや将来の時間応答モデルを用いた最適化計算を逐次的に行うことで望ましい応答を得る方式であり、産業界で広く用いられている(例えば、非特許文献1)。
【0003】
なお、モデル予測制御等で用いる線形予測モデルに関してパラメータ同定を行う手法として、逐次最小2乗法、すなわち、RLS(Recursive Least Squares)法やカルマンフィルタ法等が知られている(例えば、非特許文献2及び非特許文献3)。また、収束する数列の収束速度を加速する手法として、エイトケンのΔ2乗加速法が知られている(非特許文献4)。
【先行技術文献】
【非特許文献】
【0004】
【非特許文献1】ヤン M. マチエヨフスキー,「モデル予測制御 制約のもとでの最適制御」,東京電機大学出版局,2005
【非特許文献2】足立 修一,「MATLABによる制御のためのシステム同定」,東京電機大学出版局,1996
【非特許文献3】足立 修一,「MATLABによる制御のための上級システム同定」,東京電機大学出版局,2004
【非特許文献4】長田直樹,「収束の加速法」,数理解析研究所講究録,京都大学 (1994), 880: 28-43
【発明の概要】
【発明が解決しようとする課題】
【0005】
ところで、プラントの制御に用いる応答モデルのパラメータを計測データから推定するために、操作量に対して或る形状の信号(例えば、ステップ信号やランプ信号、M系列信号等)を印加するモデル同定試験が従来から一般的に行われている。
【0006】
しかしながら、例えば、設備の規模が大きかったり、長時間の熱的・化学的反応を伴ったりするような時定数の長いプラントを対象とする場合、ステップ応答が収束するのに数時間以上を要する場合がある。
【0007】
また、操作量が複数存在する場合に、操作量毎にステップ応答試験を行うと、(操作変数の数)×(1操作変数あたりのステップ応答試験時間)のモデル同定試験時間が必要となる。
【0008】
更に、プラントのプロセスに投入する材料や外気温等の条件毎にモデルを分けて作成する必要がある場合、(操作変数の数)×(1操作変数あたりのステップ応答試験時間)×(条件数)のモデル同定試験時間が必要となる。
【0009】
このように、プラントの応答モデルを得るための同定試験には非常に多くの時間が掛かり、運用上のコストや負荷が大きい。
【0010】
本発明の一実施形態は、上記の点に鑑みてなされたもので、モデル同定試験に要する時間を削減することを目的とする。
【課題を解決するための手段】
【0011】
上記目的を達成するため、一実施形態に係るモデル作成支援装置10は、制御対象となるプラントのモデルを同定するモデル作成支援装置であって、前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を計算するパラメータ推定部と、前記パラメータの推定値を加速法により加速した加速推定値を計算する加速計算部と、を有することを特徴とする。
【発明の効果】
【0012】
モデル同定試験に要する時間を削減することができる。
【図面の簡単な説明】
【0013】
図1】本実施形態に係るモデル作成支援システムの全体構成の一例を示す図である。
図2】本実施形態に係るモデル作成支援装置のハードウェア構成の一例を示す図である。
図3】本実施形態に係るモデルパラメータ推定処理の一例を示すフローチャートである。
図4】本実施形態に係る加速値推定処理の一例を示すフローチャートである。
図5】実施例1における計測データを示す図である。
図6】実施例1におけるモデルパラメータ推定値及び加速推定値を示す図である。
図7】実施例1におけるモデルパラメータ推定値及び加速推定値を同時表示した場合を示す図である。
図8】実施例2における計測データを示す図である。
図9】実施例2におけるモデルパラメータ推定値及び加速推定値を示す図である。
図10】実施例2におけるモデルパラメータ推定値及び加速推定値を同時表示した場合を示す図である。
【発明を実施するための形態】
【0014】
以下、本発明の一実施形態について説明する。本実施形態では、モデル同定試験(つまり、各種プラント等の制御対象の制御に用いる応答モデルのパラメータを推定するための試験)に要する時間を削減することができるモデル作成支援システム1について説明する。
【0015】
<全体構成>
まず、本実施形態に係るモデル作成支援システム1の全体構成について、図1を参照しながら説明する。図1は、本実施形態に係るモデル作成支援システム1の全体構成の一例を示す図である。
【0016】
図1に示すように、本実施形態に係るモデル作成支援システム1には、モデル作成支援装置10と、操作装置20と、プラント30と、計測装置40と、表示装置50とが含まれる。
【0017】
モデル作成支援装置10は、所定の制御周期毎に計測装置40から計測データ(つまり、制御量y(t)、操作量u(t)及び外乱v(t))を受信して、これらの計測データを用いてプラント30の応答モデルのパラメータ(以下、「モデルパラメータ」ともいう。)を推定するコンピュータ又はコンピュータシステムである。なお、プラント30の応答モデルを表す関数を「プラント応答関数」ともいう(つまり、モデルパラメータとはプラント応答関数のパラメータのことを意味する。)。ここで、tは時刻ステップを表す。例えば、制御周期をT、現在の時刻ステップ(つまり、最新の時刻ステップ)をtとすれば、時刻ステップは、t,t−T,t−2T,・・・と表される。
【0018】
このとき、モデル作成支援装置10は、モデルパラメータの推定値(以下、「モデルパラメータ推定値」ともいう。)を加速法により加速させることで、通常よりも早くモデルパラメータを推定する。なお、以下、モデルパラメータ推定値を加速法により加速させた値を「加速推定値」ともいう。この加速推定値をモデルパラメータ値と扱うことで、通常のモデル同定試験よりも早期にモデル同定が可能となり、その試験時間を削減することが可能となる。
【0019】
操作装置20は、操作信号に応じて、所定の制御周期毎にプラント30に対して操作量uを出力する。プラント30は、制御対象の各種設備や工場、発電所、機器等である。計測装置40は、所定の制御周期毎にプラント30の制御量y(t)、当該プラント30に対する操作量u(t)及び当該プラント30に対する外乱v(t)を計測し、これら制御量y(t)、操作量u(t)及び外乱v(t)を計測データとしてモデル作成支援装置10に送信するセンサ等である。表示装置50は、モデル作成支援装置10で推定されたモデルパラメータ推定や加速推定値を表示する。なお、操作信号は、例えば、プラント30を制御する制御装置等から送信される。
【0020】
ここで、図1に示すように、本実施形態に係るモデル作成支援装置10は、フィルタ部101と、モデルパラメータ推定部102と、加速計算部103とを有する。
【0021】
フィルタ部101は、所定の制御周期毎に計測装置40から計測データ(制御量y(t)、操作量u(t)及び外乱v(t))を受信し、フィルタパラメータを用いて、これらの制御量y(t)、操作量u(t)及び外乱v(t)の時系列データ(つまり、{y(t)}、{u(t)}及び{v(t)})に対して所定のフィルタ処理を行う。このようなフィルタ処理を行うことで、計測データがデータ処理に適した形式に変換される。なお、フィルタ処理としては、例えば、移動平均やダウンサンプリング等が挙げられる。フィルタパラメータは、例えば、移動平均の区間やダウンサンプリングの間隔等の種々のパラメータである。以下、制御量y(t)、操作量u(t)及び外乱v(t)はフィルタ処理後の計測データであるものとする。また、時刻ステップt=t(つまり、現在の時刻ステップ)における制御量y(t)、操作量u(t)及び外乱v(t)をそれぞれy、u及びvとも表す。
【0022】
モデルパラメータ推定部102は、制御量y(t)と操作量u(t)と外乱v(t)とを入力として、モデルパラメータ推定値θest(k)を逐次的に算出する。ここで、モデルパラメータ推定部102には、モデルパラメータ推定値の初期値θが与えられる。以下、モデルパラメータ推定値の初期値を「モデルパラメータ初期値」ともいう。なお、kは初期値として1を取り、モデルパラメータ推定値θest(k)が算出される毎に1ずつ加算されるインデックスである。
【0023】
加速計算部103は、モデルパラメータ推定値θest(k)を入力として、加速法により、モデルパラメータの加速推定値θacc(k)を逐次的に算出する。
【0024】
ここで、モデルパラメータ推定値θest(k)及び加速推定値θacc(k)は表示装置50に出力され、当該表示装置50上に表示される。また、必要に応じて、当該表示装置50上でアラームが発報される。アラームが発報される場合としては、例えば、モデルパラメータ推定値θest(k)と加速推定値θacc(k)とを比較結果、所定の閾値以上の乖離が存在する場合等である。これは、モデルパラメータ推定値θest(k)と加速推定値θacc(k)とで乖離が大きい場合は、モデル同定試験のデータが不十分であることを示唆するため、当該試験を不適切に早期に打ち切ってしまうことを防止するためである。これにより、モデル同定試験の担当者等は、その試験時間を削減(短縮)することができると共に、短縮し過ぎない適切な試験時間の長さを知ることができる。なお、モデルパラメータ推定値θest(k)と加速推定値θacc(k)との間に所定の閾値以上の乖離が存在するか否かの比較判定は、モデル作成支援装置10で行われてもよいし、表示装置50が演算装置を備える場合には当該表示装置50で行われてもよい。
【0025】
また、上記以外にも、例えば、加速推定値θacc(k)の変化が一定期間の間、微小である場合(例えば、当該加速推定値θacc(k)の全ての要素の変化が予め設定された範囲内に収まっている場合等)に当該表示装置50上にアラームが発報されてもよい。これは、この場合、加速推定値θacc(k)が収束したと考えられるためである。
【0026】
<ハードウェア構成>
次に、本実施形態に係るモデル作成支援装置10のハードウェア構成について、図2を参照しながら説明する。図2は、本実施形態に係るモデル作成支援装置10のハードウェア構成の一例を示す図である。
【0027】
図2に示すように、本実施形態に係るモデル作成支援装置10は、外部I/F201と、通信I/F202と、RAM(Random Access Memory)203と、ROM(Read Only Memory)204と、CPU(Central Processing Unit)205と、補助記憶装置206とを有する。これら各ハードウェアは、バス207を介して相互に通信可能に接続されている。
【0028】
外部I/F201は、例えば記録媒体201a等の外部装置とのインタフェースである。モデル作成支援装置10は、外部I/F201を介して、記録媒体201aの読み取りや書き込み等を行うことができる。記録媒体201aには、例えば、SDメモリカード(SD memory card)やUSBメモリ、CD(Compact Disk)、DVD(Digital Versatile Disk)等がある。なお、記録媒体201aには、モデル作成支援装置10が有する各機能部(フィルタ部101、モデルパラメータ推定部102及び加速計算部103)を実現する1以上のプログラムが格納されていてもよい。
【0029】
通信I/F202は、モデル作成支援装置10が他の機器や装置等とデータ通信を行うためのインタフェースである。なお、モデル作成支援装置10が有する各機能部を実現する1以上のプログラムは、通信I/F202を介して、所定のサーバ等から取得(ダウンロード)されてもよい。
【0030】
RAM203は、プログラムやデータを一時保持する揮発性の半導体メモリである。ROM204は、電源を切ってもデータを保持することができる不揮発性の半導体メモリである。CPU205は、例えば補助記憶装置206やROM204等からプログラムやデータをRAM203上に読み出して、各種処理を実行する演算装置である。補助記憶装置206は、例えばHDD(Hard Disk Drive)やSSD(Solid State Drive)等の不揮発性の記憶装置である。
【0031】
本実施形態に係るモデル作成支援装置10は、図2に示すハードウェア構成を有することにより、後述する各種処理を実現することができる。なお、図2に示すハードウェア構成は一例であって、モデル作成支援装置10は、他のハードウェア構成を有していてもよい。例えば、モデル作成支援装置10は、GPU(Graphics Processing Unit)等の演算装置を有していてもよいし、複数のCPU205や複数の補助記憶装置206を有していてもよい。
【0032】
<モデルパラメータ推定処理>
次に、モデルパラメータ推定部102がモデルパラメータ推定値θest(k)を逐次的に算出する場合について説明する。モデルパラメータ推定部102は、k=1,2,・・・に対してモデルパラメータ推定処理を繰り返し実行することで、モデルパラメータ推定値θest(k)を逐次的に算出する。以下、プラント応答関数としてはARMA(autoregressive and moving average)モデルを想定した上で、或るkにおけるモデルパラメータ推定処理について、図3を参照しながら説明する。図3は、本実施形態に係るモデルパラメータ推定処理の一例を示すフローチャートである。なお、モデルパラメータ推定部102は、例えば、所定の制御周期(又は、予め決められた時間幅)毎にモデルパラメータ推定処理を繰り返し実行すればよい。
【0033】
まず、モデルパラメータ推定部102は、モデルパラメータ推定値θestと共分散行列Pとを初期化するか否かを判定する(ステップS101)。なお、初期化が必要な場合としては、例えば、モデルパラメータ推定処理を初回に実行する場合(つまり、k=1のとき)、外部(例えばオペレータ等のユーザ)からの初期化指示があった場合等が挙げられる。
【0034】
上記のステップS101で初期化すると判定された場合、モデルパラメータ推定部102は、モデルパラメータ推定値θest(k−1)=θ、共分散行列P(k−1)=Iと初期化する(ステップS102)。特に、k=1のとき、θest(0)=θ、P(0)=Iである。
【0035】
なお、モデルパラメータ初期値θとしては、例えば、プラント応答関数に既に設定されているパラメータを用いてもよいし、固定の初期値(例えば、零ベクトル等)を用いてもよいし、その他の任意の方法で決定された初期値を用いてもよい。また、Iとしては、例えば、単位行列を用いてもよいし、予め決められた任意の行列を用いてもよいし、その他の任意の方法で決定された行列を用いてもよい。
【0036】
上記のステップS101で初期化すると判定されなかった場合又は上記のステップS102に続いて、モデルパラメータ推定部102は、以下により予測誤差ε(k)を計算する(ステップS103)。
【0037】
y(k)=y
ε(k)=y(k)−φ(k)θest(k−1)
ここで、Tは行列又はベクトルの転置を表す記号である。また、φ(k)は状態ベクトルであり、例えば、以下で表される。
【0038】
【数1】

すなわち、状態ベクトルφ(k)は、過去のN個の制御量y(k−1)〜y(k−N)と過去のM個の操作量u(k−1)〜u(k−M)とを順に並べたN+M次元ベクトルである。なお、N及びMは予め設定される1以上の整数である。
【0039】
ただし、上記の数1に示す状態ベクトルφ(k)は外乱vを考慮しない場合の状態ベクトルφ(k)であって、上記の数1に示す状態ベクトルφ(k)の代わりに、外乱vを考慮した状態ベクトルφ(k)が用いられてもよい。外乱vを考慮する場合、状態ベクトルφ(k)は、過去のN個の制御量y(k−1)〜y(k−N)と過去のM個の操作量u(k−1)〜u(k−M)と過去のK個の外乱v(k−1)〜v(k−K)とを順に並べたN+M+K次元ベクトルとなる。なお、Kは予め設定される1以上の整数である。
【0040】
なお、状態ベクトルφ(k)の次元数をL(L=N+M又はL=N+M+K)とすれば、モデルパラメータ推定値θest(k)はL次元ベクトル、共分散行列P(k)はL×L行列である。
【0041】
次に、モデルパラメータ推定部102は、状態ベクトルφ(k)を用いて、以下により共分散行列P(k)を更新する(ステップS104)。
【0042】
【数2】
次に、モデルパラメータ推定部102は、モデルパラメータ推定値θestを更新するか否かを判定する(ステップS105)。なお、モデルパラメータ推定値θestを更新する場合としては、例えば、上記のステップS102の初期化が行われた後に所定の時間が経過した場合(つまり、当該初期化後に予め決められた時刻ステップが経過した場合)、外部(例えばオペレータ等のユーザ)からの更新指示があった場合等が挙げられる。
【0043】
上記のステップS105でモデルパラメータ推定値θestを更新すると判定されなかった場合、モデルパラメータ推定部102は、モデルパラメータ推定値θest(k)を前回の値のままとする(ステップS106)。すなわち、モデルパラメータ推定部102は、θest(k)=θest(k−1)とする。
【0044】
一方で、上記のステップS105でモデルパラメータ推定値θestを更新すると判定された場合、モデルパラメータ推定部102は、状態ベクトルφ(k)と共分散行列P(k−1)と予測誤差ε(k)と前回のモデルパラメータ推定値θest(k−1)とを用いて、以下によりモデルパラメータ推定値θest(k)を更新(算出)する(ステップS107)。
【0045】
【数3】
そして、モデルパラメータ推定部102は、状態ベクトルφ(k)を状態ベクトルφ(k+1)に更新する(ステップS108)。すなわち、モデルパラメータ推定部102は、以下の状態ベクトルφ(k+1)を作成する。
【0046】
【数4】
ここで、状態ベクトルφ(k+1)に含まれるy(k)にはyが設定される。また、状態ベクトルφ(k+1)に含まれるy(k−1)〜y(k−N+1)には、状態ベクトルφ(k)に含まれるy(k−1)〜y(k−N+1)がそれぞれ設定される。
【0047】
同様に、状態ベクトルφ(k+1)に含まれるu(k)にはuが設定される。また、状態ベクトルφ(k+1)に含まれるu(k−1)〜u(k−M+1)には、状態ベクトルφ(k)に含まれるu(k−1)〜u(k−M+1)がそれぞれ設定される。
【0048】
なお、外乱vを考慮する場合、状態ベクトルφ(k+1)に含まれるv(k)にはvが設定される。また、状態ベクトルφ(k+1)に含まれるv(k−1)〜v(k−K+1)には、状態ベクトルφ(k)に含まれるv(k−1)〜v(k−K+1)がそれぞれ設定される。
【0049】
以上により、本実施形態に係るモデル作成支援装置10は、k=1,2,・・・に対してモデルパラメータ推定値θest(k)を逐次的に算出することができる。これらのモデルパラメータ推定値θest(k)は逐次的に算出される都度、加速計算部103と表示装置50にそれぞれ出力される。
【0050】
<加速値推定処理>
次に、加速計算部103が加速推定値θacc(k)を逐次的に算出する場合について説明する。加速計算部103は、k=3,4,・・・に対して加速値推定処理を繰り返し実行することで、加速推定値θacc(k)を逐次的に算出する。以下、或るk(ただし、kは3以上)における加速値推定処理について、図4を参照しながら説明する。図4は、本実施形態に係る加速値推定処理の一例を示すフローチャートである。なお、加速計算部103は、例えば、モデルパラメータ推定部102からモデルパラメータ推定値θest(k)(ただし、kは3以上)を受信する毎に加速値推定処理を繰り返し実行すればよい。
【0051】
ここで、以下では、少なくとも過去3点のモデルパラメータ推定値θest(k−1)、θest(k−2)及びθest(k−3)がモデルパラメータ履歴として補助記憶装置206等に格納されているものとする。ただし、加速値推定処理の初回実行時(つまり、k=3のとき)は、θest(2)及びθest(1)の過去2点のみであってもよい。
【0052】
まず、加速計算部103は、以下によりモデルパラメータ履歴[θest(k−1),θest(k−2),θest(k−3)]を[θest(k),θest(k−1),θest(k−2)]に更新する(ステップS201)。
【0053】
次に、加速計算部103は、モデルパラメータ推定値θest(k)の各要素の番号を表すインデックスiを1に初期化する(ステップS202)。
【0054】
次に、加速計算部103は、モデルパラメータ推定値θest(k)のi番目の要素の加速値を計算する(ステップS203)ここで、加速計算部103は、以下のStep2−1〜Step2−3により、モデルパラメータ推定値θest(k)のi番目の要素の加速値を計算する。なお、以下、モデルパラメータ推定値θest(k)のi番目の要素を「θest,i(k)」、その加速値を「θacc,i(k)」とも表す。
【0055】
Step2−1:加速計算部103は、モデルパラメータ履歴[θest(k),θest(k−1),θest(k−2)]を用いて、以下のx,x,x,d及びnを定義する。
【0056】
:=θest,i(k)
:=θest,i(k−1)
:=θest,i(k−2)
d:=x−2x+x
n:=(x−x
Step2−2:加速計算部103は、エイトケンのΔ2乗加速法により加速値xを計算する。このとき、加速計算部103は、|d|<εを満たす場合はx:=xとし、そうでない場合はx:=x−n/dとする。ここで、εは予め設定された調整パラメータである。なお、本Stepで、絶対値|d|が調整パラメータεよりも小さい場合にx:=xとしたのは、ゼロ割の回避と加速値xが過剰に大きくなるのを防ぐためである。
【0057】
Step2−3:加速計算部103は、加速値θacc,i(k)を計算する。このとき、加速計算部103は、加速値xがx(つまり、θest,i(k))と所定の比率α以上の乖離がある場合はxをθacc,i(k)とし、そうでない場合はxをθacc,i(k)とする。すなわち、加速計算部103は、|(x−x)/x|<αを満たす場合はθacc,i(k):=xとし、そうでない場合はθacc,i(k):=xとする。なお、αは予め設定されたパラメータである。本Stepは、加速法により計算された加速値xとモデルパラメータ推定処理で推定された値x=θest,i(k)との乖離が大きい場合は、当該加速値xをθacc,i(k)として採用しないことを意味している。
【0058】
なお、本実施形態では過去3点のモデルパラメータ推定値θest(k−1)、θest(k−2)及びθest(k−3)を用いてエイトケンのΔ2乗加速法により加速値を計算したが、これは一例であって、加速法はエイトケンのΔ2乗加速法に限られない。エイトケンのΔ2乗加速法以外にも任意の数列の加速法を用いることが可能である。
【0059】
上記のステップS203に続いて、加速計算部103は、モデルパラメータ推定値θest(k)の全要素に対して加速値θacc,i(k)を計算したか否かを判定する(ステップS204)。
【0060】
上記のステップS204でモデルパラメータ推定値θest(k)の全要素に対して加速値θacc,i(k)を計算したと判定されなかった場合、加速計算部103は、iに1を加算して(ステップS205)、上記のステップS203に戻る。これにより、モデルパラメータ推定値θest(k)のi番目の要素に対して加速値θacc,i(k)が繰り返し計算される。
【0061】
一方で、上記のステップS204でモデルパラメータ推定値θest(k)の全要素に対して加速値θacc,i(k)を計算したと判定された場合、加速計算部103は、各加速値θacc,i(k)をまとめて以下の加速推定値θacc(k)を得る(ステップS206)。
【0062】
【数5】
なお、Lはモデルパラメータ推定値θest(k)の要素数である。
【0063】
以上により、本実施形態に係るモデル作成支援装置10は、k=3,4,・・・に対してモデルパラメータ推定値θest(k)を加速法により加速させた加速推定値θacc(k)を逐次的に算出することができる。これらの加速推定値θacc(k)は逐次的に算出される都度、表示装置50に出力される。
【0064】
<実施例1>
次に、実施例1について説明する。実施例1では、外乱vを考慮せずに、モデルパラメータ推定値θest(k)と加速推定値θacc(k)を得る場合について説明する。
【0065】
制御対象であるプラント30の真のプラント応答は、ΔTを時間メッシュとして、以下で表されるものとする。
【0066】
y(t)=ay(t−ΔT)+ay(n−2ΔT)+bu(n−ΔT)+bu(n−2ΔT)+δ(t)
ここで、δ(t)はノイズを表す。
【0067】
このとき、ΔT=T(ただし、Tは制御周期)とすれば、プラント応答関数は2次のARMAモデルとなり、
y(k)=ay(k−1)+ay(k−2)+by(k−1)+by(k−2)と表せる。
【0068】
このとき、状態ベクトルは、
【0069】
【数6】
となる。
【0070】
また、モデルパラメータ推定値は、
【0071】
【数7】
となる。
【0072】
また、モデルパラメータの真値は、
【0073】
【数8】
となる。
【0074】
更に、制御量の推定値は、yest(k)=φ(k)θest(k−1)で計算される。
【0075】
実施例1における計測データ(つまり、制御量y、操作量u及び外乱v)を図5に示す。図5(a)が制御量y、図5(b)が操作量u及び外乱vである。図5(b)に示すように、実施例1では時刻0で操作量uをステップ状に変化させるステップ応答試験を行っている。このとき、図5(a)に示すように、制御量yは操作量uがステップ入力された後、徐々に上昇しているが、時刻17500時点では未だ収束していない。なお、実施例1では外乱vを考慮しないため、v=0で固定としている。また、時刻の単位は分とし、時刻0を基準に過去の時刻は負の時刻とする。
【0076】
このとき、本実施形態に係るモデル作成支援装置10で推定したモデルパラメータ推定値θestと加速推定値θacc図6に示す。図6(a)がθest,1〜θest,4図6(b)がθacc,1〜θacc,4である。
【0077】
従来より行われてきたステップ応答モデルを作成するステップ応答試験では、制御量yが収束するまでの波形が必要であるため、時刻17500時点では不十分で、更にモデル同定試験を継続する必要がある。
【0078】
一方で、本実施形態に係るモデル作成支援装置10で推定したモデルパラメータ推定値θestは時刻5000付近でほぼ収束し、更に、加速推定値θaccは時刻3000付近で収束している。
【0079】
ここで、本実施形態に係るモデル作成支援装置10で推定したモデルパラメータ推定値θestと加速推定値θaccとを表示装置50に同時表示した結果を図7に示す。実施例1では、4つのパラメータ(つまり、その真値がそれぞれa,a,b,bである4つのパラメータ)を推定しているが、いずれも一定時間後に収束している。なお、モデルパラメータ推定値θestの1番目の要素θest,1〜4番目の要素θest,4がそれぞれa,a,b及びbの推定値であり、加速推定値θaccの1番目の要素θacc,1〜4番目の要素θacc,4がそれぞれθest,1〜θest,4の加速値である。
【0080】
特に、パラメータθest,2及びθacc,2は収束までに要する時間が長いが、パラメータθest,2がおよそ時刻6500で収束するのに対して、パラメータθacc,2はおよそ時刻3500で収束している。すなわち、加速値θacc,2は、θest,2よりもおよそ3000分程度早く収束しており、しかもその収束値は一致している。
【0081】
したがって、本実施形態に係るモデル作成支援装置10により推定された加速推定値θaccを用いることで、モデルパラメータ推定値θestを用いる場合と比較しておよそ3000分、比率にして46%もの時間を削減することができる。この時間削減効果は、通常のステップ応答試験と比較すれば、更に顕著であることを言うまでもない。
【0082】
<実施例2>
次に、実施例2について説明する。実施例2では、外乱vを考慮して、モデルパラメータ推定値θest(k)と加速推定値θacc(k)を得る場合について説明する。なお、特に言及しない点については実施例1と同様であるものとする。
【0083】
実施例2における計測データ(つまり、制御量y、操作量u及び外乱v)を図8に示す。図8(a)が制御量y、図8(b)が操作量u及び外乱vである。図8(b)に示すように、実施例2では操作量uが波状に変動し、更に外乱vも波状に変動して同時に加わっている。このとき、図8(a)に示すように、制御量yは波状に変動しつつ徐々に上昇し、時刻17500では収束傾向すら見えていない。したがって、波形のみではモデル同定試験を終了して良いのか否かを判断することができない状況と言える。
【0084】
このとき、本実施形態に係るモデル作成支援装置10で推定したモデルパラメータ推定値θestと加速推定値θacc図9に示す。図9(a)がθest,1〜θest,4図9(b)がθacc,1〜θacc,4である。図9に示すように、モデルパラメータ推定値θestは時刻7000付近でほぼ収束し、更に加速推定値θaccは時刻5000で収束している。
【0085】
ここで、本実施形態に係るモデル作成支援装置10で推定したモデルパラメータ推定値θestと加速推定値θaccとを表示装置50に同時表示した結果を図10に示す。実施例2でも4つのパラメータを推定しているが、いずれも一定時間後に収束している。
【0086】
特に、パラメータθest,2及びθacc,2は収束までに要する時間が長いが、パラメータθest,2がおよそ時刻7000で収束するのに対して、パラメータθacc,2はおよそ時刻5000で収束している。すなわち、加速値θacc,2は、θest,2よりもおよそ2000分程度早く収束しており、しかもその収束値は一致している。
【0087】
したがって、本実施形態に係るモデル作成支援装置10により推定された加速推定値θaccを用いることで、モデルパラメータ推定値θestを用いる場合と比較しておよそ2000分、比率にして40%もの時間を削減することができる。
【0088】
ここで、実施例2では、時刻3500付近でモデルパラメータ推定値θestと加速推定値θaccとが一致するものの、時刻3700〜時刻5000にかけてその値の乖離が生じている様子が確認できる。これは、パラメータθest,2及びθacc,2が一旦0.25付近として推定されつつあったものの、実はデータが不十分であり、真値とは遠い値に落ち着きつつあったことを示している。
【0089】
したがって、モデルパラメータ推定値θestと加速推定値θaccとの間で所定の閾値以上の乖離が存在する場合、表示装置50上等にアラームを発報することで、モデル同定試験を継続する必要があることをユーザ等に知らせることができる。このため、例えば、アラーム発報によってモデル同定試験のデータが不十分であること等が報知されるため、当該試験を早期に打ち切ってしまうような誤りを防止することが可能となる。実施例2の場合は閾値の値にもよるが、例えば、時刻4000付近程度でアラームが発報されることで、モデル同定試験を継続する必要があることをユーザ等に知らせることができる。
【0090】
なお、上述したように、モデルパラメータ推定値θestと加速推定値θaccとの間に所定の閾値以上の乖離が存在するか否かの比較判定は、モデル作成支援装置10で行われてもよいし、表示装置50が演算装置を備える場合には当該表示装置50で行われてもよい。
【0091】
本発明は、具体的に開示された上記の実施形態に限定されるものではなく、特許請求の範囲の記載から逸脱することなく、種々の変形や変更、既知の技術との組み合わせ等が可能である。
【符号の説明】
【0092】
1 モデル作成支援システム
10 モデル作成支援装置
20 操作装置
30 プラント
40 計測装置
50 表示装置
101 フィルタ部
102 モデルパラメータ推定部
103 加速計算部
201 外部I/F
201a 記録媒体
202 通信I/F
203 RAM
204 ROM
205 CPU
206 補助記憶装置
207 バス
【要約】
【課題】モデル同定試験に要する時間を削減すること。
【解決手段】一実施形態に係るモデル作成支援装置10は、制御対象となるプラントのモデルを同定するモデル作成支援装置であって、前記プラントの計測データに含まれる制御量と操作量とに少なくとも基づいて、前記プラントの応答関数が持つパラメータの推定値を計算するパラメータ推定部と、前記パラメータの推定値を加速法により加速した加速推定値を計算する加速計算部と、を有することを特徴とする。
【選択図】図1
図1
図2
図3
図4
図5
図6
図7
図8
図9
図10