JP2006195543A - モデル同定装置およびモデル同定プログラム - Google Patents

Abstract

【課題】 データ数や数値演算量の増大を伴うことなく、周波数帯域を分割しながら周波数応答特性を計算することが可能となるとともに、間引きを行うためのサンプリング周波数と周波数帯域の分割位置を客観的に一意に設定した上で、分割して得られた周波数応答特性を結合させたり、モードの数を考慮した帯域分割を行ったりすることなく、伝達関数を同定する。
【解決手段】 サンプリング手段１３は、複数の入力信号ｕに対する応答信号ｙの組に対し、データ数がすべての組に対して等しくなるように各組ごとに異なるサンプリング周期で間引いてサンプリングし、パルス伝達関数同定手段１６は、複数個の周波数応答特性に基づいて、そのすべての周波数応答特性を表すことができる測定時のサンプリング周期での１つのパルス伝達関数を同定する。
【選択図】 図２

Description

本発明はモデル同定装置およびモデル同定プログラムに関し、特に、特性不明システムに対する入力信号と応答信号とから、システムの特性を示すモデルを同定する方法に適用して好適なものである。

従来の特性不明システムのモデル同定では、そのシステムの周波数応答特性を測定し、その周波数応答特性に基づいてシステムの伝達関数を推定する方法が知られている。ここで、周波数応答特性を測定する方法としては、掃引正弦波を用いた周波数応答法（非特許文献１）が知られている。一方、伝達関数を推定する方法としては、偏分反復法（非特許文献２、特許文献１）やモード円適合法（非特許文献３）が知られている。

しかし、広域に渡る周波数応答特性を精度よく測定するためには、測定用信号を長時間に渡ってシステムに入力する必要があり、システムに負担がかかる。このため、Ｍ系列信号などの擬似白色信号を、同定対象の周波数帯域下限周波数で決まる期間だけ入力し、その入力信号とそれに対する応答信号とから最小二乗法を用いて差分方程式を同定して、その差分方程式に対する周波数応答から伝達関数を推定する方法が開示されている（特許文献２）。

そして、擬似白色信号を用いることで広域に渡る周波数応答特性を精度よく測定するために、ある一定のサンプリング周期で測定した入出力信号を、各々複数の低域通過フィルタにより複数の周波数帯域成分に分割し、その帯域分割された複数組の入出力信号を各々の周波数帯域に対応したサンプリング周期で間引いて再度サンプリングし、その間引かれたデータから最小二乗法を用いて差分方程式を複数個同定し、その複数個の差分方程式に対して周波数応答を計算し、得られた複数個の周波数応答特性を結合させ、その結合させた周波数応答から伝達関数を推定する方法が開示されている（特許文献２、特許文献３）。

ここで、低域通過フィルタの通過周波数帯域の決定に関しては、同定に用いる応答信号とは別の同定対象のデータに対するスペクトル解析結果から決定するか（特許文献３）、各通過周波数帯域の上限を同定対象の各反共振周波数と一致させる方法が用いられている（非特許文献３）。
また、伝達関数を推定する時においても、周波数応答特性を精度よく同定するために周波数帯域を分割する方法が用いられており、分割した帯域ごとにモード特性を同定する方法が知られている。ここで、偏分反復法において、分割した周波数帯域の各帯域内に含まれるモードの数を限定し、各帯域内の伝達関数を同定する方法が開示されている（特許文献１）。また、モード円適合法においては、モードが帯域内に１つだけ存在するように周波数帯域を分割し、分割した帯域ごとに２次伝達関数を同定する方法が開示されている（非特許文献３）。

木村英紀著、「ディジタル信号処理と制御」、株式会社昭栄堂、昭和６２年９月、Ｐ２３３−２３６ 安田仁彦著、「モード解析と動的設計」、株式会社コロナ社、昭１９９３年１１月、Ｐ１８２−１８４ 橋本誠司、石川赳夫、足立修一、船渡寛人、神山健三、磯島彰宏氏著、「マルチデシメーション同定法を用いたロバスト制御のための柔軟構造物のモデリング手法」、電気学会論文誌Ｄ部門、１２４巻５号、２００４年、Ｐ４７１−４７８６ 特開平１０−３８７４８号公報 特開平３−２８２７１８号公報 特開平３−２１７９０１号公報

しかしながら、システムの広域に渡る周波数応答特性を示す伝達関数を精度よく同定するために、高いサンプリング周波数で低い周波数帯域を同定すると、扱うデータ数が増大するという問題があった。例えば、１Ｈｚ〜１０００Ｈｚまでの周波数帯域を同定するとして、サンプリング周波数を４０００Ｈｚで測定すると、入出力信号それぞれ４０００個のデータを測定する必要がある。さらに、これらのデータを用いて最小二乗法にて差分方程式を同定する場合、４０００個の方程式に対する最小二乗解法を行わなければならず、計算量が膨大になる。

