Counterexample Generation for Probabilistic Model Checking Micro-Scale Cyber-Physical Systems

Micro-scale Cyber-Physical Systems (MCPSs) can be automatically and formally estimated by probabilistic model checking, on the level of system model MDPs (Markov Decision Processes) against desired requirements in PCTL (Probabilistic Computation Tree Logic). The counterexamples in probabilistic model checking are witnesses of requirements violation, which can provide the meaningful information for debugging, control, and synthesis of MCPSs. Solving the smallest counterexample for probabilistic model checking MDP has been proven to be an NPC (Non-deterministic Polynomial complete) problem. Although some heuristic methods are designed for this, it is usually difficult to fix the heuristic functions. In this paper, the Genetic algorithm optimized with heuristic, i.e., the heuristic Genetic algorithm, is firstly proposed to generate a counterexample for the probabilistic model checking MDP model of MCPSs. The diagnostic subgraph serves as a compact counterexample, and diagnostic paths of MDP constitute an AND/OR tree for constructing a diagnostic subgraph. Indirect path coding of the Genetic algorithm is used to extend the search range of the state space, and a heuristic crossover operator is used to generate more effective diagnostic paths. A prototype tool based on the probabilistic model checker PAT is developed, and some cases (dynamic power management and some communication protocols) are used to illustrate its feasibility and efficiency.


Introduction
Micro-scale Cyber-Physical Systems (MCPSs) are a special kind of CPS in micromachines, which are composed of micro/nanoscale components, and the integration of computation with physical processes in micro-/nano-assembly operations. MCPSs are about the intersections, not the union, of the physical and the cyber, and the behaviors of them are defined by both cyber and physical parts of the system [1][2][3]. Most of constituent elements of MCPSs or MCPSs themselves are usually accompanied with stochastic behaviors. The reasons for this can be classified as three aspects: (1) MCPSs contain the randomized algorithms, e.g., leader election algorithm, consensus algorithms; (2) unreliable and unpredictable system behaviors incurred by execution environment, e.g., message loss, processor failure; (3) performance evaluation by random variables assigned artificially, e.g., reliability, availability [4,5]. As an automated and complete formal verification technique at the level of system models, probabilistic model checking can be used to estimate whether the MCPSs satisfy a desired requirement property, or avoid an undesired outcome. As shown in Figure 1, MCPSs are modeled as DTMCs (Discrete-time Markov Chains), MDPs (Markov Decision Processes), PTA (Probabilistic tic Timed Automata), etc. [5]. The achieved requirement properties, e.g., function, reliability, robust, etc., are specified by PCTL (Probabilistic Computation Tree Logic), LTL (Liner Temporal logic) with probability bounds, This work belongs to the vertical dimension, which propose a counterexample generation method for probabilistic model checking MCPSs with nondeterministic and discrete-time stochastic behaviors. Up to now, existing probabilistic model checking tools cannot provide a counterexample directly. Counterexamples play a very important role in estimating or verifying MCPSs: (1) Counterexamples can feedback which parts of the system violate the requirements and provide diagnostic information; (2) Counterexamples are very effective in model-based testing, which can provide a reference for the design of test cases; (3) In the process of abstract refinement, counterexamples can provide guidance information for the refinement of rough abstract models [13,14]; (4) Counterexamples can be used to obtain the core of feasible plans in planning, such as task scheduling [15]; (5) Counterexamples are recently used to synthesize attacks for showing how confidentiality of systems can be broken, and the quality assurance of multi-agent systems [16].

Related Works
The counterexample in probabilistic model checking may be a set of diagnostic paths that satisfy a given property, and the cumulative probability mass does not satisfy a given bound. It cannot be generated during the process of probabilistic model checking and needs the dedicated algorithm [17]. We divide the existing works into two kinds of This work belongs to the vertical dimension, which propose a counterexample generation method for probabilistic model checking MCPSs with nondeterministic and discretetime stochastic behaviors. Up to now, existing probabilistic model checking tools cannot provide a counterexample directly. Counterexamples play a very important role in estimating or verifying MCPSs: (1) Counterexamples can feedback which parts of the system violate the requirements and provide diagnostic information; (2) Counterexamples are very effective in model-based testing, which can provide a reference for the design of test cases; (3) In the process of abstract refinement, counterexamples can provide guidance information for the refinement of rough abstract models [13,14]; (4) Counterexamples can be used to obtain the core of feasible plans in planning, such as task scheduling [15]; (5) Counterexamples are recently used to synthesize attacks for showing how confidentiality of systems can be broken, and the quality assurance of multi-agent systems [16].

Related Works
The counterexample in probabilistic model checking may be a set of diagnostic paths that satisfy a given property, and the cumulative probability mass does not satisfy a given bound. It cannot be generated during the process of probabilistic model checking and needs the dedicated algorithm [17]. We divide the existing works into two kinds of methods for generating a counterexample in probabilistic model checking: the accurate and approximate approach.

