ABSTRACT

Nonlinear programming has always been a hot topic in the research of optimization theory. In this paper, three commonly used methods for solving unconstrained nonlinear programming (gradient method, Newton method and quasi-Newton method) are compared and analyzed. Two new algorithms are given by improving the original method, and the convergence of the improved algorithm is illustrated. On this basis, this paper makes an empirical study on the location of distribution center based on competition.

Keywords:Gradient, Convergence, Location Problem

基于牛顿法及拟牛顿法的非线性规划算法改进及实证研究

付巍巍^1,2，张立卫^1*，董玉¹

¹大连理工大学数学科学学院，辽宁 大连

²辽宁工程技术大学理学院，辽宁 阜新

收稿日期：2020年7月15日；录用日期：2020年7月28日；发布日期：2020年8月4日

摘 要

非线性规划一直是优化理论研究的热点问题。本文将求解无约束非线性规划的三种常用方法(梯度法，牛顿法和拟牛顿法)进行对比分析，改进原有方法从而给出两种新算法，并说明改进算法的收敛性。在此基础上对基于竞争的配送中心选址问题进行了实证研究。

关键词 :梯度，收敛，选址问题

Copyright © 2020 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

非线性规划是最优化问题中的重要一类。一般的最优化问题可以表达成一个二元组 $\left(D,f\right)$，D为任一集合，是问题的可行区域；f为目标函数，是由集合D到一维欧式空间 $R^1$ 的一个映射。问题是求出D中的一个元素 $X^{*}$，使对 $\forall X^{*}\in D$，有 $f\left(X^{*}\right)\le f\left(X\right)$，此时 $X^{*}$ 称为全局最优解。

对于无约束非线性规划问题：

$\min f(\; X\; )$

其中 $X={\left(x_1,x_2,\cdots ,x_n\right)}^{\text{T}}$ 是n维空间中的向量点。

2. 梯度法

梯度法是1847年法国著名数学家Cauchy最早提出 [1]。梯度法的基本原理：若目标函数有一阶连续偏导数，在 $X^{*}$ 取得极小点，通过在迭代点 $X^{\left(k\right)}$ 负梯度方向上的一维搜索，来确定使得 $f\left(X^{\left(k+1\right)}\right)=f\left(X^{\left(k\right)}+t_k{d}^{\left(k\right)}\right)$ 最小的 $t_k$ (即最佳步长)，得到逐渐趋近于极小点的序列 $\left\{X^{\left(k\right)}\right\}$，这种梯度法就是所谓的最速下降法。

梯度法方法简单，每迭代一次的工作量较小，所需要的存储量也小，即使从一个不好的初始点出发，也能保证算法的收敛性。而缺点是在极小点附近收敛的很慢。由于梯度是函数的局部性质，从局部看，在一点附近下降或上升的快，但从总体看可能走许多弯路。事实上，若 $X^{\left(k+1\right)}$ 是沿方向 $d^{\left(k\right)}=-\nabla f\left(X^{\left(k\right)}\right)$ 用一维搜索求得的，即 $t_k$ 为 $\phi \left(t\right)=f\left(X^{\left(k\right)}+t{d}^{\left(k\right)}\right)$ 的极小点，那么 ${\phi }^{\prime }\left(t_k\right)=0$，即

$\nabla f{\left(X^{\left(k+1\right)}\right)}^{\text{T}}\cdot \nabla f\left(X^{\left(k\right)}\right)=0$

说明相邻的两个迭代方向正交，迭代路线呈锯齿状，尤其在最优解附近，锯齿现象尤为严重，从而影响迭代速度，如图1。另外，梯度法产生的序列是线性收敛的，而且收敛性质与最优解处的Hessian阵的特征值关系很大，且关于小的扰动是不稳定的，任何小的扰动都可能破坏算法的收敛性。

图1. 梯度法迭代路线图

例1. 用梯度法求

$\min f\left(X\right)=-2x_1{x}_2-2x_2+x_1^2+2x_2^2$

