[回目錄]

 

 

在上一篇 [pso] C 語言第一個 pso 程式 - 架構改善 文章裡面,

已提到了關於一些 pso 架構改善之部份,原本程式碼  [pso] C 語言第一個 pso 程式

欲改善之原因主要有三

 

1.  程式碼內只適用少數之變數個數

 

這部份並不完全是改善的重大原因,因在一些問題裡,

實際上 fitness function 不一定用 array 去做紀錄,

甚至可能不要用 array 方式紀錄可能會較佳。

如 double fitness(double x, double y, double z) ;

如 double fitness(int a, double x, int y, double z);

會考慮改用 array 方式做編碼,乃在於不少問題之自變數確實以 array 紀錄為佳,

此處應視問題,或 coder 習慣而定。

 

2. 另建 pso.c / pso.h

 

這是個人習慣問題。一般在寫稍大之程式碼時,較建議將相關 function、作用之程式碼,

分開進行編譯,也就是這次程式碼與 pso 相關,就另建 pso.h / pso.c 之類的程式碼。

若不習慣以此方式進行,可全都擠到 main.c 檔案裡 (較不建議),

但建議還是要分函式寫,免得到時維護不易。

 

3.  去除抽像化與封裝

 

對程式經驗不夠之 coder 而言,維護抽像化與封裝性之程式碼反而會帶來相關之困擾,

同時做完記憶體分配後,所撰之程式碼皆稍嫌麻煩,故此份程式碼將 struct 全拿掉,

全以動態陣列方式進行。

 

4.  fitness function

 

也由於上述架構有所改變,fitness function 也稍換「難一點點」的適應函式,

一般常用來測試 meta-huristic algorithm 之測試函式,為 De Jong  之測試平台二號函式,

圖形與函式如下

 

fit4.png  

 

 

其 fitness function 為 max f(x, y),

但筆者嫌此份 fitness function 實在是測不出什麼東西出來,太弱,

要測試該 algorithm 所撰之品質如何,實有難度,

所以再換「難一點點」的 function fitness,

函式與其圖形繪制如下

 

fit5.png  

 

要找 max f(x, y)。這隻函式圖形畫出來已算是相當嚇人,

不斷的有 local solution,若還能找到 global solution,

代表其品質已不錯。

但這隻函式還不算是最好的鑑別測試函式,

更難、更好鑑別力之測式函式將於他於補充並撰之,程式碼如下展示。

 

pso_math_main.c

 

/*******************************************************************/
/*                                                                 */
/*     filename : pso_math_main.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.                    */
/*                                                                 */
/*                                                                 */
/*******************************************************************/

#include <stdio.h>
#include "pso_math.h"

#define W        1.2
#define C1       1.5
#define C2       1.5
#define DIM      2        /*
維度    */
#define PCNT     50       /* 粒子數  */
#define ITERATOR 100000   /* 迭代次數*/

int main()
{
     unsigned i = ITERATOR;
     double best_fv;                        /*
最佳適應值*/
     double  *min_pos = Allocate1Dim(DIM);  /* 維度最小值*/
     double  *max_pos = Allocate1Dim(DIM);  /* 維度最大值*/
     double *best_pos = Allocate1Dim(DIM);  /* 維度最佳值*/
     double    *max_v = Allocate1Dim(DIM);  /* 速度限制  */

     // ----------------------------------------------------
     //
設定參數
     min_pos[0]=-3.0   , max_pos[0]=12.1;
     min_pos[1]=4.1    , max_pos[1]=5.8;

     max_v[0] = 0.06 * ( max_pos[0] - min_pos[0] );
     max_v[1] = 0.06 * ( max_pos[1] - min_pos[1] ); 

     SetParam(W, C1, C2, DIM, PCNT);
     SetSolRange(min_pos, max_pos, max_v);
    
     // ----------------------------------------------------
     //
開始進行pso algorithm
    
     ParticleInit();                       /*
初始化設定      */
     while(i--)                            /* 開始進行迭代    */
           ParticleSwarm();                /* 粒子移動與更新  */
     GetOptimzation(best_pos, &best_fv);   /* 取得最佳值      */

     for(i=0; i!=DIM; ++i)
           printf("#x%u = %10.7lf\n", i+1, best_pos[i]);
     printf("find best fitness : %.16e\n", best_fv);

     // ----------------------------------------------------
     //
最後釋放記憶體
     ParticleRelease();
     Release1Dim(min_pos);
     Release1Dim(max_pos);
     Release1Dim(best_pos);
     return 0;
}

 

