Help us understand the problem. What is going on with this article?

Ensemble Kalman Filterを用いた非線形システムの推定

More than 3 years have passed since last update.

非線形システムの推定

世の中の多くの動的な現象は以下のように表現されます.

\frac{d}{dt}x(t) = f(x, t), \quad x(0) = x_0

このような現象において, 一般にセンサで測るなどして得られる出力は$x$の一部であり, しかもセンサノイズ等が含まれる場合が多いです.
この場合, よく使われるオンライン推定手法としてカルマンフィルタが用いられますが, 一般のカルマンフィルタは線形システムにしか適用できず, 非線形システムに対しては非線形カルマンフィルタと呼ばれるものを適用する必要があります.

非線形Kalman Filter

非線形Kalman Filterには大きくわけて4種類(厳密には3種類?)あります.

拡張カルマンフィルタ(Extended Kalman Filter)

毎時刻, 推定値周りで線形近似を行うことで線形カルマンフィルタの手法を適用します. 線形近似を行っているため, 非線形性が強いとすぐ発散してしまいます.

無香料カルマンフィルタ(Unscented Kalman Filter)

状態空間上に分布する確率分布のうち, 代表的な点を選んできて, その点の状態遷移を求めることで予測し, 出力フィードバックにより更新するカルマンフィルタです. 拡張カルマンフィルタに比べて精度がよく, 状態ベクトルの次元が$n$の場合は代表点として$2n+1$個の点だけで済むので, 計算コスト的にも有利です.

パーティクルフィルタ(Particle Filter)

状態空間上にパーティクルを撒き, そのパーティクルの状態遷移を計算することで予測を行います. その後, 現実の出力をフィードバックして尤度を計算し, 尤度の低いパーティクルを尤度の高いパーティクルの位置に撒き直し(リサンプリング)をすることで推定を行います. 一般にはパーティクルの平均値を推定値として, 実際に値として用いる場合が多いです. 上記2つの手法は確率分布をガウス分布に近似しているので, この手法は非常に精度が高いことが利点ですが, パーティクルの数が非常に多く必要であり, 計算コストは非常に大きくなります.

アンサンブルカルマンフィルタ(Ensemble Kalman Filter)

手法はパーティクルフィルタと似ていて, 状態空間上にいくつかのサンプル状態を用意し, 状態アンサンブル(パーティクルの集まりと同じ意味)をそれぞれ状態遷移させて予測します. ここまではパーティクルフィルタと同じですが, 現実の出力をフィードバックした際に, 状態アンサンブルをより尤度の高い方へかき集めるような更新を行います. こちらもパーティクルフィルタと似たような利点と欠点を持ちますが, パーティクルフィルタほど精度は高くないものの, サンプル数は比較的少なめで済むので, 計算コスト的にはパーティクルフィルタよりも優れています. 特に大規模系の推定では威力を発揮し, 実際気象系では用いられることが多いです(というか, もともと気象系で使うために出てきたはずだったと思います).

Ensemble Kalman Filter

パーティクルフィルタ等は比較的文献を見かけますが, アンサンブルカルマンフィルタはほとんど見かけないので, ここではアンサンブルカルマンフィルタを取り上げたいと思います.

推定対象

ここでは, 以下の離散時間非線形システムの推定を考えます.

x_{k+1} = f(x_k, k) + w_k \\
y_k = Cx_k + v_k

また, 状態空間上のサンプルを以下のように表現します.

x^{(i)}_{k+1} = f(x^{(i)}_k, k) + w^{(i)}_k\\
y^{(i)}_k = Cx^{(i)}_k + v^{(i)}_k

ただし, $i \in \{ 1, \ldots, M \}$とし, $M$が総サンプル数となります.
以下では, $\bar{x}$を予測値, $\hat{x}$を更新値, $\tilde{x}$を誤差とします.

更新ステップ

1.予測平均値との誤差アンサンブル行列を作成
予測値の平均値との誤差を以下のように求めます.

\tilde{x}^{(i)}_k = \bar{x}^{(i)}_k - \frac{1}{M}\sum_{i=1}^{M}\bar{x}^{(i)}_k

このとき, 誤差アンサンブル行列を以下のように作ります.

