Maurizio Falcone
email: falcone@mat.uniroma1.it

Dante Kalise
email: dante.kalise@oeaw.ac.at

A high-order semi-Lagrangian/finite volume scheme for Hamilton-Jacobi-Isaacs equations

Keywords: Hamilton-Jacobi-Isaacs equations, high-order schemes, semi-Lagrangian schemes, Finite Volume Methods, optimal control, differential games

We present a numerical scheme for the approximation of Hamilton-Jacobi-Isaacs equations related to optimal control problems and differential games. In the first case, the Hamiltonian is convex with respect to the gradient of the solution, whereas the second case corresponds to a non convex (minmax) operator. We introduce a scheme based on the combination of semi-Lagrangian time discretization with a high-order finite volume spatial reconstruction. The high-order character of the scheme provides an efficient way towards accurate approximations with coarse grids. We assess the performance of the scheme with a set of problems arising in minimum time optimal control and pursuit-evasion games.

Introduction

The numerical approximation of Hamilton-Jacobi-Isaacs (henceforth HJI) equations appears as a crucial step in many fields of applications, including optimal control, image processing, fluid dynamics, robotics and geophysics. In general, these equations do not have have regular solutions even if the data and the coefficients are regular, and therefore many efforts have been devoted to the development and the analysis of approximation schemes for such problems. The convergence of the schemes is understood in the sense of viscosity solutions; it is well known (see e.g. [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Falcone | Semi-Lagrangian schemes for linear and Hamilton-Jacobi equations[END_REF]) that viscosity solutions are typically Lipschitz continuous, and therefore the main difficulty is to have a good resolution around the singularities, and a good accuracy in the parts of the domain where the solution is regular. The theory of approximation schemes for viscosity solutions has been developed starting from the huge literature existing for the numerical solution of conservation laws in one dimension, exploiting the relation between entropy solutions and viscosity solutions. More precisely, the viscosity solution can be written as the space integral of the corresponding entropy solution and this relation can be also applied to the construction of numerical schemes, by simply integrating in space the schemes for conservation laws. At the very beginning, these techniques were successfully applied to the study of the class of monotone schemes; in this framework, the rate of convergence is limited to first order. Some of these schemes, like finite differences for instance, are used over structured grids and are strictly related to the above mentioned methods for conservation laws. Other approximation schemes, like the Finite Volume Method and semi-Lagrangian schemes, can easily work on unstructured grids and are based on different ideas, e.g. on the Hopf-Lax representation formula. In all these cases, the role of monotonicity is important to guarantee the convergence to the viscosity solution, and a general result for monotone schemes applied to second order fully nonlinear equations has been proved by Barles and Souganidis in [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF]. Although a complete list of the contributions to numerical methods for HJI equations goes beyond the scopes of this paper, let us quote the application of Godunov/central schemes [START_REF] Abgrall | Numerical discretization of the first-order hamilton-jacobi equation on triangular meshes[END_REF][START_REF] Abgrall | Numerical discretization of boundary conditions for first order Hamilton-Jacobi equations[END_REF], antidissipative and SuperBee/UltraBee [START_REF] Bryson | Semi-discrete central-upwind schemes with reduced dissipation for Hamilton-Jacobi equations[END_REF][START_REF] Bokanowski | Anti-dissipative schemes for advection and application to Hamilton-Jacobi-Bellman equations[END_REF], MUSCL [START_REF] Lions | Convergence of MUSCL and filtered schemes for scalar conservation laws and Hamilton-Jacobi equations[END_REF], and WENO schemes [START_REF] Zhang | High-order WENO schemes for Hamilton-Jacobi equations on triangular meshes[END_REF][START_REF] Bryson | Mapped WENO and weighted power ENO reconstructions in semi-discrete central schemes for Hamilton-Jacobi equations[END_REF]. A natural way to overcome the limitations of monotone schemes is by the application of high-order approximations. For a given accuracy, these methods can achieve acceptable error levels in coarser grids, with a considerably reduced number of nodes in comparison with low-order, monotone schemes. This can be a crucial point when the dimension of the problem is high or when complex computations are required at every grid node; both situations naturally arise in the context of HJI equations stemming from optimal control and differential games.

In this paper we propose the coupling between a semi-Lagrangian (SL) time discretization with a finite volume reconstruction in space. High-order SL schemes for HJI equations have been first considered for a semi-discretization in time in [START_REF] Falcone | Discrete time high-order schemes for viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF], and for the fully discrete scheme in [START_REF] Falcone | Semi-Lagrangian schemes for Hamilton-Jacobi equations, discrete representation formulae and Godunov methods[END_REF]. A convergence analysis based on the condition ∆x = O(∆t 2 ) is carried out in [START_REF] Ferretti | Convergence of semi-Lagrangian approximations to convex Hamilton-Jacobi equations under (very) large Courant numbers[END_REF]. The adaptation of the theory to weighted ENO reconstructions is presented in [START_REF] Carlini | A weighted essentially nonoscillatory, large time-step scheme for Hamilton-Jacobi equations[END_REF], along with a number of numerical tests comparing the various high-order versions of the scheme. Other numerical tests, mostly in higher dimension and concerned with applications to front propagation and optimal control, are presented in [START_REF] Carlini | An efficient algorithm for Hamilton-Jacobi equations in high dimension[END_REF]. Let us mention that a first convergence result for a class of motonone Finite Volume schemes has been proved in [START_REF] Kossioris | Finite volume schemes for Hamilton-Jacobi equations[END_REF]. The paper is organized as follows.

In Section 2, we illustrate our ideas with a setting related to minimum time optimal control and differential games, leading to stationary HJI equation. In Section 3, we deal with a high-order approximation scheme based on a coupling between a semi-Lagrangian discretization in time and a Finite Volume spatial reconstruction. Finally, in Section 4 we present some numerical experiments assessing the performance and accuracy of the proposed scheme.

HJI equations arising in optimal control and differential games

As we mentioned in the introduction, HJI equations often arise in optimal control and differential games; whenever a feedback controller is sought, the application of the Dynamic Programming Principle (DPP) leads to HJI equations, which can be time-dependent or stationary. Among a wide class of problems, in this section we illustrate our ideas by means of minimum time optimal control and pursuit-evasion games.

Let us start by considering consider system dynamics of the form

ẏ(t) = f (y, α(t)) for t > 0 , y(0) = x , (1)
where y ∈ R n is the state, α : [0, +∞) → A is the control and f : R n × A → R n is the controlled vector field. To get a unique trajectory for every initial condition and a given control function, we will always assume that f is continuous with respect to both variables, and Lipschitz continuous with respect to the state space (uniformly in α). Moreover, we will assume that the controls are measurable functions of time so that we can apply the Carathéodory theorem for the Cauchy problem (1).

In the minimum time optimal control problem, we want to minimize the time of arrival to a given target T . The cost will be given by t(x, α) := inf{t :

y x (t; α) ∈ T } +∞ if y x (t; α) / ∈ T ∀t. (2)
By the application of the DPP one can prove that the minimum time function

T (x) := inf α∈A t(x, α)
satisfies the Bellman equation

max a∈A {-f (x, a) • DT } = 1 ,
in the domain where T is finite (the so-called reachable set). Introducing the change of variable

v(x) := 1/µ if t x (a, b) = +∞ , 1/µ(1 -e -µtx(a,b) ) elsewhere , (3)
where µ is a free positive parameter to be suitably chosen, one can characterize T as the unique viscosity solution of the following Dirichlet problem

µv(x) + max a∈A {-f (x, a) • Dv} = 1 for x ∈ R n \ T , v(x) = 0 for x ∈ ∂T . (4)
Another interesting example comes from the DPP approximation of the Hamilton-Jacobi-Isaacs equations related to pursuit-evasion games (see [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF] for more details). Player-a (the pursuer) wants to catch player-b (the evader) who is escaping, and the controlled dynamics for each player are known. To simplify the notations, we will denote by y(t) = (y P (t), y E (t)) the state of the system, where y P (t) and y E (t) are the positions at time t of the pursuer and of the evader, both belonging to R n , and by f : R 2n × A × B → R 2n the dynamics of the system (clearly, here the dynamics depend on the controls for both players). The payoff is defined as the time of capture but, in order to have a fair game, we need to restrict the strategies of the players to the so-called non-anticipating strategies (i.e. they cannot exploit the knowledge of the future strategy of the opponent). These strategies will be denoted respectively by α and β. Given the strategies α(•) and β(•) for the first and the second player, we can define the corresponding time of capture as

t x (α[β], β) = inf t > 0 : y P (t) = y E (t) .
If there is no capture for those strategies we set t x (α[β], β) = +∞. Then we can define the lower time of capture as

T (x) = inf α∈A sup β∈B t x (α[β], β),
and again T can be infinite if there is no way to catch the evader from the initial position of the system x. In order to get a fixed point problem and to deal with finite values, it is useful to scale time by the change of variable (3), which corresponds to the payoff

J x (α, β) = tx(α,β) 0 e -µt dt
The rescaled minimal time will be given by

v(x) = inf α∈A sup b∈B J x (α[β], β).
Assuming v to be continuous, the application of the DPP leads to the following characterization of the value function

v + min b∈B max a∈A {-Dv • f (x, a, b)} = 1 on R n \ T , v(x) = 0 on ∂T . (5)
Note that the equation is complemented by the natural homogeneous boundary condition on the target

T (x) = v(x) = 0. If v( • ) is continuous, then v is a viscosity solution in R n \ T of the Dirichlet problem (5).
3 Semi-Lagrangian schemes for HJI equations

In this section we introduce the main building blocks for the construction of semi-Lagrangian/finite volume schemes for HJI equations of the form (4)-( 5). The general procedure is decomposed into a time discretization step, and a space discretization procedure. In the time discretization step, the system dynamics (1) are approximated by a suitable integration rule, and the DPP is applied on its discrete-time version. In the space discretization procedure, the resulting HJI equation is then approximated over a finite set of elements. The same procedure holds for all the problems presented in the previous section, however, for the sake of simplicity, this section is illustrated by means of the minimum time problem and its associated HJB equation ( 4) (which is also a particular case of ( 5) for a single player setting).

Time discretization

For the implementation of a time disretization procedure, we follow the ideas presented in [START_REF] Falcone | Discrete time high-order schemes for viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]. The first step towards the construction of a high-order scheme for the equations ( 4)-( 5) is to consider discrete time approximation of the system dynamics (1) of the form

y n+1 = y n + h Φ(y n , A n , h) , for n > 0 , (6)
y 0 = x ,
where h > 0 corresponds to a time discretization parameter, Φ = Φ(y n , A n , h) is the Henrici function of a one-step approximation of the dynamical system, and A n stands for a multidimensional control defined accordingly to the order q of the numerical integrator,

A n = (a 0 n , a 1 n , . . . , a q n ) , A n ∈ A q+1 , q ≥ 0.
Particular cases of the aforementioned setting are

(i) Explicit Euler's method: Φ(y n , A n , h) = f (y n , a n ), with A n = a 0 n . (ii) Midpoint rule: Φ(y n , A n , h) = f (y n + hf (y n , a 0 n )/2, a 1 n ), with A n = (a 0 n , a 1 n ). (iii) Fourth-order Runge-Kutta scheme: Φ(y n , A n , h) = 1 6 (K 0 + 2K 1 + 2K 2 + K 3 ) , A n = (a 0 n , a 1 n , a 2 n , a 3 n ) , K 0 = f (y n , a 0 n ) , K 1 = f (y n + h K 0 2 , a 1 n ) , K 2 = f (y n + h K 1 2 , a 2 n ) , K 3 = f (y n + hK 2 , a 3 n ) .
The application of the DPP for the discrete-time dynamics leads to an approximation of eq. ( 4) of the form

v h (x) = min An∈A q {βv h (x + hΦ(x, A n , h))} + 1 -β for x ∈ R n \ T , v(x) = 0 for x ∈ ∂T , (7)
where β = e -h appears as a consequence of the application of the Kruzkhov transform (3) to the discrete version of the minimum time solution (in this case with µ = 1).

Space discretization

Note that, despite having introduced an approximation in time for eq. ( 4), the resulting semi-discrete version ( 7) is still continuously defined over the state space. In order to implement a fully-discrete computational scheme, it is necessary to realize this expression in a bounded domain with a finite set of elements. Classical schemes for HJI equations of this form are based upon finite difference discretizations, where the domain Ω ⊂ R n over which the solution is sought, is discretized into a set of grid points, and the approximation is understood in a pointwise sense. A natural problem in this setting arises from the fact that the r.h.s. of eq. ( 7) requires the evaluation of v h at the arrival points x+hΦ(x, A n , h), which are not necessarily part of the grid. In the low-order version of the SL scheme, this evaluation is performed via piecewise linear interpolation from the grid values, whereas in this work we focus on a high-order definition of such an operation. We follow an approach based on a Finite Volume approximation of the problem. For a given mesh parameter k, and a set of central nodes

