Inner-approximate Reachability Computation via Zonotopic Boundary Analysis

The notations and operations concerning space, vectors, matrices, and sets utilized in this paper are presented in Table 1, where the symbols and descriptions for operations on vectors, matrices, and sets are mainly illustrated with specific examples of a vector $\bm{x}$, a matrix $\bm{M}$, and a set $\Delta$.

This paper considers nonlinear systems which are modelled by ordinary differential equations of the following form:

where $\bm{x}\in\mathbb{R}^{n}$ and $\bm{f}$ is a locally Lipschitz continuous function. Thus, given an initial state $\bm{x}_{0}$, there exists an unique solution $\phi(\cdot;\bm{x}_{0}):[0,T_{\bm{x}_{0}})\rightarrow\mathbb{R}^{n}$ to system (1), where $[0,T_{\bm{x}_{0}})$ is the maximal time interval on which $\phi(\cdot;\bm{x}_{0})$ is defined.

Given a set $\mathcal{X}_{0}$ of initial states, the reachable set is defined as follows:

The exact reachable set $\Phi(t;\mathcal{X}_{0})$ is usually impossible to be computed, especially for nonlinear systems. Outer-approximations and inner-approximations are often computed for formal reasoning on the system.

In this paper, we focus on the computation of an inner-approximation represented by zonotopes. Zonotope is a special class of convex polytopes with the centrally symmetric nature. It can be viewed as a Minkowski sum of a finite set of line segments, known as G-representation, which is defined as the following.

For a zonotope $Z=\left\langle\bm{c},\bm{G}\right\rangle$ in space $\mathbb{R}^{n}$, it is called $k$-dimensional if $\operatorname*{rank}(\bm{G})=k,k\leq n$. A $k$-dimensional zonotope can be reduced into space $\mathbb{R}^{k}$ without altering its shape. Furthermore, the facets of a $k$-dimensional ($k\geq 1$) zonotope are $(k-1)$-dimensional zonotopes. If an $n$-dimensional zonotope has $n$ independent generators, then it’s called parallelotope. Additionally, If there is a zonotope $Z^{\prime}$ such that $Z^{\prime}\subsetneqq Z$, $Z^{\prime}$ is called a sub-zonotope of $Z$.

The framework of computing inner-approximations in this paper follows the one proposed in [35], but with minor modifications.

Given system (1) with an initial set $\mathcal{X}_{0}$, represented by a zonotope, and a time duration $T=Nh$, where $h>0$ is the time step and $N$ is a non-negative integer, we compute a zonotopic inner-approximation $U((k+1)h;\mathcal{X}_{0})$ of the reachable set $\Phi((k+1)h;\mathcal{X}_{0})$ for $k\in\{0,1,\cdots,N\}$. The inner-approximation $U_{k+1}=U((k+1)h;\mathcal{X}_{0})$ is computed based on $U_{k}=U(kh;\mathcal{X}_{0})$ ($U_{0}:=\mathcal{X}_{0})$ with the following procedures:

1. 1.

   extract and refine the boundary $\partial U_{k}$ of $U_{k}$ ;

2. 2.

   compute a zonotopic outer-approximation $O(h;U_{k})$ of reachable set $\Phi(h;U_{k})$, and an outer-approximation $O(h;\partial U_{k})$ of reachable set $\Phi(h;\partial U_{k})$. These outer-approximations can be computed using existing zonotope-based approaches such as [3].

3. 3.

   contract $O(h;U_{k})$ to obtain a zonotopic inner-approximation candidate $U^{\prime}_{k+1}$ by excluding the set $O(h;\partial U_{k})$, i.e., let $U^{\prime}_{k+1}\cap O(h;\partial U_{k})=\emptyset$.

4. 4.

   compute an outer-approximation of the reachable set $O(h;\bm{c})$ of the time-inverted system $\dot{\bm{x}}=-\bm{f}(\bm{x})$ with the single initial state $\bm{c}$, where $\bm{c}$ is the center of the zonotope $U^{\prime}_{k+1}$. If the computed outer-approximation $O(h;\bm{c})$ is included in the set $U_{k}$, then $U_{k+1}:=U^{\prime}_{k+1}$ is an inner-approximation of the reachable set $\Phi((k+1)h;\mathcal{X}_{0})$.

The overall computational workflow is visualized in Fig. 1. There are three computational procedures that affect the efficacy (i.e., accuracy and efficiency) of inner-approximation computation in the aforementioned framework: the extraction and refinement of the boundary $\partial U_{k}$, reachability analysis for computing outer-approximations $O(h;U_{k}),O(h;\partial U_{k})$, and contraction of $O(h;U_{k})$ to obtain an inner-approximation candidate $U_{k+1}^{\prime}$. Since there are well-developed reachability algorithms in existing literature for computing outer-approximations such as [12, 3], we in the following focus on other two computational procedures. For the first one, as the outer-approximation computed $O(h;\partial U_{k})$ would be excluded from $O(h;U_{k})$, the accuracy of $O(h;\partial U_{k})$ significantly affects the one of $U_{k+1}$. Additionally, the accuracy of $O(h;\partial U_{k})$ strongly correlates with the size of $\partial U_{k}$. To improve the accuracy of $U_{k+1}$, two algorithms are proposed: one for extracting and the other for tiling the boundary of a zonotope (i.e., splitting the boundary into sub-zonotopes without overlaps). As for the third one, an adaptive strategy is developed to make the inner-approximation $U_{k+1}$ much tighter. This is achieved by contracting $O(h;U_{k})$ in a flexible way, deviating from the proportional reduction of the size of $O(h;U_{k})$ in the existing methods [35].

In this subsection we introduce the algorithm for extracting the exact boundary of a zonotope. The concept of cross product of a matrix provided by [26] will be utilized herein, which is formulated below.

The cross product of $\bm{M}\in\mathbb{R}^{n,n-1}$ can be viewed as the normal vector of the hyperplane spanned by $n-1$ linearly independent column vectors in $\bm{M}$.

The boundary extraction algorithm is established on the fact that a zonotope is centrally symmetric and each facet, which is a zonotope, has an congruent facet on the opposite side of the center (e.g., two dark blue facets in Fig. 2).

Given an $n$-dimensional zonotope $Z=\left\langle\bm{c},\bm{G}\right\rangle$, where $\bm{c}\in\mathbb{R}^{n}$ and $\bm{G}=(\bm{g}_{i})_{1\leq i\leq p}\in\mathbb{R}^{n,p}$, for each two symmetric facets, they lie in parallel hyperplanes and share the same generators. The two parallel hyperplanes are spanned by a part of generators of $Z$, which can form a submatrix of $\bm{G}$ with rank $n-1$. In boundary extraction algorithm, i.e., Alg. 1, we firstly enumerate all potential $n\times(n-1)$ submatrices of $\bm{G}$ which are able to span a hyperplane. For a certain hyperplane spanned by a submatrix $\bm{B}_{b}$, to confirm the center and generators of its corresponding facets, we compute its normal vector by the cross product operator $\mathtt{CP}(\cdot)$, then the center of the two symmetric facets can be respectively determined by moving the center $\bm{c}$ along the positive and negative directions of generators which are not perpendicular to $\mathtt{CP}(\bm{B}_{b})$, and the generator matrix of these corresponding facets can be represented by $\bm{B}_{b}$ appending generators parallel to the hyperplane. The visible operations stated above are shown in Fig. 2.

The computation of a zonotope’s boundary is summarized in Alg. 1. Its soundness, i.e., the set computed by Alg. 1 is equal to the boundary $\partial Z$ of the zonotope $Z$, is justified in Theorem 3.1, whose proof is available in Appendix A. In order to enhance the understanding of Alg. 1, we provide a simple example, Example 1 in Appendix B, to illustrate the computational process of Alg. 1.

This subsection introduces our tiling algorithm which can split a zonotope into sub-zonotopes without overlaps and then elaborates how this tiling algorithm is employed to refine the boundaries of zonotopes.

The boundary matrix, which is constructed according to Alg. 1, plays an important role in our tiling algorithm. Its entries are able to characterize the centers and generators for all facets of a zonotope.

From the boundary matrix of a zonotope, one can obtain all its facets. Appendix B provides an example (Example 2) to illustrate this claim.

Another matrix, tiling matrix, is constructed to store the outcomes of the tiling algorithm, i.e., all the non-overlapping sub-zonotopes whose union covers the original zonotope. Similar to the boundary matrix, a row of tiling matrix represents a sub-zonotope.

Our tiling algorithm is based an intuitive observation: for a zonotope, moving its one-sided facets towards to the opposite side along the direction of a generator results in a new zonotope with this generator removed, simultaneously, several sub-zonotopes are generated by adding this generator to all these facets. This process, which is visualized in Fig. 3, can be iteratively conducted, until a parallelotope remains. At this point, the tiling algorithm terminates, yielding a collection of tiles denoted as zonotopes that tile the original zonotope.

The tiling algorithm leverages operations on boundary matrix $\bm{B}$ to implement the facets’ movement and sub-zonotopes generation aforementioned. The results of each step, namely the sub-zonotopes after one-step tiling, are recorded in the tiling matrix $\bm{T}$.

