深入理解粒子群优化算法（PSO）：详解、三大变种及Java实践实例
引言：
粒子群优化算法简称PSO（Particle Swarm Optimization），PSO算法是在1995 年由美国社会心理学家Kennedy 和电气工程师Eberhart提出的 ,主要思想来源于对鸟类群体行为的研究,他们的模型和仿真算法主要利用了生物学家Heppner提出的模型。
文章目录
- 一。粒子群优化算法简介
- 1.PSO概念
- 2.群体智能概念
- 3.场景模型
- 二。基本PSO算法
- 1.公式
- 2.解析
- 三。带惯性权重w的pso算法
- 1.公式
- 2.解析
- 3.惯性权重调整策略
- 四。带压缩因子k的pso算法
- 1.公式
- 2.解析
- 五。Java实现PSO算法
- 1.算法流程图
- 2.代码部分
- 3.基准优化函数
- 参考文献
一。粒子群优化算法简介
1.PSO概念
什么是粒子群优化算法？
PSO算法与其他进化类算法类似，具有进化和群体智能的特点，算法模拟鸟群飞行觅食的行为，通过鸟之间的集体协作和竞争达到群体智能的目的。在PSO算法中，每个候选解称为一个“粒子”，若干个候选解就构成了鸟的群体。每个粒子没有重量和体积，通过目标函数确定他的适应值。每个粒子在解空间中运动，并由速度决定它的运动方向和距离，粒子通过追随自身的个体最好位置与群体的全局最好位置来动态调整自己的位置信息。
2.群体智能概念
什么是群体智能？
动物群比独行的动物更容易找到食物，因为发现猎物的机会大大增加，此外如果一群动物围住猎物，它们捕获猎物的机会也大大增加。一群动物比一只独处的动物看起来具有更大的威胁，更能降低遭遇攻击的可能性。一群个体合作不仅能改进它们在某件任务上的集合性能而且还能提高每一个个体的性能，这正是粒子群优化的原理。
我们不仅在动物的行为中清楚看到粒子群优化的原理，在人类的行为中也能看到，当我们试图改进某个任务的性能时，会根据一些基本观点调整我们的方法：
-
惯性：我们往往会保留过去已证明是成功的那些旧的方式。“我总是这样做，所以我还会继续这样做” -
受社会的影响：我们听到他人的成功后会试图仿照他人的办法，我们可能从书籍，或者互联网，
或者报纸读到他人成功的事迹。“如果那样做对他们管用，对我可能也管用。” -
受邻居的影响：我们从与自己亲近的人那里学到的最多，受朋友的影响会比社会影响的更多。
3.场景模型
鸟类捕食模型
鸟类会根据自己的意识以及群体的意识来改变自己的捕食线路，使之更有可能在该路线上发现食物，如下图所示：
粒子群算法就是根据鸟群来的，那我们就谈谈鸟群吧。正如上图所示：
鸟儿变化后的速度跟三个方面有关，首先就是与前一刻的速度有关，变化需要基准和代价的，不能随意变向和变大小（鸟儿不可能把自己速度一下提升为无穷大吧）肯定要考虑之前速度。当然，鸟儿肯定很想吃东西，它的印象中好像上次在某个方向上发现很好吃的虫子，这次那个方向会不会也有好吃的，哈哈，是要考虑一下。最后，它的队友们曾经在另一个方向上找到最好吃的东西，还很多，这次会不会还向上次一样运气那么好。当然它一部分队友说的也可以（局部最优）。它考虑一番后，作出了如图决定。
书上那个图也很经典，可以拿来对应对应，是不是更好理解它了：
注：本块区域来源于https://blog.****.net/weixin_43820463/article/details/89075133
二。基本PSO算法
1.公式
2.解析
其中 i 表示第 i 个粒子，t 表示第 t 次迭代，j 表示 第 j 维度。
第一部分表示粒子先前的自身速度：
第二部分表示个体认知部分：
第三部分表示社会认知部分：
c1和c2称为学习因子，c1调节粒子飞向自身最好位置方向的步长，c2调节粒子飞向全局最好位置的步长。r1和r2是（0,1）之间均匀分布的相互独立的随机数序列。
三。带惯性权重w的pso算法
1.公式
为了提高基本粒子群算法的收敛性能和避免算法陷入局部最优，提出了惯性权重这一概念，在进化方程中引入惯性权重因子w：
2.解析
各部分组成以及含义与基本PSO算法一致，只有w不同：
w起到了一个平衡全局搜索能力和局部搜索能力的作用
- w值较大时，全局搜索能力较强，局部搜索能力弱
- w值较小时，全局搜索能力较弱，局部搜索能力强
恰当的w值可以提高算法性能，提高寻优能力，减少迭代次数
3.惯性权重调整策略
线性递减策略
①典型线性递减策略
将w设置为0.9到0.4的线性下降，使得pso在开始时，搜索较大的区域，较快的定位最优解的大致位置，随着w逐渐减小，粒子速度减慢，开始精细的局部搜索。这种典型的惯性权重线性递减策略在目前是应用最为广泛，但是由于这种策略下，迭代初期的全局搜索能力强，如果在初期搜索不到最好点，那么随着w的减小，局部搜索能力加强，就易陷于局部极值。
如下图所示：
②线性微分递减策略
初始迭代的时候，w变化缓慢，有利于迭代是寻找满足条件的局部最优解，在接近最大迭代次数时，w变化较快，在寻找局部最优值之后能够快速地收敛逼近与全局最优解，提高运算效率。
如下图所示：
先增后减策略
这种先增后减的惯性权重，在前期有较快的收敛速度，而后期的局部搜索能力也不错，在一定程度上保持了递减和递增策略的优点，同时克服了一些缺点，相对提高算法性能。
如下图所示：
四。带压缩因子k的pso算法
1.公式
标准PSO算法的另-种形式是由Clerc于1999年提出的带压缩因子X的PSO算法。这是考虑到即使基于社会相互作用的算法表现出好的性能,但是研究者没能充分解释它是如何工作的;另外,算法的传统版本中有一些令人不满意的动态特性,特别是需要限制粒子的速度以控制粒子的轨迹。由此,Clerc提出了压缩因子的概念来控制PSO算法的收敛趋势。
2.解析
k为压缩因子
五。Java实现PSO算法
1.算法流程图
2.代码部分
Particle：
/**
* 粒子类
* @author mrpod2g
*
*/
public class Particle {
//维数
public static int dimension = 2;
//粒子的位置
public double[] X = new double[dimension];
//局部最好位置
public double[] pbest = new double[dimension];
//粒子的速度
public double[] V = new double[dimension];
//最大速度
public double Vmax = 2;
//适应值
public double fitness;
/**
* 根据当前位置计算适应值
* @return newFitness
*/
public double calculateFitness() {
//1.Ackley's function:Ackley函数广泛用于测试优化算法。如上图所示，其二维形式的特征是外部区域几乎平坦，中心处有一个大孔。该函数存在使优化算法（尤其是爬坡算法）陷入其许多局部最小值中的风险。
double newFitness = -20*Math.pow(Math.E,(-0.2*Math.sqrt(0.5*(X[0]*X[0]+X[1]*X[1]))))-Math.pow(Math.E,(0.5*(Math.cos(2*Math.PI*X[0])+Math.cos(2*Math.PI*X[1]))))+Math.E+20;
//2.Sphere function:球形函数具有d局部最小值，全局变量除外。它是连续的，凸的和单峰的。该图显示了其二维形式。
//double newFitness = X[0]*X[0]+X[1]*X[1];
//3.Rosenbrock function:Rosenbrock函数，也称为Valley或Banana函数，是基于梯度的优化算法的常见测试问题。在上面的图中以二维形式显示了它。
//该函数是单峰函数，全局最小值位于狭窄的抛物线形山谷中。但是，即使容易找到该山谷，也很难将其收敛到最小（Picheny等，2012）。
//double newFitness = 100*(Math.pow((X[1]-X[0]*X[0]),2))+Math.pow((X[0]-1), 2);
return newFitness;
}
/**
* 初始化自己的位置和pbest
*/
public void initialX() {
for(int i=0;i<dimension;i++) {
//该方法的作用是生成一个随机的int值，该值介于[0,100)的区间，也就是0到100之间的随机int值，包含0而不包含100。
X[i] = new Random().nextInt(100);
pbest[i] = X[i];
}
}
/**
* 初始化自己的速度
*/
public void initialV() {
for(int i=0;i<dimension;i++) {
//随机产生一个0~1的随机小数
double tmp = new Random().nextDouble();
V[i] = tmp*4+(-2);
}
}
}
PSO
/**
* PSO算法实现
* @author mrpod2g
*
*/
public class PSO {
private static double[] gbest;//全局最优位置
private static double gbest_fitness = Double.MAX_VALUE;//全局最优位置对应的fitness
private static int particle_num = 20;//粒子数
private static int N = 50;//迭代次数
private static int c1,c2 = 2;
private static double w = 1.4;//惯性因子
private static List<Particle> particles = new ArrayList<Particle>();//粒子群
/**
* 初始化所有粒子
*/
public static void initialParticles() {
for(int i=0;i<particle_num;i++) {
Particle particle = new Particle();
particle.initialX();
particle.initialV();
particle.fitness = particle.calculateFitness();
particles.add(particle);
}
}
/**
* 更新群组最优位置gbest 以及适应值fitness
*/
public static void updateGbest() {
// Double.MAX_VALUE: 最大的double值
double fitness = Double.MAX_VALUE;
int index = 0;
for(int i=0;i<particle_num;i++) {
if(particles.get(i).fitness<fitness) {
index = i;
fitness = particles.get(i).fitness;
}
}
if(fitness<gbest_fitness) {
gbest = particles.get(index).pbest.clone();
gbest_fitness = fitness;
}
}
/**
* 更新每个粒子的速度并返回多维度速度存贮列表
* @return
*/
public static List<String> updateV() {
//多维度速度存贮vList
List<String> vList=new ArrayList<>();
for(Particle particle:particles) {
String vString="";
for(int i=0;i<particle.dimension;i++) {
double v = w*particle.V[i]+c1*rand()*(particle.pbest[i]-particle.X[i])+c2*rand()*(gbest[i]-particle.X[i]);
/*//个人惯性
double inertia=w*particle.V[i];
//个人认知
double personal=c1*rand()*(particle.pbest[i]-particle.X[i]);
//社会认知
double society=c2*rand()*(gbest[i]-particle.X[i]);
//速度更新
double v=inertia+personal+society;*/
if(v>particle.Vmax)
v = particle.Vmax;
else if(v<-particle.Vmax)
v = -particle.Vmax;
particle.V[i] = v;//更新Vi
//拼接维度速度
vString=vString+v+",";
}
vList.add(vString);
}
return vList;
}
/**
* 更新每个粒子的位置和pbest
*/
public static void updateX() {
for(Particle particle:particles) {
for(int i=0;i<particle.dimension;i++) {
//根据速度更新粒子所在位置
particle.X[i] = particle.X[i] + particle.V[i];
}
//新的适应值
double newFitness = particle.calculateFitness();
//如果新的适应值比原来的小则更新fitness和pbest
if(newFitness<particle.fitness) {
particle.pbest = particle.X.clone();
particle.fitness = newFitness;
}
}
}
/**
* 算法主要流程
*/
public static void process() {
int n = 0;
initialParticles();
updateGbest();
while(n++<N) {
List<String> vList=updateV();
updateX();
updateGbest();
System.out.println(n+".当前gbest:("+gbest[0]+","+gbest[1]+") 更新后适应值fitness="+gbest_fitness);
for(int i=0;i<particle_num;i++){
String vString=vList.get(i);
//分割各维度速度
String[] vArray=vString.split(",");
for(int j=0;j<vArray.length;j++){
System.out.println("第 "+(i+1)+" 粒子第 "+(j+1)+" 维速度为："+vArray[j]);
}
System.out.println();
}
System.out.println("-------------------------------------------------------------------------------");
}
}
/**
* 返回一个0~1的随机数
* @return
*/
public static double rand() {
return new Random().nextDouble();
}
public static void main(String[] args) {
process();
}
}
注：以上代码源码思路来自于https://www.cnblogs.com/mrpod2g/p/4575185.html
3.基准优化函数
以上代码所用的测试函数有Ackley、Sphere和Rosenbrock，都属于无约束基准函数，有兴趣了解其他测试函数的小伙伴可以访问这个网址：http://www.sfu.ca/~ssurjano/optimization.html
Ackley
Ackley函数广泛用于测试优化算法。如上图所示，其二维形式的特征是外部区域几乎平坦，中心处有一个大孔。该函数存在使优化算法（尤其是爬坡算法）陷入其许多局部最小值中的风险。
Sphere
球形函数具有d局部最小值，全局变量除外。它是连续的，凸的和单峰的。该图显示了其二维形式。
Rosenbrock
Rosenbrock函数，也称为Valley或Banana函数，是基于梯度的优化算法的常见测试问题。在上面的图中以二维形式显示了它。
该函数是单峰函数，全局最小值位于狭窄的抛物线形山谷中。但是，即使容易找到该山谷，也很难将其收敛到最小（Picheny等，2012）。
参考文献
1.《量子行为粒子群优化：原理及其应用》 作者：孙俊、方伟、吴小俊、须文波
2.《进化优化算法——基于仿生和种群的计算机智能方法》 作者：【美】丹.西蒙 翻译：陈曦