Python退火算法在高次方程的應(yīng)用
一,簡介
退火算法不言而喻,就是鋼鐵在淬煉過程中失溫而成穩(wěn)定態(tài)時(shí)的過程,熱力學(xué)上溫度(內(nèi)能)越高原子態(tài)越不穩(wěn)定,而溫度有一個(gè)向低溫區(qū)輻射降溫的物理過程,當(dāng)物質(zhì)內(nèi)能不再降低時(shí)候該物質(zhì)原子態(tài)逐漸成為穩(wěn)定有序態(tài),這對我們從隨機(jī)復(fù)雜問題中找出最優(yōu)解有一定借鑒意義,將這個(gè)過程化為算法,具體參見其他資料。
二,計(jì)算方程
我們所要計(jì)算的方程是f(x) = (x - 2) * (x + 3) * (x + 8) * (x - 9),是一個(gè)一元四次方程,我們稱為高次方程,當(dāng)然這個(gè)函數(shù)的開口是向上的,那么在一個(gè)無限長的區(qū)間內(nèi)我們可能找不出最大值點(diǎn),因此我們嘗試在較短區(qū)間內(nèi)解最小值點(diǎn),我們成為最優(yōu)解。
解法1:
毫無疑問,數(shù)學(xué)方法多次求導(dǎo)基本可以解出,但是這個(gè)過程較復(fù)雜,還容易算錯(cuò),我就不贅述了,讀者有時(shí)間自己可以嘗試解一下。
解法二:
這個(gè)解法就是暴力解決了,我們這里只求解區(qū)間[-10,10]上的最優(yōu)解,直接隨機(jī)200個(gè)點(diǎn),再除以10(這樣可以得到非整數(shù)橫坐標(biāo)),再依此計(jì)算其縱坐標(biāo)f(x),min{f(x)}一下,用list的index方法找出最小值對應(yīng)位置就行了,然后畫出圖形大致瞄一瞄。
直接貼代碼:
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("橫坐標(biāo)為:",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("縱坐標(biāo)為:",list_y)
#經(jīng)驗(yàn)證,這里算出來的結(jié)果6.5和最優(yōu)解1549都是對的
print("最小值為:",min(list_y))
num = min(list_y)
print("最優(yōu)解:",list_y.index(num)/10)
print("第",list_y.index(num)/10-10,"個(gè)位置取得最小值")
plt.plot(list_x, list_y, label='NM')
#plt.plot(x2, y2, label='Second Line')
plt.xlabel('X') #橫坐標(biāo)標(biāo)題
plt.ylabel('Y') #縱坐標(biāo)標(biāo)題
#plt.title('Interesting Graph\nCheck it out',loc="right") #圖像標(biāo)題
#plt.title('Interesting Graph\nCheck it out')
plt.legend() #顯示Fisrt Line和Second Line(label)的設(shè)置
plt.savefig('C:/Users/zhengyong/Desktop/1.png')
plt.show()
得到如下結(jié)果:
那么我們得出最優(yōu)解的坐標(biāo)是(6.5,-1549.6875),結(jié)果先放這里,接下來用退火算法看能不能解出。
解法三:
我們看一張圖(解法二中的方法得出的圖),然后講講退火算法的最核心的思想。
首先,先隨機(jī)一個(gè)[-10.10]之間的隨機(jī)解,作為初始解空間,比方說隨機(jī)了一個(gè)位于[-2.5.2.5]中最高的那個(gè)點(diǎn)就是點(diǎn)1(橫坐標(biāo)為x1),他有對于的縱坐標(biāo)的值y1,這時(shí)候我們把這個(gè)點(diǎn)的橫坐標(biāo)隨機(jī)加或者減去一個(gè)值(注意這個(gè)值的大小很重要,我們先叫他隨機(jī)移動(dòng)值),加或者減后得到新的橫坐標(biāo)的值x2,再算出這個(gè)橫坐標(biāo)的對應(yīng)縱坐標(biāo)(y2),對比之前的縱坐標(biāo)的大小,這里設(shè)置
delta = y2-y1,發(fā)現(xiàn)無論怎樣都是小于原先的縱坐標(biāo)(前提是隨機(jī)移動(dòng)值足夠?。@時(shí)候我們把新得到的x2賦值給x1,這時(shí)候現(xiàn)在的x2的值傳給x1,x1是原先隨機(jī)的值,這個(gè)過程可以重復(fù)iter_num 次,大小就根據(jù)自己的區(qū)間來。
上述的整個(gè)過程是在一個(gè)溫度下進(jìn)行的,這個(gè)過程結(jié)束后我們用溫度更新公式再次的更新溫度,再去重復(fù)上述步驟。
溫度更新我是用的常用的公式是T(t)=aT0(t-1),其中0.85≦a≦0.99。也可用相應(yīng)的熱能衰減公式來計(jì)算,T(t)=T0/(1+lnt),t=1,2,3,...,這都是簡單的狀態(tài)更新方法。
也就是說,不管你隨機(jī)的是幾我都能朝著優(yōu)化的方向前進(jìn)(前提是非最優(yōu)點(diǎn))。
其次,點(diǎn)2 是同理的,區(qū)別在于他是局部最優(yōu)解,那么跳出這個(gè)局部最優(yōu)解的機(jī)制是什么呢?
若初始點(diǎn)是(x3,y3),然后用上述方法得出(x4,y4),在點(diǎn)二處得到的delta肯定是大于0的,那么怎么辦呢?當(dāng)大于0的時(shí)候我們每次都有一定的概率來接受這個(gè)看起來不是最優(yōu)的點(diǎn),叫Metropolis準(zhǔn)則,具體是這樣的:
這里的E就是y,T就是當(dāng)前溫度,delta小于0就是百分百接受新值,否者就是按照這個(gè)概率接受,當(dāng)?shù)啻蔚臅r(shí)候,每次向右移動(dòng)的步長累加到點(diǎn)1 時(shí)候他就有可能找到最終的最優(yōu)解了,步長是累加的但是概率是累成的,意味著這個(gè)概率很小,但是一旦迭代次數(shù)多久一定會(huì)跑出來到最優(yōu)解處。
最優(yōu),點(diǎn)3不解釋了哈,和上面一樣。
那么我們上代碼:
#自己改寫的退火算法計(jì)算方程(x - 2) * (x + 3) * (x + 8) * (x - 9)的計(jì)算方法
#class沒啥用
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import pyplot as plt
#設(shè)置基本參數(shù)
#T初始溫度,T_stop,iter_num每個(gè)溫度的迭代次數(shù),Q溫度衰減次數(shù)
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)
#這個(gè)區(qū)間(2 * np.random.rand() - 1)本身是(-1,1),所以加上就是一個(gè)隨機(jī)加或者減過程
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("當(dāng)前溫度:",T_start)
T_start = T_start * Q
print("最優(yōu)解x是:", init_x) #這里最初寫的是new_x,所以結(jié)果一直不對
print("最優(yōu)解是:", init_y)
#比如我加上new_x,真假之間的誤差實(shí)際就是最后一次的賦值“init_x = new_x”
print("假最優(yōu)解x是:", new_x) #這里最初寫的是new_x,所以結(jié)果一直不對
print("假最優(yōu)解是:", 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') #橫坐標(biāo)標(biāo)題
plt.ylabel('Y for tuihuo') #縱坐標(biāo)標(biāo)題
#plt.title('Interesting Graph\nCheck it out',loc="right") #圖像標(biāo)題
#plt.title('Interesting Graph\nCheck it out')
plt.legend() #顯示Fisrt Line和Second Line(label)的設(shè)置
plt.savefig('C:/Users/zhengyong/Desktop/1.png')
plt.show()
這里用了class,發(fā)現(xiàn)并不需要,但是不想改了,就這樣吧。
最優(yōu)結(jié)果為:
得出的示意圖為:
三,總結(jié)
退火算法的具體思想我沒怎么講,但是核心的點(diǎn)我都寫出來了,經(jīng)過驗(yàn)證發(fā)現(xiàn)退火算法得出了(6.551677228904226,-1548.933671426107)的最優(yōu)解,看看解法二的(6.5,-1549.6875),我們發(fā)現(xiàn),呵呵,差不多,誤差來講的話,能接受,當(dāng)然讀者也可以多跑幾個(gè)數(shù)據(jù)出來驗(yàn)證。
我的實(shí)驗(yàn)環(huán)境是Python3.6,Numpy1.14.3,matplotlib2.2.2,64位win10,1709教育版,OS內(nèi)核16299.547,就這樣吧,盡量講詳細(xì)點(diǎn)。
總結(jié)
以上所述是小編給大家介紹的Python退火算法在高次方程的應(yīng)用,希望對大家有所幫助,如果大家有任何疑問請給我留言,小編會(huì)及時(shí)回復(fù)大家的。在此也非常感謝大家對腳本之家網(wǎng)站的支持!
相關(guān)文章
python將數(shù)據(jù)插入數(shù)據(jù)庫的代碼分享
在本篇文章里小編給大家整理的是關(guān)于python將數(shù)據(jù)插入數(shù)據(jù)庫的代碼內(nèi)容,有興趣的朋友們可以參考下。2020-08-08
python+matplotlib實(shí)現(xiàn)禮盒柱狀圖實(shí)例代碼
這篇文章主要介紹了python+matplotlib實(shí)現(xiàn)禮盒柱狀圖實(shí)例代碼,具有一定借鑒價(jià)值,需要的朋友可以參考下2018-01-01
使用SQLAlchemy操作數(shù)據(jù)庫表過程解析
這篇文章主要介紹了使用SQLAlchemy操作數(shù)據(jù)庫表過程解析,文中通過示例代碼介紹的非常詳細(xì),對大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友可以參考下2020-06-06
使用Pandas實(shí)現(xiàn)MySQL窗口函數(shù)的解決方法
本文主要介紹 MySQL 中的窗口函數(shù)row_number()、lead()/lag()、rank()/dense_rank()、first_value()、count()、sum()如何使用pandas實(shí)現(xiàn),同時(shí)二者又有什么區(qū)別,感興趣的朋友一起看看吧2023-02-02
python playwright 自動(dòng)等待和斷言詳解
這篇文章主要為大家介紹了python playwright 自動(dòng)等待和斷言,具有一定的參考價(jià)值,感興趣的小伙伴們可以參考一下,希望能夠給你帶來幫助2021-11-11
Python爬蟲實(shí)現(xiàn)爬取京東手機(jī)頁面的圖片(實(shí)例代碼)
下面小編就為大家分享一篇Python爬蟲實(shí)現(xiàn)爬取京東手機(jī)頁面的圖片實(shí)例,具有很好的參考價(jià)值,希望對大家有所幫助。一起跟隨小編過來看看吧2017-11-11
Python greenlet和gevent使用代碼示例解析
這篇文章主要介紹了Python greenlet和gevent使用代碼示例解析,文中通過示例代碼介紹的非常詳細(xì),對大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值,需要的朋友可以參考下2020-04-04

