1. 程式人生 > >模擬退火演算法解旅行商(TSP)問題

模擬退火演算法解旅行商(TSP)問題

該帖子的程式碼主要轉自模擬退火演算法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 = 0; i < n; i++) for (int j = i + 1; 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 = 0; i < n; i++) scanf("%lf %lf", &p[i].x, &p[i].y); } void Init(int n) { nCase = 0; path.len = 0; for (int i = 0; i < n; i++) { path.citys[i] = i; if (i != n - 1) { printf("%d--->", i); path.len += w[i][i + 1]; } 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 = 0; i < n; i++) { if (i != n - 1) 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 = 0, y = 0; while (x == y) { x = (int)(n * (rand() / (RAND_MAX + 1.0))); y = (int)(n * (rand() / (RAND_MAX + 1.0))); } Path ans = p; swap(ans.citys[x], ans.citys[y]); ans.len = 0; for (int i = 0; i < n - 1; i++) ans.len += w[ans.citys[i]][ans.citys[i + 1]]; //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 = 0; // 連續找到更差結果的次數 int P_F = 0; // while迴圈次數 while (1) //外迴圈,主要更新引數t,模擬退火過程 { for (int i = 0; i < ILOOP; i++) //內迴圈,尋找在一定溫度下的最優值 { newPath = GetNext(curPath, n); double dE = newPath.len - curPath.len; if (dE < 0) //如果找到更優值,直接更新 { curPath = newPath; P_L = 0; P_F = 0; } else { double rd = rand() / (RAND_MAX + 1.0); // Modify by DD: dE取負數才有可能接受更差解,否則e>1 double e = exp(-dE / t); if ( e > rd && e < 1) //如果找到比當前更差的解,以一定概率接受該解,並且這個概率會越來越小 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 0; }

資料檔案如下:

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

修改過的地方,程式碼中寫了註釋。
其中幾個重要的、需要解釋的地方有如下幾個:

  1. 記錄全域性最優解

    記錄最優解時,原文使用的是

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.