Fully Automated Verification of Linear Time-Invariant Systems against Signal Temporal Logic Specifications via Reachability Analysis

Nonlinear Analysis: Hybrid Systems

SBU]Niklas Kochdumper
SBU]Stanley Bak

[SBU]organization=Stony Brook University,
city=Stony Brook (NY),

While reachability analysis is one of the most promising approaches for formal verification of dynamic systems, a major disadvantage preventing a more widespread application is the requirement to manually tune algorithm parameters such as the time step size. Manual tuning is especially problematic if one aims to verify that the system satisfies complicated specifications described by signal temporal logic formulas since the effect the tightness of the reachable set has on the satisfaction of the specification is often non-trivial to see for humans. We address this problem with a fully-automated verifier for linear systems, which automatically refines all parameters for reachability analysis until it can either prove or disprove that the system satisfies a signal temporal logic formula for all initial states and all uncertain inputs. Our verifier combines reachset temporal logic with dependency preservation to obtain a model checking approach whose over-approximation error converges to zero for adequately tuned parameters.

While we in this work focus on linear systems for simplicity, the general concept we present can equivalently be applied for nonlinear and hybrid systems.

Reachability analysis linear systems signal temporal logic formal verification model checking automation

§ INTRODUCTION

Temporal logic enjoys an increasing popularity in science and engineering, where it is for example used to specify desired behaviors for robots <cit.> and power systems <cit.>,

or to formalize traffic rules for road <cit.> and marine <cit.> traffic. However, while the expressiveness of temporal logic is on the one hand advantageous since it allows to model complex behaviors, it unfortunately also makes it very hard to check if a system satisfies a temporal logic formula.

Consequently, while automated verification of linear systems is already possible for simple specifications given by unsafe sets <cit.>, such an approach does not yet exist for the more challenging case of temporal logic specifications.
In this work we address this shortcoming with an automated verifier for linear systems, which decides whether or not the system satisfies a signal temporal logic (STL) <cit.> formula for all initial states and uncertain inputs.

§.§ State of the Art

Approaches that check if a system satisfies a temporal logic specification can be divided into the two groups runtime verification and static verification. Runtime verification, which is often realized via monitors <cit.>, treats the system as a black box and checks if the observations obtained from the system satisfy the specification. Static verification, on the other hand, considers the case where a model of the system is available and checks if all executions of the model satisfy the temporal logic specification. Since our method falls into the static verification category, we focus on this group from now on.
A standard approach for formal verification against temporal logic specifications is to convert the temporal logic formula into an equivalent acceptance automaton <cit.>.

If the system itself can be represented by a finite state automaton, formal verification reduces to checking if there exists an accepting trace for the automaton obtained by taking the automaton product of the system and the acceptance automaton for the negated temporal logic formula <cit.>, which can be realized with standard automaton analysis tools <cit.>.

Exemplary verification problem, where the initial set is shown in white with a black border, the reachable set is depicted in gray, the unsafe sets defined by the STL formula $\varphi$ are visualized in red, and exemplary trajectories of the system are shown in black. Even though the system satisfies the temporal logic formula $\varphi$, the reachset temporal logic approach <cit.> classifies the system as unsafe since it checks if the whole reachable set satisfies the predicates $x_2 < 0.5$ and $x_2 > -0.5$. This procedure produces so-called spurious traces (orange), which are not consistent with the system dynamics.

This concept carries over to dynamic systems where the behavior is described by a differential equation. Many approaches <cit.> focus on discrete-time systems since they can be transformed into an equivalent transition system: By adequately choosing the control input, discrete-time linear systems can be represented by a finite state transition system <cit.>, and therefore be verified with the same methods as finite state automata. Another approach <cit.> considers piecewise linear discrete-time systems and uses reachability analysis to refine a finite state automaton abstraction of the system. Reachability analysis is also used by a method for general discrete-time systems <cit.>, which transforms the temporal logic specification into a temporal logic tree rather than an acceptance automaton.

Also methods for continuous-time systems apply the acceptance automaton concept <cit.>. They consider temporal logic specifications given in HyLTL <cit.>, an extension of linear temporal logic <cit.> to hybrid systems, or as pattern templates <cit.>, which are close to natural language and therefore very intuitive for humans. One major disadvantage of these methods is that reachability analysis for hybrid systems, which are obtained from the automaton product of the continuous dynamics with the acceptance automaton, is quite challenging and often provides very conservative results. Apart from the acceptance automaton concept, there also exist other strategies <cit.>: An early approach <cit.> introduced TLA$^+$, a special logic for describing the behavior of hybrid systems. This enables the formal verification using a theorem prover, which however requires an analytical solution for the differential equation and is therefore only applicable for very simple systems. Other approaches <cit.> apply SMT solvers to verify hybrid systems against STL specifications. However, this technique is restricted to systems without inputs and suffers from the curse-of-dimensionality since SMT solvers split the state space. Yet another method <cit.> exploits connections between temporal logic operators and reachability problems to compute an approximation for the set of states from which an STL formula can be satisfied based on the Hamilton-Jacobi framework <cit.>. Finally, it is also possible to convert an STL formula to reachset temporal logic (RTL) <cit.>, a special logic that can be directly evaluated on the reachable set. All existing approaches have the disadvantage that parameters such as the time step size have to be tuned manually by the user for the verification to succeed.

§.§ Contribution

In this work we present the first fully-automated verifier for linear systems and signal temporal logic specifications, which is based on the reachset temporal logic approach <cit.>. One major drawback of the previous RTL method is demonstrated by the example shown in Fig. <ref>: RTL checks if the atomic predicates are satisfied by the whole reachable set, which yields very conservative results, especially if the reachable set becomes large. We solve this problem by keeping track which parts of the initial set and input set satisfy or violate the predicates at specific points in time, as it is visualized in Fig. <ref>. In summary, our verifier has the following features:

$\bullet$ The algorithm is guaranteed to converge to the solution in finite time for decidable problem instances.

$\bullet$ Our verifier is fully-automated, so no manual parameter tuning is required.

$\bullet$ In case the system violates the specification we return a counterexample in form of a falsifying trajectory, which might provide valuable insights for a system engineer.

$\bullet$ Since our verifier naturally divides the initial set and input set into parts that satisfy or violate the temporal logic formula, it can also be used to determine a safe set of states or cut away all states that yield a violation, which is very useful for controller synthesis and set-based prediction.

We will demonstrate all of these features on several interesting benchmarks.

§.§ Overview

The remainder of the paper is structured as follows:

We first introduce some preliminaries including notations, set representations as well as operations, reachable sets, and temporal logic in Sec. <ref>.
Afterward, we specify the formal problem definition in Sec. <ref>. To keep track of dependencies between reachable states and initial states as well as inputs, we require a special type of reachability analysis, which is presented in Sec. <ref>. Next, we introduce our model checking approach in Sec. <ref>, before we describe the overall verifier in Sec. <ref>. Finally, in Sec. <ref>, we demonstrate the performance of our automated verifier on several challenging benchmarks, before we provide an outlook to future directions in Sec. <ref>.

§ PRELIMINARIES

Let us first introduce the notation as well as some preliminaries and definitions.

§.§ Notation

Scalars and vectors are denoted by lowercase letters, whereas matrices are denoted by uppercase letters.
Given a vector $v \in \R^{n}$, $v_{(i)}$ represents the $i$-th entry and $\| v\|_p$ is the p-norm. Moreover, all vectors considered in this work are column vectors.
For a matrix $M \in \R^{w \times n}$, $M_{(i,\cdot)}$ refers to the $i$-th row and $M_{(\cdot,j)}$ to the $j$-th column.
The identity matrix of dimension $n$ is denoted by $I_n$ and we use $\mathbf{0}$ and $\mathbf{1}$ to represent vectors and matrices of proper dimension containing only zeros or ones. Moreover, we denote the horizontal concatenation of two matrices $M_1, M_2$ by $[M_1~M_2]$ and the vertical concatenation by $[M_1 \concat M_2]$.

The floor function $\lfloor x \rfloor$ with $x \in \R$ returns the largest integer that is smaller than $x$, and $x~\text{mod}~y$ with $x,y \in \R$ denotes the modulo operator.
While sets are denoted by standard calligraphic letters $\mathcal{S}$, interval matrices are denoted by bold calligraphic letters: $\mathbfcal{M} = [\underline{M},\overline{M}] = \{ M \in \R^{w \times n} ~ | ~ \underline{M} \leq M \leq \overline{M} \}$, where the inequality is evaluated element-wise. $||\mathbfcal{M}||_F$ denotes the Frobenius norm of an interval matrix.
Intervals are a special case of interval matrices, where the lower and upper bounds are vectors.

Tuples are denoted by bold uppercase letters: Given a $n$-tuple $\mathbf{H} = (h_1,\dots,h_n)$, $|\mathbf{H}| = n$ denotes the cardinality of the tuple and $\mathbf{H}_{(i)} = h_i$ refers to the $i$-th entry of tuple $\mathbf{H}$. Moreover, given two tuples $\mathbf{H} = (h_1, \dots ,h_n)$ and $\mathbf{K} = (k_1, \dots, k_w)$, $( \mathbf{H}, \mathbf{K} ) = (h_1, \dots ,h_n,k_1,\dots,k_w)$ denotes their concatenation, operation $\mathbf{H} \setminus \mathbf{K}$ removes all elements in $\mathbf{K}$ from $\mathbf{H}$, and the empty tuple is denoted by $\emptyset$. Finally, given a matrix $M \in \R^{w \times n}$ and a tuple $\mathbf{H} = (h_1,\dots,h_m)$ with $h_1,\dots,h_m \in \mathbb{N}_{\leq n}$, we use the shorthand notation $M_{(\cdot,\mathbf{H})} = [M_{(\cdot,h_1)} ~ \dots ~ M_{(\cdot,h_m)}]$.

§.§ Set Representations and Operations

Given sets $\mathcal{S}_1,\mathcal{S}_2 \subset \Rn{}$ and a matrix $M \in \R^{w \times n}$, we require the set operations linear map $M \mathcal{S}_{1}$, Minkowski sum $\mathcal{S}_{1} \oplus \mathcal{S}_{2}$, Cartesian product $\mathcal{S}_1 \times \mathcal{S}_2$, intersection $\mathcal{S}_1 \cap \mathcal{S}_2$, and convex hull $conv(\mathcal{S}_{1},\mathcal{S}_{2})$, which are defined as

\begin{align}
& M \mathcal{S}_{1} := \{ M s ~|~ s \in \mathcal{S}_{1} \} , \label{eq:defLinTrans} \\
& \mathcal{S}_{1} \oplus \mathcal{S}_{2} := \{ s_1 + s_2 ~|~ s_1 \in \mathcal{S}_{1}, s_2 \in \mathcal{S}_{2} \} , \label{eq:defMinSum} \\
& \mathcal{S}_{1} \times \mathcal{S}_{2} := \{ [s_1 \concat s_2] ~|~ s_1 \in \mathcal{S}_{1}, s_2 \in \mathcal{S}_{2} \} , \label{eq:defCartProd} \\
& \mathcal{S}_1 \cap \mathcal{S}_2 := \{s ~|~ s \in \mathcal{S}_1 \wedge s \in \mathcal{S}_2 \}, \label{eq:defIntersection} \\
\begin{split} \label{eq:defConvHull}
& conv(\mathcal{S}_{1},\mathcal{S}_{2}) := \bigg \{\sum_{i=1}^{n+1} \lambda_i\,s_i~\bigg|~s_i \in \mathcal{S}_1 \cup \mathcal{S}_2,~\lambda_i \geq 0,~\sum_{i=1}^{n+1} \lambda_i = 1 \bigg\}.
\end{split}
\end{align}

In this paper, we represent reachable sets by zonotopes:

(Zonotope) Given a center vector $c \in \Rn{}$ and a generator matrix $G \in \R^{n \times \gens{}}$, a zonotope $\mathcal{Z} \subset \Rn{}$ is

\begin{equation*}
\mathcal{Z} := \bigg\{ c + \sum_{i = 1}^{\gens} G_{(\cdot,i)} \, \alpha_i ~\bigg| ~ \alpha_i \in [-1,1] \bigg\},
\end{equation*}

where the scalars $\alpha_i$ are called factors. We use the shorthand $\mathcal{Z} = \langle c,G \rangle_Z$.

For zonotopes $\mathcal{Z}_{1} = \langle c_1,G_1 \rangle_Z,\mathcal{Z}_{2} = \langle c_2,G_2 \rangle_Z \subset \Rn{}$ with $\gamma_1$ and $\gamma_2$ generators, respectively, linear map, Minkowski sum, and convex hull can be computed as <cit.>

\begin{align}
& M \mathcal{Z}_{1} = \langle M c_1, M G_1 \rangle_Z , \label{eq:zonoLinTrans} \\
& \mathcal{Z}_{1} \oplus \mathcal{Z}_{2} = \langle c_1 + c_2, [G_1~G_2] \rangle_Z , \label{eq:zonoMinSum} \\
& conv(\mathcal{Z}_{1},\mathcal{Z}_{2}) \subseteq \big \langle 0.5(c_1 + c_2),~\big[ 0.5 (G_1 + G^{(1)}_2) \nonumber \\
& \qquad \qquad \qquad \qquad  0.5 (G_1 - G_2^{(1)})~ 0.5(c_1 - c_2)~ G_2^{(2)}  \big] \big \rangle_Z \label{eq:zonoConvHull}
\end{align}

\begin{equation*}
G_2^{(1)} = [G_{2(\cdot,1)} \, \dots \, G_{2(\cdot,\gens{}_1)}], ~ G_2^{(2)} = [G_{2(\cdot,\gens{}_1+1)} \, \dots \, G_{2(\cdot,\gens{}_2)}],
\end{equation*}

where we assume without loss of generality that $\mathcal{Z}_{2}$ has more generators than $\mathcal{Z}_{1}$. In addition, the multiplication $\mathbfcal{I} \, \mathcal{Z}$ of an interval matrix $\mathbfcal{I}$ with a zonotope $\mathcal{Z}$ can be enclosed as specified in <cit.>, and the operation $\texttt{interval}(\mathcal{Z})$ returns
the interval enclosure according to <cit.>.

Another set representation we require are polytopes, for which we consider the halfspace representation:

(Polytope) Given a constraint matrix $C \in \R^{\poly \times n}$ and a constraint offset $d \in \R^{\poly}$, a polytope $\mathcal{P} \subseteq \R^{n}$ is defined as

\begin{equation*}
\mathcal{P} := \big \{ x \in \Rn ~\big|~ C \,x \leq d \big \}.
\end{equation*}

We use the shorthand $\mathcal{P} = \langle C,d \rangle_P$.

The intersection of two polytopes $\mathcal{P}_1 = \langle C_1,d_1\rangle_P,\mathcal{P}_2 = \langle C_2,d_2 \rangle_P \subseteq \Rn$ can be computed as

\begin{equation}
\mathcal{P}_1 \cap \mathcal{P}_2 = \big \langle [C_1 \concat C_2],[d_1 \concat d_2] \big \rangle_P.
\end{equation}

Computation of the intersection might result in redundant halfspaces, which can be removed using linear programming.

§.§ Reachability Analysis

The reachable set of a dynamic system is defined as follows:

(Reachable Set) We consider a dynamic system whose behavior is described by the differential equation

\begin{equation}
\dot x(t) = f\big(x(t),u(t)\big),
\label{eq:ode}
\end{equation}

where $x(t) \in \Rn$ is the system state and $u(t) \in \R^m$ is the input. Given an initial set $\initSet \subset \Rn$ and an input set $\mathcal{U} \subset \R^m$, the reachable set at time $t \geq 0$ is defined as

\begin{equation*}
\mathcal{R}^{\text{\normalfont e}}(t) := \big \{ \xi(t,x_0,u(\cdot))~\big |~ x_0 \in \initSet, \forall \theta \in [0,t]:~ u(\theta) \in \mathcal{U} \big \},
\end{equation*}

where $\xi(t,x_0,u(\cdot))$ denotes the solution to (<ref>) for the initial state $x_0 = x(0)$ and the input signal $u(\cdot)$.

Since the exact reachable set $\mathcal{R}^{\text{e}}(t)$ as defined in Def. <ref> cannot be computed in general, the goal of reachability analysis is to calculate a tight enclosure $\mathcal{R}(t) \supseteq \mathcal{R}^{\text{\normalfont e}}(t)$ instead. Moreover, it is common practice to compute the reachable set for consecutive time intervals $\tau_i = [t_i,t_{i+1}]$ with $t_i = i \cdot \Delta t$, $i \in \{0,\dots,t_{\text{end}}/\Delta t\}$, where $\Delta t$ is the time step size and $t_{\text{end}}$ is the final time. We assume without loss of generality that time starts at $t_0 = 0$ and $t_{\text{end}}$ is a multiple of $\Delta t$. Consequently, the reachable set for the whole time horizon $\mathcal{R}([0,t_{\text{end}}])$ is given by a sequence

\begin{equation}
\mathcal{R}(t_0),\mathcal{R}(\tau_0),\mathcal{R}(t_1),\mathcal{R}(\tau_1),\dots,\mathcal{R}(\tau_{t_{\text{end}}/\Delta t -1}),\mathcal{R}(t_{\text{end}})
\label{eq:reachSeq}
\end{equation}

of time point reachable sets $\mathcal{R}(t_i)$ and time interval reachable sets $\mathcal{R}(\tau_i)$. For the remainder of the paper, we use the shorthand notation $\forall u(\cdot) \in \mathcal{U}$ to refer to the set $\{u(\cdot)\,|\, \forall \theta \in [0,t_{\text{end}}]:\,u(\theta) \in \mathcal{U} \}$ of all input signals that are contained in the input set $\mathcal{U}$ at all times.

Some reachability algorithms preserve dependencies between initial states and inputs and the corresponding reachable states <cit.>. For those algorithms, the computed enclosure of the reachable set can be used to construct a symbolic function that approximates the solution to the differential equation within a guaranteed error bound:
(Dependency Preservation)
Let $\mathcal{R}(t)$ be an enclosure of the reachable set for the initial set $\mathcal{X}_0$ and the set of uncertain inputs $\mathcal{U}$ computed with a reachability algorithm $\mathcal{A}$. The reachability algorithm $\mathcal{A}$ is dependency preserving if the computed reachable set $\mathcal{R}(t)$ can be used to construct a function $\mu(t,x_0,u(\cdot))$ that approximates the solution $\xi(t,x_0,u(\cdot))$ of the differential equation in (<ref>) within some time-varying error bound $\mathcal{E}(t)$ for all initial states $x_0 \in \mathcal{X}_0$ and input signals $u(\cdot) \in \mathcal{U}$:

\begin{equation*}
\forall t \in [0,t_{\text{\normalfont end}}], \forall x_0 \in \mathcal{X}_0, \forall u(\cdot) \in \mathcal{U}: ~\xi(t,x_0,u(\cdot)) \in \mu(t,x_0,u(\cdot)) \oplus \mathcal{E}(t) \, \subseteq \, \mathcal{R}(t),
\end{equation*}

where the trivial solution $\mu(t,x_0,u(\cdot)) = \mathbf{0}$ and $\mathcal{E}(t) = \mathcal{R}(t)$ is only allowed if the reachable set enclosure $\mathcal{R}(t)$ is singleton.

The main advantage of dependency preservation is that with the analytical relation $\mu(t,x_0,u(\cdot))$ between initial states and inputs and the corresponding reachable states enclosures of reachable sets for single initial states and inputs can be computed very efficiently in terms of a simple function evaluation.

§.§ Temporal Logic

We consider specifications in signal temporal logic <cit.>:

(Signal Temporal Logic) The syntax of a signal temporal logic formula over a finite set of atomic predicates $\pred \in \mathcal{AP}$ is

\begin{equation*}
\varphi := \pred \,|\, \neg \varphi \,|\, \varphi_1 \vee \varphi_2 \,|\, \varphi_1 \wedge \varphi_2 \,|\, \varphi_1 \, \until_{[a,b]} \, \varphi_2 \,|\, \varphi_1 \, \release_{[a,b]} \, \varphi_2 \,|\, \finally_{[a,b]} \, \varphi  \,|\, \globally_{[a,b]} \, \varphi \,|\, \nextOp_a \, \varphi
\end{equation*}

with $a,b \in \R_{\geq 0}$ and $b \geq a$, where we in addition to the until operator $\varphi_1 \, \until_{[a,b]} \, \varphi_2$ consider the operators release $\varphi_1 \, \release_{[a,b]} \, \varphi_2 := \neg( \neg \varphi_1 \, \until_{[a,b]} \, \neg \varphi_2)$, finally $\finally_{[a,b]} \, \varphi := \mathrm{true}~ \until_{[a,b]} \, \varphi$, globally $\globally_{[a,b]} \, \varphi := \neg \finally_{[a,b]} \, \neg \varphi$, and next $\nextOp_a \, \varphi := \finally_{[a,a]} \, \varphi$. For a trace $\xi(t)$, the semantics of a signal temporal logic formula is defined as follows:

\begin{equation*}
\begin{split}
\xi \vDash \pred & ~~ \Leftrightarrow ~~ \pi_{\pred}(\xi(0)) = \mathrm{true} \\
\xi \vDash \neg \varphi & ~~ \Leftrightarrow ~~ \neg (\xi \vDash \varphi) \\
\xi \vDash \varphi_1 \vee \varphi_2 & ~~ \Leftrightarrow ~~ (\xi \vDash \varphi_1) \vee (\xi \vDash \varphi_2) \\
\xi \vDash \varphi_1 \, \until_{[a,b]} \, \varphi_2 & ~~\Leftrightarrow ~~ \exists t \in [a,b]: \xi_t \vDash \varphi_2 ~\wedge~ \forall t' \in [0,t): \xi_{t'} \vDash \varphi_1
\end{split}
\end{equation*}

using a predicate evaluation function $\pi_{\pred}$ and notation $\xi_a(t) := \xi(t+a)$. The semantics for the conjunction follows from the equality $\varphi_1 \wedge \varphi_2 = \neg(\neg \varphi_1 \vee \neg \varphi_2)$.

In this work we consider atomic predicates given by linear inequality constraints such as $2 \, x_1 - 3 \, x_2 \leq 1$.

Please note that this also includes polytope containment since $x \in \mathcal{P}$ with $\mathcal{P} = \langle C,d \rangle_P$ can be represented as $C_{(1,\cdot)} x \leq d_{(1)} \wedge \dots \wedge C_{(\poly,\cdot)} x \leq d_{(\poly)}$. Reachset temporal logic <cit.> is a special type of logic that can be directly evaluated on a reach sequence (<ref>):

(Reachset Temporal Logic) The syntax of a reachset temporal logic formula[In contrast to <cit.>, we use a shortened definition of RTL that is restricted to the relevant parts.] over a finite set of atomic predicates $\pred \in \mathcal{AP}$ is

\begin{equation*}
\varphi := \all \,\pred ~|~ \varphi_1 \vee \varphi_2 ~|~ \varphi_1 \wedge \varphi_2 ~|~  \nextOp_a \, \varphi
\end{equation*}

with $a = i \cdot \Delta t/2$ and $i \in \mathbb{N}_{\geq 0}$, where $\Delta t$ is the time step size. The all operator $\all \, \pred$ specifies that a predicate is satisfied for all states inside the reachable set, and the next operator $\nextOp_a \, \varphi$ refers to a time point reachable set if $a$ is an integer multiple of $\Delta t$ and to a time interval reachable set otherwise. For a reachable set $\mathcal{R}(t)$, the semantics of a reachset temporal logic formula is defined as follows:

\begin{equation*}
\begin{split}
\mathcal{R} \vDash \all \, \pred & ~~ \Leftrightarrow ~~ \forall r \in \mathcal{R}(0): ~\pi_{\pred}(r) = \mathrm{true} \\
\mathcal{R} \vDash \varphi_1 \vee \varphi_2 & ~~ \Leftrightarrow ~~ (\mathcal{R} \vDash \varphi_1) \vee (\mathcal{R} \vDash \varphi_2) \\
\mathcal{R} \vDash \varphi_1 \wedge \varphi_2 & ~~ \Leftrightarrow ~~ (\mathcal{R} \vDash \varphi_1) \wedge (\mathcal{R} \vDash \varphi_2) \\
\mathcal{R} \vDash \nextOp_a \, \varphi & ~~\Leftrightarrow ~~ \begin{cases}\mathcal{R}_a \vDash \varphi, & a ~ \text{\normalfont mod} ~ \Delta t = 0\\ \forall t \in \left \lfloor \frac{a}{\Delta t} \right \rfloor + [0, \Delta t]: ~ \mathcal{R}_t \vDash \varphi, & \text{\normalfont otherwise}\end{cases}
\end{split}
\end{equation*}

using a predicate evaluation function $\pi_{\pred}$ and notation $\mathcal{R}_a(t) := \mathcal{R}(t+a)$.

The paper <cit.> that introduced RTL also provides an approach for converting a temporal logic formula in STL to RTL. Given an STL formula $\varphi$, the resulting RTL formula $\varphi_{\text{rtl}}$ has according to <cit.> the following structure:

\begin{equation*}
\varphi_{\text{rtl}} = \bigwedge_{h = 0}^H \bigvee_{j = 0}^J \nextOp_{j \, \Delta t/2} \bigvee_{k=0}^K \all \, \pred_{hjk},
\end{equation*}

where $H$, $J$, and $K$ denote the number of conjunctions and disjunctions.

Since we consider the case where the atomic predicates $\pred_{hjk}$ are all given by linear inequality constraints, the entailment check $\mathcal{R} \vDash \varphi_{\text{rtl}}$ can be equivalently formulated in terms of intersection checks between the reachable set $\mathcal{R}(t)$ and polytopes $\mathcal{P}_{hjkl}$ in halfspace representation <cit.>:

\begin{equation}
\mathcal{R} \vDash \varphi_{\text{rtl}} ~~ \Leftrightarrow ~~\bigwedge_{h=0}^H \bigvee_{j=0}^J \bigvee_{k=0}^K \bigg( \widetilde{\mathcal{R}}(j \, \Delta t/2) \cap \bigcup_{l=0}^L \mathcal{P}_{hjkl} = \emptyset \bigg),
\label{eq:RTLfinal}
\end{equation}

where $L$ is the number of polytopes required to represent the atomic predicates $\pred_{hjk}$ and the auxiliary variable $\widetilde{\mathcal{R}}(t)$ is used to distinguish between time point and time interval reachable sets:

\begin{equation*}
\widetilde{\mathcal{R}}(t) = \begin{cases} \mathcal{R}(t_j), & t ~ \text{\normalfont mod} ~ \Delta t = 0  \\ \mathcal{R}(\tau_j), & \text{\normalfont otherwise} \end{cases} ~~~~\text{with} ~~~~ j = \lfloor t / \Delta t \rfloor.
\end{equation*}

Since we use the same time step size $\Delta t$ for reachability analysis and for the time-discretization of reachset temporal logic, (<ref>) can be directly evaluated on a reach sequence (<ref>).
An example demonstrating the conversion from STL to RTL is provided in <ref>.

According to <cit.>, it holds that if the reachable set satisfies the RTL formula $\varphi_{\text{rtl}}$, then it is guaranteed that all traces contained in the reachable set satisfy the corresponding STL formula $\varphi$. Due to the conservatism introduced by the time-discretization for RTL with $\Delta t$, the contrary does not hold, meaning that satisfaction of the STL formula does not automatically imply satisfaction of the RTL formula. However, since the conservatism introduced by the time-discretization converges to zero for $\Delta t \to 0$ <cit.>, satisfaction of STL and RTL becomes equivalent in the limit case $\Delta t \to 0$.

§ PROBLEM FORMULATION

We consider linear time-invariant systems

\begin{align}
\dot{x}(t)	&= A \, x(t) + B \, u(t) \label{eq:linsys}
\end{align}

with $A \in \R^{n \times n}$, $B \in \R^{n \times m}$, where $x(t) \in \R^{n}$ is the state and $u(t) \in \R^{m}$ is the input.
The initial state $x(t_0)$ is uncertain within the initial set $\initSet \subset \R^{n}$ and the input $u(t)$ is uncertain within the input set $\mathcal{U} \subset \R^{m}$. In this work, we assume that $\initSet$ and $\mathcal{U}$ are represented by zonotopes.

While we consider the set of uncertain inputs $\mathcal{U}$ to be constant over time for simplicity, the extension to a time-varying set of uncertain inputs $\mathcal{U}(t)$ is straightforward. Time-varying sets of uncertain inputs can for example be used to pass external signals such as reference trajectories or a sequence of control commands to the model.

Given an STL formula $\varphi$, our goal is to decide whether the system (<ref>) satisfies $\varphi$ for all initial states $x_0 \in \initSet$ and all input signals $u(\cdot) \in \mathcal{U}$:

\begin{equation*}
\forall t, \forall x_0 \in \initSet, \forall u(\cdot) \in \mathcal{U}: ~ \xi(t,x_0,u(\cdot)) \vDash \varphi,
\end{equation*}

where $\xi(t,x_0,u(\cdot))$ is the solution to (<ref>) for the initial state $x_0 = x(t_0)$ and the input signal $u(\cdot)$.

§ REACHABILITY ANALYSIS

For the approach presented in this paper we require a reachability algorithm that is dependency preserving according to Def. <ref>. To construct such an algorithm, we modify the reachability algorithm in <cit.> by using a different enclosure for the reachable set due to uncertain inputs $\mathcal{P}(\Delta t)$. In particular, to preserve dependencies on the inputs, we approximate the reachable set due to time-varying inputs with the reachable set due to constant inputs $\mathcal{P}_c(\Delta t) = T \, \mathcal{U}_0$ in each time step, where

\begin{equation}
T = A^{-1} (e^{A \Delta t} - I_n)
\label{eq:propMat}
\end{equation}

is the propagation matrix for constant inputs and $\mathcal{U}_0 = B \big (\mathcal{U} \oplus (-c_u)\big) \subset \R^{n}$ is an auxiliary variable defined using the geometric center $c_u$ of the input set $\mathcal{U} = \langle c_u,G_u\rangle_Z$. If the matrix $A$ is not invertible, $A^{-1}$ in (<ref>) can be integrated into the power-series for the exponential matrix <cit.>. Since $\mathcal{P}_c(\Delta t) = T \, \mathcal{U}_0$ is computed using a linear map, it preserves dependencies according to <cit.>. We account for the approximation error with the bloating term

\begin{equation}
\mathcal{D} = \bigg( \sum_{j=1}^\kappa \frac{A^j \Delta t^{j+1}}{(j+1)!} \bigg) \mathcal{U}_0 \oplus \bigoplus_{j=1}^\kappa \frac{A^j \Delta t^{j+1}}{(j+1)!}\, \mathcal{U}_0 \oplus 2 \,  \Delta t \, \mathbfcal{E} \, \mathcal{U}_0,
\label{eq:inputDiff}
\end{equation}

which according to <cit.> encloses the difference between constant and time-varying inputs. Consequently, we overall obtain the following enclosure for the reachable set due to time-varying inputs:

\begin{equation}
\mathcal{P}(\Delta t) \subseteq T \, \mathcal{U}_0 \oplus \mathcal{D}.
\label{eq:inputReach}
\end{equation}

Another essential part of the reachability algorithm in  <cit.> is the set $\mathcal{C}_i$ that accounts for the curvature of trajectories and is according to <cit.> given as

\begin{equation} \label{eq:Cu}
\mathcal{C}_i = (e^{A \Delta t})^{i} \big( \mathbfcal{F} \, \mathcal{X}_0 \oplus \mathbfcal{G} \, \widetilde{u} \big)
\end{equation}

using the auxiliary variable $\widetilde u = B c_u \in \R^{n}$ and the interval matrices <cit.>

\begin{equation}
\mathbfcal{F} = \mathbfcal{T}^{(\kappa)} A \oplus \mathbfcal{E}, ~ \mathbfcal{G} = \mathbfcal{T}^{(\kappa + 1)} \oplus \mathbfcal{E} \Delta t, ~ \mathbfcal{T}^{(o)} = \bigoplus_{j=2}^o \big[\big( j^{\frac{-j}{j-1}} - j^{\frac{-1}{j-1}} \big) \Delta t^j,0\big] \frac{A^{j-1}}{j!},
\label{eq:curv}
\end{equation}

where the parameter $\kappa$ and the interval matrix $\mathbfcal{E}$ are the truncation order and the remainder of the Taylor series for the exponential matrix <cit.>:

\begin{equation*} %\label{eq:E}
\mathbfcal{E} = [-E,E] , \quad
E = e^{|A|\Delta t} - \sum_{j=0}^{\kappa} \frac{\big( |A|\Delta t \big)^j}{j!} .
\end{equation*}

Using the enclosure of the reachable set due to time-varying inputs $\mathcal{P}(\Delta t)$ in (<ref>) as well as the curvature enclosure $\mathcal{C}_i$ in (<ref>), the reachability algorithm in <cit.> finally computes tight enclosures of the time point reachable set $\mathcal{R}(t_i)$ and the time interval reachable set $\mathcal{R}(\tau_i)$ as follows:

\begin{equation}
\begin{split}
& \mathcal{R}(t_{i}) = \mathcal{H}(t_{i}) \oplus \mathcal{P}(t_{i}), \\
& \mathcal{R}(\tau_{i}) = conv\big(\mathcal{H}(t_i), \mathcal{H}(t_{i+1})\big) \oplus \mathcal{P}(t_{i+1}) \oplus \mathcal{C}_i,
\end{split}
\label{eq:reach}
\end{equation}

where the homogeneous solution $\mathcal{H}(t_i)$ resulting from the propagation of the initial set and constant inputs as well as the particular solution $\mathcal{P}(t_i)$ due to uncertain time-varying inputs are computed using the following propagation scheme:

\begin{equation}
\begin{split}
\mathcal{H}(t_{i+1}) &= e^{A \Delta t} \, \mathcal{H}(t_{i}) \oplus T \, \widetilde{u}, \\
\mathcal{P}(t_{i+1}) &= e^{A \Delta t} \, \mathcal{P}(t_i) \oplus \mathcal{P}(\Delta t),
\end{split}
\label{eq:propScheme}
\end{equation}

where the initial values are $\mathcal{H}(0) = \mathcal{X}_0$ and $\mathcal{P}(0) = \mathbf{0}$. The resulting reachability algorithm is summarized in Alg. <ref>, where we enclose $\mathcal{C}_i$ and $\mathcal{D}$ by intervals to reduce the number of generators of the reachable set.

Reachability Analysis

Input: Linear system $\dot x = A x + B u$, initial set $\initSet = \langle c_x,G_x\rangle_Z$, input set $\mathcal{U}= \langle c_u,G_u \rangle_Z$, final time $t_{\text{end}}$, time step size $\Delta t$, truncation order $\kappa$.

Output: Reach sequence $\mathcal{R}(t_0),\mathcal{R}(\tau_0),\mathcal{R}(t_1),\mathcal{R}(\tau_1),\dots,\mathcal{R}(\tau_{t_{\text{end}}/\Delta t-1}),\mathcal{R}(t_{\text{end}})$ with time point reachable sets $\mathcal{R}(t_i)$ and time interval reachable sets $\mathcal{R}(\tau_i)$.

[1]
$\widetilde{u} \gets B \, c_u$, $\mathcal{U}_0 \gets B \big( \mathcal{U} \oplus (-c_u) \big) $, $T \gets \eqref{eq:propMat}$, $\mathbfcal{F},\mathbfcal{G} \gets \eqref{eq:curv}$, $\mathcal{D} \gets \eqref{eq:inputDiff}$
$t_0 \gets 0$, $\mathcal{H}(t_0),\mathcal{R}(t_0) \gets \initSet$, $\mathcal{P}(t_0),\mathcal{P}_c(t_0),\mathcal{D}_0 \gets \mathbf{0}$
$\mathcal{C}_0 \gets \texttt{interval}(\mathbfcal{F} \, \mathcal{X}_0 \oplus \mathbfcal{G} \, \widetilde{u})$
$i \gets 0$ to $t_{\text{end}}/\Delta t -1$
$t_{i+1} \gets t_i + \Delta t$, $\tau_i \gets [t_i,t_{i+1}]$, $\mathcal{C}_{i+1} \gets e^{A \Delta t} \mathcal{C}_i$
$\mathcal{H}(t_{i+1}) \gets e^{A \Delta t} \, \mathcal{H}(t_{i}) \oplus T \, \widetilde{u}$, $\mathcal{D}_{i+1} \gets \texttt{interval}(e^{A \Delta t} \, \mathcal{D}_i \oplus \mathcal{D})$
$\mathcal{P}_c(t_{i+1}) \gets e^{A \Delta t} \, \mathcal{P}_c(t_i) \oplus T \, \mathcal{U}_0$, $\mathcal{P}(t_{i+1}) \gets \mathcal{P}_c(t_{i+1}) \oplus \mathcal{D}_{i+1}$
$\mathcal{R}(t_{i+1}) \gets \mathcal{H}(t_{i+1}) \oplus \mathcal{P}(t_{i+1})$
$\mathcal{R}(\tau_{i}) \gets conv\big(\mathcal{H}(t_i), \mathcal{H}(t_{i+1})\big) \oplus \mathcal{P}(t_{i+1}) \oplus \mathcal{C}_i$

An important property of Alg. <ref> is that it is dependency preserving for piecewise constant inputs. Please note that it is sufficient to consider piecewise constant inputs in our case since those approximate time-varying inputs arbitrary well for $\Delta t \to 0$. We therefore now show how to construct the approximate solution $\mu(t,x_0,u(\cdot))$ to the differential equation as well as the error bound $\mathcal{E}(t)$ in Def. <ref> from the computed reachable set:

We consider the initial set $\initSet = \langle c_x,G_x \rangle_Z \subset \Rn$ with $\gens_x$ generators, the input set $\mathcal{U} = \langle c_u,G_u \rangle_Z \subset \R^m$ with $\gens_u$ generators, and the corresponding enclosure of the reachable set $\mathcal{R}(t)$ computed with Alg. <ref> using time step size $\Delta t$. Given an initial state $x_0 \in \initSet$ as well as the piecewise constant input signal $u(\cdot) \in \mathcal{U}$ defined as

\begin{equation}
x_0 = c_x + G_x \alpha_x,\quad \forall t \in [(j-1) \Delta t,j \Delta t]:~ u(t) = c_u + G_u \alpha_{u,j}
\label{eq:alpha}
\end{equation}

with $j = 1,\dots,t_{\text{\normalfont end}}/\Delta t$, the approximate solution $\mu(t,x_0,u(\cdot))$ and the time-varying error bound $\mathcal{E}(t)$ required for Alg. <ref> to be dependency preserving according to Def. <ref> are given as

\begin{equation}
\begin{split}
\mu(t,x_0,u(\cdot)) &= \begin{cases} [G_{i1(\cdot,\mathbf{H}_{i})}~\mathbf{0}] \, \widetilde{\alpha}, & t~\text{\normalfont mod}~\Delta t = 0 \\ [G_{i2(\cdot,\mathbf{N}_{i})}~\mathbf{0}] \, \widetilde{\alpha}, & \text{\normalfont otherwise} \end{cases} \\
\mathcal{E}(t) &= \begin{cases} \langle c_{i1},G_{i1(\cdot,\mathbf{K}_{i})}\rangle_Z, & t~\text{\normalfont mod}~\Delta t = 0 \\ \langle c_{i2},G_{i2(\cdot,\mathbf{M}_{i})}\rangle_Z, & \text{\normalfont otherwise}\end{cases}
\end{split}
\label{eq:dependency}
\end{equation}

\begin{equation*}
i = \lfloor t/\Delta t  \rfloor, ~~~~\mathcal{R}(t_i) = \langle c_{i1},G_{i1} \rangle_Z, ~~~~\mathcal{R}(\tau_i) = \langle c_{i2},G_{i2} \rangle_Z,
\end{equation*}

where the vector

\begin{equation}
\widetilde{\alpha} = [\alpha_x \concat \alpha_{u,1} \concat \dots \concat \alpha_{u,{t_{\text{\normalfont end}}/\Delta t}}]
\label{eq:alphaVec}
\end{equation}

stores the factor values that define the initial state and input signal in (<ref>) and the tuples

\begin{equation}
\begin{split}
& \mathbf{H}_{i} = (1,\dots,\gamma_x,\gamma_x + 1,\dots,\gamma_x + \gamma_u \, i ), ~ \mathbf{K}_{i} = (1,\dots,\gamma_{i1}) \setminus \mathbf{H}_{i} \\
\mathbf{N}_{i} = &( 1, \dots,\gamma_x,2 \, \gamma_x + 2,\dots,2 \, \gamma_x + 2 + \gamma_u \, (i+1) ), ~ \mathbf{M}_{i} = (1,\dots,\gamma_{i2}) \setminus \mathbf{N}_{i}  \\
\end{split}
\label{eq:indices}
\end{equation}

store the indices of specific zonotope generators.

Instead of directly mapping initial states and inputs to reachable states, our approximate solution $\mu(t,x_0,u(\cdot))$ relies on an intermediate representation by zonotope factors $\widetilde{\alpha}$. We therefore first represent the initial state $x_0$ as well as the input signal $u(\cdot)$ by the corresponding zonotope factors $\widetilde{\alpha}$ using the zonotopes $\mathcal{X}_0$ and $\mathcal{U}$, and then map these zonotope factors to reachable states using the zonotopes that represent the reachable set $\mathcal{R}(t)$. In particular, according to Def. <ref>, every point $x \in \mathcal{Z}$ inside a zonotope $\mathcal{Z} = \langle c,G \rangle_Z \subset \Rn$ can equivalently be represented by the corresponding zonotope factors $\alpha \in \R^\gamma$ using the relation $x = c + G\, \alpha$. This allows us to equivalently represent the initial state $x_0 \in \mathcal{X}_0$ and the piecewise constant input signal $u(\cdot) \in \mathcal{U}$ by the vector of zonotope factors $\widetilde{\alpha}$ in (<ref>) according to (<ref>).

Since Alg. <ref> is composed of the set operations linear map, Minkowski sum, and convex hull, which are all dependency preserving for zonotopes according to <cit.>, it holds that Alg. <ref> is dependency preserving, too <cit.>. This means that we can directly insert $\widetilde{\alpha}$ into the zonotopes that represent the reachable set $\mathcal{R}(t)$ in (<ref>) to map from factors to reachable states, where we use the time point reachable set $\mathcal{R}(t_i)$ if the time $t$ is equivalent to $t_i$ and the time interval reachable set $\mathcal{R}(\tau_i)$ if the time $t$ is between two time points. When mapping from factors to reachable states we have to distinguish between generators that correspond to initial states and inputs, and generators that represent uncertainty arising from the difference between constant and time-varying inputs $\mathcal{D}$, the curvature enclosure $\mathcal{C}_i$, or the over-approximation from the enclosure of the convex hull of two zonotopes computed using (<ref>). We achieve this with the tuples $\mathbf{H}_i,\mathbf{N}_i,\mathbf{K}_i,\mathbf{M}_i$ in (<ref>), where $\mathbf{H}_i,\mathbf{N}_i$ store the indices of the generators that correspond to initial states as well as inputs and the tuples $\mathbf{K}_i,\mathbf{M}_i$ store the indices of the generators that represent uncertainty.

In addition to dependency preservation, another important property of Alg. <ref> is that the computed enclosure of the reachable set converges to the exact reachable set for $\Delta t \to 0$. The corresponding proof is provided in <ref>. As a direct consequence, also the size of the error $\mathcal{E}(t)$ in Prop. <ref> converges to 0 for $\Delta t \to 0$, so that the approximate solution $\mu(t,x_0,u(\cdot))$ converges to the exact solution $\xi(t,x_0,u(\cdot))$. While a large over-approximation error for reachability analysis might prevent the successful verification of the temporal logic specification, improving the accuracy to achieve a smaller error increases the computation time. To resolve this trade-off, our strategy for automated verification is to start with a quite inaccurate reachable set enclosure that can be computed very efficiently, and then iteratively refine the accuracy until the temporal logic specification can be either verified or falsified.

§ MODEL CHECKING

The main disadvantage of model checking using RTL is that the intersection checks in (<ref>) are evaluated in a pure true/false manner, losing all information about which initial states or inputs result in an intersection. As visualized in Fig. <ref>, the model checking procedure therefore also considers so-called spurious traces, which are contained in the reachable set but are not consistent with the system dynamics. In this work we solve this problem by explicitly keeping track which initial states and which inputs result in an intersection.

In particular, we approximate the set of all possible time-varying inputs $u(\cdot) \in \mathcal{U}$ with the set of all possible piecewise constant inputs, because our reachability algorithm Alg. <ref> is dependency preserving for piecewise constant inputs. This approximation does not negatively affect the completeness of our model checking algorithm, since in the limit case $\Delta t \to 0$ piecewise constant inputs approximate time-varying inputs arbitrary well. For our model checking algorithm, we have to consider the set containing all possible combinations of initial states and piecewise constant inputs, which is given as $\mathcal{X}_0 \times \mathcal{U} \times \dots \times \mathcal{U}$. As shown in (<ref>), since both the initial set and the input set are represented by zonotopes, we can equivalently represent every initial state and every piecewise constant input signal by the corresponding zonotope factors $\alpha_x$ and $\alpha_{u,i}$. Therefore, the set $\mathcal{X}_0 \times \mathcal{U} \times \dots \times \mathcal{U}$ containing all possible combinations of initial states and inputs can equivalently be represented in the space of zonotope factors by the hypercube $\widetilde \alpha = [-\mathbf{1},\mathbf{1}]$ with $\widetilde \alpha$ defined as in (<ref>). Since the hypercube $\widetilde \alpha = [-\mathbf{1},\mathbf{1}]$ is computationally much easier to handle than the zonotope $\mathcal{X}_0 \times \mathcal{U} \times \dots \times \mathcal{U}$, our model checking algorithm represents initial states and inputs in factor space $\widetilde \alpha = [-\mathbf{1},\mathbf{1}]$.

As visualized in Fig. <ref>, for each intersection between a reachable set $\widetilde{\mathcal{R}}(j \, \Delta t/2)$ and a polytope $\mathcal{P}_{hjkl}$ we can divide the domain $\widetilde{\alpha} \in [-\mathbf{1},\mathbf{1}]$ into factors that potentially lead to an intersection and factors that are guaranteed to not result in an intersection with the polytope $\mathcal{P}_{hjkl}$. The corresponding polytopes $\fac$ representing factors that potentially result in a violation can be computed as follows:

Schematic visualization of the model checking approach presented in Sec. <ref>, where the reachable sets $\widetilde{\mathcal{R}}(j \, \Delta t/2)$ are visualized in blue and the polytopes $\mathcal{P}_{hjkl}$ as well as the corresponding sets of unsafe factors $\fac$ are depicted in red. Moreover, the green arrows represent the mapping from the state space to the space of zonotope factors $\widetilde{\alpha} \in [-\mathbf{1},\mathbf{1}]$.

Model Checking

Input: Linear system $\dot x = A x + B u$, STL specification $\varphi$, initial set $\initSet$, input set $\mathcal{U}$, final time $t_{\text{end}}$, time step size $\Delta t$, truncation order $\kappa$.

Output: List $\mathbf{L}$ containing polytopes that represent zonotope factors of the initial and input set which violate the specification.

[1]
$\mathcal{R}(t_0),\mathcal{R}(\tau_0),\dots,\mathcal{R}(t_{\text{end}}) \gets $ comp. reach. set with $\Delta t$ and $\kappa$ using Alg. <ref>
$\bigwedge_h \bigvee_j \bigvee_k (\widetilde{\mathcal{R}}(j \,\Delta t/2) \cap \bigcup_l \mathcal{P}_{hjkl} = \emptyset) \gets$ convert $\varphi$ to RTL using $\Delta t$
[12]  (see (<ref>) and <cit.>)
$\mathbf{L} \gets \emptyset$
$h \gets 0$ to $H$
$\mathbf{V} \gets [-\mathbf{1},\mathbf{1}]$
$j \gets 0$ to $J$ and $k \gets 0$ to $K$
$\mathbf{W} \gets \emptyset$
$l \gets 0$ to $L$
$\widetilde{\mathcal{R}}(j \,\Delta t/2) \cap \mathcal{P}_{hjkl} \neq \emptyset$
$\fac \gets $ unsafe factors from $ \widetilde{\mathcal{R}}(j \,\Delta t/2) \cap \mathcal{P}_{hjkl}$ using Prop. <ref>
$\mathbf{W} \gets (\mathbf{W},\fac)$
$\mathbf{V} \gets (\mathbf{V}_{(1)} \cap \mathbf{W}_{(1)}, \dots, \mathbf{V}_{(|\mathbf{V}|)} \cap \mathbf{W}_{(|\mathbf{W}|)})$
$\mathbf{V} \gets $ remove empty polytopes from $\mathbf{V}$
$\mathbf{L} \gets (\mathbf{L}, \mathbf{V})$

Given the time point reachable set $\mathcal{R}(t_i) = \langle c_{i1},G_{i1}\rangle_Z$ and time interval reachable set $\mathcal{R}(\tau_i) = \langle c_{i2},G_{i2} \rangle_Z \subset \mathbb{R}^n$ computed using Alg. <ref> as well as the polytope $\mathcal{P} = \langle C,d \rangle_P \subseteq \mathbb{R}^n$, it holds that for all zonotope factors $\widetilde \alpha$ that are not located in the polytopes

\begin{align*}
\fac_{i1} & = \bigg \langle \big[C G_{i1(\cdot,\mathbf{H}_i)}~\mathbf{0} \big],~d - C c_{i1} + \sum_{j \in \mathbf{K}_i} | C G_{i1(\cdot,j)} | \bigg \rangle_P \\
\fac_{i2} & = \bigg \langle \big[C G_{i2(\cdot,\mathbf{N}_i)}~\mathbf{0} \big],~d - C c_{i2} + \sum_{j \in \mathbf{M}_i} | C G_{i2(\cdot,j)} | \bigg \rangle_P
\end{align*}

there is no intersection between the corresponding reachable states and $\mathcal{P}$:

\begin{equation*}
\begin{split}
\forall \widetilde{\alpha} \in [-\mathbf{1},\mathbf{1}]: ~~ &\big( \widetilde{\alpha} \not \in \fac_{i1} \big) \Rightarrow \big( \xi(t_i,x_0,u(\cdot)) \not \in \mathcal{P} \big) \\
\forall \widetilde{\alpha} \in [-\mathbf{1},\mathbf{1}],\, \forall t \in \tau_i: ~~ &\big( \widetilde{\alpha} \not \in \fac_{i2} \big) \Rightarrow \big( \xi(t,x_0,u(\cdot)) \not \in \mathcal{P} \big),
\end{split}
\end{equation*}

where $\widetilde{\alpha} = [\alpha_x \concat \alpha_{u,1} \concat \dots \concat \alpha_{u,t_{\text{\normalfont end}}/\Delta t}]$, $x_0$ and $u(\cdot)$ are defined as in (<ref>), and the tuples $\mathbf{H}_i$,$\mathbf{K}_i$,$\mathbf{N}_i$ and $\mathbf{M}_i$ are defined as in (<ref>).

According to Prop. <ref> and Def. <ref> it holds that

\begin{equation}
\xi(t_i,x_0,u(\cdot)) \in \underbrace{[G_{i1(\cdot,\mathbf{H}_{i})}~\mathbf{0}] \, \widetilde{\alpha}}_{\mu(t,x_0,u(\cdot))}  \oplus \underbrace{\langle c_{i1},G_{i1(\cdot,\mathbf{K}_{i})}\rangle_Z}_{\mathcal{E}(t)}.
\label{eq:proof1}
\end{equation}

Moreover, a reachable state $\xi(t_i,x_0,u(\cdot))$ is according to Def. <ref> contained in the polytope $\mathcal{P} = \langle C,d\rangle_P$ if $C \, \xi(t_i,x_0,u(\cdot)) \leq d$ holds. Combining this with (<ref>) yields the condition

\begin{equation*}
\exists x \in \langle C c_{i1}, C G_{i1(\cdot,\mathbf{K}_i)} \rangle_Z:~ [C G_{i1(\cdot,\mathbf{H}_i)} ~ \mathbf{0}] \, \widetilde \alpha + x \leq d.
\end{equation*}

Finally, we bring the zonotope $\langle C c_{i1}, C G_{i1(\cdot,\mathbf{K}_i)} \rangle_Z$ that represents uncertainty in the computed enclosure of the reachable set to the other side of the inequality using the formula for the interval enclosure of a zonotope <cit.>, which yields the polytope $\fac_{i1}$. In a similar way we obtain the result for the time interval reachable set.

As shown in Fig. <ref>, the conjunctions in (<ref>) correspond to intersections of the sets containing safe factors, which is equivalent to a union of the polytopes $\fac$ representing potentially unsafe factors. Similarly, the disjunctions in (<ref>) correspond to a union of safe factors which is equivalent to an intersection of unsafe factors.
Alg. <ref> summarizes the corresponding model checking procedure. The system satisfies the STL specification if the conjunctions and disjunctions cancel out all unsafe sets $\fac$, which corresponds to an empty list $\mathbf{L} = \emptyset$.

Automated Verification

Input: Linear system $\dot x = A x + B u$, STL specification $\varphi$, initial set $\initSet$, input set $\mathcal{U}$.

Output: Safe ($\mathcal{R}^{\text{\normalfont e}}(t) \vDash \varphi$) or unsafe ($\mathcal{R}^{\text{\normalfont e}}(t) \nvDash \varphi$).

[1]
$t_{\text{end}} \gets $ final time for temporal logic formula $\varphi$
$\Delta t \gets t_{\text{end}}$

$\kappa \gets$ increase $\kappa$ until $1 - ||\mathbfcal{T}^{(\kappa)}||_F / ||\mathbfcal{T}^{(\kappa+1)}||_F \leq 10^{-10}$ (see (<ref>))
$\mathbf{L} \gets$ model check $\varphi$ with parameters $\Delta t$ and $\kappa$ using Alg. <ref>
$\mathbf{L} = \emptyset$
$\textbf{L} \gets$ model check $\neg \varphi$ with parameters $\Delta t$ and $\kappa$ using Alg. <ref>
$\exists \widetilde \alpha \in [-\mathbf{1},\mathbf{1}]: ~ \widetilde \alpha \not \in \fac_j ~ \forall \fac_j \in \mathbf{L}$
$\Delta t \gets \Delta t /2$

§ AUTOMATED VERIFICATION

The overall verifier is summarized in Alg. <ref>. Since all temporal operators are bounded in time, the STL formula $\varphi$ is restricted to a finite time interval, which we determine in Line <ref>. We then first initialize the time step size $\Delta t$ with the overall time horizon $t_{\text{end}}$ (see Line <ref>), and afterward refine it in each iteration of the main loop (see Line <ref>). For tuning the truncation order $\kappa$ we apply the strategy from <cit.>, which is based on the interval matrices $\mathbfcal{T}^{(o)}$ in (<ref>).

Starting from $\kappa = 2$, we increase $\kappa$ until the relative change between the Frobenius norms $|| \mathbfcal{T}^{(\kappa)} ||_F$ and $|| \mathbfcal{T}^{(\kappa+1)} ||_F$ computed according to <cit.> is smaller than $10^{-10}$.
As explained in Sec. <ref>, the system satisfies the STL formula if the list $\mathbf{L}$ containing the sets of unsafe factors is empty, which we check in Line <ref>. On the other hand, the system violates the STL formula $\varphi$ if there exists a single initial state and input signal that satisfy the negated formula $\neg \varphi$. Consequently, to prove $\mathcal{R}(t) \not \vDash \varphi$, we first run Alg. <ref> to check if $\mathcal{R}(t)$ satisfies $\neg \varphi$ in Line <ref>, which yields a list of safe sets $\mathbf{L}$. Then, we need to show that there exists a vector of factors $\widetilde \alpha \in [-\mathbf{1},\mathbf{1}]$ that does not intersect any of the safe sets $\fac_j = \langle C_j,d_j\rangle_P \in \mathbf{L}$ (see Line <ref>), and therefore corresponds to an unsafe initial state and input signal. This can be realized by solving the following mixed-integer linear program:

\begin{equation}
\begin{split}
\min_{\widetilde \alpha \in [-\mathbf{1},\mathbf{1}]} \| \widetilde \alpha \|_1  ~~ \text{subject to} ~~ & \forall j \in \{1,\dots,|\mathbf{L}|\}, ~ \forall k \in \{ 1,\dots,\poly_j\}: \\
& C_{j(\cdot,k)} \, \widehat{\alpha}_{jk} > \lambda_{jk} \, d_{j(k)}, ~ - \mathbf{1} \, \lambda_{jk} \leq \widehat{\alpha}_{jk} \leq \mathbf{1} \, \lambda_{jk}, \\
& \textstyle \lambda_{jk} \in \{0,1\}, ~ \widetilde \alpha = \sum_{k=1}^{s_j} \widehat{\alpha}_{jk}, ~\sum_{k=1}^{s_j} \lambda_{jk} = 1.
\end{split}
\label{eq:mixedInteger}
\end{equation}

If (<ref>) has a feasible solution, the system violates the specification $\varphi$. In addition, the optimal solution $\widetilde \alpha = [\alpha_x \concat \alpha_{u,1} \concat \dots \concat \alpha_{u,t_{\text{end}}/\Delta t}]$ for (<ref>) defines a falsifying trajectory

\begin{equation*}
x(t_{i+1}) = e^{A \Delta t} x(t_i) + T(c_u + G_u \alpha_{u,i}), \quad \quad i = 0,\dots, t_{\text{end}}/\Delta t
\end{equation*}

with initial state $x(t_0) = c_x + G_x \alpha_x$.

The verification problem for linear systems against STL specifications is undecidable in general since it is not possible to compute the exact reachable set <cit.>. However, undecidable verification problems only occur if the reachable set is located exactly on the decision boundary between safe and unsafe, and are therefore quite rare in practice. Consequently, most verification problem are decidable, and Alg. <ref> returns the correct result in finite time in these cases:

If the verification problem defined by the linear system $\dot x = Ax + Bu$, the STL specification $\varphi$, the initial set $\mathcal{X}_0$, and the input set $\mathcal{U}$ represents a decidable problem instance, Alg. <ref> will terminate in finite time.

Our verification approach contains the following three sources of over-approximation errors:

* The error from the conservative enclosure of the reachable set.

* The error from the time-discretization of the STL formula during the conversion to RTL.

* The error through the approximation of time-varying inputs with piecewise constant inputs.

As we show in <ref>, the over-approximation error from reachability analysis converges to 0 for $\Delta t \to 0$. Similarly, the error from the time-discretization during conversion to RTL becomes 0 for $\Delta t \to 0$ <cit.>, and the space of piecewise constant inputs contained inside the input set $\mathcal{U}$ converges to the space of time-varying input signals contained in $\mathcal{U}$ for $\Delta t \to 0$. Consequently, all errors converge to 0 for $\Delta t \to 0$, so that Alg. <ref> will always be able to find a $\Delta t$ small enough to either verify or falsify the specification.

If the verification problem instance is not decidable, Alg. <ref> does not terminate. To accelerate the computation, we additionally implemented two improvements for our basic verification approach presented in this section:

$\bullet$ We already evaluate all predicates prior to the conversion to RTL and substitute predicates that are either satisfied or violated for the whole reachable set by true/false.

$\bullet$ We omit the conversion to conjunctive normal form performed during conversion to RTL for cases where it does not bring a benefit. In particular, this is the case if the conversion does not result in a disjunction of atomic predicates.

Finally, Alg. <ref> naturally splits the domain $\widetilde \alpha \in [-\mathbf{1},\mathbf{1}]$ corresponding to the initial set and space of input signals into safe and unsafe subsets. We can therefore directly use the resulting list of unsafe sets $\mathbf{L}$ to identify initial states and input signals that are safe or to cut away all reachable states that violate the temporal logic specification.

§ NUMERICAL EXAMPLES

We now demonstrate the performance of our automated verifier on several challenging benchmark systems. All computations are carried out in MATLAB on a 2.9GHz quad-core i7 processor with 32GB memory. Moreover, we integrated the implementation of our automated verifier into the open-source reachability toolbox CORA <cit.>, and published a repeatability package that reproduces all results presented in this section on CodeOcean[<https://codeocean.com/capsule/4858011/tree/v1>].

§.§ Comparison with Other Approaches

We first compare our method with the tool SMTmc <cit.>, which implements formal verification using SMT solving <cit.>. For the comparison we consider the system and STL specification in Fig. <ref>. Since SMTmc does not support systems with inputs, we slightly modify the system by keeping the input $u$ constant over time, so that it can be modeled as an additional state with dynamics $\dot u = 0$. Even for this very simple system SMTmc needs 1462 seconds to verify that the system is safe, where we used a time step size of $\Delta t = 0.1$. Our approach, on the other hand, verifies safety in only 0.25 seconds and is in addition fully automatic, so that we do not have to manually select a suitable time step size. Moreover, for the more complex systems discussed in Sec. <ref> and Sec. <ref>, SMTmc is not able at all to provide a verification result in reasonable time, where we aborted the analysis after one hour.

Reachable set for the mobile robot, where the initial set is shown in white with a black border and the sets $\mathcal{S}_1$, $\mathcal{S}_2$, $\mathcal{S}_3$ are visualized in green and red.

§.§ Mobile Robot

One typical application for temporal logic is the formulation of tasks for mobile robots. Given the sets $\mathcal{S}_1$, $\mathcal{S}_2$, $\mathcal{S}_3$ visualized in Fig. <ref>, we consider the following reach-avoid task specified in natural language:

"Visit $\mathcal{S}_1$ during the first 4 seconds and afterward visit  $\mathcal{S}_2$
$~~~~~~$exactly 6 seconds after $\mathcal{S}_1$ while avoiding $\mathcal{S}_3$ at all times."

This corresponds to the following temporal logic specification:

\begin{equation*}
\varphi = \finally_{[0,4]} \big( x \in \mathcal{S}_1 \wedge \nextOp_{6} \, x \in \mathcal{S}_2 \big) \wedge \globally_{[0,10]} \, x \not \in \mathcal{S}_3.
\end{equation*}

We model the dynamics of the mobile robot by a double-integrator for the x- and y-position, which yields the system matrices $A = [\mathbf{0} ~ [I_2 \concat\mathbf{0}]]$ and $B = [\mathbf{0} \concat I_2]$. Moreover, we consider that the mobile robot tracks a reference trajectory $x_{\text{ref}}(t)$ that corresponds to the piecewise constant control input $u_{\text{ref}}(t)$ using the feedback control law $u_{\text{ctr}}(t) = u_{\text{ref}}(t) + K(x(t) - x_{\text{ref}}(t))$. The corresponding feedback matrix $K \in \R^{2 \times 4}$ is determined by applying an LQR control approach with state weighting matrix $Q = I_4$ and input weighting matrix $R = 0.1 \cdot I_2$ to the open-loop system. Overall, the dynamics of the controlled system is given as

\begin{equation*}
\begin{bmatrix} \dot x(t) \\ \dot x_{\text{ref}}(t) \end{bmatrix}
= \begin{bmatrix} A \hspace{-2pt} + \hspace{-2pt} BK & -BK \\ \mathbf{0} & A \end{bmatrix} \begin{bmatrix} x(t) \\ x_{\text{ref}}(t) \end{bmatrix}
+ \begin{bmatrix} B & B \\ B & \mathbf{0} \end{bmatrix} u(t).
\end{equation*}

Moreover, the initial set is $\initSet = [-0.1,0.1] \si{\meter} \times [-0.1,0.1]\si{\metre} \times \mathbf{0}$ and the set of uncertain inputs is $\mathcal{U} = u_{\text{ref}}(t) \times [-0.1,0.1]\si{\meter \per \square \second} \times [-0.1,0.1] \si{\meter \per \square \second}$.

The resulting reachable set shown in Fig. <ref> demonstrates that the mobile robot satisfies the temporal logic specification $\varphi$. However, since at no point in time the reachable set is fully contained inside the set $\mathcal{S}_1$, the original RTL approach <cit.> always conservatively classifies the system as unsafe, no matter how often we refine the tightness of the reachable set enclosure. On the other hand, our verifier summarized in Alg. <ref> fully automatically proves that the system is safe in only $2.89$ seconds.

Computation time in seconds for our automated verifier on different benchmarks from the ARCH competition, where $n$ is the number of states and $m$ the number of inputs of the system.

Benchmark      $n$     $m$     Specification     Safe?     Time

2*BLDC01     2*49     2*0     $\finally_{[0,0.2]}\big( \globally_{[0,0.18]}\, x_{40} < 0.0046 \vee \nextOp_{0.2} \,x_{40} > 0.006 \big)$         57.7

$\finally_{[0,0.2]}\big( \globally_{[0,0.19]}\, x_{40} < 0.0046 \vee \nextOp_{0.2} \, x_{40} > 0.006 \big)$         25.2
2*BLDF01     2*48     2*1     $\finally_{[0,0.2]}\big( \globally_{[0,0.18]}\, x_{40} < 0.0046 \vee \nextOp_{0.2} \,x_{40} > 0.006 \big)$         59.8

$\finally_{[0,0.2]}\big( \globally_{[0,0.19]}\, x_{40} < 0.0046 \vee \nextOp_{0.2} \, x_{40} > 0.006 \big)$         32.5
2*CBC01     2*201     2*0     $x_{110} < 69 ~ \until_{[0,0.004]} \, (x_{140} > 69 \vee x_{140} < -60)$         5.2

$x_{110} < 68 ~ \until_{[0,0.004]} \, (x_{140} > 69 \vee x_{140} < -60)$         9.1
2*CBF01     2*200     2*1     $x_{110} < 69 ~ \until_{[0,0.004]} \, (x_{140} > 69 \vee x_{140} < -60)$         17.9

$x_{110} < 68 ~ \until_{[0,0.004]} \, (x_{140} > 69 \vee x_{140} < -60)$         30.8
2*HEAT01     2*125     2*0     $\globally_{[4,8]} \, x_{10} > 0.1  \vee \finally_{[8,14]} \, x_{10} < 0.1$         9.1

$\globally_{[4,8]} \, x_{10} > 0.1  \vee \finally_{[8,13]} \, x_{10} < 0.1$          20.1
2*HEAT02     2*1000     2*0     $\globally_{[3,4]} \, x_{151} > 0.07  \vee \finally_{[4,6]} \, x_{151} < 0.07$         52.8

$\globally_{[3,4]} \, x_{151} > 0.07  \vee \finally_{[4,5]} \, x_{151} < 0.07$         50.8

§.§ ARCH Benchmarks

To demonstrate that our approach also scales to high-dimensional systems we consider the benchmarks from the 2022 ARCH competition <cit.>, which represent the limit of complexity that can be handled by state-of-the-art reachability tools.

For each benchmark, we introduce one safe and one unsafe temporal logic specification. To make the problems more challenging all specifications are close to the decision boundary, meaning that already slight modifications in the specification change the verification result from safe to unsafe or the other way round.

The original RTL approach <cit.> conservatively classifies all problem instances as unsafe, even those who satisfy the specification. As shown in Tab. <ref>, our automated verifier on the other hand correctly verifies and falsifies all benchmarks in under one minute, even though we consider complex temporal logic specifications with nested temporal operators and high-dimensional systems with up to 1000 states. Moreover, if the system is unsafe our verifier returns a falsifying trajectory that demonstrates the safety violation, as exemplary shown for the HEAT02 benchmark on the right side of Fig. <ref>. The left side of Fig. <ref> highlights the difference between our method and the original RTL approach <cit.>: The overall reachable set does neither satisfy the specification $\globally_{[3,4]}\, x_{151} > 0.07$ nor $\finally_{[4,6]}\, x_{151} < 0.07$, so that RTL conservatively classifies the system as unsafe. However, since all initial states that violate the specification $\globally_{[3,4]}\, x_{151} > 0.07$ satisfy the specification $\finally_{[4,6]}\, x_{151} < 0.07$, our method is able to prove that the system is safe. Finally, the decomposition of the computation time into the time spend on the different parts of Alg. <ref> shown in Fig. <ref> demonstrates that the percentage of time spend on each part heavily depends on the benchmark, and that no part clearly dominates the computation time for all benchmarks.

Reachable set (left) and falsifying trajectory (right) for the safe and the unsafe specification of the HEAT02 benchmark, where the dashed red line visualizes the boundary of the safe region and the dashed blue lines mark the points in time relevant for the temporal logic specification.

Computation time for our automated verifier on the ARCH benchmarks from Tab. <ref> partitioned into the time required for reachablity analysis (see Alg. <ref>), for verification (see Line <ref>-<ref> of Alg. <ref>), and for falsification (see Line <ref>-<ref> of Alg. <ref>).

§.§ Prediction of Traffic Participants

In addition to verifying dynamic systems against temporal logic specifications, our approach can also be used to cut away the parts of the reachable set which violate a temporal logic formula.

To demonstrate this, we consider set-based predication of traffic participants <cit.> in this section instead of a verification task. For set-based prediction of traffic participants, one first computes the space potentially occupied by other cars using reachability analysis, and then cuts away the regions that violate traffic rules. The trajectory of an autonomous car can then be planned to avoid the resulting regions, which guarantees safety at all times. In particular, we examine the traffic scenario shown in Fig. <ref>, which contains a no passing traffic sign. For this scenario the traffic rules defined by the no passing sign as well as the rule that the car is not allowed to leave the road can be formalized by the following temporal logic specification:

\begin{equation*}
\varphi = \globally_{[0,1]}\, (x_1 < 22 \vee x_2 < 2) \wedge \globally_{[0,1]} \, x_2 > -2 \wedge \globally_{[0,1]} \, x_2 < 6.
\end{equation*}

We model the dynamics of the car by a double-integrator for the x- and y-position, which yields the system $\dot x = A \, x + B \, u$ with $A = [\mathbf{0} ~ [I_2 \concat \mathbf{0}]]$ and $B = [\mathbf{0} \concat I_2]$. The uncertainty in the initial set $\initSet = [-0.1,0.1]\si{\meter} \times [-0.1,0.1]\si{\meter} \times [29.9,30.1]\si{\meter \per \second} \times [-0.1,0.1] \si{\meter \per \second}$ captures measurement errors in the position and velocity of the car. Moreover, the uncertain input $u(\cdot) \in \mathcal{U} = [-9,9] \si{\meter \per \square \second} \times [-9,9] \si{\meter \per \square \second}$ represents the unknown behavior of the driver, which is bounded by the physical limits of the car. For set-based prediction we execute Alg. <ref> with time step size $\Delta t = 0.05\si{\second}$ and truncation order $\kappa = 10$ for the negated specification $\neg \varphi$. The list $\mathbf{L}$ returned by the algorithm then contains all factor values that potentially satisfy the specification $\varphi$. A guaranteed enclosure of the set of all legal behaviors of the car is therefore given by the union of all polytopes in the list $\mathbf{L}$. Please note that the set resulting from the combination of the polytopes for the safe factors and the zonotopes for the reachable sets can be represented as a constrained zonotope <cit.>.

Visualization of set-based prediction of traffic participants showing the considered traffic scenario (top left), the overall reachable set (top right), the naive approach that simply cuts away unsafe regions (bottom left), and the set of all legal behaviors obtained with our approach (bottom right).

The resulting sets are visualized in Fig. <ref>. A naive approach is to simply cut away the unsafe regions $x_1 > 22 \wedge x_2 > 2$, $x_2 < -2$, and $x_2 > 6$. However, this ignores the dynamic behavior of the car and therefore yields a quite conservative result (Fig. <ref>, bottom left). Our approach, on the other hand, explicitly takes the dynamic behavior into account and therefore also cuts away regions that are guaranteed to reach an unsafe region in the future (Fig. <ref>, bottom right). This yields a smaller region still containing all legal behaviors, which increases the chances for finding a safe trajectory for the autonomous car.

§ FUTURE WORK

While we in this work focused on linear systems for simplicity, the general framework we presented can equivalently be applied for nonlinear and hybrid systems. Therefore, the extension to nonlinear and hybrid systems is our main focus for future work, and we outline some challenges as well as potential solutions in this section.

One key requirement for our automated verification framework is a reachability algorithm that preserves dependencies and converges to the exact reachable set if the algorithm parameters are tuned adequately. While dependency preserving reachability algorithms are readily available for nonlinear systems <cit.>, it is more complicated to preserve dependencies for hybrid systems due to the intersections with the guard sets at discrete transitions. One promising strategy to handle those discrete transitions in a dependency-preserving fashion is the guard-mapping approach <cit.>, since it circumvents the explicit computation of geometric intersections. Moreover, automated parameter tuning is significantly harder for nonlinear and hybrid systems compared to linear systems due to the increased number of algorithm parameters. A promising approach for automated parameter tuning for those systems has been published recently <cit.>, but does not yet guarantee convergence to the exact reachable set.

Another challenge arises from the fact that reachable sets for nonlinear and hybrid systems are in general non-convex, so that non-convex set representations such as polynomial zonotopes <cit.> or Taylor models <cit.> are required to compute tight enclosures. In our framework, replacing the convex zonotopes we used for linear systems with those non-convex set representations yields polynomial level sets instead of polytopes for the regions of unsafe zonotope factors. With polynomial level sets the checks for empty sets in Line <ref> of Alg. <ref> are computationally more expensive, and the optimization problem in (<ref>) becomes a mixed-integer polynomial program, which is harder to solve. Potential solutions to address these issues are using contractor programming <cit.> to check if a polynomial level set is empty, and applying continuous relaxation  <cit.>

to warm-start the mixed-integer polynomial program.

§ CONCLUSION

We presented the first fully-automated verifier for linear systems that considers the very general case of specifications defined by signal temporal logic formulas. Our algorithm avoids the conservatism of the previous reachset temporal logic approach by explicitly keeping track of which initial states and which uncertain inputs satisfy or violate single parts of the formula, and is therefore guaranteed to always find the correct solution in finite time for decidable verification problem instances. As we demonstrated with numerical experiments, our automated verifier can even solve high-dimensional benchmarks with up to 1000 states as well as complex temporal logic specifications with nested temporal operators very efficiently. Other advantageous features of our algorithm are that it returns a falsifying trajectory in case of a safety violation and that it naturally divides the initial as well as the input set into parts which satisfy or violate the specification. This is for instance beneficial for set-based prediction of traffic participants, as we demonstrated with an exemplary traffic scenario.

Acknowledgements. This material is based upon work supported by the Air Force Office of Scientific Research and the Office of Naval Research under award numbers FA9550-19-1-0288, FA9550-21-1-0121, FA9550-23-1-0066 and N00014-22-1-2156. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the United States Air Force or the United States Navy.

We now demonstrate the conversion from signal temporal logic to reachset temporal logic according to <cit.> for the exemplary STL formula

\begin{equation*}
\varphi = \nextOp_{0.5} \, x_1 > 2 \vee \neg \finally_{[0,0.8]} \, x_2 \leq 3
\end{equation*}

and time step size $\Delta t = 0.5$. The first step is to convert the formula to negation normal form by moving all negations inward until they only appear in front of non-temporal expressions:

\begin{equation*}
\varphi_{\text{nnf}} = \nextOp_{0.5} \, \underbrace{x_1 > 2}_{\pred_1} \vee \, \globally_{[0,0.8]} \, \underbrace{x_2 > 3}_{\pred_2}.
\end{equation*}

Next, we convert the formula to sampled-time STL with sampling period $\Delta t$. For this, we rewrite all times and time intervals for the temporal operators as integer multiples $i \, \Delta t$, $i \in \mathbb{N}_0$ of the time step size. If the start and end times are not divisible by the time step size we can either extend or shorten the corresponding time intervals in a sound matter. In our formula, for example, we can replace $\globally_{[0,0.8]} \, \pred_2$ by $\globally_{[0,1]} \, \pred_2$ since satisfaction of $\globally_{[0,1]} \, \pred_2$ implies satisfaction of $\globally_{[0,0.8]} \, \pred_2$. For time step size $\Delta t = 0.5$ we therefore obtain

\begin{equation*}
\varphi_{\text{st}_{\text{1}}} = \nextOp_{1} \, \pred_1 \vee \globally_{[0,2]} \, \pred_2,
\end{equation*}

which after applying the rewriting rules in <cit.> results in the sampled-time STL formula

\begin{equation*}
\varphi_{\text{st}_{\text{2}}} = \nextOp_{1} \pred_1 \vee \big( \hspace{-2pt} \nextOp_0 \hspace{-1pt}\pred_2 \wedge \globally_{[0,1]} \, \pred_2 \wedge \nextOp_1  \pred_2 \wedge \nextOp_1 \, \globally_{[0,1]} \, \pred_2 \wedge \nextOp_2 \, \pred_2 \big).
\end{equation*}

Afterward, we have to convert the formula to conjunctive normal form, which yields

\begin{equation*}
\begin{split}
\varphi_{\text{cnf}} = \big( \hspace{-2pt} & \nextOp_{1} \hspace{-1pt}\pred_1 \vee \nextOp_0 \, \pred_2 \big) \wedge \big(\hspace{-2pt} \nextOp_{1} \hspace{-2pt}\pred_1 \vee \globally_{[0,1]} \, \pred_2 \big) \wedge \\
& \nextOp_{1} \hspace{-1pt}(\pred_1 \vee \pred_2 ) \wedge \big(\hspace{-2pt} \nextOp_{1} \hspace{-1pt}\pred_1 \vee \nextOp_1 \, \globally_{[0,1]} \, \pred_2 \big) \wedge \big( \hspace{-2pt} \nextOp_{1}\hspace{-1pt} \pred_1 \vee \nextOp_2 \, \pred_2 \big).
\end{split}
\end{equation*}

Finally, we can apply the conversion to reachset temporal logic in <cit.> to obtain

\begin{align*}
\varphi_{\text{rtl}} = \big( & \hspace{-2pt} \nextOp_{1} \hspace{-2pt}\all \pred_1 \vee \nextOp_0 \all \pred_2 \big) \wedge \big( \hspace{-2pt} \nextOp_{1} \hspace{-2pt}\all \pred_1 \vee \nextOp_{0.5} \all \pred_2 \big) \wedge \\
& \hspace{-2pt} \nextOp_{1} \hspace{-2pt}\all (\pred_1 \vee \pred_2 ) \wedge \big( \hspace{-2pt} \nextOp_{1} \hspace{-2pt}\all \pred_1 \vee \nextOp_{1.5} \all \pred_2 \big) \wedge \big( \hspace{-2pt} \nextOp_{1}\hspace{-2pt} \all \pred_1 \vee \nextOp_2 \all \pred_2 \big).
\end{align*}

Moreover, according to <cit.> entailment can be equivalently formulated in terms of intersection checks with polytopes as in (<ref>):

\begin{align*}
\mathcal{R} \vDash \varphi_{\text{rtl}} ~~ \Leftrightarrow ~~ & \big( \mathcal{R}(t_1)  \cap \mathcal{P}_1 = \emptyset \vee \mathcal{R}(t_0) \cap \mathcal{P}_2 = \emptyset \big) \wedge \\
& \big( \mathcal{R}(t_1) \cap \mathcal{P}_1 = \emptyset \vee \mathcal{R}(\tau_0) \cap \mathcal{P}_2 = \emptyset \big) \wedge \\
& ~ \mathcal{R}(t_1) \cap \mathcal{P}_{12} = \emptyset \wedge \\
&    \big( \mathcal{R}(t_1) \cap \mathcal{P}_1 = \emptyset \vee \mathcal{R}(\tau_1) \cap \mathcal{P}_2 = \emptyset \big) \wedge \\
&  \big( \mathcal{R}(t_1) \cap \mathcal{P}_1 = \emptyset \vee \mathcal{R}(t_2) \cap \mathcal{P}_2 = \emptyset \big),
\end{align*}

where the polytopes $\mathcal{P}_1$, $\mathcal{P}_2$, and $\mathcal{P}_{12}$ corresponding to the predicates $\pred_1$, $\pred_2$, and $\pred_{1} \vee \pred_{2}$ are defined as

\begin{equation*}
\mathcal{P}_1 = \big \langle [1~0],2 \big \rangle_P, ~~ \mathcal{P}_2 = \big \langle [0~1],3 \big \rangle_P, ~~ \mathcal{P}_{12} = \bigg \langle \begin{bmatrix} 1~ & ~0 \\ 0~ & ~1 \end{bmatrix},\begin{bmatrix} 2 \\ 3 \end{bmatrix} \bigg \rangle_P.
\end{equation*}

The resulting formula can therefore be directly evaluated on the reach sequence $\mathcal{R}(t_0),\mathcal{R}(\tau_0),\mathcal{R}(t_1),\mathcal{R}(\tau_1),\mathcal{R}(t_2)$.

We now prove that the enclosure computed with the reachability algorithm in Alg. <ref> converges to the exact reachable set for $\Delta t \to 0$. In summary, Alg. <ref> contains three sources of over-approximation:

* The enclosure of the difference between the reachable set due to constant and time-varying inputs $\mathcal{D}$.

* The curvature enclosure $\mathcal{C}_i$.

* The over-approximation in the zonotope enclosure of the convex hull
$conv(\mathcal{H}(t_i),\mathcal{H}(t_{i+1}))$ computed using (<ref>).

We therefore have to show that all these over-approximation errors converge to zero for $\Delta t \to 0$. We begin with the difference between constant and time-varying inputs:

The enclosure of the difference between the reachable set due to constant and time-varying inputs $\mathcal{D}$ according to (<ref>) satisfies

\begin{equation*}
\lim_{\Delta t \to 0} \mathcal{D} = \mathbf{0}.
\end{equation*}

With the formula for $\mathcal{D}$ in (<ref>), we obtain

\begin{align*}
\lim_{\Delta t \to 0} \mathcal{D} & \overset{\eqref{eq:inputDiff}}{=}	\lim_{\Delta t \to 0} \bigg( \sum_{j=1}^\kappa \frac{A^j \Delta t^{j+1}}{(j+1)!} \bigg) \mathcal{U}_0 \oplus \bigoplus_{j=1}^\kappa \frac{A^j \Delta t^{j+1}}{(j+1)!}\, \mathcal{U}_0 \oplus 2 \,  \Delta t \, \mathbfcal{E} \, \mathcal{U}_0 \\
& ~= \bigg( \sum_{j=1}^\kappa \frac{A^j \, 0}{(j+1)!} \bigg) \mathcal{U}_0 \oplus \bigoplus_{j=1}^\kappa \frac{A^j \, 0}{(j+1)!}\, \mathcal{U}_0 \oplus 2 \,  0 \, \mathbfcal{E} \, \mathcal{U}_0 = \mathbf{0},
\end{align*}

which concludes the proof.

Next, we consider the curvature enclosure:

The curvature enclosure $\mathcal{C}_i$ according to (<ref>) satisfies

\begin{equation*}
\lim_{\Delta t \to 0} \mathcal{C}_i = \mathbf{0}.
\end{equation*}

According to <cit.>, it holds that

\begin{equation}
\lim_{\Delta t \to 0} \mathbfcal{F} = \mathbf{0} ~~~~ \text{and} ~~~~ \lim_{\Delta t \to 0} \mathbfcal{G} = \mathbf{0},
\label{eq:proofC}
\end{equation}

with $\mathbfcal{F}$ and $\mathbfcal{G}$ defined as in (<ref>). With the formula for $\mathcal{C}_i$ in (<ref>) we therefore obtain

\begin{align*}
\lim_{\Delta t \to 0} \mathcal{C}_i & \overset{\eqref{eq:Cu}}{=} \lim_{\Delta t \to 0} (e^{A \Delta t})^{i} \big( \mathbfcal{F} \, \mathcal{X}_0 \oplus \mathbfcal{G} \, \widetilde{u} \big) = \lim_{\Delta t \to 0} \underbrace{(e^{A \, 0})^{i}}_{=I_n} \big( \mathbfcal{F} \, \mathcal{X}_0 \oplus \mathbfcal{G} \, \widetilde{u} \big) \\
& ~\, = \lim_{\Delta t \to 0}  \mathbfcal{F} \, \mathcal{X}_0 \oplus \mathbfcal{G} \, \widetilde{u} \overset{\eqref{eq:proofC}}{=} \mathbf{0},
\end{align*}

which concludes the proof.

Finally, we examine the over-approximation introduced by the enclosure of the convex hull:

Given two zonotopes $\mathcal{Z}_1,\mathcal{Z}_2 \subset \mathbb{R}^n$, let $\text{\normalfont\texttt{diff}}(\mathcal{Z}_1,\mathcal{Z}_2)$ denote the difference between the exact convex hull $conv(\mathcal{Z}_1,\mathcal{Z}_2)$ as defined in (<ref>) and the enclosure computed according to (<ref>). The homogeneous solution $\mathcal{H}(t_{i})$ in (<ref>) satisfies:

\begin{equation*}
\lim_{\Delta t \to 0} \text{\normalfont\texttt{diff}}\big(\mathcal{H}(t_i),\mathcal{H}(t_{i+1})\big) = \mathbf{0}.
\end{equation*}

For the limit case $\Delta t \to 0$, we obtain the following for the convex hull $conv \big(\mathcal{H}(t_i),\mathcal{H}(t_{i+1}) \big)$:

\begin{equation}
\begin{split}
\lim_{\Delta t \to 0} conv \big(\mathcal{H}(t_i),\mathcal{H}(t_{i+1}) \big) &\overset{\eqref{eq:propScheme}}{=} \lim_{\Delta t \to 0} conv \big(\mathcal{H}(t_{i}), e^{A \Delta t} \, \mathcal{H}(t_i) \oplus T \, \widetilde{u} \big) \\
&~= \lim_{\Delta t \to 0} conv \big(\mathcal{H}(t_{i}), e^{A \Delta t} \, \mathcal{H}(t_i) \oplus A^{-1} (e^{A \Delta t} - I_n) \, \widetilde{u} \big) \\
&~=  conv \big(\mathcal{H}(t_{i}), \underbrace{e^{A \, 0}}_{=I_n} \, \mathcal{H}(t_i) \oplus A^{-1} (\underbrace{e^{A \, 0}}_{=I_n} - I_n) \, \widetilde{u} \big) \\
&~= conv \big(\mathcal{H}(t_{i}),\mathcal{H}(t_i)\big).
\end{split}
\label{eq:proofConv1}
\end{equation}

Moreover, given a set $\mathcal{S} \subset \Rn$, the exact convex hull as defined in (<ref>) satisfies

\begin{equation}
conv(\mathcal{S},\mathcal{S}) = \bigg \{\sum_{i=1}^{n+1} \lambda_i\,s_i~\bigg|~s_i \in \underbrace{\mathcal{S} \cup \mathcal{S}}_{=\mathcal{S}},~\lambda_i \geq 0,~\sum_{i=1}^{n+1} \lambda_i = 1 \bigg\} = \mathcal{S}.
\label{eq:proofConv2}
\end{equation}

In addition, given a zonotope $\mathcal{Z} = \langle c,G \rangle_Z$, the enclosure of the convex hull according to (<ref>) satisfies

\begin{equation}
\begin{split}
conv\big(\mathcal{Z},\mathcal{Z} \big) & \overset{\eqref{eq:zonoConvHull}}{=} \big\langle 0.5(c+c),\big[0.5(G+G)~0.5(G-G)~0.5(c-c) \big]\big\rangle_Z  \\[5pt]
& = \big \langle c, \big[G~ \mathbf{0} ~ \mathbf{0}\big] \big \rangle_Z = \mathcal{Z}.
\end{split}
\label{eq:proofConv}
\end{equation}

For $conv \big(\mathcal{H}(t_{i}),\mathcal{H}(t_i)\big)$ in (<ref>), both, the exact convex hull as defined in (<ref>) as well as the zonotope enclosure computed using (<ref>) yield the result $\mathcal{H}(t_i)$ according to (<ref>) and (<ref>), which proves that the difference between the exact convex hull and the enclosure computed according to (<ref>) converges to 0 for $\Delta t \to 0$.


* Simon (2019) Simon J. D., 2019, ARA&A, 57, 375
  * Simons et al. (2020) Simons R. C., et al., 2020, ApJ, 905, 167
  * Skillman et al. (2013) Skillman E. D., et al., 2013, AJ, 146, 3
  * Sormani et al. (2018) Sormani M. C., Sobacchi E., Pezzulli G., Binney J., Klessen R. S., 2018, MNRAS, 481, 3370
  * Stanimirović et al. (2002) Stanimirović S., Dickey J. M., Krčo M., Brooks A. M., 2002, ApJ, 576, 773
  * Teklu et al. (2015) Teklu A. F., Remus R.-S., Dolag K., Beck A. M., Burkert A., Schmidt A. S., Schulze F., Steinborn L. K., 2015, ApJ, 812, 29
  * Tepper-García & Bland-Hawthorn (2018) Tepper-García T., Bland-Hawthorn J., 2018, MNRAS, 478, 5263
  * Tepper-García et al. (2015) Tepper-García T., Bland-Hawthorn J., Sutherland R. S., 2015, ApJ, 813, 94
  * Tepper-García et al. (2019) Tepper-García T., Bland-Hawthorn J., Pawlowski M. S., Fritz T. K., 2019, MNRAS, 488, 918
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * Thom et al. (2008) Thom C., Peek J. E. G., Putman M. E., Heiles C., Peek K. M. G., Wilhelm R., 2008, ApJ, 684, 364
  * Walker et al. (2007) Walker M. G., Mateo M., Olszewski E. W., Gnedin O. Y., Wang X., Sen B., Woodroofe M., 2007, ApJ, 667, L53
  * Weisz et al. (2014) Weisz D. R., Dolphin A. E., Skillman E. D., Holtzman J., Gilbert K. M., Dalcanton J. J., Williams B. F., 2014, ApJ, 789, 147
  * White & Frenk (1991) White S. D. M., Frenk C. S., 1991, ApJ, 379, 52
  * Wright et al. (2020) Wright R. J., Lagos C. d. P., Power C., Mitchell P. D., 2020, MNRAS, 498, 1668
  * Yamasaki & Totani (2020) Yamasaki S., Totani T., 2020, ApJ, 888, 105
  * Yao et al. (2009) Yao Y., Wang Q. D., Hagihara T., Mitsuda K., McCammon D., Yamasaki N. Y., 2009, ApJ, 690, 143
  * Yoshino et al. (2009) Yoshino T., et al., 2009, PASJ, 61, 805
  * Zhao et al. (2009) Zhao D. H., Jing Y. P., Mo H. J., Börner G., 2009, ApJ, 707, 354
  * de Boer et al. (2012a) de Boer T. J. L., et al., 2012a, A&A, 539, A103
  * de Boer et al. (2012b) de Boer T. J. L., et al., 2012b, A&A, 544, A73
  * de Boer et al. (2014) de Boer T. J. L., Tolstoy E., Lemasle B., Saha A., Olszewski E. W., Mateo M., Irwin M. J., Battaglia G., 2014, A&A, 572, A10
  * van Dokkum et al. (2013) van Dokkum P. G., et al., 2013, ApJ, 771, L35
  * van den Bosch et al. (2014) van den Bosch F. C., Jiang F., Hearin A., Campbell D., Watson D., Padmanabhan N., 2014, MNRAS, 445, 1713

##  Appendix A Growth of the corona

The initial corona at redshift $z\gtrsim 2$ would have been much denser and
more compact than today as it has since expanded as the virial radius of the
MW has grown. While the density is expected to have decreased, the total mass
of the corona has grown through infall of gas from the IGM that is unable to
reach the disc as well as hot gas expelled from the MW’s ISM by supernovae.

We refer to all the gas within the virial radius but outside of the ISM as
‘circumgalactic medium’ (CGM). From mass conservation it then follows that the
mass in the CGM must change according to the sum of four non-negative terms:

$\dot{M}_{\mathrm{CGM}}=\dot{M}_{\mathrm{IGM\rightarrow
CGM}}+\dot{M}_{\mathrm{ISM\rightarrow CGM}}-\dot{M}_{\mathrm{CGM\rightarrow
IGM}}-\dot{M}_{\mathrm{CGM\rightarrow ISM}}.$ (23)

The arrows indicate the direction of the mass flow. These terms represent
inflow of gas from the IGM, gas ejected from the ISM, gas ejected from the
CGM, and gas accreted onto the ISM, respectively.

While some of the gas in the CGM at any point in time will be cold, any cold
cloud or filament cannot remain in the CGM. It will either be heated or fall
onto the disc. Gas that remained cold throughout its journey onto the disc is
commonly referred to as ‘cold-mode accretion’ while gas that was heated and
then later cooled to eventually make its way to the disc is referred to as
‘hot-mode accretion’ (see e.g. Birnboim & Dekel, 2003; Kereš et al., 2005).
Thus, cold gas will not accumulate in the CGM over longer time-scales. This is
supported also by the results of cosmological zoom-in simulations which find
that the fraction of the CGM mass that is in hot gas increases substantially
with time (Hafen et al., 2019). Considering then the change in mass over a
time $\Delta t$ much longer than the time-scale for cold gas infall (which is
estimated to be of the order of $\lesssim 2$ Gyr at $z\leq 2$; see Nelson et
al. (2015)) we can equate the change in the CGM mass over that time with the
change in the hot corona mass $\Delta M_{\mathrm{cor}}$. Hence, from eq. (23)
we can derive

$\Delta M_{\mathrm{cor}}=f_{\mathrm{cor}}\Delta M_{\mathrm{IGM\rightarrow
CGM}},$ (24)

where

$f_{\mathrm{cor}}\equiv 1-\frac{\Delta M_{\mathrm{CGM\rightarrow ISM}}+\Delta
M_{\mathrm{CGM\rightarrow IGM}}-\Delta M_{\mathrm{ISM\rightarrow CGM}}}{\Delta
M_{\mathrm{IGM\rightarrow CGM}}}.$ (25)

$f_{\mathrm{cor}}<0$ would imply that the corona is losing mass while
$f_{\mathrm{cor}}>1$ would imply that $\Delta M_{\mathrm{ISM\rightarrow
CGM}}>\Delta M_{\mathrm{CGM\rightarrow ISM}}$ and thus the disc would be
ejecting more gas than it would accrete. Because the mass in both the corona
and the disc of the MW should be growing, $f_{\mathrm{cor}}$ is hence in the
range $0\leq f_{\mathrm{cor}}\leq 1$ and represents the fraction of baryons
accreted onto the halo over $\Delta t$ that end up in the corona.
$f_{\mathrm{cor}}$ is related to the fraction of baryons accreted onto the
halo that ends up in the disc (as either ISM or stars):

$f_{\mathrm{disc}}\equiv\frac{\Delta M_{\mathrm{CGM\rightarrow ISM}}-\Delta
M_{\mathrm{ISM\rightarrow CGM}}}{\Delta M_{\mathrm{IGM\rightarrow CGM}}}.$
(26)

If no gas is expelled from the halo, i.e. $\Delta M_{\mathrm{CGM\rightarrow
IGM}}=0$, then $f_{\mathrm{cor}}=1-f_{\mathrm{disc}}$. Generally,
$f_{\mathrm{cor}}=1-f_{\mathrm{disc}}-f_{\mathrm{out}}$ where

$f_{\mathrm{out}}\equiv\frac{\Delta M_{\mathrm{CGM\rightarrow IGM}}}{\Delta
M_{\mathrm{IGM\rightarrow CGM}}}.$ (27)

This is the fraction of gas accreted onto the halo that was later ejected from
it.

$f_{\mathrm{disc}}$ is related to the _instantaneous_ (i.e. in terms of
immediate mass changes $\dot{M}$ rather than long time scale changes $\Delta
M$) fraction of gas accreted onto the halo that reaches the disc $\zeta\geq
f_{\mathrm{disc}}$. $\zeta$ is generally greater than $f_{\mathrm{disc}}$
because it only concerns infall onto the disc and so does not have the outflow
term included in $f_{\mathrm{disc}}$. This parameter is an input for many
analytical ‘bath-tub’ equilibrium models of galaxy evolution (e.g. Davé et
al., 2012; Lilly et al., 2013; Dekel & Mandelker, 2014). It is generally found
to be in the range $0.3-0.5$ for MW progenitors at $z<2$ in theoretic models
and simulations (Davé et al., 2012; Lilly et al., 2013; Mitchell et al.,
2020).

A standard assumption is that the gas infall from the IGM equals the DM
accretion scaled by the cosmological baryon fraction $f_{\mathrm{b}}$. Hence,
we may write eq. (24) as

$\Delta M_{\mathrm{cor}}=f_{\mathrm{cor}}f_{\mathrm{b}}\Delta
M_{\mathrm{DM}}.$ (28)

This is convenient since the growth of DM haloes is a robust prediction of
cosmological simulations that agrees well between different simulations
regardless of differences in their gas physics. Cosmological hydrodynamical
zoom in simulations with AGN feedback suggest that eq. (28) substantially
overestimates the baryon infall at low halo masses because strong feedback is
able to affect gas outside the virial radius and largely prevent it from being
accreted onto the halo (Nelson et al., 2015; Hafen et al., 2019; Wright et
al., 2020). However, the mass below which this effect becomes important
depends on details in the feedback implementation and it is not clear if this
effect is significant at the $M_{\mathrm{vir}}\gtrsim 10^{11.5}M_{\odot}$
virial masses at $z<2$ that we consider.