The Euler scheme for state constrained ordinary differential inclusions

We propose and analyze a variation of the Euler scheme for state constrained ordinary differential inclusions under weak assumptions on the right-hand side and the state constraints. Convergence results are given for the space-continuous and the space-discrete versions of this scheme, and a numerical example illustrates in which sense these limits have to be interpreted.


Introduction
Numerical methods for ordinary differential inclusions have mainly been considered in the unconstrained case and for convex-valued right-hand sides. Summaries of early results can be found in the surveys [10] and [14]. An approach using subdivision and continuation techniques for time-independent systems with affine controls was given in [18].
The variation of the Euler scheme proposed in [12] only uses extremal points of the right-hand side, so that a fully discretized scheme is obtained for polytope-valued right-hand sides. The high complexity of this method motivated the detailed analysis of spatial discretization effects on the Euler scheme given in [7]. A further reduction of computational costs can be achieved by tracking the boundaries of the reachable sets of the fully discretized Euler scheme using only lower-dimensional data of the right-hand side, see [16].
A very different type of algorithm has been proposed in [4], where a reachable set is computed by an optimal control routine that minimizes the distance to grid points in the relevant region at the terminal time. As the optimizer may find a local instead of a global minimum, the propriety of this method cannot be guaranteed.
An implicit Euler scheme has been analyzed in [8], which is considerably more efficient than the explicit Euler scheme when applied to stiff differential inclusions. The disadvantage of this method is that it has to solve an algebraic inclusion in every time step. This problem is successfully avoided by the semi-implicit Euler schemes discussed in [15], leading to a better performance of the algorithm.
The Euler scheme for differential inclusions with nonconvex right-hand sides has been investigated in [17]. Error estimates for the Euler scheme applied to a differential inclusion with convex, but one-sided Lipschitz righthand side have been published in [9].
Up to our knowledge, the only convergence analysis of the Euler scheme for ordinary differential inclusions with state constraints has been given in the paper [3]. The analysis is based on stability theorems that quantify the relationship between the solution set of the constrained problem and the same differential inclusion without constraints. These results play an important role in optimal control theory and are still a subject of intensive research, see [5], [6] and the references therein.
The drawback of these stability theorems in the present context are the strict assumptions on the compatibility between right-hand side and constraints that are needed to obtain such a result. In applications like collision avoidance, see [11], [13] and the references therein, such assumptions do certainly not hold.
The present article is motivated by the lack of a convergence analysis for the Euler scheme in the present of state constraints under weaker conditions. In Section 2 we fix the notation, give a formal problem statement and formulate our basic assumptions, that the initial set is compact, the right-hand side is Lipschitz, the state constraints are upper semicontinuous in time and that all data are compatible in the very weak sense that there is a nontrivial time interval on which a solution of the state constrained inclusion exists.
Under these conditions, we prove some auxiliary results in Section 3, which are used in Section 4 to show that the solution sets of the inflated Euler scheme converge towards the exact solution sets in Hausdorff distance. In Section 5, we show that the same holds for a spatially discretized version of the scheme provided all constraints are relaxed in a proper way. A numerical example given in Section 6 illustrates in which sense this limit has to be interpreted. Note that convergence results for the inflated Euler scheme without spatial discretization are of interest for first-discretize-then-optimize strategies in optimal control and algorithms such as the one proposed in [4], which consider the Euler scheme as a constraint of a finite-dimensional optimization problem.
At the end of the paper, in Section 7, we give a brief explanation how stability results in the spirit of [5] and [6] automatically yield linear convergence of the inflated Euler scheme in Hausdorff distance.

Setting and notation
Throughout this paper, the Euclidean norm and the modulus will be denoted by · : Ê d → Ê + and | · | : Ê → Ê + . As some of the following concepts will be used not only in the state space Ê d , but also in the space of continuous functions, we introduce them in a general Banach space E.
Let A, B ⊂ E be compact. Then is called the Hausdorff semidistance between A and B. In the particular case, where A = {a} with a ∈ E, the Hausdorff semidistance is the usual distance dist(a, A) from a point to a set, and we denote is a metric on the compact subsets of E. The metric projection from x ∈ E to A is the set Since A is compact, the metric projection is nonempty. If, in addition, E is a Hilbert space and A is convex, then proj(x, A) is a singleton.
, we denote the Banach space of absolutely continuous functions from [0, τ ] to Ê d , i.e. the space of all functions, which possess a weak derivative in L 1 ((0, τ ), Ê d ).
Let F : [0, T ] × Ê d ⇒ Ê d be a set-valued mapping. A solution to the ordinary differential inclusionẋ is an absolutely continuous function x(·) ∈ W 1,1 ([0, T ], Ê d ) satisfying (1) almost everywhere in [0, T ]. Depending on the application, solutions may be required to satisfy an initial condition with initial set X 0 ⊂ Ê d and state constraints is called a solution of the inflated Euler scheme with step-size h N , which may as well be subject to a possibly relaxed initial condition and possibly relaxed state constraints dist(y N,n , A(t N,n )) ≤ δ N .
For ν ∈ {0, . . . , N}, we denote the solution sets of the discrete unconstrained and the discrete constrained problems with time horizon ν by Throughout this paper, we posit the following assumptions on the ordinary differential inclusion and the constraints.
Hypothesis A. The right-hand side F , the constraints A and the initial set X 0 satisfy the following conditions:

convex and compact values and
(2) The state constraints A : [0, T ] ⇒ Ê d have closed values and are upper semicontinuous in the sense that t = lim n→∞ t n , x n ∈ A(t n ) and x = lim n→∞ x n imply x ∈ A(t).
(3) The set X 0 ⊂ Ê d is compact, and the data F , A and X 0 are compatible in the sense that there exists τ ∈ (0, T ] with S c (τ ) = ∅.
We will also posit assumptions on the blowup size β N of the Euler scheme and the relaxation parameter δ N for the initial condition and the constraints.
Hypothesis B. The parameters β N and δ N satisfy the following conditions: The relaxation parameter δ N will not be needed before Section 5. We introduce it at this early stage, because we wish to provide a-priori bounds in Lemma 1, which apply to all exact and numerical trajectories in this paper.

Auxiliary results
, the bounds given below hold for the constrained problems as well.
for all x(·) ∈ S u (τ ) and Proof. The bounds for S u (t) and S u N (ν, β N , δ N ) follow by applying Gronwall's Lemma and induction, respectively, in a straight-forward manner. Boundedness of F is a consequence of hypothesis (A1).
The following lemma quantifies the error between exact trajectories and Euler trajectories with perturbed initial value. It will be used to determine an appropriate size β N of the blowup for the inflated Euler scheme.
Proof. As the metric projection to a compact and convex set is continuous according to [2,Theorem 9.3.4], we may define and v ∈ F (t, z) holds by [19,Theorem I.6.13]. Because of Lemma 1, we have so that The following result, which is proved by compactness type arguments, is the core of this paper. In the current situation, it replaces Filippov's theorem (see [ Then there exists some x(·) ∈ S c (τ ) and functionsx N ( along a subsequence.
As a consequence, we obtain the following statement about the life span of the solutions to the the state constrained differential inclusion. Proof. By hypothesis (A3), there exists τ ∈ (0, T ] with S c (τ ) = ∅, so that τ * is well-defined. By definition of the supremum, there exist x N (·) ∈ S c (τ N ), N ∈ AE, such that τ N → τ * . According to Lemma 3, they possess a limit, which is a member of S c (τ * ).

