粒子群优化算法

粒子群优化算法属于群智能(swarm intelligence)优化算法。群智能分两种，一种是粒群优化，另一种是蚁群优化。

群智能概念

       假设你和你的朋友正在寻宝，每个人有个探测器，这个探测器可以知道宝藏到探测器的距离。你们一群人在找，每个人都可以把信息共享出去，就跟打dota时你可以有你队友的视野，你可以知道其他所有人距离宝藏的距离，这样，你看谁离宝藏最近，就向谁靠近，这样会使你发现宝藏的机会变大，而且，这种方法比你单人找要快的多。

       这是一个群行为(swarm behavior)的简单实例，群中各个体交互作用，使用一个比单一个体更有效的方法求解全局目标。可以把群(swarm)定义为某种交互作用的组织或Agent之结构集合，在群智能计算研究中，群的个体组织包括蚂蚁，白蚁，蜜蜂，黄蜂，鱼群，鸟群等。在这些群体中，个体在结构上是很简单的，而它们的集体行为却可能变得相当复杂。研究人员发现，蚂蚁在鸟巢和食物之间的运输路线，不管一开始多随机，最后蚂蚁总能找到一条最短路径。

粒群优化概念

粒群优化(particle swarm optimization,PSO)算法是一种基于群体搜索的算法，它建立在模拟鸟群社会的基础上。粒群概念的最初含义是通过图形来模拟鸟群优美和不可预测的舞蹈动作，发现鸟群支配同步飞行和以最佳队形突然改变飞行方向并重新编队的能力。这个概念已经被包含在一个简单有效的优化算法中。

在粒群优化中，被称为“粒子”(particle)的个体通过超维搜索空间“流动”。粒子在搜索空间中的位置变化是以个体成功地超过其他个体的社会心理意向为基础的，因此，群中粒子的变化是受其邻近粒子（个体）的经验或知识影响的。一个粒子的搜索行为受到群中其他粒子的搜索行为的影响。由此可见，粒群优化是一种共生合作算法。

算法描述

先通过一个形象的场景来描述一下：5只鸟觅食，每个鸟都知道自己与食物的距离，并将此信息与其他鸟共享。

一开始，5只鸟分散在不同的地方，假设没只鸟每秒钟更新自己的速度和方向，问题是怎么更新呢？

每只鸟记下自己离食物最近的位置，称为pbest，pbest0,pbest1,..分别表示5只鸟的pbest，从这里面选一个gbest，组里最好的。

每过一秒钟，鸟更新自己的速度v(此处为矢量)，

v_new  = v_old + c1*rand()*(pbest-pcurrent) +c2*rand()*(gbest-pcurrent)

c1,c2一般为2，rand()产生0~1的随机数。

然后以这个速度飞行1秒，再次更新，最终离食物越来越近。

 以下伪码摘自百度百科

　程序的伪代码如下

 

　　For each particle

 

　　____Initialize particle

 

　　END

 

　　Do

 

　　____For each particle

 

　　________Calculate fitness value

 

　　________If the fitness value is better than the best fitness value (pBest) in history

 

　　____________set current value as the new pBest

 

　　____End

 

　　____Choose the particle with the best fitness value of all the particles as the gBest

 

　　____For each particle

 

　　________Calculate particle velocity according equation (a)

 

　　________Update particle position according equation (b)

 

　　____End

 

　　While maximum iterations or minimum error criteria is not attained

 

　　在每一维粒子的速度都会被限制在一个最大速度Vmax，如果某一维更新后的速度超过用户设定的Vmax，那么这一维的速度就被限定为Vmax。

程序实例

计算f(x) = x*x - 20x + 100 的最小值及取最小值时的x

#include <iostream>
#include <cmath>
#include <cstdlib>
using namespace std;

#define C1 2
#define C2 2
#define VMAX 5.0
#define MAX_ITERATIONS 100
float rand01()
{
        return (float) (rand()/(double)RAND_MAX);
}
struct particle{
        float current;
        float pbest;
};

float fitness(float x)
{
        return x*x - 20*x + 100;
}

float gbest = 10000;
struct particle p[5];
float v[5] = {0};

void init_particles()
{
        int i;
        for(i = 0; i < 5; i++)
        {
                p[i].current = -2+i;
                p[i].pbest = p[i].current;
        }
}
void find_gbest()
{
        int i;
        for(i = 0; i < 5; i++)
        {
                        if(fitness(gbest) > fitness(p[i].current))
                                gbest = p[i].current;
        }
}
void adjust_v()
{
        int i ;
        for(i = 0; i < 5; i++)
        {
                v[i] = v[i] + C1*rand01()*(p[i].pbest - p[i].current) + C2*rand01()*(gbest - p[i].current);
                if(v[i] > VMAX)
                        v[i] = VMAX;
        }
}
void pso()
{
        int i,iter_num;
        iter_num = 1;
        while(iter_num < MAX_ITERATIONS)
        {
                /*for(i = 0; i < 5; i++)
                {
                        cout <<"p"<<i<<":current "<<p[i].current<<" pbest "<<p[i].pbest<<endl;
                }
                cout <<"gbest:"<<gbest<<endl;
                cout <<endl;
                getchar();*/
                for(i = 0; i < 5; i++)
                {
                        if(fitness(p[i].current) < fitness(p[i].pbest))
                                p[i].pbest = p[i].current;
                }
                find_gbest();
                adjust_v();
                for(i = 0; i < 5; i++)
                        p[i].current += v[i];
                iter_num ++;
        }
}
int main()
{

        init_particles();
        pso();
        printf("After %d iterations,gbest is %f\n",MAX_ITERATIONS,gbest);
        return 0;
}

运行结果

下面是每次迭代后的结果，可以看到，经过5次迭代，结果就很好了。

After 1 iterations
p0:current -2 pbest -2
p1:current -1 pbest -1
p2:current 0 pbest 0
p3:current 1 pbest 1
p4:current 2 pbest 2
gbest:10000


After 2 iterations
p0:current 1.15506 pbest -2
p1:current 3.79064 pbest -1
p2:current 0.790205 pbest 0
p3:current 2.53646 pbest 1
p4:current 2 pbest 2
gbest:2


After 3 iterations
p0:current 6.15506 pbest 1.15506
p1:current 8.58128 pbest 3.79064
p2:current 5.79021 pbest 0.790205
p3:current 5.87216 pbest 2.53646
p4:current 4.17373 pbest 2
gbest:3.79064


After 4 iterations
p0:current 11.1551 pbest 6.15506
p1:current 13.3719 pbest 8.58128
p2:current 10.7902 pbest 5.79021
p3:current 9.79741 pbest 5.87216
p4:current 8.27141 pbest 4.17373
gbest:8.58128


After 5 iterations
p0:current 13.8766 pbest 11.1551
p1:current 10.1764 pbest 8.58128
p2:current 14.7492 pbest 10.7902
p3:current 13.7227 pbest 9.79741
p4:current 13.2714 pbest 8.27141
gbest:9.79741


After 6 iterations
p0:current 8.03327 pbest 11.1551
p1:current 6.98078 pbest 10.1764
p2:current 13.2414 pbest 10.7902
p3:current 4.78856 pbest 9.79741
p4:current 11.6974 pbest 8.27141
gbest:10.1764


After 7 iterations
p0:current 5.84287 pbest 11.1551
p1:current 9.25245 pbest 10.1764
p2:current 5.23059 pbest 10.7902
p3:current -3.28694 pbest 9.79741
p4:current 9.93147 pbest 11.6974
gbest:10.1764