\tilde{X}_k = \begin{bmatrix}\tilde{x}^{(1)}_k & \cdots & \tilde{x}^{(M)}_k \end{bmatrix}

2.予測出力も誤差アンサンブル行列を作成
予測出力を以下のように計算します.

\bar{y}^{(i)}_k = C\bar{x}^{(i)}_k + v^{(i)}_k

予測出力も同様に誤差アンサンブル行列を作ります.

\tilde{y}^{(i)}_k = \bar{y}^{(i)}_k - \frac{1}{M}\sum_{i=1}^{M}\bar{y}^{(i)}_k \\
\tilde{Y}_k = \begin{bmatrix}\tilde{y}^{(1)}_k & \cdots & \tilde{y}^{(M)}_k \end{bmatrix}

3.共分散行列を計算
共分散行列をそれぞれ以下のように計算します.

U_k = \frac{1}{M-1} \tilde{X}_k \tilde{X}^{\sf T}_k\\
V_k = \frac{1}{M-1} \tilde{Y}_k \tilde{Y}^{\sf T}_k

4.カルマンゲインを計算
カルマンゲイン$H_k$を以下のように計算します.

H_k = U_k(V_k)^{-1}

ここで, カルマンゲインの計算に逆行列が入っていますが, $V_k$のサイズが(出力のサイズ)×(出力のサイズ)であることに注意すると, 比較的低コストな逆行列計算であることに注意しておきます.

5.カルマンゲインを用いて更新
先ほど計算したカルマンゲインを用いて予測値を更新します.

\hat{x}^{(i)}_k = \bar{x}^{(i)}_k + H_k\left(y_k - \bar{y}^{(i)}_k \right)

6.推定値をアンサンブル平均とする
推定値を利用する場合は, アンサンブルの平均値を用いて推定値とします.

\hat{x}_k = \frac{1}{M}\sum_{i=1}^{M}\hat{x}^{(i)}_k

予測ステップ

状態の予測値は, システムの状態遷移関数を用いて以下のように計算します.

\bar{x}^{(i)}_{k+1} = f(\hat{x}^{(i)}_k, k) + w^{(i)}_k

数値シミュレーション

ここでは, Lorenz方程式と呼ばれる非線形ダイナミクスを利用して推定シミュレーションを行います.

Lorenz方程式

Lorenz方程式は以下の様な形で表されます.

\dot{x} = \sigma(y-x)\\
\dot{y} = rx -y-z\\
\dot{z} = -bz + xy

ただし, $\sigma = 10, r = 28, b = \frac{8}{3}$とします.

Lorenz方程式を表す関数は以下のようになります.

lorenz.m
function dx = lorenz(x, t)
    s = 10;
    R = 28;
    b= 8/3;

    dx = [-s*(x(1)-x(2))  R*x(1)-x(2)-x(1)*x(3)  x(1)*x(2) - b*x(3)];   
end

なお, ここでは出力値も好みに変更できるよう関数化しました.

myoutput.m
function y = myoutput(x,t)
    y = x';
end

Ensemble Kalman Filterのシミュレーション

サンプル数を100として, 以下のスクリプトによりシミュレーションを行いました.
今回は4次のルンゲクッタ法により離散化を行っています.

enkfsim.m
clear
close
clc

%% 分散の初期設定
Rtr = 0.3;
R = diag([0.5, 0.2, 0.1]);
Q = diag([2, 4, 3]);
dt = 0.01;

%% 真値の計算
t = [0:dt:30]; %時間
xtr0(1,:) = [5 3 3]; %真値の初期値
ytr(:,1) = myoutput(xtr0(1,:),t(1)); %真値の出力
w(:,1) =  Q*randn(length(xtr0(1,:)),1);
xtr(1,:) = xtr0;
for i = 2:length(t)
    k1 = lorenz(xtr(i-1,:), t(i-1));
    k2 = lorenz(xtr(i-1,:)+dt*k1/2 ,t(i-1)+dt/2);
    k3 = lorenz(xtr(i-1,:)+dt*k2/2 ,t(i-1)+dt/2);
    k4 = lorenz(xtr(i-1,:)+dt*k3 ,t(i-1)+dt);
    w(:,i) =  Q*randn(length(xtr(1,:)),1);
    xtr(i,:) = xtr(i-1,:) + dt* (k1+2*k2+2*k3+k4)/6 + w(:,i-1)'*dt;
    ytr(:,i) = myoutput(xtr(i,:),t(i));