{x i } N i=1 , the domain is discretized into a set of cells Ω i = [x i -k/2, x i + k/2].
Instead of considering pointwise nodal values of v h , the solution will be represented by a set of cell-averaged values V := {v i } N i=1 defined as

v i := 1 k xi+k/2 xi-k/2 v h (x) dx , i = 1, . . . , N.
Its is straightforward to see that the exact expression for the averaged values of the solution of ( 7) is given by

v i = T k,i (v h ) for i = 1, . . . , N , T k,i (v h ) : = 1 k xi+k/2 xi-k/2 min An∈A q {βv h (x + hΦ(x, A n , h))} + 1 -β dx , (8)
v(x) = 0 for x ∈ ∂T .

A first approximation is introduced when the integral in ( 8) is replaced by a suitable Gaussian quadrature rule

T k,i (v h ) ≈ 1 k i w i min An∈A q {βv h (x i + hΦ(x i , A n , h))} + 1 -β , (9)
where x i and w i are Gauss points and weights inside the i-th cell, respectively. This expression requires the evaluation of the exact v h at a set of arrival points, which is not available. Analogously to the grid-based schemes, we approximate this evaluation with an interpolation operator I = I[V ] defined upon the set of cell averages, i.e.

v h (x i + hΦ(x i , A n , h)) ≈ I[V ](x i + hΦ(x i , A n , h)).
where I : R N →→ S k corresponds to a WENO (weighted essentially nonoscillatory) interpolation routine performed over the averaged dataset V . The WENO reconstruction procedure and related numerical schemes date back to the work of [START_REF] Liu | Weighted Essentially Non-oscillatory Schemes[END_REF], in the context of numerical methods for conservation laws, as a way of circumventing Godunov's barrier theorem by considering nonlinear (on the data) reconstruction procedures for the implementation of high-order accurate schemes. As it has been shown in [START_REF] Ferretti | Convergence of semi-Lagrangian approximations to convex Hamilton-Jacobi equations under (very) large Courant numbers[END_REF], the use of a WENO interpolation procedure can be considered as a building block in high-order, semi-Lagrangian schemes for time-dependent HJB equations, whereas here we introduce an application to static HJI equations. We now briefly describe the main ideas for a 1D WENO reconstruction.

From a set of cell values V and a polynomial degree r, the WENO reconstruction procedure yields a set of polynomials P = {p i (x)} N i=1 of degree r, holding standard interpolation properties

v i = 1 k Ωi p i (x) dx , v(x) = p i (x) + o(∆x r ) , ∀x ∈ Ω i , (10)
and an essentially non-oscillatory condition [START_REF] Harten | Uniformly high order accurate essentially non-oscillatory schemes[END_REF]; in general, such an interpolant is built by considering a set of stencils per cell, and weighting them according to a smoothness indicator. Several variations of this procedure can be found in the literature; for illustration purposes, we restrict ourselves to the reconstruction procedure presented in [START_REF] Balsara | Efficient, high accuracy ADER-WENO schemes for hydrodynamics and divergence-free magnetohydrodynamics[END_REF], on its 1D version, and reconstruction degree 2. In this case, given a set of averaged values V , the reconstruction procedure seeks, for every cell, a local quadratic expansion upon a linear combination of Legendre polynomials rescaled in local coordinates ξ = [-1/2, 1/2], expressed in the form

p(ξ) = v 0 + v ξ p 1 (ξ) + v ξξ p 2 (ξ), with p 1 (ξ) = ξ , p 2 (ξ) = ξ 2 - 1 12 .
We assign the subscript "0" to the cell where we compute the coefficients, other values indicating location and direction with respect to v 0 (note that the notation is coherent with the fact that the first coefficient in the expansion v 0 , holds v 0 = v i , i.e., the centered value). Next, for this particular problem we define three stencils

S 1 = {v -2 , v -1 , v 0 } , S 2 = {v -1 , v 0 , v 1 } , S 3 = {v 0 , v 1 , v 2 } ,
and in every stencil we compute a polynomial of the form

v (i) (ξ) = v (i) 0 + v (i) ξ p 1 (ξ) + v (i) ξξ p 2 (ξ) i = 1, 2, 3.
Imposing the conservation condition [START_REF] Bokanowski | Anti-dissipative schemes for advection and application to Hamilton-Jacobi-Bellman equations[END_REF], the coefficients are given by

S 1 : v (1) ξ = -2v -1 + v -2 /2 + 3v 0 /2, v (1)
ξξ = (v -2 -2v -1 + v 0 )/2 , S 2 : v (2) ξ = (v 1 -v -1 )/2, v (2) ξξ = (v -1 -2v 0 + v 1 )/2 , S 3 : v (3) ξ = -3v 0 /2 + 2v 1 -v 2 /2, v (3) ξξ = (v 0 -2v -1 + v 2 )/2.
For every polynomial we calculate a smoothness indicator defined as

IS (i) = r l=1 Ω0 k 2l-1 ∂ l p (i) ∂x l 2 dx ,
where r is the polynomial reconstruction degree (in our case r = 2), and which in our case yields to

IS (i) = p (i) ξ 2 + 13 3 p (i) ξξ 2 .
This leads to the following WENO weights:

ω (i) = α (i) 3 i=1 α (i) , α (i) = λ (i) ( + IS (i) ) r ,
where is a parameter introduced in order to avoid division by zero; usually = 10 -12 . The scheme is rather insensitive to the parameter r, which we set r = 5. The parameter λ is usually computed in an optimal way to increase the accuracy of the reconstruction at certain points; we opt for a centered approach instead, thus λ (1) = λ (3) = 1, while λ (2) = 100. Finally, the expression for the 1D reconstructed polynomial at the i -th cell is given by

p i (ξ) = ω (1) p (1) (ξ) + ω (2) p (2) (ξ) + ω (3) p (3) (ξ).
Having defined all the buildings blocks for a fully-discrete, high-order approximation of eq. ( 4), we need to solve the following nonlinear system

v i = [T k,i (V )] i for i = 1, . . . , N , (11)
[T k,i (V )] i : = 1 k i w i min An∈A q {βI[V ](x i + hΦ(x i , A n , h))} + 1 -β , v(x) = 0 for x ∈ T , v(x) = 1 for x ∈ ∂Ω c \T .
Note that we added an additional boundary condition related to the external part of the computational domain which is not the target, computationally equivalent to setting a high value which is neglected in the minimization procedure for the interior elements. The computational domain must be set accordingly to this condition, in order to generate a consistent result. With respect to the solution of the nonlinear system [START_REF] Bryson | Semi-discrete central-upwind schemes with reduced dissipation for Hamilton-Jacobi equations[END_REF], the approach which we follow is motivated by the standard approach undertaken in the low-order setting, which is to solve the system by some variation of a fixed point iteration

V n+1 = T (V n ) , (12)
which, in the low-order monotone scheme, is well-justified since T is a contraction mapping. A key point is the fact that the corresponding linear interpolation operator is monotone, which is lost in the high-order scheme. However, it is still possible to develop a convergence theory for interpolation operators that are not monotone but have additional properties such as the WENO operator.

In [START_REF] Ferretti | Convergence of semi-Lagrangian approximations to convex Hamilton-Jacobi equations under (very) large Courant numbers[END_REF], a convergence framework has been developed for time-dependent HJB equations, and recently in [START_REF] Bokanowski | Value iteration convergence of -monotone schemes for stationary Hamilton-Jacobi equations[END_REF], convergence results have been obtained for the stationary case. One of the advantages of this setting is the vast amount of available literature dealing with acceleration techniques for HJI iterative solvers (we refer the reader to [START_REF] Alla | An efficient policy iteration algorithm for dynamic programming equations[END_REF] and references therein for a recent update on such methods).

Numerical examples

We now present two numerical examples assessing the performance of the proposed scheme. We recall that, although we will present examples dealing with minimum time optimal control and pursuit-evasion games, the presented ideas can be applied in a straightforward manner to infinite/finite horizon optimal control, reachability analysis and differential games.

A two-dimensional minimum time problem

