Matlab: ошибка в обратной операции при реализации фильтра Калмана

Я пытаюсь реализовать базовые уравнения для фильтра Калмана для следующей одномерной модели AR:

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;
      0 1 ];
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);

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

Код выдает ошибку

Error using inv
Matrix must be square.

Error in innovation_update (line 2)
K = Ppred*C'*inv(S); %% Kalman gain

Error in Kalman_run (line 65)
[xnew, P] = innovation_update(xpred, Ppred, nu, S, C);

Это потому, что матрица S не квадратная. Как сделать его квадратным? Есть ли проблема с шагами до этого?


person SKM    schedule 09.03.2015    source источник
comment
Что касается ряда других ваших вопросов, я откатил редактирование здесь - ответ ниже был дан на основе вопроса с кодом, и возможно, что ответ не был бы дан, если бы код не был там.   -  person halfer    schedule 20.06.2016


Ответы (1)


Насколько я понимаю, матрица R должна быть ковариационной матрицей для шума измерений.

Следующие строки:

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

Измените R с диагональной матрицы 2x2 на вектор-столбец 2x1.

Поскольку ваше наблюдение y является скаляром, шум наблюдения v также должен быть скаляром. Это означает, что R представляет собой ковариационную матрицу 1x1 или просто дисперсию случайной величины v.

Вы должны сделать R скаляром, чтобы ваш код работал правильно.

person eigenchris    schedule 09.03.2015
comment
Спасибо за быстрый ответ. Причина, по которой я сделал диагонализацию, заключается в том, что в исходной реализации github.com /cswetenham/pmr/blob/master/toolboxes/lds/ , в строках 49-51 авторы используют диагональ ковариационной матрицы R. My inv(S) — это их temp1. Я не понимаю, почему они это сделали, но поскольку Авторы являются первопроходцами в этой области, я следовал их процедуре. Их код не выдает эту ошибку. Я пытаюсь реализовать это сам, так как считаю, что это лучший способ научиться, а не копировать. Если можно, не могли бы вы взглянуть на это? - person SKM; 09.03.2015
comment
@SKM Я пытаюсь сказать, что проблема в строке R=diag(R);, и вы должны избавиться от нее. Перед этой строкой R уже является диагональной диагональю, ее можно оставить как есть. Когда diag применяется к матрице, матрица превращается в вектор-столбец, поэтому, используя diag, вы берете R, который вам нужен, и превращаете его во что-то, чего вы не хотите. Имеет ли это смысл? Что касается логики всего этого, то, боюсь, я знаком с фильтром Калмана только мимоходом и не понимаю математических деталей. Я просто смотрел на форму уравнений. - person eigenchris; 09.03.2015
comment
Поскольку наблюдение y является одномерным вектором, R будет скаляром. - person SKM; 09.03.2015
comment
@SKM Это имеет смысл. Тогда вам придется изменить свое определение R, чтобы оно было скаляром, верно? Это будет дисперсия v? - person eigenchris; 09.03.2015
comment
Я сгенерировал зашумленные наблюдения из дисперсии v = 0,01. Если R = дисперсия v, то в чем разница? Или они должны быть одинаковыми (эту часть я не уловил из теории). Я запускал код с R = 1 и с R = 0,01. Создание R как скаляров не возвращает предыдущее. ошибка. Но алгоритм не работает, оценки состояния ужасно неверны. Что вы думаете об этом? Спасибо. - person SKM; 09.03.2015
comment
Код не работает должным образом, график неверный. xnew или xhat — это матрица 2 на 2. Разве это не должен быть вектор временного ряда, иначе как мне построить? - person SKM; 09.03.2015
comment
@SKM Если y представляет собой вектор-столбец m x 1, то R должен быть ковариационной матрицей m x m. Поскольку здесь m=1, R действительно должно быть скаляром, поскольку это ковариационная матрица v. Боюсь, я не смогу вам помочь, не прочитав больше о фильтрах Калмана. Тем временем вы можете попробовать прочитать это. Я не уверен, что не так с вашей моделью. - person eigenchris; 09.03.2015