Achref Bachouch

Côme Huré

Nicolas Langrené

Huyen Pham

Huyên Pham

Deep neural networks algorithms for stochastic control problems on finite horizon: numerical applications *

Keywords: Deep learning, policy learning, performance iteration, value iteration, Monte Carlo, quantization

published or not. The documents may come

Deep neural networks algorithms for stochastic control problems on finite horizon: numerical applications

Introduction

This paper is devoted to the numerical resolution of discrete-time stochastic control problem over a finite horizon. The dynamics of the controlled state process X = (X n ) n valued in R d is given by X n+1 = F (X n , α n , ε n+1 ), n = 0, . . . , N -1,

X 0 = x 0 ∈ R d , (1.1)
where (ε n ) n is a sequence of i.i.d. random variables valued in some Borel space (E, B(E)), and defined on some probability space (Ω, F, P) equipped with the filtration F = (F n ) n generated by the noise (ε n ) n (F 0 is the trivial σ-algebra), the control α = (α n ) n is an F-adapted process valued in A ⊂ R q , and F is a measurable function from R d × R q × E into R d which is known by the agent. Given a running cost function f defined on R d × R q and a terminal cost function g defined on R d , the cost functional associated with a control process α is

J(α) = E N -1 n=0 f (X n , α n ) + g(X N ) . (1.2)
In this framework, we assume f and g to be known by the agent. The set A of admissible controls is the set of control processes α satisfying some integrability conditions ensuring that the cost functional J(α) is well-defined and finite. The control problem, also called Markov decision process (MDP), is formulated as

V 0 (x 0 ) := inf α∈A J(α), (1.3)
and the goal is to find an optimal control α * ∈ A, i.e., attaining the optimal value: V 0 (x 0 ) = J(α * ). Notice that problem (1.1)-(1.3) may also be viewed as the time discretization of a continuous time stochastic control problem, in which case, F is typically the Euler scheme for a controlled diffusion process. It is well-known that the global dynamic optimization problem (1.3) can be reduced to local optimization problems via the dynamic programming (DP) approach, which allows to determine the value function in a backward recursion by

V N (x) = g(x), x ∈ R d , V n (x) = inf a∈A Q n (x, a), (1.4) with Q n (x, a) = f (x, a) + E V n+1 (X n+1 ) X n = x, α n = a , (x, a) ∈ R d × A.
Moreover, when the infimum is attained in the DP formula (1.4) at any time n by a * n (x) ∈ arg min a∈A Q n (x, a), we get an optimal control in feedback form (policy) given by: α * = (a * n (X * n )) n where X * is the Markov process defined by

X * n+1 = F (X * n , a * n (X * n ), ε n+1 ), n = 0, . . . , N -1, X * 0 = x 0 .
The practical implementation of the DP formula may suffer from the curse of dimensionality and large complexity when the state space dimension d and the control space dimension are high. In [START_REF] Huré | Deep neural networks algorithms for Stochastic Control Problems on finite horizon, part I: convergence analysis[END_REF], we proposed algorithms relying on deep neural networks for approximating/learning the optimal policy and then eventually the value function by performance/policy iteration or hybrid iteration with Monte Carlo regressions now or later. This research led to three algorithms, namely algorithms NNcontPI, Hybrid-Now and Hybrid-LaterQ that are recalled in Section 2, and which can be seen as a natural extension of actor-critic methods, developed in the reinforcement learning community for stationary stochastic problem ( [START_REF] Sutton | Reinforcement Learning[END_REF]), to finite-horizon control problems. Note that for stationary control problem, it is usual to use techniques such as temporal difference learning, which relies on the fact that the value function and the optimal control do not depend on time, to improve the learning of the latter. Such techniques do not apply to finite horizon control problems. In Section 3, we perform some numerical and comparative tests to illustrate the efficiency of our different algorithms, on 100-dimensional nonlinear PDEs examples as in [START_REF] Weinan | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF] and quadratic Backward Stochastic Differential equations as in [CR16], as well as on high-dimensional linear quadratic stochastic control problems. We present numerical results for an option hedging problem in finance, and energy storage problems arising in the valuation of gas storage and in microgrid management. Numerical results and comparisons to quantization-type algorithms Qknn, introduced in this paper as an efficient algorithm to numerically solve low-dimensional control problems, are also provided. Finally, we conclude in Section 4 with some comments about possible extensions and improvements of our algorithms.

Algorithms

We introduce in this section four neural network-based algorithms for solving the discretetime stochastic control problem (1.1)-(1.3). The convergence of these algorithms have been analyzed in detail in our companion paper [START_REF] Huré | Deep neural networks algorithms for Stochastic Control Problems on finite horizon, part I: convergence analysis[END_REF], and for self-contained purpose, we recall in this section the description of these algorithms and the convergence results. We also introduce at the end of this section a quantization and k-nearest-neighbor-based algorithm (Qknn) that will be used as benchmark when testing our algorithms on low-dimensional control problems.

We are given a class of deep neural networks (DNN) for the control policy represented by the parametric functions x ∈ R d → A(x; β) ∈ A, with parameters β ∈ R q , and a class of DNN for the value function represented by the parametric functions: x ∈ R d → Φ(x; θ) ∈ R, with parameters θ ∈ R p . Recall that these DNN functions A and Φ are compositions of linear combinations and nonlinear activation functions, see [START_REF] Goodfellow | Deep learning[END_REF].

Additionally, we shall be given a sequence of probability measures on the state space R d , that we call training measure and denoted (µ n ) N -1 n=0 , which should be seen as dataset providers to learn the optimal strategies and the value functions at time n = 0, . . . , N -1.

Remark 2.1 (Training sets design)

The choice of the training sets is critical for numerical efficiency. This problem has been largely investigated in the reinforcement learning community, notably with multi-armed bandits algorithms [START_REF] Auer | Finite-time Analysis of the Multiarmed Bandit Problem[END_REF], and more recently in the numerical probability literature, see [START_REF] Ludkovski | Simulation methods for stochastic storage problems: a statistical learning perspective[END_REF], but remains a challenging issue. Here, two cases are considered for the choice of the training measure µ n used to generate the training sets on which the estimates at time n will be computed. The first one is a knowledge-based selection, relevant when the controller knows with a certain degree of confidence where the process has to be driven in order to optimize her cost functional. The second case is when the controller has no idea where or how to drive the process to optimize the cost functional.

(1) Exploitation only strategy In the knowledge-based setting, there is no need for exhaustive and expensive (in time mainly) exploration of the state space, and the controller can take a training measure µ n that assigns more points in the region of the state space that is likely to be visited by the optimally-driven process.

In practice, at time n, assuming we know that the optimal process is likely to lie in a region D, we choose a training measure in which the density assigns a lot of weight to the points of D, for example U(D), the uniform distribution in D.

(2) Explore first, exploit later When the controller has no idea where or how to drive the process to optimize the cost functional, we suggest to build the training measures as empirical measures of the process, driven by estimates of the optimal control computed using alternative methods.

(i) Explore first: Use an alternative method to obtain good estimates of the optimal strategy.

In high-dimension: one can for example think of approximating the control at all time by neural network, and obtain a good estimate of the optimal control by performing a global optimization of the function:

J(θ 0 , . . . , θ N +1 ) := E N -1 n=0 f (X n , A(X n ; θ n )) + g(X N ) ,
where X is the process controlled by the feedback control A(.; θ n ) at time n.

(ii) Exploit later: Take the training measures µ n := P Xn , for n = 0, . . . , N -1, where X is driven using the optimal control estimated in step (i); and apply the procedure (1). Such an idea has been recently exploited in [START_REF] Kou | A general Monte Carlo algorithm with monotonicity for stochastic control problems[END_REF].

Remark 2.2 (Choice of Neural Networks) Unless otherwise specified, we use feedforward Neural Networks with two or three hidden layers and d+10 neurons per hidden layer, since we noticed empirically that these parameters were enough to approximate the relatively smooth objective functions considered here. We tried sigmoid, tanh, ReLU and ELU activation functions and noticed that ELU is most often the one providing the best results in our applications. We normalize the input data of each neural network in order to speed up the training of the latter. 2

Remark 2.3 (Neural Networks Training) We use the Adam optimizer, as implemented in TensorFlow, with initial learning-rate set to 0.001 or 0.005, which are the default values in TensorFlow, to train by gradient-descent the optimal strategy and the value function defined in the algorithms described later. TensorFlow takes care of the Adam gradient-descent procedure by automatic differentiation when the function to optimize is an expectation of TensorFlow functions, such as the usual differentiable activation functions sin, log, exp but also popular non-differentiable activation functions such as ReLu:

x → max(0, x).

In order to force the weights and biases of the neurons to stay small, we use an L 2 regularization with parameter mainly set to 0.01, but the value can change in order to make sure that the regularization term is neither too strong or too weak when added to the loss when training neural networks.

We consider a large enough number of mini-batches of size 64 or 128 for the training, depending essentially empirically on the dimension of the problem. We use at least 10 epochs a and stop the training when the loss computed on a validation set of size 100 stops decreasing. We noticed that taking more than one epoch really improves the quality of the estimates.

2

Remark 2.4 (Constraints) The proposed algorithms can deal with state and control constraints at any time, which is useful in several applications:

(X α n , α n ) ∈ S a.s., n ∈ N,
where S is some given subset of R d × R q . In this case, in order to ensure that the set of admissible controls is not empty, we assume that the sets A(x) := a ∈ R q : (F (x, a, ε 1 ), a) ∈ S a.s.

are non empty for all x ∈ S, and the DP formula now reads

V n (x) = inf a∈A(x) f (x, a) + P a V n+1 (x) , x ∈ S.
From a computational point of view, it may be more convenient to work with unconstrained state/control variables, hence by relaxing the state/control constraint and introducing into the running cost a penalty function L(x, a): f (x, a) ← f (x, a) + L(x, a), and g(x) ← g(x) + L(x, a). For example, if the constraint set S is in the form:

S = {(x, a) ∈ R d × R q : h k (x, a) = 0, k = 1, . . . , p, h k (x, a) ≥ 0, k = p + 1, .
. . , q}, for some functions h k , then one can take as penalty functions:

L(x, a) = p k=1 µ k |h k (x, a)| 2 + q k=p+1 µ k max(0, -h k (x, a)).
where µ k > 0 are penalization coefficients (large in practice). 2

Control Learning by Performance Iteration

We present in this section Algorithm 1, which combines an optimal policy estimation by neural networks and the dynamic programming principle. We rely on the performance iteration procedure, i.e. paths are always recomputed up to the terminal time N .

