粒子群算法（PSO）的Matlab实现

1、粒子群算法的概念

粒子群优化算法(PSO：Particle swarm optimization) 源于对鸟群捕食的行为研究。粒子群优化算法的基本思想：是通过群体中个体之间的协作和信息共享来寻找最优解，每个个体对比最佳位置，得出群体最佳位置。

2、算法分析

粒子群算法通过无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度v=(v1,v2,....vn)和位置x(x1,x2,....xn),速度代表移动的快慢。每个粒子在搜索空间中单独的搜寻最优解，并将其标记为当前个体极值，并将个体极值与整个粒子群里的其他粒子共享，在所有粒子中找到最优的那个个体极值作为整个粒子群的当前全局最优解，并返回拥有全局最优解的粒子位置，粒子群中的所有粒子根据自己找到的当前个体极值的位置和整个粒子群共享的当前全局最优解位置来调整自己的移动方向，在迭代中根据方向和速度调整位置。

3、更新规则

PSO初始化为一群随机粒子的位置(随机解)。然后通过迭代目标函数找到最优解的位置。在每一次的迭代中，粒子通过跟踪两个“极值”(个体极值pbest，全局极值gbest)来更新自己。在找到这两个最优值后，粒子通过下面的公式来更新自己的速度和位置。

公式(1)的第一部分称为【惯性项】，表示上次速度大小和方向的影响；第二部分称为【自身认知项】，是从当前点指向粒子自身最好点的一个矢量，表示粒子的动作来源于寻找自身的最优值；公式(1)的第三部分称为【群体认知项】，是一个从当前点指向种群最好点的矢量，反映了粒子间的协同合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。以上面两个公式为基础，形成了PSO的标准形式。

g为当前迭代次数。其算法流程如下：

4、算例

y=x1^2+x2^2-10cos(2*pi*x1)-10cos(2*pi*x2)+20

clc
clear 
close all

E=0.000001;
maxnum=800;%最大迭代次数
narvs=2;%目标函数的自变量个数，空间维度 N 这里面为X1和X2，二维变量
particlesize=50;%粒子群规模
c1=2;%每个粒子的个体学习因子，加速度常数
c2=2;%每个粒子的社会学习因子，加速度常数
w=0.6;%惯性因子
vmax=5;%粒子的最大飞翔速度
limit = [-5.12,5.12]; 
v=2*rand(particlesize,narvs);%初始化粒子飞翔速度，50行2列的速度矩阵
x=limit(1) + (  limit( 2 ) -  limit( 1)  ).*rand(particlesize,narvs);%初始化粒子位置，50行2列

%定义适应度函数
fitness=@(x,y)   20 +  x.^2 + y.^2 - 10*cos(2*pi.*x)  - 10*cos(2*pi.*y);%即求解的函数
x0=[-5.12:0.05:5.12];
y0=x0;
[X,Y]=meshgrid(x0,y0);
Z=fitness(X,Y);
figure(1);mesh(X,Y,Z);
for i=1:particlesize
	f(i)=fitness(x(i,1),x(i,2));	
end
personalbest_x=x;
personalbest_faval=f;
[globalbest_faval,i]=min(personalbest_faval);%求取适应度最小的值和粒子
globalbest_x=personalbest_x(i,:); %最佳粒子位置
k=1;
%%
while (k<=maxnum)
	for i=1:particlesize
			f(i)=fitness(x(i,1),x(i,2));
		if f(i)<personalbest_faval(i)
			personalbest_faval(i)=f(i);
			personalbest_x(i,:)=x(i,:);
		end
	end
	[globalbest_faval,i]=min(personalbest_faval);
	globalbest_x=personalbest_x(i,:);
	for i=1:particlesize
		v(i,:)=w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:))...
			+c2*rand*(globalbest_x-x(i,:));
		for j=1:narvs
			if v(i,j)>vmax
				v(i,j)=vmax;
			elseif v(i,j)<-vmax
				v(i,j)=-vmax;
			end		
		end
		x(i,:)=x(i,:)+v(i,:);
    end
    ff(k)=globalbest_faval;
    if globalbest_faval<E
        break
    end
	k=k+1;
end
%xbest=globalbest_x;
figure(2)
plot(1:length(ff),ff)
title("自适应度")
disp(['最优值：',num2str(globalbest_faval)]);  
disp(['变量取值：',num2str(globalbest_x)]);