We begin by considering a two-dimensional minimum time problem. In this first example, system dynamics are given by f (x, y, (a 1 , a 2 )) = a 1 a 2 , the domain is Ω =] -1, 1[ 2 , the target is T = ∂Ω, h = 0.8k and A = {(1, 0), (0, 1), (-1, 0), (0, -1)} is the set of 4 directions pointing to the facets of the unit square. The exact solution for this problem is the distance function to the unit square. As the characteristic curves for this problem correspond to straight lines moving towards the boundary, integration in time can be achieved exactly with a solver of any order. We consider then a Euler discretization in time, and a two-dimensional WENO reconstruction of order 2 in space. The multidimensional WENO reconstruction is based on a product of unidimensional reconstructions along every direction (we refer the reader to [START_REF] Balsara | Efficient, high accuracy ADER-WENO schemes for hydrodynamics and divergence-free magnetohydrodynamics[END_REF] for the specific version used in this test). Convergence rates and errors are shown in Table 1.

Note that the second order of the space interpolation is achieved for the • 1 norm, while a lower order is observed for the • ∞ norm. This is expected from the fact that the solution is not differentiable across the kinks of the solution. Note that if error computation is performed over a restricted zone excluding non-differentiable points, as in Table 2, higher order of accuracy and convergence are achieved for both norms, as it is generally expected for WENO-based schemes. However, this might not be the case for any high-order scheme. As an example, in Figure 1, it can be seen that if a generic quadratic interpolation is used to build a similar scheme, spurious oscillations arise in non-differentiable areas, degenerating the high-order accuracy of the scheme. This latter justifies the use of WENO reconstruction operators in space, as they are accordingly designed in order to detect and penalize highly oscillatory stencils. A reduced-coordinate pursuit-evasion game

In a second example, we consider a 1D pursuit-evasion game with dynamics given by ẋP

= v P a , ẋE = v E b ,
where v P and v E denote the velocity of the pursuer and the evader respectively; a ∈ [0, 1] and b ∈ [-1, 1] are control variables. By defining the reduced coordinate x = x E -x p , the game is written as If we consider the target set T = B(0, R), the exact solution is given by

ẋ = v E b -v P a.
v(x) =      1 -exp(-|x + R|) if x < R 0 if x ∈ (-R, R) 1 if x > R .
We implement our SL/FV scheme with a fourth-order RK scheme in time and a WENO reconstruction in space of degree 2; results are shown in Figure 2. A natural advantage of high-order methods is the level of accuracy that can be reached with a reduced number of elements, which is particularly relevant when fixed point iterations of the form [START_REF] Bryson | Mapped WENO and weighted power ENO reconstructions in semi-discrete central schemes for Hamilton-Jacobi equations[END_REF] involving min or minmax operators are considered. However, at it has been previously discussed, for high-order schemes the fixed point operator is not a contraction anymore, and convergence has to be understood in a different sense. In [START_REF] Bokanowski | Value iteration convergence of -monotone schemes for stationary Hamilton-Jacobi equations[END_REF], the -monotonicity concept has been introduced in order to characterize the convergence behavior of such high-order schemes. Figure 2 illustrates this situation, as for different values of h and k, convergence of the fixed point iteration is achieved in an oscillatory way, whereas the oscillation behavior decreases when h = h(k) is reduced.

Concluding remarks

We have introduced a semi-Lagrangian/finite volume scheme for the approximation of HJI equations. The main building blocks are a high-order approximation of the system dynamics, combined with high order of accuracy in space via a WENO interpolation operator. The resulting fully-discrete scheme is then solved by means of a fixed point iteration. High-order of accuracy is observed in smooth regions, and the convergence of the fixed point iteration is achieved as long as in the spatial-resolution building block, non-oscillatory interpolation operators are considered. Further developments in the directions of this paper will include the implementation of high-dimensional interpolation routines, as well as the construction of adaptive schemes with an ad-hoc refinement criterion.

Fig. 1 .

 1 Fig.1. 2D minimum time problem: value function for different schemes. Left: highorder scheme with WENO reconstruction in space. Right: high-order scheme using quadratic interpolation.

Fig. 2 .

 2 Fig.2. SL/FV scheme for a 1D differential game. Left: exact and approximated solution for 100 elements. Right: oscillatory, but convergent behavior is achieved for the fixed point iteration[START_REF] Bryson | Mapped WENO and weighted power ENO reconstructions in semi-discrete central schemes for Hamilton-Jacobi equations[END_REF], for different values of h.

Table 1 .

 1 Errors for the 2D minimum time problem with a second-order WENO reconstruction.

	k	• ∞-error • ∞-order • 1-error • 1-order #iterations
	0.08 0.0023	-		3.667e-4	-	21
	0.04 0.0011	1.0641	1.068e-4	1.7791	30
	0.02 5.711e-4	0.946		2.887e-5	1.8878	54
	0.01 2.966e-4	0.945		7.51e-6	1.9427	104
		0.08 7.763e-9	-		1.051e-9	-
		0.04 5.362e-10	3.8558	1.438e-10	2.8696
		0.02 3.523e-11	3.9279	1.0758e-11 3.7406
		0.01 2.257e-12	3.9643	7.1681e-13 3.9076

k • ∞-error • ∞-order • 1-error • 1-order

Table 2 .

 2 Error computation (as in Table1) over a restricted zone excluding nondifferentiable kinks.

This research was supported by the following grants: AFOSR Grant no. FA9550-10-1-0029, ITN-Marie Curie Grant no. 264735-SADCO, and the FWF-START project Sparse Approximation and Optimization in High Dimensions.


Elena K Kostousova

N N Krasovskii

On target control synthesis under set-membership uncertainties using polyhedral techniques

Keywords: Differential systems, uncertain systems, control synthesis, polyhedral estimates, parallelotopes, interval analysis

Problems of feedback terminal target control for linear and bilinear uncertain systems are considered. We continue the development of control synthesis using polyhedral (parallelotope-valued) solvability tubes. The paper deals with two types of problems, where controls appear either additively or in the system matrix. For both problems, the cases without uncertainties, with additive parallelotope-bounded uncertainties, and also with interval uncertainties in coefficients of the system (a bilinear uncertainty) are considered. Ordinary differential equations, which describe the mentioned polyhedral solvability tubes, are presented for each of these cases. New control strategies, which can be calculated by explicit formulas on the base of the mentioned tubes, are proposed. Results of computer simulations are presented.

Introduction

Problems of feedback terminal target control for linear and bilinear differential uncertain systems are considered. There are known approaches for solving problems like these, in particular, based on constructing solvability tubes and the extremal aiming strategies of N.N. Krasovskii [START_REF] Krasovskii | Positional Differential Games[END_REF][START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF]. The problem statement for linear systems, approaches for solving, and the tight interconnections between solvability tubes, the Pontriagin alternated integral, Hamilton-Jacobi-Bellman equations, and funnel equations can be found, for example, in [START_REF] Kurzhanskii | On the Problem of the Synthesis of Controls: the Pontryagin Alternative Integral and the Hamilton-Jacobi Equation[END_REF][START_REF] Kurzhanski | On the Problem of Synthesizing Control Strategies. Evolution Equations and Set-Valued Integration[END_REF][START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF].

Since practical construction of the mentioned tubes can be cumbersome, different numerical methods are devised, in particular, methods for approximating the set-valued integrals and for numerical solving the mentioned equations, including methods based on approximations of sets by arbitrary polytopes with a large number of vertices [START_REF] Baier | Computing Aumanns Integral[END_REF][START_REF] Chernousko | State Estimation for Dynamic Systems[END_REF][START_REF] Taras'yev | Approximation Schemas and Finite-Difference Operators for Constructing Generalized Solutions of Hamilton-Jacobi Equations[END_REF][START_REF] Veliov | Second Order Discrete Approximations to Strongly Convex Differential Inclusions[END_REF] (here and below we mention, as examples, only some references from numerous publications; see also references therein). Such methods are devised to obtain approximations as accurate as possible. But they can require much calculations, especially for large dimensional systems. Other techniques are based on estimates of sets by domains of some fixed shape such as ellipsoids and parallelepipeds, including boxes aligned with coordinate axes as in interval analysis [4, 5, 7-12, 14, 15, 17, 18, 20]. The main advantage of such techniques is that they enable to obtain approximate/particular solutions using relatively simple tools (up to explicit formulas). More accurate approximations may be obtained by using the whole families (varieties) of such simple estimates (as was proposed by A.B. Kurzhanski) [START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF][START_REF] Kurzhanskii | On the Problem of the Synthesis of Controls: the Pontryagin Alternative Integral and the Hamilton-Jacobi Equation[END_REF][START_REF] Kurzhanski | On Ellipsoidal Techniques for Reachability Analysis. Part I: External Approximations. Part II: Internal Approximations. Box-valued Constraints[END_REF][START_REF] Kostousova | Guaranteed Estimates of Accuracy of Computations in Problems of Control and Estimation[END_REF][START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF].

For linear differential systems, the constructive computation schemes for solving the feedback target control problems by means of ellipsoidal techniques were proposed [START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF][START_REF] Kurzhanskii | On the Problem of the Synthesis of Controls: the Pontryagin Alternative Integral and the Hamilton-Jacobi Equation[END_REF] and then expanded to a polyhedral technique [START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF]. Here we continue the development of the polyhedral control synthesis using polyhedral (parallelotope-valued) solvability tubes. The paper deals with two types of problems, where the controls appear either additively or in the system matrix. For both problems, the cases without uncertainties, with additive uncertainties, and also with interval uncertainties in coefficients of the system (the bilinear uncertainty) are considered. Ordinary differential equations (ODE) for the mentioned polyhedral solvability tubes are presented. New control strategies, which can be calculated by explicit formulas on the base of the mentioned tubes, are proposed. In opposite to [START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF][START_REF] Kurzhanskii | On the Problem of the Synthesis of Controls: the Pontryagin Alternative Integral and the Hamilton-Jacobi Equation[END_REF][START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF], they are concretized by explicit formulas when the state belongs to a tube. Also the polyhedral control synthesis for discrete-time systems is considered. The results of computer simulations are presented.

Note that there are also some works devoted to other approaches for solving different control problems under uncertainty and works concerning systems with bilinear uncertainties (see, for example, [START_REF] Anan'evskii | Synthesis of a Bounded Control for Linear Dynamical Systems Using the General Lyapunov Function[END_REF][START_REF] Chernousko | State Estimation for Dynamic Systems[END_REF][START_REF] Filippova | On the Estimation of Trajectory Tubes of Differential Inclusions[END_REF][START_REF] Mazurenko | A Differential Equation for the Gauge Function of the Star-Shaped Attainability Set of a Differential Inclusion[END_REF][START_REF] Polyak | Robust Stability and Control[END_REF]).

The following notation is used below: R n is the n-dimensional vector space; is the transposition symbol; x 2 = (x x) 1/2 , x ∞ = max 1≤i≤n |x i | are vector norms for x=(x 1 , x 2 , . . . , x n ) ∈ R n ; e i =(0, . . . , 0, 1, 0, . . . , 0) is the unit vector oriented along the axis 0x i (the unit stands at position i); e = (1, 1, . . . , 1) ; R n×m is the space of real n×m-matrices A = {a j i } = {a j } (with columns a j ); I is the identity matrix; 0 is the zero matrix (vector); Abs A = {|a j i |} for A = {a j i }; diag π, diag {π i } are the diagonal matrix A with a i i = π i (π i are the components of the vector π); det A is the determinant of A; tr A= n i=1 a i i is the trace of A; A = max 1≤i≤n m j=1 |a j i | for A∈R n×m ; int X is the set of interior points of the set X ⊂ R n ; the notation of the type k=1, . . . , N is used instead of k=1, 2, . . . , N .

Problems formulation

Consider the controlled system with a given terminal set M (x∈R n is the state):

ẋ = (A(t) + U (t) + V (t)) x + u(t) + v(t), t ∈ T = [0, θ]. (1)
Here A(t) ∈ R n×n is a given matrix function; Lebesgue measurable functions U (t) ∈ R n×n and u(t) ∈ R n serve as controls and satisfy either (2) or (3):

U (t) ≡ 0, u(t) ∈ R(t), a.e. t ∈ T, (2)
U (t) ∈ U(t) = {U ∈ R n×n | Abs (U -Ũ (t))≤ Û (t)}, u(t) ≡ 0, a.e. t ∈ T ; (3)
V (t) ∈ R n×n and v(t) ∈ R n stand for unknown disturbances and satisfy

V (t) ∈ V(t) = {V ∈ R n×n | Abs (V -Ṽ (t))≤ V (t)}, v(t) ∈ Q(t), a.e. t ∈ T. (4)
Matrix and vector inequalities (≤, <, ≥, >) here and below are understood componentwise. We presume the sets R(t), Q(t), and M to be parallelotopes and a parallelepiped respectively:

R(t) = P[r(t), R(t)], R(t) ∈ R n×n1 , Q(t) = P[q(t), Q(t)], Q(t) ∈ R n×n2 , M = P(p f , P f , π f ) = P[p f , Pf ],
Pf ∈ R n×n , det Pf = 0;

(5) r(t), R(t), q(t), Q(t), as well as A(t), Ũ (t), Û (t) ≥ 0, Ṽ (t), V (t) ≥ 0, are known continuous vector and matrix functions; the parallelepiped M is nondegenerate.

By a parallelepiped P(p, P , π)⊂R n we mean a set such that P =P(p,

P , π) = {x ∈ R n | x=p+ n i=1 p i π i ξ i , ξ ∞ ≤1}, where p∈R n ; P ={p i }∈R n×n is such that det P =0, p i 2 =1 1 ; π ∈ R n , π ≥ 0.
It may be said that p determines the center of the parallelepiped, P is the orientation matrix, p i are the "directions" and π i are the values of its "semi-axes". We call a parallelepiped nondegenerate if π>0.

By a parallelotope P[p, P ] ⊂ R n we mean a set

P = P[p, P ] = {x ∈ R n | x = p+ P ζ, ζ ∞ ≤ 1}
, where p ∈ R n and the matrix P = {p i } ∈ R n×m , m ≤ n, may be singular. We call a parallelotope P nondegenerate if m = n and det P = 0.

Each parallelepiped P(p, P , π) is a parallelotope P[p, P ] with P = P diag π; each nondegenerate parallelotope is a parallelepiped with

P = P diag { pi -1 2 }, π i = pi
2 or, in a different way, with P = P , π = e, where e = (1, 1, . . . , 1) . We can consider the above system for the following cases: (I) without uncertainty when v and V ≡ 0 are given functions, i.e., Q ≡ 0, Ṽ ≡ V ≡ 0; (II) under uncertainty including the following three subcases: (II,i) only additive uncertainty (V ≡ 0); (II,ii) only matrix uncertainty ( Q ≡ 0); (II,iii) both ones.

In [START_REF] Kurzhanskii | On the Problem of the Synthesis of Controls: the Pontryagin Alternative Integral and the Hamilton-Jacobi Equation[END_REF][START_REF] Kurzhanski | On the Problem of Synthesizing Control Strategies. Evolution Equations and Set-Valued Integration[END_REF][START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF], for cases (I) and (II,i) with controls (2), the following problem of terminal target control synthesis under uncertainty was investigated.

Problem 1. For the system (1),( 2),(4), case (I) or (II,i), specify a solvability set W(τ, θ, M) = W(τ ) and a set-valued feedback control strategy 1 The normality condition p i 2=1 may be omitted to simplify formulas. 2 Here the class U c R of feasible control strategies is taken to consist of all convex compact-valued multifunctions u(t, x) that are measurable in t, upper semicontinuous in x, being restricted by u(t, x) ⊆ R(t), t ∈ T . The condition u(•, •) ∈ U c R ensures that the corresponding differential inclusion does have a solution.

2 u = u(t, x), u(•, •) ∈ U c R , such that all solutions to the differential inclusion ẋ ∈ A(t)x + u(t, x) + Q(t), t ∈ T , that start from any given position {τ, x τ }, x τ = x(τ ) ∈ W(τ, θ, M), τ ∈ [0, θ), would reach the terminal set M at time θ: x(θ) ∈ M.
The multivalued function W(t), t ∈ T , is known as a solvability tube W(•). The ellipsoidal synthesis was elaborated in [START_REF] Kurzhanskii | On the Problem of the Synthesis of Controls: the Pontryagin Alternative Integral and the Hamilton-Jacobi Equation[END_REF][START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF] for solving Problem 1. In [START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF], the families of external P + (•) and internal P -(•) parallelotope-valued (shorter, polyhedral ) estimates for W(•) were introduced. The extremal aiming strategies of N.N. Krasovskii were used there. They were constructed in an analytical form on the base of a solution of some specific mathematical programming problem. Now let us consider two following problems, which concern all above cases of uncertainties. Unlike Problem 1, they involve single-valued control strategies. This is possible because our strategies will be continuous and even linear with respect to x. Moreover, they will be constructed in an explicit form.

Problem 2. For the system (1), ( 2), ( 4), ( 5), find a polyhedral tube P -(t) = P[p -(t), P -(t)], t ∈ T , with P -(θ) = M, and find a corresponding feedback control strategy u = u(t, x) such that u(t, x) ∈ R(t) for x ∈ P -(t), t ∈ T , and each solution x(•) to the differential equation ẋ = (A(t) + V (t))x + u(t, x) + v(t), t ∈ T , with x(0) = x 0 ∈ int P -(0) would be defined on T and would satisfy x(t) ∈ P -(t), t ∈ T , whatever are v(•) and V (•) subjected to (4). Moreover, introduce a whole family of such tubes P -(•). Problem 3. For the system (1), ( 3), ( 4), ( 5), find a polyhedral tube P -(t) = P[p -(t), P -(t)], t ∈ T , with P -(θ) = M, and find a corresponding feedback control strategy U = U (t, x) such that U (t, x) ∈ U(t) for x ∈ P -(t), t ∈ T , and each solution x(•) to the differential equation

ẋ = (A(t) + U (t, x) + V (t))x + v(t), t ∈ T, (6)
with x(0)=x 0 ∈int P -(0) would be defined on T and would satisfy x(t)∈P -(t), t∈T , whatever are V (•), v(•) subjected to [START_REF] Chernousko | State Estimation for Dynamic Systems[END_REF]. Introduce a family of such tubes.

Solutions to Problem 2

First, let us consider the following ODE system for P -(t) = P[p -(t), P -(t)]:

dp - dt = (A(t) + Ṽ (t))p -+ r(t) + q(t), p -(θ) = p f ; (7)
d P - dt = (A(t) + Ṽ (t)) P -+ P -diag β(t, P -) + R(t) Γ (t) + P -diag γ(t, P -), β(t, P -) = max{Abs (( P -) -1 ) V (t) Abs (p -(t) + P -ξ) | ξ ∈ E(C)}, γ(t, P -) = Abs (( P -) -1 Q(t)) e, P -(θ) = Pf . (8)
Here (and below) the operation of maximum is understood componentwise, E(C) denotes the set of all vertices of C = P(0, I, e) (i.e., points ξ ∈ R n with ξ j ∈ {-1, 1}); Γ (t) ∈ R n1×n is an arbitrary Lebesgue measurable matrix function satisfying

Γ (t) ∈ G, a.e. t ∈ T , where G = {Γ = {γ j i } ∈ R n1×n | Γ ≤ 1}, Γ = max 1≤i≤n1 n j=1 |γ j i |.
Let G be the set of all such functions Γ (•). Let us consider the following control strategy, being connected with P -(•) from ( 7), [START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF]:

u(t, x) = r(t) + R(t)Γ (t) P -(t) -1 (x -p -(t)). ( 9
)
Theorem 1. We consider the system (1),( 2),( 4), [START_REF] Filippova | Trajectory Tubes of Nonlinear Differential Inclusions and State Estimation Problems[END_REF], where det Pf = 0. Let Γ (•) ∈ G. Then the system ( 7), ( 8) has a unique solution (p -(•), P -(•)) at least on some subinterval T 1 =[τ 1 , θ]⊆T , where 0≤τ 1 <θ. If T 1 =T and we have det P -(t) =0, t∈T , then the tube P -(•) and the control strategy [START_REF] Kostousova | State Estimation for Control Systems with a Multiplicative Uncertainty through Polyhedral Techniques[END_REF] give a particular solution to Problem 2; in cases (I), (II,i), all solutions x(•) with x(0) ∈ P -(0) (not only with x(0) ∈ int P -(0)) generated by ( 9) satisfy x(t) ∈ P -(t), t ∈ T .

The scheme of the proof is similar to the proof of Theorem 2 (see below). Theorem 1 describes the whole family of tubes P -(•), where Γ (•) serves as a parameter. Thus the set W 0 = {int P -(0)| Γ (•)∈G such that det P -(t) =0, t∈T } (or, in cases (I), (II,i), the analogous set W 0 = P -(0)) provides the set of initial positions which can be steered to the terminal set M during the time θ by solving Problem 2. But, generally speaking, it is not true that det P -(0) = 0 or even P -(0) = ∅ for each Γ (•)∈G. For cases (I), (II,i), the above family of the tubes P -(•) coincides with the family of internal estimates for W(•) introduced in [START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF]. It follows from [START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF][START_REF] Kostousova | On Tight Polyhedral Estimates for Reachable Sets of Linear Differential Systems[END_REF] that for the case (I) we have T 1 = T for each Γ (•) ∈ G and W(0)= {P -(0)| Γ (•)∈G}. But we can not conclude from here that W 0 = W(0). The attractive property of the control strategies ( 9) is their explicit form.

Remark 1. One of the heuristic ways to construct the parameter Γ (•) is to apply arguments of a "local" volume optimization similarly to [START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF] (see also Remark 2 below). Namely, assuming det Pf > 0 and introducing a grid T N of times τ k = kh N , k=0, . . . , N , h N = θN -1 , we can construct the piecewise constant function

Γ (t) ≡ Γ (τ k ) ∈ Argmin Γ ∈G tr ( P -(τ k ) -1 R(τ k )Γ ), t ∈ (τ k-1 , τ k ], k=N, . . . , 1.
For case (I), we can use, similarly to [START_REF] Kostousova | On Tight Polyhedral Estimates for Reachable Sets of Linear Differential Systems[END_REF], minimization over Γ that satisfy Γ ∈G and some constraints introduced to produce tight estimates P -(t) for W(t). Solutions of both optimization problems are known in the explicit form [START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF][START_REF] Kostousova | On Tight Polyhedral Estimates for Reachable Sets of Linear Differential Systems[END_REF].

Solutions to Problem 3

Let us consider the following ODE system for P -(t) = P[p -(t), P -(t)]:

dp - dt = (A(t) + Ũ (t) + Ṽ (t))p -+ q(t), p -(θ) = p f ; (10)
d Pdt = (A(t)+ Ũ (t)+ Ṽ (t)) P --diag α(t, P -) P -+ Pdiag (β(t, P -)+γ(t, P -)), α i (t, P -)=α i (t, P -; J(t))=û ji i (t) η ji (t, P -)(e i (Abs P -) e) -1 , i=1, . . . , n, η(t, P -) = max{0, Abs p -(t)-(Abs P -)e},

β(t, P -) = max{Abs (( P -) -1 ) V (t) Abs (p -(t) + P -ξ) | ξ ∈ E(C)}, γ(t, P -) = Abs (( P -) -1 Q(t)) e, P -(θ) = Pf , (11)
where ûj i stand for elements of Û . Here J = {j 1 , . . . , j n } is an arbitrary permutation of numbers {1, . . . , n} or even a measurable vector function with values J(t) being arbitrary permutations. Let J be the set of all such functions J(•). Let us consider the control strategy connected with P -(•) from ( 10), ( 11):

e i U (t, x)= e i Ũ (t)-α i (t, P -(t))(x i -p - i (t))(x ji ) -1 e ji if x ji = 0, e i Ũ (t) if x ji = 0, i=1, . . . , n. (12)
Theorem 2. We consider the system (1), ( 3)-( 5), where det Pf = 0. Let J(•) ∈ J. Then the system (10), ( 11) has a unique solution (p -(•), P -(•)) at least on some subinterval T 1 = [τ 1 , θ] ⊆ T , where 0 ≤ τ 1 < θ. If T 1 = T and we have det P -(t) = 0, t ∈ T , then the tube P -(•) and the control strategy [START_REF] Kostousova | Guaranteed Estimates of Accuracy of Computations in Problems of Control and Estimation[END_REF] give a particular solution to Problem 3; in cases (I), (II,i), all solutions x(•) to ( 6) with x(0) ∈ P -(0) (not only with x(0) ∈ int P -(0)) satisfy x(t) ∈ P -(t), t ∈ T .

Proof. Here we give a sketch. First, it can be checked that the strategy [START_REF] Kostousova | Guaranteed Estimates of Accuracy of Computations in Problems of Control and Estimation[END_REF] acts for x ∈ P -(t) according to the rule

U (t, x) x = Ũ (t) x -diag α(t, P -(t)) • (x-p -(t)).
Existence and uniqueness of the solution follow from the known results similarly to [START_REF] Kostousova | Control Synthesis via Parallelotopes: Optimization and Parallel Computations[END_REF][START_REF] Kostousova | On Polyhedral Estimates for Reachable Sets of Differential Systems with Bilinear Uncertainty[END_REF]. Let x 0 ∈ int P -(0) (x 0 ∈ P -(0) for cases (I) and (II,i)). Let x(•) be the solution of ( 6) that corresponds to x(0) = x 0 (i.e., x(0

) = p -(0) + P -(0)ζ 0 , where ζ 0 ∞ < 1 (respectively, ζ 0 ∞ ≤ 1)
), to the control U (t, x) from ( 12), and to arbitrary admissible functions v(•) 6), ( 10), [START_REF] Kostousova | Guaranteed Estimates of Accuracy of Computations in Problems of Control and Estimation[END_REF], it is not difficult to see that

(such that v(t) = q(t) + Q(t)χ(t), χ(t) ∞ ≤ 1)
(x -p -) = (A + Ũ + Ṽ )(x -p -) - (diag α)(x -p -) + (V -Ṽ )x + v -q, which follows from (
ζ = -(diag β + diag γ)ζ + ( P -) -1 ((V -Ṽ )x + v -q). Let us denote b(t) = β(t, P -(t)) + γ(t, P -(t)), c(t, ζ) = P -(t) -1 ((V (t) -Ṽ (t)) • (p -(t) + P -(t)ζ) + Q(t)χ(t)
). Then, using ( 4), ( 5), we have

ζi = -b i (t)ζ i + c i (t, ζ), i = 1, . . . , n, ζ(0) = ζ 0 ; b(t) ≥ 0, Abs c(t, ζ) ≤ b(t) for ζ ∈ C = P(0, I, e). ( 13
)
It is not difficult to check that if ζ(•) satisfies ( 13) and

ζ 0 ∈ int C, then ζ(t) ∈ int C, t ∈ T ; if ζ 0 ∈ C
(x i -p - i )(x ji ) -1 | ≤ ûji i , i=1, . .

. , n (this can be obtained by simple estimates).

Theorem 2 describes the family of tubes P -(•), where J(•) serves as a parameter. Thus the set W 0 = {int P -(0)| J(•)∈J such that det P -(t) =0, t∈T } provides the set of x 0 that can be steered to M during the time θ by solving Problem 3. However it is not true that det P -(0) =0 or P -(0) =∅ for each J(•)∈J.

Remark 2. One of the ways of constructing J(•) is to apply arguments of a "local" volume optimization similarly to [START_REF] Kostousova | State Estimation for Control Systems with a Multiplicative Uncertainty through Polyhedral Techniques[END_REF][START_REF] Kostousova | On Polyhedral Estimates for Reachable Sets of Differential Systems with Bilinear Uncertainty[END_REF]. Namely, assume, without loss of generality, that det Pf >0. Fix a natural number N and introduce a grid T N of times τ k =kh N , k=0, . . . , N , h N =θN -1 . Integrating the system (10), [START_REF] Kostousova | On Tight Polyhedral Estimates for Reachable Sets of Linear Differential Systems[END_REF] from right to left, let us, for each τ ∈ T N , solve the optimization problem which is to maximize n i=1 α i (τ, P -(τ ); J) over all possible permutations J = {j 1 , . . . , j n }. This is equivalent to finding the maximal possible velocity of increasing (from right to left) det P -(τ ) (therefore vol P -(τ )) at time τ , by the choice of the value J, when the value P -(τ ) has already been found. Thus we can sequentially construct the piecewise constant function

J(t) ≡ J(τ k ) ∈ Argmax J n i=1 α i (τ k , P -(τ k ); J), t ∈ (τ k-1 , τ k ], k = N, .
. . , 1, and find P -(•).

Control synthesis for discrete-time systems

Now let us briefly consider a problem of control synthesis, similar to Problem 3, for discrete-time systems. This is of independent interest and also may be useful for constructing difference schemes for solving the system [START_REF] Kostousova | On Polyhedral Estimates for Reachable Sets of Differential Systems with Bilinear Uncertainty[END_REF], [START_REF] Kostousova | On Tight Polyhedral Estimates for Reachable Sets of Linear Differential Systems[END_REF]. The analog of Problem 2 can be considered in a similar way.

Consider the controlled discrete-time system with a given terminal set M:

x[k] = (A[k] + U [k] + V [k]) x[k-1] + v[k], k = 1, . . . , N, x[N ] ∈ M = P[p f , Pf ], det Pf = 0, (14)
U [k]∈U[k]={U | Abs (U -Ũ [k])≤ Û [k]}, V [k]∈{V | Abs (V -Ṽ [k])≤ V [k]}, (15) v[k] ∈ Q[k] = P[q[k], Q[k]], k = 1, . . . , N. (16
[k] = (A[k] + U [k, x[k-1]] + V [k]) • x[k-1] + v[k], k = 1, . . . , N , with x[0] = x 0 ∈ P -[0] would satisfy x[k] ∈ P -[k], k = 1, . . . , N , whatever are V [•] and v[•]
subjected to [START_REF] Kurzhanskii | On the Problem of the Synthesis of Controls: the Pontryagin Alternative Integral and the Hamilton-Jacobi Equation[END_REF], [START_REF] Kurzhanski | On the Problem of Synthesizing Control Strategies. Evolution Equations and Set-Valued Integration[END_REF]. Introduce a family of such tubes

P -[•].
Let us consider the following system of relations for

P -[k]=P[p -[k], P -[k]]: p -[k-1]=B[k] -1 (p -[k]-q[k]), B[k]=A[k]+ Ũ [k]+ Ṽ [k], k=N, . . . , 1, p -[N ]=p f , (17)
P -[k-1] = H[k, P -[k-1]], k = N, . . . , 1, P -[N ] = Pf , (18)
[k] -1 ) V [k] Abs (p -[k-1] + P ξ) | ξ ∈ E(C)}, γ[k] = (Abs ( P -[k] -1 Q[k]))e, k = N, . . . , 1. (19)
Note that (17) is the system of explicit recurrent relations while ( 18)-( 19) is the system of implicit ones, i.e., for any time step k∈{N, . . . , 1}, we need to solve the system of nonlinear equations with respect to the unknown matrix P =P -[k-1].

Theorem 3. In the system ( 14)-( 16

U [k, x] x = Ũ [k] x -diag α[k, P -[k-1]; J[k]](x-p -[k-1]), k = 1, . . . , N, (20)
(a formula similar to (12) is true), gives a particular solution to Problem 4.

Proof. We give a sketch following the scheme of the proof of Theorem 2 and keeping the similar notation.

Let x[•] corresponds to x[0] = x 0 ∈ P -[0], i.e., x[0] = p -[0] + P -[0]ζ 0 , where ζ 0 ∞ ≤1. Let us represent x[k] in the form x[k] = p -[k] + P -[k]ζ[k], k = 0, . . . , N .
The proof is by induction on the time step k. Let we already have x[k-1] ∈ P -[k-1]. Then it follows from ( 14), ( 17) that

x[k]=p -[k]+B[k] P -[k-1]ζ[k-1]+(U [k, x[k-1]]-Ũ [k]+∆V [k])x[k-1]+v[k]-q[k],
where

∆V [k] = V [k] -Ṽ [k].
Taking into account [START_REF] Polyak | Robust Stability and Control[END_REF], we obtain 18),( 19),( 15),( 16), we have

ζ[k] = P -[k] -1 (B[k] -diag α[k, P -[k-1]]) P -[k-1]ζ[k-1] + c[k, x[k-1]], c[k, x]= P -[k] -1 ∆V [k] x+ P -[k] -1 Q[k]χ[k]. Using (
ζ[k] = diag (e -β[k, P -[k-1]] -γ[k])ζ[k-1] + c[k, x[k-1]]; Abs c[k, x] ≤ β[k, P -[k-1]] + γ[k] for x ∈ P -[k-1]. It is not difficult to see that if ζ[k-1] ∞ ≤ 1 and e -β[k, P -[k-1]] -γ[k] ≥ 0, then ζ[k] ∞ ≤ 1. Thus we obtain the desired inclusion x[k] ∈ P -[k]. Also it is not difficult to see that Abs (U [k, x] -Ũ [k]) ≤ Û [k] for x ∈ P -[k-1].
Remark 3. Let the system ( 14)-( 16) be obtained by the Euler approximations of ( 1), ( 3)-( 5) with the same M, A

[k]=I+h N A(t k-1 ), Ũ [k]=h N Ũ (t k-1 ), Û [k] = h N Û (t k-1 ), Ṽ [k]=h N Ṽ (t k-1 ), V [k]=h N V (t k-1 ), Q[k]=h N Q(t k-1 ), t k = kh N , h N = θN -1 .
Let, for a fixed k, det P -[k] = 0 and the time step h N be sufficient small. Then the operator H[k, P ] is contractive in some domain D

[k] = {P | P - P -[k] ≤ δ[k]}, i.e., H[k, P 1 ]-H[k, P 2 ] ≤ L P 1 -P 2 for any P 1 , P 2 ∈ D[k],
where L=L[k]∈(0, 1), and therefore [3, p. 319] the equation P = H[k, P ] from ( 18), ( 19) has a solution P = P -[k-1], which can be found by the simple iteration P l+1 =H[k, P l ], l=0, 1, . . ., starting from P 0 = P -[k], and we have P l -P ≤ L l (1-L) -1 P 1 -P 0 . Also, the relation γ

[k] + β[k, P -[k-1]
] < e is satisfied. But certainly we can not derive from here the existence of nonsingular matrices P -[k] for all k = N, . . . , 1 because the value of such "small" h N depends on k.

Examples

We consider model examples for Problem 3. For computations we use the Euler approximations (see Remark 3) with N =200. Let M = P((1, 1) , I, (0.1, 0.1) ), A(t)≡ -0.5 0 0 -0.5 , Ũ (t)≡ 0 2 0 0 , Û (t)≡ 0 1.5 0 0 , Ṽ (t)≡ 0 0 2 0 , V (t)≡ 0 0 0.2 0 or V (t)≡0, Q(t) ≡ P(0, I, 0) or Q(t) ≡ P(0, I, (0.05, 0.05) ), θ=0.25. We consider 3 cases: (I), (II,ii), and (II,iii). The results are presented in Fig. 1. In the second example, we put the realization V (t) ≡ Ṽ (t) + V (t); in the third one we presume V (•) to be the same and v(•) to be some extremal bang-bang type disturbance [17, p.234], where the length of intervals of constancy of v(t) is equal to θ/4. All the presented trajectories reach the target set M including the one with a small violation of the inclusion x 0 ∈ P -[0], though if x 0 / ∈ P -[0], then there is not guarantee that the trajectory can be steered into the target set M by using the control strategy (12) under any disturbances.

  The research was partially supported by the Program of Basic Research of the Ural Branch of RAS (Project 12-P-1-1019), by the Russian Foundation for Basic Research (Grant 12-01-00043), and by the Program for the State Support of Leading Scientific Schools of Russian Federation (Grant 2692.2014.1).

) Problem 4 .

 4 Find a polyhedral tube P-[k] = P[p -[k], P -[k]], k = 1, . . . , N , with P -[N ] = M, and find a corresponding feedback control strategy U = U [k, x] such that U [k, x] ∈ U[k] for x ∈ P -[k-1], k = 1, . . ., N , and each solution x[•]to the equation x

H

  [k, P ] = (B[k] -diag α[k, P ]) -1 P -[k] diag (e -β[k, P ] -γ[k]), α i [k, P ] = α i [k, P ; J[k]] = ûji i [k] η ji [k, P ] (e i (Abs P ) e) -1 , i=1, . . . , n, η[k, P ] = max{0, Abs p -[k-1] -(Abs P )e}, β[k, P ] = max{Abs ( P -

  ), let det Pf = 0 and all detB[k] = 0. Let J[k] = {j 1 [k], . . . , j n [k]} be arbitrary permutations of numbers {1, . . . , n}, k=N, . . . , 1, and the system (17)-(19) has a solution (p-[•], P -[•]) such that we obtain det P -[k] = 0 and e -β[k, P -[k-1]] -γ[k] > 0, k = N, . . . , 1.Then the tube P -[•] and the control strategy which acts according to the following rule

Fig. 1 .

 1 Fig. 1. Examples of polyhedral control synthesis for Problem 3 (n = 2). (a) Case (I): the set M (dash line), two parallelotopes P -[0] and the controlled trajectory for x0 = (0.5, 0.8) . (b) Case (I): the tube P -[•] and the controlled trajectory. (c) Case (II,ii): M, P -[0] corresponding to J[•] from Remark 2, and the controlled trajectory for x0 = (0.5, 0.8) . (d) Case (II,iii): M, P -[0] corresponding to J[•] from Remark 2, and controlled trajectories for two initial points x0 = (0.5, 0.8) and x0 = (1, 0.6) .

  Taking into account[START_REF] Kostousova | On Tight Polyhedral Estimates for Reachable Sets of Linear Differential Systems[END_REF] and the relation d dt

and V (•) (which satisfies (4)). Let us represent x(t) -p -(t) in the form x(t) -p -(t) = P -(t)ζ(t). Then we have d dt ζ = -( P -) -1 ( d dt P -)ζ + ( P -) -1 d dt (x -p -) for the above function ζ.


Diana R Merluşcȃ
email: dianam1985@yahoo.com

Application of the Fenchel theorem to the obstacle problem

Keywords: obstacle problem, duality, approximate problem

In this paper we apply a duality algorithm to the general obstacle problem for second order operators. We reduce the problem to the null obstacle case and we solve it by using an algorithm based on a dual approximate problem. This method generates a quadratic minimization problem, which is easy to implement numerically. The convergence properties and the numerical results show that the algorithm is working properly for any admissible obstacle.

Introduction

The obstacle problem is a very well studied subject. Many methods have been applied for solving it. For instance, Glowinski [START_REF] Glowinski | Numerical Methods for Nonlinear Variational Problems[END_REF] used the finite element method for solving the null obstacle problem, while Barbu and Precupanu [START_REF] Barbu | Precupanu T. Convexity and optimization in Banach spaces[END_REF] studied the problem from the duality point of view. The general obstacle problem is also intensively studied for its wide range of applications in mechanics and physics. We refer here to Rodrigues [START_REF] Rodrigues | Obstacle problems in mathematical physics[END_REF], Caffarelli and Friedman [START_REF] Caffarelli | The obstacle problem for the biharmonic operator[END_REF], Duvaut and Lions [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF] and Ciarlet [START_REF] Ciarlet | Numerical analysis of the finite element method[END_REF].

We extend here the the duality method developed in the articles Merluşcȃ [START_REF] Merluşcȃ | A duality algorithm for the obstacle problem[END_REF][START_REF] Merluşcȃ | A duality-type method for the obstacle problem[END_REF], by the application of the Fenchel theorem to the obstacle problem. We discuss the general obstacle problem in Section 2. We reduce the problem to the null obstacle case and we compute the solutions using the duality method (Merluşcȃ [9]). In Merluşcȃ [START_REF] Merluşcȃ | A Duality-Type Method For The Fourth Order Obstacle Problem[END_REF], the case of the fourth order obstacle problem was analysed. In section 3, we apply this technique in numerical examples for one dimensional problems and the obtained results are very accurate. Finally, we mention the works of Neittaanmaki, Sprekels and Tiba [START_REF] Neittaanmaki | Optimization of elliptic systems. Theory and applications[END_REF] and Sprekels and Tiba [START_REF] Sprekels | Sur les arches lipschitziennes (French)[END_REF], where a related duality approach was used in the study of Kirchhoff-Love arches and explicit solutions were obtained.

The general obstacle problem for second order equations

We consider the following obstacle problem

min 1 2 Ω |∇y| 2 - Ω f y : y ∈ K ψ , (1)
where

K ψ = {y ∈ H 1 0 (Ω) : y ≥ ψ}, ψ ∈ H 1 (Ω) is such that ψ| ∂Ω ≤ 0 and f ∈ L 2 (Ω).
It is known that the unique solution of problem ( 1) is an element in H 2 (Ω) (Th. 2.5, [START_REF] Barbu | Optimal control of variational inequalities[END_REF]).

Lemma 1. Let y ψ be the solution of the problem (1) and ŷ the solution of the problem

-∆ŷ = f, on Ω , ŷ = 0, on ∂Ω , (2)
then y ψ ≥ ŷ almost everywhere on Ω.

The problem (1) in which we replace ψ by ψ = max{ŷ, ψ} ∈ H 1 0 (Ω) has the same solution y ψ .

Proof. Denoting β ⊂ R × R a maximal monotone operator defined by

β(z) =    ] -∞, 0], z = 0 , 0, z > 0 , ∅, z < 0 .
we rewrite (1) as

