遗传算法的C语言实现(一):以非线性函数求极值为例
以前搞数学建模的时候，研究过(其实也不算是研究，只是大概了解)一些人工智能算法，比如前面已经说过的粒子群算法(PSO),还有著名的遗传算法(GA),模拟退火算法(SA),蚁群算法(ACA）等。当时懂得非常浅，只会copy别人的代码(一般是MATLAB),改一改值和参数，东拼西凑就拿过来用了，根本没有搞懂的其内在的原理到底是什么。这一段时间，我又重新翻了一下当时买的那本《MATLAB智能算法30个案例分析》，重读一遍，发现这本书讲的还是非常不错的，不仅有现成算法的MATLAB实现，而且把每一种算法的实现原理都基本讲清楚了，我再次看的时候，以前不懂的(也许以前太急功近利了)现在居然基本全部弄懂了。真是让我感到非常开心~~~所以我萌发了一个冲动，想用现在自己比较熟悉的C语言，把这本书上的案例重新实现一遍，彻底弄懂每个算法的原理。我打算不定期地更新几篇来介绍其中几种常见的智能算法和其应用(当然我也不止会看这本书，还会参考别的书籍和论文),尽量把每种算法的原理给讲明白，也算是对我接触这些智能算法的一个总结吧。
今天首先介绍遗传算法(genetic algorithm,GA)。
遗传算法是一种进化算法，其基本原理是模仿自然界中的生物“物竞天择，适者生存”的进化法则，把问题参数编码为染色体，再利用迭代的方式进行选择、交叉、变异等运算法则来交换种群中染色体的信息，最终生成符合优化目标的染色体。
在遗传算法中，染色体对应的是数据或者数组，通常是由一维的串结构数据来表示，串上各个位置对应基因的的取值。基因组成的串就是染色体，或者称为基因型个体。一定数量的个体组成了种群，种群中个体数目的大小称为种群大小，也称为种群规模，而各个个体对环境的适应程度称为适应度。
标准遗传算法的步骤如下：
(1) 编码：遗传算法在搜索解空间之前需要将解数据表示成遗传空间的基因型串结构数据，这些串结构数据的不同组合构成了不同的染色体。
(2)初始化：即生成初始种群。具体的做法是随机生成N个初始的染色体(解空间的解),每一个染色体其实就相当于一个个体，N个个体构成了一个初始种群。遗传算法以这N个个体作为初始值开始进化。
(3) 适应度评估:适应度表明个体或者解的优劣性。不同的问题，定义适应度函数的方式也不同。比如如果是求一个函数的最大值，那么一般某个解对应的函数值越大，那么这个个体(解)的适应度也就越高。
（4）选择:选择的目的是为了从当前种群中选出优良的个体，使它们有机会成为父代繁殖下一代子孙。遗传算法通过选择过程体现这一思想，进行选择的原则是适应性强的个体为下一代贡献一个或者多个后代的概率大。这体现了达尔文的适者生存原则。
(5)交叉：交叉操作是遗传算法中最主要的遗传操作。通过交叉可以得到新一代个体，新个体组合了其父代的个体特征。交叉体现了信息交换的思想。
(6)变异:变异首先在群体中随机选择一个个体，对于选中的个体以一定的概率（通常是比较小的概率，这与自然界一致，自然界的变异都是小概率事件）随机改变染色体中某个基因的值。
有的时候除了选择选择、交叉、变异这三种操作之外，我们还会针对具体的问题加入其它的操作（比如逆转之类），但是选择、交叉、变异是所有的遗传算法都共同的拥有的遗传操作。
本文以遗传算法常见的一个应用领域——求解复杂的非线性函数极值问题为例，来说明如何用代码(这里是用C语言，当然你可以换做别的任何一种语言)来具体实现上述的遗传算法操作。求解的函数来自《MATLAB智能算法30个案例分析》，即：
f = -5*sin(x1)*sin(x2)*sin(x3)*sin(x4)*sin(x5) - sin(5*x1)*sin(5*x2)*sin(5*x3)*sin(5*x4)*sin(5*x5) + 8
其中x1,x2,x3,x4,x5是0~0.9*PI之间的实数。该函数的最小值为2,当x1,x2,x3,x4,x5都取PI/2时得到。
使用C语言实现的遗传算法求解代码如下：
/*
* 遗传算法求解函数的最优值问题
* 参考自《MATLAB智能算法30个案例分析》
* 本例的寻优函数为:f = -5*sin(x1)*sin(x2)*sin(x3)*sin(x4)*sin(x5) - sin(5*x1)*sin(5*x2)*sin(5*x3)*sin(5*x4)*sin(5*x5) + 8
* 其中x1,x2,x3,x4,x5是0~0.9*PI之间的实数。该函数的最小值为2,当x1,x2,x3,x4,x5都取PI/2时得到。
* update in 16/12/3
* author:Lyrichu
* email:919987476@qq.com
*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#define PI 3.1415926 //圆周率
#define sizepop 50 // 种群数目
#define maxgen 500 // 进化代数
#define pcross 0.6 // 交叉概率
#define pmutation 0.1 // 变异概率
#define lenchrom 5 // 染色体长度
#define bound_down 0 // 变量下界，这里因为都相同，所以就用单个值去代替了，如果每个变量上下界不同，也许需要用数组定义
#define bound_up (0.9*3.1415926) // 上界
double chrom[sizepop][lenchrom]; // 种群数组
double fitness[sizepop]; //种群每个个体的适应度
double fitness_prob[sizepop]; // 每个个体被选中的概率(使用轮盘赌法确定)
double bestfitness[maxgen]; // 每一代最优值
double gbest_pos[lenchrom]; // 取最优值的位置
double average_best[maxgen+1]; // 每一代平均最优值
double gbest; // 所有进化中的最优值
int gbest_index; // 取得最优值的迭代次数序号
// 目标函数
double fit_func(double * arr)
{
double x1 = *arr;
double x2 = *(arr+1);
double x3 = *(arr+2);
double x4 = *(arr+3);
double x5 = *(arr+4);
double func_result = -5*sin(x1)*sin(x2)*sin(x3)*sin(x4)*sin(x5) - sin(5*x1)*sin(5*x2)*sin(5*x3)*sin(5*x4)*sin(5*x5) + 8;
return func_result;
}
// 求和函数
double sum(double * fitness)
{
double sum_fit = 0;
for(int i=0;i<sizepop;i++)
sum_fit += *(fitness+i);
return sum_fit;
}
// 求最小值函数
double * min(double * fitness)
{
double min_fit = *fitness;
double min_index = 0;
static double arr[2];
for(int i=1;i<sizepop;i++)
{
if(*(fitness+i) < min_fit)
{
min_fit = *(fitness+i);
min_index = i;
}
}
arr[0] = min_index;
arr[1] = min_fit;
return arr;
}
// 种群初始化
void init_chrom()
{
for(int i=0;i<sizepop;i++)
{
for(int j=0;j<lenchrom;j++)
{
chrom[i][j] = bound_up*(((double)rand())/RAND_MAX);
}
fitness[i] = fit_func(chrom[i]); // 初始化适应度
}
}
// 选择操作
void Select(double chrom[sizepop][lenchrom])
{
int index[sizepop];
for(int i=0;i<sizepop;i++)
{
double * arr = chrom[i];
fitness[i] = 1/(fit_func(arr)); // 本例是求最小值，适应度为目标函数的倒数，即函数值越小，适应度越大
}
double sum_fitness = 0;
for(int i=0;i<sizepop;i++)
{
sum_fitness += fitness[i]; // 适应度求和
}
for(int i=0;i<sizepop;i++)
{
fitness_prob[i] = fitness[i]/sum_fitness;
}
for(int i=0;i<sizepop;i++)
{
fitness[i] = 1/fitness[i]; // 恢复函数值
}
for(int i=0;i<sizepop;i++) // sizepop 次轮盘赌
{
double pick = ((double)rand())/RAND_MAX;
while(pick < 0.0001)
pick = ((double)rand())/RAND_MAX;
for(int j=0;j<sizepop;j++)
{
pick = pick - fitness_prob[j];
if(pick<=0)
{
index[i] = j;
break;
}
}
}
for(int i=0;i<sizepop;i++)
{
for(int j=0;j<lenchrom;j++)
{
chrom[i][j] = chrom[index[i]][j];
}
fitness[i] = fitness[index[i]];
}
}
// 交叉操作
void Cross(double chrom[sizepop][lenchrom])
{
for(int i=0;i<sizepop;i++)
{
// 随机选择两个染色体进行交叉
double pick1 = ((double)rand())/RAND_MAX;
double pick2 = ((double)rand())/RAND_MAX;
int choice1 = (int)(pick1*sizepop);// 第一个随机选取的染色体序号
int choice2 = (int)(pick2*sizepop);// 第二个染色体序号
while(choice1 > sizepop-1)
{
pick1 = ((double)rand())/RAND_MAX;
choice1 = (int)(pick1*sizepop); //防止选取位置过界
}
while(choice2 > sizepop-1)
{
pick2 = ((double)rand())/RAND_MAX;
choice2 = (int)(pick2*sizepop);
}
double pick = ((double)rand())/RAND_MAX;// 用于判断是否进行交叉操作
if(pick>pcross)
continue;
int flag = 0; // 判断交叉是否有效的标记
while(flag == 0)
{
double pick = ((double)rand())/RAND_MAX;
int pos = (int)(pick*lenchrom);
while(pos > lenchrom-1)
{
double pick = ((double)rand())/RAND_MAX;
int pos = (int)(pick*lenchrom); // 处理越界
}
// 交叉开始
double r = ((double)rand())/RAND_MAX;
double v1 = chrom[choice1][pos];
double v2 = chrom[choice2][pos];
chrom[choice1][pos] = r*v2 + (1-r)*v1;
chrom[choice2][pos] = r*v1 + (1-r)*v2; // 交叉结束
if(chrom[choice1][pos] >=bound_down && chrom[choice1][pos]<=bound_up && chrom[choice2][pos] >=bound_down && chrom[choice2][pos]<=bound_up)
flag = 1; // 交叉有效，退出交叉，否则继续下一次交叉
}
}
}
// 变异操作
void Mutation(double chrom[sizepop][lenchrom])
{
for(int i=0;i<sizepop;i++)
{
double pick = ((double)rand())/RAND_MAX; // 随机选择一个染色体进行变异
int choice = (int)(pick*sizepop);
while(choice > sizepop-1)
{
pick = ((double)rand())/RAND_MAX;
int choice = (int)(pick*sizepop); // 处理下标越界
}
pick = ((double)rand())/RAND_MAX; // 用于决定该轮是否进行变异
if(pick>pmutation)
continue;
pick = ((double)rand())/RAND_MAX;
int pos = (int)(pick*lenchrom);
while(pos > lenchrom-1)
{
pick = ((double)rand())/RAND_MAX;
pos = (int)(pick*lenchrom);
}
double v = chrom[i][pos];
double v1 = v - bound_up;
double v2 = bound_down - v;
double r = ((double)rand())/RAND_MAX;
double r1 = ((double)rand())/RAND_MAX;
if(r >= 0.5)
chrom[i][pos] = v - v1*r1*(1-((double)i)/maxgen)*(1-((double)i)/maxgen);
else
chrom[i][pos] = v + v2*r1*(1-((double)i)/maxgen)*(1-((double)i)/maxgen); // 注意这里写法是不会越界的，所以只用一次变异就可以了
}
}
// 主函数
int main(void)
{
time_t start,finish;
start = clock(); // 程序开始计时
srand((unsigned)time(NULL)); // 初始化随机数种子
init_chrom(); // 初始化种群
double * best_fit_index = min(fitness);
int best_index =(int)(*best_fit_index);
gbest = *(best_fit_index+1); // 最优值
gbest_index = 0;
average_best[0] = sum(fitness)/sizepop; //初始平均最优值
for(int i=0;i<lenchrom;i++)
gbest_pos[i] = chrom[best_index][i];
// 进化开始
for(int i=0;i<maxgen;i++)
{
Select(chrom); // 选择
Cross(chrom); //交叉
Mutation(chrom); //变异
for(int j=0;j<sizepop;j++)
{
fitness[j] = fit_func(chrom[j]);
}
double sum_fit = sum(fitness);
average_best[i+1] = sum_fit/sizepop; // 平均最优值
double * arr = min(fitness);
double new_best = *(arr+1);
int new_best_index = (int)(*arr); // 新最优值序号
if(new_best < gbest)
{
gbest = new_best;
for(int j=0;j<lenchrom;j++)
{
gbest_pos[j] = chrom[new_best_index][j];
}
gbest_index = i+1;
}
}
finish = clock(); // 程序计算结束
double duration = ((double)(finish-start))/CLOCKS_PER_SEC;
printf("程序计算耗时:%lf秒\n.",duration);
printf("遗传算法进化了%d次，最优值为:%lf,最优值在第%d代取得,此代的平均最优值为%lf.\n",maxgen,gbest,gbest_index,average_best[gbest_index]);
printf("取得最优值的地方为(%lf,%lf,%lf,%lf,%lf).\n",gbest_pos[0],gbest_pos[1],gbest_pos[2],gbest_pos[3],gbest_pos[4]);
return 0;
}
说明：
（1）这里是求函数的最小值，函数值越小的个体适应度越大，所以这里采用倒数变换，即把适应度函数F(x)定义为1/f(x)(f(x)为目标函数值),这样就保证了适应度函数F(x)越大，个体是越优的。(当然我们已经保证了f（x）始终是一个正数)
(2) 选择操作采用的是经典的轮盘赌法，即每个个体被选为父代的概率与其适应度是成正比的，适应度越高，被选中的概率越大，设pi为第i个个体被选中的概率，则
pi=Fi/(∑Fi),Fi为第i个个体的适应度函数。
(3)交叉操作：由于个体采用实数编码，所以交叉操作也采用实数交叉法，第k个染色体ak和第l个染色体al在第j位的交叉操作方法为:
akj=aij(1-b)+aijb; alj=alj(1-b)+akjb,其中b是[0,1]之间的随机数。k，l，j都是随机产生的，即随机选取两个父代，再随机选取一个基因，按照公式完成交叉操作。
(4)变异操作：第i个个体的第j个基因进行变异的操作方法是：
r>=0.5时，aij=aij+(amax-aij)*f(g)
r<0.5时，aij=aij+(amin-aij)*f(g)
其中amax是基因aij的上界，amin是aij的下界，f(g)=r2*(1-g/Gmax）2，r2是一个[0,1]之间的随机数，g是当前迭代次数，Gmax是进化最大次数，r是[0,1]之间的随机数。
可以看出，当r>=0.5时，aij是在aij到amax之间变化的；r<0.5时，aij是在amin到aij之间变化的。这样的函数其实可以构造不止一种，来完成变异的操作。只要符合随机性，以及保证基因在有效的变化范围内即可。
我在ubuntu16.04下，使用gcc编译器运行得到的结果如下：
从结果上来看，种群数目为50，进化代数为500时，得到的最优解为2.014102已经非常接近实际的最优解2了，当然，如果把种群数目再变大比如500，进化代数变多比如1000代，那么得到的结果将会更精确，种群数目为500，进化代数为1000时的最优解如下：
可以看出，增加种群数目以及增加进化代数之后，最优值从2.014102变为2.000189,比原来的精度提高了两位，但是相应地，程序的运行时间也显著增加了(从0.08秒左右增加到2秒左右)。所以在具体求解复杂问题的时候(比如有的时候会碰到函数参数有几十个之多),我们就必须综合考虑最优解的精度以及运算复杂度(运算时间)这两个方面，权衡之后，取合适的参数进行问题的求解。
当然，这里只是以求解多元非线性函数的极值问题为例来说明遗传算法的具体实现步骤，遗传算法的实际运用非常广泛，不仅仅可以求解复杂函数极值，而且在运筹学等众多领域有着非常广泛的应用，也经常会和其他的智能算法或者优化算法结合来求解问题。这个后面我还会再较为详细地论述。