a We denote by epoch one pass of the full training set.

Algorithm NNContPI

Our first algorithm, referred to as NNContPI, is well-designed for control problems with continuous control space such as R q or a ball in R q . The main idea is:

1. Represent the controls at time n = 0, . . . , N -1 by neural networks in which the activation function for the output layers takes values in the control space. For example, one can take the identity function as activation function for the output layer if the control space is R q ; or the sigmod function if the control space is [0, 1].

2. Learn sequentially in time, and in a backward way, the optimal parameters βn for the representation of the optimal control. In particular, notice that the learning of the optimal control at time n highly relies on the accuracy of the estimates of the optimal controls at time k = n + 1, . . . , N -1, computed previously.

Algorithm 1: NNContPI

Input: the training distributions (µ n ) N -1 n=0 ; Output: estimates of the optimal strategy (â n ) N -1 n=0 ; for n = N -1, . . . , 0 do Compute βn ∈ argmin β∈R q E f X n , A(X n ; β) + N -1 k=n+1 f X β k , âk X β k + g X β N (2.1)
where X n ∼ µ n and where X β k N k=n+1 is defined by induction as:

X β n+1 = F X n , A X n ; β , ε n+1 X β k+1 = F X β k , âk X β k , ε k+1 , for k = n + 1, . . . , N -1.
Set ân = A(.; βn ). ân is the estimate of the optimal policy at time n

Algorithm ClassifPI

In the special case where the control space A is finite, i.e., Card(A) = L < ∞ with A = {a 1 , . . . , a L }, a classification method can be used: consider a DNN that takes state x as input and returns a probability vector p(x; β) = (p (x; β)) L =1 with parameters β. Such a usual DNN can be build using k hidden layers with ReLu activation functions, an output layer with L neurons, and a Softmax b activation function for the output layer. Algorithm 2, presented below, is based on this idea, and is called ClassifPI.

b The Softmax function is defined as follows: x → e β 1 x L k=1 e β k x , . . . ,

e β 1 x L k=1 e β k x
where β1, . . . , βL are part of the parameters that will be learned by gradient-descent.

Algorithm 2: ClassifPI Input: the training distributions (µ n ) N -1 n=0 ; Output: estimates of optimal strategies (â n ) N -1 n=0 and probabilities p l (.; βn ); for n = N -1, . . . , 0 do Represent the discrete control at time n by neural network with parameter β n :

a n (x) = a n(x) with n (x) ∈ argmax =1,...,L p (x; β n ),
and compute the optimal parameter: βn ∈ argmin

β∈R q E L =1 p (X n ; β) f (X n , a ) + N -1 k=n+1 f X k , âk (X k ) + g(X N ) , (2.2)
where

X n ∼ µ n on R d , X n+1 = F (X n , a , ε n+1 ), X k+1 = F (X k , âk (X k ), ε k+1 ), for k = n + 1, .
. . , N -1 and = 1, . . . , L; Set ân (.) = a ˆ n(.) with ˆ n (x) ∈ argmax =1,...,L p (x; βn ); ân is the estimate of the optimal policy at time n Note that, when using Algorithms 1 and 2, the estimate of the optimal strategy at time n highly relies on the estimates of the optimal strategy at time n + 1, . . . , N -1, that have been computed previously. In particular, the practitioner who wants to use Algorithms 1 and 2 needs to keep track of the estimates of the optimal strategy at time n + 1, . . . , N -1 in order to compute the estimate of the optimal strategy at time n.

Remark 2.5 In practice, for n = N -1, ..., 0, one should minimize the expectations (2.1) and (2.2) by stochastic gradient-descent, where mini-batches of finite number of paths (X β k ) N k=n+1 are generated by drawing independent samples under µ n for the initial position at time n, and independent samples under ε k , for k = n + 1, . . . , N . The convergence of Algorithms 1 and 2 is analyzed in [START_REF] Huré | Deep neural networks algorithms for Stochastic Control Problems on finite horizon, part I: convergence analysis[END_REF] in terms of the error approximation of the optimal control by neural networks, and in terms of the estimation error by stochastic gradient descent methods, see their Theorem 4.7. 2

Control and value function learning by double DNN

We present in this section two algorithms, which in contrast with Algorithms 1 or 2, only keep track of the estimates of the value function and optimal control at time n + 1 in order to build an estimate of the value function and optimal control at time n.

Regress Now (Hybrid-Now)

The Algorithm 3, refereed to as Hybrid-Now, combines optimal policy estimation by neural networks and dynamic programming principle, and relies on an hybrid procedure between value and performance iteration.

Algorithm 3: Hybrid-Now

Input: the training distributions (µ n ) N -1 n=0 ; Output:
-estimate of the optimal strategy (â n ) N -1 n=0 ; -estimate of the value function ( Vn ) N -1 n=0 ; Set VN = g; for n = N -1, . . . , 0 do Compute:

βn ∈ argmin

β∈R q E f X n , A(X n ; β) + Vn+1 (X β n+1 ) (2.3)
where X n ∼ µ, and X β n+1 = F X n , A(X n ; β) , ε n+1 ); Set ân = A(.; βn ); ân is the estimate of the optimal policy at time n Compute θn ∈ argmin The Algorithm 4, called Hybrid-LaterQ, combines regress-later and quantization methods to build estimates of the value function. The main idea behind Algorithm 4 is to first interpolate the value function at time n + 1 by a set of basis functions, which is in the spirit of the regress-later-based algorithms, and secondly regress the interpolation at time n using quantization. The usual regress-later approach requires the ability to compute closed-form conditional expectations, which limits the stochastic dynamics and regression bases that can be considered. The use of quantization avoids this limitation and makes the regress-later algorithm more generally applicable. Let us first recall the basic ingredients of quantization. We denote by ε a K-quantizer of the R d -valued random variable ε n+1 ∼ ε 1 (typically a Gaussian random variable), that is a discrete random variable on a grid Γ = {e 1 , . . . , e K } ⊂ (R d ) K defined by

θ∈R p E (f (X n , ân (X n )) + Vn+1 (X βn n+1 ) -Φ(X n ; θ) 2 . (2.
ε = Proj Γ (ε 1 ) := K =1 e 1 ε 1 ∈C (Γ) ,
where C 1 (Γ), . . ., C K (Γ) are Voronoi tesselations of Γ, i.e., Borel partitions of the Euclidian space

(R d , |.|) satisfying C (Γ) ⊂ e ∈ R d : |e -e | = min j=1,...,K |e -e j | .
The discrete law of ε is then characterized by

p := P[ε = e ] = P[ε 1 ∈ C (Γ)], = 1, . . . , K.
The grid points (e ) which minimize the L 2 -quantization error ε 1ε 2 lead to the so-called optimal K-quantizer, and can be obtained by a stochastic gradient descent method, known as Kohonen algorithm or competitive learning vector quantization (CLVQ) algorithm, which also provides as a byproduct an estimation of the associated weights (p ). We refer to [START_REF] Pagès | Optimal quantization methods and applications to numerical problems in finance[END_REF] for a description of the algorithm, and mention that for the normal distribution, the optimal grids and the weights of the Voronoi tesselations are precomputed on the website http://www.quantize.maths-fi.com.

Algorithm 4: Hybrid-LaterQ Input:

-the training distributions (µ n ) N -1 n=0 ; -The grid {e 1 , . . . , e K } of K points in R d , with weights p 1 , . . . , p K for the quantization of the noise ε n ; Output:

-estimate of the optimal strategy (â

n ) N -1 n=0 ; -estimate of the value function ( Vn ) N -1 n=0 ; Set VN = g; for n = N -1, . . . , 0 do Compute: βn ∈ argmin β∈R q E f X n , A(X n ; β) + Vn+1 (X β n+1 ) (2.5)
where X n ∼ µ n , and X β n+1 = F X n , A(X n ; β) , ε n+1 ); Set ân = A(.; βn ); ân is the estimate of the optimal policy at time n Compute θn+1 ∈ argmin

θ∈R p E Vn+1 (X βn n+1 ) -Φ(X n+1 ; θ) 2 (2.6)
and set Ṽn+1 = Φ(.; θn+1 ); interpolation at time n + 1 Set

Vn (x) = f (x, ân (x)) + K =1 p Ṽn+1 F (x, ân (x), e ) ;

Vn is the estimate by quantization of the value function at time n

Quantization is mainly used in Algorithm 4 to efficiently approximate the expectations: recalling the dynamics (1.1), the conditional expectation operator for any functional W is equal to

P âM n (x) W (x) = E W (X âM n n+1 )|X n = x = E W (F (x, âM n (x), ε 1 )) , x ∈ R d ,
that we shall approximate analytically by quantization via:

P âM n (x) W (x) := E W (F (x, âM n (x), ε)) = K =1 p W F (x, âM n (x), e ) .
Observe that the solution to (2.6) actually provides a neural network Φ(.; θn+1 ) that interpolates Vn+1 . Hence the Algorithm 4 contains an interpolation step, and moreover, any kind of distance in R d can be chosen as a loss to compute θn+1 . In (2.6), we decide to take the L 2 -loss, mainly because it is the one that worked the best in our applications.

Remark 2.7 (Quantization) In dimension 1, we used the optimal grids and weights with K = 21 points, to quantize the reduced and centered normal law N (0, 1); and took 100 points to quantize the reduced and centered normal law in dimension 2, i.e. N 2 (0, 1). All the grids and weights for the optimal quantization of the normal law in dimension d are available in http://www.quantize.maths-fi.com for d = 1, . . . , 100. 2

Some remarks on Algorithms 3 and 4

As in Remark 2.5, all the expectations written in our pseudo-codes in Algorithm 3 and 4 should be approximated by empirical mean using a finite training set. The convergence of these algorithms has been analyzed in [START_REF] Huré | Deep neural networks algorithms for Stochastic Control Problems on finite horizon, part I: convergence analysis[END_REF] in terms of the approximation error of the optimal control and value function by neural networks, in terms of the estimation error by stochastic gradient descent methods, and in terms of the quantization error (for Algorithm 4, see their Theorems 4.14 and 4.19). Algorithms 3 or 4 are quite efficient to use in the usual case where the value function and the optimal control at time n are very close to the value function and the optimal control at time n + 1, which happens e.g. when the value function and the optimal control are approximations of the time discretization of a continuous in time value function and an optimal control. In this case, it is recommended to follow this two-step procedure:

(i) initialize the parameters (i.e. weights and bias) of the neural network approximations of the value function and the optimal control at time n to the ones of the neural network approximations of the value function and the optimal control at time n + 1.

(ii) take a very small learning rate parameter, for the Adam optimizer, that guarantees the stability of the parameters' updates from the gradient-descent based learning procedure.

Doing so, one obtains stable estimates of the value function and optimal control, which is desirable. We highlight the fact that this stability procedure is applicable here since the stochastic gradient descent method benefits from good initial guesses of the parameters to be optimized. It is an advantage compared to alternative methods proposed in the literature, such as classical polynomial regressions.

Quantization with k-nearest-neighbors (Qknn-algorithm)

Algorithm 5 presents the pseudo-code of an algorithm based on the quantization and k-nearest neighbors methods, called Qknn, which will be the benchmark in all the low-dimensional control problems that will be considered in Section 3 to test NNContPI, Clas-sifPI, Hybrid-Now and Hybrid-Later. Also, comparisons of Algorithm 5 to other well-known algorithms on various control problems in low-dimension are performed in [START_REF] Balata | A Class of Finite-Dimensional Numerically Solvable McKean-Vlasov Control Problems[END_REF], which show in particular that Algorithm 5 works very well to solve low-dimensional control problems. Actually, in our experiments, Algorithm 5 always outperforms the other algorithms based either on regress-now or regress-later methods whenever the dimension of the problem is low enough for Algorithm 5 to be feasible.

As done in Section 2.2.2, we consider a K-optimal quantizer of the noise ε n , i.e. a discrete random variable εn valued in a grid {e 1 , . . . , e K } of K points in E, and with weights p 1 , . . . , p K . We also consider grids Γ n , n = 0, . . . , N of points in R d , which are assumed to properly cover the region of R d that is likely to be visited by the optimally driven process X at time n = 0, . . . , N -1. These grids can be viewed as samples of wellchosen training distributions where more points are taken in the region that is likely to be visited by the optimally driven controlled process (see Remark 2.1 for details on the choice of the training measure).

Algorithm 5: Qknn Input: -Grids Γ k , k = 0, . . . , N in R d ; -Grid {e 1 , . . . , e K } of K points in E, with weights p 1 , . . . , p K for the quantization of ε n Output:
-estimate of the optimal strategy (â n ) N -1 n=0 ; -estimate of the value function ( Vn ) N -1 n=0 ; Set VN = g;

for n = N -1, . . . , 0 do Compute for (z, a) ∈ Γ n × A, Qn (z, a) = f (z, a) + K =1 p V n+1 Proj Γ n+1 F (z, a, e ) , (2.7)
where Proj Γ n+1 is the Euclidean projection over Γ n+1 ; Qn is the approximated Q-value c at time n Compute the optimal control at time n

Ân (z) ∈ argmin a∈A Qn (z, a) , ∀z ∈ Γ n ; (2.8)
use classical optimization algorithms of deterministic functions for this step Set V n (z) = Qn z, Ân (z) , ∀z ∈ Γ n ;

V n is the estimate by quantization of the value function c The Q-value at time n, denoted by Qn, is defined as the function that takes the couple state-action (x, a) as argument, and returns the expected optimal reward earned from time n to time N when the process X is at state x and action a is chosen at time n; i.e. Qn :

R d × R q ∈ (x, a) → f (x, a) + E a n,x [Vn+1(Xn+1)].
Remark 2.8 The estimate of the Q-value at time n given by (2.7) is not continuous w.r.t. the control variable a, which might cause some stability issues when running Qknn, especially during the optimization procedure (2.8). We refer to Section 3.2.2. in [START_REF] Balata | A Class of Finite-Dimensional Numerically Solvable McKean-Vlasov Control Problems[END_REF] for a detailed presentation of an extension of Algorithm 5 where the estimates of the Q value function Q n is continuous w.r.t. the control variable. 2

Numerical applications

In this section, we test the Neural-Networks-based algorithms presented in Section 2 on different examples. In high-dimension, we first took the same example as already considered in [START_REF] Weinan | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF] so that we can directly compare our results to theirs, and take another example from linear quadratic control problem with explicit analytic solution that is served as reference value. In low-dimension, we compared the results of our algorithms to the ones provided by Qknn, which has been introduced in Section 2 as an excellent benchmark for low-dimensional control problems.