-∆y ψ + β(y ψ -ψ) f in Ω . (3)
Then, since y ψ ∈ H 2 (Ω), β(y ψ -ψ) ∈ L 2 (Ω) and β(y ψ -ψ) ≤ 0 a.e. on Ω. By a comparison of (2) and (3), we obtain that y ψ ≥ ŷ a. e. on Ω.

We denote K = {y ∈ H 1 0 (Ω) : y ≥ ψ}. Then y ψ ∈ K, ∆y ψ + f ≤ 0 a.e. on Ω. For every v ∈ K, we compute

Ω (∆y ψ + f )(v -y ψ ) = Ω (∆y ψ + f )( ψ -y ψ ) + Ω (∆y ψ + f )(v -ψ) ≤ Ω (∆y ψ + f )( ψ -y ψ ) = 0 .
The last equality is due to the classical formulation of the obstacle problem (the complementarity property)

-∆y ψ = f, in Ω + = {y ψ ∈ Ω : y ψ (x) > ψ(x)} , -∆y ψ ≥ f, in Ω \ Ω + = {y ψ ∈ Ω : y ψ (x) = ψ(x)} , y ψ = ψ and ∂y ψ ∂n = ∂ψ ∂n , on ∂Ω + ∩ Ω , y ψ = 0 on ∂Ω .
and to the fact that y ψ (x) = ψ(x) means that ŷ(x) ≤ ψ(x) and that yields that y ψ (x) = ψ(x).

