TSP problem
TSP (Traveling Salesman Problem) is a typical NP complete problem, that is, its worst-case time complexity increases exponentially with the increase of problem scale. TSP problem can be described as: given the mutual distance between n cities, a traveler starts from a city, visits each city once and only once, and finally returns to the starting city. How to arrange the shortest route.
The mathematical model can be expressed as: searching for natural subsets X = { 1 , 2 , . . . , n } X=\{1,2,...,n\} X={1,2,...,n} (n indicates the number of cities) P ( X ) = { V 1 , V 2 , . . . . , V n } P(X)=\{V1,V2,....,Vn\} P (x) = {V1, V2,..., VN}, so that 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) takes the minimum value, where d ( V i , V i + 1 ) d(V_i,V_{i+1}) d(Vi, Vi+1) indicates the city V i V_i Vi. to V i + 1 V_{i+1} The distance between Vi+1 +.
Genetic algorithm steps
-
Coding problem: since this is a discrete problem, we use integer coding 1 1 1 ~ n n n to represent n n n cities, 1 1 1 ~ n n Any permutation of n constitutes a solution to the problem. You know, for n n TSP problems in n cities, a total of n ! n! n! Two different routes.
-
Population initialization: for N N The population of N individuals is given randomly N N The solution of N problems (equivalent to chromosome) is used as the initial population. The specific methods adopted here are: 1 , 2 ， . . . , n 1,2，...,n 1,2，..., n as the first individual, and then 2 , 3 ， . . n 2,3，..n 2,3，.. n and 1 1 1 exchange position obtained n − 1 n-1 n − 1 solution, from 2 2 2 start, 3 , 4 , . . . , n 3,4,...,n 3,4,...,n and 2 2 2. The exchange position is obtained n − 2 n-2 n − 2 solutions, and so on. (if this is not enough for the initial population, it can be considered again n , n − 1 , . . . , 1 n,n-1,...,1 n,n−1,..., 1 this sequence, and then generate it in the same way)
-
Fitness function: let a solution traverse the total distance of the initial walk as D, then the fitness f i t n e s s = 1 D fitness=\frac{1}{D} fitness=D1 . That is, the higher the total distance, the lower the fitness, the lower the total distance (the better the solution), and the higher the fitness. The fitness value here can take the opposite number of the total distance.
-
Selection operation: the probability of an individual being selected is directly proportional to the fitness. The higher the fitness, the greater the probability of an individual being selected. Roulette is still used here.
-
Crossover operation: crossover operation is the most important operation of genetic algorithm and the main source of generating new individuals. It is directly related to the global optimization ability of the algorithm. Here, Partial Mapping Crossover (PMX) is adopted. For example, for n = 10 n=10 In the case of n=10, for two paths:
$$
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: randomly generate two$ ̲ Random number between [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
purple colour Department branch surface show hand over fork too come of base because ， this individual Time Wait meeting hair present can can hand over fork too come of base because And primary come his he position Set upper of base because have heavy complex ， Allow easy straight reach ， The first one individual individual body heavy complex base because of number order And The first two individual individual body heavy complex base because of number order yes mutually with of ( this in all yes 3 individual ) ， only need want hold The first one individual individual body heavy complex of base because ( use green colour mark knowledge ) And The first two individual individual body heavy complex of base because do hand over change ， Namely can with eliminate except punching Outburst . eliminate except punching Outburst of after of solution as lower : The purple part indicates the crossed genes. At this time, it will be found that the possibly crossed genes are duplicated with the original genes in other positions. It is easy to find that the number of duplicated genes in the first individual is the same as that in the second individual (there are three here), The conflict can be eliminated by simply exchanging the duplicated gene of the first individual (marked with green) with the duplicated gene of the second individual. The solution after eliminating the conflict is as follows: The purple part indicates the crossed genes. At this time, it will be found that the possible crossed genes are duplicated with the original genes in other positions. It is easy to find that the number of repeated genes in the first individual is the same as that in the second individual (there are three here), The conflict can be eliminated by simply exchanging the duplicated gene of the first individual (marked with green) with the duplicated gene of the second individual. The solution after eliminating the conflict is as follows:
1 \rightarrow9 \rightarrow7 \rightarrow6 \rightarrow5 \rightarrow3 \rightarrow2 \rightarrow10\rightarrow 8 \rightarrow4
\
3\rightarrow2 \rightarrow4 \rightarrow5 \rightarrow8 \rightarrow10 \rightarrow6 \rightarrow1 \rightarrow9\rightarrow7
$$ -
Mutation operation: Exchange Mutation (EM) is adopted. Randomly select two gene positions from the parent individuals, and then exchange the genes at these two positions.
-
Evolution reversal operation: This is not available in the standard genetic algorithm. It is an operation added to accelerate evolution. Evolution here means that the reversal operation is unidirectional, that is, the reversal operation will be performed only when the individual becomes better after reversal, otherwise the reversal is invalid. The specific method is random generation [ 1 , 10 ] [1,10] Two random numbers between [1,10] (still taking 10 cities as an example) r 1 r_1 r1 ＾ and r 2 r_2 r2 (in fact, the same is allowed, but r 1 r_1 r1， r 2 r_2 r2 ＾ the reversal is ineffective after the same, but this does not often occur), and then the genes between r1 and r2 are reverse sequenced. For example, for chromosomes:
KaTeX parse error: Expected '{', got '[' at position 38: ...arrow\textcolor[̲rgb]{1,0,0}{7 \...
Random generation r 1 = 3 r_1=3 r1=3， r 2 = 5 r_2=5 r2 = 5, the chromosomes obtained after the reverse arrangement of genes (marked in red) between them are as follows:
KaTeX parse error: Expected '{', got '[' at position 38: ...arrow\textcolor[̲rgb]{1,0,0}{5 \...
Source code
#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]; // Swap conflicting locations
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; // Used to control whether this reversal is effective
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;
//
//}