# An Online Method for the Data Driven Stochastic Optimal Control Problem with
Unknown Model Parameters

Richard Archibald  Computational Science and Mathematics Division, Oak Ridge
National Laboratory, Oak Ridge, Tennessee. Feng Bao  Department of
Mathematics, Florida State University, Tallahassee, Florida,
(bao@math.fsu.edu). Jiongmin Yong  Department of Mathematics, University of
Central Florida, Orlando, Florida.

###### Abstract

In this work, an efficient sample-wise data driven control solver will be
developed to solve the stochastic optimal control problem with unknown model
parameters. A direct filter method will be applied as an online parameter
estimation method that dynamically estimates the target model parameters upon
receiving the data, and a sample-wise optimal control solver will be provided
to efficiently search for the optimal control. Then, an effective overarching
algorithm will be introduced to combine the parameter estimator and the
optimal control solver. Numerical experiments will be carried out to
demonstrate the effectiveness and the efficiency of the sample-wise data
driven control method.

##  1 Introduction

In this paper, we develop an efficient online method for solving the data
driven stochastic optimal control problem with unknown model parameters. In
the classic stochastic optimal control problem, we want to find an optimal
“control process” that controls a stochastic dynamical system (called the
“state process”) to meet certain optimality conditions – such like minimizing
a cost function. In many practical applications (e.g., smart grids [31, 36],
nano-phase materials [26, 18, 34], and mathematical biology [11]), the
controlled state model often contains unknown parameters. Therefore, we need
to estimate the parameters first and then use the estimated parameters to
design the optimal control. When parameter estimation and stochastic optimal
control are combined together, the classic optimal control problem becomes the
“data driven optimal control”, which is an important research topic that
requires uncertainty quantification. The concept of “data driven” in our
control problem refers to an extra parameter estimation procedure [29], and
the estimated parameters will be the “driver” that guides the design of the
optimal control. Extensive studies have been carried out to solve data driven
optimal control problems with time-invariant linear controlled systems [40,
46, 1]. In this work, we will focus on data driven methods for more general
nonlinear control problems.

As a key research topic in mathematical modeling, parameter estimation has
been extensively studied. The efforts in solving the parameter estimation
problem aim to collect observational data about model states and use the data
information to estimate parameters. There are two major approaches to estimate
unknown parameters: the optimization based point estimation and the Bayesian
estimation [45, 35, 33, 24], and both approaches can be implemented in either
the offline manner or the online manner. For the offline parameter estimation,
we use all the observational data that we have ever collected to estimate the
unknown parameters. On the other hand, online parameter estimation methods
dynamically estimate and update the parameters based on the newly received
data. Offline parameter estimation methods are well-studied, and they
typically provide very accurate estimation results. However, the offline
parameter estimation also has several disadvantages. First of all, when
observations form a massive data set, large-scale data storage and analytics
are needed, which often makes the offline parameter estimation prohibitive in
practice. Especially, when a set of new observational data are received,
offline methods need to go through all the historic data again to incorporate
the incoming data, which can be very expensive. Secondly, although we often
assume that the unknown parameters are static, it is also possible that the
true values of the parameters change unexpectedly. For example, when
estimating parameters that determine physical models in scientific
experiments, the true values of the parameters may change unexpectedly due to
sudden environmental changes or some unaware difference in the experimental
set-up [26, 18]. In this case, offline parameter estimation methods are not
suitable since the full data set may be generated by different true
parameters.

The most widely used online parameter estimation method is the augmented
filter [20, 28, 30]. The main idea of the augmented filter is to construct an
augmented state-parameter process and let state estimation guide parameter
estimation under the optimal filtering framework, which is a mathematical
model that provides the best estimate for the state of a hidden stochastic
dynamical system based on observational data [23, 32]. The drawback of the
augmented filter is that the indirect parameter estimation strategy may not
always be efficient enough in terms of usage of data. In fact, most effort in
the optimal filtering procedure in the augmented filter would contribute to
state estimation since the observations only describe the state, and the data
only give direct evidence on the accuracy of state estimation.

In this work, we introduce a novel direct filter method to estimate the
unknown parameters in the controlled state model. As an optimal filtering type
approach, the direct filter defines a parameter process to model the unknown
parameters. The novelty of the direct filter method is that we treat the state
model as a dynamically designed “observation function”, which provides partial
noisy observations on parameters. Then, we apply the particle filter method,
which is a state-of-the-art optimal filtering method, to approximate the
conditional probability density function (PDF) of the target parameters by
using a set of random samples (called “particles”). Then the estimated
parameters are characterized by the parameter particles that describe the PDF
of the parameter variables [3]. The major difference between the direct filter
and the augmented filter is that the direct filter does not use state
estimation to drive parameter estimation. Instead, it utilizes the dynamical
state model to “project” observational information to the parameter space,
which would directly use observational data to estimate the unknown
parameters.

The main contribution of this paper is to construct an effective overarching
algorithm that dynamically combines the direct filter as an online parameter
estimation method with an efficient stochastic optimal control solver and let
the estimated parameters drive the design of the optimal control. Therefore,
besides the direct filter, the stochastic optimal control solver is another
important component in our data driven control algorithm. There are two well-
known approaches to solve the stochastic optimal control problem: the dynamic
programing principle and the stochastic maximum principle [16, 21, 38]. 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 [47]. 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 [22, 44], 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 [49, 48]. Since we aim
to dynamically design control actions based on estimated parameters,
efficiency of the algorithm is very important.

To achieve higher efficiency, we adopt a novel concept that the random samples
used to describe the solution of the adjoint BSDE can be considered as “data
samples”, and we shall provide an efficient sample-wise optimal control
solver. The main theme of the sample-wise optimal control solver is that we
only use one realization of random sample (or a small batch of samples) to
represent the solution of the BSDE and derive numerical approximation for the
solution on the single-realization of random sample. Although the solution of
the BSDE is a random variable, which should be characterized by a set of
statistically selected samples, the rationale that allows us to solve the BSDE
by using only one realization of random sample is that the solution of the
adjoint BSDE is used to formulate the gradient process, which will be used to
carry out gradient descent optimization for the optimal control, and obtaining
a complete approximation for the solution of the adjoint BSDE is not needed in
the stochastic optimal control problem [4]. In this connection, the
methodology of applying stochastic approximation in the stochastic gradient
descent optimization [41, 42] can justify the application of sample-wise
approximation for the solution of the BSDE in the stochastic optimal control
solver. As a result, we avoid solving the adjoint BSDE on a very large number
of random samples, and we transfer the computational cost of completely
solving BSDEs to the optimization procedure of searching for the desired
optimal control, which would improve the efficiency of the stochastic optimal
control solver.

By combining the direct filter with the sample-wise stochastic optimal control
solver, in this paper we shall construct an efficient sample-wise data driven
control solver. Upon receiving a set of new observational data, we first use
the direct filter to estimate/update the unknown parameters in the state
model. Then, we use the sample-wise stochastic optimal control solver to
design control actions based on the state model determined by the estimated
model parameters. To incorporate more data information, which is carried by
the particles that we use to describe the parameter variables in the direct
filter, we treat those particles as “data samples” as well. Therefore, instead
of using the mean value of the particles as a point estimate for the
parameter, at each iteration step in the stochastic gradient descent
optimization, we randomly select one particle in the parameter particle cloud
to represent the corresponding estimated model parameter. In this way, all the
particles are considered in the iteration procedure for the optimal control,
which makes our optimal control solver sufficiently data-driven.

The rest of this paper is organized as follows: In Section 2, we introduce the
data driven stochastic optimal control problem with unknown parameters in the
state model, and we will provide preliminary knowledge that is needed in the
derivation of numerical algorithms. In Section 3, we shall provide detailed
discussions to derive our sample-wise data driven control solver. Numerical
experiments that demonstrate the effectiveness and efficiency of our algorithm
will be presented in Section 4. Finally, we will conclude our research in
Section 5.

##  2 Problem setting

In this section, we will first introduce the data driven stochastic optimal
control problem with unknown parameters. Then, we shall provide the
mathematical foundations that we need to derive the numerical methods for
solving the data driven control problem.

In the complete filtered probability space
$(\Omega,\mathcal{F},\mathbb{F}^{W},\mathbb{P})$ with
$\mathbb{F}^{W}:=\\{\mathcal{F}_{s}^{W}\\}_{s\geq 0}$ being its natural
filtration augmented by all the $\mathbb{P}$-null sets in $\mathcal{F}$, we
consider the following controlled dynamical system in the form of a stochastic
differential equation (SDE)

$dX^{\alpha}_{t}=b(t,X^{\alpha}_{t},u_{t},\alpha)dt+\sigma(t,X^{\alpha}_{t},u_{t})dW_{t},\hskip
10.00002ptt\in[0,T],$ (1)

where the vector $\alpha\in{\cal A}$ with ${\cal A}\subseteq\mathbb{R}^{q}$
being a set of unknown model parameters; $u_{t}$ is the control process valued
in some set $U\subset\mathbb{R}^{m}$; $X^{\alpha}_{t}\in\mathbb{R}^{d}$ is the
controlled process driven by the parameter $\alpha$ and the control process
$u$, which is usually called the “state process”;
$b:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\times\mathbb{R}^{q}\rightarrow\mathbb{R}^{d}$
and
$\sigma:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{d\times
d}$ are suitable maps called drift and diffusion, respectively;
$W:=\\{W_{t}\\}_{t\geq 0}$ is a standard $d$-dimensional Brownian motion in
$(\Omega,\mathcal{F},\mathbb{F}^{W},\mathbb{P})$.

Denote

$\mathcal{U}[0,T]:=\Big{\\{}u:[0,T]\times\Omega\rightarrow U\big{|}\text{ $u$
is $\mathbb{F}^{W}$-progressively measurable}\Big{\\}}.$

Under some mild conditions [47], for every choice of parameter $\alpha\in{\cal
A}$ and control $u\in\mathcal{U}[0,T]$, SDE (1) admits a unique solution
$X^{\alpha}$. The performance of the control $u$ is measured by the following
cost functional

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

where $f$ is the running cost rate and $h$ is the terminal cost at $T$. The
goal of the stochastic optimal control problem is to find the “optimal
control” $\bar{u}^{\alpha}$ that minimizes the cost $J$, i.e.

$J(\bar{u}^{\alpha})=\inf_{u\in\mathcal{U}[0,T]}J(u).$ (3)

In the case that the parameter $\alpha$ is given, the classic stochastic
optimal control problem can be solved via the following gradient descent type
optimization.

Assume that the optimal control $\bar{u}^{\alpha}$ is in the interior of
$\mathcal{U}$111In general, this might not be the case. But we only consider
such a case, and leave the general case to our future investigation.. Then we
can deduce by using the Gâteaux derivative of $\bar{u}^{\alpha}$ and the
maximum principle that the gradient process of the cost functional $J$ with
respect to the control process on time interval $t\in[0,T]$ has the following
form [47]

$\nabla
J_{u}(\bar{u}^{\alpha}_{t})={\mathbb{E}}\left[b_{u}(t,\bar{X}^{\alpha}_{t},\bar{u}^{\alpha}_{t},\alpha)^{\top}\bar{Y}^{\alpha}_{t}+\sigma_{u}(t,\bar{X}^{\alpha}_{t},\bar{u}^{\alpha}_{t})^{\top}\bar{Z}^{\alpha}_{t}+f_{u}(t,\bar{X}^{\alpha}_{t},\bar{u}^{\alpha}_{t})^{\top}\right],$
(4)

where subscripts are used to denote partial derivatives of functions. The
stochastic process $(\bar{X}^{\alpha},\bar{Y}^{\alpha},\bar{Z}^{\alpha})$ is
the adapted solution of the following forward backward stochastic differential
equations (FBSDEs) system

