Gabriel Velho

Jean Auriol

Riccardo Bonalli

A Gradient Descent-Ascent Method for Continuous-Time Risk-Averse Optimal Control

Keywords: Gradient descent-ascent, Coherent risk measures, Stochastic optimal control, Non-linear control, Risk-averse control

In this paper, we consider continuous-time stochastic optimal control problems where the cost is evaluated through a coherent risk measure. We provide an explicit gradient descent-ascent algorithm which applies to problems subject to non-linear stochastic differential equations. More specifically, we leverage duality properties of coherent risk measures to relax the problem via a smooth min-max reformulation which induces artificial strong concavity in the max subproblem. We then formulate necessary conditions of optimality for this relaxed problem which we leverage to prove convergence of the gradient descent-ascent algorithm to candidate solutions of the original problem. Finally, we showcase the efficiency of our algorithm through numerical simulations involving trajectory tracking problems and highlight the benefit of favoring risk measures over classical expectation.

Introduction

Risk-averse stochastic optimal control is a powerful tool for designing control laws that feature robustness against uncertainties. Relevant applications of this theory range from safe financial investment to safe control of autonomous systems, as evidenced by the monographs [START_REF] Shapiro | Lectures on Stochastic Programming: Modeling and Theory[END_REF][START_REF] Chapman | Risk-Sensitive Safety Analysis Using Conditional Valueat-Risk[END_REF] and references therein. For example, when controlling non-linear stochastic systems in robotics or space applications, robustness against uncertainties is often leveraged to ensure safety during motion. The problem also arises in finance, where, to mitigate unrecoverable losses, an agent may want to invest safe assets uniquely, even in exchange for a smaller average gain. This is in contrast with "classical" stochastic optimal control, which rather aims at minimizing the expectation of a stochastic cost. Although high reward on average might be attained, these risk-neutral methods result in optimal controls that may generate trajectories that may cause catastrophic losses with small, though non-zero, probability. As such, these solutions are considered unacceptable in terms of safety in many applications [START_REF] Wang | Risk-averse autonomous systems: A brief history and recent developments from the perspective of optimal control[END_REF].

To safely control systems whose dynamics are subject to uncertainties, various methods have been developed. A first possible approach is robust control, where optimal control problems are directly solved in the worst-case scenarios [START_REF] Kothare | Robust constrained model predictive control using linear matrix inequalities[END_REF]. Although robust control theoretically guarantees resilience to unfavorable situations, it can be over-conservative and does not apply to dynamics with unbounded uncertainties, e.g., stochastic differential equations [START_REF] Wang | Risk-averse autonomous systems: A brief history and recent developments from the perspective of optimal control[END_REF]. A more flexible approach is risk-averse control, which aims at minimizing a suitable risk measure of the stochastic cost, offering a good trade-off between minimizing the expected cost and mitigating catastrophic scenarios. Many risk-averse control methods involve discretizing the dynamics before proceeding with optimization [START_REF] Chapman | Risk-Sensitive Safety Analysis Using Conditional Valueat-Risk[END_REF]. For discrete problems, risk-averse Model Predictive Control (MPC) is a popular choice [START_REF] Chow | A framework for time-consistent, riskaverse model predictive control: Theory and algorithms[END_REF][START_REF] Sopasakis | Risk-averse riskconstrained optimal control[END_REF]. It operates similarly to standard MPC but optimizes a risk measure instead of the expectation of the cost. Continuous-time risk-averse control methods are less common. Although risk-averse analogs of the Hamilton-Jacobi-Bellman (HJB) equation [START_REF] Ruszczyński | A Risk-Averse Analog of the Hamilton-Jacobi-Bellman Equation[END_REF] and the Pontryagin Maximum Principle (PMP) [8] have been established, these methods do not directly provide tractable numerical schemes.

Designing a risk-averse method that scales with many different risk measures is a challenging task. Therefore, the focus has mostly been on developing methods that minimize specific risk measures. A significant hindrance of risk measures is their lack of time consistency, meaning that the optimal solution at the current time may not solely depend on the current state's values [START_REF] Shapiro | On a time consistency concept in risk averse multistage stochastic programming[END_REF]. A class of time-consistent risk measures exists for which dynamic programming can be employed to solve the associated Markov Decision Process (MDP) [START_REF] Ruszczyński | Risk-averse dynamic programming for Markov decision processes[END_REF][START_REF] Ruszczyński | Risk-Averse Control of Undiscounted Transient Markov Models[END_REF]. Nevertheless, such methods do not apply to the widely adopted Conditional Value at Risk (CV@R), a risk measure that is not time consistent [START_REF] Shapiro | Lectures on Stochastic Programming: Modeling and Theory[END_REF]. To address this issue, a state-space lifting approach has been proposed in [START_REF] Bäuerle | Markov Decision Processes with Average-Valueat-Risk criteria[END_REF]. Although effective, these approaches are limited to discretetime and finite discrete space MDP. A continuous-time and continuous space approach for the CV@R has been proposed in [START_REF] Miller | Optimal Control of Conditional Value-at-Risk in Continuous Time[END_REF], where a representation of the CV@R is used to re-write the problem into a double minimization that is solved with a gradient descent approach. Another important class of risk measures is the family of risk measures that can be expressed as a linear function of the state's probability distribution [START_REF] Shapiro | Lectures on Stochastic Programming: Modeling and Theory[END_REF]Chapter 6]. To minimize these risk measures, a mean field game approach can be applied. In [START_REF] Pfeiffer | Risk-averse Merton's Portfolio Problem[END_REF] and [START_REF] Laurent | Optimality conditions in variational form for non-linear constrained stochastic control problems[END_REF], the Frank-Wolf algorithm is used to solve the problem, with proof of convergence properties to a necessary optimality condition. Yet, these approaches work with very specific risk measures.

In this paper, we propose a gradient-based method to solve continuoustime, non-convex risk-averse control problems where the corresponding risk measure is only coherent, a fairly general requirement encompassing all the aforementioned risk measures. We will only consider the case of a diffusion independent of the control to keep the computation of gradients tractable. This is a reasonable assumption for applications such as aerospace and robotics, where the diffusion process often models external perturbations (e.g., gusts of wind on a robot), which are then uncorrelated from the control. The motivation to pursue this research direction came from the fact that gradient descent methods have already seen a lot of success in the more classical setting of average-oriented minimization [START_REF] Archibald | A Stochastic Gradient Descent Approach for Stochastic Optimal Control[END_REF], non-stochastic control [START_REF] Konolige | A gradient method for realtime robot control[END_REF], as well as some risk-averse methods [START_REF] Miller | Optimal Control of Conditional Value-at-Risk in Continuous Time[END_REF][START_REF] Laurent | Optimality conditions in variational form for non-linear constrained stochastic control problems[END_REF][START_REF] Garreis | An Interior-Point Approach for Solving Risk-Averse PDE-Constrained Optimization Problems with Coherent Risk Measures[END_REF]. Minimizing directly the risk measure is challenging. We, therefore, leverage the duality properties of coherent risk measures to transform the minimization problem into a min-max problem of a smooth, non-convex/linear function. Min-max optimization has been widely studied in the convex-concave setting [START_REF] Nedić | Subgradient Methods for Saddle-Point Problems[END_REF]. Still, recent applications to machine learning (particularly Generative Adversarial Networks) have pushed studies toward the non-convex / strongly concave setting [START_REF] Mokhtari | A Unified Analysis of Extra-gradient and Optimistic Gradient Methods for Saddle Point Problems: Proximal Point Approach[END_REF]. The gradient descent-ascent method has proven to be a successful technique for addressing these problems in finite-dimensional spaces. In [START_REF] Lin | On gradient descent ascent for nonconvex-concave minimax problems[END_REF], the proposed algorithm reaches an ϵ distance from a critical point within a finite time frame. In [START_REF] Chen | Proximal Gradient Descent-Ascent: Variable Convergence under K L Geometry[END_REF], the algorithm's convergence to the critical point is guaranteed, given a specific geometric property of the cost function.

Our contribution consists of applying gradient descent ascent to solve the min-max problem, which stems from dual representations of coherent risk measures, by extending the work in [START_REF] Chen | Proximal Gradient Descent-Ascent: Variable Convergence under K L Geometry[END_REF]. This is not an easy task in that our problem is not smooth nor strongly concave. To address such an issue, we propose to apply gradient descent-ascent to a modified version of the problem which is smooth, non-convex / strongly concave, and arbitrarily "close" to the original formulation. Via appropriate proofs of convergence, we then show that the gradient of the cost converges to 0. By requiring more regularity in the dynamics of the system, we can also prove that the sequence of controls generated by our algorithm converges (in an appropriate sense) to the set of critical points of the original problem. We finally implement the algorithm on a realistic robotic application to showcase the efficiency of our easy-to-use method.

The paper is organized as follows. In Section 2, we begin by providing a brief summary of some basic concepts of stochastic calculus and risk measures. Section 3 outlines the specific settings of dynamics and control strategies under consideration. Next, we present in Section 4 our algorithm, providing details for the computation of the gradient. In Section 5, we prove various convergence properties under different assumptions on the dynamics. Finally, we analyze in Section 6 numerical convergence, and we showcase the effectiveness of our approach through a series of simulations.

Notations and Preliminary Results

In this section, we briefly present some basic concepts and results of stochastic calculus we need in our work. More details can be found in [START_REF] Gall | Brownian Motion, Martingales, and Stochastic Calculus[END_REF] and [START_REF] Yong | Stochastic controls: Hamiltonian systems and HJB equations[END_REF]Chapter 1.6]. From now on, we fix n, m, d ∈ N, a finite time horizon T > 0, and p ∈ [1, +∞). Let E be a vector space. We denote E * its dual. As a result of Riesz theorem, if E is also a Hilbert space, then it is isomorphic to its dual space E * . Consequently, throughout the paper, elements in (L 2 ) * are implicitly assumed as elements in L 2 . We denote ∥ • ∥ as the Euclidian norm. Other norms will be properly introduced as they are needed. We consider random variables in a probability space (Ω, G, P). For any subsigma algebra S ⊂ G, we denote by L p S (Ω, R n ) the Banach space of random variables z : Ω → R n that are S-measurable and satisfy:

∥z∥ L p ≜ E [∥z∥ p ] 1/p < ∞. It is a standard consequence of Riesz theorem that L p S (Ω, R n ) * ∼ = L q S (Ω, R n ), with q the conjugate exponent of p. Let (W t ) t∈[0,T ] be a d-dimensional Wiener process, generating a filtration F ≜ (F t ) t∈[0,T ] = (σ(W s , 0 ≤ s ≤ t)) t∈[0,T ] .
We denote by L p F ([0, T ]×Ω, R n ) the Banach space of processes x : [0, T ]×Ω → R n which are progressive with respect to F and such that:

∥x∥ L p F ≜ E T 0 ∥x∥ p (t)dt 1/p < ∞.
In addition, we denote by C p F ([0, T ] × Ω, R n ) the Banach space of F-adapted processes x which have continuous sample paths and finite sup norm

∥x∥ C p F ≜ E sup 0≤t≤T ∥x∥ p (t) 1/p < ∞.
In particular,

C p F ([0, T ] × Ω, R n ) ⊂ L p F ([0, T ] × Ω, R n ).
In what follows, we often use the notation

x(t), t ∈ [0, T ] to denote progressive processes. A F- adapted process x : [0, T ] × Ω → R n such that x(t) ∈ L 1 F T for every t ∈ [0, T ] is called a martingale provided that E[x(t)|F s ] = x(s),
for all 0 ≤ s ≤ t ≤ T . We then say that a martingale x :

[0, T ] × Ω → R n is uniformly bounded in L p F ([0, T ] × Ω, R n ) if there exists a constant C > 0 such that sup t∈[0,T ] ∥x(t)∥ L p ≤ C. For every x = (x 1 | . . . |x d ) ∈ L p F ([0, T ] × Ω, R n×d ), we write
y(t) = T 0 x(t)dW t ≜ d i=1 T 0 x i (t)dW i t ,
for the Itô integral of x with respect to W , and recall that y is then a martingale in

C p F ([0, T ] × Ω, R n ).
We also recall the Burkholder-Davis-Gundy inequality for stochastic integrals:

E sup 0≤t≤T ∥y(t)∥ p ≤ C p E   T 0 ∥x(t)∥ 2 dt p/2   , ( 1
)
where C p is a constant which depends on p uniquely. As outlined in the introduction, our work is centered on controlling the extreme values or tail of the distribution of a stochastic cost rather than controlling its expectation uniquely. To estimate these extreme values, we leverage risk measures. Unlike expectation, risk measures are not always linear or Fréchet-differentiable. However, they possess several interesting properties we discuss below. We refer to [START_REF] Shapiro | Lectures on Stochastic Programming: Modeling and Theory[END_REF]Chapter 6] for all the results presented below on risk measures. Definition 2.1. A finite risk measure ρ is a mapping from a space Z = L p (Ω, R) with p ∈ [1, +∞) to R. A finite risk measure is said to be coherent if it verifies the following properties: Convexity:

∀(Z, Z ′ ) ∈ Z ∀λ ∈ [0, 1], ρ(λZ + (1 -λ)Z ′ ) ≤ λρ(Z) + (1 -λ)ρ(Z ′ ). Monotonicity 1 : ∀(Z, Z ′ ) ∈ Z, Z ≤ Z ′ ⇒ ρ(Z) ≤ ρ(Z ′ ).
Translation equivariance:

∀Z ∈ Z, ∀a ∈ R, ρ(Z + a) = ρ(Z) + a.
Positive homogeneity:

∀Z ∈ Z, ∀t ∈ R + , ρ(tZ) = tρ(Z).
Definition 2.1 yields some interesting properties that are used throughout our work.

Theorem 2.2 ([1]

). A finite coherent risk measure ρ : Z → R verifies the following properties:

1. The convex subdifferential of ρ at 0, denoted by ∂ρ(0), is a nonempty closed and bounded subset of Z * . It is, in particular, weakly-* compact. 2. For every Z ∈ Z, the risk measure can be represented as

ρ(Z) = max ζ∈∂ρ(0) E[ζZ].

For every Z ∈ Z, ∂ρ(Z) is a nonempty closed and bounded subset of

Z * , and it can be expressed as

∂ρ(Z) = arg max ζ∈∂ρ(0) E[ζZ]. 1 for (Z, Z ′ ) ∈ Z, it is said that Z ≤ Z ′ if P({ω : Z(ω) ≤ Z ′ (ω)}) = 1.
4. For every Z, H ∈ Z, the mapping ρ is Hadamard directionally (sub-) differentiable at Z along H, and the differential can be expressed as

ρ ′ (Z, H) = max ζ∈∂ρ(Z) E[ζH].
The most commonly adopted risk measure is the Conditional-Value-at-Risk, denoted by CV@R, which we will regularly consider in our examples. Intuitively, CV@R represents the average of the highest values that a random variable can attain. Consequently, it is a highly valuable metric for assessing potential losses in a worst-case scenario. Definition 2.3. For Z in L 2 (Ω, R) and α ∈ (0, 1], the Conditional-Valueat-Risk of parameter α of Z is given by:

CV@R α [Z] = inf t + E[max(0, Z -t)] α : t ∈ R .
The parameter α sets up the number of "worst-case scenarios". For instance, setting α = 0.05 means that CV@R α [Z] corresponds to the average of the top 5% values that Z can attain. Note that the subdifferential at 0 of the Conditional-Value-at-Risk can be explicitly calculated, which is extremely helpful for numerical simulations.

Proposition 2.4. The subdifferential of CV@R α at 0 is given by:

∂(CV@R α )(0) = {ζ ∈ L 2 (Ω, R) : E[ζ] = 1, ζ ∈ [0, α -1 ] a.s.}.
(2)

Problem Formulation

In what follows, we detail the settings in which we study the controlled stochastic dynamics. We then devise the optimization problem and present some necessary conditions for optimality an optimal solution must satisfy. As already mentioned, we focus on two different stochastic dynamics:

1. Setting A. A general non-linear stochastic system with no control term in the diffusion, which is coherent with our framework of modeling exterior perturbations affecting a mechanical system. The proposed algorithm can be run when considering such dynamics, although only partial theoretical guarantees of convergence may be derived. Nevertheless, through numerical simulations, we show that reliable solutions are often generated.

Setting B.

A control affine stochastic system with no control term in the diffusion. This system is a particular case of the general dynamics which are often encountered in many applications, such as aerospace and robotics [START_REF] Trélat | Contrôle optimal : théorie & applications[END_REF][START_REF] Bonalli | Solving optimal control problems for delayed control-affine systems with quadratic cost by numerical continuation[END_REF][START_REF] Bonalli | Optimal Control of Endoatmospheric Launch Vehicle Systems: Geometric and Computational Issues[END_REF][START_REF] Bonalli | Optimal Control of Aerospace Systems with Control-State Constraints and Delays[END_REF]. Importantly, under such dynamics, it is possible to theoretically prove the convergence of our algorithm.

From now on, we assume state variables take values in R n , whereas control variables take values in R m . For the sake of clarity in the exposition and without loss of generality, from now on, we assume the dynamics are perturbed through a weighted one-dimensional Wiener process W , as below.

Setting A: The general non-linear control system

In this general setting, we consider non-linear stochastic dynamics in the form

dX u (t) = b(t, X u (t), u(t))dt + σ(t, X u (t))dW t , X u (0) = X 0 , ( 3
)
where X 0 ∈ R n is a deterministic initial condition (note that extending our results where the initial condition is a random variable is straightforward).

In particular, we consider diffusion terms σ independent from the control input. This is generally the case in many applications such as aerospace and robotics [START_REF] Bonalli | Sequential convex programming for non-linear stochastic optimal control[END_REF]. We also consider the following classical assumptions on the mappings b :

[0, T ] × R n × R m → R n and σ : [0, T ] × R n → R n to hold true: 1. The mappings b(•, x, u) : [0, T ] → R n σ(•, x) : [0, T ] → R n are continuous for every (x, u) ∈ R n × R m . 2. For all t in [0,T], the mappings b(t, •, •) : R n × R m → R n σ(t, •) : R n → R n are differentiable with bounded L-Lipschitz gradient. That is, there exists a constant L > 0 such that ∂b ∂x (t, x, u) + ∂b ∂u (t, x, u) + ∂σ ∂x (t, x) ≤ L, ∂b ∂x (t, x, u) - ∂b ∂x (t, y, v) ≤ L(|x -y| + |u -v|), ∂b ∂u (t, x, u) - ∂b ∂u (t, y, v) ≤ L(|x -y| + |u -v|), ∂σ ∂x (t, x) - ∂σ ∂x (t, y) ≤ L|x -y|.
In particular, given a progressive control input u ∈ L 2 F ([0, T ] × Ω, R m ), the assumptions above guarantee the existence and uniqueness of a solution 3), see, e.g., [START_REF] Yong | Stochastic controls: Hamiltonian systems and HJB equations[END_REF]Chapter 1.6]. The additional assumption on Lipschitz gradients, in particular, will be required to compute the Fréchet derivative of X u with respect to u (see Section 4.2).

