pso.cpp

粒子群优化算法

//微粒群算法；
#include<iostream.h>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<malloc.h>
#include<time.h>

#define FALSE 0
#define TRUE 1
#define POPSIZE 40 //种群大小
#define MAXGEN 1000  //最大迭代次数
#define MAXRUN 5 //最多运行次数
#define MAXW 1.2 //最大惯性权重
#define MINW 0.1 //最小惯性权重
#define MAXV 5 //最大的速度
#define C1 2 
#define C2 2
#define DIM 2 //维数
#define ERRGOAL 1e-14
#define E 1e-7
#define L 5  //不等式约束数目
#define K 6  //等式约束数目
#define	N L+K //约束数目
#define PI 3.1415926

typedef double (*Tfun)(double pos[DIM]);//定义指向函数的指针类型；
typedef double Particle[DIM];
Tfun g[L];
Tfun h[K];
Particle pos[POPSIZE],pbestpos[POPSIZE];//种群数组，每个个体历史最优位置数组

//变量定义；
int i,j,counter;
int gen,success,gbestpar;
double w,wgen,wdec,lastbest;
double pbestvalue[POPSIZE],parvalue[POPSIZE],weight[POPSIZE];
double upper[DIM]={10.0,10.0},lower[DIM]={-10.0,-10.0};


double vec[POPSIZE][DIM],gbestpos[DIM];

//函数申明；
double fun_randval(double low,double high);//在low-high之间取随机数
void evaluate(int popsize,double pop[][DIM],double parvalue[],double w[]);//评估种群
double fun_object(double *x);//求适应值
double fun_w(double pos[]);//求罚值
int best(double pos[][DIM],double pbestvalue[],int size);//求种群中的最好个体
int better(double w1,double w2,double funcvalue1,double funcvalue2);//比较两个个体的好坏

double g1(double x[DIM]);
double g2(double x[DIM]);
double g3(double x[DIM]);
double g4(double x[DIM]);
double g5(double x[DIM]);

double h1(double x[DIM]);
double h2(double x[DIM]);
double h3(double x[DIM]);
double h4(double x[DIM]);
double h5(double x[DIM]);
double h6(double x[DIM]);



void main()
{
	int curun=0,sumgen=0;
	double sum=0.0;
	double avgbest=0.0;
	double thebestvalue;
	double thebestpos[DIM];
	for(curun=0;curun<MAXRUN;curun++)
	{

//初始化;
	g[0]=g1;
	g[1]=g2;
	g[2]=g3;
	g[3]=g4;
	g[4]=g5;

    
	h[0]=h1;
	h[1]=h2;
	h[2]=h3;
	h[3]=h4;
	h[4]=h5;
	h[5]=h6;

srand(time(NULL));//随机过程初始化
success=FALSE;
counter=0;
gen=0;
wgen=0.75*(double)MAXGEN;//求惯性权重w变化的最大代数；
wdec=((double)MAXW-double(MINW))/wgen;//每代中w的递减直；


//初始化pos and pbestpos
for(j=0;j<DIM;j++)
{
	for(i=0;i<POPSIZE;i++)
	{
		pos[i][j]=fun_randval(lower[j],upper[j]);
		pbestpos[i][j]=pos[i][j];
	}
	    
}

//初始化速度
for(i=0;i<POPSIZE;i++)
{
	for(j=0;j<DIM;j++)
		vec[i][j]=fun_randval(double(MAXV),double(-MAXV));
}

evaluate(POPSIZE,pos,parvalue,weight);//评估种群，适应值存放在parvalue，罚值存放在weight中

//初始化历史最好位置的适应值
for(i=0;i<POPSIZE;i++)
{
	for(j=0;j<DIM;j++)
		pbestvalue[i]=parvalue[i];


}


//初始化全局最优微粒；
gbestpar=best(pos,pbestvalue,POPSIZE);


while((gen<MAXGEN)&&(success==FALSE))
{
	lastbest=pbestvalue[gbestpar];//暂时存放前一次进化之后的全局最好位置的相关信息；
	gen=gen+1;

	//更新惯性权重；
	if(gen<=wgen)
		w=MAXW-(gen-1)*wdec;

	//更新速度；
	for(i=0;i<POPSIZE;i++)
	{
	    for(j=0;j<DIM;j++)
		{
			vec[i][j]=w*vec[i][j]+(double)C1*fun_randval(0,1)*(pbestpos[i][j]-pos[i][j])+
				  (double)C2*fun_randval(0,1)*(pbestpos[gbestpar][j]-pos[i][j]);



		    if(vec[i][j]>MAXV)//速度的取值不能过大
			    vec[i][j]=MAXV;
		    else
			    if(vec[i][j]<-MAXV)
				    vec[i][j]=-MAXV;
		}

	}


	//更新位置；
	for(i=0;i<POPSIZE;i++)
	{
	    for(j=0;j<DIM;j++)
		{
			pos[i][j]=vec[i][j]+pos[i][j];
		    if(pos[i][j]>upper[j])
			    pos[i][j]=upper[j];
		    else
			    if(pos[i][j]<lower[j])
			        pos[i][j]=lower[j];
		}
	}

	evaluate(POPSIZE,pos,parvalue,weight);


        //更新每个个体的历史最好位置
	for(i=0;i<POPSIZE;i++)
		if(better(weight[i],fun_w(pbestpos[i]),parvalue[i],pbestvalue[i]))
		{
			pbestvalue[i]=parvalue[i];
			for(j=0;j<DIM;j++)
				pbestpos[i][j]=pos[i][j];
		}

	gbestpar=best(pos,pbestvalue,POPSIZE);//求取全局最好位置的下标；
	
	//计数连续多少次最好值没有发生改变
	if(fabs(lastbest-pbestvalue[gbestpar])<ERRGOAL)
			counter++;
	else
		counter=0;

        //如果连续5000次未找到更好位置就认为种群已经收敛
	if(counter>=5000)
		success=TRUE;
		
		

}

//打印相关信息
printf("run %d ,minimum:%-16.14f\n",curun,pbestvalue[gbestpar]);

for(j=0;j<DIM;j++)
printf("x%d = %-16.14f  ",j+1,pbestpos[gbestpar][j]);
printf("\n");

sum=sum+pbestvalue[gbestpar];//多次运行
sumgen=sumgen+gen;

//thebestvalue记录多次运行得到的最好位置
if((thebestvalue>pbestvalue[gbestpar])||(curun==0))
{
  thebestvalue=pbestvalue[gbestpar];
  for(i=0;i<DIM;i++)
	  thebestpos[i]=pbestpos[gbestpar][i];
}


}

//输出信息；
printf("The average minium : %-16.14f \n",sum/curun);
printf("It takes avgen:%d \n",sumgen/curun);

printf("The best position is:\n");
for(j=0;j<DIM;j++)
printf("x%d = %-16.14f  ",j+1,thebestpos[j]);
printf("\n");
printf("The minium is %-16.14f\n",thebestvalue);
printf("success!\n");

}



