Underapproximation of Reach-Avoid Sets for Discrete-Time Stochastic Systems via Lagrangian Methods

TL;DR

This paper introduces a grid-free, Lagrangian-based method to efficiently compute conservative underapproximations of stochastic reach-avoid sets for discrete-time nonlinear systems, improving scalability and computational efficiency.

Contribution

It unifies existing Lagrangian approaches, characterizes the disturbance subset for guaranteed underapproximation, and demonstrates significant efficiency gains over traditional gridding methods.

Findings

Method provides scalable, grid-free underapproximations.

Approach is conservative but computationally efficient.

Validated on space vehicle rendezvous and 2D integrator examples.

Abstract

We examine Lagrangian techniques for computing underapproximations of finite-time horizon, stochastic reach-avoid level-sets for discrete-time, nonlinear systems. We use the concept of reachability of a target tube in the control literature to define robust reach-avoid sets which are parameterized by the target set, safe set, and the set in which the disturbance is drawn from. We unify two existing Lagrangian approaches to compute these sets and establish that there exists an optimal control policy of the robust reach-avoid sets which is a Markov policy. Based on these results, we characterize the subset of the disturbance space whose corresponding robust reach-avoid set for the given target and safe set is a guaranteed underapproximation of the stochastic reach-avoid level-set of interest. The proposed approach dramatically improves the computational efficiency for obtaining an underapproximation of stochastic reach-avoid level-sets when compared to the traditional approaches based on gridding. Our method, while conservative, does not rely on a grid, implying scalability as permitted by the known computational geometry constraints. We demonstrate the method on two examples: a simple two-dimensional integrator, and a space vehicle rendezvous-docking problem.

Tables

Table 1. TABLE I : Computation times, in seconds, for double integrator problem with dynamic programming and with Algorithm 1, and the ratio of computation times (dynamic programming by Algorithm 1). Algorithm 1 does not require a grid.

Grid Size	Dynamic Programming	Algorithm 1	Ratio
$41\times 41$	8.16	0.98	8.3
$82\times 82$	59.76	0.98	60.9

Equations

$$x_{k+1}=f(x_k,u_k)+w_k$$

$$x_{k+1}$$

$$D_t(\mathcal{T},\mathcal{K},E)$$

$$\forall k\in {\mathbb{Z}}_{[0,t-1]},x_k\in \mathcal{K},x_t\in \mathcal{T}\}.$$

$$x_{k+1}$$

$$V_t(\mathcal{K})$$

$$=D_t(\mathcal{K},\mathcal{K},\{0_n\}).$$

$$\underset{\pi \in \mathcal{M}}{\sup }{\mathbb{E}}_{\bar{x}}^{N,\pi }\left[\left(\prod _{i=0}^{N-1}{1}_{\mathcal{K}}(x_i)\right)1_{\mathcal{T}}(x_N)\right].$$

$$V_k^{\ast }(x)$$

$$=1_{\mathcal{K}}(x){\mathbb{P}}_{{\bar{x}}_k}^{N-k,{\pi }^{\ast }}(x_N\in \mathcal{T},x_{N-1}\in \mathcal{K},\dots ,$$

$$x_{k+1}\in \mathcal{K}|x)$$

$$V_N^{\ast }(x)$$

$${\mathcal{L}}_k(\beta )=\left\{x\in \mathcal{X}:V_{N-k}^{\ast }(x)\ge \beta \right\},$$

$${\mathcal{F}}_1(x)$$

$${\mathcal{R}}_1(\mathcal{S})$$

$$=\left\{x^{-}\in \mathcal{X}:{\mathcal{F}}_1(x^{-})\cap \mathcal{S}\ne \varnothing \right\}$$

$${\mathcal{F}}_1(x)$$

$${\mathcal{R}}_1(\mathcal{S})$$

$$D_0(\mathcal{T},\mathcal{K},E)$$

$$D_k(\mathcal{T},\mathcal{K},E)$$

$${\mathcal{F}}_1(x_0)\cap (D_{k-1}(\mathcal{T},\mathcal{K})\ominus E)\ne \varnothing \}.$$

$$D_1(\mathcal{T},\mathcal{K},E)$$

$$\exists x^{+}\in \mathcal{T},x^{+}=f(x_0,{\nu }_0(x_0))+w_0\}$$

$$=\{x_0\in \mathcal{X}:x_0\in \mathcal{K},\exists u\in \mathcal{U},\exists {\nu }_0(\cdot )\in \mathcal{F},$$

$$\exists y\in (\mathcal{T}\ominus E),y=f(x_0,{\nu }_0(x_0)),u={\nu }_0(x_0)\}$$

$$=\left\{x_0\in \mathcal{K}:\exists y,y\in {\mathcal{F}}_1(x_0)\wedge y\in (\mathcal{T}\ominus E)\right\}$$

$$=\left\{x_0\in \mathcal{K}:{\mathcal{F}}_1(x_0)\cap (D_0(\mathcal{T},\mathcal{K},E)\ominus E)\ne \varnothing \right\}$$

$$D_{t-1}(\mathcal{T},\mathcal{K},E)$$

$$\forall k\in {\mathbb{Z}}_{[0,t-2]},x_k\in \mathcal{K},x_{t-1}\in \mathcal{T}\}.$$

