最优化大作业 （一）：一维搜索法

• 问题描述

对以下优化问题

（1）黄金分割法分别取初始搜索区间[-2,0]和[0,3]；

（2）牛顿法分别取初始点或；

（3）二次插值法分别取初始点 ，，和，，。

 

• 基本原理

（1）黄金分割法

 

图1  黄金分割法流程图

（2）牛顿法

图2  牛顿法流程图

 

（3）二次插值法

 

图3  二次插值法流程图

 

• 实验结果

（1）黄金分割法

迭代 次数	1	2	3	4	5	6	7	8	9
	-0.76	-0.47	-0.76	-0.94	-0.88	-0.94	-0.99	-0.97	-0.99
	-1.24	-0.76	-0.94	-1.06	-0.94	-0.99	-1.01	-0.99	-1.00
	-1.00	-0.62	-0.85	-1.00	-0.91	-0.97	-1.00	-0.98	-0.99
	-5.00	-3.20	-4.67	-5.00	-4.87	-4.98	-5.00	-4.99	-5.00

表1  (a,b)为[-2, 0]时的迭代过程

 

迭代次数	1	2	3	4	5	6	7	8	9	10	11
	1.85	2.29	1.85	2.02	2.12	2.02	2.06	2.02	2.00	2.01	2.00
	1.15	1.85	1.58	1.85	2.02	1.96	2.02	2.00	1.98	2.00	1.99
	1.50	2.07	1.72	1.94	2.07	1.99	2.04	2.01	1.99	2.00	1.99
	-25.3	-31.8	-29.5	-31.8	-31.8	-32.0	-31.9	-32.0	-32.0	-32.0	-32.0

表2  (a,b)为[0, 3]时的迭代过程

 

图4  黄金分割法迭代过程图

 

（2）牛顿法

迭代次数	1(初始点)	2	3	4	5
x	-1.20	-1.03	-1.00	-1.00	-1.00
y	-4.14	-4.97	-4.99	-4.99	-5.00

表3  x0为-1.2时的迭代过程

 

迭代次数	1(初始点)	2	3	4	5
x	2.50	2.12	2.01	2.00	2.00
y	-20.312	-31.37	-31.99	-31.99	-32.00

表4  x0为2.5时的迭代过程

图5  牛顿法迭代过程图

 

（3）二次插值法

迭代次数	1	2
x	-0.98	-5.00
y	-0.98	-4.99

表5  x1=-1.2，x2=-1.1，x3=-0.8时的迭代过程

 

迭代次数	1	2	3	4	5	6	7	8	……	36
x	1.03	7.62	-0.21	0.44	0.79	0.97	1.07	1.12	……	1.91
y	-13.65	7646.09	-0.50	-2.58	-8.34	-12.40	-14.66	1.16	……	-31.71

表6  -1.5，x2=1.7，x3=2.5时的迭代过程

 

 

图6  插值法迭代过程图

 

由于设置ε=0.01，部分结果受精度影响未能迭代到最优值点。

 

• 代码展示

import matplotlib.pyplot as plt
import numpy as np

# 定义搜索函数
def f(t):
    y = 3*np.power(t, 4) - 4*np.power(t, 3) - 12*np.power(t, 2)
    return y

#定义优化函数的一阶导数
def f1(t):
    y = 12*np.power(t, 3) - 12*np.power(t, 2) - 24*t
    return y

#定义优化函数的二阶导数
def f2(t):
    y = 36*np.power(t, 2) - 24*t - 24
    return y

# 黄金分割法
def gold_ratio(X):
    print("(a,b)为 ", X)
    beta = 0.382
    t2 = X[0]+0.382*(X[1]-X[0])
    t1 = X[0] + X[1] - t2

    count = 1
    while True:
        plt.scatter((t1+t2)/2, f((t1+t2)/2), c='b')
        print("迭代次数：", count, " (t1,t2,(t1+t2)/2,f): (%0.2f,%0.2f,%0.2f,%0.2f)" % (t1, t2, (t1+t2)/2, f((t1+t2)/2)))
        if abs(t1-t2) < 0.01:
            print("迭代次数： ", count, " 最终结果:", "t*=", (t1+t2)/2, ", f(t*)=", f((t1+t2)/2))
            return
        else:
            count += 1
            if f(t1) < f(t2):
                X[0] = t2
                t2 = t1
                t1 = X[0] + X[1] - t2

            else:
                X[1] = t1
                t1 = t2
                t2 = X[0] + beta*(X[1]-X[0])

#牛顿法
def newton(t0):
    count = 1
    plt.scatter(t0, f(t0), c='b')
    print("迭代次数：", count, " 结果:", "t*=%0.2f, f(t*)=%0.2f"%(t0, f(t0)))
    while True:
        count += 1
        t = t0 - f1(t0) / f2(t0)
        plt.scatter(t, f(t), c='b')
        print("迭代次数：", count, " 结果:", "t*=%0.2f, f(t*)=%0.2f"%(t, f(t)))
        if abs(t-t0) < 0.01:
            return
        else:
            t0 = t

def draw():
    plt.figure()
    x = np.array([i/10 for i in range(-20,30)])
    y = f(x)
    plt.plot(x, y)

def gold_ratio_main():
    a = [-2, 0]
    b = [0, 3]
    draw()
    gold_ratio(a)
    gold_ratio(b)
    plt.show()

def newton_main():
    draw()
    print("当x0=-1.2：")
    newton(-1.2)
    print("当x0=2.5：")
    newton(2.5)
    plt.show()

#定义求二次拟合函数极小值点的函数。即书中4.10式
def cal(t0, t1, t2):
    a = t0**2
    b = t1**2
    c = t2**2
    molecular = (a-c)*f(t1) + (c-b)*f(t0) + (b-a)*f(t2)
    denominator = (t0-t2)*f(t1) + (t2-t1)*f(t0) + (t1-t0)*f(t2)
    mini = molecular/(2*denominator)
    return mini

#二次插值法
def interpolation(t0,t1,t2):
    count = 1
    while True:
        t = cal(t0, t1, t2)
        if f(t)<0:
            plt.scatter(t, f(t), c='b')
        print("迭代次数：", count, " 结果:", "t*=%0.2f, f(t*)=%0.2f" % (t, f(t)))
        count += 1
        if abs(t-t0) < 0.01:
            return
        else:
            if t > t0:
                if f(t) <= f(t0):
                    t1=t0;t0=t
                else:
                    t2 = t
            else:
                if f(t) <= f(t0):
                    t2=t0;t0=t
                else:
                    t1 = t

def interpolation_main():
    draw()
    print("当x1=-1.2，x2=-1.1，x3=-0.8时，结果为：")
    interpolation(-1.2, -1.1, -0.8)
    print("当x1=-1.5，x2=1.7，x3=2.5时，结果为：")
    interpolation(-1.5, 1.7, 2.5)
    plt.show()

if __name__ == "__main__":
    gold_ratio_main()
    newton_main()
    interpolation_main()