A semilinear PDE

We consider the following semilinear PDE with quadratic growth in the gradient:

   ∂v ∂t + ∆ x v -|D x v| 2 = 0, (t, x) ∈ [0, T ) × R d , v(T, x) = g(x), x ∈ R d . (3.1)
By observing that for any p ∈ R d , -|p| 2 = inf a∈R d [|a| 2 + 2a.p], the PDE (3.1) can be written as a Hamilton-Jacobi-Bellman equation

   ∂v ∂t + ∆ x v + inf a∈R d |a| 2 + 2a.D x v] = 0, (t, x) ∈ [0, T ) × R d , v(T, x) = g(x), x ∈ R d , (3.2)
hence associated with the stochastic control problem

v(t, x) = inf α∈A E T t |α s | 2 ds + g(X t,x,α T ) , (3.3)
where X = X t,x,α is the controlled process governed by

dX s = 2α s ds + √ 2dW s , t ≤ s ≤ T, X t = x, W is a d-dimensional
X α n+1 = X α n + 2α n h + √ 2hε n+1 =: F (X α n , α n , ε n+1 ), n = 0, . . . , N -1,
where (ε n ) n is a sequence of i.i.d. random variables with law N (0, I d ), and the cost functional

J(α) = E N -1 n=0 h|α n | 2 + g(X α N ) .
On the other hand, it is known that an explicit solution to (3.1) (or equivalently (3.2)) can be obtained via a Hopf-Cole transformation (see e.g. [CR16]), and is given by

v(t, x) = -ln E exp -g(x + √ 2W T -t ) , (t, x) ∈ [0, T ] × R d . (3.4)
We choose to run tests on two different examples that have already been considered in the literature:

Test 1 Some recent numerical results have been obtained in [START_REF] Weinan | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF] (see Section 4.3 in [EHJ17]) when T = 1 and g(x) = ln( 1 2 (1 + |x| 2 )) in dimension d = 100 (see Table 2 and Figure 3 in [START_REF] Weinan | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF]). Their method is based on neural network regression to solve the BSDE representation associated with the PDE (3.1), and provide estimates of the value function at time 0 and state 0 for different values of a coefficient γ. We plotted the results of the Hybrid-Now algorithm in Figure 1. Hybrid-Now took one hour to achieves a relative error of 0.11%, using a 4-cores 3GHz intel Core i7 CPU. We want to highlight the fact that the algorithm presented in [START_REF] Weinan | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF] only needed 330 seconds to provide a relative error of 0.17%. However, in our experience, it is difficult to reduce the relative error from 0.17% to 0.11% using their algorithm. Also, we believe that the computation time of our algorithm can easily be reduced; some ideas in this direction are discussed in Section 4. The main trick that can be used is the transfer learning (also referred to as pre-training in the literature): we rely on the continuity of the value function and the optimal control w.r.t. time to claim that the value function and the optimal control at time n are very close to the ones at time n + 1. Hence, one can initialize the weights of the value function and optimal control at time n with the optimal ones estimated at step n + 1, reduce the learning rate of the optimizer algorithm, and reduce the number of steps for the gradient descent algorithm. All this procedure really speeds up the learning of the value function and the optimal control, and insures stability of the estimates. Doing so, we were able to reduce the computation time from one hour to twenty minutes.

We also considered the same problem in dimension d = 2, for which we plotted the first component of X w.r.t. time in Figure 2, for five different paths of the Brownian motion, where for each ω, the agent follows either the naive (α = 0) or the Hybrid-Now strategy. One can see that both strategies are very similar when the terminal time is far; but the Hybrid-Now strategy clearly forces X to get closer to 0 when the terminal time gets closer, in order to reduce the terminal cost.

Let us provide further implementation details on the algorithms presented in Test 1:

• As one can guess from the representation of v in (3.3), it is probably optimal to drive the process X around 0. Hence we decided to take µ n := ( nT N ) 1/2 N d (0, I d ) as a training measure at time n to learn the optimal strategy and value function at time n, for n = 0, . . . , N -1.

• We tested the algorithm with 1, 2 and 3 layers for the representation of the value function and the optimal control by neural networks, and noticed that the quality of the estimate significantly improves when using more than one layer, but does not vary significantly when considering more than 3 layers. the optimal strategy estimated by the Hybrid-Now (opt in blue) and the naive strategy α = 0 (bench in red). We consider the problem in dimension d=2. Observe that the optimal strategy (estimated by Hybrid-Now) is to do nothing when the terminal time is far in order to avoid any running cost, i.e. α opt = 0; and push X toward 0 when the terminal time is close, in order to minimize the terminal cost.

Test 2 Tests of the algorithms are proposed in dimension 1 with the terminal cost g(x) = -x γ 1 0≤x≤1 -1 1≤x and γ ∈ (0, 1). This problem was already considered in [START_REF] Richou | Etude théorique et numérique des équations différentielles stochastiques rétrogrades[END_REF], where the author proposed an algorithm based on a smart temporal discretization of the BSDE representation of the PDE (3.1) in order to deal with the quadratic growth of the driver of the BSDE, and usual projection on basis functions techniques for the approximation of conditional expectations that appear in the dynamic programming equation associated with the BSDE. We refer to equations (13),( 14),(15) in [START_REF] Richou | Numerical Simulation of BSDEs with Drivers of Quadratic Growth[END_REF] for details on the proposed algorithm, and its Theorem 4.14 for the convergence result. Their estimates of the value function at time 0 and state 0, when γ = 1, 0.5, 0.1, 0, are available in [START_REF] Richou | Etude théorique et numérique des équations différentielles stochastiques rétrogrades[END_REF], and have been reported in the column Y &R of Table 1. Also, the exact values for the value function have been computed for these values of γ by Monte Carlo using the closed-form formula (3.4), and are reported in the column Bench of Table 1. Tests of the Hybrid-Now and Hybrid-LaterQ algorithms have been run, and the estimates of the value function at time 0 and state x = 0 are reported in the Hybrid-Now and Hybrid-LaterQ columns. We also tested Qknn and reported its results in column Qknn. Note that Qknn is particularly wellsuited to 1-dimensional control problems. In particular, it is not time-consuming since the dimension of the state space is d=1. Actually, it provides the fastest results, which is not surprising since the other algorithms need time to learn the optimal strategy and value function through gradient-descent method at each time step n = 0, . . . , N -1. Moreover, Table 1 reveals that Qknn is the most accurate algorithm on this example, probably because it uses local methods in space to estimate the conditional expectation that appears in the expression of the Q-value.

