1 /*2  * 使用C语言实现粒子群算法(PSO)3  * 参考自《MATLAB智能算法30个案例分析》4  * update: 16/12/35 * 本例的寻优非线性函数为6  * 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.712897  * 该函数有很多局部极大值点，而极限位置为(0,0),在(0,0)附近取得极大值8  */9 #include<stdio.h>10 #include<stdlib.h>11 #include<math.h>12 #include<time.h>13 #define c1 1.49445 //加速度因子一般是根据大量实验所得14 #define c2 1.4944515 #define maxgen 300  // 迭代次数16 #define sizepop 20 // 种群规模17 #define popmax 2 // 个体最大取值18 #define popmin -2 // 个体最小取值19 #define Vmax 0.5 // 速度最大值20 #define Vmin -0.5 //速度最小值21 #define dim 2 // 粒子的维数22 #define PI 3.1415926 //圆周率23 24 double pop[sizepop][dim]; // 定义种群数组25 double V[sizepop][dim]; // 定义种群速度数组26 double fitness[sizepop]; // 定义种群的适应度数组27 double result[maxgen];  //定义存放每次迭代种群最优值的数组28 double pbest[sizepop][dim];  // 个体极值的位置29 double gbest[dim]; //群体极值的位置30 double fitnesspbest[sizepop]; //个体极值适应度的值31 double fitnessgbest; // 群体极值适应度值32 double genbest[maxgen][dim]; //每一代最优值取值粒子33 34 //适应度函数35 double func(double * arr)36 {37     double x = *arr; //x 的值38     double y = *(arr+1); //y的值39     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;40     return fitness;41 42 }    43 // 种群初始化44 void pop_init(void)45 {46     for(int i=0;i<sizepop;i++)47     {48         for(int j=0;j<dim;j++)49         {50             pop[i][j] = (((double)rand())/RAND_MAX-0.5)*4; //-2到2之间的随机数51             V[i][j] = ((double)rand())/RAND_MAX-0.5; //-0.5到0.5之间52         }53         fitness[i] = func(pop[i]); //计算适应度函数值54     }55 }56 // max()函数定义57 double * max(double * fit,int size)58 {59     int index = 0; // 初始化序号60     double max = *fit; // 初始化最大值为数组第一个元素61     static double best_fit_index[2];62     for(int i=1;i<size;i++)63     {64         if(*(fit+i) > max)65             max = *(fit+i);66             index = i;67     }68     best_fit_index[0] = index;69     best_fit_index[1] = max;70     return best_fit_index;71 72 }73 // 迭代寻优74 void PSO_func(void)75 {76     pop_init();77     double * best_fit_index; // 用于存放群体极值和其位置(序号)78     best_fit_index = max(fitness,sizepop); //求群体极值79     int index = (int)(*best_fit_index);80     // 群体极值位置81     for(int i=0;i<dim;i++)82     {83         gbest[i] = pop[index][i];84     }85     // 个体极值位置86     for(int i=0;i<sizepop;i++)87     {88         for(int j=0;j<dim;j++)89         {90             pbest[i][j] = pop[i][j];91         }92     }93     // 个体极值适应度值94     for(int i=0;i<sizepop;i++)95     {96         fitnesspbest[i] = fitness[i];97     }98     //群体极值适应度值99     double bestfitness = *(best_fit_index+1);
100     fitnessgbest = bestfitness;
101
102     //迭代寻优
103     for(int i=0;i<maxgen;i++)
104     {
105         for(int j=0;j<sizepop;j++)
106         {
107             //速度更新及粒子更新
108             for(int k=0;k<dim;k++)
109             {
110                 // 速度更新
111                 double rand1 = (double)rand()/RAND_MAX; //0到1之间的随机数
112                 double rand2 = (double)rand()/RAND_MAX;
113                 V[j][k] = V[j][k] + c1*rand1*(pbest[j][k]-pop[j][k]) + c2*rand2*(gbest[k]-pop[j][k]);
114                 if(V[j][k] > Vmax)
115                     V[j][k] = Vmax;
116                 if(V[j][k] < Vmin)
117                     V[j][k] = Vmin;
118                 // 粒子更新
119                 pop[j][k] = pop[j][k] + V[j][k];
120                 if(pop[j][k] > popmax)
121                     pop[j][k] = popmax;
122                 if(pop[j][k] < popmin)
123                     pop[j][k] = popmin;
124             }
125             fitness[j] = func(pop[j]); //新粒子的适应度值
126         }
127         for(int j=0;j<sizepop;j++)
128         {
129             // 个体极值更新
130             if(fitness[j] > fitnesspbest[j])
131             {
132                 for(int k=0;k<dim;k++)
133                 {
134                     pbest[j][k] = pop[j][k];
135                 }
136                 fitnesspbest[j] = fitness[j];
137             }
138             // 群体极值更新
139             if(fitness[j] > fitnessgbest)
140             {
141                 for(int k=0;k<dim;k++)
142                     gbest[k] = pop[j][k];
143                 fitnessgbest = fitness[j];
144             }
145         }
146         for(int k=0;k<dim;k++)
147         {
148             genbest[i][k] = gbest[k]; // 每一代最优值取值粒子位置记录
149         }
150         result[i] = fitnessgbest; // 每代的最优值记录到数组
151     }
152 }
153
154 // 主函数
155 int main(void)
156 {
157     clock_t start,finish; //程序开始和结束时间
158     start = clock(); //开始计时
159     srand((unsigned)time(NULL)); // 初始化随机数种子
160     PSO_func();
161     double * best_arr;
162     best_arr = max(result,maxgen);
163     int best_gen_number = *best_arr; // 最优值所处的代数
164     double best = *(best_arr+1); //最优值
165     printf("迭代了%d次，在第%d次取到最优值，最优值为:%lf.\n",maxgen,best_gen_number+1,best);
166     printf("取到最优值的位置为(%lf,%lf).\n",genbest[best_gen_number][0],genbest[best_gen_number][1]);
167     finish = clock(); //结束时间
168     double duration = (double)(finish - start)/CLOCKS_PER_SEC; // 程序运行时间
169     printf("程序运行耗时:%lf\n",duration);
170     return 0;
171 }