## 1 Introduction

Reachability analysis for continuous and hybrid systems has been an attractive research topic for the last two decades since it is an essential problem for verification of safety-critical CPS. In this context, numerous techniques and tools have been proposed. Reachability analysis using zonotopes [2, 14] and support functions [15, 12] are efficient approaches when dealing with linear, continuous and hybrid systems. For nonlinear, continuous and hybrid systems, dReal[18] using reachability analysis and Flow [7] using Taylor model are well-known and efficient approaches.

State-space explosion is the main challenge that limits the applicability of over-approximation based approaches [2, 12, 7] to small and medium scale systems. Therefore, some promising approaches have been proposed recently to deal with large scale systems. For linear cases, the simulation-equivalent reachability analysis [3, 10] utilizing the *generalized star set* as the state-set representation has shown an impressive result by successfully dealing with linear systems up to state variables. In this approach, the *discrete simulation-equivalent* reachable set of a linear ODE system can be computed efficiently using standard ODE solvers by taking advantage of the superposition property. This simulation-based method is practically useful for falsification of large systems. Another method utilizes order-reduction abstraction [22, 16] in which a large system can be abstracted to a smaller system with bounded error. The smaller system and the error are then used to construct the reachable set or for checking the safety of the original, large system. For nonlinear cases, C2E2 [9, 11] utilizing simulation has shown great improvement on time performance and scalability in comparison with other methods.

Although many methods have been developed for reachability analysis of CPS, most of these methods focus on CPS with ODE dynamics. There is a lack of methodology in analyzing systems with DAE dynamics. To the best of the authors’ knowledge, there exists only one reachability analysis approach for nonlinear, semi-explicit index-1 DAE systems [1] using zonotopes as state-set representation that has been proposed recently. Dealing with index-1 DAE is slightly different from dealing with pure ODE because, with a consistent initial condition, a semi-explicit index-1 DAE can be converted to an ODE. As CPS involving high-index DAE dynamics appear extensively in engineering and science such as multi-body mechanics, electrical circuit design, heat and gas transfer, chemical process, atmospheric physics, thermodynamic systems, and water distribution network [5], there is an urgent need for novel reachability analysis methods and tools that can verify (or falsify) the safety property of such CPS. Solving this challenging problem is the main contribution of this research.

In this paper, we investigate the reachability analysis for large linear DAE systems with index up to , which appear widely in practice. There are a variety of definitions for the index of a linear DAE. However, throughout the paper, we use the concept of *tractability index* proposed in [19] to determine *the index of a linear DAE*. Our approach consists of three main steps: decoupling and consistency checking, reachable set computation, and safety verification or falsification. The novelty of our approach comes from its objective in dealing with high-index DAE which is a popular class of dynamics that has not been addressed in existing literature.

In the first step, we use the Marz decoupling method [19, 4] to decouple a high-index DAE into one ODE subsystem and one or several algebraic constraint (AC) subsystems. The core step in decoupling is constructing a set of admissible projectors which has not previously been discussed deeply in existing literature. In this paper, we propose a novel algorithm that can construct such admissible projectors for a linear DAE system with index up to . Additionally, we define *a consistent space* for the DAE because, unlike ODE reachability analysis where the initial set of states can be freely defined by a user, in order to guarantee a numerical solution for the DAE, the initial state and inputs of the system must be consistent and satisfy certain constraints. It is interesting to emphasized that the decoupling and consistency checking methods used in our approach can be combined with existing over-approximation reachability analysis methods [2, 12] to compute an over-approximated reachable set for high-index, linear DAE systems with small and medium dimensions.

The second step in our approach is reachable set computation. Since our main objective is to verify or falsify large linear DAEs, we extend ODE simulation-based reachability analysis to DAEs. In particular, we modify the generalized star-set proposed in [3] to enhance the efficiency in checking the initial condition consistency and safety for DAEs. From a consistent initial set of states and inputs, the reachable set of a DAE can be constructed by combining the reachable sets of its subsystems. It is also worth pointing out that the piecewise constant inputs assumption for ODE with inputs used in [3] may lead a DAE system to *impulsive behavior*. Therefore, in this paper, we assume the inputs applied to the system are *smooth functions*. These kinds of inputs can be obtained by *smoothing* piecewise constant inputs with filters.

The last step in our approach is to verify or falsify the safety property of the DAE system using the constructed reachable set. In this paper, we consider linear safety specifications. We are interested in checking the safety of the system in a specific direction defined using a directional matrix. Using the modified star-set and the directional matrix, checking the safety property can be solved efficiently as a low-dimensional feasibility linear programming problem. In the case of violation, our approach generates a trace that falsifies the system safety.

## 2 Preliminaries

### 2.1 Linear DAE system

We are interested in the reachability analysis of a high-index large linear DAE system described as follows:

