Deep Learning Methods for Mean Field Control Problems with Delay

We consider a general class of mean field control problems described by stochastic delayed differential equations of McKean-Vlasov type. Two numerical algorithms are provided based on deep learning techniques, one is to directly parameterize the optimal control using neural networks, the other is based on numerically solving the McKean-Vlasov forward anticipated backward stochastic differential equation (MV-FABSDE) system. In addition, we establish the necessary and sufficient stochastic maximum principle of this class of mean field control problems with delay based on the differential calculus on function of measures, and the existence and uniqueness results are proved for the associated MV-FABSDE system under suitable conditions.

Both mean field control and mean field games are used to characterize the asymptotic behavior of a stochastic game as the number of players grows to infinity under the assumption that all the agents behave similarly, but with different notation of equilibrium. The mean field games consist of solving a standard control problem, where the flow of measures is fixed, and solving a fixed point problem such that this flow of measures matches the distribution of the dynamic of a representative agent. Whereas, the mean field control problem is a nonstandard control problem in the sense that the law of state is present in the McKean-Vlasov dynamic, and optimization is performed while imposing the constraint of distribution of the state. More details can be found in [3] and [1].
In this paper, we considered a general class of mean field control problem with delay effect in the McKean-Vlasov dynamic. We derived the adjoint process associated with the McKean-Vlasov stochastic delayed differential equation, which is an anticipated backward stochastic differential equation of McKean-Vlasov type due to the fact that the conditional expectation of the future of adjoint process as well as the distribution of the state dynamic are involved. This type of anticipated backward stochastic differential equations (BSDE) was introduced in [15], and for the general theory of BSDE, we refer [18]. The necessary and sufficient part of stochastic maximum principle for control problem with delay in state and control can be found in [7]. Here, we also establish a necessary and sufficient stochastic maximum principle based on differential calculus on functions of measures as we consider the delay in the distribution. In the meantime, we also prove the existence and uniqueness of the system of McKean-Vlasov forward anticipated backward stochastic differential equations (MV-FABSDE) under some suitable conditions using the method of continuation, which can be found in [18], [14], [2] and [5]. For a comprehensive study of FBSDE theory, we refer to [13].
When there was no delay effect in the dynamic, [12] proved the relation between the solution to the FBSDE and quasi-linear partial differential equation (PDE) via "Four Step Scheme". [8] and [16] explored the use of deep learning for solving these PDEs in high dimensions. However, the class of fully coupled MV-FABSDE considered in our paper has no explicit solution. Here, we present one algorithm to tackle the above problem by means of deep learning techniques. Due to the non-Markovian nature of the state dynamic, we apply the long short-term memory (LSTM) network, which is able to capture the arbitrary long-term dependencies in the data sequence. It also partially solves the vanishing gradient problem in vanilla recurrent neural networks (RNNs), as was shown in [11]. The idea of our algorithm is to approximate the solution to the adjoint process and the conditional expectations of the adjoint process. The optimal control was readily obtained after the MV-FABSDE being solved. We may also emphasis that the way that we present here for numerically compute conditional expectation may have a wide range of applications, and it is simple to implement. We also present another algorithm solving the mean field control problem by directly parameterizing the optimal control. Similar idea can be found in the policy gradient method in the regime of reinforcement learning [17] as well as in [10]. Numerically, the two algorithms that we propose in this paper yield the same results. Besides, our approaches are benchmarked to the case with no delay for which we have explicit solutions.. The paper is organized as follows. We start with a N player game with delay, and let number of players goes to infinity to introduce a mean field control problem in Section 2. Next, in Section 3, we mathematically formulate the feedforward neural networks and LSTM networks, and we propose two algorithms to numerically solve the mean field control problem with delay using deep learning techniques. This is illustrated on a simple linear-quadratic toy model, however with delay in the control. One algorithm is based on directly parameterizing the control, and the other depends on numerically solving the MV-FABSDE system. In addition, we also provide an example of solving a linear quadratic mean field control problem with no delay both analytically, and numerically. The adjoint process associated with the delayed dynamic is derived, as well as the stochastic maximum principle is proved in Section 4. Finally, the uniqueness and existence solution for this class of MV-FABSDE are proved under suitable assumptions via continuation method in Section 5.

Formulation of the Problem
We consider a N -player game with delay in both state and control. The dynamic (X i t ) 0≤t≤T for player i ∈ {1, . . . , N } is given by a stochastic delayed differential equation (SDDE), for T > 0, τ > 0 be given constants, where X t = (X 1 t , · · · , X N t ), and where ((W i t ) t∈[0,T ] ) i=1,···,N are N independent Brownian motions defined on the space (Ω, F, P), (F t ) 0≤t≤T being the natural filtration of Brownian motions.
are progressively measurable functions with values in R. We denote A a closed convex subset of R, the set of actions that player i can take, and denote A the set of admissible control processes. For each i ∈ {1, . . . , N }, A-valued measurable processes (α i t ) 0≤t≤T satisfy an integrability condition such that E T −τ |α i t | 2 dt < +∞. Given an initial condition x 0 = (x 1 0 , · · · , x N 0 ) ∈ R N , each player would like to minimize his objective functional: In order to study the mean-field limit of (X t ) t∈[0,T ] , we assume that the system (2.1) satisfy a symmetric property, that is to say, for each player i, the other players are indistinguishable. Therefore, drift b i and volatility σ i in (2.1) take the form of and the running cost f i and terminal cost g i are of the form where we use the notation µ N t for the empirical distribution of X = (X 1 , · · · , X N ) at time t, which is defined as Next, we let the number of players N goes to +∞ before we perform the optimization. According to symmetry property and the theory of propagation of chaos, the joint distribution of the N dimensional process (X t ) 0≤t≤T = (X 1 t , . . . , X N t ) 0≤t≤T converges to a product distribution, and the distribution of each single marginal process converges to the distribution of (X t ) 0≤t≤T of the following Mckean-Vlasov stochastic delayed differential equation (MV-SDDE). For more detail on the argument without delay, we refer [3] and [4].
3) We then optimize after taking the limit. The objective for each player of (2.2) now becomes where we denote µ t := L(X t ) the law of X t .

Solving Mean-Field Control Problems Using Deep Learning Techniques
Due to the non-Markovian structure, the above mean-field optimal control problem (2.3)-(2.4) is difficult to solve either analytically or numerically. Here we propose two algorithms together with four approaches to tackle the above problem based on deep learning techniques. We would like to use two types of neural networks, one is called the feedforward neural network, and the other one is called Long Short-Term Memory (LSTM) network. For a feedforward neural network, we first define the set of layers M ρ d,h , for x ∈ R d , as d is called input dimension, h is known as the number of hidden neurons, A ∈ R h×d is the weight matrix, b ∈ R h is the bias vector, and ρ is called the activation function. The following activation functions will be used in this paper, for some x ∈ R, Then feedforward neural network is defined as as a composition of layers, so that the set of feedforward neural network with l hidden layers we used in this paper is defined as The LSTM network is one of recurrent neural network(RNN) architectures, which are powerful for capturing long-range dependence of the data. It is proposed in [11], and it is designed to solve the shrinking gradient effects which basic RNN often suffers from. The LSTM network is a chain of cells. Each LSTM cell composes of a cell state, which contains information, and three gates, which regulate the flow of information. Mathematically, the rule inside tth cell follows, where the operator denotes the Hadamard product. (Γ ft , Γ it , Γ ot ) ∈ R h × R h × R h represents forget gate, input gate and output gate respectively, h refers the number of hidden neurons. x t ∈ R d is the input vector with d features. a t ∈ R h is known as the output vector with initial value a 0 = 0, and c t ∈ R h is known as the cell state with initial value c 0 = 0. A · ∈ R h×d are the weight matrices connecting input and hidden layers, U · ∈ R h×h are the weight matrix connecting hidden and output layers, and b ∈ R h represents bias vector. The weight matrices and bias vectors are shared through all time steps, and are going to be learned during training process by back-propagation through time (BPTT, which can be implemented in Tensorflow platform. Here we define the set of LSTM network up to time t as In particular, we specify the model in a linear-quadratic form, which is inspired by [5] and [9]. The objective function is defined as where σ, c f , c t > 0 are given constants, and m t := R xdµ t (x) denotes the mean of X at time t. In the following subsections, we solve the above problem numerically using two algorithms together with four approaches. The first two approaches are to directly approximate the control by either a LSTM network or a feedforward neural network, and minimize the objective (3.5) using stochastic gradient descent algorithm. The third and fourth approaches are to introduce the adjoint process associated with (3.6), and approximate the adjoint process, conditional expectation of adjoint process using neural networks.

Approximating the Optimal Control Using Neural Networks
We first set ∆t = T /N = τ /D for some positive integer N . The time discretization becomes The discretizing the SDDE (3.6) according to EulerMaruyama scheme now reads where (∆W t i ) 0≤i≤N −1 are independent, normal distributed sequence of random variables with mean 0 and variance 1. First, from the definition of open loop control, and due to non-Markovian feature of (3.6), the open-loop optimal control is a function of the path of the Brownian motions up to time t, i.e., α(t, (W s ) 0≤s≤t ). We are able to describe this dependency by a LSTM network by parametrizing the control as a function of current time and the discretized increments of Brownian motion path, i.e., 8) for some h 1 ∈ Z + . We remark that the last dense layer is used to match the desired output dimension.
The second approach is again directly approximate the control but with a feedforward neural network. Due to the special structure of our model, where the mean of dynamic in (3.6) is constant, the mean field control problem coincides with the mean field game problem. In [9], authors solved the associated mean field game problem using infinite dimensional PDE approach, and found that the optimal control is a function of current state and the past of control. Therefore, the feedforward neural network with l layers, which we use to approximate the optimal control, is defined as From Monte Carlo algorithm, and trapezoidal rule, the objective function (3.5) now becomes where M denotes the number of realizations andX := 1 M M j=1 X (j) denotes the sample mean. After plugging in the neural network either given by (3.8) or (3.9), the optimization problem becomes to find the best set of parameters either (Φ 1 , Ψ 1 ) or Ψ 2 such that the objective J(Φ 1 , Ψ 1 ) or J(Ψ 2 ) is minimized with respect to those parameters.
The algorithm works as follows: Algorithm 1: Algorithms for solving mean field control problem with delay by directly approximating the optimal control using neural networks Initialization of parameters Θ 0 = (Φ 1 , Ψ 1 ) for approach 1 (3.8) for some network ϕ given by (3.8) or by (3.9) at t 0 with proper inputs;    On the left, we compare one representative optimal trajectory of (X t i ) t 0 ≤t i ≤t N . The plot on the right show the comparison of one representative optimal trajectory of (α t i ) t 0 ≤t i ≤t N between approach 1 and approach 2. On the left, we compare the sample mean of optimal trajectories of (X t i ) t 0 ≤t i ≤t N . The plot on the right show the comparison of sample mean of trajectories of optimal control (α t i ) t 0 ≤t i ≤t N between approach 1 and approach 2.

