# Data-Driven Reachability Analysis from Noisy Data

We consider the problem of computing reachable sets directly from noisy data without a given system model. Several reachability algorithms are presented, and their accuracy is shown to depend on the underlying system generating the data. First, an algorithm for computing over-approximated reachable sets based on matrix zonotopes is proposed for linear systems. Constrained matrix zonotopes are introduced to provide less conservative reachable sets at the cost of increased computational expenses and utilized to incorporate prior knowledge about the unknown system model. Then we extend the approach to polynomial systems and under the assumption of Lipschitz continuity to nonlinear systems. Theoretical guarantees are given for these algorithms in that they give a proper over-approximative reachable set containing the true reachable set. Multiple numerical examples show the applicability of the introduced algorithms, and accuracy comparisons are made between algorithms.

There are no comments yet.

## Authors

• 11 publications
• 4 publications
• 9 publications
• 11 publications
• ### Data-Driven Reachability Analysis Using Matrix Zonotopes

In this paper, we propose a data-driven reachability analysis approach f...
11/17/2020 ∙ by Amr Alanwar, et al. ∙ 0

• ### Enhancing Data-Driven Reachability Analysis using Temporal Logic Side Information

This paper presents algorithms for performing data-driven reachability a...
09/15/2021 ∙ by Amr Alanwar, et al. ∙ 0

• ### Sampling-based Reachability Analysis: A Random Set Theory Approach with Adversarial Sampling

Reachability analysis is at the core of many applications, from neural n...
08/24/2020 ∙ by Thomas Lew, et al. ∙ 0

• ### On the Decidability of Reachability in Linear Time-Invariant Systems

We consider the decidability of state-to-state reachability in linear ti...
02/19/2018 ∙ by Nathanaël Fijalkow, et al. ∙ 0

• ### Adaptive Parameter Tuning for Reachability Analysis of Linear Systems

Despite the possibility to quickly compute reachable sets of large-scale...
06/22/2020 ∙ by Mark Wetzlinger, et al. ∙ 0

• ### A Physics-Constrained Data-Driven Approach Based on Locally Convex Reconstruction for Noisy Database

Physics-constrained data-driven computing is a hybrid approach that inte...
07/26/2019 ∙ by Qizhi He, et al. ∙ 0

• ### Functional sets with typed symbols: Framework and mixed Polynotopes for hybrid nonlinear reachability and filtering

Verification and synthesis of Cyber-Physical Systems (CPS) are challengi...
09/15/2020 ∙ by Christophe Combastel, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

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 [30]. 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 [9]. The resulting enclosure is then inflated by a bound on the truncation error. Other set representations are polyhedra [41], zonotopes [22], (sparse) polynomial zonotopes [26], ellipsoids [29], and support functions [32]. 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 [16] 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 [33]. 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 [13]

, 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

[14] 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 [15] 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 [19], 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 available.

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:

• An algorithm is proposed to compute the reachable sets of an unknown control system from noise-corrupted input-state measurements using matrix zonotopes (Algorithm 1). The resultant reachable sets are shown to over-approximate the model-based reachable sets for LTI systems (Theorem 1).

• A new set representation named constrained matrix zonotope (Definition 5) and its essential operations are proposed.

• The constrained matrix zonotope is exploited in Algorithm 2 by providing a less conservative reachable sets using the exact noise description. The computed reachable sets over-approximate the model-based reachable sets for LTI systems (Theorem 2).

• 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).

• An algorithm (Algorithm 5) is proposed for computing the reachable sets that are guaranteed to over-approximate the exact reachable set for polynomial systems (Theorem 4). A variant of the LTI side information framework can be used for the polynomial systems.

• An algorithm is proposed for computing the reachable set (Algorithm 6), which results in reachable sets that are guaranteed to over-approximate the exact reachable sets of nonlinear systems under a Lipschitz continuity assumption (Theorem 5).

• 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 [37]

