✅作者简介：热爱科研的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
end
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.
