説明

カルマン・フィルタの処理方法、プログラム及びシステム

【課題】 Unscentedカルマン・フィルタにおいて、プラント特性を考慮した非一様更新タイミングに基づき、推定精度を改善する技法を提供すること。
【解決手段】 カルマン・フィルタへの入力を、分子と分母の多項式の関係が鏡像多項式の関係によって記述される、オールパスフィルタを通すことで非一様更新を実現し、ポリフェーズフィルタによるリサンプリングとは異なった演算量を抑える構成とする。次に、カルマン・フィルタによって得られた推定状態量はそのままでは位相が非一様間隔のものであるため、上記非一様化の逆変換となるフィルタリングを行う。

【発明の詳細な説明】
【技術分野】
【0001】
この発明は、カルマン・フィルタの処理技法に関し、より詳細には、カルマン・フィルタの推定誤差の改善技法に関するものである。
【背景技術】
【0002】
従来より、エンジン、衛星の姿勢制御、リチウム・イオン・バッテリ、燃料電池のように、全体としては複雑であっても、個別のサブシステムの物理・化学法則がある程度分かっているプラントに対する制御器を設計するためのモデルとして、非線形状態空間モデルが用いられている。
【0003】
このような非線形状態空間モデルの内部状態量を非線形関数の線形近似で推定するために、拡張カルマン・フィルタ(Extended Kalman Filter, EKF)に基づく構成が採用されてきた。
【0004】
しかし、モデル・ベース開発の浸透によって、プラントの物理・化学的記述が増えていき、そのことから、動作性能に無駄のない制御器開発への志向が強まるなか、非線形システムの同定や状態量推定を、線形近似を回避してプラントの性質をより一層反映した非線形システムの同定や状態量推定の一般的枠組みを提供すること、及び精度の向上により高い性能をもった制御システムの設計技法を提供することが課題となってきた。
【0005】
このような課題に対し、ヤコビ行列式やヘッセ行列式を計算できない不連続系や一次近似誤差が大きくなる場合でも非線形関数を近似することなくそのまま扱える枠組みとしてUnscentedカルマン・フィルタ(UKF)が提案されており、S.J.Julier, J.K.Uhlmann, “Unscented Filtering and Nonlinear Estimation,” Proc. IEEE Vol.92, No.3, pp.401-422, March 2004などの論文で、EKFよりも優れた推定を行えることが示されている。
【0006】
また、Y.Wu, et al., “Unscented Kalman Filtering for Additive Noise Case: Augmented versus Nonaugmented,” IEEE Signal Processing Letters, Vol.12, No.5, pp.357-360, May 2005は、拡大系UKFが、通常の状態空間モデル形式のものよりも精度が向上すると報告する。
【先行技術文献】
【非特許文献】
【0007】
【非特許文献1】S.J.Julier, J.K.Uhlmann, “Unscented Filtering and Nonlinear Estimation,” Proc. IEEE Vol.92, No.3, pp.401-422, March 2004
【非特許文献2】Y.Wu, et al., “Unscented Kalman Filtering for Additive Noise Case: Augmented versus Nonaugmented,” IEEE Signal Processing Letters, Vol.12, No.5, pp.357-360, May 2005
【発明の概要】
【発明が解決しようとする課題】
【0008】
従来の拡大系UKFは、プラントの特性を表す非線形関数をそのまま用いており、マクロな特性がわかる(周波数成分の相対的多寡等を把握可能)にも関わらず、フィルタの更新サンプリング時間は一様であることが暗黙的に仮定されており、非線形関数によって生じると考えられる時間的・周波数的解像度の部分的違いを考慮している事例は見受けられなかった。
【0009】
そこで、この発明の目的は、拡大系UKFにおいて、プラント特性を考慮した非一様更新タイミングに基づき、推定精度を改善する技法を提供することにある。
【課題を解決するための手段】
【0010】
この発明に従うシステムは先ず、好適には、システム同定対象の区間を定め、その観測信号(ノイズを含んでいてもかまわない)について周波数分析を行う。そして、分析結果に応じて、位相の非一様化を高周波数成分中心に行うか、低周波数成分中心に行うかをワープパラメータとして反映する。
【0011】
次にこの発明に従うシステムは、カルマン・フィルタへの入力を、分子と分母の多項式の関係が鏡像多項式の関係によって記述される、オールパス・フィルタを通すことで非一様更新を実現し、ポリフェーズフィルタによるリサンプリングとは異なった演算量を抑える構成とする。
【0012】
次にこの発明に従うシステムは、カルマン・フィルタによって得られた推定状態量はそのままでは位相が非一様間隔のものであるため、上記非一様化の逆変換となるフィルタリングを行う。
【発明の効果】
【0013】
本願発明者の実験によれば、従来の拡大系UKFでも高い精度を達成することが難しいとされた非線形モデルに、本発明に従う非一様更新タイミングによる拡大系UKFを適用することにより、平均二乗誤差及び平均誤差分布の両方で、従来の拡大系UKFよりも改善がみられた。すなわち、非線形モデル対する拡大系UKFで、推定精度を向上することが可能となる。
【図面の簡単な説明】
【0014】
【図1】この発明を実施するためのハードウェア構成の概要ブロック図である。
【図2】この発明を実施するための機能構成の概要ブロック図である。
【図3】この発明の処理のフローチャートを示す図である。
【発明を実施するための形態】
【0015】
以下、図面に基づき、この発明の実施例を説明する。特に断わらない限り、同一の参照番号は、図面を通して、同一の対象を指すものとする。尚、以下で説明するのは、本発明の一実施形態であり、この発明を、この実施例で説明する内容に限定する意図はないことを理解されたい。
【0016】
図1は、この発明を実施するためのハードウェア構成の概要ブロック図である。この実施例は、エンジン、人工衛星、リチウム・イオン・バッテリ、燃料電池などのプラントからの状態を信号として測定し、それらの測定信号から測定信号を推定するカルマン・フィルタの処理を実行するシステムである。
【0017】
特にこの発明は、非線形関数を近似することなくそのまま扱える枠組みとしてUnscentedカルマン・フィルタ(UKF)、より詳しくは、拡大系UKFを実行するシステムに関する。
【0018】
図1において、システム・バス102には、CPU104と、主記憶(RAM)106と、フラッシュ・メモリ108が接続されている。システム・バス102にはさらに、PCI、USBなどの所定のインターフェース・カード110を介して、プラント112及び各種センサ(図示しない)が接続されている。
【0019】
CPU104は、Intel(R) 386EX、PowerPC、ARM(R)アーキテクチャのプロセッサなどの組み込みシステムに適合する任意のプロセッサを用いることができる。
【0020】
フラッシュ・メモリ108には、組み込み用Linux(R)などのオペレーティング・システム、拡大系UKFを実行するモジュール(図2のブロック204に対応)、オールパス・フィルタとして動作するモジュール(図2のブロック202、206に対応)、プラント112に対するデータの入出力をコントロールするモジュール(図示しない)、全体の処理を制御するモジュール(図示しない)などが保存されており、これらのモジュールは、システムの起動時、あるいは必要時に主記憶106にロードされて実行される。
【0021】
図2は、本発明の処理を実行する機能ブロック図を示す図である。機能ブロック204は、下記の式で記述される処理に従い、出力推定処理を行う。
【数1】


ここで、xkは内部状態、ykは出力、vkはシステムを駆動するプロセス・ノイズ、nkは観測ノイズで、
【数2】


は、プラント112の物理的あるいは化学的性質等を考慮して制御器を構成するための、所与の非線形関数である。
【0022】
本発明を特徴づけるのは、オールパス・フィルタ202とオールパス・フィルタ206の存在である。
【0023】
オールパス・フィルタ202は、そのz変換表示が分子と分母の多項式の関係が鏡像多項式の関係によって記述されるものであり、プラント112からの観測出力信号の位相間隔を非一様化してUKFブロック204に入力する。
【0024】
オールパス・フィルタ206は、オールパス・フィルタ202と実質的に同一でよく、但し、オールパス・フィルタ202と逆の動作を行うようにパラメータをセットされ、すなわち、UKFブロック204からの、位相を非一様化された出力の位相を元に戻すように働く。
【0025】
次に、図3のフローチャートを参照して、図1に示す、この実施例のシステムが実行する処理について説明する。
【0026】
ステップ302では、システムは、インターフェース110を介して、プラント112の特性把握のための入出力信号の計測を行う。
【0027】
ステップ304でシステムは、計測結果を用いて周波数分析を行い、高周波成分Hpと、低周波成分Lpの値を得る。このとき、システム同定の区間を決め、その観測信号(ノイズを含んでいてもかまわない)について周波数分析を行う。
【0028】
ステップ306でシステムは、下記の式でワープ・パラメータλを計算する。ここでcは、|λ| < 1となるように、適当に決めた正の定数である。
【数3】


