Mahdi Achache
email: mahdi.achache@math.u-bordeaux.fr

El Maati Ouhabaz
email: elmaati.ouhabaz@math.u-bordeaux.fr.

Non-autonomous right and left multiplicative perturbations and maximal regularity

Keywords: Maximal regularity, non-autonomous evolution equations, multiplicative and additive perturbations. Mathematics Subject Classification (2010): 35K90, 35K45, 47D06

We consider the problem of maximal regularity for non-autonomous Cauchy problems

In both cases, the time dependent operators A(t) are associated with a family of sesquilinear forms and the multiplicative left or right perturbations B(t) as well as the additive perturbation P (t) are families of bounded operators on the considered Hilbert space. We prove maximal L p -regularity results and other regularity properties for the solutions of the previous problems under minimal regularity assumptions on the forms and perturbations.

Introduction

The present paper deals with maximal L p -regularity for non-autonomous evolution equations in the setting of Hilbert spaces. Before explaining our results we introduce some notations and assumptions. Let (H, (•, •), • ) be a Hilbert space over R or C. We consider another Hilbert space V which is densely and continuously embedded into H. We denote by V ′ the (anti-) dual space of V so that

V ֒→ d H ֒→ d V ′ .
We denote by , the duality V-V ′ and note that ψ, v = (ψ, v) if ψ, v ∈ H. We consider a family of sesquilinear forms

a : [0, τ ] × V × V → C such that • [H1]: D(a(t)) = V (constant form domain), • [H2]: |a(t, u, v)| ≤ M u V v V (uniform boundedness),
• [H3]: Re a(t, u, u)+ν u 2 ≥ δ u 2 V (∀u ∈ V) for some δ > 0 and some ν ∈ R (uniform quasi-coercivity).

Here and throughout this paper, • V denotes the norm of V.

To each form a(t) we can associate two operators A(t) and A(t) on H and V ′ , respectively. Recall that u ∈ H is in the domain D(A(t)) if there exists h ∈ H such that for all v ∈ V: a(t, u, v) = (h, v). We then set A(t)u := h. The operator A(t) is a bounded operator from V into V ′ such that A(t)u = a(t, u, •). The operator A(t) is the part of A(t) on H. It is a classical fact that -A(t) and -A(t) are both generators of holomorphic semigroups (e -rA(t) ) r≥0 and (e -rA(t) ) r≥0 on H and V ′ , respectively. The semigroup e -rA(t) is the restriction of e -rA(t) to H. In addition, e -rA(t) induces a holomorphic semigroup on V (see, e.g., Ouhabaz [START_REF] Ouhabaz | Analysis of Heat Equations on Domains[END_REF]Chapter 1]).

A well known result by J.L. Lions asserts that the Cauchy problem

u ′ (t) + A(t)u(t) = f (t), u(0) = u 0 ∈ H (1.1)
has maximal L 2 -regularity in V ′ , that is, for every f ∈ L 2 (0, τ ; V ′ ) there exists a unique u ∈ W 1 2 (0, τ ; V ′ ) which satisfies (1.1) in the L 2 -sense. The maximal regularity in H is however more interesting since when dealing with boundary value problems one cannot identify the boundary conditions if the Cauchy problem is considered in V ′ . The maximal regularity in H is more difficult to prove. J.L. Lions has proved that this is the case for initial data u 0 ∈ D(A(0)) under a quite restrictive regularity condition, namely t → a(t, g, h) is C 2 (or C 1 if u 0 = 0). It was a question by him in 1961 (see [START_REF] Lions | Équations Différentielles Opérationnelles et Problèmes aux Limites[END_REF] p. 68) whether maximal L 2 -regularity holds in general in H.

A lot of progress have been made in recent years on this problem. It was proved by Ouhabaz and Spina [START_REF] Ouhabaz | Maximal regularity for non-autonomous Schrödinger type equations[END_REF] that maximal L p -regularity holds in H if t → a(t, g, h) is C α for some α > 1/2 (for all g, h ∈ V). This result is however proved for the case u 0 = 0 only. In Haak and Ouhabaz [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF], it is proved that for u 0 ∈ (H, D(A(0))) 1-1 p ,p and |a(t, g, h)a(s, g, h

)| ≤ ω(|t -s|) h V g V (1.2)
for some non-decreasing function ω such that

τ 0 ω(t) t 3 2
dt < ∞ and

τ 0 ω(t) t p dt < ∞, (1.3)
then the Cauchy problem (1.1) has maximal L p -regularity in H. The condition (1.3) can be improved if (1.2) holds with norms in some complex interpolation spaces (see Arendt and Monniaux [START_REF] Arendt | Maximal regularity for non-autonomous Robin boundary conditions[END_REF] and Ouhabaz [START_REF] Ouhabaz | Maximal regularity for non-autonomous evolution equations governed by forms having less regularity[END_REF]). It was observed by Dier [7] that the answer to Lions' problem is negative in general. His example is based on non-symmetric forms for which the Kato square root property D(A(t)) 1/2 ) = V is not satisfied. Recently, Fackler [START_REF] Fackler | Lions' problem concerning maximal regularity of equations governed by non-autonomous forms[END_REF] proved a negative answer to the maximal regularity problem for forms which are C α for any α < 1/2 (even symmetric ones). Let us also mention a recent positive result of Dier and Zacher [START_REF] Dier | Non-autonomous maximal regularity in Hilbert spaces[END_REF] on maximal L 2 -regularity in which the condition (1.3) is replaced by a norm in a Sobolev space of order > 1 2 . For forms associated with divergence form elliptic operators, Auscher and Egert [START_REF] Auscher | On non-autonomous maximal regularity for elliptic operators in divergence form[END_REF] proved that the order of this Sobolev space can be 1 2 . One of the aims of the present paper is to study the same problem for multiplicative perturbations. More precisely, we study maximal L pregularity for

u ′ (t) + B(t)A(t)u(t) + P (t)u(t) = f (t), u(0) = u 0 (1.4)
and also for

u ′ (t) + A(t)B(t)u(t) + P (t)u(t) = f (t), u(0) = u 0 , (1.5)
where B(t) and P (t) are bounded operators on H such that Re (B(t) -1 g, g) ≥ δ g 2 for some δ > 0 and all g ∈ H. The left perturbation problem (1.4) was already considered by Arendt et al. [START_REF] Arendt | Maximal regularity for evolution equations governed by non-autonomous forms[END_REF] and the right perturbation one (1.5) by Augner et al. [START_REF] Augner | On the right multiplicative perturbation of non-autonomous L p -maximal regularity[END_REF]. The two problems are motivated by applications to semi-linear evolution equations and boundary value problems. We extend the results in [START_REF] Arendt | Maximal regularity for evolution equations governed by non-autonomous forms[END_REF] and [START_REF] Augner | On the right multiplicative perturbation of non-autonomous L p -maximal regularity[END_REF] in three directions. The first one is to consider general forms which may not satisfy the Kato square root property, a condition which was used in an essential way in the previous two papers. The second direction is to deal with maximal L p -regularity, whereas in the mentioned papers only the maximal L 2 -regularity is considered. The third direction, which is our main motivation, is to assume less regularity on the forms a(t) with respect to t. In both papers [START_REF] Arendt | Maximal regularity for evolution equations governed by non-autonomous forms[END_REF] and [START_REF] Augner | On the right multiplicative perturbation of non-autonomous L p -maximal regularity[END_REF] it is assumed that t → a(t, g, h) is Lipschitz continuous on [0, τ ]. In applications to elliptic operators with time dependent coefficients, the regularity assumption on the forms reflects the regularity needed for coefficients with respect to t.

Our main results can be summarized as follows (see Theorems 3.6 and 5.1 for more general and precise statements). Suppose that for some β, γ ∈ [0, 1],

|a(t, g, h) -a(s, g, h)| ≤ ω(|t -s|) g [H,V] β h [H,V] γ , u, v ∈ V where ω : [0, τ ] → [0, ∞) is a non-decreasing function such that τ 0 ω(t) t 1+ γ 2 dt < ∞.
Suppose also that t → B(t) is continuous on [0, τ ] with values in L(H). Then the Cauchy problem (1.4) has maximal L p -regularity in H for all p ∈ (1, ∞)

when u 0 = 0. If in addition, τ 0 ω(t) p t 1 2 (β+pγ) dt < ∞ (1.6) then (1.4) has maximal L p -regularity in H provided u 0 ∈ (H, D(A(0))) 1-1 p ,p . We also prove that if ω(t) ≤ Ct ε for some ε > 0 and D(A(t) 1/2 ) = V for all t ∈ [0, τ ], then the solution u ∈ C([0, τ ]; V) and s → A(s) 1/2 u(s) ∈ C([0, τ ]; H).
Concerning (1.5), we assume as in [START_REF] Augner | On the right multiplicative perturbation of non-autonomous L p -maximal regularity[END_REF] that t → B(t) is Lipschitz continuous on [0, τ ] with values in L(H). The assumptions on a(t) are the same as above. The maximal L p -regularity results we prove are the same as previously. We could also consider both left and right perturbations, see the end of Section 5.

We point out in passing that condition (1.6) is slightly better than the second condition in (1.3) which was assumed in [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF] and [START_REF] Ouhabaz | Maximal regularity for non-autonomous evolution equations governed by forms having less regularity[END_REF] (for the unperturbed problem). In the natural case ω(t) ∼ t α , one sees immediately that for large p, (1.3) requires larger α (and then more regularity) than (1.6).

In order to prove our results we follow similar ideas as in [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF] and [START_REF] Ouhabaz | Maximal regularity for non-autonomous evolution equations governed by forms having less regularity[END_REF]. However, several modifications are needed in order to deal with multiplicative perturbations. Also, at several places we appeal to classical tools from harmonic analysis such as square function estimates or Hörmander type conditions for singular integral operators with vector-valued kernels.

Our results on maximal L p -regularity could be applied to boundary values problems as well as to some semi-linear evolution equations. Such applications have been already considered in [START_REF] Arendt | Maximal regularity for evolution equations governed by non-autonomous forms[END_REF] and [START_REF] Augner | On the right multiplicative perturbation of non-autonomous L p -maximal regularity[END_REF]. The gain here is that we are able to assume less regularity with respect to the variable t. We shall not write these applications explicitly in this paper since the ideas are the same as in [START_REF] Arendt | Maximal regularity for evolution equations governed by non-autonomous forms[END_REF] and [START_REF] Augner | On the right multiplicative perturbation of non-autonomous L p -maximal regularity[END_REF], one has just to insert our new results on maximal regularity. The reader interested in applications of non-autonomous maximal regularity is referred to the previous articles and the references therein.

The maximal regularity for the unperturbed problem

Let H and V be as in the introduction. We consider a family of sesquilinear forms a(t) :

V × V → C, t ∈ [0, τ ]
which satisfy the classical assumptions [H1]-[H3]. We denote again by A(t) and A(t) the operators associated with a(t) on H and V ′ , respectively. Note that by adding a positive constant to A(t) we may assume that [H3] holds with ν = 0. Therefore, there exists

w 0 ∈ [0, π 2 ) such that a(t, u, u) ∈ Σ(w 0 ), ∀t ∈ [0, τ ], u ∈ V. (2.1) Here Σ(w 0 ) := {z ∈ C * , | arg(z)| ≤ w 0 }.
In (2.1) we take w 0 to be the smallest possible value for which the inclusion holds.

Definition 2.1. Fix u 0 ∈ H. We say that the problem

u ′ (t) + A(t)u(t) = f (t) (t ∈ [0, τ ]), u(0) = u 0 (2.2)
has maximal L p -regularity in H if for each f ∈ L p (0, τ ; H), there exists a unique u ∈ W 1 p (0, τ ; H) such that u(t) ∈ D(A(t)) for almost all t and satisfies (2.2) in the L p -sense.

We denote by V β := [H, V] β the classical complex interpolation space. Its usual norm is denoted • V β . We start with the following result on maximal L p -regularity of (2.2). Theorem 2.2. Suppose that the forms (a(t)) t∈[0,τ ] satisfy the standing hypotheses [H1]- [H3]. Suppose that for some β, γ ∈ [0, 1]

|a(t, u, v) -a(s, u, v)| ≤ ω(|t -s|) u V β v Vγ , u, v ∈ V, (2.3
)

where ω : [0, τ ] → [0, ∞) is a non-decreasing function such that τ 0 w(t) t 1+ γ 2 dt < ∞.
Then The first part of the theorem (i.e., the case u 0 = 0) was proved in [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF] when

β = γ = 1 (and hence [H, V] β = [H, V] γ = V).
The case with different values β and γ was proved in [START_REF] Ouhabaz | Maximal regularity for non-autonomous evolution equations governed by forms having less regularity[END_REF]. See also [START_REF] Arendt | Maximal regularity for non-autonomous Robin boundary conditions[END_REF] for a related result. In order to treat the case of a non-trivial initial data u 0 ∈ (H, D(A(0))) 1-1 p ,p , the assumption required on ω in [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF]

