(57)【特許請求の範囲】
【請求項１】
コンピュータの処理による、Unscentedカルマン・フィルタの処理方法であって、
【数１】
0005818608.tif 000025
を与えるステップであって、  ここで、xkは内部状態、ykは出力、vkはシステムを駆動するプロセス・ノイズ、nkは観測ノイズで、
【数２】
0005818608.tif 000026
は所与の非線形関数である、ステップと、
下記ステップ(a)及びステップ(b)を実行する非一様Unscentedカルマン・フィルタの処理ステップであって、
前記ステップ(a)は、 z変換の伝達関数の分子と分母が鏡像多項式の関係によって記述される実係数オールパス・フィルタによって、xkの位相間隔を非一様化するステップであり及び、
前記ステップ(b)は、 前記非一様化されたxkを
【数３】
0005818608.tif 000027
に適用する計算ステップである、処理ステップと、
前記実係数オールパス・フィルタにおけるパラメータの符号を反転させた実係数オールパス・フィルタによって、前記非一様Unscentedカルマン・フィルタで計算されたxkの推定値の位相間隔を元に戻すステップを有する、
Unscentedカルマン・フィルタの処理方法。
【請求項２】
前記実係数オールパス・フィルタが、１次のオールパス・フィルタである、請求項１に記載の処理方法。
【請求項３】
前記実係数オールパス・フィルタが、２次のオールパス・フィルタである、請求項１に記載の処理方法。
【請求項４】
プラントの観測信号について周波数分析を行い、高周波成分と低周波成分の値に基づきワープ・パラメータを計算するステップと、該計算されたワープ・パラメータを前記実係数オールパス・フィルタにセットするステップをさらに有する、請求項１に記載の処理方法。
【請求項５】
コンピュータの処理による、Unscentedカルマン・フィルタの処理プログラムであって、
前記コンピュータに、
【数４】
0005818608.tif 000028
を与えるステップであって、  ここで、xkは内部状態、ykは出力、vkはシステムを駆動するプロセス・ノイズ、nkは観測ノイズで、
【数５】
0005818608.tif 000029
は所与の非線形関数である、ステップと、
下記ステップ(a)及びステップ(b)を実行する非一様Unscentedカルマン・フィルタの処理ステップであって、
前記ステップ(a)は、z変換の伝達関数の分子と分母が鏡像多項式の関係によって記述される実係数オールパス・フィルタによって、xkの位相間隔を非一様化するステップであり及び、
前記ステップ(b)は、 (b) 前記非一様化されたxkを
【数６】
0005818608.tif 000030
に適用する計算ステップである、処理ステップと、
前記実係数オールパス・フィルタにおけるパラメータの符号を反転させた実係数オールパス・フィルタによって、前記非一様Unscentedカルマン・フィルタで計算されたxkの推定値の位相間隔を元に戻すステップを実行させる、
Unscentedカルマン・フィルタの処理プログラム。
【請求項６】
前記実係数オールパス・フィルタが、１次のオールパス・フィルタである、請求項５に記載の処理プログラム。
【請求項７】
前記実係数オールパス・フィルタが、２次のオールパス・フィルタである、請求項５に記載の処理プログラム。
【請求項８】
プラントの観測信号について周波数分析を行い、高周波成分と低周波成分の値に基づきワープ・パラメータを計算するステップと、該計算されたワープ・パラメータを前記実係数オールパス・フィルタにセットするステップをさらに有する、請求項５に記載の処理プログラム。
【請求項９】
コンピュータの処理により動作する制御システムであって、
Unscentedカルマン・フィルタと、
プラントの出力信号を、位相間隔を非一様化して前記Unscentedカルマン・フィルタに入力する手段と、
前記Unscentedカルマン・フィルタの推定値の位相間隔を元に戻す手段を有する、
制御システム。
【請求項１０】
前記位相間隔を非一様化して前記Unscentedカルマン・フィルタに入力する手段と、前記Unscentedカルマン・フィルタの推定値の位相間隔を元に戻す手段が、z変換の伝達関数の分子と分母が鏡像多項式の関係によって記述される実係数オールパス・フィルタである、請求項９に記載の制御システム。
【請求項１１】
前記実係数オールパス・フィルタが、１次のオールパス・フィルタである、請求項１０に記載の制御システム。
【請求項１２】
前記実係数オールパス・フィルタが、２次のオールパス・フィルタである、請求項１０に記載の制御システム。
【請求項１３】
プラントの観測信号について周波数分析を行い、高周波成分と低周波成分の値に基づきワープ・パラメータを計算する手段と、該計算されたワープ・パラメータを前記実係数オールパス・フィルタにセットする手段をさらに有する、請求項１０に記載の制御システム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、カルマン・フィルタの処理技法に関し、より詳細には、カルマン・フィルタの推定誤差の改善技法に関するものである。
【背景技術】
【０００２】
従来より、エンジン、衛星の姿勢制御、リチウム・イオン・バッテリ、燃料電池のように、全体としては複雑であっても、個別のサブシステムの物理・化学法則がある程度分かっているプラントに対する制御器を設計するためのモデルとして、非線形状態空間モデルが用いられている。
【０００３】
このような非線形状態空間モデルの内部状態量を非線形関数の線形近似で推定するために、拡張カルマン・フィルタ(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を実行するモジュール（図２のブロック２０４に対応）、オールパス・フィルタとして動作するモジュール（図２のブロック２０２、２０６に対応）、プラント１１２に対するデータの入出力をコントロールするモジュール（図示しない）、全体の処理を制御するモジュール（図示しない）などが保存されており、これらのモジュールは、システムの起動時、あるいは必要時に主記憶１０６にロードされて実行される。
【００２１】
図２は、本発明の処理を実行する機能ブロック図を示す図である。機能ブロック２０４は、下記の式で記述される処理に従い、出力推定処理を行う。
【数１】
0005818608.tif 000002
ここで、xkは内部状態、ykは出力、vkはシステムを駆動するプロセス・ノイズ、nkは観測ノイズで、
【数２】
0005818608.tif 000003
は、プラント１１２の物理的あるいは化学的性質等を考慮して制御器を構成するための、所与の非線形関数である。
【００２２】
本発明を特徴づけるのは、オールパス・フィルタ２０２とオールパス・フィルタ２０６の存在である。
【００２３】
オールパス・フィルタ２０２は、そのz変換表示が分子と分母の多項式の関係が鏡像多項式の関係によって記述されるものであり、プラント１１２からの観測出力信号の位相間隔を非一様化してUKFブロック２０４に入力する。
【００２４】
オールパス・フィルタ２０６は、オールパス・フィルタ２０２と実質的に同一でよく、但し、オールパス・フィルタ２０２と逆の動作を行うようにパラメータをセットされ、すなわち、UKFブロック２０４からの、位相を非一様化された出力の位相を元に戻すように働く。
【００２５】
次に、図３のフローチャートを参照して、図１に示す、この実施例のシステムが実行する処理について説明する。
【００２６】
ステップ３０２では、システムは、インターフェース１１０を介して、プラント１１２の特性把握のための入出力信号の計測を行う。
【００２７】
ステップ３０４でシステムは、計測結果を用いて周波数分析を行い、高周波成分Hpと、低周波成分Lpの値を得る。このとき、システム同定の区間を決め、その観測信号（ノイズを含んでいてもかまわない）について周波数分析を行う。
【００２８】
ステップ３０６でシステムは、下記の式でワープ・パラメータλを計算する。ここでcは、|λ| < 1となるように、適当に決めた正の定数である。
【数３】
0005818608.tif 000004
このとき、分析結果に応じて、位相の非一様化を高周波成分中心に行うか、低周波成分中心に行うかを、ワープ・パラメータλとして反映する。例えば、ＦＦＴによるパワースペクトルを求め、規格化角周波数で直流成分を除くπ/2以下の成分をLpとし、π/2より大きくπ以下の成分総和をHpとして、上記の式を用いる。周波数分析は観測信号の傾向を知る目的で行うので、ウェーブレット解析を用いてもよい。
【００２９】
ステップ３０８でシステムは、インターフェース１１０を介して、プラント１１２の入出力信号の計測を行う。
【００３０】
ステップ３１０でシステムは、変数k = 0とセットする。ステップ３１２ではシステムは、サンプルが最後かどうか判断し、そうなら処理を終了する。そうでないならステップ３１４で、k+1番目のサンプル信号を用意し、kをk+1にインクリメントする。
【００３１】
ステップ３１６でシステムは、用意されたサンプル信号に対して、下記のz変換の式で記述されるオールパス・フィルタリングを適用して、位相の非一様化を行う。
【数４】
0005818608.tif 000005
【００３２】
上記のオールパス・フィルタリングの式は、z変換表示が分子と分母の多項式の関係が鏡像多項式の関係によって記述されるものであり、特に上記の例は１次の場合を示す。
【００３３】
より一般にM次の因果的実係数オールパス・フィルタは、次のように記述される。
【数５】
0005818608.tif 000006
ここで、
【数６】
0005818608.tif 000007
この式で、d1やdMはフィルタ係数である。
【００３４】
また、２次の例は下記のとおりである。
【数７】
0005818608.tif 000008
この実施例では、１次の式を用いるものとする。１次のオールパス・フィルタの場合、λが正（負）であれば高（低)周波数成分が多い箇所で解像度を増やす効果が得られる。
【００３５】
より具体的には、下記の式により観測信号の位相の非一様化を行う。
【数８】
0005818608.tif 000009
ここで、係数
【数９】
0005818608.tif 000010
の定義は、ステップ３２４に関連して後述する。
【００３６】
尚、以下では、上添字や下添字に含まれるwは非一様サンプリングによってワープ（非一様化）した領域の変数をあらわし、時間に関する項の下添字wも非一様サンプリングによってワープした領域での更新タイミングをあらわすこととする。また、現時刻kで必要な状態量や共分散行列等を時刻k-1までの観測で推定したものを添字"k|k-1"であらわすが、これを変数の肩に(-)を付けて代用表現することもある。また、肩のawは、状態ベクトルxに対して、プロセスノイズと観測ノイズをあわせたベクトルを拡張状態ベクトルとして再定義したものを意味する。
【００３７】
ステップ３１８でシステムは、下記の式により、非一様位相UKFの初期値を設定する。
【数１０】
0005818608.tif 000011
【数１１】
0005818608.tif 000012
【数１２】
0005818608.tif 000013
ここで、
【数１３】
0005818608.tif 000014
また、
【数１４】
0005818608.tif 000015
は非一様サンプリング上でのプロセスノイズ共分散行列、
【数１５】
0005818608.tif 000016
は非一様サンプリング上での観測ノイズ共分散行列である。
【００３８】
ステップ３２０でシステムは、下記の式により、各kw ∈ {1,...,∞}に対して、非一様位相UKFシグマポイント計算を行う。なお、従来のUKFでは、真値の周辺に妥当なシグマポイントをどのように与えるかという側面での提案がなされてきたが、本発明によれば、真値のダイナミックスに応じたシグマポイントを与えることができるという点でも特徴的である。すなわち、真値の変化が速いと想定される観測信号に対しては、更新速度を上げてシグマポイントのより多くのパターンを生成することに相当し、一方、真値の変化が緩やかであると想定される観測信号に対しては、過更新を避けるためのシグマポイントのパターン変化を抑えることに相当するオペレーションで、推定精度を向上させるものである。
【数１６】
0005818608.tif 000017
ここで、(P)iはPの第i列、
【数１７】
0005818608.tif 000018
で、ηはスケーリング・パラメータで、Lは下記で示す拡大状態ベクトルの次数である。
η = α2(L + μ) - Lで、1つの例ではα = 0.5, μ = 3 - Lである。
【数１８】
0005818608.tif 000019
【００３９】
ステップ３２２でシステムは、下記の式により、非一様位相UKF時間更新計算を行う。下記の式で、fとhは、拡大系において、
xk+1 = f(xk,νk,k), yk = h(xk,nk,k)の式で与えられているものである。
【数１９】
0005818608.tif 000020
【数２０】
0005818608.tif 000021
αとβは、任意の正の数である。
さらに、
【数２１】
0005818608.tif 000022
【００４０】
ステップ３２４でシステムは、下記の式により、非一様位相UKF観測更新計算を行う。
【数２２】
0005818608.tif 000023
【００４１】
ステップ３２６でシステムは、下記の式により、オールパス・フィルタリングで元の位相間隔信号量に戻す。
【数２３】
0005818608.tif 000024
【００４２】
システムはステップ３１４に戻り、サンプルが最後かどうか判断し、そうなら処理を終了する。そうでないならステップ３１４で、k+1番目のサンプル信号を用意し、kをk+1にインクリメントし、ステップ３１６に進む。
【００４３】
以上、特定の実施例に基づき本発明を説明してきたが、本発明はこの実施例に限定されず、拡大系UKFが適用可能な任意のプラントの制御系で、信号の推定精度を向上させることが可能である。
【符号の説明】
【００４４】
１１２    プラント
２０２    オールパス・フィルタ
２０４    Unscentedカルマン・フィルタ
２０６    オールパス・フィルタ
【図１】
0005818608.tif 000031
【図２】
0005818608.tif 000032
【図３】
0005818608.tif 000033