粒子群（PSO）算法的理解与应用

最近在学习粒子群算法，看了很多资料都有点摸不清头脑，直到看了一篇博客中超级简洁的粒子群C++实现代码，才明白粒子群算法的原理，真心感谢博主，在此贴出博主的博客地址：

http://blog.sina.com.cn/s/blog_4ed027020100c9lx.html

在前面的文章中，我们提到的梯度下降法、牛顿法、高斯-牛顿法，以及LM算法等都属于最优化算法，这些最优化算法的共同点是优化一组可能的解，使该解尽可能接近真实解。粒子群算法（也称为PSO算法）也属于最优化算法，但该算法与上述提到算法的主要区别在于，该算法同时优化多组可能的解，最后在多组可能的解中选择最接近真实解的一组作为最终解。每一组可能的解就是一个粒子，因此多组可能解组成了粒子群，这就是该算法名称的由来。下面将讲解我对粒子群算法的理解，并将该算法应用于一个简单例子的最优化求解。

1. 粒子群算法的核心思路

粒子群算法被提出的灵感来源于鸟群觅食。如下图，鸟群觅食过程中，每只鸟沿着各个方向飞行去寻找食物，每只鸟儿都能记住到目前为止自己在飞行过程中最接近食物的位置，同时每只鸟儿之间也有信息共享，它们会比较到目前为止各自与食物之间的最近距离，从各自的最近距离中，选择并记忆整体的一个最近距离位置。由于每只鸟儿都是随机往各个方向飞行的，各自的最近距离位置与整体最近距离位置不断被更新，也即它们记忆中的最近位置越来越接近食物，当它们飞行到达的位置足够多之后，它们记忆的整体最近位置也就达到了食物的位置。

抽象成数学模型，每只鸟儿就是一个粒子，食物的位置也就是问题的最优解，鸟儿与食物的距离也即当前粒子的目标函数值。比如要求z=x*x+y*y的最优解，设置粒子数为6个，如下图所示，相当于各个粒子在曲面上滚动，每个粒子i（0≤i≤5）记住自己滚动过程中（迭代）的最低位置（位置越低z值越小）Zi（0≤i≤5），同时不同粒子之间也有交流，它们会比较Z0~Z5，从中取得一个最小值作为当前的全局最低位置。

2. 粒子群算法的实现

假设目标函数有n个输入参数：f(x1, x2, ..., xn)。

对于每个粒子，它包含的元素：

(1) 当前时刻位置：(x1, x2, ..., xn)

(2) 当前时刻的目标函数值（也称为适应度）：f(x1, x2, ..., xn)

(3) 该粒子的历史最优位置：(x1_pbest, x2_pbest, ..., xn_pbest)

(4) 该粒子的历史最优目标函数值：f_pbest = f(x1_pbest, x2_pbest, ..., xn_pbest)

对于全部粒子，它们共有的元素：

(1) 粒子总数：num

(2) 迭代总次数：cnt。粒子群算法也是一个迭代的过程，需要多次迭代才能获取到理想的最优解。

(3) 全局最优位置：(x1_gbest, x2_gbest, ..., xn_gbest)

(4) 全局最优目标函数值：f_gbest = f(x1_gbest, x2_gbest, ..., xn_gbest)

(5) 位置随机化的上下限：xmin，xmax。迭代开始的时候以及迭代的过程中，均需要对粒子的位置进行随机分布，需要设置随机分布的上下限，不然随机分布偏离得太远，会严重影响优化结果。

(6) 速度的上下限：Vmin，Vmax。迭代过程中，速度也具有一定的随机性，需要限制速度的大小在一定范围内，不然如果速度值太大也会严重影响优化结果。

(7) 速度计算参数：c1、c2。通常取1.0~1.8的值。

粒子位置的更新公式：

参考上述提到的博客，对于每一个粒子的每一个位置参数xi，其位置更新公式如下，其中randf为0~1的随机数，c1与c2通常取1.0~1.8之间的固定值。

为了加快收敛速度，根据速度具有惯性的原理，后来人们提出了在速度计算中增加惯性，也即加上上一轮迭代时该位置参数的速度。如下式，其中w为0~1的参数，通常随着迭代次数的增加，逐渐减小w，因为越到后面可能解就越接近真实解，迭代收敛就越慢，所以需要减小w来放慢速度，否则容易错过最优解。

C++代码实现（本代码参考了上述提到的博客的代码）：

定义全局参数：

#define DATA_SIZE 500
const int NUM = 300;  //粒子总数
const float c1 = 1.65;  //速度计算参数1
const float c2 = 1.65;  //速度计算参数2
const float xmin = -150;   //位置下限
const float xmax = 150;    //位置上限


