[回目錄]


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

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

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

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

 

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

假設一個函式為 | x3 - 0.8x2 - 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 試試結果,

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

 

[回目錄]

/*******************************************************************/
/*                                                                 */
/*     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 <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 - 1000x + 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);    
}


edisonx 發表在 痞客邦 PIXNET 留言(3) 人氣()


留言列表 (3)

發表留言
  • Shih-Che Ku
  • 不好意思
    您在最上面的公式部分少打一個零
    原本是1000x
    應該是10000x
    感謝您的詳細解說讓我學會PSO^^
  • 客氣了 :) 也謝謝您的指正, 已修正。

    edisonx 於 2012/07/11 14:44 回覆

  • sky
  • 不好意思 我想請問 step 這參數代表啥意思 是 迭代次數嗎?
  • step 參數在 pso 裡面並沒有用到 (這篇文章之 pso 版本沒用到) ,是「暴力法」使用的參數。
    | x3 - 0.8x2 - 10000x + 8000 |, x 從 -100 開始,每次 x + step,到 +100 為止,求出最大值。

    edisonx 於 2012/11/05 13:01 回覆

  • 訪客
  • 不好意思 還有個問題就是

    particle gbest; /* 全域最佳值 */

    這邊 particle 前面是不是少了型別 如 double? 之類的?
  • 自己回去翻 C 語言,這篇文章在第一點就講了什麼是 particle,程式語言太基礎的東西不講。

    edisonx 於 2012/11/06 13:22 回覆

您尚未登入,將以訪客身份留言。亦可以上方服務帳號登入留言

請輸入暱稱 ( 最多顯示 6 個中文字元 )

請輸入標題 ( 最多顯示 9 個中文字元 )

請輸入內容 ( 最多 140 個中文字元 )

請輸入左方認證碼:

看不懂,換張圖

請輸入驗證碼

【 X 關閉 】

【PIXNET 痞客邦】國外旅遊調查
您是我們挑選到的讀者!

填完問卷將有機會獲得心動好禮哦(注意:關閉此視窗將不再出現)

立即填寫取消