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

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

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

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

🔥 内容介绍 

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

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

二、分数阶微积分基础

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

1. 粒子表示

   ：FPF 基于蒙特卡罗方法，通过一组随机采样的粒子来近似系统状态的概率分布。每个粒子表示一个可能的系统状态，粒子的权值反映了该状态出现的概率。在非线性离散时间分数阶系统中，根据系统的动态方程和观测方程，对粒子进行采样和更新。

2. 分数阶动态模型

   ：考虑分数阶系统的动态特性，在粒子的状态更新过程中引入分数阶导数运算。根据分数阶系统的状态转移方程，利用分数阶微积分计算粒子状态的变化。例如，根据 Caputo 定义的分数阶导数，更新粒子的位置、速度等状态信息。

3. 权值更新与估计

   ：根据观测值，利用贝叶斯准则更新粒子的权值。观测值与粒子预测状态之间的差异越小，粒子的权值越大。通过不断迭代，粒子逐渐集中在真实状态附近，最终通过对粒子及其权值的统计计算得到系统状态的估计值，如加权平均等方法。

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.