Probabilistic Reachability of Discrete-Time Nonlinear Stochastic Systems

Computing the reachable sets for deterministic systems is a fundamental problem in dynamical systems and control theory. Consider the deterministic system

which can be viewed as a noiseless version of the stochastic system (1). The deterministic reachable set (DRS) at time $t$ is the set in the state space that the system (2) can reach, starting from an initial configuration, under all possible inputs within a specified time.

In general, computing the exact DRS of a dynamic system is challenging [26]. Therefore, most methods in reachability analysis focus on providing over-approximation of DRS [1]. A set $\overline{\mathcal{R}}_{t}\subseteq\mathbb{R}^{n}$ is an over-approximation of the DRS $\mathcal{R}_{t}$ if

The calculation of $\overline{\mathcal{R}}_{t}$ depends on the properties of system dynamics $f$ [27], and types of input $u_{t}$ [5]. In Section 5, we review two approaches to compute $\overline{\mathcal{R}}_{t}$: Lipschitz bound reachability and interval-based reachability.

We are interested in characterizing the reachable set of (1). We focus on the setting where the disturbance $w_{t}$ is sub-Gaussian. The sub-Gaussian distribution is quite general and includes a wide range of distributions such as Gaussian, uniform, and any zero-mean distributions with bounded support.

The theory for deterministic reachability analysis is not effective in handling stochastic disturbances. When $w_{t}$ is unbounded sub-Gaussian noise, the DRS in the sense of (3) is often trivial. To see this, consider the system $X_{t+1}=X_{t}+w_{t}$ where $w_{t}$ is a zero-mean Gaussian noise. The state $X_{t}$ is a Gaussian vector whose possible range is the entire $\mathbb{R}^{n}$ space and thus $\mathcal{R}_{t}=\mathbb{R}^{n}$. When $w_{t}$ is bounded, the deterministic reachability analysis ignores the statistics of the disturbance and treats it in a worst-case manner, leading to conservative results. In reality, the worst-case realization of stochastic disturbance occurs with extremely low probability and is thus not representative. For these reasons, we turn to the concept of probabilistic reachable set, which refers to the set that any possible trajectory starting from the initial state set $x_{0}\in\mathcal{X}_{0}$ and with the input $u_{\tau}\in\mathcal{U}$ can reach with high probability.

Apparently, as can be seen from the definition, $\delta$-PRS is not unique. Indeed, if $\mathcal{R}_{\delta,t}$ is a $\delta$-PRS, then any $\mathcal{R}_{\delta,t}^{\prime}\supseteq\mathcal{R}_{\delta,t}$ is also a $\delta$-PRS. We say $\mathcal{R}_{\delta,t}$ is a tighter $\delta$-PRS than $\mathcal{R}_{\delta,t}^{\prime}$ if $\mathcal{R}_{\delta,t}\subseteq\mathcal{R}_{\delta,t}^{\prime}$. An illustration of DRS and $\delta$-PRS is in Figure 1. In many applications like safety-critical control, it is desirable to have a tight $\delta$-PRS, and a loose $\delta$-PRS can result in very conservative strategies. Therefore, the main goal in probabilistic reachability is to find a tight $\delta$-PRS.

The trajectories of the stochastic system (1) are driven by both deterministic terms ($f$,$u_{t}$) and stochastic disturbances ($w_{t}$). The impacts of these two parts on the trajectories are relatively independent and can be addressed separately. Building on this intuition, we proposed a strategy termed the separation strategy in our recent work [24] for probabilistic reachability analysis for continuous-time systems.

This separation strategy turns out to be extendable to the discrete-time setting. We say $X_{t}$ and $x_{t}$ are associated trajectories if they have the same initial state $x_{0}$ and the same input $u_{t}$. The separation strategy states that the $\delta$-PRS can be separated into a deterministic part encoded by the DRS of (2), and a stochastic part, which is characterized as the distance $\|X_{t}-x_{t}\|$ between the associated trajectories, as formalized below.

Let $X_{t}$ be any trajectory of (1) associated with a trajectory $x_{t}$ of (2), then, by (6), $X_{t}\in\mathcal{B}^{n}(r_{\delta,t},x_{t})$ with probability at least $1-\delta$. By the definition of $\overline{\mathcal{R}}_{t}$, $x_{t}\in\overline{\mathcal{R}}_{t}$. Therefore, by the definition of the Minkowski sum [28], with probability at least $1-\delta$,

which completes the proof. We term the difference $\|X_{t}-x_{t}\|$ between associated trajectories stochastic deviation. This separation strategy decomposes the probabilistic reachability analysis problem into two parts: approximate the DRS of (2) and estimate the probabilistic bound $r_{\delta,t}$ of the stochastic deviation. Once a bound $r_{\delta,t}$ of the stochastic deviation is obtained, one can combine it with any existing deterministic reachability method to approximate the $\delta$-PRS.

