遗传算法的并行实现.doc

北京工业大学并行计算课程设计（论文）

遗
传
算
法
（基于遗传算法求函数最大值）

指导老师：刘建丽
学 号：S201007156
姓 名：杨平
班 级：研10级1班

遗传算法
一、 遗传算法的基本描述
遗传算法（Genetic Algorithm，GA）是通过模拟自然界生物进化过程来求解优化问题的一类自组织、自适应的人工智能技术。它主要基于达尔文的自然进化论和孟德尔的遗传变异理论。多数遗传算法的应用是处理一个由许多个体组成的群体，其中每个个体表示问题的一个潜在解。对个体存在一个评估函数来评判其对环境的适应度。为反映适者生存的思想，算法中设计一个选择机制，使得：适应度好的个体有更多的机会生存。在种群的进化过程中，主要存在两种类型的遗传算子：杂交和变异。这些算子作用于个体对应的染色体，产生新的染色体，从而构成下一代种群中的个体。该过程不断进行，直到找到满足精度要求的解，或者达到设定的进化代数。显然，这样的思想适合于现实世界中的一大类问题，因而具有广泛的应用价值。遗传算法的每一次进化过程中的，各个体之间的操作大多可以并列进行，因此，一个非常自然的想法就是将遗传算法并行化，以提高计算速度。本报告中试图得到一个并行遗传算法的框架，并考察并行化之后的一些特性。为简单起见（本来应该考虑更复杂的问题，如TSP。因时间有些紧张，做如TSP等复杂问题怕时间不够，做不出来，请老师原谅），考虑的具有问题是：对给定的正整数n、n元函数f，以及定义域D，求函数f在D内的最大值。
二、 串行遗传算法
1． 染色体与适应度函数
对函数优化问题，一个潜在的解就是定义域D中的一个点，因此，我们只需用一个长度为n的实数数组来表示一个个体的染色体。由于问题中要求求函数f的最大值，我们可以以个体所代表点在f函数下的值来判断该个体的好坏。因此，我们直接用函数f作为个体的适应度函数。
2． 选择机制
选择是遗传算法中最主要的机制，也是影响遗传算法性能最主要的因素。若选择过程中适应度好的个体生存的概率过大，会造成几个较好的可行解迅速占据种群，从而收敛于局部最优解；反之，若适应度对生存概率的影响过小，则会使算法呈现出纯粹的随机徘徊行为，算法无法收敛。下面我们介绍在实验中所使用的选择机制。
我们定义为当前种群内所有个体的集合，为中所有个体的一个固定排列。若为某一个体，表示该个体的适应度，则种群的适应度定义为：

对任意个体，的相对适应度定义为。相对适应度反映了个体的适应度在整个适应度总和中所占的比例。个体适应度越高，被选中的概率越高。累积适应度定义为:

进行选择之前，先产生一个0到1之间的随机实数，若满足，则第k+1个个体被选中。循环以上过程，即得到生成下一代种群的母体。
具体实现见如下函数：
void pop_select(void)
{
int mem, i, j, k;
double sum = 0;
double p;

/* 计算种群适应度之和 */
for (mem = 0; mem < POPSIZE; mem++) {

/* 按照累积适应度概率选取母体种群 */
for (i = 0; i < POPSIZE; i++) {
p = rand()%1000/1000.0;
if (p < population[0].cfitness)
newpopulation[i] = population[0];
else {
for (j = 0; j < POPSIZE;j++)
if (p >= population[j].cfitness &&
p < population[j+1].cfitness)
newpopulation[i] = population[j+1];
}
}
/*计算种群的总适应度*/
for (i = 0; i < POPSIZE; i++)
population[i] = newpopulation[i];
}
sum += (population[mem].fitness - lower_fitness);
}
/* 计算相对适应度 */
for (mem = 0; mem < POPSIZE; mem++) {
population[mem].rfitness = (population[mem].fitness - lower_fitness)/sum;
}
population[0].cfitness = population[0].rfitness;

