(57)【特許請求の範囲】
【請求項１】
同定対象とする動的システムのインパルス応答を入力とするシステム同定装置であって、
前記インパルス応答から前記動的システムを記述する線形離散時間システムの直達項を同定して出力する直達項同定部と、
前記インパルス応答からブロックハンケル行列を生成して出力するブロックハンケル行列生成部と、
前記ブロックハンケル行列生成部から出力される前記ブロックハンケル行列の特異値分解により、前記ブロックハンケル行列の左特異ベクトルを列ベクトルとする第１の直交行列、前記ブロックハンケル行列の右特異ベクトルを列ベクトルとする第２の直交行列及び前記ブロックハンケル行列の特異値を出力する特異値分解部と、
前記第１の直交行列、前記第２の直交行列、前記特異値及び指定されたシステム次元の探索範囲に基づいて、当該探索範囲に属する各次元に対する前記線形離散時間システムのシステム行列のうち前記直達項以外のシステム行列を同定して、当該システム行列及び前記直達項に基づいて算出した前記線形離散時間システムのシステム特性と前記動的システムの実際のシステム特性との比較から、システム次元を決定して出力するシステム次元決定部と、
前記第１の直交行列、前記第２の直交行列、前記特異値及び前記システム次元決定部から出力される前記システム次元に基づいて、前記線形離散時間システムのシステム行列のうち前記直達項以外のシステム行列を同定するシステム行列同定部と、を備え、
前記直達項同定部で同定した前記直達項及び前記システム行列同定部で同定した前記システム行列を前記線形離散時間システムとして出力することを特徴とするシステム同定装置。
【請求項２】
前記システム次元決定部が、
前記探索範囲に属する各次元に対して、前記線形離散時間システムのシステム特性を推定インパルス応答として算出して出力するシステム特性推定部と、
前記システム特性推定部から出力される前記線形離散時間システムの前記推定インパルス応答と前記動的システムの実際のインパルス応答との時間領域における誤差２乗和ノルムがしきい値以下となる次元のうち最小の次元をシステム次元として決定して出力するシステム次元推定部と、を備えることを特徴とする請求項１に記載のシステム同定装置。
【請求項３】
前記システム次元決定部が、
前記探索範囲に属する各次元に対して、前記線形離散時間システムのシステム特性を推定周波数応答として算出して出力するシステム特性推定部と、
前記システム特性推定部から出力される前記線形離散時間システムの推定周波数応答と前記インパルス応答をフーリエ変換した実際の周波数応答との周波数領域における誤差２乗和ノルムがしきい値以下となる次元のうち最小の次元をシステム次元として決定して出力するシステム次元推定部と、を備えることを特徴とする請求項１に記載のシステム同定装置。
【請求項４】
前記システム次元推定部が、
前記フーリエ変換した前記実際の周波数応答に基づいて重み関数を決定し、
前記システム特性推定部から出力される前記線形離散時間システムの推定周波数応答と前記動的システムの前記実際の周波数応答との周波数領域における誤差２乗値に対して、当該重み関数を乗じた値の加算値を算出して、当該加算値のノルムがしきい値以下となる次元のうち最小の次元をシステム次元として決定して出力することを特徴とする請求項３に記載のシステム同定装置。
【請求項５】
前記しきい値は、誤差２乗和許容値と誤差２乗和ノルムの最小値の積により規定されることを特徴とする請求項２から４のいずれか一項に記載のシステム同定装置。
【請求項６】
前記システム次元決定部が、
前記探索範囲の中で第１の次元よりも１段階低い第２の次元に対応したシステム行列の同定結果と、前記第１の直交行列、前記第２の直交行列及び前記特異値からそれぞれ前記第２の次元よりも大きく前記第１の次元以下となる左特異ベクトル、右特異ベクトル及び特異値と、を用いて、再帰的手法により前記第１の次元に対応するシステム行列を同定する再帰的システム行列推定部を備えることを特徴とする請求項１に記載のシステム同定装置。
【請求項７】
前記インパルス応答に対応する情報が前記動的システムを打撃加振したときのシステム入力及びシステム出力であり、
前記システム入力が有意な値を持つ複数の時刻を特定するシステム入力印加時刻特定部と、
前記システム入力印加時刻特定部から出力される複数の時刻に対応するシステム入力の加算値をインパルス入力振幅とし、前記システム入力印加時刻特定部から出力される複数の時刻の最大値をインパルス入力印加時刻として、当該インパルス入力印加時刻以降のシステム出力を前記インパルス入力振幅で除した信号を前記インパルス応答として出力するインパルス応答変換部と、を備え、
前記同定対象とする前記動的システムを打撃加振した場合のシステム入力及びシステム出力並びに前記探索範囲を入力として、前記線形離散時間システムを同定することを特徴とする請求項１に記載のシステム同定装置。
【請求項８】
前記システム入力印加時刻特定部が、
比率しきい値とシステム入力の最大値とを乗じた値をシステム入力しきい値とし、システム入力がシステム入力しきい値以上となる時刻を、システム入力が有意な値を持つ複数の時刻として特定することを特徴とする請求項７に記載のシステム同定装置。
【請求項９】
前記インパルス応答に対応する情報が前記動的システムの周波数応答であり、
前記周波数応答の逆フーリエ変換により対応するインパルス応答を出力する逆フーリエ変換部を備え、
前記同定対象とする前記動的システムの周波数応答及び前記探索範囲を入力として、前記線形離散時間システムを同定することを特徴とする請求項１に記載のシステム同定装置。
【請求項１０】
前記システム次元決定部が、
前記探索範囲に属する各次元に対して前記線形離散時間システムの安定性を評価するシステム安定性評価部を備え、
安定システムとなる次元に対応する前記線形離散時間システムのシステム特性からシステム次元を決定することを特徴とする請求項１に記載のシステム同定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、システム同定装置に関するものである。
【背景技術】
【０００２】
従来のインパルス応答によるシステム同定装置として、例えば非特許文献１にはＨｏ−Ｋａｌｍａｎの方法に基づくものが開示されている。この方法では、線形離散時間システム（Ａｄ，Ｂｄ，Ｃｄ，Ｄｄ）で記述される動的システムのインパルス応答（Ｇ０，Ｇ１，Ｇ２，…）に基づいて、Ｇ０からシステムの直達項Ｄｄを決定すると共に、Ｇ１，Ｇ２，…からブロックハンケル行列Ｈｋｌを生成する。次に、ブロックハンケル行列Ｈｋｌを特異値分解し、有意な値を持つ特異値の個数をシステム次元として決定して、特異値分解の結果と決定したシステム次元から、拡大可観測行列Ｏｋ及び拡大可到達行列Ｃｌを算出する。最後に、拡大可観測行列Ｏｋ及び拡大可到達行列Ｃｌに基づいてシステム行列Ａｄ，Ｂｄ，Ｃｄを算出することで、動的システムを記述する線形離散時間システム（Ａｄ，Ｂｄ，Ｃｄ，Ｄｄ）を同定している。【０００３】
また、従来のインパルス応答によるシステム同定装置の他の例として、例えば特許文献１にはプラント・モデリング装置が開示されている。このプラント・モデリング装置では、動的システムのインパルス応答（Ｇ０，Ｇ１，Ｇ２，…）に対して上記Ｈｏ−Ｋａｌｍａｎの方法を適用するが、システム次元の決定方法として次の２種類を選択的に適用している。第１のシステム次元の決定方法は、特異値と当該特異値に対応する次数の関係をグラフィック端末に対数スケールで表示し、有意な値を持つ特異値の個数、つまりシステム次元を作業者が決定する方法である。第２のシステム次元の決定方法は、特異値の変化率と観測ノイズに基づく評価関数を適用し、当該評価関数を最小とする次元をシステム次元として自動的に決定する方法である。【０００４】
このようなインパルス応答によるシステム同定装置では、有意な値を持つ特異値の個数からシステム次元を決定している。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開昭６１−２６７１０２号公報
【非特許文献】
【０００６】
【非特許文献１】片山徹著、「システム同定−部分空間法からのアプローチ−」、朝倉書店、２００４年２月、ｐ．１０２−１０７
【発明の概要】
【発明が解決しようとする課題】
【０００７】
しかしながら、上記従来の技術によれば、現実のインパルス応答から算出したブロックハンケル行列の特異値は、なだらかな単調減少となる場合が多数存在し、この場合は有意な値を持つ特異値と、無視できる微小な値となる特異値との境界が不明確となる。そのため、第１のシステム次元の決定方法による従来のシステム同定装置では、システム次元の決定が作業者の判断に依存し、常に最適なシステム次元を決定しているとは限らず、またはシステム次元の決定に試行錯誤が必要になる、という問題があった。【０００８】
また、このような問題に対応する方法として、特許文献１では第２のシステム次元の決定方法を適用しているが、本手法でも評価関数の与え方によって決定されるシステム次元が変化するため、システム次元の決定が評価関数の与え方に依存し、本手法でも常に最適なシステム次元を決定しているとは限らず、またはシステム次元の決定に関して試行錯誤が必要であり、問題を解決するには至っていない。【０００９】
本発明は、上記に鑑みてなされたものであって、現実のインパルス応答から算出したブロックハンケル行列の特異値がなだらかな単調減少となり、有意な値を持つ特異値と無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、常に最適なシステム次元を決定し、動的システムを記述する線形離散時間システムの同定が可能なシステム同定装置を得ることを目的とする。
【課題を解決するための手段】
【００１０】
上述した課題を解決し、目的を達成するために、本発明は、同定対象とする動的システムのインパルス応答を入力とするシステム同定装置であって、前記インパルス応答から前記動的システムを記述する線形離散時間システムの直達項を同定して出力する直達項同定部と、前記インパルス応答からブロックハンケル行列を生成して出力するブロックハンケル行列生成部と、前記ブロックハンケル行列生成部から出力される前記ブロックハンケル行列の特異値分解により、前記ブロックハンケル行列の左特異ベクトルを列ベクトルとする第１の直交行列、前記ブロックハンケル行列の右特異ベクトルを列ベクトルとする第２の直交行列及び前記ブロックハンケル行列の特異値を出力する特異値分解部と、前記第１の直交行列、前記第２の直交行列、前記特異値及び指定されたシステム次元の探索範囲に基づいて、当該探索範囲に属する各次元に対する前記線形離散時間システムのシステム行列のうち前記直達項以外のシステム行列を同定して、当該システム行列及び前記直達項に基づいて算出した前記線形離散時間システムのシステム特性と前記動的システムの実際のシステム特性との比較から、システム次元を決定して出力するシステム次元決定部と、前記第１の直交行列、前記第２の直交行列、前記特異値及び前記システム次元決定部から出力される前記システム次元に基づいて、前記線形離散時間システムのシステム行列のうち前記直達項以外のシステム行列を同定するシステム行列同定部と、を備え、前記直達項同定部で同定した前記直達項及び前記システム行列同定部で同定した前記システム行列を前記線形離散時間システムとして出力することを特徴とする。
【発明の効果】
【００１１】
この発明によれば、同定対象とする動的システムにおいて、現実のインパルス応答から算出したブロックハンケル行列の特異値がなだらかな単調減少となり、有意な値を持つ特異値と無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、常に最適なシステム次元の決定及び動的システムを記述する線形離散時間システムの同定が可能となる、という効果を奏する。
【図面の簡単な説明】
【００１２】
【図１】図１は、実施の形態１にかかるシステム同定装置の全体構成を示すブロック線図である。
【図２】図２は、実施の形態１にかかるシステム同定装置において、システム入出力の時間波形を示す概略図である。
【図３】図３は、実施の形態１にかかるシステム同定装置において、ブロックハンケル行列の特異値と次元の関係を示す概略図である。
【図４】図４は、実施の形態１にかかるシステム同定装置において、システム次元決定部の内部構成と動作を示すブロック線図である。
【図５】図５は、実施の形態１にかかるシステム同定装置において、同定した線形離散時間システムの推定インパルス応答とシステムの実際のインパルス応答との誤差２乗和ノルムと、次元との関係を示す概略図である。
【図６】図６は、実施の形態２にかかるシステム同定装置において、動的システムを打撃加振した場合のシステム入出力の時間波形を示す概略図である。
【図７】図７は、実施の形態２にかかるシステム同定装置において、動的システムを打撃加振した場合のシステム入出力を、当該システムのインパルス応答に変換する手順を示すブロック線図である。
【図８】図８は、実施の形態２にかかるシステム同定装置において、システム次元決定部の内部構成と動作を示すブロック線図である。
【図９】図９は、実施の形態３にかかるシステム同定装置において、動的システムの周波数応答を、当該システムのインパルス応答に変換する手順を示すブロック線図である。
【図１０】図１０は、実施の形態４にかかるシステム同定装置の全体構成を示すブロック線図である。
【発明を実施するための形態】
【００１３】
以下に、本発明にかかるシステム同定装置の実施の形態を図面に基づいて詳細に説明する。なお、この実施の形態によりこの発明が限定されるものではない。
【００１４】
実施の形態１．
図１は、本発明にかかるシステム同定装置の実施の形態１の全体構成を示すブロック線図である。図２は、本実施の形態にかかるシステム同定装置において、システム入出力の時間波形を示す概略図である。図１に示すシステム同定装置１０では、同定対象とする動的システムのインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）を入力とする。このインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）は、図２に示すように、時刻０では振幅が１となり、当該時刻以外では振幅が０となる理想的なインパルス入力で動的システムへのシステム入力１１を構成した場合における時刻０以降のシステム出力１２で与えられる。
【００１５】
システム同定装置１０では、こうして得られた動的システムのインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）から、直達項同定部１によって動的システムを記述する線形離散時間システムの直達項Ｄｄを同定し、更にブロックハンケル行列生成部２によってブロックハンケル行列Ｈｋｌを生成する。
【００１６】
特異値分解部３では、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｈｋｌを特異値分解し、ブロックハンケル行列Ｈｋｌの左特異ベクトルを列ベクトルとする第１の直交行列Ｕ、ブロックハンケル行列Ｈｋｌの右特異ベクトルを列ベクトルとする第２の直交行列Ｖ及びブロックハンケル行列Ｈｋｌの特異値σｉ（ｉ＝１，２，３，…）を出力する。
【００１７】
システム次元決定部４では、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σｉ（ｉ＝１，２，３，…）と、作業者が指定したシステム次元の探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に基づいて、当該探索範囲に属する次元ｎｉ（ｉ＝１，２，…，ａ）のそれぞれに対して動的システムを記述する線形離散時間システムの直達項Ｄｄを除くシステム行列を同定する。更に、当該システム行列及び直達項同定部１から出力される直達項Ｄｄに基づいて、当該探索範囲に属する次元ｎｉ（ｉ＝１，２，…，ａ）のそれぞれに対して線形離散時間システムの推定インパルス応答を算出し、動的システムの実際のインパルス応答（図１における動的システムのシステム特性）との比較から、システム次元ｎを決定する。
【００１８】
システム行列同定部５では、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σｉ（ｉ＝１，２，３，…）と、システム次元決定部４から出力されるシステム次元ｎに基づいて、線形離散時間システムの直達項Ｄｄを除くシステム行列Ａｄ，Ｂｄ，Ｃｄを同定する。
【００１９】
システム同定装置１０では、直達項同定部１で同定した直達項Ｄｄ及びシステム行列同定部５で同定したシステム行列Ａｄ，Ｂｄ，Ｃｄを、動的システムを記述する線形離散時間システムとして最終的に出力する。
【００２０】
図３は、本実施の形態にかかるシステム同定装置１０において、ブロックハンケル行列Ｈｋｌの特異値σｉと次元（ｉ＝１，２，３，…）の関係を示す概略図である。図４は、本実施の形態にかかるシステム同定装置１０において、システム次元決定部４の内部構成と動作を示すブロック線図である。図５は、本実施の形態にかかるシステム同定装置１０において、同定した線形離散時間システムの推定インパルス応答と動的システムの実際のインパルス応答との時間領域における誤差２乗和ノルム||ｅｎｉ||と、次元ｎｉ（ｉ＝１，２，…，ａ）との関係を示す概略図である。なお、図４において、「ｉ＝１」の処理、「ｉ＋＋」の処理及び「ｉ≦ａ」の判断は、システム特性推定部３２またはシステム次元推定部３３が担ってもよいし、図示しない他の構成が担ってもよい。
【００２１】
図３に示すように、動的システムのインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）から生成したブロックハンケル行列Ｈｋｌの特異値σｉは、理想的なインパルス応答では次元（ｉ＝１，２，３，…）に対して、例えば特異値分布２１に示す関係となる。
【００２２】
この場合、有意な値を持つ特異値の個数を明確に規定することができ、当該個数が動的システムのシステム次元ｎに対応する（図３の場合はシステム次元ｎ＝４）。一方、観測ノイズ等の影響を受ける現実のインパルス応答に基づいて算出した特異値σｉは、次元（ｉ＝１，２，３，…）に対して例えば特異値分布２２に示す関係となるため、有意な値を持つ特異値と無視できる微小な値となる特異値との境界が不明確となり、常に最適なシステム次元ｎを決定しているとは限らず、またはシステム次元ｎの決定に関して試行錯誤が必要になるという問題がある。
【００２３】
そこで、システム次元決定部４は、図４に示すように再帰的システム行列推定部３１により、作業者が指定したシステム次元の探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に属する第１の次元ｎｉに関して、第１の次元ｎｉよりも１段階低い第２の次元ｎｉ−１に対応したシステム行列Ａｄ，ｎｉ−１，Ｂｄ，ｎｉ−１，Ｃｄ，ｎｉ−１の同定結果と、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ、及び特異値σｉ（ｉ＝１，２，３，…）のうち、それぞれ第２の次元ｎｉ−１よりも大きく、第１の次元ｎｉ以下となる左特異ベクトルｕｊ、右特異ベクトルｖｊ及び特異値σｊ（ｊ＝ｎｉ−１＋１，ｎｉ−１＋２，…，ｎｉ）を用いて、再帰的手法によって第１の次元ｎｉに対応するシステム行列Ａｄ，ｎｉ，Ｂｄ，ｎｉ，Ｃｄ，ｎｉを同定する。
【００２４】
システム次元決定部４が、探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）の中で第１の次元ｎｉよりも１段階低い第２の次元ｎｉ−１に対応したシステム行列Ａｄ，ｎｉ−１，Ｂｄ，ｎｉ−１，Ｃｄ，ｎｉ−１の同定結果と、第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σｉ（ｉ＝１，２，３，…）からそれぞれ第２の次元ｎｉ−１よりも大きく第１の次元ｎｉ以下となる左特異ベクトルｕｊ、右特異ベクトルｖｊ及び特異値σｊ（ｊ＝ｎｉ−１＋１，ｎｉ−１＋２，…，ｎｉ）と、を用いて、再帰的手法により第１の次元ｎｉに対応するシステム行列Ａｄ，ｎｉ，Ｂｄ，ｎｉ，Ｃｄ，ｎｉを同定する再帰的システム行列推定部３１を備えるシステム同定装置１０は、現実の動的システムに対して一致度の高いシステム次元を決定するための演算量を低減することが可能である。
【００２５】
次に、システム特性推定部３２は、システム次元の探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に属するそれぞれの次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａｄ，ｎｉ，Ｂｄ，ｎｉ，Ｃｄ，ｎｉ及び直達項同定部１から出力される直達項Ｄｄに基づいて、同定した線形離散時間システムの推定インパルス応答を算出する。そして、ｉに１を加算し、ｉ≦ａの場合には再帰的システム行列推定部３１により、システム行列Ａｄ，ｎｉ，Ｂｄ，ｎｉ，Ｃｄ，ｎｉを同定し、ｉ＞ａの場合にはシステム次元推定部３３の処理に進む。
【００２６】
システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムの推定インパルス応答と動的システムの実際のインパルス応答（図４における動的システムのシステム特性）との時間領域における誤差２乗和ｅｎｉ（ｉ＝１，２，…，ａ）を算出し、図５に示すように誤差２乗和ノルム分布４１が誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ６）。
【００２７】
すなわち、システム次元決定部４が、探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）に属する各次元に対して、線形離散時間システムのシステム特性を推定インパルス応答として算出して出力するシステム特性推定部３２と、システム特性推定部３２から出力される線形離散時間システムの推定インパルス応答と動的システムの実際のインパルス応答との時間領域における誤差２乗和ノルムがしきい値以下となる次元のうち最小の次元をシステム次元として決定して出力するシステム次元推定部３３と、を備えるシステム同定装置１０は、現実のインパルス応答から算出したブロックハンケル行列の特異値がなだらかな単調減少となり、有意な値を持つ特異値と無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、現実の動的システムに対して時間領域で一致度の高いシステム次元の決定と、動的システムを記述する線形離散時間システムの同定が可能となる。
【００２８】
次に、動作について説明する。同定対象とする動的システムが下記の式（１）に示す１入力Ｐ出力のｎ次元線形離散時間システムで記述できるとする。上記の動的システムへのシステム入力ｕ（ｊＴｓ）を、図２に示すシステム入力１１、つまり理想的なインパルス入力となるように下記の式（２）で構成すると、下記の式（１）に対応するシステム出力ｙ（ｊＴｓ）、つまり図２に示すシステム出力１２は下記の式（３）となり、動的システムのインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）は、時刻０以降のシステム出力１２で与えられる。
【００２９】
【数１】
0006009106.tif 000002
【００３０】
ただし、状態ベクトルｘ∈Ｒｎであり、システム入力ｕ∈Ｒであり、システム出力ｙ∈ＲＰであり、システム行列Ａｄ∈Ｒｎ×ｎ，Ｂｄ∈Ｒｎ，Ｃｄ∈ＲＰ×ｎ，Ｄｄ∈ＲＰであり、インパルス応答列ｇ（ｊＴｓ）（ｊ＝０，１，２，…）である。
【００３１】
【数２】
0006009106.tif 000003
【００３２】
【数３】
0006009106.tif 000004
【００３３】
図１に示すように、本実施の形態にかかるシステム同定装置１０は、上記で得られた動的システムのインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）を入力とし、直達項同定部１において上記の式（１）にて動的システムを記述する線形離散時間システムの直達項Ｄｄを下記の式（４）で同定し、更にブロックハンケル行列生成部２によって下記の式（５）で与えられるブロックハンケル行列Ｈｋｌを生成する。
【００３４】
【数４】
0006009106.tif 000005
【００３５】
【数５】
0006009106.tif 000006
【００３６】
次に、特異値分解部３によって、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｈｋｌを特異値分解し、下記の式（６）で与えられるブロックハンケル行列Ｈｋｌの左特異ベクトルｕｊを列ベクトルとする第１の直交行列Ｕ、ブロックハンケル行列Ｈｋｌの右特異ベクトルｖｊを列ベクトルとする第２の直交行列Ｖ及びブロックハンケル行列Ｈｋｌの特異値σｉ（ｉ＝１，２，３，…）を出力する。
【００３７】
【数６】
0006009106.tif 000007
【００３８】
ただし、第１の直交行列Ｕ＝［ｕ１  ｕ２  …  ｕｋＰ］∈ＲｋＰ×ｋＰであり、第２の直交行列Ｖ＝［ｖ１  ｖ２  …  ｕｌ］∈Ｒｌ×ｌであり、ブロックハンケル行列Ｈｋｌの特異値σ１≧σ２≧…≧σｎ≧σｎ＋１≧σｎ＋２≧…であり、Σは下記の式（７）で表される。
【００３９】
【数７】
0006009106.tif 000008
【００４０】
対象とする動的システムのシステム次元ｎは、ブロックハンケル行列Ｈｋｌの特異値において有意な値を持つものはｎ個であり、それらに比してｎ＋１個目以降は十分小さく、無視できる微小な値であるとする下記の式（８）の関係に基づいて決定することができる。
【００４１】
【数８】
0006009106.tif 000009
【００４２】
図３に示すように、動的システムのインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）から生成したブロックハンケル行列Ｈｋｌの特異値σｉは、理想的なインパルス応答では次元（ｉ＝１，２，３，…）に対して例えば特異値分布２１に示す関係となる。この場合、有意な値を持つ特異値の個数を明確に規定することができ、当該個数から動的システムのシステム次元ｎを決定することができる（図３の場合はシステム次元ｎ＝４）。
【００４３】
一方、観測ノイズ等の影響を受ける現実のインパルス応答に基づいて算出した特異値σｉは、次元（ｉ＝１，２，３，…）に対して例えば特異値分布２２に示す関係となるため、有意な値を持つ特異値と、無視できる微小な値となる特異値との境界が不明確となる。したがって、従来の手法では常に最適なシステム次元ｎを決定しているとは限らず、またはシステム次元ｎの決定に関して試行錯誤が必要になるという課題があった。そこで、本実施の形態のシステム同定装置１０では、「時間領域において実際のインパルス応答に最も適合する」という前提の下、システム次元決定部４において最適なシステム次元ｎを決定する。
【００４４】
システム次元決定部４は、図１に示すように特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σｉ（ｉ＝１，２，３，…）と、作業者が指定したシステム次元の探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に基づいて、当該探索範囲に属する次元ｎｉ（ｉ＝１，２，…，ａ）のそれぞれに対して動的システムを記述する線形離散時間システムの直達項Ｄｄを除くシステム行列を同定する。
【００４５】
更に、当該システム行列及び直達項同定部１から出力される直達項Ｄｄに基づいて、当該探索範囲に属する次元ｎｉ（ｉ＝１，２，…，ａ）のそれぞれに対して同定した線形離散時間システムの推定インパルス応答を算出し、動的システムの実際のインパルス応答（図１における動的システムのシステム特性）との比較から、システム次元ｎを決定する。具体的には図４に示すように、再帰的システム行列推定部３１によって、作業者が指定したシステム次元の探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に属する第１の次元ｎｉに関して、第１の次元ｎｉよりも１段階低い第２の次元ｎｉ−１に対応したシステム行列Ａｄ，ｎｉ−１，Ｂｄ，ｎｉ−１，Ｃｄ，ｎｉ−１の同定結果と、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σｉ（ｉ＝１，２，３，…）のうち、それぞれ第２の次元ｎｉ−１よりも大きく、第１の次元ｎｉ以下となる左特異ベクトルｕｊ、右特異ベクトルｖｊ及び特異値σｊ（ｊ＝ｎｉ−１＋１，ｎｉ−１＋２，…，ｎｉ）とを用いて、次式に示す再帰的手法によって第１の次元ｎｉに対応するシステム行列Ａｄ，ｎｉ，Ｂｄ，ｎｉ，Ｃｄ，ｎｉを同定する。まず、第１の次元ｎｉに対応する拡大可観測行列は下記の式（９）にて表される。
【００４６】
【数９】
0006009106.tif 000010
【００４７】
第１の次元ｎｉに対応する拡大可到達行列は下記の式（１０）にて表される。
【００４８】
【数１０】
0006009106.tif 000011
【００４９】
そして、第１の次元ｎｉに対応するシステム行列Ａｄ，ｎｉ，Ｂｄ，ｎｉ，Ｃｄ，ｎｉは下記の式（１１）にて表される。
【００５０】
【数１１】
0006009106.tif 000012
【００５１】
次にシステム特性推定部３２において、システム次元の探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に属する次元ｎｉのそれぞれに対して、再帰的システム行列推定部３１から出力されるシステム行列Ａｄ，ｎｉ，Ｂｄ，ｎｉ，Ｃｄ，ｎｉ及び直達項同定部１から出力される直達項Ｄｄに基づいて、上記の式（１）に対するシステム入力ｕ（ｊＴｓ）として上記の式（２）を適用することで、同定した線形離散時間システムの推定インパルス応答ｇ＾ｎｉ（ｊＴｓ）（ｊ＝０，１，２，…）を算出する。
【００５２】
なお、本明細書中で、「ｇ＾」は「ｇ」の上に「＾（サーカムフレックス）」が配置された文字の代替表記である。
【００５３】
システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムの推定インパルス応答ｇ＾ｎｉ（ｊＴｓ）（ｊ＝０，１，２，…）と動的システムの実際のインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）（図４における動的システムのシステム特性）との時間領域における誤差２乗和ｅｎｉ（ｉ＝１，２，…，ａ）を下記の式（１２）で算出する。
【００５４】
【数１２】
0006009106.tif 000013
【００５５】
この誤差２乗和ノルム||ｅｎｉ||を最小とする次元ｎｉが「時間領域において実際のインパルス応答に最も適合する」システム次元ｎとなる。一方、実際のノルム||ｅｎｉ||は、観測ノイズが白色性であれば、そのノイズレベルによらず次元ｎｉの増加に伴って単調減少し、図５に示すように所定の次元（図５の場合はｎ６）以上でほぼ一定値となる。そこで、ここでは、システム次元ｎの推定値が必要以上に高次元となることを避けるため、下記の式（１３）で与えられる誤差２乗和ノルムしきい値４２を規定し、誤差２乗和ノルム分布４１が誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ６）。
【００５６】
【数１３】
0006009106.tif 000014
【００５７】
最後に、システム行列同定部５において、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σｉ（ｉ＝１，２，３，…）と、システム次元決定部４から出力されるシステム次元ｎと、に基づき、線形離散時間システムの直達項Ｄｄを除くシステム行列Ａｄ，Ｂｄ，Ｃｄを下記の式（１４）で同定する。
【００５８】
【数１４】
0006009106.tif 000015
【００５９】
システム同定装置１０は、システム行列同定部５で同定したシステム行列Ａｄ，Ｂｄ，Ｃｄ及び直達項同定部１で同定した直達項Ｄｄを、動的システムを記述する線形離散時間システムとして最終的に出力する。
【００６０】
このように、本実施の形態のシステム同定装置１０により、現実のインパルス応答から算出したブロックハンケル行列Ｈｋｌの特異値σｉ（ｉ＝１，２，３，…）がなだらかな単調減少となり、有意な値を持つ特異値と同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して時間領域で一致度の高いシステム次元ｎの決定及び動的システムを記述する線形離散時間システムの同定が可能となる。
【００６１】
また、再帰的システム行列推定部３１により、現実の動的システムに対して一致度の高いシステム次元ｎを決定するための演算量を低減することが可能となる。
【００６２】
なお、本実施の形態のシステム同定装置１０では、線形離散時間システムのシステム特性を推定インパルス応答として算出し、当該インパルス応答と動的システムの実際のインパルス応答との時間領域における誤差２乗和ノルム分布４１が誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定している。しかしながら、本発明はこれに限定されず、線形離散時間システムのシステム特性を推定周波数応答として算出し、当該周波数応答と動的システムのインパルス応答をフーリエ変換した実際の周波数応答との周波数領域における誤差２乗和に基づいて、システム次元ｎを決定してもよい。この場合、更に動的システムの実際の周波数応答に基づいて重み関数を決定し、線形離散時間システムの推定周波数応答と動的システムの実際の周波数応答との周波数領域における誤差２乗値に当該重み関数を乗じた値の加算値に基づいて、システム次元ｎを決定してもよい。
【００６３】
実施の形態２．
本実施の形態では、動的システムの打撃加振によって得られる現実のシステム入出力も適用可能にすると共に、現実の動的システムが安定であることが明確である場合に安定システムに限定して同定することができるシステム同定装置について説明する。
【００６４】
動的システムのインパルス応答を得るための方法としては、対象とする動的システムをインパルスハンマによって打撃加振し、そのときのシステム出力を計測する方法が最も一般的である。しかしながら、打撃加振によって与えられる現実のシステム入力は正弦半波状となるため、真のインパルス入力とは異なる。したがって、当該入力によって得られるシステム出力も、動的システムの真のインパルス応答とは異なるものとなるため、当該システム出力をそのままＨｏ−Ｋａｌｍａｎの方法に適用することはできない、という問題があった。
【００６５】
加えて、従来のインパルス応答によるシステム同定装置では、同定結果として得られる線形離散時間システム（Ａｄ，Ｂｄ，Ｃｄ，Ｄｄ）の安定性は全く考慮されていないため、現実の動的システムが安定であるにも関わらず、不安定システムとして同定される場合がある、という問題もあった。
【００６６】
本実施の形態では、上記問題に鑑み、動的システムの打撃加振によって得られる現実のシステム入出力も適用可能にすると共に、現実の動的システムが安定であることが明確である場合に、安定システムに限定して同定することができるシステム同定装置を得ることを目的とする。
【００６７】
本実施の形態において、システム同定装置の全体構成を示すブロック線図は実施の形態１の図１と同一であり、ブロックハンケル行列Ｈｋｌの特異値σｉと次元（ｉ＝１，２，３，…）の関係を示す概略図は実施の形態１の図３と同一であり、同定した線形離散時間システムの推定システム出力と動的システムの実際のシステム出力との時間領域における誤差２乗和ノルム||ｅｎｉ||と、次元ｎｉ（ｉ＝１，２，…，ａ）との関係を示す概略図は実施の形態１の図５と同一である。
【００６８】
図６は、本実施の形態のシステム同定装置において、動的システムを打撃加振した場合のシステム入出力の時間波形を示す概略図である。図７は、本実施の形態のシステム同定装置において、動的システムを打撃加振した場合のシステム入出力を、当該システムのインパルス応答に変換する手順を示すブロック線図である。このとき、図７のインパルス応答変換部５２の出力である動的システムのインパルス応答が、図１の直達項同定部１及びブロックハンケル行列生成部２の入力となる。同定対象とする動的システムを打撃加振した場合のシステム入力１１及びシステム出力１２を入力とし、図７に示すシステム入力印加時刻特定部５１及びインパルス応答変換部５２並びに図１のシステム同定装置１０全体を併せて、打撃加振によるシステム同定装置５０とする。
【００６９】
図６に示すように、本実施の形態の打撃加振によるシステム同定装置５０では、同定対象とする動的システムを打撃加振した場合のシステム入力１１及びシステム出力１２に基づいて、当該システムを記述する線形離散時間システムを同定する。このシステム入力１１及びシステム出力１２は、図７に示すようにシステム入力印加時刻特定部５１によって、システム入力１１が有意な値を持つ複数の時刻が特定される。次に、インパルス応答変換部５２によって、システム入力印加時刻特定部５１から出力される複数の時刻に対応するシステム入力の加算値をインパルス入力振幅とし、システム入力印加時刻特定部５１から出力される複数の時刻の最大値をインパルス入力印加時刻として、インパルス入力印加時刻以降のシステム出力をインパルス入力振幅で除した信号を、動的システムのインパルス応答として出力する。
【００７０】
インパルス応答に対応する情報が動的システムを打撃加振したときのシステム入力１１及びシステム出力１２であり、システム入力１１が有意な値を持つ複数の時刻を特定するシステム入力印加時刻特定部５１と、システム入力印加時刻特定部５１から出力される複数の時刻に対応するシステム入力の加算値をインパルス入力振幅とし、システム入力印加時刻特定部５１から出力される複数の時刻の最大値をインパルス入力印加時刻として、インパルス入力印加時刻以降のシステム出力を前記インパルス入力振幅で除した信号を前記インパルス応答として出力するインパルス応答変換部５２と、を備え、同定対象とする動的システムを打撃加振した場合のシステム入力１１及びシステム出力１２並びに探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）を入力として、線形離散時間システムを同定する打撃加振によるシステム同定装置５０は、動的システムを打撃加振して得られる現実のシステム入出力から、当該システムを記述する線形離散時間システムを同定することが可能である。
【００７１】
図８は、本実施の形態のシステム同定装置において、システム次元決定部４ａの内部構成と動作を示すブロック線図である。図８において、図４と同一符号を付したものは、実施の形態１と同一または同等の構成要素である。図８に示すように、本実施の形態のシステム次元決定部４ａでは、作業者が指定したシステム次元の探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に属する次元ｎｉのそれぞれに対して、再帰的システム行列推定部３１で同定したシステム行列Ａｄ，ｎｉに基づいて、システム安定性評価部３４が線形離散時間システムの安定性を評価する。なお、図８において、「ｉ＝１」の処理、「ｉ＋＋」の処理及び「ｉ≦ａ」の判断は、システム特性推定部３２またはシステム次元推定部３３が担ってもよいし、図示しない他の構成が担ってもよい。
【００７２】
システム特性推定部３２は、システム安定性評価部３４において安定システムと判断された次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａｄ，ｎｉ，Ｂｄ，ｎｉ，Ｃｄ，ｎｉ及び直達項同定部１から出力される直達項Ｄｄに基づいて、同定した線形離散時間システムに関して現実のシステム入力１１に対応する推定システム出力を算出する。
【００７３】
システム次元決定部４ａが、探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に属する各次元に対して線形離散時間システムの安定性を評価するシステム安定性評価部３４を備え、安定システムとなる次元に対応する線形離散時間システムのシステム特性からシステム次元を決定するシステム同定装置１０は、現実の動的システムが安定システムであることが明確である場合に、安定システムに限定してシステムを記述する線形離散時間システムを同定することが可能である。
【００７４】
システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムの推定システム出力と動的システムの実際のシステム出力１２（図８における動的システムのシステム特性）との時間領域における誤差２乗和ｅｎｉ（ｎｉ：安定システムとなる次元）を算出し、図５に示すように誤差２乗和ノルム分布４１が誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ６）。
【００７５】
次に、動作について説明する。同定対象とする動的システムが、１入力Ｐ出力のｎ次元線形離散時間システムとして上記の式（１）で記述できるとする。この動的システムを打撃加振した場合のシステム入力ｕ（ｊＴｓ）は、例えば図６に示すシステム入力１１のような時間波形となり、この場合のシステム入力ｕ（ｊＴｓ）は下記の式（１５）で与えられる。
【００７６】
【数１５】
0006009106.tif 000016
【００７７】
このとき、上記の式（１）に対応するシステム出力ｙ（ｊＴｓ）、つまり図６に示すシステム出力１２は下記の式（１６）となるため、システム出力ｙ（ｊＴｓ）から直接的にインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）を求めることが困難となる。
【００７８】
【数１６】
0006009106.tif 000017
【００７９】
そこで、本実施の形態の打撃加振によるシステム同定装置５０では、打撃加振における加振時間は微小であるという前提の下、図６に示す現実のシステム入力１１が図２のような理想的なインパルス入力で近似される。具体的には、図７に示すシステム入力印加時刻特定部５１によって、システム入力ｕ（ｊＴｓ）に重畳するノイズレベルによって規定される比率しきい値とシステム入力ｕ（ｊＴｓ）の最大値とを乗じた下記の式（１７）をシステム入力しきい値１３とし、システム入力ｕ（ｊＴｓ）がシステム入力しきい値１３以上となる時刻を、システム入力ｕ（ｊＴｓ）が有意な値を持つ複数の時刻ｊ’Ｔｓ（図６の場合はｊ’＝−２，−１，０）として特定する。
【００８０】
【数１７】
0006009106.tif 000018
【００８１】
システム入力印加時刻特定部５１が、比率しきい値とシステム入力の最大値とを乗じた値をシステム入力しきい値１３とし、システム入力がシステム入力しきい値１３以上となる時刻を、システム入力が有意な値を持つ複数の時刻として特定する打撃加振によるシステム同定装置５０は、動的システムを打撃加振して得られる現実のシステム入力から、システム入力が有意な値を持つ加振時刻を正確に抽出することが可能である。
【００８２】
次にインパルス応答変換部５２によって、システム入力印加時刻特定部５１から出力される複数の時刻ｊ’Ｔｓに対応して、システム入力ｕ（ｊＴｓ）の加算値をインパルス入力振幅ａ〜として下記の式（１８）で、システム入力印加時刻特定部５１から出力される複数の時刻ｊ’Ｔｓの最大値をインパルス入力印加時刻ｊ”Ｔｓ（図６の場合はｊ”＝０）として下記の式（１９）で導出し、動的システムのインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）を、インパルス入力印加時刻ｊ”Ｔｓ以降のシステム出力ｙ（ｊＴｓ）をインパルス入力振幅ａ〜で除した信号として下記の式（２０）で算出する。
【００８３】
なお、本明細書中で、「ａ〜」は「ａ」の上に「〜（チルダ）」が配置された文字の代替表記である。
【００８４】
【数１８】
0006009106.tif 000019
【００８５】
【数１９】
0006009106.tif 000020
【００８６】
【数２０】
0006009106.tif 000021
【００８７】
本実施の形態では、動的システムを打撃加振した場合のシステム入出力に基づいて、上記の式（２０）で得られた動的システムのインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）を、図１に示すシステム同定装置１０への入力とする。本実施の形態においては、図１に示すように、システム同定装置１０では、上記の式（２０）で得られた動的システムのインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）を入力とし、直達項同定部１において動的システムを記述する上記の式（１）の線形離散時間システムの直達項Ｄｄを上記の式（４）で同定すると共に、ブロックハンケル行列生成部２によって、下記の式（２１）で与えられるブロックハンケル行列Ｈｋｌを生成する。次に特異値分解部３によって、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｈｋｌを特異値分解し、上記の式（６）で与えられるブロックハンケル行列Ｈｋｌの第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σｉ（ｉ＝１，２，３，…）を出力する。
【００８８】
【数２１】
0006009106.tif 000022
【００８９】
システム次元決定部４ａは、図８に示すように再帰的システム行列推定部３１によって、作業者が指定したシステム次元の探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に属する次元ｎｉのそれぞれに関して、上記の式（９）〜（１１）に示す再帰的手法によって対応するシステム行列Ａｄ，ｎｉ，Ｂｄ，ｎｉ，Ｃｄ，ｎｉを同定する。
【００９０】
システム安定性評価部３４は、作業者が指定したシステム次元の探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に属する次元ｎｉのそれぞれに対して、再帰的システム行列推定部３１で同定したシステム行列Ａｄ，ｎｉに基づいて、線形離散時間システムの安定性を下記の数式（２２）により評価する。
【００９１】
【数２２】
0006009106.tif 000023
【００９２】
次に、システム特性推定部３２において、システム安定性評価部３４で安定システムと判断された次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａｄ，ｎｉ，Ｂｄ，ｎｉ，Ｃｄ，ｎｉ及び直達項同定部１から出力される直達項Ｄｄに基づいて、上記の式（１）に対して現実のシステム入力ｕ（ｊＴｓ）（例えば式（１５））を適用した場合の推定システム出力ｙ＾ｎｉ（ｊＴｓ）を算出する。
【００９３】
なお、本明細書中で、「ｙ＾」は「ｙ」の上に「＾（サーカムフレックス）」が配置された文字の代替表記である。
【００９４】
システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムの推定システム出力ｙ＾ｎｉ（ｊＴｓ）と動的システムの実際のシステム出力ｙ（ｊＴｓ）（図８における動的システムのシステム特性）との時間領域における誤差２乗和ｅｎｉ（ｎｉ：安定システムとなる次元）を下記の式（２３）で算出する。
【００９５】
【数２３】
0006009106.tif 000024
【００９６】
この誤差２乗和ノルム||ｅｎｉ||を最小とする次元ｎｉが「時間領域において実際のシステム出力に最も適合する」安定なシステム次元ｎとなる。ここでは、図５に示すように誤差２乗和ノルム分布４１が上記の式（１３）で与えられる誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ６）。
【００９７】
最後に、システム行列同定部５において、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σｉ（ｉ＝１，２，３，…）と、システム次元決定部４ａから出力されるシステム次元ｎと、に基づき、線形離散時間システムの直達項Ｄｄを除くシステム行列Ａｄ，Ｂｄ，Ｃｄを上記の式（１４）で同定して、直達項同定部１で同定した直達項Ｄｄ及びシステム行列同定部５で同定したシステム行列Ａｄ，Ｂｄ，Ｃｄを、動的システムを記述する線形離散時間システムとして出力する。
【００９８】
このように、本実施の形態の打撃加振によるシステム同定装置５０により、現実のシステム入出力をインパルス応答に変換し、当該インパルス応答から算出したブロックハンケル行列Ｈｋｌの特異値σｉ（ｉ＝１，２，３，…）がなだらかな単調減少となり、有意な値を持つ特異値と同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して時間領域で一致度の高いシステム次元ｎの決定及び動的システムを記述する線形離散時間システムの同定が可能となる。
【００９９】
また、再帰的システム行列推定部３１により、現実の動的システムに対して一致度の高いシステム次元ｎを決定するための演算量を低減することが可能となる。
【０１００】
更に、動的システムを打撃加振して得られる現実のシステム入出力から、システム入力が有意な値を持つ加振時刻を正確に抽出し、これに基づいて動的システムのインパルス応答に変換することで、動的システムを打撃加振した場合の現実のシステム入出力から、当該システムを記述する線形離散時間システムを正確に同定することが可能となる。
【０１０１】
加えて、システム安定性評価部３４により、現実の動的システムが安定システムであることが明確である場合に、安定システムに限定した線形離散時間システムの同定が可能となる。
【０１０２】
なお、本実施の形態の打撃加振によるシステム同定装置５０では、線形離散時間システムのシステム特性を推定システム出力として算出し、当該システム出力と動的システムの実際のシステム出力との時間領域における誤差２乗和ノルム分布４１が誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定している。しかしながら、本発明はこれに限定されるものではなく、線形離散時間システムのシステム特性を推定周波数応答として算出し、当該周波数応答と動的システムのシステム入出力を変換して得られるインパルス応答に対して当該インパルス応答をフーリエ変換した実際の周波数応答との周波数領域における誤差２乗和に基づいて、システム次元ｎを決定してもよい。
【０１０３】
この場合、更に動的システムの実際の周波数応答に基づいて重み関数を決定し、線形離散時間システムの推定周波数応答と動的システムの実際の周波数応答との周波数領域における誤差２乗値に当該重み関数を乗じた値の加算値に基づいて、システム次元ｎを決定してもよい。
【０１０４】
本実施の形態によれば、動的システムの打撃加振によって得られる現実のシステム入出力も適用可能にすると共に、現実の動的システムが安定であることが明確である場合に、安定システムに限定して同定することができるシステム同定装置を得ることができる。
【０１０５】
実施の形態３．
本実施の形態では、インパルス応答に対応する情報が動的システムの周波数応答である場合について説明する。本実施の形態において、システム同定装置の全体構成を示すブロック線図は実施の形態１の図１と同一であり、ブロックハンケル行列Ｈｋｌの特異値σｉと次元（ｉ＝１，２，３，…）の関係を示す概略図は実施の形態１の図３と同一であり、同定した線形離散時間システムの推定周波数応答と動的システムの実際の周波数応答との周波数領域における誤差２乗和ノルム||ｅｎｉ||と次元ｎｉ（ｉ＝１，２，…，ａ）の関係を示す概略図は、実施の形態１の図５と同一である。また、本実施の形態において、システム次元決定部は実施の形態２の図８に示すシステム次元決定部４ａと同一である。
【０１０６】
図９は、本実施の形態のシステム同定装置において、動的システムの周波数応答を、当該システムのインパルス応答に変換する手順を示すブロック線図である。このとき、図９の逆フーリエ変換部６１の出力である動的システムのインパルス応答が、図１の直達項同定部１及びブロックハンケル行列生成部２の入力となる。動的システムの周波数応答を入力とし、図９に示す逆フーリエ変換部６１及び図１のシステム同定装置１０からなるシステム同定装置を周波数応答によるシステム同定装置６０とする。本実施の形態の周波数応答によるシステム同定装置６０では、同定対象とする動的システムの周波数応答に基づいて、当該システムを記述する線形離散時間システムを同定する。この動的システムの周波数応答は、図９に示すように逆フーリエ変換部６１によって逆フーリエ変換され、動的システムのインパルス応答として出力される。
【０１０７】
インパルス応答に対応する情報が動的システムの周波数応答であり、周波数応答の逆フーリエ変換により対応するインパルス応答を出力する逆フーリエ変換部６１を備え、同定対象とする動的システムの周波数応答及び探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）を入力として、線形離散時間システムを同定する周波数応答によるシステム同定装置６０は、動的システムの周波数応答から、当該システムを記述する線形離散時間システムを同定することが可能である。
【０１０８】
次に、動作について説明する。同定対象とする動的システムは、１入力Ｐ出力のｎ次元線形離散時間システムとして上記の式（１）で記述できるとする。この動的システムの（複素）周波数応答Ｈ（ｋΔｆ）（ｋ＝０，１，２，…，Ｎ−１）が与えられた場合、図９に示す逆フーリエ変換部６１によって下記の式（２４）で逆フーリエ変換し、動的システムのインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）を得る。
【０１０９】
【数２４】
0006009106.tif 000025
【０１１０】
ただし、サンプリング周期Ｔｓ＝Ｔ／Ｎであり、サンプリング周波数ｆｓ＝１／Ｔｓ＝Ｎ／Ｔであり、周波数分解能Δｆ＝１／Ｔであり、時刻ｔ＝ｊＴｓ＝ｊＴ／Ｎであり、周波数ｆ＝ｋΔｆ＝ｋ／Ｔである。
【０１１１】
本実施の形態では、動的システムの周波数応答Ｈ（ｋΔｆ）（ｋ＝０，１，２，…，Ｎ−１）に基づいて、上記の式（２４）で得られた動的システムのインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）を、システム同定装置１０への入力とする。本実施の形態においては、図１に示すように、システム同定装置１０では、上記の式（２４）で得られた動的システムのインパルス応答ｇ（ｊＴｓ）（ｊ＝０，１，２，…）を入力とし、直達項同定部１において動的システムを記述する上記の式（１）の線形離散時間システムの直達項Ｄｄを式（４）で同定すると共に、ブロックハンケル行列生成部２によって下記の式（２５）で与えられるブロックハンケル行列Ｈｋｌを生成する。
【０１１２】
【数２５】
0006009106.tif 000026
【０１１３】
次に、特異値分解部３によって、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｈｋｌを特異値分解し、上記の式（６）で与えられるブロックハンケル行列Ｈｋｌの第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σｉ（ｉ＝１，２，３，…）を出力する。
【０１１４】
システム次元決定部４ａは、図８に示すように再帰的システム行列推定部３１によって、作業者が指定したシステム次元の探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に属する次元ｎｉのそれぞれに関して、上記の式（９）〜（１１）に示す再帰的手法によって対応するシステム行列Ａｄ，ｎｉ，Ｂｄ，ｎｉ，Ｃｄ，ｎｉを同定する。
【０１１５】
システム安定性評価部３４は、作業者が指定したシステム次元の探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に属する次元ｎｉのそれぞれに対して、再帰的システム行列推定部３１で同定したシステム行列Ａｄ，ｎｉに基づいて、線形離散時間システムの安定性を上記の式（２２）で評価する。
【０１１６】
次に、システム特性推定部３２において、システム安定性評価部３４で安定システムと判断された次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａｄ，ｎｉ，Ｂｄ，ｎｉ，Ｃｄ，ｎｉ及び直達項同定部１から出力される直達項Ｄｄに基づいて、線形離散時間システムの推定周波数応答Ｈ＾ｎｉ（ｋΔｆ）（ｋ＝０，１，２，…，Ｎ−１）を算出する。
【０１１７】
なお、本明細書中で、「Ｈ＾ｎｉ」は「Ｈｎｉ」の上に「＾（サーカムフレックス）」が配置された文字の代替表記である。
【０１１８】
本実施の形態の周波数応答によるシステム同定装置６０では、「周波数領域において実際の周波数応答に最も適合する」という前提の下、システム次元推定部３３において最適なシステム次元ｎを決定する。具体的には、動的システムの実際の周波数応答Ｈ（ｋΔｆ）（ｋ＝０，１，２，…，Ｎ−１）に基づいて、例えば高ゲインかつ低周波領域に重みを持たせた下記の式（２６）のような重み関数Ｗ（ｋΔｆ）（ｋ＝０，１，２，…，Ｎ−１）を決定し、システム特性推定部３２から出力される線形離散時間システムの推定周波数応答Ｈ＾ｎｉ（ｋΔｆ）と動的システムの実際の周波数応答Ｈ（ｋΔｆ）（図８における動的システムのシステム特性）との周波数領域における誤差２乗値に対して、重み関数Ｗ（ｋΔｆ）を乗じた値の加算値ｅｎｉ（ｎｉ：安定システムとなる次元）を下記の式（２７）で算出する。
【０１１９】
【数２６】
0006009106.tif 000027
【０１２０】
【数２７】
0006009106.tif 000028
【０１２１】
この重み付き誤差２乗和ノルム||ｅｎｉ||を最小とする次元ｎｉが、「重み関数に応じて、周波数領域において実際の周波数応答に最も適合する」安定なシステム次元ｎとなる。ここでは、図５に示すように誤差２乗和ノルム分布４１が、上記の式（１３）で与えられる誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ６）。
【０１２２】
最後に、システム行列同定部５において、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ、及び特異値σｉ（ｉ＝１，２，３，…）と、システム次元決定部４ａから出力されるシステム次元ｎに基づき、線形離散時間システムの直達項Ｄｄを除くシステム行列Ａｄ，Ｂｄ，Ｃｄを、上記の式（１４）で同定して、直達項同定部１で同定した直達項Ｄｄ及びシステム行列同定部５で同定したシステム行列Ａｄ，Ｂｄ，Ｃｄを、動的システムを記述する線形離散時間システムとして出力する。
【０１２３】
このように、本実施の形態の周波数応答によるシステム同定装置６０により、現実の周波数応答をインパルス応答に変換し、当該インパルス応答から算出したブロックハンケル行列Ｈｋｌの特異値σｉ（ｉ＝１，２，３，…）がなだらかな単調減少となり、有意な値を持つ特異値と同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して周波数領域で重み関数に応じて一致度の高いシステム次元ｎの決定及び動的システムを記述する線形離散時間システムの同定が可能となる。
【０１２４】
また、再帰的システム行列推定部３１により、現実の動的システムに対して一致度の高いシステム次元ｎを決定するための演算量を低減することが可能となる。
【０１２５】
更に逆フーリエ変換に基づいて、動的システムの周波数応答をインパルス応答に変換することで、動的システムの現実の周波数応答から、当該システムを記述する線形離散時間システムを正確に同定することが可能となる。
【０１２６】
加えて、システム安定性評価部３４により、現実の動的システムが安定システムであることが明確である場合に、安定システムに限定した線形離散時間システムの同定が可能となる。
【０１２７】
なお、本実施の形態の周波数応答によるシステム同定装置６０では、線形離散時間システムのシステム特性を推定周波数応答として算出し、当該周波数応答と動的システムの実際の周波数応答との周波数領域における重み付き誤差２乗和ノルム分布４１が誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定している。しかしながら、本発明はこれに限定されるものではなく、線形離散時間システムのシステム特性を推定インパルス応答として算出し、当該インパルス応答と動的システムの周波数応答を逆フーリエ変換して得られる実際のインパルス応答との時間領域における誤差２乗和に基づいて、システム次元ｎを決定してもよい。
【０１２８】
システム次元決定部４ａが、探索範囲ｎｉ＝（ｎ１，ｎ２，…，ｎａ）（ただしｎ１＜ｎ２＜…＜ｎａ）に属する各次元に対して、線形離散時間システムのシステム特性を推定周波数応答として算出して出力するシステム特性推定部３２と、システム特性推定部３２から出力される線形離散時間システムの推定周波数応答とインパルス応答をフーリエ変換した実際の周波数応答との周波数領域における誤差２乗和ノルムがしきい値４２以下となる次元のうち最小の次元をシステム次元として決定して出力するシステム次元推定部３３と、を備えるシステム同定装置は、現実のインパルス応答から算出したブロックハンケル行列の特異値がなだらかな単調減少となり、有意な値を持つ特異値と無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、現実の動的システムに対して周波数領域で一致度の高いシステム次元の決定と、動的システムを記述する線形離散時間システムの同定が可能となる。
【０１２９】
システム次元推定部３３が、フーリエ変換した実際の周波数応答に基づいて重み関数を決定し、システム特性推定部３２から出力される線形離散時間システムの推定周波数応答と動的システムの実際の周波数応答との周波数領域における誤差２乗値に対して、重み関数を乗じた値の加算値を算出して、加算値のノルムがしきい値以下となる次元のうち最小の次元をシステム次元として決定して出力するシステム同定装置は、現実のインパルス応答から算出したブロックハンケル行列の特異値がなだらかな単調減少となり、有意な値を持つ特異値と無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、現実の動的システムに対して周波数領域で重み関数に応じて一致度の高いシステム次元の決定と、動的システムを記述する線形離散時間システムの同定が可能となる。
【０１３０】
実施の形態４．
本実施の形態では、同定対象とする動的システムがＤＣサーボモータである場合について説明する。
【０１３１】
図１０は、本発明にかかるシステム同定装置の実施の形態４の全体構成を示すブロック線図である。本実施の形態のシステム同定装置７０は、ＤＣサーボモータ７１と、周波数応答変換部７２と、周波数応答によるシステム同定装置６０と、制御系またはフィルタパラメータ設計部７３とを備える。周波数応答によるシステム同定装置６０は、実施の形態３において図１及び図９を用いて説明したものと同一である。
【０１３２】
本実施の形態では、ＤＣサーボモータ７１の入力電流［Ａｒｍｓ］として、例えばサインスイープを入力し、これをシステム入力とする。また、角速度［ｒａｄ／ｓ］をシステム出力として取得する。システム入出力は、周波数応答変換部７２により、周波数応答Ｈ（ｋΔｆ）に変換される。周波数応答変換部７２は、例えばＦＦＴ（Ｆａｓｔ  Ｆｏｕｒｉｅｒ  Ｔｒａｎｓｆｏｒｍ）などにより下記の式(２８)に基づいて周波数応答を出力する。
【０１３３】
【数２８】
0006009106.tif 000029
【０１３４】
周波数応答によるシステム同定装置６０は、この動的システムの周波数応答及びシステム次元の探索範囲を入力とし、ＤＣサーボモータ７１の入力電流［Ａｒｍｓ］から角速度［ｒａｄ／ｓ］までの動的システムを記述する線形離散時間システムを同定して制御系またはフィルタパラメータ設計部７３に出力する。このとき、システム次元の探索範囲は、ｎi＝（１，２，…，５０）のように、予測されるシステム次元に対して十分な幅をとって設定すればよい。
【０１３５】
周波数応答によるシステム同定装置６０では、現実の動的システムに対して一致度の高いシステム次元の決定と、動的システムを記述する線形離散時間システムの同定とを可能とするため、当該線形離散時間システムをサーボモータの制御系またはフィルタのパラメータ設計に利用することができる。
【産業上の利用可能性】
【０１３６】
以上のように、本発明にかかるシステム同定装置は、インパルス応答によるシステム同定装置として有用である。
【符号の説明】
【０１３７】
１  直達項同定部、２  ブロックハンケル行列生成部、３  特異値分解部、４，４ａ  システム次元決定部、５  システム行列同定部、１０，７０  システム同定装置、１１  システム入力、１２  システム出力、１３  システム入力しきい値、２１  （理想的なインパルス応答における）特異値分布、２２  （現実のインパルス応答における）特異値分布、３１  再帰的システム行列推定部、３２  システム特性推定部、３３  システム次元推定部、３４  システム安定性評価部、４１  誤差２乗和ノルム分布、４２  誤差２乗和ノルムしきい値、５０  打撃加振によるシステム同定装置、５１  システム入力印加時刻特定部、５２  インパルス応答変換部、６０  周波数応答によるシステム同定装置、６１  逆フーリエ変換部、７１  ＤＣサーボモータ、７２  周波数応答変換部、７３  制御系またはフィルタパラメータ設計部。
【図１】
0006009106.tif 000030
【図２】
0006009106.tif 000031
【図３】
0006009106.tif 000032
【図４】
0006009106.tif 000033
【図５】
0006009106.tif 000034
【図６】
0006009106.tif 000035
【図７】
0006009106.tif 000036
【図８】
0006009106.tif 000037
【図９】
0006009106.tif 000038
【図１０】
0006009106.tif 000039