JP5818608B2

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

Info

Publication number: JP5818608B2
Application number: JP2011210265A
Authority: JP
Inventors: 淳也清水
Original assignee: International Business Machines Corp
Current assignee: International Business Machines Corp
Priority date: 2011-09-27
Filing date: 2011-09-27
Publication date: 2015-11-18
Anticipated expiration: 2031-09-27
Also published as: US20130080492A1; JP2013074365A; US8903879B2; CN103023459B; CN103023459A

Description

この発明は、カルマン・フィルタの処理技法に関し、より詳細には、カルマン・フィルタの推定誤差の改善技法に関するものである。

従来より、エンジン、衛星の姿勢制御、リチウム・イオン・バッテリ、燃料電池のように、全体としては複雑であっても、個別のサブシステムの物理・化学法則がある程度分かっているプラントに対する制御器を設計するためのモデルとして、非線形状態空間モデルが用いられている。

このような非線形状態空間モデルの内部状態量を非線形関数の線形近似で推定するために、拡張カルマン・フィルタ(Extended Kalman Filter, EKF)に基づく構成が採用されてきた。

しかし、モデル・ベース開発の浸透によって、プラントの物理・化学的記述が増えていき、そのことから、動作性能に無駄のない制御器開発への志向が強まるなか、非線形システムの同定や状態量推定を、線形近似を回避してプラントの性質をより一層反映した非線形システムの同定や状態量推定の一般的枠組みを提供すること、及び精度の向上により高い性能をもった制御システムの設計技法を提供することが課題となってきた。

このような課題に対し、ヤコビ行列式やヘッセ行列式を計算できない不連続系や一次近似誤差が大きくなる場合でも非線形関数を近似することなくそのまま扱える枠組みとして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よりも優れた推定を行えることが示されている。

また、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が、通常の状態空間モデル形式のものよりも精度が向上すると報告する。

S.J.Julier, J.K.Uhlmann, "Unscented Filtering and Nonlinear Estimation," Proc. IEEE Vol.92, No.3, pp.401-422, March 2004 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は、プラントの特性を表す非線形関数をそのまま用いており、マクロな特性がわかる（周波数成分の相対的多寡等を把握可能）にも関わらず、フィルタの更新サンプリング時間は一様であることが暗黙的に仮定されており、非線形関数によって生じると考えられる時間的・周波数的解像度の部分的違いを考慮している事例は見受けられなかった。

そこで、この発明の目的は、拡大系UKFにおいて、プラント特性を考慮した非一様更新タイミングに基づき、推定精度を改善する技法を提供することにある。

この発明に従うシステムは先ず、好適には、システム同定対象の区間を定め、その観測信号（ノイズを含んでいてもかまわない）について周波数分析を行う。そして、分析結果に応じて、位相の非一様化を高周波数成分中心に行うか、低周波数成分中心に行うかをワープパラメータとして反映する。

次にこの発明に従うシステムは、カルマン・フィルタへの入力を、分子と分母の多項式の関係が鏡像多項式の関係によって記述される、オールパス・フィルタを通すことで非一様更新を実現し、ポリフェーズフィルタによるリサンプリングとは異なった演算量を抑える構成とする。

次にこの発明に従うシステムは、カルマン・フィルタによって得られた推定状態量はそのままでは位相が非一様間隔のものであるため、上記非一様化の逆変換となるフィルタリングを行う。

本願発明者の実験によれば、従来の拡大系UKFでも高い精度を達成することが難しいとされた非線形モデルに、本発明に従う非一様更新タイミングによる拡大系UKFを適用することにより、平均二乗誤差及び平均誤差分布の両方で、従来の拡大系UKFよりも改善がみられた。すなわち、非線形モデル対する拡大系UKFで、推定精度を向上することが可能となる。

この発明を実施するためのハードウェア構成の概要ブロック図である。この発明を実施するための機能構成の概要ブロック図である。この発明の処理のフローチャートを示す図である。

以下、図面に基づき、この発明の実施例を説明する。特に断わらない限り、同一の参照番号は、図面を通して、同一の対象を指すものとする。尚、以下で説明するのは、本発明の一実施形態であり、この発明を、この実施例で説明する内容に限定する意図はないことを理解されたい。

図１は、この発明を実施するためのハードウェア構成の概要ブロック図である。この実施例は、エンジン、人工衛星、リチウム・イオン・バッテリ、燃料電池などのプラントからの状態を信号として測定し、それらの測定信号から測定信号を推定するカルマン・フィルタの処理を実行するシステムである。

特にこの発明は、非線形関数を近似することなくそのまま扱える枠組みとしてUnscentedカルマン・フィルタ(UKF)、より詳しくは、拡大系UKFを実行するシステムに関する。

図１において、システム・バス１０２には、ＣＰＵ１０４と、主記憶（ＲＡＭ）１０６と、フラッシュ・メモリ１０８が接続されている。システム・バス１０２にはさらに、ＰＣＩ、ＵＳＢなどの所定のインターフェース・カード１１０を介して、プラント１１２及び各種センサ（図示しない）が接続されている。

ＣＰＵ１０４は、Intel(R) 386EX、PowerPC、ARM(R)アーキテクチャのプロセッサなどの組み込みシステムに適合する任意のプロセッサを用いることができる。

フラッシュ・メモリ１０８には、組み込み用Linux(R)などのオペレーティング・システム、拡大系UKFを実行するモジュール（図２のブロック２０４に対応）、オールパス・フィルタとして動作するモジュール（図２のブロック２０２、２０６に対応）、プラント１１２に対するデータの入出力をコントロールするモジュール（図示しない）、全体の処理を制御するモジュール（図示しない）などが保存されており、これらのモジュールは、システムの起動時、あるいは必要時に主記憶１０６にロードされて実行される。

図２は、本発明の処理を実行する機能ブロック図を示す図である。機能ブロック２０４は、下記の式で記述される処理に従い、出力推定処理を行う。

ここで、x_kは内部状態、y_kは出力、v_kはシステムを駆動するプロセス・ノイズ、n_kは観測ノイズで、

は、プラント１１２の物理的あるいは化学的性質等を考慮して制御器を構成するための、所与の非線形関数である。

本発明を特徴づけるのは、オールパス・フィルタ２０２とオールパス・フィルタ２０６の存在である。

オールパス・フィルタ２０２は、そのz変換表示が分子と分母の多項式の関係が鏡像多項式の関係によって記述されるものであり、プラント１１２からの観測出力信号の位相間隔を非一様化してUKFブロック２０４に入力する。

オールパス・フィルタ２０６は、オールパス・フィルタ２０２と実質的に同一でよく、但し、オールパス・フィルタ２０２と逆の動作を行うようにパラメータをセットされ、すなわち、UKFブロック２０４からの、位相を非一様化された出力の位相を元に戻すように働く。

次に、図３のフローチャートを参照して、図１に示す、この実施例のシステムが実行する処理について説明する。

ステップ３０２では、システムは、インターフェース１１０を介して、プラント１１２の特性把握のための入出力信号の計測を行う。

ステップ３０４でシステムは、計測結果を用いて周波数分析を行い、高周波成分Hpと、低周波成分Lpの値を得る。このとき、システム同定の区間を決め、その観測信号（ノイズを含んでいてもかまわない）について周波数分析を行う。

ステップ３０６でシステムは、下記の式でワープ・パラメータλを計算する。ここでcは、|λ| < 1となるように、適当に決めた正の定数である。

このとき、分析結果に応じて、位相の非一様化を高周波成分中心に行うか、低周波成分中心に行うかを、ワープ・パラメータλとして反映する。例えば、ＦＦＴによるパワースペクトルを求め、規格化角周波数で直流成分を除くπ/2以下の成分をLpとし、π/2より大きくπ以下の成分総和をHpとして、上記の式を用いる。周波数分析は観測信号の傾向を知る目的で行うので、ウェーブレット解析を用いてもよい。

ステップ３０８でシステムは、インターフェース１１０を介して、プラント１１２の入出力信号の計測を行う。

ステップ３１０でシステムは、変数k = 0とセットする。ステップ３１２ではシステムは、サンプルが最後かどうか判断し、そうなら処理を終了する。そうでないならステップ３１４で、k+1番目のサンプル信号を用意し、kをk+1にインクリメントする。

ステップ３１６でシステムは、用意されたサンプル信号に対して、下記のz変換の式で記述されるオールパス・フィルタリングを適用して、位相の非一様化を行う。

上記のオールパス・フィルタリングの式は、z変換表示が分子と分母の多項式の関係が鏡像多項式の関係によって記述されるものであり、特に上記の例は１次の場合を示す。

より一般にM次の因果的実係数オールパス・フィルタは、次のように記述される。

ここで、

この式で、d₁やd_Mはフィルタ係数である。

また、２次の例は下記のとおりである。

この実施例では、１次の式を用いるものとする。１次のオールパス・フィルタの場合、λが正（負）であれば高（低)周波数成分が多い箇所で解像度を増やす効果が得られる。

より具体的には、下記の式により観測信号の位相の非一様化を行う。

ここで、係数

の定義は、ステップ３２４に関連して後述する。

尚、以下では、上添字や下添字に含まれるwは非一様サンプリングによってワープ（非一様化）した領域の変数をあらわし、時間に関する項の下添字wも非一様サンプリングによってワープした領域での更新タイミングをあらわすこととする。また、現時刻kで必要な状態量や共分散行列等を時刻k-1までの観測で推定したものを添字"k|k-1"であらわすが、これを変数の肩に(-)を付けて代用表現することもある。また、肩のa_wは、状態ベクトルxに対して、プロセスノイズと観測ノイズをあわせたベクトルを拡張状態ベクトルとして再定義したものを意味する。

ステップ３１８でシステムは、下記の式により、非一様位相UKFの初期値を設定する。

ここで、

また、

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

は非一様サンプリング上での観測ノイズ共分散行列である。

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

ここで、(P)_iはPの第i列、

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

ステップ３２２でシステムは、下記の式により、非一様位相UKF時間更新計算を行う。下記の式で、fとhは、拡大系において、
x_k+1 = f(x_k,ν_k,k), y_k = h(x_k,n_k,k)の式で与えられているものである。

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

ステップ３２４でシステムは、下記の式により、非一様位相UKF観測更新計算を行う。

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

システムはステップ３１４に戻り、サンプルが最後かどうか判断し、そうなら処理を終了する。そうでないならステップ３１４で、k+1番目のサンプル信号を用意し、kをk+1にインクリメントし、ステップ３１６に進む。

以上、特定の実施例に基づき本発明を説明してきたが、本発明はこの実施例に限定されず、拡大系UKFが適用可能な任意のプラントの制御系で、信号の推定精度を向上させることが可能である。

１１２プラント
２０２オールパス・フィルタ
２０４ Unscentedカルマン・フィルタ
２０６オールパス・フィルタ

Claims

コンピュータの処理による、Unscentedカルマン・フィルタの処理方法であって、

を与えるステップであって、ここで、x_kは内部状態、y_kは出力、v_kはシステムを駆動するプロセス・ノイズ、n_kは観測ノイズで、

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

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

を与えるステップであって、ここで、x_kは内部状態、y_kは出力、v_kはシステムを駆動するプロセス・ノイズ、n_kは観測ノイズで、

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

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