# A novel control method for solving high-dimensional Hamiltonian systems
through deep neural networks

Shaolin Ji  Zhongtai Securities Institute for Financial Studies, Shandong
University, 250100, China  Shige Peng  School of Mathematics, Shandong
University, 250100, China  Ying Peng  Zhongtai Securities Institute for
Financial Studies, Shandong University, 250100, China  Xichuan Zhang  School
of Mathematics, Shandong University, 250100, China

###### Abstract

In this paper, we mainly focus on solving high-dimensional stochastic
Hamiltonian systems with boundary condition, which is essentially a Forward
Backward Stochastic Differential Equation (FBSDE in short), and propose a
novel method from the view of the stochastic control. In order to obtain the
approximated solution of the Hamiltonian system, we first introduce a
corresponding stochastic optimal control problem such that the extended
Hamiltonian system of the control problem is exactly what we need to solve,
then we develop two different algorithms suitable for different cases of the
control problem and approximate the stochastic control via deep neural
networks. From the numerical results, comparing with the Deep FBSDE method
developed previously from the view of solving FBSDEs, the novel algorithms
converge faster, which means that they require fewer training steps, and
demonstrate more stable convergences for different Hamiltonian systems.

Keywords stochastic Hamiltonian system, FBSDE, optimal control, PDE

##  1 Introduction

The theory of the Hamiltonian system is known as one of the dominant tools for
the description of dynamic phenomenons in the field of physics and economics
[1]. For example, in physics, the mechanical and electrical systems are
usually represented as energy functions, which are at the same time
Hamiltonian systems. Actually, the Hamiltonian system could reflect the laws
of energy conservation and dissipation [1, 2].

A determined Hamiltonian system can be given as

$\begin{cases}\mathrm{d}x_{t}=H_{y}(x_{t},y_{t})\mathrm{d}t,\vspace{1ex}\\\
\mathrm{d}y_{t}=H_{x}(x_{t},y_{t})\mathrm{d}t,\end{cases}$ (1.1)

where $H(x,y):\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a
given real function called the Hamiltonian, $H_{x}(\cdot)$ and $H_{y}(\cdot)$
are partial derivatives of $H(\cdot)$ with respect to $x$ and $y$,
respectively. When considering a terminal condition $y_{T}=\Phi_{x}(x_{T})$
for a given function $\Phi:\mathbb{R}^{n}\rightarrow\mathbb{R}$, (1.1) becomes
a boundary problem.

For more complex environments where the physical system can not be represented
with deterministic form, the Hamiltonian system is usually combined with a
stochastic process. Here we consider a boundary problem of stochastic
Hamiltonian system, as shown in the following,

$\begin{cases}\mathrm{d}x_{t}=H_{y}(t,x_{t},y_{t},z_{t})\mathrm{d}t+H_{z}(t,x_{t},y_{t},z_{t})\mathrm{d}B_{t},\vspace{1ex}\\\
-\mathrm{d}y_{t}=H_{x}(t,x_{t},y_{t},z_{t})\mathrm{d}t-z_{t}\mathrm{d}B_{t},\vspace{1ex}\\\
x_{0}=a,\qquad y_{T}=-\Phi_{x}(x_{T}),\end{cases}$ (1.2)

which is essentially a fully coupled forward-backward stochastic differential
equation (FBSDE in short). Many research work have studied the solutions of
FBSDEs and the eigenvalue of the Hamiltonian systems [3, 4, 5, 6, 7, 8, 9, 10,
11]. The significance of studying this kind of Hamiltonian system is that on
the one hand it can be applied in solving the stochastic optimal control
problems via the well-known stochastic maximum principle [12, 13]; on the
other hand, it helps to obtain the solutions of nonlinear partial differential
equations (PDEs in short) according to the connection between the FBSDEs and
the PDEs [11].

In most cases, it is difficult to obtain the explicit solution of the
Hamiltonian system (1.2), thus numerical methods should be studied. As (1.2)
is essentially a FBSDE, an intuitive way is to solve (1.2) from the
perspective of FBSDEs. Therefore, numerical methods for solving the FBSDEs can
be applied [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 8] , such as the PDE
methods, the probabilistic methods, etc. However, most of the traditional
numerical methods can not deal with high-dimensional problems. Moreover, it is
worth to point out that solving the fully coupled FBSDEs numerically has been
a much more challenging problem than the general FBSDEs, even for low
dimensional cases. Recently, with the application of the deep learning
technique in a wide range of areas, numerical methods based on deep neural
networks have been proposed to solve high-dimensional Backward Stochastic
Differential Equations (BSDEs in short) and FBSDEs and achieved remarkable
success. Among them, a breakthrough work was developed by [24, 25], the main
idea is to reformulate the BSDE into a stochastic optimal control problem by
rewriting the backward process into a forward form and taking the terminal
error as the cost functional, then the solution of the BSDE is approximated by
deep neural network. Other different deep learning algorithms are proposed to
solve the BSDEs and related PDEs [26, 27, 28, 29], where they also focus on
approximating the solution of the BSDE (or PDE) with the deep neural network.
For solving coupled and fully coupled FBSDEs, [30, 31] developed numerical
algorithms which are also inspired by [24, 25].

In this paper, we propose a novel method to solve the Hamiltonian system (1.2)
via deep learning. As equation (1.2) is at the same time a fully coupled
FBSDE, this method is also suitable for solving fully coupled FBSDEs. However,
different from the above mentioned deep learning methods which aim to solve
the FBSDEs directly, we first look for the corresponding stochastic optimal
control problem of the Hamiltonian system, such that the extended Hamiltonian
system of the stochastic control problem is exactly what we need to solve.
Then we approximate the optimal control with deep neural networks. In order to
solve the optimal control problem, two different cases are considered which
correspond to two different algorithms. The first algorithm (Algorithm 1)
deals with the case where the function $f(t,x,u,v)$ defined in (2.5) has an
explicit form. For the case that $f(t,x,u,v)$ cannot be expressed explicitly,
the original control problem is transformed to a double objective optimization
problem and we develop the second algorithm (Algorithm 2) to solve it.
Finally, the numerical solutions $(y_{t},z_{t})$ of (1.2) are obtained by
calculating the solution of the extended Hamiltonian system for the optimal
control according to the stochastic maximum principle.

We also compare the results of our novel proposed algorithms with that of the
algorithm developed in our previous work [31] ( called the Deep FBSDE method
here), which can be used to solve the Hamiltonian system from the view of the
FBSDEs. Comparing with the Deep FBSDE method, the novel algorithms have two
advantages. The first advantage is that less numbers of iteration steps are
required to achieve convergent results. When the Deep FBSDE method converges,
it needs more iterations to achieve a convergent result. The second advantage
is that our proposed algorithms have more stable convergences. For some
Hamiltonian systems, the Deep FBSDE method is easier to diverge with the same
piecewise decay learning rate as our novel proposed algorithms. The details
can be referred to the numerical results in section 4.

This paper is organized as follows. In section 2, we describe the Hamiltonian
system that we aim to solve, and introduce its corresponding stochastic
optimal control problem. In section 3, we introduce two schemes to solve the
stochastic control problem according to whether the function $f(t,x,u,v)$
defined in (2.5) has an explicit form, and then give the corresponding neural
network architectures. The numerical results for different examples are shown
in section 4, and a brief conclusion is made in section 5.

##  2 Problem formulation

In this section, we first introduce the stochastic Hamiltonian system that we
aim to solve, then we show that solving this kind of stochastic Hamiltonian
systems is equivalent to solving a stochastic optimal control problem.

###  2.1 The stochastic Hamiltonian system

Let $T>0$, $(\Omega,\mathcal{F},\mathbb{F},\mathbb{P})$ be a filtered
probability space, in which $B:[0,T]\times\Omega\rightarrow\mathbb{R}^{d}$ is
a $d$-dimensional standard $\mathbb{F}$-Brownian motion;
$\mathbb{F}=\\{\mathcal{F}_{t}\\}_{0\leq t\leq T}$ is the natural filtration
generated by the Brownian motion $B$ Suppose that
$(\Omega,\mathcal{F},\mathbb{P})$ is complete, $\mathcal{F}_{0}$ contains all
the $\mathbb{P}$-null sets in $\mathcal{F}$ and $\mathbb{F}$ is right
continuous.

For $z^{1},z^{2}\in\mathbb{R}^{n\times d}$, define $\left\langle
z^{1},z^{2}\right\rangle=\mbox{tr}(z^{1}(z^{2})^{\operatorname{T}}{})$ and
$|z|^{2}=\left\langle z,z\right\rangle$. The space of all mean square-
integrable $\mathcal{F}_{t}$-adapted and $\mathbb{R}^{n}$-valued processes
will be denoted by $M^{2}(0,T;\mathbb{R}^{n})$, which is a Hilbert space with
the norm

$\|v(\cdot)\|=\Big{(}\mathbb{E}\big{[}\int_{0}^{T}|v(t)|^{2}dt\big{]}\Big{)}^{1/2}$

and

$L^{2}(\Omega,\mathcal{F}_{t},\mathbb{P})\triangleq\left\\{\xi|\xi\in\mathbb{R}^{n}\mbox{
is }\mathcal{F}_{t}\mbox{-measurable and
}\mathbb{E}\left[|\xi|^{2}\right]<\infty\right\\}.$

Let

$H:[0,T]\times\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{n\times
d}\rightarrow\mathbb{R},$ (2.1)

be a $C^{2}$ real function of $(x,y,z)$, called a Hamiltonian and let

$\Phi:\mathbb{R}^{n}\rightarrow\mathbb{R},$ (2.2)

be a $C^{1}$ real function of $x$. In our context, unless otherwise stated, we
always assume that the Hamiltonian $H$ is strictly convex with respect to $y$
and $z$.

Consider the following stochastic Hamiltonian system:

