Taming Delays in Dynamical Systems Unbounded Veriﬁcation of Delay Differential Equations

. Delayed coupling between state variables occurs regularly in technical dynamical systems, especially embedded control. As it consequently is omnipresent in safety-critical domains, there is an increasing interest in the safety veriﬁcation of systems modelled by Delay Differential Equations (DDEs). In this paper, we leverage qualitative guarantees for the existence of an exponentially decreasing estimation on the solutions to DDEs as established in classical stability theory, and present a quantitative method for constructing such delay-dependent estimations, thereby facilitating a reduction of the veriﬁcation problem over an unbounded temporal horizon to a bounded one. Our technique builds on the linearization technique of nonlinear dynamics and spectral analysis of the linearized counterparts. We show experimentally on a set of representative benchmarks from the literature that our technique indeed extends the scope of bounded veriﬁcation techniques to unbounded veriﬁcation tasks. Moreover, our technique is easy to implement and can be combined with any automatic tool dedicated to bounded veriﬁcation of DDEs.


Introduction
The theory of dynamical systems featuring delayed coupling between state variables dates back to the 1920s, when Volterra [41,42], in his research on predator-prey models and viscoelasticity, formulated some rather general differential equations incorporating the past states of the system.This formulation, now known as delay differential equations (DDEs), was developed further by, e.g., Mishkis [30] and Bellman and Cooke [2], and has witnessed numerous applications in many domains.Prominent examples include population dynamics [25], where birth rate follows changes in population size with a delay related to reproductive age; spreading of infectious diseases [5], where delay is induced by the incubation period; or networked control systems [21] with their associated transport delays when forwarding data through the communication network.These applications range further to models in optics [23], economics [38], and ecology [13], to name just a few.Albeit resulting in more accurate models, the presence of time delays in feedback dynamics often induces considerable extra complexity when one attempts to design or even verify such dynamical systems.This stems from the fact that the presence of feedback delays reduces controllability due to the impossibility of immediate reaction and enhances the likelihood of transient overshoot or even oscillation in the feedback system, thus violating safety or stability certificates obtained on idealized, delay-free models of systems prone to delayed coupling.
Though established automated methods addressing ordinary differential equations (ODEs) and their derived models, like hybrid automata, have been extensively studied in the verification literature, techniques pertaining to ODEs do not generalize straightforwardly to delayed dynamical systems described by DDEs.The reason is that the future evolution of a DDE is no longer governed by the current state instant only, but depends on a chunk of its historical trajectory, such that introducing even a single constant delay immediately renders a system with finite-dimensional states into an infinite-dimensional dynamical system.There are approximation methods, say the Padé approximation [39], that approximate DDEs with finite-dimensional models, which however may hide fundamental behaviors, e.g.(in-)stability, of the original delayed dynamics, as remarked in Sect.5.2.2.8.1 of [26].Consequently, despite well-developed numerical methods for solving DDEs as well as methods for stability analysis in the realm of control theory, hitherto in automatic verification, only a few approaches address the effects of delays due to the immediate impact of delays on the structure of the state spaces to be traversed by state-exploratory methods.
In this paper, we present a constructive approach dedicated to verifying safety properties of delayed dynamical systems encoded by DDEs, where the safety properties pertain to an infinite time domain.This problem is of particular interests when one pursues correctness guarantees concerning dynamics of safety-critical systems over a long run.Our approach builds on the linearization technique of potentially nonlinear dynamics and spectral analysis of the linearized counterparts.We leverage qualitative guarantees for the existence of an exponentially decreasing estimation on the solutions to DDEs as established in classical stability theory (see, e.g., [2,19,24]), and present a quantitative method to construct such estimations, thereby reducing the temporally unbounded verification problems to their bounded counterparts.
The class of systems we consider features delayed differential dynamics governed by DDEs of the form ẋ (t) = f (x (t) , x (t − r 1 ) , . . ., x (t − r k )) with initial states specified by a continuous function φ (t) on [−r max , 0] where r max = max{r 1 , . . ., r k }.It thus involves a combination of ODE and DDE with multiple constant delays r i > 0, and has been successfully used to model various real-world systems in the aforementioned fields.In general, formal verification of unbounded safety or, dually, reachability properties of such systems inherits undecidability from similar properties for ODEs (cf.e.g., [14]).We therefore tackle this unbounded verification problem by leveraging a stability criterion of the system under investigation.

