遗传算法求解多元函数极值点-C++实现

遗传算法是根据大自然中生物体进化规律而设计提出的。是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。该算法通过数学的方式,利用计算机仿真运算,将问题的求解过程转换成类似生物进化中的染色体基因的交叉、变异等过程。在求解较为复杂的组合优化问题时,相对一些常规的优化算法,通常能够较快地获得较好的优化结果。

本文将通过一个简单的例子：求解的多元函数的极值点，初步理解遗传算法。

题目

函数f(x, y) = cos(x / 40) + cos(7x / 40) / 2 + cos(17x / 40) / 3 + sin(y / 40) + sin(7y / 40) / 2 + sin(17y / 40) / 3, x, y∈N.
求函数f(x, y)的最大值点.

解题思路

该二元函数所确定的空间曲面，类似于一片地势起伏的“山地”，求该二元函数的极值点的过程，相当于求这片山地的“最高峰”。
固然，我们可以用数学知识直接求得答案，但这只是一个简单的例子，如果情况更加复杂，繁冗的数学计算就显得相对困难。
在本例中，我们用遗传算法间接求解。假设在这片“山地”的不同海拔随机投放大量的“生物”，“生物”在这片“山地”繁衍生息。每隔一段时间，我们“猎杀”一些低海拔地区的“生物”，不难想到，高海拔地区的生物在“生存斗争”中将更容易存活下来。在经历足够多代的“进化”后，“生物”将聚集在山地的“相对高峰”。

建立模型

基因型 决定“生物”所处海拔的因素有2个：“生物”的x坐标值和y坐标值。把x坐标和y坐标抽象成“生物”的基因型，即“生物”有且只有x基因和y基因。易知该二元函数对x和y的最小正周期都为80π≈251<255，故可以用8位二进制数编码x基因和y基因：
x基因：00000000~11111111（对应十进制0 ~255）
y基因：00000000~11111111（对应十进制0 ~255）
“生物”的基因型：00000000 00000000~11111111 11111111
表现型 忽略环境因素的影响，基因型决定表现型。在本例中，“生物”的表现型即生物所处的海拔高度。
例：某“生物”个体的基因型为01101111 00110010（对应十进制x0=111，y0=50），则该生物的表现型为f(x0,y0)=1.36246。
环境适应度 衡量生物在一定环境条件下的生存能力，在本例中，可以简单的认为环境适应度与表现型成线性关系，即“生物”所处海拔越高，越容易生存。
基因重组 父本和母本产生后代的过程中，不同基因随机组合的过程。注意：基因重组的对象为“基因”，故不应该把基因片段的一部分“重组”。
例：亲本基因型：01010101 000011111x00001111 10101111
子代基因型：01010101 10101111（随机组合亲本的不同基因，是基因重组）
子代基因型：01011111 00001111（把亲本的基因拆开组合，不是基因重组）
基因突变 小概率使“生物”个体的基因序列的某一位或几位发生改变（0变成1,1变成0），不考虑染色体变异。
例：某“生物”个体的基因型为01010101 000011111
基因突变后该个体的基因型为11010101 000011111（第一位基因发生了改变）

※自然选择 遗传算法的核心，是“猎杀”低环境适应度个体，同时尽量让高环境适应度个体产生后代的过程。对于我们预先投放的“生物”种群，每隔一段时间计算一次种群平均环境适应度，随机选出一定量的“生物”个体，如果该个体的环境适应度低于种群平均环境适应度，就“猎杀”该个体。同时随机选择两个不同的个体交配，产生后代，维持“生物”种群的平衡。为了使环境适应度高的个体更容易获得交配的机会，我们将使用轮盘赌算法选择获得交配权的个体。本文对轮盘赌算法不做赘述，有兴趣的读者可以自行查询。

进化在本例中，一次完整的自然选择、种群信息统计、种群信息输出过程，称为进化。一次进化后，新一代的“生物”种群将拥有更好的基因和环境适应度。

遗传算法求最优解的可行性