$\begin{cases}\mathrm{d}x_{t}=H_{y}(t,x_{t},y_{t},z_{t})\mathrm{d}t+H_{z}(t,x_{t},y_{t},z_{t})\mathrm{d}B_{t},\vspace{1ex}\\\
-\mathrm{d}y_{t}=H_{x}(t,x_{t},y_{t},z_{t})\mathrm{d}t-z_{t}\mathrm{d}B_{t},\vspace{1ex}\\\
x_{0}=a,\qquad y_{T}=-\Phi_{x}(x_{T}),\end{cases}$ (2.3)

where $H_{x}$, $H_{y}$, $H_{z}$ are derivatives of $H$ with respect to $x$,
$y$, $z$, respectively. And the above system is essentially a special kind of
fully coupled FBSDEs.

Set

$w=\begin{pmatrix}x\\\ y\\\
z\end{pmatrix}\in\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{n\times
d},\qquad A(t,w)=\begin{pmatrix}-H_{x}\\\ H_{y}\\\ H_{z}\end{pmatrix}(t,w),$

and

$\left\langle w^{1},w^{2}\right\rangle=\left\langle
x^{1},x^{2}\right\rangle+\left\langle y^{1},y^{2}\right\rangle+\left\langle
z^{1},z^{2}\right\rangle.$

######  Definition 2.1.

A triple of process
$(x(\cdot),y(\cdot),z(\cdot)):[0,T]\times\Omega\rightarrow\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{n\times
d}$ is called an adapted solution of (2.3), if
$(x(\cdot),y(\cdot),z(\cdot))\in
M^{2}(0,T;\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{n\times d})$,
and it satisfies (2.3).

######  Assumption 1.

For any
$w,w^{\prime}\in\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{n\times
d}$ and $x,x^{\prime}\in\mathbb{R}^{n}$,

  * (i)

there exists a constant $\beta_{1}>0$, such that

$\begin{array}[]{c}|A(t,w)-A(t,w^{\prime})|\leq\beta_{1}|w-w^{\prime}|,\vspace{1ex}\\\
\end{array}$

and

$|\Phi_{x}(x)-\Phi_{x}(x^{\prime})|\leq\beta_{1}|x-x^{\prime}|.$

  * (ii)

there exists a constant $\beta_{2}>0$, such that the following monotonic
conditions hold.

$\displaystyle\langle A(t,w)-A(t,w^{\prime}),w-w^{\prime}\rangle$
$\displaystyle\leq-\beta_{2}|w-w^{\prime}|^{2}$
$\displaystyle\langle-\Phi_{x}(x)+\Phi_{x}(x^{\prime}),x-x^{\prime}\rangle$
$\displaystyle\geq\beta_{2}|x-x^{\prime}|^{2}.$

As equation (2.3) is at the same time a fully coupled FBSDE, we recall the
following existence and uniqueness theorem in [3, 4].

######  Theorem 1 (Theorem 3.1 in [3]).

Let Assumption 1 hold. Then there exists a unique adapted solution
$(x(\cdot),y(\cdot),z(\cdot))$ for (2.3).

Recently, numerical algorithms for solving the BSDEs and FBSDEs with deep
learning method [24, 25, 30, 31] have been proposed and demonstrated
remarkable performance. The main idea is to reformulate the BSDE into a
stochastic optimal control problem, where the solution $z$ of the BSDE is
regarded as a control and approximated with deep neural network, and the
terminal error is taken as the cost functional. Other different numerical
algorithms have also been developed for solving the FBSDEs and the related
PDEs [26, 27, 28, 29], which also regard the solution of the FBSDE ($y$ or
$z$) as a control and approximate it with appropriate loss function.

###  2.2 A novel method to solve the stochastic Hamiltonian system

As noted in the previous sections, the stochastic Hamiltonian system (2.3) is
essentially a fully coupled FBSDE and can be solved through the methods for
solving the FBSDEs. In this paper, we develop a novel method for solving the
Hamiltonian system(2.3). Different from the above mentioned algorithms for
solving the BSDEs or FBSDEs, our main idea is to find the corresponding
stochastic optimal control problem of the stochastic Hamiltonian system, and
then directly apply the deep learning method to solve the control problem.

$\forall x,y,u\in\mathbb{R}^{n},z,v\in\mathbb{R}^{n\times d}$, set

$F(t,x,u,v,y,z)=\langle y,u\rangle+\langle z,v\rangle-H(t,x,y,z)$ (2.4)

and

$f(t,x,u,v)=\max_{y,z}F(t,x,u,v,y,z).$ (2.5)

Here $f(t,x,u,v)$ is the Legendre-Fenchel transform of $H(t,x,y,z)$ with
respect to $(y,z)$. Due to the differentiability and strict convexity of the
Hamiltonian $H$, $f$ is also differential and strictly convex with respect to
$u$ and $v$ [32].

Consider the following control system,

$\begin{cases}\mathrm{d}x_{t}=u_{t}\mathrm{d}t+v_{t}\mathrm{d}B_{t},\vspace{1ex}\\\
x_{0}=a,\end{cases}$ (2.6)

and the cost functional is given as

$J(u(\cdot),v(\cdot))=\mathbb{E}\left[\int_{0}^{T}f(t,x_{t},u_{t},v_{t})\mathrm{d}t+\Phi(x_{T})\right],$
(2.7)

where the controls $u(\cdot)$ and $v(\cdot)$ belong to
$M^{2}(0,T;\mathbb{R}^{n})$ and $M^{2}(0,T;\mathbb{R}^{n\times d})$,
respectively. The set of all admissible controls is denoted by
$\mathcal{U}_{ad}[0,T]$. Any
$(u^{*}(\cdot),v^{*}(\cdot))\in\mathcal{U}_{ad}[0,T]$ satisfying

$J(u^{*}(\cdot),v^{*}(\cdot))=\inf_{(u(\cdot),v(\cdot))\in\mathcal{U}_{ad}[0,T]}J(u(\cdot),v(\cdot))$
(2.8)

is called an optimal control. The corresponding state trajectory
$x^{*}(\cdot)$ is called an optimal trajectory and the corresponding triple
$(x^{*}(\cdot),u^{*}(\cdot),v^{*}(\cdot))$ is called an optimal triple.

In the following we prove that solving the stochastic Hamiltonian system (2.3)
is equivalent to solving the stochastic optimal control problem (2.6)-(2.7).

We need the following assumption.

######  Assumption 2.

$f(t,x,u,v)$ is continuously differentiable with respect to $x$, $u$, $v$, and

$\displaystyle|f_{x}(t,x,u,v)|$ $\displaystyle\leq
C(|x|+|u|+|v|+1),\vspace{1ex}$ $\displaystyle|f_{u}(t,x,u,v)|$
$\displaystyle\leq C(|x|+|u|+|v|+1),\vspace{1ex}$
$\displaystyle|f_{v}(t,x,u,v)|$ $\displaystyle\leq
C(|x|+|u|+|v|+1),\vspace{1ex}$ $\displaystyle|f(t,0,0,0)|$ $\displaystyle\leq
C$

for some given $C>0$.

######  Theorem 2.

Let $H$ be a given $C^{2}$ real function and strictly convex with respect to
$y$ and $z$. The derivatives of $H$ and $\Phi$ satisfy Assumption 1; $f$
satisfies Assumption 2. Suppose that
$(x^{*}(\cdot),u^{*}(\cdot),v^{*}(\cdot))$ is the optimal triple of the
optimal control problem (2.6)-(2.7). Then
$(x^{*}(\cdot),y^{*}(\cdot),z^{*}(\cdot))$ uniquely solves the Hamiltonian
system (2.3), where $y^{*}(\cdot)$ can be given as

$\displaystyle y_{t}^{*}$
$\displaystyle=\mathbb{E}\left[\int_{t}^{T}-f_{x}(s,x_{s}^{*},u^{*}_{s},v^{*}_{s})\mathrm{d}s-\Phi_{x}(x^{*}_{T})\Big{|}\mathcal{F}_{t}\right],$
(2.9)

and $z^{*}(\cdot)$ can be solved by this following BSDE

$\left\\{\begin{array}[]{l}-\mathrm{d}y_{t}^{*}=-f_{x}(t,x^{*}_{t},u^{*}_{t},v^{*}_{t})\mathrm{d}t-z_{t}^{*}\mathrm{d}B_{t},\vspace{1ex}\\\
y_{T}^{*}=-\Phi_{x}(x^{*}_{T}),\qquad t\in[0,T].\end{array}\right.$ (2.10)

###### Proof.

Set

$h(t,x,u,v,y,z)=\langle y,u\rangle+\langle z,v\rangle-f(t,x,u,v).$ (2.11)

Under our assumptions, for any optimal triple
$(x^{*}(\cdot),u^{*}(\cdot),v^{*}(\cdot))$ of the optimal control problem
(2.6)-(2.7), we have the following extended stochastic Hamiltonian system
through the well-known stochastic maximum principle (SMP in short) ( e.g.
Theorem 4.1 in [12]),

$\begin{cases}\mathrm{d}x_{t}^{*}=u_{t}^{*}\mathrm{d}t+v_{t}^{*}\mathrm{d}B_{t},\vspace{1ex}\\\
-\mathrm{d}y_{t}^{*}=h_{x}(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})\mathrm{d}t-z_{t}^{*}\mathrm{d}B_{t},\vspace{1ex}\\\
x_{0}^{*}=a,\qquad y_{T}^{*}=-\Phi_{x}(x^{*}_{T}),\qquad
t\in[0,T],\end{cases}$ (2.12)

and

$h(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})=\max_{u\in\mathbb{R}^{n},\\\
v\in\mathbb{R}^{n\times d}}h(t,x^{*}_{t},u,v,y^{*}_{t},z^{*}_{t}),\ a.e.\
t\in[0,T].$ (2.13)

The solution of the extended stochastic Hamiltonian system (2.12)-(2.13) is a
5-tuple $(x^{*}(\cdot),u^{*}(\cdot),v^{*}(\cdot),y^{*}(\cdot),z^{*}(\cdot))$.