X u ∈ C 2 F ([0, T ] × Ω, R n ) to (
We consider cost functions of the form

J(u) ≜ T 0 f (t, X u (t), u(t))dt + g(T, X u (T )), with X u solution of (3),
where f and g are positive continuous functions, such that f (t, •, •) and g(t, .) are C 1 with Lipschitz gradient in the state and control variables for all t in [0, T ]. Note that these functions cover a variety of costs that are commonly used in many applications. For instance, in trajectory tracking f consists of the quadratic difference between a reference trajectory and the output trajectory [START_REF] Trélat | Contrôle optimal : théorie & applications[END_REF].

Setting B: control affine non-linear systems

In this setting, we consider the particular case of dynamics of the form

dX u (t) = [b 0 (t, X u (t)) + m i=0 u i (t)b i (t, X u (t))] dt + σ(t, X u (t))dW t , X u (0) = X 0 , (4) where the mappings b i : [0, T ] × R n → R n and σ : [0, T ] × R n → R n are C 1
with Lipschitz gradient in the state variables for all t in [0, T ]. Again, these assumptions ensure existence and uniqueness of [START_REF] Kothare | Robust constrained model predictive control using linear matrix inequalities[END_REF]. We additionally require the functions b i to be uniformly bounded, which still allows for a wide variety of dynamics. Finally, for this setting, we consider slightly less generic stochastic cost functions of the form

X u ∈ C 2 F ([0, T ] × Ω, R n ) solution to
J(u) ≜ C u ∥u∥ 2 U + T 0 f (t, X u (t))dt + g(T, X u (T )),
with X u solution of ( 4),

where f and g follow the same assumptions as above. Additionally, we require that f and g satisfy

∂(f, g) ∂x (t, x) - ∂(f, g) ∂x (t, y) ≤ L∥x -y∥, ∂(f, g) ∂x (t, 0) ≤ L,
for some constant L > 0. These assumptions, for instance, enable the use of quadratic functionals. Note that such cost functions are still relevant for many applications, including trajectory tracking-type problems [START_REF] Trélat | Contrôle optimal : théorie & applications[END_REF].

Optimization problem

We seek optimal controls in the space

U = L 2 ([0, T ], R m ) × R q for some fixed q ∈ N, which parametrizes feedback controls in L 2 F ([0, T ] × Ω, R m
). More specifically, to each control parameter u = (v, λ) ∈ U we associate a feedback control in the following manner:

u v,λ (t, x) = v(t) + K λ (x), (5)
where K λ : R n → R m is a vector-valued function that is smoothly parametrized by the vector λ. We assume that

∀x, y ∈ R n , ∀λ ∈ R q , ∥K λ (x) -K λ (y)∥ ≤ C(1 + ∥λ∥)∥x -y∥ ∀x ∈ R n , ∀λ, µ ∈ R q , ∥K λ (x) -K µ (x)∥ ≤ C(1 + ∥x∥)∥λ -µ∥. ( 6
)
We also assume the function (x, λ) → K λ (x) has a lipschitz gradient, meaning

∥∇ λ,x K λ (x) -∇ λ,x K µ (y)∥ ≤ L(∥x -y∥ + ∥λ -µ∥) (7)
for some constant L > 0. In our implementations, we use a feedback law that is linear in λ and x, which therefore automatically verifies both assumptions, but more generic choices are clearly possible. We denote by

X v,λ ∈ C 2 F ([0, T ] × Ω, R n ) the unique solution to (3) (or ( 4
)
) which is generated by the control u v,λ . These control laws will be composed of a deterministic term to steer the trajectory and a stochastic term to compensate for the diffusion. This choice of controls is motivated by the fact that feedback-type controls are easy to implement and are widely used in practice for risk-averse control [START_REF] Wang | Risk-averse autonomous systems: A brief history and recent developments from the perspective of optimal control[END_REF] 2 . One easily verifies that under [START_REF] Sopasakis | Risk-averse riskconstrained optimal control[END_REF], the controlled systems under these controls still satisfy the assumptions of Setting A or Setting B. In what follows, we endow the space U with the scalar product

⟨(v 1 , λ 1 ), (v 2 , λ 2 )⟩ U ≜ T 0 v 1 (t) • v 2 (t)dt + λ 1 • λ 2 .
From now on, since clear from the context we equivalently denote u

= (v, λ) = u v,λ and X u = X v,λ = X u v,λ for simplicity.
Let us now state and reformulate our optimal control problem as an infinite-dimensional min-max optimization problem using the coherent measure properties described in Theorem 2.2. In particular, let ρ be a finite coherent risk measure defined on L 1 F T (Ω, R). The main objective of this paper is to find a local optimum of the problem

min u∈U ρ(J(u)). ( 8
)
This problem is not easily solvable as ρ is not always Fréchet differentiable. Given the properties of coherent risk measures we listed previously, we can rewrite (8) via the following min-max formulation:

min u∈U max ζ∈∂ρ(0) E[ζJ(u)]. (9)
Let us remark that ∂ρ(0

) is a subset of L ∞ F T (Ω, R) ⊂ L 2 F T (Ω, R).
Using this formulation, we can leverage saddle point algorithms to solve problem (8). Saddle point algorithms have been extensively investigated for min-max optimization, motivating our approach [START_REF] Nedić | Subgradient Methods for Saddle-Point Problems[END_REF][START_REF] Mokhtari | A Unified Analysis of Extra-gradient and Optimistic Gradient Methods for Saddle Point Problems: Proximal Point Approach[END_REF][START_REF] Lin | On gradient descent ascent for nonconvex-concave minimax problems[END_REF][START_REF] Chen | Proximal Gradient Descent-Ascent: Variable Convergence under K L Geometry[END_REF].

Before introducing our method to solve (9) (and therefore (8)), we characterize its optimal solutions through appropriate necessary conditions for optimality.

Necessary conditions for optimality

In what follows, we assume our problem is in the most general Setting A. As of common use in min-max optimization settings, we propose to solve Problem (9) by finding extremal points [START_REF] Benzi | Numerical solution of saddle point problems[END_REF]. For this, let us first introduce appropriate first-order necessary conditions for optimality which are satisfied by any optimal solution u * to (8). These conditions involve the differential of the cost function J that we first need to establish clearly.

U denoted ∇J(u, ω) such that ∀h ∈ U, D u J(u, ω) • h = ⟨∇J(u, ω), h⟩ U .
The proof can be found in Appendix A.1, along with an explicit formula of ∇J.

Additional notations:

In what follows, we call gradient of J the function ∇J(•, •), that we see as a mapping ∇J : U → L 2 F T (Ω, U). We may need to differentiate the two components for clarity, in which case we will write

∇J = (∇ v J, ∇ λ J), with ∇ v J(u) ∈ L 2 F T (Ω, L 2 ([0, T ], R m )) and ∇ λ J(u) ∈ L 2 F T (Ω, R q ). Additionally, for an element z in L 2 F T (Ω, U) and u in U, we will denote by ⟨z, u⟩ U the random variable in L 2 F T (Ω, R) such that ∀ω ∈ Ω, ⟨z, u⟩ U (ω) = ⟨z(ω), u⟩ U .
In particular, we can "invert" the expectation and the scalar product as follows:

E[⟨z, u⟩ U ] = ⟨E[z], u⟩ U . ( 10
)
With these properties in mind, we can establish a necessary condition any optimal solution to (9) must verify as follows:

Proposition 3.2. Let ρ be a finite coherent risk measure and let u * be an optimal solution to (8). There exists

ζ * ∈ ∂ρ(0) such that ⟨E[ζ * ∇J(u * )], h⟩ U = 0, for every h ∈ U. ( 11
)
The proof of Proposition 3.2 is quite similar to the proof given in [START_REF] Kouri | Existence and Optimality Conditions for Risk-Averse PDE-Constrained Optimization[END_REF] or [8]. Nevertheless, we report the full proof of Proposition 3.2 in Appendix A.2 for the sake of completeness.

The modified problem

Seeking extremals of min-max problems is generally a hard task, and not that many algorithms which enjoy convergence guarantees exist in the literature. Notably, for finite-dimensional problems with strongly concave costs, some convergence results are provided in [START_REF] Chen | Proximal Gradient Descent-Ascent: Variable Convergence under K L Geometry[END_REF]. Inspired by this latter work, we propose an algorithm that works on a modified version of problem [START_REF] Shapiro | On a time consistency concept in risk averse multistage stochastic programming[END_REF] whose cost is strongly concave and whose solutions are arbitrarily close to the solutions of the original problem. More specifically, the modified problem takes the following form:

min u∈U max ζ∈∂ρ(0) E[(ζ -γζ 2 )J(u)], (12)
where γ > 0 is a "concavifying" parameter, which controls how close ( 12) is to the original problem. Note that the term

E[ζ 2 J(u)] is well defined as ζ ∈ ∂ρ(0) ⊂ L ∞ F T (Ω, R).
In Section 5, we discuss in what sense the solutions to [START_REF] Bäuerle | Markov Decision Processes with Average-Valueat-Risk criteria[END_REF] converge to solutions to (9) for γ → 0, endowing our methodology with guarantees of success. Also, from a numerical point of view, we show later that γ > 0 may be chosen in the order of the machine precision.

Ascent-Descent Algorithm

To solve the min-max problem (9), we first introduce in this section a gradient descent-ascent-based approach to solve the modified problem [START_REF] Bäuerle | Markov Decision Processes with Average-Valueat-Risk criteria[END_REF]. A gradient descent-ascent is an iterative optimization algorithm that can be used to find saddle points. Descent-ascent algorithms represent arguably sound paradigms to seek extremals for general min-max problems. Importantly, under appropriate assumptions, convergence to extremals may be theoretically proven, as we will show later. Until further notice, we start by considering the setting of general non-linear dynamics (Setting A).

Update rule

The gradient descent-ascent algorithm is an iterative optimization algorithm for min-max problems. At each iteration, it computes the Riesz representation of the gradient of the function being maximized and the gradient of the function being minimized. It then updates the estimates by moving in the direction of their respective gradients. With obvious notations, our gradient descent-ascent algorithm updates at each iteration n our estimates ζ n and (v n , K λn ) with the following rule:

v n+1 = v n -ηE[ζ n ∇ v J(v n , K λn )], λ n+1 = λ n -ηE [ζ n ∇ λ J(v n , K λn )] , ζ n+1 = P ∂ρ(0) (ζ n + βJ(v n , K λn )), (13)
where P ∂ρ(0) is the projection on the closed convex set ∂ρ(0), and η and β are two positive constants that represent the step size of the gradient. To increase convergence speed, one could make these steps vary from iteration to iteration based on the change of direction of the gradient [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF]. We will keep them constant for the sake of simplicity. All convergence results should be valid for a range of step sizes that are both upper and lower bounded, although proof of this improvement is left for future work.We recall that ∂ρ(0) can be viewed as a subset of L ∞ F T (Ω, R), which justifies the above mathematical operations between ζ n and the other terms. Note that it is important to project ζ n + βJ(u n ) onto ∂ρ(0), as there is no a priori guarantee that ζ n + βJ(u n ) stays inside ∂ρ(0).

Unfortunately, through the scheme [START_REF] Miller | Optimal Control of Conditional Value-at-Risk in Continuous Time[END_REF], no clear convergence result is available. To develop theoretical guarantees of convergence, we propose to rather leverage the following alternative rule, which essentially represents an update rule for Problem [START_REF] Bäuerle | Markov Decision Processes with Average-Valueat-Risk criteria[END_REF]:

v n+1 = v n -ηE[ζ n (1 -γζ n )∇ v J(v n , K λn )], λ n+1 = λ n -ηE [ζ n (1 -γζ n )∇ λ J(v n , K λn )] , ζ n+1 = P ∂ρ(0) (ζ n + β(1 -2γζ n )J(v n , K λn )). ( 14
)
The algorithms described by ( 13) and ( 14) offer the advantage of being easily implementable, provided explicit expressions of the gradient ∇J are available. We provide such expressions in the following section.

Computation of gradients

The practical implementation of ( 13) and ( 14) requires computing the gradient of the cost, and in particular of the trajectory, with respect to both the variables ζ and u. Below, we summarize the expressions of this gradient, referring the reader to Appendix A.1 for their proof. Let us clarify that we denote K ′ λ (x) the derivative of K λ (x) with respect to x.

∇ v J(v, λ)(t) = I u (t)ϕ -1 u (t)∂ u b(t, X u (t), u(t)) ∇ λ J(v, λ) = T 0 I u (s)ϕ -1 u (s)∂ u b(s, X u (s), u(s))∇ λ K λ (X u (s))ds.
where

I u (s) ≜ ∂ x g(T, X u (T ))ϕ u (T ) + T s [∂ x f (t, X u (t), u(t)) + ∂ u f (t, X u (t), u(t))K ′ λ (X u (t))]ϕ u (t)dt.
Here, ϕ u (.) is the fundamental matrix solution to the linearized problem, i.e., the unique solution to

     dϕ(t) = [∂ x b(t, X u (t), u(t)) + ∂ u b(t, X u (t), u(t))K ′ λ (X u (t))]ϕ(t)dt + ∂ x σ(t, X u (t))ϕ(t)dW t , ϕ(0) = Id n . (15)
Thus, the primary challenge of implementing this method lies in simulating the trajectories of X u and ϕ u and computing expectations of random variables that involve terms from those trajectories. To deal with these technical details, we used a classical approach presented in [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]Chapter 10] and refer to it for more extensive elaboration. The contents of this section furnish the necessary components to implement the algorithm and evaluate its performance for the general Setting A.

In the next section, we establish some convergence results that may be of particular significance for Setting B, as they provide insight into the convergence of the algorithm towards a critical point that satisfies the necessary optimality condition [START_REF] Ruszczyński | Risk-Averse Control of Undiscounted Transient Markov Models[END_REF].

Convergence Results

We leverage the method introduced in [START_REF] Chen | Proximal Gradient Descent-Ascent: Variable Convergence under K L Geometry[END_REF] to examine the convergence of the smooth non-convex/strong concave gradient descent-ascent algorithm in the most general Setting A. In what follows, we adapt the approach to an infinite-dimensional space. With some minor adjustments, we show that the gradient of the cost converges to 0 and that the cost function converges to a finite value. While the study in [START_REF] Chen | Proximal Gradient Descent-Ascent: Variable Convergence under K L Geometry[END_REF] also establishes that in finitedimensional spaces, every converging subsequence of the control estimates tends to a critical point, we were unable to prove this result in our general setting. Nevertheless, in the slightly more restrictive Setting B, any weakly converging subsequence of our control estimates converges to a point that is a "γ-critical point." That is, a control u γ that verifies

∃ζ, ⟨E[ζ ∇J(u)], h⟩ U = ∥h∥ U O(γ) ∀h ∈ U. ( 16
)
This important result shows that by progressively selecting smaller values of γ, we can obtain a control that verifies the necessary optimality condition [START_REF] Ruszczyński | Risk-Averse Control of Undiscounted Transient Markov Models[END_REF]. In addition, taking γ → 0 yields critical points for the original problem, as we state at the end of this section.

Convergence of the gradient in Setting A

Let (u n , ζ n ) be the estimates given by algorithm [START_REF] Pfeiffer | Risk-averse Merton's Portfolio Problem[END_REF], in what follows, we prove the convergence of ∥u n+1 -u n ∥ U to 0 in Setting A, as well as the convergence of the γ-approximated cost, denoted Φ γ : U → R and defined as

Φ γ (u) ≜ max ζ∈∂ρ(0) E[ζ(1 -γζ)J(u n )].
The proof of the convergence closely follows the one presented in [START_REF] Chen | Proximal Gradient Descent-Ascent: Variable Convergence under K L Geometry[END_REF], and we, therefore, do not explicit all of the technical details. We indicate where changes were needed and develop the modifications in the Appendix for completeness. For the sake of clarity, we denote

ζ u ≜ arg max ζ∈∂ρ(0) E[ζ(1 -γζ)J(u)] which is well defined as ζ → E[ζ(1 -γζ)J(u)
] is a weakly continuous and strongly concave function over a weakly compact set. Proposition 5.1. Let L > 0 be such that J has L-Lipschitz gradient, and

µ > 0 such that ζ → E[ζ(1 -γζ)∇J(u n )] is µ-strongly concave. We define K = L
µ . If the gradient rates η and β verify η < 1 K 3 (L+3) 2 and β < 1 L , then the sequences (u n , ζ n ) given by the update rule [START_REF] Pfeiffer | Risk-averse Merton's Portfolio Problem[END_REF] verify:

∥u n+1 -u n ∥ U → 0, ∥ζ n+1 -ζ n ∥ 2 L 2 F T → 0, ∥ζ n -ζ un ∥ 2 L 2 F T → 0,
Proof. As mentioned above, we generalise the proof from [22, Proposition 2] with f (u, ζ) = E[ζ(1 -γζ)J(u)], g = 0 and h = δ ∂ρ(0) , the convex characteristic function of ∂ρ(0). The key argument is to consider the function

H(u, ζ) = Φ γ (u) + 1 - 1 4K 2 ∥ζ -ζ u ∥ 2 L 2 F T ( 17
)
and show that H(u n , ζ n ) is decreasing over the iterations of the algorithm.

To prove it, we must show the Lipschitzianity of the functions u → ζ u and u → ∇Φ γ (u). The first one is straightforward, as shown in [22, Proposition 1].

