Stochastic Reachability of Uncontrolled Systems via Probability
  Measures: Approximation via Deep Neural Networks
This paper poses a theoretical characterization of the stochastic
reachability problem in terms of probability measures, capturing the
probability measure of the state of the system that satisfies the reachability
specification for all probabilities over a finite horizon. We achieve this by
constructing the level sets of the probability measure for all probability
values and, since our approach is only for autonomous systems, we can determine
the level sets via forward simulations of the system from a point in the state
space at some time step in the finite horizon to estimate the reach
probability. We devise a training procedure which exploits this forward
simulation and employ it to design a deep neural network (DNN) to predict the
reach probability provided the current state and time step. We validate the
effectiveness of our approach through three examples.
Karthik Sivaramakrishnan, Vignesh Sivaramakrishnan, Rosalyn Alex Devonport, Meeko M. K. Oishi
# Stochastic Reachability of Discrete-Time Stochastic Systems via Probability Measures

###### Abstract

We develop a framework for stochastic reachability for discrete-time systems in terms of probability measures. We reframe the stochastic reachability problem in terms of reaching and avoiding known probability measures, and seek to characterize a) the initial set of distributions which satisfy the reachability specifications, and b) a Markov based feedback control that ensures these properties for a given initial distribution. That is, we seek to find a Markov feedback policy which ensures that the state probability measure reaches a target probability measure while avoiding an undesirable probability measure over a finite time horizon. We establish a backward recursion for propagation of state probability measures, as well as sufficient conditions for the existence of a Markov feedback policy via the 1-Wasserstein distance. Lastly, we characterize the set of initial distributions which satisfy the reachability specifications. We demonstrate our approach on an analytical example.

## I Introduction

Stochastic reachability analysis is an established method for probabilistic safety, that provides an assurance that states that start within some initial set can reach a desired target set while avoiding a "bad" set, with at least some known likelihood. Such an assurance is important for stochastic dynamical systems in which satisfaction of state and input constraints is paramount, including problems in space vehicle rendezvous [1] and robotics [2, 3]. A theoretical foundation based in dynamic programming has been developed for stochastic reachability of discrete-time stochastic hybrid systems [4, 5]. However, the computational hurdles associated with dynamic programming are significant, as it requires gridding the state space and evaluating a value function at all points on the grid. Alternative approaches have been developed, based in abstraction [6, 7], sampling [8], and underapproximative methods [9] as well as approximative methods [10] in convex optimization.

We propose an alternative technique for stochastic reachability that is based in probability measures. Probability measures can be used to represent not only the distributions of the state, but also whether a state lies within some region of interest. Probability measures have been employed in covariance steering [11, 12, 13] and distribution steering [14, 15] problems, in which the optimization problems minimize the distance between probability measures, as well as in data-driven reachability via Christoffel functions [16]. One advantage of employing probability measures is that we can exploit tools from measure theory to enforce guarantees that hold almost surely, meaning that properties over a probability measure hold with probability one. Further, recent work in optimal transport has yielded significant gains in computing with probability measures [17].

In this paper, we develop a theoretical foundation for propagation of probability measures that represent the distribution of the state of a discrete-time stochastic system over a finite horizon, and synthesize a Markov feedback policy that satisfies the reachability specifications. The Markov feedback policy minimizes the Wasserstein distance between distributions at each time step. In contrast to standard methods (Figure 1), which seek to synthesize the largest set of initial conditions which satisfy reachability specifications with at least a minimum likelihood, our approach assures reachability specifications almost surely. We obtain the set of initial probability measures, as well as a corresponding Markov feedback policy, in lieu of an initial set. Our approach is

Fig. 1: Top: The standard stochastic reach-avoid problem finds the set of initial states for which there exists a Markov feedback policy that ensures that the state avoids an unsafe set (red) and reaches a target set (purple). Bottom: We re-frame the stochastic reach-avoid problem in terms of probability measures, to find the set of initial probability measures (orange) for which a Markov feedback policy ensures that the state probability measure avoids an unsafe probability measure (red) and reaches a desired target probability measure (purple) almost surely.
enabled by the representation of standard target and avoid sets as probability measures.

We describe mathematical preliminaries from probability and measure theory in Section II. In Section III, we describe the recursive evolution of probability measures over a finite-time horizon, presuming a known Markov feedback policy, and in Section IV we synthesize a controller which is assured to satisfy the reachability specifications almost surely. We demonstrate our approach on a small, analytic example in Section V, and conclude in Section VI.

## II Preliminaries and Problem Formulation

Let \(\mathbb{N}\) denote the natural numbers and \(\mathbb{N}_{[a,b]}\) denote the set of natural numbers from \(a\) to \(b\) inclusively, for \(a,b\in\mathbb{N}\) and \(a\leq b\)[18]. Let \(\mathbb{R}^{p}\) describe the extended real numbers \(\mathbb{R}^{p}\cup\{-\infty,\infty\}^{p}\) of dimension \(p\in\mathbb{N}\). Matrices and vectors are denoted with uppercase \(D\in\mathbb{R}^{l\times p}\) and lowercase \(d\in\mathbb{R}^{l}\), respectively. Random vectors are bolded, \(\mathbf{y}\in\mathbb{R}^{p}\).

A probability space is a tuple \((\Omega,\mathscr{F},\mu)\), where \(\Omega\) is the sample space and \(\mathscr{F}\) is the \(\sigma\)-algebra containing measurable subsets of \(\Omega\), i.e. events in the sample space, and a probability measure \(\mu:\mathscr{F}\to[0,1]\). We consider the set of reals \(\mathbb{R}^{p}\) that is separable and metrizable, that is, a Polish space [19, Conventions]. A continuous random vector \(\mathbf{y}=[\mathbf{y}_{1}\cdots\mathbf{y}_{p}]^{\intercal}\) is a measurable function \(\mathbf{y}:\Omega\to\mathbb{R}^{p}\) which maps from the probability space to a Borel measurable space \((\mathbb{R}^{p},\mathscr{B}(\mathbb{R}^{p}))\), where \(\mathscr{B}(\mathbb{R}^{p})\) denotes the Borel \(\sigma\)-algebra and the inverse of the random vector is a subset of the \(\sigma\)-algebra, i.e., \(\mathbf{y}^{-1}(\mathscr{B}(\mathbb{R}^{p}))=\{\omega:\omega\in\Omega, \mathbf{y}(\omega)\in\mathscr{B}(\mathbb{R}^{p})\}\).

Since we assume Borel measurable probability spaces on the set of reals, we utilize the cumulative distribution function (cdf) as the probability measure, \(\mu_{\mathbf{y}}\), to define a probability space, \((\mathbb{R}^{p},\mathscr{B}(\mathbb{R}^{p}),\mu_{\mathbf{y}})\). The following Definition and Theorem establish this relationship.

**Definition 1** (Cumulative Distribution Function [20, Sec. 6.3]).: _A function \(\Phi:\mathbb{R}^{p}\to[0,1]\subset\mathbb{R}\) is a cdf when it satisfies the following conditions:_

1. _The left and right limits of_ \(\Phi\) _are zero and one respectively, that is,_ \(\lim_{y\to-\infty}\Phi(y_{1},\ldots,y_{p})=0\) _and_ \(\lim_{y\to\infty}\Phi(y_{1},\ldots,y_{p})=1\)_,_
2. _For a product of intervals,_ \(\times_{i=1}^{p}(a_{i},b_{i}]\subset\mathbb{R}^{p}\)_, the mixed monotonicity of the function is non-negative, i.e._ \[\Delta_{p}^{(a,b)}\Phi= \Phi(b_{1},\ldots,b_{p})\] \[-\sum_{j=1}^{p}\Phi(b_{1},\ldots,b_{j-1},a_{j},b_{j+1},\ldots,b_{ p})\] \[+\sum_{1\leq j<k\leq p}\Phi(b_{1},\ldots,b_{j-1},a_{j},b_{j+1},\ldots,\] \[b_{k-1},a_{k},b_{k+1},\ldots,b_{p})\] \[+(-1)^{p}\Phi(a_{1},\ldots,a_{p})\geq 0\] (1)

_where \(y=[y_{1}\cdots y_{p}]^{\intercal}\in\mathbb{R}^{p}\). Further, we denote for ease of notation \(\Phi(y_{1},\ldots,y_{p})=\Phi(y)\)._

**Theorem 1** (Defining a cdf for a real random vector on a probability space).: _For any random vector \(\mathbf{y}\) on the probability space \((\mathbb{R}^{p},\mathscr{B}(\mathbb{R}^{p}),\mu_{\mathbf{y}})\), there exists a unique cdf, \(\Phi_{\mathbf{y}}\), that satisfies Definition 1, where_

\[\mu_{\mathbf{y}}((-\infty,y_{1}]\times\cdots\times(-\infty,y_{p}]) \tag{2a}\] \[=\mathbb{P}\{\mathbf{y}_{1}(\omega)\leq y_{1}\cap\cdots\cap \mathbf{y}_{p}(\omega)\leq y_{p}\}\] (2b) \[=\lim_{a\to\infty}\Delta_{p}^{(a,b)}\Phi_{\mathbf{y}}=\Phi_{ \mathbf{y}}(y) \tag{2c}\]

Proof.: From [20, Sec. 6.1 and 6.3]. 

Since cdfs are unique, the converse of Theorem 1 allows one to define a random vector and its probability space, \((\mathbb{R}^{p},\mathscr{B}(\mathbb{R}^{p}),\mu_{\mathbf{y}})\), given a cdf.

**Corollary 1** (Defining a probability space via cdf).: _If a function \(\Phi\) satisfies Definition 1, then there always exists a random vector \(\mathbf{y}\) on a probability space \((\mathbb{R}^{p},\mathscr{B}(\mathbb{R}^{p}),\mu_{\mathbf{y}})\) such that \(\Phi\) is its cdf, that is, \(\Phi=\Phi_{\mathbf{y}}\)._

Proof.: From [20, Sec. 6.1 and 6.3]. 

Thus, we can either define a random vector \(\mathbf{y}\) such that there exists a corresponding cdf, \(\Phi_{\mathbf{y}}(y)\), or define a cdf, \(\Phi(y)\), such that it is the cdf of the random vector \(\mathbf{y}\), that is, \(\Phi_{\mathbf{y}}(y)=\Phi(y)\).

Finally, we say a function \(g:\mathbb{R}^{p}\to\mathbb{R}^{l}\) is a Borel measurable function as long as it is continuous almost everywhere, meaning that continuity holds everywhere except a measurable subset of measure zero [21, Sec. 2.2].

### _Discrete-time Stochastic Dynamical System_

We define a discrete-time, stochastic system that evolves over a time horizon of \(N\) steps, for time steps \(k\in\mathbb{N}_{[0,N-1]}\).

**Definition 2**.: _A discrete-time, stochastic system is a time-varying, Borel measurable mapping, \(f_{k}:\mathbb{R}^{n}\times\mathcal{U}\times\mathcal{W}_{k}\to\mathbb{R}^{n}\),_

\[x_{k+1}=f_{k}(x_{k},u_{k},w_{k})=[f_{1,k}(x_{k},u_{k},w_{k})\cdots\] \[f_{n,k}(x_{k},u_{k},w_{k})]^{\intercal}, \tag{3}\]

_where:_