$\displaystyle d\bar{X}^{\alpha}_{t}$
$\displaystyle=b(t,\bar{X}^{\alpha}_{t},\bar{u}^{\alpha}_{t},\alpha)dt+\sigma(t,\bar{X}^{\alpha}_{t},\bar{u}^{\alpha}_{t})dW_{t}\qquad$
(FSDE) (5) $\displaystyle d\bar{Y}^{\alpha}_{t}$
$\displaystyle=\big{(}-b_{x}(t,\bar{X}^{\alpha}_{t},\bar{u}^{\alpha}_{t},\alpha)^{\top}\bar{Y}^{\alpha}_{t}-\sigma_{x}(t,\bar{X}^{\alpha}_{t},\bar{u}^{\alpha}_{t})^{\top}\bar{Z}^{\alpha}_{t}-f_{x}(t,\bar{X}^{\alpha}_{t},\bar{u}^{\alpha}_{t})^{\top}\big{)}dt$
$\displaystyle\qquad+\bar{Z}^{\alpha}_{t}dW_{t},\qquad\bar{Y}^{\alpha}_{T}=h_{x}(\bar{X}^{\alpha}_{T})^{\top},$
(BSDE)

where the first equation in (5) is a forward SDE (FSDE, for short), which has
the same form as the state equation (1) in the stochastic optimal control
problem; the second equation is a backward stochastic differential equation
(BSDE) with given terminal condition at time $T$, and we call it the adjoint
equation of the forward state equation [38]. Note that $\bar{Z}^{\alpha}$ is
the integrand in the Itô’ integral of the martingale representation for
$\bar{Y}^{\alpha}$ with respect to $W$. The existence and uniqueness of the
FBSDEs (5) have been well studied under some standard assumptions, and it is
worthy of mentioning that the value of solution $\bar{Y}^{\alpha}$ of the BSDE
depends on state $\bar{X}^{\alpha}$ [37].

Then, with the gradient $\nabla J_{u}$ introduced in Eq. (4) and the FBSDEs
(5) that provides solutions $(Y,Z)$ to describe $\nabla J_{u}$, we introduce
the following gradient descent scheme

$u_{t}^{l+1,\alpha}=u_{t}^{l,\alpha}-\rho_{l}\nabla
J_{u}(u_{t}^{l,\alpha}),\hskip 15.00002pt0\leq t\leq T,\hskip
10.00002ptl=0,1,2,\cdots,L-1,$ (6)

to find the optimal control, where $L$ is the total number of gradient descent
iterations that satisfies certain stopping criteria, $\rho_{l}$ is the step-
size for optimization, and we let $\\{u_{t}^{0,\alpha}\\}_{0\leq t\leq T}$ be
our initial guess for the optimal control.

The above mathematical framework for solving the stochastic optimal control
problem requires the exact value for state $X_{t}^{\alpha}$ in Eq. (1).
However, in reality $X^{\alpha}_{t}$ may not always be explicitly given.
Instead, we can often receive observations on $X^{\alpha}_{t}$ that provide
information of the state. In this work, we let

$M_{t}=X^{\alpha}_{t}+\xi_{t}$ (7)

be noise perturbed direct observational data about $X^{\alpha}_{t}$ with a
Gaussian error term $\xi_{t}\sim N(0,\Sigma)$, where $\Sigma$ is the
covariance matrix, and we assume that the observations are accurate, i.e. the
covariance $\Sigma$ of the error-size is small. Since the observational data
are accurate, when $M_{t}$ is available we can use it as our estimate for
$X_{t}^{\alpha}$ at time $t$. Then, we calculate the gradient process beyond
time $t$ based on the estimate of the current state $X_{t}^{\alpha}$, i.e. the
observation $M_{t}$. This allows us to adjust the control actions based on
observations, which can provide real-time feedback on the controlled state. In
the case that the observations are partial/indirect with relatively large
errors, an extra filtering procedure is needed [5, 8, 10] for the feedback
control, which is out of scope of this work.

On the other hand, in this paper we assume that the parameter vector $\alpha$
in the data driven stochastic optimal control problem (1) - (3) is unknown.
Therefore, in addition to implementing the gradient descent optimization
procedure (6), the primary challenge that we need to address is to design the
optimal control based on the unknown model parameter $\alpha$, and we need an
effective parameter estimation method to determine $\alpha$.

There are extensive studies on estimating model parameters by using
observational data. Most existing methods focus on offline estimation, which
needs to analyze all the historic data in order to carry out parameter
estimation. This may be inefficient when large amount of data are collected,
which will require large-scale data storage and analytics. Moreover, in many
practical scenarios in real-world applications, the true value of the
parameter may change unpredictably. In this case, traditional offline
parameter estimation methods are not suitable since the mixture of static data
may be generated by different parameter values, and it may neglect the
sequential order of data reception, which also contains important information
about the parameter.

The standard technique to dynamically estimate unknown parameters in the
online manner is the optimal filtering, which aims to find the best estimate
for the state of a (hidden) dynamical system by using partial noisy
observational data. The most important approach that uses the optimal
filtering to estimate parameters is to construct an augmented state-parameter
vector, and then estimate the augmented state vector by using the
observational information. Such an approach is often called the “augmented
filter”.

Specifically, the augmented filter considers the following state process

$dS_{t}=B(S_{t})dt+\overline{\sigma}_{t}d\overline{W}_{t},$ (8)

where $S_{t}:=[X^{\alpha}_{t},\alpha_{t}]^{\top}\in\mathbb{R}^{d+q}$ is a
$(d+q)$-dimensional augmented vector that combines the controlled state
$X^{\alpha}_{t}\in\mathbb{R}^{d}$ and the parameter $\alpha\in\mathbb{R}^{q}$.
Note that we treat the unknown parameter $\alpha$ as a stochastic process with
“zero” dynamics in Eq. (8). The function
$B(S_{t}):=[b(t,X^{\alpha}_{t},u^{\alpha}_{t},\alpha),\bm{0}_{q}]^{\top}$
augments the controlled state dynamics $b$ with the “zero” dynamics of the
parameter process $\alpha_{t}$, and $\overline{\sigma}_{t}$,
$\overline{W}_{t}$ represent the augmented diffusion coefficient and the
augmented Brownian motion, respectively. Then, by using the same measurement
process $M_{t}$ introduced in Eq. (7), the augmented filter aims to find the
best estimate for the augmented state $S$ by using the observational
information $\mathcal{M}_{t}:=\sigma(M_{s},0\leq s\leq t)$. The mathematical
expression for the “best estimate” under the optimal filtering framework is
the conditional expectation of $S_{t}$ conditioning on $\mathcal{M}_{t}$, i.e.
${\mathbb{E}}[S_{t}\big{|}\mathcal{M}_{t}]$. In this way, the observational
data $M$ about $X$ give partial noisy observations on $S$, which includes the
unknown parameter. The estimation procedure for state $S$ (based on data about
$X$) will drive the estimation procedure for parameter $\alpha$, which would
eventually generate an estimate for the parameter. Such an augmented optimal
filtering problem can be solved by several effective optimal filtering
methods. The most important approach for solving the optimal filtering problem
is the recursive Bayesian filter, which includes the Kalman type filters [19,
25], the particle filters [2, 27], and our recently developed backward SDE
filter [15, 13, 12, 6, 7, 9, 14].

Although the idea of combining a model and its parameters in an augmented
state process in the augmented filter method provides a feasible framework to
dynamically estimate the unknown parameters, the strategy of using the
accuracy of state estimation to determine wether estimated parameters would
fit the model only provides indirect evidence to the accuracy of parameter
estimation. This would make the augmented filter less accurate especially when
estimating parameters in more complicated dynamical models.

In the following section, we shall introduce an effective novel sample-wise
data driven control method for solving the data driven stochastic optimal
control problem (1) - (3).

##  3 Numerical derivation for the sample-wise data driven control method

Our sample-wise data driven control solver consists two major computational
components: (i) a direct filter method to implement online parameter
estimation, and (ii) an efficient sample-wise optimal control solver, which is
driven by parameter estimation. In what follows, we first introduce the direct
filter algorithm in Section 3.1. Then, in Section 3.2 we introduce a fully
calculated stochastic optimal control solver for the stochastic optimal
control problem with estimated parameters. Finally, in Section 3.3 we
introduce a sample-wise optimization procedure to efficiently calculate the
optimal control, and we let the direct filter based parameter estimation guide
the optimal control solver.

In most situations, it’s very difficult to control the diffusion coefficient
$\sigma$. So we let $\sigma$ be independent from control $u$. Also, in this
work we assume that $\sigma$ does not depend on the state $X$ either (i.e.
$\sigma(t,X^{\alpha}_{t},u_{t})=\sigma(t)$), which is a common assumption in
the parameter estimation problem. Despite of the aforementioned assumptions,
the main difficulties that challenge the online parameter estimation problem
and the stochastic optimal control problem still remain in the computational
framework that we shall introduce.

In order to derive our numerical algorithm, we introduce the following
temporal partition over the time interval $[0,T]$, i.e.

$\Pi^{T}_{N}=\\{t_{n}:0=t_{0}<t_{1}<t_{2}<\cdots<t_{N_{T}}=T\\},$

where $N_{T}$ is the total number of time steps, which also represents how
many sets of observational data $\\{M_{t_{n}}\\}_{n=1}^{N_{T}}$ that will be
received, and we let $\Delta t_{n}:=t_{n+1}-t_{n}$ be the temporal step-size.
For convenience of presentation, in this work we assume that the temporal
partition $\Pi^{T}_{N}$ is uniform and denote $\Delta t=\frac{T}{N}$.

###  3.1 Direct filter for parameter estimation

In the direct filter method, instead of constructing an augmented state-
parameter dynamics Eq. (8) as in the augmented filter approach, we only
consider a parameter process. With the temporal partition $\Pi_{N}^{T}$, we
introduce the following zero dynamics to model the unknown parameter as a
stochastic process:

$\alpha_{t_{n+1}}=\alpha_{t_{n}}+\gamma_{n},\quad n=0,1,2,\cdots,N_{T}-1,$ (9)

where $\gamma_{n}$ is a Gaussian random variable with a pre-determined
covariance, which is an artificial noise added to the parameter process. In
this work, we present our direct filter method following the particle filter
implementation [23]. Other effective nonlinear filtering methods may also be
applied to implement the direct filter following the similar methodology.

At time instant $t_{n}$, assuming we have a set of $M$ samples (“particles”),
denoted by $\\{\zeta_{n}^{(m)}\\}_{m=1}^{M}$, that describe the conditional
probability distribution of the parameter $\alpha_{t_{n}}$ giving the
observational information $\mathcal{M}_{t_{n}}$, we introduce the empirical
distribution $\tilde{P}(\alpha_{t_{n}}\big{|}\mathcal{M}_{t_{n}})$ as follows:

$P(\alpha_{t_{n}}\big{|}\mathcal{M}_{t_{n}})\approx\tilde{P}(\alpha_{t_{n}}\big{|}\mathcal{M}_{t_{n}}):=\frac{1}{M}\sum_{m=1}^{M}\delta_{\zeta_{n}^{(m)}}(\alpha_{t_{n}}),$

where $\delta_{\cdot}$ is a delta function.

From time $t_{n}$ to $t_{n+1}$, the direct filter algorithm is composed of two
steps: a prediction step and an update step:

Prediction step: In the prediction step, we use the Chapman-Kolmogorov
formula, i.e.