局部最优解 即小范围的最优解，从函数图像中，可以看出有很多“小山峰”，每一个小山峰都是一个局部最优解。遗传算法通过基因重组实现优良基因个体对劣质基因个体的替换，保证了求解局部最优解的可行性。
全局最优解 即所有局部最优解的最优解。在函数图像中，全局最优解对应“最高的山峰”。遗传算法通过基因突变实现在全局空间内的大幅跳跃，保证了求解全局最优解的可行性。

代码实现

类和函数声明

struct Genotype;     //基因型
struct Individual;      //个体
//初始化种群
void init();
//统计种群信息
void statistics();
//输出种群信息
void print();
//自然选择函数：返回值：0-正常结束，1-求得局部最优解，选择非正常结束
bool select(int replace, int _variate = 0);
//基因重组：返回子代个体
Individual recombination(Individual const&mather, Individual const&father);
//变异函数：使输入个体基因突变
void variate(Individual&individual);
//进化函数：进化max_generation代
void evolve(int max_generation = 1);
//退出程序
void quit();

宏定义

//定义种群规模的宏
#define MAX_INDIVIDUALS 1000
//定义目标二元函数的宏
#define F(x,y) (cos(x/40)+cos(7*x/40)/2+cos(17*x/40)/3+sin(y/40)+sin(7*y/40)/2+sin(17*y/40)/3)

全局变量

Individual*population;       //种群
int generation;             //当前进化代数
unsigned aver_adaptability; //种群平均环境适应度
double aver_phenotype;      //种群平均表现型
Individual*best;            //最优个体

定义基因型

//_GENOTYPE_BEGINstruct Genotype {unsigned x : 8;    //x基因unsigned y : 8;    //y基因Genotype() = default;Genotype(unsigned geno);operator unsigned()const;
};Genotype::Genotype(unsigned geno) {x = geno >> 8;y = geno ^ (x << 8);
}Genotype::operator unsigned()const {unsigned gene = x;gene <<= 8;gene |= y;return gene;
}ostream&operator<<(ostream&output, Genotype const&genotype) {unsigned gene, ribosome = 1 << 7;gene = genotype.x;for (int count = 0; count != 8; count++) {output << (gene&ribosome ? 1 : 0);gene <<= 1;}output.put(32);gene = genotype.y;for (int count = 0; count != 8; count++) {output << (gene&ribosome ? 1 : 0);gene <<= 1;}return output;
}//_GENOTYPE_END

定义生物个体

//_INDIVIDUAL_BEGINstruct Individual {Genotype genotype;     //基因型double phenotype;      //表现型unsigned adaptability; //环境适应度void set();              //根据基因型，设置表现型和环境适应度
};void Individual::set() {unsigned x, y;x = genotype.x;y = genotype.y;phenotype = F(x, y);if (phenotype < 0)phenotype = 0;adaptability = static_cast<unsigned>(100 * phenotype);
}ostream&operator<<(ostream&output, Individual const&individual) {output << individual.genotype;output.put(32);output << setiosflags(ios::fixed | ios::left) <<setprecision(2) << setw(4) << individual.phenotype;output.put(32);output << individual.adaptability;return output;
}//_INDIVIDUAL_END

初始化种群

void init() {population = new Individual[MAX_INDIVIDUALS];best = new Individual;srand(static_cast<unsigned>(time(NULL)));auto it = population;for (int count = 0; count != MAX_INDIVIDUALS; count++) {it->genotype.x = rand() % 256;it->genotype.y = rand() % 256;it++->set();}statistics();
}

统计种群信息

void statistics() {double aver[2] = {};auto it = population;for (int count = 0; count != MAX_INDIVIDUALS; count++) {aver[0] = (aver[0] * count + it->adaptability) / (count + 1);aver[1] = (aver[1] * count + it->phenotype) / (count + 1);it++;}aver_adaptability = static_cast<unsigned>(aver[0]);aver_phenotype = aver[1];
}

输出种群信息

void print() {cout << "[Generation:" << generation << ']' << endl;cout << "[Aver:Adaptability:" << aver_adaptability << ']' << endl;cout << "[Aver:Phenotype:" << setiosflags(ios::fixed) << setprecision(2) <<aver_phenotype << ']' << endl;
}

自然选择：轮盘赌算法实现

