Towards a Statistical Model Checking Method for Safety-Critical Cyber-Physical System Verification

College of Computer Science and Technology, Nanjing University of Aeronautics and Astronautics, Nanjing, China Key Laboratory of Safety-Critical Software, Nanjing University of Aeronautics and Astronautics, Nanjing, China Collaborative Innovation Center of Novel Software Technology and Industrialization, Nanjing, China College of Information Engineering, Anhui Finance and Trade Vocational College, Hefei, China


Introduction
Safety-Critical Cyber-Physical System (SCCPS) is characterized with high safety and high reliability and are widely used in fields closely related to the national economy and people's livelihoods, such as aerospace, nuclear industry, public transportation, finance, and medical care. Once the execution of such system fails, it will deeply threaten the safety of human's life and property [1][2][3]. erefore, it is vital to analyze and verify the safety and reliability of safety-critical systems, and it is of great significance to the design and development of safety-critical systems. Indeed, it has attracted wide attention from researchers and has extensively grown as a prominent research topic in the community [4][5][6][7].
Essentially, SCCPS is a kind of complex cyber-physical fusion system [8][9][10]. For this kind of systems, the continuously changing behavior in their physical layer is intertwined with the discrete changing behavior in their decision control layer. eir state spaces are infinite as well. It increases the difficulty and brings severe challenges to the safety analysis and verification of SCCPS. However, the traditional model checking has the problem of state space explosion, and it is difficult to effectively verify it [11].
With the execution path of the sampling system, Statistical Model Checking (SMC) uses statistical analysis techniques to approximate the probability that the target system meets the sequential logic attributes and can provide arbitrarily small error limits [12][13][14]. Because SMC does not need to analyze the complex logic inside the target system to verify the timing logic properties of the system, it can effectively avoid the complexity of the system and the explosion of the state space [15,16]. erefore, SMC is the most effective solution to verify the timing properties of complex SCCPS [12,[17][18][19]. However, for SCCPS requiring extremely high safety, the probability of occurrence of the negative events of its safety attributes and the probability of system failures are extremely low. It is infeasible for SMC to sample extremely low probability events. us, how to use SMC to verify the extremely secure SCCPS is an urgent problem to be solved [20,21].
To date, verification of the SMC rare attributes mainly relies on the importance sampling method. For CTMC and DTMC random models, Reijsbergen et al. [22] and Barbot et al. [23] utilized the heuristic methods to obtain an importance sampling distribution to complete the attribute verification of the two models, respectively. Clarke and Zuliani [24] proposed the cross-entropy minimization importance sampling-based SMC method to verify the safety properties of the Stateflow/Simulink model system. Zuliani et al. [17] used the SMC method in his study [24] to verify the secure attribute of the discrete-time SHS. e methods proposed by Clarke and Zuliani assume that the distribution of the system path space is an exponential distribution. By simply increasing the failure rate of the system parameters, several paths that satisfy the rare attributes are extracted at one time to calculate the optimal parameters for the exponential distribution to obtain an importance sampling distribution [25]. J´egourel et al. [26] leveraged the cross-entropy minimum optimization method in the random model of a random guardian command system, which can approximate the path distribution of the system by increasing the number of commands (number of parameters), to obtain an importance sampling distribution in the random model. However, the optimal importance sampling distribution obtained with the aforementioned methods is not from the distribution family of the system path space, but essentially is a heuristic importance sampling method. us, the verification results are only rough approximation.
In this paper, we propose a method with the SCCPS path space to construct a cross-entropy optimization model and use an iterative learning method to obtain an optimal importance sampling distribution from the parameterized distribution cluster of the path space. It can ensure that the optimal importance sampling distribution is from the spatial distribution family in the SCCPS path, and the iterative learning method can ensure that the distribution evenly covers the unsafe path distribution area. As evaluated in our experiments, the accuracy and efficiency of the rare attribute verification are significantly improved.

Statistical Model
Checking. Statistical Model Checking (SMC) can be simply described as follows: given a system model M and system properties φ described by the bounded linear temporal logic (BLTL) [18], it uses the Monte Carlo sampling, model checking, and statistical analysis techniques to qualitatively/quantitatively verify the following two questions: (i) e probability that M satisfies the attribute φ: Pr(M⊨π) (ii) Whether the probability of M satisfying the attribute φ is higher than or equal to the threshold θ: In SMC, it first simulates the execution of the system model M to extract a random execution path ω. en, the BLTL model detector is used to determine whether ω satisfies the attribute φ, and a certain number of samples will be generated after multiple simulations. It further leverages the statistical method to perform statistical analysis on the samples to assess the probability of the system model M satisfying the attribute φ, as well as give the confidence interval or the estimated error margin. Let I(ω) represent the output result of the BLTL model detector. If ω⊨π, I(ω) � 1; otherwise, it is 0. I(ω) is a Bernoulli random variable, so the behavior of M can be modeled by the Bernoulli distribution with a parameter p: (1) e parameter p represents the probability that the model M satisfies the BLTL attribute φ. With the Bernoulli distribution, we note that p � E[I(ω)], var[I(ω)] � p × (12212p). Since the value of p is unknown, the goal of SMC is to estimate the value of p.
SMC can be divided into two categories: hypothesis testing and parameter estimation. e hypothesis testing is used to determine whether the probability of the system satisfying the temporal logic attribute is greater than or equal to a given threshold, which is a qualitative result, while the parameter estimation is a quantitative result to represent the approximate probability of the system satisfying the temporal logic attribute. SMC qualitative algorithms include the single sampling plan (SSP) algorithm [27], the sequential probability ratio test (SPRT) algorithm [27], and the Bayesian hypothesis test (BHT) algorithm [18]. SMC quantitative algorithms mainly include the approximate probabilistic model checking (APMC) [28] algorithm and the Bayesian interval estimation testing (BIET) algorithm [18]. Kim et al. [29] conducted an empirical evaluation on the performance and applicability of the four algorithms (i.e., SSP, SPRT, BHT, and BIET).

Safety Requirement Specification.
In this paper, we use Bounded Linear Temporal Logic (BTCL) as our specification language. BLTL restricts Linear Temporal Logic (LTL) with time bounds on the temporal operators. Formally, the syntax of BLTL is given as (the finite set of state variables), v ∈ R, t ∈ R ≥0 , and ∨, ∧, and are the usual Boolean connectives. e formulas x ∼ v is called the atomic propositions (AP). e formula φ _1 { } is true and φ _1 { } will hold within the time t. e operators ◇ t and □ t can be defined as follows by using the ∪ t operator: ◇tφ � True ∪ tφ, which required φ to hold true within time t (true). □tφ �¬◇t¬φ requires φ to hold true up to time t. e semantics of BLTL formulas [28,30,31] is defined with respect to system traces (or executions). A trace is a sequence σ � (s 0 , t 0 ), (s 1 , t 1 ), . . ., where s i is the state of the system at the represented time t i . e pair (s i , t i ) expresses the fact that the system moved to state s i+1 after having spent t i time units in state s i . If the trace σ satisfies the property φ, we write σ⊨φ. e trace suffix of σ starting at k ∈ N is denoted by σ k , and σ 0 denotes the full trace σ. e semantics of BLTL for a trace σ k is defined as follows: Lemma 1 (Bounded sampling). e problem "σ⊨φ" is welldefined and can be checked for BLTL formulas φ and traces σ based on only a finite prefix of σ of bounded duration.
Proof. According to Lemma 1, the decision "σ⊨φ" is uniquely determined (and well-defined) by considering only a prefix of σ of duration #(φ)∈∈Q ≥ 0. By divergence of time, σ reaches or exceeds this duration #(φ) in some finite number of steps n. Let σ 0 denote a finite prefix of σ of length n, such that 0≤l<n tl ≥ #(φ). Again by Lemma 3, the semantics of σ 0 ⊨φ is well-defined because any extension σ" of σ' satisfies σ"⊨φ if and only if σ'⊨φ. Consequently, the semantics of σ'⊨φ coincides with the semantics of σ ⊨φ. On the finite trace σ 0 , it is easy to see that BLTL is decidable by evaluating the atomic formulas x ∼ v at each state s i of the system simulation. □ Lemma 2 (BLTL on bounded simulation traces). Let φ be a BLTL formula, k ∈ N. en, for any two infinite traces, σ � (s 0 , t 0 ), (s 1 , t 1 ), . . . and σ � (s 0 , t 0 ), (s 1 , t 1 ), . . . with s k+I � s k+I and t k+I � t k+I ∀I ∈ N with t k+I ≤ #(ϕ) [17]. We have that σ k ⊨φ if σ k ⊨φ.
Proof. IH is short for induction hypothesis.
(1) If φ is of the form x ∼ v, σ k ⊨φ if σ k ⊨φ since s k+I � s k+I and t k+I � t k+I by using [17] for i � 0. (2) If φ is of the form φ 1 ∨φ 2 , by induction hypothesis as #(φ 1 ∨φ 2 ) ≥ #(φ 1 ) and #(φ 1 ∨φ 2 ) ≥ #(φ 2 ). e proof is similar to φ 1 and lowing three conditions are satisfied: ) ≥ t such that the durations of trace σ and σ are t k+l � t ∼ k+l for each index l with 0 ≤ l < i by the assumption [17]. (b ′ ). σ k+i ⊨φ 2 by induction hypothesis as follows: we know that the traces σ and σ match at k for duration #(φ 1 ∪ t φ 2 ) and need to show that the semantics of φ 1 ∪ t φ 2 matches at k. By IH, we know that φ 2 has the same semantics at k + i (that is, k + i⊨φ 2 if k + i⊨φ 2 ) provided that we can show that the traces σ and σ match at k + i for duration #(φ 2 ). For this case, it considers any I ∈ N with 0≤l<I t k+i+l ≤ #(φ 2 ).
As I ∈ N was arbitrary, we conclude from this with assumption [17] that, indeed s I � s ∼ I and t I � t ∼ I for all I ∈ N with 0≤l<I t k+i+l ≤ #(φ 2 ). us, the IH for φ 2 yields the equivalence of σ k+i ⊨φ 2 and σ k+i ⊨φ 2 when using the equivalence of (a) and (a'). (c ′ ). For each 0 ≤ j < i, σ k+i ⊨φ 1 . e proof of equivalence to (c) is similar to that for (b') using j < i. e existence of an i ∈ N for which these conditions (a ′ ), (b ′ ), and (c ′ ) hold is equivalent to k⊨φ 1 ∪ t φ 2 . (ii) X⊆R n is a finite set of continuous variables; (iii) E ⊂ L × L is a collection of discrete changes; (iv) Inv: L ⟶ 2 X represents the mapping from the discrete state set L to the continuous state space. For ∀l ∈ L, Inv (l) is the invariant-position set of l; (v) D: L ⟶ (X ⟶ X) is a mapping of a vector domain, which assigns a set of Stochastic Differential Equations (SDE) to each control mode l ∈ L to describe the continuous random dynamic behavior with respect to the different control modes l, Security and Communication Networks standard Wiener process defined in the real number field. It assumes that ∀l ∈ L, f(l, ·), and g(l, ·) are bounded and Lipschitz continuous; (vi) G: E ⟶ 2 X is to assign a guardian condition to each discrete transition, satisfying the following conditions: represents a set of probability measures defined on X, and continuous variables are reset according to the probability distribution.
According to the definition, the SCCPS hybrid state space is L × X, and (l, x) ∈ L × X represents the hybrid state. e continuous dynamics of SCCPS evolves according to the SDE in the current control mode. However, the discrete dynamics refers to migrating one control mode to another control mode with the guardian condition on the discrete transition, when the continuous variable cannot reach the boundary of the invariant.
Let x l (t) be the SDE solution of the initial state x l (0); τ(l) � inf t ∈ R >0 , x l (t) ∉ Inv(l) means that, in the control mode l, the first time that the evolution of a continuous variable violates the invariant, that is, the first time of exiting the control mode l.
SCCPS execution semantics: a random execution of SCCPS is defined as a random process (l(t),

Our Approach
In this section, we present our proposed method with the SCCPS path space to construct a cross-entropy optimization model and use an iterative learning method to obtain an optimal importance sampling distribution from the parameterized distribution cluster of the path space.

Model Representation.
To avoid the complexity of the dynamic evolution of SCCPS, SMC does not pay attention to the structure of SCCPS, but focus on sampling the execution path of SCCPS. e behavior of SCCPS evolving over time can be characterized by the path of the system. According to the execution semantics of SCCPS, the execution path generation process of SCCPS can be described as follows: in the current control mode l i , the continuous variable x i evolves according to the SDE. When the evolution of x i satisfies the guardian condition (x i ∈ G(l i , l i+1 )), it migrates to the next control mode l i+1 and the initial value of x i+1 is determined by the random reset kernel R. e residence time of . t i is a random variable, and its value depends on the SDE of l i and the initial values x i (0) and Inv (l i ). According to the generation process of the SCCPS execution path, the next state of SCCPS depends on the current state and the related residence time of the current state. erefore, the execution path of the SCCPS can be regarded as that it is generated in the continuous-time Markov process in the hybrid state space. As the residence time of l i is longer, the probability of migration from l i is higher. It can further presume that the residence time of l i obeys the exponential distribution, and the continuous-time Markov process then becomes CTMC.
Let G l denote the guard condition set of all edges starting from l: where G(e) ∈z Inv (l) and G(e i ) ∩ G(e j ) � ∅, i ≠ j. In l, the time for the continuous variable evolving to satisfying the conditions of each guard is τ 1 , τ 2 , . . . , τ |G l | . en, the residence time in l is t l � min τ 1 , τ 2 , . . . , τ |G l | . Supposing τ 1 , τ 2 , . . . , τ |G l | , respectively, obey the exponential distribution of parameters λ l,l′ , l ′ ∈ L, (l, l ′ ) ∈ E , then the residence time t l in l obeys the exponential distribution of parameters l′∈Loc,(l,l′)∈E λ l,l′ . With this assumption, the execution path of SCCPS can be generated by the CTMC random process. It can be seen from this definition that when the CTMC structure is known, its behavior is controlled by the migration rate matrix λ, whose value comes from SCCPS. e value of λ is estimated with the maximum likelihood method according to simulating the execution of SCCPS to obtain the time samples of the state transition. 4 Security and Communication Networks

Algorithm of Learning Model Parameters.
e rarity of the path does not necessarily imply that the conversion rate between two adjacent discrete states is low, and the rarity of the safety attributes in the path space does not necessarily imply that the optimal parameters in the parameter space are rare. Based on this observation, this section introduces our approach of leveraging the maximum likelihood estimation method to estimate the migration rate of two adjacent discrete states of SCCPS and obtain the migration rate matrix λ. With the simulation operation of each discrete state of SCCPS, the discrete state is sampled to migrate to the next discrete state time; we then use the maximum likelihood estimation to obtain an estimate of λ.
For the state s i ∈ S, we simulate executing the SDE in the running state s i to obtain the migration time t k (k � 1, . . . , N) samples of the adjacent state s j . Assuming that the migration time between s i and s j obeys the exponential distribution of the parameter λ ij , then the likelihood function of λ ij can be obtained: and its log likelihood function is as follows: We further take the derivative of λ ij with the loglikelihood function and make it equal to 0, and its estimated value can be resolved, , it can be seen that the estimated value is an unbiased estimate of λ ij . e estimated variance is but the estimated variance is biased, and the variance will be decreased as the samples increase.
In most cases, it is difficult to obtain a clear expression for the random execution of SCCPS. However, what the safety concerned is the accessibility analysis of discrete states. e discrete state set S and its transitions can capture all necessary information. erefore, we derive the DTMC from the SCCPS path generation model to represent the path space of SCCPS. e value of DTMC's migration probability matrix P: S × S ⟶ [0, 1] can be obtained from the migration rate matrix λ of the SCCPS path generation model. For two states s i and s j ∈ S, where λ i � s j ∈S λ ij .

Method of Sampling Rare Attributes.
In the path space of the high-safety SCCPS, it is difficult to obtain samples satisfying the rare attributes, which makes the SMC infeasible. To address this challenge, we propose a method for sampling the rare attributes. It uses the cross-entropy method to learn an optimal-importance sample distribution from the path space of the SCCPS. With this sample distribution, it is easy to obtain the samples that satisfy the rare attributes. us, the convergence of the SMC can be accelerated. e importance sampling distribution is corrected by the likelihood ratio weighting to ensure that the SMC verification result is unbiased.

Zero-Variance Importance Sampling Distribution.
e basic idea of the importance sampling method [33,34] is to change the probability density distribution of random variables, so as to obtain the samples of extremely small probability events with a higher probability. We now present the SMC method based on the importance sampling. Let f(ω) be the true distribution of path ω, and let g(ω) be the importance sampling distribution, and g(ω) can obtain the samples of the extremely small probability events with a higher probability when g(ω) ≠ 0 and f(ω) ≠ 0. In the case of verifying the extremely small probability events, it is difficult to sample from f(ω) to meet the requirements, but the importance sampling method is to sample from g(ω). e probability p � E f [I(ω)] satisfying the system attribute can be described as where W(ω) � (f(ω)/g(ω)) is the likelihood ratio, and g(ω) is for the importance sampling. We leverage the likelihood ratio to correct the weighting to ensure that the estimated value of p is unbiased. We then randomly sample N independent execution paths ω i , i ∈ 1, . . . , N { } from the importance distribution g(ω) and obtain the unbiased estimate: and estimated variance for p, respectively. e efficiency and accuracy of importance sampling rely on the selection of the distribution g(ω). If the selection is inadequate, the importance sampling method is difficult to effectively achieve the acceleration effect and may play a decelerating effect. e key problem of importance sampling is to find a density function for the optimal sampling Security and Communication Networks 5 probability to minimize the estimated variance. With formula (10) returning 0, it can obtain the following formula: where g * (ω) is a zero-variance importance sampling distribution, which means that extracting only one sample from the zero-variance importance sampling distribution can be used to calculate its estimated value, that is, any sample is an unbiased estimate of its mean. However, the zero-variance importance sampling distribution depends on the true value p, and the value of p is unknown. erefore, it is impossible to sample from g * (ω). is paper proposes to use the crossentropy method to find an approximate optimal importance sampling distribution closest to g * (ω) from the parameterized distribution family of the sample path space, so as to reduce the SMC variance and accelerate the convergence of the SMC algorithm.

Cross-Entropy Optimization Model.
is section is to obtain the optimal importance sampling distribution by minimizing the cross entropy between the two probability distributions. According to the definition of cross entropy [35], this section provides the definition of cross entropy for the SCCPS path space.

Definition 2.
Cross entropy for the SCCPS path space: the cross entropy between two probability measures f(ω) and f ′ (ω) for the SCCPS path space Ω is as follows: e cross entropy is used to assess the similarity of two probability distributions. e value of cross entropy is smaller, and f(ω) and f ′ (ω) are more similar, i.e., CE(f(ω), f ′ (ω)) � 0 if and only if f(ω) � f ′ (ω).
According to Definition 2, the construction of the crossentropy optimization model on the SCCPS path space is given below. Assume that the original distribution f(ω) of the SCCPS path ω comes from the parameterized distribution family f(ω, θ) , e cross-entropy optimization method is used to select a distribution f(ω, λ * ), λ * ∈ θ in the parameterized distribution family, λ * ∈ θ and the optimal distribution g * (ω) have the smallest cross-entropy.
is optimization problem can be described for e first term of formula (13) has nothing to do with λ and minimizing cross entropy is equivalent to maximizing the second term. Let D(λ) � Ω g * (ω)Inf(ω, λ)d ω ; the minimization problem of formula (13) is equivalent to the maximization problem of formula (14): Solving the optimization problem of formula (14) requires sampling from the true distribution f(ω). However, in the case of rare attribute verification, it is difficult to sample from f(ω) to the path sample that satisfies the rare attribute. By using importance again, the sampling method samples from the distribution f(ω, μ) and the selection of parameter μ should be able to increase the probability of the path that meets the rare attribute. erefore, the optimization problem of formula (14) can be re-formed as Among them, the likelihood ratio function W(ω, μ) � (f(ω)/f(ω, μ)). In formula (16), the optimal solution of its optimization problem λ * can be estimated by the path sample, and the sample mean is replaced by the expectation Get the estimated value of λ * where ω 1 , ω 2 , . . . , ω N is a sample from the distribution f(ω, μ).

Algorithm of Verifying the Cross-Entropy Safety.
In Section 3.1, we provide a DTMC-based method to approximate the SCCPS path space. SMC mainly considers the system execution path ω � s 0 , s 1 , . . . , s k (k > 0) within a bounded time T, where k is a random variable to represent the number of state transitions, and its value varies with ω. Let 〈l, m〉 denote two adjacent and ordered state pairs in ω, S(ω) represent the set of ordered state pairs in ω, n (ω) lm represent the number of transitions from state l to state m in ω, and n (ω) l represent the number of occurrences of the state l in ω; then, the probability measure function of path ω under system parameter p can be formulated as Substituting f(ω i , λ) of formulas (16) with (17), we obtain and formula (18) can be transformed by the Lagrangian multiplier method into the following optimization problem: where ] i is the Lagrangian multiplier. Taking the derivative of formula (19) to p lm and making it equal to 0, the solution can be where ω i (1 ≤ i ≤ N) is the sample path from the distribution f(ω, μ), and f(ω i ) represents the true probability distribution of the SCCPS path.
With formula (20), it indicates that the estimated value of the optimal solution relies on the initial distribution f(ω, μ). However, the distribution of f(ω, μ) is generally far from the optimal distribution. erefore, in order to reduce the influence of the initial distribution f(ω, μ) on the optimal importance sampling distribution, this paper proposes the iterative solution in the path space. rough the iteration, the algorithm can explore a wider path space, so as to obtain a better approximate optimal solution.
Let the initial distribution parameter be u � p (0) , and an iterative formula can be obtained from formula (20): (22) where N is the number of samples per iteration, i , p (j) )) represents the likelihood ratio of the nth iteration, and ω Usually, only a few state transitions can be seen in each simulated execution. During each iteration, some parameters do not work in the path that satisfies the extremely small probability event. Formula (21) will set these parameter values to zero so that these parameters will not work in all subsequent iterations. As a result, the iterative algorithm converges too prematurely to detect a wider parameter space. To avoid this situation, this paper adopts a smoothing strategy to temporarily reduce the importance of inoperative parameters in the iteration instead of simply setting them to zero. e smoothing strategy is to weight current iteration value and the parameters of the previous iteration: e smoothing strategy can retain important but not yet effective parameters. Iterative formula (21) and smoothing formula (22) can jointly ensure that approximately uniform sampling is obtained from the path set of events satisfying the minimal probability.
e selected initial distribution f(·; p (0) ) should be able to produce some paths that satisfy the event with minimal probability in the first iteration, that is, the selected parameter p (0) should be able to increase the probability of occurrence of the extremely small probability events. erefore, in this paper, we set the initial parameter p (0) to a uniform distribution, and the uniform distribution can quickly obtain the sample path that satisfies the extremely small probability event. e condition for stopping the iteration can be that the coefficient of variance or the distance between two iteration parameter vectors are not higher than a certain constant or the maximum number of iterations. For example, given any small positive number ϵ > 0, if ‖p (j) − p(j − 1)‖ < ϵ is satisfied, the iteration will be stopped. To facilitate the comparison, we limit the maximum number of iterations in the experiment. To sum up, Algorithm 1 presents the description of the importance sampling distribution learning algorithm, which iteratively solves the approximate optimal importance sampling distribution in the SCCPS path space of the attributes for being verified.
Regardless of sample acquisition time and BLTL model checking time, the time complexity of Algorithm 1 is O(j max |p|N). Since the optimized objective function is convex, there is a unique optimal solution. If Algorithm 1 can converge, it must converge to the vicinity of the unique optimal solution [36]. Since the number of samples in each iteration is limited, the convergence is probabilistic but not necessarily monotonic. By simply limiting the maximum number of iterations j max , the algorithm can be guaranteed to be terminated with 100% probability. For the proof of convergence of cross-entropy optimization, please refer to [37]; thus, a formal proof of convergence is not provided in this paper. In experiments, we observe that the parameters Security and Communication Networks are convergent. Once the parameters converge, the last set of simulated samples is used to estimate the probability p that SCCPS satisfies the safety attribute with the optimal importance sampling distribution. Algorithm 2 describes the verification process of the safety verification algorithm.

