# An Online Algorithm for Solving Feedback Optimal Control Problems with
Partial Observations

Siming Liang , Ruoyu Hu, Feng Bao , Richard Archibald , and Guannan Zhang
Department of Mathematics, Florida State University, Tallahassee, Florida,
USADepartment of Mathematics, Florida State University, Tallahassee, Florida,
USADepartment of Mathematics, Florida State University, Tallahassee, Florida,
USAComputer Science and Mathematics Division, Oak Ridge National Laboratory,
Oak Ridge, Tennessee, USA Computer Science and Mathematics Division, Oak Ridge
National Laboratory, Oak Ridge, Tennessee, USA

###### Abstract

This paper presents a novel methodology to tackle feedback optimal control
problems in scenarios where the exact state of the controlled process is
unknown. It integrates data assimilation techniques and optimal control
solvers to manage partial observation of the state process, a common
occurrence in practical scenarios. Traditional stochastic optimal control
methods assume full state observation, which is often not feasible in real-
world applications. Our approach underscores the significance of utilizing
observational data to inform control policy design. Specifically, we introduce
a kernel learning backward stochastic differential equation (SDE) filter to
enhance data assimilation efficiency and propose a sample-wise stochastic
optimization method within the stochastic maximum principle framework.
Numerical experiments validate the efficacy and accuracy of our algorithm,
showcasing its high efficiency in solving feedback optimal control problems
with partial observation.

##  1 Introduction

In this paper, we consider a feedback optimal control problem, in which the
exact state of the controlled process is not available. Stochastic optimal
control is an important research subject that attracts scientists and
engineers in various fields from theoretical scientific research to practical
industrial production. In the stochastic optimal control problem, there is a
control process (also called control policy), which controls a stochastic
dynamical system whose solution called state process, and the goal of the
stochastic optimal control problem is to find an optimal control process to
meet some optimality conditions. For the classic stochastic optimal control
problem with full observation of the state, both theoretical results and
numerical methods are extensively studied. However, in practice the full
observation of the state process is often not available. Instead, we have
detectors/observation facilities to collect partial observational data, which
provide indirect information about the state process. The theoretical
formulation for partially observable stochastic optimal control is derived
analytically [17, 22, 26, 41, 43], and some preliminary numerical method for
the corresponding feedback control problem are developed [3, 5]. To highlight
the influence of data in designing control policies, in this paper we call the
stochastic optimal control with partially observed controlled processes the
data driven feedback control.

Since the controlled state needs to be inferred from observations, the
procedure of finding the optimal control requires data analysis of
observational data. The design of control actions is therefore guided by the
information contained within this data. In this connection, the numerical
approach for solving data-driven feedback control involves both a data
assimilation procedure, which addresses the optimal filtering problem for
estimating the state of the controlled dynamical system with partial noisy
observations, and an optimal control solver. In this research, we present a
novel approach by introducing a kernel learning backward stochastic
differential equation (SDE) filter, which serves as an efficient and accurate
data assimilation technique. For the stochastic optimal control problem, we
propose a sample-wise stochastic optimization method within the stochastic
maximum principle framework.

The optimal filtering problem for data assimilation is typically solved by
recursive implementation of Bayesian inference. The well-known approached
include Kalman type filters and particle filters. The state-of-the-art Kalman
fype filter is the ensemble Kalman filter (EnKF) [20, 42]. The main idea of
EnKF is to use an ensemble of Kalman filter samples to characterize the
probability distribution of the target state as a Gaussian distribution. As a
Kalman type filter, the ensemble Kalman filter stores the information of the
state variable as the mean and the covariance of Kalman filter samples. In
this way, the probability density function (PDF) of the target state, which is
often called the “filtering density”, is approximated as a Gaussian
distribution. However, in nonlinear filtering problems, the filtering density
is usually non-Gaussian. Therefore, the ensemble Kalman filter, which still
relies on the Gaussian assumption, is not the ideal approach to solve the
nonlinear filtering problem.

In addition to the ensemble Kalman filter, various effective nonlinear
filtering methods have been developed to tackle nonlinearity, such as the
particle filter [1, 18, 19, 25, 28, 40, 31], the Zakai filter [12, 45], and
the backward SDE filter [7, 9, 10, 14]. Among these methods, the particle
filter stands out as the most widely utilized approach for addressing
nonlinear filtering problems. This method utilizes a collection of Monte Carlo
samples, referred to as “particles” to construct an empirical distribution
that represents the filtering density of the target state. Upon receiving
observational data, a Bayesian inference procedure is applied to assign
likelihood weights to these particles, and a resampling procedure is
repeatedly executed to generate additional duplicates of particles with higher
weights while discarding particles with lower weights. The particle filter can
handle the nonlinearity inherent in optimal filtering problems by
incorporating nonlinear state dynamics into the filtering density through
particle simulations. However, a significant drawback of the particle filter
is the phenomenon known as “degeneracy issue”, where particles residing in
high probability regions may not adequately capture highly probable features
in filtering densities. This issue becomes even more challenging when dealing
with high-dimensional problems due to the “curse of dimensionality”. Although
advanced resampling methods are proposed to address the degeneracy issue by
relocating particles to high probability regions [1, 18, 36, 40], the finite
particle representation for filtering density inherent in the sequential Monte
Carlo framework restricts the extent of improvements that can be achieved for
the particle filter in solving high-dimensional nonlinear filtering problems.

In this work, we introduce a backward SDE filter approach for solving the
optimal filtering problem. Specifically, a forward backward SDE system is
employed to approximate the prediction of the (controlled) state dynamics.
Then, we apply Bayesian inference to incorporate the observational information
and update the prediction. Numerical experiments have been conducted to show
the advantageous performance of the backward SDE filter over Kalman type
filters and particle filters [14].

For solving the stochastic optimal control problem, there are two well-known
approaches: the dynamic programing principle and the stochastic maximum
principle [15, 21, 34]. Both approaches involve numerical simulations for
large (stochastic) differential systems. In this work, we choose the
stochastic maximum principle approach due to its its advantage in solving
control problems with random coefficients in the state model and the fact that
it could solve problems with finite dimensional terminal state constraints
[44]. The mathematical foundation of the stochastic maximum principle approach
is to derive a system of backward stochastic differential equation (BSDE) as
the adjoint equation of the (forward) controlled state process, and we can use
the solution of the adjoint BSDE to formulate the gradient of the cost
function with respect to the control process. With the gradient process, the
desired optimal control can be solved via gradient descent type optimization.
There are several stochastic maximum principle based numerical methods to
solve the stochastic optimal control problem [23, 30, 41], and the primary
computational cost in those methods lies in obtaining the numerical solution
for the BSDE. Especially, when the dimension of the problem is high, a very
large number of random samples are needed to describe the solution of the
BSDE, which makes the numerical approximation for the gradient process very
expensive [8, 46, 47]. Since we aim to dynamically design control actions
based on estimated parameters, efficiency of the algorithm is very important.

To enhance efficiency, we introduce a novel approach wherein the random
samples used to describe the solution of the adjoint BSDE are treated as “data
samples”. This concept underpins the development of an efficient sample-wise
optimal control solver. The essence of this solver lies in utilizing only one
realization of a random sample (or a small batch of samples) to represent the
solution of the BSDE and derive numerical approximations for the solution
based on this single realization [2, 3, 5]. Additionally, the application of
stochastic approximation in stochastic gradient descent optimization
methodologies [38, 39] justifies the use of sample-wise approximation for
solving the BSDE in the stochastic optimal control solver [4]. In this way,
our approach does not require solving the adjoint BSDE on a large number of
random samples, shifting the computational burden from solving BSDEs to the
optimization procedure for identifying the optimal control. This strategy
ultimately enhances the efficiency of the stochastic optimal control solver.

The rest of this paper is organized as follows. In Section 2, we introduce the
formulation of the data driven feedback control problem and provide our
backward SDE filter and stochastic maximum principle based approach. Section
3, we describe our efficient algorithm for solving the data driven feedback
control in the online manner. Numerical experiments will be presented in
Section 4.

##  2 Problem formulation and the computational framework for the data driven
feedback control problem

###  2.1 Problem setting for the data driven optimal control problem.

In probability space $(\Omega,\mathcal{F},\mathbb{P})$, we consider the
following augmented system on time interval $[0,T]$

$\displaystyle d\begin{pmatrix}S_{t}\\\
M_{t}\end{pmatrix}=\begin{pmatrix}b(t,S_{t},u_{t})\\\
g(S_{t})\end{pmatrix}dt+\begin{pmatrix}\sigma(t,S_{t},u_{t})&0\\\
0&I\end{pmatrix}d\begin{pmatrix}W_{t}\\\
B_{t}\end{pmatrix},\qquad\begin{pmatrix}S_{0}=\xi\\\ M_{0}=0\end{pmatrix},$
(1)

where $S:=\\{S_{t}\\}_{t=0}^{T}$ is the $\mathbb{R}^{d}$-dimensional
controlled state process with dynamics
$b:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{d}$,
$\sigma:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{d\times
q}$ is the diffusion coefficient for a $d$-dimensional Brownian motion $W$
that perturbs the the state $S$, and $u$ is an $m$-dimensional control process
valued in some set $U$ that controls the state process $S$. In the case that
the state $S$ is not directly observable, we have an observation process $M$
that collects partial noisy observations on $S$ with observation function
$g:\mathbb{R}^{d}\rightarrow\mathbb{R}^{p}$, and $B$ is a $p$-dimensional
Brownian motion that is independent from $W$.

Let $\mathbb{F}^{B}=\\{\mathcal{F}_{t}^{B}\\}_{t\geq 0}$ be the filtration of
$B$ augmented by all the $\mathbb{P}$-null sets in $\mathcal{F}$, and
$\mathbb{F}^{W,B}\equiv\\{\mathcal{F}^{W,B}_{t}\\}_{t\geq 0}$ be the
filtration generated by $W$ and $B$ (augmented by $\mathbb{P}$-null sets in
$\mathcal{F}$). Under mild conditions, for any square integrable random
variable $\xi$ independent of $W$ and $B$, and any
$\mathbb{F}^{W,B}$-progressively measurable process $u$ (valued in $U$), Eq.
(1) admits a unique solution $(S,M)$ which is $\mathbb{F}^{W,B}$-adapted.
Next, we let $\mathbb{F}^{M}=\\{\mathcal{F}^{M}_{t}\\}_{t\geq 0}$ be the
filtration generated by $M$ (augmented by all the $\mathbb{P}$-null sets in
$\mathcal{F}$). Clearly, $\mathbb{F}^{M}\subset\mathbb{F}^{W,B}$, and
$\mathbb{F}^{M}\neq\mathbb{F}^{W}$, $\mathbb{F}^{M}\neq\mathbb{F}^{B}$, in
general. The $\mathbb{F}^{M}$ progressively measurable control process,
denoted by $u^{M}$, are control actions driven by the information in
observational data.

We introduce the set of data driven admissible controls as

$\mathcal{U}_{ad}[0,T]=\left\\{u^{M}:[0,T]\times\Omega\rightarrow
U\subset\mathbb{R}^{m}\big{|}u^{M}\text{ is
}\mathbb{F}^{M}-\text{progressively measurable}\right\\},$

and the cost functional that measures the performance of data driven control
$u^{M}$ is defined as

$J(u^{M})={\mathbb{E}}\left[\int_{0}^{T}f(t,S_{t},u_{t}^{M})dt+h(S_{T})\right],$
(2)

where $f$ is the running cost, and $h$ is the terminal cost.

The goal of the data driven feedback control problem is to find the optimal
data driven control $u^{*}\in\mathcal{U}_{ad}[0,T]$ such that

$J(u^{\ast})=\inf_{u^{M}\in\mathcal{U}_{ad}[0,T]}J(u^{M}).$ (3)

###  2.2 Optimization for optimal control

To solve the data driven feedback control problem, we will use the following
optimization procedure derived from the stochastic maximum principle.

When the optimal control $u^{\ast}$ is in the interior of $\mathcal{U}_{ad}$,
the gradient process of the cost functional $J^{\ast}$ with respect to the
control process on time interval $t\in[0,T]$ can be derived using the Gâteaux
derivative of $u^{\ast}$ and the stochastic maximum principle in the following
form(see [6] for details):

$(J^{\ast})_{u}^{{}^{\prime}}(u_{t}^{\ast})=E\left[b_{u}(t,S_{t}^{\ast},u_{t}^{\ast})^{\top}\bar{Y}_{t}+\sigma_{u}(t,S_{t}^{\ast},u_{t}^{\ast})^{\top}\bar{Z}_{t}+f_{u}(t,S_{t}^{\ast},u_{t}^{\ast})^{\top}|\mathcal{F}_{t}^{M}\right]$
(4)

stochastic processes $\bar{Y}$ and $\bar{Z}$ are solutions of the following
forward-backward stochastic differential equations (FBSDEs) system

$\begin{cases}dS_{t}^{\ast}=b(t,S_{t}^{\ast},u_{t}^{\ast})dt+\sigma(t,S_{t}^{\ast},u_{t}^{\ast})dW_{t},\hfill
S_{0}=\xi\\\ dM_{t}^{\ast}=\phi(S_{t})dt+dB_{t},\hfill M_{0}=0\\\
d\bar{Y}_{t}=(-b_{x}(t,S_{t}^{\ast},u_{t}^{\ast})^{\top}\bar{Y}_{t}-\sigma_{x}(t,S_{t}^{\ast},u_{t}^{\ast})^{\top}\bar{Z}_{t}-f_{x}(t,S_{t}^{\ast},u_{t}^{\ast})^{\top})dt\\\
\hskip
50.00008pt+\bar{Z}_{t}dW_{t}+\zeta_{t}dB_{t},\hfill\bar{Y}_{T}=(g_{x}(S_{T}))^{\top}\end{cases}$
(5)

where $\bar{Z}$ is the martingale representation of $\bar{Y}$ with respect to
$W$ and $\zeta$ is the martingale representation of $\bar{Y}$ with respect to
$B$.

To solve the data driven feedback optimal control problem, we also use
gradient descent type optimization and the gradient process
$(J^{\ast})_{u}^{{}^{\prime}}$ is defined in (4). Then, we can use the
following gradient descent iteration to find the optimal control
$u_{t}^{\ast}$ at any time instant $t\in[0,T]$

$u_{t}^{l+1,M}=u_{t}^{l,M}-r(J^{\ast})_{u}^{{}^{\prime}}(u_{t}^{l,M}),\;\;\;\;l=0,1,2,\dots$
(6)

where $r$ is the step size for the gradient. We know that the observational
information $\mathcal{F}_{t}^{M}$ is gradually increased as we collect more
and more data over time. Therefore, at a certain time instant $t$, we target
on finding the optimal control $u_{t}^{\ast}$ with accessible information
$\mathcal{F}_{t}^{M}$. Since the evaluation for
$(J^{\ast})_{u}^{{}^{\prime}}(u_{t}^{l,M})$ requires trajectories
${(\bar{Y}_{\tau},\bar{Z}_{\tau})}_{t\leq\tau\leq T}$ as $\bar{Y}_{t}$ and
$\bar{Z}_{t}$ are solved backward from T to t, we take conditional expectation
$E[\cdot|\mathcal{F}_{t}^{M}]$ to the gradient process
$\\{(J^{\ast})_{u}^{{}^{\prime}}(u_{t}^{l,M})\\}_{t\leq\tau\leq T}$ i.e.

$\displaystyle
E[(J^{\ast})_{u}^{{}^{\prime}}(u_{\tau}^{l,M})|\mathcal{F}_{t}^{M}]$
$\displaystyle=E\Big{[}b_{u}(\tau,S_{\tau},u_{\tau}^{l,M})^{\top}\bar{Y}_{\tau}+\sigma_{u}(\tau,S_{\tau},u_{\tau}^{l,M})^{\top}\bar{Z}_{\tau}$
(7)
$\displaystyle+f_{u}(\tau,S_{\tau},u_{\tau}^{l,M})^{\top}|\mathcal{F}_{t}^{M}\Big{]},\;\;\;\;\;\;\;\;\tau\in[t,T]$

where $S_{\tau}$, $\bar{Y}_{\tau}$ and $\bar{Z}_{\tau}$ are corresponding to
the estimated control $u_{\tau}^{l,M}$. for the gradient descent iteration (6)
on the time interval $[t,T]$, by taking conditional expectation
$E[\cdot|\mathcal{F}_{t}^{M}]$, we obtain

$E[u_{\tau}^{l+1,M}|\mathcal{F}_{t}^{M}]=E[u_{\tau}^{l,M}|\mathcal{F}_{t}^{M}]-rE[(J^{\ast})_{u}^{{}^{\prime}}(u_{\tau}^{l,M})|\mathcal{F}_{t}^{M}],\;\;\;\;l=0,1,2,\dots\;\;\;\tau\in[t,T]$
(8)

When $\tau>t$, the observational information
$\\{\mathcal{F}_{\tau}^{M}\\}_{t\leq\tau\leq T}$ is not available at time t.
We use conditional expectation $E[u_{\tau}^{l+1,M}|\mathcal{F}_{t}^{M}]$ to
replace $u_{\tau}^{l,M}$ since it provides the best approximation. We denote

$u_{\tau}^{l,M}|_{t}:=E[u_{\tau}^{l+1,M}|\mathcal{F}_{t}^{M}]$

and then the gradient descent iteration is

$u_{\tau}^{l+1,M}|_{t}=u_{\tau}^{l,M}|_{t}-rE[(J^{\ast})_{u}^{{}^{\prime}}(u_{\tau}^{l,M}|_{t})|\mathcal{F}_{t}^{M}],\;\;\;\;l=0,1,2,\dots\;\;\;\tau\in[t,T]$
(9)

where
$E[(J^{\ast})_{u}^{{}^{\prime}}(u_{\tau}^{l,M}|_{t})|\mathcal{F}_{t}^{M}]$ can
be obtained by solving the FBSDEs

$\displaystyle dS_{t}$
$\displaystyle=b(\tau,S_{\tau},u_{\tau}^{l,M}|_{t})d\tau+\sigma(\tau,S_{\tau},u_{\tau}^{l,M}|_{t})dW_{\tau},\hskip
130.0002pt\tau\in[t,T]$ (10) $\displaystyle d\bar{Y}_{t}$
$\displaystyle=(-b_{x}(\tau,S_{\tau},u_{\tau}^{l,M}|_{t})^{\top}\bar{Y}_{\tau}-\sigma_{x}(\tau,S_{\tau},u_{\tau}^{l,M}|_{t})^{\top}\bar{Z}_{\tau}-f_{x}(\tau,S_{\tau},u_{\tau}^{l,M}|_{t})^{\top})d\tau$
$\displaystyle\hskip
50.00008pt+\bar{Z}_{\tau}dW_{\tau}+\zeta_{\tau}dB_{\tau},\hskip
150.00023pt\bar{Y}_{T}=(g_{x}(S_{T}))^{\top}$

and evaluated effectively using the numerical algorithm introduced later.

When the controlled dynamics and the observation function $\phi$ are
nonlinear, we will use the optimal filtering framework to obtain the
conditional expectation $E[\Psi(\tau)|\mathcal{F}_{t}^{M}]$. For simplicity,
we define

$\Psi(\tau,S_{\tau},u_{\tau}^{l,M}|_{t}):=b_{u}^{\top}(\tau,S_{\tau},u_{\tau}^{l,M}|_{t})\bar{Y}_{\tau}+\sigma_{u}^{\top}(\tau,S_{\tau},u_{\tau}^{l,M}|_{t})\bar{Z}_{\tau}+f_{u}^{\top}(\tau,S_{\tau},u_{\tau}^{l,M}|_{t})$
(11)

for $\tau\in[t,T]$. With the conditional pdf $p(S_{t}|\mathcal{F}_{t}^{M})$
that we obtain through optimal filtering methods and the fact that
$\Psi(\tau,S_{\tau},u_{\tau}^{l,M}|_{t})$ is a stochastic process depending on
the state of random variable $S_{t}$, the conditional gradient process
$E[(J^{\ast})_{u}^{{}^{\prime}}(u_{t}^{l,M}|_{t})|\mathcal{F}_{t}^{M}]$ in (7)
can be obtained by the following integral

$E[(J^{\ast})_{u}^{{}^{\prime}}(u_{\tau}^{l,M}|_{t})|\mathcal{F}_{t}^{M}]=\int_{\mathbb{R}^{d}}E[\Psi(\tau,S_{\tau},u_{\tau}^{l,M}|_{t})|S_{t}]\cdot
p(S_{t}|\mathcal{F}_{t}^{M})dx,\hskip 20.00003pt\tau\in[t,T]$ (12)

where $(\bar{Y}_{\tau},\bar{Z}_{\tau})$ in $\Psi$ can be approximated by the
numerical scheme of the FBSDE system. $p(S_{t}|\mathcal{F}_{t}^{M})$ can be
approximated by the distribution of the backward SDE filter method introduced
in the next section.

###  2.3 Estimation for the state - nonlinear filtering

In this work, we focus on the Bayes filter approach, which is carried out by
recursive Bayesian estimations. In the Bayes filter, we estimate the target
state on a sequence of discrete time instants
$0=t_{0}<t_{1}<\cdots<t_{N_{T}}=T$ over the time interval $[0,T]$, where
$N_{T}\in\mathbb{N}$ is the total number of time steps. The general framework
of recursive Bayesian estimations is composed of two steps: a prediction step
and an update step.

In the prediction step, assuming that the filtering density
$p(S_{t_{n}}|\mathcal{F}_{t_{n}}^{M})$ is available at the time instant
$t_{n}$, we use the Chapman-Kolmogorov formula to propagate the dynamical
model in the state equation as follows:

$p(S_{t_{n+1}}|\mathcal{F}_{t_{n}}^{M})=\int
p(S_{t_{n+1}}|S_{t_{n}})p(S_{t_{n}}|\mathcal{F}_{t_{n}}^{M})dS_{t_{n}}$ (13)

where $p(S_{t_{n+1}}|S_{t_{n}})$ is the transition probability of the state
equation in (1), and the predicted filtering density
$p(S_{t_{n+1}}|\mathcal{F}_{t_{n}}^{M})$, which is the prior distribution in
the Bayesian inference, describes the state $S$ at the time instant $t_{n+1}$
before receiving the new observational data $M_{t_{n+1}}$.

In the update step, we use the following Bayesian inference formula to
incorporate the observational data into the state estimation:

$p(S_{t_{n+1}}|\mathcal{F}_{t_{n+1}}^{M})=\frac{p(M_{t_{n+1}}|S_{t_{n+1}})p(S_{t_{n+1}}|\mathcal{F}_{t_{n}}^{M})}{p(M_{t_{n+1}}|\mathcal{F}_{t_{n}}^{M})}$
(14)

where, $p(M_{t_{n+1}}|S_{t_{n+1}})$ is the likelihood function, and the
denominator $p(M_{t_{n+1}}|\mathcal{F}_{t_{n}}^{M})$ is a normalization
factor.

Then, we introduce the backward SDE filter as the theoretical preparation for
our kernel learning method.

####  2.3.1 The backward SDE filter

The backward SDE filter adopts the recursive Bayesian estimations framework.
The central idea is to use a system of (time-inverse) forward-backward
stochastic differential equations to propagate the filtering density, and we
also use Bayesian inference to incorporate the observational information into
the predicted filtering density.

To proceed, we first introduce the forward-backward stochastic differential
equations (FBSDEs) corresponding to the nonlinear filtering problem of the
observation function in (1), and we consider the following FBSDEs system:

$\displaystyle S_{t}$
$\displaystyle=S_{0}+\int_{0}^{t}b(S_{s})\,ds+\int_{0}^{t}\sigma_{s}\,dW_{s},$
(15) $\displaystyle\tilde{Y}_{0}$
$\displaystyle=\tilde{Y}_{t}-\int_{0}^{t}\tilde{Z}_{s}\,dW_{s},\hskip
56.9055pt\tilde{Y}_{t}=\Phi(S_{t}),$

where the first equation coincides with the state equation in the nonlinear
filtering problem (a standard forward SDE), and the second equation is a
backward SDE. The solution of the above FBSDEs system is the pair
$(\tilde{Y},\tilde{Z})$, which is adapted to the Brownian motion $W$, i.e.,
$\tilde{Y}_{t},\tilde{Z}_{t}\in\mathcal{F}_{t}^{W}$, and the solution
$\tilde{Z}$ is the martingale representation of $\tilde{Y}$ with respect to
$W$ [35]. $\Phi$ is a test function representing the quantity of interest.
Note that the side condition of the backward SDE, i.e.,
$\tilde{Y}_{t}=\Phi(S_{t})$, is given at the time instant $t$, and the
solution pair $(\tilde{Y},\tilde{Z})$ propagates backward from $t$ to $0$.

For a fixed initial state $S_{0}=x\in\mathbb{R}^{d}$, we take the conditional
expectation $E[\cdot|S_{0}=x]$ on both sides of the backward SDE in (15) and
obtain $\tilde{Y}_{0}(x)=E[\Phi(S_{t})|S_{0}=x]$, which is a simplified
version of the Feynman-Kac formula. Here, we note that the value of
$\tilde{Y}_{0}$ is determined by the value of the state $S_{0}$. In addition,
the solution $\tilde{Y}$ of (15) is equivalent to the solution of the
Kolmogorov backward equation. In other words, for the following backward
parabolic type partial differential equation (PDE),

$-\frac{du_{s}}{ds}=\sum_{i=1}^{d}b_{i}\frac{\partial u_{s}}{\partial
x_{i}}+\frac{1}{2}\sum_{i,j=1}^{d}(\sigma_{s}\sigma_{s}^{T})_{i,j}\frac{\partial^{2}u_{s}}{\partial
x_{i}\partial x_{j}},\quad\quad u_{t}(x)=\Phi(x)$ (16)

where $b_{i}$ is the $i$-th component of the vector function $b$ in the state
equation in (1), and we have $\tilde{Y}_{0}(S_{0}=x)=u_{0}(x)$ [33].

In the nonlinear filtering problem, we need to propagate the filtering density
forward. The PDE that propagates the PDF of the state $S$ driven by (1) is the
following Fokker-Planck equation:

$\frac{dp_{t}}{dt}=-\sum_{i=1}^{d}\frac{\partial b_{i}p_{t}}{\partial
x_{i}}+\frac{1}{2}\sum_{i,j=1}^{d}(\sigma_{t}\sigma_{t}^{T})_{i,j}\frac{\partial^{2}p_{t}}{\partial
x_{i}\partial x_{j}}$ (17)

where the initial condition $p_{0}$ is the distribution of the state $S_{0}$.
Equation (17) is the adjoint equation of the Kolmogorov backward equation, and
thus, the Fokker-Planck equation is also called the Kolmogorov forward
equation.

Following the analysis in [33] that establishes the equivalence between the
FBSDEs (15) and the Kolmogorov backward equation (16), one can derive that the
solution $p_{t}$ of the above Fokker-Planck equation is equivalent to the
solution $Y_{t}$ of the following FBSDEs:

$\displaystyle X_{0}$
$\displaystyle=X_{t}-\int_{0}^{t}b(X_{s})\,ds+\int_{0}^{t}\sigma_{s}\,d\overleftarrow{W_{s}},$
(18a) $\displaystyle Y_{t}$
$\displaystyle=Y_{0}-\int_{0}^{t}\sum_{i=1}^{d}\frac{\partial b_{i}}{\partial
x_{i}}(X_{s})Y_{s}\,ds-\int_{0}^{t}Y_{s}\,d\overleftarrow{W_{s}},\quad
Y_{0}=p_{0}.$ (18b)

where the integral $\int_{0}^{t}\cdot d\overleftarrow{W_{s}}$ is a backward
Itô integral, integrated backward [32]. The equation (18a) is a ”backward SDE”
as it propagates backward from $t$ to $0$, and the equation (18b) is a
”forward SDE” since its side condition is given at the time instant $0$. On
the other hand, the propagation direction of (18a) is the same as the
integration direction of the backward Itô integral. At the same time, the
equation (18b) propagates forwards given the side condition $Y_{0}$ = $p_{0}$
with a backward Itô integral. In this way,the equations in (18) compose a
time-inverse FBSDEs system.

Due to the equivalence $Y_{t}=p_{t}$, the solution $Y_{t}$ also propagates the
PDF of the state $S_{t}$ forward. In this connection, we use the time-inverse
FBSDEs system (18) to predict the filtering density in the nonlinear filtering
problem. Specifically, assuming we have the filtering density
$p(S_{t_{n}}|\mathcal{F}_{t}^{M})$ at the time step $t_{n}$, we solve the
following time-inverse FBSDEs system:

$\displaystyle X_{t_{n}}$
$\displaystyle=X_{t_{n+1}}-\int_{t_{n}}^{t_{n+1}}b\left(X_{s}\right)ds+\int_{t_{n}}^{t_{n+1}}\sigma_{s}d\overleftarrow{W}_{s}$
(19a) $\displaystyle Y_{t_{n+1}}^{M_{t_{n}}}$
$\displaystyle=Y_{t_{n}}^{M_{t_{n}}}-\int_{t_{n}}^{t_{n+1}}\sum_{i=1}^{d}\frac{\partial
b_{i}}{\partial
x_{i}}\left(X_{s}\right)Y_{s}^{M_{t_{n}}}ds-\int_{t_{n}}^{t_{n+1}}Z_{s}^{M_{t_{n}}}d\overleftarrow{W}_{s}$
(19b) $\displaystyle X_{t_{n+1}}$ $\displaystyle=S_{t_{n+1}},\quad\quad\quad
Y_{t_{n}}^{M_{t_{n}}}=p\left(S_{t_{n}}\mid\mathcal{F}_{t_{n}}^{M}\right)$
(19c)

we obtain the solution $Y_{t_{n+1}}^{M_{t_{n}}}$, which is the predicted
filtering density, i.e.,
$Y_{t_{n+1}}^{M_{t_{n}}}=p\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n}}^{M}\right)$.
In other words, the time-inverse FBSDEs system (19) provides a mechanism to
carry out the Chapman-Kolmogorov formula (13) in the prediction step of the
Bayes filter approach. Here, we use the superscript $M_{t_{n}}$ in
$Y_{t_{n+1}}^{M_{t_{n}}}$ to emphasize that the solution $Y$ depends on the
observational information $\mathcal{F}_{t_{n}}^{M}$.

Then, we substitute the prior distribution in the Bayesian inference (14) by
the solution $Y_{t_{n+1}}^{M_{t_{n}}}$ obtained in (19) to get the posterior
distribution $p\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n+1}}^{M}\right)$, i.e.

$p\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n+1}}^{M}\right)=\frac{p\left(M_{t_{n+1}}\mid
S_{t_{n+1}}\right)Y_{t_{n+1}}^{M_{t_{n}}}}{p\left(M_{t_{n+1}}\mid\mathcal{F}_{t_{n}}^{M}\right)}$
(20)

which is the updated filtering density at time step $t_{n+1}$ that will be
used for the next recursive stage.

From the above discussion, we can see that the central idea of the backward
SDE filter is to use the time-inverse FBSDEs system to predict the filtering
density, and then use Bayesian inference to update the predicted filtering
density. As a method that carries out recursive Bayesian estimations, the
backward SDE filter is also composed of a prediction step and an update step.
In most practical situations, FBSDEs are not explicitly solvable. Therefore,
numerical solutions for FBSDEs are needed. In the following section, we
introduce a numerical algorithm to implement the above backward SDE filter
framework and develop our efficient kernel learning method for the backward
SDE filter.

##  3 Numerical approach for data driven feedback control by BSDE-SGD

For the numerical framework, we need the temporal partition $\Pi_{N_{T}}$

$\Pi_{N_{T}}=\\{t_{n},0=t_{0}<t_{1}<\cdots<t_{N_{T}}=T\\}$

and we use the control sequence $\\{u_{t_{n}}^{\ast}\\}_{n=1}^{N_{T}}$ to
represent the control process $u^{\ast}$.

###  3.1 Numerical schemes for the FBSDE in optimal control

For the FBSDE system (5), we get the following schemes:

$\displaystyle S_{i+1}=S_{i}+b(t_{i},S_{i},u_{t_{i}}^{l,M}|_{t_{n}})\triangle
t_{i}+\sigma(t_{i},S_{i},u_{t_{i}}^{l,M}|_{t_{n}})\triangle W_{t_{i}}$ (21)
$\displaystyle\bar{Y}_{i}=E_{i}[\bar{Y}_{i+1}]+E_{i}[b_{x}(t_{i+1},S_{i+1},u_{t_{i+1}}^{l,M}|_{t_{n}})^{\top}\bar{Y}_{i+1}$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+\sigma_{x}(t_{i+1},S_{i+1},u_{t_{i+1}}^{l,M}|_{t_{n}})^{\top}\bar{Z}_{i+1}+f_{x}(t_{i+1},S_{i+1},u_{t_{i+1}}^{l,M}|_{t_{n}})^{\top}]\triangle
t_{i}$ $\displaystyle\bar{Z}_{i}=\frac{1}{\triangle
t_{i}}E_{i}[\bar{Y}_{i+1}\triangle W_{t_{i}}]$

where $S_{i+1},\bar{Y}_{i},\;and\;\bar{Z}_{i}$ are numerical approximations
for $S_{t_{i+1}},\bar{Y}_{t_{i}},\;and\;\bar{Z}_{t_{i}}$, respectively.

Then, the standard Monte Carlo method can approximate expectations with K
random samples:

$\displaystyle
S_{i+1}^{k}=S_{i}+b(t_{i},S_{i},u_{t_{i}}^{l,M}|_{t_{n}})\triangle
t_{i}+\sigma(t_{i},S_{i},u_{t_{i}}^{l,M}|_{t_{n}})\sqrt{\triangle
t_{i}}\omega_{i}^{k},k=1,2,\dots,K,$ (22)
$\displaystyle\bar{Y}_{i}=\sum_{k=1}^{K}\frac{\bar{Y}_{i+1}^{k}}{K}+\frac{\triangle
t_{i}}{K}\sum_{k=1}^{K}[b_{x}(t_{i+1},S_{i+1}^{k},u_{t_{i+1}}^{l,M}|_{t_{n}})^{\top}\bar{Y}_{i+1}^{k}$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+\sigma_{x}(t_{i+1},S_{i+1}^{k},u_{t_{i+1}}^{l,M}|_{t_{n}})^{\top}\bar{Z}_{i+1}^{k}+f_{x}(t_{i+1},S_{i+1}^{k},u_{t_{i+1}}^{l,M}|_{t_{n}})^{\top}]$
$\displaystyle\bar{Z}_{i}=\frac{1}{\triangle
t_{i}}\sum_{k=1}^{K}\frac{\bar{Y}_{i+1}^{k}\sqrt{\triangle
t_{i}}\omega_{i}^{k}}{K}$

where $\\{\omega_{i}^{k}\\}_{k=1}^{K}$ is a set of random samples following
the standard Gaussian distribution that we use to describe the randomness of
$\triangle W_{t_{i}}$.

The above schemes solve the FBSDE system (5) as a recursive algorithm, and the
convergence of the schemes is well studied — cf. ([13],[11], and [48]).

###  3.2 Numerical algorithms for kernel learning backward SDE filter

We first provide numerical schemes to solve the time-inverse FBSDEs system
(19) and then give Bayesian inference based on the numerical solution of the
FBSDEs system in Subsection 3.2.1. To approximate the entire filtering density
that carries the information contained in the state dynamical model and the
observational data, in Subsection 3.2.2 we introduce a kernel learning method
to ”learn” the filtering density from discrete density values. In Subsection
3.3 we summarize our kernel learning backward SDE filter and apply it to the
stochastic optimization for the control process.

####  3.2.1 Numerical schemes for time-inverse FBSDEs

Since the equation (19a) is essentially an SDE with inverse propagation
direction, we apply the Euler-Maruyama scheme [29] and get

$X_{t_{n}}=X_{t_{n+1}}-b\left(X_{t_{n+1}}\right)\Delta
t_{n}+\sigma_{t_{n+1}}\Delta W_{t_{n}}+R_{n}^{X}$

where $\Delta t_{n}:=t_{n+1}-t_{n},\Delta W_{t_{n}}:=W_{t_{n+1}}-W_{t_{n}}$,
and $R_{n}^{X}:=b\left(X_{t_{n+1}}\right)\Delta
t_{n}-\int_{t_{n}}^{t_{n+1}}b\left(X_{s}\right)ds$ is the approximation error.
By dropping the error term $R_{n}^{X}$, we obtain the following numerical
scheme for (19a):

$X_{n}=X_{n+1}-b\left(X_{n+1}\right)\Delta t_{n}+\sigma_{t_{n+1}}\Delta
W_{t_{n}}$ (23)

where $X_{n}$ is the numerical approximation for $X_{t_{n}}$, and $X_{n+1}$ is
a representation for the state variable $S_{t_{n+1}}$.

To solve the backward SDE (19b), we take conditional expectation
$\mathbb{E}_{n+1}^{X}[\cdot]:=\mathbb{E}\left[\cdot\mid
X_{t_{n+1}},\mathcal{F}_{t_{n}}^{M}\right]$ on both sides of the equation and
obtain

$\mathbb{E}_{n+1}^{X}\left[Y_{t_{n+1}}^{M_{t_{n}}}\right]=\mathbb{E}_{n+1}^{X}\left[Y_{t_{n}}^{M_{t_{n}}}\right]-\int_{t_{n}}^{t_{n+1}}\mathbb{E}_{n+1}^{X}\left[\sum_{i=1}^{d}\frac{\partial
b_{i}}{\partial x_{i}}\left(X_{s}\right)Y_{s}^{M_{t_{n}}}\right]ds$ (24)

where the backward Itô integral
$\int_{t_{n}}^{t_{n+1}}Z_{s}d\overleftarrow{W}_{s}$ is eliminated due to the
martingale property of Itô integrals. For the left hand side of the above
equation, we have
$Y_{t_{n+1}}^{M_{t_{n}}}=\mathbb{E}_{n+1}^{X}\left[Y_{t_{n+1}}^{M_{t_{n}}}\right]$
since $Y_{t_{n+1}}^{M_{t_{n}}}$ is adapted to $X_{t_{n+1}}$ and is
$\mathcal{F}_{t_{n}}^{M}$ measurable. In this work, we approximate the
deterministic integral on the right hand side of (24) by using the right-point
formula and get

$Y_{t_{n+1}}^{M_{t_{n}}}=\mathbb{E}_{n+1}^{X}\left[Y_{t_{n}}^{M_{t_{n}}}\right]-\Delta
t_{n}\cdot\sum_{i=1}^{d}\frac{\partial b_{i}}{\partial
x_{i}}\left(X_{t_{n+1}}\right)Y_{t_{n+1}}^{M_{t_{n}}}+R_{n}^{Y}$ (25)

where $R_{n}^{Y}:=\Delta t_{n}\cdot\sum_{i=1}^{d}\frac{\partial
b_{i}}{\partial
x_{i}}\left(X_{t_{n+1}}\right)Y_{t_{n+1}}^{M_{t_{n}}}-\int_{t_{n}}^{t_{n+1}}\mathbb{E}_{n+1}^{X}\left[\sum_{i=1}^{d}\frac{\partial
b_{i}}{\partial x_{i}}\left(X_{s}\right)Y_{s}^{M_{t_{n}}}\right]ds$ is the
approximation error for the integral, and we have used the fact
$\sum_{i=1}^{d}\frac{\partial b_{i}}{\partial
x_{i}}\left(X_{t_{n+1}}\right)Y_{t_{n+1}}^{M_{t_{n}}}=$
$\mathbb{E}_{n+1}^{X}\left[\sum_{i=1}^{d}\frac{\partial b_{i}}{\partial
x_{i}}\left(X_{t_{n+1}}\right)Y_{t_{n+1}}^{M_{t_{n}}}\right]$.

Then, we drop the approximation error term $R_{n}^{Y}$ in (25) and obtain the
following approxi- mation scheme for $Y_{t_{n+1}}^{M_{t_{n}}}$ :

$Y_{n+1}^{M_{t_{n}}}=\mathbb{E}_{n+1}^{X}\left[Y_{n}^{M_{t_{n}}}\right]-\Delta
t_{n}\cdot\sum_{i=1}^{d}\frac{\partial b_{i}}{\partial
x_{i}}\left(X_{t_{n+1}}\right)Y_{n+1}^{M_{t_{n}}}$ (26)

where $Y_{n+1}^{M_{t_{n}}}$ is the approximated solution and
$Y_{n}^{M_{t_{n}}}$ is an approximation of the filtering density
$p\left(S_{t_{n}}\mid\mathcal{F}_{t_{n}}^{M}\right)$ that we obtained in the
previous recursive step. We can see that the above approximation scheme is an
implicit scheme. In order to calculate $Y_{n+1}^{M_{t_{n}}}$, we introduce the
following fixed-point iteration procedure:

$Y_{n+1}^{M_{t_{n}},l+1}=\mathbb{E}_{n+1}^{X}\left[Y_{n}^{M_{t_{n}}}\right]-\Delta
t_{n}\cdot\sum_{i=1}^{d}\frac{\partial b_{i}}{\partial
x_{i}}\left(X_{t_{n+1}}\right)Y_{n+1}^{M_{t_{n}},l},\quad l=0,1,2,\cdots,L-1$
(27)

and we let the approximated solution be
$Y_{n+1}^{M_{t_{n}}}=Y_{n+1}^{M_{t_{n}},L}$, where $L$ is a number that
satisfies certain stopping criteria for the iterations, and we let the initial
guess for the solution $Y_{n+1}^{M_{t_{n}}}$ be
$Y_{n+1}^{M_{t_{n}},0}=Y_{n}^{M_{t_{n}}}$.

From the equivalence
$Y_{t_{n+1}}^{M_{t_{n}}}=p\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n}}^{M}\right)$,
the approximated solution $Y_{n+1}^{M_{t_{n}}}$ gives an approximation for the
predicted filtering density. Hence the iterative scheme (27) accomplishes the
prediction step in Bayesian estimation. To incorporate the observational
information, we carry out the update step through Bayesian inference as
follows

$\tilde{p}\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n+1}}^{M}\right)=\frac{p\left(M_{t_{n+1}}\mid
S_{t_{n+1}}\right)Y_{n+1}^{M_{t_{n}}}}{C}$ (28)

where $C$ is a normalization factor, and the prior distribution is replaced by
the approximation of the predicted filtering density $Y_{n+1}^{M_{t_{n}}}$. As
a result, we obtain the approximated filtering density
$\tilde{p}\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n+1}}^{M}\right)$ as desired in
the backward SDE filter.

The numerical schemes (27)-(28) compose a general computational framework for
the backward SDE filter, which provides a recursive prediction-update
mechanism that formulates the temporal propagation of the filtering density.
On the other hand, the filtering density is a function that connects state
positions to the probability density values at those positions. Therefore,
spatial dimension approximation for the filtering density function with
respect to the state variable is needed. In the following subsection, we
introduce a kernel learning method to generate a continuous global
approximation for the filtering density over the state space.

####  3.2.2 Efficient kernel learning in the backward SDE filter

In the scheme (27), we need to evaluate the conditional expectation
$\mathbb{E}_{n+1}^{X}\left[Y_{n}^{M_{t_{n}}}\right]$ to calculate the
predicted filtering density $Y_{n+1}^{M_{t_{n}}}$. So, we first discuss the
approximation method for $\mathbb{E}_{n+1}^{x}\left[Y_{n}^{M_{t_{n}}}\right]$
given that $X_{t_{n+1}}=x$.

Approximating the conditional expectation

When the target state is a high dimensional variable, Monte Carlo simulation
is usually applied to evaluate high dimensional expectations. Specifically,
for a given point $x\in\mathbb{R}^{d}$ in the state space, we let $X_{n+1}=x$
in the scheme (23). Then, the conditional expectation in (27) is approximated
by

$\mathbb{E}_{n+1}^{x}\left[Y_{n}^{M_{t_{n}}}\right]\approx\hat{\mathbb{E}}_{n+1}^{x}\left[Y_{n}^{M_{t_{n}}}\right]:=\frac{\sum_{m=1}^{M}Y_{n}^{M_{t_{n}}}\left(\tilde{X}_{n}^{x,m}\right)}{M}$
(29)

where $M$ is the total number of Monte Carlo samples that we use to
approximate the conditional expectation. The random variable
$\tilde{X}_{n}^{x,m}$ in the above approximation is a Monte Carlo sample
simulated by using the scheme (23) as follows:

$\tilde{X}_{n}^{x,m}=x-b(x)\Delta t_{n}+\sigma_{t_{n+1}}\sqrt{\Delta
t_{n}}\omega_{m}$ (30)

where $\omega_{m}$ is a sample drawn from the $d$-dimensional standard
Gaussian distribution, and $\sqrt{\Delta t_{n}}\omega_{m}$ is the $m$-th
realization of $\Delta W_{t_{n}}$. However, when the dimension of the state
variable is high, the number $M$ of Monte Carlo samples needs to be very
large, hence evaluating the conditional expectation
$\mathbb{E}_{n+1}^{x}\left[Y_{n}^{M_{t_{n}}}\right]$ by using the Monte Carlo
simulation (29) is a computationally expensive task.

Inspired by the stochastic approximation method and its application in
stochastic optimization [37], [16], we treat the large number of Monte Carlo
samples that we use to approximate the conditional expectation
$\mathbb{E}_{n+1}^{x}\left[Y_{n}^{M_{t_{n}}}\right]$ as a ”large data set”.
Then, we adopt the methodology of stochastic approximation and use a single-
sample (or a small batch of samples) to represent conditional expectations.
Specifically, in each fixed-point iteration step (27), instead of using the
fully-calculated Monte Carlo simulation (29) to compute the conditional
expectation, one may use one realization of the simulated sample
$\tilde{X}_{n}^{x,l}$ to represent the entire set of Monte Carlo samples
$\left\\{\tilde{X}_{n}^{x,m}\right\\}_{m=1}^{M}$ at each iteration step, where
$\tilde{X}_{n}^{x,l}$ is also simulated through (30) indexed by the iteration
step $l$. In this way, the fixed-point iteration scheme for the approximated
solution $Y_{n+1}^{M_{t_{n}}}$ at the spatial point $X_{n+1}=x$ can be carried
out as follows:

$Y_{n+1}^{M_{t_{n}},l+1}(x)=\tilde{\mathbb{E}}_{n+1}^{x,l}\left[Y_{n}^{M_{t_{n}}}\right]-\Delta
t_{n}\cdot\sum_{i=1}^{d}\frac{\partial b_{i}}{\partial
x_{i}}(x)Y_{n+1}^{M_{t_{n}},l},\quad l=0,1,2,\cdots,L-1$ (31)

where the conditional expectation
$\mathbb{E}_{n+1}^{x}\left[Y_{n}^{M_{t_{n}}}\right]$ is represented by a
single-sample of $Y_{n}^{M_{t_{n}}}$ corresponding to $\tilde{X}_{n}^{x,l}$,
and we have

$\tilde{\mathbb{E}}_{n+1}^{x,l}\left[Y_{n}^{M_{t_{n}}}\right]=Y_{n}^{M_{t_{n}}}\left(\tilde{X}_{n}^{x,l}\right).$
(32)

As a result, in each fixed-point iteration step, we only need to generate one
sample of $X_{n}$ and evaluate the function value of the previous filtering
density $Y_{n}^{M_{t_{n}}}$ at the spatial point $\tilde{X}_{n}^{x,l}$. In
this way, we transfer the cost of simulating a large number of Monte Carlo
samples to carrying out fixed-point iterations. Although the single-sample
representation does not provide an accurate approximation for the conditional
expectation, every simulated sample $\tilde{X}_{n}^{x,l}$ is effectively used
to improve the estimate of the desired predicted filtering density
$Y_{n+1}^{M_{t_{n}}}(x)$, which makes the overall fixed-point iteration
procedure more efficient.

To use the simulated samples more effectively and to make the approximation
for the conditional expectation more accurate, we modify the single-sample
representation (32) and use the batch of samples
$\left\\{Y_{n}^{M_{t_{n}}}\left(\tilde{X}_{n}^{x,l}\right)\right\\}_{l=0}^{L-1}$
to approximate the expectation. Precisely, we let the approximated expectation
in (31) at each iteration step be

$\tilde{\mathbb{E}}_{n+1}^{x,l}\left[Y_{n}^{M_{t_{n}}}\right]=\frac{\sum_{i=1}^{l}Y_{n}^{M_{t_{n}}}\left(\tilde{X}_{n}^{x,i}\right)}{l}$
(33)

In this way, all the samples previously generated are used to evaluate the
expectation. Note that the expectation
$\mathbb{E}_{n+1}^{x}\left[Y_{n}^{M_{t_{n}}}\right]$ in the fixed-point
iteration is independent of the estimation for the solution
$Y_{n+1}^{M_{t_{n}}}$. Hence using more samples to approximate the expectation
at each iteration step only makes the iterative scheme more accurate.

Approximating the filtering density on random spatial sample points

By using the iterative scheme (31), we can calculate the predicted filtering
density $Y_{n+1}^{M_{t_{n}}}$ on the state point $x$. Then, through the
Bayesian inference scheme (28), we can obtain an approximation for the updated
filtering density
$\tilde{p}\left(S_{t_{n+1}}=x\mid\mathcal{F}_{t_{n+1}}^{M}\right)$ at the time
step $t_{n+1}$.

To provide a complete description for the filtering density, we need to
approximate the conditional PDF of the target state as a mapping from the
state variable to PDF values. Standard function approximation methods use
tensor-product grid points, on which we approximate function values, and then
use polynomial interpolation to construct an interpolatory approximation for
the entire function. If the dimension of the problem is moderately high,
sparse-grid methods are often adopted as efficient alternatives to tensor-
product grid interpolations. However, even advanced adaptive sparse-grid
methods suffer from the ”curse of dimensionality” problem. When the dimension
of the problem is higher, i.e., $d\geq 10$, the cost of implementing sparse-
grid approximation becomes extremely high. In many practical nonlinear
filtering problems, the state dimensions are very high. Hence, applying
traditional grid-based function approximation methods for filtering density is
infeasible in solving high dimensional real-world problems.

An advantage of the backward SDE filter is that it allows us to approximate
the filtering density on any point $x$ in the state space. In this way, we
don’t have to solve the problem on pre-determined meshes. Instead, in this
work we use randomly generated state samples as our spatial points, and we
generate spatial sample points so that they adaptively follow the conditional
distribution of the target state. Specifically, assuming that we have a set of
spatial points $\left\\{x_{i}^{n}\right\\}_{i=1}^{N}$ that follow the previous
approximated filtering density
$\tilde{p}\left(S_{t_{n}}\mid\mathcal{F}_{t_{n}}^{M}\right)$, where $N$ is the
total number of spatial points, we propagate those spatial samples through the
state dynamics in (1) by using the following Euler-Maruyama scheme:

$\tilde{x}_{i}^{n+1}=x_{i}^{n}+b\left(x_{i}^{n}\right)\Delta
t_{n}+\sigma_{t_{n}}\sqrt{\Delta t_{n}}\omega^{i},\quad i=1,2,\cdots,N$ (34)

where $\left\\{\omega^{i}\right\\}_{i=1}^{N}$ is a sequence of i.i.d. standard
$d$-dimensional Gaussian random variables. As a result, the sample set
$\left\\{\tilde{x}_{i}^{n+1}\right\\}_{i=1}^{N}$ forms an empirical
distribution for the prior distribution
$p\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n}}^{M}\right)$. Then, we solve the
time-inverse FBSDEs (19) on those spatial sample points through the scheme
(31) to get
$\left\\{Y_{n+1}^{M_{t_{n}}}\left(\tilde{x}_{i}^{n+1}\right)\right\\}_{i=1}^{N}$.
When the new observational data $M_{t_{n+1}}$ is available, we use Bayesian
inference to update the predicted filtering density and get
$\left\\{\tilde{p}\left(S_{t_{n+1}}=\right.\right.$
$\left.\left.\tilde{x}_{i}^{n+1}\mid\mathcal{F}_{t_{n+1}}^{M}\right)\right\\}_{i=1}^{N}$
as our approximations for the updated filtering density. In this way, the
approximations on the scattered sample points
$\left\\{\tilde{x}_{i}^{n+1}\right\\}_{i=1}^{N}$ provide a partial description
for the desired filtering density.

On the other hand, in our iterative scheme (31) we need the function value of
$Y_{n}^{M_{t_{n}}}$ on the spatial point
$\tilde{X}_{n}^{\tilde{x}_{i}^{n+1},l}$, which is calculated from the scheme
(30) by choosing $x=\tilde{x}_{i}^{n+1}$. Apparently,
$\tilde{X}_{n}^{\tilde{x}_{i}^{n+1},l}$ is unlikely to be one of the existing
sample points, on which we have approximated the function values for
$Y_{n}^{M_{t_{n}}}$. Therefore, we need to derive an approximation for the
filtering density over the entire state space. Since the filtering density is
approximated on random spatial points, mesh-free methods are needed. Although
traditional mesh-free interpolation methods, such as the moving least square
method and the radial basis function interpolation method, could compute
interpolatory approximation based on density function values at nearby spatial
points, calculating the filtering density at each point separately is
computationally expensive. Especially, when the dimension of the problem is
high, we have to approximate the filtering density on a very large number of
spatial points, which makes local approximation methods very difficult to
implement.

Then, we introduce a kernel learning method to ”learn” a global approximation
for the entire filtering density.

Kernel learning for the filtering density

The kernel machine utilizes the combination of a set of pre-chosen kernels to
represent a model [27]. For a target model in the form of a function $F$, the
kernel learning method approximates the function as $F(x)\approx
g\left(\sum_{k=1}^{K}\alpha_{k}\phi_{k}(x)+\beta\right),x\in\mathcal{R}^{d}$,
where $\left\\{\phi_{k}\right\\}_{k=1}^{K}$ is a set of $K$ kernels, $g$ is an
optional nonlinearity, $\left\\{\alpha_{k}\right\\}_{k=1}^{K}$ are weights of
kernels, and $\beta$ is a bias parameter. Then we can apply the kernel
learning method to approximate the filtering density under the backward SDE
filter framework. Since our target function in the nonlinear filtering problem
is a probability distribution, which is nonnegative and often bell-shaped, we
drop the nonlinear function $g$ and the bias $\beta$ in the kernel learning
model, and we choose Gaussian type functions as our kernels. Specifically, at
the time step $t_{n+1}$, we use the following kernel learning scheme to
formulate a global approximation for the filtering density
$p\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n+1}}^{M}\right)$ :

$p_{n+1}(x):=\sum_{k=1}^{K}\alpha_{k}^{n+1}\phi_{k}^{n+1}(x),\quad
x\in\mathcal{R}^{d}$ (35)

where $p_{n+1}$ is the kernel learned filtering density at the time step
$t_{n+1}$, the Gaussian type kernel is chosen as
$\phi_{k}^{n+1}(x)=\exp\left(-\left(\hat{x}_{k}^{n+1}-x\right)^{2}/\left(\lambda_{k}^{n+1}\right)^{2}\right)$
with center $\hat{x}_{k}^{n+1}$ and covariance $\lambda_{k}^{n+1}$, and the
coefficient $\alpha_{k}^{n+1}>0$ is the weight of the $k$-th Gaussian kernel
$\phi_{k}^{n+1}$. We can see from the scheme (35) that the features of the
kernel learned density $p_{n+1}$ depend on the choice of kernel centers
$\left\\{\hat{x}_{k}^{n+1}\right\\}_{k=1}^{K}$ and the parameters
$\left(\alpha^{n+1},\lambda^{n+1}\right)$, where
$\alpha^{n+1}:=\left(\alpha_{1}^{n+1},\alpha_{2}^{n+1},\cdots,\alpha_{K}^{n+1}\right)^{T}$
and
$\lambda^{n+1}:=\left(\lambda_{1}^{n+1},\lambda_{2}^{n+1},\cdots,\lambda_{K}^{n+1}\right)^{T}$
denote all the weights and covariances of kernels.