The second one however needs some adjustments to our infinite-dimensional setting. Using the version of the Danskin theorem found in [START_REF] Bernhard | On a theorem of Danskin with an application to a theorem of Von Neumann-Sion[END_REF] with the weak topology, we can first prove the following lemma.

Lemma 5.2. Φ γ is Fréchet differentiable and its gradient is given by

∇Φ γ = E[ζ u (1 -γζ u )∇J(u)].
The Lipschitzianity of ∇Φ γ (u) can be deduced as follows:

∥Φ(u 2 ) -Φ(u 1 )∥ = ∥∇ u J γ (u 2 , ζ u 2 ) -∇ u J γ (u 1 , ζ u 1 )∥ ≤ L(∥u 2 -u 1 ∥ + ∥ζ u 2 -ζ u 1 ∥) ≤ L(1 + K)∥u 2 -u 1 ∥
With these two Lipschitzianity results, we can execute the same computations as in [22, Proposition 2] and prove the inequality

H(u n+1 , ζ n+1 ) ≤ H(u n , ζ n )-2∥u n+1 -u n ∥ 2 U - 1 4K 2 (∥ζ n+1 -ζ u n+1 ∥ 2 L 2 F T +∥ζ n -ζ un ∥ 2 L 2 F T
).

From this inequality we can deduce the sequence H(u n , ζ n ) is decreasing and therefore converges to H * . By summing the inequality over the iterations, we obtain a telescoping series that yields

N n=0 2∥u n+1 -u n ∥ 2 U - 1 4K 2 (∥ζ n+1 -ζ u n+1 ∥ 2 L 2 F T + ∥ζ n -ζ un ∥ 2 L 2 F T ) ≤ N n=0 H(u n , ζ n ) -H(u n+1 , ζ n+1 ) ≤ H(u 0 , ζ 0 ) -H * .
This shows that the series on the left is convergent, meaning that its terms tend to zero, which proves our proposition.

This result additionally shows some other interesting properties of the algorithm, which are of straightforward proof.

Corollary 5.3.

1. Since ∥ζ n -ζ un ∥ 2 L 2 F T
tends to zero and H decreases, the sequence Φ γ (u n ) converges to a finite limit. 2. If the cost J is coercive, the sequence u n is bounded in U and is thus in a weakly compact set.

To go further, one would wish to find a generalization of [22, Theorem 1] by showing that all weak limit points of our sequence u n are "gamma-critical points" and verify equation [START_REF] Archibald | A Stochastic Gradient Descent Approach for Stochastic Optimal Control[END_REF]. Unfortunately, we could not prove this in the general case, as a weakly converging sequence of controls u n does not necessarily yield a converging sequence of states x un . Nonetheless, we can demonstrate the result in the slightly more restrictive setting of linear affine controls, as shown in the next subsection.

Convergence of the gradient in Setting B

As mentioned previously, restricting to control affine systems allows us to use important theoretical properties and improve convergence guarantees for our gradient ascent-descent algorithm. The most important property we take advantage of is the compacity of the input-output mapping u → X u , as shown in the following lemma: Lemma 5.4. Let u n be a sequence controls in U, and p ≥ 2. Under the hypothesis of Setting B, we have that

u n U ⇀ u ⇒ X un C p F → X u .

Consequently, if we decompose the cost J of Setting B as

J(u) = C u ∥u∥ 2 U + J X (X u ), we have u n U ⇀ u ⇒ ∇ u (J X (X un )) C 2 F → ∇ u (J X (X u )).
The proof of this lemma requires some long and technical computations that we detailed in Appendix A.3. This compacity result, in turn, allows us to prove the following: Proposition 5.5. Let u be a weak limit of a subsequence of u n , then for any ζ ∈ ∂ρ(0) that is a weak limit of a subsequence of ζ n we have that

⟨E[ζ ∇J(u)], h⟩ U = ∥h∥ U O(γ) ∀v ∈ U. ( 18
)
Proof. Let us first use our result on the convergence of the general algorithm. We can see that

v n+1 -v n = -ηE[(ζ n -γζ 2 n )∇ v J(u n )], = -ηE[ζ n ∇ v J(u n )] + ηγE[ζ 2 n ∇ v J(u n )]. Since ζ n is bounded in L ∞ (Ω, R), ∇ v J(u n ) is bounded in L 2 (Ω, R), and (v n+1 , λ n+1 ) -(v n , λ n ) tends to 0, there exists a constant C > 0 such that ∥E[ζ n ∇ v J(u n )]∥ L 2 (0,T ) ≤ Cγ = O(γ). ( 19
)
Below, we will implicitly overload the constant C. Similarly, we have

|E[ζ n ∇ λ J(u n )]| ≤ Cγ = O(γ). ( 20
)
We now go back to the equality we want to prove. For this, let h be a control in U, and let u and ζ be weak limits of respectively u n and ζ n . We may compute

⟨E[ζ∇J(u)], h⟩ U = ⟨E[(ζ -ζ n )∇J(u)], h⟩ U + ⟨E[ζ n (∇J(u) -∇J(u n ))], h⟩ U + ⟨E[ζ n ∇J(u n )], h⟩ U .
Using equations ( 19) and ( 20), we have that

|⟨E[ζ n ∇J(u n )], h⟩ U | = ∥h∥ U O(γ) ∀v ∈ U.
and since ζ n weakly converges to ζ, with an abuse of notation we may write

⟨E[ζ∇J(u)], h⟩ U = ⟨E[ζ n (∇J(u) -∇J(u n ))], h⟩ U + ∥h∥ U O(γ).
By further decomposing J(u) = C u ∥u∥ 2 U + J X (X u ), we additionally compute

⟨E[ζ∇J(u)], h⟩ U = 2C u ⟨E[ζ n (u -u n )], h⟩ U + ⟨E[ζ n (∇ u J X (X u ) -∇ u J X (X un ))], h⟩ U + ∥h∥ U O(γ).
Using Lemma 5.4 and the fact that u = (v, λ) is deterministic, again with an abuse of notation, we can write

⟨E[ζ∇J(u)], h⟩ U = 2C u E[ζ n ]⟨u -u n , h⟩ U + ∥h∥ U O(γ),
and thanks to the weak convergence of u n to u * we finally obtain that

⟨E[ζ∇J(u)], h⟩ U = ∥h∥ U O(γ).
This result implies that, as we make γ tend to 0, the solutions obtained will get closer to our original problem, as summarized in the following corollary.

Corollary 5.6. Let u γ be a weak limit of a sequence of controls u γ n obtained with the update rule [START_REF] Pfeiffer | Risk-averse Merton's Portfolio Problem[END_REF] with parameter γ. For any decreasing sequence (γ k ) k∈N tending to 0, we can define a sequence (u γ k ) k∈N of controls whose weak limits all verify the original optimality condition [START_REF] Ruszczyński | Risk-Averse Control of Undiscounted Transient Markov Models[END_REF].

This corollary indicates that if we find with our algorithm a γ-critical point for smaller and smaller values of γ, we can find a control law that verifies the necessary optimality condition. Theoretically, lowering γ should imply lowering the gradient step η to satisfy the assumptions in Proposition 5.1, implying a slower convergence in practice. However, the numerical results we present in the next section show that the two convergences, i.e., with γ = 0 or γ ∼ 0, are essentially equivalent, thus there is no downside in taking γ > 0 very small, e.g., in the order of the epsilon machine.

Numerical Simulations

In this section, we present some simulation results. Our algorithm was implemented and tested on two systems, one in Setting A, and the other in Setting B. We respectively denote them as System A and System B. By testing our algorithm on these two different systems, we can evaluate its effectiveness and robustness in solving optimal control problems in both settings. System A is an academic non-linear system that simulates a steering problem, where a vehicle is controlled through its angular velocity while moving at a constant speed. It is subject to perturbations in the state which, for instance, may model measurement errors. On the other hand, System B is a control affine system that models a real-world robot called AstroBee [START_REF] Smith | Astrobee: A New Platform for Free-Flying Robotics on the International Space Station[END_REF], designed to operate in a zero-gravity environment, such as a space station. The diffusion term may, for instance, model sudden gusts of air randomly generated by onboard air purification systems, which destabilizes the robot. For simplicity, we designed and leveraged a very simple stochastic model to model the disturbances; note that other models exist in the literature, e.g., [START_REF] Siurna | Transient solutions to a stochastic model of ventilation[END_REF][START_REF] Thavlov | A non-linear stochastic model for an office building with air infiltration[END_REF].

System presentation

For both systems, the objective is to move from one point to another in a corridor while following a reference trajectory and avoiding obstacles. The reference trajectory is given by a reference control that provides a satisfying solution to the deterministic problem. However, when we add perturbations to the system, the trajectories resulting from the reference control often collide with obstacles. To overcome this challenge, our algorithm searches for a control strategy that minimizes a risk measure of the stochastic cost function. The cost function is composed of two terms: a trajectory-tracking term that tracks the reference trajectory and a term J obstacles (X) that heavily penalizes trajectories that collide with obstacles. This penalization method as presented in [START_REF] Bonnans | Numerical Optimization: Theoretical and Practical Aspects[END_REF] is a classic optimal control approach to get rid of the constraints and is widely used in the literature [START_REF] Gilbert | Distance functions and their application to robot path planning in the presence of obstacles[END_REF][START_REF] Rasekhipour | A Potential Field-Based Model Predictive Path-Planning Controller for Autonomous Road Vehicles[END_REF][START_REF] Bonalli | GuSTO: Guaranteed Sequential Trajectory optimization via Sequential Convex Programming[END_REF].

J(u) = ∥u∥ 2 U + T 0 ∥X u (t) -X ref (t)∥ 2 dt + J obstacles (X u ).
The term J obstacles (X) acts as a potential that is extremely high in the obstacles and thus forces the trajectories away. The risk measure significantly weighs the worst-case scenarios, which means that the cost of a control strategy will increase substantially if some random trajectories produced by the control collide with an obstacle. A control minimizing the risk measure should therefore avoid collisions with high probability.

Results

As shown in Figure 1, our algorithm converges arbitrarily close to solutions for both System A and System B, as long as we use a sufficiently small step size. Additionally, implementing the algorithm with an infinitesimal value of γ results in the same level of convergence as when using a non-negligible value of γ. In particular, although theoretical proofs typically require γ to be strictly greater than 0, our numerical results suggest that it may not always be necessary to achieve convergence.

In Table 1, we evaluate the effectiveness of our algorithm under three different values of α: α = 0.05, α = 0.2, and α = 1 (the latter value corresponds to the minimization of the mean). By testing our algorithm under these different values of α, we can assess the impact of different risk levels on the control strategy. Our main metric is the percentage of trajectories produced by the control strategy that resulted in collisions with obstacles. We compare these percentages between strategies obtained with different values of α to assess its impact on the level of risk in the control strategy. Systems α = 1 α = 0.2 α = 0.05 System A 0.353% 0.112% 0.011% System B 0.049% 0.007% 0.004% To evaluate the robustness of the control strategy, we can increase the size of the obstacles to see how much the probability of collision increases. We rapidly re-trained the controls in this modified problem, with obstacles of a bigger size. This allows us to assess the ability of the control strategy to handle more challenging and uncertain conditions. From Table 4, for both System A and System B, we observe the control solutions that minimized the risk measure with a low value of α are the most conservative. This results in trajectories that keep a larger distance from obstacles but also deviate more from the reference trajectory.

To further assess the impact of the parameter α on the trade-off between mean scores and worst-case scenarios, we computed the probability densities of the costs of the obtained control strategies. By comparing the probability densities of the costs under different values of α, we can identify the optimal level of risk for each system and obtain a control strategy that balances the trade-off between mean scores and worst-case scenarios. As expected, density distributions have higher means but smaller tails for lower values of α. This suggests that a more conservative risk estimate leads to a more robust control strategy with lower worst-case scenarios but also a lower mean score. Obstacles bonus size α = 1 α = 0.2 α = 0.05 +0% 0.353% 0.112% 0.011% +5% 0.367% 0.111% 0.011% +10% 0.470% 0.145% 0.013% +15% 0.600% 0.177% 0.015% +20% 0.735% 0.225% 0.020%

Obstacles bonus size