Then integrating by parts we get

Ω ∇y ψ ∇(v -y ψ ) ≤ Ω f (v -y ψ ), ∀v ≥ ψ .
Remark 1. Lemma 1 is known (see Murea and Tiba [START_REF] Murea | A direct algorithm in some free boundary problems[END_REF]) and we indicate its proof for easy reference.

The null obstacle problem that we use is

min y∈K0 1 2 Ω |∇y| 2 - Ω f y + Ω ∇ ψ∇y . (4)
Let y 0 be the unique solution of problem (4).

Proposition 1. Then the solution of the problem (1) can be computed by just adding ψ, i.e.

y ψ = y 0 + ψ . (5)
Proof. The weak formulation of problem ( 4) is given by the form

Ω ∇y 0 ∇(y 0 -v) ≤ Ω f (y 0 -v) - Ω ∇ ψ∇(y 0 -v), ∀v ∈ K 0 .
We translate the problem by adding ψ. Then, for every v ∈ K 0 , we get v + ψ ≥ ψ ≥ ψ. With this translations from the variational inequality we obtain

Ω ∇(y 0 + ψ)(∇(y 0 + ψ) -∇( ψ + v)) ≤ Ω f (y 0 + ψ -v -ψ) .
Using Lemma 1 it yields that, by a translation with ψ, the problem (4) is equivalent to problem (1). Then we conclude that y 0 + ψ = y ψ .