pso_math.h

 

/*******************************************************************/
/*                                                                 */
/*     filename : pso_math.h                                       */
/*     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.                    */
/*                                                                 */
/*                                                                 */
/*******************************************************************/

#ifndef _PSO_MATH_H_
#define _PSO_MATH_H_

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

double fitness(double *x, int dim) ;  /*
適應函式                 */


double** Allocate2Dim(int m, int n);   /* 配置double arr[m][n]     */
void     Release2Dim(double** ptr);    /* 釋放二維記憶體           */
doubleAllocate1Dim(int x);          /* 配置double arr[x]        */
void     Release1Dim(double* ptr);     /* 釋放一維記憶體           */

void   SetParam(                       /* 設置PSO 參數,配置記憶體 */
                 double w,             /* 自身權重參數w            */
                 double c1,            /* 經驗常數c1               */
                 double c2,            /* 整體經驗常數c2           */
                 unsigned dim,         /* 解空間之維度             */  
                 unsigned particle_cnt /*
粒子個數                 */
               );

void   SetSolRange(                    /*
設定解空間限制           */
                    double* min_pos,   /* 解空間最小限制           */
                    double* max_pos,   /* 解空間最大限制           */
                    double* max_v      /* 粒子移動最大速度限制     */
                  );

void ParticleInit();                   /*
粒子群初始化             */
void ParticleSwarm();                  /* 各粒子開始移動           */
void ParticleDisplay();                /* 顯示粒子詳細資訊         */

void GetOptimzation(                   /* 取得最優解               */
                     double* x,        /* 最優之解                 */
                     double* fitv      /* 最優解之適應值           */
                        );           

void ParticleRelease();                /*
粒子釋放記憶體           */

//////////////////////////////////////////////////////////////////////////
#endif

 

pso_math.c

 

/*******************************************************************/
/*                                                                 */
/*     filename : pso_math.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.                    */
/*                                                                 */
/*                                                                 */
/*******************************************************************/

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

//////////////////////////////////////////////////////////////////////////
//  macro
設定

/* 產生[0,1] 亂數*/
#define RND() ((double)rand() / RAND_MAX)
/*
產生[low,up] 亂數*/
#define random(low, up) ( ((up)-(low)) * RND()  + (low) )
/*
定義PI */
#define PI 3.1415926535897932


//////////////////////////////////////////////////////////////////////////
//
靜態變數

static double w;
static double c1;
static double c2;

static unsigned pcnt;          /*
粒子群個數                       */
static unsigned  dim;          /* 解空間維度                       */

static double **particle_pos/* 各粒子目前位置[pcnt][dim]         */
static double **particle_v;    /* 各粒子目前速度[pcnt][dim]         */
static double  *particle_fv;   /* 各粒子目前適應值[pcnt]            */
static double **pbest_pos;     /* 各粒子目前找到最佳位置[pcnt][dim] */
static double  *pbest_fv ;     /* 各粒子目前最佳適應值[pcnt]        */

static double  *gbest_pos;     /* 全域目前找到最佳位置[dim]         */
static double   gbest_fv;      /* 全域目前最佳適應值               */

static double *max_pos;        /* 各維度之最大值[dim]               */
static double *min_pos;        /* 各維度之最小值[dim]               */
static double *max_v  ;        /* 各維度最大速度限制[dim]           */

void ParticleDisplay()
{
     int i, j;
     for(i=0; i!=pcnt; ++i){
           for(j=0; j!=dim; ++j)
                printf("%+6.3lf(%+6.3lf) ", particle_pos[i][j],  particle_v[i][j]);
           printf("[%+10.5lf]\n",particle_fv[i]);
     }
}

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

// ------------------------------------------------------------------
//
適應函式, De Jong, Goldberg
double fitness(double *x, int dim)
{
     return 21.5 + x[0]*sin(4.0*PI*x[0]) + x[1]*sin(20.0*PI*x[1]);
}


