PSO粒子群算法（鸟群算法）计算二元函数极值（C语言实现、matlab工具箱实现）

  PSO算法是经典的智能优化算法，在数学建模等比赛中非常常用，求解时的效果不错。对于智能优化算法，个人倾向于matlab实现，因为计算起来非常方便。但是这次因为老师的要求，准备C语言实现PSO算法求二元函数极值。说实话，第一次听到老师的要求时吓了一跳，因为从来没想过尝试过C语言实现这个算法，感觉这操作多少是有点阴间了。不过静下心来写写发现其实还是很容易的，可能是因为粒子群优化算法比较简单。

算法

下面是一个不太标准的流程图，个人认为认真读代码3分钟就能看懂（可以看C语言代码中的PSO_search）。

  这里浅浅讲一下自己对PSO算法的理解，个人认为PSO算法实际上就是有方向的暴力搜索，一开始瞎搜，然后粒子们逐渐朝表现最好的粒子的方向靠拢着搜索，最后会收敛于某个点。

  粒子群优化算法里面有粒子的加速度和速度的最大最小值，个人认为这个东西有点像深度学习里的学习率，用来控制收敛的速度。大了可能在最优值附近摇摆不定，小了可能陷入局部最优出不来。

实例

  待求极值函数如下：

求解极值 C 语言实现

  C语言代码如下，主要的流程在 PSO_search 函数里写的比较清楚；两个题目可以在 fitting 函数里面切换，修改的时候记得改 #define 里的最大最小值。

#include<bits/stdc++.h>
using namespace std;

// 种群大小 
#define population_size 300
// 维度=> 二元的，二维的 
#define dimension 2 
// 迭代次数
#define max_epoch 15
// 最大值最小值
#define population_max 10
#define population_min 0
// 定义速度
#define velocity_max 5
#define velocity_min -5
// 定义维数
#define dimension 2 
// 网上说取这个很好 
#define c1 1.49445 
#define c2 1.49445


// 鸟群 
double population[population_size][dimension];
// 鸟群速度 
double velocity[population_size][dimension];
// 适应度，其实就是每个粒子对应的解 
double fitness[population_size];

// 个体极值的位置
double pbest[population_size][dimension];  
// 个体极值适应度的值
double fitnesspbest[population_size]; 
// 群体极值的位置
double gbest[dimension];
// 群体极值适应度值
double fitnessgbest; 

// 存每次迭代种群后储存本次迭代最优值的数组=>存的是适应度 
double result[max_epoch];
// 存种群每次迭代后最优值（适应度）时的x,y => 当然不是二维也可以，不一定时x,  y 
double fitnessepoch[population_size][dimension];

// 不让用库，只能手写它俩 
double max(double x1, double x2){
	if(x1>x2)
		return x1;
	else 
		return x2;
}

double min(double x1, double x2){
	if(x1<x2)
		return x1;
	else 
		return x2;
}

// 适应度函数，其实就是目标函数 
double fitting(double x, double y){ 
//	第一题：   
//	double fitness = sin(x)/x+sin(y)/y;
//	第二题 
	double fitness = (6.452*(x+0.125*y)*pow((cos(x)-cos(2*y)),2))/sqrt(0.8+pow((x-4.2),2)+2*pow((y-7),2))+3.226*y;
	return fitness;
}

// 初始化函数，一开始还是随机初始化 
void initialization(void){
	for(int i = 0; i < population_size; i++){
		for(int j = 0; j < dimension; j++){
			// 初始化第一批粒子的位置和速度 
			population[i][j] = ((double)rand()/(RAND_MAX))*(population_max-population_min)+population_min;
			velocity[i][j] = ((double)rand()/(RAND_MAX))*(velocity_max-velocity_min)+velocity_min;
		}
		// 顺便算出来第一批粒子的适应度函数 
		fitness[i] = fitting(population[i][0], population[i][1]);
//		printf("%lf\t%lf\tfitness:%lf\n", population[i][0], population[i][1], fitness[i]);
	}
}

// 作用其实就是比较所有散点对应的适应度，输出最大的适应度和对应的位置 
double* max_in_community(double fit[], int size){
	static double result[2];
	result[0] = *fit;
	result[1] = 0;
	for(int i = 0; i < size; i++)
		if(fit[i]>result[0]){
			result[0] = fit[i];
			result[1] = i;
		}
//	printf("%lf\n%lf\n",result[0], result[1]);
	return result;
}

