JPWO2015118736A1 - システム同定装置 - Google Patents

Abstract

システム同定装置は、指定されたシステム次元の探索範囲に属するそれぞれの次元に対して再帰的手法によって線形離散時間システムを同定し、各次元に対応する線形離散時間システムに対して、実際の同定用入力データを適用した場合のシステム出力をシステム特性として算出し、当該システム出力と、動的システムの実際の同定用出力データとの時間領域における誤差２乗和のノルム分布（４１）が閾値（４２）以下となる次元のうち、最小の次元をシステム次元ｎとして決定すると共に、動的システムの入出力ベクトルおよび、決定されたシステム次元を用いて生成された状態ベクトルに基づいて線形離散時間システムのシステム行列を同定する。

Description

この発明は、対象とする動的システムに対して擬似ランダム入力を印加した場合のシステム入出力に基づいて、当該システムの数学モデルを構築するシステム同定装置に関する。

従来の擬似ランダム入力によるシステム同定装置として、例えば非特許文献１に記載されたＮ４ＳＩＤ法に基づくものがある。このＮ４ＳＩＤ法では、線形離散時間システム（Ａ_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を結合したデータ行列を生成し、当該データ行列をＬＱ分解した部分行列Ｌ₂₂,Ｌ₃₂を生成して出力する。

平行射影生成部５は、ＬＱ分解部４から出力される部分行列Ｌ₂₂,Ｌ₃₂と、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｕ_p,Ｙ_pに基づいて、動的システムの平行射影Θを生成する。

特異値分解部６は、平行射影生成部５から出力される平行射影Θを特異値分解し、平行射影Θの左特異ベクトルを列ベクトルとする第１の直交行列Ｕ、平行射影Θの右特異ベクトルを列ベクトルとする第２の直交行列Ｖおよび、平行射影Θの特異値σ_i(i=1,2,3…)を出力する。

