C语言实现标准PSO算法

http://blog.csdn.net/fovwin/article/details/8154123



C语言实现标准PSO算法

  3999人阅读  评论(12)  收藏  举报
  分类:
 

简介以及基本PSO的实现:http://blog.csdn.net/fovwin/article/category/1256709

相对于基本PSO,标准PSO加入了惯性权重系数W,W代表者粒子群对全局空间的搜索能力和局部收敛速度的权衡。

也就是说,若惯性权重W越大,速度更新式子的第一项所占的比重,即速度惯性项,比较大。粒子就会比较不受约束,可以冲来冲去。不受世俗的约束,也就不容易变“俗”(陷入局部极小点出不来);若惯性权重W越小,速度更新式子的后两项所占的比重,即个体历史极值和全局最优解,比较大。粒子被自己的经验和社会的世俗所约束了,听“妈妈”的话,当然人生走的顺畅,但是你一般在你20岁就知道你后辈子的路,即容易陷入局部极小点出不来——过早收敛早熟。

所以一般在一开始一般设置较大的W值,让它充分搜索全局空间(就像小时候不听话,有时候老爸老妈说了这样不行,我们也会自己试了再说,呵呵),到粒子迭代的后半阶段,设置W较小,可以提高收敛速度(后半辈子大部分按经验活着)。

本文的W按照如下设置:

[cpp]  view plain  copy
  1. /*惯性权重系数*/  
  2. #define W_START 1.4  
  3. #define W_END   0.4  
  4. #define W_FUN   (W_START-(W_START-W_END)*pow((double)cur_n/ITE_N,2))  

当然,我们可以有的别的设置公式,但是大体呈现下降的趋势,也可以用模糊控制规则(下下步实现目标)来设置~~~

至于作用效果么,比基本的能更快早到最优解~~~不过木有pp,1.X版再说~~~

相对与Verson 0.0版减少好几个函数,第一版初级阶段,主要为了学习原理,Version 1.0版代码更加简洁,将pso.c和pso.h独立出来,好以后扩展,通过define来设置你的待优化函数(下一步尝试多种诡异的函数,有的还蛮漂亮的),但是还要完善,但是加入了扩展适应度函数的框架~~~将初始化第一次计算P,X和之后的合并为一体~~~

main.c

[cpp]  view plain  copy
  1. #include "pso.h"  
  2. #include "stdio.h"  
  3.   
  4.   
  5. int main(int argc, const char *argv[])  
  6. {  
  7.     cur_n=0;  
  8.     RandInitofSwarm();             //初始化粒子群  
  9.     while(cur_n++ != ITE_N)  
  10.     {  
  11.         UpdatePandGbest();    //更新个体历史最优解P和全局最优解GBest  
  12.         UpdateofVandX();        //速度和位置更新,即飞翔  
  13.     }  
  14.   
  15.     getchar();  
  16.     return 0;  
  17. }  
pso.h

[cpp]  view plain  copy
  1. #ifndef _PSO_H_  
  2. #define _PSO_H_  
  3.   
  4. /*各种适应度函数选择,要用哪个,就设置为1,但只能有一个为1*/  
  5. #define TEST 0    
  6. #define GOLDSTEIN_PRICE 0  
  7. #define SCHAFFER 1  
  8. #define HANSEN  0  
  9. #define NEEDLE  0  
  10.   
  11. #define Dim 2       //粒子维度  
  12. #define PNum 20     //种群规模  
  13. #define ITE_N  50   //最大迭代次数  
  14. int cur_n;          //当前迭代次数  
  15.   
  16. /*惯性权重函数*/  
  17. #define W_START 1.4  
  18. #define W_END   0.4  
  19. #define W_FUN   (W_START-(W_START-W_END)*pow((double)cur_n/ITE_N,2))  
  20. /*个体和种群结构体*/  
  21. struct PARTICLE  
  22. {  
  23.     double X[Dim];  
  24.     double P[Dim];  
  25.     double V[Dim];  
  26.     double Fitness;  
  27. } particle;  
  28.   
  29. struct SWARM  
  30. {  
  31.     struct PARTICLE Particle[PNum];  
  32.     int GBestIndex;  
  33.     double GBest[Dim];  
  34.     double C1;  
  35.     double C2;  
  36.     double Xup[Dim];  
  37.     double Xdown[Dim];  
  38.     double Vmax[Dim];  
  39. } swarm;  
  40. /*是的,只要三个就好了,更好理解点感觉*/  
  41. void RandInitofSwarm(void);  
  42. void UpdateofVandX(void);  
  43. void UpdatePandGbest(void);  
  44.   
  45. #endif  