/* 计算累积适应度 */
for (mem = 1; mem < POPSIZE; mem++) {
population[mem].cfitness = population[mem-1].cfitness +
population[mem].rfitness;
}
3． 杂交算子
杂交算子的流程一般如下：
（1） 按杂交概率选择一对参与进化的个体；
（2） 随机确定一个截断点；
（3） 将两个个体的染色体从截断点处截断，并交换，从而得到新的染色体。
具体算法见如下函数：
void crossover(void)
{
int i, j, k, m, point;
int first = 0;
double x;
for (k = 0; k < POPSIZE; k++) {
x = rand()%1000/1000.0; //产生随机交叉概率
if (x < PXOVER) /*如果随机交叉概率小于交叉概率，则进行交叉*/
{
first++;
if (first % 2 == 0) {
if (NVARS == 2) point = 1; //得到一个交叉点
else point = (rand() % (NVARS - 1)) + 1;

for (j = 0; j < point; j++)
//交叉运算，两个个体的交叉点前的基因进行交换
swap(&population[m].gene[j], &population[k].gene[j]);
}
else m = k;
}
}
}
4． 变异算子
在遗传算法中使用变异算子有两个目的：改善遗传算法的局部搜索能力。维持群体的多样性，防止出现早熟现象。变异操作的实现相当简单，只需遍历各染色体的各个单元，按某一变异概率将该单元变成一个随机的合法值。
其执行过程是：
（1）对个体的每一个基因组，依变异概率Pm指定为变异点。
（2）对每一个指定的变异点，对其基因取非或者用其他等位基因值来代替，从而产生一个新的个体。实现代码如下：
void mutate(void)
{
int i, j;
double lbound, hbound;
double p; //定义p为随机变异概率

for (i = 0; i < POPSIZE; i++)
for (j = 0; j < NVARS; j++) {
p = rand()%1000/1000.0;
if (p < PMUTATION) {
population[i].gene[j] = randval(lower[j], upper[j]); }
}
}
串行遗传算法的主要流程如图1所示。在每一次进化过程中，总是找出种群中的最优解与最差解，并将最优解保存，将本次最差解用上次保存的最优解替换，这样保证了各次进化的最优解的适应度不会降低，从而增快收敛的速度。

图1 串行遗传算法基本流程

三、 算法设计
分析图1中的串行算法，容易看出，在选择函数中，计算相对适应度需要用到全局种群的适应度之和，计算个体xk+1的累积适应度依赖于xk的累积适应度，如果在并行算法中要原封不动地模拟串行算法的运算，这些数据依赖关系都将产生通讯。更为不幸的是，选择后的个体需在各进程中作大量数据迁移。杂交算子中，一次杂交需要用到母体中的两个个体，若在这两个个体分配在不同进程，则需要进行一次通讯。此后的变异和评估都可以非常容易的实现并行，并且完全不需要任何通讯。但最后一步求最优个体和最差个体需要对各进程进行归约。由这些分析可以看出，完全地模拟串行情形将使算法变得相当低效。
幸运地是，遗传算法本身是一个概率算法，我们完全可以对串行算法作些必要的改变。如图2所示，我们将整个种群分成p个子种群，每一子种群由一个单一的进程负责。各进程独立地完成串行遗传算法的整个过程，唯一不同的是选择函数。各进程作选择操作时，首先计算各子种群内的局部累积适应度，然后根据局部累积适应度选择若干（本算法实现中使用的是常数3，也可以设为子种群大小的一个函数）个体按一固定规则轮流发送到其他进程；同时，按照该规则相应地从其他进程获取若干用来进行交流的个体。获取到个体后，先将其暂存；然后按串行算法中的选择机制从原子种群中选择进行进化的母体；最后再用之前暂存的个体完成进程间的种群交流。对每一个待交流的个体，具体策略如下：
（1） 随机地从本地的待进化母体种群内抽取与之进行交流的母体；
（2） 比较本地个体与传送过来的待交流个体，选取适应度高者作为最终母体。
各进程在每一次进化过程中，均分别保留各自的局部最优解，用来在下一次进化中替换局部最差的个体。各进程均完成所预定的进化迭代后，最后对各进程的局部最优解进行归约，从而得到整个算法的全局最优解。算法的主要流程详见图2。

图2 并行遗传算法基本流程

四、 算法实现
该算法实现的最关键部分为选择中的种群交流，该功能有如下函数实现

