C语言实现标准PSO算法

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

C语言实现标准PSO算法

2012-11-06 18:48  3999人阅读  评论(12)  收藏  举报

  分类：

语言学习篇：C/C++（42）   智能算法篇：Soft Computing（14） 

版权声明：本文为博主原创文章，未经博主允许不得转载。

简介以及基本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. }