Experiment and Analysis
To evaluate the effectiveness and performance of the Cross-Entropy Safety Verification Algorithm (CESVA) method proposed in this paper, we apply CESVA to a fault-tolerant controller for an aircraft elevator system (FTC4AE), that is, a Stateflow/Simulink hybrid system modeling case from MATLAB. It introduces the randomness in terms of the fault injection and simulates with MATLAB to obtain the system execution path. Path checking is realized by the BLTL model detector of Plasma-Lab [38]. In the experiment, the rare attributes of FTC4AE is verified with the CESVA method, which is further compared with the Heuristic Importance Sampling (HIS) method [17].

Validity Measurement of Experimental Results.
In the case of nonrare attribute verification, the confidence interval is used to assess the accuracy of various methods, while in the case of rare attribute verification, the relative error of sampling is used to assess the accuracy of the estimation: where E[p] is replaced by the current estimated value p, Skewness is a measure of assessing the skewing direction and degree of data distribution and is the characteristic number that characterizes the degree of asymmetry of the probability distribution density curve with respect to the average. Skewness is defined as the third-order standardized moment of the sample, and the skewness of the normal distribution is 0, and its estimator is evenly distributed around the mean: e negative skewness means that the distribution is lefttailed. At this time, the data on the left of the mean are less than the data on the right. Intuitively, the tail on the left is longer than the tail on the right. In contrast, the positive skewness means that the distribution is right-tailed. e data on the right of the mean is less than the left. Intuitively, the tail on the right is longer than the tail on the left.