$P(\alpha_{t_{n+1}}|\mathcal{M}_{t_{n}})=\int
P(\alpha_{t_{n}}|\mathcal{M}_{t_{n}})P(\alpha_{t_{n+1}}|\alpha_{t_{n}})d\alpha_{t_{n}},$

to generate a set of proposed particles and predict possible parameter values,
where $P(\alpha_{t_{n+1}}|\alpha_{t_{n}})$ is the transition probability that
propagates the particles from time $t_{n}$ to $t_{n+1}$. Since there’s no
prior knowledge about the parameter and we use the zero dynamics for the
parameter process, the transition probability just produces random walks. In
this way, for each particle $\zeta_{n}^{(m)}$ we just add a version of
Gaussian random variable, i.e. $\gamma_{n}^{(m)}$, to inflate
$\zeta_{n}^{(m)}$ with some artificial noise and get

$\tilde{\zeta}_{n+1}^{(m)}=\zeta_{n}^{(m)}+\gamma_{n}^{(m)}.$ (10)

Then, we obtain the empirical distribution that approximates the (prior)
conditional PDF $P(\alpha_{t_{n+1}}\big{|}\mathcal{M}_{t_{n}})$ as follows:

$P(\alpha_{t_{n+1}}\big{|}\mathcal{M}_{t_{n}})\approx\tilde{\pi}(\alpha_{t_{n+1}}\big{|}\mathcal{M}_{t_{n}})=\frac{1}{M}\sum_{m=1}^{M}\delta_{\tilde{\zeta}_{n+1}^{(m)}}(\alpha_{t_{n+1}}).$
(11)

Update step: In the update step, one needs to incorporate the observational
information into the prediction and update the prior distribution to obtain a
posterior distribution. This procedure is typically achieved by Bayesian
inference. In the online parameter estimation problem that we consider, the
observational data only provide observations on the state, which don’t
directly reflect parameter states. On the other hand, since the parameter
variable is not augmented with the state variable in the direct filter method,
we don’t have an estimated state to be compared with the data.

The main idea of the direct filter is to treat the dynamical model as a pseudo
“observation function”. Specifically, we let the state dynamics map a proposed
parameter sample $\tilde{\zeta}_{n+1}^{(m)}$ to a state point
$X_{n+1}^{l,\tilde{\zeta}_{n+1}^{(m)}}$ that describes $X_{t_{n+1}}^{\alpha}$,
i.e.

$X_{n+1}^{l,\tilde{\zeta}_{n+1}^{(m)}}:=X^{\alpha}_{n}+b(t_{n},X^{l,\alpha}_{t_{n}},u_{n}^{l,\alpha},\tilde{\zeta}_{n+1}^{(m)})\Delta
t,$ (12)

where $u_{n}^{l,\alpha}$ is an estimated optimal control at the $l$-th
iteration step and $X^{\alpha}_{n}$ is an estimate for the state variable
$X_{t_{n}}^{\alpha}$. Since data $M_{t_{n+1}}$ directly observe
$X_{t_{n+1}}^{\alpha}$ in this paper, we introduce the following likelihood
function

$p(M_{t_{n+1}}|\tilde{\zeta}_{n+1}^{(m)}):=\exp\left\\{-\frac{\big{(}M_{t_{n+1}}-X_{n+1}^{l,\tilde{\zeta}_{n+1}^{(m)}}\big{)}^{2}}{(\sigma(t_{n}))^{2}\Delta
t+\|\Sigma\|_{L^{2}}^{2}}\right\\}$ (13)

to measure the likelihood of the sample $\tilde{\zeta}_{n+1}^{(m)}$ in fitting
the observational data, where the dynamics Eq. (12) that we use to generate
$X_{n+1}^{l,\tilde{\zeta}_{n+1}^{(m)}}$ contains both the parameter and the
state, and it connects the proposed parameter $\tilde{\zeta}_{n+1}^{(m)}$ to
the estimated state $X_{n+1}^{l,\tilde{\zeta}_{n+1}^{(m)}}$. Note that the
observational error (with covariance $\Sigma$) is incorporated into the
likelihood together with the diffusion term of the model through function
(13).

Then, we apply Bayesian inference with likelihood
$p(M_{t_{n+1}}|\tilde{\zeta}_{n+1}^{(m)})$ and empirical prior
$\tilde{\pi}(\alpha_{t_{n+1}}\big{|}\mathcal{M}_{t_{n}})$ (introduced in Eq.
(11)) to obtain an empirical posterior distribution
$\tilde{\pi}(\alpha_{t_{n+1}}\big{|}\mathcal{M}_{t_{n+1}})$ as follows

$\tilde{\pi}(\alpha_{t_{n+1}}\big{|}\mathcal{M}_{t_{n+1}}):=\frac{p(M_{t_{n+1}}|\\{\tilde{\zeta}_{n+1}^{(m)}\\}_{m})\tilde{\pi}(\alpha_{t_{n+1}}\big{|}\mathcal{M}_{t_{n}})}{C},$

where $C$ is a normalization factor. The practical implementation of the above
Bayesian inference scheme is to assign a weight $w_{n+1}^{(m)}$ to each sample
$\tilde{\zeta}_{n+1}^{(m)}$ with the likelihood value, i.e.
$w_{n+1}^{(m)}:=p(M_{t_{n+1}}|\tilde{\zeta}_{n+1}^{(m)})/C$, and the empirical
distribution $\tilde{\pi}(\alpha_{t_{n+1}}\big{|}\mathcal{M}_{t_{n+1}})$
becomes

$\tilde{\pi}(\alpha_{t_{n+1}}\big{|}\mathcal{M}_{t_{n+1}})=\sum_{m=1}^{M}w_{n+1}^{(m)}\delta_{\tilde{\zeta}_{n+1}^{(m)}}(\alpha_{t_{n+1}}).$
(14)

To stabilize the particle filter algorithm and to alleviate the so-called
“degeneracy issue” [43], we resample the particle cloud
$\\{\tilde{\zeta}_{n+1}^{(m)}\\}_{m}$ and get an equally weighted set of
samples $\\{\zeta_{n+1}^{(m)}\\}_{m}$. There are several well-designed
resampling methods, and we refer to [2, 39, 17] for some examples. As a result
of resampling, we obtain an empirical distribution
$\tilde{P}(\alpha_{t_{n+1}}\big{|}\mathcal{M}_{t_{n+1}})$ described by
$\\{\zeta_{n+1}^{(m)}\\}_{m}$ as

$\tilde{P}(\alpha_{t_{n+1}}\big{|}\mathcal{M}_{t_{n+1}}):=\frac{1}{M}\sum_{m=1}^{M}\delta_{\zeta_{n+1}^{(m)}}(\alpha_{t_{n+1}}),$

which will be used in the next recursive prediction-update procedure. With the
empirical distribution
$\tilde{P}(\alpha_{t_{n+1}}\big{|}\mathcal{M}_{t_{n+1}})$ and the particles
$\\{\zeta_{n+1}^{(m)}\\}_{m}$, we let the direct filter estimated parameter at
the time instant $t_{n+1}$ be $\hat{\alpha}_{t_{n+1}}$, which is defined as

$\hat{\alpha}_{t_{n+1}}:=\tilde{{\mathbb{E}}}[\alpha_{t_{n+1}}|\mathcal{M}_{t_{n+1}}]=\frac{1}{M}\sum_{m=1}^{M}\zeta_{n+1}^{(m)}.$

It’s worthy to mention that the estimated state $X_{n}^{\alpha}$, which
describes the state $X_{t_{n}}^{\alpha}$ in Eq. (12), is needed to implement
the above direct filter method. Since we do not know the exact value of
$X_{t_{n}}^{\alpha}$ and our observational data
$M_{t_{n}}=X_{t_{n}}^{\alpha}+\xi_{t_{n}}$ has high accuracy (with low noise
level for $\xi_{t_{n}}$), in the direct filter method we use the data, i.e.
$M_{t_{n}}$, as our “best” estimate for $X_{t_{n}}^{\alpha}$. In this way, the
estimated state $X_{n}^{\alpha}$ in Eq. (12) is replaced by $M_{t_{n}}$ in
this work.

In the following subsection, we introduce the standard fully calculated
stochastic optimal control solver to implement the gradient descent
optimization scheme (6).

###  3.2 Fully calculated stochastic optimal control solver

We first assume that we have an estimated parameter $\hat{\alpha}$. The
gradient descent optimization procedure will be implemented on discrete
temporal points over the temporal partition $\Pi_{N}^{T}$, i.e.

$u_{t_{n}}^{l+1,\hat{\alpha}}=u_{t_{n}}^{l,\hat{\alpha}}-\rho_{l}\nabla
J_{u}(u_{t_{n}}^{l,\hat{\alpha}}),\hskip 15.00002ptn=0,1,2,\cdots,N_{T}-1,$

and the gradient $\nabla J_{u}(u_{t_{n}}^{l})$ at time step $t_{n}$ is given
by

$\nabla
J_{u}(u^{l,\hat{\alpha}}_{t_{n}})={\mathbb{E}}\left[b_{u}(t_{n},X^{l,\hat{\alpha}}_{t_{n}},u^{l,\hat{\alpha}}_{t_{n}},\hat{\alpha})^{\top}Y^{l,\hat{\alpha}}_{t_{n}}+f_{u}(t_{n},X^{l,\hat{\alpha}}_{t_{n}},u^{l,\hat{\alpha}}_{t_{n}})^{\top}\right],$
(15)

where $Y^{l,\hat{\alpha}}_{t_{n}}$ and $Z^{l,\hat{\alpha}}_{t_{n}}$ are
solutions in the FBSDEs system Eq. (5) driven by the estimated parameter
$\hat{\alpha}$ and the control $u^{l,\hat{\alpha}}$. Hence, numerical
approximations for $X^{l,\hat{\alpha}}_{t}$ and $Y^{l,\hat{\alpha}}_{t}$ are
needed over the time partition $\Pi^{T}_{N}$. In what follows, we shall
introduce a standard numerical scheme to solve the FBSDEs.

Consider the FBSDEs system on the sub-interval $[t_{n},t_{n+1}]$, i.e.

$\displaystyle X^{l,\hat{\alpha}}_{t_{n+1}}$
$\displaystyle=X^{l,\hat{\alpha}}_{t_{n}}+\int_{t_{n}}^{t_{n+1}}b(t,X^{l,\hat{\alpha}}_{t},u^{l,\hat{\alpha}}_{t},\hat{\alpha})dt+\int_{t_{n}}^{t_{n+1}}\sigma(t)dW_{t},\qquad$
(FSDE) (16) $\displaystyle Y^{l,\hat{\alpha}}_{t_{n}}$
$\displaystyle=Y^{l,\hat{\alpha}}_{t_{n+1}}+\int_{t_{n}}^{t_{n+1}}\big{(}b_{x}(t,X^{l,\hat{\alpha}}_{t},u^{l,\hat{\alpha}}_{t},\hat{\alpha})^{\top}Y^{l,\hat{\alpha}}_{t}+f_{u}(t,X^{l,\alpha}_{t},u^{l,\alpha}_{t})^{\top}\big{)}dt$
$\displaystyle\hskip
80.00012pt-\int_{t_{n}}^{t_{n+1}}Z^{l,\hat{\alpha}}_{t}dW_{t}.$ (BSDE)

The forward SDE in Eq. (16) can be solved by the following classic Euler-
Maruyama scheme

$X^{l,\hat{\alpha}}_{t_{n+1}}=X^{l,\hat{\alpha}}_{t_{n}}+b(t_{n},X^{l,\hat{\alpha}}_{t_{n}},u^{l,\hat{\alpha}}_{t_{n}},\hat{\alpha})\Delta
t+\sigma(t_{n})\Delta W_{t_{n}}+R_{n}^{X},$ (17)

where $\Delta W_{t_{n}}:=W_{t_{n+1}}-W_{t_{n}}$ and

$R_{n}^{X}:=\int_{t_{n}}^{t_{n+1}}b(t,X^{l,\hat{\alpha}}_{t},u^{l,\alpha}_{t},\hat{\alpha})dt-b(t_{n},X^{l,\hat{\alpha}}_{t_{n}},u^{l,\hat{\alpha}}_{t_{n}},\hat{\alpha})\Delta
t+\int_{t_{n}}^{t_{n+1}}\sigma(t)dW_{t}-\sigma(t_{n})\Delta W_{t_{n}}$

is the approximation error term.

To solve the BSDE and obtain a numerical solution for $Y$, we use the right-
point formula to approximate the deterministic integral on the right hand side
of the BSDE in Eq. (16) and take the conditional expectation
${\mathbb{E}}_{n}[\cdot]:={\mathbb{E}}[\cdot|X^{l,\alpha}_{t_{n}}]$ on both
side of the equation to get

$Y^{l,\hat{\alpha}}_{t_{n}}={\mathbb{E}}_{n}[Y^{l,\hat{\alpha}}_{t_{n+1}}]+{\mathbb{E}}_{n}\left[b_{x}(t_{n+1},X^{l,\hat{\alpha}}_{t_{n+1}},u^{l,\hat{\alpha}}_{t_{n+1}},\hat{\alpha})^{\top}Y^{l,\hat{\alpha}}_{t_{n+1}}+f_{u}(t_{n+1},X^{l,\hat{\alpha}}_{t_{n+1}},u^{l,\hat{\alpha}}_{t_{n+1}})^{\top}\right]\Delta
t+R_{n}^{Y},$ (18)

where

$\displaystyle R_{n}^{Y}:=$
$\displaystyle{\mathbb{E}}_{n}\left[\int_{t_{n}}^{t_{n+1}}\big{(}b_{x}(t,X^{l,\hat{\alpha}}_{t},u^{l,\alpha}_{t},\hat{\alpha})^{\top}Y^{l,\hat{\alpha}}_{t}+f_{u}(t,X^{l,\hat{\alpha}}_{t},u^{l,\hat{\alpha}}_{t})^{\top}\big{)}dt\right]$
$\displaystyle\quad-{\mathbb{E}}_{n}\left[b_{x}(t_{n+1},X^{l,\hat{\alpha}}_{t_{n+1}},u^{l,\hat{\alpha}}_{t_{n+1}},\hat{\alpha})^{\top}Y^{l,\hat{\alpha}}_{t_{n+1}}+f_{u}(t_{n+1},X^{l,\hat{\alpha}}_{t_{n+1}},u^{l,\hat{\alpha}}_{t_{n+1}})^{\top}\right]\Delta
t$

is the approximation error for the deterministic integral, and the stochastic
integral $\int_{t_{n}}^{t_{n+1}}Z^{l,\hat{\alpha}}_{t}dW_{t}$ is eliminated
due to the martingale property of Itô integral, i.e.
${\mathbb{E}}_{n}\left[\int_{t_{n}}^{t_{n+1}}Z^{l,\hat{\alpha}}_{t}dW_{t}\right]=0$.
Note that we have also used that fact that
$Y^{l,\hat{\alpha}}_{t_{n}}={\mathbb{E}}_{n}\left[Y^{l,\hat{\alpha}}_{t_{n}}\right]$.

By dropping the error terms $R_{n}^{X}$ and $R_{n}^{Y}$ in (17) and (18),
respectively, we have the following numerical scheme for solving the FBSDEs
system (16):

$\displaystyle X^{l,\hat{\alpha}}_{n+1}$
$\displaystyle=X^{l,\hat{\alpha}}_{n}+b(t_{n},X^{l,\hat{\alpha}}_{n},u^{l,\hat{\alpha}}_{t_{n}},\hat{\alpha})\Delta
t+\sigma(t_{n})\Delta W_{t_{n}},$ (19) $\displaystyle Y^{l,\hat{\alpha}}_{n}$
$\displaystyle={\mathbb{E}}_{n}[Y^{l,\hat{\alpha}}_{n+1}]+{\mathbb{E}}_{n}\left[b_{x}(t_{n+1},X^{l,\hat{\alpha}}_{n+1},u^{l,\hat{\alpha}}_{t_{n+1}},\hat{\alpha})^{\top}Y^{l,\hat{\alpha}}_{n+1}+f_{u}(t_{n+1},X^{l,\hat{\alpha}}_{n+1},u^{l,\hat{\alpha}}_{t_{n+1}})^{\top}\right]\Delta
t,$

where $X^{l,\hat{\alpha}}_{n+1}$ and $Y^{l,\hat{\alpha}}_{n+1}$ are
approximations to solutions $X$ and $Y$ in Eq. (16) with a given estimated
parameter $\hat{\alpha}$ and an estimated optimal control $u^{l,\alpha}$.
Although $Z^{l,\hat{\alpha}}$ is also a solution in Eq. (16), since it does
not appear in the expression (15) for the gradient $\nabla
J_{u}(u^{l,\hat{\alpha}}_{t_{n}})$ and $Z$ is not included in the numerical
schemes for $X$ and $Y$, we don’t need to obtain a numerical solution for $Z$
in the data driven control problem that we try to solve in this work. Then,
with approximation schemes (19), the gradient $\nabla
J_{u}(u^{l,\hat{\alpha}}_{t_{n}})$ described in Eq. (15) can be approximated
as

$\nabla
J_{u}(u^{l,\hat{\alpha}}_{t_{n}})\approx{\mathbb{E}}\left[b_{u}(t_{n},X^{l,\hat{\alpha}}_{n},u^{l,\hat{\alpha}}_{t_{n}},\hat{\alpha})^{\top}Y^{l,\hat{\alpha}}_{n}+f_{u}(t_{n},X^{l,\hat{\alpha}}_{n},u^{l,\hat{\alpha}}_{t_{n}})^{\top}\right].$
(20)

Since we apply online methods to estimate the unknown parameters in this
paper, the estimated parameter $\hat{\alpha}$ in Eq. (20) becomes
$\hat{\alpha}_{t_{n}}$, which is the dynamically estimated parameter generated
at time $t_{n}$. Note that other online parameter estimation methods can also
be applied in the above stochastic optimal control solver to drive the design
of the optimal control. As a result, with the estimated parameter
$\hat{\alpha}_{t_{n}}$ at time $t_{n}$, we estimate the gradient $\nabla
J_{u}$ for a future time $t_{n_{s}}$ ($n_{s}\geq n$) by

$\nabla
J_{u}(u^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}})\approx{\mathbb{E}}\left[b_{u}(t_{n_{s}},X^{l,\hat{\alpha}_{t_{n}}}_{n_{s}},u^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}},\hat{\alpha}_{t_{n}})^{\top}Y^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}+f_{u}(t_{n_{s}},X^{l,\hat{\alpha}}_{n_{s}},u^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}})^{\top}\right],\quad
n\leq n_{s}\leq N_{T}.$ (21)

When $n=0$, the estimated parameter $\hat{\alpha}_{t_{0}}$ is chosen as our
initial guess for the parameter.

In order to implement the gradient descent optimization procedure, we need to
approximate the conditional expectation ${\mathbb{E}}_{n}[\cdot]$ in Eq. (19)
and the expectation ${\mathbb{E}}[\cdot]$ in Eq. (21) with respect to the
state process $X$. The standard approach to approximate expectations is the
Monte Carlo method.

To be specific, we generate $P$ realizations of simulated paths for
$X^{l,\hat{\alpha}_{t_{n}}}_{n_{s}+1}$ as

$X^{l,\hat{\alpha}_{t_{n}},(p)}_{n_{s}+1}=X^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}+b(t_{n_{s}},X^{l,\hat{\alpha}_{t_{n}}}_{n_{s}},u^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}},\hat{\alpha}_{t_{n}})\Delta
t+\sigma(t_{n_{s}})\Delta W^{(p)}_{t_{n_{s}}},\quad n_{s}\geq n,\
p=1,2,\cdots,P,$ (22)

where $\Delta W^{(p)}_{t_{n}}\sim N(0,\Delta t\cdot I_{d})$ is the $p$-th
realization of $\Delta W_{t_{n}}$, and $X^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}$ is
a “point of interest” in the state space at which we approximate solution
$Y_{n_{s}}^{l,\hat{\alpha}_{t_{n}}}$. Then, the Monte Carlo approximation
based numerical solution $\bar{Y}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}$
corresponding to the state point $X^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}$ becomes

$\displaystyle\hat{Y}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}(X^{l,\hat{\alpha}_{t_{n}}}_{n_{s}})=$
$\displaystyle\sum_{p=1}^{P}\frac{\hat{Y}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}+1}(X^{l,\hat{\alpha}_{t_{n}},(p)}_{n_{s}+1})}{P}+\sum_{p=1}^{P}\frac{b_{x}(t_{n_{s}+1},X^{l,\hat{\alpha}_{t_{n}},(p)}_{n_{s}+1},u^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}+1}})^{\top}\hat{Y}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}+1}(X^{l,\hat{\alpha}_{t_{n}},(p)}_{n_{s}+1})}{P}\Delta
t$ (23)
$\displaystyle\quad+\sum_{p=1}^{P}\frac{f_{x}(t_{n_{s}+1},X^{l,\hat{\alpha}_{t_{n}},(p)}_{n_{s}+1},u^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}+1}})^{\top}}{P}\Delta
t.$

In order to approximate the expectation ${\mathbb{E}}[\cdot]$ in Eq. (21), we
need to simulate another set of sample paths for $X$. To this end, we simulate
altogether $Q$ trajectories of $X$ as

$X^{l,\hat{\alpha}_{t_{n}}}_{n_{s}+1,(q)}=X^{l,\hat{\alpha}_{t_{n}}}_{n_{s},(q)}+b(t_{n_{s}},X^{l,\hat{\alpha}_{t_{n}}}_{n_{s},(q)},u^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}},\hat{\alpha}_{t_{n}})\Delta
t+\sigma(t_{n_{s}})\Delta W_{t_{n_{s}},(q)},\ \ n_{s}\geq n,\ q=1,2,\cdots,Q.$
(24)

Then, the fully calculated Monte Carlo approximation for the gradient process,
denoted by $\nabla\hat{J}_{u}(u^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}})$, is
computed as follows:

$\nabla\hat{J}_{u}(u^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}})=\frac{1}{Q}\sum_{q=1}^{Q}\left(b_{u}(t_{n_{s}},X^{l,\hat{\alpha}_{t_{n}}}_{n_{s},(q)},u^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}},\hat{\alpha}_{t_{n}})^{\top}\hat{Y}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}(X^{l,\hat{\alpha}_{t_{n}}}_{n_{s},(q)})+f_{u}(t_{n_{s}},X^{l,\hat{\alpha}_{t_{n}}}_{n_{s},(q)},u^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}})^{\top}\right),$
(25)

and the gradient descent procedure for the optimal control becomes

$\hat{u}_{n_{s}}^{l+1,\hat{\alpha}_{t_{n}}}=\hat{u}_{n_{s}}^{l,\hat{\alpha}_{t_{n}}}-\rho_{l}\nabla\hat{J}_{u}(\hat{u}_{n_{s}}^{l,\hat{\alpha}_{t_{n}}}),\hskip
15.00002ptn\leq n_{s}\leq N_{T}.$ (26)

We refer to [22] for more detailed discussions of the above fully calculated
stochastic optimal control solver.