$$D_t(\mathcal{T},\mathcal{K},E)$$

$$=\{x_0\in \mathcal{X}:\exists {\rho }_t\in {\mathcal{P}}_t,\forall {\bar{w}}_t\in E^t,\forall k\in {\mathbb{Z}}_{[0,t-1]},$$

$$x_k\in \mathcal{K},x_t\in \mathcal{T}\}$$

$$=\{x_0\in \mathcal{X}:x_0\in \mathcal{K},\exists {\nu }_0(\cdot )\in \mathcal{F},\exists {\rho }_{t-1}\in {\mathcal{P}}_{t-1},$$

$${\rho }_t=[{\nu }_0,{\rho }_{t-1}],\forall w_0\in E,\forall {\bar{w}}_{t-1}\in E^{t-1},$$

$$\forall k\in {\mathbb{Z}}_{[1,t-1]},x_k\in \mathcal{K},x_t\in \mathcal{T}\}$$

$$=\{x_0\in \mathcal{K}:\exists {\nu }_0(\cdot )\in \mathcal{F},\forall w_0\in E,\exists x_0^{+}\in D_{t-1}(\mathcal{T},\mathcal{K},E)$$

$$f(x_0,{\nu }_0(x_0))+w_0=x_0^{+}\}$$

$$=\{x_0\in \mathcal{K}:\exists {\nu }_0(\cdot )\in \mathcal{F},\exists y\in D_{t-1}(\mathcal{T},\mathcal{K})\ominus E,$$

$$y=f(x_0,{\nu }_0(x_0))\}$$

$$=\left\{x_0\in \mathcal{K}:{\mathcal{F}}_1(x_0)\cap (D_{t-1}(\mathcal{T},\mathcal{K},E)\ominus E)\ne \varnothing \right\}.$$

$$D_0(\mathcal{T},\mathcal{K},E)$$

$$D_k(\mathcal{T},\mathcal{K},E)$$

$${\mathcal{S}}_2\ominus \{0_n\}$$

Full text

I Introduction

Reach-avoid analysis is an established verification tool that provides formal guarantees of both safety (via avoiding unsafe regions) and performance (via reaching a target set). It has been used in safety-critical or expensive systems, for example, with application to space systems [1], aviation [2, 3], biomedical systems [4], and other domains [5, 6, 7]. The reach-avoid set is the set of initial states for which there exists control that enables the state to reach a target within some finite time horizon, while remaining within a safe set (avoiding an unsafe set) for all instants in the time horizon. In a probabilistic system, satisfaction of the reach-avoid objective is accomplished stochastically. The stochastic reach-avoid level-set for a given likelihood is the set of states for which probabilistic success of the reach-avoid objective is assured with at least the given likelihood.

The theoretical framework for the probabilistic reach-avoid calculation is based on dynamic programming [8, 7], and, hence, is computationally infeasible for even moderate-sized systems due to the gridding of not only the state-space, but also of the input and disturbance spaces [9]. Recent work has focused on alternatives to dynamic programming, including approximate dynamic programming [10, 11, 6], Gaussian mixtures [11], particle filters [1, 6], and convex chance-constrained optimization [1, 5]. These methods have been applied to systems that are at most 10-dimensional, at high memory and computational costs [6]. Further, since an analytical expression of the value function is not accessible, stochastic reach-avoid level-sets can be computed only up to the accuracy of the gridding.

We propose a method to compute an underapproximation of probabilistic reach-avoid sets via robust reach-avoid sets, the set of states assured to reach the target set and remain in the safe region despite any disturbance input. Robust reach-avoid sets can be theoretically posed as the solution to the reachability of a target tube problem [12, 13, 14], originally framed to compute reachable sets of discrete-time controlled systems with bounded disturbance sets. Motivated by the scalability of the Lagrangian method proposed in [4, 15] for viability analysis in deterministic systems (that is, systems without a disturbance input but with a control input), we seek a similar approach to compute the robust reach-avoid sets via tractable set theoretic operations. Lagrangian methods rely on computational geometry, whose scalability depends on the representation and the operation used [16], including polyhedrons (implementable using Model Parametric Toolbox (MPT) [17]), support functions [18], and ellipsoids (implementable via the Ellipsoidal Toolbox [19]).

In this paper, we unify these two approaches to create an efficient algorithm for underapproximation of the stochastic reach-avoid set, and demonstrate our approach on practical examples. Our main contributions are: a) synthesis of the approaches presented in [4, 15] and [12, 13, 14] to compute the robust reach-avoid sets, b) sufficient conditions under which an optimal control policy for a given robust reach-avoid set is a Markov policy, and c) an algorithm to compute an underapproximation of the stochastic reach-avoid level-sets using the robust reach-avoid sets. Specifically, we establish the sufficient conditions under which an optimal control policy is comprised of universally measurable state-feedback laws. For these conditions, we characterize the subset of the disturbance space whose corresponding robust reach-avoid set is a guaranteed underapproximation of the desired stochastic reach-avoid level-set. Leveraging established Lagrangian methods, we demonstrate that our approach dramatically reduces the computation time required for computing a conservative underapproximation of the desired stochastic reach-avoid level-set. Further, the Lagrangian methods does not rely on grids, freeing the underapproximated sets from any numerical artifacts arising due to the discretization.