γ

Y&R Hybrid-LaterQ Hybrid-Now Qknn Bench 1.0 -0.402 -0.456 -0.460 -0.461 -0.464 0.5 -0.466 -0.495 -0.507 -0.508 -0.509 0.1 -0.573 -0.572 -0.579 -0.581 -0.586 0.0 -0.620 -1.000 -1.000 -1.000 -1.000

We end this paragraph by giving some implementation details for the different algorithms as part of Test 2:

• Y&R: The algorithm Y&R converged only when using a Lipschitz version of g. The following approximation was used to obtain the results in Table 1:

g N (x) = g(x) if x ∈ [0, N -1 1-γ ] -N x otherwise.
• Hybrid-Now: We used N = 40 time steps for the time-discretization of [0, T ]. The value functions and optimal controls at time n = 0, . . . , N -1 are estimated using neural networks with 3 hidden layers and 10+5+5 neurons.

• Hybrid-LaterQ: We used N = 40 time steps for the time-discretization of [0, T ]. The value functions and optimal controls at time n = 0, . . . , N -1 are estimated using neural networks with 3 hidden layers containing 10+5+5 neurons; and 51 points for the quantization of the exogenous noise.

• Qknn: We used N = 40 time steps for the time-discretization of [0, T ]. We take 51 points to quantize the exogenous noise, ε n ∼ N (0, 1), for n = 0, . . . , N ; and decided to use the 200 points of the optimal grid of N 2 (0, 1) for the state space discretization.

The main conclusion regarding the results in this semilinear PDE problem is that Hybrid-Now provides better estimates of the solution to the PDE in dimension d=100 than the previous results available in [START_REF] Weinan | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF] but requires more time to do so.

Hybrid-Now and Hybrid-Later provide better results than those available in [START_REF] Richou | Numerical Simulation of BSDEs with Drivers of Quadratic Growth[END_REF] to solve the PDE in dimension 2; but are outperformed by Qknn, which is arguably very accurate.

A linear quadratic stochastic test case

We consider a linear controlled process with dynamics in R d according to

dX t = (BX t + Cα t )dt + p j=1 D j α t dW j t , (3.5)
where W j , j = 1, . . . , p, are independent real-Brownian motion, the control process α ∈ A is valued in R m , and the constant coefficients B ∈ R d×d , C, D j ∈ R d×m , j = 1, . . . , p. The value function of the linear quadratic stochastic control problem is

v(t, x) = inf α∈A E T t (X t,x,α s .QX t,x,α s + λ|α t | 2 )dt + X t,x,α T .P X t,x,α T , (t, x) ∈ [0, T ] × R d ,
where X t,x,α is the solution to (3.5) starting from x at time t, given a control process α ∈ A, P, Q are nonnegative symmetric d × d matrices, and λ > 0. The Bellman equation associated with this stochastic control problem is a fully nonlinear equation in the form

∂v ∂t + x.Qx + inf a∈R (Bx + Ca).D x v + a λI m + p j=1 D j D 2 x vD j 2 a = 0, on [0, T ) × R d , v(T, x) = x.P x, x ∈ R d ,
and it is well-known, see e.g. [START_REF] Yong | Stochastic Controls Hamiltonian Systems and HJB Equations[END_REF], that an explicit solution is given by

v(t, x) = x.K(t)x, (3.6)
where K(t) is a nonnegative symmetric d × d matrix, solution to the Riccati equation

K + B K + KB + Q -KC(λI m + p j=1 D j KD j ) -1 C K = 0, K(T ) = P, (3.7)
while an optimal feedback control is equal to

a * (t, x) = -λI m + p j=1 D j K(t)D j -1 C K(t)x, (t, x) ∈ [0, T ) × R d . (3.8)
We numerically solve this problem by considering a time discretization (with time step h = T /N ), which leads to the discrete-time control problem with dynamics

X α n+1 = X α n + (BX α n + Cα n )h + Dα n √ hε n+1 =: F (X α n , α n , ε n+1 ), n = 0, . . . , N -1,
where (ε n ) n is a sequence of i.i.d. random variables with law N (0, 1), and cost functional

J(α) = E N -1 n=0 X α n .QX α n + λ|α n | 2 )h + X α N .P X α N .
For the numerical tests, we take m = 1, p = d, and the following parameters:

T = 1, N = 20, B = I d , C = 1 d , D j = (0, . . . , 1 j-th term
, . . . , 0) , j = 1, . . . , p,

Q = P = I d , λ = 1,
where we denote

1 d := ( 1, . . . , 1 d times
) .

Numerical results

We implement our algorithms in dimension d = 1, 10, 100, and compare our solutions with the analytic solution via the Riccati equation (3.7) solved by Matlab d .

• For d = 1, we plotted the estimates of the optimal control at time n = 0, . . . , N -1 in Figure 3 and the value function in Figure 4. Observe that, as expected, the estimated optimal control is linear and the estimated value function is quadratic at each time.

• For d = 10, we reported in Table 2 the estimates of v(0, X 0 ), computed by running forward simulations of X using the estimated optimal strategy. "Riccati" is v(0, X 0 ) computed by solving (3.7) with Matlab. We set the initial position to X 0 = 1 d . We also plotted in Figure 5 a forward simulation of the components of X optimally controlled.

Observe that NNContPI is more accurate than Hybrid-Now. Notice that the estimates provided by the algorithms are biased, which is due to the time discretization.

• For d = 100, we reported in Table 3 the estimates of the value function, computed by running forward simulations of X using the estimated optimal strategy. "Riccati" is v(0, X 0 ) computed by solving (3.7) with Matlab. We set the initial position to X 0 = 0.11 d and X 0 = 0.51 d . Once again, NNContPI is slightly more accurate than Hybrid-Now, and the estimates provided by the latter are biased due to the time discretization.

Implementation details:

We implemented Hybrid-Now and NNContPI using training sets from the distribution µ n := N d (0, 1) for n = 0, . . . , N -1. We represented the value function and optimal control at time n, n = 0, . . . , N -1 using two hidden layers with d+20 and d+10 neurons, and 1 neuron for the output layers. We used Elu as activation function for the hidden layers, and identity for the output layer. Comments on the algorithms: Hybrid-Now behaved similarly as for the SemiLinear PDE example, and we can make the same remarks. NNContPI is much slower than Hybrid-Now, because the data have to go through the N -n-1 neural networks that represent the optimal controls at time n + 1, . . . , N -1, in order to estimate the optimal control at time n. The bench curve represents one of the identical component of X when it is driven using the strategy α = 0. One can see that the optimal control tends to reduce the norm of each component of X.

Table 2: Estimate of v(0, X 0 ) obtained by forward simulation of the process controlled by the optimal strategy estimated by Hybrid-Now and NNContPI. "Riccati" is v(0, X 0 ) computed by solving (3.7) with Matlab. We took d = 10, and X 0 = 1 d . Mean and standard deviation are computed on 10 sets of 10,000 simulations each.

Mean std

Option hedging

Our third example comes from a classical hedging problem in finance. We consider an investor who trades in q stocks with (positive) price process (P n ) n , and we denote by (α n ) valued in A ⊂ R q the amount held in these assets over the period (n, n + 1]. We assume for simplicity that the price of the riskless asset is constant equal to 1 (zero interest rate).

It is convenient to introduce the return process as: R n+1 = diag(P n ) -1 (P n+1 -P n ), n = 0, . . . , N -1, so that the self-financed wealth process of the investor with a portfolio strategy α, and starting from some capital w 0 , is governed by

W α n+1 = W α n + α n .R n+1 , n = 0, . . . , N -1, W α 0 = w 0 .
Given an option payoff h(P N ), the objective of the agent is to minimize over her portfolio strategies α her expected square replication error

V 0 = inf α∈A E h(P N ) -W α N ,
where is a convex function on R. Assuming that the returns R n , n = 1, . . . , N are i.i.d, we are in a (q + 1)-dimensional framework of Section 1 with X α = (W α , P ) with ε n = R n valued in E ⊂ R q , with the dynamics function

F (w, p, a, r) = w + a.r p + diag(p)r, x = (w, p) ∈ R × R q , a ∈ R q , r ∈ E,
the running cost function f = 0 and the terminal cost g(w, p) = (h(p)w). We test our algorithm in the case of a square loss function, i.e. (w) = w 2 , and when there is no portfolio constraints A = R q , and compare our numerical results with the explicit solution derived in [START_REF] Bertsimas | Hedging derivative securities and incomplete markets: an ε-arbitrage approach[END_REF]: denote by ν(dr) the distribution of R n , by ν = E[R n ] = rν(dr) its mean, and by M2 = E[R n R n ] assumed to be invertible; we then have

V n (w, p) = K n w 2 -2Z n (p)w + C n (p)
where the functions K n > 0, Z n (p) and C n (p) are given in backward induction, starting from the terminal condition

K N = 1, Z N (p) = h(p), C N (p) = h 2 (p),
and for n = N -1, . . . , 0, by

K n = K n+1 1 -ν M -1 2 ν , Z n (p) = Z n+1 (p + diag(p)r)ν(dr) -ν M -1 2 Z n+1 (p + diag(p)r)rν(dr), C n (p) = C n+1 (p + diag(p)r)ν(dr) - 1 K n+1 Z n+1 (p + diag(p)r)rν(dr) M -1 2 Z n+1 (p + diag(p)r)rν(dr) , so that V 0 = K 0 w 2 0 -2Z 0 (p 0 )w 0 + C 0 (p 0 )
, where p 0 is the initial stock price. Moreover, the optimal portfolio strategy is given in feedback form by α * n = a * n (W * n , P n ), where a * n (w, s) is the function

a * n (w, p) = M -1 2 Z n+1 (p + diag(p)r)rν(dr) K n+1 -νw ,
and W * is the optimal wealth associated with α * , i.e., W * n = W α * n . Moreover, the initial capital w * 0 that minimizes V 0 = V 0 (w 0 , p 0 ), and called (quadratic) hedging price is given by

w * 0 = Z 0 (p 0 ) K 0 .
Test Take N = 6, and consider one asset q = 1 with returns modeled by a trinomial tree:

ν(dr) = π + δ r + + π 0 δ 0 + π -δ r -, π 0 + π + + π -= 1,
with r + = 5%, r -= -5%, π + = 60%, π -= 30%. Take p 0 = 100, and consider the call option h(p) = (pκ) + with κ = 100. The price of this option is defined as the initial value of the portfolio that minimizes the terminal quadratic loss of the agent when the latter follows the optimal strategy associated with the initial value of the portfolio. In this test, we want to determine the price of the call and the associated optimal strategy using different algorithms.

