prog.c

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

/***************************************************************/
/* This is a simple genetic algorithm implementation where the */
/* evaluation function takes positive values only and the      */
/* fitness of an individual is the same as the value of the    */
/* objective function                                          */
/* author:Denis Cormier(North Carolina State University)       */
/*   Sita S.Raghavan(University of North Carolina at Charlotte)*/
/* Modified by:Bruce Chiang(China University of Geosciences)   */
/* Mayjor:Geophysics                                           */
/***************************************************************/
#include <stdio.h >
#include <stdlib.h>
#include <math.h  >
#include <time.h  >

/* Change any of these parameters to match your needs */

#define POPSIZE 50           /* population size */
#define MAXGEN 150          /* max. number of generations */
#define NVARS 2              /* no. of problem variables */
#define PXOVER 0.8           /* probability of crossover */
#define PMUTATION 0.1       /* probability of mutation */
#define TRUE 1
#define FALSE 0
#define PI 3.1415926
#define B   8                /*非均勻變異中的形狀因子*/
#define Lenth 100            /* length of generate Gauss distribution*/
int generation;              /* current generation no. */
int cur_best;                /* best individual */
FILE *galog;                 /* an output file */

struct genotype              /* genotype (GT), a member of the population */
{
  double gene[NVARS];        /* a string of variables */
  double fitness;            /* GT's fitness */
  double upper[NVARS];       /* GT's variables upper bound */
  double lower[NVARS];       /* GT's variables lower bound */
  double rfitness;           /* relative fitness */
  double cfitness;           /* cumulative fitness */
};

struct genotype population[POPSIZE+1];    /* population */
struct genotype newpopulation[POPSIZE+1]; /* new population; */
                                          /* replaces the */
                                          /* old generation */

/* Declaration of procedures used by this genetic algorithm */
void   initialize(void);
void   evaluate(void);
void   keep_the_best(void);
void   elitist(void);
void   select(void);
void   crossover(void);
void   Xover(int,int);
void   swap(double*,double*);
void   mutate(void);
void   report(void);
double randval(double, double);
double randvalG(double,double);
double randvalNu(int,int,double,double);
double obFun(double x[]);
/***************************************************************/
/* Initialization function: Initializes the values of genes    */
/* within the variables bounds. It also initializes (to zero)  */
/* all fitness values for each member of the population. It    */
/* reads upper and lower bounds of each variable from the      */
/* input file `gadata.txt'. It randomly generates values       */
/* between these bounds for each gene of each genotype in the  */
/* population. The format of the input file `gadata.txt' is    */
/* var1_lower_bound var1_upper bound                           */
/* var2_lower_bound var2_upper bound ...                       */
/***************************************************************/

void initialize(void)
{
	FILE *infile;
    int i, j;
    double lbound, ubound;
	if ((infile = fopen("gadata.txt","r"))==NULL)
	{
		fprintf(galog,"\nCannot open input file!\n");
        exit(1);
    }

/* initialize variables within the bounds */
	for (i=0; i<NVARS; i++)
	{
		fscanf(infile, "%lf",&lbound);
        fscanf(infile, "%lf",&ubound);
		for (j=0; j<POPSIZE; j++)
		{
			population[j].fitness = 0;
            population[j].rfitness= 0;
            population[j].cfitness= 0;
            population[j].lower[i]= lbound;
            population[j].upper[i]= ubound;
            population[j].gene[i] = randval(population[j].lower[i],population[j].upper[i]);
        }
	}
	fclose(infile);
}

/***********************************************************/
/* Random value generator: Generates a value within bounds */
/* Generate uniform random numbers,the mutation is uniform */
/* And I changed the form of the mutation,I used the Gauss */
/* and Non-uniform mutation ,also I modified some place    */
/* But I don't change it's structure,because I love it     */
/***********************************************************/

