Reachability analysis computes the union of all possible trajectories that a system can reach within a finite or infinite time when starting from a bounded set of initial states . The most popular approaches for computing reachable sets are set-propagation and simulation-based techniques. Set-propagation techniques propagate reachable sets for consecutive time steps. The efficiency depends on the set representation and the used technique. For instance, Taylor series methods propagate enclosures over discrete time by constructing a Taylor expansion of the states with respect to time and bounding the coefficients . The resulting enclosure is then inflated by a bound on the truncation error. Other set representations are polyhedra , zonotopes , (sparse) polynomial zonotopes , ellipsoids , and support functions . Zonotopes have favorable properties as they can be represented compactly, and they are closed under Minkowski sum and linear mapping. The simulation-based approach in  over-approximates the reachable set by a collection of tubes around trajectories such that the union of these tubes provides an over-approximation of the reachable set. Sampling-based reachability analysis utilizes random set theory in . Other simulation-based techniques are proposed in [23, 18, 17, 4, 35].
While there is a considerable amount of literature on computing reachable sets for different model classes, these methods assume an a priori given model. Obtaining a model that adequately describes the system at hand from first principles or noisy data is usually a challenging and time-consuming task. Simultaneously, system data in the form of measured trajectories is often readily available in many applications. Therefore, we are interested in reachability analysis directly from noisy data of an unknown system model. One recent contribution in this direction can be found in 
, where the authors introduce two data-driven methods for computing the reachable sets with probabilistic guarantees. The first method represents the reachability problem as a binary classification problem, using a Gaussian process classifier. The second method makes use of a Monte Carlo sampling approach to compute the reachable set. A probabilistic reachability analysis is proposed for general nonlinear systems using level sets of Christoffel functions in where they provide a guarantee that the output of the algorithm is an accurate reachable set approximation in a probabilistic sense. Over-approximating reachable sets from data are considered in  based on interval Taylor-based methods applied to systems with dynamics described as differential inclusions; however, the proposed approach only works under the assumption of noise-free data. Another interesting method is introduced in , where the model is assumed to be partially known and data is used to learn an additional Lipschitz continuous state-dependent uncertainty, where the unknown part is assumed to be bounded by a known set. We believe that computing guaranteed reachable sets from noisy data is still an open problem.
The main idea underlying the introduced data-driven reachability framework is visualized in Fig. 1 where we compute data-driven reachable sets based on matrix zonotope recursion. In order to guarantee that the reachable set encloses all system trajectories from finite noisy data, we compute a matrix zonotope that encloses all models that are consistent with the noisy data. We then propagate the initial set forward, utilizing this matrix zonotope to compute the reachable set. We also provide an improved reachability analysis algorithm that provides a less conservative over-approximation of the reachable set at the cost of additional computational expenses. This improved analysis scheme is enabled by introducing a new set representation, namely constrained matrix zonotopes. We utilize this novel set representation and the corresponding operations additionally to provide a systematic approach on how supplementary prior knowledge about the unknown model like states decoupling, partial model knowledge, or given bounds on certain entries in the system matrices can be incorporated into the reachability analysis. We then extend our approach to two classes of nonlinear systems: polynomial systems and Lipschitz systems. All used codes to recreate our findings are publicly available111https://github.com/aalanwar/Data-Driven-Reachability-Analysis.
We will specifically build upon ideas used in [45, 25, 24] and [12, 5, 6, 44] for data-driven analysis and data-driven controller design, respectively. In these works, data is used to characterize all models that are consistent with the data. This characterization enables a computational approach for direct systems analysis and design without explicitly identifying a model.
The main contributions of this paper are as follows:
A new set representation named constrained matrix zonotope (Definition 5) and its essential operations are proposed.
A general framework is introduced for incorporating side information like states decoupling, partial model knowledge, or given bounds on certain entries in the system matrices about the unknown model (Algorithm 3) into the reachability analysis, which further decreases the conservatism of the resulting reachable sets. The resultant reachable sets over-approximate the model-based reachable set for LTI systems (Theorem 3).
An algorithm is proposed for reachability analysis in the existence of process and measurement noise for LTI systems (Algorithm 4).
A comparison between the proposed algorithms and the alternative direction of system identification and model-based reachability analysis is provided.
As discussed in more detail before, there have been a few approaches to infer the reachable sets directly from data, mostly without providing guarantees in the case of noisy data. An alternative approach is to apply well-known system identification approaches and consecutively do model-based reachability analysis. Thus, we include a comparison to a standard system identification method in Section 6 while considering uncertainty bound in the model-based reachability analysis. Interesting recent results on system identification with probabilistic guarantees from finite noisy data include concentration bounds 
. Another very related approach is set membership estimation (see, e.g.,), where the trade-off between conservatism and computational expenses is usually of central importance, which we will also encounter in the course of this paper.
This paper is an extension to our previous work in  where we introduced the basic idea of computing the reachable set by using matrix zonotopes over-approximating the set of models consistent with the data. In this work, we significantly extend and improve the ideas in  by introducing constrained matrix zonotopes and its essential set of operations which allow the incorporation of side information about the unknown model and provide less conservative reachable sets. We also enhance the proposed nonlinear Lipschitz reachability analysis method in . Furthermore, we propose a new approach for computing the reachable sets of polynomial systems, and we consider measurement noise for linear systems. Our initial ideas in  have been utilized in different applications like data-driven predictive control  and set-based estimation .
The rest of the paper is organized as follows: the problem statement and preliminaries are defined in Section 2. The new set representation (constrained matrix zonotope) is proposed in Section 3. Data-driven reachability analysis for LTI systems is proposed in Section 4 including a framework to include prior knowledge into the analysis approach. Then, we extend the proposed approach to nonlinear systems in Section 5. The introduced approaches are applied to multiple numerical examples in Section 6, and Section 7 concludes the paper.
2 Problem Statement and Preliminaries
We start by defining some set representations that are used in the reachability analysis.
2.1 Set Representations
We define the following sets:
A linear map is defined as . Given two zonotopes and , the Minkowski sum is
For simplicity, we use the notation instead of to denote the Minkowski sum as the type can be determined from the context. Similarly, we use to denote . We define the Cartesian product of two zonotopes and by
(Matrix Zonotope [3, p.52]) Given a center matrix and generator matrices , a matrix zonotope is defined as
We use the shorthand notation for matrix zonotopes.
Zonotopes have been extended in  to represent arbitrary convex shapes by applying constraints on the factors multiplied with the generators.
(Constrained Zonotope [42, Prop. 1]) An -dimensional constrained zonotope is defined by
where is the center, is the generator matrix and and denote the constraints. In short, we use the shorthand notation for constrained zonotopes.
2.2 Problem Statement
We consider a discrete-time system
where denotes the noise bounded by a noise zonotope , the input bounded by an input zonotope , the measured state that is additionally corrupted by measurement noise bounded by measurement noise zonotope , and the initial state of the system bounded by the initial set .
Reachability analysis computes the set of states which can be reached given a set of uncertain initial states and a set of uncertain inputs . More formally, it can be defined as follows:
(Exact Reachable Set) The exact reachable set after time steps subject to inputs , noise , is the set of all states trajectories starting from initial set after steps:
|Input zonotope at time .|
|Initial reachable set.|
|Zonotope bounding the process noise.|
|Matrix zonotope bounding the process noise matrix .|
|Constrained matrix zonotope as a subset of given additional constraints from data for LTI system.|
|Matrix zonotope consistent with the data and noise bound for LTI system.|
|Constrained matrix zonotope as a subset of given additional constraints from data for LTI system.|
|Exact reachable set.|
|Reachable set computed using zonotopes for LTI system.|
|Reachable set computed given exact noise description using constrained zonotopes for LTI system.|
|Reachable set computed given side information using constrained zonotopes for LTI system.|
|Reachable set computed using zonotopes in case of measurement noise for LTI system.|
|Reachable set computed using constrained zonotopes in case of measurement noise for LTI system.|
|Reachable set computed using a data-based approximation in case of measurement noise for LTI system.|
|Reachable set computed using zonotopes for polynomial systems.|
|Reachable set computed using constrained zonotopes for polynomial systems.|
|Reachable set computed given side information using constrained zonotopes for polynomial systems.|
|Reachable set for Lipschitz nonlinear systems.|
We aim to compute an over-approximation of the exact reachable sets when the model of the system in (6) is unknown, but input and noisy state trajectories are available. More specifically, we aim to compute data-driven reachable sets in the following cases:
Instead of having access to a mathematical model of the system, we consider input-state trajectories of different lengths , denoted by , , . We collect the set of all data sequences in the following matrices
Let us further denote
The total amount of data points from all available trajectories is denoted by , and we denote the set of all available data by . Note that when dealing with measurement noise, we will consider the trajectories instead of .
We denote the unknown process noise in state trajectory by . It follows directly that the stacked matrix of the noise in the collected data:
is an element of the set where , with
Note that is the matrix zonotope resulting from the concatenation of multiple noise zonotopes as follows:
, . We denote the Kronecker product by . We denote the element at row and column of matrix by and column of by . For vectors, we denote the element of vector by . For a given matrix , denotes a matrix of same size as A with zero entries everywhere except for the value at row i and column j. The vectorization of a matrix is defined by . The element-wise multiplication of two matrices is denoted by . The right inverse is denoted by . We denote the over-approximation of a reachable set by an interval by . We define also for time steps
Finally, we denote all system matrices that are consistent with the data by :
By assumption, . Throughout the paper, we keep record of the used set descriptions in Table 1.
3 Constrained Matrix Zonotope
Inspired by constrained zonotopes, we propose to extend the notion of matrix zonotopes to constrained matrix zonotopes as follows:
(Constrained Matrix Zonotope) Given a center matrix and a number of generator matrices , as well as matrices and constraining the factors , a constrained matrix zonotope is defined by
Furthermore, we define the shorthand notation for the constrained matrix zonotope.
The constrained matrix zonotopes are closed under Minkowski sum and multiplication by a scalar which can be done as follows:
For every , , and the following identities hold
A proof of Proposition 1 is provided in the Appendix. During the propagation of the reachable sets, we additionally need to multiply the constrained matrix zonotope by a zonotope or constrained zonotope. The result of both operations can be over-approximated by a constrained zonotope. We provide these operations in the next two propositions.
For every , and the following identity holds
A proof of Proposition 2 is provided in the Appendix.
For every , and the following identity holds
A proof of Proposition 3 is provided in the Appendix.
4 Data-Driven Reachability for Linear Systems
We consider in this section LTI systems given (i) data corrupted by process noise, (ii) data corrupted by process noise while having additional prior information on the system matrices, and (iii) data corrupted by process noise and measurement noise.
4.1 Linear Systems with Process Noise
Consider a discrete-time linear system
where , and . Due to the presence of noise, there generally exist multiple matrices that are consistent with the data. To provide reachability analysis guarantees, we need to consider all models that are consistent with the data. Therefore, we are interested in computing a set that contains all possible that are consistent with the input-state measurements and the given noise bound. We build upon ideas from  to our zonotopic noise descriptions, which yields a matrix zonotope paving the way to a computationally simple reachability analysis.
Given input-state trajectories of the system (26) and a matrix such that
then the matrix zonotope
contains all matrices that are consistent with the data and the noise bound, i.e., .
For any , we know that there exists a such that
Every can be represented by a specific choice , , , that results in a matrix inside the matrix zonotope :
Multiplying from the right to both sides in (29) yields
The condition (27) in Lemma 1 requires that there exists a right-inverse of the matrix . This is equivalent to requiring this matrix to have full row rank, i.e. . This condition can be easily checked given the data. Note that for noise-free measurements this rank condition can also be enforced by choosing the input persistently exciting of order (compare to [46, Cor. 2]).
To guarantee an over-approximation of the reachable sets for the unknown system, we need to consider the union of reachable sets of all that are consistent with the data. We apply the results of Lemma 1 and do reachability analysis to all systems in the set . Let denotes the reachable set computed based on the noisy data using matrix zonotopes. We propose Algorithm 1 to compute as an over-approximation of the exact reachable set . The set of models that is consistent with data is computed in line 2 which is then utilized in the recursion of computing the reachable set in line 4. The following theorem proves that .
The reachable set computed based on the model can be found using
Since according to Lemma 1 and both and start from the same initial set , it holds that .
Lemma 1 provides a matrix zonotope which comprises all that are consistent with the data and the noise bound. However, not all elements of the matrix zonotope correspond to a system in (26) that can explain the data given the noise bound, i.e., is in fact a superset of all that are consistent with the data (). As discussed in [5, 25], might not be explainable by for all possible . More precisely, there might not exists a solution to the system of linear equations
for all . An exact description for all systems consistent with the data and the noise bound would therefore be the set
where denotes a matrix containing a basis of the kernel of . Representing and is not possible using state of the art set representations. Therefore, we propose the constrained matrix zonotope introduced in Section 3 as a new set representation that can represent the sets and thereby to compute a less conservative reachable set at the cost of increasing the computational complexity in Algorithm 2. Due to adding constraints, is a constrained zonotope different from which is a zonotope. We first compute the exact noise description in line 2 to line 5. Then, we compute the set of models that is consistent with the exact noise description in line 6 which is further utilized in the recursion of computing the reachable set in line 8. The following theorem proves that .
As pointed out in , the condition for the existence of a solution to the system of linear equations
can be reformulated via the Fredholm alternative as
which means that any vector in the kernel of must also lie in the kernel of . Since contains a basis of the kernel of , another equivalent condition for the existence of a solution in (34) is hence
Considering the constraint (35) together with the bounding matrix zonotope , we find:
Rearranging (36) results in
Algorithm 2 provides a less conservative description of the data-driven reachable set compared to Algorithm 1 by utilizing a less conservative description of the set of systems consistent with the data. To be more precise, in (32) with (33) is an equivalent description of all systems consistent with the data and the noise bound (compare [25, Lemma 8]). However, applying the reachability analysis in line 8 of Algorithm 2 requires multiplying constrained matrix zonotopes by zonotopes and constrained zonotopes. For this multiplication, we introduced a guaranteed over-approximation in Proposition 2 and Proposition 3, respectively, which hence introduces conservatism into the proposed reachability analysis approach.
Note that initial zonotope captures all the uncertainty in the initial state. Next, we provide a general framework for incorporating side information about the unknown model.
4.2 Linear Systems with Side information
Consider a scenario in which we have prior side information about the unknown model from the physics of the problem or any other source. It would be beneficial to make use of this side information to have less conservative reachable sets. In the following, we propose a framework to incorporate side information about the unknown model like decoupled dynamics, partial model knowledge, or prior bounds on entries in the system matrices. More specifically, we consider any side information that can be formulated as
where , , and are matrices defining the side information which is known to hold for the true system matrices . Here, the absolute and are element-wise operators. To incorporate such side information into the reachability analysis, we utilize once again the newly introduced concept of constrained matrix zonotopes. We introduce a reachability analysis in Algorithm 3 on the basis of the set of system matrices that are consistent with the data (including the less conservative noise handling in ) as well as the a priori known side information in (37). We denote the reachable set computed based on the side information by . Algorithm 3 summarizes the required computation to incorporate the side information. After setting in line 1, we compute the exact noise description and exact set of models consistent with the noisy data in lines 2:6 similar to Algorithm 2. Next, we compute the set of models consistent with the side of the information in line 7 to line 14. Finally, we compute the recursion of the reachable sets in line 16. The following theorem proves that .