このとき、分析結果に応じて、位相の非一様化を高周波成分中心に行うか、低周波成分中心に行うかを、ワープ・パラメータλとして反映する。例えば、FFTによるパワースペクトルを求め、規格化角周波数で直流成分を除くπ/2以下の成分をLpとし、π/2より大きくπ以下の成分総和をHpとして、上記の式を用いる。周波数分析は観測信号の傾向を知る目的で行うので、ウェーブレット解析を用いてもよい。
【0029】
ステップ308でシステムは、インターフェース110を介して、プラント112の入出力信号の計測を行う。
【0030】
ステップ310でシステムは、変数k = 0とセットする。ステップ312ではシステムは、サンプルが最後かどうか判断し、そうなら処理を終了する。そうでないならステップ314で、k+1番目のサンプル信号を用意し、kをk+1にインクリメントする。
【0031】
ステップ316でシステムは、用意されたサンプル信号に対して、下記のz変換の式で記述されるオールパス・フィルタリングを適用して、位相の非一様化を行う。
【数4】

【0032】
上記のオールパス・フィルタリングの式は、z変換表示が分子と分母の多項式の関係が鏡像多項式の関係によって記述されるものであり、特に上記の例は1次の場合を示す。
【0033】
より一般にM次の因果的実係数オールパス・フィルタは、次のように記述される。
【数5】


ここで、
【数6】


この式で、d1やdMはフィルタ係数である。
【0034】
また、2次の例は下記のとおりである。
【数7】


この実施例では、1次の式を用いるものとする。1次のオールパス・フィルタの場合、λが正(負)であれば高(低)周波数成分が多い箇所で解像度を増やす効果が得られる。
【0035】
より具体的には、下記の式により観測信号の位相の非一様化を行う。
【数8】


ここで、係数
【数9】


の定義は、ステップ324に関連して後述する。
【0036】
尚、以下では、上添字や下添字に含まれるwは非一様サンプリングによってワープ(非一様化)した領域の変数をあらわし、時間に関する項の下添字wも非一様サンプリングによってワープした領域での更新タイミングをあらわすこととする。また、現時刻kで必要な状態量や共分散行列等を時刻k-1までの観測で推定したものを添字"k|k-1"であらわすが、これを変数の肩に(-)を付けて代用表現することもある。また、肩のawは、状態ベクトルxに対して、プロセスノイズと観測ノイズをあわせたベクトルを拡張状態ベクトルとして再定義したものを意味する。
【0037】
ステップ318でシステムは、下記の式により、非一様位相UKFの初期値を設定する。
【数10】


【数11】


【数12】


ここで、
【数13】


また、
【数14】


は非一様サンプリング上でのプロセスノイズ共分散行列、
【数15】


は非一様サンプリング上での観測ノイズ共分散行列である。
【0038】
ステップ320でシステムは、下記の式により、各kw ∈ {1,...,∞}に対して、非一様位相UKFシグマポイント計算を行う。なお、従来のUKFでは、真値の周辺に妥当なシグマポイントをどのように与えるかという側面での提案がなされてきたが、本発明によれば、真値のダイナミックスに応じたシグマポイントを与えることができるという点でも特徴的である。すなわち、真値の変化が速いと想定される観測信号に対しては、更新速度を上げてシグマポイントのより多くのパターンを生成することに相当し、一方、真値の変化が緩やかであると想定される観測信号に対しては、過更新を避けるためのシグマポイントのパターン変化を抑えることに相当するオペレーションで、推定精度を向上させるものである。
【数16】


ここで、(P)iはPの第i列、
【数17】


で、ηはスケーリング・パラメータで、Lは下記で示す拡大状態ベクトルの次数である。
η = α2(L + μ) - Lで、1つの例ではα = 0.5, μ = 3 - Lである。
【数18】

【0039】
ステップ322でシステムは、下記の式により、非一様位相UKF時間更新計算を行う。下記の式で、fとhは、拡大系において、
xk+1 = f(xkk,k), yk = h(xk,nk,k)の式で与えられているものである。
【数19】


【数20】


αとβは、任意の正の数である。
さらに、
【数21】

【0040】
ステップ324でシステムは、下記の式により、非一様位相UKF観測更新計算を行う。
【数22】

【0041】
ステップ326でシステムは、下記の式により、オールパス・フィルタリングで元の位相間隔信号量に戻す。
【数23】

【0042】
システムはステップ314に戻り、サンプルが最後かどうか判断し、そうなら処理を終了する。そうでないならステップ314で、k+1番目のサンプル信号を用意し、kをk+1にインクリメントし、ステップ316に進む。
【0043】
以上、特定の実施例に基づき本発明を説明してきたが、本発明はこの実施例に限定されず、拡大系UKFが適用可能な任意のプラントの制御系で、信号の推定精度を向上させることが可能である。
【符号の説明】
【0044】
112 プラント
202 オールパス・フィルタ
204 Unscentedカルマン・フィルタ
206 オールパス・フィルタ

【特許請求の範囲】
【請求項1】
コンピュータの処理による、Unscentedカルマン・フィルタの処理方法であって、
【数1】

