0618法、最速下降法、牛顿法解方程最优解的python实现

最优化理论（二）0618法、最速下降法、牛顿法（python实现）

一、0618法

0618法代码如下：

# @Time : 2023/5/4
# @title: 0618法
# 写传入的函数
def hanshu(x):
    return 2*x**2-x-1

def golden_section(f,a,b,epsilon):
    lamuda = a+0.382*(b-a)
    mu = a+ 0.618 *(b-a)
    k = 0
    result = 0

    while True:
        if hanshu(lamuda) < hanshu(mu):
            if b - lamuda <= epsilon:
                result = mu
                break
            else:
                a = lamuda
                lamuda = mu
                mu = a + 0.618 * (b-a)
        else:
            if mu - a < epsilon:
                result = lamuda
                break
            else:
                b = mu
                mu = lamuda
                lamuda = a + 0.382 * (b-a)
        print(a,b)
        k+=1
    return result

if __name__ == "__main__":
    x = golden_section(hanshu,-1,1,0.16)
    print("%3.2f" % x)

二、最速下降法（二元方程）

代码如下

from sympy import *

x_1 = symbols('x_1')
x_2 = symbols('x_2')

# 函数求解
fun = 2*x_1**2 + x_2**2
print(fun)

# 求导
grad_1 = diff(fun, x_1)
grad_2 = diff(fun, x_2)

# 精度0.1
MaxIter = 100
epsilon = 0.1

# 定义起始点
x_1_value = 1
x_2_value = 1

iter_cnt = 0
current_step_size = 10000

grad_1_value = (float)(grad_1.subs({x_1: x_1_value, x_2: x_2_value}).evalf())
grad_2_value = (float)(grad_2.subs({x_1: x_1_value, x_2: x_2_value}).evalf())

current_obj = fun.subs({x_1: x_1_value, x_2: x_2_value}).evalf()

print(
    '迭代次数: %2d  当前点 (%3.2f, %3.2f)   当前数值: %5.4f     梯度x1: %5.4f     梯度x2 : %5.4f     步长 : %5.4f' % (
    iter_cnt, x_1_value, x_2_value, current_obj, grad_1_value, grad_2_value, current_step_size))

while (abs(grad_1_value) + abs(grad_2_value) >= epsilon):
    iter_cnt += 1
    # find the step size
    t = symbols('t')
    x_1_updated = x_1_value + grad_1_value * t
    x_2_updated = x_2_value + grad_2_value * t
    Fun_updated = fun.subs({x_1: x_1_updated, x_2: x_2_updated})
    grad_t = diff(Fun_updated, t)
    t_value = solve(grad_t, t)[0]  # solve grad_t == 0

    # 更新点的位置
    grad_1_value = (float)(grad_1.subs({x_1: x_1_value, x_2: x_2_value}).evalf())
    grad_2_value = (float)(grad_2.subs({x_1: x_1_value, x_2: x_2_value}).evalf())

    x_1_value = (float)(x_1_value + t_value * grad_1_value)
    x_2_value = (float)(x_2_value + t_value * grad_2_value)

    current_obj = fun.subs({x_1: x_1_value, x_2: x_2_value}).evalf()
    current_step_size = t_value

    print(
        '迭代次数: %2d  当前点 (%3.2f, %3.2f)   当前数值: %5.4f     梯度x1: %5.4f     梯度x2 : %5.4f     步长 : %5.4f' % (
        iter_cnt, x_1_value, x_2_value, current_obj, grad_1_value, grad_2_value, current_step_size))
print("实际上，问题的最优解x*=(%3.2f, %3.2f)^T" % (int(x_1_value), int(x_2_value)))

运行结果：

2*x_1**2 + x_2**2
迭代次数:  0  当前点 (1.00, 1.00)   当前数值: 3.0000     梯度x1: 4.0000     梯度x2 : 2.0000     步长 : 10000.0000
迭代次数:  1  当前点 (-0.11, 0.44)   当前数值: 0.2222     梯度x1: 4.0000     梯度x2 : 2.0000     步长 : -0.2778
迭代次数:  2  当前点 (-0.11, 0.44)   当前数值: 0.2222     梯度x1: -0.4444     梯度x2 : 0.8889     步长 : 0.0000
迭代次数:  3  当前点 (0.07, 0.07)   当前数值: 0.0165     梯度x1: -0.4444     梯度x2 : 0.8889     步长 : -0.4167
迭代次数:  4  当前点 (0.07, 0.07)   当前数值: 0.0165     梯度x1: 0.2963     梯度x2 : 0.1481     步长 : 0.0000
迭代次数:  5  当前点 (-0.01, 0.03)   当前数值: 0.0012     梯度x1: 0.2963     梯度x2 : 0.1481     步长 : -0.2778
迭代次数:  6  当前点 (-0.01, 0.03)   当前数值: 0.0012     梯度x1: -0.0329     梯度x2 : 0.0658     步长 : 0.0000
实际上，问题的最优解x*=(0.00, 0.00)^T

三、牛顿法（二元方程）

代码如下：

针对上面的例题，由于求导后的数值在程序的Hesse matrix中无法表示为浮点数，而是未知数。所以导致不能对hesse矩阵进行取逆，这里我直接输入了一个求导后的数值。

from sympy import *
import numpy as np
# 定义符号
x1,x2 = symbols('x1, x2')
# 定义所求函数
f = (x1-1)**4 + x2**2
# 求解梯度值
def get_grad(f, X):
    # 计算一阶导数
    f1 = diff(f, x1)
    f2 = diff(f, x2)
    # 代入具体数值计算
    grad = np.array([[f1.subs([(x1, X[0]), (x2, X[1])])],
                   [f2.subs([(x1, X[0]), (x2, X[1])])]])
    return grad
# 求解Hession矩阵
def get_hess(f, X):
    # 计算二次偏导
    f1 = diff(f, x1)
    f2 = diff(f, x2)
    # f11 = diff(f,x1,2)
    f22 = diff(f,x2,2)
    f12 = diff(f1,x2)
    f21 = diff(f2,x1)
    # 计算具体数值计算
    hess = np.array([[12,
                        f12.subs([(x1,X[0]), (x2,X[1])])],

                        [f21.subs([(x1,X[0]), (x2,X[1])]),
                        f22.subs([(x1,X[0]), (x2,X[1])])]])
    # 转换数值类型为了后续求逆矩阵
    hess = np.array(hess, dtype = 'float')
    return hess
# 牛顿迭代
def newton_iter(X0, err):
    # 记录迭代次数
    count = 0
    X1 = np.array([[0],[0]])
    while True:
        # 得到两次迭代结果差值
        X2 = X0 - X1
        if sqrt(X2[0]**2 + X2[1]**2) <= err:
            break
        else:
            hess = get_hess(f, X0)
            # 得到Hession矩阵的逆
            hess_inv = np.linalg.inv(hess)
            grad = get_grad(f, X0)
            X1 = X0 #保存上次迭代结果
            # 迭代公式
            X0 = X1 - np.dot(hess_inv, grad)
            count += 1
            print('第',count,'次迭代:','x1=',X0[0],'x2=',X0[1])
    print('迭代次数为:',count)
    print('迭代结果为',X0)
# 设置初始点
X0 = np.array([[1],[1]])
err = 1e-10
newton_iter(X0, err)

运行结果

第 1 次迭代: x1= [1] x2= [0]
第 2 次迭代: x1= [1] x2= [0]
迭代次数为: 2
迭代结果为 [[1]
 [0]]