The remainder of the paper is as follows: Section II describes the problem and the necessary notation. In Section III, we describe the relationship between the the recursion established in [12] for the robust reach-avoid set and the Lagrangian approach in [4], and establish the desired measurability properties of the optimal controller. We present an algorithm for underapproximation of stochastic reach-avoid level-sets in Section IV. We demonstrate our algorithm on two examples—a simple two-dimensional integrator and a space vehicle rendezvous-docking problem—in Section V and provide conclusions and directions of future work in Section VI.

II Problem Statement

The following notation will be used throughout the paper: we denote discrete-time time intervals by $\mathbb{Z}_{[a,b]}=\mathbb{Z}\cap\{a,a+1,\dots,b-1,b\}$ for $a,b\in\mathbb{Z}$; the set of natural numbers (including zero) as $\mathbb{N}$; the Minkowski sum of two sets $\mathcal{S}_{1},\mathcal{S}_{2}$ as $\mathcal{S}_{1}\oplus\mathcal{S}_{2}=\{s_{1}+s_{2}:s_{1}\in\mathcal{S}_{1},s_{2}\in\mathcal{S}_{2}\}$; the Minkowski difference (or Pontryagin difference) of two sets $\mathcal{S}_{1},\mathcal{S}_{2}$ as $\mathcal{S}_{2}\operatorname{\ominus}\mathcal{S}_{1}=\{s:s+s_{1}\in\mathcal{S}_{2}\ \forall s_{1}\in\mathcal{S}_{1}\}$; and the indicator function corresponding to a set $\mathcal{S}$ as $\boldsymbol{1}_{\mathcal{S}}:\mathcal{X}\rightarrow\{0,1\}$ where $\boldsymbol{1}_{\mathcal{S}}(x)=1$ if $x\in\mathcal{S}$ and is zero otherwise.

II-A System formulation

We consider a discrete-time, nonlinear, time-invariant system with an affine disturbance,

[TABLE]

with state $x_{k}\in\mathcal{X}\subseteq\mathbb{R}^{n}$, input $u_{k}\in\mathcal{U}\subseteq\mathbb{R}^{m}$, disturbance $w_{k}\in\mathcal{W}\subseteq\mathbb{R}^{n}$, and a function $f:\mathcal{X}\times\mathcal{U}\rightarrow\mathcal{X}$. Without loss of generality, we assume $\mathcal{W}$ contains $0_{n}$, the zero vector of $\mathbb{R}^{n}$. We will also consider the discrete, LTI system of form

[TABLE]

for some matrices $A\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{n\times m}$. We assume $A$ is non-singular, which holds true for discrete-time systems that arise from the discretization of continuous-time systems.

II-B Robust reach-avoid sets

Let $t\in\mathbb{N}$ and $\mathcal{F}$ denote the set of admissible state-feedback laws, $\nu:\mathcal{X}\rightarrow\mathcal{U}$. We define a control policy as a sequence of state-feedback laws, $\rho_{t}=[\nu_{0}(\cdot),\dots,\nu_{t-1}(\cdot)]$ with $\nu_{k}\in\mathcal{F}$ for $k\in\mathbb{Z}_{[0,t-1]}$. We denote the corresponding set of admissible control policies as $\mathcal{P}_{t}$.

Let $\operatorname{\mathcal{E}}\subseteq\mathcal{W}$ be a subset of the disturbance set. We define the $t$-time robust reach-avoid set corresponding to $\operatorname{\mathcal{E}}$ as the set of initial states $x_{0}\in\mathcal{X}$ such that there exists an admissible control policy $\rho_{t}\in\mathcal{P}_{t}$ that ensures $x_{k}$ remains in a safe set $\mathcal{K}\subseteq\mathcal{X}$ for $k\in\mathbb{Z}_{[0,t-1]}$ and $x_{t}$ lies in a target set $\mathcal{T}\subseteq\mathcal{X}$ (reach-avoid objective) despite the presence of the disturbance $w_{k}\in\operatorname{\mathcal{E}}$ at each instant. Denoting $\bar{w}_{t}={[w_{0}^{\top},\dots,w_{t-1}^{\top}]}^{\top}\in\mathcal{W}^{t}$, the $t$-time robust reach-avoid set is

[TABLE]

Note that for $\operatorname{\mathcal{E}}=\{0_{n}\}$, the system (1) is equivalent to a deterministic, discrete-time, nonlinear system

[TABLE]

when $w_{k}\in\operatorname{\mathcal{E}}$. The $t$-time viable set of the system (4) is the set of initial states $x_{0}\in\mathcal{X}$ such that there exists an admissible control policy $\rho_{t}\in\mathcal{P}_{t}$ such that $x_{k}$ remains in a safe set $\mathcal{K}$ for $k\in\mathbb{Z}_{[0,t]}$. That is,

[TABLE]

The authors in [4] presented a Lagrangian formulation to compute $\operatorname{\mathcal{V}}_{t}(\mathcal{K})$ and discussed the scalability of the viability analysis using MPT, ET, and support functions.

II-C Stochastic reach-avoid level-sets

In this subsection, we further assume the disturbance $w_{k}$ in (1) is an $n$-dimensional random vector defined in the probability space $(\mathcal{W},\sigma(\mathcal{W}),\mathbb{P}_{w})$. Here, $\sigma(\mathcal{W})$ denotes the minimal $\sigma$-algebra associated with the random vector $w_{k}$. We assume the disturbance $w_{k}$ is absolutely continuous with a probability density function (PDF) $\psi_{w}$, the disturbance process ${\{w_{k}\}}_{k=0}^{N-1}$ is an independent and identically distributed (i.i.d.) random process, and $N\in\mathbb{N}$ is a finite time horizon. We assume that $f$ is Borel-measurable, $\mathcal{U}$ is compact, the sets $\mathcal{K},\mathcal{T}$ are Borel, and $\psi_{w}$ is continuous.

We denote the set of universally measurable state-feedback laws $\mu(\cdot):\mathcal{X}\rightarrow\mathcal{U}$ as $\mathcal{F}_{u}$. We define the Markov control policy as $\pi=[\mu_{0}(\cdot),\dots,\mu_{N-1}(\cdot)]$ where $\mu_{k}\in\mathcal{F}_{u}\ \forall k\in\mathbb{Z}_{[0,N-1]}$, and $\mathcal{M}$ is the set of admissible Markov policies. Since no measurability restrictions were imposed on the feedback laws in Section II-B, $\mathcal{F}_{u}\subseteq\mathcal{F}$ and $\mathcal{M}\subseteq\mathcal{P}_{N}$.

Given a Markov policy $\pi$ and initial state $x_{0}\in\mathcal{X}$, the concatenated state vector $\bar{x}=[x_{1},\dots,x_{N}]$ for the system (1) is a random vector defined in the probability space $(\mathcal{X}^{N},\sigma(\mathcal{X}^{N}),\mathbb{P}^{N,\pi}_{\bar{x}})$. The probability measure $\mathbb{P}^{N,\pi}_{\bar{x}}$ is induced from the probability measure $\mathbb{P}_{w}$ via (1) [7]. We will denote the probability space associated with the random vector $\bar{x}_{k}=[x_{k+1},\dots,x_{N}]$ as $(\mathcal{X}^{N-k},\sigma(\mathcal{X}^{N-k}),\mathbb{P}^{N-k,\pi}_{\bar{x}_{k}})$ for $k\in\mathbb{Z}_{[0,N-1]}$.

For stochastic reachability analysis, we are interested in the maximum likelihood that the system (1) starting at an initial state $x_{0}\in\mathcal{X}$ will achieve the reach-avoid objective using a Markov policy. The maximum likelihood and the optimal Markov policy can be determined as the solution to the optimization problem, [7]

[TABLE]

A dynamic programming approach was presented in [7] to solve problem (7). Let the optimal solution to problem (7) be $\pi^{\ast}=[\mu_{0}^{\ast}(\cdot)\ \ldots\ \mu_{N-1}^{\ast}(\cdot)]$, the maximal Markov policy in the terminal sense [7, Def. 10]. The existence of a Markov policy is guaranteed for a continuous $\psi_{w}$ and compact $\mathcal{U}$ [20, Thm. 1]. The approach in [7] generates value functions $V_{k}^{\ast}:\mathcal{X}\rightarrow[0,1]$ for $k\in[0,N]$,

[TABLE]

initialized with

[TABLE]

By definition, the optimal value function $V_{0}^{\ast}(x_{0})$ provides the maximum likelihood, optimal value of problem (7), of achieving the reach-avoid objective by the system (1) for the time horizon $N$ and the initial state $x_{0}\in\mathcal{X}$.

For $\beta\in[0,1]$ and $k\in\mathbb{Z}_{[0,N]}$, the stochastic reach-avoid $\beta$-level-set,

[TABLE]

is the set of states $x$ that achieve the reach-avoid objective by the time horizon with a probability of, at minimum, $\beta$, in the time interval $\mathbb{Z}_{[0,k]}$.

II-D Problem statements

The following problems are addressed in this paper:

Problem 1**.**

Construct a recursion for exact computation of the robust reach-avoid sets (3) for the system (1).

Problem 2**.**

Given a set $\operatorname{\mathcal{E}}\subseteq\mathcal{W}$ and the corresponding robust reach-avoid set (3), characterize the sufficient conditions under which there exists an optimal control policy that is a Markov control policy for the system (1).

Problem 3**.**

Given $\beta\in[0,1]$, characterize $\operatorname{\mathcal{E}}\subseteq\mathcal{W}$ whose corresponding robust reach-avoid set (3) underapproximates the stochastic reach-avoid $\beta$-level-set (10).

Problem 3a**.**

For a given $\beta\in[0,1]$, characterize an algorithm to compute $\operatorname{\mathcal{E}}$ for Problem 3 when the disturbance in (1) is Gaussian.

III Robust Reach-Avoid Set Computation

