算法介紹
模拟退火算法（Simulate Anneal，SA）是一種通用機率演算法，用來在一個大的搜尋空間内找尋命題的最優解。美國實體學家 N.Metropolis 和同仁在1953年發表研究複雜系統、計算其中能量分布的文章，他們使用蒙特卡羅模拟法計算多分子系統中分子的能量分布。
模拟退火的出發點是基于實體中固體物質的退火過程與一般組合優化問題之間的相似性。模拟退火算法是一種通用的優化算法，其實體退火過程由加溫過程、等溫過程、冷卻過程這三部分組成。
啟發來源
在熱力學上，退火（annealing）現象指物體逐漸降溫的實體現象，溫度愈低，物體的能量狀态會低；夠低後，液體開始冷凝與結晶，在結晶狀态時，系統的能量狀态最低。大自然在緩慢降溫（亦即，退火）時，可“找到”最低能量狀态：結晶。但是，如果過程過急過快，快速降溫（亦稱「淬煉」，quenching）時，會導緻不是最低能态的非晶形。
如下圖所示，首先（左圖）物體處于非晶體狀态。我們将固體加溫至充分高（中圖），再讓其徐徐冷卻，也就退火（右圖）。加溫時，固體内部粒子随溫升變為無序狀，内能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達到平衡态，最後在常溫時達到基态，内能減為最小（此時物體以晶體形态呈現）。
似乎，大自然知道慢工出細活：緩緩降溫，使得物體分子在每一溫度時，能夠有足夠時間找到安頓位置，則逐漸地，到最後可得到最低能态，系統最安穩。
爬山算法
爬山算法是一種Greegy算法，該算法每次從目前解的臨近解空間中選擇一個最優解作為目前解，直到達到一個局部最優解。
如上圖：從A開始疊代搜尋解空間，尋找最優解，但是該算法采用的是貪心政策，是以當搜尋到B點時，發現再往後小幅度移動搜尋時，不能得到比B點更好的結果了，至此停止了搜尋，陷入到局部最優解的問題中。
模拟退火
介紹了以上的爬山算法，這種貪心算法往往得到的解并不是解空間中的全局最優解，如何有效的避免陷入局部最優解中，模拟退火算法應運而生。
其實模拟退火算法也是一種Greedy算法，但是在搜尋最優解過程中，加入了随機的因素，一定的機率接受一個比目前解要差的解，故此，存在一定的機率跳出局部最優解，進而找到真正的全局最優解。
- 算法流程
- 由一個産生函數從目前的解空間 { x 1 , x 2 , x 3 , . . . , x n } \{x_1, x_2, x_3, ..., x_n\} {x1,x2,x3,...,xn}中産生一個新解（一般采用随機方法生成新解，比如交換順序等）。
- 使用目标函數 f ( x ) f(x) f(x)，計算出新解對應的目标值 f ( x n e w ) f(x_{new}) f(xnew)。
-
判斷新解是否可以采納，變成目前的解。采納采用的準則——Metropolis準則：
p = { 1 , E ( x n e w ) < E ( x o l d ) e x p ( − E ( x n e w ) − E ( x o l d ) T ) , E ( x n e w ) ≥ E ( x o l d ) p = \begin{cases} 1, \quad \quad \quad \quad \quad \quad \quad \quad \quad E(x_{new}) < E(x_{old}) \\ exp(-\frac{E(x_{new}) - E(x_{old})}{T}), \quad E(x_{new}) \ge E(x_{old}) \end{cases} p={1,E(xnew)<E(xold)exp(−TE(xnew)−E(xold)),E(xnew)≥E(xold)
x n e w : 新 解 、 x o l d : 之 前 的 解 、 E ( x n e w ) : 當 前 解 的 能 量 、 T ： 當 前 溫 度 x_{new}: 新解、x_{old}: 之前的解、E(x_{new}): 目前解的能量、T：目前溫度 xnew:新解、xold:之前的解、E(xnew):目前解的能量、T：目前溫度根據以上的公式可以發現：随着溫度的降低， e x p ( − E ( x n e w ) − E ( x o l d ) T ) exp(-\frac{E(x_{new}) - E(x_{old})}{T}) exp(−TE(xnew)−E(xold))即 P P P采納比目前解差的新解的機率在不斷的降低。其中的溫度T，随着時間的流逝存在一個降溫的過程，降溫的速率可以自己定義，定義降溫速率 t _ a l p h a t\_alpha t_alpha。則：溫度的變化為： T = T ∗ t _ a l p h a T = T * t\_alpha T=T∗t_alpha
- 新解被采納後，使用新解替代目前解。
以上為一次算法的疊代。
- 算法參數
- 溫度初始值T。（影響模拟退火算法全局搜尋性能的重要因素之一、初始溫度高，則搜尋到全局最優解的可能性大，但是以要花費大量的計算時間；反之，則可節約計算時間，但全局搜尋性能可能受到影響。）
- 溫度降低速率，即退火速率 t _ a l p h a t\_alpha t_alpha。（同上，退火速率快，可節約計算時間，但全局搜尋性能可能受到影響。反之）
與普通的貪心算法對比(參考一個比較有趣的解釋)：
- 普通貪心算法：兔子朝着比現在低的地方跳去。它找到了不遠處的最低的山谷。但是這座山谷不一定最低的。
- 模拟退火：兔子喝醉了。它随機地跳了很長時間。這期間，它可能走向低處，也可能踏入平地。但是，它漸漸清醒了并朝最低的方向跳去。這就是模拟退火。
TSP-SA方法
TSP（旅行商問題）：一個商品推銷員要去若幹個城市推銷商品，該推銷員從一個城市出發，需要經過所有城市後，回到出發地。應如何選擇行進路線，以使總的行程最短。從圖論的角度來看，該問題實質是在一個帶權完全無向圖中，找一個權值最小的Hamilton回路。由于該問題的可行解是所有頂點的全排列，随着頂點數的增加，會産生組合爆炸，它是一個NP完全問題。
使用模拟退火算法求解TSP：
class Cities:
"""
storage the information:
- the location dataframe of cities
- the number of cities
- the distance matrix of cities
- the total distance of solution
"""
def __init__(self, coordinates):
self.coordinates = coordinates.copy()
self.city_number = coordinates.shape[0]
self.getDistance()
self.getSolutionDistance()
def getDistance(self):
"""
get the distance matrix of cities
"""
self.distances = np.zeros((self.city_number, self.city_number))
for i in range(self.city_number):
for j in range(i, self.city_number):
self.distances[i][j] = self.distances[j][i] = np.linalg.norm(
self.coordinates.iloc[i] - self.coordinates.iloc[j])
def getSolutionDistance(self):
"""
get the total distance of the array of cities by order
"""
self.solution_distance = 0
for i in range(self.city_number):
self.solution_distance += self.distances[i][(i + 1) % self.city_number]
return self.solution_distance
def getNewSolution(self):
"""
get the new sequence of cities array.
calculate the distance matrix and solution distance
"""
while True: # exchange the city index randomly
index_1 = np.random.randint(self.city_number - 1)
index_2 = np.random.randint(self.city_number - 1)
if index_1 != index_2:
if index_1 > index_2:
index_1, index_2 = index_2, index_1
break
city_index_1 = self.coordinates[index_1: index_1 + 1]
city_index_2 = self.coordinates[index_2: index_2 + 1]
above = self.coordinates[: index_1]
mid = self.coordinates[index_1 + 1: index_2]
below = self.coordinates[index_2 + 1:]
self.coordinates = above.append(city_index_2).append(mid).append(city_index_1).append(below)
self.city_number = self.coordinates.shape[0]
self.getDistance()
self.getSolutionDistance()
def SA(temperature, IN_Loop, temperature_min, best_path, t_alpha, result):
new_path = copy.deepcopy(best_path) # new solution
cur_path = copy.deepcopy(best_path) # current solution
while temperature > temperature_min:
for i in range(IN_Loop): # inner loop times
new_path.getNewSolution()
delta_distance = new_path.solution_distance - cur_path.solution_distance
if delta_distance < 0: # the better path
cur_path = copy.deepcopy(new_path)
else: # replace with random probability
if np.random.rand() < np.exp(- delta_distance / temperature):
cur_path = copy.deepcopy(new_path)
if best_path.solution_distance > cur_path.solution_distance:
# update the best solution
best_path = copy.deepcopy(cur_path)
result.append(best_path.solution_distance)
temperature *= t_alpha
return result
由于考慮到時間和運算，測試使用了15個城市，資料如下：
cities_location = [[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]]
cities = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O']
b = pd.DataFrame(cities_location, columns=['x', 'y'], index=cities)
city_path = Cities(b)
temperature = 100
IN_Loop = 1000
temperature_min = 1
alpha = 0.9
result = list()
result = pd.Series(SA(temperature, IN_Loop, temperature_min, city_path, alpha, result))
plt.plot(result)
plt.ylabel("best_distance")
plt.xlabel("iterators")
plt.show()
尋找最優解的過程如下：
參考資料
- 算法大全
- 模拟退火算法 - ranjiewen
- Simulated annealing - Wikipedia