(57)【特許請求の範囲】
【請求項１】
同定対象とする動的システムに対して擬似ランダム入力を印加した場合のシステム入出力を入力とするシステム同定装置において、
前記動的システムのシステム入出力から同定に適用する同定用入出力データを抽出するシステム入出力抽出部と、
前記同定用入出力データに基づいて、ブロックハンケル行列を生成するブロックハンケル行列生成部と、
前記ブロックハンケル行列に基づいて、前記動的システムの入力ベクトルおよび出力ベクトルを生成する入出力ベクトル生成部と、
前記ブロックハンケル行列を結合してデータ行列を生成し、当該データ行列をＬＱ分解した部分行列を出力するＬＱ分解部と、
前記部分行列と前記ブロックハンケル行列に基づいて、平行射影を生成する平行射影生成部と、
前記平行射影の特異値分解により、前記平行射影の左特異ベクトルを列ベクトルとする第１の直交行列、当該平行射影の右特異ベクトルを列ベクトルとする第２の直交行列および当該平行射影の特異値を出力する特異値分解部と、
前記第２の直交行列および前記特異値と、前記動的システムの入力ベクトルおよび出力ベクトルと、指定されたシステム次元の探索範囲に基づいて、当該探索範囲に属するそれぞれの次元に対して動的システムを記述する線形離散時間システムのシステム行列を同定し、さらに当該システム行列に基づいて算出した線形離散時間システムのシステム特性と、動的システムの実際のシステム特性との比較から、システム次元を決定するシステム次元決定部と、
前記第２の直交行列および特異値と、前記決定されたシステム次元に基づいて、前記動的システムの状態ベクトルを生成する状態ベクトル生成部と、
前記動的システムの入力ベクトルおよび出力ベクトルならびに前記動的システムの状態ベクトルに基づいて、当該動的システムを記述する線形離散時間システムのシステム行列を同定するシステム行列同定部と、を備え、
前記同定されたシステム行列を、前記動的システムを記述する線形離散時間システムとして出力することを特徴とするシステム同定装置。
【請求項２】
前記システム次元決定部は、
前記探索範囲に属するそれぞれの次元に対して、
同定した線形離散時間システムに関して実際の同定用入力データを適用した場合のシステム出力を算出し、当該システム出力を線形離散時間システムのシステム特性として出力するシステム特性推定部と、
前記線形離散時間システムのシステム出力と、前記動的システムの実際の同定用出力データとの時間領域における誤差２乗和のノルムが設定閾値以下となる次元のうち、最小の次元をシステム次元として決定して出力するシステム次元推定部と
を備えたことを特徴とする請求項１に記載のシステム同定装置。
【請求項３】
前記システム次元決定部が、
前記探索範囲に属するそれぞれの次元に対して、
同定した線形離散時間システムの周波数応答を算出し、当該周波数応答を線形離散時間システムのシステム特性として出力するシステム特性推定部と、
前記線形離散時間システムの周波数応答と、前記動的システムのシステム入出力から得られる実際の周波数応答との周波数領域における誤差２乗和のノルムが設定閾値以下となる次元のうち、最小の次元をシステム次元として決定して出力するシステム次元推定部と、
を備えたことを特徴とする請求項１に記載のシステム同定装置。
【請求項４】
前記システム次元推定部は、
前記動的システムのシステム入出力から得られる実際の周波数応答に基づいて重み関数を決定し、
前記線形離散時間システムの周波数応答と、前記動的システムの実際の周波数応答との周波数領域における誤差２乗値に対して、当該重み関数を乗じた値の加算値を算出し、
当該加算値のノルムが設定閾値以下となる次元のうち、最小の次元をシステム次元として決定して出力する
ことを特徴とする請求項３に記載のシステム同定装置。
【請求項５】
前記システム次元決定部は、
前記探索範囲に属する第１の次元に対応するシステム行列の同定に関して、
当該探索範囲の中で第１の次元よりも１段階低い第２の次元に対応したシステム行列の同定結果と、
前記第２の直交行列および特異値のうち、それぞれ前記第２の次元よりも大きく、前記第１の次元以下となる右特異ベクトルおよび特異値、ならびに、前記動的システムの入力ベクトルおよび出力ベクトルを用いて、再帰的手法により前記第１の次元に対応するシステム行列を同定する再帰的システム行列推定部と、
を備えたことを特徴とする請求項１に記載のシステム同定装置。
【請求項６】
前記システム入出力抽出部は、
設定された比率閾値とシステム入力の最大値とを乗じた値をシステム入力閾値とし、システム入力の絶対値がシステム入力閾値以上となる時刻の最小値を擬似ランダム入力印加時刻として、擬似ランダム入力印加時刻以降のシステム入力およびシステム出力を、それぞれ同定用入力データおよび同定用出力データとして抽出することを特徴とする請求項１に記載のシステム同定装置。
【請求項７】
前記システム次元決定部は、
前記探索範囲に属するそれぞれの次元に対して線形離散時間システムの安定性を評価するシステム安定性評価部を備え、
安定システムとなる次元に対応する線形離散時間システムのシステム特性から、システム次元を決定することを特徴とする請求項１に記載のシステム同定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、対象とする動的システムに対して擬似ランダム入力を印加した場合のシステム入出力に基づいて、当該システムの数学モデルを構築するシステム同定装置に関する。
【背景技術】
【０００２】
従来の擬似ランダム入力によるシステム同定装置として、例えば非特許文献１に記載されたＮ４ＳＩＤ法に基づくものがある。このＮ４ＳＩＤ法では、線形離散時間システム（Ａd,Ｂd,Ｃd,Ｄd）で記述される動的システムに擬似ランダム入力を印加した場合のシステム入出力に基づいて、システム入力に関するブロックハンケル行列（Ｕp,Ｕf）、およびシステム出力に関するブロックハンケル行列（Ｙp,Ｙf）を生成し、さらにブロックハンケル行列（Ｕf,Ｙf）に基づいて、入出力ベクトル（~ＵK|K,~ＹK|K）を生成する。なお、「~」の表記は、本来であれば「Ｕ」の文字の上部に横線（オーバーバー）を引くべきところであるが、その表記ができない。このため、本明細書では、イメージで挿入する数式部分を除き、横線（オーバーバー）を「~」で代替する。【０００３】
次に、上記ブロックハンケル行列を結合したデータ行列をＬＱ分解し、ＬＱ分解により得られた部分行列と、ブロックハンケル行列Ｕp,Ｙpから平行射影Θを生成する。この平行射影Θを特異値分解し、有意な値を持つ特異値の個数をシステム次元として決定して、特異値分解の結果と決定したシステム次元から、動的システムの状態ベクトル（~ＸK,~ＸK+1）を算出する。最後に、入出力ベクトル（~ＵK|K,~ＹK|K）および状態ベクトル（~ＸK,~ＸK+1）に対して最小２乗法を適用することで、動的システムを記述する線形離散時間システム（Ａd,Ｂd,Ｃd,Ｄd）を同定している。【０００４】
また、従来の擬似ランダム入力によるシステム同定装置の他の例として、例えば特許文献１に記載された露光装置及び除振装置、システム同定装置及びその方法がある。【０００５】
この露光装置及び除振装置、システム同定装置及びその方法では、対象とする動的システムに擬似ランダム入力を印加した場合のシステム入出力に基づいて、Ｎ４ＳＩＤ法に代表される部分空間法により、動的システムの状態方程式を同定する。このとき、同定する状態方程式のシステム次元を、動的システムの運動方程式から定まるシステム次元に一致させることで、運動方程式に基づく特性方程式と、同定した状態方程式に基づく特性方程式との比較により、運動方程式に含まれる未知の物理パラメータを同定している。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２０００−８２６６２号公報
【非特許文献】
【０００７】
【非特許文献１】システム同定−部分空間法からのアプローチ、朝倉書店  ｐｐ．１１７−１２０
【発明の概要】
【発明が解決しようとする課題】
【０００８】
このような擬似ランダム入力によるシステム同定装置では、有意な値を持つ特異値の個数、または動的システムの運動方程式から定まるシステム次元から、対象とする動的システムのシステム次元を決定している。【０００９】
しかしながら、現実のシステム入出力から算出した平行射影Θの特異値は、なだらかな単調減少となる場合が多数存在し、この場合は有意な値を持つ特異値と、無視できる微小な値となる特異値との境界が不明確となる。したがって、非特許文献１に記載された従来のシステム同定装置では、システム次元の決定が作業者の判断に依存し、常に最適なシステム次元を決定しているとは限らないか、もしくはシステム次元の決定に関して試行錯誤が必要になるという問題があった。【００１０】
また、動的システムのモデリングによって得られる運動方程式では、動的システムの実際の動特性を全て記述することは困難であり、一般に、“運動方程式から定まるシステム次元＜動的システムの実際のシステム次元”となる。したがって、特許文献１に記載された従来のシステム同定装置では、そもそも動的システムを記述する最適なシステム次元を決定することができないという問題があった。【００１１】
加えて、従来の擬似ランダム入力によるシステム同定装置では、同定結果として得られる線形離散時間システム（Ａd,Ｂd,Ｃd,Ｄd）の安定性は全く考慮されていないため、現実の動的システムが安定であるにも関わらず、不安定システムとして同定される場合があるという問題もあった。【００１２】
本発明は、上記に鑑みてなされたものであって、現実のシステム入出力から算出した平行射影Θの特異値がなだらかな単調減少となり、したがって有意な値を持つ特異値と、無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、最適なシステム次元を決定することができるシステム同定装置を得ることを目的とする。【００１３】
また、本発明は、現実の動的システムが安定であることが明確である場合に、安定システムに限定して同定することができるシステム同定装置を得ることを目的とする。
【課題を解決するための手段】
【００１４】
上述した課題を解決し、目的を達成するため、本発明に係るシステム同定装置は、同定対象とする動的システムに対して擬似ランダム入力を印加した場合のシステム入出力を入力とするシステム同定装置において、前記動的システムのシステム入出力から同定に適用する同定用入出力データを抽出するシステム入出力抽出部と、前記同定用入出力データに基づいて、ブロックハンケル行列を生成するブロックハンケル行列生成部と、前記ブロックハンケル行列に基づいて、前記動的システムの入力ベクトルおよび出力ベクトルを生成する入出力ベクトル生成部と、前記ブロックハンケル行列を結合してデータ行列を生成し、当該データ行列をＬＱ分解した部分行列を出力するＬＱ分解部と、前記部分行列と前記ブロックハンケル行列に基づいて、平行射影を生成する平行射影生成部と、前記平行射影の特異値分解により、前記平行射影の左特異ベクトルを列ベクトルとする第１の直交行列、当該平行射影の右特異ベクトルを列ベクトルとする第２の直交行列および当該平行射影の特異値を出力する特異値分解部と、前記第２の直交行列および前記特異値と、前記動的システムの入力ベクトルおよび出力ベクトルと、指定されたシステム次元の探索範囲に基づいて、当該探索範囲に属するそれぞれの次元に対して動的システムを記述する線形離散時間システムのシステム行列を同定し、さらに当該システム行列に基づいて算出した線形離散時間システムのシステム特性と、動的システムの実際のシステム特性との比較から、システム次元を決定するシステム次元決定部と、前記第２の直交行列および特異値と、前記決定されたシステム次元に基づいて、前記動的システムの状態ベクトルを生成する状態ベクトル生成部と、前記動的システムの入力ベクトルおよび出力ベクトルならびに前記動的システムの状態ベクトルに基づいて、当該動的システムを記述する線形離散時間システムのシステム行列を同定するシステム行列同定部と、を備え、前記同定されたシステム行列を、前記動的システムを記述する線形離散時間システムとして出力することを特徴とする。
【発明の効果】
【００１５】
この発明によれば、同定対象とする動的システムにおいて、現実のシステム入出力から算出した平行射影の特異値がなだらかな単調減少となり、したがって有意な値を持つ特異値と、無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元の決定から試行錯誤を排除し、常に最適なシステム次元の決定と、動的システムを記述する線形離散時間システムの同定が可能となる。
【図面の簡単な説明】
【００１６】
【図１】図１は、実施の形態１および実施の形態２に係るシステム同定装置の全体構成を示すブロック線図である。
【図２】図２は、実施の形態１のシステム同定装置におけるシステム入出力の時間波形を示す概略図である。
【図３】図３は、実施の形態１および実施の形態２に係るシステム同定装置における平行射影の特異値と次元との関係を示す概略図である。
【図４】図４は、実施の形態１のシステム同定装置におけるシステム次元決定部の内部構成を示すブロック線図である。
【図５】図５は、実施の形態１および実施の形態２に係るシステム同定装置において、同定した線形離散時間システムの時間領域または周波数領域における誤差２乗和のノルムと次元との関係を示す概略図である。
【図６】図６は、実施の形態２のシステム同定装置における動的システムをＭ系列加振した場合のシステム入出力の時間波形を示す概略図である。
【図７】図７は、実施の形態２のシステム同定装置におけるシステム次元決定部の内部構成を示すブロック線図である。
【図８】図８は、実施の形態３に係る全体構成を示すブロック線図である。
【発明を実施するための形態】
【００１７】
以下、添付図面を参照し、本発明の実施の形態に係るシステム同定装置について説明する。なお、以下に示す実施の形態により本発明が限定されるものではない。
【００１８】
実施の形態１．
図１は、実施の形態１に係るシステム同定装置の全体構成を示すブロック線図であり、図２は、実施の形態１のシステム同定装置におけるシステム入出力の時間波形を示す概略図である。
【００１９】
実施の形態１に係るシステム同定装置１０では、図１および図２に示すように、同定対象とする動的システムに対して擬似ランダム入力を印加した場合のシステム入力１１（ｕ(jＴS)(ｊ=0,1,2,…)）および、システム出力１２（ｙ(jＴS)(ｊ=0,1,2,…)）を入力とする。
【００２０】
システム入出力抽出部１は、予め設定された比率閾値とシステム入力１１の最大値とを乗じた値で決定されるシステム入力閾値１３に対して、システム入力１１の絶対値がシステム入力閾値１３以上となる時刻の最小値を擬似ランダム入力印加時刻（図２では３ＴS)として、擬似ランダム入力印加時刻以降のシステム入力１１およびシステム出力１２を、それぞれ同定用入力データ（ｕid(jＴS)(ｊ=0,1,2,…)）および、同定用出力データ（ｙid(jＴS)(ｊ=0,1,2,…)）として抽出して出力する。
【００２１】
ブロックハンケル行列生成部２は、システム入出力抽出部１から出力される同定用入力データｕid(jＴS)(ｊ=0,1,2,…)および、同定用出力データｙid(jＴS)(ｊ=0,1,2,…）に基づいて、ブロックハンケル行列Ｕp,ＵfおよびＹp,Ｙfを生成する。
【００２２】
入出力ベクトル生成部３は、ブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfに基づいて、動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kを生成する。
【００２３】
ＬＱ分解部４は、ブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfを結合したデータ行列を生成し、当該データ行列をＬＱ分解した部分行列Ｌ22,Ｌ32を生成して出力する。
【００２４】
平行射影生成部５は、ＬＱ分解部４から出力される部分行列Ｌ22,Ｌ32と、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｕp,Ｙpに基づいて、動的システムの平行射影Θを生成する。
【００２５】
特異値分解部６は、平行射影生成部５から出力される平行射影Θを特異値分解し、平行射影Θの左特異ベクトルを列ベクトルとする第１の直交行列Ｕ、平行射影Θの右特異ベクトルを列ベクトルとする第２の直交行列Ｖおよび、平行射影Θの特異値σi(i=1,2,3, …)を出力する。
【００２６】
システム次元決定部７は、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3, …)と、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kならびに、作業者が指定したシステム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に基づいて、当該探索範囲に属する次元ｎi(i=1,2,…,a)のそれぞれに対して動的システムを記述する線形離散時間システムのシステム行列を同定する。さらに当該システム行列に基づいて、当該探索範囲に属する次元ｎi(i=1,2,…,a)のそれぞれに対応する線形離散時間システムに対して実際の同定用入力データｕid(jＴS)(ｊ=0,1,2,…)を適用した場合のシステム出力を算出し、動的システムの実際の同定用出力データｙid(jＴS)(ｊ=0,1,2,…）（図１では動的システムのシステム特性と記載）との比較から、システム次元ｎを決定する。
【００２７】
状態ベクトル生成部８では、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3, …)および、システム次元決定部７から出力されるシステム次元ｎsに基づいて、動的システムの状態ベクトル~ＸK+1,~ＸKを生成する。
【００２８】
システム行列同定部９は、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|Kおよび、出力ベクトル~ＹK|Kと、状態ベクトル生成部８から出力される動的システムの状態ベクトル~ＸK+1,~ＸKに基づいて、動的システムを記述する線形離散時間システムのシステム行列Ａd,Ｂd,Ｃd,Ｄdを同定して出力する。
【００２９】
図３は、実施の形態１のシステム同定装置１０における平行射影Θの特異値σiと次元(i=1,2,3, …)の関係を示す概略図であり、図４は、実施の形態１のシステム同定装置１０におけるシステム次元決定部７の内部構成を示すブロック線図であり、図５は、実施の形態１のシステム同定装置１０における同定した線形離散時間システムのシステム出力と、動的システムの実際のシステム出力との時間領域における誤差２乗和のノルム||ｅn||と次元ｎi(i=1,2,…,a)との関係を示す概略図である。
【００３０】
図３に示すように、動的システムのシステム入出力から算出した平行射影Θの特異値σi(i=1,2,3, …)は、理想的には次元(i=1,2,3, …)に対して、例えば特異値分布２１に示す関係となる。この場合、有意な値を持つ特異値の個数を明確に規定することができ、当該個数が動的システムのシステム次元ｎに対応する（図３の場合はシステム次元ｎ＝４）。
【００３１】
一方、観測ノイズ等の影響を受ける現実のシステム入出力に基づいて算出した特異値σiは、次元(i=1,2,3, …)に対して例えば特異値分布２２に示す関係となるため、有意な値を持つ特異値と、無視できる微小な値となる特異値との境界が不明確となり、常に最適なシステム次元ｎを決定しているとは限らない。よって、システム次元ｎの決定に関して試行錯誤が必要になるという問題が発生する。
【００３２】
そこで、実施の形態１に係るシステム同定装置１０では、システム次元決定部７により、図４に示す処理を実行する。具体的には、以下の通りである。
【００３３】
システム次元決定部７には、再帰的システム行列推定部３１、システム特性推定部３２およびシステム次元推定部３３が設けられる。
【００３４】
再帰的システム行列推定部３１は、作業者が予め指定したシステム次元の探索範囲ｎi(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属する第１の次元ｎiに対応するシステム行列の同定に関して、第１の次元ｎiよりも１段階低い第２の次元ｎi-1に対応したシステム行列Ａd,n(i-1),Ｂd,n(i-1),Ｃd,n(i-1),Ｄd,n(i-1)［ここでのｎ（ｉ−１）はｎｉ−１を指す］の同定結果と、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3, …)のうち、それぞれ第２の次元ｎi-1よりも大きく、第１の次元ｎi以下となる右特異ベクトルｖjおよび特異値σj(j=ｎi-1+1,ｎi-1+2,…,ｎi)と、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|K、および出力ベクトル~ＹK|Kとを用いて、再帰的手法によって第１の次元ｎiに対応するシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niを同定する。
【００３５】
次にシステム特性推定部３２は、システム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属するそれぞれの次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niに基づいて、同定した線形離散時間システムに対して、実際の同定用入力データｕid(jＴS)(ｊ=0,1,2,…)を適用した場合のシステム出力を算出する。
【００３６】
なお、再帰的システム行列推定部３１およびシステム特性推定部３２の処理は、ｉをインクリメントし、ｉ＝ａとなるまで実行される。
【００３７】
システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムのシステム出力と、動的システムの実際の同定用出力データｙid(jＴS)(ｊ=0,1,2,…)（図４では動的システムのシステム特性と記載）との時間領域における誤差２乗和ｅni(i=1,2,…,ａ)を算出し、図５に示したように誤差２乗和のノルム||ｅni||の分布４１が、予め設定された誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する構成としている（図５の場合はシステム次元ｎ＝ｎ6）。
【００３８】
次に、実施の形態１に係るシステム同定装置の動作について説明する。
【００３９】
同定対象とする動的システムが、次式のように、１入力Ｐ出力のｎ次元線形離散時間システムで記述できるとする。
【００４０】
【数１】
0006009105.tif 000002
【００４１】
上記の動的システムへのシステム入力ｕ(jＴS)を擬似ランダム入力で構成すると、当該システム入力ｕ(jＴS)および、上記［数１］式に対応するシステム出力ｙ(jＴS)は、例えば図２に示したシステム入力１１、およびシステム出力１２のような時間波形となる。
【００４２】
ここで、図１および図２の説明で述べたように、予め設定された比率閾値とシステム入力１１（ｕ(jＴS)）の最大値とを乗じた次式をシステム入力閾値１３として使用する。
【００４３】
【数２】
0006009105.tif 000003
【００４４】
システム入出力抽出部１は、システム入力１１の絶対値がシステム入力閾値１３以上となる時刻の最小値を擬似ランダム入力印加時刻ｊminＴSとして特定する（図２の場合はｊminＴS＝３ＴS）。
【００４５】
また、システム入出力抽出部１は、擬似ランダム入力印加時刻ｊminＴS以降のシステム入力１１およびシステム出力１２を次式を用いて抽出する。
【００４６】
【数３】
0006009105.tif 000004
【００４７】
さらに、システム入出力抽出部１は、上記［数３］を用いて抽出したそれぞれの値を同定用入力データｕid(ｊＴS)および同定用出力データｙid(ｊＴS)とすることで、対象とする動的システムのシステム入出力から擬似ランダム入力印加前のシステム静止時間領域データを除去する。
【００４８】
ブロックハンケル行列生成部２は、システム入出力抽出部１から出力される同定用入力データｕid(jＴS)(ｊ=0,1,2,…)および同定用出力データｙid(jＴS)(ｊ=0,1,2,…)に基づいて、次式で与えられるブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfを生成する。
【００４９】
【数４】
0006009105.tif 000005
【００５０】
入出力ベクトル生成部３は、このブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfに基づいて、次式で与えられる動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kを生成する。
【００５１】
【数５】
0006009105.tif 000006
【００５２】
ＬＱ分解部４は、ブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfを結合した次式で与えられるデータ行列を生成する。
【００５３】
【数６】
0006009105.tif 000007
【００５４】
また、ＬＱ分解部４は、上記データ行列を次式のようにＬＱ分解すると共に、ＬＱ分解した行列の要素から、部分行列Ｌ22,Ｌ32を出力する。
【００５５】
【数７】
0006009105.tif 000008
【００５６】
平行射影生成部５は、ＬＱ分解部４から出力される部分行列Ｌ22,Ｌ32と、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｕp,Ｙpに基づいて、次式で定義される動的システムの平行射影Θを生成する。
【００５７】
【数８】
0006009105.tif 000009
【００５８】
特異値分解部６は、上式で示される平行射影Θを特異値分解することで、次式で与えられる平行射影Θの左特異ベクトルｕjを列ベクトルとする第１の直交行列Ｕ、平行射影Θの右特異ベクトルｖjを列ベクトルとする第２の直交行列Ｖおよび、平行射影Θの特異値σi(i=1,2,3, …)を出力する。
【００５９】
【数９】
0006009105.tif 000010
【００６０】
対象とする動的システムのシステム次元ｎは、平行射影Θの特異値において有意な値を持つものはｎ個であり、それらに比してｎ＋１個目以降は十分小さいとする以下の関係に基づいて決定することができる。
【００６１】
【数１０】
0006009105.tif 000011
【００６２】
図３に示したように、動的システムのシステム入出力から算出した平行射影Θの特異値σiは、理想的には次元(i=1,2,3, …)に対して、例えば特異値分布２１に示す関係となる。この場合、有意な値を持つ特異値の個数を明確に規定することができ、当該個数から動的システムのシステム次元ｎを決定することができる（図３の場合はシステム次元ｎ＝４）。一方、観測ノイズ等の影響を受ける現実のシステム入出力に基づいて算出した特異値σiは、次元(i=1,2,3, …)に対して、例えば特異値分布２２に示す関係となるため、有意な値を持つ特異値と、無視できる微小な値となる特異値との境界σn>>σn+1が不明確となる。したがって、従来の手法では、常に最適なシステム次元ｎを決定しているとは限らず、最適なシステム次元ｎを決定するためには、試行錯誤が必要になるという課題があった。
【００６３】
そこで、実施の形態１に係るシステム同定装置１０では、“時間領域において実際のシステム入出力に最も適合する”という前提の下、システム次元決定部７において最適なシステム次元ｎを決定する。システム次元決定部７では、図１に示したように特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3, …)と、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kならびに、作業者が指定したシステム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）とに基づいて、当該探索範囲に属する次元ｎi(i=1,2,…,a)のそれぞれに対して動的システムを記述する線形離散時間システムのシステム行列を同定する。さらに当該システム行列に基づいて、当該探索範囲に属する次元ｎi(i=1,2,…,a)のそれぞれに対応する線形離散時間システムに対して、実際の同定用入力データｕid(jＴS)(ｊ=0,1,2,…)を適用した場合のシステム出力を算出し、動的システムの実際の同定用出力データｙid(jＴS)(ｊ=0,1,2,…）（図１では動的システムのシステム特性と記載）との比較から、システム次元ｎを決定する。
【００６４】
具体的には図４に示したように、再帰的システム行列推定部３１は、作業者が指定したシステム次元の探索範囲ｎi(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属する第１の次元ｎiに対応するシステム行列の同定に関して、第１の次元ｎiよりも１段階低い第２の次元ｎi-1に対応したシステム行列Ａd,n(i-1),Ｂd,n(i-1),Ｃd,n(i-1),Ｄd,n(i-1)［ここでのｎ（ｉ−１）はｎｉ−１を指す］の同定結果と、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3, …)のうち、それぞれ第２の次元ｎi-1よりも大きく、第１の次元ｎi以下となる右特異ベクトルｖjおよび特異値σj(j=ｎi-1+1,  ｎi-1+2,…,ｎi)と、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kとを用いて、次式に示す再帰的手法によって第１の次元ｎiに対応するシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niを同定する。
【００６５】
【数１１】
0006009105.tif 000012
【００６６】
システム特性推定部３２は、システム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属する次元ｎiのそれぞれに対して、再帰的システム行列推定部３１から出力されるシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niに基づいて、同定した線形離散時間システムに対して実際の同定用入力データｕid(jＴS)(ｊ=0,1,2,…)（［数３］参照）を適用した場合のシステム出力^ｙid,ni(jＴS)(ｊ=0,1,2,…）を算出する。なお、「^ｙ」の表記は、「ｙ」という文字の上部に「^」の記号を付したことを意味する代替表記である。
【００６７】
また、システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムのシステム出力^ｙid,ni(jＴS)(ｊ=0,1,2,…）と、動的システムの実際の同定用出力データｙid(jＴS)(ｊ=0,1,2,…）（図４では動的システムのシステム特性と記載）との時間領域における誤差２乗和ｅni(i=1,2,…,a）を次式を用いて算出する。
【００６８】
【数１２】
0006009105.tif 000013
【００６９】
上式で示される誤差２乗和のノルム||ｅni||を最小とする次元ｎiが、“時間領域において実際のシステム入出力に最も適合する”システム次元ｎとなる。一方、実際のノルム||ｅni||は、観測ノイズが白色性であれば、そのノイズレベルによらず次元ｎiの増加に伴って単調減少し、図５に示すようにある次元以上でほぼ一定値となる。そこでここでは、システム次元ｎの推定値が必要以上に高次元となることを避けるため、次式で与えられる誤差２乗和ノルム閾値４２を規定する。
【００７０】
【数１３】
0006009105.tif 000014
【００７１】
システム次元推定部３３は、誤差２乗和のノルム||ｅni||の分布４１が上記誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の例では、システム次元ｎ＝ｎ6）。
【００７２】
状態ベクトル生成部８は、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3, …)と、システム次元決定部７から出力されるシステム次元ｎに基づいて、動的システムの状態ベクトル~ＸK,~ＸK+1を次式に従って生成する。
【００７３】
【数１４】
0006009105.tif 000015
【００７４】
最後に、システム行列同定部９は、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kならびに、状態ベクトル生成部８から出力される動的システムの状態ベクトル~ＸK,~ＸK+1に基づいて、動的システムを記述する線形離散時間システムのシステム行列Ａd,Ｂd,Ｃd,Ｄdを次式を用いて同定して出力する。
【００７５】
【数１５】
0006009105.tif 000016
【００７６】
このように、実施の形態１に係るシステム同定装置１０によれば、現実のシステム入出力から算出した平行射影Θの特異値σi(i=1,2,3, …)がなだらかな単調減少となり、したがって有意な値を持つ特異値と、同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して時間領域で一致度の高いシステム次元ｎを決定することができ、動的システムを記述する線形離散時間システムの同定が可能となる。
【００７７】
また、動的システムの現実のシステム入出力から、擬似ランダム入力印加前のシステム静止時間領域データを除去することで、同定精度を向上させることも可能となる。
【００７８】
さらに、再帰的システム行列推定部３１の存在により、現実の動的システムに対して一致度の高いシステム次元ｎを決定するための演算量を低減することが可能となる。
【００７９】
なお、実施の形態１のシステム同定装置１０では、線形離散時間システムに対して実際の同定用入力データを適用した場合のシステム出力をシステム特性として算出し、当該システム出力と、動的システムの実際の同定用出力データとの時間領域における誤差２乗和のノルム分布４１が閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定している。しかしながら、本発明はこれに限定するものではなく、線形離散時間システムのシステム特性を周波数応答として算出し、当該周波数応答と、動的システムの同定用入出力データから得られる実際の周波数応答との周波数領域における誤差２乗和に基づいて、システム次元ｎを決定してもよい。この場合、さらに動的システムの実際の周波数応答に基づいて重み関数を決定し、線形離散時間システムの周波数応答と、動的システムの実際の周波数応答との周波数領域における誤差２乗値に当該重み関数を乗じた値の加算値に基づいて、システム次元ｎを決定してもよい。
【００８０】
実施の形態２．
つぎに、実施の形態２に係るシステム同定装置について説明する。なお、実施の形態２に係るシステム同定装置の全体構成を示すブロック線図、平行射影Θの特異値σiと次元(i=1,2,3, …)の関係を示す概略図および、同定した線形離散時間システムの周波数応答と動的システムの実際の周波数応答との周波数領域における誤差２乗和のノルム||ｅni||と次元ｎi(i=1,2,…,a)との関係を示す概略図は、それぞれ実施の形態１の説明において使用した図１、図３および図５と同一である。
【００８１】
図６は、実施の形態２のシステム同定装置における動的システムをＭ系列加振した場合のシステム入出力の時間波形を示す概略図である。
【００８２】
図６に示したように、実施の形態２に係るシステム同定装置１０では、同定対象とする動的システムに対してＭ系列信号を印加した場合のシステム入力１１（ｕ(jＴS)(ｊ=0,1,2,…)）およびシステム出力１２（ｙ(jＴS)(ｊ=0,1,2,…））に基づいて、当該システムを記述する線形離散時間システムを同定する。
【００８３】
図７は、実施の形態２のシステム同定装置におけるシステム次元決定部７の内部構成を示すブロック線図である。図７において、図４と同一符号を付したものは、実施の形態１と同一または同等の構成要素であり、システム安定性評価部３４が追加されている。
【００８４】
実施の形態２によるシステム次元決定部７では、図７に示すように、作業者が指定したシステム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属する次元ｎiのそれぞれに対して、再帰的システム行列推定部３１で同定したシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niに基づき、システム安定性評価部３４にて線形離散時間システムの安定性を評価する。
【００８５】
システム特性推定部３２は、システム安定性評価部３４にて安定システムと判断された次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niに基づいて、同定した線形離散時間システムに関する周波数応答を算出する。
【００８６】
システム次元推定部３３は、動的システムのシステム入出力から得られる実際の周波数応答（図７では動的システムのシステム特性と記載）に基づいて重み関数を決定し、システム特性推定部３２から出力される線形離散時間システムの周波数応答と、動的システムの実際の周波数応答との周波数領域における誤差２乗値に対して、当該重み関数を乗じた値の加算値ｅni（ｎi：安定システムとなる次元）を算出し、図５に示したように当該加算値のノルム||ｅni||の分布４１が予め設定された誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ6）。
【００８７】
次に、実施の形態２に係るシステム同定装置の動作について説明する。
【００８８】
同定対象とする動的システムが、１入力Ｐ出力のｎ次元線形離散時間システムとして［数１］で記述できるとする。この動的システムへのシステム入力ｕ(jＴS)をＭ系列信号で構成すると、当該システム入力ｕ(jＴS)および［数１］に対応するシステム出力ｙ(jＴS)は、例えば図６に示したシステム入力１１およびシステム出力１２のような時間波形となる。
【００８９】
実施の形態２に係るシステム同定装置１０では、図１および図６に示すように、システム入出力抽出部１は、予め設定された比率閾値とシステム入力１１（ｕ(jＴS)）の最大値とを乗じた［数２］をシステム入力閾値１３とし、システム入力１１の絶対値がシステム入力閾値１３以上となる時刻の最小値をＭ系列信号印加時刻ｊminＴsとして特定する（図６の例はｊminＴs＝２Ｔs）。
【００９０】
また、システム入出力抽出部１は、Ｍ系列信号印加時刻ｊminＴs以降のシステム入力１１およびシステム出力１２を［数３］で抽出し、抽出したそれぞれを同定用入力データｕid(jＴS)および同定用出力データｙid(jＴS)とすることで、対象とする動的システムのシステム入出力からＭ系列信号印加前のシステム静止時間領域データを除去する。
【００９１】
次に実施の形態１と同様に、ブロックハンケル行列生成部２は［数４］で与えられるブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfを生成し、入出力ベクトル生成部３は［数５］で与えられる動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kを生成し、ＬＱ分解部４はブロックハンケル行列Ｕp,Ｕf,Ｙp,Ｙfを結合したデータ行列（［数６］）をＬＱ分解し（［数７］）、部分行列Ｌ22,Ｌ32を出力する。
【００９２】
平行射影生成部５は、［数８］で定義される動的システムの平行射影Θを生成し、特異値分解部６は、生成された平行射影Θを特異値分解することで、［数９］で与えられる第１の直交行列Ｕ、第２の直交行列Ｖおよび、特異値σi(i=1,2,3, …)を出力する。
【００９３】
システム次元決定部７では、図７に示す処理が実行される。まず、再帰的システム行列推定部３１は、作業者が指定したシステム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属する次元ｎiのそれぞれに関して、［数１１］に示される再帰的手法によって、対応するシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niを同定する。
【００９４】
次にシステム安定性評価部３４は、作業者が指定したシステム次元の探索範囲ｎi=(ｎ1,ｎ2,…,ｎa)（ただし、ｎ1<ｎ2<…<ｎa）に属する次元ｎiのそれぞれに対して、再帰的システム行列推定部３１が同定したシステム行列Ａd,niに基づいて、線形離散時間システムの安定性を、以下の内容にて評価する。
【００９５】
【数１６】
0006009105.tif 000017
【００９６】
システム特性推定部３２は、システム安定性評価部３４で安定システムと判断された次元に対して、再帰的システム行列推定部３１が生成したシステム行列Ａd,ni,Ｂd,ni,Ｃd,ni,Ｄd,niに基づいて、同定された線形離散時間システムの周波数応答^Ｈni(kΔf)(k=0,1,2,…,N/2-1）を算出する。
【００９７】
実施の形態２のシステム同定装置１０では、“周波数領域において実際の周波数応答に最も適合する”という前提の下、システム次元推定部３３にて最適なシステム次元ｎを決定する。具体的には、以下の通りである。
【００９８】
まず、次式で与えられる同定用入出力データｕid(jＴS)，ｙid(jＴS)の有限離散フーリエ変換Ｕid(kΔf),Ｙid(kΔf)(k=0,1,2,…,N/2-1）から、次々式で得られる動的システムの実際の周波数応答Ｈ(kΔf)(k=0,1,2,…,N/2-1）（図７では動的システムのシステム特性と記載）を算出する。
【００９９】
【数１７】
0006009105.tif 000018
【０１００】
【数１８】
0006009105.tif 000019
【０１０１】
次に周波数応答Ｈ(kΔf)(k=0,1,2,…,N/2-1）に基づいて、例えば高ゲインかつ低周波領域に重みを持たせた次式で示されるような重み関数Ｗ(kΔf)(k=0,1,2,…,N/2-1）を決定する。
【０１０２】
【数１９】
0006009105.tif 000020
【０１０３】
そして、システム特性推定部３２から出力される線形離散時間システムの周波数応答^Ｈni(kΔf)と、動的システムの実際の周波数応答Ｈ(kΔf)との周波数領域における誤差２乗値に対して、重み関数Ｗ(kΔf)を乗じた値の加算値ｅni（ｎi：安定システムとなる次元）を次式にて算出する。
【０１０４】
【数２０】
0006009105.tif 000021
【０１０５】
この重み付き誤差２乗和のノルム||ｅni||を最小とする次元ｎiが、“重み関数に応じて、周波数領域において実際の周波数応答に最も適合する”安定なシステム次元ｎとなる。ここでは、図５に示すように重み付き誤差２乗和のノルム||ｅni||の分布４１が、［数１３］で与えられる誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の例は、システム次元ｎ＝ｎ6）。
【０１０６】
状態ベクトル生成部８は、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σi(i=1,2,3, …)と、システム次元決定部７から出力されるシステム次元ｎに基づいて、動的システムの状態ベクトル~ＸK,~ＸK+1を［数１４］を用いて生成する。
【０１０７】
最後にシステム行列同定部９は、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~ＵK|Kおよび出力ベクトル~ＹK|Kならびに、状態ベクトル生成部８から出力される動的システムの状態ベクトル~ＸK,~ＸK+1に基づいて、動的システムを記述する線形離散時間システムのシステム行列Ａd,Ｂd,Ｃd,Ｄdを［数１５］を用いて同定して出力する。
【０１０８】
このように、実施の形態２に係るシステム同定装置１０によれば、現実のシステム入出力から算出した平行射影Θの特異値σi(i=1,2,3, …)がなだらかな単調減少となり、したがって有意な値を持つ特異値と、同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して、周波数領域での重み関数に応じて一致度の高いシステム次元ｎを決定することができ、動的システムを記述する線形離散時間システムの同定が可能となる。
【０１０９】
また、動的システムの現実のシステム入出力から、Ｍ系列信号印加前のシステム静止時間領域データを除去することで、同定精度を向上させることも可能となる。
【０１１０】
さらに、再帰的システム行列推定部３１の存在により、現実の動的システムに対して一致度の高いシステム次元ｎを決定するための演算量を低減することが可能となる。
【０１１１】
加えて、システム安定性評価部３４の存在により、現実の動的システムが安定システムであることが明確である場合に、安定システムに限定した線形離散時間システムの同定が可能となる。
【０１１２】
なお、実施の形態２のシステム同定装置１０では、線形離散時間システムのシステム特性を周波数応答として算出し、当該周波数応答と、動的システムの同定用入出力データから得られる実際の周波数応答との周波数領域における誤差２乗和のノルム分布４１が予め設定された閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定している。しかしながら、本発明はこれに限定するものではなく、線形離散時間システムに対して実際の同定用入力データを適用した場合のシステム出力をシステム特性として算出し、当該システム出力と、動的システムの実際の同定用出力データとの時間領域における誤差２乗和に基づいて、システム次元ｎを決定してもよい。
【０１１３】
実施の形態３．
実施の形態３では、同定対象とする動的システムがＤＣサーボモータである場合について説明する。図８は、実施の形態３に係る全体構成を示すブロック線図である。本実施の形態において、図８に示すシステム同定装置１０は、図１に示した実施の形態１に係るシステム同定装置１０と同一または同等の構成をとる。本実施の形態では、ＤＣサーボモータ５１の入力電流[Ａｒｍｓ]として、例えばＭ系列信号などの疑似ランダム信号を入力し、当該疑似ランダム信号を動的システムのシステム入力１１（ｕ(jＴS)(ｊ=0,1,2,…)）とする。また、角速度[rad/s]を動的システムのシステム出力１２（ｙ(jＴS)(ｊ=0,1,2,…)）として取得する。システム同定装置１０は、これらのシステム入出力および、システム次元の探索範囲を入力とし、ＤＣサーボモータ５１を記述する線形離散時間システムを同定する。このとき、システム次元の探索範囲は、ｎi＝(1,2,…,50)のように、予測されるシステム次元に対して十分な幅を取って設定すればよい。システム同定装置１０では、現実の動的システムに対して一致度の高いシステム次元の決定と、動的システムを記述する線形離散時間システムの同定とを可能とするため、当該線形離散時間システムをサーボモータの制御系におけるパラメータ設計、フィルタのパラメータ設計などに利用することができる。
【符号の説明】
【０１１４】
１  システム入出力抽出部、２  ブロックハンケル行列生成部、３  入出力ベクトル生成部、４  ＬＱ分解部、５  平行射影生成部、６  特異値分解部、７  システム次元決定部、８  状態ベクトル生成部、９  システム行列同定部、１０  システム同定装置、１１  システム入力、１２  システム出力、１３  システム入力閾値、２１  （理想的なシステム入出力における平行射影の）特異値分布、２２  （現実のシステム入出力における平行射影の）特異値分布、３１  再帰的システム行列推定部、３２  システム特性推定部、３３  システム次元推定部、３４  システム安定性評価部、４１  （時間領域または周波数領域における）誤差２乗和のノルム分布、４２  （時間領域または周波数領域における）誤差２乗和ノルム閾値、５１  ＤＣサーボモータ。
【図１】
0006009105.tif 000022
【図２】
0006009105.tif 000023
【図３】
0006009105.tif 000024
【図４】
0006009105.tif 000025
【図５】
0006009105.tif 000026
【図６】
0006009105.tif 000027
【図７】
0006009105.tif 000028
【図８】
0006009105.tif 000029