nichegamodify.c

采用FORTRAN编制的小生境遗传算法反演程序

    } 
}
/**************************************************************/
/* Selection function: Standard proportional selection for    */
/* maximization problems incorporating elitist model - makes  */
/* sure that the best member survives                         */
/**********??????????????????????????????????******************/

void select(void)
{
	int mem, i, j;
    double sum = 0;
    double p;
	int expectN,COUNT=0;

/* find total fitness of the population */
    for (mem=0; mem<POPSIZE;mem++)
	{
      sum+=population[mem].fitness;
    }

/* calculate cumulation fitness */
    for (mem=0; mem<POPSIZE;mem++)
	{
		population[mem].cfitness =  population[mem].fitness/sum;
    }
	for(mem=0;mem<POPSIZE;mem++)
	{
		expectN=(int)(POPSIZE*population[mem].cfitness);
		if(expectN==1)     //也就是当population[mem].fitness大于全部适应值的平均值时
		{
						newpopulation[COUNT]=population[mem]; //大于适应值的平均值的个体全部保留
		    COUNT++;
		}//???????????????疑问?????????????????????????????
	}
    
/* calculate cumulative fitness */
    for (mem=1;mem<POPSIZE;mem++)
	{
		population[mem].cfitness=population[mem-1].cfitness+population[mem].cfitness;
    }
	//調試程序階段
	/*
	for(i=0;i<POPSIZE;i++)
			printf("第%4d個種群 相對適應度輸出=%10.7f\n",(i+1),population[i].rfitness);
        for(i=0;i<POPSIZE;i++)
            printf("第%4d個種群 纍計適應度輸出=%10.7f\n",(i+1),population[i].cfitness);
*/
	//将前面N个个体直接参与下代
	
//	for (i=0; i<N;i++)
//		newpopulation[i] =population[i] ;
	for (i=COUNT;i<POPSIZE;i++)
	{
		p = rand()%1000/1000.0;
      	for (j=0;j<POPSIZE;j++)
				if (p>=population[j].cfitness&&p<population[j+1].cfitness)
				{
					newpopulation[i] = population[j+1];
				//	printf("第=%4d个新群体为原来第j=%4d个群体\n",(i+1),(j+1));
				}		
	}
/* finally select survivors using cumulative fitness. */
/*
    for (i = 0; i < POPSIZE; i++)
	{
		p = rand()%1000/1000.0;
        if (p<population[0].cfitness)
		{
			newpopulation[i] = population[0];		
			printf("第1个群体被选中\n");
		}
		else
		{
			for (j=0;j<POPSIZE;j++)
				if (p>=population[j].cfitness&&p<population[j+1].cfitness)
				{newpopulation[i] = population[j+1];
				printf("第=%4d个新群体为原来第j=%4d个群体\n",(i+1),(j+1));}
		}
		
	}
	*/
/* once a new population is created, copy it back */

    for (i=0; i<POPSIZE;i++)
		population[i] = newpopulation[i];      
}

/***************************************************************/
/* Crossover selection: selects two parents that take part in  */
/* the crossover. Implements a single point crossover          */
/***************************************************************/

void crossover(void)
{
	int mem, one=0;
    int first=0; /* count of the number of members chosen */
    double x;

    for (mem=0;mem<POPSIZE;++mem)
	{
		x=rand()%1000/1000.0;
        if (x<PXOVER)
		{
            ++first;
            if (first%2==0)
			{
				Xover(one, mem);
			//	printf("第%4d个和第%4d群体选择了交叉\n",one,mem);
			}
            else
                  one = mem;
         }
      }
}
/**************************************************************/
/* Crossover: performs crossover of the two selected parents. */
/**************************************************************/

void Xover(int one, int two)
{
int i=0;
int point; /* crossover point */
double x1=0,x2=0; 
double alph=(rand()%1000/1000.0);

/* select crossover point */
if(NVARS>1)
   {
	if(NVARS==2)
         point=1;
	else
         point=(rand()%(NVARS-1))+1;
	for (i=0;i<point;i++)//swap(&population[one].gene[i], &population[two].gene[i]);
	{
		x1=population[one].gene[i];
	    x2=population[two].gene[i];
	    population[one].gene[i]=alph*x2+(1.0-alph)*x1;
	    population[two].gene[i]=alph*x1+(1.0-alph)*x2;  //算术交叉
	}

   }
}

