Monitoring Bounded LTL Properties Using Interval Analysis

Verification of temporal logic properties plays a crucial role in proving the desired behaviors of hybrid systems. In this paper, we propose an interval method for verifying the properties described by a bounded linear temporal logic. We relax the problem to allow outputting an inconclusive result when verification process cannot succeed with a prescribed precision, and present an efficient and rigorous monitoring algorithm that demonstrates that the problem is decidable. This algorithm performs a forward simulation of a hybrid automaton, detects a set of time intervals in which the atomic propositions hold, and validates the property by propagating the time intervals. A continuous state at a certain time computed in each step is enclosed by an interval vector that is proven to contain a unique solution. In the experiments, we show that the proposed method provides a useful tool for formal analysis of nonlinear and complex hybrid systems.


Introduction
Reasoning of the temporal logic properties in a hybrid system is a challenging and important task that lies in the intersection among computer science, numerical analysis, and control theory.Various methods for falsification of hybrid systems with temporal properties have been developed, e.g., [21,20,6,27], and these methods enable verification of various properties (e.g., safety, stability, and robustness) of large and complex systems.The state-of-the-art tools are based on numerical simulations whose numerical errors often produce a qualitatively wrong result and become problematic even in a statistical evaluation.
A fundamental process in formal methods for hybrid systems is computation of rigorously approximated reachable sets.The techniques based on interval analysis

Related Work
Many previous studies have applied interval methods to reachability analysis of hybrid systems [8,5,24,15,3,11,12].The outcome of these methods is an overapproximation of a set of reachable states with a set of boxes.In interval analysis, a computation often provides a proof of unique existence of a solution within a resulting interval.This technique also applies in interval-based reachability analysis [15,14], but it is not considered in most of the methods for hybrid systems.Our method enforces the use of the proof to verify more generic temporal properties.
Reasoning of real-time temporal logic has been a research topic of interest [2,25].Numerical method for falsification of a temporal property is straightforward [16].It simulates a trajectory of a bounded length and checks the satisfiability of the negation of the property described by a bounded temporal logic.This paper presents an interval extension of this falsification method.
A tree-search method for searching witness trajectories [21], a falsification method based on a Monte-Carlo optimization technique [20], and statistical model checking methods [6,27] have been proposed.These methods have been shown their practicality in the verification of realistic nonlinear models; however, their implementations are based on numerical simulations and might suffer from numerical error.Applications of our interval method include an integration with these statistical methods to achieve both reliability and practicality.An integrated statistical and interval method was also proposed in [26] for reachability analysis.

Ishii Yonezaki Goldsztejn
Notions of robustness have been proposed to facilitate the simulation-based verification of temporal properties [9,7,20].In these works, the degree of robustness is represented as a distance between a trajectory and a region where a proposition holds.A non-robust trajectory, which is computed numerically, is likely to be inconsistent with the considered model due to numerical errors.Our method ensures a robustness rigorously by verifying that a trajectory intersects with each boundary in the state space.
There exist a few methods for model checking of temporal logic properties [23,4].[23] proposed a method specialized in stability properties, which is described as a specific form of temporal logic formula.[4] proposed a method that translates a verification problem into a reachability problem with the k-Liveness scheme, which is incomplete in general settings.Our method can be viewed as a bounded model checking method that validates a bounded temporal property by ensuring that all trajectories that emerge from an initial interval value satisfy the property.

Interval Analysis
This section introduces selected topics and techniques based on interval analysis [17,19].The techniques are used in the proposed method in Section 6.All of these definitions are naturally extended to interval vectors; an n-dimensional box (or interval vector) a is a tuple of n intervals (a 1 , . . ., a n ), and I n denotes the set of n-dimensional boxes.For a ∈ R n and a ∈ I n , we use the notation a ∈ a, which is interpreted as ∀i ∈ {1, . . ., n} a i ∈ a i .In an actual implementation, the bounds of intervals should be machine-representable floating-point numbers and other real values are rounded in the appropriate directions.

