Python實現二分法和黃金分割法


  運籌學課上,首先介紹了非線性規划算法中的無約束規划算法。二分法和黃金分割法是屬於無約束規划算法的一維搜索法中的代表。

  二分法:$$x_{1}^{(k+1)}=\frac{1}{2}(x_{R}^{(k)}+x_{L}^{(k)}-\Delta)$$$$x_{2}^{(k+1)}=\frac{1}{2}(x_{R}^{(k)}+x_{L}^{(k)}+\Delta)$$

  黃金分割法:$$x_{1}^{(k+1)}=x_{R}^{(k)}-(\frac{\sqrt{5}-1}{2})(x_{R}^{(k)}-x_{L}^{(k)})$$$$x_{2}^{(k+1)}=x_{L}^{(k)}+(\frac{\sqrt{5}-1}{2})(x_{R}^{(k)}-x_{L}^{(k)})$$

  選擇的$x_{1}^{(k+1)}$和$x_{2}^{(k+1)}$一定滿足$$x_{L}^{(k)}<x_{1}^{(k+1)}<x_{2}^{(k+1)}<x_{R}^{(k)}$$

  下面確定新的不確定空間$I^{(k+1)}$

  情況1:若$f(x_{1}^{(k+1)})>f(x_{2}^{(k+1)})$,則$I^{(k+1)}=\left[x_{L}^{(k)},x_{2}^{(k+1)}\right]$

  情況2:若$f(x_{1}^{(k+1)})<f(x_{2}^{(k+1)})$,則$I^{(k+1)}=\left[x_{1}^{(k+1)},x_{R}^{(k)}\right]$

  情況3:若$f(x_{1}^{(k+1)})=f(x_{2}^{(k+1)})$,則$I^{(k+1)}=\left[x_{1}^{(k+1)},x_{2}^{(k+1)}\right]$

  下面記錄下用Python實現二分法和黃金分割法的代碼。

  二分法:

 1 import math
 2 import numpy as np
 3 
 4 
 5 def anyfunction(x):  # 在這里我們定義任意一個指定初始區間內的單峰函數,以x*cos(x)為例
 6     return x*math.cos(x)
 7 
 8 
 9 Low = float(input("Please enter the lowbound: "))
10 High = float(input("Please enter the highbound: "))
11 High = np.pi  # 在這里我們取初始上界為π,如果可以輸入則注釋掉這一行
12 echos = int(input("Please enter the echos: "))  # 迭代次數
13 small = float(input("Please enter the smallvalue: "))  # 公式中的Delta
14 
15 for i in range(1, echos + 1):
16     Lowvalue = anyfunction(0.5*(Low + High - small))
17     Highvalue = anyfunction(0.5*(Low + High + small))
18     print("echos: " + str(i))
19     print('before ' + "Lowbound: " + str(0.5*(Low + High - small)) + " Highbound: " + str(0.5*(Low + High + small)))
20     print('Lowvalue: ' + str(Lowvalue) + ' ' + 'Highvalue: ' + str(Highvalue))
21     if(Lowvalue == Highvalue):
22         Low = 0.5*(Low + High - small)
23         High = 0.5*(Low + High + small)
24     elif(Lowvalue < Highvalue):
25         Low = 0.5*(Low + High - small)
26     else:
27         High = 0.5*(Low + High + small)
28     print("Lowbound: " + str(Low) + " Highbound: " + str(High))

  輸出結果如下:

  5次循環后極值點被限制在[0.7828981633974482,0.8907604338221292]內。

  黃金分割法:

 1 from math import sqrt, cos
 2 import numpy as np
 3 
 4 
 5 def anyfunction(x):  # 同上以函數x*cos(x)為例
 6     return x*cos(x)
 7 
 8 
 9 Low = float(input("Please enter the lowbound: "))
10 High = float(input("Please enter the highbound: "))
11 High = np.pi  # 同上,使用時應該注釋掉
12 echos = int(input("Please enter the echos: "))
13 
14 # 初始化,第一次運算不存在運算簡化
15 uniquevalue = ((sqrt(5)-1)/2)*(High-Low)
16 value1 = anyfunction(High - uniquevalue)
17 value2 = anyfunction(Low + uniquevalue)
18 
19 for i in range(1, echos + 1):
20     print("echos: " + str(i))
21     print('before ' + "Lowbound: " + str(High - uniquevalue) + " Highbound: " + str(Low + uniquevalue))
22     print('value1: ' + str(value1) + ' ' + 'value2: ' + str(value2))
23     # 利用黃金分割法的性質減少一半的運算量
24     if(value1 == value2):
25         Low = High - uniquevalue
26         High = Low + uniquevalue
27         uniquevalue = ((sqrt(5)-1)/2)*(High-Low)
28         value1 = anyfunction(High - uniquevalue)
29         value2 = anyfunction(Low + uniquevalue)
30     elif(value1 < value2):
31         Low = High - uniquevalue
32         uniquevalue = ((sqrt(5)-1)/2)*(High-Low)
33         value1 = value2
34         value2 = anyfunction(Low + uniquevalue)
35     else:
36         High = Low + uniquevalue
37         uniquevalue = ((sqrt(5)-1)/2)*(High-Low)
38         value2 = value1
39         value1 = anyfunction(High - uniquevalue)
40     print("Lowbound: " + str(Low) + " Highbound: " + str(High))

  輸出結果如下:

   5次循環后極值點被限制在[0.7416294238611398,1.0249066567190932]


免責聲明!

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



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