Edison.X. Blog

[回目錄]

註1 : 由於 Online High Light 會把我原本的縮排打亂，本文後面不再特別包 code。

註2：這篇文章前面看起來很亂，建議從頭到尾看兩遍，才會對整體設計架構有所了解。

註3：程式碼很短、說明本身沒什麼，但請尊重原創，請別整份 copy 就走。

註4：IE 無法正常進行顯示或留言者，請換 firefox 瀏覽。

要實作 pso alogrithm 時，通常第一個程式，我習慣放上比較無腦的程式。

假設一個函式為 | x³ - 0.8x² - 10000x + 8000 |， x 範圍定在 [-100 , +100] 之間，

可能是小數，要找最大值。一開始先用暴力破解法，去得到它的解落在哪，C 語言暴力去做大概長這樣，

#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

接下來再以 C 語言開始設計整個 pso 流程，去筆對暴力法找出來的解，看設計的好不好。聲明，這份 C language 是寫給初學者看的，避免在 pso 不明、C 語言不熟，馬上看一些難以理解的 code ，故一些變數刻意安排成 global variable，且也沒進行函式分開編譯，一定要強調，它不是個好的佈局，(但它可以動)，若對 C language 很不熟又要馬上 coding pso，建議這份 code 不懂的地方拿去問別人問清楚 (用留言問語法實在是很沒效率，用留言問演算法或意見可能還好點)。

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] 亂數         */

整份完整程式碼如下，註解已盡可能完整，這份程式碼扣掉註解，實際上只有 150 行左右，即使不會 C 語言，但有了上一篇之介紹 後，花時間看懂它 (程式碼語法不會，記得去問人)，最多應不花超過一星期。

/*******************************************************************/
/*                                                                 */
/*     filename : simple_pso.c                                     */
/*     author   : edison.shih/edisonx                              */
/*     compiler : Visual C++ 2008                                  */
/*     date     : 2011.03.07                                       */
/*                                                                 */
/*         A.L.L.  R.I.G.H.T.S.  R.E.S.E.R.V.E.                    */
/*                                                                 */
/*                                                                 */
/*******************************************************************/


/************************************************************************/
/* 重要變數與參數說明                                                   */
/*                                                                      */
/* 1.particle_cnt: 粒子數, 通常於-40,複雜設-200                         */
/*                                                                      */
/* 2.粒子長度: 即解空間維度,或使用變數個數,由於此問題之適應函式只用到   */
/*             x 變數, 故粒子長度為，將再補程式碼做為考慮粒子長度為n    */
/*             之情形                                                   */
/*                                                                      */
/* 3.max_pos, min_pos: 粒子範圍,即解空間中,各維度(變數)之範圍限制,      */
/*             普遍性應考慮n 維之情形                                   */
/*                                                                      */
/* 4.max_v : 即最大速度限制, 通常設定成粒子之寬度,即解空間範圍,普遍性考 */
/*           慮n 維情形,即每維之解空間範圍                              */
/* 5.c1,c2 : 學習常數, c1,c2多設, 一般c1=c2=[0,4]                       */
/*                                                                      */
/* 6.w     : 慣性權重, 此常數為後來由Yuhui Shi and Russell Eberhart 提  */
/*           出,多設在[0,1.5] 之間, 設1.0 為保守值                      */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>

//////////////////////////////////////////////////////////////////////////
// -----------------------------------------------------------------------
// fitness value
// best value on [-100, 100] , aboue fit(-57.469) = 390245.791738

double fit(double x)
{
    
     // x**3 - 0.8x**2 - 1000x + 8000
     return fabs(8000.0 + x*(-10000.0+x*(-0.8+x)));
}


//////////////////////////////////////////////////////////////////////////
// -----------------------------------------------------------------------
// * 使用pso 求x**3 - 0.8x**2 - 10000x + 8000 於[-100,100] 間之最大值    *
// -----------------------------------------------------------------------
//////////////////////////////////////////////////////////////////////////

/* 定義結構體particle */
typedef struct tag_particle{
     double position;  /* 目前位置, 即x value    */
     double velocity;  /* 目前粒子速度           */
     double fitness ;  /* 適應函式值             */
     double pbest_pos; /* particle 目前最好位置  */
     double pbest_fit; /* particle 目前最佳適應值*/
}particle;

// -----------------------------------------------------------------------
// * 全域變數宣告, 較好佈局為拆pos.h / pos.c , 放在pos.c 宣告成static    *
// * 再做set_param function 為入口, 設定所有static param variable        *
// * 或直接用#define 方式, 也需get_gbest 做入口取得全域最佳值            *
// -----------------------------------------------------------------------

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

// -----------------------------------------------------------------------
// * pso 相關函式宣告                                                    *
// -----------------------------------------------------------------------