is τ 0 ω(t) t p dt < ∞, (2.5)
and in [START_REF] Ouhabaz | Maximal regularity for non-autonomous evolution equations governed by forms having less regularity[END_REF],

τ 0 ω(t) t β+γ 2 p dt < ∞. (2.6)
In the previous theorem we replace these conditions by the weaker condition (2.4). The important example ω(t) = t α shows that (2.5) and (2.6)) require a large α (and hence more regularity) in the case p > 2, whereas (2.4) does not require any additional regularity than α > γ 2 which is already needed for the first condition

τ 0 w(t) t 1+ γ 2 dt < ∞.
Proof. As explained above the sole novelty here is the treatment of the case u 0 ∈ (H, D(A(0))) 1-1 p ,p under the condition (2.4). Following [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF] and [START_REF] Ouhabaz | Maximal regularity for non-autonomous evolution equations governed by forms having less regularity[END_REF], we have to prove that

t → A(t)e -tA(t) u 0 ∈ L p (0, τ ; H). (2.7)
Since we can assume without loss of generality that A(0) is invertible, then

u 0 ∈ (H, D(A(0))) 1-1 p ,p is equivalent to (see [17, Theorem 1.14]) t → A(0)e -tA(0) u 0 ∈ L p (0, τ ; H). (2.8)
For g ∈ H and a chosen contour Γ in the positive half-plane we write by the holomorphic functional calculus (A(t)e -tA(t) u 0 -A(0)e -tA(0) u 0 , g)

= 1 2πi Γ (ze -tz (zI -A(t)) -1 -(zI -A(0)) -1 u 0 , g) dz = 1 2πi Γ (ze -tz A(0) -A(t) (zI -A(0)) -1 u 0 , (zI -A(0) * ) -1 g)dz = 1 2πi Γ ze -tz a(0, (zI -A(0)) -1 u 0 , (zI -A(0) * ) -1 g)- a(t, (zI -A(0)) -1 u 0 , (zI -A(0)) -1 * g) dz.
Hence by (2.3), the modulus is bounded by

Cω(t) ∞ 0 |z|e -ct|z| (zI -A(0)) -1 u 0 V β (zI -A(t)) -1 * g Vγ d|z|.
Note that by interpolation (see e.g. [START_REF] Ouhabaz | Maximal regularity for non-autonomous evolution equations governed by forms having less regularity[END_REF])

(zI -A(t) * ) -1 L(H,Vγ) ≤ C |z| 1-γ 2 .
(2.9)

On the other hand for

f ∈ D(A(0)), δ (zI -A(0)) -1 f 2 V ≤ Re (A(0)(zI -A(0)) -1 f, (zI -A(0)) -1 f ) ≤ (zI -A(0)) -1 A(0)f (zI -A(0)) -1 f ≤ C |z| A(0)f (zI -A(0)) -1 f V . The embedding V ֒→ V β gives (zI -A(0)) -1 L(D(A(0)),V β ) ≤ C |z| .
Hence, by (2.9) and interpolation

(zI -A(0)) -1 L((H,D(A(0))) 1-1 p ,p ,V β ) ≤ C |z| 1-β 2p .
(2.10)

Using these estimates we obtain

|(A(t)e -tA(t) u 0 -A(0)e -tA(0) u 0 , g)| ≤ Cω(t) ∞ 0 e -ct|z| |z| 1-1 2 (γ+ β p ) d|z| g u 0 (H,D(A(0))) 1-1 p ,p ≤ C ′ ω(t) t 1 2 (γ+ β p ) g u 0 (H,D(A(0))) 1-1 p ,p . Hence, t → A(t)e -tA(t) u 0 ∈ L p (0, τ, H) for u 0 ∈ (H, D(A(0))) 1-1 p ,p if ω(t) satisfies (2.4).

Maximal regularity for left perturbations

This section is devoted to the main subject of this paper in which we are interested in maximal regularity for operators B(t)A(t) for a wide class of operators B(t) and A(t). We will consider in another section the same problem for right multiplicative perturbations A(t)B(t).

Single left multiplicative pertubation-Resolvent estimates

Let H and V be as above. We denote again by • and • V their associated norms, respectively. Let a : V × V → C be a closed, coercive and continuous sesquilinear form. We denote by A and A its associated operators on H and V ′ , respectively. Let b : H × H → C be a bounded sesquilinear form. We assume that b is coercive, that is there exists a constant δ > 0 such that

Re b(u, u) ≥ δ u 2 , u ∈ H. (3.1)
There exists a unique bounded operator associated with b. We denote temporarily this operator by C. Note that by coercivity, it is obvious that C is invertible on H. Now we introduce another operator A b which we call the operator associated with a with respect to b. It is defined as follows

D(A b ) = {u ∈ V, ∃v H : a(u, φ) = b(v, φ) ∀φ ∈ V}, A b u := v.
The difference with A is that we take the form b instead of the scalar product of H in the equality a(u, φ) = b(v, φ). The operator A b is well defined. Indeed, if b(v 1 , φ) = b(v 2 , φ) for all φ ∈ V then by density this equality holds for all φ ∈ H. Therefore, taking φ = v 2 -v 1 and using (3.1), we obtain

v 2 = v 1 . Proposition 3.1. Let B := C -1 . Then A b = BA with domain D(A b ) = D(A). Proof. Let u ∈ D(A b ) and v = A b u. Then a(u, φ) = b(v, φ) = (Cv, φ) ∀φ ∈ V.

Thus, u ∈ D(A) and Au

= Cv = B -1 v. This gives, u ∈ D(A) and A b u = v = BAu.
For the converse, we write for u ∈ D(A) and φ ∈ V

a(u, φ) = (Au, φ) = (CBAu, φ) = b(BAu, φ).

This gives u ∈ D(A b ) and BAu

= A b u.
It is obvious that BA is a closed operator on H. In order to continue we assume that a is coercive (i.e., it satisfies [H3] with ν = 0) and define w 0 and w 1 to be the angles of the numerical ranges of A and B, respectively. That is

(Au, u) ∈ Σ(w 0 ) := {z ∈ C * , | arg(z)| ≤ w 0 } and b(u, u) = (B -1 u, u) ∈ Σ(w 1 )
where w 0 and w 1 are the smallest possible values for which these two properties hold for all u ∈ V. Note that w 0 , w 1 ∈ [0, π 2 ) because of the coercivity property.

Proposition 3.2. For all λ /

∈ Σ(w 0 + w 1 ), the operator λI-BA is invertible on H and

(λI -BA) -1 L(H) ≤ δ -1 B -1 L(H) dist(λ, Σ(w 0 + w 1 )) . Proof. Let u ∈ D(A). We write (λI -BA)u u = B(λB -1 -A)u u ≥ 1 B -1 L(H) (λB -1 I -A)u u ≥ 1 B -1 L(H) |(λB -1 u -Au, u)| = |(B -1 u, u)| B -1 L(H) |λ - (Au, u) (B -1 u, u) |. Since (Au,u) (B -1 u,u) = a(u,u) b(u,u) ∈ Σ(w 0 + w 1 ) it follows that |(λ - (Au, u) (B -1 u, u) | ≥ dist(λ, Σ(w 0 + w 1 )).
On the other hand, by (3.1), |(B -1 u, u)| ≥ δ u 2 and so

(λI -BA)u u ≥ δ B -1 L(H) u 2 dist(λ, Σ(w 0 + w 1 )).
Hence,

(λI -BA)u ≥ δ B -1 L(H) u dist(λ, Σ(w 0 + w 1 )) ∀u ∈ D(A). (3.2)
This implies that λI-BA is injective and has closed range for λ ∈ Σ(w 0 + w 1 ).

In order to prove that λI -BA is invertible it remains to prove that it has dense range. By duality, one has to prove that the adjoint is injective.The adjoint operator is λI -A * B * . We write

λI -A * B * = (λB * -1 -A * )B * .
The previous arguments show that λB -1 -A is injective. This also applies to λB * -1 -A * . Since B * is invertible, we obtain λI -A * B * is injective and hence λI -BA is invertible . Now (3.2) gives

(λI -BA) -1 ≤ B -1 L(H) δ.dist(λ, Σ(w 0 + w 1 ))
for all λ ∈ Σ(w 0 + w 1 )).

Corollary 3.3. Suppose that w 0 + w 1 < π 2 .Then -BA is the generator of a bounded holomorphic semigroup on H.

Proof. By Proposition 3.2, (λI -BA) -1 ≤ c |λ| , ∀λ ∈ Σ(w 0 + w 1 ))
In other words, λI + BA is invertible for λ ∈ Σ(π -(w 0 + w 1 )) and :

(λI + BA) -1 ≤ c |λ| , ∀λ ∈ Σ(π -(w 0 + w 1 )).
It is a classical fact that the latter estimate implies that -BA generates a bounded holomorphic semigroup of angle π 2 -(w 0 + w 1 ). Obviously, one cannot remove the assumption w 0 + w 1 < π 2 in the previous result. Indeed, let A = -e i π 3 ∆ on L 2 (R d ) and B be the multiplication by e i π 3 . Then -BA = e i 2π 3 ∆ is not a generator of a C 0 -semigroup.

Single pertubation-Maximal regularity

Let (a(t)) t∈[0,τ ] , A(t), A(t) and b be as in the previous sub-section. We assume that [H3] holds with ν = 0. In particular, (2.1) holds. We also have

b(u, u) ∈ Σ(w 1 ) (3.3)
for some w 1 ∈ [0, π 2 ) by coercivity of b. We make the assumption w 0 + w 1 < π 2 . By Corollary 3.3, for each t ∈ [0, τ ], the operator -BA(t) generates a holomorphic semigroup (e -sBA(t) ) s≥0 on H.

Our aim in this section is to prove maximal regularity in H for the Cauchy problem associated with BA(t), t ∈ [0, τ ]. The definition of maximal L p -regularity in this context is the same as in Definition 2.1. Set R(λ, BA(t)) := (λI + BA(t)) -1 for λ ∈ ρ(-BA(t)).

Proposition 3.4. Assume that w

0 + w 1 < π 2 . Then 1-(λB -1 + A(t)) -1 L(H) ≤ C |λ|+1 , λ ∈ Σ(π -(w 0 + w 1 )), 2-R(λ, BA(t))B L(V ′ ,H) ≤ C (|λ|+1) 1 2 , λ ∈ Σ(π -(w 0 + w 1 )), 3-e -(t-s)BA(t) B L(V ′ ,H) ≤ C (t-s) 1 2 , 4-e -(t-s)BA(t) B L(V ′ ,V) ≤ C (t-s)
. The constant C is independent of t and λ.

Proof. We have (λB -1 +A(t)) -1 = (λ+BA(t)) -1 B, then we obtain assertion 1-from Proposition 3.2. Note that

(λB -1 +A(t)) -1 = (λ+A(t)) -1 +(λB -1 +A(t)) -1 (λ(-B -1 +I))(λ+A(t)) -1 . (3.4) Then R(λ, BA(t))B L(V ′ ,H) = (λB -1 + A(t)) -1 L(V ′ ,H) ≤ (λI + A(t)) -1 L(V ′ ,H) + (λB -1 + A(t)) -1 (λ(-B -1 + I)) L(H) (λ + A(t)) -1 L(V ′ ,H) . Since (λI + A(t)) -1 L(V ′ ,H) ≤ C (|λ| + 1) 1 2
(see e.g. [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF]), we obtain

R(λ, BA(t))B L(V ′ ,H) ≤ C (|λ| + 1) 1 2
, which proves assertion 2. Now we choose an appropriate contour Γ = ∂Σ(θ) with θ < π 2 and write by the functional calculus

e -(t-s)BA(t) B = 1 2πi Γ e -(t-s)λ (λ -BA(t)) -1 Bdλ. Then e -(t-s)BA(t) B L(V ′ ,H) ≤ 1 2π ∞ 0 e -(t-s)Re λ (λ -BA(t)) -1 B L(V ′ ,H) d|λ| ≤ C ∞ 0 e -(t-s)Re λ 1 (|λ| + 1) 1 2 d|λ| ≤ C ′ (t -s) 1 2
.

In order to prove assertion 4-we write

e -(t-s)BA(t) B L(V ′ ,V) ≤ e -(t-s) 2 BA(t) BB -1 L(H,V) e -(t-s) 2 BA(t) B L(V ′ ,H) and e -(t-s) 2 BA(t) BB -1 L(H,V) ≤ B -1 L(H) e -(t-s) 2 BA(t) B L(H,V) .
We use the equality