1. \(\mathcal{U}\subseteq\mathbb{R}^{m}\) _is the set of admissible inputs that is compact, where_ \(u_{k}\in\mathcal{U}\) _is the input vector._
2. \(\mathcal{W}=\{(\mathcal{W}_{k},\mathscr{B}(\mathcal{W}_{k}))\}_{k=0}^{N-1}\)_, where_ \(\mathcal{W}_{k}\subseteq\mathbb{R}^{p}\)_, is the set of Borel measurable spaces denoting the disturbance space, with random vector_ \(\mathbf{w}_{k}\) _on a probability space_ \((\mathcal{W}_{k},\mathscr{B}(\mathcal{W}_{k}),\mu_{\mathbf{w}_{k}})\) _where_ \(\Phi_{\mathbf{w}_{k}}(w_{k})\) _is its cdf._
3. \(\Phi_{\mathbf{x}_{k+1}}:\mathscr{B}(\mathbb{R}^{n})\times\mathbb{R}^{n}\times \mathbb{R}^{n}\times\mathcal{W}_{k}\times\mathcal{U}\to[0,1]\) _is a cdf representing the time-varying transition kernel for a random vector_ \(\mathbf{x}_{k+1}\) _on the probability space_ \((\mathbb{R}^{n},\mathscr{B}(\mathbb{R}^{n}),\mu_{\mathbf{x}_{k+1}})\)_, given the current state_ \(x_{k}\in\mathbb{R}^{n}\)_, disturbance_ \(\mathbf{w}_{k}\)_, and input_ \(u_{k}\in\mathcal{U}\)_._

_We define the transition kernel in terms of the probability measure of the random vector \(\mathbf{w}_{k}\)[22, Ch. 8],_

\[\Phi_{\mathbf{x}_{k+1}}(x_{k+1},x_{k},u_{k})\] \[=\mu_{\mathbf{w}}(\{x_{k+1}\in\mathbb{R}^{n}:f_{1,k}(x_{k},u_{k}, \mathbf{w}_{k})\leq x_{1,k+1}\times\] \[\qquad\qquad\qquad\ldots\times f_{n,k}(x_{k},u_{k},\mathbf{w}_{k}) \leq x_{n,k+1}\}|x_{k},u_{k}) \tag{4a}\] \[=\Phi_{f_{k}(x_{k},u_{k},\mathbf{w}_{k})}(x_{k+1},x_{k},u_{k}). \tag{4b}\]
We also presume a Markov feedback policy \(\pi_{k}:\mathbb{R}^{n}\to\mathcal{U}\) that depends on the current state, \(x_{k}\).

**Definition 3** ([22,, Sec. 8.1, Def. 8.2]).: _A Markov feedback policy is a sequence \(\pi=\{\pi_{0},\ldots,\pi_{k},\ldots,\pi_{N-1}\}\) such that each Markov feedback map, \(\pi_{k}:\mathbb{R}^{n}\to\mathcal{U}\) for \(k\in\mathbb{N}_{[0,N-1]}\), is a time-varying, non-random, Borel measurable function, that is, \(u_{k}=\pi_{k}(x_{k})\)._

### _Set Representations via Measure Theory_

We use tools from measure theory to represent a change of measure and probabilities of a random vector \(\mathbf{y}\) residing in a Borel measurable set via probability measures.

**Theorem 2** (Radon-Nikodym Theorem, [20]).: _If \(\nu\) and \(\mu\) are \(\sigma\)-finite measures (e.g., probability measures) on the Borel measurable space \((\mathbb{R}^{p},\mathscr{B}(\mathbb{R}^{p}))\), where \(\mu\) is absolutely continuous with respect to \(\nu\), that is, \(\mu<<\nu\), then there exists a measurable function \(g(y)\geq 0\) such that_

\[\mu(G)=\int_{G}g(y)\mathrm{d}\nu(y),\ \forall G\in\mathscr{B}(\mathbb{R}^{p}). \tag{5}\]

_The function \(g(y)=\frac{\mathrm{d}\mu}{\mathrm{d}\nu}\) is the Radon-Nikodym derivative of \(\mu\) with respect to \(\nu\)._

Theorem 2 imposes a change of measure via a non-negative measurable function. We make use of the the Dirac measure, a probability measure over a set, as such a non-negative measurable function.

**Definition 4** (Dirac Measures [21, Sec. 1.3]).: _Given a probability space \((\mathbb{R}^{n},\mathscr{B}(\mathbb{R}^{n}))\), a Dirac measure is \(\delta_{y}:G\to\{0,1\}\) where_

\[\delta_{y}(G)=\begin{cases}1,\ y\in\mathcal{G}\\ 0,\ y\notin\mathcal{G}\end{cases}. \tag{6}\]

Note that the Dirac measure functions similarly to an indicator function (\(1_{G}(y)=1\) if \(y\in G\), and is 0 otherwise).

**Corollary 2** (Random Vector Residing in a Set via Dirac Measure).: _Given the probability measure of a random vector residing within a set, \(\mu(G)=\mathbb{P}\{\mathbf{y}\in G\}\), and the cdf of \(\mathbf{y}\), \(\Phi_{\mathbf{y}}\), such that \(\mu(G)<<\Phi_{\mathbf{y}}\), then the Dirac measure of a point \(y\) residing in a set, \(\delta_{y}(G)\), is a non-negative measurable function such that_

\[\mu(G)=\int_{\mathbb{R}^{p}}\delta_{y}(G)\mathrm{d}\Phi_{\mathbf{y}}(y). \tag{7}\]

Proof.: This follows directly from Theorem 2. 

### _Problem Formulation_

We seek to find the set of initial probability measures for which the state of the stochastic system in Definition 2, with some Markov feedback policy \(\pi\), will satisfy

\[\mathbb{P}\left\{\mathbf{x}_{N}\in\mathcal{T}\cap\left(\cap_{k=0}^{N-1} \mathbf{x}_{k}\not\in\mathcal{A}\right)\right\} \tag{8}\]

almost surely. In words, we wish for the system to reside in a target region, \(\mathcal{T}\), at final time \(N\), while avoiding an unsafe region, \(\mathcal{A}\), for \(k\in\mathbb{N}_{[0,N-1]}\) with probability 1 almost everywhere. Equation (8) imposes reaching a target region, \(\mathcal{T}\), while avoiding an undesired region, \(\mathcal{A}\), by staying within a target tube that is its complement, \(\mathcal{A}^{\mathcal{G}}\)[9].

We first transform these set representations to probability measures of the target set and the target tube via Definition 4.