Approximating the Adjoint Process Using Neural Networks
The third and fourth approaches are based on numerically solving the MV-FABSDE system using LSTM network and feedforward neural networks. From Section 4, we derive the adjoint process, and prove the sufficient and necessary parts of stochastic maximum principle. From (4.7), we are able to write the backward stochastic differential equation associated to (3.6) as, with terminal condition Y T = c t (X T − m T ), and Y s = 0 for s ∈ (T, T + τ ]. The optimal control (α t ) 0≤t≤T can be obtained in terms of the adjoint process Y t from the maximum principle, and it is given by From the Euler-Maruyama scheme, the discretized version of (3.6) and (3.11) now reads, where we use the sample averageX t = 1 t to approximate the expectation of X t . In order to solve the the above MV-FABSDE system, we need to approximate The third approach consists of approximating , Z t i ) 0≤t i ≤t N using three LSTM networks as functions of current time and the discretized path of Brownian motions respectively, i.e., Again, the last dense layers are used to match the desired output dimension.
Since approach 3 consist three neural networks with large number of parameters, which is hard to train in general, we would like to make the following simplification in approach 4 for approximating (3.15) In words, the algorithm works as follows. We first initialize the parameters given by either in (3.14) or in (3.15), we update X t i+1 and Y t i+1 according to (3.12), and the solution to the backward equation at t i+1 is denoted byỸ t i+1 . In the mean time, Y t i+1 is also approximated by a neural network. In such case, we referỸ · as the label, and Y · given by the neural network as the prediction. We would like to minimize the mean square error between these two. At time T , Y t N is also supposed to match c(X t N −X t N ), from the terminal condition of (3.11). In addition, the conditional expectation given by a neural network should be the best predictor ofỸ t i+D , which implies that we would like to find the set of parameters Θ EY such that E[(Ỹ t i+D − ϕ EY t i (Θ EY )] is minimized for all t i ∈ {t 0 , . . . , t N −D }. Therefore, for M samples, we would like to minimize two objective functions L 1 and L 2 defined as (3.16) The algorithm works as follows: Algorithm 2: Algorithms for solving mean field control problem with delay according to MV-FABSDE Initialization of parameters (Θ Y , Θ EY , Θ Z ) for approach 3 as in (3.14) or approach 4 as in (3.15); given by either (3.14) or by (3.15) ; α given by (3.14) or by(3.15) at t i+1 ; given by (3.14) or by (3.15) at t i+1 ; • Compute the gradient ∇L 1 (Θ Y , Θ Z )by backpropagation through time; • Update Θ Y e+1 and Θ Z e+1 according to SGD algorithm; • Compute the gradient ∇L 2 (Θ EY )by backpropagation through time; • Update Θ EY e+1 according to SGD algorithm; end Again, in the following graphics, we choose x 0 = 0, c f = c t = 1, σ = 1, T = 10, τ = 4, ∆t = 0.1, M = 4000. For a specific representative path, the underlying Brownian motion paths are the same for different approaches. Figure 3.3 compares one representative optimal trajectory of the state dynamic and the control via two approaches, and they coincide. Figure 3.4 plots the sample average of optimal trajectory of the dynamic and the control, which are trajectories of 0, which is the same as the theoretical mean. Comparing to Figure 3.1 and Figure 3.2, as well as based on numerous experiments, we find that given a path of Brownian motion, the two algorithms would yield the same optimal trajectory of state dynamic and the same path for the optimal control. From Figure 3.6, the loss L 1 as defined in (3.16) are minimized to 0 for both approach 3 and approach 4. This can also be observed from Figure 3.5, since the red dash line and the blue solid line coincide for both left and right graphs. In addition, from the righthand side of Figure 3.6, we observe the loss L 2 as defined in (3.16) converges. This is due to the fact that the conditional expectation can be understand as an orthogonal projection. Figure 3.7 plots 64 sample paths of the process (Z t i ) t 0 ≤t i ≤t N , which seems to be a deterministic function since σ is constant in this example. Finally, Figure 3.8 shows the convergence of the value function as number of epoches increases. Although, two algorithms would arrive approximately at the same optimal value with is around 6. This confirms that the out control problem has a unique solution. In section 5, we show that the MV-FASBDE system is uniquely solvable. It is also observable that the first algorithm converges faster than the second one, since it directly paramerizes the control using one neural network, instead of solving the MV-FABSDE system, which uses three neural networks. On the left, we compare one representative optimal trajectory of (X t i ) t 0 ≤t i ≤t N . The plot on the right shows the comparison of one representative optimal trajectory of (α t i ) t 0 ≤t i ≤t N between approach 3 and approach 4.
The plot on the right shows the comparison of sample mean of trajectories of optimal control (α t i ) t 0 ≤t i ≤t N between approach 3 and approach 4.
from approach 3 (one the left) and from approach 4 (on the right).

Numerically Solving the Optimal Control Problem with No Delay
Since the algorithms we proposed embrace the case with no delay, we illustrate the comparison between numerical results and the analytical results. By letting τ > T we obtain α t−τ = 0 in (3.6), and we aim at solving the following linear-quadratic mean-field control problem by minimizing  subject to dX t =α t dt + σdW t , t ∈ [0, T ], Again, from Section 4, we find the optimal control where (Y t , Z t ) is the solultion of the following adjoint process, Next, we make the ansatz for some deterministic function φ t , satisfying the terminal condition φ T = c t . Differentiating the ansatz, the backward equation should satisfy whereφ t denotes the time derivative of φ t . Comparing with (3.19), and identifying the drift and volatility term, φ t must satisfy the scalar Riccati equation, Numerically, we apply the two deep learning algorithms proposed in the previous section. The first algorithm directly approximates the control. According to the open loop formulation, we set 25) for some h ∈ Z + . We remark that the last dense layer is used to match the desired output dimension. The second algorithm numerically solves the forward backward system as in (3.18) and (3.19). From the ansatz (3.20) and the Markovian feature, we approximate (Y t , Z t ) 0≤t≤T using two feedforward neural networks, i.e., (3.26) Figure 3.9 shows the representative optimal trajectory of (X t i −X t i ) t 0 ≤t i ≤t N and (α t i ) t 0 ≤t i ≤t N from both algorithms, which are exactly the same. The symmetry feature can be seen from the computation (3.24). On the left of Figure 3.10 confirms the mean of the processes (X t i −X t i ) t 0 ≤t i ≤t N and (α t i ) t 0 ≤t i ≤t N are 0 from both algorithms. The right picture of Figure 3.10 plots the data points of x against α, and we can observe that the optimal control α is linear in x as a result of (3.24), and the slope tends to be -1, since φ = 1 solves the scalar Riccati equation (3.22). Finally, Figure 3.11 plots representative optimal trajectories of the solution to the adjoint equations (Y t , Z t ). On the left, we observe that the ajoint process (Y t ) 0≤t≤T matches the terminal condtion, and on the right, (Z t ) 0≤≤T appears to be a deterministic process of value 1, and this matches the result we compute previously. Figure 3.9: Representative optimal trajectory of(X t i −X t i ) t 0 ≤t i ≤t N and (α t i ) t 0 ≤t i ≤t N from algorithm 1 (on the left) and from algorithm 2 (on the right).
The picture on the right plots 64 trajectories of (Z t i ) t 0 ≤t i ≤t N .