The inflated Euler scheme
The following example shows that it is necessary to modify the Euler scheme to treat state constrained problems.
Motivated by Lemma 2 and the above example, we choose a blowup size In this context, there is no need to relax the initial condition or the state constraints, so that we have δ N = 0. A relaxation will become necessary in Section 5, where not only the time interval, but also the phase space is discretized.
The following proposition states that with this choice of β N , the exact dynamics are overapproximated by the inflated Euler scheme. In particular, Euler trajectories exist at least as long as exact trajectories.
Proof. The existence of a sequence (y N ) N with (14) can be obtained by successively applying Lemma 2 and taking (13) into account.
The statement of Proposition 6 can be reformulated in terms of the Hausdorff semidistance as The following proposition guarantees that the opposite Hausdorff semidistance between numerical and exact trajectories converges to zero. In particular, Euler trajectories do not live longer asymptotically than exact trajectories. Assume that statement (15) is false. Then there exist ε > 0 and a sequence y N ∈ S c The linearly interpolated Euler trajectories y N (·) satisfy the assumptions of Lemma 3, because we have t N,⌊τ /h N ⌋ → τ , the estimate ẏ N (s) ≤ P ∀s ∈ [0, t N,⌊τ /h N ⌋ ] holds according to Lemma 1, we have dist((s, y N (s),ẏ N (s)), graph(F )) by Hypothesis (B1), and the initial condition as well as the state constraints are satisfied by definition. Hence, by Lemma 3, there exist x(·) ∈ S c (τ ), a subsequence AE ′′ ⊂ AE ′ and functionsỹ N ( This contradicts assumption (16), and hence statement (15) is correct.
The second assertion is proved in a very similar way. Assume that τ > τ * and that there exist y N (·) ∈ S c N (⌊τ /h N ⌋, β N , 0), N ∈ AE ′ , for a subsequence AE ′ ⊂ AE. As above, we can extract another subsequence converging to a solution x(·) ∈ S c (τ ), which contradicts S c (τ ) = ∅.

Spatial discretization
When the solution sets are computed in practice, it is necessary to discretize the underlying state space. Simple pathological examples show that in this situation, it is necessary to relax the state constraints in order to obtain convergence of the numerical scheme. The results obtained in this context as well as their proofs are similar to those in Section 4.
In addition to Hypotheses A and B, we posit the following assumptions on the discretization of the state space.
Hypothesis C. The spatial nets ∆ N ⊂ Ê d satisfy the following conditions: (1) For every x ∈ Ê d , there exists some y # ∈ ∆ N with x − y # ≤ δ N .

Consider the discrete solution sets
The following computations show that β # N > 0 is a good choice for the blowup parameter. Hypothesis (C2) strengthens Hypothesis (B2) and ensures that β # N satisfies Hypothesis (B1).
The next proposition implies that the fully discrete Euler trajectories are as close to the exact ones as the grid size allows. In particular, they do not cease to exist before the exact solutions violate the constraints.
Proof. We construct the desired trajectory recursively. By Hypothesis (C1), there exists an element which satisfies the relaxed initial condition (5) as well as the relaxed state constraints (6) and statement (17). Assume that y # N,0 , . . . , y # N,n with (y # N,k ) n k=0 ∈ S c,# N (n, β # N , δ N ) and have been constructed and that n < ν. Then Lemma 2 yields By Hypothesis (C1), there exists a point satisfying the relaxed state constraint (6). Because of (18), we have so that y # N,0 , . . . , y # N,n+1 is a trajectory of the fully discrete Euler scheme.
The following proposition shows that the opposite Hausdorff semi-distance between fully discrete and exact trajectories converges to zero. In particular, fully discrete Euler trajectories do not exist longer asymptotically than exact trajectories.
Proposition 9. For any τ ∈ (0, τ * ], we have Proof. The proof is the same as that of Proposition 7, provided the assumptions of Lemma 3 can be verified in the present setting. But for any sequence and by the relaxed initial condition and state constraints, we have so that Lemma 3 is applicable.

Numerical example
We continue the simple, but instructive Example 5 from Section 4. As it is difficult to visualize sets of trajectories, Figures 1, 2 and 3 show the reachable sets In the figures, the blue semiarc depicts the admissible set which is constant in time. The black line drawn around the semiarc represents the boundary of the admissible set A + B δ N (0) of the numerical scheme. Blue dots are elements of the discrete reachable sets R c,# N (ν, β # N , δ N ), and the centers of the large red dots indicate the position of the exact solution at time t = t N,ν . Whenever dots exceed their corresponding admissible sets, this is due to limited plotting precision.
The step-sizes have been chosen in such a way that we can observe the typical behavior of the inflated Euler scheme, as predicted by the analysis in the sections above. The exact solution is always approximated by a numerical trajectory up to the mesh-size of the spatial grid, which implies that some approximate trajectories live at least as long as the exact solutions. The approximate reachable sets, however, can be substantially larger than that of the exact dynamics. Moreover, numerical trajectories may live much longer than their exact counterparts. Figure 1 shows an extreme case in which the reachable set of the numerical scheme becomes stationary after time t = 1.5, and hence some trajectories are defined on Ê + , while the exact solution ceases to exist at time τ * = π. In Figures 2 and 3, we see that the discrete reachable sets converge towards the true solution as h N → 0 and that the life span τ * N := max{t N,ν : S c,# N (ν, β # N , δ N ) = ∅} of the numerical scheme converges to the exact value τ * = π as h N → 0. This behavior is clearly recognizable in the data given in Figures 4 and 5, where the inflation and relaxation parameters β # N and δ N as well as the life span τ * N are plotted as functions of h N .

Exploiting stability theorems
Many papers in optimal control theory are concerned with the relationship between the solution sets S u (τ ) and S c (τ ) for τ ∈ [0, T ]. The specification of conditions, under which property (D1) below holds, is subject of ongoing research. An overview over the relevant literature as well as instructive counterexamples can be found in [5]. It is certain that conditions, which ensure property (D1), are not satisfied by the class of examples we are aiming at, because property (D1) implies S c (τ ) = ∅ for arbitrarily large τ . Nevertheless, in this section, we will assume this property as given and describe briefly its consequences for the inflated Euler scheme.
The ordinary Euler scheme (without inflation) has been analyzed in this context in [3], so that similar results should also be made available for the inflated scheme. We would also like to provide a template, how property (D1) can be exploited to obtain short and simple proofs for linear convergence of a numerical scheme.
To this end, we suppose the following hypothesis in addition to Hypotheses A and B.
Hypothesis D. The solution sets S u (τ ) and S c (τ ) enjoy a linear stability property, and the constraints are Lipschitz: (1) For any τ ∈ (0, τ * ], there exists some K = K τ > 0 such that for any dist(x u (t), A(t)).
(2) The state constraints are L A -Lipschitz, i.e. there exists As Proposition 6 yields that the first Hausdorff semidistance between the exact and the numerical solution sets is zero, we only need to take care of the second semidistance. The proof of Proposition 10 is remarkably short compared to the machinery set up in [3] to prove a similar estimate under similar assumptions.
It is not difficult to see that the above result persists under a spatial discretization with parameter δ N = h 2 N . In this case, the first Hausdorff semidistance is still covered by Proposition 8.

Conclusion
The aim of the present paper was to prove convergence of the inflated Euler scheme in the presence of state constraints and under weak assumptions, not excluding applications such as the collision avoidance problems discussed in [11] and [13]. This aim is achieved, but the convergence statements are not as strong as one might wish, because we only showed pure convergence. It would be desirable to know if there were any conditions on the right-hand side and the state constraints that are substantially weaker than those imposed in [3], [5] and [6], but strong enough to specify a rate or a speed of convergence of some Euler-like scheme. In that sense, we leave a gap in the literature, which we are not able to fill at the moment.