In this section, we characterize the robust reach-avoid set for the system described in (1). To solve Problem 1, we first extend the approach presented in [4, 15] to reproduce the results presented in [12]. The authors in [4] demonstrated scalability of the Lagrangian methods for viability analysis in deterministic systems. By unifying these approaches, we aim for a tractable and efficient Lagrangian computation of the robust reach-avoid set with established scalability properties. We also demonstrate that the recursion presented for the viable set computation in deterministic system [4] is a special case of the proposed Lagrangian approach. Finally, we solve Problem 2 and establish that there is an optimal control policy for the robust reach-avoid set that is also a Markov policy.

III-A Iterative computation for robust reach-avoid sets

Similar to the work in [4], for the system (1), we define the unperturbed, one-step forward reach set from a point $x\in\mathcal{X}$ as $\operatorname{\mathscr{F}_{1}}(x)$, and the unperturbed, one-step backward reach set from a set $\mathcal{S}\subseteq\mathcal{X}$ as $\operatorname{\mathscr{R}_{1}}(\mathcal{S})$. Formally, for the system (1),

[TABLE]

where (12) follows from (11). For the system (2),

[TABLE]

Proposition 1**.**

Given a set $\operatorname{\mathcal{E}}\subseteq\mathcal{W}$, the finite horizon robust reach-avoid sets for the system (1) can be computed recursively as follows for $k\geq 1,k\in\mathbb{N}$:

[TABLE]

*Proof: * We first show the case $k=1$, which differs slightly from other cases. From (1) and (3),

[TABLE]

For any $t\in\mathbb{N},t>1$, from (3),

[TABLE]

Using (17), we construct $\operatorname{\mathcal{D}}_{t}(\mathcal{T},\mathcal{K},\operatorname{\mathcal{E}})$ in the form of (16).

[TABLE]

Since the choice of $w_{0}$ depends only $(x_{0},\nu(x_{0}))$, the terms $\exists\rho_{t-1}$ and $\forall w_{0}$ can be exchanged in (18). We obtain (19) after exchanging the terms and applying (17).

Theorem 1**.**

For the system given in (1), the finite-time robust reach-avoid sets $\operatorname{\mathcal{D}}_{k}$ can be computed using the recursion for $k\geq 1,k\in\mathbb{N}$:

[TABLE]

*Proof: * Follows from Proposition 1 and (12).

Figure 1 depicts the recursion in Theorem 1 (21) graphically. From (21), we have the following corollary.

Corollary 1**.**

$\operatorname{\mathcal{D}}_{k}(\mathcal{T},\mathcal{K},\operatorname{\mathcal{E}})\subseteq K\ \forall k\geq 1,k\in\mathbb{N}$.

For completeness, we establish that the viability analysis presented in [4] is a special case of Theorem 1. From [21, Theorem 2.1], for any $\mathcal{S}_{1},\mathcal{S}_{2}\subseteq\mathcal{X}$, $\mathcal{S}_{2}\operatorname{\ominus}\mathcal{S}_{1}=\cap_{s\in\mathcal{S}_{1}}(\mathcal{S}_{2}\oplus\{-s\})$. Hence,

[TABLE]

Corollary 2**.**

[4, Theorem 1]** The finite horizon viable sets for (4) can be computed recursively as follows:

[TABLE]

*Proof: * Follows from Theorem 1, (4), and (22).

A similar recursion can be provided for computing the reach-avoid sets for a deterministic system.

Lemma 1**.**

[12, Proposition 3]** For the dynamics (2), if $\mathcal{U},\mathcal{K},\mathcal{T}$ are convex and compact sets, $\operatorname{\mathcal{E}}$ is a compact set, and $A$ in the dynamics (2) is non-singular, then $\operatorname{\mathcal{D}}_{k}(\mathcal{T},\mathcal{K},\operatorname{\mathcal{E}})$ is convex and compact $\forall k\in\mathbb{N}$.

Note that convexity of $\operatorname{\mathcal{D}}_{k}(\mathcal{T},\mathcal{K},\operatorname{\mathcal{E}})$ does not require convexity of $\operatorname{\mathcal{E}}$. Further, for polyhedral $\mathcal{U},\mathcal{K},\mathcal{T}$, the robust reach-avoid set $\operatorname{\mathcal{D}}_{k}(\mathcal{T},\mathcal{K},\operatorname{\mathcal{E}})$ is polyhedral for $k\in\mathbb{N}$. Note that the same can not be said be for ellipsoids [12, Sec. 4]. A detailed discussion for the implementation of Theorem 1 for polyhedral sets using support functions is given in [12, App. A].

III-B Minmax problem for robust reach-avoid set computation

We will now address Problem 2. A minmax optimization problem was presented in [12, Sec. 1], [22, Sec. 4.6.2] to compute the robust reach-avoid sets (3) for the system (1). The optimization problem is:

[TABLE]

where the decision variables are $\rho_{t}$ and $\bar{w}_{t}$. Here, $g_{k}(\cdot)=1-\boldsymbol{1}_{\mathcal{K}}(\cdot)$ for $k\in\mathbb{Z}_{[0,t-1]}$ and $g_{t}(\cdot)=1-\boldsymbol{1}_{\mathcal{T}}(\cdot)$. The objective function $J(\cdot)$ is parameterized by the initial state $x_{0}\in\mathcal{X}$ and the time horizon $t\in\mathbb{N}$. Problem (29) can be solved using dynamic programming [22, Sec 1.6] to generate the value functions $J^{\ast}_{k}(x;t):\mathcal{X}\rightarrow\mathbb{Z}_{[0,t-k+1]}$ for $k\in\mathbb{Z}_{[0,t]}$