(1) |

where

is the state vector of the system;

, are the system’s matrices in which is*singular*; and is the input of the system. Let be the

-dimensional identity matrix. The regularity, the tractability index, the admissible projectors, the fixed-step bounded-time simulation, and the bounded-time simulation-equivalent reachable set of the system are defined below.

###### Definition 1 (Regularity [8])

The pair is said to be regular if is not identically zero.

###### Remark 1

For *any specified initial conditions*, the regularity of the pair guarantees the existence and uniqueness of a solution of the system (1).

###### Definition 2 (Tractability index [19])

Assume that the DAE system (1) is solvable, i.e., the matrix pair is regular. A matrix chain is defined by:

(2) |

where are projectors onto , i.e., , and . Then, there exists an index such that is nonsingular and all are singular for . It is said that the system (1) has tractability index-. In the rest of the paper, we use the term “index” to state for the “tractability index” of the system.

###### Definition 3 (Admissible projectors [19])

Given a DAE with tractability index-, the projectors in Definition 2 are called admissible if and only if they satisfy the following property: , .

###### Definition 4 (Fixed-step, bounded-time simulation)

Given consistent initial state and input , a time bound , and a time step , the finite sequence:

is a -simulation of the DAE system (1) if and only if for all , is the state of the system trajectory starting from when provided with input function for . If there is no input, .

The consistent condition for the initial state and input will be discussed in detail in Section 4. From the fixed-step, bounded-time simulation of a DAE system, we define the following bounded-time, simulation-equivalent reachable set of the DAE system.

###### Definition 5 (Bounded-time, simulation-equivalent reachable set)

Given sets of consistent initial state and input , the bounded-time, simulation-equivalent reachable set of the system (1) is the set of all states that can be encountered by any -simulation starting from any and input .

Let be the unsafe set of the DAE system (1) in which is the state vector of the system, is the *unsafe matrix* and is the *unsafe vector*. Given sets of consistent initial state and input , the simulation-based safety verification and falsification problem is defined in the following.

###### Definition 6 (Simulation-based safety verification and falsification)

The DAE system (1) is said to be “simulationally safe” up to time if and only if its simulation-equivalent reachable set, , and the unsafe set, , are disjoint, i.e., . Otherwise, it is simulationally unsafe.

The DAE system is said to be “simulationally falsifiable” if and only if it is simulationally unsafe and there exists a simulation, , that leads the initial state, , of the system to an unsafe state, .

The main objective of the paper is to compute the simulation-equivalent reachable set, , of the DAE system and use it to verify or falsify the safety property of the system. In the rest of the paper, we use the term *reachable set* to stand for *simulation-equivalent reachable set*. Next, we define a *modified star set* which is used as the state-set representation of the DAE system.

### 2.2 Modified star set

In our approach, we use a modified star set to represent the reachable set of the DAE system. The modified star set is slightly different from the generalized star set [3] because it does not have a *center vector* and is only defined on a star’s basis matrix.

###### Definition 7 (Modified star set)

A modified star set (or simply star) is a tuple where is a star basis matrix and is a linear predicate. The set of states represented by the star is given by:

(3) |

where , , and is the number of linear constraints.

One can see that if we let be a *center vector* and , then the modified star set becomes the generalized star set proposed in [3]. The benefit of the modified star set come from its form given as a *matrix-vector product* which is convenient for checking initial condition consistency and safety properties. In the rest of the paper, we will refer to both the tuple and the set of states as .

To construct the reachable set of the DAE system (1), we decouple the system into subsystems where is the index of the DAE system. The underlining technique used in our approach is the Marz decoupling method utilizing admissible projectors which is presented in detail in the following section.

## 3 Decoupling

In this section, we discuss how to decouple a high-index DAE system into one ODE subsystem and one or several AC subsystems using the matrix chain and admissible projectors defined in the previous section. Since we are particularly interested in DAE systems with index up to 3, the proofs of decoupling process for index-1, -2, and -3 are given in detail. A generalization of decoupling for a DAE with arbitrary index is presented in [19]. As the construction of admissible projectors used in decoupling has not been discussed clearly in existing literature, in this section, we propose a method and an algorithm to solve this problem.

Proof is given in Appendix 0.A.1.

Proof is given in Appendix 0.A.2.

Proof is given in Appendix 0.A.3.

It should be noted that the AC subsystems and in Lemma 2 and 3 are called algebraic constraints, though they contain the derivatives of and . This is because the explicit forms of these algebraic constraints can be obtained if we further extend the derivatives using the corresponding ODE subsystems. In addition, one can see that for a DAE system with index- or -, a set of admissible projectors need to be constructed for decoupling. In the following, we give a Proposition and Lemmas that are used to construct such admissible projectors.

###### Proposition 1 (Orthogonal projector on a subspace)

