在上一篇 [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 之測試平台二號函式，
圖形與函式如下
其 fitness function 為 max f(x, y)，
但筆者嫌此份 fitness function 實在是測不出什麼東西出來，太弱，
要測試該 algorithm 所撰之品質如何，實有難度，
所以再換「難一點點」的 function fitness，
函式與其圖形繪制如下
要找 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); /* 釋放二維記憶體 */
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_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]
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];
}
}
}
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) 參數設定技巧 便針對這些問題做討論。
留言列表
開發手札 (2)
{{ article.title }}