%基本粒子群优化算法（Particle Swarm Optimization）-----------
%名称：基本粒子群优化算法（PSO）
%作用：求解优化问题
%说明：是2维变量的求解,程序中做了少量修改，可以准确迭代出结果
%作者：周楠
%时间：2009-12-18-2:00
clear all;
clc;
%横坐标长度百分比
x_rate=0.02;
%纵坐标长度百分比
y_rate=0.1;
%字体显示位置横坐标长度百分比
x1_show=0.5;
x2_show=0.5;
%字体显示位置纵坐标长度百分比
y1_show=0.2;
y2_show=0.9;
%x1的范围[a1 b1];
a1=-3.0;
b1=12.1;
%x2的范围[a2 b2];
%X2[a2,b2]
a2=4.1;
b2=5.8;
c1=2; %学习因子1
c2=2; %学习因子2
%惯性权重
w_0=0.9;
w_end=0.4;
MaxDT=30; %最大迭代次数
D=2; %搜索空间维数（未知数个数）
N=50; %初始化群体个体数目
eps=10^(-16); %设置精度(在已知最小值时候用)，即运算精度
rate=0.2; %如果步长超出边界，则缩短为原来的rate
%初始化种群的个体(位置和速度)
%rand是0-1的均匀分布
rr=rand(2,N);
r1=rr(1,:);
r2=rr(2,:);
%二维初始位置变量生成
xx=rand(N,2);
x1=xx(:,1);
x2=xx(:,2);
x1=(b1-a1)*x1+a1;
x2=(b2-a2)*x2+a2;
x_x=[x1,x2];
%速度
v1=zeros(N,N);
v2=v1;
%某次迭代的最优解初始化
pbest=zeros(MaxDT,N);
uu=fx1Tox2(x_x(1,:));
pbest(1,1)=uu(1,1);
%某次迭代的最优解的位置初始化
x_pbest=zeros(MaxDT,N,D);
x_pbest(1,1,:)=x_x(1,:);
%gbest为全局最优解的初始化
gbest=pbest(1,1);
draw_gbest=pbest(1,1); %draw_gbest为画出全局最优解的初始化
%gbest为全局最优解的位置初始化
x_gbest=zeros(1,D);
x_gbest(1,:)=x_x(1,:);
xx_x=zeros(1,D);%局部变量，用来存储临时变量
%进入主要循环，按照公式依次迭代，直到满足精度要求
for t=1:MaxDT
for i=1:N
for j=1:N
w(i)=((w_0-w_end)*(MaxDT-i))/(MaxDT)+w_end;
v1(i,j)=w(i)*v1(i,j)+c1*r1(j)*(x_pbest(t,j)-x_x(i,1))+c2*r2(j)*(x_gbest(1)-x_x(i,1));
v2(i,j)=w(i)*v2(i,j)+c1*r1(j)*(x_pbest(t,j)-x_x(i,2))+c2*r2(j)*(x_gbest(2)-x_x(j,2));
x_x(i,1)=x_x(i,1)+v1(i,j);
x_x(i,2)=x_x(i,2)+v2(i,j);
%对于步长太长，导致变量越界后的处理
while (((x_x(i,1))>b1)||((x_x(i,2))>b2)||((x_x(i,1))<a1)||((x_x(i,2))<a2))
%方法一：如果步长超出边界，则缩短为原来的rate
% x_x(i,1)=x_x(i,1)-v1(i,j);
% x_x(i,2)=x_x(i,2)-v2(i,j);
% v1(i,j)=v1(i,j)*rate;
% v2(i,j)=v2(i,j)*rate;
% x_x(i,1)=x_x(i,1)+v1(i,j);
% x_x(i,2)=x_x(i,2)+v2(i,j);
%方法二：把如果步长超出边界，则取边界附近某 个值
if ((x_x(i,1))>b1)
x_x(i,1)=x_x(i,1)-rate;
end
if ((x_x(i,2))>b2)
x_x(i,2)=x_x(i,2)-rate;
end
if (( x_x(i,1)<a1))
x_x(i,1)=x_x(i,1)+rate;
end
if (( x_x(i,2)<a2))
x_x(i,2)=x_x(i,2)+rate;
end
end
%保存D个值，保存第t次迭代的第i次局部最优解
xx_x(:,1)=x_x(i,1);
xx_x(:,2)=x_x(j,2);
if (fx1Tox2(xx_x)>pbest(t,i))
pbest(t,i)=fx1Tox2(xx_x);
x_pbest(t,i,:)=xx_x;
end
end
%保存全局最优解
if (pbest(t,i)>gbest)
gbest=pbest(t,i);
x_gbest=x_pbest(t,i,:);
end
end
end
%draw_gbest为画出全局最优解
draw_gbest=max(pbest');
mindle=draw_gbest(1);
for t=1:MaxDT
if (draw_gbest(t)>mindle)
mindle=draw_gbest(t);
else
draw_gbest(t)=mindle;
end
end
%最后给出计算结果
disp('函数的全局最优位置为：')
Solution=x_gbest
disp('最后得到的优化极值为：')
Result=gbest
figure
x=(1:MaxDT)';
plot(x,pbest,'b:');
title('粒子群优化算法')
xlabel('迭代次数')
ylabel('优化极值')
aa=min(min(pbest));
bb=max(gbest);
string_cell{1}='/fontsize{10}/fontname{arial}绿实线：到目前为止遇到的最好的解';
text((MaxDT+MaxDT*x_rate)*x1_show,((aa+aa*y_rate)+((bb+bb*y_rate)-(aa+aa*y_rate))*y1_show),string_cell,'color','m')
string_cell{1}='/fontsize{10}/fontname{arial}红虚线：当前代中最好的个体';
text((MaxDT+MaxDT*x_rate)*x1_show,((aa+aa*y_rate)+((bb+bb*y_rate)-(aa+aa*y_rate))*y1_show*y1_show),string_cell,'color','m')
cell_string{4}=['/fontsize{15}/rm全局最优位置为= ' num2str(x_gbest,6)];
cell_string{5}='/fontsize{15}';
cell_string{6}=['/fontsize{15}/rm全局优化极值= ' num2str(gbest,6)];
cell_string{7}='/fontsize{15}';
cell_string{8}=['/fontsize{15}/rm迭代次数= ' num2str(MaxDT,6)];
text((MaxDT+MaxDT*x_rate)*x2_show,((aa+aa*y_rate)+((bb+bb*y_rate)-(aa+aa*y_rate))*y2_show),cell_string,'color','b','HorizontalAlignment','center')
hold on
plot(x,draw_gbest,'r-')
axis([0,fix(MaxDT+MaxDT*x_rate),fix(aa+aa*y_rate),fix(bb+bb*y_rate)])
hold off
grid on