随机微粒群算法.cpp

随机粒子群算法对解答NP问题、复杂问题

/************************************************************************/
/* 该程序用于计算Goldstein-Price函数的最小值							*/
/*  随机微粒群算法源程序						                        */
/************************************************************************/
#include "stdlib.h"
#include "time.h"
#include "stdio.h"
#include "math.h"

#define rdint(i) (rand()%(int)i)
#define rdft() (float)((double)rdint(16384)/(16383.0))
#define rnd(a,b) (rdint((int)(b)-(int)(a)+1) + (int)(a))

#define DIMENSION 2
#define POPSIZE 20

float W=0.0;
float C1=1.8;
float C2=1.8;
float VMAX=2;
float XMIN = -2.0;
float XMAX = 2.0;
float P[DIMENSION];
float PBEST;

struct indi
{
	float number[DIMENSION];
	float best[DIMENSION];
	float bestfitness;
	float fitness;
	float speed[DIMENSION];
}	individual[POPSIZE];	//粒子数目

void initiate();
void calculation(int number);
void globalbest(int number);
void localbest(int number);
void generate(int number);

//这个函数起什么作用呢？
//对第num个粒子....
void generate(int num)
{
	int j;
	for(j=0;j<DIMENSION;j++)
		individual[num].number[j] = rdft()*(XMAX-XMIN) + XMIN;	//修改它的位置
	calculation(num);	//计算什么东西？猜测是对粒子里面的参数进行更新吧！
	individual[num].bestfitness=individual[num].fitness;	//从此可以看出，generate应该调用得比较早，这么简单就将bestfitness得到，应该是初始化吧！
	for (j=0;j<DIMENSION;j++)
		individual[num].best[j]=individual[num].number[j];	//这里同上面一样...证实我的猜想——如果在计算更优解过程中，应该不会这么草率地赋值呢！
}

void initiate()
{
	int i,j;
	for(i=0;i<POPSIZE;i++)
		for(j=0;j<DIMENSION;j++)
			individual[i].number[j] = rdft()*(XMAX-XMIN) + XMIN;
		for(i=0;i<POPSIZE;i++)
			for(j=0;j<DIMENSION;j++)
				individual[i].speed[j]=VMAX*rdft();
			for(i=0;i<POPSIZE;i++)
				for(j=0;j<DIMENSION;j++)
					individual[i].best[j] = individual[i].number[j];
				for(i=0;i<POPSIZE;i++)
					calculation(i);		//先计算，再得出适应值
				for(i=0;i<POPSIZE;i++)
					individual[i].bestfitness=individual[i].fitness;
				globalbest(0);		//在初始的时候就已经有全局最优解了，我现在很想知道：全局最优是以什么样的形式保存的？位置？适应值？
}
void localbest(int number)
{
	int i;
	//注意，这里是希望求得函数的最小值，所以适应值越小越好！
	//所谓的localbest，指的是单个粒子本身不断调整所得到的最优解；
	if(individual[number].bestfitness > individual[number].fitness)
		for(i=0;i<DIMENSION;i++)
			individual[number].best[i] = individual[number].number[i];
		individual[number].bestfitness = individual[number].fitness;
}
void globalbest(int number)
{
	//所谓的globalbest，指的是所有粒子中得到的最优解。
	int i,j;
	float s = 0;
	int flag=0;
	if(number==0)
	{
		s=individual[0].fitness;
		flag=0;
		for(i=1;i<POPSIZE;i++)
			if(individual[i].fitness<s)
			{
				s=individual[i].fitness;
				flag = i;
			}
			for(i=0;i<DIMENSION;i++)
				P[i]=individual[flag].number[i];
			PBEST=individual[flag].fitness;
	}
	else
	{
		for(i=0;i<POPSIZE;i++)
			if(individual[i].bestfitness<PBEST)
			{
				for(j=0;j<DIMENSION;j++)
					P[j]=individual[i].best[j];
				PBEST=individual[i].bestfitness;
			}
	}
}
void calculation(int num)
{
	//计算函数值
	int i;
	float s=0.0,h=0.0;
	s=pow(individual[num].number[0] + individual[num].number[1] + 1.0, 2.0)*(19.0 -14.0*individual[num].number[0]
			+ 3.0*pow(individual[num].number[0],2.0)-14.0*individual[num].number[1] + 6.0*individual[num].number[0]*
			individual[num].number[1] + 3.0*pow(individual[num].number[1],2.0))+1.0;
	h=30.0 + pow(2.0*individual[num].number[0] - 3.0*individual[num].number[1], 2.0)*(18.0-32.0*individual[num].number[0]
			+12.0*pow(individual[num].number[0],2.0) + 48.0*individual[num].number[1] - 36.0*individual[num].number[0]*individual[num].number[1]
			+27.0*pow(individual[num].number[1],2.0));
	individual[num].fitness = s*h;
}
void main()
{
	int i,j,k,t,total=0,flag=-1,fflag;
	float sum=0,temp[DIMENSION],s;
	srand(time(NULL));
	for (j=0;j<50;j++)
	{
		initiate();
		for (i=0;i<500;i++)
		{
			flag = -1;
			//寻找是否存在该最优微粒,yathing:什么是最优、个体微粒如何利用最优来修正自己?
			for(k=0;k<POPSIZE;k++)
			{
				if((individual[k].fitness == individual[k].bestfitness) && (individual[k].fitness == PBEST))
				{
					flag=k;break;	
				}
			}
			if(flag != -1)
				/************************************************************************/
				/* 为什么在这里调用generate?                                            */
				/************************************************************************/
				generate(flag);
			for(k=0;k<POPSIZE;k++)
				if(k!=flag)
				{
					for(t=0;t<DIMENSION;t++)
					{
						individual[k].speed[t] = W*individual[k].speed[t] + C1*rdft()*(individual[k].best[t] - 
											individual[k].number[t])+C2*rdft()*(P[t]-individual[k].number[t]);
						if(individual[k].speed[t] > VMAX)
							individual[k].speed[t] = VMAX;
						individual[k].number[t]=individual[k].number[t] + individual[k].speed[t];
						if(individual[k].number[t] < XMIN)
							individual[k].number[t] = 2*XMIN - individual[k].number[t];
						if(individual[k].number[t]>XMAX)
							individual[k].number[t]=2*XMAX-individual[k].number[t];
					}
					calculation(k);
					localbest(k);
				}
				globalbest(1);
				if ((PBEST-3.0)<0.00001)
				{
					total++;
					sum=sum+i;
					break;
				}	
		}
		printf("%d,%f\t",i,PBEST);
		printf("Sum is %f\n",sum);
		printf("Average is %f\n",sum/total);
	}
}