α = 1 α = 0.2 α = 0.05 +0% 0.049% 0.007% 0.004% +5% 0.044% 0.006% 0.004% +10% 0.048% 0.006% 0.006% +15% 0.054% 0.006% 0.006% +20% 0.066% 0.008% 0.006%

Conclusion and future directions

In this paper, we develop an algorithm to solve risk-averse stochastic optimal control problems subject to non-linear stochastic differential equations, where a risk measure replaces expectation in the cost. By leveraging duality results for coherent risk measures, we recast our minimization problem into a min-max problem that we solve with a gradient descent-ascent approach. We prove convergence properties for this algorithm, which we showcase through appropriate numerical implementation on non-trivial control systems. Our results show our algorithm genuinely converges under very general assumptions and settings. Importantly, our approach of minimizing the cost under a risk measure instead of the expectation yields more robust control strategies, which are capable of considerably better mitigating the uncertainties generated by the diffusion term in the stochastic differential equation.

Several exciting avenues for further exploration and testing are listed hereafter. One direction would consist in extending our algorithm to sys-tems that are modeled through more complex stochastic processes, such as processes with control variables in the diffusion or processes with jumps. Control-dependent diffusion would allow us to model systems where the uncertainties are caused by inputs as well. While putting a control term in the diffusion is already theoretically feasible, we found it to be computationally very expensive as it requires computing conditional expectations, calling for deeper investigation. In particular, recent techniques found in the literature [START_REF] Chessari | Numerical methods for backward stochastic differential equations: A survey[END_REF] may allow us to efficiently compute these conditional expectations. Another direction consists of testing our approach for other risk measures than just the CV@R. Additionally, we are interested in extending our approach to more sophisticated agents, e.g., soft robots. Soft robots have recently received a particular surge of interest as they appear in an increasing number of applications. An emerging model for the dynamics of soft robots hinges upon infinite-dimensional dynamics ruled by Partial Differential Equations (PDEs) [START_REF] Della Santina | Model-Based Control of Soft Robots: A Survey of the State of the Art and Open Challenges[END_REF]. It would be beneficial to extend our approach to handle stochastic PDEs: we look forward to future developments in this area. in the two trajectories X u and X h is given by

X u (s) -X h (s) = s 0 [b(r, X u (r), u(r)) -b(r, X h (r), h(r))] dr + s 0 [σ(r, X u (r)) -σ(r, X h (r))] dW r .
where we denote the control value u(r) = v u (r) + K λu (X u (r)). If we note C > 0 a positive constant (that will be overloaded in the rest of the proof), we can bound the C 2 F norm of the difference above by

E sup 0≤s≤t ∥X u (s) -X h (s)∥ 2 ≤ CE sup 0≤s≤t s 0 [σ(r, X u (r)) -σ(r, X h (r))] dW r 2 + CE sup 0≤s≤t s 0 [b(r, X u (r), u(r)) -b(r, X h (r), h(r))] dr 2 .
Using Cauchy-Schwartz in the deterministic integral and the Burkholder-Davis-Gundy inequality in the stochastic integral, we get:

E sup 0≤s≤t ∥X u (s) -X h (s)∥ 2 ≤ CE t 0 ∥σ(r, X u (r)) -σ(r, X h (r))∥ 2 dr + CE t 0 ∥b(r, X u (r), u(r)) -b(r, X h (r), h(r))∥ 2 dr .
By using the lipschitzianity of the dynamics and taking the sup in the integral, we get

E sup 0≤s≤t ∥X u (s) -X h (s)∥ 2 ≤ CE t 0 sup 0≤r≤s ∥X u (r) -X h (r)∥ 2 ds + CE t 0 ∥u(s) -h(s)∥ 2 ds .
We can write

u(s) -h(s) = v u (s) -v h (s) + K λu (X u (s)) -K λ h (X u (s)) + K λ h (X u (s)) -K λ h (X h (s)),
and therefore have

E sup 0≤s≤t ∥X u (s) -X h (s)∥ 2 ≤ CE t 0 sup 0≤r≤s ∥X u (r) -X h (r)∥ 2 ds + CE t 0 ∥v u (s) -v h (s)∥ 2 ds + CE t 0 ∥K λu (X u (s)) -K λ h (X u (s))∥ 2 ds + CE t 0 ∥K λ h (X u (s)) -K λ h (X h (s))∥ 2 ds .
Using the assumptions on K λ (X) in ( 6), we obtain that

E sup 0≤s≤t ∥X u (s) -X h (s)∥ 2 ≤ CE t 0 sup 0≤r≤s ∥X u (r) -X h (r)∥ 2 ds + C t 0 ∥v u (s) -v h (s)∥ 2 ds + C 1 + E t 0 ∥X u (s)∥ 2 ds ∥λ u -λ h ∥ 2 + C(1 + ∥λ h ∥ 2 )E t 0 ∥X u (s) -X h (s)∥ 2 ds .
The dynamics of the SDE (3) are continuous in u = (v, λ) and the solution X u changes continuously with the dynamics [START_REF] Yong | Stochastic controls: Hamiltonian systems and HJB equations[END_REF], implying that E t 0 ∥X u (s)∥ 2 ds is a continuous function of u. We may therefore introduce C u,h , an overloaded constant that depends continuously on u and h. From this remark, by taking the sup in the integral above, we may infer that

E sup 0≤s≤t ∥X u (s) -X h (s)∥ 2 ≤C u,h E t 0 sup 0≤r≤s ∥X u (r) -X h (r)∥ 2 ds + C u,h t 0 ∥v u (s) -v h (s)∥ 2 ds + ∥λ u -λ h ∥ 2 .
Finally, the Gronwall lemma applied to t → sup

0≤s≤t ∥X u (s) -X h (s)∥ 2 yields E sup 0≤s≤t ∥X u (s) -X h (s)∥ 2 ≤ C u,h ∥u(s) -h(s)∥ 2 U .
In what follows, we denote y u,h the unique element of L 2 F (Ω × [0, T ], R n ) that verifies the linearized SDE Lemma Appendix A.2. Let u = (v u , λ u ) and h = (v h , λ h ) be two controls in U, let us consider the following linear SDE with stochastic dynamics:

         dY (t) = [∂ x b(t, X u (t), u(t)) + ∂ u b(t, X u (t), u(t))K ′ λu (X u (t))]Y (t) +∂ u b(t, X u (t), u(t))[v h (t) + ∇ λ K λu (X u (t))λ h ]dt +∂ x σ(t, X u (t))Y (t)dW t , Y (0) = 0.

This equation has a unique solution denoted y

u,h ∈ L 2 F (Ω × [0, T ], R n )
. This solution can be written explicitly as

y u,h (t) = ϕ u (t) t 0 ϕ -1 u (s)∂ u b(s, X u (s), u(s))[v h (s) + ∇ λ K λ (X u (s))λ h ]ds
where ϕ u , defined in [START_REF] Laurent | Optimality conditions in variational form for non-linear constrained stochastic control problems[END_REF], is the resolvant of the linear SDE [START_REF] Laurent | Optimality conditions in variational form for non-linear constrained stochastic control problems[END_REF]. Consider ϵ a positive constant and denote δX u,ϵ,h ≜ X u+ϵh -X u . Then,

∥δX u,ϵ,h -ϵy u,h ∥ C 2 F ≤ C(ϵ)ϵ∥h∥ L 2 F = o(ϵ)
where C(ϵ) > 0 is a real positive function of ϵ that tends to zero as ϵ → 0.

Proof. The process Y is the solution of a linear equation of the form

dY (t) = (A(t)Y (t) + r(t))dt + C(t)Y (t)dW t ,
where A(t), C(t) and r(t) are F t -adapted processes. If A(t) and C(t) are L ∞ F (Ω×[0, T ]), and r(t) is L 2 F (Ω×[0, T ]), then [START_REF] Yong | Linear forward-backward stochastic differential equations with random coefficients[END_REF] ensures us that there exists a unique strong solution that can be written explicitly. The assumptions made on K λ (x) ensure us that

∥K ′ λ (x)∥ ≤ C(1 + ∥λ∥), ∥∇ λ K λ (x)∥ ≤ C(1 + ∥x∥), which in turn implies ∥A∥ ∞ ≤ ∥∂ x b∥ ∞ + C∥∂ u b∥ ∞ (1 + ∥λ u ∥) ∥C∥ ∞ ≤ ∥∂ x σ∥ ∞ ∥r∥ L 2 ≤ ∥∂ u b∥ ∞ (∥v h ∥ L 2 + C(1 + ∥X u ∥ C 2 F )∥λ h ∥)
The solution therefore exists and can be written explicitly.

As for the approximation, the functions δX u,ϵ,h and y u,h follow the same SDE up to a o(ϵ) term. Therefore, the proof follows the exact same steps as the proof of Lemma Appendix A.1, which we then avoid detailing here.

We can now move to the proof of Lemma 3.1. To prove differentiability, we start by estimating the difference

δJ u,ϵ,h = J(u + ϵh) -J(u) = T 0 [f (t, X u+ϵh (t), v u (t) + ϵv h (t) + K λu+ϵλ h (X u+ϵh (t))) -f (t, X u (t), v u (t) + K λu (X u (t)))]dt + g(T, X u+ϵh (T )) -g(T, X u (T )) = T 0 ∂ x f (t, X u (t), u(t))δX u,ϵ,h (t)dt + T 0 ∂ u f (t, X u (t), u(t))[ϵv h (t) + K λu+ϵλ h (X u+ϵh (t)) -K λu (X u (t))]dt + ∂ x g(T, X u (T ))δX u,ϵ,h (T ) + o(ϵ) + o(δX u,ϵ,h ).
where δX u,ϵ,h ≜ X u+ϵh -X u . Due to Lemma Appendix A.1, we have

∥X u+ϵh -X u ∥ 2 C 2 F ≤ Cϵ 2 ∥h∥ 2 U ,
where C is a positive constant, meaning that o(δX u,ϵ,h ) = o(ϵ). Using this latter fact, and the assumption that (λ, x) → K λ (x) has lipschitz gradient, meaning it verifies [START_REF] Ruszczyński | A Risk-Averse Analog of the Hamilton-Jacobi-Bellman Equation[END_REF], we may use the Taylor expansion of K λ (x) and write

K λu+ϵλ h (X u+ϵh (t))-K λu (X u (t)) = K ′ λu (X u (t))δX u,ϵ,h (t)+∇ λ K λu (X u (t)
)ϵλ h +o(ϵ), and then finally compute

δJ u,ϵ,h = T 0 [∂ x f (t, X u (t), u(t)) + ∂ u f (t, X u (t), u(t))K ′ λ (X u (t))]δX u,ϵ,h (t)dt + T 0 ∂ u f (t, X u (t), u(t))ϵ[v h (t) + ∇ λ K λ (X u (t))λ h ]dt + ∂ x g(T, X u (T ))δX u,ϵ,h (T ) + o(ϵ).
We now need an estimate for δX u,ϵ,h that is linear in ϵh. In order to obtain this latter, we first notice that δX u,ϵ,h solves the following SDE

dδX u,ϵ,h (t) = [∂ x b(t, X u (t), u(t)) + ∂ u b(t, X u (t), u(t))K ′ λ (X u (t))]δX u,ϵ,h (t)dt + ∂ u b(t, X u (t), u(t))ϵ[v h (t) + ∇ λ K λ (X u (t))λ h ]dt + o(ϵ)dt + [∂ x σ(t, X u (t))δX u,ϵ,h (t) + o(ϵ)] dW t .
Lemma Appendix A.2 yields that the solution of the linearized SDE is indeed a good approximation. The gradients of f and g are bounded, we can therefore replace δX u,ϵ,h by the new estimate ϵy u,h and obtain

δJ u,ϵ,h = T 0 [∂ x f (t, X u (t), u(t)) + ∂ u f (t, X u (t), u(t))K ′ λ (X u (t))]ϵy u,h (t)dt + T 0 ∂ u f (t, X u (t), u(t))ϵ[v h (t) + ∇ λ K λ (X u (t))λ h ]dt + ∂ x g(T, X u (T ))ϵy u,h (T ) + o(ϵ).
By injecting the explicit formula of y u,h and inverting the integral signs, we finally obtain