[TABLE]

initialized with $J^{\ast}_{t}(x;t)=g_{t}(x)$. The optimal value of problem (29) when starting at $x_{0}$ is $J^{\ast}_{0}(x_{0};t)$. Further,

[TABLE]

Recall that lower semi-continuous functions are functions whose sublevel-sets are closed and upper-semicontinuous functions are functions whose negative is a lower semi-continuous function [23, Definition 7.13]. Also, the supremum of a lower-semicontinuous function is the negative of the infimum of an upper-semicontinuous function. Let an optimal control policy for problem (29) be $\rho_{t}^{\ast}=[\nu_{0}^{\ast}(\cdot)\ \ldots\ \nu_{t-1}^{\ast}(\cdot)]$. Note that $\rho_{t}^{\ast}$ need not be unique.

Theorem 2**.**

For closed sets $\mathcal{K},\mathcal{T}$ and compact set $\mathcal{U}$, there exists an optimal policy $\rho^{\ast}_{t}$ for problem (29) which is also a Markov policy.

*Proof: * We show by induction that the statement S: the optimal value functions $J_{k}^{\ast}$ of (29) are lower-semicontinuous and there exists a Borel-measurable state-feedback law $\nu_{k}^{\ast}(\cdot)$ for every $k\in\mathbb{Z}_{[0,N-1]}$. Since Borel-measurability implies universal measurability [23, Definition 7.20], the proof of Theorem 2 follows from S and the definition of a Markov policy.

Proof of S: The closedness property of $\mathcal{K},\mathcal{T}$ imply $g_{k}(\cdot)$ is lower semi-continuous for $k\in\mathbb{Z}_{[0,t]}$. Hence, $J_{t}^{\ast}(x;t)$ is lower-semicontinuous.

Consider the base case $k=t-1$. From [23, Prop. 7.32 (b)], we can see that $H_{t-1}^{\ast}(u,x;t)$ is lower semi-continuous. From [23, Prop. 7.33], we conclude that $J_{t-1}^{\ast}(x;t)$ is lower semi-continuous and an optimal state-feedback policy $\nu_{t-1}^{\ast}(\cdot)$ exists which is also Borel-measurable.

Let $\tau\in\mathbb{Z}_{[1,t-2]}$. Assume, for induction, the case $k=\tau$ is true, i.e, $J_{\tau}^{\ast}$ is lower semi-continuous. The proof that $J_{\tau-1}^{\ast}$ is lower semi-continuous and the existence of a Borel-measurable $\nu_{\tau-1}^{\ast}(\cdot)$ follows from [23, Prop. 7.32(b) and 7.33]. This completes the induction.

IV Conservative Approximation of Stochastic Reach-Avoid Level-Set

We will now focus on the stochastic system described in Section II-C and use the theory developed in Section III to solve Problems 3 and 3a.

Theorem 3**.**

Given closed sets $\mathcal{K},\mathcal{T}$, compact set $\mathcal{U}$ and a set $\operatorname{\mathcal{E}}\subseteq\mathcal{W}$. For every $x\in\operatorname{\mathcal{D}}_{t}(\mathcal{T},\mathcal{K},\operatorname{\mathcal{E}})$ with $t\in\mathbb{Z}_{[1,N]}$,

[TABLE]

*Proof: * Follows from Theorem 2 and the definition of $\operatorname{\mathcal{D}}_{t}(\mathcal{T},\mathcal{K},\operatorname{\mathcal{E}})$ (3).

Theorem 4**.**

Given $\beta\in[0,1]$, closed sets $\mathcal{K},\mathcal{T}$, and a compact set $\mathcal{U}$, if for any $t\in\mathbb{Z}_{[0,N]}$, $\operatorname{\mathcal{E}}\subseteq\mathcal{W}$ such that $\mathbb{P}_{w}(w_{k}\in\operatorname{\mathcal{E}})=\beta^{\frac{1}{t}}$ for all $k\in\mathbb{Z}_{[0,t-1]}$, then $\operatorname{\mathcal{D}}_{t}(\mathcal{T},\mathcal{K},\operatorname{\mathcal{E}})\subseteq\mathcal{L}_{t}(\beta)$ .

*Proof: * The case for $t=0$ follows trivially from (9), (10), and (15). Let $t>0$ and $x\in\operatorname{\mathcal{D}}_{t}(\mathcal{T},\mathcal{K},\operatorname{\mathcal{E}})$. We are interested in underapproximating $\mathcal{L}_{t}(\beta)=\{x:V_{N-t}^{\ast}(x)\geq\beta\}$ as defined in (10). From (8),

[TABLE]

Equation (34) follows from the law of total probability and Corollary 1 which implies $\boldsymbol{1}_{\mathcal{K}}(x)=1$. Equation (35) follows from (34) after ignoring the second term (which is non-negative). Simplifying (35) using Theorem 3 and the i.i.d. assumption of the disturbance process, we obtain

[TABLE]