/*************************************************************/
/* Swap: A swap procedure that helps in swapping 2 variables */
/*************************************************************/

void swap(double *x, double *y)
{
	double temp;
	temp = *x;
	*x = *y;
	*y = temp;
}

/**************************************************************/
/* Mutation: Random uniform mutation. A variable selected for */
/* mutation is replaced by a random value between lower and   */
/* upper bounds of this variable                              */
/**************************************************************/

void mutate(void)
{
	int i, j;
    double lbound, hbound;
    double x;
    for (i=0;i<POPSIZE;i++)
		for (j=0;j<NVARS;j++)
		{
			x = rand()%500/500.0;
            if (x<PMUTATION)
			{
				/* find the bounds on the variable to be mutated */
				lbound = population[i].lower[j];
				hbound = population[i].upper[j];
				population[i].gene[j] = randvalNu(i,j,lbound, hbound);
            }
		}
}
/***************************************************************/
/* nicheGA function: calculate the niche value                 */
/***************************************************************/
/*我觉得这里实现了小生境的算法思想，今天中午突然有所悟而想到的 */
void nicheGA(void)
{
	int i=0,j=0,k=0;
	double sum=0.0;
	struct genotype tempGA;
	for(i=0;i<POPSIZE;i++)
		nicheP[i]=population[i];
    for(i=POPSIZE;i<(POPSIZE+N);i++)
		nicheP[i]=popMemory[i-POPSIZE];

	for(i=0;i<(POPSIZE+N-1);i++)
	{
		for(j=i+1;j<(POPSIZE+N);j++)
		{
			if(hamming1(i,j)<HD)
			{
	
				if(nicheP[j].fitness>nicheP[i].fitness)
					nicheP[i].fitness=Penalty*nicheP[i].fitness;
				else
					nicheP[j].fitness=Penalty*nicheP[j].fitness;
				
			}
		}
                                                
	}
 
	//對(POPSIZE+N)個群體進行排序	
	for (i=0;i<(POPSIZE+N-1);i++)
	{
		for(j=i+1;j<(POPSIZE+N);j++)
		if (nicheP[i].fitness<nicheP[j].fitness)
		{
			tempGA=nicheP[i];
			nicheP[i]=nicheP[j];
			nicheP[j]=tempGA;
		}
	}
	//記憶前N個個体
    for(i=0;i<N;i++)
		popMemory[i]=nicheP[i];
	//把前M個個体給作爲下一代的群體
	for(i=0;i<POPSIZE;i++)
		population[i]=nicheP[i];
    for(i=0;i<N;i++)
		population[POPSIZE-1-i]=popMemory[i];
	if (population[0].fitness>=population[POPSIZE].fitness)
    {
		population[POPSIZE]=population[0];
    }
	population[POPSIZE-1]=population[POPSIZE];

	//else
    //{
    //for (i=0; i<NVARS; i++)
     //  population[worst_mem].gene[i] = population[POPSIZE].gene[i];
    //population[worst_mem].fitness = population[POPSIZE].fitness;
  //  } 




}

/***************************************************************/
/* simulate function: simulates progress of the simulation. Data   */
/* dumped into the  output file are separated by commas        */
/***************************************************************/

void simulate(void)
{
	int i;
    double best_val;            /* best population fitness */
    double avg;                 /* avg population fitness */
    double stddev;              /* std. deviation of population fitness */
    double sum_square;          /* sum of square for std. calc */
    double square_sum;          /* square of sum for std. calc */
    double sum;                 /* total population fitness */
	sum = 0.0;
    sum_square = 0.0;
	for (i=0; i<POPSIZE; i++)
	{
		sum += population[i].fitness;
        sum_square += population[i].fitness * population[i].fitness;
    }            
	avg = sum/(double)POPSIZE;
	square_sum = avg * avg * POPSIZE;           // sum_square》=square_sum
	stddev     =sqrt((sum_square - square_sum)/(POPSIZE - 1));
	best_val=-population[POPSIZE].fitness+5000.0;
//	if(fabs(best_val+186.731)<0.0009)
//	{
//		printf("Gneration=%3d\n",generation);
//		exit(0);
//	}
	fprintf(galog, "\n%6.2f\n", best_val);
}