It is apparent that the size of the $\delta$-PRS $\overline{\mathcal{R}}_{t}\oplus\mathcal{B}^{n}(r_{\delta,t},0)$ increases with $r_{\delta,t}$. Therefore, to ensure that $\overline{\mathcal{R}}_{t}\oplus\mathcal{B}^{n}(r_{\delta,t},0)$ is not an overly-conservative $\delta$-PRS of (1), it is crucial to establish a probabilistic bound $r_{\delta,t}$ for the stochastic deviation that is as tight as possible. This is the main challenge addressed in this paper.

We employ a similar technique to that in [24, 29] to analyze the stochastic deviation $\|X_{t}-x_{t}\|$ and establish bounds on $\mathbb{E}(\|X_{t}-x_{t}\|^{2})$. The result is presented in the following proposition.

Define $v_{t}=X_{t}-x_{t}$, $\beta_{t}=f(X_{t},u_{t},t)-f(x_{t},u_{t},t)$ and $V(\cdot)=\|\cdot\|^{2}$, then by (1) and (2) we have

It follows that

Take the expectation of (11). By Assumption 1, $\mathbb{E}\left(V(\beta_{t})\right)\leq L_{t}^{2}\mathbb{E}\left(v_{t}\right)$. Since $w_{t}$ is independent from $\beta_{t}$, $\mathbb{E}\left(\langle\beta_{t},w_{t}\rangle\right)=0$. Let $w_{t}=\begin{bmatrix}w_{t,1}&\cdots&w_{t,n}\end{bmatrix}^{\top}$, then by Definition 2.2, each entry $w_{t,i}$ of $w_{t}$ is sub-Gaussian with variance proxy $\sigma_{t}^{2}$. Therefore,

where the last “$\leq$” uses the classical result that $\text{Var}(X)\leq\sigma^{2}$ if $X\in\mathbb{R}$ and $X\sim subG(\sigma^{2})$ [30]. Solving the linear difference inequality (12) we obtain

where $\Psi_{t}$ is defined in (9). This completes the proof.

The probabilistic bound directly derived from the expectation bound (8) can be loose. By Markov’s Inequality, for any $r>0$,

Let $r=\sqrt{\frac{n}{\delta}\Psi_{t}}$ and apply Proposition 2 to (14), then

Given a probability level $1-\delta$, (15) provides a probabilistic bound on stochastic deviation that has $\mathcal{O}(\sqrt{1/\delta})$ dependence on $\delta$. This is much worse than that for linear systems. To see this, consider a linear dynamics $X_{t+1}=AX_{t}+Bu_{t}+w_{t}$ which corresponds to (1) with $f(x,u,t)=Ax+Bu$. It satisfies Assumption 1 with $L_{t}\equiv\|A\|$, for every $t\geq 0$. For associated $X_{t}$ and $x_{t}$, by linearity,

By the sum-up property of independent sub-Gaussian variables [31], (16) implies that $X_{t}-x_{t}$ is sub-Gaussian with variance proxy

Sub-Gaussian random variables enjoy a strong norm concentration property stated below. Note that the dependence $\mathcal{O}(\sqrt{n})$ and $\mathcal{O}(\sqrt{\log(1/\delta)})$ in Lemma 3.1 are tight for sub-Gaussian distributions but the coefficients can be improved. We refer to [31, Chapter 1.4][24] for more details.

Plugging (17) into Lemma 3.1, we conclude that when the system (1) is linear, the probabilistic bound on the stochastic deviation becomes

which has an $\mathcal{O}(\sqrt{\log(1/\delta))}$ dependence on $\delta$. When $\delta$ is sufficiently small (e.g., $\delta=10^{-6}$), which is crucial for safety-critical systems, $\sqrt{\log(1/\delta))}$ is significantly smaller than the term $\sqrt{1/\delta}$ in (15) (e.g., $2.6$ v.s. $10^{3}$).

Thus, a significant gap between the result (15) for nonlinear dynamics and the probabilistic bound (20) for linear dynamics exists. This points to the question: is the gap fundamental or an artifact of the analysis?

The limitation of the expectation bound (8) primarily lies in the quadratic Lyapunov function $V_{t}=\|X_{t}-x_{t}\|^{2}$. The analysis focuses only on the evolution of the second order moment $\mathbb{E}(\|X_{t}-x_{t}\|^{2})$, so its upper bound at best guarantees a probabilistic bound for $\|X_{t}-x_{t}\|$ of order $\mathcal{O}(\sqrt{1/\delta})$ via Markov inequality. In contrast, in Definition 2.2, the function on the left-hand side of (4) captures arbitrarily high-order moments of a random variable, and its boundness leads to the norm concentration property of sub-Gaussian variables. This function is known as Moment Generating Function (MGF) [31, Chapter 1.1]