Given a real matrix such that

, the Singular-Value Decomposition (SVD) of

has the form:(4) |

where and . Then, the matrix is an orthogonal projector on , i.e., , and .

Proof is given in Appendix 0.A.4.

For an index- or - DAE system, using Proposition 1, we can construct a set of projectors of the matrix chain defined in Equation (2). However, these projectors are not yet admissible, because . Instead, the admissible projectors can be constructed based on these inadmissible projectors using the following Lemmas.

###### Lemma 4 (Admissible projectors for an index-2 DAE system)

Proof is given in Appendix 0.A.5.

###### Lemma 5 (Admissible projectors for an index-3 DAE system)

Given an index-3 DAE system described by (1), let , and respectively be the orthogonal projectors of , and of the matrix chain defined in Equation (2). We define the following projectors and the corresponding new matrices for the matrix chain as:

where and . Let be the orthogonal projector on and , then the following projectors and are admissible: .

Proof is given in Appendix 0.A.6.

Lemmas 4 and 5 are the constructions of admissible projectors for index- and - DAE systems. Algorithm 3.1 clarifies how to construct admissible projectors for a DAE system with an index up to . Next, based on the decoupled DAE system, we discuss the consistent condition of the system and analyze the system behavior under the effect of input functions.

## 4 Consistency

In this section, we discuss the consistent condition for a DAE system. Using the decoupled DAE system, the consistent condition for the initial state and inputs is derived. Additionally, the piecewise constant assumption on the inputs used in [3] for ODE systems may lead to *impulsive behavior* in high-index DAE systems. To avoid this, we limit our problem to *smooth* and *specific-user-defined* inputs. As a result, DAE systems with inputs can be converted to autonomous DAE systems, where *consistent spaces* for the initial states and inputs can be conveniently defined and checked. Furthermore, the reachable set computation is executed efficiently using a decoupled autonomous DAE system.

Using Lemmas 1, 2, and 3, to guarantee a solution for the DAE system, the initial states and inputs must satisfy the following conditions:

(5) |

Assuming that the consistent condition is satisfied, Lemmas 2 and 3 indicate the solution of the system involves the derivatives of the input functions and . In cases where we apply piecewise constant inputs to a high-index DAE system, the impulsive behavior may appear in the system at an exact discrete time point . For example, let be a step function in , then , where is the Dirac function describing an impulse. To avoid such impulsive behavior and do reachability analysis for high-index DAE systems, we limit our approach to smooth inputs which are governed by the following ODE: , where is the user-defined input matrix, and is the set of initial inputs.

###### Remark 2

By introducing the input matrix , we limit the safety verification and falsification of a high-index DAE system to a class of *specific-user-defined* inputs. If , then the input set is a set of constant inputs.

Given a user-defined input matrix , a DAE system described by (1) can be converted to an equivalent autonomous DAE system of the following form:

(6) |

where , and the state of the original DAE is: .

Similar to the original DAE system, the autonomous DAE system (6) can be decoupled to form one autonomous ODE subsystem and one or several AC subsystems. It should be noted that the autonomous DAE system has the same index as the original one.

We have discussed the conversion of a DAE system with user-defined input to an autonomous DAE system. Next, we derive the *consistent space* for the initial condition of an autonomous DAE system. All previous results apply to these systems given that .

###### Definition 8 (Consistent Space for an autonomous DAE system)

Consider an autonomous DAE system () defined in Equation (1) by letting . From this, we define in the following a “consistent matrix” as:

then, is the consistent space of the system , where denotes the null space of the matrix .

An initial state is consistent if it is in the consistent space, i.e., . The consistent matrix and consistent space is introduced because it is useful and convenient for checking the consistency of an initial set of states represented using a star set. For example, assume that the initial set of states is defined by , then this set is consistent for all satisfying the predicate if . This means that we require consistency for all points in the initial set. With a consistent initial set of states, we investigate the reachable set computation and safety verification/falsification of an autonomous DAE system in the next section.

## 5 Reachability analysis

### 5.1 Reachable set computation

The reachable set of an autonomous DAE system is constructed by combining the reachable set of all of its decoupled subsystems. The reachable set of all AC subsystems can be derived from the reachable set of the ODE subsystem, which can be computed efficiently using existing ODE solvers. We first discuss the reachable set computation of the ODE subsystem by exploiting its *superposition property*. Then, the reachable set of the autonomous DAE system is constructed conveniently using only matrix addition and multiplication.

Let be the initial set of states of an autonomous DAE system defined in (1) by letting . Assume that the initial set of states, , satisfies the consistent condition. After decoupling, the initial set of states of the ODE subsystem is obtained as follows: where , is the index of the DAE system, and are defined in Lemma 1 or 2 or 3 corresponding to the index .

