【卡尔曼滤波】粗略模型和过滤技术在模型不确定情况下的应用研究（Matlab代码实现）

 💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码及文献

---

💥1 概述

摘要—本文讨论了在只有粗略模型可用时系统的估计和滤波问题。基于经典正则化最小二乘问题的修改版本，提出了一种新的设计准则，考虑了测量和创新作为不确定性的可能来源。在高斯假设下，它表现为最大后验贝叶斯估计的上界。通过利用非光滑分析工具获得最优解，最优解揭示了残差空间中的一个区域，其中估计的非变化是最优的。该方法以递归形式从随机角度提供了鲁棒估计器。为了说明，导出了类似卡尔曼滤波器，并与经典的最坏情况鲁棒设计滤波器进行了比较。
索引词—估计，不确定系统，卡尔曼滤波。
一、引言
从噪声观测中估计一组未知参数的问题在科学和工程中有着广泛的应用。然而，通常情况下，数据生成的精确模型可能不为人所知，因此需要一个鲁棒的估计器。
在文献中，主要有两种主要的鲁棒状态估计方法。第一种方法是H∞估计，试图最小化输入噪声对估计误差的最坏情况能量增益，参见[1]，[2]。第二种方法是所谓的保证成本估计，其中设计了一个满足上限误差的估计算法，参见[3]，[4]。此外，还有一种基于确定性正则化最小二乘问题的不同方法已被探索用于设计鲁棒估计器，参见[5]，[6]，[7]。然而，这些方法涉及基于最坏情况的优化，如果模型和测量的不确定性相当大，解决方案中可能存在高度保守性。
为了获得针对模型不确定性的鲁棒估计器，本文采用了较少保守性的随机观点，并依赖于经典正则化最小二乘问题的修改版本[8]。这使我们能够获得最小二乘鲁棒估计和卡尔曼滤波器的更新阶段。
该立场的主要贡献是一种用于模型不确定性的静态估计器和滤波器模型，采用了一种新颖的随机建模工具。结果表明，在高斯假设下，估计准则表现为最大后验贝叶斯估计的上界。
由于模型中嵌入了由估计变化的绝对值调制的随机噪声，这种滤波器的行为在某些情况下，最佳估计仍然是先验，而不是后验，考虑到系统动态提供的先验的信心和新测量的信心之间的权衡，这在现有的鲁棒滤波的文献中是不可见的，参见[5]，[6]，[9]，[10]。与这些不同的是，得到的滤波器增益是非线性的。

📚2 运行结果

主函数代码：

% This Matlab code implements the numerical example presented in the paper 
% M. R. Fernandes, J. B. R. do Val and R. F. Souto, "Robust Estimation and 
% Filtering for Poorly Known Models" in IEEE Control Systems Letters.
% doi: 10.1109/LCSYS.2019.2951611
% https://ieeexplore.ieee.org/document/8891731

close all
clear all
clc
rng('default'); rng(0)
%% System Model
n=2; %number of states
p=1; %number of outputs
A0=[0.9802 0.0196;0 0.9802]; %Nominal Dynamic Matrix
A1=[0 0.099;0 0];
H0=[1 -1]; %Measurement Matrix
G=eye(n);
Ef=[0 5];
Eg=[0 0];
M=[0.0198;0];

%% Noise Covariances
sigma=sqrtm([1.9608 0.0195;0.0195 1.9605]);
nu=eye(p);
Q = sigma*sigma';
R = nu*nu';

%% EVIU Parameters
sigmas.x=[0 0;0 0];
sigmas.bx=0.1*eye(n);
sigmas.y=[0 0.1];
sigmas.by=[0 0.1];
sigmas.v=[0 0.1];
sigmas.bv=sigmas.v;

%% Simulation
N=1000; %time horizon
a=-1:0.1:1; % grid for uncertain parameter
H=H0; %output matrix is considered exactly
gamma=99; %gain for the Hinf Filter

