模擬退火算法(附代碼)

上一篇我們講了旅行商問題及遺傳算法來解決此類問題：遺傳算法(附代碼)

今天介紹另外一種解決此類NP問題的方法是模擬退火算法(Simulated Annealing, SA)

模擬退火算法的思想借鑑於固體的退火原理，當固體的溫度很高的時候，內能比較大，固體的內部粒子處於快速無序運動，當溫度慢慢降低的過程中，固體的內能減小，粒子的慢慢趨於有序，最終，當固體處於常溫時，內能達到最小，此時，粒子最為穩定。模擬退火算法便是基於這樣的原理設計而成。

那麼為什麼在算法開始的時候要處於不穩定的狀態呢？我們先看來一下爬山法吧

爬山法的問題

爬山法是一種貪婪算法，在當前位置的附近尋找一個更大的值，不斷迭代這個過程直到處於穩定狀態

如圖中要尋找這個函數的最大值，採用爬山法的話，如果初始點選擇D點，那麼爬山法能夠找到最大值B點，但是如果初始值選擇C或者E點，那麼爬山法找到的最大值就是局部最優A點，而不是全局最優點B點

這就是爬山法的局限性，如果選擇的初始點不好，那麼算法很有可能會陷入局部最優解，而模擬退火引入的溫度，在溫度越高的時候越有機會跳到一個相對不那麼優的解，從而跳出局部最優。

模擬退火的過程

1.初始化一個溫度T, 如T=100

2.老的解損失為: old_cost, 生成一個新的解,損失為: new_cost

3.計算當前是否採納新的解概率: P = math.exp(-(new_cost-old_cost)/T)

4.如果損失new_cost<old_cost或者隨機一個(0,1)之間數值<P, 則接受這個新解，否則還是使用老的解

5.當前的溫度T下降一定的值(溫度降低)

6.一直迭代2~5步，直到穩定狀態

可以看到一開始溫度T很高的時候P也很大，那麼採納新的解(損失比老的解更大)的概率也更高，從而跳出了局部最優

隨著溫度T的下降P也越來越小，會越來越多的採用損失大小來評價是否採納新的解

代碼實現

import randomimport mathimport matplotlib.pyplot as plt

class SA:    def __init__(self):        self.T = 10000.        self.cool = 0.95
    def new_route(self, tour, data):        c1 = random.randint(0, len(tour['tour']) - 1)        c2 = random.randint(0, len(tour['tour']) - 1)
        while c1 == c2:            c2 = random.randint(0, len(tour['tour']) - 1)
        new_tour = []        for i in range(len(tour['tour'])):            if i == c1:                new_tour.append(tour['tour'][c2])            elif i == c2:                new_tour.append(tour['tour'][c1])            else:                new_tour.append(tour['tour'][i])
        return self.get_tour_detail(new_tour, data)
    def get_distance(self, c1, c2):                xd = abs(c1['x'] - c2['x'])        yd = abs(c1['y'] - c2['y'])        distance = math.sqrt(xd * xd + yd * yd)
        return distance
    def get_cost(self, distance):        return distance
    def get_tour(self, data):                tour = []        for key, value in data.items():            tour.append(key)
        random.shuffle(tour)
        return self.get_tour_detail(tour, data)
    def get_tour_detail(self, tour, data):        tmp = None        distance = 0        for item in tour:            if tmp is not None:                distance_tmp = self.get_distance(data[tmp], data[item])                distance += distance_tmp
            tmp = item
        return {'tour': tour, 'distance': distance, 'cost': self.get_cost(distance)}
    def run(self, data):        route = self.get_tour(data)        print 'before route:'        print route        i = 0        while self.T > 0.1:            new_route = self.new_route(route, data)            old_cost, new_cost = route['cost'], new_route['cost']
            if new_cost < old_cost or random.random() < math.exp(-(new_cost - old_cost) / self.T):                route = new_route
            self.T = self.T * self.cool            i += 1
        print 'total gen:', i        print route        return route

def init_data(num):    data = {}    def get_random_point():                return {            'x': random.randint(1, 99),            'y': random.randint(1, 99)        }
    for i in range(num):        data[i + 1] = get_random_point()
    return data

def show(citys, tour):        plt.bar(left=0, height=100, width=100, color=(0, 0, 0, 0), edgecolor=(0, 0, 0, 0))    plt.title(u'tsp')    plt.xlabel('total distance: %s m' % tour['distance'])
    x = []    y = []    i = 0    for item in tour['tour']:        city = citys[item]        x.append(city['x'])        y.append(city['y'])
        i += 1        if i == 1:            plt.plot(city['x'], city['y'], 'or')        else:            plt.plot(city['x'], city['y'], 'bo')
    plt.plot(x, y, 'g')    plt.xlim(0.0, 100)    plt.ylim(0.0, 100)    plt.show()

if __name__ == '__main__':    scale, num, max_gen, pc, elite = 50, 10, 1000, 0.8, 0.2    data = init_data(num)    sa = SA()    new_fittest = sa.run(data)    show(data, new_fittest)
運行效果：

參考
[1].<<集體智慧編程>>
[2].https://www.cnblogs.com/heaad/archive/2010/12/20/1911614.html
[3].https://zhuanlan.zhihu.com/p/33184423

長按二維碼, 關注本公眾號