Contributions.
In this paper, we present a quantitative method for constructing a delaydependent, exponentially decreasing upper bound, if existent, that encloses trajectories of a DDE originating from a certain set of initial functions.This method consequently yields a temporal bound T * such that for any T > T * , the system is safe over [−r max , T ] iff it is safe over [−r max , ∞).For linear dynamics, such an equivalence of safety applies to any initial set of functions drawn from a compact subspace in R n ; while for nonlinear dynamics, our approach produces (a subset of) the basin of attraction around a steady state, and therefore a certificate (by bounded verification in finitely many steps) that guarantees the reachable set being contained in this basin suffices to claim safety/unsafety of the system over an infinite time horizon.Our technique is easy to implement and can be combined with any automatic tool for bounded verification of DDEs.We show experimentally on a set of representative benchmarks from the literature that our technique effectively extends the scope of bounded verification techniques to unbounded verification tasks.

Related Work.
As surveyed in [14], the research community has over the past three decades vividly addressed automatic verification of hybrid discrete-continuous systems in a safety-critical context.The almost universal undecidability of the unbounded reachability problem, however, confines the sound key-press routines to either semi-decision procedures or even approximation schemes, most of which address bounded verification by computing the finite-time image of a set of initial states.It should be obvious that the functional rather than state-based nature of the initial condition of DDEs prevents a straightforward generalization of this approach.
Prompted by actual engineering problems, the interest in safety verification of continuous or hybrid systems featuring delayed coupling is increasing recently.We classify these contributions into two tracks.The first track pursues propagation-based bounded verification: Huang et al. presented in [21] a technique for simulation-based timebounded invariant verification of nonlinear networked dynamical systems with delayed interconnections, by computing bounds on the sensitivity of trajectories to changes in initial states and inputs of the system.A method adopting the paradigm of verificationby-simulation (see, e.g., [9,16,31]) was proposed in [4], which integrates rigorous error analysis of the numeric solving and the sensitivity-related state bloating algorithms (cf.[7]) to obtain safe enclosures of time-bounded reachable sets for systems modelled by DDEs.In [46], the authors identified a class of DDEs featuring a local homeomorphism property which facilitates construction of over-and under-approximations of reachable sets by performing reachability analysis on the boundaries of the initial sets.Goubault et al. presented in [17] a scheme to compute inner-and outer-approximating flowpipes for DDEs with uncertain initial states and parameters using Taylor models combined with space abstraction in the shape of zonotopes.The other track of the literature tackles unbounded reachability problem of DDEs by taking into account the asymptotic behavior of the dynamics under investigation, captured by, e.g., Lyapunov functions in [32,47] and barrier certificates in [35].These approaches however share a common limitation that a polynomial template has to be specified either for the interval Taylor models exploited in [47] (and its extension [29] to cater for properties specified as bounded metric interval temporal logic (MITL) formulae), for Lyapunov functionals in [32], or for barrier certificates in [35].Our approach drops this limitation by resorting to the linearization technique followed by spectral analysis of the linearized counterparts, and furthermore extends over [47] by allowing immediate feedback (i.e.x(t)) as well as multiple delays in the dynamics), to which their technique does not generalize immediately.In contrast to the absolute stability exploited in [32], namely a criterion that ensures stability for arbitrarily large delays, we give the construction of a delaydependent stability certificate thereby substantially increasing the scope of dynamics amenable to stability criteria, for instance, the famous Wright's equation (cf.[44]).Finally, we refer the readers to [34] and [33] for related contributions in showing the existence of abstract symbolic models for nonlinear control systems with time-varying and unknown time-delay signals via approximate bisimulations.

