WO2015118737A1

WO2015118737A1 - システム同定装置

Info

Publication number: WO2015118737A1
Application number: PCT/JP2014/079339
Authority: WO
Inventors: 光伯齊藤; 百合夏金井
Original assignee: 三菱電機株式会社
Priority date: 2014-02-07
Filing date: 2014-11-05
Publication date: 2015-08-13
Also published as: US10387116B2; DE112014006345T5; US20170010861A1; JP6009106B2; CN105960613A; JPWO2015118737A1

Abstract

　システム次元の決定から試行錯誤を排除し、常に最適なシステム次元を決定し、動的システムを記述する線形離散時間システムの同定を可能とするために、同定対象とする動的システムのインパルス応答と、システム次元の探索範囲と、を入力とし、直達項同定部と、ブロックハンケル行列生成部と、ブロックハンケル行列の特異値分解により、第１の直交行列、第２の直交行列及び特異値を出力する特異値分解部と、第１の直交行列、第２の直交行列、特異値及び探索範囲に基づいて、探索範囲に属する各次元に対して直達項以外のシステム行列を同定して、システム行列及び直達項に基づいて算出したシステム特性と実際のシステム特性との比較から、システム次元を決定するシステム次元決定部と、第１の直交行列、第２の直交行列、特異値及びシステム次元に基づいて、直達項以外のシステム行列を同定するシステム行列同定部と、を備えるシステム同定装置とする。

Description

システム同定装置

　本発明は、システム同定装置に関するものである。

　従来のインパルス応答によるシステム同定装置として、例えば非特許文献１にはＨｏ－Ｋａｌｍａｎの方法に基づくものが開示されている。この方法では、線形離散時間システム（Ａ_ｄ，Ｂ_ｄ，Ｃ_ｄ，Ｄ_ｄ）で記述される動的システムのインパルス応答（Ｇ_０，Ｇ_１，Ｇ_２，…）に基づいて、Ｇ_０からシステムの直達項Ｄ_ｄを決定すると共に、Ｇ_１，Ｇ_２，…からブロックハンケル行列Ｈ_ｋｌを生成する。次に、ブロックハンケル行列Ｈ_ｋｌを特異値分解し、有意な値を持つ特異値の個数をシステム次元として決定して、特異値分解の結果と決定したシステム次元から、拡大可観測行列Ｏ_ｋ及び拡大可到達行列Ｃ_ｌを算出する。最後に、拡大可観測行列Ｏ_ｋ及び拡大可到達行列Ｃ_ｌに基づいてシステム行列Ａ_ｄ，Ｂ_ｄ，Ｃ_ｄを算出することで、動的システムを記述する線形離散時間システム（Ａ_ｄ，Ｂ_ｄ，Ｃ_ｄ，Ｄ_ｄ）を同定している。

　また、従来のインパルス応答によるシステム同定装置の他の例として、例えば特許文献１にはプラント・モデリング装置が開示されている。このプラント・モデリング装置では、動的システムのインパルス応答（Ｇ_０，Ｇ_１，Ｇ_２，…）に対して上記Ｈｏ－Ｋａｌｍａｎの方法を適用するが、システム次元の決定方法として次の２種類を選択的に適用している。第１のシステム次元の決定方法は、特異値と当該特異値に対応する次数の関係をグラフィック端末に対数スケールで表示し、有意な値を持つ特異値の個数、つまりシステム次元を作業者が決定する方法である。第２のシステム次元の決定方法は、特異値の変化率と観測ノイズに基づく評価関数を適用し、当該評価関数を最小とする次元をシステム次元として自動的に決定する方法である。

　このようなインパルス応答によるシステム同定装置では、有意な値を持つ特異値の個数からシステム次元を決定している。

特開昭６１－２６７１０２号公報

片山徹著、「システム同定－部分空間法からのアプローチ－」、朝倉書店、２００４年２月、ｐ．１０２－１０７

　しかしながら、上記従来の技術によれば、現実のインパルス応答から算出したブロックハンケル行列の特異値は、なだらかな単調減少となる場合が多数存在し、この場合は有意な値を持つ特異値と、無視できる微小な値となる特異値との境界が不明確となる。そのため、第１のシステム次元の決定方法による従来のシステム同定装置では、システム次元の決定が作業者の判断に依存し、常に最適なシステム次元を決定しているとは限らず、またはシステム次元の決定に試行錯誤が必要になる、という問題があった。

　また、このような問題に対応する方法として、特許文献１では第２のシステム次元の決定方法を適用しているが、本手法でも評価関数の与え方によって決定されるシステム次元が変化するため、システム次元の決定が評価関数の与え方に依存し、本手法でも常に最適なシステム次元を決定しているとは限らず、またはシステム次元の決定に関して試行錯誤が必要であり、問題を解決するには至っていない。

　本発明は、上記に鑑みてなされたものであって、現実のインパルス応答から算出したブロックハンケル行列の特異値がなだらかな単調減少となり、有意な値を持つ特異値と無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、常に最適なシステム次元を決定し、動的システムを記述する線形離散時間システムの同定が可能なシステム同定装置を得ることを目的とする。

　上述した課題を解決し、目的を達成するために、本発明は、同定対象とする動的システムのインパルス応答と、指定されたシステム次元の探索範囲と、を入力とするシステム同定装置であって、前記インパルス応答から前記動的システムを記述する線形離散時間システムの直達項を同定して出力する直達項同定部と、前記インパルス応答からブロックハンケル行列を生成して出力するブロックハンケル行列生成部と、前記ブロックハンケル行列生成部から出力される前記ブロックハンケル行列の特異値分解により、前記ブロックハンケル行列の左特異ベクトルを列ベクトルとする第１の直交行列、前記ブロックハンケル行列の右特異ベクトルを列ベクトルとする第２の直交行列及び前記ブロックハンケル行列の特異値を出力する特異値分解部と、前記第１の直交行列、前記第２の直交行列、前記特異値及び前記探索範囲に基づいて、当該探索範囲に属する各次元に対する前記線形離散時間システムのシステム行列のうち前記直達項以外のシステム行列を同定して、当該システム行列及び前記直達項に基づいて算出した前記線形離散時間システムのシステム特性と前記動的システムの実際のシステム特性との比較から、システム次元を決定して出力するシステム次元決定部と、前記第１の直交行列、前記第２の直交行列、前記特異値及び前記システム次元決定部から出力される前記システム次元に基づいて、前記線形離散時間システムのシステム行列のうち前記直達項以外のシステム行列を同定するシステム行列同定部と、を備え、前記直達項同定部で同定した前記直達項及び前記システム行列同定部で同定した前記システム行列を前記線形離散時間システムとして出力することを特徴とする。

　この発明によれば、同定対象とする動的システムにおいて、現実のインパルス応答から算出したブロックハンケル行列の特異値がなだらかな単調減少となり、有意な値を持つ特異値と無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、常に最適なシステム次元の決定及び動的システムを記述する線形離散時間システムの同定が可能となる、という効果を奏する。