\[\mu_{\mathcal{T}}(x_{N}) =\delta_{x_{N}}(\mathcal{T})\] (9a) \[\mu_{\mathcal{A}^{\mathcal{G}}}(x_{k}) =\delta_{x_{k}}(\mathcal{A}^{\mathcal{\bar{\mathbf{\mathbf{\mathbf{ \mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{ \mathbf{      \mathbf{      }}}}}}}}}}}})\] (\ref{\mu_{\mathcal{A}}(x_{k})=\delta_{x_{k}}(\mathcal{A}^{\mathcal{ \bar{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{ \mathbf{\mathbf{\mathbf{\mathbf{\mathbf{ }}}}}}}}}}}}})\)\] (\ref{\mu_{\mathcal{A}}(x_{k})=1-\mu_{\mathcal{A}^{\mathcal{\bar{ \mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{ \mathbf{\mathbf{ }}}}}}}}}}}}}})\] (\ref{\mu_{\mathcal{A}}(x_{k})=1}\). Note that this transformation is readily extendable to random sets [23].

**Problem 1** (Stochastic Reach-Avoid via Probability Measures).: _Given a target probability measure, \(\mu_{\mathcal{T}}\), and an avoid probability measure, \(\mu_{\mathcal{A}}\), we seek to find \(P_{\mathbf{x}_{0}}\), the set of initial probability measures, for which there exists a Markov feedback policy such that the probability measure of state of the system in Definition 2 reaches \(\mu_{\mathcal{T}}\) at the final time \(N\) and avoids \(\mu_{\mathcal{A}}\) for all time steps \(k\in\mathbb{N}_{[0,N-1]}\), almost surely._

We solve Problem 1 in two steps. First, in Section III, we presume that a Markov feedback policy exists and we recursively construct the state probability measures which satisfy the reach-avoid specifications almost surely backwards-in-time. Second, in Section IV, we construct Markov feedback policies for every initial probability measure that satisfies the reach-avoid specification almost surely, then construct the set of initial probability measures, \(P_{\mathbf{x}_{0}}\).

## III Backward Recursion of the State Probability Measure

We seek to characterize the backward recursion that enables propagation of probability measures (Figure 2) that satisfy the target and avoid probability measures. To do so, we first describe the state probability measure \(\mathbf{x}_{k+1}\) in terms of the cdf of the prior state, \(\mathbf{x}_{k}\), using the relationship in (4).

\[\Phi_{\mathbf{x}_{k+1}}(x_{k+1})\\ =\int_{\mathbb{R}^{n}}\Phi_{f_{k}(x_{k},u_{k},\mathbf{w}_{k})}(x_{ k+1},x_{k},\pi_{k}(x_{k}))\mathrm{d}\Phi_{\mathbf{x}_{k}}(x_{k}) \tag{10}\]

**Lemma 1** (Uniqueness of State Probability Measure at Time \(k\)).: _For a probability measure \(\Phi_{\mathbf{x}_{k+1}}\) and a known Markov feedback map \(\pi_{k}(x_{k})\), there exists a unique probability measure \(\Phi_{\mathbf{x}_{k}}\) such that (10) holds._

Proof.: From [21, Lemma 10.4.3]. 

Satisfaction of the avoid and target probability measures requires additional constraints to be placed on this propagation. To constrain the state probability measure \(\Phi_{\mathbf{x}_{k}}\) so that it avoids the avoid probability measure \(\mu_{\mathcal{A}}\), we seek to satisfy

