求解双层规划问题的光滑化牛顿方法

A Smoothing Newton Method for Bilevel Programming Problems

摘要: 考虑一类具有特殊结构的双层规划问题，其上层问题无约束，而下层问题为凸优化问题。首先，通过引入下层问题的值函数，将原问题转化为单层优化问题(VP)。在部分平稳性条件下，通过将值函数约束惩罚到目标函数，进一步得到问题(VP)的部分精确罚问题。随后，在一定的约束规范条件下，推导出该罚问题的最优性条件，并利用光滑障碍增广拉格朗日函数方法处理其中的互补约束条件，从而将问题转化为超定方程组进行求解。针对该方程组，采用经典的高斯牛顿法进行求解，并理论证明了当障碍参数趋近于零时算法的收敛性。最后，基于BOLIB库中的小规模问题开展数值实验，将所提算法与LM算法、拟牛顿法等方法的计算结果进行对比分析，实验结果表明了该算法的有效性。

Abstract: Consider a class of bilevel programming problems with a special structure, where the upper-level problem is unconstrained and the lower-level problem is a convex optimization problem. First, by introducing the value function of the lower-level problem, the original problem is reformulated into a single-level optimization problem (VP). Under the partial calmness condition, an exact penalty problem for (VP) is further derived by penalizing the value function constraint into the objective function. Subsequently, under certain constraint qualifications, optimality conditions for the penalized problem are established. The complementary constraints in the optimality conditions are handled using a smoothing barrier augmented Lagrangian function method, transforming the problem into an overdetermined system of equations. To solve this system, the classical Gauss-Newton method is employed, and the convergence of the algorithm is theoretically proven as the barrier parameter tends to zero. Finally, numerical experiments are conducted on small-scale problems from the BOLIB library. The proposed algorithm is compared with the LM algorithm, the Quasi-Newton method, and other approaches. The experimental results demonstrate the effectiveness of the proposed algorithm.

1. 引言

双层规划(Bilevel Programming)是指一类具有两级层次结构的系优化问题，其特点是上层问题和下层问题均具有独立的决策变量，且两层的约束条件与目标函数相互耦合。两层决策者在不同层次上进行决策。上层决策者制定某目标函数并对下层决策者做出一些限制条件，但不直接干涉下层的决策；而下层决策者则在上层给定的约束条件下，优化自身目标函数以确定最优策略。值得注意的是，下层决策会进一步影响上层目标函数的值，从而形成双向交互的双层决策机制。目前，双层规划在管理科学，工程优化及经济决策等领域具有广泛的应用(见[1]-[4])。同时，其理论研究促进了博弈论，最优化以及计算机等研究领域的发展，成为了解决复杂现实问题的重要工具。因此，开发高效的双层规划求解算法具有重要的理论与应用价值。

我们考虑以下双层规划模型：

$\left(\text{BP}\right)\text{ min }F\left(x,y\right)\text{ s}\text{.t}\text{. }y\in S\left(x\right),$ (1)

其中$S\left(x\right)$ 是下层问题的最优解集

$\left({\text{P}}_{\text{x}}\right)\underset{y\in Y\left(x\right)}{\min }f\left(x,y\right),$ (2)

即

$S\left(x\right):=\arg \underset{y}{\min }\left\{f\left(x,y\right)|g\left(x,y\right)\le 0\right\}.$ (3)

这里我们假设$Y\left(x\right)=\left\{y\in {ℝ}^m:g\left(x,y\right)\le 0\right\}$ ，函数$F:{ℝ}^n\times {ℝ}^m\to ℝ$ 是连续可微的，$f:{ℝ}^n\times {ℝ}^m\to ℝ$ ，$g:{ℝ}^n\times {ℝ}^m\to {ℝ}^p$ 是连续可微函数，且对变量$y$ 是二阶连续可微的。定义$ℱ$ 是问题(BP)的可行域。我们

的结果很容易推广到标准的双层规划模型，因此为方便分析结果，我们忽略了上层问题的约束条件。

众所周知，由于双层规划问题复杂的结构使得求解问题(1)较为困难。即使对于最简单的双层规划问题，也已被证明属于NP-hard问题(参见文献[5])，这主要是因为双层规划模型的解通常是近似最优解而非全局最优解。目前有多种方法可以求解双层规划问题(1)，其中最常用的方法就是将其转化为一个单层问题。

一种方法是基于下层问题KKT条件，将双层规划问题转化为带有互补约束的优化问题：

$\begin{array}{l}\left(\text{MPEC}\right)\underset{x,y,\lambda }{\min }F\left(x,y\right)\\ \text{ s}\text{.t}\text{. }{\nabla }_{y}f\left(x,y\right)+{\lambda }^{\text{T}}{\nabla }_{y}g\left(x,y\right)=0,\\ \lambda \ge 0,g\left(x,y\right)\le 0,\\ {\lambda }^{\text{T}}g\left(x,y\right)=0.\end{array}$

该方法为最常用的转化方法之一，上述问题被称为互补约束数学规划(Mathematical Program with Equilibrium Constraint，简称MPEC)。需要指出的是，即便在下层问题为凸的情形下，双层规划问题与对应的MPEC模型也未必等价[6]。为确保二者的等价性，需要添加一些额外的约束条件。Dempe等人[7] [8]研究了当下层问题满足Slater约束规范时，双层规划的全局解(局部解)与MPEC模型全局解(局部解)之间的关系。另一方面，转化之后的问题是非凸的，且理论上MPEC在任何可行点均不满足MFCQ [9]。这会导致MPEC可能存在局部最优解不满足KKT条件的情况，使得一般的非线性规划理论不能直接应用于MPEC问题。为此，学者们引入了一些较弱的稳定点和约束规范(见[10] [11])。此外，在数值求解方面，众多学者提出近似算法求解MPEC模型。[12]中设计了求解MPEC问题的序列二次规划方法，并建立了该方法的强收敛性理论。[13]中提出了求解一般约束优化问题的对偶内点方法，在不依赖正则性条件的假设下，证明了强收敛理论，并结合不等式松弛技术，用该方法求解MPEC模型。[14]中考虑一类上下层问题具有相同变量的问题，研究MPEC模型的最优性理论，通过构造对偶间隙函数将MPEC模型转化为非光滑问题进行研究。关于MPEC理论和算法的最新发展，请读者参阅文献[12]-及其引用的相关研究工作。

另一种常用的方法是利用下层问题的值函数将原问题转化为单层问题。对于任意的$x$ ，假设下层问题的最优解集$S\left(x\right)\ne \varnothing$ ，定义下层问题的值函数为

$V\left(x\right):=\underset{y\in Y\left(x\right)}{\inf }f\left(x,y\right).$ (4)

