优化方法上机实验

最新推荐文章于 2022-04-22 14:48:32 发布

红帽小生

最新推荐文章于 2022-04-22 14:48:32 发布

阅读量893

点赞数 3

分类专栏：科学计算文章标签：算法 python

本文链接：https://blog.csdn.net/abraham76/article/details/115695732

版权

科学计算专栏收录该内容

6 篇文章 0 订阅

订阅专栏

本文介绍了数值优化中的几种经典方法，包括最速下降法、牛顿法、共轭梯度法和拟牛顿法（BFGS）。通过Python实现展示了这些方法在无约束非线性规划问题中的应用，并给出了约束优化问题的增广拉格朗日方法。文章强调了这些算法在确定步长和方向选择上的差异，以及它们在不同问题上的效率和适用性。

摘要由CSDN通过智能技术生成

引言

求解无约束非线性规划问题，一般是求目标函数 $f (x)$ 的极小值。至于求极大值问题可以转化为求 $- f (x)$ 的极小值问题。
一般情况下，不能用解析方法求得准确解，只能用数值方法逐步求其近似解。其思想很简单，从一个初始估计点 $x^{(0)}$ 出发，计算一系列点 $x^{(k)}$ ，这个点列的极限就是 $f (x)$ 的一个极小点。
那么由 $x^{(k)}$ 如何求得 $x^{(k+1)}$ 呢？因为 $x^{(k+1)}-x^{(k)}$ 是一个向量，而向量是由它的方向和长度确定的。那么有：
$x^{(k+1)}-x^{(k)}=\lambda_{k}s^{k}$
或写成这样：
$x^{(k+1)}=x^{(k)}+\lambda_{k}s^{k}$
只要确定了步长 $\lambda_{k}$ 和方向 $s^{k}$ 就可以进行迭代。不同的算法仅仅是在方向和步长的选取方法上不同。
来看下题目：

最速下降法

最速下降法也叫梯度法。是用来求多个变量函数极值问题的最早的一种方法。最速下降的方向为负梯度方向，步长采用线性搜索来确定。有精确线性搜索(一维搜索)，也有不精确线性搜索。后者更常用。主要有Goldstein不精确线性搜索和Wolfe法线性搜索。书上较多采用了精确一维搜索。
使用python编程，代码如下：

import random
import numpy as np
import matplotlib.pyplot as plt
"""
最速下降法
目标函数 f(x)=10*(x(1)-1).^2+(x(2)+1).^4
梯度 g(x)=( 20*(x(1)-1) ,4*(x(2)+1).^3)^(T)
"""
def goldsteinsearch(f,df,d,x,alpham,rho,t):
    '''
    线性搜索子函数
    数f，导数df，当前迭代点x和当前搜索方向d
    '''
    flag = 0
    a = 0
    b = alpham
    fk = f(x)
    gk = df(x)
    phi0 = fk
    dphi0 = np.dot(gk, d)
    # print(dphi0)
    alpha=b*random.uniform(0,1)
 
    while(flag==0):
        newfk = f(x + alpha * d)
        phi = newfk
        # print(phi,phi0,rho,alpha ,dphi0)
        if (phi - phi0 )<= (rho * alpha * dphi0):
            if (phi - phi0) >= ((1 - rho) * alpha * dphi0):
                flag = 1
            else:
                a = alpha
                b = b
                if (b < alpham):
                    alpha = (a + b) / 2
                else:
                    alpha = t * alpha
        else:
            a = a
            b = alpha
            alpha = (a + b) / 2
    return alpha

def objective_function(x):
    return 10*(x[0]-1)**2+(x[1]+1)**4

def jacobian(x):
    return np.array([20*(x[0]-1),4*(x[1]+1)**3])

def steepest(x0):
    print('初始点为:')
    print(x0,'\n')
    imax = 20000
    W = np.zeros((2, imax))
    epo=np.zeros((2, imax))
    W[:, 0] = x0
    i = 1
    x = x0
    grad = jacobian(x)
    delta = sum(grad ** 2)  # 初始误差
 
    f=open("最速.txt",'w')
 
    while i < imax and delta > 10 ** (-4):
        p = -jacobian(x)
        x0 = x
        alpha = goldsteinsearch(objective_function, jacobian, p, x, 1, 0.1, 2)
        x = x + alpha * p
        W[:, i] = x
        if i % 5 == 0:
 
            epo[:,i] =np.array((i,delta))
            f.write(str(i)+"        "+str(delta)+"\n")
            print(i,np.array((i,delta)))
        grad = jacobian(x)
        delta = sum(grad ** 2)
        i = i + 1
    print("迭代次数为:", i)
    print("近似最优解为:")
    print(x, '\n')
    W = W[:, 0:i]  # 记录迭代点
    return [W,epo]
 
if __name__=="__main__":
    X1 = np.arange(-1.5, 1.55, 0.05)
    X2 = np.arange(-3.5, 3.55, 0.05)
    [x1, x2] = np.meshgrid(X1, X2)
    f = 10 *(x1 - 1 ) ** 2 + (x2+1) ** 4  # 目标函数
    plt.contour(x1, x2, f, 20)  # 画出函数的20条轮廓线
    x0 = np.array([0, 0])   #初始点
    list_out = steepest(x0)
    W=list_out[0]
    epo=list_out[1]
    plt.plot(W[0, :], W[1, :], 'g*-')  # 画出迭代点收敛的轨迹
    plt.show()

牛顿法

一元函数和多元函数的牛顿法思想都是一样的，不再赘述。直接上代码：

import random
import numpy as np
import matplotlib.pyplot as plt
 
'''
牛顿法求解无约束最优化问题
'''
def dampnm(fun,gfun,hess,x0):
    #x0是初始点，fun，gfun和hess分别是目标函数值，梯度，海森矩阵的函数
    maxk = 500
    rho = 0.55
    sigma = 0.4
    k = 0
    epsilon = 1e-4     
    f=open("牛顿.txt","w")
    W = np.zeros((2, 20000))
    while k < maxk:
        W[:, k] = x0
        gk = gfun(x0)
        Gk = hess(x0)
        dk = -1.0*np.linalg.solve(Gk,gk)
        print(k, np.linalg.norm(dk))
        f.write(str(k)+'   '+str(np.linalg.norm(gk))+"\n")
        if np.linalg.norm(dk) < epsilon:
            break
        x0 =x0+ dk
        k += 1
    W = W[:, 0:k + 1]  # 记录迭代点
    return x0,fun(x0),k,W
# 函数表达式fun
fun = lambda x:10*(x[0]-1)**2 + (x[1]+1)**4
# 梯度向量 gfun
gfun = lambda x:np.array([20*(x[0]-1), 4*(x[1]+1)**3])
# 海森矩阵 hess
hess = lambda x:np.array([[ 20 , 0 ],[0,12*(x[1]+1)**2]])
if __name__=="__main__":
    X1 = np.arange(-1.5, 1.55, 0.05)
    X2 = np.arange(-3.5, 2.05, 0.05)
    [x1, x2] = np.meshgrid(X1, X2)
    f = 10 *(x1 - 1 ) ** 2 + (x2+1) ** 4   # 目标函数
    plt.contour(x1, x2, f, 40)  # 画出函数的20条轮廓线
    x0 = np.array([0, 0])
    out=dampnm(fun, gfun, hess, x0)
    print(out[2],out[0])
    W = out[3]
    print(W[:,:])
    plt.plot(W[0, :], W[1, :], 'g*-')
    plt.show()

共轭梯度法

最速下降法每次都是在局部最优的方向上选取了最优的步长，但是在其后面的迭代过程中可能会出现相同的方向，这样使我们在迭代过程中多次对一个方向进行迭代，在极值点附近的迭代轨迹为锯齿波，影响了收敛速度。
共轭梯度法解决了这个问题，每一次迭代都将该方向优化为最优方向，理论上对n维问题，只需求出n个方向上的极小值。

