查看: 467|回复: 4|关注: 0

[已答复] 粒子群优化粒子滤波的matlab仿真程序

[复制链接]

新手

11 麦片

财富积分


050


9

主题

14

帖子

0

最佳答案
  • 关注者: 1
发表于 2019-1-7 16:56:53 | 显示全部楼层 |阅读模式
粒子群优化粒子滤波的matlab仿真程序,求共享一份啊。
PSO-PF程序,感谢!!!!!!

新手

11 麦片

财富积分


050


9

主题

14

帖子

0

最佳答案
  • 关注者: 1
 楼主| 发表于 2019-1-7 17:01:32 | 显示全部楼层
clear all
clc
x = 0.1; % initial state 初始状态
Q = 1; % process noise covariance 过程噪声方差
R = 1; % measurement noise covariance 测量噪声方差
tf = 50; % simulation length 模拟长度
N = 100; % number of particles in the particle filter 颗粒过滤器中的粒子数
xhat = x; %xhat=x=0.1
P = 2;
xhatPart = x; %xhatPart=x=0.1
%粒子群初始化部分
nn=100; %粒子群规模,即粒子数 %粒子的维数,即寻优参数的个数
X_pso=rand(nn); %初始化粒子群的位置(因为x的取值范围在-100到100之间)
V=rand(nn); %初始化粒子群的速度
Pbest=rand(nn); %粒子群的个体最优解
Gbest=rand(1); %粒子群的全体最优解
Wmax=0.9;
Wmin=0.4;
ITERmax=100; %进化代数
iter=1;
C1=2; %学习因子
C2=2; %学习因子
Xmax=1;
%以上程序到上一个%号为粒子群的初始化部分