从值函数的定义易得$y\in S\left(x\right)$ 当且仅当对于$y\in Y\left(x\right)$ ，$f\left(x,y\right)-V\left(x\right)=0$ ，令$\epsilon >0$ ，我们考虑以下单层问题

$\begin{array}{l}\left({\text{VP}}_{\epsilon }\right)\underset{x,y}{\min }F\left(x,y\right)\\ \text{ s}\text{.t}\text{. }g\left(x,y\right)\le 0,\\ f\left(x,y\right)\le V\left(x\right)+\epsilon ,\end{array}$ (5)

若$\epsilon =0$ ，问题$\left({\text{VP}}_{\epsilon }\right)$ 称为问题(VP)，此时为问题(BP)的一种重构；若$\epsilon >0$ ，问题$\left({\text{VP}}_{\epsilon }\right)$ 则构成问题(BP)的近似模型。这一转化问题最初是由Outrata [15]提出用于数值计算，随着被Ye和Zhu [16]将其应用于推导双层规划问题的最优必要性条件。通常情况下，值函数是不可微的，且值函数约束的存在使得该问题很难满足一般的约束规范。针对这一挑战，文献[17]提出了积分熵函数来近似值函数，证明了积分熵函数是值函数的光滑逼近函数，并且在平稳性条件下，设计了一种光滑投影梯度法来求解该问题。除了用光滑函数近似值函数外，将原问题转化为方程组形式也是一种重要的求解方法。文献[18]-[20]中在部分平稳性条件和一定的约束规范下，推导出问题(VP)部分精确罚问题的最优性条件，并用光滑FB函数处理其中的互补约束条件，此时原问题就可转化为一个方程组。针对该系统，学者们运用了高斯牛顿法，半光滑牛顿法，LM算法进行求解，从而近似求解了原问题。

本文中同样基于问题(VP)，在部分平稳性条件和一定的约束规范下，系统推导了该优化问题的最优性条件。针对其中的互补约束条件，我们采用光滑障碍增广拉格朗日函数方法进行处理，从而将其转化为超定方程组。在数值求解方面，文中提出了一种结合线搜索策略的高斯牛顿算法，并证明了当障碍参数趋近于零时算法的收敛性。

第二章系统介绍了相关约束规范的定义以及光滑障碍增广拉格朗日函数方法。第三章建立了在一定条件下将原问题转化为超定方程组的框架以及设计了高斯牛顿法和自适应线搜索法结合的算法，并证明了其收敛性。第四章中基于BOLIB库，进行了小规模的数值实验，并与一些现有典型算法进行了对比分析。

2. 预备知识

2.1. 变分分析

首先介绍本文用到的变分分析[21]中的一些基础知识。给定集合$C\subseteq {ℝ}^d$ ，对于$\forall z\in C$ ，它的正则法锥定义为

${\widehat{N}}_C\left(z\right):=\left\{v\in {ℝ}^d:\langle v,z^{\prime }-z\rangle \le o\left(\Vert z-z^{\prime }\Vert \right)\forall z\in C\right\},$

在该点的极限法锥定义为

$N_C\left(z\right):=\left\{\underset{k\to \infty }{\lim }{v}_k:v_k\in {\widehat{N}}_C\left(z_k\right),z_k\in C,z_k\to z\right\}.$

若集合$C\subseteq {ℝ}^d$ 在$z$ 处局部闭且满足$N_C\left(z\right)={\widehat{N}}_C\left(z\right)$ ，则称该集合在$z$ 处是Clarke正则的。若$f:{ℝ}^d\to \overline{ℝ}$ 是一个下半连续的函数且在$x\in {ℝ}^d$ 处是有限的，可定义该函数在点$x$ 的正则次微分为

$\widehat{\partial }f\left(x\right):=\left\{\zeta \in {ℝ}^d:f\left(x^{\prime }\right)-f\left(x\right)-\langle \zeta ,x-x^{\prime }\rangle \le o\left(\Vert x-x^{\prime }\Vert \right)\right\},$

在该点的极限次微分为

$\partial f\left(x\right):=\left\{\underset{k\to \infty }{\lim }{\xi }_k:{\xi }_k\in \widehat{\partial }f\left(x_k\right),x_k\to x,f\left(x_k\right)\to f\left(x\right)\right\},$

并且函数$f$ 在点$x$ 的Clarke次微分为

${\partial }^{c}f\left(x\right):=co\partial f\left(x\right).$

令$f:{ℝ}^d\to \overline{ℝ}$ 在点$x$ 是Lipschitz连续的，若$\widehat{\partial }f\left(x\right)=\partial f\left(x\right)$ ，我们称函数$f$ 在该点是Clarke正则的。

2.2. 约束规范

若优化问题的目标函数和约束条件均为凸函数，则称该问题为完全凸的。我们考虑以下系统

$\begin{array}{l}\left(\text{P}\right)g_i\left(x\right)\le 0,i=1,\cdots ,p,\\ h_j\left(x\right)=0,j=p+1,\cdots ,q,\end{array}$

其中函数$g_i\left(i=1,\cdots ,p\right)$ 和$h_j\left(j=p+1,\cdots ,q\right)$ 是连续可微的。

对于问题(P)的可行点$\overline{x}$ ，定义$I\left(\overline{x}\right):=\left\{i=1,\cdots ,p,g_i\left(\overline{x}\right)=0\right\}$ 为$\overline{x}$ 处的积极集。下面介绍一些常用的约束规范[21]。

定义2.1. (MFCQ) 在问题(P)的可行点$\overline{x}$ 处满足MFCQ，如果存在向量$d\in {ℝ}^n$ 使得

$\begin{array}{l}\nabla g_i{\left(\overline{x}\right)}^{\text{T}}d<0,\forall i\in I\left(\overline{x}\right),\\ \nabla h_j{\left(\overline{x}\right)}^{\text{T}}d=0,j\in \left\{p+1,\cdots ,q\right\}\end{array}$

成立，并且等式约束梯度集$\left\{\nabla h_j\left(\overline{x}\right),j\in \left\{p+1,\cdots ,q\right\}\right\}$ 是线性无关的。

2.3. 光滑障碍增广拉格朗日函数方法

考虑以下只有不等式约束的非线性优化问题

$\begin{array}{l}\underset{x,y}{\min }f\left(x,y\right)\\ \text{ s}\text{.t}\text{. }g\left(x,y\right)\le 0.\end{array}$

若$g\left(x,y\right)<0$ 严格成立，用文献[22]中提出的对数障碍法重构该问题

$\begin{array}{l}\underset{x,y}{\min }f\left(x,y\right)-r{\displaystyle \sum _{i=1}^p\ln z_i}\\ \text{ s}\text{.t}\text{. }z+g\left(x,y\right)=0.\end{array}$