Given an $n$-dimensional zonotope $Z=\left\langle\bm{c},\bm{G}\right\rangle,$ where $\bm{G}=(\bm{g}_{i})_{1\leq i\leq p}\in\mathbb{R}^{n,p}$, we require that the right-most $n\times n$ submatrix of $\bm{G}$ is full rank to ensure that the sub-zonotopes with one generator removed after one-step tiling remain $n$-dimensional. For the specific $j$-th column of boundary matrix $\bm{B}$, where $1\leq j\leq(p-n)$, we process the following operations to its entries:

1. 1.

   if there exist $i$’s such that $\bm{B}(i,j)=-1$, we add these rows in the boundary matrix $\bm{B}$ into the tiling matrix $\bm{T}$ as new rows, but change their $j$-th entry to $0$ in the tiling matrix $\bm{T}$. Meanwhile, the $j$-th entries of these rows in the boundary matrix $\bm{B}$ are modified into $1$, i.e., $\bm{B}(i,j)=1$;

2. 2.

   if there exist $i$’s such that $\bm{B}(i,j)=0$, we delete these rows from the boundary matrix $\bm{B}$.

After the $j$-th iteration, the updated boundary matrix $\bm{B}$ characterizes the boundary of a new zonotope. This new zonotope is derived by removing the first through the $j$-th generators from the original zonotope $Z$. Simultaneously, the sub-zonotopes generated by adding the generator $\bm{g}_{j}$ to the facets are incorporated into the tiling matrix $\bm{T}$. Finally, after $p-n$ iterations, there remains one parallelotope, whose generator matrix is the right-most $n\times n$ submatrix of $\bm{G}$, we put this parallelotope into the tiling matrix $\bm{T}$ and then output the result.

The above computational procedures are summarized in Alg. 2. Its soundness is justified by Theorem 3.2, whose proof is available in Appendix A. Moreover, Appendix B supplements an example (Example 3) to illustrate the main steps tiling a zonotope using Alg. 2.

In this subsection we present our contraction method, yielding the inner approximation candidate $U_{k+1}^{\prime}$ by contracting $O(h;U_{k})$. In contrast to the approaches in [35], which contracts $O(h;U_{k})$ by reducing size proportionally, our contraction method offers a more flexible way. Specifically, the length of each generator of $O(h;U_{k})$ can be adjusted and some generators can be removed. The incorporation of this adaptive contraction method enhances the tightness of the computed inner-approximation.

By extracting and refining of boundary $\partial U_{k}$ of $U_{k}$, we get a collection of sub-zonotopes, i.e., $\{\partial U_{k}^{(i)}\}_{i\in\mathbb{N}_{[1,s]}}$, where $\bigcup_{i=1}^{s}\partial U_{k}^{(i)}=\partial U_{k}$. Then $O(h;\partial U_{k})$ can be obtained by uniting all the out-approximations $\partial O_{k+1}^{(i)}:=O(h;\partial U_{k}^{(i)}),i\in\mathbb{N}_{[1,s]}$, i.e., $O(h;\partial U_{k})=\bigcup_{i=1}^{s}\partial O_{k+1}^{(i)}$.

Noticing that the shape of every outer-approximation $\partial O_{k+1}^{(i)}=\left\langle\bm{c}_{o},\bm{G}_{o}\right\rangle$ is usually long and narrow (refer to Fig. 1), we choose the top $n-1$ independent generators by norm (such as Euclidean norm) to span a hyperplane, which can be seen as an $(n-1)$-dimensional form approximating $\partial O_{k+1}^{(i)}$. Then we compute the cross product $\mathtt{CP}(\cdot)$ of this hyperplane as its normal vector to represent the attitude of $\partial O_{k+1}^{(i)}$, denoted by $\mathtt{AT}\left(\partial O_{k+1}^{(i)}\right)$ (i.e. $\mathtt{CP}(\hat{\bm{G}}_{o})$, where $\hat{\bm{G}}_{o}$ contains top $n-1$ indenpent generators of $\bm{G}_{o}$ by norm).

Initially, we set the inner-approximation candidate $U_{k+1}^{\prime}:=O(h;U_{k})$. Subsequently, we iteratively reduce the length of generators and adjust the position (by changing the center) of $U_{k+1}^{\prime}$ until the intersections between $U_{k+1}^{\prime}$ and all outer-approximations $\partial O_{k+1}^{(i)}$ become empty sets. For each outer-approximation $\partial O_{k+1}^{(i)}$, we begin by shortening the length of generators that are most likely to yield collisions between $U_{k+1}^{\prime}$ and $\partial O_{k+1}^{(i)}$, which would prevent the unnecessary contraction of $U_{k+1}^{\prime}$ and make the result tighter. Heuristically, the generators with directions closest to $\mathtt{AT}\left(\partial O_{k+1}^{(i)}\right)$, or in other words, those most likely to “perpendicular” to $\partial O_{k+1}^{(i)}$ (precisely, perpendicular to hyperplane spanned by column vectors of $\hat{\bm{G}}_{o}$) should be given priority considerations. When encountering a generator that dose not need to be shortened, indicating that $U_{k+1}^{\prime}$ and $\partial O_{k+1}^{(i)}$ have no overlapping parts, we turn to the next outer-approximation $\partial O_{k+1}^{(i+1)}$. The details of contraction method proposed is summarized below.