Problem Formulation
Notations.Let N, R and C be the set of natural, real and complex numbers, respectively.Vectors will be denoted by boldface letters.For z = a + ib ∈ C with a, b ∈ R, the real and imaginary parts of z are denoted respectively by R(z) = a and is the modulus of z.For a vector x ∈ R n , x i refers to its i-th component, and its maximum norm is denoted by x = max 1≤i≤n |x i |.We define for δ > 0, B(x, δ) = {x ∈ R n | x − x ≤ δ} as the δ-closed ball around x.The notation • extends to a set X ⊆ R n as X = sup x∈X x , and to an m × n complex-valued matrix A as A = max 1≤i≤m n j=1 |a ij |.X is the closure of X and ∂X denotes the boundary of X.For a ≤ b, let C 0 ([a, b], R n ) denote the space of continuous functions from [a, b] to R n , which is associated with the maximum norm f = max t∈[a,b] f (t) .We abbreviate C 0 ([−r, 0], R n ) as C r for a fixed positive constant r, and let C 1 consist of all continuously differentiable functions.Given f : [0, ∞) → R a measurable function such that f (t) ≤ ae bt for some constants a and b, then the Laplace transform L{f } defined by L{f }(z) = ∞ 0 e −zt f (t) dt exists and is an analytic function of z for R(z) > b.
Delayed Differential Dynamics.We consider a class of dynamical systems featuring delayed differential dynamics governed by DDEs of autonomous type: where x is the time-dependent state vector in R n , ẋ denotes its temporal derivative dx/dt, and t is a real variable modelling time.The discrete delays are assumed to be ordered as r k > . . .> r 1 > 0, and the initial states are specified by a vector-valued function φ ∈ C r k .Suppose f is a Lipschitz-continuous vector-valued function in C 1 R (k+1)n , R n , which implies that the system has a unique maximal solution (or trajectory) from a given initial condition φ ∈ C r k , denoted as ξ φ : [−r k , ∞) → R n .We denote in the sequel by f x = ∂f ∂x1 • • • ∂f ∂xn the Jacobian matrix (i.e., matrix consisting of all firstorder partial derivatives) of f w.r.t. the component x (t).Similar notations apply to components x (t − r i ), for i = 1, . . ., k.
Example 1 (Gene regulation [12,36]).The control of gene expression in cells is often modelled with time delays in equations of the form where the gene is transcribed producing mRNA (x 1 ), which is translated into enzyme x 2 that in turn produces another enzyme x 3 and so on.The end product x n acts to repress the transcription of the gene by ġ < 0. Time delays are introduced to account for time involved in transcription, translation, and transport.The positive β j 's represent decay rates of the species.The dynamic described in Eq. ( 2) falls exactly into the scope of systems considered in this paper, and in fact, it instantiates a more general family of systems known as monotone cyclic feedback systems (MCFS) [28], which includes neural networks, testosterone control, and many other effects in systems biology.
Lyapunov Stability.Given a system of DDEs in Eq. ( 1), suppose f has a steady state (a.k.a., equilibrium) at x e such that f (x e , . . ., x e ) = 0 then x e is said to be Lyapunov stable, if for every > 0, there exists δ > 0 such that, if φ − x e < δ, then for every t ≥ 0 we have ξ φ (t) − x e < .-x e is said to be asymptotically stable, if it is Lyapunov stable and there exists δ > 0 such that, if φ − x e < δ, then lim t→∞ ξ φ (t) − x e = 0. -x e is said to be exponentially stable, if it is asymptotically stable and there exist α, β, δ > 0 such that, if φ − x e < δ, then ξ φ (t) − x e ≤ α φ − x e e −βt , for all t ≥ 0. The constant β is called the rate of convergence.
Here x e can be generalized to a constant function in C r k when employing the supremum norm φ − x e over functions.This norm further yields the locality of the above definitions, i.e., they describe the behavior of a system near an equilibrium, rather than of all initial conditions φ ∈ C r k , in which case it is termed the global stability.W.l.o.g., we assume f (0, . . ., 0) = 0 in the sequel and investigate the stability of the zero equilibrium thereof.Any nonzero equilibrium can be straightforwardly shifted to a zero one by coordinate transformation while preserving the stability properties, see e.g., [19].
Safety Verification Problem.Given X ⊆ R n a compact set of initial states and U ⊆ R n a set of unsafe or otherwise bad states, a delayed dynamical system of the form ( 1) is said to be T -safe iff all trajectories originating from any φ(t) satisfying φ and T -unsafe otherwise.In particular, we distinguish unbounded verification with T = ∞ from bounded verification with T < ∞.
In subsequent sections, we first present our approach to tackling the safety verification problem of delayed differential dynamics coupled with one single constant delay (i.e., k = 1 in Eq. ( 1)) in an unbounded time domain, by leveraging a quantitative stability criterion, if existent, for the linearized counterpart of the potentially nonlinear dynamics in question.A natural extension of this approach to cater for dynamics with multiple delay terms will be remarked thereafter.In what follows, we start the elaboration of the method from DDEs of linear dynamics that admit spectral analysis, and move to nonlinear cases afterwards and show how the linearization technique can be exploited therein.