初始点 $X^{\left(0\right)}={\left(0,0\right)}^{\text{T}}$，$\epsilon =0.1$。

图2. 例1中目标函数的函数图像

解：图2为目标函数的函数图像，从图上可以看出二次函数的精确最优解 $X^{\ast }={\left(1,1\right)}^{\text{T}}$，目标函数的

梯度为 $\nabla f\left(X\right)=-{\left(2x_2-2x_1,2x_1+2-4x_2\right)}^{\text{T}}$，所以 $\nabla f\left(X^{\left(0\right)}\right)={\left(0,-2\right)}^{\text{T}}$，$\Vert \nabla f\left(X^{\left(0\right)}\right)\Vert =2>0.1$。

第一次迭代：

从 $X^{\left(0\right)}={\left(0,0\right)}^{\text{T}}$ 出发，沿方向 $d^{\left(0\right)}=-\nabla f\left(X^{\left(0\right)}\right)=\left(\begin{array}{c}0\\ 2\end{array}\right)$ 进行一维搜索， $t_0=\frac{1}{4}$，在直线上的极大点

$X^{\left(1\right)}=X^{\left(0\right)}+t_0{d}^{\left(0\right)}=\left(\begin{array}{c}0\\ 0\end{array}\right)+\frac{1}{4}\left(\begin{array}{c}0\\ 2\end{array}\right)=\left(\begin{array}{c}0\\ \frac{1}{2}\end{array}\right)$

如此重复进行迭代，会如图所示，得到点列 $\left(\frac{1}{2}，\frac{3}{4}\right),\left(\frac{3}{4},\frac{3}{4}\right),\left(\frac{3}{4},\frac{7}{8}\right),\cdots$ 逐渐趋近于问题的最优解 $X^{*}=\left(\begin{array}{c}1\\ 1\end{array}\right)$，可以验证 $\nabla f\left(X^{*}\right)=0$。具体迭代情况如图3及表1所示。

图3. 梯度法求解例1时迭代点及迭代路线图

表1. 梯度法求解例1的初值点和迭代次数表

3. 牛顿法

牛顿法最初由艾萨克·牛顿于1736年在 Method of Fluxions中公开提出 [2]。若非线性目标函数 $f\left(X\right)$ 具有二阶连续偏导数，又设 $X^{\left(k\right)}$ 为其极值点的某一近似，在这一点做 $f\left(X\right)$ 的二阶泰勒展开，即：

$f\left(X\right)\approx g\left(X\right)=f\left(X^{\left(k\right)}\right)+\nabla f{\left(X^{\left(k\right)}\right)}^{\text{T}}\left(X-X^{\left(k\right)}\right)+\frac{1}{2}{\left(X-X^{\left(k\right)}\right)}^{\text{T}}{\nabla }^{2}f\left(X^{\left(k\right)}\right)\left(X-X^{(\; k\; )}\right)$

其中 ${\nabla }^{2}f\left(X^{\left(k\right)}\right)$ 为 $f\left(X\right)$ 在 $X^{\left(k\right)}$ 处的Hessian矩阵。

为求 $g\left(X\right)$ 的平稳点，令 $\nabla g\left(X\right)=0$，即

$\nabla f\left(X^{\left(k\right)}\right)+{\nabla }^{2}f\left(X^{\left(k\right)}\right)\left(X-X^{\left(k\right)}\right)=0$

设矩阵 ${\nabla }^{2}f\left(X^{\left(k\right)}\right)$ 可逆，则得到牛顿法的算法迭代公式：

$X^{\left(k+1\right)}=X^{\left(k\right)}-{\left[{\nabla }^{2}f\left(X^{\left(k\right)}\right)\right]}^{-1}\nabla f\left(X^{(\; k\; )}\right)$

如果 $f\left(X\right)$ 是二次函数，则其Hessian矩阵 ${\nabla }^{2}f\left(X\right)$ 为常数阵。在这种情况下，从任意一点出发，只要一步即可求出 $f\left(X\right)$ 的极值点(假设Hessian矩阵正定)。