Accurate Approach
Han, Katoen, et al. [18] provide the theoretical and algorithmic foundations for counterexample generation in probabilistic model checking, which is the earliest research on the counterexample in probabilistic model checking. They define the counterexample for DTMC and translate the counterexample solving problem into the shortest path problem in the corresponding directed weight graph. Thus, the counterexample generation can explicitly enumerate the paths according to the probability they carry. To obtain the smallest counterexample, the enumeration stops when the cumulative probability of all generated paths exceeds the probability bound set in advance. Algorithmically, this can be done efficiently by translating this problem into a k-Shortest-Paths problem, where k is not fixed in advanced and is determined during the calculation. Eppstein and REA algorithms are applied to generate a counterexample for DTMC. On the basis of [19], the approach of calculating and representing a counterexample in a succinct way using a regular expression is presented in [20]. According to [20], the regular expressions that represent a counterexample can be calculated by the state elimination method of the automaton theory, which is guided by the k shortest paths search. The counterexample can also be compacted based on the strongly connected components (SCCs), and the obtained acyclic model is helpful in reducing efforts to determine the counterexample.
A hierarchical counterexample generated by performing SCC reduction is presented in [21]. Considering the unbounded-until property formula P ≤p ΦU ≤h Ψ , where h = ∞, counterexample generation can be carried out by applying k shortest paths algorithms such as Epstein algorithm, and the time complexity of which is O(a + blogb + c), where a is the number of states and b is the number of transitions of DTMC. For bound-until property P ≤p ΦU ≤h Ψ , where h ∈ N, a recursive enumeration algorithm to generate counterexample is presented in [20]. Reference [22] defines the counterexample of MDP as the general DTMC. Lal and Prabhakar [14] use a counterexample generated by this method to guide the abstraction-refinement for the polyhedral probabilistic hybrid systems. Taking PCTL path formula Φ 1 UΦ 2 as an example, the process of generating a counterexample for MDP by the Eppstein algorithm can be summarized as follows: (1) make Φ 2 -states and all ¬Φ 1 ∧ ¬Φ 2 -states absorbing, (2) insert a sink state and redirect all outgoing edges of Φ 2 -states to it, (3) turn it into a weighted digraph, (4) implicit representation of paths, (5) represent paths by a heap, (6) find the k shortest paths as the counterexample.

Approximate Approach
A critical subsystem, as a subsystem of DTMC, can represent a counterexample for probabilistic model checking. It is called a minimal critical subsystem if the number of states and transitions included is minimal compared to other subsystems, and the smallest subsystem if it is minimal and carries a maximal probability to reach a target state. Chadha et al. [23] further define the counterexample as the general MDP and give the corresponding algorithm for solving the minimal counterexample. This is the first work of counterexample generation for MDP.
Generating the smallest critical subsystems to represent the counterexample is an NP-complete problem which has been proved in [23,24], and it is hard to solve with exact and complete algorithms. A heuristic search method like best first search to obtain a counterexample is presented in [25,26]. However, it may not be a smallest, or even a minimal counterexample. Symbolic methods like bounded model checking (BMC) for finding a small critical subsystem have been developed in [25,26]; it uses a symbolic representation to reduce the size of the counterexample. Another option to determine the smallest critical subsystem is to use mixed integer linear programming techniques in [27,28]. Aljzzar et al. [29] propose a directed state space search method called XBF (XZ, XUZ) to generate a counterexample, which is the first work of counterexample generation with a heuristic. There are some works in this direction, such as [13,30,31], which optimizes the heuristics to generate a counterexample for DTMC. XBF extends the search strategy BF (Best-First) to generate a counterexample for MDP. In the framework of generating a counterexample for MDP by the Eppstein algorithm, XBF makes the following three modifications: (1) XBF records all parent states for each state; (2) XBF maintains a graph which contains, at any point in the search, the currently selected diagnostic subgraph; (3) when XBF finds the first target state, it continues the search for further target states until the whole state space has been processed, or termination is explicitly requested.

Our Contribution
This paper deals with counterexample generation for probabilistic model checking MCPSs with nondeterministic and discrete-time stochastic behaviors, in which MCPSs are modeled as MDP, and the requirements are specified as PCTL. Solving the smallest counterexample for probabilistic model checking MDP has been proven to be an NPC problem. Although some heuristic methods are designed to generate the counterexample for MDP, it is usually difficult to fix the global heuristic functions. Instead, we design the local heuristic function integrated into the Genetic algorithm (heuristic Genetic algorithm, HGA) to explore the state space, and propose a counterexample generation method for probabilistic model checking MDP against requirement property in PCTL. We use a diagnostic subgraph to represent the compact counterexample for MDP, and exploit diagnostic paths to constitute an AND/OR tree for constructing the diagnostic subgraph. We adopt the indirect path coding of the Genetic algorithm to extend the search range of state space, and the heuristic crossover operator to generate more effective diagnostic paths. A corresponding prototype tool is implemented based on PAT [16], and some MCPSs cases are used to illustrate its feasibility and efficiency.

Outline of the Paper
The rest of this paper is organized as follows. In Section 2, we introduce some preliminaries and definitions. Section 3 describes how to optimize the Genetic algorithm with a heuristic to generate a counterexample for MCPSs model MDPs. In Section 4, we show an experimental evaluation of some MCPSs cases, in order to illustrate the feasibility and efficiency. Conclusion and future work are in Section 5.

MDP
An MDP can be viewed as an extension of DTMC, which permits both probabilistic and non-deterministic choices. In an MDP, each transition includes a non-deterministic choice of actions in state s. Formally, an MDP is defined as follows.
Definition 1 (MDP). Let AP be a set of atomic propositions. A Markov decision process M is a tuple (S, s init , A, P, L), where S is a finite set of states, s init ∈ S is the initial state, A is a set of actions, P : S × A × S → [0, 1] is a probability transition function such that for every state, s ∈ S and an action α ∈ A: ∑ s ∈ S, P(s, α, s ) ∈ {0, 1}, and L : S → 2 AP is a labeling function of atomic propositions.
For a state s ∈ S, the probability of s to its successor state s by action a is given by P(s, α, s ). If and only if ∑ s ∈S P(s, α, s ) = 1, we call the action α enabled in the state s, otherwise, the action α is disabled. We use A(s) to represent a set of actions that are enabled at state s.
A path represents the execution of the system modeled by MDP; that is, the possible event behavior of the system, which can be described as an infinite sequence of states.