図１は、実施の形態１にかかるシステム同定装置の全体構成を示すブロック線図である。図２は、実施の形態１にかかるシステム同定装置において、システム入出力の時間波形を示す概略図である。図３は、実施の形態１にかかるシステム同定装置において、ブロックハンケル行列の特異値と次元の関係を示す概略図である。図４は、実施の形態１にかかるシステム同定装置において、システム次元決定部の内部構成と動作を示すブロック線図である。図５は、実施の形態１にかかるシステム同定装置において、同定した線形離散時間システムの推定インパルス応答とシステムの実際のインパルス応答との誤差２乗和ノルムと、次元との関係を示す概略図である。図６は、実施の形態２にかかるシステム同定装置において、動的システムを打撃加振した場合のシステム入出力の時間波形を示す概略図である。図７は、実施の形態２にかかるシステム同定装置において、動的システムを打撃加振した場合のシステム入出力を、当該システムのインパルス応答に変換する手順を示すブロック線図である。図８は、実施の形態２にかかるシステム同定装置において、システム次元決定部の内部構成と動作を示すブロック線図である。図９は、実施の形態３にかかるシステム同定装置において、動的システムの周波数応答を、当該システムのインパルス応答に変換する手順を示すブロック線図である。図１０は、実施の形態４にかかるシステム同定装置の全体構成を示すブロック線図である。

　以下に、本発明にかかるシステム同定装置の実施の形態を図面に基づいて詳細に説明する。なお、この実施の形態によりこの発明が限定されるものではない。

実施の形態１．
　図１は、本発明にかかるシステム同定装置の実施の形態１の全体構成を示すブロック線図である。図２は、本実施の形態にかかるシステム同定装置において、システム入出力の時間波形を示す概略図である。図１に示すシステム同定装置１０では、同定対象とする動的システムのインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）を入力とする。このインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）は、図２に示すように、時刻０では振幅が１となり、当該時刻以外では振幅が０となる理想的なインパルス入力で動的システムへのシステム入力１１を構成した場合における時刻０以降のシステム出力１２で与えられる。

　システム同定装置１０では、こうして得られた動的システムのインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）から、直達項同定部１によって動的システムを記述する線形離散時間システムの直達項Ｄ_ｄを同定し、更にブロックハンケル行列生成部２によってブロックハンケル行列Ｈ_ｋｌを生成する。

　特異値分解部３では、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｈ_ｋｌを特異値分解し、ブロックハンケル行列Ｈ_ｋｌの左特異ベクトルを列ベクトルとする第１の直交行列Ｕ、ブロックハンケル行列Ｈ_ｋｌの右特異ベクトルを列ベクトルとする第２の直交行列Ｖ及びブロックハンケル行列Ｈ_ｋｌの特異値σ_ｉ（ｉ＝１，２，３，…）を出力する。

　システム次元決定部４では、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σ_ｉ（ｉ＝１，２，３，…）と、作業者が指定したシステム次元の探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）に基づいて、当該探索範囲に属する次元ｎ_ｉ（ｉ＝１，２，…，ａ）のそれぞれに対して動的システムを記述する線形離散時間システムの直達項Ｄ_ｄを除くシステム行列を同定する。更に、当該システム行列及び直達項同定部１から出力される直達項Ｄ_ｄに基づいて、当該探索範囲に属する次元ｎ_ｉ（ｉ＝１，２，…，ａ）のそれぞれに対して線形離散時間システムの推定インパルス応答を算出し、動的システムの実際のインパルス応答（図１における動的システムのシステム特性）との比較から、システム次元ｎを決定する。

　システム行列同定部５では、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σ_ｉ（ｉ＝１，２，３，…）と、システム次元決定部４から出力されるシステム次元ｎに基づいて、線形離散時間システムの直達項Ｄ_ｄを除くシステム行列Ａ_ｄ，Ｂ_ｄ，Ｃ_ｄを同定する。

　システム同定装置１０では、直達項同定部１で同定した直達項Ｄ_ｄ及びシステム行列同定部５で同定したシステム行列Ａ_ｄ，Ｂ_ｄ，Ｃ_ｄを、動的システムを記述する線形離散時間システムとして最終的に出力する。

　図３は、本実施の形態にかかるシステム同定装置１０において、ブロックハンケル行列Ｈ_ｋｌの特異値σ_ｉと次元（ｉ＝１，２，３，…）の関係を示す概略図である。図４は、本実施の形態にかかるシステム同定装置１０において、システム次元決定部４の内部構成と動作を示すブロック線図である。図５は、本実施の形態にかかるシステム同定装置１０において、同定した線形離散時間システムの推定インパルス応答と動的システムの実際のインパルス応答との時間領域における誤差２乗和ノルム||ｅ_ｎｉ||と、次元ｎ_ｉ（ｉ＝１，２，…，ａ）との関係を示す概略図である。なお、図４において、「ｉ＝１」の処理、「ｉ＋＋」の処理及び「ｉ≦ａ」の判断は、システム特性推定部３２またはシステム次元推定部３３が担ってもよいし、図示しない他の構成が担ってもよい。

　図３に示すように、動的システムのインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）から生成したブロックハンケル行列Ｈ_ｋｌの特異値σ_ｉは、理想的なインパルス応答では次元（ｉ＝１，２，３，…）に対して、例えば特異値分布２１に示す関係となる。

　この場合、有意な値を持つ特異値の個数を明確に規定することができ、当該個数が動的システムのシステム次元ｎに対応する（図３の場合はシステム次元ｎ＝４）。一方、観測ノイズ等の影響を受ける現実のインパルス応答に基づいて算出した特異値σ_ｉは、次元（ｉ＝１，２，３，…）に対して例えば特異値分布２２に示す関係となるため、有意な値を持つ特異値と無視できる微小な値となる特異値との境界が不明確となり、常に最適なシステム次元ｎを決定しているとは限らず、またはシステム次元ｎの決定に関して試行錯誤が必要になるという問題がある。

　そこで、システム次元決定部４は、図４に示すように再帰的システム行列推定部３１により、作業者が指定したシステム次元の探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）に属する第１の次元ｎ_ｉに関して、第１の次元ｎ_ｉよりも１段階低い第２の次元ｎ_ｉ－１に対応したシステム行列Ａ_{ｄ，ｎｉ－１}，Ｂ_{ｄ，ｎｉ－１}，Ｃ_{ｄ，ｎｉ－１}の同定結果と、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ、及び特異値σ_ｉ（ｉ＝１，２，３，…）のうち、それぞれ第２の次元ｎ_ｉ－１よりも大きく、第１の次元ｎ_ｉ以下となる左特異ベクトルｕ_ｊ、右特異ベクトルｖ_ｊ及び特異値σ_ｊ（ｊ＝ｎ_ｉ－１＋１，ｎ_ｉ－１＋２，…，ｎ_ｉ）を用いて、再帰的手法によって第１の次元ｎ_ｉに対応するシステム行列Ａ_ｄ，ｎｉ，Ｂ_ｄ，ｎｉ，Ｃ_ｄ，ｎｉを同定する。

　システム次元決定部４が、探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）の中で第１の次元ｎ_ｉよりも１段階低い第２の次元ｎ_ｉ－１に対応したシステム行列Ａ_{ｄ，ｎｉ－１}，Ｂ_{ｄ，ｎｉ－１}，Ｃ_{ｄ，ｎｉ－１}の同定結果と、第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σ_ｉ（ｉ＝１，２，３，…）からそれぞれ第２の次元ｎ_ｉ－１よりも大きく第１の次元ｎ_ｉ以下となる左特異ベクトルｕ_ｊ、右特異ベクトルｖ_ｊ及び特異値σ_ｊ（ｊ＝ｎ_ｉ－１＋１，ｎ_ｉ－１＋２，…，ｎ_ｉ）と、を用いて、再帰的手法により第１の次元ｎ_ｉに対応するシステム行列Ａ_ｄ，ｎｉ，Ｂ_ｄ，ｎｉ，Ｃ_ｄ，ｎｉを同定する再帰的システム行列推定部３１を備えるシステム同定装置１０は、現実の動的システムに対して一致度の高いシステム次元を決定するための演算量を低減することが可能である。