对于$\rho >0$ ，$s\in {ℝ}^p$ ，我们考虑上述问题的增广拉格朗日函数

$L_r^{\rho }\left(x,y,z,s\right):=f\left(x,y\right)+{\displaystyle \sum _{i=1}^p\left[-r\ln z_i+s_i\left(z_i+g_i\left(x,y\right)\right)+\frac{1}{2\rho }{\left(z_i+g_i\left(x,y\right)\right)}^2\right]}.$

为确保乘子的可行性，我们需要关于$s$ 对上述函数取最大值，进而有以下问题：

$\underset{x,y,z}{\min }\underset{s}{\max }{L}_r^{\rho }\left(x,y,z,s\right),$

对于任意解$\left(x,y,s\right)$ ，一定有${\nabla }_z{L}_r^{\rho }\left(x,y,z,s\right)=0$ ，从而推导出$z$ 是关于$x,y,s,r,\rho$ 的函数：

$z_i\left(x,y,s_i,r,\rho \right):=\frac{1}{2}\left(\sqrt{{\left(\rho s_i+g_i\left(x,y\right)\right)}^2+4r\rho }-\left(\rho s_i+g_i\left(x,y\right)\right)\right),$ (6)

同时定义函数

$\begin{array}{l}{\kappa }_i\left(x,y,s_i,r,\rho \right):=z_i\left(x,y,s_i,r,\rho \right)+g_i\left(x,y\right)+\rho s_i\\ =\frac{1}{2}\left(\sqrt{{\left(\rho s_i+g_i\left(x,y\right)\right)}^2+4r\rho }-\left(\rho s_i+g_i\left(x,y\right)\right)\right),\end{array}$ (7)

其中$i=1,\cdots ,p$ 。为方便叙述，我们将$z_i\left(x,y,s_i,r,\rho \right)$ 和${\kappa }_i\left(x,y,s_i,r,\rho \right)$ 简写为$z_i$ 和${\kappa }_i$ 。

类似于刘等人[23]，我们列出函数$z_i$ 和${\kappa }_i$ 的相关性质。

引理2.1. 对于$r\ge 0$ ，$\rho >0$ ，$i=1,\cdots ,p$ ，下列结论成立。

a) $z_i\ge 0$ ，${\kappa }_i\ge 0$ ，$z_i+g_i\left(x,y\right)={\kappa }_i-\rho s_i$ 且有$z_i{\kappa }_i=r\rho$ ；

b) $z_i+g_i\left(x,y\right)=0$ 当且仅当$g_i\left(x,y\right)\le 0$ ，$s_i\ge 0$ ，$-s_i{g}_i\left(x,y\right)=r$ ；

c) 对于$r>0$ ，$z_i$ 和${\kappa }_i$ 关于变量$x,y,s_i$ 是可微的，且有

${\nabla }_{\left(x,y\right)}{z}_i=\frac{-z_i}{z_i+{\kappa }_i}\nabla g_i\left(x,y\right),{\nabla }_{\left(x,y\right)}{\kappa }_i=\frac{{\kappa }_i}{z_i+{\kappa }_i}\nabla g_i\left(x,y\right),$

${\nabla }_s{z}_i=\frac{-\rho z_i}{z_i+{\kappa }_i}{e}_i,{\nabla }_s{\kappa }_i=\frac{\rho {\kappa }_i}{z_i+{\kappa }_i}{e}_i.$

其中$e_i$ 是第$i$ 个元素为1，其余元素为0的单位向量。

3. 最优性条件

为求得问题(VP)的最优性条件，我们还需引入部分平稳性条件。

定义3.1. (部分平稳性条件) 设$\left(\overline{x},\overline{y}\right)$ 是(VP)的局部最优解，该问题在$\left(\overline{x},\overline{y}\right)$ 满足部分平稳条件，如果存在$\lambda >0$ 及$\left(\overline{x},\overline{y},0\right)$ 的邻域$\mathcal{U}$ 使得

根据文献[24]，如果问题(VP)在局部最优解$\left(\overline{x},\overline{y}\right)$ 满足部分平稳条件，那么存在参数$\lambda >0$ 使得也是以下问题的局部最优解：

$\underset{x,y}{\min }F\left(x,y\right)+\lambda \left(f\left(x,y\right)-V\left(x\right)\right)\text{ s}\text{.t}\text{. }g\left(x,y\right)\le 0.$ (8)

显然，上述问题是将值函数约束作为目标函数的惩罚项。因此，我们通常称它为问题(VP)的部分精确罚问题。

基于上述问题，可应用约束规范推导最优性条件，因此，我们有以下结论。

定理3.1. 令$\left(x,y\right)$ 是问题(VP)的局部最优解。假设该问题中的所有函数均为二阶连续可微的，值函数$V\left(x\right)$ 在$x$ 附近有限且下层问题完全凸。此外，问题(VP)在$\left(x,y\right)$ 处满足部分平稳性条件，下层问题在$\left(x,y\right)$ 处满足MFCQ。那么，存在$\lambda \in \left(0,+\infty \right)$ 和拉格朗日乘子$s$ ，$w$ 使得

${\nabla }_{x}F\left(x,y\right)+\left(w-\lambda s\right){\nabla }_{x}g\left(x,y\right)=0,$ (9)

${\nabla }_{y}F\left(x,y\right)+\lambda {\nabla }_{y}f\left(x,y\right)+w{\nabla }_{y}g\left(x,y\right)=0,$ (10)

${\nabla }_{y}f\left(x,y\right)+s{\nabla }_{y}g\left(x,y\right)=0,$ (11)

$s_i\ge 0,g_i\left(x,y\right)\le 0,s_i{g}_i\left(x,y\right)=0,$ (12)

$w_i\ge 0,g_i\left(x,y\right)\le 0,w_i{g}_i\left(x,y\right)=0,$ (13)

其中$i=1,\cdots ,p$ 。

在上述的最优性条件中，互补约束条件(9)~(13)的存在会使得常见的约束规范不成立从而导致KKT条件无法直接应用。同时，考虑到该系统的可微性，基于引理2.1 (2)，我们用光滑函数(6)处理互补约束条件。对于$\rho >0$ 和充分小的$r>0$ ，我们构造如下光滑系统

${\psi }_{\left(r,\rho ,\lambda \right)}\left(x,y,s,w\right)=\left(\begin{array}{c}{\nabla }_{x}F\left(x,y\right)+\left(w-\lambda s\right){\nabla }_{x}g\left(x,y\right)\\ {\nabla }_{y}F\left(x,y\right)+\lambda {\nabla }_{y}f(x,y)+w{\nabla }_{y}g\left(x,y\right)\\ {\nabla }_{y}f\left(x,y\right)+s{\nabla }_{y}g\left(x,y\right)\\ \begin{array}{l}z\left(x,y,s,r,\rho \right)+g\left(x,y\right)\\ z\left(x,y,w,r,\rho \right)+g\left(x,y\right)\end{array}\end{array}\right)=0.$ (14)