bool select(int replace, int _variate) {//创建整型数组：储存亲本下标int*m_subscript = new int[replace];int*f_subscript = new int[replace];//选择母本和父本：环境适应性决定unsigned long long SUM = MAX_INDIVIDUALS * aver_adaptability;auto it = population;for (int i = 0; i != replace;) {//根据总环境适应性，创建轮盘unsigned long long m_path, f_path;m_path = SUM * rand() / RAND_MAX;do f_path = SUM * rand() / RAND_MAX;while (m_path == f_path);   //母本和父本不能相同//计算母本下标int subscript[2];it = population;for (int j = 0; j != MAX_INDIVIDUALS; j++, it++)if (m_path <= it->adaptability) {subscript[0] = j;break;}else m_path -= it->adaptability;//计算父本下标it = population;for (int j = 0; j != MAX_INDIVIDUALS; j++, it++)if (f_path <= it->adaptability) {subscript[1] = j;break;}else f_path -= it->adaptability;bool flag = false;for (int j = 0; j != i; j++) {if (subscript[0] == m_subscript[j] || subscript[0] == f_subscript[j] ||subscript[1] == f_subscript[j] || subscript[1] == f_subscript[j]) {flag = true;break;}}if (flag)continue;else {m_subscript[i] = subscript[0];f_subscript[i] = subscript[1];i++;}}//淘汰环境适应性低于平均值的个体it = population;int succeed = 1000;    //标识符：进化成功for (int count = 0; count != replace && succeed;) {int temp = rand() % MAX_INDIVIDUALS;if (it[temp].adaptability < aver_adaptability) {it[temp] = recombination(it[m_subscript[count]], it[f_subscript[count]]);count++;succeed = 1000;}else succeed--;//如果连续1000次找不到环境适应性低于平均值的个体//则说明该种群的个体已经高度相似，进化提前完成}//变异：环境适应性低于平均值的个体发生基因突变if(succeed)for (int count = 0; count != _variate && succeed;) {int temp = rand() % MAX_INDIVIDUALS;if (it[temp].adaptability < aver_adaptability) {variate(it[temp]);count++;succeed = 1000;//同理，如果连续1000次找不到环境适应性低于平均值的个体，进化提前完成}else succeed--;}return succeed ? false : true;delete[]m_subscript;delete[]f_subscript;
}

基因重组：x基因和y基因的自由组合

Individual recombination(Individual const&mather, Individual const&father) {unsigned gene;gene = rand() % 2 ? mather.genotype.x : father.genotype.x;gene <<= 8;gene |= rand() % 2 ? mather.genotype.y : father.genotype.y;Individual child;child.genotype = Genotype(gene);child.set();return child;
}

变异：改变基因的单个结点

void variate(Individual&individual) {int position = rand() % 16;    //基因突变的位置unsigned temp;if (position < 8) {temp = 1 << position;individual.genotype.y ^= temp;}else {position -= 8;temp = 1 << position;individual.genotype.x ^= temp;}
}

进化

void evolve(int max_generation) {init(); //初始化种群信息//最多进行max_generation代进化for (int count = 0; count != max_generation; count++) {statistics();    //统计种群信息print();        //输出种群信息if (select(50))break;   //更新50个个体，如果进化提前完成，则退出循环generation++;   //进化的代数+1}//进化结束：筛选最优个体auto it = population;*best = *population;for (int count = 0; count != MAX_INDIVIDUALS; count++) {if (it->adaptability > best->adaptability)*best = *it;it++;}//输出最优个体信息：遗传算法结束cout << "\n[Best]" << endl;cout << "[x:" << best->genotype.x << "]:[y:" << best->genotype.y << ']' << endl;cout << *best << endl;
}

退出程序：释放申请的内存空间

void quit() {delete[]population;delete best;
}

结果与分析

创建一个规模为1000的种群，最大进化1000代，观察输出结果。

int main() {evolve(1000);quit();return 0;
}

输出结果：
初始种群信息：

进化完成后种群信息：

分析：
种群在第75代就完成了进化，平均环境适应度提高了约6倍，最优个体的基因型为x=2，y=48。
结论：
函数f(x, y) 的最大值点为(2, 48)，最大值为f(2, 48)=3.47。