// ------------------------------------------------------------------
//
配置二維記憶體double arr[m][n]
double** Allocate2Dim(int m, int n)
{
     int i;
     double *trunk = (double*)malloc(sizeof(double)*m*n);
     double **p = (double**)malloc(sizeof(double*)*m);
     for(i=0; i<m; ++i){
           p[i]= trunk;
           trunk+=n;
     }
     return p;
}

// ------------------------------------------------------------------
//
釋放二維記憶體
void Release2Dim(double **ptr)
{
     free(*ptr), free(ptr);
}

// ------------------------------------------------------------------
//
配置一維記憶體double arr[x]
doubleAllocate1Dim(int x)
{
     return (double*)malloc(sizeof(double)*x);
}
// ------------------------------------------------------------------
//
配置一維記憶體
void Release1Dim(double* ptr)
{
     free(ptr);
}

// ------------------------------------------------------------------
//
設置pso 相關參數
void   SetParam( double _w,             /* 自身權重參數w         */
                 double _c1,            /* 經驗常數c1            */
                 double _c2,            /* 整體經驗常數c2        */
                 unsigned _dim,         /* 解空間之維度          */  
                 unsigned _particle_cnt /*
粒子個數              */
               )
{
     w=_w, c1=_c1, c2=_c2;
     dim=_dim, pcnt=_particle_cnt;

     //
配置相關記憶體空間
     particle_pos = Allocate2Dim(pcnt, dim); /* 各粒子目前位置        */
     particle_v   = Allocate2Dim(pcnt, dim); /* 各粒子目前速度        */
     particle_fv  = Allocate1Dim(pcnt);      /* 各粒子目前適應值      */
     pbest_pos    = Allocate2Dim(pcnt, dim); /* 各粒子目前找到最佳位置*/
     pbest_fv     = Allocate1Dim(pcnt);      /* 各粒子目前最佳適應值  */
     gbest_pos    = Allocate1Dim(dim);       /* 全域目前找到最佳位置  */
     max_pos      = Allocate1Dim(dim);       /* 各維度之最大值        */
     min_pos      = Allocate1Dim(dim);       /* 各維度之最小值        */
     max_v        = Allocate1Dim(dim);       /* 各維度最大速度限制    */
}

// ------------------------------------------------------------------
//
設置pso 解空間相關參數
void   SetSolRange(                    /* 設定解空間限制        */
                    double* _min_pos,   /* 解空間最小限制        */
                    double* _max_pos,   /* 解空間最大限制        */
                    double* _max_v      /* 粒子移動最大速度限制  */
                  )
{
     memcpy( (void*)min_pos, (void*)_min_pos, dim*sizeof(double));
     memcpy( (void*)max_pos, (void*)_max_pos, dim*sizeof(double));
     memcpy( (void*)max_v  , (void*)_max_v  , dim*sizeof(double));
}

// ------------------------------------------------------------------
//
粒子群初始化
void ParticleInit()
{
     int i, j;
     srand((unsigned)time(NULL));

     for(i=0; i!=pcnt; ++i){
           for(j=0; j!=dim; ++j){
                /*
隨機產生各粒子位置*/
                particle_pos[i][j] = random(min_pos[j], max_pos[j]);
                /*
隨機產生各粒子速度*/
                particle_v[i][j]   = random(0.0, max_v[j]);
           }
           /*
計算各粒子適應值*/
           pbest_fv[i]=particle_fv[i] = fitness(particle_pos[i], dim);

          /*
更新全域最佳解*/
           if(i==0 || pbest_fv[i] > gbest_fv) {
                memcpy((void*)gbest_pos, (void*)pbest_pos[i], sizeof(double)*dim);
                gbest_fv = pbest_fv[i];
           }         
     }
}
static int c=0;
// ------------------------------------------------------------------
//
各粒子開始移動
void ParticleSwarm()
{
     int i, j;
     double tmp;
    
     for(i=0; i!=pcnt; ++i){
           for(j=0; j!=dim; ++j){
                /*
更新各粒子速度*/
                particle_v[i][j] = \
                     w*particle_v[i][j] + \
                     c1 * RND() * (pbest_pos[i][j]-particle_pos[i][j]) + \
                     c2 * RND() * (gbest_pos[j]   -particle_pos[i][j]) ;

                /*
判斷速度是否超過上下界*/
                if(particle_v[i][j] > max_v[j]) particle_v[i][j]=max_v[j];
                else if(particle_v[i][j] < -max_v[j]) particle_v[i][j]=-max_v[j];

                /*
更新各粒子位置*/
                particle_pos[i][j] += particle_v[i][j];
                if(particle_pos[i][j] > max_pos[j]) particle_pos[i][j] = max_pos[j];
                else if(particle_pos[i][j] < min_pos[j]) particle_pos[i][j] = min_pos[j];
           }
           /*
更新各粒子適應值*/
           tmp = particle_fv[i] = fitness(particle_pos[i], dim);
           /*
更新各粒子最佳適應值*/
           if(tmp > pbest_fv[i]) {
                pbest_fv[i] = tmp;
                memcpy((void*)pbest_pos[i], (void*)particle_pos[i], sizeof(double)*dim);
           }
           /*
更新全域最佳適應值 */
           if(tmp > gbest_fv) {
                gbest_fv = tmp;
                memcpy((void*)gbest_pos, (void*)particle_pos[i], sizeof(double)*dim);
           }
     }
}