Definition 2 (Path). An (infinite) path of MDP M is an infinite sequence
We use len(σ) to denote the length of σ, which is determined by the number of states. For an infinite path σ, len(σ) is ∞. For a natural number i such that 0 ≤ i < len(σ), σ[i] refers to the ith state of σ, i.e., s i . Using σ (i) to indicate the ith prefix of σ, formally, → s i , which represents the i + 1 state of the path σ. For 0 ≤ i ≤ len(σ), A σ (i) denotes the ith action in σ, namely α i . For a finite path, last(σ) denotes the last state of σ. Paths M and Paths M f in represent the set of infinite paths and the set of all finite paths in M, respectively. Paths M (s) denotes the set of infinite paths which start at s and Paths M f in (s) denotes the set of finite paths starting from s. Non-determinism in an MDP is resolved by schedulers (also called adversary, policy, or strategy). A scheduler for an MDP M = (S, s init , A, P, L) is a function mapping every finite path σ in M onto an action d(σ) ∈ A(last(σ)). According to [32], we consider the deterministic scheduler that can deterministically select actions, which induces the maximal and minimal probability measures. The scheduler converts MDP into DTMC for which the probability of paths is measurable. We refer to the set of infinite paths under this schedule as Paths d (s 0 ). Paths in an MDP are called valid if the paths are allowed under a given scheduler d. We give the definition of the valid path as follows.
Otherwise, the path σ is invalid under scheduler d.
The underlying σ-algebra is formed by the cylinder sets which are induced by finite paths under the scheduler denoted FinitePaths d (s 0 ). The probability of this cylinder set is computed by using the following equation:

Probabilistic Computation Tree Logic
We consider the calculation of the (constrained) reachability probability in MDP. Let a finite MDP , , , , and ⊆ a set of a target. The maximal probability of reaching a state in B from the initial state of MDP is equivalent to determining Pr ⊨ ◊ Pr ⊨◊ . Note that the supremum ranges over all, potentially

Probabilistic Computation Tree Logic
We consider the calculation of the (constrained) reachability probability in MDP. Let a finite MDP M = (S, s init , A, P, L) and B ⊆ S a set of a target. The maximal probability of reaching a state in B from the initial state of MDP is equivalent to determining Pr max (s M ♦B) = sup D Pr(s ♦B). Note that the supremum ranges over all, potentially infinitely many, schedulers for M. Probabilistic computation tree logic (PCTL) is a probabilistic branching-time temporal logic. PCTL state formulae over the atomic propositions set AP are formed according to the following grammar: where a ∈ AP, φ is a path formula, ∼∈ {< , ≤, >, ≥}, and p ∈ [0, 1]. PCTL path formula is formed according to the following grammar: where Φ 1 and Φ 2 are state formulae, and n ∈ IN ≥0 U{∞}. The formula Φ 1 UΦ 2 means that Φ 2 is satisfied and all preceding states satisfy Φ 1 . The path formula Φ 1 U ≤n Φ 2 is the step-bounded variant of Φ 1 UΦ 2 . n = ∞ in a path formula means that it is unbounded until formula.
Let s ∈ S be a state and σ be an infinite path under D, where D denotes the set of all schedulers in MDP; the semantics of PCTL formulas over MDP are defined by satisfaction relation D : We abbreviate Pr D s ( σ ∈ Path D (s) σ D φ ) as Pr D s (φ) for convenience. It needs to calculate the probability of each path under D starting from s and satisfying φ, and determine which state satisfies the formula P ∼p (φ). We use the set Sat P ∼p (φ) to represent all states that satisfy P ∼p (φ). For MDP M = (S, s init , A, P, L), a formula P ∼p (φ) is satisfied if and only if for every d ∈ D : Pr d (φ) ∼ p, where Pr d (φ) denotes the probability of the set of all finite paths satisfying φ under scheduler d. The probability of paths in MDP is only defined in a given scheduler d. Probabilistic model checking MDP should consider the maximizing or minimizing probability values incurred by the different schedulers. Let Pr max (φ) denote the maximal probability mass at which a MDP M satisfies φ, Pr max (φ) = max(Pr d (φ)|d ∈ D) , and dually, the minimal probability Pr min (φ) = min(Pr d (φ) d ∈ D) . For properties in the upper bound, it is obvious that

Genetic Algorithm
The Genetic algorithm is a heuristic search, inspired by Charles Darwin's theory of natural evolution. This algorithm reflects the process of natural selection where the fittest individuals are selected for reproduction in order to produce offspring of the next generation. The Genetic algorithm has been applied to a board range of learning and optimization problems [33]. It begins with a set of a random population of coded candidate solutions which are called the chromosome. Then, the fitness is evaluated, measured by the fitness function of each candidate solution in the current population, and the fittest candidate solution is selected as parent of the next generation. The next step is to generate a second-generation population of solutions from those selected through crossover and mutation. The fittest parents and the new offspring form a new population. We give the standard Genetic algorithm in pseudocode, as shown in Algorithm 1.