如果 $f\left(X\right)$ 不是二次函数，二阶泰勒展式仅是一个近似表达式。此时，求得的极值点只是 $f\left(X\right)$ 的近似极值点。在这种情况下，选取搜索方向：

$d^{\left(k\right)}=-{\left[{\nabla }^{2}f\left(X^{\left(k\right)}\right)\right]}^{-1}\nabla f\left(X^{(\; k\; )}\right)$

$X^{\left(k+1\right)}=X^{\left(k\right)}+d^{(\; k\; )}$

按照这种方式求函数 $f\left(X\right)$ 极小点的方法称为牛顿法，称方向 $d^{\left(k\right)}=-{\left[{\nabla }^{2}f\left(X^{\left(k\right)}\right)\right]}^{-1}\nabla f\left(X^{\left(k\right)}\right)$ 为牛顿

方向。牛顿法收敛的速度很快，当 $f\left(X\right)$ 的二阶导数及Hessian矩阵的逆矩阵便于计算时，这一方法非常有效。

例2用牛顿法求解例1问题

$\min f\left(X\right)=-2x_1{x}_2-2x_2+x_1^2+2x_2^2$

初始点 $X^{\left(0\right)}={\left(0,0\right)}^{\text{T}}$，$\epsilon =0.1$。

解：目标函数 $f\left(X\right)$ 的梯度为 $\nabla f\left(X\right)={\left(-2x_2+2x_1,-2x_1-2+4x_2\right)}^{\text{T}}$。

Hessian阵 $H\left(X\right)={\nabla }^{2}f\left(X\right)=\left(\begin{array}{cc}2& -2\\ -2& 4\end{array}\right)$ 为正定矩阵(即目标函数为二次凸函数)。

初始点 $X^{\left(0\right)}={\left(0,0\right)}^{\text{T}}$，于是有 $\nabla f\left(X^{\left(0\right)}\right)={\left(0,-2\right)}^{\text{T}}$，此时 $\Vert \nabla f\left(X^{\left(0\right)}\right)\Vert =2>0.1$，开始第一次迭代

$H\left(X^{\left(0\right)}\right)={\nabla }^{2}f\left(X^{\left(0\right)}\right)=\left(\begin{array}{cc}2& -2\\ -2& 4\end{array}\right)$，$H{\left(X^{\left(0\right)}\right)}^{-1}=\left(\begin{array}{cc}1& \frac{1}{2}\\ \frac{1}{2}& \frac{1}{2}\end{array}\right)$

根据牛顿法的迭代公式有：

$X^{\left(1\right)}=X^{\left(0\right)}-H{\left(X^{\left(0\right)}\right)}^{-1}\nabla f\left(X^{\left(0\right)}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)-\left(\begin{array}{cc}1& \frac{1}{2}\\ \frac{1}{2}& \frac{1}{2}\end{array}\right)\left(\begin{array}{c}0\\ -2\end{array}\right)=(\; 1\; 1\; )$

于是 $\nabla f\left(X^{\left(1\right)}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)$，所以 $X^{*}=X^{\left(1\right)}=\left(\begin{array}{c}1\\ 1\end{array}\right)$ 为问题精确的极小点，极小值为 $f\left(X^{\left(1\right)}\right)=-1$。

显然，与例1对比可以知道，牛顿法的收敛速度要远远大于梯度法。

牛顿法每一步的迭代步长都是 $t_k=1$，这样虽然算法简单，但是会在一定的程度上影响收敛速度，所以每次迭代时，我们依然可以采用一维搜索的方法计算出最优步长，即找到搜索方向上的最优迭代点，该种方法被称为阻尼牛顿法。

选择非线性规划常用的测试函数 [3]： $\min f=100{\left(x_2-x_1^2\right)}^2+{\left(1-x_1\right)}^2$，该函数图像如图4所示，此问题解析解为(1, 1)，最优值为0，设误差为0.1进行迭代求解，结果见表2。

