黄金分割法的python实现

在准备学优化方法期末考试的时候拿到习题没有答案,就想写一个程序然后和自己算出来的答案对一下就写了个程序,考完试了留着也没用索性写一篇出来造福学弟学妹,黄金分割法通俗的讲就是给一个区间给一个函数在这个区间里求函数的最小值

1.一维搜索黄金分割法

这里搬一下书上的概念书上的概念严谨但确实不好懂做两道例题就好了其实
考察一维搜索问题:

{ m i n φ ( λ ) s . t . a 1 ≤ λ ≤ b 1 \begin{cases} min&\varphi(\lambda)\\ s.t.&a_1\le\lambda\le b_1 \end{cases} {mins.t.φ(λ)a1λb1
其中 φ ( λ ) \varphi(\lambda) φ(λ)为凸函数

1 。 1^。 1 ε > 0 \varepsilon>0 ε>0为允许的搜索区间最后的长度令 α \alpha α=0.618则黄金分割法计算步骤如下

λ 1 = a 1 + ( 1 − α ) ( b 1 − a 1 ) φ ( λ 1 ) μ 1 = a 1 + α ( b 1 − a 1 ) φ ( μ 1 ) \lambda_1=a_1+(1-\alpha)(b_1-a_1) \quad \varphi(\lambda_1)\\ \mu_1=a_1+\alpha(b_1-a_1) \qquad \varphi(\mu_1) λ1=a1+(1α)(b1a1)φ(λ1)μ1=a1+α(b1a1)φ(μ1)
k = 1 k=1 k=1

2 。 2^。 2 b k − a k < ε b_k-a_k<\varepsilon bkak<ε,则计算结束,最优解 λ ∗ ∈ [ a k , b k ] \lambda^*\in[a_k,b_k] λ[ak,bk],可取 λ ∗ ≈ ( 1 / 2 ) ( a k + b k ) \lambda^*\approx(1/2)(a_k+b_k) λ(1/2)(ak+bk);否则 φ ( λ k ) > φ ( μ k ) \varphi(\lambda_k)>\varphi(\mu_k) φ(λk)>φ(μk), 则转 3 。 3^。 3;若 φ ( λ k ) ≤ φ ( μ k ) \varphi(\lambda_k)\le\varphi(\mu_k) φ(λk)φ(μk),则转 4 。 4^。 4

3 。 3^。 3 a k + 1 = λ k , b k + 1 = b k a_{k+1}=\lambda_k,b_{k+1}=b_k ak+1=λk,bk+1=bk,再令 λ k + 1 = μ k , μ k + 1 = a k + 1 + α ( b k + 1 − a k + 1 ) \lambda_{k+1}=\mu_k, \mu_{k+1}=a_{k+1}+\alpha(b_{k+1}-a_{k+1}) λk+1=μk,μk+1=ak+1+α(bk+1ak+1),计算 φ ( μ k + 1 ) \varphi(\mu_{k+1}) φ(μk+1) 5 。 5^。 5

4 。 4^。 4 a k + 1 = a k , b k + 1 = μ k a_{k+1}=a_k,b_{k+1}=\mu_k ak+1=ak,bk+1=μk,再令 μ k + 1 = λ k , λ k + 1 = a k + 1 + ( 1 − α ) ( b k + 1 − a k + 1 ) \mu_{k+1}=\lambda_k, \lambda_{k+1}=a_{k+1}+(1-\alpha)(b_{k+1}-a_{k+1}) μk+1=λk,λk+1=ak+1+(1α)(bk+1ak+1),计算 φ ( μ k + 1 ) \varphi(\mu_{k+1}) φ(μk+1) 5 。 5^。 5

5 。 5^。 5 k = k + 1 k=k+1 k=k+1,返回 2 。 2^。 2

2.python代码实现

例:求解 min ⁡ φ ( λ ) = λ 2 + 2 λ , s . t . − 3 ≤ λ ≤ 5 \min\varphi(\lambda)=\lambda^2+2\lambda,s.t. -3\le\lambda\le 5 minφ(λ)=λ2+2λ,s.t.3λ5
ε = 0.2 \varepsilon=0.2 ε=0.2取精度为0.001
取三位小数用python的round(x,3)实现

# 只需改如下四个参数即可
f = lambda x: x**2 + 2*x
# 定义变量
e = 0.2       # 终止线
a = -3           # 左边
b = 5          # 右边


# 黄金分割法计算法则
lanb = lambda a, b : a + 0.382* (b -a)
mui = lambda a, b : a + 0.618* (b -a)

Mui = round(mui(a, b),3)
Mui_value = round(f(Mui), 3)

Lamb = round(lanb(a, b) ,3)
Lamb_value= round(f(Lamb), 3)

print(f'a={a} b={b}')
print(f'lanb:{Lamb}')
print(f'mui:{Mui}')
print(f'f(lanb)= {Lamb_value}')
print(f'f(mui)= {Mui_value}')
print() #换行最后看着舒服

while b - a > e: #判断是否到达终止线
    if Lamb_value < Mui_value:
        print('f(lanb) < f(mui)')
        print()
        b = Mui
        Mui = Lamb
        Lamb = round(lanb(a, b), 3)
    else:
        print('f(mui) < f(lanb) ')
        print()
        a = Lamb
        Lamb = Mui
        Mui = round(mui(a, b), 3)

    Mui_value = round(f(Mui), 3)
    Lamb_value = round(f(Lamb), 3)
    print(f'a={a} b={b}')
    print(f'lanb:{Lamb}')
    print(f'mui:{Mui}')
    print(f'f(lanb)= {Lamb_value}')
    print(f'f(mui)= {Mui_value}')


print(f'b-a绝对值为{abs(b -a)} 小于终止线 算法停止')
print(f'X⭐包含于[{a}, {b}]')
print(f'X⭐ = {round((a + b)/ 2, 3)}')

3.运行结果

a=-3 b=5
lanb:0.056
mui:1.944
f(lanb)= 0.115
f(mui)= 7.667
f(lanb) < f(mui)

a=-3 b=1.944
lanb:-1.111
mui:0.056
f(lanb)= -0.988
f(mui)= 0.115
f(lanb) < f(mui)

a=-3 b=0.056
lanb:-1.833
mui:-1.111
f(lanb)= -0.306
f(mui)= -0.988
f(mui) < f(lanb) 

a=-1.833 b=0.056
lanb:-1.111
mui:-0.666
f(lanb)= -0.988
f(mui)= -0.888
f(lanb) < f(mui)

a=-1.833 b=-0.666
lanb:-1.387
mui:-1.111
f(lanb)= -0.85
f(mui)= -0.988
f(mui) < f(lanb) 

a=-1.387 b=-0.666
lanb:-1.111
mui:-0.941
f(lanb)= -0.988
f(mui)= -0.997
f(mui) < f(lanb) 

a=-1.111 b=-0.666
lanb:-0.941
mui:-0.836
f(lanb)= -0.997
f(mui)= -0.973
f(lanb) < f(mui)

a=-1.111 b=-0.836
lanb:-1.006
mui:-0.941
f(lanb)= -1.0
f(mui)= -0.997
f(lanb) < f(mui)

a=-1.111 b=-0.941
lanb:-1.046
mui:-1.006
f(lanb)= -0.998
f(mui)= -1.0
b-a绝对值为0.17 小于终止线 算法停止
X⭐包含于[-1.111, -0.941]
X⭐ = -1.026

结束

用的书上的例子,书上的 λ 2 \lambda_2 λ2算成-1.112算错了应该是-1.111所以后面的结果不一样
在这里插入图片描述

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