Arvind Adimoolam

Thao Dang
email: thao.dang@univ-grenoble-alpes.fr

Safety Verification of Networked Control Systems by Complex Zonotopes

Keywords: 2012 ACM Subject Classification Computer systems organization, Dependable and fault-tolerant systems and networks, Reliability Safety Verification, Networked Control System, Reachability Analysis, Invariant, Complex Zonotope Digital Object Identifier

Networked control systems (NCS) are widely used in real world applications because of their advantages, such as remote operability and reduced installation costs. However, they are prone to various inaccuracies in execution like delays, packet dropouts, inaccurate sensing and quantization errors. To ensure safety of NCS, their models have to be verified under the consideration of aforementioned uncertainties. In this paper, we tackle the problem of verifying safety of models of NCS under uncertain sampling time, inaccurate output measurement or estimation, and unknown disturbance input. Unbounded-time safety verification requires approximation of reachable sets by invariants, whose computation involves set operations. For uncertain linear dynamics, two important set operations for invariant computation are linear transformation and Minkowski sum operations. Zono-topes have the advantage that linear transformation and Minkowski sum operations can be efficiently approximated. However, they can not encode directions of convergence of trajectories along complex eigenvectors, which is closely related to encoding invariants. Therefore, we extend zonotopes to the complex valued domain by a representation called complex zonotope, which can capture contraction along complex eigenvectors for determining invariants. We prove a related mathematical result that in case of accurate feedback sampling, a complex zonotope will represent an invariant for a stable NCS. In addition, we propose an algorithm to verify the general case based on complex zonotopes, when there is uncertainty in sampling time and in input. We demonstrate the efficiency of our algorithm on benchmark examples and compare it with a stateof-the-art verification tool.

Introduction

The computerized control of spatially distributed components through communication channels is called networked control. The emergence of fast and reliable communication networks has enabled efficient networked control of many applications in aerospace, automation, manufacturing and robotics [START_REF] Ding | Distributed cooperative optimal control of dc microgrids with communication delays[END_REF][START_REF] Kim | Cyber-physical systems: A perspective at the centennial[END_REF][START_REF] Sandberg | Cyberphysical security in networked control systems: An introduction to the issue[END_REF][START_REF] Wang | Network-based modelling and dynamic output feedback control for unmanned marine vehicles in network environments[END_REF][START_REF] Zhang | Networkbased h∞ filtering using a logic jumping-like trigger[END_REF]. Advantages of networked control systems (NCS) include reduced installation cost due to absence of wiring and remote operability. However, they are prone to inaccuracies in execution such as delays, packet dropouts and errors in sensing and quantization.

To mitigate the risk of system failure, we need to verify the safety requirements of NCS in the presence of such inaccuracies. In this paper, we tackle the problem of verifying unbounded time safety of linear networked control systems with a uncertain sampling period for feedback input, inaccurate sensors for output estimation and disturbance input.

Safety verification involves proving that the set of reachable states of the system are contained within a specified safe set. However, exactly computing the set of reachable states of models containing linear differential equations with discontinuous switching between states or vector fields, called affine hybrid systems, is computationally intractable [START_REF] Alur | Formal verification of hybrid systems[END_REF]. Instead, a common approach is to overapproximate this set of reachable states to prove safety. There has been a lot of research on overapproximating the set of reachable states of affine hybrid systems in both bounded and unbounded time [START_REF] Frehse | Spaceex: Scalable verification of hybrid systems[END_REF][START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF][START_REF] Sankaranarayanan | Symbolic model checking of hybrid systems using template polyhedra[END_REF]. But networked control systems have time triggered switching which can be handled accurately by these methods. The reason is explained below.

Challenge of safety verification of NCS.

In an NCS, there is time triggered switching of states with uncertainty in the switching time (sampling time). The approach followed by most techniques [START_REF] Frehse | Spaceex: Scalable verification of hybrid systems[END_REF][START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF][START_REF] Sankaranarayanan | Symbolic model checking of hybrid systems using template polyhedra[END_REF] to overapproximate reachable sets of hybrid systems involves computing intersection between reachable sets and guard conditions by some set representation like polytopes [START_REF] Sankaranarayanan | Symbolic model checking of hybrid systems using template polyhedra[END_REF], zonotopes [START_REF] Girard | Zonotope/hyperplane intersection for hybrid systems reachability analysis[END_REF], ellipsoids [START_REF] Xavier Allamigeon | A fast method to compute disjunctive quadratic invariants of numerical programs[END_REF],

