Second-order differential equation method for solving stochastic variational inequality

Huiting ZHUANG; Li WANG; Juhe SUN; Danna JIA; Yanhong YUAN

doi:10.3969/j.issn.2095-1248.2023.04.012

Journal of Shenyang Aerospace University >

2023 , Vol. 40 >Issue 4: 88 - 96

DOI: https://doi.org/10.3969/j.issn.2095-1248.2023.04.012

Fundamental Science and Engineering

Second-order differential equation method for solving stochastic variational inequality

Huiting ZHUANG ^,¹ ,
Li WANG ^,¹ ,
Juhe SUN ¹ ,
Danna JIA ¹ ,
Yanhong YUAN ²

Expand

^1. College of Sciences，Shenyang Aerospace University，Shenyang 110136，China
^2. College of Economics and Management，Taiyuan University of Technology，Taiyuan 030024，China

Received date: 2023-03-31

Online published: 2023-11-09

Fold

Abstract

A system of second-order differential equation with positive viscous damping coefficients and time-scale coefficients was applied to solve the stochastic variational inequality problem （SVIP）. Firstly， the complementary function and the sample average approximation （SAA） method were applied to equate the original problem， and the stochastic variational inequality problem was transformed into a system of equations. Based on this， a second-order differential equation system with positive viscous damping coefficients $γ t$ and time-scale coefficients $β t$ was established. Secondly， the convergence and convergence rate of the trajectory of the second-order differential equation system were obtained. Finally， two numerical experiments were presented to demonstrate the effectiveness of the second-order differential equation system in solving stochastic variational inequality problems.

Key words： stochastic variational inequality; second-order differential equation; complementary function; sample average approximation method; convex optimization problem

Cite this article

Huiting ZHUANG , Li WANG , Juhe SUN , Danna JIA , Yanhong YUAN . Second-order differential equation method for solving stochastic variational inequality[J]. Journal of Shenyang Aerospace University, 2023 , 40(4) : 88 -96 . DOI: 10.3969/j.issn.2095-1248.2023.04.012

变分不等式是数学领域的一个重要分支。经过几十年的研究与发展，变分不等式的理论和算法研究已经趋于成熟，并在多个学科中有所应用。变分不等式的研究主要集中于解的存在条件^［1-2］、迭代算法及收敛条件^［3-5］，并在经济、交通、工程等领域有广泛应用^［6-8］。基于此，该问题得到了越来越多研究者的关注。

Ferguson等^［9］在1956年提出了随机规划问题，而随机变分不等式问题作为随机规划的一个重要分支，在经济分析、交通均衡、能源建模等领域^［10-11］都有着广泛的应用，更多关于随机规划知识的介绍见文献［12］。基于此，随机变分不等式问题成为近年来的研究热点之一^［13-15］。本文将研究随机变分不等式问题，求解

x ∈ C

，使得

E f x, ξ, y - x ≥ 0

，

∀ y ∈ C

（1）

式中：

E f x, ξ = ∫ Ω f x, ξ d P ξ

是随机变量

ξ

的数学期望。约束集合

C

如式（2）所示

C = x ∈ R n E h x, η = 0, - E g x, ς ≥ 0

（2）

式中：

R n

为

n

维实数空间；

f ⋅, ξ

对任意

ξ

是局部Lipschitz连续。对每一个固定的

η

和

ς

，函数

h ⋅, η : R n × Ω → R l

和

g ⋅, ς : R n × Ω

→ R m

均是连续可微的凸函数。

ξ

、

η

和

ς

是定义在概率空间

Ω, F, P

上的随机变量，

F

是Sigma代数，

P

是定义在

Ω

上的概率测度。

显然，随机变分不等式问题中存在着随机变量，即该问题中可能存在不确定因素。随着时代的发展，人们越来越重视对实际问题中一些不确定性风险的预测和处理，特别是在经济金融、工程应用和决策科学等领域中所涉及的很多因素都是不确定的，因此研究包含的随机因素优化问题更符合实际意义。近年来，随机变分不等式作为处理不确定环境下最优化问题和博弈论的重要方法之一，受到了广泛的关注和研究。

Facchinei等^［16］对变分不等式问题和互补问题的基本理论和应用进行了全面的介绍和总结，该书被认为是最优化领域中最有价值的工作，记录了变分不等式和互补问题很大一部分重要结果，但是却没有记录微分方程方法求解变分不等式的理论。因此，应用微分方程方法求解变分不等式这一问题的研究还远未深入，这部分内容还需要进一步补充，这也是本文将要讨论的主要问题。

受到Attouch等^［17］的启发，利用二阶微分方程系统求解随机变分不等式问题（1）是一个新的思想。在上述研究成果的基础之上，本文将研究求解随机变分不等式问题（1）的二阶微分方程方法。许多文献的数值结果表明SAA方法是一种有效的求解随机问题的方法^［18-19］，所以本文将运用互补函数和SAA方法将原始问题等价转换成方程组，利用该方程组建立具有正黏性阻尼系数

γ t

和时间尺度系数

β t

的二阶微分方程系统。然后进一步研究该二阶微分方程系统的收敛性和收敛速度。最后，给出两个数值实验说明该二阶微分方程系统求解随机变分不等式问题的有效性。

为了建立随机变分不等式的二阶微分方程系统，首先介绍下面的概念与相关性质。

定理1 ^［17］设

ℋ

是希尔伯特（Hilbert）空间，

Φ : ℋ → R

是凸的连续可微的函数，其梯度

∇ Φ

是Lipschitz连续的。

若

β, γ : t 0, + ∞ → R +

，且

g : t 0, + ∞ → ℋ

是局部可积的，则二阶微分方程系统

x ¨ t + γ t x ˙ t + β t ∇ Φ x t = g t

在初始条件

x t 0, x ˙ t 0 ∈ ℋ × ℋ

下，存在唯一的全局解

x ∈ t 0, + ∞ → ℋ

。

定义1 ^［17］如果映射

ϕ : R n × R n → R n

，满足

x ∈ R n, y ∈ R n, x T y = 0 ⇔ ϕ x, y = 0,

则称

ϕ x, y

为

R n

上的互补函数。特别地，NR函数

φ ε, u = 12 u + ε 2 e + u 2

就是一类互补函数，将NR函数光滑化后可得光滑的NR函数为

ϕ N R ε x, y = x - φ ε, x - y

。

1 随机变分不等式的转化

由优化理论不难得到SVIP问题可以等价为

- E f x, ξ ≥ 0

。因此，SVIP问题的KKT条件为

L x, μ, λ = 0, E h x, η = 0, E g x, ς T λ = 0, 且 λ, - E g x, ς ≥ 0,

（3）

L x, μ, λ = E f x, ξ + J x E h x, η T μ

+ J x E g x, ς T λ

（

J x F x, y T

表示对函数

F x, y

中的变量

x

所求得的Jacobian矩阵的转置表示SVIP问题的拉格朗日函数。由SAA方法，问题（1）可以转化为如下问题：找到

x ∈ C N

，使得

F N x, y - x ≥ 0

，

∀ y ∈ C N

（4）

其中约束集合

C N

表示为

C N = x ∈ R N h N x = 0, - g N x ≥ 0

（5）

F N = 1 N ∑ i = 1 N f x, ξ i

，

h N = 1 N ∑ i = 1 N h x, η i

，

g N = 1 N ∑ i = 1 N g x, ς i

，

ξ i

、

η i

和

ς i

分别是随机变量

ξ

、

η

和

ς

的

N

个独立同分布的样本。本文称问题（4）为SAA问题，问题（1）为原问题。类似地，可以得到SAA问题（4）的KKT条件为

L N x, μ, λ = 0, h N x = 0, g N x T λ = 0, 且 λ, - g N x ≥ 0

（6）

其中

L N x, μ, λ = F N x + J x h N x T μ + J x g N x T

根据互补函数的定义1，设

S : R × R n × R l × R m → R × R n × R l × R m

，

z = ε, x, μ, λ ∈ R × R n × R l × R m