当$r=0$ ，该系统与原系统(9)~(13)等价。根据文献[25]中提出的光滑化方法，对于固定参数$\lambda$ ，我们需要构造一个递减到0的序列$\left\{r_k\right\}$ 来近似求解系统(9)~(13)，从而使得当$k\to \infty$ ，

${\psi }_{\left(r_k,\rho ,\lambda \right)}\left(x_k,y_k,s_k,w_k\right)\to 0,$

显然，式(14)包含$n+m+2p$ 个变量，而方程组个数为$n+2m+2p$ ，这表明该系统是一个超定方程组。接下来，我们采用经典的高斯牛顿法求解该方程组：

算法1

步骤0：给定$\left(x_0,y_0,s_0,w_0\right)$ ，$\lambda >0$ ，${\rho }_0>0$ ，$\epsilon >0$ ，$K>0$ ，$r,\upsilon ,\omega \in \left(0,1\right)$ ，令$k:=0$ 。

步骤1：若$\Vert {\psi }_{\left(r,{\rho }_k,\lambda \right)}\left(x_k,y_k,s_k,w_k\right)\Vert \le \epsilon$ 或$k\ge K$ ，算法结束。

步骤2：计算雅可比矩阵$\nabla {\psi }_{\left(r,{\rho }_k,\lambda \right)}\left(x_k,y_k,s_k,w_k\right)$ 和迭代方向$d_k$

$d_k=-{\left(\nabla {\psi }_{\left(r,{\rho }_k,\lambda \right)}^{\text{T}}\nabla {\psi }_{\left(r,{\rho }_k,\lambda \right)}\right)}^{-1}\nabla {\psi }_{\left(r,{\rho }_k,\lambda \right)}^{\text{T}}{\psi }_{\left(r,{\rho }_k,\lambda \right)}.$ (15)

步骤3：在$\left\{1,\nu ,{\nu }^2,\cdots \right\}$ 中选择最大值作为迭代步长${\alpha }_k\in \left(0,1\right]$ ，并且该步长需要满足不等式

$\begin{array}{l}{\Vert {\psi }_{\left(r,{\rho }_k,\lambda \right)}\left(x_k+{\alpha }_k{d}_{xk},y_k+{\alpha }_k{d}_{yk},s_k+{\alpha }_k{d}_{sk},w_k+{\alpha }_k{d}_{wk}\right)\Vert }^2\\ \le {\Vert {\psi }_{\left(r,{\rho }_k,\lambda \right)}\left(x_k,y_k,s_k,w_k\right)\Vert }^2+\omega {\alpha }_k\nabla {\psi }_{\left(r,{\rho }_k,\lambda \right)}^{\text{T}}\left(x_k,y_k,s_k,w_k\right){\psi }_{\left(r,{\rho }_k,\lambda \right)}\left(x_k,y_k,s_k,w_k\right)d_k.\end{array}$ (16)

步骤4：更新$x_{k+1}=x_k+{\alpha }_k{d}_{xk}$ ，$y_{k+1}=y_k+{\alpha }_k{d}_{yk}$ ，$s_{k+1}=s_k+{\alpha }_k{d}_{sk}$ ，$w_{k+1}=w_k+{\alpha }_k{d}_{wk}$ ，${\rho }_k\ge {\rho }_{k+1}$ 。

步骤5：令$k:=k+1$ ，返回步骤1。

在算法1中，$\epsilon$ 表示容差，K是最大迭代次数。对于步长${\alpha }_k$ 的选取，我们采用Armijo方法。由式(15)可知，要使该算法良定，矩阵$\nabla {\psi }_{\left(r,{\rho }_k,\lambda \right)}^{\text{T}}\nabla {\psi }_{\left(r,{\rho }_k,\lambda \right)}$ 必须是非奇异的。又因为$\nabla {\psi }_{\left(r,\rho ,\lambda \right)}$ 的行数大于列数，该矩阵$\nabla {\psi }_{\left(r,\rho ,\lambda \right)}^{\text{T}}\nabla {\psi }_{\left(r,\rho ,\lambda \right)}$ 的列满秩可保证的非奇异性。接下来，我们将给出确保这一条件成立的可行条件。

为方便叙述，我们定义上下层问题的Lagrange函数

$L^{\lambda }\left(x,y,s,w\right):=F\left(x,y\right)+\left(w-\lambda s\right)g\left(x,y\right),{ℒ}^{\lambda }\left(x,y,s\right):=f\left(x,y\right)+sg\left(x,y\right),$

以及这两个函数关于变量$\left(x,y\right)$ 的Hessian矩阵

${\nabla }^2{L}^{\lambda }:=\left[\begin{array}{cc}{\nabla }_{xx}^2{L}^{\lambda }& {\nabla }_{yx}^2{L}^{\lambda }\\ {\nabla }_{xy}^2{L}^{\lambda }& {\nabla }_{\text{yy}}^2{L}^{\lambda }\end{array}\right],\nabla \left({\nabla }_y{ℒ}^{\lambda }\right):=\left[\begin{array}{cc}{\nabla }_{xy}^2{ℒ}^{\lambda }& {\nabla }_{yy}^2{ℒ}^{\lambda }\end{array}\right].$

此外，令$\nabla g{\left(x,y\right)}^{\text{T}}:=\left[\begin{array}{c}{\nabla }_{x}g{\left(x,y\right)}^{\text{T}}\\ {\nabla }_{y}g{\left(x,y\right)}^{\text{T}}\end{array}\right]$ ，进而计算${\psi }_{\left(r,\rho ,\lambda \right)}$ 的雅可比矩阵：

$\nabla {\psi }_{\left(r,\rho ,\lambda \right)}\left(x,y,s,w\right)=\left(\begin{array}{ccc}{\nabla }^2{L}^{\lambda }& -\lambda \nabla g\left(x,y\right)& \nabla g\left(x,y\right)\\ \nabla \left({\nabla }_yℒ\right)& {\nabla }_{y}g\left(x,y\right)& 0\\ \mathcal{J}\nabla g{\left(x,y\right)}^{\text{T}}& \Gamma & 0\\ \Theta \nabla g{\left(x,y\right)}^{\text{T}}& 0& \mathfrak{B}\end{array}\right),$

其中$\mathcal{J}:=diag\left\{{\tau }_1,\cdots ,{\tau }_p\right\}$ ，$\Gamma :=diag\left\{{\gamma }_1,\cdots ,{\gamma }_p\right\}$ ，$\Theta :=diag\left\{{\theta }_1,\cdots ,{\theta }_p\right\}$ ，$\mathfrak{B}:=diag\left\{{\beta }_1,\cdots ,{\beta }_p\right\}$ ，这些矩阵的元素分别定义为

