A comparison of numerical splitting-based methods for Markovian dependability and performability models

A Comparison of Numerical Splitting-based Methods for Markovian Dependability and Performability Models ? V´ıctor Su˜n´e and Juan A. Carrasco Departament d’Enginyeria Electr`onica, Universitat Polit`ecnica de Catalunya, Diagonal 647, plta. 9, 08028 Barcelona, Spain Abstract. Iterative numerical methods are an important ingredient for the solution of continuous time Markov dependability models of faulttolerant systems. In this paper we make a numerical comparison of several splitting-based iterative methods. We consider the computation of steady-state reward rate on rewarded models. This measure requires the solution of a singular linear system. We consider two classes of models. The first class includes failure/repair models. The second class is more general and includes the modeling of periodic preventive test of spare components to reduce the probability of latent failures in inactive components. The periodic preventive test is approximated by an Erlang distribution with enough number of stages. We show that for each class of model there is a splitting-based method which is significantly more efficient than the other methods. 1 Introduction Continuous time Markov chains (CTMCs) are widely used for dependability modeling. For these models, several measures of interest can be computed from the solution vector of a linear system of equations. Typically, such a system is sparse and may have hundreds of thousands of unknowns, so it must, in general, be solved numerically using an iterative method. Several currently available tools allow us to solve dependability models. These are, among others, SAVE [7], SPNP [3], UltraSAN [5] and SURF-2 [2]. SPNP uses Successive Overrelaxation (SOR) with dynamic adjustment of the relaxation parameter ω[4]. SAVE uses SOR for the computation of the steady-state probability vector and SOR combined with an acceleration technique [10] for computation of mean time to failure (MTTF) like measures. UltraSAN offers a direct method with techniques to reduce the degree of fill-in and SOR, being ω selected by the user. Finally, SURF-2 uses the gradient-conjugate method (see, for instance, [14]). ?This work has been supported by the Comisi´on Interministerial de Ciencia y Tecnolog´ıa (CICYT) under the research grant TIC95–0707–C02–02. R. Puigjaner et al. (Eds.): Tools’98, LNCS 1469, pp. 154–164, 1998. c Springer-Verlag Berlin Heidelberg 1998 A Comparison of Numerical Splitting-based Methods for Markovian Models 155 Several papers have compared numerical methods for solving the linear systems of equations which arise when solving CTMC models. In an early paper [8], performance models are considered and several iterative methods are compared for the computation of the stationary probability vector of an ergodic Markov chain. These methods include Gauss-Seidel (GS), SOR, block SOR and Chebyshev acceleration with GS preconditioning. For SOR, an algorithm based on the theory of p-cyclic matrices [17] is used to select a value for ω. In [14], failure/repair models are considered and SOR with dynamic adjustment of ωalso based on the theory of p-cyclic matrices is compared with GS and the power methods, showing that SOR is considerably more efficient specially for the linear systems arising in MTTF computations. In [11] a number of direct and iterative methods are reviewed in the context of performance models. Among others, three spliting-base methods are considered: GS, SOR and symmetric SOR. In [6] the generalized minimal residual method and two variants of the quasi-minimal residual algorithm are compared. In [9], direct and splitting-based iterative methods are considered for solving CTMC models arising in communication systems and the authors suggest to use SOR with suitable values for ωin combination with some aggregation/disaggregation steps. In this paper we compare splitting-based iterative methods for the solution of linear systems which arise in the computation of the steady-state reward rate (SSRR) defined over rewarded CTMC models. We start by defining formally the measure and establishing the linear system which has to be solved. Let X={X(t); t≥0}be a finite irreducible CTMC. Xhas state space Ωand infinitesimal generator Q=(qij)i,j∈Ω. Let ri,i∈Ωbe a reward rate structure defined over X. The steady-state reward rate is defined as: SSRR = lim t→∞ E[rX(t)] and can be computed as SSRR =X i∈Ω riπi, where π=(πi)i∈Ωis the steady-state probability distribution vector of X, which is the only normalized (kπk1= 1) solution of: QTπ=0,(1) where matrix QTis singular and the superscript Tindicates transpose. The steady-state unavailability is a particular case of SSRR obtained by defining a reward rate structure ri=0,i∈U,ri=1,i∈D, where Uis the subset of Ωincluding the up (operational) states and Dis the subset of Ωincluding the down states. In this paper we are concerned with numerical iterative methods to solve the linear system (1). Two classes of models will be considered. The first class include failure/repair models like those which can be specified by the SAVE modeling language [7]. Basically, these models correspond to fault-tolerant systems 156 V. Su˜n´e and J.A. Carrasco made up of components which fail and are repaired with exponential distributions. There is an state in which all components are unfailed having only outgoing failure transitions. The remaining states have at least an outgoing repair transition. Note that in this class of models the detection of the failure of a component is assumed to be instantaneous, i.e. all failed components are immediately scheduled for repair. In the second class of models which we will consider, failures of “spare” (inactive) components will be detected only when they are tested. Test of spare components will be assumed to be performed periodically with deterministic intertests times. To be able to use CTMCs to represent such systems the deterministic intertests time will be approximated by a K-Erlang distribution with Klarge enough to obtain convergence in SSRR as Kis incremented. We will analyze and compare GS, SOR and block Gauss-Seidel (BGS). A more complete comparison including an efficient implementation of GMRES (see, for instance, [13]) and MTTF like measures can be found in [16]. We will show that GS is guaranteed to converge for (1) when the chain Xis generated breadth-first and the first generated state is exchanged with the last generated state. Also, an algorithm to select dynamically ωin SOR will be briefly reviewed. The rest of the paper is organized as follows. Section 2 reviews the iterative methods. Section 3 analyzes convergence issues. Section 4 presents examples and numerical results and Sect. 5 presents the conclusions. 2 Numerical Methods We are interested in solving a linear system of the form Ax =b,(2) where A=QTand b=0. In the following we will let nbe the dimension of A. We next review splitting-based iterative numerical methods which can be used to solve (2). 2.1 Gauss-Seidel, SOR and Block Gauss-Seidel Splitting-based methods are based on the decomposition of matrix Ain the form A=M−N, where Mis nonsingular. The iterative method is then: x(k+1) =M−1Nx(k)+M−1b, where x(k)is the k-th iterate for x. Both GS and SOR are easily derived by considering the decomposition A= D−E−F, where Dis the diagonal of Aand −Eand −Fare, respectively, the strict lower and upper part of A. GS is obtained by taking M=D−Eand N=F. The iterative step of GS can then be described as: x(k+1) =(D−E)−1Fx (k)+(D−E)−1b, A Comparison of Numerical Splitting-based Methods for Markovian Models 157 or in terms of the components of Aas: x(k+1) i=1 ai,i − i−1 X j=1 ai,jx(k+1) j− n X j=i+1 ai,jx(k) j+bi,i=1,...,n . (3) SOR is obtained by taking M=(D−ωE)/ω and N=(1 −ω)D+ωF/ω. The iterative step of SOR can then be described as: x(k+1) =(D−ωE)−1(1 −ω)D+ωFx(k)+(D−ωE)−1ωb, or in terms of the components of Aas: x(k+1) i=ωx GS i+(1−ω)x(k) i,i=1,...,n , where xGS iis the right-hand side of (3). BGS is the straightforward generalization of GS when the coefficient matrix, the right-hand side and the solution vector of (2) are partitioned in pblocks as follows: A=   A1,1... A1,p . . ..... . . Ap,1... Ap,p   ,x=   x1 . . . xp   ,b=   b1 . . . bp   . The iterative step of BGS is: x(k+1) i=A−1 i,i − i−1 X j=1 Ai,jx(k+1) j− p X j=i+1 Ai,jx(k) j+bi,i=1,...,p . (4) Hence, each iteration of BGS requires to solve psystems of linear equations of the form Ai,ixi=zi. Depending on the sizes of the matrices Ai,i, such systems may be solved using direct or iterative methods. 2.2 An Algorithm for the Optimization of ωin SOR In this section we briefly describe an algorithm for the optimization of the relaxation parameter ωof SOR. The algorithm does not assume any special property on the matrix of the linear system and searches the optimum ωin the interval [0,2]. The algorithm is based on estimations of the convergence factor (modulus of the sub-dominant eigenvalue of the iteration matrix). After each iteration k such that the last two iterations have been performed with the same value of ω, the convergence factor ηis estimated as: ˜η=kx(k)−x(k−1)k∞ kx(k−1) −x(k−2)k∞ . 158 V. Su˜n´e and J.A. Carrasco Stabilization of ˜ηis monitored and it is assumed that a good estimate has been achieved when the relative difference in ˜η/(1 −˜η) is smaller than or equal to a given threshold parameter TOLNI three consecutive times. For a given ω, except ω= 1, for which no limit is imposed, a maximum of M= max{MAXITEST, est/RATIOETAST}iterations are allocated for the stabilization of ˜η/(1 −˜η), where est is the number of iterations required for the stabilization of ˜η/(1 −˜η) for ω= 1. If after Miterations ˜η/(1 −˜η) has not been stabilized, SOR is assumed not to converge for the current ω. Selection of appropriate values for TOLNI is a delicate matter. If TOLNI is chosen too large an erroneous estimate of the convergence factor may result and the optimization method may become confused. If TOLNI is chosen too small non-convergence may be assumed when the method converges but ˜η/(1−˜η) takes a large number of iterations to stabilize. Selection of values for MAXITEST and RATIOETAST also involves a tradeoff. If the resulting Mis too small, non-convergence may be assumed erroneously. If the resulting Mis too large, iterations may be wasted for a bad ω. After some experimentation we found TOLNI =0.0001, MAXITEST = 150, and RATIOETAST = 5 to be appropriate choices. The algorithm starts with ω= 1 and, while the estimate for ηdecreases, makes a scanning in the interval [1,2] taking increments for ωof 0.1. If a minimum for ηis bracketed, a golden search (see, for instance, [12]) is initiated. If for an ωit is found that the method does not converge, the increment for ωis divided by 10 and the search continues to the right starting from the last ωfor which the method converged. This process is repeated till the increment for ωis 0.001 (the minimum allowed). If the estimate for ηfor ω>1 is found to increase a similar scanning is made to the left in the interval [0,1]. At any point of the algorithm, the best ωis recorded and used till convergence or the maximum number of allowed iterations is reached when the algorithm becomes “lost”. The complete algorithm is implemented using an automaton with 13 different states corresponding to different states of the search. Due to lack of space we cannot give a precise description of the automaton, but only highlight the main ideas on which it is based. 3 Convergence −QTis a singular M-matrix and it is well known that SOR converges for −QTπ=0(and, therefore, for (1)) if 0 <ω<1 [15, Theorem 3.17]. We prove next that if Xis generated breadth-first and the first state and the last one are exchanged, convergence of GS when solving (1) is also guaranteed. First we briefly describe breadth-first generation. The initial state is put in an empty FIFO queue. From that point, the generation process continues by taking a state from the queue, generating all its successors and putting in the queue the successors not previously generated. The generation process finishes when the queue becomes empty. Given the n×nmatrix A, its associated directed graph Γ(A)=(V,E)is defined by a set of vertices V={1,...,n}and a set of edges E={(i, j)∈ A Comparison of Numerical Splitting-based Methods for Markovian Models 159 V|ai,j 6=0}. A sequence of vertices α=(i0,i 1,...,i l,i 0) is called a cycle of Γ(A)ifij6=ik,j6=k,0≤j, k ≤land, (ik,i (k+1) mod (l+1))∈E,0≤k≤l. The cycle is said to be monotone increasing if i0<i 1<...<i land it is said to be monotone decreasing if i0>i 1>...>i l[1]. Theorem 1. Let Qbe the infinitesimal generator of a finite and irreducible CTMC obtained by generating the CTMC breadth-first and exchanging the first and last states. Forward1Gauss-Seidel converges for the linear system QTπ=0 for each initial guess π(0). Proof. Let Γ(Q)=(V,E) be the directed graph associated to Q, infinitesimal generator of X, and let in=nbe the index of the last state of X. Also, let in−1be any state with (in−1,i n)∈E. Because of irreducibility of Q, some state in−1exists. More generally, since Xis generated breadth-first, each state ilhas a predecessor which was generated before it. Such a precedessor may be in(the state from which Xis generated) or il−1<i l. Then, it is clear that a cycle (in,i k,...,i n−2,i n−1,i n) with il−1<i l,1≤k<l≤ncan be formed in Γ(Q). Such a cycle becomes (in,i n−1,i n−2,...,i k,i n)inΓ(QT), which is monotone decreasing. Then [1, Corollary 1], forward Gauss-Seidel converges for solving QTπ=0for each π(0).ut Of course, with π(0) =0the method will converge to 0, a trivial solution of the linear system which is not of interest. Thus, we should start with any π(0) >0. Regarding BGS when applied to solve (1), it is known that there always exists a convergent block splitting for QTprovided that an appropriate ordering of the states is used [8]. As convergence test we require the relative variation on SSRR to be smaller than or equal to a specified tolerance three consecutive times. The rationale for this test is that it takes into account only “important” components of the solution vector. 4 Numerical Results In this section we compare the performance of the numerical methods to solve (1) using two examples. The first example is a model with failure and repair transitions and immediate detection of component failures; the second example is a model with failure and repair transitions and K-Erlang intertests time of spare components. For all methods, the relative tolerance for convergence is taken =1×10−8 and a maximum of 100,000 iterations is allowed. In all cases, the CTMC is generated breadth-first and, when the GS and SOR methods (SOR reverts to GS when it cannot find an appropriate ω6= 1) are used for the solution of (1), 1The method we have called Gauss-Seidel should be more properly called forward Gauss-Seidel. 160 V. Su˜n´e and J.A. Carrasco the first state (state with all components unfailed) is exchanged with the last generated state so that convergence of GS is guaranteed (Theorem 1). CPU times have been all measured on a Ultra 1 SPARC workstation. The first example is the distributed fault-tolerant database system depicted in Fig. 1. The system includes two processors, two controllers and three disk clusters, each with four disks. When both processors are unfailed, one of them is in the active state and the other in the spare state. Similarly, when both controllers are unfailed, one of them is active and the other spare. The system is operational if at least one processor, one controller and three disks of each cluster are unfailed. Processors, controllers and disks fail with constant rates 2 ×10−5, 2×10−4and 3×10−5, respectively. Spare components fail with rate 0.2λ, where λis the failure rate of the active component. There are two failure modes for processors: “soft” mode, which occurs with probability 0.8, and “hard” mode, which occurs with probability 0.2. Soft failures are recovered by an operator restart, while hard failures require hardware repair. Coverage is assumed perfect for all failures except those of the controllers, for which the coverage probability is C. Uncovered controller failures are propagated to two failure free disks of a randomly chosen cluster. Processor restarts are performed by an unlimited number of repairmen. There is only one repairman who gives preemptive priority first to disks, next to controllers, and last to processors in hard failure mode. Failed components with the same priority are taken at random for repair. Repair rates for processors in soft and hard failure mode are, respectively, 0.5 and 0.2. Controllers and disks are repaired with rates 0.5 and 1, respectively. Components continue to fail when the system is down. The measure of interest is the steadystate unavailability (UA), a particular case of the SSRR generic measure. The generated CTMC has 2,250 states and 19,290 transitions. Four values for the coverage probability are considered: C=0.9, 0.99, 0.999, and 0.9999. For this example we only experimented with GS and SOR. The CPU time required for the generation of the model was 0.263 s. processors controllers disks D1 D2 D3 CO PP Fig. 1. Distributed fault-tolerant database system A Comparison of Numerical Splitting-based Methods for Markovian Models 161 In Table 1 we show the number of iterations and CPU times in seconds for the first example. The GS method is the faster. SOR requires the same number of iterations to achieve convergence as GS because the convergence is so fast that it is achieved before any adjustment on ωcan be done. The time per iteration for SOR is slightly greater than for GS. Table 1. Number of iterations under, respectively, GS and SOR, itGS,itSOR, CPU time in seconds under, respectively, GS and SOR, tGS,tSOR, and UA for the first example CitGS tGS itSOR tSOR UA 0.9 20 0.097 20 0.12 4.054 ×10−5 0.99 19 0.094 19 0.12 4.461 ×10−6 0.999 19 0.094 19 0.12 8.537 ×10−7 0.9999 19 0.094 19 0.12 4.929 ×10−7 The system considered in the second example is exactly the same as the system of the first example, with the only difference that spare processors and controllers are tested with deterministic intertests times Tapproximated by a K-Erlang distribution with expected value Twith Klarge enough to make the approximation error small. The measure of interest is again the steady-state unavailability UA. It is clear that the greater Tthe greater UA. Intuitively, the fact that faulty spare units are not immediately scheduled for repair “increases” the repair times of such units and so increases UA. Then, only values for the intertests time not much greater than the average repair times of the components are reasonable choices. Since the minimum repair rate is 0.2, we consider the following five values for T: 100, 10, 1, 0.1, and 0.01. For the sake of brevity, we will give only results for a coverage probability Cequal to 0.99. The value of K is chosen as the minimum value which makes the relative difference between UA for two consecutive K’s smaller than or equal to 5 ×10−4.InTable2weshow, for each value of T, the number of Erlang stages K, the number of states and transitions of the CTMC Xand its generation time. Table 2. Number of Erlang stages K, number of states, number of transitions and generation time tgin seconds for the second example and C=0.99 TKstates transitions tg 0.01 3 12,000 125,766 1.98 0.1 3 12,000 125,766 1.98 1 6 24,000 251,532 4.14 10 9 36,000 377,298 6.36 100 25 100,000 1,048,050 18.7 162 V. Su˜n´e and J.A. Carrasco The state descriptions of the second example have a component φ,1≤φ≤K used to indicate the phase of the K-Erlang distribution. For BGS, the blocks are chosen to include all states which only differ in the value of the state variable φ. In addition, states within each block are sorted following increasing values of φ (from 1 to K). With that ordering, the diagonal matrices Aii of BGS have the form:         qm,m 0... 0qn,m qm,m+1 qm+1,m+1 0... 0 qm+1,m+2 qm+2,m+2 ... 0 ... 0... q n−2,n−1qn−1,n−10 0... 0qn−1,n qn,n         . Taking advantage of this form, we solve efficiently the linear systems (4) of BGS using Gaussian elimination with fill-in only in the last column. The iterative methods considered for the second example are GS, SOR and BGS. Table 3 shows the results obtained. Notice first that although UA tends fast to the value corresponding to instantaneous detection of failed spare components (4.461×10−6), its dependence on Tis significant, at least for moderate values of the intertests time. The performance of the numerical methods is also affected by T. For large values of T, the GS method performs very well, but its performance degrades quickly as Tdecreases. The same type of comments can be made for the SOR algorithm. Note, however, that as the number of iterations required by GS increases, the relative reduction in the number of iterations achieved by SOR is greater. This means that the algorithm used for selecting the relaxation parameter ωis efficient. BGS is the method which requires less iterations. For T= 100 it requires more CPU time than GS and SOR. This is due to the time required to sort the states as explained before. For T= 10 BGS is as fast as GS and SOR, and for smaller values of Tit should be clearly considered as the method of choice. Overall, BGS seems to be the method of choice for the second example. Table 3. Number of iterations under, respectively, GS, SOR and BGS, itGS,itSOR, itBGS, CPU time in seconds under, respectively, GS, SOR and BGS, tGS,tSOR,tBGS, and UA for the second example, C=0.99 and several values of T TitGS tGS itSOR tSOR itBGS tBGS UA 100 21 10.5 21 11.9 10 15.6 6.129 ×10−6 10 31 5.01 31 5.74 11 4.92 4.641 ×10−6 1 162 14.8 162 17.5 12 3.15 4.480 ×10−6 0.1 1,391 58.4 1,045 51.4 12 1.43 4.463 ×10−6 0.01 12,632 527 5,953 284 12 1.45 4.461 ×10−6