const float vmin = -250.5;   //速度下限
const float vmax = 250.5;    //速度上线

定义粒子结构体：

struct particle
{
  float x[DATA_SIZE];   //该粒子的当前时刻位置
  float bestx[DATA_SIZE];//该粒子的历史最优位置
  float f;       //该粒子的当前适应度
  float bestf;   //该粒子的历史最优适应度
}swarm[NUM];   //总共有num个粒子

获取0~1的随机数：

#define randf ((rand()%10000+rand()%10000*10000)/100000000.0) //产生0-1随机浮点数

目标函数，显然目标函数的真实解为(0, 1, 2，..., dim-1)：

float f1(float x[], int dim)
{
  float z = 0;
  for (int i = 0; i < dim; i++)
  {
    z += (i+20)*(x[i] - i)*(x[i] - i);
  }
  return z;
}

迭代优化过程：

/*
best为全局最优位置
DIM为目标函数的输入参数总个数，也即每个粒子的位置参数总个数
*/
void PSO(float *best, int DIM)
{
  for (int i = 0; i < DIM; i++)//初始化全局最优
  {
    best[i] = randf*(xmax - xmin) + xmin; //随机生成xmin~xmax的随机数
  }
  //初始化全局最优目标函数值为一个非常大的值
  double gbestf = 100000000000000.0;  
  
  for (int i = 0; i < NUM; i++)  //初始化粒子群
  {
    particle* p1 = &swarm[i];   //获取每一个粒子结构体的地址
    for (int j = 0; j < DIM; j++)
    {
      //随机生成xmin~xmax的随机数作为该粒子的位置参数
      p1->x[j] = randf*(xmax - xmin) + xmin;
    }
    p1->f = f1(p1->x, DIM);   //计算当前时刻的目标函数值
    p1->bestf = 100000000000000.0;  //初始化每个粒子的历史最优目标函数值为一个非常大的值
  }


  float *V = (float *)calloc(DIM, sizeof(float));
  const int cnt = 2000;   //总共2000次迭代
  float w = 0.0025/(cnt-1);  //设置w初值，后面逐渐减小
 
  for (int t = 0; t < cnt; t++)   //cnt次迭代
  {
    for (int i = 0; i < NUM; i++)  //NUM个粒子
    {
      particle* p1 = &swarm[i];
      for (int j = 0; j < DIM; j++)   //计算速度，并更新位置
      {
        //w*(cnt-1-t)随着t的增加逐渐减小
        V[j] = w*(cnt-1-t)*V[j] + c1*randf*(p1->bestx[j] - p1->x[j]) + c2*randf*(best[j] - p1->x[j]);
        //V[j] = c1*randf*(p1->bestx[j] - p1->x[j]) + c2*randf*(best[j] - p1->x[j]);
        //限制速度的上下限
        V[j] = (V[j] < vmin) ? vmin : ((V[j] > vmax) ? vmax : V[j]);
        p1->x[j] = p1->x[j] + V[j];  //更新位置参数
      }


      p1->f = f1(p1->x, DIM);   //计算当前粒子的当前时刻目标函数值
      //如果当前粒子的当前时刻目标函数值小于其历史最优值，则替换该粒子的历史最优
      if (p1->f < p1->bestf)   
      {
        for (int j = 0; j < DIM; j++)
        {
          p1->bestx[j] = p1->x[j];
        }
        p1->bestf = p1->f;
      }


      //如果当前粒子的历史最优值小于全局最优值，则替换全局最优值
      if (p1->bestf < gbestf)   
      {
        for (int j = 0; j < DIM; j++)
        {
          best[j] = p1->bestx[j];
        }
        
        //这里有点不理解，为什么替换全局最优值之后要重新随机分布该粒子的位置参数？
        //如果没有这部分代码，整体优化效果就会很差
        for (int j = 0; j < DIM; j++) 
        {
          p1->x[j] = randf*(xmax - xmin) + xmin;
        }


        gbestf = p1->bestf;
        printf("t = %d, gbestf = %lf\n", t, gbestf);
      }
    }
  }
  
  free(V);
}

主函数：

void main(void)
{
  int DIM = 400;//维数
  float gbestx[DATA_SIZE];//全局最优位置


  PSO(gbestx, DIM);


  printf("全局最优位置:\n");
  for(int i = 0; i < DIM; i++)
  {
    printf("%f ", gbestx[i]);
    if (i > 0 && i % 7 == 0)
      printf("\n");
  }
  printf("\n");
}

运行上述代码，得到结果如下，可知获取的最优解还是非常接近真实解的。

计算速度时，分别加入惯性与不加入惯性，记录每轮迭代的全局最优目标函数值，如下图所示，可以看到加入惯性之后全局最优目标函数值减小得更快。