求解线性互补约束问题的一个SQP方法

A SQP Method for Linear Complementary Constraints Problem

摘要: 对线性互补约束优化问题，利用一个连续可微的光滑互补函数，光滑系数趋于零时，原问题转化为光滑非线性规划问题，再利用SQP算法来求解该光滑非线性规划问题。在适当的条件下，该算法具有全局收敛性，最后给出了数值实验，也证明了该算法是有效的。

Abstract: In this paper, equilibrium problem with linear equilibrium constrains is studied. By using a smoothing complimentarily function which is differentiable everywhere, when smoothing parameter tends to zero, the original problem is equivalently transformed to a smoothing nonlinear optimization, then the smooth nonlinear optimization is solved by SQP algorithm. The algorithm converges globally under certain conditions. Finally, numerical experiments are given to prove the effectiveness of the algorithm.

1. 引言

近年来，带均衡约束的数学规划问题因为在经济、工程技术、对策决策等领域的广泛应用而引起关注 [1] [2] 。本文求解一个线性互补均衡约束问题，可以描写为以下形式：

$\begin{array}{l}\min f\left(x,y\right)\\ \text{s}\text{.t}\text{.}Ax\le b,\\ w=Nx+My+q,\\ 0\le w\perp y\ge 0,\end{array}$ (1.1)

其中 $f:R^{n+m}\to R$ ，是连续可微的函数， $A\in R^{p\times n},N\in R^{m\times n},M\in R^{m\times m},b\in R^p,q\in R^m$ 。

如果将互补条件 $w\perp y$ 写成内积的形式： $w^{\text{T}}y=0$ ，那么均衡约束问题等价于一个光滑的非线性规划问题。但众所周知，对于此类非线性划归问题，即使较弱的MFCQ规格,在任何可行点处都不成立，故传统的优化算法一般不能直接用来求解此类问题。

对线性互补均衡约束问题，先介绍一些常见的光滑函数：

广义Fischer-Burmeister互补函数：

${\varphi }_1\left(a,b\right)=a+b-\sqrt{a^2+b^2+\mu },$

Zang在文献 [3] 提出一个如下形式的光滑函数：