由下面的分析可以知道，梯度法对于初值的依赖性不强(初值点从 $\left(0,0\right)\to \left(10,10\right)\to \left(20,20\right)\to \left(100,100\right)$，在满足精度要求的前提下，迭代次数从 $156\to 18374\to 57561\to 619505$ 次)，越接近极值点收敛速度越慢，而恰好牛顿法对于初值的依赖性很强(初值点从 $\left(0,0\right)\to \left(10,10\right)\to \left(20,20\right)\to \left(40,40\right)$，迭代次数从 $3\to 12598\to 88854\to 167330$ 次，初值点在(0, 0)时，只需迭代3次即可满足精度要求，而初值点在(10, 10)时，却需要迭代12598次才能满足精度要求)，在极值点附近收敛速度很快。再对比牛顿法和阻尼牛顿法，显然可以看出阻尼牛顿法由于每一次迭代都在搜索方向上寻求了最优的迭代点，所以他的收敛速度相比牛顿法会有明显的提升，但牛顿法固有的对初值依赖性强的问题并没有得到解决，所以我们想到结合两种算法，首先使用梯度法进行迭代，当迭代点进入极值点的某个充分小邻域时，再使用阻尼牛顿法进行搜索，从而提升问题的收敛速度，本文将此种改进算法称作GN算法。

图4. 测试函数的函数图像

表2. 三种方法求解测试函数的初值点和迭代次数表

GN算法的算法步骤：

(1) 给出初始点 $X^{\left(0\right)}\in E^n$，精度，令 $k=0$ ；

(2) 计算 $\nabla f\left(X^{\left(k\right)}\right)$ ；

(3) 若 $\Vert \nabla f\left(X^{\left(k\right)}\right)\Vert <{\epsilon }_1$，则迭代结束，取 $\bar{X}=X^{\left(k\right)}$，否则转(4)；

(4) 若 $\Vert \nabla f\left(X^{\left(k\right)}\right)\Vert \ge {\epsilon }_1$。

(I)若 $\Vert \nabla f\left(X^{\left(k\right)}\right)\Vert \ge {\epsilon }_2$，用一维搜索求 $\phi \left(t\right)=f\left(X^{\left(k\right)}+t{d}^{\left(k\right)}\right)$ 的一个极小点 $t_k$，使 $f\left(X^{\left(k\right)}+t_k{d}^{\left(k\right)}\right)<f\left(X^{\left(k\right)}\right)$，其中 $d^{\left(k\right)}=-\nabla f\left(X^{\left(k\right)}\right)$。

令 $X^{\left(k+1\right)}=X^{\left(k\right)}+t_k{d}^{\left(k\right)},k=k+1$，转(2)。

(II)若 $\Vert \nabla f\left(X^{\left(k\right)}\right)\Vert <{\epsilon }_2$，用一维搜索求 $\phi \left(t\right)=f\left(X^{\left(k\right)}+t{d}^{\left(k\right)}\right)$ 的一个极小点 $t_k$，使 $f\left(X^{\left(k\right)}+t_k{d}^{\left(k\right)}\right)<f\left(X^{\left(k\right)}\right)$，其中 $d^{\left(k\right)}=-{\left[{\nabla }^{2}f\left(X^{\left(k\right)}\right)\right]}^{-1}\cdot \nabla f\left(X^{\left(k\right)}\right)$。

令 $X^{\left(k+1\right)}=X^{\left(k\right)}+t_k{d}^{\left(k\right)},k=k+1$，转(2)。

表3. GN算法求解测试函数的初值点及迭代次数表

由表3可以看出，当初值点位于(100, 100)时，G-N算法经18,000次迭代，即可得到非常精确(误差仅为3.9211e−14)的最优解，而此时无论是牛顿法还是阻尼牛顿法甚至都是不能得到近似解的，因为牛顿法每一次迭代都严格要求迭代点处的Hessian矩阵可逆，而满足这点并不容易。我们改进后的算法只在极值点的充分小的邻域内使用阻尼牛顿法，既发挥了牛顿法迭代速度快的优点，又最大限度的避免了初值对其影响大的缺点，在局部最优解的充分小的邻域内依然可以保证至少二阶的收敛速度，是一种有效的改进算法。

4. 拟牛顿法

牛顿法，它的突出优点是收敛很快，但是运用牛顿法需要计算二阶偏导数，而目标函数的Hessian矩阵可能非正定。由于实际问题的目标函数往往相当复杂，计算二阶导数矩阵及其逆矩阵相当困难或根本不可能，因此，确定牛顿方向存在很大的障碍。为避免计算二阶导数矩阵 ${\nabla }^{2}f\left(X^{\left(k\right)}\right)$ 及其逆矩阵 ${\left[{\nabla }^{2}f\left(X^{\left(k\right)}\right)\right]}^{-1}$，提出了拟牛顿法。

拟牛顿法的基本思想：使用不含二阶导数的矩阵 $H_k$ 近似替代牛顿法中的Hessian阵的逆矩阵 ${\left[{\nabla }^{2}f\left(X^{\left(k\right)}\right)\right]}^{-1}$，构造 $H_k$ 应遵循如下三条原则：

(1) 每一步迭代均能按已有的信息确定下一个搜索方向；

(2) 每一步迭代均能使目标函数有所改进；

(3) 近似矩阵 $H_k$ 最终应收敛于极值点处的Hessian矩阵的逆矩阵。

由于构造近似矩阵的方法不同，因而出现不同的拟牛顿法。拟牛顿法不是一种算法，而是一族公认的比较有效的算法。

前面已经给出了牛顿法的迭代方向，即牛顿方向为：

$d^{\left(k\right)}=-{\left[{\nabla }^{2}f\left(X^{\left(k\right)}\right)\right]}^{-1}\nabla f\left(X^{(\; k\; )}\right)$

构造 ${\left[{\nabla }^{2}f\left(X^{\left(k\right)}\right)\right]}^{-1}$ 的近似矩阵 $H_k$，得 $d^{\left(k\right)}=-H_k\nabla f\left(X^{\left(k\right)}\right)$，称之为点 $X^{\left(k\right)}$ 处的拟牛顿方向。

设在第k次迭代后，得到点 $X^{\left(k+1\right)}$，将目标函数 $f\left(X\right)$ 在 $X^{\left(k+1\right)}$ 点处展开成Taylor级数，并取二阶近似，记 $p^{\left(k\right)}=X^{\left(k+1\right)}-X^{\left(k\right)}$，$q^{\left(k\right)}=\nabla f\left(X^{\left(k+1\right)}\right)-\nabla f\left(X^{\left(k\right)}\right)$。

为了用不包含二阶导数的矩阵 $H_{k+1}$ 取代牛顿法中的Hessian矩阵的逆矩阵 ${\left[{\nabla }^{2}f\left(X^{\left(k\right)}\right)\right]}^{-1}$，有理由令 $H_{k+1}$ 满足：

$p^{\left(k\right)}=H_{k+1}\cdot q^{(\; k\; )}$

这个条件称为拟牛顿条件。

显然，若想使拟牛顿方向 $d^{\left(k\right)}=-H_k\nabla f\left(X^{\left(k\right)}\right)$ 较好的近似于牛顿方向，下面基本条件应当得到满足：

(1) $d^{\left(k\right)}=-H_k\nabla f\left(X^{\left(k\right)}\right)$ 是目标函数 $f\left(X\right)$ 在点 $X^{\left(k\right)}$ 的下降方向；

(2) $H_{k+1}$ 满足 $p^{\left(k\right)}=H_{k+1}\cdot q^{\left(k\right)}$ ；

(3) $H_{k+1}$ 与之间应具有某种简单的迭代关系

$H_{k+1}=H_k+E_k$

其中 $E_k$ 是一个容易计算的矩阵，称之为校正矩阵。即 $p^{\left(k\right)}=\left(H_k+E_k\right)\cdot q^{\left(k\right)}$。

