巡回セールスマン問題を焼きなまし法で解いてみた

巡回セールスマン問題を焼きなまし法によりC++で解いてみました。
図1は初期経路、図2は評価関数にExchange（巡回経路中の二つの町の番号の入れ替え）を使った場合の最終経路です。
図3は評価関数に2-Optを使った場合の最終経路です。
2-OptはExchangeに比べて14.6%ぐらい巡回経路が短くなっているようです。
実行するたびに最終経路長も変わってくるし、初期経路や町の数によっても変わってくるので、あくまで参考程度ですが。
図4は評価関数に2-Optを使った場合の途中経過です。
SimulatedAnnealing1.cppは評価関数にExchangeを使ったプログラムであり、SimulatedAnnealing2.cppは評価関数に2-Optを使ったプログラムです。
このような方法で最適化ができてしまうのがちょっと不思議ですが、焼きなまし法はけっこう強力な手法ですね。

こちらが参考にさせていただいたサイトです。
https://github.com/chncyhn/simulated-annealing-tsp

図1　初期経路（n=100、L=5446.6）
![_or-opt_240603_003.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/222722/9f304aae-1dfe-0fca-c50b-c1028b1acfc8.jpeg)

図2　Exchangeの最終経路 (n=100、L=931.84)
![exchange_001.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/222722/c59c9999-394b-7993-bc22-438d2f09f777.jpeg)

図3　2-Optの最終経路 (n=100、L=795.6)
![__20240606-205412.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/222722/89950cb0-242e-532f-b99f-5a03bc70cbbb.jpeg)

図4　途中経過 (2-Opt)
![__20240606-205333_001.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/222722/ce3f91b5-f932-cade-8d06-2dd91c0aacf0.gif)

```c++:SimulatedAnnealing1.cpp
// Simulated Annealing (Exchange)

#include <iostream>
using namespace std;

#define NUM_TOWN 20    // 町の数
#define AREA_SCALE 100.0    // 町が存在するエリアの大きさ（正方形の一辺の長さ）

#define ITERATION 1000000 // 繰り返し回数

double town_map[NUM_TOWN][2];    // 町のx, y座標
int route[NUM_TOWN];    // 巡回経路
double length; // 巡回経路の長さ

int route_min[NUM_TOWN];    // 最小の巡回経路
double length_min;    // 最小の巡回経路の経路長

bool b_auto_cool_rate = true;    // クーリングレートを自動的に決定するかどうかのフラグ

// ルートの中の二つの町の行く順番を交換
void swap(int* route1, int i, int j)
{
	int tmp = route1[i];
	route1[i] = route1[j];
	route1[j] = tmp;
}

// 町の座標をランダムに生成
void create_random_map()
{
	for (int i = 0; i < NUM_TOWN; i++)
		for (int j = 0; j < 2; j++)
			town_map[i][j] = AREA_SCALE * (double)rand() / ((double)RAND_MAX + 1);
}

// 都市間の距離を計算
double calc_dist(int i, int j)
{
	double dx = town_map[i][0] - town_map[j][0];
	double dy = town_map[i][1] - town_map[j][1];
	double dist = sqrt(dx * dx + dy * dy);
	return dist;
}

// 指定された巡回経路の長さを計算
double calc_length(int route1[])
{
	double sum_dist = 0.0;
	for (int i = 0; i < NUM_TOWN; i++)
		sum_dist += calc_dist(route1[i], route1[(i + 1) % NUM_TOWN]);
	return sum_dist;
}

// 巡回経路の初期化
void init_route()
{
	for (int i = 0; i < NUM_TOWN; i++)
		route[i] = i;

	length = calc_length(route);
}

void print_route()
{
	for (int i = 0; i < NUM_TOWN; i++)
		cout << route[i] + 1 << " ";
//		cout << route[i] << " ";
	cout << endl;

	length = calc_length(route);
	cout << "Total Distance: " << length << endl;
}

// 町の座標を出力
void print_route_min()
{
	for (int i = 0; i < NUM_TOWN; i++)
		cout << route_min[i] + 1 << " ";
//		cout << route[i] << " ";
	cout << endl;

	length_min = calc_length(route_min);
	cout << "Total Distance: " << length_min << endl;
}

// 町の座標を出力
void print_town_map()
{
	for (int i = 0; i < NUM_TOWN; i++) {
		cout << route[i] + 1 << " ";
		cout << town_map[route[i]][0] << " " << town_map[route[i]][1] << endl;
	}
}

double get_random()
{
	return ((double)rand() + 0.5) / ((double)RAND_MAX + 1.0);
}

void copy_route(int route_dest[], int route_source[])
{
	for (int i = 0; i < NUM_TOWN; i++)
		route_dest[i] = route_source[i];
}

#define DIST(n,m)  calc_dist(route[(n) % NUM_TOWN], route[(m) % NUM_TOWN])

//int repetition = 100;
double initialT = 1.0;
double finalT = 0.0;
double coolingRate = 0.995;

double drand()
{
	return (double)rand()/(double)RAND_MAX;
}

int shouldChange(double delta, double t)
{
	if (delta <= 0.0)
		return 1;
	if (t == 0.0)
		return 0;
	if (drand() < exp(-delta / t))
		return 1;
	return 0;
}

double sa_exchange(double t, int route1[], double current_length)
{
	int town1 = rand() % NUM_TOWN;
	int town2 = rand() % NUM_TOWN;

	swap(route1, town1, town2);

	double new_length = calc_length(route1);

	if (shouldChange(new_length - current_length, t))
		current_length = new_length;
	else
		swap(route1, town1, town2);

	return current_length;
} 

void calc_sa()
{
	if (b_auto_cool_rate == true) {
		initialT = AREA_SCALE * AREA_SCALE;
		coolingRate = pow(pow((double)AREA_SCALE / NUM_TOWN, 2)
			/ (100.0 * initialT), 1.0 / ITERATION);
	}

	copy_route(route_min, route);

	double length_new = length;
	length_min = length;

	double t = initialT;
	for (int i = 0; i < ITERATION; i++) {
		length_new = sa_exchange(t, route, length);

		if (length_new != length) {
			length = length_new;

			if (length_new < length_min) {
				length_min = length_new;
				copy_route(route_min, route);
			}
		}

//		if (b_auto_cool_rate == false)
			if (t < finalT)
				break;

		t *= coolingRate;
	}

//	if (b_auto_cool_rate == true)
//		finalT  = t;
}

int main()
{
	unsigned int seed = (unsigned)time(NULL);
	cout << "Seed = " << seed << endl;
	srand(seed);

	create_random_map();
	init_route();

	cout << "Number of town: " << NUM_TOWN << endl;

	cout << endl << "Town map" << endl;
	print_town_map();

	cout << endl << "First root" << endl;
	print_route();

	calc_sa();

//	print_route();

	cout << endl << "Last root" << endl;
	print_route_min();

	return 0;
}

```