void pop_select(void)
{
MPI_Status status;
MPI_Request handle;
int mem, i, j, k;
double sum = 0;
double p;
static struct genotype ex_member[EX_NUM];

/* 计算子种群的总适应度 */
for (mem = 0; mem < TASK_NUM(pid); mem++) {
sum += (population[mem].fitness - lower_fitness);
}

/* 计算各个体相应适应度 */
for (mem = 0; mem < TASK_NUM(pid); mem++) {
population[mem].rfitness =(population[mem].fitness - lower_fitness)/sum;
}
population[0].cfitness = population[0].rfitness;

/* 计算各个体累积适应度 */
for (mem = 1; mem < TASK_NUM(pid); mem++) {
population[mem].cfitness=population[mem-1].cfitness+population[mem].rfitness;
}

/* 按照累积适应度概率选取种群交流个体，并发送和接收 */
for (i = 1; i <= EX_NUM; i++) {
p = rand()%1000/1000.0;
if (p < population[0].cfitness) {
MPI_Isend(&population[0], sizeof(struct genotype)/sizeof(char),
MPI_CHAR, (pid+i*generation)%pnum, 0, MPI_COMM_WORLD, &handle);
}
else {
for (j = 0; j < TASK_NUM(pid);j++)
if (p >= population[j].cfitness && p < population[j+1].cfitness){
MPI_Isend(&population[j+1], sizeof(struct genotype)/sizeof(char), MPI_CHAR, (pid+i*generation)%pnum, 0, MPI_COMM_WORLD, &handle);
break;
}
}

MPI_Recv (&ex_member[i-1], sizeof(struct genotype)/sizeof(char), MPI_CHAR, (pid+(pnum-i)*generation)%pnum, 0, MPI_COMM_WORLD, &status);
}

/* 按照累积适应度概率选取母体种群 */
for (i = 0; i < TASK_NUM(pid); i++)
{
p = rand()%1000/1000.0;
if (p < population[0].cfitness)
newpopulation[i] = population[0];
else {
for (j = 0; j < TASK_NUM(pid); j++)
if (p >= population[j].cfitness &&
p < population[j+1].cfitness)
newpopulation[i] = population[j+1];
}
}

for (i = 0; i < TASK_NUM(pid); i++)
population[i] = newpopulation[i];

/* 按优胜劣汰的原则完成种群交流 */
for (i = 0; i < EX_NUM; i++) {
j = rand()%TASK_NUM(pid);
if (population[j].fitness < ex_member[i].fitness) {
for (k = 0; k < NVARS; k++) {
population[j].gene[k] = ex_member[i].gene[k];
}
population[j].rfitness = 0;
population[j].cfitness = 0;
population[j].fitness = ex_member[i].fitness;
}
}
}
另外，全局最优解的归约由如下代码实现：
MPI_Op_create((MPI_User_function *)gene_max, 1, &my_op);
MPI_Reduce( local_best_individual, best_individual, NVARS+1,
MPI_DOUBLE, my_op, pnum-1, MPI_COMM_WORLD );
其中，具体的归约操作由如下函数实现：
void gene_max(double *in, double *inout, int *len, MPI_Datatype *dptr)
{
int i;

if (inout[0] < in[0]) { /* 比较适应度 */
for (i=0; i < *len; ++i) {
inout[i] = in[i]; /* 复制适应度较高的个体 */
}
}
}
五、 算法分析与实验结果
下面的实验结果是在192.169.129.47上利用结点vC0—168-129-48（slaver）和vC0—168-129-49（slaver）和vC0—168-129-46（slaver）获得的。用来计算最大值的函数为

