一类变分不等式问题的二阶微分方程方法

The Second Order Differential Equation Method for a Class of Variational Inequality Problem

摘要: 本文运用光滑化的自然残差函数建立了具有不等式约束条件的变分不等式问题的光滑化KKT方程组，并建立了与其等价的无约束优化问题。建立了具有阻尼惯性参数和时间尺度参数的二阶微分方程系统来求解该无约束优化问题，并证明了该二阶微分方程系统的稳定性，从而得到了具有不等式约束的变分不等式问题的KKT点的收敛性。并将二阶微分方程方法与已有的一阶微分方程方法进行了理论条件和数值结果的对比。在理论条件的要求上，二阶微分方程方法的条件要更容易实现，而在数值结果上，一阶微分方程方法的收敛速度要快，但是两种方法的差距可以忽略不计。

Abstract: In this paper, a smoothing KKT equation system with inequality constraints is established by using the smoothing natural residual function, and the unconstrained optimization problem is established. A second order differential equation system with damping coefficient and time scale coefficient is established to solve the unconstrained optimization problem, and the stability of the second-order differential equation system is proved, then the convergence of the KKT points of the variational inequality problem with inequality constraints is obtained. The theoretical conditions and numerical results of the second-order differential equation method are compared with the existing first order differential equation method. In terms of theoretical conditions, the conditions of the second-order differential equation method are easier to implement, while in the numerical results, the convergence speed of the first order differential equation method is faster, but the difference between the two methods is negligible.

1. 引言

变分不等式(Variational Inequality，简称VI)是构成非线性分析的重要组成部分，它广泛地应用在微分方程、力学、数理经济、控制论、优化理论、对策理论及非线性规划等理论和应用学科。随着经典变分不等式的发展，产生了各类变分不等式问题，带约束的变分不等式问题是将经典变分问题的约束性条件放宽到某一些单边约束的一种变分方法，多被用来研究最优控制、偏微分方程等，是变分学研究中的一个重大发现。本文将要研究的变分不等式如下：

$\langle F\left(x^{*}\right),y-x^{*}\rangle \ge 0,\forall y\in K,$ (1.1)

其中，约束集合K如下表示

$K=\left\{x\in R^n|g\left(x\right)\le 0\right\},$

$\langle \cdot ,\cdot \rangle$ 表示欧式空间的内积，$F:R^n\to R^n,$ $g:R^n\to R^m$ 是映射，且都是连续可微的。$R^n$ 表示为n维实向量空间。

变分不等式问题的研究可以分为理论研究与算法设计两大类。本文是对具有不等式约束条件的变分不等式问题(1.1)的微分方程方法的研究。Facchinei和Pang专著[1]收集了n维欧氏空间中的互补问题和变分不等式问题的几乎所有重要成果，却没有记录微分方程方法求解变分不等式问题。用微分方程方法求解变分不等式问题的工作还值得继续深入开展。在投影算子的基础上，Gao、Liao和Qi [2]提出了求解具有线性和非线性约束条件的变分不等式问题的微分方程系统，给出了三个充分条件以确保所提出的具有非对称映射的微分方程系统是Lyapunov意义下稳定的，且收敛于原始问题的一个确定的解。其余结果还可以参见文献He和Yang [3]，Zhou、Zhang和He [4]，Liao、Qi和Qi [5]。不必构造Lyapunov函数，Antipin [6]提出了具有双约束条件的变分不等式问题，运用投影算子，建立了微分方程系统。同时，证明了该微分方程系统的轨迹收敛于具有双约束条件的变分不等式问题的解。王莉等学者在Antipin学者工作的基础之上，继续研究了不同种类的变分不等式问题的微分方程方法，详见参考文献[7]-[9]。近几年，Attouch [10]-[13]运用二阶微分方程系统求解了不同类型的优化问题。运用这一思想，在此基础上，将构造具有阻尼惯性参数和时间尺度参数的二阶微分方程系统来求解具有不等式约束条件的变分不等式的问题(1.1)，并证明它们的收敛性定理。这与传统的微分方程方法相比，所需条件比较简洁。又由于微分方程系统的求解有很成熟的数值解法，比如Matlab中有Runge-Kutta龙格－库塔算法等。将这些数值软件用于求解该二阶微分方程系统会有良好的效果。

2. 预备知识

设集合C是一个凸闭集合，对于任意$x\in R^n$ ，都存在唯一的$\overline{x}\in C$ ，使得

$\Vert x-\overline{x}\Vert =\min \left\{\Vert x-y\Vert ,y\in C\right\}$ ,

则称$\overline{x}$ 为点x在集合C上的投影点。

对任意$x,y\in R^n$ ，有$\Vert {\Pi }_C\left(y\right)-{\Pi }_C\left(x\right)\Vert <\Vert y-x\Vert$ 成立。其中${\Pi }_C:R^n\to C$ 是定义好的映射，称之为投影算子。

引理2.1 [14]：设H为Hilbert空间，集合$C\subset H$ 为一闭凸集，对任意的$z\in H$ ，则存在$u\in C$ 满足

$\langle u-z,v-u\rangle \ge 0,\forall v\in C,$

当且仅当$u-{\Pi }_C\left(z\right)=0$ 。

受到Attouch [11]-[13]的启示，为了建立二阶微分方程系统，需要对原始的变分不等式问题(1.1)进行相应的转化，得到一个无约束优化问题。光滑化互补函数在该转化过程中起着重要的作用，具体介绍如下：

所谓互补函数是指对于$\phi :R^m\times R^m$ ，满足$\phi \left(x,y\right)=0$ 当且仅当$x\perp y$ 成立。

本文将主要运用的互补函数为自然残差函数(NR函数)

${\phi }_{NR}\left(x,y\right)=x-{\Pi }_K\left(x-y\right),\forall x,y\in R^n,$ (2.1)

其中${\Pi }_K\left(x-y\right)$ 表示$x-y$ 到集合K上的投影算子。根据互补函数的定义可知，对任意的$x,y\in R^n$ ，有${\phi }_{NR}\left(x,y\right)=0$ 当且仅当$\langle x,y\rangle =0$ 。光滑化的NR函数表示为

${\phi }_{NR}^{\epsilon }\left(x,y\right)=x-\varphi \left(\epsilon ,x-y\right),$

其中，对于$\left(\epsilon ,x-y\right)\in R\times R^n$ ，关于函数$\varphi :R\times R^n\to R^n$ 有

$\varphi \left(\epsilon ,x-y\right)=\frac{1}{2}\left(x-y+\sqrt{{\epsilon }^{2}e+{\left(x-y\right)}^2}\right),$

其中e表示$R^n$ 的单位向量。

具有不等式约束条件的变分不等式问题(1.1)的KKT条件

$\{\begin{array}{l}L\left(x,\lambda \right)=0,\\ \lambda \perp g\left(x\right),\\ -g\left(x\right)\ge 0,\lambda \ge 0,\end{array}$ (2.2)

其中

$L\left(x,\lambda \right)=F\left(x\right)+Jg{\left(x\right)}^{\text{T}}\lambda ,$

$L\left(x,\lambda \right)$ 为变分不等式的拉格朗日函数。对于$g\left(x\right)\perp \lambda$ , 其等价于

$-g_i\left(x\right)\perp {\lambda }_i,-g_i\left(x\right)\ge ,{\lambda }_i\ge ,i=1,\cdots ,p.$

对此，运用光滑化的NR函数(2.7)，将$g\left(x\right)\perp \lambda$ 表示为

${\phi }_{NR}^{\epsilon }\left(-g_i\left(x\right),{\lambda }_i\right)=0.$

对于$z=\left(\epsilon ,x,\lambda \right)\in R^{1+n+m}$ ，定义如下的光滑化函数

$S\left(z\right)=S\left(\epsilon ,x,\lambda \right)=\left\{\begin{array}{c}\epsilon \\ L\left(x,\lambda \right)\\ {\phi }_{NR}^{\epsilon }\left(-g\left(x\right),\lambda \right)\end{array}\right\},$ (2.3)

其中

${\phi }_{NR}^{\epsilon }\left(-g\left(x\right),\lambda \right)=\left\{\begin{array}{c}{\phi }_{NR}^{\epsilon }\left(-g_1\left(x\right),{\lambda }_1\right)\\ {\phi }_{NR}^{\epsilon }\left(-g_2\left(x\right),{\lambda }_2\right)\\ ⋮\\ {\phi }_{NR}^{\epsilon }\left(-g_m\left(x\right),{\lambda }_m\right)\end{array}\right\}.$ (2.4)

则求解变分不等式问题(1.1)的KKT条件的解转化为求解光滑化方程组

$S\left(z\right)=S\left(\epsilon ,x,\lambda \right)=\left\{\begin{array}{c}\epsilon \\ L\left(x,\lambda \right)\\ {\phi }_{NR}^{\epsilon }\left(-g\left(x\right),\lambda \right)\end{array}\right\}=0.$ (2.5)

接下来，引入效益函数$\Phi \left(z\right)=\frac{1}{2}{\Vert S\left(z\right)\Vert }^2$ ，则解决光滑化方程组(2.5)等价于求解下面的无约束优化问题

$\min \Phi \left(z\right)=\min \frac{1}{2}{\Vert S\left(z\right)\Vert }^2.$ (2.6)

而且容易知道，若益函数$z_0=\left({\epsilon }_0,x_0,{\lambda }_0\right)$ 是(2.6)的解，则$\left(x_0,{\lambda }_0\right)$ 就是具有不等式约束条件的变分不等式问题(1.1)的KKT点。

3. 变分不等式的二阶微分方程方法

本章将建立具有阻尼惯性参数和时间尺度参数的二阶微分方程系统来求解无约束凸优化问题(2.6)，从而实现对不等式约束的变分不等式问题(1.1)的KKT点的求解。

3.1. 二阶微分方程系统的建立

为求解无约束优化问题(2.6)，建立二阶微分方程系统如下

$\left(\begin{array}{c}\ddot{\epsilon }\left(t\right)\\ \ddot{x}\left(t\right)\\ \ddot{\lambda }\left(t\right)\end{array}\right)+\gamma \left(t\right)\left(\begin{array}{c}\ddot{\epsilon }\left(t\right)\\ \ddot{x}\left(t\right)\\ \ddot{\lambda }\left(t\right)\end{array}\right)+\beta \left(t\right)\left(\begin{array}{c}{\nabla }_{\epsilon }S{\left(z\right)}^{\text{T}}S\left(z\right)\\ {\nabla }_{x}S{\left(z\right)}^{\text{T}}S\left(z\right)\\ {\nabla }_{\lambda }S{\left(z\right)}^{\text{T}}S\left(z\right)\end{array}\right)=0,$ (3.1)

其中，$\gamma \left(t\right)$ 表示为阻尼惯性参数，$\beta \left(t\right)$ 为时间尺度参数，且$\gamma ,\beta$ 在$\left[t_0,+\infty \right)$ 上是非负连续函数。由Attouch等[7]-[10]的研究结果知道该系统的解是存在的，并且系统对应的平衡点$z^{*}=\left({\epsilon }^{*},x^{*},{\lambda }^{*}\right)$ 是无约束优化问题(2.6)的解。基于对阻尼惯性参数和时间尺度参数的调整，将分析该二阶微分方程系统解轨迹的全局收敛性。

在初始情形$\Phi =\left(\epsilon ,x,\lambda \right)$ 的情况下，对二阶微分方程系统(3.1)直接积分，可以得到

$p\left(t\right)={\text{e}}^{{\displaystyle {\int }_{t_0}^t\gamma \left(u\right)\text{d}u}}.$

假设其满足$H_0$ 条件

${\displaystyle {\int }_{t_0}^{+\infty }\frac{\text{d}u}{p\left(u\right)}}<+\infty .$

在$H_0$ 条件下，定义函数$\Gamma :\left[t_0,+\infty \right)\to R$ 如下

$\Gamma \left(t\right)=p\left(t\right){\displaystyle {\int }_t^{+\infty }\frac{\text{d}u}{p\left(u\right)}}.$ (3.2)

通过对式(3.2)进行求导，可以得到

$\dot{\Gamma }\left(t\right)=\gamma \left(t\right)\Gamma \left(t\right)-1.$ (3.3)

定义全局能量函数

$W\left(t\right)=\frac{1}{2}{\Vert \begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\Vert }^2+\beta \left(t\right)\left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right).$ (3.4)

以及锚函数$h\left(t\right)$

$h\left(t\right)=\frac{1}{2}{\Vert \left(\begin{array}{c}\epsilon \left(t\right)\\ x\left(t\right)\\ \lambda \left(t\right)\end{array}\right)-\left(\begin{array}{c}{\epsilon }^{*}\\ x^{*}\\ {\lambda }^{*}\end{array}\right)\Vert }^2,$ (3.5)

其中，$\left(\begin{array}{c}{\epsilon }^{*}\\ x^{*}\\ {\lambda }^{*}\end{array}\right)\in \arg \min \Phi$ ，解集$\arg \min \Phi$ 非空。

定义能量函数$\xi :\left[t_0,+\infty \right)\to R^{+}$

$\xi \left(t\right)=\Gamma {\left(t\right)}^{2}W\left(t\right)+h\left(t\right)+\Gamma \left(t\right)h\left(t\right).$ (3.6)

该能量函数在二阶微分方程系统的稳定性分析中起着非常重要的作用。结合(3.4)和(3.5)，可以得到

$\xi \left(t\right)=\Gamma {\left(t\right)}^2\beta \left(t\right)\left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right)+\frac{1}{2}{\Vert \left(\begin{array}{c}\epsilon \left(t\right)\\ x\left(t\right)\\ \lambda \left(t\right)\end{array}\right)-\left(\begin{array}{c}{\epsilon }^{*}\\ x^{*}\\ {\lambda }^{*}\end{array}\right)+\Gamma \left(t\right)\left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right)\Vert }^2.$

3.2. 二阶微分方程系统的稳定性分析

下面给出二阶微分方程系统的稳定性定理。

定理3.1：时间尺度参数$\beta \left(t\right)$ 是一个连续非负函数，阻尼惯性参数$\gamma \left(t\right)$ 是一个满足假设条件$H_0$ 的连续函数且二者满足$H_{\gamma ,\beta }$ 增长条件

$\Gamma \left(t\right)\beta \left(t\right)\le \beta \left(t\right)\left(3-2y\left(t\right)\Gamma \left(t\right)\right)$

则当$t\to +\infty$ 时，二阶微分方程系统(3.1)的解的轨迹$z=\left(\epsilon ,x,\lambda \right)\in R\times R^n\times R^m$ 在$\left[t_0,+\infty \right)$ 上对应的函数值的收敛速度满足

$\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi =\text{O}\left(\frac{1}{\beta \left(t\right)\Gamma {\left(t\right)}^2}\right),$

而且解的轨迹$z=\left(\epsilon ,x,\lambda \right)$ 在$\left[t_0,+\infty \right)$ 上是有界的。

证明：为了对能量函数$\xi \left(t\right)$ 求导，需要对能量函数$\xi \left(t\right)$ 中的全局能量函数$W\left(t\right)$ 和锚函数$h\left(t\right)$ 求导。全局能量函数$W\left(t\right)$ 的导数计算如下

$\begin{array}{c}\dot{W}\left(t\right)=\langle \left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right),\left(\begin{array}{c}\ddot{\epsilon }\left(t\right)\\ \ddot{x}\left(t\right)\\ \ddot{\lambda }\left(t\right)\end{array}\right)\rangle +\beta \left(t\right)\langle \left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right),\left(\begin{array}{c}{\nabla }_{\epsilon }S{\left(z\right)}^{\text{T}}S\left(z\right)\\ {\nabla }_{x}S{\left(z\right)}^{\text{T}}S\left(z\right)\\ {\nabla }_{\lambda }S{\left(z\right)}^{\text{T}}S\left(z\right)\end{array}\right)\rangle +\dot{\beta }\left(t\right)\left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right)\\ =\langle \left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right),\left(\begin{array}{c}\ddot{\epsilon }\left(t\right)\\ \ddot{x}\left(t\right)\\ \ddot{\lambda }\left(t\right)\end{array}\right)+\beta \left(t\right)\left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right),\left(\begin{array}{c}{\nabla }_{\epsilon }S{\left(z\right)}^{\text{T}}S\left(z\right)\\ {\nabla }_{x}S{\left(z\right)}^{\text{T}}S\left(z\right)\\ {\nabla }_{\lambda }S{\left(z\right)}^{\text{T}}S\left(z\right)\end{array}\right)\rangle +\dot{\beta }\left(t\right)\left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right)\\ =\langle \left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right),-\gamma \left(t\right)\left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right)\rangle +\dot{\beta }\left(t\right)\left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right),\end{array}$

继续整理为

$\dot{W}\left(t\right)=-\gamma \left(t\right){\Vert \left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right)\Vert }^2+\dot{\beta }\left(t\right)\left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right).$

锚函数求一阶导数可得

$\dot{h}\left(t\right)=\langle \left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right),\left(\begin{array}{c}\epsilon \left(t\right)-{\epsilon }^{*}\\ x\left(t\right)-x^{*}\\ \lambda \left(t\right)-{\lambda }^{*}\end{array}\right)\rangle ,$

继续求二阶导数可得

$\ddot{h}\left(t\right)={\Vert \left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right)\Vert }^2+\langle \left(\begin{array}{c}\ddot{\epsilon }\left(t\right)\\ \ddot{x}\left(t\right)\\ \ddot{\lambda }\left(t\right)\end{array}\right),\left(\begin{array}{c}\epsilon \left(t\right)-{\epsilon }^{*}\\ x\left(t\right)-x^{*}\\ \lambda \left(t\right)-{\lambda }^{*}\end{array}\right)\rangle .$

综合以上结果及(3.6)计算能量函数导数如下

$\begin{array}{c}\ddot{\xi }\left(t\right)=2\Gamma \left(t\right)\dot{\Gamma }\left(t\right)W\left(t\right)+\Gamma {\left(t\right)}^2\dot{W}\left(t\right)+\dot{h}\left(t\right)+\dot{\Gamma }\left(t\right)\dot{h}\left(t\right)+\Gamma \left(t\right)\ddot{h}\left(t\right)\\ =2\Gamma \left(t\right)\dot{\Gamma }\left(t\right)W\left(t\right)+\Gamma {\left(t\right)}^2\dot{W}\left(t\right)+\gamma \left(t\right)\Gamma \left(t\right)\dot{h}\left(t\right)+\Gamma \left(t\right)\ddot{h}\left(t\right)\\ =2\Gamma \left(t\right)\dot{\Gamma }\left(t\right)W\left(t\right)+\Gamma {\left(t\right)}^2\dot{W}\left(t\right)+\Gamma \left(t\right)\left(\gamma \left(t\right)\dot{h}\left(t\right)+\ddot{h}\left(t\right)\right)\\ =-\Gamma \left(t\right){\Vert \dot{z}\left(t\right)\Vert }^2+2\Gamma \left(t\right)\dot{\Gamma }\left(t\right)\beta \left(t\right)\left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right)\\ +\Gamma {\left(t\right)}^2\dot{\beta }\left(t\right)\left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right)+\Gamma \left(t\right)\left(\gamma \left(t\right)\dot{h}\left(t\right)+\ddot{h}\left(t\right)\right),\end{array}$

在上式中${\Vert \dot{z}\left(t\right)\Vert }^2={\Vert \left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right)\Vert }^2$ 。根据二阶微分方程系统(3.1)及$\Phi \left(z\right)$ 的凸性，

继续计算

$\begin{array}{l}\gamma \left(t\right)\dot{h}\left(t\right)+\ddot{h}\left(t\right)\\ =\langle \gamma \left(t\right)\left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right),\left(\begin{array}{c}\epsilon \left(t\right)-{\epsilon }^{*}\\ x\left(t\right)-x^{*}\\ \lambda \left(t\right)-{\lambda }^{*}\end{array}\right)\rangle +\langle \left(\begin{array}{c}\ddot{\epsilon }\left(t\right)\\ \ddot{x}\left(t\right)\\ \ddot{\lambda }\left(t\right)\end{array}\right),\left(\begin{array}{c}\epsilon \left(t\right)-{\epsilon }^{*}\\ x\left(t\right)-x^{*}\\ \lambda \left(t\right)-{\lambda }^{*}\end{array}\right)\rangle +{\Vert \left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right)\Vert }^2\\ =\langle \gamma \left(t\right)\left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right)+\left(\begin{array}{c}\ddot{\epsilon }\left(t\right)\\ \ddot{x}\left(t\right)\\ \ddot{\lambda }\left(t\right)\end{array}\right),\left(\begin{array}{c}\epsilon \left(t\right)-{\epsilon }^{*}\\ x\left(t\right)-x^{*}\\ \lambda \left(t\right)-{\lambda }^{*}\end{array}\right)\rangle +{\Vert \left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right)\Vert }^2\\ =-\beta \left(t\right)\langle \left(\begin{array}{c}{\nabla }_{\epsilon }S{\left(z\right)}^{T}S\left(z\right)\\ {\nabla }_{x}S{\left(z\right)}^{T}S\left(z\right)\\ {\nabla }_{\lambda }S{\left(z\right)}^{T}S\left(z\right)\end{array}\right),\left(\begin{array}{c}\epsilon \left(t\right)-{\epsilon }^{*}\\ x\left(t\right)-x^{*}\\ \lambda \left(t\right)-{\lambda }^{*}\end{array}\right)\rangle +{\Vert \left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right)\Vert }^2\\ \le -\beta \left(t\right)\left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right)+{\Vert \left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right)\Vert }^2\end{array}$

进而可得

$\begin{array}{c}\dot{\xi }\left(t\right)=-\Gamma \left(t\right){\Vert \dot{z}\left(t\right)\Vert }^2+2\Gamma \left(t\right)\dot{\Gamma }\left(t\right)\beta \left(t\right)\left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right)\\ +\Gamma {\left(t\right)}^2\dot{\beta }\left(t\right)\left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right)+\Gamma \left(t\right)\left(\gamma \left(t\right)\dot{h}\left(t\right)+\ddot{h}\left(t\right)\right)\\ \le \left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right)\Gamma \left(t\right)\left(2\dot{\Gamma }\left(t\right)\beta \left(t\right)+\Gamma \left(t\right)\dot{\beta }\left(t\right)-\beta \left(t\right)\right).\end{array}$

由式(3.3)，上式继续化简

$\dot{\epsilon }\left(t\right)\le \left(\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \right)\Gamma \left(t\right)\left(\Gamma \left(t\right)\dot{\beta }\left(t\right)+\beta \left(t\right)\left(2\gamma \left(t\right)\Gamma \left(t\right)-3\right)\right).$

根据条件$H_{\gamma ,\beta }$ 得$\dot{\xi }\left(t\right)\le 0$ ，因此能量函数$\xi \left(t\right)$ 在$\left[t_0,+\infty \right)$ 上是递减的，从而当$t\in \left[t_0,+\infty \right)$ 时有$\xi \left(t\right)\le \xi \left(t_0\right)$ 成立。

根据能量函数的定义(3.6)可以得到，对于任意的$t\ge t_0$ 都有下式成立

$\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi \le \frac{\xi \left(t_0\right)}{\beta \left(t\right)\Gamma {\left(t\right)}^2}$

等价于

$\Phi \left(\epsilon \left(t\right),x\left(t\right),\lambda \left(t\right)\right)-\min \Phi =\text{O}\left(\frac{1}{\beta \left(t\right)\Gamma {\left(t\right)}^2}\right).$

下面讨论二阶微分方程系统(3.1)的解的轨迹的有界性，已知函数$\xi \left(t\right)$ 在$\left[t_0,+\infty \right)$ 上是递减的，对于任意$t\in \left[t_0,+\infty \right)$ 时，有以下不等式成立

${\Vert \left(\begin{array}{c}\epsilon \left(t\right)\\ x\left(t\right)\\ \lambda \left(t\right)\end{array}\right)-\left(\begin{array}{c}{\epsilon }^{*}\\ x^{*}\\ {\lambda }^{*}\end{array}\right)+\Gamma \left(t\right)\left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right)\Vert }^2\le 2\xi \left(t\right)\le 2\xi \left(t_0\right).$

在展开上述不等式后可以得到

${\Vert \left(\begin{array}{c}\epsilon \left(t\right)\\ x\left(t\right)\\ \lambda \left(t\right)\end{array}\right)-\left(\begin{array}{c}{\epsilon }^{*}\\ x^{*}\\ {\lambda }^{*}\end{array}\right)\Vert }^2+2\Gamma \left(t\right)\langle \left(\begin{array}{c}\epsilon \left(t\right)\\ x\left(t\right)\\ \lambda \left(t\right)\end{array}\right)-\left(\begin{array}{c}{\epsilon }^{*}\\ x^{*}\\ {\kappa }^{*}\end{array}\right)，\left(\begin{array}{c}\dot{\epsilon }\left(t\right)\\ \dot{x}\left(t\right)\\ \dot{\lambda }\left(t\right)\end{array}\right)\rangle \le 2\xi \left(t_0\right).$ (3.7)

假设$q\left(t\right)={\displaystyle {\int }_t^{+\infty }\frac{\text{d}s}{p\left(s\right)}}$ ，当满足假设条件$H_0$ 时，则函数$q\left(t\right)$ 在$t\ge t_0$ 时是有界的，而且满足

$\dot{q}\left(t\right)=-\frac{1}{p\left(t\right)},q\left(t\right)=\frac{\Gamma \left(t\right)}{p\left(t\right)}.$ (3.8)

将(3.7)式除以$p\left(t\right)$ ，并令$C=\xi \left(t_0\right)$ ，结合(3.8)式可以得到

$\frac{1}{p\left(t\right)}h\left(t\right)+q\left(t\right)\dot{h}\left(t\right)\le \frac{C}{p\left(t\right)},\forall t\in \left[t_0,+\infty \right).$

进一步整理得

$q\left(t\right)\dot{h}\left(t\right)-\dot{q}\left(t\right)\left(h\left(t\right)-C\right)\le 0,\forall t\in \left[t_0,+\infty \right),$

其中，$h\left(t\right)$ 为锚函数。

对上式除以$q{\left(t\right)}^2$ 得到

$\frac{1}{q{\left(t\right)}^2}\left[q\left(t\right)\dot{h}\left(t\right)-\dot{q}\left(t\right)\left(h\left(t\right)-C\right)\right]=\frac{\text{d}}{\text{d}t}\left(\frac{h\left(t\right)-C}{q\left(t\right)}\right)\le 0,\forall t\in \left[t_0,+\infty \right).$

对上述不等式$\frac{\text{d}}{\text{d}t}\left(\frac{h\left(t\right)-C}{q\left(t\right)}\right)\le 0$ 积分，则得$h\left(t\right)\le C_1\left(1+q\left(t\right)\right)$ 及$C_1>0$ 。根据锚函数$h\left(t\right)$ 的定义(3.5)，因此二阶微分方程系统(3.1)的解的轨迹是有界的，定理3.1得证。

4. 数值实验

Sun等[15]中建立了一阶微分方程系统求解一类互补问题，现将该一阶微分方程系统与本文的二阶微分方程系统所需的条件进行对比。其一阶微分方程系统如下

$\{\begin{array}{l}\frac{\text{d}z\left(t\right)}{\text{d}t}=-\rho \nabla \Phi \left(z\left(t\right)\right),\hfill \\ z\left(t_0\right)=z_0,\hfill \end{array}$ (4.1)

其中，$\rho >0$ 是该模型中的缩放因子。对于系统稳定性的研究中，Sun等[15]要求$\nabla S\left(z\right)$ 为非奇异的，则一阶微分方程系统(4.1)的平衡点为指数稳定的。下面对一阶微分方程系统(4.1)与二阶微分方程系统(3.1)的理论条件作对比。具体如下表1：

Table 1. Comparison of theoretical conditions

表1. 理论条件对比

理论条件	一阶微分方程系统(1.8)	二阶微分方程系统(3.1)
解集$\arg \min \Phi$ 非空	√	√
$\nabla S\left(z\right)$ 的非奇异性	√
$JF\left(x^{*}\right)$ 为半正定的	√
$\arg \min \Phi$ 是紧的	√
$\gamma \left(t\right)$ ，$\beta \left(t\right)$ 满足假设条件$H_{\gamma ,\beta }$		√

下面将用二阶微分方程系统(3.1)来求解两个具有不等式约束条件的变分不等式问题(1.1)，说明二阶微分方程方法的有效性，同时对该变分不等式也运用一阶微分方程系统(4.1)进行求解，将其结果进行对比。

例 4.1：考虑以下变分不等式

$\langle F\left(x^{*}\right),y-x^{*}\rangle \ge 0,\forall y\in K,$

其中

$K=\left\{x\in R^5|-g\left(x\right)=x\ge 0\right\}$

及

$F\left(x\right)=\left(\begin{array}{c}2x_1-4\\ {\text{e}}^{x_2}-1\\ 3x_3-4\\ -\sin \left(x_4\right)\\ x_5\end{array}\right),$

该问题的解为$x^{*}=\left(2,0,1.333,0,0\right)$ 。计算$Jg{\left(x\right)}^{\text{T}}=-I_5$ ，其中$I_5$ 表示5阶单位矩阵，其KKT条件如下

$\{\begin{array}{l}L\left(x,\lambda \right)=F\left(x\right)-\lambda I_5=0,\\ g\left(x\right)\perp \lambda .\end{array}$

下面写出其光滑化的函数为

$S\left(\epsilon ,x,\lambda \right)=\left(\begin{array}{c}\epsilon \\ L\left(x,\lambda \right)\\ {\phi }_{NR}^{\epsilon }\left(x,\lambda \right)\end{array}\right)=0,$

其中

${\phi }_{NR}^{\epsilon }\left(x,\lambda \right)=\left(\begin{array}{c}{x}_1-\frac{1}{2}\left(x_1-{\lambda }_1+\sqrt{{\epsilon }^2+{\left(x_1-{\lambda }_1\right)}^2}\right)\\ x_2-\frac{1}{2}\left(x_2-{\lambda }_2+\sqrt{{\epsilon }^2+{\left(x_2-{\lambda }_2\right)}^2}\right)\\ x_3-\frac{1}{2}\left(x_3-{\lambda }_3+\sqrt{{\epsilon }^2+{\left(x_3-{\lambda }_3\right)}^2}\right)\\ x_4-\frac{1}{2}\left(x_4-{\lambda }_4+\sqrt{{\epsilon }^2+{\left(x_4-{\lambda }_4\right)}^2}\right)\\ x_5-\frac{1}{2}\left(x_5-{\lambda }_5+\sqrt{{\epsilon }^2+{\left(x_5-{\lambda }_5\right)}^2}\right)\end{array}\right).$

则$\Phi \left(\epsilon ,x.\lambda \right)=\frac{1}{2}{\Vert S\left(\epsilon ,x,\lambda \right)\Vert }^2$ 。运用Matlab软件对微分方程系统(3.1)和(4.1)求解。如下表2。

例 4.2：考虑以下变分不等式

$\langle F\left(x^{*}\right),y-x^{*}\rangle \ge 0,\forall y\in K,$

其中

$K=\left\{x\in R^3|x_1^2+0.4x_2^2+0.6x_3^2-1\le 0\right\}$

Table 2. Numerical results of Example 4.1

表2. 例4.1的数值结果

	一阶微分方程系统(4.1)	二阶微分方程系统(3.1)
CPU时间(秒)	0.26187	1.46514
$x^{\ast }$	$\left(\begin{array}{c}2\\ 1.1078\text{e}-16\\ 1.33335\\ 1.2608\text{e}-17\\ 1.0759\text{e}-16\end{array}\right)$	$\left(\begin{array}{c}2\\ 3.4092\text{e}-5\\ 1.3333\\ 4.0908\text{e}-6\\ 3.5492\text{e}-5\end{array}\right)$

及

$F\left(x\right)=\left(\begin{array}{c}2x_1+0.2x_1^3-0.5x_2+0.1x_3-4\\ -0.5x_1+x{}_2+0.1x_2^3+0.5\\ 0.5x_1-0.2x_2+2x_3-0.5\end{array}\right).$

该问题的解$x^{*}={\left(1,0,0\right)}^{\text{T}}$ 。计算$Jg{\left(x\right)}^{\text{T}}\lambda =\left(\begin{array}{c}2x_1\\ 0.8x_2\\ 1.8x_3\end{array}\right)\lambda$ ，其KKT条件如下

$\{\begin{array}{l}L\left(x,\lambda \right)=F\left(x\right)-\left(\begin{array}{c}2x_1\\ 0.8x_2\\ 1.8x_3\end{array}\right)\lambda =0,\\ g\left(x\right)\perp \lambda .\end{array}$

其光滑化函数为

$S\left(\epsilon ,x,\lambda \right)=\left(\begin{array}{c}\epsilon \\ L\left(x,\lambda \right)\\ {\phi }_{NR}^{\epsilon }\left(-g\left(x\right),\lambda \right)\end{array}\right)=0,$

其中

${\phi }_{NR}^{\epsilon }\left(-g\left(x\right),\lambda \right)=-g\left(x\right)-\frac{1}{2}\left(-g\left(x\right)-\lambda +\sqrt{{\epsilon }^2+{\left(-g\left(x\right)-\lambda \right)}^2}\right).$

则$\Phi \left(\epsilon ,x,\lambda \right)=\frac{1}{2}{\Vert S\left(\epsilon ,x,\lambda \right)\Vert }^2$ 。运用Matlab软件对微分方程系统(3.1)和(4.1)求解。如下表3。

Table 3. The numerical results of Example 4.2

表3. 例4.2的数值结果

	一阶微分方程系统(4.1)	二阶微分方程系统(3.1)
CPU时间(秒)	0.0687	0.5943
$x^{*}$	$\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)$	$\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right)$

5. 小结

运用Attouch等学者的微分方程系统求解优化问题的思想，本文研究了具有不等式约束的变分不等式(1.1)的KKT解的收敛性问题，这一方法的运用丰富了变分不等式的微分方程方法的研究。首先运用光滑化的自然残差函数建立了具有不等式约束条件的变分不等式问题(1.1)的KKT条件的光滑化方程，并建立了无约束优化问题。建立了具有阻尼惯性参数和时间尺度参数的二阶微分方程系统，并研究了该二阶微分方程系统的稳定性定理，从而得到了具有不等式约束条件的变分不等式问题(1.1)的KKT点的收敛性。并与已有的一阶微分方程系统进行了理论条件和数值结果的对比，在理论条件上，二阶微分方程系统要更容易实现。但是在数值计算上，一阶微分方程系统要更快一些，但是其差值不是特别大，可忽略不计。

项目支持

沈阳航空航天大学2024年大学生创新创业项目，编号X202410143190。