多目标优化的遗传算法及其改进(浮点数编码)，对多个函数进行测试

本文链接：https://blog.csdn.net/yanguilaiwuwei/article/details/46699801

上一篇博客主要写了遗传算法的基本操作，主要是对单目标优化的算法，经过测试函数，可以知道算法的准确度十分高，但是仍然会存在陷入局部最优的情况。想了解上一篇博客的网友可以点击：http://blog.csdn.net/yanguilaiwuwei/article/details/46670805

为了解决算法陷入局部最优的现象，本文主要采用以下改进算法：把每一代种群中最优的一定数量的个体，无条件的遗传到下一代中，所以种群的最优适应度一定会随着遗传代数的增加不断升高或者不变(达到最大值时不再变化)，通过这种方法可以大大减小遗传算法陷入局部最优的概率。

本文在上一篇的基础上，把单目标优化，扩展到了多目标优化，C语言代码如下：

#include<stdio.h>
#include<stdlib.h> 
#include<math.h>
#include<time.h> 

#define M 100				//种群数量
#define N 2					//变量个数
#define PI 3.1415926
#define PC 0.8				//交叉概率
#define PM 0.4				//变异概率
#define PA 0.01				//交叉因子

#define COPYNUM 4


struct Node
{
	double Pmember[N];
	double Myfitness;		//Myfitness是适应度
	double Myfitsum;		//Myfitsum是适应度占总体适应度的百分比，然后从第一个个体往后累加，主要用于选择操作
}Nownode[M],Nextnode[M],TempNode[COPYNUM];	//本代群体和下一代群体,TempNode存放每一代最优的COPYNUM个个体

int nodeindex[M];			//交叉时随机配对，存放配对的群体下标
int T=0;
double XMIN[N]={-10,-10};//每一个参数的最小值
double XMAX[N]={10,10};//每一个参数的最大值


double fx(double *x);
int sortnode();
int copym(int n);
int copyn(int n);
int calfitsum();
int calfitness();
double randn(int index);
int initpopulation();
int assignment(struct Node *node1,struct Node *node2);
int copypopulation();
int isrepeat(int temp,int n);
int crossover();
int mutation();
int findmaxfit();
int diplaynode();


double fx(double *x)			//根据x计算fx
{
	double y;
	int i=0,px[N]={0};
	y=0.5-(pow(sin(sqrt(pow(x[0],2)+pow(x[1],2))),2)-0.5)/pow((1+0.001*(pow(x[0],2)+pow(x[1],2))),2);//schaffer函数,通过测试
	//y=3+pow(x[0],2)-cos(18*x[0])+pow(x[1],2)-cos(18*x[1]);//通过测试
	//y=4*pow(x[0],2)-2.1*pow(x[0],4)+pow(x[0],6)/3+x[0]*x[1]-4*pow(x[1],2)+pow(x[1],4)+5;//通过测试
	//y=pow(x[1]-5.1*pow(x[0],2)/(4*PI*PI)+5*x[0]/PI-6,2)+10*(1-1/(8*PI))*cos(x[0])+10;//通过测试
	//y=x[0]*exp(-pow(x[0],2)-pow(x[1],2))+1;//通过测试
	//y=sin(sqrt(pow(x[0],2)+pow(x[1],2)))/sqrt(pow(x[0],2)+pow(x[1],2))+1;//通过测试
	//y=1.0/4000*(pow(x[0],2)+pow(x[1],2))-cos(x[0])*cos(x[1]/sqrt(2.0))+1;//Griewank函数通过测试
	//y=pow(x[0],2)-10*cos(2*PI*x[0])+10+pow(x[1],2)-10*cos(2*PI*x[1])+10;//Rastrigin函数通过测试
	
	/*for(i=0;i<N;i++)
	{
		if(x[i]>=0)
			px[i]=1;
		else
			px[i]=0;
	}
	y=px[1]*(1+px[0])+fabs(x[0]+50*px[1]*(1-2*px[0]))+fabs(x[1]+50*(1-2*px[1]));*///Tripod函数通过测试
	//y=20+exp(1)-20*exp(-0.2*sqrt((0.5)*(pow(x[0],2)+pow(x[1],2))))-exp(0.5*(cos(2*PI*x[0])+cos(2*PI*x[1])));//Ackley函数通过测试
	//y=100*pow((pow(x[0],2)-x[1]),2)+pow((1-x[0]),2)+1;//Rosenbrock函数范围缩小时可找到,或者把y=1/y改成y=exp(1/y)扩大个体适应度之间的差别
	//y=pow(x[0],2)+pow(x[1],2)+1;//NDparabola函数通过测试
	y=1/y;
	return y;
}