Remark 3.1 The option hedging problem is linear-quadratic, hence belongs to the class of problems where the agent has ansatzes on the optimal control and the value function. Indeed, we expect here the optimal control to be affine w.r.t. w and the value function to be quadratic w.r.t. w. For these kind of problems, the algorithms presented in Section 2 can easily be adapted so that the expressions of the estimators satisfy the ansatzes. See (3.9) and (3.10) for the option hedging problem. 2

Numerical results

In Figure 6, we plot the value function at time 0 w.r.t w 0 , the initial value of the portfolio, when the agent follows the theoretical optimal strategy (benchmark), and the optimal strategy estimated by the Hybrid-Now or Hybrid-LaterQ algorithms. We perform forward Monte Carlo using 10,000 samples to approximate the lower bound of the One can observe that all the algorithms estimate the price to be 4.5, but Hybrid-LaterQ is better than Hybrid-Now at reducing the quadratic risk.

value function at time 0 (see [START_REF] Henry-Labordere | Deep Primal-Dual Algorithm for BSDEs: Applications of Machine Learning to CVA and IM[END_REF] for details on how to get an approximation of the upperbound of the value function via duality). One can observe that while all the algorithms give a call option price approximately equal to 4.5, Hybrid-LaterQ clearly provides a better strategy than Hybrid-Now to reduce the quadratic risk of the terminal loss. We plot in Figure 7 three different paths of the value of the portfolio w.r.t the time n, when the agent follows either the theoretical optimal strategy (red), or the estimated one using Hybrid-Now (blue) or Hybrid-LaterQ (green). We set w 0 = 100 for these simulations.

Comments on Hybrid-Now and Hybrid-LaterQ

The Option Hedging problem belongs to the class of linear-quadratic control problems for which we expect the optimal control to be affine w.r.t. w and the value function to be quadratic w.r.t. w. It is then natural to consider the following classes of controls A M and functions F M to properly approximate the optimal controls and the values functions at time n=0, . . . , N -1:

A M := (w, p) → A(x; β) • 1, w ; β ∈ R p , (3.9)
F M := (w, p) → Φ(x; θ) • 1, w, w 2 ; θ ∈ R p , (3.10)
where β describes the parameters (weights+bias) associated with the neural network A and θ describes those associated with the neural network Φ. The notation stands for the transposition, and • for the inner product. Note that there are 2 (resp. 3) neurons in the output layer of A (resp. Φ), so that the inner product is well-defined in (3.10) and (3.9).

Hybrid-Now

Hybrid-LaterQ

Figure 7: Three simulations of the agent's wealth w.r.t. time n when, for each ω, the latter follows the theoretical optimal strategy (red), the estimated one using Hybrid-Now (blue) and the one using Hybrid-LaterQ (green). We took w 0 = 100. Observe that the process is driven similarly to the optimally controlled process, when the agent follows the estimated optimal strategy using Hybrid-LaterQ or Hybrid-Now.

Valuation of energy storage

We present a discrete-time version of the energy storage valuation problem studied in [START_REF] Carmona | Valuation of energy storage: an optimal switching approach[END_REF]. We consider a commodity (gas) that has to be stored in a cave, e.g. salt domes or aquifers. The manager of such a cave aims to maximize the real options value by optimizing over a finite horizon N the dynamic decisions to inject or withdraw gas as time and market conditions evolve. We denote by (P n ) the gas price, which is an exogenous real-valued Markov process modeled by the following mean-reverting process:

P n+1 = p(1 -β) + βP n + ξ n+1 , (3.11)
where β < 1, and p > 0 is the stationary value of the gas price. The current inventory in the gas storage is denoted by (C α n ) n and depends on the manager's decisions represented by a control process α = (α n ) valued in {-1, 0, 1}: α n = 1 (resp. -1) means that she injects (resp. withdraws) gas with an injection (resp. withdrawal) rate

a in (C α n ) (resp. a out (C α n )) requiring (causing) a purchase (resp. sale) of b in (C α n ) ≥ a in (C α n ) (resp. b out (C α n ) ≤ a out (C α n ))
, and α n = 0 means that she is doing nothing. The difference between b in and a in (resp. b out and a out ) indicates gas loss during injection/withdrawal. The evolution of the inventory is then governed by

C α n+1 = C α n + h(C α n , α n ), n = 0, . . . , N -1, C α 0 = c 0 , (3.12)
where we set

h(c, a) =     
a in (c) for a = 1 0 for a = 0 -a out (c) for a = -1, and we have the physical inventory constraint:

C α n ∈ [C min , C max ], n = 0, . . . , N.
The running gain of the manager at time n is f (P n , C α n , α t ) given by

f (p, c, a) =      -b in (c)p -K 1 (c) for a = 1 -K 0 (c) for a = 0 b out (c)p -K -1 (c) for a = -1,
and K i (c) represents the storage cost in each regime i = -1, 0, 1. The problem of the manager is then to maximize over α the expected total profit

J(α) = E N -1 n=0 f (P n , C α n , α n ) + g(P N , C α N ) , (3.13)
where a common choice for the terminal condition is

g(p, c) = -µp(c 0 -c) + ,
which penalizes for having less gas than originally, and makes this penalty proportional to the current price of gas (µ > 0). We are then in the 2-dimensional framework of Section 1 with X α = (P, C α ), and the set of admissible controls in the dynamic programming loop is given by:

A n (c) = a ∈ {-1, 0, 1} : c + h(c, a) ∈ [C min , C max ], c ∈ [C min , C max ] , n = 0, . . . , N -1.
Test We fixed the parameters as follows, to run our numerical tests:

a in (c) = b in (c) = 0.06, a out (c) = b out (c) = 0.25 K i (c) = 0.01c
C max = 8, C min = 0, c 0 = 4, p = 5, β = 0.5, ξ n+1 ; N (0, σ 2 ) with σ 2 = 0.05, and µ = 2 in the terminal penalty function, N =30.

Numerical results

We plotted in Figure 8 the estimates of the value function at time 0 w.r.t. a in using Qknn, as well as the reward function (3.13) associated with the naive do-nothing strategy α = 0 (see Bench in figure 8). As expected, the naive strategy performs well when a in is small compared to a out , since, in this case, it takes time to fill the cave, so that the agent is likely to do nothing in order to avoid any penalization at terminal time. When a in is of the same order as a out , it is easy to fill up and empty the cave, so the agent has more freedom to buy and sell gas in the market without worrying about the terminal cost. Observe that the value function is not monotone, due to the fact that the C component in the state space takes its value in a bounded and discrete set (see (3.12)). We also plotted the cost functional associated with the naive passive strategy α = 0 (Bench). See that for small values of a in such as 0.06, doing nothing is a reasonable strategy. Observe also that the value function is not monotone w.r.t. a in which is due to the dynamics of C (3.12).

Table 4 provides the estimates of the value function using the ClassifPI, Hybrid-Now and Qknn algorithms. Observe first that the estimates provided by Qknn are larger than those provided by the other algorithms, meaning that Qknn outperforms the other algorithms. The second best algorithm is ClassifPI, while Hybrid-Now performs poorly and clearly suffers from instability, due to the discontinuity of the running rewards w.r.t. the control variable. As expected, one can observe on each plot that the optimal strategy is to inject gas when the price is low, to sell gas when the price is high, and to make sure to have a volume of gas greater than c 0 in the cave when the terminal time is getting closer to minimize the terminal cost.

Let us now comment on the implementation of the algorithms:

• Qknn: Table 4 shows that once again, due to the low-dimensionality of the problem, Qknn provides the best value function estimates. The estimated optimal strategies, shown on Figure 9, are very good estimates of the theoretical ones. The three decision regions on Figure 9 are natural and easy to interpret: basically it is optimal to sell when the price is high, and to buy when it is low. However, a closer look reveals that the waiting region (where it is optimal to do nothing) has an unusual triangular-based shape, due essentially to the discreteness of the space on which the C component of the state space takes its values. We expect this shape to be very hard to reproduce with the DNN-based algorithms proposed in Section 2.

• ClassifPI: As shown on Figure 10, the ClassifPI algorithm manages to provide accurate estimates for the optimal controls at time n = 0, . . . , N -1. However, the latter is not able to catch the particular triangular-based shape of the waiting region, which explains why Qknn performs better.

• Hybrid-Now: As shown on Figure 11, Hybrid-Now only manages to provide relatively poor estimates, compared to ClassifPI and Qknn, of the three different regions at time n = 0, . . . , N -1. In particular, the regions suffer from instability.

We end this paragraph by providing some implementation details for the different algorithms we tested.

• Qknn: We used the extension of Algorithm 5 introduced in the paragraph "semi-linear interpolation" of the Section 3.2.2. in [START_REF] Balata | A Class of Finite-Dimensional Numerically Solvable McKean-Vlasov Control Problems[END_REF] and used a projection of each state on its k=2-nearest neighbors to get an estimate of the value function which is continuous w.r.t. the control variable at each time n = 0, . . . , N -1. The optimal control is computed at each point of the grids using the Brent algorithm, which is a deterministic function optimizer already implemented in Python e .

• Implementation details for the neural network-based algorithms: We use neural networks with two hidden layers, ELU activation functions f and 20+20 neurons . The output layer contains 3 neurons with softmax activation function for the ClassifPI algorithm and no activation function for the Hybrid-Now one. We use a training set of size M=60,000 at each time step. Note that given the expression of the terminal cost, the ReLU activation functions (Rectified Linear Units) could have been deemed a better choice to capture the shape of the value functions, but our tests revealed that ELU activation functions e We could have chosen other algorithms to optimize the Q-value, but, in our tests, Brent was faster than the other choices that we tried, such as GoldenSearch, and always provided accurate estimates of the optimal controls. provide better results. At time n = 0, . . . , N -1, we took µ n = U(C min , C max ) as training measure.

We did not use the pre-train trick discussed in Section 2.2.3, which explains the instability in the decisions that can be observed in Figure 11.

The main conclusion of our numerical comparisons on this energy storage example is that ClassifPI, the DNN-based classification algorithm designed for stochastic control problems with discrete control space, appears to be more accurate than the more general Hybrid-Now. Nevertheless, ClassifPI was not able to capture the unusual triangle-based shape of the optimal control as well as Qknn did.

Microgrid management