pso.c

[cpp]  view plain  copy
  1. #include "stdio.h"  
  2. #include "stdlib.h"  
  3. #include "time.h"  
  4. #include "math.h"  
  5. #include "pso.h"  
  6.   
  7.   
  8. //初始化种群  
  9. void RandInitofSwarm(void)  
  10. {  
  11.     int i, j;  
  12.         //学习因子C1,C2  
  13.     swarm.C1 = 2.0;  
  14.     swarm.C2 = 2.0;  
  15.     for(j = 0; j < Dim; j++)  
  16.     {  
  17.         swarm.Xdown[j] = -10;    //搜索空间范围  
  18.         swarm.Xup[j] = 10;  
  19.         swarm.Vmax[j] = 0.1;       //粒子飞翔速度最大值  
  20.     }  
  21.   
  22.     srand((unsigned)time(NULL));  
  23.     for(i = 0; i < PNum; i++)  
  24.     {  
  25.         for(j = 0; j < Dim; j++)  
  26.         {  
  27.             swarm.Particle[i].X[j] = rand() / (double)RAND_MAX * (swarm.Xup[j] - swarm.Xdown[j]) + swarm.Xdown[j];  //-Xdown~Xup  
  28.             swarm.Particle[i].V[j] = rand() / (double)RAND_MAX * swarm.Vmax[j] * 2 - swarm.Vmax[j];                 //-Vmax~Vmax  
  29.         }  
  30.     }  
  31. }  
  32. /*计算适应度函数的适应度,待进一步完善*/  
  33. static double ComputAFitness(double X[])  
  34. {  
  35.     /* 
  36.         OK 
  37.         min-f(0,0)=3 
  38.     */  
  39. #if TEST  
  40.     return X[0]*X[0]+X[1]*X[1]+3;  
  41. #endif  
  42.   
  43.     /* 
  44.          
  45.         min-f(0,-1)=3.x,y-(-2,2) 
  46.     */  
  47. #if GOLDSTEIN_PRICE  
  48.     return (1+pow(X[0]+X[1]+1,2)*(19-14*X[0]+3*pow(X[0],2)-14*X[1]+6*X[0]*X[1]+ 3*pow(X[1],2))*(30+pow((2*X[0]-3*X[1]),2)*\  
  49.         (18-32*X[0]+12*pow(X[0],2)+48*X[1]-36*X[0]*X[1] + 27*pow(X[1],2))));  
  50. #endif  
  51.   
  52.     /* 
  53.          
  54.         min-f(0,0)=-1.x,y-(-4,4) 
  55.     */  
  56. #if SCHAFFER  
  57.     //函数:Schaffer's F6  
  58.   
  59.     return -1*(0.5-(sin(sqrt(X[0]*X[0]+X[1]*X[1]))*\  
  60.         sin(sqrt(X[0]*X[0]+X[1]*X[1]))-0.5)/pow(( 1+0.001*(X[0]*X[0]+X[1]*X[1])),2));  
  61.   
  62. #endif  
  63.   
  64.     /* 
  65.         此函数局部极值点有760个。 
  66.         min-f(x,y)=-176.541793.x,y-(-10,10) 
  67.         (-7.589893,-7.708314),(-7.589893,-1.425128) 
  68.         (-7.5898893,4.858057),(-1.306708,-7.708314) 
  69.         (-1.306707,-1.425128),(-1.306708,4.858057) 
  70.         (4.976478,-7.708314),(4.976478,-1.425128) 
  71.         (4.976478,4.858057) 
  72.     */  
  73. #if HANSEN  
  74.     int i;  
  75.     double temp1=0,temp2=0;  
  76.     double hansenf=0;  
  77.     for (i=1;i <= 5;i++)  
  78.     {  
  79.         temp1+=i*cos((i-1)*X[0]+i);  
  80.         temp2+=i*cos((i-1)*X[1]+i);  
  81.     }  
  82.     hansenf=-1*temp1*temp2;  
  83.   
  84.     return hansenf;  
  85. #endif  
  86.   
  87.     /* 
  88.          
  89.         min-f(0,0)=-3600,x,y-(-5.12,5.12) 
  90.         该问题的全局最优解被最差解包围, 
  91.         局部极值点:(-5.12,5.12),(-5.12,-5.12) 
  92.                     (5.12,-5.12),(5.12,5.12) 
  93.                     f=-2748.78 
  94.     */  
  95. #if NEEDLE  
  96.     return -1*pow((3/(0.05+pow(X[0]-1,2)+pow(X[1]-1,2))),2);  
  97. #endif  
  98.   
  99. }  
  100.   
  101. //update  V and X  
  102. void UpdateofVandX(void)  
  103. {  
  104.     int i, j;  
  105.     srand((unsigned)time(NULL));  
  106.     for(i = 0; i < PNum; i++)  
  107.     {  
  108.         for(j = 0; j < Dim; j++)  
  109.             swarm.Particle[i].V[j] = W_FUN * swarm.Particle[i].V[j] +  
  110.             rand() / (double)RAND_MAX * swarm.C1 * (swarm.Particle[i].P[j] - swarm.Particle[i].X[j]) +  
  111.             rand() / (double)RAND_MAX * swarm.C2 * (swarm.GBest[j] - swarm.Particle[i].X[j]);  
  112.         for(j = 0; j < Dim; j++)  
  113.         {  
  114.             if(swarm.Particle[i].V[j] > swarm.Vmax[j])  
  115.                 swarm.Particle[i].V[j] = swarm.Vmax[j];  
  116.             if(swarm.Particle[i].V[j] < -swarm.Vmax[j])  
  117.                 swarm.Particle[i].V[j] = -swarm.Vmax[j];  
  118.         }  
  119.   
  120.         for(j = 0; j < Dim; j++)  
  121.         {  
  122.             swarm.Particle[i].X[j] += swarm.Particle[i].V[j];  
  123.             if(swarm.Particle[i].X[j] > swarm.Xup[j])  
  124.                 swarm.Particle[i].X[j] = swarm.Xup[j];  
  125.             if(swarm.Particle[i].X[j] < swarm.Xdown[j])  
  126.                 swarm.Particle[i].X[j] = swarm.Xdown[j];  
  127.         }  
  128.     }  
  129. }  
  130.   
  131. /*更新个体极值P和全局极值GBest*/  
  132. void UpdatePandGbest(void)  
  133. {  
  134.     int i, j;  
  135.     //update of P if the X is bigger than current P  
  136.     for (i = 0; i < PNum; i++)  
  137.     {  
  138.         if (swarm.Particle[i].Fitness < ComputAFitness(swarm.Particle[i].P))  
  139.         {  
  140.             for(j = 0; j < Dim; j++)  
  141.             {  
  142.                 swarm.Particle[i].P[j] = swarm.Particle[i].X[j];  
  143.             }  
  144.         }  
  145.     }  
  146.     //update of GBest  
  147.     for(i = 0; i < PNum; i++)  
  148.         if(ComputAFitness(swarm.Particle[i].P) < ComputAFitness(swarm.Particle[swarm.GBestIndex].P))  
  149.             swarm.GBestIndex = i;  
  150.     for(j = 0; j < Dim; j++)  
  151.     {  
  152.         swarm.GBest[j] = swarm.Particle[swarm.GBestIndex].P[j];  
  153.     }  
  154.     printf("The %dth iteraction.\n",cur_n);  
  155.     printf("GBestIndex:%d \n",swarm.GBestIndex );  
  156.     printf("GBest:" );  
  157.     for(j=0;j<Dim;j++)  
  158.     {  
  159.         printf("%.4f ,",swarm.GBest[j]);  
  160.     }  
  161.     printf("\n" );  
  162.     printf("Fitness of GBest: %f \n\n",ComputAFitness(swarm.Particle[swarm.GBestIndex].P));  
  163. }  
  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值