python進階學習：Python退火演算法在高次方程的應用

　　退火演算法就是鋼鐵在淬鍊過程中失溫而成穩定態時的過程，熱力學上溫度(內能)越高原子態越不穩定。這篇文章主要介紹了Python退火演算法在高次方程的應用,需要的朋友可以參考下。

　　一、簡介

　　退火演算法不言而喻，就是鋼鐵在淬鍊過程中失溫而成穩定態時的過程，熱力學上溫度(內能)越高原子態越不穩定，而溫度有一個向低溫區輻射降溫的物理過程，當物質內能不再降低時候該物質原子態逐漸成為穩定有序態，這對我們從隨機複雜問題中找出最優解有一定借鑑意義，將這個過程化為演算法，具體參見其他資料。

　　二、計算方程

　　我們所要計算的方程是f(x) = (x - 2) * (x + 3) * (x + 8) * (x - 9),是一個一元四次方程，我們稱為高次方程，當然這個函式的開口是向上的，那麼在一個無限長的區間內我們可能找不出最大值點，因此我們嘗試在較短區間內解最小值點，我們成為最優解。

　　解法一：

　　毫無疑問，數學方法多次求導基本可以解出，但是這個過程較複雜，還容易算錯，我就不贅述了，讀者有時間自己可以嘗試解一下。

　　解法二：

　　這個解法就是暴力解決了，我們這裡只求解區間[-10,10]上的最優解，直接隨機200個點，再除以10(這樣可以得到非整數橫座標)，再依此計算其縱座標f(x)，min{f(x)}一下，用list的index方法找出最小值對應位置就行了，然後畫出圖形大致瞄一瞄。

　　直接貼程式碼：

　　import random

　　import matplotlib.pyplot as plt

　　list_x = []

　　# for i in range(1):

　　# #print(random.randint(0,100))

　　# for i in range(0,100):

　　# print("sss",i)

　　#

　　# list_x.append(random.randint(0,100))

　　for i in range(-100,100):

　　list_x.append(i/10)

　　print("橫座標為：",list_x)

　　print(len(list_x))

　　list_y = []

　　for x in list_x:

　　# print(x)

　　#y = x*x*x - 60*x*x -4*x +6

　　y = (x - 2) * (x + 3) * (x + 8) * (x - 9)

　　list_y.append(y)

　　print("縱座標為：",list_y)

　　#經驗證，這裡算出來的結果6.5和最優解1549都是對的

　　print("最小值為：",min(list_y))

　　num = min(list_y)

　　print("最優解：",list_y.index(num)/10)

　　print("第",list_y.index(num)/10-10,"個位置取得最小值")

　　plt.plot(list_x, list_y, label='NM')

　　#plt.plot(x2, y2, label='Second Line')

　　plt.xlabel('X') #橫座標標題

　　plt.ylabel('Y') #縱座標標題

　　#plt.title('Interesting Graph\nCheck it out',loc="right") #影像標題

　　#plt.title('Interesting Graph\nCheck it out')

　　plt.legend() #顯示Fisrt Line和Second Line(label)的設定

　　plt.savefig('C:/Users/zhengyong/Desktop/1.png')

　　plt.show()

　　得到如下結果：

　　

　　那麼我們得出最優解的座標是(6.5，-1549.6875)，結果先放這裡，接下來用退火演算法看能不能解出。

　　解法三：

　　我們看一張圖(解法二中的方法得出的圖)，然後講講退火演算法的最核心的思想。

　　首先，先隨機一個[-10.10]之間的隨機解，作為初始解空間，比方說隨機了一個位於[-2.5.2.5]中最高的那個點就是點1(橫座標為x1)，他有對於的縱座標的值y1，這時候我們把這個點的橫座標隨機加或者減去一個值(注意這個值的大小很重要，我們先叫他隨機移動值)，加或者減後得到新的橫座標的值x2，再算出這個橫座標的對應縱座標(y2)，對比之前的縱座標的大小，這裡設定

　　delta = y2-y1，發現無論怎樣都是小於原先的縱座標(前提是隨機移動值足夠小)，這時候我們把新得到的x2賦值給x1，這時候現在的x2的值傳給x1，x1是原先隨機的值，這個過程可以重複iter_num 次，大小就根據自己的區間來。

　　上述的整個過程是在一個溫度下進行的，這個過程結束後我們用溫度更新公式再次的更新溫度，再去重複上述步驟。

　　溫度更新我是用的常用的公式是T(t)=aT0(t-1)，其中0.85≦a≦0.99。也可用相應的熱能衰減公式來計算，T(t)=T0/(1+lnt),t=1,2,3,…，這都是簡單的狀態更新方法。

