Matlab:在通过卡尔曼滤波器进行状态估计后，如何模拟模型 [英] Matlab: How do I simulate the model after state estimation from Kalman filter

问题描述

我正在尝试为以下一维AR模型实现卡尔曼滤波器的基本方程式:

I am trying to implement the basic Equations for Kalman filter for the following 1 dimensional AR model:

x(t) = a_1x(t-1) + a_2x(t-2) + w(t)  

y(t) = Cx(t) + v(t);

KF状态空间模型:

x(t+1) = Ax(t) + w(t)

y(t) = Cx(t) + v(t)

w(t) = N(0,Q)

v(t) = N(0,R)

其中

 % A - state transition matrix
% C - observation (output) matrix
% Q - state noise covariance
% R - observation noise covariance
% x0 - initial state mean
% P0 - initial state covariance

%%% Matlab script to simulate data and process usiung Kalman for the state
%%% estimation of AR(2) time series.
% Linear system representation
% x_n+1 = A x_n + Bw_n
% y_n = Cx_n + v_n
% w = N(0,Q); v = N(0,R)
clc
clear all

T = 100; % number of data samples
order = 2;
% True coefficients of AR model
  a1 = 0.195;
  a2 = -0.95;

A = [ a1  a2;
      1  0 ];
C = [ 1 0 ];
B = [1;
      0];

 x =[ rand(order,1) zeros(order,T-1)];



sigma_2_w =1;  % variance of the excitation signal for driving the AR model(process noise)
sigma_2_v = 0.01; % variance of measure noise


Q=eye(order);
P=Q;




%Simulate AR model time series, x;



sqrtW=sqrtm(sigma_2_w);
%simulation of the system
for t = 1:T-1
    x(:,t+1) = A*x(:,t) + B*sqrtW*randn(1,1);
end

%noisy observation

y = C*x + sqrt(sigma_2_v)*randn(1,T);

%R=sigma_2_v*diag(diag(x));
%R = diag(R);

R = var(y);
z = zeros(1,length(y));
z = y;

 x0=mean(y);
for i=1:T-1
[xpred, Ppred] = predict(x0,P, A, Q);
[nu, S] = innovation(xpred, Ppred, z(i), C, R);
[xnew, P] = innovation_update(xpred, Ppred, nu, S, C);
end

%plot
xhat = xnew';


plot(xhat(:,1),'red');
hold on;
plot(x(:,1));



function [xpred, Ppred] = predict(x0,P, A, Q)
xpred = A*x0;
Ppred = A*P*A' + Q;
end

function [nu, S] = innovation(xpred, Ppred, y, C, R)
nu = y - C*xpred; %% innovation
S = R + C*Ppred*C'; %% innovation covariance
end

function [xnew, Pnew] = innovation_update(xpred, Ppred, nu, S, C)
K = Ppred*C'*inv(S); %% Kalman gain
xnew = xpred + K*nu; %% new state
Pnew = Ppred - Ppred*K*C; %% new covariance
end

问题:我想通过绘图查看估计状态xnew与实际状态x有多接近.但是，函数innovation_update返回的xnew是2by2矩阵！如何用估计值模拟时间序列?

Question: I want to see how close the estimated state xnew is to the actual state x by a plot. But, the xnew returned by the function innovation_update is a 2by2 matrix! How do I simulate a time series with the estimated values?

推荐答案

您无需将x初始化为任何东西，只需设置初始状态x(:,1)和系统仿真"循环将填充其余部分. 糟糕，我知道您已经在这样做了！

You don't need to initialise x to anything, just set the initial state x(:,1) and the "simulation of the system" loop will fill in the rest. Oops, I see you were already doing that!

稍后，在从嘈杂的观测结果y推断状态xnew的循环中，您可以添加以下行:

Later, in the loop that infers the state xnew from the noisy observations y you can add the line:

[xnew, P] = innovation_update(xpred, Ppred, nu, S, C);
yhat(i) = C*xnew; % Observed value at time step i, assuming inferred state xnew

最后，您应该绘制yhat和y进行比较.

Finally, you should plot yhat and y for comparison.

如果要为估计的不确定性添加误差线，则还应该存储Phat(i) = sqrt(C*P*C')并调用errorbar(yhat,Phat)而不是plot(yhat).

If you want to add error bars for the uncertainty on the estimate then you should also store Phat(i) = sqrt(C*P*C') and call errorbar(yhat,Phat) instead of plot(yhat).

这篇关于Matlab:在通过卡尔曼滤波器进行状态估计后，如何模拟模型的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！