具有增强状态向量的离散时间卡尔曼滤波器

Posted 2023-02-21

【中文标题】具有增强状态向量的离散时间卡尔曼滤波器

【英文标题】：Discrete time kalman filter with augmented state vector

【发布时间】：2015-04-26 19:04:12

【问题描述】：

我正在尝试为具有增强状态向量的状态空间模型实现离散时间卡尔曼滤波器。也就是说，状态方程由下式给出：

x(t) = 0.80x(t-1) + 噪声

在我的示例中只有潜在过程，因此 x(t) 的维度为 1x1。测量方程为：

Y(t) = AX(t) + 噪声

其中大写 Y(t) 的尺寸为 3x1。它包含三个观测序列，Y(t) = [y1(t),y2(t),y3(t)]'。观察到的两个系列偶尔会以系统的方式缺失：

y1

(t) 在所有时间点都被观察到。

y2

(t) 在每四个时间点观察一次。

y3

(t) 在每十二个时间点观察一次。

度量矩阵A的维度是3x12，大写的维度X(t)是12x1。这是一个与潜在过程的最后十二个实现堆叠在一起的向量：

X(t) = [x(t),x(t-1),...,x （t-11）]'。

大写X(t)叠加的原因是y1(t)与x(t)和x(t-1)。此外，y2(t) 与 x(t)、x(t-1)、x 相关(t-2) 和 x(t-3)。最后，y3(t) 与潜在过程的所有最后十二个实现有关。

随附的 matlab 代码模拟来自该状态空间模型的数据，随后通过具有增强状态空间向量 X(t) 的卡尔曼滤波器运行。