int sortnode()
{
	int i,j;
	struct Node temp;
	for(i=0;i<M;i++)
	{
		for(j=i;j<M;j++)
		{
			if(Nownode[i].Myfitness<Nownode[j].Myfitness)
			{
				assignment(&temp,&Nownode[i]);
				assignment(&Nownode[i],&Nownode[j]);
				assignment(&Nownode[j],&temp);
			}
		}
	}
	calfitsum();
	return 0;
}

int copyn(int n)
{
	int i;
	sortnode();
	for(i=0;i<n;i++)
	{
		assignment(&Nownode[M-1-i],&TempNode[i]);
	}
	calfitsum();
	return 0;
}


int copym(int n)
{
	int i;
	sortnode();
	for(i=0;i<n;i++)
	{
		assignment(&TempNode[i],&Nownode[i]);
	}
	return 0;
}

int calfitsum()
{
	int i;
	Nownode[0].Myfitsum=Nownode[0].Myfitness;
	for(i=1;i<M;i++)
	{
		Nownode[i].Myfitsum=Nownode[i].Myfitness+Nownode[i-1].Myfitsum;//每一个Myfitsum都是自己的适应度加上前一个的Myfitsum
	}


	for(i=0;i<M;i++)
	{
		Nownode[i].Myfitsum=Nownode[i].Myfitsum/Nownode[M-1].Myfitsum;//每一个Myfitsum除以所有适应度之和，使Myfitsum为0~1之间
	}
	return 0;
}

int calfitness()			//计算适应度值
{
	int i;
	double minfitness,maxfitness,avefitness=0;
	double temp;
	minfitness=Nownode[0].Myfitness=fx(Nownode[0].Pmember);
	maxfitness=minfitness;
	avefitness=maxfitness;
	for(i=1;i<M;i++)
	{	
		Nownode[i].Myfitness=fx(Nownode[i].Pmember);
		avefitness+=Nownode[i].Myfitness;
		if(minfitness>Nownode[i].Myfitness)
		{
			minfitness=Nownode[i].Myfitness;
		}
		if(maxfitness<Nownode[i].Myfitness)
		{
			maxfitness=Nownode[i].Myfitness;
		}
	}
	
	if(minfitness<0)//如果有负的适应度值，就把所以的适应度都加上一个数，使适应度全都为正数
	{
		temp=minfitness;
		Nownode[0].Myfitness+=-temp;
		avefitness=Nownode[0].Myfitness;
		maxfitness=Nownode[0].Myfitness;
		minfitness=Nownode[0].Myfitness;
		for(i=1;i<M;i++)
		{
			Nownode[i].Myfitness+=-temp;
			avefitness+=Nownode[i].Myfitness;
			if(minfitness>Nownode[i].Myfitness)
			{
				minfitness=Nownode[i].Myfitness;
			}
			if(maxfitness<Nownode[i].Myfitness)
			{
				maxfitness=Nownode[i].Myfitness;
			}	
		}
	}

	calfitsum();
	return 0;
}

double randn(int index)		//产生XMIN到XMAX之间的随机数
{
	return XMIN[index]+1.0*rand()/RAND_MAX*(XMAX[index]-XMIN[index]);
}

int initpopulation()	//初始化种群
{
	int i,j;
	for(i=0;i<M;i++)
	{
		for(j=0;j<N;j++)
		{
			Nownode[i].Pmember[j]=randn(j);
		}
	}
	calfitness();		//计算适应度
	return 0;
}

int assignment(struct Node *node1,struct Node *node2)//把node2的值赋值给node1
{
	int j;
	for(j=0;j<N;j++)
	{
		node1->Pmember[j]=node2->Pmember[j];
	}
	node1->Myfitness=node2->Myfitness;
	node1->Myfitsum=node2->Myfitsum;
	return 0;
}

int copypopulation()		//复制操作
{
	int i,num=0;
	double temp;
	while(num<M)
	{
		temp=1.0*rand()/RAND_MAX;
		for(i=1;i<M;i++)
		{
			if(temp<=Nownode[0].Myfitsum)
			{
				assignment(&Nextnode[num++],&Nownode[0]);//把第一个个体复制到下一代
				break;
			}
			if(temp>=Nownode[i-1].Myfitsum&&temp<=Nownode[i].Myfitsum)//把第i个个体复制到下一代
			{
				assignment(&Nextnode[num++],&Nownode[i]);
				break;
			}
		}
	}
	for(i=0;i<M;i++)
	{
		assignment(&Nownode[i],&Nextnode[i]);	//更新本代个体
	}
	calfitsum();
	//calfitness();		//计算适应度
	return 0;
}

