Matlab: Comment puis-je simuler le modèle après estimation de l'état de filtre de Kalman

Question

Je suis en train de mettre en œuvre les équations de base pour le filtre de Kalman pour le modèle 1 dimensions AR suivant:

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

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

Le modèle espace d'état KF:

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

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

w(t) = N(0,Q) 

v(t) = N(0,R) 

où

% 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 

question: Je veux voir à quel point l'estime état xnew est à l'état actuel x par un complot. Mais, le xnew renvoyé par la fonction innovation_update est une matrice 2by2! Comment simuler une série temporelle avec les valeurs estimées?

Answer 1

Vous n'avez pas besoin d'initialiser x à quoi que ce soit, il suffit de définir l'état initial x(:,1) et la boucle "Simulation du système" remplira le reste. Oups, je vois que tu faisais déjà ça!

Plus tard, dans la boucle qui infère l'état xnew des observations bruyantes y vous pouvez ajouter la ligne:

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

Enfin, vous devez tracer yhat et y à titre de comparaison.

Si vous souhaitez ajouter des barres d'erreur pour l'incertitude sur l'estimation, vous devez également enregistrer Phat(i) = sqrt(C*P*C') et appeler errorbar(yhat,Phat) au lieu de plot(yhat).