We present a neural network approach for approximating the value function of high-dimensional stochastic control problems. Our training process simultaneously updates our value function estimate and identifies the part of the state space likely to be visited by optimal trajectories.
Our approach leverages insights from optimal control theory and the fundamental relation between semi-linear parabolic partial differential equations and forward-backward stochastic differential equations.
To focus the sampling on relevant states during neural network training, we use the stochastic Pontryagin maximum principle (PMP) to obtain the optimal controls for the current value function estimate.
By design, our approach coincides with the method of characteristics for the non-viscous Hamilton-Jacobi-Bellman equation as the dynamics become deterministic.
Our training loss consists of a weighted sum of the objective functional of the control problem and penalty terms that enforce the HJB equations along the sampled trajectories.
Importantly, training is unsupervised in that it does not require solutions of the control problem.

In our numerical experiments, we compare our method to existing solvers for a more general class of semi-linear PDEs.
Using a two-dimensional model problem, we demonstrate the importance of the stochastic PMP to inform the sampling.
For a 100-dimensional benchmark problem, we demonstrate that our approach improves accuracy and time-to-solution.

Hamilton-Jacobi-Bellman equation, high-dimensional stochastic optimal control, forward-backward stochastic differential equations, neural networks, stochastic maximum principle

35F21, 49M99, 68T07

§ INTRODUCTION

We consider Stochastic Optimal Control (SOC) problems that require finding a policy to control randomly perturbed dynamical systems so as to optimize a given objective functional.
Problems of this type arise in many fields including finance, biology, robotics and many other engineering applications; see, e.g., [10, 34] for references and extensive theoretical discussion on SOC problems.

Dynamic programming (DP) is a prominent framework for solving SOC problems.
At its core, DP seeks to find the value function, which assigns every state of the system the optimal cost-to-go.
For many problems, recovering the optimal control from the value function is straightforward.
One way to compute the value function is by solving the Hamilton-Jacobi-Bellman (HJB) equation [2, 37], which is a semi-linear parabolic Partial Differential Equation (PDE).
The two main challenges in solving the HJB equations is the forward-backward structure and the high dimensionality.

The space dimension of the HJB equation equals the state dimension of the dynamical system to be controlled.
Hence, many numerical schemes for solving PDEs cannot be applied to realistic problem instances due to the Curse-of-Dimensionality (CoD); for example, the computational costs of approaches that discretize the state space with a mesh typically grows exponentially in the dimensions of the system.
Hence, the applicability of existing approaches that employ spatial discretizations (for example, [4, 6, 16, 20, 32, 22, 23, 32, 35]) is limited to state dimensions $d \leq 3$. Mesh-free methods based on radial basis function methods (for example, [36]) can be effective for slightly larger $d$ but ultimately also suffer from CoD for $d = \mathcal{O}(100)$, as we consider here.

A first requirement for mitigating the CoD in SOC problems is to alleviate the need for a spatial discretization.
This can be achieved by using a nonlinear version of the Feynman-Kac lemma and replacing the HJB equation by a system of Forward-Backward Stochastic Differential Equations (FBSDEs); see, for example, [5, 28, 37].
Discretizing the FBSDE system using the classical Euler Maruyama scheme [19] yields mesh-free numerical schemes that have gained importance for SOC problems and non-linear second-order PDEs; see, for example, [3, 8, 9, 13, 25] and references therein.

A second requirement for mitigating the CoD is to parameterize the value function effectively in high-dimensions.
Due to their universal approximation property and many advances in deep learning, the use of neural networks (NNs) as function approximators has been gathering significant attention recently.
The idea of combining the BSDE approach with NNs has been pioneered by the seminal works [7, 12] that use neural networks to solve high-dimensional semi-linear parabolic PDEs and include an example for HJB equations as a special case.
Around the same time,  [33] achieves similarly promising results for neural-network based solution of the HJB equation.
These works developed learning algorithms for approximating the value function or its gradient using the FBSDE system.
Their success has sparked several extensions that improve the NN architectures and consider other types of PDEs; see, for example, [15, 17, 29].

Even after combining FBSDEs and deep learning, solving the HJB equation globally remains affected by the curse of dimensionality as it would require sampling the entire state-space.
Fortunately, the optimal state trajectories of many SOC problems do not cover the entire space. This motivates us to develop semi-global approaches, that is, we seek to reliably estimate the value function in those parts of the state space that are likely to be visited when following an optimal policy.
Clearly, these tasks are interrelated: Finding the relevant parts of the state space depends on the value function and, vice versa, the value function needs to be trained using samples from the state space.

The key idea of our approach is to use the stochastic Pontryagin maximum principle (PMP) [31] to define the forward system in the FBSDE approach in terms of the value function.
The PMP provides a set of necessary optimality conditions, known as the Hamiltonian system, and a feedback form that allows recovering the optimal control from the value function; see, for example, [37].
As we discuss in more detail later, the choice of the forward SDE is crucial as it is used to sample the parts of the state space along which the value function estimate is improved using a loss that includes the backward SDE.
It is worth noting that the FBSDE framework allows choosing the forward system almost arbitrarily.
Thus, our choice is theoretically on par with the standard Brownian motion used for exploring the state space in [12, 33], but as our numerical experiments suggest enable solving a wider class of practical problems; see <ref>.

A theoretical advantage of our proposed scheme over existing approaches is its consistency with deterministic optimal control problems.
In the absence of uncertainty in the dynamics, our approach reduces to the method of characteristics for the HJB equation.
Hence our approach can also be seen as an extension of the neural network approaches for deterministic optimal control proposed in [27, 21].

The rest of the paper is organised as follows:
In <ref> we introduce the parts of SOC theory used to obtain our proposed forward SDE.
In <ref>, we provide the NN approach to learn the value function using the FBSDE reformulation from <ref>. Using an intuitive 2D example, we show the importance of sampling.  To illustrate the potential of our method, we consider the 100-dimensional benchmark problem also used in [12, 33], propose a modification that leads to a more challenging sampling problem, and investigate the scaling of our approach with respect to problem dimension.
Finally we conclude the paper and discuss future directions.

§ STOCHASTIC OPTIMAL CONTROL BACKGROUND

In this section, we describe the Stochastic Optimal Control (SOC) problems considered in this work and review the key results from SOC theory that motivates our approach; our discussion follows [37] and we refer to this text book for a more comprehensive background. We first introduce the SOC problem and then review its underlying theory; to be specific, we review the stochastic Pontryagin Maximum Principle (PMP), the Hamilton-Jacobi-Bellman (HJB) equation, and its reformulation into a system of forward backward stochastic differential equations (FBSDEs) obtained from a nonlinear version of the Feynman-Kac formula.

