在算法里面有一類問(wèn)題叫做最優(yōu)化問(wèn)題欠肾，這類問(wèn)題通常的特點(diǎn)是問(wèn)題的解空間很大忱详，用傳統(tǒng)的算法沒有辦法窮舉虚汛。
比如時(shí)間復(fù)雜度是O（NＸ依恕）的問(wèn)題，這類問(wèn)題解空間會(huì)非常非常的大卷哩，比如在N=15的時(shí)候蛋辈，就有1307674368000種可能性。

在用自己的筆記本上面計(jì)算了一下将谊，求解15個(gè)元素的全排列冷溶，跑了沒多久電腦就卡死了，然后在實(shí)驗(yàn)室的服務(wù)器上面試了一下尊浓，計(jì)算的用時(shí)大概是300s左右逞频，內(nèi)存占了60G左右。這已經(jīng)是一個(gè)很大的開銷了眠砾。

#如果你感興趣的話虏劲，可以在本地跑一下，改變‘a(chǎn)bcd’的大小褒颈，超過(guò)10以上，程序運(yùn)行時(shí)間開始明顯的上升
import itertools
for i in itertools.permutations('abcd',4):
        print (''.join(i))

對(duì)于這種復(fù)雜的組合優(yōu)化問(wèn)題励堡，解空間往往非常大谷丸，如果暴力計(jì)算的話，往往是需要計(jì)算幾十甚至是幾百個(gè)元素的全排列应结。剛剛在計(jì)算15的全排列的時(shí)候刨疼，已經(jīng)開銷已經(jīng)非常大了泉唁。對(duì)于，幾十幾百的全排列揩慕，是一個(gè)不可思議的計(jì)算量亭畜。

so，對(duì)于這類問(wèn)題迎卤，一個(gè)常見的解法是利用組合優(yōu)化的算法拴鸵，比如爬山法，模擬退火算法蜗搔，蟻群算法等等劲藐。這類算法的核心思想是：放棄遍歷所有解的嘗試，而是通過(guò)遍歷一部分解空間樟凄，來(lái)找到一個(gè)接近全局最優(yōu)解的一個(gè)近似最優(yōu)解聘芜。

模擬退火算法

模擬退火算法是一種啟發(fā)式尋優(yōu)算法，其出發(fā)點(diǎn)是基于物理中固體物質(zhì)的退火過(guò)程與一般組合優(yōu)化問(wèn)題之間的相似性汰现。模擬退火算法從某一較高初溫出發(fā)，伴隨溫度參數(shù)的不斷下降叔壤，結(jié)合概率突跳特性在解空間中隨機(jī)尋找目標(biāo)函數(shù)的全局最優(yōu)解瞎饲，即在局部最優(yōu)解能概率性地跳出并最終趨于全局最優(yōu)。

模擬退火算法是通過(guò)賦予搜索過(guò)程一種時(shí)變且最終趨于零的概率突跳性百新，從而可有效避免陷入局部極小并最終趨于全局最優(yōu)的串行結(jié)構(gòu)的優(yōu)化算法

image.png

如果我們現(xiàn)在有上面這樣一個(gè)函數(shù)圖企软，如圖二所示，現(xiàn)在想求函數(shù)的全局最優(yōu)解饭望。如果采用貪心策略仗哨，那么從A點(diǎn)開始試探，如果函數(shù)值繼續(xù)減少铅辞，那么試探過(guò)程就會(huì)繼續(xù)厌漂。而當(dāng)?shù)竭_(dá)點(diǎn)B時(shí)，顯然我們的探求過(guò)程就結(jié)束了斟珊。最終苇倡，只能找打一個(gè)局部最后解B。

模擬退火算法的搜索過(guò)程引入了隨機(jī)因素囤踩。模擬退火算法以一定的概率來(lái)接受一個(gè)比當(dāng)前解要差的解旨椒，因此有可能會(huì)跳出這個(gè)局部的最優(yōu)解，達(dá)到全局的最優(yōu)解堵漱。以上圖為例综慎，模擬退火算法在搜索到局部最優(yōu)解B后，會(huì)以一定的概率接受向右繼續(xù)移動(dòng)勤庐。也許經(jīng)過(guò)幾次這樣的不是局部最優(yōu)的移動(dòng)后會(huì)到達(dá)B 和C之間的峰點(diǎn)示惊，于是就跳出了局部最小值B好港。

根據(jù)Metropolis準(zhǔn)則，粒子在溫度T時(shí)趨于平衡的概率為exp(-ΔE/(kT))米罚，其中E為溫度T時(shí)的內(nèi)能钧汹，ΔE為其改變數(shù)，k為Boltzmann常數(shù)录择。Metropolis準(zhǔn)則常表示為

image.png

Metropolis準(zhǔn)則表明拔莱，在溫度為T時(shí)，出現(xiàn)能量差為dE的降溫的概率為P(dE)糊肠，表示為：P(dE) = exp( dE/(kT) )辨宠。其中k是一個(gè)常數(shù)，exp表示自然指數(shù)货裹，且dE<0嗤形。所以P和T正相關(guān)。這條公式就表示：溫度越高弧圆，出現(xiàn)一次能量差為dE的降溫的概率就越大赋兵；溫度越低，則出現(xiàn)降溫的概率就越小搔预。又由于dE總是小于0（因?yàn)橥嘶鸬倪^(guò)程是溫度逐漸下降的過(guò)程）霹期，因此dE/kT < 0 ，所以P(dE)的函數(shù)取值范圍是(0拯田，1) 历造。隨著溫度T的降低，P(dE)會(huì)逐漸降低船庇。

