[数据结构与算法]C语言实现粒子群算法

C语言实现粒子群算法

一、粒子群算法介绍

??粒子群算法是一种进化算法，其思想来源是模仿自然界中的鸟类觅食。

??假设有50只鸟随机出现在一个位置，并且他们有随机的初始速度，假设单位时间内初始速度不变，单位时间后，他们会到达一个新的位置，并且会判断自己这个位置的好坏程度(可以理解成离食物的远近)，其他的鸟儿下一次选择速度的时候会学习在好坏程度上的自己的历史最优位置和所有鸟的历史最优位置，调整自己的速度，然后重复上述过程，直到轮次结束。

??我们可以用这种思想来寻找全局最优解，假设我们要找一个二维函数f(x,y)的最小值，(x,y)∈D，我们可以先生成50个粒子，他们有随机的初始二维坐标，然后随机生成50个粒子的初始速度，然后经过单位时间后，50个粒子到达了新的位置，计算他们在这个位置的函数值，然后把自己的值放到局部最优里，把所有值中最小的值对应的粒子的信息(位置，坐标，值)放到全局最优数组里。

??接下来的每一轮行为都是相似的，利用循环完成：根据到上一轮以前的全局最优和自己的历史最优和自己上一轮的速度调整自己的速度，然后再求值，再放到局部最优数组的第t个位置里，然后对局部最优数组求最小值放到局部最优数组的第t个位置，然后求此轮全局的最小值，然后比较过往的全局最小值和此轮全局最小值，取较小者放到全局最小值数组的第t个位置中。

??最后，当t大于循环轮次的时候，跳出循环，输出全局最小值数组中的第t-1个位置。

??负责速度更新的式子可以这样写。
$$v_i(t)=\alpha v_i(t?1)+\beta (jub{u}_i(t?1)?r_i(t?1))+\gamma (quanju(t?1)?r_i(t?1))$$
??其中，a是表示之前速度对这次速度选择的影响程度(惯性)的参数,b表示历史最优位置对这次速度选择的影响程度的参数,y是表示历史全局最优对这次速度选择影响程度的参数，vi(t)表示i粒子第t轮的速度,jubui(t-1)表示i粒子第t-1轮以前的历史最优位置，quanjui(t-1)表示t到t-1轮为止的全局最优位置,ri(t-1)表示i粒子第t-1轮的位置。

??负责位置更新的式子可以这样写
$$r_i(t)=r_i(t?1)+v_i(t)$$

二、C语言实现粒子群算法

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define EPS 1e-9
typedef struct {
	double x;
	double y;
}Point;
//二维函数，点是二维坐标。

typedef struct {
	Point point;
	double z;
}Abird;
//Abird结构体存放粒子的点和对应函数值的信息

typedef struct {
	Abird a[500];
}zuiyou;
//局部最优则建立一个zuiyou数组，全局最优建立一个zuiyou就可
//这里的500代表会放500个轮次的全局最优进去，所以500是轮次数。

double f(double x, double y)//目标函数
{
	//return 2 * x * x * x + 3 * y * y * y + 4 * x * y * y - 5 * x * x * y + 10;
	return 2 * x * x + 3 * y * y + 4 * x * y;
}


Abird mymin(Abird bird[], int sz)//选择Abird作为返回值目的是返回整个带有完整信息(点，值)的最小值。
{
	Abird ret = bird[0];
	for (int i = 0; i < sz; i++)
	{
		if (ret.z - bird[i].z > EPS && bird[i].point.x - 0.0 > EPS && 5.0 - bird[i].point.x > EPS && bird[i].point.y - 3.0 > EPS && 5.0 - bird[i].point.y > EPS)
        //控制粒子不要飞出界，飞出界的粒子我们就舍弃了
		{
			ret = bird[i];
		}
	}
	return ret;
}





int main()
{
	Abird bird[50] = { 0 };//50只鸟
	zuiyou jubu[50] = { 0 };//jubu[i].a[t]表示i鸟第t轮的局部最优
	zuiyou quanju = { 0 };//quanju.a[t]表示第t轮的全局最优
	Point v[50] = { 0 };//50只鸟速度矢量的创建
	srand((unsigned)time(NULL));//随机种子
	for (int i = 0; i < 50; i++)
	{
		bird[i].point.x = 0.0 + 1.0 * (rand() % RAND_MAX) / RAND_MAX * (5.0 - 0.0);
		bird[i].point.y = 3.0 + 1.0 * (rand() % RAND_MAX) / RAND_MAX * (5.0 - 3.0);
		v[i].x = 0.0 + 0.1 * (rand() % RAND_MAX) / RAND_MAX * (0.1 - 0.0);
		v[i].y = 0.0 + 0.1 * (rand() % RAND_MAX) / RAND_MAX * (0.1 - 0.0);
	}//初始化50个粒子的速度和位置。
//随机生成[x, y)的浮点数，x ＋ 1.0 * ( rand() % RAND_MAX ) / RAND_MAX * ( y - x )
//随机生成[x,y]的浮点数，x ＋ 1.0 * rand() / RAND_MAX * ( y - x )

	for (int i = 0; i < 50; i++)
	{
		bird[i].point.x += v[i].x;
		bird[i].point.y += v[i].y;
		bird[i].z = f(bird[i].point.x, bird[i].point.y);
		jubu[i].a[0] = bird[i];//第0轮的局部最优就是这次的位置
	}
	quanju.a[0] = mymin(bird, 50);
    //第0轮的全局最优
	int t = 1;
	while (t < 500)//进行499轮
	{
		for (int i = 0; i < 50; i++)
		{
            //更新速度
			v[i].x = 0.5 * v[i].x + 0.2 *  (jubu[i].a[t - 1].point.x - bird[i].point.x) + 0.3 *  (quanju.a[t - 1].point.x - bird[i].point.x);
			v[i].y = 0.5 * v[i].y + 0.2 *  (jubu[i].a[t - 1].point.y - bird[i].point.y) + 0.3 *  (quanju.a[t - 1].point.y - bird[i].point.y);
            //更新位置
			bird[i].point.x += v[i].x;
			bird[i].point.y += v[i].y;
			if (bird[i].point.x - 0.0 > EPS && 5.0 - bird[i].point.x > EPS && bird[i].point.y - 3.0 > EPS && 5.0 - bird[i].point.y > EPS)
             //保证经过第t轮后出界的粒子被舍弃
			{
				bird[i].z = f(bird[i].point.x, bird[i].point.y);
			}
			else
			{
				bird[i].z = 100000.0;
			}
			jubu[i].a[t] = bird[i];//先把此次的局部最优储存起来
			jubu[i].a[t] = mymin(jubu[i].a, t + 1);//与历史过往的局部最优比较
		}
		quanju.a[t] = mymin(bird, 50);//求第t轮的全局最优
		if (quanju.a[t].z - quanju.a[t - 1].z > EPS)
		{
			quanju.a[t] = quanju.a[t - 1];
		}//与上一轮的全局最优相比，若比上一轮差则交换。
		t++;
	}

	printf("粒子群算法的解是%lf，对应坐标为(%lf,%lf)", quanju.a[t - 1].z, quanju.a[t - 1].point.x, quanju.a[t - 1].point.y);
	return 0;
}