```c++:SimulatedAnnealing2.cpp
// Simulated Annealing (2-Opt)

#include <iostream>
using namespace std;

#define NUM_TOWN 20    // 町の数
#define AREA_SCALE 100.0    // 町が存在するエリアの大きさ（正方形の一辺の長さ）

#define ITERATION 1000000 // 繰り返し回数

double town_map[NUM_TOWN][2];    // 町のx, y座標
int route[NUM_TOWN];    // 巡回経路
double length; // 巡回経路の長さ

int route_min[NUM_TOWN];    // 最小の巡回経路
double length_min;    // 最小の巡回経路の経路長

bool b_auto_cool_rate = true;    // クーリングレートを自動的に決定するかどうかのフラグ

// ルートの中の二つの町の行く順番を交換
void swap(int* route1, int i, int j)
{
	int tmp = route1[i];
	route1[i] = route1[j];
	route1[j] = tmp;
}

// 町の座標をランダムに生成
void create_random_map()
{
	for (int i = 0; i < NUM_TOWN; i++)
		for (int j = 0; j < 2; j++)
			town_map[i][j] = AREA_SCALE * (double)rand() / ((double)RAND_MAX + 1);
}

// 都市間の距離を計算
double calc_dist(int i, int j)
{
	double dx = town_map[i][0] - town_map[j][0];
	double dy = town_map[i][1] - town_map[j][1];
	double dist = sqrt(dx * dx + dy * dy);
	return dist;
}

// 指定された巡回経路の長さを計算
double calc_length(int route1[])
{
	double sum_dist = 0.0;
	for (int i = 0; i < NUM_TOWN; i++)
		sum_dist += calc_dist(route1[i], route1[(i + 1) % NUM_TOWN]);
	return sum_dist;
}

// 巡回経路の初期化
void init_route()
{
	for (int i = 0; i < NUM_TOWN; i++)
		route[i] = i;

	length = calc_length(route);
}

void print_route()
{
	for (int i = 0; i < NUM_TOWN; i++)
		cout << route[i] + 1 << " ";
//		cout << route[i] << " ";
	cout << endl;

	length = calc_length(route);
	cout << "Total Distance: " << length << endl;
}

// 町の座標を出力
void print_route_min()
{
	for (int i = 0; i < NUM_TOWN; i++)
		cout << route_min[i] + 1 << " ";
//		cout << route[i] << " ";
	cout << endl;

	length_min = calc_length(route_min);
	cout << "Total Distance: " << length_min << endl;
}

// 町の座標を出力
void print_town_map()
{
	for (int i = 0; i < NUM_TOWN; i++) {
		cout << route[i] + 1 << " ";
		cout << town_map[route[i]][0] << " " << town_map[route[i]][1] << endl;
	}
}

double get_random()
{
	return ((double)rand() + 0.5) / ((double)RAND_MAX + 1.0);
}

void copy_route(int route_dest[], int route_source[])
{
	for (int i = 0; i < NUM_TOWN; i++)
		route_dest[i] = route_source[i];
}

#define DIST(n,m)  calc_dist(route[(n) % NUM_TOWN], route[(m) % NUM_TOWN])

//int repetition = 100;
double initialT = 1.0;
double finalT = 0.0;
double coolingRate = 0.995;

double drand()
{
	return (double)rand()/(double)RAND_MAX;
}

int shouldChange(double delta, double t)
{
	if (delta <= 0.0)
		return 1;
	if (t == 0.0)
		return 0;
	if (drand() < exp(-delta / t))
		return 1;
	return 0;
}

double sa_2opt(double t, int route1[], double current_length)
{
	int i = rand() % NUM_TOWN;
	int j = rand() % (NUM_TOWN - 3) + i + 2;
	double diff = DIST(i, i + 1) + DIST(j, j + 1) - (DIST(i, j) + DIST(i + 1, j + 1));

	if (shouldChange(-diff, t)) {
		// 部分経路を逆順に並び替える
		for (int k = 0; k < (j - i) / 2; k++)
			swap(route1, (i + 1 + k) % NUM_TOWN, (j - k) % NUM_TOWN);
		current_length = calc_length(route1);
	}

	return current_length;
} 

void calc_sa()
{
	if (b_auto_cool_rate == true) {
		initialT = AREA_SCALE * AREA_SCALE;
		coolingRate = pow(pow((double)AREA_SCALE / NUM_TOWN, 2)
			/ (100.0 * initialT), 1.0 / ITERATION);
	}

	copy_route(route_min, route);

	double length_new = length;
	length_min = length;

	double t = initialT;
	for (int i = 0; i < ITERATION; i++) {
		length_new = sa_2opt(t, route, length);

		if (length_new != length) {
			length = length_new;

			if (length_new < length_min) {
				length_min = length_new;
				copy_route(route_min, route);
			}
		}

//		if (b_auto_cool_rate == false)
			if (t < finalT)
				break;

		t *= coolingRate;
	}

//	if (b_auto_cool_rate == true)
//		finalT  = t;
}

int main()
{
	unsigned int seed = (unsigned)time(NULL);
	cout << "Seed = " << seed << endl;
	srand(seed);

	create_random_map();
	init_route();

	cout << "Number of town: " << NUM_TOWN << endl;

	cout << endl << "Town map" << endl;
	print_town_map();

	cout << endl << "First root" << endl;
	print_route();

	calc_sa();

//	print_route();

	cout << endl << "Last root" << endl;
	print_route_min();

	return 0;
}

```

```title:実行結果 (sa + 2-Opt)
Seed = 1722772335
Number of town: 20

Town map
1 30.4138 4.86755
2 98.941 36.8439
3 95.5261 54.0039
4 69.8853 13.681
5 61.8195 67.2607
6 86.087 20.0012
7 9.53064 29.7546
8 64.7186 36.7432
9 48.4497 4.36096
10 90.4327 40.097
11 9.61609 69.104
12 65.7379 62.8021
13 23.9502 93.7653
14 43.6462 8.77075
15 40.9241 96.3898
16 20.6482 86.618
17 90.3809 70.3064
18 18.1732 77.8168
19 40.6555 22.2748
20 30.6427 5.00488

First root
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Total Distance: 1088.51

Last root
4 8 19 14 9 20 1 7 11 18 16 13 15 5 12 17 3 10 2 6
Total Distance: 355.527
```