. Another very related approach is set membership estimation (see, e.g.,

[38]), 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 [1] 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 [1] 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 [1]. 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 [1] have been utilized in different applications like data-driven predictive control [2] and set-based estimation [7].

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:

###### Definition 1

(Zonotope [28]) Given a center and

generator vectors in a generator matrix

, a zonotope is defined as

 Z={x∈Rn∣∣x=cZ+γZ∑i=1β(i)g(i)Z,−1≤β(i)≤1}. (1)

We use the shorthand notation for zonotopes.

A linear map is defined as . Given two zonotopes and , the Minkowski sum is

 Z1+Z2=⟨cZ1+cZ2,[GZ1,GZ2]⟩. (2)

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

 Z1×Z2 =⟨[cZ1cZ2],[GZ100GZ2]⟩. (3)
###### Definition 2

(Matrix Zonotope [3, p.52]) Given a center matrix and generator matrices , a matrix zonotope is defined as

 M={X∈Rn×p∣∣X=CM+γM∑i=1β(i)G(i)M,−1≤β(i)≤1}. (4)

We use the shorthand notation for matrix zonotopes.

Zonotopes have been extended in [42] to represent arbitrary convex shapes by applying constraints on the factors multiplied with the generators.

###### Definition 3

(Constrained Zonotope [42, Prop. 1]) An -dimensional constrained zonotope is defined by

 (5)

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

 x(k+1)=f(x(k),u(k))+w(k),y(k)=x(k)+v(k). (6)

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:

###### Definition 4