§.§ Stochastic Optimal Control Problem

Let $(\Omega, \Fc, \bbF = \{\Fc_t\}_{t\ge 0}, \bbP)$ be a given complete probability space, $W(s)$ be a $d$-dimensional Brownian motion on $(\Omega, \Fc, \bbF, \bbP)$ where $s$ denotes the time.
For a fixed initial state $\bfx$ at some time $0<t<T<\infty$, we seek to control the randomly perturbed dynamical system
\begin{equation}\label{eq:charOC_stoc}
\begin{cases}
\du \bfz_{t,\bfx}(s) = f(s,\bfz_{t,\bfx}(s), \bfu_{t,\bfx}(s, \bfz_{t,\bfx}(s)))\du s + \sigma(s, \bfz_{t,\bfx}(s)) \du W(s), \quad t \leq  s \leq  T, \\
\bfz_{t,\bfx}(t)= \bfx.
\end{cases}
\end{equation}

Here, $\bfz_{t,\bfx} : [t,T]\ \to \R^d$ describes the state and $\bfu_{t,\bfx}: [t,T]\times  \R^d \to U$ describes the control of the system,
the function $\sigma: [t,T]\times \R^d \to \R^{d\times d}$ represents the diffusion coefficient, and $f:[t,T]\times \R^d \times U\to \R^d$ represents the drift of the system.
We assume that the set of admissible controls $U\subset \R^k$ is closed. We seek to minimize the objective functional

\begin{equation}\label{eq:Joc_stoc}
J_{t,\bfx}(\bfu_{t,\bfx}) = \bbE \left\{G\big(\bfz_{t,\bfx}(T)\big) + \int_{t}^T L\big(s,\bfz_{t,\bfx}(s), \bfu_{t,\bfx}(s, \bfz_{t,\bfx}(s))\big) \, \du s\right\},
\end{equation}

which is comprised of the running cost $L:[t,T]\times \R^d \times U\to \R$ and the terminal cost $G:\R^d\to \R$.
Here, the expectation is taken with respect to perturbances of the dynamics <ref> that is described by the probability measure $\bbP$. We assume sufficient regularity conditions on $f$, $\sigma$, $G$ and $L$, see <cit.> for a list of assumptions.

The value function assigns the optimal cost-to-go to any initial state, that is,

\begin{equation}\label{eq:OC_stoc}
\Phi(t,\bfx)=\inf_{\bfu_{t,\bfx}}  J_{t,\bfx}(\bfu_{t,\bfx}),
\end{equation}

and a solution $\bfu_{t,\bfx}^*$ to <ref> incurring this minimum value is called an optimal control.

The generalized Hamiltonian, $H :[t,T] \times \R^d \times \R^d \times \R^{d\times d} \to \R \cup \{\infty\}$, is a key ingredient for the SOC theory in the following sections and the backbone of our numerical scheme. For the problem defined in <ref> and <ref> it reads

\begin{equation}\label{eq:H_stoc}
\begin{split}
H(s,\bfz,\bfp,\bfM)=\sup_{\bfu\in U} \Hc(s,\bfz,\bfp,\bfM,\bfu),
\end{split}
\end{equation}

where $\bfp$ and $\bfM$ are called adjoint variables and

\[ \Hc(s,\bfz,\bfp,\bfM,\bfu) =  \frac{1}{2}\text{tr}\left(\sigma(s,\bfz)\bfM\right) + \bfp \cdot f(s,\bfz,\bfu)-L(s,\bfz,\bfu).  \]

To make it notationally convenient, in the rest of the paper, we drop the second argument for the controls and denote controls by $\bfu_{t,\bfx}(s)$.

§.§ Stochastic Pontryagin Maximum Principle

The PMP provides first-order necessary conditions for the SOC problem and also states that the optimal control $\bfu_{t,\bfx}^*$ must satisfy an (extended) Hamiltonian system along the optimal state and adjoint trajectory.
This is made precise by the following result from <cit.>.
Assume that $(\bfz^*_{t,\bfx},\bfu^*_{t,\bfx})$ is an optimal pair that solves <ref> and <ref>. Then there exist adjoint states $\bfp_{t,\bfx} \colon [t,T]\to \R^d$ and $\bfM_{t,\bfx} \colon [t,T]\to \R^{d \times d}$ satisfying the adjoint equation

\begin{equation}\label{eq:adjoint}
\begin{cases}
\du \bfp_{t,\bfx}(s)= & \bfM_{t,\bfx}(s) \du W(s) - \nabla_{\bfz} \mathcal{H} \big( s,\bfz^*_{t,\bfx}(s), \bfp_{t,\bfx}(s), \bfM_{t,\bfx}(s),\bfu^*_{t,\bfx}(s)) \du s \\
\bfp_{t,\bfx}(T)= & - \nabla_{\bfz} G \big(\bfz_{t,\bfx}^*(T) \big),
\end{cases}
\end{equation}

where $s \in [t,T]$ and the optimal control satisfies

\begin{equation}\label{eq:PMP}
\begin{split}
\bfu^*_{t,\bfx}(s) \in \argmax_{\bfu \in U}~ \mathcal{H} \big( s,\bfz_{t,\bfx}^*(s), \bfp_{t,\bfx}(s), \bfM_{t,\bfx}(s),\bfu(s) \big)
\end{split}
\end{equation}

for almost all $s\in [t,T]$, $\bbP$-almost surely.

We note that the optimal control defined in <ref> only depends on the adjoint variable $\bfp_{t,\bfx}$ but not on $\bfM_{t,\bfx}$ since the diffusion $\sigma(\cdot,\cdot)$ does not depend on the control in our case.

We further assume that there exists a unique continuous closed form solution to <ref>. Although not demonstrated in this work, this assumption can be weakened to include implicitly defined functions as long as they can be obtained efficiently; this allows, for example, modeling more general convex running costs.

We note that when the control satisfies <ref>, the dynamics in <ref> is equal to

