1. 程式人生 > >模擬退火演算法解決TSP問題+Python實現

模擬退火演算法解決TSP問題+Python實現

網上看見的比喻:

爬山演算法:兔子朝著比現在高的地方跳去。它找到了不遠處的最高山峰。但是這座山不一定是珠穆朗瑪峰。這就是爬山演算法,它不能保證區域性最優值就是全域性最優值。

模擬退火:兔子喝醉了。它隨機地跳了很長時間。這期間,它可能走向高處,也可能踏入平地。但是,它漸漸清醒了並朝最高方向跳去。這就是模擬退火。

旅行商問題 ( TSP,Traveling Salesman Problem ) :有N個城市,要求從其中某個問題出發,唯一遍歷所有城市,再回到出發的城市,求最短的路線。

       旅行商問題屬於所謂的NP完全問題(是世界七大數學難題之一。NP的英文全稱是Non-deterministicPolynomial的問題,即多項式複雜程度的非確定性問題),精確的解決TSP只能通過窮舉所有的路徑組合,其時間複雜度是O(N!) 。

求解旅行商問題的已有演算法很多:

  • 精確演算法如線性規劃方法、動態規劃方法、分支定界方法;
  • 近似方法有插入法、最近鄰演算法、Clark&Wright演算法、生成樹法、Chrisrofides演算法、r-opt演算法、混合演算法、概率演算法等。
  • 近年來的一些嘗試:緊急搜尋方法、遺傳演算法、模擬退火演算法、神經網路方法、蟻群演算法等。

使用模擬退火演算法可以比較快的求出TSP的一條近似最優路徑。模擬退火解決TSP的思路:

  1. 產生一條新的遍歷路徑P(i+1),計算路徑P(i+1)的長度L( P(i+1))
  2. 若L(P(i+1))< L(P(i)),則接受P(i+1)為新的路徑,否則以模擬退火的那個概率接受P(i+1) ,然後降溫
  3. 重複步驟1,2直到滿足退出條件

產生新的遍歷路徑的方法有很多,下面列舉其中3種:

  1. 隨機選擇2個節點,交換路徑中的這2個節點的順序。
  2. 隨機選擇2個節點,將路徑中這2個節點間的節點順序逆轉。
  3. 隨機選擇3個節點m,n,k,然後將節點m與n間的節點移位到節點k後面。

python 3.6.1 程式碼如下

import numpy as np
import matplotlib.pyplot as plt 
import pdb

"旅行商問題 ( TSP , Traveling Salesman Problem )"
coordinates = np.array([[565.0,575.0],[25.0,185.0],[345.0,750.0],[945.0,685.0],[845.0,655.0],
                        [880.0,660.0],[25.0,230.0],[525.0,1000.0],[580.0,1175.0],[650.0,1130.0],
                        [1605.0,620.0],[1220.0,580.0],[1465.0,200.0],[1530.0,  5.0],[845.0,680.0],
                        [725.0,370.0],[145.0,665.0],[415.0,635.0],[510.0,875.0],[560.0,365.0],
                        [300.0,465.0],[520.0,585.0],[480.0,415.0],[835.0,625.0],[975.0,580.0],
                        [1215.0,245.0],[1320.0,315.0],[1250.0,400.0],[660.0,180.0],[410.0,250.0],
                        [420.0,555.0],[575.0,665.0],[1150.0,1160.0],[700.0,580.0],[685.0,595.0],
                        [685.0,610.0],[770.0,610.0],[795.0,645.0],[720.0,635.0],[760.0,650.0],
                        [475.0,960.0],[95.0,260.0],[875.0,920.0],[700.0,500.0],[555.0,815.0],
                        [830.0,485.0],[1170.0, 65.0],[830.0,610.0],[605.0,625.0],[595.0,360.0],
                        [1340.0,725.0],[1740.0,245.0]])

#得到距離矩陣的函式
def getdistmat(coordinates):
    num = coordinates.shape[0] #52個座標點
    distmat = np.zeros((52,52)) #52X52距離矩陣
    for i in range(num):
        for j in range(i,num):
            distmat[i][j] = distmat[j][i]=np.linalg.norm(coordinates[i]-coordinates[j])
    return distmat

def initpara():
    alpha = 0.99
    t = (1,100)
    markovlen = 1000

    return alpha,t,markovlen
num = coordinates.shape[0]
distmat = getdistmat(coordinates) #得到距離矩陣


solutionnew = np.arange(num)
#valuenew = np.max(num)

solutioncurrent = solutionnew.copy()
valuecurrent =99000  #np.max這樣的原始碼可能同樣是因為版本問題被當做函式不能正確使用,應取一個較大值作為初始值
#print(valuecurrent)

solutionbest = solutionnew.copy()
valuebest = 99000 #np.max

alpha,t2,markovlen = initpara()
t = t2[1]

result = [] #記錄迭代過程中的最優解
while t > t2[0]:
    for i in np.arange(markovlen):

        #下面的兩交換和三角換是兩種擾動方式,用於產生新解
        if np.random.rand() > 0.5:# 交換路徑中的這2個節點的順序
            # np.random.rand()產生[0, 1)區間的均勻隨機數
            while True:#產生兩個不同的隨機數
                loc1 = np.int(np.ceil(np.random.rand()*(num-1)))
                loc2 = np.int(np.ceil(np.random.rand()*(num-1)))
                ## print(loc1,loc2)
                if loc1 != loc2:
                    break
            solutionnew[loc1],solutionnew[loc2] = solutionnew[loc2],solutionnew[loc1]
        else: #三交換
            while True:
                loc1 = np.int(np.ceil(np.random.rand()*(num-1)))
                loc2 = np.int(np.ceil(np.random.rand()*(num-1))) 
                loc3 = np.int(np.ceil(np.random.rand()*(num-1)))

                if((loc1 != loc2)&(loc2 != loc3)&(loc1 != loc3)):
                    break

            # 下面的三個判斷語句使得loc1<loc2<loc3
            if loc1 > loc2:
                loc1,loc2 = loc2,loc1
            if loc2 > loc3:
                loc2,loc3 = loc3,loc2
            if loc1 > loc2:
                loc1,loc2 = loc2,loc1

            #下面的三行程式碼將[loc1,loc2)區間的資料插入到loc3之後
            tmplist = solutionnew[loc1:loc2].copy()
            solutionnew[loc1:loc3-loc2+1+loc1] = solutionnew[loc2:loc3+1].copy()
            solutionnew[loc3-loc2+1+loc1:loc3+1] = tmplist.copy()  

        valuenew = 0
        for i in range(num-1):
            valuenew += distmat[solutionnew[i]][solutionnew[i+1]]
        valuenew += distmat[solutionnew[0]][solutionnew[51]]
       # print (valuenew)
        if valuenew<valuecurrent: #接受該解
           
            #更新solutioncurrent 和solutionbest
            valuecurrent = valuenew
            solutioncurrent = solutionnew.copy()

            if valuenew < valuebest:
                valuebest = valuenew
                solutionbest = solutionnew.copy()
        else:#按一定的概率接受該解
            if np.random.rand() < np.exp(-(valuenew-valuecurrent)/t):
                valuecurrent = valuenew
                solutioncurrent = solutionnew.copy()
            else:
                solutionnew = solutioncurrent.copy()
    t = alpha*t
    result.append(valuebest)
    print (t) #程式執行時間較長,列印t來監視程式進展速度
#用來顯示結果
plt.plot(np.array(result))
plt.ylabel("bestvalue")
plt.xlabel("t")
plt.show()

實驗結果: