当前位置：首页 > news >正文

分数阶扩展卡尔曼滤波器、分数阶中心差分卡尔曼滤波器、分数阶无迹卡尔曼滤波器和分数阶粒子滤波器的状态估计附matlab代码

news 2026/3/26 15:55:30

✅作者简介：热爱科研的Matlab仿真开发者，擅长毕业设计辅导、数学建模、数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。

🍎 往期回顾关注个人主页：Matlab科研工作室

👇 关注我领取海量matlab电子书和数学建模资料

🍊个人信条：格物致知,完整Matlab代码获取及仿真咨询内容私信。

🔥 内容介绍

一、状态估计的重要性与挑战

在众多工程和科学领域，如机器人运动控制、电力系统监测、生物医学信号处理等，准确估计系统的状态至关重要。对于非线性离散时间分数阶系统，其状态演化不仅涉及复杂的非线性关系，还具有分数阶导数的特性，这使得状态估计面临巨大挑战。传统的状态估计方法往往基于整数阶模型，难以准确描述此类系统的动态行为，因此需要专门针对非线性离散时间分数阶系统的估计方法。

二、分数阶微积分基础

六、FPF 分数阶粒子滤波器原理

粒子表示
：FPF 基于蒙特卡罗方法，通过一组随机采样的粒子来近似系统状态的概率分布。每个粒子表示一个可能的系统状态，粒子的权值反映了该状态出现的概率。在非线性离散时间分数阶系统中，根据系统的动态方程和观测方程，对粒子进行采样和更新。
分数阶动态模型
：考虑分数阶系统的动态特性，在粒子的状态更新过程中引入分数阶导数运算。根据分数阶系统的状态转移方程，利用分数阶微积分计算粒子状态的变化。例如，根据 Caputo 定义的分数阶导数，更新粒子的位置、速度等状态信息。
权值更新与估计
：根据观测值，利用贝叶斯准则更新粒子的权值。观测值与粒子预测状态之间的差异越小，粒子的权值越大。通过不断迭代，粒子逐渐集中在真实状态附近，最终通过对粒子及其权值的统计计算得到系统状态的估计值，如加权平均等方法。

FEKF、FCDKF、FUKF 和 FPF 分别从不同角度针对非线性离散时间分数阶系统进行状态估计，每种方法都有其独特的优势和适用场景，为解决此类系统的状态估计问题提供了多样化的选择。

⛳️ 运行结果

📣 部分代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 分数阶粒子滤波仿真复现

% 论文： fractional order PF

% 目的：分数阶粒子滤波算法测试

% 对系统噪声均值进行估计

% 函数实验: D^{0.7} x_k = 3*sin(2*x_{k-1}) -x_{k-1} + w_k

% y_k = x_k + v_k

% 结果：较好的对状态进行估计，常值系统噪声均值收敛

% 备注：分数阶粒子滤波的算法测试

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

clear all;

LineWidth = 1.5;

h_con = sqrt(3);

tf = 100; % 仿真时长

N = 100; % 粒子个数

%系统矩阵设置

% A = [0,1; -0.1,-0.2]; %系统矩阵

% B = [0; 1]; %

% C = [0.1,0.3]; %

I = eye(1,1); %生成单位阵

%I(3,3) = 0;

%噪声

q = 0; %系统噪声均值

r = 0; %测量噪声均值

Q = 0.81; %系统噪声方差矩阵

R = 0.25; %测量噪声方差矩阵

W_noise = sqrt(Q)*randn(1,N) + q; %系统噪声

V_noise = sqrt(R)*randn(1,N) + r; %测量噪声

x = zeros(1,tf); % 系统状态真实值初始值0

y = zeros(1,tf); % 系统状态真实值初始值0

y(1,1) = x(1,1) + sqrt(R) * randn;

P = zeros(1,tf); % 采样方差

P(1,1) = 2; % 初始采样分布的方差

xhatPart = zeros(1,tf); %状态估计值

xpart = zeros(N,tf);

for i = 1 : N

xpart(i,1) = x(1,1) + sqrt(P(1,1)) * randn; %初始状态服从x=0均值，方差为sqrt(P)的高斯分布

end

% xArr = [x];

% yArr = [];

% xhatArr = [x];