システム次元決定部７は、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σ_i(i=1,2,3…)と、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~Ｕ_K|Kおよび出力ベクトル~Ｙ_K|Kならびに、作業者が指定したシステム次元の探索範囲ｎ_i=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_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…)および、システム次元決定部７から出力されるシステム次元ｎに基づいて、動的システムの状態ベクトル~Ｘ_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(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_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,n_i,Ｂ_d,n_i,Ｃ_d,n_i,Ｄ_d,n_iを同定する。

次にシステム特性推定部３２は、システム次元の探索範囲ｎ_i=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_a）に属するそれぞれの次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａ_d,n_i,Ｂ_d,n_i,Ｃ_d,n_i,Ｄ_d,n_iに基づいて、同定した線形離散時間システムに対して、実際の同定用入力データｕ_id(jＴ_S)(ｊ=0,1,2,…)を適用した場合のシステム出力を算出する。

なお、再帰的システム行列推定部３１およびシステム特性推定部３２の処理は、ｉをインクリメントし、ｉ＝ａとなるまで実行される。

システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムのシステム出力と、動的システムの実際の同定用出力データｙ_id(jＴ_S)(ｊ=0,1,2,…)（図４では動的システムのシステム特性と記載）との時間領域における誤差２乗和ｅ_ni(i=1,2,…,ａ)を算出し、図５に示したように誤差２乗和のノルム||ｅ_ni||の分布４１が、予め設定された誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する構成としている（図５の場合はシステム次元ｎ＝ｎ₆）。

次に、実施の形態１に係るシステム同定装置の動作について説明する。

同定対象とする動的システムが、次式のように、１入力Ｐ出力のｎ次元線形離散時間システムで記述できるとする。

上記の動的システムへのシステム入力ｕ(jＴ_S)を擬似ランダム入力で構成すると、当該システム入力ｕ(jＴ_S)および、上記［数１］式に対応するシステム出力ｙ(jＴ_S)は、例えば図２に示したシステム入力１１、およびシステム出力１２のような時間波形となる。

ここで、図１および図２の説明で述べたように、予め設定された比率閾値とシステム入力１１（ｕ(jＴ_S)）の最大値とを乗じた次式をシステム入力閾値１３として使用する。

システム入出力抽出部１は、システム入力１１の絶対値がシステム入力閾値１３以上となる時刻の最小値を擬似ランダム入力印加時刻ｊ_minＴ_Sとして特定する（図２の場合はｊ_minＴ_S＝３Ｔ_S）。

また、システム入出力抽出部１は、擬似ランダム入力印加時刻ｊ_minＴ_S以降のシステム入力１１およびシステム出力１２を次式を用いて抽出する。

さらに、システム入出力抽出部１は、上記［数３］を用いて抽出したそれぞれの値を同定用入力データｕ_id(ｊＴ_S)および同定用出力データｙ_id(ｊＴ_S)とすることで、対象とする動的システムのシステム入出力から擬似ランダム入力印加前のシステム静止時間領域データを除去する。

ブロックハンケル行列生成部２は、システム入出力抽出部１から出力される同定用入力データｕ_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を結合した次式で与えられるデータ行列を生成する。

また、ＬＱ分解部４は、上記データ行列を次式のようにＬＱ分解すると共に、ＬＱ分解した行列の要素から、部分行列Ｌ₂₂,Ｌ₃₂を出力する。

平行射影生成部５は、ＬＱ分解部４から出力される部分行列Ｌ₂₂,Ｌ₃₂と、ブロックハンケル行列生成部２から出力されるブロックハンケル行列Ｕ_p,Ｙ_pに基づいて、次式で定義される動的システムの平行射影Θを生成する。

特異値分解部６は、上式で示される平行射影Θを特異値分解することで、次式で与えられる平行射影Θの左特異ベクトルｕ_jを列ベクトルとする第１の直交行列Ｕ、平行射影Θの右特異ベクトルｖ_jを列ベクトルとする第２の直交行列Ｖおよび、平行射影Θの特異値σ_i(i=1,2,3…)を出力する。

対象とする動的システムのシステム次元ｎは、平行射影Θの特異値において有意な値を持つものはｎ個であり、それらに比してｎ＋１個目以降は十分小さいとする以下の関係に基づいて決定することができる。

図３に示したように、動的システムのシステム入出力から算出した平行射影Θの特異値σ_iは、理想的には次元(i=1,2,3…)に対して、例えば特異値分布２１に示す関係となる。この場合、有意な値を持つ特異値の個数を明確に規定することができ、当該個数から動的システムのシステム次元ｎを決定することができる（図３の場合はシステム次元ｎ＝４）。一方、観測ノイズ等の影響を受ける現実のシステム入出力に基づいて算出した特異値σ_iは、次元(i=1,2,3…)に対して、例えば特異値分布２２に示す関係となるため、有意な値を持つ特異値と、無視できる微小な値となる特異値との境界σ_n>>σ_n+1が不明確となる。したがって、従来の手法では、常に最適なシステム次元ｎを決定しているとは限らず、最適なシステム次元ｎを決定するためには、試行錯誤が必要になるという課題があった。

そこで、実施の形態１に係るシステム同定装置１０では、“時間領域において実際のシステム入出力に最も適合する”という前提の下、システム次元決定部７において最適なシステム次元ｎを決定する。システム次元決定部７では、図１に示したように特異値分解部６から出力される第２の直交行列Ｖおよび特異値σ_i(i=1,2,3…)と、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~Ｕ_K|Kおよび出力ベクトル~Ｙ_K|Kならびに、作業者が指定したシステム次元の探索範囲ｎ_i=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_a）に基づいて、当該探索範囲に属する次元ｎ_i(i=1,2,…,a)のそれぞれに対して動的システムを記述する線形離散時間システムのシステム行列を同定する。さらに当該システム行列に基づいて、当該探索範囲に属する次元ｎ_i(i=1,2,…,a)のそれぞれに対応する線形離散時間システムに対して、実際の同定用入力データｕ_id(jＴ_S)(ｊ=0,1,2,…)を適用した場合のシステム出力を算出し、動的システムの実際の同定用出力データｙ_id(jＴ_S)(ｊ=0,1,2,…）（図１では動的システムのシステム特性と記載）との比較から、システム次元ｎを決定する。

具体的には図４に示したように、再帰的システム行列推定部３１は、作業者が指定したシステム次元の探索範囲ｎ_i(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_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,n_i,Ｂ_d,n_i,Ｃ_d,n_i,Ｄ_d,n_iを同定する。

システム特性推定部３２は、システム次元の探索範囲ｎ_i=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_a）に属する次元ｎ_iのそれぞれに対して、再帰的システム行列推定部３１から出力されるシステム行列Ａ_d,n_i,Ｂ_d,n_i,Ｃ_d,n_i,Ｄ_d,n_iに基づいて、同定した線形離散時間システムに対して実際の同定用入力データｕ_id(jＴ_S)(ｊ=0,1,2,…)（［数３］参照）を適用した場合のシステム出力^ｙ_id,n_i(jＴ_S)(ｊ=0,1,2,…）を算出する。なお、「^ｙ」の表記は、「ｙ」という文字の上部に「^」の記号を付したことを意味する代替表記である。

また、システム次元推定部３３は、システム特性推定部３２から出力される線形離散時間システムのシステム出力^ｙ_id,n_i(jＴ_S)(ｊ=0,1,2,…）と、動的システムの実際の同定用出力データｙ_id(jＴ_S)(ｊ=0,1,2,…）（図４では動的システムのシステム特性と記載）との時間領域における誤差２乗和ｅn_i(i=1,2,…,a）を次式を用いて算出する。

上式で示される誤差２乗和のノルム||ｅn_i||を最小とする次元ｎ_iが、“時間領域において実際のシステム入出力に最も適合する”システム次元ｎとなる。一方、実際のノルム||ｅn_i||は、観測ノイズが白色性であれば、そのノイズレベルによらず次元ｎ_iの増加に伴って単調減少し、図５に示すようにある次元以上でほぼ一定値となる。そこでここでは、システム次元ｎの推定値が必要以上に高次元となることを避けるため、次式で与えられる誤差２乗和ノルム閾値４２を規定する。

システム次元推定部３３は、誤差２乗和のノルム||ｅn_i||の分布４１が上記誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の例では、システム次元ｎ＝ｎ₆）。

状態ベクトル生成部８は、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σ_i(i=1,2,3…)と、システム次元決定部７から出力されるシステム次元ｎに基づいて、動的システムの状態ベクトル~Ｘ_K,~Ｘ_K+1を次式に従って生成する。

最後に、システム行列同定部９は、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~Ｕ_K|Kおよび出力ベクトル~Ｙ_K|Kならびに、状態ベクトル生成部８から出力される動的システムの状態ベクトル~Ｘ_K,~Ｘ_K+1に基づいて、動的システムを記述する線形離散時間システムのシステム行列Ａ_d,Ｂ_d,Ｃ_d,Ｄ_dを次式を用いて同定して出力する。

このように、実施の形態１に係るシステム同定装置１０によれば、現実のシステム入出力から算出した平行射影Θの特異値σ_i(i=1,2,3…)がなだらかな単調減少となり、したがって有意な値を持つ特異値と、同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して時間領域で一致度の高いシステム次元ｎを決定することができ、動的システムを記述する線形離散時間システムの同定が可能となる。

また、動的システムの現実のシステム入出力から、擬似ランダム入力印加前のシステム静止時間領域データを除去することで、同定精度を向上させることも可能となる。

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

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

実施の形態２．
つぎに、実施の形態２に係るシステム同定装置について説明する。なお、実施の形態２に係るシステム同定装置の全体構成を示すブロック線図、平行射影Θの特異値σ_iと次元(i=1,2,3…)の関係を示す概略図および、同定した線形離散時間システムの周波数応答と動的システムの実際の周波数応答との周波数領域における誤差２乗和のノルム||ｅn_i||と次元ｎ_i(i=1,2,…,a)との関係を示す概略図は、それぞれ実施の形態１の説明において使用した図１、図３および図５と同一である。

図６は、実施の形態２のシステム同定装置における動的システムをＭ系列加振した場合のシステム入出力の時間波形を示す概略図である。

図６に示したように、実施の形態２に係るシステム同定装置１０では、同定対象とする動的システムに対してＭ系列信号を印加した場合のシステム入力１１（ｕ(jＴ_S)(ｊ=0,1,2,…)）およびシステム出力１２（ｙ(jＴ_S)(ｊ=0,1,2,…））に基づいて、当該システムを記述する線形離散時間システムを同定する。

図７は、実施の形態２のシステム同定装置におけるシステム次元決定部７の内部構成を示すブロック線図である。図７において、図４と同一符号を付したものは、実施の形態１と同一または同等の構成要素であり、システム安定性評価部３４が追加されている。

実施の形態２によるシステム次元決定部７では、図７に示すように、作業者が指定したシステム次元の探索範囲ｎ_i=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_a）に属する次元ｎ_iのそれぞれに対して、再帰的システム行列推定部３１で同定したシステム行列Ａ_d,n_i,Ｂ_d,n_i,Ｃ_d,n_i,Ｄ_d,n_iに基づき、システム安定性評価部３４にて線形離散時間システムの安定性を評価する。

システム特性推定部３２は、システム安定性評価部３４にて安定システムと判断された次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａ_d,n_i,Ｂ_d,n_i,Ｃ_d,n_i,Ｄ_d,n_iに基づいて、同定した線形離散時間システムに関する周波数応答を算出する。

システム次元推定部３３は、動的システムのシステム入出力から得られる実際の周波数応答（図７では動的システムのシステム特性と記載）に基づいて重み関数を決定し、システム特性推定部３２から出力される線形離散時間システムの周波数応答と、動的システムの実際の周波数応答との周波数領域における誤差２乗値に対して、当該重み関数を乗じた値の加算値ｅn_i（ｎ_i：安定システムとなる次元）を算出し、図５に示したように当該加算値のノルム||ｅn_i||の分布４１が予め設定された誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ₆）。

次に、実施の形態２に係るシステム同定装置の動作について説明する。

同定対象とする動的システムが、１入力Ｐ出力のｎ次元線形離散時間システムとして［数１］で記述できるとする。この動的システムへのシステム入力ｕ(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を結合したデータ行列（［数６］）をＬＱ分解し（［数７］）、部分行列Ｌ₂₂,Ｌ₃₂を出力する。

平行射影生成部５は、［数８］で定義される動的システムの平行射影Θを生成し、特異値分解部６は、生成された平行射影Θを特異値分解することで、［数９］で与えられる第１の直交行列Ｕ、第２の直交行列Ｖおよび、特異値σ_i(i=1,2,3…)を出力する。

システム次元決定部７では、図７に示す処理が実行される。まず、再帰的システム行列推定部３１は、作業者が指定したシステム次元の探索範囲ｎ_i=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_a）に属する次元ｎ_iのそれぞれに関して、［数１１］に示される再帰的手法によって、対応するシステム行列Ａ_d,n_i,Ｂ_d,n_i,Ｃ_d,n_i,Ｄ_d,n_iを同定する。

次にシステム安定性評価部３４は、作業者が指定したシステム次元の探索範囲ｎ_i=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_a）に属する次元ｎ_iのそれぞれに対して、再帰的システム行列推定部３１が同定したシステム行列Ａ_d,n_iに基づいて、線形離散時間システムの安定性を、以下の内容にて評価する。

システム特性推定部３２は、システム安定性評価部３４で安定システムと判断された次元に対して、再帰的システム行列推定部３１が生成したシステム行列Ａ_d,n_i,Ｂ_d,n_i,Ｃ_d,n_i,Ｄ_d,n_iに基づいて、同定された線形離散時間システムの周波数応答^Ｈn_i(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）（図７では動的システムのシステム特性と記載）を算出する。

次に周波数応答Ｈ(kΔf)(k=0,1,2,…,N/2-1）に基づいて、例えば高ゲインかつ低周波領域に重みを持たせた次式で示されるような重み関数Ｗ(kΔf)(k=0,1,2,…,N/2-1）を決定する。

そして、システム特性推定部３２から出力される線形離散時間システムの周波数応答^Ｈn_i(kΔf)と、動的システムの実際の周波数応答Ｈ(kΔf)との周波数領域における誤差２乗値に対して、重み関数Ｗ(kΔf)を乗じた値の加算値ｅn_i（ｎ_i：安定システムとなる次元）を次式にて算出する。

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

状態ベクトル生成部８は、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σ_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)のように、予測されるシステム次元に対して十分な幅を取って設定すればよい。システム同定装置１０では、現実の動的システムに対して一致度の高いシステム次元の決定と、動的システムを記述する線形離散時間システムの同定とを可能とするため、当該線形離散時間システムをサーボモータの制御系におけるパラメータ設計、フィルタのパラメータ設計などに利用することができる。

１ システム入出力抽出部、２ ブロックハンケル行列生成部、３ 入出力ベクトル生成部、４ ＬＱ分解部、５ 平行射影生成部、６ 特異値分解部、７ システム次元決定部、８ 状態ベクトル生成部、９ システム行列同定部、１０ システム同定装置、１１ システム入力、１２ システム出力、１３ システム入力閾値、２１ （理想的なシステム入出力における平行射影の）特異値分布、２２ （現実のシステム入出力における平行射影の）特異値分布、３１ 再帰的システム行列推定部、３２ システム特性推定部、３３ システム次元推定部、３４ システム安定性評価部、４１ （時間領域または周波数領域における）誤差２乗和のノルム分布、４２ （時間領域または周波数領域における）誤差２乗和ノルム閾値、５１ ＤＣサーボモータ。

上述した課題を解決し、目的を達成するため、本発明に係るシステム同定装置は、同定対象とする動的システムに対して擬似ランダム入力を印加した場合のシステム入出力を入力とするシステム同定装置において、前記動的システムのシステム入出力から同定に適用する同定用入出力データを抽出するシステム入出力抽出部と、前記同定用入出力データに基づいて、ブロックハンケル行列を生成するブロックハンケル行列生成部と、前記ブロックハンケル行列に基づいて、前記動的システムの入力ベクトルおよび出力ベクトルを生成する入出力ベクトル生成部と、前記ブロックハンケル行列を結合してデータ行列を生成し、当該データ行列をＬＱ分解した部分行列を出力するＬＱ分解部と、前記部分行列と前記ブロックハンケル行列に基づいて、平行射影を生成する平行射影生成部と、前記平行射影の特異値分解により、前記平行射影の左特異ベクトルを列ベクトルとする第１の直交行列、当該平行射影の右特異ベクトルを列ベクトルとする第２の直交行列および当該平行射影の特異値を出力する特異値分解部と、前記第２の直交行列および前記特異値と、前記動的システムの入力ベクトルおよび出力ベクトルと、指定されたシステム次元の探索範囲に基づいて、当該探索範囲に属するそれぞれの次元に対して動的システムを記述する線形離散時間システムのシステム行列を同定し、さらに当該システム行列に基づいて算出した線形離散時間システムのシステム特性と、動的システムの実際のシステム特性との比較から、システム次元を決定するシステム次元決定部と、前記第２の直交行列および特異値と、前記決定されたシステム次元に基づいて、前記動的システムの状態ベクトルを生成する状態ベクトル生成部と、前記動的システムの入力ベクトルおよび出力ベクトルならびに前記動的システムの状態ベクトルに基づいて、当該動的システムを記述する線形離散時間システムのシステム行列を同定するシステム行列同定部と、を備え、前記同定されたシステム行列を、前記動的システムを記述する線形離散時間システムとして出力することを特徴とする。

特異値分解部６は、平行射影生成部５から出力される平行射影Θを特異値分解し、平行射影Θの左特異ベクトルを列ベクトルとする第１の直交行列Ｕ、平行射影Θの右特異ベクトルを列ベクトルとする第２の直交行列Ｖおよび、平行射影Θの特異値σ_i(i=1,2,3, …)を出力する。

システム次元決定部７は、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σ_i(i=1,2,3, …)と、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~Ｕ_K|Kおよび出力ベクトル~Ｙ_K|Kならびに、作業者が指定したシステム次元の探索範囲ｎ_i=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_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を生成する。

図３は、実施の形態１のシステム同定装置１０における平行射影Θの特異値σ_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(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_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=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_a）に属するそれぞれの次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａ _d,ni ,Ｂ _d,ni ,Ｃ _d,ni ,Ｄ _d,ni に基づいて、同定した線形離散時間システムに対して、実際の同定用入力データｕ_id(jＴ_S)(ｊ=0,1,2,…)を適用した場合のシステム出力を算出する。

特異値分解部６は、上式で示される平行射影Θを特異値分解することで、次式で与えられる平行射影Θの左特異ベクトルｕ_jを列ベクトルとする第１の直交行列Ｕ、平行射影Θの右特異ベクトルｖ_jを列ベクトルとする第２の直交行列Ｖおよび、平行射影Θの特異値σ_i(i=1,2,3, …)を出力する。

図３に示したように、動的システムのシステム入出力から算出した平行射影Θの特異値σ_iは、理想的には次元(i=1,2,3, …)に対して、例えば特異値分布２１に示す関係となる。この場合、有意な値を持つ特異値の個数を明確に規定することができ、当該個数から動的システムのシステム次元ｎを決定することができる（図３の場合はシステム次元ｎ＝４）。一方、観測ノイズ等の影響を受ける現実のシステム入出力に基づいて算出した特異値σ_iは、次元(i=1,2,3, …)に対して、例えば特異値分布２２に示す関係となるため、有意な値を持つ特異値と、無視できる微小な値となる特異値との境界σ_n>>σ_n+1が不明確となる。したがって、従来の手法では、常に最適なシステム次元ｎを決定しているとは限らず、最適なシステム次元ｎを決定するためには、試行錯誤が必要になるという課題があった。

そこで、実施の形態１に係るシステム同定装置１０では、“時間領域において実際のシステム入出力に最も適合する”という前提の下、システム次元決定部７において最適なシステム次元ｎを決定する。システム次元決定部７では、図１に示したように特異値分解部６から出力される第２の直交行列Ｖおよび特異値σ_i(i=1,2,3, …)と、入出力ベクトル生成部３から出力される動的システムの入力ベクトル~Ｕ_K|Kおよび出力ベクトル~Ｙ_K|Kならびに、作業者が指定したシステム次元の探索範囲ｎ_i=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_a）とに基づいて、当該探索範囲に属する次元ｎ_i(i=1,2,…,a)のそれぞれに対して動的システムを記述する線形離散時間システムのシステム行列を同定する。さらに当該システム行列に基づいて、当該探索範囲に属する次元ｎ_i(i=1,2,…,a)のそれぞれに対応する線形離散時間システムに対して、実際の同定用入力データｕ_id(jＴ_S)(ｊ=0,1,2,…)を適用した場合のシステム出力を算出し、動的システムの実際の同定用出力データｙ_id(jＴ_S)(ｊ=0,1,2,…）（図１では動的システムのシステム特性と記載）との比較から、システム次元ｎを決定する。

具体的には図４に示したように、再帰的システム行列推定部３１は、作業者が指定したシステム次元の探索範囲ｎ_i(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_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=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_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）を次式を用いて算出する。

上式で示される誤差２乗和のノルム||ｅ _ni ||を最小とする次元ｎ_iが、“時間領域において実際のシステム入出力に最も適合する”システム次元ｎとなる。一方、実際のノルム||ｅ _ni ||は、観測ノイズが白色性であれば、そのノイズレベルによらず次元ｎ_iの増加に伴って単調減少し、図５に示すようにある次元以上でほぼ一定値となる。そこでここでは、システム次元ｎの推定値が必要以上に高次元となることを避けるため、次式で与えられる誤差２乗和ノルム閾値４２を規定する。

システム次元推定部３３は、誤差２乗和のノルム||ｅ _ni ||の分布４１が上記誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の例では、システム次元ｎ＝ｎ₆）。

状態ベクトル生成部８は、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σ_i(i=1,2,3, …)と、システム次元決定部７から出力されるシステム次元ｎに基づいて、動的システムの状態ベクトル~Ｘ_K,~Ｘ_K+1を次式に従って生成する。

このように、実施の形態１に係るシステム同定装置１０によれば、現実のシステム入出力から算出した平行射影Θの特異値σ_i(i=1,2,3, …)がなだらかな単調減少となり、したがって有意な値を持つ特異値と、同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して時間領域で一致度の高いシステム次元ｎを決定することができ、動的システムを記述する線形離散時間システムの同定が可能となる。

実施の形態２．
つぎに、実施の形態２に係るシステム同定装置について説明する。なお、実施の形態２に係るシステム同定装置の全体構成を示すブロック線図、平行射影Θの特異値σ_iと次元(i=1,2,3, …)の関係を示す概略図および、同定した線形離散時間システムの周波数応答と動的システムの実際の周波数応答との周波数領域における誤差２乗和のノルム||ｅ _ni ||と次元ｎ_i(i=1,2,…,a)との関係を示す概略図は、それぞれ実施の形態１の説明において使用した図１、図３および図５と同一である。

実施の形態２によるシステム次元決定部７では、図７に示すように、作業者が指定したシステム次元の探索範囲ｎ_i=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_a）に属する次元ｎ_iのそれぞれに対して、再帰的システム行列推定部３１で同定したシステム行列Ａ _d,ni ,Ｂ _d,ni ,Ｃ _d,ni ,Ｄ _d,ni に基づき、システム安定性評価部３４にて線形離散時間システムの安定性を評価する。

システム特性推定部３２は、システム安定性評価部３４にて安定システムと判断された次元に対して、再帰的システム行列推定部３１から出力されるシステム行列Ａ _d,ni ,Ｂ _d,ni ,Ｃ _d,ni ,Ｄ _d,ni に基づいて、同定した線形離散時間システムに関する周波数応答を算出する。

システム次元推定部３３は、動的システムのシステム入出力から得られる実際の周波数応答（図７では動的システムのシステム特性と記載）に基づいて重み関数を決定し、システム特性推定部３２から出力される線形離散時間システムの周波数応答と、動的システムの実際の周波数応答との周波数領域における誤差２乗値に対して、当該重み関数を乗じた値の加算値ｅ _ni （ｎ_i：安定システムとなる次元）を算出し、図５に示したように当該加算値のノルム||ｅ _ni ||の分布４１が予め設定された誤差２乗和ノルム閾値４２以下となる次元のうち、最小の次元をシステム次元ｎとして決定して出力する（図５の場合はシステム次元ｎ＝ｎ₆）。

平行射影生成部５は、［数８］で定義される動的システムの平行射影Θを生成し、特異値分解部６は、生成された平行射影Θを特異値分解することで、［数９］で与えられる第１の直交行列Ｕ、第２の直交行列Ｖおよび、特異値σ_i(i=1,2,3, …)を出力する。

システム次元決定部７では、図７に示す処理が実行される。まず、再帰的システム行列推定部３１は、作業者が指定したシステム次元の探索範囲ｎ_i=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_a）に属する次元ｎ_iのそれぞれに関して、［数１１］に示される再帰的手法によって、対応するシステム行列Ａ _d,ni ,Ｂ _d,ni ,Ｃ _d,ni ,Ｄ _d,ni を同定する。

次にシステム安定性評価部３４は、作業者が指定したシステム次元の探索範囲ｎ_i=(ｎ₁,ｎ₂,…,ｎ_a)（ただし、ｎ₁<ｎ₂<…<ｎ_a）に属する次元ｎ_iのそれぞれに対して、再帰的システム行列推定部３１が同定したシステム行列Ａ _d,ni に基づいて、線形離散時間システムの安定性を、以下の内容にて評価する。

システム特性推定部３２は、システム安定性評価部３４で安定システムと判断された次元に対して、再帰的システム行列推定部３１が生成したシステム行列Ａ _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）（図７では動的システムのシステム特性と記載）を算出する。

そして、システム特性推定部３２から出力される線形離散時間システムの周波数応答^Ｈ _ni (kΔf)と、動的システムの実際の周波数応答Ｈ(kΔf)との周波数領域における誤差２乗値に対して、重み関数Ｗ(kΔf)を乗じた値の加算値ｅ _ni （ｎ_i：安定システムとなる次元）を次式にて算出する。

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

状態ベクトル生成部８は、特異値分解部６から出力される第２の直交行列Ｖおよび特異値σ_i(i=1,2,3, …)と、システム次元決定部７から出力されるシステム次元ｎに基づいて、動的システムの状態ベクトル~Ｘ_K,~Ｘ_K+1を［数１４］を用いて生成する。

このように、実施の形態２に係るシステム同定装置１０によれば、現実のシステム入出力から算出した平行射影Θの特異値σ_i(i=1,2,3, …)がなだらかな単調減少となり、したがって有意な値を持つ特異値と、同定において無視できる微小な値となる特異値との境界が不明確となる場合においても、システム次元ｎの決定から試行錯誤を排除し、現実の動的システムに対して、周波数領域での重み関数に応じて一致度の高いシステム次元ｎを決定することができ、動的システムを記述する線形離散時間システムの同定が可能となる。

Claims (7)

1. 同定対象とする動的システムに対して擬似ランダム入力を印加した場合のシステム入出力および指定されたシステム次元の探索範囲を入力とするシステム同定装置において、
   前記動的システムのシステム入出力から同定に適用する同定用入出力データを抽出するシステム入出力抽出部と、
   前記同定用入出力データに基づいて、ブロックハンケル行列を生成するブロックハンケル行列生成部と、
   前記ブロックハンケル行列に基づいて、前記動的システムの入力ベクトルおよび出力ベクトルを生成する入出力ベクトル生成部と、
   前記ブロックハンケル行列を結合してデータ行列を生成し、当該データ行列をＬＱ分解した部分行列を出力するＬＱ分解部と、
   前記部分行列と前記ブロックハンケル行列に基づいて、平行射影を生成する平行射影生成部と、
   前記平行射影の特異値分解により、前記平行射影の左特異ベクトルを列ベクトルとする第１の直交行列、当該平行射影の右特異ベクトルを列ベクトルとする第２の直交行列および当該平行射影の特異値を出力する特異値分解部と、
   前記第２の直交行列および前記特異値と、前記動的システムの入力ベクトルおよび出力ベクトルと、前記探索範囲に基づいて、当該探索範囲に属するそれぞれの次元に対して動的システムを記述する線形離散時間システムのシステム行列を同定し、さらに当該システム行列に基づいて算出した線形離散時間システムのシステム特性と、動的システムの実際のシステム特性との比較から、システム次元を決定するシステム次元決定部と、
   前記第２の直交行列および特異値と、前記決定されたシステム次元に基づいて、前記動的システムの状態ベクトルを生成する状態ベクトル生成部と、
   前記動的システムの入力ベクトルおよび出力ベクトルならびに前記動的システムの状態ベクトルに基づいて、当該動的システムを記述する線形離散時間システムのシステム行列を同定するシステム行列同定部と、を備え、
   前記同定されたシステム行列を、前記動的システムを記述する線形離散時間システムとして出力することを特徴とするシステム同定装置。
2. 前記システム次元決定部は、
   前記探索範囲に属するそれぞれの次元に対して、
   同定した線形離散時間システムに関して実際の同定用入力データを適用した場合のシステム出力を算出し、当該システム出力を線形離散時間システムのシステム特性として出力するシステム特性推定部と、
   前記線形離散時間システムのシステム出力と、前記動的システムの実際の同定用出力データとの時間領域における誤差２乗和のノルムが設定閾値以下となる次元のうち、最小の次元をシステム次元として決定して出力するシステム次元推定部と
   を備えたことを特徴とする請求項１に記載のシステム同定装置。
3. 前記システム次元決定部が、
   前記探索範囲に属するそれぞれの次元に対して、
   同定した線形離散時間システムの周波数応答を算出し、当該周波数応答を線形離散時間システムのシステム特性として出力するシステム特性推定部と、
   前記線形離散時間システムの周波数応答と、前記動的システムのシステム入出力から得られる実際の周波数応答との周波数領域における誤差２乗和のノルムが設定閾値以下となる次元のうち、最小の次元をシステム次元として決定して出力するシステム次元推定部と、
   を備えたことを特徴とする請求項１に記載のシステム同定装置。
4. 前記システム次元推定部は、
   前記動的システムのシステム入出力から得られる実際の周波数応答に基づいて重み関数を決定し、
   前記線形離散時間システムの周波数応答と、前記動的システムの実際の周波数応答との周波数領域における誤差２乗値に対して、当該重み関数を乗じた値の加算値を算出し、
   当該加算値のノルムが設定閾値以下となる次元のうち、最小の次元をシステム次元として決定して出力する
   ことを特徴とする請求項３に記載のシステム同定装置。
5. 前記システム次元決定部は、
   前記探索範囲に属する第１の次元に対応するシステム行列の同定に関して、
   当該探索範囲の中で第１の次元よりも１段階低い第２の次元に対応したシステム行列の同定結果と、
   前記第２の直交行列および特異値のうち、それぞれ前記第２の次元よりも大きく、前記第１の次元以下となる右特異ベクトルおよび特異値、ならびに、前記動的システムの入力ベクトルおよび出力ベクトルを用いて、再帰的手法により前記第１の次元に対応するシステム行列を同定する再帰的システム行列推定部と、
   を備えたことを特徴とする請求項１に記載のシステム同定装置。
6. 前記システム入出力抽出部は、
   設定された比率閾値とシステム入力の最大値とを乗じた値をシステム入力閾値とし、システム入力の絶対値がシステム入力閾値以上となる時刻の最小値を擬似ランダム入力印加時刻として、擬似ランダム入力印加時刻以降のシステム入力およびシステム出力を、それぞれ同定用入力データおよび同定用出力データとして抽出することを特徴とする請求項１に記載のシステム同定装置。
7. 前記システム次元決定部は、
   前記探索範囲に属するそれぞれの次元に対して線形離散時間システムの安定性を評価するシステム安定性評価部を備え、
   安定システムとなる次元に対応する線形離散時間システムのシステム特性から、システム次元を決定することを特徴とする請求項１に記載のシステム同定装置。

Images