STAMINA: STochastic Approximate Model-Checker for INfinite-State Analysis

Stochastic model checking is a formal method that designers and engineers can use to determine the likelihood of safety and liveness properties. Checking properties using numerical model checking techniques requires enumerating the state space of the system to determine the probability that the system is in any given state at a desired time [17]. Real-world applications often have very large or even infinite state spaces.

Numerous state representation, reduction, and approximation methods have been proposed. Symbolic model checking based on multi-terminal binary decision diagrams (MTBDDs) [23] has achieved success in representing large Markov Decision Process (MDP) models with a few distinct probabilistic choices at each state, e.g., the shared coin protocol [3]. MTBDDs, however, are often inefficient for models with many different and distinct probability/rate values due to the inefficient representation of solution vectors. Continuous-time Markov chain (CTMC) models, whose state transition rate is a function of state variables, generally contain many distinct rate values. As a result, symbolic model checkers can run out of memory while verifying a typical CTMC model with as few as 73,000 states [23]. State reduction techniques, such as bisimulation minimization [7, 8, 14], abstraction [6, 12, 14, 20], symmetry reduction [5, 16], and partial order reduction [9] have been mainly extended to discrete-time, finite-state probabilistic systems. The three-valued abstraction [14] can reduce large, finite-state CTMCs. It may, however, provide inconclusive verification results due to abstraction.

To the best of our knowledge, only a few tools can analyze infinite-state probabilistic models, namely, STAR [19] and INFAMY [10]. The STAR tool primarily analyzes biochemical reaction networks. It approximates solutions to the chemical master equation (CME) using the method of conditional moments (MCM) [11] that combines moment-based and state-based representations of probability distributions. This hybrid approach represents species with low concentrations using a discrete stochastic description and numerically integrates a small master equation using the fourth order Runge-Kutta method over a small time interval [2]; and solves a system of conditional moment equations for higher concentration species, conditioned on the low concentration species. This method has been optimized to drop unlikely states and add likely states on-the-fly. STAR relies on a well-structured underlying Markov process with small sensitivity on the transient distribution. Also, it mainly reports state reachability probabilities, instead of checking a given probabilistic property. INFAMY is a truncation-based approach that explores the model’s state space up to a certain finite depth k. The truncated state space still grows exponentially with respect to exploration depth. Starting from the initial state, breadth-first state search is performed up to a certain finite depth. The error probability computed during the model checking depends on the depth of state exploration. Therefore, higher exploration depth generally incurs lower error probability.

This paper presents a new infinite-state stochastic model checker, STochastic Approximate Model-checker for INfinite-state Analysis (STAMINA). Our tool also takes a truncation-based approach. In particular, it maintains a probability estimate of each path being explored in the state space, and when the currently explored path probability drops below a specified threshold, it halts exploration of this path. All transitions exiting this state are redirected to an absorbing state. After all paths have been explored or truncated, transient Markov chain analysis is applied to determine the probability of a transient property of interest specified using Continuous Stochastic Logic (CSL) [4]. The calculated probability forms a lower bound on the probability, while the upper bound also includes the probability of the absorbing state. The actual probability of the CSL property is guaranteed to be within this range. An initial version of our tool and preliminary results are reported in [22]. Since that paper, our tool has been tightly integrated within the PRISM model checker [18] to improve performance, and we have also developed a new property-guided state expansion technique to expand the state space to tighten the reported probability range incrementally. This paper reports our results, which show significant improvement on both efficiency and verification accuracy over several non-trivial case studies from various application domains.

Figure 1 presents the architecture of STAMINA. Based on a user-specified probability threshold \(\varkappa \) (kappa), it first constructs a finite-state CTMC model \(\mathcal {C}\rfloor _{\varkappa }\) from the original infinite-state CTMC model \(\mathcal {C}\) using the state space approximation method presented in Sect. 2.1. \(\mathcal {C}\rfloor _{\varkappa }\) is then checked using the PRISM explicit-state model checker against a given CSL property \(\mathsf {P}_{\thicksim p} (\phi )\), where \(\thicksim \in \{<, \; >, \; \leqslant , \; \geqslant \}\) and \(p \in [0, 1]\) (for cases where it is desired that a predicate be true within a certain probability bound) or \(\mathsf {P}_{=?} (\phi )\) (for cases where it is desired that the exact probability of the predicate being true be calculated). Lower- and upper-bound probabilities that \(\phi \) holds, namely, \(P_{min}\) and \(P_{max}\), are then obtained, and their difference, i.e., \((P_{max} - P_{min})\), is the probability accumulated in the absorbing state \(\mathbf {x}_{abs}\) which abstracts all the states not included in the current state space. If \(p \in [P_{min}, P_{max}]\), it is not known whether \(\mathsf {P}_{\thicksim p} (\phi )\) holds. If exact probability is of interest and the probability range is larger than the user-defined precision \(\epsilon \), i.e., \((P_{max} - P_{min}) > \epsilon \), then the method does not give a meaningful result.