一方、同定に用いる入出力信号を各々複数の低域通過フィルタにより複数の周波数帯域成分に分割して各周波数帯域の周波数応答特性を計算する方法では、分割される周波数帯域の数に比例してデータ量が増えるとともに、最小二乗法による差分方程式の導出時に、データ量の増加に比例して数値演算量が増えるという問題があった。
例えば、表１に示すように、１Ｈｚ〜１０００Ｈｚまでの周波数帯域を１Ｈｚ〜２５０Ｈｚ、２５０Ｈｚ〜５００Ｈｚ、５００Ｈｚ〜７５０Ｈｚ、７５０Ｈｚ〜１０００Ｈｚの４個に分割した場合、各周波数帯域での入力信号のデータ数および最小二乗解法の方程式の数はそれぞれ１０００、２０００、３０００、４０００になる。このため、１Ｈｚ〜１０００Ｈｚまでの周波数帯域を４個に分割すると、合計で１００００個のデータを用意する必要があるだけでなく、最小二乗解法の方程式の数も合計で１００００個となり、データ数および最小二乗解法の方程式の数は、周波数帯域を分割しない場合に比べて２．５倍に増加する。

また、周波数帯域ごとに得られた周波数応答特性を結合させる際には、周波数帯域の上限周波数と下限周波数付近は間引きサンプリングによる影響を考慮しなければならないが、周波数帯域の上限周波数と下限周波数付近では互いに隣り合う周波数帯域が重なるように設定し、重なり合う部分の平均値を求める方法が取られる。このため、周波数帯域の結合後の精度は、周波数帯域の上限周波数と下限周波数の取り方に依存するという問題があった。

また、低域通過フィルタの通過周波数帯域の決定に関し、同定に用いる応答信号とは別の同定対象のデータに対するスペクトル解析結果から決定する方法（特許文献３）では、同定に用いる応答信号とは別個に信号を測定する必要があるだけでなく、その信号に対するスペクトル解析を別途行う必要があり、作業工程の増大を招くという問題があった。また、各通過周波数帯域の上限を同定対象の各反共振周波数と一致させる方法（非特許文献３）では、特性不明システムの同定を行う際に、反共振点がどの周波数点にあるか不明であるため、反共振点を容易に設定することができないという問題があった。また、擬似白色信号の周波数応答から反共振点を探し出す方法もあるが、図３に示すような周波数応答解析結果から視覚的に探索する必要があるため、恣意的な探索となり、同定以前に作業工程が増加するという問題があった。

さらに、伝達関数を推定する場合にも、周波数応答特性を同定する時と同様な問題がある。例えば、偏分反復法では、偏分反復法の演算の前にモードの数を考慮した帯域分割のための工程が必要となる。また、モード円適合法では、モードが帯域内に１つだけ存在するように周波数帯域を分割しなければならず、特性不明システムの同定を行うために周波数応答特性を測定する際に、モードの数やモードの周波数帯域を精度よく解析しなければならないという問題があった。

そこで、本発明の目的は、データ数や数値演算量の増大を伴うことなく、周波数帯域を分割しながら周波数応答特性を計算することが可能となるとともに、間引きを行うためのサンプリング周波数と周波数帯域の分割位置を客観的に一意に設定した上で、分割して得られた周波数応答特性を結合させたり、モードの数を考慮した帯域分割を行ったりすることなく、伝達関数を同定することが可能なモデル同定装置およびモデル同定プログラムを提供することである。

上述した課題を解決するために、請求項１記載のモデル同定装置によれば、一定のサンプリング周期で期間の異なる複数の入力信号を入力する入力信号入力手段と、前記入力信号に対する応答信号を前記サンプリング周期で測定する応答信号測定手段と、前記複数の入力信号に対する応答信号の組に対し、データ数がすべての組に対して等しくなるように各組ごとに異なるサンプリング周期で間引いてサンプリングするサンプリング手段と、前記間引かれたデータを用いることにより、サンプリング周期が異なる応答信号ごとに複数個の差分方程式を同定する差分方程式同定手段と、前記複数個の差分方程式について、各サンプリング周期と応答信号の測定期間で決まる周波数帯域の周波数応答特性を計算する周波数応答特性計算手段と、前記計算された複数個の周波数応答特性に基づいて、そのすべての周波数応答特性を表すことができる測定時のサンプリング周期での１つのパルス伝達関数を同定するパルス伝達関数同定手段とを備えることを特徴とする。

