Under-Approximate Reachability Analysis for a Class of Linear Systems with Inputs

Reachable sets and tubes of dynamical systems are central in control synthesis and verification applications, especially in the presence of uncertainties and constraints [1, 2]. Mere approximations of reachable sets and tubes are not sufficient in such frameworks. Instead, conservative estimations, i.e., over (outer)-approximations, are typically utilized to ensure all possible behaviors of a given control system are accounted for, which explains the sheer number of over-approximation methods in the literature [3, 4].

In the last few years, there has been a growing interest in additionally under-approximating reachable sets and tubes for synthesis and verification (see, e.g., [5, 6, 7, 8]), because under-approximations can be used in estimating subsets of the states that are attainable under given control constraints [9], obtaining subsets of the initial states from which all trajectories fulfill safety and reachability specifications [10], solving falsification problems: verifying if reachable sets/tubes intersect with unsafe sets [11], and measuring the accuracy of computed over-approximations. Motivated by the aforementioned applications, we aim in this paper at investigating under-approximations of forward reachable sets and tubes for continuous-time LTI systems with uncertainties or constraints on the initial and input values, where we attempt to overcome some of the limitations associated with available methods in the literature.

Optimal control-based polytopic approaches [12, 13] were proposed for linear systems with uncertain inputs and initial conditions. These methods rely on obtaining boundary points of reachable sets, associated with specified direction vectors, and then computing the convex hull of the obtained boundary points. Given a reachable set, the convergence of the polytopic approaches requires computing an increasing number of boundary points until the whole reachable set boundary is obtained, which is computationally expensive, especially if the dimension of the reachable set is high. An ellipsoidal method was proposed in [14] for controllable linear systems with ellipsoidal initial and input sets. The method relies on solving initial value problems, derived from maximal principles similar to those presented in [13], to obtain ellipsoidal subsets that touch a given reachable set at some boundary points, depending on specified direction vectors. The accuracy of the ellipsoidal method of [14] in under-approximating a given reachable set is increased by evaluating an increasing number of ellipsoids, which necessitates considering an increasing number of direction vectors, and then taking their union. The aforementioned optimal control-based approaches [12, 13, 14] assume the ability to compute transition matrices and their integrals exactly, and that is not feasible in general. In addition, when under-approximating a reachable tube, the mentioned approaches use non-convex representations for the under-approximations, which are challenging to analyze in the contexts of verification and synthesis (see the introduction in [15]).

In [16], a set-propagation technique was proposed, yielding convergent under-approximations of forward reachable sets and tubes for a general class of linear systems with uncertainties or constraints on the initial and input values, where the under-approximations of reachable sets and tubes are given as convex sets and finite unions of convex sets, respectively, which can be analysed with relative ease. However, the method in [16], like the approaches in [17, 12, 13, 14], assumes the ability to compute transition matrices and their integrals exactly. A similar set-propagation approach was proposed in [18] for LTI systems, which relies on computing Minkowski differences to under-approximate reachable tubes. The method in [18] also suffers from the issue of the method in [16], in addition to the computational hurdle of evaluating Minkowski differences (see, e.g., [19]).

The issue of evaluating transition matrices and their integrals exactly can be solved by adopting formally correct under-approximation methods that are designed for general nonlinear systems. For example, a formally correct interval arithmetic under-approximation approach was proposed in [8]; however, such a method lacks convergence guarantees, and may produce empty under-approximations. In [20], a novel method was proposed for nonlinear systems with uncertain initial conditions that depends on computing over-approximations and then scaling them down to obtain under-approximations. A drawback of the approach in [20] is that the scaling necessitates solving (sub-optimal) optimization problems that involve enclosures of boundaries of reachable sets, which can be computationally expensive (see the comparison in Section VII-C1). Finally, a recent approach has been proposed in [21] to under-approximate reachable sets when the system parameters (e.g. system matrices) are not known exactly , where collected trajectories (i.e., data) are utilized to estimate the system dynamics. Such an approach is highly valuable in applications when system identification cannot be attained, due, e.g., to failure or damage mid-operation. However, this approach is conservative (i.e., convergence cannot be attained in general) as it considers the set of all systems that can generate the collected trajectories, while fulfilling some specified assumptions.

In this work, we present a novel efficient approach that results in under-approximations of forward finite-time reachable sets and tubes for a class of LTI systems with inputs, where approximations of the matrix exponential and its integral are used, first-order convergence guarantees are provided, and approximations of reachable sets and tubes are given as convex sets and finite unions of convex sets, respectively. Our approach is fundamentally based on set-based recursive relations (see, e.g., [22, 16]), where truncation errors are accounted for in an under-approximating manner, utilizing matrix lower bounds [23].

