巡回セールスマン問題を蟻コロニー最適化によりC++で解いてみました。
図1は初期経路、図2は蟻コロニー最適化による最終経路です。
蟻コロニー最適化は局所解に陥りやすく、局所解から抜け出すことが難しかったため、2-Optを併用するアルゴリズムも実装してみました。
図3が2-Optを併用したアルゴリズムの最終経路です。
蟻コロニー最適化の単独のものに比べて、11%ほど巡回経路長を短縮できました。
AntColony1.cppは蟻コロニー最適化の単独のプログラムであり、AntColony2.cppは蟻コロニー最適化と2-Optを併用したプログラムです。
なお、最初から2-Optを併用するよりは、蟻コロニー最適化の単独で十分に最適化を実行してから2-Optを併用したほうが、巡回経路長は短くなる傾向にあるようです。
こちらが参考にさせていただいたサイトです。アルゴリズム的には参考にさせていただいたサイトとほぼ同じです。
巡回セールスマン問題を蟻コロニー最適化と遺伝的アルゴリズムで解いてみる[Python]
https://tony-mooori.blogspot.com/2016/02/tsppython.html
図1 初期経路（n=100、L=5446.6）
図2 最終経路（Ant Colony、n=100、L=879.4）
図3 最終経路（Ant Colony + 2-Opt、n=100、L=780.0）
AntColony1.cpp
// Ant Colony Optimization
#include <iostream>
using namespace std;
#define NUM_TOWN 20 // 町の数
#define AREA_SCALE 100.0 // 町が存在するエリアの大きさ（正方形の一辺の長さ）
#define ITERATION 1000 // 繰り返し回数
double town_map[NUM_TOWN][2]; // 町のx, y座標
int route[NUM_TOWN]; // 巡回経路
double length; // 巡回経路の長さ
// ルートの中の二つの町の行く順番を交換
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_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;
}
}
// 0より大きく1より小さな乱数の生成
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];
}
double ALPHA = 1.0; // フェロモンの優先度
double BETA = 5.0; // ヒューリスティック情報(距離)の優先度
double Q = 1.0; // フェロモン変化量の係数
double VANISH_RATIO = 0.95; // 蒸発率
class TSP_ANT{
protected:
double** weight; // フェロモンの量、[NUM_TOWN][NUM_TOWN]
double** delta; // フェロモン変化量、[NUM_TOWN][NUM_TOWN]
int* route_new1; // 巡回経路、[NUM_TOWN]
public:
TSP_ANT(); // 初期化を行う関数
~TSP_ANT();
void initialize();
// 任意の確率分布に従って乱数を生成する関数
int random_index(int n_percentage, double percentage[]);
void solve_new_route(int route_new[]);
double solve_pheromone(int route_new[]);
double solve(); // 巡回セールスマン問題を蟻コロニー最適化で解く
int* get_route_new() { return route_new1; }
};
TSP_ANT* pTspAnt = nullptr;
// 初期化
TSP_ANT::TSP_ANT()
{
weight = new double*[NUM_TOWN];
for (int i = 0; i < NUM_TOWN; i++)
weight[i] = new double[NUM_TOWN];
delta = new double*[NUM_TOWN];
for (int i = 0; i < NUM_TOWN; i++)
delta[i] = new double[NUM_TOWN];
route_new1 = new int[NUM_TOWN];
initialize();
}
TSP_ANT::~TSP_ANT()
{
for (int i = 0; i < NUM_TOWN; i++)
delete []weight[i];
delete []weight;
for (int i = 0; i < NUM_TOWN; i++)
delete []delta[i];
delete []delta;
delete []route_new1;
}
void TSP_ANT::initialize()
{
// フェロモン量の初期化
for (int i = 0; i < NUM_TOWN; i++)
for (int j = 0; j < NUM_TOWN; j++)
// weight[i][j] = get_random();
weight[i][j] = 0.0;
// フェロモンの変化量を初期化
for (int j = 0; j < NUM_TOWN; j++)
for (int k = 0; k < NUM_TOWN; k++)
delta[j][k] = 0.0;
// 巡回経路の初期化
for (int i = 0; i < NUM_TOWN; i++)
route_new1[i] = i;
}
// 任意の確率分布に従って乱数を生成
int TSP_ANT::random_index(int n_percentage, double percentage[])
{
int index;
while (true) {
index = rand() % n_percentage;
double y = get_random();
if (y < percentage[index])
break;
}
return index;
}
void TSP_ANT::solve_new_route(int route_new[])
{
int* city = new int[NUM_TOWN];
double* upper = new double[NUM_TOWN];
double* evaluation = new double[NUM_TOWN]; // 評価関数
double* percentage = new double[NUM_TOWN]; // 移動確率
for (int j = 0; j < NUM_TOWN; j++)
city[j] = j;
int now_city = rand() % NUM_TOWN; // 現在居る都市番号
int other_city = NUM_TOWN - 1; // 現在都市を除いた残りの都市の数
for (int j = 0; j < other_city; j++)
if (city[j] >= now_city) {
city[j] = city[j + 1];
}
route_new[0] = now_city;
for (int j = 1; j < NUM_TOWN; j++) {
double sum_upper = 0.0;
for (int k = 0; k < other_city; k++) {
upper[k] = pow(weight[now_city][city[k]], ALPHA)
* pow(calc_dist(now_city, city[k]), -1.0 * BETA);
sum_upper += upper[k];
}
if (sum_upper == 0.0)
sum_upper = DBL_MIN; // ゼロ除算回避のためDBL_MINを代入
double sum_eva = 0.0;
for (int k = 0; k < other_city; k++) {
evaluation[k] = upper[k] / sum_upper; // 評価関数
if (evaluation[k] == 0.0)
evaluation[k] = DBL_MIN; // ゼロ除算回避のためDBL_MINを代入
sum_eva += evaluation[k];
}
for (int k = 0; k < other_city; k++) {
percentage[k] = evaluation[k] / sum_eva; // 移動確率
}
int index = random_index(other_city, percentage); // 移動先の要素番号取得
// 状態の更新
now_city = city[index];
other_city--;
for (int k = 0; k < other_city; k++)
if (city[k] >= now_city) {
city[k] = city[k + 1];
}
route_new[j] = now_city;
}
delete []city;
delete []upper;
delete []evaluation;
delete []percentage;
}
double TSP_ANT::solve_pheromone(int route_new[])
{
double L = calc_length(route_new); // 経路のコストを計算
// フェロモンの変化量を初期化
for (int j = 0; j < NUM_TOWN; j++)
for (int k = 0; k < NUM_TOWN; k++)
delta[j][k] = 0.0;
double c = Q / L;
for (int j = 0; j < (NUM_TOWN - 1); j++) {
delta[route_new[j]][route_new[j + 1]] = c;
delta[route_new[j + 1]][route_new[j]] = c;
}
// フェロモン更新
for (int j = 0; j < NUM_TOWN; j++) {
for (int k = 0; k < NUM_TOWN; k++) {
weight[j][k] = weight[j][k] * VANISH_RATIO + delta[j][k];
}
}
return L;
}
// 巡回セールスマン問題を蟻コロニー最適化で解く
double TSP_ANT::solve()
{
solve_new_route(route_new1);
double L = solve_pheromone(route_new1);
return L;
}
void calc_ant()
{
for (int i = 0; i < ITERATION; i++) {
double length_new = pTspAnt->solve();
if (length_new < length) {
length = length_new;
copy_route(route, pTspAnt->get_route_new());
}
}
}
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();
pTspAnt = new TSP_ANT;
calc_ant();
cout << endl << "Last root" << endl;
print_route();
delete pTspAnt;
return 0;
}
AntColony2.cpp
// Ant Colony Optimization + 2-Opt
#include <iostream>
using namespace std;
#define NUM_TOWN 20 // 町の数
#define AREA_SCALE 100.0 // 町が存在するエリアの大きさ（正方形の一辺の長さ）
#define ITERATION 1000 // 繰り返し回数
double town_map[NUM_TOWN][2]; // 町のx, y座標
int route[NUM_TOWN]; // 巡回経路
double length; // 巡回経路の長さ
// ルートの中の二つの町の行く順番を交換
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_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;
}
}
// 0より大きく1より小さな乱数の生成
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];
}
double ALPHA = 1.0; // フェロモンの優先度
double BETA = 5.0; // ヒューリスティック情報(距離)の優先度
double Q = 1.0; // フェロモン変化量の係数
double VANISH_RATIO = 0.95; // 蒸発率
class TSP_ANT{
protected:
double** weight; // フェロモンの量、[NUM_TOWN][NUM_TOWN]
double** delta; // フェロモン変化量、[NUM_TOWN][NUM_TOWN]
int* route_new1; // 巡回経路、[NUM_TOWN]
public:
TSP_ANT(); // 初期化を行う関数
~TSP_ANT();
void initialize();
// 任意の確率分布に従って乱数を生成する関数
int random_index(int n_percentage, double percentage[]);
void solve_new_route(int route_new[]);
double solve_pheromone(int route_new[]);
double solve(); // 巡回セールスマン問題を蟻コロニー最適化で解く
bool ant_two_opt(int route[]);
int* get_route_new() { return route_new1; }
};
TSP_ANT* pTspAnt = nullptr;
// 初期化
TSP_ANT::TSP_ANT()
{
weight = new double*[NUM_TOWN];
for (int i = 0; i < NUM_TOWN; i++)
weight[i] = new double[NUM_TOWN];
delta = new double*[NUM_TOWN];
for (int i = 0; i < NUM_TOWN; i++)
delta[i] = new double[NUM_TOWN];
route_new1 = new int[NUM_TOWN];
initialize();
}
TSP_ANT::~TSP_ANT()
{
for (int i = 0; i < NUM_TOWN; i++)
delete []weight[i];
delete []weight;
for (int i = 0; i < NUM_TOWN; i++)
delete []delta[i];
delete []delta;
delete []route_new1;
}
void TSP_ANT::initialize()
{
// フェロモン量の初期化
for (int i = 0; i < NUM_TOWN; i++)
for (int j = 0; j < NUM_TOWN; j++)
// weight[i][j] = get_random();
weight[i][j] = 0.0;
// フェロモンの変化量を初期化
for (int j = 0; j < NUM_TOWN; j++)
for (int k = 0; k < NUM_TOWN; k++)
delta[j][k] = 0.0;
// 巡回経路の初期化
for (int i = 0; i < NUM_TOWN; i++)
route_new1[i] = i;
}
// 任意の確率分布に従って乱数を生成する関数
int TSP_ANT::random_index(int n_percentage, double percentage[])
{
int index;
while (true) {
index = rand() % n_percentage;
double y = get_random();
if (y < percentage[index])
break;
}
return index;
}
void TSP_ANT::solve_new_route(int route_new[])
{
int* city = new int[NUM_TOWN];
double* upper = new double[NUM_TOWN];
double* evaluation = new double[NUM_TOWN]; // 評価関数
double* percentage = new double[NUM_TOWN]; // 移動確率
for (int j = 0; j < NUM_TOWN; j++)
city[j] = j;
int now_city = rand() % NUM_TOWN; // 現在居る都市番号
int other_city = NUM_TOWN - 1; // 現在都市を除いた残りの都市の数
for (int j = 0; j < other_city; j++)
if (city[j] >= now_city) {
city[j] = city[j + 1];
}
route_new[0] = now_city;
for (int j = 1; j < NUM_TOWN; j++) {
double sum_upper = 0.0;
for (int k = 0; k < other_city; k++) {
upper[k] = pow(weight[now_city][city[k]], ALPHA)
* pow(calc_dist(now_city, city[k]), -1.0 * BETA);
sum_upper += upper[k];
}
if (sum_upper == 0.0)
sum_upper = DBL_MIN; // ゼロ除算回避のためDBL_MINを代入
double sum_eva = 0.0;
for (int k = 0; k < other_city; k++) {
evaluation[k] = upper[k] / sum_upper; // 評価関数
if (evaluation[k] == 0.0)
evaluation[k] = DBL_MIN; // ゼロ除算回避のためDBL_MINを代入
sum_eva += evaluation[k];
}
for (int k = 0; k < other_city; k++) {
percentage[k] = evaluation[k] / sum_eva; // 移動確率
}
int index = random_index(other_city, percentage); // 移動先の要素番号取得
// 状態の更新
now_city = city[index];
other_city--;
for (int k = 0; k < other_city; k++)
if (city[k] >= now_city) {
city[k] = city[k + 1];
}
route_new[j] = now_city;
}
delete []city;
delete []upper;
delete []evaluation;
delete []percentage;
}
double TSP_ANT::solve_pheromone(int route_new[])
{
double L = calc_length(route_new); // 経路のコストを計算
// フェロモンの変化量を初期化
for (int j = 0; j < NUM_TOWN; j++)
for (int k = 0; k < NUM_TOWN; k++)
delta[j][k] = 0.0;
double c = Q / L;
for (int j = 0; j < (NUM_TOWN - 1); j++) {
delta[route_new[j]][route_new[j + 1]] = c;
delta[route_new[j + 1]][route_new[j]] = c;
}
// フェロモン更新
for (int j = 0; j < NUM_TOWN; j++) {
for (int k = 0; k < NUM_TOWN; k++) {
weight[j][k] = weight[j][k] * VANISH_RATIO + delta[j][k];
}
}
return L;
}
// 巡回セールスマン問題を蟻コロニー最適化で解く
double TSP_ANT::solve()
{
solve_new_route(route_new1);
ant_two_opt(route_new1); // 2-Optを実行
double L = solve_pheromone(route_new1);
return L;
}
#define DIST(n,m) calc_dist(route[(n) % NUM_TOWN], route[(m) % NUM_TOWN])
// 2-OPTで経路修正
bool TSP_ANT::ant_two_opt(int route[])
{
bool changed = false;
bool fImproved = true;
while (fImproved == true) {
fImproved = false;
for (int i = 0; i < NUM_TOWN; i++) {
for (int j = i + 2; j < (i + NUM_TOWN - 1); j++) {
double diff = DIST(i, i + 1) + DIST(j, j + 1)
- (DIST(i, j) + DIST(i + 1, j + 1));
if (diff > 0) {
// 部分経路を逆順に並び替える
for (int k = 0; k < ((j - i) / 2); k++)
swap(route, (i + 1 + k) % NUM_TOWN, (j - k) % NUM_TOWN);
fImproved = true; // 改善した
changed = true;
}
}
}
}
return changed;
}
void calc_ant()
{
for (int i = 0; i < ITERATION; i++) {
double length_new = pTspAnt->solve();
if (length_new < length) {
length = length_new;
copy_route(route, pTspAnt->get_route_new());
}
}
}
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();
pTspAnt = new TSP_ANT;
calc_ant();
cout << endl << "Last root" << endl;
print_route();
delete pTspAnt;
return 0;
}
実行結果 (Ant Colony + 2-Opt)
Seed = 1722745481
Number of town: 20
Town map
1 62.793 52.4963
2 96.0693 2.45972
3 96.5118 51.0132
4 98.5992 9.43298
5 5.55115 13.1195
6 7.27844 10.2966
7 86.5112 26.8219
8 74.7498 47.1832
9 62.8448 22.1863
10 59.4971 48.5779
11 87.2375 76.767
12 72.0734 3.2074
13 82.4188 5.61218
14 67.5659 95.0623
15 98.6115 1.46484
16 87.5183 90.1855
17 15.686 40.7837
18 4.37927 15.8112
19 68.1732 21.1945
20 71.6064 22.5647
First root
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Total Distance: 1022.94
Last root
19 9 6 5 18 17 10 1 8 14 16 11 3 7 4 15 2 13 12 20
Total Distance: 374.098