Genetic Algorithm 遗传算法

遗传算法是计算数学中用于解决最优化的搜索算法，是进化算法的一种。进化算法最初是借鉴了进化生物学中的一些现象而发展起来的，这些现象包括遗传、突变、自然选择以及杂交等。


学习遗传算法的时候,看到资料是这一页:
http://math.ecnu.edu.cn/sxsykc/jxnr/sy10.htm
讲解的很详细,最后的练习是在[-1,2]上求函数f(x)=-x^2+2x+0.5的最大值.二次函数求极值,有公式可以直接使用,结果是x=1时,f(x)=1.5.主要是学习利用遗传算法来逼近的方法.

主要步骤:
1. 设计编码规则,因为所有可行解都是在解空间内用二进制编码表示,因此必须有一种可行的映射规则,从解空间到一定长度的二进制字符串进行映射.

2. 设计适应度函数(fitness).求极值问题函数的设计是g(x)=f(x)-N(这个N可以是任意值).

3. 随机产生初始种群,计算适应度函数和各样本的被选择概率;

4. 设置迭代次数,交叉概率Pc,变异概率Pm,开始遗传算法演化.

4.1 判断多样化指标,如果多样性较低,直接跳出循环.(这是我看完其他资料后,加入的步骤,可以加速收敛过程,否则必须老老实实迭代指定的次数才会结束.唯一的问题是,初始种群随机产生的样本比较单一,进化一次就完成.不过这种情况基本不会出现,只是有理论上的可能,但是要注意,因为当这种情况发生时,有可能得到不正确的结果).

4.2 产生新一代种群,包含几个步骤:

4.2.1 选择:有各种策略可用,较常用的是轮盘赌法,具体计算方式是: 根据被选择概率数组P得到一个新的累加数组A,其中 (A(i)=\Sigma_0^i(P_i)).然后随机产生一个数r,如果r<A(d),则第d个样本被选中.另外一种策略是最简单的,扔掉适应度最差的m个样本,复制适应度最好的样本m次,这样可以保证种群样本数不变.这种策略也叫精英(elitism method)演化策略.实际中可以在轮盘赌的基础上结合精英演化策略来保证演化总是朝着适应度高的方向进行.(轮盘赌法是否有存在的必要?)

4.2.2 产生随机数rc,如果rc<Pc,则交叉演化(crossover),继续精英策略,最好的某几个样本不参与交叉,其余样本两两配对随机交叉(一般取60%的样本进行交叉,为什么这么取?)

4.2.3 产生随机数rm,如果rm<Pm,则变异演化(mutation),还是精英策略,最好的某些样本永不变异,其余的一定比例的样本在某个随机位置变异(这个比例不知怎么选择,我选择的10%,有没有什么要求?)

4.3 计算新种群的适应度数组

4.4 计算新种群的被选择概率数组,准备进入下次演化(回到4)

最终可以输出结果.