(4) 若有 $\nabla f{\left(X^{\left(k\right)}\right)}^{\text{T}}\left[-H_k\nabla f\left(X^{\left(k\right)}\right)\right]>0$，即 $\nabla f{\left(X^{\left(k\right)}\right)}^{\text{T}}{H}_k\nabla f\left(X^{\left(k\right)}\right)<0$ 成立，则拟牛顿方向 $d^{\left(k\right)}=-H_k\nabla f\left(X^{\left(k\right)}\right)$ 必是下降方向。

构造校正矩阵 $E_k$ 的方法有很多，每一种特殊的构造方法对应一种具体的算法，因而拟牛顿法是一族算法，本文选择有代表性的变尺度算法，其主要的优点是：既避免了计算二阶导数矩阵及其求逆过程，又比梯度法收敛的速度快，特别是对高维问题具有显著的优越性。变尺度法最早由Davidon于1959年提出，后经Fletcher和Powell二人改进，因此变尺度法也被称为DFP法 [2]。

变尺度法中，校正矩阵 $E_k$ 的构造方式为：

$E_k={\alpha }_k{u}^{\left(k\right)}{u}^{\left(k\right)\text{T}}+{\beta }_k{v}^{\left(k\right)}{v}^{\left(k\right)\text{T}}$

其中 ${\alpha }_k,{\beta }_k$ 是待定的数，是待定的向量。

将上式代入拟牛顿条件 $p^{\left(k\right)}=\left(H_k+E_k\right)\cdot q^{\left(k\right)}$ 中，可以得到：

$\begin{array}{l}{\alpha }_k=\frac{1}{u^{\left(k\right)\text{T}}{q}^{\left(k\right)}}=\frac{1}{p^{\left(k\right)\text{T}}{q}^{\left(k\right)}}\\ {\beta }_k=\frac{-1}{v^{\left(k\right)\text{T}}{q}^{\left(k\right)}}=\frac{-1}{q^{\left(k\right)\text{T}}{H}_k{q}^{(\; k\; )}}\end{array}$

于是

$E_k=\frac{p^{\left(k\right)}{p}^{\left(k\right)\text{T}}}{p^{\left(k\right)\text{T}}{q}^{\left(k\right)}}-\frac{H_k{q}^{\left(k\right)}{q}^{\left(k\right)\text{T}}{H}_k}{q^{\left(k\right)\text{T}}{H}_k{q}^{(\; k\; )}}$

下面通过例3来展示变尺度法的计算过程。

例3用变尺度法(DFP)求例1

$\min f\left(X\right)=-2x_1{x}_2-2x_2+x_1^2+2x_2^2$

初始点 $X^{\left(1\right)}={\left(0,0\right)}^{\text{T}}$，初始矩阵 $H_1=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$

解：由题可知 $f\left(X\right)=-2x_1{x}_2-2x_2+x_1^2+2x_2^2$

$\nabla f\left(X\right)=\left(\begin{array}{c}2x_1-2x_2\\ -2x_1-2+4x_2\end{array}\right)$

在点 $X^{\left(1\right)}={\left(0,0\right)}^{\text{T}}$ 处的梯度 $g_1=\left(\begin{array}{c}0\\ -2\end{array}\right)$，$\Vert g_1\Vert =2\ne 0$。

第一次迭代：

令搜索方向

$d^{\left(1\right)}=-H_1{g}_1=-\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}0\\ -2\end{array}\right)=(\; 0\; 2\; )$

从 $X^{\left(1\right)}={\left(0,0\right)}^{\text{T}}$ 出发，沿方向 $d^{\left(1\right)}$ 进行一维搜索：

$\underset{t\ge 0}{\min }f\left(X^{\left(1\right)}+t{d}^{(\; 1\; )}\right)$

求得 $t_1=\frac{1}{4}$，令