import random
import numpy as np
import matplotlib.pyplot as plt
'''
共轭梯度法
  线性搜索子函数
  数f，导数df，当前迭代点x和当前搜索方向d，t试探系数>1，
'''
def goldsteinsearch(f,df,d,x,alpham,rho,t):
    flag = 0
    a = 0
    b = alpham
    fk = f(x)
    gk = df(x)
    phi0 = fk
    dphi0 = np.dot(gk, d)
    alpha=b*random.uniform(0,1)
    while(flag==0):
        newfk = f(x + alpha * d)
        phi = newfk
        if (phi - phi0 )<= (rho * alpha * dphi0):
            if (phi - phi0) >= ((1 - rho) * alpha * dphi0):
                flag = 1
            else:
                a = alpha
                b = b
                if (b < alpham):
                    alpha = (a + b) / 2
                else:
                    alpha = t * alpha
        else:
            a = a
            b = alpha
            alpha = (a + b) / 2
    return alpha
'''
 线性搜索子函数
 数f，导数df，当前迭代点x和当前搜索方向d
σ∈(ρ,1)=0.75
'''
def Wolfesearch(f,df,d,x,alpham,rho,t):
    sigma=0.75
    flag = 0
    a = 0
    b = alpham
    fk = f(x)
    gk = df(x)
    phi0 = fk
    dphi0 = np.dot(gk, d)
    alpha=b*random.uniform(0,1)
    while(flag==0):
        newfk = f(x + alpha * d)
        phi = newfk
        if (phi - phi0 )<= (rho * alpha * dphi0):
            if (phi - phi0) >= ((1 - rho) * alpha * dphi0):
                flag = 1
            else:
                a = alpha
                b = b
                if (b < alpham):
                    alpha = (a + b) / 2
                else:
                    alpha = t * alpha
        else:
            a = a
            b = alpha
            alpha = (a + b) / 2
    return alpha
 
def frcg(fun,gfun,x0):
 # x0是初始点，fun和gfun分别是目标函数和梯度
 # x,val分别是近似最优点和最优值，k是迭代次数
 # dk是搜索方向，gk是梯度方向
 # epsilon是预设精度，np.linalg.norm(gk)求取向量的二范数
    maxk = 5000
    rho = 0.6
    sigma = 0.4
    k = 0
    epsilon = 1e-4
    n = np.shape(x0)[0]
    itern = 0
    W = np.zeros((2, 20000))
    f = open("共轭.txt", 'w')
    while k < maxk:
        W[:, k] = x0
        gk = gfun(x0)
        itern += 1
        itern %= n
        if itern == 1:
            dk = -gk
        else:
            beta = 1.0 * np.dot(gk, gk) / np.dot(g0, g0)
            dk = -gk + beta * d0
            gd = np.dot(gk, dk)
            if gd >= 0.0:
                dk = -gk
        if np.linalg.norm(gk) < epsilon:
            break
        alpha=goldsteinsearch(fun,gfun,dk,x0,1,0.1,2)
        x0=x0+alpha*dk
        f.write(str(k)+' '+str(np.linalg.norm(gk))+"\n")
        print(k,alpha)
        g0 = gk
        d0 = dk
        k += 1
    W = W[:, 0:k+1] # 记录迭代点
    return [x0, fun(x0), k,W]
 
def fun(x):
    return 10*(x[0]-1)**2+(x[1]+1)**4
def gfun(x):
    return np.array([20*(x[0]-1),4*(x[1]+1)**3])

if __name__=="__main__":
    X1 = np.arange(-1.5, 1.55, 0.05)
    X2 = np.arange(-3.5, 3.55, 0.05)
    [x1, x2] = np.meshgrid(X1, X2)
    f = 10 *(x1 - 1 ) ** 2 + (x2+1) ** 4   # 目标函数
    plt.contour(x1, x2, f, 20) # 画出函数的20条轮廓线
    x0 = np.array([0, 0])
    x=frcg(fun,gfun,x0)
    print(x[0],x[2])
    W=x[3]
    plt.plot(W[0, :], W[1, :], 'g*-') # 画出迭代点收敛的轨迹
    plt.show()