Counterexample Generation with Heuristic Genetic Algorithm
In this section, we adopt a diagnostic subgraph to represent the counterexample, and optimize the Genetic algorithm with heuristic (heuristic Genetic algorithm, HGA) to generate a counterexample for probabilistic model checking MCPSs. The system model of MCPSs is MDP, and the requirement property is expressed by PCTL, i.e., P ∼p (φ). We will focus on the formula P ≤p (φ) in PCTL; the lower bounded formulas P >p (φ) and P ≥p (φ) can be converted to upper bound formulas.

Counterexample Represented by Diagnostic Subgraph
In classical model checking, the checked property will determine the shape of a counterexample. For linear temporal logic such as LTL, a failure path serves as a counterexample; and for CTL, the path as a counterexample is only for a subclass of global quantifier formulas. In probabilistic model checking, the situation is very complex. Let does not satisfy at state s, then the probability sum of all paths which satisfy formula φ with initial state s exceeds the probability bound p. It is obvious that the counterexample can be represented by finite paths. For a low-bounded property formula P ≥p (φ), the definition of the counterexample does not make sense, since the empty set satisfies this condition; diagnostic path Pr D s (φ) < P does not carry useful information. Therefore, in the following, we only consider the upper bound.
The PCTL property P ≤p (φ) is refused in an MDP, if there exists at least one scheduler d such that the probability mass of the paths FinitePaths d (s 0 ) satisfying φ under d exceeds p, where FinitePaths d (s 0 ) means a set of paths starting from state s 0 and satisfying φ under a scheduler d. We denote this set by FinitePaths d (s 0 ϕ). These finite paths are also called diagnostic path.

Definition 4 (Diagnostic paths).
Let M be an MDP model of a MCPS, and upper bounded formula P ≤p (φ) be a requirement property under consideration. A counterexample of the property P ≤p (φ) for MDP M is a set X of diagnostic paths in M such that Pr max (φ) > p holds. Figure 3 and the property P ≤0.4 (aUb). This property is violated in this model since there exists a scheduler that induces a set of diagnostic paths that satisfy formula aUb and their probability mass greater than 0.4. We have the following diagnostic paths:  We argue that a counterexample in probabilistic model checking should following three conditions: (1) explain why the probabilistic system does not property; (2) represent violation of a class of requirement properties; (3) identif the system simply and specifically. A set X of diagnostic paths can show the v a property, but it contains too much redundancy information. Thus, we diagnostic subgraph to represent a counterexample for MDP.

Example 2. Let us consider the example of MDP shown in
An MDP can be regarded as a weighted directed graph where the vertice represent states of MDP and the edges of graph represent transitions of MDP transition probability between states in the MDP can be converted into t between the vertices in the graph.  We argue that a counterexample in probabilistic model checking should satisfy the following three conditions: (1) explain why the probabilistic system does not satisfy the property; (2) represent violation of a class of requirement properties; (3) identify errors in the system simply and specifically. A set X of diagnostic paths can show the violation of a property, but it contains too much redundancy information. Thus, we define the diagnostic subgraph to represent a counterexample for MDP.
An MDP can be regarded as a weighted directed graph where the vertices of graph represent states of MDP and the edges of graph represent transitions of MDP. Thus, the transition probability between states in the MDP can be converted into the weight between the vertices in the graph.
The subgraph of MDP not only contains states and transitions of paths, but also all transitions connecting them in the original MDP.
Definition 6 (Diagnostic subgraph). Let M = (S, s init , A, P, L) be an MDP; the counterexample for s D P ≤p (φ) (s D P <p (φ)) is a diagnostic subgraph such that Pr D s (φ) > p (s D P ≥p (φ)) does not hold. Theorem 1. For MDPs, a Counterexample represented by the diagnostic sub-graph, has less than or equal to the state space of the counterexample represented by diagnostic paths.
Proof of Theorem 1. If a counterexample for MDP is one diagnostic path, it also can be represented with a diagnostic subgraph with this diagnostic path. If a counterexample for MDP is a measurable subset X = {X 1 , X 2 . . . X n ,} ⊆ {π ∈ FinitePaths d (s)|π D φ} such that Pr D s (φ) > p (s D P ≥p (φ)), each path in the measurable subset can be derived from a diagnostic loop path. We can get that the probability of a diagnostic loop path is greater than or equal to the diagnostic path. Thus, the state space of the counterexample represented by a diagnostic subgraph is less than or equal to the counterexample represented by the diagnostic paths.

Genetic Algorithm with Heuristic (HGA)
The most critical issue of counterexample generation with HGA is how to code paths in the weighted directed graph for MDP as a chromosome. In order to represent all possible paths in the graph, we adopt a priority-based coding method. Then, we use heuristic crossover operators to generate better individuals than the parents, which can improve the search speed. We calculate the shortest path by the fitness function.

Coding the Path
According to [34,35], there are two main coding schemas of the Genetic algorithm for solving the shortest path, namely, direct and indirect coding. Direct coding is based on the identification number of the node. The disadvantage of this coding is that an invalid path is generated because the random sequence of node identification numbers may not correspond to a valid path. Therefore, direct coding is not a good choice. We choose priority-based indirect coding, which is used in [35,36].
Indirect coding has significant advantages in generating finite paths compared to direct coding schemes. The indirect coding scheme used in the Genetic algorithm is shown in Figure 4. The initial node identification number appears directly on the path; it is different from the path that uses the instruction information about forming the path node. The guidance information is the priority of each node in a weighted directed graph. In the initialization phase, the priorities are randomly distributed. The path starting from the initial node and terminating at the target node is generated by a sequential node attach procedure. In each step of the path construction, if there are several successor nodes to consider, then the node with the highest path priority is selected. Repeat the above process until reaching the target node.
represented by the diagnostic paths. □

