最优化理论 无约束规划 FR算法 共轭梯度法求极小点

前言

调用库
sympy，符号计算库，可以用来求偏导、带值计算、求解方程等。

import sympy as sp
import numpy as np

针对规划问题

取初始点x0=(3, -1, 0, 1) 设置精度范围 e = 0.05
这里精度小了会发现迭代次数非常多，我设置了10^(-3) 能迭代两百多次，这里设置为0.05迭代大概四十八次结束，但是会发现两次计算结果不一样，我估计是因为精度要求不同下降得到的极小点不同，因而选择了其他的局部解。

共轭梯度法基本步骤

①取初始点 x0，设置精度，令迭代次数 k=1；
②递归终点（循环终点）：||▽f(x(k)) || <= e，停止计算，x(k)作为无约束规划问题的解；否则
置（下面两个式子等价)

再计算
③一维搜索求解问题

然后，置 x(k+1) = x(k) + α(k)d(k)。
注意，这里括号内的k为迭代次序，并非变量；
④置 k += 1，转步骤②。

代码

import sympy as sp
import numpy as np

x1 = sp.Symbol('x1')
x2 = sp.Symbol('x2')
x3 = sp.Symbol('x3')
x4 = sp.Symbol('x4')
f = (x1 + 10 * x2) ** 2 + 5 * ((x3 - x4) ** 2) + (x2 - 2 * x3) ** 4 + 10 * ((x1 - x4) ** 4)
x0 = np.array([3, -1, 0, 1])
dfx1 = sp.diff(f, x1)
dfx2 = sp.diff(f, x2)
dfx3 = sp.diff(f, x3)
dfx4 = sp.diff(f, x4)
grad_f = np.array([dfx1, dfx2, dfx3, dfx4])  # 梯度
d2f1 = [sp.diff(dfx1, x1), sp.diff(dfx2, x1), sp.diff(dfx3, x1), sp.diff(dfx4, x1)]
d2f2 = [sp.diff(dfx1, x2), sp.diff(dfx2, x2), sp.diff(dfx3, x2), sp.diff(dfx4, x2)]
d2f3 = [sp.diff(dfx1, x3), sp.diff(dfx2, x3), sp.diff(dfx3, x3), sp.diff(dfx4, x3)]
d2f4 = [sp.diff(dfx1, x4), sp.diff(dfx2, x4), sp.diff(dfx3, x4), sp.diff(dfx4, x4)]
hesse = np.array([d2f1, d2f2, d2f3, d2f4])  # Hesse矩阵
e = 0.05
k = 1


def endow_grad(x):  # 给梯度赋值
    dicx = {'x1': x[0], 'x2': x[1], 'x3': x[2], 'x4': x[3]}
    grad_fx = np.array([dfx1.evalf(subs=dicx),
                        dfx2.evalf(subs=dicx),
                        dfx3.evalf(subs=dicx),
                        dfx4.evalf(subs=dicx)],
                       dtype=np.float64)
    return grad_fx


def endow_hesse(x):  # 给汉斯矩阵赋值
    dicx = {'x1': x[0], 'x2': x[1], 'x3': x[2], 'x4': x[3]}
    hessex = np.array(
        [[d2f1[0].evalf(subs=dicx), d2f1[1].evalf(subs=dicx), d2f1[2].evalf(subs=dicx), d2f1[3].evalf(subs=dicx)],
         [d2f2[0].evalf(subs=dicx), d2f2[1].evalf(subs=dicx), d2f2[2].evalf(subs=dicx), d2f2[3].evalf(subs=dicx)],
         [d2f3[0].evalf(subs=dicx), d2f3[1].evalf(subs=dicx), d2f3[2].evalf(subs=dicx), d2f3[3].evalf(subs=dicx)],
         [d2f4[0].evalf(subs=dicx), d2f4[1].evalf(subs=dicx), d2f4[2].evalf(subs=dicx), d2f4[3].evalf(subs=dicx)]],
        dtype=np.float64)
    return hessex


def FR(pre_x, cur_x, pre_d):
    global k
    grad_fx1, grad_fx2 = endow_grad(pre_x), endow_grad(cur_x)  # 获取前一个x对应的梯度值和当前x对应的梯度值
    hesse_x = endow_hesse(cur_x)  # 获取当前x对应的已代入值得汉斯矩阵
    b = ((np.linalg.norm(grad_fx2)) ** 2) / (np.linalg.norm(grad_fx1) ** 2)  # 求 β
    dk = -grad_fx2 + b * pre_d  # 求当前迭代次数对应的方向向量
    a = -(np.dot(dk.T, grad_fx2) / np.dot(np.dot(dk.T, hesse_x), dk))  # 求最优步长
    new_x = cur_x + a * dk  # 求出新一轮的 x
    grad_fx3 = endow_grad(new_x)
    k += 1
    print("(%d) x =" % k, new_x)
    # print(grad_fx1, grad_fx2)
    # print(hesse_x)
    print("范数值为：", np.linalg.norm(grad_fx3))
    if np.linalg.norm(grad_fx3) <= e:
        print("最终解为：", new_x)
        return
    else:
        FR(cur_x, new_x, dk)


grad_fx0 = endow_grad(x0)
hessex0 = endow_hesse(x0)
d0 = -grad_fx0
a0 = -(np.dot(d0.T, grad_fx0) / np.dot(np.dot(d0.T, hessex0), d0))
x1 = x0 + a0 * d0
print("(%d)  x1 =" % k, x1)
print("d0 =", d0)
print("a0 =", a0)
FR(x0, x1, d0)

附言

值得一提的是，当时写共轭梯度法，因为不太熟悉sympy库内函数的调用，结果一直报这个错误，当然，上面贴的代码是正确的。

我就很疑惑，为什么写Newton法求解无约束规划问题时这样写没有报错，但是后来将报错部分改成这段就没事了。没错就是加个单引号，但是Newton法并没有加单引号也并没有报错，有人告诉我原因吗。
dicx = {'x1': x[0], 'x2': x[1], 'x3': x[2], 'x4': x[3]}
改个bug改了三天可还行，网络上都搜不到类似问题。

此外，最优步长的求解使用一维搜索是需要目标函数满足特定条件的，不能任何时候都使用，否则会落入局部解甚至得到不正确解。