% PArr = [P];

%xhatPartArr = [xhatPart]; %

%计算alpha阶次对应的GL定义系数 binomial coefficient

bino_fir = zeros(1,N); %微分阶次为0.7时GL定义下的系数

alpha = 1;

bino_fir(1,1) = 0.7;

for i = 2:1:N

bino_fir(1,i) = (1-(alpha+1)/(i-1))*bino_fir(1,i-1);

end

% diff_X_real 表示k时刻状态的微分

diff_X_real = 0;

%% 开始循环

for k = 2 : tf

diff_X_real = 3*sin(2*x(1,k-1)) -x(1,k-1) + W_noise(1,k-1);

rema = 0;

for i = 2:1:k

rema = rema + bino_fir(1,i)*x(1,k+1-i);

end

x(1,k) = diff_X_real - rema;

%k时刻真实值

y(1,k) = x(1,k) + V_noise(1,k); %k时刻观测值

%% 采样N个粒子

for i = 1 : N

%采样获得N个粒子

xpartminus(i) = 3*sin(2*xpart(i,k-1)) - xpart(i,k-1) + sqrt(Q) * randn;

temp = 0;

for j = 2 : 1 : k

temp = temp + bino_fir(1,j)*xpart(i,k+1-j);

end

xpartminus(i) = xpartminus(i) - temp;

ypart = xpartminus(i); %每个粒子对应观测值

vhat = y(1,k) - ypart; %与真实观测之间的似然

q(i) = (1 / sqrt(R) / sqrt(2*pi)) * exp(-vhat^2 / 2 / R);

%每个粒子的似然即相似度

end

%权值归一化

qsum = sum(q);

for i = 1 : N

q(i) = q(i) / qsum; %归一化后的权值 q

end

%根据权值重新采样

for i = 1 : N

u = rand;

qtempsum = 0;

for j = 1 : N

qtempsum = qtempsum + q(j);

if qtempsum >= u

xpart(i,k) = xpartminus(j);

break;

else

xpart(i,k) = xpart(i,k-1);

end

xhatPart(1,k) = mean(xpart(:,k));

%最后的状态估计值即为N个粒子的平均值，这里经过重新采样后各个粒子的权值相同

% xArr = [xArr x];

% yArr = [yArr y];

% % xhatArr = [xhatArr xhat];

% PArr = [PArr P];

% xhatPartArr = [xhatPartArr xhatPart];

end

t = 1 : tf;

figure;

plot(t, x, 'r', t, xhatPart, 'b--','linewidth',LineWidth);

Esitimated_state = legend('Real Value','Estimated Value','Location','best');

set(Esitimated_state,'Interpreter','latex')

set(gcf,'Position',[200 200 400 300]);

axis([0 50 -6 6]) %设置坐标轴在指定的区间

axis normal

set(gca,'FontSize',10);

xlabel('time step','FontSize',7);

ylabel('state','FontSize',7);

%设置坐标轴刻度字体名称，大小

set(gca,'FontName','Helvetica','FontSize',8)

%title('Fractional particle filter')

%xhatRMS = sqrt((norm(x - xhat))^2 / tf);

%xhatPartRMS = sqrt((norm(xArr - xhatPartArr))^2 / tf);

% figure;

% plot(t,abs(x-xhatPart),'b');

% title('The error of FPF')

% t = 0 : tf;

% figure;

% plot(t, xArr, 'b-.', t, xhatPartArr, 'k-');

% legend('Real Value','Estimated Value');

% set(gca,'FontSize',10);

% xlabel('time step');

% ylabel('state');

% title('Particle filter')

% xhatRMS = sqrt((norm(xArr - xhatArr))^2 / tf);

% xhatPartRMS = sqrt((norm(xArr - xhatPartArr))^2 / tf);

% figure;

% plot(t,abs(xArr-xhatPartArr),'b');

% title('The error of PF')

🔗 参考文献

[1] Liu T , Wei Y , Yin W ,et al.State estimation for nonlinear discrete-time fractional systems: A Bayesian perspective[J].Signal processing, 2019, 165(Dec.):250-261.DOI:10.1016/j.sigpro.2019.06.037.

查看全文

http://www.jsqmd.com/news/529883/