これにより、分割された各周波数帯域ごとに入力信号が入力される期間を異ならせることが可能となり、一定のサンプリング周期で入力信号を入力した場合においても、分割された各周波数帯域の入出力信号のデータ数を互いに一致させることが可能となる。このため、周波数帯域を分割した場合においても、高いサンプリング周波数で低い周波数帯域を同定する必要がなくなり、扱う総データ数と最小二乗解法による数値演算量の増大を伴うことなく、周波数応答特性を同定することが可能となるとともに、特性不明システムのモードの数に依存することなく、間引きを行うためのサンプリング周波数や周波数帯域の分割位置を設定することが可能となり、周波数帯域の分割するためのデータを別途測定したり、スペクトル解析による試行錯誤をなくすことを可能として、作業工程の増大を抑制することができる。

また、請求項２記載のモデル同定装置によれば、前記差分方程式同定手段は、未知変数の最大数が前記差分方程式を構成する連立方程式の個数を超えないように次数を変化させた時に、複数個の差分方程式の入力信号に対する出力と測定した応答信号との差の平方和が最小となるように差分方程式の次数を設定することを特徴とする。
これにより、差分方程式の次数を適切に決定することが可能となり、差分方程式から周波数応答特性を精度よく計算することができる。

また、請求項３記載のモデル同定装置によれば、前記パルス伝達関数同定手段は、全周波数応答特性に対する線形化最小二乗法に基づいてパルス伝達関数を同定することを特徴とする。
これにより、複数の周波数帯域の周波数応答特性から１つの伝達関数を推定する際に、分割して得られた周波数応答特性を結合させることなく、複数組の周波数応答特性をそのまま用いることで、その特性を数学的に近似したパルス伝達関数を同定することが可能となるとともに、周波数応答特性を分割してパルス伝達関数を同定する必要がなくなり、周波数応答特性に対する視覚的解析を基にした帯域分割のための領域設定などの事前の作業を不要とすることができる。

また、請求項４記載のモデル同定装置によれば、前記パルス伝達関数同定手段は、線形化されたパルス伝達関数に対して反復重み付き最小二乗法を適用し、導出されたパルス伝達関数の周波数応答特性と基準となる差分方程式の周波数応答特性との相対偏差が閾値を下回るように、前記パルス伝達関数の係数を計算することを特徴とする。
これにより、非線形最小二乗法を適用することなく、パルス伝達関数を同定することが可能となり、収束解が得られなくなる懸念を払拭することが可能となるとともに、パルス伝達関数の線形化に伴う偏差を低減することができる。

また、請求項５記載のモデル同定装置によれば、前記パルス伝達関数同定手段は、前記パルス伝達関数の分母多項式および分子多項式の各係数で構成されるシルベスタ行列の階数がそのシルベスタ行列の列数と等しくなるように前記パルス伝達関数の次数を決定することを特徴とする。
これにより、伝達関数の係数で構成されるシルベスタ行列の行列解析によりパルス伝達関数の次数を一意に決定することができる。このため、パルス伝達関数の次数をモードの数から決定する必要がなくなり、モードの探索やモード探索のための周波数応答解析などを不要となることから、パルス伝達関数の次数を決定するための作業工程の負担を軽減することができる。

また、請求項６記載のモデル同定プログラムによれば、一定のサンプリング周期で期間の異なる複数の入力信号を入力するステップと、前記入力信号に対する応答信号を前記サンプリング周期で測定するステップと、前記複数の入力信号に対する応答信号の組に対し、データ数がすべての組に対して等しくなるように各組ごとに異なるサンプリング周期で間引いてサンプリングするステップと、前記間引かれたデータを用いることにより、サンプリング周期が異なる応答信号ごと複数個の差分方程式を同定するステップと、前記複数個の差分方程式について、各サンプリング周期と応答信号の測定期間で決まる周波数帯域の周波数応答特性を計算するステップと、前記計算された複数個の周波数応答特性に基づいて、そのすべての周波数応答特性を表すことができる測定時のサンプリング周期での１つのパルス伝達関数を同定するステップとをコンピュータに実行させることを特徴とする。

これにより、モデル同定プログラムをコンピュータに実行させることで、周波数帯域を分割した場合においても、扱う総データ数と最小二乗解法による数値演算量の増大を伴うことなく周波数応答特性を同定することが可能となるとともに、特性不明システムのモードの数に依存することなく、間引きを行うためのサンプリング周波数や周波数帯域の分割位置を設定することが可能となり、周波数帯域の分割するためのデータを別途測定したり、スペクトル解析による試行錯誤をなくすことができ、作業工程の増大を抑制することができる。

また、請求項７記載のモデル同定プログラムによれば、未知変数の最大数が前記差分方程式を構成する連立方程式の個数を超えないように次数を変化させた時に、複数個の差分方程式の入力信号に対する出力と測定した応答信号との差の平方和が最小となるように差分方程式の次数を設定するステップをコンピュータに実行させることを特徴とする請求項６記載のモデル同定プログラム。
これにより、モデル同定プログラムをコンピュータに実行させることで、差分方程式の次数を適切に決定することが可能となり、差分方程式から周波数応答特性を精度よく計算することができる。

また、請求項８記載のモデル同定プログラムによれば、全周波数応答特性に対する線形化最小二乗法に基づいてパルス伝達関数を同定するステップをコンピュータに実行させることを特徴とする。
これにより、モデル同定プログラムをコンピュータに実行させることで、複数の周波数帯域の周波数応答特性から１つの伝達関数を推定する際に、分割して得られた周波数応答特性を結合させることなく、複数組の周波数応答特性をそのまま用いることで、その特性を数学的に近似したパルス伝達関数を同定することが可能となるとともに、周波数応答特性を分割してパルス伝達関数を同定する必要がなくなり、周波数応答特性に対する視覚的解析を基にした帯域分割のための領域設定などの事前の作業を不要とすることができる。

また、請求項９記載のモデル同定プログラムによれば、線形化されたパルス伝達関数に対して反復重み付き最小二乗法を適用し、導出されたパルス伝達関数の周波数応答特性と基準となる差分方程式の周波数応答特性との相対偏差が閾値を下回るように、前記パルス伝達関数の係数を計算するステップをコンピュータに実行させることを特徴とする。
これにより、モデル同定プログラムをコンピュータに実行させることで、非線形最小二乗法を適用することなく、パルス伝達関数を同定することが可能となり、収束解が得られなくなる懸念を払拭することが可能となるとともに、パルス伝達関数の線形化に伴う偏差を低減することができる。

また、請求項１０記載のモデル同定プログラムによれば、前記パルス伝達関数の分母多項式および分子多項式の各係数で構成されるシルベスタ行列の階数がそのシルベスタ行列の列数と等しくなるように前記パルス伝達関数の次数を決定するステップをコンピュータに実行させることを特徴とする。
これにより、モデル同定プログラムをコンピュータに実行させることで、伝達関数の係数で構成されるシルベスタ行列の行列解析によりパルス伝達関数の次数を一意に決定することができる。このため、パルス伝達関数の次数をモードの数から決定する必要がなくなり、モードの探索やモード探索のための周波数応答解析などを不要となることから、パルス伝達関数の次数を決定するための作業工程の負担を軽減することができる。

以上説明したように、本発明によれば、データ数や数値演算量の増大を伴うことなく、周波数帯域を分割しながら周波数応答特性を計算することが可能となるとともに、間引きを行うためのサンプリング周波数と周波数帯域の分割位置を客観的に一意に設定した上で、分割して得られた周波数応答特性を結合させたり、モードの数を考慮した帯域分割を行ったりすることなく、伝達関数を同定することが可能となる。

以下、本発明の実施形態に係るモデル同定装置およびその方法について図面を参照しながら説明する。
図１は、本発明の一実施形態に係るモデル同定方法を示すブロック図である。
図１において、特性不明システム１に入力信号ｕが入力されると、その入力信号ｕに対する応答信号ｙが出力される。また、特性不明システム１にはシステム同定装置２が接続され、システム同定装置２は、特性不明システム１に入力される入力信号ｕおよび特性不明システム１から出力される応答信号ｙのうち、ある一定期間のサンプリングデータｕ［ｎ］、ｙ［ｎ］を測定する。そして、システム同定装置２は、これらのサンプリングデータｕ［ｎ］、ｙ［ｎ］に基づいて、特性不明システム１の特性を示すモデルを同定する。なお、特性不明システム１の特性を示すモデルは、例えば、特性不明システム１の伝達関数で表すことができる。

図２は、図１のシステム同定装置２の概略構成を示すブロック図である。
図２において、システム同定装置２には、入力信号入力手段１１、応答信号測定手段１２、サンプリング手段１３、差分方程式同定手段１４、周波数応答特性計算手段１５およびパルス伝達関数同定手段１６が設けられている。ここで、入力信号入力手段１１は、一定のサンプリング周期で期間の異なる複数の入力信号ｕを入力する。なお、複数の入力信号ｕとしては、例えば、複数の擬似白色信号または多重正弦波を用いることができる。応答信号測定手段１２は、一定のサンプリング周期で入力された入力信号ｕに対する応答信号ｙを測定する。サンプリング手段１３は、複数の入力信号ｕに対する応答信号ｙの組に対し、データ数がすべての組に対して等しくなるように各組ごとに異なるサンプリング周期で間引いてサンプリングする。差分方程式同定手段１４は、サンプリング手段１３にて間引かれたデータを用いることにより、サンプリング周期が異なる応答信号ｙごとに複数個の差分方程式を同定する。周波数応答特性計算手段１５は、差分方程式同定手段１４にて同定された複数個の差分方程式について、各サンプリング周期と応答信号ｙの測定期間で決まる周波数帯域の周波数応答特性を計算する。パルス伝達関数同定手段１６は、周波数応答特性計算手段１５にて計算された複数個の周波数応答特性に基づいて、そのすべての周波数応答特性を表すことができる測定時のサンプリング周期での１つのパルス伝達関数を同定する。