${\tau }_i=\frac{{\kappa }_i}{z_i+{\kappa }_i},{\gamma }_i=\frac{-\rho z_i}{z_i+{\kappa }_i},{\theta }_i=\frac{{\tilde{\kappa }}_i}{{\tilde{z}}_i+{\tilde{\kappa }}_i},{\beta }_i=\frac{-\rho {\tilde{z}}_i}{{\tilde{z}}_i+{\tilde{\kappa }}_i}\text{, }i=1,\cdots ,p.$ (17)

其中$z_i=z_i\left(x,y,s_i,r,\rho \right)$ ，${\tilde{z}}_i={\tilde{z}}_i\left(x,y,w_i,r,\rho \right)$ ，${\kappa }_i$ 和${\tilde{\kappa }}_i$ 同理定义。

类似于文献[25]，我们首先讨论上述元素的相关性质。

引理3.1. 对于每一$r>0$ ，$\rho >0$ 和使得${\psi }_{\left(r,\rho ,\lambda \right)}\left(x,y,s,w\right)=0$ 成立的点$\left(x,y,s,w\right)$ ，以下结论成立：

$\begin{array}{l}{\tau }_i>0,{\gamma }_i<0,i=1,\cdots ,p,\\ {\theta }_i>0,{\beta }_i<0,i=1,\cdots ,p.\end{array}$

Proof. 根据引理2.1 (1)，对于任意$r>0$ ，$\rho >0$ ，该引理显然成立。 ■

基于上述引理，我们可以推导出保证矩阵$\nabla {\psi }_{\left(r,{\rho }_k,\lambda \right)}^{\text{T}}\nabla {\psi }_{\left(r,{\rho }_k,\lambda \right)}$ 非奇异的充分条件。这将确保高斯牛顿迭

代格式(15)是良定的。如同文献[25]所述，我们只需证明是列满秩的即可。

定理3.2. 对于某些$r>0$ ，$\rho >0$ ，$0<\lambda <\frac{{\gamma }_i}{{\tau }_i}\frac{{\theta }_i}{{\beta }_i}$ ，$i=1,\cdots ,p$ 和满足(14)的点$\left(\overline{x},\overline{y},\overline{s},\overline{w}\right)$ ，假设${\nabla }^2{L}^{\lambda }\left(\overline{x},\overline{y},\overline{s},\overline{w}\right)$ 是正定的。那么，$\nabla {\psi }_{\left(r,\rho ,\lambda \right)}\left(\overline{x},\overline{y},\overline{s},\overline{w}\right)$ 的列向量是线性无关的。

Proof. 考虑任意向量$d={\left(d_1^{\text{T}},d_2^{\text{T}},d_3^{\text{T}}\right)}^{\text{T}}$ 使得$\nabla {\psi }_{\left(r,\rho ,\lambda \right)}\left(\overline{x},\overline{y},\overline{s},\overline{w}\right)d=0$ ，其中$d_1\in {ℝ}^{n+m}$ ，$d_2\in {ℝ}^p$ ，$d_3\in {ℝ}^p$ 。因此我们有

${\nabla }^2{L}^{\lambda }\left(\overline{x},\overline{y},\overline{s},\overline{w}\right)d_1-\lambda \nabla g{\left(\overline{x},\overline{y}\right)}^{\text{T}}{d}_2+\nabla g{\left(\overline{x},\overline{y}\right)}^{\text{T}}{d}_3=0,$ (18)

$\nabla \left({\nabla }_{{}^y}{ℒ}^{\lambda }\left(\overline{x},\overline{y},\overline{s},\overline{w}\right)\right)d_1+{\nabla }_{y}g{\left(\overline{x},\overline{y}\right)}^{\text{T}}{d}_2=0,$ (19)

${\tau }_i\nabla g_i{\left(\overline{x},\overline{y}\right)}^{\text{T}}{d}_1+{\gamma }_i{d}_{2i}=0,$ (20)

${\theta }_i\nabla g_i{\left(\overline{x},\overline{y}\right)}^{\text{T}}{d}_1+{\beta }_i{d}_{3i}=0,$ (21)

其中$i=1,\cdots ,p$ 。根据引理3.1，我们可将式(20)和式(21)写为

$\nabla g_i{\left(\overline{x},\overline{y}\right)}^{\text{T}}{d}_1=-\frac{{\gamma }_i}{{\tau }_i}{d}_{2i},\nabla g_i{\left(\overline{x},\overline{y}\right)}^{\text{T}}{d}_1=-\frac{{\beta }_i}{{\theta }_i}{d}_{3i},$(22)

其中$-\frac{{\gamma }_i}{{\tau }_i}>0$ ，$-\frac{{\beta }_i}{{\theta }_i}>0$ 。因此，我们对式(18)左乘$d_1^{\text{T}}$ ：

$d_1^{\text{T}}{\nabla }^2{L}^{\lambda }\left(\overline{x},\overline{y},\overline{s},\overline{w}\right)d_1-\lambda d_1^{\text{T}}\nabla g{\left(\overline{x},\overline{y}\right)}^{\text{T}}{d}_2+d_1^{\text{T}}\nabla g{\left(\overline{x},\overline{y}\right)}^{\text{T}}{d}_3=0.$ (23)

将式(22)带入上式有

$d_1^{\text{T}}{\nabla }^2{L}^{\lambda }\left(\overline{x},\overline{y},\overline{s},\overline{w}\right)d_1-\lambda {\displaystyle \sum _{i=1}^p\left(-\frac{{\gamma }_i}{{\tau }_i}\right)d_{2i}^2}+{\displaystyle \sum _{i=1}^p\left(-\frac{{\beta }_i}{{\theta }_i}\right)d_{3i}^2}=0.$ (24)

更多地，由式(22)得

$d_{2i}=-\frac{{\tau }_i}{{\gamma }_i}\nabla g_i{\left(\overline{x},\overline{y}\right)}^{\text{T}}{d}_1=\frac{{\tau }_i}{{\gamma }_i}\frac{{\beta }_i}{{\theta }_i}{d}_{3i},$ (25)

并将该式带入式(24)有

$d_1^{\text{T}}{\nabla }^2{L}^{\lambda }\left(\overline{x},\overline{y},\overline{s},\overline{w}\right)d_1+{\displaystyle \sum _{i=1}^p\left(1-\lambda \frac{{\tau }_i}{{\gamma }_i}\frac{{\beta }_i}{{\theta }_i}\right)}\left(-\frac{{\beta }_i}{{\theta }_i}\right)d_{3i}^2=0.$(26)