　次に、システム特性推定部３２は、システム次元の探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）に属するそれぞれの次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａ_ｄ，ｎｉ，Ｂ_ｄ，ｎｉ，Ｃ_ｄ，ｎｉ及び直達項同定部１から出力される直達項Ｄ_ｄに基づいて、同定した線形離散時間システムの推定インパルス応答を算出する。そして、ｉに１を加算し、ｉ≦ａの場合には再帰的システム行列推定部３１により、システム行列Ａ_ｄ，ｎｉ，Ｂ_ｄ，ｎｉ，Ｃ_ｄ，ｎｉを同定し、ｉ＞ａの場合にはシステム次元推定部３３の処理に進む。

　システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムの推定インパルス応答と動的システムの実際のインパルス応答（図４における動的システムのシステム特性）との時間領域における誤差２乗和ｅ_ｎｉ（ｉ＝１，２，…，ａ）を算出し、図５に示すように誤差２乗和ノルム分布４１が誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ_６）。

　すなわち、システム次元決定部４が、探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）に属する各次元に対して、線形離散時間システムのシステム特性を推定インパルス応答として算出して出力するシステム特性推定部３２と、システム特性推定部３２から出力される線形離散時間システムの推定インパルス応答と動的システムの実際のインパルス応答との時間領域における誤差２乗和ノルムがしきい値以下となる次元のうち最小の次元をシステム次元として決定して出力するシステム次元推定部３３と、を備えるシステム同定装置１０は、現実のインパルス応答から算出したブロックハンケル行列の特異値がなだらかな単調減少となり、有意な値を持つ特異値と無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、現実の動的システムに対して時間領域で一致度の高いシステム次元の決定と、動的システムを記述する線形離散時間システムの同定が可能となる。

　次に、動作について説明する。同定対象とする動的システムが下記の式（１）に示す１入力Ｐ出力のｎ次元線形離散時間システムで記述できるとする。上記の動的システムへのシステム入力ｕ（ｊＴ_ｓ）を、図２に示すシステム入力１１、つまり理想的なインパルス入力となるように下記の式（２）で構成すると、下記の式（１）に対応するシステム出力ｙ（ｊＴ_ｓ）、つまり図２に示すシステム出力１２は下記の式（３）となり、動的システムのインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）は、時刻０以降のシステム出力１２で与えられる。

　ただし、状態ベクトルｘ∈Ｒ^ｎであり、システム入力ｕ∈Ｒであり、システム出力ｙ∈Ｒ^Ｐであり、システム行列Ａ_ｄ∈Ｒ^ｎ×ｎ，Ｂ_ｄ∈Ｒ^ｎ，Ｃ_ｄ∈Ｒ^Ｐ×ｎ，Ｄ_ｄ∈Ｒ^Ｐであり、インパルス応答列ｇ（ｊＴｓ）（ｊ＝０，１，２，…）である。

　図１に示すように、本実施の形態にかかるシステム同定装置１０は、上記で得られた動的システムのインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）を入力とし、直達項同定部１において上記の式（１）にて動的システムを記述する線形離散時間システムの直達項Ｄ_ｄを下記の式（４）で同定し、更にブロックハンケル行列生成部２によって下記の式（５）で与えられるブロックハンケル行列Ｈ_ｋｌを生成する。

　次に、特異値分解部３によって、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｈ_ｋｌを特異値分解し、下記の式（６）で与えられるブロックハンケル行列Ｈ_ｋｌの左特異ベクトルｕ_ｊを列ベクトルとする第１の直交行列Ｕ、ブロックハンケル行列Ｈ_ｋｌの右特異ベクトルｖ_ｊを列ベクトルとする第２の直交行列Ｖ及びブロックハンケル行列Ｈ_ｋｌの特異値σ_ｉ（ｉ＝１，２，３，…）を出力する。

　ただし、第１の直交行列Ｕ＝［ｕ_１　ｕ_２　…　ｕ_ｋＰ］∈Ｒ^{ｋＰ×ｋＰ}であり、第２の直交行列Ｖ＝［ｖ_１　ｖ_２　…　ｕ_ｌ］∈Ｒ^ｌ×ｌであり、ブロックハンケル行列Ｈｋｌの特異値σ_１≧σ_２≧…≧σ_ｎ≧σ_ｎ＋１≧σ_ｎ＋２≧…であり、Σは下記の式（７）で表される。

　対象とする動的システムのシステム次元ｎは、ブロックハンケル行列Ｈ_ｋｌの特異値において有意な値を持つものはｎ個であり、それらに比してｎ＋１個目以降は十分小さく、無視できる微小な値であるとする下記の式（８）の関係に基づいて決定することができる。

　図３に示すように、動的システムのインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）から生成したブロックハンケル行列Ｈ_ｋｌの特異値σ_ｉは、理想的なインパルス応答では次元（ｉ＝１，２，３，…）に対して例えば特異値分布２１に示す関係となる。この場合、有意な値を持つ特異値の個数を明確に規定することができ、当該個数から動的システムのシステム次元ｎを決定することができる（図３の場合はシステム次元ｎ＝４）。

　一方、観測ノイズ等の影響を受ける現実のインパルス応答に基づいて算出した特異値σ_ｉは、次元（ｉ＝１，２，３，…）に対して例えば特異値分布２２に示す関係となるため、有意な値を持つ特異値と、無視できる微小な値となる特異値との境界が不明確となる。したがって、従来の手法では常に最適なシステム次元ｎを決定しているとは限らず、またはシステム次元ｎの決定に関して試行錯誤が必要になるという課題があった。そこで、本実施の形態のシステム同定装置１０では、「時間領域において実際のインパルス応答に最も適合する」という前提の下、システム次元決定部４において最適なシステム次元ｎを決定する。

　システム次元決定部４は、図１に示すように特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σ_ｉ（ｉ＝１，２，３，…）と、作業者が指定したシステム次元の探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）に基づいて、当該探索範囲に属する次元ｎ_ｉ（ｉ＝１，２，…，ａ）のそれぞれに対して動的システムを記述する線形離散時間システムの直達項Ｄ_ｄを除くシステム行列を同定する。

　更に、当該システム行列及び直達項同定部１から出力される直達項Ｄ_ｄに基づいて、当該探索範囲に属する次元ｎ_ｉ（ｉ＝１，２，…，ａ）のそれぞれに対して同定した線形離散時間システムの推定インパルス応答を算出し、動的システムの実際のインパルス応答（図１における動的システムのシステム特性）との比較から、システム次元ｎを決定する。具体的には図４に示すように、再帰的システム行列推定部３１によって、作業者が指定したシステム次元の探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）に属する第１の次元ｎ_ｉに関して、第１の次元ｎ_ｉよりも１段階低い第２の次元ｎ_ｉ－１に対応したシステム行列Ａ_{ｄ，ｎｉ－１}，Ｂ_{ｄ，ｎｉ－１}，Ｃ_{ｄ，ｎｉ－１}の同定結果と、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σ_ｉ（ｉ＝１，２，３，…）のうち、それぞれ第２の次元ｎ_ｉ－１よりも大きく、第１の次元ｎ_ｉ以下となる左特異ベクトルｕ_ｊ、右特異ベクトルｖ_ｊ及び特異値σ_ｊ（ｊ＝ｎ_ｉ－１＋１，ｎ_ｉ－１＋２，…，ｎ_ｉ）とを用いて、次式に示す再帰的手法によって第１の次元ｎ_ｉに対応するシステム行列Ａ_ｄ，ｎｉ，Ｂ_ｄ，ｎｉ，Ｃ_ｄ，ｎｉを同定する。まず、第１の次元ｎ_ｉに対応する拡大可観測行列は下記の式（９）にて表される。