これにより、システム同定装置２は、入力信号ｕが入力される期間を異ならせることで周波数帯域を分割しながら、特性不明システム１の周波数応答特性を計算することが可能となるとともに、一定のサンプリング周期で特性不明システム１に入力信号ｕを入力した場合においても、分割された各周波数帯域の入力信号ｕのデータ数を互いに一致させることが可能となる。このため、システム同定装置２は、特性不明システム１の周波数帯域を分割した場合においても、高いサンプリング周波数で低い周波数帯域を同定する必要がなくなり、扱う総データ数と最小二乗解法による数値演算量の増大を伴うことなく、周波数応答特性を同定することが可能となるとともに、特性不明システム１のモードの数に依存することなく、間引きを行うためのサンプリング周波数や周波数帯域の分割位置を設定することが可能となり、周波数帯域の分割するためのデータを別途測定したり、スペクトル解析による試行錯誤をなくすことを可能として、作業工程の増大を抑制することができる。

なお、入力信号入力手段１１、応答信号測定手段１２、サンプリング手段１３、差分方程式同定手段１４、周波数応答特性計算手段１５およびパルス伝達関数同定手段１６は、これらの手段で行われる処理を遂行させる命令が記述されたプログラムをコンピュータに実行させることにより実現することができる。
そして、このプログラムをＣＤ−ＲＯＭなどの記憶媒体に記憶しておけば、その記憶媒体をコンピュータに装着し、そのプログラムをコンピュータにインストールすることにより、入力信号入力手段１１、応答信号測定手段１２、サンプリング手段１３、差分方程式同定手段１４、周波数応答特性計算手段１５およびパルス伝達関数同定手段１６で行われる処理を実現することができる。また、このプログラムをインターネットを介してダウンロードすることにより、このプログラムを容易に普及させることができる。

また、入力信号入力手段１１、応答信号測定手段１２、サンプリング手段１３、差分方程式同定手段１４、周波数応答特性計算手段１５およびパルス伝達関数同定手段１６で行われる処理を遂行させる命令が記述されたプログラムをコンピュータに実行させる場合、スタンドアロン型コンピュータで実行させるようにしてもよく、ネットワークに接続された複数のコンピュータに分散処理させるようにしてもよい。

以下、図１のシステム同定装置２で実施される処理について数式を用いながら説明する。
図１において、特性不明システム１が線形動的システムであると仮定すると、あるサンプリング点ｎでの特性不明システム１の応答信号ｙは、そのサンプリング点ｎでの入力信号ｕと、過去にサンプリングされた入力信号ｕ［ｎ−ｉ］および応答信号ｙ［ｎ−ｉ］とを用いることで、以下の（１）式の差分方程式で表すことができる。

ただし、ａ_i、ｂ_iは差分方程式の未知係数、ｑは差分方程式の未知の次数である。
（１）式は、特性不明システム１の応答信号ｙを測定した期間Ｔの全てのサンプリング点ｎで成立し、差分方程式の次数ｑを設定すると、未知係数ａ_i、ｂ_iに対するｍ個の連立方程式が得られる。ただし、ｍ＝Ｔ／Δｔ、Δｔはサンプリング周期である。従って、（２）式に示すように、未知係数ａ_i、ｂ_iを変数ｘとした行列方程式Ａｘ＝ｂが（１）式から得られる。

例えば、１Ｈｚ〜１０００Ｈｚまでの周波数帯域を同定するとして、サンプリング周波数を４０００Ｈｚで測定すると、図２の入力信号入力手段１１は、表２に示すような期間の異なるＭ系列信号を入力信号ｕとして用いることができる。ここで、サンプリング周波数を４０００Ｈｚで一定の場合、入力信号ｕの各データ数は、表２に示すように、入力信号ｕが入力される期間に比例して増加する。例えば、サンプリング周波数が４０００Ｈｚで入力信号ｕを０．０６４秒、０．２５６秒、１．０２４秒の期間だけそれぞれサンプリングした場合、データ数はそれぞれ２５６個、１０２４個、４０９６個となる。