By Theorem 12.2 in [32], we have the inverse Legendre-Fenchel transform of
(2.5):

$H(t,x,y,z)=\max_{u,v}h(t,x,u,v,y,z).$ (2.14)

Because $h(t,x,u,v,y,z)$ is strictly concave in $u,v$, it yields that the
maximum point $(u^{*}_{t},v^{*}_{t})$ of (2.13) is uniquely determined by
$(x^{*}_{t},y^{*}_{t},z^{*}_{t})$ due to the implicit function existence
theorem:

$\displaystyle u^{*}_{t}$
$\displaystyle=\bar{u}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}),$ (2.15)
$\displaystyle v^{*}_{t}$
$\displaystyle=\bar{v}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}),$

and $\bar{u},\bar{v}$ are differentiable functions. By (2.14),

$\displaystyle H(t,x^{*}_{t},y^{*}_{t},z^{*}_{t})$
$\displaystyle=h(t,x^{*}_{t},\bar{u}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}),\bar{v}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}),y^{*}_{t},z^{*}_{t}),$
(2.16)

which leads to

$\displaystyle
f(t,x^{*}_{t},\bar{u}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}),\bar{v}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}))$
(2.17) $\displaystyle=$ $\displaystyle\langle
y^{*}_{t},\bar{u}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t})\rangle+\langle
z^{*}_{t},\bar{v}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t})\rangle-H(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}).$

Thus, the derivatives of $f(t,x,u,v)$ with respect to $u,v$ are

$\displaystyle
f_{u}(t,x^{*}_{t},\bar{u}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}),\bar{v}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}))=y^{*}_{t},$
(2.18) $\displaystyle
f_{v}(t,x^{*}_{t},\bar{u}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}),\bar{v}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}))=z^{*}_{t}.$

It can be verified that

$\displaystyle H_{x}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t})=$ $\displaystyle-
f_{x}(t,x^{*}_{t},\bar{u}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}),\bar{v}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}))$
(2.19) $\displaystyle=$ $\displaystyle-f_{x}(t,x^{*}_{t},u^{*}_{t},v^{*}_{t})$
$\displaystyle=$ $\displaystyle
h_{x}(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t}),$

which implies that $(y^{*}(\cdot),z^{*}(\cdot))$ solves the BSDE (2.10).
Taking the conditional expectation in the backward equation of (2.10), we have
(2.9) hold.

Similarly, we have

$\displaystyle
H_{y}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t})=\bar{u}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t})=u^{*}_{t},$
(2.20) $\displaystyle
H_{z}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t})=\bar{v}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t})=v^{*}_{t}.$

Thus, the extended stochastic Hamiltonian system (2.12) is just the
Hamiltonian system (2.3) and $(x^{*}(\cdot),y^{*}(\cdot),z^{*}(\cdot))$ solves
(2.3). Because $H$ satisfies the monotonicity condition in Assumption 1, the
uniqueness of the solution is proved. ∎

The following proposition can help us to construct our algorithms in the next
section.

######  Proposition 3.

Under the same assumptions as in Theorem 2, we have

$\displaystyle J(u^{*}(\cdot),v^{*}(\cdot))$
$\displaystyle=\mathbb{E}\left[\int_{0}^{T}f(t,x^{*}_{t},u^{*}_{t},v^{*}_{t})\mathrm{d}t+\Phi(x_{T}^{*})\right]$
(2.21)
$\displaystyle=\mathbb{E}\left[\int_{0}^{T}F(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})dt+\Phi(x_{T}^{*})\right],$

and $(y^{*}(\cdot),z^{*}(\cdot))$ of (2.3) can also be obtained by solving the
following BSDE:

$\left\\{\begin{array}[]{l}-\mathrm{d}y_{t}^{*}=-F_{x}(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})\mathrm{d}t-z_{t}^{*}\mathrm{d}B_{t},\vspace{1ex}\\\
y_{T}^{*}=-\Phi_{x}(x^{*}_{T}),\qquad t\in[0,T],\end{array}\right.$ (2.22)

where

$F(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})=\max_{y,z}F(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y,z).$
(2.23)

Then $y^{*}(\cdot)$ can be expressed as

$\displaystyle y_{t}^{*}$
$\displaystyle=\mathbb{E}\left[\int_{t}^{T}-F_{x}(s,x_{s}^{*},u^{*}_{s},v^{*}_{s},y^{*}_{s},z^{*}_{s})\mathrm{d}s-\Phi_{x}(x^{*}_{T})\Big{|}\mathcal{F}_{t}\right].$
(2.24)

###### Proof.

By the definition of $H$ and the SMP,

$\displaystyle H(t,x^{*}_{t},y^{*}_{t},z^{*}_{t})$
$\displaystyle=\max_{u,v}h(t,x^{*}_{t},u,v,y^{*}_{t},z^{*}_{t}),$ (2.25)
$\displaystyle=\langle y^{*}_{t},u^{*}_{t}\rangle+\langle
z^{*}_{t},v^{*}_{t}\rangle-f(t,x^{*}_{t},u^{*}_{t},v^{*}_{t}).$

Then, we have

$f(t,x^{*}_{t},u^{*}_{t},v^{*}_{t})=\langle y^{*}_{t},u^{*}_{t}\rangle+\langle
z^{*}_{t},v^{*}_{t}\rangle-H(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}).$ (2.26)

The definition of $f(t,x,u,v)$ shows that

$\displaystyle f(t,x^{*}_{t},u^{*}_{t},v^{*}_{t})$
$\displaystyle=\max_{y,z}F(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y,z),$ (2.27)
$\displaystyle=\langle\bar{y},u^{*}_{t}\rangle+\langle\bar{z},v^{*}_{t}\rangle-H(t,x^{*}_{t},\bar{y},\bar{z}),$

where $(\bar{y},\bar{z})$ is the minimum point.

Notice that $F(t,x,u,v,y,z)$ is strictly concave with respect to $y,z$. If
$(y^{*}_{t},z^{*}_{t})\neq(\bar{y},\bar{z})$, then

$\displaystyle F(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},\bar{y},\bar{z})$
$\displaystyle>F(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})$ (2.28)
$\displaystyle=\langle y^{*}_{t},u^{*}_{t}\rangle+\langle
z^{*}_{t},v^{*}_{t}\rangle-H(t,x^{*}_{t},y^{*}_{t},z^{*}_{t})$
$\displaystyle=f(t,x^{*}_{t},u^{*}_{t},v^{*}_{t}),$

which contradicts the formula (2.27). So we have
$(y^{*}_{t},z^{*}_{t})=(\bar{y},\bar{z})$ and (2.21) holds.

Because the strict concavity of $F(t,x,u,v,y,z)$ with respect to $y,z$, we
have

$\displaystyle F_{y}(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})$
$\displaystyle=u^{*}_{t}-H_{y}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t})=0$ (2.29)
$\displaystyle F_{z}(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})$
$\displaystyle=v^{*}_{t}-H_{z}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t})=0.$

Similar to the proof of Theorem 2, according to the implicit function
existence theorem, we can easily check that

$F_{x}(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})=-H_{x}(t,x^{*}_{t},y^{*}_{t},z^{*}_{t}),$
(2.30)

holds, and then (2.22) holds. Taking the conditional expectation on (2.22), we
have (2.24). ∎

In Theorem 2 and Proposition 3, we choose the stochastic optimal control
problem (2.6)-(2.7), whose diffusion term $b$ and drift term $\sigma$ of the
state equation are simplely $u$ and $v$. In fact, to simplify the linear terms
of $y$ and $z$ in the Hamiltonian $H$, we can also choose other forms of $b$
and $\sigma$, such as $\alpha_{1}(x)+\alpha_{2}(x)u$ and
$\beta_{1}(x)+\beta_{2}(x)v$, which are linear with respect to $u$ and $v$. In
this cases, the transformations (2.5) and (2.14) also hold. We show an example
of this form in subsection 4.2 for the numerical results.

Besides, we can still solve the Hamiltonian system (2.3) even if the
coefficients $H_{x},H_{y},H_{z}$ do not satisfy the monotonic conditions
(Assumption 1 (ii)) in Theorem 2 and Proposition 3. For example, the articles
[8, 9] studied the solvability of FBSDEs under relatively loose conditions. In
this situation, as long as the optimal controls reach the optimal values, the
Hamiltonian system (2.3) can be solved, however the solution is not
necessarily unique.

##  3 Numerical method for solving Hamiltonian systems

In Section 2, we present the idea of the stochastic optimal control method to
solve the Hamiltonian system (2.3). According to Theorem 2, we only need to
find the optimal control triple $(x^{*}(\cdot),u^{*}(\cdot),v^{*}(\cdot))$ of
the stochastic control problem (2.6) and (2.7), then the solution
$(y^{*}(\cdot),z^{*}(\cdot))$ can be obtained by taking the conditional
expectation on the backward SDE of (2.10). Therefore, effective approximation
method should be used to obtain the optimal triple of (2.6)-(2.7), especially
for high dimensional cases.

Deep neural networks are usually used to approximate functions defined on
finite-dimensional space, and the approximation relies on the composition of
layers with simple functions. On the basis of the universal approximation
theorem [33, 34], the neural networks have shown to be an effective tool and
gained great successes in many practical applications. In this paper, inspired
by [35], we simulate the stochastic optimal control problem (2.6)-(2.7) from a
direct way with the deep neural network and develop two different numerical
algorithms suitable for different cases.

Let $\pi$ be a partition of the time interval,
$0=t_{0}<t_{1}<t_{2}<\cdots<t_{N-1}<t_{N}=T$ of $[0,T]$. Define $\Delta
t_{i}=t_{i+1}-t_{i}$ and $\Delta B_{t_{i}}=B_{t_{i+1}}-B_{t_{i}}$, where
$B_{t_{i}}\sim\mathcal{N}(0,t_{i})$, for $i=0,1,2,\cdots,N-1$. We also denote