${\varphi }_2\left(a,\mu \right)=\{\begin{array}{l}0,\text{}\text{}\text{}\text{}a\le -\frac{\mu }{2},\\ \frac{1}{2\mu }{\left(a+\frac{\mu }{2}\right)}^2,-\frac{\mu }{2}<a<\frac{\mu }{2},\\ a,a\ge -\frac{\mu }{2}.\end{array}$

Song Xu在文献 [4] 中提出一个如下形式的新的光滑函数：

${\varphi }_3\left(a,\mu \right)=\mu \ln {\displaystyle \underset{i=1}{\overset{m}{\sum }}{\text{e}}^{\frac{a_i}{\mu }}}.$

本文尝试利用互补函数 $\varphi \left(y_i,w_i,\mu \right)=-\mu \ln \left(1+{\text{e}}^{-\frac{y_i}{\mu }}\right)-\mu \ln \left(1+{\text{e}}^{-\frac{w_i}{\mu }}\right)$ ，将问题(1.1)中的互补约束条件

$w\perp y$ 光滑化，将问题(1.1)转化为一光滑非线性规划的一般约束优化问题，然后借助罚函数型SQP方法思想，通过求解后者逼近前者的解，并证明该算法是全局收敛的.

为了方便，本文采用如下记号：

$\begin{array}{l}z=\left(x,y,w\right),s=\left(x,y\right),t=\left(y,w\right),z^k=\left(x^k,y^k,w^k\right),s^k=\left(x^k,y^k\right),t^k=\left(y^k,w^k\right),\\ t_i=\left(y_i,w_i\right),\text{d}z=\left(\text{d}x,\text{d}y,\text{d}w\right),\text{d}{z}^k=\left(\text{d}{x}^k,\text{d}{y}^k,\text{d}{w}^k\right),X=\left\{z\in X_0|0\le w\perp y\ge 0\right\},\\ X_0=\left\{z|Ax\le b,w=Nx+My+q\right\},A^{\text{T}}=\left(a_1^{\text{T}},\cdots ,a_p^{\text{T}}\right),b^{\text{T}}=\left(b_1,\cdots ,b_p\right).\end{array}$

本文主要符号解释如下：

$R^n$ 表示一个实n维Euclide空间。

e 表示各分量均为1的列向量。

$\nabla f\left(x\right)$ 表示函数 $f\left(x\right)$ 的一阶导数或梯度。

${\nabla }^{2}f\left(x\right)$ 表示函数 $f\left(x\right)$ 的二阶导数或 Hessian 矩阵。

$u>0$ 表示向量u的每一个分量 $u_j>0$ 。

$u\ge 0$ 表示向量u的每一个分量 $u_j\ge 0$ 。

$x^{\text{T}}$ 和 $A^{\text{T}}$ 分别表示向量x和矩阵A的转置。

$diag\left(a_j,j=1,2,~,m\right)$ 表示以 $a_j$ 为对角元素的m阶对角矩阵。

2. 问题光滑化

记F是问题(1.1)的可行域，定义在向量 $z\in F$ 处F的切锥如下：

$L\left(z,F\right)=\left\{d=\underset{k\to \infty }{\lim }\frac{z^k-z}{\tau }:z^k\in F,\underset{k\to \infty }{\lim }{z}^k=z,{\tau }_k\downarrow 0\right\},$

对问题(1.1)，需做如下基本假设：

H 2.1 问题(1.1)中的矩阵M是一个 $P_0$ 矩阵，即所有的主子式均非负。

H 2.2 对任意的 $z\in X_0$ ，矩阵A的行向量组 $\left\{a_j|j\in L\left(z\right)\right\}$ 线性无关。

H 2.3 问题(1.1)的可行域非空。

H 2.4 函数 $f\left(x,y\right)$ 为二阶连续可微的。

定义 2.1 如果对于可行点 $z^{\ast }=\left(x^{\ast },y^{\ast },w^{\ast }\right)$ 满足下面的条件，那么称 $z^{\ast }$ 是问题(1.1)的稳定点：

$\text{d}z=\left(\text{d}x,\text{d}y,\text{d}w\right)\in \Gamma \left(X;z^{\ast }\right)\Rightarrow \nabla f{\left(x^{\ast },y^{\ast }\right)}^{\text{T}}\text{d}s\ge 0,$

其中 $\Gamma \left(X;z^{\ast }\right)$ 为X在 $z^{\ast }$ 处的切锥。

问题(1.1)的稳定点与其最优解两者之间存在密切的联系，求解问题(1.1)的算法多半集中在稳定点的求解。

定理 2.1 [5] 假设 $z^{\ast }=\left(x^{\ast },y^{\ast },w^{\ast }\right)\in X$ 满足下层非退化条件：

$\left(y_i,w_i\right)\ne \left(0,0\right),\forall i=1,\cdots ,m,$

则 $z^{\ast }$ 为问题(1.1)的一个稳定点的充分必要条件是存在乘子 $\left({\lambda }^{\ast },u^{\ast },v^{\ast }\right)$ ，使得下面的方程组成立：

$\begin{array}{l}\left(\begin{array}{c}\nabla f_x\left(x^{\ast },y^{\ast }\right)\\ \nabla f_y\left(x^{\ast },y^{\ast }\right)\\ 0\end{array}\right)+\left(\begin{array}{c}{A}^{\text{T}}\\ 0\\ 0\end{array}\right){\lambda }^{\ast }+\left(\begin{array}{c}{N}^{\text{T}}\\ M^{\text{T}}\\ -I\end{array}\right)u^{\ast }+\left(\begin{array}{c}0\\ W^{\ast }\\ Y^{\ast }\end{array}\right)v^{\ast }=0,\\ {\left(A{x}^{\ast }-b\right)}^{\text{T}}{\lambda }^{\ast }=0,A{x}^{\ast }-b\ge 0,{\lambda }^{\ast }\ge 0,\end{array}$ (2.1)

其中 $W^{\ast }=diag\left(w_i^{\ast },i=1,\cdots ,m\right),Y^{\ast }=diag\left(y_i^{\ast },i=1,\cdots ,m\right)$ 。

本文为问题(1.1)提出一个新的光滑函数.

令 $\phi \left(y,w\right)=\min \left\{y,w\right\}$ ，其中 $\min \left\{y,w\right\}$ 定义为 $\min \left\{y_i,w_i\right\}e,i=1,\cdots ,m$ ， $e={\left(1,\cdots ,1\right)}^{\text{T}}\in R^m$ 。

有

$y^{\text{T}}w=0,y\ge 0,w\ge 0\iff \phi \left(y,w\right)=\min \left\{y,w\right\}=0.$ (2.2)

定义函数 $\Phi \in R^{2m}\times \left[0,+\infty \right)$ ：

$\Phi \left(a,b,\mu \right)=\left(\begin{array}{c}\varphi \left(a_1,b_1,\mu \right)\\ ⋮\\ \varphi \left(a_m,b_m,\mu \right)\end{array}\right),$ (2.3)

其中 $\varphi \in R^{2m}\times \left[0,+\infty \right)$ ， $\mu$ 是一个非负的参数，

$\varphi \left(a,b,\mu \right)=-\mu \ln \left(1+{\text{e}}^{-\frac{a}{\mu }}\right)-\mu \ln \left(1+{\text{e}}^{-\frac{b}{\mu }}\right),$ (2.4)

易知对于固定的 $\mu >0$ ，函数 $\varphi \left(\cdot ,\cdot ,\mu \right)$ 处处可微.

$\begin{array}{l}\nabla \varphi \left(y,w,\mu \right)\\ =\left(\begin{array}{l}{\nabla }_y\varphi \left(y,w,\mu \right)\\ \nabla w\varphi \left(y,w,\mu \right)\end{array}\right)\\ ={\left(\frac{\partial \varphi \left(y_1,w_1,\mu \right)}{\partial y_1},\cdots ,\frac{\partial \varphi \left(y_m,w_m,\mu \right)}{\partial y_m},\frac{\partial \varphi \left(y_1,w_1,\mu \right)}{\partial w_1},\cdots ,\frac{\partial \varphi \left(y_m,w_m,\mu \right)}{\partial w_m}\right)}^{\text{T}}\\ ={\left(\frac{{\text{e}}^{-\frac{y_1}{\mu }}}{1+{\text{e}}^{-\frac{y_1}{\mu }}},\cdots ,\frac{{\text{e}}^{-\frac{y_m}{\mu }}}{1+{\text{e}}^{-\frac{y_m}{\mu }}},\frac{{\text{e}}^{-\frac{w_1}{\mu }}}{1+{\text{e}}^{-\frac{w_1}{\mu }}},\cdots ,\frac{{\text{e}}^{-\frac{w_m}{\mu }}}{1+{\text{e}}^{-\frac{w_m}{\mu }}}\right)}^{\text{T}}.\end{array}$ (2.5)

$\varphi \left(a,b,\mu \right)$ 有一些很好的性质：

$\underset{\mu \to 0}{\lim }\varphi \left(a,b,\mu \right)=\underset{\mu \to 0}{\lim }-\mu \ln \left(1+{\text{e}}^{-\frac{a}{\mu }}\right)-\mu \ln \left(1+{\text{e}}^{-\frac{b}{\mu }}\right)=0,$

且当 $\mu =0$ 时，

$\underset{\mu \to 0}{\lim }\Phi \left(y,w,\mu \right)=\phi \left(y,w\right)=0.$ (2.6)

由上式，做如下定义：

$\Phi \left(y,w,0\right)=\phi \left(y,w\right),$ (2.7)

从而有

$\underset{\mu \to 0}{\lim }\varphi \left(y,w,\mu \right)=\varphi \left(y,w,0\right).$ (2.8)

如将问题(1.1)的互补条件转化为线性方程组 $\varphi \left(y,w,0\right)=0,i=1,\cdots ,m$ 。则提出以下标准非线性规划问题：

$\begin{array}{l}\min f\left(x,y\right)\\ \text{s}\text{.t}\text{.}Ax\le b,\\ w=Nx+My+q,\\ \Phi \left(y,w,\mu \right)=0,\end{array}$ (2.9)

当 $\mu \ge 0$ 且充分逼近于0时，问题(2.9)是原问题(1.1)的一个近似逼近；当 $\mu =0$ 时，问题(2.9)是原问题(1.1)的等价；当 $\mu >0$ 时，问题(2.9)是一个一般的可微非线性规划问题。

对于 $z=\left(x,y,w\right)\in R^{n+2m}$ ，可以定义如下罚函数 $\phi :R^{n+2m}\times \left(0,+\infty \right)\times \left(0,+\infty \right)\to R$ ：

$\phi \left(z,\epsilon ,\mu \right)=f\left(x,y\right)+\epsilon {\displaystyle \underset{i=1}{\overset{m}{\sum }}\left|\varphi \left(y_i,w_i,\mu \right)\right|},$ (2.10)

其中 $\epsilon >0$ 是罚参数.

罚函数 $\phi$ 在z处沿方向 $\text{d}z$ 的方向导数为 ${\phi }^{\prime }\left(z,\epsilon ,\mu ,\text{d}z\right)$ ：

$\begin{array}{l}{\phi }^{\prime }\left(z,\epsilon ,\mu ,\text{d}z\right)\\ =\nabla f{\left(s\right)}^{\text{T}}\text{d}s+\epsilon {\displaystyle \underset{i\in I_{+}}{\sum }\nabla \varphi {\left(t_i,\mu \right)}^{\text{T}}\text{d}{t}_i}+\epsilon {\displaystyle \underset{i\in I_0}{\sum }\left|\nabla \varphi {\left(t_i,\mu \right)}^{\text{T}}\text{d}{t}_i\right|}-\epsilon {\displaystyle \underset{i\in I_{-}}{\sum }\left|\nabla \varphi {\left(t_i,\mu \right)}^{\text{T}}\text{d}{t}_i\right|},\end{array}$ (2.11)

其中下指标集定义如下：

$\begin{array}{l}{I}_{+}\equiv I_{+}\left(y,w,\mu \right)=\left\{i|\varphi \left(y_i,w_i,\mu \right)>0\right\},\\ I_0\equiv I_0\left(y,w,\mu \right)=\left\{i|\varphi \left(y_i,w_i,\mu \right)=0\right\},\\ I_{-}\equiv I_{-}\left(y,w,\mu \right)=\left\{i|\varphi \left(y_i,w_i,\mu \right)<0\right\}.\end{array}$

3. 预备知识与算法

引理3.1 对于任意的 $\left(a,b,\mu \right)\in R^2\times \left(0,+\infty \right)$ ，且 $\left(a,b,\mu \right)\ne \left(0,0,0\right)$ ，有

$0<{\partial }_a\varphi \left(a,b,{\mu }_k\right)<1,0<{\partial }_b\varphi \left(a,b,{\mu }_k\right)<1,i=1,\cdots ,m,$ (3.1)

$0<{\left(\frac{\partial \varphi \left(a,b,{\mu }_k\right)}{\partial a}\right)}^2+{\left(\frac{\partial \varphi \left(a,b,{\mu }_k\right)}{\partial b}\right)}^2\le 2$ (3.2)

成立。

证明：由(2.5)，有(3.1)成立。

接下来证明(3.2)式如下：

$\begin{array}{l}{\left(\frac{\partial \varphi \left(y,w,{\mu }_k\right)}{\partial y}\right)}^2+{\left(\frac{\partial \varphi \left(y,w,{\mu }_k\right)}{\partial w}\right)}^2\\ ={\left(\frac{{\text{e}}^{-\frac{y}{\mu }}}{1+{\text{e}}^{-\frac{y}{\mu }}}\right)}^2+{\left(\frac{{\text{e}}^{-\frac{w}{\mu }}}{1+{\text{e}}^{-\frac{w}{\mu }}}\right)}^2\\ \le {\left(\frac{{\text{e}}^{-\frac{y}{\mu }}}{{\text{e}}^{-\frac{y}{\mu }}}\right)}^2+{\left(\frac{{\text{e}}^{-\frac{w}{\mu }}}{{\text{e}}^{-\frac{w}{\mu }}}\right)}^2\\ =2.\end{array}$

设 $z^k=\left(x^k,y^k,w^k\right)\in X_0$ ，是算法第k步的迭代点， $\mu >0$ ， $B_k$ 是一个 $n+2m$ 阶的对称正定矩阵，在算法的第k步迭代中，求解以下二次序列子问题

$\begin{array}{l}\min \nabla f{\left(x^k{y}^k\right)}^{\text{T}}\left(\begin{array}{c}\text{d}x\\ \text{d}y\end{array}\right)+\frac{1}{2}\left(\text{d}{x}^{\text{T}}\text{d}{y}^{\text{T}}\text{d}{w}^{\text{T}}\right)B_k\left(\begin{array}{c}\text{d}x\\ \text{d}y\\ \text{d}w\end{array}\right),\\ \text{s}\text{.t}\text{.}A{x}^k+A\text{d}x-b\le 0,\\ \text{d}w-N\text{d}x-M\text{d}y=0,\\ \Phi \left(t^k,{\mu }^k\right)+{\nabla }_t\Phi {\left(t^k,\mu \right)}^{\text{T}}\text{d}t=0.\end{array}$ (3.3)

求解(3.3)式，得解 $\text{d}{z}^k=\left(\text{d}{x}^k\text{d}{y}^k\text{d}{w}^k\right)$ ，及相应的KKT乘子 $\left({\lambda }^k,l^k,m^k\right)$ ，定义 $\nabla \Phi \left(y^k,w^k,{\mu }^k\right)\in R^{2m}\times R^m\times R\to R^m$ 。

$\nabla \Phi \left(y^k,w^k,{\mu }^k\right)=\left(\begin{array}{c}{D}_y^k\\ D_w^k\end{array}\right),$ (3.4)

其中

$\begin{array}{l}{D}_y^k=diag\left({\partial }_{y_i}\varphi \left(y_i^k,w_i^k,{\mu }_k\right)\right),i=1,\cdots ,m\\ D_w^k=diag\left({\partial }_{w_i}\varphi \left(y_i^k,w_i^k,{\mu }_k\right)\right),i=1,\cdots ,m\end{array}$

从而QP子问题(3.3)也可以写成如下形式：

$\begin{array}{l}\min \nabla f{\left(x^k{y}^k\right)}^{\text{T}}\left(\begin{array}{c}\text{d}x\\ \text{d}y\end{array}\right)+\frac{1}{2}\left(\text{d}{x}^{\text{T}}\text{d}{y}^{\text{T}}\text{d}{w}^{\text{T}}\right)B_k\left(\begin{array}{c}\text{d}x\\ \text{d}y\\ \text{d}w\end{array}\right),\\ \text{s}\text{.t}\text{.}A{x}^k+A\text{d}x-b\le 0,\\ \left(\begin{array}{lll}N\hfill & M\hfill & -I\hfill \\ 0\hfill & D_y^k\hfill & D_w^k\hfill \end{array}\right)\left(\begin{array}{c}\text{d}x\\ \text{d}y\\ \text{d}w\end{array}\right)+\left(\begin{array}{c}0\\ \Phi \left(y_i,w_i,\mu \right)\end{array}\right)=0.\end{array}$ (3.5)

QP子问题(3.3)是一个带线性约束的凸规划问题，所以问题(3.5)的最优解与KKT点等价，进而其KKT系统可以写成：

$\left(\begin{array}{c}\nabla f\left(x^k\right)\\ \nabla f\left(y^k\right)\\ 0\end{array}\right)+B_k\left(\begin{array}{c}\text{d}x\\ \text{d}y\\ \text{d}w\end{array}\right)+\left(\begin{array}{cc}{N}^{\text{T}}& 0\\ M^{\text{T}}& D_y^k\end{array}\right)\left(\begin{array}{c}l\\ m\end{array}\right)+\left(\begin{array}{c}{A}^{\text{T}}\\ 0\\ 0\end{array}\right)\lambda =0,$ (3.6)

$0\le \left(b-A{x}^k-A\text{d}x\right)\perp \lambda \ge 0,$ (3.7)

$\left(\begin{array}{lll}N\hfill & M\hfill & -I\hfill \\ 0\hfill & D_y^k\hfill & D_w^k\hfill \end{array}\right)\text{d}z+\left(\begin{array}{c}0\\ \Phi \left(y_i,w_i,\mu \right)\end{array}\right)=\text{0}，$ (3.8)

其中 $\left(\lambda ,l,m\right)\in R^p\times R^m\times R^m$ 是KKT乘子。

命题3.1 [6] 设假设H 2.1成立， $z^k=\left(x^k,y^k,w^k\right)\in X$ ，罚参数 ${\mu }_k$ 是个给定的正数， $B_k$ 是一个对称正定矩阵，那么QP子问题(3.3)有唯一最优解。

现假设问题(3.3)的唯一最优解为 $\text{d}{z}^k=\left(\text{d}{x}^k,\text{d}{y}^k,\text{d}{w}^k\right)$ ，通过对罚函数 $\phi$ 作线搜索，得到下一个迭代点 $z^{k+1}=\left(x^{k+1},y^{k+1},w^{k+1}\right)$ ，也就是说对于给定的参数 $\rho ,\sigma$ ，满足 $0<\rho <1,0<\sigma <1$ ，第 $k+1$ 次迭代点 $z^{k+1}=z^k+{\alpha }_k\text{d}{z}^k$ ，其中 ${\alpha }_k$ 是步长因子，令 ${\alpha }_k={\delta }^{m_k}$ ， $m_k$ 是使得下列不等式成立的最小的非负整数：

$\phi \left(z^k+{\delta }^{m_k}\text{d}{z}^k,\epsilon ,{\mu }_k\right)\le \phi \left(z^k,\epsilon ,{\mu }_k\right)+\sigma {\delta }^{m_k}{\phi }^{\prime }\left(z^k,\epsilon ,{\mu }_k;\text{d}{z}^k\right).$ (3.9)

下面给出求解问题(1.1)的SQP算法。

算法3.1

步骤0. 选取初始点 $z^0=\left(x^0,y^0,z^0\right)\in X_0$ ，参数 $\xi ,{\mu }_0,{\epsilon }_{-1}\in \left(0,+\infty \right)$ ， $B_0\in R^{n+2m}\times R^{n+2m}$ 。

选取序列 ${\left\{{\mu }_k\right\}}_{k=0}^{\infty }$ 满足 ${\mu }_k>0,{\mu }_{k+1}<{\mu }_k,\underset{k\to \infty }{\lim }\frac{{\mu }_{k+1}}{{\mu }_k^{\zeta }}=\eta$ ，其中 $0<\eta <1,1<\zeta <2$ ，令 $k=0$ 。

步骤1. 求解QP子问题(3.3)得唯一最优解 $\text{d}{z}^k=\left(\text{d}{x}^k,\text{d}{y}^k,\text{d}{w}^k\right)$ 以及相应的乘子 $\left({\lambda }^k,l^k,m^k\right)$ ，如果 $\text{d}{z}^k=0$ ，则令 $z^{k+1}=z^k$ ， ${\epsilon }_k={\epsilon }_{k-1}$ ，转到步骤4。

步骤2. 调整罚参数，计算

${\epsilon }_k=\{\begin{array}{l}{\epsilon }_{k-1},{\epsilon }_{k-1}\ge {\Vert c^k\Vert }_{\infty }+\xi ,\\ \max \left\{{\Vert c^k\Vert }_{\infty }+\xi ,{\epsilon }_{k-1}+2\xi \right\},{\epsilon }_{k-1}<{\Vert c^k\Vert }_{\infty }+\xi ,\end{array}$ (3.10)

其中 ${\Vert c^k\Vert }_{\infty }=\underset{1\le i\le m}{\max }\left|c_i^k\right|$ 。

步骤3. 对函数 $\phi \left(\cdot ,{\epsilon }_k,{\mu }_k\right)$ 按(3.10)做Armijo先搜索求得步长 ${\alpha }_k$ ,令 $z^{k+1}=z^k+{\alpha }_k\text{d}{z}^k$ 。

步骤4. 如果 $\text{d}{z}^k=0$ ，停止；否则，令 ${\mu }_{k+1}=\beta {\mu }_k$ ，更新修正 $B_k$ ，得新的对称正定阵， $B_{k+1}\in R^{\left(n+2m\right)\times \left(n+2m\right)}$ ，令 $k=k+1$ ，返回步骤1。

引理3.2 若假设H2.1成立， $\left(x^k,y^k,w^k\right)\in X_0$ ， ${\mu }_k>0$ ， $B_k$ 是对称正定阵，假设 $\text{d}{z}^k=\left(\text{d}{x}^k,\text{d}{y}^k,\text{d}{w}^k\right)$ 是子问题(3.3)的唯一最优解，相应的KKT乘子为 $\left({\lambda }^k,l^k,m^k\right)$ ，若 ${\epsilon }_k\ge {\Vert c^k\Vert }_{\infty }$ ，有

${\phi }^{\prime }\left(z^k,{\epsilon }^k,{\mu }_k;\text{d}{z}^k\right)\le {\left(-\text{d}{z}^k\right)}^{\text{T}}{B}_k\text{d}{z}^k<0$ (3.11)

成立，且对任意的 $\rho \in \left(0,1\right)$ ，存在 $\bar{\rho }>0$ ，使得对于所有的 $\rho \in \left(0,\bar{\rho }\right)$ ，有

$\phi \left(z^k+\rho \text{d}z,\epsilon ,{\mu }_k\right)\le \phi \left(z^k,\epsilon ,{\mu }_k\right)+\sigma \rho {\phi }^{\prime }\left(z^k,\epsilon ,{\mu }_k;\text{d}{z}^k\right).$ (3.12)

证明：由(3.5)有

$D_y^k\text{d}{y}^k+D_w^k\text{d}{w}^k+\Phi \left(y^k,w^k,{\mu }_k\right)=0,$

${\left(\begin{array}{c}{\nabla }_y\varphi \left(y_i,w_i,\mu \right)\\ {\nabla }_w\varphi \left(y_i,w_i,\mu \right)\end{array}\right)}^{\text{T}}\left(\begin{array}{c}\text{d}{y}_i^k\\ \text{d}{w}_i^k\end{array}\right)+\Phi \left(y_i^k,w_i^k,\mu \right)=0.$ (3.13)

由(2.10)和(2.11)有

${\phi }^{\prime }\left(z^k,{\epsilon }_k,{\mu }_k;\text{d}{z}^k\right)=\nabla f{\left(x^k{y}^k\right)}^{\text{T}}\left(\begin{array}{c}\text{d}{x}^k\\ \text{d}{y}^k\end{array}\right)-{\epsilon }_k{\displaystyle \sum \left|\varphi \left(y^k,w^k,{\mu }_k\right)\right|}.$ (3.14)

由(3.6)得

$\begin{array}{l}\nabla f{\left(x^k{y}^k\right)}^{\text{T}}\left(\begin{array}{c}\text{d}{x}^k\\ \text{d}{y}^k\end{array}\right)+{\left({\lambda }^k\right)}^{\text{T}}A\text{d}{x}^k+{\left(\text{d}{z}^k\right)}^{\text{T}}{B}_k\text{d}{z}^k\\ +{\left(l^k\right)}^{\text{T}}\left(N\text{d}x+M\text{d}y-\text{d}{w}^k\right)+{\left(m^k\right)}^{\text{T}}\left(D_y^k\text{d}{y}^k+D_w^k\text{d}{w}^k\right)=0,\end{array}$

故

$\begin{array}{l}\nabla f{\left(x^k{y}^k\right)}^{\text{T}}\left(\begin{array}{c}\text{d}{x}^k\\ \text{d}{y}^k\end{array}\right)+{\left(\text{d}{z}^k\right)}^{\text{T}}{B}_k\text{d}{z}^k-{\left(m^k\right)}^{\text{T}}\Phi \left(y^k,w^k,{\mu }_k\right)\\ ={\left({\lambda }^k\right)}^{\text{T}}\left(Ax-b\right)\le 0.\end{array}$ (3.15)

由(3.14)，(3.15)有

$\begin{array}{l}{\phi }^{\prime }\left(z^k,{\epsilon }_k,{\mu }_k;\text{d}{z}^k\right)\\ \le -{\left(\text{d}{z}^k\right)}^{\text{T}}{B}_k\text{d}{z}^k+{\displaystyle \underset{i\in I_{+}^k}{\sum }\left(c_i^k-{\epsilon }_k\right)\varphi \left(y_i^k,w_i^k,{\mu }_k\right)}+{\displaystyle \underset{i\in I_{-}^k}{\sum }\left(c_i^k+{\epsilon }_k\right)\varphi \left(y_i^k,w_i^k,{\mu }_k\right)},\end{array}$

其中 $I_{+}^k,I_{-}^k$ 分别定义如下：

$I_{+}^k=\left\{i:\varphi \left(y_i^k,w_i^k,{\mu }_k\right)>0\right\},$

$I_{-}^k=\left\{i:\varphi \left(y_i^k,w_i^k,{\mu }_k\right)<0\right\}.$

又由(3.10)及 ${\epsilon }_k\ge {\Vert c^k\Vert }_{\infty }$ ，有

${\phi }^{\prime }\left(z^k,{\epsilon }^k,{\mu }_k;\text{d}{z}^k\right)\le {\left(-\text{d}{z}^k\right)}^{\text{T}}{B}_k\text{d}{z}^k<0.$

由(2.10)，(3.10) (3.11)和(3.14)易得(3.13)成立.

4. 收敛性分析

本节分析算法的收敛性，为此，我们对对称正定矩阵 $B_k$ 以及算法产生的无穷点列 $\left\{z_k\right\}$ 作如下基本假设：

H4.1 无穷点列 $\left\{z_k\right\}$ 有界。

H4.2 若H4.1成立，序列 $\left\{z_k\right\}$ 的每个极限点 $z^{\ast }=\left(x^{\ast },y^{\ast },w^{\ast }\right)$ ，满足下层非退化条件： $\left(y^{\ast },w^{\ast }\right)\ne \left(0,0\right),i=1,\cdots ,m.$

H4.3 存在常数 $a,b>0$ ，有

$a{\Vert z\Vert }^2\le z^{\text{T}}{B}_{k}z\le b\Vert z\Vert {|}^2,\forall z\in R^{n+2m},k=1,2,\cdots .$

根据上面的几个假设，有以下引理成立.

引理4.1 [5] 如果假设H4.1~H4.3成立，则

1) 存在一个常数 $C>0$ ，有

$\Vert {\left(\begin{array}{cc}M& -I\\ D_y^k& D_W^k\end{array}\right)}^{-1}\Vert \le C,\forall k\in K.$

2) 序列 $\left\{\left({\lambda }^k,l^k,m^k\right)\right\}$ ， $\left\{\text{d}{z}^k\right\}$ 有界。

3) 存在一个正整数 $k_0$ ，使得 ${\epsilon }_k={\epsilon }_{k_0}=\epsilon ,\forall k\ge k_0$ 。

命题4.2 [6] 若数列 $\left\{v_k\right\},\left\{{\gamma }_k\right\}$ 满足

$v_k\ge 0,{\displaystyle \underset{k=1}{\overset{\infty }{\sum }}{v}_k<\infty },{\gamma }_{k+1}\le {\gamma }_k+v_k,k=1,2,\cdots .$ (4.1)

则

1) 数列 $\left\{{\gamma }_k\right\}$ 有上界。即

$\underset{k\to \infty }{\overset{\_\_\_\_}{\lim }}{\gamma }_k<+\infty ;$

2) 数列 $\left\{{\gamma }_k\right\}$ 整列收敛。

命题4.3 若H4.1成立，则对任意的 ${\delta }_1>{\delta }_2>0$ ，点列 $\left(y^k,w^k\right)\in R^2$ ，有不等式

$\left|\varphi \left(y^k,w^k,{\delta }_2\right)\right|\le \left|\varphi \left(y^k,w^k,{\delta }_1\right)\right|+2\left(\ln 2+\frac{\tilde{M}}{{\delta }_2}\right){\delta }_1$

成立。

证明：存在一个 $\bar{\delta }>0$ ，且 ${\delta }_1<\bar{\delta }<{\delta }_2$ ，由中值定理，有

$\begin{array}{l}\left|\varphi \left(y^k,w^k,{\delta }_2\right)\right|\\ \le \left|\varphi \left(y^k,w^k,{\delta }_1\right)+{{\varphi }^{\prime }}_{\mu }\left(y^k,w^k,\bar{\delta }\right)\left({\delta }_2-{\delta }_1\right)\right|\\ \le \left|\varphi \left(y,w,{\delta }_1\right)\right|+\left[\ln \left(1+{\text{e}}^{-\frac{y}{\mu }}\right)+\frac{\frac{y}{\mu }\cdot {\text{e}}^{-\frac{y}{\mu }}}{1+{\text{e}}^{-\frac{y}{\mu }}}+\ln \left(1+{\text{e}}^{-\frac{w}{\mu }}\right)+\frac{\frac{w}{\mu }\cdot {\text{e}}^{-\frac{w}{\mu }}}{1+{\text{e}}^{-\frac{w}{\mu }}}\right]\left({\delta }_2-{\delta }_1\right)\\ \le \left|\varphi \left(y,w,{\delta }_1\right)\right|+\left[2\ln 2+\frac{1}{\mu }\left(\frac{y}{\mu }+\frac{w}{\mu }\right)\right]\left({\delta }_2-{\delta }_1\right)\end{array}$

$\begin{array}{l}\le \left|\varphi \left(y,w,{\delta }_1\right)\right|+\left[2\ln 2+\frac{2}{{\mu }^2}\max \left\{y,w\right\}\right]\left({\delta }_2-{\delta }_1\right)\\ \le \left|\varphi \left(y,w,{\delta }_1\right)\right|+2\left[\ln 2+\frac{1}{\delta }\max \left\{y,w\right\}\right]\left({\delta }_2-{\delta }_1\right)\\ \le \left|\varphi \left(y^k,w^k,{\delta }_1\right)\right|+2\left(\ln 2+\frac{\tilde{M}}{{\delta }_2}\right){\delta }_1.\end{array}$

由引理3.2，命题4.2及命题4.3易知有以下结论成立.

引理4.2 若H4.1成立，则

$\underset{k\to \infty }{\lim }\phi \left(z^k,{\epsilon }^k,{\mu }_k\right)=\underset{k\to \infty }{\lim }\phi \left(z^{k+1},{\epsilon }^k,{\mu }_k\right)=\phi \left(z^{\ast },\epsilon ,0\right).$

由引理4.1(2)，不妨做如下假设：

$\text{d}{z}^k=\text{d}{z}^{\ast }=\left(\text{d}{x}^{\ast },\text{d}{y}^{\ast },\text{d}{w}^{\ast }\right),x^k\to x^{\ast },B_k\to B_{\ast },l^k\to l^{\ast },m^k\to m^{\ast },k\in K.$ (4.2)

利用上述引理及假设,来证明本文的算法具有全局收敛性。

引理 4.3 若假设H 4.1-H4.3成立，则序列 $\left\{\text{d}{z}^k\right\}$ 收敛到0，即

$\underset{k\in K}{\lim }\text{d}{z}^k=\text{d}{z}^{\ast }=0.$ (4.3)

证明：我们用反证法来证明，假设 $\text{d}{z}^{\ast }\ne 0$ 。

由(3.12)有

${\phi }^{\prime }\left(z^{\ast },\epsilon ,0;\text{d}{z}^k\right)=-{\left(\text{d}{z}^{\ast }\right)}^{\text{T}}{B}_{\ast }\text{d}{z}^{\ast }<0.$

当 $k\to \infty$ ，序列 $\left\{z^k\right\}$ 收敛于 $z^{\ast }$ ，记

$\underset{k\in K}{\lim }{\phi }^{\prime }\left(z^k,\epsilon ,{\mu }_k;\text{d}{z}^k\right)={\phi }^{\prime }\left(z^{\ast },\epsilon ,0;\text{d}{z}^{\ast }\right).$

考虑算法步骤3的线搜索，步长 $\alpha$ 在K上有界，且远大于0，设

${\alpha }_k\ge {\alpha }_{\ast }=\inf \left({\alpha }_k,k\in K\right)>0,k\in K.$

又由序列 $\left\{\phi \left(z^k,\epsilon ,{\mu }_k\right)\right\}$ 为单调递减序列，当 ${\mu }_k$ 固定时，由序列 $\left\{z^k\right\}$ 是收敛到 $z^{\ast }$ 的,及函数 $\phi$ 的连续性，当 $k\to \infty$ 时，有 $\phi \left(z^k,\epsilon ,{\mu }_k\right)\to \phi \left(z^{\ast },\epsilon ,0\right)$ 。

由(3.9)及引理2.4.1(3)，有

$0=\underset{k\in K}{\lim }\left(\phi \left(z^{k+1},\epsilon ,{\mu }_k\right)-\phi \left(z^k,\epsilon ,{\mu }_k\right)\right)\le \underset{k\in K}{\lim }\delta \alpha {\phi }^{\prime }\left(z^k,\epsilon ,{\mu }_k;\text{d}{z}^k\right)\le \frac{1}{2}\delta {\alpha }_{\ast }{\phi }^{\prime }\left(z^k,\epsilon ,0;\text{d}{z}^k\right),$

当 $k\in K,k\to \infty$ 时。这与(3.12)矛盾，从而假设 $\text{d}{z}^{\ast }\ne 0$ 不成立，即(4.3)成立.

定理4.1 H4.1~H4.3成立，则算法3.1产生的点列 $\left\{z^k\right\}$ ，它的每一个聚点 $z^{\ast }=\left(x^{\ast },y^{\ast },w^{\ast }\right)$ 都是问题(1.1)

的稳定点，即本文的算法是全局收敛的。

证明：由已知，序列 $\left\{z^k\right\}$ 是个无穷序列， $z^k\to z^{\ast }$ ，有 $z^{\ast }\in X_0$ ,即

$A{x}^{\ast }\le b,w^{\ast }=N{x}^{\ast }+M{y}^{\ast }+q,$

由(2.7)和(3.3)有

$\begin{array}{l}\nabla f\left(x^{\ast },y^{\ast }\right)+A^{\text{T}}{\lambda }^{\ast }+N^{\text{T}}{l}^{\ast }=0,\\ M^{\text{T}}{l}^{\ast }+D_y^{\ast }{m}^{\ast }=0,\\ l^{\ast }=D_w^{\ast }{m}^{\ast },\\ 0\le {\lambda }^{\ast }\perp \left(b-A{x}^{\ast }\right)\ge 0,\\ \underset{k\to \infty }{\lim }\Phi \left(a,b,{\mu }_k\right)=\min \left(a,b\right)=0,\end{array}$ (4.4)

其中 $D_y^{\ast }$ ， $D_w^{\ast }$ 定义为

$D_y^{\ast }=diag\left({\nabla }_{y_i}\varphi \left(y_i^{\ast },w_i^{\ast },0\right),i=1,\cdots ,m\right),$

$D_w^{\ast }=diag\left({\nabla }_{w_i}\varphi \left(y_i^{\ast },w_i^{\ast },0\right),i=1,\cdots ,m\right).$

由(2.2)和(3.11)有

$0\le w^{\ast }\perp y^{\ast }\ge 0,$

又 $z^{\ast }\in X_0$ ，从而有 $z^{\ast }\in X$ 。

综上所述， $z^{\ast }$ 是问题(1.1)的稳定点。

5. 数值实验

本文提出了求解均衡约束优化问题的算法，在这一节我们将对这个算法进行数值试验，利用MATLAB软件算法进行了数值实验，通过实验说明这个算法具有有效性和可行性.

算法中的矩阵 $B_k$ 是Lagrange函数 $L\left(z,\mu \right)=f\left(z\right)+{\displaystyle \underset{j=1}{\overset{m}{\sum }}{\mu }_j{g}_j\left(z,\mu \right)}$ 的二阶Hessian阵的近似估计。矩阵

$B_k$ 采用BFGS公式修正：

$B_{k+1}=B_k-\frac{B_k{\hat{s}}_k{\left(B_k{\hat{s}}_k\right)}^{\text{T}}}{{\hat{s}}_k^{\text{T}}{B}_k{\hat{s}}_k}+\frac{{\hat{\eta }}_k{\hat{\eta }}_k^{\text{T}}}{{\hat{s}}_k^{\text{T}}{\hat{\eta }}_k},$

其中 ${\hat{s}}_k=z^{k+1}-z^k,{\hat{\eta }}_k={\theta }_k{\hat{\gamma }}_k+\left(1-\theta \right)B_k{\hat{s}}_k$ ，这里

${\hat{\gamma }}_k=\nabla L_z\left(z^k+{\hat{s}}_k,{\mu }_k\right)-\nabla L_z\left(z^k,{\mu }_k\right),$

$\theta =\{\begin{array}{l}1,{\hat{s}}_k{\hat{\gamma }}_k\ge 0.2{\hat{s}}_k^{\text{T}}{B}_k{\hat{s}}_k\\ \frac{0.8{\hat{s}}_k^{\text{T}}{B}_k{\hat{s}}_k}{{\hat{s}}_k^{\text{T}}{B}_k{\hat{s}}_k-{\hat{s}}_k{\hat{\gamma }}_k},否则.\end{array}$

分别选取参数 $\xi =10$ ， ${\epsilon }_0=10$ ， $B_0$ 为 $\left(n+2m\right)\times \left(n+2m\right)$ 的单位阵.

问题1

$\begin{array}{l}\min f\left(x\right)=\frac{1}{2}{x}^2+\frac{1}{2}xy-95x\\ \text{s}\text{.t}.0\le x\le 200,\\ 0\le w\perp y\ge 0.\end{array}$

问题2

$\begin{array}{l}\min f\left(x\right)=\frac{1}{2}\left[{\left(x_1+x_2+y_1-15\right)}^2+{\left(x_1+x_2+y_2-15\right)}^2\right]\\ \text{s}\text{.t}\text{.}0\le x\le 10,\\ w=Nx+My+q,\\ 0\le w\perp y\ge 0,\\ \left(x,y,w\right)\in R^2\times R^2\times R^2.\end{array}$

其中

$q=\left(\begin{array}{c}-36\\ -25\end{array}\right),M=\left(\begin{array}{cc}2& 8/3\\ 5/4& 2\end{array}\right),N=\left(\begin{array}{cc}8/3& 2\\ 2& 5/4\end{array}\right).$

问题3

$\begin{array}{l}\min f\left(x,y\right)=\frac{1}{2}{\left(x_1+3y_1-4y_2-5\right)}^2+{\left(x_2-\frac{15}{8}{y}_1+3y_2-9\right)}^2\\ \text{s}\text{.t}\text{.}0\le x_1\le 10,\\ 0\le x_2\le 10,\\ 0\le y\perp \left[{\left(\begin{array}{cc}0& -1\\ -1& 0\end{array}\right)}^{\text{T}}x+{\left(\begin{array}{cc}3& -\frac{15}{8}\\ -4& 3\end{array}\right)}^{\text{T}}y+\left(\begin{array}{c}10\\ 6\end{array}\right)\right]\ge 0.\end{array}$

问题4

$\begin{array}{l}\min f\left(x\right)=\frac{1}{2}{x}^{\text{T}}x+e^{\text{T}}y\\ \text{s}\text{.t}\text{.}Ax\le b,\\ w=Nx+My+q,\\ 0\le w\perp y\ge 0.\end{array}$

其中， $M\in R^{m\times m},N\in R^{m\times n},A\in R^{p\times n}$ 分别为随机产生的矩阵，矩阵M严格对角占优的，是P矩阵，其非对角线上的非负向量 $q\in R^m,b\in R^p$ 是随机产生的，对于单位向量 $e\in R^n,A,b$ 满足 $Ae<b$ 。

表1和表2是算法3.1的数值实验结果。

Table 1. The numerical experiment results of problem 1 - problem 3

表1. 问题1~问题3的数值实验结果

Table 2. The numerical experiment results of problem 4

表2. 问题4的数值实验结果

在表1和表2中，NO.表示问题的编号， $\left(p,m,n\right)$ 分别表示x的维数，m表示y的维数，n表示约束的个数。 $\Vert \text{d}{z}^k\Vert$ 表示终止时， $\Vert \text{d}{z}^k\Vert$ 的对应值。数值实验表明本文的算法3.1是有效的。

基金项目

河池学院2019年高层次人才科研启动费项目“比例时滞神经网络的周期解与稳定状态研究”(2019GCC005)。