那么，根据引理3.1，由于${\nabla }^2{L}^{\lambda }\left(\overline{x},\overline{y},\overline{s},\overline{w}\right)$ 是正定的，$0<\lambda <\frac{{\gamma }_i}{{\tau }_i}\frac{{\theta }_i}{{\beta }_i}$ ，$i=1,\cdots ,p$ ，式(26)是非负项的求和，则$d_1$ 和$d_3$ 只能是零向量。根据式(25)，进而推出$d_{2i}=0$ ，$i=1,\cdots ,p$ 。证毕。 ■

需要注意的是，$\lambda <\frac{{\gamma }_i}{{\tau }_i}\frac{{\theta }_i}{{\beta }_i}$ 这一假设并不与$\lambda$ 严格为正相冲突，因为根据引理3.1，有$\frac{{\gamma }_i}{{\tau }_i}\frac{{\theta }_i}{{\beta }_i}>0$ 。

基于上述定理，我们还需做出以下假设，以证明$\nabla {\psi }_{\left(r,\rho ,\lambda \right)}^{\text{T}}\nabla {\psi }_{\left(r,\rho ,\lambda \right)}$ 的非奇异性。

假设3.1. 以下矩阵的每一行都是非零向量：

$\left[{\nabla }^2{L}^{\lambda }{\left(x,y,s,w\right)}^{\text{T}}\nabla {\left({\nabla }_yℒ\left(x,y,s,w\right)\right)}^{\text{T}}\nabla g{\left(x,y\right)}^{\text{T}}\right].$

假设3.2. 对于$\lambda >0$ ，$r>0$ ，$\rho >0$ ，矩阵$\nabla {\psi }_{\left(r,\rho ,\lambda \right)}^{\text{T}}\nabla {\psi }_{\left(r,\rho ,\lambda \right)}$ 是严格对角占优的。

引理3.2. 令假设3.1在点$\left(x,y,s,w\right)$ 成立，那么对于任意的$\lambda >0$ ，$r>0$ ，$\rho >0$ ，矩阵$\nabla {\psi }_{\left(r,\rho ,\lambda \right)}^{\text{T}}\nabla {\psi }_{\left(r,\rho ,\lambda \right)}$ 的对角线元素是严格正的。

Proof. 定义$a_{ij}$ ，$i,j=1,2,3$ 为矩阵的组成元素，那么该矩阵写为

$\nabla {\psi }_{\left(r,\rho ,\lambda \right)}^{\text{T}}\nabla {\psi }_{\left(r,\rho ,\lambda \right)}=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right].$

显然，我们只需考虑元素$a_{11}$ ，$a_{22}$ ，$a_{33}$ 我们将这三个元素的表达式列出：

$\begin{array}{l}{\left(a_{11}\right)}_{jj}={\displaystyle \sum _{k=1}^{n+m}{\left({\nabla }_{j,k}^2{L}^{\lambda }\right)}^{\text{T}}}{\nabla }_{j,k}^2{L}^{\lambda }+{\displaystyle \sum _{k=1}^m{\nabla }_j{\left({\nabla }_{y_k}ℒ\right)}^{\text{T}}}{\nabla }_j\left({\nabla }_{y_k}ℒ\right)\\ \text{ +}{\displaystyle \sum _{k=1}^p{\nabla }_j{g}_k{\left(x,y\right)}^{\text{T}}{\nabla }_j{g}_k\left(x,y\right){\left({\tau }_k\right)}^2}+{\displaystyle \sum _{k=1}^p{\nabla }_j{g}_k{\left(x,y\right)}^{\text{T}}{\nabla }_j{g}_k\left(x,y\right){\left({\beta }_k\right)}^2,}\end{array}$

其中${\nabla }_j$ 是$\nabla :=\left({\nabla }_{x_1},\cdots ,{\nabla }_{x_n},{\nabla }_{y_1},\cdots ,{\nabla }_{y_m}\right)$ 的第$j$ 个元素和${\nabla }_{j,k}^2$ 是以下矩阵的第$j$ 行第$k$ 列元素：

$\left(\begin{array}{cccccc}{\nabla }_{x_1{x}_1}& \cdots & {\nabla }_{x_1{x}_n}& {\nabla }_{x_1{y}_1}& \cdots & {\nabla }_{x_1{y}_m}\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ {\nabla }_{x_n{x}_1}& \cdots & {\nabla }_{x_n{x}_n}& {\nabla }_{x_n{y}_1}& \cdots & {\nabla }_{x_n{y}_m}\\ {\nabla }_{y_1{x}_1}& \cdots & {\nabla }_{y_1{x}_n}& {\nabla }_{y_1{y}_1}& \cdots & {\nabla }_{y_1{y}_m}\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ {\nabla }_{y_m{x}_1}& \cdots & {\nabla }_{y_m{x}_n}& {\nabla }_{y_m{y}_1}& \cdots & {\nabla }_{y_m{y}_m}\end{array}\right).$

结合假设3.1和引理3.1，易得${\left(a_{11}\right)}_{jj}>0$ ，$j=1,\cdots ,n+m$ 。类似地，其他对角线元素为

$\begin{array}{l}{\left(a_{22}\right)}_{ii}=\nabla g_i\left(x,y\right)\nabla g_i{\left(x,y\right)}^{\text{T}}+{\left({\gamma }_i\right)}^2,\\ {\left(a_{33}\right)}_{ii}=\nabla g_i\left(x,y\right)\nabla g_i{\left(x,y\right)}^{\text{T}}+{\left({\beta }_i\right)}^2,i=1,\cdots ,p.\end{array}$

根据引理4.1，易得这两个元素也均大于0。证毕。 ■

定理3.3. 对于某些$\lambda >0$ ，$r>0$ ，$\rho >0$ ，令$\left(x,y,s,w\right)$ 是满足系统(14)的稳定点。若假设3.1和假设3.2在该点成立，则$\nabla {\psi }_{\left(r,\rho ,\lambda \right)}^{\text{T}}\nabla {\psi }_{\left(r,\rho ,\lambda \right)}$ 是非奇异的。

Proof. 该结果可参考文献[13]中定理4.8类似证明。 ■

为证明算法1的收敛性，我们还需证明步长的存在性。令$\left(x_k,y_k,s_k,w_k\right)$ 是第k次迭代点，则我们有以下结果。

引理3.3. 对于$\lambda >0$ ，$r>0$ ，$\rho >0$ ，假设函数$F$ ，$f$ 和$g$ 是二阶连续可微的，那么存在步长${\alpha }_k\in \left(0,1\right]$ 使得式(16)成立。

