遗传算法求解旅行商问题

1.遗传算法

        遗传算法是受大自然的启发，模拟生物在自然环境中的遗传和进化过程而形成的一种自适应、具有全局优化能力的随机搜索算法。

自然界的进化包括3个原则：

（1）适者生存原则，这意味着适应能力强的物种，会在残酷的竞争中生存下来，而适应能力差的物种会逐渐地消亡。

(2) 两性繁殖。这意味着种群中性别不同的个体，生活在一起，产生新的个体。

(3) 变异。 由于环境的变化，新物种的出现，以及不同物种的交互都会引起种群的变异。

        遗传算法的思路是通过从给定一个初始群体出发，利用选择算子、杂交算子以及变异算子来模拟自然进化的三种原则，逐步改进种群，一步步逼近最优解，以达到求解最优华问题的目的。

GA算法的计算步骤：

        记住遗传算法的过程很重要，首先是初始化一群解，然后再在这些解中选择较优的一部分，将选择的这一部分解进行交叉，且以一定概率变异，（交叉一般能得到比当前解更好的解，而变异则很可能让结果变差，所以变异的概率一般不是很大，但是这样有助于我们跳出局部最优）。交叉变异以后进行群体更新，对于TSP问题，群体更新时保存这一次迭代产生的最好的解，然后继续进行下一次迭代，直到满足终结条件为止。

GA的算法过程：

        初始化t，t代表while循环已经迭代了多少次。其中f(pop(t))是指这个解的适应度，对于TSP问题，适应度就是它的代价，第8行是按一定的概率选择较有的解。第9行以Pc概率进行交叉，第10行以Pm概率进行变异。

2. 问题建模

        遗传算法其实很简单，就是初始化一群解，然后选择这一群里面较优的解，在较优的解里面，让其中的个体交叉，使得交叉后得到更好的解，再按一定概率进行变异，希望变异能跳出局部最优。对于遗传算法求解TSP问题，最难的地方在于问题建模，刚开始根本不知道如何用遗传算法来求解旅行商问题，如何交叉，如何变异。

        首先初始化一群解，可以通过C++提供的库函数来产生一个城市的随机排列，每一个排列代表一个解random_shuffle(temp.path, temp.path + nCities)。然后以一定概率选择较优的解，选择的方法有很多，我们不一定非要按照上面伪代码的方式来选择，比如我们希望每次保存当前这群解中的前60%,则我们可以按解的适应度排序，然后取前60%的解，对于后40%的解，我们可以用前40%的解去覆盖它，则前40%的解就有2个副本，只要我们交叉的时候不要让相同的两个副本交叉就行了，因为相同的两个解交叉，不会让结果变得更好。

        变异也很简单，只需要在一个解中随机的选择两个城市，然后交换它们即可。注意，变异的概率不宜太大。

        最难的部分是交叉，我们要如何用两个解得到一个更好的解？这就是交叉，让一代比一代强，我们才可能慢慢接近最优解。交叉的方法有很多，可以参考http://blog.csdn.net/xuyuanfan77/article/details/6726477

源码中采用类似于三交换启发交叉(THGA)，我把它改成了二交叉的。

三交换启发交叉方法的基本思想如下:

选3个参加交配的染色体作为父代,以8个城市为例来说明这一过程,其中dij由前面的表1给出,父代染色体为

A　＝　3　2　1　4　8　7　6　5

B　＝　2　4　6　8　1　3　5　7

C　＝　8　7　5　6　4　3　2　1

SUM1＝42,SUM2＝40,SUM3＝46(SUM1，SUM2，SUM3分别为这3种排法所走的距离总和数).

随机选出初始城市j=1,Sj=3右转动，使3成为3父代的第1位置.

A　＝　3　2　1　4　8　7　6　5

B　＝　3　5　7　2　4　6　8　1

C　＝　3　2　1　8　7　5　6　4

由于d(3,2)＞d(3,5),所以有:

A　＝　×　5　2　1　4　8　7　6

B　＝　×　5　7　2　4　6　8　1

C　＝　×　5　6　4　2　1　8　7

由此规则计算可得:

O　＝　3　5　7　6　8　4　2　1

我们本来是3个不同的解，现在得到了一个比三个解都优的解，总不能让原来的三个解都等于现在的这个局部最优解吧，这样不利于下次交叉，我们可以用如下的方法改变另外两个解的路径：