经测试,在精度为6位小数(字符串长度22),样本多样化参数(size(set(P))/size(P)设置为小于0.15时break循环的情况下的一些测试结果: (实际值x=1,f(x)=1.5)

样本6,50次迭代输出结果为x=1.0069,f(x)=1.4995.
样本10,100次迭代在第15次跳出,x=1.00067,f(x)=1.5.
样本30,100次迭代在第78次跳出,x=1(有2个样本的x是0.999994),f(x)=1.5

#include <iostream>
#include <string>
#include <vector>
#include <ctime>
#include <cmath> // ceil
#include <numeric> // accumulate
#include <map> // for quick sort to get new generation by probability
using namespace std;

// random number generator in [a,b]
double randRange(int a, int b)
{
 return (double)rand()/RAND_MAX*(b-a)+a;
}

// m^n
double powint(int m, int n)
{
 if (n==0)
 {
  return 1;
 }
 else
 {
  return m*powint(m,n-1);
 }
}

// from bin string to decimal number
double bin2dec(string binNum)
{
 double temp = 0.0;
 int numDigit;
 int len = (int)binNum.size();
 for (int i=0;i<len;++i)
 {
  numDigit = binNum.substr(len-1-i,1).compare("1")==0 ? 1 : 0;
  temp += numDigit*powint(2,i);
 }

 return temp;
}

// from decimal number to bin string
string dec2bin(int num, int len)
{
 string binNum = "";

 while (num>=2)
 {
  binNum.insert(0,num%2==0 ? "0" : "1");
  num /= 2;
 }

 binNum.insert(0,num==0 ? "0" : "1");

 while(binNum.size()<len)
 {
  binNum.insert(0,"0");
 }

 return binNum;
}

// double number in range [a,b] to represent val
double rangeValue(double val, int a, int b, int len)
{
 return a+val*(b-a)/(powint(2,len)-1);
}

// real value in [0, RAND_MAX] to represents val in [a,b]
int realValue(double val, int a, int b, int len)
{
 return floor((val-a)*(powint(2,len)-1)/(b-a));
}

// input the number of digits needed within range [a,b]
// 2 means 0.01 or 10^(-2)
int getLength(int a, int b, int n)
{
 int maxab, minab, nDigits=1;
 maxab = a>=b ? a : b;
 minab = a>=b ? b : a;

 if (maxab==minab)
 {
  cout << "It's a point or zero length range!" << endl;
   return 0;
 }

 double temp = (maxab-minab)*powint(10,n);
 while(temp>powint(2,nDigits)) nDigits++;

 return nDigits;
}

// generate num binary string vector with length=len represent the numbers in [a,b]
vector<string> randSample(int num, int len, int a, int b)
{
 vector<string> samples;

 for (int i=0;i<num;++i)
 {
  samples.push_back(dec2bin(realValue(randRange(a,b),a,b,len),len));
 }

 return samples;
}

// target function: f(x) = -x*x+2*x+0.5
double targetFunc(double x)
{
 return -x*x+2*x+0.5;
}

// adaptive values: g(x) = f(x)-Fmin
vector<double> Fitness(vector<string> samples, int a, int b, int len, double fmin=-2.0)
{
 vector<double> fitness;
 for (int i=0;i<(int)samples.size();++i)
 {
  fitness.push_back(targetFunc(rangeValue(bin2dec(samples[i]),a,b,len))-fmin);
 }

 return fitness;
}

// chosen probability, p(i) = g(i)/sum_of_g
vector<double> choseProbability(vector<double> fitness)
{
 vector<double> chosenProbs;
 double sum = accumulate(fitness.begin(),fitness.end(),0.0);

 for (int i=0;i<(int)fitness.size();++i)
 {
  chosenProbs.push_back(fitness[i]/sum);
 }

 //// found the smallest one
 //double minNum = fitness[0];
 //for (int i=1;i<(int)fitness.size();++i)
 //{
 // minNum = minNum>fitness[i] ? fitness[i] : minNum;
 //}

 //for (int i=0;i<(int)fitness.size();++i)
 //{
 // chosenProbs.push_back((fitness[i]-minNum)/sum);
 //}
 return chosenProbs;
}

// roulette wheel selection to choose one sample
// the probability of x_i to be chosen is probs[i]
// if r>sigma_0^n[probs(i)], the n-th sample is chosen
string rouletteSelection(vector<double> probs, vector<string> samples)
{
 double accSum = 0.0;

 // roulette choose
 double randChose = randRange(0,1);

 int i;

 // find the one
 for (i=0;i<(int)probs.size();++i)
 {
  accSum += probs[i];
  if (randChose<=accSum) { break; }
 }

 return samples[i];

}

// cross: only corssRate% samples will be usef for crossover
vector<string> Crossover(vector<string> samples, double crossRate=0.6)
{
 vector<string> newSamples;

 int nums = (int)samples.size();

 if (nums<=0)
 {
  cout << "Empty generation!\n";
  return newSamples;
 }
 else
 {
  // crossed number, must be even
  int crossNum = floor(nums*crossRate);
  crossNum = crossNum%2==0 ? crossNum : crossNum+1;

  // length of the string
  int len = (int)samples[0].size();

  // crossover: i and crossNum-1-i
  for (int i=0;i<crossNum/2;++i)
  {
   // where to start cross
   int crossPos = floor(randRange(0,len-1));

   // crossover
   if (0==samples[i].compare(samples[crossNum-1-i]))
   {
    newSamples.push_back(samples[i]);
    newSamples.push_back(samples[i]);
   }
   else
   {
    string firsthalf = samples[i].substr(0,crossPos);
    string secondhalf = samples[crossNum-1-i].substr(crossPos);
    newSamples.push_back(samples[i].substr(0,crossPos)+samples[crossNum-1-i].substr(crossPos));

    firsthalf = samples[crossNum-1-i].substr(0,crossPos);
    secondhalf = samples[i].substr(crossPos);
    newSamples.push_back(samples[crossNum-1-i].substr(0,crossPos)+samples[i].substr(crossPos));
   }
  }

  // add elite samples
  int ind = 0;
  while((int)newSamples.size()<nums)
  {
   newSamples.push_back(samples[nums-1-ind]);
   ind++;
  }

  return newSamples;
 }
}

// mutation: only the mutateRate% samples will be used for mutation
vector<string> Mutation(vector<string> samples, double mutateRate=0.1)
{
 vector<string> newSamples;

 int nums = (int)samples.size();

 if (nums<=0)
 {
  cout << "Empty generation!\n";
  return newSamples;
 }
 else
 {
  // number of samples that will mutate
  int mutateNum = ceil(nums*mutateRate);

  // length of the string
  int len = (int)samples[0].size();

  for (int i=0;i<nums;++i)
  {
   if (i<mutateNum)
   {
    // position
    int mutatePos = floor(randRange(0,len-1));
    mutatePos = mutatePos>=0 ? mutatePos : 0;
    mutatePos = mutatePos<(int)samples[0].size() ? mutatePos : (int)samples[0].size()-1;

    // mutate
    string first = samples[i].substr(0,mutatePos-1);
    string theOne = samples[i].substr(mutatePos,1);
    string second = samples[i].substr(mutatePos+1);
    theOne = theOne.compare("0")==0 ? "1" : "0";
    newSamples.push_back(first+theOne+second);
   }
   else
   {
    newSamples.push_back(samples[i]);
   }
  }

  return newSamples;
 }
}

// get new generation
// this function always drop the sample with lowest probability
// it should use roulette wheel selection method
// pc: corss probability, if random<pc, crossover
// pm: mutation probability, if random<pm, mutation
vector<string> newGeneration(vector<string> samples, vector<double> probs, double pc=0.6, double pm=0.01)
{
 // map double to string, use double values as keys
 // why? because the map will keep the keys sorted automatically,
 // it's easier to get the string indices to be changed at the same time
 map<double,string> probSamplePair;
 for (int i=0;i<(int)probs.size();++i)
 {
  probSamplePair[probs[i]] = samples[i];
 }

 // return value
 vector<string> newSamples;

 // associated probabilities
 vector<double> newProbs;

 // reproduce the samples
 map<double,string>::iterator it;
 for (it=probSamplePair.begin();it!=probSamplePair.end();it++)
 {
  //cout << (*it).first << "->" << (*it).second.c_str() << endl;
  newSamples.push_back((*it).second);
  newProbs.push_back((*it).first);
 }

 // replace the one with minimal probability with the one with the biggest
 // this is the genetic strategy which can be modified
 for (int i=0;i<(int)newSamples.size()-1;++i)
 {
  newSamples[i] = newSamples[i+1];
  newProbs[i] = newProbs[i+1];
 }

 if ((int)newSamples.size()<(int)samples.size())
 {
  string lastString = newSamples[(int)newSamples.size()-1];
  double lastDouble = newProbs[(int)newProbs.size()-1];

  // map will keep the keys identical, here to check the length
  while((int)newSamples.size()<(int)samples.size())
  {
   newSamples.push_back(lastString);
   newProbs.push_back(lastDouble);
  }
 }

 // crossover
 if (randRange(0,1)<pc)
 {
  newSamples = Crossover(newSamples);
 }

 // mutation
 if (randRange(0,1)<pm)
 {
  newSamples = Mutation(newSamples);
 }

 return newSamples;
}

// variation
double Variation(vector<string> samples)
{
 map<double,string> pairs;
 for (int i=0;i<(int)samples.size();++i)
 {
  pairs[bin2dec(samples[i])] = samples[i];
 }

 return (double)pairs.size()/(int)samples.size();
}

int main(int argc, char* argv[])
{
 srand((unsigned int)time(0));

 // lowerbound, upperbound, precision and number of samples
 int a = -1, b = 2, precision = 6, nums = 30;

 // number of string length
 int len = getLength(a,b,precision);

 // generate initial population
 vector<string> samples = randSample(nums,len,a,b);

 // fitness
 vector<double> fitness = Fitness(samples,a,b,len);

 // probability distribution
 vector<double> probs = choseProbability(fitness);

 // number of generation of the evolving process
 int numGeneration = 100;

 // current generation number
 int curGeneration = 0;

 // evolving
 while (curGeneration<numGeneration)
 {
  // generate new generation
  samples = newGeneration(samples,probs);

  // output
  cout << "The " << curGeneration << "-th generation: ";
  for (int i=0;i<(int)samples.size();++i)
  {
   cout << rangeValue(bin2dec(samples[i]),a,b,len) << " ";
  }
  cout << endl;

  // check variation
  if (Variation(samples)<0.15)
  {
   break;
  }

  // fitness
  fitness = Fitness(samples,a,b,len);

  // probability
  probs = choseProbability(fitness);

  curGeneration++;
 }

 // print results
 vector<double> doubleValue;
 for (int i=0;i<(int)samples.size();++i)
 {
  doubleValue.push_back(rangeValue(bin2dec(samples[i]),a,b,len));
 }

 cout << "\nThe result:\n";
 for (int i=0;i<(int)doubleValue.size();++i)
 {
  cout << "f(" << doubleValue[i] << ")" << " = " << targetFunc(doubleValue[i]) << "\n";
 }

 system("PAUSE");
 return 0;
}