著者:佐藤 楓真
今回は、「焼きなまし法」と呼ばれる手法で、すべての都市を1回訪問して出発地点に帰るまでの最短経路を求める「巡回セールスマン問題」を近似的に解いてみます。焼きなまし法のベースとなる「山登り法」についても解説します。使用するプログラミング言語は、C++です。
シェルスクリプトマガジン Vol.91は以下のリンク先でご購入できます。![]()
![]()
図3 図1の巡回セールスマン問題を山登り法で解くC++プログラムの例
#include <iostream>
#include <algorithm>
using namespace std;
#define TL 3 //制限時間
#define N 4
int d[N][N] = {{0, 2, 10, 5},
{2, 0, 5, 3},
{10, 5, 0, 4},
{5, 3, 4, 0}};
int ans[N+1] = {0, 1, 3, 2, 0};
int rand_num(void) {
return rand() % N + 1;
}
int calc_dist(void) {
int sum = 0;
for (int i = 0; i < N; i++) {
sum += d[ans[i]][ans[i+1]];
}
return sum;
}
int main(void) {
srand((unsigned)time(NULL));
int st_time = clock();
while ((clock() - st_time) / CLOCKS_PER_SEC < TL) {
int L = rand_num(), R = rand_num();
if (L > R) swap(L, R);
int bf_dist = calc_dist();
reverse(ans + L, ans + R);
int af_dist = calc_dist();
if (bf_dist < af_dist) reverse(ans + L, ans + R);
}
for (int i = 0; i < N+1; i++) {
cout << ans[i] + 1 << " ";
}
cout << endl << calc_dist() << endl;
return 0;
}
図4 都市数「5000」のテスト用データを作成するC++プログラム
#include <iostream>
#include <vector>
#include <fstream>
#include <string>
using namespace std;
#define N 5000
int main(void) {
srand((unsigned)time(NULL));
for (int num = 1; num <= 100; num++) {
vector<vector<int>> d(N, vector<int>(N));
for (int i = 0; i < N; i++) {
for (int j = i; j < N; j++) {
if (i == j) {
d[i][j] = 0;
} else {
d[i][j] = rand() % 5000 + 1;
}
}
}
for (int i = 0; i < N; i++) {
for (int j = 0; j < i; j++) {
d[i][j] = d[j][i];
}
}
string name = "case" + to_string(num) + ".txt";
ofstream out(name);
out << N << endl;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
out << d[i][j] << " ";
}
out << endl;
}
}
return 0;
}
図5 テスト用データに基づく巡回セールスマン問題を山登り法で解くC++プログラムの例
#include <iostream>
#include <algorithm>
#include <vector>
#include <string>
#include <fstream>
using namespace std;
#define TL 3 //制限時間
int calc_dist(vector<vector<int>>& dist, vector<int>& ans, int N) {
int sum = 0;
for (int i = 0; i < N; i++) {
sum += dist[ans[i]][ans[i+1]];
}
return sum;
}
int rand_num(int N){
return rand() % N + 1;
}
int main(void) {
long long ave = 0;
for (int num = 1; num <= 100; num++) {
srand((unsigned)time(NULL));
string name = "case" + to_string(num) + ".txt";
ifstream infile(name);
int N;
infile >> N;
vector<vector<int>> d(N, vector<int>(N));
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
infile >> d[i][j];
}
}
vector<int> ans(N + 1);
for (int i = 0; i < N; i++) ans[i] = i;
ans[N] = 0;
int st_time = clock();
while ((clock() - st_time) / CLOCKS_PER_SEC < TL) {
int bf_dist = calc_dist(d, ans, N);
int L = rand_num(N), R = rand_num(N);
if (L > R) swap(L, R);
reverse(ans.begin() + L, ans.begin() + R);
int af_dist = calc_dist(d, ans, N);
if (af_dist > bf_dist) reverse(ans.begin() + L, ans.begin() + R);
}
ave += calc_dist(d, ans, N);
}
cout << ave / 100.0 << endl;
return 0;
}
図7 テスト用データに基づく巡回セールスマン問題を焼きなまし法で解くC++プログラムの例
#include <iostream>
#include <algorithm>
#include <vector>
#include <string>
#include <fstream>
#include <cmath>
using namespace std;
#define MAX_TMP 100.0 //初期温度
#define TL 3 //制限時間
int calc_dist(vector<vector<int>>& dist, vector<int>& ans, int N) {
int sum = 0;
for (int i = 0; i < N; i++) {
sum += dist[ans[i]][ans[i+1]];
}
return sum;
}
int rand_num(int N) {
return rand() % N + 1;
}
double rand_p(void) {
return (double)rand() / RAND_MAX;
}
int main(void) {
srand((unsigned)time(NULL));
long long ave = 0;
for (int num = 1; num <= 100; num++) {
string name = "case" + to_string(num) + ".txt";
ifstream infile(name);
int N;
infile >> N;
vector<vector<int>> d(N, vector<int>(N));
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
infile >> d[i][j];
}
}
vector<int> ans(N + 1);
for (int i = 0; i < N; i++) ans[i] = i;
ans[N] = 0;
int st_time = clock(); //焼きなまし開始時間
while ((clock() - st_time) / CLOCKS_PER_SEC < TL ) {
int bf_dist = calc_dist(d, ans, N);
int L = rand_num(N), R = rand_num(N);
if (L > R) swap(L, R);
reverse(ans.begin() + L, ans.begin() + R);
int af_dist = calc_dist(d, ans, N);
double t = (double)(clock() - st_time) / CLOCKS_PER_SEC; //経過時間
double tmp = MAX_TMP - MAX_TMP * (t / TL); //現在の温度
double p = exp((bf_dist - af_dist) / tmp); //採用確率
if (rand_p() > p) reverse(ans.begin() + L, ans.begin() + R);
}
ave += calc_dist(d, ans, N);
}
cout << ave / 100.0 <<endl;
return 0;
}