【状态估计】FEKF分数阶扩展卡尔曼滤波器、FCDKF分数阶中心差分卡尔曼滤波器、FUKF分数阶无迹卡尔曼滤波器和 FPF分数阶粒子滤波器的非线性离散时间分数阶系统状态估计附matlab代码

✅作者简介热爱科研的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); % 系统状态真实值 初始值0y zeros(1,tf); % 系统状态真实值 初始值0y(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 : Nxpart(i,1) x(1,1) sqrt(P(1,1)) * randn; %初始状态服从x0均值方差为sqrt(P)的高斯分布end% xArr [x];% yArr [];% xhatArr [x];% PArr [P];%xhatPartArr [xhatPart]; %%计算alpha阶次对应的GL定义系数 binomial coefficientbino_fir zeros(1,N); %微分阶次为0.7时GL定义下的系数alpha 1;bino_fir(1,1) 0.7;for i 2:1:Nbino_fir(1,i) (1-(alpha1)/(i-1))*bino_fir(1,i-1);end%%% diff_X_real 表示k时刻状态的微分diff_X_real 0;%% 开始循环for k 2 : tfdiff_X_real 3*sin(2*x(1,k-1)) -x(1,k-1) W_noise(1,k-1);rema 0;for i 2:1:krema rema bino_fir(1,i)*x(1,k1-i);endx(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 : ktemp temp bino_fir(1,j)*xpart(i,k1-j);endxpartminus(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 : Nq(i) q(i) / qsum; %归一化后的权值 qend%%%根据权值重新采样for i 1 : Nu rand;qtempsum 0;for j 1 : Nqtempsum qtempsum q(j);if qtempsum uxpart(i,k) xpartminus(j);break;elsexpart(i,k) xpart(i,k-1);endendendxhatPart(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 normalset(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.往期回顾扫扫下方二维码