Since Ω ∇y∇ ψ = -Ω ∆ ψy and ∆ ψ ∈ H -1 (Ω), then we can consider the approximate problem

min 1 2 Ω |∇y| 2 - Ω f + ∆ ψ y : y ∈ C k , (6)
where

C k = {y ∈ H 1 0 (Ω) : y(x i ) ≥ 0, ∀i = 1, 2, . . . , k} and {x i } i is a dense set in Ω.
In [START_REF] Glowinski | Numerical Methods for Nonlinear Variational Problems[END_REF] , we assume that dim Ω = 1 (for the numerical applications in Section 3), but the result can be extended in higher dimension by using non hilbertian Sobolev spaces. Using the Sobolev imbedding theorem and the weak lower semicontinuity of the norm, then we can prove the following approximation result (see Merluşcȃ [START_REF] Merluşcȃ | A duality-type method for the obstacle problem[END_REF]) Theorem 1. The sequence {ȳ k } k of the solutions of problems (6), for k ∈ N, is a strongly convergent sequence in H 1 0 (Ω) to the unique solution ȳ of the problem (4).

We denote f = f + ∆ ψ ∈ H -1 (Ω). Applying the Fenchel duality theorem to problem [START_REF] Glowinski | Numerical Methods for Nonlinear Variational Problems[END_REF] we obtain the dual problem

min 1 2 |y * + f | 2 H -1 (Ω) : y ∈ C * k , (7)
where

C * k = {y * ∈ H -1 (Ω) : y * = k i=1 α i δ xi , α i ≥ 0} is the dual cone.
Remark 2. Let y * k be the solution of the dual approximate problem [START_REF] Hatch | Freefem documetation, Laboratoire Jacques-Louis Lions[END_REF]. Since y * k ∈ C * k , it is sufficient to compute the coefficients α * i , due to the formula

y * k = k i=1 α * i δ xi .
The solution y k of the approximate problem ( 6) is computed using the equality y k = J -1 (y * k + f ), where J is the duality mapping J : H 1 0 (Ω) → H -1 (Ω) and we also have

α * i y k (x i ) = 0, ∀i = 1, k.
We obtain the formula for the solution of the approximate problem, denoted by y 0 k ,

y 0 k = k i=1 α * i J -1 (δ i ) + J -1 ( f )
using the fact that the duality mapping J : H 1 0 (Ω) → H -1 (Ω) is defined by J(y) = -∆y.

Then applying [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF] we find the approximate solution of the general obstacle problem (1).

Numerical applications

In this section we discuss two examples in one dimension for the general obstacle problem for second order operators.

We consider the obstacle problem [START_REF] Barbu | Optimal control of variational inequalities[END_REF] with Ω =]-1, 1[ the domain, ψ ≡ -1/18 the obstacle and

f (x) = -1, |x| > 1/4 , 1 -32x 2 , |x| ≤ 1/4 .
The solution of this problem is, Ockendon and Elliott [START_REF] Ockendon | Weak and variational methods for moving boundary problems[END_REF], (pp. 93 -94)

u(x) =        -1 18 + 1 2 x ± 2 3 2 , 2 3 < |x| ≤ 1 , -1 18 , 1 3 ≤ |x| ≤ 2 3 , -1 18 + 1 2 x ± 1 3 2 , 1 4 ≤ |x| < 1 3 , -1 32 + 8 3 x 2 x 2 -3 16 , |x| < 1 4 .
The duality mapping J : H 1 0 (Ω) → H -1 (Ω) is defined as J(y) = -y . It is a linear bounded operator.

We consider the discretization x i = 2ih -1, where h = 1/k, for all i ∈ {0, 1, 2, . . . , k}.

We denote by d i = J -1 (δ xi ). Computing them we get

d i (x) = 0.5(1 -x i )(x + 1), x ≤ x i , 0.5(1 -x i )(x + 1) -x + x i , x > x i ,
∀i ∈ {0, 1, . . . , k} .

We consider y f = J -1 (f + ∆ ψ) the solution of the problem

-y f = f + ∆ ψ, on ] -1, 1[ , y f (-1) = y f (1) = 0 . Then |y * + f + ∆ ψ| 2 H -1 (Ω) = k i=1 α i d i + y f 2 H 1 0 (Ω)
and

k i=1 α i d i + y f 2 H 1 0 (Ω) = k i,j=1 α i α j Ω d i d j + 2 k i=1 α i Ω d i y f + Ω (y f ) 2 .
Denoting Computing the components, we get b i = y f (x i ) and

a ij = Ω d i d j , b i = Ω d i y f ,
a ij = 0.5(1 + x i )(1 -x j ), j > i , 0.5(1 + x j )(1 -x i ), j ≤ i .
In Figure 1 we represent the coefficients {α * i } i=1,100 , the solution of the problem [START_REF] Merluşcȃ | A duality algorithm for the obstacle problem[END_REF].

We construct the solution of the problem (6) using Remark 2, then we apply formula [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF] to obtain the approximate solution of (1).

In Figure 2 we represent three solutions: the one computed by the duality method, the one computed with the IPOPT optimizer (Freefem++ script included in the ff-Ipopt dynamic library, Hatch [START_REF] Hatch | Freefem documetation, Laboratoire Jacques-Louis Lions[END_REF]; for details about the method see Wächter and Biegler [START_REF] Wächter | On the Implementation of a Primal-Dual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming[END_REF]) and the exact solution given by Ockendon and Elliot [START_REF] Ockendon | Weak and variational methods for moving boundary problems[END_REF]. They coincide graphically. We now consider an example with general obstacle: After solving the quadratic minimization problem (8), the solution of which is represented in Figure 3, we compute the solution of the approximate problem (6) using the Remark 2. We represent in Figure 4 the obstacle ψ and the solutions, one computed by the duality method and the other one computed by the IPOPT method [START_REF] Hatch | Freefem documetation, Laboratoire Jacques-Louis Lions[END_REF]. The two solutions are very close and coincide graphically.

min 1 2 Ω |∇y| 2 - Ω f y : y ∈ K ψ , (9)

1 2 α

 12 we solve the dual problem[START_REF] Hatch | Freefem documetation, Laboratoire Jacques-Louis Lions[END_REF] which is equivalent to the quadratic minimization problemmin α≥0 T Aα + b T α ,(8)where A = [a ij ] and b = [b i ].

Fig. 1 :Fig. 2 :

 12 Fig. 1: The coefficients {α * i } i .

Fig. 3 :

 3 Fig. 3: The coefficients {α * i } i .

Fig. 4 :

 4 Fig. 4: The dual solution and the IPOPT solution coincide graphically.

This paper is supported by the Sectorial Operational Programme Human Resources Development (SOP HRD), financed from the European Social Fund and by the Romanian Government under contract number SOP HRD/107/1.5/S/82514


Marius Cornel
email: cornel.murea@uha.fr

Murea

Dan Tiba
email: dan.tiba@imar.ro

A penalization method for the elliptic bilateral obstacle problem

Keywords: obstacle problem, free boundary problems, penalization

In this paper we propose a new algorithm for the wellknown elliptic bilateral obstacle problem. Our approach enters the category of fixed domain methods and solves just linear elliptic equations at each iteration. The approximating coincidence set is explicitly computed. In the numerical examples, the algorithm has a fast convergence.

Introduction

The obstacle problem may be formulated as an elliptic variational inequality. Detailed theoretical discussions of various variational inequalities may be found in [START_REF] Brézis | Problèmes unilatéraux[END_REF][START_REF] Kinderlehrer | An introduction to variational inequalities and their applications[END_REF][START_REF] Rodrigues | Obstacle Problems in Mathematical Physics[END_REF]. Applications, including optimal control problems are investigated in the books [START_REF] Barbu | Optimal control of variational inequalities[END_REF][START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF][START_REF] Elliott | Weak and variational methods for moving boundary problems[END_REF][START_REF] Tiba | Optimal Control of Nonsmooth Distributed Parameter Systems[END_REF]. From the point of view of the numerical approximation, we quote just the monographs [START_REF] Elliott | Weak and variational methods for moving boundary problems[END_REF][START_REF] Glowinski | Numerical analysis of variational inequalities[END_REF][START_REF] Neittaanmaki | Optimal Control of Nonlinear Parabolic Systems[END_REF].

In this paper we propose an algorithm for the elliptic bilateral obstacle problem which is of fixed domain type in the sense that the finite element discretization is given in the whole domain, independently of the position of the unknown free boundary. In each iteration a linear elliptic equation has to be solved in the whole domain and the corresponding stiffness matrix is common for all iterations. This is a clear advantage from the point of view of the implementation and the approximating coincidence set is explicitly computed in each iteration and it converges in the Hausdorff-Pompeiu sense [START_REF] Neittaanmaki | Optimization of elliptic systems[END_REF] to the searched geometry. Moreover, we need just a scalar penalization parameter in our method. A similar strategy was employed in [START_REF] Murea | A direct algorithm in some free boundary problems, submitted to[END_REF] for the elliptic unilateral obstacle problem and for parabolic variational inequalities.

Supported by Grant 145/2011 CNCS, Romania.

Our approach is inspired from shape optimization techniques, but no shape optimization problem is used here although this is a known method in free boundary problems, [START_REF] Burger | A level set based shape optimization method for an elliptic obstacle problem[END_REF]. One may compare the present approach to the recent works [START_REF] Halanay | Existence and approximation for a steady fluid-structure interaction problem using fictitious domain approach with penalization[END_REF][START_REF] Neittaanmaki | Fixed domain approaches in shape optimization problems with Dirichlet boundary conditions[END_REF][START_REF] Neittaanmaki | Fixed domain approaches in shape optimization problems[END_REF]. An efficient Lagrangian method together with a primal-dual active set strategy with regularization is studied in [START_REF] Ito | Semi-smooth Newton methods for variational inequalities of the first kind[END_REF]. But our approach and arguments are certainly diferent. We also quote the multi grid method employed in [START_REF] Hoppe | Adaptive multilevel methods for obstacle problems[END_REF], the path-following method for semi-smooth Newton schemes [START_REF] Hintermuller | Path-following methods for a class of constrained minimization problems in function space[END_REF] and a duality-type method [START_REF] Merlusca | A duality algorithm for the obstacle problem[END_REF][START_REF] Merlusca | A duality-type method for the fourth order obstacle problem[END_REF].

Formulation of the problem and the algorithm

Let D be a smooth domain in R d , d ∈ N * and f ∈ L 2 (D) be given. We denote the obstacles by

ψ 1 , ψ 2 : D → R, such that ψ 1 , ψ 2 ∈ H 2 (D), ψ 1 ≤ ψ 2 in D, ψ 1|∂D ≤ 0, ψ 2|∂D ≥ 0. The admisible set K = v ∈ H 1 0 (D); ψ 1 (x) ≤ v(x) ≤ ψ 2 (x) a.e. in D
is a nonvoid closed convex subset of H 1 0 (D). To K, the following variational inequality, may be associated:

D ∇y • (∇y -∇v) dx ≤ D f (y -v) dx, ∀v ∈ K. (1)
The existence of a unique solution y ∈ K is wellknown. We introduce β ⊂ R × R the maximal monotone graph given by

β(r) =    ∅, r < 0, ] -∞, 0], r = 0, 0, r > 0, (2)
the maximal monotone graph γ ⊂ R × R given by

γ(r) =    0, r < 0, [0, +∞[ , r = 0, ∅, r > 0 (3)
and denote by β , γ , > 0, their Yosida approximations. We have

β (r) = 1 r, r ≤ 0, 0, r > 0, γ (r) = 0, r < 0, 1 r, r ≥ 0.
Notice that β is a concave and γ is a convex function in R.

In the case when

f ∈ L 2 (D), ψ 1 , ψ 2 ∈ H 2 (D) with the compatibility condi- tion ψ 1 ≤ ψ 2 in D, ψ 1|∂D ≤ 0, ψ 2|∂D ≥ 0, it
is known that the solution of (1) satisfies the regularity property y ∈ H 2 (D). Moreover, in this case, the obstacle problem may be written as a multivalued equation

-∆y + β(y -ψ 1 ) + γ(y -ψ 2 ) f in D. (4)
One can define two coincidence sets, corresponding to the two obstacles:

D 1 = {x ∈ D; y(x) = ψ 1 (x)} D 2 = {x ∈ D; y(x) = ψ 2 (x)}
and associated to [START_REF] Barbu | Optimal control of variational inequalities[END_REF].

We state now our algorithm.

Algorithm 1) Choose n = 0, 0 > 0, Ω 0 1 ⊂ D, Ω 0 2 ⊂ D open subsets such that (D \ Ω 0 1 ) ∩ (D \ Ω 0 2 ) = ∅, ỹ-1 = 0; 2) Compute y n ∈ H 1 0 (D)
as solution of the linear elliptic equation