Experiment and Analysis on a Fault-Tolerant Controller for the Aircraft Elevator System.
e fault-tolerant controller for an aircraft elevator system is a part of a large Simulink model of HL-20 rescuers developed by the National Aeronautics and Space Administration [39]. e two horizontal tails on the two side of the aircraft's fuselage are controlled by two elevators, respectively. Each elevator has two independent hydraulic actuators. In the normal operation process, each elevator is positioned by its corresponding external actuator, and its internal actuator can be used when the external actuator does not work. e two external actuators are driven by two independent hydraulic circuits, and the two internal actuators are both connected to the third hydraulic circuit. e system should ensure that only one set of actuators (i.e., external or internal) locates the elevator at any given time. If the external actuator or its corresponding hydraulic circuit fails, the system will activate the internal actuator. If the fault still exists, the external actuator will be shut down and eventually isolated. e fault in the hydraulic circuit may be temporary, and if the fault is cleared, the hydraulic circuit can always be restored to the online state. e control logic of the system is implemented in the form of a state flow diagram, while the hydraulic actuators and elevators are modeled by using Simulink.
According to modifying the Stateflow/Simulink model, we add random faults into three hydraulic circuits. Setting the fault model with an out-of-bounds' reading of circuit pressure, we model the fault injection as three independent Poisson processes. When the hydraulic circuit fails, the circuit will stay in the fault state for one second. en, the pressure reading will restore to its normal value, and the fault state will be terminated. In our experiments, the being estimated safety attribute is the probability that, within 25 seconds, the horizontal tails will not respond to the control inputs in the duration of 1 second.
We estimated the probability of the BLTL formula φ: where H 1 and H 3 represent the hydraulic circuit that drives the external actuator, while H 2 represents the hydraulic circuit that drives the internal actuator.
In the experiment, the failure rate of the three hydraulic circuits is set to 0.001, and the failure repair rate is 1. With the two parameters, the parameter v in Algorithm 1 can be calculated. It still is difficult to obtain samples that satisfy the attribute φ with the previous parameters. erefore, to ensure that the obtained samples can satisfy the attribute φ, the initial failure rate is set as 0.1 and the fault repair rate is set as 1. According to these two parameters, the initial parameter of iteration p (0) in Algorithm 1 can be calculated. In order to assess the performance of verifying the rare attributes with the CESVA method, 20 iterations of Algorithm 1 are performed. In each iteration, the number of samples is N � 104, the smoothing factor α � 0.2, and the total number of required samples is 2.0 × 10 5 . Figure 1 shows the change trend of the failure rate parameters during the 20 iterations of the CESVA method. At the beginning of the iteration, the parameters converge rapidly. When the parameters are close to their optimal values, the convergences of their values slow down with random fluctuations. From the 16th iteration, the failure rate parameters start to converge to the stable values. From the perspective of the parameter convergence trend, it seems that the value of the failure rate parameter increases with the increasing iteration times. It indicates that the proportion of sampling the paths satisfying the rare attribute is gradually increasing. Figure 2 illustrates the distribution of the estimated values of the CESVA method during the iterations. e estimated value gradually converges from the 17th iteration. Figure 3 presents the distribution of the relative error of the CESVA method during the iterations. e relative error gradually converges from the 16th iteration. Finally, the probability estimated value of the security attribute φ is 1.682 × 10 − 12 , and the value of the relative error is 0.01.
In order to verify the statistical performance of the CESVA method, 100 experiments were carried out under the above parameters, and 2.0 × 10 5 samples were used in each experiment. Compared with the performance of the HIS method under the same sample size, Table 1 shows the mean, skewness, and statistical indicators such as standard deviation (likelihood ratio standard deviation), relative error, and sample size for each experiment. As presented in Table 1, with the same sample size, the estimated values of the CESVA method are more closely distributed around the mean value, and the likelihood is over 10 times less than the standard deviation and relative error, when comparing against the HIS method. Although the true probability is unknown, statistical indicators such as the standard deviation, skewness, and relative error of the likelihood ratio illustrate that the true probability and the mean are very close.

