gas.cpp

非线性参数估计的遗传算法解。《非线性模型参数估计理论与方法》

//基本遗传算法源程序（二进制编码）
/*****************************************
    Simple Genetic Algorithm
	Programed by Tian yugang
******************************************/

//#include "stdafx.h"
#include "stdio.h"
#include "stdlib.h"
#include "time.h"
//#include "mfc.h"
#include "math.h"
#include "string.h"
#include "afx.h"
#include <vector>
//#include <graphics.h>
//#include "graph.h"

//The Definition of Constant

#define POPSIZE        500    // population size
#define MAXIMIZATION   1      // maximization flag
#define MINIMIZATION   2      // minimization flag

// The Definition of User Data
// (For different problem,there are some difference)

#define Cmax     100         // certain maximal value
#define Cmin     0           // certain minimal value

#define LENGTH1  14          // the chromosome length of 1st variable
#define LENGTH2  14          // the chromosome length of 2nd variable
                                       
#define CHROMLENGTH  LENGTH1 + LENGTH2 // the total length of chromosome
int FunctionMode = MINIMIZATION;       // optimization type(求最大值或最小值)

int PopSize      = 120;                 // population size

int MaxGeneration= 150;                // max number of generation

double Pc = 0.591;                      // probability of crossover

double Pm = 0.0089;                    // probability of mutation

// The Definition of Data Structure
 
struct  individual                     // data structure of individual                    
{
  char chrom[CHROMLENGTH + 1] ;        // a string of code representing individual

  double value;                        // object value of this individual
  double fitness;                      // fitness value of this individual
  
};

// The Definition of Global Variables   

int generation;                            // number of generation
int best_index;                            // index of best individual
int worst_index;                           // index of worst individual
struct individual  bestindividual;         // best individual of current generation
struct individual  worstindividual;        // worst individual of current generation
struct individual  currentbest;            // best individual by now
struct individual  population[POPSIZE];    // population

//Declaration of Prototype

void GenerateInitialPopulation(void);
void GenerateNextPopulation(void);
void EvaluatePopulation(void);
long DecodeChromosome(char *,int,int);
void CalculateObjectValue(void);
void CalculateFitnessValue(void);
void FindBestAndWorstIndividual(void);
void PerformEvolution(void);
void SelectionOperator(void);
void CrossoverOperator(void);
void MutationOperator(void);
void OutputTextReport(void);
int  randomize(int nMax);

//main program

void main(void)
{
   generation = 0;

   GenerateInitialPopulation();

   EvaluatePopulation();

   while(generation < MaxGeneration)
   {
     generation ++ ;
     
	 GenerateNextPopulation();

     EvaluatePopulation();

     PerformEvolution();

     OutputTextReport();
      
   }
}

//Function: Generate the first population.
//variable: none.

void GenerateInitialPopulation(void)
{
  int i,j;
  
  for(i=0; i <PopSize; i++)
  {
    for(j=0; j <CHROMLENGTH; j++)
	{
	   population[i].chrom[j] = (randomize(10) < 5)?  '0':'1';
	
	}
	population[i].chrom[CHROMLENGTH] = '\0' ;
  }
}


// Function: Initial the first generation.
// Variable:none.

void GenerateNextPopulation(void)
{
    SelectionOperator();
    CrossoverOperator();
    MutationOperator();
}

//Function: Evaluate population according to certain formula.
//variable:none.
void EvaluatePopulation(void)
{
   CalculateObjectValue();
   CalculateFitnessValue();
   FindBestAndWorstIndividual();
}

//Function: To decode a binary chromosome into a decimal integer
//variable: none.
//note: The returned value may be plus,or minus.For difference coding method,
//this value may be changed into "unsigned int"

long DecodeChromosome(char *string,int point,int length)
{
   int i;
   long decimal = 0L;
   char *pointer;

   for(i=0,pointer = string + point; i<length;i++,pointer++)
   {
      decimal  = decimal + ((*pointer - '0') << (length - 1 - i));
   }
   return(decimal);
}

//Function:To calculate object value.
//variable: none.
//note:For different problem,user must change these code.

void CalculateObjectValue(void)
{
  int i,j;
  long temp1,temp2;
  double x1,x2;
  static double l[5] = {4.2,3.25,2.52,1.95,1.51};           //观测值

  //rosen brock function
  for(i=0; i <PopSize; i++)
  {
     temp1 = DecodeChromosome(population[i].chrom,0,LENGTH1);
	 temp2 = DecodeChromosome(population[i].chrom,LENGTH1,LENGTH2);
	 x1 = 0.5 * temp1/(pow(2.0,14)-1) + 5.1;
	 x2 = 0.5 * temp2/(pow(2.0,14)-1) - 0.6;
     double sum = 0.0;
	 double t =0.0;
	 for(j=1;j<6;j++)
	 {
	   //sum = sum + (x1*x1*exp(j*2*x2) - 2*x1*l[j-1]*exp(j*x2));
		 t=x1*exp(j*x2);
		 sum = sum + fabs(t - l[j-1]);
	 }
	 population[i].value = sum;
  
  }
}

//Function: To calculate fitness value.
//variable: None.

void CalculateFitnessValue(void)
{
   int i ;
   double temp;

   for(i=0; i <PopSize; i++)
   {
      if(FunctionMode == MAXIMIZATION)    //     acquire the maximization value of object
	  {
	     if((population[i].value + Cmin) > 0.0) 
		 {
		    temp = Cmin + population[i].value;
		 }
		 else
		 {
		    temp = 0.0;
		 }
	  }
	  else if(FunctionMode == MINIMIZATION)  // acquire the minimization value of object  
	  {
	     if(population[i].value < Cmax)
		 {
		    temp = Cmax - population[i].value;
		 }
		 else
		 {
		    temp = 0.0;
		 }
	  }
	   population[i].fitness = temp;
     // population[i].fitness = 1.0/(sqrt(population[i].fitness)+0.0001); 
   }
}