$X^{\left(2\right)}=X^{\left(1\right)}+t_1{d}^{\left(1\right)}=\left(\begin{array}{c}0\\ 0\end{array}\right)+\frac{1}{4}\left(\begin{array}{c}0\\ 2\end{array}\right)=\left(\begin{array}{c}0\\ \frac{1}{2}\end{array}\right)$

此时 $g_2=\nabla f\left(X^{\left(2\right)}\right)=\left(\begin{array}{c}-1\\ 0\end{array}\right)$，$\Vert g_2\Vert =1\ne 0$

第二次迭代：

令 $p^{\left(1\right)}=X^{\left(2\right)}-X^{\left(1\right)}=t_1{d}^{\left(1\right)}=\left(\begin{array}{c}0\\ \frac{1}{2}\end{array}\right)$

$q^{\left(1\right)}=g_2-g_1=\left(\begin{array}{c}-1\\ 0\end{array}\right)-\left(\begin{array}{c}0\\ -2\end{array}\right)=\left(\begin{array}{c}-1\\ 2\end{array}\right)$

于是

$H_2=H_1+E_1=H_1+\frac{p^{\left(1\right)}{p}^{\left(1\right)\text{T}}}{p^{\left(1\right)\text{T}}{q}^{\left(1\right)}}-\frac{H_1{q}^{\left(1\right)}{q}^{\left(1\right)\text{T}}{H}_1}{q^{\left(1\right)\text{T}}{H}_1{q}^{\left(1\right)}}=\frac{1}{20}\left(\begin{array}{cc}16& 8\\ 8& 9\end{array}\right)$

令

$d^{\left(2\right)}=-H_2{g}_2=-\frac{1}{20}\left(\begin{array}{cc}16& 8\\ 8& 9\end{array}\right)\left(\begin{array}{c}-1\\ 0\end{array}\right)=\frac{1}{5}(\; 4\; 2\; )$

从 $X^{\left(2\right)}={\left(0,\frac{1}{2}\right)}^{\text{T}}$ 出发，沿方向 $d^{\left(2\right)}$ 进行一维搜索：

$\underset{t\ge 0}{\min }f\left(X^{\left(2\right)}+t{d}^{(\; 2\; )}\right)$

求得 $t_2=\frac{5}{4}$，令

$X^{\left(3\right)}=X^{\left(2\right)}+t_2{d}^{\left(2\right)}=(\; 1\; 1\; )$

此时 $g_3=\nabla f\left(X^{\left(3\right)}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)$。

所以 $X^{*}=\left(\begin{array}{c}1\\ 1\end{array}\right)$ 为此问题的最优解。

在例3的讨论中，第一个尺度矩阵 $H_1$ 为单位矩阵(对称正定阵)，以后的尺度矩阵 $H_k$ 按算法步骤逐步形成。可以证明，如此构造的尺度矩阵均为对称正定阵，进而确保搜索方向为下降方向。

DFP方法由于不需要在迭代点处计算Hessin矩阵的逆，所以在极值点附近的收敛速度比牛顿法更快，效率更高。可由表4列出的数值实验可以看出，当初始点距离极值点很远时，受迭代误差的影响，DFP方法很难收敛到极值点。因此，为了提升收敛速度和求解的准确性我们同样可以在迭代点距离极值点很近时再使用拟牛顿方法求解问题。本文将此种改进方法称作GNN算法。

表4. DFP算法求解测试函数的初值点及迭代次数表

由表4和表5，明显可以看出GNN算法的精确度和收敛速度对比DFP方法均有提升，是一种行之有效的优化算法。

表5. GNN算法求解测试函数的初值点及迭代次数表

5. 物流配送中心选址问题

假设在某区域内有8个需求点，需求点位置坐标如表6所示，选址区域D内已有2个配送中心，配送中心坐标分别为 $\left(x_j,y_j\right)$，(4, 4)，(6, 5)，现欲在该区域内选择一个配送中心，问该如何选择配送中心的位置才能使总收益达到最大?

表6. 需求点信息表