% Initialize the particle filter. 初始化粒子滤波,xpart值用来在不同时刻生成粒子
for i = 1 : N
xpart(i) = x + sqrt(P) * randn; % randn产生标准正态分布的随机数或矩阵的函数。
end %初始化xpart(i)为生成的100个随机粒子
xArr = [x]; %xArr=x=0.1
xhatPartArr = [xhatPart]; %xhatPartArr = [xhatPart]=0.1
close all;
for k = 1 : tf %tf为时间长度,k可以理解为时间轴上的k时刻
% System simulation系统仿真
% x数据为时刻k的真实状态值
x = 0.5 * x + 25 * x / (1 + x^2) + 8 * cos(1.2*(k-1)) + sqrt(Q) * randn; %状态方程(1)
y = x^2 / 20 + sqrt(R) * randn;%观测方程(2)。观测方程是在观测值和待估参数之间建立的函数关系式。
% Particle filter 生成100个粒子并根据预测和观测值差值计算各个粒子的权重
for i = 1 : N
xpartminus(i) = 0.5 * xpart(i) + 25 * xpart(i) / (1 + xpart(i)^2) + 8 * cos(1.2*(k-1)) + sqrt(Q) * randn;
ypart = xpartminus(i)^2 / 20;
X_pso(i)=ypart; %这一步实现了将每一个预测值赋给粒子群,作为粒子群优化的基本粒子。
end
%粒子群对于100个观察值进行i代,使得这100个预测值向最优方向移动(一个预测值与一个粒子群的粒子对应)
while(iter<=ITERmax)
for i=1:nn
R1=rand;
R2=rand;
% W_pso=Wmax-(Wmax-Wmin)/ITERmax*iter; %惯性因子
W_pso = Wmin + (Wmax-Wmin) * exp(-3*(iter/ITERmax).^2);
Vinter=abs(randn)*(Pbest(i)-X_pso(i))+abs(randn)*(Gbest(1)-X_pso(i));
Xinter=X_pso(i)+Vinter;
if Xinter>Xmax
X_pso(i)=Xmax;
elseif Xinter<-Xmax
X_pso(i)=-Xmax;
else
X_pso(i)=Xinter;
end
end
for i=1:nn
pres_X=fit_m(X_pso(i),y);
pres_Pbest=fit_m(Pbest(i),y);
pres_Gbest=fit_m(Gbest(1),y);
if pres_X>pres_Pbest %粒子的适应值越高,那么粒子的分布就越靠近真实状态,因此粒子群的优化方向是让粒子朝一个较大适应值方向走。
Pbest(i)=X_pso(i);
end
if pres_X>pres_Gbest
Gbest(1)=X_pso(i);
end
end
% err(iter)=pres_Gbest; %%%%%%%构造一个iter维的行向量,并将里面的每个元素值取为pres-Gbest
if pres_Gbest>=10^(-10) %粒子群最优值符合某个阀值的时候,退出循环,认为粒子已经分布在真实状态附近。
iter=ITERmax;
else
iter=iter+1;
end
end %粒子群循环结束。
%将粒子群优化后的粒子重新赋给100个观测值。利用每一个观测值通过下面for循环可得到较优的100个权重值
for i = 1 : N
ypart= X_pso(i);
vhat = y - ypart; %观测和预测的差
q(i) = (1 / sqrt(R) / sqrt(2*pi)) * exp(-vhat^2 / 2 / R); %根据差值给出100个粒子对应的权重(正态分布函数)重要性概率密度函数
end
% Normalize the likelihood of each a priori estimate. 每一个先验估计正常化的可能性
qsum = sum(q);
for i = 1 : N
q(i) = q(i) / qsum;%归一化权重,较大的权重除以qsum不为零,大部分较小的权重除以qsum后为零
end
% Resample.重新取样
for i = 1 : N
u = rand; % uniform random number between 0 and 1 0和1之间的均匀随机数
qtempsum = 0;
for j = 1 : N
qtempsum = qtempsum + q(j);
if qtempsum >= u
%重采样对低权重进行剔除,同时保留高权重,防止退化的办法
xpart(i) = xpartminus(j);
break;
end
end
end
% The particle filter estimate is the mean of the particles. 粒子滤波的估计是颗粒的平均值
xhatPart = mean(xpart); %经过粒子滤波处理后的均值
% Plot the estimated pdf's at a specific time. 绘制在特定的时间估计的概率密度函数
% Save data in arrays for later plotting
xArr = [xArr x];
xhatPartArr = [xhatPartArr xhatPart];
end
t = 0 : tf;
toc
figure;
plot(t, xArr, 'b.', t, xhatPartArr, 'g'); %此图2对应xArr为真值,xhatPartArr为粒子滤波值
xlabel('time step'); ylabel('state');
legend('True state', 'Particle filter estimate');
xhatPartRMS = sqrt((norm(xArr - xhatPartArr))^2 / tf);
disp(['Particle filter RMS error = ', num2str(xhatPartRMS)]);
各位大神这个程序里缺少适应度函数,但是pres_X=fit_m(X_pso(i),y);
pres_Pbest=fit_m(Pbest(i),y);
pres_Gbest=fit_m(Gbest(1),y);这个程勋里fit却没定义,算不出来求解怎么弄呢?

新手

5 麦片

财富积分


050


2

主题

17

帖子

0

最佳答案
发表于 2019-1-8 16:48:27 | 显示全部楼层
luluiszero 发表于 2019-1-7 17:01
clear all
clc
x = 0.1; % initial state 初始状态

还以为你发了个完整的

新手

11 麦片

财富积分


050


9

主题

14

帖子

0

最佳答案
  • 关注者: 1
 楼主| 发表于 2019-1-8 20:26:50 | 显示全部楼层
我也想发完整的呢

新手

5 麦片

财富积分


050


2

主题

5

帖子

0

最佳答案
发表于 2019-3-14 09:25:51 | 显示全部楼层
楼主,我用了你的方法,也把适应度函数加进去了,但是效果不对,你做出来了吗可以交流一下 我Q422721657  微信18840822252
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

快速回复 返回顶部 返回列表