Then, for any , we have . The solution of the ODE subsystem at time is given by: , where and . Therefore, the reachable set of the ODE subsystem at anytime is also a star set defined by . Using existing ode solvers, we can construct the matrix at anytime . From , the reachable set of the autonomous DAE system can be obtained using the following Lemma.

###### Lemma 6 (Reachable set construction)

Given an autonomous DAE system defined in Equation (1) where and a consistent initial set of states , let be the reachable set at time of the corresponding ODE subsystem after decoupling. Then, the reachable set at time of the system is given by , where is a “reachable set projector” defined below.

(7) |

Proof is given in Appendix 0.B.1.

The Algorithm 5.1 describes how to construct a reachable set of an autonomous DAE system. Next, from the constructed reachable set, we discuss how to verify or falsify the system safety property.

### 5.2 Safety verification and falsification

By utilizing the star set to represent the reachable set of a DAE system, the safety verification and falsification problem is solved in the following manner. Let be the unsafe set of an autonomous DAE system and assume that we want to check the safety of the system at the time step . This is equivalent to checking subject to , where is the basic matrix of the reachable set of the system at time computed using Algorithm 5.1. Combining these constraints, the problem changes to checking the feasibility of the following linear predicate: , where and . This can be solved efficiently using existing linear programming toolboxes. Algorithm 5.2 illustrates how to verify or falsify the safety property of an autonomous DAE system. In the next section, we evaluate our approach using a set of DAE benchmarks with several to thousands of states.

## 6 Experimental Results

The detailed description of how our approach works using the index-2, interconnected rotating masses system [21] is presented in Appendix 0.C. In this section, we first analyze the time performance of our approach using the index-2, two-dimensional semidiscretized Stokes Equation benchmark [20]. Then, we test and compare the scalability of our approach with the well-known tool SpaceEx [12] using a set of benchmarks. Using SpaceEx for DAE systems verification is non-trivial. To do this, we first decouple a DAE and check the consistency condition of the initial set of states and inputs. Then, we construct an automaton with the decoupled ODE subsystem. Finally, the state of the DAE system is declared as a set of invariants of the automaton using the “reachable set projector” in Lemma 6. Our approach is implemented in a tool called *Daev*^{1}^{1}1https://www.dropbox.com/s/w1mrbl3nw31n2gs/daev.zip?dl=0 using python and its standard packages numpy, scipy, and mathplotlib. All experiments were done on a computer with the following configuration: Intel Core i7-6700 CPU @ 3.4GHz × 8 Processor, 62.8 GiB Memory, 64-bit Ubuntu 16.04.3 LTS OS.

###### Example 1 (Semidiscretized Stokes Equation [20])

This example studies the safety of a Stokes equation that describes the flow of an incompressible fluid in a two-dimensional spatial domain . The mathematical description of the Stokes-equation is given in Appendix 0.D. An index-2 DAE system is derived from the Stokes-equation by discretizing the domain by a number of uniform square cells. Let be the number of discretized segments of the domain on the x- or y-axes, then the dimension of the DAE system is . Additionally, we are interested in the velocity along the x- and y- axes, and , of the fluid in the *central cell* of the domain . The unsafe set of the system is defined: .

By increasing the number of cells used to discretize the domain , we can produce an index-2 DAE system with arbitrarily large dimension. We evaluate the time performance of our approach via three scenarios. First, we discuss how the times for decoupling, reachable set computation, and safety checking are affected by changes in the system dimension. Second, we analyze the reachable set computation time along with the *width* of the basic matrix of the initial set , i.e., the number of the initial basic vectors. Finally, because the reachable set of the system is constructed from the reachable set of its corresponding ODE subsystem, which is computed using ODE solvers as shown in Algorithm 5.1, we investigate the time performance of reachable set computation using different ODE solving schemes.

Table 1 presents the verification time, *V-T*, for the Stokes-equation benchmark with different dimensions. The verification time is broken into three components measured in seconds: decoupling time *D-T*, reachable set computation time *RSC-T*, and checking safety time *CS-T*. Table 1 shows the decoupling and reachable set computation times dominate the time for verification process. In addition, these times increase as the system size grows. The time for checking safety is almost unchanged and very small. This happens because the size of the feasibility problem defined in Algorithm 5.2 is unchanged and usually small when we only check the safety in some specific directions defined by the unsafe matrix in the Algorithm.

n | ||||||

D-T | ||||||

RSC-T | ||||||

CS-T | ||||||

V-T |

Since the reachable set the Stokes-equation benchmark is constructed by simulating its corresponding ODE subsystem with each initial vector of its initial basic matrix, the time for computing the reachable set of the Stokes-equation depends linearly on the number of the initial basic vectors . Table 2 shows the reachable set computation time, , for the Stokes-equation of dimension versus the number of the initial basic vectors .

k |

Comments

There are no comments yet.