//函数体定义；
double fun_randval(double low,double high)
{
	double val;
	val=((double)(rand()%1001)/1000.0)*(high-low)+low;
	return(val);

}
void evaluate(int popsize,Particle pos[],double parvalue[],double w[])
{
	for(i=0;i<popsize;i++)
	{
	  parvalue[i]=fun_object(pos[i]);
	  weight[i]=fun_w(pos[i]);

	}

}


int best(double pos[][DIM],double pbestvalue[],int size)
{
	int g=0,i;
	for(i=1;i<size;i++)
		if(better(weight[i],weight[g],pbestvalue[i],pbestvalue[g]))
			g=i;
	return g;
}



double fun_w(double pos[])
{
	double result=0.0;
	for(int i=0;i<K;i++)
	{
		double temp=(*g[i])(pos);
		if(temp>0)
			result=result+pow(temp,2);
	}
	for(i=0;i<L;i++)
	{
		double temp=(*h[i])(pos);
		if(temp>0)
			result=result+pow(temp,2);
	}
	return result;
}

int better(double w1,double w2,double funcvalue1,double funcvalue2)
{
	if(w1<w2)
		return 1;
	else
		if(w1>w2)
			return 0;
		else
			if((funcvalue1)<=(funcvalue2))
				return 1;
			else
				return 0;
}

double fun_object(double *x)
{
	
	int i;
	double sum1=0.0;
	double sum2=0.0;
	double sum=0.0;
	for(i=1;i<6;i++)
		sum1+=i*cos((i+1)*x[0]+i);
	for(i=1;i<6;i++)
		sum2+=i*cos((i+1)*x[1]+i);
	sum=sum1*sum2;
	return sum;

}

double g1(double x[DIM])
{
	//return(x[0]*x[0]+x[1]*x[1]-9*9);
	return 0.0;
}
double g2(double x[DIM])
{
	return 0.0;
}
double g3(double x[DIM])
{
	return 0.0;
}
double g4(double x[DIM])
{
	return 0.0;
}
double g5(double x[DIM])
{
    return 0.0;
}

double h1(double x[DIM])
{
	return 0.0;
}
double h2(double x[DIM])
{
	return 0.0;
}
double h3(double x[DIM])
{
	return 0.0;
}
double h4(double x[DIM])
{
	return 0.0;
}
double h5(double x[DIM])
{
	return 0.0;
}
double h6(double x[DIM])
{
	return 0.0;
}