そこで、サンプリング手段１３は、入力される期間の最も短い入力信号ｕのデータ数と、その他の入力される期間の長い入力信号ｕのデータ数が互いに一致するように、入力される期間の長い入力信号ｕのデータ数を間引いてサンプリングする。この結果、例えば、表２に示すように、０．０６４秒、０．２５６秒、１．０２４秒の期間に対し、間引きが行われた後のサンプリング周波数として、４０００Ｈｚ、１０００Ｈｚ、２５０Ｈｚという互いに異なるサンプリング周波数がそれぞれ得られる。なお、間引きが行う前に、折り返し雑音を除去するためのフィルタリングを行うようにしてもよい。

この間引かれたサンプリング周波数の異なる３組のデータを用いて、（２）式に示すような２５６個の連立方程式を３組だけ作成する。そして、差分方程式同定手段１４は、これらの３組の連立方程式を最小二乗法によって解くことにより、サンプリング周波数の異なる応答信号ｙごとの差分方程式の未知係数ａ_i、ｂ_iを導出することができる。
ここで、周波数応答特性を精度よく計算するためには、差分方程式の次数ｑを適切に決定する必要がある。このため、差分方程式の次数ｑを２から１２７まで設定し、各次数ｑに対する未知係数ａ_i、ｂ_iを（２）式により導出する。そして、複数個の差分方程式の入力信号ｕに対する各出力を（３）式により計算する。

ただし、
ｏ［ｋ］：差分方程式のｋ番目の出力
ｕ［ｋ］：ｋ（ｋ＝０，・・・，２５５）番目の入力信号
である。
そして、複数個の差分方程式の入力信号ｕに対する各出力が求まると、測定した特性不明システム１の応答信号ｙとの差の平方和Ｅｑを（４）式にてそれぞれ計算し、平方和Ｅｑが最小となる次数ｑを最適値とすることができる。

そして、差分方程式の未知係数ａ_i、ｂ_iが導出されると、周波数応答特性計算手段１５は、３組の差分方程式を用いることにより、特性不明システム１の周波数応答特性を計算する。ここで、３組の差分方程式に対する周波数応答特性Ｇ_fs（ｚ）は、（５）式に示すように、これらに差分方程式をｚ演算子で構成されるパルス伝達関数に変換することにより計算することができる。ここで、ｚ^-iはｉサンプル遅延を表すため、（２）式により導出された差分方程式の係数ａ_i、ｂ_iをそのまま用いることができる。

ただし、
ｚ_k＝ｅｘｐ（ｊ２πｆ_k／ｆｓ）
ｆｓ：３組の差分方程式に対応する各サンプリング周波数
ｆ_k：各サンプリング周波数に対応する周波数帯域内のｋ番目の周波数点
である。
ここで、３組の差分方程式に対する周波数帯域と周波数点は、それぞれのサンプリング周波数ｆｓと測定期間Ｔにより、以下の（６）式、（７）式で表すことができる。
１／Ｔ≦ｆ_k≦ｆｓ／５ （６）
ｆ_k＝ｆｓ／Ｎ×ｋ （７）
ただし、
Ｎ：各サンプリング周波数での入力信号ｕのデータ数（表２の場合、２５６個）
ｋ：１以上の正の整数
である。

そして、パルス伝達関数同定手段１６は、周波数応答特性計算手段１５にて得られた３組の周波数応答データに対して、そのすべての周波数応答特性を表すことができる測定時のサンプリング周期での１つのパルス伝達関数を同定する。ここで、パルス伝達関数同定手段１６は、（８）式に示すように、全周波数応答特性に対する最小二乗法に基づいて、α_i、β_iを未知係数とした測定時のサンプリング周期でのパルス伝達関数を同定することができる。

これにより、複数の周波数帯域の周波数応答特性から１つのパルス伝達関数を推定する際に、分割して得られた周波数応答特性を結合させることなく、複数組の周波数応答特性をそのまま用いることで、その特性を数学的に近似したパルス伝達関数を同定することが可能となるとともに、周波数応答特性を分割してパルス伝達関数を同定する必要がなくなり、周波数応答特性に対する視覚的解析を基にした帯域分割のための領域設定などの事前の作業を不要とすることができる。
ただし、（８）式のままでは、非線形最小二乗法を適用することとなり、収束解が得られない場合が懸念される。このため、（８）式の分母を払うことにより（９）式を生成する。

そして、（１０）式に示すように、線形化された（９）式の全周波数応答特性に対する連立方程式を生成し、（１０）式の連立方程式を線形化最小二乗法により解いて、未知係数α_i、β_iを導出する。
Ａｘ＝ｄ （１０）
ただし、Ａ、ｄ、ｘはそれぞれ以下のように表すことができる。

