Statistical Verification of Hyperproperties for Cyber-Physical System

TL;DR

This paper introduces a new framework for verifying hyperproperties in cyber-physical systems using a novel temporal logic and statistical model checking, enabling analysis of complex probabilistic relationships across multiple system executions.

Contribution

It proposes HyperPSTL, a new logic for specifying hyperproperties, and develops SMC algorithms for verifying probabilistic hyperproperties in CPS modeled as PUSs.

Findings

Effective verification of hyperproperties demonstrated on CPS benchmarks.

New SMC algorithms handle nested probabilistic operators and multiple paths.

Validated on Toyota Powertrain Control System benchmark.

Abstract

Many important properties of cyber-physical systems (CPS) are defined upon the relationship between multiple executions simultaneously in continuous time. Examples include probabilistic fairness and sensitivity to modeling errors (i.e., parameters changes) for real-valued signals. These requirements can only be specified by hyperproperties. In this work, we focus on verifying probabilistic hyperproperties for CPS. To cover a wide range of modeling formalisms, we first propose a general model of probabilistic uncertain systems (PUSs) that unify commonly studied CPS models such as continuous-time Markov chains (CTMCs) and probabilistically parametrized Hybrid I/O Automata. To formally specify hyperproperties, we propose a new temporal logic, hyper probabilistic signal temporal logic (HyperPSTL) that serves as a hyper and probabilistic version of the conventional signal temporal logic (STL). Considering complexity of real-world systems that can be captured as PUSs, we adopt a statistical model checking (SMC) approach for their verification. We develop a new SMC technique based on the direct computation of the significance levels of statistical assertions for HyperPSTL specifications, which requires no a priori knowledge on the indifference margin. Then, we introduce SMC algorithms for HyperPSTL specifications on the joint probabilistic distribution of multiple paths, as well as specifications with nested probabilistic operators quantifying different paths, which cannot be handled by existing SMC algorithms. Finally, we show the effectiveness of our SMC algorithms on CPS benchmarks with varying levels of complexity, including the Toyota Powertrain Control~System.

Tables

Table 1. Table 1 . Accuracy (Acc.), average number of samples (Sam.), average execution time (Time), and SMC results (Ans.) for sensitivity ( 5 ) with parameters δ 𝛿 \delta and ε 𝜀 \varepsilon of Thermostat under significance level α 𝛼 \alpha (note the 1 − α 1 𝛼 1-\alpha accuracy guarantee by our SMC method).

$\delta$	$\epsilon$	$\alpha$	Acc.	Sam.	Time (s)	Ans.
0.9	0.05	0.05	1.00	1.8e+02	2.1e+02	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$
0.9	0.05	0.01	1.00	5.0e+02	6.0e+02	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$
0.9	0.01	0.05	1.00	2.8e+01	3.4e+01	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$
0.9	0.01	0.01	1.00	4.6e+01	5.5e+01	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$
1.1	0.05	0.05	0.99	3.0e+02	3.5e+02	$\mathrm{𝚃𝚛𝚞𝚎}$
1.1	0.05	0.01	0.99	6.1e+02	7.3e+02	$\mathrm{𝚃𝚛𝚞𝚎}$
1.1	0.01	0.05	1.00	1.3e+02	1.6e+02	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$
1.1	0.01	0.01	1.00	2.2e+02	2.6e+02	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$

Table 2. Table 2 . Accuracy (Acc.), average number of samples (Sam.), average execution time (Time), and SMC results (Ans.) for sensitivity ( 5 ) with parameters δ 𝛿 \delta and ε 𝜀 \varepsilon of Toyota Powertrain under significance level α 𝛼 \alpha (note the 1 − α 1 𝛼 1-\alpha accuracy guarantee by our SMC method).

$\delta$	$\epsilon$	$\alpha$	Acc.	Sam.	Time (s)	Ans.
0.15	0.95	0.05	1.00	5.9e+01	8.1e+00	$\mathrm{𝚃𝚛𝚞𝚎}$
0.15	0.95	0.01	1.00	9.0e+01	1.3e+01	$\mathrm{𝚃𝚛𝚞𝚎}$
0.15	0.99	0.05	0.99	6.6e+01	9.1e+00	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$
0.15	0.99	0.01	1.00	9.7e+01	1.4e+01	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$
0.20	0.95	0.05	0.98	5.9e+01	8.1e+00	$\mathrm{𝚃𝚛𝚞𝚎}$
0.20	0.95	0.01	1.00	9.0e+01	1.2e+01	$\mathrm{𝚃𝚛𝚞𝚎}$
0.20	0.99	0.05	1.00	3.0e+02	4.2e+01	$\mathrm{𝚃𝚛𝚞𝚎}$
0.20	0.99	0.01	0.99	4.6e+02	1.8e+02	$\mathrm{𝚃𝚛𝚞𝚎}$

Table 3. Table 3 . Accuracy (Acc.), average number of samples (Sam.), average execution time (Time), and SMC results (Ans.) for fairness ( 7 ) with parameters t 𝑡 t , δ 𝛿 \delta and ε 𝜀 \varepsilon of a Small Queueing Network (with 1 front-end server and 2 back-end servers) under significance level α 𝛼 \alpha .

$t$	$\delta$	$\epsilon$	$\alpha$	Acc.	Sam.	Time (s)	Ans.
0.1	0.1	0.1	0.05	1.00	1.5e+02	8.0e+01	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$
0.1	0.5	0.5	0.05	1.00	1.4e+02	8.2e+01	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$
5.0	0.1	0.1	0.05	1.00	7.3e+02	3.9e+02	$\mathrm{𝚃𝚛𝚞𝚎}$
5.0	0.5	0.5	0.05	1.00	3.1e+01	1.9e+01	$\mathrm{𝚃𝚛𝚞𝚎}$

Table 4. Table 4 . Accuracy (Acc.), average number of samples (Sam.), average execution time (Time), and SMC results (Ans.) for fairness of a Large Queueing Network (with 25 front-end servers and 20 back-end servers) with parameters defined as in Table   3 .

$t$	$\delta$	$\epsilon$	$\alpha$	Acc.	Sam.	Time (s)	Ans.
0.1	0.1	0.1	0.05	1.00	2.1e+02	3.2e+02	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$
0.1	0.5	0.5	0.05	1.00	3.3e+02	4.9e+02	$\mathrm{𝙵𝚊𝚕𝚜𝚎}$
5.0	0.1	0.1	0.05	1.00	6.8e+02	1.1e+03	$\mathrm{𝚃𝚛𝚞𝚎}$
5.0	0.5	0.5	0.05	1.00	4.2e+01	6.6e+01	$\mathrm{𝚃𝚛𝚞𝚎}$

Equations

$$\frac{\mathrm{d}\mathtt{x_{\mathcal{A}}}}{\mathrm{d}t}=\mathtt{F_{\mathcal{A}}}\big{(}\mathtt{m_{\mathcal{A}}},\mathtt{x_{\mathcal{A}}},\mathtt{u_{\mathcal{A}}}\big{)};$$

$$X(t)=X_i\text{ if }t\in [T_i,T_{i+1}),$$

$$X_{i+1}=j\text{ when }{X}_i=l\text{ if }{\sum }_{k=1}^{j-1}{M}_{lk}\le d_{1i}<{\sum }_{k=1}^j{M}_{lk},$$

$$T_{i+1}=T_{i}+\mathbf{F}^{-1}_{\mathrm{Exp}\big{(}\sum_{k\in[n],k\neq l}M_{lk}\big{)}}\big{(}d_{2i}\big{)},\ T_{0}=0,$$

$$\phi ::={\mathsf{a}}^{\pi }|{\phi }^{\pi }|\neg \phi |\phi \wedge \phi |\phi {\mathcal{U}}_{[t_1,t_2]}\phi |p\bowtie p$$

$$p::={\mathbb{P}}^{\Pi }\phi |{\mathbb{P}}^{\Pi }p|f(p,\dots ,p)$$

$$\begin{array}{rl}& \mathsf{fv}({\mathsf{a}}^{\pi })=\{\pi \},\mathsf{fv}({\phi }^{\pi })=\{\pi \},\mathsf{fv}(\neg \phi )=\mathsf{fv}(\phi ),\\ & \mathsf{fv}({\phi }_1\wedge {\phi }_2)=\mathsf{fv}({\phi }_1)\cup \mathsf{fv}({\phi }_2),\mathsf{fv}({\phi }_1{\mathcal{U}}_{[t_1,t_2]}{\phi }_2)=\mathsf{fv}({\phi }_1)\cup \mathsf{fv}({\phi }_2),\\ & \mathsf{fv}(p_1\bowtie p_2)=\mathsf{fv}(p_1)\cup \mathsf{fv}(p_2),\mathsf{fv}(f(p_1,\dots ,p_n))={\bigcup }_{i\in [n]}\mathsf{fv}(p_i),\\ & \mathsf{fv}({\mathbb{P}}^{\Pi }(\phi ))=\mathsf{fv}(\phi )\backslash \Pi ,\mathsf{fv}({\mathbb{P}}^{\Pi }(p))=\mathsf{fv}(p)\backslash \Pi .\end{array}$$

$$\begin{array}[]{l@{\hspace{1em}}c@{\hspace{1em}}l}(\mathcal{S},X)\models\varphi\hfil\hskip 10.00002pt&\Leftrightarrow\hfil\hskip 10.00002pt&\mathcal{S}\models\llbracket\varphi\rrbracket_{V_{X}},\end{array}$$

$$\begin{array}[]{l@{\hspace{1em}}c@{\hspace{1em}}l}\mathcal{S}\models\llbracket p\Join p\rrbracket_{V_{X}}\hfil\hskip 10.00002pt&\Leftrightarrow\hfil\hskip 10.00002pt&\mathcal{S}\models\llbracket p\rrbracket_{V_{X}}\Join\llbracket p\rrbracket_{V_{X}}\\ \mathcal{S}\models\llbracket f(p,\ldots,p)\rrbracket_{V_{X}}\hfil\hskip 10.00002pt&\Leftrightarrow\hfil\hskip 10.00002pt&\mathcal{S}\models f\big{(}\llbracket p\rrbracket_{V_{X}},\ldots,\llbracket p\rrbracket_{V_{X}}\big{)}\\ \mathcal{S}\models\llbracket\mathbin{\mathbb{P}}^{\Pi}(\varphi)\rrbracket_{V_{X}}\hfil\hskip 10.00002pt&\Leftrightarrow\hfil\hskip 10.00002pt&\mathcal{S}\models\mathbf{Pr}_{\boldsymbol{\sigma}\sim\mathrm{Path}^{|\Pi|}(X)}\Big{(}\big{(}\mathcal{S},V_{X}[\Pi\rightarrow\boldsymbol{\sigma}]\big{)}\models\varphi\Big{)}\\ (\mathcal{S},V_{X})\models\mathsf{a}^{\pi}\hfil\hskip 10.00002pt&\Leftrightarrow\hfil\hskip 10.00002pt&\mathsf{a}\in\mathtt{L}(V_{X}(\pi)(0))\\ (\mathcal{S},V_{X})\models\Phi^{\pi}\hfil\hskip 10.00002pt&\Leftrightarrow\hfil\hskip 10.00002pt&(\mathcal{S},V_{X}(\pi))\models\Phi\\ (\mathcal{S},V_{X})\models\neg\varphi\hfil\hskip 10.00002pt&\Leftrightarrow\hfil\hskip 10.00002pt&(\mathcal{S},V_{X})\not\models\varphi\\ (\mathcal{S},V_{X})\models\varphi_{1}\land\varphi_{2}\hfil\hskip 10.00002pt&\Leftrightarrow\hfil\hskip 10.00002pt&(\mathcal{S},V_{X})\models\varphi_{1}\textrm{ and }(\mathcal{S},V_{X})\models\varphi_{2}\\ (\mathcal{S},V_{X})\models\varphi_{1}\mathbin{\mathcal{U}}_{[t_{1},t_{2}]}\varphi_{2}\hfil\hskip 10.00002pt&\Leftrightarrow\hfil\hskip 10.00002pt&\exists t\in[t_{1},t_{2}].\ \Big{(}\forall t^{\prime}<t.\big{(}\mathcal{S},V_{X}^{(t^{\prime})}\big{)}\models\varphi_{1}\Big{)}\land\big{(}\mathcal{S},V_{X}^{(t)}\big{)}\models\varphi_{2}\end{array}$$

$$\left[\begin{array}{ccc}-1& 1& 0\\ 2& -3& 1\\ 0& 2& -2\end{array}\right]$$

$$\varphi=\mathbin{\mathbb{P}}^{\pi_{1}}\big{(}(\neg s_{1}^{\pi_{1}})\mathbin{\mathcal{U}}(s_{1}^{\pi_{1}}\mathbin{\mathcal{U}}_{[0,1]}s_{0}^{\pi_{1}})\big{)}-\mathbin{\mathbb{P}}^{\pi_{2}}\big{(}(\neg s_{2}^{\pi_{2}})\mathbin{\mathcal{U}}(s_{2}^{\pi_{2}}\mathbin{\mathcal{U}}_{[0,1]}s_{1}^{\pi_{2}})\big{)}>0.05,$$

$$\varphi=\Big{(}\frac{\mathbin{\mathbb{P}}^{\pi}\big{(}\mathrm{Init}^{\pi}\Rightarrow\Diamond(\mathsf{a}_{1}^{\pi}\land\mathsf{a}_{2}^{\pi})\big{)}}{\mathbin{\mathbb{P}}^{\pi}\big{(}\mathrm{Init}^{\pi}\Rightarrow\Diamond\mathsf{a}_{2}^{\pi}\big{)}}=\frac{1}{2}\Big{)},$$

$$\mathbf{Pr}_{\pi_{1},\pi_{2}}\big{(}|\tau^{\pi_{1}}-\tau^{\pi_{2}}|\leq\delta\big{)}>1-\varepsilon,$$

$$\mathbin{\mathbb{P}}^{\{\pi_{1},\pi_{2}\}}\Big{(}(\neg\mathtt{Q}^{\pi_{1}}\land\neg\mathtt{Q}^{\pi_{2}})\,\mathbin{\mathcal{U}}\big{(}(\mathtt{Q}^{\pi_{1}}\land\Diamond_{[0,\delta]}\mathtt{Q}^{\pi_{2}})\lor(\mathtt{Q}^{\pi_{2}}\land\Diamond_{[0,\delta]}\mathtt{Q}^{\pi_{1}})\big{)}\Big{)}\geq 1-\varepsilon.$$

$$\mathbin{\mathbb{P}}^{\{\pi,\pi^{\prime}\}}\bigg{(}\Big{(}\Box\big{(}\mathit{step}^{\pi}\Rightarrow\Box_{I}(x^{\pi}<c)\big{)}\wedge\Diamond\big{(}\mathit{step}^{\pi^{\prime}}\wedge\Diamond_{I}(x^{\pi^{\prime}}>c)\big{)}\Big{)}\Rightarrow\ \big{(}\Diamond_{I}d(y^{\pi},y^{\pi^{\prime}})>c^{\prime}\big{)}\bigg{)}>1-\varepsilon,$$

$$\mathbf{Pr}_{\pi_{1}}\Big{(}\big{|}\mathbf{Pr}_{\pi_{2}}(\tau_{i}^{\pi_{1}}-\tau_{j}^{\pi_{2}}>t)-\mathbf{Pr}_{\pi_{2}}(\tau_{j}^{\pi_{2}}-\tau_{i}^{\pi_{1}}>t)\big{|}<\delta\Big{)}>1-\varepsilon.$$

$$\begin{split}&\mathbin{\mathbb{P}}^{\pi_{1}}\Big{(}\big{|}\mathbin{\mathbb{P}}^{\pi_{2}}\big{(}(\neg\mathtt{Q}_{i}^{\pi_{1}}\land\neg\mathtt{Q}_{j}^{\pi_{2}})\mathbin{\mathcal{U}}(\mathtt{Q}_{i}^{\pi_{1}}\land\Diamond_{[\tau,\infty)}\mathtt{Q}_{j}^{\pi_{2}})\big{)}\\ &\hskip 6.0pt-\mathbin{\mathbb{P}}^{\pi_{2}}\big{(}(\neg\mathtt{Q}_{i}^{\pi_{1}}\land\neg\mathtt{Q}_{j}^{\pi_{2}})\mathbin{\mathcal{U}}(\mathtt{Q}_{j}^{\pi_{2}}\land\Diamond_{[\tau,\infty)}\mathtt{Q}_{i}^{\pi_{1}})\big{)}\big{|}\leq\delta\Big{)}\geq 1-\varepsilon.\end{split}$$

$$\quad p_{\varphi}=\mathbf{Pr}_{\boldsymbol{\sigma}\sim\mathrm{Path}^{|\Pi|}(X)}\Big{(}\big{(}\mathcal{S},V_{X}[\Pi\rightarrow\boldsymbol{\sigma}]\big{)}\models\varphi\Big{)}<p.$$

$$\phi (\mathbf{\sigma })=\{\begin{array}{ll}1,& \text{ if φ is true on }\mathbf{\sigma },\\ 0,& \text{ otherwise.}\end{array}$$

$$|p_{\phi }-p|>\delta$$

$$\begin{array}{r}{H}_0:p_{\phi }=p-\delta ,H_1:p_{\phi }=p+\delta ,\end{array}$$

$$T=\sum\nolimits_{i\in[N]}\varphi\big{(}\sigma^{(i)}_{1},\ldots,\sigma^{(i)}_{K}\big{)}$$

$$\mathcal{A}((\mathcal{S},X)\models \Phi )=\{\begin{array}{ll}1,& \text{ if }T/N<p\\ 0,& \text{ otherwise.}\end{array}$$

$$p_{\phi }\ne p,$$

$${\alpha }_{\text{CP}}(a,b|T,N)=1-\{\begin{array}{ll}(1-a{)}^N-(1-b{)}^N& \text{ if }T=0\\ b^N-a^N& \text{ if }T=N\\ F_{\text{Beta}}(b|T+1,N-T)-F_{\text{Beta}}(a|T,N-T+1)& \text{ otherwise.}\end{array}$$

$$[a,b]=\{\begin{array}{ll}[0,p],& \text{ if }T/N<p,\\ [p,1],& \text{ if }T/N>p.\end{array}$$

$${\alpha }_{\text{CP}}(a,b|T,N)\approx 1-F_{\text{Beta}}(b|T,N-T)+F_{\text{Beta}}(a|T,N-T).$$

$$(\mathcal{S},X)\models(p_{1},\ldots,p_{n})\in D\Leftrightarrow\mathcal{S}\models\big{(}\llbracket p_{1}\rrbracket_{V_{X}},\ldots,\llbracket p_{n}\rrbracket_{V_{X}}\big{)}\in D,$$

$$\begin{split}&\Phi_{1}=\big{(}f_{1}(\mathbin{\mathbb{P}}^{\Pi_{1}}\varphi_{1},f_{2}(\mathbin{\mathbb{P}}^{\Pi_{2}}\varphi_{2},\mathbin{\mathbb{P}}^{\Pi_{3}}\varphi_{3}))>c_{1}\big{)}\\ &\Phi_{2}=\big{(}f_{3}(\mathbin{\mathbb{P}}^{\Pi_{2}}\varphi_{2}))<c_{2}\big{)}\end{split}$$

$$D=\big{\{}(x_{1},x_{2},x_{3})\in[0,1]^{3}\,|\,f_{1}(x_{1},f_{2}(x_{2},x_{3}))>c_{1},f_{3}(x_{2})<c_{2}\big{\}}.$$

$$(\mathcal{S},X)\models\big{(}\allowbreak\mathbin{\mathbb{P}}^{\Pi_{1}}\varphi_{1},\allowbreak\ldots,\allowbreak\mathbin{\mathbb{P}}^{\Pi_{n}}\varphi_{n}\big{)}\in D,$$

Full text

Cyber-physical systems, hyperproperties, statistical model checking, embedded control software.

††copyright: acmcopyright††journalyear: 2019††journal: TECS††journalvolume: 1††journalnumber: 1††article: 1††publicationmonth: 1††ccs: Computer systems organization Embedded and cyber-physical systems††ccs: Theory of computation Logic and verification††ccs: Security and privacy Formal security models††ccs: Mathematics of computing Hypothesis testing and confidence interval computation††ccs: Software and its engineering Software verification and validation

1. Introduction

Ensuring safety of controllers in embedded and cyber-physical systems (CPS) using closed-loop system verification is a challenging problem, due to the inherent uncertainties in system dynamics and the environment. In systems where the uncertainties can be captured in a probabilistic manner, two prominent verification approaches are exhaustive (Baier and Katoen, 2008) and statistical (Legay and Viswanathan, 2015). The exhaustive approach, with full knowledge of a system model, computes the satisfying probability of the desired properties arithmetically; on the other hand, the statistical approach estimates the probabilities from sampling, and makes assertions with a certain significance level (an upper bound of the probability of returning a wrong answer). Accordingly, statistical model checking (SMC) is more capable of handling “black-box”, high-dimension or large-scale system models.

Properties of interest for such systems, from ‘simple’ probabilistic safety or reachability, to the ones that capture dynamical responses under complex conditions, are usually formally defined by probabilistic temporal logic specifications. Conventional probabilistic temporal logics, such as the probabilistic computational tree logic (PCTL) (Hansson and Jonsson, 1994), as well as its extension PCTL∗ (Baier and Katoen, 2008), can only specify probabilistic properties without explicitly quantifying over different paths of the system; that is, they cannot simultaneously and explicitly quantify over multiple distinctive paths to fully express their inter-relations. This prevents them from capturing important safety/performance hyperproperties (Clarkson and Schneider, 2008; Ábrahám and Bonakdarpour, 2018) that involve multiple execution paths, such as sensitivity to model errors, detectability of system anomalies, and fairness when more than once process/client are controlled/serviced.

For example, consider embedded controllers such as the Toyota Powertrain control system benchmark (Jin et al., 2014b), for which both exhaustive (e.g., (Duggirala et al., 2015)) and statistical (e.g., (Roohi et al., 2017)) verification techniques have been introduced. However, all these techniques are restricted to the use on a dynamical system model obtained for fixed system parameters. In general, such parameters are experimentally derived and thus, should be considered as random variables with unknown probability distributions, instead of the fixed values. In addition, some of the system parameters might change to a degree, due to system ‘wear-and-tear’. Hence, it is critical to enable analysis of system sensitivity to model errors, by providing a formal logic to capture such properties, as well as methods to verify how system evolution (execution) changes for different parameters of the system model.

Specifying these properties, such as sensitivity to change of parameters, involves probabilistic quantification over multiple paths (a path and its deviation), and thus can only be captured as hyperproperties (Clarkson and Schneider, 2008). For example, as illustrated in Figure 1, for sensitivity analysis we can check whether the deviation $\pi_{2}$ of a path $\pi_{1}$ under probabilistic uncertainty stays close probabilistically, such that there is limited variation in the hitting time $\tau$ to a desired working region. Although the sensitivity may be expressed as a non-hyperproperty if ‘expected’ hitting times are known in advance (as a reference) for any entry into the desired operating region, such information is usually unavailable for complex systems.

Consequently, in this work, we first introduce a probabilistic temporal logic for hyperproperties expressed on real-valued continuous-time signals, referred to as Hyper Probabilistic Signal Temporal Logic (HyperPSTL). HyperPSTL can be viewed as a hyper extension and generalization of the probabilistic signal/metric-interval temporal logic (Sadigh and Kapoor, 2016; Wang et al., 2016), a probabilistic version of HyperSTL (Nguyen et al., 2017), and a continuous-time extension and generalization of HyperPCTL (Ábrahám and Bonakdarpour, 2018). HyperPSTL extends those logics by enabling (1) reasoning about the probability of paths by adding a probability operator, and (2) reasoning about multiple paths simultaneously, i.e., hyperproperties specifying the relationship between different paths by associating path variables to the atomic propositions.

To allow us to cover a range of modeling formalisms, we introduce a very general system model – probabilistic uncertain systems (PUS) – and define the semantics of HyperPSTL on it. Generally, they are ‘black-box’ probabilistic dynamical systems with unknown dynamics on a given state space. A PUS incorporates nondeterminism as its input and probabilism as its parameters, both of which are time functions of values of general types, including real, integer or categorical/Boolean. Given the values of the input and the parameters, we can draw a time-dependent sample path from the PUS, which can also be of general types. We note that this general model captures commonly studied models such as CTMCs and hybrid I/O automata with probabilistic parameters – referred to as probabilistically parameterized hybrid I/O automata (P2HIOA).

We define the semantics of HyperPSTL on a set of paths of the PUS through pre-defined labels or predicates on the state spaces. We show by concrete examples that this allows for the capturing of properties such as (i) anomaly detectability, or (ii) sensitivity to model errors due to probabilistic uncertainty of parameters of hybrid I/O automata, and (iii) workload fairness in queueing networks modeled as continuous-time Markov chains (CTMC). To verify HyperPSTL specifications on a PUS for a given input, we introduce new SMC algorithms using Clopper-Pearson (CP) bounds (Clopper and Pearson, 1934).

Unlike most previous SMC methods based on the sequential probability ratio test (SPRT) (Legay and Viswanathan, 2015; Larsen and Legay, 2016), this approach requires no a priori knowledge on the indifference margin. Our SMC algorithm can verify a HyperPSTL specification to arbitrarily small (non-zero) significance levels, which is also different from (Barbot et al., 2017) on using CP bounds to estimate satisfying probability for given samples. The conservativeness of the CP significance level ensures that the desired significance level is strictly achieved even in the worst case. Specifically, we iteratively draw samples from the objective model and compute the significance level derived from the CP confidence interval. The algorithm stops if the CP significance level is smaller than the desired significance levels. Further improvement in sample efficiency is also possible using the upper bounds derived from Wilson scores, Jeffreys interval or Agresti-Coull interval, at the cost of yielding only asymptotic correctness (see (Brown et al., 2001) for a review of these confidence intervals).

To the best of our knowledge, this work is the first to enable SMC on probabilistic temporal logic for hyperproperties in continuous time. Compared to common non-hyper temporal logics, HyperPSTL allows uniquely for (1) probability quantification over multiple paths, (2) defining specifications on their joint probabilistic distribution, such as comparison of probabilities, and (3) nesting of probabilistic operators that quantify over different sets of paths. We address these challenges by deriving SMC methods capable of handling these cases, while providing provable significance levels. We also show that HyperPSTL subsumes several existing signal temporal logics for probabilistic properties (Roohi et al., 2017; Wang et al., 2016) (see Section 3), and thus, the SMC algorithms for HyperPSTL introduced in this work, also apply to them.

Finally, we show the effectiveness of our SMC methods on embedded case studies with different complexity and modeling formalism. Specifically, we statistically verified sensitivity of a thermostat and the Toyota Powertrain System with uncertain parameters, as well as fairness in queueing networks of different sizes; we generated probabilistic guarantees with high significance levels. For example, for the Toyota Powertrain and the large queueing networks, where exhaustive verification is not possible, we construct upper bounds on the $0.95$ and $0.99$ percentiles of their sensitivity and fairness, respectively, with a confidence level of $0.99$ using only a few hundred samples. This achieves, for the first time, statistical verification of probabilistic hyperproperties on continuous signals, which can be used even for real-world sized CPS.

Organization

We introduce probabilistic uncertain systems in Section 2, followed by the syntax and semantics of HyperPSTL in Section 3. In Section 4, we show how HyperPSTL can be used to capture desired properties of CPS. Our SMC algorithms are presented in Section 5 with an emphasis on the use of CP significance levels and the handling of probabilities involving multiple paths. We apply the SMC techniques to three embedded case studies (e.g., Toyota Powertrain System) in Section 6, before concluding in Section 7.

Notation

We denote the sets of integers, rational, real, and non-negative real numbers by $\mathbb{N}$, $\mathbb{Q}$, , and ≥0, respectively. The domain and image of a function is denoted by $\mathrm{Dom}(\cdot)$ and $\mathrm{Image}(\cdot)$, respectively. Let $\mathbb{N}_{\infty}=\mathbb{N}\cup\{\infty\}$ and $\mathbb{Q}_{\infty}=\mathbb{Q}\cup\{\infty\}$. For $n\in\mathbb{N}$, let $[n]=\{1,\ldots,n\}$. The indicator function is denoted by $\mathbf{I}$. We denote the (Borel) measure of a measurable set by $\mu_{\mathrm{Borel}}(\cdot)$. The cardinality and the power set of a set are denoted by $|\cdot|$ and $2^{\cdot}$. For any set $I\subseteq^{n}$, we denote its boundary, interior, and closure by $\partial I$, $I^{\circ}$, and $\bar{I}$, respectively. The empty set is denoted by $\emptyset$. With a slight abuse of notation, given a map $V:A\rightarrow B$ and $A^{\prime}\subseteq A$, let $V(A^{\prime})={\bigcup}_{a\in A^{\prime}}V(a)$.

We refer to a function of time $\sigma:_{\geq 0}\rightarrow^{n}$ as a signal, and denote by $\sigma^{(t_{1})}$ its $t_{1}$-time shift defined as $\sigma^{(t_{1})}(t_{2})=\sigma(t_{1}+t_{2})$ for $t_{1}\in_{\geq 0}$. We denote the binomial distribution by $\mathrm{Binom}(n,p)$, the exponential distribution parametrized by rate by $\mathrm{Exp}(\lambda)$, and the beta distributions with shape parameters by $\mathrm{Beta}(\alpha,\beta)$. A random variable $X$ drawn from probability distribution $\mu$ is denoted by $X\sim\mu$.

2. Probabilistic Uncertain Systems

To allow for defining HyperPSTL for a large class of commonly used modeling formalisms in a way that facilitates design of SMC techniques, we introduce a very general system model, which we refer to as probabilistic uncertain systems (PUS) (Figure 2). Such model can be viewed as a “black-box” probabilistic dynamical system with an explicit state space, where the randomness only comes from the time-dependent parameters drawn from random processes. The PUS generalizes models such as CTMCs and P2HIOA, and provides a unified framework for the SMC of probabilistic hyperproperties.

A PUS is a ’black-box’ probabilistic dynamical system on a given state space, which we denote by $\mathcal{X}$. Its probabilistic uncertainty comes from a set of $n$ time-dependent parameters D(t)=\big{(}d_{1}(t),...,d_{n}(t)\big{)} where $t\in_{\geq 0}$, drawn from some probability distribution $\mu(\mathcal{D})$ on its domain $\mathcal{D}$. That is, the parameters of the system are drawn from an $n$-dimensional random process. The input to the system I(t)=\big{(}i_{1}(t),\dots,i_{m}(t)\big{)}\in\mathcal{I} is an $m$-dimensional function of time $t$. Given the input $I(t)$, the values of the parameters $D(t)$, and an initial state X^{\mathrm{init}}=\big{(}x^{\mathrm{init}}_{1},\dots,x^{\mathrm{init}}_{l}\big{)}\in\mathcal{X}, the system deterministically generates a time-dependent path $X:_{\geq 0}\rightarrow\mathcal{X}$ with X(t)=\big{(}x_{1}(t),\dots,x_{l}(t)\big{)}. That is, the randomness in system evolution only comes from system parameters.

The values of the system inputs, parameters, and states can be of mixed types: real, discrete or categorical, depending on the system formulation. Here, we make no assumption on the system dynamics (Markovian, causal, etc); rather, the system should be viewed just as a general deterministic map from the functions $I(t)$ and $D(t)$ to the function $X(t)$. Without knowing the value of the system parameters and given the initial state, the map from the input $I(t)$ to the path $X(t)$ is only probabilistic, i.e., the probability distribution for all uncertain variables are given.

Given the finite set of atomic propositions $\mathsf{AP}$, the labeling function on the state space $\mathtt{L}:\mathcal{X}\rightarrow 2^{\mathsf{AP}}$ defines for each state, a set of atomic propositions that hold. Alternatively, the labeling of each atomic proposition in $\mathsf{AP}$ can be represented by a predicate/Boolean function that indicates the subset of states $\mathcal{X}$ that are labeled by that atomic proposition. Then, a path of the system induces a signal $\sigma(t)=\mathtt{L}(X(t)):_{\geq 0}\rightarrow 2^{\mathsf{AP}}$. This signal indicates the set of properties that is satisfied by the system $\mathcal{S}$ at each time $t\in_{\geq 0}$ during system evolution. In Section 3, we will define hyper temporal properties on sets of such signals.

A PUS $\mathcal{S}$ is both probabilistic in the value of its parameters and nondeterministic in its input. We note that the distinction between the parameters and the input of the system is only mathematical, not physical. In practice, we can model the probabilistic input of a real system as a parameter of its PUS model, and a nondeterministic parameter of the system as an input of the PUS. In this work, we focus on SMC related to the probabilistic nature of a PUS; i.e., we consider the cases, where the input of the PUS is given such as when it is determined by the actions of the controller used to control the system (PUS). Before we introduce HyperPSTL and SMC methods for PUS, in the rest of this section, we show how our abstract notion of PUS captures two prominent computing models: continuous-time Markov chains (CTMCs) that mathematically model queuing networks and probabilistic hybrid I/O automata that model the Toyota Powertrain system mentioned in Section 1.

2.1. Probabilistically Parameterized Hybrid I/O Automata

Probabilistically Parameterized Hybrid I/O Automata (P2HIOA) are extensions of Hybrid I/O Automata (HIOA) (Lynch et al., 2003; Henzinger, 2000; Jin et al., 2014b) with probabilistic system parameters. Given the variety of mathematical models for HIOA, we build the one used in (Jin et al., 2014b) to describe the dynamics of hybrid systems with fully observable states. Specifically, an HIOA is a tuple $\mathcal{A}=(\mathtt{M_{\mathcal{A}}},\mathtt{X_{\mathcal{A}}},\mathrm{I_{\mathcal{A}}},\mathtt{U_{\mathcal{A}}},\allowbreak\mathtt{Inv_{\mathcal{A}}},\mathtt{T_{\mathcal{A}}},\mathtt{Q})$ where

• •

$\mathtt{M_{\mathcal{A}}}$ is a finite set of modes;

• •

$\mathtt{X_{\mathcal{A}}}$ is a set of $n$ state variables of real values, i.e., $\mathtt{X_{\mathcal{A}}}\subseteq^{n}$;

• •

$\mathtt{U_{\mathcal{A}}}$ is a set of $m$ input variables; their valuations can be of different types, such as $\mathbb{N}$, , and Boolean;

• •

$\mathtt{F_{\mathcal{A}}}:\mathtt{M_{\mathcal{A}}}\times\mathtt{X_{\mathcal{A}}}\times\mathtt{U_{\mathcal{A}}}\rightarrow^{n}$ defines a deterministic flow capturing system evolution in each mode, i.e., a differential equation:

[TABLE]

• •

For any mode $\mathtt{m_{\mathcal{A}}}\in\mathtt{M_{\mathcal{A}}}$, $\mathtt{Inv_{\mathcal{A}}}(\mathtt{m_{\mathcal{A}}})=\mathrm{Dom}(\mathtt{F_{\mathcal{A}}}(\mathtt{m_{\mathcal{A}}},\cdot,\cdot))$ defines the invariant of the mode;

• •

$\mathtt{T_{\mathcal{A}}}:\mathtt{M_{\mathcal{A}}}\times\mathtt{X_{\mathcal{A}}}\times\mathtt{U_{\mathcal{A}}}\rightarrow\mathtt{M_{\mathcal{A}}}\times\mathtt{X_{\mathcal{A}}}$ defines deterministic jumps triggered by $\mathtt{X_{\mathcal{A}}}\in\partial\mathtt{Inv_{\mathcal{A}}}(\mathtt{m_{\mathcal{A}}})$, i.e., when the flow hits the boundary of the invariant.

• •

$\mathrm{I_{\mathcal{A}}}\in\mathtt{M_{\mathcal{A}}}\times\mathtt{X_{\mathcal{A}}}$ is the initial condition;

• •

$\mathtt{Q}\subseteq\mathtt{U_{\mathcal{A}}}\times\mathtt{X_{\mathcal{A}}}\rightarrow\{0,1\}$ is a finite set of predicates.

The P2HIOA extend Hybrid I/O Automata by allowing system parameters, used to capture $\mathtt{F_{\mathcal{A}}},\mathtt{T_{\mathcal{A}}},\mathtt{Inv_{\mathcal{A}}}$, to be probabilistic instead of fixed. This differs from the Probabilistic Hybrid Automata introduced in (Sproston, 2000; Zhang et al., 2010), where the randomness comes only from the probabilistic jumps. In standard hybrid models, the deterministic and nonlinear flow, invariant, and jump functions $\mathtt{F_{\mathcal{A}}},\mathtt{T_{\mathcal{A}}},\mathtt{Inv_{\mathcal{A}}}$, which define its dynamics, are parameterized by quantities that are estimated from physical experiments. For example, in the Toyota Powertrain model, the mass flow rate of intake air is determined by the RPM of the engine and the pressure of the intake manifold through a polynomial, whose five parameters are fit from experimental data. Due to experimental errors, these parameters are better represented as random variables with unknown probability distributions (e.g., Gaussian or uniform with means and variances inferred from experimental data) than real numbers as in (Jin et al., 2014b). Thus, we denote the parameters for $\mathtt{F_{\mathcal{A}}}$, $\mathtt{T_{\mathcal{A}}}$, and $\mathtt{Inv_{\mathcal{A}}}$ by $D=(d_{\mathtt{F_{\mathcal{A}}}},d_{\mathtt{T_{\mathcal{A}}}},d_{\mathtt{Inv_{\mathcal{A}}}})$. To simplify our presentation, in this work, we assume that these parameters for P2HIOA are time-invariant.

Consequently, a P2HIOA can be represented by a PUS by treating (i) $I=\mathtt{U_{\mathcal{A}}}$ as the input, (ii) $X=\mathtt{M_{\mathcal{A}}}\times\mathtt{X_{\mathcal{A}}}$ as the state (the input variables $\mathtt{U_{\mathcal{A}}}$ are encoded as part of the state), (iii) $X^{\mathrm{init}}=\mathrm{I_{\mathcal{A}}}$ as the initial state, (iv) $D=(d_{\mathtt{F_{\mathcal{A}}}},d_{\mathtt{T_{\mathcal{A}}}},d_{\mathtt{Inv_{\mathcal{A}}}})$ as the parameters, and (v) the predicate $\mathtt{Q}$ as the labeling function. The probabilistic hyperproperties related to P2HIOA, as the one discussed in Sections 1 and 4, will be statistically verified on this PUS.

2.2. Continuous-Time Markov Chains

Another example of PUS are CTMCs, which are commonly used to model queuing and task scheduling in embedded computing and communication systems with uncertainties. Consider a CTMC with

• •

the states $[n]$,

• •

the initial state $X_{0}\in[n]$,

• •

the labeling function on the states $\mathtt{L}:[n]\rightarrow 2^{\mathsf{AP}}$ for a given set of labels $\mathsf{AP}$,

• •

the probability transition rate matrix $M\in^{n\times n}$, such that $\sum_{j\in[n]}\allowbreak M_{ij}=0$, where $M_{ij}$ is the transition rate from a state $i$ to a state $j$.

The CTMC can be represented by a PUS with (i) the states $\mathcal{X}=[n]$, (ii) the initial state $X_{0}$, (iii) the empty inputs (i.e., the PUS has no input), (iv) the labels $\mathsf{AP}$, and (v) the labeling function $\mathtt{L}:[n]\rightarrow 2^{\mathsf{AP}}$. For $i\in\mathbb{N}$, we draw $d_{1i}$ and $d_{2i}$ from $[0,1]$ uniformly and independently. Using Gillespie algorithm (Gillespie, 1976), with the convention $\sum_{k=1}^{0}\cdot=0$, the state $X$ of the PUS representing the CTMC evolves by:

[TABLE]

where the state jumps are determined by $d_{1i}$ as

[TABLE]

and the time lapses are determined by $d_{2i}$ using

[TABLE]

with $\mathrm{Exp}(\cdot)$ denoting the exponential distribution parameterized by rate and $\mathbf{F}^{-1}_{\cdot}$ the inverse function of the cumulative distribution function. Treating the two sequences $\{d_{1i}\}_{i\in\mathbb{N}}$ and $\{d_{2i}\}_{i\in\mathbb{N}}$ as the two parameters of the PUS, provides a PUS presentation of the CTMC. Although they are discrete sequences, we can easily represent them as continuous functions of time to fit into the definition of the parameters of the PUS. Then, the probabilistic hyperproperties (e.g., fairness) for the initial CTMC, should be statistically verified on the aforementioned PUS.

3. Hyper Probabilistic Signal-Temporal Logic

To formally express and reason about probabilistic hyperproperties on real-valued signals, we introduce the logic HyperPSTL, which can be viewed as a probabilistic extension of the signal temporal logic (STL) for hyperproperties. We introduce the syntax and semantics of HyperPSTL for PUS in Sections 3.1 and 3.2, before presenting its use to capture relevant properties of embedded systems in Section 4.

3.1. Syntax

We define HyperPSTL formulas inductively as:

[TABLE]

where

• •

$\mathsf{a}\in\mathsf{AP}$, and $\mathsf{AP}$ is the finite set of atomic propositions,

• •

$t_{1}<t_{2}$ with $t_{1},t_{2}\in\mathbb{Q}_{\infty}$,

• •

$\pi$ is a path variable, and $\Pi$ is a set of path variables,

• •

$\mathbin{\mathbb{P}}$ is the probability operator,

• •

$\Join\ \in\{<,>,=,\leq,\geq\}$,

• •

$f:^{n}\rightarrow$ is a $n$-ary measurable function, constants are viewed as [math]-ary functions,

• •

$\mathsf{fv}(\cdot)$ denotes the set of free path variables in $\varphi$, i.e., the path variables not quantified by a probability operator through (2), and $\mathsf{fv}(\varphi)=\emptyset$ in the $2^{\mathrm{nd}}$ rule of (1). This is recursively defined by:

[TABLE]

Other common logic operators can be derived as follows: $\varphi\lor\varphi^{\prime}\equiv\neg(\neg\varphi\land\neg\varphi^{\prime})$, $\mathtt{True}\equiv\varphi\vee\neg\varphi$, $\varphi\Rightarrow\varphi^{\prime}\equiv\neg\varphi\lor\varphi^{\prime}$, $\Diamond_{[t_{1},t_{2}]}\varphi\equiv\mathtt{True}\,\mathbin{\mathcal{U}}_{[t_{1},t_{2}]}\,\varphi$, and $\Box_{[t_{1},t_{2}]}\varphi\equiv\neg\Diamond_{[t_{1},t_{2}]}\neg\varphi$. In addition, we denote $\mathbin{\mathcal{U}}_{[0,\infty)}$, $\Diamond_{[0,\infty)}$, and $\Box_{[0,\infty)}$ by $\mathbin{\mathcal{U}}$, $\Diamond$, $\Box$, respectively.

A formula with $\mathsf{fv}(\varphi)=\emptyset$ is called a state formula and requires no instantiation of free path variables, thus can be evaluated on a state of the PUS. Therefore, we can associate another path variable $\pi$ to it. All other formulas are referred to as path formulas, as their correctness depends on the instantiation of their free path variables. It specifies on which path a state formula should be satisfied, so that a HyperPSTL formula can reason simultaneously on multiple paths. Finally, for the first two rules of (2), we assume that $\Pi$ is contained in $\mathsf{fv}(\varphi)$ or $\mathsf{fv}(p)$ for the probability quantification to be non-trivial. Observe that HyperPSTL can be viewed as the probabilistic version of HyperSTL (Nguyen et al., 2017) by replacing the existential and universal quantifiers over signals with probabilistic quantifiers over paths in (2).

HyperPSTL has the following unique features. It allows for the simultaneous probability quantification over several paths (as we show for the sensitivity analysis of powertrain controllers in Section 4). It also allows for the arithmetics and comparison of probabilities, and the nesting of probability operators quantifying different paths (as shown in the queueing fairness analysis described in Section 4).

HyperPSTL reduces to a non-hyper probabilistic signal temporal logic (PSTL) if it only has one path variable in it. PSTL can still define probability satisfaction of single atomic propositions, so it subsumes the MITL for probability distributions from (Wang et al., 2016). But, PSTL (and thus HyperPSTL) does not subsume PrSTL in (Sadigh and Kapoor, 2016), since it does not allow a time-varying probability threshold, as is allowed in PrSTL. Still, augmenting HyperPSTL syntax to allow time-varying functions is straight-forward.

Finally, note that to simplify our presentation of HyperPSTL syntax and semantics, while allowing for verification of complex systems such as P2HIOA, we only include simultaneous or consecutive probabilistic quantification (e.g., $\mathbin{\mathbb{P}}^{\{\pi_{1},\pi_{2}\}}$ or $\mathbin{\mathbb{P}}^{\pi_{1}}\mathbin{\mathbb{P}}^{\pi_{2}}$) over the paths from a single initial state. HyperPSTL can be augmented by allowing nested existential and universal quantification over multiple states in the same way as (Wang et al., 2019; Ábrahám and Bonakdarpour, 2018). Specifically, in addition to the probabilistic quantification over the paths, one can add extra state quantification of these paths to specify from which state the path starts, like $\exists X_{1}^{\pi_{1}}.\forall X_{2}^{\pi_{2}}.\mathbin{\mathbb{P}}^{\{\pi_{1},\pi_{2}\}}$. However, verifying such formulas generally requires exhaustive iteration over all the states, which is challenging, if not impossible, on systems with infinite state spaces like P2HIOA. Therefore, in this paper, we do not include state quantification which in our logic can be done as presented in (Wang et al., 2019; Ábrahám and Bonakdarpour, 2018).

3.2. Semantics

We define the satisfaction relation for HyperPSTL state formulas on a PUS $\mathcal{S}$ by

[TABLE]

where $V_{X}$ is an assignment of path variables to the paths of the PUS $\mathcal{S}$ starting from a state $X$, and $\llbracket\varphi\rrbracket_{V_{X}}$ is the instantiation of the assignment $V_{X}$ on $\varphi$. Here, the instantiation specifies the initial state for the formula $\varphi$.

The satisfaction relation for the HyperPSTL path formulas is defined with respect to the assignment $V_{X}$ by:

[TABLE]

where

• •

$\mathtt{L}(X)$ is the set of labels of the PUS state $X$,

• •

$\mathrm{Path}^{|\Pi|}(X)$ is the collection of all $|\Pi|$-tuples of paths starting from state $X$,

• •

$V_{X}[\Pi\rightarrow\boldsymbol{\sigma}]$ is a revision of the assignment $V_{X}$ by assigning the set $\boldsymbol{\sigma}$ of paths to the set $\Pi$ of path variables, respectively,

• •

$V_{X}^{(t)}$ is the $t$-shift of the assignment $V_{X}$, defined by \big{(}V_{X}^{(t)}(\pi)\big{)}\allowbreak=(V_{X}(\pi))^{(t)} for all path variables $\pi$ in assignment $V$.

Finally, we note the equivalence $\mathsf{a}\in\mathtt{L}(V_{X}(\pi)(0))\ \Leftrightarrow\ \mathsf{a}\in\mathtt{L}(X)$ and $(\mathcal{S},V_{X}(\pi))\models\Phi\ \Leftrightarrow\ (\mathcal{S},X)\models\Phi$ for the $4^{\mathrm{th}}$ and $5^{\mathrm{th}}$ rules in (4). The following example illustrates the semantics of HyperPSTL.

Example 0.

As shown in Figure 3, consider a CTMC $\mathcal{S}$ that models a queue that is initially empty and is of buffer size $2$. The transition rate matrix of $\mathcal{S}$ is the following:

[TABLE]

The CTMC satisfies the following HyperPSTL formula:

[TABLE]

The formula asserts that the probability difference between (i) finishing a task within $1$ time delay after first having $1$ in queue and (ii) finishing a task within $1$ time delay after first having $2$ in queue, is greater than $0.05$. The is because the probability for (i) is $\frac{2}{3}(1-e^{-3})\approx 0.633$ and that for (ii) is $1-e^{-2}\approx 0.865$. $\lhd$

3.3. Expressivity of HyperPSTL

Theorem 3.2.

HyperPSTL subsumes PSTL on CTMCs.

Proof.

It suffices to show that there exist formulas in HyperPSTL that cannot be expressed in PSTL. The general idea is to show that PSTL cannot express conditional probabilities, while HyperPSTL can. Consider the CTMC $\mathcal{S}$ given in Figure 4, and the HyperPSTL state formula

[TABLE]

expressing a conditional probability for a path $\pi$ of $\mathcal{S}$. Clearly, it is satisfied for the state $\mathrm{Init}$.

We claim that $\varphi$ cannot be expressed in PSTL. By the syntax and semantics of PSTL (Sadigh and Kapoor, 2016), it suffices to show that $\varphi$ cannot be expressed by a formula $\mathbin{\mathbb{P}}(\psi)$, where $\psi$ is a PSTL path formula derived by concatenating a set of PSTL state formulas $\varphi_{1},\ldots,\varphi_{n}$ with $\land,\neg$, or temporal operators. These state formulas are either $\mathtt{True}$ or $\mathtt{False}$ on the states $X_{0}$, $X_{1}$, $X_{2}$, and $X_{3}$. Thus, whether a path satisfies $\psi$, defines a subset of the paths from the state $X_{0}$ in the CTMC, Since the probability of a path ended up in any $X_{i}$ is $1/3$ for $i=1,2,3$, the formula $\mathbin{\mathbb{P}}(\psi)$ can only take values in $\{0,1/3,2/3,1\}$. However, by the semantics of HyperPSTL, the fractional probability on the right side of the equation has value $1/2$, thus $\varphi$ evaluates to true and cannot be expressed by $\mathbin{\mathbb{P}}(\psi)$ in PSTL. ∎

4. HyperPSTL in Action

In this section, we demonstrate how HyperPSTL can be used to capture relevant properties of CPS.

4.1. Sensitivity to Modeling Errors in P2HIOA

A typical example of a CPS that can be modeled as a P2HIOA is the automotive powertrain, where the response to the change in driving behaviors is of key interest. Note that the dynamical system parameters, even for the same type of powertrains, vary across different systems. Thus, it is critical to analyze if in most cases, the change in dynamical response of the controlled system stays within permitted amount $\delta$ when the system parameters change; i.e., if the dynamical response deviation is within $\delta$ with probability of at least $1-\varepsilon$.

Consider as an example, the sensitivity of the Toyota Powertrain Controller (Jin et al., 2014b) under probabilistic uncertainty in its dynamical parameters (i.e., the system model). As shown in Figure 1, we consider statistically verifying the probabilistic boundedness of the sensitivity of the (first) hitting time $\tau$ to a desired working region (where the error to the desired AF ratio is less than $5\%$) under the probabilistic uncertainty in RPM. Mathematically, this can be represented by

[TABLE]

for some given values $\delta,\varepsilon>0$, where $\pi_{1}$ and $\pi_{2}$ are two statistically independent sample paths of the system. To express this formally using HyperPSTL, we introduce a predicate $\mathtt{Q}$ for the desired working region of the system (i.e., under the dashed line in Figure 1). If $0\leq\tau^{\pi_{2}}-\tau^{\pi_{1}}\leq\delta$, the switch from $\neg\mathtt{Q}^{\pi_{2}}$ to $\mathtt{Q}^{\pi_{2}}$ for the path $\pi_{2}$ happens within time $\delta$ after $\neg\mathtt{Q}^{\pi_{1}}$ changes to $\mathtt{Q}^{\pi_{1}}$ for the path $\pi_{1}$.111Meanwhile, $\mathtt{Q}^{\pi_{2}}$ may switch back to $\neg\mathtt{Q}^{\pi_{2}}$ for $\pi_{2}$, but the first hitting time $\tau^{\pi_{2}}$ will not change. This can be equivalently expressed as $(\neg\mathtt{Q}^{\pi_{1}}\land\neg\mathtt{Q}^{\pi_{2}})\mathbin{\mathcal{U}}(\mathtt{Q}^{\pi_{1}}\land\Diamond_{[0,\delta]}\mathtt{Q}^{\pi_{2}})$, and accordingly, the probabilistic boundedness of sensitivity is expressed as

[TABLE]

Note that (5) involves probability quantification over a set of paths, which cannot be expressed in non-hyper temporal logics.

4.2. Probabilistic Anomaly Detectability

An important feature of CPS is detectability of system anomalies, independently of the type of used sound detector; this can be captured as probabilistic overshoot observability on system outputs, where the input overshoot captures that an anomaly has occurred. Specifically, we require that with probability of at least $1-\varepsilon$ it holds that: if (i) in one execution, a signal $\pi$ steps (i.e., anomaly starts) and then stays bounded (e.g., within the modeled noise bound) for some time interval $I$; and (ii) in another execution, signal $\pi^{\prime}$ steps and then overshoots (i.e., beyond the noise bound); then (iii) the distance between the two signals is greater than a predefined threshold (i.e., the anomaly overshoot can be observed by a detector on system output). This is captured as the following HyperPSTL formula

[TABLE]

where $x$ is the input and $y$ is the output. As (5), the formula in (6) also involves probability quantification over a set of paths.

4.3. Workload Fairness in Queueing Networks

As shown in Figure 5, consider an embedded processing system with $n$ front servers and $m$ back servers (as in e.g., (Tang and Bergmann, 2015)). The requests (e.g., task, packets) arrive at each front-end queue, probabilistically over time, into buffers of different sizes. For each queue, the requests are preprocessed with probabilistic execution times, and delivered to back servers with different buffer sizes, following some scheduling policy. We can probabilistically model the arrival and processing of requests by Markov Modulated Poisson Processes (MMPP) of different parameters across all the servers (Bolch et al., 2006). In the general case, this setup yields no easy exhaustive solution.

Our goal is to check if a request-delivering policy between the front and back servers is fair (Eryilmaz and Srikant, 2006; Georgiadis et al., 2006); i.e., if a back server $i$ is more likely to be overloaded than another back server $j$. Let $\tau_{i}$ and $\tau_{j}$ be the overload times of the back server $i$ and $j$, respectively. We define a fairness property that with probability of at least $1-\varepsilon$, given the overload time for the back server $i$, the back server $j$ is overloaded much earlier or later (more than some $t>0$ than that time), with approximately equal probability (i.e., difference less than some $\delta$)

[TABLE]

To express this property in HyperPSTL, for $i\in[1,m]$, let $\mathtt{Q}_{i}$ be the predicate of the overload of the buffer of the $i^{th}$ back server. The event that the back server $i$ is overloaded earlier than the back server $j$ more than time $\tau>0$, can be expressed as $(\neg\mathtt{Q}_{i}^{\pi_{1}}\land\neg\mathtt{Q}_{j}^{\pi_{2}})\mathbin{\mathcal{U}}(\mathtt{Q}_{i}^{\pi_{1}}\land\Diamond_{[\tau,\infty)}\mathtt{Q}_{j}^{\pi_{2}})$. Hence, a fairness policy requirement is captured by a HyperPSTL formula

[TABLE]

Note that (7) involves both comparison of probabilities and nesting of probability operators quantifying different paths, which are not allowable in common non-hyper temporal logics.

5. Statistical Verification of HyperPSTL Properties

In this section, we study the SMC of HyperPSTL on a PUS with given inputs. As with previous works, we focus on handling the probability operators in HyperPSTL by sampling, which is the main issue for the SMC of probabilistic temporal logic. Through Sections 5.1 – 5.3, we propose SMC algorithms based on Clopper-Pearson (CP) significance level calculation for all the ways the probability operators can be used or nested in HyperPSTL. Accordingly, any nested HyperPSTL can be verified recursively by applying these SMC algorithms. The handling of temporal operators is similar to that of hyperSTL (Nguyen et al., 2017) and is thus not fully discussed due to the space limitations. For a PUS, verifying bounded-time properties is straightforward; verifying unbounded-time properties is more involving, and will be part of the future work.

5.1. SMC via CP Significance Level

A new feature of HyperPSTL compared to common temporal logics is the simultaneous probability quantification over multiple path variables. We now illustrate the idea of SMC for such formulas. Consider HyperPSTL formula $\Phi=(\mathbin{\mathbb{P}}^{\Pi}\varphi<p)$, where (i) $\Pi=\{\pi_{1},\ldots,\pi_{K}\}=\mathsf{fv}(\varphi)$ is the set of free path variables of $\varphi$, (ii) $p\in[0,1]$ is a probability threshold, and (iii) $\varphi$ contains no probability operator. The semantics of $\Phi$ is

[TABLE]

The truth value of $\varphi$ can be evaluated on a set of concrete sample paths $\boldsymbol{\sigma}=\{\sigma_{1},\ldots,\sigma_{K}\}$ by assigning the concrete sample path $\sigma_{i}$ to the free path variable $\sigma_{i}$ for $i\in[K]$. Hence, with a slight abuse of notation, we denote

[TABLE]

Previous SMC approaches have used the sequential probability ratio test (SPRT) to evaluate $\Phi$ with the specification of an indifference margin (Legay and Viswanathan, 2015). Specifically, assuming that

[TABLE]

for some $\delta>0$, to evaluate $\Phi$, it suffices to test the two most indistinguishable cases, i.e., a Simple Hypothesis Testing (SHT) problem with two hypothesis,

[TABLE]

which can then be solved by SPRT (Hogg et al., 2005).

Since the choice of indifference margin is somewhat arbitrary, we propose an indifferent margin-free SMC approach via significance level calculation. For $i\in[N]$ and $K=|\Pi|$, let $(\sigma^{(i)}_{1},...,\sigma^{(i)}_{K})$ be a tuple of i.i.d. sample paths drawn from the PUS $\mathcal{S}$ starting from the state $X$. Checking the correctness of $\varphi$ by (8) on each tuple gives the sum statistic

[TABLE]

that obeys the binomial distribution $\mathrm{Binom}(n,p_{\varphi})$. The average statistics $T/N$ is a unbiased estimator for $p_{\varphi}$. Intuitively, when $T/N<p$, it is more likely that $p_{\varphi}<p$; and the same for the other case. Hence, we define the following statistical asserting function based on the samples

[TABLE]

To ensure the asymptomatic correctness of the SMC algorithm, we assume that

[TABLE]

which is a weaker assumption than (9). When (12) holds, as the number of samples increases, the samples will be increasingly concentrated on one side of $p$ by the central limit theorem. Therefore, a statistical analysis based on the majority of the samples has an increasing accuracy. When (12) is violated, the samples would be evenly distributed on the two sides of $p$, regardless of the sample size. Thus, no matter how the sample size increases, the accuracy of any statistical test would not increase. This will be illustrated later in the proof of Theorem 5.1.

In general, the significance level for claiming $p_{\varphi}\in[a,b]\subseteq[0,1]$ (i.e., an upper bound of the probability of making a wrong claim), when $T/N\in[a,b]$, can be computed using a method from Clopper and Pearson (Clopper and Pearson, 1934) by

[TABLE]

where $F_{\text{Beta}}(\cdot\,|\,T_{1},T_{2})$ is the cumulative probability function (CDF) of the beta distribution $\text{Beta}(T_{1},T_{2})$ with the shape parameters $(T_{1},T_{2})$. For computing the significance level of the assertion $\mathcal{A}((\mathcal{S},X)\models\Phi)$, we utilize

[TABLE]

Remark 1.

The CP significance levels (13) have the following properties. First, $\alpha_{\text{CP}}(a,b\,|\,T,N)$ increases as the interval $[a,b]$ shrinks. That is, for any $[a^{\prime},b^{\prime}]\subseteq[a,b]\subseteq[0,1]$, we have $\alpha_{\text{CP}}(a^{\prime},b^{\prime}\,|\,T,N)\leq\alpha_{\text{CP}}(a,b\,|\,T,N)$. In addition, for $T\notin\{0,N\}$ and $N\gg 1$, we have that

[TABLE]

Since the beta distribution $\text{Beta}(T,N-T)$ has the mean $T/N$ and the variance $T(N-T)/N^{2}(N+1)$, for fixed $T/N$, as the number of samples $N\rightarrow\infty$, the beta distribution becomes increasingly concentrated at $T/N$, and thus $\alpha_{\text{CP}}(a,b\,|\,T,N)\rightarrow 0$. This implies that the probability of making the wrong claim decreases as more samples are available.

Given a desired significance level $\alpha$, we can design a new SMC algorithm by, at each iteration, collecting $B$ new samples, computing the CP significance interval and stopping when the result is less than $\alpha$, as summarized in Algorithm 1. Correctness of Algorithm 1 follows directly from the definition of significance level.

Theorem 5.1.

Algorithm 1* terminates with probability $1$ and gives the correct statistical assertion with probability at least $1-\alpha_{d}$.*

Proof.

Termination: For $p_{\varphi}\in\{0,1\}$, the proof is trivial. For $p_{\varphi}\in(0,1)$, from assumption (12) and without loss of generality, let $p_{\varphi}<p$ and $\delta=p-p_{\varphi}$. Recalling the second half of Remark 1, for any $T/N\in[p_{\varphi}-\delta/2,p_{\varphi}+\delta/2]\subseteq(0,p)$, the variance of $\text{Beta}(T,N-T)$ is lower bounded by $\min_{x\in[p_{\varphi}-\delta/2,p_{\varphi}+\delta/2]}x(1-x)/(N+1)$. Therefore, as $N\rightarrow\infty$, it uniformly converges to [math]. This implies that $\alpha_{\textrm{CP}}(0,p\,|\,T,N)$ uniformly converges to [math] – i.e., for any given $\alpha_{d}>0$, there exists $N_{0}(p_{\varphi},\delta)\in\mathbb{N}$, such that $\alpha_{\text{CP}}(0,p\,|\,T,N)<\alpha_{d}$ for any $N\geq N_{0}(p_{\varphi},\delta)$ and $T/N\in[p_{\varphi}-\delta/2,p_{\varphi}+\delta/2]\subseteq(0,p)$.

With $N\geq N_{0}(p_{\varphi},\delta)$, by the law of large numbers (or central limit theorem), we have \mathbf{Pr}\big{(}T/N\in[p_{\varphi}-\delta/2,p_{\varphi}+\delta/2]\big{)}\rightarrow 1, as the number of samples $N\rightarrow\infty$. Therefore, Algorithm 1 terminates with probability $1$.

Correctness: Let $\tau$ be the step Algorithm 1 terminates and $A$ be “the assertion $\mathcal{A}$ in (11) is correct”, then $\mathbf{Pr}(A)=\sum\nolimits_{i\in\mathbb{N}}\mathbf{Pr}(A|\,\tau=i)\mathbf{Pr}(\tau=i)$. By construction of the significance intervals, for any $i\in\mathbb{N}$, we have $\mathbf{Pr}(A\allowbreak\,|\,\tau=i)>1-\alpha_{d}$. In addition by Termination, we have $\sum_{i\in\mathbb{N}}\mathbf{Pr}(\tau=i)=1$, Thus, $\Pr(A)\geq 1-\alpha_{d}$. ∎

5.2. SMC of Joint Probabilities

Another new feature of HyperPSTL is the arithmetics and comparisons of the probabilities of multiple sub-properties. For example, we can compare the satisfaction probability of $\varphi_{1}$ and $\varphi_{2}$ by the HyperPSTL formula $p_{1}<p_{2}$, where $p_{1}=\mathbin{\mathbb{P}}^{\Pi_{1}}\varphi_{1}$ and $p_{2}=\mathbin{\mathbb{P}}^{\Pi_{2}}\varphi_{2}$, according to the syntax (1), (2). For simplicity, let $\Pi_{1}=\mathsf{fv}(\varphi_{1})$ and $\Pi_{2}=\mathsf{fv}(\varphi_{2})$. Mathematically, this is equivalent to reasoning over the joint probabilities of these properties. Specifically, $p_{1}<p_{2}$ can be equivalently expressed as a specification on the joint probability $(p_{1},p_{2})\in D$, where D=\big{\{}(x_{1},x_{2})\in[0,1]^{2}\ |\ x_{1}<x_{2}\big{\}}.

To formally capture this, we introduce an additional syntactic rule $(p_{1},\ldots,p_{n})\in D$ in (1), whose semantics is given by

[TABLE]

where $D\subseteq[0,1]^{n}$ is measurable.

While the expressiveness of HyperPSTL is unchanged with the new rule (16), the conjunction and disjunction of several HyperPSTL formula can be simplified. For example, the HyperPSTL formula $\Phi_{1}\land\Phi_{2}$ with

[TABLE]

can be equivalently written as $(\mathbin{\mathbb{P}}^{\Pi_{1}}\varphi_{1},\mathbin{\mathbb{P}}^{\Pi_{2}}\varphi_{2},\mathbin{\mathbb{P}}^{\Pi_{3}}\varphi_{3})\in D$, where

[TABLE]

In addition, the new rule simplifies the SMC of the specification. Previously, it requires checking both $\Phi_{1}$ and $\Phi_{2}$ separately, and then a probabilistic composition of the two results. With the new rule, a single procedure of checking whether the joint probability $(\mathbin{\mathbb{P}}^{\Pi_{1}}\varphi_{1},\mathbin{\mathbb{P}}^{\Pi_{2}}\varphi_{2},\mathbin{\mathbb{P}}^{\Pi_{3}}\varphi_{3})$ is in $D$ is sufficient.

We now demonstrate the idea of verifying the joint probability by checking a non-nested state formula

[TABLE]

where for $i\in[n]$, and $\varphi_{i}$ contains no probability operator and $\Pi_{i}=\mathsf{fv}(\varphi_{i})$. As a statistical approach is adopted, we assume that the exact probability of satisfying $\varphi$ does not lie within the boundary of the test region $D$, as stated in Assumption 1.

Assumption 1.

To check (\mathcal{S},X)\models\big{(}\allowbreak\mathbin{\mathbb{P}}^{\Pi_{1}}\varphi_{1},\allowbreak\ldots,\allowbreak\mathbin{\mathbb{P}}^{\Pi_{n}}\varphi_{n}\big{)}\in D, we assume that (i) the test region $D$ is a simply connected domain with $\mu_{\mathrm{Borel}}(D)\neq 0$, and (ii)

[TABLE]

Assumption 1 can be viewed as the multidimensional generalization of (12), and is necessary for asymptomatic correctness of the SMC algorithm, as discussed in Section 5.1.

Remark 2.

Compared to previous studies on SMC using sequential probability ratio tests (SPRT) (Zuliani, 2015; Wang et al., 2018), Assumption 1 is weaker as it requires no a priori knowledge on the indifference margin.

Remark 3.

From Assumption 1, we have that \big{(}\allowbreak\mathbin{\mathbb{P}}^{\Pi_{1}}\varphi_{1},\allowbreak\ldots,\allowbreak\mathbin{\mathbb{P}}^{\Pi_{n}}\varphi_{n}\big{)}\allowbreak\in D_{1} and \big{(}\allowbreak\mathbin{\mathbb{P}}^{\Pi_{1}}\varphi_{1},\allowbreak\ldots,\allowbreak\mathbin{\mathbb{P}}^{\Pi_{n}}\varphi_{n}\big{)}\in D_{2} are semantically equivalent if $\overline{D_{1}}=\overline{D_{2}}$. In addition, $\mathbin{\mathbb{P}}^{\Pi}\in D$ and $\mathbin{\mathbb{P}}^{\Pi}(\neg\varphi)\in D^{c}$ are semantically equivalent when $D^{c}=[0,1]\backslash D$.

By the semantic rule (16), the SMC problem is converted to a composite hypothesis testing problem

[TABLE]

For each $i\in[n]$, let $\{\boldsymbol{\sigma}_{j}\}_{j\in[N_{i}]}$ be $N_{i}$ tuples of i.i.d sample paths of the PUS $\mathcal{S}$ starting from the same state $X$ that are used to estimate $p_{\varphi_{i}}$. Similar to Section 5.1 approach, we consider the statistics

[TABLE]

where $T_{i}\sim\mathrm{Binom}(N_{i},p_{\varphi_{i}})$. We define the assertion using the average statistics $T_{i}/N_{i}$ for $i\in[n]$ as

[TABLE]

For this multi-dimensional case, to compute the exact CP significance level for a general domain $D$ involves multi-dimensional integrations on it, which can be computationally intensive. Therefore, we compute an upper bound $\alpha_{\mathrm{CP}}$ on the significance level of the assertion (20) by finding a hypercube $\prod_{i\in[n]}[a_{i},b_{i}]$ such that

[TABLE]

Due to the monotonicity of significance levels, the significance level of $(T_{1}/N_{1},...,T_{n}/N_{n})\in D$ is upper bounded by that of $(T_{1}/N_{1},...,\allowbreak T_{n}/N_{n})\in\prod_{i\in[n]}[a_{i},b_{i}]$, which can be computed directly by compositing the significance level of $T_{i}/N_{i}\in[a_{i},b_{i}]$ for $i\in[n]$, using the results in Section 5.1. Thus, we compute an upper bound $\bar{\alpha}_{\mathrm{CP}}$ of the exact CP significance level as

[TABLE]

where $\alpha_{\mathrm{CP}}(a_{i},b_{i}\,|\,T_{i},N_{i})$ is defined in (13).

When implementing this verification approach, for each iteration, we look for a hypercube $\prod_{i\in[n]}[a_{i},b_{i}]$ satisfying (21) with $a_{i}<b_{i}$ for $i\in[n]$. Although finding such a hypercube is only possible if $(T_{1}/N_{1},...,T_{n}/N_{n})\notin\partial D$, this is guaranteed with probability $1$ for large samples when Assumption 1 holds.

To minimize the upper bound of the significance level $\bar{\alpha}_{\mathrm{CP}}$ in (22), the hypercube should be preferably as large as possible. For a simple domain $D$, the analytic solutions of such a largest hypercube can be derived directly as a function of $T_{i},N_{i}$ for $i\in[n]$ and the functions defining $\partial D$ More generally, especially if $D$ is convex, the largest hypercube can be derived by solving the optimization problem of maximizing its volume while keeping it inside $D$. Admittedly, solving the optimization problem at every iteration can still be inefficient for some cases. To remedy for this, we can (1) reduce the frequency of computing the significance level by drawing samples in batches; and (2) search only for approximate maxima in optimization.

Finally, we note that the upper bound $\bar{\alpha}_{\mathrm{CP}}$ is asymptotically tight if a largest hypercube is used to compute it. This holds because, by the law of large numbers, as the number of samples increases, $(T_{1}/N_{1},...,T_{n}/N_{n})$ concentrates near $(p_{\varphi_{1}},\ldots,p_{\varphi_{n}})$ and the largest hypercube converges to a constant one strictly containing $(p_{\varphi_{1}},\ldots,\allowbreak p_{\varphi_{n}})$. Thus, the probability of $(T_{1}/N_{1},...,\allowbreak T_{n}/N_{n})\in\prod_{i\in[n]}[a_{i},b_{i}]$ converges to that of that of $(T_{1}/N_{1},...,\allowbreak T_{n}/N_{n})\in D$.

Based on the previous discussions, we derive Algorithm 2. Correctness of Algorithm 2 is given by Theorem 5.2 that can be proved in the same way as Theorem 5.1.

Theorem 5.2.

Algorithm 2* terminates with probability $1$ and gives the correct statistical assertion with probability at least $1-\alpha$.*

5.3. SMC of Nested Probability Operators

Finally, we consider SMC of nested HyperPSTL formulas on the PUS $\mathcal{S}$. By the syntax of HyperPSTL in Section 3.1, a nested HyperPSTL formula is constructed iteratively in two ways: (i) replacing an atomic proposition with a general state formula, and (ii) consecutively nested quantification of free path variables in a non-nested formula. The former also appears in common temporal logics and the latter is unique to HyperPSTL.

For (i), we show the idea of SMC by checking the satisfaction on a state $X$ of the nested formula $\Psi=\mathbin{\mathbb{P}}^{\Pi}\psi[\rho]\in D_{1}$ derived by replacing an atomic proposition of $\psi$ with a non-nested state formula \rho=\big{(}\allowbreak\mathbin{\mathbb{P}}^{\Pi_{1}}\varphi_{1},\allowbreak\ldots,\allowbreak\mathbin{\mathbb{P}}^{\Pi_{n}}\varphi_{n}\big{)}\in D_{2}. In $\Psi$, we have $D_{1}\subseteq[0,1]$, $D_{2}\subseteq[0,1]^{n}$. If $\rho$ is treated as an atomic proposition, $\Psi$ becomes a non-nested state HyperPSTL formula as discussed in Section 5.1. The HyperPSTL state formulas nesting more probability operators in this fashion can be statistically verified in the same way.

To verify $\Psi$, we follow a compositional analysis similar to (Sen et al., 2005a, b). If the state space $\mathcal{X}$ of the PUS is finite (e.g., a CTMC from Section 2.2), we can statistically verify the sub-formula $\rho$ on each state $X$ of the PUS $\mathcal{S}$ with significance level $\alpha_{X}$, using Theorem 5.2 and Algorithm 2, and label the state with $\rho$ if the assertion given by Algorithm 2 is \mathcal{A}\big{(}V\models\rho\big{)}=1. Then, we can statistically verify the full formula $\Psi$ on the relabeled PUS as a non-nested formula with significance level $\alpha_{0}$ using Theorem 5.2. The overall significance level is $\alpha=\sum_{X\in\mathcal{X}}\alpha_{X}+\alpha_{0}$, which is only bounded when $\mathcal{X}$ is finite. The SMC for $\Psi$ on infinite-state PUS provides an avenue for future work. For this work, it brings no limitation as we only verify such nested formulas on finite-state PUS (e.g., queueing in Section 6), while for P2HIOA such formulas are not used to capture properties of interest.

To implement the SMC algorithm for $\Psi$ on finite-state PUS, given the overall significance level $\alpha$, we can split it into the summation $\sum_{X\in\mathcal{X}}\alpha_{X}+\alpha_{0}$. The simplest way is $\alpha_{X}=\alpha_{0}=\alpha/(|\mathcal{X}|+1)$ for $X\in\mathcal{X}$, where $\mathcal{X}$ is the number of states of the PUS $\mathcal{S}$. Then, we can employ Algorithm 2 to verify $\rho$ on each state $X$ with significance level $\alpha_{X}$ and then assert $\Psi$ with significance level $\alpha_{0}$ using Algorithm 2 again. This is summarized by Algorithm 3.

For (ii), we show the idea of SMC by checking the satisfaction on a state $X$ of a HyperPSTL formula $\Psi=\mathbin{\mathbb{P}}^{\Pi_{1}}(\mathbin{\mathbb{P}}^{\Pi_{2}}\varphi<p_{2})<p_{1}$, where the sub-formula $\varphi$ contains no probability operator and all its path variables are probabilistically quantified by $\mathbin{\mathbb{P}}^{\Pi_{1}}\mathbin{\mathbb{P}}^{\Pi_{2}}$, i.e., $p_{1},p_{2}\in[0,1]$, $\Pi_{1}\cap\Pi_{2}=\emptyset$, and $\mathsf{fv}(\varphi)=\Pi_{1}\cup\Pi_{2}$.

The SMC can be similarly done for HyperPSTL state formulas that nest more probability operators in this fashion. The formula $\Psi$ says that with probability at most $p_{1}$, we can find a set of paths $\Pi_{1}$ such that the probability to find another set of paths $\Pi_{2}$ to satisfy $\varphi$ is at most $p_{2}$. This formula can be equivalently expressed by \Psi=\mathbin{\mathbb{P}}^{\Pi_{1}}\big{(}\mathbin{\mathbb{P}}^{\Pi_{2}}\in[0,p_{2}]\big{)}\in[0,p_{1}] using the rule (16). We note that for (ii), unlike (i), the state space $\mathcal{X}$ of the PUS can be infinite.

The main idea is as follows. For $i\in[N]$, let $\boldsymbol{\sigma}_{i}$ be a $|\Pi_{1}|$-tuple of i.i.d. sample paths starting from $X$ in the PUS; define the indicator

[TABLE]

of whether the partly instantiated formula $(\mathcal{S},X)\models\llbracket\mathbin{\mathbb{P}}^{\Pi_{2}}\varphi<p_{2}\rrbracket_{V[\Pi_{1}\rightarrow\boldsymbol{\sigma}_{i}]}$ is true under this instantiating. Although $T_{i}$ is not directly accessible, we can estimate it statistically using the assertion $A_{i}$ of Algorithm 2 for any given significance level $\alpha_{1}>0$.

To verify the full formula $\Psi$, we only need to estimate the total number of positive instantiation $T=\sum_{i\in[N]}T_{i}$. We estimate it using $A=\sum_{i\in[N]}A_{i}$. Since $A_{i}\neq T_{i}$ with probability at most $\alpha_{1}$ for all $i\in[N]$, we have $|T-A|<\Delta$ with the significance level

[TABLE]

where $F_{\mathrm{Binom}}$ is the Binomial cumulative distribution function. Since

[TABLE]

we can check the full formula $\Psi$ using these minimal and maximal estimations of $T$. Intuitively, if $T_{2}/N<p$ (hence $T_{1}/N<p$), it is more likely that $\Psi$ is true; if $T_{1}/N>p$ (hence $T_{2}/N>p$), it is more likely that $\Psi$ is false; otherwise, further sampling is needed. Thus, we define the following statistical asserting function by

[TABLE]

When a final assertion is possible, its significance level is the larger one between plugging $T_{1}$ and $T_{2}$ into (13) and (14). Accordingly, the overall significance level $\alpha$ is

[TABLE]

where $\alpha_{\mathrm{CP}}$ is given by (13).

To implement the SMC algorithm, given the overall significance level $\alpha$, we need to simultaneously decrease both the significance level $\alpha_{1}$ for making assertions on the partly instantiation of $\psi$ for given values of $\Pi_{1}$, and the significance level $\alpha_{2}$ for estimating the sum of those assertions. We start from $\Delta=c\alpha_{1}N$ with $c=1$ and $\alpha_{1}=\alpha_{d}$. When $\alpha_{2}$ is the main source of statistical error, i.e., $\alpha_{2}>\alpha/2$, we decrease $\alpha_{2}$ by increasing the parameter $c$ by $1$; otherwise, we decrease $\alpha_{1}$ by reducing it by half. This is summarized in Algorithm 4.

6. Evaluation

We numerically evaluate our SMC algorithms on several benchmarks with different complexity levels. Other probabilistic hyperproperties on different systems are handled in a similar manner, but due to space constraints here we focus on the discussed properties/systems. All benchmarks are implemented in Matlab/Simulink and are available in (CPSL@Duke, 2019). Specifically, the Toyota powertrain model is derived from (Jin et al., 2014a); and the queueing networks are implemented in Simulink using the SimEvents Toolbox (The MathWorks, Inc., 2019). Evaluations are performed on a laptop with 16 GB RAM and Intel Xeon E-2176 CPU.

For each benchmark, we evaluate the proposed SMC algorithms in different setups by changing the desired significance level $\alpha$, as well as the parameters (e.g., $\delta$, $\varepsilon$, and $t$) in the objective HyperPSTL specifications. The proposed SMC algorithms are executed repeately on each setup for $100$ times. This is to check whether the probability for the proposed SMC algorithms to return the correct assertion is at least $1-\alpha$: i.e., we repeat the SMC algorithm for each setup for $100$ times, and check if it makes the correct assertion for at least $100(1-\alpha_{d})$ times.

For each setup, to compare it with the assertions of the proposed SMC algorithms, the truth value of the HyperPSTL specification of interest is derived by estimating the probabilities involved in it using an analytic solution or numerous sampling. Specifically, for the thermostat, we derive an analytic solution for the left hand side of (5) from its dynamics; this can be done due to its simplicity. For the powertrain, we estimate the left hand side of (5) by sampling $10^{5}$ pairs of $(\pi_{1},\pi_{2})$, for which the standard error is less than $0.01$. For the queueing networks, we estimate the left hand side of (7) by drawing $500$ samples for $\pi_{1}$. For each sample of $\pi_{1}$, we draw $500$ samples for $\pi_{2}$ to evaluate the truth value of \big{|}\mathbf{Pr}_{\pi_{2}}(\tau_{i}^{\pi_{1}}-\tau_{j}^{\pi_{2}}>t)-\mathbf{Pr}_{\pi_{2}}(\tau_{j}^{\pi_{2}}-\tau_{i}^{\pi_{1}}>t)\big{|}<\delta. The total standard error for estimating the left hand side of (7) is less than $0.05$.

Results for all considered setups are shown in Tables 1, 2, 3 and 4; as can be seen, the estimated accuracy of our SMC algorithms are very close to $1$, showing the conservativeness of the CP significance level. We also report the average number of sample paths and the execution times for each setup based on the $100$ repetitions, rounded up to their standard errors. In all simulations, the number increases when the desired significance level decreases, showing the trade-off between accuracy and the sampling cost. In addition, the execution time is mainly consumed by drawing samples from the Simulink models, and is approximately propositional to the number of samples.

6.1. Thermostat

A thermostat can be modeled as a simple P2HIOA with two modes $\mathsf{Heat}$ and $\mathsf{Cool}$ (Figure 6), and one state variable $T$ that varies within the temperature interval $[T_{l},T_{h}]=[15,40]\subseteq$. The mean heating and cooling rates are $c_{1}=c_{2}=5$; they are subject to time-invariant but random Gaussian error $n_{1},n_{2}\sim\mathbf{N}(0,0.5^{2})$. The thermostat starts from \big{(}T=T_{l},\mathsf{Heat}\big{)}. We verify the sensitivity of the running period of a heat and cool cycle under the noise, which is represented in HyperPSTL by (5) with \mathtt{Q}:=\big{(}T=T_{l},\mathsf{Cool}\big{)}. We statistically verified the sensitivity specification using Algorithm 2, with specification (5) parameters $\delta\in\{0.9,1.1\}$ and $\varepsilon\in\{0.05,0.01\}$, under the desired significance levels $\alpha\in\{0.01,0.05\}$.

The derived results in Table 1 give accurate estimations (with significance level as low as $\alpha=0.01$) on the probability distribution of the sensitivity, with a relatively small number of samples (at most a few hundred samples for each setup). We verified that the sensitivity of the running period of the thermostat is less than $\delta=1.1$ with probability $1-\varepsilon=0.95$, but not less than $\delta=0.9$ with the same probability, showing that the $\mathbf{0.95}$ percentile is between $\mathbf{[0.9,1.1]}$. Also, as the sensitivity specification is false for $1-\varepsilon=0.99$ for both $\delta=0.9$ and $\delta=1.1$, showing that the $\mathbf{0.99}$ percentile is in $\mathbf{[1.1,\infty]}$.

6.2. Toyota Powertrain Control System

We use the Simulink model for the Toyota Powertrain with a four-mode embedded controller from (Jin et al., 2014b). It can be considered as a P2HIOA  with four modes and $15$ state variables. We consider the sensitivity of the recovery time after start of the deviation percentage of the air/fuel (A/F) ratio $\mu$ to the level $|\mu|<0.05$, under the mean RPM $2500$ subjecting to Gaussian noise $\mathbf{N}(0,25^{2})$. The sensitivity specification is formally expressed by HyperPSTL formula (5) with \mathtt{Q}=\big{(}|\mu|<0.05\big{)}; we statistically verified it using Algorithm 2 with parameters $\varepsilon\in\{0.05,0.01\}$ and $\delta\in\{0.06,0.07\}$, under the desired significance levels $\alpha\in\{0.01,0.05\}$.

The results shown in Table 2 give accurate estimations (with significance level as low as $0.01$) on the probability distribution of the sensitivity under the given embedded controller; this is achieved with a relatively small number of samples (at most a few hundred samples for each setup). We verified that the sensitivity of the recovery time of the powertrain is less than $\delta=0.15s$ with probability $1-\varepsilon=0.99$, but not less than $\delta=0.20s$ with the same probability, showing that the $\mathbf{0.99}$ percentile is between $\mathbf{[0.15s,0.20s]}$. Also, as the sensitivity specification is true for $1-\varepsilon=0.95$ for both $\delta=0.15s$ and $\delta=0.20s$, showing that the $\mathbf{0.95}$ percentile is in $\mathbf{[0,0.15s]}$.

6.3. Queueing Networks

We consider two queuing networks with different sizes $m$ and $n$, as shown in Figure 5: the Small has $1$ front-end servers and $2$ back-end servers, the Large has $25$ front-end servers and $20$ back-end servers. For the small model, the package arrival and processing are modeled by exponential distribution, while in the second model, the package arrival and processing are modeled by Markov Modulated Poisson Processes with different parameters over all queues. Note that the parameters for different front-end servers are not identical, so deriving an exhaustive solution is virtually impossible for large sizes. We consider the fairness of workloads in queuing networks under the policy that the front servers deliver to the back server with the shortest queue. The fairness specification is formally defined in HyperPSTL by (7) with $\mathtt{Q}_{i}$ representing the overloading of the back server $i$. We set $i=1$, $j=2$, and statistically verify the specification with Algorithm 4 with $t\in\{0.1,5.0\}$, $\delta\in\{0.1,0.5\}$ and $\varepsilon\in\{0.1,0.5\}$, for the significance level $\alpha=0.95$.

From the first row of Tables 3 and 4, for both the small and large queueing network, with probability $1-\varepsilon=0.9$ for the back server $1$, we have that $|p_{1}-p_{2}|<0.1$, where $p_{1}$ and $p_{2}$ are the probabilities that back server $2$ overloads $t=5$ earlier/later than back server $1$, respectively. This shows that for $90\%$ cases, the shortest-queue-delivery policy is roughly fair for back-end server $1$, in the sense that its overload time is not significantly sooner or later (for $t=5$) than back-end server $2$. On the other hand, from the third row of Tables 3 and 4, with probability $1-\varepsilon=0.9$ for back-end server $1$, we have $|p_{1}^{\prime}-p_{2}^{\prime}|>0.5$, where $p_{1}^{\prime}$ and $p_{2}^{\prime}$ are the probabilities that back server $2$ overloads $t=0.1$ earlier/later than back server $1$, respectively. Thus, for $90\%$ cases, the shortest-queue-delivery policy is not exactly fair for back server $1$, in the sense that its overload time can be moderately sooner or later (for $t=0.1$) than for back-end server $2$.

7. Conclusion

In this work, we have studied statistical verification of hyperproperties for cyber-physical systems (CPS). We have first defined a general model of probabilistic uncertain systems (PUS) that unify commonly studied modeling formalisms such as continuous-time Markov chains and hybrid I/O automata with probabilistic dynamical parameters. We have then introduced hyperproperties, such as fairness and sensitivity, that involve relationships between multiple paths simultaneously in continuous time. To formally specify such hyperproperties, we have introduced Hyper Probabilistic Signal Temporal Logic (HyperPSTL), which is a hyper and probabilistic version of the conventional the signal temporal logic (STL). To simplify our presentation of the HyperPSTL syntax and semantics, while considering SMC of specific hyperproperties of complex CPS, in the proposed HyperPSTL, we have disallowed nested existential and universal quantifies on states to avoid exhaustive iteration on the possibly infinite state space. Still, this logic can be augmented (as done in detail in (Wang et al., 2019)) by allowing nested existential and universal quantifications over multiple states.

To verify HyperPSTL specifications on the PUS, we have developed statistical model checking (SMC) algorithms with three new features: (1) the significance level of the HyperPSTL specifications are computed directly using the Clopper-Pearson significance level; (2) statistically verifying HyperPSTL specifications on the joint probabilistic distribution of multiple paths, and (3) HyperPSTL specifications with nested probabilistic operators quantifying different paths are allowed. Finally, we have evaluated the introduced SMC algorithms on different CPS benchmarks with varying levels of complexity.

The work in this paper opens many new avenues for further research. Our work has direct application in software doping (D’Argenio et al., 2017) and in particular, identifying whether two black-box systems are statistically identical. This problem is in particular challenging in the CPS domain, since the physical environment may make the behavior of cyber components more unpredictable. Another application area of our work is in conformance-based testing (Abbas et al., 2014). We also plan to expand our work to analyzing information-flow security and in particular differential privacy.

Acknowledgements.

This work is sponsored in part by the ONR under agreements N00014-17-1-2504, AFOSR under award number FA9550-19-1-0169, as well as the NSF CNS-1652544 and NSF SaTC-1813388 grants.