There are two major computational challenges to carry out the fully calculated
gradient descent optimization (26). First of all, there are two Monte Carlo
simulations involved in the approximation of
$\nabla\hat{J}_{u}(u^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}})$: (i) Monte Carlo
simulation for the conditional expectation ${\mathbb{E}}_{n}[\cdot]$ as
described in Eq. (23) (to be implemented on each spatial-state point
$X^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}$ at each time step $t_{n_{s}}$), and (ii)
Monte Carlo simulation for the expectation ${\mathbb{E}}[\cdot]$ as described
in Eq. (25). Hence it requires a very large number of Monte Carlo samples to
approximate the gradient at each iteration step. Secondly, in order to use
$\hat{Y}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}$ to calculate the gradient
$\nabla\hat{J}_{u}(\hat{u}_{n_{s}}^{l,\hat{\alpha}_{t_{n}}})$, a complete
numerical approximation for $\hat{Y}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}$ in the
entire state space is needed since
$\\{X^{l,\hat{\alpha}_{t_{n}}}_{n_{s}+1,(q)}\\}_{q=1}^{Q}$ are randomly
generated sample points, which may appear anywhere in the state space. Those
two computational challenges make the fully calculated stochastic optimal
control solver (22) - (26) extremely computational expensive – especially when
solving high dimensional problems.

In what follows, we shall introduce a sample-wise data driven control solver
to efficiently implement the gradient descent optimization scheme (26) and to
incorporate the probability distribution
$\tilde{P}(\alpha_{t_{n}}|\mathcal{M}_{t_{n}})$ for the unknown parameter
(obtained by the direct filter) into the optimization procedure.

###  3.3 Sample-wise data driven optimal control solver

The novel concept that we introduce in the sample-wise data driven control
solver is to treat the random samples generated by the state variable as “data
samples”. Then, we apply stochastic approximation in the gradient descent
optimization procedure to improve the efficiency of numerical implementation
of the stochastic optimal control solver.

Recall that in the fully calculated stochastic optimal control solver, we use
a set of $P$ samples, i.e.
$\\{X^{l,\hat{\alpha}_{t_{n}},(p)}_{n_{s}+1}\\}_{p=1}^{P}$, to represent the
state variable $X_{t_{n_{s}+1}}^{l,\hat{\alpha}_{t_{n}}}$, and we approximate
the conditional expectation ${\mathbb{E}}_{n}[\cdot]$ in scheme (19) by using
the Monte Carlo average of all $P$ samples. In addition, we simulate $Q$
realizations of sample paths, i.e.
$\\{X^{l,\hat{\alpha}_{t_{n}}}_{n_{s},(q)}\\}_{q=1,n_{s}\geq n}^{Q}$, to
approximate the expectation ${\mathbb{E}}[\cdot]$ in the gradient process.
Note that both sample sets
$\\{X^{l,\hat{\alpha}_{t_{n}},(p)}_{n_{s}+1}\\}_{p=1}^{P}$ and
$\\{X^{l,\hat{\alpha}_{t_{n}}}_{n_{s},(q)}\\}_{q=1,n_{s}\geq n}^{Q}$ are
generated to approximate expectations, which appear in the gradient process.
Therefore, we can apply the stochastic approximation method to the gradient
descent optimization procedure and approximate the expectations by using a
single realization of state sample (or a small batch of samples) to represent
the state process.

Specifically, we generate a sample path for the state process, denoted by
$\\{\tilde{X}_{n_{s}}^{l,\hat{\alpha}_{t_{n}}}\\}_{n_{s}=n}^{N_{T}}$,
corresponding to the $l$-th gradient descent iteration step as follows

$\tilde{X}_{n_{s}+1}^{l,\hat{\alpha}_{t_{n}}}=\tilde{X}_{n_{s}}^{l,\hat{\alpha}_{t_{n}}}+b(t_{n_{s}},\tilde{X}_{n_{s}}^{l,\hat{\alpha}_{t_{n}}},\tilde{u}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}},\hat{\alpha}_{t_{n}})\Delta
t+\sigma(t_{n_{s}})\Delta W^{(l)}_{t_{n_{s}}},\qquad n\leq n_{s}\leq N_{T}-1,$
(27)

where $\tilde{u}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}$ is an estimated optimal
control and $W^{(l)}_{t_{n_{s}}}$ is a sample of Brownian motion corresponding
to the iteration index $l$. With the simulated sample path representing the
state process, we approximate the conditional expectations in the scheme (19)
by using only one sample of $X$, i.e.
$\\{\tilde{X}_{n_{s}}^{l,\hat{\alpha}_{t_{n}}}\\}_{n_{s}=n}^{N_{T}}$, instead
of using altogether $P$ samples in the fully calculated scheme (23). As a
result, we have the following sample-wise approximation for solution $Y$ of
the BSDE in Eq. (16):

$\displaystyle\tilde{Y}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}(\tilde{X}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}})=$
$\displaystyle\tilde{Y}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}+1}(\tilde{X}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}+1})+\Big{(}b_{x}(t_{n_{s}+1},\tilde{X}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}+1},\tilde{u}^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}+1}})^{\top}\tilde{Y}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}+1}(\tilde{X}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}+1})$
(28) $\displaystyle\hskip
30.00005pt+f_{x}(t_{n_{s}+1},\tilde{X}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}+1},\tilde{u}^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}+1}})^{\top}\Big{)}\Delta
t,\qquad n\leq n_{s}\leq N_{T}-1.$

Moreover, we also use the sample path
$\\{\tilde{X}_{n_{s}}^{l,\hat{\alpha}_{t_{n}}}\\}_{n_{s}=n}^{N_{T}}$ generated
in Eq. (27) to represent the simulated sample paths
$\\{X^{l,\hat{\alpha}_{t_{n}}}_{n_{s},(q)}\\}_{q=1,n_{s}\geq n}^{Q}$ that we
use in the Monte Carlo approximation for the gradient process as decribed in
Eq. (25). In this way, we introduce the following sample-wise approximation
for the gradient process

$\displaystyle\nabla\tilde{J}_{u}(\tilde{u}^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}})=$
$\displaystyle
b_{u}(t_{n_{s}},\tilde{X}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}},\tilde{u}^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}},\hat{\alpha}_{t_{n}})^{\top}\tilde{Y}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}(\tilde{X}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}})$
(29)
$\displaystyle\quad+f_{u}(t_{n_{s}},\tilde{X}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}},\tilde{u}^{l,\hat{\alpha}_{t_{n}}}_{t_{n_{s}}})^{\top},\hskip
40.00006pt\qquad n\leq n_{s}\leq N_{T},$

where $\tilde{X}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}$ and
$\tilde{Y}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}}(\tilde{X}^{l,\hat{\alpha}_{t_{n}}}_{n_{s}})$
are calculated via Eq. (27) and Eq. (28), respectively. Then, we implement the
following stochastic gradient descent optimization procedure to find the
optimal control

$\tilde{u}_{n_{s}}^{l+1,\hat{\alpha}_{t_{n}}}=\tilde{u}_{n_{s}}^{l,\hat{\alpha}_{t_{n}}}-\rho_{l}\nabla\tilde{J}_{u}(\tilde{u}_{n_{s}}^{l,\hat{\alpha}_{t_{n}}}),\qquad\qquad
n\leq n_{s}\leq N_{T},\quad l=0,1,2,\cdots,L.$ (30)

The numerical schemes (27) - (30) constitute the sample-wise data driven
control method, and the dynamically estimated parameter $\hat{\alpha}_{t_{n}}$
is the driver that guides the sample-wise optimal control solver.

It’s worthy to point out that the scheme (28) only provides a sample-wise
approximation to $Y$. It does not calculate the numerical solution for $Y$
completely, and accuracy of the sample-wise approximation $\tilde{Y}$ is not
guaranteed. Actually, approximation for $Y$ in the entire state space of $X$
is needed if we want to obtain accurate numerical solutions for BSDEs, and
Monte Carlo simulation for conditional expectations can not be avoided (unless
other numerical integration methods are applied). However, since the BSDE that
we consider in the stochastic optimal control problem is the adjoint equation
of the (forward) state equation and $Y$ only appears in the gradient process
under expectation, obtaining a complete numerical solution for $Y$ is not
needed in the optimization procedure for finding the optimal control. In this
connection, the justification for applying stochastic approximation in
stochastic gradient descent type optimization can be used to explain the
methodology of our sample-wise optimal control solver.

In this work, we use the direct filter introduced in Section 3.1 to be the
online parameter estimation method to calculate the estimated parameter
$\hat{\alpha}_{t_{n}}$. In order to sufficiently incorporate parameter
distribution $\tilde{P}(\alpha_{t_{n}}|\mathcal{M}_{t_{n}})$, which contains
more data information then the average value of parameter particles, into the
stochastic optimal control solver, we randomly select
$\zeta_{n}^{(l)}\sim\\{\zeta_{n}^{(m)}\\}_{m=1}^{M}$ to replace the mean value
estimate $\hat{\alpha}_{t_{n}}$ in each stochastic gradient descent iteration
step. In this way, we rewrite the approximation scheme (29) for the gradient
process as

$\displaystyle\nabla\tilde{J}_{u}(\tilde{u}^{l,\zeta_{n}^{(l)}}_{t_{n_{s}}})=$
$\displaystyle
b_{u}(t_{n_{s}},\tilde{X}^{l,\zeta_{n}^{(l)}}_{n_{s}},\tilde{u}^{l,\zeta_{n}^{(l)}}_{t_{n_{s}}},\hat{\alpha}_{t_{n}})^{\top}\tilde{Y}^{l,\zeta_{n}^{(l)}}_{n_{s}}(\tilde{X}^{l,\zeta_{n}^{(l)}}_{n_{s}})$
(31)
$\displaystyle\quad+f_{u}(t_{n_{s}},\tilde{X}^{l,\zeta_{n}^{(l)}}_{n_{s}},\tilde{u}^{l,\zeta_{n}^{(l)}}_{t_{n_{s}}})^{\top},\hskip
40.00006ptn\leq n_{s}\leq N_{T}.$

As a result, when we run a large enough number of iteration steps, the
distribution of parameter particles will be well incorporated into the
stochastic optimal control solver.

With the sample-wise gradient approximator $\nabla\tilde{J}_{u}$ introduced in
Eq. (31), we carry out the following iteration scheme to implement the
optimization procedure for determining the optimal control:

$\tilde{u}_{n_{s}}^{l+1,\zeta_{n}^{(l)}}=\tilde{u}_{n_{s}}^{l,\zeta_{n}^{(l)}}-\rho_{l}\nabla\tilde{J}_{u}(\tilde{u}_{n_{s}}^{l,\zeta_{n}^{(l)}}),\qquad
n\leq n_{s}\leq N_{T},\quad l=0,1,2,\cdots,L.$ (32)

###  3.4 Summary of the algorithm

In Table 1, we provide a pseudo algorithm to summarize the computational
framework of our online method for data driven stochastic optimal control with
unknown model parameters.

Table 1: Summary of the algorithm Algorithm: Sample-wise data driven control
(SW-DDC)
---
0: the number of particles $M$, the learning rate $\rho_{l}$ and the number of
iterations $L\in\mathbb{N}$. Generate parameter particles
$\\{\zeta_{0}^{(m)}\\}_{m=1}^{M}$ to represent the initial guess of the
unknown parameter and introduce an initial guess for the optimal control.  0:
$n=0,1,2,\cdots,N_{T}-1$, do $\bullet$ Implement the sample-wise optimal
control solver with $\\{\zeta_{n}^{(m)}\\}_{m=1}^{M}$: for $l=1,2,\cdots L$ -
Generate a sample path
$\\{\tilde{X}_{n_{s}}^{l,\zeta_{n}^{(l)}}\\}_{n_{s}=n}^{N_{T}}$ by using the
scheme (27), where $\zeta_{n}^{(l)}\in\\{\zeta_{n}^{(m)}\\}_{m=1}^{M}$ is a
representation for $\hat{\alpha}_{t_{n}}$ corresponding to the $l$-th
iteration. - Generate a sample-wise approximation
$\\{\tilde{Y}^{l,\zeta_{n}^{(l)}}_{n_{s}}(\tilde{X}^{l,\zeta_{n}^{(l)}}_{n_{s}})\\}_{n_{s}=n}^{N_{T}}$
corresponding to
$\\{\tilde{X}_{n_{s}}^{l,\zeta_{n}^{(l)}}\\}_{n_{s}=n}^{N_{T}}$ by using the
scheme (28) for solution $Y$ of the adjoint BSDE. - Generate a sample-wise
gradient process
$\\{\nabla\tilde{J}_{u}(\tilde{u}^{l,\zeta_{n}^{(l)}}_{t_{n_{s}}})\\}_{n_{s}=n}^{N_{T}}$
according to Eq. (31). - Carry out gradient descent iteration Eq. (32) to
obtain an updated optimal control
$\\{\tilde{u}_{n_{s}}^{l+1,\zeta_{n}^{(l)}}\\}_{n_{s}=n}^{N_{T}}$. end
Implement the direct filter method to generate
$\\{\zeta_{n+1}^{(m)}\\}_{m=1}^{M}$: - Prediction step: Generate a set of
predicted parameter particles $\\{\tilde{\zeta}_{n+1}^{(m)}\\}_{m=1}^{M}$ via
Eq. (10). - Update step: Assign weight $w_{n+1}^{(m)}$ to each particle
according to the likelihood function defined by Eq. (13) to generate a
posterior distribution that incorporates the new observational data
$M_{t_{n+1}}$. - Resampling: Resample parameter particles to generate an
equally weighted set of particles $\\{\zeta_{n+1}^{(m)}\\}_{m=1}^{M}$ to
formulate the empirical distribution of the estimated parameter
$\hat{\alpha}_{t_{n+1}}$. 0:

As a summary, the sample-wise data driven control (SW-DDC) utilizes the direct
filter as an online parameter estimation method to drive the design of the
optimal control. Instead of implementing the fully calculated stochastic
optimal control solver (introduced in Section 3.2), which requires
approximating the gradient process with respect to the optimal control in the
entire state space, we treat the random samples that describe the state
variable (as well as the parameter variable) as “data samples” and adopt the
methodology of stochastic approximation to approximate expectations in the
gradient process by using a single simulated state (parameter) sample or a
small batch of samples. As a result, the SW-DDC solver can efficiently solve
the data driven stochastic optimal control problem with unknown parameters.

##  4 Numerical Experiments

In this section, we use two numerical examples to examine the performance of
our SW-DDC method. In the first example, we solve a classic linear-quadratic
(LQ) stochastic optimal control problem with unknown parameters in the state
process. Since the LQ problem can be solved analytically, we use this example
as a benchmark scenario to verify the efficiency and effectiveness of our
method. In the second example, we solve an unmanned aerial vehicle (drone)
maneuvering problem. In order to make the control problem data driven, we
assume that there’s an unknown parameter that can influence the control of the
drone, and we use this example to demonstrate the applicability of our SW-DDC
method in solving practical problems.

###  4.1 Example 1. Linear-quadratic data driven stochastic optimal control
problem

Consider the following linear controlled dynamical system

$dX_{t}=A(t)X_{t}dt+BU_{t}dt+CdW_{t},\quad X_{0}=x_{0},$ (33)

where $A(t)$, $B$ and $C$ are given coefficients and $A(t)$ contains unknown
parameters, which will influence the behavior of the dynamics, $X_{t}$ is the
multi-dimensional state process with the given initial position $X_{0}$, and
$U_{t}$ is the multi-dimensional control process. The cost function $J$ is
defined by

$J(U^{M})=E\Big{[}\frac{1}{2}\int_{0}^{T}\Big{(}\langle
QX_{t},X_{t}\rangle+\langle RU_{t},U_{t}\rangle\Big{)}dt+\frac{1}{2}\langle
FX_{T},X_{T}\rangle\Big{]}$ (34)

The data that we use to estimate the unknown parameters is noise perturbed
direct observations on $X$, i.e.

$M_{t}=X_{t}+\xi_{t},$ (35)

where $\xi_{t}\sim N(0,\Sigma)$ is a Gaussian noise with covariance $\Sigma$.

The optimal control $\bar{U}_{t}$ for the above LQ stochastic optimal control
problem can be derived analytically as

$\bar{U}_{t}=-R^{-1}B^{T}P(t)X_{t},$ (36)

where $P(t)$ is the unique solution of the following Riccati equation

$\frac{dP(t)}{dt}=-P(t)A(t)-A^{T}(t)P(t)+P(t)BR^{-1}B^{T}P(t)-Q,\quad P(T)=F.$

#### Case 1: 1D state dynamics with 1 unknown parameter

We first carry out an experiment to demonstrate the necessity of using online
parameter estimation methods. To this end, we solve the data driven optimal
control problem over the time interval $t\in[0,1]$ with $A(t)=\alpha$, where
$\alpha$ is the unknown parameter to be estimated, and we let the true value
for parameter $\alpha$ be $1$. Although we don’t assume that $\alpha$ is a
time dependent parameter, in practice the true parameter may change
unexpectedly in some circumstances. In this experiment, we let $\alpha$ change
from $1$ to $5$ at time instant $t=0.5$. Since the switch of the true
parameter at the time instant $t=0.5$ is unexpected, the conventional offline
parameter estimation methods cannot adjust the change of the true parameter
due to the mixture data from two different true parameters. The other
coefficients in the data-driven LQ control problem (33)-(35) are chosen as
$B=0.5$, $C=0.01$, $Q=1$, $R=0.1$, $F=1$, and $\Sigma=0.001$.

In our numerical experiment, we let $\Delta t=0.02$, i.e. we receive $50$ sets
of observational data, and we use altogether $M=200$ particles in the SW-DDC
algorithm.

(a) Parameter estimation (b) Estimation for the optimal control

Figure 1: Comparison between fully calculated control solver with sample-wise
control solver

In Figure 1, we present the estimated parameter and the estimated optimal
control by using the SW-DDC algorithm. The black straight lines in Figure 1
(a) give the real parameter $\alpha$ over the time interval $[0,1]$. The black
curve marked by crosses in Figure 1 (b) gives the optimal control computed via
Eq. (36) by using the real parameter $\alpha$. Note that the first piece of
the true optimal control in subplot Figure 1 (b) is obtained based on the
assumption that the true parameter is $\alpha=1$ over the entire time interval
$[0,1]$ since the switch of the parameter occurs unexpected at the time
instant $t=0.5$. The red curves marked by circles in subplots Figure 1 (a) and
(b) are the estimated parameter and the estimated optimal control,
respectively. From this figure, we can see that the direct filter quickly
captured the real model parameter, and it can detect and modify its estimation
to adjust the change of true parameter at the time instant $t=0.5$. On the
other hand, the estimated optimal control presented in subplot Figure 1 (b)
shows that our SW-DDC algorithm can provide accurate estimates for the optimal
control, and it only takes a few steps to re-capture the analytical optimal
control driven by the switched real parameter.

In the second experiment for the 1D case, we demonstrate the efficiency of our
sample-wise optimal control solver by comparing with the conventional fully-
calculated stochastic optimal control solver described by (22) - (26) (see
[22] for details). The spatial points that we use to represent the state
variable $X_{n_{s}}^{l,\hat{\alpha}_{t_{n}}}$ in the fully-calculated
stochastic optimal control solver is the uniform mesh in the state space.
Since we aim to compare the efficiency between two optimal control methods, we
use the direct filter as our parameter estimation method to drive both optimal
control solvers.

In this experiment, we let $A(t)=2\alpha\sin(t)$, and all the other
coefficients remain the same as in the first experiment. The true parameter
$\alpha$ that we choose is $\alpha=1$.

(a) Fully calculated control solver with coarse mesh (b) Fully calculated
control solver with medium mesh (c) Fully calculated control solver with fine
mesh (d) Sample-wise control solver

Figure 2: Comparison between fully calculated control solver with sample-wise
optimal control solver

In Figure 2, we show the comparison results in estimating the optimal control.
In subplots Figure 2 (a), (b), (c), we present the estimated optimal control
obtained by the fully-calculated stochastic optimal control solver with coarse
mesh ($\Delta x=0.8$, $\Delta t=0.1$), medium mesh ($\Delta x=0.6$, $\Delta
t=0.075$) and fine mesh ($\Delta x=0.4$, $\Delta t=0.05$). The black curve
marked by crosses in each subplot is the true optimal control (calculated
based on its corresponding temporal discretization), and the green curve
marked by diamonds is the corresponding estimated optimal control. In subplot
Figure 2 (d), we present the estimated optimal control obtained by using the
sample-wise optimal control solver in the SW-DDC framework with temporal step-
size $\Delta t=0.05$, where the red curve marked by circles is our estimated
control, and the estimated optimal control is obtained through $10^{4}$
iteration steps in the optimization procedure. We can see from this figure
that by using finer meshes, the fully-calculated stochastic optimal control
solver can provide better estimated optimal control, while our sample-wise
optimal control solver outperforms the fully-calculated stochastic optimal
control solver with the finest mesh.

In Table 2, we present the computational cost for each method, and we carry
out our numerical experiments on an Apple M1 CPU with 8 cores and 16 GB
memory. We can see from Table 2 and the comparison results in Figure 2 that
the SW-DDC has much lower computational cost while is still more accurate
compared with the the fully-calculated stochastic optimal control solver with
the finest mesh.

Table 2: Example 1. 1D computational cost comparison Method | Coarse mesh | Medium mesh | Fine mesh | SW-DDC
---|---|---|---|---
CPU time (seconds) | $40.6$ | $84.2$ | $287.2$ | $4.49$

Through the numerical experiments that we carried out in solving the 1D LQ
data driven stochastic optimal control problems, we have verified that the
direct filter in the SW-DDC algorithm works well for unexpected changes in the
true parameter, which indicates the necessity of online parameter estimation
in data driven control. Moreover, the comparison experiments in this 1D case
show that the sample-wise optimal control solver in our SW-DDC framework
outperforms conventional fully-calculated stochastic optimal control solver in
efficiency. Due to the fact that our sample-wise optimal control solver does
not solve the adjoint BSDE in the entire state space, the efficiency feature
of our method would be more advantageous in solving higher dimensional
problems. In this way, we have demonstrated the necessity of using sample-wise
optimal control solver to solve data driven control problems.

In what follows, we shall demonstrate the effectiveness and efficiency of the
direct filter, and we are going to compare the direct filter with other state-
of-the-art online parameter estimation methods as a “driver” to drive the
sample-wise optimal control solver in the SW-DDC framework.

#### Case 2: 2D state dynamics with 1 unknown parameter

We consider a 2D case of the the LQ data-driven stochastic optimal control
problem (33)-(35), and we let $A(t)=\big{(}\alpha\sin(t),\cos(t)\big{)}\
I_{2}$, $B=(0.5,0.5)\ I_{2}$, and $C=0.1\ I_{2}$ in (33), where $\alpha$ in
$A(t)$ is a scaler unknown parameter. The control term $U_{t}$ is chosen as a
scaler in this 2D problem as well. For the cost functional, we let $Q=I_{2}$,
$R=1$ and $F=I_{2}$. In the numerical experiments, we solve this data driven
stochastic optimal control problem over the time interval $[0,1]$ by choosing
the temporal step-size $\Delta t=0.02$, i.e. $N_{T}=50$, and we let the
covariance matrix $\Sigma$ of the observational noise be
$\sigma=(0.01)^{2}I_{2}$. The initial state of the controlled process is
chosen as $X_{0}=(2,-2)^{\top}$.

(a) Parameter estimation by the direct filter (b) Parameter estimation by the
AugPF

Figure 3: Example 1. Comparison between the the direct filter and the AugPF in
parameter estimation (2D case)

The goal of the numerical experiments for this 2D problem is to show that the
direct filter method in the SW-DDC framework will outperform the augmented
particle filter (AugPF), which is a particle filter implementation of the
augmented filter method, and both the direct filter and the AugPF will drive
the same sample-wise optimal control solver (due to its high efficiency
advantage) to determine the data driven optimal control. In Figure 3, we
present the comparison between the direct filter (subplot (a)) and the AugPF
(subplot (b)) in estimating the unknown parameter $\alpha$, where we have used
$100$ particles in the direct filter and $1000$ particle in the AugPF. The
black solid line in each subplot indicates the true parameter value. We can
see from this figure that the direct filter captured the real parameter near
the time instant $t=0.4$, and the AugPF cannot find the target parameter in
the entire time interval $[0,1]$.

(a) Optimal control estimation based on direct filter parameter estimation (b)
Optimal control estimation based on AugPF parameter estimation

Figure 4: Example 1. Comparison of estimated optimal controls corresponding to
each parameter estimation method (2D case)

In Figure 4, we present the estimated optimal control obtained by the sample-
wise optimal control solver, which is driven by estimated parameters using the
direct filter (subplot (a)) and the AugPF (subplot (b)). We can see that with
accurate estimation for the parameter, the direct filter driven SW-DDC can
provide accurate estimation for the optimal control. On the other hand, since
the AugPF failed to find the true parameter, the sample-wise optimal control
solver cannot provide accurate estimates for the optimal control.

We also want to mention that the failure of the AugPF in the above experiment
could be caused by inefficiency of the particle filter in estimating the
augmented state process. For a relatively low dimensional problem, the
parameter estimation performance of the AugPF may improve by increasing the
particle-size.

(a) Parameter estimation by AugPF (b) Optimal control estimation based on
AugPF estimated parameter

Figure 5: Example1. Performance of the AugPF with $20,000$ particles (2D case)

In Figure 5, we present the parameter estimation performance (subplot (a) )
and the optimal control estimation performance (subplot (b) ) of the AugPF
driven SW-DDC by using $20000$ particles in the AugPF procedure. We can see
that in this experiment the AugPF driven SW-DDC could provide much better
estimation for both the parameter and the optimal control when significantly
more particles are used.

(a) Comparison of RMSEs in parameter estimation (b) Comparison of RMSEs in
optimal control estimation

Figure 6: Example1. Comparison of RMSEs (2D case)

To further compare the performance of the direct filter and the AugPF in
driving the sample-wise optimal control solver under our SW-DDC framework, we
carry out the above comparison experiments repeatedly $20$ times (with
different particle-sizes in the AugPF) and present the root mean square errors
(RMSEs) of each implementation in Figure 6, where Figure 6 (a) shows the RMSEs
of parameter estimation and Figure 6 (b) shows the RMSEs of optimal control
estimation. We can see from this figure that by using more and more particles,
the AugPF can provide more and more accurate parameter estimation results,
which will result better and better estimated optimal control. On the other
hand, the direct filter can provide the best estimated parameter, which
resulted the best estimated optimal control. In Table 3, we present the
comparison of (average) computational cost between the AugPF (with $1,000$
particles, $20,000$ particles and $40,000$ particles) driven SW-DDC and the
direct filter (with only $100$ particles) driven SW-DDC.

Table 3: Example 1. 2D computational cost comparison Method | AugPF (1K) | AugPF (20K) | AugPF (40K) | Direct filter (100)
---|---|---|---|---
CPU time (seconds) | $73.0$ | $126.9$ | $161.8$ | $67.9$

We can see from this table that the computational cost for the direct filter
is the lowest since it only used $100$ particles to estimate the unknown
parameter.

#### Case 3: 4D state dynamics with 2 unknown parameters

The third case that we consider in the LQ data driven stochastic control
example has 4D state dynamics with 2 unknown parameters. Specifically, we let
$A(t)=\big{(}\sin(t),\cos(t),\alpha,\beta\big{)}\ I_{4}$,
$B=\big{(}0.5,0;0.5,0;1,0;1,0\big{)}^{\top}$, $C=0.1\ I_{4}$, where $\alpha$
and $\beta$ in $A(t)$ are 2 unknown parameters to be estimated, and $U_{t}$ is
a 2D control vector. For the cost functional, we let $Q=I_{4}$, $R=I_{2}$ and
$F=I_{4}$. In the numerical experiments, we solve this data driven control
problem over the time interval $[0,1]$ by choosing the temporal step-size
$\Delta t=0.025$, i.e. $N_{T}=40$, and we let the covariance matrix $\Sigma$
of the observational noise be $\sigma=(0.01)^{2}I_{2}$. The initial state of
the controlled process is chosen as $X_{0}=(1,2,-1,2)^{\top}$.

(a) Parameter estimation for $\alpha_{1}$ (b) Parameter estimation for
$\alpha_{1}$ (c) Estimated optimal control

Figure 7: Example 1. Performance of the SW-DDC algorithm driven by the direct
filter (4D case)

In Figure 7, we present the performance of the SW-DDC algorithm driven by the
direct filter, and we used $500$ particles in the direct filter method to
estimate the parameters. Figure 7 (a) and (b) provide the performance of the
direct filter in estimating the unknown parameters $\alpha$ and $\beta$,
respectively, and Figure 7 (c) shows the estimation accuracy for the optimal
control. We can see from this figure that the direct filter can capture the
true target parameters quickly, and once accurate estimates for parameters are
obtained, the sample-wise optimal control solver starts to provide good
estimates for the optimal control.

As a higher dimensional problem (with 4D state and 2D parameter), we observe
that the AugPF suffers more severe degeneracy problem.

(a) RMSEs for parameter estimation (b) RMSEs for optimal control

Figure 8: Example 1. Comparison of RMSEs (4D case)

In Figure 8, we show the comparison of RMSEs between the direct filter driven
SW-DDC and the AugPF driven SW-DDC, where subplot (a) shows comparison of
accuracy in estimating the unknown parameters and subplot (b) shows comparison
of accuracy in estimating the optimal control. We can see from this figure
that although we have used $20,000$, $50,000$ and $80,000$ particles, the
AugPF provided similar parameter estimation performance, which resulted
similar performance in estimating the optimal control. This indicates that
increasing particles from $20,000$ to $80,000$ would not be enough to let
AugPF accurately capture the true parameters in this experiment, and this is
caused by the well-known “degeneracy problem” in the particle filter when
solving higher dimensional problems (see [43, 2, 39, 17]). On the other hand,
by using only $500$ particles, the direct filter provided very accurate
estimated parameters, which resulted very accurate estimates for the optimal
control in this experiment.

Table 4: Example 1. 4D computational cost comparison Method | AugPF (20K) | AugPF (50K) | AugPF (80K) | Direct filter (500)
---|---|---|---|---
CPU time (seconds) | $81.3$ | $159.8$ | $283.1$ | $58.3$

In Table 4, we present comparison of computational costs corresponding to each
experiment in Figure 8, and we can see that the direct filter method can
provide the most accurate parameter/control estimation while using the lowest
amount of CPU time.

Moreover, we want to mention that although the dimension of the problem was
increased from 2D state/1D parameter to 4D state/2D parameter, the
computational cost for the direct filter driven SW-DDC method only increased
from $1.35$ seconds/step to $1.46$ seconds/step. It indicates that our SW-DDC
algorithm is not very sensitive to dimensions, and it can alleviate the so-
called “curse of dimensionality” in moderately high dimensional problems.

###  4.2 Example 2: Data driven control of an unmanned aerial vehicle

In this example, we consider a more practical scenario and solve an unmanned
aerial vehicle maneuvering problem. The controlled state process is given by
the following nonlinear dynamical system

$\displaystyle dX_{t}$ $\displaystyle=\sin(\theta_{t})dt+\sigma dW_{t},$ (37)
$\displaystyle dY_{t}$ $\displaystyle=\cos(\theta_{t})dt+\sigma dW_{t},$
$\displaystyle dZ_{t}$ $\displaystyle=(F_{t}-mg)dt+\sigma dW_{t},$
$\displaystyle d\theta_{t}$ $\displaystyle=(\alpha_{t}-\mu
m)dt+\sigma^{2}dW_{t},$

where $(X,Y,Z)^{\top}$ describes the position of an unmanned aerial vehicle
(drone) moving in the 3D space, $\theta$ is the steering angle that controls
the horizontal moving direction of the drone, which is governed by a steering
control factor $\alpha_{t}$, and $\sigma$ is the noise that perturbs both the
motions and control actions. The unknown parameter $m$ is the total mass
(weight) of the drone with its load, $g=9.8$ is the standard gravity constant.
Then, we let $F_{t}$ be another control factor that controls the elevating
force. Therefore, the control term $u_{t}:=(\alpha_{t},F_{t})$ becomes a 2D
process. The term $\mu m$ in the equation of $d\theta_{t}$ (the fourth
equation in Eq. (37)) can be interpreted as the resistance for the steering
control factor $\alpha_{t}$, which is influenced by the weight, and we choose
$\mu=0.1$ in this example.

The cost functional that we aim to minimize is defined as

$J^{*}(u)=E\Big{[}\int_{0}^{T}\frac{1}{2}((F_{t})^{2}+(\alpha_{t})^{2})dt+10\Big{(}(X_{T}-X_{P})^{2}+(Y_{T}-Y_{P})^{2}+(Z_{T}-Z_{P})^{2}+(Y_{T}-Y_{P})^{2}\Big{)}\Big{]},$
(38)

where $(X_{P},Y_{P},Z_{P})^{\top}$ is the pre-designated target platform, and
we let $(X_{P},Y_{P},Z_{P})^{\top}=(6,7,8)^{\top}$ in our numerical
experiments. The goal of the optimal control problem is to let the controlled
drone be as close as to the target at terminal time $T$ while trying to keep
the cost of control actions low.

The observational data that we use to design data driven control actions
$\alpha_{t}$ and $F_{t}$ is

$M_{t}=(X_{t},Y_{t},Z_{t},\theta_{t})^{\top}+(R,R,R,0.1R)\xi_{t},$

where $\xi_{t}\sim N(0,I_{4})$ is a standard multi-variate Gaussian random
variable and $R$ determines the size of observation errors.

In our numerical experiments, we choose $\sigma=0.2$, $R=0.01$, and we let the
designated terminal time be $T=1$ with each observation/estimation step-size
$\Delta t=0.02$, i.e. $N_{T}=50$. The initial condition of the drone is
$(X_{0},Y_{0},Z_{0},\theta_{0})^{\top}=(0,0,5,0)^{\top}$.

(a) $X_{1}-X_{2}$ location (b) $X_{3}$ location

Figure 9: Example 2. Location of the controlled drone controlled by the SW-DDC
algorithm

In Figure 9, we present the controlled trajectory of the drone in the
$X_{1}-X_{2}$ plane (subplot (a)) and in the $X_{3}$ direction (subplot (b))
based on the control actions designed by the SW-DDC algorithm, and the
parameter estimation method was chosen as the direct filter method with $200$
particles in estimating the unknown parameter. The red curves marked by
circles are drone trajectories, and the black box (in subplot (a)) and the
black solid line (in subplot (b)) indicate the target platform in the
$X_{1}-X_{2}$ plane and the height ($X_{3}$-coordinate), respectively.

Figure 10: Example 2. 3D demonstration of the controlled drone

Also, we plot the 3D trajectory of the drone in Figure 10 to provide a more
direct viewing for the performance of our algorithm, where the black box is
the 3D position of the target platform. We can see from Figure 9 and Figure 10
that the SW-DDC can effectively guide the drone to the target platform at the
pre-designated time. The overall computational time to design altogether $50$
data driven control actions in this experiment is $23.99$ seconds, which shows
that our algorithm is potentially feasible to implement real-time data
analysis and optimal control when solving practical application problems.

Figure 11: Example 2. Parameter estimation performance of the direct filter

In Figure 11, we present the performance of parameter estimation by using the
direct filter method. We can see that after receiving only a few sets of
observational data, the direct filter (by using only $200$ particles) obtained
accurate estimation for the unknown parameter.

To further demonstrate the advantageous performance of our SW-DDC method, we
use the augmented ensemble Kalman filter (AugEnKF) (the ensemble Kalman filter
implementation for the augmented filter method) to drive the sample-wise
optimal control solver, where the AugEnKF is well known as the most widely
used method in solving the online parameter estimation problem.

(a) Parameter estimation (AugEnKF) (b) 3D trajectory of the controlled state
(AugEnKF)

Figure 12: Example 2. Performance of SW-DDC driven by AugEnKF estimated
parameter

In Figure 12, we present the parameter estimation performance of the AugEnKF
with $50$ realizations of Kalman filter samples (in subplot (a)) and the 3D
trajectory of the controlled drone (in subplot (b)), where the black solid
line in subplot (a) is the real parameter and the black box in subplot (b)
shows the target platform in the 3D space. Note that $50$ realizations of
Kalman filter samples is already a very large number of samples for the
ensemble Kalman filter in this relatively low dimensional problem. From this
figure, we can see that the AugEnKF made unreliable estimation for the unknown
parameter, which is mainly because of the nonlinearity of the dynamical model
and the nonlinear parameter-state relation. As a result, the control actions
designed for the drone are based on inaccurate parameter, which made the drone
fail to arrive at the target platform (as presented in Figure 12 (b)).

(a) RMSEs in parameter estimation (b) Average distance to the target platform

Figure 13: Example 2. Performance comparison with repeated experiments

To further compare the performance between the direct filter and the AugEnKF
in driving the sample-wise optimal control solver and show the advantage of
our SW-DDC algorithm introduced in Algorithm 1, we repeat the above
experiments $20$ times. In Figure 13 (a), we show the comparison of RMSEs in
parameter estimation, and in Figure 13 (b) we plot the average distance to the
target platform. Specifically, the red curves marked by circles are RMSEs of
parameter estimation (in subplot (a)) and average distance to the target
platform (in subplot (b)) obtained by the direct filter and the direct filter
driven SW-DDC, respectively. The blue curves marked by crosses are RMSEs of
parameter estimation (in subplot (a)) and average distance to the target
platform (in subplot (b)) obtained by the AugEnKF and the AugEnKF driven
optimal control solver, respectively. From this figure, we can see that the
direct filter significantly outperforms the AugEnKF in this experiment, and
therefore the controlled state designed by using the direct filter estimated
parameter also outperforms the controlled state driven by the AugEnKF.

##  5 Concluding remarks

In this paper, we introduced a novel sample-wise data driven control (SW-DDC)
solver to solve the stochastic optimal control problem with unknown model
parameters. In the computational framework for SW-DDC, we utilize the direct
filter as an effective online parameter estimation method to drive the design
of the optimal control through an efficient sample-wise optimal control
solver. The main novelty in our methodology is to treat samples that
characterize state variables and parameter variables in the controlled
dynamical model as pseudo-data, and then we can use stochastic approximation
to accelerate the gradient descent optimization procedure for finding the data
driven optimal control. The primary contribution of the SW-DDC framework is
that it allows us to efficiently generate accurate estimates for both the
unknown parameters and the optimal control in the online manner. Numerical
experiments are carried out to demonstrate the efficiency advantage of the
sample-wise optimal control solver and the accuracy advantage of the direct
filter, and both components combining together resulted the advantageous
overall performance of the SW-DDC solver.

## References

  * [1] Shahrokh Akhlaghi, Ning Zhou, and Zhenyu Huang.  A multi-step adaptive interpolation approach to mitigating the impact of nonlinearity on dynamic state estimation.  IEEE Transactions on Smart Grid, 9(4), 7 2018.
  * [2] C. Andrieu, A. Doucet, and R. Holenstein.  Particle markov chain monte carlo methods.  J. R. Statist. Soc. B, 72(3):269–342, 2010.
  * [3] R. Archibald, F. Bao, and X. Tu.  A direct filter method for parameter estimation.  J. Comput. Phys., 398:108871, 17, 2019.
  * [4] R. Archibald, F. Bao, and J. Yong.  A stochastic gradient descent approach for stochastic optimal control.  East Asian Journal on Applied Mathematics, to appear, 2020.
  * [5] 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.
  * [6] F. Bao, Y. Cao, and H. Chi.  Adjoint forward backward stochastic differential equations driven by jump processes and its application to nonlinear filtering problems.  International Journal of Uncertainty Quantification, 9(2):143–159, 2019.
  * [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] F. Bao, Y. Cao, and W. Zhao.  Numerical solutions for forward backward doubly stochastic differential equations and zakai equations.  International Journal for Uncertainty Quantification, 1(4):351–367, 2011.
  * [10] F. Bao, Y. Cao, and W. Zhao.  A backward doubly stochastic differential equation approach for nonlinear filtering problems.  Commun. Comput. Phys., 23(5):1573–1601, 2018.
  * [11] F. Bao, N. Cogan, A. Dobreva, and R. Paus.  Data assimilation of synthetic data as a novel strategy for predicting disease progression in alopecia areata.  Mathematical Medicine and Biology: A Journal of the IMA, 2021.
  * [12] 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.
  * [13] 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.
  * [14] Feng Bao, Yanzhao Cao, and Weidong Zhao.  A first order semi-discrete algorithm for backward doubly stochastic differential equations.  Discrete and Continuous Dynamical Systems-Series B, 5(2):1297 – 1313, 2015.
  * [15] Feng Bao and Vasileios Maroulas.  Adaptive meshfree backward SDE filter.  SIAM J. Sci. Comput., 39(6):A2664–A2683, 2017.
  * [16] R. Bellman.  Dynamic Programming.  1957\.
  * [17] A. J. Chorin and X. Tu.  Implicit sampling for particle filters.  Proc. Nat. Acad. Sc. USA, 106:17249–17254, 2009.
  * [18] O. Dyck, M. Ziatdinov, S. Jesse, F. Bao, A. Yousefzadi Nobakht, A. Maksov, B.G. Sumpter, R. Archibald, K.J.H. Law, and S.V. Kalinin.  Probing potential energy landscapes via electron-beam-induced single atom dynamics.  Acta Materialia, 203:116508, 2021.
  * [19] G. Evensen.  Data assimilation: the ensemble Kalman filter.  Springer, 2006.
  * [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] 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.
  * [23] 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.
  * [24] Ajay Jasra, Kengo Kamatani, Kody Law, and Yan Zhou.  Bayesian static parameter estimation for partially observed diffusions via multilevel Monte Carlo.  SIAM J. Sci. Comput., 40(2):A887–A902, 2018.
  * [25] S.J. Julier and J.K. Uhlmann.  Unscented filtering and nonlinear estimation.  Proceedings of the IEEE, 92:401–422, 2004.
  * [26] S. Kalinin, A. Borisevich, and S. Jesse.  Fire up the atom forge.  Nature, 22 November 2016.
  * [27] 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.
  * [28] N. Kantas, A. Doucet, S. Singh, J. Maciejowski, and N. Chopin.  On particle methods for parameter estimation in state-space models.  Statistical Science, 30(3):328 – 351, 2015.
  * [29] Dzmitry Kastsiukevich and Natalia Dmitruk.  Data-driven optimal control of linear time-invariant systems.  IFAC-PapersOnLine, 53(2):7191–7196, 2020.  21st IFAC World Congress.
  * [30] G. Kitagawa.  A self-organizing state-space model.  J. Amer. Statist. Assoc., 93:1203–1215, 1998.
  * [31] R. V. Gamkrelidze E. F. Mishchenko L. S. Pontryagin, V. G. Boltyanskii.  Power system stability and control.  1994\.
  * [32] J. Liu and M. West.  Combined parameter and state estimation in simulation-based filtering.  In Sequential Monte Carlo methods in practice, Stat. Eng. Inf. Sci., pages 197–223. Springer, New York, 2001.
  * [33] Matthias Morzfeld, Marcus S. Day, Ray W. Grout, George Shu Heng Pau, Stefan A. Finsterle, and John B. Bell.  Iterative importance sampling algorithms for parameter estimation.  SIAM J. Sci. Comput., 40(2):B329–B352, 2018.
  * [34] Ali Yousefzadi Nobakht, Ondrej Dyck, David B Lingerfelt, Feng Bao, Maxim Ziatdinov, Artem Maksov, Bobby G Sumpter, Richard Archibald, and Sergei V Kalinin Stephen Jesse, and Kody JH Law.  Reconstruction of effective potential from statistical analysis of dynamic trajectories.  AIP Advances, 10:065034, 2020.
  * [35] D. Oliver, A. Reynolds, and N. Liu.  Inverse theory for petroleum reservoir characterization and history matching.  Cambridge University Press, 2008.
  * [36] N. J. Balu P. Kundur and M. G. Lauby.  The Mathematical Theory of Optimal Processes.  1962\.
  * [37] É. Pardoux and S. Peng.  Backward stochastic differential equations and quasilinear parabolic partial differential equations.  In Stochastic partial differential equations and their applications (Charlotte, NC, 1991), volume 176 of Lecture Notes in Control and Inform. Sci., pages 200–217. Springer, Berlin, 1992.
  * [38] S. Peng.  A general stochastic maximum principle for optimal control problems.  SIAM J. Control Optim., pages 966–979, 1990.
  * [39] Michael K. Pitt and Neil Shephard.  Filtering via simulation: auxiliary particle filters.  J. Amer. Statist. Assoc., 94(446):590–599, 1999.
  * [40] Mario Sassano and Alessandro Astolfi.  A local separation principle via dynamic approximate feedback and observer linearization for a class of nonlinear systems.  IEEE Trans. Automat. Control, 64(1):111–126, 2019.
  * [41] 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.
  * [42] A. Shapiro and Y. Wardi.  Convergence analysis of gradient descent stochastic algorithms.  Journal of Optimization Theory and Applications, pages 439–454, 1996.
  * [43] C. Snyder, T. Bengtsson, P. Bickel, and J. Anderson.  Obstacles to high-dimensional particle filtering.  Mon. Wea. Rev., 136:4629–4640, 2008.
  * [44] Shanjian Tang.  The maximum principle for partially observed optimal control of stochastic differential equations.  SIAM J. Control Optim., 36(5):1596–1617, 1998.
  * [45] S. Wilks.  Mathematical Statistics.  New York: John Wiley & Sons, 1962.
  * [46] W. M. Wonham.  On the separation theorem of stochastic control.  SIAM J. Control, 6:312–326, 1968.
  * [47] 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.
  * [48] Jianfeng Zhang.  A numerical scheme for BSDEs.  Ann. Appl. Probab., 14(1):459–488, 2004.
  * [49] 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.