[回目錄]

 

 

在上一篇 [pso] 粒子移動(pso) 參數設定技巧  文章裡面,

筆者提到了幾個 pso 參數設定方式,及一般 coder 較少注意之地方,

這次要實做之程式碼,乃根據其一般 coding 時之缺點加以改善,

同時 fitness function 測試函式,強度再加強 (更難找到 global solution),

以下為這份程式碼一些細節更改之部份。

 

1.  更新速度時使用之亂數

 

vid = w * vid + c1 * rnd() * (pbest_xid - xid)  + c2 * rnd() * (gbest_xd-xid)

原本 rnd () 範圍為 [0,1] ,此處改用 [-1, 1],同時定成 macro RND1()

 

2.  逾界時之處理

 

在 simple_pos 與 pso_math 裡面,不論是速度或位置,

當逾界 (比 max 大或比 min 小) 時,是直接設成邊界值,

此處改成,逾界時,隨機再產生一初始值。

 

3. fitness function

 

在上一篇   文章時,筆者曾提到,該 fitness function 並非最佳鑑別,

目前筆者所知較佳鑑別函式如下

fit7.png  

到處都是極值,容易陷入 local solution,要找到 global solution 不易。

而這份 fitness function,並不只適用於二維變數,若要擴到 n 維變數時,

其 fitness function 變為

f(x1, x2,...xn) = sin(x1)*sin(x2)*...*sin(xn) * sqrt ( fabs (x1*x2*...*xn) )

且在 x1, x2, ... xn 之值域範圍沒太限制,故認為適合拿來當作大多

meta-heuristic algorithm 品質測試函式。

 

架構乃不採用 struct 封裝、抽像化方式,此份程式碼參考如下

 

pso_quality_main.c

 

/*******************************************************************/
/*                                                                 */
/*     filename : pso_quality_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_quality.h"

#define W        0.9
#define C1       2.0
#define C2       2.0
#define DIM      2          /*
維度    */
#define PCNT     30         /* 粒子數  */
#define ITERATOR 10000      /* 迭代次數*/

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]=-10.0   , max_pos[0]=10.0;
     min_pos[1]=-10.0   , max_pos[1]=10.0;

     max_v[0] = 0.0006 * ( max_pos[0] - min_pos[0] );
     max_v[1] = 0.0006 * ( 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_quality.h

 

/*******************************************************************/
/*                                                                 */
/*     filename : pso_quality.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_QUALITY_H_
#define _PSO_QUALITY_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_quality.c

 

/*******************************************************************/
/*                                                                 */
/*     filename : pso_quality.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_quality.h"

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

/* 產生[0,1] 亂數*/
// #define RND() ((double) ((rand()<<15)|rand() )/( (RAND_MAX << 15) | RAND_MAX))
#define RND() ((double)rand()/RAND_MAX )
/*
產生[low,up] 亂數*/
#define random(low, up) ( ((up)-(low)) * RND()  + (low) )
/*
產生[-1,1] 亂數*/
#define RND1() (random(-1.0, 1.0))

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

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)
{
     int i;
     double SinMult=1.0, mult=1.0;
     for(i=0; i!=dim; ++i){
           mult*=x[i];
           SinMult*=sin(x[i]);
     }
     return sqrt(fabs(mult))*SinMult;
}


// ------------------------------------------------------------------
//
配置二維記憶體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];
           }         
     }
}

// ------------------------------------------------------------------
//
各粒子開始移動
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 * RND1() * (pbest_pos[i][j]-particle_pos[i][j]) + \
                     c2 * RND1() * (gbest_pos[j]   -particle_pos[i][j]) ;

                /*
判斷速度是否超過上下界*/
                if(particle_v[i][j] > max_v[j]  ||
                     particle_v[i][j] < -max_v[j]) /*
重新產生一速度*/
                     particle_v[i][j] = random(0.0, max_v[j]);

                /*
更新各粒子位置*/
                particle_pos[i][j] += particle_v[i][j];

                /*
超過上下界,重新產生一新位置*/
                if(particle_pos[i][j] > max_pos[j] ||
                     particle_pos[i][j] < min_pos[j]
                     )
                     particle_pos[i][j] = random(min_pos[j], max_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);
}                   
//////////////////////////////////////////////////////////////////////////


結果說明

 

此處值域邊界均設 [-10.0, 10.0],

以暴力法,step = 5E-3 情況下 (用 5E-4 實在跑太久了),

得到最佳之 global solution 如下

 

(7.915000, 7.915000) 7.8855671033511063e+000

 

以這份 pso 跑出來的結果大致如下 (截取自其中十次平均)

 

#x1 =  7.9157735
#x2 =  7.9170704
find best fitness : 7.8855941946765666e+000

 

每次跑出來之結果並不會差太多,速度比暴力法快,

且得到之值較暴力法略優。若不針對上述第一、二點細節做修正,

事實上得到之解常比暴力法來得差,此處可由讀者自行測試驗證。

筆者所提供改善之方案,似乎還是有點成效。

 

 

[回目錄]

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


留言列表 (2)

發表留言
  • edisonx
  • 今天好像把一篇留言,當作廣告誤刪了...
    我不是故意刪您的留言的, 方便的話歡迎再回來留言 @@
  • 悄悄話

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

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

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

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

請輸入左方認證碼:

看不懂,換張圖

請輸入驗證碼