δJ u,ϵ,h = ϵ T 0 I u (s)ϕ -1 u (s)∂ u b(s, X u (s), u(s))[v h (s) + ∇ λ K λ (X u (s))λ h ]ds + o(ϵ), (A.2)
By using the Riesz representation of D u J(u * ) and using equation [START_REF] Ruszczyński | Risk-averse dynamic programming for Markov decision processes[END_REF] we get inf

ζ∈∂ρ(J(u * )) ⟨E [-ζ ∇J(u * )] , u -u * ⟩ U ≤ 0, ∀u ∈ U. (A.3)
By taking the supremum of (A.3) over u ∈ U, we infer that

sup u∈U inf ζ∈∂ρ(J(u * )) ⟨E [-ζ ∇J(u * )] , u -u * ⟩ U ≤ 0. (A.4)
Now, we need to invert the sup with the inf in the above expression, and for this, we use the following version of Sion min-max theorem [START_REF] Simons | Minimax Theorems and Their Proofs[END_REF]:

Theorem Appendix A.4.
Let X and Y be two convex subsets of Haussdorff topological spaces with X being compact, and consider a continuous map f such that

• x → f (x, y) is convex in X for all y ∈ Y , • y → f (x, y) is concave in Y for all x ∈ X.
Then, it holds that

sup y∈Y inf x∈X f (x, y) = inf x∈X sup y∈Y f (x, y).
In our setting, we select Y = U, which is convex, and X = ∂ρ(J(u * )), which is convex, weakly-* compact and therefore weakly compact since it is a subset of a Hilbert space. Also, by definition, the mapping

(ζ, u) → ⟨E [-ζ ∇J(u * )] , u -u * ⟩ U (A.5)
is weakly continuous and concave in u, as well as weakly continuous and convex in ζ. Therefore, Theorem Appendix A.4 can be applied to equation (A.4), yielding inf

ζ∈∂ρ(J(u * )) sup u∈U ⟨E [-ζ ∇J(u * )] , u -u * ⟩ U ≤ 0. (A.6)
Since the supremum of a family of lower semi-continuous functions is also lower semi-continuous, the mapping

ζ → sup u∈U ⟨E [-ζ ∇J(u * )] , u -u * ⟩ U
is lower semi-continuous, and therefore it attains a minimum which we denote by ζ * ∈ ∂ρ(J(u * )). Thanks to this latter remark, from (A.6) we infer that

⟨E [ζ * ∇J(u * )] , u -u * ⟩ U ≥ 0, ∀u ∈ U,
and the conclusion follows from the fact that U is a vector space.

Appendix A.3. Proof of Lemma 5.4 (compacity of the output of a control affine system) We recall we are under the hypothesis of Setting B. Proof of the compacity of u → X u : Let u n be a sequence of controls in U weakly converging to u. Let p ≥ 2, we prove below the more general convergence

∥X un -X u ∥ p C p F → 0
The theorem of existence of strong solutions of SDEs in [START_REF] Yong | Stochastic controls: Hamiltonian systems and HJB equations[END_REF] ensures us that, since our initial condition is deterministic, X u ∈ C p F . Below, we denote by C > 0 a constant that will be implicitly overloaded throughout the proof. By definition, for every n ∈ N we may compute ∥K λn (X un (s)) (i) -K λn (X u (s)) (i) ∥ p ∥b i (s, X u (s))∥ p ds

∥X un -X u ∥ p C p F ≤ CE sup
+ CE t 0 m i=1
∥K λn (X u (s)) (i) -K λ (X u (s)) To conclude, we now need to prove that E sup 0≤s≤t ∥h n (s, ω)∥ p tends to 0. The following lemma gives us the convergence of sup 0≤s≤t ∥h n (s, ω)∥ p for ω fixed in Ω (see [START_REF] Trélat | Contrôle optimal : théorie & applications[END_REF] for proof). If the sequence f n converges simply to an application f , then it converges uniformly.

Lemma

For ω ∈ Ω, h n (s, ω) converges to 0 by weak convergence of u n . Let us recall that the functions b i are bounded by assumption. We can now show that the sequence h n (•, ω) is uniformly Hölder thanks to Cauchy Schwartz and to the fact that v n is uniformly bounded in L 2 (0, T ) :

|h n (t 2 , ω)-h n (t 1 , ω)| = t 2 t 1 m i=1 (v n ) (i) (s) -(v) (i) (s) b i (s, X u (s, ω))ds ≤ C m i=1   t 2 t 1 (|(v n ) (i) (s)| 2 + |(v) (i) (s)| 2 )ds 1 2 t 2 t 1 ∥b i ∥ ∞ ds 1 2   ≤ C|t 2 -t 1 | 1 2
Using Lemma Appendix A.5, we have that sequence of random variables g n (ω) := sup 0≤s≤t ∥h n (s, ω)∥ p converges pointwise to 0. The random variable g n can be uniformly bounded by the deterministic value

|g n (ω)| ≤ T m T 0 m i=1 (M u + ∥(v) (i) (s)∥ 2 )∥b i ∥ 2 ∞ p 2
with M u a uniform bound of the sequence of controls. We therefore have by dominated convergence that E[g n ] tends to 0 , which in turn gives

∥X un -X u ∥ p C p F → 0.

Lemma 3 . 1 .

 31 The cost function J : U → L 2 F T is Fréchet differentiable. Moreover, let ω ∈ Ω and u ∈ U. Thanks to the Riesz representation theorem in the Hilbert space U, we can represent the differential D u J(u, ω) by an element of

Proposition 4 . 1 .

 41 Fix a control u = (v, λ) ∈ U. Under the assumptions of Setting A, the Riesz representation of the gradient of J at u = (v, λ) is

Figure 1 :

 1 Figure 1: Convergence rates obtained for α = 0.05 and different values of γ, for both the steering system (left) and the AstroBee system (right).

  (a) System A: 100 trajectories, α = 0.05, and γ = 0. (b) System B: 100 trajectories, α = 0.05, and γ = 0.

Figure 2 :

 2 Figure 2: 100 trajectories given by a trained control for both systems.

  (a) Density of the cost for System A. (b) Density of the cost for System B.

Figure 3 :

 3 Figure 3: Densities of the cost for different α and γ = 0.

0≤r≤t r 0 (K.r 0 (b 0 K

 000 b 0 (s, X un (s)) -b 0 (s, X u (s)))dsp + CE sup 0≤r≤t r 0 (σ(s, X un (s)) -σ(s, X u (s)))dW s p ) (i) (s) -(v) (i) (s) b i (s, X u (s))ds p λn (X un (s)) (i) -K λ (X u (s)) (i) b i (s, X u (s))ds p n (s)(b i (s, X un (s)) -b i (s, X u (s)))ds p Let us denote h n (t, ω) ≜ t 0 m i=1 (v n ) (i) (s) -(v) (i) (s) b i (s, X u (s, ω))ds.We now appropriately bound the different integral terms. Hölder's inequality yields sup0≤r≤t (s, X un (s)) -b 0 (s, X u (s)))ds p s, X un (s)) -b 0 (s, X u (s))∥ p ds ≤ C t 0 ∥X un (s) -X u (s)∥ p ds ≤ C t 0 sup 0≤s ′ ≤s ∥X un (s ′ ) -X u (s ′ )∥ p ds.Using that the functions b i are bounded and the assumptions on the feedback (6) paired with Hölder's inequality, we have that λn (X un (s)) (i) -K λ (X u (s)) (i) b i (s, X u (s))ds p

K 2 ≤ CE t 0 ∥X 0 E

 200 λn (X un (s)) (i) -K λ (X u (s)) (i) b i (s, X u (s))ds p ∥X un (s ′ ) -X u (s ′ )∥ p ds + C∥λ n -λ∥ p .For the stochastic integral, thanks to Burkholder-Davis-Gundy inequality (1) combined with Hölder we may computeE sup 0≤r≤t r 0 (σ(s, X un (s)) -σ(s, X u (s)))dW s p , X un (s)) -σ(s, X u (s))∥ 2 ds p un (s) -X u (s)∥ p ds ≤ C t 0 E sup 0≤s ′ ≤s ∥X un (s ′ ) -X u (s ′ )∥ p ds.Gathering all the previous bounds yieldsE sup 0≤s≤t ∥X un (s) -X u (s)∥ p ≤ C t sup 0≤s ′ ≤s ∥X un (s ′ ) -X u (s ′ )∥ p ds + C λ n -λ p + CE sup 0≤s≤t ∥h n (s, ω)∥ p ,and an application of Grönwall inequality allows us to infer that E sup 0≤s≤t ∥X un (s) -X u (s)∥ p ≤ C E sup 0≤s≤t ∥h n (s, ω)∥ p + λ n -λ p .

Appendix A. 5 .

 5 Let a, b ∈ R and let E be a normed vector space. For all n ∈ N, letf n : [a, b] → E be uniformly α-Hölder, that is ∃α, K > 0, ∀n ∈ N, ∀x, y ∈ [a, b] ∥f n (x) -f n (y)∥ ≤ K∥x -y∥ α (A.7)

Table 1 :

 1 Percentage of collisions with obstacles for different α and γ = 0. Control strategies trained with higher values of α yield significantly more trajectories colliding with obstacles.

Table 2 :

 2 Percentage of collisions with obstacles in the steering problem.

Table 3 :

 3 Percentage of collisions with obstacles in the AstroBee system.

Table 4 :

 4 Percentage of collisions with obstacles for different α in the two systems.

  (i) ∥ p ∥b i (s, X u (s))∥ p ds ∥X u (s)∥ p )∥λ n -λ∥ p ds .Since λ n converges, it is bounded. Using the fact that X u is in C p F , by taking the sup we finally have that

	≤ CE
	E sup

