Quantifying Conformance Using the Skorokhod Metric

A fundamental question in model-based design is conformance testing: whether two models of a system display similar behavior. For discrete systems, this question is well-studied [19, 20, 28, 29], and there is a rich theory of process equivalences based, e.g., on bisimilarity. For continuous and hybrid systems, however, the state of the art is somewhat unsatisfactory. While there is a straightforward generalization of process equivalences to the continuous case, in practice, equivalence notions such as bisimilarity are always too strong and most systems are not bisimilar. Since equivalence is a Boolean notion, one gets no additional information about the systems other than they are “not bisimilar.” Further, even if two dynamical systems are bisimilar, they may still differ in many control-theoretic properties. Thus, classical notions for equivalence and conformance have been of limited use in industrial practice.

In recent years, the notion of bisimulation has therefore been generalized to metrics on systems, which quantify the distance between them. For example, one approach is that of \(\epsilon \)-bisimulation, which requires that the states of the two systems remain “close” forever (within an \(\epsilon \)-ball), rather than coincide exactly. Under suitable stability assumptions on the dynamics, one can construct \(\epsilon \)-bisimulations [17, 18]. Unfortunately, proving the pre-requisites for the existence of \(\epsilon \)-bisimulations for complex dynamical models, or coming up with suitable and practically tractable bisimulation functions is extremely difficult in practice. In addition, establishing \(\epsilon \)-bisimulation requires full knowledge of the system dynamics making the scheme inapplicable where one system is an actual physical component with unknown dynamics. So, these notions have also been of limited industrial use so far.

In addition, there is the more informal requirement of usability: the metric should classify systems which the engineers consider close as being close, and conversely.

The simplest candidate metric is a pointwise metric that computes the maximum pointwise difference between two trajectories, sometimes generalized to apply a constant time-shift to one trajectory [15]. Unfortunately, for many practical models, two trajectories may be close only under variable time-shifts. This is the case, for example, for two dynamical models that may use different numerical integration techniques (e.g., fixed step versus adaptive step) or when some component in the implementation has some jitter. Thus, the pointwise metric spuriously reports large distances for “close” models. More nuanched hybrid distances have been proposed [1], but the transference properties of these metrics w.r.t. common temporal logics are not yet clear.

Transference. We show that the Skorokhod metric gives a robust quantification of system conformance by relating the metric to TLTL (timed LTL) enriched with (i) predicates of the form \(f(x_1,\dots , x_n) \ge 0\), as in Signal Temporal Logic [15], for specifying constraints on trace values; and (ii) freeze quantifiers, as in TPTL  [4], for specifying temporal constraints (freeze quantifiers can express more complex timing constraints than bounded timing constraints, e.g., of MTL). TLTL subsumes MTL and STL  [15]. We prove a transference theorem: flows (and propositional traces) which are close under the Skorokhod metric satisfy “close” TLTL formulae for a rich class of temporal and spatial predicates, where the untimed structure of the formulae remains unchanged, only the predicates are enlarged.

Tractability. We improve on recent polynomial-time algorithms for the Skorokhod metric [25] by taking advantage of the fact that, in practice, only retimings that map the times in one trace to “close” times in the other are of interest. This enables us to obtain a streaming sliding-window based monitoring procedure which takes only O(W) time per sample, where W is the window size (assuming the dimension n of the system to be a constant).

Usability. Using the Skorokhod distance checking procedure as a subroutine, we have implemented a Simulink toolbox for conformance testing. Our tool integrates with Simulink’s model-based design flow for control systems, and provides a stochastic search-based approach to find inputs which maximize the Skorokhod distance between systems under these inputs.

We present three case studies from the control domain, including industrial challenge problems; our empirical evaluation shows that our tool computes sharp estimates of the conformance distance reasonably fast on each of them. Our input models were complex enough that techniques such as \(\epsilon \)-bisimulation functions are inapplicable. We conclude that the Skorokhod metric can be an effective foundation for semi-formal conformance testing for complex dynamical models. Proofs of the theorems are given in the accompanying technical report [13].

Related Work. The work of [1, 2] is closely related to ours. In it, robustness properties of hybrid state sequences are derived with respect to a trace metric which also quantifies temporal and spatial variations. Our work differs in the following ways. First, we guarantee robustness properties over flows rather than only over (discrete) sequences. Second, the Skorokhod metric is a stronger form of the \((T,J,(\tau , \epsilon ))\)-closeness degree^1,2 (for systems which do not have hybrid time); and allows us to give stronger robustness transference guarantees. The Skorokhod metric requires order preservation of the timeline, which the \((T,J,(\tau , \epsilon ))\)-closeness function does not. Preservation of the timeline order allows us to (i) keep the untimed structure of the formulae the same (unlike in the transference theorem of [1]); (ii) show transference of a rich class of global timing constraints using freeze quantifiers (rather than only for the standard bounded time quantifiers of MTL/MITL). However, for implementations where the timeline order is not preserved, we have to settle for the less stronger guarantees provided by [1]. The work of [15] deals with spatial robustness of STL; the only temporal disturbances considered are constant time-shifts for the entire signal where the entire signal is moved to the past, or to the future by the same amount. In contrast, the Skorokhod metric incorporates variable time-shifts.

In this section, we demonstrate a transference result for the Skorokhod metric for a version of the timed linear time logic TLTL  [4]. The logic we consider generalizes MTL and STL. We show that if the Skorokhod distance between two traces is small, they satisfy close TLTL formulae. Given a formula \(\phi \) of TLTL satisfied by trace \(\pi _1\), we can compute a “relaxation” of \(\phi \) that will be satisfied by the “close” trace \(\pi _2\). We first present the results in a propositional framework, and then extend to \({\mathbb {R}}^n\)-valued spaces for a logic generalizing STL.

3.3 Transference of TLTL Properties for \({\mathbb {R}}^n\)-valued Signals

A timed\({\mathbb {R}}^n\)-valued trace\(\pi \) is a function from a closed interval I of \({\mathbb {R}}_+\) to \({\mathbb {R}}^n\). For \( \overline{\alpha } = (\alpha ^0, \dots , \alpha ^n)\in {\mathbb {R}}^n\), we denote the k-th dimensional value \(\alpha ^k\) as \(\overline{\alpha }[k]\). The \(\pi \) projected function onto the k-th \({\mathbb {R}}\) dimension is denoted by \(\pi _{k}: I \mapsto {\mathbb {R}}\).

To define the semantics of TLTL formulae over timed \({\mathbb {R}}^n\)-valued sequences, we use booleanizing predicates \(\mu : {\mathbb {R}}^n \mapsto {\mathbb {B}}\), as in STL  [15], to transform \({\mathbb {R}}^n\)-valued sequences into timed propositional sequences. These predicates are part of the logical specification. In this work, we restrict our attention to traces and predicates such that each predicate varies only finitely often on the finite time traces under consideration. Since we also have freeze variables, TLTL with predicates is strictly more expressive than STL⁵ (as in the propositional case [6]).

Definition 5

(TLTL\((\mathcal {F}_{{\mathsf {T}}}, \mathcal {F}_{{\mathsf {S}}})\)Syntax). Given a set of variables \(V_{\mathsf {T}}\) (the freeze variables), a set of ordered variables \(V_{\mathsf {S}}\) (the signal variables), and two sets \(\mathcal {F}_{{\mathsf {T}}}, \mathcal {F}_{{\mathsf {S}}}\) of functions, the formulae of TLTL\((\mathcal {F}_{{\mathsf {T}}}, \mathcal {F}_{{\mathsf {S}}})\) are defined by the grammar:

$$\begin{aligned} \phi := {true}\mid f_{\mathsf {T}}(\overline{x}) \sim 0 \mid f_{\mathsf {S}}(\overline{y}) \sim 0 \mid \lnot \phi \mid \phi _1 \wedge \phi _2 \mid \phi _1\vee \phi _2\mid \phi _1 {\,\mathcal{U}\,}\phi _2 \mid x.\phi \quad \text { where} \end{aligned}$$

• \(x\in V_{\mathsf {T}}\), and \(\overline{x} = (x_1, \dots , x_l)\) with \(x_i\in V_{\mathsf {T}}\) for all \(1\le i\le l\);

• \(\overline{z} = (z_1, \dots , z_d)\) with \(z_j\in V_{\mathsf {S}}\) for all \(1\le j \le d\) (with \(d\le n\));

• \(V_{\mathsf {T}} \) and \(V_{\mathsf {S}}\) are disjoint;

• \(f_{\mathsf {T}} \in \mathcal {F}_{{\mathsf {T}}}\) and \(f_{\mathsf {S}} \in \mathcal {F}_{{\mathsf {S}}}\) are real-valued functions, and \(\sim \) is \( \le , <, \ge , \) or \(>\). \(\quad \square \)

The semantics of TLTL\((\mathcal {F}_{{\mathsf {T}}}, \mathcal {F}_{{\mathsf {S}}})\) is straightforward and similar to the propositional case (Definition 3). The only new ingredients are the booleanizing predicates \(f_{\mathsf {S}}(\overline{z}) \sim 0 \): we define \(\pi \models _{\mathcal {E}} f_{\mathsf {S}}(z_1, \dots , z_d) \sim 0\) iff \( f_{\mathsf {S}}( \pi _{j_1}[t_0], \dots , \pi _{j_d}[t_0]) \sim 0\) for any freeze variable environment \(\mathcal {E}\), where \(t_0 = \min {\mathsf {tdom}}(\pi )\), and \(z_i\) is the \(j_i\)-th variable in \(V_{\mathsf {S}}\) (i.e., \(z_i\) refers to the \(j_i\)-th dimension in the signal trace). We require that for a timed \({\mathbb {R}}^n\)-valued trace \(\pi \) to satisfy \(\phi \), the arity of the functions in \( \mathcal {F}_{{\mathsf {S}}}\) occurring in \(\phi \) should not be more than n, that is, functions should not refer to dimensions greater than n for an \({\mathbb {R}}^n\) trace.

\(\delta \)relaxation ofTLTL\((\mathcal {F}_{{\mathsf {T}}}, \mathcal {F}_{{\mathsf {S}}})\). Let \( \mathbf {J}_{V_{\mathsf {S}}}\) be a mapping from \(V_{\mathsf {S}}\) to closed intervals of \({\mathbb {R}}\), thus \(\mathbf {J}_{V_{\mathsf {S}}}(z)\) denotes a sub-domain of \(z\in V_{\mathsf {S}}\). The relaxation function \({\mathsf {rx}}_{{J}, \mathbf {J}_{V_{\mathsf {S}}}}^{\delta }\) which operates on TLTL\((\mathcal {F}_{{\mathsf {T}}}, \mathcal {F}_{{\mathsf {S}}})\) formulae is defined analogous to the relaxation function \({\mathsf {rx}}_{J}^{\delta }\) in Definition 4. We omit the similar cases, and only present the new case for the predicates formed from \( \mathcal {F}_{{\mathsf {S}}}\).

$$\begin{aligned} {\mathsf {rx}}_{J, \mathbf {J}_{V_{\mathsf {S}}}}^{\delta }\left( f_{\mathsf {S}}(z_1, \dots , z_l) \right) \sim 0) \ = {\left\{ \begin{array}{ll} f_{\mathsf {S}}(z_1, \dots , z_l) \, +\, K^{f_{\mathsf {S}}}_{J, \mathbf {J}_{V_{\mathsf {S}}}}(\delta ) \ \sim \, 0 &{} \text { if } \sim \, \in \{ >, \ge \};\\ f_{\mathsf {S}}(z_1, \dots , z_l) \, -\, K^{f_{\mathsf {S}}}_{J, \mathbf {J}_{V_{\mathsf {S}}}}(\delta ) \ \sim \, 0 &{} \text { if } \sim \,\in \{ <, \le \} \end{array}\right. } \end{aligned}$$

where \(K^{f_{\mathsf {S}}}_{J, \mathbf {J}_{V_{\mathsf {S}}}}: \big [0,\ \max _{z\in V_{\mathsf {S}}} |\max \mathbf {J}_{V_{\mathsf {S}}}(z) \, - \, \min \mathbf {J}_{V_{\mathsf {S}}}(z)| \big ] \mapsto {\mathbb {R}}_+\) is a function s.t.

$$\begin{aligned} K^{f_{\mathsf {S}}}_{J, \mathbf {J}_{V_{\mathsf {S}}}}(\delta ) = \! \sup _{\begin{array}{c} u_i\in \mathbf {J}_{V_{\mathsf {S}}}(z_i); \, u'_i\in \mathbf {J}_{V_{\mathsf {S}}}(z'_i)\\ \text { for all } i \end{array}} \!\! \left\{ \! \left| \begin{array}{c} f_{\mathsf {S}}(u_1,\dots , u_l) \\ -\\ f_{\mathsf {S}}(u_1',\dots ,u_l') \end{array}\right| \text { s.t. } |u_i \!-\!u_i'| \le \delta \text { for all } i \right\} \!. \end{aligned}$$

The functions \(K^{f_{\mathsf {S}}}_{J, \mathbf {J}_{V_{\mathsf {S}}}}(\delta )\) define the maximal change in the value of \(f_{\mathsf {S}}\) that can occur when the input variables can vary by \(\delta \) over the intervals in \(\mathbf {J}_{V_{\mathsf {S}}}(z)\) and J. The role of \(\mathbf {J}_{V_{\mathsf {S}}}\) in the above definition is to restrict the domain of the signal variables in order to obtain the least possible relaxation bounds on the signal constraints; as was done in Definition 4 for the freeze variables.

Theorem 3

(Transference for\({\mathbb {R}}^n\)-valued Traces). Let \(\pi , \pi '\) be two \({\mathbb {R}}^n\)-valued traces such the Skorokhod distance between them is less than \(\delta \) for some finite \(\delta \). Let \(\phi \) be a closed TLTL\((\mathcal {F}_{{\mathsf {T}}}, \mathcal {F}_{{\mathsf {S}}})\) formula in negation-normal form. If \(\pi \models \phi \), then \(\pi '\models {\mathsf {rx}}_{I_{\pi , \pi '}, \mathbf {I}_{V_{\mathsf {S}}}}^{\delta }(\phi )\), where

• \({I_{\pi ,\pi '}}\) is the convex hull of \({\mathsf {tdom}}(\pi ) \cup {\mathsf {tdom}}(\pi ')\); and

• \(\mathbf {I}_{V_{\mathsf {S}}}(z)\) is the convex hull of \(\{ \pi (t)[k] \mid t\in {\mathsf {tdom}}(\pi ) \} \cup \{ \pi '(t)[k] \mid t\in {\mathsf {tdom}}(\pi ') \} \); where z is the k-th variable in the ordered set \(V_{\mathsf {S}}\). \(\quad \square \)

Theorem 3 can be strengthened similar to the strengthening mentioned for Theorem 2 by relaxing the variables over smaller intervals obtained from assignments to variables which witness \(\pi \models \phi \).

Example 5

(Spatial Constraints and Transference). Recall Example 2, suppose that the events p, q, and r are defined by the following predicates over real variables \(\alpha _1\) and \(\alpha _2\). Let \(p \equiv \alpha _1 + 10\!\cdot \!\alpha _2 \ge 3\); the predicate \(q \equiv |\alpha _1| + |\alpha _2| \le 20\); and \(r \equiv |\alpha _1| + |\alpha _2| \le 15\). Let \(\pi \) satisfy this formula with these predicates, and let \(\pi '\) be \(\delta \) close to \(\pi \), for a finite \(\delta \) under the Skorokhod metric for \({\mathbb {R}}^2\). Our robustness theorem ensures that \(\pi '\) will satisfy the relaxed formula

$$ x. \left( p^{\delta } \rightarrow \Diamond \big (y. \left( q^{\delta } \wedge \Diamond \left[ z. \left( r^{\delta } \wedge \left( (y-x)^2 + (z-y)^2 + (z-x)^2 \le d+12\!\cdot \!\delta ^2 \right) \right) \right] \right) \big )\right) \!. $$

where the relaxed predicates \(p^{\delta },q^{\delta },r^{\delta }\) are defined as follows: \(p^{\delta } \equiv \alpha _1 + 10\!\cdot \!\alpha _2 \ge 3- 22\!\cdot \!\delta \); \(q^{\delta } \equiv |\alpha _1| + |\alpha _2| \le 20+4\!\cdot \!\delta \); and \(r^{\delta } \equiv |\alpha _1| + |\alpha _2| \le 15+4\!\cdot \!\delta \). \(\quad \square \)

Open image in new window

We have implemented a streaming, sliding window-based monitoring routine which checks, given a fixed \(\delta \), whether the linear interpolations of two time-sampled traces are at Skorokhod distance at most \(\delta \) away from each other. The least such \(\delta \) value is then computed by binary search over the monitoring routine. The upper limit of the search range is set to the pointwise metric (i.e., assuming the identity retiming) between the two traces. The traces to the monitoring routine are pre-scaled, each dimension (and the time-stamp) is scaled by a different constant. The constants are chosen so that after scaling, one unit of deviation in one dimension is as undesirable as one unit of jitter in other dimensions.

We have integrated the monitoring routine in an adaptive testing procedure for Simulink blocks based on Algorithm 1. The output of Algorithm 1 is compared against tolerance levels (e.g., maximum allowed jitter) given by the engineering requirements. In the following, we evaluate the effectiveness of the Skorokhod metric in conformance testing of Simulink applications.

Recall that the Skorokhod distance computation involves a sequence of monitoring calls with different \(\delta \) values picked by a binary-search procedure. Thus, the total time to compute \({\mathscr {D}}_{\mathscr {S}}\) is the sum over the computation times for individual monitoring calls plus some bookkeeping. In Table 1, we make a distinction between the average time to monitor traces (given a \(\delta \) value), and the average time to compute \({\mathscr {D}}_{\mathscr {S}}\). There are an average of 6 monitoring calls per \({\mathscr {D}}_{\mathscr {S}}\) computation. We ran 64 simulations by choosing different input values, and then computing \({\mathscr {D}}_{\mathscr {S}}\) for increasing window sizes. As the window size increases, the average \({\mathscr {D}}_{\mathscr {S}}\) decreases and the computation time increases linearly, as expected from Theorem 1. Finally, the Skorokhod distance can be significantly smaller than the simpler metric \({\mathscr {D}}_2\) (defined as the maximum of the pointwise \(L_2\) norm). This discrepancy becomes more prominent with increased window size. With a window size of 100, the variation between \({\mathscr {D}}_{\mathscr {S}}\) and \({\mathscr {D}}_2\) was up to \(38\,\%\) (mean difference of \(10\,\%\) with std. deviation of \(11\,\%\)).

Case Study I: LQR-Based Controller. The first case study for conformance testing is an aircraft pitch control application taken from the openly accessible control tutorials for Matlab and Simulink [27]. The authors describe a linear dynamical system of the form: \(\dot{\mathbf {x}} = (A-BK)\mathbf {x} + B\theta _{des}\). Here, \(\mathbf {x}\) describes the vector of continuous state variables and \(\theta _{des}\) is the desired reference provided as an external input. One of the states in the \(\mathbf {x}\) vector is the pitch angle \(\theta \), which is also the system output. The controller gain matrix K is computed using the linear quadratic regulator method [5], a standard technique from optimal control. We are interested in studying a digital implementation of the continuous-time controller obtained using the LQR method. To do so, we consider sampled-data control where the controller samples the plant output, computes, and provides the control input to the plant every \(\varDelta \) s. To model sensor delay, we add a fixed delay element to the system; thus, the overall system now represents a delay-differential equation.

Control engineers are typically interested in the step response of a system. In particular, quantities such as the overshoot/undershoot of the output signal (maximum positive/negative deviation from a reference value) and the settling time (time it takes for transient behaviors to converge to some small region around the reference value) are of interest. Given a settling time and overshoot for the first system, we would like the second system to display similar characteristics. We remark that both of these properties can be expressed in STL (and hence in TLTL\((\mathcal {F}_{{\mathsf {T}}}, \mathcal {F}_{{\mathsf {S}}})\)), see [21] for details. We quantify system conformance (and thereby adherence to requirements) in terms of the Skorokhod distance, or, in other words, maximum permitted time/space-jitter value \(\delta \). For this system, we know that at nominal conditions, the settling time is approximately 2.5 s, and that we can tolerate an increase in settling time of about 0.5 s. Thus, we chose a time-scaling factor of \(2 = \frac{1}{0.5}\). We observe that the range of \(\theta \) is about 0.4 radians, and specify an overshoot of \(20\,\%\) of this range as being permissible. Thus, we pick a scaling factor of \(\frac{1}{0.08}\) for the signal domain. In other words, Skorokhod distance \(\delta = 1\) corresponds to either a time-jitter of 0.5 s, or a space-discrepancy of 0.08 radians.

Case Study II: Air-Fuel Ratio Controller. In [21], the authors present three systems representing an air-fuel ratio (\(\lambda \)) controller for a gasoline engine, that regulate \(\lambda \) to a given reference value of \(\lambda _{\text {ref}} = 14.7\). Of interest to us are the second and the third systems. The former has a continuous-time plant model with highly nonlinear dynamics, and a discrete-time controller model. In [22], the authors present a version of this system where the controller is also continuous. We take this to be \(\mathfrak {A}_1\). The third system in [21] is a continuous-time closed-loop system where all the system differential equations have right-hand-sides that are polynomial approximations of the nonlinear dynamics in \(\mathfrak {A}_1\). We call this polynomial dynamical system \(\mathfrak {A}_2\). The rationale for these system versions is as follows: existing formal methods tools cannot reason about highly nonlinear dynamical systems, but tools such as Flow* [10], C2E2 [16], and CORA [3] demonstrate good capabilities for polynomial dynamical systems. Thus, the hope is to analyze the simpler systems instead. In [21], the authors comment that the system transformations are not accompanied by formal guarantees. By quantifying the difference in the system behaviors, we hope to show that if the system \(\mathfrak {A}_2\) satisfies the temporal requirements \(\varphi \) presented in [21], then \(\mathfrak {A}_1\) satisfies a relaxation of \(\varphi \). We pick a scaling factor of 2 for the time domain, as a time-jitter of 0.5 s is the maximum deviation we wish to tolerate in the settling time, and pick \(0.68 = \frac{1}{0.1*\lambda _{\text {ref}}}\) as the scaling factor for \(\lambda \) (which corresponds to the worst case tolerated discrepancy in the overshoot).

Open image in new window

Case Study III: Engine Timing Model. The Simulink demo palette from Mathworks [26] contains a system representing a four-cylinder spark ignition internal combustion engine based on a model by Crossley and Cook [11]. This system is then enhanced by adding a proportional plus integral (P+I) control law. The integrator is used to adjust the steady-state throttle as the desired engine speed set-point changes, and the proportional term compensates for phase lag introduced by the integrator. In an actual implementation of such a system, such a P+I controller is implemented using a discrete-time integrator. Such integrator blocks are typically associated with a particular numerical integration technique, e.g., forward-Euler, backward-Euler, trapezoidal, etc. It is expected that different numerical techniques will produce slight variation in the results. We wish to quantify the effect of using different numerical integrators in a closed-loop setting. We checked if the user-provided tolerance of \(\delta = 1.0\) is satisfied by systems \(\mathfrak {A}_1\) and \(\mathfrak {A}_2\), where \(\mathfrak {A}_1\) is the original system provided in [26] and \(\mathfrak {A}_2\) is a modified system that uses the backward Euler method to compute the discrete-time integral in the controller. We scale the outputs in such a way that a value discrepancy of \(1\,\%\) of the the output range (\(\sim 1000\)) is equivalent to a time discrepancy of 0.1 s. These values are chosen to bias the search towards finding signals that have a small time jitter. This is an interesting scenario for this case study where the two systems are equivalent except for the underlying numerical integration solver. We find the signal shown in Fig. 2, for which we find output traces with Skorokhod distance 1.04. The experiment uses 296 simulations and the total time taken to find the counterexample is 677 s.

We argue that the Skorokhod metric provides a robust basis for checking conformance between dynamical systems. We showed that it provides transference of a rich class of temporal logic properties and that it can be computed efficiently, both in theory and in practice. Our experiments indicate that conformance checking using the Skorokhod metric can be integrated into a testing flow for Simulink models and can find non-conformant behaviors effectively, allowing for independent weighing of time and value distortions.