Linear Dynamics
Consider the linear sub-class of dynamics given in Eq. ( 1): where A, B ∈ R n×n , φ ∈ C r , and the system is associated with the characteristic equation where I is the n×n identity matrix.Denote by h(z) = zI−A−Be −rz the characteristic matrix in the sequel.Notice that the characteristic equation can be obtained by seeking nontrivial solutions to Eq. ( 3) of the form ξ φ (t) = ce zt , where c is an n-dimensional nonzero constant vector.The roots λ ∈ C of Eq. ( 4) are called characteristic roots or eigenvalues and the set of all eigenvalues is referred to as the spectrum, denoted by σ = {λ | det (h(λ)) = 0}.Due to the exponentiation in the characteristic equation, the DDE has, in line with its infinite-dimensional nature, infinitely many eigenvalues possibly, making a spectral analysis more involved.The spectrum does however enjoy some elementary properties that can be exploited in the analysis.For instance, the spectrum has no finite accumulation point in C and therefore for each positive γ ∈ R, the number of roots satisfying |λ| ≤ γ is finite.It follows that the spectrum is a countable (albeit possibly infinite) set: Lemma 1 (Accumulation freedom [6,19]).Given γ ∈ R, there are at most finitely many characteristic roots satisfying R(λ) > γ.If there is a sequence {λ n } of roots of Eq. (4) Lemma 1 suggests that there are only a finite number of solutions in any vertical strip in the complex plane, and there thus exists an upper bound α ∈ R such that every characteristic root λ in the spectrum satisfies R(λ) < α.This upper bound captures essentially the asymptotic behavior of the linear dynamics: Theorem 1 (Globally exponential stability [6,36]).Suppose R(λ) < α for every characteristic root λ.Then there exists K > 0 such that where ξ φ (t) is the solution to Eq. (3).In particular, x = 0 is a globally exponentially stable equilibrium of Eq. (3) if R(λ) < 0 for every characteristic root; it is unstable if there is a root satisfying R(λ) > 0.
Theorem 1 establishes an existential guarantee that the solution to the linear delayed dynamics approaches the zero equilibrium exponentially for any initial conditions in C r .To achieve automatic safety verification, however, we ought to find a constructive means of estimating the (signed) rate of convergence α and the coefficient K in Eq. ( 5).This motivates the introduction of the so-called fundamental solution ξ φ (t) to Eq. ( 3), whose Laplace transform will later be shown to be h −1 (z), the inverse characteristic matrix, which always exists for z satisfying R(z) > max λ∈σ R(λ).
Lemma 2 (Variation-of-constants [19,36]).Let ξ φ (t) be the solution to Eq. (3).Denote by ξ φ (t) the solution that satisfies Eq. (3) for t ≥ 0 and satisfies a variation of the initial condition as φ (0) = I and φ (t) = O for all t ∈ [−r, 0), where O is the n × n zero matrix, then for t ≥ 0, Note that in Eq. ( 6), φ(t) is extended to [−r, ∞) by making it zero for t > 0. In spite of the discontinuity of φ at zero, the existence of the solution ξ φ (t) can be proven by the well-known method of steps [8].[19]).The solution ξ φ (t) to Eq. (3) with initial data φ is the fundamental solution; that is for z s.t.R(z) > max λ∈σ R(λ),
Following Theorem 2, an exponentially decreasing bound on the solution ξ φ (t) to linear DDEs of the form (3) can be assembled by computing α satisfying μ < α < 0 and the coefficient K > 0.

Identifying the Rightmost Roots
Due to the significance of characteristic roots in the context of stability and bifurcation analysis, numerical methods on identifying-particularly the rightmost-roots of linear (or linearized) DDEs have been extensively studied in the past few decades, see e.g., [3,11,43,45].There are indeed complete methods on isolating real roots of polynomial exponential functions, for instances [37] and [15] based on cylindrical algebraic decomposition (CAD).Nevertheless, as soon as non-trivial exponential functions arise in the characteristic equation, there appear to be few, if any, symbolic approaches to detecting complex roots of the equation.
In this paper, we find α that bounds the spectrum from the right of the complex plane, by resorting to the numerical approach developed in [11].The computation therein employs discretization of the solution operator using linear multistep (LMS) methods to approximate eigenvalues of linear DDEs with multiple constant delays, under an absolute error of O (τ p ) for sufficiently small stepsize τ , where O (•) is the big Omicron notation and p depends on the order of the LMS-methods.A well-developed MATLAB package called DDE-BIFTOOL [10] is furthermore available to mechanize the computation, which will be demonstrated in our forthcoming examples.

Constructing K
By the inverse Laplace transform (cf.Theorem 5.2 in [19] for a detailed proof), we have where α is the exponent associated with the bound on ξ φ (t) in Eq. ( 7), and hence by substituting z = α + iν, we have , together with the fact that an integral over a quadratic integrand is convergent, it follows that By taking the norm while observing that e iνt = 1, we get .
(8) For the integral (8-a), the fact1 that implies Notice that the second integral (8-b) is computable, since it is convergent and independent of t.The underlying computation of the improper integral, however, can be rather time-consuming.We therefore detour by computing an upper bound of (8-b) in the form of a definite integral, due to Lemma 4, which suffices to constitute an exponential estimation of ξ φ (t) while reducing computational efforts pertinent to the integration.Lemma 4.There exists M > 0 such that inequation (11) below holds for any α > μ.
Proof.The proof depends essentially on constructing a threshold M > 0 such that the integral over |ν| > M can be bounded, thus transforming the improper integral in question to a definite one.To find such an M , observe that Without loss of generality, suppose the entry of h −1 (z) at (i, j) takes the form where p ij k (•) and q k (•) are polynomials in e −rz as coefficients of z k .Since e −rz is bounded by e −rα along the vertical line z = α + iν, we can conclude that there exist

We can thus choose |ν|
where the third inequality holds since |ν| > M.
In contrast to the existential estimation guarantee established in Theorem 2, exploiting the construction of α and K gives a constructive quantitative criterion permitting to reduce an unbounded safety verification problem to its bounded counterpart: Theorem 3 (Equivalence of bounded and unbounded safety).Given X ⊆ R n a set of initial states and U ⊆ R n a set of bad states satisfying 0 / ∈ U, suppose we have α satisfying μ < α < 0 and K from Eq. (12).Let K = K 1 + B r 0 e −ατ dτ X , then there exists T * < ∞, defined as such that for any T > T * , the system (3) is ∞-safe iff it is T -safe.
Proof.The "only if" part is for free, as ∞-safety subsumes by definition T -safety.
For the "if" direction, the constructed K in Eq. ( 12) suffices as an upper bound of e −αt ξ φ (t) , and hence by Theorem 2, ξ φ (t) ≤ Ke αt for any t ≥ 0 and φ constrained by X .As a consequence, it suffices to show that T * given by Eq. ( 13) is finite, which then by definition implies that system (3) is safe over t > T * .Note that the assumption 0 / ∈ U implies that there exists a ball B(0, δ) such that B(0, δ) ∩ U = ∅.Moreover, Ke αt is strictly monotonically decreasing w.r.t.t, and thus T = max{0, ln(δ/ K)/α} is an upper bound3 of T * , which further implies T * < ∞.
Example 2 (PD-controller [17]).Consider a PD-controller with linear dynamics defined, for t ≥ 0, as which controls the position y and velocity v of an autonomous vehicle by adjusting its acceleration according to the current distance to a reference position y * .A constant time delay r is introduced to model the time lag due to sensing, computation, transmission, and/or actuation.We instantiate the parameters following [17] as κ p = 2, κ d = 3, y * = 1, and r = 0.35.The system described by Eq. ( 14) then has one equilibrium at (1; 0), which shares equivalent stability with the zero equilibrium of the following system, with ŷ = y − 1 and v = v: Suppose we are interested in exploiting the safety property of the system (15) in an unbounded time domain, relative to the set of initial states X = [−0.1,0.1] × [0, 0.1] and the set of unsafe states U = {(ŷ; v) | |ŷ| > 0.2}.Following our construction process, we obtain automatically some key arguments (depicted in Fig. 1) as α = −0.5,M = 11.9125,K = 7.59162 and K = 2.21103, which consequently yield T * = 4.80579 s.By Theorem 3, the unbounded safety verification problem thus is reduced to a T -bounded one for any T > T * , inasmuch as ∞-safety is equivalent to T -safety for the underlying dynamics.
[− Ke αt , Ke αt ] n in Eq. ( 13) can be viewed as an overapproximation of all trajectories originating from X .As shown in the right part of Fig. 1, this overapproximation, however, is obviously too conservative to be utilized in proving or disproving almost any safety specifications of practical interest.The contribution of our approach lies in the reduction of unbounded verification problems to their bounded counterparts, thereby yielding a quantitative time bound T * that substantially "trims off" the verification efforts pertaining to t > T * .The derived T -safety verification task can be tackled effectively by methods dedicated to bounded verification of DDEs of the form (3), or more generally, (1), e.g., approaches in [17] and [4].