Basic Notions and Techniques
For a function f : R n → R, f : I n → I is known as an interval extension of f if and only if it satisfies the containment condition ∀a ∈ For arbitrary intervals a, b, d ∈ I, the extended division 2{d ∈ d | ∃a ∈ a ∃b ∈ b a = bd} can be implemented as follows (see Section 4.3 of [19]): In the second and third cases, when b = 0 (resp.b = 0), we set a/b and a/b as −∞ and ∞ (resp.a/b and a/b as ∞ and −∞).Given a differentiable function f (a) : R → R and a domain interval a, a root ã ∈ a of f such that f (ã) = 0 is included in the result of an interval Newton operator where â ∈ a, and f and d da f are interval extensions of f and the derivative of f .Iterative applications of the operator will converge.Let a be the result of applying the operator to a.If a ⊆ int a holds, a unique root exists in a .

ODE Integration
An initial value problem (IVP) for an ordinary differential equation (ODE) is specified by a triple (t 0 , x 0 , F ) consisting of an initial value x 0 ∈ R n at time t 0 ∈ R and a flow function F : R n → R n (assume Lipschitz continuity).Given a time interval t ∈ I and a continuous trajectory x(t) : t → R n , the satisfaction for IVP-ODEs is defined as Given t 0 ∈ I and x 0 ∈ I n , we can consider a parametric IVP-ODE (t 0 , x 0 , F ), where the initial condition is parameterized, and its satisfaction relation is defined as x, t |= (t 0 , x 0 , F ) iff ∃t 0 ∈ t 0 ∃x 0 ∈ x 0 x, t |= (t 0 , x 0 , F ).
TS t (t 0 , x 0 , F ) denotes the set of satisfied trajectories on t.
Using the tools based on the interval Taylor methods, e.g., CAPD4 and VN-ODE [18], we can obtain an interval extension X : I → I n of solution trajectories in TS t (t 0 , x 0 , F ).Given t ∈ I, such tools compute a value X(t ) by performing the stepwise integration of the flow function F from the initial time t 0 to time t .In the stepwise computation of the interval Taylor methods, the unique existence of a solution is verified for a box enclosure computed in each step based on the Picard-Lindelöf operator and Banach's fixpoint theorem.Accordingly, when an interval enclosure X(t ) (assume t ≥ t 0 ) is computed with an interval Taylor method, the following property holds: In principle, if F is Lipschitz continuous and we can assume an arbitrary precision, we obtain an arbitrary narrow interval enclosure X([t]) for t ∈ R.However, since the implementations use machine-representable real numbers, it may fail to compute an enclosure in the process that verifies the unique existence property, even with the smallest step size.
We model a hybrid system as a hybrid automaton [1].For simplicity in this paper, we consider deterministic systems, i.e., the location invariant is the negation of guard conditions and two guards do not overlap in a location.The proposed method can be extended to handle non-deterministic systems, e.g., by enumerating possible paths and computing a trajectory enclosure for each path.Definition 4.1 A hybrid automaton is a septet HA := Q, x, X, Init, {F q } q∈Q , {G q,q } q∈Q,q ∈Q , {R q,q } q∈Q,q ∈Q , that consists of the following components: • A finite set of locations Q = {q 1 , . . ., q nq }.
• A domain X ⊆ I n for the valuation of the variables.
• A set of initial values Init ⊆ {q}×X where q ∈ Q.
• A set of vector fields F q : X → X (assume Lipschitz continuity).
• A set of guards G q,q ⊆ X described by a condition of the form g(x) = 0 ∧ h(x) < 0 where g, h : Behaviors of the states σ ∈ Q×X over the timeline are formalized as trajectories.
In this work, we assume that there are no consecutive multiple discrete changes.Definition 4.2 Given an HA, an initial state (q 0 , s 0 ) ∈ Init, and a time interval t = [0, t max ] (t max ∈ R ≥0 ), a state at each time t ∈ t is determined as a trajectory (q, s), which consists of a location trajectory q : t → Q and a continuous state trajectory s : t → X.The value of the trajectory is defined recursively as follows: (q(0), s(0)) := (q 0 , s 0 ), where the relation σ 1 t − → σ 2 ∈ (Q×X)×t×(Q×X) is given by the following rules: Note that the second rule also applies for t = 0.The set of trajectories of length t max is denoted by TS tmax (HA).
When a discrete change 0 − → is applied at time t, the state (q(t), s(t)) overwrites the state σ before the discrete change.An HA has a unique trajectory (q, s) starting from an initial state (q 0 , s 0 ) ∈ Init because we have assumed that two guards do not hold simultaneously.
Variables x 1 , x 2 , and x 3 represent the height and velocity of the ball, and the (global) time, respectively.Air resistance is considered in the dynamics.The height of the table sinusoidally oscillates within [−1, 1] and is represented as sin x 3 .The second proposition of the guard is to forbid the guard to hold right after a discrete change.Possible trajectories of x 1 and x 2 are illustrated in Figure 1.

Bounded Linear Temporal Logic
We consider a fragment [16] of the real-time metric temporal logic [2] such that the temporal modalities are bounded by an interval t = [t, t] such that the bounds t, t are in Q.We refer to the logic bounded linear temporal logic (BLTL) as in [27].Definition 5.1 We consider constraints in the real domain as atomic propositions.The syntax of the BLTL formulae is defined by the grammar where p belongs to a set of atomic propositions AP ϕ , U t is the "until" operator bounded with a non-empty positive time interval t ∈ I, x = (x 1 , . . ., x n ) is a vector of variables, and f : R n → R. We use the standard abbreviations, e.g., ϕ 1 ∧ ϕ 2 := ¬(¬ϕ 1 ∨ ¬ϕ 2 ), F t ϕ := true U t ϕ ("eventually"), and G t ϕ := ¬F t ¬ϕ ("always").An equation f (x) = 0 can be encoded as f (x) ≤ 0 ∧ −f (x) ≤ 0.

Semantics
The necessary length ||ϕ|| of trajectories for checking a formula ϕ is inductively defined by the structure of the formula: Let (q, s) be a trajectory in TS ||ϕ|| (HA) and ϕ be a BLTL property.We have a satisfaction relation defined as follows: ϕ 1 U t ϕ 2 intuitively means that (assuming we are at time t) ϕ 2 will hold within the time interval t+t and ϕ 1 always hold until then.We also have a validation relation defined as:

Method for Monitoring BLTL Formulae
Our interval method is based on the method proposed in [16] that decides whether a trajectory satisfies a BLTL property.In this section, we explain this basic method.First, we introduce the notion of consistent time intervals against BLTL formulae.
Definition 5.2 Let (q, s) be a trajectory of length t max and ϕ be a BLTL formula.
We say that a left-closed and right-open interval The satisfiability of a property ϕ by a trajectory is checked as follows: (i) For each atomic proposition p in ϕ, monitor the trajectory of length ||ϕ|| and identify a non-overlapping set of consistent time intervals T p = {t 1 , . . ., t np }.
(ii) Following the parse tree of ϕ in a bottom-up fashion, compute a set of consistent time intervals of ϕ.For each construct of BLTL, compute as follows: (iii) Check whether the smallest time interval in T ϕ contains time 0. If yes, ϕ is satisfied; otherwise, it is not satisfied.

Example 5.3 We verify the property
for the model in Example 4.3.Computation with the monitoring method (which is extended to an interval method) is illustrated in Figure 1.

