Particle Swarm Optimization (III)

1. 暴力法找出最佳解

假設一個函式為 | x3 - 0.8x2 - 10000x + 8000 |， x 範圍定在 [-100 , +100] 之間，
可能是小數，要找最大值。一開始先用暴力破解法，去得到它的解落在哪。

#include <stdio.h>
#include <math.h>

double fit(double x) { return fabs(8000.0 + x*(-10000.0+x*(-0.8+x)));}

int main()
{
    const double low=-100.0, up=100.0, step=5E-4;
    double x, max=fit(low), max_x=low;
    for(x=low+step; x<=up; x+=step)
        if(fit(x)>max)
            max_x=x, max=fit(x);
    printf("max fitness is : \n x = %.6lf\nf(x) = %.16e\n", max_x, max);
    return 0;
}

執行結果

max fitness is :
x = -57.469000
f(x) = 3.9024579173849098e+005

接下來是利用PSO演算法求最佳解去和上述暴力法結果做比較。

1. 定義結構 particle

回想一下，一個粒子有哪些資訊要存下來的？它要存所在的位置、目前的速度、目前的適應值，還有本身曾找到最佳位置、本身曾找到最佳的適應值。 struct 便根據這些資訊予以定義

/* 定義結構體particle */

typedef struct tag_particle{ 

    double position; /* 目前位置, 即x value */

    double velocity; /* 目前粒子速度 */ 

    double fitness ; /* 適應函式值 */ 

    double pbest_pos; /* particle 目前最好位置 */ 

    double pbest_fit; /* particle 目前最佳適應值 */ 

}particle;

2. 相關參數宣告

權重相關參數 w, c1, c2 ;

最大速度限制 max_v ;

最大位置限制，也就是上面設定 x 的範圍，這裡命名成 max_pos , min_pos；

粒子數量 particle_cnt ;

還有一個，要多設一個粒子，它是紀錄目前所有粒子，找到的最佳解，gbest。

上面這些變數設成全域變數，整體宣告如下

double w, c1, c2; /* 相關權重參數 */
double max_v; /* 最大速度限制 */
double max_pos, min_pos; /* 最大,小位置限制*/
unsigned particle_cnt; /* 粒子數量 */
particle gbest; /* 全域最佳值 */

3. 主程式的長相

基於 pso 演算法的流程，主程式 main 的長相，一開始是去設定可調試的參數值，如下所示。

int main()
{
    /* 設定參數*/
    min_pos = -100.0 , max_pos=+100.0; /* 位置限制, 即解空間限制 */
    w = 1.00, c1=2.0, c2=2.0 ; /* 慣性權重與加速常數設定 */
    particle_cnt = 20; /* 設粒子個數 */
    max_v = (max_pos-min_pos) * 1.0; /* 設最大速限 */

    /* 做其它事*/
    return 0;
}

設完參數值後，後面的整個流程大致就長如下

    動態配置粒子群記憶體
    初始化各個粒子
    for(i=0; i< 最大演化代數; i++) { // 進行迭代
    移動每個粒子,更新速度, 更新位置
    更新最佳之粒子 gbest
    }
    顯示最後結果
    釋放粒子群記憶體

大致只有幾個函數要做

particle* AllocateParticle(); /* 配置particle_cnt 個粒子*/
void ParticleInit(particle *p); /* 粒子初始化 */
void ParticleMove(particle *p); /* 開始移動 */
void ParticleRelease(particle* p); /* 釋放粒子記憶體 */
void ParticleDisplay(particle* p); /* 顯示所有粒子資訊 */

然而在做速度更新時，會不斷地重覆取 [0,1] 之間的亂數，於是定一個 macro 出來

#define RND() ((double)rand()/RAND_MAX) /* 產生[0,1] 亂數 */