Since the state sample points that we generate through the scheme (34)
adaptively follow the filtering density, which provide a good representation
for the target PDF, we choose the kernel centers as a subset of the state
samples $\left\\{\tilde{x}_{i}^{n+1}\right\\}_{i=1}^{N}$. Note that we have
the value of the approximated filtering density
$\tilde{p}\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n+1}}^{M}\right)$ on each
sample point $S_{t_{n+1}}=\tilde{x}_{i}^{n+1}$. Therefore, we can choose state
samples with high density values as kernel centers. To avoid using too many
samples in the mode of the distribution and to capture more features of the
filtering density, we use importance sampling to choose
$\left\\{\hat{x}_{k}^{n+1}\right\\}_{k=1}^{K}$ [24] instead of only using
samples with highest density values.

To determine the parameters of the kernels in kernel learning, we use the
approximated filtering density values
$\left\\{\tilde{p}\left(S_{t_{n+1}}=\tilde{x}_{i}^{n+1}\mid\mathcal{F}_{t_{n+1}}^{M}\right)\right\\}_{i=1}^{N}$
as simulated ”training data”, and we aim to find kernel parameters so that the
kernel learned filtering density $p_{n+1}$ matches the training data. Then, we
implement stochastic gradient descent optimization to determine the parameters
$\alpha^{n+1}$ and $\lambda^{n+1}$. Specifically, we define the loss function
to be minimized as

$\displaystyle F_{\alpha,\lambda}^{n+1}$
$\displaystyle:=\mathbb{E}\left[\left(p_{n+1}\left(S_{t_{n+1}}\right)-\tilde{p}\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n+1}}^{M}\right)\right)^{2}\right]$
$\displaystyle=\mathbb{E}\left[\left(\sum_{k=1}^{K}\alpha_{k}^{n+1}\phi_{k}^{n+1}\left(S_{t_{n+1}}\right)-\tilde{p}\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n+1}}^{M}\right)\right)^{2}\right]$

Since the state $S_{t_{n+1}}$ is represented by spatial samples
$\left\\{\tilde{x}_{i}^{n+1}\right\\}_{i=1}^{N}$, the fully calculated Monte
Carlo simulation for the above loss function is given as

$F_{\alpha,\lambda}^{n+1}\approx\frac{1}{N}\sum_{i=1}^{N}\left(p_{n+1}\left(S_{t_{n+1}}=\tilde{x}_{i}^{n+1}\right)-\tilde{p}\left(S_{t_{n+1}}=\tilde{x}_{i}^{n+1}\mid\mathcal{F}_{t_{n+1}}^{M}\right)\right)^{2}$

Then, for pre-chosen initial estimates $\alpha^{n+1}(0)$ and
$\lambda^{n+1}(0)$, we have the following stochastic gradient descent
iteration to search for the parameters $\alpha^{n+1}$ and $\lambda^{n+1}$ :

$\displaystyle\alpha^{n+1}(j+1)$
$\displaystyle=\alpha^{n+1}(j)-\left.\rho_{\alpha}^{j}\nabla_{\alpha}F_{\alpha,\lambda}^{n+1}\right|_{S_{t_{n+1}}=\bar{x}_{j}},j=0,1,2,\cdots,J-1,$
(36a) $\displaystyle\lambda^{n+1}(j+1)$
$\displaystyle=\lambda^{n+1}(j)-\left.\rho_{\lambda}^{j}\nabla_{\lambda}F_{\alpha,\lambda}^{n+1}\right|_{S_{t_{n+1}}=\bar{x}_{j}},j=0,1,2,\cdots,J-1,$
(36b)

where $\rho_{\alpha}^{j}$ and $\rho_{\lambda}^{j}$ are learning rates for the
parameters $\alpha$ and $\lambda$, respectively, and $J$ is the total number
of iterations corresponding to a stopping criteria. The gradients
$\left.\nabla_{\alpha}F_{\alpha,\lambda}^{n+1}\right|_{S_{t_{n+1}}=\bar{x}_{j}}$
and
$\left.\nabla_{\lambda}F_{\alpha,\lambda}^{n+1}\right|_{S_{t_{n+1}}=\bar{x}_{j}}$
of the cost function $F_{\alpha,\lambda}^{n+1}$ are single-sample
representations of the gradients $\nabla_{\alpha}F_{\alpha,\lambda}^{n+1}$ and
$\nabla_{\lambda}F_{\alpha,\lambda}^{n+1}$ by choosing a specific state sample
$\bar{x}_{j}$ for the state variable $S_{t_{n+1}}$, where the sample
$\bar{x}_{j}$ is picked among the sample set
$\left\\{\tilde{x}_{i}^{n+1}\right\\}_{i=1}^{N}$. As a result, we improve the
estimates for $\alpha^{n+1}$ and $\lambda^{n+1}$ gradually by comparing the
kernel learned filtering density $p_{n+1}$ with approximated filtering density
values on samples $\left\\{\tilde{x}_{i}^{n+1}\right\\}_{i=1}^{N}$.

To use the approximated filtering density values more effectively, instead of
picking samples uniformly from the sample set, we use importance sampling to
choose samples according to their density values. In this way, it is more
likely to consider higher filtering density values in the optimization
procedure, which makes the stochastic gradient descent procedure more
efficient. In this work, we let the covariance matrices for Gaussian type
kernels be diagonal to reduce the dimension of optimization, and note that the
scattered kernel centers also provide covariant features of the target
filtering density.

Resampling random spatial points

In the nonlinear filtering problem, the state equation is a diffusion process.
Thus the random spatial samples propagated through the scheme (34) diffuse
after several estimation steps. As a result, fewer and fewer spatial samples
will remain in high probability regions of the filtering density as we
estimate the target state step-by-step. So, we need to resample spatial points
to make them better represent the filtering density by using the kernel
learned updated filtering density $p_{n+1}$ to generate a set of new spatial
samples, denoted by $\left\\{x_{i}^{n+1}\right\\}_{i=1}^{N}$, to replace the
samples $\left\\{\tilde{x}_{i}^{n+1}\right\\}_{i=1}^{N}$, which follow the
predicted filtering density. We want to point out that the kernel learned
filtering density provides the conditional PDF for the target state over the
entire state space. Therefore, this resampling procedure also allows us to
consider probabilistically insignificant regions. On the other hand, the
filtering density $p_{n+1}$ is essentially a combination of Gaussian kernels.
Hence drawing samples from $p_{n+1}$ is very efficient. For example, to
generate the sample $x_{i}^{n+1}$, we first use importance sampling to pick a
kernel $\phi_{k}^{n+1}$ based on weights of kernels. Since $\phi_{k}^{n+1}$ is
a Gaussian kernel, we can simply draw the sample $x_{i}^{n+1}$ from the
Gaussian distribution $N\left(\hat{x}_{k}^{n+1},\lambda_{k}^{n+1}\right)$.

Remark The state spatial samples used in the backward SDE filter have similar
behavior to the particles in the particle filter method, which roughly
characterizes the filtering density in the nonlinear filtering problem. We
want to emphasize that the backward SDE filter can also approximate the
filtering density values on spatial samples. On the other hand, the particle
filter only utilizes particle positions to construct an empirical distribution
for the target state. Therefore, each spatial sample in the backward SDE
filter carries more information about the state distribution than a particle
in the particle filter. Moreover, since the kernel learned filtering density
is a global continuous approximation for the state distribution, it covers a
wide range in the state space, which can provide more robust/stable
performance for the backward SDE filter.

###  3.3 Stochastic optimization for control process

In this subsection, we combine the numerical schemes for the adjoint FBSDEs
system (10) and the BSDE filter algorithm to formulate an efficient stochastic
optimization algorithm to solve for the optimal control process $u^{\ast}$.

On a time instant $t_{n}\in\Pi_{N_{T}}$, we have

$E[(J^{\ast})_{u}^{{}^{\prime}}(u_{t_{i}}^{l,M}|_{t_{n}})|\mathcal{F}_{t_{n}}^{M}]=\int_{\mathbb{R}^{d}}E[\Psi({t_{i}},S_{t_{i}},u_{t_{i}}^{l,M}|_{t_{n}})|S_{t_{n}}=x]\cdot
p(x|\mathcal{F}_{t_{n}}^{M})dx,$ (37)

where $t_{i}\geq t_{n}$ is a time instant after $t_{n}$.

Then, we use the approximate solutions $(\bar{Y}_{i},\bar{Z}_{i})$ of FBSDEs
from schemes (LABEL:PF-FBSDE:EQ15) to replace
$(\bar{Y}_{t_{i}},\bar{Z}_{t_{i}})$ and the conditional distribution
$p(S_{t_{n}}|\mathcal{F}_{t_{n}}^{M})$ is approximated by the filtering
density $p(S_{t_{n}}|\mathcal{F}_{t_{n}}^{M})$ obtained from the BSDE filter
through the kernel learning method. Then, we can generate S samples from
$p(S_{t_{n}}|\mathcal{F}_{t_{n}}^{M})$ and solve the optimal control
$u_{t_{n}}^{\ast}$ through the following gradient descent optimization
iteration

$\displaystyle u_{t_{i}}^{l+1,M}|_{t_{n}}=u_{t_{i}}^{l,M}|_{t_{n}}$
$\displaystyle-r\frac{1}{S}\sum_{s=1}^{S}E\Big{[}b_{u}^{\top}(t_{i},S_{t_{i}},u_{t_{i}}^{l,M}|_{t_{n}})\bar{Y}_{{i}}+\sigma_{u}^{\top}(t_{i},S_{t_{i}},u_{t_{i}}^{l,M}|_{t_{n}})\bar{Z}_{{i}}$
(38)
$\displaystyle+f_{u}^{\top}(t_{i},S_{t_{i}},u_{t_{i}}^{l,M}|_{t_{n}})|S_{t_{n}}=x_{n}^{(s)}\Big{]}$

Then, the standard Monte Carlo method can approximate expectation
$E\Big{[}\cdot|S_{t_{n}}=x_{n}^{(s)}\Big{]}$ by $\Lambda$ samples:

$\displaystyle u_{t_{i}}^{l+1,M}|_{t_{n}}\approx u_{t_{i}}^{l,M}|_{t_{n}}$
$\displaystyle-r\frac{1}{S}\frac{1}{\Lambda}\sum_{s=1}^{S}\sum_{\lambda=1}^{\Lambda}\Big{[}b_{u}^{\top}(t_{i},S_{t_{i}}^{(\lambda,s)},u_{t_{i}}^{l,M}|_{t_{n}})\bar{Y}_{{i}}+\sigma_{u}^{\top}(t_{i},S_{t_{i}}^{(\lambda,s)},u_{t_{i}}^{l,M}|_{t_{n}})\bar{Z}_{{i}}$
(39)
$\displaystyle+f_{u}^{\top}(t_{i},S_{t_{i}}^{(\lambda,s)},u_{t_{i}}^{l,M}|_{t_{n}})|S_{t_{n}}=x_{n}^{(s)}\Big{]}$

We can see from the above Monte Carlo approximation that in order to
approximate the expectation in one gradient descent iteration step, we need to
generate $S\times\Lambda$ samples. This is even more computationally expensive
when the controlled system is a high-dimensional process.

So, we can also apply the idea of stochastic gradient descent (SGD) to improve
the efficiency of classic gradient descent optimization as we used in the
kernel learning for the BSDE filter. Instead of using the fully calculated
Monte Carlo simulation to approximate the conditional expectation, we use only
one realization of $S_{t_{i}}$ to represent the expectation. For the
conditional distribution of the controlled process, we can use the spatial
samples to describe. So, we have

$\displaystyle
E[(J^{\ast})_{u}^{{}^{\prime}}(u_{t_{i}}^{l,M}|_{t_{n}})|\mathcal{F}_{t_{n}}^{M}]\approx
b_{u}^{\top}(t_{i},S_{t_{i}}^{(\hat{l},\hat{s})},u_{t_{i}}^{l,M}|_{t_{n}})\bar{Y}_{{i}}^{(\hat{l},\hat{s})}+\sigma_{u}^{\top}(t_{i},S_{t_{i}}^{(\hat{l},\hat{s})},u_{t_{i}}^{l,M}|_{t_{n}})\bar{Z}_{{i}}^{(\hat{l},\hat{s})}+f_{u}^{\top}(t_{i},S_{t_{i}}^{(\hat{l},\hat{s})},u_{t_{i}}^{l,M}|_{t_{n}})$
(40)

where $l$ is the iteration index, the index $\hat{l}$ indicates that the
random generation of the controlled process varies among the gradient descent
iteration steps .$S_{t_{i}}^{(\hat{l},\hat{s})}$ indicates a randomly
generated realization of the controlled process with a randomly selected
initial state $S_{t_{n}}^{(\hat{l},\hat{s})}=x_{n}^{\hat{s}}$ from the spatial
sample cloud $\\{x_{n}^{(s)}\\}_{s=1}^{S}$.

Then, we have the following SGD schemes:

$\displaystyle u_{t_{i}}^{l+1,M}|_{t_{n}}=u_{t_{i}}^{l,M}|_{t_{n}}$
$\displaystyle-r\Big{(}b_{u}^{\top}(t_{i},S_{t_{i}}^{(\hat{l},\hat{s})},u_{t_{i}}^{l,M}|_{t_{n}})\bar{Y}_{{i}}^{(\hat{l},\hat{s})}+\sigma_{u}^{\top}(t_{i},S_{t_{i}}^{(\hat{l},\hat{s})},u_{t_{i}}^{l,M}|_{t_{n}})\bar{Z}_{{i}}^{(\hat{l},\hat{s})}$
(41)
$\displaystyle+f_{u}^{\top}(t_{i},S_{t_{i}}^{(\hat{l},\hat{s})},u_{t_{i}}^{l,M}|_{t_{n}})\Big{)}$

$\bar{Y}_{{i}}^{(\hat{l},\hat{s})}$ is the approximate solution $\bar{Y}_{i}$
corresponding to the random sample $S_{t_{i}}^{(\hat{l},\hat{s})}$. And the
path of $S_{t_{i}}^{(\hat{l},\hat{s})}$ is generated as following