Finally, we consider a discrete-time model for power microgrid inspired by the continuoustime models developed in [START_REF] Benjamin Heymann | Continuous optimal control approaches to microgrid energy management[END_REF] and [START_REF] Daniel | An approximate dynamic programming algorithm for monotone value functions[END_REF]; see also [START_REF] Alasseur | Regression Monte Carlo for Microgrid Management[END_REF]. The microgrid consists of a photovoltaic (PV) power plant, a diesel generator and a battery energy storage system (BES), hence using a mix of fuel and renewable energy sources. These generation units are decentralized, i.e., installed at a rather small scale (a few kW power), and physically close to electricity consumers. The PV produces electricity from solar panels with a generation pattern (P n ) n depending on the weather conditions. The diesel generator has two modes: on and off. Turning it on consumes fuel, and produces an amount of power α n . The BES can store energy for later use but has limited capacity and power. The aim of the microgrid management is to find the optimal planning that meets the power demand, denoted by (D n ) n , while minimizing the operational costs due to the diesel generator. We denote by R n = D n -P n , the residual demand of power: when R n > 0, one should provide power through diesel or battery, and when R n < 0, one can store the surplus power in the battery.

The optimal control problem over a fixed horizon N is formulated as follows. At any time n = 0, . . . , N -1, the microgrid manager decides the power production of the diesel generator, either by turning it off: α n = 0, or by turning it on, hence generating a power α n valued in [A min , A max ] with 0 < A min < A max < ∞. There is a fixed cost κ > 0 associated with switching from the on/off mode to the other one off/on, and we denote by M α n the mode valued in {0 = off, 1 = on} of the generator right before time n, i.e., M α n+1 = 1 αn =0 . When the diesel generator and renewable provide a surplus of power, the excess can be stored into the battery (up to its limited capacity) for later use, and in case of power insufficiency, the battery is discharged for satisfying the power demand. The input power process I α for charging the battery is then given by

I α n = (α n -R n ) + ∧ (C max -C α n ),
where C max is the maximum capacity of the battery with current charge C α , while the output power process O α for discharging the battery is given by

O α n = (R n -α n ) + ∧ C α n .
Remark 3.2 The state/space constraint is managed in our NN-based algorithm by introducing a penalty function into the running cost (see Remark 2.4): f (x, a) ← f (x, a)+L(x, a)

L(x, a) = Q + rca + with large Q + taken much larger than Q -. Doing so, the NN-based estimate of the optimal control learns not to take any forbidden decision. 2

The control space {0} ∪ [A min , A max ] is a mix between a discrete space and a continuous space, which is challenging for algorithms with neural networks. We actually use a mixture of classification and standard DNN for the control: (p 0 (x; θ), π(x; β)) valued in [0, 1] × [A min , A max ], where p 0 (x; θ) is the probability of turning off in state x, and π(x; β) is the amount of power when turning on with probability 1p 0 (x; θ). In other words,

X n+1 = F (X n , 0, ε n+1 ) with probability p 0 (X n ; θ n ) F (X n , π(X n ; β n ), ε n+1 ) with probability 1 -p 0 (X n ; θ n )
The pseudo-code of this approach, specifically designed for this problem, is written in Algorithm 6, and we henceforth refer to it as ClassifHybrid. Note in particular that it is an Hybrid version of ClassifPI.