我們將一次向較差解的移動(dòng)看做一次溫度跳變過(guò)程吭产，我們以概率P(dE)來(lái)接受這樣的移動(dòng)。也就是說(shuō)鸭轮，在用固體退火模擬組合優(yōu)化問(wèn)題臣淤，將內(nèi)能E模擬為目標(biāo)函數(shù)值 f，溫度T演化成控制參數(shù) t窃爷，即得到解組合優(yōu)化問(wèn)題的模擬退火演算法：由初始解 i 和控制參數(shù)初值 t 開始邑蒋，對(duì)當(dāng)前解重復(fù)“產(chǎn)生新解→計(jì)算目標(biāo)函數(shù)差→接受或丟棄”的迭代，并逐步衰減 t 值按厘，算法終止時(shí)的當(dāng)前解即為所得近似最優(yōu)解医吊，這是基于蒙特卡羅迭代求解法的一種啟發(fā)式隨機(jī)搜索過(guò)程。

用模擬算法求解旅行商問(wèn)題

有 N ( <=20 ) 臺(tái) PC 放在機(jī)房?jī)?nèi)逮京，現(xiàn)在要求由你選定一臺(tái) PC遮咖，用共 N-1 條網(wǎng)線從這臺(tái)機(jī)器開始一臺(tái)接一臺(tái)地依次連接他們，最后接到哪個(gè)以及連接的順序也是由你選定的造虏，為了節(jié)省材料御吞，網(wǎng)線都拉直。求最少需要一次性購(gòu)買多長(zhǎng)的網(wǎng)線漓藕。

"""
模擬退火算法
"""
import numpy as np
import math
import matplotlib.pyplot as plt
import random
import copy

#calu_city_length 計(jì)算的是城市0到城市1陶珠，城市1到城市2，... ,城市18到城市19享钞，城市19到城市0揍诽，組成一個(gè)環(huán)的距離
#為什么要這么連接？因?yàn)檫@里做了一個(gè)coding上面的調(diào)整栗竖，假設(shè)城市之間的連接關(guān)系是不能調(diào)整的暑脆，只能按照city0到city1...city18到city19，city19到city0的關(guān)系連接狐肢。但是添吗，在后面調(diào)參的時(shí)候，需要調(diào)整的是城市之間的連接順序
def calu_city_length(city):
    sum_distance=0
    for i in range(len(city)-1):
        x_diff=city[i][0][0]-city[i+1][0][0]
        y_diff=city[i][0][1]-city[i+1][0][1]
        distance= math.sqrt(x_diff*x_diff+y_diff*y_diff)
        sum_distance += distance
    #計(jì)算city[n-1]-city[0]的距離
    distance=math.sqrt( (city[len(city)-1][0][0]-city[0][0][0])**2 + (city[len(city)-1][0][1]-city[0][0][1])**2 )
    sum_distance+=distance
    return sum_distance

def show(city):
    x_list=[]
    y_list=[]
    for i in range(len(city)):
        x=city[i][0][0]
        y=city[i][0][1]
        x_list.append(x)
        y_list.append(y)
    plt.plot(x_list,y_list)
    plt.show()


def perturb(city):
    a=range(len(city))
    c1,c2=random.sample(a,2)
    temp_x=city[c1][0][0]
    temp_y=city[c1][0][1]
    city[c1][0][0]=city[c2][0][0]
    city[c1][0][1]=city[c2][0][1]   
    city[c2][0][0]=temp_x
    city[c2][0][1]=temp_y
    return city


def sa(city,temperature):
    cost=[]
    while temperature >0.001:
        for i in range(100):
            len1=calu_city_length(city)
            #擾動(dòng)
            orignal_city=copy.deepcopy(city)
            #perturb 函數(shù)里面把city 已經(jīng)修改了份名，結(jié)果是city與perturb_city是相等的
            perturb_city=perturb(city)
            len2=calu_city_length(perturb_city)

            delta=len2-len1
            if (delta <0):   #新路線好與舊路線
                city=perturb_city
            else:
                #以一定的概率接收較差的解
                r=random.random()
                if math.exp(((-delta)/temperature)) > r:
                    city=perturb_city
                else:
                    city=orignal_city
        temperature=float(temperature) * 0.99
        cost.append(calu_city_length(city))
    return city,cost            


if __name__=="__main__":
    #城市的個(gè)數(shù)
    n=20
    temperature=2000

    #隨機(jī)初始化城市的坐標(biāo)
    #所有的城市都是連通的碟联，城市之間的距離就是坐標(biāo)的距離
    city=np.random.randint(1,100,size=(n,1,2))

    #訪問(wèn)第20城市的坐標(biāo)
    #坐標(biāo)x：city[19][0][0]
    #坐標(biāo)y：city[19][0][1]

    city,cost=sa(city,temperature)
    plt.plot(cost)
    plt.show()
    print (calu_city_length(city))  

運(yùn)行結(jié)果：

如果把每一次的代價(jià)輸出的話，那么結(jié)果如下僵腺，橫坐標(biāo)的迭代的次數(shù)鲤孵，縱坐標(biāo)是將所有城市連通的總代價(jià)。cost很明顯的隨著程序的運(yùn)行的次數(shù)在下降辰如，然后收斂在某個(gè)值普监。

2017-09-24 20-53-21屏幕截圖.png

PS:
在今年的數(shù)學(xué)建模里面，我們team就用的是這個(gè)算法琉兜，但是用之前一直不是很確定模擬退化的算法會(huì)不會(huì)在模型的算法里面生效凯正。直到程序運(yùn)行的結(jié)果出來(lái)，算法成功收斂！！淮菠！開心到爆炸Ｆ优馈！沥寥！

參考資料：
https://zhuanlan.zhihu.com/p/23968011
http://blog.csdn.net/kai940325/article/details/43267917