Genetic Algorithm with Heuristic (HGA)
The most critical issue of counterexample generation with HGA is how to cod in the weighted directed graph for MDP as a chromosome. In order to repre possible paths in the graph, we adopt a priority-based coding method. Then, heuristic crossover operators to generate better individuals than the parents, wh improve the search speed. We calculate the shortest path by the fitness function.

Coding the Path
According to [34,35], there are two main coding schemas of the Genetic algor solving the shortest path, namely, direct and indirect coding. Direct coding is b the identification number of the node. The disadvantage of this coding is that an path is generated because the random sequence of node identification numbers m correspond to a valid path. Therefore, direct coding is not a good choice. We priority-based indirect coding, which is used in [35,36].
Indirect coding has significant advantages in generating finite paths comp direct coding schemes. The indirect coding scheme used in the Genetic algorithm i in Figure 4. The initial node identification number appears directly on the pa different from the path that uses the instruction information about forming the pa The guidance information is the priority of each node in a weighted directed grap initialization phase, the priorities are randomly distributed. The path starting f initial node and terminating at the target node is generated by a sequential nod procedure. In each step of the path construction, if there are several successor n consider, then the node with the highest path priority is selected. Repeat the above until reaching the target node.

Fitness Function
The path quality is measured by the fitness function. In order to make the counterexample as small as possible, the goal of the fitness function is to minimize the search for the candidate path. The fitness function is described as follows: where C is the path set in the graph, that is, each chromosome in the Genetic algorithm. e is an edge of path C, and ω(e) is the weight of the edge.
The construction process of diagnostic paths is selecting individuals with higher fitness from the current population and eliminating individuals with lower fitness. Assuming that the population size is N and fitness of the ith chromosome is p i , then the probability of which the ith chromosome is selected is

Heuristic Crossover Operator
The crossover operation replaces the reorganization of the partial structure of two parent individuals to generate a new individual. This paper uses a heuristic crossover approach. For the cross of individuals c 1 and c 2 , we firstly generate a cross break-point K randomly, then copy the Kth to Nth gene string c 11 in c 1 to the back of c 2 and delete the same gene in c 2 as in c 11 ; the same process is performed to c 2 , the Kth to Nth gene string c 22 of c 2 is copied to the back of c 1 , and the same gene in c 11 as in gene string c 22 is deleted. This produces two legal intermediate individuals c 3 , c 4 then selects two individuals with high fitness as crossover offspring from c 1 , c 2 , c 3 , c 4 and puts the selected descendants into the mating pool for the next crossover.

Mutation Operator
The mutation operation is to change the gene values of certain gene spots on chromosome individuals in the population. Its purpose is to enable Genetic algorithms to have local random search capabilities and maintain group diversity. Suppose we perform mutation operations on individual A: [s a 1 a 2 a 3 a 4 a 5 t]; select a gene block Y: [a 2 a 3 a 4 a 5 ] on A first, and then randomly generate a path X: [a 2 R a 5 ] from a 2 to a 5 in graph, where R represents all genes on the path X, i.e., diagnostic paths. Then, individual A after the mutation operation is [s a 1 a 2 R a 5 t].

Generating Counterexample with HGA
The diagnostic paths are added in the set R, and we need to filter out the invalid paths in R. We call the set of all diagnostic paths from R which are valid under d the maximum of R. The goal is to provide a counterexample that contains only valid paths.

Definition 7 (Maximum of R). Let d be a maximizing scheduler of R, and the set of all diagnostic paths which are valid under d is called a maximum of R.
If R is a counterexample, then each maximum of R only contains the valid paths. Now, we need to compute a maximum X of R and to compute its probability Pr max (X). To this end, we define the compatible paths as follows.