其定义域如文件yangping.txt中所示，总种群大小为500，最大进化次数为2000。
进程个数
1
2
3
4
5
6
7
运算时间
32.308594
9.132812
4.335938
2.777344
3.699219
2.949219
2.621094
加速比
1
3.537639
7.451351
11.632910
8.733896
10.954966
12.326377
运
行
结
果
x1
22.455050
22.455050
22.455050
22.455050
22.455050
22.455050
22.455050
x2
7.205500
7.279000
7.289500
7.237000
7.289500
7.289500
7.289500
x3
19.269500
19.239000
19.269500
19.269500
19.269500
19.269500
19.269500
x4
4.944000
4.992000
4.072000
4.984000
4.888000
4.984000
4.968000
x5
-1.193000
-1.196500
-1.200000
-1.200000
-1.200000
-1.196500
-1.200000
x6
-9.120000
-9.113180
-9.072260
-9.120000
-9.120000
-9.120000
-9.120000
f
157521400.960884
157373629.694664
157325373.684378
157701606.886265
157673921.515628
157623544.425141
157702393.783168
进程个数
8
9
10
11
12
13
14
运算时间
2.226562
2.574219
2.449219
2.617188
2.289062
2.664062
2.597656
加速比
14.510530
12.550833
13.191386
12.344774
14.114338
12.127568
12.437595
运
行
结
果
x1
22.455050
22.455050
22.455050
22.455050
22.455050
22.455050
22.455050
x2
7.289500
7.289500
7.279000
7.289500
7.289500
7.289500
7.279000
x3
19.269500
19.269500
19.269500
19.269500
19.269500
19.269500
19.269500
x4
4.920000
4.976000
4.992000
4.984000
4.992000
4.968000
4.992000
x5
-1.200000
-1.200000
-1.200000
-1.200000
-1.200000
-1.196500
-1.200000
x6
-9.120000
-9.120000
-9.120000
-9.120000
-9.120000
-9.120000
-9.120000
f
157689427.608764
157704566.391334
157624693.663186
157706742.306205
157708921.527178
157619195.230127
157707891.211178
表1 实验结果
结果分析：
表1中最为有趣的现象是，当进程数小于5时，该算法的加速比似乎与进程数p存在一个平方关系，也就是说，存在一个超线性加速的关系。进程数大于等于5时，这种超线性加速实际也应该存在，只是由于节点数的限制，被进程管理的开销所限制。下面我们通过估计时间复杂性来分析造成这种超线性加速的原因。
如果将对染色体中每一变元上的一个计算看作一个基本计算，并设变元数为k，总种群中个体数为n，进程数则对每一进程，分析容易得到：pop_select函数最坏情形的时间复杂性为O((kn/p)2)，crossover函数最坏情形的时间复杂性为O(kn/p)，mutate函数最坏情形的时间复杂性为O(kn/p)，评估函数最坏情形的时间复杂性为O(kn/p)，elitist函数最坏情形的时间复杂性为O(n/p+k)。此外，按照算法的设计，在选择过程中的通讯所耗费的时间为O(kn/p)。综合可知，一次进化的时间复杂性为O((kn/p)2)。因此，所有进程总的计算时间最坏情形的渐近上界为O((kn)2/p)。而串行遗传算法一次进化的时间复杂性为O((kn)2)，这就解释了为什么p小于5的情形会具有超线性加速。当然，这并不能说明并行计算真能产生超线性加速比，因为我们可以非常有效地用一个进程来模拟p个进程的计算，也就是说在串行的情形下也能达到这样的加速。真正值得研究的问题是分析上述建立并行遗传算法收敛速度与串行遗传算法的收敛速度之间的关系。不过从表1可以看出，进程增加时，解得质量并没有任何降低。

六、 源程序清单
// genetic.cpp : 定义控制台应用程序的入口点。
//
//#include<stdafx.h>
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<time.h>

#define POPSIZE 500 /*种群大小*/
#define NVARS 6 /*染色体的长度*/
#define MAXGENS 2000
#define PXOVER 0.8 /* 杂交概率*/
#define PMUTATION 0.15 /* 变异概率*/

#define TASK_NUM(i) ((POPSIZE+pnum-1)/pnum)
//#define TASK_NUM(i) ((i)==pnum-1 && POPSIZE%((POPSIZE+pnum-1)/pnum)? POPSIZE%((POPSIZE+pnum-1)/pnum):(POPSIZE+pnum-1)/pnum)
#define EX_NUM 3

struct genotype
{
double gene[NVARS]; /*记录染色体的基因*/
double fitness; /* 适应度*/
double rfitness; /* 相对适应度*/
double cfitness; /* 累积适应度*/
} *population, *newpopulation;

double lowerbound[NVARS];
double upperbound[NVARS];
double local_best_individual[NVARS+1]; /* 局部最优解*/
double best_individual[NVARS+1]; /* 全局最优解*/

int generation;
double lower_fitness; /*最小的适应度*/
FILE *galog;

int pid, pnum; /*记录进程编号和进程数*/

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

void initialize(void) /*初始化函数*/
{
FILE *infile;
int i, j, r;
double lbound, ubound;
MPI_Status status;

population = (struct genotype*)malloc(sizeof(struct genotype)*(TASK_NUM(pid)+1)); //定义一个genotype类型指针，用于存放种群
newpopulation = (struct genotype*)malloc(sizeof(struct genotype)*(TASK_NUM(pid)+1)); //定义一个genotype类型指针newpopulation，用于存放产生新的种群

if (pid !=pnum - 1) {
if ((infile = fopen("yangping.txt","r"))==NULL) //打开文件，读入每个自变量的定义域
{
fprintf(galog,"\nCannot open input file!\n");
exit(1);
}

srand(time(0));
for (i=0; i<pnum; i++) {
r = rand();
MPI_Send(&r, 1, MPI_INT, i, 1, MPI_COMM_WORLD); //初始化每个进程的种群
}

for (i = 0; i < NVARS; i++) {
fscanf(infile, "%lf",&lowerbound[i]); //lowerbound记录各个自变量的定义域的下限
fscanf(infile, "%lf",&upperbound[i]); //upperbound记录各个自变量的定义域的上限
}
}
else {
MPI_Recv(&r, 1, MPI_INT, pnum-1, 1, MPI_COMM_WORLD, &status); //进程号为pnum-1的作为主进程，负责从其他子进程收集信息
srand(r);
}

MPI_Bcast(lowerbound, NVARS, MPI_DOUBLE, pnum-1, MPI_COMM_WORLD); //MPI_Bcast 是从一序号为pnum-1的主进程将一条消息广播发送到进程组内的所有进程
MPI_Bcast(upperbound, NVARS, MPI_DOUBLE, pnum-1, MPI_COMM_WORLD);

for (j = 0; j < TASK_NUM(pid); j++) {
population[j].fitness = 0; //初始化适应度为
population[j].rfitness = 0; //初始化相对适应度为
population[j].cfitness = 0; // 初始化累计适应度为

//printf("\nGene of member %d in %d is:\n", j, pid);
for (i = 0; i < NVARS; i++) {
population[j].gene[i] = randval(lowerbound[i], upperbound[i]); //将编码后的个体染色体放入种群数组population中
//printf(" %lf ", population[j].gene[i]);
}
}
if (pid == pnum - 1) fclose(infile);
}