　第１の次元ｎ_ｉに対応する拡大可到達行列は下記の式（１０）にて表される。

　そして、第１の次元ｎ_ｉに対応するシステム行列Ａ_ｄ，ｎｉ，Ｂ_ｄ，ｎｉ，Ｃ_ｄ，ｎｉは下記の式（１１）にて表される。

　次にシステム特性推定部３２において、システム次元の探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）に属する次元ｎ_ｉのそれぞれに対して、再帰的システム行列推定部３１から出力されるシステム行列Ａ_ｄ，ｎｉ，Ｂ_ｄ，ｎｉ，Ｃ_ｄ，ｎｉ及び直達項同定部１から出力される直達項Ｄ_ｄに基づいて、上記の式（１）に対するシステム入力ｕ（ｊＴ_ｓ）として上記の式（２）を適用することで、同定した線形離散時間システムの推定インパルス応答ｇ^＾ _ｎｉ（ｊＴ_ｓ）（ｊ＝０，１，２，…）を算出する。

　なお、本明細書中で、「ｇ^＾」は「ｇ」の上に「＾（サーカムフレックス）」が配置された文字の代替表記である。

　システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムの推定インパルス応答ｇ^＾ _ｎｉ（ｊＴ_ｓ）（ｊ＝０，１，２，…）と動的システムの実際のインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）（図４における動的システムのシステム特性）との時間領域における誤差２乗和ｅ_ｎｉ（ｉ＝１，２，…，ａ）を下記の式（１２）で算出する。

　この誤差２乗和ノルム||ｅ_ｎｉ||を最小とする次元ｎ_ｉが「時間領域において実際のインパルス応答に最も適合する」システム次元ｎとなる。一方、実際のノルム||ｅ_ｎｉ||は、観測ノイズが白色性であれば、そのノイズレベルによらず次元ｎ_ｉの増加に伴って単調減少し、図５に示すように所定の次元（図５の場合はｎ_６）以上でほぼ一定値となる。そこで、ここでは、システム次元ｎの推定値が必要以上に高次元となることを避けるため、下記の式（１３）で与えられる誤差２乗和ノルムしきい値４２を規定し、誤差２乗和ノルム分布４１が誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ_６）。

　最後に、システム行列同定部５において、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σ_ｉ（ｉ＝１，２，３，…）と、システム次元決定部４から出力されるシステム次元ｎと、に基づき、線形離散時間システムの直達項Ｄ_ｄを除くシステム行列Ａ_ｄ，Ｂ_ｄ，Ｃ_ｄを下記の式（１４）で同定する。

　システム同定装置１０は、システム行列同定部５で同定したシステム行列Ａ_ｄ，Ｂ_ｄ，Ｃ_ｄ及び直達項同定部１で同定した直達項Ｄ_ｄを、動的システムを記述する線形離散時間システムとして最終的に出力する。

　このように、本実施の形態のシステム同定装置１０により、現実のインパルス応答から算出したブロックハンケル行列Ｈ_ｋｌの特異値σ_ｉ（ｉ＝１，２，３，…）がなだらかな単調減少となり、有意な値を持つ特異値と同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して時間領域で一致度の高いシステム次元ｎの決定及び動的システムを記述する線形離散時間システムの同定が可能となる。

　また、再帰的システム行列推定部３１により、現実の動的システムに対して一致度の高いシステム次元ｎを決定するための演算量を低減することが可能となる。

　なお、本実施の形態のシステム同定装置１０では、線形離散時間システムのシステム特性を推定インパルス応答として算出し、当該インパルス応答と動的システムの実際のインパルス応答との時間領域における誤差２乗和ノルム分布４１が誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定している。しかしながら、本発明はこれに限定されず、線形離散時間システムのシステム特性を推定周波数応答として算出し、当該周波数応答と動的システムのインパルス応答をフーリエ変換した実際の周波数応答との周波数領域における誤差２乗和に基づいて、システム次元ｎを決定してもよい。この場合、更に動的システムの実際の周波数応答に基づいて重み関数を決定し、線形離散時間システムの推定周波数応答と動的システムの実際の周波数応答との周波数領域における誤差２乗値に当該重み関数を乗じた値の加算値に基づいて、システム次元ｎを決定してもよい。

実施の形態２．
　本実施の形態では、動的システムの打撃加振によって得られる現実のシステム入出力も適用可能にすると共に、現実の動的システムが安定であることが明確である場合に安定システムに限定して同定することができるシステム同定装置について説明する。

　動的システムのインパルス応答を得るための方法としては、対象とする動的システムをインパルスハンマによって打撃加振し、そのときのシステム出力を計測する方法が最も一般的である。しかしながら、打撃加振によって与えられる現実のシステム入力は正弦半波状となるため、真のインパルス入力とは異なる。したがって、当該入力によって得られるシステム出力も、動的システムの真のインパルス応答とは異なるものとなるため、当該システム出力をそのままＨｏ－Ｋａｌｍａｎの方法に適用することはできない、という問題があった。

　加えて、従来のインパルス応答によるシステム同定装置では、同定結果として得られる線形離散時間システム（Ａ_ｄ，Ｂ_ｄ，Ｃ_ｄ，Ｄ_ｄ）の安定性は全く考慮されていないため、現実の動的システムが安定であるにも関わらず、不安定システムとして同定される場合がある、という問題もあった。

　本実施の形態では、上記問題に鑑み、動的システムの打撃加振によって得られる現実のシステム入出力も適用可能にすると共に、現実の動的システムが安定であることが明確である場合に、安定システムに限定して同定することができるシステム同定装置を得ることを目的とする。

　本実施の形態において、システム同定装置の全体構成を示すブロック線図は実施の形態１の図１と同一であり、ブロックハンケル行列Ｈ_ｋｌの特異値σ_ｉと次元（ｉ＝１，２，３，…）の関係を示す概略図は実施の形態１の図３と同一であり、同定した線形離散時間システムの推定システム出力と動的システムの実際のシステム出力との時間領域における誤差２乗和ノルム||ｅ_ｎｉ||と、次元ｎ_ｉ（ｉ＝１，２，…，ａ）との関係を示す概略図は実施の形態１の図５と同一である。

　図６は、本実施の形態のシステム同定装置において、動的システムを打撃加振した場合のシステム入出力の時間波形を示す概略図である。図７は、本実施の形態のシステム同定装置において、動的システムを打撃加振した場合のシステム入出力を、当該システムのインパルス応答に変換する手順を示すブロック線図である。このとき、図７のインパルス応答変換部５２の出力である動的システムのインパルス応答が、図１の直達項同定部１及びブロックハンケル行列生成部２の入力となる。同定対象とする動的システムを打撃加振した場合のシステム入力１１及びシステム出力１２を入力とし、図７に示すシステム入力印加時刻特定部５１及びインパルス応答変換部５２並びに図１のシステム同定装置１０全体を併せて、打撃加振によるシステム同定装置５０とする。

　図６に示すように、本実施の形態の打撃加振によるシステム同定装置５０では、同定対象とする動的システムを打撃加振した場合のシステム入力１１及びシステム出力１２に基づいて、当該システムを記述する線形離散時間システムを同定する。このシステム入力１１及びシステム出力１２は、図７に示すようにシステム入力印加時刻特定部５１によって、システム入力１１が有意な値を持つ複数の時刻が特定される。次に、インパルス応答変換部５２によって、システム入力印加時刻特定部５１から出力される複数の時刻に対応するシステム入力の加算値をインパルス入力振幅とし、システム入力印加時刻特定部５１から出力される複数の時刻の最大値をインパルス入力印加時刻として、インパルス入力印加時刻以降のシステム出力をインパルス入力振幅で除した信号を、動的システムのインパルス応答として出力する。