拟牛顿法(使用BFGS方法)

拟牛顿法属于变尺度算法的一种。拟牛顿法的本质思想是改善牛顿法每次需要求解复杂的Hessian矩阵的逆矩阵的缺陷，它使用正定矩阵来近似Hessian矩阵的逆，从而简化了运算的复杂度。拟牛顿法和最速下降法一样只要求每一步迭代时知道目标函数的梯度。通过测量梯度的变化，构造一个目标函数的模型使之足以产生超线性收敛性。这类方法大大优于最速下降法，尤其对于困难的问题。另外，因为拟牛顿法不需要二阶导数的信息，所以有时比牛顿法更为有效。如今，优化软件中包含了大量的拟牛顿算法用来解决无约束，约束，和大规模的优化问题。

"""
拟牛顿法
"""
import numpy as np
from sympy import symbols
# 首先定义二维向量内的元素
x1 = symbols("x1")
x2 = symbols("x2")
def jacobian(f, x):
    """
    求函数一阶导
    :param f: 原函数
    :param x: 初始值
    :return: 函数一阶导的值
    """
    grandient = np.array([20*(x[0]-1),4*(x[1]+1)**3], dtype=float)
    return grandient


def bfgs_newton(f, x, iters):
    """
    实现BFGS拟牛顿法
    :param f: 原函数
    :param x: 初始值
    :param iters: 遍历的最大epoch
    :return: 最终更新完毕的x值
    """
    # 步长。设为1才能收敛，小于1不能收敛
    learning_rate = 1
    # 初始化B正定矩阵
    B = np.eye(2)
    x_len = x.shape[0]
    # 一阶导g的第二范式的最小值（阈值）
    epsilon = 1e-4
    for i in range(1, iters):
        g = jacobian(f, x)
        if np.linalg.norm(g) < epsilon:
            break
        p = np.linalg.solve(B, g)
        # 更新x值
        x_new = x - p*learning_rate
        print("第" + str(i) + "次迭代后的结果为:", x_new)
        g_new = jacobian(f, x_new)
        y = g_new - g
        k = x_new - x
        y_t = y.reshape([x_len, 1])
        Bk = np.dot(B, k)
        k_t_B = np.dot(k, B)
        kBk = np.dot(np.dot(k, B), k)
        # 更新B正定矩阵。完全按照公式来计算
        B = B + y_t*y/np.dot(y, k) - Bk.reshape([x_len, 1]) * k_t_B / kBk
        x = x_new
    return x
if __name__ == "__main__":
    x0 = np.array([0, 0], dtype=float)
    f = 10 *(x1 - 1 ) ** 2 + (x2+1) ** 4  # 目标函数
    print(bfgs_newton(f, x0, 2000))

约束优化方法

前面讨论的都是无约束的优化问题。但更常见的数学规划问题是约束优化问题。目前尚没有一种对一切问题都普遍有效的算法，而且现有的算法求得的解多是局部最优解。
下面看一下题目：

增广拉格朗日方法

增广拉格朗日法是惩罚函数法的一种。增广拉格朗日方法在拉格朗日方法的基础上添加了二次惩罚项，从而使得转换后的问题能够更容易求解，不至于因条件数变大不好求。代码如下：

from scipy.optimize import minimize
import numpy as np 
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt 
'''
增广拉格朗日法
'''
#目标函数：
def func(args):
    fun = lambda x: (x[0]-2)**4 + (x[0]-2*x[1])**2
    return fun
#约束条件，包括等式约束和不等式约束
def con(args):
    cons = ({'type': 'ineq', 'fun': lambda x: 4-x[0]**2-x[1]**2},
    {'type': 'eq', 'fun': lambda x: x[0]**2-x[1]})
    return cons 
