前言
通过代码解释遗传算法解决旅行商问题的思路(CUDA版本，使用GPU并行加速)。
首先，定义常量参数
const int nCities = 99; //城市数量
const int GEN_OUT = 15; //最大代数
const int GEN_IN = 1; //最大代数
const int UNIT_NUM = 20000; //群体规模为20000
const int BlockDIm = 5;//开五个Block
const int ThreadDIm = 200;//每个Blockl里开两百个线程
const int EachBlock = 20;//每个线程处理20个个体
const int Exchage = 10;
const float PM = 0.1; //变异概率
const float PS = 0.8;//保留概率
//以PS 的概率选择前 num2 个解，抛弃其后的num1 个解。
const int NUM1 = EachBlock * ( 1 - PS);
const int NUM2 = EachBlock * PS;
//城市间的距离
__constant__ float data[nCities * nCities]
同时，定义两个结构体，分别存放城市信息和个体信息
//城市
struct node
{
int num;//城市的编号
float x;//横坐标
float y;//纵坐标
}nodes[nCities]; //定义99维度的数组，用来存放99座城市的信息
//个体
struct unit
{
float length;//代价，总长度
int path[nCities];//路径 （定义99维度的数组，用来存放每一种不重复遍历全部99座城市的方法）
};
然后从主函数开始看起
int main(int argc, char* argv[])
{
struct timeval start;//记录开始时间
struct timeval end;//记录结束时间
unsigned long timer;//用于记录总用时的变量（timeval end - timeval start）
struct unit bestone;//bestone 用于保存最好的一个解
//群体规模（群体规模是指有不同的解，而bestone 用于保存最好的一个解)
struct unit group[UNIT_NUM]; //建立两万个个体
srand(time(NULL));
init_dis(); //读入99座城市的信息，包括城市的编号，横纵坐标，存放在node结构体中
gettimeofday(&start,NULL);
//个体初始化,进化,选择,交叉,变异，交换
bestone = GAStart( group );
gettimeofday(&end,NULL);
timer = 1000000 * (end.tv_sec-start.tv_sec)+ end.tv_usec-start.tv_usec;
printf("timer = %ld us\n",timer);
cout << "变异概率PM = " << PM << endl;
cout << "保留概率PS = " << PS << endl;
cout << "各子群交换次数 ＝ " << GEN_OUT << endl;
cout << "子群迭代次数 ＝ " << GEN_IN << endl;
cout << "群体规模 ＝ " << UNIT_NUM << endl;
print(bestone);
}
init_dis()函数的主要作用是读入99座城市的信息，包括城市的编号，横纵坐标，存放在node结构体中，具体函数如下：
void init_dis()
{
int i, j;
ifstream in("source.txt");
if(in.is_open())//打开source.txt文件
{
for (i = 0; i < nCities; i++)
{
in >> nodes[i].num >> nodes[i].x >> nodes[i].y; //读入城市序号以及相应的横纵坐标
}
}
else
{
cout<<"open source.txt error"<<endl;//如何没有正确读入则报错
}
float *temp_s = (float *)malloc( sizeof(float) * nCities * nCities );//开一片99*99的内存空间，用来存放任意两座城市之间的距离，并且使用temp_s指针指向这片内存
//下面给城市距离赋值，相当于定义了一个99*99的矩阵，这是为了解释方便，记该矩阵为矩阵A，第i行对应第i个城市到每一个城市的距离
for(int i=0;i<nCities;i++)
{
temp_s[i * nCities + i] = (float)INT_MAX; //自己到自己的距离，随意定义距离值，没有什么具体意义
for (j = i + 1; j < nCities; j++)
{
temp_s[i * nCities + j] = temp_s[j * nCities + i] = sqrt(
(nodes[i].x - nodes[j].x) * (nodes[i].x - nodes[j].x) +
(nodes[i].y - nodes[j].y) * (nodes[i].y - nodes[j].y) );
}//根据每座城市的x，y坐标，运用两点之间的距离公式求解任意两座城市之间的距离，显然，矩阵A为对称矩阵（A(i,j)和A(j,i)都表示i,j两座城市之间的距离）
}
cudaMemcpyToSymbol( data, temp_s, sizeof(float) * nCities * nCities);//放入到专门存在任意城市之间距离的数组变量data中
free( temp_s ); //释放指针
}
回到main函数中，函数gettimeofday(&start,NULL)表示开始计时，然后运行bestone = GAStart( group )，从而跳转到GAStart中，GAStart()是这里的核心函数，代码如下：
struct unit GAStart( unit group[])
{
unit temp;//定义的用来存放个体的结构体，一个个体即不重复的包含1-99座城市序列及其总长度的结构体
//1, 2, 3, 4 ...... nCities
generate(temp.path, temp.path + nCities, GenByOne(0)); //temp里的path数组（99维）依次赋值1到99
//初始化UNIT_NUM个解，也就是群体
for (int i = 0; i < UNIT_NUM; i++)
{
random_shuffle(temp.path, temp.path + nCities); //对path里的1-99数字进行随机打乱，构造个体
memcpy(&group[i], &temp, sizeof(temp)); //将构造的个体存入到group中，一共构造了20000个个体
}
//保存最好的一个解
struct unit bestone[1];
struct unit *d_in;
struct unit *b_tmp;
struct unit *b_arr;
cudaMalloc((void**)&d_in, UNIT_NUM * sizeof(struct unit));//开放了两万个个体存贮所需显存
cudaMalloc((void**)(&b_tmp), sizeof(struct unit));//存放最优个体
cudaMalloc((void**)(&b_arr), ThreadDIm * BlockDIm * sizeof(struct unit));// Block数目乘上每个Block里所开的线程数，等于总的Block数目
cudaMemcpy(d_in, group, UNIT_NUM * sizeof(struct unit), cudaMemcpyHostToDevice);//将group中的数据（20000个个体信息）复制到d_in中
//这里需要注意group在计算机的内存中，而d_in在显卡的显存中，通过这一步将内存中的数据传入到显存中，下面可以用cuda 进行GPU并行计算
InitLength <<<BlockDIm, ThreadDIm >>>(d_in);
for (int i1 = 0 ; i1 < GEN_OUT ; i1++)
{ cout<<"round: "<<i1<<endl;
GA <<<BlockDIm, ThreadDIm >>>(d_in);
Exchange <<<BlockDIm, ThreadDIm >>>(d_in);
}
FinalSort <<<BlockDIm, ThreadDIm >>>(d_in,b_arr);
BestOne <<<1, 1 >>>(b_tmp ,b_arr);
//cudaMemcpy(group, d_in, UNIT_NUM * sizeof(struct unit), cudaMemcpyDeviceToHost);
cudaMemcpy(bestone, b_tmp, sizeof(struct unit), cudaMemcpyDeviceToHost);
cudaFree(d_in);
cudaFree(b_arr);
cudaFree(b_tmp);
//for (int i = 0; i < UNIT_NUM; i++)
//{
// print(group[i]);
//}
return bestone[0];
}
执行InitLength <<<BlockDIm, ThreadDIm >>>(d_in)，接着来看InitLength()函数
__global__ void InitLength(unit *dataIn)
{
int thread = blockIdx.x * ThreadDIm + threadIdx.x;//当前线程号=当前block序号*每个block中的线程数+当前block中的线程号序号
//一个线程处理EachBlock组距离
int index = thread * EachBlock;//EachBlock=20，即每个线程处理二十个个体（5个block，每个block中有200个线程，每个线程20个个体，共20000个个体）
for (int i = 0;i < EachBlock;i++)
{
CalCulate_length(dataIn[index + i]); //计算每个线程中对应的二十个个体里的城市序列的总距离，dataIn总维度为20000对应每个个体对应的总距离
}
}
CalCulate_length()函数用来遍历每个个体中城市序列，求得总距离，data中包含了每两个城市间的距离如下：
__device__ void CalCulate_length(unit &p)
{
p.length = 0;
for (int j = 1; j < nCities; j++)
{
p.length += data[ (p.path[j-1]) * nCities + p.path[j] ];
}
p.length += data[ (p.path[nCities - 1]) * nCities + p.path[0] ];
}
回到GAStart()中，开始下列迭代，一共迭代十五次(GEN_OUT是迭代总次数)
for (int i1 = 0 ; i1 < GEN_OUT ; i1++)
{ cout<<"round: "<<i1<<endl;
GA <<<BlockDIm, ThreadDIm >>>(d_in);
Exchange <<<BlockDIm, ThreadDIm >>>(d_in);
}
先进入到GA()函数中，函数代码如下：
__global__ void GA(unit *dataIn)
{
int thread = blockIdx.x * ThreadDIm + threadIdx.x;//和上述同理，计算当前线程数
int index = thread * EachBlock;//计算index
for (int i2 = 0 ; i2 < GEN_IN ; i2++)//GEN_IN为1，可以省略
{
quicksort(dataIn + index , EachBlock);//将每个线程中对应的二十个个体对应的总距离进行升序排序(总目标为求得总距离最短的对应个体)
for (int j = 0; j <= NUM1 - 1; j++) //NUM1 = EachBlock * ( 1 - PS)，PS = 0.8，故NUM=4
{
//复制操作
memcpy(&dataIn[index + NUM2 + j], &dataIn[index + j], sizeof(unit));//NUM2为16，即用距离最短的4个个体（前4个个体）替换距离最长的4个个体（后4个个体）
}
for (int jq = 0; jq < EachBlock / 2; jq++)
{
Cross_group(dataIn[index + jq], dataIn[index + EachBlock - jq -1]); //这里为交叉操作每个线程对应的二十个个体，第i个个体和第19-i个个体做交叉（数组元素索引从0开始）
}
//变异
Varation_group(&dataIn[index]);
}
quicksort(dataIn + index , EachBlock);
// printf("%f\t",(*(dataIn+index)).length);
}
复制操作用距离最短的4个个体（前4个个体）替换距离最长的4个个体（后4个个体），下面进入到交叉函数中：
__device__ void Cross_group( unit &p, unit &q)
{
unsigned int seed = blockIdx.x * ThreadDIm + threadIdx.x;
__shared__ curandState s;
// seed a random number generator
curand_init(seed, 0, 0, &s);
int i = 0;
int pos1, pos2; //分别用来记录p,q个体中的具体第几位的索引
int len = nCities; //len值为99
int first;
float len1 = data[p.path[0] * nCities + p.path[1] ];//计算p个体里第一座城市到第二座城市距离
float len2 = data[q.path[0] * nCities + q.path[1] ];//计算q个体里第一座城市到第二座城市距离
if (len1 <= len2)
{
first = p.path[0]; //如果len1小，记录p个体中第一个城市的序号
}
else
{
first = q.path[0]; //如果len2小，记录q个体中第一个城市的序号
}
pos1 = Search_son( p.path + i, len, first); //记录first中对应的城市序号在p个体中的索引，即p中的第几个城市为first所对应的城市（本程序中城市名序号即为城市名）
pos2 = Search_son( q.path + i, len, first); //记录first中对应的城市序号在q个体中的索引，即q中的第几个城市为first所对应的城市
Rotate(p.path + i, len, pos1); //翻转
Rotate(q.path + i, len, pos2);
while ( --len > 1)
{
i++;
float span1 = data[ p.path[i - 1] * nCities + p.path[i] ];
float span2 = data[ q.path[i - 1] * nCities + q.path[i] ];
if ( span1 <= span2 )
{
pos2 = Search_son( q.path + i, len, p.path[i]);
Rotate(q.path + i, len, pos2);
}
else
{
pos1 = Search_son( p.path + i, len, q.path[i]);
Rotate(p.path + i, len, pos1);
}
}
//Rotate(q.path, nCities, 44);
Rotate(q.path, nCities, curand(&s) % nCities);
CalCulate_length(p);
CalCulate_length(q);
}
这里，我们就Rotate()函数做进一步说明，Rotate()主要调用了三次Reverse()函数，这里为了解释方便，将两个函数放在一起：
//旋转 m 位
__device__ void Rotate(int path[],int len, int m)
{
if( m < 0 )
{
return;
}
if (m > len)
{
m %= len;
}
Reverse(path, 0, m -1); //将path中第m个对应城市前面的城市顺序进行倒置
Reverse(path, m, len -1); //将path中第m个对应城市及其后面的城市顺序进行倒置
Reverse(path, 0, len -1); //将path中所有城市顺序进行倒置
}
__device__ void Reverse(int path[], int b, int e) //将path中从path[b]到path[e]做倒置
{
int temp;
while (b < e)
{
temp = path[b];
path[b] = path[e];
path[e] = temp;
b++;
e--;
}
}
经过Rotate()函数处理后，实际相当于将第m个城市到最后一个城市移到到了原始的第一个城市到第m-1个城市的前面，再来看整个交叉函数中的交叉处理
float len1 = data[p.path[0] * nCities + p.path[1] ];
float len2 = data[q.path[0] * nCities + q.path[1] ];
if (len1 <= len2)
{
first = p.path[0];
}
else
{
first = q.path[0];
}
pos1 = Search_son( p.path + i, len, first);
pos2 = Search_son( q.path + i, len, first);
Rotate(p.path + i, len, pos1);
Rotate(q.path + i, len, pos2);
while ( --len > 1)
{
i++;
float span1 = data[ p.path[i - 1] * nCities + p.path[i] ];
float span2 = data[ q.path[i - 1] * nCities + q.path[i] ];
if ( span1 <= span2 )
{
pos2 = Search_son( q.path + i, len, p.path[i]);
Rotate(q.path + i, len, pos2);
}
else
{
pos1 = Search_son( p.path + i, len, q.path[i]);
Rotate(p.path + i, len, pos1);
}
}
实际相当于p,q个体对应的第i座城市和第i+1座城市之间的距离依次对比，假设p对应的距离短，则对q个体的第i座城市到最后一座城市的顺序以p中第i个城市在q个体中的位置为依据进行旋转处理，最终p和q个体中应该包含一致的城市序列，即两个个体经过层层比较得出的最优排列，为了保证种群的多样性，我们q个体进行一定的处理，用随机数对q个体进行旋转
Rotate(q.path, nCities, curand(&s) % nCities);
然后分别计算p,q个体的总长度，至此，交叉操作完成。
接着我们来看变异操作：
__device__ void Varation_group(unit group[])
{
int i, j, k;
int temp;
//变异的数量，即，群体中的个体以PM的概率变异,变异概率不宜太大
int num = EachBlock * PM;
unsigned int seed = blockIdx.x * ThreadDIm + threadIdx.x;
__shared__ curandState s;
// seed a random number generator
curand_init(seed, 0, 0, &s);
while (num--)
{
//确定发生变异的个体
k = curand(&s) % EachBlock;
//确定发生变异的位
i = curand(&s) % nCities;
j = curand(&s) % nCities;
temp = group[k].path[i];
group[k].path[i] = group[k].path[j];
group[k].path[j] = temp;
CalCulate_length(group[k]);
}
}
变异的操作也是通过生成随机数的方式，随机生成变异的个体号，对应变异个体中需要交换顺序的两个城市，然后进行交换处理即可。变异操作完成后对每个线程对应的二十个个体的总距离在进行一次升序的快速排序，至此GA()函数执行完毕，然后回到GAStart()函数中，接一下进入GAStart()函数中的Exchange()函数中：
__global__ void Exchange(unit *dataIn)
{
int thread = blockIdx.x * ThreadDIm + threadIdx.x;
int index = thread * EachBlock;
for (int ee = 0 ; ee < Exchage ; ee++)
{
if (dataIn[index].length < dataIn[((thread + 1 + ThreadDIm * BlockDIm/Exchage * ee)%(ThreadDIm * BlockDIm)) * EachBlock + EachBlock - ee - 1].length)
{
memcpy(&dataIn[((thread + 1 + ThreadDIm * BlockDIm/Exchage * ee)%(ThreadDIm * BlockDIm)) * EachBlock + EachBlock - ee - 1], &dataIn[index], sizeof(unit));
}
}
}
Exchange()函相当于适应度函数（相当于自然选择，选择更好的个体类型），对每个线程的第一个个体（最优距离个体）和附近的Exchage（此处为10）个线程中最优个体对比，如果距离更小，则替代附近线程中的最优个体。
经历十五次迭代之后，在经过FinalSort()和BestOne()函数，得到最优个体。
2223
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