#define RND() ((double)rand()/RAND_MAX) /* 產生[0,1] 亂數         */
particle* AllocateParticle();           /* 配置particle_cnt 個粒子*/
void ParticleInit(particle *p);         /* 粒子初始化             */
void ParticleMove(particle *p);         /* 開始移動               */
void ParticleRelease(particle* p);      /* 釋放粒子記憶體         */
void ParticleDisplay(particle* p);      /* 顯示所有粒子資訊       */

//////////////////////////////////////////////////////////////////////////

int main()
{
     /* 變數宣告*/
     unsigned i;
     unsigned max_itera=100000;           /* max_itera : 最大演化代數*/
     particle* p;                         /* p         : 粒子群      */

     /* 設定參數*/
     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;    /* 設最大速限               */

     /* 開始進行*/
     p = AllocateParticle();      // 配置記憶體
     ParticleInit(p);             // 粒子初始化
     for(i=0; i<max_itera; i++)   // 進行迭代
           ParticleMove(p);       // 粒子移動
     ParticleDisplay(p);          // 顯示最後結果
     ParticleRelease(p);          // 釋放記憶體

     // 暴力取得之較佳值
     printf("know2: %10.6lf , %lf\n", -57.469, fit(-57.469));
     return 0;
}

//////////////////////////////////////////////////////////////////////////

// 配置particle_cnt 個粒子
particle* AllocateParticle()
{
     return (particle*)malloc(sizeof(particle)*particle_cnt);
}

// 粒子初始化
void ParticleInit(particle *p)
{
     unsigned i;
     const double pos_range = max_pos - min_pos; // 解寬度
     srand((unsigned)time(NULL));

     // 以下程式碼效率不佳, 但較易懂一點
     for(i=0; i<particle_cnt; i++) {
           // 隨機取得粒子位置, 並設為該粒子目前最佳適應值
           p[i].pbest_pos = p[i].position = RND() * pos_range + min_pos;
           // 隨機取得粒子速度
           p[i].velocity = RND() * max_v;
           // 計算該粒子適應值, 並設為該粒子目前最佳適應值
           p[i].pbest_fit = p[i].fitness = fit(p[i].position);

           // 全域最佳設定
           if(i==0 || p[i].pbest_fit > gbest.fitness)
                memcpy((void*)&gbest, (void*)&p[i], sizeof(particle));
     }
}

// 開始移動
void ParticleMove(particle *p)
{
     unsigned i;
     double v, pos;     // 暫存每個粒子之速度, 位置用
     double ppos, gpos; // 暫存區域及全域最佳位置用
     gpos = gbest.position;

     // 更新速度與位置
     for(i=0; i<particle_cnt; i++){
           v = p[i].velocity; // 粒子目前速度
           pos=p[i].position; // 粒子目前位置
           ppos=p[i].pbest_pos; // 粒子目前曾到到最好位置
          
           v = w*v + c1*RND()*(ppos-pos) + c2*RND()*(gpos-pos); // 更新速度
           if(v<-max_v) v=-max_v;    // 限制最大速度
           else if(v>max_v) v=max_v; // 限制最大速度
          
           pos = pos + v;               // 更新位置
           if(pos>max_pos) pos=max_pos; // 限制最大位置
           else if(pos<min_pos) pos=min_pos; // 限制最小位置

           p[i].velocity= v;        // 更新粒子速度     
           p[i].position=pos;       // 更新粒子位置
           p[i].fitness = fit(pos); // 更新粒子適應值

           // 更新該粒子目前找過之最佳值
           if(p[i].fitness > p[i].pbest_fit) {
                p[i].pbest_fit = p[i].fitness ;
                p[i].pbest_pos = p[i].position;
           }

           // 更新全域最佳值
           if(p[i].fitness > gbest.fitness)
                memcpy((void*)&gbest, (void*)&p[i], sizeof(particle));
     }
}

// 釋放粒子記憶體
void ParticleRelease(particle* p)
{
     free(p);
}

// 顯示所有粒子資訊
void ParticleDisplay(particle* p)
{
     unsigned i;

     /* 若想看完整的粒子資料，可把下面三行註解拿掉，這裡只顯示最佳解。*/
//   for(i=0; i<particle_cnt; i++)
//         printf("#%d : %lf , %lf . %lf\n", i+1, p[i].position, p[i].fitness, p[i].velocity);
//   puts("------------------------------\n");
     printf("best : %10.6lf , %lf\n", gbest.position, gbest.fitness);    
}

單次參考執行結果

best : -57.469015 , 3.9024579173833411e+005

pso 每次找出來的結果可能都不大一樣，但這問題多測幾次，值是差不了多少的

(廢話，因這問題太容易了) 。

而上面用暴力法每次增 5E-4，得到的結果是

x = -57.469000 , 答案為 3.9024579173849098e+005

明顯，若用暴力法要得到像 PSO 那種精度的解答， step 至少要設 5E-6，

到時執行時間，會比 pso 還要久許多。

若有興趣，可再調那些 w, c1, c2, particle_cnt,  max_v , max_itera 試試結果，

通常一些較偏激的值，得到的結果會讓人意外 (可能是較佳，也可能是較差)。

[回目錄]