Interval-Based Simulation and Monitoring Method
In this section, we propose an interval extension of the monitoring method in Section 5.2.
In Step (i) of the method, given an HA and a BLTL property ϕ, we first simulate the HA for length ||ϕ||.Our simulation method computes an over-approximation of a trajectory of HA in which the existence of a unique trajectory is verified, i.e., it verifies the property ∀(q 0 , x 0 ) ∈ Init ∃unique (q, s) ∈ TS ||ϕ|| (HA) q(0) = q 0 ∧ s(0) = x 0 meaning that for each initial value, there exists a unique trajectory of length ||ϕ||.
Next, our method identifies the time intervals that are consistent with an atomic proposition in ϕ (Definition 5.2).A consistent time interval [t, t) is, in general, not representable in an actual implementation; therefore, we approximate it by a pair of intervals u and u such that each of them encloses the boundaries t or t.Given an atomic proposition f (x) • 0 (• ∈ {<, ≤}) and a trajectory s, we search for a boundary enclosure u such that f (s(u)) 0 where f is an interval extension of f .
In Step (ii), the set of consistent time intervals is updated to be consistent with ϕ.This computation requires u to contain a unique boundary point, and thus a naive over-approximation ũ does not suffice because an interval extension f (ũ) may contain 0 although f (ũ) does not contain a boundary or contains several boundaries.Thanks to interval techniques, our method verifies the unique existence of a boundary in u.Finally, as an extension of the above property, our method verifies ∀(q 0 , x 0 ) ∈ Init ∃unique (q, s) ∈ TS ||ϕ|| (HA) q(0) = q 0 ∧ s(0) = x 0 ∧ s, 0 |= ϕ.
The proposed method has some limitations.First, it is a semi-decision procedure that may output an inconclusive result (unknown) because of a failure in the verification of unique existence; both the procedures for enclosing a continuous trajectory and enclosing a time where a discrete change occurs may cause errors.However, this mechanism is valuable in terms of reliability and complexity of the problem; a non-robust trajectory and a zeno HA will be rejected as an error in the verification process.In practice, when addressing a nonlinear HAs, the method may only work successfully with a sufficiently small subset Init ⊂ Init of initial values.In this way, the method can be still used for sat/invalid checking.Second, the method is a bounded model-checking method in the sense that the domain X of the variables is bounded, and it assumes a bounded length and a number of discrete changes in a trajectory.

Main Algorithm
Given an HA and a BLTL property ϕ, the proposed Monitor algorithm (Figure 2) outputs the following results: valid that implies HA |= ϕ; unsat that implies HA |= ¬ϕ; or unknown when the computation is inconclusive.
An iteration of the outmost loop corresponds to a continuous phase of the trajectory and a discrete change.At Lines 4-11, each guard for a possible transition is evaluated.The SearchZero algorithm is described in Section 6.4 and will return a time interval t c within which the guard is satisfied or ∅ if no state satisfies the guard.Next, the algorithm attempts to decide the earliest time interval by checking whether t c is strongly less than t c .If two guard crossings are too close, so two crossing time intervals overlap, the algorithm returns unknown.At Lines 12-19, for each atomic proposition of ϕ of the form f • 0 (• ∈ {<, ≤}), the algorithm searches for a boundary where the sign of f changes.Because several boundaries can exist in a continuous phase, the inner loop searches for all of them.The detected time intervals are saved in the set T associated with the atomic propositions.At Lines 20-21, the discrete change between the locations q and q is computed by evaluating an Ishii Yonezaki Goldsztejn interval extension of R q,q (X(t c )).A jump of state might switch the state of an atomic proposition; if such a switch exists, the Jump procedure should detect and record it in T .Finally, boundary points (which are enclosed by intervals) of the consistent time intervals saved in T are analyzed by the AnalIntervals procedure (Line 24, Section 6.2).The procedures X (Section 3.2) and SearchZero (Section 6.4) may results in errors.These errors are caught by the catch clause at Line 22.

Evaluation of BLTL Properties
The BLTL evaluation explained in Section 5.2 can be implemented as a rigorously approximated procedure AnalIntervals.
First, we approximate a set of consistent time intervals Next, the evaluation on the set of time intervals in Step (ii) is extended to address the approximated sets.The set of the inverted time intervals for ¬ϕ can be represented as ; otherwise, the evaluation results in unknown if u 1 0. The union of two sets of time intervals for ϕ 1 ∨ ϕ 2 can be implemented as a merge and sort process of the two approximated sets.The Shift t procedure for ϕ 1 U t ϕ 2 can be implemented as translations of the time intervals u i and u i for t and t, respectively.Some more case analyses should be applied, e.g., when the intervals in T ϕ become redundant or when they overlap.
Finally, we obtain T ϕ and conclude that ϕ is valid if

Computation of a Continuous Trajectory
If the system is in a location q ∈ Q and the value of the state variables is x ∈ I n at time t ∈ I, then the subsequent continuous evolution is specified by an IVP-ODE (t, x, F q ).Next, we can obtain an interval extension X : I → I n of the continuous trajectories in TS [t,||ϕ||] (t, x, F q ) as described in Section 3.2.

Evaluation of Boundary Conditions
Guards and atomic propositions can be treated as boundary conditions in the state space X ⊆ R n of the form where g : X → R and h : X → R. We propose the SearchZero algorithm shown in Figure 3 for searching the intersection between a trajectory and a boundary condition.Inputs to the algorithm consist of an interval extension of the continuous trajectory X : I → I n , a vector field of the current location F : X → X, the boundary condition B(x), and a time interval t init ∈ I to be searched.SearchZero searches for Next, at Lines 5 and 6, the interval Newton is applied.The extended division (Section 3) is used to implement the interval Newton to handle the numerator intervals d g , d h containing zero.Because we expand the interval Newton on the lower bound t and the extended division encloses the values in the domain t − t, the resulting t is filtered its inconsistent portion, without losing the solutions or being expanded.When the interval Newton results in ∅, SearchZero also returns ∅ to signal the unsatisfiability.At Line 9, because t may contain several solutions, t is reset to the lower bound as a starting value to compute an enclosure of the earliest solution.At Lines 10-18, SearchZero applies the interval Newton method with the inclusion test to prove the unique existence of a solution within the contracted interval t .This interval Newton verification is repeated with an inflation process of the time interval (see [13] for a detailed implementation).When reaching Line 19 with no error, the time interval t is a sharp enclosure of the first zero of g(s(t)) = 0.It remains to check that the inequality constraint h(s(t)) < 0 is satisfied inside t.
When SearchZero is implemented with machine representable real numbers, or when there is a tangency between the trajectory and the guard constraint, a computation may result in an error.Line 12 of SearchZero may give rise to error if the derivative on an (inflated) time interval contains zero.At Line 17, we limit the number of iterations according to whether the inflation ratio reaches the threshold as proposed in [13].

Experiments
We have implemented the proposed method and experimented on two examples to confirm the effectiveness of the method.Experiments were run using a 2.4GHz Intel Core i5 processor with 16GB of RAM.

Implementation
We have implemented the proposed algorithms in Figures 2 and 3 in OCaml and C/C++.The CAPD library was used for solving ODEs.Parameters t min , , and θ should be configured.Each parameter corresponds to the smallest integration step size that CAPD can take, the threshold used in Figure 3, or a threshold used in Inflate.In the experiments, these parameters were set as t min := 10 −14 , := 10 −14 , and θ := 0.01.
For most of models, our implementation only accepts small interval values, as reported in the next section, because an interval enclosure of the state after a continuous evolution and a jump expands by quite a large amount, and thus, the verification process in SearchZero or the solving process of CAPD will fail.