#画三维模式图
def draw3D():
    fig = plt.figure()
    ax = Axes3D(fig)
    x_arange = np.arange(-5.0, 5.0)
    y_arange = np.arange(-5.0, 5.0)
    X, Y = np.meshgrid(x_arange, y_arange)
    Z1 = (X-2)**4 - (X-2*Y)**2
    Z2 = 4 - X**2 - Y**2
    Z3 = X**2 - Y
    plt.xlabel('x')
    plt.ylabel('y')
    ax.plot_surface(X, Y, Z1, rstride=1, cstride=1, cmap='rainbow')
    ax.plot_surface(X, Y, Z2, rstride=1, cstride=1, cmap='rainbow')
    ax.plot_surface(X, Y, Z3, rstride=1, cstride=1, cmap='rainbow')
    plt.show()
#画等高线图
def drawContour():
    x_arange = np.linspace(-3.0, 4.0, 256)
    y_arange = np.linspace(-3.0, 4.0, 256)
    X, Y = np.meshgrid(x_arange, y_arange)
    Z1 = (X-2)**4 - (X-2*Y)**2
    Z2 = 4 - X**2 - Y**2
    Z3 = X**2 - Y
    plt.xlabel('x')
    plt.ylabel('y')
    plt.contourf(X, Y, Z1, 8, alpha=0.75, cmap='rainbow')
    plt.contourf(X, Y, Z2, 8, alpha=0.75, cmap='rainbow')
    plt.contourf(X, Y, Z3, 8, alpha=0.75, cmap='rainbow')
    C1 = plt.contour(X, Y, Z1, 8, colors='black')
    C2 = plt.contour(X, Y, Z2, 8, colors='blue')
    C3 = plt.contour(X, Y, Z3, 8, colors='red')
    plt.clabel(C1, inline=1, fontsize=10)
    plt.clabel(C2, inline=1, fontsize=10)
    plt.clabel(C3, inline=1, fontsize=10)
    plt.show()
if __name__ == "__main__":
    args = ()
    args1 = ()
    cons = con(args1)
    x0 = np.array((0, 0)) #设置初始值，初始值的设置很重要，很容易收敛到另外的极值点中，建议多试几个值
    #求解#
    res = minimize(func(args), x0, method='SLSQP', constraints=cons)
    print("最优值：")
    print(res.fun)
    print("最优点：")
    print(res.x)
    draw3D()
    drawContour()

总结

最优化方法（也称做运筹学方法）是近几十年形成的，它主要运用数学方法研究各种系统的优化途径及方案，为决策者提供科学决策的依据。其目的在于针对所研究的系统，求得一个合理运用人力、物力和财力的最佳方案，发挥和提高系统的效能及效益，最终达到系统的最优目标。
无论是约束优化问题还是无约束优化问题，目前都没有一种算法能适用于普遍的场景。不同类型的最优化问题有不同的最优化方法，即使同一类型的问题也可有多种最优化方法。反之,某些最优化方法可适用于不同类型的模型。所以还有很大得改进空间，等待大家去研究。

红帽小生

关注

3
点赞
踩
13

收藏

觉得还不错? 一键收藏
4
评论
优化方法上机实验

优化方法上机实验引言最速下降法牛顿法共轭梯度法拟牛顿法(使用BFGS方法)约束优化方法增广拉格朗日方法总结引言求解无约束非线性规划问题，一般是求目标函数f(x)f(x)f(x)的极小值。至于求极大值问题可以转化为求−f(x)-f(x)−f(x)的极小值问题。一般情况下，不能用解析方法求得准确解，只能用数值方法逐步求其近似解。其思想很简单，从一个初始估计点x(0)x^{(0)}x(0)出发，计算一系列点x(k)x^{(k)}x(k)，这个点列的极限就是f(x)f(x)f(x)的一个极小点。那么由
复制链接

扫一扫

专栏目录