Jean-Baptiste Caillau

Roberto Ferretti

Emmanuel Trélat
email: emmanuel.trelat@sorbonne-universite.fr

Hasnaa Zidani
email: hasnaa.zidani@insa-rouen.fr

Numerics for finite-dimensional optimal control problems

Keywords: Optimal control, Pontryagin maximum principle, Hamilton-Jacobi-Bellman equation, direct and indirect methods, optimistic planning

We survey on numerics for finite-dimensional nonlinear optimal control. The chapter is written as a guide to practitioners who wish to get rapidly acquainted with the main numerical methods used to efficiently solve an optimal control problem. We consider throughout two classical examples, quite simple but representative enough to be complexified and generalized to other problems: the Zermelo and the Goddard problems. We provide their solving codes that are available on the web and make the point on the most up-to-date and efficient methods existing nowadays. We range on direct and indirect methods, on Hamilton-Jacobi approaches and we end with optimistic planning. Our examples illustrate the pros and cons of the methods and we also show how those various approaches can be combined in view of augmenting the efficiency of the numerical solving.

1 Introduction and statement of the problem

Brief overview

Optimal control theory has been widely developed since many decades. The theoretical and numerical achievements, motivated by a body of diverse applications in various domains, provide valuable insights into the nature of optimal controls and the corresponding optimal trajectories. In this chapter our objective is to provide to practitioners a guid to get rapidly acquainted with the most powerful but however easy-to-use methods and algorithms to solve efficiently a given nonlinear optimal control problem in finite dimension.

The most intuitive and popular numerical methods for solving an optimal control problem, called direct methods, consist of first discretizing then optimizing. Such approaches have been investigated in a number of contributions (see, e.g., [START_REF] Betts | Practical methods for optimal control using nonlinear programming, volume 36 of Advances in Design and Control[END_REF]). From the theoretical point of view, the efficiency of these methods has been established for some classes of control problems (see [START_REF] Bonnans | Computation of order conditions for symplectic partitioned Runge-Kutta schemes with application to optimal control[END_REF][START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF][START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF][START_REF] Ross | Legendre pseudospectral approximations of optimal control problems[END_REF][START_REF] Sanz-Serna | Symplectic Runge-Kutta schemes for adjoint equations, automatic differentiation, optimal control, and more[END_REF]). From the numerical point of view, the direct approach benefits from the tremendous advances in numerical optimization methods achieved in the last decades. As shown 1 in Section 2, the direct approach is extremely easy to implement for general control problems with constraints on the control variable and on the state variable as well. However, in general, direct methods may provide only locally optimal solutions and may require a fair initialization of the iterative process in optimization algorithms. Also, they may lack of numerical accuracy, or may become computationally demanding in high dimension.

A major breakthrough in optimal control theory were achieved in the 1950's by Pontryagin's research group, who successfully generalized and extended to a general nonlinear optimal control setting the classical Euler-Lagrange and Weierstrass conditions of the Calculus of Variations. The early optimality conditions, called Pontryagin Maximum principle (in short, PMP), were subsequently strengthened and extended using methods of convex and non-smooth analysis, or methods of differential geometry (see [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF]). The PMP inspired effective computational schemes, like the shooting method presented in Section 3. As illustrated on some classical examples, the shooting method provides very accurate numerical solutions. The major drawback of this method is that it requires an a priori knowledge of the structure of the solution as well as a good approximation of the adjoint state. In particular, shooting methods are not much suitable for control problems in presence of state constraints or at least, are much more difficult to implement in particular if there are many state constraints.

Another major advance in optimal control was achieved by Richard Bellman in the 1950's, who provided a description of how the minimum cost depends on initial conditions, as the solution to a Hamilton-Jacobi partial differential equation. The Hamilton-Jacobi approach, when it can be applied, yields a global solution to the optimal control problem. It also provides the optimal control in a feedback form suitable for many engineering applications. Despite of these advantages, the approach suffers from some difficulties of computing the solution to the Hamilton-Jacobi equation in higher dimensions. The numerical simulations presented in Section 4 show that the approximation of the HJB equations, even on coarse grids, provides solutions which depict well the qualitative structure of the optimal trajectories and may thus be used to guess the structure of optimal solutions. However, to get an accurate computation of the optimal trajectories, an important numerical effort is necessary to approach the solutions of the HJB equations on very fine grids.

More recently, other global methods have been developed for control problems, such as the optimistic planning (OP) algorithms introduced in Section 5. These methods are based on a discretization process in the control space and do not require any discretization of the state space. As a consequence, OP methods are particularly suitable for control problems with lowdimensional control space. Moreover, these approaches have the potential to be applied to problems where the dynamics and the cost are given by learning models. A preliminary analysis of the complexity of the OP methods is now well established, but it remains to further develop the approach to make it more accurate, to integrate the knowledge of the structure of trajectories when available and also for a more efficient implementation.

Throughout the chapter we consider two well known but representative optimal control examples: the Zermelo and the Goddard problem. Both are quite simple but can be complexified and adapted to more intricate models. We use these two examples to illustrate the numerical methods and show how they can be rapidly and efficiently implemented, with the most up-to-date existing solutions. We also show how the various approaches can be combined to give the best of themselves. For instance, HJB approach and direct methods can be used to obtain a rough solution that successfully provides an initial guess for the (more accurate) shooting method. We hope that the codes that we provide, which are also available on the web (as indicated further), can serve as templates to readers interested in adapting them to their specific setting.

Notations. Throughout the chapter, R denotes the set of real numbers, •, • and • denote respectively the Euclidean inner product and the norm on R N (for any N ≥ 1), B N = {x ∈ R N : x ≤ 1} is the closed unit ball (also denoted B if there is no ambiguity) and B(x; r) = x + rB.

For any set S ⊆ R N ,

• S, S, ∂S, coS denote its interior, closure, boundary, and convex envelope, respectively. For any a, b ∈ R, we define a b := max(a, b). Similarly, for a 1 , • • • , a m ∈ R, we define m i=1 a i := max(a 1 , • • • , a m ). The notation W 1,1 ([a, b]) stands for the usual Sobolev space {f ∈ L 1 (a, b), f ∈ L 1 (a, b)}. Finally, the abbrevation "w.r.t." stands for "with respect to", and "a.e." means "almost everywhere".

Formulation of the optimal control problem

Let d, r ∈ N * , let T > 0 be a fixed final time horizon and let U be a compact subset of R r (with r ≥ 1). We consider the finite-dimensional control system in R d (for 0 ≤ t < T ) ẋ(s) = f (s, x(s), u(s)), a.e. s ∈ (t, T ),

where the control input u : [0, T ] -→ R r (r ≥ 1) is a measurable function such that u(s) ∈ U for almost every s ∈ [0, T ]. Throughout the paper, we assume that

(H 0 ) U is a compact subset of R r .
We denote by U the set of all admissible controls U := {u : [0, T ] -→ R r measurable, and u(s) ∈ U a.e.}.

The dynamics f : [0, T ] × R d × U → R d and the distributed cost : [0, T ] × R d × U → R satisfy:

(H 1a )                (i) f is continuous on [0, T ] × R d × U ;
(ii) x → f (s, x, u) is locally Lipschitz continuous in the following sense:

∀R > 0, ∃k R ≥ 0, ∀(x, y) ∈ (B d (0; R)) 2 , ∀(s, u) ∈ [0, T ] × U f (s, x, u) -f (s, y, u) ≤ k R x -y ; (iii) ∃c f > 0 such that f (s, x, u) ≤ c f (1 + x ) ∀(s, x, u) ∈ [0, T ] × R d × U ; and (H 1b )                (i) is continuous on [0, T ] × R d × U ;
(ii) x → (s, x, u) is locally Lipschitz continuous in the following sense:

∀R > 0, ∃k R ≥ 0, ∀(x, y) ∈ (B d (0; R)) 2 , ∀(s, u) ∈ [0, T ] × U | (s, x, u) -(s, y, u)| ≤ k R x -y ; (iii) ∃c > 0 such that | (s, x, u)| ≤ c (1 + x ) ∀(s, x, u) ∈ [0, T ] × R d × U.
The assumptions (ii)-(iii) of (H 1a ) guarantee, for every u ∈ U and (t, x) ∈ [0, T ] × R d , the existence of an absolutely continuous curve x : [t, T ] → R N which satisfies [START_REF] Agrachev | Any sub-Riemannian metric has points of smoothness[END_REF] and the initial data x(t) = x. The Gronwall lemma gives

1 + x(s) ≤ (1 + x )e c f (s-t) ∀s ∈ [t, T ]. (2)
Given any t ∈ [0, T ] and any x ∈ R d , the set of all admissible pairs control-and-trajectories starting at x at time t is denoted by X [t,T ] (x) := {(x, u) ∈ W 1,1 (t, T ) × U | (x, u) satisfies (1) with x(t) = x}.

Throughout the chapter, we consider the optimal control problem

           minimize ϕ(x(T )) + T t
(s, x(s), u(s)) ds, (x, u) ∈ X [t,T ] (x), g(x(s)) ≤ 0 for s ∈ [t, T ], g f (x(T )) ≤ 0, [START_REF] Alla | An adaptive POD approximation method for the control of advectiondiffusion equations[END_REF] with the agreement that inf ∅ = +∞. The final cost ϕ : R d → R and the constraint functions g and g f are given functions satisfying: (H 2 ) ϕ : R d → R is locally Lipschitz continuous and there exists a constant c ϕ ≥ 0 such that |ϕ(y)| ≤ c ϕ (1 + y ).

(H 3 ) The constraint functions g : R d → R m and g f : R d → R r are locally Lipschitz continuous and there exists a constant c g ≥ 0 such that, for every y ∈ R d ,

|g(y)| + |g f (y)| ≤ c g (1 + y ).
The so-called augmented control system ẋ(s) = f (x, x(s), u(s)) a.e. s ∈ (t, T ),

ż(s) = -(s, x(s), u(s)) a.e. s ∈ (t, T ),

is usually considered in optimal control theory to recast the Bolza problem in the Mayer form

       minimize ϕ(x(T )) -z(T ), (x, u) ∈ X [t,T ] (x 0 ), g(x(s)) ≤ 0 for s ∈ [t, T ], g f (x(T )) ≤ 0, (5)
In what follows, we will sometimes assume that the augmented dynamics satisfies the following assumption (convex epigraph):

(H 4 ) For any s ∈ [0, T ] and every x ∈ R d , f (s, x, u) -(s, x, u) + η , u ∈ U, -c (1 + x ) + (s, x, u) ≤ η ≤ 0 is a convex set. Under assumptions (H 0 ), (H 1a ), (H 1b ), by [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF], the set S [a,b] (x, z) is bounded in W 1,1 ([a, b]). Moreover, if (H 4 ) is satisfied then S [a,b] (x, z) is a compact set in W 1,1 ([a, b]) endowed with the C 0 ([a, b])-topology (see [START_REF] Aubin | Differential inclusions. Set-valued maps and viability theory[END_REF]Theorem 1.4.1]). Therefore, if there exists a trajectory x ∈ S [0,T ] (x 0 , 0) that satisfies the constraints g(x(s)) ≤ 0 for all s ∈ [t, T ] and g(x(T )) ≤ 0, then the the control problem [START_REF] Altarovici | A general Hamilton-Jacobi framework for non-linear state-constrained control problems[END_REF] has an optimal solution.

Examples

In this section, we present two well-known examples that we are going to consider throughout. In our opinion, they illustrate nicely the most classical difficulties encountered in theoretical and numerical optimal control, and lend themselves to a number of more complicated variants. They are expected to serve as "templates" to the reader who aims at getting acquainted with the main issues in numerical optimal control.

Example 1: Zermelo problem A boat with coordinates x(s) = (y 1 (s), y 2 (s)) navigates through a canal R × [a, b], starting at x(0) = x = (x 1 , x 2
). We will either take a = 0 and b = 1, or a = -2 and b = 2. One may want to reach a fixed final point, or to reach an island C, with minimal cost. The cost function may be an energy, or the final time, or the lateral deport, etc. The control system is

ẏ1 (s) = v(s) cos(u(s)) + h(y(s)), (6a)
ẏ2 (s) = v(s) sin(u(s)), (6b)
where u(s) ∈ [0, 2π] is the first control (angle), v(s) ∈ [0, V max ] is a second control (speed of the boat), and h(y(s)) is the current drift (along the x 1 -axis). Because of the drift term (which can be greater than V max ), the system may not be controllable. Consider a target C := B(c, r 0 ) that is a ball with radius r 0 ≥ 0 and centered at a given point c located in the canal. The target represented by a function ψ defined by

g f (x) := x -c -r 0 as x ∈ C ⇐⇒ g f (x) ≤ 0.
Consider also a set of constraints K := {x ∈ R 2 , g(x) ≤ 0} where g is a given function that is non-positive in a region where the boat can move and g is positive in the location of the obstacles that the boat should avoid. In this example, the cost function could be the time, or the energy, required to steer the boat from a given position x to the target C.

Example 2: Goddard Problem We consider the optimal control problem associated with a vertical ascension flight of a rocket, known as Goddard problem. The dynamics of the rocket is defined with three state variables: r, the altitude, v the relative speed and m, the total mass. In general, a dimensionless version for the motion equations is considered in the literature:

ṙ = v (7) v = 1 m (T max u -D(r, v)) - 1 r 2 (8) ṁ = -T max b u (9)
where D(r, v) is the drag force and T max •u is the trust force. The dimensionless initial state of the rocket is given by r(0) = 1.0, v(0) = 0, m(0) = 1. Here r = 1 corresponds to the Earth's ground level and m = 1 to the initial total mass of the rocket. The motion of the rocket is controlled by the trust factor u ∈ [0, 1] so that the trust force is on the interval

T = T max • u ∈ [0, T max ].
The drag force is a nonlinear function of r and v. Its expression depends on the choice of a model for the atmosphere and on the structure of the rocket: D(r, v) = C D ρ(r)|v|v, where C D is the drag coefficient of the rocket and ρ is the atmospheric density. In this section we consider the case of constant drag coefficient (it depends in general on the Mach number) and exponential model of the atmospheric density:

D(r, v) = C D v 2 e -β (r-1) .
This definition of the rocket's model has been widely studied in the literature. We take C D = 310.0, β = 500.0, T max = 3.5, b = 2.0, as in [START_REF] Tsiotras | Goddard problem with constrained time of flight[END_REF] and other references. The corresponding model is then an approximation of real flight conditions for some rockets. The optimal control problem consists of maximizing the final altitude r(t f ) at a given fixed final time t f > 0 under the following final constraint (which represents a limit on the fuel consumption): m(t f ) ≥ m * , with m * = 0.6. Finally, the optimal control problem is formulated as:

inf u(.) -r(t f ) ṙ = v r(0) = 1 v = 1 m (T max u -D(r, v)) - 1 r 2 v(0) = 0 ṁ = -T max b u m(0) = 1 m * -m(t f ) ≤ 0
2 Direct methods: nonlinear programming

Principle

We call direct methods all numerical methods consisting of numerically solving the optimal control problem as follows: without applying a priori any first-order necessary condition for optimality, we choose a discretization for the state and for the control, we choose a numerical scheme to discretize the control system (differential equation) and the cost functional (integral quadrature), so that the discretized optimal control problem is expressed as a family of nonlinear optimization problems in finite dimension, indexed by a discretization parameter N , of the form

min Z∈C F (Z) (10)
where Z = (x 1 , . . . , x N , u 1 , . . . , u N ) and

C = {Z | g i (Z) = 0, i ∈ 1, . . . , r, g j (Z) ≤ 0, j ∈ r + 1, . . . , m}. (11)
This is a classical optimization problem under constraints in finite dimension, with a dimension being as larger as the discretization is finer. Of course, there exist an infinite number of variants to discretize the problem and ending up with a problem of the form [START_REF] Augoula | High order numerical discretization for hamilton-jacobi equations on triangular meshes[END_REF]. We discuss hereafter several classes of discretizations. Once this transcription has been done, one can then implement a number of various optimization routines to solve [START_REF] Augoula | High order numerical discretization for hamilton-jacobi equations on triangular meshes[END_REF].

Let us first explain hereafter one of the simplest possible discretizations. Consider the optimal control problem (3) with t = 0 as initial time. Consider a subdivision

0 = t 0 < t 1 < • • • < t N = T of the interval [0, T ].
Controls are discretized on U -valued piecewise constant functions on this subdivision. To discretize the control system, we choose the explicit Euler method: setting

h i = t i+1 -t i , we have x i+1 = x i + h i f (t i , x i , u i ) for i = 0, . . . , N -1.
To discretize the integral cost, we choose the left rectangle method (which is equivalent to the explicit Euler method for the augmented system). We obtain the nonlinear programming problem min C(x 0 , . . . , x N , u 0 , . . . , u N ),

x i+1 = x i + h i f (t i , x i , u i ), u i ∈ U, g(x i ) ≤ 0, i = 0, . . . , N -1, g f (x N ) ≤ 0.
We note that this approach is very soft insofar it is not much sensitive to the model (contrarily to the shooting method, described further): it is very easy to add various constraints to the optimal control problem. This is why direct methods are often privileged when the model is not completely fixed. The resulting numerical simulations often give an interesting feedback that may lead to change or adapt the optimal control model under consideration.

Practical numerical implementation