Bouncing ball
Example 5.3 can be verified using the implementation by limiting the initial value of x 2 to a small interval of width at most 0.01.We verified the model with three configurations: 1) with setting a point initial value to x 2 ; 2) with an interval initial value of width 0.01; and 3) with a point initial value and the property in which the time bound for the G operator was set to [0, 100].Each experiment was run 1000 times with the initial values of x 2 randomly picked within [2,7].The results are shown in Table 1.In the table, the second column represents the width of the For each r ∈ {valid, unsat, unknown}, the column "# r" represents the number of runs that resulted in r.The column "# errors" represents that how many of unknown results are caused by an error in the SearchZero procedure.The column "time" shows average timings.The computation of each experiment was quite efficient.
Regarding the rate of inconclusive runs in each experiment, all the verifications succeeded in the first experiment.Because we set the point initial values and the bounded simulation time, considered trajectories were always enclosed with tight intervals, and thus the verification process succeeded even in a situation that was close to singular.The result also implied that we did not meet a zeno behavior.Contrastingly, around 90% and 80% of the runs in the second and third experiments, respectively, resulted in errors.Our verification process with the interval Newton failed more often if a trajectory and a guard approached or they became close to tangent.Because the model was chaotic, a coarser enclosure of states or a longer simulation time increased the possibility to meet such situations.
In the second and third experiments, most of the successful runs resulted in valid.These runs became stable (i.e., there are less difference between the continuous trajectories in each step) in the later steps, and thus satisfied the property ϕ.We confirmed that most of the inconclusive runs fell into zeno behaviors.
It was quite rare that the result of AnalIntervals became unknown because there were always a few (or no) tight boundaries in T ϕ and they rarely contained zero.
Next, we experimented with dReal (version 2.14.08), a solver for δ-weakened SMT problems, for comparison.We consider a portion of the model of (another instance of) the bouncing ball such that the initial state is (1.1, 0, 0) and the trajectories of the ball and the table become close to tangent (Figure 4 (a)).Next, ∧ F [10,20] (−(x i − cx) 2 − (y i − cy) 2 + 100 < 0) we analyzed this model by simulating the underlying HA with our method and by solving the SMT problem with dReal, respectively.Figure 4 (b) and (c) show the computed witness trajectories.Our implementation verified the occurrence of the first contact with the floor.dReal computed two enclosures for the possible trajectories of the model; they seemed corresponding to the first and last intersections with the guard.The second witness seemed wrong because the guard condition became δ-sat around t = 0.5 with δ = 0.001.The computation of a boundary can be troublesome in this way, without the unique existence verification process.

Air Traffic Maneuver
We performed a verification of a simplified model of an air traffic maneuver (ATM) [22] in which the number of aircrafts was parameterized.A trajectory of (x i , y i ) when m = 4 is illustrated in Figure 5.We verified the following property shown in Figure 5.The first line describes that the distance between each pair of aircrafts is larger than the threshold 8/m during ||ϕ|| = 30 time units.The following of the property describes that all aircrafts reach within the circle with radius 5 within the time interval [0, 10], stay there at least 10 time units, and reach outside the circle with radius 10 after another 10 time units.We verified 10 times for each instance with m = 2, 4, 6, 8.In each verification, we randomly picked a point initial value.All runs resulted in valid.The specification of the instances and the results are shown in Table 2.The columns "# vars" and "# APs" represent the numbers We have presented a sound BLTL validation method that assures that all initialized trajectories satisfy the property.The proposed method is able to detect a witness trajectory that is verified its unique existence with an interval-based ODE integration and an interval Newton method.We consider the experimental results are promising for the practical use.
In future work, our method and implementation should be improved to allow large and uncertain initial values.Examples in a realistic setting should be demonstrated with the implementation.
a} and I denotes the set of intervals.For an interval a, a and a denote the lower and upper bounds; the width is defined as a − a; and int a denotes the interior {b ∈ R | a < b < a}.[a] denotes a point interval [a, a].For intervals a and b, d(a, b) denotes the hypermetric between the two, i.e., max(|a − b|, |a − b|).For a set S ⊂ R, 2S denotes the interval [inf S, sup S].

Figure 1 .
Figure 1.Verification of the bouncing ball example Example 4.3 We model a bouncing ball on a moving table as an HA:

Figure 4 .
Figure 4. Results of boundary detection of the bouncing ball example

Figure 5 .
Figure 5.A trajectory of ATM (left) and the BLTL property to verify (right)

Table 2
Experimental results (ATM) of variables in HA and APs in the property.The average timings rose exponentially as m increased.