Thus, to establish a tight probabilistic bound on the stochastic deviation, a potential approach is to upper bound the MGF of $X_{t}-x_{t}$ so that $X_{t}-x_{t}$ is sub-Gaussian and Lemma 3.1 can be applied. Unfortunately, this is infeasible. For associated $X_{t}$ and $x_{t}$, $\mathbb{E}(X_{t})\neq x_{t}$ when the systems are nonlinear, because

where the last “$\neq$” can be improved to “$=$” only if $f$ is linear. Therefore, the upper bound on $\mathbb{E}\left(M_{\lambda,\ell}(X_{t}-x_{t})\right)$ is influenced by the value of $\mathbb{E}\left(\langle\ell,X_{t}-x_{t}\rangle\right)$ [32, Chapter 2.5], which can be large for some $\ell$ when $\mathbb{E}(\|X_{t}-x_{t}\|)$ is large.

Motivated by the advantages and limitations of MGF, we resort to a modified version of MGF termed Averaged Moment Generating Function (AMGF) in our method.

The AMGF was recently proposed in [33] in the field of Sampling Theory. It has a natural interpretation as the average of MGF over the unit sphere $\ell\in\mathcal{S}^{n-1}$. AMGF can also be viewed as a special MGF while replacing the exponential energy function $e^{\lambda\langle\ell,x\rangle}$ by its sphere-averaged version $\Phi_{n,\lambda}(x)=\mathbb{E}_{\ell\sim\mathcal{S}^{n-1}}\left(e^{\lambda\langle\ell,x\rangle}\right)$. This energy function $\Phi_{n,\lambda}(x)$ exhibits several intriguing properties.

Moreover, AMGF induces the same norm concentration property as MGF.

We emphasize that even though (24) is a weaker condition than sub-Gaussian, it still guarantees norm concentration property as MGF. To see why Lemma 4.2 holds, define an intermediate random variable $\tilde{X}=QX$ where $Q\sim\mathbb{U}^{n}$ is a random unitary matrix with $\mathbb{U}^{n}$ denoting the set of all the unitary matrices in $\mathbb{R}^{n\times n}$. Then the AMGF over $X$ is equal to the MGF over $\tilde{X}$, which means $\tilde{X}$ is sub-Gaussian with variance proxy $\sigma^{2}$ if (24) holds. Norm concentration of $X$ then follows by noticing that the transformation $\tilde{X}=QX$ does not affect the norm.

Equipped with the AMGF, we establish a probabilistic bound of order $\mathcal{O}(\sqrt{\log(1/\delta)})$ for the stochastic deviation $\|X_{t}-x_{t}\|$ by bounding the evolution of AMGF $\mathbb{E}(\Phi_{n,\lambda}(X_{t}-x_{t}))$ over time.

We begin with the special case where the Lipschitz constant $L=1$ at any time $t\geq 0$. Define $v_{t}=X_{t}-x_{t}$ and $\beta_{t}=f(X_{t},u_{t},t)-f(x_{t},u_{t},t)$. By (1) and (2) we have

By Assumption 1,

By Assumption 2 and Lemma 4.1(3), the conditional expectation of $\Phi_{n,\lambda}(v_{t+1})$ can be bounded as

Combining (27) and (28), invoking Lemma 4.1(2), we obtain

Taking the expectation over $v_{t}$ on both sides of (29) yields

which is a linear difference inequality that points to

By Lemma 4.2, it follows that for any $\delta\in(0,1)$

holds with probability at least $1-\delta$, where $\varepsilon_{1},\varepsilon_{2}$ are as in (19). This completes the proof of the special case.

Next, we consider the general cases under Assumption 1. Define $\tilde{X}_{t}=\prod_{k=0}^{t-1}L_{k}^{-1}X_{t}$ and $\tilde{x}_{t}=\prod_{k=0}^{t-1}L_{k}^{-1}x_{t}$. Then $\tilde{X}_{t}$ satisfies

and similarly, $\tilde{x}_{t}$ satisfies

which has the same drift term as (33). For any $\tilde{x}_{t},\tilde{y}_{t}\in\mathbb{R}^{n}$, it holds that

meaning the scaled systems (32) and (33) satisfy Assumption 1 with $\tilde{L}=1$. Also, $\tilde{w}_{t}$ satisfies Assumption 2 with variance proxy $\tilde{\sigma}^{2}_{t}=\prod_{k=0}^{t}L_{k}^{-2}\sigma_{t}^{2}$ due to the sum-up property of independent sub-Gaussian noise [31]. Thus the result of the special case can be applied.

By (31), with probability at least $1-\delta$,

Because of the relation $\|X_{t}-x_{t}\|=\prod_{k=0}^{t-1}L_{k}\|\tilde{X}_{t}-\tilde{x}_{t}\|$, the probabilistic bound (25) follows. This completes the proof.