end

%% 観測値の計算
[len_y, len_x] = size(ytr);
for i = 1:length(t)
    yt(:,i) = ytr(:,i) + R*randn(len_y,1);
end

%%  EnKFの初期設定
M = 100; %サンプルの数
barx0 = [2;10;1]; %状態の初期値
P0 = diag([1,3, 3]); %共分散行列の初期値
for j = 1:M
    x0(:,j) = barx0 + Q*randn(length(xtr(1,:)),1); %各サンプルの初期値
end
for j = 1:M
    y0(:,j) = myoutput(x0(:,j)',t(1)) + R*randn(len_y,1); %各サンプルの予測値計算
end

%% Ensemble Kalman Filterの計算
for i = 1:length(t)-1
    %%%%  EnKF: 観測更新ステップ  %%%%
    xhat = sum(x0(:,:)')'/M; %アンサンブル平均の計算
    yhat = sum(y0(:,:)')'/M; %予測値の平均値
    yest(:,i) = yhat;
    Pxy = ((kron(eye(1,100),xhat)) - x0)*((kron(eye(1,100),yhat)) - y0)'/(M-1);
    Pyy = ((kron(eye(1,100),yhat))-y0)*((kron(eye(1,100),yhat))-y0)'/(M-1);
    K = Pxy*inv(Pyy); %カルマンゲインの計算
    %%%%  EnKF: 予測更新ステップ  %%%%
    for j = 1:M
        x(:,j) = x0(:,j) + K*(yt(:, i) - y0(:,j));
        x0(:,j) = x(:,j) + lorenz(x(:,j),t(i))'.*dt + Q*randn(length(xtr(1,:)),1).*dt;
        y0(:,j) = myoutput(x0(:,j)',t(i+1)) + R*randn(len_y,1);
    end
end

%% 以下, 表示に関する処理
yest(:,length(t)) = sum(y0(:,:)')'/M;

figure(1);
subplot(2,2,1);
graph1 = plot(t,ytr(1,:),'r', t, yest(1,:));
hold on
grid on
legend('True', 'Estimate');
title('Ensemble Kalman Filter Simulation');
set(graph1, 'linewidth', 1.5);
xlabel('Time [s]');
ylabel('x');

subplot(2,2,2);
graph2 = plot(t,ytr(2,:),'r', t, yest(2,:));
hold on
grid on
legend('True', 'Estimate');
title('Ensemble Kalman Filter Simulation');
set(graph2, 'linewidth', 1.5);
xlabel('Time [s]');
ylabel('y');

subplot(2,2,3);
graph3 = plot(t,ytr(3,:),'r', t, yest(3,:));
hold on
grid on
legend('True', 'Estimate');
title('Ensemble Kalman Filter Simulation');
set(graph3, 'linewidth', 1.5);
xlabel('Time [s]');
ylabel('z');

subplot(2,2,4);
gg = plot3(yest(1,:), yest(2,:), yest(3,:), 'b', ytr(1,:), ytr(2,:), ytr(3,:),'r');
hold on
g2 = plot3(yest(1,1), yest(2,1), yest(3,1), 'b.', ytr(1,1), ytr(2,1), ytr(3,1),'r');
title('EnKF Simulation (Lorenz)');
grid on
set(gg, 'linewidth', 2);
set(g2, 'markersize', 10);
legend('Estimatin', 'True');

結果

enkf.png
青が推定値, 赤が真値です.
比較的いい感じに推定できてると思います.

結論

非線形カルマンフィルタのひとつであるEnsemble Kalman Filterの紹介と数値シミュレーションを行いました.
数値シミュレーションでは状態空間は3次元ですが, Ensemble Kalman Filterが威力を発揮するのは状態空間が高次元になったときです.
このような場合はParticle Filterでも十分な計算速度で高精度な推定を行えると思います.
しかし, Ensemble Kalman Filterでは100サンプルでも十分に推定できるところが利点と言えるでしょう.

g_rej55
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした