Jong-Shi Pang
email: pangj@rpi.edu

David E Stewart
email: dstewart@math.uiowa.edu

Differential variational inequalities

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

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 [START_REF] Brenan | Numerical Solution of Initial-Value Problems in Differential Algebraic Equations[END_REF][START_REF] Petzold | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF] 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 [START_REF] Anitescu | A constraint-stabilized time-stepping for multi-body dynamics with contact and friction[END_REF][START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF][START_REF] Anitescu | Time-stepping for three-dimensional rigid-body dynamics[END_REF][START_REF] Anitescu | A time-stepping method for stiff multi-body dynamics with friction and contact[END_REF][START_REF] Brogliato | Nonsmooth Mechanics. Models, Dynamics and Control[END_REF][START_REF] Brogliato | Numerical simulation of finite dimensional multibody nonsmooth mechanical systems[END_REF][START_REF] Song | Analysis of rigid-body dynamic models for simulation of systems with frictional contacts[END_REF][START_REF] Song | Semi-implicit time-stepping models for frictional compliant contact problems[END_REF][START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of Painlevé's problem[END_REF][START_REF] Stewart | Rigid-body dynamics with friction and impact[END_REF][START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction[END_REF]101,102] and on hybrid engineering systems [START_REF] Çamlibel | Complementarity methods in the analysis of piecewise linear dynamical systems[END_REF][START_REF] Heemels | Linear complementarity systems: a study in hybrid dynamics[END_REF][START_REF] Heemels | Well-posedness of linear complementarity systems[END_REF][START_REF] Heemels | Linear complementarity systems[END_REF], 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 [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF], 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 multivalued differential equations) [START_REF] Aubin | Differential Inclusions[END_REF][START_REF] Deimling | Multivalued Differential Equations[END_REF][START_REF] Monteiro Marques | Differential Inclusions in Nonsmooth Mechanical Problems. Shocks and Dry Friction[END_REF][START_REF] Smirnov | Introduction to the Theory of Differential Inclusions[END_REF] and finite-dimensional variational inequalities (VIs) and complementarity problems (CPs) [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF].

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 [START_REF] Robinson | Generalized equations and their solutions[END_REF][START_REF] Robinson | Strongly regular generalized equations[END_REF][START_REF] Robinson | Generalized equations and their solutions. II. Applications to nonlinear programming[END_REF]; 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 [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]. To this end, we should mention the second author's paper [START_REF] Stewart | A high accuracy method for solving ODEs with discontinuous right-hand side[END_REF], 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.

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.

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 Φ : m → m and a nonempty closed convex set K in m , the VI (K, Φ), is to find a vector u ∈ K such that

( u -u ) T Φ(u) ≥ 0, ∀u ∈ K.
Let SOL(K, Φ) denote the solution set of this problem. If Φ(u) ≡ r + Du is an affine function for some vector r ∈ m and matrix D ∈ m×m , we write VI (K, r, D) for VI (K, Φ) and SOL(K, r, D) for SOL(K, Φ). When K is a polyhedron, we use the adjective "linearly constrained" to describe the VI (K, Φ); we use the acronym AVI to mean the affine variational inequality when (K, Φ) is an affine pair; i.e., K is polyhedral and Φ is affine. We refer to [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] for a comprehensive study of the finite-dimensional variational inequality and the related complementarity problem (CP), where K = C is a cone:

C u ⊥ Φ(u) ∈ C * , (2.1)
with C * ≡ {v ∈ m : u T v ≥ 0, ∀u ∈ C} being the dual cone of C. When C is the nonnegative orthant m + , (2.1) becomes the nonlinear complementarity problem (NCP), which further reduces to the linear complementarity problem (LCP) [START_REF] Cottle | The Linear Complementarity Problem[END_REF] when Φ(u) is additionally an affine function. When Φ is a gradient map, say Φ ≡ ∇φ, where φ is a real-valued function, then the VI (K, Φ) is the first-order stationarity condition for the optimization problem: minimize φ(u) subject to u ∈ K. Standard notation and known results from the LCP literature [START_REF] Cottle | The Linear Complementarity Problem[END_REF] will be employed; in particular, we write LCP (r, D) and SOL(r, D) for VI ( m + , r, D) and SOL( m + , r, D), respectively. Finally, we write ẋ ≡ dx dt for the time-derivative of a function x(t).

Formal definition

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

ẋ(t) = f (t, x(t), u(t)) u(t) ∈ SOL(K, F(t, x(t), •)) (2.2) (x(0), x(T)) = 0. (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 u : [0, T] → K, T 0 ( u(t)u(t) ) T F(t, x(t), u(t)) dt ≥ 0.

(2.4)

The latter implies that for almost all t, u(t) ∈ SOL(K, F(t, x(t), •)). 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 u. Expanding the set of functions u under which (2.4) holds does not change the set of solutions.

Consider the situation where t → F(t, x(t), u(t)) is in L p ([0, T], m ) for some p ≥ 1, and let p satisfy 1/p + 1/p = 1. Since the space of continuous functions C([0, T]; m ) is dense in L p ([0, T]; m ), we can approximate any u in the latter space by continuous w ε in the former space; we can then set u ε (t) ≡ K (w ε (t)), where K is the nearest-point map to K (i.e., the Euclidean projection onto K), in order to ensure that u ε (t) ∈ K for all t. Thus if (2.4) holds for all continuous u, it also holds for all u ∈ L p . In fact, it is not necessary to require (2.4) to hold for all continuous u. We could require only, for example, that (2.4) holds for u ∈ C ∞ . Every continuous u can be approximated by convolutions ψ ε * u (if necessary, extend u beyond [0, T] with u(t) = u(0) for t < 0 and u(t) = u(T) for t > T) for ψ ε (s) = ε -1 ψ(s/ ) with ψ ≥ 0 in C ∞ , and ψ(s) ds = 1. Since K is convex, u(t) ∈ K for all t implies u ε (t) ∈ K for all t. By uniform continuity of u on [0, T], we can show that u ε → u in the uniform norm. Hence if (2.4) holds for u that are C ∞ , it holds for all continuous u, and even holds for all u ∈ L p .

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

x(t) = x(s) + t s f (τ , x(τ ), u(τ )) dτ
holding for all 0 ≤ s ≤ t ≤ 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 , 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 (x, y) ≡ xx 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 ∈ SOL(K, F(0, x 0 , •)), 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 δ > 0, a continuously differentiable function x : [0, δ] → n , and a Lipschitz continuous function u : [0, δ] → K such that x(0) = x 0 , u 0 ∈ SOL(K, F(0, x 0 , •)), and (x(t), u(t)) satisfies (2.2) for all t ∈ [0, δ]. The existence of such solutions is treated in Sects. 5.3 and 5.4.

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 [START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF][START_REF] Brenan | Numerical Solution of Initial-Value Problems in Differential Algebraic Equations[END_REF]. 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, •) is a strongly monotone function. Under a certain uniformity of the strong monotonicity and a Lipschitz continuity of F, the map (t, x) → SOL(K, F(t, x, •)) 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 ẋ = 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 SOL(K, F(t, x, •)) might have multiple solutions. In that case, substitution of solutions u into the differential equation would give not an ODE, but a DI:

ẋ ∈ F(t, x) ≡ f (t, x, SOL(K, F(t, x, •)).
(2.5)

If the set-valued map 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) ∈ 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) ≡ 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-δ 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 [START_REF] Heemels | Linear complementarity systems[END_REF]), 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-δ functions for which statements such as "u(t) ≥ 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 [START_REF] Acary | Higher order Moreau's sweeping process: Mathematical formulation and numerical simulation[END_REF].

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 [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]. Based on the natural map of the VI, we have u ∈ SOL(K, F(t, x, •)) if and only if 0 = u -K (u -F(t, x, u)), where K denotes the Euclidean projector onto the closed convex set K. The nondifferentiability of the projector K is well known (see [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]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 [START_REF] Brenan | Numerical Solution of Initial-Value Problems in Differential Algebraic Equations[END_REF] 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 [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] 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 [START_REF] Aubin | Differential Inclusions[END_REF][START_REF] Deimling | Multivalued Differential Equations[END_REF][START_REF] Monteiro Marques | Differential Inclusions in Nonsmooth Mechanical Problems. Shocks and Dry Friction[END_REF][START_REF] Smirnov | Introduction to the Theory of Differential Inclusions[END_REF] 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 ≡ m , (2.2) becomes

ẋ(t) = f (t, x(t), u(t)), 0 = F(t, x(t), u(t)).
(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):

ẋ(t) = f (t, x(t), u(t)), C u(t) ⊥ F(t, x(t), u(t)) ∈ C * . (2.7)
In turn, a proper specialization of the latter DCP yields the linear complementarity system (LCS) studied extensively in [START_REF] Çamlibel | Complementarity methods in the analysis of piecewise linear dynamical systems[END_REF][START_REF] Heemels | Linear complementarity systems: a study in hybrid dynamics[END_REF][START_REF] Heemels | Well-posedness of hybrid systems[END_REF][START_REF] Heemels | Well-posedness of linear complementarity systems[END_REF][START_REF] Heemels | Linear complementarity systems[END_REF][START_REF] Schumacher | Complementarity systems in optimization[END_REF]. Specifically, the LCS is

ẋ = p + Ax + Bu, 0 ≤ u ⊥ q + Cx + Du ≥ 0, (2.8)
which is clearly a special case of (2.7) where

f (t, x, u) ≡ p + Ax + Bu and F(t, x, u) ≡ q + Cx + Du,
and C ≡ 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 × n, n × m, m × n, and m × m, respectively. As in the study of the finite-dimensional LCP [START_REF] Cottle | The Linear Complementarity Problem[END_REF], 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 [START_REF] Çamlibel | Lyapunov stability of linear complementarity systems[END_REF][START_REF] Çamlibel | Conewise linear systems[END_REF][START_REF] Pang | Strongly regular variational systems[END_REF][START_REF] Shen | Linear complementarity systems:Zeno states[END_REF].

Pioneered by Chen and Mangasarian [START_REF] Chen | Smoothing methods for convex inequalities and linear complementarity problems[END_REF][START_REF] Chen | A class of smoothing functions for nonlinear and mixed complementarity problems[END_REF], 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.

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 n and F be a given (single-valued) mapping from 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) ∈ K for all t ∈ [0, T], and (ii) ẋ(t) ∈ -F(x(t)) -N (x(t); K) for almost all t ∈ [0, T], where, as in convex analysis [START_REF] Rockafellar | Convex Analysis[END_REF], N (z; K) denotes the normal cone of K at an element z ∈ K; i.e.,

N (z; K) ≡ { v ∈ n : v T ( y -z ) ≤ 0, ∀y ∈ 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 [START_REF] Cojocaru | Projected dynamics and evolutionary variational inequalities via Hilbert spaces with applications[END_REF], the authors chose the term "evolutionary variational inequality" to mean a VI in the Hilbert space L p ([0, T], q ) that involves no time derivative. This terminology clashes with the time-varying system ẋ(t) ∈ -F(x(t)) -N (x(t); K), that has been called an "evolution variational inequality" in previous literature; see e.g. [START_REF] Goeleven | Stability and instability matrices for linear evolution variational inequalities[END_REF][START_REF] Goeleven | On the stability of stationary solutions of evolution variational inequalities[END_REF]. 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 [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]. One motivation for considering the VIE is that its stationary solutions, i.e., those for which ẋ = 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 [START_REF] Henry | Differential equations with discontinuous right-hand side for planning procedures[END_REF][START_REF] Henry | An existence theorem for a class of differential equations with multivalued righthand side[END_REF] 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 [START_REF] Hipfel | The nonlinear differential complementarity problem[END_REF], who used the terminology DVI for the VIE. In the same year, Dupuis and Nagurney [START_REF] Dupuis | Dynamical systems and variational inequalities[END_REF] published their study on the class of "projected dynamical systems". For a general existence result for such a differential system in Hilbert space, see [START_REF] Cojocaru | Projected dynamical systems on Hilbert spaces[END_REF][START_REF] Cojocaru | Existence of solutions to projected differential equations in Hilbert space[END_REF]. The recent papers [START_REF] Adly | A stability theory for second order nonsmooth dynamical systems with applications to friction problems[END_REF][START_REF] Goeleven | Stability and instability matrices for linear evolution variational inequalities[END_REF][START_REF] Goeleven | Necessary conditions of asymptotic stability for unilateral dynamical systems[END_REF][START_REF] Goeleven | On the stability of stationary solutions of evolution variational inequalities[END_REF] have developed an extensive stability theory for the VIE.

As an alternative dynamic disequilibrium model, we could consider the ODE: ẋ = x -K (x -F(x)) whose "equilibria" (i.e., vectors x e such that 0 = x e -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 [START_REF] Çamlibel | Lyapunov stability of linear complementarity systems[END_REF].

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 -N (z; K) if and only if

K * u ⊥ z ∈ K. Thus the DI ẋ ∈ -F(x(t)) -N (x(t); K) is equivalent to ẋ(t) = -F(x(t)) + u(t), K * ∈ u(t) ⊥ x(t) ∈ K, (2.9)
which is in the form of the DCP (2.7) with f (t, x, u) ≡ -F(x) + u, F(t, x, u) ≡ x, and C ≡ 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 ≡ { x ∈ n : g(x) ≤ 0 },
where g : n → m is continuously differentiable such that each component function g i is convex, for i = 1, . . . , 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 ∈ K, a vector v belongs to N (x; K) if and only if there exist multipliers u

∈ m satisfying v = m i=1 u i ∇g i (x) and 0 ≤ u ⊥ g(x) ≤ 0.
Consequently, the DI ẋ ∈ -F(x(t)) -N (x(t); K) is equivalent to

ẋ(t) = -F(x(t)) - m i=1 u i ∇g i (x), 0 ≤ u ⊥ -g(x) ≥ 0, (2.10)
which can be identified as the DCP (2.7) with f (t, x, u) ≡ -F(x) -Jg(x) T u, C ≡ m + , and F(t, x, u) ≡ -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 [START_REF] Goeleven | Stability and instability matrices for linear evolution variational inequalities[END_REF] 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 ẋ ∈ -F(x(t)) -N (x(t); K) is equivalent to

ẋ(t) = -F(x(t)) + w(t), 0 = x(t) -y(t), y(t) ∈ K, ( v -y(t) ) T w(t) ≥ 0, ∀v ∈ K. (2.11)
This is in the form of the DVI with u ≡ (w, y), f (t, x, u) ≡ -F(x) + w, F(t, x, u) ≡ (xy, w) and K ≡ n × K. It is interesting to note that when 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:

ẋ(t) = f (x(t), y(t), w(t)), 0 = g(x(t), y(t), w(t)), y(t) ∈ K, ( v -y(t) ) T w(t) ≥ 0, ∀v ∈ K, (2.12)
we obtain a unified formulation for both the DVI and the VIE.

Extended systems

There are many generalizations of the VI (see [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]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] × m into closed convex subsets of m : i.e., a system of the kind:

ẋ(t) = f (t, x(t), u(t)), u(t) ∈ SOL(K(t, x(t)), F(t, x(t), •)).
Properly specializing the functions f and F, we obtain the sweeping processes developed by Moreau, Monteiro-Marques, Kunze and others [START_REF] Castaing | Evolution equations governed by the sweeping process[END_REF][START_REF] Kunze | Existence of solutions for degenerate sweeping processes[END_REF][START_REF] Kunze | An introduction to Moreau's sweeping process[END_REF][START_REF] Monteiro Marques | Differential Inclusions in Nonsmooth Mechanical Problems. Shocks and Dry Friction[END_REF][START_REF] Moreau | Evolution problem associated with a moving convex set in a Hilbert space[END_REF][START_REF] Moreau | Bounded variation in time[END_REF][START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Stewart | Reformulations of measure differential inclusions and their closed graph property[END_REF]. These processes can be formulated as DIs of the kind:

ẋ(t) ∈ -F(x(t)) -N (x(t); K(t, x(t))), ( 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(•, •) 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(•) 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 [START_REF] Monteiro Marques | Differential Inclusions in Nonsmooth Mechanical Problems. Shocks and Dry Friction[END_REF][START_REF] Moreau | Evolution problem associated with a moving convex set in a Hilbert space[END_REF][START_REF] Moreau | Bounded variation in time[END_REF][START_REF] Stewart | Reformulations of measure differential inclusions and their closed graph property[END_REF]. Generalizations of (2.13) to nonconvex K(t, x) can be found, for instance, in [START_REF] Bounkhel | General existence results for second order nonconvex sweeping process with unbounded perturbations[END_REF][START_REF] Faik | Differential inclusions governed by a nonconvex sweeping process[END_REF]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 [START_REF] Petrov | A pseudodifferential linear complementarity problem related to a one-dimensional viscoeleastic model with Signorini conditions[END_REF] and further analyzed by the second author [START_REF] Stewart | Convolution complementarity problems with application to impact problems[END_REF]. 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.

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 . The first issue is dealt with by time stepping, i.e., replacing the time derivative ẋ 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 ) 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 < • • • < 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, ∀i = 0, 1, . . . , N h .
We approximate the time derivative ẋ by a backward finite-difference quotient:

ẋ(t h,i+1 ) ≈ x(t h,i+1 ) -x(t h,i ) h , i = 0, 1, . . . , N h .
We employ a scalar θ ∈ [0, 1] to distinguish an explicit (θ = 1), an implicit (θ = 0), or a semi-implicit (θ ∈ (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 , . . . , x h,N h +1 } ⊂ n and {u h,1 , u h,2 , . . . , u h,N h +1 } ⊂ m . (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 ∈ x h,i + h F(t h,i+1 , θ x h,i + (1 -θ)x h,i+1 ), i = 0, 1, . . . , N h .
Unwrapping the set-valued map F, we obtain

x h,i+1 = x h,i + h f (t h,i+1 , θ x h,i + (1 -θ)x h,i+1 , u h,i+1 ) u h,i+1 ∈ SOL(K, F(t h,i+1 , θ x h,i + (1 -θ)x h,i+1 , •)) 0 = (x h,0 , x h,N h +1 ) ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ , ( 3.2)
Dontchev and Lempio [START_REF] Dontchev | Difference methods for differential inclusions: A survey[END_REF] have surveyed various difference methods for initialvalue 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 θ 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:

x h,i+1 = x h,i + h f (t h,i+1 , θ x h,i + (1 -θ)x h,i+1 , u h,i+1 ) u h,i+1 ∈ SOL(K, F(t h,i+1 , x h,i+1 , •)) 0 = (x h,0 , x h,N h +1 ) ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ , (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

x h,i+1 = x h,i + h f (t h,i+1 , θ x h,i + (1 -θ)x h,i+1 , u h,i+1 ) u h,i+1 ∈ SOL(K, F(t h,i+1 , x h,i , •)) 0 = (x h,0 , x h,N h +1 ) ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ , ( 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 h,i } 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 h,i given, the first two conditions in (3.3) define a VI ( K, H h,i+1 ) in the variable (x h,i+1 , u h,i+1 ), where K ≡ n × K and

H h,i+1 (x, u) ≡ x -x h,i -h f (t h,i+1 , θ x h,i + (1 -θ)x, u) F(t h,i+1 , x, u) , ( x, u ) ∈ 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). Namely, 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 differential 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

x h,i+1 = x h,i + t h,i+1 t h,i f (τ , θx h,i + (1 -θ)x h,i+1 ) dτ +h B(t h,i+1 , θ x h,i + (1 -θ)x h,i+1 ) u h,i+1 u h,i+1 ∈ SOL(K, F(t h,i+1 , x h,i+1 , •)) 0 = (x h,0 , x h,N h +1 ) ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ . (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.

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 [START_REF] Çamlibel | Complementarity methods in the analysis of piecewise linear dynamical systems[END_REF][START_REF] Heemels | Linear complementarity systems: a study in hybrid dynamics[END_REF][START_REF] Heemels | Well-posedness of linear complementarity systems[END_REF][START_REF] Heemels | Linear complementarity systems[END_REF][START_REF] Schumacher | Complementarity systems in optimization[END_REF]. The survey paper [START_REF] De | Modeling and simulation of genetic regulatory systems: a literature review[END_REF] 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 applications, such as the dynamic traffic assignment problem [START_REF] Ban | Quasi-variational inequality formulations and solution approaches for dynamic user equilibria[END_REF][START_REF] Liu | An analytical dynamic traffic assignment model with probabilistic travel times and perceptions[END_REF][START_REF] Ran | Modeling Dynamic Transportation Networks[END_REF], multibody frictional contact problems (see [START_REF] Pang | Convergence of time-stepping methods for initial and boundary value frictional compliant contact problems[END_REF][START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of Painlevé's problem[END_REF] 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 modeling of atmospheric aerosols [START_REF] Kakhu | Dynamic modeling of acqueous electrolyte systems[END_REF]75], as well as many applications of the VIE, which has a large body of literature by itself.

ODEs with discontinuous right-hand sides

While pioneering work on ODEs with discontinuous right-hand sides was initiated by Filippov [START_REF] Filippov | Differential equations with discontinuous right-hand side[END_REF][START_REF] Filippov | Differential Equations with Discontinuous Right-Hand Sides[END_REF], the DCP approach to these ODEs originated from the second author's Ph.D. thesis [START_REF] Stewart | A high accuracy method for solving ODEs with discontinuous right-hand side[END_REF], who used linear complementarity theory to study "piecewise active solutions" to ODEs with piecewise smooth (but discontinuous) 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 ∈ n×n and consider:

ẋ = E sgn(x), (4.1)
where sgn(x) be the signum function with components sgn(x i ) with sgn(τ )

⎧ ⎪ ⎨ ⎪ ⎩ = 1 i fτ > 0, ∈ [ -1, 1 ] if τ = 0, = -1 i fτ < 0.
The inclusion sgn(0) ∈ [-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 ± are the nonnegative and nonpositive part of x, respectively. In terms of these two nonnegative (and complementary) variables x ± , we can further write sgn(x) = 1v, where v satisfies

0 ≤ v ⊥ x + ≥ 0 and 0 ≤ 2 -v ⊥ x -≥ 0,
with 1 and 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):

ẋ = E 1 -Ev, 0 ≤ v ⊥ x + x -≥ 0, 0 ≤ x -⊥ 2 -v ≥ 0,
which corresponds to (2.8) with u ≡ (v, x -), and

p ≡ E 1, A ≡ 0, B ≡ -E 0 , q ≡ 0 2 , C ≡ I 0 , D ≡ 0 I -I 0 . (4.2)
Note that D is a positive semidefinite matrix. Moreover, if -E is positive semidefinite, then so is D + τ CB for all τ > 0. These are important properties when it comes to solving the problem.

Differential Nash games

Consider a noncooperative Nash equilibrium problem with N players where each player solves an optimal control problem. Let (x ν , u ν ) ∈ n ν +m ν denote player ν's pair of state and control variables, respectively; let U ν ⊆ m ν be the closed convex set of player ν's admissible controls; let x ≡ (x ν ) N ν=1 ∈ n and u ≡ (u ν ) N ν=1 ∈ m , where

n ≡ N ν=1 n ν and m ≡ N ν=1 m ν
be the collection of all players' variables. Dependent on all players' variables, player ν's cost function is given by

θ ν (x, u) ≡ ψ ν (x(T)) + T 0 ϕ ν (t, x(t), u(t)) dt,
where T denotes the terminal time and ψ ν : n → and ϕ ν : 1+n+m → are continuously differentiable functions. Parameterized by the rival players' strategies x -ν ≡ (x ν ) ν =ν and u -ν ≡ (u ν ) ν =ν , which are taken to be exogenous to player ν's optimization problem but endogenous to the overall model, player ν's optimization problem is to compute an optimal trajectory (x ν (t), u ν (t)), for each fixed but arbitrary (x -ν , u -ν ), in order to minimize θ ν (x, u)

subject to ẋν (t) = g ν (t, x ν (t), u ν (t)), x ν (0) = x 0,ν and u ν (t) ∈ U ν ⊆ m ν , for almost all t ∈ [0, T],
where g ν : 1+n ν +m ν → ν is continuously differentiable. A differential Nash equilibrium solution [START_REF] Başar | Dynamic Noncooperative Game Theory[END_REF][START_REF] Dockner | Differential Games in Economics and Management Science[END_REF] is a pair of trajectories (x * (t), u * (t)) such that for every ν = 1, . . . , N, (x * ,ν (t), u * ,ν (t)) solves player ν's problem, given that ν's rivals all play their Nash strategies (x * ,-ν (t), u * ,-ν (t)). To derive a set of necessary conditions for a differential Nash equilibrium trajectory, let

H ν (t, x, u, λ ν ) ≡ ϕ ν (t, x, u) + (λ ν ) T g ν (t, x ν , u ν ), λ ν ∈ ν ,
be player ν's Hamiltonian function, where λ ν is the adjoint variable of the ODE constraint in player ν's control problem. The first-order necessary optimality condition of player ν's problem is

λν (t) = -∇ x ν H ν (t, x(t), u(t), λ ν (t)), ẋν (t) = g ν (t, x ν (t), u ν (t)), u ν (t) ∈ SOL(U ν , ∇ u ν H ν (t, x(t), u -ν (t), •, λ ν (t)), x ν (0) = x 0,ν and λ ν (T) = ∇ x ν ψ ν (x(T)).
Concatenating these conditions for all players, we obtain the BV-DVI:

λ(t) ẋ(t) = H(t, x, u, λ) g(t, x, u) , u(t) ∈ SOL(U, F(t, x(t), λ(t), •)), 0 = (x(0), x(T), λ(T)),
where U ≡ N ν=1 U ν ⊆ m is the Cartesian product of lower-dimensional sets,

H(t, x, u, λ) ≡ -∇ x ν H ν (t, x, u, λ ν ) N ν=1 , g(t, x, u) ≡ g ν (t, x, u) N ν=1 , F(t, x, λ, u) ≡ ∇ u ν H ν (t, x, u, λ ν ) N ν=1 , (y, z, λ) ≡ y ν -x 0,ν λ ν -∇ z ν ψ ν (z) N ν=1 .
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 ν is replaced by the variable set U ν (x -ν , u -ν ). In this case, the above DVI becomes the DQVI where the variational condition:

u(t) ∈ SOL(U, F(t, x(t), λ(t), •)) is replaced by the quasi-variational condition: u(t) ∈ SOL(U(x(t), u(t)), F(t, x(t), λ(t), •)), where U(x(t), u(t)) ≡ N ν=1 U ν (x -ν , u -ν ).

The linear-quadratic case

The linear-quadratic differential Nash game has three specifications: (a) each

ψ ν (x(T)) is quadratic in x(T): ψ ν (x(T)) ≡ x ν (T) T c ν + ν ∈N W νν x ν (T)
for some constant vector c ν ∈ n ν and matrices W νν ∈ n ν ×n ν with W νν symmetric (note: the terms not involving x ν (T) are not included because they are exogenous to player ν's optimization problem); and (b) each ϕ ν (t, x, u) is quadratic in (x, u): for some time-dependent vectors p ν (t) ∈ n ν and q ν (t) ∈ m ν and time-dependent matrices

P νν (t) ∈ n ν ×n ν , Q νν (t) ∈ n ν ×m ν , R νν (t) ∈ m ν ×n ν , and S νν (t) ∈ m ν ×m ν : ϕ ν (t, x, u) ≡ ( x ν ) T p ν (t) + ν ∈N P νν (t)x ν + Q νν (t)u ν +( u ν ) T q ν (t) + ν ∈N R νν (t)x ν + S νν (t)u ν
and (c) each g ν (t, x ν , u ν ) is an affine function in (x ν , u ν ): for some vector r ν (t) ∈ n ν and matrices G ν (t) ∈ n ν ×n ν and H ν (t) ∈ n ν ×m ν :

g ν (t, x ν , u ν ) ≡ r ν (t) + G ν (t)x ν + H ν (t)u ν ,
Under the above specifications, we have

H(t, x, u, λ) g(t, x, u) = d(t) + A(t) λ x + B(t)u, (4.3)
where

d(t) = ( -p ν (t) ) ν∈N ( r ν (t) ) ν∈N ∈ 2n A(t) = ⎡ ⎢ ⎣ -diag(G ν (t) T ) ν∈N -P νν (t) νν ∈N -diag(P νν (t)) ν∈N 0 d i a g (G ν (t)) ν∈N ⎤ ⎥ ⎦ ∈ 2n×2n B(t) = ⎡ ⎢ ⎣ -Q νν (t) νν ∈N -diag R νν (t) T ν∈N diag(H ν (t)) ν∈N ⎤ ⎥ ⎦ ∈ 2n×m .
Moreover,

F(t, x, λ, u) = q(t) + C(t) λ x + S(t)u, (4.4)
where S(t) ≡ [S νν (t)] νν ∈N + diag(S νν (t)) ν∈N and

C(t) ≡ diag(H ν (t) T ) ν∈N R νν (t) νν ∈N + diag Q νν (t) T ν∈N ∈ m×2n .
Finally, we also have

∇ x ν (T) ψ ν (x(T)) = c ν + W νν x ν (T) + ν ∈N W νν x ν (T).
Rewriting the VI (U, F(t, x, λ, •)) as the triple (U, q(t, x, λ), S(t)), where

q(t, x, λ) ≡ q(t) + C(t) λ x ,
we see that the DNG can be formulated as:

λ ẋ = d(t) + A(t) λ x + B(t)u u(t) ∈ SOL(U, q(t, x(t), λ(t)), S(t)) x(0) = x 0 , λ(T) = c + Wx(T) ⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ , ( 4.5)
where c ≡ (c ν ) ν∈N and W ≡ [W νν ] ν,ν ∈N +diag(W νν ) ν∈N . Notice that in (4.5), the right-hand side in the ODE is an affine function of the state variables (x, λ) and algebraic variable u and that the middle VI is defined by the affine mapping u → q(t, x(t), λ(t)) + S(t)u.

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 [START_REF] Coddington | Theory of Ordinary Differential Equations[END_REF]. 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 ẋ = 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 [START_REF] Pang | Strongly regular variational systems[END_REF] where the issue of "non-Zenoness" of a strongly regular DVI is studied (see Sect. 5.3 for the definition of strong regularity).

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, •)) has a unique solution, say u(t, x), that is a Lipschitz continuous function of (t, x) for all (t, x) ∈ [0, T] × 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 ≡ N ν=1 K ν , ( 5.1)
where each K ν is a closed convex subset of n ν . We postulate the following conditions on the pair (K, F):

(A 0 ) F(t, x,
•) is a continuous, uniformly P function [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] on K with a modulus that is independent of (t, x); i.e., there exists a constant η F > 0 such that

max 1≤ν≤N ( u ν -u ν ) T ( F ν (t, x, u) -F ν (t, x, u ) ) ≥ η F u -u 2 for all (t, x) ∈ [0, T] × X and u ≡ (u ν ) N ν=1 and u ≡ (u ν ) N ν=1 in K = N ν=1 K ν ; (B 0 ) F(•, •, 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) ≤ L F [ | t -t | + x -x ],
for all (t, x) and (t , x ) in [0, T] × X and all u ∈ K.

From well-known theory of partitioned VIs [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]Sect. 3.5], we have Theorem 5.1 Let K be given by (5.1), where each K ν is a closed convex subset of n ν . Under the two assumptions (A 0 ) and (B 0 ), there exists a Lipschitz continuous function u :

[0, T] × X → K such that for each pair (t, x) ∈ [0, T] × X, u(t, x) is the unique solution of the VI (K, F(t, x, •)).
In terms of the implicitly defined function u(t, x), (2.2) is equivalent to the ODE:

ẋ = f (t, x, u(t, x)). ( 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 [START_REF] Ascher | Numerical Solution of Ordinary Value Problems for Ordinary Differential Equations[END_REF][START_REF] Keller | Numerical Methods for Two-Point Boundary-Value Problems[END_REF][START_REF] Roberts | Two-Point Boundary Problems: Shooting Methods[END_REF], which invariably utilize some Newton methods for solving smooth systems of algebraic equations. In a subsequent paper [START_REF] Pang | Solution dependence on initial conditions in differential variational inequalities Math[END_REF], we have shown that the implicit function u(t, x) is "B(ouligand)-differentiable" [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]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 [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]Chaps. 7 and 8]); this is indeed possible for problems of the affine type; for details, see [START_REF] Pang | Solution dependence on initial conditions in differential variational inequalities Math[END_REF].

Differential semi-affine systems

It is possible to relax assumption (A 0 ) 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 [START_REF] Robinson | Normal maps induced by linear transformations[END_REF][START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] 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:

ẋ = f (x) + Bu, x(0) = x 0 u ∈ SOL(K, Cx, D), ( 5
functions { f 1 , . . . , f k } such that Φ(x) ∈ { f 1 (x) . . . , f k (x)} for all x ∈ dom Φ. If dom Φ is convex, then Φ is Lipschitz continuous on dom Φ.
Proposition 5.1 Let f : n → n be Lipschitz and let K ⊆ m be polyhedral. A necessary and sufficient condition for (5.3) to have a unique C 1 solution trajectory x(t) on [0, ∞) for all x 0 ∈ n is that BSOL(K, Cx, D) is a singleton for all x ∈ n .

Proof The necessity is obvious; indeed, for any x 0 ∈ n and any u ∈ SOL (K, Cx 0 , D), we must have Bu = ẋ(0)-f (x 0 ), which is necessarily unique because the trajectory x(t) is uniquely determined by x 0 by assumption. Conversely, if BSOL(K, Cx, D) is a singleton for all x ∈ n , then since the set-valued map x → SOL(K, Cx, D) is a polyhedral multifunction [START_REF] Robinson | Some continuity properties of polyhedral multifunctions[END_REF][START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF], by Lemma 5.1, the single-valued map x → BSOL(K, Cx, D) is Lipschitz continuous. The existence and uniqueness of a C 1 solution trajectory x(t) of (5.3) on [0, ∞) follows from classical ODE theory.

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 ≡ E T GE for some (possibly asymmetric) positive definite matrix G ∈ k×k and an arbitrary matrix E ∈ k×m , and B ≡ HE some matrix H ∈ n×k . For such a matrix D, it is known [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] that ESOL(K, r, D), thus BSOL(K, r, D), is a singleton for all vectors 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 [START_REF] Çamlibel | Lyapunov stability of linear complementarity systems[END_REF][START_REF] Pang | Strongly regular variational systems[END_REF][START_REF] Shen | Linear complementarity systems:Zeno states[END_REF].

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 , •)). 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 , •)) has a strongly regular solution u 0 [START_REF] Robinson | Strongly regular generalized equations[END_REF][START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]. Under (A 0 ) and (B 0 ), it follows that there exist a scalar ε > 0, a neighborhood N of x 0 , a neighborhood U of u 0 , and a Lipschitz continuous function u :

[0, ε] × N → K ∩ U such that for each pair (t, x) ∈ [0, ε] × N , u(t, x) is the unique solution of the VI (K, F(t, x, •)) that belongs to U. Therefore in the domain [0, ε] × 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 [START_REF] Coddington | Theory of Ordinary Differential Equations[END_REF].

It is possible to refine the condition (A 0 ) in the case when the set K is the nonnegative orthant m + so that the VI (K, F(t, x, •)) is equivalent to the nonlinear complementarity problem (NCP):

0 ≤ u ⊥ F(t, x, u) ≥ 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 ≤ u 0 ⊥ F(0, x 0 , u 0 ) ≥ 0.

(5.4)

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

α ≡ { i : u 0 i > 0 = F i (0, x 0 , u 0 ) } the support of u 0 β ≡ { i : u 0 i = 0 = F i (0, x 0 , u 0 ) } the degenerate set of u 0 γ ≡ { i : u 0 i = 0 < F i (0, x 0 , u 0 ) } the inactive set of u 0 .
For any two index subsets α 1 and α 2 of {1, . . . , m}, let J u α 1 F α 2 denote the submatrix of the (partial) Jacobian matrix J u F of F with respect to the u-variable with rows and columns indexed by α 2 and α 1 , respectively; i.e.,

J u α 1 F α 2 ≡ ∂F i ∂u j : ( i, j ) ∈ α 2 × α 1 .
The strong regularity of u 0 is equivalent to two conditions [START_REF] Robinson | Strongly regular generalized equations[END_REF][START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]: (a) the matrix J u α F α (0, x 0 , u 0 ) is nonsingular, and (b) the Schur complement

J u β F β (0, u 0 , u 0 ) -J u α F β (0, x 0 , u 0 ) J u α F α (0, x 0 , u 0 ) -1 J u β F α (0, x 0 , u 0 ) of J u α F α (0, x 0 , u 0 ) in J u α F α (0, x 0 , u 0 ) J u β F α (0, x 0 , u 0 ) J u α F β (0, x 0 , u 0 ) J u β F β (0, x 0 , u 0 )
is a P-matrix (see [START_REF] Cottle | The Linear Complementarity Problem[END_REF] 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 α F α (0, x 0 , u 0 ) is singular (therefore u 0 is not strongly regular), we consider the extreme situation where F α is independent of the variable u and write

Φ α (t, x) ≡ F α (t, x, u).
(The other components of F can still depend on u.) We define the function

F α (t, x, u) ≡ ∂Φ α (t, x) ∂t + J x Φ α (t, x) f (t, x, u),
which is obtained by differentiating Φ α with respect to t and then substituting f (t, x, u) for ẋ (provided that the latter substitution is legal). Using the function F α (t, x, u) instead of the original F α (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):

ẋ = f (t, x, u), x(0) = x 0 , 0 ≤ u ⊥ F(t, x, u) ≥ 0 (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 α (t, x, u) = Φ α (t, x) is independent of u and that F α (0, x 0 , u 0 ) = 0. If the matrix J u α F α (0, x 0 , u 0 ) = J x Φ α (0, x 0 )J u α f (0, x 0 , u 0
) is nonsingular, and the following Schur complement is a P-matrix:

J u β F β (0, x 0 , u 0 ) -J u α F β (0, x 0 , u 0 ) J x Φ α (0, x 0 )J u α f (0, x 0 , u 0 ) -1 J x Φ α (0, x 0 )J u β f (0, x 0 , u 0 )
then the DNCP (5.5) has a locally smooth solution near (0, x 0 , u 0 ).

Proof Consider the modified DNCP:

ẋ = f (t, x, u), x(0) = x 0 , 0 ≤ u ⊥ F(t, x, u) ≥ 0, (5.6)
where the components of F are equal to the corresponding components of F except for those in α which are equal to the ones in F α . 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, δ] for some δ > 0. Moreover, provided that δ > 0 is sufficiently small, u * α (t) > 0 on [0, δ]. Hence it remains to show that F α (t, x * (t), u * (t)) = Φ α (t, x * (t)) = 0 on [0, δ]. We have F α (t, x * (t), u * (t)) = d dt Φ α (t, x * (t)) = 0 for all t ∈ [0, δ]. Since Φ α (0, x * (0)) = 0,
simple integration yields Φ α (t, x * (t)) = 0 on [0, δ], completing the proof.

Linear complementarity systems

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

0 ≤ u ⊥ q + Cx + Du ≥ 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 [START_REF] Çamlibel | Complementarity methods in the analysis of piecewise linear dynamical systems[END_REF][START_REF] Heemels | Linear complementarity systems: a study in hybrid dynamics[END_REF] and the subsequent paper [START_REF] Heemels | Linear complementarity systems[END_REF] 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 [START_REF] Heemels | Linear complementarity systems[END_REF]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

D αα D αβ D βα D ββ
that are linearly dependent on another set of rows in the same matrix, where α is the support of u 0 and β is the degenerate index set of u 0 as a solution of the LCP (r 0 , D). Note that β = ∅ implies in particular that u 0 is not strongly regular. For simplicity, we assume specifically that a nonsingular matrix P αα can be found such that

P αα [ D αα D αβ ] = E α 1 α E α 1 β 0 0 ∈ |α|×(|α|+|β|)
for some subset α 1 (with complement α 2 ) of α, where the two zero matrices are of order α 2 × α and α 2 × β, respectively. (For instance, if D is the zero matrix, then P αα can be taken to be the identity matrix.) Setting u γ = 0, premultiplying the equation 0

= C α• x + D αα u α + D αβ u β by P αα , we obtain 0 = P α 2 α C α• x, (5.7)
in which the variable u has been eliminated. Thus, we can differentiate the latter equation with respect to time and then substitute ẋ using the first equation in (2.8). This process results in the following modified LCS of reduced size:

ẋ = Ax + B •α u α + B •β u β , x(0) = x 0 0 ≤ ⎛ ⎜ ⎝ u α 1 u α 2 u β ⎞ ⎟ ⎠ ⊥ ⎡ ⎢ ⎣ P α 1 α C α• P α 2 α C α• A C β• ⎤ ⎥ ⎦ x + ⎡ ⎢ ⎣ E α 1 α E α 1 β P α 2 α C α• B •α P α 2 α C α• B •β D βα D ββ ⎤ ⎥ ⎦ u α u β ≥ 0.
If this modified LCS has a locally smooth solution (x * (t), u * α (t), u * β (t)) on an interval [0, δ] such that

P α 1 α C α• P α 2 α C α• A x * (t) + E α 1 α E α 1 β P α 2 α C α• B •α P α 2 α C α• B •β u * α (t) u * β (t) = 0 and C γ • x * (t) + D γ α u * α (t) + D γβ u * β (t) ≥ 0 on [0, t], with u * γ ≡ 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 * α (t), u * β (t)).
Theorem 5. [START_REF] Anitescu | A constraint-stabilized time-stepping for multi-body dynamics with contact and friction[END_REF] In the above setting, suppose that P

α 2 α C α• (Ax 0 +Bu 0 ) = 0. Assume further that the matrix E α 1 α P α 2 α C α• B •α (5.8)
is nonsingular and that the Schur complement

D ββ -D βα E α 1 α P α 2 α C α• B •α -1 E α 1 β P α 2 α C α• B •β (5.9)
is a P-matrix. The LCS (2.8) has a locally smooth solution near (0, x 0 , u 0 ).

P α 1 α C α• x * (t) + D αα u * α (t) + D αβ u * β (t) = 0, (5.10
)

P α 2 α C α• Ax * (t) + B •α u * α (t) + B •β u * β (t) = 0 (5.11)
on the domain of the trajectory. From (5.11), we obtain 0 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.

= P α 2 α C α• ẋ * (t) on [0, δ], which implies P α 2 α C α• x * (t) + D αα u * α (t) + D αβ u * β (t) = 0.

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 ≡ [0, T] × n and let B denote the closed Euclidean unit ball centered at the origin. → n be an upper semicontinuous set-valued map with nonempty closed convex values. Suppose that there exists a scalar ρ F > 0 satisfying

sup{ y : y ∈ F(t, x) } ≤ ρ F ( 1 + x ), ∀ ( t, x ) ∈ . (6.1)
For every x 0 ∈ n , the DI: ẋ ∈ 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 F holds if F is closed (that is, if the graph of F is closed). Indeed, this assumption implies that F is bounded on compact sets; therefore the closedness of F is equivalent to its upper semicontinuity.

Setting of the DVI

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

ẋ = f (t, x) + B(t, x)u, x(0) = x 0 , u ∈ SOL(K, G(t, x) + F(•)), ( 6.2)
where K ⊆ m is closed and convex, (f , B, G) :

→ n × n×m × m , and F : m → m are given functions. We let K ∞ be the recession cone of K; see [START_REF] Rockafellar | Convex Analysis[END_REF]. 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 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 with Lipschitz constants L f , L B , and L G , respectively; (B) B is bounded on with σ B ≡ sup

(t,x)∈ B(t, x) < ∞.
Notice that (A) implies that f and G have linear growth on in x; i.e., for some positive constants ρ f and ρ G and for all (t, x) ∈ ,

f (t, x) ≤ ρ f ( 1 + x ) and G(t, x) ≤ ρ G ( 1 + x ). (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) ∈ ,

F(t, x) ≡ { f (t, x) + B(t, x)u : u ∈ SOL(K, G(t, x) + F) }, (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 F(t, x) is a convex set; however, the converse is false. Proposition 6.5 identifies a case where F(t, x) is convex without F being monotone. The lemma below shows that the remaining conditions on the map 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 : m → m be continuous. Suppose that SOL(K, q + F) = ∅ for all q ∈ G( ), and there exists a constant ρ > 0 such that for all q ∈ G( ),

sup{ u : u ∈ SOL(K, q + F) } ≤ ρ ( 1 + q ). (6.5)
Then a constant ρ F > 0 exists such that (6.1) holds for the map F defined by (6.4). Hence, F is upper semicontinuous, and thus closed-valued, on .

Proof It is easy to see that

sup{ y : y ∈ F(t, x) } ≤ [ ρ f + ρ σ B ( 1 + ρ G ) ] ( 1 + x ), on .
Hence F has linear growth. To prove the upper semicontinuity of F, it suffices to show that F is closed. Suppose that {(t k , x k )} ⊂ is a sequence converging to some vector (t ∞ , x ∞ ) in and {f (t k , x k ) + B(t k , x k )u k } converges to some vector z ∞ ∈ n as k → ∞, where u k ∈ 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 ∞ , which is easily seen to be a solution of the limiting VI (K, G(t ∞ , x ∞ )+ F), by a simple limiting argument. Thus,

z ∞ = f (t ∞ , x ∞ )+B(t ∞ , x ∞ )u ∞ belongs to F(t ∞ , x ∞ ), completing the proof.
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 [START_REF] Filippov | On certain questions in the theory of optimal control[END_REF]; see also Theorem 2.3 in [START_REF] Smirnov | Introduction to the Theory of Differential Inclusions[END_REF] and Corollary 1 in [START_REF] Aubin | Differential Inclusions[END_REF]. The reader can consult these references for a proof. Lemma 6.3 Let h : × m → n be a continuous function and U : → m be a closed set-valued map such that for some constant η U > 0,

sup u∈U(t,x) u ≤ η U ( 1 + x ), ∀ ( t, x ) ∈ .

Let v : [0, T] → n be a measurable function and x : [0, T] → n be a continuous function satisfying v(t) ∈ h(t, x(t), U(t, x(t))) for almost all t ∈ [0, T]. There exists a measurable function u : [0, T] → m such that u(t) ∈ U(t, x(t)) and v(t) = h(t, x(t), u(t)) for almost all t ∈ [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 : m → m be continuous. Suppose there exists a constant ρ > 0 such that (6.5) holds for all q ∈ G( ). If the DI: ẋ ∈ 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) ≤ x 0 + t 0 ρ F ( 1 + x(s) ) ds.
By Gronwall's lemma, we obtain

x(t) ≤ ( x 0 + ρ F T ) exp(ρ F T).
( u(t) -u(t) ) T ( G(t, x(t)) + F(u(t)) ) dt ≥ 0,
because the integrand is nonnegative for almost all t. This shows that the pair (x, u) is a weak solution of the DVI.

The linear growth of VI solutions

Although there is an extensive literature on the VI, including the study of boundedness of solutions (see [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]), 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 : m → m is monotone plus on a convex set K ⊆ m if H is monotone on K, i.e., (uu ) T (H(u) -H(u )) ≥ 0, for all u and u ∈ K and the following "plus" property holds:

[ u, u ∈ K, ( u -u ) T ( H(u) -H(u ) ) = 0 ] ⇒ H(u) = H(u ).
It is known that if H is a monotone map on the closed convex set K, then 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 ∈ SOL(K, H), provided that the latter set is nonempty; see [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]. 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 ≡ E T • • E for some matrix E ∈ ×m and a mapping : → that is strongly monotone on EK ⊆ ; that is, there is a constant η > 0 where

( u -u ) T ( (u ) -(u) ) ≥ η u -u 2 , ∀ u , u ∈ 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 ⊆ m be a nonempty closed convex set. Let F : m → m be a continuous map. Suppose that there exists u ref ∈ K satisfying

lim inf u∈K u →∞ ( u -u ref ) T F(u) u 2 > 0. (6.6)
For every vector r ∈ m , SOL(K, r + F) = ∅; moreover, a constant ρ > 0 exists such that (6.5) holds for all r ∈ m . Finally, if F is monotone, then 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 ρ exists. There exist sequences {r k } and {u k } such that, for every k,

( u -u k ) T ( r k + F(u k ) ) ≥ 0, ∀ u ∈ K and u k > k(1 + r k ).
The sequence {u k } clearly is unbounded. Moreover, we have

( u ref -u k ) T ( r k + F(u k ) ) ≥ 0
Dividing by u k 2 and letting k → ∞, 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:

lim inf u∈K u →∞ max 1≤ν≤N F ν (u) T ( u -u ref ) ν 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 ⊆ m be a nonempty closed convex set. Let F ≡ E T • • E, where E ∈ ×m and :

→ is continuous and strongly monotone on EK. If K ∞ ∩ ker E = {0}, then SOL(K, r + F) is a nonempty convex set for all r ∈ m ; moreover, there exists ρ > 0 such that (6.5) holds for all r ∈ m .

Proof By Theorem 2.3.16 in [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF], the nonemptiness of SOL(K, r+F) will follow if we can show that

K ∞ ∩ [-((r + F(K)) * )] = {0}. For this purpose let v be a nonzero element of K ∞ ∩ [-((r + F(K)) * )]. We then have, for any x ∈ K and τ > 0, x + τ v ∈ K, so v T r + ( Ev ) T (Ex + τ Ev) ≤ 0, ∀ τ > 0.
By the strong monotonicity of , a constant η > 0 exists such that

η τ 2 Ev 2 ≤ τ ( Ev ) T ( (Ex + τ Ev) -(Ex) ) ≤ -τ v T r -τ ( Ev ) T (Ex).
Since this holds for all τ > 0, we must have Ev = 0, which means that 0 = v ∈ K ∞ ∩ 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 ∞ ∩ ker E = {0} instead of (6.6). Details are not repeated.

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 [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]. Let D be an m × m matrix. The solution set of the following homogeneous CP is called the VI kernel of the pair (K, D) and is denoted by K(K, D):

K ∞ v ⊥ Dv ∈ ( K ∞ ) * , (6.7)
where (K ∞ ) * is the dual cone of K ∞ . Nonzero elements of 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 K(K, D) = {0}. We denote the interior of the dual cone of the VI kernel K(K, D) by int K(K, D) * . When D is positive semidefinite, the following properties hold: 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.

(a) (D + D T )SOL(K, r, D) is a singleton for all r ∈ m [38, Lemma 2.4.14]; (b) K(K, D) = {v ∈ K ∞ : -D T v ∈ (K ∞ ) * } = K ∞ ∩ (-DK ∞ ) * [38, Lemma 2.5.2]; in particular, K(K, D) is a polyhedral cone if K is polyhedral; (c) if K is a cone (polyhedrality is not needed), then SOL(K,
( D + D T ) ( u -u ) ≤ L D r -r ,
for all u ∈ SOL(K, r, D), u ∈ SOL(K, r , D), and r and r in K(K, D).

Proof We first show that SOL(K, r, D) = ∅ for all r ∈ int K(K, D) * . By degree theory, it suffices to show that for some vector u ref ∈ K, the union t∈(0,1)

SOL(K, tr -(1 -t)u ref , tD + (1 -t)I),
i.e., the set of vectors u ∈ K such that

( u -u ) T [ ( 1 -t ) ( u -u ref ) + t ( r + Du ) ] ≥ 0, ∀ u ∈ K
for some t ∈ (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 } ⊂ (0, 1) and a sequence of vectors {u k } ⊂ K such that, for every k,

( u -u k ) T [ ( 1 -t k ) ( u k -u ref ) + t k ( r + Du k ) ] ≥ 0, ∀ u ∈ K,
and { u k } → ∞. Without loss of generality, we may assume that

lim k→∞ t k = t ∞ and lim k→∞ u k -u ref u k -u ref = v ∞ ,
for some scalar t ∞ ∈ [0, 1] and some nonzero vector v ∞ , which must belong to K ∞ . Letting u = 0, we have, by the positive semidefiniteness of D,

0 ≤ -( u k ) T [ ( 1 -t k ) ( u k -u ref ) + t k ( r + Du k ) ] ≤ -( u k ) T [ ( 1 -t k ) ( u k -u ref ) + t k r ].
Dividing by u ku ref 2 and letting k → ∞, we deduce that t ∞ = 1. Let d be an arbitrary vector in K ∞ . We have

d T [ ( 1 -t k ) ( u k -u ref ) + t k ( r + Du k ) ] ≥ 0.
Dividing by u ku ref and letting k → ∞, we deduce that d T Dv ∞ ≥ 0. Consequently,

Dv ∞ ∈ (K ∞ ) * . Since 0 ≤ ( 1 -t k ) ( u k ) T u ref -t k ( u k ) T ( r + Du k ), dividing by u k -u ref 2 and letting k → ∞, we deduce (v ∞ ) T Dv ∞ ≤ 0. Hence v ∞
is a nonzero solution of the CP (6.7). From the last inequality, we have

0 ≤ ( 1 -t k ) ( u k ) T u ref -t k ( u k ) T r.
Dividing by u k -u ref , letting k → ∞, and using t ∞ = 1, we deduce (v ∞ ) T r ≤ 0. Since v ∞ 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 ∞ > 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 ρ is similar to the previous proof. This completes the proof of statement (a). Assume next that K is polyhedral. Write

K ≡ { u ∈ m : Au ≤ 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 } ⊂ int K(K, D) * and an unbounded sequence {u k } such that u k ∈ 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 λ k such that

0 = r k + Du k + A T λ k , 0 ≤ λ k ⊥ b -Au k ≥ 0,
for every k. Without loss of generality, we may assume that

lim k→∞ u k = ∞, lim k→∞ u k u k = v ∞ ,
for some nonzero vector v ∞ which must belongs to K ∞ = {v ∈ n : Av ≤ 0}.

For each j such that (Av ∞ ) j < 0, we must have (Au k -b) j < 0 for all k sufficiently large, which implies λ k j = 0 for all such k. Let I ≡ {i :

(Av ∞ ) i = 0}. Since 0 = r k u k + D u k u k + i∈I λ k i ( A i• ) T
has a solution λ k i ≡ λ k i / u k ≥ 0 for all k, it follows that for each k sufficiently large, the above equation has a basic feasible solution { λ k i } i∈I with the property that the vectors {(A i• ) T } corresponding to the positive components of λ k i (i ∈ I) are linearly independent. Moreover, we may choose each λ k i such that λ k i > 0 ⇒ (Au kb) i = 0. Hence there exist an infinite subset κ of {1, 2, . . .} and an index set I ∞ ⊆ I such that I ∞ = {i : λ k i > 0} for all k ∈ κ. Since the vectors in {(A i• ) T } i∈I ∞ are linearly independent, it follows that the sequence { λ k I } k∈κ is bounded. Without loss of generality, assume that it converges to a vector λ ∞ I , which must be nonnegative and satisfies

0 = Dv ∞ + i∈I λ ∞ i ( A i• ) T . Furthermore, if λ ∞ i > 0, then λ k i > 0, hence (Au k -b) i = 0,
for all k ∈ κ sufficiently large. Therefore, letting λ ∞ i ≡ 0 ≡ λ k i for all i ∈ I and for all k ∈ κ, we have 0 = Dv ∞ +A T λ and ( λ ∞ ) T (Au k -b) = 0 for all k ∈ κ sufficiently large. Since 0 ∈ K, it follows that b ≥ 0. Hence ( λ ∞ ) T Au k ≥ 0 for all k ∈ κ sufficiently large, which implies

(v ∞ ) T Dv ∞ = 0; thus v ∞ ∈ K(K, D). Since D is positive semidefinite, it follows that (D+D T )v ∞ = 0; Hence, A T λ ∞ = -Dv ∞ = D T v ∞ , which yields (Du k ) T v ∞ ≥ 0 for all k ∈ κ sufficiently large. Consequently, 0 = ( v ∞ ) T r k u k + D u k u k + A T λ k = ( v ∞ ) T r k u k + D u k u k ≥ ( v ∞ ) T r k u k .
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, The LCP kernel of D is {0} × n + ; the dual of this kernel is therefore n × n + . Clearly, the vector r belongs to the interior of this dual for any x ∈ n . Therefore, part (b) of Proposition 6.4 applies. The reader can easily verify that

sup{ u : u ∈ SOL(r, D) } ≤ 2 ( 1 + r ), ∀ 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 [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF][START_REF] Luo | A decomposition property for a class of square matrices[END_REF]; i.e., D is positive semidefinite and

[ u T Du = 0 ⇒ Du = 0 ].
This is equivalent to saying the linear map u → Du is monotone-plus. The following known result [START_REF] Luo | A decomposition property for a class of square matrices[END_REF] shows among other things that for a linear map, the strongly monotone (linear) composite property is equivalent to the monotoneplus property. If D is a psd-plus matrix, then Du is a constant for all u ∈ SOL(K, r, D) [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]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 ⊆ m be a polyhedron containing either the origin or no lines. Let D ∈ m×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 ∈ m . Moreover, there is a constant L V > 0, such that Du 1 -Du 2 ≤ L V r 1r 2 , (6.9)

for all r i in m and u i ∈ 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 ∞ ∩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 • •E in Proposition 6.3.
Yet, for the latter map, while Eu is a constant for u ∈ SOL(K, r + F), the (singlevalued) map r → 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 [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]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 + Φ is not necessarily monotone. Proposition 6.5 Let K ⊆ m be a polyhedron containing either the origin or no lines. Let D ∈ m×m be a psd-plus matrix such that (K, D) is an R 0 pair. Let ρ > 0 be such that (6.8) holds. Let Φ : m → m be a continuous function such that Φ(u) ≤ L Φ u , ∀ u ∈ K, (6.10)

for some constant L Φ ∈ (0, 1/ρ). For every vector q ∈ m , SOL(K, q + F) = ∅, where F ≡ D + Φ; moreover,

sup{ u : u ∈ SOL(K, q + F) } ≤ ρ ( 1 + q ) 1 -ρ L Φ . ( 6

.11)

Assume in addition

Φ(u) -Φ(u ) ≤ L Φ Du -Du , ∀ u, u ∈ K, (6.12)
for some constant L Φ ∈ (0, 1/L V ), where L V > 0 is such that (6.9) holds for all r i in m and u i ∈ SOL(K, r i , D), i = 1, 2. It holds that

Du 1 -Du 2 ≤ L V q 1 -q 2 1 -L V L Φ , (6.13)
for all q i ∈ m and u i ∈ SOL(K, q i + F), i = 1, 2. Finally, for every q ∈ m ,

SOL(K, q + F) = D -1 (v(q)) ∩ argmin{ u T w(q) : u ∈ K }, (6.14)
where v(q) ≡ D ũ and w(q) ≡ q+F( ũ) for any ũ ∈ SOL(K, q+F), and D -1 (v(q)) is the inverse image of v(q) under the ( possibly singular) matrix D. Consequently, SOL(K, q + F) is a convex set.

Proof Choose an arbitrary scalar R > ρ ( 1 + q ) 1 -ρ L Φ such that the intersection, denoted U ≡ K ∩ RB, of K with the Euclidean ball RB in m centered at the origin and with radius R is nonempty. We define a set-valued mapping from U into subsets of U as follows. For u ∈ U, let (u) be the solution set of the VI (K, r(u), D), where r(u) ≡ q + Φ(u). Clearly (u) is a nonempty closed convex subset of K. We need to show that v ≤ R for all v ∈ (u). Indeed, we have

v ≤ ρ ( 1 + r(u) ) ≤ ρ ( 1 + q + L Φ R ) < R.
Since Φ is continuous, it is easy to see that is a closed (set-valued) map. Consequently, by Kakutani's fixed-point theorem, there exists u ∈ (u). This vector u clearly is an element of SOL(K, q + F). Moreover, the bound (6.11) clearly follows from the above bound on v .

Notice that if u ∈ SOL(K, q + F), then u ∈ SOL(K, q + Φ(u), D). Consequently, if u i ∈ SOL(K, q i + F) for i = 1, 2, then we have

Du 1 -Du 2 ≤ L V q 1 + Φ(u 1 ) -q 2 -Φ(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 ∈ SOL(K, q + F); hence, by (6.12), so is Φ(u), thus so is F(u). Consequently, the inclusion

SOL(K, q + F) ⊆ D -1 (v(q)) ∩ argmin{ u T w(q) : u ∈ 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 ũ; thus Φ(u ) = Φ( ũ). Hence w(q) = q + F( ũ) = q + F(u ). Consequently, for all u ∈ K, 0 ≤ ( u -u ) T w(q) = ( u -u ) T ( q + F(u ) ),
establishing the equality (6.14). This representation of SOL(K, q + F) clearly shows that this solution set is convex.

Although the VI (K, q + F) in Proposition 6.5 has a convex solution set, the map F ≡ D+Φ 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 × 2 matrix of all ones. Let Φ(u) ≡ Eu, where

E ≡ ε ε 0 0 .
For any ε > 0 sufficiently small, the matrix

D + E = 1 + ε 1 + ε 1 1
is column sufficient [30, Definition 3.5.1] but not positive semidefinite. Yet with such an ε, all other conditions in Proposition 6.5 are valid for the pair (D, Φ). 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( ) ⊆ int K(K, D) * ; while the latter inclusion is trivially valid when (K, D) is an R 0 pair, condition (c) allows K to be non-polyhedral. where D is a psd-plus matrix such that (K, D) is an R 0 pair, and where Φ : m → 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, F(x, t) is nonempty and convex, and by Proposition 6.4(a), F satisfies the linear growth property (6.5). (d) By Proposition 6.4, F(x, t) is nonempty and convex, and by Proposition 6.4(b), F satisfies the linear growth property (6.5). (e) By Proposition 6.5, and equation (6.14), F(x, t) is non-empty and convex; by (6.11), F satisfies the linear growth property (6.5).

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:

ẋ = p + Ax + Bu, C u ⊥ q + Cx + Du ∈ C * , (6.15)
where C is a polyhedral cone, D is positive semidefinite, and q + C n ⊆ int K(C, D) * . This corollary is related to similar results for the LCS (2.8) with p = 0 and q = 0 established in [START_REF] Çamlibel | Complementarity methods in the analysis of piecewise linear dynamical systems[END_REF][START_REF] Heemels | Linear complementarity systems: a study in hybrid dynamics[END_REF], where it is further assumed that the tuple (A, B, C, D) is "minimal" and "passive". (See also [START_REF] Goeleven | Stability and instability matrices for linear evolution variational inequalities[END_REF] which contains further results on passivity and the LCS.) There is one important difference, however.

Our result shows that if q + C n ⊆ int K(C, D) * , then (6.15) has a solution for all initial conditions x 0 ∈ n ; in contrast, a result in the cited references shows that, in the case where p = 0 and q = 0 and C = 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.

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, •) is Lipschitz with Lipschitz constant L f for all t, f (•, x) is measurable for all x ∈ n , and

f (t, 0) ≤ ρ f ,0 for all t ∈ [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 t h,i+1 t h,i f (τ , z h,i ) dτ . 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 ≡ x 0 given, we compute {x h,1 , x h,2 , . . . , x h,N h +1 } ⊂ n and {u h,1 , u h,2 , . . . , u h,N h +1 } ⊂ K (7.1) by the recursion: for i = 0, 1, . . . , N h , where N h ≡ T/h -1,

x h,i+1 = x h,i + h[f (t h,i+1 , θ x h,i + (1 -θ)x h,i+1 ) + B(t h,i , x h,i )u h,i+1 ], u h,i+1 ∈ SOL(K, G(t h,i+1 , x h,i+1 ) + F). (7.2)
Despite the similarity with some schemes developed for mechanical systems subject to unilateral constraints and Coulomb friction (see e.g. [START_REF] Anitescu | Time-stepping for three-dimensional rigid-body dynamics[END_REF][START_REF] Jean | The nonsmooth contact dynamics approach[END_REF][START_REF] Pang | Convergence of time-stepping methods for initial and boundary value frictional compliant contact problems[END_REF][START_REF] Song | Semi-implicit time-stepping models for frictional compliant contact problems[END_REF][START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of Painlevé's problem[END_REF][START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction[END_REF]), 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 , θ x h,i + (1 -θ)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 ) - m j=1 v h,i+1 j ∇g i (x h,i+1 ), 0 ≤ v h,i+1 ⊥ -g(x h,i+1 ) ≥ 0, (7.3)
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) ≡ xx 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 ) ≤ 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 [START_REF] Dirkse | The PATH solver: A non-monotone stabilization scheme for mixed complementarity problems[END_REF][START_REF] Ferris | Interfaces to PATH 3.0: Design, implementation and usage[END_REF] 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 ∈ (0, h 0 ], (x ref , u) ∈ n+m , θ ∈ [0, 1], and t and t ref in [0, T], there exists a unique vector x(u) satisfying

x(u) -x ref = h [ f (t, θx ref + (1 -θ)x(u)) + B(t ref , x ref )u ].
Moreover, for all u, u in m ,

x(u) -x(u ) ≤ h σ B 1 -h ( 1 -θ ) L f u -u , x(u) -x ref ≤ h ρ f ( 1 + x ref ) + σ B u 1 -h ( 1 -θ ) ρ f . (7.4)
Proof It suffices to choose h 0 to satisfy

0 < h 0 < min 1 ( 1 -θ ) L f , 1 ( 1 -θ ) ρ f .
(The right-hand side is taken to be 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 x h (•) be the continuous piecewise linear interpolant of the family {x h,i }; i.e., for i = 0, 1, . . . , N h ,

∞ if θ = 1.) Under
x h (t) ≡ x h,i + t -t h,i h ( x h,i+1 -x h,i ), ∀ t ∈ [ t h,i , t h,i+1 ];
let u h (•) be the (possibly discontinuous) piecewise constant interpolant of the family {u h,i }, i.e., u h (t) ≡ u h,i+1 for t ∈ (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 ≡ • E, where E ∈ ×m and : → m . We let w h (•) 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 SOL(K, G(0, x 0 ) + F), provided that the latter VI is solvable. We have, for i = 0, 1, . . . , N h ,

w h (t) ≡ Eu h,i + t -t h,i h E( u h,i+1 -u h,i ), ∀ t ∈ [ t h,i , t h,i+1 ].
Consequently, for every h > 0 sufficiently small, the functions x h : [0, T] → n , u h : [0, T] → m , and w h : [0, T] → 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 ⊆ 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 h such that for all h ∈ (0, h] and all i = 0, 1, . . . , N h , x h,i+1 ≤ c 0,x + c 1,x x 0 and u h,i+1 ≤ c 0,u + c 1,u x 0 . (7.5)

There is a sequence {h ν } ↓ 0 such that the following two limits exist: x h ν → x uniformly on [0, T] and u h ν → u weakly in L 2 (0, T). Furthermore, under either one of the following two conditions:

(a) F(u) ≡ (Eu)
, where E ∈ ×m and : → 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 ≤ h c 2,u , (7.6)
(b) F(u) ≡ Du for some positive semidefinite matrix D, all such limits ( x, u) are weak solutions of the initial-value DVI (6.2).

Proof Throughout the proof below, h ∈ (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 ≤ h ρ f ( 1 + x h,i ) + σ B u h,i+1 1 -h ( 1 -θ ) ρ 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, . . . , N h , x h,i+1x h,i ≤ L x 0 h. This implies that the piecewise interpolants 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 { 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

x h L ∞ ≡ sup t∈[0,T]
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 ν } ↓ 0 such that { x h ν } converges in the supremum (i.e., L ∞ ) norm to a Lipschitz function 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 { u h } are uniformly bounded in the L ∞ norm on (0, T); that is, for some constant ψ u > 0, u h L ∞ ≤ ψ u for all h > 0 sufficiently small. By Alaoglu's theorem [63, pp. 71-72], and by working with an appropriate subsequence of {h ν } if necessary, we may assume without loss of generality that the sequence { u h ν } has a weak* limit 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 ϕ ∈ L 2 (0, T),

lim ν→∞ T 0 ϕ(t) T u h ν (t) dt = T 0 ϕ(t) T u(t) dt.
Mazur's theorem [63, p. 88] implies that there is a sequence of convex combinations of { u h ν } that converges strongly in L 2 (0, T). Since u h (t) ∈ K for all t, and K is a closed convex set, any convex combination of { u h (t)} will also be in K. Thus, 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 < ∞) have subsequences that are pointwise convergent almost everywhere [63, p. 210], it follows that u(t) ∈ K for almost all t.

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

[0, T] → K, T 0 ( u(t) -u(t) ) T [ G(t, x(t)) + (E u(t)) ] dt ≥ 0, (7.7)
(ii) for any 0

≤ s ≤ t ≤ T, x(t) -x(s) = t s [ f (τ , x(τ )) + B(τ , x(τ )) u(τ ) ] dτ , (7.8)
and (iii) the initial condition x(0) = x 0 . The last condition (iii) holds since x h (0) = x 0 for all h > 0 sufficiently small and x h ν → x uniformly as ν → ∞. To show (i), note that the sequence {G(t,

x h ν ) + (E u h ν )} converges to G(t, x) + (E u) uniformly. This implies that lim ν→∞ T 0 ( u(t) -u h ν (t) ) T G(t, x h ν (t)) + (E u h ν (t)) dt = T 0 ( u -u ) T [ G(t, x) + (E u) ] dt.
On the other hand,

T 0 ( u -u h ) T G(t, x h ) + (E u h ) dt = N h i=0 t h,i+1 t h,i ( u(t) -u h,i+1 ) T G(t, x h (t)) + (Eu h,i+1 ) dt ≥ h N h i=0 ⎛ ⎜ ⎝ 1 h t h,i+1 t h,i u(t) dt -u h,i+1 ⎞ ⎟ ⎠ T G(t h,i+1 , x h,i+1 ) + (Eu h,i+1 ) -h 2 ( N h + 1 ) u L ∞ + u h L ∞ L G L x ≥ -h T ( u L ∞ + u h L ∞ ) L G L x ,
where the last inequality holds because 1 h t i+1 t i u(t)dt belongs to K by convexity of K and also because u h,i+1 ∈ SOL(K, G(t h,i+1 , x h,i+1 ) + F). Hence,

T 0 ( u -u ) T [ G(t, x) + (E u) ] dt ≥ 0
for any continuous u : [0, T] → K. To complete the convergence proof under (a), it remains to show (7.8). Similar to the above derivation, we have

x h,i+1 -x h,i = h [ f (t h,i+1 , θ x h,i + (1 -θ)x h,i+1 ) + B(t h,i , x h,i )u h,i+1 ] = t h,i+1 t h,i [ f (τ , x h (τ )) + B(τ , x h (τ )) u h,i+1 ] dτ + h 2 ζ i , where ζ i ≤ (L f + L f L x + σ B ψ u ). It follows that for any 0 ≤ s ≤ t ≤ T, x h (t) -x h (s) = t s [ f (τ , x h (τ )) + B(τ , x h (τ )) u h (τ ) ] dτ + O(h).
By continuity of f and uniform convergence of x h ν → x, we deduce

lim ν→∞ t s f (τ , x h ν (τ )) dτ = t s f (τ , x(τ )) dτ . Furthermore, t s [ B(τ , x h (τ )) u h (τ ) -B(τ , x(τ )) u(τ ) ] dτ ≤ t s B(τ , x h (τ )) -B(τ , x(τ )) u h (τ ) dτ + t s [B(τ , x(τ )) u h (τ ) -B(τ , x(τ )) u(τ )]dτ .
Restricted to the subsequence {h ν }, the first summand on the right-hand side converges to zero because { x h ν } → x uniformly and B is continuous; the second summand converges to zero because { u h ν } converges weakly in L 2 (0, T) to u and B(τ , x(τ )) is continuous in τ . Consequently, (7.8) holds.

Next assume condition (b). It suffices to show that for all continuous u :

[0, T] → K, T 0 ( u(t) -u(t) ) T [ G(t, x(t)) + D u(t) ] dt ≥ 0,
In turn, it suffices to show that

lim sup ν→∞ T 0 ( u(t) -u h ν (t) ) T G(t, x h ν (t)) + D u h ν (t) dt ≤ T 0 ( u(t) -u(t) ) T [ G(t, x(t)) + D u(t)
] dt because the left-hand limit is nonnegative as shown above. Since

lim ν→∞ T 0 u(t) T G(t, x h ν (t)) + D u h ν (t) dt = T 0 u(t) T [G(t, x(t)) + D u(t)] dt and lim ν→∞ T 0 ( u h ν (t) ) T G(t, x h ν (t)) dt = T 0 ( u(t)) T G(t, x(t)) dt, it remains to show lim inf ν→∞ T 0 ( u h ν (t) ) T D u h ν (t) dt ≥ T 0 ( u(t) ) T D u(t) dt. ( 7.9)
We can write

( u h ν (t) ) T D u h ν (t) = ( u(t) ) T D u(t) + ( u h ν -u(t) ) T ( D + D T ) u(t) +( u h ν (t) -u(t) ) T D( u h ν -u ) ≥ ( u(t) ) T D u(t) + ( u h ν (t) -u(t) ) T ( D + D T ) u(t).
Since u h ν u weakly in L 2 (0, T) and (D + D T ) u ∈ L 2 (0, T), it follows that

lim ν→∞ T 0 ( u h ν (t) -u(t) ) T ( D + D T ) 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 , ρ u , and ψ x such that for all h ∈ (0, h 1 ) and all nonnegative integers i with

(i + 1)h ≤ T, u h,i+1 ≤ ρ u ( 1 + 2 x h,i ), (7.10) x h,i+1 -x h,i ≤ h ψ x ( 1 + x h,i ). (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 ∈ (0, h 1 ] and all i = 0, 1, . . .

, N h ≡ T/h -1.
Proof For all i = 0, 1, . . . , N h , x h,i+1 ≤ (1

+ h ψ x ) x h,i + h ψ x , yielding x h,i+1 ≤ ( 1 + h ψ x ) i+1 x 0 + h ψ x i j=0 ( 1 + h ψ x ) j ≤ ( 1 + h ψ x ) i+1 x 0 + h ψ x ( 1 + h ψ x ) i+1 -1 h ψ x ≤ e h (i+1) ψ x x 0 + e h (i+1) ψ x -1 ≤ e T ψ x x 0 + e T ψ x -1,
where the last inequality holds because (i+1)h ≤ (N h +1)h ≤ T. Consequently, with c 1,x ≡ e Tψ x ≡ 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).

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 ρ > 0, the solutions to the VI (K, q + F) satisfies the linear growth property (6.5) for all q ∈ G( ). There exist positive scalars c 0,x , c 1,x , c 0,u , and c 1,u and h 1 such that (7.5) holds for all h ∈ (0, h 1 ] and all i = 0, 1, . . . , N h . Proof Throughout the proof below, the scalar h > 0 is taken to be sufficiently small. Letting

ρ x ≡ ρ f + σ B 1 -h 0 ( 1 -θ ) ρ f , (7.12)
it follows from (7.4) that

x h,i+1 -x h,i ≤ h ρ x ( 1 + x h,i + u h,i+1
). (

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

u h,i+1 ≤ ρ ( 1 + G(t h,i+1 , x h,i+1 ) ) ≤ ρ [ 1 + ρ G ( 1 + x h,i+1 ) ) ] ≤ ρ [ 1 + ρ G { 1 + x h,i + h ρ x ( 1 + x h,i + u h,i+1 ) } ].

Consequently, letting

ρ

IVP 1 u ≡ 2 ρ ( 1 + ρ G + ρ G ρ x )
we obtain, for all h ∈ (0, 1) satisfying h < 1/(2ρρ G ρ x ),

u h,i+1 ≤ ρ IVP 1 u ( 1 + 2 x h,i ).
Therefore, with ψ

IVP 1 x ≡ ρ x (1 + 2ρ IVP 1 u
), we deduce

x h,i+1 -x h,i ≤ h ψ IVP 1 x ( 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 [START_REF] Eilenberg | Fixed point theorems for multi-valued transformations[END_REF] and the other due to Górniewicz for admissible maps [START_REF] Górniewicz | Homological methods in fixed-point theory of multi-valued maps. Dissertationes Mathematicae[END_REF][START_REF] Górniewicz | Topological Fixed Point Theory of Multivalued Mappings[END_REF]; the recent paper [START_REF] Saveliev | Fixed points and selections of set-valued maps on spaces with convexity[END_REF] 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 [START_REF] Górniewicz | Topological Fixed Point Theory of Multivalued Mappings[END_REF].

An acyclic set is a topological space X where the rational homology groups of X are isormorphic 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 → Y is a Vietoris map if it is continuous and for each y ∈ 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 → Y with compact values which can be represented as F = r • p -1 where r : Z → Y is a continuous single-valued map and p : Z → 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 → X on a compact AR X has a fixed point: x ∈ F(x) for some x ∈ X. Theorem 7.3 (Górniewicz) Every admissible multifunction F : X → X on a compact AR X has a fixed point: x ∈ F(x) for some x ∈ X.

Returning to the IVP iteration (7.2), we note that the constant ψ 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

ψ IVP 1 x ≡ ρ f + σ B 1 -h 0 ( 1 -θ ) ρ f 1 + 4 ρ ( 1 + ρ G + ρ G ρ x ) ,
and recall that (i) h 0 is the constant asserted by Lemma 7.1, (ii) ρ f and ρ G are the respective linear growth constants of f and G implied by assumption (A); cf. (6.3); (iii) σ B is the bound of B(t, x) on in assumption (B), (iv) ρ is the linear growth constant of the solutions to the VIs (K, q + F), (v) ρ x is the constant defined in (7.12), and (vi) θ ∈ [0, 1] is the parameter used in the iterative scheme (7.2). Most importantly, ψ

IVP 1 x
is independent of the step size h and any other properties of the functions (f , B, G, F). The role of ψ IVP 1 x 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, . . . , N h . Proposition 7.1 Let K ⊆ 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 ρ > 0, SOL(K, q + F) is a nonempty convex set satisfying the linear growth property (6.5) for all q ∈ G( ). A scalar h R > 0 exists such that for all h ∈ (0, h R ], θ ∈ [0, 1], and x 0 ∈ n , a pair (x h,i+1 , u h,i+1 ) exists satisfying (7.2) for every i = 0, 1, . . . , N h .

Proof Associated with the constant ψ

IVP 1 x
, we define, for any scalar h > 0 sufficiently small, the scalars ρ 1 , ρ 2 , . . ., and ρ N h +1 recursively by

ρ i+1 ≡ ( 1 + h ψ IVP 1 x ) ρ i + h ψ IVP 1 x , i = 0, 1, . . . , N h ,
where ρ 0 is arbitrary. By the proof of Lemma 7.3, we can show that

ρ i ≤ e T ψ IVP 1 x ρ 0 + e T ψ IVP 1 x -1, ∀ i = 0, 1, . . . , N h + 1.
Let α denote the quantity on the right-hand side, which depends additionally on ρ 0 , but is also independent of h. Let h R > 0 satisfy h R < min(h 0 , h 1 ) and also

h R ρ f ( 1 + α ) + σ B ρ ρ G ( 1 + 2 α ) 1 -h R ( 1 -θ ) ρ f < α,
where h 0 and h 1 are as described in Lemmas 7.1 and 7.3, respectively. Fix an arbitrary h ∈ (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 ≤ ρ i+1 . Let B α denote the Euclidean ball in n with center at the origin and radius 2α. For any x ∈ B α , let S(t, x) denote the nonempty convex set SOL(K, G(t, x) + F). By the linear growth assumption, we have

sup{ u : u ∈ S(t, x) } ≤ ρ ρ G ( 1 + 2 α ).
Define the map S i from B α into subsets of B α as follows: for x ∈ B,

S i (x) ≡ ( I -h f (t h,i+1 , θ x h,i + (1 -θ)•) ) -1 ( x h,i + h B(t h,i , x h,i )S(t h,i+1 , x) ).
Since the map (Ihf (t h,i+1 , θ x i + (1 -θ)•)) -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 B. Let x be an arbitrary element in S i (x) and let u ∈ S(t h,i+1 , x) be such that

x = x h,i + h [ f (t h,i+1 , θ x i + (1 -θ) x) + B(t h,i , x h,i )u ] From (7.4), we have x -x h,i ≤ h ρ f ( 1 + x h,i ) + σ B u 1 -h ( 1 -θ ) ρ f ,
which yields, by induction hypothesis,

x ≤ ρ i + h ρ f ( 1 + ρ i ) + σ B ρ ρ G ( 1 + 2 α ) 1 -h ( 1 -θ ) ρ f < 2 α, because x h,i ≤ ρ i ≤ α.
Since the solution map S(t h,i+1 , x) is upper semicontinuous by Lemma 6.2, it follows that S i : B α → B α 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 ≤ ρ i+1 . But this follows from the proof of Lemma 7.3 and the definition of ρ i+1 .

We now have all the preparatory results to establish the promised corollary of Theorem 7. 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).

Like part (d) of Theorem 6.1, part (c) of Theorem 7.4 is applicable to the 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 n ⊆ int K(C, D) * .

Dividing by τ 2 and letting τ ↓ 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 ref , x ref )u) -G(t, r + B(t ref , x ref )u ) ) = 1 0 (u -u ) T J x G(t, r + B(t ref , x ref )(u + τ (u -u )))B(t ref , x ref )(u -u )dτ ≥ η 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 , θ x h,i + (1 -θ)x h,i+1 ) + B(t h,i , x h,i )u h,i+1 ], u h,i+1 ∈ SOL(K, G h,i+1 + F), ( 8.2)
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 ) + ∂G(t h,i , x h,i ) ∂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) ∈ , 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 ∅ = K ⊆ m be closed and convex. Let ( f , B, G, F) satisfy conditions (A)-(D). If 0 < h < min η G ( 1 -θ ) L f ( L G σ B + η G ) , 1 ( 1 -θ ) ρ f , (8.3)
then there exists a pair (x h , u h ) satisfying

x h = x ref + h f (t, θ x ref + (1 -θ)x h ) + B(t ref , x ref )u h , u h ∈ SOL(K, G(t, x h ) + F). (8.4)
Moreover, the pair (x h , u h ) is unique if F is additionally monotone on K.

the vector on the left-hand side is clearly equal to Zµ 1 + Wλ 2 . Thus the right-hand side of (8.6) follows. Conversely, suppose the latter implication holds. Let u i ∈ K for i = 1, 2 and write

µ i λ i ≡ P T u i = Z T u i W T u i .
Thus, ZZ T u 1 + WW T u 2 = Zµ 1 + Wλ 2 ; since (8.6) holds, the vector in the right-hand side of the last expression belongs to K.

Since is strongly monotone on E m and W has full column rank, it follows that the map ϒ ≡ W T • • W is a strongly monotone self-map from m-k onto itself. Thus, a constant η ϒ > 0 exists such that

( ϒ(λ) -ϒ(λ ) ) T ( λ -λ ) ≥ η ϒ λ -λ 2 , ∀ λ, λ ∈ m-k . ( 8.7)
(When E = 0, we let η ϒ 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+Φ h +F), where Φ h is intended to be a mapping of the kind G • x h , with x h (u) being the unique vector x satisfying the equation

x = x ref + h [ f (t, x ref + (1 -θ)(x -x ref )) + B(t ref , x ref )u ];
cf. the proof of Proposition 8.2. With the nonsingular transformation of variables: u ≡ Pv = Zµ + Wλ, where v ≡ (µ, λ) ∈ m the VI (K, q + Φ h + F) is equivalent to (since P T = P -1 ): (µ, λ) ∈ P T K and for all (µ , λ ) ∈ P T K,

µ -µ λ -λ T Z T q + Z T Φ h (Zµ + Wλ) W T q + ϒ(λ) + W T Φ h (Zµ + Wλ) ≥ 0.
Notice that Eu = Wλ. Based on the above preparations, we state and prove the following result.

Lemma 8.2 Let F satisfy condition (C ).

Let Φ h : m → m be strongly monotone on m with constant h η Φ > 0 and Lipschitz continuous with constant h L Φ > 0. Provided that h > 0 is sufficiently small satisfying

0 < h < 1 2 η ϒ L Φ W 2 1 + L Φ η Φ , ( 8.8)
the map q → ESOL(K, q + Φ h + F) is single-valued and Lipschitz continuous with a constant independent of h.

Proof For i = 1, 2, let u i ∈ SOL(K, q i + Φ h + F) and write P T u i ≡ (µ i , λ i ). By Lemma 8.1, (µ 1 , λ 2 ) and (µ 2 , λ 1 ) both belong to P T K. Hence

( µ 2 -µ 1 ) T [ Z T q 1 + Z T Φ h (Zµ 1 + Wλ 1 ) ] ≥ 0, ( λ 2 -λ 1 ) T [ W T q 1 + ϒ(λ 1 ) + W T Φ h (Zµ 1 + Wλ 1 ) ] ≥ 0, ( µ 1 -µ 2 ) T [ Z T q 2 + Z T Φ h (Zµ 2 + Wλ 2 ) ] ≥ 0, ( λ 1 -λ 2 ) T [ W T q 2 + ϒ(λ 2 ) + W T Φ h (Zµ 2 + Wλ 2 ) ] ≥ 0.
Adding the first and third inequality and rearranging terms, we obtain

( q 2 -q 1 ) T Z( µ 1 -µ 2 ) ≥ ( Z( µ 2 -µ 1 ) ) T [ Φ h (Zµ 2 + Wλ 2 ) -Φ h (Zµ 1 + Wλ 1 ) ] = ( Z( µ 2 -µ 1 ) ) T [ Φ h (Zµ 2 + Wλ 2 ) -Φ h (Zµ 2 + Wλ 1 ) + Φ h (Zµ 2 + Wλ 1 ) -Φ h (Zµ 1 + Wλ 1 ) ] ≥ η Φ h Z( µ 2 -µ 1 ) 2 -L h Z( µ 2 -µ 1 ) W( λ 1 -λ 2 ) .
Consequently,

Z( µ 2 -µ 1 ) ≤ 1 η Φ h q 1 -q 2 + L Φ η Φ W( λ 1 -λ 2 ) .
Similarly, we can deduce

( q 2 -q 1 ) T W( λ 1 -λ 2 ) ≥ η ϒ -h W 2 L 2 Φ η Φ + L Φ λ 1 -λ 2 2 -W( λ 2 -λ 1 ) L Φ η Φ q 1 -q 2 .
Consequently, if (8.8) holds, then

λ 2 -λ 1 ≤ 1 + L Φ η Φ W η ϒ -h L Φ W 2 1 + L Φ η Φ q 1 -q 2 ,
which yields

Eu 1 -Eu 2 ≤ W 1 + L Φ η Φ W η ϒ -h L Φ W 2 1 + L Φ η Φ q 1 -q 2 ,
where the multiplicative constant on the right-hand side is less than a constant that is independent of h.

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 m and θ ∈ [0, 1] be a given scalar. Let (f , B, G) satisfy the conditions (A), (B) and (D). Let F satisfy (C ) and let η ϒ satisfy (8.7). Assume (0) = 0 and condition (E). Let x ref ∈ n be an arbitrary vector in n . If h satisfies (8.3), then for all t and t ref in [0, T], a unique pair (x h , u h ) exists satisfying (8.4). If in addition there exists

u ref ∈ K such that G(t ref , x ref ) + F(u ref ) is in K * , and h > 0 is sufficiently small satisfying in particular 0 < h < η ϒ W 2 L G σ B 1 + L G σ B η G , ( 8.9)
then for some constant ξ > 0 that is independent of h and x ref , u h is bounded above by

ξ 1 + 1 h |t -t | + G(t, x ref ) + h{ f (t, x ref ) + x h -x ref } . ( 8
(µ -µ h ) T Z T G(t, x h ) ≥ 0, (λ -λ h ) T [ ϒ(λ h ) + W T G(t, x h ) ] ≥ 0, (8.12)
for all (µ, λ) ∈ P T K. Since K contains the origin, we have, by substituting µ = 0 in the first of the two inequalities in (8.12) and by (8.11),

0 ≥ (µ h ) T Z T G(t, x h ) ≥ (Zµ h ) T [G(t, x ref ) -G(t ref , x ref ) + G(t, x h ) -G(t, x ref )]. Since x h = x ref + h[f (t, x ref + (1 -θ)(x h -x ref )) + B(t ref , x ref )(Zµ h + Wλ h )], we deduce, 0 ≥ ( Zµ h ) ] T [ G(t, x ref ) -G(t ref , x ref ) + G(t, x ref + h(f (t, x ref + (1 -θ)(x h -x ref )) + B(t ref , x ref )(Zµ h + Wλ h ))) -G(t, x ref + h(f (t, x ref + (1 -θ)(x h -x ref )) + B(t ref , x ref )Wλ h )) + G(t, x ref + h(f (t, x ref + (1 -θ)(x h -x ref ))+B(t ref , x ref )Wλ h ))-G(t, x ref )],
which implies

η G h Zµ h ≤ L G | t -t | + hL G [ f (t, x ref ) + L f ( 1 -θ) x h -x ref + B(t ref , x ref ) Wλ h ].
Hence Zµ h is bounded above by

L G |t -t ref | η G h + L G [ f (t, x ref ) + L f (1 -θ) x h -x ref + σ B Wλ h ] η G . (8.13)
Substituting λ = 0 into the second inequality in (8.12) and using the assumption ϒ(0) = 0 and the strong monotonicity of ϒ, we deduce

η ϒ λ h ≤ W G(t, x h ) ≤ W [ G(t, x ref ) + h L G f (t, x ref + (1 -θ)(x h -x ref ) ] + W h L G B(t ref , x ref ) Zµ h + Wλ h , which implies ( η ϒ -h W 2 L G σ B ) λ h ≤ W h L G σ B Zµ h + W [ G(t, x ref ) + h L G { f (t, x ref ) + L f ( 1 -θ ) x h -x ref } ].
Consequently, from (8.13) and after rearranging terms, we obtain

η ϒ -h W 2 L G σ B 1 + L G σ B η G λ h ≤ W [ G(t, x ref ) ] + W L 2 G σ B η G | t -t ref | + 1 + L G σ B η G W h L G [ f (t, x ref ) + L f ( 1 -θ ) x h -x ref ].
Therefore, if h satisfies (8.9), the multiplicative factor of λ h in the left-hand side is positive. Dividing by this constant factor, we obtain

λ h ≤ c [ G(t, x ref ) + | t -t | + h [ f (t, x ref ) + x h -x 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µ h . Since u h ≤ Zµ h + W λ h , the desired upper bound (8.10) follows readily with an appropriate constant ξ > 0.

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 m and θ ∈ [0, 1] be a given scalar. Let (f , B, G) satisfy the conditions (A), (B) and (D). Let F satisfy (C ) and let η ϒ satisfy (8.7). Assume (0) = 0 and condition (E). A positive scalar h > 0 exists such that for all x 0 ∈ n for which SOL(K, G(t 0 , x 0 ) + F) = ∅ and for all h ∈ (0, h], a unique pair (x h,i+1 , u h,i+1 ) exists satisfying (7.2) for all i = 0, 1, 2, . . . , 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, . . . , 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 ∈ 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 ) ∈ K * . Thus, with t ≡ t h,i+1 , t ref ≡ t h,i , and u ref ≡ u h,i , (8.10) gives

u h,i+1 ≤ ξ [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 }] ≤ ξ [ h + 1 + ( ρ G + h ρ f ) ( 1 + x h,i ) + h x h,i+1 -x h,i } ],
which yields, for some constant ω > 0 that is independent of h and x 0 ,

u h,i+1 ≤ ω ( 1 + x h,i + h x h,i+1 -x h,i ). ( 8.14)
From (7.13), we have x h,i+1x h,i ≤ hρ x (1 + x h,i + u h,i+1 ), which yields

x h,i+1 -x h,i ≤ h ρ x 1 + ω 1 -ω ρ x h 2 ( 1 + x h,i ).
Therefore, for some positive constant ψ IVP 2 x 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 ρ IVP 2 u . 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 h such that (7.5) holds for all h ∈ (0, h] and all i = 0, 1, . . . , 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 ∈ m , let xi+1 (u) denote the unique vector x satisfying

x -x h,i = h[f (t h,i+1 , θ x h,i + (1 -θ)x) + B(t h,i , x h,i )u]. By the proof of Proposition 8.2, the map Φ h : u → G(t h,i+1 , xi+1 (u)) is strongly monotone with modulus h η G - h(1 -θ)L f L G σ B 1 -h(1 -θ)L f , which is greater than h η G /2, provi- ded that h > 0 is sufficiently small. Moreover, Φ h is Lipschitz continuous with modulus hL G σ B 1 -h(1 -θ)L f
, which is less than 2hL G σ B , provided that h > 0 is sufficiently small. Applying Lemma 8.2 to the pair (K, Φ 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 ) ≤ c u q 1q 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 + Φ 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 , xi+1 (u h,i )) + Φ 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 , xi+1 (u h,i )), we deduce Eu h,i+1 -Eu h,i ≤ c u G(t h,i , x h,i ) -G(t h,i+1 , xi+1 (u h,i )) ≤ c u L G [ h + xi+1 (u h,i ) -x h,i ] ≤ h c u L G 1 + ρ f ( 1 + x h,i ) + σ B u h,i 1 -h ρ f ( 1 -θ ) ,
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: A positive scalar h > 0 exists such that for all x 0 ∈ G(t 0 , •) -1 (C * ) and for all h ∈ (0, h], a unique pair (x h,i+1 , u h,i+1 ) exists satisfying for all i = 0, 1, 2, . . . , N h :

ẋ = f (t, x) + B(t, x)u, x(0) = x 0 , C u ⊥ G(t, x) ∈ C * . ( 8
x h,i+1 -x h,i = h[f (t h,i+1 , θ x h,i + (1 -θ)x h,i+1 ) + B(t h,i , x h,i )u h,i+1 ] C u h+1 ⊥ G(t, x h,i+1 ) ∈ C * . ( 8

.16)

The discrete-time trajectory {(x h,i , u h,i )} is uniformly bounded; moreover, there is a sequence {h ν } ↓ 0 such that x h ν → x uniformly on [0, T] and u h ν u weakly in L 2 (0, T); moreover all such limits ( x, u) are weak solutions of (8.15). Specializing the above result to the DCP (2.9) derived from the VIE: ẋ ∈ -F(x)-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 n and θ ∈ [0, 1] be a given scalar. Let F : n :→ n be a Lipschitz continuous function. A positive scalar h > 0 exists such that for all x 0 ∈ K and for all h ∈ (0, h], a unique pair (x h,i+1 , u h,i+1 ) exists satisfying for all i = 0, 1, 2, . . . , N h :

K * u h,i+1 ⊥ x h,i+1 = I + hF(θx h,i + (1 -θ)•) -1 (x h,i + hu h,i+1 ) ∈ K.
The discrete-time trajectory {(x h,i , u h,i )} is uniformly bounded; moreover, there is a sequence {h ν } ↓ 0 such that x h ν → x uniformly on [0, T] and u h ν → u weakly in L 2 (0, T); finally, any such limit x is a weak solution of the VIE ẋ ∈ -F(x) -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) ≡ -F(x), B(t, x) ≡ I, C ≡ K * , G(t, x) ≡ x. Clearly, any element of C is a solution of the initial CP (C, G(0, x 0 )).

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 ∈ K (for (7.3), and all x ∈ 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.

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:

ẋ = f (t, x) + B(t, x)u, u ∈ SOL(K, G(t, x) + F), b = Mx(0) + Nx(T), (9.1)
where b is an n-vector and M and N are n × 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 ≡ ⎛ ⎜ ⎜ ⎜ ⎝ x h,0 x h,1 . . . x h,N h +1 ⎞ ⎟ ⎟ ⎟ ⎠ ∈ (N h +2)n and u h ≡ ⎛ ⎜ ⎜ ⎜ ⎝ u h,1 u h,2 . . . u h,N h +1 ⎞ ⎟ ⎟ ⎟ ⎠ ∈ K h , (9.2)
where 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, . . . , N h ,

x h,i+1 = x h,i + h[f (t h,i+1 , θ x h,i + (1 -θ)x h,i+1 ) + B(t h,i , x h,i )u h,i+1 ], u h,i+1 ∈ SOL(K, G(t h,i+1 , x h,i+1 ) + F) (9.3)
plus the boundary condition:

b = Mx h,0 + Nx h,N h +1 . ( 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 (x h , u h ) have to be computed simultaneously. Detailed discussion of how to solve such (discretetime 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 (x h , 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 ψ

IVP 1 x
. 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 ⊆ 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 ρ > 0 exists such that solutions to the VI (K, q + F) exist and satisfy the linear growth property (6.5) for all q ∈ G( );

(b) condition (D) holds; (c) M + N is nonsingular and e Tψ IVP 1 x < 1 + 1 ( M + N ) -1 N . ( 9

.5)

There exists h > 0 such that for all h ∈ (0, 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, . . . , N h , the continuous map x i : n → n as follows: x i (x 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 ref ) = (t h,i+1 , t h,i ). The well-definedness of the map x i is justified by Proposition 8.2, and its continuity by an easy argument. Define the composite maps i , for i = 0, 1, 2, . . . , N h + 1, recursively by: 0 ≡ I, and

i+1 ≡ x i • i , i = 0, 1, . . . , N h . From Lemma 7.2 we have, for i = 0, 1, . . . , N h , i+1 (x ref ) ≤ ρ i+1 and i+1 (x ref ) -i (x ref ) ≤ ρ i+1 -ρ i , (9.6)
where ρ i+1 satisfies the recursion

ρ i+1 ≡ ( 1 + h ψ x ) ρ i + h ψ IVP 1 x , i = 0, 1, . . . , N h , (9.7) with ρ 0 ≥ x 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 ≡ ( M + N ) -1 b -( M + N ) -1 N • ( N h +1 -0 ), ( 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 ρ 0 > 0 such that H is a self-map from the closed Euclidean ball ρ 0 B in n with center at the origin and radius ρ 0 into itself. Since

H(x ref ) ≤ (M + N) -1 b + (M + N) -1 N N h +1 (x ref ) -0 (x ref ) ≤ ( M + N ) -1 b + ( M + N ) -1 N ( ρ N h +1 -ρ 0 ),
where the second inequality follows from the recursion (9.6), it follows that if ρ 0 satisfies

( M + N ) -1 b + ( M + N ) -1 N ( ρ N h +1 -ρ 0 ) < ρ 0 , (9.9)
then H(x ref ) < ρ 0 for every

x ref satisfying x ref ≤ ρ 0 ; i.e., H is a self-map as described. Since ρ N h +1 ≤ e T ψ IVP 1 x ρ 0 + e T ψ IVP 1
x -1, (9.9) holds if

( M + N ) -1 b + ( M + N ) -1 N ( e T ψ IVP 1 x -1 ) ( ρ 0 + 1 ) < ρ 0 .
Under condition (9.5), we have 1 -(M + N) -1 N (e Tψ IVP 1

x -1) > 0; therefore, with the choice of

ρ 0 > ( M + N ) -1 b + ( M + N ) -1 N ( e T ψ IVP 1 x -1 ) 1 -( M + N ) -1 N ( e T ψ IVP 1 x -1 ) , ( 9.10)
(9.9) follows.

An important remark about the above proof is that the choice of ρ 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. (i) K contains the origin and (K, D) is an R 0 pair, (ii) K is a polyhedron containing the origin and G( ) ⊆ int K(K, D) * , the following three statements hold:

(a) there exists h > 0 such that for all h ∈ (0, h], a discrete-time boundary-value trajectory (9.2) exists that satisfies (9.3) and (9.4); (b) there is a sequence {h ν } ↓ 0 such that the following limits exist: x h ν → x uniformly on [0, T] and u h ν u weakly in L 2 (0, T); (c) all such limits ( x, u) from part (b) are weak solutions of the boundary-value DVI (9.1).

Proof Let ρ 0 > 0 satisfy (9.10). Assertion (a) follows from Lemma 9.1. Since {x h,0 } is bounded in norm by ρ 0 for all h > 0 sufficiently small, we have the bounds x h,i ≤ ρ i and x h,i+1x h,i ≤ ρ i+1 -ρ i by (9.6) with ρ 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.

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 x i (x ref ) defined in Lemma 9.1 are multi-valued; thus so are the composite maps 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 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, . . . , N h ,

x h,i+1 = x h,i + h [ f (t h,i+1 , θ x h,i + (1 -θ)x h,i+1 ) + B(t h,i , x h,i )u h,i+1 ], u h,i+1 ∈ SOL(K, G(t h,i , x h,i ) + F), b = Mx h,0 + Nx h,N h +1 .
Analytically, a major simplification results from the modified scheme. Namely, corresponding to a given pair (t h,i , x h,i ), 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

x h,i+1 ∈ x i (x h,i ), i = 0, 1, . . . , N h , b = Mx h,0 + Nx h,N h +1 , (9.11)
where x i : n → n is the set-valued map defined by: for

x ref ∈ n , x i (x ref ) ≡ [ I -hf (t h,i+1 , θ x ref + (1 -θ)•) ] -1 ×[x ref + hB(t h,i , x ref )SOL(K, G(t h,i , x ref ) + F) ].
We have noted several times that the map

[I -hf (t h,i+1 , θ x ref + (1 -θ)•)] -1 is
a homeomorphism for all h > 0 sufficiently small, provided that f is Lipschitz continuous. Consequently, for such an h, the map x i is well defined. By way of relating x i to the map S i defined in Proposition 7.1, we note that x i (x h,i ) = S i (x h,i ). The two maps 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 , . . . , x h,N h +1 } are well defined. We do this by applying Górniewicz' fixed-point Theorem 7.3. Define j : 2n → n by j(x, y)

≡ ( M + N ) -1 b -( M + N ) -1 N( y -x ), ( x, y ) ∈ 2n ,
and, for a given h > 0, also the set-valued map Φ : n → n by

Φ(x ref ) ≡ j(x ref , x N h • x N h -1 • • • • • x 1 • x 0 (x ref )), x ref ∈ n .
(9.12)

Clearly, a fixed point x h,0 of Φ yields the desired iterates via the iterative set membership: x h,i+1 ∈ x i (x h,i ), for i = 0, 1, . . . , N h . Although we could formulate the following result for more general boundaryvalue DIs, we refrain from doing so and restrict our attention to the DVI on hand. In the theorem below, each function 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 ⊆ m be a nonempty closed convex set. Let (f , B, G) satisfy conditions (A) and (B). Let F : m → m be continuous. Assume that SOL(K, q + F) is a nonempty convex set for all q ∈ G( ) and that for some constant ρ > 0, elements in SOL(K, q + F) satisfy the linear growth property (6.5) for all q ∈ G( ). Assume further that M + N is nonsingular and (9.5) holds. Then the following two statements are valid:

(a) there exists h > 0 such that for all h ∈ (0, h], a discrete-time boundary-value trajectory {x h,0 , x h,1 , • • • , x h,N h +1 } exists satisfying (9.11); (b) there is a sequence {h ν } ↓ 0 such that x h ν → x uniformly on [0, T]; moreover, any such limit x along with a suitable u(t) ∈ 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 x i is an acyclic map. Proceeding as in Lemma 9.1 and choosing ρ 0 as described therein, we can show that Φ in (9.12) is a set-valued map from ρ 0 B into subsets of ρ 0 B. This map Φ, 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 Φ as a composition of acyclic maps as follows: Let : ρ 0 B → ρ 0 B × ρ 0 B be the diagonal map (z) = (z, z). Then Φ can be identified as the composition

ρ 0 B → ρ 0 B × ρ 0 B id× x 0 → ρ 0 B × ρ 1 B id× x 1 → • • • id× x N h → ρ 0 B × ρ N h +1 B j → 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 Φ is an admissible map. Since ρ 0 B is a compact AR, Φ therefore has a fixed point by Górniewicz' fixed-point theorem.

The proof of the theorem can be easily completed as before.

The generalized LCS

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

ẋ = p + Ax + Bu, C u ⊥ q + Cx + Du ∈ C * , b = Mx(0) + Nx(T), (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, . . . , N h + 1,

x h,i+1 -x h,i = h [ p + θ Ax h,i + ( 1 -θ ) Ax h,i+1 + Bu h,i+1 ], C u h,i+1 ⊥ q + Cx h,i+1 + Du h,i+1 ∈ C * , b = Mx h,0 + Nx h,N h +1 . (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 then u h must be a solution of the CP (C, q h , D h ), where

q h ≡ q + C[ I -h ( 1 -θ ) A ] -1 [ h p + ( I + h θ A ) x ref ].
By Lemma 9.2, it follows that Bu h is a single-valued Lipschitz continuous function of x ref ; hence so is x h , which is xi (x 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 n ⊆ int K(C, D) * , which becomes q + C n ⊆ 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 ψ IVP 2 x , ρ IVP 2 u , and h exist such that for all h ∈ (0, h] and all x h,0 ∈ C -1 (C * ), a unique pair (x h,i+1 , u h,i+1 ) exists satisfying, for all i = 0, 1, . . . , N h ,

x h,i+1 -x h,i = h [ p + θ Ax h,i + ( 1 -θ ) Ax h,i+1 + Bu h,i+1 ], C u h,i+1 ⊥ q + Cx h,i+1 ∈ C * , b = Mx h,0 + Nx h,N h +1 , (9.16)
x h,i+1x h,i ≤ hψ IVP 2 x

(1 + x h,i ) and u h,i+1 ≤ ρ IVP 2 u (1 + 2 x h,i ). Based on these preliminary remarks, we can establish the following result. 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 ≡ ρ 0 B ∩ C -1 (C * ) into itself, where ρ 0 is the constant satisfying (9.10) corresponding to the constant ψ IVP sufficiently small. It follows that for some constant σ > 0 independent of h and x ref , sup{ C T u h : u h ∈ SOL(C, q h , D h ) } ≤ σ ( 1 + x ref ). (9.19) Returning to the above string of expressions, we deduce that we may choose the constant σ sufficiently large so that sup{ v T u h : u h ∈ SOL(C, q h , D h ) } ≤ σ ( 1 + x ref ).

Consider the family of polyhedra with variable right-hand sides:

P(a, α) ≡ { u ∈ C : C T u = a, v T u = α }, for ( a, α ) ∈ n × + .
By the renowned Hoffman bound for linear inequality systems [START_REF] Hoffman | On approximate solutions of systems of linear inequalities[END_REF], it follows that there exists a constant ς > 0, dependent on the triple (C, C, v) only, such that for all pairs (a, α) for which the polyhedron P(a, α) is nonempty, there exists a vector u ∈ P(a, α) such that u ≤ ς( a + α). By [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]Corollary 2.3.13]

SOL(C, q h , D h ) = { u ∈ C : C T u = C T u h , v T u = v T u h }, (9.20)
for any u h ∈ SOL(C, q h , D h ), it follows that for all h > 0 sufficiently small and all x ref ∈ n , u h ∈ SOL(C, q h , D h ) exists such that u h ≤ 2σ ς(1 + x ref ).

From this bound, it is easy to derive a corresponding bound for x hx ref .

The proof that C T u h is a single-valued Lipschitz function of x ref is very similar to the proof of (9.19). Indeed, for any two vectors q + Cr i for i = 1, 2, with u i ∈ SOL(C, q + Cr i , D h ), we have C T u 1 -C T u 2 ≤ σ -1 r 1r 2 . The last assertion of the lemma follows.

Remark The CP (9.18) may have multiple solutions that are possibly unbounded; yet, one exists with the much needed linear growth property.

Based on the above constant ξ , we can define the constants ψ x > 0 and ρ 0 satisfying (9.10). Moreover, for each h > 0 sufficiently small, we can define the scalars ρ 1 , . . . , ρ N h +1 recursively by (9.7). The rest of the proof of the following result does not need to be repeated. Theorem 9.5 Let C be a polyhedral cone. Suppose that q ∈ C n + C * . Assume further that M + N is nonsingular, (9.5) holds, and C(M + N) -1 (b + N n ) ⊆ C * . A discrete-time BVP trajectory satisfying, for all i = 0, 1, . . . , N h , x h,i+1x h,i = h [ p + θ Ax h,i + ( 1 -θ ) Ax h,i+1 + C T u h,i+1 ], C u h,i+1 ⊥ q + Cx h,i+1 ∈ C * , b = Mx h,0 + Nx h,N h +1 , exists; moreover, a subsequence of this trajectory converges to a weak solution of the boundary-value generalized LCS (9.17).

The linear-quadratic Nash game revisited

A major effort in the treatment of the boundary-value DVI in the previous section has been placed on the demonstration that the discrete-time subproblems (9.3) and (9.11) have solutions. Such a treatment leads to a condition such as (9.5) that restricts the terminal time T based on several model constants.

Once a discrete-time trajectory is known to exist, the convergence of such a trajectory to a continuous-time solution of the boundary-value DVI (9.1) turns out to be not as hard to show; it fact, it can be accomplished rather simply with the aid of the IVP analysis. The linear-quadratic differential Nash game discussed in Sect. 4.2 is an example of such a problem. In what follows, we describe a time-stepping scheme for solving this game for which the solvability of the discrete-time subproblems follows from a known result for a noncooperative Nash game in finite-dimensions (thus, no advanced fixed-point theorem is needed); more importantly, the condition (9.5) can be removed completely.

To simplify the discussion, we focus on the linear-quadratic game and assume that the data functions d(t), q(t), B(t), C(t), and S(t) are all time invariant. The DVI formulation for the resulting linear-quadratic Nash game is a boundaryvalue DAVI:

λ ẋ = d + A λ x + Bu, u(t) ∈ SOL(U, q + Cz, D)
with very special boundary conditions: x(0) = x 0 and λ(T) = Wx(T) + c. The time-stepping scheme computes, for i = 0, 1, 2, . . . , N h + 1, the vectors x h,i ≡ (x ν,h,i ) ν∈N ∈ n , u h,i ≡ (u ν,h,i ) ν∈N ∈ m , and λ h,i ≡ (λ ν,h,i ) ν∈N ∈ n which are approximations of the state trajectory x(t), the control trajectory u(t), and the adjoint trajectory λ(t), respectively, at times t = t h,0 , t h,1 , t h,2 , . . . , t h,N h +1 .

Respecting the initial condition x h,0 = x 0 and the terminal condition λ(T) = c + Wx(T), (the latter suggests that λ should be treated backward in time, which leads to the term λ ν,h,i in (c) and the right-hand side of (a)), the discrete-time system is as follows: for every ν ∈ N and for all i = 0, 1, . . . , N h , (a) λ ν,h,i+1 -λ ν,h,i = -h p ν + ν ∈N P νν x ν ,h,i+1 + Q νν u ν ,h,i+1

+P νν x ν,h,i+1 + (R νν ) T u ν,h,i+1 + (G ν ) T λ ν,h,i (b) x ν,h,i+1x ν,h,i = h r ν + G ν x ν,h,i+1 + H ν u ν,h,i+1 (c) u ν,h,i+1 ∈ U ν satisfies: for all u ν ∈ U ν , u νu ν,h,i+1 T q ν + ( H ν ) T λ ν,h,i + ( Q νν ) T x ν,h,i+1 + S νν u ν,h,i+1 ν ∈N R νν x ν ,h,i+1 + S νν u ν ,h,i+1 ≥ 0, (d) λ ν,h,N h +1 = c ν + W νν x ν,h,N h +1 + ν ∈N W νν x ν ,h,N h +1 , (e) x ν,h,0 = x ν,0 .

We claim that the pair (x h , u h ), where x h ≡ (x h,i )

N h +1 i=1
and u h ≡ (u h,i )

N h +1
i=1 , constitutes an equilibrium to a discrete-time Nash game, with λ h ≡ (λ h,j ) N h j=0 being the vector of KKT multipliers of the players' dynamic constraints in discrete time. The latter static game has the same number of players; player ν's optimization problem is in terms of the variables (x ν,h,i , u ν,h,i )

N h +1
i=1 , with the rival players' variables (x -ν,h,i , u -ν,h,i ) N h +1 i=1 taken as exogenous:

minimize θ ν (x h , u h ) ≡ (x ν,h,N h +1 ) T c ν + ν ∈N W νν x ν ,h,N h +1 + h N h i=0 (x ν,h,i+1 ) T p ν + ν ∈N P νν x ν ,h,i+1 + Q νν u ν ,h,i+1
+ (u ν,h,i+1 ) T q ν + ν ∈N R νν x ν ,h,i+1 + S νν u ν ,h,i+1 (10.1) subject to x ν,h,0 = x ν,0 and for i = 0, 1, . . . , N h ,

x ν,h,i+1x ν,h,i = h r ν + G ν x ν,h,i+1 + H ν u ν,h,i+1 , u ν,h,i+1 ∈ U ν .

When each U ν is a polytope, Proposition 10.1 establishes the above claim, from which the existence of a tuple (x h , u h , λ h ) satisfying the discrete-time system (a)-(e) follows readily.

Proposition 10.1 Let each U ν be a nonempty compact polyhedron in m ν . Assume that for all ν ∈ N , W νν and

ν ≡ 2P νν Q νν + ( R νν ) T ( Q νν ) T + R νν 2S νν (10.2)
are symmetric positive semidefinite. A scalar h > 0 exists such for all h ∈ (0, h], a tuple (x h , u h , λ h ) exists satisfying the discrete-time system (a)-(e).

Proof For fixed (x -ν,h,i , u -ν,h,i )

N h +1
i=1 , (10.1) is a convex quadratic program in player ν's variables (x ν,h,i , u ν,h,i )

N h +1
i=1 , which can be reformulated in terms of (u ν,h,i )

N h +1 i=1
only by solving for x ν,h,i+1 from the difference equation in the constraints of (10.1). Indeed, for all h > 0 sufficiently small, we have x ν,h,i+1 = ( Ih G ν ) -1 x ν,h,i + h ( r ν + H ν u ν,h,i+1 ) , i = 0, 1, . . . , N h .

Lemma 6 . 1 [ 32 ,

 6132 Theorem 5.1] Let F :

Proposition 6 . 4

 64 Let K ⊆ m be a nonempty closed convex set containing the origin. Let D be an m × m positive semidefinite matrix. For every vector r ∈ int K(K, D) * , SOL(K, r, D) is a nonempty convex set on which (D + D T )u is a constant for all u ∈ SOL(K, r, D). Furthermore, (a) if (K, D) is an R 0 pair, then a constant ρ > 0 exists such that for all r ∈ m , SOL(K, r, D) = ∅ and sup{ u : u ∈ SOL(K, r, D) } ≤ ρ ( 1 + r ). (6.8) (b) if K is a polyhedron, then (6.8) holds for all r ∈ int K(K, D) * ; moreover, the (single-valued) map r → (D + D T )SOL(K, r, D) is Lipschitz continuous on K(K, D); i.e., there exists a constant L D > 0 such that

  part (b) of Theorem 2.5.15], it follows that SOL(K, r, D) = ∅ for all r ∈ K(K, D) * . Noticing that K(K, D) * is a polyhedral 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 ≡ 2n + ,

Lemma 6 . 4

 64 Let M ∈ m×m . The following three statements are equivalent: (a) There exists c > 0 such that u T Mu ≥ c Mu 2 for all u in m . (b) M = E T AE for some matrices E ∈ ×m and A ∈ × for some positive integer , where A is positive definite. (c) M is psd-plus.

Theorem 6 . 1

 61 Let K ⊆ m be a nonempty closed convex set. Let (f , B, G) satisfy conditions (A) and (B). Let F : m → 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 ref ∈ K satisfying (6.6); (b) F ≡ E T • • E for some matrix E ∈ ×m satisfying K ∞ ∩ ker E = {0} and for some continuous map : → 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 R 0 pair; (d) K is a polyhedron containing the origin, F(u) ≡ Du for some positive semidefinite matrix D, and G( ) ⊆ int K(K, D) * ; (e) K is a polyhedron containing either the origin or no lines, and F ≡ D + Φ,

1 .Theorem 7 . 4

 174 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. Let K ⊆ 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( ) ⊆ int K(K, D) * ; (d) K is a polyhedron containing either the origin or no line and F(u) ≡ Du + ϒ(Du), where ϒ : m → m is Lipschitz continuous such that with Φ(u) ≡ ϒ(Du), the pair (K, F) satisfies condition (e) in Theorem 6.1;

. 15 ) 8 . 1

 1581 Corollary Let C be a closed convex cone in m and θ ∈ [0, 1] be a given scalar. Let (f , B, G) satisfy the conditions (A), (B) and (D).

Theorem 9 . 1

 91 Let K ⊆ 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:

Theorem 9 . 4

 94 Let C be a polyhedral cone in m . Assume that CB is positive definite. Assume further that M + N is nonsingular, (9.5) holds, and C(M + N) -1 (b + N n ) ⊆ 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.

2

 2

x

  noted above. For a vector x ref ∈ S, the proof of Lemma 9.1 shows that H(x ref ) ∈ ρ 0 B. Since H(x ref ) ∈ (M + N) -1 (b + N n ), the assumption implies that H(x 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 conditionC(M + N) -1 (b + N n ) ⊆ C * becomes Cx 0 ∈ C * ,which is exactly the special case of the familiar condition x 0 ∈ G(t 0 , •) -1 (C * ) in Corollary 8.1 that corresponds to an initial-value generalized LCS with D = 0.

  .[START_REF] Anitescu | A constraint-stabilized time-stepping for multi-body dynamics with contact and friction[END_REF] where f : n → n is a Lipschitz function, K is a polyhedron in 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: ẋ = 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[START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF] Exercises 5.6.14] and also[START_REF] Billups | Solutions to affine generalized equations using proximal mappings[END_REF] Lemma 3.3]. The domain of a multifunction (i.e., set-valued function) Φ,

denoted dom Φ, is the set of vectors x such that Φ(x) = ∅; a multifunction is polyhedral if its graph, i.e., the set {(x, y) : y ∈ Φ(x)}, is the union of finitely many polyhedra. Lemma 5.1 (Gowda) A single-valued polyhedral multifunction Φ is piecewise linear on its domain; i.e., Φ is continuous on dom Φ and there exist finitely many linear

  This together with the equation (5.10) yields C α• x * (t) + D αα u * α (t) + D αβ u * β (t) = 0, by the nonsingularity of P αα . By restricting δ if necessary, we have C γ • x * (t) + D γ α u * α (t) + D γβ u * β (t) > 0 on [0, δ], by continuity of the solution trajectory. Consequently, by letting u * γ (t) ≡ 0 on [0, δ], we obtain a desired locally smooth solution of the original LCS (2.8).

  .10) Proof It suffices to prove(8.10). WriteP T u h ≡ (µ h , λ h ). Since K contains the G(t ref , x ref ) + E T (Eu ref ) ) T Zµ h = ( Zµ h ) T G(t ref , x ref ).

	origin, Lemma 8.1 implies that P	µ h 0	= Zµ h belongs to K, Therefore,
	0 ≤ ( (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

Acknowledgments

The authors are very grateful to two referees and the Associate Editor for their many insightful and constructive comments. Both referees have identified numerous technical details that have helped to clarify the proofs and improve the presentation of the paper; their diligence and great care in reviewing our lengthy work is greatly appreciated. One referee has kindly reminded us of related references on the VIE and stressed that the class of DVIs treated in the paper is of index no higher than 2. As seen in Sect. 2.4, the VIE is a particular realization of the DVI. The two problems can be unified by the common formulation (2.12), which was observed

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. The work of D. E. Stewart is supported by the National Science Foundation under a Focused Research Group grant DMS-0138708. J.-S. Pang (B)

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

(D) there exists a constant η G > 0 such that

for all r ∈ n , u and u in m , t and t ref in [0, T], and x ref ∈ n .

The liminf condition in (C) obviously holds if F is monotone on K. Specialized to the LCS (2.8) where K = m + , condition (C) holds if the matrix D is copositive (i.e., u T Du ≥ 0 for all u ≥ 0); whereas condition (D) holds if and only if the matrix CB is positive definite. In particular, D is not required to be an R 0 matrix, as required by conditions (b) and (d) in Theorem 7. [START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF].

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 ref , x ref ) is uniformly positive definite for all tuples (t, t ref , x, x 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

Proof Suppose (D) holds. Let u ≡ u + τ v for τ > 0. We then have

Hence, for every vector u, there exists a unique vector x(u) satisfying

We claim that the map u → G(t, x(u)) is strongly monotone on m , provided that h satisfies (8.3). Once this is establish, the existence of the pair (x h , u h ) satisfying (8.4) follows readily (see [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF]Proposition 2.2.7]), because the map

Moreover, if F is monotone on K, then the map u → 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

. By (8.5) and (7.4), we have

Moreover, by assumption (D),

Consequently,

By (8.3), it follows that the map u → G(t, x(u)) is strongly monotone.

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 The most interesting case of the above assumption is when ker E = {0}. Indeed, when the latter kernel is a singleton, then F is strongly monotone on m , and this case is covered by Theorem 7.4(a). Ideally, we would like to make only the additional assumption that K ∞ ∩ 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 E and E ⊥ denote, respectively, the null space of E and its orthogonal complement. Let K 1 and K 2 denote the orthogonal projection of K onto E and E ⊥ , respectively. We clearly have K ⊆ K 1 ⊕ K 2 , where ⊕ denotes the orthogonal sum. We postulate the following assumption:

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 ∈ m×k be any matrix whose columns form an orthonormal basis of E so that EZ = 0. Extend the basis Z through E ⊥ to an orthonormal basis P of m . Thus P = [Z W] ∈ m×m for some matrix W ∈ m×(m-k) whose columns form an orthonormal basis of E ⊥ . (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 W ≡ EW ∈ ×(m-k) must have full column rank; moreover, E m = W m-k . In terms of the matrices Z and W, assumption (E) is equivalent to

This follows because ZZ T and WW T are the Euclidean projectors onto E and E ⊥ , 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

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 (µ i , λ i ) satisfy the left-hand condition in (8.6).

We then have Zµ i + Wλ i ∈ K. By (E), it follows that

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

substituting this into the CP yields

where

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 [START_REF] Çamlibel | Complementarity methods in the analysis of piecewise linear dynamical systems[END_REF][START_REF] Heemels | Linear complementarity systems: a study in hybrid dynamics[END_REF]. 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

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 0 is defined to be 1):

At first sight, condition (ABCD-2b) seems fairly strong because it requires in particular that CA i B be positive semidefinite for all i (when θ < 1); nevertheless, unlike (ABCD-2a), condition (ABCD-2b) does not require any of the CA i B 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 h is 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 i B = 0 for all i ≥ 1 and

If u T CBu = 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 m . Under assumption (ABCD), the matrix D h is positive semidefinite for all h > 0 sufficiently small; moreover, for all r ∈ K(C, D h ) * , Bu is constant for all u ∈ SOL(C, r, D h ); hence the map B r : r ∈ K(C, D h ) * → 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 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 m . Assume condition (ABCD) and that q + C n ⊆ 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 xi (x ref ) remains single-valued. This is because if (x h , u h ) is any pair, which must exist, satisfying

Next we consider the following class of boundary-value generalized LCSs:

where 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 C being the nonnegative orthant, problem (9.17) provides an equivalent formulation for the VIE: ẋ ∈ p + Ax -N (x; K), where K is the polyhedron {x ∈ n : q + Cx ≥ 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

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 C be a polyhedral cone. Suppose that q ∈ C n + C * . There exist constants h > 0 and ξ > 0 such that for all h ∈ (0, h] and all x ref ∈ n , a pair (x h , u h ) exists satisfying

and max(

Lipschitz function of x ref with a constant that is independent of h > 0 sufficiently small.

Proof Let h > 0 be sufficiently small so that Ih(1 -θ)A is positive definite. Since (C, D h ) is an affine pair and q belongs to C n + C * , it follows that for all x ref ∈ n , the generalized LCP:

x ref ] has a solution. To establish the existence of the constant ξ , write q ≡ Cs + v where s ∈ n and v ∈ C * . We have, for any u h ∈ SOL(C, q h , D h ),

where σ is any positive constant not exceeding the smallest eigenvalue of the symmetric part of the positive definite matrix [Ih(1 -θ)A] -1 , for all h > 0 Since x ν,h,0 is known, it follows by induction that the variables (x ν,h,i )

i=1 can be expressed as a linear function of (u ν,h,i ) N h +1 i=1 . Hence, after substituting the former variables into the objective function of (10.1), we obtain a quadratic program in the variables (u ν,h,i ) N h +1 i=1 , which remains convex and has an optimal solution, by the compactness of each U ν . Therefore, by a well-known existence result for a standard, static Nash game with compact convex strategy sets, such as [38, Proposition 2.2.9], it follows that the discrete-time Nash problem has an equilibrium solution (x ν,h,i , u ν,h,i )

. Letting λ ν,h,i be the KKT multiplier of the constraint:

we can readily show that the tuple (x h , u h , λ h ) satisfies the desired discrete-time system (a)-(e).

Since U ν is bounded and u ν,h,i+1 ∈ U ν , it follows that the assumptions of Lemma 7.2 are readily satisfied. The convergence of the above time-stepping method for solving the linear-quadratic differential Nash game therefore follows from Theorem 7.1, provided that in addition the matrix D is positive semidefinite. The details are omitted.

Concluding remarks

In this paper, we have undertaken an extensive study on the DVI and several of its special cases, including the generalized LCS. These problems have the general form of a differential equation where a term in the right-hand side is a solution of a VI/CP. Depending on the structure of the VI and the differential equation, the regularity of the solutions and their theory can vary considerably. In the simplest case, the VI has a (locally or globally) unique solution that depends in a Lipschitz fashion on the data of the VI. This leads to Lipschitz ODEs for which there is a well-known theory. If, however, there are many solutions then we can use the theory of differential inclusions. Our theory goes beyond these cases to situations where the VI does not always have solutions, and where the solutions can be unbounded. Then our existence theory depends on the joint structure of the differential equation and the VI. To prove these results we often rely on the use of basic time-stepping schemes, whose extension to boundary-value problems is also analyzed. There remain many issues that require further investigation. These include the efficient numerical solution of the discrete-time subproblems and the application of our theory to engineering and economic problems in which dynamics is an important concern.

by Professor Hans Schumacher of Tilburg University, who has made other valuable comments that we have incorporated in the revision of the paper. We owe these individuals our sincere gratitude. Finally, we thank Professors Seetharama Gowda of the University of Maryland Baltimore County and Michael Ferris of the University of Wisconsin at Madison for providing us with the references [START_REF] Billups | Solutions to affine generalized equations using proximal mappings[END_REF] and [START_REF] Ban | Quasi-variational inequality formulations and solution approaches for dynamic user equilibria[END_REF], respectively.