該文章的代碼主要轉自模拟退火算法1.
該文對模拟退火算法作了較好的分析,不過該文中舉例的TSP的代碼有一些問題,我對此作了修正,并在文中最後做出解釋。
代碼如下:
#include <iostream>
#include <string.h>
#include <stdlib.h>
#include <algorithm>
#include <stdio.h>
#include <time.h>
#include <math.h>
#define N 30 //城市數量
#define T 3000 //初始溫度
#define EPS 1e-8 //終止溫度
#define DELTA 0.98 //溫度衰減率
#define LIMIT 10000 //機率選擇上限
#define OLOOP 1000 //外循環次數
#define ILOOP 15000 //内循環次數
using namespace std;
//定義路線結構體
struct Path
{
int citys[N];
double len;
};
//定義城市點坐标
struct Point
{
double x, y;
};
Path path; //記錄最優路徑
Point p[N]; //每個城市的坐标
double w[N][N]; //兩兩城市之間路徑長度
int nCase; //測試次數
double dist(Point A, Point B)
{
return sqrt((A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}
void GetDist(Point p[], int n)
{
for (int i = ; i < n; i++)
for (int j = i + ; j < n; j++)
w[i][j] = w[j][i] = dist(p[i], p[j]);
}
void Input(Point p[], int &n)
{
scanf("%d", &n);
for (int i = ; i < n; i++)
scanf("%lf %lf", &p[i].x, &p[i].y);
}
void Init(int n)
{
nCase = ;
path.len = ;
for (int i = ; i < n; i++)
{
path.citys[i] = i;
if (i != n - )
{
printf("%d--->", i);
path.len += w[i][i + ];
}
else
printf("%d\n", i);
}
printf("\nInit path length is : %.3lf\n", path.len);
}
void Print(Path t, int n)
{
printf("Path is : ");
for (int i = ; i < n; i++)
{
if (i != n - )
printf("%d-->", t.citys[i]);
else
printf("%d\n", t.citys[i]);
}
printf("\nThe path length is : %.3lf\n", t.len);
}
// 随機交換2個城市的位置
Path GetNext(Path p, int n)
{
// Modify by DD: 確定x!=y
int x = , y = ;
while (x == y)
{
x = (int)(n * (rand() / (RAND_MAX + )));
y = (int)(n * (rand() / (RAND_MAX + )));
}
Path ans = p;
swap(ans.citys[x], ans.citys[y]);
ans.len = ;
for (int i = ; i < n - ; i++)
ans.len += w[ans.citys[i]][ans.citys[i + ]];
//cout << "nCase = " << nCase << endl;
//Print(ans, n);
nCase++;
return ans;
}
// 模拟退火
void SA(int n)
{
double t = T;
srand(time(NULL));
Path curPath = path;
Path newPath = path;
int P_L = ; // 連續找到更差結果的次數
int P_F = ; // while循環次數
while () //外循環,主要更新參數t,模拟退火過程
{
for (int i = ; i < ILOOP; i++) //内循環,尋找在一定溫度下的最優值
{
newPath = GetNext(curPath, n);
double dE = newPath.len - curPath.len;
if (dE < ) //如果找到更優值,直接更新
{
curPath = newPath;
P_L = ;
P_F = ;
}
else
{
double rd = rand() / (RAND_MAX + );
// Modify by DD: dE取負數才有可能接受更差解,否則e>1
double e = exp(-dE / t);
if ( e > rd && e < ) //如果找到比目前更差的解,以一定機率接受該解,并且這個機率會越來越小
curPath = newPath;
P_L++;
}
if (P_L > LIMIT)
{
P_F++;
break;
}
}
// Modify by DD: 記錄全局最優解
if (curPath.len < path.len)
path = curPath;
if (P_F > OLOOP || t < EPS)
break;
t *= DELTA;
}
}
int main()
{
freopen("TSP.txt", "r", stdin);
int n;
Input(p, n);
GetDist(p, n);
Init(n);
SA(n);
Print(path, n);
printf("Total test times is : %d\n", nCase);
return ;
}
資料檔案如下:
27
41 94
37 84
53 67
25 62
7 64
2 99
68 58
71 44
54 62
83 69
64 60
18 54
22 60
83 46
91 38
25 38
24 42
58 69
71 71
74 78
87 76
18 40
13 40
82 7
62 32
58 35
45 21
修改過的地方,代碼中寫了注釋。
其中幾個重要的、需要解釋的地方有如下幾個:
-
記錄全局最優解
記錄最優解時,原文使用的是
if (curPath.len < newPath.len)
path = curPath;
這樣顯然記錄的是目前溫度中的最優解。更合理的做法是記錄全局出現過的最優解,即
if(curPath.len < path .len)
path= curPath;
(…這裡不加文字,markdown有序清單編号會出錯…)
2. 連續搜尋到最差結果的處理
代碼中的 P_L用于标示連續搜尋到比目前更差結果的次數。代碼中,如果連續發生了LIMIT次(代碼中是10000次),就意味着目前解空間已經找不到更好的結果了。
多次執行代碼測試會發現,原文代碼更換城市次數卻僅僅在1w次左右,最終路徑在450附近,這就說明了原文代碼會陷入局部最優解。即使每次把P_L複位為0,執行次數是上去了,但是仍然會陷入局部最優解。
3. 必須以一定機率接受更差結果
能夠跳出全局最優,正是模拟退火算法最大的優勢所在。
原來,原文中計算的是exp(dE / t),注意dE是為正的,是以exp(dE / t)恒大于1. 是以更差的結果永遠不會被選中。而由于交換操作中每次隻交換curPath的一對城市,是以陷入了curPath的局部最優解。
修改為exp(-dE / t)後,curPath有一定機率被替換。測試結果發現,更換城市次數在700w次左右,最終路徑在370附近。顯然獲得了全局最優解。
其他不錯的資料,如随機模拟的基本思想和常用采樣方法(sampling)2.
- http://blog.csdn.net/xianlingmao/article/details/7768833. ↩
- http://blog.csdn.net/xianlingmao/article/details/7768833. ↩