(λB -1 + A(t)) -1 = (λI + A(t)) -1 + (λI + A(t)) -1 λ(I -B -1 )(λB -1 + A(t)) -1
in place of (3.4) to estimate R(λ, BA(t))B L(H,V) and then argue as previously.

Now, let P (t) ∈ L(H) such that t → P (t) is strongly measurable and

P (t) L(H) ≤ M, t ∈ [0, τ ] (3.5)
for some constant M . We consider the Cauchy problem

u ′ (t) + BA(t)u(t) + P (t)u(t) = f (t), u(0) = u 0 . ( 3.6)
Recall that B -1 is the operator associated with b. We are interested in maximal regularity of (3.6). As explained at the beginning of the proof of the next proposition, we may assume without loss of generality that the forms a(t) are coercive and hence (2.1) is satisfied for some w 0 ∈ [0, π 2 ). Proposition 3.5. Suppose that the forms (a(t)) t∈[0,τ ] satisfy [H1]-[H3], the form b satisfies (3.1) and w 0 + w 1 < π 2 . Suppose that for some β, γ

∈ [0, 1] |a(t, u, v) -a(s, u, v)| ≤ ω(|t -s|) u V β v Vγ , u, v ∈ V where ω : [0, τ ] → [0, ∞) is a non-decreasing function such that : τ 0 w(t) t 1+ γ 2 dt < ∞.
Then the Cauchy problem (3.6)

with u 0 = 0 has maximal L p -regularity in H for all p ∈ (1, ∞). If in addition, τ 0 ω(t) p t 1 2 (β+pγ) dt < ∞ (3.7)
then (3.6) has maximal L p -regularity for all u 0 ∈ (H, D(A(0))) 1-1 p ,p . Moreover there exists a positive constant C such that :

u W 1 p (0,τ ;H) + Au Lp(0,τ ;H) ≤ C f Lp(0,τ ;H) + u 0 (H,D(A(0))) 1-1 p ,p
.

Here C depends only on the constants in

[H1]-[H3], B L(H) , B -1 L(H)
and M in (3.5).

Proof. Firstly, we note that for c ∈ R, (3.6) has maximal L p -regularity if and only if the Cauchy problem

v ′ (t) + (BA(t) + P (t) + cI)v(t) = e -ct f (t), v(0) = u 0
has maximal L p -regularity. The reason is that v(t) = u(t)e -ct and it is clear that u ∈ W 1 p (0, τ ; H) if and only if v ∈ W 1 p (0, τ ; H). Thus, by adding a large constant c we may assume that [H3] holds with ν = 0 and BA(t) + P (t) is invertible for each t ∈ [0, τ ]. Note that BA(t) = A(t) b is the operator associated with the form a(t) with respect to b (see Section 3.1). This allows us to use the same strategy of proof as for Theorem 2.2 (cf. [START_REF] Ouhabaz | Maximal regularity for non-autonomous evolution equations governed by forms having less regularity[END_REF] or [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF] in the case β = γ = 1).

Set v(s) := e -(t-s)BA(t) u(s). Writing v(t) -v(0) = t 0 v ′ (s)ds we obtain

A(t)u(t) = A(t)e -tBA(t) u 0 + A(t) t 0 e -(t-l)BA(t) B(A(t) -A(l))u(l)dl + A(t) t 0 e -(t-l)BA(t) (-P (l))u(l)dl + A(t) t 0 e -(t-l)BA(t) f (l)dl.
Note that by Proposition 3.4, the term e -(t-l)BA(t) B(A(t)-A(l))u(l) is well defined.

We first prove the proposition in the case u 0 = 0. We define (Lf )(t) := A(t) t 0 e -(t-l)BA(t) f (l)dl.

Following [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF] the operator L is a pseudo-differential operator with the vector-valued symbol σ(t, ξ) given by

σ(t, ξ) :=      A(0)(iξ + B(0)A(0)) -1 if t < 0 A(t)(iξ + B(t)A(t)) -1 if 0 ≤ t ≤ τ A(τ )(iξ + B(τ )A(τ )) -1 if t > τ.
Then we use Proposition 3.2 and argue as in the proof of Lemmas 10 and 11 in [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF] to prove the boundedness on L p (0, τ ; H), 1 < p < ∞, of the operator L. We continue as in [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF] and [START_REF] Ouhabaz | Maximal regularity for non-autonomous evolution equations governed by forms having less regularity[END_REF]. We set (Sg)(t) := A(t) t 0 e -(t-l)BA(t) (P (l))A(l) -1 g(l)dl.

By the boundedness of the operator L on L p (0, τ ; H),

Sg Lp(0,τ ;H) ≤ C A -1 g Lp(0,τ ;H) . We define (Qg)(t) := A(t) t 0 e -(t-l)BA(t) B(A(t) -A(l))A(l) -1 g(l)dl.
Then, arguing as in [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF] or [START_REF] Ouhabaz | Maximal regularity for non-autonomous evolution equations governed by forms having less regularity[END_REF] we obtain easily from Proposition 3.4

(Qg)(t) ≤ t 0 w(|t -l|) (t -l) 1+ γ 2 A -1 (l)g(l) V dl.
Thus,

Qg Lp(0,τ ;H) ≤ C τ 0 w(t) t 1+ γ 2 dt A -1 g Lp(0,τ ;V) .
From these estimates, we see that by replacing A(t) by A(t) + cI for c large enough we obtain

S L(Lp(0,τ ;H)) < 1 4 and Q L(Lp(0,τ ;H)) < 1 4 .
In particular, I -

(S + Q) is invertible. Since (Au)(t) = (I -(S + Q)) -1 (L(f ))(t)
we obtain Au ∈ L p (0, τ ; H) and hence u ∈ W 1 p (0, τ ; H). This proves maximal L p -regularity.

In order to treat the case u 0 = 0 we need to estimate the difference of the resolvents, i.e., R(λ, A(t) b ) -R(λ, A(s) b ) L(H) in terms of ω(|t -s|). Let f, g ∈ H and λ ∈ Σ(π -(w 0 + w 1 )). We write

([R(λ, A(t) b ) -R(λ, A(s) b )]f, g) = -([R(λ, BA(t))B(A(t) -A(s))R(λ, BA(t))]f, g).
Note that the RHS is well defined since R(λ, BA(t))B is a bounded operator from V ′ to V (cf. Proposition 3.4). Therefore,

([R(λ, A(t) b ) -R(λ, A(s) b )]f, g) = (A(t) -A(s))R(λ, BA(t))f, B * R(λ, BA(t)) * g = a(s, R(λ, BA(s))f, (λB * -1 + A(t) * ) -1 g) -a(t, R(λ, BA(s))f, (λB * -1 + A(t) * ) -1 g).
Hence the modulus is bounded by

ω(|t -s|) R(λ, BA(s))f V β λB * -1 + A(t) * ) -1 g Vγ .
Let w 0 be the common angle for the numerical range of a(t). By Proposition 3.2 we have for all λ / ∈ Σ(w 0 + w 1 )

δ R(λ, A(s) b )f 2 V ≤ Re a(s, R(λ, A(s) b )f, R(λ, A(s) b )f ) = Re (A(s)R(λ, A(s) b )f, R(λ, A(s) b )f ) = Re (BA(s)R(λ, A(s) b )f, (B -1 ) * R(λ, A(s) b )f ) ≤ C |λ| f 2 .
Hence, by interpolation

R(λ, A(s) b )f 2 V β ≤ C |λ| 1-β 2 . (3.8)
Putting together the previous estimates yields