　インパルス応答に対応する情報が動的システムを打撃加振したときのシステム入力１１及びシステム出力１２であり、システム入力１１が有意な値を持つ複数の時刻を特定するシステム入力印加時刻特定部５１と、システム入力印加時刻特定部５１から出力される複数の時刻に対応するシステム入力の加算値をインパルス入力振幅とし、システム入力印加時刻特定部５１から出力される複数の時刻の最大値をインパルス入力印加時刻として、インパルス入力印加時刻以降のシステム出力を前記インパルス入力振幅で除した信号を前記インパルス応答として出力するインパルス応答変換部５２と、を備え、同定対象とする動的システムを打撃加振した場合のシステム入力１１及びシステム出力１２並びに探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）を入力として、線形離散時間システムを同定する打撃加振によるシステム同定装置５０は、動的システムを打撃加振して得られる現実のシステム入出力から、当該システムを記述する線形離散時間システムを同定することが可能である。

　図８は、本実施の形態のシステム同定装置において、システム次元決定部４ａの内部構成と動作を示すブロック線図である。図８において、図４と同一符号を付したものは、実施の形態１と同一または同等の構成要素である。図８に示すように、本実施の形態のシステム次元決定部４ａでは、作業者が指定したシステム次元の探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）に属する次元ｎ_ｉのそれぞれに対して、再帰的システム行列推定部３１で同定したシステム行列Ａ_ｄ，ｎｉに基づいて、システム安定性評価部３４が線形離散時間システムの安定性を評価する。なお、図８において、「ｉ＝１」の処理、「ｉ＋＋」の処理及び「ｉ≦ａ」の判断は、システム特性推定部３２またはシステム次元推定部３３が担ってもよいし、図示しない他の構成が担ってもよい。

　システム特性推定部３２は、システム安定性評価部３４において安定システムと判断された次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａ_ｄ，ｎｉ，Ｂ_ｄ，ｎｉ，Ｃ_ｄ，ｎｉ及び直達項同定部１から出力される直達項Ｄ_ｄに基づいて、同定した線形離散時間システムに関して現実のシステム入力１１に対応する推定システム出力を算出する。

　システム次元決定部４ａが、探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）に属する各次元に対して線形離散時間システムの安定性を評価するシステム安定性評価部３４を備え、安定システムとなる次元に対応する線形離散時間システムのシステム特性からシステム次元を決定するシステム同定装置１０は、現実の動的システムが安定システムであることが明確である場合に、安定システムに限定してシステムを記述する線形離散時間システムを同定することが可能である。

　システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムの推定システム出力と動的システムの実際のシステム出力１２（図８における動的システムのシステム特性）との時間領域における誤差２乗和ｅ_ｎｉ（ｎ_ｉ：安定システムとなる次元）を算出し、図５に示すように誤差２乗和ノルム分布４１が誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ_６）。

　次に、動作について説明する。同定対象とする動的システムが、１入力Ｐ出力のｎ次元線形離散時間システムとして上記の式（１）で記述できるとする。この動的システムを打撃加振した場合のシステム入力ｕ（ｊＴ_ｓ）は、例えば図６に示すシステム入力１１のような時間波形となり、この場合のシステム入力ｕ（ｊＴ_ｓ）は下記の式（１５）で与えられる。

　このとき、上記の式（１）に対応するシステム出力ｙ（ｊＴ_ｓ）、つまり図６に示すシステム出力１２は下記の式（１６）となるため、システム出力ｙ（ｊＴ_ｓ）から直接的にインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）を求めることが困難となる。

　そこで、本実施の形態の打撃加振によるシステム同定装置５０では、打撃加振における加振時間は微小であるという前提の下、図６に示す現実のシステム入力１１が図２のような理想的なインパルス入力で近似される。具体的には、図７に示すシステム入力印加時刻特定部５１によって、システム入力ｕ（ｊＴ_ｓ）に重畳するノイズレベルによって規定される比率しきい値とシステム入力ｕ（ｊＴ_ｓ）の最大値とを乗じた下記の式（１７）をシステム入力しきい値１３とし、システム入力ｕ（ｊＴ_ｓ）がシステム入力しきい値１３以上となる時刻を、システム入力ｕ（ｊＴ_ｓ）が有意な値を持つ複数の時刻ｊ’Ｔ_ｓ（図６の場合はｊ’＝－２，－１，０）として特定する。

　システム入力印加時刻特定部５１が、比率しきい値とシステム入力の最大値とを乗じた値をシステム入力しきい値１３とし、システム入力がシステム入力しきい値１３以上となる時刻を、システム入力が有意な値を持つ複数の時刻として特定する打撃加振によるシステム同定装置５０は、動的システムを打撃加振して得られる現実のシステム入力から、システム入力が有意な値を持つ加振時刻を正確に抽出することが可能である。

　次にインパルス応答変換部５２によって、システム入力印加時刻特定部５１から出力される複数の時刻ｊ’Ｔ_ｓに対応して、システム入力ｕ（ｊＴ_ｓ）の加算値をインパルス入力振幅ａ^～として下記の式（１８）で、システム入力印加時刻特定部５１から出力される複数の時刻ｊ’Ｔ_ｓの最大値をインパルス入力印加時刻ｊ”Ｔ_ｓ（図６の場合はｊ”＝０）として下記の式（１９）で導出し、動的システムのインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）を、インパルス入力印加時刻ｊ”Ｔ_ｓ以降のシステム出力ｙ（ｊＴ_ｓ）をインパルス入力振幅ａ^～で除した信号として下記の式（２０）で算出する。

　なお、本明細書中で、「ａ^～」は「ａ」の上に「^～（チルダ）」が配置された文字の代替表記である。

　本実施の形態では、動的システムを打撃加振した場合のシステム入出力に基づいて、上記の式（２０）で得られた動的システムのインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）を、図１に示すシステム同定装置１０への入力とする。本実施の形態においては、図１に示すように、システム同定装置１０では、上記の式（２０）で得られた動的システムのインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）を入力とし、直達項同定部１において動的システムを記述する上記の式（１）の線形離散時間システムの直達項Ｄ_ｄを上記の式（４）で同定すると共に、ブロックハンケル行列生成部２によって、下記の式（２１）で与えられるブロックハンケル行列Ｈ_ｋｌを生成する。次に特異値分解部３によって、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｈ_ｋｌを特異値分解し、上記の式（６）で与えられるブロックハンケル行列Ｈ_ｋｌの第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σ_ｉ（ｉ＝１，２，３，…）を出力する。

　システム次元決定部４ａは、図８に示すように再帰的システム行列推定部３１によって、作業者が指定したシステム次元の探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）に属する次元ｎ_ｉのそれぞれに関して、上記の式（９）～（１１）に示す再帰的手法によって対応するシステム行列Ａ_ｄ，ｎｉ，Ｂ_ｄ，ｎｉ，Ｃ_ｄ，ｎｉを同定する。

　システム安定性評価部３４は、作業者が指定したシステム次元の探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）に属する次元ｎ_ｉのそれぞれに対して、再帰的システム行列推定部３１で同定したシステム行列Ａ_ｄ，ｎｉに基づいて、線形離散時間システムの安定性を下記の数式（２２）により評価する。