Stochastic Maximum Principle for Optimality
In this section, we derive the adjoint equation associated to our mean field stochastic control problem (2.3) and (2.4). The necessary and sufficient parts of stochastic maximum principle have been proved for optimality. We assume (H4.1) b, σ are differentiable with respect to (X t , µ t , X t−τ , µ t−τ , α t , α t−τ ); f is differentiable with respect to (X t , µ t , X t−τ , µ t−τ , α); g is differentiable with respect to (X T , µ T ). Their derivatives are bounded.
In order to simplify our notations, let θ t = (X t , µ t , α t ). For 0 < < 1, we denote α the admissible control defined by for any (α) 0≤t≤T and (β) 0≤t≤T ∈ A. X t := X α t is the corresponding controlled process. We define ∇X t := lim →0 X t − X α t to be the variation process, which should follow the following dynamic for t ∈ (0, T ], is a copy of (X t , ∇X t ) defined on (Ω,F,P), where we apply differential calculus on functions of measure, see [3] for detail. ∂ x b, ∂ xτ b, ∂ α b, ∂ ατ b are derivatives of b with respect to (X t , X t−τ , α t , α t−τ ) respectively, and we use the same notation for ∂ · σ.
In the meantime, the Gateaux derivative of functional α → J(α) is given by In order to determine the adjoint backward equation of (Y t , Z t ) 0≤t≤T associated to (2.3), we assume it is of the following form: Next, we apply integration by part to ∇X t and Y t . It yields We integrate from 0 to T , and take expectation to get Using the fact that Y t = Z t = 0 for t ∈ (T, T + τ ], we are able to make a change of time, and by Fubini's theorem, so that (4.4) becomes x, µ, x τ , µ τ , α, α τ , y, z) = b(t, x, µ, x τ , µ τ , α, α τ )y + σ(t, x, µ, x τ , µ τ , α, α τ )z + f (t, x, µ, x τ , µ τ , α) (4.6) Using the terminal condition of Y T , and plug (4.5) into (4.2), and set the integrand containing ∇X t to zero, we are able to obtain the adjoint equation is of the following form (4.7) Theorem 4.1. Let (α t ) 0≤t≤T ∈ A be optimal, (X t ) 0≤t≤T be the associated controlled state, and (Y t , Z t ) 0≤t≤T be the associated adjoint processes defined in (4.7). For any β ∈ A, and t ∈ [0, T ], (4.8) Proof. Given any (β t ) 0≤t≤T ∈ A, we perturbate α t by α t := α t + (β t − α t ) for 0 ≤ ≤ 1. Using the adjoint process (4.7), and apply integration by parts formula to (∇X t Y t ). Then plug the result into (4.2), and the Hamiltonian H is defined in (4.6). Also, since α is optimal, we have (4.9) Now, let C ∈ F t be an arbitrary progressively measurable set, and denote C the complement of C. We choose β t to be β t := β1 C + α t 1 C for any given β ∈ A. Then, (4.11) Remark 4.2. When we further assume that H is convex in (α t , α t−τ ), then for any β, β τ ∈ A in Theorem 4.1, we have as a direct consequence of (4.8).
Proof. Let (α t ) 0≤t≤T ∈ A be a admissible control, and let (X t ) 0≤t≤T = (X α t ) 0≤t≤T be the corresponding controlled state. From the definition of the objective function as in (2.4), we first use convexity of g, and the terminal condition of the adjoint process Y t in (4.7), then use the fact that H is convex, and because of (4.12), we have the following ≤0. (4.13)