(57)【要約】【課題】精度よく推定対象の状態を推定することができる。【解決手段】カルマンフィルタ部３２により、状態空間モデルで表わされるカルマンフィルタに従って、入力された観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、状態推定誤差を算出し、ロバストエンハンスメント部３４により、状態推定値に基づく観測値の残差の二乗の期待値を算出し、算出された期待値が予め定められた値以上であるときに、算出された期待値に基づいて、状態推定手段によって算出された状態推定値及び状態推定誤差の各々を補正し、出力部５０により、時系列データの観測値の各々について、カルマンフィルタ部３２による算出及びロバストエンハンスメント部３４による補正を繰り返したときに得られる状態推定値の各々を出力する。【選択図】図２
2014215822.tif 000002
【特許請求の範囲】
【請求項１】
カルマンフィルタ手段と、ロバストエンハンスメント手段と、出力手段とを含み、入力された推定対象の観測値の時系列データに基づいて、前記推定対象の状態を推定する状態推定装置における状態推定方法において、
前記カルマンフィルタ手段は、システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて予め定められた状態空間モデルで表わされるカルマンフィルタに従って、入力された前記観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、前記状態推定誤差を算出し、
前記ロバストエンハンスメント手段は、前記状態推定手段によって算出された前記状態推定値に基づく前記観測値の残差の二乗の期待値を算出し、前記算出された前記期待値が予め定められた値以上であるときに、前記算出された前記期待値に基づいて、前記状態推定手段によって算出された前記状態推定値及び前記状態推定誤差の各々を補正し、
前記出力手段は、前記時系列データの観測値の各々について、前記カルマンフィルタ手段による算出及び前記ロバストエンハンスメント手段による補正を繰り返したときに得られる前記状態推定値の各々を出力する
状態推定方法。
【請求項２】
前記カルマンフィルタ手段は、以下の（１）式〜（５）式に従って、前記状態推定値^Ｘtを算出すると共に、前記状態推定誤差^Ｖtを算出し、
前記ロバストエンハンスメント手段は、以下の（８）式に従って、前記観測値の残差の二乗の期待値γtを算出し、前記算出された前記期待値γtが予め定められた値τ以上であるときに、以下の（６）式、（７）式に従って、前記状態推定手段によって算出された前記状態推定値^Ｘt及び前記状態推定誤差^Ｖtの各々を状態推定値＾Ｘ´ｔ及び状態推定誤差＾Ｖ´ｔに補正する請求項１記載の状態推定方法。
【数１】
2014215822.tif 000013
ただし、Ｆは状態遷移行列であり、Ｈは観測出力行列であり、Ｒ及びＱは共分散行列であり、ｖｉは時刻ｉにおける前記観測値の残差を示すベクトルであり、βｔ，ｍｉｎは＾Ｖｔの最小の固有値であり、Ｍ´ｔは＾Ｖｔの最小の固有値に対応する固有ベクトルを長さが１となるように正規化したベクトルであり、Ｎは定数である。
【請求項３】
システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて予め定められた状態空間モデルで表わされるカルマンフィルタに従って、入力された前記観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、前記状態推定誤差を算出するカルマンフィルタ手段と、
前記状態推定手段によって算出された前記状態推定値に基づく前記観測値の残差の二乗の期待値を算出し、前記算出された前記期待値が予め定められた値以上であるときに、前記算出された前記期待値に基づいて、前記状態推定手段によって算出された前記状態推定値及び前記状態推定誤差の各々を補正するロバストエンハンスメント手段と、
前記時系列データの観測値の各々について、前記カルマンフィルタ手段による算出及び前記ロバストエンハンスメント手段による補正を繰り返したときに得られる前記状態推定値の各々を出力する出力手段と、
を含む、状態推定装置。
【請求項４】
コンピュータを、
システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて予め定められた状態空間モデルで表わされるカルマンフィルタに従って、入力された前記観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、前記状態推定誤差を算出するカルマンフィルタ手段、
前記状態推定手段によって算出された前記状態推定値に基づく前記観測値の残差の二乗の期待値を算出し、前記算出された前記期待値が予め定められた値以上であるときに、前記算出された前記期待値に基づいて、前記状態推定手段によって算出された前記状態推定値及び前記状態推定誤差の各々を補正するロバストエンハンスメント手段、及び
前記時系列データの観測値の各々について、前記カルマンフィルタ手段による算出及び前記ロバストエンハンスメント手段による補正を繰り返したときに得られる前記状態推定値の各々を出力する出力手段、
として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、状態推定装置、方法、及びプログラムに係り、特に、カルマンフィルタを用いて推定対象の状態を観測値から推定する状態推定装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
従来技術として、カルマンフィルタのロバスト性を確保するため、システムノイズ及び状態の推定誤差共分散行列を人工的に拡大させる方法がある（非特許文献１）。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】D. Simon, “Optimal State Estimation: Kalman, H∞, and Nonlinear Approaches”, John Wiley & Sons,Inc., Hoboken, New Jersey, 2006.
【発明の概要】
【発明が解決しようとする課題】
【０００４】
非特許文献１の方法では、いつ、どれ位、システムノイズ及び状態の推定誤差共分散行列を拡大させれば良いかが解析者の経験に基づいて決定されることが多く、フィルタの初心者にとってはパラメータ設定の敷居が高く扱いづらいという問題がある。
【０００５】
本発明では、上記問題点を解決するために成されたものであり、精度よく推定対象の状態を推定することができる状態推定装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００６】
上記目的を達成するために、第１の発明に係る状態推定方法は、カルマンフィルタ手段と、ロバストエンハンスメント手段と、出力手段とを含み、入力された推定対象の観測値の時系列データに基づいて、前記推定対象の状態を推定する状態推定装置における状態推定方法において、前記カルマンフィルタ手段は、システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて予め定められた状態空間モデルで表わされるカルマンフィルタに従って、入力された前記観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、前記状態推定誤差を算出し、前記ロバストエンハンスメント手段は、前記状態推定手段によって算出された前記状態推定値に基づく前記観測値の残差の二乗の期待値を算出し、前記算出された前記期待値が予め定められた値以上であるときに、前記算出された前記期待値に基づいて、前記状態推定手段によって算出された前記状態推定値及び前記状態推定誤差の各々を補正し、前記出力手段は、前記時系列データの観測値の各々について、前記カルマンフィルタ手段による算出及び前記ロバストエンハンスメント手段による補正を繰り返したときに得られる前記状態推定値の各々を出力する。
【０００７】
第２の発明に係る状態推定装置は、システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて予め定められた状態空間モデルで表わされるカルマンフィルタに従って、入力された前記観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、前記状態推定誤差を算出するカルマンフィルタ手段と、前記状態推定手段によって算出された前記状態推定値に基づく前記観測値の残差の二乗の期待値を算出し、前記算出された前記期待値が予め定められた値以上であるときに、前記算出された前記期待値に基づいて、前記状態推定手段によって算出された前記状態推定値及び前記状態推定誤差の各々を補正するロバストエンハンスメント手段と、前記時系列データの観測値の各々について、前記カルマンフィルタ手段による算出及び前記ロバストエンハンスメント手段による補正を繰り返したときに得られる前記状態推定値の各々を出力する出力手段と、を含んで構成されている。
【０００８】
第１の発明及び第２の発明によれば、カルマンフィルタ手段により、予め定められた状態空間モデルであらわされるカルマンフィルタに従って、入力された観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、状態推定誤差を算出し、ロバストエンハンスメント手段により、状態推定値に基づく観測地の残差の二乗の期待値を算出し、算出された期待値が予め定められた値以上であるときに、算出された期待値に基づいて、状態推定値及び状態推定誤差の各々を補正する。
【０００９】
そして、出力手段により、カルマンフィルタ手段による算出及びロバストエンハンスメント手段による補正を繰り返したときに得られる状態推定値の各々を出力する。
【００１０】
このように、第１の発明及び第２の発明によれば、精度よく推定対象の状態を推定することができる。
【００１１】
また、本発明のプログラムは、コンピュータを、システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて予め定められた状態空間モデルで表わされるカルマンフィルタに従って、入力された前記観測値と前回算出された状態推定値及び状態推定誤差とに基づいて、状態推定値を算出すると共に、前記状態推定誤差を算出するカルマンフィルタ手段、前記状態推定手段によって算出された前記状態推定値に基づく前記観測値の残差の二乗の期待値を算出し、前記算出された前記期待値が予め定められた値以上であるときに、前記算出された前記期待値に基づいて、前記状態推定手段によって算出された前記状態推定値及び前記状態推定誤差の各々を補正するロバストエンハンスメント手段、及び前記時系列データの観測値の各々について、前記カルマンフィルタ手段による算出及び前記ロバストエンハンスメント手段による補正を繰り返したときに得られる前記状態推定値の各々を出力する出力手段、として機能させるためのプログラムである。
【発明の効果】
【００１２】
以上説明したように、本発明の状態推定装置、方法、及びプログラムによれば、精度よく推定対象の状態を推定することができる。
【図面の簡単な説明】
【００１３】
【図１】観測データと状態推定結果のデータ構造の例を示す図である。
【図２】本発明の実施の形態に係る状態推定装置の機能的構成を示すブロック図である。
【図３】本発明の実施の形態に係る状態推定装置における状態推定処理ルーチンを示すフローチャート図である。
【図４】カルマンフィルタの結果とロバストカルマンフィルタの推定結果の例を示す図である。
【発明を実施するための形態】
【００１４】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１５】
＜発明の概要＞
本実施の形態では、入力されるデータは、観測機器によって観測された推定対象（現象）の観測データと、推定対象の背後にある支配方程式を解析者がモデル化した状態空間モデルである。そして、入力された観測データと状態空間モデルとに基づいて、カルマンフィルタ及びロバストエンハンスメントからなるロバストカルマンフィルタによる処理を行い、推定対象の状態を推定する結果を出力する。
【００１６】
本実施の形態において、フィルタが異常を検知した際、カルマンフィルタを通常の二乗誤差規範に対して準最適化する。これは、最小化する損失関数に罰則項を追加することで実現され、結果としてフィルタのモデリング誤差に対するロバスト性が向上する。一方、フィルタが正常な区間においては、通常の二乗誤差規範に対して最適であるカルマンフィルタと等価なアルゴリズムを適用する。
【００１７】
観測データは時系列データとなっており、観測データの次元をＭとする。状態空間モデルは解析対象の背後にある支配方程式を解析者がモデル化したものである。具体的には、システムノイズを用いて状態推定値を時間更新するための状態更新式、及び観測ノイズを用いて状態推定値と観測値との関係を示す観測方程式を用いて定められ、下記（１）式及び（２）式に表される線形な状態空間モデルである。
【００１８】
【数１】
2014215822.tif 000003
【００１９】
ここで、確率変数Ｘｔを時刻ｔにおけるＬ次元の状態（推定対象）、確率変数Ｙｔを時刻ｔにおける観測値と定義すると、上記（１）式及び（２）式は、上記（１）式が状態Ｘｔの時間発展となり、上記（２）式が状態Ｘｔの観測機構を表す連立方程式となる。
【００２０】
上記（１）式及び（２）式のＦ、Ｈはそれぞれ状態遷移行列、観測出力行列と定義する。μｔ−１、ωｔはシステムノイズ、及び観測ノイズと定義し、それぞれ平均０、共分散行列Ｑ及びＲの正規分布に従うノイズとする。状態空間モデルにおけるＦ，Ｈ、Ｑ、Ｒの値は解析者により設定する。
【００２１】
フィルタの問題は、観測値Ｙｔ：１＝Ｙｔ，Ｙｔ−１，・・・，Ｙ１を使用した状態Ｘｔの推定であり、カルマンフィルタはその推定誤差の二乗を最小化する最適な推定値＾Ｘｔを出力する。つまり、カルマンフィルタは下記（３）式の損失関数を最小化する推定値を出力する。
【００２２】
【数２】
2014215822.tif 000004
【００２３】
しかし、カルマンフィルタには、フィルタの回数が増加すると推定誤差の二乗が小さくなり過ぎ、推定誤差の共分散行列が小さくなり過ぎるために状態の時間変化に対する素早い追従性を失うという問題がある。推定対象は一般的に複雑な時間変化を示すので、状態空間モデルで事前に想定されたモデルの誤差もそれに伴い時間変化するが、フィルタの回数が大きくなるとその誤差に対してフィルタが敏感になってしまい、状態の時間変化に適切に追従することができなくなる。
【００２４】
＾Ｖｔは、観測値Ｙｔ：１が与えられた後の状態の推定誤差共分散行列であり、上記（３）式の二乗誤差はこの値により正規化されている。本実施の形態におけるロバストカルマンフィルタは下記（４）式及び（５）式の損失関数を最小化するように設計される。
【００２５】
【数３】
2014215822.tif 000005
【００２６】
【数４】
2014215822.tif 000006
【００２７】
は、予測残差の二乗の期待値（観測値の残差の二乗の期待値）であり、上記（４）式及び（５）式の損失関数は、状態の推定誤差は小さいが、フィルタの予測残差が大きくなっているときに通常の二乗誤差項（第一項）に罰則を与える形になっている。
【００２８】
従って、フィルタが異常を検知した場合には、状態の二乗誤差を最小化することを止め、第二項の付加によって準最適な推定値を出力するようにする。準最適な推定値の推定誤差共分散行列は、最適な推定値の推定誤差共分散行列より大きくなるため、結果として次のフィルタリングにおけるフィルタのロバスト性が自動的に回復する。
【００２９】
上記（５）式の予測残差の二乗は、￣Ｖ［ｖｔ］（観測値Ｙｔ−１：１が与えられた後の残差の予測共分散行列）で正規化されており、第一項の二乗誤差の正規化と合わせてある。
【００３０】
上記（４）式の第一項、第二項共にカイ二乗分布に従う確率変数の期待値になっており、第一項の期待値は状態の次元のＬ、第二項の期待値は観測値の次元のＭである。従って、両項のバランスをとるため、кを第二項にかける。
【００３１】
上記（４）式及び（５）式の損失関数を最小化するロバストカルマンフィルタの状態推定アルゴリズムは下記（６）〜（１３）式に示すようになる。
【００３２】
＜カルマンフィルタの状態推定アルゴリズム＞
【数５】
2014215822.tif 000007
【００３３】
＜準最適化によるロバストエンハンスメントの状態推定アルゴリズム＞
【数６】
2014215822.tif 000008
【００３４】
なお、上記の予測残差の二乗の期待値が上記（１３）式のγtに対応する。
【００３５】
カルマンフィルタには、状態空間モデルで設定された上記（１）式及び（２）式と時刻ｔにおける観測データの観測データＹｔが入力される。この情報に基づいて、上記（６）〜（１０）式を用いて、状態推定値の初期値＾Ｘ´０及び状態推定誤差の初期値＾Ｖ´０又は、ロバストエンハンスメントにおいて前回算出された状態推定値＾Ｘ´ｔ及び状態推定誤差＾Ｖ´ｔを時間更新、及び観測更新させ、カルマンフィルタにおける状態推定値＾Ｘｔと、カルマンフィルタにおける状態推定誤差＾Ｖｔを算出する。なお、状態推定値及び状態推定誤差の初期値は予め適切な値を解析者が設定しておく。
【００３６】
また、ロバストエンハンスメントにはカルマンフィルタの状態推定結果である＾Ｘｔと、＾Ｖｔとが入力され、上記（１１）式〜（１３）式を用いて、準最適化によるロバストエンハンスメントで新たに計算される状態推定値＾Ｘ´ｔと、状態推定誤差＾Ｖ´ｔが状態推定結果として出力される。
【００３７】
具体的には、上記（１３）式においてγｔを算出し、γｔ≧τの場合には、カルマンフィルタの状態推定結果である＾Ｘｔと、＾Ｖｔとに基づいて、上記（１１）式及び（１２）式を用いて、状態推定値＾Ｘｔを補正した状態推定値＾Ｘ´ｔと、状態推定誤差＾Ｖｔを補正した状態推定誤差＾Ｖ´ｔとを求める。
【００３８】
また、γｔ≦＜τの場合には、状態推定値＾Ｘｔを＾Ｘ´ｔとし、状態推定誤差＾Ｖｔを状態推定誤差＾Ｖ´ｔとする。
【００３９】
ここで、
【００４０】
【数７】
2014215822.tif 000009
【００４１】
であり、ｖｉは時刻ｉにおける観測残差（観測値の予測誤差）のベクトルを表す。βｔ，ｍｉｎは、＾Ｖｔの最小の固有値であり、Ｍ´ｔ，ｍｉｎはと＾Ｖｔの最小の固有値に対応する固有ベクトルを長さが１となるように正規化したベクトルである。τ及びＮの標準的な設定値（推奨値）は観測値が１次元の場合、
【００４２】
【数８】
2014215822.tif 000010
【００４３】
、Ｎ＝５〜１０である。観測値が２次元以上の場合もカイ二乗分布の有意水準５％の点を採用する。
【００４４】
＾Ｘ´ｔの計算における±はどちらでもよいが、＋と−を交互に使うことで、ロバストカルマンフィルタの推定量のバイアス性を軽減する効果を得る事が出来る。
【００４５】
カルマンフィルタとロバストエンハンスメントは繰返しの処理関係になっており、時刻ｔの観測値Ｙｔを得る度にカルマンフィルタ→ロバストエンハンストメントを繰り返し、ロバストエンハンスメントは状態の推定結果を出力する。この繰返しの処理体系全体をロバストカルマンフィルタと定義する。図１に観測データと状態推定結果のデータ構造を示す。
【００４６】
＜状態推定装置の構成＞
次に、本発明の実施の形態に係る状態推定装置の構成ついて説明する。図２に示すように、本発明の実施の形態に係る状態推定装置１００は、ＣＰＵとＲＡＭと後述する状態推定処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。この状態推定装置１００は、機能的には図２に示すように入力部１０と、演算部２０と、出力部５０とを備えている。
【００４７】
入力部１０は、観測機器（図示省略）によって得られた、推定対象の次元がＭの時系列データである時刻ｔの観測データＹｔの入力を受け付ける。
【００４８】
入力部１０は、解析者により入力される上記（１）式及び（２）式からなる状態空間モデルを受け付ける。また、入力部１０は、状態遷移行列Ｆ、観測出力行列Ｈ、共分散行列Ｑ及びＲ、の値の入力を受け付ける。
【００４９】
演算部２０は、カルマンフィルタ部３２と、ロバストエンハンスメント部３４とを備えている。
【００５０】
カルマンフィルタ部３２は、入力部１０により受け付けた各時刻ｔにおける観測データＹｔと、状態空間モデルと、状態推定値の初期値＾Ｘ´０及び状態推定誤差の初期値＾Ｖ´０又は、ロバストエンハンスメント部３４において前回算出された状態推定値＾Ｘ´ｔ及び状態推定誤差＾Ｖ´ｔとに基づいて、カルマンフィルタにおける状態推定値＾Ｘ´ｔと、カルマンフィルタにおける状態推定誤差＾Ｖｔを算出する。
【００５１】
具体的には、上記（６）式と、（７）式と、（１０）式とにより，￣Ｘｔと、￣Ｖｔと、Ｋｔとを求めた後に、上記（８）式及び（９）式により、＾Ｘｔと、＾Ｖｔとを求め、ロバストエンハンスメント部３４に出力する。
【００５２】
ロバストエンハンスメント部３４は、カルマンフィルタ部３２で算出された￣Ｖｔ、￣Ｘｔ、及び観測データＹｔに基づいて、上記（１３）式に従って、γｔを算出し、γｔとτの値を比較する。γｔ≧τの場合には、カルマンフィルタ部３２から入力された状態推定値＾Ｘｔと、状態推定誤差＾Ｖｔと、算出したγｔとに基づいて、上記（１１）式、及び（１２）式に従って、補正後の状態推定値＾Ｘ´ｔと、状態推定誤差＾Ｖ´ｔを求めカルマンフィルタ部３２及び出力部５０に出力する。なお、βｔ，ｍｉｎ及びＭ´ｔは、＾Ｖｔの値から求められる。
【００５３】
また、γｔ＜τの場合には、カルマンフィルタ部３２から入力された状態推定値＾Ｘｔを補正せずに状態推定値＾Ｘ´ｔとし、状態推定誤差＾Ｖｔを補正せずに状態推定誤差＾Ｖ´ｔとしてカルマンフィルタ部３２及び出力部５０に出力する。
【００５４】
＜状態推定装置の作用＞
次に、本発明の実施の形態に係る状態推定装置１００の作用について説明する。図３に示す状態推定処理ルーチンを実行する前に、入力部１０により状態空間モデルを受け付ける。そして、入力部１０により、各時刻ｔの観測データＹｔを受け付けているときに、状態推定装置１００は、図３に示す状態推定処理ルーチンを実行する。
【００５５】
まず、ステップＳ１００では、入力部１０において受け付けた上記（１）式及び（２）式からなる状態空間モデルを読み込む。なお、同時に状態遷移行列Ｆ、観測出力行列Ｈ、共分散行列Ｑ及びＲの値も読み込む。
【００５６】
次に、ステップＳ１０２では、入力部１０において受け付けた現時刻ｔの観測データＹｔを読み込む。
【００５７】
次に、ステップＳ１０４では、ステップＳ１００において読み込んだ状態空間モデルと、ステップＳ１０２において読み込んだ現時刻ｔにおける観測データＹｔと、状態推定値の初期値＾Ｘ´０及び状態推定誤差の初期値＾Ｖ´０又は、ステップＳ１１０又はステップＳ１１２において前回算出された状態推定値＾Ｘ´ｔ及び状態推定誤差＾Ｖ´ｔとに基づいて、上記（６）式〜（１０）式に従って、カルマンフィルタにおける状態推定値＾Ｘｔ及び状態推定誤差＾Ｖｔを算出する。
【００５８】
次に、ステップＳ１０６では、ステップＳ１０４において、上記（６）式及び（７）式に従って算出された￣Ｖｔ、￣Ｘｔ、及びステップＳ１０２において取得した観測データＹｔに基づいて、上記（１３）式を用いてγｔを算出する。
【００５９】
次に、ステップＳ１０８では、ステップＳ１０６において算出したγｔの値がτの値以上か否かの判定を行う。γｔの値がτの値以上である場合には、ステップＳ１１０に移行し、γｔの値がτの値未満の場合には、ステップＳ１１２に移行する。
【００６０】
ステップＳ１１０では、ステップＳ１０４において算出されたカルマンフィルタにおける状態推定値＾Ｘｔ、状態推定誤差＾Ｖｔ、及びステップＳ１０６において算出されたγｔに基づいて、上記（１１）式及び（１２）式に従って状態推定値＾Ｘｔを＾Ｘ´ｔに補正し、状態推定誤差＾Ｖｔを＾Ｖ´ｔに補正する。
【００６１】
ステップＳ１１２では、ステップＳ１０４において算出されたカルマンフィルタにおける状態推定値＾Ｘｔを補正をせずに＾Ｘ´ｔとし、状態推定誤差＾Ｖｔを補正をせずに＾Ｖ´ｔとする。
【００６２】
次に、ステップＳ１１４では、ステップＳ１１０又はステップＳ１１２において求められた状態推定値＾Ｘ´ｔ及び状態推定誤差＾Ｖ´ｔを推定結果としてカルマンフィルタ部３２及び出力部５０に出力し、ステップＳ１０２に移行し、次の時刻ｔの観測データＹｔを読み込みステップＳ１０２〜ステップＳ１０８の処理を繰り返す。
【００６３】
＜実験例＞
本発明の状態推定装置におけるロバストカルマンフィルタの有効性を、以下（１４）式に示すＡＲ（２）過程の係数推定の実験で検証した。
【００６４】
【数９】
2014215822.tif 000011
【００６５】
上記（１４）式で示されるＡＲ（２）過程はｔ≦１０００の区間は定常で、１０００＜ｔ≦２０００の区間は独立過程になり、その後また元の定常過程に戻るといった構造の変化を有するものである。このＡＲ（２）過程を状態空間モデルで表すと下記（１５）式になる。
【００６６】
【数１０】
2014215822.tif 000012
【００６７】
問題はこの２次元の状態ベクトルＸｔのフィルタ推定であり、通常のカルマンフィルタの結果と共に、ロバストカルマンフィルタの推定結果を以下に説明する。
【００６８】
図４の左上が観測値の時系列で、右上が時変ＡＲ係数の真値をプロットしたものである。左下がカルマンフィルタの推定結果で、最初ｎ区間（ｔ≦１０００）の状態は正しく推定されているが、その後の係数の変化は追従できていないことがわかる。これは最初の１０００回のフィルタ後、状態の推定誤差共分散行列が極端に小さくなり、フィルタのロバスト性がほぼ無くなったためである。右下がロバストカルマンフィルタの結果で、ＡＲ係数の時間変化を素早く追従していることがわかる。推定値の準最適化が正しく動作し、必要な時にフィルタのロバスト性が適切にエンハンスされている。
【００６９】
上述した数値実験の結果もシステムモデルの時間変化に対する優れた追従性を見せていることから、本アルゴリズムの有効性が確認できた。
【００７０】
以上説明したように、本発明の実施の形態に係る状態推定装置によれば、システムモデルの時間変化に対する追従性を向上させて、精度よく推定対象の状態を推定することができる。
【００７１】
また、推定対象（現象）の状態を観測値から推定するフィルタリングの問題において、モデリング誤差に対してロバストなカルマンフィルタを提案することができる。
【００７２】
また、フィルタが異常を感知すると自動的にロバスト性を向上させる自己回復機能を持ち、状態の推定値を観測値に基づき素早く修正することができる。
【００７３】
また、必要な時にフィルタのロバスト化が適切に行われるため、解析者の熟練した経験を必要とせず、結果として安定した状態の推定結果を出力することができる。
【００７４】
なお、本発明は、上記の実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【００７５】
上記の実施の形態では、出力部５０において、状態推定値＾Ｘ´ｔ及び状態推定誤差＾Ｖ´ｔを推定結果として出力する場合について説明したが、これに限定されるものではない。例えば、出力部５０において、状態推定値＾Ｘ´ｔのみを出力するようにしてもよい。
【符号の説明】
【００７６】
１０   入力部
２０   演算部
２２   入力部
３２   カルマンフィルタ部
３４   ロバストエンハンスメント部
５０   出力部
１００  状態推定装置
【図１】
2014215822.tif 000014
【図２】
2014215822.tif 000015
【図３】
2014215822.tif 000016
【図４】
2014215822.tif 000017