假设新建配送中心的位置在 $\left(x,y\right)$ 处，根据Huff引力模型 [4]，在选址区域D内，已存在多个配送中心，新建一个配送中心，采用配送中心利润最大原则进行，则问题模型为求新建配送中心利润的最大值 [5]。

$\underset{\left(x,y\right)\in D}{\max }F\left(x,y\right)=\omega {\displaystyle \underset{i=1}{\overset{n}{\sum }}{p}_i\left(x,y\right)\cdot v_i}=\omega {\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{1}{1+{\displaystyle \underset{j=1}{\overset{m}{\sum }}{\text{e}}^{-\lambda \left(d_{ij}-d_i\right)}}}\cdot v_i}$，其中表示新建配送中心 $\left(x,y\right)$ 分担需

求点i的需求量的比例， $\left(a_i,b_i\right),i=1,2,\cdots ,8$ 表示8个需求点的位置坐标， $v_i$ 表示需求点i的需求量， $\omega$ 表示配送单位需求量的利润， $\lambda$ 为距离调节系数， $d_i$ 表示新建配送中心到需求点i的距离， $d_{ij}$ 表示已有配送中心j到需求点i的距离，即

$\begin{array}{l}{d}_{ij}=\sqrt{{\left(a_i-x_j\right)}^2+{\left(b_i-y_j\right)}^2};\\ d_i=\sqrt{{\left(x-a_i\right)}^2+{\left(y-b_i\right)}^2}\end{array}$

取 $\omega =1,\lambda =0.5$，于是该问题可表示为如下非线性规划 [6]：

$\underset{\left(x,y\right)\in D}{\max }F\left(x,y\right)={\displaystyle \underset{i=1}{\overset{8}{\sum }}\frac{1}{1+{\displaystyle \underset{j=1}{\overset{2}{\sum }}{\text{e}}^{-\frac{1}{2}\left(d_{ij}-d_{{}_i}\right)}}}{v}_i}$，

该问题目标函数如图5所示，若利用Nelder-Mead simplex direct search方法可求出该问题的最优解为 $x^{*}=\left(6.9201,3.1403\right)$，在该处建立配送中心可获得利润的最大值为：157.6933。

图5. 物流配送中心问题目标函数图像

接下来，我们分别使用梯度法、牛顿法、拟牛顿方法以及改进的两种算法来求解此问题。

表7. 五种算法对比分析表

由表7的对比分析可以看出，五种算法中效率最差的是梯度法，一般而言梯度法只具有线性收敛速度，而阻尼牛顿法和DFP算法在极值点的充分小邻域内具有至少二阶的收敛速度。但是对于本问题而言，在极值点附近阻尼牛顿法的效率明显高于DFP算法，是因为该问题的运算量不大，各迭代点处Hessian矩阵的逆存在，阻尼牛顿法在有限步内迭代简单可逐渐收敛到极值点，而DFP算法在此邻域内因为要叠加误差，所以收敛的速度反而没有阻尼牛顿法快，但是当迭代点到达(18, 18)时，受Hessian矩阵的逆的存在性影响，阻尼牛顿法已经不能收敛到极值点，但此问题显然借助DFP算法可以解决。而毫无疑问改进的两种算法无论从迭代速度还是初值点的选择上都增加了问题的实用性，对初值点的选择没有那么强，同时迭代速度对比DFP算法有明显提升。

6. 结论

本文结合梯度法和阻尼牛顿法提出了GN算法，结合梯度法和变尺度法提出了GNN算法，从理论上解释算法的收敛性，用测试函数进行了收敛速度的验证。同时将几种算法应用于有竞争机制的物流中心选址问题，验证算法的实用性，同时给出算法的对比分析。

文章引用

付巍巍,张立卫,董 玉. 基于牛顿法及拟牛顿法的非线性规划算法改进及实证研究
The Improvement and Empirical Study of Nonlinear Programming Algorithm Based on Newton Method and Quasi-Newton Method[J]. 应用数学进展, 2020, 09(08): 1146-1158. https://doi.org/10.12677/AAM.2020.98134

参考文献