，则SAA问题的KKT条件（6）可以转化为如下的光滑化方程

S z = ε L N x, μ, λ h N x - g N x - φ ε, - g N x - λ = 0

（7）

则有

∇ S z = ∇ ε S z T ∇ x S z T ∇ μ S z T ∇ λ S z T

（8）

运用光滑化的NR函数和SAA方法将原SVIP问题（1）转化为光滑化方程组（7），下一节将运用光滑方程组（7）建立二阶微分方程系统。

2 随机变分不等式的二阶微分方程方法

由于原SVIP问题（1）转化为光滑化方程组（7），为了求解方程组（7）的值，定义函数

Φ z = 12 S z 2

，考虑下面的凸优化问题

m i n z ∈ R 1 + n + l + m Φ z = m i n z ∈ R 1 + n + l + m 12 S z 2

（9）

则有

∇ Φ z = ∇ S z T S z = ∇ ε S z T S z ∇ x S z T S z ∇ μ S z T S z ∇ λ S z T S z

（10）

首先假设以下条件成立：

（1）

Φ z = 12 S z 2

为凸函数，

P = a r g m i n Φ

非空；（2）

γ

和

β : t 0, + ∞ → R +

为非负连续函数；（3）

t 0 > 0

是初始时间。建立二阶微分方程系统如下

ε ¨ t x ¨ t μ ¨ t λ ¨ t + γ t ε ˙ t x ˙ t μ ˙ t λ ˙ t + β t ∇ ε S z T S z ∇ x S z T S z ∇ μ S z T S z ∇ λ S z T S z = 0

（11）

式中：

∇ Φ : R → R

为

Φ

的梯度，还有两个随时间变化的参数分别为：

γ t

是正黏性阻尼系数，

β t

是时间尺度系数。根据定理1可知，该方程存在唯一的全局解。在平凡情况

Φ z = 12 S z 2 ≡ 0

时，上述二阶微分方程系统（11）可直接积分得到

p t = e ∫ t 0 t γ u d u

（12）

设二阶微分方程系统（11）满足以下假设条件

H 0

：

∫ t 0 + ∞ d u p u < + ∞

在该假设条件下，可以定义函数

Γ : t 0, + ∞ → R +

，

Γ t = p t ∫ t + ∞ d u p u

Γ t

的定义不依赖于初始时间

t 0

的选择，当

t 0 → + ∞

时，可以得到

Γ ˙ t = γ t Γ t - 1

。

定义全局能量函数如下

W t = 12 ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 + β t

12 S z t 2 - m i n z ∈ R 1 + n + l + m 12 S z 2

（13）

以及锚点函数

h t = 12 ε t x t μ t λ t - ε * x * μ * λ * 2

（14）

其中

ε * x * μ * λ * ∈ a r g m i n 12 S z 2

有了上述三类函数，下面定义函数

χ : t 0, + ∞ → R +

，该函数在二阶微分方程系统（11）的收敛性定理的证明过程中起着重要的作用。具体形式如下

χ t = Γ t 2 W t + h t + Γ t h ˙ t

=

Γ t 2 12 ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 + 12 ε t x t μ t λ t - ε * x * μ * λ * 2

+

Γ t 2 β t 12 S z t 2 - m i n z ∈ R 1 + n + l + m 12 S z 2

+

Γ t ε ˙ t x ˙ t μ ˙ t λ ˙ t, ε t x t μ t λ t - ε * x * μ * λ *

根据内积与范数的关系，上式的右边可以化为

Γ t 2 β t 12 S z t 2 - m i n R 1 + n + l + m 12 S z 2

+

12 ε t x t μ t λ t - ε * x * μ * λ * + Γ t ε ˙ t x ˙ t μ ˙ t λ ˙ t 2

显然，

χ ⋅

是一个非负函数且是递减的。运用函数

χ ⋅

可以得到二阶微分方程系统式（11）的收敛性定理，该定理也给出了收敛率的估计。

定理2 设

Φ z = 12 S z 2

为凸函数，

P = a r g m i n Φ

非空，

β, γ : t 0, + ∞ → R +

为非负连续函数，且满足假设条件

H 0

，并满足以下增长条件

H γ, β

H γ, β

Γ t β ˙ t ≤ β t 3 - 2 γ t Γ t

则对于二阶微分方程系统（11）的轨迹

z t = ε t, x t, μ t, λ t ∈ R × R n × R l × R m

满足以下收敛性

12 S z t 2 - m i n z ∈ R 1 + n + l + m 12 S z 2

= O 1 β t Γ t 2, t → + ∞

此外，可以得到该轨迹

z t = ε t, x t, μ t, λ t

在

t 0, + ∞

上是有界的。

证明：设

m = m i n z ∈ R 1 + n + l + m 12 S z 2

，为了计算

χ ⋅

的导数，首先计算能量函数

W

和锚点函数

h

的导数

W ˙ t = ε ˙ t x ˙ t μ ˙ t λ ˙ t, ε ¨ t x ¨ t μ ¨ t λ ¨ t +

β t ε ˙ t x ˙ t μ ˙ t λ ˙ t, ∇ ε S z t T S z t ∇ x S z t T S z t ∇ μ S z t T S z t ∇ λ S z t T S z t

+ β ˙ t 12 S z t 2 - m

= - γ t ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 + β ˙ t 12 S z t 2 - m

和

h ˙ t = ε ˙ t x ˙ t μ ˙ t λ ˙ t, ε t x t μ t λ t - ε * x * μ * λ *

，

h ¨ t = ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 + ε ¨ t x ¨ t μ ¨ t λ ¨ t, ε t x t μ t λ t - ε * x * μ * λ *

所以有

h ¨ t + γ t h ˙ t = ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 +

ε ¨ t x ¨ t μ ¨ t λ ¨ t + γ t ε ˙ t x ˙ t μ ˙ t λ ˙ t, ε t x t μ t λ t - ε * x * μ * λ *

≤ ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 - β ˙ t 12 S z t 2 - m

≤ ε ˙ t x ˙ t μ ˙ t λ ˙ t 2

其中，上述不等式是由

12 S z t 2

的凸性得到的，综合上述结果可以得到

χ ˙ t = Γ t 2 W ˙ t + 2 Γ t Γ ˙ t W t

+ h ˙ t + Γ ˙ t h ˙ t + Γ t h ¨ t

= Γ t 2 - γ t ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 + β ˙ t 12 S z t 2 - m +

2 Γ t Γ ˙ t 12 ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 + β t 12 S z t 2 - m

+ h ˙ t + Γ ˙ t h ˙ t + Γ t h ¨ t

根据

1 + Γ ˙ t = γ t Γ t

，可以得到

h ˙ t + Γ ˙ t h ˙ t + Γ t h ¨ t

= Γ t h ¨ t + γ t h ˙ t ≤

Γ t ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 - β t 12 S z t 2 - m

结合以上结果，可以得到

χ ˙ t ≤ - Γ t 2 γ t ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 +

Γ t 2 β ˙ t 12 S z t 2 - m +

2 Γ t Γ ˙ t β t 12 S z t 2 - m +

Γ t Γ ˙ t ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 + Γ t ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 -

Γ t β t 12 S z t 2 - m

≤ Γ t ε ˙ t x ˙ t μ ˙ t λ ˙ t 2 - β t 12 S z t 2 - m

+ 12 S z t 2 - m Γ t

Γ t β ˙ t + 2 Γ ˙ t β t - β t

由

1 + Γ ˙ t = γ t Γ t

可得

χ ˙ t ≤ 12 S z t 2 - m Γ t

Γ t β ˙ t + 2 Γ ˙ t β t - β t

与

χ ˙ t ≤ 12 S z t 2 - m Γ t

Γ t β ˙ t + β t 2 γ t Γ t - 3

等价。

因此，由假设

H γ, β

可以推断出

χ ˙ t ≤ 0

。所以在

t 0, + ∞

上有

χ t ≤ χ t 0

，根据

χ t

的公式可以得到对于所有的

t > t 0

有

12 S z t 2 - m i n z ∈ R 1 + n + l + m 12 S z 2 ≤ χ t 0 β t Γ t 2

（15）

此外，证明

z t = ε t, x t, μ t, λ t

在

t 0, + ∞

上仍为有界的，由

χ ⋅

的公式和递减性，可得

ε t x t μ t λ t - ε * x * μ * λ * + Γ t ε ˙ t x ˙ t μ ˙ t λ ˙ t 2

≤ 2 χ t ≤ 2 χ t 0

（16）

由上述不等式可以得到

ε t x t μ t λ t - ε * x * μ * λ * 2

+

2 Γ t ε t x t μ t λ t - ε * x * μ * λ *, ε ˙ t x ˙ t μ ˙ t λ ˙ t ≤ 2 χ t 0

（17）

设

h t = 12 ε t x t μ t λ t - ε * x * μ * λ * 2

，且由

q t = ∫ t + ∞ d s p s

可得

q t t ≥ t 0

是有界的，因为

∫ t + ∞ d s p s < + ∞

。将（17）除以

p t

，又因为

Γ t p t = q t

，且

C = χ t 0

，可得

1 p t h t + q t h ˙ t ≤ C p t, ∀ t ∈ t 0, + ∞

。

因为

q ˙ t = - 1 p t

，所以有

q t h ˙ t - q ˙ t h t -

C ≤ 0

，再除以

q t 2

后，可得

d d t h t - C q t =

1 q t 2 q t h ˙ t - q ˙ t h t - C ≤ 0

上式积分后可得：对于某些

c 1 > 0

，有

h t ≤ c 1 1 + q t

。

因此

z t = ε t, x t, μ t, λ t

在

t 0, + ∞

上有界，证毕。

3 数值实验

数值实验采用Matlab9.0编写，常微分方程求解采用龙格-库塔ode 45。在实验中，选取系数

γ t = α t

，

α = 25

和

β t ≡ 1

。

例1 考虑随机变分不等式问题（1），其中

F x, ξ = ξ x 1 - 3 2, s i n x 2, x 3 - 1 2

和

C = x ∈ R 3 - g x, ς = ς x ≥ 0

式中：

ξ

服从均值为1、方差为1的正态分布；

ς

是

0,1

上的随机分布；该问题的最优解为

x * = 3.0000,1.5708,1.0000

。本实验中取初始点分别为

ε 0 = ε ¯ = 0.0001

，

x 0 = 3,1, 1

，

λ 0 = - 0.5, - 0.1, - 0.1

。

表1 总结了运用二阶微分方程系统（11）求解例1的数值结果。图1中给出了从所给的初始点出发，二阶微分方程系统（11）求解例1所得到的解的轨迹 $x t$ 。

时间/ms
50	0.5	1
100	0.25	1
150	0.167	1

例2 考虑随机变分不等式问题（1），其中

F x, ξ = ξ x 1 - 5 2, ξ s i n x 2, ξ x 3 - 1 2, x 4 2,

x 5 - 1 2

和

C = x ∈ R 5 - g x = ς x ≥ 0

，

式中：

ξ

服从均值为1，方差为1的正态分布；

ς

是

0,1

上的随机分布；该问题的最优解为

x * = 5,1.57,1, 0,1

。本实验中取初始点分别为

ε 0 = ε ¯ = 0.000 1

，

x 0 = 3,1, 1,0, 0

，

λ 0 = - 0.5, - 0.1,

- 0.1,0, 0

。

在表2中总结了运用二阶微分方程系统（11）求解例2的数值结果。图2中给出了从所给的初始点出发，二阶微分方程系统（11）求解例2所得到的解的轨迹

x t

。

表2 例2的二阶微分方程系统（11）的数值结果

时间/ms	γ（t）	β（t）
50	0.5	1
100	0.25	1
150	0.167	1

显示原图|下载原图ZIP|生成PPT

图2 例2的二阶微分方程系统（11）从初始点出发的 x （ t ）轨迹曲线

4 结论

本文首次采用Attouch等 ^［17］提出的基于惯性阻尼系数

γ t

和时间标度系数

β t

的二阶微分方程系统来求解随机变分不等式问题（1），得到了该二阶微分方程系统轨迹的收敛性和收敛速度，并研究了在

γ t

和

β t

影响下解的收敛性质。除了理论结果外，本文还给出两个数值算例说明二阶微分方程系统在求解随机变分不等式问题中的有效性。

References

Publishing order | Descend order by publishing year | Descend order by cited within

1	Thong D V， Hieu D V.Weak and strong convergence theorems for variational inequality problems ［J］.Numerical Algorithms，2018，78（4）：1045-1060.

2	Wu Z L， Lu Y.Minimum and maximum principle sufficiency for a nonsmooth variational inequality ［J］.Bulletin of the Malaysian Mathematical Sciences Society，2021，44（3）：1233-1257.

3	Nnakwe M O.An algorithm for approximating a common solution of variational inequality and convex minimization problems［J］.Optimization，2021，70（10）：2227-2246.

4	Jolaoso L O， Taiwo A， Alakoya T O，et al.A unified algorithm for solving variational inequality and fixed point problems with application to the split equality problem ［J］.Computational and Applied Mathematics，2020，39（1）：1-28.

5	Eslamian M.Variational inequality over the set of common solutions of a system of bilevel variational inequality problem with applcations ［J］.Revista de la Real Academia de Ciencias Exactas，Físicas y Naturales.Serie A.Matemáticas，2022，116（1）：1-18.

6	Parise F， Ozdaglar A.A variational inequality framework for network games：existence，uniqueness，convergence and sensitivity analysis ［J］.Games Economic Behavior，2019，114（13）：47-82.

7	杨振平.随机变分不等式问题及其在天然气市场中的应用研究［D］.上海：上海大学，2019.

8	陈琦琼.不确定变分不等式及其在非合作博弈中的应用［D］.南京：南京理工大学，2018.

9	Ferguson A R， Dantzig G B.The allocation of aircraft to routes-an example of linear programming under uncertain demand ［J］.Management Science，1956，3（1）：45-73.

10	R.Tyrrell Rockafellar and Roger J-B Wets.Stochastic variational inequalities：single-stage to multistage ［J］.Mathematical Programming，2017，165（1）：331-360.

11	Harker P T.A variational inequality approach for the determination of oligopolistic market equilibrium ［J］.Mathematical Programming，1984，30（1）：105-111.

12	Shapiro A， Dentcheva D， Ruszczynski A.Lectures on Stochastic Programming：Modeling and Theory ［M］.Philadeslphia：Society for Industrial and Applied Mathematics，2021.

13	Zhao Y， Zhang J， Yang X M，et al.Expected residual minimization formulation for a class of stochastic vector variational inequalities ［J］.Journal of Optimization Theory and Applications，2017，175（2）：545-566.

14	Chen S， Pang L P， Ma X F，et al.SAA method based on modified Newton method for stochastic variational inequality with second-order cone constraints and application in portfolio optimization ［J］.Mathematical Methods of Operations Research，2016，84（1）：129-154.

15	Sun J H， Chen J S， Ko C H.Neural networks for solving second-order cone constrained variational inequality problem ［J］.Computational Optimization and Applications，2010，51（2）：623-648.

16	Facchinei F， Pang J S.Finite-Dimensional Variational Inequalities and Complementarity Problems［M］.New York：Springer-Verlag，2003.

17	Attouch H， Chbani Z， Riahi H，Fast convex optimization via a third-order in time evolution equation ［J］.Optimization，2020，71（5）：1275-1304.

18	Wang C M， He S X， Wu H Y.An implementable SAA nonlinear lagrange algorithm for constrained minimax stochastic optimization problems［J］.Mathematical Problems in Engineering，2018（16）：1-13.

19	Lam H， Jiang G X， Fu M C，et al.On efficiencies of stochastic optimization procedures under importance sampling ［C］ //Winter Simulation Conference Proceedings.Gothenburg，Sweden，2018，1862-1873.

Options

Outlines

模态框（Modal）标题