residual

拡張カルマン フィルターまたはアンセンテッド カルマン フィルターを使用するとき、測定残差と残差共分散を返す

ページ内をすべて折りたたむ

構文

[Residual,ResidualCovariance] = residual(obj,y)
[Residual,ResidualCovariance] = residual(obj,y,Um1,...,Umn)

説明

residual コマンドは、extendedKalmanFilter オブジェクトと unscentedKalmanFilter オブジェクトについて実際の測定値と予測された測定値の差を返します。残差を確認することは、フィルターの性能を検証する 1 つの方法です。残差は、"イノベーション" とも呼ばれ、予測誤差を定量化し、拡張カルマン フィルターおよびアンセンテッド カルマン フィルターの更新シーケンスの補正ステップを制御します。correct および predict (System Identification Toolbox) を使用して推定されたカルマン フィルターの状態を更新するときは、correct コマンドを使用する直前に residual コマンドを使用します。

[Residual,ResidualCovariance] = residual(obj,y) は、測定値 y とカルマン フィルター obj によって生成された予測測定値の間の残差 Residual を返します。この関数は、残差の共分散 ResidualCovariance も返します。

extendedKalmanFilter コマンドまたは unscentedKalmanFilter コマンドを使用して obj を作成します。非線形システムの状態遷移関数 f と測定関数 h を obj に指定します。オブジェクトの State プロパティには最新の推定状態値が格納されます。各タイム ステップで、correct および predict を一緒に使用して、状態 x を更新します。残差 s は、そのタイム ステップにおける実際の測定値と予測された測定値の差であり、s = y - h(x) と表現されます。残差の共分散 S は、R + RP の和です。ここで、R はフィルターの MeasurementNoise プロパティによって設定された測定ノイズの行列で、RP は測定空間に投影された状態共分散の行列です。

obj.MeasurementFcn で指定した測定関数 h に次のいずれかの形式が含まれている場合、この構文を使用します。

  • y(k) = h(x(k)) (加法性の測定ノイズ向け)

  • y(k) = h(x(k),v(k)) (非加法性の測定ノイズ向け)

ここで、y(k)、x(k)、および v(k) はタイム ステップ k におけるシステムの測定出力、状態、そして測定ノイズです。h への入力のみが状態および測定ノイズです。

例

システムの測定関数がこれらの入力を必要とする場合、[Residual,ResidualCovariance] = residual(obj,y,Um1,...,Umn) は追加の入力引数を指定します。複数の引数を指定できます。

測定関数 h に次のいずれかの形式がある場合は、この構文を使用します。

  • y(k) = h(x(k),Um1,...,Umn) (加法性の測定ノイズ向け)

  • y(k) = h(x(k),v(k),Um1,...,Umn) (非加法性の測定ノイズ向け)

例

例

すべて折りたたむ

ライブ スクリプトを開く

拡張カルマン フィルター アルゴリズムと測定出力データを使用して、ファン デル ポール振動子の状態を推定します。振動子には 2 つの状態と 1 つの出力があります。

振動子の拡張カルマン フィルター オブジェクトを作成します。前に記述して保存した状態遷移関数 vdpStateFcn.m と測定関数 vdpMeasurementFcn.m を使用します。これらの関数は、1 と等しい非線形パラメーター mu を使用して、ファン デル ポール振動子への離散近似を記述します。これらの関数は、システム内の加法性プロセスと測定ノイズを仮定します。2 つの状態の初期状態の値を [1;0] と指定します。これは、時間 k-1 までのシステム出力に関する情報に基づいて、初期時間 k における状態を推定した値 xˆ[k|k-1] です。

obj = extendedKalmanFilter(@vdpStateFcn,@vdpMeasurementFcn,[1;0]);

振動子から測定出力データ y を読み込みます。この例では、説明のためにシミュレーションされた静的データを使用します。データは、vdp_data.mat ファイルに保存されています。

load vdp_data.mat y

振動子のプロセス ノイズ共分散と測定ノイズ共分散を指定します。

obj.ProcessNoise = 0.01;
obj.MeasurementNoise = 0.16;

配列を初期化して、推定の結果を取得します。

residBuf = [];
xcorBuf = [];
xpredBuf = [];

拡張カルマン フィルター アルゴリズムを実装し、correct コマンドと predict コマンドを使用して振動子の状態を推定します。最初に、時間 k における測定値を使用して xˆ[k|k-1] を修正し、xˆ[k|k] を取得します。次に、時間 k までの測定値を使用して推定されたタイム ステップ k における状態推定 xˆ[k|k] を使用して、次のタイム ステップにおける状態値 xˆ[k+1|k] を予測します。

リアルタイムのデータ測定をシミュレートするには、一度に 1 つのタイム ステップの測定データを使用します。予測値と実際の測定値の残差を計算して、フィルターの実行と収束がどの程度適切であるか評価します。残差の計算はオプションのステップです。residual を使用する場合は、correct コマンドの直前にこのコマンドを置きます。予測が測定値と一致した場合、残差はゼロです。

タイム ステップのリアルタイム コマンドを実行した後に、結果をバッファー処理して、実行が完了した後に結果をプロットできます。

for k = 1:length(y)
    [Residual,ResidualCovariance] = residual(obj,y(k));
    [CorrectedState,CorrectedStateCovariance] = correct(obj,y(k));
    [PredictedState,PredictedStateCovariance] = predict(obj);
    
     residBuf(k,:) = Residual;
     xcorBuf(k,:) = CorrectedState';
     xpredBuf(k,:) = PredictedState';
