C語言實現粒子羣算法(PSO)一
最近在溫習C語言，看的書是《C primer Plus》，突然想起來之前在參加數學建模的時候，用過的一些智能算法，好比遺傳算法、粒子羣算法、蟻羣算法等等。當時是使用MATLAB來實現的，並且有些MATLAB自帶了工具箱，當時有些只是利用工具箱求最優解問題，沒有本身動手親自去實現一遍，如今都忘的差很少了。我以爲那樣層次實在是很淺，沒有真正理解算法的核心思想。本着「紙上得來終覺淺，絕知此事要躬行」的態度，我決定如今從新複習一遍算法，而後手工用C語言從新實現一遍。說作就作，我第一個實現的算法是相對來講比較簡單的粒子羣算法(與遺傳算法等相比，至少我本身以爲實現要簡單一些)。html
首先簡單介紹一下啓發式算法和智能算法。粒子羣算法、遺傳算法等都是從傳統的搜索算法演變而來的啓發式算法。啓發式算法（heuristic algorithm)是相對於最優化算法提出的。一個問題的最優算法求得該問題每一個實例的最優解。 啓發式算法能夠這樣定義：一個基於直觀或經驗構造的算法，在可接受的花費（指計算時間和空間）下給出待解決組合優化問題每個實例的一個可行解，該可行解與最優解的偏離程度通常不能被預計，可是一般狀況下啓發式算法能夠給出接近最優解的不錯的解，可是沒法保證每次它均可以獲得很好的近似解。 啓發式算法中有一類被稱之爲智能算法，所謂"智能"二字，指的是這種算法是經過模仿大天然中的某種生物或者模擬某種現象而抽象獲得的算法，好比遺傳算法就是模擬天然界生物天然選擇，優勝劣汰，適者生存而獲得的進化算法，粒子羣是源於對於鳥類捕食行爲的研究，而模擬退火算法則是根據物理學中固體物質的退火過程抽象獲得的優化算法。智能算法興起於上個世紀80年代左右，以後就一直髮展迅速，除了傳統的智能算法以外，近幾年又涌現出了一些新的算法好比魚羣算法、蜂羣算法等。算法
言歸正傳，下面來介紹今天的主角：粒子羣算法。粒子羣算法的基本原理以下（參考《MATLAB智能算法30個案例分析》）：數組
假設在一個D維的搜索空間中，由n個粒子組成的種羣X=(X1,X2,..,Xn),其中第i個粒子表示爲一個D維的向量Xi=(xi1,xi2,xiD),表明第i個粒子在D維搜索空間中的位置，亦表明問題的一個潛在解。根據目標函數便可以計算出每一個粒子位置Xi對應的適應度值。第i個粒子的速度爲Vi = (Vi1,Vi2,...,ViD),其個體極值爲Pi=(Pi1,Pi2,...,PiD),種羣的羣體極值爲Pg=(Pg1,Pg2,...,PgD)。在每次迭代的過程當中，粒子經過個體極值和羣體極值更新自身的速度和位置，即：函數
Vid(k+1)=w*Vid(k)+c1*r1*(Pid(k)-Xid(k))+c2*r2*(Pgd(k)-Xid(k))
Xid(k+1) = Xid(k) + Vid(k+1)
其中w爲慣性權重，若是不考慮能夠默認爲1,後面還會再詳細討論w對於PSO的影響。d=1,2,..,D;i=1,2,...,n;k爲當前迭代次數;Vid爲粒子的速度;c1和c2是非負的常數，稱爲加速度因子;r1和r2是分佈於[0,1]之間的隨機數。爲了防止粒子的盲目搜索，通常建議將其位置和速度限制在必定的區間內。工具
下面是我用C語言實現的求一個二元函數最大值的粒子羣算法：優化
/*
* 使用C語言實現粒子羣算法(PSO)
* 參考自《MATLAB智能算法30個案例分析》
* update: 16/12/3
* 本例的尋優非線性函數爲
* f(x,y) = sin(sqrt(x^2+y^2))/(sqrt(x^2+y^2)) + exp((cos(2*PI*x)+cos(2*PI*y))/2) - 2.71289
* 該函數有不少局部極大值點，而極限位置爲(0,0),在(0,0)附近取得極大值
*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#define c1 1.49445 //加速度因子通常是根據大量實驗所得
#define c2 1.49445
#define maxgen 300 // 迭代次數
#define sizepop 20 // 種羣規模
#define popmax 2 // 個體最大取值
#define popmin -2 // 個體最小取值
#define Vmax 0.5 // 速度最大值
#define Vmin -0.5 //速度最小值
#define dim 2 // 粒子的維數
#define PI 3.1415926 //圓周率
double pop[sizepop][dim]; // 定義種羣數組
double V[sizepop][dim]; // 定義種羣速度數組
double fitness[sizepop]; // 定義種羣的適應度數組
double result[maxgen]; //定義存放每次迭代種羣最優值的數組
double pbest[sizepop][dim]; // 個體極值的位置
double gbest[dim]; //羣體極值的位置
double fitnesspbest[sizepop]; //個體極值適應度的值
double fitnessgbest; // 羣體極值適應度值
double genbest[maxgen][dim]; //每一代最優值取值粒子
//適應度函數
double func(double * arr)
{
double x = *arr; //x 的值
double y = *(arr+1); //y的值
double fitness = sin(sqrt(x*x+y*y))/(sqrt(x*x+y*y)) + exp((cos(2*PI*x)+cos(2*PI*y))/2) - 2.71289;
return fitness;
}
// 種羣初始化
void pop_init(void)
{
for(int i=0;i<sizepop;i++) //sizepop 種羣規模
{
for(int j=0;j<dim;j++)//dim 粒子的維數
{
pop[i][j] = (((double)rand())/RAND_MAX-0.5)*4; //-2到2之間的隨機數 //pop[i][j] 種羣數組
V[i][j] = ((double)rand())/RAND_MAX-0.5; //-0.5到0.5之間 // V[i][j] 種羣速度數組
}
fitness[i] = func(pop[i]); //計算適應度函數值
}
}
// max()函數定義
double * max(double * fit,int size)//尋找數組內的最大值和最大值所處的位置
{
int index = 0; // 初始化序號
double max = *fit; // 初始化最大值爲數組第一個元素
static double best_fit_index[2];
for(int i=1;i<size;i++)
{
if(*(fit+i) > max)
{
max = *(fit+i);
index = i;
}
}
best_fit_index[0] = index;
best_fit_index[1] = max;
return best_fit_index;
}
// 迭代尋優
void PSO_func(void)
{
pop_init();
double * best_fit_index; // 用於存放羣體極值和其位置(序號)
best_fit_index = max(fitness,sizepop); //求羣體極值
int index = (int)(*best_fit_index);
// 羣體極值位置
for(int i=0;i<dim;i++)
{
gbest[i] = pop[index][i];
}
// 個體極值位置
for(int i=0;i<sizepop;i++)
{
for(int j=0;j<dim;j++)
{
pbest[i][j] = pop[i][j];
}
}
// 個體極值適應度值
for(int i=0;i<sizepop;i++)
{
fitnesspbest[i] = fitness[i];
}
//羣體極值適應度值
double bestfitness = *(best_fit_index+1);
fitnessgbest = bestfitness;
//迭代尋優
for(int i=0;i<maxgen;i++)
{
for(int j=0;j<sizepop;j++)
{
//速度更新及粒子更新
for(int k=0;k<dim;k++)
{
// 速度更新
double rand1 = (double)rand()/RAND_MAX; //0到1之間的隨機數
double rand2 = (double)rand()/RAND_MAX;
V[j][k] = V[j][k] + c1*rand1*(pbest[j][k]-pop[j][k]) + c2*rand2*(gbest[k]-pop[j][k]);
if(V[j][k] > Vmax)
V[j][k] = Vmax;
if(V[j][k] < Vmin)
V[j][k] = Vmin;
// 粒子更新
pop[j][k] = pop[j][k] + V[j][k];
if(pop[j][k] > popmax)
pop[j][k] = popmax;
if(pop[j][k] < popmin)
pop[j][k] = popmin;
}
fitness[j] = func(pop[j]); //新粒子的適應度值
}
for(int j=0;j<sizepop;j++)
{
// 個體極值更新
if(fitness[j] > fitnesspbest[j])
{
for(int k=0;k<dim;k++)
{
pbest[j][k] = pop[j][k];
}
fitnesspbest[j] = fitness[j];
}
// 羣體極值更新
if(fitness[j] > fitnessgbest)
{
for(int k=0;k<dim;k++)
gbest[k] = pop[j][k];
fitnessgbest = fitness[j];
}
}
for(int k=0;k<dim;k++)
{
genbest[i][k] = gbest[k]; // 每一代最優值取值粒子位置記錄
}
result[i] = fitnessgbest; // 每代的最優值記錄到數組
}
}
// 主函數
int main(void)
{
clock_t start,finish; //程序開始和結束時間
start = clock(); //開始計時
srand((unsigned)time(NULL)); // 初始化隨機數種子
PSO_func();
double * best_arr;
best_arr = max(result,maxgen); //result 存放每次迭代種羣最優值的數組 //maxgen 迭代次數
int best_gen_number = *best_arr; // 最優值所處的代數
double best = *(best_arr+1); //最優值
printf("迭代了%d次，在第%d次取到最優值，最優值爲:%lf.\n",maxgen,best_gen_number+1,best);
printf("取到最優值的位置爲(%lf,%lf).\n",genbest[best_gen_number][0],genbest[best_gen_number][1]);
finish = clock(); //結束時間
double duration = (double)(finish - start)/CLOCKS_PER_SEC; // 程序運行時間
printf("程序運行耗時:%lf\n",duration);
return 0;
}
我運行C採用的是Ubuntu16 下的gcc編譯器，運行結果截圖以下：spa
屢次運行結果差很少，基本每次均可以很接近最優解。並且發現C語言運行時間要遠快於MATLAB實現（我記得MATLAB要用好幾秒，這裏就不貼MATLAB代碼進行運行時間對比了），只須要耗時0.004秒左右。這裏只討論了基本的粒子羣算法，後面一篇我還會對於粒子羣的參數w進行詳細的討論，討論不一樣的w參數的取法對於粒子羣尋優能力的影響。code
轉自：http://www.cnblogs.com/lyrichu/p/6151272.htmlhtm
- 1. C語言實現粒子羣算法(PSO)一
- 2. C語言實現粒子羣算法（PSO）二
- 3. 粒子羣算法（PSO）
- 4. 優化算法——粒子羣算法(PSO)
- 5. 粒子羣優化算法PSO
- 6. 粒子羣優化（PSO）算法
- 7. 標準粒子羣算法（PSO）
- 8. 標準粒子羣優化算法 PSO
- 9. 粒子羣優化算法(PSO)
- 10. 進化算法之粒子羣（PSO）
- 更多相關文章...
- • XSL 語言 - XSLT 教程
- • C# 基本語法 - C#教程
- • Kotlin學習（一）基本語法
- • C# 中 foreach 遍歷的用法
-
每一个你不满意的现在，都有一个你没有努力的曾经。
- 1. C語言實現粒子羣算法(PSO)一
- 2. C語言實現粒子羣算法（PSO）二
- 3. 粒子羣算法（PSO）
- 4. 優化算法——粒子羣算法(PSO)
- 5. 粒子羣優化算法PSO
- 6. 粒子羣優化（PSO）算法
- 7. 標準粒子羣算法（PSO）
- 8. 標準粒子羣優化算法 PSO
- 9. 粒子羣優化算法(PSO)
- 10. 進化算法之粒子羣（PSO）