-∆y n + 1 n χ D\Ω n 1 (y n -ψ 1 ) + 1 n χ D\Ω n 2 (y n -ψ 2 ) = f in D (5) 3) Compute y n = min {ψ 2 , max{y n , ψ 1 }}, Ω n+1 1 = {x ∈ D; y n (x) > ψ 1 (x)}, Ω n+1 2 = {x ∈ D; y n (x) < ψ 2 (x)} n+1 = n 2 ; 4) If y n -y n-1 H 1 (D) < tol then STOP else n=n+1 GO TO step 2.
Remark 1. By the classical result of [START_REF] Brézis | Équivalence de deux inéquations variationnelles et applications[END_REF], the elastic-plastic torsion problem is equivalent with a variational inequality of obstacle type and our algorithm may be applied as well.

We convene to extend the value 1 for β and γ in the origin and we can rewrite the step 2 of the Algorithm as

-∆y n + β n (y n-1 -ψ 1 ) (y n -ψ 1 ) + γ n (y n-1 -ψ 2 ) (y n -ψ 2 ) = f. (6)
Recall that the usual approximation by regularization of the variational inequality (1) is

-∆ỹ n + β n (ỹ n -ψ 1 ) + γ n (ỹ n -ψ 2 ) = f in D, (7)
plus homogeneous boundary conditions on ∂D. Notice that β (r) = β (r)r and γ (r) = γ (r)r, under the above convention, which shows that ( 6) and ( 7) have very similar structure. Clearly, ( 7) is a nonlinear elliptic equation, while the decoupling operated in (6) allows to use linear elliptic equations.

Stability

We present in this section a stability result in L 2 (D) for the algorithm introduced above applied to the bilateral obstacle problem.

Theorem 1. i) The sequence {y n } is bounded in L 2 (D).
ii) There is C > 0, independent of n, such that:

D\Ω n 2 (y n -ψ 2 ) 2 + dx + D\Ω n 1 (y n -ψ 1 ) 2 -dx ≤ C n . (8)
Proof. Using β (r) = β (r)r, the concavity of β (•) and the definition of the subdifferential of concave mapping, we obtain

β n (y n-1 -ψ) (y n -ψ) = β n (y n-1 -ψ) (y n-1 -ψ) + β n (y n-1 -ψ) (y n -ψ -y n-1 + ψ) ≥ β n (y n-1 -ψ) +β n (y n -ψ) -β n (y n-1 -ψ) = β n (y n -ψ) .
We use the above inequality in the equation ( 5). We get

-∆y n + β n (y n -ψ 1 ) + 1 n χ D\Ω n 2 (y n -ψ 2 ) ≤ f, (9)
where β is given by ( 2) and β n is its regularization. We multiply ( 9) by (y n -ψ 2 ) + and we use that

β n (y n -ψ 1 )(y n -ψ 2 ) + = 0. ( 10
)
While β n (y n -ψ 1 ) may take negative values, this happens for

y n ≤ ψ 1 , that is y n -ψ 2 ≤ 0 (since ψ 1 ≤ ψ 2 )
. Then (y n -ψ 2 ) + = 0 and (10) follows. We infer

D |∇(y n -ψ 2 ) + | 2 + 1 n D\Ω n 2 (y n -ψ 2 ) 2 + (11)
≤ D f (y n -ψ 2 ) + + D ∇ψ 2 • ∇(y n -ψ 2 ) + .
By the conditions ψ 2 |∂D ≥ 0 and y n |∂D = 0, we have (y n -ψ 2 ) + = 0 on ∂D and the Poincaré inequality shows that {(y n -ψ 2 ) + } is bounded in H 1 0 (D), by [START_REF] Hintermuller | A mesh-independence result for semismooth Newton methods[END_REF].

Equation ( 5) may be rewritten in the form

-∆y n + 1 n χ D\Ω n 1 (y n -ψ 1 ) + γ n (y n-1 -ψ 2 )(y n -ψ 2 ) = f. (12)
We compute

γ n (y n-1 -ψ 2 )(y n -ψ 2 ) = γ n (y n-1 -ψ 2 )(y n-1 -ψ 2 ) ( 13
) +γ n (y n-1 -ψ 2 )(y n -ψ 2 -y n-1 + ψ 2 ) = γ n (y n-1 -ψ 2 ) +γ n (y n-1 -ψ 2 )(y n -ψ 2 -y n-1 + ψ 2 ) ≤ γ n (y n-1 -ψ 2 ) +γ n (y n -ψ 2 ) -γ n (y n-1 -ψ 2 ) = γ n (y n -ψ 2 )
using the subdifferential of convex mappings.

By ( 12), ( 13), we obtain

-∆y n + 1 n χ D\Ω n 1 (y n -ψ 1 ) + γ n (y n -ψ 2 ) ≥ f. (14)
Multiply ( 14) by -(y n -ψ 1 ) -∈ H 1 0 (D), due to y n = 0 on ∂D, ψ 1 ≤ 0 on ∂D:

D |∇(y n -ψ 1 ) -| 2 + 1 n D\Ω n 1 (y n -ψ 1 ) 2 --γ n (y n -ψ 2 )(y n -ψ 1 ) -(15) ≤ - D f (y n -ψ 1 ) -+ D ∇ψ 1 • ∇(y n -ψ 1 ) -. Notice that -γ n (y n -ψ 2 )(y n -ψ 1 ) -= 0 (16)
since -γ n (y n -ψ 2 ) may take negative values just for y n ≥ ψ 2 ≥ ψ 1 and in this case (y n -ψ 1 ) -= 0. By ( 15), ( 16) we obtain:

D |∇(y n -ψ 1 ) -| 2 + 1 n D\Ω n 1 (y n -ψ 1 ) 2 (17) ≤ - D f (y n -ψ 1 ) -+ D ∇ψ 1 • ∇(y n -ψ 1 ) -.
Relation [START_REF] Merlusca | A duality-type method for the fourth order obstacle problem[END_REF] shows that {(y n -ψ 1 ) -} is bounded in H 1 0 (D), by the Poincaré inequality.

We use the inequality

(x -b) + ≤ (x -a) + + (a -b) +
and we have

(y n -ψ 1 ) + ≤ (y n -ψ 2 ) + + (ψ 2 -ψ 1 ) + = (y n -ψ 2 ) + + ψ 2 -ψ 1 .
Relation [START_REF] Hintermuller | A mesh-independence result for semismooth Newton methods[END_REF] shows that {(y n -ψ 1 ) + } is bounded in L 2 (D). In combination with [START_REF] Merlusca | A duality-type method for the fourth order obstacle problem[END_REF], it yields {y n } bounded in L 2 (D). Relation [START_REF] Halanay | Existence and approximation for a steady fluid-structure interaction problem using fictitious domain approach with penalization[END_REF] follows by adding ( 11) and ( 17) and using the already established boundedness of all the terms except the penalization term. This ends the proof.

Remark 2. In fact, the above proof shows that {y n } bounded in L p (D), p > 2 depending on the dimension of D. Relation (8) says that the sequence {y n } does not overpass the obstacles ψ 1 , ψ 2 , in the limit. The proof also provides partial information on {∇y n }, but it is unclear whether {y n } is bounded in H 1 0 (D).

Numerical tests

We have used the software FreeFem++ v 3.19, [9]. For all tests, we use the same initial guess for the coincidence set D \ Ω 0 1 = ∅ and D \ Ω 0 2 = ∅.

Test 1. We consider the torsion of an elastic-plastic prism studied in [START_REF] Glowinski | Numerical analysis of variational inequalities[END_REF] p. 133 and [START_REF] Xue | An algorithm for solving the obstacle problems[END_REF] We use a mesh of 39158 triangles, 19836 vertices and size h = 1 128 . The tolerance for the stopping test is tol = 10 -3 and the penalization parameter is Our algorithm stops after 6 iterations and the relative error in the H 1 norm at the last iteration is y n -y n-1 H 1 (D) = 1.5 × 10 -5 .

In [START_REF] Murea | A direct algorithm in some free boundary problems, submitted to[END_REF], we have tested numerically with positive results the stability of a similar algorithm when f the right-hand side in (1) is perturbed.

Test 2.

We solve the problem [START_REF] Barbu | Optimal control of variational inequalities[END_REF] where

K = {v ∈ H 1 0 (D); ψ 1 (x) ≤ v(x) ≤ ψ 2 (x) a.e. in D}, D = [0, 1] × [0, 1], ψ 1 (x) = -dist(x, ∂D), ψ 2 (x) = dist(x, ∂D) and f (x) = 11(x + y -1)
. Now both obstacles are active. We use a mesh of 39158 triangles, 19836 vertices, the size h = 1 128 , the tolerance for the stopping test tol = 10 -3 and the penalization parameter is n = 0.003. The algorithm stops in 4 iterations and the relative error in the H 1 norm at the last iteration is y n -y n-1 H 1 (D) = 0.000209.

The coincidence sets are presented in Figure 2 and the computed solution in Figure 3. We solved the problem on different meshes, see Table 1. We denote by u i , the solution obtained using the mesh no. i. In [START_REF] Hintermuller | A mesh-independence result for semismooth Newton methods[END_REF], for the semi-smooth Newton method, it is proved a mesh-independence result: the continuous and the discrete process, converge q-linearly with the same rate.

  . The cross-section of the prism is D = [0, 1]×[0, 1]. We solve the problem (1) with f (x) = -8 where K = {v ∈ H 1 0 (D); -1 ≤ ∇v(x) ≤ 1 a.e. in D}. For v ∈ K, the set {x ∈ D; |∇v(x)| < 1} is the elastic zone and {x ∈ D; |∇v(x)| = 1} is the plastic zone, [5, p. 264]. By the result of [3], the elastic-plastic torsion problem is equivalent with a variational inequality of obstacle type K = {v ∈ H 1 0 (D); ψ 1 (x) ≤ v(x) ≤ ψ 2 (x) a.e. in D} where ψ 1 (x) = -dist(x, ∂D) and ψ 2 (x) = dist(x, ∂D). If f < 0, then y ≤ 0, consequently the top obstacle will be inactive.

Fig. 1 .

 1 Fig. 1. Test 1. The computed coincidence set of the plastic zone for the bottom obstacle at the first (left), second (middle) and last (right) iteration.

Fig. 2 .

 2 Fig. 2. Test 2. Coincidence sets of the plastic zone for the bottom obstacle (left) and for the top obstacle (right).

Fig. 3 .

 3 Fig. 3. Test 2. Computed solution.

Fig. 4 .

 4 Fig. 4. Test 3. Computed solution.

Fig. 5 .

 5 Fig. 5. Test 3. Coincidence sets for the bottom obstacle (left) and for the top obstacle (right).

n = 0.003. The coincidence set of the solution presented in Figure1at the right is similar to the above references.

mesh no. mesh size h triangles vertices uiui-1 H 1 (D) Test 3. Now we test the algorithm for the torsion of the elastic-plastic prism discussed in [START_REF] Karkkainen | Augmented Lagrangian active set methods for obstacle problems[END_REF] and [START_REF] Wang | An algorithm for solving the double obstacle problems[END_REF]. We can put this problem in the form [START_REF] Barbu | Optimal control of variational inequalities[END_REF].

We use a mesh of 39158 triangles, 19836 vertices, the size h = 1 128 , the tolerance for the stopping test tol = 10 -3 and the penalization parameter is n = 0.03. The computed solution after 6 iterations is presented in Figure 4 and the corresponding coincidence sets in Figure 5. The relative error in the H 1 norm at the last iteration is y n -y n-1 H 1 (D) = 2.7 × 10 -5 .

In [START_REF] Karkkainen | Augmented Lagrangian active set methods for obstacle problems[END_REF], an augmented lagrangian active set strategy is employed. At each iteration, a reduced linear system associated with the inactive set is solved. In [START_REF] Wang | An algorithm for solving the double obstacle problems[END_REF], at each iteration, linear systems associated to the complementary of the coincidence sets are solved.


Salah-Eddine Rebiai

Fatima Zohra Sidi

Ali Ltm

Exponential stability of compactly coupled wave equations with delay terms in the boundary feedbacks

Keywords: Coupled wave equations, time delays, boundary stabilization

We consider a linear system of compactly coupled wave equations with Neumann feedback controllers that contain delay terms. First, we prove under some assumptions that the closed-loop system generates a C0-semigroup of contractions on an appropriate Hilbert space. Then, under further assumptions, we show that the closed-loop system is exponentially stable. This result is obtained by introducing a suitable energy function and by using an observability estimate.

Introduction

In [START_REF] Datko | Not all feedback stabilized hyperbolic systems are robust with respect to small time delays in their feedbacks[END_REF] and [START_REF] Datko | An example on the effect of time delays in boundary feedback stabilization of wave equations[END_REF], Datko et al presented examples of infinite-dimensional secondorder systems that become unstable when arbitrary small time delays occur in the damping.

