粒子群优化算法

 粒子群算法描述

  粒子群算法，也称粒子群优化算法或鸟群觅食算法，缩写为PSO，是一种新的进化算法。PSO算法属于进化算法中的一种，和模拟退火算法相似，它也是从随机解出发，通过迭代寻找最优解，通过适应度来评价解的品质，但是比遗传算法规则更为简单，并没有遗传算法的”交叉“和”变异“操作，它通过追随当前搜索到的最优质来寻找全局最优。这种算法以其容易实现、精度高、收敛快等优点引起重视，并在解决实际问题中展示了其优越性。粒子群算法是一种并行算法。

粒子仅有两个属性：速度和位置，速度代表移动的快慢，粒子代表移动的方向。

鸟被抽象为没有质量和体积的微粒（点），并延伸到N维空间，粒子i在N维空间的位置表示为矢量Xi = （x1，x2，…，xn），飞行速度表示为矢量Vi=（v1，v2，…，vn）。每个粒子都有一个由目标函数决定的适应值，并且知道自己到目前为止发现的最好位置（pbest）和现在的位置Xi。这个可以看作是粒子自己的飞行经验。除此之外，每个粒子还知道到目前为止整个群体中所有粒子发现的最好的位置（gbest），这个可以看作是粒子同伴的经验。粒子就是通过自己的经验和同伴中最好的经验来决定下一步运动。

速度和位置的更新

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

对于公式（1）：

公式（1）中的第一部分称为记忆项，表示上次速度大小和方向的影响；
公式（1）中的第二部分称为自身认知项，是从当前点指向粒子自身最好点的一个矢量，表示粒子的动作来源于自己经验的部分；
公式（1）中的第三部分称为群体认知项，是一个从当前点指向种群最好点的矢量，反映了粒子间的协调合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。

 标准PSO算法流程

1.初始化一群微粒（群体规模为N），包括随机位置和速度；
2.评价每个微粒的适应度；
3.对每个微粒，将其适应值与其经过的最好位置pbest作比较，如果较好，则将其作为当前的最好位置pbest；
4.对每个微粒，将其适应值与其经过的最好位置gbest作比较，如果较好，则将其作为当前的最好位置gbest；
5.根据公式（2）、（3）调整微粒的速度和位置；
6.未达到结束条件则转到第二步。
  迭代终止条件根据具体问题一般选为最大迭代次数Gk或微粒群迄今为止搜索到的最优位置满足预定最小适应阈值。

学习因子c1、c2分析
公式（2）、（3）中pbest和gbest分别表示微粒群的局部和全局最优位置。

当C1 = 0 时，则粒子没有了认知能力，变为只有社会的模型（social-only）：

称为全局PSO算法。粒子有扩展搜索空间的能力，具有较快的收敛速度，但由于缺少局部搜索，对复杂问题比标准PSO更易陷入局部最优。

c1 = 0 表示与单粒子本身的pbest无关了

当C2 = 0 时，则粒子之间没有社会信息，模型变为只有认知（cognition-only）模型。

此时称为局部PSO算法。由于个体之间没有信息交流，整个群体相当于多个粒子进行盲目的随机搜索，收敛速度满，因而得到最优解的可能性小。
 

操作思路

在D维空间中，有N个粒子；

粒子i的速度：xi = （xi1，xi2，…，xiD），将xi带入**适应函数f（xi）**求适应值；

粒子i的速度：vi = （vi1，vi2，…，viD）

粒子i个体经历的最好位置：pbest = （pi1，pi2，…，piD）

种群所经历的最好位置：gbest = （g1，g2，…，gD）

通常，在第d（1<=d<=D）维的位置变化范围限定在（Xmin d，Xmax d）内，速度变化范围限定在（-Vmax d，Vmax d）内（即在迭代中若超出了边界值，则该维的速度或位置被限制维该维最大速度或边界位置）

粒子i的第d维速度更新公式：

粒子i的第d维位置更新公式：

粒子其运动的过程受到三方面影响：粒子先前的速度、粒子认知部分、粒子社会部分

种群大小m
m很小容易陷入局部最优；m很大，PSO的优化能力很好，当种群数目增长到一定水平时，再增长将不再有显著的作用。

权重因子
对于粒子的速度更新的三部分：

1.惯性因子w=1表示基本的粒子群算法，w=0表示失去对粒子本身的速度的记忆；
2.自我认知部分的学习因子C1 = 0 表示无私型的粒子群算法，只有社会，没有自我，这样会使群体丧失多样性，从而导致容易陷入局部最优而无法跳出；
3.社会经验部分的学习因子C2=0表示自我型的粒子群算法，只有自我没有社会，这样导致没有信息的社会共享，算法收敛速度较慢。

最大速度

速度限制作用为：维护算法的探索能力与开发能力的平衡。

Vm较大时，探索能力强，但是粒子容易飞过最优解；

Vm较小时，开发能力强，但是粒子容易陷入局部最优解；

Vm一般设定为每维变量变化范围的10%–20%。

停止准则
最大迭代次数；

可以接受的满意解（通过fitness function判断是否满意）

粒子空间的初始化
较好地选择粒子的初始化空间，将大大缩短收敛时间。初始化空间根据具体问题的不同而不同，根据具体问题进行设定。该算法为数不多的关键参数的设置却对算法的精度和效率有着显著影响。

线性递减权值

Wmax最大惯性权重，Wmin为最小惯性权重，run为当前迭代次数，runmax为算法迭代总次数。

较大的w有较好的全局收敛能力，较小的w则有较强的局部收敛能力。因此随着迭代次数的增加，惯性权重w应不断减少，从而使得粒子群算法在初期具有较强的全局收敛能力，而晚期具有较强的局部收敛能力。

全局粒子群算法：

1.粒子自己历史最优值；
2.粒子群体的全局最优值。

局部粒子群算法：

1.粒子自己历史最优值；
2.粒子邻域内粒子的最优值。
邻域随迭代次数的增加先行增大，最后邻域扩展至整个粒子群。

实践证明：全局版本的粒子群算法收敛速度快,但是容易陷入局部最优。局部最优的粒子群算法收敛速度慢，但是很难陷入局部最优。现在粒子群算法大都在收敛速度与摆脱局部最优这两个方面下功夫。

PSO寻找函数最小值MATLAB代码

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)]);  

PSO在TSP问题中的应用

简单解释一下，m等于7即7个城市，k等于2，选取第二个，比较xi和xj，寻找到两个向量相等的位置的值，开始比较，第一轮发现3和1不等，那么将Xj中3和1的位置交换，开启下一轮。

代码后续补上