\[\int_{\mathbb{R}^{n}}\mu_{\mathcal{A}^{\mathcal{\bar{\mathbf{\mathbf{\mathbf{\mathbf{ \mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbf{\mathbfmathbfmathbfmathbf{\mathbfmathbfmathbfmathbfmathbf{\mathbfmathbfmathbf{\mathbfmathbfmathbf{ \mathbfmathbfmathbf{ \mathbfmathbfmathbfmathbfmathbfmathbfmathbfmathbfmathbfmathbfmathbf \mathbfmathbfmathbf{\leftleftleftleftleftleftleftleftleftleftleftleft \leftleftleftleftleftleftleftleftleftleft {ileftleftleftleftleftleftleftleftleftleftleftleftleftleft {i\leftleftleftleftleftleftleftleftleftleftleftleftleft {i{\leftleftleftleftleftleftleftleftleftleftleftleftleftleftleftleft {\leftleftleftleftleftleftleftleftleftleftleftleft {\leftleftleftleftleftleftleftleftleftleft {\leftleftleftleftleftleftleftleftleft {\leftleftleftleftleftleftleftleftleftleftleft {\leftleftleftleftleftleftleftleftleftleftleftleft {\leftleftleftleftleftleftleftleftleftleft \leftleftleftleftleftleftleft {\leftleftleftleftleftleftleftleftleftleft \leftleftleftleftleftleftleft { \leftleftleftleftleftleftleftleft { \leftleftleftleftleftleftleftleftleftleftleft { { \leftleftleftleftleftleftleftleftleftleftleftleftleft  { \leftleftleftleftleftleftleftleftleftleft  {    \right} \right}}\right}\right}\)\))))))))))))))))))))))))))))))))))))))))))}}}}}}}}\)))))))((((((((((((((((((((((((((((((()))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))}}}}}}}}}((())))))))))))))))))}
For the state probability measure to meet the target set requirement, we must constrain \(\Phi_{\mathbf{x}_{N-1}}\) as follows,

\[\int_{\mathbb{R}^{n}}\int_{\mathbb{R}^{n}}\mu_{\mathcal{T}}(x_{N}) \mathrm{d}\Phi_{\mathbf{\hat{x}}_{N}}(x_{N},x_{N-1},\pi_{N-1}(x_{N-1}))\\ \cdot\mathrm{d}\Phi_{\mathbf{x}_{N-1}}(x_{N-1})=1, \tag{12}\]

for \(\hat{\mathbf{x}}_{N}=f_{N-1}(x_{N-1},u_{N-1},\mathbf{w}_{N-1})\). We formalize these notions in the following Theorem.

**Theorem 3**.: _Given a target probability measure \(\mu_{\mathcal{T}}\), a target tube probability measure \(\mu_{\mathcal{A}}\), and a Markov feedback policy \(\pi\), the state probability measure must satisfy (10) for \(k\in\mathbb{N}_{[0,N-1]}\), as well as_

1. _(_11_) for_ \(k\in\mathbb{N}_{[0,N-2]}\)_,_
2. _(_12_) for_ \(k=N-1\)_,_

_in order to satisfy the reachability constraints almost surely._

Proof.: Consider the case in which \(k=N-1\): Note that (12) holds as long as \(\Phi_{\mathbf{x}_{N-1}}(x_{N-1})\) satisfies (11).

Now consider \(k\in\mathbb{N}_{[0,N-2]}\). We employ backwards induction, beginning with the base case \(k=N-2\),

\[\Phi_{\mathbf{x}_{N-1}}(x_{N-1})\\ =\int_{\mathbb{R}^{n}}\Phi_{f_{N-2}(x_{N-2},u_{N-2},\mathbf{w}_{N -2})}(x_{N-1},x_{N-2},\pi_{N-2}(x_{N-2}))\\ \cdot\mathrm{d}\Phi_{\mathbf{x}_{N-2}}(x_{N-2}), \tag{13}\]

holds almost everywhere as long as \(\Phi_{\mathbf{x}_{N-2}}(x_{N-2})\) satisfies (11). If the equality holds for the case \(k=j\) where \(j<N-2\), then it must also hold for the case \(k=j-1\). Observe that for \(k=j-1\) and \(k=j\) the recursions, respectively,

\[\Phi_{\mathbf{x}_{j+1}}(x_{j+1})=\int_{\mathbb{R}^{n}}\Phi_{f_{j} (x_{j},u_{j},\mathbf{w}_{j})}(x_{j+1},x_{j},\pi_{j}(x_{j}))\\ \cdot\mathrm{d}\Phi_{\mathbf{x}_{j}}(x_{j})\\ \Phi_{\mathbf{x}_{j}}(x_{j})\\ =\int_{\mathbb{R}^{n}}\Phi_{f_{j-1}(x_{j-1},u_{j-1},\mathbf{w}_{ j-1})}(x_{j},x_{j-1},\pi_{j-1}(x_{j-1}))\\ \cdot\mathrm{d}\Phi_{\mathbf{x}_{j-1}}(x_{j-1}), \tag{14b}\]

hold almost everywhere as long as \(\Phi_{\mathbf{x}_{j}}(x_{j})\) and \(\Phi_{\mathbf{x}_{j-1}}(x_{j-1})\) satisfy (11). Thus, the equality holds for all \(k\leq N-2\). 

With establishing the conditions for which the backward recursion exists, we also can form a _certificate_ if the Markov feedback policy does not ensure that the state probability measure satisfies the reachability specifications almost surely.

**Definition 5** (Certificate of Infeasibility).: _A backward recursion with a given Markov feedback policy provides a certificate of infeasibility at time step \(k\in\mathbb{N}_{[0,N]}\) when the recursion in Theorem 3 fails. That is, it cannot ensure a state probability measure, \(\Phi_{\mathbf{x}_{k}}\), that satisfies (11) for \(k\in\mathbb{N}_{[0,N-2]}\) or (11) and (12) for \(k=N-1\)._

As a result, when a given Markov feedback policy does not satisfy Definition 5, it may be an indication that the feedback policy should be modified to allow more control effort, for example, to ensure satisfaction.

## IV Feedback controller synthesis via probability measures

With a complete description of the machinery of propagating probability measures via a backward recursion, we can now consider synthesis of Markov feedback maps in this same framework. First, we note that one of the main theoretical challenges associated with this framework is the need to enforce equality of probability measures almost everywhere. In particular, the formulation of an optimal Markov feedback policy within this framework requires a description of distances between probability measures. For this reason, we employ the Wasserstein distance from optimal transport [19] as it allows us to match two probability measures by measuring the distance between them.

**Theorem 4** (Kantorovich-Rubinstein Theorem [19], Ch. 6).: _Let \((\mathcal{M},\mathscr{B}(\mathcal{M}))\) denote a Polish space. Let \(\mu_{1}\) and \(\mu_{2}\) represent probability measures on \((\mathcal{M},\mathscr{B}(\mathcal{M}))\), that is, \(\mu_{1},\mu_{2}\in\mathscr{B}(\mathcal{M})\). Thus, \(\forall\mu_{1},\mu_{2}\in\mathscr{B}(\mathcal{M})\), the dual of the 1-Wasserstein distance can be written as_

\[W_{1}(\mu_{1},\mu_{2})\\ =\sup_{\|g\|_{\mu}\leq 1}\left\{\int_{\mathcal{M}}g(x)\mathrm{d} \mu_{1}(x)-\int_{\mathcal{M}}g(x)\mathrm{d}\mu_{2}(x)\right\}, \tag{15}\]

_where \(g:\mathcal{M}\to\mathbb{R}_{\geq 0}\) that satisfy the Lipschitz condition \(|g(x_{1})-g(x_{2})|\leq D(x_{1},x_{2})\) for all \(x_{1},x_{2}\in\mathcal{M}\)._

Theorem 4 allows us to measure the largest point of deviation between two probability measures via this distance

Fig. 2: The backward recursion (10) and constraints (11), (12) together ensure that the target probability measure can be reached and the avoid probability measure avoided, subject to stochastic system dynamics and a known controller.
metric. In addition, the \(1\)-Wasserstein distance exhibits the property of lower semicontinuity, which is an essential property that allows us to optimize over functions.

**Lemma 2** (Lower semicontinuity of Wasserstein [19], Ch. 6).: _Let \((\mathcal{M},\mathscr{B}(\mathcal{M}))\) be a Polish space which is a complete, separable metric space with Borel \(\sigma\)-algebra \(\mathscr{B}(\mathcal{M})\) and \(W_{1}:\mathcal{M}\to\mathbb{R}_{\geq 0}\). The 1-Wasserstein distance in Theorem 4 is lower semicontinuous on \(\mathscr{B}(\mathcal{M})\), if for any sequence of probability measures \(\mu_{i}\) and \(\nu_{i}\) on \(\mathscr{B}_{1}(\mathcal{M})\) that converges weakly to probability measures \(\mu\) and \(\nu\) on \(\mathscr{B}(\mathcal{M})\), then we have that \(W_{1}(\mu,\nu)\leq\liminf\limits_{i\to\infty}\,W_{1}(\mu_{i},\nu_{i})\)._

Proof.: Given by Theorem 4.1 and Remark 6.12 in [19]. 

Lemma 2 essentially states that, if a sequence of probability measures \((\mu_{i},\nu_{i})\) weakly converges to another pair of probability measures \((\mu,\nu)\), then \(W_{1}(\mu_{i},\nu_{i})\) is an infimum to \(W_{1}(\mu,\nu)\) as \(i\in\mathbb{N}\) goes to \(\infty\). Therefore, the \(1-\)Wasserstein distance is lower semicontinuous, meaning that it is possible to optimize over this distance metric.

### _Controller Synthesis_

We now apply Theorem 4 to the left and right sides of the equalities in equations (10) and (12), to obtain (16) and (17), respectively. Equations (16) and (17) can be used to construct a Markov feedback policy, such that (10) and (12), respectively, are satisfied. Feedback policies that solve (16) and (17) are optimal, in that they minimize the 1-Wasserstein distance. (Note that this is a different notion of optimality of the control than in standard stochastic reachability, in which an optimal control minimizes a cost function that encodes a a reachability-specific cost function.)

**Theorem 5** (Existence of an Optimal Markov Feedback Policy).: _Consider a recursion of \(\Phi_{\mathbf{x}_{k}}\) that satisfies (16) for \(k\in\mathbb{N}_{[0,N-2]}\) and (17) for \(k=N-1\). If \(\pi_{k}^{*}(x_{k}):\mathbb{R}^{n}\to\mathcal{U}\) satisfies_

\[\pi_{k}^{*}(x_{k})=\operatorname*{arg\,inf}_{\pi_{k}\in\mathcal{U}}\,W_{1}\left( \Phi_{f_{k}(x_{k},u_{k},\mathbf{w}_{k})},\Phi_{\mathbf{x}_{k+1}}\right) \tag{18}\]

_for every \(\Phi_{\mathbf{x}_{k}}\) that satisfies (11) for \(k\in\mathbb{N}_{[0,N-2]}\), and if \(\pi_{N-1}^{*}(x_{N-1}):\mathbb{R}^{n}\to\mathcal{U}\) satisfies_

\[\pi_{N-1}^{*}(x_{N-1})=\operatorname*{arg\,inf}_{\pi_{N-1}\in \mathcal{U}}\,W_{1}\left(\mu_{\mathcal{T}}(x_{N})\Phi_{\mathbf{\hat{x}}_{N}}, \mu_{\mathcal{T}}(x_{N})\right), \tag{19}\]

_where \(\mathbf{\hat{x}}_{N}=f_{N-1}(x_{N-1},u_{N-1},\mathbf{w}_{N-1})\) for every \(\Phi_{\mathbf{x}_{N-1}}\) that satisfies (11) and (12) for \(k=N-1\), then \(\pi^{*}=\{\pi_{0}^{*}(x_{0}),\pi_{1}^{*}(x_{1}),\ldots,\pi_{N-1}^{*}(x_{N-1})\}\) is an optimal Markov feedback policy, where each \(\pi_{k}^{*}(x_{k})\) is an optimal Markov feedback map._

Proof.: We first consider the \(k=N-1\) time step. Through lower semicontinuity of the Wasserstein distance, provided via Lemma 2, we know that since the optimization in (19) is sound, there must exist an optimal Markov feedback map \(\pi_{N-1}^{*}(x_{N-1})\in\mathcal{U}\) as in (19) for every \(\Phi_{\mathbf{x}_{N-1}}\) that satisfies (11). Therefore, \(\Phi_{\mathbf{x}_{N}}\) and \(\Phi_{\mathbf{x}_{N-1}}\) will satisfy (12) for \(\pi_{N-1}^{*}(x_{N-1})\) that minimizes (17).

For \(k\in\mathbb{N}_{[0,N-2]}\), we again exploit lower semicontinuity of the Wasserstein distance via Lemma 2, to note that there must exist an optimal Markov feedback map \(\pi_{k}^{*}(x_{k})\in\mathcal{U}\) that optimizes (16). Hence for every \(\Phi_{\mathbf{x}_{k}}\) that satisfies (11), \(\Phi_{\mathbf{x}_{k+1}}\) and \(\Phi_{\mathbf{x}_{k}}\) will satisfy (10). 

With existence of an optimal Markov feedback policy assured, we now turn to synthesis of the policy.

**Corollary 3**.: _Given an optimal Markov feedback policy \(\pi^{*}\) as in Theorem 5, if the Wasserstein distance in (16) and (17) is zero, then \(\pi^{*}\) satisfies the state probability measure reaching \(\mu_{\mathcal{T}}\) at final time \(N\), while staying within \(\mu_{\mathcal{A}^{\mathsf{B}}}\) for \(k\in\mathbb{N}_{[0,N-1]}\) almost surely._

Proof.: Using a similar logic as in Lemma 2, note that for a sequence of state probability measures, we assume an optimal Markov feedback policy that guarantees the Wasserstein distance is zero. Thus, the state probability measures on the left and right sides of (10) and (12), respectively, are equivalent. With this equivalence, the target probability measure, \(\mu_{\mathcal{T}}\), and the target tube measure, \(\mu_{\mathcal{A}^{\mathsf{B}}}\), are satisfied almost surely. 

Corollary 3 ensures that the state probability measure reaches the target probability measure at time step \(k=N\)
and avoids the avoid probability measure for \(k\in\mathbb{N}_{[0,N-1]}\) almost surely.

### _Characterizing the Set of Initial State Probability Measures_

Recall that Problem 1 seeks to find the set of initial state probability measures, \(P_{\mathbf{x}_{0}}\), for which there exists a Markov feedback policy such that the state probability measure satisfies the reach-avoid specification almost surely. Corollary 3 provides a means to generate the set \(P_{\mathbf{x}_{0}}\), which we capture via the following definition.

**Definition 6**.: _The set of initial state probability measures, \(P_{\mathbf{x}_{0}}\), consists of the probability measures \(\Phi_{\mathbf{x}_{0}}\) that satisfy Corollary 3._

Synthesis of the set \(P_{\mathbf{x}_{0}}\) is the focus of future work.

### _Relationship to Discrete-Time Stochastic Reachability_

Canonical stochastic reachability frameworks [4, 5, 23] construct a set of initial states for which there exists a Markov feedback policy that maximizes (8). The approach described here provides a stricter interpretation, by insisting that (8) be satisfied almost surely. Additionally, in contrast to a point-based or set-based evaluation, the approach proposed here is based on entire distributions. Lastly, the controller synthesis that we propose varies considerably from that addressed via existing approaches, in that instead of maximizing a likelihood of safety, our proposed controller minimizes violation of equality constraints that describe propagation of distributions in a manner consistent with the reachability specifications. We plan to investigate these relationship more precisely in future work.

## V Examples

Consider the time-discretized single integrator with additive disturbance over a horizon \(N=1\),

\[\mathbf{x}_{k+1}=\mathbf{x}_{k}+\Delta T\mathbf{u}_{k}+\mathbf{w}_{k}, \tag{20}\]

with sampling rate \(\Delta T=1\) and \(\mathbf{w}_{k}\) that follows a normal distribution \(\mathbf{w}_{k}\sim\mathcal{N}(m,\sigma)\), with mean \(m=0\) and standard deviation \(\sigma=1\). We assume the state probability measure is described by Gaussian probability measures, which are amenable to analytical solutions. For ease of calculation, we consider an affine feedback controller (as opposed to a general Markov controller) that preserves Gaussianity of the state probability measure,

\[\mathbf{u}_{k}=h_{k}\mathbf{x}_{k}+v_{k}, \tag{21}\]

where \(h_{k}\in\mathbb{R}\) is a scaling term and \(v_{k}\in\mathbb{R}\) is a feedforward term. We presume the input lies within \(u_{\min}=-0.1\) and \(u_{\max}=0.1\) via the probabilistic constraint,

\[\mathbb{P}\left\{\cap_{k=0}^{N-1}(u_{\min}\leq\mathbf{u}_{k}\leq u_{\max}) \right\}=1. \tag{22}\]

We combine (21) with (20) to form the closed-loop dynamics,

\[\mathbf{x}_{k+1} =\mathbf{x}_{k}+\Delta T(h_{k}\mathbf{x}_{k}+v_{k})+\mathbf{w}_{k}, \tag{23a}\] \[=(1+\Delta Thh_{k})\mathbf{x}_{k}+\Delta Tv_{k}+\mathbf{w}_{k}, \tag{23b}\]

with mean and variance

\[m_{\mathbf{x}_{k+1}} =(1+\Delta Thh_{k})m_{\mathbf{x}_{k}}+\Delta Tv_{k}+m_{\mathbf{w} _{k}}, \tag{24a}\] \[\sigma_{\mathbf{x}_{k+1}}^{2} =(1+\Delta Thh_{k})\sigma_{\mathbf{x}_{k}}^{2}+\sigma_{\mathbf{w} _{k}}^{2}. \tag{24b}\]
### _Non-Random Target and Avoid Sets_

For the target set \(\mathcal{T}=\{x\in\mathbb{R}:x\leq 1\}\), we construct the target probability measure \(\mu_{\mathcal{T}}(x)=\delta_{x}(\mathcal{T})\). For the avoid set with \(\mathcal{A}=(-0.5,0.5)\subset\mathbb{R}\), we construct the target tube probability measure \(\mu_{\mathcal{A}^{\complement}}(x)=\delta_{x}(\mathcal{A}^{\complement})\). Because we only consider one time step, we must assure that \(\Phi_{\mathbf{x}_{0}}\) satisfies (11) and \(\Phi_{\mathbf{x}_{1}}\) satisfies (12), almost surely. Since \(\int_{\mathbb{R}}\mu_{\mathcal{A}^{\complement}}(x_{0})\mathrm{d}\Phi_{\mathbf{ x}_{0}}(x_{0})=1\) for \(\mathbf{x}_{0}\sim\mathcal{N}(m_{\mathbf{x}_{0}},\sigma_{\mathbf{x}_{0}})\), we have two solutions: \(\Phi_{\mathbf{x}_{0}}^{(2)}\) and \(\Phi_{\mathbf{x}_{0}}^{(1)}\), where \(\Phi_{\mathbf{x}_{0}}^{(2)}=1-\bar{\Phi}_{\mathbf{x}_{0}}^{(2)}\) is the complement cdf, as shown in Figure 4. Similarly, with the affine control, we must ensure \(\int_{\mathbb{R}}\int_{\mathbb{R}}\mu_{\mathcal{T}}(x_{1})\mathrm{d}\Phi_{( \downarrow\Delta Th_{0})x_{0}+\Delta Tu_{0}+\mathbf{w}_{0}}(x_{1})\mathrm{d} \Phi_{\mathbf{x}_{0}}(x_{0})=1\). This is possible by the use of the normal approximation of Dirac measures, since the normal distribution is a Schwartz function [24]. That is, for some mean \(m_{\mathbf{x}_{0}}\), \(\lim_{\sigma_{\mathbf{x}_{0}}\to 0}\Phi\left(\frac{x_{0}-m_{\mathbf{x}_{0}}}{ \sigma_{\mathbf{x}_{0}}}\right)=\delta_{x}(\{x\in\mathbb{R}:m_{\mathbf{x}_{0}} \leq x\})\), where \(\Phi(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-t^{2}/2}\) is a standard Gaussian cdf. Hence we note that the two cases which satisfy (11) are \(m_{x_{0}}^{(1)}=0.5\) and \(m_{x_{0}}^{(2)}=-0.5\) with \(\sigma_{\mathbf{x}_{0}}\to 0\).

Conroller synthesis is accomplished by satisfying (12), to ensure the transition from \(\Phi_{\mathbf{x}_{0}}\) to \(\Phi_{\mathbf{x}_{1}}\). However, because control authority is bounded, it is not possible to satisfy (12). With the maximal value of control, the resulting state probability measure is still far from the desired probability measure. As shown in Figure 5 (left column), to achieve the desired measure (red), which is represented by a Dirac measure, we would need to have a probability measure with standard deviation that approaches 0 in the limit (right column). This would require \(v_{0}\) to be infinite. Note that \(h_{0}=0\) since both \(\Phi_{\mathbf{x}_{0}}^{(1)}\) and \(\Phi_{\mathbf{x}_{0}}^{(2)}\) utilize the normal approximation of a Dirac measure.

The infeasibility shown in Figure 5 is due to the fact that this analysis seeks to satisfy (12) _almost surely_, which is quite strict. We anticipate that for this analysis to have utility beyond certificates of infeasibility for control, it will be important to examine relaxations that allow for probability measures within some allowable distance of the desired probability measure.

### _Random Target and Avoid Sets_

We now consider target and avoid sets that are random, meaning that the set takes a probabilistic representation. For a Bernoulli random variable \(\mathbf{y}\) with \(p=0.2\) and \(q=0.8\), the target set is \(\mathcal{T}\) is random, with the set \(\mathcal{T}_{1}\) evaluated with 20% likelihood or \(\mathcal{T}_{2}\) evaluated with 80% likelihood, where \(\mathcal{T}_{1}=\{x\in\mathbb{R}:x\leq-0.5\}\) and \(\mathcal{T}_{2}=\{x\in\mathbb{R}:x\leq-1\}\). The resulting target probability measure is \(\mu_{\mathcal{T}}(x)=0.2\cdot\delta_{x}(\mathcal{T}_{1})+0.8\cdot\delta_{x}( \mathcal{T}_{2})\). For the target tube probability measure, consider the same \(\mathbf{y}\) such that \(\mu_{\mathcal{A}^{\complement}}(x)=0.2\cdot\delta_{x}(\mathcal{A}^{\complement} _{1})+0.8\cdot\delta_{x}(\mathcal{A}^{\complement}_{2})\), for \(\mathcal{A}^{\complement}_{1}=\{x\in\mathbb{R}:\ x\leq-0.5\ \text{or}\ x\geq 0.5\}\) and \(\mathcal{A}^{\complement}_{2}=\{x\in\mathbb{R}:\ x\leq-1\ \text{or}\ x\geq 1\}\).

This case follows a similar argument to the non-random sets, in which \(\Phi_{\mathbf{x}_{0}}\) satisfies (11) and \(\Phi_{\mathbf{x}_{1}}\) satisfies (12), almost surely. Thus, we have two possible state probability measures in the form of \(\Phi_{\mathbf{x}_{0}}(x_{0})=p\cdot\Phi_{\mathbf{x}_{0,1}}(x_{0})+q\cdot\Phi_{ \mathbf{x}_{0,2}}(x_{0})\) which is a mixture of Gaussian probability measures where for \(p=0.2\) with mean \(m_{\mathbf{x}_{0,1}}\) and standard deviation \(\sigma_{\mathbf{x}_{0,1}}\) and \(q=0.8\) with \(m_{\mathbf{x}_{0,2}}\) and \(\sigma_{\mathbf{x}_{0,2}}\). We use a similar normal approximation argument as in the non-random case, resulting in \(m_{\mathbf{x}_{1}}^{(1)}=0.5\), \(m_{\mathbf{x}_{0,2}}^{(1)}=1\) and \(m_{\mathbf{x}_{0,1}}^{(2)}=-0.5\), \(m_{\mathbf{x}_{0,2}}^{(2)}=-1\) with all \(\sigma_{\mathbf{x}_{0,1}},\sigma_{\mathbf{x}_{0,2}}\to 0\), as shown in Figure 6.

Similarly to the previous analysis, it is not possible to satisfy (12) with bounded control authority, as shown in Figure 7 (left column). Satisfaction of (12) would require the standard deviation of the mixture of Gaussian probability measure, \(\Phi_{\mathbf{x}_{1}}\), to approach zero, which would require \(v_{0}\) to be infinite and \(h_{0}=0\) (right column).

As in the non-random set case, we cannot assure (12) due to the affine controller and Gaussian probability measure, as shown in Figure 7 (left) for the value of \(v_{k}\) which saturates the control (with \(h_{0}=0\)). For (12) to be satisfied, we require \(v_{0}\rightarrow-\infty\) (right), which, as in the previous case, requires infinite control authority.

Fig. 4: For a one-step time horizon for a Gaussian system with Dirac measure target and avoid sets and an affine control, satisfying (11) results in a Gaussian with standard deviation \(\sigma_{\mathbf{x}_{0}}\to 0\) (right, blue solid). For small values of \(\sigma_{\mathbf{x}_{0}}\) (left and middle) discrepancies between the Gaussian (cyan dashed) and Dirac measure (red) that captures the target tube probability measure are evident.

Fig. 5: With bounded control authority, for each of the two initial probability measures, \(\Phi_{\mathbf{x}_{1}}^{(1)}\) and \(\Phi_{\mathbf{x}_{1}}^{(2)}\), it is not possible to construct an affine controller that satisfies (12). This is seen in the fact that the state probability measure \(\Phi_{\mathbf{x}_{1}}\) (left, cyan dashed) does not attain the desired value of \(1\) before the target probability measure (red) does (as \(x_{1}\) increases). In contrast, if we were to allow infinite control authority (right), the feedforward gain of \(v_{0}\rightarrow\infty\) (respectively \(v_{0}\rightarrow-\infty\)) would satisfy (12).
## VI Conclusion

We present a framework for stochastic reachability via probability measures. Given a Markov feedback policy, we establish the existence of the backward recursion of state probability measures. Then, we establish the conditions under which there exists a Markov feedback policy that minimizes a Wasserstein distance, ensuring satisfaction of propagation of the probability measures in a manner that respects reachability specifications. Future work will focus on relationships to the standard stochastic reachability problem and relaxations that exploit computational advantages of the Wasserstein distance.

## Acknowledgements

This work has been supported in part by the NSF under awards CNS-1836900 and CMMI-2105631 as well as by NASA under the University Leadership Initiative award #80NSSC20M0163. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF or any NASA entity.