Architecture of STAMINA.

Full size image

For an inconclusive verification result from the previous step, STAMINA applies a property-guided approach, described in Sect. 2.2, to further expand \(\mathcal {C}\rfloor _{\varkappa }\), provided \(\mathsf {P}_{\thicksim p} (\phi )\) is a non-nested “until” formula; otherwise, it uses the previous method to expand the state space. Note that \(\varkappa \) also drops by the reduction factor \(\kappa _r\) to enable states that were previously ignored due to a low probability estimate to be included in the current state expansion. The expanded CTMC model \(\mathcal {C}\rfloor _{\varkappa }\) is then checked to obtain a new probability bound \([P_{min}, P_{max}]\). This iterative process repeats until one of the following conditions holds: (1) the target probability p falls outside the probability bound \([P_{min},P_{max}]\), (2) the probability bound is sufficiently small, i.e, \((P_{max} - P_{min}) < \epsilon \), or (3) a maximal number of iterations N has been reached (\(r \geqslant N\)).

2.1 State Space Approximation

The state space approximation method [22] truncates the state space based on a user-specified reachability threshold \(\varkappa \). During state exploration, the reachability-value function, \(\hat{\kappa }: \mathbf {X}{} \rightarrow \mathbb {R}^{+}\), estimates the probability of reaching a state on-the-fly, and is compared against \(\varkappa \) to determine whether the state search should terminate. Only states with a higher reachability-value than the reachability threshold are explored further.

State space approximation.

Full size image

Figure 2 illustrates the standard breadth first search (BFS) state exploration for reachability threshold \(\varkappa =0.25\). It starts from the initial state whose reachability-value i.e., \(\hat{\kappa }(\mathbf {x}_0)\), is initialized to 1.0 as shown in Fig. 2a. In the first step, two new states \(\mathbf {x}_1\) and \(\mathbf {x}_4\) are generated and their reachability-values are 0.8 and 0.2, respectively, as shown in Fig. 2b. The reachability-value in \(\mathbf {x}_0\) is distributed to its successor states, based on the probability of outgoing transitions from \(\mathbf {x}_0\) to its successor state. For the next step, only state \(\mathbf {x}_1\) is scheduled for exploration because \(\hat{\kappa }(\mathbf {x}_1) \ge \varkappa \). Note that the transition from \(\mathbf {x}_4\) to \(\mathbf {x}_0\) is executed because \(\mathbf {x}_0\) is already in the explored set. Expanding \(\mathbf {x}_1\) leads to two new states, namely \(\mathbf {x}_2\) and \(\mathbf {x}_5\) as shown in Fig. 2c, from which only \(\mathbf {x}_5\) is scheduled for further exploration. This leads to the generation of \(\mathbf {x}_6\) and \(\mathbf {x}_9\) shown in Fig. 2d. State exploration terminates after Fig. 2e since both newly generated states have reachability-values less than 0.25. States \(\mathbf {x}_2\), \(\mathbf {x}_4\), \(\mathbf {x}_6\) and \(\mathbf {x}_9\) are marked as terminal states. During state exploration, the reachability-value update is performed every time a new incoming path is added to a state because a new incoming path can add its contribution to the state, potentially bringing the reachability-value above \(\varkappa \), which in turn changes a terminal state to be non-terminal. When the truncated CTMC model \(\mathcal {C}\rfloor _{\varkappa }\) is analyzed, it introduces some error in the probability value of the property under verification, because of leakage the probability (i.e., cumulative path probabilities of reaching states not included in the explored state space) during the CTMC analysis. To account for probability loss, an abstract absorbing state \(\mathbf {x}_{abs}\) is created as the sole successor state for all terminal states on each truncated path. Figure 2e shows the addition of the absorbing state.

2.2 Property Based State Space Exploration

This paper introduces a property-guided state expansion method, in order to efficiently obtain a tightened probability bound. Since all non-nested CSL path formulas \(\phi \) (except those containing the “next” operator) derive from the “until” formula, \(\varPhi \; \mathcal {U}^I \; \varPsi \), construction of the set of terminal states for further expansion boils down to eliminating states that are known to satisfy or dissatisfy \(\varPhi \; \mathcal {U}\; \varPsi \). Given a state graph, a path starting from the initial state can never satisfy \(\varPhi \; \mathcal {U}\; \varPsi \), if it includes a state satisfying \(\lnot \varPhi \wedge \lnot \varPsi \). Also, if a path includes a state satisfying \(\varPsi \), satisfiability of \(\varPhi \; \mathcal {U}\; \varPsi \) can be determined without further expanding this path beyond the first \(\varPsi \)-state. Our property-guided state space expansion method identifies the path prefixes, from which satisfiability of \(\varPhi \; \mathcal {U}\; \varPsi \) can be determined, and shortens them by making the last state of each prefix absorbing based on the satisfiability of \((\lnot \varPhi \vee \varPsi )\). Only the non-absorbing states whose path probability is greater than the state probability estimate threshold \(\varkappa \) are expanded further. For detailed algorithms of STAMINA, readers are encouraged to read [21].

This section presents results on the following case studies to illustrate the accuracy and efficiency of STAMINA: a genetic toggle switch [20, 22]; the following examples from the PRISM benchmark suite [15]: grid world robot, cyclic server polling system, and tandem queuing network; and the Jackson queuing network from INFAMY case studies [1]. All case studies are evaluated on STAMINA and INFAMY, except the genetic toggle switch ^{Footnote 1}. Experiments are performed on a 3.2 GHz AMD Debian Linux PC with six cores and 64 GB of RAM. For all experiments, the maximal number of iterations N is set to 10, and the reduction factor \(\kappa _r\) is set to 1000. All experiments terminate due to \((P_{max} - P_{min}) < \epsilon \), where \(\epsilon =10^{-3}\), before they reach N. STAMINA is freely available at: https://github.com/formal-verification-research/stamina.

We compare the runtime, state size, and verification results between STAMINA and INFAMY using the same precision \(\epsilon =10^{-3}\). For all tables in this section, column \(\varkappa \) reports the probability estimate threshold used to terminate state generation in STAMINA. The state space size is listed in column \(|\mathcal {G}^{}|(K)\), where K indicates one thousand states. Column \(T(C{/}A)\) reports the state space construction (C) and analysis (A) time in seconds. For STAMINA, the total construction and analysis time is the cumulation of runtime for all \(\varkappa \) values for a model configuration. Columns \(P_{min}\) and \(P_{max}\) list the lower and upper probability bounds for the property under verification, and column P lists the single probability value (within the precision \(\epsilon \)) reported by INFAMY. We select the best runtime reported by three configurations of INFAMY. The improvement in state size (column \(|\mathcal {G}^{}|(X)\)) and runtime (column \(T(\%)\)) are represented by the ratio of state count generated by INFAMY to that of STAMINA (higher is better) and percentage improvement in runtime (higher is better), respectively.

Genetic Toggle Switch. The genetic toggle switch circuit model has two inputs, aTc and IPTG. It can be set to the OFF state by supplying it with aTc and can be set to the ON state by supplying it with IPTG [20]. Two important properties for a toggle switch circuit are the response time and the failure rate. The first experiments set IPTG to 100 to measure the toggle switch’s response time. It should be noted that the input value of 100 molecules of IPTG is chosen to ensure that the circuit switches to the ON state. The later experiments initialize IPTG to 0 to compute the failure rate, i.e., the probability that the circuit changes state erroneously within a cell cycle of 2, 100 s (an approximation of the cell cycle in E. coli [24]). Initially, LacI is set to 60 and TetR is set to 0 for both experiments. The CSL property used for both experiments, \(\mathsf {P}_{=?} \; [true \; \mathcal {U}^{\leqslant 2100} \; (TetR>40 \; \wedge \; LacI<20)]\), describes the probability of the circuit switching to the ON state within a cell cycle of 2, 100 s. The ON state is defined as LacI below 20 and TetR above 40 molecules.

The property-agnostic state space is generated with the probability estimate threshold \(\varkappa =10^{-3}\). Table 1 shows large probability bounds: [0, 0.999671] for \(\text {IPTG}=100\) and [0, 0.6975] for \(\text {IPTG} = 0\). It is obvious that they are significantly inaccurate w.r.t. the precision \(\epsilon \) of \(10^{-3}\). The \(\varkappa \) is then reduced to \(10^{-6}\) and state generation switches to the property-guided state expansion mode, where the CSL property is used to guide state exploration, based on the previous state graph. Each state expansion step reduces the \(\varkappa \) value by a factor of \(\kappa _r = 1000\). To measure the effectiveness of the property-guided state expansion approach, we compare state graphs generated with and without the property-guided state expansion, as indicated by the “property agnostic” and “property guided” rows in the table. Property-guided state expansion reduces the size of the state space without losing the analysis precision for the same value of \(\varkappa \). Specifically, the state expansion approach reduces the state space by almost \(20\%\) for the response rate experiment.