$\delta=\sup\limits_{0\leq i\leq N-1}\Delta t_{i},$

which is small enough. Then the Euler-Maruyama scheme of the state equation
(2.6) can be written as

$\left\\{\begin{array}[]{l}x_{t_{i+1}}^{\pi}=x_{t_{i}}^{\pi}+u_{t_{i}}^{\pi}\Delta
t_{i}+v_{t_{i}}^{\pi}\Delta B_{t_{i}},\vspace{1ex}\\\
x_{0}=a,\end{array}\right.$ (3.1)

and the corresponding cost functional is given as

$J(u^{\pi}(\cdot),v^{\pi}(\cdot))=\dfrac{1}{M}\sum_{m=1}^{M}\Big{[}\sum_{i=0}^{N-1}f(t_{i},x_{t_{i}}^{\pi,m},u_{t_{i}}^{\pi,m},v_{t_{i}}^{\pi,m})\Delta
t_{i}+\Phi(x_{t_{N}}^{\pi,m})\Big{]},$ (3.2)

where $M$ represents the number of Monte Carlo samples.

We introduce a feedforward neural network
$\varphi^{\theta}:[0,T]\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ of the
form

$\varphi^{\theta}=\mathcal{A}_{\ell}\circ\sigma_{\ell-1}\circ\mathcal{A}_{\ell-1}\circ\cdots\circ\sigma_{1}\circ\mathcal{A}_{1},$
(3.3)

where

  * •

$\ell$ is a positive integer specifying the depth of the neural network,

  * •

$\mathcal{A}_{1},\cdots,\mathcal{A}_{\ell}$ are functions of the form

$\displaystyle\mathcal{A}_{1}$
$\displaystyle=w_{1}x+b_{1}\in\mathbb{R}^{d_{1}},$
$\displaystyle\mathcal{A}_{i}$
$\displaystyle=w_{i}\mathcal{A}_{i-1}+b_{i}\in\mathbb{R}^{d_{i}},\qquad\text{for
}2\leq i\leq\ell,$

the matrix weights $w_{i}$ and bias vector $b_{i}$ are trainable parameters,
$\theta$ is the whole trainable parameters $\theta=(w_{i},b_{i})_{1\leq
i\leq\ell}$, and $d_{i}$ is the number of neurons at layer $i$,

  * •

$\sigma_{\ell-1},\cdots,\sigma_{1}$ are the nonlinear activation functions,
such as the sigmoid, the rectified linear unit (ReLU), the exponential linear
unit (ELU), etc.

We approximate the controls $u,v$ with two different neural networks, which
can be represented with (3.3) and denoted as $\varphi^{\theta^{u}}_{u}$ and
$\varphi^{\theta^{v}}_{v}$, respectively:

$\begin{cases}u=\varphi^{\theta^{u}}_{u}(t,x)=\varphi_{u}(t,x;\theta^{u})=\mathcal{A}_{\ell_{u}}^{u}\circ\sigma_{\ell_{u}-1}^{u}\circ\mathcal{A}_{\ell_{u}-1}^{u}\circ\cdots\circ\sigma_{1}^{u}\circ\mathcal{A}_{1}^{u}(t,x)\vspace{1ex}\\\
v=\varphi^{\theta^{v}}_{v}(t,x)=\varphi_{v}(t,x;\theta^{v})=\mathcal{A}_{\ell_{v}}^{v}\circ\sigma_{\ell_{v}-1}^{v}\circ\mathcal{A}_{\ell_{v}-1}^{v}\circ\cdots\circ\sigma_{1}^{v}\circ\mathcal{A}_{1}^{v}(t,x).\end{cases}$
(3.4)

The two neural networks have the same input dimension but different output
dimensions. In this paper, we use the common parameters of the neural networks
for all the time points, i.e. a single network is developed for simulating
each of the control, and the time point $t$ is regarded as an input of the
neural network.

###  3.1 Case 1: the function $f(t,x,u,v)$ has an explicit form

When the function $f(t,x,u,v)$ defined as (2.5) has an explicit form, then the
discrete cost functional (3.2) can be approximated directly with

$J(u^{\pi}(\cdot),v^{\pi}(\cdot))=\dfrac{1}{M}\sum_{m=1}^{M}\Big{[}\sum_{i=0}^{N-1}f(t_{i},x_{t_{i}}^{\pi,m},u_{t_{i}}^{\pi,m},v_{t_{i}}^{\pi,m})\Delta
t_{i}+\Phi(x_{t_{N}}^{\pi,m})\Big{]},$ (3.5)

which is also the loss function we need to minimize in the whole neural
network, and $u_{t_{i}}^{\pi},v_{t_{i}}^{\pi}$ are the outputs of the two
neural networks at time $t_{i}$. Both the neural networks approximating
$u_{t_{i}}^{\pi},v_{t_{i}}^{\pi}$ contain one $(n+1)$-dim input layer, three
$(n+10)$-dim hidden layers. The network of $u_{t_{i}}^{\pi}$ has an $n$-dim
output layer and that of $v_{t_{i}}^{\pi}$ has a $n\times d$-dim output layer.
In order to simplify the representation, here we use $\theta$ to represent the
training parameters $(\theta^{u},\theta^{v})$ for both of the neural networks.

To minimize the loss function (3.5) and learn the optimal parameters, some
basic optimization algorithms, such as stochastic gradient descent (SGD),
AdaGrad, RMSProp, and Adam which are already implemented in TensorFlow can be
used. In this paper, the Adam method [36] is adopted as the optimizer.

Once we obtain the approximations of the optimal controls $u^{*}$ and $v^{*}$,
we can calculate the numerical solution $y^{*}_{t}$ by taking the conditional
expectation on the Backward SDE of (2.10), which can be approximated with
Monte Carlo simulation:

$y_{t_{i}}^{\pi}=\dfrac{1}{M}\sum\limits_{m=1}^{M}\Big{[}\sum\limits_{j=i}^{N-1}-f_{x}(t_{j},x_{t_{j}}^{\pi,m},u_{t_{j}}^{\pi,m},v_{t_{j}}^{\pi,m})\Delta
t_{j}-\Phi_{x}(x_{t_{N}}^{\pi,m})\Big{]}.$ (3.6)

We show the whole network architecture in Figure 1, where $h^{u}$ and $h^{v}$
represent respectively the hidden layers of the neural networks
$\varphi^{\theta^{u}}_{u}$ and $\varphi^{\theta^{v}}_{v}$. For each of the
neural network, the common parameters is used for all the time points, and the
time $t$ is taken as one of the inputs of the neural network.

Figure 1: The whole network architecture for case 1: $f(t,x,u,v)$ has an
explicit representation. The boxes in purple, green and orange represent
respectively the input layer, the hidden layers and the output layers of the
neural networks $\varphi^{\theta^{u}}_{u}$ and $\varphi^{\theta^{v}}_{v}$. The
data flow of the neural networks is represented with black arrows.

The pseudo-code is shown in Algorithm 1.

Algorithm 1 Numerical algorithm for solving the Hamiltonian system of case 1

1:The Brownian motion $\Delta B_{t_{i}}$, initial state $a$, and time $t_{i}$;

2:The output controls $u_{t_{i}}^{\pi}$, $v_{t_{i}}^{\pi}$ and
$y_{t_{i}}^{\pi}$.

3:for $l=0$ to $maxstep$ do

4: $x_{0}^{l,\pi,m}=a$, $loss=0$;

5: for $i=0$ to $N-1$ do

6: $u^{l,\pi,m}_{t_{i}}=\varphi_{u}(t_{i},x^{l,\pi,m}_{t_{i}};\theta^{l});$

7: $v^{l,\pi,m}_{t_{i}}=\varphi_{v}(t_{i},x^{l,\pi,m}_{t_{i}};\theta^{l});$

8: $x^{l,\pi,m}_{t_{i+1}}=x^{l,\pi,m}_{t_{i}}+u^{l,\pi,m}_{t_{i}}\Delta
t_{i}+v^{l,\pi,m}_{t_{i}}\Delta B_{t_{i}};$

9:
$loss=loss+f(t_{i},x_{t_{i}}^{l,\pi,m},u_{t_{i}}^{l,\pi,m},v_{t_{i}}^{l,\pi,m})\Delta
t_{i};$

10: end for

11:
$loss=\dfrac{1}{M}\sum\limits_{m=1}^{M}\left[loss+\Phi(x_{t_{N}}^{l,\pi,m});\right]$

12: $\theta^{l+1}=Adam(\theta^{l},\nabla loss);$

13:end for

14:$y_{t_{i}}^{\pi}=\dfrac{1}{M}\sum\limits_{m=1}^{M}\Big{[}\sum\limits_{j=i}^{N-1}-f_{x}(t_{j},x_{t_{j}}^{l,\pi,m},u_{t_{j}}^{l,\pi,m},v_{t_{j}}^{l,\pi,m})\Delta
t_{j}-\Phi_{x}(x_{t_{N}}^{l,\pi,m})\Big{]}.$

###  3.2 Case 2: the function $f(t,x,u,v)$ does not have an explicit form

For the cases where the function $f$ does not have an explicit form, we can
still solve the Hamiltonian system (2.3) by constructing a different neural
network architecture. For any given optimal triple
$(x^{*}(\cdot),u^{*}(\cdot),v^{*}(\cdot))$ of the optimal control problem
(2.7), we assume that the solution $(y^{*}(\cdot),z^{*}(\cdot))$ of (2.22)
satisfy

$y^{*}_{t}=Y(t,x^{*}_{t}),\qquad z^{*}_{t}=Z(t,x^{*}_{t}),\qquad\forall
t\in[0,T],$ (3.7)

for some functions $Y$ and $Z$. Then according to Proposition 3, we have

$J(u^{*}(\cdot),v^{*}(\cdot))=\mathbb{E}\displaystyle\left[\int_{0}^{T}F(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})dt+\Phi(x_{T}^{*})\right],$
(3.8)