Related Work
e verification of the rare attribute for SMC mainly includes the importance sampling method, the importance splitting method, and the statistical learning method. e importance sampling method is an effective method to solve the verification of rare attributes. For the CTMC and DTMC random models, Reijsbergen et al. [40] and Barbot et al. [23] leveraged the heuristic methods to obtain an importance sampling distribution to complete the attribute verification of the two types of models. For the Stateflow/Simulink model, Clarke and Zuliani [24] proposed the SMC method of cross-entropy minimization importance sampling to verify its safety properties. Zuliani et al. [17] further used the SMC method in paper [24] to verify the safety properties of a class of discretetime SHS. e method proposed by Clarke and Zuliani [24] assumes that the distribution of the system path space is exponential distribution. By simply increasing the failure rate of the system parameters and calculating the optimal parameters of the exponential distribution with the paths satisfying the rare attributes extracted at one time, an importance sampling distribution can be obtained. J´egourel et al. [26] used a random guardian command to the importance sampling distribution. is model can approximate the path distribution of the system by increasing the number of commands (the number of parameters) and uses the minimized cross-entropy method to obtain an importance sampling distribution in the random model. However, the optimal importance sampling distribution obtained by the above method does not come from the distribution family of the system path space, and these methods actually belong to the heuristic importance sampling method. e importance segmentation method [34] is a method of reducing the estimated variance. Based on the importance segmentation method, J´egourel et al. [33] proposed the SMC algorithm for the verification of small probability events.
e key idea is to decompose the system logic  attributes into embedded attributes, which makes its probability easier to be estimated and reduces the number of sample paths required by verification. To improve the performance, the attributes need to be decomposed into multiple levels with different probabilities. During the decomposition process, copying or eliminating paths depend on their intermediate behavior. When the decomposition is over, an estimated probability that the attribute is satisfied can be obtained. e importance segmentation method is essentially heuristic and depends on the model, but lacks the support of theoretical results.
Applying statistical learning methods to SMC is also an important research direction. Du et al. [19] proposed a learning SMC framework based on support vector machinebased two classifiers. It uses cost-sensitive and resampling methods to solve the unbalanced data learning problem of support vector machines and implements predicting and assessing the probability of occurrence of small-probability events with a relatively small number of samples. However, this method cannot obtain rare attribute samples. For the low-probability attributes of hardware circuits with multiple failure regions, Kumar et al. [41] assumed that the system failure distribution is a Gaussian mixture model, thus proposed to use the variational Bayes method to learn an optimal importance sampling distribution from the Gaussian mixture model. However, the optimal importance sampling distribution is not a distribution family from the system path space. Kalajdzic et al. [42] proposed an SMC method based on the principle of feedback control. is method learns a model of a cyber-physical fusion system by Input: N, the number of samples per iteration. Input: v, the true path distribution parameter of SCCPS. Input: p (0) , the initialization parameter. Input: j max , the maximum number of iterations. Output: p * Optimal parameters.
Input: N I S, e number of samples. Input: v, the true path distribution parameter of SCCPS. Input: p * , the optimal parameters calculated by Algorithm 1. Output: p, Probability of SCCPS meeting safety attributes.
(1) Function verifyingAlg (N, v, p (0) j max ) (2) A � 0, i � 1 (3) while i ≤ N do (4) generate a path ω i according to the pdf f(., p (j) ) (5) if ω i ⊨φ then (6) W i � 〈l,m〉∈S(ω i ) (v lm /p lm ) n (ω i ) lm : (7) A � A + W i ; (8) i � i + 1 (9) return (A/N IS ) ALGORITHM 2: Safety verification algorithm. using importance sampling to estimate the system state and importance division to control the system. So it can infer the probability that the system satisfies the given attributes. e method proposed in this paper starts from the SCCPS path probability space, constructs a cross-entropy optimization model, and uses an iterative learning method to obtain an optimal importance sampling distribution from the parameterized distribution clusters of the path space. It ensures that the optimal importance sampling distribution can come from the distribution family in the path probability space of SCCPS. And, the iterative learning method ensures that the distribution can evenly cover the unsafe path distribution area. erefore, the accuracy and efficiency of the rare attribute verification can be improved significantly.

Conclusion
SMC has been successfully applied to SCCPS safety attribute verification and has become the most effective solution, but rare attribute verification is still a challenge for SMC. To be able to extract samples satisfying the rare attributes from SCCPS, CTMC is used to construct the probability space model of the execution path of SCCPS given with the probability measure of the random execution path as well as the parameterized probability distribution function family, to construct the cross-entropy iterative model. According to the iteratively learning from finding the approximate optimal importance sampling distribution in the SCCPS path probability space, the efficient sampling of rare attribute samples in SCCPS is achieved. With the evaluating experiments, the experimental results show that, for the verification of rare attributes, comparing against the heuristic importance sampling method with the same number of samples, the estimated value of our method is better distributed around the mean, and the standard deviation and relative error are reduced by more than an order of magnitude. Based on the method proposed in this paper, combining with the current mainstream SMC method to develop an adaptive SMC tool is set as the future work.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request. e authors apply CESVA to a fault-tolerant controller for an aircraft elevator system (FTC4AE) that is a State-flow/ Simulink hybrid system modeling case from MATLAB.

Conflicts of Interest
e authors declare that they have no conflicts of interest.