// ------------------------------------------------------------------
//
取得最優解
void GetOptimzation(double *x, double *fitv)
{
     memcpy((void*)x, (void*)gbest_pos, sizeof(double)*dim);
     *fitv = gbest_fv;
}


// ------------------------------------------------------------------
//
粒子釋放記憶體
void ParticleRelease()
{
     Release2Dim(particle_pos);
     Release2Dim(particle_v);
     Release1Dim(particle_fv);
     Release2Dim(pbest_pos);
     Release1Dim(pbest_fv);
     Release1Dim(gbest_pos);
     Release1Dim(max_pos);
     Release1Dim(min_pos);
}                   
//////////////////////////////////////////////////////////////////////////

 

結果與說明

 

上述之 fitness function,其定義域已定在

-3.0 <= x <= 12.1    ,    4.1 <=y <= 5.8 ,

若使用暴力去算 fitness function, step 精度抓到 5E-4 時,得到的解為

 

(11.625500, 5.725000) 3.8850270522582363e+001

 

然以上述之 pso 方式進行,大約有 5%~10% 的次數會離這個解很遠,

若有找到附近時,數值範例效果大概如下 (截取其中一結果)

#x1 = 12.1000000
#x2 =  5.7250443
find best fitness : 3.8732805969533452e+001

 

雖稍嫌比原暴力解差了點,但在執行速度上卻快了些,

然而要找到比原先暴力法更好的解,次數少之又少,

一方面原因在於電腦本身浮點數精度問題 (取亂數會跑掉一些精度),

一方面問題更大 - 這支 pso 程式,雖是以大多數 coder 習慣為主,

但進入演算法的分析後,會發現有些部份做些調整為佳,

下篇文章 [pso] 粒子移動(pso) 參數設定技巧 便針對這些問題做討論。

 

 

[回目錄]

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


留言列表 (1)

發表留言
  • 關於&quot;[pso] C 語言第一個 pso 程式 - 架構改善實作&quot; 的程式問題
  • edisonx版主你好,很感謝你的pso程式提供免費給普羅大眾參考,
    我有一個小問題不知道是我有所誤會。

    1. 在pso_math.c 的ParticleInit() 裡,
    你雖然在初始化設定給予適應值 pbest_fv[i]=particle_fv[i]
    但好像沒有設定particle_pos 給 pbest_pos

    2. 導致粒子開始更新移動void ParticleSwarm() 的第一個迴圈
    pbest_pos 與 gbest_pos 都會是呈0的狀態,
    造成更新速度一開始都容易為負值。

    3. 於是我在ParticleInit()裡/* 計算各粒子適應值*/的下一行
    加入 memcpy((void*)pbest_pos[i], (void*)particle_pos[i], sizeof(double)*dim);
    即可改善。

    如果錯誤煩請指導,很感謝你的程式,我已結合你的程式並修改成功,謝謝。
    dlnx38@hotmail.com
  • pbest_pos 確實是我忘了給,有空時我再修過,謝謝指正。

    edisonx 於 2013/03/04 12:17 回覆

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

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

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

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

請輸入左方認證碼:

看不懂,換張圖

請輸入驗證碼

【 X 關閉 】

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

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

立即填寫取消