$S_{t_{i+1}}^{(\hat{l},\hat{s})}=S_{t_{i}}^{(\hat{l},\hat{s})}+b(t_{i},S_{t_{i}}^{(\hat{l},\hat{s})},u_{t_{i}}^{l,M}|_{t_{n}})\triangle
t_{i}+\sigma(t_{i},S_{t_{i}}^{(\hat{l},\hat{s})},u_{t_{i}}^{l,M}|_{t_{n}})\sqrt{\triangle
t_{i}}\omega_{i}^{(\hat{l},\hat{s})}$ (42)

where $\omega_{i}^{(\hat{l},\hat{s})}\sim N(0,1)$. Then, an estimate for our
desired data driven optimal control at time instant $t_{n}$ is

$\hat{u}_{t_{n}}:=u_{t_{n}}^{L,M}|_{t_{n}}$

The scheme for FBSDE is

$\displaystyle\bar{Y}_{{i}}^{(\hat{l},\hat{s})}=\bar{Y}_{{i+1}}^{(\hat{l},\hat{s})}+\Big{[}b_{x}(t_{i+1},S_{t_{i+1}}^{(\hat{l},\hat{s})},u_{t_{i}}^{l,M}|_{t_{n}})^{\top}\bar{Y}_{{i+1}}^{(\hat{l},\hat{s})}+\sigma_{x}(t_{i+1},S_{t_{i+1}}^{(\hat{l},\hat{s})},u_{t_{i}}^{l,M}|_{t_{n}})^{\top}\bar{Z}_{{i+1}}^{(\hat{l},\hat{s})}$
(43)
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;+f_{x}(t_{i+1},S_{t_{i+1}}^{(\hat{l},\hat{s})},u_{t_{i}}^{l,M}|_{t_{n}})^{\top}\Big{]}\triangle
t_{i}$
$\displaystyle\bar{Z}_{{i}}^{(\hat{l},\hat{s})}=\frac{\bar{Y}_{{i+1}}^{(\hat{l},\hat{s})}\omega_{i}^{(\hat{l},\hat{s})}}{\sqrt{\triangle
t_{i}}}$

Then, we summarize the algorithm in the table below:

Algorithm 1 BSDE-SGD algorithm for data driven feedback control problem

Initialize the spatial sample cloud $\left\\{x_{i}^{0}\right\\}_{i=1}^{N}\sim
p_{0}$, the number of kernels $K$, the learning rates $\rho_{\alpha}$ and
$\rho_{\lambda}$, the number of iterations $L,J\in\mathbb{N}$, the total
number of time steps $N_{T}$.

while $n=0,1,2,\cdots,N_{T},$ do

Initialize an estimated control process
$\\{u_{t_{i}}^{(0,M)}|_{t_{n}}\\}_{i=n}^{N_{T}}$ and a step-size $\rho$

for SGD iteration steps $l=0,1,2,\cdots,L$, do

 _$\bullet$_ Simulate one realization of controlled process
$\\{S_{t_{i+1}}^{(\hat{l},\hat{s})}|_{t_{n}}\\}_{i=n}^{N_{T}-1}$ through
scheme (42) with
$\\{S_{t_{n}}^{(\hat{l},\hat{s})}\\}=x_{n}^{(\hat{s})}\in\\{x_{n}^{(s)}\\}_{s=1}^{S}$;

 _$\bullet$_ Calculate solution
$\\{\bar{Y}_{t_{i}}^{(\hat{l},\hat{s})}|_{t_{n}}\\}_{i=N_{T}}^{n}$ of the
FBSDEs system (10) corresponding to
$\\{S_{t_{i+1}}^{(\hat{l},\hat{s})}|_{t_{n}}\\}_{i=n}^{N_{T}-1}$ through
schemes (43);

 _$\bullet$_ Update the control process to obtain
$\\{u_{t_{i}}^{(l+1,M)}|_{t_{n}}\\}_{i=n}^{N_{T}}$ through scheme (41);

end for

Update PDF of the state

 _$\bullet$_ Propagate samples $\left\\{x_{i}^{n}\right\\}_{i=1}^{N}$ through
the scheme (34) to get $\left\\{\tilde{x}_{i}^{n+1}\right\\}_{i=1}^{N}$.

 _$\bullet$_ Let $Y_{n}^{M_{t_{n}}}=p_{n}$ be the previous kernel learned
filtering density. Solve the timeinverse FBSDEs system (19) for
$Y_{n+1}^{M_{t_{n}}}$ on spatial samples
$\left\\{\tilde{x}_{i}^{n+1}\right\\}_{i=1}^{N}$ through the iterative scheme
(31). $Y_{n+1}^{M_{t_{n}}}$ is the approximation for the predicted filtering
density $\tilde{p}\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n}}^{M}\right)$.

 _$\bullet$_ Incorporate the observational information through Bayesian
inference to get the updated filtering density
$\tilde{p}\left(S_{t_{n+1}}\mid\mathcal{F}_{t_{n+1}}^{M}\right)$ on the
spatial samples $\left\\{\tilde{x}_{i}^{n+1}\right\\}_{i=1}^{N}$.

 _$\bullet$_ Select kernel centers from the spatial samples
$\left\\{\tilde{x}_{i}^{n+1}\right\\}_{i=1}^{N}$ by using the updated
filtering density values on those samples.

 _$\bullet$_ Consider the approximated filtering density values
$\left\\{\tilde{p}\left(S_{t_{n+1}}=\tilde{x}_{i}^{n+1}\mid\mathcal{F}_{t_{n+1}}^{M}\right)\right\\}_{i=1}^{N}$
as training data. Use the optimization procedure (36) to obtain the kernel
learned filtering density $p_{n+1}$ introduced in (35);

 _$\bullet$_ Carry out the resampling procedure to generate new samples
$\left\\{x_{i}^{n+1}\right\\}_{i=1}^{N}$ that follow the kernel learned
filtering density $p_{n+1}$;

end while

##  4 Numerical examples

###  4.1 10D Linear Quadratic Experiment

In this experiment, we want to show that the BSDE-SGD method can be applied to
solve the classic linear quadratic control problem.

Assume A, B, R, K, Q are symmetric, positive definite. The forward process $Y$
and the observation process $M$ is given by

$\displaystyle dY(t)$ $\displaystyle=A(u(t)-r(t))dt+\sigma Bu(t)dW_{t}$ (44)
$\displaystyle dM(t)$ $\displaystyle=sin(Y(t))dt+dB_{t}$

The cost functional is given by

$J[u]=\frac{1}{2}E\left[\int_{0}^{T}\langle
R(Y_{t}-Y^{\ast}_{t}),(Y_{t}-Y^{\ast}_{t})\rangle
dt+\frac{1}{2}\int_{0}^{T}\langle Ku_{t},u_{t}\rangle dt+\frac{1}{2}\langle
QY_{T},Y_{T}\rangle\right]$ (45)

we want to find $J(u^{\ast})=\inf_{u\in\mathcal{U}_{ad}[0,T]}J(u)$.

####  4.1.1 Construction of exact solutions

An interesting fact of such example is that one can construct a time
deterministic exact solution which depend only on $x_{0}$.

By simplifying (45), we have

$J[u]=\frac{1}{2}\int_{0}^{T}(RE[Y_{t}^{T}Y_{t}]-2RY^{\ast
T}_{t}E[Y_{t}]+RY^{\ast T}Y^{\ast}+\langle Ku_{t},u_{t}\rangle
dt)+\frac{1}{2}E[\langle QY_{T},Y_{T}\rangle]$ (46)

Then, we define:

$\displaystyle S_{t}$
$\displaystyle=E[Y_{t}]=E[Y_{0}+\int_{0}^{T}A(u(s)-r(s))ds+\int_{0}^{T}\sigma
Bu(s)dW_{s}]$ (47) $\displaystyle=E[Y_{0}+\int_{0}^{T}A(u(s)-r(s)ds]$

Hence, we see that

$\displaystyle E[Y_{t}^{T}Y_{t}]$
$\displaystyle=E[(Y_{0}+\int_{0}^{t}A(u(s)-r(s))ds+\int_{0}^{t}\sigma
Bu(s)dW_{s})^{2}]$ (48)
$\displaystyle=E[Y_{0}^{T}Y_{0}+\int_{0}^{t}(u(s)-r(s))^{T}A^{T}A(u(s)-r(s)ds+Y_{0}^{T}\int_{0}^{t}A(u(s)-r(s))ds$
$\displaystyle+\int_{0}^{t}(u(s)-r(s))^{T}A^{T}Y_{0}ds]+E[\int_{0}^{t}\sigma^{2}u(s)^{T}B^{T}Bu(s)ds]$
$\displaystyle=S_{t}^{T}S_{t}+\sigma^{2}\int_{0}^{t}u(s)^{T}B^{T}Bu(s)ds$

And (48) is true because all the terms are deterministic in time given
$x_{0}$. Also, we observe that

$\displaystyle E[Y_{T}^{T}Y_{T}]$
$\displaystyle=E[(Y_{0}+\int_{0}^{T}A(u(s)-r(s))ds+\int_{0}^{T}\sigma
Bu(s)dW_{s})^{2}]$ (49)
$\displaystyle=S_{T}^{T}S_{T}+\sigma^{2}\int_{0}^{T}u(s)^{T}B^{T}Bu(s)ds$

As a result, we see that now (46) takes the form:

$\displaystyle J[u]=$
$\displaystyle\frac{1}{2}\int_{0}^{T}(S_{s}^{T}RS_{s}-2S_{s}^{T}RS_{s}^{\ast}+S_{s}^{\ast
T}RS_{s}^{\ast}+u_{s}^{T}(\sigma^{2}B^{T}QB+K)u_{s})ds$ (50)
$\displaystyle+\frac{1}{2}\sigma^{2}\int_{0}^{T}\int_{0}^{t}u_{s}^{T}B^{T}RBu_{s}dsdt+\frac{1}{2}S_{T}^{T}QS_{T}$

by doing a simple integration by part, we have

$\displaystyle J[u]=$
$\displaystyle\frac{1}{2}\int_{0}^{T}(S_{s}^{T}RS_{s}-2S_{s}^{T}RS_{s}^{\ast}+S_{s}^{\ast
T}RS_{s}^{\ast}+u_{s}^{T}(\sigma^{2}B^{T}QB+K)u_{s})ds$ (51)
$\displaystyle+\frac{1}{2}\sigma^{2}\int_{0}^{T}(T-t)u_{s}^{T}B^{T}RBu_{s}ds+\frac{1}{2}S_{T}^{T}QS_{T}$

As a result, we have the following standard deterministic control problem:

$\displaystyle J[u]=$
$\displaystyle\frac{1}{2}\int_{0}^{T}\underbrace{(S_{s}^{T}RS_{s}-2S_{s}^{T}RS_{s}^{\ast}+S_{s}^{\ast
T}RS_{s}^{\ast}+u_{s}^{T}(\sigma^{2}B^{T}QB+K)u_{s}+\sigma^{2}(T-t)u_{s}^{T}B^{T}RBu_{s})}_{2f}ds$
(52) $\displaystyle+\frac{1}{2}S_{T}^{T}QS_{T}$
$\frac{dX_{t}}{dt}=\underbrace{A(u(t)-r(t))}_{b},X_{t_{0}}=X_{0}$ (53)

Then, one can form the following Hamiltonian

$H(x,p,u)=bp+f$ (54)

where $p$ is $\nabla_{x}v$, $v$ is the value function.

Then, to find the optimal control, we have

$\displaystyle\frac{\partial}{\partial u}H=0$ (55)

which is

$\displaystyle Ap+u(\sigma^{2}B^{T}RB(T-t)+(K+\sigma^{2}B^{T}QB))=0$ (56)

So, we got

$\displaystyle
u_{t}=\frac{-Ap(t)}{\sigma^{2}B^{T}RB(T-t)+(K+\sigma^{2}B^{T}QB)}$ (57)

Also, notice that

$\displaystyle\frac{d}{dt}p(t)=-R(X_{t}-X_{t}^{\ast}),p(T)=QX_{T}$ (58)

then,

$\displaystyle p(t)=QX_{T}+\int_{t}^{T}R(X_{s}-X_{s}^{\ast})ds$ (59)

Combine (53),(57),(59) together, we can solve the control of the system.

Set

$\displaystyle\mathbf{A}=\begin{bmatrix}1&0.2&\cdots&0.2&0.2\\\
0.2&1&\cdots&0.2&0.2\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
0.2&0.2&\cdots&1&0.2\\\ 0.2&0.2&\cdots&0.2&1\end{bmatrix}_{10\times 10}$ (60)

Set B, R, K ,Q be identity matrices. To solve (59), let

$\displaystyle X_{t}^{1}-X_{t}^{\ast 1}$ $\displaystyle:=t$ $\displaystyle
X_{t}^{2}-X_{t}^{\ast 2}:=cos(t)$ $\displaystyle X_{t}^{3}-X_{t}^{\ast 3}$
$\displaystyle:=t^{2}$ $\displaystyle X_{t}^{4}-X_{t}^{\ast 4}:=-2\pi sin(2\pi
t)$ $\displaystyle X_{t}^{5}-X_{t}^{\ast 5}$
$\displaystyle:=-\frac{1}{2}e^{(T-t)}$ $\displaystyle X_{t}^{6}-X_{t}^{\ast
6}:=\frac{1}{1+t}$ $\displaystyle X_{t}^{7}-X_{t}^{\ast 7}$
$\displaystyle:=sec^{2}(t)$ $\displaystyle X_{t}^{8}-X_{t}^{\ast
8}:=\frac{1}{1+t^{2}}$ $\displaystyle X_{t}^{9}-X_{t}^{\ast 9}$
$\displaystyle:=1$ $\displaystyle X_{t}^{10}-X_{t}^{\ast{10}}:=e^{(t-T)}$

Then we have

$\displaystyle p(t)=\begin{bmatrix}X_{t}^{1}\\\ X_{t}^{2}\\\ X_{t}^{3}\\\
X_{t}^{4}\\\ X_{t}^{5}\\\ X_{t}^{6}\\\ X_{t}^{7}\\\ X_{t}^{8}\\\ X_{t}^{9}\\\
X_{t}^{10}\end{bmatrix}-\begin{bmatrix}\frac{T^{2}}{2}-\frac{t^{2}}{2}\\\
sin(T)-sin(t)\\\ \frac{T^{3}}{3}-\frac{t^{3}}{3}\\\ cos(2\pi T)-cos(2\pi t)\\\
sinh(T)-sinh(t)\\\ log(1+T)-log(1+t)\\\ tan(T)-tan(t)\\\
arctan(T)-arctan(t)\\\ T-t\\\ 1-e^{(t-T)}\end{bmatrix}$ (61)

Let $\hat{r}(t)$ as

$\displaystyle\hat{r}_{t}^{1}=\frac{-t^{2}/2}{\beta_{t}},\hat{r}_{t}^{2}=\frac{-sin(t)}{\beta_{t}},\hat{r}_{t}^{3}=\frac{-t^{3}}{3\beta_{t}},\hat{r}_{t}^{4}=\frac{-cos(2\pi
t)}{\beta_{t}},\hat{r}_{t}^{5}=\frac{-sinh(t)}{\beta_{t}}$
$\displaystyle\hat{r}_{t}^{6}=\frac{-log(1+t)}{\beta_{t}},\hat{r}_{t}^{7}=\frac{-tan(t)}{\beta_{t}},\hat{r}_{t}^{8}=\frac{-arctan(t)}{\beta_{t}},\hat{r}_{t}^{9}=\frac{-t}{\beta_{t}},\hat{r}_{t}^{10}=\frac{-e^{(t-T)}}{\beta_{t}}$

where $\beta_{t}=(1+\sigma^{2})+\sigma^{2}(T-t)$. then

$r(t)=A\hat{r}(t)$

Then, plug (61) into (57), then we solve (53)

$\displaystyle\frac{dX_{t}}{dt}$ $\displaystyle=A(u(t)-r(t))$ (62)
$\displaystyle=A(-\frac{Ap(t)}{\beta_{t}}-A\hat{r}(t)))$ (63)

we get

$\displaystyle X_{t}=\begin{bmatrix}X_{t}^{1}\\\ X_{t}^{2}\\\ X_{t}^{3}\\\
X_{t}^{4}\\\ X_{t}^{5}\\\ X_{t}^{6}\\\ X_{t}^{7}\\\ X_{t}^{8}\\\ X_{t}^{9}\\\
X_{t}^{10}\end{bmatrix}=\frac{\alpha_{t}A^{2}}{\sigma^{2}}(\begin{bmatrix}\frac{T^{2}}{2}\\\
sin(T)\\\ \frac{T^{3}}{3}\\\ cos(2\pi T)\\\ sinh(T)\\\ log(1+T)\\\ tan(T)\\\
arctan(T)\\\ T\\\ 1\end{bmatrix}-\begin{bmatrix}X_{T}^{1}\\\ X_{T}^{2}\\\
X_{T}^{3}\\\ X_{T}^{4}\\\ X_{T}^{5}\\\ X_{T}^{6}\\\ X_{T}^{7}\\\ X_{T}^{8}\\\
X_{T}^{9}\\\ X_{T}^{10}\end{bmatrix})$ (64)

So, replace $X_{t}$ with $X_{t}^{\ast}$, we get

$\displaystyle X_{t}^{\ast}=\begin{bmatrix}X_{t}^{1,\ast}\\\ X_{t}^{2,\ast}\\\
X_{t}^{3,\ast}\\\ X_{t}^{4,\ast}\\\ X_{t}^{5,\ast}\\\ X_{t}^{6,\ast}\\\
X_{t}^{7,\ast}\\\ X_{t}^{8,\ast}\\\ X_{t}^{9,\ast}\\\
X_{t}^{10,\ast}\end{bmatrix}=\begin{bmatrix}t\\\ cos(t)\\\ t^{2}\\\ -2\pi
sin(2\pi t)\\\ -\frac{1}{2}e^{(T-t)}\\\ \frac{1}{1+t}\\\ sec^{2}(t)\\\
\frac{1}{1+t^{2}}\\\ 1\\\
e^{(t-T)}\end{bmatrix}+\frac{\alpha_{t}A^{2}}{\sigma^{2}}(\begin{bmatrix}\frac{T^{2}}{2}\\\
sin(T)\\\ \frac{T^{3}}{3}\\\ cos(2\pi T)\\\ sinh(T)\\\ log(1+T)\\\ tan(T)\\\
arctan(T)\\\ T\\\ 1\end{bmatrix}-\begin{bmatrix}X_{T}^{1}\\\ X_{T}^{2}\\\
X_{T}^{3}\\\ X_{T}^{4}\\\ X_{T}^{5}\\\ X_{T}^{6}\\\ X_{T}^{7}\\\ X_{T}^{8}\\\
X_{T}^{9}\\\ X_{T}^{10}\end{bmatrix})$ (65)

where $X_{T}^{i}$ can be obtained from the system (64) by letting $t=T$, and

$\alpha_{t}=ln\frac{1+\sigma^{2}+\sigma^{2}T}{(1+\sigma^{2})+\sigma^{2}(T-t)}$

Then, to find an the exact form by following the trajectory of $y_{t}$ in this
setup, one will have to solve the following Coupled forward-backward ODE.

$\displaystyle\frac{dX_{t}}{dt}$ $\displaystyle=A(u_{t}-r_{t})\ \ \ \ \
x_{n}=y_{t_{n}}$ (66) $\displaystyle\frac{dp(t)}{dt}$
$\displaystyle=X_{t}-X^{*}_{t},\ \ \ p_{T}=x^{t_{n},y_{t_{n}}}_{T}$ (67)

with $u_{t}=-Ap_{t}\big{/}\big{(}\sigma^{2}(T-t)+(1+\sigma^{2})\big{)}$. As a
result, we have

$\displaystyle\frac{dX_{t}}{dt}$
$\displaystyle=A(-Ap_{t}\big{/}\big{(}\sigma^{2}(T-t)+(1+\sigma^{2})\big{)}-r_{t})\
\ \ \ \ x_{n}=y_{t_{n}}$ (68) $\displaystyle\frac{dp(t)}{dt}$
$\displaystyle=X_{t}-X^{*}_{t},\ \ \ p_{T}=x^{t_{n},y_{t_{n}}}_{T}$ (69)

That is, we need to solve the above coupled FBODE. Then, seeing that
$p_{t}=x^{t_{n},y_{t_{n}}}_{T}+\int^{T}_{t_{n}}X_{s}-X^{*}_{s}ds$. Writing
$a_{t}:=1/\big{(}\sigma^{2}(T-t)+(1+\sigma^{2})\big{)}$, we have

$\displaystyle\frac{dX_{t}}{dt}$
$\displaystyle=A\big{(}-a_{t}AX_{T}-a_{t}A\int^{T}_{t_{n}}(X_{s}-X^{*}_{s})ds-
r_{t}\big{)},\ \ \ \ \ x_{t_{n}}=y_{t_{n}}$ (70)

To solve (70) numerically, we do a numerical discretization:

$\displaystyle x_{t_{n+1}}-x_{t_{n}}$
$\displaystyle=-a_{t_{n}}A^{2}X_{T}\Delta t-a_{t_{n}}(\Delta
t)^{2}A^{2}\sum^{N-1}_{i=n}(X_{t_{i}}-X^{*}_{t_{i}})-Ar_{t}$
$\displaystyle\Rightarrow Ar_{t}-a_{t_{n}}(\Delta
t)^{2}A^{2}\sum^{N-1}_{i=n}X^{*}_{t_{i}}$
$\displaystyle=X_{t_{n}}-X_{t_{n+1}}-a_{t_{n}}(\Delta
t)^{2}A^{2}\sum^{N-1}_{i=n}X_{t_{i}}-a_{t_{n}}A^{2}X_{T},\ \ \
x_{t_{n}}=y_{t_{n}}$ (71)

We can put (71) into a large linear system, and solve it numerically.

####  4.1.2 Result

Set iteration $L=10^{4},\sigma=0.1$, $T=1$, $X_{0}=0$ we have the following
results in figure 1 and 2. We demonstrate the effectiveness of our data driven
feedback control algorithm by comparing with the analytical solution, where
the analytical optimal control is calculated from (57) with the exact
controlled state and we let the analytical optimal control drive its
corresponding controlled system through (44) to get the analytical controlled
process. The estimated optimal control and its corresponding controlled
process are calculated by our data driven feedback control Algorithm based on
fully nonlinear observations. figure 1 shows that the estimated control is
very close to the analytical optimal control. figure 2 shows the controlled
process governed by our data driven feedback control well aligns with the real
state of the analytical controlled process.

Figure 1: Estimated control vs True optimal control Figure 2: Estimated
trajectory vs True trajectory

####  4.1.3 Testing convergence

To test the convergence of our algorithm, we increase the number of
iterations, we got the result in (figure 3(a)), where the error between the
estimated control and the analytical control is the average of
$||u^{est}-u^{\ast}||_{2}$ of 50 repeated experiments. The result shows that
when the number of iterations in the stochastic gradient descent algorithm
increases, the error decreases and converges. This means we can get a better
result by increasing the number of iterations. At the same time, we also got
the result in (figure 3(b)), where the error between the controlled process
governed by our data driven feedback control and the analytical controlled
process is the average of $||X^{est}-X^{ture\;trajectory}||_{2}$ of the same
50 repeated experiments. It tells that the BSDE filter provides a very
accurate and stable estimation of the true state X.

(a) Control u

(b) State X

Figure 3: Error of BSDE-SGD Estimation on Control u and State X

###  4.2 PDE Experiment

In this experiment, we want to show that the BSDE-SGD method can be applied to
solve the PDE with control terms. The idea is we use a well-studied scheme to
discretize the PDE where we get a high-dimensional ODE system. Then we can add
a diffusion term to convert it into a high-dimensional SDE system. We also
assume that we do not have the true state and we can only have observation
data. Finally, we can solve this stochastic optimal feedback control problem.

First, the problem is to find the optimal control of the 1-D heat equation:

$\displaystyle
P_{t}=aP_{xx}+bP_{x}+cu(x,t)+df(x,t),\;\;\;\;(x,t)\in[0,1]\times[0,1]$ (72)
$\displaystyle P(x,0)=10x(1-x)\frac{10}{(1+0)}e^{-\frac{x-0}{(1+0)^{2}}}$
$\displaystyle f(x,t)=x(1-x)t(1-t)$ $\displaystyle P(0,t)=0$ $\displaystyle
P(1,t)=0$

The cost function J is defined by

$J[u]=\frac{1}{2}E\left[\int_{0}^{T}\langle RP(t),P(t)\rangle
dt+\int_{0}^{T}\langle Qu(t),u(t)\rangle dt+\langle KP(T),P(T)\rangle\right]$
(73)

The discretized system is given by:

$\frac{d}{dt}P=AP+Bu+Cf$ (74)

using forward-time for $P_{t}$, central-space for $P_{xx}$ and backward-space
for $P_{x}$, the matrices $A$,$B$, and $C$ are given by

$\displaystyle\mathbf{A}=\frac{a}{h^{2}}\begin{bmatrix}-2&1&&&\\\ 1&-2&1&&\\\
&\ddots&\ddots&\ddots&\\\ &&1&-2&1\\\ &&&1&-2\end{bmatrix}_{n\times
n}+\frac{b}{h}\begin{bmatrix}1&&&&\\\ -1&1&&&\\\ &\ddots&\ddots&&\\\
&&-1&1&\\\ &&&-1&1\end{bmatrix}_{n\times n}$ (75)

$\displaystyle\mathbf{B}=\begin{bmatrix}0\\\ c_{1}\\\ \vdots\\\ c_{n-1}\\\
0\end{bmatrix}_{n\times 1},\mathbf{C}=\begin{bmatrix}0\\\ d_{1}\\\ \vdots\\\
d_{n-1}\\\ 0\end{bmatrix}_{n\times 1}$ (76)

For the system of (74), the Riccati equations are given by:

$\frac{d}{dt}G(t)=-A^{T}G(t)-G(t)A-G(t)BR^{-1}B^{T}G(t)+Q$ (77)

The optimal control $u$ for the above stochastic optimal control problem can
be derived analytically as

$u^{\ast}=(-R^{-1}B^{T}G(t)P(t)-df(t))/c$ (78)

To ensure the discretization is stable, we set $a=b=5\times 10^{-5}$. We also
set $c=d=1$, $R$, $Q$, and $K$ be identity matrix.

To apply the BSDE-SGD method, we can add the diffusion term to the system (74)
and the observation operator is fully nonlinear:

$\displaystyle dP_{t}$ $\displaystyle=(AP+Bu+Cf)dt+\sigma dW_{t}$ (79)
$\displaystyle dM_{t}$ $\displaystyle=sin(P_{t})dt+\eta dB_{t}$

Then we can use the BSDE-SGD method to the stochastic optimal feedback control
problem giving the cost function J in (73).

Set iteration $L=10^{4},\sigma=0.1,\eta=0.1$, and number of nodes on space
$n=21$. We have the following results. From Figure 4, we see that the
estimated control is very accurate compared with the analytical solution of
control. From Figure 5, we also see that the estimated heat distribution
governed by our data driven feedback control well aligns with the analytical
real heat distribution. We demonstrate the effectiveness of our data driven
feedback control algorithm of solving a PDE control problem by converting it
into a high-dimensional SDE system governed by the control.

(a) t=0

(b) t=0.2

(c) t=0.4

(d) t=0.6

(e) t=0.8

(f) t=1

Figure 4: BSDE-SGD Estimation VS Optimal Control

(a) t=0

(b) t=0.2

(c) t=0.4

(d) t=0.6

(e) t=0.8

(f) t=1

Figure 5: BSDE-SGD Estimation VS True State

###  4.3 3D Dubins Airplane maneuvering problem

In this example, we solve a Dubins airplane maneuvering problem. The
controlled process is described by the following nonlinear controlled
dynamics:

$\displaystyle dS_{t}$ $\displaystyle=\begin{bmatrix}dX_{t}\\\ dY_{t}\\\
dZ_{t}\end{bmatrix}=\begin{bmatrix}v\cos(\theta_{t})\cos(\phi_{t})dt+\sigma
dW_{1t}\\\ v\cos(\theta_{t})\sin(\phi_{t})dt+\sigma dW_{2t}\\\
v\sin(\theta_{t})+\sigma dW_{3t}\end{bmatrix}$ (80) $\displaystyle
d\theta_{t}$ $\displaystyle=u_{t}+\sigma^{2}dB_{1t}$ (81) $\displaystyle
d\phi_{t}$ $\displaystyle=p_{t}+\sigma^{2}dB_{2t}$ (82) $\displaystyle dM_{t}$
$\displaystyle=[arctan(\frac{X_{t}+3}{Y_{t}+2}),arctan(\frac{X_{t}-2}{Y_{t}+2}),arctan(\frac{Z_{t}-2}{Y_{t}+2})]^{T}+\eta_{t}$
(83)

where the pair $(X,Y,Z)$ gives the position of a airplane-like robot moving in
the 3D plane. $\theta$ is the angle that controls the moving direction of the
airplane vertically, which is governed by the control action $u_{t}$. $\phi$
is the angle that controls the moving direction of the airplane horizontally
which is governed by the control action $p_{t}$. $\sigma$ is the noise that
perturbs the motion and control actions. Assume that we do not have direct
observations on the robot. Instead, we use two detectors located on different
observation platforms at $(-3,-2,2)$ and $(2,-2,-2)$ to collect the bearing
angles of the target robot as indirect observations. So, we have the
observation process $M_{t}$ Given the expected path $S^{\ast}$, the airplane
should follow it and arrive at the terminal position on time. The performance
cost functional based on observational data that we aim to minimize is defined
as:

$J[u]=E\left[\frac{1}{2}\int_{0}^{T}\langle
R(S_{t}-S^{\ast}_{t}),(S_{t}-S^{\ast}_{t})\rangle
dt+\frac{1}{2}\int_{0}^{T}\langle Ku_{t},u_{t}\rangle dt+\langle
Q(S_{T}-S^{\ast}_{T}),(S_{T}-S^{\ast}_{T})\rangle\right]$ (84)

In our numerical experiments, we let the airplane start from
$(X_{0},Y_{0},Z_{0})=(0,0.5,0)$ to $(X_{T},Y_{T},Z_{T})=(0,0.5,1)$. The
expected path $S^{\ast}_{t}$ is

$X(t)=0.5\sin(2\pi t),\;Y(t)=0.5\cos(2\pi t),\;Z(t)=t$

Other settings are $T=1$, $\Delta t=0.02$ i.e. $N_{T}=50$, $\sigma=0.1$,
$\eta_{t}\sim N(0,0.1)$, $L=1000$, $K=1$ and the initial heading direction is
$\theta=0$ and $\phi=arctan(1/2\pi)$. To emphasize the importance of following
the expected path and arriving at the target location at the terminal time,
let $Q=20$ and $R=40$.

Figure 6: BSDE-SGD Estimate VS Designed trajectory Figure 7: BSDE-SGD Estimate
VS Designed trajectory in Each Dimension (X, Y, Z from left to right)

We can see that the observation is highly non-linear. But the estimation of
the BSDE-SDG method gives the result where the airplane moves to the designed
terminal position. At the meantime, it also follows the designed trajectory
closely.

####  4.3.1 Compare with the Particle Filter Method

In this part, we will add a big jump and huge noise to the Airplane and the
system to simulate the situation where the airplane is suddenly in front of a
thundercloud. In this case, the airplane has to dive quickly to avoid crossing
the cloud. In the meantime, radars on the ground will receive huge noise due
to the weather. In the figure 8, we see after it happens, the particle filter
method fails to track to position of the airplane. But figure 9 shows that the
BSDE filter works well to track the airplane when a big jump happens. It also
shows that when the filter works well, our data driven feedback control can
”drag” the airplane back to the designed path and move to the target position.

Figure 8: PF-SGD Estimate VS Designed trajectory in Each Dimension (X, Y, Z
from left to right) Figure 9: BSDE-SGD Estimate VS Designed trajectory in Each
Dimension (X, Y, Z from left to right)

####  4.3.2 Compare with the Dynamic Programming Method

In this section, we will compare our algorithm with the Dynamic Programming
method. In figure 6, we have shown the result of our algorithm where we used
1,000 iterations which cost 67 seconds. To compare with the Dynamic
Programming method, we set the 3 by 3 grid for the control, $T=1$, $N_{T}=8$.
Note that we already know the estimated control from our algorithm, so the 3
by 3 grid for the control is very close to its best choice. But even with this
”trick”, it still takes 627.5 seconds to get the result in figure 10 below.

Figure 10: DP Estimate VS Designed trajectory in Each Dimension (X, Y, Z from
left to right)

From figure 10, we see the Dynamic Programming method works, but the result is
very rough. If we when to get a better result, we will need to increase the
time partition and set a control grid that covers more points. But by just
setting $N_{T}=50$, it will take $6\times 10^{42}$ seconds to finish.

## References

  * [1] C. Andrieu, A. Doucet, and R. Holenstein.  Particle markov chain monte carlo methods.  J. R. Statist. Soc. B, 72(3):269–342, 2010.
  * [2] R. Archibald, F. Bao, and J. Yong.  A stochastic gradient descent approach for stochastic optimal control.  East Asian Journal on Applied Mathematics, 10(4):635–658, 2020\.
  * [3] R. Archibald, F. Bao, J. Yong, and T. Zhou.  An efficient numerical algorithm for solving data driven feedback control problems.  Journal of Scientific Computing, 85(51), 2020.
  * [4] Richard Archibald, Feng Bao, Yanzhao Cao, and Hui Sun.  Numerical analysis for convergence of a sample-wise backpropagation method for training stochastic neural networks.  SIAM Journal on Numerical Analysis, 62(2):593–621, 2024.
  * [5] Richard Archibald, Feng Bao, and Jiongmin Yong.  A stochastic maximum principle approach for reinforcement learning with parameterized environment.  Journal of Computational Physics, 488:112238, 2023.
  * [6] Richard Archibald, Feng Bao, Jiongmin Yong, and Tao Zhou.  An efficient numerical algorithm for solving data driven feedback control problems.  Journal of Scientific Computing, 85(2):1–27, 2020.
  * [7] F. Bao, Y. Cao, and P. Maksymovych.  Backward sde filter for jump diffusion processes and its applications in material sciences.  Communications in Computational Physics, 27:589–618, 2020.
  * [8] F. Bao, Y. Cao, and J. Yong.  Data informed solution estimation for forward backward stochastic differential equations.  Analysis and Applications, to appear, 2020.
  * [9] Feng Bao, Yanzhao Cao, and Xiaoying Han.  Forward backward doubly stochastic differential equations and optimal filtering of diffusion processes.  Communications in Mathematical Sciences, 18(3):635–661, 2020.
  * [10] Feng Bao, Yanzhao Cao, Amnon Meir, and Weidong Zhao.  A first order scheme for backward doubly stochastic differential equations.  SIAM/ASA J. Uncertain. Quantif., 4(1):413–445, 2016.
  * [11] Feng Bao, Yanzhao Cao, Amnon Meir, and Weidong Zhao.  A first order scheme for backward doubly stochastic differential equations.  SIAM/ASA Journal on Uncertainty Quantification, 4(1):413–445, 2016\.
  * [12] Feng Bao, Yanzhao Cao, Clayton Webster, and Guannan Zhang.  A hybrid sparse-grid approach for nonlinear filtering problems based on adaptive-domain of the Zakai equation approximations.  SIAM/ASA J. Uncertain. Quantif., 2(1):784–804, 2014.
  * [13] Feng Bao, Yanzhao Cao, and Weidong Zhao.  Numerical solutions for forward backward doubly stochastic differential equations and zakai equations.  Visualization of Mechanical Processes: An International Online Journal, 1(4), 2011.
  * [14] Feng Bao and Vasileios Maroulas.  Adaptive meshfree backward SDE filter.  SIAM J. Sci. Comput., 39(6):A2664–A2683, 2017.
  * [15] R. Bellman.  Dynamic Programming.  1957\.
  * [16] Léon Bottou and Olivier Bousquet.  The tradeoffs of large scale learning.  Advances in neural information processing systems, 20, 2007.
  * [17] Charalambos D. Charalambous and Robert J. Elliott.  Classes of nonlinear partially observable stochastic optimal control problems with explicit optimal control laws.  SIAM J. Control Optim., 36(2):542–578, 1998.
  * [18] A. J. Chorin and X. Tu.  Implicit sampling for particle filters.  Proc. Nat. Acad. Sc. USA, 106:17249–17254, 2009.
  * [19] A. Doucet, N. de Freitas, and N. Gordon.  Sequential Monte Carlo Methods in Practice.  Springer, New York, 2001.
  * [20] G. Evensen.  The ensemble Kalman filter for combined state and parameter estimation: Monte Carlo techniques for data assimilation in large systems.  IEEE Control Syst. Mag., 29(3):83–104, 2009.
  * [21] Xiaobing Feng, Roland Glowinski, and Michael Neilan.  Recent developments in numerical methods for fully nonlinear second order partial differential equations.  SIAM Rev., 55(2):205–267, 2013.
  * [22] Wendell H. Fleming and Étienne Pardoux.  Optimal control for partially observed diffusions.  SIAM J. Control Optim., 20(2):261–285, 1982.
  * [23] Bo Gong, Wenbin Liu, Tao Tang, Weidong Zhao, and Tao Zhou.  An efficient gradient projection method for stochastic optimal control problems.  SIAM J. Numer. Anal., 55(6):2982–3005, 2017.
  * [24] Neil J Gordon, David J Salmond, and Adrian FM Smith.  Novel approach to nonlinear/non-gaussian bayesian state estimation.  In IEE proceedings F (radar and signal processing), volume 140, pages 107–113. IET, 1993.
  * [25] N.J Gordon, D.J Salmond, and A.F.M. Smith.  Novel approach to nonlinear/non-gaussian bayesian state estimation.  IEE PROCEEDING-F, 140(2):107–113, 1993.
  * [26] U. G. Haussmann.  On the existence of optimal controls for partially observed diffusions.  SIAM J. Control Optim., 20(3):385–407, 1982.
  * [27] Thomas Hofmann, Bernhard Schölkopf, and Alexander J Smola.  Kernel methods in machine learning.  2008\.
  * [28] Kai Kang, Vasileios Maroulas, Ioannis Schizas, and Feng Bao.  Improved distributed particle filters for tracking in a wireless sensor network.  Comput. Statist. Data Anal., 117:90–108, 2018.
  * [29] Peter E. Kloeden and Eckhard Platen.  Stochastic Differential Equations, pages 103–160.  Springer Berlin Heidelberg, Berlin, Heidelberg, 1992.
  * [30] T Kong, W Zhao, and T. Zhou.  High order numerical schemes for second order fbsdes with applications to stochastic optimal control.  Communications in Computational Physics, 21:808–834, 2017.
  * [31] Xin Li, Feng Bao, and Kyle Gallivan.  A drift homotopy implicit particle filter method for nonlinear filtering problems.  Discrete and Continuous Dynamical Systems - S, 15(4):727–746.
  * [32] Étienne Pardoux and Shige Peng.  Backward doubly stochastic differential equations and systems of quasilinear spdes.  Probability Theory and Related Fields, 98(2):209–227, 1994.
  * [33] Etienne Pardoux and Shige Peng.  Backward stochastic differential equations and quasilinear parabolic partial differential equations.  In Stochastic Partial Differential Equations and Their Applications: Proceedings of IFIP WG 7/1 International Conference University of North Carolina at Charlotte, NC June 6–8, 1991, pages 200–217. Springer, 2005.
  * [34] S. Peng.  A general stochastic maximum principle for optimal control problems.  SIAM J. Control Optim., pages 966–979, 1990.
  * [35] Shige Peng.  Backward stochastic differential equation, nonlinear expectation and their applications.  In Proceedings of the International Congress of Mathematicians 2010 (ICM 2010) (In 4 Volumes) Vol. I: Plenary Lectures and Ceremonies Vols. II–IV: Invited Lectures, pages 393–432. World Scientific, 2010.
  * [36] Michael K. Pitt and Neil Shephard.  Filtering via simulation: auxiliary particle filters.  J. Amer. Statist. Assoc., 94(446):590–599, 1999.
  * [37] Herbert Robbins and Sutton Monro.  A stochastic approximation method.  The annals of mathematical statistics, pages 400–407, 1951.
  * [38] I. Sato and H. Nakagawa.  Convergence analysis of gradient descent stochastic algorithms.  Proceedings of the 31st International Conference on Machine Learning, pages 982–990, 2014.
  * [39] A. Shapiro and Y. Wardi.  Convergence analysis of gradient descent stochastic algorithms.  Journal of Optimization Theory and Applications, pages 439–454, 1996.
  * [40] C. Snyder, T. Bengtsson, P. Bickel, and J. Anderson.  Obstacles to high-dimensional particle filtering.  Mon. Wea. Rev., 136:4629–4640, 2008.
  * [41] Shanjian Tang.  The maximum principle for partially observed optimal control of stochastic differential equations.  SIAM J. Control Optim., 36(5):1596–1617, 1998.
  * [42] Xin T. Tong, Andrew J. Majda, and David Kelly.  Nonlinear stability and ergodicity of ensemble based Kalman filters.  Nonlinearity, 29(2):657–691, 2016.
  * [43] Guangchen Wang, Zhen Wu, and Jie Xiong.  An introduction to optimal control of FBSDE with incomplete information.  SpringerBriefs in Mathematics. Springer, Cham, 2018.
  * [44] Jiongmin Yong and Xun Yu Zhou.  Stochastic controls, volume 43 of Applications of Mathematics (New York).  Springer-Verlag, New York, 1999.  Hamiltonian systems and HJB equations.
  * [45] Moshe Zakai.  On the optimal filtering of diffusion processes.  Z. Wahrscheinlichkeitstheorie und Verw. Gebiete, 11:230–243, 1969\.
  * [46] Jianfeng Zhang.  A numerical scheme for BSDEs.  Ann. Appl. Probab., 14(1):459–488, 2004.
  * [47] Weidong Zhao, Yu Fu, and Tao Zhou.  New kinds of high-order multistep schemes for coupled forward backward stochastic differential equations.  SIAM J. Sci. Comput., 36(4):A1731–A1751, 2014.
  * [48] Weidong Zhao, Tao Zhou, and Tao Kong.  High order numerical schemes for second-order fbsdes with applications to stochastic optimal control.  Communications in Computational Physics, 21(3):808–834, 2017.