end

correct コマンドを使用すると、obj.State および obj.StateCovariance は、タイム ステップ k での修正済み状態値 CorrectedState と修正済み状態推定誤差の共分散値 CorrectedStateCovariance で更新されます。predict コマンドを使用すると、obj.State および obj.StateCovariance は、タイム ステップ k+1 での予測値 PredictedState と PredictedStateCovariance で更新されます。residual コマンドを使用する場合は、obj プロパティを変更しないでください。

この例では、初期状態値が xˆ[k|k-1] (つまり、時間 k-1 までのシステム出力に基づく初期時間 k での状態推定値) であったため、predict の前に correct を使用しました。初期状態の値が xˆ[k-1|k-1] (つまり、k-1 までの測定値に基づく前の時間 k-1 での値) である場合は、最初に predict コマンドを使用します。predict と correct の使用順序の詳細については、predict コマンドと correct コマンドの使用を参照してください。

修正後の値を使用して、推定状態をプロットします。

plot(xcorBuf(:,1),xcorBuf(:,2))
title("Estimated States")

Figure contains an axes object. The axes object with title Estimated States contains an object of type line.

実際の測定値、修正された推定測定値、および残差をプロットします。vdpMeasurementFcn の測定関数の場合、測定値が最初の状態になります。

M = [y,xcorBuf(:,1),residBuf];
plot(M)
grid on
title("Actual and Estimated Measurements, Residual")
legend("Measured","Estimated","Residual")

Figure contains an axes object. The axes object with title Actual and Estimated Measurements, Residual contains 3 objects of type line. These objects represent Measured, Estimated, Residual.

推定値は測定値に近接して追跡します。初期の過渡状態の後、残差は実行全体で比較的小さく維持されます。

ライブ スクリプトを開く

次の状態遷移方程式と測定方程式に従って状態 x と測定値 y が変化するような入力 u をもつ、非線形システムについて考えます。

x[k]=x[k-1]+u[k-1]+w[k-1]

y[k]=x[k]+2*u[k]+v[k]2

システムのプロセス ノイズ w は加法性で、測定ノイズ v は非加法性です。

システムの状態遷移関数と測定関数を作成します。追加入力 u を使用して関数を指定します。

f = @(x,u)(sqrt(x+u));
h = @(x,v,u)(x+2*u+v^2);

f と h は状態遷移関数と測定関数をそれぞれ保存する無名関数に対する関数ハンドルです。測定関数では、測定ノイズが非加法性であるため、v も入力として指定されます。追加入力 u の前に、v が入力として指定されることに注意してください。

指定した関数を使用して、非線形システムの状態を推定するために拡張カルマン フィルター オブジェクトを作成します。状態の初期値を 1、測定ノイズを非加法性として指定します。

obj = extendedKalmanFilter(f,h,1,'HasAdditiveMeasurementNoise',false);

測定ノイズ共分散を指定します。

obj.MeasurementNoise = 0.01;

これで、predict コマンドと correct コマンドを使用して、システムの状態を推定できます。u の値を predict と correct に渡すと、状態遷移関数と測定関数にそれぞれ渡されます。

タイム ステップ k で測定値 y[k]=0.8 と入力 u[k]=0.2 を使用して状態推定値を修正します。

correct(obj,0.8,0.2)

次のタイム ステップで u[k]=0.2 が与えられた場合の状態を予測します。

predict(obj,0.2)

予測値と測定値の誤差、つまり "残差" を取得します。

[Residual, ResidualCovariance] = residual(obj,0.8,0.2);

入力引数

すべて折りたたむ

次のいずれかのコマンドを使用して作成された拡張カルマン フィルターまたはアンセンテッド カルマン フィルター。

現在のタイム ステップで測定されたシステム出力。N 要素のベクトルとして指定します。ここで、N は測定の回数です。

システムの測定関数への追加の入力引数。任意の型の入力引数として指定します。測定関数 h は obj の MeasurementFcn プロパティまたは MeasurementLikelihoodFcn プロパティで指定されます。状態値と測定ノイズ値の他に、関数で入力引数が必要な場合、これらの入力を residual コマンド構文で指定します。residual コマンドは、これらの入力を測定関数または測定尤度関数に渡し、推定出力を計算します。複数の引数を指定できます。

たとえば、測定関数または測定尤度関数が状態 x に加えて、システム入力 u と現在の時間 k を使用して、推定システム出力 y を計算すると仮定します。したがって、Um1 項と Um2 項は u(k) と k になります。これらの入力の結果、次の推定出力が得られます。

y(k) = h(x(k),u(k),k)

タイム ステップ k でオンライン状態推定補正を実行する前に、これらの追加の入力を residual コマンド構文で指定します。

[Residual,ResidualCovariance] = residual(obj,y,u(k),k);

追加の入力引数を使用する方法の例については、追加入力を使用した状態遷移関数と測定関数の指定を参照してください。

出力引数

すべて折りたたむ

現在の測定値と予測された測定値の間の残差。以下として返されます。

  • 単出力システムの場合、スカラー

  • 多出力システムの場合、サイズ N のベクトル。ここで N は測定出力の数

残差共分散。N 行 N 列の行列として返されます。ここで、N は測定出力の数です。

バージョン履歴

R2019b で導入

参考

| | |

トピック