The numerical implementation of such a nonlinear programming problem is standard and can be done in a number of ways, for instance using a penalty method or a SQP method (sequential quadratic programming) or dual methods (like Uzawa). A survey on the use of direct methods in optimal control, with a special interest to applications in aerospace, can be found in [START_REF] Betts | Practical methods for optimal control using nonlinear programming, volume 36 of Advances in Design and Control[END_REF].

From the point of view of practical implementation, in the last years much progress has been done in the direction of combining automatic differentiation softwares (such as the modeling language AMPL, see [START_REF] Fourer | AMPL: A modeling language for mathematical programming[END_REF], or CasADi, see [6]) with expert optimization routines (such as the opensource package Ipopt, see [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF], providing an interior point optimization algorithm for large-scale differential algebraic systems combined with a filter line-search method). With such tools it has become very simple to implement with only few lines of code difficult (nonacademic) optimal control problems, with success and within a reasonable time of computation. Websites such as NEOS (http://www.neos-server.org/neos/solvers/index.html) propose to launch online such kinds of computation: codes can be written in a modeling language such as [START_REF] Fourer | AMPL: A modeling language for mathematical programming[END_REF] (or others) and can be combined with many optimization routines (specialized either for linear problems, nonlinear, mixed, discrete, etc). The advantage of using NEOS is that one has nothing to install on his own machine, and moreover one can test a large number of possible optimization routines. Note that there exists a large number (open-source or not) of automatic differentiation softwares and of optimization routines, it is however not our aim to provide a list of them. They are easy to find on the web. Note also that AMPL, which is very easy and friendly to use, is however not free of charge (although the licence is not expensive) and that CasADi offers a very good and efficient free alternative, see https://web.casadi.org.

Variants

As alluded above, there exist many possible approaches to discretize an optimal control problem, see [START_REF] Betts | Practical methods for optimal control using nonlinear programming, volume 36 of Advances in Design and Control[END_REF] where the important sparsity issues are also discussed. Among those various approaches, we quote the following.

Collocation methods consist of choosing specific points or nodes on every subinterval of a given subdivision of the time interval. Such methods consist of approximating the trajectories and the controls by polynomials on each subinterval. The collocation conditions state that the derivatives of the approximated state match with the dynamics at each node.

Spectral and pseudospectral methods are another class in which the above nodes are chosen as the zeros of specific polynomials such as Gauss-Legendre or Gauss-Lobatto polynomials. Such polynomials are used as a basis to approximate trajectories and controls in appropriate approximation spaces. Since they share nice orthogonality properties, the collocation conditions turn into constraints that are easily tractable for numerical purposes. We refer the reader to [START_REF] Elnagar | Pseudospectral Chebyshev optimal control of constrained nonlinear dynamical systems[END_REF][START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF][START_REF] Ross | Legendre pseudospectral approximations of optimal control problems[END_REF] and to the references therein for more details.

There exist also some probabilistic approaches, such as the method described in [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI-relaxations[END_REF] which consists of first relaxing the optimal control problem in measure spaces and then of seeking the optimal control as an occupation measure, which is approximated by a finite number of its moments (see [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]). This approach relies on algebraic geometry tools and reduces the optimal control problem to some finite dimensional optimization problem involving linear matrix inequalities (LMI). On this topic involving Sums-of-Square (SoS) considerations, we refer the reader to another chapter of the Handbook, [START_REF] Lasserre | Moments and convex optimization for analysis and control of nonlinear partial differential equations[END_REF].

Remark 2.1. Direct methods consist of first discretizing and then optimizing, i.e., we apply optimality conditions in a second step, to the discretized model, in contrast to indirect methods that are discussed in the next section, which consist of applying first optimality conditions (the Pontryagin maximum principle) and then of discretizing the resulting boundary value problem.

While the latter method clearly falls in the classical Lax scheme, "consistency plus stability imply convergence", there is a serious gap there in the direct approaches: to ensure convergence using the Lax scheme, one would a priori need a uniform (with respect to N ) consistency property, which is not true in general because the optimal control problem is an optimization problem in infinite dimension. Surprisingly simple examples of divergence are provided in [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF]. In this same paper, it is shown that convergence is obtained for "smooth enough" optimal control problems without any final state constraint, discretized with Runge-Kutta methods with positive coefficients. We also refer to [START_REF] Bonnans | Computation of order conditions for symplectic partitioned Runge-Kutta schemes with application to optimal control[END_REF][START_REF] Sanz-Serna | Symplectic Runge-Kutta schemes for adjoint equations, automatic differentiation, optimal control, and more[END_REF] for further comments and considerations on symplectic integrators. Convergence has also been established for classes of Legendre pseudospectral methods (see [START_REF] Elnagar | Pseudospectral Chebyshev optimal control of constrained nonlinear dynamical systems[END_REF][START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF][START_REF] Ross | Legendre pseudospectral approximations of optimal control problems[END_REF]).

Goddard problem by a direct approach

To illustrate the method, we present a treatment of the Goddard case presented Section 1.3 using the Crank-Nicolson scheme. This scheme, which is dual to the midpoint rule (see [START_REF] Hairer | Geometric Numerical Integration[END_REF]), has the advantage of using only gridpoints (contrary to the midpoint scheme). It is very easy to implement and the code below is a basic "do it yourself" direct solver that can easily be adpated to other problems. We use the nice JuMP interface in julia to define the discretization of the problem and call the celebrated interior point solver Ipopt previously mentioned. The symbolic-numeric framework put forward by Julia allows to efficiently and transparently use AD (automatic differentiation / differentiable programming) and sparse linear algebra (which is important for structured constraints stemming from one-step like methods to discretize the dynamics). In the spirit of reproducible research [START_REF] Buckheit | Wavelab and reproducible research[END_REF], the code itself is available and executable online on the gallery of the ct: control toolbox project. 1 Discretizing the dynamics boils down to the following lines (note that nonlinear expressions encoding the right-hand side are first defined):

# Dynamics @NLexpressions(sys, begin # D = Cd v^2 exp(-β(r-1)) D[i = 1:N+1], Cd * v[i]^2 * exp(-β * (r[i] -1.0)) # r'= v dr[i = 1:N+1], v[i] # v' = (Tmax.u-D)/m -1/r^2 dv[i = 1:N+1], (Tmax*u[i]-D[i])/m[i] -1/r[i]^2 # m' = -b.Tmax.u dm[i = 1:N+1], -b*Tmax*u[i] end) # Crank-Nicolson scheme @NLconstraints(sys, begin con_dr[i = 1:N], r[i+1] == r[i] + Δt * (dr[i] + dr[i+1])/2.0 con_dv[i = 1:N], v[i+1] == v[i] + Δt * (dv[i] + dv[i+1])/2.0 con_dm[i = 1:N], m[i+1] == m[i] + Δt * (dm[i] + dm[i+1])/2.0 end)
Also note that the constraints have been labeled so that the corresponding Lagrange multipliers can be retrieved. They are indeed approximations of the costate of the continuous problem and will be used as such to initialize successfully a subsequent shooting method (see Section 3). A strong benefit of direct methods is that state constraints are very easy to handle. We may for instance add a state constraint on the velocity (note that some other constraints have been added, directly when defining the unknowns of the problem, to improve convergence of the solver; what would be a complication for indirect methods is actually an asset here): Figure 1: Goddard problem: result of the direct code (Crank-Nicolson scheme), states and Lagrange multipliers associated with the equality constraints that discretise the dynamics. These multipliers approximate the adjoint states and will used as such to initialise an indirect method. A boundary arc is observed as the state constraint on the velocity is saturated. @variables(sys, begin 0.0 ≤ Δt # time step (unknown as tf is free) r

[1:N+1] ≥ r0 # r 0 ≤ v[1:N+1] ≤ vmax # v mf ≤ m[1:N+1] ≤ m0 # m 0.0 ≤ u[1:N+1] ≤ 1.0 # u end)
Boundary constraints are obviously added in the same way, and for N = 100 gridpoints we get the results given on Figures 1 and2. The resulting optimisation problem solved has about 400 variables (the state is of dimension 3 and the control is scalar) and 300 equality constraints (the discretised dynamics), plus box constraints. One could use a finer grid (up to the price of a larger problem), but the result we obtain turns to be precise enough to trigger convergence of much more accurate solve by shooting in the next section. The bang-singular-constrained-bang structure of the solution has indeed been captured by the direct solve, and this is essentially all we need to resort to indirect methods to complete the computation.

3 Indirect approaches: the shooting method

Pontryagin maximum principle

Consider the optimal control problem (3) with t = 0 as initial time, and with g = 0, i.e., for the moment, there is no state constraint.

Recall that the end-point mapping E : R d × R + × U T,U → R d of the system is defined by E(x 0 , T, u) = x(x 0 , T, u), where t → x(x 0 , t, u) is the trajectory solution of the control system, corresponding to the control u, such that x(x 0 , 0, u) = x 0 . Here, U T,U is the set of controls , possibly of a singular arc interior to the control bounds (u ∈ (0, 1)), of a boundary arc due to the state constraint on the velocity, and of a final bang arc (u 0). This structure information, of combinatorial nature, is obtained here without any a priori knowledge on the solution. It is the key to define an appropriate shooting function and solve very accurately the problem thanks to an indirect method.

taking their values in U and such that the corresponding trajectory is well defined on [0, T ] (no blow-up).

The set U T,R r , endowed with the standard topology of L ∞ ([0, T ], R r ), is open, and the endpoint mapping is C 1 on U T,R r (and C k whenever the dynamics are C k ).

Denoting by C the cost functional, the optimal control problem can be written as the infinite dimensional minimization problem of minimizing C over the considered set of controls, under the constraint g f (E(x 0 , T, u)) ≤ 0 on the final state.

Assume for one moment that we are in the simplified situation where, in the optimal control problem, the initial point x 0 is fixed, the final point x 1 is fixed, the final time T is fixed, U = R r and there are no state constraints along the trajectory. Then, the optimal control problem is exactly equivalent to the optimization problem min

E(x 0 ,T,u)=x 1 C(T, u). ( 12
)
If u is optimal then there must exist Lagrange multipliers (ψ,

ψ 0 ) ∈ R d × R \ {0} such that ψ.dE x 0 ,T (u) = -ψ 0 dC T (u). (13)
Defining the Lagrangian L T (u, ψ, ψ 0 ) = ψE x 0 ,T (u) + ψ 0 C T (u), the first-order condition (13) is written as

∂L T ∂u (u, ψ, ψ 0 ) = 0. (14)
The Pontryagin maximum principle, recalled hereafter, is a far-reaching statement of these conditions (see [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF], see also [START_REF] Lee | Foundations of optimal control theory[END_REF]).

We recall here the statement without any state constraint (g = 0).

almost everywhere on [0, T ]. If moreover the final time T is not fixed, then

max v∈U H(T, x(T ), p(T ), p 0 , v) = 0. ( 16
)
Additionally, if g f (x(T )) = 0 then the adjoint vector can be built in order to satisfy the transversality condition p(T ) -p 0 ∇ϕ(x(T )) ⊥ ker dg f (x(T )).

If only some components of g f (x(T )) are zero then the orthogonality relation [START_REF] Betts | Practical methods for optimal control using nonlinear programming, volume 36 of Advances in Design and Control[END_REF] holds for the corresponding components.

Actually, compared with (13) (at least, in the simplified case), the adjoint vector of the Pontryagin maximum principle can be constructed so that (ψ, ψ 0 ) = (p(T ), p 0 ) up to scaling. In particular, the Lagrange multiplier ψ is unique (up to a multiplicative scalar) if and only if the trajectory x(•) admits a unique extremal lift (up to scaling).

When p 0 < 0, the extremal is said to be normal, and in this case, since the Lagrange multiplier is defined up to scaling, it is usual to normalize it so that p 0 = -1. When p 0 = 0, the extremal is said to be abnormal.

In the absence of control constraints, i.e., when U = R r , abnormal extremals project exactly onto singular trajectories. Recall that a trajectory (x, u) defined on [0, T ] is said to be singular when the linearized control system along it is not controllable in time T ; equivalently, in terms of the end-point mapping, the Fréchet differential dE x 0 ,T is not surjective.

In the normal case, ψ = p(T ) coincides up to scaling with the gradient of the value function (solution of the Hamilton-Jacobi equation); see, e.g., [START_REF] Clarke | The relationship between the maximum principle and dynamic programming[END_REF] for precise results. Remark 3.1. In some sense, generically the solution of the optimal control problem is unique, and moreover it has a unique extremal lift. This well known fact is related to the differentiability properties of the value function (see for instance [START_REF] Aubin | Set-Valued Analysis[END_REF][START_REF] Clarke | The relationship between the maximum principle and dynamic programming[END_REF], and see [START_REF] Cannarsa | Semiconcave functions, Hamilton-Jacobi equations, and optimal control[END_REF][START_REF] Rifford | On the stabilization problem for nonholonomic distributions[END_REF][START_REF] Rifford | Morse-Sard type results in sub-Riemannian geometry[END_REF][START_REF] Stefani | Regularity properties of the minimum-time map[END_REF] for results on the size of the set where the value function is differentiable). Remark 3.2. When applying the Pontryagin Maximum Principle, we must distinguish between two extremal flows: the normal flow, with p 0 < 0 (and in general we normalize to p 0 = -1, as said above), and the abnormal one, for which p 0 = 0. The scalar p 0 is a Lagrange multiplier associated with the instantaneous cost.

In many situations, where some qualification conditions hold, abnormal extremals do not exist in the problem under consideration, but in general it is difficult to guess whether, given some initial and final conditions, these qualification conditions hold or not.

Remark 3.3. The fact that p 0 ≤ 0 is a convention due to Pontryagin, which leads to the maximum principle. The choice p 0 ≥ 0 would lead to a minimum condition, instead.

The component p 0 appears in the transversality condition [START_REF] Betts | Practical methods for optimal control using nonlinear programming, volume 36 of Advances in Design and Control[END_REF], even when we consider a cost with = 0. Note that, if the final point x(T ) is let free (i.e., there is no g f ) then this condition leads to p(T ) = p 0 ∇ϕ(x(T )), and then necessarily p 0 = 0 and we can normalize to p 0 = -1.

Remark 3.4. When there are some state constraints g(x(t)) ≤ 0 along the trajectory, the Pontryagin maximum principle is modified as follows. We keep the same definition for the Hamiltonian H.

If x(•) is optimal then there must exist p 0 ≤ 0, an absolutely continuous adjoint vector p(•) and a nonnegative Radon measure µ, the triple (p, p 0 , µ) being nontrivial, such that the adjoint equation is

ṗ(t) = - ∂H ∂x (t, x(t), q(t), p 0 , u(t))
almost everywhere, with

q(t) = p(t) + [0,t) ∂g ∂x (x(s)) dµ(s),
(the interval being closed when t = T ) and the maximization condition becomes

H(t, x(t), q(t), p 0 , u(t)) = max v∈U H(t, x(t), q(t), p 0 , v)
almost everywhere. Finally, in the transversality condition on the final adjoint vector, one replaces p(T ) with q(T ).

Note that, taking q(t) as a new adjoint, we have q = ṗ+ ∂g ∂x (x(t)) dµ(t), and since ṗ = -∂H ∂x (q), we have q = -q ∂f ∂x -p 0 ∂f 0 ∂x + ν ∂g ∂x with ν ≥ 0, which is the formulation that one can also find in the existing literature (possibly with an opposite sign for the state constraint, and obvious changes in the previous expressions).

Shooting method

To compute optimal trajectories thanks to the Pontryagin Maximum Principle, the first step is to make explicit the maximization condition, at least, when this is possible (otherwise this can be done numerically). A usual assumption, often satisfied, is the strict Legendre assumption, stipulating that ∂ 2 H ∂u 2 (t, x, p, p 0 , u) is negative definite. Under this assumption, an implicit function argument gives, locally, a control u expressed in function of the state and of the adjoint. Let us assume, below, that we are in the normal case, p 0 = -1. Plugging the resulting expression of the control in the Hamiltonian equations, and defining the reduced (normal) Hamiltonian by H r (t, x, p) = H(t, x, p, -1, u(x, p)), we obtain that every normal extremal is solution of the reduced Hamiltonian system

ẋ(t) = ∂H r ∂p (t, x(t), p(t)), ṗ(t) = - ∂H r ∂x (t, x(t), p(t)). (18)
Definition 3.1. Denoting by (x(t, x 0 , p 0 ), p(t, x 0 , p 0 )) the solution of (18) starting at (x 0 , p 0 ) for t = 0, the exponential mapping is defined by exp x 0 (t, p 0 ) = x(t, x 0 , p 0 ).

The exponential mapping parametrizes the normal extremal flow. The abnormal extremal flow can be parametrized as well provided that an appropriate Legendre assumption holds in the abnormal case. Remark 3.5. Let us give an example where the Hessian of the Hamiltonian is degenerate: the minimal time problem for single-input control affine systems ẋ(t) = F 0 (x(t)) + u(t)F 1 (x(t)) without constraint on controls. In that case, the maximization condition implies that p(t), F 1 (x(t)) = 0 along the extremals, and the optimal control u is singular. To compute u, we perform two successive derivations in time of the latter relation, obtaining

p(t), [F 0 , F 1 ](x(t)) = 0 and p(t), [F 0 , [F 0 , F 1 ]](x(t)) + u(t) p(t), [F 1 , [F 0 , F 1 ]](x(t)) = 0, where [•, •]
is the Lie bracket of vector fields. We thus obtain u in function of x and p, provided that p(t),

[F 1 , [F 0 , F 1 ]](x(t)) = 0,
which is the so-called strong generalized Legendre-Clebsch condition (see, e.g., [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]). Actually, under generic conditions on the vector fields, the above computation can always be performed (see [START_REF] Chitour | Genericity results for singular curves[END_REF][START_REF] Chitour | Singular trajectories of control-affine systems[END_REF]). See Goddard example in this section for an example of this computation.

Remark 3.6. When an abnormal flow can be well defined, we then have to deal with two extremal flows (and two exponential mappings). In general, however, the abnormal flow "does not fill much space". For example, in [START_REF] Agrachev | Any sub-Riemannian metric has points of smoothness[END_REF][START_REF] Rifford | On the stabilization problem for nonholonomic distributions[END_REF] it is proved that for control-affine systems without drift (satisfying the Hörmander assumption), with quadratic cost, the image of the abnormal exponential mapping has an empty interior in the state space, and is even of Lebesgue measure zero under some additional assumptions.

Remark 3.7. The Pontryagin Maximum Principle is a first-order necessary condition for optimality, asserting that if a trajectory is optimal then it should be sought among projections of extremals joining the initial point to the final target. This yields the shooting method that is described hereafter. But, before coming to the description of that method, it is interesting to observe that, conversely, the projection of a given extremal may not be (locally or globally) optimal. This important observation has led to develop second-order optimality conditions in optimal control. In terms of the Lagrangian (in the simplified setting), considering the intrinsic second order derivative Q T of the Lagrangian, given by

Q T = ∂ 2 L T ∂ 2 u (u, ψ, ψ 0 ) | ker ∂L T ∂u ,
a second-order necessary condition for optimality is that Q T be nonpositive, and a second-order sufficient condition for local optimality is that Q T be negative definite. Such conditions admit a number of generalizations for optimal control problems involving control and/or state constraints. It happens that, given a fixed extremal starting at (x 0 , p 0 ), in the simplified context and under appropriate generic assumptions, the quadratic form Q T is not degenerate (i.e., its kernel is trivial) if and only if the exponential mapping exp x 0 (t c , •) is not an immersion at p 0 (that is, its differential is not injective). This result, coming from symplectic considerations and generalizing the Riccati theory (see [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF][START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF]), yields to algorithms for computing the first conjugate time along a given extremal (see [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF] for a survey). By definition, the first conjugate time along an extremal is the first time t c at which the quadratic form Q tc has a nontrivial kernel. This means that the trajectory x(•) under consideration is locally optimal (in L ∞ topology) on [0, t] if and only if t < t c . Computing a first conjugate time amounts to computing the first time at which some determinant along the extremal vanishes. More generally, the fact that the exponential mapping is not an immersion can be translated in terms of Jacobi fields, like in Riemannian geometry.

Let us now describe the contents of the shooting method. The Pontryagin Maximum Principle states that every optimal trajectory is the projection of an extremal. After having made explicit the maximization condition, the problem is reduced (for instance, in the normal case) to an extremal system of the form ż(t) = F (t, z(t)), where z(t) = (x(t), p(t)), and initial, final, transversality conditions, are written as R(z(0), z(T )) = 0. We thus end up with a boundary value problem (BVP) of the form

ż(t) = F (t, z(t)), R(z(0), z(T )) = 0. ( 19
)
Denote by z(t, z 0 ) the solution of the Cauchy problem ż(t) = F (t, z(t)), z(0) = z 0 , and set

G(z 0 ) = R(z 0 , z(T, z 0 )).
The boundary value problem ( 19) is then equivalent to solving G(z 0 ) = 0, i.e., to finding a zero of the function G. By definition, the (single) shooting method consists of implementing a Newton like method to find a zero of G.

The feasibility of the shooting method relies on the fact that the Jacobian of the mapping G is nonzero. According to the above remark, in the simplified case, this determinant is nonzero, i.e., the (single) shooting method is well-posed at time T , locally around p 0 , if and only if the exponential mapping exp x 0 (t, •) is an immersion at p 0 , that is, if and only if T is not a conjugate time.

Although this result admits generalizations to a number of contexts (see, e.g., [START_REF] Bonnans | Well-posedness of the shooting algorithm for state constrained optimal control problems with a single constraint and control[END_REF][START_REF] Bonnans | Second-order analysis for optimal control problems with pure state constraints and mixed control-state constraints[END_REF]), there still misses a complete conjugate time theory involving state and control constraints, in which the trajectories may have bang, singular, boundary arcs.

The single shooting method can be refined to the multiple shooting method, in which one may add a number of intermediate nodes, thus incorporating new (matching) conditions in the shooting function G. This can be useful to face with numerical instability issues for instance, or to implement the shooting method in bang-bang situations where one knows in advance the number of switchings.

Efficient shooting methods are available in the HamPath package2 , now encapsulated in the Python package nutopy of the ct: control toolbox project. 3 These open-source packages also contain implementations of conjugate time computations and of several homotopy routines that are particularly useful in a number of contexts (see [START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF] for a survey on the use of continuation methods in optimal control). Remark 3.8. Numerically, the shooting method is the combination of a numerical integration of a differential equation with a Newton method for finding a zero of a map (the shooting function). It thus inherits of the main features of a Newton method: when it converges, the convergence is extremely quick and the result is very accurate. However, it may be difficult to initialize successfully: finding a good initial guess for z 0 , in the above notations, may be a real challenge. To face with this difficulty, several possible remedies are known, such as the following, surveyed in [START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF] (see also the references therein):

• Since direct methods are less sensitive to the initialization, it is often successful to first run a direct approach (even with a quite rough mesh) so that, if can obtain its convergence, then the corresponding optimal solution and Lagrange multiplier can be used as an approximation of the searched extremal.

• Continuation and homotopy methods can be combined with the shooting method (and also, by the way, with direct methods): when a given problem happens to be difficult to solve, or quite ill-posed, because of some too restricted parameters or because of too constraining terms in the dynamics, one can try to relax the optimal control problem by adding some continuation parameters in front of those terms, then run a series of shooting methods with the continuation parameters ranging iteratively from 0 to 1 (with adaptive steps if necessary).

• Geometric control gives useful information on the local or global structure of the optimal controls, for instance one can guess in advance the number of switchings in a bang-bang strategy under some appropriate assumptions. This knowledge can then be combined with the Pontryagin maximum principle in order to drastically reduce the complexity of the shooting problem.

Much more could be said on these classical issues, but since they have already been surveyed in [START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF] we do not elaborate more. We next describe another powerful remedy that has emerged recently, although it relies on an old concept.

Turnpike property

Assume that f and do not depend on time. In few words, the turnpike property stipulates that, for optimal control problems in large time, under slight assumptions it is expected that the optimal solution, the optimal control and the associated adjoint remain essentially close to static values, except at the beginning and at the end of the time interval. Moreover these static values correspond to the optimal solution of a static optimization problem.

The idea is very easy to understand. When the time T is large, setting ε = 1/T and τ = t/T = εT , the optimal control problem consists of determining a trajectory y(τ ) = x(t) and a control v(τ ) = u(t) solution of εy (τ ) = f (y(τ ), v(τ )) and minimizing the cost 1 0 (y(τ ), v(τ )) dτ , under various constraints. At the formal level, we see that, when ε → 0, at the limit we find the static optimization problem consisting of minimizing (y, v) under the constraint f (y, v) = 0, i.e., under the constraint of being an equilibrium of the controlled dynamics.

To give a more precise insight, let us establish the so-called exponential turnpike phenomenon in the linear-quadratic case. Let x 0 , x 1 , x ∈ R d and û ∈ R r be fixed. We consider the optimal control problem in fixed final time T > 0:

ẋ(t) = Ax(t) + Bu(t) x(0) = x 0 , x(T ) = x 1 min T 0 (x(t) -x) Q(x(t) -x) + (u(t) -û) U (u(t) -û) dt (20)
where x(t) ∈ R d and u(t) ∈ R r , and where Q et U are real-valued symmetric positive definite matrices. By strict convexity, there exists a unique optimal solution (x T (•), u T (•)) of [START_REF] Bokanowski | An efficient filtered scheme for some first order timedependent hamilton-jacobi equations[END_REF]. We assume that the pair (A, B) satisfies the Kalman condition.

The Hamiltonian of the problem is

H(x, p, p 0 , u) = p, Ax + p, Bu + p 0 ((x -x) Q(x -x) + (u -û) U (u -û)).
Let us prove that p 0 = 0. By contradiction, if p 0 = 0 then the condition ∂H ∂u = 0 yields p(t), B = 0, and by successive derivations and using the fact that ṗ(t) = -A p(t), we obtain p(t), A k B = 0, which raises a contradiction with the Kalman condition since p(t) = 0.

We choose then to normalize the adjoint so that p 0 = -1/2. The condition ∂H ∂u = 0 yields u T (t) = û + U -1 B p T (t), and the extremal system is

ẋT (t) = Ax T (t) + BU -1 B p T (t) + B û ṗT (t) = Qx T (t) -A p T (t) -Qx i.e., d dt x T (t) p T (t) = M x T (t) p T (t) + B û -Qx where M = A BU -1 B Q -A
Besides, the static optimization problem is

min (x,u)∈R d ×R r Ax+Bu=0 (x -x) Q(x -x) + (u -û) U (u -û) .
This strictly convex problem has a unique solution (x, ū), associated with a normal Lagrange multiplier (p, -1) (the problem is qualified because ker(A ) ∩ ker(B ) = {0}, as a consequence of the Kalman condition). According to the Lagrange multiplier rule, which is here a necessary and sufficient condition for optimality, there exists p ∈ R d \ {0} such that ū = û + U -1 B * p and

Ax + BU -1 B * p + B û = 0 Qx -A * p -Qx = 0 (21)
i.e.,

M x p + B û -Qx = 0 0
We have the following exponential turnpike property.

Proposition 3.1. There exist constants C > 0 and ν > 0, not depending on T , such that

x T (t) -x + u T (t) -ū + p T (t) -p ≤ C(e -νt + e -ν(T -t) ) ∀t ∈ [0, T ]. (22)
Proof. Using the above optimality systems, we have

d dt x T (t) - x p T (t) -p = M x T (t) - x p T (t) -p
In order to prove [START_REF] Bokanowski | An adaptive sparse grid semi-lagrangian scheme for first order Hamilton-Jacobi Bellman equations[END_REF], the crucial observation is that the matrix M , which is Hamiltonian, is hyperbolic, i.e., all its eigenvalues have a nonzero real part (actually, the number of unstable modes is equal to the number of stable modes). To prove this hyperbolicity property, we start by noting that, as a consequence of the Kalman condition on (A, B), we have

ker(A -pI n ) ∩ ker(B ) = {0} ∀p ∈ C. ( 23
)
It follows that the matrix M has no purely imaginary eigenvalue. Indeed, let

(z 1 , z 2 ) ∈ C n × C n and let µ ∈ R be such that (M -iµ) z 1 z 2 = 0. Then (A -iµ)z 1 + BU -1 B z 2 = 0 Qz 1 -(A + iµ)z 2 = 0 hence z 1 = Q -1 (A + iµ)z 2 and thus (A -iµ)Q -1 (A + iµ)z 1 + BU -1 B z 2 = 0. Multiplying to the left by z 2 , we obtain Q -1/2 (A + iµ)z 2 2 + U -1/2 B z 2 2
= 0 and hence (A + iµ)z 2 = 0 and B z 2 = 0. We infer that z 2 = 0 by using [START_REF] Bonnans | Well-posedness of the shooting algorithm for state constrained optimal control problems with a single constraint and control[END_REF]. The claim is proved.

Since M is hyperbolic, there exists a real-valued square matrix P of size 2n such that

P -1 M P = M 1 0 0 M 2
where all eigenvalues of the matrix M 1 have a negative real part, and all eigenvalues of M 2 have a positive real part. Now, setting

x T - x p T -p = P v w , the extremal system gives d dt v(t) w(t) = P -1 M P v(t) w(t) hence v(t) = M 1 v(t) and w(t) = M 2 w(t).
All eigenvalues of M 1 have a negative real part, hence there exist constants C 1 > 0 and ν 1 > 0, not depending on T , such that v(t) ≤ C 1 e -ν 1 t . For w(t), we reverse time and we apply the same argument, hence w(t) ≤ C 2 e -ν 2 (T -t) . Finally, [START_REF] Bokanowski | An adaptive sparse grid semi-lagrangian scheme for first order Hamilton-Jacobi Bellman equations[END_REF] follows by noting that x T (t) et p T (t) are linear combinations of v(t) and w(t).

The exponential turnpike property [START_REF] Bokanowski | An adaptive sparse grid semi-lagrangian scheme for first order Hamilton-Jacobi Bellman equations[END_REF] says that, except near t = 0 and t = T , the optimal state, control and adjoint are exponentially close to static values, themselves corresponding to the solution of the associated static optimization problem.

In [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF], the proof of the above proposition is a bit different and relies on the use of the Riccati theory: actually, the matrix P is built by considering the minimal and maximal solutions of the Riccati algebraic equation, which gives an interpretation of the constants C and ν in function of these matrices. Anyway, the argument remains very easy and it withstands a number of generalizations, to infinite dimension, to nonlinear dynamics, also to situations where the turnpike is not restricted, as above, to a singleton but may even consist of a nontrivial set of trajectories (an example being the periodic turnpike).

We refer to the chapter [START_REF] Faulwasser | Turnpike properties in optimal control: An overview of discrete-time and continuous-time results[END_REF] of the Volume 1 of the present Handbook, for a recent survey on the turnpike property in optimal control, containing a number of references and commenting also on the important related notion of dissipativity.

What we want to point out here is that, when the turnpike property is satisfied, we can use it to successfully initialize a variant of the shooting method. Remark 3.9 (Variant of the shooting method). In the turnpike context, we know that, in the middle of the trajectory, x T (T /2) and p T (T /2) are exponentially close to x and p. In such conditions, if this is feasible, it is relevant to first compute the solution of the static problem and then to implement a variant of the shooting method by initializing it "at the middle", as follows. Using the notations of the previous section, the unknown is now z(T /2) ∈ R 2d , which is initialized at (x, p), the steady-state solution of the static optimal control problem. Then, integrating backwards the extremal system over [0, T /2], we compute z(0); integrating forward the extremal system over [T /2, T ], we compute z(T ). Finally, the unknown z(T /2) is tuned so that G(z(0), z(T )) = 0, thanks to a Newton method. It has been observed in [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF] that this variant of the single shooting method is very efficient.

Solving the Zermelo problem by the shooting method

We consider the navigation problem of Zermelo presented in Section 1.3, with a = 0 and b = 1. We take

V max = 1, y f = (20, 1), h(y 2 ) = 3 + 0.2y 2 (1 -y 2 )
, where the target is the point y f (a slightly simpler problem than the general previously presented), while the initial point is y 0 = (0, 0). We treat the minimum time case so that the velocity can be set to V max and the only control is the angle. An obstacle is inserted along the unconstrained optimal trajectory:

(y 1 -y 1,f /2) 2 a 2 1 + (y 2 -y 2,f /2.5) 2 a 2 2 ≤ 1
with a 1 = 2 and a 2 = 0.1. We use a logarithmic barrier to penalize internally the state constraint and we consider the augmented cost

T -α 1 0 ln (y 1 (t) -y 1,f /2) 2 a 2 1 + (y 2 (t) -y 2,f /2.5) 2 a 2 2 -1 dt → min .
For a detailed study of internal penalization, we refer to [START_REF] Malisani | An interior penalty method for optimal control problems with state and input constraints of nonlinear systems[END_REF]. For the computation below, α is set to 1e -3. (A finer approach would be to perform a continuation on α.) We apply the maximum principle which leads to the (normal) maximized Hamiltonian

H(y, p) = V max p 2 1 + p 2 2 + p 1 h(y 2 ) + α ln (y 1 -y 1,f /2) 2 a 2 1 + (y 2 -y 2,f /2.5) 2 a 2 2 -1 -1. ( 24
)
The shooting problem then consists in finding the initial value p(0) ∈ R 2 so that, integrating the flow of the maximized Hamiltonian, the target y f is reached. Moreover, the equation H = 0 is added to accomodate the fact that the final time T is free. Depending on the initialisation of the shooting method, two solutions are obtained, see Figure 3. The second one is clearly a local minimizer. An interesting approach in such a situation is to rely on an HJB solver (see Section 4) to retrieve an initial guess of the adjoint state (and of the final time) that avoid strict local minima. It can indeed be checked numerically that the HJB solution on the previous data allows shooting to converge towards the global minimizer on the internally penalized problem.

Remark 3.10. A refinement of this computation could involve a continuation (or differential homotopy) on the size of the obstacle. Using path following methods, one expects to be able to track the two branches associated with the global and strictly local minimum, then to decide which one is optimal for each size of the obstacle by comparing the associated final times.

Another relevant observation on this problem is related to the turnpike property described in Section 3.3. Although, the problem is not linear-quadratic, one can easily guess the role played by this property for "large" (in terms of y 1 ) final conditions. We keep the same data for V max and for the current, remove the obstacle, and now target y f = (200, 1). For such a distant target, one expects the optimal control to use the "fast lane" so that, most of the time, y 2 remains close to 1/2 where the current h is maximum. So the guess would be that, for a large part of the trajectory, y 2 1/2, u 2 = ẏ2 0, u 1 1, and T (y 1,f -y 1,0 )/(1 + h(1/2)). The maximized (normal) Hamiltonian is (compare [START_REF] Bonnans | Second-order analysis for optimal control problems with pure state constraints and mixed control-state constraints[END_REF])

H(y, p) = V max p 2 1 + p 2 2 + p 1 h(y 2 ) -1,
so that, on H = 0, p (1/(1 + h(1/2, 0)). These approximations can be used to initialise the variant of the shooting method described in Remark 3.9. It is straightforwardly implemented in Julia according to4 The structure of the optimal solution is as expected, see Figure 4.

Remark 3.11. Further analysis is required to tackle the case of "strong" currents. When there exist zones where the drift cannot be compensated by the control, abnormal trajectories come into play and discontinuities of the value function (minimum time) are observed. See [START_REF] Bonnard | Abnormal geodesics in 2d-zermelo navigation problems in the case of revolution and the fan shape of the small time balls[END_REF] for a detailed treatment of such cases.

Solving the Goddard problem by the shooting method

We chain the result of Section 2.4, obtained by a direct method, with a shooting method. The previous numerical solution, although not very accurate, has captured the structure of the solution (more precisely of what one may hope to be a local minimiser, at least). This knowledge allows us to define the appropriate shooting function, namely one that combines four arcs (bangsingular-constrained-bang), each one being the flow of a relevant Hamiltonian. On can then leverage the accurate knowledge (including the Hamiltonian character) gained on each subarc by means of Pontryagin maximum principle to obtain a very accurate solve. In general, knowing the structure alone is not enough to actually solve the problem. One also needs a good initial guess for the zero of the shooting function, which turns to be also provided by the previous direct solve. In order to set up our shooting, we rely on the maximum principle and observe that the control can be either bang, singular of boundary. Indeed, system (7)-( 9) is affine in the control and can be written according to

ẋ(t) = F 0 (x(t)) + u(t)F 1 (x(t)), u(t) ∈ [0, 1],
with x = (r, v, m) ∈ R 3 and vector fields that we shall define in Julia, completing the code started Section 2.4:

# Dynamics function F0(x) r, v, m = x D = Cd * v^2 * exp(-β*(r-1.0)) F = [ v, -D/m-1.0/r^2, 0.0 ] return F end function F1(x) r, v, m = x F = [ 0.0, Tmax/m, -b*Tmax ] return F end
In order to deal with the state constraint g(x(t)) := v max -v(t) ≥ 0, we follow Remark 3.4 (note that we use the opposite sign for the constraint) and introduce the Hamiltonian (where the constraint has been directly adjoined, see [START_REF] Hartl | A survey of the maximum principles for optimal control problems with state constraints[END_REF])

H(x, p, u, µ) = H 0 (x, p) + uH 1 (x, p) + µg(x)
where H i (x, p) := p, F i (x) are the Hamiltonian lifts of the aforementioned fields. The maximisation condition tells us that u is bang (0 or 1) whenever H 1 is not zero. Besides, whenever the state constraint is not active, the associated non-negative multiplier µ vanishes (complementarity condition). As a result, bang arcs are obtained by computing the flow of either H 0 (case H 1 < 0, u = 0) or H 0 + H 1 (case H 1 > 0, u = 1). Conversely, when H 1 vanishes identically, assuming the state constraint is not active, as indicated in Remark 3.5 one can differentiate a.e. two times H 1 to retrieve the singular control provided the length three Poisson bracket H 101 := {H 1 , {H 0 , H 1 }} is not zero (singular of order one; see, e.g., [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]): whose flow coincides with the extremal flow on Σ := {H 1 = H 01 = 0} (see, e.g., [START_REF] Agrachev | Control theory from the geometric viewpoint[END_REF]). Along a boundary arc where the state constraint is activated, if the control is interior (u(t) ∈ (0, 1)), H 1 must also vanish. Moreover, if the constraint is of order one (that is if the control appears when the equality g(x(t)) = 0 is differentiated once) which gives in this case (Lie derivative of the constraint along F 1 ) F 1 • g = 0, the interiority of the boundary control implies that the contact with the constraint is transverse at the exit time t 3 (where, in our case, a bang arc u = 0 is joined): ġ(t 3 +) = 0, and there is no jump on the adjoint [START_REF] Hartl | A survey of the maximum principles for optimal control problems with state constraints[END_REF]. Denoting t 2 the entry point, on (t 1 , t 2 ) one retrieves the boundary control by differentiating once g(x(t)) = 0 as

u s (x, p) = -
0 = ġ = F 0 •g +uF 1 •g, which implies that u b (x) = -F 0 •g F 1 •
g under the two previous assumptions. Similarly, differentiating once H 1 = 0 allows to compute the multiplier µ as 0

= Ḣ1 = {H 0 + uH 1 + µg, H 1 } = H 01 -µF 1 • g, so that µ b (x, p) = H 01 F 1 •g .
The relevant flow is an integral curve of the boundary Hamiltonian

H b (x, p) := H(x, p, u b (x), µ b (x, p)).
All in all, the symbolic-numeric framework allows to define everything in terms of the vector fields and of the constraint: It is straightforward to retrieve initial guesses for these unknowns from the direct code solution, most notably by retrieving the Lagrange multipliers from the optimizer through JuMP interface (note the minus sign to accomodate the convention on the adjoint state vs. the one on Lagrange multipliers):

#
p = -[ [ dual(con_dr[i]), dual(con_dv[i]), dual(con_dm[i]) ] for i in 1:N ]
Moreover, in view of these results, the constraint on the final mass is assumed to be active (which is in accordance with the final zero bang arc). Automatic differentiation can also be used to compute the derivative of the shooting function (several Julia backends to do this are available and include differentiating properly calls to ODE solvers). See Figures 5 and6 for the associated numerical simulations, while the code itself is available and executable online.

Hamilton-Jacobi-Bellman approach

The key tool for the Hamilton-Jacobi-Bellman approach (HJB in what follows) is the value function V : [0, T ] × R d -→ R, which associates to any initial condition the optimal value of the control problem, and is defined by

V(t, x) := inf ϕ(x(T )) + T t (t, x(t), u(t))dt (x, u) ∈ X [t,T ] (x), g(x(τ )) ≤ 0 for all τ ∈ [t, T ], and g f (x(T )) ≤ 0 . (25)
It is known that the value function V can be characterized as the unique solution, in a suitably weak sense, of a Hamilton-Jacobi type equation [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]. Starting from the knowledge of the value function, which is typically obtained via numerical approximation, it is possible to reconstruct the optimal solution in feedback form, i.e., with an optimal control given as a function of the current state.

Unconstrained Bolza control problems

In this section we present first some classical results of HJB approach when the optimal control problems is free of state constraints (i.e., g = g f ≡ 0). The control problem is

V(t, x) = inf ϕ(x(T )) + T t (t, x(t), u(t))dt (x, u) ∈ X [t,T ] (x) . (26)
When (H 4 ) is satisfied, then the control problem admits a solution. The Gronwall estimate on the trajectories and the Lipschitz regularity of the cost functions ensure that the value function is Lipschitz continuous.

Proposition 4.1. Under (H 0 )-(H 1a )-(H 1b )-(H 2 )-(H 3 ), the value function V is Lipschitz con- tinuous.
The systematic study of optimal control problems dates from the late 1950s, and one important method is the Dynamical Programming and HJB approach. This approach reduces the study to investigating the analytic solution to a partial differential equation of the HJB type. However, the problems are usually nonlinear and the analytic solutions do not exist. To deal with the lack of smoothness, two important tools have been developed: the theory of viscosity solutions and the non-smooth analysis. The theory of viscosity solutions for nonlinear Hamilton-Jacobi equations, introduced in the early 1980s by Crandall-Lions [START_REF] Crandall | Condition d'unicité pour les solutions generalisées des équations de Hamilton-Jacobi du premier ordre[END_REF][START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF] and Crandall-Evans-Lions [START_REF] Crandall | Some properties of viscosity solutions of Hamilton-Jacobi equations[END_REF], allows to define generalized solutions to broad classes of nonlinear partial differential equations, including the HJB equations of optimal control problems. We refer also to the books [START_REF] Barles | Solutions de viscosité des équations de hamilton-jacobi[END_REF][START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF] for a more complete introduction to this theory. Another important tool is the non-smooth analysis, which addresses to differential analysis for non-smooth functions. We refer the reader to [START_REF] Aubin | Differential inclusions. Set-valued maps and viability theory[END_REF][START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF][START_REF] Clarke | Nonsmooth analysis and control theory[END_REF][START_REF] Vinter | Optimal control[END_REF] for an introduction of the theory and its applications.

Dynamic programming and Hamilton-Jacobi-Bellman equation

The fundamental idea of Dynamic Programming is that the value function V satisfies a functional equation, often called the Dynamic Programming Principle (DPP). Proposition 4.2 (Dynamic Programming Principle). Assume (H 0 )-(H 3 ), and denote by x u (s) the solution of (1) for a given control u, and such that x u (t) = x. Then, for all x ∈ R d , t ∈ [0, T ] and h ∈ [0, T -t], the value function V satisfies the equality

V(t, x) = inf u∈U t+h t (s, x u (s), u(s))ds + V(t + h, x u (t + h)) . (27)
This principle provides two properties, called sub-and-super-optimality, which are defined as follows. For any function v : [0, T ] × R d → R, 1. we say that v satisfies the super-optimality principle if for any t ∈ [0, T ], x ∈ R d , there exists

(x u , u) ∈ X [t,T ] (x) such that v(t, x) ≥ v(t + h, x u (t + h)) + t+h t (s, x u (s), u(s)) ds, ∀ h ∈ [0, T -t];
2. we say that v satisfies the sub-optimality principle if for any t ∈ [0, T ], x ∈ R d , and

(x u , u) ∈ X [t,T ] (x), v(t, x) ≤ v(t + h, x u (t + h)) + t+h t (s, x u (s), u(s)) ds, ∀ h ∈ [0, T -t].
In principle, once chosen a "small" time increment h, the DPP allows to compute the value function at the point (t, x) by splitting the trajectories at time t+h and starting with the position of the trajectory y t,x at time t+h. As it will be seen later on, it is possible to construct numerical schemes based on this idea to compute an approximation of the value function.

If the function V is differentiable, we can derive from the DPP its infinitesimal version, the Hamilton-Jacobi-Bellman (HJB) equation

-V t (t, x) + H(t, x, DV(t, x)) = 0 for (t, x) ∈ (0, T ) × R d , (28a)
V(T, x) = ϕ(x) for x ∈ R d , (28b)
where the Hamiltonian is given by

H(t, x, q) = sup u∈U -f (t, x, u) • q -(t, x, u) . (29)
In general, V is Lipschitz continuous but not differentiable. Furthermore, the nonlinear equation ( 28) is not expected to admit a classical solution. To deal with a value function lacking smoothness, as mentioned before, the theory of viscosity solutions and the non-smooth analysis will be applied.

Theory of viscosity solutions. Hamilton-Jacobi-Bellman equations

First, we recall the definition of viscosity solution for HJB equations (see [START_REF] Crandall | Condition d'unicité pour les solutions generalisées des équations de Hamilton-Jacobi du premier ordre[END_REF][START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF][START_REF] Barles | Solutions de viscosité des équations de hamilton-jacobi[END_REF][START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]). (i) We say that V is a viscosity supersolution if V is lower semicontinuous (lsc) and for any

φ ∈ C 1 ((0, T ) × R d ) and (t 0 , x 0 ) ∈ (0, T ) × R d local minimum point of V -φ, we have -φ t (t 0 , x 0 ) + H(t 0 , x 0 , Dφ(t 0 , x 0 )) ≥ 0.
(ii) We say that V is a viscosity subsolution if V is upper semicontinuous (usc) and for any φ ∈ C 1 ((0, T ) × R d ) and (t 0 , x 0 ) ∈ (0, T ) × R d local maximum point of V -φ, we have -φ t (t 0 , x 0 ) + H(t 0 , x 0 , Dφ(t 0 , x 0 )) ≤ 0.

(iii) We say that V is a viscosity solution if it is both a viscosity supersolution and a viscosity subsolution and the final condition is satisfied:

V (T, x) = ϕ(x) in R d .
Remark 4.1. There are also some equivalent definitions which are more local using the superand sub-differentials, this meaning that the differentials of the test functions can be replaced by some weak differentials of the viscosity solution. See [START_REF] Barles | Solutions de viscosité des équations de hamilton-jacobi[END_REF][START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF] for the definition using the Dini-differentials and [START_REF] Clarke | Nonsmooth analysis and control theory[END_REF] for the definition using the proximal differentials. The first statement of this theorem is that V is a viscosity solution of [START_REF] Bonnard | Abnormal geodesics in 2d-zermelo navigation problems in the case of revolution and the fan shape of the small time balls[END_REF]. The proof of this claim relies on the regularity of the value function (continuity) and on the dynamic programming principle. Actually, one can show that any usc function v : [0, T ] × R d which verifies the suboptimality principle, v is a viscosity sub-solution of (28a). Indeed, let φ be in C 1 ([0, T ] × R d ) and (t 0 , x 0 ) ∈ [0, T [×R d local maximum of v -φ. Let u ∈ U , and consider the trajectory x associated with the constant control u = u on [0, T ], and starting from x 0 at time t 0 (i.e., x(t 0 ) = x 0 ). For τ > 0 small enough, we have v(t 0 + h, x(h)) -φ(t 0 + h, x(h)) ≤ v(t 0 , x 0 ) -φ(t 0 , x 0 ), for every 0 ≤ h ≤ τ. On another hand, the sub-optimality principle gives

v(t 0 , x 0 ) ≤ v(t 0 + h, x(h)) + t 0 +h t 0 (s, x(s), u) ds.
We conclude that the test function φ satisfies also the sub-optimality along the trajectory x(•) on [t 0 , t 0 + τ ], that is

φ(t 0 , x 0 ) ≤ φ(t 0 + h, x(h)) + t 0 +h t 0 (s, x(s), u) ds for 0 ≤ h ≤ τ.
Since φ is C 1 , the classical differential calculus leads to

-∂ t φ(t 0 , x 0 ) -f (t 0 , x 0 , u) • D x φ(t 0 , x 0 ) -(t 0 , x 0 , u) ≤ 0.
The above inequality being valid for any u ∈ U , by taking the supremum over u ∈ U , we get that φ satisfies -∂ t φ(t 0 , x 0 ) + H(t 0 , x 0 , D x φ(t 0 , x 0 )) ≤ 0, which concludes that v is a viscosity sub-solution of the HJB equation (28a). Similarly, we can prove that any lsc function that verifies the super-optimality principle is a viscosity super-solution of (28a).

The theorem claims also that the value function is the unique solution of ( 28). This is the consequence of the following equivalences that can be established by non-smooth analysis (see [START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF][START_REF] Vinter | Optimal control[END_REF] for the proof):

v : [0, T ] × R d is usc, v satisfies the sub-optimality principle ⇐⇒ v is a viscosity sub-solution of (28a),

v : [0, T ] × R d is lsc, v satisfies the super-optimality principle ⇐⇒ v is a viscosity super-solution of (28a).

From the point of view of PDEs and viscosity theory, one can also obtain uniqueness by using a general theorem on the comparison principle which can be stated as follows.

Theorem 4.2. Let v 1 be a subsolution of (28) and v 2 be a supersolution of (28

) with v 1 (T, x) ≤ ϕ(x) ≤ v 2 (T, x) for x ∈ R d . Then for any t ∈ [0, T ], x ∈ R d , v 1 (t, x) ≤ v 2 (t, x).
For the proof, we refer to [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]Theorem III.3.7]). The classical method is based on the doubling variable technique (see [START_REF] Crandall | Condition d'unicité pour les solutions generalisées des équations de Hamilton-Jacobi du premier ordre[END_REF][START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF][START_REF] Barles | Solutions de viscosité des équations de hamilton-jacobi[END_REF][START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]). It consists of considering

sup t,s,x,y v 1 (t, x) -v 2 (s, y) - |t -s| 2 + |x -y| 2 ε 2 ,
where ε > 0. It is a regularization technique, called sup/inf-convolution for sub/super-solutions.

Then the regularization of v 1 and u 2 can be considered as the viscosity test functions for v 2 and v 1 respectively, and the comparison result is deduced by the information obtained through the viscosity tests.

Other unconstrained control problems and their HJB formulation

In addition to the Bolza problem, which has been taken here as a model, various other formulations have been considered for optimal control problems, in particular without a final time, or in which the final time is itself a parameter to be optimally chosen. We briefly review some of these formulations, while a more extensive discussion can be found, for example, in [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF].

Minmax control problems. In this class of control problems, the cost function is defined as a maximum cost over the time interval (and not in integral form). More precisely, the control problem reads as

V # (t, x) := min ϕ(x(T )) max s∈[t,T ] Ψ(s, x(s)) (x, u) ∈ X [t,T ] (x) .
Here the cost function is the maximum between the final cost and a maximum running cost along the trajectory. Relation between minmax control problems and state constrained control problems have been noticed and exploited to derive necessary conditions of optimal trajectories, see [START_REF] Vinter | Optimal control[END_REF]Chapter 9]. In Hamilton-Jacobi approach, minmax control problems have been also analyzed in [START_REF] Barron | Viscosity solutions and analysis in L ∞[END_REF]. It has been shown that the value function is Lipschitz continuous and can be characterized by a HJB variational inequality

min ∂ t V # (t, x) + H # (t, x, DV # (t, x)), V # (t, x) -Ψ(t, x) = 0 on [0, T ) × R d , (30a) V # (T, x) = ϕ(x) ∨ Ψ(T, x) for x ∈ R d . (30b)
In this HJB inequality, the running cost function plays the role of an "obstacle". The value function satisfies V # (t, x) ≥ Ψ(t, x). It is also a super-solution of the equation:

∂ t v(t, x) + H # (t, x, D x v(t, x)) = 0 on [0, T ) × R d . (31)
However, V # is a sub-solution of [START_REF] Buckheit | Wavelab and reproducible research[END_REF] only in open sets where V # < Ψ.

Infinite horizon. Infinite horizon problems are intended to model optimal control strategies in the long-time behaviour. Assume that both the dynamics f and the running cost do not depend on time, so that ẋ(s) = f (x(s), u(s)),

with the initial condition x(0) = x. The discounted infinite horizon cost functional is defined as

J(x, u) = ∞ 0 (x(s), u(s))e -λs ds,
where, in addition to the basic assumptions, we require that λ > 0, and, for simplicity, that is uniformly bounded. According to the definition of the cost functional, it is possible to define a value function, which will depend in this case on x alone:

V ∞ (x) := inf J(x, u) (x, u) ∈ X [0,∞) (x) .
The value function can still be characterized as the viscosity solution of a (stationary) Hamilton-Jacobi-Bellman equation, which takes the form

λV ∞ (x) + H(x, DV ∞ (x)) = 0, (33)
in which x ∈ R d , and the Hamiltonian function H is defined by [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF]. As for the regularity of the value function, the basic assumptions imply uniform continuity of the value function. Hölder regularity of V ∞ holds under the assumptions of boundedness for , and global Lipschitz continuity for both and f ; in addition, V ∞ is itself globally Lipschitz continuous if λ is larger than the Lipschitz constant of f (see [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]).

Free final time. In free final time problems, also termed as optimal stopping time problems, the endtime θ of the control interval is itself a free parameter to be chosen in an optimal way. In the simplest case, the dynamics is set in the form [START_REF] Buşoniu | Discounted near-optimal control of general continuous-action nonlinear systems using optimistic planning[END_REF] and the cost functional in the form J(x, (u, θ)) = θ 0 (x(s), u(s))e -λs ds + e -λθ ϕ(x(θ)).

Accordingly, the value function is defined as

V f (x) := inf J(x, (u, θ)) (x, u) ∈ X [0,∞) (x), θ ≥ 0 .
In this case, the HJB equation is stationary, but comes in the form of an obstacle problem:

max λV f (x) + H(x, DV f (x)), V f (x) -ϕ(x) = 0.
The state space is then split in two (possibly overlapping) sets: in the first one the first argument of the max vanishes, so that

λV f (x) + H(x, DV f (x)) = 0,
and the optimal control requires to keep the system evolving; in the second set the second argument of the max vanishes, and therefore

V f (x) = ϕ(x).
As soon as the state of the system enters this set, the optimal strategy requires to stop the system, paying the stopping cost ϕ(x(θ)). Under the basic assumptions, the value function is uniformly continuous (see [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]).

Minimum time. For simplicity, we assume in this paragraph that the dynamics does not depend explicitly on time, that is f (s, x, u) = f (x, u). For minimum time function with timedependent dynamics, we refer to [START_REF] Bokanowski | Minimum time control problems for non autonomous differential equations[END_REF] and the references therein.

In minimum time problems, the goal is to drive the state of the system, in the shortest time, to the final closed set C := {x ∈ R d , g f (x) ≤ 0} (called the target). In general, the possibility of driving the state to the target in finite time may not be ensured for each initial state; this leads to define the so-called backward reachable set R as the set of initial states x which can be driven to the target in finite time. The minimum time control problem is formulated as

T (x) := inf t g f (x(t)) (x, u) ∈ X [0,∞) (x) , (34)
with the convention that T (x) = +∞ when there is no trajectory that starts from x and reaches the target at a finite time. It is not difficult to show that the value function satisfies a dynamic programming principle (DPP) in the form The drawback of this formulation is that the reachable set should be known in advance. An alternative formulation makes use of the so-called Kružkov transformation

T (x) = inf T (x(h)) + h, (x, u) ∈ X [0,∞) (x) .
V(x) = 1 if T (x) = +∞
1 -e -T (x) else, which solves (yet when the time function T is continuous) an auxiliary infinite horizon problem of the form [START_REF] Buşoniu | Continuous-action planning for discounted infinite-horizon nonlinear optimal control with lipschitz values[END_REF], with λ = 1, ≡ 1, and the boundary condition

V(x) = 0 for x ∈ C.
In this case, the reachable set is obtained as a byproduct of the computation of the value function

V, as R = x ∈ R d | V(x) < 1 .
Continuity of the function T is closely related to controllability properties satisfied by the system in a neighbourhood of C and more precisely, to the so-called Small-Time Local Controllability, see [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]Chapter IV]. When the target is smooth enough with a compact boundary, a necessary condition for the continuity of the minimum time function is given by the condition

min u∈U f (x, u) • ν(x) ≤ 0. ( 35
)
This condition is restrictive and excludes a large class of systems with drift. An equivalent condition to the Lipschitz continuity of the minimal time function is the Petrov condition

min u∈U f (x, u) • ν(x) < 0, (36)
which is even more restrictive than [START_REF] Carlini | A weighted essentially nonoscillatory, large time-step scheme for hamilton-jacobi equations[END_REF]. From a geometrical point of view, Petrov condition states that at every point of a neighborhood of the target there exists an admissible control such that the corresponding trajectory points towards the target.

A practical approach to compute the minimum time function and the corresponding optimal trajectories is based on the level set method, introduced by Osher and Sethian. Consider the final cost function as

Φ(x) := max g 1,f (x), • • • , g r,f (x) for x ∈ R d , where (g 1,f , • • • , g r,f
) are the components of the function g f . The level set approach consists of considering the value function V associated to the Mayer problem with final cost Φ, defined by

V (t, x) := min{Φ(x(t)), (x, u) ∈ X [0,t] (x)}. ( 37
)
The value function V can be characterized by an HJB equation as in section 5.1. This function is Lipschitz continuous, while the minimum time function may be discontinuous. Besides, one can notice immediately that when V (t, x) ≤ 0, then there exists an admissible pair (x, u) ∈ X [0,t] (x) such that Φ(x u x (t)) which means that (x, u) satisfies the final constraint. More precisely, the set of positions from where it is possible to reach the target at time t is given by

R(t) := {x ∈ R d | ∃(x, u) ∈ X [0,t] (x), g f (x(t)) ≤ 0} = {x ∈ R d | V (t, x) ≤ 0}.
Therefore, the value function V gives valuable information on the minimum time function and the corresponding trajectories, without any controllability assumption.

Theorem 4.3. Let V be the value function associated to problem [START_REF] Chitour | Singular trajectories of control-affine systems[END_REF]. Then, for every x ∈ R d ,

T (x) := min{t | V (t, x) ≤ 0}. ( 38
)
Furthermore, for x ∈ R d , any optimal trajectory for the minimal time problem [START_REF] Cannarsa | Semiconcave functions, Hamilton-Jacobi equations, and optimal control[END_REF] is also an optimal trajectory of the control problem [START_REF] Chitour | Singular trajectories of control-affine systems[END_REF] where t is fixed to the minimum time given by [START_REF] Clarke | Nonsmooth analysis and control theory[END_REF].

The level set approach provides an effective way to compute the minimum time to reach a target without assuming any specific regularity. Moreover, it has been shown that the level-set method can be generalized to minimum time problems with state constraints (for instance, the case when the trajectory should avoid some obstacles, see Example 1). In this case, the minimum time function is defined as

T (x) := inf t ∃(x, u) ∈ X [0,∞) (x) with g f (x(t)) ≤ 0, and g(x(s)) ≤ 0 on [0, t] .
To use the level set approach in this context, the definition of the value function should be adapted and defined as

V (t, x) := min Φ(x(t)) max s∈[0,t] Ψ(x(s)) | (x, u) ∈ X [0,t] (x) ,
with Ψ(x) := max(g 1 (x), . . . , g m (x)) for every x ∈ R d . With this new definition of the value function V , the statement of Theorem 4.3 remains valid. The new value function corresponds to a minmax problem, so its characterization is given by the HJB inequality (30).

Constrained Bolza problems

In this section, we consider a control problem with state constraints. We denote by K the set of constraints K := {x ∈ R d , g(x) ≤ 0}. We define the set of admissible trajectories by

X g [t,T ] (x) := (x, u) ∈ X [t,T ] (x) | g(x(s)) ≤ 0 for s ∈ [t, T ] .
We adopt the convention V(t, x) = +∞, when the set of admissible trajectories is empty, i.e.,

X g

[t,T ] (x) = ∅. Similarly to the unconstrained case, the value function V satisfies a dynamic programming principle that can be stated as follows.

i) For all x ∈ K, V(T, x) = ϕ(x).

ii) Dynamic programming principle: for all x ∈ K, τ ∈ [0, T ] and h ∈ [0, T -τ ], we have:

V(t, x) = inf (x,u)∈X g [t,T ] (x) V(t + h, x(t + h)) + t+h t (s, x(s), u(s)) ds. ( 39
)

Inward pointing condition

To analyze the properties of the value function V, it is important first to understand the structure of the set the admissible trajectories. This structure depends on an interplay between the dynamics of the state equation and the set of constraints K. Assume in this section that g : R d → R is C 1,1 function, and its zero-level set is very smooth. Consider the following controllability assumption: (HK1) Inward pointing qualification (IPQ) condition: For every R > 0, there exists β > 0 and ρ > 0 such that for every t ∈ [0, T ],

min u∈U f (t, y, u) • ∇g(y) < -β, ∀y ∈ ∂K ∩ B(0, R). ( 40
)
The IPQ condition states that the set of constraints K has a smooth structure and that on each point of the boundary K it is possible to find an admissible control that allows the trajectory to stay in the set K. So, the IPC implies that the set K is weakly invariant. Moreover, the IPC guarantees even a nicer property, called (NFB) Neighbouring feasible trajectories principle, whose proof can be found in [START_REF] Bettiol | On trajectories satisfying a state constraint: w 1,1 estimates and counterexamples[END_REF]Theorem 2.1].

Lemma 4.1. Assume (HK1). Let (t 0 , x 0 ) ∈ [0, T ] × K and let (y, v) ∈ X [t 0 ,T (x 0 ). There exists a constant C and a feasible pair

(x, u) ∈ X [t 0 ,T ] (x 0 ) such that x(t) ∈ K, ∀ t ∈ [t 0 , T ], y -x W 1,1 ([t 0 ,T ];R d ) ≤ Cg + (y(•)), where g + (y(•)) = max t∈[t 0 ,T ] {max(g(y(t)), 0)}.
A consequence of this lemma it that for any x 0 ∈ • K, there exists an admissible trajectory x on [t 0 , T ] such that

x(t 0 ) = x 0 , x(t) ∈ • K, ∀ t ∈ [t 0 , T ].
In fact, for R > 0, there exists a small enough ε > 0 such that

x 0 ∈ K ε := {x | g(x) + ε ≤ 0}, and min ν∈f (t,x,U ) ∇g(x) • ν < - β 2 ∀x ∈ ∂K ε ∩ ∩B(0, R), t ∈ [0, T ],
by the continuity of ∇g and x → f (t, x, U ). Therefore, there exists an admissible trajectory that remains in K ε ⊂ • K. By using similar arguments, we can also prove the following stability result.

Proposition 4.3. Assume (HK1) is satisfied. Consider (K ε ) ε>0 a sequence of subsets of R d such that K ε = {x ∈ R d | g(x) -ε ≤ 0},
and we denote by V ε the value function associated to the control problem with state constraints in K ε (instead of K). Then

lim ε-→0 V ε (t, x) = V(t, x) uniformly on (0, T ) × K.
These stability results lead to the continuity of the value function and to its characterization by a HJB equation on K. Theorem 4.4. Assume (HK1), the value function V(•, •) is uniformly continuous and bounded on [0, T ] × K. Moreover, it is the unique constrained viscosity solution of the following HJB equation:

-∂ t V(t, x) + H(t, x, D x V(t, x)) ≥ 0 on [0, T [×K, (41a)
-∂ t V(t, x) + H(t, x, D x V(t, x)) ≤ 0 on [0, T [× • K, (41b)
V(T, x) = ϕ(x), (41c)
with V(t, x) = +∞ for every x ∈ R d \ K.

It should be noticed that the HJB equation in the above theorem provides only partial information on the boundary of K. Moreover, the function V takes infinite values outside K (i.e., V (t, x) = +∞ for every x ∈ K). These two facts make the approximation of V on K very challenging and require some penalization techniques.

Finally, let us mention that the IPQ condition (HK1) can be weakened a bit and generalized to the case with several constraints, see [START_REF] Frankowska | On relations of the adjoint state to the value function for optimal control problems with state constraints[END_REF]. Nevertheless, this condition remains restrictive for many control problems. Besides, the stability results guaranteed by the IPQcondition ensure the qualification of any optimal trajectory. Therefore, when this assumption is satisfied, any optimal trajectory corresponds to a normal extremal (see Remark 3.2).

Case of state constraints without controllability assumptions

As we mentioned in the previous section, the Lipschitz regularity of the value function requires an interplay between the dynamics f and the set of constraints K. When the controllabity condition (HK1) is not satisfied, the value function may be discontinuous and its characterization by a HJB equation becomes very delicate. In this section, we shall introduce another alternative to deal with state-constrained control problems when the controllability assumption is not satisfied.

Auxiliary Control problem

We set G(y) := max(g 1 (y), • • • , g m (y)) and G f (y) := max(g f,1 (y), • • • , g f,r (y)) for every y ∈ R d . We introduce an auxiliary control problem and its associated value function W defined by

W(t, x, z) := inf y=(x,ζ)∈ Ŝ[t,T ] (x,z) ϕ(x(T )) -ζ(T ) max θ∈(t,T ) G(x(θ)) G f (x(T ))(42) for x ∈ R d , z ∈ R, t ∈ [0, T ],
and where a ∨ b := max(a, b).

In this auxiliary control problem, the term max θ∈[t,T ] g x(θ) is an exact penalization of the state constraints. Here, we shall use the problem [START_REF] Crandall | Condition d'unicité pour les solutions generalisées des équations de Hamilton-Jacobi du premier ordre[END_REF] to characterize the epigraph of the value function ϑ without requiring any additional controllability assumption.

(i) V(t, x) -z ≤ 0 ⇐⇒ W(t, x, z) ≤ 0.
(ii) In addition, the function V is characterized by W through the relation

V(t, x) = min z ∈ R, W(t, x, z) ≤ 0 . (43)
Proof. (i) Let us assume that V(t, x) ≤ z. So there exists a sequence (x n , u n ) n∈N of admissible pairs in X [t,T ] (x), such that

lim n→+∞ T t (s, x n (s), u n (s)) ds + ϕ(x n (T )) -z = V(t, x) -z ≤ 0.
By admissibility, we have for each n ≥ 0, max θ∈[t,T ] G(x n (θ)) ≤ 0 and G f (x n (T )) ≤ 0. Hence

W(t, x, z)≤ lim inf n→∞ T t (s, x n (s), u n (s)) ds + ϕ(y n (T )) -z max θ∈[t,T ] G(x n (θ)) G f (x n (T )) ≤ 0.
Conversely, let us assume that W(t, x, z) ≤ 0. We know that Ŝ[t,T ] (ξ) is a compact set in C 0 ([t, T ]), therefore the infimum in W(t, x, z) is achieved by some trajectory x ∈ S[t, T ]((x, z)) (with an associated control u ∈ U). Moreover,

0 ≥ W(t, x, z) = T t (s, x(s), u(s)) ds + ϕ(x(T )) -z max θ∈[t,T ] G(y(θ)) G f (x(T )).
On the one hand, max θ∈[t,T ] G(x(θ)) G f (x(T )) ≤ 0 and x satisfies the state constraints, and on the other hand,

V(t, x) -z ≤ T t (s, x(s), u(s)
) ds + ϕ(x(T )) -z ≤ 0
which is the desired result. Finally, statement (ii) is an immediate consequence of (i).

Remark 4.2. Should the convexity assumption (H 4 ) not be satisfied, the statements of the above theorem may not hold. Indeed, in general we have

inf{z ∈ R, W(t, x, z) ≤ 0} ≤ V(t, x) ≤ inf{z ∈ R, W(t, x, z) < 0}.
If inf{z ∈ R|W(t, x, z) ≤ 0} = +∞, thereby v(t, x) = +∞ which means that there is no strategies of the first player that guarantee admissibility of trajectories for any action of the second player.

However if inf{z ∈ R|W(t, x, z) ≤ 0} ∈ R, the infimum is reached by some z ∈ R. Furthermore, if inf{z ∈ R|W(t, x, z) < 0} ∈ R, then inf{z ∈ R, W(t, x, z) ≤ 0} = V(t, x) = inf{z ∈ R, W(t, x, z) < 0}.
The auxiliary control problem suggests a reformulation of the state-constrained optimal control problem in an augmented state space. In this new formulation, the constraints are integrated into the functional to be minimized. The value function W is Lipschitz continuous and it can be characterized by an HJB equation without any additional controllability assumption. Theorem 4.6. Assume that (H 0 )-(H 3 ) are satisfied. Then the auxiliary value function is Lipschitz continuous and it is the unique viscosity solution of the following HJB equation

min -∂ t W(t, x, z) + H(t, x, D x W(t, x, z), D z W(t, x, z)), W(t, x, z) -G(x) = 0, on [0, T [×R d × R, W(T, x, z) = ϕ(x) -z G(x) G f (x), on R d × R.
Here the Hamiltonian is defined as

H(t, x, p, q) = max u∈U -f (t, x, u) • p -(t, x, u)q for every (t, x, p, q) ∈ [0, T ] × R d × R d × R.
Let us also point out that any optimal trajectory for the original problem is also a solution of the auxiliary problem when z = z := V (t, x). Conversely, any solution of the auxiliary problem with z = z is an optimal solution of the original state-constrained problem. By consequence, the auxiliary problem provides the value of the original state-constrained control problem and also allow the reconstruction of the optimal trajectories (see [START_REF] Altarovici | A general Hamilton-Jacobi framework for non-linear state-constrained control problems[END_REF][START_REF] Assellaou | Value function and optimal trajectories for a maximum running cost control problem with state constraints. application to an abort landing problem[END_REF]).

Relationship between HJB and PMP

To explain the relationship between HJB and PMP, we first recall a classical result (see for example [START_REF] Fleming | Deterministic and stochastic optimal control[END_REF]), valid under the very restrictive assumption that V ∈ C 2 . Theorem 4.7. We consider the optimal control problem (3) with g = 0 (no state constraint), under the assumptions (H 0 )-(H 1a )-(H 1b ). Assume moreover that f and are continuously differentiable with respect to the space variable, that V ∈ C 2 ([0, T ] × R d ), and that there exists an optimal pair (x * , u * ). Then, the vector p(t) defined by

p(t) := -DV(t, x * (t)) (44)
satisfies the Pontryagin maximum principle.

While this relatively simple result works under unrealistically strong assumptions, more recent theory [START_REF] Vinter | New results on the relationship between dynamic programming and the maximum principle[END_REF] justifies the PMP-HJB relationship (in a suitably weakened form) under general assumptions. In the case of Mayer problems with a locally Lipschitz continuous cost, the sensitivity relations have also been studied in [START_REF] Clarke | The relationship between the maximum principle and dynamic programming[END_REF][START_REF] Vinter | New results on the relationship between dynamic programming and the maximum principle[END_REF]. In these results, the value function is only required to be Lipschitz continuous in a neighborhood of the optimal trajectory. The final cost is not differentiable and therefore, the costate function is not necessarily unique. In this context the sensitivity relations assert that there exists p verifying the PMP, the terminal conditions

-p(0) ∈ ∂ x V(0, x * (0)), -p(T ) ∈ ∂ x V(T, x * (T )), (45a)
and both a partial sensitivity relation

-p(t) ∈ ∂ x V(t, x * (t)), a.e. on (0, T ), (45b)
and a global sensitivity relation

(H(t, x * (t), p(t)), -p(t)) ∈ ∂V(t, x * (t)), for all t ∈ [0, T ]. (45c)
These relations extend (44) by using the generalized gradient of the value function (which is well defined for locally Lipschitz continuous functions). The sensitivity relation (45a) can be simply derived by noting that the optimal solution (x * , u * ) is also solution of the free initial state problem

min ϕ(x(T )) + T 0 (t, x(t), u(t)) dt -V(0, x(0)) ẋ(s) = f (s, x(s), u(s)) on (0, T ) .
Applying necessary optimality conditions to this problem yields a costate arc p 0 . Relation (45a) turns out to be nothing else than the transversality conditions at the endpoints. With the same reasoning, we can notice that (x * , u * ) is also solution of the free initial state problem on [t, T ] for every

0 ≤ t ≤ T min ϕ(x(T )) + T 0 (t, x(t), u(t)) dt -V(t, x(t)) ẋ(s) = f (s, x(s), u(s)) on (0, T ) .
Here again, the optimality condition applied to the free intial point asserts the existence of an adjoint arc p t (which depends on the initial time t). The left-endpoint transversality condition yields the relation p t (t) ∈ ∂V(t, x * (t)). When the final cost function is C 1 -regular, the costate p 0 restricted to [t, T ] is the unique solution to the costate equation on this interval satisfying the right transversality condition. It follows that p 0 (t) coincides with p t (t); the proof of (45b) is then completed. This analysis breaks down if, as is typically the case when either the final cost ϕ is non-smooth or if some state constraints are present. Indeed, in that cases, co-state trajectory may not be unique. The sensitivity relations is in essence a strengthened necessary condition, asserting that it is possible to choose a co-state trajectory to satisfy the sensitivity relation. Some examples are given in [START_REF] Clarke | The relationship between the maximum principle and dynamic programming[END_REF][START_REF] Vinter | Optimal control[END_REF] showing that, in some cases, there are a number of possible choices of co-state trajectories associated with the same optimal control problem, but not all of them satisfy the sensitivity relations.

Numerical methods for HJB

In order to present the general theory for the approximation of viscosity solutions of HJB equations, we refer to an abstract forward problem

v t + H(t, x, Dv) = 0 (t, x) ∈ (0, T ] × R d , v(0, x) = v 0 (x) x ∈ R d , (46)
(for some continuous Hamiltonian H) and set ourselves in the usual finite difference scheme framework. Time is discretized with a (fixed) time step ∆t, so that t k = k∆t; space is discretized with a fixed space step ∆x = (∆x 1 , • • • , ∆x d ). A generic node will be denoted by x j = j∆x, for j ∈ Z d . We also define ∆ = (∆x, ∆t). More general options can be considered, in particular variable time steps and unstructured space grids, but we will restrict here to the basic ideas. In the following, we denote by V n i the desidred approximation of v(t n , x i ), and by V n the set of nodal values for the numerical solution v(t n , •) at time t n . A scheme may be written in compact form as

V n+1 = S(∆; V n ), (47)
where S may be defined in terms of its components S j , for j ∈ Z d .

Monotone schemes

The first and basic convergence theory aimed at approximating HJ equations of the form [START_REF] Falcone | Discrete time high-order schemes for viscosity solutions of hamiltonjacobi-bellman equations[END_REF] uses the concept of monotone scheme. Among the various results, we quote here the Barles-Souganidis theory [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF], which applies to the widest class of schemes and models, including the possibility of treating second-order, degenerate and singular equations. Roughly speaking, this theory states that any monotone, stable and consistent scheme converges to the exact viscosity solution, provided there exists a comparison principle for the limiting equation. Consider a scheme in the general form [START_REF] Falcone | Semi-lagrangian schemes for hamilton-jacobi equations, discrete representation formulae and godunov methods[END_REF]. We recall the concepts of consistency, monotonicity and stability.

Consistency. Let ∆ m = (∆x m , ∆t m ) be a generic sequence of discretization parameters, (t jm , x jm ) be a generic sequence of nodes in the space-time grid such that, for m → ∞,

(∆x m , ∆t m ) → 0 and (t nm , x jm ) → (t, x). (48)
The scheme S is said to be consistent if for any Let φ ∈ C ∞ ((0, T ] × R d ), we have

lim inf m→∞ φ(t nm , x jm ) -S jm (∆ m ; φ(t nm-1 , •)) ∆t m ≥ φ t (t, x) + H(t, x, φ(t, x), Dφ(t, x)), (49a) lim sup m→∞ φ(t nm , x jm ) -S jm (∆ m ; φ(t nm-1 , •)) ∆t m ≤ φ t (t, x) + H(t, x, φ(t, x), Dφ(t, x)), (49b)
where H and H denote respectively the lower and upper semicontinuous envelopes of H.

In [START_REF] Faulwasser | Turnpike properties in optimal control: An overview of discrete-time and continuous-time results[END_REF], the index of the sequence is m, while j m and n m denote the corresponding node indices with respect to the mth space-time grid. Note that, When the Hamiltonian is continuous, then conditions (49) are equivalent to

lim m→∞ φ(t nm , x jm ) -S jm (∆ m ; φ(t nm-1 , •)) ∆t m = φ t (t, x) + H(t, x, φ(t, x), Dφ(t, x)),
Monotonicity. The scheme S is said to be monotone if, for any couple of vectors V and W such that

V j ≥ W j : S j (∆; V ) ≥ S j (∆; W ), (50)
for any ∆x, ∆t satisfying suitable compatibility conditions, that are typically in the form of the so-called CFL conditions.

It is also possible to give a generalized form of definition of monotonicity, to treat some case of high-order scheme. We refer the reader to the discussion carried out in [START_REF] Bokanowski | Value iteration convergence of"-monotone schemes for stationary hamilton-jacobi equations[END_REF].

Given a numerical solution V n , we define its piecewise constant (in time) interpolation V ∆ as

V ∆ (t, x) = I[V n ](x) if t ∈ [t n , t n+1 ) , v 0 (x) if t ∈ [0, ∆t), (51)
where I[V n ](x) denotes an interpolation of the node values in V n , computed at x. We remark that the interpolation operator has to satisfy itself a monotonicity property to obtain a monotone scheme (this holds, for example, for a piecewise linear reconstruction).

We can now state (in a slightly rephrased form) the convergence result given in [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF]:

Theorem 4.8. Assume that (46) satisfies a comparison principle, and let v(t, x) be the unique viscosity solution of [START_REF] Falcone | Discrete time high-order schemes for viscosity solutions of hamiltonjacobi-bellman equations[END_REF]. Assume that (49a)-( 49b) and (50) hold. Assume in addition that the family v ∆t is uniformly bounded in L ∞ . Then, V ∆ (t, x) → V(t, x) locally uniformly on R d ×[0, T ] as ∆ → 0.

This result directly applies to the most classical cases of monotone schemes, as in the examples below.

Finite difference schemes. Given a numerical Hamiltonian

H : [0, T ] × R d × R d × R d → R,
we define an explicit scheme (see [START_REF] Crandall | Two approximations of solutions of Hamilton-Jacobi equations[END_REF]) as follows:

V n+1 i = V n i -∆t H(t n , x i , D -V n i , D + V n i ), (52a)
V 0 i = v 0 (x i ). (52b)
Here, i ∈ Z d , and the terms D -V n i and D + V n i represent finite difference approximations of the gradient at x i . Some examples of such approximations are given here below.

-A monotone (first order) finite difference approximation is defined as

D ± V n i = (D ± k V n i ) 1≤k≤d with D ± k V n i := ± V n i±e k -V n i ∆x k ,
and where {e k } k=1,...,d is the canonical basis of R d ((e k ) k = 1 and (e k ) j = 0 if j = k).

-A second order ENO approximation of D ± k V n i k (see [START_REF] Osher | High essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF]) is given by

D ± k V n i := ± V n i±e k -V n i ∆x k ∓ 1 2 ∆x k m D 2 k,0 V n i D 2 k,±1 V n i (53) with D 2 k,ε V n i := (-V n i+(-1+ε)e k + 2V n i+εe k -V n i+(1+ε)e k )/(∆x k ) 2
then the numerical scheme (52) satisfies all the requirements of Theorem 4.8, and it provides a sequence of solutions V ∆ that converges to the desired solution when ∆ goes to 0.

Some classical choices for the numerical Hamiltonian are listed below.

-Example 1. If the Hamiltonian H is defined by, then a classical choice of numerical Hamiltonian uses upwind approximations as follows

H(t, x, p -, p + ) = max u∈U d i=1 max(-f i (t, x, u), 0) p - i + min(-f i (t, x, u), 0) p + i -(t, x, u) . ( 55
)
This numerical Hamiltonian satisfies all the required conditions of Remark 4.3 for a sufficiently small time step ∆t satisfying the CFL condition, which can be written here as

∆t 1≤i≤d max t,x,u |f i (t, x, u)| ∆x i ≤ 1. (56)
-Example 2. The Lax-Friedrichs scheme can be defined for a very general Hamiltonian.

H LF (t, x, p -, p + ) := H(t, x, p -+ p + 2 ) - d i=1 C i p + i -p - i 2 , (57)
The numerical Hamiltonian H LF satisfies the requirements of Remark 4.3 provided the constants C i are chosen such that C i ≥ max u∈U |f i (t, x, u)| and (∆t, ∆x) satisfies the CFL condition [START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF].

Semi-Lagrangian scheme. The Semi-Lagrangian (SL) scheme is written here directly in the form suitable for the backward dynamic programming equation ( 28)- [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF]. In fact, the SL scheme can be derived by discretizing the Dynamic Programming Principle on a single time step:

   v n-1 j = min α∈U {∆t (t, x, α) + I 1 [V n ](x j + ∆tf (t, x, α)} v N j = ϕ(x j ), (58)
in which I 1 [V ](x) denotes the P 1 (piecewise linear) interpolate of the vector V of node values, computed at the point x.

High-order schemes

In this section, we sketch some basic ideas about high-order approximations for HJ equations, using the simpler model problem

v t + H(Dv) = 0. (59)
Theoretical tools. While the framework of monotone schemes remains the most classical, in the last decades other techniques have been proved to be suitable for the convergence analysis of high-order numerical schemes for Hamilton-Jacobi-Bellman equations, and in particular:

• ε-monotonicity

As remarked in the original paper [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF], the Barles-Souganidis theorem allows for a small (more precisely, an o(∆t)) monotonicity defect. This makes it possible to prove convergence for quasi-monotone schemes, although the property of quasi-monotonicity is more subtle, and more difficult to prove. This theory has been applied to high-order SL schemes and to filtered schemes (for which the monotonicity defect can be set a priori), as in [START_REF] Augoula | High order numerical discretization for hamilton-jacobi equations on triangular meshes[END_REF][START_REF] Falcone | Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations[END_REF][START_REF] Bokanowski | Value iteration convergence of"-monotone schemes for stationary hamilton-jacobi equations[END_REF].

• Lin-Tadmor theory Lin-Tadmor convergence theory is inspired by the Lip'-stability theory for conservation laws (see [START_REF] Lin | l 1 -stability and error estimates for approximate Hamilton-Jacobi solutions[END_REF]). Here, a different concept of stability is singled out, i.e., the concept of semi-concave stability: Definition 4.2. A family of approximate solutions v of (59) is said to be semi-concave stable if there exists a function k(t) ∈ L 1 ([0, T ]) such that

D 2 u (t, x) ≤ k(t)I ( 60
)
for t ∈ [0, T ].
is the conservation law (or system of conservation laws, see [START_REF] Jin | Numerical passage from systems of conservation laws to hamilton-jacobi equations, and relaxation schemes[END_REF]) associated to the HJ equation. Following [START_REF] Hu | A discontinuous galerkin finite element method for hamilton-jacobi equations[END_REF], the approximate solution is constructed in the space

V k ∆x = w : w| I j ∈ P k (I j ) ,
in which I j denotes the j-th element of the computational domain, for example the interval

[x j-1/2 , x j+1/2 ] = [x j -∆x/2, x j + ∆x/2]
in a one-dimensional uniform grid. A DG scheme of order k for the equation ( 59) is defined by looking for the function w ∈ V k ∆x such that

I j w xt φdx - I j H(w x )φ x dx + H j+1/2 φ - j+1/2 -H j-1/2 φ + j-1/2 = 0 (64)
for all j and all φ ∈ V k-1 ∆x . In [START_REF] Kurganov | Adaptive central-upwind schemes for hamilton-jacobi equations with nonconvex hamiltonians[END_REF], the values H j-1/2 are defined by

H j±1/2 := H w x (x - j±1/2 ), w x (x + j±1/2 ) ,
and H is as usual a monotone numerical Hamiltonian. The outcome of the scheme is the derivative w x (or, in the multi-dimensional case, all the partial derivatives). The approximate solution w is recovered via suitable post-processing techniques.

Filtered schemes. The general idea of filtered schemes (which had previously appeared in the context of conservation laws as flux-limiter schemes) is to provide a clever coupling between a monotone and a high-order scheme. The high-order scheme is applied only where the solution is regular enough, this being detected by a suitable "filter" function.

The construction of a filtered scheme needs three ingredients: a monotone scheme S M , a high-order scheme S HO and a bounded filter function F : R → R. The filtered scheme S F is then defined as

v n+1 j = S F j (V n ) := S M j (V n ) + ∆tF S HO j (V n ) -S M j (V n ) ∆t , (65)
where ε = ε(∆) > 0 is a parameter vanishing for ∆t, ∆x → 0, which controls the monotonicity defect of the filtered scheme (more hints on the choice of ε and of the filter function can be found in [START_REF] Froese | Convergent filtered schemes for the Monge-Ampère partial differential equation[END_REF][START_REF] Bokanowski | An efficient filtered scheme for some first order timedependent hamilton-jacobi equations[END_REF]). In constructing the filter function, the basic idea is that "large" values of the ratio ρ = (S HO -S M )/(∆tε) indicate a singularity (where the scheme needs to be monotone), while "small" values indicate a smooth region (in which the scheme can be high-order). It can be shown that, for a suitable choice of , the filtered scheme converges to the viscosity solution by quasi-monotonicity.

Further comments. Several advances have been made to improve the numerical schemes of approximations of HJB equations, in particular in high dimension. Let us mention the resolution techniques on sparse grids [START_REF] Bokanowski | An adaptive sparse grid semi-lagrangian scheme for first order Hamilton-Jacobi Bellman equations[END_REF][START_REF] Garcke | Suboptimal feedback control of PDEs by solving HJB equations on adaptive sparse grids[END_REF], on tree structures as in [START_REF] Alla | An efficient DP algorithm on a tree-structure for finite horizon optimal control problems[END_REF], or approximation by sophisticated model reduction techniques as in [START_REF] Alla | An adaptive POD approximation method for the control of advectiondiffusion equations[END_REF]. We also mention antoher chapter [START_REF] Ballarin | Space-time pod-galerkin approach for parametric flow control[END_REF] of the Volume 1 of this Handbook, dedicated to model reduction methods. All these methods aim at providing accurate numerical approaches for solving efficiently HJB equations with reasonable numerical efforts (measured by complexity of algorithms, CPU time and accuracy).

Optimal trajectory reconstruction from the value function

From a control viewpoint, the approximation of the value function V has a relatively lesser interest with respect to the construction of the (approximate) optimal control. We propose in this section some algorithms that lead, given an approximations of the value function, to construct a quasioptimal controls in feedback form. The procedure does not depend on the specific scheme used to compute V. For simplicity, we consider the trajectory reconstruction on the time interval [0, T ]. However, all the results remain valid for a reconstruction on any sub-interval [t, T ]. For h > 0, consider an integer n h ∈ N, a partition

s 0 = 0 < s 1 < • • • < s n h = T of [0, T ], denote h k := s k+1 -s k and that max k h k ≤ h.
Consider a numerical approximation f h of the dynamics f such that, for every R > 0, we have:

|f h (t, x, u) -f (t, x, u)| ≤ C R h, ∀t ∈ [0, T ], |x| ≤ R, u ∈ U, (66)
where the constant C R is independent of h ∈ [0, 1]. An approximation scheme for the differential equation ẋ(t) = f (t, x(t), u) (for a constant control u, discrete times s k and time step h k = s k+1 -s k ) can be written as

y k+1 = y k + h k f h k (s k , y k , u), k ≥ 0. ( 67
)
Here and in the sequel we use the notation y k to denote a state at discrete times. The case of the Euler forward scheme corresponds to the choice f h := f. Higher order Runge Kutta schemes can also be written as [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI-relaxations[END_REF] and with a function f h satisfying [START_REF] Lasserre | Moments and convex optimization for analysis and control of nonlinear partial differential equations[END_REF]. For instance, the Heun scheme corresponds to the choice

f h (t, y, u) := 1 2 (f (t, y, u) + f (t + h, y + hf (t, y, u), u)).
Bolza or Mayer problems. Consider first, the case of Bolza or Mayer unconstrained problems. Let V the exact value function defined in [START_REF] Alla | An adaptive POD approximation method for the control of advectiondiffusion equations[END_REF]. Let V h be an approximation of V, and define E h as a uniform bound on the error:

|V h (t, x) -V(t, x)| ≤ E h , ∀t ∈ [0, T ], |x| ≤ R,
with R > 0 large enough. The approximate feedback is defined on the basis of the approximate value function with a discrete dynamic programming procedure. Following some arguments introduced in [START_REF] Rowland | Construction of optimal feedback controls[END_REF], it can be shown that any cluster point of (y h ) h>0 is an optimal trajectory that realizes a minimum in the definition of the original control problem.

Theorem 4.10. Assume (H 0 ), (H 1 ), (H 2 ) and (H 4 ). Assume also that the approximation (66) is valid and the error estimate E h = o(h). Let (y, z) be in R d × R and let (y h k ) be the sequence generated by Algorithm (TR).

(i) The approximate trajectories (y h k ) k=0,...,n h constitute a minimizing sequence in the following sense:

V(0, x) = lim h→0 ϕ(y h n h ) + h n h k=0 (s k , y h k , u h k ) . (70)
(ii) Moreover, the family (y h ) h>0 admits cluster points, for the L ∞ norm, when h → 0. For any such cluster point ȳ, we have ȳ ∈ S [0,T ] (x) and x is an optimal trajectory for V(0, x).

Algorithm (TR) -Trajectory reconstruction algorithm for Bolza problems

Require: First we set y h 0 := x. 1: For k = 0, . . . , n h -1, knowing the state y h k we define (i) an optimal control value u h k ∈ U such that

u h k ∈ argmin u∈U V h s k , y h k + h k f h (s k , y h k , u) + h (s k , y h k , u) (68)
(ii) a new state position y h k+1

y h k+1 := y h k + h k f h (s k , y h k u h k ). (69)
2: return a piecewise constant control u h (s) := u h k on [s k , s k+1 [, and a piecewise linear trajectory y h such that y h (s k ) = y h k .

Let us emphasize that the condition E h = o(h) indicates that the approximation V h should be provided with a given precision. Typically, a numerical scheme would provide an approximation V ∆ for ∆ = (∆t, ∆x). Under a CFL condition, the error estimate is of order 0( √ ∆t. To ensure that V ∆ -V = o(h), it suffices to take √ ∆t = o(h).

Minmax problems. The same idea as in Algorithm (TR) can be adapted for minmax control problems. Let V #,h be an approximation of the value function V # ,. Here again, the function V #,h could be a numerical approximation obtained by solving a discretized form of the HJB equation.

Algorithm (TRM) -Trajectory reconstruction for Minmax problems

Require: Set y h 0 = x. 1: Define the positions (y h k ) k=0,...,n h , and control values (u h k ) k=0,...,n h -1 , by recursion as follows. or k = 0, . . . , n h -1, knowing the state y h k we define

(i) an optimal control value u h k ∈ U such that u h k ∈ argmin u∈U V #,h s k , y h k + h k f h (s k , y h k , u), z Ψ(s k , y h k ) (71)
(ii) a new state position y h k+1

y h k+1 := y h k + h k f h (s k , y h k , u h k ).
Note that in (71) the value of u h k can also be defined as a minimizer of u → V #,h s k , y h k + h k f h (s k , y h k , u), z , since this will imply in turn to be a minimizer of (71) Theorem 4.11. Assume (H 0 ), (H 1 ), (H 2 ) and (H 4 ) hold true. Assume also that the approximation (66) is valid and that V #,h is an approximation of V # with error estimate V # -V #,h = o(h). let (y h k ) be the sequence generated by Algorithm 2. (i) The approximate trajectories (y h k ) k=0,...,n h constitute a minimizing sequence in the following sense:

V # (0, y) = lim h→0 max 0≤k≤n h Ψ(s k , y h k ) ϕ(y h n h ).
(ii) Moreover, the family (y h ) h>0 admits cluster points, for the L ∞ norm, when h → 0. For any such cluster point ȳ, we have ȳ ∈ S [0,T ] (y) and ȳ is an optimal trajectory for V # (0, x).

Numerical examples

Test 1 -Unconstrained Zermelo problem In this example, we consider the same setting as in Section 3.4. To compute the minimal time function, we use the level-set approach described in Section 4.2. We consider a domain of computation large enough to contain the initial position (0, 0) and the final target y f = (20, 1). More precisely, the computation will be performed on D := [-1, 21] × [-0.5, 1.5]. We use the finite difference scheme (52) combined with ENO approximation ( 53) and the Lax-Friedrichs numerical Hamiltonian (57). Figure 7 shows the optimal solution computed on a uniform grid with 500 × 100 nodes. The optimal time to steer the system from the initial point (0, 0) to the final state y f = (20, 1) is 4.94916. While the approximation of the optimal trajectory on Figure 7(left) seems quite accurate, the optimal control law in Figure 7(right) presents some oscillations.

In this example, the minimum time function T is continuous on its domain, as it can be seen on Figure 8, where some level-sets of the T are presented. Recall that for this version of Zermelo problem, the shooting method provides a global optimal trajectory or a local optimal trajectory depending on the initialization of the adjoint state (see Section 4). By using the HJB approach, the reconstructed trajectory is an approximation of the global trajectory even when performing the calculations on very coarse grids. Note also that an approximation of the derivative of the minimal time function at the initial position (0,0) can be An initialization of the shooting method with this vector allows the method to converge towards the global solution.

Test 2 -Unconstrained Zermelo problem -Case with a strong current. Here, we consider a variant of the Zermelo problem where the dynamics is given by

ẏ1 (s) = v(s) cos(u(s)) + 2 - 1 2 y 2 2 (s), (72a) ẏ2
(s) = v(s) sin(u(s)), (72b)
where the control functions are the speed u 1 (s) and the angle of orientation u 2 (s) of the boat (such that u(s)(u 1 (s), u 2 (s)) ∈ U := [0, 1] × [0, 2π], for a.e. s ∈ [0, T ]). In this example, the drift is strong in the middle of the channel R × [-2, 2] and is zero along the channel banks. The target is a ball centred at the origin and with radius r = 0.2. The domain of computation is

D := [-5, 2] × [-2, 2]. -Test 2-1:
The numerical simulations are performed on a uniform grid with N x = N y = 100 nodes on each axis.

-Test 2-2 : Here, the numerical simulations are performed on a uniform grid with N x = N y = 500 nodes on each axis.

Test 3 -Constrained Zermelo problem. Consider again the same dynamics for Zermelo problem as in the previous paragraph. Now, the state is required to avoid an obstacle that is a rectangle [-2.5, -1.5] × [-0.5, 0.5]. Figure 11 shows the level sets of the minimum time function and some sample optimal trajectories. Figure 11(left) displays the results obtained with 250×250 nodes on the domain of computation, while figure 11(right) corresponds to computation on grid of 500×500 nodes. Note that optimal trajectories tend to avoid the central section of the channel, in which the current is stronger; moreover, gradient singularities downstream of the obstacles indicate points at which the optimal trajectory is not unique (e.g., it can go either left or right of the obstacle).

Test 4 -Goddard Problem. Now, consider the Goddard problem as described in Sections 2.4 & 3.5. The problem is with free final time and state constraints. Since the inward pointing condition is not satisfied, we reformulate the problem as described in Section 4.4.1. We compute the auxiliary value function on D = [1.0, 1.2] × [0, 0.12] × [0, 1]. The optimal trajectory obtained from this approximation indicates that the constraint on the variable v is saturated during a time interval. In Figures 12 and13, we present the obtained approximate solution on a grid with N x 4 node points.

Figure 12 corresponds to a simulation with N x = 20, while Figure 13 corresponds to a simulation with N x = 40. In both figures, we show the three state variables in the top line, and the adjoint states in the middle line. In the third line, the control variable is displayed on the left. In this example again, the approximation of the control law presents high oscillations. We present also, in the third line, an approximation of the derivative of the Hamiltonian with respect to variable u. We notice that this derivative vanishes on come time interval, which indicates that on this interval, the control is singular. Moreover, we notice that the constraint on the velocity is saturated on another time interval.

Optimistic planning algorithms

As mentioned previously, direct and indirect methods are quite simple to implement and they provide local solutions with a great precision. These methods depend on the initialization -especially the shooting method which is particularly sensitive to the initialization and also requires an a priori knowledge of the optimal trajectory structure. On the other hand, the HJB approach always provides a global optimum, but it requires a much greater computational effort because the value function must first be computed in a quite large domain.

The approach that we present here is a global approach based on a discretization in the space of controls, combined with optimistic planning (OP) algorithms [START_REF] Buşoniu | Discounted near-optimal control of general continuous-action nonlinear systems using optimistic planning[END_REF][START_REF] Buşoniu | Continuous-action planning for discounted infinite-horizon nonlinear optimal control with lipschitz values[END_REF] (without requiring any discretization of the state space). This approach is interesting especially for applications where the control dimension r is much lower compared to the state dimension d. On a given discretization of the time interval, our approach will seek to identify the best control strategy to apply on each time sub-interval. The OP methods perform the optimal control search by branch and bound on the control set, always refining the region with the best lower bound of the optimal value (this is what justifies the term "optimistic"). An interesting feature of these algorithms is the close relationship between computational resources and quasi-optimality, which exploits some ideas of reinforcement learning [START_REF] Munos | From bandits to Monte-Carlo tree search: the optimistic principle applied to optimization and planning[END_REF]. Indeed, for given computational resources, the OP approaches provide a sub-optimal strategy whose performance is as close as possible to the optimal value (with the available resources).

First, for N ≥ 2, consider a uniform partition of [0, T ] with N + 1 time steps: t k = k∆t, k = 0, . . . , N , where ∆t = T N is the step size. For a sequence of actions u = (u k ) 0≤k≤N -1 ∈ U N , we consider the discrete-time dynamical system

y 0 = x, y k+1 = F k (y k , u k ) k = 0, ..., N -1, (73)
where F k (x, u) is an approximation of the solution to the system ẋ(s) = f (x(s), u) on (t k , t k+1 ), with the initialcondition x(t k ) = x. More precisely, we assume that

F k (x, u) -x - t k +1 t k f (s, x(s), u) ds = O(∆t) for all x ∈ R d , u ∈ U.
Consider also an instantaneous cost function L k that approximates the integral of over an interval [t k , t k+1 ], for k = 0, ..., N -1:

L k (x, u) - t k +1 t k (s, x(s), u) ds = O(∆t) for all x ∈ R d , u ∈ U.
In this section, we assume that (H 0 )-(H 4 ) are satisfied and that f and are Lipschitz continuous with respect to the control variable. The approximations F k and L k are also assumed to be Lipschitz continuous:

F k (x, u) -F k (x , u ) ≤ L F,x x -x + L F,u u -u , |L k (x, u) -L k (x , u )| ≤ L L,x x -x + L L,u u -u ,
Notice that W (x, z) converges to W(0, x, z) as N → ∞ (i.e., ∆t = T N → 0), where the continuous value function W is defined in [START_REF] Crandall | Condition d'unicité pour les solutions generalisées des équations de Hamilton-Jacobi du premier ordre[END_REF]. Under Assumption (H 4 ), the error estimate of |W (x, z) -W(0, x, z)| is bounded by O( 1N ), see [START_REF] Bokanowski | Optimistic planning algorithms for state-constrained optimal control problems[END_REF]Appendice B]. Furthermore, the sequence of discretetime optimal trajectories (for N ∈ N) provide convergent approximations of optimal trajectories of the continuous problem [START_REF] Rowland | Construction of optimal feedback controls[END_REF], see [START_REF] Assellaou | Value function and optimal trajectories for a maximum running cost control problem with state constraints. application to an abort landing problem[END_REF].

Remark 5.1. The discrete dynamics and the discrete cost can be obtained by a discretization scheme of the time-continuous function f and .

It is worth mentioning that he algorithms that will be presented in this section allow also to consider situations where the dynamics F and L are obtained by some statistical models which can be enriched during the computational process.

With similar arguments as in the proof of (43), we have

V (x) = inf{z | W (x, z) ≤ 0}. (78)
Furthermore, a pair (ȳ, ū) is optimal for problem [START_REF] Rifford | On the stabilization problem for nonholonomic distributions[END_REF] if and only if it is optimal for problem (77) with z = V (x).

For the sake of simplicity and without loss of generality, we suppose that the control is of dimension r = 1 and we denote by D its maximal diameter (∀a, a ∈ A, a -a ≤ D). Let us emphasize that the approach can be generalized to control variables with multiple dimensions.

Planning algorithms are based on the principles of optimistic optimization. In order to minimize the objective function J over the space U N , we refine, in an iterative way the search space into smaller subsets. A search space, called node and denoted by U i with i ∈ N, is a Cartesian product of sub-intervals of U , i.e., U i := U i,0 × U i,1 × • • • × U i,N -1 ⊆ U N , where U i,k represents the control interval at time step k, for k = 0, . . . , N -1. The collection of nodes will be organized into a tree Υ that will be constructed progressively by expanding the tree nodes. Expanding a node U i , with i ∈ N, consists in choosing an interval U i,k , for k = 0, . . . , N -1, and splitting it uniformly to M sub-intervals where M > 1 is a parameter of the algorithm. The order of expanded nodes and the intervals that have to be split will be chosen in such a way to minimize the cost J. For now, we introduce some useful notations related to the tree Υ: • We associate, for any node U i ∈ Υ, a sample sequence of controls u i := (u i,k ) N -1 k=0 ∈ U i such that u i,k corresponds to the midpoint of the interval U i,k for any k = 0, . . . , N -1.

• Denote d i,k , for k = 0, ..., N -1, the diameter of the interval U i,k of some node U i ∈ Υ. In particular, k) , where s i (k) indicates the number of splits needed to obtain the interval U i,k for k = 0, . . . , N -1.

d i,k = D M s i (
• The depth p i of a node U i is the total number of splits done to obtain this node:

p i := N -1 k=0 s i (k). (79)
We denote by Depth(Υ) the maximal depth in the tree Υ. • A node U i is a tree leaf if it has not been expanded. The set of tree leaves is denoted by Λ.

• Finally, we denote by Λ p := U i ∈ Υ s.t. p i = p the set of leaves of Υ of depth p ∈ N.

By selecting controls at the intervals centers and by taking M odd, we guarantee that after expanding a node U i we generate at least one node U j with J(x, z, u j ) ≤ J(x, z, u i ). Indeed, the middle child U j contains the control sequence of U i . Proposition 5.1. By the tree construction, there exists at least a leaf node U i ∈ Λ containing an optimal control sequence and satisfying J(x, z, u i ) -σ i ≤ W (x, z) ≤ J(x, z, u i ), [START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF] where u i is the sample control sequence in U i and where

σ i := N -1 k=0 β k d i,k N -1 k=0 γ k d i,k , (81)
with β k and γ k positive constants only depending on the Lipschitz constants of F, L, Φ and of Ψ.

In the optimistic planning algorithms, at each iteration, one or several optimistic nodes are chosen and split to get from each node M children (M > 1 is a fixed parameter of the algorithm). To expand a node U i , we choose an interval from U i,0 × U i,1 × • • • × U i,N -1 and we partition it uniformly to M sub-intervals. If we choose to split the interval U i,k , for some k = 0, . . . , N -1, then M nodes with will be generated and then the new error term σ + i (k) is

σ + i (k) := N -1 j=0,j =k β j d i,j + β k d i,k M N -1 j=0,j =k γ j d i,j + γ k d i,k M .
Henceforth, in order to minimize the error σ + i (k), the best choice of the interval to split, k * i , is given by:

k * i ∈ argmin 0≤k≤N -1 σ + i (k). (82)
The following result gives an upper bound on the error term σ i , of any node U i ∈ Υ.

Proposition 5.2. Assume that the number of split M > L F,x > 1. Let τ := log M log(L F,x ) . Consider a node U i at some depth p i = p. For p large enough, the error σ i (defined in [START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF]) is bounded as follows:

σ i ≤ δ p := c 1 (N ) ∆t M -p N , (83)
where c 1 (N ) := C 1 1-M 1/τ L F,x M q(N ) , q(N ) := 2 -(N -1) τ -2 2τ (τ -1) , and C ≥ 0 are independent of N and p.

Notice that, by definition, we have τ ≥ 2, and L -1 F,x M 1/τ ≤ 1. Hence ( 83) is meaningful only when the strict inequality L -1 F,x M 1/τ < 1 holds. Moreover, q(N ) ≤ 2, hence c 1 (N ) is bounded independently of N .

For a given depth p ∈ N, we define the set of nodes that will be expanded by optimistic planning algorithms:

Υ * p := {U i ∈ Υ at depth p | J(x, z, u i ) -δ p ≤ W (x, z)},
where δ p is defined as in [START_REF] Vinter | New results on the relationship between dynamic programming and the maximum principle[END_REF]. The set containing all expanded nodes : Υ * := q≥0 Υ * q is in general smaller than the whole tree Υ. The next definition introduces a notion of branching factor that is a sort of complexity of the algorithm. The precise definition of this notion is given hereafter.

Now, we will present the rules for refining the search of an optimal control strategy. In the first algorithm, at each iteration, the node U i * minimizing the lower bound (J(x, z, u i ) -σ i ) will be selected and split to M children. More precisely, we identify an interval U i * ,k * i * whose partition in M sub-intervals will produce the lowest error σ i * (k * i * ).

Algorithm Optimistic Planning (OP)

Require: The number of intervals N , the split factor M , the maximal number of expanded nodes I max 1: Initialize Υ with a root U 0 := U N and n = 0 (n := number of expanded nodes). 2: while n < I max do 3:

Select an optimistic node to expand: U i * ∈ argmin U i ∈Λ (J(x, z, u i ) -σ i ).

4:

Select k * i * , defined in [START_REF] Tsiotras | Goddard problem with constrained time of flight[END_REF], the interval to split for the node U i * .

5:

Update Υ by expanding U i * along k * i * and adding its M children.

6:

Update n = n + 1. 7: end while 8: return Control sequence u i * = (u i * ,k ) k ∈ U N of the node U i * ∈ argmin U i ∈Λ J(x, z, u i ).

Theorem 5.1. Assume that M > L F,x > 1. Let u i * and J(x, z, u i * ) be the output of the OP algorithm, and let n ≥ 1 be the corresponding number of expanded nodes. For n large enough, we have

0 ≤ J(x, z, u i * ) -W (x, z) ≤ E OP (N, n) := c 2 (N ) n - log(M ) N log(m) if m > 1, c 2 (N )(log(n)) -log(M ) N R if m = 1, (85)
where c 2 (N ) > 0 is bounded uniformly with respect to N .

In Algorithm , the number I max represents a maximal available computational resource. The number of expanded nodes corresponds to the number of iterations, since at each iteration only one node is expanded. Other optimistic planning methods can be considered. For instance, the simultaneous optimistic planning (SOP) algorithm or simultaneous optimistic planning with multiple steps (SOPMS) algorithm that expand at each iteration several nodes at every iteration.

Test 5. To show the relevance of this approach, we consider a variant of Zermelo problem where a boat targets the set C := B(0; 0.1), at time T = 1, with minimal fuel consumption. The dynamics is similar to the one considered in [START_REF] Osher | High essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF]. We consider also two rectangular obstacles with horizontal and vertical half lengths (r x , r y ). The first obstacle is centered at (-2.0, 0.5) with (r x , r y ) = (0.4, 0.4), and the second obstacle is centered at (-2.5, -1) with (r x , r y ) = (0.2, 1). To take into account the pointwise and final state constraints, we define the functions g and g f by

Remark 1 . 1 .Remark 1 . 2 .

 1112 Note that if ≡ 0 (Mayer problem), (H 4 ) reduces tof (s, x, U ) convex for all (s, x) ∈ [0, T ] × R d . For 0 ≤ a < b ≤ T and (x, z) ∈ R d × R,consider the set of all trajectories satisfying (4) on the time interval [a, b], for a control input u ∈ U, and starting at a position (x, z) at time a: S [a,b] (x, z) := {(x, z) ∈ W 1,1 ([a, b]) | ∃u ∈ U such that (x, z, u) satisfies (4), with x(a) = x, z(a) = z}.

Figure 2 :

 2 Figure 2: Goddard problem: result of the direct code, control values. A structure with four different subarcs is observed. A simple inspection indicates a concatenation of a bang arc (u 1), possibly of a singular arc interior to the control bounds (u ∈ (0, 1)), of a boundary arc due to the state constraint on the velocity, and of a final bang arc (u 0). This structure information, of combinatorial nature, is obtained here without any a priori knowledge on the solution. It is the key to define an appropriate shooting function and solve very accurately the problem thanks to an indirect method.

Figure 3 :

 3 Figure 3: Zermelo problem with obstacle: result of the shooting method. The state, costate and control are portrayed for the two solutions (penalisation parameter α = 1e -3). The first solution (in blue) is slightly shorter (T 4.98) than the second one (yellow, T 4.99).

#

  Regular maximized Hamiltonians and associated flow H(y, p, u) = -1.0 + Vmax*p'*u + p[1]*h(y[2]) ur(p) = p / sqrt(p[1]^2+p[2]^2) Hr(y, p) = H(y, p, ur(p)) fr = Flow(Hr) # Shooting function function shoot(y1, p1, tf) yy0, p0 = fr(tf/2.0, y1, p1, t0) yyf, pf = fr(tf/2.0, y1, p1, tf) s = zeros(eltype(y1), 5)

Figure 4 :

 4 Figure 4: Zermelo problem, turnpike property: result of the variant of the shooting method. The state and control are portrayed, illustrating that most of the time is spent close to y 2 = 1/2 where the current is maximum. Accordingly, u 2 0 on the corresponding control subarc.

(

  Same notation used for H 001 .) Plugging this dynamic feedback control into the original Hamiltonian (with µ = 0) defines the singular Hamiltonian H s (x, p) := H(x, p, u s (x, p), µ = 0)

  Computation of singular control of order 1 H0(x, p) = p' * F0(x) H1(x, p) = p' * F1(x) H01 = Poisson(H0, H1) H001 = Poisson(H0, H01) H101 = Poisson(H1, H01) us(x, p) = -H001(x, p)/H101(x, p) # Computation of boundary control g(x) = vmax-x[2] # vmax -v ≥ 0 ub(x) = -Lie(F0, g)(x) / Lie(F1, g)(x) µb(x, p) = H01(x, p) / Lie(F1, g)(x) # Hamiltonians (regular, singular, boundary) and associated flows H(x, p, u, µ) = H0(x, p) + u*H1(x, p) + µ*g(x) Hr(x, p) = H(x, p, 1.0, 0.0) Hs(x, p) = H(x, p, us(x, p), 0.0) Hb(x, p) = H(x, p, ub(x), µb(x, p)) Then, to integrate the Hamiltonians to obtain the flows and define the shooting function in terms of the initial adjoint p 0 , the entry point of the singular arc t 1 , the entry point of the boundary arc t 2 , the exit point t 3 , and the free final time T : f0 = Flow(H0) fr = Flow(Hr) fs = Flow(Hs) fb = Flow(Hb) # Shooting function function shoot(p0, t1, t2, t3, tf) x1, p1 = fr(t0, x0, p0, t1) x2, p2 = fs(t1, x1, p1, t2)

Figure 5 :

 5 Figure 5: Goddard problem: result of the shooting method, states and costates. The values are in line with those obtained by the direct method. They are computed by concatenating the integration of the four Hamiltonian flows involved.

5

 5

Figure 6 :

 6 Figure 6: Goddard problem: result of the shooting method, control and state constraint multiplier. The bang-singular-constrained-bang structure of the control is very accurately determined, while the state constraint multiplier is positive only along the boundary arc.

Definition 4 . 1 .

 41 (viscosity solution) Let V : [0, T ] × R d -→ R.

Theorem 4 . 1 .

 41 Suppose that (H 0 )-(H 1a )-(H 1b ) and (H 3 ) hold. Then the value function V is the unique viscosity solution of (28) in the sense of Definition 4.1.

  When the backward reachable set is open and the minimum time function is continuous, then the DPP leads to a characterization of T by the HJB equation sup u∈U -f (x, u) • DT (x) = 1, complemented with the boundary condition T (x) = 0 for x ∈ C and lim x→z T (x) = +∞ for z ∈ ∂R.

Theorem 4 . 5 .

 45 Assume that (H 0 )-(H 4 ) are satisfied. Then for any t ∈ [0, T ] and (x, z) ∈ R d ×R, there holds

Remark 4 . 3 .Friedrich

 43 , and where m(a, b) := a if ab > 0 and |a| ≤ |b|, m(a, b) = b if ab > 0 and |a| > |b|, and m(a, b) = 0 if ab ≤ 0. The numerical Hamiltonian can be defined in different ways to satisfy the requirements of the convergence Theorem 4.8. In particular, if the numerical Hamiltonian H is Lipschitz continuous on all its arguments, consistent with H (i.e. H(t, x, p, p) = H(t, x, p)) and monotone (i.e., ∂H ∂p - k (t, x, p -, p + ) ≥ 0, ∂H ∂p + k (t, x, p -, p + ) ≤ 0) together with the following Courantp -, p + ) + ∂H ∂p + k (y, p -, p + ) ≤ 1,

Figure 7 :

 7 Figure 7: (Test 1) Approximation of the optimal trajectory for Zermelo problem with HJB approach.

Figure 8 :

 8 Figure 8: (Test 1) Approximation of the minimum time function on [-0.5, 20] × [-0, 5, 1.1].

Figure 9 :

 9 Figure 9: (Test 2-1) Unconstrained Zermelo problem: minimal time function (left) and optimal trajectories (right). N x = 100

Figure 10 :

 10 Figure 10: (Test 2-2) Unconstrained Zermelo problem: minimal time function (left) and optimal trajectories (right). N x = 500

Figure 11 :

 11 Figure 11: (Test 3) Constrained Zermelo problem: minimum time function and optimal trajectories. Nx=250 (left) and Nx=500 (right)

Figure 13 :

 13 Figure 13: (Test 4) Goddard problem solved by HJB approach, with Nx=40

  for every x ∈ R d and u ∈ U . Moreover, we assume that the Lipschitz constants of F k and L k are related to the Lipschitz constants of f and by the following relationsL F,x := 1 + ∆tL f,x C L F,u := ∆tL f,u C and (74a) L L,x := ∆t 2 L ,x C, L L,a := ∆t(L ,a + C),(74b)where the constant C > 0 may depend on ∆t and the Lipschitz constants of f and . Now, we define the state-constrained optimal control problemV (x) := inf u=(u k ) k ∈U N N -1 k=0 L k (y x,u k , u k ) + ϕ(y x,u N ) | g(y x,uk ) ≤ 0 ∀k = 0, ..., N, g f (y x,u N ) ≤ 0 .

( 75 ) 1 k=0L

 751 Then, for the discrete auxiliary control problem, we define the cost functional J byJ(x, z, u) := Nk (y x,u k , u k ) + ϕ(y x,u N ) -z max 0≤k≤N g(y x,u k ) g f (y x,u N ) (76) (for (x, z, u = (u k )) ∈ R d × R × U N) and the corresponding auxiliary value is defined, for(x, z) ∈ R d × R, by W (x, z) := inf u=(u k ) k ∈U NJ(x, z, u).

Definition 5 . 1 .

 51 For a given depth p ∈ N, we define the asymptotic branching factor m, as the smallest real m ∈ [1, M ] such that ∃R ≥ 1, ∀p ≥ 0, |Υ * p | ≤ Rm p .

1 N N - 1 k=0 u 1

 111 g(x) := 0.4-x-(-2, 0.5) ∞ min 0.2-|x 1 +2.5|, 1-|x 2 +1| and g f (x) := x ∞ -0.1.For a given N , the discrete control problem becomes:V (x) = inf ,k with u = ((u 1,k , u 2,k )) k ∈ U N , g(y u k ) ≤ 0 for k = 0, • • • , N, and Ψ(y u N ) ≤ 0 ,where (y u k ) k is the discrete state variable, corresponding to the control policy u ∈ U N , and starting at the initial position x. Henceforth, the discrete auxiliary value function isW (x, z) = inf u∈U N

Figure 14 displays

 14 Figure 14 displays optimal trajectories obtained from three different initial positions. A simultaneous optimistic planning algorithm is used for this simulation with I max = 3200 and N = 40. The optimal controls are displayed on Figure 15.

Figure 14 :

 14 Figure 14: (Test 5) An optimistic planning approach with N = 40. Optimal trajectories corresponding to three different initial data.

Figure 15 :

 15 Figure 15: (Test 5) An optimistic planning approach with N = 40. Control laws corresponding to three different initial data.

ct.gitlabpages.inria.fr/gallery

hampath.org

ct.gitlabpages.inria.fr/gallery

The code is available and executable online at ct.gitlabpages.inria.fr/gallery

ct.gitlabpages.inria.fr/gallery

[6] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl. CasADi -A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 2018.

More explicitly, the condition [START_REF] Hu | A discontinuous galerkin finite element method for hamilton-jacobi equations[END_REF] means that the matrix D 2 u (t, x) -k(t)I is negative semidefinite, that is, that all eigenvalues of D 2 u are bounded from above by k(t). Numerically, the semi-concave stability is enforced as a uniform upper bound on the second directional incremental ratios. The core of the theory is an abstract result of convergence for perturbed semi-concave stable solutions. Theorem 4.9. Consider problem (59) for a semi-concave initial condition v 0 with compact support, and assume the family v is semi-concave stable. Define the truncation error associated to v as

Then, for any t ∈ [0, T ],

The second ingredient of the Lin-Tadmor theory, i.e., the L 1 estimation of the truncation error, requires a reversed approach in measuring the truncation error, that is, to plug the numerical solution into the exact equation. An easier expression can be derived for Godunov-type schemes, taking into account that in this case numerical errors are generated only by the projection step, and not by the evolution operator. In this case, a practical application of this theory is presented in [START_REF] Lin | l 1 -stability and error estimates for approximate Hamilton-Jacobi solutions[END_REF].

Higher-order FD schemes. The basic strategy for constructing high-order finite difference methods has been first proposed in [START_REF] Osher | High essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF], and uses an intermediate semi-discrete formulation. In a successive step, a TVD Runge-Kutta method is applied to obtain a fully discrete approximation. The semi-discrete scheme is constructed using a monotone numerical Hamiltonian, in the

in which H(•, •) is increasing with respect to its first argument, and decreasing with respect to the second. In ( 63), D ± j [V ] denote high order approximations of the right/left derivative at the node x j , which replace in the numerical Hamiltonian the mere right/left incremental ratios. In the most extensively studied versions of the scheme, these estimates are usually obtained via non-oscillatory (ENO/WENO) techniques (see [START_REF] Osher | High essentially nonoscillatory schemes for Hamilton-Jacobi equations[END_REF][START_REF] Bryson | High-order semi-discrete central-upwind schemes for multi-dimensional hamilton-jacobi equations[END_REF][START_REF] Jiang | Weighted eno schemes for hamilton-jacobi equations[END_REF][START_REF] Kurganov | Semidiscrete central-upwind schemes for hyperbolic conservation laws and hamilton-jacobi equations[END_REF][START_REF] Kurganov | Adaptive central-upwind schemes for hamilton-jacobi equations with nonconvex hamiltonians[END_REF]).

Higher-order SL schemes. The SL scheme ( 58) is easily extended to a higher consistency rate by replacing the P 1 space interpolation I 1 with an interpolation of higher accuracy [START_REF] Falcone | Semi-lagrangian schemes for hamilton-jacobi equations, discrete representation formulae and godunov methods[END_REF][START_REF] Carlini | A weighted essentially nonoscillatory, large time-step scheme for hamilton-jacobi equations[END_REF]. In general, since characteristics are not straight lines, a more accurate method of characteristics tracking is also desirable [START_REF] Falcone | Discrete time high-order schemes for viscosity solutions of hamiltonjacobi-bellman equations[END_REF].

In some model cases (a single space dimension, no x-dependence of the Hamiltonian) convergence of high-order SL schemes, for both the evolutive and the stationary case, can be proved by showing their quasi-monotonicity (see [START_REF] Ferretti | Convergence of semi-lagrangian approximations to convex hamilton-jacobi equations under (very) large courant numbers[END_REF][START_REF] Falcone | Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations[END_REF][START_REF] Bokanowski | Value iteration convergence of"-monotone schemes for stationary hamilton-jacobi equations[END_REF]).

Discontinuous Galerkin. The application to HJ equations of Discontinuous Galerkin (DG) methods uses the relationship with conservation laws. In fact, what is discretized in this case