粒子群算法（初步）

本文参照：http://blog.csdn.net/liangyixin19800304/article/details/12361165

给出了一种C++解决的方法，非常简单，希望与大家共同交流。目前我的疑惑是，在粒子群算法中的那些参数应该如何选取？

1.惯性系数采用了线性变化的形式，在一开始注重大范围搜索，w值比较大为0.9，最后注重局部走索，w为0.4。（论文上说0.9到0.4）。

2.通过查阅论文，c1和c2的取值应该在1.5到2之间，而且通常而言，应该一开始c1比较大，后期c2比较大。

3.惯性参数可以采用线性变化，也可以用指数函数，有相关论文的。

4.速度限幅决定了精度大小。

5.第二份代码，尝试用粒子群算法解决多峰问题，并且得到比较好的结果。

#include "iostream"
#include "algorithm"
#include "vector"
#include "random"
#include "iterator"
using namespace std;
#define MAX_COUNT 500
#define MAX_M 0.9
#define MIN_M 0.4
double fitness_function(double x)
{
	return -x*(x - 1);//适应函数
}
class seed
{
public:
	seed() = default;
	seed(seed& s)
	{
		speed = s.speed;
		w = s.w;
		c1 = s.c1;
		c2 = s.c2;
		present = s.present;
	}
	seed& operator=(const seed& s)
	{
		if (this == &s) return *this;
		speed = s.speed;
		w = s.w;
		c1 = s.c1;
		c2 = s.c2;
		present = s.present;
	}	
	double get_speed()
	{
		return speed;
	}

	~seed() {};
	
	double speed = 0.5;//速度
	double w = 0.9;//惯性系数
	double c1 = 1.5;//认知系数
	double c2 = 2;//社会系数
	double present = 1;//当前位置
	double pbest = 0;

};

double find_max(double a, double b)
{
	return a>b ? a : b;
}




int main()
{
	static uniform_real_distribution<double> u(0, 1);//随机数生成范围
	static default_random_engine e;//随机数生成引擎
	static uniform_real_distribution<double> uu(-2, 2);//随机数生成范围
	static default_random_engine ee;//随机数生成引擎
	int seed_num = 7;
	vector<seed> v(seed_num);
	for (int i = 0; i < seed_num; i++)
	{
		v[i].present = uu(ee);
		v[i].pbest = fitness_function(v[i].present);
	}
	double gbest;
	for (int j = 0; j < MAX_COUNT; j++)
	{
		for (int k = 0; k < seed_num-1; k++)
		{
			gbest = find_max(v[k].pbest, v[k + 1].pbest);//找到的是群体中的最大值
			v[k].w = MAX_M - (MAX_M - MIN_M)*j / MAX_COUNT;
		}
		v[seed_num-1].w = MAX_M - (MAX_M - MIN_M)*j / MAX_COUNT;
		for (int i = 0; i < seed_num;i++)
		{
			v[i].pbest = find_max(v[i].pbest, fitness_function(v[i].present));
			v[i].speed = v[i].w * v[i].speed + u(e) * v[i].c1 * (v[i].pbest - v[i].present) + u(e) * v[i].c2 * (gbest - v[i].present);
			if (v[i].speed < -0.5)
				v[i].speed = -0.5;
			if (v[i].speed > 0.5)
				v[i].speed = 0.5;
			v[i].present += v[i].speed;

		}
		for (int i = 0; i < seed_num; i++)
		{
			cout << v[i].present << "    ";
		}
		cout << gbest << endl;

	}

	//cout << v[0].present << "   " << v[1].present << "   " << max << endl;

	while (1);
	return 0;
}

第二个粒子群算法，用于解决多峰情况。

//粒子群算法
#include "iostream"
#include "algorithm"
#include "vector"
#include "random"
#include "iterator"
using namespace std;
#define MAX_COUNT 200
#define MAX_M 0.9
#define MIN_M 0.4
#define VMAX 0.5
#define N_LIMIT 0
#define P_LIMIT 9
double fitness_function(double x)
{
	//return -x*(x - 1);//适应函数
	//return x*x*sin(3 * x)*cos(2 * x);
	return x + 10 * sin(5 * x) + 7 * cos(4 * x);
}
class seed
{
public:
	seed() = default;
	seed(seed& s)
	{
		speed = s.speed;
		w = s.w;
		c1 = s.c1;
		c2 = s.c2;
		present = s.present;
	}
	seed& operator=(const seed& s)
	{
		if (this == &s) return *this;
		speed = s.speed;
		w = s.w;
		c1 = s.c1;
		c2 = s.c2;
		present = s.present;
	}	
	double get_speed()
	{
		return speed;
	}

	~seed() {};
	
	double speed = 0.9;//速度
	double w = MAX_M;//惯性系数
	double c1 = 1.25;//认知系数
	double c2 = 1.8;//社会系数
	double present = 1;//当前位置
	double pbest = 0;

};

double find_max(double a, double b)
{
	return a>b ? a : b;
}




int main()
{
	static uniform_real_distribution<double> u(0, 1);//随机数生成范围
	static default_random_engine e;//随机数生成引擎
	static uniform_real_distribution<double> uu(0, 9);//随机数生成范围
	static default_random_engine ee;//随机数生成引擎
	int seed_num = 20;
	vector<seed> v(seed_num);
	for (int i = 0; i < seed_num; i++)
	{
		v[i].present = uu(ee);
		v[i].pbest = fitness_function(v[i].present);
	}
	double gbest;
	for (int j = 0; j < MAX_COUNT; j++)
	{
		for (int k = 0; k < seed_num-1; k++)
		{
			gbest = find_max(v[k].pbest, v[k + 1].pbest);//找到的是群体中的最大值
			v[k].w = MAX_M - (MAX_M - MIN_M)*j / MAX_COUNT;
		}
		v[seed_num-1].w = MAX_M - (MAX_M - MIN_M)*j / MAX_COUNT;
		for (int i = 0; i < seed_num;i++)
		{
			v[i].pbest = find_max(v[i].pbest, fitness_function(v[i].present));
			v[i].speed = v[i].w * v[i].speed + u(e) * v[i].c1 * (v[i].pbest - v[i].present) + u(e) * v[i].c2 * (gbest - v[i].present);
			if (v[i].speed < -VMAX)
				v[i].speed = -VMAX;
			if (v[i].speed > VMAX)
				v[i].speed = VMAX;
			double tmp = v[i].present + v[i].speed;
			if (tmp > P_LIMIT || tmp < N_LIMIT)
				v[i].present = uu(ee);//如果过大或者过小，就重新生成这个种子。
			else
				v[i].present = tmp;
		}
	}

	cout << gbest << endl;
	while (1);
	return 0;
}