%Make a grid on [-1,1] for the uncertain parameter and apply simultaneosly
%the EVIU, BDU, Hinf and Kalman filter.
for i=1:length(a)
    A=A0+a(i)*A1; %define 'true' dynamic matrix
    % Monte Carlo Simulation
    for j=1:500
        %% Initialization
        clear x hx m hx2 hx_opt hx_s hx_hinf P P2 P_s P_hinf P_opt
        x(:,1)=1*ones(n,1);%real states
        hx(:,1)=zeros(n,1);%classic kalman states
        hx_opt(:,1)=zeros(n,1);%opt states
        m=zeros(n,1);
        S=zeros(n,1); %beta0
        hx2=hx; %eviu states
        hx_s=hx; %BDU states
        hx_hinf=hx; %hinf states
        P=eye(n); %initial covariance
        L=P;
        P_s=P; %covariance for BDU filter
        P2=P; %covariance for EVIU filter
        Phinf=P; %covariance for Hinf filter
        P_opt=P; %covariance for optimal filter
        trP(1)=trace(P); %trace of KF covariance
        trP2(j,1)=trP(1); %trace of EVIU covariance
        trP_opt(j,1)=trP(1); %trace of Optimal covariance
        trP_s(j,1)=trP(1); %trace of BDU covariance
        trP_hinf(j,1)=trP(1); %trace of Hinf covariance
        mse(j,1)=norm(hx(:,1)-x(:,1))^2; %KF error
        mse2(j,1)=mse(j,1); % EVIU error
        mse_opt(j,1)=mse(j,1); % Minimum error
        mse_hinf(j,1)=mse(j,1); % Hinf error
        mse_s(j,1)=mse(j,1); % BDU error
        %time simulation
        for k=1:N-1
            %% real system
            x(:,k+1)=A*x(:,k)+sigma*randn(n,1);
            y(:,k+1)=H*x(:,k+1)+nu*randn(p,1);
            %% Classic Kalman Filter (hx)
            [hx(:,k+1),P]=kalman_filter(A0,H0,Q,R,hx(:,k),P,y(:,k+1));
            trP(k+1)=trace(P);
            mse(j,k+1)=norm(x(:,k+1)-hx(:,k+1))^2;
            %% EVIU Filter (hx2)
        [hx2(:,k+1),P2,m,L,S(:,k+1)]=eviu_filter(A0,H0,sigmas,Q,R,hx2(:,k),P2,m,L,S(:,k),y(:,k+1));
            trP2(j,k+1)=trace(P2);
            mse2(j,k+1)=norm(x(:,k+1)-hx2(:,k+1))^2;
            %% Optimal Filter (Kalman with exactly model)
            [hx_opt(:,k+1),P_opt]=kalman_filter(A,H,Q,R,hx_opt(:,k),P_opt,y(:,k+1));
            trP_opt(j,k+1)=trace(P_opt);
            %% BDU Filter
            [hx_s(:,k+1),P_s]=BDU_filter(A0,H0,M,Ef,Eg,G,Q,R,hx_s(:,k),P_s,y(:,k+1));
            trP_s(j,k+1)=trace(P_s);
            mse_s(j,k+1)=norm(x(:,k+1)-hx_s(:,k+1))^2;
            %% Hinf filter
            [hx_hinf(:,k+1),Phinf]=hinf_filter(A0,H0,G,Q,R,hx_hinf(:,k),Phinf,y(:,k+1),eye(n),gamma);
            trP_hinf(j,k+1)=trace(Phinf);
            mse_hinf(j,k+1)=norm(x(:,k+1)-hx_hinf(:,k+1))^2;
        end
    end
    mse=mean(mse); %KF
    mse2=mean(mse2); %EVIU
    mse_s=mean(mse_s); %BDU
    trP_opt=mean(trP_opt); %Optimal
    mse_hinf=mean(mse_hinf); %Hinf
    mmse(i)=mean(mse); %KF
    mmse2(i)=mean(mse2); %EVIU
    mmse_s(i)=mean(mse_s); %BDU
    mmse_hinf(i)=mean(mse_hinf); %Hinf
    mmse_opt(i)=mean(trP_opt); % Optimal
    fprintf('%d gain=%.3f\n',i,(1-mmse2(i)/mmse(i))*100)
end
%% Plot results
figure(1)
p1=plot(a,sqrt(mmse),'*-','linewidth',2);
hold on
p2=plot(a,sqrt(mmse2),'-','linewidth',2);
p3=plot(a,sqrt(mmse_s),'o-','linewidth',2);
p5=plot(a,sqrt(mmse_hinf),'c-','linewidth',2);
p4=plot(a,sqrt(mmse_opt),'--k','linewidth',2);
legend([p1,p2,p3,p5,p4],'KF','EVIU','BDU','$\mathcal{H}_\infty$','Optimal','Interpreter','latex')
ylabel('RMSE_{total}')
xlabel('\delta')
axis([-1 1 0 22])
grid on
figure(2)
%Worst-case (delta=1)
p1=semilogx(10*log10(mse),'-','linewidth',2);
hold on
p2=semilogx(10*log10(mse2),'linewidth',2);
p3=semilogx(10*log10(mse_s),'-','linewidth',2);
p4=semilogx(10*log10(trP_opt),'--k','linewidth',2);
p5=semilogx(10*log10(mse_hinf),'c-','linewidth',2);
ylabel('MSE(dB)')
xlabel('Tempo (k)')
grid on
legend([p1,p2,p3,p5,p4],'KF','EVIU','BDU','$\mathcal{H}_\infty$','Optimal','Interpreter','latex')

🎉3 参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。

Marcos R. Fernandes, Joao Bosco R. do Val, Rafael. F. Souto.

🌈4 Matlab代码及文献