(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:

 RN={x(N)∈Rn∣∣ x(k+1)=f(x(k),u(k))+w(k), x(0)∈X0,u(k)∈Uk,w(k)∈Zw: ∀k∈{0,...,N−1}}. (7)

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:

1. LTI systems in Subsection 4.1:

 x(k+1)=Atrx(k)+Btru(k)+w(k).
2. LTI systems given additional side information about the unknown model in Subsection 4.2.

3. LTI systems with measurement noise in Subsection 4.3:

 x(k+1)=Atrx(k)+Btru(k)+w(k),y(k)=x(k)+v(k).
4. Polynomial systems in Subsection 5.1.

5. Lipschitz nonlinear systems in Subsection 5.2.

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

 X =[x(1)(0)…x(1)(T1)…x(K)(0)…x(K)(TK)], U− =[u(1)(0)…u(1)(T1−1)…u(K)(0)…u(K)(TK−1)].

Let us further denote

 X+ =[x(1)(1)…x(1)(T1)…x(K)(1)…x(K)(TK)], X− =[x(1)(0)…x(1)(T1−1)…x(K)(0)…x(K)(TK−1)].

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 .

### 2.3 Preliminaries

We denote the unknown process noise in state trajectory by . It follows directly that the stacked matrix of the noise in the collected data:

 ^W−=[^w(1)(0)…^w(1)(T1−1)…^w(K)(0)…^w(K)(TK−1)]

is an element of the set where , with

 ~GMw=[G(1)Mw…G(γMw)Mw]. (8)

Note that is the matrix zonotope resulting from the concatenation of multiple noise zonotopes as follows:

 CMw (9) (10)

, . 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

 F=∪Nk=0(Rk×Uk). (11)

Finally, we denote all system matrices that are consistent with the data by :

 NΣ={[AB]|X+=AX−+BU−+W−,W−∈Mw}.

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:

###### Definition 5

(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

 N={X∈Rn×p∣∣ X=CN+γN∑i=1β(i)G(i)N, γN∑i=1β(i)A(i)N=BN,−1≤β(i)≤1}.

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:

###### Proposition 1

For every , , and the following identities hold

 RN1 =⟨RCN1,R~GN1,~AN1,BN1⟩, (12) N1+N2 =⟨CN1+CN2,[~GN1,~GN2],~AN12,[BN100BN2]⟩, (13)

where

 ~AN12=⎡⎣[A(1)N1000]…⎡⎣A(γN1)N1000⎤⎦[000A(1)N2]…⎡⎣000A(γN2)N2⎤⎦⎤⎦.

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.

###### Proposition 2

For every , and the following identity holds

 NZ⊂ [ANZ00],vec(BN)⟩, (14)

where

 Gf=[g(1)f…g(γZγN)f], g(k)f=f(i)G(i)Ng(j)Z,∃ k ∀i={1,…,γN},and j={1,…,γZ} such that k={1,…,γZγN}, (15) f(i)=max(|β(i)L,N|,|β(i)U,N|), (16) β(i)L,N=minβ(i)β(i): γN∑j=1β(j)A(j)N=BN, ∥β∥∞≤1, (17) β(i)U,N=maxβ(i)β(i): γN∑j=1β(j)A(j)N=BN, ∥β∥∞≤1. (18)

A proof of Proposition 2 is provided in the Appendix.

###### Proposition 3

For every , and the following identity holds

 NC⊂ [ANC000AC0],[vec(BN)bC]⟩, (19)

where

 ANC =[vec(A(1)N)…vec% (A(γN)N)], Gf =[g(1)f…g(γCγN)f], g(k)f =f(k)G(i)Ng(j)C,∃ k ∀i={1,…,γN},and j={1,…,γC} such that k={1,…,γCγN}, (20) f(k) =max(|β(i)L,Nβ(j)L,C|,|β(i)L,Nβ(j)U,C|,|β(i)U,Nβ(j)L,C|,|β(i)U,Nβ(j)U,C|), (21) β(i)L,N =minβ(i)β(i): γN∑j=1β(j)A(j)N=BN, ∥β∥∞≤1, (22) β(i)U,N =maxβ(i)β(i): γN∑j=1β(j)A(j)N=BN, ∥β∥∞≤1, (23) β(j)L,C (24) β(j)U,C (25)

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

 x(k+1)=Atrx(k)+Btru(k)+w(k), (26)

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 [25] to our zonotopic noise descriptions, which yields a matrix zonotope paving the way to a computationally simple reachability analysis.

###### Lemma 1

Given input-state trajectories of the system (26) and a matrix such that

 [X−U−]H=I, (27)

then the matrix zonotope

 MΣ=(X+−Mw)H (28)

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

 AX−+BU−=X+−W−. (29)

Every can be represented by a specific choice , , , that results in a matrix inside the matrix zonotope :

 W− =CMw+γMw∑i=1^β(i)MwG(i)Mw.

Multiplying from the right to both sides in (29) yields

 [AB]=(X+−CMw+γMw∑i=1^β(i)MwG(i)Mw)H. (30)

Hence, for all , there exists , , , such that (30) holds. Therefore, for all , it also holds that as defined in (28), which concludes the proof.

###### Remark 1

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 .

###### Theorem 1

Given input-state trajectories of the system in (26) and a matrix as defined in (27), then the reachable set computed in Algorithm 1 over-approximates the exact reachable set, i.e., .

The reachable set computed based on the model can be found using

 Rk+1 (31)

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

 [AB][X−U−]=X+−W−

for all . An exact description for all systems consistent with the data and the noise bound would therefore be the set

 NΣ=(X+−Nw)H (32)

with

 (33)

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 .

###### Theorem 2

Given input-state trajectories of the system in (26) and a matrix as defined in (27), then the reachable set computed in Algorithm 2 over-approximates the exact reachable set, i.e., .

As pointed out in [5], the condition for the existence of a solution to the system of linear equations

 F[A B][X−U−]=X+−W−,

or equivalently

 [X−U−]⊤F⊤[A B]=(X+−W−)⊤ (34)

can be reformulated via the Fredholm alternative as

 [X−U−]~z=0⇒(X+−W−)~z=0,

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

 (X+−W−)[X−U−]⊥=0. (35)

Considering the constraint (35) together with the bounding matrix zonotope , we find:

 (X+−CMw−γZwT∑i=1β(i)G(i)Mw)[X−U−]⊥=0. (36)

Rearranging (36) results in

###### Remark 2

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

 |¯Q[AtrBtr]−¯Y|≤¯R, (37)

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 .