　次に、システム特性推定部３２において、システム安定性評価部３４で安定システムと判断された次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａ_ｄ，ｎｉ，Ｂ_ｄ，ｎｉ，Ｃ_ｄ，ｎｉ及び直達項同定部１から出力される直達項Ｄ_ｄに基づいて、上記の式（１）に対して現実のシステム入力ｕ（ｊＴ_ｓ）（例えば式（１５））を適用した場合の推定システム出力ｙ^＾ _ｎｉ（ｊＴ_ｓ）を算出する。

　なお、本明細書中で、「ｙ^＾」は「ｙ」の上に「＾（サーカムフレックス）」が配置された文字の代替表記である。

　システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムの推定システム出力ｙ^＾ _ｎｉ（ｊＴ_ｓ）と動的システムの実際のシステム出力ｙ（ｊＴ_ｓ）（図８における動的システムのシステム特性）との時間領域における誤差２乗和ｅ_ｎｉ（ｎ_ｉ：安定システムとなる次元）を下記の式（２３）で算出する。

　この誤差２乗和ノルム||ｅ_ｎｉ||を最小とする次元ｎ_ｉが「時間領域において実際のシステム出力に最も適合する」安定なシステム次元ｎとなる。ここでは、図５に示すように誤差２乗和ノルム分布４１が上記の式（１３）で与えられる誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ_６）。

　最後に、システム行列同定部５において、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σ_ｉ（ｉ＝１，２，３，…）と、システム次元決定部４ａから出力されるシステム次元ｎと、に基づき、線形離散時間システムの直達項Ｄ_ｄを除くシステム行列Ａ_ｄ，Ｂ_ｄ，Ｃ_ｄを上記の式（１４）で同定して、直達項同定部１で同定した直達項Ｄ_ｄ及びシステム行列同定部５で同定したシステム行列Ａ_ｄ，Ｂ_ｄ，Ｃ_ｄを、動的システムを記述する線形離散時間システムとして出力する。

　このように、本実施の形態の打撃加振によるシステム同定装置５０により、現実のシステム入出力をインパルス応答に変換し、当該インパルス応答から算出したブロックハンケル行列Ｈ_ｋｌの特異値σ_ｉ（ｉ＝１，２，３，…）がなだらかな単調減少となり、有意な値を持つ特異値と同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して時間領域で一致度の高いシステム次元ｎの決定及び動的システムを記述する線形離散時間システムの同定が可能となる。

　更に、動的システムを打撃加振して得られる現実のシステム入出力から、システム入力が有意な値を持つ加振時刻を正確に抽出し、これに基づいて動的システムのインパルス応答に変換することで、動的システムを打撃加振した場合の現実のシステム入出力から、当該システムを記述する線形離散時間システムを正確に同定することが可能となる。

　加えて、システム安定性評価部３４により、現実の動的システムが安定システムであることが明確である場合に、安定システムに限定した線形離散時間システムの同定が可能となる。

　なお、本実施の形態の打撃加振によるシステム同定装置５０では、線形離散時間システムのシステム特性を推定システム出力として算出し、当該システム出力と動的システムの実際のシステム出力との時間領域における誤差２乗和ノルム分布４１が誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定している。しかしながら、本発明はこれに限定されるものではなく、線形離散時間システムのシステム特性を推定周波数応答として算出し、当該周波数応答と動的システムのシステム入出力を変換して得られるインパルス応答に対して当該インパルス応答をフーリエ変換した実際の周波数応答との周波数領域における誤差２乗和に基づいて、システム次元ｎを決定してもよい。

　この場合、更に動的システムの実際の周波数応答に基づいて重み関数を決定し、線形離散時間システムの推定周波数応答と動的システムの実際の周波数応答との周波数領域における誤差２乗値に当該重み関数を乗じた値の加算値に基づいて、システム次元ｎを決定してもよい。

　本実施の形態によれば、動的システムの打撃加振によって得られる現実のシステム入出力も適用可能にすると共に、現実の動的システムが安定であることが明確である場合に、安定システムに限定して同定することができるシステム同定装置を得ることができる。

実施の形態３．
　本実施の形態では、インパルス応答に対応する情報が動的システムの周波数応答である場合について説明する。本実施の形態において、システム同定装置の全体構成を示すブロック線図は実施の形態１の図１と同一であり、ブロックハンケル行列Ｈ_ｋｌの特異値σ_ｉと次元（ｉ＝１，２，３，…）の関係を示す概略図は実施の形態１の図３と同一であり、同定した線形離散時間システムの推定周波数応答と動的システムの実際の周波数応答との周波数領域における誤差２乗和ノルム||ｅ_ｎｉ||と次元ｎ_ｉ（ｉ＝１，２，…，ａ）の関係を示す概略図は、実施の形態１の図５と同一である。また、本実施の形態において、システム次元決定部は実施の形態２の図８に示すシステム次元決定部４ａと同一である。

　図９は、本実施の形態のシステム同定装置において、動的システムの周波数応答を、当該システムのインパルス応答に変換する手順を示すブロック線図である。このとき、図９の逆フーリエ変換部６１の出力である動的システムのインパルス応答が、図１の直達項同定部１及びブロックハンケル行列生成部２の入力となる。動的システムの周波数応答を入力とし、図９に示す逆フーリエ変換部６１及び図１のシステム同定装置１０からなるシステム同定装置を周波数応答によるシステム同定装置６０とする。本実施の形態の周波数応答によるシステム同定装置６０では、同定対象とする動的システムの周波数応答に基づいて、当該システムを記述する線形離散時間システムを同定する。この動的システムの周波数応答は、図９に示すように逆フーリエ変換部６１によって逆フーリエ変換され、動的システムのインパルス応答として出力される。

　インパルス応答に対応する情報が動的システムの周波数応答であり、周波数応答の逆フーリエ変換により対応するインパルス応答を出力する逆フーリエ変換部６１を備え、同定対象とする動的システムの周波数応答及び探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）を入力として、線形離散時間システムを同定する周波数応答によるシステム同定装置６０は、動的システムの周波数応答から、当該システムを記述する線形離散時間システムを同定することが可能である。

　次に、動作について説明する。同定対象とする動的システムは、１入力Ｐ出力のｎ次元線形離散時間システムとして上記の式（１）で記述できるとする。この動的システムの（複素）周波数応答Ｈ（ｋΔｆ）（ｋ＝０，１，２，…，Ｎ－１）が与えられた場合、図９に示す逆フーリエ変換部６１によって下記の式（２４）で逆フーリエ変換し、動的システムのインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）を得る。

　ただし、サンプリング周期Ｔ_ｓ＝Ｔ／Ｎであり、サンプリング周波数ｆｓ＝１／Ｔ_ｓ＝Ｎ／Ｔであり、周波数分解能Δｆ＝１／Ｔであり、時刻ｔ＝ｊＴ_ｓ＝ｊＴ／Ｎであり、周波数ｆ＝ｋΔｆ＝ｋ／Ｔである。

　本実施の形態では、動的システムの周波数応答Ｈ（ｋΔｆ）（ｋ＝０，１，２，…，Ｎ－１）に基づいて、上記の式（２４）で得られた動的システムのインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）を、システム同定装置１０への入力とする。本実施の形態においては、図１に示すように、システム同定装置１０では、上記の式（２４）で得られた動的システムのインパルス応答ｇ（ｊＴ_ｓ）（ｊ＝０，１，２，…）を入力とし、直達項同定部１において動的システムを記述する上記の式（１）の線形離散時間システムの直達項Ｄ_ｄを式（４）で同定すると共に、ブロックハンケル行列生成部２によって下記の式（２５）で与えられるブロックハンケル行列Ｈ_ｋｌを生成する。