ここで、
ｚ_k＝ｅｘｐ（ｊ２πｆ_k／ｆｓ）
ｆｓ：測定時のサンプリング周波数（表２の場合、４０００Ｈｚ）
ｆ_k：全周波数帯域のｋ番目の周波数点
である。
また、ｚ演算子は複素関数であるので、（１１）式に示すように、１０）式の連立方程式の実部と虚部とを別々に計算できるように構成し、連立代数方程式から未知係数α_i、β_iを導出する。
Ｂｘ＝ｅ （１１）
ただし、Ｂ、ｅはそれぞれ以下のように表すことができる。

しかし、（１１）式で得られる解には線形化による偏差が含まれるので、偏差を低減させるために、反復重み付き最小二乗法を適用する。すなわち、重み係数行列Ｗ⁽λ⁾を（１１）式に追加して、以下の（１１）式に示される連立代数方程式を構成する。
Ｃ⁽λ⁾ｘ＝ｇ⁽λ⁾ （１２）
ただし、λは繰り返し回数である。
また、Ｃ⁽λ⁾、ｇ⁽λ⁾は以下の式で表すことができる。

ここで、重み係数行列Ｗ⁽λ⁾は以下の式で表すことができる。

ここで、重み係数ω_i ⁽λ⁾は以下の式で表すことができる。

ただし、
α_i ⁽λ⁾：λ回目に得られたパルス伝達関数を構成する分母多項式のｉ番目の係数である。
ここで、係数α_i、β_iが未知であるため、初回は重み係数ω_i ⁽λ⁾はすべて１とする。そして、重み係数ω_i ⁽λ⁾を更新しながら、重み付き最小二乗解を繰り返し計算することで、最小二乗解の偏差を小さくすることができる。そして、導出されたパルス伝達関数の周波数応答特性Ｈ（ｚ）と基準となる差分方程式の周波数応答特性Ｇ_fs（ｚ）との相対偏差が、（１３）式に示すように、閾値γを下回った時に繰り返し作業を終了する。

ただし、パルス伝達関数の周波数応答特性Ｈ（ｚ）は、以下の式で表すことができる。

ここで、
ｎ：周波数点の総数
である
基準となる差分方程式の周波数応答特性Ｇ_fs（ｚ）を精度よく同定するためには、パルス伝達関数の次数を適切に設定する必要がある。そのため、以下の（１４）式に示すように、パルス伝達関数の分母多項式および分子多項式の各係数で構成されるシルベスタ行列Ｓを求め、シルベスタ行列Ｓの階数がそのシルベスタ行列Ｓの列数と等しくなるようにパルス伝達関数の次数を決定する。

これにより、伝達関数の係数で構成されるシルベスタ行列Ｓの行列解析によりパルス伝達関数の次数を一意に決定することができる。このため、パルス伝達関数の次数をモードの数から決定する必要がなくなり、モードの探索やモード探索のための周波数応答解析などを不要となることから、パルス伝達関数の次数を決定するための作業工程の負担を軽減することができる。

ここで、最適なパルス伝達関数の次数の決定手順としては、例えば、以下に示すような方法を用いることができる。
まず、初期次数を２次と設定し、線形化最小二乗法にて２次のパルス伝達関数を導出する。そして、合計６個の係数から４行４列のシルベスタ行列Ｓを構成し階数を計算する。そして、得られたシルベスタ行列Ｓの階数が列数と等しい場合、パルス伝達関数の次数を１つ増やして、３次のパルス伝達関数を導出する。そして、合計８個の係数から６行６列のシルベスタ行列Ｓを構成し階数を計算する。

もし、得られたシルベスタ行列Ｓの階数が列数より少ない場合、２次を最適な次数として決定する。一方、得られたシルベスタ行列Ｓの階数が列数と等しい場合、パルス伝達関数の次数を１つ増やして、４次のパルス伝達関数を導出する。もし、得られたシルベスタ行列Ｓの階数が列数より少ない場合、３次を最適な次数として決定する。一方、得られたシルベスタ行列Ｓの階数が列数と等しい場合、パルス伝達関数の次数を１つ増やして、同様な操作を繰り返すことにより、最適なパルス伝達関数の次数を決定することが可能となる。

本発明は、データ量や数値演算量の増大および作業工程の負担を抑制しつつ、特性不明システムのモデルを精度よく同定することができ、人工衛星の姿勢制御系やロボットの適応制御系に利用することができる。