Robot World. This case study considers a robot moving in an n-by-n grid and a janitor moving in a larger grid Kn-by-Kn, where the constant K is used to significantly scale up the state space. The robot starts from the bottom left corner to reach the top right corner. The janitor moves around randomly. Either the robot or janitor can occupy one grid location at any given time. The robot also randomly communicates with the base station. The property of interest is the probability that the robot reaches the top right corner within 100 time units while periodically communicating with the base station, encoded as \(\mathsf {P}_{=?} \; [ \; (\mathsf {P}_{\geqslant 0.5} \; [ \; true \; \mathcal {U}^{\leqslant 7} \; communicate \;]) \; \mathcal {U}^{\leqslant 100} \; goal \;]\).

Table 2 provides a comparison of results for \(K=1024,64\) and \(n=64,32\). For smaller grid size i.e, 32-by-32, the robot can reach the goal with a high probability of \(97.56\%\). Where as for a larger value of \(n=64\) and \(K=64\), the robot is not able to reach the goal with considerable probability. STAMINA generates precise results that are similar to INFAMY, while exploring less than half of states with shorter runtime.

Jackson Queuing Network. A Jackson queuing network consists of N interconnected nodes (queues) with infinite queue capacity. Initially, all queues are considered empty. Each station is connected to a single server which distributes the arrived jobs to different stations. Customers arrive as a Poisson stream with intensity \(\lambda \) for N queues. The model is taken from [10, 13]. We compute the probability that, within 10 time units, the first queue has more that 3 jobs and the second queue has more than 5 jobs, given by \(\mathsf {P}_{=?} \; [ \; true \; \mathcal {U}^{\leqslant 10} \; (jobs\_1\geqslant 4 \; \wedge \; jobs\_2\geqslant 6)]\).

Table 2 summarizes the results for this model. STAMINA uses roughly equal time to construct and analyze the model for \(N=5\), whereas INFAMY takes significantly longer to construct the state space, making it slower in overall runtime. For \(N=4\), STAMINA is faster in generating verification results In both configurations, STAMINA only explores approximately one third of the states explored by INFAMY.

Cyclic Server Polling System. This case study is based on a cyclic server attending N stations. We consider the probability that station one is polled within 10 time units, \(\mathsf {P}_{=?} \; [ \; true \; \mathcal {U}^{\leqslant 10} \; station1\_polled \;]\). Table 2 summarizes the verification results for \(N=12,16,20\). The probability of station one being polled within 10 s is 1.0 for all configurations. Similar to previous case studies, STAMINA explores significantly smaller state space. The advantage of STAMINA in terms of runtime starts to manifest as the size of model (and hence the state space size) grows.

Tandem Queuing Network. A tandem queuing network is the simplest interconnected queuing network of two finite capacity (c) queues with one server each [18]. Customers join the first queue and enter the second queue immediately after completing the service. This paper considers the probability that the first queue becomes full in 0.25 time units, depicted by the CSL property \(\mathsf {P}_{=?} \; [ \; true \; \mathcal {U}^{\leqslant 0.25} \; queue1\_full \;]\).

As seen in Table 2, there is almost fifty percent probability that the first queue is full in 0.25 s irrespective of the queue capacity. As in the polling server, STAMINA explores significantly smaller state space. The runtime is similar for model with smaller queue capacity (\(c=2047\)). But the runtime improves as the queue capacity is increased.

This paper presents an infinite-state stochastic model checker, STAMINA, that uses path probability estimates to generate states with high probability and truncate unlikely states based on a specified threshold. Initial state construction is property agnostic, and the state space is used for stochastic model checking of a given CSL property. The calculated probability forms a lower and upper bound on the probability for the CSL property, which is guaranteed to include the actual probability. Next, if finer precision of the probability bound is required, it uses a property-guided state expansion technique to explore states to tighten the reported probability range incrementally. Implementation of STAMINA is built on top of the PRISM model checker with tight integration to its API. Performance and accuracy evaluation is performed on case studies taken from various application domains, and shows significant improvement over the state-of-art infinite-state stochastic model checker INFAMY. For future work, we plan to investigate methods to determine the reduction factor on-the-fly based on the probability bound. Another direction is to investigate heuristics to further improve the property-guided state expansion, as well as, techniques to dynamically remove unlikely states.