　次に、特異値分解部３によって、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｈ_ｋｌを特異値分解し、上記の式（６）で与えられるブロックハンケル行列Ｈ_ｋｌの第１の直交行列Ｕ、第２の直交行列Ｖ及び特異値σ_ｉ（ｉ＝１，２，３，…）を出力する。

　システム安定性評価部３４は、作業者が指定したシステム次元の探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）に属する次元ｎ_ｉのそれぞれに対して、再帰的システム行列推定部３１で同定したシステム行列Ａ_ｄ，ｎｉに基づいて、線形離散時間システムの安定性を上記の式（２２）で評価する。

　次に、システム特性推定部３２において、システム安定性評価部３４で安定システムと判断された次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａ_ｄ，ｎｉ，Ｂ_ｄ，ｎｉ，Ｃ_ｄ，ｎｉ及び直達項同定部１から出力される直達項Ｄ_ｄに基づいて、線形離散時間システムの推定周波数応答Ｈ＾_ｎｉ（ｋΔｆ）（ｋ＝０，１，２，…，Ｎ－１）を算出する。

　なお、本明細書中で、「Ｈ＾_ｎｉ」は「Ｈ_ｎｉ」の上に「＾（サーカムフレックス）」が配置された文字の代替表記である。

　本実施の形態の周波数応答によるシステム同定装置６０では、「周波数領域において実際の周波数応答に最も適合する」という前提の下、システム次元推定部３３において最適なシステム次元ｎを決定する。具体的には、動的システムの実際の周波数応答Ｈ（ｋΔｆ）（ｋ＝０，１，２，…，Ｎ－１）に基づいて、例えば高ゲインかつ低周波領域に重みを持たせた下記の式（２６）のような重み関数Ｗ（ｋΔｆ）（ｋ＝０，１，２，…，Ｎ－１）を決定し、システム特性推定部３２から出力される線形離散時間システムの推定周波数応答Ｈ＾_ｎｉ（ｋΔｆ）と動的システムの実際の周波数応答Ｈ（ｋΔｆ）（図８における動的システムのシステム特性）との周波数領域における誤差２乗値に対して、重み関数Ｗ（ｋΔｆ）を乗じた値の加算値ｅ_ｎｉ（ｎ_ｉ：安定システムとなる次元）を下記の式（２７）で算出する。

　この重み付き誤差２乗和ノルム||ｅ_ｎｉ||を最小とする次元ｎ_ｉが、「重み関数に応じて、周波数領域において実際の周波数応答に最も適合する」安定なシステム次元ｎとなる。ここでは、図５に示すように誤差２乗和ノルム分布４１が、上記の式（１３）で与えられる誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ_６）。

　最後に、システム行列同定部５において、特異値分解部３から出力される第１の直交行列Ｕ、第２の直交行列Ｖ、及び特異値σ_ｉ（ｉ＝１，２，３，…）と、システム次元決定部４ａから出力されるシステム次元ｎに基づき、線形離散時間システムの直達項Ｄ_ｄを除くシステム行列Ａ_ｄ，Ｂ_ｄ，Ｃ_ｄを、上記の式（１４）で同定して、直達項同定部１で同定した直達項Ｄ_ｄ及びシステム行列同定部５で同定したシステム行列Ａ_ｄ，Ｂ_ｄ，Ｃ_ｄを、動的システムを記述する線形離散時間システムとして出力する。

　このように、本実施の形態の周波数応答によるシステム同定装置６０により、現実の周波数応答をインパルス応答に変換し、当該インパルス応答から算出したブロックハンケル行列Ｈ_ｋｌの特異値σ_ｉ（ｉ＝１，２，３，…）がなだらかな単調減少となり、有意な値を持つ特異値と同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して周波数領域で重み関数に応じて一致度の高いシステム次元ｎの決定及び動的システムを記述する線形離散時間システムの同定が可能となる。

　更に逆フーリエ変換に基づいて、動的システムの周波数応答をインパルス応答に変換することで、動的システムの現実の周波数応答から、当該システムを記述する線形離散時間システムを正確に同定することが可能となる。

　なお、本実施の形態の周波数応答によるシステム同定装置６０では、線形離散時間システムのシステム特性を推定周波数応答として算出し、当該周波数応答と動的システムの実際の周波数応答との周波数領域における重み付き誤差２乗和ノルム分布４１が誤差２乗和ノルムしきい値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定している。しかしながら、本発明はこれに限定されるものではなく、線形離散時間システムのシステム特性を推定インパルス応答として算出し、当該インパルス応答と動的システムの周波数応答を逆フーリエ変換して得られる実際のインパルス応答との時間領域における誤差２乗和に基づいて、システム次元ｎを決定してもよい。

　システム次元決定部４ａが、探索範囲ｎ_ｉ＝（ｎ_１，ｎ_２，…，ｎ_ａ）（ただしｎ_１＜ｎ_２＜…＜ｎ_ａ）に属する各次元に対して、線形離散時間システムのシステム特性を推定周波数応答として算出して出力するシステム特性推定部３２と、システム特性推定部３２から出力される線形離散時間システムの推定周波数応答とインパルス応答をフーリエ変換した実際の周波数応答との周波数領域における誤差２乗和ノルムがしきい値４２以下となる次元のうち最小の次元をシステム次元として決定して出力するシステム次元推定部３３と、を備えるシステム同定装置は、現実のインパルス応答から算出したブロックハンケル行列の特異値がなだらかな単調減少となり、有意な値を持つ特異値と無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、現実の動的システムに対して周波数領域で一致度の高いシステム次元の決定と、動的システムを記述する線形離散時間システムの同定が可能となる。

　システム次元推定部３３が、フーリエ変換した実際の周波数応答に基づいて重み関数を決定し、システム特性推定部３２から出力される線形離散時間システムの推定周波数応答と動的システムの実際の周波数応答との周波数領域における誤差２乗値に対して、重み関数を乗じた値の加算値を算出して、加算値のノルムがしきい値以下となる次元のうち最小の次元をシステム次元として決定して出力するシステム同定装置は、現実のインパルス応答から算出したブロックハンケル行列の特異値がなだらかな単調減少となり、有意な値を持つ特異値と無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、現実の動的システムに対して周波数領域で重み関数に応じて一致度の高いシステム次元の決定と、動的システムを記述する線形離散時間システムの同定が可能となる。

実施の形態４．
　本実施の形態では、同定対象とする動的システムがＤＣサーボモータである場合について説明する。

　図１０は、本発明にかかるシステム同定装置の実施の形態４の全体構成を示すブロック線図である。本実施の形態のシステム同定装置７０は、ＤＣサーボモータ７１と、周波数応答変換部７２と、周波数応答によるシステム同定装置６０と、制御系またはフィルタパラメータ設計部７３とを備える。周波数応答によるシステム同定装置６０は、実施の形態３において図１及び図９を用いて説明したものと同一である。

　本実施の形態では、ＤＣサーボモータ７１の入力電流［Ａ_ｒｍｓ］として、例えばサインスイープを入力し、これをシステム入力とする。また、角速度［ｒａｄ／ｓ］をシステム出力として取得する。システム入出力は、周波数応答変換部７２により、周波数応答Ｈ（ｋΔｆ）に変換される。周波数応答変換部７２は、例えばＦＦＴ（Ｆａｓｔ　Ｆｏｕｒｉｅｒ　Ｔｒａｎｓｆｏｒｍ）などにより下記の式(２８)に基づいて周波数応答を出力する。