Nonlinear Dynamics
In this section, we address a more general form of dynamics featuring substantial nonlinearity, by resorting to linearization techniques and thereby establishing a quantitative stability criterion, analogous to the linear case, for nonlinear delayed dynamics.
Consider a singly delayed version of Eq. ( 1): with f being a nonlinear vector field involving possibly non-polynomial functions.Let where f x and f y are the Jacobian matrices of f in terms of x and y, respectively; g is a vector-valued, high-order term whose Jacobian matrix at (0, 0) is O.
By dropping the high-order term g in f , we get the linearized counterpart of Eq. ( 16): which falls in the scope of linear dynamics specified in Eq. ( 3), and therefore is associated with a characteristic equation of the same form as that in Eq. ( 4).Equation ( 17) will be in the sequel referred to as the linearization of Eq. ( 16) at the steady state 0, and σ is used to denote the spectrum of the characteristic equation corresponding to Eq. (17).
In light of the well-known Hartman-Grobman theorem [18,20] in the realm of dynamical systems, the local behavior of a nonlinear dynamical system near a (hyperbolic) equilibrium is qualitatively the same as that of its linearization near this equilibrium.The following statement uncovers the connection between the locally asymptotic behavior of a nonlinear system and the spectrum of its linearization: Theorem 4 (Locally exponential stability [6,36]).Suppose max λ∈σ R(λ) < α < 0. Then x = 0 is a locally exponentially stable equilibrium of the nonlinear systems (16).In fact, there exists δ > 0 and K > 0 such that where ξ φ (t) is the solution to Eq. (16).If R(λ) > 0 for some λ in σ, then x = 0 is unstable.
Akin to the linear case, Theorem 4 establishes an existential guarantee that the solution to the nonlinear delayed dynamics approaches the zero equilibrium exponentially for initial conditions within a δ-neighborhood of this equilibrium.The need of constructing α, K and δ quantitatively in Theorem 4, as essential to our automatic verification approach, invokes again the fundamental solution ξ φ (t) to the linearized dynamics in Eq. ( 17): Lemma 5 (Variation-of-constants [19,36]).Consider nonhomogeneous systems of the form Let ξ φ (t) be the solution to Eq. (18).Denote by ξ φ (t) the solution that satisfies Eq. (17) for t ≥ 0 and satisfies a variation of the initial condition as φ (0) = I and φ (t) = O for all t ∈ [−r, 0).Then for t ≥ 0, In what follows, we give a constructive quantitative estimation of the solutions to nonlinear dynamics, which admits a reduction of the problem of constructing an exponential upper bound of a nonlinear system to that of its linearization, as being immediately evident from the constructive proof.
Proof.The existence of K follows directly from Eq. ( 7) in Theorem 2. By the variationof-constants formula (19), we have, for t ≥ 0, •) being a higher-order term yields that for any > 0, there exists δ > 0 such that x φ t ≤ δ implies g (x(t), x(t − r)) ≤ x φ t .Due to the fact that ξ φ (t) ≤ Ke αt and the monotonicity of ξ φ (t) with α < 0, we have x φ t ≤ Ke α(t−r) .This, together with Eq. ( 20), leads to By the Grönwall-Bellman inequality [1] we obtain The above constructive quantitative estimation of the solutions to nonlinear dynamics gives rise to the reduction, analogous to the linear case, of unbounded verification problems to bounded ones, in the presence of a local stability criterion.
Theorem 6 (Equivalence of safety properties).Given initial state set X ⊆ R n and bad states U ⊆ R n satisfying 0 / ∈ U. Let σ denote the spectrum of the characteristic equation corresponding to Eq. (17).Suppose that max λ∈σ R(λ) < α < 0, and the fundamental solution to Eq. (17) satisfies ξ φ (t) ≤ Ke αt for any t ≥ 0. Let K = Ke −rα 1 + B r 0 e −ατ dτ X .Then there exists δ > 0 and T * < ∞, defined as Proof.The proof is analogous to that of Theorem 3, particularly following from the local stability property stated in Theorem 5.
Note that for nonlinear dynamics, the equivalence of safety claimed by Theorem 6 holds on the condition that X ≤ δ, due to the locality stemming from linearization.In fact, such a set B ⊆ R n satisfying B ≤ δ describes (a subset of) the basin of attraction around the local attractor 0, in a sense that any initial condition in B will lead the trajectory eventually into the attractor.Consequently, for verification problems where X ⊇ B, if the reachable set originating from X is guaranteed to be subsumed within B in the time interval [T − r, T ], then T + T * suffices as a bound to avoid unbounded verification, namely for any T > T + T * , the system is ∞-safe iff it is T -safe.This is furthermore demonstrated by the following example.[4,25]).Consider a slightly version of the delayed logistic equation introduced by G. Hutchinson in 1948 (cf.[22])