int isrepeat(int temp,int n)	//产生随机下标判断是否重复
{
	int i;
	for(i=0;i<n;i++)
	{
		if(nodeindex[i]==temp)
			return 1;
	}
	return 0;
}

int crossover()
{
	int i,j,temp;
	double temp_pc;
	for(j=0;j<N;j++)
	{
		for(i=0;i<M;i++)		//产生交叉点的下标
		{
			do {
				temp=rand()%M;
			} while(isrepeat(temp,i));
			nodeindex[i]=temp;
		}
		for(i=0;i<M;i=i+2)
		{
			temp_pc=1.0*rand()/RAND_MAX;	//如果满足交叉的条件，就开始交叉
			if(temp_pc<=PC)
			{
				Nownode[nodeindex[i]].Pmember[j]=PA*Nownode[nodeindex[i+1]].Pmember[j]+(1-PA)*Nownode[nodeindex[i]].Pmember[j];
				Nownode[nodeindex[i+1]].Pmember[j]=PA*Nownode[nodeindex[i]].Pmember[j]+(1-PA)*Nownode[nodeindex[i+1]].Pmember[j];
				
			}	
		}
	}
	calfitness();		//计算适应度
	return 0;
}


int mutation()			//变异操作
{
	int i,j,temp;
	double k=0.7,temp_pm;
	for(j=0;j<N;j++)
	{
		for(i=0;i<M;i++)
		{
			temp_pm=1.0*rand()/RAND_MAX;
			if(temp_pm<=PM)	//如果满足变异条件，就开始变异
			{
				temp=rand()%2;
				if(temp==0)
				{
					Nownode[i].Pmember[j]=Nownode[i].Pmember[j]+k*(XMAX[j]-Nownode[i].Pmember[j])*1.0*rand()/RAND_MAX;
				}
				else
				{
					Nownode[i].Pmember[j]=Nownode[i].Pmember[j]-k*(Nownode[i].Pmember[j]-XMIN[j])*1.0*rand()/RAND_MAX;
				}
			}
		}
	}
	calfitness();	//计算适应度
	return 0;
}

int findmaxfit()//找到适应度最大的个体
{
	int i,index=0;
	double temp=0;
	for(i=0;i<M;i++)
	{
		if(temp<Nownode[i].Myfitness)
		{
			index=i;
			temp=Nownode[i].Myfitness;
		}
	}
	return index;
}

int displaynode()
{
	int i,j;
	printf("下标\t");
	for(j=0;j<N;j++)
	{
		printf("x%d\t\t",j+1);
	}
	printf("适应度\t\t适应度占比\n");
	for(i=0;i<M;i++)
	{
		printf("%d\t",i);
		for(j=0;j<N;j++)
		{
			printf("%.4lf\t\t",Nownode[i].Pmember[j]);
		}
		printf("%.4lf\t\t%.4lf\n",Nownode[i].Myfitness,Nownode[i].Myfitsum);
	}
	return 0;
}

int main()
{
	int j,index;
	double x[2]={0.3432,0.4399};
	int num=0,num1=0,num2=0;
	T=0;
	srand(time(NULL));
	initpopulation();
	while(T++<500)
	{
		copym(COPYNUM);
		copypopulation();
		crossover();
		mutation();
		copyn(COPYNUM);
	}
	index=findmaxfit();
	printf("下标\t");
	for(j=0;j<N;j++)
	{
		printf("x%d\t\t",j+1);
	}
	printf("适应度\t\t函数值\n");
	printf("%d\t",index);
	for(j=0;j<N;j++)
	{
		printf("%.4lf\t\t",Nownode[index].Pmember[j]);
		x[j]=Nownode[index].Pmember[j];
	}
	printf("%.4lf\t\t%lf\n",(Nownode[index].Myfitness),fx(Nownode[index].Pmember));	
	return 0;
}

本文所采用的测试函数主要有以下几个：

1.Schaffer函数

10.

11.

12.

上述函数在matlab画出图像的代码如下所示：

clear;
clc;