/**************************************************************/
/* Main function: Each generation involves selecting the best */
/* members, performing crossover & mutation and then          */
/* evaluating the resulting population, until the terminating */
/* condition is satisfied                                     */
/**************************************************************/

void main(void)
{

	double bi[201],ri[RECEIVES],OBS[RECEIVES];//定义数组
	double dt=0.5,f0=60;  //数组赋值
	double t[3];//时间
	double v[4]={1500,1000,1500,2000};
	double SUM=0.0;
    int i=0,j=0,nw=100;
   
    Penalty=pow(10,-3);
	population[POPSIZE].fitness=0.0;//初始化存放最大適應度值的群體
	//初始化反射系数数组		  
	for(i=0;i<RECEIVES;i++)
		ri[i]=0.0;

	t[0]=2*h[0]/v[0];
	j=(int)(1000*t[0]/dt);
	ri[j]=(v[1]-v[0])/(v[1]+v[0]); //1层 

	
	for(i=1;i<NVARS;i++)
	{
		t[i]=t[i-1]+2*h[i]/v[i];
		j=(int)(1000*t[i]/dt);
		ri[j]=(v[i+1]-v[i])/(v[i+1]+v[i]);//2-4层
	}	

	
	for(i=-nw;i<nw+1;i++)
	{
		double a=(0.001*3.1415926*f0*i*dt);
        bi[i+nw]=(1.0-2.0*a*a)*exp(-a*a);//求取子波，0相位，
        
	}

	
	//Convolution   子波和参数做卷积
	for(i=0;i<RECEIVES;i++)//从第一道循环
	{
		double sum=0.0;
	    for(j=0;j<2*nw+1;j++)
		{
			if(i-j>=0&&i-j<=2*nw)
		    sum=sum+bi[j]*ri[i-j+1];
		}
		OBS[i]=sum;
		
	}
    for (i=0;i<NVARS;i++)
		population[POPSIZE].gene[i]=0.0;
	srand((unsigned int)time(NULL));    //激活rand()%函数
    if ((galog = fopen("galog.txt","w"))==NULL)
		exit(1);    
	generation = 0;
	fprintf(galog, "\n generation  best  average  standard \n");
    fprintf(galog, " number      value fitness  deviation \n");
	initialize();    //初始化得到群体
	evaluate();          //计算得到群体中各个个体的适应值
	searchBest();    //按适应值由大到小排序，并记录出前N个个体，
	                 //且把最大的一个记录到population[POPSIZE]
    while(generation<MAXGENS)
	{		
		generation++;
	//	Gfactor=sqrt(1-((double)generation/MAXGENS));
        select();
		crossover();		
        mutate();
        evaluate();	      //计算得到各个新个体的适应值
		simulate();	        //计算并在输出文本匡中打印出来这次叠代的结果
		nicheGA();
        elitist();          //输出本次叠代产生的最大最小函数值
		
		//break;
    }
	fprintf(galog,"\n\n Simulation completed\n");
    fprintf(galog,"\n Best member: \n");
//	for(j=0;j<POPSIZE;j++)
//	{
		fprintf (galog,"\n");
		for (i = 0; i < NVARS; i++)
	//	{
			fprintf (galog," X(%d) = %8.2f",i,population[POPSIZE].gene[i]);
	//	}
	//	fprintf(galog,"\n\n Best fitness = %6.3f",(-population[j].fitness)+0);
//	}
	fprintf(galog,"\n\n Best fitness = %6.3f",(-population[POPSIZE].fitness)+5000.0);
    fclose(galog);
    printf("Success\n");
}
/***************************************************************/