double randval(double low, double high)
{
	double val;
	val = low+((double)(rand()%1000)/1000.0)*(high-low);
	return(val);
}

/* Gauss mutation*/

double randvalG(double low,double high)
{
	int i=0;
	double val,sum=0.0;
	for(i=0;i<Lenth;i++)
	{
		sum+=(double)(rand()%1000)/1000.0;
	}
	return (val=(low+high)/2+(high-low)*(sum-Lenth/2)/(Lenth/2));
}

/***********************************************************/
/* Random value generator: Generates a value within bounds */
/* Generate Non-unform random numbers,also the mutation    */
/***********************************************************/

double randvalNu(int i,int j,double low,double high)
{
	double val;
	double r=generation/MAXGEN;
	double randN=(rand()%1000)/1000.0;
	if( (int) (randN+0.5)==0 )
		val=population[i].gene[j]+(high-population[i].gene[j])*(pow(randN*(1-r),B));
	if( (int) ( randN+0.5)==1 )
		val=population[i].gene[j]-(population[i].gene[j]-low)*(pow(randN*(1-r),B));
	return(val);
}

/*************************************************************/
/* Calculate the object function                             */
/* Each time this is changed, the code has to be recompiled. */
/* The current function is:  x[1]^2-x[1]*x[2]+x[3]           */
/* Because the principle of GA,it calculates the maxim valuse*/
/* however we will get the minimum value ,so we will change  */
/* the form of the objective fuction ,so I add the constant  */
/* I chose Shuber function as a test function,while X[-10 10]*/
/* This function has the global minimum value,F(X)=-186.731  */
/*************************************************************/

double obFun(double x[])//x[]表示反演的参数
{
	
	int i=0;
	double sum0=0.0,sum1=0.0;
	for(i=1;i<=5;i++)
	{
		sum0=sum0+i*cos((i+1)*x[0]+i);
		sum1=sum1+i*cos((i+1)*x[1]+i);
	}
	return(500-sum0*sum1);
}

/*************************************************************/
/* Evaluation function: This takes a user defined function.  */
/*************************************************************/

void evaluate(void)
{
	int mem;
	int i;
	double x[NVARS];
	for (mem = 0; mem < POPSIZE; mem++)
	{
		for (i=0;i<NVARS;i++)
            x[i] = population[mem].gene[i];
		population[mem].fitness =obFun(x);
      }
}

/***************************************************************/
/* Keep_the_best function: This function keeps track of the    */
/* best member of the population. Note that the last entry in  */
/* the array Population holds a copy of the best individual    */
/***************************************************************/

void keep_the_best()
{
	int mem;
    int i;
    cur_best = 0; /* stores the index of the best individual */
    for (mem = 0; mem < POPSIZE; mem++)
	{
		if (population[mem].fitness>population[POPSIZE].fitness)
		{
			cur_best = mem;
            population[POPSIZE].fitness=population[mem].fitness;
		}
	}
/* once the best member in the population is found, copy the genes */
    for (i = 0; i < NVARS; i++)
		population[POPSIZE].gene[i] = population[cur_best].gene[i];
}

/****************************************************************/
/* Elitist function: The best member of the previous generation */
/* is stored as the last in the array. If the best member of    */
/* the current generation is worse then the best member of the  */
/* previous generation, the latter one would replace the worst  */
/* member of the current population                             */
/****************************************************************/

void elitist()
{
	int i;
    double best, worst;             /* best and worst fitness values */
    int best_mem, worst_mem; /* indexes of the best and worst member */
    best = population[0].fitness;
    worst= population[0].fitness;
    for (i=0;i<POPSIZE-1;++i)
    {
		if(population[i].fitness>population[i+1].fitness)
		{  
            if (population[i].fitness>=best)
			{
				best = population[i].fitness;
                best_mem = i;
             }
            if (population[i+1].fitness<=worst)
            {
				worst=population[i+1].fitness;
                worst_mem=i+1;
            }
            }
		else