Thus, $\operatorname{\mathcal{D}}_{t}(\mathcal{T},\mathcal{K},\operatorname{\mathcal{E}})\subseteq\mathcal{L}_{t}(\beta)$ by (10).

Theorem 4 solves Problem 3 for an arbitrary density $\psi_{w}$. Computation of $\operatorname{\mathcal{D}}_{t}(\mathcal{T},\mathcal{K},\operatorname{\mathcal{E}})$ can be done via Theorem 1. Note that $\operatorname{\mathcal{E}}$ characterized by Theorem 4 is not unique. Recall that Corollary 2 states that the robust reach-avoid set is exact viable set for the case when $\operatorname{\mathcal{E}}=\{0_{n}\}$. We therefore prescribe $\operatorname{\mathcal{E}}$ that contains $\{0_{n}\}$ and has the least Lebesgue measure to reduce the the degree of conservativeness in Theorem 4. We also recommend the set $\operatorname{\mathcal{E}}$ be convex and compact for computational ease.

Next, we provide a method to compute $\operatorname{\mathcal{E}}\subseteq\mathcal{W}$ for any $t\in\mathbb{Z}_{[0,N-1]}$ such that $\mathbb{P}_{w}(w_{k}\in\operatorname{\mathcal{E}})=\beta^{\frac{1}{t}}$ for all $k\in\mathbb{Z}_{[0,t-1]}$ when the disturbance in the system (1) is a Gaussian random vector.

IV-A Computation of $\operatorname{\mathcal{E}}$ for Gaussian disturbance

Let the disturbance in (1) be $w_{k}=v$, an $n$-dimensional Gaussian random variable with mean vector $\mu$ and covariance matrix $\Sigma$. The probability density of a multivariate Gaussian random vector is [24, Ch. 29]

[TABLE]

Consider the $n$-dimensional ellipsoid parameterized by $R^{2}\in[0,\infty)$ for $\operatorname{\mathcal{E}}$,

[TABLE]

For $\mu=0,\Sigma=r^{2}I_{n}$, we have $\operatorname{\mathcal{E}}_{R^{2}}=\{w:w^{\top}w\leq r^{2}R^{2}\}$, an $n$-dimensional hypersphere of radius $rR$. We aim to compute the parameter $R^{2}$ such that $\mathbb{P}_{v}\{v\in\operatorname{\mathcal{E}}_{R^{2}}\}=\beta^{\frac{1}{t}}$ for application of Theorem 4.

Given a standard normal distributed $n$-dimensional random vector $\eta\sim N(0,I_{n})$, $v={\Sigma}^{\frac{1}{2}}\eta+\mu$ [24, Ch. 29]. Also, $\operatorname{\mathcal{E}}_{R^{2}}={\Sigma}^{\frac{1}{2}}\operatorname{\mathcal{E}}_{\eta,R^{2}}\oplus\{\mu\}$ with $\operatorname{\mathcal{E}}_{\eta,R^{2}}=\{s\in\mathbb{R}^{n}:s^{\top}s\leq R^{2}\}$. Since the affine transformation of $\eta$ to $v$ is deterministic, $\mathbb{P}_{v}\{v\in\operatorname{\mathcal{E}}_{R^{2}}\}=\mathbb{P}_{\eta}\{\eta\in\operatorname{\mathcal{E}}_{\eta,R^{2}}\}=\beta^{\frac{1}{t}}$. From [24, Ex. 20.16], we have

[TABLE]

where $\chi^{2}(n)$ is a chi-squared random variable with $n$ degrees of freedom and $F_{\chi^{2}(n)}(\cdot)$ denotes its cumulative distribution function. Consequently, we have

[TABLE]

Equations (37) and (38) solves Problem 3a.

IV-B Computing the stochastic level-set underapproximation

A pseudo-algorithm to compute the underapproximation of the $N$-time stochastic reach-avoid $\beta$-level-set is shown in Algorithm 1 using robust reach-avoid sets. Note that while the system dynamics permitted for Algorithm 1 is the nonlinear system given in (1), the computation of $\operatorname{\mathscr{R}_{1}}(S)$ is accessible only for linear system (2) as defined in (14). Further, Lemma 1 guarantees convexity and compactness of the robust reach-avoid sets, allowing for easy representation, only for linear system dynamics.

Since these sets are formed using Lagrangian techniques, Algorithm 1 is more computationally efficient than the dynamic programming based discretization approach. Algorithm 1 requires a number of basic geometric operations. We will focus on the implementation of Algorithm 1 in a polyhedral representation using the readily available MATLAB toolbox MPT. From Lemma 1, we note that support function methods can also be used [4, 12]. The conservativeness of the underapproximations obtained using Algorithm 1 are very problem dependent. The system dynamics, strength of the disturbance process, and size of the targe and safe sets can all have non-trivial affects on the resulting conservativeness.

The robust reach-avoid set computation requires a Minkowski difference operation as well as an intersection operation in the recursion (21). Hence both facet and vertex representations are typically required for polytopes and numerical implementations will be limited by the well-known vertex-facet enumeration problem. Support functions would not be subject to this problem but require analytic solutions to support vector calculations. Minkowski differences can be handled for polytopes using the MPT toolbox [17] but implementation using ellipses is not feasible without further underapproximation. Additional problems such as redundancy in vertices and facets also commonly arise using polytope representations.

