%% 粒子群 PSO 线性权重
clear
clc
tic
%% 绘制原图
x1=-5:0.2:5;
x2=-5:0.2:5;
[x1,x2]=meshgrid(x1,x2);
x1=x1+eps;%eps=2.2204e-16是一个小量，后续计算过程，x,y出现在分母上，为防止x,y等于0，加上一个小量
x2=x2+eps;
% y=(sin(x1).*sin(x2))./(x1.*x2);
y =3*(1-x1).^2.*exp(-(x1.^2) - (x2+1).^2) ...
- 10*(x1/5 - x1.^3 - x2.^5).*exp(-x1.^2-x2.^2) ...
- 1/3*exp(-(x1+1).^2 - x2.^2);
mesh(x1,x2,y)
axis vis3d
% shading interp
hold on
%% 预制参数
n=100;%粒子数量
d=2;%变量个数/变量维度
c1=2;
c2=2;
w=0.9;
ws=0.9;
we=0.4;
K=100;
%% 分布粒子
x=-5+10*rand(n,d);%rand(n,d)函数可随即分布范围（0,1）
v=-5+10*rand(n,d);
%% 计算适应度
fit=zeros(n,1);
for j=1:n
fit(j)=demo_11_1(x(j,:));
end
pbest=x;
ind=find(min(fit)==fit);
gbest=x(ind,:);
h=scatter3(x(:,1),x(:,2),fit,'bo');
%% 更新速度与位置
for i=1:K
%更新速度与位置
for m=1:n
v(m,:)=w*v(m,:)+c1*rand*(pbest(m,:)-x(m,:))+c2*rand*(gbest-x(m,:));
w=ws-(ws-we)*i/K;%线性递减权重
v(m,find(v(m,:)<-5))=-5;
v(m,find(v(m,:)>5))=5;
x(m,:)=x(m,:)+v(m,:);
x(m,find(x(m,:)<-5))=-5;
x(m,find(x(m,:)>5))=5;
%重新计算适应度
fit(m)=demo_11_1(x(m,:));
if fit(m)<demo_11_1(pbest(m,:))
pbest(m,:)=x(m,:);
end
if demo_11_1(pbest(m,:))<demo_11_1(gbest)
gbest=pbest(m,:);
end
end
fitnessbest(i)=demo_11_1(gbest);
pause(0.01);
h.XData=x(:,1);
h.YData=x(:,2);
h.ZData=fit;
end
demo_11_1(gbest)
% plot(fitnessbest)
% xlabel('迭代次数');
% ylabel('适应度值')
% title('迭代优化过程')
toc