t 0 (1 + ∥λ n ∥ p )∥X un (s) -X u (s)∥ p ds + CE t 0 (1 +

Note that in this work, we assume that the state variables are completely observable, leaving the investigation of a more general setting as a future research direction.

Appendix A. Technical proofs

Appendix A.1. Proof of Lemma 3.1 In this appendix, we give the proof of Lemma 3.1. We want to show the differentiability of the cost function. We assume to be dealing with the cost and dynamics in Setting A. We prove first that the function J : U → L 2 F T is Fréchet differentiable. We also retrieve a numerically tractable formula for the Riesz representation of the differential, that we denote ∇J(u) = (∇ v J(u), ∇ λ J(u)). For the sake of clarity, we recall that if u = (v, λ) is a control in U, we denote u(t) the value u(t) ≜ v(t) + K λ (X u (t)), where X u uniquely solves (3) with control [START_REF] Chow | A framework for time-consistent, riskaverse model predictive control: Theory and algorithms[END_REF]. We also denote K ′ λ (x) the differential of K λ (x) with respect to x. We first start giving two preliminary lemmas.

Lemma Appendix A.1. Under assumptions of Setting

Proof. Let u = (v u , λ u ) and h = (v h , λ h ) be two controls in U. The difference where I u is the term

The function u = (v, λ) → J is therefore Fréchet differentiable and we can compute the Riesz representation of its gradient, denoted ∇J(u). We recall that

From the expression of δJ u,ϵ,h in (A.2) we can directly identify the gradient as

Let ρ be a finite coherent risk measure, U and ∂ρ(0) are both convex subsets of a Hilbert space. As shown in 3.1, u → J(u) is also Frechet differentiable. Additionally, the risk function is Hadamard differentiable (H-differentiable) as stated in Theorem 2.2, and thus u → ρ(J(u)) is Hdifferentiable. From [45, Chapter 3.1], we obtain the following optimality condition for H-differentiable functions: Proposition Appendix A.3. Let X be a convex subset of a Banach space and f : X → R an H-differentiable function. If f has a minimum on X reached at x * , it is necessary for x * to verify:

Let u * be a solution to Problem (8). Therefore, the previous proposition yields:

-ρ ′ (J(u * ))

In particular, from the formula of the differential of ρ given in Theorem 2.2 and from the chain rule for directional differentials, we obtain that:

Proof of the compacity of u → ∇ u J X (X u ): For more clarity, let us consider a cost where g = 0. Let u be a control in U, we can write the gradient as

with ϕ u satisfying for t ∈ [0, T ] the linear system

Since functions b i , ∂ x b i , ∂ x σ i and K ′ λ are all uniformly bounded, we can prove in a similar way as previously the compacity of u → ϕ u in C p F . We can also prove the compacity of u → ϕ -1 u as it satisfies a similar SDE. We then conclude using the Cauchy-Schwartz inequality twice on the expectation of the four terms and obtaining

The lipschitzianity of ∂ x f ensures the convergence of the first term to 0 and implies that ∥∂ x f (s, x)∥ ≤ C(1 + ∥x∥) which in turn yields

, which is bounded thanks to the convergence of the sequence X un . Finally, we conclude with the lipschitzianity of ∂ x b i which implies that the last term tends to 0. We therefore have


Alain Queffelec
email: alain.queffelec@u-bordeaux.fr

Vincent Delvigne
email: vincent.delvigne@cnrs.fr

Arthur Leck
email: arthur.leck@gmail.com

Morgane De Parthenay

Ludovic Bellot- Gurlet

Céline Leandri-Bressy

Isaac Shearn
email: isaacshearn@gmail.com

Christian Stouvenot

Morgane De

Ludovic Bellot-Gurlet
email: ludovic.bellot-gurlet@sorbonne-universite.fr

Céline Leandri
email: celine.leandri@culture.gouv.fr

Lithic raw materials provenance in the Caribbean islands: flint, jasper, obsidian and so on

Keywords: Lithic raw materials, Silicites, Obsidian, Ceramic Age, Lesser Antilles, Sourcing

Lithic materials in archaeological records have always been the object of special attention by archaeologists, especially from prehistoric contexts. They draw such attention thanks to their often excellent preservation, and the information that the intrinsic characteristics of the rocks can yield about past behaviors including mobility of goods and people through chemical analysis and provenance studies. Although they do not represent the majority of the artifacts recovered from precolonial archaeological sites in the Lesser Antilles, their diversity necessitates varied knowledge, methods, and skills to be studied. Here, we present preliminary results from lithic specialists in characterization of volcanic rocks, diorite, turquoise used for bead production, and silicites.

Several volcanic glass artifacts have been identified for the first time from Saladoid sites from Guadeloupe and Dominica. However, many questions remain about the geological origin and the petrogenesis of these artifacts given their very specific chemical composition. On the other hand, detailed microtextural and micropaleontological studies of the different siliceous raw materials provide first results identifying and mapping their origin and distribution among the Caribbean islands. As a first glance of what will become a long-term project, results presented in this paper demonstrate the capacity of petrological studies carried out with appropriate technical methods to help us improve our knowledge of the Ceramic Age populations' movement and rock distribution networks.

Introduction

The mobility of Amerindians within and between islands during the Ceramic Age has been a particularly intriguing topic for Caribbean archaeologists. In this paper, we aim to contribute to the knowledge of intra-island mobility, inter-islands exchanges, and regional networks within the Caribbean archipelago. Our previous efforts to describe circulation of lapidary raw materials at a regional scale [START_REF] Queffelec | Local production and long-distance procurement of beads and pendants with high mineralogical diversity in an early Saladoid settlement of Guadeloupe (French West Indies)[END_REF][START_REF] Queffelec | Stone ornaments from Guadeloupe and Martinique Early Ceramic period sites (200 BC -AD 400), detailed analysis and comparison with a Late Ceramic period site (AD 750 -1000)[END_REF][START_REF] Queffelec | First identification of sudoite in the Caribbean Ceramic Age lapidary craftsmanship[END_REF] inspired us to pursue a new direction and investigate the provenance and circulation of lithic raw materials more broadly. Our analysis of the current state of knowledge of lithic raw materials in the Caribbean reveals a lack of thorough studies and descriptions, as well as a shortage of specialists in this domain. To address this gap in knowledge, our project aims to gather specialists in a variety of stones used as utilitarian or lapidary raw material and with this in mind, create, describe, and valorize a reference collection (i.e., lithoteca). The aim is to make it available to the community to enable robust identification using objective criteria and ultimately to contribute to an understanding of the landscape exploitation by people who lived in the Caribbean.

The current state of knowledge regarding inter-and intra-island mobility based on lithic material distribution in the Caribbean is heavily reliant on the observation of few key stone sources. These include i) the identification of Antigua's flint in numerous archaeological sites (e.g., [START_REF] Knippenberg | The provenance of flint in the Leeward region, West Indies[END_REF][START_REF] Knippenberg | Chert Sourcing in the Northern Lesser Antilles : The Use of Geochemical Techniques in Discriminating Chert Materials. Crossing the Borders: New Methods and Techniques in the Study of Archaeological Materials from the Caribbean[END_REF][START_REF] Van Gijn | Flint exploitation on Long Island, Antigua, West-Indies[END_REF], ii) the studies of the distribution of the distinctive rocks from St. Martin found in several sites on that island and on others [START_REF] Knippenberg | Distribution and exchange of lithic materials : Three-pointers and axes from St. Martin[END_REF][START_REF] Knippenberg | Stone artefact production and exchange among the Lesser Antilles[END_REF][START_REF] Stouvenot | Examen pétrographique sommaire des objets lithiques de la fouille de la Place Saint-François à Basse-Terre (Guadeloupe)[END_REF], and iii) the identification of radiolarite from La Désirade in several archaeological sites [START_REF] De Waal | On the reconstruction of prehistoric social territories : The La Désirade lithic workshops and the distribution of La Désirade chert (French West Indies)[END_REF].

To achieve our goals, we have organized our project into three distinct work packages, each focused on the different highly specialized knowledge and methods required for the analysis of i) silicites 1 , ii) volcanic rocks, and iii) greenstones. Our approach involves the combination of pre-existing lithotecas with new fieldwork to establish a comprehensive georeferenced collection. We undertake thorough characterization of the samples, using binocular observations for silicites and conducting geochemical and mineralogical analysis for volcanic rocks, in addition to documenting the criteria of identification with photography. Finally, all data will be used for the study of complete archaeological collections without subsampling, allowing the identification of both common and rare materials. By producing new knowledge on these rocks and creating a reference collection, we hope to contribute to the wider archaeological community's understanding of the Caribbean's pre-Columbian past.

Method

The petrological characterization of geological samples of silicites follows the method developed by P. Fernandes and J.-P. Raynal [START_REF] Fernandes | Pétroarchéologie du silex : Un retour aux sources[END_REF]. This method involves filling out three observation and description grids, which were developed based on the concept of the evolutionary chain of silicites [START_REF] Delvigne | Géoressources et expressions technoculturelles dans le sud du Massif central : Des déterminismes et des choix[END_REF][START_REF] Delvigne | Taphonomic method and database to establish the origin of silicites from the Middle-recent Gravettian of La Picardie open-air site (Indre-et-Loire, France)[END_REF][START_REF] Delvigne | Quand la pétroarchéologie questionne la notion de site : États de surface et taphonomie des objets lithiques de La Belle-Roche[END_REF][START_REF] Fernandes | Itinéraires et transformations du silex : Une pétroarchéologie refondée, application au Paléolithique moyen[END_REF][START_REF] Thiry | Driving forces for the weathering and alteration of silica in the regolith : Implications for studies of prehistoric flint tools[END_REF]. This concept was formalized following the observation of the evolution and transformation of this type of material during their non-anthropogenic path, particularly after their genesis. Indeed, after being brought to the surface, they undergo mechanical and/or chemical alterations that are characteristic of the different environments they crossed. The goal of petrological characterization is therefore to identify the different processes at the origin of such alteration. It should be noted that petrology (i.e., the study of rock evolution) is coupled with petrographic analysis, which aims to describe the composition and structure of silicites, allowing for the distinction of genetic types. All observations carried out using exoscopy at low magnification (50x to 200x), with the sample under a slice of water, were recorded in a data table, which was then synthesized into a summary table. This table lists the featured elements (authigenic, non-organic, and organic), as well as inherited structures from the protolith (thus pre-genetic), during siliceous diagenesis, and during different moments of alteration (thus post-genetic). This very detailed description is necessary for identification of archaeological artifacts and possibility of replication by other archaeologists. Applied to the archaeological record, this method also makes it possible to differentiate the alterations undergone before the collection by the prehistoric collectives (predepositional), from those undergone during the burial phase in the archaeological site (postdepositional). In this way, it is possible to determine the location of the outcrop collected by the prehistoric groups, also specifying the type of deposit exploited (primary or secondary).

Non-destructive elemental analysis by PIXE (Particle Induced X-Ray Emission) was performed on the AGLAE facility [START_REF] Calligaro | Ion beam analysis with external beams : Recent set-up improvements[END_REF] at the Centre de Recherche et de Restauration des Musées de France (Paris) using the procedures presented in [START_REF] Bellot-Gurlet | A PIXE/fission-track dating approach to sourcing studies of obsidian artefacts in Colombia and Ecuador[END_REF]. This analysis uses a proton beam to stimulate the X-ray fluorescence of most major elements and several trace elements in the sample. Raman spectroscopy was performed on minerals visible on some artifacts' surfaces with a Bruker SENTERRA Raman microspectrometer equipped with a 532 nm laser, and a 1200-line diffraction grating resulting in a resolution of 3 to 5 cm-1.

Portable XRF analyses were performed with a Bruker Tracer 5i, equipped with a SDD detector and a 4W rhodium X-ray tube. Measurements are carried out in two modes of tube operation: a mode at 15 kV and 100 μA to analyze light elements (aluminum to manganese) and a mode at 45 kV and 64 μA to analyze heavier elements (iron to zirconium). Analyses were carried out on the standards and on the artifacts using a collimator to reduce the beam diameter to 3 mm.

Results

Reference collection and preliminary archaeological study for silicites. The analysis of the reference collection created by one of us (CS) and the samples we collected during our field mission allowed for the identification of 24 types of raw materials. They can be grouped in several categories: orthoquartzite (Figure 1a-c), hypersiliceous volcanic rocks (Figure 1d-f), silicified cinerite, hydrothermal siliceous rocks (Figure 1g-i), radiolarite, agate and silicified coral. These types have been described thoroughly and photographed, so that future research can be based on these documents, either by our team or other researchers. Alongside our efforts to complete and describe the reference samples, we conducted a preliminary study of the archaeological site of Stade José Bade located in Marie-Galante [START_REF] Serrand | Stade José Bade : Les occupations cedrosan-saladoïde et troumassantroumassoïde de Tourlourous[END_REF]. Due to time constraints during our first mission, this first study is a first approach to a re-examination of the collection. This collection had already been studied, without extensive description and binocular observations. This initial work identified eight different types of siliceous raw materials. Upon the application of the methodology by trained specialists, and with time dedicated to binocular observations, the same collection of artifacts was robustly described as composed by raw materials belonging to four types, one of which consists of two sub-types. These four types are notably different from the eight categories previously identified, demonstrating the importance of using such a methodology (Figure 2). The detailed description of the artifacts from this site will also allow later identification when we update the reference collection with new raw materials. In the end, the largely dominant Type A was attributed to Antigua's flint, more precisely to the "distal part of the lagoon" (Type A1) and the "entrance of lagoon" (TypeA2) facies of Flinty Bay. While two other types of flint were identified, their sources are still unknown and will hopefully be recognized later in the project. Additionally, we described one type of lithic material that is likely of volcanogenic origin. This preliminary study demonstrates the significant contributions that trained specialists can make to the field.