Example 3 (Population dynamics
which is used to model a single population whose percapita rate of growth Ṅ (t)/N (t) depends on the population size r time units in the past.This would be a reasonable model for a population that features a significant minimum reproductive age or depends on a resource, like food, needing time to grow and thus to recover its availability.
If we change variables, putting u = N − 1, then Eq. ( 21) becomes the famous Wright's equation (see [44]): The steady state N = 1 is now u = 0. We instantiate the verification problem of Eq. ( 22) over [−r, ∞) as X = [−0.2,0.2], U = {u | |u| > 0.6}, under a constant delay r = 1.Note that delay-independent Lyapunov techniques, e.g.[32], cannot solve this problem, since Wright's conjecture [44], which has been recently proven in [40], together with corollaries thereof implies that there does not exist a Lyapunov functional guaranteeing absolute stability of Eq. ( 22) with arbitrary constant delays.To achieve an exponential estimation, we first linearize the dynamics by dropping the nonlinearity Following our constructive approach, we obtain automatically for Eq. ( 23) α = −0.3(see the left of Fig. 2), M = 2.69972, K = 3.28727, and thereby for Eq. ( 22) δ = 0.00351678, K = 0.0338039 and T * = 0 s.It is worth highlighting that by the bounded verification method in [17], with Taylor models of the order 5, an overapproximation Ω of the reachable set w.r.t.system (22) over the time interval [14.5, 15.5] was verified to be enclosed in the δ-neighborhood of 0, i.e., Ω ≤ δ, yet escaped from this region around t = 55.3 s, and tended to diverge soon, as depicted in the right part of Fig. 2, and thus cannot prove unbounded safety properties.However, with our result of T * = 0s and the fact that Ω over [−1, 15.5] is disjoint with U, we are able to claim safety of the underlying system over an infinite time domain.
DDEs with Multiple Different Delays.Delay differential equations with multiple fixed discrete delays are extensively used in the literature to model practical systems where components coupled with different time lags coexist and interact with each other.We remark that previous theorems on exponential estimation and equivalence of safety w.r.t.cases of single delay extend immediately to systems of the form (1) with almost no change, except for replacing B e −rα with k i=1 A i e −riα and B with k i=1 A i , where A i denotes the matrix attached to x(t − r i ) in the linearization.For a slightly modified form of the variation-of-constants formula under multiple delays, we refer the readers to Theorem 1.2 in [19].Fig. 2. Left: the identified rightmost eigenvalues of h(z) and an upper bound α = −0.5 such that max λ∈σ R(λ) < α < 0; Right: overapproximation of the reachable set of the system (22) produced by the method in [17] using Taylor models for bounded verification.Together with this overapproximation we prove the equivalence of ∞-safety and T -safety of the system, for any T > (T + T * ) = 15.5 s.

Implementation and Experimental Results
To further investigate the scalability and efficiency of our constructive approach, we have carried out a prototypical implementation4 in Wolfram MATHEMATICA, which was selected due to its built-in primitives for integration and matrix operations.By interfacing with DDE-BIFTOOL 5 (in MATLAB or GNU OCTAVE) for identifying the rightmost characteristic roots of linear (or linearized) DDEs, our implementation computes an appropriate T * that admits a reduction of unbounded verification problems to bounded ones.A set of benchmark examples from the literature has been evaluated on a 3.6 GHz Intel Core-i7 processor with 8 GB RAM running 64-bit Ubuntu 16.04.All computations of T * were safely rounded and finished within 6 s for any of the examples, including Examples 2 and 3.In what follows, we demonstrate in particular the applicability of our technique to DDEs featuring non-polynomial dynamics, high dimensionality and multiple delays.

Conclusion
We have presented a constructive method, based on linearization and spectral analysis, for computing a delay-dependent, exponentially decreasing upper bound, if existent, that encloses trajectories of a DDE originating from a certain set of initial functions.We showed that such an enclosure facilitates a reduction of the verification problem over an unbounded temporal horizon to a bounded one.Preliminary experimental results on a set of representative benchmarks from the literature demonstrate that our technique effectively extends the scope of existing bounded verification techniques to unbounded verification tasks.
Peeking into future directions, we plan to exploit a tight integration of our technique into several automatic tools dedicated to bounded verification of DDEs, as well as more permissive forms of stabilities, e.g.asymptotical stability, that may admit a similar reduction-based idea.An extension of our method to deal with more general forms of DDEs, e.g., with time-varying, or distributed (i.e., a weighted average of) delays, will also be of interest.Additionally, we expect to refine our enclosure of system trajectories by resorting to a topologically finite partition of the initial set of functions.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/),which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material.If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.