//Function: To find out the best individual so far current generation.
//Variable: None.

void FindBestAndWorstIndividual(void)
{
   int i;
   double  sum = 0.0;
   //find out the best and worst individual of this generation.
   bestindividual = population[0];
   worstindividual = population[0];

   for(i=1; i <PopSize; i++)
   {
      if(population[i].fitness > bestindividual.fitness)
	  {
	     bestindividual = population[i];
		 best_index = i;
	  }
	  else if(population[i].fitness < worstindividual.fitness)
	  {
	     worstindividual = population[i];
		 worst_index = i;
	  }
	  sum = sum + population[i].fitness;
   }
   //find out the best individual so far
   if(generation == 0)    //initial the best individual
   {
      currentbest = bestindividual;
   }
   else if(bestindividual.fitness > currentbest.fitness)
   {
      currentbest = bestindividual; 
   }
}

//Function: To perform evolution operation based on elitise model.
//elitist model is to replace the worst individual of this generation by current best one.

void PerformEvolution(void)
{
  if(bestindividual.fitness > currentbest.fitness)
  {
     currentbest = population[best_index];
  }
  else
  {
     population[worst_index] = currentbest;
  }
}

//Function: To reproduce a chromosome by proportional selection.  Variable : none

void SelectionOperator(void)
{
   int i,index;
   double p,sum = 0.0;
   double cfitness[POPSIZE];      //cumulative fitness value
   struct individual newpopulation[POPSIZE];

   //caculate relative fitness
   for(i=0; i < PopSize; i++)
   {
      sum = sum + population[i].fitness;
   }
   for(i=0; i < PopSize; i++)
   {
      cfitness[i] = population[i].fitness/sum;
   }

   //calculate cumulative fitness
   for(i=1; i < PopSize; i++)
   {
      cfitness[i] = cfitness[i-1] + cfitness[i];
   }

   //selection operator
   for(i=0; i < PopSize; i++)
   {
      p = rand()%1000 / 1000.0;
	  index = 0;
	  while(p > cfitness[index])
	  {  
		 index ++;    
	  }
	  newpopulation[i] = population[index];
   }
   for(i=0; i <PopSize; i++)
   {
      population[i] = newpopulation[i];
   }
}


//function: Crossover two chromosome by means of one-point crossover.
//variable:none

void CrossoverOperator(void)
{
  int i, j;
  int index[POPSIZE];
  int point,temp;
  double p;
  char ch;

  //make a pair of individual randomly
  for(i=0; i < PopSize; i++)
  {
     index[i] = i;
  }
  for(i=0; i < PopSize; i++)
  {
     point = randomize(PopSize - i);

	 temp = index[i];
	 index[i] = index[point + i];
     index[point + i] = temp;
  }

  //one-point crossover operation
  for(i=0; i < PopSize-1; i = i+2 )
  {
     p = rand()%1000/1000.0;
	 if(p<Pc)
	 {
	    point = randomize(CHROMLENGTH - 1) + 1;
		for(j=point; j<CHROMLENGTH; j++)
		{
		  ch = population[index[i]].chrom[j];
          population[index[i]].chrom[j] = population[index[i+1]].chrom[j];
          population[index[i+1]].chrom[j] = ch;
		}
	 }
  }

}

//function: Mutation of a chromsome.

void MutationOperator(void)
{
   int i,j;
   double p;

   //bit mutation 
   for(i=0; i<PopSize; i++)
   {
      for(j=0; j<CHROMLENGTH; j++)
	  {
	     p = rand()%1000/1000.0;
		 if(p < Pm)
		 {
			 population[i].chrom[j] = (population[i].chrom[j] == '0')?  '1' : '0' ;
		 }
	  }
   }
} 



//Function: crossover two chromosome by meansof one_point crossover.  Variable:None

void OutputTextReport(void)
{
   int i;
   double sum;      //temporary sum
   double average;  //average of population object value.

////////calculate average object value
   sum = 0.0;
   for(i=0; i <PopSize;i++)
   {
     sum  = sum + population[i].value;
   }

   average = sum/PopSize;

   //printf results of this population.
   

    long temp1,temp2;
    double x1,x2;
    temp1 = DecodeChromosome(currentbest.chrom,0,LENGTH1);
	temp2 = DecodeChromosome(currentbest.chrom,LENGTH1,LENGTH2);

	 x1 = 0.5 * temp1/(pow(2.0,14)-1) + 5.1;
	 x2 = 0.5 * temp2/(pow(2.0,14)-1) - 0.6;


	/*CString WriteStr , OutFile ;
	CStdioFile FileOut;
	OutFile = "e:\\out.txt" ;
	if(!FileOut.Open(OutFile , CFile::modeCreate | CFile::modeWrite))
	{
//		AfxMessageBox("打开输出文件错误!" , MB_ICONSTOP , NULL);
		return ;
	}
	
		   
	
	WriteStr.Format(" \ngen= %d, best= %f,x1= %f,x2= %f",generation,currentbest.value,x1,x2);
	FileOut.WriteString(WriteStr);
	
   /*printf("chromosome = ");
   for(i =0; i< CHROMLENGTH; i++)
   {
	   printf("%c", currentbest.chrom[i]);
   }*/
   printf(" \ngen= %d, best= %f,x1= %f,x2= %f",generation,currentbest.value,x1,x2);	
   printf("\n");
  // FileOut.Close();
}

int randomize(int nMax)
{
	int nCreate = rand();
	int nValue = nMax * nCreate / RAND_MAX;
	return nValue;
}