|b([R(λ, BA(t)) -R(λ, BA(s))]f, g)| ≤ C ω(|t -s|) |λ| 2-β+γ 2 f g . This shows R(λ, A(t) b ) -R(λ, A(s) b ) L(H) ≤ C ω(|t -s|) |λ| 2-β+γ 2 .
This is the estimate we need in order to obtain the proposition when u 0 ∈ (H, D(A(0)) 1-1 p ,p (see [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF] or [START_REF] Ouhabaz | Maximal regularity for non-autonomous evolution equations governed by forms having less regularity[END_REF] for the details).

Time dependent perturbations-Maximal regularity

Let a(t), A(t), V and H be as above and suppose again that the standard assumptions [H1]-[H3] are satisfied. Let (B(t)) t∈[0,τ ] be a family of bounded invertible operators on H. We assume that there exist constants δ > 0 and M > 0 independent of t such that Re (B(t

) -1 u, u) ≥ δ u 2 H ∀u ∈ H, (3.9)
and

B(t) -1 L(H) ≤ M. (3.10)
Let (P (t)) t∈[0,τ ] be a family of bounded operators on H. We assume that

P (t) L(H) ≤ M. (3.11)
As a consequence of (3.9) and (3.10) the numerical range of B(t) -1 is contained in a sector of angle w 1 for some w 1 ∈ [0, π 2 ), independent of t. Note that (3.9) implies that

B(t) -1 u ≥ δ u and hence B(t) L(H) ≤ 1 δ . (3.12)
We denote as previously by w 0 the common angle of the numerical range of forms a(t), t ∈ [0, τ ]. We assume again that

w 0 + w 1 < π 2 . (3.13)
The following is our main result.

|a(t, u, v) -a(s, u, v)| ≤ ω(|t -s|) u V β v Vγ , u, v ∈ V (3.14)
where ω : [0, τ ] → [0, ∞) is a non-decreasing function such that :

τ 0 ω(t) t 1+ γ 2 dt < ∞. ( 3
then u ′ (t) + B(t)A(t)u(t) + P (t)u(t) = f (t), u(0) = u 0 (3.18) has maximal L p -regularity in H provided u 0 ∈ (H, D(A(0))) 1-1 p ,p .
Moreover there exists a positive constant C such that :

u W 1 p (0,τ ;H) + BAu Lp(0,τ ;H) ≤ C f Lp(0,τ ;H) + u 0 (H,D(A(0))) 1-1 p ,p . ( 3

.19) The constant C depends only on the constants in [H1]-[H3]

, δ and M in (3.9)- (3.11).

Remark. As we shall see in the proof, the regularity assumption on B(t) can be weakened considerably. Indeed, continuity at finite number of appropriate points is sufficient. uniqueness of the fixed point of S. The a priori estimate of Proposition 3.5 and the fact that v = u on [0, t 0 ] give

u M R(p,H) ≤ C f Lp(0,t 0 ;H) + u 0 (H,D(A(0))) 1-1 p ,p .
(3.21)

Next, we divide the interval [0, τ ] into ∪ i∈1,....,N [t i-1 , t i ] with t i -t i-1 small enough. On each interval [t i-1 , t i ], we search for a solution u i to (3.16) with initial data u i (t i-1 ) = u i-1 (t i-1 ). The forgoing arguments prove existence and uniqueness of a solution on

[t i-1 , t i ] with maximal L p - regularity provided u i-1 (t i-1 ) ∈ (H, D(A(t i-1 ))) 1-1 p ,p .
Once we do this we glue these solutions and obtain a unique solution u of (3.16) with maximal L p -regularity on [0, τ ] for all u 0 ∈ (H, D(A(0))) 1-1 p ,p . Thus, our task now is to prove that u i-1 (t i-1 ) ∈ ( H, D(A(t i-1 ))) 1-1 p ,p . In order to make the notation simpler, we work on [0, τ ] (with τ small enough) instead of [t i-1 , t i ] and set A := A(τ ). We have to prove that the solution u to (3.16) satisfies u(τ ) ∈ (H, D(A)) 1-1 p ,p . This means that (remember we always assume w.l.o.g. that the operators A(t) are invertible, i.e. ν = 0 in [H3])

t → Ae -tA u(τ ) ∈ L p (0, τ ; H). (3.22)
We start with an expression for u(τ ). Set

v(s) := e -(τ -s)A u(s), 0 ≤ s ≤ τ.
We have v ′ (s) = Ae -(τ -s)A u(s) + e -(τ -s)A (-B(s)A(s)u(s) -P (s)u(s) + f (s)).

Hence

u(τ ) = e -τ A u 0 + τ 0 e -(τ -s)A (A(τ ) -A(s))u(s)ds (3.23) + τ 0 e -(τ -s)A [(I -B(s)A(s)u(s) -P (s)u(s) + f (s)] ds. Since Ae -(t+τ )A L(H) ≤ C t + τ it is clear that t → Ae -tA e -τ A u 0 ∈ L p (0, τ ; H).
For the second term we have

A τ 0 e -(t+τ -s)A (A(τ ) -A(s))u(s)ds ≤ τ 0 C t + τ -s e -1 2 (t+τ -s)A (A(τ ) -A(s))u(s)ds ≤ C τ 0 ω(τ -s) (t + τ -s) 1+ γ 2 u(s) V β ds ≤ C τ 0 ω(t + τ -s) (t + τ -s) 1+ γ 2 u(s) V β ds,
where ω(r) = ω(r) for r ∈ [0, τ ] and = ω(τ ) for r > τ . Therefore, using the assumption (3.15) on ω and Young's inequality we obtain

τ 0 A τ 0 e -(t+τ -s)A (A(τ ) -A(s))u(s)ds p dt ≤ C ′ u p L p (0,τ ;V) . (3.24)
Now we consider the last term in (3.23). We start with the case p = 2. Set g(s

) := (I -B(s)A(s)u(s) -P (s)u(s) + f (s). (3.25)
We have

τ 0 Ae -tA τ 0 e -(τ -s)A g(s)ds 2 dt 1/2 = τ 0 A 1/2 e -tA τ 0 A 1/2 e -(τ -s)A g(s)ds 2 dt 1/2 ≤ C τ 0 A 1/2 e -(τ -s)A g(s)ds .
In the last inequality we use the boundedness of the square function, namely

∞ 0 A(t) 1/2 e -rA(t) x 2 dr ≤ C x 2 (3.26)
for all x ∈ H. This estimate is a consequence of the fact that A(t) has a bounded holomorphic functional calculus as an accretive operator, see [START_REF] Cowling | Banach space operators with a bounded H ∞ functional calculus[END_REF]. We repeat the previous argument but since g is not necessarily constant in s we cannot use directly the square function estimate. We argue by duality. For x ∈ H we have

|( τ 0 A 1/2 e -(τ -s)A g(s)ds, x)| = | τ 0 (g(s), A * 1/2 e -(τ -s)A * x)| ≤ ( τ 0 g(s) 2 ds) 1/2 ( τ 0 A * 1/2 e -(τ -s)A * x 2 ds) 1/2 ≤ C x ( τ 0 g(s) 2 ds) 1/2 .
Since this is true for all x ∈ H, we obtain

τ 0 Ae -tA τ 0 e -(τ -s)A g(s)ds 2 dt 1/2 ≤ C( τ 0 g(s) 2 ds) 1/2 .
We define the operator T by

T g(t) = τ 0
Ae -(τ +t-s)A g(s)ds.

We have proved that T : L 2 (0, τ ; H) → L 2 (0, τ ; H) is bounded. We extend this operator to L p (0, τ ; H) for all p ∈ (1, ∞). Indeed, note that T is a singular integral operator with kernel

K(t, s) = Ae -(τ +t-s)A
and we use Hörmander's integral condition for K(t, s) and K(s, t) (see, e.g. [START_REF] Rubio De Francia | Calderón-Zygmund theory for operator-valued kernels[END_REF], Theorems III 1.2 and III 1.3). A similar argument was used in [START_REF] Haak | Maximal regularity for non-autonomous evolution equations[END_REF].

We have to prove that

|t-s|≥2|s ′ -s| K(t, s) -K(t, s ′ ) L(H) dt ≤ C (3.27)
for some constant C independent of s, s ′ ∈ (0, τ ).

Assume for example that s ≤ s ′ . Since the semigroup generated by -A is bounded holomorphic we have for some constant

C |t-s|≥2|s ′ -s| K(t, s) -K(t, s ′ ) L(H) dt = τ 2s ′ -s Ae -(τ +t-s)A -Ae -(τ +t-s ′ )A L(H) = τ 2s ′ -s s ′ s A 2 e -(τ +t-r)A dr L(H) dt ≤ C τ 2s ′ -s s ′ s 1 (τ + t -r) 2 drdt = C τ 2s ′ -s 1 τ + t -s ′ - 1 τ + t -s dt = C log τ + t -s τ + t -s ′ t=τ t=2s ′ -s ≤ C log 2.
This proves (3.27). The same arguments apply for the kernel of the adjoint T * . Hence T : L p (0, τ ; H) → L p (0, τ ; H) is a bounded operator for all p ∈ (1, ∞). We obtain

τ 0 Ae -tA τ 0 e -(τ -s)A g(s)ds p dt 1/p ≤ C( τ 0 g(s) p ds) 1/p . (3.28)
Note that for f ∈ L p (0, τ ; H) we have A(•)u(•) ∈ L p (0, τ ; H) since we have proved maximal L p -regularity for small τ . Hence g ∈ L p (0, τ ; H) (remember that g is given by (3.25)). This finishes the proof of u(τ ) ∈ (H, D(A(τ ))

u(t i ) (H;D(A(t i ))) 1-1 p ;p ≤ C f Lp(0,t i ;H) + u 0 (H,D(A(0))) 1-1 p ,p + u Lp(0,t i ;V) .
(3.29) Remember that we can replace A(t) by A(t) + cI for a given constant c ∈ (0, +∞). Since (A(t) + cI) -1 L(H,V) ≤ c 0 (c + 1)

1 2 it follows that u(t) V ≤ (A(t) + cI) -1 L(H,V) (A(t) + cI)u(t) H ≤ c 0 (c + 1) 1 2 (A(t) + cI)u(t) H .
Summing over i in (3.29) and taking c large enough we see that for some constant

C 1 u W 1 p (0,τ ;H) + A(•)u(•) Lp (0,τ ;H) ≤ C 1 f Lp(0,τ ;H) + u 0 (H,D(A(0))) 1-1 p ,p .
This proves the desired a priori estimate and finishes the proof of the theorem.

Further regularity results

We continue our investigations on the solution of the problem (3.18). We work with the same assumptions as in Theorem 3.6. For f ∈ L 2 (0, τ ; H) the solution u ∈ W 1 2 (0, τ ; H) and by the Sobolev embedding u ∈ C([0, τ ]; H). It is interesting to know whether u is also continuous for the norm of V. This is indeed the case if the forms a(t) are symmetric (or perturbations of symmetric forms) and t → a(t, x, y) is Lipschitz continuous on [0, τ ] for all x, y ∈ V. This is proved in [START_REF] Arendt | Maximal regularity for evolution equations governed by non-autonomous forms[END_REF]. Continuity in V was also proved in [START_REF] Arendt | Maximal regularity for non-autonomous Robin boundary conditions[END_REF] for the unperturbed problem (i.e., without multiplicative and additive perturbations) when γ in (3.14) is < 1. This is a rather restrictive condition but turns out to be satisfied in some cases such as time-dependent Robin boundary conditions. Here we make no restriction on γ and β and we assume less regularity for t → a(t, x, y) than what was previously known. The continuity of the solution with respect to the norm of V is used in [START_REF] Arendt | Maximal regularity for evolution equations governed by non-autonomous forms[END_REF] in applications to some semi-linear PDE's. In this section, we look at again this question in the setting of Theorem 3.6 in which we assume less regularity (than Lipschitz continuous) on t → a(t, g, h).

In the statements below we shall need the following square root property (called Kato's square root property)

D(A(t) 1/2 ) = V and c 1 A(t) 1/2 v ≤ v V ≤ c 2 A(t) 1/2 v (4.1)
for all v ∈ V and t ∈ [0, τ ], where the positive constants c 1 and c 2 are independent of t. Note that this assumption is always true for symmetric forms when ν = 0 in [H3].

We start with the following lemma which will be used later.

Lemma 4.1. Suppose (4.1). Then for all f ∈ L 2 (0, t;

H), 0 ≤ s ≤ t ≤ τ , t s e -(t-r)A(t) f (r)dr V ≤ C f L 2 (s,t;H) .
Proof. By (4.1),

t s e -(t-r)A(t) f (r)dr V ≤ c 2 t 0 A(t) 1/2 e -(t-r)A(t) f (r)dr = c 2 sup x =1 | t s (f (r), A(t) * 1/2 e -(t-r)A(t) * x)dr| ≤ c 2 sup x =1 t s A(t) * 1/2 e -(t-r)A(t) * x 2 dr 1 2 × f L 2 (s,t;H) ≤ C f L 2 (s,t;H) .
Note that in the last inequality we use again the square function estimate for A(t) * (see (3.26)). This proves the lemma.

In the next result we prove continuity of the solution to (3.18) as a function with values in V. Note that if D(A(0)

1/2 ) = V, then (D(A(0), H) 1 2 ,2 = D(A(0) 1/2 ) = V.
Theorem 4.2. Suppose (4.1) and that the assumptions of Theorem 3.6 are satisfied. Suppose also that ω(t) ≤ ct ε for some ε > 0. Let f ∈ L 2 (0, τ ; H) and u 0 ∈ V. Then the solution u to the problem (3.18) satisfies u ∈ C([0, τ ]; V).

Proof. We fix s and t in [0, τ ] such that s < t. We first derive a formula similar to (3.23). Define v(r) := e -(t-r)A(t) u(r) for r ∈ [s, t]. After derivation and integration from s to t we obtain

u(t) = e -(t-s)A(t) u(s) + t s e -(t-r)A(t) (A(t) -B(r)A(r))u(r)dr + t s e -(t-r)A(t) [-P (r)u(r) + f (r)]dr.
Hence

u(t) = e -(t-s)A(t) u(s) + t s e -(t-r)A(t) (A(t) -A(r))u(r)dr + t s e -(t-r)A(t) [(-B(r) + I)A(r)u(r) -P (r)u(r) + f (r)]dr, (4.2)
and so

u(t) -u(s) V ≤ t s e -(t-r)A(t) (A(t) -A(r))u(r)dr V + t s e -(t-r)A(t) [(I -B(r))A(r)u(r) -P (r)u(r) + f (r)]dr V + [e -(t-s)A(t) -I]u(s) V .
We estimate each term in the RHS and obtain

u(t) -u(s) V ≤ C( t s ω(t -r) t -r u(r) V dr + f L 2 (s,t;H) + Au L 2 (s,t;H) + u L 2 (s,t;H) ) + [e -(t-s)A(t) -I]u(s) V ≤ C( t-s 0 ω(r) r dr u L ∞ (0,τ ;V) + f L 2 (s,t;H) + Au L 2 (s,t;H) + u L 2 (s,t;H) ) + [e -(t-s)A(t) -I]u(s) V .
Note that M R(2, H) is continuously embedded into L ∞ (0, τ ; V) by Proposition 4.5 in [START_REF] Monniaux | The incompressible Navier-Stokes system with time-dependent Robin-type boundary conditions[END_REF]. This proposition is proved for forms which are symmetric but it remains true under the assumption (4.1). Thus, by maximal regularity result (Theorem 3.6), u ∈ M R(2, H) and hence u ∈ L ∞ (0, τ ; V). Next, since ω satisfies the assumptions of Theorem 3.6, τ 0 ω(r) r dr < ∞ and since u, Au ∈ L 2 (0, τ ; H) we see that the first four terms in the RHS converge to 0 as t → s (or as s → t). It remains to proves that [e -(t-s)A(t) -I]u(s) V → 0 as t → s (or as s → t).

(4.3)

We first prove (4.3) when t → s (for fixed s). We write

[e -(t-s)A(t) -I]u(s) V ≤ [e -(t-s)A(t) -e -(t-s)A(s) ]u(s) V + [e -(t-s)A(s) -I]u(s) V .
Using Proposition 3.4 and the functional calculus (on the sector Σ(θ) for appropriate θ ∈ (0, π 2 )) we estimate the first terms as follows.

[e -(t-s)A(t) -e -(t-s)A(s) ]u(s)

V = 1 2π Γ e -(t-s)λ [(λ -A(t)) -1 (A(t) -A(s))(λ -A(s)) -1 ]u(s)dλ V ≤ cω(t -s) ∞ 0 e -(t-s)|λ| cos θ (1 + |λ|) -1 d|λ| u(s) V ≤ c ′ ω(t -s) (t -s) ε ′ u(s) V .
Here we use u(s) ∈ D(A(s) 12 ) (see the proof of Theorem 3.6) and (4.1). Now the fact that ω(t) ≤ ct ε for some ε > 0 and the strong continuity of the semigroup e -tA(s) on V imply that [e -(t-s)A(t) -I]u(s) V → 0 as t → s. This proves that u is right continuous for the norm of V.

It remains to prove left continuity of u. We need a formula similar to (4.2) but with u(s) expressed in terms of u(t). Fix 0 ≤ s < t ≤ τ and set v(r) := e -(r-s)A(s) u(r) for r ∈ [s, t]. Then v ′ (r) = -e -(r-s)A(s) (A(s) + B(r)A(r) + P (r))u(r)) + e -(r-s)A(s) f (r), and hence

u(s) = e -(t-s)A(s) u(t) + t s e -(r-s)A(s) (A(s) + B(r)A(r))u(r)dr - t s e -(r-s)A(s) [f (r) -P (r)u(r)]dr. (4.4) Therefore u(s) -u(t) = [e -(t-s)A(s) u(t) -u(t)] + t s e -(r-s)A(s) (A(s) -A(r))u(r)dr + t s e -(r-s)A(s) ((B(r) + I)A(r)u(r))dr - t s e -(r-s)A(s) [f (r) -P (r)u(r)]dr =: I 1 (s, t) + I 2 (s, t) + I 3 (s, t) + I 4 (s, t).
By Lemma 4.1,

I 4 (s, t) V ≤ C[ u L 2 (s,t;H) + f L 2 (s,t;H) ].
By Lemma 4.1

I 3 (s, t) V ≤ C A(•)u(•) L 2 (s,t;H) . For I 2 (s, t) we have immediately, I 2 (s, t) V ≤ C t s w(r -s) r -s dr u L ∞ (0;τ,V) .
For I 1 (s, t) we proceed as before. We write [e -(t-s)A(s) u(t) -u(t)] = [e -(t-s)A(s) u(t) -e -(t-s)A(t) u(t)]

+ [e -(t-s)A(t) u(t) -u(t)]

We use again the functional calculus as above to obtain [e -(t-s)A(s) u(t) -e -(t-s)A(t) u(t)] V ≤ c ω(t -s) (t -s) ε ′ u(t) V .

The remaining term e -(t-s)A(t) u(t) -u(t) V converges to 0 as s → t by strong continuity of the semigroup on V. We have proved that u is left continuous in V and finally u ∈ C([0, τ ]; V). t s e -(t-r)A(t) (A(t) -A(r))u(r)dr

+ A(t) 1 2
t s e -(t-r)A(t) [(-B(r) + I)A(r)u(r) -P (r)u(r) + f (r)]dr. By (4.1), the norms in H of the last two terms are equivalent to the norms in V of the same terms but without A(t) 1 2 . We have seen in the proof of Theorem 4.2 that these norms in V converge to 0 as t → s or as s → t.

It remains to consider the term A(t)

1 2 e -(t-s)A(t) u(s) -A(s) 1 2 u(s). We use again the functional calculus to write

A(t)

1 2 e -(t-s)A(t) u(s) -A(s)

1 2 u(s) = A(t)
≤ (λ -A(t)) -1 L(V ′ ,H) (A(s) -A(t)) L(V,V ′ ) (λ -A(s)) -1 L(V) ≤ C|λ| -1/2 ω(|t -s|) 1 1 + |λ| .
Therefore for any ε > 0,

A(t)
1 2 e -(t-s)A(t) u(s) -A(s) Remember that u ∈ L ∞ (0, τ ; V) by Theorem 4.2. Using the assumption on ω, the latest term converges to 0 as |t-s| → 0. The term A(s)

1 2 e -(t-s)A(s) u(s)-A(s) 1 2 u(s) converges to 0 as t → s by the strong continuity of the semigroup. This proves the right continuity of s → A(s) 1 2 u(s). The left continuity is proved similarly, we use (4.4) instead of (4.2).

Right perturbations-Maximal regularity

Let B(t) and P (t) (t ∈ [0, τ ]) be bounded operators on H. We investigate the maximal L p -regularity property for right multiplicative perturbations A(t)B(t). As mentioned in the introduction, this problem has been considered in [START_REF] Augner | On the right multiplicative perturbation of non-autonomous L p -maximal regularity[END_REF] and was motivated there by several applications. We will extend the results from [START_REF] Augner | On the right multiplicative perturbation of non-autonomous L p -maximal regularity[END_REF] in the sense that we require much less regularity for t → a(t).

Let a(t) be a family of sesquilinear forms satisfying again [H1]-[H3] and denote as before A(t) the corresponding associate operators. Under the assumptions of Theorem 3.6, for each t, the operator -B(t)A(t) generates a holomorphic semigroup on H. The same is also true for -B(t) * A(t) * since the adjoint operators B(t) * and A(t) * satisfy the same properties as B(t) and A(t). Hence by duality, -A(t)B(t) generates a holomorphic semigroup on H. Here, the domain of A(t)B(t) is given by

D(A(t)B(t)) = {x ∈ H, B(t)x ∈ D(A(t))}.

For right perturbations, we say that the Cauchy problem u ′ (t) + A(t)B(t)u(t) + P (t)u(t) = f (t), u(0) = u 0

Notation.

  We denote by L(E, F ) (or L(E)) the space of bounded linear operators from E to F (from E to E). The spaces L p (a, b; E) and W 1 p (a, b; E) denote respectively the Lebesgue and Sobolev spaces of function on (a, b) with values in E. Recall that the norms of H and V are denoted by • and • V . The scalar product of H is (•, •). Finally, we denote by C, C ′ or c... all inessential constants. Their values may change from line to line.

Proposition 4 . 3 . 1 2

 431 Suppose that the assumptions of the previous theorem are satisfied. Let f ∈ L 2 (0, τ ; H), u 0 ∈ V and u be the solution of(3.18). Thens → A(s) u(s) ∈ C([0, τ ]; H).Proof. We use again (4.2) and write

1 2 e

 2 -(t-s)A(t) u(s) -A(s)

1 2 e

 2 -(t-s)A(t) u(s) -A(s)

1 2 e

 2 -(t-s)A(s) u(s) + A(s)

1 2 e 1 2

 21 -(t-s)A(s) u(s) -A(s) e -(t-s)λ [(λ -A(t)) -1 -(λ -A(s)) -1 ]u(s)dλ + A(s)

1 2 e

 2 -(t-s)A(s) u(s) -A(s) 1 2 u(s).

1 and

 1 By the resolvent equation(λ -A(t)) -1 -(λ -A(s)) -1 = (λ -A(t)) -1 (A(s) -A(t))(λ -A(s)) -Proposition 3.4 we have (λ -A(t)) -1 -(λ -A(s)) -1

1 2 e

 2 -(t-s)A(s) u(s) ≤ Cω(|t -s|) e -(t-s)r dr u(s) V ≤ C ε ω(|t -s|) |t -s| ε ′ u(s) V .

Research partially supported by the ANR project HAB: ANR-12-BS01-0013-02.

Before starting the proof, let us define the maximal regularity space M R(p, H) := {u ∈ W 1 p (0, τ ; H) : u(t) ∈ D(A(t)) a.e., A(.)u(.) ∈ L p (0, τ ; H)}.

It is a Banach space for the norm u M R(p,H) := u W 1 p (0,τ ;H) + Au Lp(0,τ ;H) .

Proof. Let f ∈ L p (0, τ ; H) and u 0 ∈ (H, D(A(0)) 1-1 p ,p . By Proposition 3.5, there exists a unique u ∈ M R(p, H) such that

Thus, by Proposition 3.5, there exists a constant C such that

By continuity at 0, for ǫ > 0 there exists t 0 > 0 such that for t ∈ [0, t 0 ] (B(0) -B(t)) L(H) < ǫ.

Hence for τ = t 0 small enough, the operator S is a contraction on M R(p, H) and so it has a fixed point u ∈ M R(p, H). Clearly, u is a solution of the Cauchy problem (3.16) on [0, t 0 ]. The uniqueness of u follows from the has maximal L p -regularity if for every f ∈ L p (0, τ, H) there exists a unique u ∈ W 1 p (0, τ ; H), B(t)u(t) ∈ D(A(t)) a.e. and u satisfies the Cauchy problem in the L p -sense.

Our main result in this section is the following.

Then the Cauchy problem

has maximal L p regularity in H provided u 0 ∈ B(0) -1 (H, D(A(0))) 1-1 p ,p . Moreover there exists a positive constant C such that :

We start with the following lemma.

Lemma 5.2. Under the above assumptions on B(t), the mapp t → B(t) -1 x is differentiable on (0, τ ) with values in L(H) and

for all x ∈ H.

Proof. We write

and since B(t + h) -1 has norm bounded with respect to h it follows that B(t + h) -1 converges uniformly to B(t) -1 . Using this and the fact that t → B(t)x is Lipschitz we obtain the lemma.

Proof of Theorem 5.1. Let f ∈ L p (0, τ ; H) and initial data u 0 such that u 0 ∈ B(0) -1 (H, D(A(0))) 1-1 p ,p . We consider the Cauchy problem with left multiplicative perturbations

Note also that t → -B ′ (t)B(t) -1 +B(t)P (t) is strongly measurable with values in L(H). Thus, we can apply Theorem 3.6. We obtain existence and uniqueness v ∈ W 1 p (0, τ ; H) such that v(t) ∈ D(A(t)) for a.e. t ∈ [0, τ ] which satisfies (5.6). We set u(t) := B(t) -1 v(t). Using Lemma 5.2 we check easily that u ∈ W 1 p (0, τ ; H), B(t)u(t) ∈ D(A(t)) for a.e. t and it is the unique solution of (5.1). Finally, (5.4) follows immediately from the a priori estimate of Theorem 3.6.

Note that we may consider both left and right multiplicative perturbations at the same time. Let B 0 (t) and B 1 (t) be bounded operators satisfying the same assumptions (3.9) and (3.10). We assume that t → B 0 (t) is continuous and t → B 1 (t) is Lipschitz continuous on [0, τ ]. We assume that the forms a(t) and P (t) are as in Theorem 5.1. We consider the Cauchy problem u ′ (t) + B 0 (t)A(t)B 1 (t)u(t) + P (t)u(t) = f (t), u(0) = u 0 .

(5.7)

Then the maximal L p -regularity results of Theorem 5.1 hold for (5.7) for initial data u 0 ∈ B 1 (0) -1 (H, D(A(0))) 1-1 p ,p . The proof is very similar to the previous one. We consider the Cauchy problem with left perturbations      v ′ (t) + B 1 (t)B 0 (t)A(t)v(t) -B ′ 1 (t)B 1 (t) -1 v(t) + B 1 (t)P (t)B 1 (t) -1 v(t) = B 1 (t)f (t) v(0) = B 1 (0)u 0 .

(5.8) We obtain the maximal L p -regularity for (5.8) by Theorem 3.6 and set as above u(t) = B 1 (t) -1 v(t).


Céline Parzani

Stéphane Puechmorel

On a Hamilton-Jacobi-Bellman approach for coordinated optimal aircraft trajectories planning

Keywords: Hamilton-Jacobi-Bellman, Aircraft trajectory coordinated planning, high dimension

In the context of future air traffic management (ATM), an increasing importance is given to environmental considerations and especially fuel consumption. It is thus advisable to make an optimal use of external conditions knowledge, like wind or temperature, to reduce the total fuel needed to complete a flight. On the other hand, safety must be guaranteed all over the trajectory: encounters below the regulatory separation minima, termed as conflicts, must never occur. In this paper, we consider the problem of optimally planning conflict free aircraft trajectories, based on a minimal time criterion and taking into account an ambient wind field. Aircraft motions are restricted to the horizontal plane only in order to reduce the dimension of the problem and to avoid costly changes of flight level. Since global optimality for the whole set of aircraft is sought after, the admissible space is modeled after a Cartesian product of the individual, two-dimensional state spaces with forbidden configurations removed. An Hamilton-Jacobi-Bellman approach in the so-called configuration space obtained that way is then applied to get a coordinated, conflict-free optimal planning. One of the main advantages of this method is the insurance of getting a global optimum without having to rely on an complex initialization procedure. However, as the size of the instances increases with the power of the number of aircraft considered, a dedicated numerical algorithm must be designed. In this work, the Ultra-Bee scheme is adapted to configuration spaces and solved in a four dimensional plus time setting.

Introduction

Aircraft trajectory optimization while ensuring flights safety is, from the operational point of view, a subject of great importance in air traffic management (ATM). It becomes more and more challenging with the expected growth of air traffic for the next 20 years.

As a consequence, many ATM applications aiming at improving the performance of the overall system involve optimal aircraft trajectories planning. This question deals with the generation of a path from departure to destination taking into account several objectives (minimizing distance, motion time or fuel) while avoiding obstacles and complying with aircraft dynamics.

A survey of mathematical models for aircraft trajectory design can be found in [START_REF] Delahaye | Mathematic Models for Aircraft Trajectory Design: A Survey, 2013 ENRI International Workshop on ATM/CNS (EIWAC2013[END_REF]. Planning several trajectories interacting with each other is a complex problem and different approaches can be considered to solve it. Among them, the simplest and most common one is the sequential approach. It consists first in computing the optimal trajectory for a given mobile. Next, the trajectory for a second mobile is computed considering the previous one as a moving obstacle and so on. It has been used for instance in [START_REF] Dougui | A light-propagation model for aircraft trajectory planning[END_REF][START_REF] Schouwenaars | Decentralized Cooperative Trajectory Planning of Multiple Aircraft with Hard Safety Guarantees[END_REF][START_REF] Zheng | A Novel Approach for Multiple Mobile Robot Path Planning in Dynamic Unknown Environment[END_REF].

The two main drawbacks are that a admissible solution may not be found even if one exists and second that a global optimum may not be reached. A second class concerns iterative method (see [START_REF] Chaimatanan | Strategic Deconfliction of Aircraft Trajectories[END_REF]) and consists in first planning individual trajectories, then identifying potential conflicts and finally using metaheuristics like genetic algorithms. This method can reach an asymptotically optimal solution. However, for a given computing time, a feasible (conflict-free) solution is not guaranteed to be reached. A third class of methods is the Model Predictive Control (see for instance [START_REF] Pannequin | Multiple aircraft deconflicted path planning with weather avoidance constraints[END_REF]) which consists in solving at each time step an optimal control problem to go from the current state to the target state satisfying the constraints and then finding the optimal control for each aircraft minimizing the sum, over the whole set of aircraft, of individual cost functions. Moreover, an approach based on Multi-Agents System has also been explored in [START_REF] Oikonomopoulos | Coordination of multiple non-holonomic agents with input constraints[END_REF][START_REF] Sislak | Agent-Based Approach to Air-Traffic Modeling, Simulation and Collision Avoidance[END_REF]. The main drawback is that emergence of the expected complex behavior is not always guaranted. Last, we find the coordinated motion approach which takes into account the whole set of mobiles simultaneously leading to more complex methods. Among deterministic methods, we can cite some works based on navigation functions (see [START_REF] Guys | Automatic conflict solving using biharmonic navigation functions[END_REF]). Moreover, in [START_REF] Alton | Dijkstra-like Ordered Upwind Methods for Solving Static Hamilton-Jacobi Equations[END_REF], Alton optimizes the trajectory of two robots without collision. Instead of solving the problem in R 2 , he worked in the configuration space of several robots. For instance, for two robots, the work space is constructed removing from R 4 the collision points between robots. In our case, the mathematical definition of the configuration space will be given in Section 2.1.

which is defined such that the collision points between robots are excluded from the space. Computed trajectories are then ensured to be without collision. Space dimension increasing exponentially is the main drawback of these algorithms.

In this context, the objective of this paper is to propose a new framework for trajectories planning with coordinated motion ensuring the following key points. First, we are looking for a global solution for simultaneous planning several aircraft trajectories. On the other hand, the model has to be able to avoid collisions and to take into account external components such as wind fields.

To reach our objective, we investigate an optimal control approach. First introduced in [START_REF] Bryson | Applied Optimal Control: Optimization, Estimation and Control[END_REF], it has been widely used in litterature for trajectories planning for which it is well-suited. The objective of such approach is to determine the control inputs such that the response of a dynamical system satisfies the physical constraints, while, at the same time, minimizes some performance criterion. Aircraft trajectory optimization can then be considered as an optimal control problem in which the control variables are time varying. Trajectory feasibility is automatically ensured and reconstruction is then performed from the optimal value of the control at each time step.

Among the methods to solve optimal control problems, we will consider in this work, the Hamilton-Jacobi-Bellman approach. Based on the Dynamic Programming Principle [START_REF] Bellman | Dynamic Programming[END_REF], it can be used to solve fully nonlinear control problems and minimum time problems, taking into account mixed state/control constraints (see [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]).

On the opposite of Direct Approach [START_REF] Betts | Practical methods for optimal control using nonlinear programming[END_REF] or Pontryagin-based optimization algorithms that rely on necessary condition for optimality, HJB approach gives a global optimum avoiding initialization and convergence questions.

For instance, according to [START_REF] Delahaye | Mathematic Models for Aircraft Trajectory Design: A Survey, 2013 ENRI International Workshop on ATM/CNS (EIWAC2013[END_REF], the optimal control formulation of a trajectory optimization problem using Pontryagin's Maximum Principle leads to a Two-point Boundary Value Problem. Numerical techniques such as shooting methods can be applied to solve accurately this Two-point Boundary Value Problem, but their convergence is very sensitive to the choice of an initial guess for the solution. Moreover, an interesting byproduct of the HJB approach is the synthesis of the optimal control in feedback form. Once the HJB equation is solved, for any starting point, the reconstruction of the optimal trajectory can be performed in real time. In return, the method has to face the so-called curse of dimensionality which implies challenging development in numerical methods to make it amenable to implementation. HJB framework has already been used for trajectories planning (see for instance [START_REF] Bokanowski | Optimization of the launcher ascent trajectory leading to the global optimum without any initialization: the breakthrough of the HJB approach[END_REF][START_REF] Bokanowski | Hamilton-Jacobi-Bellman approach for the climbing problem for multi-stage launchers[END_REF][START_REF] Cruck | Collision analysis for a UAV[END_REF]). In this paper, the aim is to demonstrate the potential of HJB approach to plan trajectories for coordinated motion in the presence of wind. Let us note that when wind is taken into account and has a non negligible speed with respect to aircraft, the problem is no longer isotropic. Compared with existing works, the main contribution of this paper is to derive a model able to take into account several mobiles moving simultaneously while the wind effect is described through the dynamic of the system using some control parameters. Moreover, we show how to take advantage of the mathematical formalism to describe conflict avoidance solving HJB system in the configuration space.

In order to evaluate the relevance of our approach, we consider in this study a scenario with only two aircraft. The model is derived following the formalism proposed in [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF] and it leads to the resolution of the HJB equations in a four dimensional space (2D for each aircraft) plus time. Since the aim is to evaluate the HJB approach in term of modeling, we focus on the precision of its numerical solution. Indeed, in order to be able to describe phenomena such as conflict avoidance or obstacles bypassing, the numerical scheme used must provide an accurate description of a front expansion problem. Therefore, we propose to solve the 4D HJB system with the anti-diffusive UltraBee scheme first developed in [START_REF] Desprès | Un schéma non linéaire anti-dissipatif pour l'équation d'advection linéaire. A non-linear anti-diffusive scheme for the linear advection equation[END_REF][START_REF] Desprès | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF]. The CPU time reduction particularly interesting when considering more than two aircraft with an operational objective will be studied in a further work taking advantage of recent numerical developments.

The paper is organized as follows. In Section 2, we first describe the optimal control problem for two aircraft in a wind field. Next, we present the formalism leading to the Hamilton-Jacobi-Bellman system for the problem we are interested in. In Section 3, we present the numerical methods used to solve the problem. First, we recall the Ultra-Bee scheme for the 4-dimensional space case with obstacles. Then, the algorithm for trajectory reconstruction is presented. Numerical experiments showing the relevance of this approach are performed in Section 4 for a two aircraft conflict avoidance scenario. Finally, conclusions of the work and directions for further works are given in Section 5.

Hamilton-Jacobi-Bellman approach for two aircraft trajectories planning

This section aims at describing the physical problem and its mathematical formulation. We are interested in solving a minimum-time path problem for two aircraft under the following assumptions. The two mobiles are moving in a two dimensional space (planar cartesian positions) and aircraft are cruising at a fixed altitude.

Their speed is assumed to be constant and they are only controlled by the heading direction. The objective is then to plan the optimal paths avoiding obstacles and taking into account wind currents. The approach developed in the following section is based on a complex control optimal problem in high dimension solved using HJB equations. This technique has already been studied for trajectories planning (see [START_REF] Bokanowski | Optimization of the launcher ascent trajectory leading to the global optimum without any initialization: the breakthrough of the HJB approach[END_REF][START_REF] Bokanowski | Hamilton-Jacobi-Bellman approach for the climbing problem for multi-stage launchers[END_REF][START_REF] Cruck | Collision analysis for a UAV[END_REF][START_REF] Margellos | Hamilton-Jacobi formulation for Reach-Avoid Problems with an application to Air Traffic Management[END_REF]) and to formulate the model we will follow the mathematical framework proposed in [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF][START_REF] Bokanowski | Minimum time control problems for non autonomous differential equations[END_REF]. The main contribution of this work is to adapt existing works to propose a model able to manage several mobiles taking into account conflict avoidance in the HJB framework. The technique presented in this paper is based on the computation of the backward reachable set.

It consists in computing the set of initial departure position from which aircraft can reach its destinationthe target -avoiding conflicts with other aircraft and obstacles within a given time t ≥ 0. The model description only concerns two aircraft in the following paragraph but it can be extended to a higher number of mobiles. This specific aspect of the problem requires a special numerical treatment and will be addressed in a separate work.

Optimal control problem formulation

Let us consider two aircraft M 1 and M 2 at a given flight level during cruise phase. For i = 1, 2, suppose that at any time t, the position of aircraft M i is a vector X i (t) ∈ R 2 . The optimal control problem consists in computing the optimal coordinated control for the system composed by the two aircraft in order to minimize a criterion.

Let us first define the state variables. In our case, they will correspond to the planar cartesian coordinates of the aircraft M 1 (respectively M 2 ) denoted X 1 (t) (respectively X 2 (t)). Let us note that in this study the initial positions of aircraft as well as the final ones are known.

Resulting from aircraft's dynamics, these state variables are subject to dynamical laws yielding to differential equations on position. For i = 1, 2 let A i be a nonempty compact subset of R. Let f i : R 2 ×A i → R 2 be Lipschitz continuous with respect to (X i , a i ). Moreover, for each aircraft i = 1, 2, we define the set of controls A i

A i := {a i : R + → A i , a i measurable}.
Given an initial state X 0 i corresponding to departure position for aircraft M i and i = 1, 2, the dynamic of the system is given by

           Ẋ1 (t) = f 1 X 1 (t), a 1 (t) t ≥ 0, Ẋ2 (t) = f 2 X 2 (t), a 2 (t) t ≥ 0, X 1 (0), X 2 (0) = X 0 1 , X 0 2 ∈ R 4 , (2.1)
where the controls a 1 ∈ A 1 , a 2 ∈ A 2 must be determined later according to conflict avoidance rules.

In particular, when we consider wind currents, for each aircraft M i with position X i (t) = (x i , y i ) ∈ R 2 at time t, the dynamic of the system is defined as follows:

Ẋi (t) = V h(a i (t)) + W i (t, X i (t)), (2.2)
where V is the constant speed of the aircraft, W is the wind field vector and a i ∈ A i is the control parameter of the optimal path corresponding in our case to the heading angle relative to North direction for aircraft M i .

Here the admissible control set is A i = [0, 2π] and the associated function h is defined by

h : [0, 2π[ → R × R θ → (cos(θ), sin(θ)).
The vector W = (u, v) is given such that u is the east component of the wind, and v is the north component of the wind.

Let us note that in this study, we assume that during cruise phase, aicrafts are moving at a given optimal flight level. Therefore, the vertical component for wind is not taken into account.

In the sequel, we first define the set of controls

A = A 1 × A 2 and the state variable Z = (X 1 , X 2 ) ∈ R 4 .
Then, for any control value a ∈ A, the function Z is solution of the system

   Ż(t) = f (Z(t), a(t)) t ≥ 0, Z(0) = z 0 := (X 0 1 , X 0 2 ). (2.3)
Furthermore, for any control value a ∈ A, we denote the set of associated trajectories on [0, T ] by

K a = Z(t) ∈ R 4 , with t ∈ [0, T ], Z satisfies system (2.3) (2.4)
leading to the whole set of trajectories denoted

K := a∈A K a .
Now we define the state-constraints set that allows to model both conflict avoidance between aircraft as well as fixed obstacles bypassing.

Definition 2.1 (conflict avoidance). We say that two aircraft M 1 and M 2 are in a "conflict avoidance" configuration if the distance between the two at any time is larger than a safety value > 0. Hence we define

C ε 1 ⊂ R 4 as C ε 1 = Z = (X 1 , X 2 ) ∈ R 4 : d(X 1 , X 2 ) > ε (2.5)
the subset of R 4 in which conflicts may not occur.

Definition 2.2 (fixed obstacles bypassing). Let us denote O ⊂ R 4 a set of fixed obstacles then we define C ε 2 the set of constraints enforcing the aircraft to bypass the set of obstacles by

C ε 2 = Z ∈ R 4 ; d(X 1 , O 1 ) > ε, d(X 2 , O 2 ) > ε with (O 1 , O 2 ) ∈ O . (2.6)
We are now able to define the set where the admissible trajectories of aircraft lie. For instance in the case where the aim is to take into account conflict avoidance and fixed obstacles in the model, we have:

Definition 2.3. A trajectory S := {Z(t); t ∈ [0, T ]} will be "admissible" if Z is solution of (2.
3) and satifies the state-constraints on the time interval [0, T ], that means

S ∈ K ∩ (C ε 1 ∪ C ε 2 ).
In this work, we aim at planning aircraft trajectories from departure to arrival positions which corresponds to a front propagation problem for the optimal control problem we consider. Therefore, we will introduce the target C a nonempty closed set of R 4 modeling arrival positions and representing the front at time t = 0. For numerical implementation, we will assume that C is bounded. We have the following definition Definition 2.4. The capture basin Cap f C (t) is the the set of initial points z 0 ∈ R 4 from which starts an admissible trajectory {Z(s); s ∈ [0, t]} solution of (2.3) associated to an admissible control a ∈ A and such that Z reaches the target C at time t ≥ 0 Remark 2.5. Let us note that for all t ≥ 0, the set Cap f C (t) contains the initial positions which can be steered to the target exactly at time t.

We also define the backward reachable set associated to C for the dynamic f by Definition 2.6. The backward reachable set R f C ([0, t]) associated to the target C is the set of all initial conditions z 0 ∈ R 4 from which it is possible to find an admissible trajectory {Z(s); s ∈ [0, t]} where Z is solution of (2.3) reaching the target C before time t ≥ 0 starting from z 0 .

In this work, the objective is to reduce costs associated to cruise phase. These costs mainly depend on two correlated parameters: consumption and flight duration. Let us denote by T the cost to optimize. From the Bolza formulation and for a free final time problem T , the cost functional can be written as

T (z 0 , a) := T 0 g t, Z(t), a(t) dt + l T, Z(T ) . (2.7)
where g represents the instantaneous cost and l represents the cost associated to the final position.

The control problem we are interested in can then be written as follows: for a given z 0 ∈ R 4 , we look for

a ∈ A, Z such that Z (0) = z 0 , Z (T ) ∈ C and {Z (t), t ∈ [0, T ]} ⊂ K a ∩ (C ε 1 ∪ C ε 2 )
where a is solution to

T z 0 , a := min a∈A T z 0 , a (2.8)
Remark 2.7. Let us note that the formulation proposed for problem (2.8) ensure to get the minimal cost for each aircraft. Indeed, the only coupling condition between the two mobiles is expressed through the admissible space K which can be viewed as a forbidden area for the two aircraft. As a consequence, the optimal trajectory for each aircraft can easily be deduced from the 4D trajectory reconstructed from the solution of problem (2.8).

Hamilton-Jacobi-Bellman equation

Assuming that the speed aircraft is constant, the consumption mainly depends on time flight. Hence, neglecting the over-consumption resulting from changes of direction, minimizing time flight leads to minimizing consumption. As a consequence, the criterion to optimize in our study will be the time flight. Therefore, letting g equal to 1 and l equal to 0 in (2.7), the cost T is now the minimal time for each aircraft needed to reach the destination C from the departure position z 0 = (X 0 1 , X 0 2 ) ∈ R 4 with an admissible trajectory obeying the state/control constraints.

The optimization problem then only consists in finding the best control a (t) ∈ A for t ∈ [0, T ] to minimize the total travel time.

In order to solve problem (2.8) on R 2 ×R 2 , we will use the Hamilton-Jacobi-Bellmann approach which consists in proving that T is a solution of a nonlinear time evolution Hamilton Jacobi partial differential equation set

in (t, Z(t)) ∈]0, ∞[×R 4 .
The form of the Hamilton-Jacobi equation is based on Bellman's optimality principle in discrete case [START_REF] Bellman | Dynamic Programming[END_REF]. The resolution of this Partial Differential Equation on the whole space gives a necessary and sufficient condition for optimality. This is a main advantage of the methodology developed in this paper compared to others methods that need suitable initialization and convergence study. It can be seen as a front expansion problem where the wave front represents the minimum time to reach the arrival point. The wave front is computed gradually from the previous front positions.

In order to compute the function T , we will focus on computation of backward reachable sets associated to the target C for the dynamic f of the problem (2.3). In [START_REF] Mitchell | A Time-Dependent Hamilton-Jacobi Formulation of Reachable Sets for Continuous Dynamic Games[END_REF], the authors proved that that the viscosity solution of a particular time-dependent Hamilton Jacobi Isaacs PDE provides an implicit surface representation of the continuous backwards reachable set. From these results (see also [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF]), we have the following proposition: We will take advantage of this property to compute the minimum time T . Let now φ : R 4 → R be Lipschitz continuous and a level set function for the target that means :

φ(Z(t)) ≤ 0 ⇐⇒ Z(t) ∈ C.
We also introduced g a level set function associated to the state-constraints K∩(C ε 1 ∪C ε 2 ) such that g : R 4 → R Lipschtiz continuous and

g(Z(s), s) ≤ 0 ⇐⇒ Z(s) ∈ K ∩ (C ε 1 ∪ C ε 2 ) ∀s ∈ [0, T ], Z(s) ∈ R 4 .
We now consider the value function v defined by

v(z 0 , t) = inf a∈A max φ(Z(t)), max s∈[0,t]
g(Z(s) .

Then it can be proven (see [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF]) that this function v is a level set function for the capture basin

Cap f C (t) such that Cap f C (t) = {z 0 ∈ R 4 |v(z 0 , t) ≤ 0}.
Moreover, we have the following result (see [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF])

Proposition 2.9. The value function v is the unique continuous viscosity solution (in the sense of [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]) of the following Hamilton-Jacobi-Bellman (HJB) equation

   min ∂ t v + H Z(t), ∇ Z v ; v -g(Z(t)) = 0 t > 0, Z(t) ∈ R 4 v(Z(t), 0) = max(φ(Z(t)), g(Z(t))), Z(t) ∈ R 4 (2.9)
where the Hamiltonian is given by

H(Z(t), ∇ Z v) := max ã∈ Ã (-f (Z(t), ã) • p) p ∈ R 4 (2.10)
Let us now detail the computation of the minimum time function T we aim to minimize and defined by (2.8). From [START_REF] Bokanowski | Reachability and minimal times for state constrained nonlinear problems without any controllability assumption[END_REF], it can be shown that

T (z 0 , a ) = inf t ∈ [0, T ], v(z 0 , t) ≤ 0. (2.11)
Afterwards, the minimal time function is used to recover the optimal control and optimal trajectories from any point of the state space.

Numerical resolution

3.1. Ultra-Bee scheme for HJB resolution

In this section, we are interested in the numerical approximation of the first order Hamilton Jacobi Bellman (HJB) equation resulting from aircraft trajectories planning given by (2.9)-(2.10). A large variety of numerical methods have been studied to determine approximations of the viscosity solutions of systems of type (2.9)-(2.10) which remains a challenging problem, in particular in high dimensions.

Some of them deal with finite difference methods as the first order monotone schemes proposed in [START_REF] Crandall | Two approximations of solutions of Hamilton Jacobi equations[END_REF] or higher ones (ENO scheme in [START_REF] Osher | High order essentially non-oscillatory schemes for Hamilton-Jacobi equations[END_REF] and WENO scheme in [START_REF] Zhang | High order WENO schemes for Hamilton-Jacobi equations on triangular meshes[END_REF]). A second class of methods concerns Discontinuous Galerkin methods, a direct method was proposed in [START_REF] Cheng | A discontinuous Galerkin finite element method for directly solving the Hamilton-Jacobi equations[END_REF] and a scheme for front propagation with obstacles in [START_REF] Bokanowski | A discontinuous Galerkin scheme for front propagation with obstacles[END_REF]. Then, we can also consider Semi-Lagrangian schemes which are based on the discretization of the dynamic programming principle as developed in [START_REF] Falcone | Semi-Lagrangian schemes for Hamilton-Jacobi equations, dis-crete representation formulae and Godunov methods[END_REF][START_REF] Falcone | Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations[END_REF] and [START_REF] Carlini | An efficient algorithm for Hamilton-Jacobi equations in high dimensions[END_REF] for instance. A brief review of different efficient techniques that have been proposed can also be found in [START_REF] Cacace | Can local single-pass methods solve any stationnary Hamilton-Jacobi-Bellman equation[END_REF].

In this paper, we use the Ultra Bee scheme (a finite difference type scheme) to solve the HJB equations (2.9)-(2.10). First developed for advection equations with constant velocity [START_REF] Desprès | Un schéma non linéaire anti-dissipatif pour l'équation d'advection linéaire. A non-linear anti-diffusive scheme for the linear advection equation[END_REF][START_REF] Desprès | Contact discontinuity capturing schemes for linear advection and compressible gas dynamics[END_REF], a generalization to velocity with changing sign is suggested in [START_REF] Bokanowski | Anti-diffusive schemes for linear advection and application to Hamilton-Jacobi-Bellman equations[END_REF] and the adaptation to solve HJB equations is done in [START_REF] Bokanowski | An anti-diffusive scheme for viability problems[END_REF]. Due to very good anti-diffusive properties, this scheme allows computations of the reachable sets with good accuracy and therefore to get a precise localization of the discontinuity. Indeed, the problem we are interested in can be viewed as a front expansion problem where the wave front represents the minimum time to reach a position, which corresponds to a discontinuity from numerical point of view.

Let us now recall the algorithm applied to discretize the problem (2.9)-(2.10

). Let D = [X 1,min , X 2,min ] × [Y 1,min , Y 2,
min ] be a bounded square box corresponding to the discrete computational domain. We consider a regular spatial grid G of D ⊂ R 4 defined by the constant mesh size ∆x in all directions. Then, for I = {i, j, k, l} ∈ ZZ 4 , let

ξ I := (i∆x, j∆x, k∆x, l∆x)
be the center of the cell

M I = [x 1,i-1 2 , x 1,i+ 1 2 ] × [y 1,j-1 2 , y 1,j+ 1 2 ] × [x 2,k-1 2 , x 2,k+ 1 2 ] × [y 2,l-1 2 , y 2,l+ 1 2 ]
We solve the HJB equation on the time interval [0, T max ]. Let us denote ∆t the time step and for n ≥ 0,

t n := n∆t.
The following steps are applied.

Let us first starts by the initialization step. The finite-volume like UltraBee scheme aims at computing a numerical approximation V n I of the mean value of v solution of (2.9)-(2.10) on cell M I at time t n . Thus,

V n I ≈ 1 m(M I ) M I v t, ζ dζ
The initialization is then done using the average values of ϕ:

V 0 I = 1 m(M I ) M I ϕ(ζ)dζ,
the average values of g on the grid are computed following

W I = 1 m(M I ) M I g(ζ)dζ
and we consider an approximation S of λ such that for ζ ∈ D

S(ζ) =    {1} if ζ / ∈ C, {0, 1} if ζ ∈ C
We recall that C is a non empty closed set of R 4 and represents aircraft destination in our application.

The next step consists in the set controls discretization.

The set of controls A = A 1 × A 2 is discretized using N 1 N 2 values and we set

A d := {a 1 1 , ..., a N1 1 } × {a 1 2 , ..., a N2 2 }.
Now, let us detail the time evolution step of the algorithm. For n ≥ 1, let us assume that V n I , solution of (2.9)-(2.10) on cell M I at time t n , is known. Moreover, for a d ∈ A d and λ d an approximation of λ(ξ I ), let us denote U n+1 I a d , λ d the numerical solution of the advection equation

∂u ∂t (t, Z(t)) -λ(Z(t))f (t, a) • ∇ Z u(t, Z(t)) = 0. (3.12)
The solution of this advection problem is obtained using the Ultra-Bee Scheme. Since 4 dimensions have to be considered (2D for each mobile) in our study, we will use a classical Trotter splitting or "alternate direction method" such that the numerical solution evolves during a time step in each direction. The resolution using this splitting technique is stable under the following CFL condition, max

I,a max ∆t ∆x |f 1 (ξ I , a 1 )|, ∆t ∆x |f 2 (ξ I , a 2 )| ≤ 1 (3.13)
For sake of simplicity, let us present the UltraBee scheme only for the first step of the splitting corresponding to the x 1 -direction. Thus, ∀i ∈ ZZ and for (j, k, l) fixed in ZZ 3 , let us note U n i,. = U n I , then we have

U n+1 i,. (a d , λ d ) = U n i,. - ∆t ∆x f 1 (ξ I , a 1 ) F L i+ 1 2 -F R i-1 2
Let us denote by ν i = ∆t ∆x f 1 (ξ I , a 1 ) and define the "fluxes"

F L i+ 1 2 and F R i-1 2
. Following the extension of the original scheme in the case where the velocity is changing of sign proposed in [START_REF] Bokanowski | Anti-diffusive schemes for linear advection and application to Hamilton-Jacobi-Bellman equations[END_REF], they are given by

1. if ν i > 0, F L i+ 1 2 = min max(U n i+1,. , b n,+ i ), B n,+ i where b n,+ i and B n,+ i are defined by        b n,+ i = max(U n i,. , U n i-1,. ) + 1 ν i U n i,. -max(U n i,. , U n i-1,. ) , B n,+ i = min(U n i,. , U n i-1,. ) + 1 ν i,. U n i,. -min(U n i,. , U n i-1,. ) 2. if ν i < 0, F R i-1 2 = min max(U n i-1,. , b n,- i ), B n,- i with        b n,- i = max(U n i,. , U n i+1,. ) + 1 |ν i | U n i,. -max(U n i,. , U n i+1,. ) , B n,- i = min(U n i,. , U n i+1,. ) + 1 |ν i | U n i,. -min(U n i,. , U n i+1,. ) 3. if ν i ≤ 0 ans ν i+1 ≥ 0, then F L i+ 1 2 = U n i,. and F R i+ 1 2 = U n i+1,. 4. last, if ν i ν i+1 > 0 then      F R i+ 1 2 = F L i+ 1 2 if ν i > 0, F L i+ 1 2 = F R i+ 1 2 if ν i+1 > 0
Last, the solution of the HJB equation (2.9)-(2.10) on the cell M I is then obtained following

V n+1 I = max W I , min λ d ∈S(ξ I ) a d ∈A d U n+1 I a d , λ d (3.14)

Optimal trajectories reconstruction

The optimal trajectory reconstruction is obtained by minimizing the time function along possible trajectories (see [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF], Appendix by Falcone). It is then necessary to know the minimal time function T on all over the domain D and not only around the front at the final time. In practice, the values of T are saved during the computation of the solution v of (2.9)-(2.10). We follow the process proposed in [START_REF] Bokanowski | Minimum time control problems for non autonomous differential equations[END_REF]. Thus, once we have solved the HJB equation (2.9)-(2.10), we get the reachability function v and we use (2.11) to reconstruct T from v. A numerical approximation T I of T (ξ I ) is computed as follows. We first set T I equal to a large numerical value everywhere on the domain D. Then, as the computation of V n I goes on, for n = 0 up to n = N -1, (N being the number of iterations in time) we consider that if

(V n I > 0, V n+1 I = 0 and T n I = ∞) then we set T n+1 I := t n+1 .
In other words, T I is the first time when the value of the cell M I goes to 0, otherwise T I = ∞.

Numerical Simulations in 4D

This section is devoted to the evaluation from a numerical point of view of the HJB approach proposed in this paper. So, we present 4-dimensional numerical simulations for different scenarios corresponding to configurations arising in air traffic management area and we show how the method allows both to take into account wind current effects, fixed obstacles and conflict avoidance.

For each scenario, we investigate optimal trajectories reconstruction for a two-aircraft configuration in the presence or not of constraints.

4D numerical method validation

This first part is devoted to the validation of the numerical method applied to discretize the HJB equation (2.9). We consider here two simple configurations for which the behavior of the solution is known a priori.

Thus, we study two aircraft M 1 and M 2 each one flying in a square domain [-1 : 1] × [-1 : 1] without wind current. The dynamics for each aircraft is such that (X i (t)) = (cos(θ i ), sin(θ i )). We use 15 points in each space direction and 10 control values in [0, 2π] while the time step satisfies the CFL condition. The results obtained are plotted in Figure 1 for two cases: the first one (left side) is such that departure positions of the aircraft do not lead to potential conflicts. This is satisfied since as expected under assumptions detailed above, each of the two aircraft follows a straight line trajectory to its respective arrival position. The second case plotted on the right side in Figure 1 illustrates a conflict configuration since departure positions are chosen such that aircraft still follow a straight line trajectory but are at the same position (0, 0) at the same time. We then consider that the 4D numerical method is validated.

Optimal trajectories under wind effect

In this part, we present a scenario inspired from the Zermelo problem with two objectives: the first one is an additional validation of the method in the presence of wind and the second one is to show the effect of wind current on trajectories behaviors. The Zermelo's problem first introduced in 1931 consists in considering a ship that has to travel a region with strong current minimizing the travel time from a point A to a point B. The ship has a constant speed and its heading can change such that the dynamics is given by ẋ = V cos(θ) + u(x, y),

ẏ = V sin(θ) + v(x, y). (4.15)
It is then similar to the one we use in our model defined by equations (2.2). Moreover, for some specific form of wind current, an analytic solution can be computed in 2D.

In this scenario, we consider two aircraft M 1 and M 2 each one moving in a square domain [-1

: 1] × [-1 : 1].
Since we are interested in validation and wind effect, we do not consider any interaction, that means no conflict, between the two aircraft neither fixed obstacle in the computational domain. The first aircraft M 1 is going from its departure position (0.7, 0.2) to the arrival one (-0.7, 0.2) and respectively for the second aircraft M 2 going from (-0.7, -0.2) to (0.7, -0.2). The wind current is such that V = 0.5, u(x, y) = -0.2 * y and v(x, y) = 0 in system (4.15). Thus, wind is going from east to west in the upper part of the domain, growing gradually from y = 0 to the North direction and symmetrically, wind is going from west to east in the lower part of the domain, growing gradually from y = 0 to the South direction. Figure 2 presents the results obtained with the HJB system (2.9) with D = [-1, 1] 2 × [-1, 1] 2 discretized with the Ultra-Bee scheme presented in Section 3.1. We use 15 points in each space direction, 10 values of control and the time step is satisfying the CFL defined by (3.13).

We plot the trajectories reconstructed following the method described in Subsection 3.2 and then smoothed.

We first note that the method leads to the optimal trajectory for each aircraft. As expected according to the analytic solution of the Zermelo problem in 2D, reconstructed trajectories follow a parabolic curve. Moreover, symmetry in wind definition leads to symmetry in reconstructed trajectories (more or less depending on the smoothing option as also expected). This allows to validate the proposed approach taking into account several mobiles. Now, looking at the effect of the wind, we can observe that trajectories are deviating from the no-wind straight lines ones (plotted in dash lined on Figure 2) to take advantage of the wind. Let us note that the longer distance is counterbalanced by the higher ground speed of the aircraft such that the arrival time is less optimal in the presence of wind.

Optimal trajectories under wind effect taking into account fixed obstacles

In this part, we will consider a more complex scenario composed of two aircraft flying through wind current with some fixed obstacles to avoid. From an operational point of view, these fixed obstacles can described for instance dangerous meteorological areas or areas for which flying over is forbidden. Here, the computational do-

main is D = [-1.2, 1.2] 2 ×[-1.2, 1.2] 2 .
The first aircraft M 1 is going from its departure position (-0.96, -0.48) to the arrival one (0.96, 0.96) and respectively for the second aircraft M 2 going from (-0.48, -0.96) to (0.96, 0.48).

The wind current we consider here is similar to the one used for scenario 1 except that it moves in the opposite direction with V = 1, W x = -0.5 * y and W y = 0 in system (2.2). The obstacles are defined by circles respectively centered in (0.4, 0.56), (0.4, -0.56) and (-0.56, -0.56) with radius equal to 0.16. We use 15 points in each space direction, 20 values of control and the time step is satisfying the CFL defined by (3.13).

Figure 3 presents the results obtained in two cases. First we plot the reconstructed trajectories taking into account obstacles without wind, this corresponds to the trajectories plotted in green line for aircraft M 1 and in red line for aircraft M 2 . Next, we plot the reconstructed trajectories with wind and obstacles in green line with circles for aircraft M 1 and respectively in red line with squares for aircraft M 2 . Let us note that once again in order to evaluate the potential of the HJB approach to take into account obstacles, departure and arrival positions for the aircraft have been chosen in order that no conflict happens in this scenario. We can first observe that the method well performs since the fixed obstacles are avoided by the aircraft. Moreover, we show that the two trajectories do not behave in the same way to get around the obstacles due to the wind direction. Indeed, the trajectory for aircraft M 1 takes more advantage of the wind while avoiding obstacles and we can note that wind-optimal trajectories can be far from the optimal trajectories computed without wind.

Optimal trajectories with conflict avoidance

In this subsection, we are now interested in evaluating the potential of the proposed approach to treat conflict avoidance. For this purpose, we will study three cases: a first one with no wind and no obstacle (see FigureTest3), a second one with wind effect and no obstacle (see FigureTest4) and the last one with no wind and obstacle (see FigureTest4).

4.4.1. Optimal trajectories with conflict avoidance, no wind and no obstacle

Once again we consider a two-aircraft configuration and in order to dissociate different effects, we consider here that there is no wind current. The first aircraft M 1 is going from its departure position (-0.7, 0.7) to the arrival one (0.7, -0.7) and respectively for the second aircraft M 2 going from (-0.7, -0.7) to (0.7, 0.7). This configuration can illustrate for instance aircraft that are just landing and are going to park in an assigned place around the airport. These positions have been chosen such that if the constraint on conflict avoidance defined by (2.5) is not included in the model then the two-aircraft will be in conflict at point (0, 0). This is illustrated with the dashed lines of Figure 4 left side. The dynamics for each aircraft is such that (X i (t)) = (cos(θ i ), sin(θ i )).

We use 15 points in each space direction and 10 control values in [0, 2π] while the time step satisfies the CFL condition. In Figure 4, we plot the global trajectories (left side) obtained taking into account the conflict avoidance conflict in the model for each aircraft (blue line with circles for aircraft M 1 and magenta line with squares for aircraft M 2 ). Each circle or square correspond to a time, the 4D trajectory has been evaluated with the reconstruction process. In order to verify that with the conflict avoidance constraint, trajectories do not cross, we plot a zoom on the conflict area in Figure 4 (right side) with the corresponding times and we show that aircraft are never at the same position in the same time in the 4D configuration space as required.

Optimal trajectories with conflict avoidance taking into account wind

In order to evaluate the relevance of the method proposed in this work, we now study a new configuration adding wind effect. Aircraft initial positions are identical with those defined in previous section 4.4.1.However, this time the dynamic for each aircraft is given by (X i (t)) = (-0.2 * x 2 + 0.5 * cos(θ i ), sin(θ i )). In Figure 5 we plot trajectories without conflict avoidance on the left side and with conflict avoidance on the right one and once again, we demonstrate the potential of the method to avoid conflict in the presence of wind.

Optimal trajectories with conflict avoidance taking into account fixed obstacles

The last test concerns conflict avoidance in the presence of fixed obstacle. We still use the computational domain and the departure and arrival positions of the two previous cases for aircraft. Wind is not consider in 335 this test but we add a squared fixed obstacle of size 0.3 centered at (0, 0). the dynamics for each aircraft is (X i (t)) = (cos(θ i ), sin(θ i )) with the same numerical parameters as in section 4.4.1. As previously, we plot in Figure 6 trajectories without conflict avoidance on the left side and with conflict avoidance on the right one and once again and show that the method is also able to avoid conflict in the presence of fixed obstacles.

Conclusion and perspectives

In this paper, we propose a new framework based on optimal control and Hamilton Jacobi Bellman approach for trajectories planning. The model we develop is able to take into account for several mobiles, conflict avoidance, wind current and fixed obstacles bypassing. We use the Ultra-Bee numerical scheme in order to discretize the 4D space domain plus time and present several numerical results that demonstrate the potential of the proposed approach for different two aircraft configurations. Compared to others existing methods, the main advantage of this approach is that it provides a global optimal solution without initialization and convergence questions. However, it has to faced the so-called curse of dimensionality. As a consequence a possible future work could be to extent the model to a configuration involving a larger number of mobiles taking advantage of recent developments of efficient algorithms to solve HJB equation. Another future direction to increase the number of aircraft taken into account in the planning, could concern the possibility to initialize PMP (Pontryagin Minimum Principle Method) using the trajectories resulting from the HJB approach. Indeed, following the method proposed in [START_REF] Cristiani | Initialization of the shooting method via the Hamilton-Jacobi-Bellman approach[END_REF], a rough approximation of the value function computed by the HJB approach can be used to obtain an initial guess for the PMP method. In this case, due to the properties of the HJB equations, the main advantage is that the initial guess will be close to the global minimum.

Proposition 2 . 8 .

 28 Introducing a new virtual control λ(.) with ∀s ∈ [0, t] λ(s) ∈ [0, 1], the backward reachable setR f C ([0, t])defined by (2.6) is exactly the capture basin Cap f C (t) associated to C and the dynamic f (z 0 , (a, λ) := λf (z 0 , a) for all (a, λ) ∈ Ã = A × [0, 1].

Figure 1 :

 1 Figure 1: Optimal trajectories for two aircraft with no wind and no obstacles: no conflict case (left side) and conflict configuration (right side)

Figure 2 :

 2 Figure 2: Optimal trajectories for two aircraft under wind effect.

Figure 3 :

 3 Figure 3: Optimal trajectories for two aircraft taking into account fixed obstacles under wind effect.

Figure 4 :

 4 Figure 4: Optimal trajectories for two aircraft taking into account conflict avoidance. Global trajectories (left side) and zoom on the conflict zone (right side).

Figure 5 :

 5 Figure 5: Optimal trajectories for two aircraft taking into account conflict avoidance under wind effect. Trajectories in conflict (left side) and trajectories when the conflict avoidance constraint is included in the model (right side).

Figure 6 :

 6 Figure 6: Optimal trajectories for two aircraft taking into account fixed obstacles and conflict avoidance. Trajectories in conflict (left side) and trajectories when the conflict avoidance constraint is included in the model (right side).