NO IMAGE

為什麼我的眼裡常含淚水?因為我有一個演算法不會。為了節約點眼淚,今天我們就來介紹著名的模擬退火演算法(Simulated Annealing),它是一種基於蒙特卡洛思想設計的近似求解最優化問題的方法。

這是本系列文章的第二篇,我們通過一個例項來程式設計演示模擬退火的執行。特別地,我們這裡所採用的例項是著名的“旅行商問題”(TSP,Traveling Salesman Problem),它是哈密爾頓迴路的一個例項化問題,也是最早被提出的NP問題之一。

在你閱讀本文前,希望你已經讀過《模擬退火演算法系列之(一):通俗理解》一文,因為本文我們主要討論實現上的一些策略,對於那些基本的概念並不會過多涉及。

歡迎關注白馬負金羈的部落格 http://blog.csdn.net/baimafujinji,為保證公式、圖表得以正確顯示,強烈建議你從該地址上檢視原版博文。本部落格主要關注方向包括:數字影象處理、演算法設計與分析、資料結構、機器學習、資料探勘、統計分析方法、自然語言處理。


TSP是一個最常被用來解釋模擬退火用法的問題,因為這個問題比較有名,我們這裡不贅言重述,下面直接給出C 實現的程式碼:


#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 1000   //概率選擇上限
#define OLOOP 20    //外迴圈次數
#define ILOOP 100   //內迴圈次數
using namespace std;
//定義路線結構體
struct Path
{
int citys[N];
double len;
};
//定義城市點座標
struct Point
{
double x, y;
};
Path bestPath;        //記錄最優路徑
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;
bestPath.len = 0;
for(int i = 0; i < n; i  )
{
bestPath.citys[i] = i;
if(i != n - 1)
{
printf("%d--->", i);
bestPath.len  = w[i][i   1];
}
else
printf("%d\n", i);
}
printf("\nInit path length is : %.3lf\n", bestPath.len);
printf("-----------------------------------\n\n");
}
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);
printf("-----------------------------------\n\n");
}
Path GetNext(Path p, int n)
{
Path ans = p;
int x = (int)(n * (rand() / (RAND_MAX   1.0)));
int y = (int)(n * (rand() / (RAND_MAX   1.0)));
while(x == y)
{
x = (int)(n * (rand() / (RAND_MAX   1.0)));
y = (int)(n * (rand() / (RAND_MAX   1.0)));
}
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((unsigned)(time(NULL)));
Path curPath = bestPath;
Path newPath = bestPath;
int P_L = 0;
int P_F = 0;
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);
//如果找到比當前更差的解,以一定概率接受該解,並且這個概率會越來越小
if(exp(dE / t) > rd && exp(dE / t) < 1)
curPath = newPath;
P_L  ;
}
if(P_L > LIMIT)
{
P_F  ;
break;
}
}
if(curPath.len < bestPath.len)
bestPath = curPath;
if(P_F > OLOOP || t < EPS)
break;
t *= DELTA;
}
}
int main(int argc, const char * argv[]) {
freopen("TSP.data", "r", stdin);
int n;
Input(p, n);
GetDist(p, n);
Init(n);
SA(n);
Print(bestPath, n);
printf("Total test times is : %d\n", nCase);
return 0;
}

注意由於是基於蒙特卡洛的方法,所以上面程式碼每次得出的結果並不完全一致。你可以通過增加迭代的次數來獲得一個更優的結果。

我們這裡需要說明的是,在之前的文章裡,我們用求最小值的例子來解釋模擬退火的執行:如果新一輪的計算結果更前一輪之結果更小,那麼我們就接受它,否則就以一個概率來拒絕或接受它,而這個拒絕的概率會隨著溫度的降低(也即是迭代次數的增加)而變大(也就是接受的概率會越來越小)。

但現在我們面對一個TSP問題,我們如何定義或者說如何獲取下一輪將要被考察的哈密爾頓路徑呢?在一元函式最小值的例子中,下一輪就是指向左或者向右移動一小段距離。而在TSP問題中,我們可以採用的方式其實是很多的。上面程式碼中GetNext()函式所採用的方式是隨機交換兩個城市在路徑中的順序。例如當前路徑為 A->B->C->D->A,那麼下一次路徑就可能是A->D->C->B->A,即交換B和D。而在文獻【3】中,作者取樣的程式碼如下(我們擷取一個片段,完整程式碼請參考原文):

public class Tour{
... ...
// Creates a random individual
public void generateIndividual() {
// Loop through all our destination cities and add them to our tour
for (int cityIndex = 0; cityIndex < TourManager.numberOfCities(); cityIndex  ) {
setCity(cityIndex, TourManager.getCity(cityIndex));
}
// Randomly reorder the tour
Collections.shuffle(tour);
}
... ...
}

可見作者的方法是把上一輪路徑做一個隨機的重排(這顯然也是一種策略)。

TSP.data的資料格式如下,第一行的數字表示一個有多少座城市,第2至最後一行,每行有兩個數字表示,城市的座標(平面直角座標系)。例如:
6
20 80
16 84
23 66
62 90
11 9
35 28
最後請讀者自行編譯執行程式並觀察分析輸出的結果。


參考文獻與推薦閱讀材料

【1】關於哈密爾頓問題和TSP問題請參考下面兩個資料以瞭解更多:

【2】上面的C 程式碼在下面這個兩個帖子中都有給出,原作者無法考證

【3】關於TSP問題的一個Java語言實現的原始碼,請參考