# Differential variational inequalities

Jong-Shi Pang, David E. Stewart

► To cite this version:

Jong-Shi Pang, David E. Stewart. Differential variational inequalities. Mathematical Programming, 2008, 113 (2), 10.1007/s10107-006-0052-x. hal-01366027

HAL Id: hal-01366027

https://hal.science/hal-01366027

Submitted on 14 Sep 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Differential variational inequalities

Jong-Shi Pang · David E. Stewart

**Abstract** This paper introduces and studies the class of *differential variational inequalities* (DVIs) in a finite-dimensional Euclidean space. The DVI provides a powerful modeling paradigm for many applied problems in which dynamics, inequalities, and discontinuities are present; examples of such problems include constrained time-dependent physical systems with unilateral constraints, differential Nash games, and hybrid engineering systems with variable structures. The DVI unifies several mathematical problem classes that include ordinary differential equations (ODEs) with smooth and discontinuous right-hand sides, differential algebraic equations (DAEs), dynamic complementarity systems, and evolutionary variational inequalities. Conditions are presented under which the DVI can be converted, either locally or globally, to an equivalent ODE with a Lipschitz continuous right-hand function. For DVIs that cannot be so converted, we consider their numerical resolution via an Euler time-stepping procedure, which involves the solution of a sequence of finite-dimensional

<footnote>The work of J.-S. Pang is supported by the National Science Foundation under grants CCR-0098013 CCR-0353074, and DMS-0508986, by a Focused Research Group Grant DMS-0139715 to the Johns Hopkins University and DMS-0353016 to Rensselaer Polytechnic Institute, and by the Office of Naval Research under grant N00014-02-1-0286.</footnote>

<footnote>The work of D. E. Stewart is supported by the National Science Foundation under a Focused Research Group grant DMS-0138708.</footnote>

<footnote>J.-S. Pang (✉)</footnote>

<footnote>Department of Mathematical Sciences and Department of Decision Science and Engineering Systems, Rensselaer Polytechnic Institute, Troy, NY 12180-3590, USA</footnote>

<footnote>e-mail: pangj@rpi.edu</footnote>

<footnote>D. E. Stewart</footnote>

<footnote>Department of Mathematics, University of Iowa, Iowa City, IA 52242, USA</footnote>

<footnote>e-mail: dstewart@math.uiowa.edu</footnote>
variational inequalities. Borrowing results from differential inclusions (DIs) with upper semicontinuous, closed and convex valued multifunctions, we establish the convergence of such a procedure for solving initial-value DVIs. We also present a class of DVIs for which the theory of DIs is not directly applicable, and yet similar convergence can be established. Finally, we extend the method to a boundary-value DVI and provide conditions for the convergence of the method. The results in this paper pertain exclusively to systems with “index” not exceeding two and which have absolutely continuous solutions.

# 1 Introduction

Ordinary differential equations (ODEs) with smooth input functions are a classical paradigm in applied mathematics that have existed for centuries. In the past couple decades, differential algebraic equations (DAEs) have become a very important generalization of ODEs and have been studied extensively; see the two excellent monographs [19,77] and the references therein. Yet, as evidenced by the growing literature that has surfaced in recent years on multi-rigid-body dynamics with frictional contacts [3–6,17,18,90,91,94,95,98,101,102] and on hybrid engineering systems [21,49,51,52], ODEs and DAEs are vastly inadequate to deal with many naturally occurring engineering problems that contain inequalities (for modeling unilateral constraints) and disjunctive conditions (for conditional phenomena such as contacts and mode switching). Built on the “static” finite-dimensional variational inequalities (VIs), for which a comprehensive reference is available [38], differential variational inequalities (DVIs), which include differential complementarity problems (DCPs), provide a fundamental generalization of DAEs that significantly extends these differential equations and opens up a broad paradigm for the enhanced modeling of complex engineering systems.

This paper has several objectives: to formally introduce the DVI and its special cases, to discuss several source problems in order to demonstrate the usefulness of the DVI as a unifying mathematical framework for a host of dynamic problems involving inequalities and discontinuities, to present solution existence results, and to investigate the numerical resolution of the DVI via a Euler time-stepping procedure. Two sets of mathematical theories provide the cornerstone for our development: differential inclusions (also known as multi-valued differential equations) [10,32,66,89] and finite-dimensional variational inequalities (VIs) and complementarity problems (CPs) [38].

As we see later, a DVI can be looked upon as a special case of a DAE (albeit with nonsmooth input functions) and also of a differential inclusion (DI). This raises the question of why it is of interest to identify the DVI as a separate problem class and give it an independent treatment. There are two principal reasons. One is the fact that many applied problems naturally lend themselves to a DVI whose special structure has only been minimally taken into account in prior studies of the DAE and the DI. To appreciate the importance of exploiting such structure one can recall the relationship between the finite-dimensional
VI and Robinson's framework of a generalized equation [80,81,83]; although the former has greatly benefited from the latter in terms of its theoretical developments, when it comes to practical computations, the rich structure of the former enables the design of a host of efficient algorithms that have no parallels for the latter; see [38, Chaps. 9–11]. To this end, we should mention the second author's paper [93], which is arguably the first work on using linear complementarity theory to study ODEs with discontinuous but "piecewise smooth" right-hand functions. The other reason for initiating a focused study on the DVI is that although the DVI can be cast as a DAE, which is a practical framework to develop computational algorithms, in such a formulation, one has to employ nonsmooth functions that lead to a differential system with a nondifferentiable algebraic equation. This immediately invalidates the established methodology of DAEs. In other words, the formulation of a DVI as a DAE is, at the present time, only of conceptual significance. In contrast, the theory of DIs with closed convex-valued multifunctions is indeed applicable to some classes of DVIs; yet, the demonstration of the applicability of the DI theory is a nontrivial task and requires extensive results for finite-dimensional VIs/CPs. Part of the contribution of this paper is to demonstrate how such results can be profitably employed for the analysis of the DVI.

In summary, the DVI occupies a significant niche between the two domains: DAE and DI, and provides a fruitful modeling ground for many applied problems. Special attention to the DVI is therefore warranted. Our work presents a systematic investigation of the DVI as a general mathematical paradigm. By briefly mentioning a differential generalized Nash game (Sect. 4.2), which leads to a differential quasi-variational inequality, we are suggesting that there remains much to be researched; our goal is that this paper will provide the motivation and impetus for the sustained development of this important class of differential problems.

## 1.1 Structure of paper

In Sect. 2 we introduce the formulation of the DVI; our solution concept is formalized in Sect. 2.1. We outline a classification of DVIs in Sect. 2.2 which serves as a guide for much of the development in this paper. Connections with other differential problem classes are noted in Sect. 2.3. In Sect. 3, we introduce two general approaches for designing iterative methods for solving DVIs. Section 4 lists a number of areas in which DVIs are likely to arise. Section 5 deals with the first of the classifications of DVIs (index) which can be treated via the theory of Lipschitz ODEs. Section 6 deals with the next of the classifications of DVIs (multivalued index one DVIs), which can be treated via the theory of convex-valued differential inclusions satisfying a linear growth property. The section ends with a main result, Theorem 6.1, that summarizes various sufficient conditions under which a DVI has a "weak solution" in the sense of Carathéodory. The convergence of time-stepping methods is analyzed in detail in Sects. 7–9, the former two treating initial-value problems that
include a class of “mixed index one and index two” DVIs; Section 9 deals with boundary-value problems. Section 10 studies the linear-quadratic differential Nash game. The paper ends with some concluding remarks.

## 2 Formulation of DVIs

The DVI comprises two major components: an ODE and a VI. We begin with a formal definition of the latter. For a generic mapping $\Phi : \mathbb{R}^m \to \mathbb{R}^m$ and a nonempty closed convex set $K$ in $\mathbb{R}^m$, the VI $(K, \Phi)$, is to find a vector $u \in K$ such that

$$ (u' - u)^T \Phi(u) \geq 0, \quad \forall u' \in K. $$

Let SOL$(K, \Phi)$ denote the solution set of this problem. If $\Phi(u) \equiv r + Du$ is an affine function for some vector $r \in \mathbb{R}^m$ and matrix $D \in \mathbb{R}^{m \times m}$, we write VI $(K, r, D)$ for VI $(K, \Phi)$ and SOL$(K, r, D)$ for SOL$(K, \Phi)$. When $K$ is a polyhedron, we use the adjective “linearly constrained” to describe the VI $(K, \Phi)$; we use the acronym AVI to mean the *affine variational inequality* when $(K, \Phi)$ is an affine pair; i.e., $K$ is polyhedral and $\Phi$ is affine. We refer to [38] for a comprehensive study of the finite-dimensional variational inequality and the related complementarity problem (CP), where $K = C$ is a cone:

$$ C \ni u \perp \Phi(u) \in C^*, \quad (2.1) $$

with $C^* \equiv \{v \in \mathbb{R}^m : u^T v \ge 0, \forall u \in C\}$ being the dual cone of $C$. When $C$ is the nonnegative orthant $\mathbb{R}_+^m$, (2.1) becomes the *nonlinear complementarity problem* (NCP), which further reduces to the *linear complementarity problem* (LCP) [30] when $\Phi(u)$ is additionally an affine function. When $\Phi$ is a gradient map, say $\Phi \equiv \nabla \phi$, where $\phi$ is a real-valued function, then the VI $(K, \Phi)$ is the first-order stationarity condition for the optimization problem: minimize $\phi(u)$ subject to $u \in K$. Standard notation and known results from the LCP literature [30] will be employed; in particular, we write LCP $(r, D)$ and SOL$(r, D)$ for VI $(\mathbb{R}_+^m, r, D)$ and SOL$(\mathbb{R}_+^m, r, D)$, respectively. Finally, we write $\dot{x} \equiv \frac{dx}{dt}$ for the time-derivative of a function $x(t)$.

### 2.1 Formal definition

Let $f : \mathbb{R}^{1+n+m} \to \mathbb{R}^n$ and $F : \mathbb{R}^{1+n+m} \to \mathbb{R}^m$ be two continuous vector functions. Let $K$ be a nonempty closed convex subset of $\mathbb{R}^m$. Let $\Gamma : \mathbb{R}^{2n} \to \mathbb{R}^n$ be a *boundary function* and $T > 0$ be a *terminal time*. The DVI defined by the triplet of functions $f$, $F$, and $\Gamma$, the set $K$, and the scalar $T$ is to find time-dependent trajectories $x(t)$ and $u(t)$ that satisfy condition (2.2) in the *weak sense*
of Carathéodory for $t \in [0, T]$ and the algebraic condition (2.3):

$$ \begin{gathered} \dot{x}(t) = f(t, x(t), u(t)) \tag{2.2} \\ u(t) \in \text{SOL}(K, F(t, x(t), \cdot)) \end{gathered} $$

$$ \Gamma(x(0), x(T)) = 0. \tag{2.3} $$

This means that $x$ is an absolutely continuous function on $[0, T]$, that $u$ is an integrable function on $[0, T]$, and that the differential equation need only be satisfied for almost all $t$; moreover, the membership for $u(t)$ means that for any continuous $\tilde{u} : [0, T] \to K$,

$$ \int_0^T (\tilde{u}(t) - u(t))^{\mathrm{T}} F(t, x(t), u(t)) \mathrm{d}t \ge 0. \tag{2.4} $$

The latter implies that for almost all $t$, $u(t) \in \text{SOL}(K, F(t, x(t), \cdot))$. Conversely, if $u(t)$ is an integrable function satisfying the latter membership for all $t$, then the integral inequality (2.4) must hold for any continuous $\tilde{u}$. Expanding the set of functions $\tilde{u}$ under which (2.4) holds does not change the set of solutions. Consider the situation where $t \mapsto F(t, x(t), u(t))$ is in $L^p([0, T], \mathfrak{R}^m)$ for some $p' \ge 1$, and let $p$ satisfy $1/p + 1/p' = 1$. Since the space of continuous functions $C([0, T]; \mathfrak{R}^m)$ is dense in $L^p([0, T]; \mathfrak{R}^m)$, we can approximate any $\tilde{u}$ in the latter space by continuous $w_\epsilon$ in the former space; we can then set $\tilde{u}_\epsilon(t) = \Pi_K(w_\epsilon(t))$, where $\Pi_K$ is the nearest-point map to $K$ (i.e., the Euclidean projection onto $K$), in order to ensure that $\tilde{u}_\epsilon(t) \in K$ for all $t$. Thus if (2.4) holds for all continuous $\tilde{u}$, it also holds for all $\tilde{u} \in L^p$.

In fact, it is not necessary to require (2.4) to hold for all continuous $\tilde{u}$. We could require only, for example, that (2.4) holds for $\tilde{u} \in C^\infty$. Every continuous $\tilde{u}$ can be approximated by convolutions $\psi_\epsilon * \tilde{u}$ (if necessary, extend $\tilde{u}$ beyond $[0, T]$ with $\tilde{u}(t) = \tilde{u}(0)$ for $t < 0$ and $\tilde{u}(t) = \tilde{u}(T)$ for $t > T$) for $\psi_\epsilon(s) = \epsilon^{-1}\psi(s/\epsilon)$ with $\psi \ge 0$ in $C^\infty$, and $\int \psi(s) \mathrm{d}s = 1$. Since $K$ is convex, $\tilde{u}(t) \in K$ for all $t$ implies $\tilde{u}_\epsilon(t) \in K$ for all $t$. By uniform continuity of $\tilde{u}$ on $[0, T]$, we can show that $\tilde{u}_\epsilon \to \tilde{u}$ in the uniform norm. Hence if (2.4) holds for $\tilde{u}$ that are $C^\infty$, it holds for all continuous $\tilde{u}$, and even holds for all $\tilde{u} \in L^p$.

Under the absolute continuity of $x(t)$, the differential equation holding for almost all $t \in [0, T]$ is equivalent to the integral equation

$$ x(t) = x(s) + \int_s^t f(\tau, x(\tau), u(\tau)) \mathrm{d}\tau $$

holding for all $0 \le s \le t \le T$ [26, p. 42]. The weak sense of a solution to the DVI is generally needed in order to deal with possible discontinuities in $u(t)$ and the problem of convergence of time-stepping methods. This will become clear subsequently when we discuss the numerical methods.
The above general framework of the DVI includes a host of important special cases. In the stated form, the problem is a *two-point boundary-value problem* (BVP) in the sense that, linked by the abstract function $\Gamma$, both the initial state $x(0)$ and the terminal state $x(T)$ are unknown variables to be computed; in particular, the former variable $x(0)$ is not completely given. The *initial-value (IVP) version of the problem* corresponds to the special case where $\Gamma(x, y) \equiv x - x^0$, with $x^0$ being given. In the latter case, $x(0)$ is of course known. As it is well known in ODEs, in an *initial-value problem*, the terminal time $T$ may not always be strictly imposed, especially in situations where the input functions $f$ and $F$ are known to have favorable properties only in a neighborhood of the initial $x^0$; in such a case, one is interested in a solution that extends beyond the initial time $t=0$ but is not bound by the terminal time. Needless to say, an IVP is easier to analyze than a BVP, which is more subtle and involves considerably more topological considerations. In this paper, we will consider both types of problems, first the IVPs, followed by the BVPs.

For an initial-value problem corresponding to a given initial state $x^0$, some classes of DVIs have “locally smooth solutions” upon a reduction locally near a triple $(0, x^0, u^0)$, where $u^0 \in \text{SOL}(K, F(0, x^0, \cdot))$, to an equivalent ODE with a Lipschitz continuous right-hand side. Specifically, we say that the DVI (2.2) has a *locally smooth solution* near the triple $(0, x^0, u^0)$ if there exist a scalar $\delta > 0$, a continuously differentiable function $x : [0, \delta] \to \mathbb{R}^n$, and a Lipschitz continuous function $u : [0, \delta] \to K$ such that $x(0) = x^0$, $u^0 \in \text{SOL}(K, F(0, x^0, \cdot))$, and $(x(t), u(t))$ satisfies (2.2) for all $t \in [0, \delta]$. The existence of such solutions is treated in Sects. 5.3 and 5.4.

## 2.2 Classification of DVIs

There can be a substantial variation of the characteristics of DVIs regarding the existence and regularity of solutions. In order to provide a coherent framework for the study of DVIs we need to provide a rough classification of them. A significant guide in this process is the notion of *index* of DAEs [9,19]. In the latter theory, the index of a DAE is essentially the number of times that algebraic equations need to be differentiated in order to obtain an ODE (possibly after some manipulation of the equations). Obviously, a system with no algebraic constraints would be “index zero”. Thus if the algebraic equations uniquely define the algebraic variables as functions of the differential variables, one differentiation will give an ODE for the algebraic variables. Such a system is called “index one”. If the system requires more than one differentiation to obtain an ODE for the algebraic variables, then we have an “index two” or higher index DAE. (One could take the point of view that our index-one systems should be classified as index-zero systems because the algebraic variable can be substituted out, thereby obtaining an ODE in the differentiable variable only without differentiating the algebraic variable. Nevertheless, throughout the paper, we count the index as just explained.)
Unlike the DAE, a straightforward differentiation of the VI condition in a DVI is not possible. In spite of this handicap, a similar classification of the DVI can be used. Consider, for example, the IVP (2.2) with $x(0) = x^0$, where $F(t, x, \cdot)$ is a strongly monotone function. Under a certain uniformity of the strong monotonicity and a Lipschitz continuity of $F$, the map $(t, x) \mapsto \text{SOL}(K, F(t, x, \cdot))$ is a single-valued Lipschitz function; and so the (algebraic) variable $u$ can be written as a Lipschitz function of $(t, x)$. Substituting this into the differential equation $\dot{x} = f(t, x, u)$ gives a Lipschitz ODE to which standard theory can be applied. This is the *index one* case. For details of this case, see Sect. 5.1.

A variation on the above case is where $\text{SOL}(K, F(t, x, \cdot))$ might have multiple solutions. In that case, substitution of solutions $u$ into the differential equation would give not an ODE, but a DI:

$$ \dot{x} \in \mathbf{F}(t, x) \equiv f(t, x, \text{SOL}(K, F(t, x, \cdot))). \qquad (2.5) $$

If the set-valued map $\mathbf{F}$ has certain desirable convexity and linear growth properties (see Sect. 6), we call this the *multivalued index one* case.

Higher index DVIs are important too. In particular, we may consider *index two* problems where $u$ does not appear in the VI: $u(t) \in \text{SOL}(K, F(t, x))$, but $J_x F(t, x) f(t, x, u)$ is a monotone function of $u$. (Here, $J_x F(t, x)$ denotes the partial Jacobian matrix of $F(t, x)$ with respect to the variable $x$.) In this case, $u$ may be discontinuous, and little more than weak convergence can be expected of numerical or other approximations. For a proper treatment of this case, we need $f(t, x, u)$ to be an affine function of $u$; i.e., $f(t, x, u) \equiv \tilde{f}(t, x) + B(t, x)u$. There are also *mixed index one and index two* problems where $u$ can be split into two (vector) components, and the VI is independent of one component but strongly monotone in the other. The latter case is treated more generally in Sect. 8.

This paper treats exclusively DVIs with index no higher than two and which have absolutely continuous solutions. As we will see, the detailed analysis of these systems already present considerable technical challenges. Even higher index DVIs can be defined. Nevertheless, their theory and computation are generally complicated by a number of issues. Numerical schemes can fail to converge, and the theory of existence of solutions relies on high-order distributions such as derivatives of Dirac-$\delta$ functions, which are typically inappropriate in nonlinear settings. Even in the case of affine time-invariant functions $f$ and $F$ (which includes the Linear Complementarity Systems (LCSs) treated in [52]), there are questions of the interpretation of the existence theory for index four or higher DVIs: even with consistent initial conditions and analytic data, the solution has to contain derivatives of Dirac-$\delta$ functions for which statements such as “$u(t) \ge 0$” are problematic. Even for index three DVIs, uniqueness is known to fail, and the set of solutions is unbounded (even in the space of distributions). In spite of these difficulties in general, there are occasional situations where solutions to higher-order systems can be obtained by reduction to lower-order systems via suitable algebraic operations and functional differentiation. Due to the length of this paper, the details of such reduced systems are omitted. For a
recent distributional treatment of higher-order measure differential inclusions in the framework of Moreau's sweeping process, see [1].

## 2.3 Connections to existing problem classes

Before discussing the special cases of (2.2), we show how the DVI is conceptually equivalent to a DAE with a nonsmooth algebraic equation. The cornerstone of this equivalence is the casting of the variational condition as a nonsmooth equation [38]. Based on the *natural map* of the VI, we have $u \in \text{SOL}(K, F(t, x, \cdot))$ if and only if $0 = u - \Pi_K(u - F(t, x, u))$, where $\Pi_K$ denotes the Euclidean projector onto the closed convex set $K$. The nondifferentiability of the projector $\Pi_K$ is well known (see [38, Chap. 4] for a detailed study on the differentiability properties of this operator and for references); therefore, the crux of the existing theory of DAEs is not applicable to the DVI, except under restrictive assumptions. Nevertheless, the notion of the *index* of a (smooth) DAE and the associated differentiation process that is involved in the determination of such an index [19] have some relevance when dealing with *differential complementarity problems* (DCPs); see Sects. 5.3 and 5.4. A potential benefit of the (nonsmooth) DAE point of view of the DVI is that one could presumably borrow from the theory of nonsmooth (algebraic) equations [38] for the study of the DVI as the former theory becomes more developed. However, a detailed investigation of this venue is beyond the scope of this paper.

We have already noted the reformulation of the DVI (2.2) as the DI (2.5). This simple conversion allows us to employ the powerful machinery of DIs for the study of the DVI. Yet a word of caution is needed. Namely, the existing DI theory as documented in such references as [10,32,66,89] is cast in a very abstract framework. Substantial work is needed to verify the assumptions made in the DI theory; our task is therefore to derive conditions on the VI that ensure the satisfaction of such DI assumptions. In addition, we are also interested in considering cases of the DVI where the DI results are not directly applicable, and in numerical methods for boundary-value DVIs, which have received to date only minimal attention in the DI literature.

### *Special cases*

We have mentioned that the DVI is conceptually equivalent to a nonsmooth DAE. With $K \equiv \mathbb{R}^m$, (2.2) becomes

$$ 
\begin{aligned}
\dot{x}(t) &= f(t, x(t), u(t)), \\
0 &= F(t, x(t), u(t)).
\end{aligned}
\qquad (2.6) $$

Therefore, we recover a standard DAE as a special case of the DVI. Another important special case of the DVI is when $K$ is a cone $C$. In this case, we obtain
the differential complementarity problem (DCP):

$$ \begin{gathered} \dot{x}(t) = f(t, x(t), u(t)), \\ C \ni u(t) \perp F(t, x(t), u(t)) \in C^*. \end{gathered} \tag{2.7} $$

In turn, a proper specialization of the latter DCP yields the linear complementarity system (LCS) studied extensively in [21,49–52,87]. Specifically, the LCS is

$$ \begin{gathered} \dot{x} = p + Ax + Bu, \\ 0 \le u \perp q + Cx + Du \ge 0, \end{gathered} \tag{2.8} $$

which is clearly a special case of (2.7) where

$$ f(t, x, u) \equiv p + Ax + Bu \quad \text{and} \quad F(t, x, u) \equiv q + Cx + Du, $$

and $C \equiv \mathbb{R}_+^m$. Here $p$ and $q$ are given $n$-vector and $m$-vector, respectively, and $A, B, C$, and $D$ are given matrices of order $n \times n, n \times m, m \times n$, and $m \times m$, respectively. As in the study of the finite-dimensional LCP [30], the symmetry of $D$ is often not assumed. The cited references contain, among other things, some basic results on the “well-posedness” of the LCS and a detailed discussion of various solution concepts for such a complementarity system and how they relate to each other and to some related solution concepts for hybrid systems. Based on different approaches, the existence results that we obtain later for the LCS will complement those in the cited references. Most recently, some systematic-theoretic issues of the DCP and the LCS have been treated in [22,23,72,88].

Pioneered by Chen and Mangasarian [24,25], a popular approach for solving finite-dimensional complementarity problems is via smoothing of the complementarity condition. Applied to the LCS (2.8), this yields a sequence of smoothed DAEs whose convergence is worthy of investigation. It is also possible to extend the smoothing approach to the DCP (2.7) for a polyhedral cone $C$. However, these topics are beyond the scope of this paper.

## 2.4 Variational inequalities of evolution

In Definition 1 in [10, Sect. 5.6], Aubin and Cellina defined the “differential variational inequality” to mean a problem of the following type. Let $K$ be a closed convex subset of $\mathbb{R}^n$ and $F$ be a given (single-valued) mapping from $\mathbb{R}^n$ into itself. (Consistent with our development throughout this paper, we take $F$ to be single-valued, although the reference allows $F$ to be set-valued.) Find an absolutely continuous function $x(t)$ such that (i) $x(t) \in K$ for all $t \in [0, T]$, and (ii) $\dot{x}(t) \in -F(x(t)) - N(x(t); K)$ for almost all $t \in [0, T]$, where, as in convex analysis [85], $N(z; K)$ denotes the normal cone of $K$ at an element $z \in K$; i.e.,

$$ N(z; K) \equiv \{ v \in \mathbb{R}^n : v^T(y - z) \le 0, \forall y \in K \}. $$
In order to distinguish this problem from our DVI, we call the above differential problem a *variational inequality of evolution* and denote it as VIE $(K, F)$, which is consistent with the notation of the VI $(K, F)$. (In the recent paper [28], the authors chose the term “evolutionary variational inequality” to mean a VI in the Hilbert space $L^p([0, T], \mathbb{R}^q)$ that involves no time derivative. This terminology clashes with the time-varying system $\dot{x}(t) \in -F(x(t)) - \mathcal{N}(x(t); K)$, that has been called an “evolution variational inequality” in previous literature; see e.g. [44,46]. Our terminology VIE avoids such conflicts.) The VIE is a subclass of the class of DIs with maximal monotone right-hand sides which has been studied extensively as early as in 1970s; see e.g. Brézis [16]. One motivation for considering the VIE is that its stationary solutions, i.e., those for which $\dot{x} = 0$, correspond to solutions of the static VI defined by the pair $(K, F)$. Thus the VIE represents the latter VI in disequilibrium/evolution.

The VIE was apparently introduced by Henry [53,54] as a class of differential inclusions known as *projected differential inclusions*; this author was interested in solving certain constrained autonomous ODEs arising from some financial models in mathematical economics. An early study of the “complementarity problem of evolution” was documented in the 1993 Ph.D. thesis of Hipfel [55], who used the terminology DVI for the VIE. In the same year, Dupuis and Nagurney [36] published their study on the class of “projected dynamical systems”. For a general existence result for such a differential system in Hilbert space, see [27,29]. The recent papers [2,44–46] have developed an extensive stability theory for the VIE.

As an alternative dynamic disequilibrium model, we could consider the ODE: $\dot{x} = x - \Pi_K(x - F(x))$ whose “equilibria” (i.e., vectors $x^e$ such that $0 = x^e - \Pi_K(x^e - F(x^e))$) are precisely the solutions of the finite-dimensional VI $(K, F)$. The solution trajectories of the latter ODE, which are necessarily $C^1$ functions of time provided that the function $F$ is Lipschitz continuous, can be thought of as smooth pathways to such equilibria. A major computational advantage of this kind of pathways compared to the trajectories derived from the above VIE, which is in the form of a DI, is obvious; namely, we deal with Lipschitz ODEs instead of DIs; for some stability results of the ODE in question, see [22].

When $K$ is a cone, the VIE $(K, F)$ is a special case of our DVI. To see this, note that a vector $u$ belongs to $-\mathcal{N}(z; K)$ if and only if $K^* \ni u \perp z \in K$. Thus the DI $\dot{x} \in -F(x(t)) - \mathcal{N}(x(t); K)$ is equivalent to

$$ 
\begin{gathered}
\dot{x}(t) = -F(x(t)) + u(t), \\
K^* \in u(t) \perp x(t) \in K,
\end{gathered}
\qquad (2.9)
$$

which is in the form of the DCP (2.7) with $f(t, x, u) \equiv -F(x) + u$, $F(t, x, u) \equiv x$, and $C \equiv K^*$. In the case where $K$ is finitely representable by convex inequalities, which is fairly common in many applications, we can rely on the explicit representation of the normal cone using multipliers to express the VIE also as a DCP. Indeed, for simplicity, suppose that

$$ K \equiv \{x \in \mathbb{R}^n : g(x) \le 0\}, $$
where $g : \mathbb{R}^n \to \mathbb{R}^m$ is continuously differentiable such that each component function $g_i$ is convex, for $i = 1, \dots, m$. Further assume that $K$ obeys some standard constraint qualifications from nonlinear programming (e.g. $K$ is a polyhedron or $K$ contains an interior point) so that for every $x \in K$, a vector $v$ belongs to $\mathcal{N}(x; K)$ if and only if there exist multipliers $u \in \mathbb{R}^m$ satisfying

$$v = \sum_{i=1}^{m} u_i \nabla g_i(x) \quad \text{and} \quad 0 \le u \perp g(x) \le 0.$$

Consequently, the DI $\dot{x} \in -F(x(t)) - \mathcal{N}(x(t); K)$ is equivalent to

$$\begin{aligned} \dot{x}(t) = -F(x(t)) - \sum_{i=1}^{m} u_i \nabla g_i(x), & \qquad (2.10) \\ 0 \le u \perp -g(x) \ge 0, & \end{aligned}$$

which can be identified as the DCP (2.7) with $f(t,x,u) \equiv -F(x) - Jg(x)^Tu$, $C \equiv \mathbb{R}_+^m$, and $F(t,x,u) \equiv -g(x)$, where $Jg(x)$ denotes the Jacobian matrix of $g$ at $x$. The derivation of the above Karush-Kuhn-Tucker (KKT) form of the VIE is a standard exercise in nonlinear programming; in the context of the linear evolution VI, the formulation (2.10) has been noted previously in [44] for a polyhedral $K$.

For a general closed convex set $K$ with no particular structure, the VIE and the DVI still bears a close connection. Specifically, by introducing an auxiliary variable, the DI $\dot{x} \in -F(x(t)) - \mathcal{N}(x(t); \tilde{K})$ is equivalent to

$$\begin{aligned} \dot{x}(t) &= -F(x(t)) + w(t), \\ 0 &= x(t) - y(t), \\ y(t) \in \tilde{K}, \quad (\nu - y(t))^T w(t) \ge 0, \quad \forall \nu \in \tilde{K}. \end{aligned} \qquad (2.11)$$

This is in the form of the DVI with $u \equiv (w, y)$, $f(t,x,u) \equiv -F(x) + w$, $F(t,x,u) \equiv (x - y, w)$ and $K \equiv \mathbb{R}^n \times \tilde{K}$. It is interesting to note that when $\tilde{K}$ is a cone, the auxiliary variable $y$ is not needed. In general, the DVI is a broader mathematical model than the VIE as defined above. Nevertheless, by extending (2.11) to the following form:

$$\begin{aligned} \dot{x}(t) &= f(x(t), y(t), w(t)), \\ 0 &= g(x(t), y(t), w(t)), \\ y(t) \in \tilde{K}, \quad (\nu - y(t))^T w(t) \ge 0, \quad \forall \nu \in \tilde{K}, \end{aligned} \qquad (2.12)$$

we obtain a unified formulation for both the DVI and the VIE.
## 2.5 Extended systems

There are many generalizations of the VI (see [38, Sect. 1.6]), each providing an extension for the DVI. One such generalization leads to the differential quasi-variational inequality (DQVI) in which the fixed set $K$ is replaced by a time-dependent set-valued map from $[0, T] \times \mathbb{R}^m$ into closed convex subsets of $\mathbb{R}^m$: i.e., a system of the kind:

$$ 
\begin{aligned}
\dot{x}(t) &= f(t, x(t), u(t)), \\
u(t) &\in \text{SOL}(K(t, x(t)), F(t, x(t), \cdot)).
\end{aligned}
 $$

Properly specializing the functions $f$ and $F$, we obtain the sweeping processes developed by Moreau, Monteiro-Marques, Kunze and others [20, 61, 62, 66–69, 96]. These processes can be formulated as DIs of the kind:

$$ \dot{x}(t) \in -F(x(t)) - \mathcal{N}(x(t); K(t, x(t))), \quad (2.13) $$

where the moving set $K(t, x(t))$ is typically closed and convex. Sweeping processes have found application in the area of mechanical systems with impact, specifically for rigid-body dynamics with unilateral constraints and are particularly useful as they allow for *discontinuous* solutions, particularly if $K(\cdot, \cdot)$ is a discontinuous function in the Hausdorff metric; this is important for handling impact problems as the velocities are discontinuous in most impact problems. Care must be taken in interpreting (2.13) in this case: the derivative must be understood in the distributional sense. Typically $x(\cdot)$ is required to have bounded variation and the derivative $dx/dt$ is understood to be a measure. The meaning of the inclusion (2.13) is as a *measure differential inclusion* [66–68, 96]. Generalizations of (2.13) to nonconvex $K(t, x)$ can be found, for instance, in [15, 39, 100]. The DQVI also arises from the class of differential generalized Nash games; see Sect. 4.2.

Another extension of the DVI is the class of convolution complementarity problems (CCPs) introduced in [76] and further analyzed by the second author [97]. The CCP is useful for the study of linear time-invariant systems, such as delays and systems obtained from partial differential equations, and for the study of impacts of elastic bodies, and is interesting in its own right. Since this paper concerns itself mainly with DVIs with absolutely continuous solutions, and since all the mentioned extended dynamical systems generally involve solutions that belong to broader function classes, the detailed study of the extended systems is beyond the scope of this paper.

## 3 Computational Methods for DVIs

Since we are particularly interested in computational methods for DVIs, we find it useful to begin early in our discussion to address this important aspect of our study. In designing computational methods for solving a general DVI, one
has to take into account two important features of this mathematical model. These are: (i) time dependence, which refers to the ODE, and (ii) nonlinearity, which refers to the input functions $f, F$, and possibly $\Gamma$. The first issue is dealt with by *time stepping*, i.e., replacing the time derivative $\dot{x}$ by a finite-difference quotient; and the second is dealt with by some kind of *linearization*. Although an explicit time-stepping scheme applied to an IVP ODE can easily handle any nonlinearity in the right-hand function $f$ (perhaps at the expense of being less suitable than an implicit scheme for some problems), the nonlinearity of $F$ (and $\Gamma$) remains to be dealt with by an algebraic method. For a BVP, the distinction between implicit and explicit time stepping disappears for the obvious reason that $x^0$ is not completely known; hence, any nonlinearity of $f$ cannot be removed in time stepping and has to be dealt with by an algebraic method in any practical implementation.

In light of the above considerations, one can therefore develop two general families of iterative methods for solving a DVI, depending on whether time stepping or linearization is the “outer” iteration and the other is the “inner” iteration. This paper focuses on the family of methods where time stepping is the outer scheme. A typical method of this kind begins with the division of the time interval $[0, T]$ into $N_h + 1$ subintervals

$$0 = t_{h,0} < t_{h,1} < t_{h,2} < \cdots < t_{h,N_h} < t_{h,N_h+1} = T,$$

each of length $h > 0$. Thus $(N_h + 1)h = T$ and

$$t_h = t_{h,i} + h, \quad \forall i = 0, 1, \dots, N_h.$$

We approximate the time derivative $\dot{x}$ by a backward finite-difference quotient:

$$\dot{x}(t_{h,i+1}) \approx \frac{x(t_{h,i+1}) - x(t_{h,i})}{h}, \quad i = 0, 1, \dots, N_h.$$

We employ a scalar $\theta \in [0,1]$ to distinguish an explicit ($\theta = 1$), an implicit ($\theta = 0$), or a semi-implicit ($\theta \in (0,1)$) discretization of the ODE. We then define a discrete-time system to compute two finite family of vectors

$$\{x^{h,1}, x^{h,2}, \dots, x^{h,N_h+1}\} \subset \mathbb{R}^n \quad \text{and} \quad \{u^{h,1}, u^{h,2}, \dots, u^{h,N_h+1}\} \subset \mathbb{R}^m. \quad (3.1)$$

The individual methods differ in how the latter system is defined.

## Sub-scheme I: The DI formulation

In this sub-scheme, we rely on the reformulation of the DVI as the DI (2.5) and consider the following discrete-time system:

$$x^{h,i+1} \in x^{h,i} + h \mathbf{F}(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1}), \quad i = 0, 1, \dots, N_h.$$
Unwrapping the set-valued map **F**, we obtain

$$
\left.
\begin{aligned}
x^{h,i+1} &= x^{h,i} + h f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1}, u^{h,i+1}) \\
u^{h,i+1} &\in \text{SOL}(K, F(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1}, \cdot)) \\
0 &= \Gamma(x^{h,0}, x^{h,N_h+1})
\end{aligned}
\right\} , \quad (3.2)
$$

Dontchev and Lempio [35] have surveyed various difference methods for initial-value DIs, including the explicit Euler method, linear multistep methods, implicit Runge-Kutta methods, among others. In principle, one could apply their results to the DVI as long as the required conditions for convergence are satisfied; among other things, we would need the solvability of the discrete-time system (3.2). The verification of these conditions is not a trivial task and some of them could be fairly restrictive in the context of the DVI. Since we are interested in identifying broad classes of DVIs for which convergent time-stepping methods can be developed, we will not pursue the DI-based approach in this paper; instead, we will present next a DVI-specific scheme, in which we respect the VI constraint exactly at each time step; see this constraint in (3.3).

Sub-scheme II: A DVI-specific approach

A distinguished feature of the iteration (3.2) is the appearance of the same
parameter $\theta$ in both the discretized differential equation and the algebraic VI.
This of course is due to the absorption of the variational condition into the
set-valued map **F**. In contrast, by not using the DI formulation, one has the
freedom of using different parameters in the first two equations in (3.2). For
instance, we may consider the following variant:

$$
\left.
\begin{aligned}
x^{h,i+1} &= x^{h,i} + h f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1}, u^{h,i+1}) \\
u^{h,i+1} &\in \text{SOL}(K, F(t_{h,i+1}, x^{h,i+1}, \cdot)) \\
0 &= \Gamma(x^{h,0}, x^{h,N_h+1})
\end{aligned}
\right\} , \quad (3.3)
$$

in which the VI is satisfied exactly by each iterate $(x^{h,i+1}, u^{h,i+1})$ at time $t_{h,i+1}$.
Another variant is

$$
\left.
\begin{aligned}
x^{h,i+1} &= x^{h,i} + h f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1}, u^{h,i+1}) \\
u^{h,i+1} &\in \text{SOL}(K, F(t_{h,i+1}, x^{h,i}, \cdot)) \\
0 &= \Gamma(x^{h,0}, x^{h,N_h+1})
\end{aligned}
\right\} , \quad (3.4)
$$

where the VI can be solved independently of the first equation (assuming $x^{h,i}$ is
known). These two schemes illustrate that by separating the algebraic from the
differential condition, one obtains more flexibility in designing finite difference
methods for solving the DVI than by including the VI within the abstract
function F. See Sects. 9.1 and 10 for situations where such flexibility is put to good use.

Besides the flexibility, another significant difference between a DI-based time-stepping method versus a DVI-specific method is that the latter approach exposes the role of the algebraic variable *u* while the former approach buries it. As a consequence, the behavior of the sequence {*u*<sup>*h*,*i*</sup>} is not of concern in the DI-based approach; whereas in the DVI-specific approach, both the computation and convergence of the algebraic *u*-iterates is as important as that of the differential *x*-iterates. From a modeling and applications point-of-view, the information in the algebraic variable *u* can be of crucial importance and therefore should not be ignored.

For an IVP, (3.3) can be solved recursively starting with *i* = 0. Indeed with *x*<sup>*h*,*i*</sup> given, the first two conditions in (3.3) define a VI (K̂, *H*<sup>*h*,*i*+1</sup>) in the variable (*x*<sup>*h*,*i*+1</sup>, *u*<sup>*h*,*i*+1</sup>), where K̂ ≡ ℝ<sup>*n*</sup> × K and

$$
H^{h,i+1}(x,u) \equiv \left( \begin{array}{c} x - x^{h,i} - h f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x, u) \\ F(t_{h,i+1}, x, u) \end{array} \right), \quad (x,u) \in \hat{K}.
$$

For the BVP, such decoupling of the subproblems at the individual time steps
is no longer valid, and we have to deal with a very large finite-dimensional VI
in the variables (3.1), whose size becomes more pronounced as the time step *h*
tends to zero.

There are three major issues associated with the above time-stepping scheme that need to be addressed: (i) existence of a solution to the discrete-time system (3.3), (ii) numerical resolution of the latter system, and (iii) convergence of the discrete-time trajectory to a continuous-time solution as *h* ↓ 0. In later sections, we shall treat the first and third issue in detail, but only give a pedestrian discussion with reference to the second computational issue. Needless to say, one could consider advanced time-stepping methods for solving the DVI; nevertheless, such an investigation is beyond the scope of this paper.

A word of caution should be said about the schemes (3.3) and (3.4). Na-
mely, their convergence relies on the the data functions *f*(*t*, *x*, *u*) and *F*(*t*, *x*, *u*).
In particular, the convergence theory developed in this paper assumes that
*f*(*t*, *x*, *u*) ≡ *f̃*(*t*, *x*) + *B*(*t*, *x*) *u*, with *f̃* and *B* Lipschitz continuous jointly in time
and state (see Sect. 6.1 and Assumption (A) therein). From the theory of diffe-
rential inclusions, such continuity seems restrictive because for time-dependent
problems, one typically assumes only continuity of the Carathéodory kind;
i.e., the functions are Lipschtiz continuous in the state for fixed times with
time-dependent Lipschitz constants that are integrable functions of the state.
For problems with *f* having such loose properties, a more appropriate scheme
could be

$$
\left.
\begin{aligned}
x^{h,i+1} &= x^{h,i} + \int_{t_{h,i}}^{t_{h,i+1}} \tilde{f}(\tau, \theta x^{h,i} + (1-\theta)x^{h,i+1}) \, d\tau \\
&\quad + h B(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1}) u^{h,i+1} \\
u^{h,i+1} &\in \text{SOL}(K, F(t_{h,i+1}, x^{h,i+1}, \cdot)) \\
0 &= \Gamma(x^{h,0}, x^{h,N_h+1})
\end{aligned}
\right\} . \quad (3.5)
$$

Our focus in this paper is on the properties of the DVI that are the result of incorporating variational inequalities into dynamics. Thus, we focus on discontinuities due to the interaction of dynamics and variational inequalities rather than discontinuities due to the external inputs and the model functions. Discontinuous and non-Lipschitz inputs can be handled, albeit at the cost of additional technical complication, such as having to deal with the integral scheme (3.5) instead of the finite-difference schemes (3.3) and (3.4). We wish to avoid overburdening the reader with this additional complication in an already very long paper, and so will restrict the data functions in the DVI to be at least Lipschitz in all variables. Nonetheless, we will indicate how our results can be extended to include functions with intermediate Carathéodory-type conditions; see Assumption (A') in Sect. 7.

4 Source problems

In the following two subsections, we discuss in detail two source problems of the
DVI: an ODE with a discontinuous right-hand-side and the differential Nash
game. For applications of the LCS to electric circuits with ideal diodes and to
piecewise linear dynamical systems, see [21, 49, 51, 52, 87]. The survey paper [31]
discusses a genetic regulatory system modeled as an ODE with a piecewise
linear right-hand side that can be cast as an LCS. There are many other appli-
cations, such as the dynamic traffic assignment problem [11, 64, 78], multibody
frictional contact problems (see [73, 94] and the references in the bibliography),
which typically lead to DVIs of index two and higher (therefore falling outside
the scope of this study), differential systems subject to optimization constraints,
which have applications in chemical equilibrium problems such as in the mode-
ling of atmospheric aerosols [59, 75], as well as many applications of the VIE,
which has a large body of literature by itself.

4.1 ODEs with discontinuous right-hand sides

While pioneering work on ODEs with discontinuous right-hand sides was ini-
tiated by Filippov [41,43], the DCP approach to these ODEs originated from
the second author’s Ph.D. thesis [93], who used linear complementarity theory
to study “piecewise active solutions” to ODEs with piecewise smooth (but dis-
continuous) right-hand sides. In what follows, we use a simple vector-ODE to
illustrate how ODEs with discontinuous right-hand sides can be formulated as DCPs. Let $E \in \mathbb{R}^{n \times n}$ and consider:

$$ \dot{x} = E \operatorname{sgn}(x), \qquad (4.1) $$

where $\operatorname{sgn}(x)$ be the signum function with components $\operatorname{sgn}(x_i)$ with

$$ \operatorname{sgn}(\tau) = \begin{cases} = 1 & \text{if } \tau > 0, \\ \in [-1, 1] & \text{if } \tau = 0, \\ = -1 & \text{if } \tau < 0. \end{cases} $$

The inclusion $\operatorname{sgn}(0) \in [-1, 1]$ is the convexification proposed by Filippov to address the discontinuity of the signum function at the origin. As in elementary linear programming, we can write $x = x^{+} - x^{-}$, where $x^{\pm}$ are the nonnegative and nonpositive part of $x$, respectively. In terms of these two nonnegative (and complementary) variables $x^{\pm}$, we can further write $\operatorname{sgn}(x) = \mathbf{1} - v$, where $v$ satisfies

$$ 0 \leq v \perp x^{+} \geq 0 \quad \text{and} \quad 0 \leq \mathbf{2} - v \perp x^{-} \geq 0, $$

with $\mathbf{1}$ and $\mathbf{2}$ being the vectors of 1's and 2's, respectively. Since $x^{+} = x + x^{-}$, we obtain the following equivalent (inhomogeneous) LCS formulation of (4.1):

$$ \begin{aligned} \dot{x} &= E\mathbf{1} - Ev, \\ 0 &\leq v \perp x + x^{-} \geq 0, \\ 0 &\leq x^{-} \perp \mathbf{2} - v \geq 0, \end{aligned} $$

which corresponds to (2.8) with $u = (v, x^{-})$, and

$$ \begin{aligned} p &\equiv E\mathbf{1}, & A &\equiv 0, & B &\equiv [-E \quad 0], \\ q &\equiv \begin{pmatrix} 0 \\ \mathbf{2} \end{pmatrix}, & C &\equiv \begin{bmatrix} I \\ 0 \end{bmatrix}, & D &\equiv \begin{bmatrix} 0 & I \\ -I & 0 \end{bmatrix}. \end{aligned} \qquad (4.2) $$

Note that $D$ is a positive semidefinite matrix. Moreover, if $-E$ is positive semi-definite, then so is $D + \tau CB$ for all $\tau > 0$. These are important properties when it comes to solving the problem.

## 4.2 Differential Nash games

Consider a noncooperative Nash equilibrium problem with $N$ players where each player solves an optimal control problem. Let $(x^{\nu}, u^{\nu}) \in \mathbb{R}^{n_{\nu}+m_{\nu}}$ denote player $\nu$'s pair of state and control variables, respectively; let $U^{\nu} \subseteq \mathbb{R}^{m_{\nu}}$ be the
closed convex set of player $v$'s admissible controls; let $\mathbf{x} \equiv (x^v)_{v=1}^N \in \mathbb{R}^n$ and
$\mathbf{u} \equiv (u^v)_{v=1}^N \in \mathbb{R}^m$, where

$$
n \equiv \sum_{v=1}^{N} n_v \quad \text{and} \quad m \equiv \sum_{v=1}^{N} m_v
$$

be the collection of all players' variables. Dependent on all players' variables,
player v's cost function is given by

$$
\theta_v(\mathbf{x}, \mathbf{u}) \equiv \psi_v(\mathbf{x}(T)) + \int_0^T \varphi_v(t, \mathbf{x}(t), \mathbf{u}(t)) \, dt,
$$

where $T$ denotes the terminal time and $\psi_v : \mathbb{R}^n \to \mathbb{R}$ and $\varphi_v : \mathbb{R}^{1+n+m} \to \mathbb{R}$
are continuously differentiable functions. Parameterized by the rival players'
strategies $x^{-v} \equiv (x^{v'})_{v' \neq v}$ and $u^{-v} \equiv (u^{v'})_{v' \neq v}$, which are taken to be exogenous
to player $v$'s optimization problem but endogenous to the overall model, player
$v$'s optimization problem is to compute an optimal trajectory $(x^v(t), u^v(t))$, for
each fixed but arbitrary $(x^{-v}, u^{-v})$, in order to

$$
\begin{array}{ll}
\text{minimize} & \theta_v(\mathbf{x}, \mathbf{u}) \\
\text{subject to} & \dot{x}^v(t) = g^v(t, x^v(t), u^v(t)), \quad x^v(0) = x^{0,v} \quad \text{and} \\
& u^v(t) \in U^v \subseteq \mathbb{R}^{m_v}, \quad \text{for almost all } t \in [0, T],
\end{array}
$$

where $g^v : \mathbb{R}^{1+n_v+m_v} \to \mathbb{R}^{l_v}$ is continuously differentiable. A differential Nash
\textit{equilibrium solution} [12,34] is a pair of trajectories $(\mathbf{x}^*(t), \mathbf{u}^*(t))$ such that for
every $v = 1, \dots, N$, $(x^{*,v}(t), u^{*,v}(t))$ solves player $v$'s problem, given that $v$'s ri-
vals all play their Nash strategies $(x^{*,-v}(t), u^{*,-v}(t))$. To derive a set of necessary
conditions for a differential Nash equilibrium trajectory, let

$$
H_v(t, \mathbf{x}, \mathbf{u}, \lambda^v) \equiv \varphi_v(t, \mathbf{x}, \mathbf{u}) + (\lambda^v)^T g^v(t, \mathbf{x}^v, \mathbf{u}^v), \quad \lambda^v \in \mathbb{R}^{l_v},
$$

be player $v$'s Hamiltonian function, where $\lambda^v$ is the adjoint variable of the ODE constraint in player $v$'s control problem. The first-order necessary optimality condition of player $v$'s problem is

$$
\begin{align*}
\dot{\lambda}^v(t) &= -\nabla_{x^v} H_v(t, \mathbf{x}(t), \mathbf{u}(t), \lambda^v(t)), \\
\dot{x}^v(t) &= g^v(t, x^v(t), u^v(t)), \\
u^v(t) &\in \text{SOL}(U^v, \nabla_{u^v} H_v(t, \mathbf{x}(t), u^{-v}(t), \cdot, \lambda^v(t)), \\
x^v(0) &= x^{0,v} \quad \text{and} \quad \lambda^v(T) = \nabla_{x^v} \psi_v(\mathbf{x}(T)).
\end{align*}
$$
Concatenating these conditions for all players, we obtain the BV-DVI:

$$ 
\begin{gathered}
\begin{pmatrix} \dot{\lambda}(t) \\ \dot{x}(t) \end{pmatrix} = \begin{pmatrix} \mathbf{H}(t, \mathbf{x}, \mathbf{u}, \boldsymbol{\lambda}) \\ \mathbf{g}(t, \mathbf{x}, \mathbf{u}) \end{pmatrix}, \\
\mathbf{u}(t) \in \mathrm{SOL}(\mathbf{U}, \mathbf{F}(t, \mathbf{x}(t), \boldsymbol{\lambda}(t), \cdot)), \\
0 = \Gamma(\mathbf{x}(0), \mathbf{x}(T), \boldsymbol{\lambda}(T)),
\end{gathered}
$$

where $\mathbf{U} = \prod_{v=1}^{N} U^v \subseteq \Re^m$ is the Cartesian product of lower-dimensional sets,

$$ 
\begin{gathered}
\mathbf{H}(t, \mathbf{x}, \mathbf{u}, \boldsymbol{\lambda}) = -(\nabla_{x^\nu} H_\nu(t, \mathbf{x}, \mathbf{u}, \boldsymbol{\lambda}^\nu))_{\nu=1}^N, \quad \mathbf{g}(t, \mathbf{x}, \mathbf{u}) = (g^\nu(t, \mathbf{x}, \mathbf{u}))_{\nu=1}^N, \\
\mathbf{F}(t, \mathbf{x}, \boldsymbol{\lambda}, \mathbf{u}) = (\nabla_{u^\nu} H_\nu(t, \mathbf{x}, \mathbf{u}, \boldsymbol{\lambda}^\nu))_{\nu=1}^N, \quad \Gamma(\mathbf{y}, \mathbf{z}, \boldsymbol{\lambda}) = \begin{pmatrix} y^\nu - x^{0,\nu} \\ \boldsymbol{\lambda}^\nu - \nabla_{z^\nu} \psi_\nu(\mathbf{z}) \end{pmatrix}_{\nu=1}^N.
\end{gathered}
$$

In a generalized differential Nash game, each player's constraints can depend on his rivals' variables. As an illustration, let us consider the case where this occurs with the set of controls; thus the constant set $U^\nu$ is replaced by the variable set $U^\nu(x^{-\nu}, u^{-\nu})$. In this case, the above DVI becomes the DQVI where the variational condition: $\mathbf{u}(t) \in \mathrm{SOL}(\mathbf{U}, \mathbf{F}(t, \mathbf{x}(t), \boldsymbol{\lambda}(t), \cdot))$ is replaced by the quasi-variational condition: $\mathbf{u}(t) \in \mathrm{SOL}(\mathbf{U}(\mathbf{x}(t), \mathbf{u}(t)), \mathbf{F}(t, \mathbf{x}(t), \boldsymbol{\lambda}(t), \cdot))$, where $\mathbf{U}(\mathbf{x}(t), \mathbf{u}(t)) = \prod_{v=1}^{N} U^\nu(x^{-\nu}, u^{-\nu})$.

### The linear-quadratic case

The linear-quadratic differential Nash game has three specifications: (a) each $\psi_\nu(\mathbf{x}(T))$ is quadratic in $\mathbf{x}(T)$: $\psi_\nu(\mathbf{x}(T)) = x^\nu(T)^T [c^\nu + \sum_{v' \in \mathcal{N}} W^{\nu v'} x^{v'}(T)]$ for some constant vector $c^\nu \in \Re^{n_\nu}$ and matrices $W^{\nu v'} \in \Re^{n_\nu \times n_{v'}}$ with $W^{\nu v}$ symmetric (note: the terms not involving $x^\nu(T)$ are not included because they are exogenous to player $v$'s optimization problem); and (b) each $\varphi_\nu(t, \mathbf{x}, \mathbf{u})$ is quadratic in $(\mathbf{x}, \mathbf{u})$: for some time-dependent vectors $p^\nu(t) \in \Re^{n_\nu}$ and $q^\nu(t) \in \Re^{m_\nu}$ and time-dependent matrices $P^{\nu v'}(t) \in \Re^{n_\nu \times n_{v'}}$, $Q^{\nu v'}(t) \in \Re^{n_\nu \times m_{v'}}$, $R^{\nu v'}(t) \in \Re^{m_\nu \times n_{v'}}$, and $S^{\nu v'}(t) \in \Re^{m_\nu \times m_{v'}}$:

$$ 
\begin{aligned}
\varphi_\nu(t, \mathbf{x}, \mathbf{u}) = (x^\nu)^T \left[ p^\nu(t) + \sum_{v' \in \mathcal{N}} \left( P^{\nu v'}(t) x^{v'} + Q^{\nu v'}(t) u^{v'} \right) \right] \\
+ (u^\nu)^T \left[ q^\nu(t) + \sum_{v' \in \mathcal{N}} \left( R^{\nu v'}(t) x^{v'} + S^{\nu v'}(t) u^{v'} \right) \right]
\end{aligned}
$$

and (c) each $g^\nu(t, x^\nu, u^\nu)$ is an affine function in $(x^\nu, u^\nu)$: for some vector $r^\nu(t) \in \Re^{n_\nu}$ and matrices $G^\nu(t) \in \Re^{n_\nu \times n_\nu}$ and $H^\nu(t) \in \Re^{n_\nu \times m_\nu}$:

$$ g^\nu(t, x^\nu, u^\nu) = r^\nu(t) + G^\nu(t) x^\nu + H^\nu(t) u^\nu, $$
Under the above specifications, we have

$$ \left( \begin{array}{c} \mathbf{H}(t, \mathbf{x}, \mathbf{u}, \boldsymbol{\lambda}) \\ \mathbf{g}(t, \mathbf{x}, \mathbf{u}) \end{array} \right) = \mathbf{d}(t) + \mathbf{A}(t) \left( \begin{array}{c} \boldsymbol{\lambda} \\ \mathbf{x} \end{array} \right) + \mathbf{B}(t)\mathbf{u}, \quad (4.3) $$

where

$$ \begin{align*} \mathbf{d}(t) &= \left( \begin{array}{c} (-p^v(t))_{v \in \mathcal{N}} \\ (r^v(t))_{v \in \mathcal{N}} \end{array} \right) \in \mathfrak{R}^{2n} \\ \mathbf{A}(t) &= \left[ \begin{array}{cc} -\mathrm{diag}(G^v(t)^T)_{v \in \mathcal{N}} & [-P^{vv'}(t)]_{v\nu' \in \mathcal{N}} - \mathrm{diag}(P^{vv}(t))_{v \in \mathcal{N}} \\ 0 & \mathrm{diag}(G^v(t))_{v \in \mathcal{N}} \end{array} \right] \in \mathfrak{R}^{2n \times 2n} \\ \mathbf{B}(t) &= \left[ \begin{array}{c} -[Q^{vv'}(t)]_{v\nu' \in \mathcal{N}} - \mathrm{diag}(R^{vv}(t)^T)_{v \in \mathcal{N}} \\ \mathrm{diag}(H^v(t))_{v \in \mathcal{N}} \end{array} \right] \in \mathfrak{R}^{2n \times m}. \end{align*} $$

Moreover,

$$ \mathbf{F}(t, \mathbf{x}, \boldsymbol{\lambda}, \mathbf{u}) = \mathbf{q}(t) + \mathbf{C}(t) \left( \begin{array}{c} \boldsymbol{\lambda} \\ \mathbf{x} \end{array} \right) + \mathbf{S}(t)\mathbf{u}, \quad (4.4) $$

where $\mathbf{S}(t) = [S^{vv'}(t)]_{v\nu' \in \mathcal{N}} + \mathrm{diag}(S^{vv}(t))_{v \in \mathcal{N}}$ and

$$ \mathbf{C}(t) = \left[ \mathrm{diag}(H^v(t)^T)_{v \in \mathcal{N}} [R^{vv'}(t)]_{v\nu' \in \mathcal{N}} + \mathrm{diag}(Q^{vv}(t)^T)_{v \in \mathcal{N}} \right] \in \mathfrak{R}^{m \times 2n}. $$

Finally, we also have $\nabla_{x^v(T)} \psi_v(\mathbf{x}(T)) = c^v + W^{vv}x^v(T) + \sum_{v' \in \mathcal{N}} W^{vv'}x^{v'}(T)$.

Rewriting the VI $(\mathbf{U}, \mathbf{F}(t, \mathbf{x}, \boldsymbol{\lambda}, \cdot))$ as the triple $(\mathbf{U}, \mathbf{q}(t, \mathbf{x}, \boldsymbol{\lambda}), \mathbf{S}(t))$, where

$$ \mathbf{q}(t, \mathbf{x}, \boldsymbol{\lambda}) = \mathbf{q}(t) + \mathbf{C}(t) \left( \begin{array}{c} \boldsymbol{\lambda} \\ \mathbf{x} \end{array} \right), $$

we see that the DNG can be formulated as:

$$ \left. \begin{array}{l} \left(\begin{array}{c} \dot{\boldsymbol{\lambda}} \\ \dot{\mathbf{x}} \end{array}\right) = \mathbf{d}(t) + \mathbf{A}(t)\left(\begin{array}{c} \boldsymbol{\lambda} \\ \mathbf{x} \end{array}\right) + \mathbf{B}(t)\mathbf{u} \\ \mathbf{u}(t) \in \mathrm{SOL}(\mathbf{U}, \mathbf{q}(t, \mathbf{x}(t), \boldsymbol{\lambda}(t)), \mathbf{S}(t)) \\ \mathbf{x}(0) = \mathbf{x}^0, \quad \boldsymbol{\lambda}(T) = \mathbf{c} + \mathbf{W}\mathbf{x}(T) \end{array} \right\} , \quad (4.5) $$

where $\mathbf{c} = (c^v)_{v \in \mathcal{N}}$ and $\mathbf{W} = [W^{vv'}]_{v,v' \in \mathcal{N}} + \mathrm{diag}(W^{vv})_{v \in \mathcal{N}}$. Notice that in (4.5), the right-hand side in the ODE is an affine function of the state variables $(\mathbf{x}, \boldsymbol{\lambda})$ and algebraic variable $\mathbf{u}$ and that the middle VI is defined by the affine mapping $\mathbf{u} \mapsto \mathbf{q}(t, \mathbf{x}(t), \boldsymbol{\lambda}(t)) + \mathbf{S}(t)\mathbf{u}$.
# 5 Reduction to Lipschitz ODEs

Beginning in this section, we undertake a systematic study of the DVI (2.2) and (2.3). We start with the case where (2.2) can be reduced to an ODE with an implicitly defined Lipschitz continuous right-hand side [26]. This situation is reminiscent of an index-one DAE (2.6), whereby under a nonsingular assumption on the partial Jacobian matrix $J_u F(0, x^0, u^0)$, where $u^0$ satisfies $F(0, x^0, u^0) = 0$, we can obtain a Lipschitz continuous implicit function $u(t, x)$, defined near $(0, x^0)$, that satisfies the algebraic equation $F(t, x, u(t, x)) = 0$ for all $(t, x)$ in the domain of the implicit function $u(t, x)$. Substituting the latter function $u(t, x)$ into the ODE $\dot{x} = f(t, x, u)$ converts the DAE into an ODE with a Lipschitz right-hand side, provided that $f$ is Lipschitz continuous.

By employing known results from VI theory, the above reduction for a DAE can be extended to the DVI. We distinguish between a global reduction and a local conversion; in both cases, we emphasize that the right-hand function of the resulting ODE is only implicitly defined, whose evaluation requires the solution of a finite-dimensional VI. While this is an important consideration in practical computations, the reduction enables the well-established theory and methods of ODE to be directly applicable to the DVI, thereby no new theoretical treatment of the DVI is needed in this case. Another noteworthy remark to emphasize is that there is a major gap between the global and local Lipschitz ODE reformulations. While a global reformulation requires restrictive assumptions, a local reformulation, if possible, can be useful for analyzing local properties near a system state. An example of such an analysis can be found in [72] where the issue of “non-Zenoness” of a strongly regular DVI is studied (see Sect. 5.3 for the definition of strong regularity).

## 5.1 Global reduction

Similar to an index-one DAE, the key idea behind the global reduction for a DVI is that for certain classes of $(K, F)$, the VI $(K, F(t, x, \cdot))$ has a unique solution, say $u(t, x)$, that is a Lipschitz continuous function of $(t, x)$ for all $(t, x) \in [0, T] \times \dot{X}$, where $X$ is a certain domain wherein $F$ has the required properties for such a solution $u(t, x)$ to exist. In what follows, we spell out the conditions on the pair $(K, F)$ that makes the global reduction possible.

To cover a wide class of problems, we take $K$ to be the Cartesian product of finitely many lower-dimensional sets:

$$ K = \prod_{v=1}^{N} K^{v}, \qquad (5.1) $$

where each $K^v$ is a closed convex subset of $\mathbb{R}^{n_v}$. We postulate the following conditions on the pair $(K, F)$:
(A₀) $F(t, x, \cdot)$ is a continuous, uniformly $P$ function [38] on $K$ with a modulus that is independent of $(t, x)$; i.e., there exists a constant $\eta_F > 0$ such that

$$ \max_{1 \le \nu \le N} (u_\nu - u'_\nu)^T (F_\nu(t, x, u) - F_\nu(t, x, u')) \ge \eta_F \|u - u'\|^2 $$

for all $(t, x) \in [0, T] \times X$ and $u = (u_\nu)_{\nu=1}^N$ and $u' = (u'_\nu)_{\nu=1}^N$ in $K = \prod_{\nu=1}^N K^\nu$; (B₀) $F(\cdot, \cdot, u)$ is Lipschitz continuous with a constant that is independent of $u$; i.e., there exists a constant $L_F > 0$ such that

$$ \| F(t, x, u) - F(t', x', u) \| \le L_F [ |t - t'| + \|x - x' \| ], $$

for all $(t, x)$ and $(t', x')$ in $[0, T] \times X$ and all $u \in K$.

From well-known theory of partitioned VIs [38, Sect. 3.5], we have

**Theorem 5.1** Let $K$ be given by (5.1), where each $K^\nu$ is a closed convex subset of $\mathbb{R}^{n_\nu}$. Under the two assumptions (A₀) and (B₀), there exists a Lipschitz continuous function $u: [0, T] \times X \to K$ such that for each pair $(t, x) \in [0, T] \times X$, $u(t, x)$ is the unique solution of the VI $(K, F(t, x, \cdot))$.

In terms of the implicitly defined function $u(t, x)$, (2.2) is equivalent to the ODE:

$$ \dot{x} = f(t, x, u(t, x)), \qquad (5.2) $$

Hence, provided that $f$ is Lipschitz continuous, we arrive at the familiar domain of an ODE with a Lipschitz right-hand side. A word of caution is needed in the case of solving boundary-value DVIs via this reduction. A class of well-known methods for solving two-point boundary ODEs is the class of shooting methods [7, 60, 79], which invariably utilize some Newton methods for solving smooth systems of algebraic equations. In a subsequent paper [74], we have shown that the implicit function $u(t, x)$ is “B(ouligand)-differentiable” [38, Chap. 5] and is F(réchet)-differentiable only under highly restrictive assumptions. Hence, caution is needed in applying the shooting methods to (5.2). In principle, we could consider applying “nonsmooth Newton methods” instead (see [38, Chaps. 7 and 8]); this is indeed possible for problems of the affine type; for details, see [74].

## 5.2 Differential semi-affine systems

It is possible to relax assumption (A₀) when the set $K$ is polyhedral and the function $F(t, x, u)$ is linear in $u$ with a constant $J_u F(t, x, u)$. In this case, under the coherent orientation property [84, 38] for finite-dimensional affine variational inequalities, the same global reduction of a DVI can be accomplished. Instead of presenting the details of this case, we consider the following autonomous “semi-affine” IVP:

$$ \begin{array}{l} \dot{x} = f(x) + Bu, \quad x(0) = x^0 \\ u \in \text{SOL}(K, Cx, D), \end{array} \qquad (5.3) $$
where $f : \mathbb{R}^n \to \mathbb{R}^n$ is a Lipschitz function, $K$ is a polyhedron in $\mathbb{R}^m$, and $B, C$, and $D$ are constant matrices of appropriate orders. If the pair $(K, D)$ is coherently oriented in the sense of Robinson (see the cited references), then the above differential system (5.3) is globally equivalent to an ODE with a Lipschitz (albeit implicitly defined) right-hand side: $\dot{x} = f(x) + Bu(x)$, where $u(x)$ is the unique solution of the AVI $(K, Cx, D)$. In what follows, we state and prove a result, Proposition 5.1, that provides a necessary and sufficient condition for (5.3) to have a unique $C^1$ $x$-trajectory for all initial conditions $x^0$, without assuming the coherent orientation of the pair $(K, D)$. This result does not require nor imply that SOL$(K, Cx, D)$ is a singleton. The proof of the result makes use of a fundamental property (Lemma 5.1) of single-valued polyhedral multifunctions that first appeared in an unpublished report of M.S. Gowda; for a sketch of the proof of the latter property, see [38, Exercises 5.6.14] and also [14, Lemma 3.3]. The domain of a multifunction (i.e., set-valued function) $\Phi$, denoted $\text{dom } \Phi$, is the set of vectors $x$ such that $\Phi(x) \neq \emptyset$; a multifunction is polyhedral if its graph, i.e., the set $\{(x, y) : y \in \Phi(x)\}$, is the union of finitely many polyhedra.

**Lemma 5.1** (Gowda) A single-valued polyhedral multifunction $\Phi$ is piecewise linear on its domain; i.e., $\Phi$ is continuous on $\text{dom } \Phi$ and there exist finitely many linear functions $\{f^1, \dots, f^k\}$ such that $\Phi(x) \in \{f^1(x) \dots, f^k(x)\}$ for all $x \in \text{dom } \Phi$. If $\text{dom } \Phi$ is convex, then $\Phi$ is Lipschitz continuous on $\text{dom } \Phi$.

**Proposition 5.1** Let $f : \mathbb{R}^n \to \mathbb{R}^n$ be Lipschitz and let $K \subseteq \mathbb{R}^m$ be polyhedral. A necessary and sufficient condition for (5.3) to have a unique $C^1$ solution trajectory $x(t)$ on $[0, \infty)$ for all $x^0 \in \mathbb{R}^n$ is that $\text{BSOL}(K, Cx, D)$ is a singleton for all $x \in \mathbb{R}^n$.

*Proof* The necessity is obvious; indeed, for any $x^0 \in \mathbb{R}^n$ and any $u \in \text{SOL}(K, Cx^0, D)$, we must have $Bu = \dot{x}(0) - f(x^0)$, which is necessarily unique because the trajectory $x(t)$ is uniquely determined by $x^0$ by assumption. Conversely, if $\text{BSOL}(K, Cx, D)$ is a singleton for all $x \in \mathbb{R}^n$, then since the set-valued map $x \mapsto \text{SOL}(K, Cx, D)$ is a polyhedral multifunction [82,38], by Lemma 5.1, the single-valued map $x \mapsto \text{BSOL}(K, Cx, D)$ is Lipschitz continuous. The existence and uniqueness of a $C^1$ solution trajectory $x(t)$ of (5.3) on $[0, \infty)$ follows from classical ODE theory. $\quad \square$

A class of problems (5.3) satisfying the assumptions of Proposition 5.1 without the pair $(K, D)$ being coherently oriented occurs when $K$ is a compact polyhedron, $D \equiv E^TGE$ for some (possibly asymmetric) positive definite matrix $G \in \mathbb{R}^{k \times k}$ and an arbitrary matrix $E \in \mathbb{R}^{k \times m}$, and $B \equiv HE$ some matrix $H \in \mathbb{R}^{n \times k}$. For such a matrix $D$, it is known [38] that $\text{ESOL}(K, r, D)$, thus $\text{BSOL}(K, r, D)$, is a singleton for all vectors $r \in \mathbb{R}^m$. The above proposition is the basis of many results for the LCS that are obtained subsequently to the completion of the present paper; see [22,72,88].
## 5.3 Locally smooth solutions

Global properties like ($A_0$) or coherent orientation for an affine pair are desirable but not always satisfied in applications. In what follows, we investigate the conversion of (2.2) into a well-behaved ODE locally in a suitable neighborhood of a triple $(0, x^0, u^0)$, where $u^0$ is a given solution of the VI $(K, F(0, x^0, \cdot))$. Like the last subsection, the discussion herein is restricted to the initial-value version of the DVI. Specifically, we are interested in identifying sufficient conditions on an initial state $x^0$ under which a locally smooth solution of the DVI (2.2) exists that satisfies the initial condition $x(0) = x^0$.

We begin by postulating a well-known property of a solution to a VI: $(A'_0)$ the VI $(K, F(0, x^0, \cdot))$ has a strongly regular solution $u^0$ [81,38]. Under $(A'_0)$ and $(B_0)$, it follows that there exist a scalar $\epsilon > 0$, a neighborhood $\mathcal{N}$ of $x^0$, a neighborhood $\mathcal{U}$ of $u^0$, and a Lipschitz continuous function $u: [0, \epsilon] \times \mathcal{N} \to K \cap \mathcal{U}$ such that for each pair $(t, x) \in [0, \epsilon] \times \mathcal{N}$, $u(t, x)$ is the unique solution of the VI $(K, F(t, x, \cdot))$ that belongs to $\mathcal{U}$. Therefore in the domain $[0, \epsilon] \times \mathcal{N}$, the DVI (2.2) is once again equivalent to an ODE with a Lipschitz right-hand side, and the existence of a locally smooth solution to the DVI satisfying the initial condition $x(0) = x^0$ therefore follows from standard ODE theory [26].

It is possible to refine the condition $(A'_0)$ in the case when the set $K$ is the nonnegative orthant $\mathbb{R}_+^m$ so that the VI $(K, F(t, x, \cdot))$ is equivalent to the nonlinear complementarity problem (NCP):

$$0 \le u \perp F(t, x, u) \ge 0.$$

(The discussion below can be generalized to mixed complementarity constraints, in particular, to the DCP (2.10) that is derived from a VIE with a finitely representable set. For simplicity, we restrict our treatment to the above NCP.) Let $u^0$ satisfy

$$0 \le u^0 \perp F(0, x^0, u^0) \ge 0. \qquad (5.4)$$

Define three fundamental index sets associated with the (initial) solution $u^0$:

$$\alpha = \{i : u_i^0 > 0 = F_i(0, x^0, u^0)\} \text{ the support of } u^0$$

$$\beta = \{i : u_i^0 = 0 = F_i(0, x^0, u^0)\} \text{ the degenerate set of } u^0$$

$$\gamma = \{i : u_i^0 = 0 < F_i(0, x^0, u^0)\} \text{ the inactive set of } u^0.$$

For any two index subsets $\alpha_1$ and $\alpha_2$ of $\{1, \dots, m\}$, let $J_{u_{\alpha_1}} F_{\alpha_2}$ denote the sub-matrix of the (partial) Jacobian matrix $J_u F$ of $F$ with respect to the $u$-variable with rows and columns indexed by $\alpha_2$ and $\alpha_1$, respectively; i.e.,

$$J_{u_{\alpha_1}} F_{\alpha_2} \equiv \left( \frac{\partial F_i}{\partial u_j} : (i, j) \in \alpha_2 \times \alpha_1 \right).$$
The strong regularity of $u^0$ is equivalent to two conditions [81, 38]: (a) the matrix $J_{u\alpha}F_\alpha(0, x^0, u^0)$ is nonsingular, and (b) the Schur complement

$$J_{u\beta}F_{\beta}(0, x^0, u^0) - J_{u\alpha}F_{\beta}(0, x^0, u^0) \left(J_{u\alpha}F_{\alpha}(0, x^0, u^0)\right)^{-1} J_{u\beta}F_{\alpha}(0, x^0, u^0)$$

of $J_{u\alpha}F_\alpha(0, x^0, u^0)$ in

$$\begin{bmatrix} J_{u\alpha}F_{\alpha}(0, x^0, u^0) & J_{u\beta}F_{\beta}(0, x^0, u^0) \\ J_{u\alpha}F_{\beta}(0, x^0, u^0) & J_{u\beta}F_{\beta}(0, x^0, u^0) \end{bmatrix}$$

is a P-matrix (see [30] for an extensive discussion of the role of a P-matrix in linear complementarity theory).

To deal with the case where the matrix $J_{u\alpha}F_\alpha(0, x^0, u^0)$ is singular (therefore $u^0$ is not strongly regular), we consider the extreme situation where $F_\alpha$ is independent of the variable $u$ and write

$$\Phi_{\alpha}(t, x) \equiv F_{\alpha}(t, x, u).$$

(The other components of $F$ can still depend on $u$.) We define the function

$$\tilde{F}_{\alpha}(t, x, u) \equiv \frac{\partial \Phi_{\alpha}(t, x)}{\partial t} + J_{x} \Phi_{\alpha}(t, x) f(t, x, u),$$

which is obtained by differentiating $\Phi_\alpha$ with respect to $t$ and then substituting $f(t, x, u)$ for $\dot{x}$ (provided that the latter substitution is legal). Using the function $\tilde{F}_\alpha(t, x, u)$ instead of the original $F_\alpha(t, x, u)$, which does not contain the variable $u$, we postulate a strong regularity hypothesis on a transformed NCP. From the latter, we obtain an equivalent ODE with a Lipschitz right-hand side whose solution can be shown to yield a smooth solution to the initial-value differential NCP (DNCP):

$$\begin{array}{l} \dot{x} = f(t, x, u), \quad x(0) = x^0, \\ 0 \le u \perp F(t, x, u) \ge 0 \end{array} \qquad (5.5)$$

locally near $(0, x^0, u^0)$. This postulate is stated in terms of two matrix-theoretic conditions analogous to the two conditions in the strongly regular case. The result below bypasses the equivalent ODE and addresses directly the existence of a locally smooth solution of (5.5).

**Theorem 5.2** Let $f$ and $F$ be continuously differentiable, and let $u^0$ satisfy (5.4). Suppose that $F_\alpha(t, x, u) = \Phi_\alpha(t, x)$ is independent of $u$ and that $\tilde{F}_\alpha(0, x^0, u^0) = 0$. If the matrix $J_{u\alpha}\tilde{F}_\alpha(0, x^0, u^0) = J_x\Phi_\alpha(0, x^0)J_{u\alpha}f(0, x^0, u^0)$ is nonsingular, and the following Schur complement is a P-matrix:

$$\begin{gathered} J_{u\beta}F_{\beta}(0, x^0, u^0) \\ -J_{u\alpha}F_{\beta}(0, x^0, u^0) \left(J_{x}\Phi_{\alpha}(0, x^0)J_{u\alpha}f(0, x^0, u^0)\right)^{-1} J_{x}\Phi_{\alpha}(0, x^0)J_{u\beta}f(0, x^0, u^0) \end{gathered}$$

then the DNCP (5.5) has a locally smooth solution near $(0, x^0, u^0)$.
Proof Consider the modified DNCP:

$$ 
\begin{array}{l}
\dot{x} = f(t, x, u), \quad x(0) = x^0, \\
0 \le u \perp \hat{F}(t, x, u) \ge 0,
\end{array}
\qquad (5.6)
$$

where the components of $\hat{F}$ are equal to the corresponding components of $F$ except for those in $\alpha$ which are equal to the ones in $\tilde{F}_\alpha$. By assumption, $u^0$ satisfies the three complementarity conditions at $(0, x^0)$. It is easy to see that the two matrix-theoretic assumptions are precisely the strong regularity property of $u^0$ as a solution of the modified NCP at $(0, x^0)$. Therefore, (5.6) is equivalent to an ODE with a Lipschitz right-hand side, and thus has a locally smooth solution $(x^*(t), u^*(t))$ with $u^*(0) = u^0$ on the interval $[0, \delta]$ for some $\delta > 0$. Moreover, provided that $\delta > 0$ is sufficiently small, $u_\alpha^*(t) > 0$ on $[0, \delta]$. Hence it remains to show that $F_\alpha(t, x^*(t), u^*(t)) = \Phi_\alpha(t, x^*(t)) = 0$ on $[0, \delta]$. We have $\tilde{F}_\alpha(t, x^*(t), u^*(t)) = \frac{d}{dt}\Phi_\alpha(t, x^*(t)) = 0$ for all $t \in [0, \delta]$. Since $\Phi_\alpha(0, x^*(0)) = 0$, simple integration yields $\Phi_\alpha(t, x^*(t)) = 0$ on $[0, \delta]$, completing the proof. $\square$

## 5.4 Linear complementarity systems

For the LCS (2.8), the treatment can be further extended by exploiting the properties of the linear complementarity condition:

$$ 0 \le u \perp q + Cx + Du \ge 0. $$

We continue to focus on the case where a strongly regular solution to the LCP $(r^0, D)$, where $r^0 = q + Cx^0$, is not readily available. Before proceeding, we should note that the two theses [21,49] and the subsequent paper [52] have investigated the existence of locally smooth solutions to an LCS in connection with a “regular initial state”. In particular, under the P-property on a “leading row coefficient matrix”, a lexicography condition is shown to be both necessary and sufficient for an initial state to be regular [52, Theorem 6.8]. In contrast, our development below is inspired by the notion of the index of a DAE.

Specifically, we describe a differentiation procedure that attempts to uncover a strongly regular solution to a certain modified LCS, from which a locally smooth solution of the original LCS (2.8), if it exists, can be determined. We begin by obtaining a solution $u^0$ of the LCP $(r^0, D)$. Assume that we can identify a set of nonzero rows in the principal submatrix

$$ \begin{bmatrix} D_{\alpha\alpha} & D_{\alpha\beta} \\ D_{\beta\alpha} & D_{\beta\beta} \end{bmatrix} $$

that are linearly dependent on another set of rows in the same matrix, where $\alpha$ is the support of $u^0$ and $\beta$ is the degenerate index set of $u^0$ as a solution of the LCP $(r^0, D)$. Note that $\beta \neq \emptyset$ implies in particular that $u^0$ is not strongly
regular. For simplicity, we assume specifically that a nonsingular matrix $P_{\alpha\alpha}$ can
be found such that

$$
P_{\alpha\alpha}[D_{\alpha\alpha} \quad D_{\alpha\beta}] = \begin{bmatrix} E_{\alpha_1\alpha} & E_{\alpha_1\beta} \\ 0 & 0 \end{bmatrix} \in \Re^{|\alpha| \times (|\alpha| + |\beta|)}
$$

for some subset $\alpha_1$ (with complement $\alpha_2$) of $\alpha$, where the two zero matrices are
of order $\alpha_2 \times \alpha$ and $\alpha_2 \times \beta$, respectively. (For instance, if $D$ is the zero matrix,
then $P_{\alpha\alpha}$ can be taken to be the identity matrix.) Setting $u_\gamma = 0$, premultiplying
the equation $0 = C_\alpha x + D_{\alpha\alpha} u_\alpha + D_{\alpha\beta} u_\beta$ by $P_{\alpha\alpha}$, we obtain

$$
0 = P_{\alpha 2 \alpha} C_{\alpha} x, \quad (5.7)
$$

in which the variable $u$ has been eliminated. Thus, we can differentiate the latter
equation with respect to time and then substitute $\dot{x}$ using the first equation in
(2.8). This process results in the following modified LCS of reduced size:

$$
\dot{x} = Ax + B_{\cdot\alpha}u_{\alpha} + B_{\cdot\beta}u_{\beta}, \quad x(0) = x^0
$$

$$
0 \le \begin{pmatrix} u_{\alpha_1} \\ u_{\alpha_2} \\ u_{\beta} \end{pmatrix} \perp \begin{bmatrix} P_{\alpha_1 \alpha} C_{\alpha} \\ P_{\alpha_2 \alpha} C_{\alpha} A \\ C_{\beta} \end{bmatrix} x + \begin{bmatrix} E_{\alpha_1 \alpha} & E_{\alpha_1 \beta} \\ P_{\alpha_2 \alpha} C_{\alpha} B_{\cdot \alpha} & P_{\alpha_2 \alpha} C_{\alpha} B_{\cdot \beta} \\ D_{\beta \alpha} & D_{\beta \beta} \end{bmatrix} \begin{pmatrix} u_{\alpha} \\ u_{\beta} \end{pmatrix} \ge 0.
$$

If this modified LCS has a locally smooth solution $(x^*(t), u_\alpha^*(t), u_\beta^*(t))$ on an interval $[0, \delta]$ such that

$$
\begin{bmatrix} P_{\alpha_1 \alpha} C_{\alpha} \\ P_{\alpha_2 \alpha} C_{\alpha} A \end{bmatrix} x^*(t) + \begin{bmatrix} E_{\alpha_1 \alpha} & E_{\alpha_1 \beta} \\ P_{\alpha_2 \alpha} C_{\alpha} B_{\cdot \alpha} & P_{\alpha_2 \alpha} C_{\alpha} B_{\cdot \beta} \end{bmatrix} \begin{pmatrix} u_{\alpha}^*(t) \\ u_{\beta}^*(t) \end{pmatrix} = 0
$$

and $C_{\gamma\alpha}x^*(t) + D_{\gamma\alpha}u_{\alpha}^*(t) + D_{\gamma\beta}u_{\beta}^*(t) \ge 0$ on $[0,t]$, with $u_{\gamma}^* \equiv 0$, $(x^*, u^*)$ is then a locally smooth solution of the original LCS (2.8). We give in the next result a sufficient condition for the successful identification of such a solution $(x^*(t), u_{\alpha}^*(t), u_{\beta}^*(t))$.

**Theorem 5.3** In the above setting, suppose that $P_{\alpha_2\alpha}C_\alpha(Ax^0+Bu^0)=0$. Assume further that the matrix

$$
\begin{bmatrix} E_{\alpha_1 \alpha} \\ P_{\alpha_2 \alpha} C_{\alpha} B_{\cdot \alpha} \end{bmatrix} \qquad (5.8)
$$

is nonsingular and that the Schur complement

$$
D_{\beta\beta} - D_{\beta\alpha} \begin{bmatrix} E_{\alpha_1\alpha} \\ P_{\alpha_2\alpha} C_{\alpha} B_{\cdot\alpha} \end{bmatrix}^{-1} \begin{bmatrix} E_{\alpha_1\beta} \\ P_{\alpha_2\alpha} C_{\alpha} B_{\cdot\beta} \end{bmatrix} \quad (5.9)
$$

is a P-matrix. The LCS (2.8) has a locally smooth solution near (0, $x^0$, $u^0$).
*Proof* The assumptions imply that the modified LCS has a locally smooth solution trajectory ($x^*(t), u_\alpha^*(t), u_\beta^*(t)$) on an interval $t \in [0, \delta]$ that emanates from the triple ($x^0, u_\alpha^0, u_\beta^0$); moreover, ($u_\alpha^*(t), u_\beta^*(t)$) is Lipschitz continuous and $u_\alpha^*(t) > 0$ on $[0, \delta]$. This implies that

$$P_{\alpha_1\alpha} \left[ C_{\alpha\cdot} x^*(t) + D_{\alpha\alpha} u_\alpha^*(t) + D_{\alpha\beta} u_\beta^*(t) \right] = 0, \quad (5.10)$$

$$P_{\alpha_2\alpha} C_{\alpha\cdot} \left[ A x^*(t) + B_{\alpha\cdot} u_\alpha^*(t) + B_{\beta\cdot} u_\beta^*(t) \right] = 0 \quad (5.11)$$

on the domain of the trajectory. From (5.11), we obtain $0 = P_{\alpha_2\alpha} C_{\alpha\cdot} \dot{x}^*(t)$ on $[0, \delta]$, which implies $P_{\alpha_2\alpha} \left[ C_{\alpha\cdot} x^*(t) + D_{\alpha\alpha} u_\alpha^*(t) + D_{\alpha\beta} u_\beta^*(t) \right] = 0$. This together with the equation (5.10) yields $C_{\alpha\cdot} x^*(t) + D_{\alpha\alpha} u_\alpha^*(t) + D_{\alpha\beta} u_\beta^*(t) = 0$, by the non-singularity of $P_{\alpha\alpha}$. By restricting $\delta$ if necessary, we have $C_{\gamma\cdot} x^*(t) + D_{\gamma\alpha} u_\alpha^*(t) + D_{\gamma\beta} u_\beta^*(t) > 0$ on $[0, \delta]$, by continuity of the solution trajectory. Consequently, by letting $u_\gamma^*(t) \equiv 0$ on $[0, \delta]$, we obtain a desired locally smooth solution of the original LCS (2.8). □

The crux of the above procedure is that we attempt to obtain, by means of elementary matrix row operations on some active constraints, an equation that contains only the state variable $x$, such as (5.7), which we can then differentiate, thereby obtaining a new LCS. If the latter LCS has a strongly regular solution that respects the active constraints employed in the row operations, then we have successfully identified a locally smooth solution of the original LCS. Otherwise, we need to repeat the procedure and try to derive another modified LCS with the desired solution property.

## 6 The approach of convex-valued DIs

When the reduction to an ODE with a Lipschitz right-hand side fails, the next best thing is to be able to apply the theory of DIs. Although there is an extensive literature on DIs, much of it requires properties on the set-valued map that translate into very restrictive requirements for the DVI. One exception is the following result, which pertains to a DI with a set-valued map that has closed and convex values and satisfies a linear growth property (6.1). This result will be the cornerstone for our development in this section. Throughout the rest of the paper, we let $\Omega = [0, T] \times \mathbb{R}^n$ and let $\mathcal{B}$ denote the closed Euclidean unit ball centered at the origin.

**Lemma 6.1** [32, Theorem 5.1] Let $\mathbf{F}: \Omega \to \mathbb{R}^n$ be an upper semicontinuous set-valued map with nonempty closed convex values. Suppose that there exists a scalar $\rho_{\mathbf{F}} > 0$ satisfying

$$\sup\{ \|y\| : y \in \mathbf{F}(t,x) \} \le \rho_{\mathbf{F}} (1 + \|x\|), \quad \forall (t,x) \in \Omega. \quad (6.1)$$
For every $x^0 \in \mathbb{R}^n$, the DI: $\dot{x} \in \mathbf{F}(t, x)$, $x(0) = x^0$ has a weak solution in the sense of Carathéodory.

We should note that under the linear growth assumption (6.1), the upper semicontinuity of $\mathbf{F}$ holds if $\mathbf{F}$ is closed (that is, if the graph of $\mathbf{F}$ is closed). Indeed, this assumption implies that $\mathbf{F}$ is bounded on compact sets; therefore the closedness of $\mathbf{F}$ is equivalent to its upper semicontinuity.

## 6.1 Setting of the DVI

We wish to apply Lemma 6.1 to the following initial-value DVI with the right-hand function in the ODE being affine in $u$ and with the VI function being separable in $u$:

$$ \begin{gathered} \dot{x} = f(t, x) + B(t, x)u, \quad x(0) = x^0, \tag{6.2} \\ u \in \text{SOL}(K, G(t, x) + F(\cdot)), \end{gathered} $$

where $K \subseteq \mathbb{R}^m$ is closed and convex, $(f, B, G) : \Omega \to \mathbb{R}^n \times \mathbb{R}^{n \times m} \times \mathbb{R}^m$, and $F : \mathbb{R}^m \to \mathbb{R}^m$ are given functions. We let $K_\infty$ be the recession cone of $K$; see [85]. From this point on, the function $F$ has only one argument.

We should say a few words about the distinguished features of the above setting, namely, the linearity in $u$ in the right-hand side of the ODE and the separability of $u$ from $(t, x)$ in the VI. First of all, linearity in $u$ is the setting commonly used in the control literature; this special feature greatly facilitates the analysis (e.g., to establish the convexity of the set-valued map $\mathbf{F}$ in (6.4) when the VIs have convex solution sets). Secondly, the special structure imposed herein indeed arises in many situations, including the LCS (2.8), the two DCPs (2.9) and (2.10) derived from a VIE, and the DNG (see (4.3) and (4.4)). Finally, by means of a sequential linearization scheme, the problem (6.2) could be used as the workhorse in an iterative procedure for solving the general case.

We postulate two blanket assumptions (A) and (B) on the triplet of functions $(f, B, G)$ that will be enforced throughout the rest of the paper. Conditions on $F$ will vary and will be stated as needed. Some results will have restrictions on $x^0$.

(A) $f, B$, and $G$ are Lipschitz continuous functions on $\Omega$ with Lipschitz constants $L_f, L_B$, and $L_G$, respectively;

(B) $B$ is bounded on $\Omega$ with $\sigma_B = \sup_{(t,x) \in \Omega} \|B(t,x)\| < \infty$.

Notice that (A) implies that $f$ and $G$ have linear growth on $\Omega$ in $x$; i.e., for some positive constants $\rho_f$ and $\rho_G$ and for all $(t, x) \in \Omega$,

$$ \|f(t, x)\| \le \rho_f (1 + \|x\|) \quad \text{and} \quad \|G(t, x)\| \le \rho_G (1 + \|x\|). \tag{6.3} $$

The assumption (A) can be replaced by weaker assumptions, such as $(A')$ in Section 7, which is more in the spirit of Carathéodory. For reasons given at the end of Section 3, we refrain from introducing such weaker assumptions at this stage.
Our goal is to show, under assumptions (A) and (B) and further assumptions to be imposed, that the set-valued map: for $(t, x) \in \Omega$,

$$ \mathbf{F}(t, x) = \{f(t, x) + B(t, x)u : u \in \mathrm{SOL}(K, G(t, x) + F)\}, \quad (6.4) $$

satisfies the assumptions in Lemma 6.1, thereby establishing the existence of a weak solution to the initial-value DVI (6.2); see Proposition 6.1. It is clear that if $F$ is monotone, then $\mathbf{F}(t, x)$ is a convex set; however, the converse is false. Proposition 6.5 identifies a case where $\mathbf{F}(t, x)$ is convex without $F$ being monotone. The lemma below shows that the remaining conditions on the map $\mathbf{F}$ are satisfied if the solutions to the VI have linear growth.

**Lemma 6.2** Let $(f, B, G)$ satisfy (A) and (B) above. Let $F: \mathbb{R}^m \to \mathbb{R}^m$ be continuous. Suppose that $\mathrm{SOL}(K, q + F) \neq \emptyset$ for all $q \in G(\Omega)$, and there exists a constant $\rho > 0$ such that for all $q \in G(\Omega)$,

$$ \sup\{\|u\| : u \in \mathrm{SOL}(K, q + F)\} \leq \rho (1 + \|q\|). \quad (6.5) $$

Then a constant $\rho_{\mathbf{F}} > 0$ exists such that (6.1) holds for the map $\mathbf{F}$ defined by (6.4). Hence, $\mathbf{F}$ is upper semicontinuous, and thus closed-valued, on $\Omega$.

*Proof* It is easy to see that

$$ \sup\{\|y\| : y \in \mathbf{F}(t, x)\} \leq [\rho_f + \rho \sigma_B (1 + \rho_G)] (1 + \|x\|), \quad \text{on } \Omega. $$

Hence $\mathbf{F}$ has linear growth. To prove the upper semicontinuity of $\mathbf{F}$, it suffices to show that $\mathbf{F}$ is closed. Suppose that $\{(t_k, x^k)\} \subset \Omega$ is a sequence converging to some vector $(t_\infty, x^\infty)$ in $\Omega$ and $\{f(t_k, x^k) + B(t_k, x^k)u^k\}$ converges to some vector $z^\infty \in \mathbb{R}^n$ as $k \to \infty$, where $u^k \in \mathrm{SOL}(K, G(t_k, x^k) + F)$ for every $k$. It follows that the sequence $\{u^k\}$ is bounded, and thus has a convergent subsequence with a limit $u^\infty$, which is easily seen to be a solution of the limiting VI $(K, G(t_\infty, x^\infty) + F)$, by a simple limiting argument. Thus, $z^\infty = f(t_\infty, x^\infty) + B(t_\infty, x^\infty)u^\infty$ belongs to $\mathbf{F}(t_\infty, x^\infty)$, completing the proof. $\square$

The linear growth condition (6.5) holds trivially for a bounded set $K$ (e.g., in a “box-constrained” VI). In the next subsection, we derive sufficient conditions for this growth property to hold, focusing on the case of an unbounded set $K$ (such as the nonnegative orthant). For now, we state an important result to show that a weak solution to the DI formulation of the DVI indeed yields a weak solution to the DVI as defined in Sect. 2.1. The following measurable selection lemma is a version of Filippov’s implicit function theorem [42]; see also Theorem 2.3 in [89] and Corollary 1 in [10]. The reader can consult these references for a proof.

**Lemma 6.3** Let $h: \Omega \times \mathbb{R}^m \to \mathbb{R}^n$ be a continuous function and $U: \Omega \to \mathbb{R}^m$ be a closed set-valued map such that for some constant $\eta_U > 0$,

$$ \sup_{u \in U(t, x)} \|u\| \leq \eta_U (1 + \|x\|), \quad \forall (t, x) \in \Omega. $$
Let $v: [0, T] \to \mathbb{R}^n$ be a measurable function and $x: [0, T] \to \mathbb{R}^n$ be a continuous function satisfying $v(t) \in h(t, x(t), U(t, x(t)))$ for almost all $t \in [0, T]$. There exists a measurable function $u: [0, T] \to \mathbb{R}^m$ such that $u(t) \in U(t, x(t))$ and $v(t) = h(t, x(t), u(t))$ for almost all $t \in [0, T]$.

Based on the above lemma, we can establish the following connection between the weak solution of a DVI and that of its DI formulation.

**Proposition 6.1** Let $(f, B, G)$ satisfy conditions (A) and (B) above. Let $F: \mathbb{R}^m \to \mathbb{R}^m$ be continuous. Suppose there exists a constant $\rho > 0$ such that (6.5) holds for all $q \in G(\Omega)$. If the DI: $\dot{x} \in F(t, x)$, $x(0) = x^0$ has a weak solution in the sense of Carathéodory, where $F$ is given by (6.4), then the initial-value DVI (6.2) has a weak solution in the same sense.

*Proof* By the integral representation of the DI [10, Lemma 1] and the linear growth of $F$, we deduce that if $x$ is a weak solution of the DI, then

$$ \|x(t)\| \le \|x^0\| + \int_0^t \rho_F (1 + \|x(s)\|) \, ds. $$

By Gronwall's lemma, we obtain

$$ \|x(t)\| \le (\|x^0\| + \rho_F T) \exp(\rho_F T). $$

Consequently, letting $U(t, x) \equiv \text{SOL}(K, G(t, x) + F)$ and $h(t, x, u) \equiv f(t, x) + B(t, x)u$, we can apply Lemma 6.3 to deduce the existence of a measurable function $u(t)$ such that $u(t) \in \text{SOL}(K, G(t, x(t)) + F)$ and $\dot{x}(t) = f(t, x(t)) + B(t, x(t))u(t)$ for almost all $t$. Since $\|u(t)\|$ is bounded on $[0, T]$ almost everywhere, it follows that $u$ is integrable on $[0, T]$. Clearly for any continuous function $\tilde{u}: [0, T] \to K$, we must have

$$ \int_0^T (\tilde{u}(t) - u(t))^T (G(t, x(t)) + F(u(t))) dt \ge 0, $$

because the integrand is nonnegative for almost all $t$. This shows that the pair $(x, u)$ is a weak solution of the DVI. □

## 6.2 The linear growth of VI solutions

Although there is an extensive literature on the VI, including the study of boundedness of solutions (see [38]), the linear growth property (6.5) has not been a focused topic of attention. Our results in this subsection, several of which are new, contribute directly to the theory of the static VI. We recall some function classes from VI theory [38, Definition 2.3.9]. We say that a mapping
$H : \mathbb{R}^m \to \mathbb{R}^m$ is monotone plus on a convex set $K \subseteq \mathbb{R}^m$ if $H$ is monotone on $K$, i.e., $(u - u')^T (H(u) - H(u')) \ge 0$, for all $u$ and $u' \in K$ and the following "plus" property holds:

$$[u, u' \in K, (u - u')^T (H(u) - H(u')) = 0] \Rightarrow H(u) = H(u').$$

It is known that if $H$ is a monotone map on the closed convex set $K$, then $\text{SOL}(K, H)$ is a convex set (possibly empty). Moreover, if $H$ is monotone plus on $K$, then $H(u)$ is a constant vector for all $u \in \text{SOL}(K, H)$, provided that the latter set is nonempty; see [38]. A large subclass of the monotone plus mappings is the class of "strongly monotone composite" maps. Specifically, $H$ is a strongly monotone composite map on $K$ if $H \equiv E^T \circ \Psi \circ E$ for some matrix $E \in \mathbb{R}^{l \times m}$ and a mapping $\Psi : \mathbb{R}^l \to \mathbb{R}^l$ that is strongly monotone on $EK \subseteq \mathbb{R}^l$; that is, there is a constant $\eta_\psi > 0$ where

$$(u' - u)^T (\Psi(u') - \Psi(u)) \ge \eta_\psi \|u' - u\|^2, \quad \forall u', u \in EK.$$

The following result pertains to the existence and linear growth of solutions to a VI under a coercivity property on its defining function. This result includes the case where $K$ is bounded so that the limit (6.6) holds vacuously.

**Proposition 6.2** Let $K \subseteq \mathbb{R}^m$ be a nonempty closed convex set. Let $F : \mathbb{R}^m \to \mathbb{R}^m$ be a continuous map. Suppose that there exists $u^{\text{ref}} \in K$ satisfying

$$\liminf_{\substack{u \in K \\ \|u\| \to \infty}} \frac{(u - u^{\text{ref}})^T F(u)}{\|u\|^2} > 0. \qquad (6.6)$$

For every vector $r \in \mathbb{R}^m$, $\text{SOL}(K, r + F) \ne \emptyset$; moreover, a constant $\rho > 0$ exists such that (6.5) holds for all $r \in \mathbb{R}^m$. Finally, if $F$ is monotone, then $\text{SOL}(K, r + F)$ is a convex set.

*Proof* It suffices to prove the second assertion. Assume for the sake of contradiction that no such scalar $\rho$ exists. There exist sequences $\{r^k\}$ and $\{u^k\}$ such that, for every $k$,

$$(u - u^k)^T (r^k + F(u^k)) \ge 0, \quad \forall u \in K$$

and $\|u^k\| > k(1 + \|r^k\|)$. The sequence $\{u^k\}$ clearly is unbounded. Moreover, we have

$$(u^{\text{ref}} - u^k)^T (r^k + F(u^k)) \ge 0$$

Dividing by $\|u^k\|^2$ and letting $k \to \infty$, we obtain a contradiction to the coercivity property (6.6). □
If $K$ is the Cartesian product of finitely many lower-dimensional closed convex sets (see (5.1)), Proposition 6.2 remains valid with the condition (6.6) weakened as follows:

$$ \liminf_{\substack{u \in K \\ \|u\| \to \infty}} \frac{\max_{1 \le v \le N} F_v(u)^T (u - u^{\text{ref}})_v}{\|u\|^2} > 0. $$

All the results in this section have a similar weakening when $K$ has the Cartesian structure (5.1); we will not repeat this remark henceforth.

The next result pertains to a VI with a strongly monotone composite map. Such a map is not necessarily coercive in the sense of the last proposition, and yet the linear growth property remains valid.

**Proposition 6.3** Let $K \subseteq \mathbb{R}^m$ be a nonempty closed convex set. Let $F \equiv E^T \circ \Psi \circ E$, where $E \in \mathbb{R}^{l \times m}$ and $\Psi: \mathbb{R}^l \to \mathbb{R}^l$ is continuous and strongly monotone on $EK$. If $K_\infty \cap \ker E = \{0\}$, then SOL($K, r + F$) is a nonempty convex set for all $r \in \mathbb{R}^m$; moreover, there exists $\rho > 0$ such that (6.5) holds for all $r \in \mathbb{R}^m$.

*Proof* By Theorem 2.3.16 in [38], the nonemptiness of SOL($K, r + F$) will follow if we can show that $K_\infty \cap [-(r + F(K))^*] = \{0\}$. For this purpose let $v$ be a nonzero element of $K_\infty \cap [-(r + F(K))^*]$. We then have, for any $x \in K$ and $\tau > 0$, $x + \tau v \in K$, so

$$ v^T r + (Ev)^T \Psi(Ex + \tau Ev) \le 0, \quad \forall \tau > 0. $$

By the strong monotonicity of $\Psi$, a constant $\eta_\psi > 0$ exists such that

$$ \eta_\psi \tau^2 \|Ev\|^2 \le \tau (Ev)^T (\Psi(Ex + \tau Ev) - \Psi(Ex)) \\ \le -\tau v^T r - \tau (Ev)^T \Psi(Ex). $$

Since this holds for all $\tau > 0$, we must have $Ev = 0$, which means that $0 \ne v \in K_\infty \cap \ker E$. This contradiction establishes the nonemptiness of SOL($K, r + F$), which must necessarily be convex because $F$ is monotone. The proof of the linear growth of the solutions to the VI ($K, r + F$) is similar to that of Proposition 6.2, using the condition $K_\infty \cap \ker E = \{0\}$ instead of (6.6). Details are not repeated. $\square$

In the case of an affine mapping $F$, the coercivity condition (6.6) basically restricts $F$ to be “strongly copositive” on the recession cone of $K$. In what follows, we derive some similar results for the VI ($K, r, D$) with a positive semidefinite matrix $D$ by relaxing such a coercivity condition. In addition to the properties of an AVI already used in the proof of Proposition 5.1, we will need some further facts that we now summarize; for details, see [38]. Let $D$ be an $m \times m$ matrix. The solution set of the following homogeneous CP is called the
VI kernel of the pair $(K, D)$ and is denoted by $\mathcal{K}(K, D)$:

$$ K_{\infty} \ni v \perp Dv \in (K_{\infty})^*, \qquad (6.7) $$

where $(K_{\infty})^*$ is the dual cone of $K_{\infty}$. Nonzero elements of $\mathcal{K}(K, D)$ are known to have an important role to play in studying the VIs associated with the pair $(K, D)$. In particular, $(K, D)$ is called an $R_0$ pair when $\mathcal{K}(K, D) = \{0\}$. We denote the interior of the dual cone of the VI kernel $\mathcal{K}(K, D)$ by $\text{int} \, \mathcal{K}(K, D)^*$. When $D$ is positive semidefinite, the following properties hold:

(a) $(D + D^T)\text{SOL}(K, r, D)$ is a singleton for all $r \in \mathbb{R}^m$ [38, Lemma 2.4.14];

(b) $\mathcal{K}(K, D) = \{v \in K_{\infty} : -D^T v \in (K_{\infty})^*\} = K_{\infty} \cap (-DK_{\infty})^*$ [38, Lemma 2.5.2]; in particular, $\mathcal{K}(K, D)$ is a polyhedral cone if $K$ is polyhedral;

(c) if $K$ is a cone (polyhedrality is not needed), then $\text{SOL}(K, r, D)$ is nonempty and bounded for every $r \in \text{int} \, \mathcal{K}(K, D)^*$ [38, Proposition 2.5.11];

(d) if $K$ is a polyhedron containing no lines and $(K, D)$ is an $R_0$ pair, then $\text{SOL}(K, r, D)$ is nonempty and bounded for every $r \in \mathbb{R}^m$ [38, Theorem 2.5.20].

Pertaining to the VI $(K, r, D)$ with a positive semidefinite $D$, the following lemma extends the above properties in several ways; most importantly, it establishes the linear growth of solutions to the VI.

**Proposition 6.4** Let $K \subseteq \mathbb{R}^m$ be a nonempty closed convex set containing the origin. Let $D$ be an $m \times m$ positive semidefinite matrix. For every vector $r \in \text{int} \, \mathcal{K}(K, D)^*$, $\text{SOL}(K, r, D)$ is a nonempty convex set on which $(D + D^T)u$ is a constant for all $u \in \text{SOL}(K, r, D)$. Furthermore,

(a) if $(K, D)$ is an $R_0$ pair, then a constant $\rho > 0$ exists such that for all $r \in \mathbb{R}^m$, $\text{SOL}(K, r, D) \neq \emptyset$ and

$$ \sup\{\|u\| : u \in \text{SOL}(K, r, D)\} \leq \rho (1 + \|r\|). \qquad (6.8) $$

(b) if $K$ is a polyhedron, then (6.8) holds for all $r \in \text{int} \, \mathcal{K}(K, D)^*$; moreover, the (single-valued) map $r \mapsto (D + D^T)\text{SOL}(K, r, D)$ is Lipschitz continuous on $\mathcal{K}(K, D)$; i.e., there exists a constant $L_D > 0$ such that

$$ \| (D + D^T) (u - u') \| \leq L_D \| r - r' \|, $$

for all $u \in \text{SOL}(K, r, D)$, $u' \in \text{SOL}(K, r', D)$, and $r$ and $r'$ in $\mathcal{K}(K, D)$.

*Proof* We first show that $\text{SOL}(K, r, D) \neq \emptyset$ for all $r \in \text{int} \, \mathcal{K}(K, D)^*$. By degree theory, it suffices to show that for some vector $u^{\text{ref}} \in K$, the union

$$ \bigcup_{t \in (0,1)} \text{SOL}(K, tr - (1-t)u^{\text{ref}}, tD + (1-t)I), $$
i.e., the set of vectors $u \in K$ such that

$$ (u' - u)^T [ (1-t)(u - u^{\text{ref}}) + t(r + Du) ] \ge 0, \quad \forall u' \in K $$

for some $t \in (0,1)$, is bounded; see, e.g., the proof of [38, Proposition 2.2.3]. Assume for the sake of contradiction that the contrary holds. Hence there exist a sequence of scalars $\{t_k\} \subset (0,1)$ and a sequence of vectors $\{u^k\} \subset K$ such that, for every $k$,

$$ (u' - u^k)^T [ (1 - t_k)(u^k - u^{\text{ref}}) + t_k (r + Du^k) ] \ge 0, \quad \forall u' \in K, $$

and $\{\|u^k\|\} \to \infty$. Without loss of generality, we may assume that

$$ \lim_{k \to \infty} t_k = t_\infty \quad \text{and} \quad \lim_{k \to \infty} \frac{u^k - u^{\text{ref}}}{\|u^k - u^{\text{ref}}\|} = v^\infty, $$

for some scalar $t_\infty \in [0,1]$ and some nonzero vector $v^\infty$, which must belong to $K_\infty$. Letting $u' = 0$, we have, by the positive semidefiniteness of $D$,

$$ \begin{aligned} 0 &\le -(u^k)^T [ (1 - t_k)(u^k - u^{\text{ref}}) + t_k (r + Du^k) ] \\ &\le -(u^k)^T [ (1 - t_k)(u^k - u^{\text{ref}}) + t_k r ]. \end{aligned} $$

Dividing by $\|u^k - u^{\text{ref}}\|^2$ and letting $k \to \infty$, we deduce that $t_\infty = 1$. Let $d$ be an arbitrary vector in $K_\infty$. We have

$$ d^T [ (1 - t_k)(u^k - u^{\text{ref}}) + t_k (r + Du^k) ] \ge 0. $$

Dividing by $\|u^k - u^{\text{ref}}\|$ and letting $k \to \infty$, we deduce that $d^T D v^\infty \ge 0$. Consequently, $D v^\infty \in (K_\infty)^*$. Since

$$ 0 \le (1 - t_k)(u^k)^T u^{\text{ref}} - t_k (u^k)^T (r + Du^k), $$

dividing by $\|u^k - u^{\text{ref}}\|^2$ and letting $k \to \infty$, we deduce $(v^\infty)^T D v^\infty \le 0$. Hence $v^\infty$ is a nonzero solution of the CP (6.7). From the last inequality, we have

$$ 0 \le (1 - t_k)(u^k)^T u^{\text{ref}} - t_k (u^k)^T r. $$

Dividing by $\|u^k - u^{\text{ref}}\|$, letting $k \to \infty$, and using $t_\infty = 1$, we deduce $(v^\infty)^T r \le 0$. Since $v^\infty$ is a nonzero vector in the VI kernel $K(K, D)$ and $r$ is an interior vector in the dual of this kernel, we must have $r^T v^\infty > 0$. We obtain a contradiction. This completes the proof of the nonemptiness of SOL($K, r, D$). The proof of the existence of the scalar $\rho$ is similar to the previous proof. This completes the proof of statement (a). Assume next that $K$ is polyhedral. Write

$$ K \equiv \{ u \in \mathbb{R}^m : Au \le b \}, $$
for some matrix $A$ and vector $b$ of compatible dimensions. Proceeding as above, we may assume for the sake of contradiction that there exist a sequence $\{r^k\} \subset \text{int } K(K, D)^*$ and an unbounded sequence $\{u^k\}$ such that $u^k \in \text{SOL}(K, r^k, D)$ for every $k$ and $\{r^k/\|u^k\|\}$ converges to zero. By the polyhedral representation of $K$, we deduce the existence of a multiplier $\lambda^k$ such that

$$ 
\begin{aligned}
0 &= r^k + Du^k + A^T \lambda^k, \\
0 &\le \lambda^k \perp b - Au^k \ge 0,
\end{aligned}
$$

for every $k$. Without loss of generality, we may assume that

$$ \lim_{k \to \infty} \|u^k\| = \infty, \quad \lim_{k \to \infty} \frac{u^k}{\|u^k\|} = v^\infty, $$

for some nonzero vector $v^\infty$ which must belongs to $K_\infty = \{v \in \mathbb{R}^n : Av \le 0\}$. For each $j$ such that $(Av^\infty)_j < 0$, we must have $(Au^k - b)_j < 0$ for all $k$ sufficiently large, which implies $\lambda_j^k = 0$ for all such $k$. Let $\mathcal{I} = \{i : (Av^\infty)_i = 0\}$. Since

$$ 0 = \frac{r^k}{\|u^k\|} + D \frac{u^k}{\|u^k\|} + \sum_{i \in \mathcal{I}} \tilde{\lambda}_i^k (A_{i\cdot})^T $$

has a solution $\tilde{\lambda}_i^k = \lambda_i^k / \|u^k\| \ge 0$ for all $k$, it follows that for each $k$ sufficiently large, the above equation has a basic feasible solution $\{\tilde{\lambda}_i^k\}_{i \in \mathcal{I}}$ with the property that the vectors $\{(A_{i\cdot})^T\}$ corresponding to the positive components of $\tilde{\lambda}_i^k$ ($i \in \mathcal{I}$) are linearly independent. Moreover, we may choose each $\tilde{\lambda}_i^k$ such that $\tilde{\lambda}_i^k > 0 \Rightarrow (Au^k - b)_i = 0$. Hence there exist an infinite subset $\kappa$ of $\{1, 2, \dots\}$ and an index set $\mathcal{I}_\infty \subseteq \mathcal{I}$ such that $\mathcal{I}_\infty = \{i : \tilde{\lambda}_i^k > 0\}$ for all $k \in \kappa$. Since the vectors in $\{(A_{i\cdot})^T\}_{i \in \mathcal{I}_\infty}$ are linearly independent, it follows that the sequence $\{\tilde{\lambda}_{\mathcal{I}}^k\}_{k \in \kappa}$ is bounded. Without loss of generality, assume that it converges to a vector $\tilde{\lambda}_{\mathcal{I}}^\infty$, which must be nonnegative and satisfies

$$ 0 = Dv^\infty + \sum_{i \in \mathcal{I}} \tilde{\lambda}_i^\infty (A_{i\cdot})^T. $$

Furthermore, if $\tilde{\lambda}_i^\infty > 0$, then $\tilde{\lambda}_i^k > 0$, hence $(Au^k - b)_i = 0$, for all $k \in \kappa$ sufficiently large. Therefore, letting $\tilde{\lambda}_i^\infty = 0 \equiv \tilde{\lambda}_i^k$ for all $i \notin \mathcal{I}$ and for all $k \in \kappa$, we have $0 = Dv^\infty + A^T \tilde{\lambda}$ and $(\tilde{\lambda}^\infty)^T (Au^k - b) = 0$ for all $k \in \kappa$ sufficiently large. Since $0 \in K$, it follows that $b \ge 0$. Hence $(\tilde{\lambda}^\infty)^T A u^k \ge 0$ for all $k \in \kappa$ sufficiently large, which implies $(v^\infty)^T D v^\infty = 0$; thus $v^\infty \in K(K, D)$. Since $D$ is positive semidefinite, it follows that $(D+D^T)v^\infty = 0$; Hence, $A^T \tilde{\lambda}^\infty = -Dv^\infty = D^T v^\infty$,
which yields $(Du^k)^T v^\infty \ge 0$ for all $k \in \kappa$ sufficiently large. Consequently,

$$ 
\begin{aligned}
0 &= (v^\infty)^T \left( \frac{r^k}{\|u^k\|} + D \frac{u^k}{\|u^k\|} + A^T \hat{\lambda}^k \right) \\
&= (v^\infty)^T \left( \frac{r^k}{\|u^k\|} + D \frac{u^k}{\|u^k\|} \right) \ge (v^\infty)^T \frac{r^k}{\|u^k\|}.
\end{aligned}
$$

This contradicts the membership of $r^k$ in int $K(K, D)^*$. Finally, since the AVI range of an affine pair is a closed set [38, part (b) of Theorem 2.5.15], it follows that $\text{SOL}(K, r, D) \ne \emptyset$ for all $r \in K(K, D)^*$. Noticing that $K(K, D)^*$ is a polyhedron set, thus convex, we deduce the existence of the Lipschitz constant $L_D$ from the aforementioned facts. □

We illustrate part (b) of Proposition 6.4 with the LCP arising from the ODE (4.1).

*Example 6.1* Let $K = \mathbb{R}_+^{2n}$,

$$ r = \begin{pmatrix} x \\ 2 \end{pmatrix} \quad \text{and} \quad D = \begin{bmatrix} 0 & I \\ -I & 0 \end{bmatrix}. $$

The LCP kernel of $D$ is $\{0\} \times \mathbb{R}_+^n$; the dual of this kernel is therefore $\mathbb{R}^n \times \mathbb{R}_+^n$. Clearly, the vector $r$ belongs to the interior of this dual for any $x \in \mathbb{R}^n$. Therefore, part (b) of Proposition 6.4 applies. The reader can easily verify that

$$ \sup\{\|u\| : u \in \text{SOL}(r, D)\} \le 2(1 + \|r\|), \quad \forall r \in \mathbb{R}^n. $$

Since $D + D^T = 0$, the Lipschitz property of $(D + D^T)u$ is obvious.

When $K$ is a polyhedron, we can introduce nonlinear perturbations of $D$ in Proposition 6.4 by restricting $D$ to be psd-plus [38,65]; i.e., $D$ is positive semidefinite and

$$ [u^T Du = 0 \Rightarrow Du = 0]. $$

This is equivalent to saying the linear map $u \mapsto Du$ is monotone-plus. The following known result [65] shows among other things that for a linear map, the strongly monotone (linear) composite property is equivalent to the monotone-plus property.

**Lemma 6.4** Let $M \in \mathbb{R}^{m \times m}$. The following three statements are equivalent:

(a) There exists $c > 0$ such that $u^T M u \ge c \|M u\|^2$ for all $u$ in $\mathbb{R}^m$.

(b) $M = E^T A E$ for some matrices $E \in \mathbb{R}^{l \times m}$ and $A \in \mathbb{R}^{l \times l}$ for some positive integer $l$, where $A$ is positive definite.

(c) $M$ is psd-plus.
If $D$ is a psd-plus matrix, then $Du$ is a constant for all $u \in \text{SOL}(K, r, D)$ [38, Corollary 2.3.10]. For ease of reference, we summarize in the lemma below further properties of an AVI with such a matrix $D$.

**Lemma 6.5** Let $K \subseteq \mathbb{R}^m$ be a polyhedron containing either the origin or no lines. Let $D \in \mathbb{R}^{m \times m}$ be a psd-plus matrix such that $(K, D)$ is an $R_0$ pair. The AVI $(K, r, D)$ has a solution for all $r \in \mathbb{R}^m$. Moreover, there is a constant $L_V > 0$, such that

$$ \| Du^1 - Du^2 \| \le L_V \| r^1 - r^2 \|, \qquad (6.9) $$

for all $r^i$ in $\mathbb{R}^m$ and $u^i \in \text{SOL}(K, r^i, D)$, $i = 1, 2$.

*Remark* If $D$ is a psd-plus matrix, $(K, D)$ is an $R_0$ pair if and only if $K_\infty \cap \ker D = \{0\}$. Hence this assumption on $D$ in Lemma 6.5 is analogous to that on the strongly monotone composite (nonlinear) map $F = E^T \circ \Psi \circ E$ in Proposition 6.3. Yet, for the latter map, while $Eu$ is a constant for $u \in \text{SOL}(K, r+F)$, the (single-valued) map $r \mapsto ESOL(K, r+F)$ is not even known to be locally Lipschitz continuous for a non-polyhedral set $K$ (for a result along these lines, see [38, Corollary 6.2.11]). In other words, the polyhedrality of $K$ plays an essential role in Lemma 6.5, and hence in Proposition 6.5.

The following result is based on the above lemma and a fixed-point argument. A major difference between the result below and the previous three propositions is that while the solution set of the VI $(K, q + F)$ remains convex, the (nonlinear) map $F = D + \Phi$ is not necessarily monotone.

**Proposition 6.5** Let $K \subseteq \mathbb{R}^m$ be a polyhedron containing either the origin or no lines. Let $D \in \mathbb{R}^{m \times m}$ be a psd-plus matrix such that $(K, D)$ is an $R_0$ pair. Let $\rho > 0$ be such that (6.8) holds. Let $\Phi : \mathbb{R}^m \to \mathbb{R}^m$ be a continuous function such that

$$ \| \Phi(u) \| \le L_{\Phi} \| u \|, \quad \forall u \in K, \qquad (6.10) $$

for some constant $L_\Phi \in (0, 1/\rho)$. For every vector $q \in \mathbb{R}^m$, $\text{SOL}(K, q + F) \neq \emptyset$, where $F = D + \Phi$; moreover,

$$ \sup\{ \|u\| : u \in \text{SOL}(K, q + F) \} \le \frac{\rho (1 + \|q\|)}{1 - \rho L_{\Phi}}. \qquad (6.11) $$

Assume in addition

$$ \| \Phi(u) - \Phi(u') \| \le L'_{\Phi} \| Du - Du' \|, \quad \forall u, u' \in K, \qquad (6.12) $$

for some constant $L'_\Phi \in (0, 1/L_V)$, where $L_V > 0$ is such that (6.9) holds for all $r^i$ in $\mathbb{R}^m$ and $u^i \in \text{SOL}(K, r^i, D)$, $i = 1, 2$. It holds that

$$ \| Du^1 - Du^2 \| \le \frac{L_V \| q^1 - q^2 \|}{1 - L_V L'_{\Phi}}, \qquad (6.13) $$
for all $q^i \in \mathbb{R}^m$ and $u^i \in \text{SOL}(K, q^i + F)$, $i = 1, 2$. Finally, for every $q \in \mathbb{R}^m$,

$$ \text{SOL}(K, q + F) = D^{-1}(v(q)) \cap \operatorname*{argmin}\{u^T w(q) : u \in K\}, \quad (6.14) $$

where $v(q) = D\tilde{u}$ and $w(q) = q + F(\tilde{u})$ for any $\tilde{u} \in \text{SOL}(K, q+F)$, and $D^{-1}(v(q))$ is the inverse image of $v(q)$ under the (possibly singular) matrix $D$. Consequently, $\text{SOL}(K, q + F)$ is a convex set.

*Proof* Choose an arbitrary scalar $R > \frac{\rho(1 + \|q\|)}{1 - \rho L_\phi}$ such that the intersection, denoted $\mathcal{U} = K \cap R\mathcal{B}$, of $K$ with the Euclidean ball $R\mathcal{B}$ in $\mathbb{R}^m$ centered at the origin and with radius $R$ is nonempty. We define a set-valued mapping $\Xi$ from $\mathcal{U}$ into subsets of $\mathcal{U}$ as follows. For $u \in \mathcal{U}$, let $\Xi(u)$ be the solution set of the VI $(K, r(u), D)$, where $r(u) = q + \Phi(u)$. Clearly $\Xi(u)$ is a nonempty closed convex subset of $K$. We need to show that $\|v\| \le R$ for all $v \in \Xi(u)$. Indeed, we have

$$ \|v\| \le \rho(1 + \|r(u)\|) \le \rho(1 + \|q\| + L_\phi R) < R. $$

Since $\Phi$ is continuous, it is easy to see that $\Xi$ is a closed (set-valued) map. Consequently, by Kakutani's fixed-point theorem, there exists $u \in \Xi(u)$. This vector $u$ clearly is an element of $\text{SOL}(K, q + F)$. Moreover, the bound (6.11) clearly follows from the above bound on $\|v\|$.

Notice that if $u \in \text{SOL}(K, q + F)$, then $u \in \text{SOL}(K, q + \Phi(u), D)$. Consequently, if $u^i \in \text{SOL}(K, q^i + F)$ for $i = 1, 2$, then we have

$$ \| Du^1 - Du^2 \| \le L_v \| q^1 + \Phi(u^1) - q^2 - \Phi(u^2) \|. $$

In view of (6.12), the bound (6.13) follows readily. Finally, to establish (6.14), we first note that, by (6.13), $Du$ is a constant vector for all $u \in \text{SOL}(K, q + F)$; hence, by (6.12), so is $\Phi(u)$, thus so is $F(u)$. Consequently, the inclusion

$$ \text{SOL}(K, q + F) \subseteq D^{-1}(v(q)) \cap \operatorname*{argmin}\{u^T w(q) : u \in K\} $$

is obvious from the constancy of $v(q)$ and $w(q)$. Conversely, if $u'$ is an element of the right-hand set, then $Du' = v(q) = D\tilde{u}$; thus $\Phi(u') = \Phi(\tilde{u})$. Hence $w(q) = q + F(\tilde{u}) = q + F(u')$. Consequently, for all $u \in K$,

$$ 0 \le (u - u')^T w(q) = (u - u')^T (q + F(u')), $$

establishing the equality (6.14). This representation of $\text{SOL}(K, q + F)$ clearly shows that this solution set is convex. $\square$

Although the VI $(K, q + F)$ in Proposition 6.5 has a convex solution set, the map $F = D + \Phi$ is not necessarily monotone. As an example, consider the linear
complementarity problem where $K$ is the nonnegative orthant in the plane. Let $D$ be the $2 \times 2$ matrix of all ones. Let $\Phi(u) \equiv Eu$, where

$$E = \begin{bmatrix} \varepsilon & \varepsilon \\ 0 & 0 \end{bmatrix}.$$

For any $\varepsilon > 0$ sufficiently small, the matrix

$$D + E = \begin{bmatrix} 1 + \varepsilon & 1 + \varepsilon \\ 1 & 1 \end{bmatrix}$$

is column sufficient [30, Definition 3.5.1] but not positive semidefinite. Yet with such an $\varepsilon$, all other conditions in Proposition 6.5 are valid for the pair $(D, \Phi)$.

Summarizing the discussion in this section, the theorem below is an existence result for the initial-value DVI under various independent assumptions on the pair $(K, F)$. In particular, each of the five cases in the theorem pertains to a different class of DVIs. For instance, for an affine pair $(K, D)$, condition (d) distinguishes itself in that the map $G$ is linked to $F$ via the interiority condition: $G(\Omega) \subseteq \text{int } \mathcal{K}(K, D)^*$; while the latter inclusion is trivially valid when $(K, D)$ is an $\mathbb{R}_0$ pair, condition (c) allows $K$ to be non-polyhedral.

**Theorem 6.1** Let $K \subseteq \mathbb{R}^m$ be a nonempty closed convex set. Let $(f, B, G)$ satisfy conditions (A) and (B). Let $F : \mathbb{R}^m \to \mathbb{R}^m$ be given. The initial-value DVI (6.2) has a weak solution under any one of the following five conditions:

(a) $F$ is continuous and monotone, and there exists $u^{\text{ref}} \in K$ satisfying (6.6);

(b) $F = E^T \circ \Psi \circ E$ for some matrix $E \in \mathbb{R}^{l \times m}$ satisfying $K_\infty \cap \ker E = \{0\}$ and for some continuous map $\Psi : \mathbb{R}^l \to \mathbb{R}^l$ that is strongly monotone on $EK$;

(c) $K$ contains the origin, $F(u) = Du$ for some positive semidefinite matrix $D$ such that $(K, D)$ is an $\mathbb{R}_0$ pair;

(d) $K$ is a polyhedron containing the origin, $F(u) = Du$ for some positive semidefinite matrix $D$, and $G(\Omega) \subseteq \text{int } \mathcal{K}(K, D)^*$;

(e) $K$ is a polyhedron containing either the origin or no lines, and $F = D + \Phi$, where $D$ is a psd-plus matrix such that $(K, D)$ is an $\mathbb{R}_0$ pair, and where $\Phi : \mathbb{R}^m \to \mathbb{R}^m$ is continuous and satisfies (6.10) and (6.12).

*Proof* All parts follow from Lemmas 6.1 and 6.2. The following proof details how the conditions of these lemmas are satisfied by virtue of the results in this subsection.

(a) As $F$ is continuous and (6.6) holds, by Proposition 6.2, the set $F(t, x)$ is nonempty and closed, and $F$ has the linear growth property (6.5); since $F$ is monotone the set $F(t, x)$ is convex.

(b) By Proposition 6.3 the set $F(x, t)$ is nonempty and convex, and furthermore satisfies the linear growth property (6.5).
(c) By Proposition 6.4, $\mathbf{F}(x,t)$ is nonempty and convex, and by Proposition 6.4(a), $\mathbf{F}$ satisfies the linear growth property (6.5).

(d) By Proposition 6.4, $\mathbf{F}(x,t)$ is nonempty and convex, and by Proposition 6.4(b), $\mathbf{F}$ satisfies the linear growth property (6.5).

(e) By Proposition 6.5, and equation (6.14), $\mathbf{F}(x,t)$ is non-empty and convex; by (6.11), $\mathbf{F}$ satisfies the linear growth property (6.5). $\square$

Continuing Example 6.1, we deduce that part (d) of Theorem 6.1 is applicable to the ODE (4.1). More broadly, the theorem is applicable to a generalized LCS:

$$ \begin{gathered} \dot{x} = p + Ax + Bu, \tag{6.15} \\ C \ni u \perp q + Cx + Du \in C^*, \end{gathered} $$

where $C$ is a polyhedral cone, $D$ is positive semidefinite, and $q + C\mathbb{R}^n \subseteq \text{int}(K(C, D)^*)$. This corollary is related to similar results for the LCS (2.8) with $p=0$ and $q=0$ established in [21,49], where it is further assumed that the tuple $(A, B, C, D)$ is “minimal” and “passive”. (See also [44] which contains further results on passivity and the LCS.) There is one important difference, however. Our result shows that if $q + C\mathbb{R}^n \subseteq \text{int}(K(C, D)^*)$, then (6.15) has a solution for all initial conditions $x^0 \in \mathbb{R}^n$; in contrast, a result in the cited references shows that, in the case where $p=0$ and $q=0$ and $C = \mathbb{R}_+^m$, the passive system (6.15) has a solution for a given $x^0$, provided that $Cx^0$ belongs to the dual of the LCP kernel of $D$.

## 7 Convergence of a time-stepping scheme: IVPs

We remind the reader that conditions (A) and (B) introduced in Sect. 6.1 will remain in force throughout this and the next two sections. At this time, it is appropriate to mention the following weaker assumption:

(A') $B$ and $G$ are satisfy the condition in (A), but $f(t, \cdot)$ is Lipschitz with Lipschitz constant $L_f$ for all $t$, $f(\cdot, x)$ is measurable for all $x \in \mathbb{R}^n$, and $\|f(t, 0)\| \le \rho_{f,0}$ for all $t \in [0, T]$.

Under Assumption (A') the theorems and lemmas in the remainder of the paper can be shown to hold. Nevertheless, suitable modifications and additional technical arguments are needed. For instance, the discrete-time equations such as (7.2) should have $h f(t_{h,i}, z_{h,i})$ and $h f(t_{h,i+1}, z_{h,i})$ replaced by $\int_{t_{h,i}}^{t_{h,i+1}} f(\tau, z_{h,i}) d\tau$. Assumption (A') is not the most general setting in which existence results for DVIs can be shown, but it suffices for most applications.

We consider a special time-stepping scheme tailored to the initial-value DVI (6.2). The scheme is as follows; cf. (3.3). With $x^{h,0} \equiv x^0$ given, we compute

$$ \{x^{h,1}, x^{h,2}, \dots, x^{h,N_h+1}\} \subset \mathbb{R}^n \quad \text{and} \quad \{u^{h,1}, u^{h,2}, \dots, u^{h,N_h+1}\} \subset K \quad (7.1) $$
by the recursion: for $i = 0, 1, \dots, N_h$, where $N_h \equiv T/h - 1$,

$$x^{h,i+1} = x^{h,i} + h[f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1}) + B(t_{h,i}, x^{h,i})u^{h,i+1}], \quad (7.2)$$
$$u^{h,i+1} \in \text{SOL}(K, G(t_{h,i+1}, x^{h,i+1}) + F).$$

Despite the similarity with some schemes developed for mechanical systems subject to unilateral constraints and Coulomb friction (see e.g. [5,57,73,91, 94,98]), the convergence of (7.2) has not been analyzed in its generality. A distinguished feature of the iteration (7.2) is that at iteration $i+1$, the arguments in $B(t_{h,i}, x^{h,i})$ are set at the previous time $t_{h,i}$. This “explicit” use of the pair $(t_{h,i}, x^{h,i})$ in the ODE allows us to conveniently solve for $x^{h,i+1}$ in terms of $u^{h,i+1}$, albeit only implicitly because the unknown $x^{h,i+1}$ is contained in the nonlinear term $f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1})$. There are many variations of the iteration (7.2), some of which may be particularly suitable for problems with special structures. For instance, for the DVI (2.10) derived from the KKT form of the VIE with a convex finitely represented set, one can consider the following fully implicit iterative scheme: (the function $F$ below refers to (2.10) and should not be confused with that in (7.2))

$$x^{h,i+1} - x^{h,i} = -h F(x^{h,i+1}) - \sum_{j=1}^{m} v_j^{h,i+1} \nabla g_i(x^{h,i+1}), \quad (7.3)$$
$$0 \le v^{h,i+1} \perp -g(x^{h,i+1}) \ge 0,$$

which can easily be seen to be the KKT conditions of the finite dimensional VI $(K, F^{h,i+1})$, where $F^{h,i+1}(x) \equiv x - x^{h,i} - hF(x)$ is strongly monotone with a modulus that is independent of $h$, provided that $h > 0$ is sufficiently small and $F$ is Lipschitz, and where $v^{h,i+1} = hu^{h,i+1}$ is the KKT multiplier of the constraint $g(x^{h,i+1}) \le 0$. As such, the existence of the pair $(x^{h,i+1}, v^{h,i+1})$ is easy to establish based on finite-dimensional VI theory.

The practical implementation of the iterative scheme (7.2) depends on the solution of finite-dimensional VIs. Relying on [38, Volume II] for a comprehensive summary of applicable methods and their convergence theory and on the PATH solver [33,40] that contains ideas of many such methods and which is publicly available on the NEOS server (http://neos.mcs.anl.gov/neos/solvers), we omit the details of the computational aspects and focus instead on the analysis of the iteration (7.2). This line of analysis can easily be modified to suit the variations of the basic scheme, such as (7.3). To begin, we establish a preliminary lemma for the time-discretized equation in (7.2) under assumptions (A) and (B).

**Lemma 7.1** Under (A) and (B), $h_0 > 0$ exists such that for all $h \in (0, h_0]$, $(x^{\text{ref}}, u) \in \mathbb{R}^{n+m}$, $\theta \in [0, 1]$, and $t$ and $t_{\text{ref}}$ in $[0, T]$, there exists a unique vector $x(u)$ satisfying

$$x(u) - x^{\text{ref}} = h[f(t, \theta x^{\text{ref}} + (1-\theta)x(u)) + B(t_{\text{ref}}, x^{\text{ref}})u].$$
Moreover, for all $u, u'$ in $\mathbb{R}^m$,

$$ 
\begin{aligned}
\|x(u) - x(u')\| &\le \frac{h \sigma_B}{1 - h(1 - \theta)L_f} \|u - u'\|, \\
\|x(u) - x^{\text{ref}}\| &\le h \frac{\rho_f (1 + \|x^{\text{ref}}\|) + \sigma_B \|u\|}{1 - h(1 - \theta)\rho_f}.
\end{aligned} 
\quad (7.4) $$

*Proof* It suffices to choose $h_0$ to satisfy

$$ 0 < h_0 < \min\left(\frac{1}{(1-\theta)L_f}, \frac{1}{(1-\theta)\rho_f}\right). $$

(The right-hand side is taken to be $\infty$ if $\theta = 1$.) Under this choice, consider any tuple $(h, x^{\text{ref}}, u, t, t_{\text{ref}})$ as specified. With $f$ being Lipschitz continuous, the map $x \mapsto hf(t, \theta x^{\text{ref}} + (1-\theta)x)$ can readily be seen to be a contraction. Consequently, $x(u)$ exists and is unique. The two inequalities in (7.4) can be established easily. $\square$

The above result does not establish the existence of $(x^{h,i+1}, u^{h,i+1})$ satisfying (7.2). Postponing this issue until later (see Propositions 7.1 and 8.2), we use the iterates (7.1) to construct functions of time by interpolation. Specifically, let $\hat{x}^h(\cdot)$ be the continuous piecewise linear interpolant of the family $\{x^{h,i}\}$; i.e., for $i = 0, 1, \dots, N_h$,

$$ \hat{x}^h(t) = x^{h,i} + \frac{t - t_{h,i}}{h} (x^{h,i+1} - x^{h,i}), \quad \forall t \in [t_{h,i}, t_{h,i+1}]; $$

let $\hat{u}^h(\cdot)$ be the (possibly discontinuous) piecewise constant interpolant of the family $\{u^{h,i}\}$, i.e., $\hat{u}^h(t) = u^{h,i+1}$ for $t \in (t_{h,i}, t_{h,i+1}]$. In the convergence analysis, we consider two main classes of functions $F$, one of which consists of composite maps of the form $F \equiv \Psi \circ E$, where $E \in \mathbb{R}^{l \times m}$ and $\Psi : \mathbb{R}^l \to \mathbb{R}^m$. We let $\hat{w}^h(\cdot)$ be the continuous piecewise linear interpolant of the family $\{Eu^{h,i}\}$, where $u^{h,0} = u^{h,1}$, or alternatively, $u^{h,0}$ can be taken to be an element of $\text{SOL}(K, G(0, x^0) + F)$, provided that the latter VI is solvable. We have, for $i = 0, 1, \dots, N_h$,

$$ \hat{w}^h(t) = Eu^{h,i} + \frac{t - t_{h,i}}{h} E(u^{h,i+1} - u^{h,i}), \quad \forall t \in [t_{h,i}, t_{h,i+1}]. $$

Consequently, for every $h > 0$ sufficiently small, the functions $\hat{x}^h : [0, T] \to \mathbb{R}^n$, $\hat{u}^h : [0, T] \to \mathbb{R}^m$, and $\hat{w}^h : [0, T] \to \mathbb{R}^l$ are well defined. We state and prove in the theorem below the main convergence of these functions as $h$ tends to zero; the theorem postulates some bounds on the iterates $\{x^{h,i+1}\}$ and $\{u^{h,i+1}\}$, which will be verified to hold under appropriate conditions on the data functions $(f, B, G, F)$ and the set $K$.
**Theorem 7.1** Let $K \subseteq \mathbb{R}^m$ be a nonempty closed convex set and let $(f, B, G)$ satisfy conditions (A) and (B). Suppose that there exist positive scalars $c_{0,x}, c_{1,x}, c_{0,u}, c_{1,u}$, and $\bar{h}$ such that for all $h \in (0, \bar{h}]$ and all $i = 0, 1, \dots, N_h$,

$$ \| x^{h,i+1} \| \le c_{0,x} + c_{1,x} \| x^0 \| \quad \text{and} \quad \| u^{h,i+1} \| \le c_{0,u} + c_{1,u} \| x^0 \|. \quad (7.5) $$

There is a sequence $\{h_v\} \downarrow 0$ such that the following two limits exist: $\hat{x}^{h_v} \to \hat{x}$ uniformly on $[0, T]$ and $\hat{u}^{h_v} \to \hat{u}$ weakly in $L^2(0, T)$. Furthermore, under either one of the following two conditions:

(a) $F(u) = \Psi(Eu)$, where $E \in \mathbb{R}^{\ell \times m}$ and $\Psi: \mathbb{R}^{\ell} \to \mathbb{R}^m$ is Lipschitz continuous, and a constant $c_{2,u} > 0$ exists such that for all $h$ sufficiently small,

$$ \| Eu^{h,i+1} - Eu^{h,i} \| \le h c_{2,u}, \quad (7.6) $$

(b) $F(u) = Du$ for some positive semidefinite matrix $D$, all such limits $(\hat{x}, \hat{u})$ are weak solutions of the initial-value DVI (6.2).

*Proof* Throughout the proof below, $h \in (0, h_0]$ is taken to be sufficiently small. By the second inequality in (7.4), we have

$$ \| x^{h,i+1} - x^{h,i} \| \le h \frac{\rho_f (1 + \| x^{h,i} \|) + \sigma_B \| u^{h,i+1} \|}{1 - h(1 - \theta) \rho_f}. $$

Consequently, for an appropriate constant $L_{x^0} > 0$ that is independent of $h$ (but depends on $\|x^0\|$), we deduce that for all $h > 0$ sufficiently small and for all $i = 0, 1, \dots, N_h$, $\|x^{h,i+1} - x^{h,i}\| \le L_{x^0}h$. This implies that the piecewise interpolants $\hat{x}^h$ are not only Lipschitz continuous on $[0, T]$, but the Lipschitz constant is independent of $h$. Hence there is a positive scalar $h'_0$, which depends only on the data of the problem, such that the family of functions $\{\hat{x}^h\}$ for $h$ in $(0, h'_0]$ is an equicontinuous family of functions. Furthermore, since $x^0$ is fixed, we have a uniform bound on

$$ \| \hat{x}^h \|_{L^\infty} \equiv \sup_{t \in [0,T]} \| \hat{x}^h(t) \|. $$

So by the Arzelá–Ascoli theorem (see, e.g., [99, p. 167] or [63, pp. 57–59]), there is a sequence $\{h_v\} \downarrow 0$ such that $\{\hat{x}^{h_v}\}$ converges in the supremum (i.e., $L^\infty$) norm to a Lipschitz function $\hat{x}$ on $[0, T]$. From the uniform boundedness of the iterates $u^{h,i+1}$, we can conclude that the family of piecewise constant interpolants $\{\hat{u}^h\}$ are uniformly bounded in the $L^\infty$ norm on $(0, T)$; that is, for some constant $\psi_u > 0$, $\|\hat{u}^h\|_{L^\infty} \le \psi_u$ for all $h > 0$ sufficiently small. By Alaoglu’s theorem [63, pp. 71–72], and by working with an appropriate subsequence of $\{h_v\}$ if necessary, we may assume without loss of generality that the sequence $\{\hat{u}^{h_v}\}$ has a weak* limit $\hat{u}$. Since $L^2(0, T)$ is a reflexive Banach space, as it is a Hilbert space, weak*
convergent sequences are also weakly convergent sequences. Consequently, for any function $\varphi \in L^2(0, T)$,

$$ \lim_{\nu \to \infty} \int_0^T \varphi(t)^T \hat{u}^{h_\nu}(t) \, dt = \int_0^T \varphi(t)^T \hat{u}(t) \, dt. $$

Mazur's theorem [63, p. 88] implies that there is a sequence of convex combinations of $\{\hat{u}^{h_\nu}\}$ that converges strongly in $L^2(0, T)$. Since $\hat{u}^h(t) \in K$ for all $t$, and $K$ is a closed convex set, any convex combination of $\{\hat{u}^h(t)\}$ will also be in $K$. Thus, $\hat{u}$ is a strong limit in $L^2(0, T)$ of functions with values in $K$. Since strongly convergent sequences in $L^p$ spaces ($p < \infty$) have subsequences that are pointwise convergent almost everywhere [63, p. 210], it follows that $\hat{u}(t) \in K$ for almost all $t$.

Suppose first that condition (a) holds. Since $\|\hat{w}^h\|_{L^\infty} \le \|E\| \|\hat{u}^h\|_{L^\infty}$, it follows that the family $\{\hat{w}^h\}$ is uniformly bounded in $L^\infty(0, T)$ as $h \downarrow 0$. By (7.6), we may demonstrate as we did for the family $\{\hat{x}^h\}$ that the sequence $\{\hat{w}^{h_\nu}\}$ contains a subsequence converging uniformly to a limit function $\hat{w}$ which is also Lipschitz continuous on $[0, T]$. Without loss of generality, we may take $\{\hat{w}^{h_\nu}\}$ itself to be such a subsequence. Since the sequence $\{\hat{u}^{h_\nu}\}$ converges weakly in $L^2(0, T)$ to $\hat{u}$, it follows that $\{E\hat{u}^{h_\nu}\}$ also converges weakly to $E\hat{u}$ in $L^2(0, T)$. Since $\hat{w}^h(t_j) = E\hat{u}^h(t_j)$ for all $t_j$, and since $E\hat{u}^h$ is Lipschitz continuous with a constant independent of $h$, a constant $L_w > 0$ exists such that $\|\hat{w}^h(t) - E\hat{u}^h(t)\| \le L_w h$ on $[0, T]$, for all $h > 0$ sufficiently small. Thus $\hat{w}^{h_\nu} - E\hat{u}^{h_\nu}$ converges uniformly to zero as $\nu \to \infty$. Since $\hat{w}^{h_\nu}$ converges uniformly to $\hat{w}$, we see that $E\hat{u}^{h_\nu}$ converges uniformly to $\hat{w}$, which is also $E\hat{u}$. We need to show that the limit functions $(\hat{x}, \hat{u})$ indeed solve the DVI. We have already shown that $\hat{u}(t) \in K$ for almost all $t$. The remaining three things to be shown are (i) for all continuous functions $\tilde{u} : [0, T] \to K$,

$$ \int_0^T (\tilde{u}(t) - \hat{u}(t))^T [G(t, \hat{x}(t)) + \Psi(E\hat{u}(t))] dt \ge 0, \quad (7.7) $$

(ii) for any $0 \le s \le t \le T$,

$$ \hat{x}(t) - \hat{x}(s) = \int_s^t [f(\tau, \hat{x}(\tau)) + B(\tau, \hat{x}(\tau))\hat{u}(\tau)] d\tau, \quad (7.8) $$

and (iii) the initial condition $\hat{x}(0) = x^0$. The last condition (iii) holds since $\hat{x}^h(0) = x^0$ for all $h > 0$ sufficiently small and $\hat{x}^{h_\nu} \to \hat{x}$ uniformly as $\nu \to \infty$. To show (i), note that the sequence $\{G(t, \hat{x}^{h_\nu}) + \Psi(E\hat{u}^{h_\nu})\}$ converges to
$G(t, \hat{x}) + \Psi(E\hat{u})$ uniformly. This implies that

$$ 
\begin{aligned}
& \lim_{v \to \infty} \int_0^T (\tilde{u}(t) - \hat{u}^{h_v}(t))^T [G(t, \hat{x}^{h_v}(t)) + \Psi(E\hat{u}^{h_v}(t))] dt \\
&= \int_0^T (\tilde{u} - \hat{u})^T [G(t, \hat{x}) + \Psi(E\hat{u})] dt.
\end{aligned}
$$

On the other hand,

$$ 
\begin{aligned}
& \int_0^T (\tilde{u} - \hat{u}^h)^T [G(t, \hat{x}^h) + \Psi(E\hat{u}^h)] dt \\
&= \sum_{i=0}^{N_h} \int_{t_{h,i}}^{t_{h,i+1}} (\tilde{u}(t) - u^{h,i+1})^T [G(t, \hat{x}^h(t)) + \Psi(Eu^{h,i+1})] dt \\
&\geq h \sum_{i=0}^{N_h} \left( \frac{1}{h} \int_{t_{h,i}}^{t_{h,i+1}} \tilde{u}(t) dt - u^{h,i+1} \right)^T [G(t_{h,i+1}, x^{h,i+1}) + \Psi(Eu^{h,i+1})] \\
&\quad - h^2 (N_h + 1) (\|\tilde{u}\|_{L^\infty} + \|\hat{u}^h\|_{L^\infty}) L_G L_x \\
&\geq -h T (\|\tilde{u}\|_{L^\infty} + \|\hat{u}^h\|_{L^\infty}) L_G L_x,
\end{aligned}
$$

where the last inequality holds because $\frac{1}{h} \int_{t_i}^{t_{i+1}} \tilde{u}(t)dt$ belongs to $K$ by convexity of $K$ and also because $u^{h,i+1} \in \text{SOL}(K, G(t_{h,i+1}, x^{h,i+1}) + F)$. Hence,

$$ \int_0^T (\tilde{u} - \hat{u})^T [G(t, \hat{x}) + \Psi(E\hat{u})] dt \geq 0 $$

for any continuous $\tilde{u} : [0, T] \to K$. To complete the convergence proof under (a), it remains to show (7.8). Similar to the above derivation, we have

$$ 
\begin{aligned}
x^{h,i+1} - x^{h,i} &= h[f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1}) + B(t_{h,i}, x^{h,i})u^{h,i+1}] \\
&= \int_{t_{h,i}}^{t_{h,i+1}} [f(\tau, \hat{x}^h(\tau)) + B(\tau, \hat{x}^h(\tau)) u^{h,i+1}] d\tau + h^2 \zeta^i,
\end{aligned}
$$
where $\|\zeta^i\| \le (L_f + L_f L_x + \sigma_B \psi_u)$. It follows that for any $0 \le s \le t \le T$,

$$ \hat{x}^h(t) - \hat{x}^h(s) = \int_s^t [f(\tau, \hat{x}^h(\tau)) + B(\tau, \hat{x}^h(\tau))\hat{u}^h(\tau)] d\tau + O(h). $$

By continuity of $f$ and uniform convergence of $\hat{x}^{h_v} \to \hat{x}$, we deduce

$$ \lim_{v \to \infty} \int_s^t f(\tau, \hat{x}^{h_v}(\tau)) d\tau = \int_s^t f(\tau, \hat{x}(\tau)) d\tau. $$

Furthermore,

$$ \left\| \int_s^t [ B(\tau, \hat{x}^h(\tau)) \hat{u}^h(\tau) - B(\tau, \hat{x}(\tau)) \hat{u}(\tau) ] d\tau \right\| \\ \le \int_s^t \| B(\tau, \hat{x}^h(\tau)) - B(\tau, \hat{x}(\tau)) \| \| \hat{u}^h(\tau) \| d\tau \\ + \left\| \int_s^t [ B(\tau, \hat{x}(\tau)) \hat{u}^h(\tau) - B(\tau, \hat{x}(\tau)) \hat{u}(\tau) ] d\tau \right\|. $$

Restricted to the subsequence $\{h_v\}$, the first summand on the right-hand side converges to zero because $\{\hat{x}^{h_v}\} \to \hat{x}$ uniformly and $B$ is continuous; the second summand converges to zero because $\{\hat{u}^{h_v}\}$ converges weakly in $L^2(0, T)$ to $\hat{u}$ and $B(\tau, \hat{x}(\tau))$ is continuous in $\tau$. Consequently, (7.8) holds.

Next assume condition (b). It suffices to show that for all continuous $\tilde{u}: [0, T] \to K$,

$$ \int_0^T (\tilde{u}(t) - \hat{u}(t))^T [G(t, \hat{x}(t)) + D\hat{u}(t)] dt \ge 0, $$

In turn, it suffices to show that

$$ \limsup_{v \to \infty} \int_0^T (\tilde{u}(t) - \hat{u}^{h_v}(t))^T [G(t, \hat{x}^{h_v}(t)) + D\hat{u}^{h_v}(t)] dt \\ \le \int_0^T (\tilde{u}(t) - \hat{u}(t))^T [G(t, \hat{x}(t)) + D\hat{u}(t)] dt $$
because the left-hand limit is nonnegative as shown above. Since

$$ \lim_{v \to \infty} \int_0^T \tilde{u}(t)^T [G(t, \hat{x}^{h_v}(t)) + D\tilde{u}^{h_v}(t)] dt = \int_0^T \tilde{u}(t)^T [G(t, \hat{x}(t)) + D\hat{u}(t)] dt $$

and

$$ \lim_{v \to \infty} \int_0^T (\hat{u}^{h_v}(t))^T G(t, \hat{x}^{h_v}(t)) dt = \int_0^T (\tilde{u}(t))^T G(t, \hat{x}(t)) dt, $$

it remains to show

$$ \liminf_{v \to \infty} \int_0^T (\hat{u}^{h_v}(t))^T D\hat{u}^{h_v}(t) dt \ge \int_0^T (\tilde{u}(t))^T D\tilde{u}(t) dt. \quad (7.9) $$

We can write

$$ \begin{aligned} (\hat{u}^{h_v}(t))^T D\hat{u}^{h_v}(t) &= (\hat{u}(t))^T D\hat{u}(t) + (\hat{u}^{h_v} - \hat{u}(t))^T (D + D^T)\hat{u}(t) \\ &\quad + (\hat{u}^{h_v}(t) - \hat{u}(t))^T D(\hat{u}^{h_v} - \hat{u}) \\ &\ge (\hat{u}(t))^T D\hat{u}(t) + (\hat{u}^{h_v}(t) - \hat{u}(t))^T (D + D^T)\hat{u}(t). \end{aligned} $$

Since $\hat{u}^{h_v} \to \hat{u}$ weakly in $L^2(0, T)$ and $(D + D^T)\hat{u} \in L^2(0, T)$, it follows that

$$ \lim_{v \to \infty} \int_0^T (\hat{u}^{h_v}(t) - \hat{u}(t))^T (D + D^T) \hat{u}(t) dt = 0, $$

and we conclude that (7.9) holds. □

Our next order of business is to identify classes of DVIs for which the desired iterates (7.1) exist and satisfy the assumptions in Theorem 7.1. By way of preparation, we state and prove a Gronwall-type lemma, whose conclusions are the main conditions for Theorem 7.1.

**Lemma 7.2** Suppose there are positive constants $h_1$, $\rho_u$, and $\psi_x$ such that for all $h \in (0, h_1)$ and all nonnegative integers $i$ with $(i+1)h \le T$,

$$ \| u^{h, i+1} \| \le \rho_u (1 + 2 \| x^{h, i} \|), \quad (7.10) $$

$$ \| x^{h, i+1} - x^{h, i} \| \le h \psi_x (1 + \| x^{h, i} \|). \quad (7.11) $$

Then there are constants $c_{0,x}$, $c_{1,x}$, $c_{0,u}$, and $c_{1,u}$ such that (7.5) holds for all $h \in (0, h_1]$ and all $i = 0, 1, \dots, N_h \equiv \lfloor T/h \rfloor - 1$.
*Proof* For all $i = 0, 1, \dots, N_h$, $\|x^{h,i+1}\| \le (1+h\psi_x)\|x^{h,i}\| + h\psi_x$, yielding

$$
\begin{align*}
\| x^{h,i+1} \| & \le (1+h \psi_x)^{i+1} \| x^0 \| + h \psi_x \sum_{j=0}^{i} (1+h \psi_x)^j \\
& \le (1+h \psi_x)^{i+1} \| x^0 \| + h \psi_x \frac{(1+h \psi_x)^{i+1}-1}{h \psi_x} \\
& \le e^{h(i+1)} \psi_x \| x^0 \| + e^{h(i+1)} \psi_x - 1 \\
& \le e^{T} \psi_x \| x^0 \| + e^{T} \psi_x - 1,
\end{align*}
$$

where the last inequality holds because $(i+1)h \le (N_h+1)h \le T$. Consequently, with $c_{1,x} \equiv e^{T\psi_x} \equiv c_{0,x} + 1$, the bound on $\|x^{h,i+1}\|$ in (7.5) holds. The desired bound on $\|u^{h,i+1}\|$ follows from (7.10). $\square$

The next lemma establishes the two bounds in (7.5) under a linear growth
property on the solutions to the VIs.

**Lemma 7.3** Let (f, B, G) satisfy conditions (A) and (B). Suppose that F is continuous and that for some constant $\rho > 0$, the solutions to the VI ($K, q + F$) satisfies the linear growth property (6.5) for all $q \in G(\Omega)$. There exist positive scalars $c_{0,x}$, $c_{1,x}$, $c_{0,u}$, and $c_{1,u}$ such that (7.5) holds for all $h \in (0, h_1]$ and all $i = 0, 1, \dots, N_h$.

*Proof* Throughout the proof below, the scalar $h > 0$ is taken to be sufficiently small. Letting

$$
\rho_x \equiv \frac{\rho_f + \sigma_B}{1 - h_0 (1 - \theta) \rho_f}, \qquad (7.12)
$$

it follows from (7.4) that

$$
\| x^{h,i+1} - x^{h,i} \| \le h \rho_x (1 + \| x^{h,i} \| + \| u^{h,i+1} \|). \quad (7.13)
$$

By the linear growth of solutions to the VI in (7.2), we have

$$
\begin{align*}
\| u^{h,i+1} \| &\le \rho (1 + \| G(t_{h,i+1}, x^{h,i+1}) \|) \\
&\le \rho [1 + \rho_G (1 + \| x^{h,i+1} \|)] \\
&\le \rho [1 + \rho_G \{1 + \| x^{h,i} \| + h \rho_x (1 + \| x^{h,i} \| + \| u^{h,i+1} \|) \}].
\end{align*}
$$

Consequently, letting

$$
\rho_u^{\text{IVP}_1} \equiv 2 \rho (1 + \rho_G + \rho_G \rho_x)
$$

we obtain, for all $h \in (0, 1)$ satisfying $h < 1/(2\rho\rho_G\rho_x)$,

$$
\| u^{h,i+1} \| \le \rho_u^{\text{IVP}_1} (1 + 2 \| x^{h,i} \|).
$$
Therefore, with $\psi_x^{\text{IVP}_1} \equiv \rho_x(1 + 2\rho_u^{\text{IVP}_1})$, we deduce

$$ \| x^{h,i+1} - x^{h,i} \| \le h \psi_x^{\text{IVP}_1} (1 + \| x^{h,i} \|). $$

Since the conditions of Lemma 7.2 hold, the desired conclusions follow. □

The proof of the existence of a pair $(x^{h,i+1}, u^{h,i+1})$ satisfying (7.2) turns out to be a nontrivial task. For this purpose and also for later use, we review two advanced topological fixed-point theorems for set-valued maps: one due to Eilenberg-Montgomery for acyclic maps [37] and the other due to Górniewicz for admissible maps [47,48]; the recent paper [86] contains a good summary of these results and their extensions. To state the two fixed-point theorems, we first review some relevant concepts, which can be found in the text [48].

An acyclic set is a topological space $X$ where the rational homology groups of $X$ are isomorphic to those of a singleton. An absolute retract (AR) is a topological space $X$ such that if $X$ is embedded as a closed subset $X'$ of a space $Y$, then $X'$ is a retract of $Y$. Every compact convex set is an AR; every homeomorphic image of a compact convex set is acyclic. An acyclic map is an upper semicontinuous set-valued map which has compact acyclic values. A single-valued map $p: X \to Y$ is a Vietoris map if it is continuous and for each $y \in Y$, $p^{-1}(y)$ is a nonempty compact acyclic set. In Górniewicz' terminology an admissible map is an upper semicontinuous set-valued map $F: X \to Y$ with compact values which can be represented as $F = r \circ p^{-1}$ where $r: Z \to Y$ is a continuous single-valued map and $p: Z \to X$ is a Vietoris map. It can be easily shown that acyclic maps are necessarily admissible maps; furthermore, compositions of admissible maps are admissible.

The following are the two fixed-point theorems that we need, with the first one being a special case of the second.

**Theorem 7.2 (Eilenberg-Montgomery)** Every acyclic multifunction $F: X \to X$ on a compact AR $X$ has a fixed point: $x \in F(x)$ for some $x \in X$.

**Theorem 7.3 (Górniewicz)** Every admissible multifunction $F: X \to X$ on a compact AR $X$ has a fixed point: $x \in F(x)$ for some $x \in X$.

Returning to the IVP iteration (7.2), we note that the constant $\psi_x$ in the proof of Lemma 7.3 plays an important role in the subsequent analysis; for this reason, we display again this important constant

$$ \psi_x^{\text{IVP}_1} \equiv \frac{\rho_f + \sigma_B}{1 - h_0 (1 - \theta) \rho_f} [1 + 4 \rho (1 + \rho_G + \rho_G \rho_x)], $$

and recall that (i) $h_0$ is the constant asserted by Lemma 7.1, (ii) $\rho_f$ and $\rho_G$ are the respective linear growth constants of $f$ and $G$ implied by assumption (A); cf. (6.3); (iii) $\sigma_B$ is the bound of $\|B(t, x)\|$ on $\Omega$ in assumption (B), (iv) $\rho$ is the linear growth constant of the solutions to the VIs $(K, q + F)$, (v) $\rho_x$ is the constant defined in (7.12), and (vi) $\theta \in [0, 1]$ is the parameter used in the iterative scheme
(7.2). Most importantly, $\psi_x^{\text{IVP}_1}$ is independent of the step size $h$ and any other properties of the functions $(f, B, G, F)$. The role of $\psi_x^{\text{IVP}_1}$ can be seen in the next result, which establishes the existence of a pair $(x^{h,i+1}, u^{h,i+1})$ satisfying (7.2) for all $h > 0$ sufficiently small and all $i = 0, 1, \dots, N_h$.

**Proposition 7.1** Let $K \subseteq \mathbb{R}^m$ be a nonempty closed convex set. Let $(f, B, G)$ satisfy conditions (A) and (B). Suppose that $F$ is continuous and that for some constant $\rho > 0$, $\text{SOL}(K, q + F)$ is a nonempty convex set satisfying the linear growth property (6.5) for all $q \in G(\Omega)$. A scalar $h_R > 0$ exists such that for all $h \in (0, h_R]$, $\theta \in [0, 1]$, and $x^0 \in \mathbb{R}^n$, a pair $(x^{h,i+1}, u^{h,i+1})$ exists satisfying (7.2) for every $i = 0, 1, \dots, N_h$.

*Proof Associated with the constant* $\psi_x^{\text{IVP}_1}$, we define, for any scalar $h > 0$ sufficiently small, the scalars $\rho_1, \rho_2, \dots$, and $\rho_{N_h+1}$ recursively by

$$ \rho_{i+1} \equiv (1 + h \psi_x^{\text{IVP}_1}) \rho_i + h \psi_x^{\text{IVP}_1}, \quad i = 0, 1, \dots, N_h, $$

where $\rho_0$ is arbitrary. By the proof of Lemma 7.3, we can show that

$$ \rho_i \le e^T \psi_x^{\text{IVP}_1} \rho_0 + e^T \psi_x^{\text{IVP}_1} - 1, \quad \forall i = 0, 1, \dots, N_h + 1. $$

Let $\alpha$ denote the quantity on the right-hand side, which depends additionally on $\rho_0$, but is also independent of $h$. Let $h_R > 0$ satisfy $h_R < \min(h_0, h_1)$ and also

$$ h_R \frac{\rho_f (1 + \alpha) + \sigma_B \rho \rho_G (1 + 2\alpha)}{1 - h_R (1 - \theta) \rho_f} < \alpha, $$

where $h_0$ and $h_1$ are as described in Lemmas 7.1 and 7.3, respectively. Fix an arbitrary $h \in (0, h_R]$. We use induction on $i$ and a fixed-point argument to show that (7.2) has a solution $(x^{h,i+1}, u^{h,i+1})$ and any such solution must satisfy $\|x^{h,i+1}\| \le \rho_{i+1}$. Let $\mathcal{B}_\alpha$ denote the Euclidean ball in $\mathbb{R}^n$ with center at the origin and radius $2\alpha$. For any $x \in \mathcal{B}_\alpha$, let $S(t, x)$ denote the nonempty convex set $\text{SOL}(K, G(t, x) + F)$. By the linear growth assumption, we have

$$ \sup\{ \|u\| : u \in S(t, x) \} \le \rho \rho_G (1 + 2\alpha). $$

Define the map $S^i$ from $\mathcal{B}_\alpha$ into subsets of $\mathcal{B}_\alpha$ as follows: for $x \in \mathcal{B}$,

$$ S^i(x) = (I - h f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)\cdot))^{-1} (x^{h,i} + h B(t_{h,i}, x^{h,i}) S(t_{h,i+1}, x)). $$

Since the map $(I - h f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)\cdot))^{-1}$ is a homeomorphism for all $h > 0$ sufficiently small, it follows that $S^i(x)$ is a compact acyclic set. We need to show that $S^i(x)$ is a subset of $\mathcal{B}$. Let $\tilde{x}$ be an arbitrary element in $S^i(x)$ and let $u \in S(t_{h,i+1}, x)$ be such that

$$ \tilde{x} = x^{h,i} + h[f(t_{h,i+1}, \theta x^{i} + (1-\theta)\tilde{x}) + B(t_{h,i}, x^{h,i})u] $$
From (7.4), we have

$$ \| \tilde{x} - x^{h,i} \| \le h \frac{\rho_f (1 + \| x^{h,i} \|) + \sigma_B \| u \|}{1 - h(1 - \theta) \rho_f}, $$

which yields, by induction hypothesis,

$$ \| \tilde{x} \| \le \rho_i + h \frac{\rho_f (1 + \rho_i) + \sigma_B \rho \rho_G (1 + 2\alpha)}{1 - h(1 - \theta) \rho_f} < 2\alpha, $$

because $\|x^{h,i}\| \le \rho_i \le \alpha$. Since the solution map $S(t_{h,i+1}, x)$ is upper semicontinuous by Lemma 6.2, it follows that $S^i : B_\alpha \to B_\alpha$ is a closed set-valued map with compact acyclic values. By Theorem 7.2, $S^i$ has a fixed point, which easily produces a solution pair $(x^{h,i+1}, u^{h,i+1})$ of (7.2). It remains to show that $\|x^{h,i+1}\| \le \rho_{i+1}$. But this follows from the proof of Lemma 7.3 and the definition of $\rho_{i+1}$. $\square$

We now have all the preparatory results to establish the promised corollary of Theorem 7.1. The assumptions (b), (c), and (d) below correspond to (c), (d), and (e) in Theorem 6.1, respectively. Although part (a) is covered by the global reduction to a Lipschitz ODE, we include it here for completeness.

**Theorem 7.4** Let $K \subseteq \mathbb{R}^m$ be a nonempty closed convex set. Let $(f, B, G)$ satisfy conditions (A) and (B). Under any one of the four conditions below:

(a) $F$ is Lipschitz continuous and strongly monotone on $K$;

(b) $K$ contains the origin, $F(u) = Du$ for some positive semidefinite matrix $D$ such that $(K, D)$ is an $R_0$ pair;

(c) $K$ is a polyhedron containing the origin, $F(u) = Du$ for some positive semidefinite matrix $D$, and $G(\Omega) \subseteq \text{int} K(K, D)^*$;

(d) $K$ is a polyhedron containing either the origin or no line and $F(u) = Du + \Upsilon(Du)$, where $\Upsilon: \mathbb{R}^m \to \mathbb{R}^m$ is Lipschitz continuous such that with $\Phi(u) = \Upsilon(Du)$, the pair $(K, F)$ satisfies condition (e) in Theorem 6.1;

the iterates (7.1) are well defined for all $h > 0$ sufficiently small; moreover, the conclusion of Theorem 7.1 holds.

*Proof* In view of Propositions 6.4 and 6.5, the well-definedness of the iterates (7.1) follows from Proposition 7.1. Moreover, by Lemma 7.3, it is easy to see that the remaining assumptions in Theorem 7.1 are satisfied. In particular, one can verify that if $F$ is Lipschitz continuous and strongly monotone on $K$, then (7.6) holds with $E$ being the identity matrix. For the remaining parts in the present theorem, we apply Theorem 7.1(b). $\square$

Like part (d) of Theorem 6.1, part (c) of Theorem 7.4 is applicable to the generalized LCS (6.15), and thus to ODE (4.1) in particular, provided that $C$ is polyhedral, $D$ is positive semidefinite, and $q + C\mathbb{R}^n \subseteq \text{int } K(C, D)^*$.
## 8 A class of nonlinear DVIs yielding convergence

The class of DVIs in Theorem 7.4 is restricted by a uniform linear growth property of the solutions to the VIs. In turn, this property is ensured by restricting the pair $(K, F)$. In what follows, we extend the convergence analysis by relaxing the assumptions on this pair and shift the restriction to the function $G$. A distinguished feature of the development in this section is that $(K, F)$ is a non-affine pair. It turns out, however, that we need $K$ to be a cone; see Theorem 8.1. Specifically, we begin by assuming

(C) $F$ is continuous and there exists $u^{\text{ref}} \in K$,

$$ \liminf_{\substack{u \in K \\ \|u\| \to \infty}} \frac{(u - u^{\text{ref}})^T F(u)}{\|u\|^2} \geq 0; $$

(D) there exists a constant $\eta_G > 0$ such that

$$ (u - u')^T (G(t, r + B(t_{\text{ref}}, x^{\text{ref}})u) - G(t, r + B(t_{\text{ref}}, x^{\text{ref}})u')) \\ \geq \eta_G \|u - u'\|^2, $$

for all $r \in \mathbb{R}^n$, $u$ and $u'$ in $\mathbb{R}^m$, $t$ and $t_{\text{ref}}$ in $[0, T]$, and $x^{\text{ref}} \in \mathbb{R}^n$.

The liminf condition in (C) obviously holds if $F$ is monotone on $K$. Specialized to the LCS (2.8) where $K = \mathbb{R}_+^m$, condition (C) holds if the matrix $D$ is copositive (i.e., $u^T Du \geq 0$ for all $u \geq 0$); whereas condition (D) holds if and only if the matrix $CB$ is positive definite. In particular, $D$ is not required to be an $\mathbb{R}_0$ matrix, as required by conditions (b) and (d) in Theorem 7.4.

Before analyzing the convergence of the scheme (7.2), we show that for a differentiable function $G$, (D) holds if and only if the matrix product $J_x G(t, x)$ $B(t_{\text{ref}}, x^{\text{ref}})$ is uniformly positive definite for all tuples $(t, t_{\text{ref}}, x, x^{\text{ref}})$. This extends the above remark for the LCS.

**Proposition 8.1** Suppose that $G$ is continuously differentiable. Condition (D) holds if and only if for all $(t, t_{\text{ref}}, x, x^{\text{ref}}, v) \in [0, T]^2 \times \mathbb{R}^{2n+m}$,

$$ v^T J_x G(t, x) B(t_{\text{ref}}, x^{\text{ref}}) v \geq \eta_G \|v\|^2. \quad (8.1) $$

*Proof* Suppose (D) holds. Let $u = u' + \tau v$ for $\tau > 0$. We then have

$$ \begin{aligned} \eta_G \tau^2 \|v\|^2 &\geq \tau v^T \left[ G(t, r + B(t_{\text{ref}}, x^{\text{ref}})(u' + \tau v)) - G(t, r + B(t_{\text{ref}}, x^{\text{ref}})u') \right] \\ &= \tau^2 \left[ v^T J_x G(t, r + B(t_{\text{ref}}, x^{\text{ref}})u') B(t_{\text{ref}}, x^{\text{ref}})v + o(\tau)/\tau \right]. \end{aligned} $$
Dividing by $\tau^2$ and letting $\tau \downarrow 0$ easily yields (8.1). Conversely, suppose (8.1) holds. By the mean-value theorem, we can write

$$ (u - u')^T (G(t, r + B(t_{\text{ref}}, x^{\text{ref}})u) - G(t, r + B(t_{\text{ref}}, x^{\text{ref}})u')) \\ = \int_0^1 (u - u')^T J_x G(t, r + B(t_{\text{ref}}, x^{\text{ref}})(u' + \tau(u - u')) B(t_{\text{ref}}, x^{\text{ref}})(u - u') d\tau \\ \geq \eta_G \|u - u'\|^2, $$

where the inequality follows from (8.1). □

Proposition 8.1 suggests a variation of the iteration (7.2) that is perhaps more suitable for a DVI with a differentiable function G. Specifically, consider

$$ x^{h,i+1} = x^{h,i} + h[f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1}) + B(t_{h,i}, x^{h,i})u^{h,i+1}], \quad (8.2) \\ u^{h,i+1} \in \text{SOL}(K, G^{h,i+1} + F), $$

where $G^{h,i+1}$ is the linearization of G at the pair $(t_{h,i}, x^{h,i})$:

$$ G^{h,i+1}(x) = G(t_{h,i}, x^{h,i}) + \frac{\partial G(t_{h,i}, x^{h,i})}{\partial t} h + J_x G(t_{h,i}, x^{h,i})(x - x^{h,i}). $$

With this modified scheme, condition (D) can be replaced by the uniform positive definiteness of $J_x G(t, x)B(t, x)$ for all $(t, x) \in \Omega$, which is a weaker condition than (8.1). For the sake of uniformity, we continue to analyze the convergence of the basic iteration (7.2) under condition (D), which does not require the differentiability of G. It is easy to adopt the same line of analysis to the modified scheme (8.2) for a differentiable G using the weaker positive definiteness condition.

The next result establishes the existence of a pair of iterates $(x^{h,i+1}, u^{h,i+1})$ satisfying (7.2) under assumptions (A)–(D), and uniqueness of this pair under the additional monotonicity of F.

**Proposition 8.2** Let $\emptyset \neq K \subseteq \mathbb{R}^m$ be closed and convex. Let $(f, B, G, F)$ satisfy conditions (A)–(D). If

$$ 0 < h < \min \left( \frac{\eta_G}{(1-\theta) L_f (L_G \sigma_B + \eta_G)}, \frac{1}{(1-\theta) \rho_f} \right), \quad (8.3) $$

then there exists a pair $(x^h, u^h)$ satisfying

$$ x^h = x^{\text{ref}} + h \left[ f(t, \theta x^{\text{ref}} + (1-\theta)x^h) + B(t_{\text{ref}}, x^{\text{ref}})u^h \right], \quad (8.4) \\ u^h \in \text{SOL}(K, G(t, x^h) + F). $$

Moreover, the pair $(x^h, u^h)$ is unique if F is additionally monotone on K.
Proof Condition (8.3) clearly implies $1 - h(1 - \theta)L_f > 0$. Hence, for every vector $u$, there exists a unique vector $x(u)$ satisfying

$$x(u) = x^{\text{ref}} + h \left[ f(t, \theta x^{\text{ref}} + (1-\theta)x(u)) + B(t_{\text{ref}}, x^{\text{ref}})u \right]. \quad (8.5)$$

We claim that the map $u \mapsto G(t, x(u))$ is strongly monotone on $\mathbb{R}^m$, provided that $h$ satisfies (8.3). Once this is established, the existence of the pair $(x^h, u^h)$ satisfying (8.4) follows readily (see [38, Proposition 2.2.7]), because the map $u \mapsto G(t, x(u)) + F(u)$ is then coercive on $K$, meaning

$$\liminf_{\substack{u \in K \\ \|u\| \to \infty}} \frac{(u - u^{\text{ref}})^T (G(t, x(u)) + F(u))}{\|u\|^2} > 0.$$

Moreover, if $F$ is monotone on $K$, then the map $u \mapsto G(t, x(u)) + F(u)$ is strongly monotone; the uniqueness of $(x^h, u^h)$ is then obvious. To establish the desired claim, we note that $G(t, x(u)) = G(t, g_h(u) + hB(t_{\text{ref}}, x^{\text{ref}})u)$, where $g_h(u) = x^{\text{ref}} + hf(t, \theta x^{\text{ref}} + (1-\theta)x(u))$. By (8.5) and (7.4), we have

$$\|g_h(u) - g_h(u')\| \le h(1-\theta)L_f \|x(u) - x(u')\| \le \frac{h^2(1-\theta)L_f \sigma_B \|u - u'\|}{1 - h(1-\theta)L_f}.$$

Moreover, by assumption (D),

$$(u - u')^T [ G(t, g_h(u) + hB(t_{\text{ref}}, x^{\text{ref}})u) - G(t, g_h(u) + hB(t_{\text{ref}}, x^{\text{ref}})u') ] \\ \geq h \eta_G \|u - u'\|^2.$$

Consequently,

$$\begin{align*}
(u - u')^T [ G(t, x(u)) - G(t, x(u')) ]
&= (u - u')^T [ G(t, g_h(u) + hB(t_{\text{ref}}, x^{\text{ref}})u) - G(t, g_h(u) + hB(t_{\text{ref}}, x^{\text{ref}})u') ] \\
&= (u - u')^T [ G(t, g_h(u) + hB(t_{\text{ref}}, x^{\text{ref}})u) - G(t, g_h(u) + hB(t_{\text{ref}}, x^{\text{ref}})u') ] \\
&\quad + (u - u')^T [ G(t, g_h(u) + hB(t_{\text{ref}}, x^{\text{ref}})u') - G(t, g_h(u) + hB(t_{\text{ref}}, x^{\text{ref}})u') ] \\
&\geq h \eta_G \|u - u'\|^2 - L_G \|g_h(u) - g_h(u')\| \|u - u'\| \\
&\geq h \left[ \eta_G - \frac{h(1-\theta)L_f L_G \sigma_B}{1 - h(1-\theta)L_f} \right] \|u - u'\|^2.
\end{align*}$$

By (8.3), it follows that the map $u \mapsto G(t, x(u))$ is strongly monotone. $\square$

<figure><img src="image_2.png" /></figure>

Our next task is to establish that the pair of solutions $(x^h, u^h)$ to (8.4) is uniformly bounded for all $h > 0$ sufficiently small. For this purpose, we need to restrict $F$ to be a strongly monotone composite map. Specifically, we assume that
(C') $F \equiv E^T \circ \Psi \circ E$, where $E \in \mathbb{R}^{\ell \times m}$ and $\Psi : \mathbb{R}^{\ell} \to \mathbb{R}^{\ell}$ is Lipschitz continuous and strongly monotone on $E\mathbb{R}^m$.

The most interesting case of the above assumption is when $\ker E \neq \{0\}$. Indeed, when the latter kernel is a singleton, then $F$ is strongly monotone on $\mathbb{R}^m$, and this case is covered by Theorem 7.4(a). Ideally, we would like to make only the additional assumption that $K_\infty \cap \ker E = \{0\}$, as in part (b) of Theorem 6.1. Nevertheless, the latter assumption does not seem to be enough to ensure the key bound (7.6). Consequently, we need to further restrict $K$ and $E$. Specifically, let $\mathcal{E}$ and $\mathcal{E}^\perp$ denote, respectively, the null space of $E$ and its orthogonal complement. Let $K_1$ and $K_2$ denote the orthogonal projection of $K$ onto $\mathcal{E}$ and $\mathcal{E}^\perp$, respectively. We clearly have $K \subseteq K_1 \oplus K_2$, where $\oplus$ denotes the orthogonal sum. We postulate the following assumption:

(E) $K_1 \oplus K_2 \subseteq K$, or equivalently, $K_1 \oplus K_2 = K$.

Assumption (E) is trivially valid in the case where $E$ is the zero matrix; i.e., when $F$ is identically equal to zero. In this case, much of the discussion below holds vacuously; most importantly, the main result of this subsection, Theorem 8.1, is valid. In general, let $Z \in \mathbb{R}^{m \times k}$ be any matrix whose columns form an orthonormal basis of $\mathcal{E}$ so that $EZ = 0$. Extend the basis $Z$ through $\mathcal{E}^\perp$ to an orthonormal basis $P$ of $\mathbb{R}^m$. Thus $P = [Z \ W] \in \mathbb{R}^{m \times m}$ for some matrix $W \in \mathbb{R}^{m \times (m-k)}$ whose columns form an orthonormal basis of $\mathcal{E}^\perp$. (When $E = 0$, we have $k = m$; the matrix $W$ and all matrices and vectors related to it that appear subsequently are vacuous.) Hence $P$ is an orthogonal matrix; i.e., $P^T P$ is the identity matrix of order $m$. It follows that the matrix $\tilde{W} \equiv EW \in \mathbb{R}^{\ell \times (m-k)}$ must have full column rank; moreover, $E\mathbb{R}^m = \tilde{W}\mathbb{R}^{m-k}$. In terms of the matrices $Z$ and $W$, assumption (E) is equivalent to

$$ZZ^T K \oplus WW^T K = K.$$

This follows because $ZZ^T$ and $WW^T$ are the Euclidean projectors onto $\mathcal{E}$ and $\mathcal{E}^\perp$, respectively. The following lemma rephrases assumption (E) in a form that is more conducive for the subsequent analysis.

**Lemma 8.1** *Condition (E) is equivalent to*

$$\left[ \begin{pmatrix} \mu^i \\ \lambda^i \end{pmatrix} \in P^T K, \quad i = 1, 2 \right] \Rightarrow \begin{pmatrix} \mu^1 \\ \lambda^2 \end{pmatrix} \in P^T K, \qquad (8.6)$$

where the partition of the vectors is according to that of the matrix $P$ in terms of the two submatrices $Z$ and $W$.

*Proof* Suppose (E) holds. Let $(\mu^i, \lambda^i)$ satisfy the left-hand condition in (8.6). We then have $Z\mu^i + W\lambda^i \in K$. By (E), it follows that

$$ZZ^T(Z\mu^1 + W\lambda^1) + WW^T(Z\mu^2 + W\lambda^2) \in K,$$
the vector on the left-hand side is clearly equal to $Z\mu^1 + W\lambda^2$. Thus the right-hand side of (8.6) follows. Conversely, suppose the latter implication holds. Let $u^i \in K$ for $i = 1, 2$ and write

$$ \begin{pmatrix} \mu^i \\ \lambda^i \end{pmatrix} \equiv P^T u^i = \begin{pmatrix} Z^T u^i \\ W^T u^i \end{pmatrix}. $$

Thus, $ZZ^T u^1 + WW^T u^2 = Z\mu^1 + W\lambda^2$; since (8.6) holds, the vector in the right-hand side of the last expression belongs to $K$. $\square$

Since $\Psi$ is strongly monotone on $E\mathfrak{R}^m$ and $\tilde{W}$ has full column rank, it follows that the map $\Upsilon \equiv \tilde{W}^T \circ \Psi \circ \tilde{W}$ is a strongly monotone self-map from $\mathfrak{R}^{m-k}$ onto itself. Thus, a constant $\eta_\Upsilon > 0$ exists such that

$$ (\Upsilon(\lambda) - \Upsilon(\lambda'))^T (\lambda - \lambda') \geq \eta_\Upsilon \|\lambda - \lambda'\|^2, \quad \forall \lambda, \lambda' \in \mathfrak{R}^{m-k}. \quad (8.7) $$

(When $E = 0$, we let $\eta_\Upsilon$ be an arbitrary constant.) In what follows, we analyze the pair $(x^h, u^h)$ under the assumptions made so far. For this purpose, we establish a preliminary lemma that pertains to a VI $(K, q + \Phi_h + F)$, where $\Phi_h$ is intended to be a mapping of the kind $G \circ x^h$, with $x^h(u)$ being the unique vector $x$ satisfying the equation

$$ x = x^{\text{ref}} + h[f(t, x^{\text{ref}} + (1-\theta)(x - x^{\text{ref}})) + B(t_{\text{ref}}, x^{\text{ref}})u]; $$

cf. the proof of Proposition 8.2. With the nonsingular transformation of variables: $u \equiv Pv = Z\mu + W\lambda$, where $v \equiv (\mu, \lambda) \in \mathfrak{R}^m$ the VI $(K, q + \Phi_h + F)$ is equivalent to (since $P^T = P^{-1}$): $(\mu, \lambda) \in P^T K$ and for all $(\mu', \lambda') \in P^T K$,

$$ \begin{pmatrix} \mu' - \mu \\ \lambda' - \lambda \end{pmatrix}^T \begin{pmatrix} Z^T q + Z^T \Phi_h (Z\mu + W\lambda) \\ W^T q + \Upsilon(\lambda) + W^T \Phi_h (Z\mu + W\lambda) \end{pmatrix} \geq 0. $$

Notice that $Eu = \tilde{W}\lambda$. Based on the above preparations, we state and prove the following result.

**Lemma 8.2** Let $F$ satisfy condition (C'). Let $\Phi_h: \mathfrak{R}^m \to \mathfrak{R}^m$ be strongly monotone on $\mathfrak{R}^m$ with constant $h \eta_\phi > 0$ and Lipschitz continuous with constant $h L_\phi > 0$. Provided that $h > 0$ is sufficiently small satisfying

$$ 0 < h < \frac{\eta_\Upsilon}{L_\phi \|W\|^2 \left(1 + \frac{L_\phi}{\eta_\phi}\right)}, \quad (8.8) $$

the map $q \mapsto \text{ESOL}(K, q + \Phi_h + F)$ is single-valued and Lipschitz continuous with a constant independent of $h$.
Proof For $i = 1, 2$, let $u^i \in \text{SOL}(K, q^i + \Phi_h + F)$ and write $P^Tu^i \equiv (\mu^i, \lambda^i)$. By Lemma 8.1, $(\mu^1, \lambda^2)$ and $(\mu^2, \lambda^1)$ both belong to $P^TK$. Hence

$$ 
\begin{align*}
(\mu^2 - \mu^1)^T [Z^T q^1 + Z^T \Phi_h (Z\mu^1 + W\lambda^1)] &\geq 0, \\
(\lambda^2 - \lambda^1)^T [W^T q^1 + \Upsilon(\lambda^1) + W^T \Phi_h (Z\mu^1 + W\lambda^1)] &\geq 0, \\
(\mu^1 - \mu^2)^T [Z^T q^2 + Z^T \Phi_h (Z\mu^2 + W\lambda^2)] &\geq 0, \\
(\lambda^1 - \lambda^2)^T [W^T q^2 + \Upsilon(\lambda^2) + W^T \Phi_h (Z\mu^2 + W\lambda^2)] &\geq 0.
\end{align*}
$$

Adding the first and third inequality and rearranging terms, we obtain

$$ 
\begin{align*}
& (q^2 - q^1)^T Z(\mu^1 - \mu^2) \\
& \geq (Z(\mu^2 - \mu^1))^T [\Phi_h(Z\mu^2 + W\lambda^2) - \Phi_h(Z\mu^1 + W\lambda^1)] \\
& = (Z(\mu^2 - \mu^1))^T [\Phi_h(Z\mu^2 + W\lambda^2) - \Phi_h(Z\mu^2 + W\lambda^1) \\
& \qquad + \Phi_h(Z\mu^2 + W\lambda^1) - \Phi_h(Z\mu^1 + W\lambda^1)] \\
& \geq \eta_\phi h \|Z(\mu^2 - \mu^1)\|^2 - L_\psi h \|Z(\mu^2 - \mu^1)\| \|W(\lambda^1 - \lambda^2)\|.
\end{align*}
$$

Consequently,

$$ \| Z(\mu^2 - \mu^1) \| \leq \frac{1}{\eta_\phi h} \| q^1 - q^2 \| + \frac{L_\phi}{\eta_\phi} \| W(\lambda^1 - \lambda^2) \|. $$

Similarly, we can deduce

$$ 
\begin{align*}
(q^2 - q^1)^T W(\lambda^1 - \lambda^2) &\geq \left[ \eta_\gamma - h \|W\|^2 \left( \frac{L_\phi^2}{\eta_\phi} + L_\phi \right) \right] \|\lambda^1 - \lambda^2\|^2 \\
&\quad - \|W(\lambda^2 - \lambda^1)\| \frac{L_\phi}{\eta_\phi} \|q^1 - q^2\|.
\end{align*}
$$

Consequently, if (8.8) holds, then

$$ \| \lambda^2 - \lambda^1 \| \leq \frac{\left(1 + \frac{L_\phi}{\eta_\phi}\right) \| W \|}{\eta_\gamma - h L_\phi \| W \|^2 \left(1 + \frac{L_\phi}{\eta_\phi}\right)} \| q^1 - q^2 \|, $$

which yields

$$ \| Eu^1 - Eu^2 \| \leq \frac{\| \tilde{W} \| \left(1 + \frac{L_\phi}{\eta_\phi}\right) \| W \|}{\eta_\gamma - h L_\phi \| W \|^2 \left(1 + \frac{L_\phi}{\eta_\phi}\right)} \| q^1 - q^2 \|, $$
where the multiplicative constant on the right-hand side is less than a constant that is independent of $h$. $\square$

Our next result, which establishes a uniform bound on the family $\{u^h\}$, requires $K$ to be a cone.

**Proposition 8.3** Let $K$ be a closed convex cone in $\mathbb{R}^m$ and $\theta \in [0, 1]$ be a given scalar. Let $(f, B, G)$ satisfy the conditions (A), (B) and (D). Let $F$ satisfy (C') and let $\eta_\tau$ satisfy (8.7). Assume $\Psi(0) = 0$ and condition (E). Let $x^{\text{ref}} \in \mathbb{R}^n$ be an arbitrary vector in $\mathbb{R}^n$. If $h$ satisfies (8.3), then for all $t$ and $t_{\text{ref}}$ in $[0, T]$, a unique pair $(x^h, u^h)$ exists satisfying (8.4). If in addition there exists $u^{\text{ref}} \in K$ such that $G(t_{\text{ref}}, x^{\text{ref}}) + F(u^{\text{ref}})$ is in $K^*$, and $h > 0$ is sufficiently small satisfying in particular

$$0 < h < \frac{\eta_{\tau}}{\|W\|^2 L_G \sigma_B \left(1 + \frac{L_G \sigma_B}{\eta_G}\right)}, \quad (8.9)$$

then for some constant $\xi > 0$ that is independent of $h$ and $x^{\text{ref}}$, $\|u^h\|$ is bounded above by

$$\xi \left[ \left(1 + \frac{1}{h}\right) |t - t'| + \|G(t, x^{\text{ref}})\| + h \left\{ \|f(t, x^{\text{ref}})\| + \|x^h - x^{\text{ref}}\| \right\} \right]. \quad (8.10)$$

*Proof* It suffices to prove (8.10). Write $P^Tu^h = (\mu^h, \lambda^h)$. Since $K$ contains the origin, Lemma 8.1 implies that $P\begin{pmatrix} \mu^h \\ 0 \end{pmatrix} = Z\mu^h$ belongs to $K$, Therefore,

$$0 \le (G(t_{\text{ref}}, x^{\text{ref}}) + E^T \Psi(Eu^{\text{ref}}))^T Z \mu^h = (Z \mu^h)^T G(t_{\text{ref}}, x^{\text{ref}}). \quad (8.11)$$

The remainder of the proof is very similar to that of the Lipschitzian proof in Lemma 8.2. By the proof of Proposition 8.2, we have

$$\begin{aligned} & (\mu - \mu^h)^T Z^T G(t, x^h) \ge 0, \\ & (\lambda - \lambda^h)^T [\Upsilon(\lambda^h) + W^T G(t, x^h)] \ge 0, \end{aligned} \quad (8.12)$$

for all $(\mu, \lambda) \in P^T K$. Since $K$ contains the origin, we have, by substituting $\mu = 0$ in the first of the two inequalities in (8.12) and by (8.11),

$$0 \ge (\mu^h)^T Z^T G(t, x^h) \ge (Z\mu^h)^T [G(t, x^{\text{ref}}) - G(t_{\text{ref}}, x^{\text{ref}}) + G(t, x^h) - G(t, x^{\text{ref}})].$$

Since $x^h = x^{\text{ref}} + h[f(t, x^{\text{ref}} + (1-\theta)(x^h - x^{\text{ref}})) + B(t_{\text{ref}}, x^{\text{ref}})(Z\mu^h + W\lambda^h)]$, we deduce,

$$\begin{aligned} 0 \ge (Z\mu^h)^T [ & G(t, x^{\text{ref}}) - G(t_{\text{ref}}, x^{\text{ref}}) \\ & + G(t, x^{\text{ref}}) + h(f(t, x^{\text{ref}} + (1-\theta)(x^h - x^{\text{ref}})) + B(t_{\text{ref}}, x^{\text{ref}})(Z\mu^h + W\lambda^h))) \\ & - G(t, x^{\text{ref}}) + h(f(t, x^{\text{ref}} + (1-\theta)(x^h - x^{\text{ref}})) + B(t_{\text{ref}}, x^{\text{ref}})W\lambda^h)) \\ & + G(t, x^{\text{ref}}) + h(f(t, x^{\text{ref}} + (1-\theta)(x^h - x^{\text{ref}})) + B(t_{\text{ref}}, x^{\text{ref}})W\lambda^h)) - G(t, x^{\text{ref}})] \end{aligned}$$
which implies

$$
\eta_G h \| Z\mu^h \| \le L_G |t-t'| + h L_G [ \|f(t,x^{\text{ref}})\| + L_f (1-\theta) \|x^h - x^{\text{ref}}\| + \|B(t_{\text{ref}}, x^{\text{ref}})\| \|W\lambda^h\| ]
$$

Hence $\|Z\mu^h\|$ is bounded above by

$$
\frac{L_G |t - t_{\text{ref}}|}{\eta_G h} + \frac{L_G [\|f(t, x^{\text{ref}})\| + L_f (1 - \theta) \|x^h - x^{\text{ref}}\| + \sigma_B \|W \lambda^h\|]}{\eta_G}. \quad (8.13)
$$

Substituting $\lambda = 0$ into the second inequality in (8.12) and using the assumption
$\Upsilon(0) = 0$ and the strong monotonicity of $\Upsilon$, we deduce

$$
\begin{align*}
\eta_{\Upsilon} \| \lambda^h \| &\le \| W \| \| G(t, x^h) \| \\
&\le \| W \| [ \| G(t, x^{\text{ref}}) \| + h L_G \| f(t, x^{\text{ref}} + (1-\theta)(x^h - x^{\text{ref}})) \| ] \\
&\quad + \| W \| h L_G \| B(t_{\text{ref}}, x^{\text{ref}}) \| \| Z\mu^h + W\lambda^h \|,
\end{align*}
$$

which implies

$$
( \eta_{\Upsilon} - h \|W\|^2 L_G \sigma_B ) \|\lambda^h\| \le \|W\| h L_G \sigma_B \|Z\mu^h\| \\
+ \|W\| [ \|G(t,x^{\text{ref}})\| + h L_G \{ \|f(t,x^{\text{ref}})\| + L_f (1-\theta) \|x^h - x^{\text{ref}}\| \} ].
$$

Consequently, from (8.13) and after rearranging terms, we obtain

$$
\begin{align*}
\left[ \eta_{\Upsilon} - h \|W\|^2 L_G \sigma_B \left( 1 + \frac{L_G \sigma_B}{\eta_G} \right) \right] \|\lambda^h\| \\
&\le \|W\| [ \|G(t,x^{\text{ref}})\| ] + \|W\| \frac{L_G^2 \sigma_B}{\eta_G} |t - t_{\text{ref}}| \\
&\quad + \left( 1 + \frac{L_G \sigma_B}{\eta_G} \right) \|W\| h L_G [ \|f(t,x^{\text{ref}})\| + L_f (1-\theta) \|x^h - x^{\text{ref}}\| ].
\end{align*}
$$

Therefore, if $h$ satisfies (8.9), the multiplicative factor of $\|\lambda^h\|$ in the left-hand
side is positive. Dividing by this constant factor, we obtain

$$
\| \lambda^h \| \le c [ \| G(t, x^{\text{ref}}) \| + |t - t'| + h [ \| f(t, x^{\text{ref}}) \| + \| x^h - x^{\text{ref}} \| ] ]
$$

for an appropriate constant $c > 0$ that is independent of $h$. Substituting the
above bound in (8.13), we obtain a corresponding bound for $\|Z\mu^h\|$. Since
$\|u^h\| \le \|Z\mu^h\| + \|W\| \|\lambda^h\|$, the desired upper bound (8.10) follows readily with
an appropriate constant $\xi > 0$. $\square$

Applying Proposition 8.3 to the discrete DVI (7.2), we obtain the following
convergence theorem that complements Theorem 7.4.
**Theorem 8.1** Let K be a closed convex cone in $\mathbb{R}^m$ and $\theta \in [0,1]$ be a given scalar. Let $(f, B, G)$ satisfy the conditions (A), (B) and (D). Let $F$ satisfy (C') and let $\eta_\gamma$ satisfy (8.7). Assume $\Psi(0) = 0$ and condition (E). A positive scalar $\bar{h} > 0$ exists such that for all $x^0 \in \mathbb{R}^n$ for which $\text{SOL}(K, G(t_0, x^0) + F) \neq \emptyset$ and for all $h \in (0, \bar{h}]$, a unique pair $(x^{h,i+1}, u^{h,i+1})$ exists satisfying (7.2) for all $i = 0, 1, 2, \dots, N_h$; moreover positive constants $c_{0,x}, c_{1,x}, c_{0,u}, c_{1,u}$, and $c_{2,u}$ exist such that (7.5) and (7.6) hold for all $h > 0$ sufficiently small and for all $i = 0, 1, 2, \dots, N_h$. Consequently, the conclusion of Theorem 7.1 holds for the discrete-time trajectories $\{(x^{h,i}, u^{h,i})\}$.

*Proof* Since K is a cone and $u^{h,i} \in \text{SOL}(K, G(t_{h,i}, x^{h,i}) + F)$ (for $i = 0$, the assumption yields the existence of $u^{h,0}$), it follows that $G(t_{h,i}, x^{h,i}) + F(u^{h,i}) \in K^*$. Thus, with $t = t_{h,i+1}$, $t_{\text{ref}} = t_{h,i}$, and $u^{\text{ref}} = u^{h,i}$, (8.10) gives

$$ \begin{aligned} \| u^{h,i+1} \| & \le \xi [1 + h + \|G(t_{h,i+1}, x^{h,i})\| + h \{\|f(t_{h,i+1}, x^{h,i})\| + \|x^{h,i+1} - x^{h,i}\|\}] \\ & \le \xi [h + 1 + (\rho_G + h \rho_f)(1 + \|x^{h,i}\|) + h \|x^{h,i+1} - x^{h,i}\|] \end{aligned} $$

which yields, for some constant $\omega > 0$ that is independent of $h$ and $x^0$,

$$ \| u^{h,i+1} \| \le \omega (1 + \| x^{h,i} \| + h \| x^{h,i+1} - x^{h,i} \|). \quad (8.14) $$

From (7.13), we have $\|x^{h,i+1} - x^{h,i}\| \le h\rho_x(1 + \|x^{h,i}\| + \|u^{h,i+1}\|)$, which yields

$$ \| x^{h,i+1} - x^{h,i} \| \le h \rho_x \frac{1+\omega}{1-\omega \rho_x h^2} (1 + \| x^{h,i} \|). $$

Therefore, for some positive constant $\psi_x^{\text{IVP}_2}$ that is independent of $h$, (7.11) holds for all $h > 0$ sufficiently small. Moreover, (8.14) shows that (7.10) holds for some positive constant $\rho_u^{\text{IVP}_2}$. By Lemma 7.3, it therefore follows that there exist positive scalars $c_{0,x}, c_{1,x}, c_{0,u}$, and $c_{1,u}$ and $\bar{h}$ such that (7.5) holds for all $h \in (0, \bar{h}]$ and all $i = 0, 1, \dots, N_h$. To complete the proof of the theorem, it remains to show the existence of a constant $c_{2,u} > 0$ such that (7.6) holds. For every $u \in \mathbb{R}^m$, let $\tilde{x}^{i+1}(u)$ denote the unique vector $x$ satisfying $x - x^{h,i} = h[f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x) + B(t_{h,i}, x^{h,i})u]$. By the proof of Proposition 8.2, the map $\Phi_h : u \mapsto G(t_{h,i+1}, \tilde{x}^{i+1}(u))$ is strongly monotone with modulus $h \left[ \eta_G - \frac{h(1-\theta)L_f L_G \sigma_B}{1-h(1-\theta)L_f} \right]$, which is greater than $h\eta_G/2$, provided that $h > 0$ is sufficiently small. Moreover, $\Phi_h$ is Lipschitz continuous with modulus $\frac{hL_G\sigma_B}{1-h(1-\theta)L_f}$, which is less than $2hL_G\sigma_B$, provided that $h > 0$ is sufficiently small. Applying Lemma 8.2 to the pair $(K, \Phi_h + F)$, we deduce the existence of a constant $c_u > 0$, which is independent of $h$, such that $\|Eu(q^1) - Eu(q^2)\| \le c_u \|q^1 - q^2\|$ for all $q^1$ and $q^2$, where for $i = 1,2$, $u(q^i)$ is
the unique solution of the VI $(K, q^i + \Phi_h + F)$. Since $u^{h,i}$ is a solution of the VI $(K, G(t_{h,i}, x^{h,i}) + F)$ and

$$G(t_{h,i}, x^{h,i}) = G(t_{h,i}, x^{h,i}) - G(t_{h,i+1}, \tilde{x}^{i+1}(u^{h,i})) + \Phi_h(u^{h,i}),$$

with $q^1$ being the zero vector and $q^2 = G(t_{h,i}, x^{h,i}) - G(t_{h,i+1}, \tilde{x}^{i+1}(u^{h,i}))$, we deduce

$$\begin{aligned}
\| Eu^{h,i+1} - Eu^{h,i} \| &\le c_u \| G(t_{h,i}, x^{h,i}) - G(t_{h,i+1}, \tilde{x}^{i+1}(u^{h,i})) \| \\
&\le c_u L_G [h + \| \tilde{x}^{i+1}(u^{h,i}) - x^{h,i} \|] \\
&\le h c_u L_G \left[ 1 + \frac{\rho_f (1 + \| x^{h,i} \|) + \sigma_B \| u^{h,i} \|}{1 - h \rho_f (1 - \theta)} \right],
\end{aligned}$$

where the last inequality follows from (7.4) of Lemma 7.1. The desired bound (7.6) now follows readily. □

Due to its importance, we give a corollary of Theorem 8.1 when $F$ is the zero map, which yields the following initial-value DCP:

$$\begin{array}{l}
\dot{x} = f(t,x) + B(t,x)u, \quad x(0) = x^0, \\
C \ni u \perp G(t,x) \in C^*.
\end{array}
\qquad (8.15)$$

**Corollary 8.1** Let $C$ be a closed convex cone in $\mathbb{R}^m$ and $\theta \in [0,1]$ be a given scalar. Let $(f, B, G)$ satisfy the conditions (A), (B) and (D). A positive scalar $\tilde{h} > 0$ exists such that for all $x^0 \in G(t_0, \cdot)^{-1}(C^*)$ and for all $h \in (0, \tilde{h}]$, a unique pair $(x^{h,i+1}, u^{h,i+1})$ exists satisfying for all $i = 0, 1, 2, \dots, N_h$:

$$\begin{array}{l}
x^{h,i+1} - x^{h,i} = h[f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1}) + B(t_{h,i}, x^{h,i})u^{h,i+1}] \quad (8.16) \\
C \ni u^{h+1} \perp G(t, x^{h,i+1}) \in C^*.
\end{array}$$

The discrete-time trajectory ${\{(x^{h,i}, u^{h,i})\}}$ is uniformly bounded; moreover, there is a sequence ${\{h_v\}} \downarrow 0$ such that $\tilde{x}^{h_v} \to \tilde{x}$ uniformly on $[0, T]$ and $\hat{u}^{h_v} \to \hat{u}$ weakly in $L^2(0, T)$; moreover all such limits $(\tilde{x}, \hat{u})$ are weak solutions of (8.15).

Specializing the above result to the DCP (2.9) derived from the VIE: $\dot{x} \in -F(x) - \mathcal{N}(x; K)$, we can establish the convergence of a certain iterative method for solving the VIE in the case of a cone $K$.

**Corollary 8.2** Let $K$ be a closed convex cone in $\mathbb{R}^n$ and $\theta \in [0,1]$ be a given scalar. Let $F: \mathbb{R}^n \to \mathbb{R}^n$ be a Lipschitz continuous function. A positive scalar $\tilde{h} > 0$ exists such that for all $x^0 \in K$ and for all $h \in (0, \tilde{h}]$, a unique pair $(x^{h,i+1}, u^{h,i+1})$ exists satisfying for all $i = 0, 1, 2, \dots, N_h$:

$$K^* \ni u^{h,i+1} \perp x^{h,i+1} = [I + hF(\theta x^{h,i} + (1-\theta)\cdot)]^{-1} (x^{h,i} + hu^{h,i+1}) \in K.$$
The discrete-time trajectory ${\{ (x^{h,i}, u^{h,i}) \}}$ is uniformly bounded; moreover, there is a sequence $\{h_v\} \downarrow 0$ such that $\hat{x}^{h_v} \to \hat{x}$ uniformly on $[0, T]$ and $\hat{u}^{h_v} \to \hat{u}$ weakly in $L^2(0, T)$; finally, any such limit $\hat{x}$ is a weak solution of the VIE $\dot{x} \in -F(x) - \mathcal{N}(x; K)$ satisfying the initial condition $x(0) = x^0$.

*Proof* This is an immediate consequence of Corollary 8.1 under the following identifications: $f(t,x) \equiv -F(x)$, $B(t,x) \equiv I$, $C \equiv K^*$, $G(t,x) \equiv x$. Clearly, any element of $C$ is a solution of the initial CP $(C, G(0, x^0))$. $\square$

Condition (D) turns out to be quite restrictive when specialized to the DCP (2.10) with a nonlinear (differentiable) constraint function $g$. Nevertheless, as we mentioned before, the proper way to solve the latter DCP is via the iteration (7.3), or alternatively, (8.2). Following a similar analysis, we can establish the convergence of both modified schemes for solving (2.10), provided that the gradients of the constraint functions are “uniformly linearly independent”; i.e., $Jg(x)Jg(x)^T$ is uniformly positive definite for all $x \in K$ (for (7.3), and all $x \in \mathbb{R}^n$ for (8.2), respectively). We omit the details. The case where the pair $(K, F)$ is affine is discussed in Sect. 9.2; see Theorem 9.5.

## 9 Time stepping for boundary-value DVIs

Based on the iterative scheme (3.3), we introduce a time-stepping method for solving a boundary-value DVI of the following type:

$$ 
\begin{aligned}
\dot{x} &= f(t, x) + B(t, x)u, \\
u &\in \text{SOL}(K, G(t, x) + F), \\
b &= Mx(0) + Nx(T),
\end{aligned} 
\quad (9.1) $$

where $b$ is an $n$-vector and $M$ and $N$ are $n \times n$ matrices. Specifically, for a given step size $h > 0$ and an integer $N_h$ satisfying $T = (N_h + 1)h$, we compute

$$ x^h = \begin{pmatrix} x^{h,0} \\ x^{h,1} \\ \vdots \\ x^{h,N_h+1} \end{pmatrix} \in \mathbb{R}^{(N_h+2)n} \quad \text{and} \quad u^h = \begin{pmatrix} u^{h,1} \\ u^{h,2} \\ \vdots \\ u^{h,N_h+1} \end{pmatrix} \in \mathbf{K}^h, \quad (9.2) $$

where $\mathbf{K}^h$ is the $(N_h + 1)$-fold Cartesian product of the set $K$, by solving a large-scale, aggregate VI that is defined by the following conditions: for $i = 0, 1, \dots, N_h$,

$$ 
\begin{aligned}
x^{h,i+1} &= x^{h,i} + h[f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1}) + B(t_{h,i}, x^{h,i})u^{h,i+1}], \\
u^{h,i+1} &\in \text{SOL}(K, G(t_{h,i+1}, x^{h,i+1}) + F)
\end{aligned} 
\quad (9.3) $$

plus the boundary condition:

$$ b = Mx^{h,0} + Nx^{h,N_h+1}. \quad (9.4) $$
Unlike an initial-value system, the above discrete-time VI does not decouple into sub-VIs pertaining to the individual time steps, due to the coupling equation (9.4); all components of the two unknown vectors $(\mathbf{x}^h, \mathbf{u}^h)$ have to be computed simultaneously. Detailed discussion of how to solve such (discrete-time boundary-value) VIs most effectively is beyond the scope of this paper. Our focus is on the existence and convergence of a resulting discrete-time boundary-value trajectory $(\mathbf{x}^h, \mathbf{u}^h)$.

Besides the blanket conditions (A) and (B), the assumptions needed here are of three kinds: (a) monotonicity of $F$ and linear growth of the elements in SOL($K, q + F$), (b) the strong monotonicity assumption (D) on $G$, and (c) a condition on the boundary matrices $M$ and $N$ that involves the time $T$ and the constant $\psi_x^{\text{IVP}_1}$. Together, these assumptions yield the existence of a solution to the discrete-time boundary-value system (9.3) and (9.4).

**Lemma 9.1** Let $K \subseteq \mathbb{R}^m$ be a nonempty closed convex set. Let $(f, B, G)$ satisfy conditions (A) and (B). Suppose that

(a) $F$ is continuous and monotone; a constant $\rho > 0$ exists such that solutions to the VI $(K, q + F)$ exist and satisfy the linear growth property (6.5) for all $q \in G(\Omega)$;

(b) condition (D) holds;

(c) $M + N$ is nonsingular and

$$e^{T \psi_x^{\text{IVP}_1}} < 1 + \frac{1}{\| (M+N)^{-1} N \|} \quad (9.5)$$

There exists $\bar{h} > 0$ such that for all $h \in (0, \bar{h}]$, a discrete-time boundary-value trajectory (9.2) exists that satisfies (9.3) and (9.4).

*Proof* Let $h$ satisfy (8.3). Define, for $i = 0, 1, \dots, N_h$, the continuous map $\tilde{x}^i$: $\mathbb{R}^n \to \mathbb{R}^n$ as follows: $\tilde{x}^i(x^{\text{ref}})$ is the vector $x^h$ such that for some $u^h$, $(x^h, u^h)$ is the unique pair satisfying (8.4) for $(t, t^{\text{ref}}) = (t_{h,i+1}, t_{h,i})$. The well-definedness of the map $\tilde{x}^i$ is justified by Proposition 8.2, and its continuity by an easy argument. Define the composite maps $\Gamma^i$, for $i = 0, 1, 2, \dots, N_h + 1$, recursively by: $\Gamma^0 = I$, and

$$\Gamma^{i+1} \equiv \tilde{x}^i \circ \Gamma^i, \quad i = 0, 1, \dots, N_h.$$

From Lemma 7.2 we have, for $i = 0, 1, \dots, N_h$,

$$\| \Gamma^{i+1}(x^{\text{ref}}) \| \le \rho_{i+1} \quad \text{and} \quad \| \Gamma^{i+1}(x^{\text{ref}}) - \Gamma^i(x^{\text{ref}}) \| \le \rho_{i+1} - \rho_i, \quad (9.6)$$

where $\rho_{i+1}$ satisfies the recursion

$$\rho_{i+1} \equiv (1 + h \psi_x) \rho_i + h \psi_x^{\text{IVP}_1}, \quad i = 0, 1, \dots, N_h, \quad (9.7)$$
with $\rho_0 \ge \|x^{\text{ref}}\|$. Since $Mx + Ny = b$ if and only if $x = (M+N)^{-1}b - (M+N)^{-1}N(y-x)$, defining the continuous map

$$H \equiv (M+N)^{-1}b - (M+N)^{-1}N \circ (\Gamma^{N_h+1} - \Gamma^0), \quad (9.8)$$

we see that a sufficient condition for the existence for a discrete-time trajectory (9.2) satisfying (9.3) and (9.4) is that the map $H$ has a fixed point. In turn, the latter will be true if there exists a scalar $\rho_0 > 0$ such that $H$ is a self-map from the closed Euclidean ball $\rho_0 B$ in $\mathbb{R}^n$ with center at the origin and radius $\rho_0$ into itself. Since

$$\|H(x^{\text{ref}})\| \le \|(M+N)^{-1}b\| + \|(M+N)^{-1}N\| \|\Gamma^{N_h+1}(x^{\text{ref}}) - \Gamma^0(x^{\text{ref}})\| \\ \le \|(M+N)^{-1}b\| + \|(M+N)^{-1}N\| (\rho_{N_h+1} - \rho_0),$$

where the second inequality follows from the recursion (9.6), it follows that if $\rho_0$ satisfies

$$\|(M+N)^{-1}b\| + \|(M+N)^{-1}N\| (\rho_{N_h+1} - \rho_0) < \rho_0, \quad (9.9)$$

then $\|H(x^{\text{ref}})\| < \rho_0$ for every $x^{\text{ref}}$ satisfying $\|x^{\text{ref}}\| \le \rho_0$; i.e., $H$ is a self-map as described. Since $\rho_{N_h+1} \le e^T \psi_x^{\text{IVP}_1} \rho_0 + e^T \psi_x^{\text{IVP}_1} - 1$, (9.9) holds if

$$\|(M+N)^{-1}b\| + \|(M+N)^{-1}N\| (e^T \psi_x^{\text{IVP}_1} - 1)(\rho_0 + 1) < \rho_0.$$

Under condition (9.5), we have $1 - \|(M+N)^{-1}N\| (e^T \psi_x^{\text{IVP}_1} - 1) > 0$; therefore, with the choice of

$$\rho_0 > \frac{\|(M+N)^{-1}b\| + \|(M+N)^{-1}N\| (e^T \psi_x^{\text{IVP}_1} - 1)}{1 - \|(M+N)^{-1}N\| (e^T \psi_x^{\text{IVP}_1} - 1)}, \quad (9.10)$$

(9.9) follows. $\square$

An important remark about the above proof is that the choice of $\rho_0$ is independent of $h$. Based on this observation, we can combine Theorem 7.4 and Lemma 9.1 to obtain our first result for the boundary-value DVI (9.1). In stating this result, we construct from the discrete-time boundary-value trajectory (9.2) a continuous-time trajectory similar to that in Sect. 7.

**Theorem 9.1** Let $K \subseteq \mathbb{R}^m$ be a nonempty closed convex set. Let $(f, B, G)$ satisfy conditions (A), (B), and (D). Let $F(u) = Du$ for some positive semidefinite matrix $D$. Assume further that $M+N$ is nonsingular and (9.5) holds. Under either one of the following conditions:

(i) $K$ contains the origin and $(K, D)$ is an $R_0$ pair,
(ii) K is a polyhedron containing the origin and $G(\Omega) \subseteq \text{int } K(K, D)^*$,

the following three statements hold:

(a) there exists $\bar{h} > 0$ such that for all $h \in (0, \bar{h}]$, a discrete-time boundary-value trajectory (9.2) exists that satisfies (9.3) and (9.4);

(b) there is a sequence $\{h_v\} \downarrow 0$ such that the following limits exist: $\tilde{x}^{h_v} \to \hat{x}$ uniformly on $[0, T]$ and $\hat{u}^{h_v} \to \hat{u}$ weakly in $L^2(0, T)$;

(c) all such limits $(\tilde{x}, \hat{u})$ from part (b) are weak solutions of the boundary-value DVI (9.1).

*Proof* Let $\rho_0 > 0$ satisfy (9.10). Assertion (a) follows from Lemma 9.1. Since $\{x^{h,0}\}$ is bounded in norm by $\rho_0$ for all $h > 0$ sufficiently small, we have the bounds $\|x^{h,i}\| \le \rho_i$ and $\|x^{h,i+1} - x^{h,i}\| \le \rho_{i+1} - \rho_i$ by (9.6) with $\rho_i$ given by (9.7). The uniform boundedness of $\|u^{h,i}\|$ follows from the linear growth property. Part (b) and (c) both follow readily from Theorem 7.1. $\square$

## 9.1 An implicit-explicit scheme for BVPs

Ideally, it would be desirable to relax/remove the strong monotonicity condition (D). It turns out that this is not an easy task because there are many technical challenges that need to be overcome. To mention one such challenge, consider the case where the maps $\tilde{x}^i(x^{\text{ref}})$ defined in Lemma 9.1 are multi-valued; thus so are the composite maps $\Gamma^{i+1}$. Most importantly, the final boundary map $H$ is set-valued. Hence, we need to apply a fixed-point theorem for set-valued maps. For this purpose, we need to establish some suitable topological properties of the images $H(x^{\text{ref}})$ in order to ensure the applicability of such a theorem. This is where the difficulty begins. In fact, we are so far not able to analyze the iteration (9.3)-(9.4) without condition (D), except in the case where the matrix $M$ is nonsingular. (See also the case of the LCS in Sect. 9.2.) Since the case of a nonsingular $M$ is somewhat restrictive in applications, we consider in what follows a variation of the above iteration and establish the convergence of this alternative method by applying Górniewicz' fixed-point Theorem 7.3.

Derived from (3.4), the alternative implicit-explicit time-stepping scheme is as follows: for $i = 0, 1, \dots, N_h$,

$$ 
\begin{align*}
x^{h,i+1} &= x^{h,i} + h[f(t_{h,i+1}, \theta x^{h,i} + (1-\theta)x^{h,i+1}) + B(t_{h,i}, x^{h,i})u^{h,i+1}], \\
u^{h,i+1} &\in \text{SOL}(K, G(t_{h,i}, x^{h,i}) + F), \\
b &= Mx^{h,0} + Nx^{h,N_h+1}.
\end{align*} 
$$

Analytically, a major simplification results from the modified scheme. Namely, corresponding to a given pair $(t_{h,i}, x^{h,i})$, $\text{SOL}(K, G(t_{h,i}, x^{h,i}) + F)$ is a fixed set whose properties depend only on the triple $(K, G, F)$, and in particular, are independent of the functions $f$ and $B$. As a matter of fact, we can totally bypass
the treatment of the algebraic variable $u$, and write the scheme very simply as

$$ 
\begin{aligned}
x^{h,i+1} &\in \tilde{x}^i(x^{h,i}), \quad i = 0, 1, \dots, N_h, \\
b &= Mx^{h,0} + N x^{h,N_h+1},
\end{aligned} 
\quad (9.11) 
$$

where $\tilde{x}^i : \mathbb{R}^n \to \mathbb{R}^n$ is the set-valued map defined by: for $x^{\text{ref}} \in \mathbb{R}^n$,

$$ 
\tilde{x}^i(x^{\text{ref}}) \equiv [I - hf(t_{h,i+1}, \theta x^{\text{ref}} + (1-\theta)\cdot)]^{-1} \times [x^{\text{ref}} + hB(t_{h,i}, x^{\text{ref}})\text{SOL}(K, G(t_{h,i}, x^{\text{ref}}) + F)]. 
$$

We have noted several times that the map $[I - hf(t_{h,i+1}, \theta x^{\text{ref}} + (1-\theta)\cdot)]^{-1}$ is a homeomorphism for all $h > 0$ sufficiently small, provided that $f$ is Lipschitz continuous. Consequently, for such an $h$, the map $\tilde{x}^i$ is well defined. By way of relating $\tilde{x}^i$ to the map $S^i$ defined in Proposition 7.1, we note that $\tilde{x}^i(x^{h,i}) = S^i(x^{h,i})$. The two maps $\tilde{x}^i$ and $S^i$ are different, however; they happen to have the same value at the vector $x^{h,i}$, which is employed in defining the map $S^i$.

In the form (9.11), the alternative time-stepping scheme becomes a DI-based method for treating a BVP, which is different from (3.2) that is derived the DI formulation (2.5). Not aware of a readily available convergence result in the literature that is applicable to (9.11), we give below a detailed treatment of the convergence of the latter scheme. Most importantly, we need to show that for all $h > 0$ sufficiently small, the iterates $\{x^{h,0}, x^{h,1}, \dots, x^{h,N_h+1}\}$ are well defined. We do this by applying Górniewicz' fixed-point Theorem 7.3. Define $j : \mathbb{R}^{2n} \to \mathbb{R}^n$ by

$$ j(x, y) \equiv (M + N)^{-1}b - (M + N)^{-1}N(y - x), \quad (x, y) \in \mathbb{R}^{2n}, $$

and, for a given $h > 0$, also the set-valued map $\Phi : \mathbb{R}^n \to \mathbb{R}^n$ by

$$ \Phi(x^{\text{ref}}) \equiv j(x^{\text{ref}}, \tilde{x}^{N_h} \circ \tilde{x}^{N_h-1} \circ \dots \circ \tilde{x}^1 \circ \tilde{x}^0(x^{\text{ref}})), \quad x^{\text{ref}} \in \mathbb{R}^n. \quad (9.12) $$

Clearly, a fixed point $x^{h,0}$ of $\Phi$ yields the desired iterates via the iterative set membership: $x^{h,i+1} \in \tilde{x}^i(x^{h,i})$, for $i = 0, 1, \dots, N_h$.

Although we could formulate the following result for more general boundary-value DIS, we refrain from doing so and restrict our attention to the DVI on hand. In the theorem below, each function $\tilde{x}^h$ is constructed from the iterates $x^{h,i}$ in the same way as before, i.e., by interpolation.

**Theorem 9.2** Let $K \subseteq \mathbb{R}^m$ be a nonempty closed convex set. Let $(f, B, G)$ satisfy conditions (A) and (B). Let $F : \mathbb{R}^m \to \mathbb{R}^m$ be continuous. Assume that SOL$(K, q + F)$ is a nonempty convex set for all $q \in G(\Omega)$ and that for some constant $\rho > 0$, elements in SOL$(K, q + F)$ satisfy the linear growth property (6.5) for all $q \in G(\Omega)$. Assume further that $M + N$ is nonsingular and (9.5) holds. Then the following two statements are valid:
(a) there exists $\bar{h} > 0$ such that for all $h \in (0, \bar{h}]$, a discrete-time boundary-value trajectory $\{x^{h,0}, x^{h,1}, \dots, x^{h,N_h+1}\}$ exists satisfying (9.11);

(b) there is a sequence $\{h_v\} \downarrow 0$ such that $\hat{x}^{h_v} \to \hat{x}$ uniformly on $[0, T]$; moreover, any such limit $\hat{x}$ along with a suitable $\hat{u}(t) \in \text{SOL}(K, G(t, x(t)) + F)$ constitutes a weak solution of the boundary-value DVI (9.1).

*Proof* As in the proof of Proposition 7.1, it can be shown that $\tilde{x}^i$ is an acyclic map. Proceeding as in Lemma 9.1 and choosing $\rho_0$ as described therein, we can show that $\Phi$ in (9.12) is a set-valued map from $\rho_0\mathcal{B}$ into subsets of $\rho_0\mathcal{B}$. This map $\Phi$, however, is not acyclic in general. But it is sufficient to show that it is an admissible map in the sense of Górniewicz. We can do this by identifying $\Phi$ as a composition of acyclic maps as follows: Let $\Delta : \rho_0\mathcal{B} \to \rho_0\mathcal{B} \times \rho_0\mathcal{B}$ be the diagonal map $\Delta(z) = (z, z)$. Then $\Phi$ can be identified as the composition

$$ \rho_0\mathcal{B} \xrightarrow{\Delta} \rho_0\mathcal{B} \times \rho_0\mathcal{B} \xrightarrow{\text{id} \times \tilde{x}^0} \rho_0\mathcal{B} \times \rho_1\mathcal{B} \xrightarrow{\text{id} \times \tilde{x}^1} \dots \xrightarrow{\text{id} \times \tilde{x}^{N_h}} \rho_0\mathcal{B} \times \rho_{N_h+1}\mathcal{B} \xrightarrow{j} \mathfrak{R}^n. $$

Each “factor” in this composition is an acyclic map, because each factor map is either a continuous single-valued function, or is a Cartesian product of an acyclic map with the identity map. Thus $\Phi$ is an admissible map. Since $\rho_0\mathcal{B}$ is a compact AR, $\Phi$ therefore has a fixed point by Górniewicz’ fixed-point theorem. The proof of the theorem can be easily completed as before. $\square$

## 9.2 The generalized LCS

In this subsection, we develop several specialized results for the boundary-value generalized LCS of the form:

$$ \begin{aligned} \dot{x} &= p + Ax + Bu, \\ C &\ni u \perp q + Cx + Du \in C^*, \\ b &= Mx(0) + Nx(T), \end{aligned} \quad (9.13) $$

where $C$ is a polyhedral cone and the data: $(p, q, A, B, C, D, b, M, N)$ are constant vectors and matrices. We focus on the following time-stepping scheme (same as (9.3) and (9.4)): for $h > 0$ and $i = 0, 1, \dots, N_h + 1$,

$$ \begin{aligned} x^{h,i+1} - x^{h,i} &= h[p + \theta Ax^{h,i} + (1 - \theta)Ax^{h,i+1} + Bu^{h,i+1}], \\ C &\ni u^{h,i+1} \perp q + Cx^{h,i+1} + Du^{h,i+1} \in C^*, \\ b &= Mx^{h,0} + Nx^{h,N_h+1}. \end{aligned} \quad (9.14) $$

Our development herein is in two independent directions: (i) to relax condition (D) which translates into the positive definiteness of the matrix $CB$, and (ii) to retain condition (D) but to consider the case where $D$ is the zero matrix. There are three main results in this subsection, Theorems 9.3, 9.4, and 9.5: the first one accommodates a nonzero $D$, the second one requires $CB$ to be positive definite but has $D = 0$, and the third one pertains to the case where $D = 0$ and $B$ is
the transpose of C. Needless to say, the linear structure in (9.13) provides the
cornerstone for these strengthened results.

We begin by considering the relaxation of condition (D) with a nonzero
matrix D. For h > 0 sufficiently small, we can use the finite difference equation
to solve for x^h[i+1] in terms of x^h[i] and u^h[i+1], obtaining

$$
x^{h,i+1} = [I - h(1-\theta)A]^{-1}[hp + (I + h\theta A)x^{h,i} + hBu^{h,i+1}];
$$

substituting this into the CP yields

$$
C \ni u^{h,i+1} \perp q^h + D^h u^{h,i+1} \in C^*,
$$

where $q^h \equiv q + C[I - h(1 - \theta)A]^{-1}[hp + (I + h\theta A)x^{h,i}]$ and $D^h \equiv D + hC[I - h(1 - \theta)A]^{-1}B$. The matrix $D^h$ is called the transfer matrix in linear control theory. The fundamental role of this matrix in the LCS is well explained in [21,49]. For our purpose, we postulate an assumption on the tuple $(A, B, C, D)$ that ensures some desirable properties of $D^h$.

(ABCD) ≡ (ABCD-1) + (ABCD-2(a or b)) + (ABCD-3), where

(ABCD-1) the matrix D is positive semidefinite;

(ABCD-2a) the first nonzero matrix in the infinite family

$$
\{ CB, (1-\theta) CAB, (1-\theta)^2 CA^2 B, (1-\theta)^3 CA^3 B, \dots \} \quad (9.15)
$$

is positive definite (to avoid triviality, we assume that the family
contains at least one nonzero member);

(ABCD-2b) every matrix in the family (9.15) is positive semidefinite;

(ABCD-3) the implication below holds (0º is defined to be 1):

$$
\begin{align*}
& [u^T Du = u^T C[(1-\theta)A]^i Bu = 0, \forall i = 0, 1, 2, \dots] \\
& \Rightarrow Bu = 0.
\end{align*}
$$

At first sight, condition (ABCD-2b) seems fairly strong because it requires
in particular that CAIB be positive semidefinite for all i (when θ < 1); never-
theless, unlike (ABCD-2a), condition (ABCD-2b) does not require any of the
CAIB to be positive definite. Notice that if condition (ABCD-2a) holds, then
(ABCD-3) holds trivially because the only vector u satisfying the left-hand side
of the implication is the zero vector. In the example below, we illustrate the
satisfaction of condition (ABCD) for a class of LCSs derived from an ODE
with a discontinuous right-hand side, for which the matrix D has not positive
definite for any scalar h.

*Example 6.1 (cont.)* Consider the LCS with data (4.2). Let $-E$ be any psd-plus
matrix. Then $CA^iB = 0$ for all $i \ge 1$ and

$$
CB = \begin{bmatrix} -E & 0 \\ 0 & 0 \end{bmatrix}.
$$
If $u^T C B u = 0$, where $u = (v, x^{-})$, then $Ev = 0$ and thus $Bu = 0$. Therefore condition (ABCD) holds. Since $D^h$ has a zero diagonal block for any scalar $h$, $D^h$ can never be positive definite.

Condition (ABCD) yields the following consequence.

**Lemma 9.2** Let $C$ be a polyhedral cone in $\mathbb{R}^m$. Under assumption (ABCD), the matrix $D^h$ is positive semidefinite for all $h > 0$ sufficiently small; moreover, for all $r \in K(C, D^h)^*$, $Bu$ is constant for all $u \in SOL(C, r, D^h)$; hence the map $B_r : r \in K(C, D^h)^* \mapsto BSOL(C, r, D^h)$ is piecewise linear. If in addition $D$ is psd-plus, then so is $D^h$ for all $h > 0$ sufficiently small.

*Proof* For all $h > 0$ sufficiently small, we have

$$D^h = D + \sum_{i=0}^{\infty} h^{i+1} (1-\theta)^i CA^i B,$$

which establishes the positive semidefiniteness of $D^h$. Hence, $SOL(C, r, D^h)$ is nonempty if and only if $r \in K(C, D^h)^*$. Once we have demonstrated the single-valuedness of the map $B_r$, the claimed properties of this map follows from Lemma 5.1, by the polyhedrality, and thus convexity, of $K(C, D^h)^*$. Take any two solutions $u$ and $u'$ of the CP $(C, r, D^h)$, we must have $(u-u')^T D^h (u-u') = 0$. If condition (ABCD-2a) holds, then $D^h$ is positive definite for all $h > 0$ sufficiently small. Hence $u = u'$. If (ABCD-2b) holds, then $\tilde{u}^T D \tilde{u} = \tilde{u}^T C[(1-\theta)A]^i B \tilde{u} = 0$ for all nonnegative integers $i$, where $\tilde{u} = u - u'$. Consequently, $Bu = Bu'$ by (ABCD-3). Finally, suppose $D$ is psd-plus. If $u^T D^h u = 0$, then $u^T D u = u^T C[(1-\theta)A]^i Bu = 0$ for all nonnegative integers $i$. We then have $Du = 0$ and $Bu = 0$, which then implies $D^h u = 0$. Hence $D^h$ is psd-plus. $\square$

Based on the above lemma, we have the following convergence theorem for the time-stepping scheme (9.14) for solving the boundary-value generalized LCS; this result supplements Theorems 9.1 and 9.2.

**Theorem 9.3** Let $C$ be a polyhedral cone in $\mathbb{R}^m$. Assume condition (ABCD) and that $q + C\mathbb{R}^n \subseteq \text{int } K(C, D)^*$. Assume further that $M+N$ is nonsingular and that (9.5) holds. A discrete-time BV trajectory satisfying (9.14) exists and is uniformly bounded; moreover, a subsequence of this trajectory converges to a weak solution of the boundary-value generalized LCS (9.13).

*Proof* We can follow the proof of Lemma 9.1 and its continuation in Theorem 9.1. The only thing to note is that the map $\tilde{x}^i(x^{\text{ref}})$ remains single-valued. This is because if $(x^h, u^h)$ is any pair, which must exist, satisfying

$$x^h - x^{\text{ref}} = h[p + \theta Ax^{\text{ref}} + (1-\theta)Ax^h + Bu^h], \\ C \ni u^h \perp q + Cx^h + Du^h \in C^*,$$
then $u^h$ must be a solution of the CP $(C, q^h, D^h)$, where

$$q^h = q + C[I - h(1 - \theta)A]^{-1}[hp + (I + h\theta A)x^{\text{ref}}].$$

By Lemma 9.2, it follows that $Bu^h$ is a single-valued Lipschitz continuous function of $x^{\text{ref}}$; hence so is $x^h$, which is $\tilde{x}^i(x^{\text{ref}})$. The rest of the proof is the same as before. □

Parallel to Corollary 8.1, we can state a convergence result for the BV generalized LCS (9.13) with $D$ identically equal to zero. We wish to relax the condition $q + C\mathbb{R}^n \subseteq \text{int } K(C, D)^*$, which becomes $q + C\mathbb{R}^n \subseteq \text{int } C^*$ in this case and is very restrictive. Without the latter condition, however, we can no longer guarantee the linear growth of the solutions to the discrete-time generalized LCPs; therefore, Lemma 7.3 is not applicable. Instead, we need to rely on the proof of Theorem 8.1, which shows that, if CB is positive definite, positive constants $\psi_x^{\text{IVP}_2}$, $\rho_u^{\text{IVP}_2}$, and $\tilde{h}$ exist such that for all $h \in (0, \tilde{h}]$ and all $x^{h,0} \in C^{-1}(C^*)$, a unique pair $(x^{h,i+1}, u^{h,i+1})$ exists satisfying, for all $i = 0, 1, \dots, N_h$,

$$ \begin{aligned} x^{h,i+1} - x^{h,i} &= h[p + \theta Ax^{h,i} + (1 - \theta)Ax^{h,i+1} + Bu^{h,i+1}], \\ C \ni u^{h,i+1} &\perp q + Cx^{h,i+1} \in C^*, \\ b &= Mx^{h,0} + Nx^{h,N_h+1}, \end{aligned} \quad (9.16) $$

$\|x^{h,i+1} - x^{h,i}\| \le h\psi_x^{\text{IVP}_2}(1 + \|x^{h,i}\|)$ and $\|u^{h,i+1}\| \le \rho_u^{\text{IVP}_2}(1 + 2\|x^{h,i}\|)$. Based on these preliminary remarks, we can establish the following result.

**Theorem 9.4** Let $C$ be a polyhedral cone in $\mathbb{R}^m$. Assume that CB is positive definite. Assume further that $M + N$ is nonsingular, (9.5) holds, and $C(M + N)^{-1}(b + N\mathbb{R}^n) \subseteq C^*$. A discrete-time BVP trajectory satisfying (9.16) exists and is uniformly bounded; moreover, a subsequence of this trajectory converges to a weak solution of the boundary-value generalized LCS (9.13) with $D = 0$.

*Proof* We proceed as in the proof of Lemma 9.1. Consider the map $H$ in (9.8). This is a single-valued, continuous map, which we wish to show is a self-map from $S = \rho_0 B \cap C^{-1}(C^*)$ into itself, where $\rho_0$ is the constant satisfying (9.10) corresponding to the constant $\psi_x^{\text{IVP}_2}$ noted above. For a vector $x^{\text{ref}} \in S$, the proof of Lemma 9.1 shows that $H(x^{\text{ref}}) \in \rho_0 B$. Since $H(x^{\text{ref}}) \in (M + N)^{-1}(b + N\mathbb{R}^n)$, the assumption implies that $H(x^{\text{ref}})$ is an element of $C^{-1}(C^*)$. Since $S$ is clearly a compact convex set, it follows that $H$ has a fixed point. The rest of the proof is now familiar. □

*Remark* With $(M, N) = (I, 0)$ and $b = x^0$, the condition $C(M + N)^{-1}(b + N\mathbb{R}^n) \subseteq C^*$ becomes $Cx^0 \in C^*$, which is exactly the special case of the familiar condition $x^0 \in G(t_0, \cdot)^{-1}(C^*)$ in Corollary 8.1 that corresponds to an initial-value generalized LCS with $D = 0$.
Next we consider the following class of boundary-value generalized LCSs:

$$
\begin{align}
\dot{x} &= p + Ax + C^T u, \nonumber \\
C &\ni u \perp q + Cx \in C^*, \tag{9.17} \\
b &= Mx(0) + Nx(T), \nonumber
\end{align}
$$

where $\mathcal{C}$ is a polyhedral cone. It should be noted that no assumption is made on the matrix $A$ throughout the discussion below. As we see from (2.10), with $\mathcal{C}$ being the nonnegative orthant, problem (9.17) provides an equivalent formulation for the VIE: $\dot{x} \in p + Ax - \mathcal{N}(x; K)$, where $K$ is the polyhedron $\{x \in \mathbb{R}^n : q + Cx \ge 0\}$. The latter problem provides an important motivation to consider (9.17).

It is easy to see that for all $h > 0$ sufficiently small, the matrix

$$
D^h = C[I - h(1 - \theta)A]^{-1}C^T
$$

is psd-plus (recall Lemma 6.4) because the matrix within the square bracket is positive definite. The special structure of $D^h$ allows us to establish the following lemma, which is analogous to Proposition 8.3 but without assuming the positive definiteness of the matrix $CC^T$. This lemma is key to a similar convergence result for the problem (9.17) that will follow.

**Lemma 9.3** Let $\mathcal{C}$ be a polyhedral cone. Suppose that $q \in \mathcal{C}\mathbb{R}^n + \mathcal{C}^*$. There exist constants $\bar{h} > 0$ and $\xi > 0$ such that for all $h \in (0, \bar{h}]$ and all $x^{\text{ref}} \in \mathbb{R}^n$, a pair $(x^h, u^h)$ exists satisfying

$$
\begin{equation}
\begin{aligned}
x^h - x^{\text{ref}} &= h[p + \theta Ax^{\text{ref}} + (1-\theta)Ax^h + C^T u^h], \\
C \ni u^h \perp q + Cx^h \in \mathcal{C}^*,
\end{aligned}
\tag{9.18}
\end{equation}
$$

and $\max(\|x^h - x^{\text{ref}}\|, \|u^h\|) \le \xi(1 + \|x^{\text{ref}}\|)$. Moreover, $C^T u^h$ is a (single-valued) Lipschitz function of $x^{\text{ref}}$ with a constant that is independent of $h > 0$ sufficiently small.

*Proof* Let $h > 0$ be sufficiently small so that $I - h(1 - \theta)A$ is positive definite. Since $(C, D^h)$ is an affine pair and $q$ belongs to $C\mathbb{R}^n + \mathcal{C}^*$, it follows that for all $x^{\text{ref}} \in \mathbb{R}^n$, the generalized LCP: $C \ni u^h \perp q^h + D^h u^h \in \mathcal{C}^*$, where $q^h = q + C[I - h(1 - \theta)A]^{-1}[hp + (I + h\theta A)x^{\text{ref}}]$ has a solution. To establish the existence of the constant $\xi$, write $q = Cs + v$ where $s \in \mathbb{R}^n$ and $v \in \mathcal{C}^*$. We have, for any $u^h \in \text{SOL}(C, q^h, D^h)$,

$$
0 \ge (C^T u^h)^T \left\{ s + [I - h(1 - \theta)A]^{-1} [hp + (I + h\theta A)x^{\text{ref}}] \right\} + \sigma \|C^T u^h\|^2,
$$

where $\sigma$ is any positive constant not exceeding the smallest eigenvalue of the symmetric part of the positive definite matrix $[I - h(1 - \theta)A]^{-1}$, for all $h > 0$