xx=-10:0.1:10;
yy=-10:0.1:10;
[x,y]=meshgrid(xx,yy);
z1=0.5-(sin(sqrt(x.^2+y.^2)).^2-0.5)./((1+0.001*(x.^2+y.^2)).^2);%Schaffer函数
mesh(x,y,z1);
set(gca,'fontsize',20);
title('Schaffer函数','fontsize',20);
z2=x.^2-cos(18*x)+y.^2-cos(18*y)+3;%后面加了3
mesh(x,y,z2);
title('z2=x^2-cos(18*x)+y^2-cos(18*y)+3','fontsize',20);
z3=4*x.^2-2.1*x.^4+x.^6/3+x.*y-4*y.^2+y.^4+5;%后面加了5
mesh(x,y,z3);
title('z3=4*x^2-2.1*x^4+x^6/3+x*y-4*y^2+y^4+5','fontsize',20);
z4=(y-5.1*x.^2/4/pi/pi+5*x/pi-6).^2+10*(1-1/8/pi)*cos(x)+10;
mesh(x,y,z4);
title('z4=(y-5.1*x^2/4/pi/pi+5*x/pi-6)^2+10*(1-1/8/pi)*cos(x)+10','fontsize',20);
z5=x.*exp(-x.^2-y.^2)+1;
mesh(x,y,z5);
title('z5=x*exp(-x^2-y^2)+1','fontsize',20);
z6=sin(sqrt(x.^2+y.^2)+1e-6)./(sqrt(x.^2+y.^2)+1e-6)+1;
mesh(x,y,z6);
title('z6=sin(sqrt(x^2+y^2))/(sqrt(x^2+y^2))+1','fontsize',20);
z7=1/4000*(x.^2+y.^2)-cos(x).*cos(y/sqrt(2))+1;%Griewank函数
mesh(x,y,z7);
title('Griewank函数','fontsize',20);
z8=(x.^2-10*cos(2*pi*x)+10)+(y.^2-10*cos(2*pi*y)+10);%Rastrigin函数
mesh(x,y,z8);
title('Rastringin函数','fontsize',20);
px1=((x)>=0);
px2=((y)>=0);
z9=px2.*(1+px1)+abs(x+50*px2.*(1-2*px1))+abs(y+50*(1-2*px2));%Tripod函数
mesh(x,y,z9);
title('Tripod函数','fontsize',20);
z10=20+exp(1)-20*exp(-0.2*sqrt((0.5)*(x.^2+y.^2)))-exp(0.5*(cos(2*pi*x)+cos(2*pi*y)));%Ackley函数
mesh(x,y,z10);
title('Ackley函数','fontsize',20);
z11=100*(x.^2-y).^2+(1-x).^2;%Rosenbrock函数
mesh(x,y,z11);
title('Rosenbrock函数','fontsize',20);
z12=x.^2+y.^2;
mesh(x,y,z12);
title('NDparabola函数','fontsize',20);

在matlab中可以不仅可以画出上述函数的图像，还可以得到函数的最值点，然后用来检测遗传算法得到系统最值点的准确性，以Schaffer函数为例，matlab得到最值点的代码如下所示：

f=@(x)-1*(0.5-(sin(sqrt(x(1).^2+x(2).^2)).^2-0.5)./((1+0.001*(x(1).^2+x(2).^2)).^2));
[minx fmin]=fminsearch(f,[1 1]);

fminsearch是寻找多变量函数的最小值，这里要寻找最大值，所以把Schaffer函数取相反数。得到的minx是函数的最小值点，fmin是最小值点的函数值，后面的[1,1]是对应x的初始值。得到的结果是minx为[0,0]，fmin为-1,即Schaffer的最大值在[0,0]取到，并且函数的最大值为1。

遗传算法要求最大值，所以适应度不需要变换，直接可用函数值作为适应度的值，当x0和x1的取值范围都是[-5,5]时，得到的结果如下：

可以看到算法得到了最优值。

经验证，当取值范围是[-5,5]时，算法准确度为100%，但是当取值范围[-10,10]时，准确度为96%以上。

当函数要求最小值时，如果函数值均大于零，就把适应度取倒数，如果函数中含有小于零的值，就把函数值都加上一个正数，使函数值全都大于零，再取倒数。

经过验证，此算法对所有的测试函数都有效。其中，当取值范围为[-10,10]时，对于Rosenbrock函数效果不是很好，因为此函数在最小值附近，函数值下降的很慢，所以收敛效果不佳，但是当取值范围为[-2,2]时，准确度能达到100%，另外，还可以改变适应度的值，按照上面的方法，可以在Rosenbrock函数值加上1，然后取倒数作为适应度值，但是由于函数值在最小值附近下降的太慢，所以适应度值区分度不大，这里可以用y=exp(1/y)作为适应度值，取指数之后，适应度之间的差别会增大，所以加快收敛速度。此外，经过测试，当取值范围为[-10,10]时，假如把种群数量设置为1000，也是能收敛到最小值点的，但是算法运行时间明显增长。这些说明：适应度的取法对算法的收敛速度有影响，此外，当种群数量足够大，也能够收敛到最小值点，但是运行时间会增大，所以可以采用折中的办法：先用大的取值范围得到最值点的大概范围，然后再把算法的范围缩小，得到精确地最值点。