Python退火算法在高次方程的應用


一,簡介

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

二,計算方程

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

解法1:

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

解法二:

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

直接貼代碼:

 1 import random
 2 import matplotlib.pyplot as plt
 3 
 4 list_x = []
 5 # for i in range(1):
 6 #     #print(random.randint(0,100))
 7 #     for i in range(0,100):
 8 #         print("sss",i)
 9 #
10 #     list_x.append(random.randint(0,100))
11 for i in range(-100,100):
12     list_x.append(i/10)
13 
14 print("橫坐標為:",list_x)
15 print(len(list_x))
16 
17 
18 list_y = []
19 for x in list_x:
20     # print(x)
21     #y = x*x*x - 60*x*x -4*x +6
22     y = (x - 2) * (x + 3) * (x + 8) * (x - 9)
23     list_y.append(y)
24 print("縱坐標為:",list_y)
25 
26 #經驗證,這里算出來的結果6.5和最優解1549都是對的
27 print("最小值為:",min(list_y))
28 num = min(list_y)
29 print("最優解:",list_y.index(num)/10)
30 print("",list_y.index(num)/10-10,"個位置取得最小值")
31 
32 plt.plot(list_x, list_y, label='NM')
33 #plt.plot(x2, y2, label='Second Line')
34 plt.xlabel('X')  #橫坐標標題
35 plt.ylabel('Y')  #縱坐標標題
36 #plt.title('Interesting Graph\nCheck it out',loc="right")   #圖像標題
37 #plt.title('Interesting Graph\nCheck it out')
38 plt.legend()    #顯示Fisrt Line和Second Line(label)的設置
39 plt.savefig('C:/Users/zhengyong/Desktop/1.png')
40 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不解釋了哈,和上面一樣。

那么我們上代碼:

 1 #自己改寫的退火算法計算方程(x - 2) * (x + 3) * (x + 8) * (x - 9)的計算方法
 2 #class沒啥用
 3 import numpy as np
 4 import matplotlib.pyplot as plt
 5 from matplotlib import pyplot as plt
 6 
 7 
 8 #設置基本參數
 9 #T初始溫度,T_stop,iter_num每個溫度的迭代次數,Q溫度衰減次數
10 class Tuihuo_alg():
11     def __init__(self,T_start,iter_num,T_stop,Q,xx,init_x):
12         self.T_start = T_start
13         self.iter =iter_num
14         self.T_stop = T_stop
15         self.Q = Q
16         self.xx = xx
17         self.init_x = init_x
18     # def cal_x2y(self):
19     #     return (x - 2) * (x + 3) * (x + 8) * (x - 9)
20 
21 
22 if __name__ == '__main__':
23 
24     def cal_x2y(x):
25         #print((x - 2) * (x + 3) * (x + 8) * (x - 9))
26         return (x - 2) * (x + 3) * (x + 8) * (x - 9)
27     T_start = 1000
28     iter_num = 1000
29     T_stop = 1
30     Q = 0.95
31     K = 1
32     l_boundary = -10
33     r_boundary = 10
34     #初始值
35     xx = np.linspace(l_boundary, r_boundary, 300)
36     yy = cal_x2y(xx)
37     init_x =10 * ( 2 * np.random.rand() - 1)
38     print("init_x:",init_x)
39 
40     t = Tuihuo_alg(T_start,iter_num,T_stop,Q,xx,init_x)
41 
42     val_list = [init_x]
43     while T_start>T_stop:
44         for i in range(iter_num):
45             init_y = cal_x2y(init_x)
46             #這個區間(2 * np.random.rand() - 1)本身是(-1,1),所以加上就是一個隨機加或者減過程
47             new_x = init_x + (2 * np.random.rand() - 1)
48             if l_boundary <= new_x <= r_boundary:
49                 new_y = cal_x2y(new_x)
50         #print("new_x:",new_x)
51         #print('new_y:',new_y)
52                 delta = new_y - init_y  #新減舊
53                 if delta < 0:
54                     init_x = new_x
55                 else:
56                     p = np.exp(-delta / (K * T_start))
57                     if np.random.rand() < p:
58                         init_x = new_x
59             #print("new_x:",new_x)
60             #print("當前溫度:",T_start)
61         T_start = T_start * Q
62 
63 print("最優解x是:", init_x)   #這里最初寫的是new_x,所以結果一直不對
64 print("最優解是:", init_y)
65 #比如我加上new_x,真假之間的誤差實際就是最后一次的賦值“init_x = new_x”
66 print("假最優解x是:", new_x)   #這里最初寫的是new_x,所以結果一直不對
67 print("假最優解是:", new_y)
68 
69 xx = np.linspace(l_boundary,r_boundary,300)
70 yy = cal_x2y(xx)
71 plt.plot(xx, yy, label='Tuihuo')
72 #plt.plot(x2, y2, label='Second Line')
73 plt.xlabel('X for tuihuo')  #橫坐標標題
74 plt.ylabel('Y for tuihuo')  #縱坐標標題
75 #plt.title('Interesting Graph\nCheck it out',loc="right")   #圖像標題
76 #plt.title('Interesting Graph\nCheck it out')
77 plt.legend()    #顯示Fisrt Line和Second Line(label)的設置
78 plt.savefig('C:/Users/zhengyong/Desktop/1.png')
79 plt.show()

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

最優結果為:

得出的示意圖為:

 

三,總結

退火算法的具體思想我沒怎么講,但是核心的點我都寫出來了,經過驗證發現退火算法得出了(6.551677228904226,-1548.933671426107)的最優解,看看解法二的(6.5,-1549.6875),我們發現,呵呵,差不多,誤差來講的話,能接受,當然讀者也可以多跑幾個數據出來驗證。

我的實驗環境是Python3.6,Numpy1.14.3,matplotlib2.2.2,64位win10,1709教育版,OS內核16299.547,就這樣吧,盡量講詳細點。

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM