be the unique $\mathbb{B}$-path such that $r$ corresponds to $\rho$.
Similarly, recall that $\Lambda^{\hat{\pi}}({\rho})$ is the set of finite
paths of $\hat{\pi}$ that correspond to $\rho$ (where the notion of
correspondence of finite paths of $\hat{\pi}$ to a $\mathbb{B}$-path is given
in the proof of Lemma 5). Consider the extension of $\rho$ with
$(\tilde{B},p)$ such that $\mathit{source}({p})=l_{n}$,
$\tilde{B}\in\mathcal{I}_{\mathbb{B}}$, $\tilde{B}\geq B_{n}$, and
$\tilde{B}\models\mathit{guard}({p})\wedge\mathit{inv}(l_{n})$.

We recall an assumption that we can make without loss of generality, as
already described in the proof of Lemma 2: given
$r^{\prime}\in\Lambda^{\hat{\pi}}({\rho})$, we assume that there exists at
most one
$\alpha\in\mathfrak{G}({\mathfrak{D}_{\mathfrak{M}[{\mathcal{P}}]}(\mathit{last}(r^{\prime}),(\tilde{B},p))})$
such that $((\tilde{B},p),\alpha)\in\mathsf{support}(\hat{\pi}(r^{\prime}))$
(this follows from the fact that the set of assignments
$\mathfrak{G}({\mathfrak{D}_{\mathfrak{M}[{\mathcal{P}}]}(\mathit{last}(r^{\prime}),(\tilde{B},p))})$
is closed under convex combinations). If such
$\alpha\in\mathfrak{G}({\mathfrak{D}_{\mathfrak{M}[{\mathcal{P}}]}(\mathit{last}(r^{\prime}),(\tilde{B},p))})$
exists, we denote it by $\alpha^{r^{\prime}}_{({\tilde{B}},{p})}$.

As in the proof of Lemma 4, to define $\pi(r)$, for each pair $(\tilde{B},p)$
used to extend $\rho$, we identify a clock valuation
$\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})}$, and then let:

$\displaystyle\pi(r)(\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})},p)$
$\displaystyle=$
$\displaystyle\frac{\sum_{r^{\prime}\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r^{\prime})\cdot\hat{\pi}(r^{\prime})((\tilde{B},p),\alpha^{r^{\prime}}_{({\tilde{B}},{p})})}{\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(\Lambda^{\hat{\pi}}({\rho}))}\;,$

where
$\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(\Lambda^{\hat{\pi}}({\rho}))=\sum_{r^{\prime}\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r^{\prime})$,
and where
$\hat{\pi}(r^{\prime})((\tilde{B},p),\alpha^{r^{\prime}}_{({\tilde{B}},{p})})=0$
if there does not exist any
$\alpha\in\mathfrak{G}({\mathfrak{D}_{\mathfrak{M}[{\mathcal{P}}]}(\mathit{last}(r^{\prime}),(\tilde{B},p))})$
such that $((\tilde{B},p),\alpha)\in\mathsf{support}(\hat{\pi}(r^{\prime}))$.

For the case in which $\tilde{B}$ is a closed interval, i.e.,
$\tilde{B}=[b,b]$ for $b\in\mathbb{B}$, we let
$\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})}=b$.

We now consider the case in which $\tilde{B}$ is an open interval. The
approach we take is analogous to that taken in the proof of Lemma 4:
intuitively, if $p$ is a constant probabilistic edge, and no non-constant
probabilistic edge has been taken while the clock value remains in
$\tilde{B}$, then the scheduler $\pi$ maintains the clock as having a value
that is a (small) amount greater than the left endpoint of $\tilde{B}$; if $p$
is instead a non-constant probabilistic edge, then the value of the clock is
advanced to a particular value in the interval $\tilde{B}$, in order to mimic
the choice of particular assignments in $\hat{\pi}$; finally, if $p$ is a
constant probabilistic edge, and a non-constant probabilistic edge _has_ been
taken previously while the clock value remains in $\tilde{B}$, then the
scheduler $\pi$ can choose the value of the clock from $\tilde{B}$ in an
arbitrary manner. As in the proof of Lemma 4, the fact that we need to
consider only these three cases relies on the assumption of initialisation.
The construction is more complicated than that of the proof of Lemma 4,
because the clock valuations must be chosen by $\pi$ to reflect _assignments_
chosen by $\hat{\pi}$ (instead, for Lemma 4, the clock valuations were chosen
in the construction of $\pi$ to reflect _clock valuations_ chosen by $\sigma$,
which was simpler conceptually). As in the proof of Lemma 4, let $k\leq n$ be
the maximum index for which $B_{k}\neq\tilde{B}$ (with $k=0$ if no such index
exists), noting that
$\tilde{B}_{k}=B_{k+1}=\tilde{B}_{k+1}=\cdots=\tilde{B}_{n-1}=B_{n}=\tilde{B}$.
We consider the following function
$h_{\tilde{B}}:\bigcup_{p^{\prime}\in\mathit{prob}}\mathsf{Dist}(\\{{(\tilde{B},p^{\prime},\mathsf{le}),(\tilde{B},p^{\prime},\mathsf{re})}\\})\rightarrow\tilde{B}$
that maps an assignment (over the target endpoint indicators
$(\tilde{B},p^{\prime},\mathsf{le}),(\tilde{B},p^{\prime},\mathsf{re})$ for
some $p^{\prime}\in\mathit{prob}$) to a clock valuation in $\tilde{B}$. The
function $h_{\tilde{B}}$ is defined as follows: given
$\alpha\in\mathsf{Dist}(\\{{(\tilde{B},p^{\prime},\mathsf{le}),(\tilde{B},p^{\prime},\mathsf{re})}\\})$,
we let
$h_{\tilde{B}}(\alpha)=\mathsf{re}(\tilde{B})-\alpha(\mathsf{le})\cdot(\mathsf{re}(\tilde{B})-\mathsf{le}(\tilde{B}))$
(note that, as in the proof of Lemma 5, we write $\alpha(\mathsf{le})$ rather
than $\alpha(\tilde{B},p^{\prime},\mathsf{le})$ for simplicity; in the sequel,
we will also usually write $\alpha(\mathsf{re})$ rather than
$\alpha(\tilde{B},p^{\prime},\mathsf{re})$). The intuition underlying the
definition of $h_{\tilde{B}}$ is that $h_{\tilde{B}}(\alpha)$ is a clock
valuation that represents faithfully the assignment $\alpha$, in the sense
that the greater the probability of $\alpha(\mathsf{le})$, the closer the
clock valuation is to the left endpoint of $\tilde{B}$.

First we consider the subcase in which all probabilistic edges in $p_{k}\cdots
p_{n}$ are constant, and $p$ is also constant. While the clock must be set to
some value in $\tilde{B}$, it is important that the value chosen is not
greater than the value needed in the future to replicate the probabilities
corresponding to a non-constant probabilistic edge. We now identify such
values. We first define the set $\Omega_{{\rho},({\tilde{B}},{p})}$, which
contains finite suffixes of paths in $\Lambda^{\hat{\pi}}({\rho})$, where the
suffixes (1) start with a transition derived from $(\tilde{B},p)$, (2) feature
only $\tilde{B}$ (apart from the suffixes’ first state), and (3) features only
constant probabilistic edges until the final transition, which features a non-
constant probabilistic edge. We now define formally
$\Omega_{{\rho},({\tilde{B}},{p})}$. For a finite path
$r^{\prime}\in\mathit{Paths}^{\hat{\pi}}_{*,{S_{\mathfrak{M}[{\mathcal{P}}]}^{\mathsf{reg}}}}(\mathsf{reg}({l,0}))$,
let $\mathit{SuffFrom}_{*}^{\hat{\pi}}({r^{\prime}})$ be the set of (finite)
suffixes of $r^{\prime}$ generated by $\hat{\pi}$ that terminate in a state in
$S_{\mathfrak{M}[{\mathcal{P}}]}^{\mathsf{end}}$; formally,
$\mathit{SuffFrom}_{*}^{\hat{\pi}}({r^{\prime}})$ is the smallest set such
that, if
$r^{\prime\prime}\in\mathit{Paths}^{\hat{\pi}}_{*,{S_{\mathfrak{M}[{\mathcal{P}}]}^{\mathsf{end}}}}(\mathsf{reg}({l,0}))$,
where $r^{\prime\prime}=r^{\prime}r^{\prime\prime\prime}$, then
$r^{\prime\prime\prime}\in\mathit{SuffFrom}_{*}^{\hat{\pi}}({r^{\prime}})$. We
let $\mathit{Suffixes}_{*}^{\hat{\pi}}(\mathsf{reg}({l,0}))$ be the set of all
suffixes of finite paths of $\hat{\pi}$, starting in a state of
$S_{\mathfrak{M}[{\mathcal{P}}]}^{\mathsf{reg}}$ and terminating in a state in
$S_{\mathfrak{M}[{\mathcal{P}}]}^{\mathsf{end}}$, that is:

$\mathit{Suffixes}_{*}^{\hat{\pi}}(\mathsf{reg}({l,0}))=\bigcup_{r^{\prime}\in\mathit{Paths}^{\hat{\pi}}_{*,{S_{\mathfrak{M}[{\mathcal{P}}]}^{\mathsf{reg}}}}(\mathsf{reg}({l,0}))}\mathit{SuffFrom}_{*}^{\hat{\pi}}({r^{\prime}})\;.$

Consider
$r^{\prime}\in\mathit{Suffixes}_{*}^{\hat{\pi}}(\mathsf{reg}({l,0}))$, where

$r^{\prime}=(l_{0},B_{0})((\tilde{B}_{0},p_{0}),\alpha_{0})(\mathsf{ep}_{0})\cdots\linebreak[0]((\tilde{B}_{m},p_{m}),\alpha_{m})(\mathsf{ep}_{m}).$

We let $\mathsf{FA}({r^{\prime}})=\alpha_{m}$ be the assignment featured in
the final transition of $r^{\prime}$. Furthermore, we say that:

  * •

$r^{\prime}$ _starts with $(\tilde{B},p)$_ if $\tilde{B}_{0}=\tilde{B}$ and
$p_{0}=p$;

  * •

$r^{\prime}$ _remains in $\tilde{B}$_ if $\tilde{B}_{i}=\tilde{B}$ for all
$0\leq i\leq m$ and $B_{i}=\tilde{B}$ for all $0<i\leq m$;

  * •

$r^{\prime}$ _terminates with a non-constant probabilistic edge_ if $p_{m}$ is
a non-constant probabilistic edge and $p_{i}$ is a constant probabilistic edge
for all $0\leq i<m$.

We say that
$r^{\prime}\in\mathit{Suffixes}_{*}^{\hat{\pi}}(\mathsf{reg}({l,0}))$ is a
_$(\tilde{B},p)$ -critical path_ if $r^{\prime}$ starts with $(\tilde{B},p)$,
remains in $\tilde{B}$, and terminates with a non-constant probabilistic edge.
Let $\mathit{Critical}_{(\tilde{B},p)}^{\hat{\pi}}$ be the set of
$(\tilde{B},p)$-critical paths. We can now let:

$\Omega_{{\rho},({\tilde{B}},{p})}=\mathit{Critical}_{(\tilde{B},p)}^{\hat{\pi}}\cap\bigcup_{r^{\prime}\in\Lambda^{\hat{\pi}}({\rho})}\mathit{SuffFrom}_{*}^{\hat{\pi}}({r^{\prime}})\;.$

Then we let $\Psi_{{\rho},({\tilde{B}},{p})}$ be the set of assignments
featured in the final transition of finite paths in
$\Omega_{{\rho},({\tilde{B}},{p})}$: formally
$\Psi_{{\rho},({\tilde{B}},{p})}=\\{{\mathsf{FA}({r^{\prime}})\mid
r^{\prime}\in\Omega_{{\rho},({\tilde{B}},{p})}}\\}$. Finally we let
$\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})}=\min\\{{h_{\tilde{B}}(\alpha)\mid\alpha\in\Psi_{{\rho},({\tilde{B}},{p})}}\\}$.

Now consider the subcase in which $\tilde{B}$ is an open interval, and all
probabilistic edges in the sequence $p_{k}\cdots p_{n}$ are constant, but $p$
is non-constant. As in the analogous case of the proof of Lemma 4, the choice
of clock valuation by $\pi$ corresponds to a weighted average from a set of
clock valuations. In the case of this proof, the set of clock valuations used
is that obtained from assignments used by $\hat{\pi}$ after finite paths
corresponding to $\rho$ extended with $(\tilde{B},p)$, using a similar
“transformation from assignment to clock valuation” approach as used in the
previous paragraph. We now give the formal details. Recall that we have
assumed w.l.o.g. that, given $r^{\prime}\in\Lambda^{\hat{\pi}}({\rho})$, there
exists at most one
$\alpha\in\mathfrak{G}({\mathfrak{D}_{\mathfrak{M}[{\mathcal{P}}]}(\mathit{last}(r^{\prime}),(\tilde{B},p))})$
such that $((\tilde{B},p),\alpha)\in\mathsf{support}(\hat{\pi}(r^{\prime}))$,
and that we denote such an assignment $\alpha$ by
$\alpha^{r^{\prime}}_{({\tilde{B}},{p})}$. Then we let:

$\displaystyle\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})}$
$\displaystyle=$
$\displaystyle\frac{\sum_{r^{\prime}\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,(\mathsf{reg}({l,0}))}(r^{\prime})\cdot\hat{\pi}(r^{\prime})((\tilde{B},p),\alpha^{r^{\prime}}_{({\tilde{B}},{p})})\cdot
h_{\tilde{B}}(\alpha^{r^{\prime}}_{({\tilde{B}},{p})})}{\sum_{r^{\prime}\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,(\mathsf{reg}({l,0}))}(r^{\prime})\cdot\hat{\pi}(r^{\prime})((\tilde{B},p),\alpha^{r^{\prime}}_{({\tilde{B}},{p})})}\;.$

Finally, we consider the subcase in which $\tilde{B}$ is an open interval,
there exists $i\in\\{{k,\ldots,n}\\}$ such that $p_{i}$ is non-constant, and
for all other $j\in\\{{k,\ldots,n}\\}$ we have that $p_{i}$ is constant, and
also $p$ is constant. As in the proof of Lemma 4, for this subcase, the clock
valuation can be arbitrary; for simplicity we retain the same clock valuation
as was used in the final state of $r$, i.e., if $\mathit{last}(r)$ is equal to
$(l^{\prime},v^{\prime})$, we let
$\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})}=v^{\prime}$.

Repeating this overall process for all
$(\tilde{B},p)\in\mathfrak{A}_{\mathfrak{M}[{\mathcal{P}}]}({l,B})$ suffices
to define the distribution $\pi(r)$. Following similar reasoning to the
analogous part of Lemma 4, we can show that $\pi(r)$ is indeed a distribution,
i.e.,
$\sum_{(\tilde{B},p)\in\mathfrak{A}_{\mathfrak{M}[{\mathcal{P}}]}({\mathsf{reg}({\mathit{last}(r)})})}\pi(r)(\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})},p)=1$
(recall that $\mathsf{reg}({\mathit{last}(r)})=(l_{n},B_{n})$). From the
construction of $\pi$, and from the fact that
$\sum_{(\tilde{B},p)\in\mathfrak{A}_{\mathfrak{M}[{\mathcal{P}}]}({l_{n},B_{n}})}\hat{\pi}(r^{\prime})((\tilde{B},p),\alpha^{r^{\prime}}_{({\tilde{B}},{p})})=1$,
we have the following:

$\displaystyle\sum_{(\tilde{B},p)\in\mathfrak{A}_{\mathfrak{M}[{\mathcal{P}}]}({l_{n},B_{n}})}\pi(r)(\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})},p)$
$\displaystyle=$
$\displaystyle\sum_{(\tilde{B},p)\in\mathfrak{A}_{\mathfrak{M}[{\mathcal{P}}]}({l_{n},B_{n}})}\left(\frac{\sum_{r^{\prime}\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r^{\prime})\cdot\hat{\pi}(r^{\prime})((\tilde{B},p),\alpha^{r^{\prime}}_{({\tilde{B}},{p})})}{\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(\Lambda^{\hat{\pi}}({\rho}))}\right)$
$\displaystyle=$
$\displaystyle\frac{\sum_{r^{\prime}\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r^{\prime})\cdot\sum_{(\tilde{B},p)\in\mathfrak{A}_{\mathfrak{M}[{\mathcal{P}}]}({l_{n},B_{n}})}\hat{\pi}(r^{\prime})((\tilde{B},p),\alpha^{r^{\prime}}_{({\tilde{B}},{p})})}{\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(\Lambda^{\hat{\pi}}({\rho}))}$
$\displaystyle=$
$\displaystyle\frac{\sum_{r^{\prime}\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r^{\prime})}{\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(\Lambda^{\hat{\pi}}({\rho}))}$
$\displaystyle=$ $\displaystyle 1\;.$

The next step of the proof is to show that $\mathrm{Pr}^{\pi}_{(l,0)}(\Diamond
T_{F})=\mathrm{Pr}^{\hat{\pi}}_{\mathsf{reg}({l,0})}(\Diamond\mathfrak{T}_{F})$.
As in previous proofs, it is sufficient to show that
$\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho}))=\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(\Lambda^{\hat{\pi}}({\rho}))$
for all $\mathbb{B}$-paths $\rho$. Proceeding by induction on the length of
$\mathbb{B}$-paths, consider the $\mathbb{B}$-path
$\rho(\tilde{B},p)(l^{\prime},B^{\prime})$, and assume that
$\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho}))=\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(\Lambda^{\hat{\pi}}({\rho}))$.
Our aim is to show that

$\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho(\tilde{B},p)(l^{\prime},B^{\prime})}))=\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(\Lambda^{\hat{\pi}}({\rho(\tilde{B},p)(l^{\prime},B^{\prime})}))\;.$

We consider the case in which $\tilde{B}$ is open, $p$ is non-constant and
$B^{\prime}=\tilde{B}$, i.e., the target region $(l^{\prime},B^{\prime})$ is
obtained using the outcome $(\emptyset,l^{\prime})$ (which does not reset the
clock, and hence the interval $B^{\prime}$ of the target region is equal to
the interval $\tilde{B}$ obtained just before the probabilistic edge $p$ is
taken). Other cases are dealt with similarly; in particular, the cases for
constant probabilistic edges are more straightforward. Our approach is
slightly different from that used in the proof of Lemma 5 because, in the
scheduler of ${[[\mathfrak{M}[{\mathcal{P}}]]]}$ constructed in the proof of
Lemma 5, the choice made after all finite paths corresponding to a particular
$\mathbb{B}$-path was the same; instead, this property does not necessarily
hold for $\hat{\pi}$. Note that in this final part of the proof, the fact that
clock dependencies are affine will be used at multiple points. Let
$r\in\Lambda^{\hat{\pi}}({\rho})$. Now consider the finite paths in
$\Lambda^{\hat{\pi}}({\rho(\tilde{B},p)(l^{\prime},B^{\prime})})$ that have
$r$ as a prefix: these finite paths are
$r((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})(\mathsf{le})(l^{\prime},B^{\prime})$
and
$r((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})(\mathsf{re})(l^{\prime},B^{\prime})$.
From the construction of $\hat{\pi}$ and by definition (where $\wp$ denotes
the distribution template of probabilistic edge $p$, and $c$ and $d$ denote
$c^{p}_{(\emptyset,l^{\prime})}$ and $d^{p}_{(\emptyset,l^{\prime})}$,
respectively), we have the following:

$\displaystyle\sum_{\mathsf{ep}\in\\{{\mathsf{le},\mathsf{re}}\\}}\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})(\mathsf{ep})(l^{\prime},B^{\prime}))$
$\displaystyle=$
$\displaystyle\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})$
$\displaystyle\times\sum_{\mathsf{ep}\in\\{{\mathsf{le},\mathsf{re}}\\}}\alpha^{r}_{({\tilde{B}},{p})}(\tilde{B},p,\mathsf{ep})\cdot\Delta_{\mathfrak{M}[{\mathcal{P}}]}((\tilde{B},p,\mathsf{ep}),(\tau,\alpha^{\tilde{B},p,\mathsf{ep}}))(l^{\prime},B^{\prime})$
$\displaystyle=$
$\displaystyle\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})\cdot\sum_{\mathsf{ep}\in\\{{\mathsf{le},\mathsf{re}}\\}}\alpha^{r}_{({\tilde{B}},{p})}(\tilde{B},p,\mathsf{ep})\cdot{\wp}[{\mathsf{ep}(\tilde{B})}](\emptyset,l^{\prime})$
$\displaystyle=$
$\displaystyle\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})\cdot\sum_{\mathsf{ep}\in\\{{\mathsf{le},\mathsf{re}}\\}}\alpha^{r}_{({\tilde{B}},{p})}(\tilde{B},p,\mathsf{ep})\cdot(c+d\cdot\mathsf{ep}(\tilde{B}))$
$\displaystyle=$ $\displaystyle
c\cdot\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})\cdot\sum_{\mathsf{ep}\in\\{{\mathsf{le},\mathsf{re}}\\}}\alpha^{r}_{({\tilde{B}},{p})}(\tilde{B},p,\mathsf{ep})$
$\displaystyle+\,d\cdot\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})\cdot\sum_{\mathsf{ep}\in\\{{\mathsf{le},\mathsf{re}}\\}}\alpha^{r}_{({\tilde{B}},{p})}(\tilde{B},p,\mathsf{ep})\cdot\mathsf{ep}(\tilde{B})\;.$

Note that
$\sum_{\mathsf{ep}\in\\{{\mathsf{le},\mathsf{re}}\\}}\alpha^{r}_{({\tilde{B}},{p})}(\tilde{B},p,\mathsf{ep})=1$
(because $\alpha^{r}_{({\tilde{B}},{p})}$ is a distribution over the set
$\\{{(\tilde{B},p,\mathsf{le}),(\tilde{B},p,\mathsf{re})}\\}$). Furthermore,
recalling the definition of $h_{\tilde{B}}$, we observe that:

$\displaystyle\alpha^{r}_{({\tilde{B}},{p})}(\tilde{B},p,\mathsf{le})$
$\displaystyle=$
$\displaystyle\frac{\mathsf{re}(\tilde{B})-h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})}{\mathsf{re}(\tilde{B})-\mathsf{le}(\tilde{B})}$
$\displaystyle\alpha^{r}_{({\tilde{B}},{p})}(\tilde{B},p,\mathsf{re})$
$\displaystyle=$
$\displaystyle\frac{h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})-\mathsf{le}(\tilde{B})}{\mathsf{re}(\tilde{B})-\mathsf{le}(\tilde{B})}\;.$

Hence we have:

$\displaystyle\sum_{\mathsf{ep}\in\\{{\mathsf{le},\mathsf{re}}\\}}\alpha^{r}_{({\tilde{B}},{p})}(\tilde{B},p,\mathsf{ep})\cdot\mathsf{ep}(\tilde{B})$
$\displaystyle=$
$\displaystyle\frac{\mathsf{re}(\tilde{B})-h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})}{\mathsf{re}(\tilde{B})-\mathsf{le}(\tilde{B})}\cdot\mathsf{le}(\tilde{B})+\frac{h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})-\mathsf{le}(\tilde{B})}{\mathsf{re}(\tilde{B})-\mathsf{le}(\tilde{B})}\cdot\mathsf{re}(\tilde{B})$
$\displaystyle=$
$\displaystyle\frac{\mathsf{re}(\tilde{B})\cdot\mathsf{le}(\tilde{B})-h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})\cdot\mathsf{le}(\tilde{B})+h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})\cdot\mathsf{re}(\tilde{B})-\mathsf{le}(\tilde{B})\cdot\mathsf{re}(\tilde{B})}{\mathsf{re}(\tilde{B})-\mathsf{le}(\tilde{B})}$
$\displaystyle=$
$\displaystyle\frac{h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})\cdot\mathsf{re}(\tilde{B})-h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})\cdot\mathsf{le}(\tilde{B})}{\mathsf{re}(\tilde{B})-\mathsf{le}(\tilde{B})}$
$\displaystyle=$ $\displaystyle
h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})\;.$

Therefore we have:

$\displaystyle
c\cdot\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})\cdot\sum_{\mathsf{ep}\in\\{{\mathsf{le},\mathsf{re}}\\}}\alpha^{r}_{({\tilde{B}},{p})}(\tilde{B},p,\mathsf{ep})$
$\displaystyle+\,d\cdot\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})\cdot\sum_{\mathsf{ep}\in\\{{\mathsf{le},\mathsf{re}}\\}}\alpha^{r}_{({\tilde{B}},{p})}(\tilde{B},p,\mathsf{ep})\cdot\mathsf{ep}(\tilde{B})$
$\displaystyle=$ $\displaystyle
c\cdot\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})$
$\displaystyle+\,d\cdot\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})\cdot
h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})\;.$

We now use this equality to establish that
$\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho}))=\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(\Lambda^{\hat{\pi}}({\rho}))$.
First observe the following:

$\displaystyle\mathrm{Pr}^{\hat{\pi}}_{*,(l,0)}(\Lambda^{\pi}({\rho(\tilde{B},p)})(l^{\prime},B^{\prime}))$
$\displaystyle=$
$\displaystyle\sum_{r\in\Lambda^{\hat{\pi}}({\rho(\tilde{B},p)})(l^{\prime},B^{\prime})}\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)$
$\displaystyle=$
$\displaystyle\sum_{r\in\Lambda^{\hat{\pi}}({\rho})}\sum_{\mathsf{ep}\in\\{{\mathsf{le},\mathsf{re}}\\}}\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})(\mathsf{ep})(l^{\prime},B^{\prime}))$
$\displaystyle=$
$\displaystyle\sum_{r\in\Lambda^{\hat{\pi}}({\rho})}c\cdot\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})$
$\displaystyle+\,\sum_{r\in\Lambda^{\hat{\pi}}({\rho})}d\cdot\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})\cdot
h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})$

Recalling the definition of
$\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})}$ (for the case of open
$\tilde{B}$ and non-constant $p$), we have that
$\sum_{r\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,(\mathsf{reg}({l,0}))}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})\cdot
h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})=\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})}\cdot\sum_{r\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,(\mathsf{reg}({l,0}))}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})$.
Hence we then obtain:

$\displaystyle\sum_{r\in\Lambda^{\hat{\pi}}({\rho})}c\cdot\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})$
$\displaystyle+\,\sum_{r\in\Lambda^{\hat{\pi}}({\rho})}d\cdot\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})\cdot
h_{\tilde{B}}(\alpha^{r}_{({\tilde{B}},{p})})$ $\displaystyle=$ $\displaystyle
c\cdot\sum_{r\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})$
$\displaystyle+\,d\cdot\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})}\cdot\sum_{r\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,(\mathsf{reg}({l,0}))}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})\;.$

Let $\check{r}$ denote the unique path in $\Lambda^{\pi}({\rho})$. By the
definition of $\pi(\check{r})$, and from the fact that
$\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho}))=\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(\Lambda^{\hat{\pi}}({\rho}))$
by induction, we have:

$\displaystyle\sum_{r\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})$
$\displaystyle=$
$\displaystyle\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho}))\cdot\pi(\check{r})(\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})},p)\;.$

We then obtain:

$\displaystyle
c\cdot\sum_{r\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})$
$\displaystyle+\,d\cdot\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})}\cdot\sum_{r\in\Lambda^{\hat{\pi}}({\rho})}\mathrm{Pr}^{\hat{\pi}}_{*,(\mathsf{reg}({l,0}))}(r)\cdot\hat{\pi}(r)((\tilde{B},p),\alpha^{r}_{({\tilde{B}},{p})})$
$\displaystyle=$ $\displaystyle
c\cdot\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho}))\cdot\pi(\check{r})(\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})},p)$
$\displaystyle+\,d\cdot\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})}\cdot\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho}))\cdot\pi(\check{r})(\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})},p)$
$\displaystyle=$
$\displaystyle\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho}))\cdot\pi(\check{r})(\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})},p)\cdot(c+d\cdot\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})})$
$\displaystyle=$
$\displaystyle\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho}))\cdot\pi(\check{r})(\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})},p)\cdot{\wp}[{\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})}}](\emptyset,l^{\prime})$
$\displaystyle=$
$\displaystyle\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho}))\cdot\pi(\check{r})(\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})},p)\cdot\Delta_{\mathcal{P}}(\mathit{last}(\check{r}),(\tilde{v}^{\mathrm{*}}_{{\rho},({\tilde{B}},{p})},p))(l^{\prime},B^{\prime})$
$\displaystyle=$
$\displaystyle\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho(\tilde{B},p)(l^{\prime},B^{\prime})}))\;.$

This establishes that
$\mathrm{Pr}^{\pi}_{*,(l,0)}(\Lambda^{\pi}({\rho(\tilde{B},p)(l^{\prime},B^{\prime})}))=\mathrm{Pr}^{\hat{\pi}}_{*,\mathsf{reg}({l,0})}(\Lambda^{\hat{\pi}}({\rho(\tilde{B},p)(l^{\prime},B^{\prime})}))$.
We note that the cases for other outcomes are similar, and for constant
probabilistic edges and/or closed intervals are more straightforward.

We characterise the size of a 1c-cdPTA as the sum of the number of its
locations, the size of the binary encoding of the clock constraints used in
invariant conditions and guards, and the size of the binary encoding of the
constants used in the distribution templates of the probabilistic edges (i.e.,
$c^{p}_{e}$ and $d^{p}_{e}$ for each $p\in\mathit{prob}$ and $e\in
2^{\\{{x}\\}}\times L$).

######  Theorem 7.

Quantitative and qualitative problems for 1c-cdPTA can be solved in polynomial
time.

The theorem follows from Lemma 4, Lemma 5, Lemma 6, Proposition 1, the fact
that the IMDP defined in this section can be constructed in polynomial time,
and the fact that quantitative and qualitative problems for IMDPs can be
solved in polynomial time, given that there exist polynomial-time algorithms
for analogous problems on IMCs with the semantics adopted in this paper
[PLSS13, CHK13, CK15, Spr18]. We add that the quantitative and qualitative
problems for 1c-cdPTAs are PTIME-hard, following from the PTIME-hardness of
the corresponding problems for MDPs [PT87, CDH10].

##  5\. Conclusion

We have presented a method for the transformation of a class of 1c-cdPTAs to
IMDPs such that there is a precise relationship between the schedulers of the
1c-cdPTA and the IMDP, allowing us to use established polynomial-time
algorithms for IMDPs to decide quantitative and qualitative reachability
problems on the 1c-cdPTA. Overall, the results establish that such problems
are in PTIME. The techniques rely on the initialisation requirement, which
ensures that optimal choices for non-constant probabilistic edges correspond
to the left or right endpoints of intervals that are derived from the
syntactic description of the 1c-cdPTA. The initialisation requirement
restricts dependencies between non-constant probabilistic edges: while this
necessarily restricts the expressiveness of the formalism, the resulting model
nevertheless retains the expressive power to represent basic situations in
which the probability of certain events depends on the exact amount of time
elapsed, such as those described in the introduction.

The IMDP construction can be simplified in a number of cases: for example, in
the case in which at most two outcomes $e_{1},e_{2}$ of every probabilistic
edge $p$ are non-constant, i.e., for which $d^{p}_{e_{1}}\neq 0$ and
$d^{p}_{e_{2}}\neq 0$, endpoint indicators are unnecessary; instead, when a
probabilistic edge is taken from an open interval $\tilde{B}$, each of $e_{1}$
and $e_{2}$ are associated with (non-singleton) intervals (other outcomes are
associated with singleton intervals), and the choice of probability to assign
between the two intervals represents the choice of clock valuation in
$\tilde{B}$. This construction is also polynomial in the size of the 1c-cdPTA.

Future work could consider time-bounded reachability problems for 1c-cdPTAs,
or lifting one of the two restrictions that we have applied on the original
clock-dependent probabilistic timed automata formalism as presented in
[Spr21], namely initialisation and the restriction to one clock. We discuss
these restrictions in turn. Consider first the case of initialisation. Observe
that initialisation enforces that the choice of clock value at the point at
which a non-constant probabilistic edge is taken is independent of the choice
of clock value at the point at which any other non-constant probabilistic edge
is taken. For example, for the 1c-cdPTA of Figure 1, the clock values chosen
when leaving location $\mathrm{W}$ and (in the case in which the outcome to
location $\mathrm{F}$ is taken) when leaving location $\mathrm{F}$ are
independent from each other, because they belong to different intervals of
$\mathcal{I}_{\mathbb{B}}$. Instead, in the example of Figure 5 (taken from
[Spr21]), in which initialisation does not hold, the choice of clock value
when leaving location $\mathrm{B}$ ($\mathrm{C}$, respectively) depends on
(i.e., cannot be less than) the clock value when leaving location $\mathrm{A}$
($\mathrm{B}$, respectively). The characteristic of independence of the choice
of clock values for different non-constant probabilistic edges permits the
reduction from 1c-cdPTAs to IMDPs presented in this paper: the IMDP
construction encodes the choice of clock value of the 1c-cdPTA by the chosen
assignment of an interval distribution (i.e., the interval distributions with
intervals $(0,1)$ available from regions). If the choice of clock values for
different non-constant probabilistic edges are _not_ independent, as is the
case when the 1c-cdPTA is not initialised, then we would require some
mechanism in the IMDP to enforce such dependencies between the choices of
assignment in different IMDP states, which does not exist in the classical
IMDP formalism. Solutions to the qualitative problem for non-initialised
1c-cdPTAs could potentially utilise connections with parametric MDPs [HHZ11,
WJPK21], in which dependence between parameters on probabilities of
transitions from different states is an inherent part of the formalism.

$\mathrm{A}$$x<1$$\mathrm{B}$$x<1$$\mathrm{C}$$x<1$$\mathrm{D}$$\mathrm{E}$$x>0$$x$$x>0$$1-x$$x$$x>0$$1-\frac{x}{2}$$1-x$$\frac{x}{2}$
Figure 5. A non-initialised 1c-cdPTA (taken from [Spr21]).

For the case of lifting the restriction to one clock, while retaining
initialisation (i.e., a clock has a natural-numbered value between any two
non-constant probabilistic edges with a dependence on that clock), a natural
way of generalising the results of this paper would be to consider classical
multidimensional regions (as presented for timed automata in [AD94] and used
for approximate analysis of clock-dependent probabilistic timed automata in
[Spr21]) as states of a finite-state MDP (whether an IMDP or a standard MDP).
Now consider a clock-dependent probabilistic timed automaton that has the same
overall structure as the 1c-cdPTA of Figure 5, but which has three clocks,
$x$, $y$ and $z$, and for which the clock dependencies of the three
probabilistic edges are as follows: for the probabilistic edge from location
$\mathrm{A}$ ($\mathrm{B}$, $\mathrm{C}$, respectively), the outcome to
location $\mathrm{B}$ ($\mathrm{C}$, $\mathrm{D}$, respectively) has clock
dependence $x$ ($1-y$, $1-\frac{z}{2}$, respectively). This example satisfies
the notion of initialisation in the case of multiple clocks, but features the
same dependence between the times at which the probabilistic edges can be
taken as in the case of the previous paragraph, and hence a reduction to IMDPs
(or even standard MDPs) is likely to be challenging if not impossible. Future
work could explore the applications of a subclass of this formalism, for
example with multiple clocks but only one clock used for clock dependencies,
or a more strict definition of initialisation than that described above: such
subclasses may yield models that have more significance in practice than
1c-cdPTAs while offering the possibility of exact rather than approximate
analysis.

Finally, we also discuss lifting the restriction to affine clock dependencies
in the case of initialised 1c-cdPTAs. We recall that the IMDP construction
presented in this paper depends on the fact that the distributions
corresponding to the endpoints of intervals in $\mathcal{I}_{\mathbb{B}}$ are
“extremal”, in the sense that, for a given interval, non-constant
probabilistic edge and outcome, the maximum or minimum probability assigned to
that outcome according to the probabilistic edge’s distribution template will
be assigned at one of the endpoints of the interval. This property allows any
scheduler of the 1c-cdPTA to be mimicked by a IMDP scheduler by using
probabilistic edges available from regions to obtain linear combinations of
the distributions corresponding to an interval’s endpoints. In the case of
non-affine clock dependencies, this property does not necessarily hold, i.e.,
the maximum or minimum probability assigned to an outcome may fall _within_ an
interval, rather than at one of its endpoints. Hence, dealing with non-affine
clock dependencies requires non-trivial developments beyond those featured in
this paper, or requires extra assumptions in order to handle the greater
expressiveness of the clock dependencies, for example restrictions on the non-
affine functions used (for example, so that their maximum or minimum is always
at the endpoint of any interval in $\mathcal{I}_{\mathbb{B}}$) or, for any
probabilistic edge, on the number of outcomes that can have non-affine clock
dependencies.

## References

  * [ABK+16] S. Akshay, P. Bouyer, S. N. Krishna, L. Manasa, and A. Trivedi.  Stochastic timed games revisited.  In P. Faliszewski, A. Muscholl, and R. Niedermeier, editors, Proc. MFCS 2016, volume 58 of LIPIcs, pages 8:1–8:14. Leibniz-Zentrum für Informatik, 2016.
  * [AD94] R. Alur and D. L. Dill.  A theory of timed automata.  Theoretical Computer Science, 126(2):183–235, 1994.
  * [BBB+14] N. Bertrand, P. Bouyer, T. Brihaye, Q. Menet, C. Baier, M. Größer, and M. Jurdzinski.  Stochastic timed automata.  Logical Methods in Computer Science, 10(4), 2014.
  * [BBG14] N. Bertrand, T. Brihaye, and B. Genest.  Deciding the value 1 problem for reachability in 1-clock decision stochastic timed automata.  In G. Norman and W. H. Sanders, editors, Proc. QEST 2004, volume 8657 of LNCS, pages 313–328. Springer, 2014.
  * [BBH+19] N. Bertrand, B. Bordais, L. Hélouët, T. Mari, J. Parreaux, and O. Sankur.  Performance evaluation of metro regulations using probabilistic model-checking.  In S. Collart Dutilleul, T. Lecomte, and A. B. Romanovsky, editors, Proc. RSSRail 2019, volume 11495 of LNCS, pages 59–76. Springer, 2019\.
  * [BFL+18] P. Bouyer, U. Fahrenberg, K. G. Larsen, N. Markey, J. Ouaknine, and J. Worrell.  Model checking real-time systems.  In E. M. Clarke, T. A. Henzinger, H. Veith, and R. Bloem, editors, Handbook of Model Checking, pages 1001–1046. Springer, 2018.
  * [BK08] C. Baier and J.-P. Katoen.  Principles of model checking.  MIT Press, 2008.
  * [BLM08] P. Bouyer, K. G. Larsen, and N. Markey.  Model checking one-clock priced timed automata.  Logical Methods in Computer Science, 4(2), 2008.
  * [CDH10] K. Chatterjee, L. Doyen, and T. A. Henzinger.  Qualitative analysis of partially-observable Markov decision processes.  In Petr Hlinený and Antonín Kučera, editors, Proc. MFCS 2010, volume 6281 of LNCS, pages 258–269. Springer, 2010.
  * [CGP01] E. M. Clarke, O. Grumberg, and D. A. Peled.  Model checking.  MIT Press, 2001.
  * [CHK13] T. Chen, T. Han, and M. Kwiatkowska.  On the complexity of model checking interval-valued discrete time Markov chains.  Information Processing Letters, 113(7):210–216, 2013.
  * [CK15] S. Chakraborty and J.-P. Katoen.  Model checking of open interval Markov chains.  In M. Gribaudo, D. Manini, and A. Remke, editors, Proc. ASMTA 2015, volume 9081 of LNCS, pages 30–42. Springer, 2015.
  * [CSH08] K. Chatterjee, K. Sen, and T. A. Henzinger.  Model-checking $\omega$-regular properties of interval Markov chains.  In R. Amadio, editor, Proc. FOSSACS 2008, volume 4962 of LNCS, pages 302–317. Springer, 2008.
  * [FKNP11] V. Forejt, M. Kwiatkowska, G. Norman, and D. Parker.  Automated verification techniques for probabilistic systems.  In M. Bernardo and V. Issarny, editors, Formal Methods for Eternal Networked Software Systems (SFM 2011), volume 6659 of LNCS, pages 53–113. Springer, 2011.
  * [FWHT16] L. Feng, C. Wiltsche, L. R. Humphrey, and U. Topcu.  Synthesis of human-in-the-loop control protocols for autonomous systems.  IEEE Trans. Automation Science and Engineering, 13(2):450–462, 2016.
  * [GJ95] H. Gregersen and H. E. Jensen.  Formal design of reliable real time systems.  Master’s thesis, Department of Mathematics and Computer Science, Aalborg University, 1995.
  * [GLD00] R. Givan, S. M. Leach, and T. L. Dean.  Bounded-parameter Markov decision processes.  Artificial Intelligence, 122(1-2):71–109, 2000.
  * [Hah13] E. M. Hahn.  Model checking stochastic hybrid systems.  PhD thesis, Universität des Saarlandes, 2013.
  * [HHK14] V. Hashemi, H. Hatefi, and J. Krcál.  Probabilistic bisimulations for PCTL model checking of interval MDPs.  In É. André and G. Frehse, editors, Proc. SynCoP 2014, volume 145 of EPTCS, pages 19–33, 2014.
  * [HHZ11] E. M. Hahn, T. Han, and L. Zhang.  Synthesis for PCTL in parametric Markov decision processes.  In M. G. Bobaru, K. Havelund, G. J. Holzmann, and R. Joshi, editors, Proc. NFM 2011, volume 6617 of LNCS, pages 146–161. Springer, 2011\.
  * [HKPV98] T. A. Henzinger, P. W. Kopke, A. Puri, and P. Varaiya.  What’s decidable about hybrid automata?  Journal of Computer and System Sciences, 57(1):94–124, 1998.
  * [HM18] S. Haddad and B. Monmege.  Interval iteration algorithm for MDPs and IMDPs.  Theoretical Computer Science, 735:111–131, 2018.
  * [JL91] B. Jonsson and K. G. Larsen.  Specification and refinement of probabilistic processes.  In Proc. LICS 1991, pages 266–277. IEEE Computer Society, 1991\.
  * [JLS08] M. Jurdziński, F. Laroussinie, and J. Sproston.  Model checking probabilistic timed automata with one or two clocks.  Logical Methods in Computer Science, 4(3):1–28, 2008.
  * [KNSS02] M. Kwiatkowska, G. Norman, R. Segala, and J. Sproston.  Automatic verification of real-time systems with discrete probability distributions.  Theoretical Computer Science, 286:101–150, 2002.
  * [KU02] I. O. Kozine and L. V. Utkin.  Interval-valued finite Markov chains.  Reliable Computing, 8(2):97–113, 2002.
  * [LMS04] F. Laroussinie, N. Markey, and P. Schnoebelen.  Model checking timed automata with one or two clocks.  In P. Gardner and N. Yoshida, editors, Proc. CONCUR 2004, volume 3170 of LNCS, pages 387–401. Springer, 2004.
  * [NE05] A. Nilim and L. El Ghaoui.  Robust control of Markov decision processes with uncertain transition matrices.  Operations Research, 53(5):780–798, 2005.
  * [NPS13] G. Norman, D. Parker, and J. Sproston.  Model checking for probabilistic timed automata.  Formal Methods in System Design, 43(2):164–190, 2013.
  * [PLSS13] A. Puggelli, W. Li, A. L. Sangiovanni-Vincentelli, and S. A. Seshia.  Polynomial-time verification of PCTL properties of MDPs with convex uncertainties.  In N. Sharygina and H. Veith, editors, Proc. CAV 2013, volume 8044 of LNCS, pages 527–542. Springer, 2013.
  * [PT87] C. H. Papadimitriou and J. N. Tsitsiklis.  The complexity of Markov decision processes.  Mathematics of Operations Research, 12(3):441–450, 1987.
  * [Spr18] J. Sproston.  Qualitative reachability for open interval Markov chains.  In I. Potapov and P.-A. Reynier, editors, Proc. RP 2018, volume 11123 of LNCS, pages 146–160. Springer, 2018.
  * [Spr20] J. Sproston.  Probabilistic timed automata with one clock and initialised clock-dependent probabilities.  In A. Gotsman and A. Sokolova, editors, Proc. FORTE 2020, volume 12136 of LNCS, pages 150–168. Springer, 2020.
  * [Spr21] J. Sproston.  Probabilistic timed automata with clock-dependent probabilities.  Fundamenta Informaticae, 178(1–2):101–138, 2021.
  * [SVA06] K. Sen, M. Viswanathan, and G. Agha.  Model-checking Markov chains in the presence of uncertainties.  In H. Hermanns and J. Palsberg, editors, Proc. TACAS 2006, volume 3920 of LNCS, pages 394–410, 2006.
  * [WJPK21] T. Winkler, S. Junges, G. A. Pérez, and J.-P. Katoen.  On the complexity of reachability in parametric Markov decision processes.  Journal of Computer and System Sciences, 119:183–210, 2021.
  * [WK08] D. Wu and X. D. Koutsoukos.  Reachability analysis of uncertain systems using bounded-parameter Markov decision processes.  Artificial Intelligence, 172(8-9):945–954, 2008.


# Time-Reversal of Stochastic Maximum Principle

Amirhossein Taghvaei A. Taghvaei is with the William E. Boeing Department of
Aeronautics and Astronautics at University of Washington, Seattle
<EMAIL_ADDRESS>
###### Abstract

Stochastic maximum principle (SMP) specifies a necessary condition for the
solution of a stochastic optimal control problem. The condition involves a
coupled system of forward and backward stochastic differential equations
(FBSDE) for the state and the adjoint processes. Numerical solution of the
FBSDE is challenging because the boundary condition of the adjoint process is
specified at the terminal time, while the solution should be adaptable to the
forward in time filtration of a Wiener process. In this paper, a “time-
reversal” of the FBSDE system is proposed that involves integration with
respect to a backward in time Wiener process. The time-reversal is used to
propose an iterative Monte-Carlo procedure to solves the FBSDE system and its
time-reversal simultaneously. The procedure involves approximating the
Föllmer’s drift and solving a regression problem between the state and its
adjoint at each time. The procedure is illustrated for the linear quadratic
(LQ) optimal control problem with a numerical example.

##  I Introduction

This paper is concerned with the continuous-time stochastic optimal control
problems that involve the controlled stochastic differential equation (SDE)

$\mathrm{d}X_{t}=a(X_{t},U_{t})\mathrm{d}t+\sigma(X_{t})\mathrm{d}W_{t},\quad
X_{0}\sim p_{0},$ (1)

where $X:=(X_{t})_{t\geq 0}$ is the $n$-dimensional state trajectory,
$U:=(U_{t})_{t\geq 0}$ is the $m$-dimensional control input,
$W:=(W_{t})_{t\geq 0}$ is the standard $n$-dimensional Wiener process, $p_{0}$
is a probability distribution on $\mathbb{R}^{n}$, and
$a:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n}$ and
$\sigma:\mathbb{R}^{n}\to\mathbb{R}^{n\times n}$ are smooth and globally
Lipschitz functions111In this article, the model is restricted to the case
that the diffusion term is not affected by the control input.. The control
input $U_{t}$ is assumed to be adapted to the filtration
$\mathcal{F}_{t}=\sigma\\{W_{s};0\leq s\leq t\\}$, generated by the Wiener
process. The objective of the optimal control problem is

$\min_{U_{t}\in\mathcal{F}_{t}}\,J(U):=\mathbb{E}\left[\int_{0}^{T}\ell(X_{t},U_{t})\mathrm{d}t+\ell_{f}(X_{T})\right],$
(2)

where $\ell:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}$ and
$\ell_{f}:\mathbb{R}^{n}\to\mathbb{R}$ are smooth functions, and the notation
$U_{t}\in\mathcal{F}_{t}$ means $U_{t}$ is $\mathcal{F}_{t}$-measurable.

The stochastic maximum principle (SMP) provides a necessary condition for the
optimal control input. In order to express SMP, define the Hamiltonian
function

$H(x,u,y,z):=\ell(x,u)+y^{\top}a(x,u)+\text{Tr}(z^{\top}\sigma(x)),$ (3a) for
$(x,u,y,z)\in\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{n}\times\mathbb{R}^{n\times
n}$ and the Hamiltonian system of forward-backward SDE (FBSDE)
$\displaystyle\mathrm{d}X_{t}$ $\displaystyle=\frac{\partial H}{\partial
y}(X_{t},U_{t},Y_{t},Z_{t})\mathrm{d}t+\sigma(X_{t})\mathrm{d}{W}_{t},~{}X_{0}\sim
p_{0},$ (3b) $\displaystyle-\mathrm{d}Y_{t}$ $\displaystyle=\frac{\partial
H}{\partial
x}(X_{t},U_{t},Y_{t},Z_{t})\mathrm{d}t-Z_{t}\mathrm{d}{W}_{t},~{}Y_{T}=\frac{\partial\ell_{f}}{\partial
x}(X_{T}).$ (3c) The solution to the FBSDE is a $\mathcal{F}_{t}$-adapted
triple $(X,Y,Z)=(X_{t},Y_{t},Z_{t})_{0\leq t\leq T}$. SMP states that if
$\hat{U}$ is the optimal control input for the objective (2) subject to the
constraint (1), then $\displaystyle\hat{U}_{t}\in\mathop{\text{\rm
arg\,min}}_{u}\,H(\hat{X}_{t},u,\hat{Y}_{t},\hat{Z}_{t}),\quad\text{for}\quad
t\in[0,T],$ (3d) where $\hat{X}$, $\hat{Y}$, and $\hat{Z}$ solve the FBSDE
(3b)-(3c) with $U=\hat{U}$ [1, Thm. 3.2, Case 1].

Numerical implementation of the FBSDE is challenging and it has been subject
of research since its introduction. The challenges include (i) obtaining a
$\mathcal{F}_{t}$-adapted solution $Y_{t}$ while satisfying a terminal
condition, (ii) and obtaining the adjoint diffusion coefficient $Z_{t}$. There
are generally two approaches to numerically solve the FBSDE, the PDE-based
approaches that use the fact that $Y_{t}=\phi(t,X_{t})$, where $\phi$ solves a
parabolic pde [2, 3, 4], and direct numerical approxiamtion of the FBSDE that
involve approximating conditional expectations with respect to the filtration
[5, 6, 7, 8, 9].

In this paper, we explore an alternative approach to numerically solve the
FBSDE system that is based on time-reversal of diffusions [10, 11, 12, 13, 14,
15]. The starting point is to formulate the optimal control problem for the
time-reversal of the SDE (1) where the solution $X_{t}$ and the control
$U_{t}$ become adaptable to a backward filtration of a Wiener process (see
(10)). Compared to SDE (1), the time-reversal involves additional drift term,
referred to as the Föllmer’s drift [13], that depends on the probability
density of $X_{t}$. As a result, the reversed SDE is interpreted as a
controlled McKean-Vlasov SDE leading to a mean-field type optimal control
problem [16, 17, 18]. The appropriate SMP for a mean-field optimal control
problem leads to a “time-reversal” of the FBSDE (3b)-(3c), where the solution
is required to be adapted to a backward filtration (see Sec. III-A). The time-
reversal procedure is justified in Prop. 1 where it is shown that the solution
to (3b)-(3c) and its time-reversal admit the same probability law.

The time-reversed FBSDE is numerically useful as it facilitates the numerical
integration of the adjoint process. Based on this, an iterative numerical
procedure is proposed in Sec. V that simultaneously solves the FBSDE and its
time-reversal. Each iteration of this algorithm involves (i) the forward
Monte-Carlo simulations of the state process; (ii) approximation of the
Föllmer’s drift; (iii) backward simulation of the state and adjoint processes
according to the time-reversed FBSDE; and (iv) a gradient descent update of
the control input to minimize the Hamiltonian function. The third step
involves solving a regression problem between the state and adjoint for all
time $t\in[0,T]$. The proposed procedure is presented for the linear quadratic
(LQ) setting, where the approximation of the Föllmer’s drift simplifies to
mean and covariance estimation, and the regression problem becomes linear. A
two-dimensional numerical example is provided to illustrates the performance
of the algorithm, while the extension to the nonlinear setup and extensive
numerical experiment is subject of future work.

Relation to the other approaches: This paper is closely related to the recent
works that solve stochastic optimal control problem by approximating the
solution to the FBSDE system [19, 20, 21, 22, 23]. In particular, similar to
[19], the proposed approach involves solving a regression type problem at each
time. However, the numerical approximation of the adjoint process is different
in comparison to [19, Eq. 18] as it involves an additional term due to the
time-reversal effect.

##  II Preliminaries

This section contains a brief introduction to the technical concepts that are
essential to the proposed procedure.

###  II-A Backward stochastic integration

We follow the notation introduced in [24, Sec. 2]. We use $W_{t}$ to denote a
standard $n$-dimensional Wiener process which introduces the forward and
backward filtrations

$\displaystyle{\mathcal{F}}_{t}:=$ $\displaystyle\sigma(W_{s};0\leq s\leq
t),\quad\cev{\mathcal{F}}_{t}:=\sigma(W_{T}-W_{s};t\leq s\leq T),$

for $t\in[0,T]$ (the backward filtration is increasing as $t\downarrow 0$).
For any continuous ${\cev{\mathcal{F}}}_{t}$-measurable process $\cev{V}_{t}$,
the backward stochastic integral

$\displaystyle\int_{T-t}^{T}\cev{V}_{s}\cev{\mathrm{d}}W_{s}:=-\int_{0}^{t}\cev{V}_{T-s}\mathrm{d}\tilde{W}_{s}$

where the right-hand-side is the usual Itô integral with respect to the Wiener
process $\tilde{W}_{t}:=W_{T-t}-W_{T}$. Numerically, the backward stochastic
integral may be approximated as

$\int_{T-t}^{T}\cev{V}_{s}\cev{\mathrm{d}}W_{s}=\lim_{n\to\infty}\sum_{k=0}^{n-1}\cev{V}_{T-k\Delta
t}(W_{T-k\Delta t}-W_{T-(k+1)\Delta t})$

where $\Delta t=\frac{t}{n}$.

We use the notation $\cev{\mathrm{d}}W_{t}$ to express stochastic differential
equations (SDEs) that involve backward stochastic integration. For example, a
SDE of the form

$\mathrm{d}\cev{X}_{t}=a(\cev{X}_{t})\mathrm{d}t+\sigma(\cev{X}_{t})\cev{\mathrm{d}}W_{t},\quad\cev{X}_{T}=\xi,$
(4)

is a short-hand notation for

$\displaystyle\cev{X}_{T-t}$
$\displaystyle=\xi-\int_{T-t}^{T}a(\cev{X}_{s})\mathrm{d}s-\int_{T-t}^{T}\sigma(\cev{X}_{s})\cev{\mathrm{d}}W_{s},~{}\forall
t\in[0,T]$ (5)

which defines a ${\cev{\mathcal{F}}}_{t}$-adapted solution $\cev{X}_{t}$.

###### Lemma 1

The SDE (4) is equivalent to the SDE

$\mathrm{d}\tilde{X}_{t}=-a(\tilde{X}_{t})\mathrm{d}t+\sigma(\tilde{X}_{t})\mathrm{d}\tilde{W}_{t},\quad\tilde{X}_{0}=\xi,$
(6)

where $\tilde{W}_{t}:=W_{T-t}-W_{T}$ and $\tilde{X}_{t}:=\cev{X}_{T-t}$ is
adapted to the forward filtration
$\tilde{\mathcal{F}}_{t}:=\sigma\\{\tilde{W}_{s};0\leq s\leq t\\}$, for
$t\in[0,T]$.

###### Proof:

Using $\tilde{X}_{t}=\cev{X}_{T-t}$ in (5) yields

$\displaystyle\tilde{X}_{t}$
$\displaystyle=\xi-\int_{T-t}^{T}a(\tilde{X}_{T-s})\mathrm{d}s-\int_{T-t}^{T}\sigma(\cev{X}_{s})\cev{\mathrm{d}}W_{s}$
$\displaystyle=\xi-\int_{0}^{t}a(\tilde{X}_{s})\mathrm{d}s+\int_{0}^{t}\sigma(\tilde{X}_{s})\mathrm{d}\tilde{W}_{s},$

concluding the forward in time SDE (6). ∎

###  II-B Time reversal of SDE

Consider the two stochastic processes $X=\\{X_{t};\,0\leq t\leq T\\}$ and
$\cev{X}=\\{\cev{X}_{t};\,0\leq t\leq T\\}$ that are governed by the following
SDEs

$\displaystyle\mathrm{d}X_{t}$
$\displaystyle=\tilde{a}(t,X_{t})\mathrm{d}t+\sigma(X_{t})\mathrm{d}W_{t},~{}X_{0}\sim
p_{0},$ (7a) $\displaystyle\mathrm{d}\cev{X}_{t}$
$\displaystyle=\tilde{a}(t,\cev{X}_{t})\mathrm{d}t+b(\cev{X}_{t},p_{X_{t}})\mathrm{d}t+\sigma(\cev{X}_{t})\cev{\mathrm{d}}{W}_{t},~{}\cev{X}_{T}\sim
p_{X_{T}},$ (7b) where, $b(x,p):=-\frac{1}{p(x)}\nabla\cdot(D(x)p(x)),$ (7c)
$D(x):=\sigma(x)\sigma(x)^{\top}$, and $p_{X_{t}}$ is the probability density
function of $X_{t}$. The Föllmer’s drift term (7c) is also interpreted as the
score function when $D(x)=I$.

The following result from [15, Thm. 1.12] relates the two processes.

###### Theorem 1

Assume $D(x)\succ 0$ for all $x\in\mathbb{R}^{n}$ and the probability law of
$X$ satisfies a certain finite entropy condition (see [15, Thm. 1.12]). Then
the probability law of $X$ and $\cev{X}$, as stochastic processes, are equal,
i.e. $X\overset{d}{=}\cev{X}$.

###### Remark 1

The standard statement of the time reversal result involves writing the SDE
for the time-reversed process $\tilde{X}_{t}:=X_{T-t}$. Theorem 1 is an
equivalent statement: application of Lemma 1 to (7b) gives the SDE for
$\cev{X}_{T-t}\overset{d}{=}X_{T-t}=\tilde{X}_{t}$.

In this paper, we like to apply the result of Thm. 1 to the case of controlled
SDEs of the form (1). This is possible whenever the control input $U_{t}$
follows a deterministic state feedback control law of the form
$U_{t}=k(t,X_{t})$ for a time-varying function
$k:[0,T]\times\mathbb{R}^{n}\to\mathbb{R}^{m}$. Under this assumption, we can
define $\tilde{a}(t,x):=a(x,k(t,x))$ and apply Thm. 1. This assumption on the
control is not restrictive because the optimal control law of the stochastic
optimal control problem (2) is known to be of this form.

###  II-C Wasserstein derivatives

Let $\mathcal{P}_{2}(\mathbb{R}^{n})$ be the space of probability
distributions with finite-second moments. For any
$\mu\in\mathcal{P}_{2}(\mathbb{R}^{n})$, we consider the inner-product of two
vector-fields on $\mathbb{R}^{n}$ as

$\langle
v_{1},v_{2}\rangle_{\mu}=\int_{\mathbb{R}^{n}}v_{1}(x)^{\top}v_{2}(x)\mathrm{d}\mu(x),$

which induces the norm $\|v\|_{\mu}:=\sqrt{\langle v,v\rangle}_{\mu}$. The
tangent space $\text{Tan}_{\mu}(\mathcal{P}_{2}(\mathbb{R}^{n}))$ is
identified with the closure of the space of vector-fields of gradient form
$\nabla\phi$ with respect to $\|\cdot\|_{\mu}$ norm [25, Lemma 8.4.2]. The
tangent space is equipped with the Riemannian metric
$\langle\nabla\phi_{1},\nabla\phi_{2}\rangle_{\mu}$. The Wasserstein
derivative of a real-valued function
$F:\mathcal{P}_{2}(\mathbb{R}^{n})\to\mathbb{R}$ at $\mu$, denoted by
$\frac{\partial
F}{\partial\mu}(\mu)\in\text{Tan}_{\mu}(\mathcal{P}_{2}(\mathbb{R}^{n}))$,
satisfies the identity

$\displaystyle\langle\nabla
F(\mu),\nabla\phi\rangle=\frac{\mathrm{d}}{\mathrm{d}t}F((\text{Id}+t\nabla\phi)_{\\#}\mu)|_{t=0},$

for all $\nabla\phi\in\text{Tan}_{\mu}(\mathcal{P}_{2}(\mathbb{R}^{n}))$,
where Id denotes the identity map and $\\#$ denotes the push-forward operator.
See [18, Def. 5.62] for a rigorous definition of the Wasserstein derivative.
In this article, we use the following property of the Wasserstein derivatives,
shown in [18, Example 3 on page 387]: If
$F(\mu)=\int_{\mathbb{R}^{n}}f(x^{\prime},\mu)\mathrm{d}\mu(x^{\prime})$ and
$f$ is differentiable, then

$\displaystyle\frac{\partial F}{\partial\mu}(\mu)(x)=\frac{\partial
f}{\partial x}(x,\mu)+\int\frac{\partial
f}{\partial\mu}(x^{\prime},\mu)(x)\mathrm{d}\mu(x^{\prime}).$ (8)

###  II-D PDE solution of FBSDE

Consider the FBSDE (3b)-(3c) assuming the control
$U_{t}=k(t,X_{t},Y_{t},Z_{t})$ for some smooth function
$k:[0,T]\times\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{n\times
n}\to\mathbb{R}^{m}$. The PDE solution of the FBSDE is based on the hypothesis
that $Y_{t}=\phi(t,X_{t})$ for some smooth function
$\phi:[0,T]\times\mathbb{R}^{n}\to\mathbb{R}^{n}$. Then, by Itô’s formula

$\displaystyle\mathrm{d}Y_{t}=$ $\displaystyle\frac{\partial\phi}{\partial
t}(t,X_{t})\mathrm{d}t+\frac{\partial\phi}{\partial x}(t,X_{t})\frac{\partial
H}{\partial y}(X_{t},k(t,X_{t},Y_{t},Z_{t}),Y_{t},Z_{t})\mathrm{d}t$
$\displaystyle+\frac{\partial\phi}{\partial
x}(t,X_{t})\sigma(X_{t})\mathrm{d}W_{t}+\frac{1}{2}\text{Tr}(D(X_{t})\frac{\partial^{2}\phi}{\partial
x^{2}}(t,X_{t}))\mathrm{d}t,$

where $\frac{\partial\phi}{\partial x}$ and $\frac{\partial^{2}\phi}{\partial
x^{2}}$ are the Jacobian and Hessian respectively. Comparing to (3c) conlcudes
that $Z_{t}=\frac{\partial\phi}{\partial x}(t,X_{t})\sigma(X_{t})$ and $\phi$
satisfies the PDE

$\displaystyle\frac{\partial\phi}{\partial t}+\frac{\partial\phi}{\partial
x}\frac{\partial H}{\partial y}(t,x,k,\phi,\frac{\partial\phi}{\partial
x}\sigma)+\frac{1}{2}\text{Tr}(D\frac{\partial^{2}\phi}{\partial x^{2}})$
$\displaystyle+\frac{\partial H}{\partial
x}(t,x,k,\phi,\frac{\partial\phi}{\partial
x}\sigma)=0,~{}\phi(T,x)=\frac{\partial\ell_{f}}{\partial x}(x),$ (9)

for all $(x,t)\in\mathbb{R}^{n}\times[0,T]$, where, for convenience, $k$ and
$\sigma$ are used to denote $k(t,x,y,z)$ and $\sigma(x)$. The existence of the
solution to this parabolic PDE is investigated in [3] under certain regularity
assumptions. In this paper, we make the assumption that the solution to the
PDE (9) exists.

##  III Time-reversal procedure

In this section, we introduce the “time-reversal” formulation of the
stochastic control problem (1)-(2). Based on the discussion presented in Sec.
II-B, we consider the time-reversal of (1) according to

$\displaystyle\mathrm{d}\cev{X}_{t}=$ $\displaystyle
a(\cev{X}_{t},\cev{U}_{t})\mathrm{d}t+b(\cev{X}_{t},p_{\cev{X}_{t}})\mathrm{d}t+\sigma(\cev{X}_{t})\cev{\mathrm{d}}{W}_{t},\quad\cev{X}_{T}\sim
p_{T},$ (10)

where $\cev{\mathrm{d}}W_{t}$ denotes a backward stochastic integration,
$b(x,p)$ is defined in (7c), and $p_{\cev{X}_{t}}$ is the probability density
function of $\cev{X}_{t}$. We consider control inputs $\cev{U}_{t}$ that are
adapted to the backward filtration $\cev{\mathcal{F}}_{t}$ that lead to an
${\cev{\mathcal{F}}}_{t}$-adapted solution $\cev{X}_{t}$. The time-reversal of
the objective (2) becomes

$\min_{\cev{U}_{t}\in{\cev{\mathcal{F}}}_{t}}\,\mathbb{E}\left[\int_{0}^{t_{f}}\ell(\cev{X}_{t},\cev{U}_{t})\mathrm{d}t+\ell_{f}(\cev{X}_{T})\right].$
(11)

The SDE (10) involves an additional drift term
$b(\cev{X}_{t},p_{\cev{X}_{t}})$ that depends on the probability distribution
of ${\cev{X}_{t}}$. As a result, it belongs to the class of McKean-Vlasov
SDEs, which requires an appropriate extension of the SMP [18, Sec. 6.3].

###  III-A SMP for the time-reversed problem

In order to describe the maximum principle for the McKean-Vlasov SDE (10) with
the objective function (11), define the (backward) Hamiltonian

$\displaystyle\cev{H}(x,p,u,y,z):=\ell(x,u)+y^{\top}(a(x,u)+b(x,p))-\text{Tr}(z^{\top}\sigma(x)).$
(12a) for
$(x,p,u,y,z)\in\mathbb{R}^{n}\times\mathcal{P}_{2}(\mathbb{R}^{n})\times\mathbb{R}^{m}\times\mathbb{R}^{n}\times\mathbb{R}^{n\times
n}$. Note that, in comparison to the Hamiltonian (3a), the sign of $z$ is
changed to account for the backward stochastic integration. The backward
Hamiltonian leads to the FBSDE $\displaystyle\mathrm{d}\cev{X}_{t}$
$\displaystyle=\\!\frac{\partial\cev{H}}{\partial
y}(\cev{X}_{t},p_{\cev{X}_{t}},\cev{U}_{t},\cev{Y}_{t},\cev{Z}_{t})\mathrm{d}t\\!+\\!\sigma(\cev{X}_{t})\mathrm{d}\cev{W}_{t},~{}\cev{X}_{T}\\!\sim\\!p_{T}$
(12b) $\displaystyle-\mathrm{d}\cev{Y}_{t}$
$\displaystyle=\frac{\partial\cev{H}}{\partial
x}(\cev{X}_{t},p_{\cev{X}_{t}},\cev{U}_{t},\cev{Y}_{t},\cev{Z}_{t})\mathrm{d}t-\cev{Z}_{t}\mathrm{d}\cev{W}_{t}$
(12c) $\displaystyle+\tilde{\mathbb{E}}\left[\frac{\partial\cev{H}}{\partial
p}(\tilde{X}_{t},p_{\cev{X}_{t}},\tilde{U}_{t},\tilde{Y}_{t},\tilde{Z}_{t})(\cev{X}_{t})\right]\mathrm{d}t,\quad\cev{Y}_{t_{f}}=\frac{\partial\ell_{f}}{\partial
x}(\cev{X}_{t_{f}}),$ where $\frac{\partial\cev{H}}{\partial
p}(x,p,u,y,z):\mathbb{R}^{n}\to\mathbb{R}^{n}$ is the Wasserstein derivative
of $\cev{H}$ with respect to $p$, and the expectation $\mathbb{\tilde{E}}$ is
with respect to
$\\{\tilde{X}_{t},\tilde{Y}_{t},\tilde{Z}_{t},\tilde{U}_{t}\\}$ which are
independent copies of $\\{\cev{X}_{t},\cev{Y}_{t},\cev{Z}_{t},\cev{U}_{t}\\}$.
The SMP states that if $\hat{U}$ is the optimal control input, then
$\displaystyle\hat{U}_{t}\in\mathop{\text{\rm
arg\,min}}_{u\in\mathbb{R}^{m}}\,\cev{H}(\hat{X}_{t},p_{\hat{X}_{t}},u,\hat{Y}_{t},\hat{Z}_{t})$
(12d) where $\hat{X}$, $\hat{Y}$, and $\hat{Z}$ solve the FBSDE (12b)-(12c)
with $\cev{U}=\hat{U}$ [18, Thm. 6.14].

In the next section, we show that the new FBSDE (12b)-(12c) and the original
FBSDE (3b)-(3c) both admit the same probability law for the state and adjoint
processes.

###  III-B Time-reversal relationship

We present the result under the assumption that the control inputs $U_{t}$ and
$\cev{U}_{t}$ follow the same feedback control law

$U_{t}=k(t,X_{t},Y_{t},Z_{t}),\quad\cev{U}_{t}=k(t,\cev{X}_{t},\cev{Y}_{t},\cev{Z}_{t}),$
(13)

for some smooth function $k$.

###### Proposition 1

Let $(\cev{X},\cev{Y},\cev{Z})$ be the triple that solves the FBSDE
(12c)-(12b) with $\cev{U}_{t}$ as in (13). Assume the PDE (9) has a smooth
solution $\phi$. Then,

$\displaystyle(\cev{Y}_{t},\cev{Z}_{t})=(\phi(t,\cev{X}_{t}),\frac{\partial\phi}{\partial
x}(t,\cev{X}_{t})\sigma(\cev{X}_{t})),\quad\text{for}\quad t\in[0,T].$ (14)

Moreover, assume $(X,Y,Z)$ solves the FBSDE (3b)-(3c) with $U_{t}$ as in (13),
and $p_{T}$ is equal to the probability density function of $X_{T}$. Then, the
probability law of the quadruple $(\cev{X},\cev{U},\cev{Y},\cev{Z})$ is equal
to the probability law of $(X,U,Y,Z)$, i.e.

$(X,U,Y,Z)\overset{d}{=}(\cev{X},\cev{U},\cev{Y},\cev{Z}).$ (15)

###### Remark 2

This result establishes that the new FBSDE (12c)-(12b) is indeed the time-
reversal of the FBSDE (3b)-(3c). The key element that allows this time-
reversal is the restriction of the control $U_{t}$ to feedback control laws of
the form (13). This assumption ensures that the control has a Markov property,
i.e. it depends only on the current value of the state, conditionally
independent of its past and future. As a result, the control is adaptable to
both forward and backward filtrations generated by the state process. This
assumption is appropriate for stochastic optimal control problems since the
optimal control law is of this form.

Before presenting the proof of this proposition, we present the following
Lemma that gives the expression for the Wasserstein derivative of the
Hamiltonian.

###### Lemma 2

If $\cev{Y}_{t}=\phi(t,\cev{X}_{t})$, then

$\displaystyle{\mathbb{E}}$
$\displaystyle\left[\frac{\partial\cev{H}}{\partial
p}(\cev{X}_{t},p_{t},\cev{U}_{t},\cev{Y}_{t},\cev{Z}_{t})(x)\right]\\!=$
$\displaystyle\frac{\partial}{\partial x}(D(x)\frac{\partial\phi}{\partial
x}(t,x)-b(x,p_{t})^{\top}\phi(t,x)).$ (16)

where $p_{t}$ is the probability density function of ${\cev{X}_{t}}$.

###### Proof:

By the definition of the Hamiltonian

$\displaystyle{\mathbb{E}}$
$\displaystyle\left[\frac{\partial\cev{H}}{\partial
p}(\cev{X}_{t},p_{t},\cev{U}_{t},\cev{Y}_{t},\cev{Z}_{t})(x)\right]={\mathbb{E}}\left[\frac{\partial
b}{\partial p}(\cev{X}_{t},p_{t})(x)^{\top}\cev{Y}_{t}\right]$
$\displaystyle=\\!\int\frac{\partial b}{\partial
p}(x^{\prime},p_{t})(x)^{\top}\phi(t,x^{\prime})p_{t}(x^{\prime})\mathrm{d}x^{\prime}.$

Define $\beta(p):=\int
b(x^{\prime},p)^{\top}\phi(t,x^{\prime})p(x^{\prime})\mathrm{d}x^{\prime}$.
Then, the identity (16) implies

$\displaystyle\int\frac{\partial b}{\partial p}(x^{\prime},p_{t})(x)^{\top}$
$\displaystyle\phi(t,x^{\prime})p_{t}(x^{\prime})\mathrm{d}x^{\prime}=$
$\displaystyle\frac{\partial\beta}{\partial p}(x)-\frac{\partial}{\partial
x}(b(x,p_{t})^{\top}\phi(t,x)).$

Direct evaluation of $\frac{\partial\beta}{\partial p}$, using the definition
of $b(x,p)$, gives

$\displaystyle\frac{\partial\beta}{\partial p}$
$\displaystyle=-\frac{\partial}{\partial
p}\int\nabla\cdot(p(x^{\prime})D(x^{\prime}))\phi(t,x^{\prime})\mathrm{d}x^{\prime}$
$\displaystyle=\frac{\partial}{\partial p}\int
D(x^{\prime})\frac{\partial\phi}{\partial
x}(t,x^{\prime})p(x^{\prime})\mathrm{d}x^{\prime}=\frac{\partial}{\partial
x}(D(x)\frac{\partial\phi}{\partial x}(t,x)),$

concluding the claimed identity(16). ∎

Proof of Prop. 1: Application of the Itô rule to
$\cev{Y}_{t}=\phi(t,\cev{X}_{t})$ implies

$\displaystyle\mathrm{d}\cev{Y}_{t}=$
$\displaystyle\frac{\partial\phi}{\partial
t}(t,\cev{X}_{t})\mathrm{d}t+\frac{\partial\phi}{\partial
x}(t,\cev{X}_{t})\frac{\partial\cev{H}}{\partial
y}(\cev{X}_{t},\cev{U}_{t},\cev{Y}_{t},\cev{Z}_{t})\mathrm{d}t$
$\displaystyle+\frac{\partial\phi}{\partial
x}(t,\cev{X}_{t})\sigma(\cev{X}_{t})\cev{\mathrm{d}}W_{t}-\frac{1}{2}\text{Tr}(D(\cev{X}_{t})\frac{\partial^{2}\phi}{\partial
x^{2}}(t,\cev{X}_{t}))\mathrm{d}t,$

where the Itô correction term appears with the negative sign due to the
backward stochastic integration, and we used $p_{t}$ to denote
$p_{\cev{X}_{t}}$. Using $\cev{Z}_{t}=\frac{\partial\phi}{\partial
x}(t,\cev{X}_{t})\sigma(\cev{X}_{t})$, $\frac{\partial\cev{H}}{\partial
y}=\frac{\partial H}{\partial y}+b$, and the fact that $\phi$ solves (9),

$\displaystyle\mathrm{d}\cev{Y}_{t}=$
$\displaystyle\frac{\partial\phi}{\partial
x}(t,\cev{X}_{t})b(\cev{X}_{t},p_{t})\mathrm{d}t-\frac{\partial H}{\partial
x}(\cev{X}_{t},\cev{U}_{t},\cev{Y}_{t},\cev{Z}_{t})\mathrm{d}t$
$\displaystyle+\cev{Z}_{t}\cev{\mathrm{d}}W_{t}-\text{Tr}(D(\cev{X}_{t})\frac{\partial^{2}\phi}{\partial
x^{2}}(t,\cev{X}_{t}))\mathrm{d}t.$

On the other-hand, starting from (12c), using the definition (12a), and the
identity (16)

$\displaystyle\mathrm{d}\cev{Y}_{t}=$ $\displaystyle-\frac{\partial
H}{\partial
x}(\cev{X}_{t},\cev{U}_{t},\cev{Y}_{t},\cev{Z}_{t})\mathrm{d}t-\frac{\partial
b}{\partial x}(\cev{X}_{t},p_{t})\cev{Y}_{t}\mathrm{d}t$
$\displaystyle-2\frac{\partial}{\partial
x}\text{Tr}(\sigma(\cev{X}_{t})^{\top}\cev{Z}_{t})\mathrm{d}t-\frac{\partial}{\partial
x}(D(\cev{X}_{t})\frac{\partial\phi}{\partial x}(t,\cev{X}_{t}))$
$\displaystyle+\frac{\partial}{\partial
x}(b(\cev{X}_{t},p_{t})\phi(t,\cev{X}_{t}))+\cev{Z}_{t}\cev{\mathrm{d}}W_{t}.$

Comparing the two derivations and using the following two identities

$\displaystyle\frac{\partial}{\partial
x}(b(x,p)\phi(t,x))=\frac{\partial\phi}{\partial x}(t,x)b(x,p)+\frac{\partial
b}{\partial x}(x,p)\phi(t,x)$ $\displaystyle\frac{\partial}{\partial
x}(D(x)\frac{\partial\phi}{\partial x}(t,x))=\text{Tr}(\frac{\partial
D}{\partial x}(x)\frac{\partial\phi}{\partial
x}(t,x))+\text{Tr}(D(x)\frac{\partial^{2}\phi}{\partial x^{2}}(t,x)).$

concludes (14). Finally, the equality (15) follows by application of Thm. 1
with $\tilde{a}(t,x)=a(x,k(t,x,\phi(t,x)),\frac{\partial\phi}{\partial
x}(t,x)\sigma(x))$. ∎

##  IV Linear quadratic setting

In this section, we consider the special case where the dynamic model is
linear, the initial distribution is Gaussian, and the control cost is
quadratic, i.e.

$\displaystyle a(x,u)=Ax+Bu,\quad\sigma(x)=\sigma,\quad
p_{0}=\mathcal{N}(m_{0},\Sigma_{0})$
$\displaystyle\ell(x,u)=\frac{1}{2}\left[x^{\top}Qx+u^{\top}Ru\right],\quad\ell_{f}(x)=\frac{1}{2}x^{\top}Q_{f}x$

where $Q,Q_{f}\succeq$ are assumed to be positive semi-definite, $R\succ 0$ is
assumed to be positive definite, and $\mathcal{N}(m,\Sigma)$ denotes a
Gaussian distribution with mean $m$ and covariance $\Sigma$. Under the LQ
setting, the Hamiltonian (3a) becomes

$H(x,u,y,z):=\frac{1}{2}\left[x^{\top}Qx+u^{\top}Ru\right]+y^{\top}(Ax+Bu)+\text{Tr}(\sigma^{\top}z),$

concluding the FBSDE

$\displaystyle\mathrm{d}X_{t}$
$\displaystyle=AX_{t}\mathrm{d}t+BU_{t}\mathrm{d}t+\mathrm{d}W_{t},\quad
X_{0}\sim p_{0}$ (17a) $\displaystyle-\mathrm{d}Y_{t}$
$\displaystyle=A^{\top}Y_{t}\mathrm{d}t+QX_{t}\mathrm{d}t-Z_{t}\mathrm{d}W_{t},\quad
Y_{t_{f}}=Q_{f}X_{t_{f}}.$ (17b) Moreover, the control input that minimizes
the Hamiltonian, according to (3d), is equal to $U_{t}=-R^{-1}B^{\top}Y_{t}$.

We present the solution of the FBSDE system above for an affine form of the
control according to

$\displaystyle U_{t}=K_{1}(t)X_{t}+K_{2}(t)Y_{t}+K_{3}(t)Z_{t}+K_{4}(t)$ (17c)

for some time-varying parameters $K_{1}$, $K_{2}$,$K_{3}$, and $K_{4}$ with
appropriate dimensions.

###### Proposition 2

Consider the FBSDE (17a)-(17b) with affine form of control as in (17c). Then,
the FBSDE is solved with $Y_{t}=G_{1}(t)X_{t}+G_{2}(t)$ and
$Z_{t}=G_{1}(t)\sigma$ where $G_{1}$ and $G_{2}$ solve matrix differential
equations

$\displaystyle\dot{G}_{1}(t)+G_{1}(t)A+A^{\top}G_{1}(t)+G_{1}(t)BK_{1}(t)$
$\displaystyle+G_{1}(t)BK_{2}(t)G_{1}(t)+Q=0,\quad G_{1}(T)=Q_{f},$ (18a)
$\displaystyle\dot{G}_{2}(t)+G_{1}(t)B(K_{2}(t)G_{2}(t)+K_{3}(t)G_{1}(t)\sigma+K_{4}(t))$
$\displaystyle+A^{\top}G_{2}(t)=0,\quad G_{2}(T)=0.$ (18b) For the case with
the optimal control law, i.e. $K_{2}(t)=-R^{-1}B^{\top}$ and
$K_{1}(t)=K_{3}(t)=K_{4}(t)=0$, we have $G_{2}(t)=0$ and $G_{1}(t)$ solving
the Riccati equation $\displaystyle\dot{G}_{1}(t)$
$\displaystyle+G_{1}(t)A+A^{\top}G_{1}(t)+Q$ $\displaystyle-
G_{1}(t)BR^{-1}B^{\top}G_{1}(t),\quad G_{1}(T)=Q_{f}.$ (18c)

###### Proof:

The proof follows upon application of the Itô rule to
$Y(t)=G_{1}(t)X_{t}+G_{2}(t)$ and using the formula for $\dot{G}_{1}(t)$ and
$\dot{G}_{2}(t)$. ∎

###  IV-A Time reversal in LQ setting

We like to express the time-reversed FBSDE (12b)-(12c) in the LQ setting. We
assume the control input $\cev{U}_{t}$ has the affine form (17c). Then,
according to Prop. 1, the probability law of
$(\cev{X}_{t},\cev{Y}_{t},\cev{X}_{t},\cev{U}_{t})$ is equal to
$(X_{t},Y_{t},X_{t},U_{t})$. As a result, the probability law of
$\cev{X}_{t}$, denoted by $p_{t}$, is Gaussian, similar to $X_{t}$, and the
drift term $b(x,p_{t})=D\Sigma_{t}^{-1}(x-m_{t})$, where $m_{t}$ and
$\Sigma_{t}$ denote the mean and covariance of $p_{t}$, respectively. This
conlcudes the following SDE for $\cev{X}_{t}$

$\displaystyle\mathrm{d}\cev{X}_{t}=A\cev{X}_{t}\mathrm{d}t+B\cev{U}_{t}\mathrm{d}t+\sigma\cev{\mathrm{d}}{W}_{t}+$
$\displaystyle D\Sigma_{t}^{-1}(\cev{X}_{t}-m_{t})\mathrm{d}t,$ (19a) where
$D:=\sigma\sigma^{\top}$. In order express the adjoint process, we evaluate
the Wasserstein derivative of the Hamiltonian according to Lemma 2, using
$\cev{Y}_{t}=G_{1}(t)\cev{X}_{t}+G_{2}(t)$, $\displaystyle{\mathbb{E}}$
$\displaystyle\left[\frac{\partial\cev{H}}{\partial
p}(\cev{X}_{t},p_{t},\cev{U}_{t},\cev{Y}_{t},\cev{Z}_{t})(x)\right]\\!=-D\Sigma_{t}^{-1}(G_{1}(t)x+G_{2}(t))$
$\displaystyle\quad-G_{1}(t)D\Sigma_{t}^{-1}(x-m_{t}),$ concluding the
following expression for the adjoint process:
$\displaystyle-\mathrm{d}\cev{Y}_{t}$
$\displaystyle=A^{\top}\cev{Y}_{t}\mathrm{d}t+Q\cev{X}_{t}\mathrm{d}t-\cev{Z}_{t}\mathrm{d}\cev{W}_{t}$
(19b) $\displaystyle-
G_{1}(t)D\Sigma_{t}^{-1}(\cev{X}_{t}-m_{t})\mathrm{d}t,\quad\cev{Y}_{t_{f}}=Q_{f}\cev{X}_{t_{f}}.$
Application of the Itô rule verifies that (19b) is indeed solved with
$\cev{Y}_{t}=G_{1}(t)\cev{X}_{t}+G_{2}(t)$ and $\cev{Z}_{t}=G_{1}(t)\sigma$
where $G_{1}$ and $G_{2}$ solve (18a)-(18b). In the proposed numerical scheme,
we do not solve (18a), but compute the matrix $G_{1}(t)$ by solving the
regression problem $G_{1}(t)=\mathop{\text{\rm
arg\,min}}_{G}\,\mathbb{E}\left[\|\cev{Y}_{t}-G\cev{X}_{t}\|^{2}\right]=\text{Cov}(\cev{X}_{t},\cev{Y}_{t})\text{Cov}(\cev{Y}_{t},\cev{Y}_{t})^{-1}$
where $\text{Cov}(V_{1},V_{2})$ denotes the covariance matrix of two random
variables $V_{1}$ and $V_{2}$.

(a) (b) (c)

Figure 1: Numerical results for a two-dimensional linear quadratic stochastic
optimal control problem presented in Sec. V-A: (a) Sampled trajectories of the
state (20a), its time-reversal (20b), and the adjoint process (20c); (b)
numerical approxiamtion of the matrix $G(t)$ according to (20d) in comparison
to its exact value according to (18c); (c) value of the objective function (2)
as the number of algorithm iterations increases.

##  V Numerical Algorithm

A numerical procedure is proposed to approximate the solution to the FBSDE
(17a)-(17b), its time-reversal (19a)-(19b), and the optimal control input.
While the procedure is presented for the LQ setting, it is generalizable to
the nonlinear case as explained in Remark 3.

Initialization: The algorithm begins with generating $N$ independent samples
$\\{X^{i}_{0}\\}_{i=1}^{N}$ from the initial distribution
$p_{0}=\mathcal{N}(m_{0},\Sigma_{0})$. Each sample $X^{i}_{0}$, for
$i=1,\ldots,N$, is then paired with a zero initialization of the control
inputs $U^{i,0}_{t}=\cev{U}^{i,0}_{t}=0$, and a random realization of the
Weiner processes $W^{i}_{t}$ and $\tilde{W}^{i}_{t}$, that will be used for
the simulation of the FBSDE and its time-reversal, respectively.

Iterations: The $k$-th iteration, for $k=1,2,3,\ldots,k_{f}$, involves the
following steps:

  1. 1.

Forward simulations of (17a) according to

$\displaystyle\mathrm{d}X^{i,k}_{t}\\!=\\!AX^{i,k}_{t}\mathrm{d}t+BU^{i,k-1}_{t}\mathrm{d}t+\sigma\mathrm{d}W^{i}_{t},~{}X^{i,k}_{0}=X^{i}_{0}$
(20a)

for $i=1,\ldots,N$ and $t\in[0,T]$.

  2. 2.

Approximation of the Fölmer’s drift. For the LQ setting, this amounts to
approximating the mean and covariance according to

$\displaystyle
m^{k}_{t}\approx\frac{1}{N}\sum_{i=1}^{N}X^{i,k}_{t},~{}\Sigma^{k}_{t}\approx\frac{1}{N}\sum_{i=1}^{N}(X^{i,k}_{t}-m^{k}_{t})(X^{i,k}_{t}-m^{k}_{t})^{\top}.$

  3. 3.

Simulation of time-reversed FBSDE (19a)-(19b) according to

$\displaystyle\mathrm{d}\cev{X}^{k,i}_{t}=A\cev{X}^{k,i}_{t}\mathrm{d}t+B\cev{U}^{i,k-1}_{t}\mathrm{d}t+\sigma\cev{\mathrm{d}}{\tilde{W}}^{i}_{t}$
(20b)
$\displaystyle+D(\Sigma^{k}_{t})^{-1}(\cev{X}^{i,k}_{t}-m^{k}_{t})\mathrm{d}t,~{}X^{i,k}_{T}\sim\mathcal{N}(m^{k}_{T},\Sigma^{k}_{T}),$
$\displaystyle-$
$\displaystyle\mathrm{d}\cev{Y}^{i,k}_{t}=A^{\top}\cev{Y}^{i,k}_{t}\mathrm{d}t+Q\cev{X}^{i,k}_{t}\mathrm{d}t-G_{1}^{k}(t)\sigma\cev{\mathrm{d}}\tilde{W}^{i}_{t}$
(20c) $\displaystyle-
G_{1}^{k}(t)D(\Sigma^{k}_{t})^{-1}(\cev{X}^{i,k}_{t}-m^{k}_{t})\mathrm{d}t,\quad
Y^{i,k}_{T}=Q_{f}X^{i,k}_{T},$

for $i=1,\ldots,N$ and $t\in[0,T]$, where

$\displaystyle G_{1}^{k}(t)$ $\displaystyle=\mathop{\text{\rm
arg\,min}}_{G}\,\frac{1}{N}\sum_{i=1}^{N}\|\cev{Y}^{i,k}_{t}-G\cev{X}^{i,k}_{t}\|^{2}.$
(20d)

  4. 4.

Gradient descent update of the control to minimize the Hamiltonian according
to

$\displaystyle U^{i,k}_{t}$ $\displaystyle=U^{i,k-1}_{t}-\eta\frac{\partial
H}{\partial u}(X^{i,k}_{t},U^{i,k-1}_{t},Y^{i,k}_{t},Z^{k}_{t}),$
$\displaystyle\cev{U}^{i,k}_{t}$
$\displaystyle=\cev{U}^{i,k-1}_{t}-\eta\frac{\partial H}{\partial
u}(\cev{X}^{i,k}_{t},\cev{U}^{i,k-1}_{t},\cev{Y}^{i,k}_{t},\cev{Z}^{k}_{t}),$

for a small step-size $\eta>0$, where
$Z^{k}_{t}=\cev{Z}^{k}_{t}=G_{1}^{k}(t)\sigma$,
$Y^{i,k}_{t}=G_{1}^{k}(t)X^{i,k}_{t}$, and $\frac{\partial H}{\partial
u}(x,u,y,z)=Ru+B^{\top}y$.

Output: At the end of the iterations, the samples
$(X^{i,k_{f}}_{t},U^{i,k_{f}}_{t},G^{k_{f}}_{1}(t)X^{i,k_{f}}_{t},G^{k_{f}}_{1}(t)\sigma)_{i=1}^{N}$
are expected to approximate the probability law of the quadruple
$(\hat{X},\hat{U},\hat{Y},\hat{Z})$ that satisfy the SMP conditions
(3b)-(3c)-(3d) in the LQ setting. A feedback control law is obtained according
to $U_{t}=\mathop{\text{\rm
arg\,min}}_{u}H(X_{t},u,Y_{t},Z_{t})=-R^{-1}B^{\top}Y_{t}\approx-R^{-1}B^{\top}G^{k_{f}}_{1}(t)X_{t}$,
given a realization of the state process $X_{t}$.

###### Remark 3

In order to extend the procedure to the nonlinear setting, it is necessary to
use a score function approximation algorithm in step 2, and solve a nonlinear
regression problem
$\min_{\phi}\,\frac{1}{N}\sum_{i=1}^{N}\|\cev{Y}^{i,k}_{t}-\phi(\cev{X}^{i,k}_{t})\|^{2}$
to find $\phi^{k}(t,x)$ for all $t\in[0,T]$ in step 3. When $\phi^{k}(t,x)$ is
available, we can evaluate the Wasserstein derivative according to Lemma 2 and
use $\cev{Z}_{t}=\frac{\partial\phi}{\partial
x}(t,\cev{X}_{t})\sigma(\cev{X}_{t})$ to find $\cev{Z}_{t}$.

###  V-A Numerical Example

The algorithm is illustrated for a two-dimensional mass-spring example with
the following model parameters

$\displaystyle A=\begin{bmatrix}0&1\\\
-1&0\end{bmatrix},~{}B=\begin{bmatrix}0\\\
1\end{bmatrix},~{}\sigma=\begin{bmatrix}1&0\\\ 0&1\end{bmatrix},$
$\displaystyle R=1,~{}m_{0}=\begin{bmatrix}0\\\ 0\end{bmatrix},\quad
Q=Q_{f}=\Sigma_{0}=\begin{bmatrix}1&0\\\ 0&1\end{bmatrix}.$

The SDEs are simulated using the Euler-Maruyama method with $\Delta t=0.02$.
The number of samples $N=1000$, iterations $k_{f}=75$, and optimization step-
size $\eta=0.02$. The numerical results are shown in Fig. 1. Panel (a)
presents the first component of the trajectories $X^{i,k_{f}}_{t}$,
$\cev{X}^{i,k_{f}}_{t}$, and $\cev{Y}^{i,k_{f}}_{t}$ for $20$ randomly
selected samples. Panel (b) shows the components of the matrix $G^{k_{f}}(t)$,
averaged over $10$ independent simulations, in comparison with their exact
value, derived by solving the Riccati equation (18c). Finally, panel (c) shows
the convergence of the value of the objective function (2) as the number of
iterations increases.

##  VI Conclusion

This paper presents the time-reversal of the FBSDE that appears in stochastic
maximum principle. The time-reversal is used to propose an iterative numerical
algorithm to approximately solve the FBSDE and obtain the optimal control
input. The algorithm is presented and numerically illustrated for the LQ
setting. The extension to the nonlinear setting and subsequent comparison to
the related algorithms is the subject of future research.

## References

  * [1] J. Yong and X. Y. Zhou, _Stochastic controls: Hamiltonian systems and HJB equations_. Springer Science & Business Media, 1999, vol. 43.
  * [2] S. Peng, “Probabilistic interpretation for systems of quasilinear parabolic partial differential equations,” _Stochastics and stochastics reports (Print)_ , vol. 37, no. 1-2, pp. 61–74, 1991.
  * [3] J. Ma, P. Protter, and J. Yong, “Solving forward-backward stochastic differential equations explicitly—a four step scheme,” _Probability theory and related fields_ , vol. 98, no. 3, pp. 339–359, 1994.
  * [4] J. Ma, P. Protter, J. San Martín, and S. Torres, “Numerical method for backward stochastic differential equations,” _Annals of Applied Probability_ , pp. 302–316, 2002.
  * [5] J. Zhang, “A numerical scheme for bsdes,” _The annals of applied probability_ , vol. 14, no. 1, pp. 459–488, 2004.
  * [6] B. Bouchard and N. Touzi, “Discrete-time approximation and monte-carlo simulation of backward stochastic differential equations,” _Stochastic Processes and their applications_ , vol. 111, no. 2, pp. 175–206, 2004.
  * [7] W. Zhao, L. Chen, and S. Peng, “A new kind of accurate numerical method for backward stochastic differential equations,” _SIAM Journal on Scientific Computing_ , vol. 28, no. 4, pp. 1563–1581, 2006.
  * [8] C. Bender and R. Denk, “A forward scheme for backward sdes,” _Stochastic processes and their applications_ , vol. 117, no. 12, pp. 1793–1812, 2007.
  * [9] G. Zhang, M. Gunzburger, and W. Zhao, “A sparse-grid method for multi-dimensional backward stochastic differential equations,” _Journal of Computational Mathematics_ , pp. 221–248, 2013.
  * [10] B. D. Anderson, “Reverse-time diffusion equation models,” _Stochastic Processes and their Applications_ , vol. 12, no. 3, pp. 313–326, 1982.
  * [11] U. G. Haussmann and E. Pardoux, “Time reversal of diffusions,” _The Annals of Probability_ , pp. 1188–1205, 1986.
  * [12] A. Millet, D. Nualart, and M. Sanz, “Integration by parts and time reversal for diffusion processes,” _The Annals of Probability_ , pp. 208–238, 1989\.
  * [13] H. Föllmer, “An entropy approach to the time reversal of diffusion processes,” in _Stochastic Differential Systems Filtering and Control: Proceedings of the IFIP-WG 7/1 Working Conference Marseille-Luminy, France, March 12–17, 1984_. Springer, 2005, pp. 156–163.
  * [14] ——, “Time reversal on wiener space,” in _Stochastic Processes—Mathematics and Physics: Proceedings of the 1st BiBoS-Symposium held in Bielefeld, West Germany, September 10–15, 1984_. Springer, 2006, pp. 119–129.
  * [15] P. Cattiaux, G. Conforti, I. Gentil, and C. Léonard, “Time reversal of diffusion processes under a finite entropy condition,” _arXiv preprint arXiv:2104.07708_ , 2021.
  * [16] A. Bensoussan, J. Frehse, P. Yam, _et al._ , _Mean field games and mean field type control theory_. Springer, 2013, vol. 101.
  * [17] M. Fornasier and F. Solombrino, “Mean-field optimal control,” _ESAIM: Control, Optimisation and Calculus of Variations_ , vol. 20, no. 4, pp. 1123–1152, 2014.
  * [18] R. Carmona, F. Delarue, _et al._ , _Probabilistic theory of mean field games with applications I-II_. Springer, 2018.
  * [19] I. Exarchos and E. A. Theodorou, “Learning optimal control via forward and backward stochastic differential equations,” in _2016 American Control Conference (ACC)_. IEEE, 2016, pp. 2155–2161.
  * [20] J. Han, A. Jentzen, and W. E, “Solving high-dimensional partial differential equations using deep learning,” _Proceedings of the National Academy of Sciences_ , vol. 115, no. 34, pp. 8505–8510, 2018.
  * [21] I. Exarchos and E. A. Theodorou, “Stochastic optimal control via forward and backward stochastic differential equations and importance sampling,” _Automatica_ , vol. 87, pp. 159–165, 2018.
  * [22] R. Archibald, “A stochastic gradient descent approach for stochastic optimal control,” _East Asian Journal on Applied Mathematics_ , vol. 10, no. 4, 2020\.
  * [23] M. Pereira, Z. Wang, T. Chen, E. Reed, and E. Theodorou, “Feynman-kac neural network architectures for stochastic control using second-order fbsde theory,” in _Learning for Dynamics and Control_. PMLR, 2020, pp. 728–738.
  * [24] E. Pardoux and P. Protter, “A two-sided stochastic integral and its calculus,” _Probability theory and related fields_ , vol. 76, no. 1, pp. 15–49, 1987.
  * [25] L. Ambrosio, N. Gigli, and G. Savaré, _Gradient flows: in metric spaces and in the space of probability measures_. Springer Science & Business Media, 2008.