and

$\mathrm{d}x_{t}^{*}=u^{*}_{t}\mathrm{d}t+v_{t}^{*}\mathrm{d}B_{t},\quad
x_{0}^{*}=a,$ (3.9)

where $y^{*}_{t},z^{*}_{t}$ satisfy

$F(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})=\max\limits_{y,z}F(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y,z),$
(3.10)

and $F$ is defined by (2.4). Because of the strict concavity and
differentiable properties of $F$ with respect to $y,z$, the constraint
condition (3.10) can be rewritten as

$\begin{cases}F_{y}(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})=0,\vspace{1ex}\\\
F_{z}(t,x^{*}_{t},u^{*}_{t},v^{*}_{t},y^{*}_{t},z^{*}_{t})=0.\end{cases}$
(3.11)

In this way, the Hamiltonian system (2.3) can be solved by solving the
stochastic optimal control

$J(u(\cdot),v(\cdot))=\mathbb{E}\displaystyle\left[\int_{0}^{T}F(t,x_{t},u_{t},v_{t},y_{t},z_{t})\mathrm{d}t+\Phi(x_{T})\right],$
(3.12)

with the state constraint

$\begin{cases}\mathrm{d}x_{t}=u_{t}\mathrm{d}t+v_{t}\mathrm{d}B_{t},\quad
x_{0}=a,\vspace{1ex}\\\
F_{y}(t,x_{t},u_{t},v_{t},y_{t},z_{t})=0,\vspace{1ex}\\\
F_{z}(t,x_{t},u_{t},v_{t},y_{t},z_{t})=0.\end{cases}$ (3.13)

Now we focus on solving (3.12)-(3.13) with a new neural network architecture.
Firstly, the Euler-Maruyama scheme (3.1) is used to obtain the discrete form
of the optimal control problem. In addition to the neural networks for
simulating the controls $u$ and $v$, we need to construct two more neural
networks for simulating the functions $Y$ and $Z$,

$\begin{cases}y=\varphi^{\theta^{y}}_{y}(t,x)=\varphi_{y}(t,x;\theta^{y})=\mathcal{A}_{\ell_{y}}^{y}\circ\sigma_{\ell_{y}-1}^{y}\circ\mathcal{A}_{\ell_{y}-1}^{y}\circ\cdots\circ\sigma_{1}^{y}\circ\mathcal{A}_{1}^{y}(t,x)\vspace{1ex}\\\
z=\varphi^{\theta^{z}}_{z}(t,x)=\varphi_{z}(t,x;\theta^{z})=\mathcal{A}_{\ell_{z}}^{z}\circ\sigma_{\ell_{z}-1}^{z}\circ\mathcal{A}_{\ell_{z}-1}^{z}\circ\cdots\circ\sigma_{1}^{z}\circ\mathcal{A}_{1}^{z}(t,x).\end{cases}$
(3.14)

We also use the common parameters on all the time points for each of the four
neural networks, and the inputs of each neural network are $(t,x)$. All of the
four neural networks contain one $(n+1)$-dim input layer and three
$(n+10)$-dim hidden layers. The dimensions of the output layer for each neural
network are different, that of $y$ and $u$ are $n$-dim, and that of $z$ and
$v$ are $(n\times d)$-dim. We still adopt Adam as the optimizer.

We denote $\theta=(\theta_{uv},\theta_{yz})$ as all the parameters of the
neural networks, where $\theta_{uv}$ are the parameters of the neural networks
$\varphi_{u}$ and $\varphi_{v}$ (for simulating $u$ and $v$), and
$\theta_{yz}$ are the parameters of the neural networks $\varphi_{y}$ and
$\varphi_{z}$ (for simulating $y$ and $z$).

Then the cost functional of the control problem is approximated by

$J(u^{\pi}(\cdot),v^{\pi}(\cdot))=\dfrac{1}{M}\sum_{m=1}^{M}\Big{[}\sum_{i=0}^{N-1}F(t_{i},x_{t_{i}}^{\pi,m},u_{t_{i}}^{\pi,m},v_{t_{i}}^{\pi,m},y_{t_{i}}^{\pi,m},z_{t_{i}}^{\pi,m})\Delta
t_{i}+\Phi(x_{t_{N}}^{\pi,m})\Big{]},$ (3.15)

which is the first loss function we need to minimize, and
$u_{t_{i}}^{\pi},v_{t_{i}}^{\pi},y_{t_{i}}^{\pi},z_{t_{i}}^{\pi}$ are the
outputs of the whole neural networks at time $t_{i}$. In addition, in order to
guarantee that the conditions (3.11) hold, we introduce the other cost
functional

$\displaystyle J(y^{\pi}(\cdot),z^{\pi}(\cdot))$
$\displaystyle:=\dfrac{1}{M}\sum_{m=1}^{M}\sum_{i=0}^{N-1}\Big{[}|F_{y}(t_{i},x_{t_{i}}^{\pi,m},u_{t_{i}}^{\pi,m},v_{t_{i}}^{\pi,m},y_{t_{i}}^{\pi,m},z_{t_{i}}^{\pi,m})|^{2}$
(3.16)
$\displaystyle\qquad\qquad+|F_{z}(t_{i},x_{t_{i}}^{\pi,m},u_{t_{i}}^{\pi,m},v_{t_{i}}^{\pi,m},y_{t_{i}}^{\pi,m},z_{t_{i}}^{\pi,m})|^{2}\Big{]}.$

which is at the same time the second loss function we need to minimize in the
neural networks.

The update of the neural network parameters is carried out as follows. Suppose
that we have finished the update at the iteration step $l$ and obtain the
parameters $\theta^{l}=(\theta_{uv}^{l},\theta_{yz}^{l})$. We first calculate
the values
$(x^{\pi}_{t_{i}},u^{\pi}_{t_{i}},v^{\pi}_{t_{i}},y^{\pi}_{t_{i}},z^{\pi}_{t_{i}})$
with the parameters $\theta^{l}$. Then the parameters $\theta_{uv}^{l}$ are
updated to $\theta_{uv}^{l+1}$ by using one step Adam optimization with the
first loss function (3.15). In the following, the parameters
$(\theta_{uv}^{l+1},\theta_{yz}^{l})$ are used to calculate the values
$(x^{\pi}_{t_{i}},u^{\pi}_{t_{i}},v^{\pi}_{t_{i}},y^{\pi}_{t_{i}},z^{\pi}_{t_{i}})$
in (3.16). Then the parameters $\theta_{yz}^{l}$ are updated with the second
loss function (3.16) by Adam optimization. In each iteration step, the update
of parameters $\theta_{yz}^{l}$ can be performed multiple times, for example
$\kappa$ times, to ensure the loss function (3.16) is enough small. And after
$\kappa$ times update, the parameters of the neural networks are denoted as
$\theta_{yz}^{l+1}$. Finally the solution $y$ is obtained by taking the
conditional expectation of the backward SDE of (2.22) which can be calculated
with Monte Carlo simulation:

$y_{t_{i}}^{\pi}=\dfrac{1}{M}\sum\limits_{m=1}^{M}\Big{[}\sum\limits_{j=i}^{N-1}-F_{x}(t_{j},x_{t_{j}}^{\pi,m},u_{t_{j}}^{\pi,m},v_{t_{j}}^{\pi,m},y_{t_{j}}^{\pi,m},z_{t_{j}}^{\pi,m})\Delta
t_{j}-\Phi_{x}(x_{t_{N}}^{\pi,m})\Big{]}.$ (3.17)

The pseudo-code is given in Algorithm 2.

In fact,we can also deal with the maximum condition (3.10) directly. in this
situation, we need to maximize the second objective functional in addition to
(3.15),

$J(y(\cdot),z(\cdot))=\mathbb{E}\left[\int_{0}^{T}F(t,x_{t},u_{t},v_{t},y_{t},z_{t})\right]\mathrm{d}t.$
(3.18)

and a similar scheme can be given. The advantage for using the conditions
(3.11) instead of (3.18) is that we can determine the influence of the
constraint conditions by the value of the cost functional (3.16), as the
optimal value of (3.16) should be 0.

##  4 Numerical results

In this section, we show the numerical results for solving the Hamiltonian
system with our proposed algorithms. If not specifically mentioned, we use
6-layer fully connected neural networks for the approximation in these
examples, the number of time divisions is set to be $N=25$, and we mainly use
a piecewise constant decay learning rate which decreases from $3\times
10^{-3}$ to $1\times 10^{-3}$ with the increase of the number of iteration
steps. We adopt ELU as the activation function. In order to show the
performance of the proposed algorithms, we compare the results among the two
proposed algorithms and the Deep FBSDE method (briefly noted as DFBSDE in the
figures and tables of this section) which was developed as Algorithm 1 in our
previous work [31]. But different from [31], for better comparison with the
novel proposed algorithms, here we use the ELU activation function and remove
the batch normalization layer in the Deep FBSDE method. For each algorithm of
the examples, we perform ten independent runs to show more accurate results.

Algorithm 2 Numerical algorithm for solving the Hamiltonian system of case 2

1:The Brownian motion $\Delta B_{t_{i}}$, initial state $a$, and time $t_{i}$;

2:The processes
$(x^{l,\pi}_{t_{i}},u^{l,\pi}_{t_{i}},v^{l,\pi}_{t_{i}},y^{l,\pi}_{t_{i}},z^{l,\pi}_{t_{i}})$.

3:for $l=0$ to $maxstep$ do

4: for $k=0$ to $\kappa+1$ do

5: $x_{0}^{l,\pi,m}=a$, $loss_{1}=0$, $loss_{2}=0$;

6: for $i=0$ to $N-1$ do