1. 1.

   Initialize inner-approximation candidate $U_{k+1}^{\prime}:=\left\langle\bm{c}_{u},\bm{G}_{u}\right\rangle=O(h;U_{k})$.

2. 2.

   For every boundary outer-approximations $\partial O_{k+1}^{(i)}=\left\langle\bm{c}_{o},\bm{G}_{o}\right\rangle,i\in\mathbb{N}_{[1,s]}$, carry out the following processing steps.

   • 2a.

     Sort the generators $\{\bm{g}_{l}\}_{1\leq l\leq\operatorname*{cols}(\bm{G}_{u})}$ of $U_{k+1}^{\prime}$ according the angle with $\mathtt{AT}\left(\partial O_{k+1}^{(i)}\right)$ from small to large (i.e., $\|{\rm cos}\theta\|=\frac{\|\bm{g}_{l}\cdot\mathtt{AT}\left(\partial O_{k+1}^{(i)}\right)\|}{\|\bm{g}_{l}\|\|\mathtt{AT}\left(\partial O_{k+1}^{(i)}\right)\|}$ from large to small).

   • 2b.

     Loop all the generators according to the sorted order, for the generator $\bm{g}_{l}$, compute its domain $[\underline{\alpha_{l}},\overline{\alpha_{l}}]$ which intersects with $\partial O_{k+1}^{(i)}$ by LPs (2) and (3) (using approach in [17, Chapter $4.2.5$]), where $\bm{\alpha}=(\alpha_{1},\cdots,\alpha_{l},\cdots,\alpha_{\operatorname*{cols}(\bm{G}_{u})})^{\intercal}$, $\bm{\beta}=(\beta_{1},\cdots,\beta_{\operatorname*{cols}(\bm{G}_{o})})^{\intercal}$.

     	$\displaystyle\min$	$\displaystyle\alpha_{l}$		(2)
     	$\displaystyle s.t.$	$\displaystyle\bm{c}_{u}+\bm{G}_{u}\bm{\alpha}=\bm{c}_{o}+\bm{G}_{o}\bm{\beta}$
     		$\displaystyle-\bm{1}\leq\bm{\alpha}\leq\bm{1},-\bm{1}\leq\bm{\beta}\leq\bm{1}$

     	$\displaystyle\max$	$\displaystyle\alpha_{l}$		(3)
     	$\displaystyle s.t.$	$\displaystyle\bm{c}_{u}+\bm{G}_{u}\bm{\alpha}=\bm{c}_{o}+\bm{G}_{o}\bm{\beta}$
     		$\displaystyle-\bm{1}\leq\bm{\alpha}\leq\bm{1},-\bm{1}\leq\bm{\beta}\leq\bm{1}$

     When the optimal value of (2) or (3) can’t be found, then terminate this loop and continue for the next boundary outer-approximation $\partial O_{k+1}^{(i+1)}$.

   • 2c.

     If $[\underline{\alpha_{l}},\overline{\alpha_{l}}]=[-1,1]$, then delete $\bm{g}_{l}$ from generator matrix $\bm{G}_{u}$. Else, update the range of $a_{l}\in\max\{[-1,\underline{a_{l}}-\epsilon],[\overline{a_{l}}+\epsilon,1]\}\triangleq[\underline{\gamma},\overline{\gamma}]$, where the operation $\max\{\cdot,\cdot\}$ means choosing the interval with maximum length and $\epsilon$ is a user-defined small positive number.

   • 2d.

     Update $\bm{c}_{u}:=\bm{c}_{u}+0.5(\overline{\gamma}+\underline{\gamma})\bm{g}_{l}$ and $\bm{g}_{l}:=0.5(\overline{\gamma}-\underline{\gamma})\bm{g}_{l}$.

Appendix B provides an example (Example 4) to illustrate the procedure of the contraction method and why is necessary to sort the generators $\{\bm{g}_{l}\}_{1\leq l\leq\operatorname*{cols}(\bm{G}_{u})}$ of $U_{k+1}^{\prime}$ according to the angle with $\mathtt{AT}\left(\partial O_{k+1}^{(i)}\right)$.