在上一篇 [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 並非最佳鑑別，
目前筆者所知較佳鑑別函式如下
到處都是極值，容易陷入 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); /* 釋放二維記憶體 */
double* Allocate1Dim(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]
double* Allocate1Dim(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
每次跑出來之結果並不會差太多，速度比暴力法快，
且得到之值較暴力法略優。若不針對上述第一、二點細節做修正，
事實上得到之解常比暴力法來得差，此處可由讀者自行測試驗證。
筆者所提供改善之方案，似乎還是有點成效。
留言列表
開發手札 (2)
{{ article.title }}