7:
$u^{l,\pi,m}_{t_{i}}=\varphi_{u}(t_{i},x^{l,\pi,m}_{t_{i}};\theta_{uv}^{l});$

8:
$v^{l,\pi,m}_{t_{i}}=\varphi_{v}(t_{i},x^{l,\pi,m}_{t_{i}};\theta_{uv}^{l});$

9:
$y^{l,\pi,m}_{t_{i}}=\varphi_{y}(t_{i},x^{l,\pi,m}_{t_{i}};\theta_{yz}^{l});$

10:
$z^{l,\pi,m}_{t_{i}}=\varphi_{z}(t_{i},x^{l,\pi,m}_{t_{i}};\theta_{yz}^{l});$

11: $x^{l,\pi,m}_{t_{i+1}}=x^{l,\pi,m}_{t_{i}}+u^{l,\pi,m}_{t_{i}}\Delta
t_{i}+v^{l,\pi,m}_{t_{i}}\Delta B_{t_{i}};$

12:
$loss_{1}=loss_{1}+F(t_{i},x_{t_{i}}^{l,\pi,m},u_{t_{i}}^{l,\pi,m},v_{t_{i}}^{l,\pi,m},y_{t_{i}}^{l,\pi,m},z_{t_{i}}^{l,\pi,m})\Delta
t_{i};$

13:
$loss_{2}=loss_{2}+\dfrac{1}{M}\sum\limits_{m=1}^{M}\left[|F_{y}(t_{i},x_{t_{i}}^{l,\pi,m},u_{t_{i}}^{l,\pi,m},v_{t_{i}}^{l,\pi,m},y_{t_{i}}^{l,\pi,m},z_{t_{i}}^{l,\pi,m})|^{2}\right];$

14:
$loss_{2}=loss_{2}+\dfrac{1}{M}\sum\limits_{m=1}^{M}\left[|F_{z}(t_{i},x_{t_{i}}^{l,\pi,m},u_{t_{i}}^{l,\pi,m},v_{t_{i}}^{l,\pi,m},y_{t_{i}}^{l,\pi,m},z_{t_{i}}^{l,\pi,m})|^{2}\right];$

15: end for

16:
$loss_{1}=\dfrac{1}{M}\sum\limits_{m=1}^{M}\left[loss_{1}+\Phi(x_{t_{N}}^{l,\pi,m})\right];$

17: if $k=0$ then

18: $\theta_{uv}^{l}=Adam(\theta_{uv}^{l},\nabla loss_{1});$

19: $\theta^{l}=(\theta_{uv}^{l},\theta_{yz}^{l});$

20: else

21: $\theta_{yz}^{l}=Adam(\theta_{yz}^{l},\nabla loss_{2});$

22: $\theta^{l}=(\theta_{uv}^{l},\theta_{yz}^{l});$

23: end if

24: end for

25: $\theta^{l+1}=(\theta_{uv}^{l},\theta_{yz}^{l});$

26:end for

27:$y_{t_{i}}^{\pi}=\dfrac{1}{M}\sum\limits_{m=1}^{M}\Big{[}\sum\limits_{j=i}^{N-1}-F_{x}(t_{j},x_{t_{j}}^{l,\pi,m},u_{t_{j}}^{l,\pi,m},v_{t_{j}}^{l,\pi,m},y_{t_{j}}^{l,\pi,m},z_{t_{j}}^{l,\pi,m})\Delta
t_{j}-\Phi_{x}(x_{t_{N}}^{l,\pi,m})\Big{]}.$

###  4.1 Example 1: a linear Hamiltonian system

Firstly, we consider the following linear quadratic Hamiltonian

$H(t,x,y,z)=\langle x,y\rangle+\dfrac{1}{4}\langle y,y\rangle+\langle
z,z\rangle,\qquad\Phi(x)=\langle\dfrac{1}{2}Qx,x\rangle,$ (4.1)

where $(x,y,z)\in\mathbb{R}^{n+n+n}$ and $Q$ is a given matrix valued in
$\mathbb{R}^{n\times n}$, and the corresponding Hamiltonian system is given as