定义残差函数${\varphi }_{\left(r,\rho ,\lambda \right)}\left(x,y,s,w\right)=\frac{1}{2}{\Vert {\psi }_{\left(r,\rho ,\lambda \right)}\left(x,y,s,w\right)\Vert }^2$ ，下列结论成立。

引理3.4. 令$z_{k+1}\left(\rho \right)=z\left(x_{k+1},y_{k+1},s_{k+1},r,\rho \right)$ ，${\tilde{z}}_{k+1}\left(\rho \right)=z\left(x_{k+1},y_{k+1},w_{k+1},r,\rho \right)$ 。假设 $\Vert z_{k+1}\left({\rho }_k\right)+g\left(x_{k+1},y_{k+1}\right)\Vert \ne 0$ ，$\Vert {\tilde{z}}_{k+1}\left({\rho }_k\right)+g\left(x_{k+1},y_{k+1}\right)\Vert \ne 0$ 。若${\rho }_k\ge {\rho }_{k+1}>0$ 则有

${\varphi }_{\left(r,{\rho }_{{}_{k+1}},\lambda \right)}\left(x_{k+1},y_{k+1},s_{k+1},w_{k+1}\right)\le {\varphi }_{\left(r,{\rho }_{{}_k},\lambda \right)}\left(x_k,y_k,s_k,w_k\right).$

Proof. 根据$z$ 和$\tilde{z}$ 表达式和引理2.1，易计算出

$\frac{\partial z_i}{\partial \rho }=\frac{1}{\rho }\frac{z_i}{z_i+{\kappa }_i}\left(z_i+g_i\right),\frac{\partial {\tilde{z}}_i}{\partial \rho }=\frac{1}{\rho }\frac{{\tilde{z}}_i}{{\tilde{z}}_i+{\tilde{\kappa }}_i}\left({\tilde{z}}_i+g_i\right),$

其中$i=1,\cdots ,p$ ，进一步可得

$\begin{array}{l}{\frac{D{\varphi }_{\left(r,\rho ,\lambda \right)}\left(x_{k+1},y_{k+1},s_{k+1},w_{k+1}\right)}{D\rho }|}_{\rho ={\rho }_k}\\ =\frac{1}{{\rho }_k}{\left({\widehat{z}}_{k+1}+g\left(x_{k+1},y_{k+1}\right)\right)}^{\text{T}}{\left(Z_{k+1}+K_{k+1}\right)}^{-1}{Z}_{k+1}\left({\widehat{z}}_{k+1}+g\left(x_{k+1},y_{k+1}\right)\right)\\ +\frac{1}{{\rho }_k}{\left({\tilde{z}}_{k+1}+g\left(x_{k+1},y_{k+1}\right)\right)}^{\text{T}}{\left({\tilde{Z}}_{k+1}+{\tilde{K}}_{k+1}\right)}^{-1}{Z}_{k+1}\left({\tilde{z}}_{k+1}+g\left(x_{k+1},y_{k+1}\right)\right)\\ >0,\end{array}$

其中${\widehat{z}}_{k+1}=z\left(x_{k+1},y_{k+1},s_{k+1},r,{\rho }_k\right)$ ，${\tilde{z}}_{k+1}=z\left(x_{k+1},y_{k+1},w_{k+1},r,{\rho }_k\right)$ ，$Z_{k+1}=diag\left({\widehat{z}}_{k+1}\right)$ ，$K_{k+1}=diag\left({\widehat{\kappa }}_{k+1}\right)$ ，${\tilde{Z}}_{k+1}=diag\left({\tilde{z}}_{k+1}\right)$ ，${\tilde{K}}_{k+1}=diag\left({\tilde{\kappa }}_{k+1}\right)$ 。因此，根据上式，可推出当$\rho >0$ 时，${\varphi }_{\left(r,\rho ,\lambda \right)}\left(x_{k+1},y_{k+1},s_{k+1},w_{k+1}\right)$ 关于$\rho$ 是单调递增的，从而得出最终结论。证毕。 ■

以下是算法1收敛性定理。

定理3.4. 若定理3.1中的假设和假设3.1~3.2在$\left(x_k,y_k,s_k,w_k\right)$ 处成立，那么

$\underset{k\to \infty }{\lim }{\varphi }_{\left(r,{\rho }_{{}_k},\lambda \right)}\left(x_k,y_k,s_k,w_k\right)=0.$

备注 在上述假设下，我们提出的光滑化的高斯牛顿法在理论上表现出良好的收敛性。但在实际生活中，这些关键假设可能难以完全满足，从而对算法的适用性带来挑战。接下来，我们对这些假设的局限性以及潜在影响进行系统分析：关于惩罚参数$\lambda$ 的后验条件对于算法的收敛以及后续数值实验的结果具有重要的影响，若惩罚参数过大，可能会导致惩罚过度从而数值不稳定，若惩罚参数过小，则可能无法保证约束从而导致解不可行。严格对角占优假设的成立保证了高斯牛顿法中雅可比矩阵的非奇异性，从而保证了算法的良定性，但在实际问题中该条件可能过于严格，在高维或非凸问题中可能难以验证，所以我们数值实验考虑解决小规模问题，若该条件不成立，可能导致迭代矩阵病态，进而导致算法收敛到非稳定点，对于该问题实际中可通过正则化方法解决，但这会增加算法的复杂度。

4. 数值实验

本节中数值实验聚焦于求解系统(14)。基于该系统，我们比较了高斯牛顿法[25]，LM (Levenberg-Marquardt)算法[25]，信赖域算法[26]，拟牛顿法[27]以及Nelder-Mead单纯形法的测试结果。实验所用案例均来自双层优化测试问题库(BOLIB) [28]，该库包含124个非线性测试案例。本实验是在配置Intel(R) Core(TM) Ultra 5 125H 3.60 GHz处理器，Windows11系统的笔记本电脑上进行，实验平台为Matlab2023b。

在算法1的第0步及其他算法中，我们将容差设为$\epsilon =1\text{e}-6$ ，最大迭代次数设为$K=1000$ 。所有实验均采用六个不同的惩罚参数值，即$\lambda \in \left\{0.01,0.1,1,10,100,1000\right\}$ 。采取不同的$\lambda$ 值的原因是为了避免对下层目标值与最优值的过度惩罚或惩罚不足。本节后续部分将展示各算法在CPU时间，相对误差以及残差等方面的性能比较。默认设置下，我们以$x_0=1_n$ ，$y_0=1_m$ 作为初始点。若默认起始点不满足至少一个约束条件时，则选择可行起始点(详见文献[28])。对于拉格朗日乘子，按以下方式进行初始化：$s_0=w_0=1_p$ 。

4.1. 时间性能曲线