Definition 8 (Compatible paths set). A set of paths C is called compatible, if and only if there exists a scheduler d such that all paths in C are valid under d.
We proposed this concept because each scheduler compatible subset X of R with maximal probability is a maximum of R. For diagnostic paths, the first prefix after the common prefix is decisive for scheduler compatibility. When a branch exists in a state, the diagnostic path is not compatible. We can always define a scheduler that allows a set of diagnostic paths to branch in an action. Therefore, the diagnostic path of the branch in the action is compatible. Checking the scheduler compatibility can be done as follows. R can be implemented as an AND/OR tree in order to check scheduler compatibility. Each diagnostic path is stored in R by mapping its states and actions to nodes. The OR nodes map to the state nodes, in which the scheduler makes a decision. The AND nodes map to the probability decisions after an action has been chosen by the scheduler. For a searched diagnostic path σ = s 0 → s k , we interpret σ as an alternating sequence of OR and AND nodes, i.e., S σ = (s 0 , α 0 . . . α k−1 , s k ). If R is empty, we add all nodes and edges in diagnostic path σ directly to R. If R is not empty, determine whether the longest prefix of S σ is already included in R, which starts from root s 0 , i.e., whether the longest prefix is shared by any path in S σ and R. When the longest prefix is included in the tree, insert diagnostic path σ into the tree. The remainder of diagnostic path σ becomes a new subtree, with the node ending with the longest prefix as root. Then, assign a probability to each node of this diagnostic path σ in a bottom-up manner. For AND root node in diagnostic path σ, its probability value is the sum of the values of its child nodes. The value marked on the leaf is the probability of the path from the root to it. The value marked in the internal AND node is the sum of values of its child nodes, and the value of the OR node is the maximum value of all child node values. Algorithm 2 illustrates the flows to a counterexample generation for MDP with HGA. Step1: Initialize R and X as an empty set, respectively, that is, without edges and nodes.
Step2: Generate a diagnostic path σ by HGA in Section 3.2 Step3: Replace diagnostic path σ with the sequence S σ = < s 0 , α 0 . . . s i α i > in the AND/OR tree Step4: If R is empty, then directly add the nodes and edges in the σ and go to step 8.
Step5: Determine whether the longest prefix in S σ already exists in R, starting from the initial node s 0 .
Step6: Add the remainder of S σ to the last node in the longest prefix.
Step7: Assign mark M p to each node in R from the bottom up, where M p is the probability of each node Step8: If a node n in S σ is an AND node, then its probability value M p is the sum of the child-nodes' value.
Step9: If a node n in S σ is an OR node, then its probability value M p is the maximum of the child-nodes' value.
Step11: Go to Step2 Theorem 2. The Algorithm 2 will terminate after generating a counterexample.
Proof of Theorem 2. Assume that Algorithm 2 does not terminate, only an infinite counterexample for s = D P ≤p (φ) is generated, i.e., the counterexample consists of an infinite path. Let C = {σ 1 , σ 2 . . .} be an infinite counterexample; we have Pr(σ i ) = a j . According to the characteristics of the limit, this means that: Let 0 < ε < L − p. According to (3), for some n ≥ N, |a n − L| < ε ⇒ |a n − L|< L − p ; but the finite set C = {σ 1 , σ 2 . . . σ n } is also a counterexample as Pr(C ) > p. This is a contradiction. Therefore, the algorithm will terminate after generating a counterexample.

An Example
We consider the MDP as shown in Figure 3 and the property P ≤0.5 (aUb). The AND/OR tree of it is shown in Figure 5. Searching the state space by HGA and generating the most probable diagnostic paths, then add it to the AND/OR tree, which is the implementation of R. s 2 is the AND node, the probability of which is the sum of its child-nodes' value, i.e., 0.35. s 4 is the OR node, the probability of which is the maximum of child-nodes' value, i.e., 0.25. The diagnostic paths we are interested in are marked with thick lines. The number identified next to each node is based on the probability calcu-lated in Algorithm 2, the probability of the root node is 0.6, which means the maximum probability of identification is 0.6. Since 0.6 > 0.5, the maximum that was obtained is a counterexample, as shown in Figure 5.
AND/OR tree of it is shown in Figure 5. Searching the state space by HGA and gen the most probable diagnostic paths, then add it to the AND/OR tree, which implementation of R.
is the AND node, the probability of which is the sum of it nodes' value, i.e., 0.35.
is the OR node, the probability of which is the maxim child-nodes' value, i.e., 0.25. The diagnostic paths we are interested in are marke thick lines. The number identified next to each node is based on the probability cal in Algorithm 2, the probability of the root node is 0.6, which means the ma probability of identification is 0.6. Since 0.6 > 0.5, the maximum that was obtain counterexample, as shown in Figure 5.

Experimentation
We develop a prototype tool CX-HGA for counterexample generation with based on PAT [7,16]. It is developed in the Java language with explicit-state data str (sparse matrices, bit-sets, etc.). MDPs are described with PAT language, a probability of the diagnostic subgraph is calculated by the PAT engine. All expe are performed on a machine with Intel Pentium CPU 3.2GHz speed and 1GB mem our algorithms, population size is set to 100 and the evolution generation i Crossover probability is 0.90 and mutation probability is 0.10. By some PRISM Ben cases, we compare HGA with the XBF-based (XZ and XUZ) algorithm and E algorithm (Eppstein) for counterexample generation, which are implemented in [37]. We also compare HGA and Genetic algorithm (GA) with directed coding. N DTMCs are a proper subset of MDPs. Any Markov chain is an MDP in which, for a s, A (s) is just a singleton set. Thus, the experiments also include two DTMC ca comparing with the existing works. We select 4 cases for demonstrating the effect

Experimentation
We develop a prototype tool CX-HGA for counterexample generation with HGA, based on PAT [7,16]. It is developed in the Java language with explicit-state data structures (sparse matrices, bit-sets, etc.). MDPs are described with PAT language, and the probability of the diagnostic subgraph is calculated by the PAT engine. All experiments are performed on a machine with Intel Pentium CPU 3.2 GHz speed and 1 GB memory. In our algorithms, population size is set to 100 and the evolution generation is 3000. Crossover probability is 0.90 and mutation probability is 0.10. By some PRISM Benchmark cases, we compare HGA with the XBF-based (XZ and XUZ) algorithm and Eppstein algorithm (Eppstein) for counterexample generation, which are implemented in DIPRO [37]. We also compare HGA and Genetic algorithm (GA) with directed coding. Note that DTMCs are a proper subset of MDPs. Any Markov chain is an MDP in which, for any state s, A (s) is just a singleton set. Thus, the experiments also include two DTMC cases for comparing with the existing works. We select 4 cases for demonstrating the effectiveness of the proposed method, which are either a typical MCPS or some important constituent protocols of MCPSs.

Dynamic Power Management
This case is the power dynamics management of the IBM disk drive [38,39]. It is a typical MCPS, and can be used for various areas, e.g., multi-robot collaboration. Its purpose is to minimize the power consumption while minimizing the impact on performance. It is modeled as a DTMC in PAT with approximately 140,000 states and transitions. We are interested in the event that is 100 or more requests getting lost within 400 milliseconds. The time consumption of transition in DTMC is 0.5 milliseconds, and the step number of property is 800. Therefore, the state formula corresponding to this event is true U ≤800 lost = 100 . The probability calculated by PAT is 0.014725297, if the timebound property of P <0.014725297 true U ≤800 lost = 100 is violated, which is a time-bound property of DTMC. We generate the counterexample for it by HGA, GA, XZ, XUZ, and Eppstein algorithms. As shown in Figure 6, it is obvious that the Eppstein algorithm is almost parallel to the y-axis for a counterexample probability close to zero. HGA, XZ, and XUZ can successfully provide counterexamples for all possible upper bounds with a reasonable computational effort. HGA finds the first increment of the counterexample after exploring a smaller fraction of the state space than XZ and XUZ. The Probability value of the first increment which is generated by HGA is greater than 0.011, and the reason for this is that we use a diagnostic subgraph to represent counterexample. After exploring about 2000 states and transitions, all algorithms can implement the second increment of the counterexample, and generate a counterexample after finding the second increment, because the total probability value exceeds the given probability. In the process of finding a counterexample, HGA explores fewer states and transitions than XZ and XUZ algorithms, except for the range between 0.006 and 0.011. HGA reaches the probability bound within about 3 s and consumes 250 KB memory. Meanwhile, GA, XZ, and XUZ algorithms reach the probability bound within about 3 s and 300 KB. GA generates a counterexample exploring more states and transitions, because it uses directed coding methods to generate many invalid paths and crossover operators used by GA cannot generate better offspring.
time-bound property of . ( = 100) is violated, which is a timebound property of DTMC. We generate the counterexample for it by HGA, GA, XZ, XUZ, and Eppstein algorithms. As shown in Figure 6, it is obvious that the Eppstein algorithm is almost parallel to the y-axis for a counterexample probability close to zero. HGA, XZ, and XUZ can successfully provide counterexamples for all possible upper bounds with a reasonable computational effort. HGA finds the first increment of the counterexample after exploring a smaller fraction of the state space than XZ and XUZ. The Probability value of the first increment which is generated by HGA is greater than 0.011, and the reason for this is that we use a diagnostic subgraph to represent counterexample. After exploring about 2000 states and transitions, all algorithms can implement the second increment of the counterexample, and generate a counterexample after finding the second increment, because the total probability value exceeds the given probability. In the process of finding a counterexample, HGA explores fewer states and transitions than XZ and XUZ algorithms, except for the range between 0.006 and 0.011. HGA reaches the probability bound within about 3s and consumes 250 KB memory. Meanwhile, GA, XZ, and XUZ algorithms reach the probability bound within about 3 s and 300 KB. GA generates a counterexample exploring more states and transitions, because it uses directed coding methods to generate many invalid paths and crossover operators used by GA cannot generate better offspring.

Synchronous Leader Election Protocol
This case is the synchronous leader election protocol [40,41]. It selects a unique leader node from N identical network nodes of the MCPSs. At each round, every node chooses a random number from {0, . . . , } as its ID. A node with the maximum unique ID will be elected as leader; otherwise, a new round begins between these nodes. The protocol system is modeled in PAT as a DTMC with 12,400 states and 16,495 transitions. Here, we set N and K as 4 and 8, respectively. This process will not stop until the leader is elected.

Synchronous Leader Election Protocol
This case is the synchronous leader election protocol [40,41]. It selects a unique leader node from N identical network nodes of the MCPSs. At each round, every node chooses a random number from {0, . . . , K} as its ID. A node with the maximum unique ID will be elected as leader; otherwise, a new round begins between these nodes. The protocol system is modeled in PAT as a DTMC with 12,400 states and 16,495 transitions. Here, we set N and K as 4 and 8, respectively. This process will not stop until the leader is elected. The probability of the formula (F <= (L * (N + 1) "elected")) is 0.95703125 by PAT, where "elected" indicates that a node has been selected as the leader, and L is set to be 1. We generate a counterexample of property P ≤0.95703125 (F <= (L * (N + 1) "elected")) by HGA, GA, XZ, XUZ, and Eppstein algorithms.
As shown in Figure 7, we can see that the Eppstein and XUZ algorithms grow almost parallel to the y-axis, so the probability of the counterexample is zero. HGA and XZ can generate a counterexample for probability upper bound by reasonable computational effort. The HGA finds the first increment after exploring about half of the entire state space. However, the probability of the first increment of counterexample is less than 0.5. After exploring all the states and transitions, the HGA and XZ algorithms can find the second increment. HGA explores fewer states and transitions than XZ; and HGA finds the first increment with consumption of about 150 s and 360 KB memory. The GA algorithm consumes too much memory to find a counterexample, the reason being that direct coding produces more invalid paths.
space. However, the probability of the first increment of counterexample is less than 0.5 After exploring all the states and transitions, the HGA and XZ algorithms can find the second increment. HGA explores fewer states and transitions than XZ; and HGA finds the first increment with consumption of about 150 s and 360 KB memory. The GA algorithm consumes too much memory to find a counterexample, the reason being that direct coding produces more invalid paths.

Zeroconf Protocol
This case is a Zeroconf protocol [42]. It offers a distributed "plug-and-play" solution in which the address configuration is managed by an individual MCPSs device when it i connected to the network. Each station has a single-message buffer and is cyclically attended by the server. The buffer could store the messages that it wants to send; in such cases, messages are not relevant after reconfiguring, and thus keeping these messages can slow down the network and make hosts reconfigure when they do not need to. We only consider version networks where the host does not do anything (No Reset) on these messages.
It is modeled as an MDP in PAT. The number of abstract hosts is denoted by N, the number of probes to send is denoted by K, and the probability of message loss is denoted by a loss. In this experiment, N and loss are set to 1000 and 0.1, respectively. K is set to 2 or 4. We are interested in the probability that a host picks an IP address already in use. In Table 1, "Time" and "Memory" represent the run time (in seconds) and memory (in KB) We compare Eppstein, XUZ, GA, and HGA on the two k-value variants of the model, a shown in Table 2 and Table 3. When the searched state space is very big, Eppstein and XUZ algorithms cannot get a counterexample. They are recorded as "Failed", which denotes out of memory. The time consumption of HGA is less than Eppstein, XUZ, and GA, at the expense of very little memory, especially when the searched state and transitions are very large.

Zeroconf Protocol
This case is a Zeroconf protocol [42]. It offers a distributed "plug-and-play" solution, in which the address configuration is managed by an individual MCPSs device when it is connected to the network. Each station has a single-message buffer and is cyclically attended by the server. The buffer could store the messages that it wants to send; in such cases, messages are not relevant after reconfiguring, and thus keeping these messages can slow down the network and make hosts reconfigure when they do not need to. We only consider version networks where the host does not do anything (No Reset) on these messages.
It is modeled as an MDP in PAT. The number of abstract hosts is denoted by N, the number of probes to send is denoted by K, and the probability of message loss is denoted by a loss. In this experiment, N and loss are set to 1000 and 0.1, respectively. K is set to 2 or 4. We are interested in the probability that a host picks an IP address already in use. In Table 1, "Time" and "Memory" represent the run time (in seconds) and memory (in KB). We compare Eppstein, XUZ, GA, and HGA on the two k-value variants of the model, as shown in Tables 2 and 3. When the searched state space is very big, Eppstein and XUZ algorithms cannot get a counterexample. They are recorded as "Failed", which denotes out of memory. The time consumption of HGA is less than Eppstein, XUZ, and GA, at the expense of very little memory, especially when the searched state and transitions are very large.

Bounded Retransmission Protocol
BRP (bounded retransmission protocol) is a variant of the alternating bit protocol [43] for the compositional MCPSs. It sends a file in several chunks, but allows only a bounded number of retransmissions of each chunk. It is modeled as an MDP in PAT, where N represents the number of blocks, and MAX represents the maximum number of retransmits allowed for each block. In the experiment, MAX is set to 5 and the size of state space is scaled by N. It is modeled as an MDP in PAT; we are interested in the probability that the sender does not report a successful transmission. Model information is shown in Table 4, and experimental results are shown in Tables 5 and 6, the fields of which have the same  meaning as Tables 1-3, respectively. The non-deterministic analysis of MDP results in 137,313 transitions at N = 640 and 685,153 transitions at N = 3200. As can be seen from Table 4, the transitions in DTMC are 98% of that in MDP, which indicates that the MDP model is low uncertainty. Tables 5 and 6 list very low probabilities because these methods fail to provide counterexamples for large probability boundaries in real time. This is caused by the very low probability of a single diagnostic path being carried in this model, which means that many diagnostic paths must be found and processed to increase the probability of counterexamples. The run time and memory consumption of HGA performs best, especially for large probability bounds. HGA has a greater advantage in memory consumption because more efficient paths are found with a heuristic, and for large probability boundaries, the number of diagnostic paths found is very large, so there is a great advantage in storing them in an AND/OR tree.

Analysis
Compared with the existing works, we can see that HGA can generate a counterexample effectively for all 4 cases, even if the Eppstein algorithm, XUZ algorithm, and Genetic algorithm cannot get the counterexample for large state space of MCPSs, which are noted as "failed" in the corresponding tables. The bigger the state space of MCPSs, the more the effectiveness of our method is apparent, compared with the related works. This is rooted in the advantages of HGA: firstly, under the action of the genetic operator, HGA has a strong search ability, which can find the global optimal solution of the counterexample with a large probability; secondly, inherent in the parallelism of HGA, it can effectively deal with the large state space of the MCPSs model MDP. Compared with GA, HGA also performs better in counterexample generation for MDP. It is because that we adopt the indirect path coding in HGA to extend the search range of state space, and the heuristic crossover operator to generate more effective diagnostic paths.

Conclusions
We propose an HGA-based counterexample generation for PCTL probabilistic model checking MCPSs with nondeterministic and discrete-time stochastic behaviors, i.e., MDP.
As far as we know, it is the first counterexample generation method that employs the Genetic algorithm to search the diagnostic subgraph of MDP. We encode paths into chromosomes to generate more efficient paths through priority-based indirect coding. The heuristic crossover operator guarantees that a number of populations are excellent individuals, which helps to improve the average value of the population and speed up the convergence. We use the diagnostic subgraph to represent the counterexample, which is more compact. Experimental results show that the time and memory consumption are less than existing works. In the future, we will continue to optimize HGA to generate counterexamples for MCPSs against requirements in PCTL* which is a super set of PCTL. At the same time, we would like to extend HGA to generate counterexamples for MCPSs models with continuous-time or hybrid-time stochastic behaviors, such as CTMDPs (Continuoustime Markov Decision Processes), PHA (Probabilistic Hybrid Automata), against CSL (Continuous-time Stochastic Logic) stochastic communication logics [44,45].