我的问题是，过滤（和预测）过程与真正的潜在过程有很大不同。这也可以从附图中看出。显然，我做错了什么，但我看不到它是什么？我的卡尔曼滤波器工作得很好，当状态向量没有堆叠时...希望有人能帮忙，谢谢！

    %-------------------------------------- %
    % AUGMENTED KALMAN FILTER WITH MISSINGS
    % ------------------------------------- %

    clear; close all; clc;

    % 1) SET PARAMETER VALUES %
    % ----------------------- %

    % Number of observations, number of latent processes and observed
    % processes.
    T = 2000;
    NumberOfLatent = 1;
    NumberOfObs = 3;

    % Measurement Matrix, A.
    A = [-0.40, -0.15, zeros(1,10);
          0.25,  0.25, 0.25, 0.25, zeros(1,8);
         ones(1,12)];

    % State transition matrix, Phi. 
    Phi = zeros(12,12);
    Phi(1,1) = 0.80;
    Phi(2:end,1:end-1) = eye(11);

    % Covariance matrix (for the measurement equation).            
    R = diag([2.25; 1.00; 1.00]);

    % State noise and measurement noise.
    MeasNoise = mvnrnd(zeros(1,NumberOfObs),R,T)';
    StateNoise = randn(NumberOfLatent,T+11);

    % One latent process (X) and three observed processes (Y).
    X = zeros(NumberOfObs,T+11);  
    Y = zeros(NumberOfObs,T); 

    % Set initial state.
    X0 = 0;

    % 2) SIMULATE DATA %
    % ---------------- %
    % Before Y (the observed processes) can be simulated we need to simulate 11
    % realisations of X (the latent process).  
    for t = 1:T+11
       if (t == 1)
            X(:,t) = Phi(1,1)*X0 + StateNoise(:,t);
       else
            X(:,t) = Phi(1,1)*X(1,t-1) + StateNoise(:,t);
       end
       if (t>=12)
            Y(:,t-11) = A*X(1,t:-1:t-11)' + MeasNoise(:,t-11);
       end
    end


    % 3) DELETE DATA SUCH THAT THERE ARE MISSINGS %
    % ------------------------------------------- %
    % First row is observed at all time points. The second rows is observed  at
    % every 4th. time point and the third row is observed at every 12th. time
    % point.
    for j = 1:3
        Y(2,j:4:end) = NaN;
    end
    for j = 1:11
        Y(3,j:12:end) = NaN;
    end


    % 4) RUN THE KALMAN FILTER WITH AUGMENTED STATE VECTOR %
    % ---------------------------------------------------- %

    % Initializing matrices.
    AugmentedStates = 12;
    X_predicted = zeros(AugmentedStates,T);                   % Predicted  states.
    X_filtered = zeros(AugmentedStates,T);            % Updated states.
    P_predicted = zeros(AugmentedStates,AugmentedStates,T);   % Predicted variances.
    P_filtered = zeros(AugmentedStates,AugmentedStates,T);   % Filtered variances.
    KG = zeros(AugmentedStates,NumberOfObs,T);       % Kalman gains.
    Q  = eye(AugmentedStates);                       % State noise for  augmented states.
    Res         = zeros(NumberOfObs,T);                       % Residuals.

    % Initial values:
    X_zero = ones(AugmentedStates,1);
    P_zero = eye(AugmentedStates);

    for t = 1:T

        % Initialisation of the Kalman filter.
        if (t == 1)
            X_predicted(:,t) = Phi * X_zero;                          %  Initial predicted state, dimension (AugmentedStates x 1).
            P_predicted(:,:,t) = Phi * P_zero * Phi'+ Q;              % Initial variance, dimension (AugmentedStates x AugmentedStates).
        else
         % #1 Prediction step.
           X_predicted(:,t) = Phi * X_filtered(:,t-1);                    %   Predicted state, dimension (AugmentedStates x 1).
           P_predicted(:,:,t) = Phi * P_filtered(:,:,t-1) * Phi'+ Q;      % Variance of prediction, dimension (AugmentedStates x AugmentedStates).
        end

        % If there is missings the corresponding entries in A and Y is set  to
        % zero.
        if sum(isnan(Y(:,t)))>0
            temp = find(isnan(Y(:,t)));                             
            Y(temp,t) = 0;                                         
            A(temp,:) = 0;              
        end

        % #3 Calculate the Kalman Gains and save them in the matrix KG.
        K = (P_predicted(:,:,t) * A')/(A * P_predicted(:,:,t) * A' + R);         % Kalman Gains, dimension (AugmentedStates x ObservedProcesses).
        KG(:,:,t) = K;  

        % Residuals.
        Res(:,t) = Y(:,t) - A*X_predicted(:,t);                

        % #3 Update step.
        X_filtered(:,t) =  X_predicted(:,t) + K * Res(:,t);                           % Updated state (AugmentedStates x 1).
        P_filtered(:,:,t) = (eye(AugmentedStates) - K * A) * P_predicted(:,:,t);    % Updated variance, dimension (AugmentedStatesstates x AugmentedStates).

         % Measurement Matrix, A, is recreated.
         A = [-0.40, -0.15, zeros(1,10);
               0.25,  0.25, 0.25, 0.25, zeros(1,8);
               ones(1,12)];

     end

【问题讨论】：

我看到的一个错误：当你删除 A 的一部分时，由于缺少观察，它们永远不会恢复。

我已经更正了我的问题。谢谢@Ben Jackson！不幸的是，这并不能解决问题。大事没有出现在我的原始脚本中，只是在我为这个问题创建的这个脚本中（不幸的是......）。

这一行K = (P_predicted(:,:,t) * A')/(A * P_predicted(:,:,t) * A' + R) ; % 卡尔曼增益，维度（AugmentedStates x ObservedProcesses）。尝试将大小为 2x2 的矩阵与大小为 3x3 的矩阵相加。这是不允许的。

【参考方案1】：

啊哈。我意识到错误是什么。 matlab代码中的状态转移矩阵（Phi）不正确。也就是说，在当前的matlab代码中是：

% State transition matrix, Phi. 
Phi = zeros(12,12);
Phi(1,1) = 0.80;
Phi(2:end,1) = 1

它实际上应该是这样的：

% State transition matrix, Phi. 
Phi = zeros(12,12);
Phi(1,1) = 0.80;
Phi(2:end,1:end-1) = eye(11)

如果没有这种调整，过滤器中的状态方程就没有任何意义！我已经在代码中实现了这一点。为了证明它确实有效，情节现在看起来像这样：

...或者至少：过滤器的性能得到了改进！