Let $\mathbb{R}$, $\mathbb{R}_{+}$, $\mathbb{Z}$, $\mathbb{Z}_{+}$, and $\mathbb{C}$ denote the sets of real numbers, non-negative real numbers, integers, non-negative integers, and complex numbers, respectively, and $\mathbb{N}\mathrel{:=}\mathbb{Z}_{+}\setminus\{0\}$. Let $[a,b]$, $\orbrack a,b\clbrack$, $[a,b\clbrack$, and $\orbrack a,b]$ denote closed, open, and half-open intervals, respectively, with end points $a$ and $b$, and $[a;b]$, $\orbrack a;b\clbrack$, $[a;b\clbrack$, and $\orbrack a;b]$ stand for their discrete counterparts, e.g., $[a;b]=[a,b]\cap\mathbb{Z}$, and $[1;4\clbrack=\{1,2,3\}$. Given any map $f\colon A\to B$, the image of a subset $C\subseteq A$ under $f$ is given by $f(C)=\left\{f(c)\,\middle|\,c\in C\right\}$. The identity map $X\to X\colon x\mapsto x$ is denoted by $\operatorname{id}$, where the domain of definition $X$ will always be clear form the context. The Minkowski sum of $M,N\subseteq\mathbb{R}^{n}$ is defined as $M+N\mathrel{:=}\left\{y+z\,\middle|\,y\in M,z\in N\right\}$. By $\|\cdot\|$ we denote any norm on $\mathbb{R}^{n}$, the norm of a non-empty subset $M\subseteq\mathbb{R}^{n}$ is defined by $\|M\|\mathrel{:=}\sup_{x\in M}\|x\|$, $\mathbf{B}_{n}\subseteq\mathbb{R}^{n}$ is the closed unit ball w.r.t. $\|\cdot\|$, and the maximum norm on $\mathbb{R}^{n}$ is denoted by $\|\cdot\|_{\infty}$ ($\|x\|_{\infty}=\max\left\{|x_{i}|\,\middle|\,i\in[1;n]\right\},~{}x\in\mathbb{R}^{n}$). Given norms on $\mathbb{R}^{n}$ and $\mathbb{R}^{m}$, $\mathbb{R}^{n\times m}$ is endowed with the usual matrix norm, $\|A\|=\sup_{\|x\|\leq 1}\|Ax\|$ for $A\in\mathbb{R}^{n\times m}$, e.g., the matrix norm of $A$ induced by the maximum norm is $\|A\|_{\infty}=\max_{i\in[1;n]}\sum_{j=1}^{m}\left|A_{i,j}\right|$. Given a square matrix $A\in\mathbb{R}^{n\times n}$, $\rho(A)$ denotes the spectral radius of $A$, i.e., $\rho(A)\mathrel{:=}\max\{\left|\lambda\right|,~{}\lambda\textrm{ is an eigenvalue of }A\}$. The spectral radius satisfies the property below, which follows from [24, proof of Lemma 5.6.10, p. 348]).

Given $A\in\mathbb{R}^{n\times m}$, $\mathrm{rank}(A)$, $\mathrm{col}(A)$, and $A^{\dagger}\in\mathbb{R}^{m\times n}$ denote the rank, the column space, and the Moore–Penrose inverse of $A$, respectively ($A^{\dagger}=A^{-1}$ if $A$ is invertible). The following lemma provides a sufficient condition to check invertibility.

Given $A\in\mathbb{R}^{n\times m}\setminus\{0\}$, $\|A\|_{l}$ denotes the matrix lower bound of $A$ w.r.t. $\|\cdot\|$, which is defined as $\|A\|_{l}\mathrel{:=}\max\{m\in\mathbb{R}~{}|~{}\forall y\in\mathrm{col}(A),\exists x\in\mathbb{R}^{m}~{}\mathrm{s.t.}~{}Ax=y~{}\mathrm{and}~{}m\|x\|\leq\|y\|\},$ see [23]. Matrix lower bounds satisfy the following properties, which are essential in our derivation of the proposed method.