/*进化函数*/
void evaluate(void)
{
int mem;
int i;
double x[NVARS+1];

lower_fitness = 0.0;

for (mem = 0; mem < TASK_NUM(pid); mem++) { //解码过程
for (i = 0; i < NVARS; i++)
x[i+1] = population[mem].gene[i]; //计算每个染色体中自变量的值

population[mem].fitness = (x[1]*x[1]*x[1]*x[1]*x[1]*x[1]) //定义每个染色体的适应度函数
- (x[1]*x[2]*x[3]*x[4]*x[4]*x[6]) + (x[3]*x[3]*x[3]*x[3]*x[3]*x[5]*x[6])
- (x[2]*x[4]*x[5]*x[5]*x[6]) - (x[1]*x[1]*x[3]*x[3]*x[5]*x[5])
+ (x[2]*x[2]*x[3]*x[5]*x[6]) + (x[1]*x[2]*x[4]*x[4]*x[4]*x[5]);
//printf(" The fitness of %d at %lf, %lf, %lf, %lf, %lf, %lf is %lf\n", mem, x[0], x[1], x[2], x[3], x[4], x[5], population[mem].fitness);
if (population[mem].fitness < lower_fitness) lower_fitness = population[mem].fitness; //记录最小的适应度
}
}

/*保存最优个体函数*/
void keep_the_best()
{
int mem, i, cur_best = 0;

for (mem = 0; mem < TASK_NUM(pid); mem++) {
if (population[mem].fitness > population[TASK_NUM(pid)].fitness) {
cur_best = mem; //用Cur_best记录最大适应度的染色体
population[TASK_NUM(pid)].fitness = population[mem].fitness;
}
}
//printf("\nG%4d: the best individual in Process %d is %d, the fitness is %lf.\n", generation, pid, cur_best, population[cur_best].fitness);
for (i = 0; i < NVARS; i++) {
population[TASK_NUM(pid)].gene[i] = population[cur_best].gene[i];
local_best_individual[i+1] = population[cur_best].gene[i]; //保存最有个体
}
local_best_individual[0] = population[cur_best].fitness;
}

/*选择函数*/
void pop_select(void)
{
MPI_Status status; //通信状态
MPI_Request handle; //定义一个通信请求对象
int mem, i, j, k;
double sum = 0;
double p;
static struct genotype ex_member[EX_NUM]; //定义一个genotype的数组用以记录交流个体

/*计算子种群的总适应度*/
for (mem = 0; mem < TASK_NUM(pid); mem++) {
sum += (population[mem].fitness - lower_fitness);
//printf(" Sum %d is %lf\n", mem, sum);
}

/* 计算各个体相应适应度*/
for (mem = 0; mem < TASK_NUM(pid); mem++) {
population[mem].rfitness = (population[mem].fitness - lower_fitness)/sum;
//printf("The fitness of %d is %lf, the rfitness is %lf\n", mem, population[mem].fitness, population[mem].rfitness);
}
population[0].cfitness = population[0].rfitness;

/* 计算各个体累积适应度*/
for (mem = 1; mem < TASK_NUM(pid); mem++) {
population[mem].cfitness = population[mem-1].cfitness + population[mem].rfitness;
}
/* 按照累积适应度概率选取种群交流个体，并发送和接收*/
for (i = 1; i <= EX_NUM; i++) {
p = rand()%1000/1000.0;
if (p < population[0].cfitness) {
MPI_Isend(&population[0], sizeof(struct genotype)/sizeof(char),
MP