Xu et al established in [START_REF] Xu | Stabilization of wave systems with input delay in the boundary control[END_REF] sufficient conditions that guarantee the exponential stability of the one-dimensional wave equation with a delay term in the boundary feedback. Nicaise and Pignotti [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF] extended this result to the multi-dimensional wave equation with a delay term in the boundary or internal feedbacks. The same type of result was obtained by Nicaise and Rebiai [START_REF] Nicaise | Stabilization of the Schrödinger equation with a delay term in boundary feedback or internal feedback[END_REF] for the Schrödinger equation.

Motivated by the references [START_REF] Xu | Stabilization of wave systems with input delay in the boundary control[END_REF], [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF], [START_REF] Komornik | Boundary stabilization of compactly wave equations[END_REF] and [START_REF] Najafi | Stabilizability of coupled wave equations in parallel under various boundary conditions[END_REF], we investigate in this paper the problem of exponential stability for a linear system of compactly coupled wave equations with delay terms in the boundary feedbacks.

Let Ω be an open bounded domain of R n with a boundary Γ of class C 2 which consists of two non-empty parts Γ 1 and Γ 2 such that Γ 1 ∩ Γ 2 = ∅. Furthermore, assume that there exists a real vector field h ∈ (C 2 (Ω)) n such that: (H.1) The Jacobian matrix J of h satisfies

Ω J(x)ζ(x).ζ(x)dΩ ≥ c Ω |ζ(x)| 2 dΩ, for some constant c > 0 and for all ζ ∈ L 2 (Ω; R n ), (H.2) h(x).ν(x) ≤ 0 on Γ 1 ,
where ν is the unit normal on Γ pointing towards the exterior of Ω.

Consider the following coupled system of two wave equations with delay terms in the boundary conditions:

∂ 2 u(x, t) ∂t 2 -∆u(x, t) + l(u(x, t) -v(x, t)) = 0 in Ω × (0, +∞), (1)
∂ 2 v(x, t) ∂t 2 -∆v(x, t) + l(v(x, t) -u(x, t)) = 0 in Ω × (0, +∞), (2)
u(x, 0) = u 0 (x), ∂u(x, 0) ∂t = u 1 (x) in Ω, (3)
v(x, 0) = v 0 (x), ∂v(x, 0) ∂t = v 1 (x) in Ω, (4)
u(x, t) = v(x, t) = 0 on Γ 1 × (0, +∞), ( 5
) ∂u(x, t) ∂ν = -α 1 ∂u(x, t) ∂t -α 2 ∂u(x, t -τ ) ∂t on Γ 2 × (0, +∞), (6)
∂v(x, t) ∂ν = -β 1 ∂v(x, t) ∂t -β 2 ∂v(x, t -τ ) ∂t on Γ 2 × (0, +∞), (7)
∂u(x, t -τ ) ∂t = g(x, t -τ ) on Γ 2 × (0, τ ), ( 8)

∂v(x, t -τ ) ∂t = h(x, t -τ ) on Γ 2 × (0, τ ). (9)
Physically, u and v may represent the displacements of two vibratings objects measured from their equilibrium positions, the coupling terms ±l(u -v) are the distributed springs linking the two vibrating objects. l, α 1 , α 2 , β 1 , β 2 are positive constants, τ is the time delay, u 0 , u 1 , v 0 , v 1 , g and h are the initial data.

It is well known that in the absence of delay (i.e. α 2 = β 2 = 0), the solution of ( 1)-( 9) with α 1 and β 1 positive, decays exponentially to zero in the energy space [START_REF] Najafi | Stabilizability of coupled wave equations in parallel under various boundary conditions[END_REF] and [START_REF] Komornik | Boundary stabilization of compactly wave equations[END_REF]). The purpose of this paper is to investigate the uniform exponential stability of system (1)-( 9) in the case where all the boundary damping coefficients α 1 , α 2 , β 1 and β 2 are positive. To this end, assume as in [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF] that

H 1 Γ1 (Ω) × L 2 (Ω)× H 1 Γ1 (Ω) × L 2 (Ω) (see
α 1 > α 2 , β 1 > β 2 (10)
and define the energy of a solution of ( 1) -( 9) by

E(t) = 1 2 Ω [|∇u(x, t)| 2 + ∂u(x, t) ∂t 2 + |∇v(x, t)| 2 + ∂v(x, t) ∂t 2 + l |u(x, t) -v(x, t)| 2 ]dx + 1 2 Γ2 1 0 [µ ∂u(x, t -τ ρ) ∂t 2 + ξ ∂v(x, t -τ ρ) ∂t 2 ]dρ dΓ (11) where τ α 2 < µ < τ (2α 1 -α 2 ) (12) and τ β 2 < ξ < τ (2β 1 -β 2 ) (13)
We show that if {Ω, Γ 1 , Γ 2 } satisfies (H.1) and (H.2), then there is an exponential decay rate for E(t).The proof of this result is based on Carleman estimates for a system of coupled nonconservative hyperbolic systems established by Lasieka and Triggiani in [START_REF] Lasiecka | Carleman estimates and exact boundary controllability for a system of coupled non-conservative second-order hyperbolic equations[END_REF] and on compactness-uniqueness arguments.

The main result of this paper can be stated as follows.

Theorem 1. Assume (H1), (H.2), ( 10), ( 12) and (13). Then there exist constants M ≥ 1 and ω > 0 such that

E(t) ≤ M e -ωt E(0).
Theorem 1 is proved in Section 3. In Section 2, we study the well-posedness of system (1) -(9) using semigroup theory.

2 Well-posedness of system (1) -( 9)

Inspired from [START_REF] Nicaise | Stability and instability results of the wave equation with a delay term in the boundary or internal feedbacks[END_REF] and [START_REF] Nicaise | Stabilization of the Schrödinger equation with a delay term in boundary feedback or internal feedback[END_REF], we introduce the auxilliary variables

y(x, ρ, t) = ∂u(x, t -τ ρ) ∂t z(x, ρ, t) = ∂v(x, t -τ ρ)
∂t With these new unknowns, system (1) -( 9) is equivalent to

∂ 2 u(x, t) ∂t 2 -∆u(x, t) + l(u(x, t) -v(x, t)) = 0 in Ω × (0, +∞), (14)
∂y(x, ρ, t) ∂t + 1 τ ∂y(x, ρ, t) ∂ρ = 0 on Γ 2 × (0, 1) × (0, +∞), (15)

∂ 2 v(x, t) ∂t 2 -∆v(x, t) + l(v(x, t) -u(x, t)) = 0 in Ω × (0, +∞), ( 16
) ∂z(x, ρ, t) ∂t + 1 τ ∂z(x, ρ, t) ∂ρ = 0 on Γ 2 × (0, 1) × (0, +∞), (17) u(x, t) = v(x, t) = 0 on Γ 1 × (0, +∞), (18) ∂u(x, t) ∂ν = -α 1 ∂u(x, t) ∂t -α 2 y(x, 1, t) on Γ 2 × (0, +∞), ( 19
) ∂v(x, t) ∂ν = -β 1 ∂u(x, t) ∂t -β 2 z(x, 1, t) on Γ 2 × (0, +∞), (20) y
(x, 0, t) = ∂u(x, t) ∂t , z(x, 0, t) = ∂v(x, t) ∂t on Γ 2 × (0, +∞), ( 21
) u(x, 0) = u 0 (x), ∂u(x, 0) ∂t = u 1 (x) in Ω, ( 22
) v(x, 0) = v 0 (x), ∂v(x, 0) ∂t = v 1 (x) in Ω, (23)
y(x, ρ, 0) = g(x, -τ ρ), z(x, ρ, 0) = h(x, -τ ρ) on Γ 2 × (0, 1). ( 24
)
Denote by H the Hilbert space

H =H 1 Γ1 (Ω) × L 2 (Ω) × L 2 (Γ 2 × L 2 (0, 1)) × H 1 Γ1 (Ω) × L 2 (Ω) × L 2 (Γ 2 × L 2 (0, 1))
where

H 1 Γ1 (Ω) = {u ∈ H 1 (Ω) : u = 0 on Γ 1 } We equip H with the inner product         ζ η θ φ χ ψ         ;          ζ η θ φ χ ψ          = Ω (∇ζ(x).∇ ζ(x) + η(x) η(x)) dx+ µ Γ2 1 0 θ(x, ρ) θ(x, ρ)dρ dΓ + Ω (∇φ(x).∇ φ(x) + χ(x) χ(x)) dx+ ξ Γ2 1 0 ψ(x, ρ) ψ(x, ρ)dρ dΓ + l Ω (ζ(x) -φ(x))( ζ(x) -φ(x))dx
Define in H a linear operator A by

D(A) = {(ζ, η, θ, φ, χ, ψ) T ∈ H 2 (Ω) × H 1 Γ1 (Ω) × L 2 (Γ 2 × H 1 (0, 1))× H 2 (Ω) × H 1 Γ1 (Ω) × L 2 (Γ 2 × H 1 (0, 1)); ∂ζ ∂ν = -α 1 η -α 2 θ(., 1), η = θ(., 0) on Γ 2 ; ∂φ ∂ν = -β 1 χ -β 2 ψ(., 1), χ = ψ(., 0) on Γ 2 } (25) A(ζ, η, θ, φ, χ, ψ) T = (η, ∆ζ + lφ -lζ, -τ -1 ∂θ ∂ρ , χ, ∆φ -lφ + lζ, -τ -1 ∂ψ ∂ρ ) T (26)
Then we can rewrite ( 14) -(24) as an abstract Cauchy problem in

H d dt W (t) = AW (t) W (0) = W 0 (27)
where

W (t) = (u(x, t), ∂u(x, t) ∂t , y(x, ρ, t), v(x, t), ∂v(x, t) ∂t , z(x, ρ, t)) T ,
and

W 0 = (u 0 , u 1 , g(., -.τ ), v 0 , v 1 , h(., -.τ )) T
We verify that A is dissipative and that λI-A is onto for a fixed λ > 0. Thus, by the Lumer-Phillips Theorem (see for instance [START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF]) A generates a strongly continuous semigroup on H and consequently we have Proposition 1. For every W 0 ∈ H, problem (27) has a unique solution W whose regularity depends on the the initial datum W 0 as follows:

W (.) ∈ C([0, +∞); H) if W 0 ∈ H, W (.) ∈ C 1 ([0, +∞); H) ∩ C([0, +∞); D(A)) if W 0 ∈ D(A).
3 Proof of Theorem 1

We prove Theorem 1 for smooth initial data. The general case follows by a standard density argument. We proceed in several steps.

Step 1.

Differentiating E(t) with respect to time, we obtain

d dt E(t) ≤ -k Γ2 { ∂u(x, t) ∂t 2 + ∂u(x, t -τ ) ∂t 2 + ∂v(x, t) ∂t 2 + ∂v(x, t -τ ) ∂t 2 }dΓ (28) where k = min{α 1 - α 2 2 - µ 2τ , µ 2τ - α 2 2 , β 1 - β 2 2 - ξ 2τ , ξ 2τ - β 2 2 }
Step 2.

We rewrite

E(t) = E(t) + E d (t)
where

E(t) = 1 2 Ω {|∇u(x,
From (29), we obtain (here and throughout the rest of the paper C is some positive constant different at different occurences)

E d (t) ≤ C T 0 Γ2 { ∂u(x, s -τ ) ∂t 2 + ∂v(x, s -τ ) ∂t 2 }dΓ ds (30)
for 0 ≤ t + τ ≤ T and T large enough.

Step 3.

From Poincaré inequality and Proposition 3.5 of [START_REF] Lasiecka | Carleman estimates and exact boundary controllability for a system of coupled non-conservative second-order hyperbolic equations[END_REF], we have for T sufficiently large and for any > 0

E(0) ≤ C T 0 Γ2 { ∂u(x, t) ∂ν 2 + ∂u(x, t) ∂t 2 + ∂v(x, t) ∂ν 2 + ∂v(x, t) ∂t 2 }dΓ dt + C{ u 2 L 2 (0,T ;H 1/2+ (Ω)) + v 2 L 2 (0,T ;H 1/2+ (Ω)) } (31)
Inserting the boundary conditions ( 6) and ( 7) into (31), we obtain

E(0) ≤ C T 0 Γ2 { ∂u(x, t) ∂t 2 + ∂u(x, t -τ ) ∂t 2 + ∂v(x, t) ∂t 2 + ∂v(x, t -τ ) ∂t 2 }dΓ dt+ C{ u 2 L 2 (0,T ;H 1/2+ (Ω)) + v 2 L 2 (0,T ;H 1/2+ (Ω)) } (32)
Step

Step 5.

We drop the lower order terms on the right-hand side of (33) by a compactnessuniqueness argument to obtain

E(0) ≤ C T 0 Γ2 { ∂u(x,
E(T ) ≤ Ck -1 1 + Ck -1 E(0) (35)
The desired conclusion follows now from (35).