　周波数応答によるシステム同定装置６０は、この動的システムの周波数応答及びシステム次元の探索範囲を入力とし、ＤＣサーボモータ７１の入力電流［Ａ_ｒｍｓ］から角速度［ｒａｄ／ｓ］までの動的システムを記述する線形離散時間システムを同定して制御系またはフィルタパラメータ設計部７３に出力する。このとき、システム次元の探索範囲は、ｎ_i＝（１，２，…，５０）のように、予測されるシステム次元に対して十分な幅をとって設定すればよい。

　周波数応答によるシステム同定装置６０では、現実の動的システムに対して一致度の高いシステム次元の決定と、動的システムを記述する線形離散時間システムの同定とを可能とするため、当該線形離散時間システムをサーボモータの制御系またはフィルタのパラメータ設計に利用することができる。

　以上のように、本発明にかかるシステム同定装置は、インパルス応答によるシステム同定装置として有用である。

　１　直達項同定部、２　ブロックハンケル行列生成部、３　特異値分解部、４，４ａ　システム次元決定部、５　システム行列同定部、１０，７０　システム同定装置、１１　システム入力、１２　システム出力、１３　システム入力しきい値、２１　（理想的なインパルス応答における）特異値分布、２２　（現実のインパルス応答における）特異値分布、３１　再帰的システム行列推定部、３２　システム特性推定部、３３　システム次元推定部、３４　システム安定性評価部、４１　誤差２乗和ノルム分布、４２　誤差２乗和ノルムしきい値、５０　打撃加振によるシステム同定装置、５１　システム入力印加時刻特定部、５２　インパルス応答変換部、６０　周波数応答によるシステム同定装置、６１　逆フーリエ変換部、７１　ＤＣサーボモータ、７２　周波数応答変換部、７３　制御系またはフィルタパラメータ設計部。

Claims

　同定対象とする動的システムのインパルス応答と、指定されたシステム次元の探索範囲と、を入力とするシステム同定装置であって、
　前記インパルス応答から前記動的システムを記述する線形離散時間システムの直達項を同定して出力する直達項同定部と、
　前記インパルス応答からブロックハンケル行列を生成して出力するブロックハンケル行列生成部と、
　前記ブロックハンケル行列生成部から出力される前記ブロックハンケル行列の特異値分解により、前記ブロックハンケル行列の左特異ベクトルを列ベクトルとする第１の直交行列、前記ブロックハンケル行列の右特異ベクトルを列ベクトルとする第２の直交行列及び前記ブロックハンケル行列の特異値を出力する特異値分解部と、
　前記第１の直交行列、前記第２の直交行列、前記特異値及び前記探索範囲に基づいて、当該探索範囲に属する各次元に対する前記線形離散時間システムのシステム行列のうち前記直達項以外のシステム行列を同定して、当該システム行列及び前記直達項に基づいて算出した前記線形離散時間システムのシステム特性と前記動的システムの実際のシステム特性との比較から、システム次元を決定して出力するシステム次元決定部と、
　前記第１の直交行列、前記第２の直交行列、前記特異値及び前記システム次元決定部から出力される前記システム次元に基づいて、前記線形離散時間システムのシステム行列のうち前記直達項以外のシステム行列を同定するシステム行列同定部と、を備え、
　前記直達項同定部で同定した前記直達項及び前記システム行列同定部で同定した前記システム行列を前記線形離散時間システムとして出力することを特徴とするシステム同定装置。
　前記システム次元決定部が、
　前記探索範囲に属する各次元に対して、前記線形離散時間システムのシステム特性を推定インパルス応答として算出して出力するシステム特性推定部と、
　前記システム特性推定部から出力される前記線形離散時間システムの前記推定インパルス応答と前記動的システムの実際のインパルス応答との時間領域における誤差２乗和ノルムがしきい値以下となる次元のうち最小の次元をシステム次元として決定して出力するシステム次元推定部と、を備えることを特徴とする請求項１に記載のシステム同定装置。
　前記システム次元決定部が、
　前記探索範囲に属する各次元に対して、前記線形離散時間システムのシステム特性を推定周波数応答として算出して出力するシステム特性推定部と、
　前記システム特性推定部から出力される前記線形離散時間システムの推定周波数応答と前記インパルス応答をフーリエ変換した実際の周波数応答との周波数領域における誤差２乗和ノルムがしきい値以下となる次元のうち最小の次元をシステム次元として決定して出力するシステム次元推定部と、を備えることを特徴とする請求項１に記載のシステム同定装置。
　前記システム次元推定部が、
　前記フーリエ変換した前記実際の周波数応答に基づいて重み関数を決定し、
　前記システム特性推定部から出力される前記線形離散時間システムの推定周波数応答と前記動的システムの前記実際の周波数応答との周波数領域における誤差２乗値に対して、当該重み関数を乗じた値の加算値を算出して、当該加算値のノルムがしきい値以下となる次元のうち最小の次元をシステム次元として決定して出力することを特徴とする請求項３に記載のシステム同定装置。
　前記しきい値は、誤差２乗和許容値と誤差２乗和ノルムの最小値の積により規定されることを特徴とする請求項２から４のいずれか一項に記載のシステム同定装置。
　前記システム次元決定部が、
　前記探索範囲の中で第１の次元よりも１段階低い第２の次元に対応したシステム行列の同定結果と、前記第１の直交行列、前記第２の直交行列及び前記特異値からそれぞれ前記第２の次元よりも大きく前記第１の次元以下となる左特異ベクトル、右特異ベクトル及び特異値と、を用いて、再帰的手法により前記第１の次元に対応するシステム行列を同定する再帰的システム行列推定部を備えることを特徴とする請求項１に記載のシステム同定装置。
　前記インパルス応答に対応する情報が前記動的システムを打撃加振したときのシステム入力及びシステム出力であり、
　前記システム入力が有意な値を持つ複数の時刻を特定するシステム入力印加時刻特定部と、
　前記システム入力印加時刻特定部から出力される複数の時刻に対応するシステム入力の加算値をインパルス入力振幅とし、前記システム入力印加時刻特定部から出力される複数の時刻の最大値をインパルス入力印加時刻として、当該インパルス入力印加時刻以降のシステム出力を前記インパルス入力振幅で除した信号を前記インパルス応答として出力するインパルス応答変換部と、を備え、
　前記同定対象とする前記動的システムを打撃加振した場合のシステム入力及びシステム出力並びに前記探索範囲を入力として、前記線形離散時間システムを同定することを特徴とする請求項１に記載のシステム同定装置。
　前記システム入力印加時刻特定部が、
　比率しきい値とシステム入力の最大値とを乗じた値をシステム入力しきい値とし、システム入力がシステム入力しきい値以上となる時刻を、システム入力が有意な値を持つ複数の時刻として特定することを特徴とする請求項７に記載のシステム同定装置。
　前記インパルス応答に対応する情報が前記動的システムの周波数応答であり、
　前記周波数応答の逆フーリエ変換により対応するインパルス応答を出力する逆フーリエ変換部を備え、
　前記同定対象とする前記動的システムの周波数応答及び前記探索範囲を入力として、前記線形離散時間システムを同定することを特徴とする請求項１に記載のシステム同定装置。
　前記システム次元決定部が、
　前記探索範囲に属する各次元に対して前記線形離散時間システムの安定性を評価するシステム安定性評価部を備え、
　安定システムとなる次元に対応する前記線形離散時間システムのシステム特性からシステム次元を決定することを特徴とする請求項１に記載のシステム同定装置。