粒子羣優化算法

粒子羣優化算法屬於羣智能(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