V Examples

All results were obtained using the MPT toolbox [17] with MATLAB R2016a running on Windows 7 computer with and Intel Core i7-2600 CPU, 3.6 GHz, and 8 GB RAM. We focus on examples in which clear comparisons of conservativeness can be made since the ability to handle high-dimensional systems is established in [4].

V-A 2-Dimensional Double Integrator

The first example considered is the stochastic viability analysis of a 2-dimensional discrete-time double integrator model. This example can be solved with both the proposed Lagrangian methods as well as with dynamic programming, allowing for direct comparisons of conservativeness and speed.

The discretized double integrator dynamics are

[TABLE]

The state $x_{k}\in\mathcal{X}\subseteq\mathbb{R}^{2}$, input $u_{k}\in\mathcal{U}\subseteq\mathbb{R}$, $T=0.25$, and the disturbance is assumed to be i.i.d. Gaussian $w_{k}\sim N(0,0.005I_{2})$.

Figure 2 compares the underapproximation via Algorithm 1 and the level-sets computed using dynamic programming techniques, as in [7].The underapproximation is closest and the approximation become progressively more conservative as $N$ increases, as is expected. For $N_{1},N_{2}\in\mathbb{N}$, $N_{2}>N_{1}$, $\beta^{\frac{1}{N_{2}}}>\beta^{\frac{1}{N_{1}}}$, and hence, from Section IV-A, $R_{N_{2}}>R_{N_{1}}$, indicating that $\operatorname{\mathcal{E}}_{N_{1}}\subset\operatorname{\mathcal{E}}_{N_{2}}$.

A comparison between the total computation time for the dynamic programming method and Algorithm 1 is provided in Table I. The accuracy of dynamic programming relies on its grid size, resulting in a trade-off between accuracy and computation speed, from which Algorithm 1 does not suffer.

For systems with Gaussian disturbance processes that have a very low variance, the underapproximation obtained through the Lagrangian methods tightly approximates the stochastic level-set and is computed significantly faster—over $7$ times faster for a $41\times 41$ grid. Figure 3 shows a comparison of the stochastic level-set and the Lagrangian underapproximation when the Gaussian disturbance process is of the form $w_{k}\sim N(0,10^{-5}I_{2})$. The bumps on the exterior of the stochastic level-set are a numerical artifact from the state-space gridding.

V-B Application to space-vehicle dynamics

In this section, we compute an underapproximation of the stochastic reach-avoid level-set for a spacecraft rendezvous docking problem using Algorithm 1. The goal is for a spacecraft, referred to as the deputy, to approach and dock to an orbiting satellite, referred to as the chief, while remaining in a predefined line-of-sight cone. The dynamics are described by the Clohessy-Wiltshire-Hill (CWH) equations [25]

[TABLE]

The chief is located at the origin, the position of the deputy is at $x,y\in\mathbb{R}$, $\omega=\sqrt{\mu/R_{0}^{3}}$ is the orbital frequency, $\mu$ is the gravitational constant, and $R_{0}$ is the orbital radius of the spacecraft.

We define the state vector $z=[x,y,\dot{x},\dot{y}]\in\mathbb{R}^{4}$ and input vector $u=[F_{x},F_{y}]\in\mathcal{U}\subseteq\mathbb{R}^{2}$. We discretize the dynamics (40) in time to obtain the discrete-time LTI system,

[TABLE]

where $w_{k}\in\mathbb{R}^{4}$ is assumed to be a Gaussian i.i.d. disturbance with $\mathbb{E}[w_{k}]=0$, $\mathbb{E}[w_{k}w_{k}^{\top}]=10^{-4}\times\mbox{diag}(1,1,0.0005,0.0005)$.

We define the target set and the constraint set as in [1]

[TABLE]

Figure 4 shows a cross-section at $\dot{x}=\dot{y}=0$ of the resulting underapproximation of the stochastic reach-avoid level set. The computation time for the $N=5$ level-set was 14.5 seconds. Because of the extreme computational requirements of solving a 4-dimensional problem via dynamic programming methods we cannot make a direct computational comparison between dynamic programming and Algorithm 1. In [1, Figure 2], a cross-section of $\dot{x}=\dot{y}=0.9$ of the stochastic reach-avoid set was approximated using convex, chance-constrained optimization and particle approximation methods. Since both of these methods require gridding, the computation time is slower, reported to be about 20 minutes, for a subset of the state space.

VI Conclusion

In this work, we provide a Lagrangian method to compute an underapproximation of a stochastic reach-avoid level-set using robust reach-avoid sets. We synthesize approaches in [4, 15] and [12, 13, 14], and characterize the sufficient conditions under which a optimal control policy for the robust reach-avoid set is also a Markov policy. We demonstrate that our Lagrangian approach to compute the underapproximation is significantly faster when compared to the dynamic programming approach. The utility of this method is problem-dependent, as the conservativeness of the underapproximations are affected by the system dynamics and noise processes.

In future, we intend to examine methods to reduce the conservativeness of the underapproximation and extend the computation of the disturbance set $\operatorname{\mathcal{E}}$ in Theorem 4 for disturbances other than a Gaussian random vector.