In this paper, we focus on verifying Cyber-Physical System via Statistical Model Checking and numerical simulation. Our verification problem is to establish whether the expected value of a given Key Performance Indicator exceeds the desired threshold. We propose to use an optimal approximation algorithm to get an (, )-approximation of a desired expected value through Monte Carlo experiments. We prove the feasibility of our approach on the Pumping System, an example derived from the Modelica Standard Library. We were able to verify the system in a reasonable time, taking advantage of a parallel implementation of the approximation algorithm.
hypotheses with a given statistical confidence, whereas quantitative properties provide a measurable interpretation of some stochastic behaviour.
The literature (see, e.g., [1, 49]) shows a plethora of tools (e.g., MultiVeSta [50], PRISM [26], UPPAAL-SMC [14], COSMOS [6], Ymer [60]) that enable SMC via numerical simulation. Many SMC approaches, even based on such tools, require the system modelling through some kind of structure (e.g., Discrete Time Markov Chain, Continuous Time Markov Chain [61, 4], Probabilistic Timed Automata [48]), and a specification through a class of logic (e.g., Probabilistic Computational Tree Logic [23], Continuous Stochastic Logic [52]) to carry out the verification. We note that specific modelling of the system is needed to generate trajectories on demand.
In this paper, instead, we focus on verifying CPSs by relying on a purely black-box approach where the System Under Verification (SUV) is available only via a simulator (e.g., OpenModelica, Dymola, Simulink). Our setting requires establishing if the expected value of a given Key Performance Indicator (KPI) does not exceed a given threshold. We propose to estimate such an expected value using an optimal approximation algorithm (namely, the Approximation Algorithm, [13]) that interfaces in a black-box fashion [55] with the SUV. Mainly, the approximation algorithm asks for a new sample when needed and outputs the estimate when enough samples have been collected. So, our approach aims to be compliant with any simulation driver and SUV model. Furthermore, we prove the feasibility of such an approach on a Modelica model. Note that, to the best of our knowledge, no other contributions exist that apply an estimation algorithm on a Modelica system for purely verification purposes.
Black-box approaches distant from ours have been proposed in, e.g., [51, 59]; in those works, the authors use hypothesis testing to verify qualitative properties when trajectories can not be generated on demand.
The paper is organised as follows. We provide a formalisation of the verification problem in Section 2 and a description of the tool in Section 3. In Section 4, we present the case study used to prove the feasibility of our approach, whereas, in Section 5, we summarise and discuss the results obtained. Finally, we draw our conclusions in Section 6.
We denote with R, R0+, and R+, respectively, the sets of all, non-negative and positive reals. N and N+ are the sets of non-negative and strictly-positive integers.
The system verification via SMC consists of three components: the SUV, the specifications, and the SUV environment. In our setting, a SUV is typically a simulable CPS model, i.e., a system described by continuous as well as discrete variables. The continuous part describes the physical component, whereas the discrete part models the dynamics of the software component. We treat (see Definition 1) our SUV as a black-box input-output deterministic causal dynamical system (see, e.g., [56, 32, 38]) that depends on its inputs and outputs and can be evaluated only through simulation. are the input and output spaces and :
A black-box SUV is a tuple ( , , , ) where is the time-set,
and → is the output (simulation) function that defines the time evolution of the system for each time point in
A black-box SUV is strictly causal if, for any ∈ and for any u1, u2 ∈ such that u1( ′) = u2( ′) for all ′ < , it holds ( ; u1) = ( ; u2) where ( ; u) denotes the . value (u) at time .
The SUV environment defines all operational scenarios the SUV must withstand (e.g., see [31, 36, 45]). Such inputs are uncontrollable, i.e., not under the system’s control (e.g., faults or disturbances). Since we need to sample the operational scenarios needed for the verification, we focus of finitely parametrisable environments, and sample within such parameter spaces. with ⊆ R ,
∈ N.
A finitely parametrisable SUV environment is a function : →
Of course, not all scenarios have the same occurrence. We model this feature through a probability density function ( ) on the set of environment parameters . This seamlessly includes both discrete and continuous parameter spaces.
The SUV must satisfy several requirements for all possible scenarios. The specifications of such requirements model KPIs on the safety properties of the system. Note that in our framework, an output function y ∈ of the SUV can be a system trajectory (possibly needed to compute a KPI) and the KPI itself. In particular, we are interested in quantitative specifications that give a metric interpretation of the properties of interest.
A quantitative KPI is a function : → [0, 1].
We want to know whether the expected value of a given KPI over the set of operational scenarios is less than or equal to a given threshold . We need to compute the expected value of ( (u)) for any u = ( ) for some probability density function ( ) on . Formally, ∈ , where ranges in [0, 1] with So we can define the verification problem as follows: ︁∫
=
( ( ( ))) ( )
︁∫
( ( ( ))) ( )
−
≤ 0
(
may be, and usually is, infinite and computing all operational scenarios is impossible and, even if they are a finite number, typically prohibitive. To avoid such a scenario explosion, we can use Monte Carlo based approaches that provide 1 − , an (, )-approximation ^ of ∫︀ ( ( ( ))) ( ) , that is, with probability at least (1 − ) ≤ ^ ≤ (1 + ) with 0 < ≤ 1 and 0 < ≤ 1.
We focus on verifying CPSs via Statistical Model Checking. Our goal is to perform Monte Carlo experiments to estimate the expected value of a given KPI and establish whether it exceeds the desired threshold. We may empirically compute the expected value through a Monte Carlo approach that uses an arbitrary large number of samples. Unfortunately, such an approach does not provide any formal guarantees, in terms of error bound and statistical confidence, about the accuracy of the expected value estimate.
To overcome this obstacle, we adopt the optimal Approximation Algorithm ( as our estimator.
computes an ( , )-approximation ˜ of the expected value > 0 of a univariate non-negative random variable , i.e., the computed estimate ˜ expected value of is guaranteed to be within (1 − ) and (1 + ) with probability at least 1 − . The algorithm optimally uses Monte Carlo experiments for the approximation, that is, within a constant factor, it requires the minimum number of experiments to output the estimate. Note that in our setting, we get samples via numerical simulation and a non-optimal estimator could be very expensive. Note that in case an absolute (and not relative) error is sought, other algorithms might seamlessly be used, e.g., based on the Bernstein Inequality [13].
[13].
in [0, 1].
Below we briefly summarise the main concepts of the . For a detailed explanation see is the three-phase algorithm illustrated in Algorithm 1. It computes an approximation ˜ of the mean > 0 of a univariate non-negative random variable distributed an initial estimate ^ within a constant factor
Given the user-defined input parameters ∈ (0, 1] and ∈ (0, 1], the following values are defined: = ( − 2) ≈ .72, ϒ = 4 (2/ )/ 2 and = { 2 , }. uses the ifrst two phases to decide the number of samples required to output ˜ in the third phase. In the first phase, it computes an initial √ ( , / 3)-approximation ^ of the mean using the Stopping Rule Algorithm (SRA). Note that in this phase the number of samples cannot be determined in advance. The second phase (Estimate Variance, EV) provides of the variance 2 using 2 · 2 samples where 2 = ϒ 2 · / ^. ^ is with probability of at least 1 − .
Estimate Mean (EM) is the third phase of the
approximation ˜ of by taking advantage from ^ and ^
samples.
algorithm. It outputs an (,
)with 3 = ϒ 2 · ^
(
( 32 )/ ( 2 ))ϒ 2: ^ ← SRA( 3: ^ ← EV(^ , , ϒ 2) 4: ˜ 5: return ˜
← EM(^ , ^ , ϒ 2) function: SRA(, ) 1: ϒ 1 ← 1 + (1 + )ϒ (1/2, √ ), / 3) function: EV(^ , , ϒ 2) ← 0 ← ϒ 2 · / ^ 3: for i=1 to N do 6: return ˜
function: EM(^, ^, ϒ 2) ← ϒ 2 · ^ ← 0 3: for i=1 to N do
/^2 ← + ← / 6: return ˜ 1: 2: ← +
is implemented as explained in the pseudocode of Section 3.1. The orchestrator updates the approximation algorithm only if new samples have been produced; otherwise, it freezes the algorithm. In such a context, the computational time required to produce new samples via simulation of the SUV model is the main bottleneck of the addressed problem, and the simulation time afects its performance.
Note that computing the expected value of multiple KPIs can be addressed seamlessly by using the same random samples, where, for each sample, the SUV simulator computes the value of all KPIs. Hence, the approximation algorithm terminates when enough samples are simulated so that (, )-approximations of all KPIs are computed.
The workflow of each We developed a parallel implementation of
The (orchestrator) and slaves (workers). The orchestrator handles both the and the input for all workers, whereas each worker waits for new input and simulates the algorithm system to produce a new sample.
The orchestrator keeps the computation status of the algorithm and sends new input to idle workers as long as new samples are required. Our implementation takes into account whether too many samples have been produced at the end of the first phase so that we can use them during the second phase. In such a way, we eficiently use the computational power without wasting any model simulation. The dispatching of the inputs is centralised to sort the obtained samples deterministically. Such a feature is not strictly required, but ensures the full reproducibility of all experiments.
We evaluated the efectiveness of our approach on the Pumping System (PS) model provided by the Modelica Standard Library. PS represents a pumping control system for drinking water. The water is pumped from a source by a pump, fitted with check valves, into a reservoir. The reservoir is described by a (nominal) height of 3 and an area
of 50 2 and it is connected to another valve that models the users, i.e., the outgoing water sink. The source and water sink work at an environmental pressure (101 325 Pa) and temperature (20 ℃) to get the ingoing and outgoing water dynamics where is the average value of atmospheric pressure in standard conditions. In addition,
models a constraint on the system dynamics, that is, the water level must be within [0, ] during the course of the whole simulation.
The control component outputs the rotational speed of the pump engine through a ifrst-order system. The component sends a small non-zero value to put the pump in a standby state when other water is not required, whereas other values indicate that the pump must inject the water into the reservoir.
Each experiment of our Monte Carlo-based approach simulates the system with an horizon ℎ of 2000 seconds. After = 200 seconds, PS opens the sink valve, and the controller starts turning on and of the pump to keep the level water value around 2.2 meters so that the gauge pressure is around 200 mbar.
We added a stochastic environment that injects a disturbance on the environmental pressure so that the ingoing and the outgoing water pressure are unbalanced, i.e., the sink and source stochastic pressure take diferent values. Our stochastic pressure ranges from (1 − ) to (1 + ) where is the percentage of disturbance. We chose a disturbance of = 0.07 (7%) as it is the maximum value that still yields plausible trajectories (for instance, with = 0.08, we witnessed trajectories in which the water level exceeded the tank height).
To verify the correctness of the PS, we modelled the Mean Relative Absolute Error (MRAE) as desired KPI on a current water level to quantify the ability of the system to keep it around 2.2 meters in presence of a stochastic environmental pressure. Given the current water level and reference water level = 2.2, we define the Relative Absolute Error (RAE) at time as follows:
Thus, we define the MRAE at time as follows:
MRAE( ) =
RAE( ) = | ( ) − |
︀∫ 0 RAE( )
(
In this section, we report the experimental results obtained on the PS. The experiments were run on a computer equipped with AMD EPYC 7301 16-Core CPU and 256 GB RAM. We estimated the expected value of the MRAE on the water level on several configurations for and . From preliminary experiments, we have observed that the simulation time drastically increases when we run more than 16 copies of the Modelica system in parallel. In our setting, the simulator is particularly expensive and requires so many computational resources that the simulation time degrades. This motivates the fact that we limited our experiment settings up to 16 workers. 16 workers. It reports the average number of samples required by of the estimated percentage error (E[ ]) and the average value of the time needed to , the average value carry out the experiments.
In each configuration, on average, the PS has a percentage error of 14.7%. This indicates that our safety property is satisfied since our KPI can range on average in [2.2 − 0.323, 2.2 + 0.323] where 0.323 is 14.7% of the
The system simulation time and hyperparameters and impact the experiments. Values for and (mainly ) afect the number of samples needed to obtain the desired estimate. In fact, accurate estimations with higher statistical confidence, i.e., smaller values for and , require more samples. The simulation time, instead, heavily dominates the execution time of our experiments in which we can get each sample only via numerical simulation of the SUV. Such an aspect is not negligible and becomes increasingly relevant when non-optimal approximation algorithms are used. In [47], the authors explain how approximation algorithms that do not take advantage of the variance have worse performance. Note that the number of required samples is the criterion to use to compare diferent Monte Carlo approaches. In such a context, having worse performance with the same and also increases the number of required simulations to estimate the expected value of the target KPI. For example, if a system takes 4-5 seconds for every simulation as in our case, performing too many of simulations could make the problem intractable even when using Statistical Model Checking.
Another observation is about the performance of the developed parallel tool. Figure 2a and Figure 2b illustrate, respectively, speed-up and eficiency by varying the number of workers used for the simulation with = = 0.05. Values for other configurations, even with smaller and , follow the same trend. The speed-up is the ratio between the time required using a single worker and the time required using workers, whereas the eficiency is the ratio between the speed-up with workers and . Both figures show that we are able to eficiently parallelise up to 16 workers with an eficiency greater than 85% and a speed-up grater than 13.5x.
In this work, we performed formal verification of CPSs via Statistical Model Checking. In our setting, solving the verification problem means establishing whether the expected value of a desired KPI does not exceed a given threshold . We developed a parallel implementation of and evaluated our approach’s viability on an example provided by the Modelica Standard Library, the Pumping System. We observed that, even in limited presence of disturbances, the system, on average, is able to satisfy the desired specification. An interesting aspect is the need for parallel approaches when dealing with simulation-based verification. Although and heavily impact the number of required samples by , the parallelisation makes the verification feasible.
Several directions may be pursued in future work. Among those, we mention (
This work was partially supported by: Italian Ministry of University and Research under grant “Dipartimenti di eccellenza 2018–2022” of the Department of Computer Science of Sapienza University of Rome; INdAM “GNCS Project 2020”; Sapienza University projects RG12117A8B393BDC, RG11816436BD4F21, RG11916B892E54DB, RP11916B8665242F; Lazio POR FESR projects E84G20000150006, F83G17000830007. [11] S. Chakraborty, et al. Automotive cyber–physical aystems: A tutorial introduction.
IEEE Des&Test, 33(
CEUR, 2021. [18] M. Esposito and L. Picchiami. Simulation-based synthesis of personalised therapies for colorectal cancer. In OVERLAY 2021, CEUR W.P., vol. 2987. CEUR, 2021. [19] M. Esposito and L. Picchiami. A comparative study of AI search methods for personalised cancer therapy synthesis in copasi. In AI*IA 2022, LNCS, vol. 13196.
Springer, 2022. [20] D. Goswami, et al. Challenges in automotive cyber-physical systems design. In
SAMOS 2012. IEEE, 2012. [21] G. Gottlob, et al. Conditional constraint satisfaction: Logical foundations and complexity. In IJCAI 2007. 2007. [22] R. Grosu and S. Smolka. Monte Carlo model checking. In TACAS 2005, LNCS, vol.
3440. Springer, 2005. [23] H. Hansson and B. Jonsson. A logic for reasoning about time and reliability. Form
Asp Comp, 6(
CAV 2011, LNCS, vol. 6806. Springer, 2011.
[27] B. Leeners, et al. Associations between natural physiological and supraphysiological
estradiol levels and stress perception. Front Psycol, 10, 2019.
[28] F. Maggioli, et al. SBML2Modelica: Integrating biochemical models within
openstandard simulation ecosystems. Bioinformatics, 36(
IEEE, 2015. [37] T. Mancini, et al. User flexibility aware price policy synthesis for smart grids. In
DSD 2015. IEEE, 2015. [38] T. Mancini, et al. Anytime system level verification via parallel random exhaustive hardware in the loop simulation. Microprocessors and Microsystems, 41, 2016. [39] T. Mancini, et al. SyLVaaS: System level formal verification as a service. Fundam
Inform, 149(
Lett, 122, 2017. [41] T. Mancini, et al. Computing personalised treatments through in silico clinical trials.
A case study on downregulation in assisted reproduction. In RCRA 2018, CEUR
W.P., vol. 2271. CEUR, 2018.
[42] T. Mancini, et al. An eficient algorithm for network vulnerability analysis under
malicious attacks. In ISMIS 2018. Springer, 2018.
[43] T. Mancini, et al. Optimal fault-tolerant placement of relay nodes in a mission
critical wireless network. In RCRA 2018, CEUR W.P., vol. 2271. CEUR, 2018.
[44] T. Mancini, et al. Parallel statistical model checking for safety verification in smart
grids. In SmartGridComm 2018. IEEE, 2018.
[45] T. Mancini, et al. Any-horizon uniform random sampling and enumeration of
constrained scenarios for simulation-based formal verification. IEEE TSE, 2021.
[46] I. Melatti, et al. A two-layer near-optimal strategy for substation constraint
management via home batteries. IEEE Trans Ind Elect, 69(
Des, 43(
Systems: a Review. Inf, 11(12), 2020. [50] S. Sebastio and A. Vandin. MultiVeStA: Statistical model checking for discrete event simulators. In ValueTools 2013. ICST/ACM, 2013. [51] K. Sen, et al. "statistical model checking of black-box probabilistic systems". In
CAV 2004, LNCS, vol. 3114. Springer, 2004. [52] K. Sen, et al. On statistical model checking of stochastic systems. In CAV 2005,
LNCS, vol. 3576. Springer, 2005. [53] S. Sinisi, et al. Complete populations of virtual patients for in silico clinical trials.
Bioinformatics, 36(22–23), 2020.
[54] S. Sinisi, et al. Optimal personalised treatment computation through in silico clinical
trials on patient digital twins. Fundam Inform, 174(
Springer, 2005.
[61] H. Younes and R. Simmons. Probabilistic verification of discrete event systems using
acceptance sampling. In CAV 2002, LNCS, vol. 2404. Springer, 2002.
[62] L. Zhang. Modeling automotive cyber physical systems. In DCABES 2013. IEEE,
2013.
[63] P. Zuliani, et al. Bayesian statistical model checking with application to
Statelfow/Simulink verification. Form Meth Sys Des, 43(