void PSO_search(void){
//——————————Step1——————————// 
	// 随机初始化种群 
	initialization();
	//先初始化下第一轮的变量 
	double *best_fitness;
	// 该轮所有鸟的最大值和index 
	best_fitness = max_in_community(fitness, population_size); 
	// 存群体的极值在哪里 
	for(int i = 0; i<dimension; i++)
		gbest[i] = population[int(best_fitness[1])][i];
		
//——————————Step2——————————// 
	// 个体极值位置 
	for(int i = 0; i < population_size; i++)
		for(int j = 0; j < dimension; j++)
			pbest[i][j] = population[i][j];
	// 个体极值
	for(int i = 0; i < population_size; i++)
		fitnesspbest[i] = fitness[i]; 
	// 群体极值
//	printf("%lf\n%lf\n", best_fitness[0], best_fitness[1]);
	fitnessgbest = best_fitness[0];

//——————————Step3——————————// 
	// 开始迭代
	for(int i=0; i<max_epoch; i++){
		for(int j = 0; j<population_size; j++){
			for (int k = 0; k<dimension; k++){
				// 更新速度 
				// 这个是PSO的公式
				double r1 = (double)rand()/RAND_MAX;
				double r2 = (double)rand()/RAND_MAX;
				velocity[j][k] =  velocity[j][k] + c1*r1*(pbest[j][k]-population[j][k]) + c2*r2*(gbest[k]-population[j][k]); 
				// 防止速度超出范围 
				velocity[j][k] = min(velocity[j][k], velocity_max);
				velocity[j][k] = max(velocity[j][k], velocity_min);
				// 更新位置
				population[j][k] += velocity[j][k];
				population[j][k] = min(population[j][k], population_max);
				population[j][k] = max(population[j][k], population_min); 
			}
//			printf("(%lf,%lf)\n", population[j][0], population[j][1]);
			// 计算更新后的fitness
			fitness[j] = fitting(population[j][0], population[j][1]); 
		}
		// 对于每一次迭代，都要更新极值(个体 群体)
		for(int j = 0; j<population_size; j++){
			if(fitness[j] > fitnesspbest[j]){
				for(int k=0;k<dimension;k++)
					pbest[j][k] = population[j][k];
				fitnesspbest[j] = fitness[j];
			}
		// 对于每一次迭代，同样也需要更新群体极值
			if(fitness[j] > fitnessgbest){
				for(int k=0;k<dimension;k++)
					gbest[k] = population[j][k];
				fitnessgbest = fitness[j];
			} 
		} 
		// 迭代完了， 记录下每次的结果
		for(int k = 0; k<dimension; k++)
			fitnessepoch[i][k] = gbest[k];
		// 记录一下 
		result[i] = fitnessgbest;
		printf("第%d次迭代：(%lf,%lf)%lf\n", i, fitnessepoch[i][0], fitnessepoch[i][1], result[i]);
	}	
} 

int main(){
	//程序开始和结束时间
	clock_t begin, end;
	begin = clock(); //开始计时 
	srand((unsigned)time(NULL)); // 初始化随机数种子
	// 搜索 
	PSO_search();
	double * solve;
	solve = max_in_community(result,max_epoch);
	printf("迭代了%d次，在第%d次收敛到最优值，最优解为：%lf\n", max_epoch, solve[1], solve[0]);
	// 这次迭代中的最优解中找位置 
	printf("取到最优值的位置为(%lf,%lf).\n", fitnessepoch[int(solve[1])][0], fitnessepoch[int(solve[1])][1]);
	end = clock();
	printf("计算总耗时：%lf秒", end - begin); 
	return 0;
} 

求解极值 matlab 语言实现

画图

  不是很确定最后的结果对不对，拿matlab试试结果。先浅浅画俩图：

	fsurf(@(x,y)sin(x)/x+sin(y)/y)
	fsurf(@(x,y)(6.452*(x+0.125*y)*(cos(x)-cos(2*y))^2/(0.8+(x-4.2)^2+2*(y-7)^2)^(1/2))+3.226*y)

  例题1：

  例题2：

PSO 求解

  matlab自带PSO求解工具计算极值，因为这个东西只能算最小值，所以对原函数取负号，得到的最小值再取负就是我们要的最大值：

	[x,fval] = particleswarm(@(x)-(sin(x(1))/x(1)+sin(x(2))/x(2)),2,[-10,-10],[10,10])
	[x,fval] = particleswarm(@(x)-(6.452*(x(1)+0.125*x(2))*(cos(x(1))-cos(2*x(2)))^2/(0.8+(x(1)-4.2)^2+2*(x(2)-7)^2)^(1/2)+3.226*x(2)),2,[0,0],[10,10])

  运行结果如下：

例题1

例题2

结论

  经过证明，我们C语言的实现是没什么大问题的。