\begin{align}\label{eq:charOC_stoc_red}
\left\{ \begin{aligned}
\du \bfz_{t,\bfx}^*(s)& = \nabla_{\bfp} \Hc \big( s,\bfz_{t,\bfx}^*(s), \bfp_{t,\bfx}(s), \bfM_{t,\bfx}(s), \bfu^*_{t,\bfx}(s) \big) \du s	+ \sigma( s,\bfz_{t,\bfx}^*(s)) \du W(s), \\ %\; t \leq  s \leq  T, \\ %
\bfz_{t,\bfx}^*(t) &= \bfx.
\end{aligned}
\right.
\end{align}

The system of equations <ref> is called the stochastic Hamiltonian system, where the maximum condition <ref> corresponds to the variational inequality for the control.

Finding a tuple
$(\bfz^*_{t,\bfx},\bfu^*_{t,\bfx},\bfp_{t,\bfx},\bfM_{t,\bfx})$ that satisfies the PMP can be extremely difficult.
However, when the value function $\Phi$ is differentiable, $(\bfp_{t,\bfx},\bfM_{t,\bfx})$ satisfying <ref> can be obtained from $\Phi$; this is formalized in the following theorem that, with weaker assumptions, can be found in  <cit.>.
Assume that $\bfu^*_{t,\bfx}$ is an optimal control and $\Phi\in C^{1,3}([t,T]\times\R^d)$. Then

\begin{equation}\label{eq:gradPhi}
\bfp_{t,\bfx}(s) = -\nabla_{\bfz} \Phi \big(s,\bfz_{t,\bfx}^*(s)\big)
\quad \text{ and } \quad
\bfM_{t,\bfx}(s))= - \sigma(s,\bfz_{t,\bfx}^*(s))^\top \nabla^2_{\bfz} \Phi \left(s,\bfz_{t,\bfx}^*(s)\right)
\end{equation}

solve <ref>.

<Ref> simplifies <ref> to

\begin{equation}\label{eq:opt_control_nec}
\begin{split}
\bfu^*_{t,\bfx}(s)=\bfu_{t,\bfx}^*\left(s,\bfz^*_{t,\bfx}(s),-\nabla_{\bfz} \Phi \big(s,\bfz_{t,\bfx}^*(s)\big)\right).
\end{split}
\end{equation}

This relation along with <ref> is one of the key ingredients of our numerical solution approach.
<Ref> characterizes optimal control in a feedback or closed-loop form, which is of utmost importance in many real-life applications. This is because optimal control can be quickly computed at any given point in time and space,  when the value function is known and its gradient is readily available. Thereby, avoiding re-computation of optimal controls for cases when multiple evaluations are needed for different time or states.

§.§ Hamilton-Jacobi-Bellman Equation

To help approximate  the value function $\Phi$, we also use the fact that  $\Phi$ satisfies the Hamilton-Jacobi-Bellman (HJB) equation, which is a result of the Dynamic Programming (DP) method or Bellman's principle. We state the following result taken from [37] under suitable assumptions, see also <cit.>.
Assume that the value function $\Phi \in C^{1,2}([t,T] \times \R^d)$. Then $\Phi$ satisfies the HJB equation

-   ∂_s Φ(s,)  +  H(s,,-∇_ Φ(s,),  -σ(s,)^⊤∇_^2 Φ(s,))=0,  ∀ (s,)∈[t,T)×^d,

Φ(T,) = G().

The smoothness of $\Phi$ can be relaxed to continuity in the weaker sense of viscosity solutions <cit.>.

Using the relation between PMP and DP from <ref> the HJB equation, <ref>, can now be re-written as the following second order parabolic PDE:

-∂_s Φ(s,) -12tr(σ(s,)σ(s,)^⊤∇^2 Φ(s,)) - ∇Φ(s,)·f(s,,^*)

-∂_s Φ(t,) -σ^22ΔΦ(t,) + ∇Φ(t,)--L(s,,^*)=0,  ∀ (s,)∈[t,T)×^d,

Φ(T,) = G((T)).

In addition, by the envelope theorem, it follows that $\nabla_{\bfp} \Hc = \nabla_{\bfp} H$ and $\nabla_{\bfM} \Hc = \nabla_{\bfM} H$. This simplifies the computation of optimal trajectories, which can now be expressed via the value function as:

\begin{align}\label{eq:HJB_FSDE}
\left\{ \begin{aligned}
\du \bfz_{t,\bfx}^*(s)& = \nabla_{\bfp} H \big( s,\bfz_{t,\bfx}^*(s), -\nabla_{\bfz} \Phi \big(s,\bfz_{t,\bfx}^*(s)\big), -\sigma(s,\bfz_{t,\bfx}^*(s))^\top \nabla^2_{\bfz} \Phi \big(s,\bfz_{t,\bfx}^*(s)\big) \big) \du s \\  & \phantom{\nabla_{\bfp} H \big( s,\bfz_{t,\bfx}^*(s), -\nabla_{\bfz} \Phi \big(s,\bfz_{t,\bfx}^*(s)\big)\big)\sigma(s,\bfz_{t,\bfx}^*(s))} + \sigma(s,\bfz_{t,\bfx}^*(s))  \du W(s), \\
\bfz_{t,\bfx}^*(t) &= \bfx.
\end{aligned}
\right.
\end{align}

These modified dynamics do not explicitly involve the control, which reduces the problem solely to state variable.
This idea was proposed for deterministic optimal control problems in [27].

§.§ FBSDE Formulation

One way to avoid the need for a spatial discretization of the HJB equation <ref> is to use a non-linear version of Feynman-Kac formula and obtain an equivalent system of stochastic differential equations.
This idea has been applied to a variety of nonlinear parabolic/elliptic PDEs; see,  e.g., [1, 24, 5, 28, 30].

Our FBSDE system uses  <ref> as the forward system to sample trajectories.
Along those trajectories, we note that the solution to the HJB equation <ref> must satisfy the backward SDE

Φ(s,(s))    = ∇_ Φ(s,(s))^⊤σ(s,(s))  W(s)
- L(s,(s), ^*(s))s,

Φ(T,(T))    = G((T)).

It is important to stress that our choice of the forward system <ref> is the key difference to other existing works.
For example,  [12] and [33] use the standard Brownian motion.
While both of these choices lead to a valid FBSDE system for <ref>, we advocate for including the control in the dynamics as motivated by stochastic PMP <ref>. As our numerical experiments demonstrate, focusing the sampling along optimal trajectories can lead to more accurate and efficient value function approximations.

§ NEURAL NETWORK APPROACH

In this section, we present a neural network framework for approximating the value function of the stochastic optimal control problem defined by the objective functional <ref> and dynamics <ref>.
The theoretical foundation of our framework are given by the PMP, FBSDE, and Dynamic Programming as presented in the previous section.
The key idea is to approximate the value function $\Phi$ in <ref> by a neural network and compute the control using the feedback form <ref>.
What sets our framework apart is the use of the feedback form to guide the sampling during training.
Thereby we seek to learn to explore the relevant part of the state space.
We also derive and experiment with various loss functions that are based on the control objective <ref>, the BSDE <ref>, and the  HJB <ref>.

§.§ Neural Network Approximation

The first building block of our framework is to parameterize the value function using a neural network.
Since finding an effective network architecture for any learning task is both crucial and an open research topic, we treat this as a modular component.
Our framework can be used with any  scalar-valued neural network that takes inputs in $\R^{d+1}$ as long as it is twice differentiable with respect to its inputs.

Among the networks we use in our numerical experiments is the Multi-Layer Perceptron (MLP) model used in [33].
As an alternative, which also satisfies the regularity needed, we propose the residual network also used for deterministic control in [27]. The network is given by

\begin{equation}
\label{eq:NNArchitecture}
% \begin{split}
\Phi(\bfy ; \bfth) = \bfw^\top \mathcal{NN}(\bfy;\bfth_\mathcal{NN})
+ \frac{1}{2} \bfy^\top (\bfA^\top
\bfA)\bfy + \bfb^\top \bfy + c,
% \end{split}
\end{equation}

with trainable weights $\bfth = (\bfw, \bfth_\mathcal{NN}, \bfA, \bfb, c)$.
Here the inputs $\bfy=(s,\bfz(s))\in \R^{d+1}$ correspond to time-space,  $\mathcal{NN}(\bfy;\bfth_\mathcal{NN}) \colon \R^{d+1} \to \R^m$ is a neural network,
and $\bfth$ contains the trainable weights: $\bfw\,\,{\in}\,\,\R^m$, $\bfth_\mathcal{NN}\,\,{\in}\,\,\R^p$
, $\bfA\,\,{\in}\,\,\R^{\gamma \times (d+1)}$, $\bfb\,\,{\in}\,\,\R^{d+1}$, $c\,{\in}\,\R$, where rank $\gamma{=}\min (10,d+1)$ limits the number of parameters in $\bfA^\top \bfA$.
Here, $\bfA$, $\bfb$, and $c$ model quadratic potentials, i.e., linear dynamics; $\mathcal{NN}$ models nonlinear dynamics.

In our experiments, for $\mathcal{NN}$, we either use a MLP [11]

\begin{equation} \label{eq:MLP}
\begin{split}
\bfa_0 & = \rm act(\bfK_0 \bfy + \bfb_0) \\
\bfa_{i+1} & = \rm act(\bfK_{i+1} \bfa_i + \bfb_{i+1}),\quad 0\le i \le M-2 \\
\mathcal{NN}(\bfy;\bfth_\mathcal{NN}) & =  \rm act(\bfK_M \bfa_{M-1} + \bfb_M),
\end{split}
\end{equation}

or a residual neural network (ResNet) [14]

\begin{equation} \label{eq:ResNet}
\begin{split}
\bfa_0 & = \rm act(\bfK_0 \bfy + \bfb_0) \\
\bfa_{i+1} & = \bfa_i + \rm act(\bfK_{i+1} \bfa_i + \bfb_{i+1}),\quad 0\le i \le M-2 \\
\mathcal{NN}(\bfy;\bfth_\mathcal{NN}) & = \bfa_{M-1} + \rm act(\bfK_M \bfa_{M-1} + \bfb_M),
\end{split}
\end{equation}

with neural network weights $\bfth_\mathcal{NN}{=}(\bfK_0,\hdots,\bfK_M, \bfb_0,\hdots,\bfb_M)$ where $\bfb_i\in \R^{m}\; \forall i$, $\bfK_0 \in \R^{m \times (d+1)}$, and $\{ \bfK_1,\hdots,\bfK_M \}\in \R^{m \times m}$ with $M$ being the depth of the network.
The choice of the element-wise nonlinearity $\rm act(\cdot)$ is discussed in the respective  experiments.

§.§ Training Problem

Ideally, we would choose $\bftheta$ such that $\Phi(s,\bfz;\bftheta)$ is equal to the value function of the control problem  globally, that is, for all $(s,\bfz) \in [t,T]\times \R^d$.
Since this is known to be cursed by the dimensionality for reasonable problem sizes, we resort to a semi-global approach which enforces this property at randomly sampled points in the space time domain.

To generate samples, we first obtain initial states $\bfx\sim \rho$ from some (possibly Dirac)  distribution $\rho$ and then use an Euler Maruyama scheme with  $N+1$ equidistant time points $s_0, \ldots, s_N$ and step size $\du s = (T-t)/N$.
This yields a state trajectory starting at $\bfz_0=\bfx$ via

\begin{equation}\label{eq:sampling}
% \begin{align}
\bfz_{i+1} = \bfz_i + f(s_i,\bfz_i,\bfu_i) \du s + \sigma(s_i,\bfz_i) \du \bfW_i, \quad i=0,\ldots,N-1
% \end{align}
\end{equation}

where $\du\bfW_i \sim \mathcal{N}({\bf 0},\du s \cdot \bfI_d)$,  and $\bfu_i = \bfu_{t,x}^*(s_i,\bfz_i)$ is the optimal control obtained from the feedback, that is, form <ref>

\begin{equation*}
\bfu_i^* \in \argmax_{\bfu \in U} \Hc\left(s_i,\bfz_i,-\nabla \Phi(s_i,\bfz_i;\bftheta),-\sigma(s_i,\bfz_i)^\top\nabla^2 \Phi(s_i,\bfz_i;\bftheta), \bfu\right).
\end{equation*}

A few comments are in place. First, it is important to note that due to the feedback form, the sampled trajectories depend on the parameters of the value function.
Second, the addition of this drift term, motivated by control theory, is the key difference to neural network solvers for the more general class of semi-linear elliptic PDEs [12, 33].
Third, the drift term can also be motivated by the fact that for $\sigma \to 0$, the trajectories above approximate the characteristic curves of the non-viscous HJB equation; thereby our SOC approach coincides with that for deterministic OC in [27].

To simplify our notation below, we collect the states, control, and noise along the discrete trajectories in <ref> column-wise in the matrices

\begin{equation*}
\bfZ \in \R^{d\times N}, \quad \bfU \in \R^{k\times N}, \quad \mathbf{dW} \in \R^{d\times N}.
\end{equation*}

To learn the parameters of the neural networks in an unsupervised way (that is,  assuming neither analytic values of $\Phi$ nor optimal control trajectories), we approximately solve the minimization problem

\begin{equation}\label{eq:full-opt}
\begin{split}
\min_{\bfth} \E_{\bfx \sim \rho}  & \left\{\E_{\bfZ,\bfU,\mathbf{dW} | \bfx} \left\{
\beta_1 P_{\rm BSDE}(\bfZ,\bfU,\mathbf{dW}) + \beta_2 P_{\rm HJB}(\bfZ) + \beta_3 J(\bfZ,\bfU) \right.\right.\\
+\beta_4 |G(\bfz_N) - \Phi(s_N, \bfz_N; \bftheta)| +  \beta_5 |\nabla G(\bfz_N) - \nabla\Phi(s_N,\bfz_N; \bftheta)|
\right\} \right\},
\end{split}
\end{equation}

where the  terms in the objective function consist of penalty functions for violations of the BSDE system and the HJB PDE, the control objective, and penalty terms for the terminal condition, respectively, and are defined below. The relative influence of each term is controlled by the components of $\beta \in \R^5_{+}$.
Different choices of $\beta$ allow us to experiment with different learning approaches; for example, setting $\beta_1 = \beta_4 = \beta_5=1$ and $\beta_2=\beta_3=0$ provides the same loss function as in [33] while $\beta_1=0$ and $\beta_i>0, \; i\in \{2,3,4,5\}$ gives the loss function used for deterministic OC problems in [27].

We penalize the violation of the BSDE <ref> via

\begin{equation}
P_{\rm BSDE}(\bfZ,\bfU,\mathbf{dW}) = \sum_{i=0}^{N-1} |\Phi_{i+1}(\bftheta)-\Phi_{i}(\bftheta) + L(s_i,\bfz_i,\bfu_i)\du s -  \nabla\Phi_i(\bftheta)^\top \sigma(s_i,\bfz_i)  \du\bfW_i |^p
\end{equation}

where we use the abbreviations $\Phi_i(\bftheta) := \Phi(s_i,\bfz_i;\bftheta)$ and $\nabla\Phi_i(\bftheta) := \nabla\Phi(s_i,\bfz_i;\bftheta)$ and the exponent $p\in\{1,2\}$.
Similarly, the HJB penalty term reads

\begin{equation}
% \begin{split}
P_{\rm HJB}(\bfZ) = \du s \sum_{i=1}^{N} |H(s_i,\bfz_i,-\nabla \Phi_i(\theta), -\sigma(s_i,\bfz_i)^\top \nabla^2 \Phi_i(\theta)) - \partial_s \Phi_i(\bftheta)|^p,
% \end{split}
\end{equation}

where $\nabla^2 \Phi_i(\theta) := \nabla^2\Phi(s_i,\bfz_i;\theta)$,  $\partial_s \Phi_i(\theta) := \partial_s \Phi(s_i,\bfz_i;\theta)$.
Finally, we approximate the objective functional via

\begin{equation*}
J(\bfZ,\bfU) = G(\bfz_N) + \du s \sum_{i=1}^N L(s_i,\bfz_i,\bfu_i).
\end{equation*}

In principle, any stochastic approximation approach can be used to approximately solve the above optimization problem.
Here, we use Adam and sample a minibatch of trajectories originating in $n$ i.i.d. samples from $\rho$.

§.§ Implementation

We implement and test our proposed approach in two software environments.

To obtain a direct comparison with [33] we modify the FBSNN code accompanying the paper.
To this end, we created a publicly available fork at <https://github.com/EmoryMLIP/FBSNNs>.
Our main modifications are adding the control objective in the training loss and also adding the proposed drift to the forward dynamics.
The experiments carried out here use a standard MLP as the neural network without the quadratic term added in <ref>.

In order to further simplify the experimentation, we also implement our own PyTorch code available at <https://github.com/EmoryMLIP/NeuralHJB>. Our implementation contains all loss terms in <ref>.
We implement both sampling techniques: pure random walk and the proposed one informed by PMP.
This facilitates comparisons of  our approach with other available methods and simplifies developing new examples.

We tested most of our examples using either Intel Xeon E5-4627 CPU or Nvidia P100 GPU.

§ NUMERICAL EXPERIMENTS

We test the efficacy of our proposed algorithm on two different SOC problems.
In <ref>, we  introduce a two-dimensional trajectory planning problem to visualize the difference between purely random exploration and our proposed sampling scheme.
In <ref>, we compare our approach to those in [7, 12] using a  100-dimensional benchmark problem. For the original version of this problem, our method shows faster initial convergence and time-to-solution but comparable accuracy.
We modify the terminal cost of this problem to further highlight the importance of the feedback form in the sampling  (see <ref>) and experimentally show that the growth of network size is not exponential with increasing dimensions.

§.§ 2D Trajectory Planning Problem

To visualize the behaviour of our PMP-based sampling approach, we consider a two-dimensional test problem.

The problem consists of planning an optimal trajectory from the initial state $\bfx \sim \rho = \mathcal{N}((-1.5,-1.5)^\top, \bfI_2)$ to the target $\bfx_{\rm target} = (1.5,1.5)^\top$.
To make the problem interesting, a hill is placed at the origin, which adds height-dependent cost for traveling around that region given by

\begin{equation*}
Q(\bfz) = 50 \exp\left(- \frac{\|\bfz\|^2}{0.8} \right).
\end{equation*}

The dynamics for the problem read

\begin{equation}
f(s,\bfz,\bfu) = \bfu \quad \text{and} \quad \sigma = 0.5.
\end{equation}

The running cost and terminal cost of the problem are given respectively, by

\begin{equation}
L(s,\bfz,\bfu) = \frac{1}{2} \| \bfu \|^2 + Q(\bfz) \quad \text{and} \quad  G(\bfz) = 50 \cdot \| \bfz - \bfx_{\rm target} \|^2.
\end{equation}

The corresponding HJB equation is

\begin{equation}
\partial_s\Phi (s, \bfz) + \frac{1}{8} \Delta \Phi(s, \bfz) - \frac{1}{2}\| \nabla \Phi(s, \bfz) \|^2 + Q(\bfz) = 0.
\end{equation}

With the above calculations, the sampling SDE <ref> simplifies to

\begin{equation}
\label{eq: 2D_FSDE}
\bfz_{i+1} = \bfz_{i}- \nabla \Phi(s_i,\bfz_{i})\du s + \frac{1}{2} \du \bfW_i.
\end{equation}

We approximate the value function using the model in <ref> where $\mathcal{NN}$ is a two-layer residual neural network with 32  neurons per layer and the rank of the quadratic term is $\gamma=3$.
Overall, the model has 1229 trainable weights.
As activation function we chose the anti-derivative of the $\tanh$ function as in [27].
We use the penalty parameters $\beta=(1,1,1,0.1,0.1)$, that is, we enable both the penalty terms, $P_{\rm BSDE}$ and $P_{\rm HJB}$.
To approximately solve <ref> we use 3000 steps of the Adam optimizer and batch size of 64. After 1500 iterations, we reduce the learning rate from 0.01 to  0.001. The average cost per iteration is about  $ 0.1$s using the CPU.

In <ref>, we visualize the estimate of the value function and the optimal control policy at $t=0$ as well as example trajectories for randomly chosen initial states after training the model.
As expected, the trajectories travel from the initial points to the target while avoiding the obstacle in the center of the domain.

Results of the two-dimensional test problem. Left: Value function function estimate $\Phi(0,\cdot;\bftheta^*)$ after training. Middle: Quiver plot of optimal controls  at $s=0$. Right:  Trajectories from randomly chosen initial states.

To gain more insight into the sampling, we store all states visited during training and plot them as two-dimensional histograms for different time points (left to right) in  <ref>.
We compare the proposed PMP-based sampling (<ref>) to the pure random walk (i.e., without drift) used in [12, 33] (<ref>).
As to be expected, sampling without drift leads to visiting points around the initial states in all (even suboptimal) directions.
It can also be seen that there are almost no samples close to the target.
With the drift term and after training, the sampled states cover more likely paths between the initial and target states.

Training samples with pure random walk as FSDE  as also used in [12] and [33].

Training samples with PMP-based drift term.

We visualize training samples of a pure random walk sampler (top row) and our proposed PMP-based sampler (bottom row) for the two-dimensional test problem.
At six time points (left to right), we visualize the sampled states as two-dimensional histrograms.
As expected, the pure random walk explores the area around the initial state in all (even suboptimal) directions, while the proposed approach   learns to sample around approximately optimal trajectories.

Another way to interpret the histogram plots in <ref> is that they show the semi-global nature of neural network approaches for SOC problems and, crucially, that the regions of the state space  visited in training can vary a lot.
Since the loss function in <ref> penalizes the HJB and BSDE loss only around points sampled from the forward SDE, one would expect that the trained model is more reliable in regions that are frequently visited.

§.§ 100-dimensional example

We consider the 100-dimensional benchmark SOC problem also used in [7, 12] with  initial state $\bfx = (0,0,\dots,0)^\top\in \R^{100}$ corresponding to time $t=0$.
The drift and diffusion of the system are given by

\[ f(s,\bfz,\bfu) = 2\bfu\quad \text{and} \quad \sigma = \sqrt{2}, \]

respectively. The terminal and Lagrangian cost  are

\begin{align}\label{eq:100D_term}
G(\bfz) = \ln \left(\frac{1 + \|\bfz \|^2 }{2}\right),\quad \text{and}\quad L(s,\bfz,\bfu ) =   \|\bfu\|^2,
\end{align}

We compute the Hamiltonian <ref> as

\begin{align*}
\begin{split}
H(s,\bfz,\bfp,\bfM) & =\sup_{\bfu\in U}\left\{\frac{\sigma}{2}\text{tr}\left(\bfM\right) + \bfp \cdot f(s,\bfz,\bfu)-L(s,\bfz,\bfu) \right\}\\
& = \sup_{\bfu\in U}\left\{\frac{1}{\sqrt{2}}\text{tr}\left(\bfM\right) + \bfp \cdot 2\bfu - \|\bfu\|^2 \right\}
\end{split}
\end{align*}

Using the first-order necessary condition we get

\begin{align*}
0 = 2\bfu - 2\bfp\implies \bfu = \bfp,
\end{align*}

and using this closed form for $\bfu$, the Hamiltonian is given by

\[ H(s,\bfz,\bfp,\bfM) =\frac{1}{\sqrt{2}} \text{tr}\left(\bfM\right) + \|\bfp\|^2 .\]

Hence, the HJB equation satisfied by the value function, $\Phi(\cdot,\cdot)$, reads

\begin{equation}\label{eq:HJB_100D}
\frac{\partial }{\partial s} \Phi(s,\bfz) + \Delta \Phi(s,\bfz) - \| \nabla \Phi(s,\bfz) \|^2 = 0, \quad \Phi(T,\bfz) = G(\bfz).
\end{equation}

and its solution is given by

\begin{equation}\label{eq:HJB_100D_sol}
\Phi(s,\bfz) = -\ln\left(\bbE \left( \exp\left(-G\left(\bfz + \sqrt{2}\,\du W(T-s)\right)\right) \right)\right),
\end{equation}

which we use to test the performance of our method.

Finally, we note that the forward SDE <ref> we propose to use for sampling the state space simplifies to

\begin{equation}\label{eq:FBSDE100d_ours}
\bfz_{i+1} = \bfz_{i}- 2\, \nabla_{\bfz}\Phi(s_i,\bfz_i)\du s + \sqrt{2} \du \bfW_{i}.
\end{equation}

§.§.§ The importance of sampling

To demonstrate the impact of using the feedback form to sample the state space, we use the same neural network model as in [33], which is given by a 5-layer feed forward neural network with 256 neurons per hidden layer to approximate the solution $\Phi(s,\bfz)$. We partition the time interval $[0,1]$ using 50-uniformly spaced points.
We use the same penalty parameters as in the original code, that is, $\beta=(1,0,20,1,1)$.
We use the Adam optimizer  [18] to update the parameters of the network with a batch size of 64 using 50k iterations. The average cost per 100 iterations was  $27$s using the CPU. For the following experiments, we use <ref> excluding $P_{\rm HJB}$ penalty and compare our method with FBSNNs in [33].

Solution to <ref> obtained using our method and method in [33]

Method       2c20k iterations

RE     RE$_0$

2c50k iterations

RE     RE$_0$

0.54%     0.12%

0.39%     0.045%

0.48%     0.0083%

0.39%     0.012%

Relative errors for <ref>  obtained using our method and method in [33]

In <Ref> we plot the exact solution (black-dashed line) <ref>, the learned solution using our approach (blue-solid line) and the solution learned using FBSNNs (red-solid line) along five random trajectories. In the top row, we present the results obtained after training the networks for 20K iterations with a learning rate of $10^{-3}$ and the bottom row presents the results after training the networks for 20K and 30K iterations with learning rates $10^{-3}$ and $10^{-4}$, respectively. These results suggest that our approach is approximating the value function better, especially in early iterations, as compared to the FBSNNs.

In <ref>, we also compare the learned solutions to the exact solution $\Phi$ in <ref> by computing the average relative errors,

\[ RE = \frac{\|\Phi(\cdot,\cdot;\bfth) - \Phi(\cdot,\cdot) \|_{2}}{\|\Phi\|_2},\quad   RE_0 = \frac{|\Phi(0,\bfz(0);\bfth) - \Phi(0,\bfz(0)) |}{|\Phi(0,\bfz(0))|},\]

for ten random trajectories.
Our method attains lower errors, especially for the initial values and at the earlier iterations.

§.§.§ Shifted target

In the example above, the minimizer of the terminal function coincides with the initial state $\bfx=(0,0,\dots,0)^\top$.
Therefore, even a random walk without drift (as used in  [33, 12]) will sample around the optimal terminal state, which is critical to accurately approximate the value function.
This also means that after training using our approach, the drift term in the sampler is relatively small and that the above experiment does not fully show the advantages of our method.

To shed more light on the importance of sampling, we modify the terminal cost to

$$G(z) = 1000\ln \left(\frac{1 + \|\bfz - \bfz_{\rm target}\|^2 }{2}\right),$$

with $\bfz_{\rm target} = [3,3,\dots, 3]^{T}$, so that the target for the state variable $\bfz$ at final time $T$ no longer coincides with the initial state.
Similar to the two-dimensional test problem in <ref>, solving the modified problem now requires sampling around the target and we expect to benefit from the added drift term.

We compare our method to FBSNNs on the modified problem keeping the same network structure and hyper-parameters. We use a smaller $\sigma = \frac{2\sqrt{2}}{5}$ to improve training speed. We evaluate the performance of the methods using the objective functional $J$ defined in  <ref> at the control obtained from the feedback form via the respective value function approximations. For this experiment we use GPU to train and the results of this comparison are shown in <ref>.

Computational results for the modified 100-dimensional benchmark problem in <ref>. Left: Control objective for both methods given the same initial state, blue line represents results using FBSNNs in [33], and orange line denotes our method.  Right: Trajectory examples generated using learned value functions on two randomly selected dimensions. Orange line represents our method and blue line FBSNNs.

To reduce the effect of the Brownian motion, we run the experiments for each method on the same problem five times and plot the average values corresponding to training iterations, we also further reduced the effect of the random walk in evaluation for clearer visualization.

As can be seen in <ref> (left), our method not only yields faster initial convergence but also achieves a considerably lower control objective. This indicates that the controls obtained from our approach are  more effective, that is, they are closer to optimal.
Since FBSNNs use a Brownian motion with no drift, the sampling is unlikely to discover the target.
Consequently, the generated trajectories in <ref> (right)  from our method approximately reach the target, while the trajectories obtained from the FBSNN method stay closer to the initial state.

§.§.§ Scalability with problem dimension

We compare the necessary number of weights required to train a working model for problem with different sizes. Each problem follows the same setup used in <ref> with the only difference being the dimensionality of the problem. For a model to be considered working we require that under the feedback form the expected value for the final state of the agent to be within 0.02 range of the target.

As some evidence that the proposed NN approach can mitigate the curse of dimensionality, we vary the dimension in $\{20,40,60,80,100\}$ in <ref>.
We use the same network architecture from the previous session.
For each dimension, we increase the width of the network until the training yields a network whose trajectories end sufficiently close to the target.

With our method we notice that the number of weights required to solve the problem does not increase exponentially with dimensionality, suggesting that our method does ease the CoD. It is however worth mentioning that this experiment does not aim to determine the minimum number of weights required to solve each problem.

§ DISCUSSION

We propose a neural network approach for approximately solving Hamilton-Jacobi-Bellman PDEs arising in high-dimensional stochastic optimal control.
Similar to existing approaches [33, 12], we parameterize the value function with a neural network and we experiment with different losses to train the weights of the network.
What sets our work apart from existing approaches is the use of feedback form given by the stochastic Pontryagin maximum principle to design the forward SDE used to explore the state space during training.

Using an intuitive two-dimensional test problem, we visualize that the improved sampling strategy allows us to effectively learn the value function and determine the relevant regions of the state space; see <ref>.
Based on this insight, we modify the  $100$-dimensional test problem also used in [33, 12] by shifting the minimizer of the terminal costs; see <ref>.
Thereby, we demonstrate that our proposed method dramatically improves the quality of the obtained control.

In additional numerical experiments, we show that the modified forward SDE and the control objective can lead to faster initial convergence compared to the approaches in [33, 12] (see <ref>) and show that, for the problem at hand, the network size needed to achieve a given tolerance does not need to grow exponentially.

We refer to our method as a semi-global method for solving the HJB equation since we do not aim at estimating the value function accurately globally.
Instead, we seek to approximate the value function well for states likely to be visited by optimal trajectories of the SOC problem.
In theory, any point in the state space has a positive probability to be visited  due to the stochasticity in the dynamics.
However, the histogram plots in <ref> suggest that the density of the optimal trajectory is concentrated in a small subset of the state space.
Therefore, focussing the exploration to this subsets may have practical advantages over FSDEs that are pure random walks.

Another benefit of our proposed forward SDE compared to purely random exploration is that it coincides with the characteristic curves of the HJB equation as the stochasticity of the system is reduced.
Therefore, our work can be seen as an extension of the neural network approaches for deterministic control problems in [27].

Compared to neural network approaches for semi-linear elliptic/parabolic PDEs such as [33, 12] it is important to highlight that our approach is only applicable to HJB equations arising in stochastic optimal control.
Since our forward SDE is derived from optimality principles, extending it to other high-dimensional PDEs  (e.g., Black Scholes and Allen Cahan equations) is not obvious and may be impossible.

In future work, we will apply our approach to problems with non-constant and non-scalar diffusion coefficients.
Although the theoretical framework supports this, we are not aware of any practical algorithms for this case.
To overcome the challenge of efficiently computing Hessians of the value function needed for HJB penalties, we will use the hessQuik package [26].

[1]
F. Antonelli, Backward-forward stochastic differential equations,
Ann. Appl. Probab., 3 (1993), pp. 777–793,
<http://links.jstor.org/sici?sici=1050-5164(199308)3:3<777:BSDE>2.0.CO;2-5   origin=MSN>.

[2]
R. Bellman, Dynamic programming and stochastic control processes,
Information and Control, 1 (1958), pp. 228–239.

[3]
C. Bender and R. Denk, A forward scheme for backward SDEs,
Stochastic Process. Appl., 117 (2007), pp. 1793–1812,

[4]
D. Bertsekas and J. Tsitsiklis, Neuro-dynamic programming: an
overview, in Proceedings of 1995 34th IEEE Conference on Decision and
Control, vol. 1, 1995, pp. 560–564 vol.1,

[5]
P. Cheridito, H. M. Soner, N. Touzi, and N. Victoir, Second-order
backward stochastic differential equations and fully nonlinear parabolic
PDEs, Comm. Pure Appl. Math., 60 (2007), pp. 1081–1110,

[6]
H. Dong and N. V. Krylov, The rate of convergence of
finite-difference approximations for parabolic Bellman equations with
Lipschitz coefficients in cylindrical domains, Appl. Math. Optim., 56
(2007), pp. 37–66, <https://doi.org/10.1007/s00245-007-0879-4>,

[7]
W. E, J. Han, and A. Jentzen, Deep learning-based numerical methods
for high-dimensional parabolic partial differential equations and backward
stochastic differential equations, Commun. Math. Stat., 5 (2017),
pp. 349–380, <https://doi.org/10.1007/s40304-017-0117-6>,

[8]
I. Exarchos, E. Theodorou, and P. Tsiotras, Stochastic differential
games: a sampling approach via FBSDEs, Dyn. Games Appl., 9 (2019),
pp. 486–505, <https://doi.org/10.1007/s13235-018-0268-4>,

[9]
I. Exarchos and E. A. Theodorou, Stochastic optimal control via
forward and backward stochastic differential equations and importance
sampling, Automatica J. IFAC, 87 (2018), pp. 159–165,

[10]
W. H. Fleming and H. M. Soner, Controlled Markov Processes and
Viscosity Solutions, vol. 25 of Stochastic Modelling and Applied
Probability, Springer, New York, second ed., 2006.

[11]
I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning, MIT
Press, 2016.

[12]
J. Han, A. Jentzen, and W. E, Solving high-dimensional partial
differential equations using deep learning, Proc. Natl. Acad. Sci. USA, 115
(2018), pp. 8505–8510, <https://doi.org/10.1073/pnas.1718942115>,

[13]
K. P. Hawkins, A. Pakniyat, E. Theodorou, and P. Tsiotras, Forward-backward rapidly-exploring random trees for stochastic optimal
control, 2020, <https://doi.org/10.48550/ARXIV.2006.12444>,

[14]
K. He, X. Zhang, S. Ren, and J. Sun, Deep residual learning for
image recognition, in Proceedings of the IEEE Conference on Computer Vision
and Pattern Recognition (CVPR), June 2016.

[15]
C. Huré, H. Pham, and X. Warin, Deep backward schemes for
high-dimensional nonlinear PDEs, Math. Comp., 89 (2020), pp. 1547–1579,

[16]
E. R. Jakobsen, On the rate of convergence of approximation schemes
for Bellman equations associated with optimal stopping time problems,
Math. Models Methods Appl. Sci., 13 (2003), pp. 613–644,

[17]
S. Ji, S. Peng, Y. Peng, and X. Zhang, Three algorithms for solving
high-dimensional fully coupled fbsdes through deep learning, IEEE
Intelligent Systems, 35 (2020), pp. 71–84,

[18]
D. P. Kingma and J. Ba, Adam: A method for stochastic optimization,
arXiv preprint arXiv:1412.6980,  (2014).

[19]
P. E. Kloeden and E. Platen, Numerical solution of stochastic
differential equations, vol. 23 of Applications of Mathematics (New York),
Springer-Verlag, Berlin, 1992,

[20]
N. V. Krylov, The rate of convergence of finite-difference
approximations for Bellman equations with Lipschitz coefficients, Appl.
Math. Optim., 52 (2005), pp. 365–399,

[21]
K. Kunisch and D. Walter, Semiglobal optimal feedback stabilization
of autonomous systems via deep neural network approximation, ESAIM Control
Optim. Calc. Var., 27 (2021), pp. Paper No. 16, 59,

[22]
H. J. Kushner, Numerical methods for stochastic control problems in
continuous time, SIAM J. Control Optim., 28 (1990), pp. 999–1048,
<https://doi.org/10.1137/0328056>, <https://doi.org/10.1137/0328056>.

[23]
H. J. Kushner and P. Dupuis, Numerical methods for stochastic
control problems in continuous time, vol. 24 of Applications of Mathematics
(New York), Springer-Verlag, New York, second ed., 2001,
Stochastic Modelling and Applied Probability.

[24]
J. Ma, P. Protter, and J. M. Yong, Solving forward-backward
stochastic differential equations explicitly—a four step scheme, Probab.
Theory Related Fields, 98 (1994), pp. 339–359,

[25]
J. Ma and J. Yong, Forward-backward stochastic differential
equations and their applications, vol. 1702 of Lecture Notes in Mathematics,
Springer-Verlag, Berlin, 1999.

[26]
E. Newman and L. Ruthotto, `hessquik`: Fast hessian computation of
composite functions, Journal of Open Source Software, 7 (2022), p. 4171,

[27]
D. Onken, L. Nurbekyan, X. Li, S. W. Fung, S. Osher, and L. Ruthotto,
A neural network approach for high-dimensional optimal control, 2021,

[28]
E. Pardoux and S. Tang, Forward-backward stochastic differential
equations and quasilinear parabolic PDEs, Probab. Theory Related Fields,
114 (1999), pp. 123–150, <https://doi.org/10.1007/s004409970001>,

[29]
M. Pereira, Z. Wang, T. Chen, E. Reed, and E. Theodorou, Feynman-kac
neural network architectures for stochastic control using second-order fbsde
theory, in Proceedings of the 2nd Conference on Learning for Dynamics and
Control, A. M. Bayen, A. Jadbabaie, G. Pappas, P. A. Parrilo, B. Recht,
C. Tomlin, and M. Zeilinger, eds., vol. 120 of Proceedings of Machine
Learning Research, PMLR, 10–11 Jun 2020, pp. 728–738,

[30]
H. Pham, Continuous-time stochastic control and optimization with
financial applications, vol. 61 of Stochastic Modelling and Applied
Probability, Springer-Verlag, Berlin, 2009,

[31]
L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F.
Mishchenko, The Mathematical Theory of Optimal Processes, Translated
by K. N. Trirogoff; edited by L. W. Neustadt, Interscience Publishers John
Wiley & Sons, Inc.  New York-London, 1962.

[32]
W. B. Powell, Approximate dynamic programming, Wiley Series in
Probability and Statistics, Wiley-Interscience [John Wiley & Sons], Hoboken,
NJ, 2007, <https://doi.org/10.1002/9780470182963>,
Solving the curses of dimensionality.

[33]
M. Raissi, Forward-backward stochastic neural networks: Deep
learning of high-dimensional partial differential equations, arXiv preprint
arXiv:1804.07010,  (2018).

[34]
R. F. Stengel, Optimal control and estimation, Dover Publications,
Inc., New York, 1994.
Corrected reprint of the 1986 original.

[35]
J. Wang and P. A. Forsyth, Maximal use of central differencing for
Hamilton-Jacobi-Bellman PDEs in finance, SIAM J. Numer. Anal., 46
(2008), pp. 1580–1601, <https://doi.org/10.1137/060675186>,

[36]
S. Yensiri and R. J. Skulkhu, An investigation of radial basis
function-finite difference (rbf-fd) method for numerical solution of elliptic
partial differential equations, Mathematics, 5 (2017),

[37]
J. Yong and X. Y. Zhou, Stochastic controls, vol. 43 of
Applications of Mathematics (New York), Springer-Verlag, New York, 1999,
Hamiltonian systems and HJB equations.