を与えるステップであって、 ここで、xkは内部状態、ykは出力、vkはシステムを駆動するプロセス・ノイズ、nkは観測ノイズで、
【数2】

は所与の非線形関数である、ステップと、
下記ステップ(a)及びステップ(b)を実行する非一様Unscentedカルマン・フィルタの処理ステップと、
(a) z変換の伝達関数の分子と分母が鏡像多項式の関係によって記述される実係数オールパス・フィルタによって、xkの位相間隔を非一様化するステップ及び、
(b) 前記非一様化されたxk
【数3】

に適用する計算ステップ、
前記オールパス・フィルタにおけるパラメータの符号を反転させた実係数オールパス・フィルタによって、前記非一様Unscentedカルマン・フィルタで計算されたxkの推定値の位相間隔を元に戻すステップを有する、
Unscentedカルマン・フィルタの処理方法。
【請求項2】
前記オールパス・フィルタが、1次のオールパス・フィルタである、請求項1に記載の方法。
【請求項3】
前記オールパス・フィルタが、2次のオールパス・フィルタである、請求項1に記載の方法。
【請求項4】
プラントの観測信号について周波数分析を行い、高周波成分と低周波成分の値に基づきワープ・パラメータを計算するステップと、該計算されたワープ・パラメータを前記オールパス・フィルタにセットするステップをさらに有する、請求項1に記載の方法。
【請求項5】
コンピュータの処理による、Unscentedカルマン・フィルタの処理プログラムであって、
前記コンピュータに、
【数4】

を与えるステップであって、 ここで、xkは内部状態、ykは出力、vkはシステムを駆動するプロセス・ノイズ、nkは観測ノイズで、
【数5】

は所与の非線形関数である、ステップと、
下記ステップ(a)及びステップ(b)を実行する非一様Unscentedカルマン・フィルタの処理ステップと、
(a) z変換の伝達関数の分子と分母が鏡像多項式の関係によって記述される実係数オールパス・フィルタによって、xkの位相間隔を非一様化するステップ及び、
(b) 前記非一様化されたxk
【数6】

に適用する計算ステップ、
前記オールパス・フィルタにおけるパラメータの符号を反転させた実係数オールパス・フィルタによって、前記非一様Unscentedカルマン・フィルタで計算されたxkの推定値の位相間隔を元に戻すステップを実行させる、
Unscentedカルマン・フィルタの処理プログラム。
【請求項6】
前記オールパス・フィルタが、1次のオールパス・フィルタである、請求項5に記載のプログラム。
【請求項7】
前記オールパス・フィルタが、2次のオールパス・フィルタである、請求項5に記載のプログラム。
【請求項8】
プラントの観測信号について周波数分析を行い、高周波成分と低周波成分の値に基づきワープ・パラメータを計算するステップと、該計算されたワープ・パラメータを前記オールパス・フィルタにセットするステップをさらに有する、請求項5に記載のプログラム。
【請求項9】
コンピュータの処理により動作する制御システムであって、
Unscentedカルマン・フィルタと、
プラントの出力信号を、位相間隔を非一様化して前記Unscentedカルマン・フィルタに入力する手段と、
前記Unscentedカルマン・フィルタの推定値の位相間隔を元に戻す手段を有する、
制御システム。
【請求項10】
前記位相間隔を非一様化して前記Unscentedカルマン・フィルタに入力する手段と、前記Unscentedカルマン・フィルタの推定値の位相間隔を元に戻す手段が、z変換の伝達関数の分子と分母が鏡像多項式の関係によって記述される実係数オールパス・フィルタである、請求項9に記載の制御システム。
【請求項11】
前記オールパス・フィルタが、1次のオールパス・フィルタである、請求項9に記載の制御システム。
【請求項12】
前記オールパス・フィルタが、2次のオールパス・フィルタである、請求項9に記載の制御システム。
【請求項13】
プラントの観測信号について周波数分析を行い、高周波成分と低周波成分の値に基づきワープ・パラメータを計算する手段と、該計算されたワープ・パラメータを前記オールパス・フィルタにセットする手段をさらに有する、請求項9に記載の制御システム。

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate


【公開番号】特開2013−74365(P2013−74365A)
【公開日】平成25年4月22日(2013.4.22)
【国際特許分類】
【出願番号】特願2011−210265(P2011−210265)
【出願日】平成23年9月27日(2011.9.27)
【出願人】(390009531)インターナショナル・ビジネス・マシーンズ・コーポレーション (4,084)
【氏名又は名称原語表記】INTERNATIONAL BUSINESS MACHINES CORPORATION
【Fターム(参考)】