本発明の一実施形態に係るモデル同定方法を示すブロック図である。 図１のシステム同定装置の概略構成を示すブロック図である。 従来のM系列入力信号に対する応答信号の周波数応答解析結果を示す図である。

符号の説明

１ 特性不明システム
２ システム同定装置
１１ 入力信号入力手段
１２ 応答信号測定手段
１３ サンプリング手段
１４ 差分方程式同定手段
１５ 周波数応答特性計算手段
１６ パルス伝達関数同定手段

Claims (10)

1. 一定のサンプリング周期で期間の異なる複数の入力信号を入力する入力信号入力手段と、
   前記入力信号に対する応答信号を前記サンプリング周期で測定する応答信号測定手段と、
   前記複数の入力信号に対する応答信号の組に対し、データ数がすべての組に対して等しくなるように各組ごとに異なるサンプリング周期で間引いてサンプリングするサンプリング手段と、
   前記間引かれたデータを用いることにより、サンプリング周期が異なる応答信号ごとに複数個の差分方程式を同定する差分方程式同定手段と、
   前記複数個の差分方程式について、各サンプリング周期と応答信号の測定期間で決まる周波数帯域の周波数応答特性を計算する周波数応答特性計算手段と、
   前記計算された複数個の周波数応答特性に基づいて、そのすべての周波数応答特性を表すことができる測定時のサンプリング周期での１つのパルス伝達関数を同定するパルス伝達関数同定手段とを備えることを特徴とするモデル同定装置。
2. 前記差分方程式同定手段は、未知変数の最大数が前記差分方程式を構成する連立方程式の個数を超えないように次数を変化させた時に、複数個の差分方程式の入力信号に対する出力と測定した応答信号との差の平方和が最小となるように差分方程式の次数を設定することを特徴とする請求項１記載のモデル同定装置。
3. 前記パルス伝達関数同定手段は、全周波数応答特性に対する線形化最小二乗法に基づいてパルス伝達関数を同定することを特徴とする請求項１または２記載のモデル同定装置。
4. 前記パルス伝達関数同定手段は、線形化されたパルス伝達関数に対して反復重み付き最小二乗法を適用し、導出されたパルス伝達関数の周波数応答特性と基準となる差分方程式の周波数応答特性との相対偏差が閾値を下回るように、前記パルス伝達関数の係数を計算することを特徴とする請求項３記載のモデル同定装置。
5. 前記パルス伝達関数同定手段は、前記パルス伝達関数の分母多項式および分子多項式の各係数で構成されるシルベスタ行列の階数がそのシルベスタ行列の列数と等しくなるように前記パルス伝達関数の次数を決定することを特徴とする請求項３または４記載のモデル同定装置。
6. 一定のサンプリング周期で期間の異なる複数の入力信号を入力するステップと、
   前記入力信号に対する応答信号を前記サンプリング周期で測定するステップと、
   前記複数の入力信号に対する応答信号の組に対し、データ数がすべての組に対して等しくなるように各組ごとに異なるサンプリング周期で間引いてサンプリングするステップと、
   前記間引かれたデータを用いることにより、サンプリング周期が異なる応答信号ごと複数個の差分方程式を同定するステップと、
   前記複数個の差分方程式について、各サンプリング周期と応答信号の測定期間で決まる周波数帯域の周波数応答特性を計算するステップと、
   前記計算された複数個の周波数応答特性に基づいて、そのすべての周波数応答特性を表すことができる測定時のサンプリング周期での１つのパルス伝達関数を同定するステップとをコンピュータに実行させることを特徴とするモデル同定プログラム。
7. 未知変数の最大数が前記差分方程式を構成する連立方程式の個数を超えないように次数を変化させた時に、複数個の差分方程式の入力信号に対する出力と測定した応答信号との差の平方和が最小となるように差分方程式の次数を設定するステップをコンピュータに実行させることを特徴とする請求項６記載のモデル同定プログラム。
8. 全周波数応答特性に対する線形化最小二乗法に基づいてパルス伝達関数を同定するステップをコンピュータに実行させることを特徴とする請求項６または７記載のモデル同定プログラム。
9. 線形化されたパルス伝達関数に対して反復重み付き最小二乗法を適用し、導出されたパルス伝達関数の周波数応答特性と基準となる差分方程式の周波数応答特性との相対偏差が閾値を下回るように、前記パルス伝達関数の係数を計算するステップをコンピュータに実行させることを特徴とする請求項８記載のモデル同定プログラム。
10. 前記パルス伝達関数の分母多項式および分子多項式の各係数で構成されるシルベスタ行列の階数がそのシルベスタ行列の列数と等しくなるように前記パルス伝達関数の次数を決定するステップをコンピュータに実行させることを特徴とする請求項８または９記載のモデル同定プログラム。

Images