Theorem 1 establishes $r_{\delta,t}=\sqrt{\Psi_{t}(\varepsilon_{1}n+\varepsilon_{2}\log(1/\delta))}$ as a probabilistic bound on stochastic deviation. The bound is comprised of a scaling term $\sqrt{\Psi_{t}}$, a dimensional term $\sqrt{n}$ and a probabilistic term $\sqrt{\log(1/\delta)}$. The scaling term is determined by the Lipschitz constants $L_{t}$ and the disturbance variance proxy $\sigma_{t}^{2}$. It reveals the fluctuation of the system. Especially, when the system is linear, $\Psi_{t}$ is exactly the variance proxy of $X_{t}-x_{t}$. The term $\sqrt{n}$ reveals the dimensional dependence of $r_{\delta,t}$, which is the same as the dimensional dependence of sub-Gaussian vectors shown in Lemma 3.1. The term $\sqrt{\log(1/\delta)}$ essentially represents the fast-decay feature of the distribution of $X_{t}-x_{t}$. Notice that this term means the probability level $\delta$ has an $\mathcal{O}(e^{-r_{\delta,t}^{2}})$ dependence on the probabilistic bound $r_{\delta,t}$, so $\mathbb{P}\left(\|X_{t}-x_{t}\|>r_{\delta,t}\right)$ decreases with $r_{\delta,t}$ even faster than the exponential rate. This phenomenon leads to the following question: Is this the best rate one can obtain? In other words, is $r_{\delta,t}$ derived from Theorem 1 tight?

We next show that the probabilistic bound in Theorem 1 is tight under Assumptions 1 and 2 and it is impossible to achieve better probabilistic bounds than (25) without additional assumptions. In particular, we show that the tightness is precisely captured by linear systems satisfying Assumptions 1 and 2.

Consider the linear stochastic system

that satisfies Assumptions 1 and 2 with $L\equiv\|A\|$ and $\sigma_{t}\equiv\sigma$. Applying the time-invariant $L$ and $\sigma$ into (9), we get that $\Psi_{t}=\frac{\sigma^{2}(L^{2t}-1)}{L^{2}-1}$. Therefore, by Theorem 1, we get that with probability at least $1-\delta$,

This is exactly the same as the $r_{\delta,t}^{(lin)}$ in (20). From [32, Chapter 4.6], we know that at least for linear systems satisfying Assumptions 1 and 2, the dependence on $\delta$ and $n$ can not be improved further without additional assumptions. Therefore, $r_{\delta,t}$ is tight in our problem setting.

The worst-case analysis is a popular method for understanding the reachability when the disturbance is bounded [4]. We point out that it can also be applied to stochastic systems to obtain a probabilistic bound on the stochastic deviation under unbounded stochastic disturbance. As explained below, this is achieved by viewing unbounded sub-Gaussian disturbance as a bounded disturbance with high probability, reminiscent of conformal prediction [34, 35]. However, the result is still significantly more conservative than our bound in Theorem 1.

By norm concentration (Lemma 3.1) of sub-Gaussian random variables, for any $\delta\in(0,1)$,

Applying union bound to (36) over $\tau=0,\dots,t-1$, we obtain

This means that the disturbance is bounded with

for all $\tau<t$ with probability at least $1-\delta$. A worst-case probabilistic bound on $\|X_{t}-x_{t}\|$ can be established by assuming this bound. More specifically, by the system dynamics (1)-(2) and the triangular inequality,

for all $\tau<t$ with probability at least $1-\delta$. It follows that

where $\psi_{t}$ is as in (9). Plugging (38) into (40), we conclude that, for any $\delta\in(0,1)$,

holds with probability at least $1-\delta$.

This bound (41) derived using worst-case analysis is substantially more conservative than that in Theorem 1. First, since $\sqrt{\Psi_{t}}\leq\sqrt{\psi_{t-1}}\sum_{k=0}^{t-1}\sigma_{k}\sqrt{\psi_{k}^{-1}}$, (41) is always worse than (25). To see more clearly the gap, consider the case when $L_{t}\equiv L$ and $\sigma_{t}\equiv\sigma$. In this case, (41) reduces to $\|X_{t}-x_{t}\|\leq\frac{L^{t}-1}{L-1}\sqrt{\sigma^{2}(\varepsilon_{1}n+\varepsilon_{2}\log\frac{t}{\delta})}$, which is much worse than (34) when $L$ is close to 1. In particular, when $L=1$, $\sqrt{\psi_{t}}\sum_{k=1}^{t}\sigma_{k-1}\sqrt{\psi_{t}^{-1}}=\sigma t$, significantly larger than $\sqrt{\Psi_{t}}=\sigma\sqrt{t}$ as $t$ increases.