Several islands in the Caribbean (Figure 3), if not most of them, display chert, flint, jaspers or other siliceous rocks described in the archaeological or geological literature [START_REF] Bérard | Les premières occupations agricoles de l'Arc Antillais, Migrations et insularité[END_REF][START_REF] Donovan | The geology of Barbados : A field guide[END_REF][START_REF] Knippenberg | Stone artefact production and exchange among the Lesser Antilles[END_REF][START_REF] Pantel | Precolumbian flaked stone assemblages in the West Indies[END_REF][START_REF] Ray | Procurement and use of chert from localized sources in Trinidad[END_REF][START_REF] Stouvenot | Matières premières[END_REF]. Our project aims to study these materials as well as archaeological collections. Another future goal for this long-term project is to expand into the Greater Antilles, a region with a complex geological history with great potential to yield a variety of unique materials. Volcanic glass identification in the Caribbean archaeological record Despite the widely accepted belief that obsidian does not occur naturally in the Caribbean nor is found in the archaeological record, we present here evidence of the contrary. We identified twenty artifacts in the literature, recovered from seven archaeological sites spanning Early to Late Ceramic age contexts and located on four islands (Table 1). We were able to analyze sixteen of them (Figure 4). In the Lesser Antilles, only one potential sources of glassy material (i.e., obsidian), is known: a limited outcrop integrated in the "Volcan du Tuf" [START_REF] Boudon | Notice explicative de la carte géologique à 1/20 000e du Massif Volcanique de la Soufrière (Département de la Guadeloupe, Petites Antilles)[END_REF][START_REF] Gadalia | La chaîne volcanique de Bouillante, Basse-Terre de Guadeloupe, (Petites Antilles) : Identitié pétrographique, volcanologique et géodynamique[END_REF][START_REF] Gadalia | Prospection géothermique de la région de Bouillante-Vieux Habitants Guadeloupe[END_REF] In 2007, a preliminary non-invasive PIXE analysis was conducted on 16 artifacts and one sample from the "Volcan du Tuf" [START_REF] Bressy | Caractérisation et recherche de provenance de matériaux lithiques de sites préhistoriques de Guadeloupe « Analyses Obsidiennes Guadeloupe 2007[END_REF]. Subsequently, additional analytical techniques including photographs, microphotographs, pXRF analysis, and Raman spectroscopy were employed in 2019 and 2020 to investigate further most of these artifacts and some discovered after the initial study. The results indicated that nearly all the samples had a comparable composition, which was not consistent with the typical rhyolitic composition of obsidian. Rather, their low silica content classified them as andesitic volcanic glass (Figure 5), Antilles origin given that most of the recent volcanic arc consists of andesitic rocks. However, our global scale literature review did not yield a single reference to natural andesitic glass suitable for knapping, either in the geological or archaeological domains, except for rocks originating from asteroid impacts [START_REF] Rochette | Impact glasses from Belize represent tektites from the Pleistocene Pantasma impact crater in Nicaragua[END_REF]. Raman spectroscopy analysis of minerals inclusions identified some plagioclase feldspars and augite (Figure 6), both standard minerals of andesites but that can be also present in impactites [START_REF] French | Traces of catastrophe : A handbook of shock-metamorphic effects in terrestrial meteorite impact structures[END_REF]. Additional analyses are currently underway to further investigate the cosmic hypothesis. Although the provenance of the glassy material in question remains unknown, we conducted a comprehensive survey of its geographic distribution throughout the Caribbean islands. This inquiry led us to investigate relevant literature, notwithstanding the widely held notion that obsidian is not present in this region. The samples listed in the literature from Maruca and Maisabel, in Puerto Rico, are probably not obsidian (Pers. Comm Rodriguez-Ramos 2021, Siegel 2021). One sample is listed in R. Murphy's dissertation [START_REF] Murphy | The prehistory of Antigua, ceramic age : Subsistence, settlement, culture and adaptation within an insular environment[END_REF], and we aspire to examine it in the future. Some samples are listed in the site of Fond Brûlé in Martinique [START_REF] Mattioni | Salvage Excavations at the Fond-Brule Site, Martinique[END_REF], but these were mistakenly identified 2 . Despite the scrutiny of existing literature, discussions with numerous Caribbean archaeologists, and exchanges held during the IACA congress in Cuba, the only volcanic glass artifact that we could identify outside French islands was the flake from Delices (Del-3) in Dominica [START_REF] Shearn | Pre-Columbian Regional Community Integration in Dominica, West Indies[END_REF].

The other magmatic rock that we aim to study for provenance studies is the diorite. This distinctive black and white rock was widely used for lapidary production during Saladoid and Post-Saladoid periods, all over the archipelago (Figure 7). Few sources are mentioned in the literature, mainly at Tobago [START_REF] Mones | Amerindian lithics of Tobago : A preliminary report[END_REF] and Puerto Rico [START_REF] Turner | Investigation of in situ weathering of quartz diorite bedrock in the Rio Icacos basin, Luquillo Experimental Forest, Puerto Rico[END_REF], but other samples from La Désirade, Saint-Martin, British Virgin Islands, and French Guyana will be also integrated as references. Despite its granular texture, which induces a certain analytical challenges, petrographic observations coupled with geochemical analyses using X-ray fluorescence and X-ray spectroscopy in a SEM have been successfully employed in recent studies (e.g., [START_REF] Burg | Advances in geochemical sourcing of granite ground stone : Ancient Maya artifacts from the middle Belize Valley[END_REF][START_REF] Turbay | Quarry identification and characterization of 2nd century AD Roman granite columns from Ecija (Spain)[END_REF]. Sourcing of turquoise by non-destructive isotopic analysis Turquoise, another rock used for lapidary production during the Early and Middle Ceramic periods, is also studied in this project. This gemstone is well distributed among the Lesser Antilles, even if always as small quantity in the sites (Figure 8). The source of this material, on the other hand, is unknown. The well-known sources of turquoise are located far in the American Cordillera, both in North and South America [START_REF] Queffelec | Provenancing turquoise-A difficult task[END_REF], while a source is mentioned in coastal Brazil [START_REF] Costa | A turquesa de Itacupim, Pará[END_REF], and a single source is known so far in the Caribbean, in the US Virgin Islands [START_REF] Alminas | Geochemistry, mineralogy, and geochronology of the US Virgin Islands[END_REF][START_REF] Tucker | Geochemical Evidence for Metallization on St[END_REF]. The application of non-destructive isotopic analysis of copper, as proposed by [START_REF] Hull | A new approach to determining the geological provenance of turquoise artifacts using hydrogen and copper stable isotopes[END_REF], and of lead, as measured destructively by [START_REF] Thibodeau | Lead isotope analysis as a new method for identifying material culture belonging to the Vázquez de Coronado expedition[END_REF], is not straightforward to implement without risk for the analyzed artifacts and with good reproducibility. We are presently in the process of creating internal standards for isotopic analysis by laser ablation mass spectrometry, to avoid the rather invasive sample preparation of SIMS analysis (Secondary Ion Mass Spectrometry).

Discussion and conclusion

The creation of a reference collection of silicites, and its characterization through microscopic observation, started recently in Guadeloupe and will continue in Martinique, hopefully gathering raw materials from other islands as well. Microscopic observation will be used to characterize the archaeological collections too, building upon a methodology established in continental France over the past decade that has successfully identified and described dozens of different types of lithic raw materials. This approach will prove invaluable in describing the range of raw materials encountered in archaeological contexts. A very preliminary study of an archaeological collection from the Saladoid period (Stade José Bade) already demonstrated the applicability and utility of this method. Future applications to archaeological collections from different islands and periods will rely on these detailed descriptions, enabling the identification of raw materials beyond the sole source of flint on Antigua.

This study also sheds light on the presence of volcanic glass artifacts in the region, which was previously unknown. The study reveals two significant observations: i) Caribbean people from the Saladoid period, living in present-day Basse-Terre, used a local volcanic glass, found today in the "Volcan du Tuf", and ii) Caribbean people from Basse-Terre, Marie-Galante, and Dominica, during this same period, used a unique and hitherto unknown andesitic glass to produce small flakes through bipolar debitage. Two of these flakes have also been found in Troumassoid context on the Ilets de Petite Terre, in an unfortunately unpublished excavation.

The creation of a reference collection of diorites -a rock used for lapidary production on many islands spanning all Ceramic age -will be also pursued in the coming years. Development of the analytical methodology for sourcing such coarse-grained rocks will be needed.

The project is also working towards sourcing turquoise, with a focus on developing a nondestructive method for analyzing the isotopic composition of this mineral. The goal is to create a methodology that can be used on archaeological artifacts such as small beads and pendants from the Early and Middle Ceramic periods in the Lesser Antilles.

The creation of lithic raw material reference collection for all the Caribbean archipelago is a long-term project, and so will be its use for accurate provenance studies of lithic artifacts. The development of such a database for silicites, and its comprehensive and detailed description, will also be helpful for the community of Caribbean archaeologists, who can hopefully apply this method to their own studies. By combining the different lithic raw materials, we aim at determining the litho-space of Caribbean communities during the Ceramic age. This will guide a more general reflection on the circulation of materials in the Caribbean islands.

Figure 1 .

 1 Figure 1. Examples of microphotographs from several types of silicites from Guadeloupe, showing internal texture, presence of featured elements like manages dendrites (dend) or exoclasts (exo).

Figure 2 .

 2 Figure 2. Comparison of silicite identification on the same archaeological collection by two separate specialists with different methodologies and experience.

Figure 3 .

 3 Figure 3. Natural occurrence of silicites mentioned in the literature. Only one point has been located on each island of the Greater Antilles even if more occurrences exist on each of these islands.

Figure 4 .

 4 Figure 4. Macrophotographs of the sixteen glassy artifacts recovered from six archaeological sites.

  and their highly similar composition for minor and trace element suggests that the samples likely share a common origin. The two most recent discoveries -Stade José Bade and Delices (Del-3) -showed the same "non-rhyolitic" chemical composition based on pXRF analysis. The unique artifact made of obsidian sensu stricto (RRC-01) originated directly from the "Volcan du Tuf", as underlined by their common chemical composition. It was excavated at the Ruelle Roches Caraïbes site, only a few kilometers away from the source.

Figure 5 .

 5 Figure 5. Total Alcali Silica (TAS) diagram (mod. after Le Maitre et al., 2005) from the PIXE analysis, showing that most samples are andesitic glasses, while one of the artifacts (RRC-01) and the geological sample of glass from "Volcan du Tuf" are of rhyolitic composition.

Figure 6 .

 6 Figure 6. Mineral identification by Raman spectroscopy of the crystals included in the Del-3 sample (top two spectra). The minerals are identified as labradorite (a plagioclase feldspar) and augite (a clinopyroxene) by comparison with the RRUFF database (bottom two spectra).

Figure 7 .

 7 Figure 7. Map of the archaeological distribution of diorite beads and pendants during the Ceramic age in the Caribbean islands (mod. after Queffelec, 2022).

Figure 8 .

 8 Figure 8. Distribution of the turquoise artifacts in the Caribbean archipelago (mod. after Queffelec, 2022).

Table 1

 1 . It is located a few kilometers from the archaeological site of Ruelle Roches Caraïbes.

	Site	Island	Longitude	Latitude	Nb artifacts Period	Reference
	Bisdary	Basse-Terre	-61.709	15.988	5*	Early Ceramic	(Romon et al., 2006)
	Ruelle Roches Caraïbes Basse-Terre	-61.736	16.038	2	Early Ceramic	(Bonnissent & Chancerel, 2001)
	Tourlourous	Marie-Galante	-61.224	15.889	7*	Early Ceramic	(Colas, 2002)
	Stade José Bade	Marie-Galante	-61.224	15.889	1	Early Ceramic	(Serrand, 2016)
	Place St. François	Basse-Terre	-61.73	15.995	1*	Early Ceramic	(Sellier-Ségard, 2001)
	Caille à Bélasse	Ilets de Petite Terre	-61.11	16.17	2	Late Ceramic	(Gagnepain, 2007)
	Del-3	Dominica	-61.256	15.275	1	Early Ceramic	(Shearn, 2014)
	Mill Reef+	Antigua	-61.679	17.057	1	Early Ceramic	(Murphy, 1999)

. List of archaeological sites with volcanic glass artifacts (+: yet to be confirmed) and their location. *: one of the samples analyzed in 2007 was missing in 2021.

The term silicites designates all the rocks having undergone a diagenetic silicification of chemical or biochemical origin (e.g., flint, chert, silcrete, hydrothermal silica). It allows to avoid the term of "flint", the meaning of which is subject to controversy and misunderstandings depending on whether it is used by geologists, petrographers, sedimentologists or archaeologists. We therefore use the term flint only to designate rocks having undergone an initial diagenesis followed by a siliceous epigenesis and bearing a cortex.

After careful examination of the complete collection from this site, we can conclude that these flakes and nucleus are not made of volcanic glass.

Acknowledgments

Our team want to thank Sebastiaan Knippenberg for his pioneering work and all the information and samples he accepted to share with us. This work was funded by the French Ministry of Culture (DAC Guadeloupe), and the Guadeloupe Régional Council. This study received financial support from the French government in the framework of the University of Bordeaux's IdEx "Investments for the Future" program / GPR "Human Past". We thank the TraceoLab of Lièges University for providing us the equipment for microphotography of the silicites.