　　也就是說，不管你隨機的是幾我都能朝著最佳化的方向前進(前提是非最優點)。

　　其次，點2 是同理的，區別在於他是區域性最優解，那麼跳出這個區域性最優解的機制是什麼呢?

　　若初始點是(x3,y3),然後用上述方法得出(x4,y4)，在點二處得到的delta肯定是大於0的，那麼怎麼辦呢?當大於0的時候我們每次都有一定的機率來接受這個看起來不是最優的點，叫Metropolis準則，具體是這樣的：

　　

　　這裡的E就是y，T就是當前溫度，delta小於0就是百分百接受新值，否者就是按照這個機率接受，當迭代多次的時候，每次向右移動的步長累加到點1 時候他就有可能找到最終的最優解了，步長是累加的但是機率是累成的，意味著這個機率很小，但是一旦迭代次數多久一定會跑出來到最優解處。

　　最優，點3不解釋了哈，和上面一樣。

　　那麼我們上程式碼：

　　#自己改寫的退火演算法計算方程(x - 2) * (x + 3) * (x + 8) * (x - 9)的計算方法

　　#class沒啥用

　　import numpy as np

　　import matplotlib.pyplot as plt

　　from matplotlib import pyplot as plt

　　#設定基本引數

　　#T初始溫度，T_stop，iter_num每個溫度的迭代次數，Q溫度衰減次數

　　class Tuihuo_alg():

　　def __init__(self,T_start,iter_num,T_stop,Q,xx,init_x):

　　self.T_start = T_start

　　self.iter =iter_num

　　self.T_stop = T_stop

　　self.Q = Q

　　self.xx = xx

　　self.init_x = init_x

　　# def cal_x2y(self):

　　# return (x - 2) * (x + 3) * (x + 8) * (x - 9)

　　if __name__ == '__main__':

　　def cal_x2y(x):

　　#print((x - 2) * (x + 3) * (x + 8) * (x - 9))

　　return (x - 2) * (x + 3) * (x + 8) * (x - 9)

　　T_start = 1000

　　iter_num = 1000

　　T_stop = 1

　　Q = 0.95

　　K = 1

　　l_boundary = -10

　　r_boundary = 10

　　#初始值 鄭州人流醫院

　　xx = np.linspace(l_boundary, r_boundary, 300)

　　yy = cal_x2y(xx)

　　init_x =10 * ( 2 * np.random.rand() - 1)

　　print("init_x:",init_x)

　　t = Tuihuo_alg(T_start,iter_num,T_stop,Q,xx,init_x)

　　val_list = [init_x]

　　while T_start>T_stop:

　　for i in range(iter_num):

　　init_y = cal_x2y(init_x)

　　#這個區間(2 * np.random.rand() - 1)本身是(-1,1)，所以加上就是一個隨機加或者減過程

　　new_x = init_x + (2 * np.random.rand() - 1)

　　if l_boundary <= new_x <= r_boundary:

　　new_y = cal_x2y(new_x)

　　#print("new_x:",new_x)

　　#print('new_y:',new_y)

　　delta = new_y - init_y #新減舊

　　if delta < 0:

　　init_x = new_x

　　else:

　　p = np.exp(-delta / (K * T_start))

　　if np.random.rand() < p:

　　init_x = new_x

　　#print("new_x:",new_x)

　　#print("當前溫度：",T_start)

　　T_start = T_start * Q

　　print("最優解x是：", init_x) #這裡最初寫的是new_x,所以結果一直不對

　　print("最優解是：", init_y)

　　#比如我加上new_x，真假之間的誤差實際就是最後一次的賦值“init_x = new_x”

　　print("假最優解x是：", new_x) #這裡最初寫的是new_x,所以結果一直不對

　　print("假最優解是：", new_y)

　　xx = np.linspace(l_boundary,r_boundary,300)

　　yy = cal_x2y(xx)

　　plt.plot(xx, yy, label='Tuihuo')

　　#plt.plot(x2, y2, label='Second Line')

　　plt.xlabel('X for tuihuo') #橫座標標題

　　plt.ylabel('Y for tuihuo') #縱座標標題

　　#plt.title('Interesting Graph\nCheck it out',loc="right") #影像標題

　　#plt.title('Interesting Graph\nCheck it out')

　　plt.legend() #顯示Fisrt Line和Second Line(label)的設定

　　plt.savefig('C:/Users/zhengyong/Desktop/1.png')

　　plt.show()

　　這裡用了class，發現並不需要，但是不想改了，就這樣吧。

　　最優結果為：

　　

　　得出示意圖