性能曲线图[29]被广泛用于比较不同算法的特性。本小节我们用性能曲线图分析各算法的CPU时间，$t_{i,j}$ 表示算法$j$ 求解问题$i$ 所需的平均时间。若问题$i$ 的最优解已知但无法被算法$j$ 求解(即上层目标函数误差 > 60%)，设定$t_{i,j}=\infty$ 。定义时间性能比率：

$r_{i,j}:=\frac{t_{i,j}}{\min \left\{t_{i,j}:j\in J\right\}},$

其中$J$ 表示算法集合。性能比是指算法$j$ 在解决问题$i$ 时的时间与集合$J$ 中表现最佳的算法所用时间的比值。定义算法$j$ 性能比的累积分布函数：

${\rho }_j\left(\tau \right):=\frac{\left|\left\{i\in I:r_{i,j}\le \tau \right\}\right|}{n_i},$

其中$I$ 表示问题集合。性能剖面图${\rho }_j\left(\tau \right)$ 用于统计算法$j$ 的性能比优于$\tau$ 的问题数量。该性能剖面函数${\rho }_j:ℝ\to \left[0,1\right]$ 关于$\tau$ 是非递减函数，其中${\rho }_j\left(1\right)$ 表示算法$j$ 表现最优的问题数与全部问题的比例。

Figure 1. CPU time performance curves of different algorithms

图1. 不同算法的CPU时间性能曲线图

上图1为不同算法性能比的折线图，曲线越高表示算法性能越优。纵轴表示算法CPU时间性能比率小于$\tau$ 的问题占比。由图1可得，除单纯形法，高斯牛顿法的表现弱于其他算法。根据${\rho }_j\left(1\right)$ 的数值可知：高斯牛顿法在约10%的问题最快，LM算法约15%的问题最快，拟牛顿法约22%的问题最快，而信赖域算法约在60%的问题中最快。图中还显示信赖域算法的${\rho }_j\left(2\right)$ 就已经达到了90%，而其他算法最多的也就只有约60%，这说明在CPU时间指标上，信赖域算法的性能明显优于其他算法。但是我们的实验针对的是小规模问题，还需通过以下表1进一步分析：

Table 1. Average execution time of algorithms

表1. 各算法的平均运行时间

算法	高斯牛顿法	LM	信赖域	拟牛顿法	单纯形法
平均时间(秒)	6.68E−02	1.61E−02	8.16E−03	1.04E−02	4.02E+00

上述表格，高斯牛顿法，LM算法和拟牛顿法在时间上量级相同，说明这三种算法的时间性能相当。除此之外，无论从性能曲线还是平均时间上，单纯形法显著差于其他算法，而信赖域算法优于其他算法。

4.2. 目标函数的相对误差

本小节我们将比较不同算法计算所得的上下层目标函数值。为进行此项比较，我们仅关注116个BOLIB中的案例，因为有部分算法在其他7个案例中未求出明确的解，例如高斯牛顿法在NieEtal2017e案例中可能因矩阵的奇异性而导致发散(详见文献中[30])。

定义$F$ ，$f$ 为通过某算法求得的上下层目标函数值，$F_k$ ，$f_k$ 为文献[28]中BOLIB库中记录的最优解对应的最优值，以及相对误差的公式定义为

$uppererror=\frac{\left|F-F_k\right|}{1+\left|F_k\right|},lowererror=\frac{\left|f-f_k\right|}{1+\left|f_k\right|}.$

以下图像分别展示了不同参数下不同算法的相对误差对比结果，图像按误差升序排列。

Figure 2. Comparison chart of the relative errors in the upper-level and lower-level objective functions across algorithms under different parameters

图2. 不同参数下各算法上下层目标函数相对误差对比图

从图2中可以看出，对于约80%的问题，高斯牛顿法产生的误差小于其他算法，说明该算法比其他算法表现好，尤其是在$\lambda =100$ 时明显能看出高斯牛顿法的上层目标函数的误差曲线低于其他曲线。具体而言，在约50%测试问题中，高斯牛顿法所得解产生的误差可忽略不计(误差小于5%)。在6%~30%的误差范围内，当罚参数分别为0.01，0.1，1时，除单纯性法外，该算法与其他算法差别不大，还原了约70%的解。在其他参数时，高斯牛顿法还原了近60%的解，而其他算法表现最佳的也仅还原了约50%的解。

4.3. 停止准则中的残差变化

在本小节中我们将评估算法1在放宽停止条件中容差时的表现。具体而言，在实际的数值实验中，当容差为1e−8时，多数案例无法实现该精度。因此，当算法达到最大迭代次数或步长过小时，算法终止。将不同算法在不同参数条件下终止时的最小残差画图如下：

Figure 3. Residual performance comparison of algorithms

图3. 各算法残差性能对比图

上图3显示的是各算法最终生成的KKT残差数值按升序排列结果。具体来看，从图3可得高斯牛顿法在95个案例中表现较好，在这些案例中残差值均小于1e−8。这表明在这些案例中，当问题可被近乎求到精确解时，该算法的停止条件更为严格，从而能得到较优的数值解。整体来看，高斯牛顿法明显优于其他算法，在同一容差下，其他算法仅能求解约50%~60%的问题，因此，该算法更易达到稳定点。

5. 结论

本文中，基于对数障碍法与增广拉格朗日函数，构造了下层问题近似值函数，进而求得问题(VP)的近似问题。基于光滑化方法，我们将该近似问题重新表述为一个非线性方程组。理论分析表明，本文中所采用的高斯牛顿法具有良定性，并在一定条件下可证明该算法的收敛性。在数值实验，我们分别从时间，相对误差和残差三方面对比了我们的算法与LM算法，信赖域算法以及拟牛顿法的结果。实验结果表明，在时间性能方面，当$\nabla {\psi }_{\left(r,\rho ,\lambda \right)}^{\text{T}}\nabla {\psi }_{\left(r,\rho ,\lambda \right)}$ 在迭代中不存在病态条件时，除单纯形法外，四种算法均展现出极快的运算速度，在不同的罚参数下平均运行时间均低于0.1秒，但根据时间性能曲线，我们的算法在运算速度上慢于其他三种算法，这可以作为我们今后的改进工作；在相对误差方面，该模型适用于多数双层规划问题，高斯牛顿法在不同参数下可还原60%~70%的问题的已知最优解，尤其当惩罚参数为100时，我们的算法的精度明显比其他算法要高，但是我们的实验只列举了六种参数下的结果，对于惩罚参数$\lambda$ 的选取规则一直是双层规划中的一大难题，后续我们将其作为一个重要的研究方向；在残差方面，高斯牛顿法明显优于其他算法，这表明我们的算法在求解双层规划问题上相较于几种经典算法具有更高的精度。综上，在BOLIB测试库的124个非线性案例中，高斯牛顿法展现出显著优势：不仅成功求解率最高，其残差下降速度也明显优于其他对比算法。