$\left\\{\begin{array}[]{l}\mathrm{d}x_{t}=(x_{t}+\dfrac{1}{2}y_{t})\mathrm{d}t+2z_{t}\mathrm{d}B_{t},\vspace{1ex}\\\
-\mathrm{d}y_{t}=y_{t}\mathrm{d}t-z_{t}\mathrm{d}B_{t},\vspace{1ex}\\\
x_{0}=a,\qquad y_{T}=-\Phi_{x}(x_{T})=-Qx_{T},\end{array}\right.$ (4.2)

which is a linear FBSDE and $B$ is a $1$-dimensional standard Brownian motion.
It can be easily check that this linear FBSDE has a unique solution [3, 4,
37].

As is known, the linear FBSDE connects with a Riccati equation. Suppose the
solution of FBSDE (4.2) is in the following form:

$\displaystyle y_{t}=-K_{t}x_{t},\qquad z_{t}=-M_{t}x_{t}.$

Combing it with (4.2), we then obtain a Riccati equation

$\begin{cases}\dot{K}_{t}-\frac{1}{2}K_{t}^{2}+2K_{t}=0,\\\ M_{t}=0,\qquad
K_{T}=Q,\end{cases}$ (4.3)

where $K_{t}$ is a matrix function and $\dot{K}_{t}$ is the derivative of
$K_{t}$ with respect to $t$. And the solution of (4.3) can be approximated
with the ODE45 method in MATLAB (ODE45 in short), which solves determined
ordinary differential equations with the four-order Runge-Kutta method. In
order to show the performance of our novel proposed algorithms, we take the
numerical solution of the Ricatti equation (4.3) with ODE45 as a benchmark.

Now we consider the corresponding optimal control problem,

$\left\\{\begin{array}[]{l}\mathrm{d}x_{t}=u_{t}\mathrm{d}t+v_{t}\mathrm{d}B_{t},\vspace{1ex}\\\
x_{0}=a,\end{array}\right.$

with cost functional

$J(u(\cdot),v(\cdot))=\mathbb{E}\Big{\\{}\int_{0}^{T}f(t,x_{t},u_{t},v_{t})+\Phi(x_{T})\Big{\\}},$

where

$\displaystyle f(t,x,u,v)$
$\displaystyle=\max_{y,z}F(t,x,u,v,y,z)\vspace{1ex}$
$\displaystyle=|u-x|^{2}+\dfrac{1}{4}|v|^{2},$

and

$F(t,x,u,v,y,z)=\langle y,u\rangle+\langle z,v\rangle-\langle
x,y\rangle-\dfrac{1}{4}\langle y,y\rangle-\langle z,z\rangle.$

In this example, we set $a=(1.0,\cdots,1.0)\in\mathbb{R}^{n},T=0.1$ and give
the form of $Q$ as

$\begin{bmatrix}1&\lambda&\lambda&\cdots&\lambda\\\
\lambda&1&\lambda&\cdots&\lambda\\\ \lambda&\lambda&1&\cdots&\lambda\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\lambda&\lambda&\lambda&\cdots&1\end{bmatrix}$

where $\lambda$ is a given constant between 0 and 1. For example, if
$\lambda=0.0$, $Q=E_{n}$ is a $n$-order unit matrix, The numerical result of
(4.3) with ODE45 can be solved with $K_{0}=1.1573E_{n}$ for $n=100$, then the
value of $y_{0}$ can be obtained by

$y_{0}=-K_{0}x_{0}=-1.1573a.$

which is taken as the benchmark results in this example.

Even though the function $f$ can be solved explicitly in this example, we
calculate the numerical results through both the two proposed algorithms and
regard that $f$ does not have an explicit form in Algorithm 2. The comparison
results on the approximated solution $y_{0}$ between our proposed stochastic
control methods (Algorithm 1 and 2) and the Deep FBSDE method are shown in
Table 1, and the solution with ODE45 is regarded as the benchmark. Note that
when the initial state $x=a$ , the solution of $y_{0}$ is a vector and all its
elements are equal, thus we show the value of the first element of $y_{0}$ in
Table 1. Moreover, we show the relative errors between our approximated
solution of $y_{0}$ and that of ODE45, and compute the variances among ten
independent runs of the approximated solution $y_{0}$. We also change the
parameter $\lambda$, and study the corresponding approximation results.

Table 1: The implementations of different terminal $Q$ with $n=100$ | Riccati | DEEP FBSDE | Alg 1 | Alg 2
---|---|---|---|---
Mean | Rela. Error | Var. | Mean | Rela. Error | Var. | Mean | Rela. Error | Var.
$\lambda=0.0$ | -1.1573 | -1.15733 | 2.907e-05 | 6.209e-08 | -1.15751 | 1.797e-04 | 3.141e-07 | -1.15651 | 6.839e-04 | 3.065e-06
$\lambda=0.2$ | -11.8093 | -11.6314 | 1.506e-02 | 2.392e-01 | -11.8113 | 1.733e-04 | 2.441e-03 | -11.8222 | 1.095e-03 | 5.997e-03
$\lambda=0.4$ | -15.2711 | -14.1265 | 7.495e-02 | 9.101e-01 | -15.3120 | 2.681e-03 | 1.739e-03 | -15.2264 | 2.930e-03 | 1.622e-01
$\lambda=0.6$ | -16.9860 | -10.4284 | 3.861e-01 | 2.753e+01 | -17.0461 | 3.859e-03 | 3.563e-03 | -17.0516 | 3.860e-03 | 4.633e-02
$\lambda=0.8$ | -18.0101 | -9.2087 | 4.887e-01 | 2.279e+01 | -17.9844 | 1.426e-03 | 2.528e-02 | -18.1494 | 7.733e-03 | 2.726e-02
$\lambda=1.0$ | -18.6920 | -9.9503 | 4.677e-01 | 2.300e+01 | -18.7162 | 1.297e-03 | 2.000e-02 | -18.8837 | 1.026e-02 | 8.245e-02

From Table 1, we can see that the novel proposed algorithms show much more
stable performance than the Deep FBSDE method. For different terminals with
different parameters $\lambda$, the novel proposed algorithms demonstrate more
stable relative errors and variances. The Deep FBSDE method perform well when
the terminal is a unit matrix ($\lambda=0.0$), but when we change the terminal
to other forms ($\lambda\not=0.0$), the results of the Deep FBSDE method
diverge. As we know, the learning rate is one of the important factors
affecting the approximation results. When we choose smaller learning rate, the
Deep FBSDE method can converge, but it need much more iteration steps than our
proposed algorithms. As an example for the case $\lambda\not=0.0$ with smaller
learning rate, we show the approximation results in the following Figure 3 for
$\lambda=0.8$.

In Figure 2, we show the curves and the variances of the approximated results
$y_{0}$ with different iteration steps for $\lambda=0.0$ and $\lambda=0.8$
respectively, and the black lines represent the results with ODE45 which is
taken as the benchmark. The upper two figures in Figure 2 exhibit the results
for $\lambda=0.0$. From the upper left figure, we can see that when the number
of iteration steps is close to 10000, the approximated solution $y_{0}$ with
Algorithm 1, 2 and the Deep FBSDE method are all very close to the results
with ODE45. Moreover, that of Algorithm 1 and 2 have smaller variation scopes
among ten independent runs and converge within less iteration steps than that
of the Deep FBSDE method. The upper right figure shows that when the number of
iteration steps tends to be 10000, the variance curve of $y_{0}$ with
Algorithm 1 and 2 are also close to that of the Deep FBSDE method. The lower
two figures in Figure 2 exhibit the results for $\lambda=0.8$. We can see that
when the number of iteration steps tends to be 10000, the approximation
results of Algorithm 1 and 2 are close to the benchmark. However, that of the
Deep FBSDE method is far from the benchmark, and the variation scope and
variance increase with the increase of iteration steps.

Figure 2: Approximation results with a piecewise decay learning rate from
$3\times 10^{-3}$ to $1\times 10^{-3}$ for $\lambda=0.0$ in the upper figures
and $\lambda=0.8$ in the lower figures. The left figures show the mean and
variation scopes of the approximated solution $y_{0}$ among 10 independent
runs, and the right figures exhibit the variance curves of $y_{0}$ among 10
independent runs. The black lines in the left figure represent the results
with ODE45 which are taken as the benchmarks. We can see that our novel
proposed algorithms(Algorithms 1 and 2) have more stable convergence and
converge within less iteration steps.

Figure 3: Approximation results with a constant learning rate of $1\times
10^{-3}$. We can see from the left figure that comparing with the novel
proposed algorithms(Algorithms 1 and 2), the Deep FBSDE method needs more
iteration steps to achieve stable convergence, and it has a much smaller
variance at the end of the training from the right figure.

###  4.2 Example 2: a nonlinear Hamiltonian system

Given the Hamiltonian $H$ as

$H(t,x,y,z)=\dfrac{1}{2}\langle y,y\circ\cos^{2}x\rangle+\dfrac{1}{2}\langle
z,z\circ\sin^{2}x\rangle+\langle y,\cos x\rangle+\langle z,\sin
x\rangle-\dfrac{1}{2}\langle x,x\rangle$ (4.4)

and

$\Phi(x)=\dfrac{1}{2}\langle x,x\rangle,$ (4.5)

where $(x,y,z)\in\mathbb{R}^{n+n+n}$. Here $\circ$ represents the Hadamard
product,

$x\circ
y=(x_{1}y_{1},x_{2}y_{2},\cdots,x_{n}y_{n})\in\mathbb{R}^{n},\qquad\forall
x,y\in\mathbb{R}^{n}$

The corresponding Hamiltonian system is

$\left\\{\begin{array}[]{l}\mathrm{d}x_{t}=\cos x_{t}\circ(y_{t}\circ\cos
x_{t}+1)\mathrm{d}t+\sin x_{t}\circ(z_{t}\circ\sin
x_{t}+1)\circ\mathrm{d}B_{t},\vspace{1ex}\\\
-\mathrm{d}y_{t}=\Big{[}-y_{t}\circ\sin x_{t}\circ(y_{t}\circ\cos
x_{t}+1)+z_{t}\circ\cos x_{t}\circ(z_{t}\circ\sin
x_{t}+1)-x_{t}\Big{]}\mathrm{d}t\\\ \qquad\qquad-
z_{t}\circ\mathrm{d}B_{t},\vspace{1ex}\\\ x_{0}=a,\qquad
y_{T}=-\Phi_{x}(x_{T}),\end{array}\right.$ (4.6)

where $B$ is a $n$-dimensional Brownian motion.

Here we introduce a stochastic optimal control problem which is different from
(2.6):

$\left\\{\begin{array}[]{l}\mathrm{d}x_{t}=\cos
x_{t}\circ(u_{t}+1)\mathrm{d}t+\sin
x_{t}\circ(v_{t}+1)\circ\mathrm{d}B_{t}\vspace{1ex}\\\
x_{0}=a,\end{array}\right.$ (4.7)

with the cost functional

$J(u(\cdot),v(\cdot))=\mathbb{E}\left\\{\int_{0}^{T}f(t,x_{t},u_{t},v_{t})dt+\Phi(x_{T})\right\\},$

where $f(t,x,u,v)$ is given as

$\displaystyle f(t,x,u,v)$
$\displaystyle=\max_{y,z}F(t,x,u,v,y,z)\vspace{1ex}$
$\displaystyle=\max_{y,z}\left\\{\langle y,\cos x\circ(u+1)\rangle+\langle
z,\sin x\circ(v+1)\rangle-H(t,x,y,z)\right\\}.$

Then $f(t,x,u,v)$ can be solved as

$f(t,x,u,v)=\dfrac{1}{2}\big{[}\langle x,x\rangle+\langle u,u\rangle+\langle
v,v\rangle\big{]},$

and

$\displaystyle y$ $\displaystyle=u\circ\dfrac{1}{\cos x},$ $\displaystyle z$
$\displaystyle=v\circ\dfrac{1}{\sin x}.$

Define $h(t,x,y,v,y,z)$ as

$\displaystyle h(t,x,u,v,y,z)=\langle y,\cos x\circ(u+1)\rangle+\langle z,\sin
x\circ(v+1)\rangle-f(t,x,u,v),$

then we get the value of $y_{0}$ by

$\displaystyle y_{0}$
$\displaystyle=\mathbb{E}\left\\{\int_{0}^{T}h_{x}(t,x_{t},u_{t},v_{t},y_{t},z_{t})dt-\Phi_{x}(x_{T})\right\\}$
$\displaystyle=\mathbb{E}\left\\{\int_{0}^{T}(-u_{t}\circ\tan
x_{t}\circ(u_{t}+1)+v_{t}\circ\cot
x_{t}\circ(v_{t}+1)-x_{t})-\Phi_{x}(x_{T})\right\\}.$

We set $a=(1.0,\cdots,1.0)\in\mathbb{R}^{n}$ and $T=0.1$. The comparison
results among Algorithm 1, Algorithm 2 and the Deep FBSDE method are shown in
Figure 4. We can see from the left figure that when the number of iteration
steps is 5000, The approximated value of $Y_{0}$ with Algorithm 1 and the Deep
FBSDE method are very close, which are $-1.0835$ and $-1.0834$, respectively.
And the approximated solution with Algorithm 2 is $-1.1208$, which is slightly
smaller than that with Algorithm 1 and the Deep FBSDE method. Similar to the
first example, the variation scope of $y_{0}$ for Algorithm 1, 2 are much
smaller than that of the Deep FBSDE method. Besides, at the end of the
training, the variances of $y_{0}$ for Algorithm 1, 2 and the Deep FBSDE
method are very close.

Figure 4: We can see from the figure that when the number of iteration steps
tends to be 5000, the approximated values of $y_{0}$ for all the three methods
are close. However, the Deep FBSDE method shows larger variation scopes than
Algorithm 1 and 2 during the training. And the variances of these three
methods are also close at the end of the training.

###  4.3 Example 3: a Hamiltonian system with exponents in the drift term

In the third example, we solve a Hamiltonian system with exponents in the
drift term. Consider the following Hamiltonian

$H(t,x,y,z)=\log\left(\sum_{i=1}^{n}\exp(y_{i})\right)+\dfrac{1}{2}\langle
y,y\rangle+\langle z,z\rangle+\langle z,x\rangle+\dfrac{1}{5}\langle
x,x\rangle,$ (4.8)

where $(x,y,z)\in\mathbb{R}^{n+n+n}$ and $y=(y_{1},\cdots,y_{n})$. The
terminal function is given as $\Phi(x)=\dfrac{1}{2}\langle x,x\rangle$. Then
the Hamiltonian system we need to solve is given as following,

$\left\\{\begin{array}[]{l}\mathrm{d}x_{t}=\left[\exp(y_{t})\left(\sum_{i=1}^{n}\exp(y_{it})\right)^{-1}+y_{t}\right]\mathrm{d}t+(x_{t}+2z_{t})\circ\mathrm{d}B_{t},\vspace{1ex}\\\
-\mathrm{d}y_{t}=(\dfrac{2}{5}x_{t}+z_{t})\mathrm{d}t-z_{t}\circ\mathrm{d}B_{t},\vspace{1ex}\\\
x_{0}=a,\qquad y_{T}=-\Phi_{x}(x_{T}),\end{array}\right.$ (4.9)

where $\exp(y_{t})=(\exp(y_{it}),\cdots,\exp(y_{nt}))$ and $B$ is a
$n$-dimensional Brownian motion.

The corresponding stochastic optimal control problem is

$\begin{cases}\mathrm{d}x_{t}=u_{t}\mathrm{d}t+v_{t}\circ\mathrm{d}B_{t},\vspace{1ex}\\\
x_{0}=a,\end{cases}$

with the cost functional

$\displaystyle
J(u(\cdot),v(\cdot))=\mathbb{E}\left\\{\int_{0}^{T}f(t,x_{t},u_{t},v_{t})\mathrm{d}t+\Phi(x_{T})\right\\},$

where

$\displaystyle f(t,x,u,v)$ $\displaystyle=\max_{y,z}F(t,x,u,v,y,z),$ (4.10)
$\displaystyle=\max_{y,z}\left\\{\langle y,u\rangle+\langle
z,v\rangle-H(t,x,y,z)\right\\}.$

Different from the previous two examples, in this example, the function $f$
defined as (4.10) does not have an explicit representation. In this situation,
Algorithm 1 is not applicable, thus we mainly make the comparison between the
results of Algorithm 2 and the Deep FBSDE method.

We set $T=0.2$, $a=(0.5,\cdots,0.5)\in\mathbb{R}^{n}$ and $n=100$ in this
example. The comparison results between Algorithm 2 and the Deep FBSDE method
are shown in Figure 5. We can see that at the end of the training, the
approximated solution $y_{0}$ of the two methods are close, and both the
variances are small enough. Similar with the previous two examples, Algorithm
2 shows smaller variation scope and converges within less iteration steps than
the Deep FBSDE method.

Figure 5: We can see from the above figure that at the end of training, the
approximated solutions of $y_{0}$ with both Algorithm 2 and the Deep FBSDE
method are very close, and the mean values of $y_{0}$ among 10 independent
runs are $-0.41297$ for Algorithm 2 and $-0.41211$ for the Deep FBSDE method.
The variance of Algorithm 2 is slightly larger than that of the Deep FBSDE
method, but it converges within less iteration steps.

In the implementations of Algorithm 2 for all the three examples, we minimize
the norm of the derivatives of the function $F$ with respect to $y$ and $z$
according to (3.16) so that they are equal to $0$. As an alternative, we can
also maximize the cost functional defined as (3.18) in the implementations.

##  5 Conclusion

In this paper, different from the general way of solving the FBSDEs, we
propose a novel method which solve the Hamiltonian system from the view of the
stochastic optimal control via deep learning. Two different algorithms
suitable for different cases are developed. From the numerical results, the
novel proposed Algorithms 1 and 2 demonstrate faster convergence rate and more
stable performance than the Deep FBSDE method. In some cases, Algorithm 1 and
2 shows higher accuracy.

## References

  * [1] R. Ortega, A. J. Van Der Schaft, I. Mareels, and B. Maschke, “Putting energy back in control,” IEEE Control Systems Magazine, vol. 21, no. 2, pp. 18–33, 2001.
  * [2] R. Ortega, A. van der Schaft, B. Maschke, and G. Escobar, “Interconnection and damping assignment passivity-based control of port-controlled hamiltonian systems,” Automatica, vol. 38, no. 4, pp. 585–596, 2002.
  * [3] Y. Hu and S. Peng, “Solution of forward-backward stochastic differential equations,” Probability Theory and Related Fields, vol. 103, no. 2, pp. 273–283, 1995.
  * [4] S. Peng and Z. Wu, “Fully coupled forward-backward stochastic differential equations and applications to optimal control,” Siam Journal on Control and Optimization, vol. 37, no. 3, pp. 825–843, 1999.
  * [5] S. Peng, “Problem of eigenvalues of stochastic hamiltonian systems with boundary conditions,” Stochastic Processes and Their Applications, vol. 88, no. 2, pp. 259–290, 2000.
  * [6] Y. Hu and J. Yong, “Forward–backward stochastic differential equations with nonsmooth coefficients,” Stochastic Processes and their Applications, vol. 87, no. 1, pp. 93–106, 2000.
  * [7] J. Ma, P. Protter, and J. Yong, “Solving forward-backward stochastic differential equations explicitly — a four step scheme,” Probability Theory and Related Fields, vol. 98, no. 3, pp. 339–359, 1994.
  * [8] J. Ma and J. Yong, “Solvability of forward-backward sdes and the nodal set of hamilton-jacobi-bellman equations,” Chinese Annals of Mathematics, 1993\.
  * [9] J. Ma and J. Yong, “Approximate solvability of forward—backward stochastic differential equations,” Applied Mathematics and Optimization, vol. 45, no. 1, pp. 1–22, 2002.
  * [10] N. E. Karoui, S. Peng, and M. C. Quenez, “Backward Stochastic Differential Equations in Finance,” Mathematical Finance, vol. 7, no. 1, pp. 1–71, 1997.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/1467-9965.00022.
  * [11] E. Pardoux and S. Tang, “Forward-backward stochastic differential equations and quasilinear parabolic PDEs,” Probability Theory and Related Fields, vol. 114, pp. 123–150, may 1999.
  * [12] A. Bensoussan, “Lectures on stochastic control,” in Nonlinear filtering and stochastic control, pp. 1–62, Springer, 1982.
  * [13] J. Yong and X. Zhou, Stochastic Controls-Hamiltonian System and HJB Equations.  Springer, 1999.
  * [14] E. Tadmor, “A review of numerical methods for nonlinear partial differential equations,” Bulletin of the American Mathematical Society, vol. 49, no. 4, pp. 507–554, 2012.
  * [15] B. Bouchard and N. Touzi, “Discrete-time approximation and monte-carlo simulation of backward stochastic differential equations,” Stochastic Processes & Their Applications, vol. 111, no. 2, pp. 175–206, 2004.
  * [16] V. Bally and G. Pages, “A quantization algorithm for solving multidimensional discrete-time optimal stopping problems,” Bernoulli, vol. 9, no. 6, pp. 1003–1049, 2003.
  * [17] M. A. Jin, S. Jie, and Y. Zhao, “On numerical approximations of forward-backward stochastic differential equations,” Siam Journal on Numerical Analysis, vol. 46, no. 5, pp. 2636–2661, 2008.
  * [18] C. Bender and J. Zhang, “Time discretization and markovian iteration for coupled fbsdes,” Annals of Applied Probability, vol. 18, no. 1, pp. 143–177, 2008.
  * [19] F. Yu, W. Zhao, and T. Zhou, “Multistep schemes for forward backward stochastic differential equations with jumps,” Journal of Scientific Computing, vol. 69, no. 2, pp. 1–22, 2016.
  * [20] E. Weinan, M. Hutzenthaler, A. Jentzen, and T. Kruse, “On multilevel picard numerical approximations for high-dimensional nonlinear parabolic partial differential equations and high-dimensional nonlinear backward stochastic differential equations,” arXiv: Numerical Analysis, 2017.
  * [21] G. N. Milstein and M. V. Tretyakov, “Numerical algorithms for forward-backward stochastic differential equations connected with semilinear parabolic equations,” vol. 28, pp. 561–582, 2004.
  * [22] F. Yu, W. Zhao, and Z. Tao, “Efficient spectral sparse grid approximations for solving multi-dimensional forward backward sdes,” Discrete and Continuous Dynamical Systems - Series B, vol. 22, no. 9, 2016.
  * [23] T. P. Huijskens, M. Ruijter, and C. W. Oosterlee, “Efficient numerical fourier methods for coupled forward-backward sdes,” Journal of Computational and Applied Mathematics, vol. 296, pp. 593–612, 2016.
  * [24] E. Weinan, J. Han, and A. Jentzen, “Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations,” Communications in Mathematics & Statistics, vol. 5, no. 4, pp. 349–380, 2017.
  * [25] J. Han, A. Jentzen, and E. Weinan, “Solving high-dimensional partial differential equations using deep learning,” Proceedings of the National Academy of Sciences of the United States of America, vol. 115, no. 34, pp. 8505–8510, 2018.
  * [26] C. Huré, H. Pham, and X. Warin, “Deep backward schemes for high-dimensional nonlinear PDEs,” arXiv:1902.01599 [cs, math, stat], June 2020.
  * [27] H. Pham, X. Warin, and M. Germain, “Neural networks-based backward scheme for fully nonlinear PDEs,” SN Partial Differential Equations and Applications, vol. 2, p. 16, Feb. 2021.
  * [28] M. Raissi, “Forward-backward stochastic neural networks: Deep learning of high-dimensional partial differential equations,” arXiv: 1804.07010, 2018\.
  * [29] H. Wang, H. Chen, A. Sudjianto, R. Liu, and Q. Shen, “Deep learning-based bsde solver for libor market model with application to bermudan swaption pricing and hedging,” 2018.
  * [30] J. Han and J. Long, “Convergence of the deep bsde method for coupled fbsdes,” arXiv:1811.01165, 2018.
  * [31] S. Ji, S. Peng, Y. Peng, and X. Zhang, “Three algorithms for solving high-dimensional fully coupled fbsdes through deep learning,” IEEE Intelligent Systems, vol. 35, no. 3, pp. 71–84, 2020.
  * [32] R. T. ROCKAFELLAR, Convex Analysis.  Princeton University Press, 1970.
  * [33] G. Cybenko, “Approximation by superpositions of a sigmoidal function,” Mathematics of Control, Signals and Systems, vol. 2, no. 4, pp. 303–314, 1989\.
  * [34] K. Hornik, “Approximation capabilities of multilayer feedforward networks,” Neural Networks, vol. 4, no. 2, pp. 251–257, 1991.
  * [35] J. Han and W. E, “Deep learning approximation for stochastic control problems,” NIPS Workshop on Deep Reinforcement Learning, 2016.
  * [36] A. C. Wilson, R. Roelofs, R. Stern, N. Srebro, and B. Recht, “The marginal value of adaptive gradient methods in machine learning,” 31st Conference on Neural Information Processing Systems (NIPS 2017), 2017.
  * [37] J. Ma, Z. Wu, D. Zhang, and J. Zhang, “On well-posedness of forward-backward SDEs–A unified approach,” The Annals of Applied Probability, vol. 25, no. 4, pp. 2168 – 2214, 2015.