TSP问题
TSP（Traveling Salesman Problem）是典型的NP完全问题，即其最坏情况下的时间复杂度随着问题规模的增大按指数方式增长。TSP问题可以描述为：已知n个城市之间的相互距离，某一旅行商从某一个城市出发，访问每个城市一次且仅一次，最后回到出发的城市，如何安排才能使其所走的路线最短。
换成数学模型可表述为：搜索自然子集 X = { 1 , 2 , . . . , n } X=\{1,2,...,n\} X={1,2,...,n}（n表示对城市的编号）的一个排列 P ( X ) = { V 1 , V 2 , . . . . , V n } P(X)=\{V1,V2,....,Vn\} P(X)={V1,V2,....,Vn}，使得 D = ∑ d ( V i , V i + 1 ) + d ( V n , V 1 ) D=∑d(V_i,V_{i+1})+d(Vn,V1) D=∑d(Vi,Vi+1)+d(Vn,V1)取最小值，其中 d ( V i , V i + 1 ) d(V_i,V_{i+1}) d(Vi,Vi+1) 表示城市 V i V_i Vi 到 V i + 1 V_{i+1} Vi+1 的距离。
遗传算法步骤
编码问题：由于这是一个离散型的问题，我们采用整数编码的方式，用 1 1 1 ~ n n n 来表示 n n n 个城市， 1 1 1 ~ n n n 的任意一个排列就构成了问题的一个解。可以知道，对于 n n n 个城市的TSP问题，一共有 n ! n! n! 种不同的路线。
种群初始化：对于 N N N 个个体的种群，随机给出 N N N 个问题的解（相当于是染色体）作为初始种群。这里具体采用的方法是： 1 , 2 ， . . . , n 1,2，...,n 1,2，...,n 作为第一个个体，然后 2 , 3 ， . . n 2,3，..n 2,3，..n 分别与 1 1 1 交换位置得到 n − 1 n-1 n−1 个解，从 2 2 2 开始， 3 , 4 , . . . , n 3,4,...,n 3,4,...,n 分别与 2 2 2 交换位置得到 n − 2 n-2 n−2 个解，依次类推。（如果这样还不够初始种群的数量，可以再考虑 n , n − 1 , . . . , 1 n,n-1,...,1 n,n−1,...,1 这个序列，然后再按照相同的方法生成）
适应度函数：设一个解遍历初始行走的总距离为D，则适应度 f i t n e s s = 1 D fitness=\frac{1}{D} fitness=D1 。即总距离越高，适应度越低，总距离越低（解越好），适应度越高，这里的适应值可以取总距离的相反数。
选择操作：个体被选中的概率与适应度成正比，适应度越高，个体被选中的概率越大。这里仍然采用轮盘赌法。
交叉操作：交叉操作是遗传算法最重要的操作，是产生新个体的主要来源，直接关系到算法的全局寻优能力，这里采用部分映射交叉(Partial Mapping Crossover，PMX) 。比如对于 n = 10 n=10 n=10 的情况，对于两个路径:
$$
1\rightarrow\textcolor[rgb]{1,0,0}{2}\rightarrow4\rightarrow\textcolor[rgb]{1,0,0}5\rightarrow6\rightarrow3\rightarrow9\rightarrow10\rightarrow8\rightarrow7\3\rightarrow\textcolor[rgb]{1,0,0}{9}\rightarrow7\rightarrow\textcolor[rgb]{1,0,0}6\rightarrow8\rightarrow10\rightarrow5\rightarrow1\rightarrow2\rightarrow4
KaTeX parse error: Can't use function '$' in math mode at position 12: 随机产生两个 $̲[1,10]$ 之间的随机数 …
1\rightarrow\textcolor[rgb]{1,0,1}{9\rightarrow7\rightarrow6}\rightarrow\textcolor[rgb]{0,1,0}6\rightarrow3\rightarrow\textcolor[rgb]{0,1,0}9\rightarrow10\rightarrow8\rightarrow\textcolor[rgb]{0,1,0}7\3\rightarrow\textcolor[rgb]{1,0,1}{2\rightarrow4\rightarrow5}\rightarrow8\rightarrow10\rightarrow\textcolor[rgb]{0,1,0}5\rightarrow1\rightarrow\textcolor[rgb]{0,1,0}2\rightarrow\textcolor[rgb]{0,1,0}4
紫 色 部 分 表 示 交 叉 过 来 的 基 因 ， 这 个 时 候 会 发 现 可 能 交 叉 过 来 的 基 因 与 原 来 其 他 位 置 上 的 基 因 有 重 复 ， 容 易 直 到 ， 第 一 个 个 体 重 复 基 因 的 数 目 与 第 二 个 个 体 重 复 基 因 的 数 目 是 相 同 的 （ 这 里 都 是 3 个 ） ， 只 需 要 把 第 一 个 个 体 重 复 的 基 因 ( 用 绿 色 标 识 ) 与 第 二 个 个 体 重 复 的 基 因 做 交 换 ， 即 可 以 消 除 冲 突 。 消 除 冲 突 之 后 的 解 如 下 : 紫色部分表示交叉过来的基因，这个时候会发现可能交叉过来的基因与原来其他位置上的基因有重复，容易直到，第一个个体重复基因的数目与第二个个体重复基因的数目是相同的（这里都是3个），只需要把第一个个体重复的基因(用绿色标识)与第二个个体重复的基因做交换，即可以消除冲突。消除冲突之后的解如下: 紫色部分表示交叉过来的基因，这个时候会发现可能交叉过来的基因与原来其他位置上的基因有重复，容易直到，第一个个体重复基因的数目与第二个个体重复基因的数目是相同的（这里都是3个），只需要把第一个个体重复的基因(用绿色标识)与第二个个体重复的基因做交换，即可以消除冲突。消除冲突之后的解如下:
1 \rightarrow9 \rightarrow7 \rightarrow6 \rightarrow5 \rightarrow3 \rightarrow2 \rightarrow10\rightarrow 8 \rightarrow4
\
3\rightarrow2 \rightarrow4 \rightarrow5 \rightarrow8 \rightarrow10 \rightarrow6 \rightarrow1 \rightarrow9\rightarrow7
$$变异操作：采取交换变异（Exchange Mutation，EM）。随机从父代个体中选择两个基因位置，然后互换这两个位置的基因。
进化逆转操作：这个是标准的遗传算法没有的，是为了加速进化而加入的一个操作。这里的进化是指逆转操作具有单向性，即只有逆转之后个体变得更优才会执行逆转操作，否则逆转无效。具体的方法是，随机产生 [ 1 , 10 ] [1,10] [1,10] （这里仍然以10个城市为例）之间的两个随机数 r 1 r_1 r1 和 r 2 r_2 r2 （其实也是允许相同的，只是 r 1 r_1 r1， r 2 r_2 r2相同之后逆转无效，但是这不会经常发生），然后将r1和r2之间的基因进行反向排序。比如对于染色体:
KaTeX parse error: Expected '{', got '[' at position 38: …arrow\textcolor[̲rgb]{1,0,0}{7 \…
随机产生 r 1 = 3 r_1=3 r1=3， r 2 = 5 r_2=5 r2=5 ，它们之间的基因（标红位置）反向排列之后得到的染色体如下:
KaTeX parse error: Expected '{', got '[' at position 38: …arrow\textcolor[̲rgb]{1,0,0}{5 \…
源码
#pragma once
#include <random>
#include <vector>
#include <iostream>
#include <cmath>
#include <ctime>
#include <cstdlib>
#include <bitset>
#include <iomanip>
#include <string>
#include <sstream>
#include <fstream>
#include <time.h>
using namespace std;
int max_gen = 6000; // max generation
const int pop_size = 20; // size of colony
double px = 0.1; // probability of crossover
double pm = 0.01; // probability of mutation
const int cities_num = 280; // number of cities
int generation;
int cur_best; // the index with the best fitness in the colony
struct COORD
{
double x;
double y;
};
struct City
{
int number;
COORD coord;
};
City cityList[cities_num]; // the initial city list read from file
struct Route // a route visits each city once and returns to the origin city
{
City city[cities_num]{}; // city sequence in the route
double length = 0; // length of the route
double fitness = 0;
double rFitness = 0;
};
Route randomRoute; // a random routes creates form the initial city list
Route colony[pop_size+1]; // a colony is an array of routes
Route newColony[pop_size+1];
double best_fitness = 0;
long fitness_ratio = 10000;
long long min_distance = 0;
double fitness_results[10][10][10] = {0};
double min_distance_results[10][10][10] = {0};
long long progress_step = 0; // count the progress steps
/* Declaration of procedures used by this genetic algorithm */
void readCityList(char *);
void init_colony();
double distance(City, City);
double path_len(Route);
void select();
void crossover();
void mutate();
void reverse();
void knuth_durstenfeld_shuffle (City [], int);
void copy_route(Route a, Route b);
void print_route(const Route);
void ga_progress();
void evaluate();
void keepTheBest();
void elitist();
void xOver(Route one, Route two);
void report();
void swap_city(City a, City b);
void cross(Route one, Route two);
void saveToTxt();
/* Implementation of procedures used by this genetic algorithm */
/* read cities from a file */
void readCityList(char *fileName) {
FILE * fp;
fp = fopen(fileName, "r");
if (fp == NULL) {
cout<< "Open file " << fileName << " failed!"<<endl;
exit(1);
}
for (auto & city : cityList) {
fscanf(fp, "%d", &city.number);
fscanf(fp, "%lf", &city.coord.x);
fscanf(fp, "%lf", &city.coord.y);
}
fclose(fp);
}
/* init the colony with pop_size routes */
void init_colony() {
int listLength = sizeof(cityList)/sizeof(cityList[0]);
int num = 0;
while (num < pop_size) {
knuth_durstenfeld_shuffle(cityList, listLength);
for (int i = 0; i < cities_num; ++i) {
randomRoute.city[i] = cityList[i];
colony[num].city[i] = randomRoute.city[i];
colony[num].length = path_len(colony[num]);
colony[num].fitness = 1/colony[num].length;
}
num++;
}
best_fitness = colony[0].fitness;
min_distance = colony[0].length;
// for (const auto & j : colony) {
// print_route(j);
// }
}
/* print the route with length and fitness */
void print_route(const Route route) {
cout<<"route: ";
for (int i = 0; i < cities_num; ++i) {
cout<<route.city[i].number;
if (i!=cities_num-1) {
cout<< "->";
}
}
cout<<" length: "<<route.length<<" ";
cout<<" fitness: "<<route.fitness<<endl;
}
/* knuth_durstenfeld_shuffle: gives me a random sequence of cities */
void knuth_durstenfeld_shuffle (City cityList[], int size)
{
srand(time(0));
for(int i = 0; i < size; i++)
{
// guarantee the i'th value does not interfere with it's front value
int index = i + rand()%(size-i);
swap(cityList[index], cityList[i]);
}
}
/* figure out the distance between tow cities */
double distance(City a, City b) {
return sqrt((a.coord.x-b.coord.x)*(a.coord.x-b.coord.x)+(a.coord.y-b.coord.y)*(a.coord.y-b.coord.y));
}
/* figure out the length of the route */
double path_len(Route route) {
double path_length = 0;
for (int i = 0; i < cities_num-1; ++i) {
path_length += distance(route.city[i], route.city[i+1]);
}
path_length += distance(route.city[0], route.city[cities_num-1]);
return path_length;
}
/* evaluate the fitness of the colony, figure out all the fitness of routes in this colony */
void evaluate() {
int num = 0;
while (num < pop_size) {
colony[num].length = path_len(colony[num]);
colony[num].fitness = 1/colony[num].length;
// cout<<colony[num].fitness<<endl;
num++;
}
for (int i = 0; i < pop_size; ++i) {
if (colony[i].fitness > best_fitness) {
best_fitness = colony[i].fitness;
min_distance = colony[i].length;
}
}
// cout<<"best_fitness: "<<best_fitness<<endl;
}
/**
* Keep_the_best function: This function keeps track of the
* best member of the colony. Note that the last entry in
* the array colony holds a copy of the best individual
*/
void keepTheBest() {
cur_best = 0; /* stores the index of the best individual */
colony[pop_size].fitness = colony[cur_best].fitness;
for (int mem = 0; mem < pop_size; ++mem) {
if (colony[mem].fitness > colony[pop_size].fitness) {
cur_best = mem;
colony[pop_size].fitness = colony[mem].fitness;
}
}
/* once the best member in the colony is found, copy the genes */
for (int i = 0; i < cities_num; i++) {
colony[pop_size].city[i] = colony[cur_best].city[i];
}
}
/****************************************************************/
/* Elitist function: The best member of the previous generation */
/* is stored as the last in the array. If the best member of */
/* the current generation is worse then the best member of the */
/* previous generation, the latter one would replace the worst */
/* member of the current colony */
/****************************************************************/
void elitist()
{
int i;
double best, worst;
int best_mem, worst_mem;
best = colony[0].fitness;
worst = colony[0].fitness;
for (i = 0; i < pop_size ; ++i) // find the best and the worst
{
if (colony[i].fitness >= best)
{
best = colony[i].fitness;
best_mem = i;
}
if (colony[i].fitness <= worst)
{
worst = colony[i].fitness;
worst_mem = i;
}
}
if (best >= colony[pop_size].fitness) // store the best individual as the last in the array.
{
for (i = 0; i < cities_num; i++)
colony[pop_size].city[i] = colony[best_mem].city[i];
colony[pop_size].fitness = colony[best_mem].fitness;
}
else // the
{
for (i = 0; i < cities_num; i++)
colony[worst_mem].city[i] = colony[pop_size].city[i];
colony[worst_mem].fitness = colony[pop_size].fitness;
}
}
/**************************************************************/
/* Selection function: Standard proportional selection for */
/* maximization problems incorporating elitist model - makes */
/* sure that the best member survives */
/**************************************************************/
void select()
{
int mem, i, j;
double sum = 0;
double cFitness;
double p;
for (mem = 0; mem < pop_size; mem++) // counting sum fitness
sum += colony[mem].fitness;
for (mem = 0; mem < pop_size; mem++)
colony[mem].rFitness = colony[mem].fitness / sum; // counting relative fitness
for (i = 0; i < pop_size; i++)
{
cFitness = 0;
p = rand() % 1000 / 1000.0; // roll the roulette : 0 < p < 1
for (j = 0; j < pop_size; j++)
{
cFitness += colony[j].rFitness;
if (p <= cFitness) // roulette
{
newColony[i] = colony[j];
break;
}
}
}
for (i = 0; i < pop_size; i++) // replace the old colony
colony[i] = newColony[i];
}
/***************************************************************/
/* Crossover selection: selects two parents that take part in */
/* the crossover. Implements a single point crossover */
/***************************************************************/
void crossover()
{
int jack, rose;
int couples = 0;
double x;
for (int i = 0; i < pop_size; ++i)
{
x = rand() % 1000 / 1000.0;
if (x < px) // the chance of crossover
{
++couples;
jack = i; // we got jack
if (couples % 2 == 0)
// xOver(colony[jack], colony[rose]);
cross(colony[jack], colony[rose]);
else
rose = i; // we got rose
}
}
}
void xOver(Route one, Route two)
{
int point_1;
int point_2;
/* select crossover point */
if(cities_num > 1)
{
point_1 = (rand() % (cities_num - 1)); // generate random int from 0 to cities_num
point_2 = (rand() % (cities_num - 1));
if (point_1 > point_2) { // make sure point_1 < point_2
int temp = point_1;
point_1 = point_2;
point_2 = temp;
}
if (point_1 = point_2) {
return;
}
for (int i = point_1; i <= point_2; i++) // reverse the gene between point_1 and point_2
{
// City middle[point_2-point_1];
City middle[100];
for (int j = 0; j < point_2-point_1; ++j) {
middle[j] = one.city[i];
}
one.city[i] = two.city[i];
for (int k = 0; k < point_2-point_1; ++k) {
two.city[i] = middle[k];
}
}
/* find the duplicate gene, swap them */
vector<int> conflict_index_in_one;
vector<int> conflict_index_in_two;
for (int l = 0; l < point_1; ++l) {
for (int i = point_1; i <= point_2; ++i) {
if (one.city[l].number == one.city[i].number) {
conflict_index_in_one.push_back(l);
}
if (two.city[l].number == two.city[i].number) {
conflict_index_in_two.push_back(l);
}
}
}
for (int m = point_2+1; m < cities_num; ++m) {
for (int i = point_1; i <= point_2; ++i) {
if (one.city[m].number == one.city[i].number) {
conflict_index_in_one.push_back(m);
}
if (two.city[m].number == two.city[i].number) {
conflict_index_in_two.push_back(m);
}
}
}
for (auto item_1 : conflict_index_in_one) {
for (auto item_2 : conflict_index_in_two) {
int temp = one.city[item_1].number;
one.city[item_1].number = two.city[item_2].number;
two.city[item_2].number = temp;
}
}
}
}
void cross(Route one, Route two)
{
int pos1,pos2;
City temp;
int conflict1[cities_num];
int conflict2[cities_num];
int num1,num2;
int index1,index2;
pos1 = (rand() % (cities_num - 1));
pos2 = (rand() % (cities_num - 1));
if(pos1 > pos2)
{
int t = pos1;
pos1 = pos2;
pos2 = t;
}
for(int j=pos1; j<=pos2; j++)
{
City city_temp = one.city[j];
one.city[j] = two.city[j];
two.city[j] = city_temp;
}
num1 = 0;
num2 = 0;
for(int j =0; j<=pos1-1; j++)
{
for(int k=pos1; k<=pos2; k++)
{
if(one.city[j].number == one.city[k].number)
{
conflict1[num1] = j;
num1++;
}
if(two.city[j].number == two.city[k].number)
{
conflict2[num2] = j;
num2++;
}
}
}
for(int j=pos2+1; j<cities_num; j++)
{
for(int k=pos1; k<=pos2; k++)
{
if(one.city[j].number == one.city[k].number)
{
conflict1[num1] = j;
num1++;
}
if(two.city[j].number == two.city[k].number)
{
conflict2[num2] = j;
num2++;
}
}
}
if((num1 == num2) && num1 > 0)
{
for(int j=0; j<num1; j++)
{
index1 = conflict1[j];
index2 = conflict2[j];
temp = one.city[index1]; // 交换冲突的位置
one.city[index1] = two.city[index2];
two.city[index2] = temp;
}
}
}
void swap_city(City a, City b) {
City t = a;
a = b;
b = t;
}
/**************************************************************/
/* Mutation: Random uniform mutation. A variable selected for */
/* mutation is replaced by a random value between lower and */
/* upper bounds of this variable */
/**************************************************************/
void mutate()
{
double x;
int point_1;
int point_2;
for (int i = 0; i < pop_size; i++) {
x = rand() % 1000 / 1000.0;
if (x < pm) {
point_1 = (rand() % (cities_num - 1)); // generate random int from 0 to cities_num
point_2 = (rand() % (cities_num - 1));
swap_city(colony[i].city[point_1], colony[i].city[point_2]); // swap two city info
}
}
}
void report() {
cout<<"generation: "<<generation<<"\t best_fitness: "<<best_fitness<<"\t min_distance: "<<min_distance<<endl;
}
// reverse
void reverse()
{
int point_1,point_2;
double dis,reverse_dis;
int n;
int flag;
Route reverse_route;
for(int i=0; i<pop_size; i++)
{
flag = 0; // 用于控制本次逆转是否有效
while(flag == 0)
{
point_1 = (rand() % (cities_num - 1)); // generate random int from 0 to cities_num
point_2 = (rand() % (cities_num - 1));
if(point_1 > point_2)
{
int temp = point_1;
point_1 = point_2;
point_2 = temp;
}
if(point_1 < point_2)
{
for(int j=0; j<cities_num; j++)
reverse_route = colony[i];
n = 0;
for(int j=point_1; j<=point_2; j++)
{
reverse_route.city[j] = colony[i].city[point_2-n];
n++;
}
reverse_dis = path_len(reverse_route); // distance after reversal
dis = path_len(colony[i]); // original distance
if(reverse_dis < dis)
{
for(int j=0; j<cities_num; j++)
colony[i] = reverse_route;
}
}
flag = 1;
}
}
}
void saveToTxt() {
for(int gen = 500; gen <= 5000; gen = gen + 500)
{
stringstream ss;
ss<<gen<<".txt";
cout << ss.str() <<endl;
ofstream f_out(ss.str());
if (!f_out.is_open())
{
cerr << "cannot create output file!" <<endl;
exit(0);
}
for(int i=0; i<10; i++)
{
for(int j=0; j<10; j++)
{
f_out<<(fitness_results[i][j][gen/500])<<"\t";
}
f_out<<endl;
}
f_out.close();
}
}
void ga_progress(int px_count, int pm_count) {
srand(time(NULL));
generation = 0;
init_colony();
evaluate();
keepTheBest();
int gen_count = 0; // count the generation 500, 1000, ..., 5000
while(generation <= max_gen) {
generation++;
select();
crossover();
mutate();
reverse();
evaluate();
elitist();
// report();
if (generation % 500 == 0) {
fitness_results[px_count][pm_count][gen_count] += colony[pop_size].fitness;
min_distance_results[px_count][pm_count][gen_count] += colony[pop_size].length;
gen_count++;
}
}
}
int main() {
clock_t begin, end;
double cost;
begin = clock();
readCityList("a280.tsp.txt");
int px_count, pm_count;
px_count = 0;
for (int i = 0; i < 10; ++i) {
pm_count = 0;
for (int j = 0; j < 10; ++j) {
progress_step++;
ga_progress(px_count, pm_count);
pm_count++;
pm = pm + 0.01;
}
px_count++;
px = px + 0.1;
}
saveToTxt();
printf("Success \n");
end = clock();
cost = (double)(end - begin)/CLOCKS_PER_SEC;
printf("constant CLOCKS_PER_SEC is: %ld, time cost is: %lf secs", CLOCKS_PER_SEC, cost);
return 0;
}
//void report()
//{
// int i;
// double best_val; /* best colony fitness */
// double avg; /* avg colony fitness */
// double stddev; /* std. deviation of colony fitness */
// double sum_square; /* sum of square for std. calc */
// double square_sum; /* square of sum for std. calc */
// double sum; /* total colony fitness */
//
// sum = 0.0;
// sum_square = 0.0;
//
// for (i = 0; i < pop_size; i++)
// {
// sum += colony[i].fitness;
// sum_square += colony[i].fitness * colony[i].fitness;
// }
//
// avg = sum/(double)pop_size;
// square_sum = avg * avg * pop_size;
// stddev = sqrt((sum_square - square_sum)/(pop_size - 1));
// best_val = colony[pop_size].fitness;
//
//}