Simulation of Markovian models using Bootstrap method Ricardo M. Czekster, Paulo Fernandes, Afonso Sales, Dione Taschetto and Thais Webber Faculdade de Inform´atica – Pontif´ıcia Universidade Cat´olica do Rio Grande do Sul Av. Ipiranga, 6681, Pr´edio 32 – 90619-900 – Porto Alegre – RS – Brasil {ricardo.czekster, paulo.fernandes, afonso.sales, dione.taschetto, thais.webber}@pucrs.br Keywords: Markovian models, Discrete event simulation, Bootstrap method
Abstract Simulation is an interesting alternative to solve Markovian models. However, when compared to analytical and numerical solutions it suffers from a lack of precision in the results due to the very nature of simulation, which is the choice of samples through pseudorandom generation. This paper proposes a different way to simulate Markovian models by using a Bootstrap-based statistical method to minimize the effect of sample choices. The effectiveness of the proposed method, called Bootstrap simulation, is compared to the numerical solution results for a set of examples described using Stochastic Automata Networks modeling formalism.
1. INTRODUCTION Modeling formalisms are usually employed to describe real systems, capturing their behavior. An example of a wellknown modeling formalism used for such purpose is Markov Chains (MC) [20]. MC models use simple primitives such as states and labeled transitions to exemplify system’s evolution and operational semantics. Besides Computer Science abundant applications, Markov modeling examples are present in bioinformatics, economics, engineering and chemistry, to name a few domains [20]. However, the major limitation of MC is the state space explosion problem, observed when the amount of possible configurations is massively huge, rendering solution intractable and easily depleting computational resources [20]. In order to mitigate this problem, structured formalisms are proposed such as Stochastic Automata Networks (SAN) [9], Generalized Stochastic Petri Nets [6] or Stochastic Process Algebras (for example, Extended Markovian Process Algebra [3] or Performance Evaluation Process Algebra [12]), among several other approaches. Every structured formalism has a correspondent underlying MC and, in a particular way, SAN presents an inner representation based on tensor operations, producing a highly memory efficient fashion to map its underlying transition system. SAN also enables complex modeling due to sophisticated primitives (such as events) to define the firing of transitions, locally or synchronously, through constant or functional rates (rates dependent on the other component states).
Beyond model mapping, its numerical solution plays an equally important role when examining complex systems. There are distinct ways to solve any given model, e.g., using numerical methods or simulation. Numerical methods often rely on iterative mathematical techniques such as the Power Method [20], Arnoldi [1] or GMRES [16] to calculate measures of interest. On the other hand, simulation techniques consider the model events, firing transitions according to pseudorandom numbers generation. Simulation precision is strongly related to the number of samples that are produced. Normally, a simulation study comprehends large amounts of time in order to obtain significant performance indices, usually in function of the model state space size. Although numerical solution produces reliable results, it is bounded by the number of states of the model. Current numerical methods are directed to research on the acceleration of iterative methods [9]. However, methods based on Perfect Samplings for SAN have been used with satisfactory results, generating unbiased samples [10]. Another alternative is to use pure traditional simulation methods for Markovian models [11]. Such methods are based on simulating sets of trajectories and counting the amount of visits for each state, where every step produces a sample (saved for further analysis). When using traditional simulation methods one must be concerned about result quality, i.e., precision. Since every simulation trajectory only approximates the numerical solution, a thoroughly data inspection must be conducted, verifying results and searching for imprecise results. Succinctly, traditional simulation starts from an arbitrarily chosen initial state and randomly walks on a predefined state space during a fixed quantity of steps (trajectory length). Ideal length for a trajectory can be estimated using confidence intervals, among other approaches. Since every visited state generates a sample, dividing the number of times a state has been visited by the trajectory length gives an approximation of the mean permanence probability for each state. Bootstrap methods were firstly proposed by Efron [8] to perform estimations applied to different areas (for instance, machine learning algorithms [2]). The main idea is to reduce the “noise” that were observed from the samples. In the context of simulation, this noise corresponds to the noticed error in relation to the numerical solution. The method consists on resamplings, where each sample maps not only one state (the common approach in the traditional simulation) but a set of
states. In comparison, these sets of results, in average, are better than the average of a single sampling. Our aim in this work is to employ Bootstrap methods as a valid alternative to the solution of Markovian models by controlling the precision when producing large sets of samples. We focus our attention on proposing a method, called Bootstrap simulation, compared against a traditional simulation in order to measure the corresponding error of both methods. The remainder of this paper is presented as follows. We show basic modeling concepts based on Markov Chains and SAN, followed by principles of simulation in Section 2. We discuss Bootstrap method and how it can be applied in the context of Markovian simulation in Section 3. We present the main results and comparisons with traditional simulation in Section 4. Finally, the conclusion draws some final considerations and future works.
2. MARKOVIAN MODEL SIMULATION Markovian modeling are commonly used to describe complex system interactions using mathematical methods for solving the linear system of equations. The results produce quantitative performance indices used to subject models for analysis and interpretation. These performance indices normally signalize bottlenecks and configurations where the system will degrade or malfunction [7, 4, 18]. A common modeling formalism is Markov Chains [20], used for system modeling through simple primitives such as a set of states S and transitions among states. The transitions can be modeled by stochastic processes that vary according to a Continuous Time Markov Chains (CTMC) or Discrete Time Markov Chains (DTMC) mapping, respectively, exponential or geometric probability distributions. A DTMC is represented by a probability matrix P, where every row i and column j indicate the probability to change from state i to j, where the matrix dimension corresponds to the total number of model states S, i.e., the cardinality of S given by |S|. The matrix is mapped to a linear system of equations for solution purposes, i.e., the determination of the permanence probability of every state, solving πP = π. For CTMC, instead of a probability matrix, a transition matrix, called infinitesimal generator Q [20], represents the MC. In an infinitesimal generator all non-diagonal elements i, j represent the occurrence rate of a transition from state i to state j, and the diagonal elements are set with negative values to let the sum of every row equal to zero. In order to solve a CTMC it is common to discretize (or uniformize) the infinitesimal generator Q to its analogous DTMC version [20], obtaining a probability matrix P. When a MC model has a huge set of states S, usually, a larger amount of memory is needed to store the transition system Q and, consequently, P as well. This problem, known as state space explosion, is always on evidence in Marko-
vian models and more compositional approaches have been defined throughout the years to mitigate the hazardous effects present on massively sized models. Among other structured formalisms, Stochastic Automata Networks (SAN) [15] are based on MC and use a high-level description language to build modular systems in a memory efficient manner due mainly to a tensor representation of Q instead of the full transition matrix storage [9]. The main idea behind SAN is to model any reality with several subsystems, or automata, acting independently and, occasionally, synchronizing activities. Each automaton has a set of states with local transitions and/or synchronizing events having constant or functional rates associated. To compose the model state space, one must combine all local states from all automata, forming the model global states, i.e., the Product State Space (PSS). However, when performing this combination, depending on initial states and the model transitions, only a possible set of states can be valid (reached), defining thus the Reachable State Space (RSS) of the model. In order to inspect the performance indices for the models, Vector-Descriptor Product methods are usually used to multiply the probability vector π by the tensor representation of Q, named the Markovian Descriptor. The operation is currently bounded by the available memory needed to entirely store the probability vector. For example, in a 4GB machine, the maximum number for the PSS is set to ≈65 million of states. Recent efforts are based on solution approximations, particularly based on simulation of the underlying MC as a valid option to handle massive state spaces [13, 19]. Simulation methods, among many other advantages, are extremely useful to easily manipulate the components of Markovian models and use the probability vector to estimate reliable performance indices for huge state spaces. Obviously, the number of collected samples should be large enough for a comprehensive statistical analysis, e.g., the chosen amount should estimate according to the RSS size and the overall variance of the model in order to deliver the desired confidence interval. Due to the fact that simulation only produces samples belonging to the reachable set RSS, simulation approaches can potentially produce interesting results within memory bounds. In SAN, a probability vector of size |RSS| stores the number of occurrences of each reachable state within S = {s0 , s1 , ..., sr } and a stochastic matrix P maps transition probabilities Pi j for any given row i and column j. The simulation method needs a pseudorandom generation function U uniformly distributed in (0..1) [11]. Every transition between two states occurs according to a transition function given by φ(si ,U), dictating how the random walking takes place for the set of valid states defined by the SAN model. We are interested in defining an initial state and randomly visit the model states, generating a valid trajectory. The transition function φ is generically described by the following expression:
φ(si ,U) =
s0 s1 .. .
for U ∈ [0, Pi0 ) for U ∈ [Pi0 , Pi0 + Pi1 ) .. . ! " j
j−1
sj .. .
∑ Pil , ∑ Pil
for U ∈
l=0
l=0
r−1
#
(1)
.. .
sr
for U ∈
"
∑ Pir , 1 l=0
3. BOOTSTRAP SIMULATION Bootstrap method is relevant in statistics to derive estimations about standard deviations for complex distributions. In essence, this method works with previously unknown data ˜ and discovers the pertaining an infinite sized population Λ distribution in random samplings Λ of size n (extracted from ˜ [14]. Intuitively, the n samples this population, i.e., Λ ⊂ Λ) obtained from the observations of the infinite population Λ, each one with probability 1n , are used to model the unknown population as real as possible. The set of values of a new sample K, also having size equal to n (|K| = |Λ|) is exclusively obtained from the first sam˜ ple Λ without generating new values for the population Λ, as showed by Figure 1. The main feature of the Bootstrap method is to sample with replacement, i.e., the values that constitutes the sample could be repeated during the process. To exemplify this method, lets consider the problem of finding the average height of the whole world population. Since it is impractical to measure every individual within the available ˜ only a subset of this population should be population (set Λ), taken (subset Λ). Using Λ, a quantity of z resamplings must be performed, corresponding to the number of bootstraps as demonstrated by Figure 1, where Ki = {x ∈ Ki | x ∈ Λ} and i ∈ [1..z]. Each bootstrap Ki has a set of n values obtained from Λ through pseudorandom samplings that could repeat or not the values (i.e., x1 ∈ Ki , x2 ∈ Ki | x1 = x2 ). The mean x¯K , or the average population height, is calculated using the averages x¯K1 , x¯K2 , . . . , x¯Kz of every bootstrap. ˜ Λ
Λ
x¯Λ
The objective of Bootstrap method is to generate more precise values for the average x¯K than the average x¯Λ . Moreover, as the number of bootstraps is increased, it reduces the effects of bias samples (noise), generally associated with pseudorandom generators. The adaptation of the method to be used in Markovian simulation concerns the transition function φ. As mentioned before, a uniformly distributed sample is generated between 0 and 1 and the transition function φ(si ,U) is applied to determine the next state to be visited, forming a simulation trajectory that maps the Λ sample. Figure 2 illustrates the Bootstrap simulation schema, where each visited state produces pseudorandom samplings that define the states accounted in every bootstrap K. The number of trials is equal to the amount of present values in Λ, which is equivalent to the trajectory length (n). In this figure, Tb (where b ∈ [1..z]) represents n samplings for each bootstrap K. Thus, n × z samplings take place for every trajectory step, comparing each pseudorandom value against an arbitrarily chosen constant value α between 0 and n − 1. If the pseudorandom generated value is equal to α, the state is accounted in the correspondent bootstrap K. This procedure is repeated for all bootstraps until a predefined trajectory length n is reached. ¯ that represents this execution. When Figure 2 also shows |b| the trajectory is finished, the probability vector π contains the permanence probabilities of model states. However, it is computationally too expensive to randomly generate n samples for every trajectory step. Since the number of times each step will actually be randomly sampled is very likely to be much smaller than n, we propose, instead of generating n samples for each bootstrap, only execute n¯ trials, where n¯ ≪ n. To evaluate the impact of reducing this number of trials, from n to n, ¯ we calculated the probability of a same value to be sampled more than n¯ times during n trials. If this probability is negligible, it means that the number of samples n¯ is sufficient. To compute such a value, we observe that: " # Probability of a given n¯ n − 1 n−i 1 i number to appear at = ∑ × (2) n n i=0 the most n¯ times Hence, simplifying (2) and considering its complement, i.e., the probability of having more than n¯ becomes: n¯ Probability of a given ∑i=0 (n − 1)n−i number appear more = 1 − nn than n¯ times
K1
K2
x¯K1
x¯K2
...
x¯K
Figure 1. Bootstrap method illustration
Kz x¯Kz
(3)
For example, assuming n¯ = 15 and n = 104 , the probability to obtain a same value more than 15 times with n samples is inferior to 10−16. Algorithm 1 presents the application of the Bootstrap in the context of Markovian simulation. In this algorithm, variables and vectors used in the method are initialized between lines
States ¯ |b| s2
s2
s2
s2
s2
s2
¯ |b|
¯ |b| s1
s1
s1
s1
s1
s0
s0
s0
s0
1
2
3
4
¯ |b| Initial state
s0
...
s1
...
s0
...
¯ |b|
Time 0 U = 0.08
U = 0.87
K1
U = 0.32
K2
s0
s0
s1
s1
s2
s2
U = 0.06
...
where:
s0
s1
s1
s2
s2
Transition Matrix
s1 s2
x¯2
s0
n
Kz s0
P= x¯1
5 ... ... U = 0.56
x¯z
s0
s1
s2
s0
0.10
0.65
0.25
s1
0.25
0.55
0.20
s2
0.30
0.25
0.45
s0 ...
s1
and T1
s2
U1 (0..n¯ − 1) .. . π π0
x¯1 [0] + x¯2 [0] + ··· + x¯z [0] z
π1
x¯1 [1] + x¯2 [1] + ··· + x¯z [1] z x¯1 [2] + x¯2 [2] + ··· + x¯z [2] z
π2
¯ = |b|
Un¯ (0..n¯ − 1) .. . Tz U1 (0..n¯ − 1) ... Un¯ (0..n¯ − 1)
Figure 2. Bootstrap simulation schema 1 and 4. The samples counted in the bootstraps for all trajectories of length n are shown from line 6 to line 16. Between lines 17 and 25, it is calculated the permanence probabilities of states for each bootstrap. For the rest part of the algorithm (lines 26-31), the average probabilities of every state in vector π are adjusted according the information contained in the bootstraps.
4. NUMERICAL ANALYSIS The numerical analysis of Bootstrap simulation method is conducted using three SAN models with different characteristics [17, 7, 4], such as: Alternate Service Patterns (ASP), First Available Server (FAS) and Resource Sharing (RS). ASP model describes an Open Queueing Network [20] with servers that map different service patterns. It has four queues represented by four automata, where its PSS is equal to (K + 1)4 × P, K is the queues’ capacity and P the number of service patterns. All states of this model are reachable. FAS model indicates the availability of N servers, where every server is composed of a two state automaton, represent-
ing the two possible server conditions: available or busy. In the model, tasks are firstly assigned to the first server, case it is available. If the server is busy, the task must be sent to the second server and so on. The first available server processes the task. All states in this model are reachable and its PSS is equal to 2N . The classical RS model maps R shared resources to P processes. Each process is represented by an automaton with two states: idle or busy. The number of available resources is represented by an automaton that counts the number of resources being used. The PSS for this model is equal to 2P × (R + 1) and, due to the nature of the model transitions, not all states are reachable. The parameters for these models adopted in this paper experiments are the following: ASP model - every queue with capacity two (K = 2) and two service patterns (P = 2); FAS model - with nine servers (N = 9); and RS model - with 10 processes (P = 10) and five resources (R = 5). Algorithm 1 Bootstrap simulation method 1: α ← U(0..n¯ − 1) { constant value α initialization with a pseudorandom value chosen between 0 and n¯ − 1 } 2: π ← 0 { initialization of the probability vector π } 3: K ← 0 { initialization of all z bootstraps K } 4: sc ← s0 { set sc as initial state s0 } 5: { walk on a trajectory of length n } 6: for t = 1 to n do 7: sd ← φ(sc ,U(0..1)) { finds destination state sd from sc according to U(0..1) } 8: for b = 1 to z do 9: for c = 1 to n¯ do 10: if (U(0..n¯ − 1) == α) then 11: Kb [sd ] ← Kb [sd ] + 1 { counts in Kb [sd ] every time the sample equals to α } 12: end if 13: end for 14: end for 15: sc ← sd { current state sc is updated to sd } 16: end for 17: for b = 1 to z do 18: ω←0 19: for i = 1 to |RSS| do 20: ω ← ω + Kb [i] { calculates in ω the total sum of accumulated values in Kb } 21: end for 22: for i = 1 to |RSS| do K [i]
23: x¯b [i] ← ωb { calculates the probability of i-th state in Kb } 24: end for 25: end for 26: for i = 1 to |RSS| do 27: for b = 1 to z do 28: π[i] ← π[i] + x¯b [i] 29: end for π[i] 30: π[i] ← z {calculates the average probability from the bootstraps} 31: end for
(a) ASP model .9
47
5
Mean (Trad.) Mean (Boot.) Maximum (Trad.) Maximum (Boot.)
96
14
54 105
106
05 08
08
0.
0. 70 01 108
0.0009
0.
00
51
0. 107
0.0028 0.0061
0.001
104
40
0 0 . .1 0 05 3 73 4 0.0082 0.0069
0.0902 0.0285
0.01
0.0272 0.0110
0. 0.3740
0.1
0.2608 0.1017
Relative error
1
0. 0. 14 20 38 38
67
45
1.
2.
91
93
10
109
Trajectory length (n)
(b) FAS model
21.7341
50.0041
7.2506
10
11.7229 1.6371
3.6133 0.4680 0.1552
0.1809
105
106
107
0.1024
108
0.0146
0.01
104
0.3522 0.2895
0.1429 0.0582
0.1
1.2137
0.0485 0.0383
0.8672 0.4214
1.5182 0.9986
1
1.2876
Relative error
Mean (Trad.) Mean (Boot.) Maximum (Trad.) Maximum (Boot.)
101.2809
100
last tested trajectory (109), the mean relative error is 0.0061 in traditional simulation and plummets to 0.0009 when the Bootstrap simulation is employed. Comparing the results for FAS model - Figure 3 (b) - we also observed a significant reduction in the mean relative error for Bootstrap when the trajectory length exceeds 106 . RS model - Figure 3 (c) - presents similar results for Bootstrap and traditional simulation (for instance, when the trajectory length is equal to 109 , the mean relative error in traditional simulation is 0.0010, whereas in Bootstrap the error is 0.0011). However, even not presenting significant precision gains, the mean relative error is very low. In simulation, the quantity of samples is a direct indication of the solution approximation quality, i.e., the precision. In the context of this paper, the measure of quality used is the simulation mean relative error that decreases, sometimes dramatically, as the trajectory length is increased. Figure 3 shows that using the Bootstrap simulation is possible to generate samples that greatly approximate the numerical solution, when compared to the traditional simulation approaches, considering trajectory lengths greater than 107 . For the worst case (RS model) presented in this work, we obtained the same precision of traditional simulation, which is a clear indication that our proposed method could be applied to greater models without significantly impairing precision.
109
Trajectory length (n)
5. CONCLUSION
(c) RS model
12
105
106
107
108
0.0011
0.0034 0.0010
0.
0. 006 00 5 51
02 0.
0.0108 0.0033
0.0135
0.001
104
06
0. 0.0466
0.0307 0.0110
0.01
06
0. 0. 160 17 4 07
0.1071 0.0325
0.3487 0.1063
0.1
0.3151
Relative error
1
0. 0. 652 59 9 81
2.
2. 408 14 0 59
10
Mean (Trad.) Mean (Boot.) Maximum (Trad.) Maximum (Boot.)
109
Trajectory length (n)
Figure 3. Traditional and Bootstrap simulation results The results obtained from traditional simulation are compared to the Bootstrap simulation, where the mean relative error is used to compare the approximate solution to the numerical solution produced by the SAN models solver [5]. The main simulation results are presented in Figure 3 showing 95% confidence intervals for 50 executions, where the x-axis corresponds to the trajectory length and the y-axis indicates the relative error. The dotted lines represent the maximum relative error and the bars indicate the mean relative error. ASP model - Figure 3 (a) - presents the best precision results for the Bootstrap simulation, since the mean relative error is reduced for the set of tested trajectory lengths. For the
The main contribution of this paper is the proposition of the Bootstrap method usage adapted for a Markovian simulation context. Research on such algorithmic adaptation is extremely important when analytical modeling fails to produce relevant performance indices. Simulation, for these cases, emerges as a valid alternative to deliver analysis of huge models, allowing result interpretations and helping decision making process as a whole. Our aim is to present the Bootstrap simulation and inspect the mean relative error when compared to traditional simulation. Our results produce evidence that Bootstrap simulation method has the potential to reduce the observed noise related to sampling from large state spaces. This fact was detected for trajectory lengths greater than 106 for the available models. However, more experiments must be done to validate such findings and report the classes of models that are more suitable for Bootstrap utilization. In regards of the traditional simulation, trajectory length plays a relevant role in terms of the associated error. We observed that small trajectory lengths induce high error levels, which is interesting only for executions requiring few precision. Future works are directed to the following studies: (i) models classification - determination of classes of models more suitable for Bootstrap simulation; (ii) method parallelization and distribution - we shall take advantage from the fact that every sample is independent, so the computation can be con-
ducted separately for each model, improving the overall quality of samples for further analysis; and (iii) rare events - we plan to study the impact of Bootstrap in concealing (or not) rare events and its influence on precision of results. The work in progress herein described clearly indicates that further research efforts must be done to increase the understanding in regard to state-based systems simulation applied with different statistical methods. However, the preliminary results reveal that Bootstrap simulation is a worthy subject of research.
ACKNOWLEDGMENTS This work is funded by Petrobras (0050.0048664.09.9). Afonso Sales receives grants from CAPES-Brazil (02388/090). Paulo Fernandes is also funded by CNPq-Brazil (307272/2007-9). The authors thank Prof. Jo˜ao Batista de Oliveira for the helpful discussion on this paper subject.
REFERENCES [1] W. E. Arnoldi. The principle of minimized iterations in the solution of the matrix eigenvalue problem. Quarterly of Applied Mathematics, 9:17–29, 1951.
on Modeling and Optimization in Mobile, Ad Hoc, and Wireless Networks (WiOpt’05), pages 164–173, Trentino, Italy, April 2005. [8] B. Efron. Bootstrap Methods: Another Look at the Jackknife. The Annals of Statistics, 7(1):1–26, 1979. [9] P. Fernandes, B. Plateau, and W. J. Stewart. Efficient descriptor-vector multiplication in Stochastic Automata Networks. Journal of the ACM (JACM), 45(3):381–414, 1998. [10] P. Fernandes, J. M. Vincent, and T. Webber. Perfect Simulation of Stochastic Automata Networks. In International Conference on Analytical and Stochastic Modelling Techniques and Applications (ASMTA’08), volume 5055 of LNCS, pages 249–263, 2008. [11] O. H¨aggstr¨om. Finite Markov Chains and Algorithmic Applications. Cambridge University Press, 2002. [12] J. Hillston. A compositional approach to performance modelling. Cambridge University Press, USA, 1996. [13] A. M. Law and D. W. Kelton. Simulation Modeling and Analysis. McGraw-Hill Higher Education, 2000.
[2] E. Bauer and R. Kohavi. An Empirical Comparison of Voting Classification Algorithms: Bagging, Boosting, and Variants. Machine Learning, 36(1-2):105–139, 1999.
[14] B. F. J. Manly. Randomization, Bootstrap and Monte Carlo Methods in Biology. Chapman & Hall/CRC, second edition, 1997.
[3] M. Bernardo and R. Gorrieri. A tutorial on EMPA: a theory of concurrent processes with nondeterminism, priorities, probabilities and time. Theoretical Computer Science, 202(1-2):1–54, 1998.
[15] B. Plateau. On the stochastic structure of parallelism and synchronization models for distributed algorithms. In Proceedings of the 1985 ACM SIGMETRICS Conference on Measurements and Modeling of Computer Systems, pages 147–154, Austin, Texas, 1985. ACM Press.
[4] C. Bertolini, L. Brenner, P. Fernandes, A. Sales, and A. F. Zorzo. Structured Stochastic Modeling of Fault-Tolerant Systems. In International Symposium on Modelling, Analysis and Simulation on Computer and Telecommunication Systems (MASCOTS’04), pages 139–146, Volendam, The Netherlands, October 2004. [5] L. Brenner, P. Fernandes, B. Plateau, and I. Sbeity. PEPS2007 - Stochastic Automata Networks Software Tool. In International Conference on the Quantitative Evaluation of Systems (QEST’07), pages 163–164, Edinburgh, UK, 2007. IEEE Press. [6] S. Donatelli. Superposed stochastic automata: a class of stochastic Petri nets with parallel solution and distributed state space. Performance Evaluation (PEVA), 18(1):21–36, 1993. [7] F. L. Dotti, P. Fernandes, A. Sales, and O. M. Santos. Modular Analytical Performance Models for Ad Hoc Wireless Networks. In International Symposium
[16] Y. Saad and M. H. Schultz. GMRES: A Generalized Minimal RESidual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing, 7(3):856–869, 1986. [17] A. Sales and B. Plateau. Reachable state space generation for structured models which use functional transitions. In International Conference on the Quantitative Evaluation of Systems (QEST’09), pages 269–278, Budapest, Hungary, September 2009. [18] K. Sayre. Improved techniques for software testing based on Markov chain usage models. PhD thesis, University of Tennessee, Knoxville, USA, 1999. [19] R. E. Shannon. Introduction to the art and science of simulation. In Conference on Winter Simulation (WSC’98), pages 7–14, Los Alamitos, USA, 1998. [20] W. J. Stewart. Probability, Markov Chains, Queues, and Simulation. Princeton University Press, USA, 2009.