The collection of full-dimensional subsets of $\mathbb{R}^{n}$ that are affine transformations of closed unit balls is denoted by $\mathbb{A}_{n}$, where, in this work, saying $\Omega=c+G\mathbf{B}_{p}\in\mathbb{A}_{n}$ implies that $c\in\mathbb{R}^{n}$, $G\in\mathbb{R}^{n\times p}$, and $\mathrm{rank}(G)=n$. Integration of single-valued functions presented herein is always understood in the sense of Bochner. Given a non-empty subset $X\subseteq\mathbb{R}^{n}$ and a measurable subset $S\subseteq\mathbb{R}$, $X^{S}$ denotes the set of Lebesgue measurable maps with domain $S$ and values in $X$. Given a non-empty compact subset $W\subseteq\mathbb{R}^{m}$, a compact interval $[a,b]\subset\mathbb{R}$, and an integrable matrix-valued function $F:[a,b]\to\mathbb{R}^{n\times m}$, we define the set-valued integral $\int_{a}^{b}F(t)W\mathrm{d}t\mathrel{:=}\bigcup_{w\in W^{[a,b]}}\int_{a}^{b}F(t)w(t)\mathrm{d}t.$ The Hausdorff distance $\mathfrak{d}(\Omega,\Gamma)$ of two non-empty bounded subsets $\Omega,\Gamma\subseteq\mathbb{R}^{n}$ w.r.t. $\|\cdot\|$ is defined as $\mathfrak{d}(\Omega,\Gamma)\mathrel{:=}\inf\left\{\varepsilon>0\,\middle|\,\Omega\subseteq\Gamma+\varepsilon\mathbf{B}_{n},\Gamma\subseteq\Omega+\varepsilon\mathbf{B}_{n}\right\}.$ The Hausdorff distance satisfies the triangle inequality, in addition to the following set of properties.

In this paper, we consider the LTI system

over the time interval $[0,T]$, where $x(t)\in\mathbb{R}^{n}$ is the system state, $u(t)\in\mathbb{R}^{n}$ is the input, and $A\in\mathbb{R}^{n\times n}$ is the system matrix. The initial value $x(0)$ and the input $u(t),~{}t\in[0,T]$ belong to known sets $X_{0}$ and $U$, respectively. Let $T$, $A$, $X_{0}$, and $U$ be fixed and assume that:

1. 1.

   The time interval $[0,T]$ is compact and $T>0$.

2. 2.

   $X_{0}=c_{x}+G_{x}\mathbf{B}_{p_{x}}\in\mathbb{A}_{n}$ and $U=c_{u}+G_{u}\mathbf{B}_{p_{u}}\in\mathbb{A}_{n}$, where $c_{x},c_{u},G_{x}$, and $G_{u}$ are known.

Given an initial value $x(0)=x_{0}$ and an integrable input signal $v:[0,T]\rightarrow\mathbb{R}^{n}$, the unique solution, $\varphi(\cdot,x_{0},v)$, to system (1), generated by $x_{0}$ and $v(\cdot)$, on $[0,T]$ is given by [26, Theorem 6.5.1, p. 114]

Herein, $\mathrm{e}^{(\cdot)A}$ (or $\exp((\cdot)A)$) is the matrix exponential function, which has the Taylor series expansion $\exp(tA)=\sum_{j=0}^{\infty}{(tA)^{j}}/{j!}.$ Define

where $\mathcal{L}(t,k)$ is the truncated $(k-1)$th-order Taylor expansion of $\exp(tA)$ and $\mathcal{T}(t,k)$ is its definite integral. It is easy to verify that, for all $t\in\mathbb{R}_{+}$ and $k\in\mathbb{N}$,

where $\theta\colon\mathbb{R}_{+}\times\mathbb{N}\rightarrow\mathbb{R}_{+}$ is defined as

The function $\theta$ is infinitely differentiable and monotonically increasing in its first argument, and monotonically decreasing in its second argument, with a greatest lower bound of zero.

Let $\mathcal{R}(t)$ denote the forward reachable set of system $\eqref{eq:LinearSystem}$ at time $t\in[0,T]$, with initial values in $X_{0}$, and input signals with values in $U$. In other words,

The set $\exp(tA)X_{0}$ is referred to as the homogeneous reachable set at time $t$ and is denoted by $\mathcal{R}_{h}(t)$, and the set $\int_{0}^{t}\exp(sA)U\mathrm{d}s$ is referred to as the input reachable set at time $t$ and is denoted by $\mathcal{R}_{u}(t)$. Furthermore, let $[a,b]\subseteq[0,T]$. Then, $\mathcal{R}([a,b])=\bigcup_{t\in[a,b]}\mathcal{R}(t)$ is the reachable tube over the time interval $[a,b]$. In this paper, we aim to compute arbitrarily precise under-approximations of $\mathcal{R}_{h}(T)$, $\mathcal{R}_{u}(T)$, $\mathcal{R}(T)$, and $\mathcal{R}([0,T])$, utilizing the approximations $\mathcal{L}$ and $\mathcal{T}$.

In this section, we thoroughly derive the proposed method. The convergence guarantees of the method are discussed in Section VI.

We start with the following theoretical recursive relation, which is algorithmically similar to efficient over-approximation methods in the literature [9], and is the basis of the proposed method of this work.

The algorithm in Lemma 5, in theory, addresses the under-approximation problem of this work; however, this algorithm cannot be implemented exactly in general. In the next sections, we address the challenges in implementing this theoretical algorithm and propose a practically implementable method, which is the main contribution of this work.