Rotate(q.path, nCities, rand() % nCities);

        

        上行代码执行以后，它的代价还是和原来一样的，路径也是一样，只是起点变了，这样有什么好处呢？有利于下次交叉的时候，原来的两个相同代价，不同路径的解能和其他解交叉出不同的结果，这样有利于找到更好的解。

3. 代码实现

/* * * * * * * Copyright(c) Computer Science Department of XiaMen University  * * * * * * Authored by lalor on: 2012年 06月 29日 星期五 23:49:57 CST * * * * * * Email: mingxinglai(at)gmail.com * * * * * * @desc: * * * * * * @history * * * * * * 说明：本程序使用的测试数据来自权威的benchmark,其最优解是1211.数据保存在source.txt * * 本例的测试数据来自http://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/tsp/ * * rat99.tsp.gz * * 数据如下 * ﹡格式：（城市编号，横坐标，纵坐标） 1  6  4  2 15 15  3 24 18  4 33 12  5 48 12  6 57 14  7 67 10  8 77 10  9 86 15 10  6 21 11 17 26 12 23 25 13 32 35 14 43 23 15 55 35 16 65 36 17 78 39 18 87 35 19  3 53 20 12 44 21 28 53 22 33 49 23 47 46 24 55 52 25 64 50 26 71 57 27 87 57 28  4 72 29 15 78 30 22 70 31 34 71 32 42 79 33 54 77 34 66 79 35 78 67 36 87 73 37  7 81 38 17 95 39 26 98 40 32 97 41 43 88 42 57 89 43 64 85 44 78 83 45 83 98 46  5 109 47 13 111 48 25 102 49 38 119 50 46 107 51 58 110 52 67 110 53 74 113 54 88 110 55  2 124 56 17 134 57 23 129 58 36 131 59 42 137 60 53 123 61 63 135 62 72 134 63 87 129 64  2 146 65 16 147 66 25 153 67 38 155 68 42 158 69 57 154 70 66 151 71 73 151 72 86 149 73  5 177 74 13 162 75 25 169 76 35 177 77 46 172 78 54 166 79 65 174 80 73 161 81 86 162 82  2 195 83 14 196 84 28 189 85 38 187 86 46 195 87 57 194 88 63 188 89 77 193 90 85 194 91  8 211 92 12 217 93 22 210 94 34 216 95 47 203 96 58 213 97 66 206 98 78 210 99 85 204 ﹡* * * * */#include <iostream>#include <string.h>#include <fstream>#include <iterator>#include <algorithm>#include <limits.h>#include <math.h>#include <stdlib.h>using namespace std;const int nCities = 99; //No. of node//const double PC = 0.9; //交叉概率double PM = 0.1; //变异概率double PS = 0.8;//保留概率int GEN_MAX = 50; //最大代数const int UNIT_NUM = 5000; //群体规模为50double length_table[nCities][nCities];//distance//城市struct node{int num;//城市的编号int x;//横坐标int y;//纵坐标}nodes[nCities];struct unit{double length;//代价，总长度int path[nCities];//路径bool operator < ( const struct unit &other) const //用于群体的排序{return length < other.length;}};//群体规模（群体规模是指有 UNIT_NUM 个不同的解，而bestone 用于保存最好的一个解)struct unit group[UNIT_NUM];//保存最好的一个解unit bestone = {INT_MAX, {0} };// create matrix to storage the Distance each cityvoid init_dis(); //计算 unit 中的length， 也就是群体的一个个体（一个解）的长度void CalCulate_length(unit &p);//查找id （代表城市） 在当前解中的位置，用于两个解的交叉int search_son(unit &p, int id);//打印一个解void print( unit &p);//初始化群体，由C++ 中的 random_shuff 产生一个随机排列void Initial_group( unit group[]);//开始进化，在本函数中执行群体中个体的交叉和变异void Evolution_group(unit group[]);//变异，随机的选择一个群体，然后随机选择两个点，交换它们的位置void Varation_group(unit group[]);//交叉void Cross_group( unit &p, unit &q);int main(int argc, char* argv[]){srand(time(NULL));init_dis();//初始化种群Initial_group( group );//种群进化：选择，交叉，变异Evolution_group( group );cout << "变异概率PM = " << PM << endl;cout << "保留概率PS = " << PS << endl;cout << "最大代数 ＝ " << GEN_MAX << endl;cout << "群体规模 ＝ " << UNIT_NUM  << endl;cout << "代价是： = " << bestone.length << endl;print(bestone);}void init_dis() // create matrix to storage the Distance each city{int i, j;ifstream in("source.txt");for (i = 0; i < nCities; i++){in >> nodes[i].num >> nodes[i].x >> nodes[i].y;}for (i = 0; i < nCities; i++){length_table[i][i] = (double)INT_MAX;for (j = i + 1; j < nCities; j++){length_table [i][j] = length_table[j][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) );}}}void CalCulate_length(unit &p){int j = 0;p.length = 0;for (j = 1; j < nCities; j++) {p.length += length_table[ p.path[j-1] ][ p.path[j] ];}p.length += length_table[ p.path[nCities - 1] ][ p.path[0] ];}void print( unit &p){int i;cout << "代价是：" << p.length << endl << "路径是：";//for (i = 0; i < nCities; i++) //{//cout << p.path[i] << " ";//}copy(p.path, p.path + nCities, ostream_iterator<int>(cout, " -> "));cout << p.path[0] << endl;}//函数对象，给generate 调用class GenByOne{public:GenByOne (int _seed = -1): seed(_seed){}int operator() (){return seed += 1;}private:int seed;};//随机产生 UNIT_NUM 个解空间void Initial_group( unit group[]){int i, j;unit temp;//1, 2, 3, 4 ...... nCitiesgenerate(temp.path, temp.path + nCities, GenByOne(0));// 产生 UNIT_NUM 个解，也就是群体for (i = 0; i < UNIT_NUM; i++) {//产生一个随机排列,也就是初始化一个解random_shuffle(temp.path, temp.path + nCities);memcpy(&group[i], &temp, sizeof(temp));CalCulate_length(group[i]);}}void Evolution_group(unit group[]){int i, j;int n = GEN_MAX;int num1, num2;//以PS 的概率选择前 num2 个解，抛弃其后的num1 个解。num1 = UNIT_NUM * ( 1 - PS);num2 = UNIT_NUM * PS;//迭代几次，即繁衍多少代while (n-- ) //循环GEN-MAX次{//选择部分优秀的种群sort(group, group + UNIT_NUM);if (group[0].length < bestone.length) {memcpy(&bestone, &group[0], sizeof(unit));}for (j = 0; j <=  num1 - 1; j++) {memcpy(&group[ num2 + j], &group[j], sizeof(unit));}//交叉for (j = 0; j < UNIT_NUM / 2; j+= 1) {Cross_group(group[j], group[ UNIT_NUM - j -1]);}//变异Varation_group(group);}//保存已找最好的解sort(group, group + UNIT_NUM);if (group[0].length < bestone.length) {memcpy(&bestone, &group[0], sizeof(unit));}}void Varation_group(unit group[]){int i, j, k;double temp;//变异的数量，即，群体中的个体以PM的概率变异,变异概率不宜太大int num = UNIT_NUM * PM;while (num--) {//确定发生变异的个体k = rand() % UNIT_NUM;//确定发生变异的位i = rand() % nCities;j = rand() % nCities;//exchangetemp  = group[k].path[i];group[k].path[i] = group[k].path[j]; group[k].path[j] = temp;CalCulate_length(group[k]);}}int Search_son( int path[], int len, int city){if (city <= 0 || city > nCities) {cout << "city outfiled, city = " << city << endl;return -1;}int i = 0;for (i = 0; i < len; i++) {if (path[i] == city) {return i;}}return -1;}//reverse a array//it's a auxiliary function for Rotate() void Reverse(int path[], int b, int e){int temp;while (b < e) {temp = path[b];path[b] = path[e];path[e] = temp;b++;e--;}}//旋转 m 位void Rotate(int path[],int len, int m){if( m < 0 ){return;}if (m > len) {m %= len;}Reverse(path, 0, m -1);Reverse(path, m, len -1);Reverse(path, 0, len -1);}void Cross_group( unit &p, unit &q){int i = 0, j ,k;int pos1, pos2;int len = nCities;int first;double len1 = length_table[p.path[0] ][ p.path[1] ];double len2 = length_table[q.path[0] ][ 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++;double span1  = length_table[ p.path[i - 1] ][ p.path[i] ];double span2  = length_table[ q.path[i - 1] ][ 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, rand() % nCities);CalCulate_length(p);CalCulate_length(q);}

4.参考资料：

[1] http://blog.csdn.net/xuyuanfan77/article/details/6726477

[2] 算法设计与分析（高级教程），国防工业出版社