or polynomial sub-level sets like barrier certificates [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF]. In the case of NCS, the guard is on the time of switching, i.e., the clock variable. For very simple hybrid systems having constant vector fields, there is linear relationship between clock variables and other state variables. Then, we can expect to reliably overapproximate the intersection of guard on the clock variable and reachable states containing value of clock variable. But an NCS is more complex involving linear differential equations with time triggered switching. There is an exponential relationship between the clock variable and the reachable states in NCS. Such exponential relationship can not be captured using the well known set representations mentioned above. Instead, a more effective way is proposed in [START_REF] Al Khatib | Stability verification of nearly periodic im-pulsive linear systems using reachability analysis[END_REF][START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched lyapunov function approach[END_REF][START_REF] Fiacchini | Set theory conditions for stability of linear impulsive systems[END_REF][START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF][START_REF] Pettersson | Hybrid system stability and robustness verification using linear matrix inequalities[END_REF] to handle time triggered switching by dividing the switching time into very small sub-intervals and using matrix exponentials to map the reachable states in different sub-intervals. The above methods for stability verification of NCS do not consider additive input in the dynamics. However, our dynamics is more general where we consider consider additive disturbance input. The above methods [START_REF] Al Khatib | Stability verification of nearly periodic im-pulsive linear systems using reachability analysis[END_REF][START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched lyapunov function approach[END_REF][START_REF] Fiacchini | Set theory conditions for stability of linear impulsive systems[END_REF][START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF][START_REF] Pettersson | Hybrid system stability and robustness verification using linear matrix inequalities[END_REF] use ellipsoids and H-polytopes which do not efficiently handle additive input in high dimensions. Ellipsoids provide poor accuracy for approximating Minkowski sum with input sets, which is required in reachability analysis.

Similarly, the complexity of H-polytopes exactly representing Minkowski sums blows up at least exponentially in the dimension of state space [START_REF] Kvasnica | Minkowski addition of convex polytopes[END_REF]. Therefore, we need better set representations to accurately overapproximate the reachable sets at various switching times.

Solution: Generalizing zonotope to complex zonotope

In this context, an effective set representation is a zonotope [START_REF] Girard | Reachability of uncertain linear systems using zonotopes[END_REF], described as a linear combination of real vectors called generators, whose combining coefficients are bounded in absolute value.

Its advantage is that linear transformation and the Minkowski sum can be efficiently computed.

However, we are required to overapproximate the unbounded time reachable set of NCS that typically involves computing invariants. Computing invariants by a set representation requires encoding the directions of convergence of reachable sets towards equilibrium. In NCS, some of the directions for convergence can be along the eigenvectors of the dynamics, which can be complex valued vectors (Theorem 7). But (simple) zonotopes, which are confined to the real valued domain, can fail to capture such complex valued directions of convergence of trajectories. Therefore, our goal is to extend them to the complex number domain to obtain a new set representation called complex zonotope. Complex zonotopes retain the merit of usual zonotopes that linear transformation and the Minkowski sum operations can be computed efficiently. Additionally, they can capture the contraction of reachable sets based on the complex eigenstructure, which is not possible using a real valued zonotope. We provide mathematical evidence (Theorem 7) to support the latter claim. Furthermore, their real projections are geometrically more expressive than usual zonotopes and can represent some non-polytopic sets.

Using complex zonotopes, we propose an algorithm to verify linear safety constraints of a linear NCS with uncertain sampling period, inaccurate output estimation and disturbance input.

The algorithm is a semi-decision procedure which can either terminate when the system is verified to be safe or fail to terminate. In the latter case, the user can set a threshold number of iterations for termination to get a bounded time result. Using the algorithm, we successfully applied our algorithm to verify benchmark examples of NCS [START_REF] Heinz | Benchmark: Reachability on a model with holes[END_REF][START_REF] Ben | Networked cooperative platoon of vehicles for testing methods and verification tools[END_REF] with high dimensions (≥ 12 dimensions (including state and controller input variables), while the simple zonotope based version of our algorithm and a state-of-the-art verification tool [START_REF] Frehse | Spaceex: Scalable verification of hybrid systems[END_REF] both failed to verify them.

In summary, we make the following contributions in this paper.

1.

We introduce the complex zonotope set representation, a geometrically more expressive representation than zonotope, to handle invariant computation in the presence of additive input and time triggered switching.

2.

We propose a theoretical result (Theorem 7) about the existence of complex zonotopic invariants based on eigenvectors of the NCS dynamics.

3.

We extend the previous algorithms for stability verification of NCS [START_REF] Al Khatib | Stability verification of nearly periodic im-pulsive linear systems using reachability analysis[END_REF][START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched lyapunov function approach[END_REF][START_REF] Fiacchini | Set theory conditions for stability of linear impulsive systems[END_REF][START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF][START_REF] Pettersson | Hybrid system stability and robustness verification using linear matrix inequalities[END_REF] without additive input to safety verification in the presence of additive input. We use complex zonotopes containing eigenvectors for computing invariants of NCS more effectively when there is additive input.

4.

As a proof-of-concept, we compare the performance of our complex zonotope with a real zonotope containing concatenation of real and imaginary parts of the complex template. We compare the performance on high dimensional (>9 state space) benchmarks examples in literature. We show that while our complex zonotope is successful in verification, the real zonotope either fails to compute an invariant or computes one with very large bounds above the safety threshold. We also compare with another state-of-the-art tool [START_REF] Frehse | Spaceex: Scalable verification of hybrid systems[END_REF], which also failed to verify the benchmarks.

This paper is an extended version of part of our work presented in the conferences [START_REF] Adimoolam | Augmented complex zonotopes for computing invariants of affine hybrid systems[END_REF][START_REF] Adimoolam | Template complex zonotopes for stability and invariant verification[END_REF][START_REF] Arvind | Template complex zonotopes: a new set representation for verification of hybrid systems[END_REF][START_REF] Arvind | Using complex zonotopes for stability verification[END_REF] and the PhD thesis [START_REF] Santosh | A Calculus of Complex Zonotopes for Invariance and Stability Verification of Hybrid Systems[END_REF]. The extensions and modifications made in this paper are explained in Section 2.

Organization. In Section 2, we review previous research related to our work and draw some comparison. In Section 3, we formalize NCS with uncertainty in sampling time, inaccurate output estimation and unknown open loop input. We explain the relation between safety verification and invariant computation at sampling times for an NCS. In Section 4, we introduce the complex zonotope representation as a generalization of usual zonotopes to complex valued domain. We describe a result that shows how a complex zonotope can specify invariants based on eigenstructure.

We discuss operations on complex zonotopes that are later used to verify safety properties. In Section 5, we describe the procedure for verification using complex zonotopes. The experiments on some benchmark examples and results are discussed in Section 6. We begin by describing in the following the mathematical notation used in this paper.

Notation

The set of real numbers is represented by R, integers by Z, and their positive subsets by R ≥0 and Z ≥0 , respectively. The set of complex numbers is C. Given a subset S of real or complex numbers we denote the set of n-dimensional vectors from S as S n and n × m matrices from S as S n×m .

The i th component of a vector x is x i while the element of the i th row and the j th column of a matrix X is X ij . The numbers of rows of matrix X and size of vector x are rows (X) and rows (x), respectively. The number of columns of X is cols (X). Given any two real vectors x, y such that rows (x) = rows (y), we say x ≤ y, if ∀i ∈ {1, . . . , rows (x)}, x i ≤ y i . The diagonal matrix containing a vector x along its diagonal is denoted by D (x). The identity matrix of size n × n is denoted by I n and a vector of ones of size n is denoted by [START_REF] Adimoolam | Augmented complex zonotopes for computing invariants of affine hybrid systems[END_REF] n×1 . The real part of a complex number x is Re (x) and the imaginary part is Im (x). The absolute value of a complex number x is |x| = Re (x) 2 + Im (x) 2 . The infinity norm of a complex vector x is ∥x∥ ∞ = max rows(x) i=1

|x i |.

For a complex matrix X, |X| is the matrix containing the absolute values of the elements of X, i.e., |X| ij = |X i j|. For any a, b ∈ R {∞, -∞} we use the following notation for intervals.

(a, b) = {x ∈ R | a < x < b} , (a, b] = {x ∈ R | a < x ≤ b} [a, b) = {x ∈ R | a ≤ x < b} , [a, b] = {x ∈ R | a ≤ x ≤ b}
Given two same dimensional subsets of complex vector space U ∈ C n and V ∈ C n , we define their Minkowski sum as follows.

U ⊕ V = {u + v | u ∈ U, v ∈ V }
Given a complex matrix M ∈ C n×n and a non-negative real matrix Υ ∈ R n×n ≥0 , we define a neighborhood of M whose difference with M is less than Υ, as follows.

M ▽ Υ = M + M | M ∈ C n×n , M ≤ Υ .

Related work

This paper is a journal extension of part our work part on safety and stability verification using complex zonotope presented in the conferences [START_REF] Adimoolam | Augmented complex zonotopes for computing invariants of affine hybrid systems[END_REF][START_REF] Adimoolam | Template complex zonotopes for stability and invariant verification[END_REF][START_REF] Arvind | Template complex zonotopes: a new set representation for verification of hybrid systems[END_REF][START_REF] Arvind | Using complex zonotopes for stability verification[END_REF] and a PhD thesis [START_REF] Santosh | A Calculus of Complex Zonotopes for Invariance and Stability Verification of Hybrid Systems[END_REF]. The present work adapts the algorithm used for stability verification [START_REF] Adimoolam | Template complex zonotopes for stability and invariant verification[END_REF][START_REF] Arvind | Using complex zonotopes for stability verification[END_REF] and tackles the problem of safety verification of NCS. In this context, we provide a sufficient condition for checking inclusion of a set of complex zonotopes whose templates are in the neighborhood of a complex zonotope, inside another complex zonotope. This generalizes inclusion checking between two complex zonotopes from the conference papers [START_REF] Adimoolam | Augmented complex zonotopes for computing invariants of affine hybrid systems[END_REF][START_REF] Santosh | A Calculus of Complex Zonotopes for Invariance and Stability Verification of Hybrid Systems[END_REF] to inclusion checking between a set of complex zonotopes and another complex zonotope. Our conference paper on safety verification [START_REF] Adimoolam | Augmented complex zonotopes for computing invariants of affine hybrid systems[END_REF] proposes an approximation of intersection of complex zonotope with linear constraints which can be coarse because complex zonotopes are not closed under intersection with linear constraints. If using this method for an NCS with uncertain sampling time, intersection of complex zonotopes with sampling time constraints can be inaccurately approximated. Instead, in this paper we make use of matrix exponential maps to compute an invariant.

A number of set representations have also been developed for verification of hybrid systems and also numerical programs. Linear relationships between state variables can be represented using polytopes [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF] and their variants, such as template polyhedra [START_REF] Sankaranarayanan | Symbolic model checking of hybrid systems using template polyhedra[END_REF], hypercubes [START_REF] Tiwari | Generating box invariants[END_REF], octagons [START_REF] Miné | The octagon abstract domain[END_REF],

zonotopes [START_REF] Girard | Reachability of uncertain linear systems using zonotopes[END_REF] and tropical polyhedra [START_REF] Xavier Allamigeon | Inferring min and max invariants using max-plus polyhedra[END_REF]. Non-linear relationships between state variables can be encoded by ellipsoids [START_REF] Xavier Allamigeon | A scalable algebraic method to infer quadratic invariants of switched systems[END_REF][START_REF] Alexander | Ellipsoidal techniques for reachability analysis[END_REF], polynomial templates [START_REF] Adjé | Coupling policy iteration with semi-definite relaxation to compute accurate numerical invariants in static analysis[END_REF][START_REF] Rodríguez | Automatic generation of polynomial invariants of bounded degree using abstract interpretation[END_REF] and polynomial zonotopes [START_REF] Adjé | Quadratic zonotopes[END_REF][START_REF] Althoff | Reachability analysis of nonlinear systems using conservative polynomialization and non-convex sets[END_REF].

Similarly, a barrier certificate [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF] can be used to separate the reachable set from an unsafe set using a sub-level set of a suitable function, such as a polynomial.

The efficiency of a set representation in verification depends on the efficiency of computing common operations used in reachability analysis. For a linear networked control system with additive uncertainty, linear transformation and the Minkowski sum are the main operations involved in the computation of the reachable sets. The zonotope representation [START_REF] Girard | Reachability of uncertain linear systems using zonotopes[END_REF] can be very efficient for bounded time reachability of a linear networked control system because the linear transformation and Minkowski sum over zonotopes can be computed efficiently. However, approximating unbounded time reachable sets typically requires computing invariants. Efficiently computing invariants requires encoding directions for convergence of state trajectories, which can depend on complex valued eigenvectors. But real valued zonotopes can not capture the convergence along complex valued eigenvectors. Therefore, we generalize zonotopes to complex zonotopes that can incorporate possibly complex eigenvectors among its generators in a way that the reachable set approximation can converge along the complex eigenvectors. This generalization of zonotope to complex zonotope can express some non-polytopic sets in addition to polytopic zonotopes in the real domain, and hence more expressive. This generalization is similar in spirit to quadratic zonotope [START_REF] Adjé | Quadratic zonotopes[END_REF] or more generally polynomial zonotope [START_REF] Althoff | Reachability analysis of nonlinear systems using conservative polynomialization and non-convex sets[END_REF] since both represent constraints on variables. However, while complex zonotope involves a linear combination of some ellipsoids and line segments, polynomial zonotope involves a non-linear combination of only line segments. Therefore, both represent very different classes of non-polyhedral sets. Moreover, we show in this paper that an infinite parametrized family of invariant complex zonotopes can be represented efficiently for any stable linear system using eigenvectors. However, there is no known guarantee of existence of invariant polynomial zonotopes, except trivially the equilibrium, for stable linear systems.

We remark that not many set representations can handle additive disturbance efficiently. Ellipsoids are not closed under the Minkowski sum [START_REF] Xavier Allamigeon | A fast method to compute disjunctive quadratic invariants of numerical programs[END_REF], which can result in reduction of approximation accuracy in the presence of additive disturbance. The work of Allamigeon et. al [START_REF] Xavier Allamigeon | A fast method to compute disjunctive quadratic invariants of numerical programs[END_REF] proposed an over-approximation of the Minkowski sum of ellipsoids based on the Löwner order, still the exact Minkowski sum can not be represented by an ellipsoid. Although polytopes are closed under Minkowski sum, there can be exponential blowup of complexity of the resulting half-space representation [START_REF] Kvasnica | Minkowski addition of convex polytopes[END_REF]. In contrast, zonotope and its extension to complex zonotope can exactly represent the Minkowski sum and its computation is also very efficient.

We draw inspirations from the algorithms for stability verification of sampled data systems which compute invariants [START_REF] Al Khatib | Stability verification of nearly periodic im-pulsive linear systems using reachability analysis[END_REF][START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched lyapunov function approach[END_REF][START_REF] Fiacchini | Set theory conditions for stability of linear impulsive systems[END_REF][START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF][START_REF] Pettersson | Hybrid system stability and robustness verification using linear matrix inequalities[END_REF] either as sub-level sets of Lyapunov functions or polytopes. In our work, we extend the use invariants to verify safety in addition to stability.

Furthermore, our complex zonotope representation can efficiently handle additive disturbance due to efficient computation of the Minskowski sum. This is an advantage over sub-level sets of Lyapunov functions where the Minkowski sum can not be represented accurately, and also over polytopes whose representation is of exponential complexity [START_REF] Kvasnica | Minkowski addition of convex polytopes[END_REF].

Networked control systems

In a networked control system (NCS), a controller input is exchanged over a network between different components as information packets. The controller input is sampled at discrete time instants, while it remains constant between successive sampling times. But the sensors which estimate the output may be inaccurate and the sampling period can be uncertain, possibly due to packet dropouts. In this paper, we consider systems with linear dynamics, linear feedback input from the controller, uncertainty in sampling period and inaccurate output estimation modeled by additive error and additive disturbance input. This system can also be categorized as a hybrid system because of periodic reset of feedback input and constraints on its sampling period.

The system is modeled as follows. The state of the entire system at time

t ∈ [0, ∞) is x t ∈ R n , the feedback input exchanged between components is u t ∈ R m , disturbance input is v t ∈ V ⊆ R m , output is y t ∈ R p , additive error estimation is w t ∈ W ⊆ R p ,
τ min is a lower bound on feedback sampling period, τ max is an upper bound on sampling time, the set of initial states is Ω and

A ∈ R n×n , B ∈ R m×n , C ∈ R p×n , D ∈ R p×m , F ∈ R m×p
are real matrices related to the system dynamics described below.

∃ (t k ) ∞ k=0 , ∀k ∈ Z ≥0 t k ∈ R ≥0 ∀t / ∈ ∞ k=0 {t k } , ∂x t ∂t = Ax t + Bu t , ∂u t ∂t = 0 (1) ∀k ∈ Z ≥0 , y t k = Cx t k + Du t k + w t k ( 2
)
u t k = F y t k + v t k (3) (t k+1 -t k ) ∈ [τ min , τ max ], t 0 = 0, x 0 ∈ Ω (4)
We denote the combined vector of the state of the plant and controller input at any time t as

z t = x T t , u T t T .
We call a sequence of combined states at sampling times (z t k ) ∞ k=0 as a sampling time trajectory which, by simple manipulation of the above equations, can be shown to be equivalently governed by the following dynamics.

z t k+1 = R (t k+1 -t k ) z t k + J 1 v t k+1 + J 2 w t k+1 where ∀τ ∈ [τ min , τ max ], R τ = A r exp (A c τ ) , A r = I n 0 F C F D , A c = A B 0 0 , J 1 = 0 I m , J 2 = 0 F (5)
▶ Example 1 (Damped harmonic oscillator). We consider a damped harmonic oscillator where the feedback driving force is communicated over a network. The negative feedback driving force is -x 1 where x 1 is the position of the oscillator. In this case, we we have the following matrices for the dynamics.

A = 0 1 0 -0.5 B = 0 1 C = 1 0 D = 0 F = -1
Then the equivalent matrices for the dynamics of combined states, i.e., A c , A r , J 1 , J 2 and R t for any t ∈ [0, ∞) are the following.

A c =   0 1 0 0 -0.5 1 0 0 0   A r =   1 0 0 0 1 0 -1 0 0   J 1 =   0 0 1   J 2 =   0 0 -1   R t =   1 0 0 0 1 0 -1 0 0   exp     0 t 0 0 -0.5t t 0 0 0    

Safety of NCS

Given a set of safe states S ⊆ R n+m , we say that an NCS is safe if the state of plant, controller and feedback input of every trajectory at all times lies within S. In other words, the set Γ of all reachable states defined as follows should lie within S.

Γ = x τ u τ | τ ∈ [0, ∞), (x t , u t ) t∈[0,∞) satisfies (1)-(4)
We shall show below that the safety of an NCS is guaranteed by the existence of an invariant set at sampling times that obeys certain conditions. We will later use this result to develop an algorithm for NCS safety verification.

▶ Definition 2 (Sampling time invariant set). A set Ψ of states is called a sampling time invariant if ∀t ∈ [τ min , τ max ], R t Ψ ⊕ J 1 V ⊕ J 2 W ⊆ Ψ.
▶ Theorem 3 (Relation between safety and sampling time invariant). For any S ⊆ R n+m , we have the reachable set Γ is included in S, if there exists a sampling time invariant Ψ such that the following is true

∀τ ∈ [0, τ max ], exp (A c τ ) Ψ ⊆ S (6)
Ω ⊆ Ψ (7)

Proof. Let us consider that at some sampling time t k , a state z t k ∈ Ψ. Then at the next sampling time t k+1 , we get the following based on (5):

z t k+1 ∈ R (t k+1 -t k ) Ψ ⊕ J 1 V ⊕ J 2 W
We have (t k+1 -t k ) ∈ [τ min , τ max ] according to the dynamics of NCS. Therefore, according to the Definition 2 of sampling time invariant, we get z t k+1 ∈ Ψ. This means that any state originating inside Ψ will remain inside Ψ at every sampling time. Since the initial set is contained in Ψ, i.e., Ω ⊆ Ψ, we get that for all possible trajectories of the system, the state remains within Ψ at the sampling time. Now, we have to prove that between any two sampling times, the state remains within the safe set.

Let t = t k + τ be any time point where t k is the latest sampling time before t, that is,

τ ∈ [0, τ max ].
Then the combined state reached at t k + τ is given by

z t = exp (A c τ ) z t k .
As we have shown that z t k ∈ Ψ, we get z t ∈ exp (A c τ ) Ψ. It follows from ( 6), we get z t ∈ S, which proves the theorem.

◀ ◀

According to the above theorem, the safety of NCS can be verified by finding a sampling time invariant satisfying [START_REF] Adjé | Quadratic zonotopes[END_REF]. The eigenstructure of the dynamics is closely related to sampling time invariants, which will be explained later in Theorem 7. Therefore, we introduce complex zonotope as a set representation in the next section, which enables us to use eigenvectors of the dynamics to find invariants. Complex zonotope also has other advantages, such as efficient computation of linear transformation and the Minkowski sum.

Complex zonotope

A simple zonotope is a set of points which are linear combinations of real vectors such that the combination coefficients are bounded in absolute values. Under this representation, the linear transformation and Minkowski sum can be exactly and quickly computed. This is useful for efficient bounded time reachability of a linear system as discussed in Girard et al. [START_REF] Girard | Reachability of uncertain linear systems using zonotopes[END_REF]. The mathematical definition of an simple zonotope or a real zonotope is as follows:

▶ Definition 4 (Simple/Real zonotope). Let P be a real valued matrix and c ∈ R rows(P ) be a real vector. The following is a real zonotope centered at c.

Z (P, c) = P ζ + c | ζ ∈ R cols(P ) , ∀i ∈ {1, . . . , cols (P )} |ζ i | ≤ 1
We have discussed in the previous section that unbounded time safety verification of NCS is related to computing sampling time invariants. For a stable NCS with real eigenvectors, a set contracts along its eigenvectors. By incorporating such eigenvectors among the generators of a zonotope, we can find sampling time invariants. However, when the eigenvectors are complex valued, real zonotopes may not be able to represent the directions for convergence among its generators. Therefore, we extend real zonotopes to the complex number domain to represent complex valued eigenvectors among their generators to find positive invariants. We show in a later lemma that incorporating complex eigenvectors allows representing invariants for a stable linear system with possibly complex eigenvalues.

▶ Definition 5 (Complex zonotope). Let P ∈ C rows(P )×cols(P ) be a complex matrix, c ∈ R rows(P ) be a real vector and s ∈ R cols(P ) ≥0

be a non-negative real vector. The following is a complex zonotope centered at c with template P and scale vector s.

Z (P, c, s) = P ζ + c | ζ ∈ C cols(P ) , |ζ| ≤ s
While real zonotopes are a subclass of polytopes, real projections of complex zonotopes are more general and include non-polytopic sets in addition to polytopic zonotopes. Geometrically speaking, real projections of complex zonotopes are the Minkowski sum of some embedded ellipses and line segments.

▶ Example 6. Let us consider a complex zonotope Z (P, c, s) where

P = 1 + 2ι 1 2 + ι 1 -2ι 1 2 -ι , c = 0 0 , s = 1 1 .
The real projection generated of The following theorem states that when the sampling time period is certain, there is no additive disturbance input and the sampling time transformation matrix R τmin is stable, we can find different sampling time invariant complex zonotopes by incorporating eigenvectors of the sampling time dynamics in its template and arbitrarily varying the scale factors.

Z 1 + 2ι 1 -2ι , 0 0 , 1 is the ellipse (x, y) ∈ R 2 | (x + y) 2 + (x-y) 2 4 ≤ 4 , real projection of Z 2 + ι 2 -ι , 0 0 , 1 is the ellipse (x, y) ∈ R 2 | (x -y) 2 + (x+y
▶ Theorem 7 (Invariance based on eigenstructure). Let us consider that τ min = τ max = t, V = {0}

and W = {0}. Let E ∈ C (n+m)×(n+m) contain the eigenvectors of R t as its columns whose corresponding complex eigenvalues are e ∈ C (n+m) . If 0 < ∥e∥ ∞ ≤ 1, then for any scale vector s ∈ R n ≥0 , the complex zonotope Z (E, 0, s) is a sampling time invariant.

Proof. Let x ∈ Z (E, 0, s). Based on the complex zonotope representation, there exists ζ ∈

C n , |ζ| < s such that x = Eζ. It then follows that

R t x = R t Eζ = ED (e) ζ ( 8
)
Let α = D (e) ζ. So R t x = Eα where we derive the following bound on th magnitude of α i for any i ∈ {1, . . . , n}.

|α i | = |e i ζ i | = |e i | |ζ i |
and since ∥e∥ ∞ ≤ 1, we have

|α i | ≤ |ζ i | ≤ s i .
As |α| ≤ s, we get R t x = Eα ∈ Z (E, 0, s). Since this is true for all x ∈ Z (E, 0, s), we can establish that R t Z (E, 0, s) ⊆ Z (E, 0, s). ◀ ◀ ▶ Remark (Advantage of complex zonotope over real zonotope). The above theorem is true for complex zonotopes but not for real zonotopes because real zonotopes can not have complex valued vectors as generators. Let us consider that the dynamics with constant sampling time has complex eigenstructure where the rank of eigenvectors is the total dimension n + m. The above theorem means that using complex zonotopes we can find a sampling time invariant set containing any initial set. However, using real zonotopes, we may not be able to find such a sampling time invariant set because real zonotopes capture contraction along complex eigenvectors. This is one main advantage of using complex zonotopes apart from being geometrically more expressive than real zonotopes. The below example illustrates how a complex zonotope sampling time invariant can be found by using the complex eigenvectors as generators. But a real zonotope invariant containing an initial set can not be found by using imaginary and real parts of eigenvectors and the coordinate directions as generators.

▶ as real generators and also the identity matrix. We used convex optimization (Algorithm 1) to search for a scaling factor h and center c such that the real zonotope Z (G, c, h) is a sampling time invariant containing the initial set. But our search failed to find it. Similarly, we tested with 5 uniformly randomly generated real valued templates, but all those real templates also could not be successful in finding an invariant. This illustrates that using complex zonotope containing eigenvectors as generators can let us find invariants in examples where real zonotopes may fail.

The result in Theorem 7 is only true when there is no uncertainty in sampling period. To handle a more general case where there is uncertainty in sampling period, we can incorporate eigenvectors of multiple matrices R t for different values of t ∈ [τ min , τ max ]. Then optimization can be used to find an appropriate scale factor that guarantees invariance, which we shall explain in Section 5.

As mentioned earlier, like real zonotopes, complex zonotopes are also closed under linear transformation and the Minkowski sum and they can be computed efficiently. This property is useful for efficiently representing sets of successor states of an NCS and is exploited in computing the operations used in our proposed NCS safety verification algorithm, namely Algorithm 1.

▶ Lemma 9 (Linear transformation and the Minkowski sum). Let us consider complex zonotopes Z (P, c, s) ⊆ C n and Z (Q, e, r) ⊆ C n and real matrices A, B ∈ R n×n . Then the following is true:

AZ (P, c, s) ⊕ BZ (Q, e, r) = Z AP BQ , Ac + Be, s r . ( 9
)
Proof. Let us consider x ∈ Z (P, c, s) and y ∈ Z (Q, e, r). Then there exist ζ, α ∈ C n such that |ζ| ≤ s, |α| ≤ r, x = P ζ + c and y = Qα + e. We derive the following:

Ax + By = AP ζ + Ac + BQα + Be AP BQ ζ α + (Ac + Be) .
We have

ζ α ≤ s r
. Therefore, Ax + By ∈ Z AP BQ , Ac + Be, s r . Since this is true for all x ∈ Z (P, c, s) and y ∈ Z (Q, e, r), we obtain [START_REF] Xavier Allamigeon | Inferring min and max invariants using max-plus polyhedra[END_REF].

◀ ◀

In order to verify sampling time invariance of a complex zonotope given an interval of uncertainty, we need to check inclusion of a set of complex zonotopes (obtained by applying a sequence of transformations) inside the original complex zonotope, based on the Definition 2 of invariance. In this regard, we define a set of complex zonotopes whose templates lie in the neighborhood of a given template as follows. Let us consider a real matrix with positive entries Υ ∈ R

rows(Q)×cols(Q) ≥0
where Q is a complex matrix and a real vector ρ ∈ R cols(e) where e is a real vector.

Z (Q ▽ Υ, e ▽ ρ, r) = Z Q + Q, e + u, r | |Q| ≤ Υ, |u| ≤ ρ
The following relation is a sufficient condition for checking the required inclusion which is proved later in Lemma 11.

▶ Definition 10 (Relation for checking inclusion). Let P be a complex matrix such that P T P is a square invertible matrix. Let Υ > 0 be a real matrix with only positive elements. We define the relation

Z (Q ▽ Υ, e ▽ ρ, r) ⊑ Z (P, c, s)
if all of the following conditions are verified:

∃X, ∆ ∈ C cols(P )×cols(Q) , y ∈ C cols(P ) :

P X = QD (r) , ∆ = P † ΥD (r) , (e -c) = P y, δ = P † ρ (10) ∀i ∈ {1, . . . , rows (X)} |y i | + δ i + cols(X) j=1 |X ij | + ∆ ij ≤ s i (11) ▶ Lemma 11 (Checking inclusion). If Z (Q ▽ Υ, e ▽ ρ, r) ⊑ Z (P, c, s) is true for Υ, ρ > 0,
then the subset inclusion Z (Q ▽ Υ, e ▽ ρ, r) ⊆ Z (P, c, s) is true.

Proof. Let us assume that Z (Q ▽ Υ, e ▽ ρ, r) ⊑ Z (P, c, s) is true. Hence, there exist matrices X, ∆ and complex vectors y, ρ such that all the equations in (11) are true. Let us consider any

x ∈ Z (Q ▽ Υ, e, R). Based on the definition of a complex zonotope, there exists ζ ∈ C cols (P ) such that |ζ| ≤ s and x = Q + Q ζ + e + u, Q ≤ Υ and |u| ≤ ρ. We now have to show that

x ∈ Z (P, c, s).

Let us consider a vector α ∈ C cols(ζ) such that for any i ∈ {1, . . . , cols (ζ)}, the following is true:

α i = ζ i /r i if r i > 0, α i = 0 if r i = 0 (12)
Since |ζ| ≤ r, it follows from the above definition that |α| ≤ 1 and ζ = D (r) α. Then we derive the following.

x = Q + Q ζ + e + u = Q + Q D (r) α + e + u = Q + Q D (r) α + (e -c) + c + u
and from ( 10)

x = P Xα + P † QD (r) α + y + P † u + c ( 13
)
We derive the following for any i ∈ {1, . . . , rows (X)}

Xα + P † QD (r) α + y + P † u i ≤ |y| i + P † |u| + cols(X) j=1 |X| ij + P † Q D (r) ij |α| j ≤ |y| i + δ i + cols(X) j=1 |X| ij + ∆ ij ≤ s i . ( 14
)
The above second inequality is true because |α| ≤ 1, Q ≤ Υ, ∆ = P † ΥD (r), |u| ≤ ρ and P † ρ = δ, and the last inequality is deduced from [START_REF] Xavier Allamigeon | A scalable algebraic method to infer quadratic invariants of switched systems[END_REF].

From ( 13) and ( 14), we get that x ∈ Z (P, c, s). As this is true for any x ∈ Z (Q, e, r), the inclusion Z (Q, e, r) ⊆ Z (P, c, s) is true.

◀ ◀

We can algebraically compute bounds on the real projection of a complex zonotope along any direction, that is, the support function, as follows.

▶ Lemma 12 (Computing support function). Let us consider a complex zonotope Z (P, c, s) and a vector w ∈ R cols(c) . We have the following equality:

max x∈Z(P, c, s) Re w T x = w T c + w T P s (15)
Proof. First we prove that max x∈Z(P, c, s)

w T x ≤ w T c + w T P s ( 16
)
Let us consider x ∈ Z (P, c, s). So, there exists ζ ∈ C cols(P ) such that |ζ| ≤ s and x = P ζ + c.

We derive the following:

Re w T x = w T c + Re w T P ζ ≤ w T c + w T P s (17)
The inequality in the above formula is deduced from the fact that |ζ| ≤ s. As the above is true for all x ∈ Z (P, c, s), it proves [START_REF] Daafouz | Stability analysis and control synthesis for switched systems: a switched lyapunov function approach[END_REF].

Next we prove the following:

max x∈Z(P, c, s) w T x ≥ w T c + w T P s ( 18
)
Let us consider x = P ζ + c where ζ ∈ C cols(P ) is defined as follows.

ζ i = s i if Re w T P i ≥ 0, ζ i = -s i otherwise (19)
Then we get the following:

Re w T x = Re w T P ζ + w T c = w T P s + w T c
The second equality in the above is obtained by using [START_REF] Frehse | Spaceex: Scalable verification of hybrid systems[END_REF]. This proves [START_REF] Fiacchini | Set theory conditions for stability of linear impulsive systems[END_REF]. From the inequalities ( 16) and ( 18), we get [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF].

◀ ◀

The relation for checking inclusion (( 10) and ( 11)) consists of a set of convex constraints on the variables X, ∆, s, r, e and c when the templates Q and P are fixed (constants). In fact, they constitute a class of convex constraints called second order conic constraints (SOCC) [START_REF] Boyd | Convex optimization[END_REF].

The SOCC constraints can be solved efficiently up to high numerical precision using convex optimization techniques and many solvers are available for the same [START_REF] Grant | Cvx: Matlab software for disciplined convex programming[END_REF].

Using complex zonotopes for verification

In this section, we describe an algorithm based on operations on complex zonotopes to verify safety of NCS. Our algorithm finds a complex zonotope which is a sampling time invariant and satisfies the other required condition (6) for safety. The algorithm has two parts.

1. We find a complex zonotope that is a sampling time invariant, i.e., is invariant with respect to the transformation at all sampling times t ∈ [τ min , τ max ].

2.

Next we verify that the the reachable set of the complex zonotope by continuous evolution, without reset, within the interval [0, τ max ] remains within the safe set [START_REF] Adjé | Quadratic zonotopes[END_REF].

The detailed procedure is explained in Algorithm 1 and its correctness is proved in Theorem 17. A safe set is specified by linear constraints and an open input set, disturbance input set and initial set bounded by complex zonotopes, as follows.

S = x ∈ R n+m | Hx ≤ d , V ⊆ Z (Q V , c V , s V ) , W ⊆ Z (Q W , c W , s W ) Ω ⊆ Z (Q init , c init , s init )

Finding sampling time invariant

We first fix the template of a complex zonotope based on the eigenvectors of the dynamics and then the templates of the input sets. We can also add arbitrary vectors to the template to increase precision. Next we derive a set of convex constraints on center, scale factor and other auxiliary variables of a complex zonotope such that the complex zonotope is a sampling time invariant. We solve for the scale factor and center using convex optimization.

Let us choose a positive integer k > 0 and denote by E k the matrix containing all independent unit eigenvectors of the matrices R t for all t ∈ T k where

T k = t = τ min + i τ max -τ min k | t ∈ [τ min , τ max ]
The template P of the complex zonotope consists of these eigenvectors, the templates of inputs

sets J 1 Q V , J 2 Q W , their transformations R τmax J 1 Q V , R τmax J 2 Q W . P = E k I n+m R τmax J 1 Q V J 2 Q W ( 20
)
If P contains repeated columns or zero columns, they are removed.

▶ Remark (Choice of template). Firstly, we add the eigenvectors E K of the transformation operators at various sampling times. This heuristic is based on Theorem 7 which says that incorporating eigenvectors in the complex zonotope is closely related to finding sampling time invariants. We also include the identity matrix inside the zonotope to capture bounds along the coordinate directions. Next, we get from Lemma 9 that the resulting complex zonotopes from applying the

linear transformation R t at
Q V J 2 Q W in a
transformed complex zonotope after the switching at sampling time, by Lemma 9 . Therefore, this template is also concatenated to P k . Furthermore, adding any arbitrary generator will only increase the accuracy because our optimization will adjust the scaling factors corresponding to the generators. So, we can increase the value of K to increase the accuracy of verification. Next we have to find the center and scale factor, for the given template P k such that we get a sampling time invariant. Before we describe the algorithm for this, we derive the prerequisite mathematical results.

We can expand a transformation operator R t+δ with δ > 0 in a neighborhood of sampling time t using the Taylor expansion as follows.

R t+δ = A r exp (A c t) exp (A c δ) = R t + R t A c δ + R t A c δ 2 /2 + Λ δ ( 21
)
where

|Λ δ | ≤ E t δ = |R t | |exp (A c t)| |Ac|δ 3 /3.
We use the following lemma, which is proved in [START_REF] Hetel | Lmi control design for a class of exponential uncertain systems with application to network controlled switched systems[END_REF] to bound all matrices

R t + R t A c δ + R t A c δ 2 /2 | δ ≤ ϵ
for any ϵ > 0 by the convex hull of a finite set of matrices.

▶ Lemma 14 ( [25] ). Let L 0 , L 1 , ..., L r be a finite sequence of real matrices and

U j (δ) = j i=1 L i δ i . If 0 ≤ δ < ϵ, then U r (δ) ∈ Conv (U 0 (ϵ) , . . . , U r (ϵ)).
Proof. This lemma is proved in [START_REF] Hetel | Lmi control design for a class of exponential uncertain systems with application to network controlled switched systems[END_REF]. ◀ Let us denote the set of finite matrices Θ (t, ϵ)

= j i=1 R t A c ϵ i | j = 0, 1, 2 .
Then using Lemma 14, we have the following set inclusion:

R t + R t A c δ + R t A c δ 2 /2 | δ ≤ ϵ ⊆ Conv (Θ (t, ϵ)) (22)
Therefore, using the expansion [START_REF] Girard | Zonotope/hyperplane intersection for hybrid systems reachability analysis[END_REF] and the above result [START_REF] Grant | Cvx: Matlab software for disciplined convex programming[END_REF], we get that for any ϵ > 0 and

0 ≤ δ ≤ ϵ, R t+δ ∈ Conv M ▽ E t ϵ | M ∈ Θ (t, ϵ) . ( 23
)
In order to verify sampling time invariance of Z (P, c, s), it is sufficient to divide the time interval [τ min , τ max ] into small intervals of a chosen size and verify invariance within each of the time intervals. This is explained in the following lemma.

▶ Lemma 15. If there exists N > 0 such that ∀t ∈ T N , ∀M ∈ Θ t, τ max -τ min N all of the following is true, then Z (P, c, s) is sampling time invariant. ∃ρ ∈ R cols(c) ≥0 : E t τmax-τ min N |c| ≤ ρ (24) Z   M P J 1 J 2 ▽ E t τmax-τ min N |P | 0 0 , (M c + J 1 c v + J 2 c w ) ▽ ρ,   s c v c w     ⊑ Z (P, c, s) (25)
Proof. Let us consider any τ ∈ [τ min , τ max ]. There exists t ∈ T N , δ ∈ 0, τ max -τ min N such that τ = t + δ. Then we derive the following. [START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF] and Lemma 11 for inclusion checking ⊆ Z (P, c, s) .

R τ Z (P, c, s) ⊕ J 1 Z (Q V , c V , s V ) ⊕ J 2 Z (Q W , c W , s W ) = Z R τ P J 1 Q V J 2 Q W , R τ c + J 1 c v + J 2 c W , s s v s W T % By (23), ∃M ∈ Θ t, τ max -τ min N , ∃ M ∈ 0 ▽ E t τmax-τ min N : R τ = M + M = Z M P J 1 J 2 + M 0 0 , M c + M c + J 1 c v + J 2 c w , s s v s W T % As M ∈ 0 ▽ E t τmax -τ min N and E t τmax-τ min N |c| ≤ ρ where E t τmax-τ min N ≥ 0 ⊆ Z   M P J 1 J 2 ▽ E t τmax-τ min N |P | 0 0 , (M c + J 1 c v + J 2 c w ) ▽ ρ,   s c v c w     % By
Since the above is true for any τ ∈ [τ min , τ max ], we get that Z (P, c, s) is a sampling time invariant set. ◀ ◀

Safety verification

After finding a sampling time invariant containing the initial set, we have to verify that the condition (6) that the reachable set of Z (P, c, s) by continuous evolution within [0, τ max ] is contained within the safe set. This can be verified by checking a sequence of linear inequalities as described in the following.

▶ Lemma 16. For all t ∈ [0, τ max ], we get exp (A c t) Z (P, c, s) ⊆ S if there exists ϵ > 0 such that the following are true:

∀k ∈ Z ≥0 : k ≤ τ max ϵ , ∀i ∈ {1, . . . , rows (H)} , β i (P, c, s) ≤ d where β i (P, c, s) = |H i exp (A c kϵ) P | s + H i exp (A c kϵ) c + ∥H i ∥ exp (∥A c ∥ ϵ) ∥exp (A c kϵ)∥ ϵ (∥c∥ + ∥P ∥ ∥s∥) (26)
Proof. Let us consider any t ∈ [0, τ max ]. There exists some k ∈ Z, t = kϵ + ρ, ρ ≤ ϵ. We can write

for ρ ∈ [0, ϵ], exp (A c ρ) = n i=0 A i c ρ i i = I n+m + A c ρM : M ∈ R (n+m)×(n+m) , ∥M ∥ ≤ exp (∥A c ∥ ϵ) (27)
We derive the following:

H i exp (A c t) Z (P, c, s) = H i Z (exp (A c t) P, exp (A c t) c, s)
% By Lemma 12

≤ |H i exp (A c t) P | s + H i exp (A c t) c
% By ( 27)

= |H i exp (A c kϵ) P | s + H i exp (A c kϵ) c + |H i M P | s + M c
% Substituting the bound from ( 27)

≤ |H i exp (A c kϵ) P | s + H i exp (A c kϵ) c + ∥H i ∥ exp (∥A c ∥ ϵ) ∥exp (A c kϵ)∥ ∥A c ∥ ϵ (∥c∥ + ∥P ∥ ∥s∥) ≤ d i Therefore, exp (A c t) Z (P, c, s) ⊆ {x ∈ R n+m | Hx ≤ d} = S, which proves the lemma. ◀ ◀
The following theorem summarizes the overall sufficient condition for verifying safety based on complex zonotopes, which can be checked by the procedure described in Algorithm 1.

▶ Theorem 17. We have Γ ⊆ S if there exist c ∈ R n , s ∈ R cols(P ) ≥0
, ϵ > 0 and N ∈ Z ≥0 such that all of the following conditions are true.

1. Z (Q init , c init , s init ) ⊑ Z (P, c, s).
2. (24)-( 25) are true ∀t ∈ T N , ∀M ∈ Θ (t, ϵ).

(26) is true.

Proof. By (( 24)-( 25)) and Lemma 15, we prove that Ψ = Z (P, c, s) is a sampling time invariant. By the first condition Z (Q init , c init , s init ) ⊑ Z (P, c, s), we get that the initial set Ω is contained inside the sampling time invariant Ψ. By [START_REF] Kim | Cyber-physical systems: A perspective at the centennial[END_REF] and Lemma 16, we prove that for all t ∈ [0, τ max ], exp (A c t) Z (P, c, s) ⊆ S. Then based on Theorem 3, we get Γ ⊆ S.

◀ ◀

Based on Theorem 17, we propose a semi-decision procedure in Algorithm 1 to verify safety. The algorithm is a semi-decision procedure because if it returns that the NCS is safe, then the NCS is indeed safe as proved in Theorem 17. However, the algorithm is not guaranteed to terminate. So, the user can choose to terminate the algorithm in any threshold number of iterations possibly with inconclusive result. In this context, we note NCS is a hybrid system and it is known that verification of reachability of very simple classes of hybrid systems is undecidable [START_REF] Alur | Formal verification of hybrid systems[END_REF]. So, it is unlikely that we can not come up with a sure shot decision procedure to verify safety of NCS.

▶ Example 18. Let us the NCS in Example 13. We took K = 3 and ran our algorithm using the template matrix containing complex eigenvectors as shown in Equation 20, without any random matrix, i.e., Π = 0. Then we could verify that the x 1 bounds on the unbounded time reachable

set is |x 1 | ≤ 2.22.
Comparison with real zonotope. Next, we concatenated the real and imaginary parts of P k into a real matrix G where repeated column vectors were removed. We ran our algorithm to find a real zonotope invariant containing the real template G, but our algorithm failed. We also considered 5 uniformly randomly generated real valued templates to run our algorithm, but the search for sampling time invariant failed. Thus, complex zonotope based on eigenstructure is shown to be the better choice for verification on this example.

Experimental results

We implemented our algorithms and tested them on benchmark examples of NCS. We drew comparison with simple zonotope, i.e., having real valued generators and also the state-of-the-art tool SpaceEx [START_REF] Frehse | Spaceex: Scalable verification of hybrid systems[END_REF]. For comparison with simple zonotope, we took the real template as the concatenation of real and imaginary parts of our complex template without repeated columns, and ran the same algorithm. In SpaceEx [START_REF] Frehse | Spaceex: Scalable verification of hybrid systems[END_REF], the verification is performed by step-by-step forward reachability computation. In SpaceEx, we modeled NCS with uncertain sampling time as a hybrid system with linear guards and linear transitions. For convex optimization, we used CVX (version 2.2) with MOSEK solver (version 7.1) and Matlab 2020a on a computer with 1.4 GHz Intel Core i5 processor and 4 GB 1600 MHz DDR3. The precision of the solver is set to the default precision of CVX.

Networked platoon of vehicles

This example is adapted from a model of a networked cooperative platoon of vehicles, which is presented as a benchmark in the ARCH workshop [START_REF] Ben | Networked cooperative platoon of vehicles for testing methods and verification tools[END_REF]. The platoon consists of three follower For a given uncertainty of sampling period, the goal is to find the minimum possible reference

distances d ref ∈ R 3 such that the vehicles should not collide, that is, -e i ≤ d ref i ∀i ∈ {1, 2, 3}.
For our experiment, we fix the lower bound on sampling period as 0.01(s) and verify the minimum possible bounds for various values of τ max ∈ {0.012, 0.014, 0.016, 0.018s, 0.02} s. We used our Algorithm 1 to verify the bounds taking K = 3. We then used repeated the same algorithm with simple real valued zonotope containing concatenation of real and imaginary parts of the complex template. We also used SpaceEx to verify the bounds and draw a comparison based on the smallest value of the bound verified. in sampling time interval. Also, the simple zonotope version of our algorithm could not find a sampling time invariant. On the other hand, our complex zonotope containing eigenvectors in its template is able to represent a finite invariant for various levels of uncertainty in sampling time. Furthermore, the fact that use of simple zonotope failed to compute invariant but complex zonotope containing complex eigenvectors is successful shows that use of complex eigenvectors increases the chance of finding invariant in the presence of complex eigenstructure.

Self-balancing two wheeled robot

This example concerns a verification problem for the model of a self-balancing two wheeled robot called NXTway-GS1 1 by Yorihisa Yamamoto, which was presented in the ARCH workshop [START_REF] Heinz | Benchmark: Reachability on a model with holes[END_REF].

We consider the linearized networked control system model from the paper. The state of the plant is represented by a 6-dimensional vector x p = ( θ, θ, ρ, ρ, φ, ϕ) T , where θ is the average angle of the left and right wheel, ρ is the body pitch angle, ϕ is the body yaw angle, and the rest coordinates are their respective angular velocities. The output of the plant is represented by a 3-dimensional vector ( ρout , θ m1 , θ mr ) T such that y p = C p x p . The input to the plant u p is a 2dimensional vector. The plant gets feedback input from a controller whose state is a 5-dimensional vector x c = θ err , θ ref , θref_lpf , ρ, θ lpf . The variable θ err is integration of error between θ and θ ref , and integration of θref is θ ref .

The low pass filter applied to θref is θref_lpf , body pitch angle is ρ, and the low pass filter applied to average valued of left/right motor angle is θ lpf . The output of the controller is a 2-dimensional vector which is used to compute the feedback input.

There is also a 2-dimensional unknown disturbance input.

In the benchmark paper [START_REF] Heinz | Benchmark: Reachability on a model with holes[END_REF], the output of the controller is sampled at 4ms to update the feedback input. In our experiment, we also consider the case where there is uncertainty in sampling period a possible disturbance in estimated output due to inaccurate sensors. The original model has an eleven dimensional continuous state of the plant (6) and controller (5) and 5-dimensional input. The trajectories of the system were unbounded along a 3-dimensional subspace of the system which can be found by diagonalization. We decoupled these unbounded directions and performed model reduction to obtain a lower dimensional system. For performing this decoupling, we used linear transformation of the state space based on block diagonalization. The transformed system has an eight dimensional continuous state, 2-dimensional controller output and four dimensional input. The matrices of the transformed NCS, disturbance input and output error sets are given below. The safety requirement is that the body pitch angle ρ of the robot should be bounded within

-π 2 , π 2 .
The verification problem for NCS we consider in this paper is to find the largest value of upper bound τ max on sampling time such that the system is safe, given a lower bound τ min = 0.1(ms) and the safe set ρ ∈ [-π 2 , π 2 ]. We used Algorithm 1 to verify safety bounds and subsequently find results for the above challenges based on complex zonotope. The same algorithm is repeated with a simple zonotope containing the concatenation of real and imaginary parts of our complex template. Concerning the experiment with SpaceEx using support functions, we tested with the octagon template and a template with 400 uniformly sampled support vectors distributed uniformly.

Results.

Our algorithm based on complex zonotope could verify safety for τ max = 2ms, given τ min = 0.1ms.

But SpaceEx could not find finite bounds on the reachable set for any value of uncertainty in sampling time. The simple zonotope based algorithm found a sampling time invariant, but its bounds were far over the threshold of safety, i.e., ρ ∈ [-5.1π/2, 5.1π/2]. The computation time for complex zonotope based algorithm is 153s. We increased K to 7, but could not find any larger verified sampling time interval.

Remarks.

The step-by-step reachability algorithm in SpaceEx possibly could not overapproximate resulting set of states from switching over the sampling time interval, due to which it failed to find bounds on reachable set. On the other hand, our complex zonotope containing eigenvectors in its template found bounded invariant which is within the limits of the safe set. Regarding the simple zonotope, although a sampling time invariant was found, the bounds were far over the threshold of safety.

The reason complex zonotope is far more accurate than simple zonotope on this example is possibly that complex zonotope is geometrically more expressive, being able to encode nonlinear boundaries of invariants.

We saw that increasing K did not allow verification in a larger time interval. So, the strategy of sampling more time stamps for eigenvector template may not be the best strategy for improving the accuracy of verification. The problem remains open how to select sub-matrices for concatenation to the template to improve the accuracy of verification.

Conclusion

Given the pervasiveness of networked control systems with safety-critical applications, it is essential to develop verification algorithms for such systems in the presence of various possible inaccuracies in their execution. In this paper, we developed an algorithm to verify unbounded time safety of NCS with uncertain feedback sampling period, inaccurate output sensing and disturbance input. Our algorithm uses a novel set representation called complex zonotope that can capture convergence of forward reachable sets along eigenvectors and represent invariants.

Complex zonotope is essentially an extension of simple zonotopes to the complex domain so as to efficiently compute invariants required by safety verification. Geometrically, their real projections represent a wider class of sets including some non-polytopic sets, while they retain the advantage of usual zonotope that the Minkowski sum and linear transformation can be computed efficiently.

The practicality of our algorithm is demonstrated by successfully verifying benchmark examples with high dimensions (≥ 12 state+controller input variables), which the simple zonotope and another state-of-the-art tool failed to verify. An important direction for extension of this research is verifying NCS with non-linear differential equations and feedback. In this context, we need to find complex zonotope approximation of non-linear transformations and also conditions for checking invariance under non-linear transformation.

) 2 4 ≤ 4

 44 and the real projection of Z the line segment (x, y) ∈ R 2 | -2 ≤ x + y ≤ 2 . So, the non-polytopic real projection of Z (P, c, s) in this example is the Minkowski sum of the two ellipses and one line segment as shown in Figure1.

Figure 1

 1 Figure 1 Complex zonotope as the Minkowski sum of two ellipsoids and a line segment

Example 8 .

 8 Let us consider the damped harmonic oscillator NCS in Example 1. Let us consider an initial set Ω = [-0.85, 0.85] × [0, 0]. Let us consider a complex matrix P = time invariant of the NCS such that Ω ⊆ Z (P, c, s). Now let us take the real matrix G = Re (P ) Im (P ) I 3 , which contains the real and imaginary parts of the complex matrix P

▶ Example 13 .

 13 Let us consider the NCS in Example 1. We consider an lower bound 0.1s and upper bound 0.3s on the sampling time. We consider the safe set S = z ∈ R 3 | [100]z ≤ 1 , initial set Ω = [0.85, 0.85] × [0, 0] × [0, 0], open input set V = [-0.2, 0.2] and bounds on disturbance input set W = [-0.2, 0.2]. For k = 3, we compute E k as the concatenation of eigenvectors of R 0.1 = A r exp (0.1A c ), R 0.2 = A r exp (0.2A c ) and R 0.3 = A r exp (0.3A c ) where A r and A c are given in Example 1.

vehicles M 1 ,| 1 T,

 11 M 2 and M 3 along with a leader board ahead M 4 . Each of the vehicles receives feedback input added to the acceleration, which depends on the communication of their relative distances, velocities and accelerations over a WLAN. The distance between a vehicle M i and its next vehicle M i+1 , relative to a reference distances d ref i , is denoted by e i . The acceleration of the leader vehicle is a L which ranges between [-9, 1](m/s). The state of the system is denoted by a 9-dimensional vector x = [e 1 , ė1 , ë1 , e 2 , ė2 , ë2 , e 3 , ė3 , ë3 ], which is the reference distance minus relative distances, the relative velocities and relative accelerations of the vehicles. The disturbance input is the acceleration a L ∈ [-9, 1](m/s 2 ) of the leader board. In our adaptation, we consider that there can be uncertainty in the sampling period of the feedback input. The matrices of the NCS model are given below. 3.5730 -0.0964 0.8472 3.2568 -0.0876 1.2726 3.0720 -3x ∈ [-9, 1] , i.e., Re Z 0 0

V

  6895 0.0000 0.1537 0.6846 0.0151 -0.0135 1.1696 -0.6538 -0.6895 0.0000 0.1537 0.6846 0.0151 -0.0135 1.1696 -0= [-100, 100] × [-100, 100] W = [-1, 1] 6.2.0.1 Verification problem.

6.1.0.2 Results.

	For the different values of τ max ∈ {0.012, 0.014, 0.016, 0.018, 0.02} s, we could verify finite values
	for d ref 1 , d ref 2 , d ref

3

such that -e ≤ d ref using our complex algorithm with complex zonotopes.

The verified bounds are given in

Table 1 .

 1 The computation time is less than 310s for every bound and sampling time interval below. On the other hand, the simple zonotope version of our algorithm was unsuccessful in finding finite bounds. SpaceEx terminated unsuccessfully without being able to find any bounds on the reachable set.

	Effect of increasing number of sampling times for eigenvectors. We increased K to 5 and got
	the following smaller bounds for τ max = 0.2s: d ref 1	= 40m, d ref 2	= 29m and d ref 3	= 19m. But the
	computation time increased to 656s.
	6.1.0.3 Remarks.
	In a continuous feedback model without switching based on digital sampling, the SpaceEx
	verification found the bounds d ref = (30, 30, 16) as reported in [29]. However, when there is
	intermittent switching of feedback input with uncertainty in the sampling time, SpaceEx could

not find bounds in our experiment. Possibly the inductive step-by-step reachability algorithm used in SpaceEx could not overapproximate the reachable set resulting from switching of states

Table 1

 1 Verified bounds for vehicle platoon

Algorithm 1: Verifying safety of NCS k ← 3;