( β0 n , β1 n ) ∈ argmax β 0 ,β 1 E p 0 (X n ; β 0 ) f (X n , 0) + Vn+1 f ( X0 n+1 + (1 -p 0 (X n ; β 0 )) f (X n , π(X n ; β 1 )) + Vn+1 X1,β 1 n+1 , where X n ; µ n , X0 n+1 = F (X n , 0, ε n+1 ), and X1,β 1 n+1 = F (X n , π(X n ; β 1 ), ε n+1 ); Compute θn ∈ argmin θ E p 0 X n ; β0 n f (X n , 0) + Vn+1 f ( X0 n+1 -Φ(.; θ) 2 + 1 -p 0 X n ; β0 n f (X n , π(X n ; β 1 n )) + Vn+1 X1, β1
C min = 0, C max = 1 or 4, C 0 = 0, K = 2, γ = 2, κ = 0.2, Q -= 10, R 0 = 0.1, A min = 0.05, A max = 10 Q + = 1000.
Results Figure 12 shows the Qknn-estimated optimal decisions to take at times n = 1, 10, 28 in the cases where m = M n = 0 and m = M n = 1. If the generator is off at time n, i.e. m = 0, the blue curve separates the region where it is optimal to keep it off and the one where it is optimal to generate power. If the generator is on at time n, i.e. m = 1, the blue curve separates the region where it is optimal to turn it off and the one where it is optimal to generate power. A colorscale is available on the right to inform how much power it is optimal to generate in both cases. Observe that the optimal decisions are quite intuitive: for example, if the demand is high and the battery is empty, then it is optimal to generate a lot of energy. Moreover, it is optimal to turn the generator off if the demand is negative or if the battery is charged enough to meet the demand. We plot in Figure 13 the estimated optimal decisions at times n = 1, 10, 28, using the Hybrid-Now algorithm, with N = 30 time steps. See that the decisions are similar to the ones given using Qknn.

Note that the plots in Figure 12 and 13 look much better than the ones obtained in [START_REF] Alasseur | Regression Monte Carlo for Microgrid Management[END_REF] in which algorithms based on regress-now or regress-later are used (see in particular Figure 4 in [START_REF] Alasseur | Regression Monte Carlo for Microgrid Management[END_REF]); hence Qknn and ClassifHybrid seem more stable than the algorithms proposed in [START_REF] Alasseur | Regression Monte Carlo for Microgrid Management[END_REF]. We report in Table 5 the result for the estimates of the value function with N =30 time steps, obtained by running 10 times a forward Monte Carlo with 10,000 simulations using the optimal strategy estimated using Qknn and ClassifHybrid algorithms. Observe that Hybrid-Now performs better than Qknn. However, Qknn run in less than a minute whereas Hybrid-Now needed seven minutes to run. We also report in Table 6 the value function estimates with N =200 time steps, obtained by running 20 times a forward Monte Carlo with 10,000 simulations using the Qknn-estimated optimal strategy. Table 5: Estimates of the value function at time 0 and state (C 0 = 0, M 0 = 0, R 0 = 0.1), for N = 30 and C max = 1, using Qknn and ClassifHybrid algorithms. Note that ClassifHybrid achieved better results than Qknn on this problem. The region under the blue line is the one where it is optimal to turn the generator off if m=1 (i.e. the generator was on at time n-1), or keep it off if m = 0 (i.e. the generator was off at time n-1). Figure 14 shows two simulations of (C, M, R) controlled using the Qknn-estimated optimal strategy, where N = 200 has been chosen. Observe in particular the natural behavior of the Qknn-decisions which consists in turning the generator on when the demand cannot be met by the battery, and turn it off when the demand is negative or when the battery is charged enough to meet the demand. Note that the plots are similar to the ones plotted in Figure 9 of [START_REF] Alasseur | Regression Monte Carlo for Microgrid Management[END_REF].

Comments on Qknn: Note that there is no need to use a penalization method with the Qknn-algorithm to constrain the control to stay in A n (x), where x is the state at time n, since, for all state x, we can simply search for the optimal control associated in A n (x), using e.g. the Brent algorithm. For n = 0, . . . , N -1, we took the training set as follows: Γ n := Γ C × {0, 1} × Γ n R ; where Γ C := {C min + i 50 (C max -C min ), i = 0, . . . , 50}, Γ n R := ρ n R 0 + σ 1-ρ n 1-ρ Γ 1 and where Γ 1 is the optimal grid for the quantization of N (0, 1), available in http://www.quantize.maths-fi.com, with 51 points. This choice of training points for the C component corresponds to the exploration procedure discussed in Remark 2.1, whereas we chose the best grid with 51 points for the (uncontrolled) R component. Comments on ClassifHybrid: We took 100 mini-batches of size 300 and took 100 epochs to run the algorithm. We chose the following training distribution at time n: µ n = U(C min , C max ) × U({0, 1}) × P Rn , where P Rn is the law of the (uncontrolled) residual demand at time t n . Note that such a choice of training distribution means that we want to explore all the available states for the controlled components of the controlled process (C, M, R) in order to learn the optimal strategy globally.

The microgrid management problem is very challenging for our algorithms because the control space {0} ∪ [a min , a max ] is a mix of discrete and continuous space, moreover the choice of the optimal control is subject to constraints. We designed ClassifHybrid, an Hybrid version of ClassifPI, to solve this problem. ClassifHybrid provided very good estimates and actually managed to perform better than Qknn.

Discussion and conclusion

Our proposed algorithms are well-designed and provide accurate estimates of optimal control and value function associated with various high-dimensional control problems. Also, when tested on low-dimensional problems, they performed as well as the Monte Carlo-based or quantization-based methods, which have shown their efficiency in low dimension, see e.g. [START_REF] Balata | A Class of Finite-Dimensional Numerically Solvable McKean-Vlasov Control Problems[END_REF] and [START_REF] Alasseur | Regression Monte Carlo for Microgrid Management[END_REF].

The presented algorithms suffer from a rather high time-consuming cost due to the expensive training of 2(N -1) neural networks to learn the value functions and optimal controls at times n = 0, . . . , N -1. However, the agent can easily alleviate the computation time. A first trick consists in reducing the number of neural networks by partially or totally ignoring the dynamic programming principle (DPP), as it has been done e.g. in [START_REF] Weinan | Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations[END_REF]. The use of one unique Recurrent Neural Networks (RNN) (in the case where the DPP is totally ignored) or a few of them (in the partial-ignored case) can also be considered to learn the optimal controls, either all at the same time (first case), or group by group in a backward way (second case). We refer to [START_REF] Chan Wai-Nam | Machine Learning for semi linear PDEs[END_REF] for algorithms in this spirit. Another trick consists in learning faster the value functions and optimal controls at times n = 0, . . . , N -1 by pre-training the neural networks. The way to proceed in that direction is to initialize at time n the weights and bias of the value function estimator Vn to the ones of Vn+1 . We then rely on the continuity of the value function w.r.t. the time n to expect that the weights will not change much from time n to n + 1, hence trainable very quickly by reducing the learning rate of the Adam algorithm for the gradient descent, and using an early-stop procedure as implemented in Keras g . Another benefit from the pre-training task is to get the stability of the estimates w.r.t. time, which is also a pleasant feature.

  Brownian motion, and the control process α is valued in A = R d . The time discretization (with time step h = T /N ) of the control problem (3.3) leads to the discrete-time control problem (1.1)-(1.2)-(1.3) with

Figure 1 :Figure 2 :

 12 Figure 1: Relative error of the Hybrid-Now estimate of the value function at time 0 w.r.t the number of mini-batches used to build the Hybrid-Now estimators of the optimal strategy. The value functions have been computed running three times a forward Monte Carlo with a sample of size 10,000, following the optimal strategy estimated by the Hybrid-Now algorithm.

Table 1 :

 1 Value function at time 0 and state 0 w.r.t. γ, computed with the Y&R, Hybrid-Now, Hybrid-Later and Qknn algorithms. Bench reports the MC estimates of the closed-form formula (3.4).

d

  We solved (3.7) with the Matlab method ode45.

Figure 3 :

 3 Figure 3: Optimal decision estimated by Hybrid-Now at time n = 0, . . . , N -1. We took d = 1, N = 20. We observe that the estimates are linear, as expected given the closed-form formula (3.8) for the optimal control.

Figure 4 :

 4 Figure 4: Value function w.r.t. x, estimated by Hybrid-Now at time n = 0, . . . , N -1. We took d = 1, N = 20. We observe that the estimates are quadratic, as expected given the closed-form formula (3.6) for the value function.

Figure 5 :

 5 Figure5: Forward simulation of X w.r.t. time, when X 0 = 1 d and d = 10, driven optimally using Hybrid estimates. The first ten curves represent the ten components of X. The bench curve represents one of the identical component of X when it is driven using the strategy α = 0. One can see that the optimal control tends to reduce the norm of each component of X.

  Estimate of v(0, X 0 ) obtained by forward simulation of the process controlled by the optimal strategy estimated by Hybrid-Now and NNContPI. "Riccati" is v(0, X 0 ) computed by solving (3.7) with Matlab. We took d = 100, and initial position X 0 = 0.51 d and X 0 = 0.11 d . Mean and standard deviation are computed on 10 sets of 10,000 simulations each.

Figure 6 :

 6 Figure6: Estimates of the value function at time 0 w.r.t. w 0 using Hybrid-Now (blue line) or Hybrid-LaterQ (green dashes). We draw the value function in red for comparison. One can observe that all the algorithms estimate the price to be 4.5, but Hybrid-LaterQ is better than Hybrid-Now at reducing the quadratic risk.

Figure 8 :

 8 Figure8: Estimate of the value function at time 0 w.r.t. a in , when the agent follows the optimal strategy estimated by Qknn, by running a forward Monte Carlo with a sample of size 100,000 (blue). We also plotted the cost functional associated with the naive passive strategy α = 0 (Bench). See that for small values of a in such as 0.06, doing nothing is a reasonable strategy. Observe also that the value function is not monotone w.r.t. a in which is due to the dynamics of C (3.12).

Table 4 :Finally,

 4 V (0, P 0 , C 0 ) estimates for different values of a in , using the optimal strategy provided by the ClassifPI , Hybrid-Now and Qknn algorithms, with a out = 0.25, P 0 = 4 and C 0 = 4. Figures 9, 10, 11 provide the optimal decisions w.r.t. (P, C) at times 5, 10, 15, 20, 25, 29 estimated respectively by the Qknn, ClassifPI and Hybrid-Now algorithms.

Figure 9 :

 9 Figure 9: Estimated optimal decisions at times 5, 10, 15, 20, 25, 29 w.r.t. (P,C) for the energy storage valuation problem using Qknn. Injection (a=-1) in red, store (a=0) in black and withdraw (a=1) in blue.

Figure 10 :

 10 Figure 10: Estimated optimal decisions at times 5, 10, 15, 20, 25, 29 w.r.t. (P,C) for the energy storage valuation problem using ClassifPI. Injection (a=-1) in purple, store (a=0) in blue and withdraw (a=1) in yellow.

Figure 11 :

 11 Figure 11: Estimated optimal decisions at times 5, 10, 15, 20, 25, 29 w.r.t. (P,C) for the energy storage valuation problem using Hybrid-Now. Injection (a=-1) in purple, store (a=0) in blue and withdraw (a=1) in yellow. Observe the instability in the decisions which come from the fact that we did not pre-train the neural networks (see Section 2.2.3)

Algorithm 6 :

 6 ClassifHybridInput: the training distributions (µ n ) N -1 n=0 ; Output: -estimate of the optimal strategy (â n ) N -1 n=0 ; -estimate of the value function ( Vn ) N -1 n=0 ; Set VN = g; for n = N -1, . . . , 0 do Compute

n

  n+1 -Φ(.; θ) 2 ; Set Vn = Φ(.; θn ); Vn is the estimate of the value function at time n Test We set the parameters to the following values to compare Qknn and ClassifHybrid:

Figure 12 :

 12 Figure 12: Estimated optimal decisions at time 1, 10 and 28, using Qknn, with N = 30 time steps.

Figure 13 :Figure 14 :

 1314 Figure 13: Estimated optimal decisions at time 1, 10 and 28, using ClassifHybrid, with N = 30 time steps.

f The Exponential Linear Unit (ELU) activation function is defined asx → exp(x) -1 if x ≤ 0 x if x > 0 .

The author's research is carried out with support of the Norwegian Research Council, within the research project Challenges in Stochastic Control, Information and Applications (STOCONINF), project The work of this author is supported by the ANR project CAESARS (ANR-15-CE05-0024), and also by FiME and the "Finance and Sustainable Development" EDF -CACIB Chair 1

Here, we denote p + = max(p, 0). Assuming for simplicity that the battery is fully efficient, the capacity charge (C α n ) n of the BES, valued in [0, C max ], evolves according to the dynamics

The imbalance process defined by

represents how well we are doing for satisfying electricity supply: the ideal situation occurs when S α n = 0, i.e., perfect balance between demand and generation. When S α n > 0, this means that demand is not satisfied, i.e., there is missing power in the microgrid, and when S α n < 0, there is an excess of electricity. In order to ensure that there is no missing power, we impose the following constraint on the admissible control:

but penalize the excess of electricity when S α n < 0 with a proportional cost Q -> 0. We model the residual demand as a mean-reverting process:

where (ε n ) n are i.i.d., R ∈ R, and < 1. The goal of the microgrid manager is to find the optimal (admissible) decision α that minimizes the functional cost

where (.) is the cost function for fuel consumption: (0) = 0, and e.g. (a) = Ka γ , with K > 0, γ > 0. This stochastic control problem fits into the 3-dimensional framework of Section 1 (see also Remark 2.4) with control α valued in


G Boudebs

V Besse

C Cassagne

H Leblond

F Sanchez

Why optical nonlinear characterization using imaging technique is a better choice?

Keywords: nonlinear optics, Z-scan, diffraction, image processing, Fourier optics

We show that the direct measurement of the beam radius in Z-scan experiments using a CCD camera at the output of a 4f-imaging system allows a higher sensitivity and a better accuracy than Baryscan. One of the advantages is to be insensitive to pointing instability of the pulsed laser because no hard aperture is employed as in the usual Z-scan. In addition, the numerical calculations involved here and the measurement of the beam radius are simplified since we do not measure the transmittance through an aperture and it is not subject to mathematical artifacts related to a normalization process, especially when the diffracted light is very low.

INTRODUCTION

The inspection of linear and nonlinear (NL) optical parameters such as refractive index and absorption coefficient is essential for most of the applications in optics (optical limiting, lasers, optical amplifiers, photovoltaic cells, ablation, surgery, all-optical switching...). Measuring the NL refraction and absorption using the far field beam diffraction as in Z-scan method [START_REF] Sheik-Bahae | Sensitive measurement of optical nonlinearities using a single beam[END_REF] is a widely used method nowadays. More recently a variant of this method called Baryscan [START_REF] Godin | Baryscan: a sensitive and user-friendly alternative to Z-scan for weak nonlinearities measurements[END_REF] reporting an increased sensitivity has been published. The authors check the displacement of the beam centroid caused by the focusing (or defocusing) effect induced by the nonlinearity and achieved by a razor blade positioned close to a Position Sensitive Detector (PSD). The razor blade truncates half of the probe beam indicating an increase (or a decrease) of the beam size. Then the PSD shows a change in the position of the beam barycenter. The authors of Ref [START_REF] Godin | Baryscan: a sensitive and user-friendly alternative to Z-scan for weak nonlinearities measurements[END_REF] claim sensitivity up to 50 000 λ using stable CW solidstate laser input. Recently van Stryland's group reported a new dual-arm Z-scan technique [START_REF] Ferdinandus | Dual-arm Z-scan technique to extract dilute solute nonlinearities from solution measurements[END_REF] improving the sensitivity up to 1 000 λ using pulsed lasers. This new method allowed the determination of NL refractive signals considerably smaller than usually found in the single-arm Z-scan method. These studies show that the sensitivity is still relevant and an active research field especially when it was understood that the normalization in the Z-scan method plays an important role. This problem in eclipsing experiments and its consequence in the normalisation process was pointed out first in Ref. [START_REF] Fedus | Sensitivity of the 4f coherent imaging system used in degenerate multi-wave mixing experiments[END_REF]. Therefore, the consideration of the diffraction efficiency allowed to compare the sensitivity of different measurement techniques inside a Z-scan 4f-system [START_REF] Fedus | Experimental techniques using 4f coherent imaging system for measuring nonlinear refraction[END_REF] and had the precious advantage to avoid a "division by zero" with non existing energy. The advantage of using the diffraction efficiency allowed us to increase both the sensitivity of the measurement and the signal to noise ratio by matching properly the object with the field-stop. This optimization is not subject to artefacts linked to a normalization process when dividing by very low intensity profiles. It was shown that the Z-scan is the most sensitive techniques with the highest signal-to-noise ratio allowing measurements of NL induced phase distortions as small as 1 000 λ . But still remain open questions to fully understand the physical phenomenon improving the sensitivity as claimed in [START_REF] Godin | Baryscan: a sensitive and user-friendly alternative to Z-scan for weak nonlinearities measurements[END_REF]. The use of a two-dimensional CCD detector is a more accurate and easier method [START_REF] Chen | Two-dimensional Z scan for arbitrary beam shape and sample thickness[END_REF] and accordingly, Tsigaridas et al presented a modification of Z-scan technique by directly measuring the beam radius as a function of the peak irradiance [START_REF] Tsigaridas | Z-scan technique through beam radius measurements[END_REF].

We will show here that the use of hard-physical apertures (as razor blade) in the Baryscan technique under real pulsed laser with pointing fluctuations could be a severe limitation leading to both lower sensitivity and signal-to-noise ratio when compared to the direct measurement of the output beam waist variation using a CCD sensor.

THE THEORETICAL MODEL

The NL image formation inside the 4f coherent imaging system (Fig. 1) can be described using a simple model based on Fourier optics (see for example Refs. [START_REF] Boudebs | Third-order susceptibility measurements by nonlinear image processing[END_REF], [START_REF] Boudebs | Contrast increasing by third-order nonlinear image processing: a numerical study for microscopic rectangular objects[END_REF]). In this section we recall this approach. The general scheme of beam propagation inside the 4f imaging system is described in details in [START_REF] Fedus | Experimental techniques using 4f coherent imaging system for measuring nonlinear refraction[END_REF].

Let us suppose that the amplitude of electric field distribution at the object plane is Gaussian, ( )

2 2 2 0 e E(x, y) E exp x y /   = - + ω 
 , where x, y are the spatial coordinates, 0 E denotes the on-axis amplitude and ωe is the beam waist at the entry of the setup. Let S(u,v) be the spatial spectrum of E(x,y) :

( ) ( ) S(u,v) E x,y E x,y exp[ j2 (ux vy)]dxdy +∞ +∞ -∞ -∞ = = -π +     ∫ ∫ ɶ F
, where ɶ F denotes the Fourier transform operation, u and v are the normalized spatial frequencies.

The spectrum of the object is propagated over a distance z ' , by taking into account the transfer function of the wave propagation phenomenon ( ) ( ) ( )

2 2 H u, v exp j2 z ' 1 u v   = π -λ -λ λ    
, where λ is the wavelength.

The field amplitude at z ' after free propagation is obtained by computing the inverse Fourier transform:

[ ] 1 E(x,y,z') S(u,v)H(u,v) - = ɶ F .
To calculate the output beam after passing through a lens with focal length f, we apply the phase transformation related to their thickness variations in paraxial approximation:

( ) ( ) 2 2 L t x, y exp j x y f   = -π + λ   .
The first propagation is performed on a distance z ' f = from the object plane to the lens L1. Then we propagate the beam up to the sample located at z using z ' f z = + in H the optical transfer function ( z 0 = at the focus of the lens 1 L ). After that the NL response of the material is taken into account. We assume cubic nonlinearity and a thin NL medium of thickness L exhibiting (i) linear absorption defined by α (m -1 ), (ii) two-photon absorption defined by β (m/W) and (iii) NL refraction defined by 2 n (m 2 /W). In these conditions the transmittance of the sample is described as:

( ) ( ) ( ) ( ) { } ( ) NL 1 2 L L eff S u, v, z T u, v, z e 1 q u, v, z exp j u, v, z S u, v, z - α   = = + ∆φ       , (1)
where ( )

( ) eff q(u, v, z) L I u, v, z = β with ( ) L eff
NL eff 2 eff eff (u, v, z) 2 n L I u, v, z ∆φ = π λ where ( ) ( ) ( ) eff I u, v, z I(u, v, z) log 1 q u, v, z q u, v, z = +    
is the effective intensity seen by the material for Gaussian beams. Let 0 eff 0 q L I = β and NL0 NL eff eff (0, 0, 0) ϕ = ∆φ denote the on-axis NL absorption and NL phase shift at the focus, respectively, I0 is the focal on-axis intensity.

Next, the measurement of the beam waist can be performed using the ISO standard definition. Based on the second moment of ( ) I x, y the D4σ method gives 4 times the standard deviation of the intensity distribution. For example, the beam radius in the x direction is:

( )( ) ( ) 2 x I x, y x x dxdy 2 I x, y dxdy +∞ +∞ -∞ -∞ +∞ +∞ -∞ -∞ - ω = ∫ ∫ ∫ ∫
, where ( )

( ) I x, y xdxdy x I x, y dxdy +∞ +∞ -∞ -∞ +∞ +∞ -∞ -∞ = ∫ ∫ ∫ ∫ (2)
is the centroid of the beam profile in the x-direction. For Gaussian beams, the D4σ method gives the same result as the 2 1 e method, whereas for other beam shapes there can be significant deviations. The 2 1 e width measurements are noisier than D4σ width measurements depending on the integral of all the pixels contained in each frame. Let us add that for multimodal marginal distributions (a beam profile with multiple peaks), the D4σ width would be a better choice.

To calculate energy inside the aperture in the Z-scan configuration we must integrate over a radius

( ) { } 1/ 2 a i r 0.5log 1/ 1 S = ω -   
 , where S is the closed aperture linear transmittance and i ω is the beam waist in the 0.34 ∆ω = × ϕ ). Note that this formula allows to measure n2 when the other parameters in eff NL0 ϕ are known. Moreover the linearity remains valid in presence of relatively high NL absorption (here 0 q 0.58 = ).

EXPERIMENTAL RESULTS

The detailed description of our experimental setup is given in [START_REF] Fedus | Degenerate multi-wave mixing inside a 4f imaging system in presence of nonlinear absorption[END_REF]. Excitation is provided by a Nd:YAG laser delivering linealy polarized 17 ps single pulses at 1064 λ = µm with 10 Hz repetition rate. The input intensity is varied using a half-wave plate and a Glan prism. A beam splitter at the entry of the setup (Fig. 1) allows to monitor any fluctuation occurring in the incident laser beam. The focal lengths of lens L1 and L2 are equal to 20 cm. To detect the light in the image plane we use a 1000 ×1018 pixels cooled (-30°C) CCD camera (12x12µm 2 ). Open and closed aperture signals are numerically processed from the acquired images. The CCD camera allows performing Z-scan, Baryscan and D4σ profiles out from the same acquired images at each z-position by numerically changing the soft aperture. Two sets of acquisitions are performed in the linear and the NL regimes in order to correct for the diffraction effects due to sample inhomogeneities and imperfections by subtracting the data related to the scans in each regime. When using the D4σ method one have to calculate the centroid of the beam before processing the beam waist using eq. 2. This processing inside the images is equivalent to follow the pointing instability of the pulsed laser and thus reduces the noise that could appear with hard apertures as the ones used with Z-scan and Baryscan. Fig. 3 (a) shows the experimental results for As2Se3, a chalcogenide glass with relatively high NL absorption. The D4σ results (filled squares) can be compared with those of Baryscan (circles) using pulsed laser beams. Notice that the useful Baryscan signal is totally buried inside the noise due to laser pointing instability while the beam waist measurement method gives a cleaner signal with a higher signal to noise ratio. The same behaviour could be found for Z-scan using hard aperture. The calculated profiles (solid line for D4σ and dashed line for Baryscan) are shown in the same figure. The agreement is very good. The NL parameters obtained were: 3b). Other parameters are: 0 I 4.4 = GW/cm 2 (absolute measurement obtained as in [START_REF] Boudebs | Absolute measurement of the nonlinear refractive indices of reference materials[END_REF]), α = 0, 1.064 µm λ = , i 1.9 ω = mm and L 0.64 mm = . Moreover the D4σ method is unambiguously two times more sensitive than the Baryscan method. Physically we consider the totality of the pixels in order to measure the beam waist. So the quantity of information is twice that obtained with the razor blade in Baryscan. This is illustrated in Fig. 3 (a) by the difference between the peaks and the valleys inside the two scans showing approximately a factor 2 in the sensitivity. Fig. °3 (b) shows the comparison of the beam waist relative variation (squares) and the usual normalized Z-scan transmittance (stars) which is the result of the division between the closed and the open aperture Z-scan traces shown in the inset. The radius of the soft circular closed aperture is calculated for S = 0.73 to optimize the signal to noise ratio [START_REF] Fedus | Experimental techniques using 4f coherent imaging system for measuring nonlinear refraction[END_REF]. The NL absorption profile is obtained by integrating over all the pixels. To obtain a good comparison the centroid method is used for both Z-scan profiles in order to monitor the beam transversally for pointing laser fluctuations (Fig. 3b). The calculations have been made using the same acquired images showing unambiguously approximately the same sensitivity. However, the advantage in the D4σ method is that there is no need to divide two different Z-scan traces. The direct measurement of the beam waist gives directly the signal of the NL refraction.

( ) -
Generally, the hard aperture at the output of an experimental setup intercepts the laser profile at different position due to the laser pointing instability. This leads to a variation of the transmitted intensity showing a fluctuating noise. The use of the CCD camera allows to follow the centroid of the image. Then, it is possible to apply a numerical soft apertures of different shapes. For arbitrary (possibly not Gaussian) beam, it is recommended to use the D4σ or second moment beam width which is the ISO standard definition [START_REF]ISO Standard 11146[END_REF]. To improve the signal-to-noise ratio, the centroid of the beam must be calculated for the Z-scan traces as well as in the Baryscan method.

L 1 e 2 I

 12 -α = -α representing the effective length and I(u, v, z) denotes the intensity of the laser beam within the sample. The latter quantity is related to the input electric field through the expression (u, v, z) S(u, v, z) ∝ . Here NL eff (u, v, z) ∆φ can be considered as the effective NL phase shift:

Figure 1 :

 1 Figure 1: The 4f imaging system. The nonlinear material (NLM) is moved around the focal plane. The labels refer to: lenses (L1-L2). The object at the entry is a Gaussian beam, while a circular aperture or a razor blade can be considered to as soft apertures spatially filtering the image to obtain Z-scan or Baryscan configurations.

Figure 2 :

 2 Figure 2: (a) Relative variation of the beam waist versus z. (b) Calculated Δωpv, versus the effective phase shift at the focus for Gaussian input beam showing a linear dependence. z(mm)

Figure 3 )

 3 Figure 3) (a) Comparison between D4σ measurements method (filled squares) with that of Baryscan (circles). The solid and dashed lines show the numerical simulations. (b): Comparison of the beam waist relative variation (squares; vertically shifted to 1) and the usual normalized Z-scan transmittance (stars) (division of the closed and the open aperture Z-scan traces) for highly NLA material. The inset shows the open aperture Z-scan transmittance. The processing has been made using the same acquired images.

CONCLUSION

The D4σ method consisting of measuring the beam radius at the output of a 4f-Z-scan method using a CCD camera provides a signal having a peak-valley shape. One of the advantages of this technique is to be insensitive to pointing instability of the pulsed laser because no hard aperture is employed as in the usual Z-scan or Baryscan methods. Numerical calculations allow to obtain simple relations that can be used for the measurements simplifying the procedure when a NL absorbing material is used. Moreover the normalization process does not lead to any ambiguity especially when the intensity of diffracted light is very low.