以下、本発明にかかる実施の一形態を図面に基づいて説明する。なお、各図において同一の符号を付した構成は、同一の構成であることを示し、適宜、その説明を省略する。
図1は、実施形態におけるモデリング装置の構成を示すブロック図である。実施形態のモデリング装置Sは、モデリングすべきシステムのモデルを前記モデルのパラメータの値を決定することによって構築する装置であって、前記モデルは、連続した時間と見なされる時間変数を含む数式であり、予め得られている前記システムに関する先験的情報を制約条件として設定する制約条件設定部と、所定の入力データを実績入力データとして前記システムに入力した場合に前記システムから出力された実績出力データと、前記実績入力データを前記モデルに入力した場合に前記モデルから出力されたモデル出力データとの誤差を評価する評価関数の評価結果を求める評価演算部と、前記制約条件設定部で設定された制約条件を満たし、かつ、前記評価演算部で演算された評価結果が最良となるように、前記モデルのパラメータの値を求めるモデル演算部とを備える装置である。
このようなモデリング装置Sは、例えば、図1に示すように、演算部1と、入力設定部2と、出力部3と、記憶部4と、これら演算部1、入力設定部2、出力部3および記憶部4の間で相互にデータを交換することができるようにこれらを互いに接続するバス5とを備えて構成される。
入力設定部2は、モデリング装置Sの演算開始指示等の各種コマンドや、例えば実績入出力データ等のシステムのモデル化に必要な各種データをモデリング装置Sに入力する機器であり、例えば、キーボードやマウス等である。そして、本実施形態では、入力設定部2は、制約条件を入力設定する制約条件入力設定部21と、モデル次数および同定パラメータを入力設定するモデル次数同定パラメータ入力設定部22と、最適化計算パラメータを入力設定する最適化計算パラメータ入力設定部23と、評価関数を入力設定する評価関数入力設定部24とを備えている。これら制約条件、モデル次数、同定パラメータおよび最適化計算パラメータは、システムのモデル化の際に利用されるように、モデリング装置Sの記憶部4に格納されることで、モデリング装置Sに設定される。
制約条件入力設定部21は、モデルの制約条件(先験的情報)を入力してモデリング装置Sに設定するためのものである。制約条件は、システムのモデル化の際に、このモデル化によって得られたモデルが満たすべき条件であり、本実施形態では、モデリング対象のシステムに関する先験的知識によって与えられる。先験的知識は、任意の手段によってシステムから得られる情報であり、例えば、システムの設計値や理論式等から計算によって推定される情報であって良く、また例えば、システムを設計した場合に過去の事例に照らして経験的に推測される情報であって良く、また例えば、システムを実際に組み、この実際のシステムを稼動させてこれを観察することによって得られる情報であって良い。この得られた先験的知識は、モデリング装置Sが取り扱うことができる形式で表現(表記、記述)され、先験的情報とされる。
本実施形態では、予め決められた規範(基準)となるデータであって時間に依存する規範入力を入力データとしてシステムに与えた場合に、システムから出力される出力データ(計算によって得られる出力データ、および、システムで観測される実績出力データを含む)が制約条件(時間応答制約条件)とされる。規範入力は、例えば、ステップ状の入力(ステップ入力)、ランプ状の入力(ランプ入力)、三角関数波(例えばsin波)や三角波等の周期入力およびそれらの周波数掃引波(スイープ波、sweep波)等の、時間に関する所定の関数である。時間応答制約条件は、例えば前記出力データのそのものであって良く、また例えば前記出力データの上限値および/またはその下限値であって良く、さらにこれらに加えてその時間範囲が設定されて良い。すなわち、規範入力の入力時点を時刻0として、どの時刻tSからどの時刻tEまでの前記出力データが制約条件とされるのか、あるいは、どの時刻tSからどの時刻tEまでの前記出力データにおける上限値および/または下限値が制約条件とされるのか、が設定される(0≦tS<tE)。この時間応答制約条件が入力設定される場合、制約条件入力設定部21は、時間応答制約条件入力設定部として機能する。
なお、Aおよび/またはBは、AおよびBのうちの少なくとも一方の意味であり、AあるいはB、または、AとBとの両方の意味である。
そして、本実施形態では、システムの定常ゲインの符号や、その範囲(定常ゲインの上限値および下限値)が制約条件(定常ゲイン制約条件)とされる。この定常ゲイン制約条件が入力設定される場合、制約条件入力設定部21は、定常ゲイン制約条件入力設定部として機能する。
モデル次数未知パラメータ入力設定部22は、モデルの次数および前記モデルの同定パラメータを入力してモデリング装置Sに設定するためのものである。モデルの同定パラメータは、前記モデルにおけるパラメータのうち未知のパラメータとして同定すべきパラメータである。モデルの同定パラメータは、後述のシステムのモデルを作成する演算によってその値が決定され、これによってシステムのモデルが構築される。これらモデルの次数および同定パラメータは、例えば制御系設計に使用するためやシミュレーションに使用するため等の、モデリングの目的によって適宜に選定される。また、モデルの精度、モデル化に使用される入出力データの質や量、および、モデル化の情報処理量(計算量)等も勘案されて良い。
最適化計算パラメータ入力設定部23は、最適化計算パラメータを入力してモデリング装置Sに設定するためのものである。最適化計算パラメータは、最適化の計算を実行する際に必要とされる調整パラメータである。最適化計算パラメータは、本実施形態では、最適化計算手法として、後述するようにPSO(Particle Swarm Optimization)法が用いられていることから、このPSO法の演算に必要なパラメータであり、粒子pの個数P、解の更新回数(繰り返し回数)K、初期値の値域、後述の式におけるw、c1、c2である。最適化計算パラメータの入力設定の際に、モデルの次数や同定パラメータの個数が勘案されて良い。例えば、モデルの次数が大きくなるに従って、あるいは、同定パラメータの個数が多くなるに従って、探索すべき解空間が広くなるので、粒子pの個数Pや解の更新回数Kが大きく設定されることが好ましい。
評価関数入力設定部24は、評価関数を入力してモデリング装置Sに設定するためのものである。評価関数は、同定パラメータのパラメータ値を決めた場合に、このパラメータ値によって規定されるモデルがシステムをどの程度よく表現しているか否かの度合い(モデルの精度)を示す値(評価値)を求める関数である。評価関数は、予め1つだけ設定され、プログラムに予め組み込まれていても良いが、本実施形態では、より適切にモデルを評価するために、評価関数入力設定部24から入力可能とされている。さらに、オペレータの便宜を図るために、本実施形態では、予め複数の評価関数が用意されており、これら複数の評価関数の中から評価関数入力設定部24を介して選定することで、評価関数が入力設定される。
出力部3は、入力部2から入力されたコマンドやデータおよびモデリング装置Sの演算結果(システムのモデル)等を出力する機器であり、例えばCRTディスプレイ、LCD、有機ELディスプレイ又はプラズマディスプレイ等の表示装置やプリンタ等の印字装置等である。
記憶部4は、データやプログラムを記憶する例えばROM(Read Only Memory)およびEEPROM(Electrically Erasable Programmable Read Only Memory)等の不揮発性の記憶素子やハードディスク等を備えるとともに、演算部1が実行するモデリングプログラムや制御プログラム等の実行中における各データを一時的に記憶するいわゆるワーキングメモリであり、例えば揮発性の記憶素子であるRAM(Random Access Memory)を備えて構成される。記憶部4は、機能的に、制約条件記憶部41と、モデル次数同定パラメータ記憶部42と、最適化計算パラメータ記憶部43と、評価関数記憶部44と、モデル出力等記憶部45と、実績データ記憶部46とを備える。
制約条件記憶部41は、制約条件を記憶するものであり、制約条件入力設定部21から入力された制約条件が制約条件記憶部41に格納され、記憶される。制約条件記憶部41は、本実施形態では、時間応答制約条件を記憶する時間応答制約条件記憶部411と、定常ゲイン制約条件を記憶する定常ゲイン制約条件記憶部412とを備えている。
モデル次数同定パラメータ記憶部42は、モデルの次数および前記モデルの同定パラメータを記憶するものであり、モデル次数同定パラメータ入力設定部22から入力されたモデルの次数および同定パラメータがモデル次数同定パラメータ記憶部42に格納され、記憶される。
最適化計算パラメータ記憶部43は、最適化計算パラメータを記憶するものであり、最適化計算パラメータ入力設定部23から入力された最適化計算パラメータが最適化計算パラメータ記憶部43に格納され、記憶される。
評価関数記憶部44は、評価関数を記憶するものであり、評価関数入力設定部24から入力された評価関数が評価関数記憶部44に格納され、記憶される。本実施形態では、上述したように、複数の所定の評価関数が予め記憶されており、これら複数の評価関数の中から1つの評価関数が評価関数入力設定部24から選定され、この選定されたことを示す情報も評価関数記憶部44に記憶される。
モデル出力等記憶部45は、モデルを規定する所定のパラメータ値に対し、演算部1で演算されたモデル出力、評価結果および制約条件を満たすか否かの情報を互いに対応付けて記憶するものである。前記所定のパラメータ値は、本実施形態では、PSO法によって与えられる値である。そして、本実施形態では、制約条件を満たすか否かの情報は、評価結果を示す評価値に変えられ、評価結果としてモデル出力等記憶部45に記憶される。このように制約条件を満たすか否かの情報が評価結果を示す評価値に変換されて評価結果としてモデル出力等記憶部45に記憶されることで、モデルのパラメータの値を決定する際に、評価結果を参酌するだけで済み、このモデルのパラメータの値を決定する処理が簡便化される。
実績データ記憶部46は、実績データを記憶するものである。実績データは、モデリング対象のシステムに入力(印加)される所定の入力データ(実績入力データ)、および、この入力データの入力によってシステムから出力される出力データ(実績出力データ)であり、これら実績入力データおよび実績出力データ(実績入出力データ)は、互いに対応付けて実績データ記憶部43に記憶される。前記所定の入力データは、任意のデータであって良く、例えば、ステップ状の入力データ、ランダム入力データおよび擬似白色系列(例えばM系列等)データ等の時系列データである。出力データもこれに対応して時系列データとなる。前記所定の入力データは、制御工学の分野では、同定入力と呼ばれる。時系列入力データおよび時系列出力データは、1組であっても良く、また複数の組であっても良い。複数の組である場合には、時系列入出力データの各組は、同日に得られたデータであっても良く、また、異なる日に得られたデータであって良い。そして、これら複数の組の時系列入出力データは、モデリングの際に、これらの中から適宜に選択されるようにモデリング装置Sが構成されてもよい。
演算部1は、例えば、マイクロプロセッサおよびその周辺回路等を備えて構成され、プログラムに従い入力設定部2、出力部3および記憶部4を当該機能に応じてそれぞれ制御することによってモデリング装置S全体の制御を司るものであって、モデリング対象のシステムのモデルを前記モデルのパラメータの値を決定することによって構築するものである。演算部1は、機能的に、モデル出力演算部11と、制約条件評価部12と、評価演算部13と、モデル更新決定部14とを備える。
モデル出力演算部11は、モデル次数、同定パラメータおよび最適化計算パラメータを設定することによって設定された前記システムのモデルに、入力データを入力(印加)することで前記モデルの出力データ(モデル出力データ)を求めるものである。この求めたモデル出力データは、モデル出力等記憶部45に記憶され、格納される。
制約条件評価部12は、モデル出力演算部11によって得られたモデル出力データが制約条件を満たしているか否かを評価するものである。この評価したモデル出力データが制約条件を満たしているか否かの情報は、モデル出力等記憶部45に記憶され、格納される。
評価演算部13は、予め決められた所定の入力データを実績入力データとしてシステムに入力した場合にこのシステムから出力された実績出力データと、前記実績入力データを前記システムのモデルに入力した場合にこのモデルから出力されたモデル出力データとの誤差を評価する評価関数の評価結果を求めるものである。この求めた評価結果は、モデル出力等記憶部45に記憶され、格納される。
モデル更新決定部14は、前記システムのモデルを予め決められた所定の更新則に従って更新するとともに、制約条件入力設定部21で設定された制約条件を満たし、かつ、評価演算部13で演算された評価結果が最良となるように、前記システムのモデルにおけるそのパラメータの値を求め、前記システムのモデルを決定するものである。
このようなモデリング装置Sは、モデリング方法のプログラム(ソフトウェア)を実装することで、例えば、デスクトップ型パソコンやノート型パソコン等のコンピュータによって構成可能である。
なお、必要に応じてモデリング装置Sは、図1に破線で示す外部記憶部6や通信インターフェース部7をさらに備えてもよい。
外部記憶部6は、例えば、フレキシブルディスク、CD−ROM(Compact Disc Read Only Memory)、CD−R(Compact Disc Recordable)およびDVD−R(Digital Versatile Disc Recordable)等の記録媒体との間でデータを読み込みおよび/または書き込みを行う装置であり、例えば、フレキシブルディスクドライブ、CD−ROMドライブ、CD−Rドライブ及びDVD−Rドライブ等である。通信インターフェース部7は、通信網(ネットワーク)に接続され、この通信網を介して他の装置との間で通信信号を送受信するための機器である。
各プログラムが格納されていない場合には、これらを記録した記録媒体から外部記憶部6を介して記憶部4にインストールされるように構成してもよく、また、これらプログラムを管理するサーバコンピュータ(不図示)から通信網および通信インターフェース部7を介して各プログラムがダウンロードされるように構成してもよい。また、システムのモデル化の演算に当たってモデリング装置Sに入力すべき例えば実績入出力データ等のデータは、このデータを記憶した記録媒体によって外部記憶部6を介してモデリング装置Sに入力されるように構成してもよく、また、このデータを管理するサーバコンピュータから通信網および通信インターフェース部7を介してモデリング装置Sに入力されるように構成してもよい。
次に、本実施形態におけるモデリング装置の動作について説明する。図2は、実施形態におけるモデリング装置の概略動作を示すフローチャートである。図3は、実施形態のモデリング装置における制約条件の設定動作を示すフローチャートである。図4は、実施形態のモデリング装置におけるモデルの演算動作を示すフローチャートである。
本実施形態のモデリング装置Sでは、概略、図2に示すように、まず、ステップS11では、制約条件が設定され、続いて、ステップS12では、モデルを構築すべく、制約条件を満たし、かつ、評価結果が最良となるように、モデルのパラメータの値を決定することによって、モデルが演算される。続いて、ステップS13では、この演算されたモデルの精度が予め決められた許容範囲以内であるか否かが判定され、この判定の結果、モデルの精度が許容範囲以内ではない場合(NO)には処理がステップS11に戻され、一方、モデルの精度が許容範囲内である場合(YES)にはステップS14が実行され、この演算されたモデルがシステムのモデルとして出力部3に出力され、出力部3で出力される。
次に、前記ステップS11についてより具体的に説明する。図3において、まず、ステップS21では、入力設定部2の制約条件入力設定部21から、時間応答波形における制約条件(前記時間応答制約条件)が入力され、記憶部4の制約条件記憶部41における時間応答制約条件記憶部411に格納され、記憶される。これによって、後述するように、演算部1の制約条件評価部12が時間応答制約条件記憶部411に記憶されている時間応答制約条件を用いることで、この時間応答制約条件を満たすように、システムのモデルが構築される。
以下に、説明の簡単化のため、1入力1出力系のシステムの場合について、より具体的に説明する。もちろん、以下の説明は、一般性を損なうことなく、多入力多出力系のシステムに拡張可能である。
まず、例えば燃焼炉プラント等の実際のシステムから入出力データ(入力データUdataおよび出力データYdata)が実績入出力データ(実績入力データUdataおよび実績出力データYdata)として取得され、実績データ記憶部46に格納され、記憶されているものとする。システムのモデルが連続した時間と見なされる時間変数tを含む連続時間状態空間モデル式または連続時間伝達関数モデル式であることから、実績入出力データUdata、Ydataは、時間的に連続しているデータに見えるように、より短いサンプリング間隔で取得されることが好ましい。ここでは、表記の簡単化のために、1秒間隔(1秒周期)でサンプリングされているものとする。サンプリング開始の時刻0から時刻tmaxまでのtmax秒間における実績入出力データUdata、Ydataは、時刻t秒の入力値をu(t)とし、時刻t秒の出力値をy(t)とする場合には、それぞれ、次の式1−1、式1−2となる。
Udata=[u(0) u(1) ・・・ u(tmax)]T ・・・(1−1)
Ydata=[y(0) y(1) ・・・ y(tmax)]T ・・・(1−2)
ここで、上付き添え字Tは、いわゆる転置を表す記号である。
時刻0から時刻tEまでのtE秒間における式2の規範入力Urをシステムに加えた場合に得られる出力データに関し、時刻tS1から時刻tE1までの間における上限Ymaxおよび時刻tS2から時刻tE2までの間における下限Yminを式3−1および式3−2とすると、これら式3−1によって記述される上限Ymaxおよび式3−2によって記述される下限Yminが制約条件入力設定部21からモデリング装置Sに与えられる。
Ur=[ur(0) ur(1) ・・・ ur(tE)]T ・・・(2)
Ymax=[ymax(ts1) ymax(ts1+1) ・・・ ymax(tE1)]T ・・・(3−1)
Ymin=[ymin(ts2) ymin(ts2+1) ・・・ ymin(tE2)]T ・・・(3−2)
ここで、ur(t)は、時刻t秒における入力値であり、ymax(t)は、時刻tにおける出力値の上限値であり、そして、ymin(t)は、時刻tにおける出力値の下限値である。また、0≦tS1<tE1≦tEであり、0≦tS2<tE2≦tEである。
なお、規範入力の一具体例として、単位ステップ入力は、Ur=[0 1 1 ・・・ 1]Tであり、単位インパルス入力は、Ur=[0 1 0 ・・・ 0]Tであり、またランプ入力は、Ur=[0 1 2 3 ・・・]Tである。また、規範入力は、1種類であって良く、また複数の種類であっても良い。
続いて、ステップS22では、入力設定部2の制約条件入力設定部21から、定常ゲインにおける制約条件(前記定常ゲイン制約条件)が入力され、記憶部4の制約条件記憶部41における定常ゲイン制約条件記憶部412に格納され、記憶される。例えば、システムの定常ゲインGがGmin≦G≦Gmaxとすると、これら下限値Gminおよび上限値Gmaxが制約条件入力設定部21からモデリング装置Sに与えられる。これによって、後述するように、演算部1の制約条件評価部12が定常ゲイン制約条件記憶部412に記憶されている定常ゲイン制約条件を用いることで、この定常ゲイン制約条件を満たすように、システムのモデルが構築される。
次に、前記ステップS12についてより具体的に説明する。図4において、まず、ステップS31では、入力設定部2のモデル次数同定パラメータ入力設定部22から、モデルの次数および前記モデルの同定パラメータが入力され、記憶部4のモデル次数同定パラメータ記憶部42に格納され、記憶される。
上述の1入力1出力系のシステムの場合において、より具体的に説明すると、例えば、モデリング対象のシステムに対する連続時間状態空間モデル式が式4(式4−1、式4−2)であるとする。
dx(t)/dt=Ax(t)+Bu(t) ・・・(4−1)
ym(t)=Cx(t)+Du(t) ・・・(4−2)
ここで、x(t)は、n次の状態ベクトルであり、ym(t)は、モデルの出力(モデル出力、モデルによる計算値)であり、A、B、CおよびDは、それぞれn行n列、n行1列、1行n列および1行1列の定数行列である。ここでは、説明の簡単化のため、モデリング対象のシステムが厳密にプロパーである、すなわち、D=0とする。
モデルの次数は、このnの値であり、この値nがモデルの次数としてモデル次数同定パラメータ入力設定部22から入力される。この例では、n=2、すなわち、モデルの次数が2である2次のモデルを求めるものとする。また、定数行列A、BおよびCは、例えば可制御正準形式等の所定の形式でよいが、冗長なパラメータが無い方が好ましいことから、この例では、可観測正準形式であるとし、それぞれ式5−1、式5−2および式5−3となる。
ここで、a0、a1、b0およびb1は、モデルを規定する未知パラメータである。これらa0、a1、b0およびb1の未知パラメータが、モデル次数同定パラメータ入力設定部22から入力され、同定パラメータとして指定される。すなわち、モデリング対象のシステムのモデルを構築することは、これら同定パラメータとして指定された未知パラメータa0、a1、b0およびb1の各値を求めることである。
なお、式4(式4−1、式4−2)および式5(式5−1〜式5−3)の求解は、式6の連続時間伝達関数モデル式G(s)を求めること、すなわち、式6の連続時間伝達関数モデル式G(s)における未知パラメータa0、a1、b0およびb1を求めることと等価である。
G(s)=(b1s+b0)/(s2+a1s+a0) ・・・(6)
続いて、ステップS32では、入力設定部2の最適化計算パラメータ入力設定部23から、最適化計算パラメータが入力され、記憶部4の最適化計算パラメータ記憶部43に格納され、記憶される。
上述の1入力1出力系のシステムの場合において、より具体的に説明すると、未知パラメータ(同定パラメータ)a0、a1、b0およびb1の同定は、本実施形態では、例えば、PSO(Particle Swarm Optimization)法が用いられる。このPSO法は、解空間に複数P個の粒子p(粒子群)をランダム(無作為)にばらまき、これら解空間にばらまかれた粒子群が所定の更新則に従って解空間内を移動しながら、最適解を探索する手法である。このPSO法を本実施形態に適用する場合では、解空間は、未知パラメータ(同定パラメータ)a0、a1、b0およびb1を各座標軸とする座標空間(この例では4次元空間)であり、粒子pは、4次元ベクトルで表される前記座標空間上の点となり、更新回数k番目におけるP個のうちのi番目の粒子p(i,k)は、式7によって表現されることになる。
p(i,k)=[a0 a1 b0 b1]T ・・・(7)
以下、適宜、このp(i,k)の[a0 a1 b0 b1]Tをパラメータセットと呼称することとする。
このような最適化手法では、粒子pの個数P(0<i≦P)、更新回数K(0<k≦K)および所定の更新則におけるパラメータ(後述のw、c1、c2)が最適化計算パラメータであり、最適化計算パラメータ入力設定部23から入力される。
続いて、ステップS33では、入力設定部2から評価関数が入力され、設定される。
上述の1入力1出力系のシステムの場合において、より具体的に説明すると、評価関数は、例えば、本実施形態では、実績入力データUdataに対応する実績出力データYdata(=[y(0) y(1) ・・・ y(tmax)]T)と、前記実績入力データUdataをモデル(上述の例では、式4(式4−1、式4−2)および式5(式5−1〜式5−3)または式6)に入力した場合にこのモデルから出力されるモデル出力データYi(=[yi(0) yi(2) ・・・ yi(tmax)]T)との誤差を評価する関数である。このような評価関数Jiは、例えば、式8−1によって与えられる。
Ji=||W・(Yi−Ydata)|| ・・・(8−1)
ここで、||v||は、いわゆる、ベクトルvのユークリッドノルムを表し、vが縦ベクトルである場合では、||v||=√(vTv)である。また、Wは、出力の予測誤差に対していずれの部分を重要視して評価するかを定める重み行列であり、例えば対角行列とする。
また例えば、評価関数Jiは、例えば、式8−2によって与えられる無限ノルムである。
Ji=||W・(Yi−Ydata)||∞ ・・・(8−2)
ここで、||v||∞は、いわゆる、ベクトルvの無限ノルム(最大値ノルム)を表し、v=[v1 v2 ・・・ vn]Tである場合では、||v||∞=max(|v1|、|v2|、・・・、|vn|)である。
また例えば、評価関数Jiは、例えば、式8−3によって与えられる。
Ji=||W・(Yi−Ydata)||1 ・・・(8−3)
ここで、||v||1は、いわゆる、ベクトルvの1ノルムを表し、v=[v1 v2 ・・・ vn]Tである場合では、||v||1=Σ|vi|である。ただし、Σは、i=1からi=nまでの和である。
評価関数記憶部44には、例えば、このような複数の評価関数Ji(式8−1〜式8−3)が記憶されており、これら複数の評価関数Jiの中から、評価関数入力設定部24を介して1つの評価関数Jiが選定され、入力設定される。
なお、このような評価関数Jが、実績出力データとモデル出力データとの誤差に重み付けを行う重みWを有し、この重みWが、その値が大きいほど評価関数Jによる評価結果に大きく影響を与えるものである場合に、この重みWの値は、所定の閾値以上のノイズを含む実績出力データに対応する区間に対し、前記区間を除く他の区間よりも相対的に小さくなるように設定されてよい。このような前記区間およびその重みWならびに前記他の区間およびその重みWは、例えば、入力設定部2の最適化計算パラメータ入力設定部23から入力設定される。このように構成することによって、所定の閾値以上のノイズを含む実績出力データの、評価結果に与える影響が低減され、システムのモデルをより高精度に構築することができる。
続いて、ステップS34では、演算部1のモデル出力演算部11によって、モデル出力が演算され、記憶部4のモデル出力等記憶部45に記憶される。
上述の1入力1出力系のシステムの場合において、より具体的に説明すると、まず、P個の粒子p(i,0)(更新回数k=0、i=0〜P)、例えば100個の粒子p(i,0)(i=0〜100)の初期値が例えば乱数を生成することによってランダムに設定される。これによってP個の各パラメータセットp(i,0)の値が設定される。そして、i番目の粒子p(i,0)の移動速度v(i,0)の初期値が例えばv(i,0)=[0 0 0 0]Tに設定される。
ここで、制約条件として、モデリング対象のシステムが安定な系であることが、先験的知識として判明している場合には、a0>0かつa1>0の範囲で、i番目の粒子p(i,0)の初期値がランダムに設定されることが好ましい。
また、制約条件として、モデリング対象のシステムにおける定常ゲインが正の値であることが、先験的知識として判明している場合には、b0の符号とa0の符号が等しくなるように、すなわち、b0×a0>0となるように、i番目の粒子p(i,0)の初期値がランダムに設定されることが好ましい。
また、制約条件として、モデリング対象のシステムに不安定なゼロ点がないことが、先験的知識として判明している場合には、b1の符号とb0の符号が等しくなるように、すなわち、b1×b0>0となるように、i番目の粒子p(i,0)の初期値がランダムに設定されることが好ましい。
また、制約条件として、モデリング対象のシステムにおける定常ゲインがGmin以上かつGmax以下であることが、先験的知識として判明している場合には、a0×Gmin≦b0≦a0×Gmaxを満たすように、i番目の粒子p(i,0)の初期値がランダムに設定されることが好ましい。
このようにi番目の粒子p(i,0)の初期値を設定する際に、制約条件を加味することによって、解が存在しない領域を探索する無駄を省くことが可能となる。
このようにP個の粒子pの値をそれぞれ設定することで、1個の粒子pに対応して1個のモデルが規定される。すなわち、i番目の粒子p(i,0)に対応して1個のモデルが規定される。上述の例では、100個の各粒子pのそれぞれに対応して100個のモデルが規定されることになる。
このi番目の粒子p(i,0)に対応して規定されるモデル(この例では式4および式5によって表現されるモデル、あるいは、式6によって表現されるモデル)に対し、実績データ記憶部46に記憶されている実績入力データUdataが代入され、モデル出力データYiが得られ、この得られたモデル出力データYiが記憶部4のモデル出力等記憶部45にi番目の粒子p(i,0)に対応付けて記憶される。このモデル出力データYiの計算には、公知の制御系シミュレーションソフトウェア、例えば、The Math Work社のMATLAB/Simulink等が用いられる。このようにP個の各パラメータセットp(i,k)のそれぞれに対し、モデル出力が演算され、記憶部4のモデル出力等記憶部45にパラメータセットp(i,k)に対応付けて記憶される。
続いて、ステップS35では、制約条件評価部12によって、制約条件が評価され、この評価結果が記憶部4のモデル出力等記憶部45にパラメータセットp(i,k)に対応付けて記憶される。
上述の1入力1出力系のシステムの場合において、より具体的に説明すると、ステップS34と同様に、このi番目の粒子p(i,0)に対応して規定されるモデルに対し、所定の規範入力Urが代入され、前記所定の規範入力Urに対応するモデル出力データYri(=[yri(0) yri(1) ・・・ yri(tmax)]T)が得られる。そして、この得られた前記所定の規範入力Urに対応するモデル出力データYriが規範入力に対する制約条件を満たしているか否かが判断される。例えば、制約条件として、時刻tS1から時刻tE1までの間における上限値Ymaxおよび時刻tS2から時刻tE2までの間における下限値Yminが設定されている場合には、この前記所定の規範入力Urに対応するモデル出力データYriが時刻tS1から時刻tE1までの間で上限Ymax以下で、かつ、時刻tS2から時刻tE2までの間で下限Ymin以上であるか否かが判断される。この判断の結果、この前記所定の規範入力Urに対応するモデル出力データYriが時刻tS1から時刻tE1までの間で上限Ymax以下で、かつ、時刻tS2から時刻tE2までの間で下限Ymin以上ではない場合には、この前記所定の規範入力Urに対応するモデル出力データYriが制約条件を満たさないものとされ、この制約条件の評価結果が記憶部4のモデル出力等記憶部45にi番目の粒子p(i,0)に対応付けて記憶される。このようにP個の各パラメータセットp(i,k)のそれぞれに対し、制約条件が評価され、この評価結果が記憶部4のモデル出力等記憶部45にパラメータセットp(i,k)に対応付けて記憶される。
なお、上述では、粒子p(i,0)の初期値の設定の際に、定常ゲインGの制約条件を考慮したが、定常ゲインGの下限値Gminおよび上限値Gmaxが制約条件として設定されている場合に、実績入力データUdataに対応するモデル出力データYiが下限値Gmin以上であって、かつ、上限値Gmax以下であるか否かが判断されることで、定常ゲインGの制約条件が評価されても良い。
続いて、ステップS36では、評価演算部13によって、評価関数による評価値が演算され、記憶部4のモデル出力等記憶部45に記憶される。
上述の1入力1出力系のシステムの場合において、より具体的に説明すると、i番目の粒子p(i,0)に対応して規定されるモデルに対し、前記評価関数Jiによる評価値が演算され、この評価関数Jiの評価値が記憶部4のモデル出力等記憶部45にi番目の粒子p(i,0)に対応付けて記憶される。このようにP個の各パラメータセットp(i,k)のそれぞれに対し、評価関数による評価値が演算され、記憶部4のモデル出力等記憶部45にパラメータセットp(i,k)に対応付けて記憶される。
ここで、上記制約条件の評価結果を記憶部4のモデル出力等記憶部45にi番目の粒子p(i,0)に対応付けて記憶する際に、上述のように、この制約条件の評価結果そのものを記憶部4のモデル出力等記憶部45にi番目の粒子p(i,0)に対応付けて記憶してもよいが、本実施形態では、この制約条件の評価結果が評価関数Jiの評価値で記憶部4のモデル出力等記憶部45にi番目の粒子p(i,0)に対応付けて記憶される。例えば、制約条件を満たさない場合に、評価が悪いことを示す評価値が記憶部4のモデル出力等記憶部45にi番目の粒子p(i,0)に対応付けて記憶される。より具体的には、例えば、評価関数Jiが、評価が高いほどその評価値が小さくなるように構成されている場合では、制約条件を満たさない場合に、例えばモデル出力データYiがその値として取り得ない大きな評価値が、記憶部4のモデル出力等記憶部45にi番目の粒子p(i,0)に対応付けて記憶される。また例えば、評価関数Jiが、その取り得る値の範囲が負であるように構成されている場合では、制約条件を満たさない場合に、正の値(例えば固定値や制約条件の外れ度合いに従った値等)が記憶部4のモデル出力等記憶部45にi番目の粒子p(i,0)に対応付けて記憶される。このような評価関数Jiは、例えば、式8−4によって与えられる。
Ji=tan−1(||W・(Yi−Ydata)||)−π/2 ・・・(8−4)
続いて、ステップS37では、モデル更新決定部14によって、最適化方法の繰り返し計算における更新が終了しているか否かが判断され、この判断の結果、更新が終了していない場合(NO)では、ステップS38の処理が実行され、一方、更新が終了している場合(YES)では、ステップS39の処理が実行される。更新回数Kは、任意の回数で良いが、PSO法では、解の精度や演算処理量等を考慮することによって適宜な回数に設定され、例えば、K=100回とされる。
ステップS38では、モデル更新決定部14によって、所定の更新則に従って更新され、処理がステップS34に戻される。
上述の1入力1出力系のシステムの場合において、より具体的に説明すると、k回目(k=0,1,2,・・・、K)の更新時点において、P個の各粒子pの中からその評価関数Jの値が最も小さい粒子pが探索され、第1最良評価粒子p(best,k)とされる。例えば、k回目の更新時点において、q番目の粒子p(q、k)の評価関数Jqの値が最も小さい場合には、p(best,k)=p(q,k)である。
また、i番目の粒子pにおいて、0回目(初期値)からk−1回目までの中で評価関数Jの値が最も小さい粒子pが探索され、第2最良評価粒子p(i,best)とされる。例えば、r回目の更新時点において、i番目の粒子p(i、r)の評価関数Jiの値が最も小さい場合には、p(i,best)=p(i,r)である。なお、k=0の場合では、それぞれの初期値そのものが各粒子pのp(i,best)である。
そして、例えば、式9(式9−1、式9−2)によって与えられる更新則に従って、各粒子pの速度および位置がそれぞれ決定され、更新される。
v(i,k+1)=wv(i,k)
+c1R1ik(p(i,best)−p(i,k))
+c2R2ik(p(best,k)−p(i,k)) ・・・(9−1)
p(i,k+1)=p(i,k)+v(i,k+1) ・・・(9−2)
ここで、w、c1、c2は、それぞれ、調整パラメータであり、最適化計算パラメータとして上述したように最適化計算パラメータ入力設定部23から入力され、設定される。例えば、w=0.9、c1=0.8およびc2=0.8等である。また、R1ikおよびR2ikは、それぞれ対角行列であり、対角要素のそれぞれは、1〜0の範囲で一様に分布する乱数である。この乱数値は、固定値ではなく、各粒子pごとに、かつ、各更新kごとに、それぞれランダムな値が使用される。
ここで、制約条件として、モデリング対象のシステムが安定な系であることが、先験的知識として判明している場合には、a0>0かつa1>0の範囲に入るように更新することが好ましい。より具体的には、更新後の粒子pがa0<0となった場合には、式10−1に示すように、更新前のa0と0とを内分する点に更新される。すなわち、0〜1の乱数αにおいて、
a0(i,k+1)=(1−α)・a0(i,k)+α・0 ・・・(10−1)
ここで、i番目の粒子pについて、a0(i,k+1)は、更新後のa0であり、a0(i,k)は、更新前のa0である。また、粒子pの速度vは、a0に対応する部分v(i,k+1){a0}について、式10−2となる。
v(i,k+1){a0}=a0(k+1)−a0(k) ・・・(10−2)
また、制約条件として、モデリング対象のシステムにおける定常ゲインが正の値であることが、先験的知識として判明している場合には、b0の符号とa0の符号が等しくなるように、すなわち、b0×a0>0となるように、更新することが好ましい。
また、制約条件として、モデリング対象のシステムに不安定なゼロ点がないことが、先験的知識として判明している場合には、b1の符号とb0の符号が等しくなるように、すなわち、b1×b0>0となるように、更新することが好ましい。
また、制約条件として、モデリング対象のシステムにおける定常ゲインがGmin以上かつGmax以下の範囲であることが、先験的知識として判明している場合には、a0×Gmin≦b0≦a0×Gmaxを満たすように、更新することが好ましい。
このように粒子pを更新する際に、制約条件を加味することによって、解が存在しない領域を探索する無駄を省くことが可能となり、よりよい解をより早く探索することができる。
さらに、制約条件として、同定パラメータの符号a0、a1、b0、b01が先験的知識として判明している場合には、これら同定パラメータが指数形式で表現されても良い。例えば、式6の伝達関数が安定かつ不安定零点なしかつゲインが負であることが先験的知識として判明している場合では、a0の符号およびa1の符号が正であって、かつ、b0の符号およびb1の符号が負でなければならないことから、a0=exp(λ0)、a1=exp(λ1)、b0=−exp(μ0)、b1=exp(μ1)とおき、式7のパラメータセットがp(i,k)=[λ0 λ1 μ0 μ1]Tとなる。したがって、同定パラメータ[a0 a1 b0 b1]Tの座標空間ではなく、このパラメータ[λ0 λ1 μ0 μ1]Tの座標空間で解が探索される。このように構成することによって、−∞<λ0<+∞、−∞<λ1<+∞、−∞<μ0<+∞、−∞<μ1<+∞の実数全体での演算処理が可能となり、符号や範囲をチェックすることなく演算処理が可能となる。
このようにPSO法による更新を行う場合に、定常ゲインの符号、定常ゲインの範囲、前記システムの安定性に関する情報および前記システムの不安定なゼロ点に関する情報のうちの少なくとも1つが用いられてよい。このようにモデリング装置Sを構成することによって、不必要な探索を行うことなく、これらのうちの少なくとも1つを制約条件として満たすシステムのモデルを構築することができる。
こうしてP個の各パラメータセットp(i,k)が所定の更新則に従って更新され、処理がステップS34に戻される。
一方、ステップS39では、モデル更新決定部14によって、制約条件を満たし、かつ最良の評価値を与えるように、同定パラメータセットの値が決定され、このシステムのモデルが構築される。例えば、上述のPSO法では、P×K個の各パラメータセットp(i,k)の中から、制約条件を満たし、かつ最良の評価値を与えるパラメータセットp(i,k)が選択され、この選択されたパラメータセットp(i,k)の値が、モデリング対象のシステムのモデルを規定する同定パラメータの値として決定され、このシステムのモデルが決定される。
このように本実施形態におけるモデリング装置Sおよびモデリング方法では、実績出力データとモデル出力データとの誤差が評価関数Jで評価され、この評価結果が最良となるようにモデルのパラメータの値が求められ、モデルが構築される。そして、先験的知識を表す先験的情報が制約条件として設定され、このモデルを構築する際に、制約条件も満たすように、モデルのパラメータの値が求められる。このため、本実施形態におけるモデリング装置Sおよびモデリング方法は、上述の非特許文献1や特許文献1とは異なる手法で、先験的知識を利用することによって、前記質および量の不充分な入出力データであっても、より精度の高いモデルを構築することができる。
また、本実施形態におけるモデリング装置Sおよびモデリング方法では、システムにおける、定常ゲインGに関する上限値Gmaxおよび下限値Gminのうちの少なくとも一方が、本実施形態では定常ゲインGに関する上限値Gmaxおよび下限値Gminが制約条件とされている。このため,本実施形態におけるモデリング装置Sおよびモデリング方法は、定常ゲインGに関する先験的情報も制約条件として加味され、定常ゲインGに関する先験的情報を満たすシステムのモデルを構築することができる。
次に、一実施例として、実施形態におけるモデリング装置Sを用いて例えばごみ焼却炉(ガス化溶融炉)のモデルを作成した場合について、説明する。
図5は、ゴミ焼却炉において、給塵機の速度をステップ状に変化させた場合における蒸気発生量の変化を示す図である。図5(A)は、給塵機の速度を示し、その横軸は、秒単位で表す経過時間を示し、その縦軸は、速度を示す。図5(B)は、蒸気発生量を示し、その横軸は、秒単位で表す経過時間を示し、その縦軸は、蒸気発生量を示す。図6は、制約条件としてのステップ応答の上限および下限を示す図である。図6の横軸は、秒単位で表す経過時間を示し、その縦軸は、蒸気発生量を示す。図7は、実施形態のモデリング装置でモデリングされたゴミ焼却炉のモデルにおけるステップ応答を説明するための図である。図7(A)は、給塵機の速度を示し、その太い実線は、モデル出力としての蒸気発生量を示し、その細い実線は、図5(B)に示す実績データとしての蒸気発生量を示す。図7(B)は、モデルのステップ応答(実線)、時刻tS1から時刻tE1までの間における上限のステップ応答Ymax(破線)および時刻tS2から時刻tE2までの間における下限のステップ応答Ymin(破線)が示されている。図8は、比較例としての従来の手法によりモデリングされたゴミ焼却炉のモデルにおけるステップ応答を説明するための図である。図8(A)は、給塵機の速度を示し、その太い実線は、モデル出力としての蒸気発生量を示し、その細い実線は、図5(B)に示す実績データとしての蒸気発生量を示す。図8(B)は、モデルのステップ応答(実線)、時刻tS1から時刻tE1までの間における上限のステップ応答Ymax(破線)および時刻tS2から時刻tE2までの間における下限のステップ応答Ymin(破線)が示されている。図9は、定常ゲインの制約条件を考慮した場合であって、実施形態のモデリング装置でモデリングされたゴミ焼却炉のモデルにおけるステップ応答を説明するための図である。図9(A)は、給塵機の速度を示し、その太い実線は、モデル出力としての蒸気発生量を示し、その細い実線は、図5(B)に示す実績データとしての蒸気発生量を示す。図9(B)は、モデルのステップ応答(実線)、時刻tS1から時刻tE1までの間における上限のステップ応答Ymax(破線)および時刻tS2から時刻tE2までの間における下限のステップ応答Ymin(破線)が示されている。
ゴミ焼却炉では、ゴミを焼却炉に供給する給塵機を入力とし、焼却炉で発生する蒸気発生量を出力とすると、図5(A)に示すように、給塵機の速度がステップ状に変化させた場合に、蒸気発生量は、図5(B)に示すように、基本的に、給塵機の稼動開始に伴って徐々に増大するとともに給塵機の稼動停止に伴って徐々に減少するプロファイルであるが、このプロファイルに外乱の影響によりノイズが重畳しており、その波形が大きく乱れている。この原因として、例えば、給塵機の速度が一定であっても給塵機によって実際に焼却炉に投入されるゴミの量は、必ずしも一定ではなく、変動し、蒸気発生量の外乱となると考えられる。しかしながら、実際に稼働中のゴミ焼却炉(生産設備やプラント等の制御対象)では、ゴミ焼却炉のモデル化のために実績入出力データを取得する実験機会や、ゴミ焼却炉に投入可能なゴミ量等が限定され、このようなノイズを含むデータしか取得できない場合もある。
このようなゴミ焼却炉において、先験的情報として、時刻0において給塵機の速度を1単位だけステップ状に変化させる単位ステップ入力で変化させた場合に、蒸気発生量(ステップ応答)は、図6に示すように、時刻tS1から時刻tE1までの間における上限のステップ応答Ymax(破線、ステップ応答上限)と、時刻tS2から時刻tE2までの間における下限のステップ応答Ymin(破線、ステップ応答下限)との間に、存在することが分かっているものとする。このような先験的知識は、通常、ゴミ焼却炉を実際に稼動させるオペレータの経験や、ゴミ焼却炉の過去実績データ等が理解されるものであるが、ゴミ焼却炉の設計値等からも理解される場合もある。例えば、ゴミのカロリー(熱量)における変動範囲と、給塵機を1単位だけ変化させた場合におけるゴミの供給量変化に基づき、蒸気発生量の変化は、概略計算が可能である。
なお、このような推定計算は、このゴミ焼却炉のようなプラントだけでなく、種々の装置で可能である。例えば、リンク機構として、マニピュレータでは、入力トルクを1単位だけ変化させた場合に、アームの回転速度がどのくらいの時間でどのくらいの範囲で変化するか、その慣性や粘性等のパラメータから、時間応答の上限および下限が概略計算することができる。
図5に示す実績データおよび図6に示す制約条件の下に、ゴミ燃焼炉を本実施形態のモデリング装置Sによってモデル化することによって作成されたゴミ燃焼炉のモデルにおけるステップ応答を図7に示す。この本実施形態のモデリング装置Sによるゴミ燃焼炉のモデルでは、そのステップ応答(蒸気発生量)は、図7(A)に示すように、図5(B)に示す実際の蒸気発生量の変化によく一致(フィッティング)しており、そして、図7(B)に示すように、ステップ応答の上限および下限の制約条件を満たしている。
一方、図5に示す実績データおよび図6に示す制約条件の下に、ゴミ燃焼炉を従来の手法によってモデル化することによって作成されたゴミ燃焼炉のモデルにおけるステップ応答を図8に示す。この従来の手法によるゴミ燃焼炉のモデルでは、図7と比較すると、図7に示すステップ応答(蒸気発生量)ほど、実際の蒸気発生量の変化に一致していない。特に、図8(A)に破線の○で囲んだ外乱の大きな部分の実績出力データにモデルが合わせ込もうとされており、そのため、図8に示すステップ応答では、図7に示すステップ応答に較べて緩やかな立ち上がりとなってしまってしい、図7に示すステップ応答に較べてフィッティングの精度が落ちている。さらに、図8(B)に示すように、ステップ応答の上限および下限の制約条件を満たしていない部分がある。この結果、従来の手法では、再度、モデル化に使用するデータ区間を変更したり、フィルタリング等のデータの前処理を行ったり、あるいは、実績入出力データを取り直ししたり等、先験的情報の制約条件を満たすモデルが作成されるまで、試行錯誤的にモデル化作業を繰り返す必要がある。
しかしながら、本実施形態のモデリング装置Sは、図7と図8とを比較すると分かるように、このようなモデル化作業が軽減される。
ここで、上述のゴミ燃焼炉のモデル化では、図6(B)に示すように、ステップ応答が1200秒で略整定(飽和)していることから、制約条件の時間範囲をこの1200秒で区切っている。このため、1200秒を越える部分で、モデルの実用上問題の無いレベルであるが、ステップ応答がステップ応答の上限を越えてしまっている。そこで、厳密に、ステップ応答がステップ応答の上限および下限を逸脱しないようにするために、定常ゲインGの上下限値が制約条件とされてもよい。すなわち、上述したように、制約条件として、モデリング対象のシステムにおける定常ゲインがGmin以上かつGmax以下であることが、先験的知識として判明している場合に、a0×Gmin≦b0≦a0×Gmaxを満たすように、i番目の粒子p(i,0)の初期値がランダムに設定され、そして、a0×Gmin≦b0≦a0×Gmaxを満たすように、更新される。このような制約条件を課すことにより、ステップ応答(蒸気発生量)は、図9(A)に示すように、図5(B)に示す実際の蒸気発生量の変化によく一致(フィッティング)しているだけでなく、図9(B)に示すように、1200秒を越えた範囲も含めて、ステップ応答の上限および下限の制約条件を満たしている。
なお、上述では、給塵機の速度を入力とし、蒸気発生量を出力としたが、入出力の関係は、これに限定されるものではなく、種々の関係であって良い。例えば、入力としては、押込空気量や、空気供給量や、バーナ油量や、蒸気弁開度や、触媒投入量等が挙げられ、出力としては、砂層温度や、炉内温度や、ボイラドラム圧力や、炉内圧力や、酸素濃度や、一酸化炭素濃度や、NOx濃度等が挙げられる。
本実施形態のモデリング装置Sの適用例としてゴミ焼却炉のモデル化について説明したが、本実施形態のモデリング装置Sは、操作量を入力とするとともに制御量を出力とすることで、様々な制御対象のモデル化に適用可能である。本実施形態のモデリング装置Sは、例えば、圧延機のモデル化に適用可能であり、例えば、圧下位置やロール速度等を入力とし、スタンド間張力や板幅・板厚等を出力とすることができる。また例えば、本実施形態のモデリング装置Sは、マニピュレータのモデル化に適用可能であり、トルク等を入力とし、アーム速度や位置等を出力とすることができる。
なお、上述の実施形態において、モデルのパラメータにおける初期値を設定してから前記モデルのパラメータの値を求める最適化方法によって前記モデルのパラメータの値を求める場合に、前記初期値自体を前記モデルのパラメータとして扱うように、モデリング装置Sが構成されてもよい。一般に、モデルのパラメータにおける初期値を設定してから前記モデルのパラメータの値を求める最適化方法では、この最適化手法によって構築されるシステムのモデルの精度は、初期値に左右されることが多い。このため、このように構成されることによって、この初期値自体も未知パラメータとして取り扱われるので、初期値に左右されることなく、システムのモデルをより高精度に構築することが可能となる。
例えば、上述の1入力1出力系のシステムの場合において、初期状態が0(ゼロベクトル)であることを前提にモデルが構築されている。すなわち、モデリング対象のシステムが整定している状態(平衡状態)において入力データが印加されたものと見なされており、モデルの構築に用いた例えば実績入出力データ等の入出力データは、初期状態が0であるとされている。
式4において、状態x(t)の挙動は、一般に、式11となることが知られている。
x(t)=Φ(t)x(0)+∫Φ(t−τ)Bu(τ)dτ ・・・(11)
ここで、積分区間は、0からtまでであり、Φ(t)は、遷移マトリックスと呼ばれ、式12によって表される。
Φ(t)=exp(At)
=I+At+A2t2/2!+A3t3/3!+・・・ ・・・(12)
この式11で表される状態x(t)を式4−2に代入すると、式13が得られ、初期状態x(0)=0でなければ、出力ym(t)は、初期状態x(0)の影響を受けることが分かる。
ym(t)=CΦ(t)x(0)+C∫Φ(t−τ)Bu(τ)dτ+Du(t) ・・・(13)
モデルを構築するための入出力データを得る実際の実験において、例えばマニピュレータ等の機械系のシステムでは、システムが静止した状態から実験を開始することが可能であり、この場合では初期状態を0と見なせるが、例えばゴミ燃焼炉や化学プラント等のプラント系のシステムでは、稼動中に実験を行わなければならず、状態が絶えず変動しており、通常、初期状態が0とはならない。
このため、上述のように、初期状態自体もパラメータとして扱い、モデルを構築することが好ましい。
より具体的には、例えば、上述の1入力1出力系のシステムの場合において、式4および式5、あるいは、式6では、モデリング対象のシステムのモデル次数が2次であるから、初期状態x(0)=[x1 x2]Tとおくことができる。このため、P個のうちのi番目の粒子p(i,k)は、式7に代えて、式14によって表現されることになる。
p(i,k)=[a0 a1 b0 b1 x1 x2]T ・・・(14)
また、この場合におけるi番目の粒子pの移動速度v(i,0)の初期値は、v(i,0)=[0 0 0 0 0 0]Tに設定される。
このような初期値の設定を行って、以下、上述と同様の処理によってシステムのモデルが構築される。
また、上述の実施形態において、予め決められた所定の最適化方法によって前記モデルのパラメータの値を求め、この求めた前記モデルのパラメータの値をPSO法の初期値に加えて用いるように、モデリング装置Sが構成されてもよい。前記所定の最適化方法として、PSO法、最小二乗法および部分空間法等の各種の最適化手法が挙げられる。なお、最小二乗法や部分空間法が用いられる場合には、離散時間で一旦状態空間モデル式や伝達関数モデル式を求め、これを連続時間のモデルに変換される。また、PSO法を用いる場合では、評価関数のみを単に最小化するモデルが求められても良いし、その際に、Ji=||W(Yi−Ydata)||<ε(εは所定の閾値である)となる粒子pを保存しておいても良い。このPSO法を用いる場合では、初期値の決定およびモデルの構築の2回、PSO法が使われることになる。一般に、モデルのパラメータにおける初期値を設定してから前記モデルのパラメータの値を求める最適化方法では、この最適化手法によって構築されるシステムのモデルの精度は、初期値に左右されることが多い。このため、PSO法による場合でもランダムに振る範囲が試行錯誤によって決定される。このため、この試行錯誤の手間を軽減し、最適なモデルを求めるためには、最初のモデルが入出力データにある程度フィッティングしていることが望ましい。この構成では、PSOに用いる初期値として、予め決められた所定の最適化方法によって求められたモデルのパラメータの値も用いるので、より真に近い値を含む初期値からPSO法の演算処理を開始することができるから、不必要な探索を行うことなく、システムのモデルをより高精度に構築することができる。
本発明を表現するために、上述において図面を参照しながら実施形態を通して本発明を適切且つ十分に説明したが、当業者であれば上述の実施形態を変更および/または改良することは容易に為し得ることであると認識すべきである。したがって、当業者が実施する変更形態または改良形態が、請求の範囲に記載された請求項の権利範囲を離脱するレベルのものでない限り、当該変更形態または当該改良形態は、当該請求項の権利範囲に包括されると解釈される。