572 zyxwvutsrqponmlkji zyxw zyxw zyxwv IEEE TRANSACTIONS ON RELIABILITY, VOL. 40, NO. 5 , 1991 DECEMBER Estimation of Network Reliability Using Graph Evolution Models T. Elperin Ben-Gurion University of the Negev, Beer Sheva I. Gertsbakh Ben-Gurion University of the Negev, Beer Sheva zyx zyxwvu zyxwvut zyxw The usual network reliability problem is: Suppose that edges of G are s-independently erased with probabilities q ( e ) , e E E. Let F denote the set of non-erased edges; the subgraph ( V J ) of G then appears with probability: M. Lomonosov < ( e ) neE Ben-Gurion University of the Negev, Beer Sheva Pr{F) = Key words - Terminal reliability, Spanning tree, Lifetime In practice this static model describes: i) systems without edge renewal, or fi) stationary regimes Of renewable systems. h Case ii, q ( e ) is the equilibrium probability that a renewable edge e is down. simulation, Kruskal algorithm, Graph evolution, Monte Carlo method Reader Aids Purpose: Widen and advance state of the art. Special math needed for explanations: Basic graph theory, probability and reliability theory. Results useful to: System reliability analysts. neE E-F q(e) (1-1) Notation R ( 3 t , q ) reliability of Tc the m-vector with coordinates q ( e ) E E. 4 T a fixed subset of V; the members of T are called terminals Abstract - Monte Carlo techniques for estimating various net- Other, standard notation is given in “Information for Readers work reliability characteristics, including terminal C O M ~ C ~are ~ V ~& ~ YAuthors” , at the rear of each issue. developed by assuming that edges are subject to failures with arThe reliability of 3t is defined as the probability R ( Z , q ) bitrary probabilities and nodes are absolutely reliable. The core of our approach is introducing network time-evolutionprocesses that the random subgraph ( V , F ) is operating. As a principal and using certain graph-theoretic machinery resulting in a con- example, we consider terminal reliability defined by: ( V J ) is siderable increase in accuracy for Monte Carlo estimates, especially up if Tlies in one component of ( V J ) . This network is denoted for highly reliable networks. Simulationstrategies and numerical by (G9T). results are presented and discussed. We consider Monte Carlo simulation of network reliabili- 1. INTRODUCTION Reliability of networks with randomly failing edges is a subject of extensive research. Several directions prevail in this research: algorithms for reliability computation [1-3,2 1,221 reliabdity estimation by means of simulation [6,8-11,14,16,17] constructing tractable lower and upper bounds on the network reliability [2-4,18-201. This paper describes an approach to network reliability simulation based on an artificial time-evolution formulation of network failure reliability characteristics. This approach incorporates both simulation and analytic methods and has no difficulty when the edge failure probabilities are distinct. The aim of this work is to develop efficient simulation procedures. In particular, we are interested in reducing the relative error of network failure estimation for highly reliable networks. A network 3t is an undirected graph G = ( V , E ) ,with nodeset V, J V J= n , and edge-set E, J E J= m , whose spanning subgraphs ( V J ) are classified as up (operational) and down (non-operational), subject to the reasonable monotonicity condition: If ( V,F) is up then all subgraphs ( V J ‘ ) are up, where F E Ff. ty based on time evolution modification of the static model. This modification leads to an important acceleration (variance reduction) of the Monte Carlo procedures, and guarantees boundedness of the relative error, irrespective of the value of graph unreliability (claim 6.2). Section 2 presents a general statistical framework for reliability evaluation. Section 3 expresses the network reliability R ( T c , q ) in terms of two different Markov processes: i) a destruction process (DP), and ii) a creation process (CP). Section 4 presents combinatorics for analyzing DP for the network ( G, T ). The central role belongs to the notion of maximal spanning tree, the Kruskal algorithm, and lemma 4.1. In the case of equal edge-failure probabilities, the DP leads to an efficient Monte Carlo sampling scheme studied earlier by Fishman [lo]. Sections 5 & 6 describe a modification of CP based on the graphtheoretical notion of closure. Section 7 presents numerical results for a family of networks, with a comparison of several Monte Carlo approaches. zyxwvutsrqpon 2. MONTE CARLO SAMPLING SCHEME By the Monte Carlo method for evaluating a sum, z z(u) = U€ U 00 18-9529/91/ 1200-0572$01.WO1991 IEEE zyxwvutsrqponmlkjih ~ zyxwvu I zyxwvut zyxwvuts zyxwvutsrq ELPERIN ET AL.: ESTIMATION OF NETWORK RELIABILITY USING GRAPH EVOLUTION MODELS over a very large set Uof “outcomes”, we mean the following method. Introduce the probability distribution p ( U ) on U and consider U as an um from which a ball u can be drawn with probability p (U). Also let Z denote the mean value of random variable Y(u) = z(u)/p(u) by: The variance and coefficient of variation of Y are: 513 balls u in this scheme are the trajectories of a certain Markov process on the state space 2“ or its proper reduction, and the value Y( U ) of the random variable is the conditional probability of the up state for a given trajectory. Except for a special choice of the urn we suggest no changes to the above basic sampling scheme. This, however, provides performance of the Monte Carlo which is in general comparable with the existing advanced sampling techniques. In certain cases, such as highly reliable networks and dense graphs, the suggested method is definitely better. 3. GRAPH DESTRUCTION AND CREATION PROCESSES zyx zyxwvutsrq zyx zyxwv zyxwvuts From basic statistics we have Introduce an artificial time t and let F ( t ) denote the set of edges existing at the instant t. Consider two types of graph evolution processes G ( t ) = ( V , F ( t ) ) , F ( t ) E E, t 1 0. - Ckim 2.1. Let S = (U’, u2,...,UN) be the result of N s-independent choices from U,with probabilities p (U). Then N Ps = N-’ Y(Ui) i=l is an unbiased estimate of E { Y) , with variance and coefficient of variation equal: Destruction Process (DP) Initially, at t =0, all edges are up: F ( 0) =E. Edges leave the set F ( t ) s-independently, at random moments T( e ) , with the Cdf, Pr{.r(e) It} = 1 -exp( -X(e)t). Let 7(X) denote the random moment when G( t ) goes down. The Cdf of T( X ) is Pr{7(X) It} = R ( X , q ) , where q is an m-vector with the components q(e) = 1 -exp( -X(e)t), e E E. The static model in the Introduction agrees with DP when the edge failure rates X(e) are chosen so that - Var{E} = N-’ * Var{Y> 6s = 0 and is realized at t = 1. Crude Monte Carlo Creation Process (CP) Consider a network X with a graph G = ( V,E ) and some operational (up) criterion. Realization of the above scheme for the set of subsets of E as the urn U,with p (F), F E E, given by (1-1), and Y(F) = 1 when (VJ) is up, and 0 otherwise, are referred to as crude Monte Carlo (CMC) for evaluating R ( X , q ) , or equivalently, for evaluating Q(X , q ) = l-R(X,q). The variance of CMC is: Initially, at t =0, all edges are down: F( 0) = 0. The edges of G join F( t) s-independently, at random moments 7 ( e ) , with the Cdf, Pr{.r(e) It} = 1-exp( -X(e)t), e E E, and operate forever. Let [ ( X ) denote the moment when G ( t) goes up. TheCdfof[(X) i s P r { t ( X ] It} = R ( X , q ) ,withq(e) = exp( -X(e)t). The static model agrees with CP when the edge birth rates h ( e ) are chosen so that, and the relative error in evaluating Q on the basis of N s-independent experiments is: The main deficiency of CMC is the unbounded growth of BCMC as Q approaches 0 (viz, for highly reliable networks). Various improvements of CMC have been suggested in order to reduce or eliminate this effect [6,8-11,14,17]. In this paper we offer another urn scheme for evaluating network reliability which guarantees finite relative error. The and is realised at t = 1. For each of these processes consider an ordering (permutation) w = (el,...,e,,,) of E specifying the order in which the edges are erased (in DP) or created (in CP). The probability of w is given by the well-known expression [7,17]: (3-3) where Eo = E, Ei = E - e l - ...- ei, 1 Ii s m-1, and h(Ei) = L E E , h(e)* For a given w ,an edge e is called DP-critical if erasing it causes G( t) to go down, and CP-critical if its creation causes 574 zyx zyxwvutsrqponmlkj IEEE TRANSACTIONS ON RELIABILITY, VOL. 40, NO. 5, 1 9 9 1 DECEMBER zyxwvutsrqp zyxwvutsrq zyxwvutsr zyxwvutsr zyxwvut zyxw zyxwvu zyxwvu zyxwvutsr G(t) to go up. The ordinal number of the critical edge in w is called the critical number of w and denoted by [w], so that min{i:G-el- ... -ei is down}, in DP (el,e2,...,ei)) is up), in CP. PutP(t(w) = Pr{7(32) It I (3-4) tJw} for DP, and Pr(t(32) I w) for CP. By a well-known property of Markov processes The following obvious variance decomposition reveals the gain in accuracy provided by DP or CP with respect to CMC (at the expense of more complex computations). VarcMc = Varp + Pr{x} P(tlx) P(tlx). (3-9) X Similarly, a tail y can be identified with the bundle of permutations w satisfying tl( w) =y. Its probability is: [7], P ( t ( w ) is a convolution of exponential r.v.’s: i=r h(ei)+ ... +h(e,) P ( t ( w ) does not depend on the order of e[wl+l,...,em in the permutation w. The following notions are therefore reasonable. Trajectories and Tails An ordered subset x = (el, e 2 , ...,er) of E is called a rrajectory of DP if G-el-e2- ... -ei is up for i < rand down for i = r; x is called a trajectory ofCP if ( V , {el,...,ei)) is down for i < r and up for i = r . An ordered subset y = (er, e,, l,. ..,em)of E is called a tail of DP if ( V, {ei,ei+1,. . .,e,} is down for i > r and up for i = r; y is called a tail of CP if G-ei-ei+l-... -em is up for i > r a n d down for i = r . The Cdf of the critical moment, given y , is: (3-10‘) 0 The critical number r=[w] divides a permutation w = (el,...,em) into the trajectory tr ( w ) = ( e1,e2,...,e,.) and the tail tl(w) = (er,er+l,...,em).A trajectory x can be identified with the set (bundle) of permutations w satisfying tr (w) = x, thus - In the important case of equal edge-failure probabilities (X(e) = X for all e E E) we have - P(tlx) = P(tlw) for w satisfying tr(w) = x , n r Pr{4 = Pr{x} = w:tr(w) = x i=l X(ei) . h ( E - e l -... -ei-l) (3-6) By the total probability formula, 1 (m-r)! ( r - l)! Pr{w} = -, Pr{x} = ____ , Pr{y} = -. m! m! m! (3-11) Now we describe generating permutations, trajectories and tails with their “natural” probabilities, as they appear in the corresponding process. Generating permutations q = {q(e ) ; e E E } is the vector of edge failure probabilities, de) = i 1-exp[-X(e)t], exp[-h(e)t], for DP for CP. For each edge e generate a value b(e) of r.v. 7(e), the lifetime of e. Then the desired permutation w = ( el,e2,...,em) is induced by the inequalities: b(el) The sum at the left in (3-7) is over all trajectories of the corresponding process. The Monte Carlo scheme based on generating trajectories x and exactly computing P ( t Ix) is characterised by the variance: (3-8) < b(e2) <...< b(e,). (3-12) This method is equivalent to drawing a permutation w from the um (see section 2) of all m! possible permutations of E with 0 probability p ( w ) given by (3-3). Generating trajectories A trajectory x = (e1,e2,...,er)is generated by sequentially choosing el from E with probability X ( e l ) /X( E ) , e2 from E-el with probability h(e2)/h(E-el), etc, until the critical zyxwvutsrqponmlkji zyxwvu zyxwvutsrq zyxwvuts zyxwvu zyxwvutsr I ELPERIN ET AL.: ESTIMATION OF NETWORK RELIABILITY USING GRAPH EVOLUTION MODELS edge is generated. This is equivalent to drawing a trajectory from an urn with probability (3-6). Sequentially generating tails in the reverse order e,, e, - l,. .., with natural probabilities (3- 10) is practically intractable, except for the case of equal edge-failure probabilities. In the latter case, tails are generated in the same sequential manner as trajectories. 0 4. IDENTIFYING THE CRITICAL EDGE OF DP FOR 3Z=(G,T) 515 the subset C of E connecting D ' & D " becomes empty. If it did not, then there would exist ej in C, with j >k, so that the spanning tree D, - t?k ej is lexicographically greater than D,-which is a contradiction. + Example Figure l a presents the graph called dodecahedron [12]. It has 20 nodes and 30 edges. The double-circled nodes 1,3,17 are the terminals. The numbers near edges specify the edge lifetimes, b ( e ) .The corresponding tree D, is shown on figure lb; its bold part is D,( T ) . The network lifetime is 19, while the overall connectivity is lost earlier, at t=6. zyxwvutsrqpo zyxwvutsr Notation component of G - el - e2- ... - ek- that contains T some spanning tree of Gk the minimal subtree of D that contains T the lexicographically maximal spanning tree with respect to w e (w,T ) the junior edge of the subtree D, ( T ) . Gk D D(T) D, Consider generating a trajectory x = ( e1,e2,. .. ) of DP. After a current edge ek is erased, we need only to check if the terminal-set T is connected by G -el - e2- ... - ek. Surely, the choice of ek can always be restricted to the edge-set of Gk. The possibilities are: i. ek 4 D. Then Gk+l=G-ek, and obviously Gk+l contains T; D might be preserved. ii. otherwise. Find an edge e of Gk-ek connecting the components of D-ek. If e exists, then Gk+l =G-ek, Gk+l connects T, and we can put D: =D - ek e; otherwise (e does not exist), consider the components of D - ek. If one of them, say D ' contains the entire T, then Gk+ = Gk( D ' ) which is the subgraph of Gk induced by the vertices of D ' , Gk+ connects T, and we can put D: =D'; otherwise T is disconnected in Gk-ek, so that the critical number is k, and (el,e2,...,ek) is the desired trajectory. 0 + a. An important fact is that when an entire permutation w is available, the critical edge can be determined by using exactly one special spanning tree D, as shown in lemma 4.1 below. For a spanning tree D,D(T ) is the union of the chains of D between all pairs of terminals. An edge permutation w induces the following lexicographic order among the spanning trees: D >D " when the senior edge of D ' -D " is greater than that of D " -D Then D, is exactly the tree constructed by the famous Kruskal algorithm [15], with the input w. I . Lemma 4.1. Let b (e) be the lifetime of edge e in DP, and w = ( el,e2,...,e,) be the edge permutation induced by inequalities b ( e l ) < b ( e 2 )e . . .< b(6,). Then ( G , T ) fails when e(w,T) fails, so that e ( w , T ) is the critical edge of w. Proof. Let k be the ordinal number of e ( w,T)in w ,and putFi=E-el- ...- ei, i = l , ...,k. For i < k , (V,Fi)contains D,( T ) and thus connects T. We show that T is disconnected in ( v,Fk). Consider the components D ' , D " of D, - ek By definition of ek, both D ' , D" contain terminals. When ek fails, b. Figure 1. The Dodecahedronand Its Maximal Spanning Tree Simulation Strategy The following Monte Carlo simulation algorithm is suggested by lemma 4.1. 576 zyxwvutsrqponml zyxwvutsrqpo zyxwvutsrq zyxwvutsrq IEEE TRANSACTIONS ON RELIABILITY, VOL. 40, NO. 5, 1991 DECEMBER 1. Generate N permutations wj (j= 1,..., N ) as described in (3-12). 2. For each j = 1,. ..,N, determine [wj] using the Kruskal algorithm; the trajectory tr ( w j ) is thus identified. 3. Compute P(tlwj) by (3-5). 4. Compute the estimate of network failure probability as the sample average of P ( t 1 w j ) (j= 1,. .., N ) , see (2-2); compute the corresponding sample variance. 0 The Kruskal algorithm used in the above simulation starts from e, and terminates with e,, ie, it deals only with the tail of the given permutation. In principle, this property could be exploited by sequentially generating the tail of a permutation wj, applying (3-10) & (3-10'). Unfortunately, however, these formulas are too complicated for straightforward calculations. The situation is, however, much easier in the important case of equal edge-failure probabilities, as in (3-1 1). The reliability simulation of 32 = ( G,Z') with equal edge-failure probabilities is considered by Fishman [8,10]. In [lo], the simulation method is based on generating tails, and lemma 4.1 was used, but without being formulated explicitly. The case of equal edge-failure probabilities is extremely favorable for the network Monte Carlo. Indeed, in this case all permutations have the same probability, 1/ (m!). On the other hand, assuming, without loss of generality, h ( e )= 1 for all e E E, the convolution (3-5) can be considerably simplified. Indeed, in the theory of order statistics [5] it is well-known that 1-exp( - ( m - i + l ) t ) istheCdf0fther.v. ~ = T ( ~ ) - T ( ~ . where T ( ~ is) order statistic i for the sample of m i.i.d. r.v.'s T~ exp(l), q,=O. Moreover, ther.v.'s V,, i = l , ...,m, are s-independent. Thus the r.h.s. of (3-5) is: TABLE 1 Simulated ID'Sof Complete Graphs, N = lo6 zyxwvu A ( r ) s for ~ 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 ~ ) , 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 0.023810 0.095238 0.002050 0.285714 0.010155 0.000127 0.595238 0.029820 0.000761 0.000006 0.071937 0.002688 0.000043 0.155212 0.007294 0.000165 0.298834 0.016465 0.ooO529 0.431992 0.033110 0.001451 0.061712 0.003070 0.108629 0.006247 0.181284 0.011650 0.277702 0.020218 0.310228 0.034246 0.055661 0.086821 0.130692 0.186388 0.241062 0.221751 0.000002 0.000002 O.ooOo16 0.oooD29 O.ooOo87 0.000259 0.000549 0.001019 0.001869 0.003360 0.005759 0.009628 0.014864 0.023444 0.035028 0.051427 0.073990 0.103387 0.138312 0.177942 0.201261 0.157776 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 zyxwvuts zyxwvu zyxwvutsrq zyxwvutsrq - The Cdf Hr(t) of the r-th order statistic is given by a wellknown formula Hr(t) = Pr{T(,)st) = binfc(r;l-e-',m) (4-2) Thus (3-7) acquires the form: m R(=,q) = A(r)Hr(t), (4-3) 0.000006 0.000005 O.ooOo13 0.000041 0.000065 0.000117 0.000216 0.000342 0.000499 0.000757 0.001205 0.001785 0.002634 0.003678 0.005431 0.007685 0.010677 0.014367 0.019461 0.026335 0.035125 0.045961 0.059512 0.075767 0.093734 0.112722 0.130512 0.140335 0.130172 0.080841 r= 1 A ( r ) = (number of permutations with [w] = r)/m!, This presentation is used in [lo] as a basis for a Monte Carlo sampling scheme. A remarkable property of (4-3) is that the relevant combinatorics of 32 = ( G ,Z') , expressed by the numbers { A ( r ), r = 1,...,m}, are totally separated from the probabilities contained in the functions Hr( t) , These functions are standard and always available for any value of t. So, the Monte Carlo simulation efforts should be turned to obtaining the distribution { A ( r ), r = 1,...,m}.It is reasonable to c d it the internal distribution (ID) of the network ( G , n. As an illustration, tables 1 & 2 present the IDS of several complete graphs with T= V, and for the dodecahedredon with various Ts. obtained by simulation. 5 . CREATION AND MERGING PROCESSES FOR 32= (G,Z') A closer look at the performance of the Kruskal algorithm reveals that on each step of constructing a maximal spanning tree, there can be identified a set of irrelevant edges whose future appearance does not affect the time t (32). These are exactly the edges complementing the existing part of the tree to its graphtheoketical closure. The closure of a subset F of E consists of F and all edges of G whose ends lie in the same component of the spanning subgraph ( V , F ) . A subset F is closed if it coincides with its closure. zyxwvu zyxwvutsrqpon ELPERIN ET AL.: ESTIMATION OF NETWORK RELIABILITY USING GRAPH EVOLUTION MODELS TABLE 2 Simulated ID’S for the Dodecahedron, N = IO5 AV) zyxwvutsrqponmlkjih zyxwvu I r T = V T = (1, 7, 8, 11, 16) T = (1, 20) 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 0.00476 0.01637 0.03366 0.05991 0.09559 0.13872 0.17746 0.20031 0.17794 0.09528 0.00130 0 .OO429 0.00959 0.01928 0.03486 0.05867 0.09259 0.13270 0.17232 0.17733 0.13625 0.08407 0.00039 0.00170 0.00374 0.00752 0.01464 0.02600 0.04515 0.07215 0.10635 0.13207 0.13642 0.12568 0.10235 0.07963 0.05561 0.03731 0.02361 0.01428 0.00836 o.oO401 0.00183 0.00089 0.00029 0.m2 577 in one of direct successors of g, say g’, obtained by merging exactly two super-nodes of g, and chosen with probability (A (g) - X(g ’ ) ) /X(g) . The above is summarized in the following. Claim 5.1. i. g(t) = ( F ( t ) ) is a Markov process on L ( G ) ; ii. the time spent by g(t) in a state g is distributed as exp ( (g) 1; iii. thetransitiong g’hastheprobabilityPr{g’Ig} = (X(g) - X(g’))/X(g), when g’ is a direct successor of g, and 0 otherwise. - For additional details see [3,19]. In the following, g ( t ) is referred to as Merging Process, MP. zyxwvutsrqponm zyxwvutsrq 0.04464 0.02017 0.00784 0.00284 o.oO091 0.00028 O.ooOo5 0 . m 2 Example Figure 2 presents L (K4), the set of all regular partitions g of the complete 4-node graph, “naturally” stratified into 4 levels according to the number of super-nodes in g. The arrows show the direct successions in L ( K4), thus forming the transition graph of the Markov process g(t). The members of L ( K 4 ) are represented by circles; the corresponding closed spanning subgraph is drawn in each circle. Let T = {2,3}. The double circles correspond to the partitions for which T lies in one super-node. zyxwvutsrqp zyxwvut For example, closing the set of bold edges of the dodecahedron shown in figure l b adds the edge e = ( 1,3). The closure operation enables us to deal simultaneously with thicker bundles of permutations than in the original DP or CP. A serious obstacle for using this approach in the DP for the network ( G,T ) is the stochastic properties of tails, expressed by (3-lo), (3-10’). This obstacle never appears for (G,T) in CP, as is shown below [3,19]. In what follows, the notion of regularpartition of Vplays the central role. Given a graph G = ( V,E), a partition g = {Xl,X2 ,...,X,} of V, where Xi f l J = 0 for i # j , and U [= 1 Xi = V, is called regular (with respect to G) if each induced subgraph G ( X i ) is connected. Arbitrary set F of edges generates a regular partition (F) = {Xl,X2,. ..,X,} where Xi are the components of the spanning subgraph ( V , F ) (including isolated nodes, if any). Subsets F’ and F“ are equivalent if (F’) = ( F ” ) , and identify every regular partition g with the class of subsets F of E satisfying (F) =g. Clearly, each such class is the collection of subsets of edges with a common closure. For every regular g, let its components be referred to as super-nodes and E(g) denote the set of exteml edges (the edges between distinct super-nodes). Put A (8) = (e). Consider the set L( G ) of all regular partitions of V, partially ordered by the relation: g ’ < g ” when g ” is obtained by merging components of g’. Suppose that a state F ( t ) of CP (see section 3) belongs to an equivalence class g ( t ) = g. Clearly, the time F( t) spends inghasthecdf l-exp[-X(g).t]. Onleavingg,F(t) jumps Figure 2. Transition Diagram for a Markov Process Whose States Are All Regular Partitions of a Complete 4-Node Graph Given a network ( G,T ), we say that g E L( G ) is up if all terminals lie in one super-node of g. A trajectory of g ( t ) is a sequence U = ( go,gl,...,g,) of regular partitions where go is the trivial partition into singletons; gi is a direct successor of gi-l for i = l , ...,r; and r is the first i such that gi is up. In general, trajectories have distinct lengths, so that r depends on U. The probability of U is: 578 zyxwvuts zyxwv zyxwvutsrqponmlk zyx IEEE TRANSACTIONS ON RELIABILITY, VOL. 40, NO. 5 , 1991 DECEMBER r-1 Pr{u} = zyxwvu zyxwvu zyxwvu zyxwvut zyxwvutsrq (si) - (gi+ 1 ) (5-1) E, of the edges between distinct components Xi.If g, is up, then stop; otherwise draw an edge from E, and form g,+1 by merging the two components connected by this edge. After a trajectory u = (go,gl,...,gr} is formed, the condiThe conditional distribution function of 4 (37.1 along u is tional Cdf ( l is computed as the convolution of functions (l-exp[-Ait]) where Ai= CeEErX(e), i=O, ...,r-1. (5-2) P (t 1 U ) = Convo, is 1 { 1 - exp[ - X (gi)t]} . the Since all Ai are distinct (in fact, & > A1 > above convolution is a linear combination of the exponents Finally, the Cdf of [ (31) is exp[ - A$], i =0,..., r - 1, whose coefficients are homogenous (5-3) functions of Ai, of order 0. The following recurrent procedure has the complexity 0 ( r2): i=O h(gi) where U is the set of all trajectories of g(t). Returning to the initial creation process F ( t) E E, we see that a trajectory x = (el,e2, ...) of CP produces a uniquely defined trajectory of MP which we denote by (x). We say that trajectories n’ , x ” of CP are equivalent if ( x ’ ) = ( x ” ) . Thus, a trajectory u of MP represents the class of trajectories of CP satisfying (x) = U ; we write it as x E U. For x E u one has - Assume Conv,-kcicr-l (1 -exp[-Aitl} k = 1- Ak,i exp[-A,-itl i=l k With Ak,j = 1. i=l (To start with, we haveAl,l = 1). Pr {x 1 U } = Pr { x } /Pr (U} Then - The Monte Carlo scheme based on generating trajectories of MP and exactly computing P ( t ( u ) using (5-2) has k Ak+l,k+l = l- Ak+l,i. i= 1 zyxwvu 6. COMPLEXITY OF THE MP-MONTE CARLO Comparison with CMC and CP is based on the expansions: VarcMc = VarMp + Pr(u}P(tJu) P ( t l u ) , (5-5) A. Estimating Q = E( X , q ) for particular value of the vector q with a given mean relative error 6. B. Estimating Q = R(%,q) for a l-parameter family qr oftheformq,(e) = ql(e)‘ = exp[-h(e)t], t E (tl,t2, ...,tk) (5-6) With a given mean relative error 6. The following statements can be easily established. U Varcp = VarMp + Pr{u} U ( Pr(x)u} P(tlx)2-P(tlx)2 . xEu ) For evaluating the computationalcomplexity of the Monte Carlo scheme based on MP consider two problems. The second term in the r.h.s. Of (5-6) is the Part of Varcp Claim 6.1. The complexity of one simulation run for A is by the state space reductionwhen cp was transform- 0(n2). For B, the values of Q ( tl 1 U),...,Q ( tk1 U ) are available ed into MP. in one simulation run: its complexity is 0 ( n 2 ) O(k-n). In the example of figure 2 the trajectories of g ( t ) are the Let 6Mp denote the coefficient of variation in the MP paths starting in 86 and terminating in doubled circles. Monte Carlo scheme: Simulation strategy For estimating the sum (5-3), the sampling scheme of sec- 6$p(Q) = V ~ ~ M - P & U Pr{u}Q(tlu)2 -1, Q2 Q(02 tion 2 is applied, with “natural” probabilities p ( u ) = Pr{u} given by (5-1). Then the complexity of the MP Monte Carlo for both A and B is: Generating a trajectory of g ( t ) + (e) Start from go = {XI,...,Xn},(Xil = l . At step r one has a sequence g0,gl,g2,...,grr with g,= {XI,...,Xn-,},and a list zyxwvu sz,,(Q) . 0(~2). 62 U zyxwvu zyx zyxwvut zy zyxwvutsrqpo zyxwvutsrqp zyxwvutsr zyxwvutsrqp ELPERIN ET AL.: ESTIMATION OF NETWORK RELIABILITY USING GRAPH EVOLUTION MODELS A pleasant feature of MP is given by the following statement. Claim 6.2. For a given n and a given operational criterion the coefficient of variation 8hp(Q) is bounded uniformly for all values of A(e), e E E and 0 5 t IW . Proof. Since the A’s appear only in the products A(e) .t and in 0-homogeneous form, it is sufficient to prove that limt-o, 8Lp(Q) exists and is bounded in the unit ball E, A2 (e) I 1. Consider the lattice L of all partitions of the node-set I/ (states), and let UP and DN denote the sets of up and down states respectively. A trajectory of MP is a sequence U = (go,gl,... ,gr), where go,. ..,gr- E DN, g, E UP (in this formulation trajectories with zero probabilities are permitted). Then - W = (Variance) x (CPU-time in sec for lo00 replications). (7-1) The values WDp & WMp for DP & MP, respectively, were compared with the corresponding WcMc for the crude Monte Carlo (CMC) and for some methods in [6,9,14,16]. The CMC was based on erasing edge e with probability q( e), e E E , and on checking the terminal connectivity of the resulting subnetwork. When no edges fail, the CMC simulation program skips the terminal connectivity check [6]. The set union algorithm of Hopcroft & Ullman [13] was applied for the connectivity check. For both DP & CP, the ratios WCMCIWDp & WCMC/WMp were computed. From the accuracy point of view, the principal parameter is the relative error: 6 = (Variance) “/(Network failure probability). - as t oo (since A(g,) >...>A(gr-l)). Define p ( % ) = min{h(g) : g E DN with a direct successor in UP} and U, = {U = (go,...,gr-l,gr) E U :A(g,-l)=p(%)}. Then- [for t large enough] U€ U, (7-2) The variance reduction factor with respect to CMC was computed for MP as, qMp = 6&C/6hp. The simulation results are based on N = lo4 replications. The followingnetworks were chosen for numerical experiments: 1. The Easton-Wong [6] communicationnetwork with 105 nodes and 127 edges, with the all-terminal connectivity as the operational criterion. 2. The dodecahedron (figure la) with the s-t connectivity as the operational criterion for s = 1 and t =20 [9]. 3. The dodecahedron with the all-terminal connectivity as the operational criterion. 4. A family of complete graphs Klo, K15, K20, K25,K30 with the all-terminal connectivity as the operational criterion. For networks #1 & #2, the performance of our methods is compared with that of alternative methods [6,9]. Table 3 presents simulation results for the Easton-Wong network, table 4 for the dodecahedrons, and table 5 for complete graphs. The following conclusions can be drawn from analyzing tables 3 - 5. zyxwvutsr Pr{u} ~ ( u ) . A(%) = 519 This limit is a continuous function of A( e), e E E , and the assertion follows. 0 In the particular case of complete graph G=K,,, with equal edge-failure rates, and all-terminal connectivity as an operational criterion, it can be shown that A&( Q) I1 for all n . The following seemingly non-trivial question is then reasonable and important. Question. Is there a universal constant A such that 6MP(Q) A for all n, all possible A’s and 0 It e oo? I 7. SIMULATION RESULTS In order to evaluate the performance of the Monte Carlo schemesbased on edge destruction & creation processes, a series of experiments has been done for several networks. Network failure probability Q=R{%,q} was estimated along the lines of the simulation strategies for DP & MP described in sections 4 & 5. As a performance measure of the simulation method we take: zyxwv 1. DP is not competitivewith MP, in terms of both variance reduction factor and in the Wperformance measure. Relative to CMC, DP has good performance parameters for nondense and very reliable graphs, as shown in table 4. The reasonable application field of DP is networks with equal edge-failure probabilities. In that case, one simulation run results in estimating the ID of (G, T ) and serves for any q value; see section 5. 2. MP is very efficient for highly reliable networks and dense graphs; see lines 4,5 and 9,lO in table 4, and lines 3,4,5 in table 5. The performance of MP increases when the network reliability approaches 1. 3. The suggested MP algorithm needs no extra modifications to include the cases of distinct edge failure probabilities and various partition-dependent operational criteria. 4. The complexity of MP-evaluation of Q (t) for several ti ( i = 1,. .., k ) is essentially the same as for one value of t; see section 6. In particular, all 5 values of Q in table 4 can be obtained in a single simulation experiment of io4 replications; this would increase the W-performance ratio by a factor of 5. 580 zyxwvutsrqponm zyx zyxwvuts zyxwvut IEEE TRANSACTIONS ON RELIABILITY, VOL. 40, NO. 5, 1991 DECEMBER TABLE 3 Simulation Results of the Easton - Wong Network IS] (all-terminal connectivity) zyxwvutsrqponmlkj Network failure probability Edge failure probabilities pH pV WCMC e*’ PS 6?fP VMP ~ WCMC** zyxwvutsrqponmlkjihgfedcbaZ WSD 0.02 0.005 0.01 0.01 0.001 0.0005 0.0438 0.00538 0.9 1.6 <1 1.5 0.6 1.2 67 135 3.4 5.4 2.4 7.6 *)Estimated by MP **)WcMc/WsD is the performance ratio of the CMC versus the Easton-Wong’s sequential destruction method [6]. zyxwv TABLE 4 Simulation Results for the Dodecahedron (s- t connectivity; s = “1 t = “20”,see Fig. 1,a.) ‘I, (1)’ 0.5 0.2 0.1 0.05 0.02 0.71023 0.0358 0.00282 0.000288 O.oooO167 34 43 0.2 0.1 0.05 0.02 0.01 0.1876 0.0226 0.002650 0.000158 O.oooO2 1.2 1.5 1.6 1.7 <1 <1 <I 1.3 8.2 0.38 1.8 3.1 3.8 4.1 2.8 8.3 37.7 246 3472 0.67 2.0 8.8 55.7 495 1.56 (2)’ {3)* 1.81 1.91 - 0.68 1.40 2.71 - 0.56 12.3 136 - - - - - {41* 0.05 - 70.3 3714.4 - All-terminal connectivity <1 4.0 10.2 23.0 109 0.43 0.72 0.86 0.95 0.98 23 85 504 6960 51590 5.4 20 108 1227 7635 - *) Computed by MP. {I)’ Dagger sampling [16],source [9]. {2}* Sequential Construction [6],source [9]. (3)* Method of Bounds [9]. (4)’ Methods of failure sets [14],source [9].The failure set method produced two estimates, the table presents the better one, with larger value of WcMc/WFs. zyxwvutsrqponm 5. In our experiments, sparse networks were represented by the dodecahedron. Based on the results in tables 3 & 4, we suggest using MP for 0 IQ < 0.05. 6. The comparison with the results from the literature reveals that MP (for a particular value of the q-vector) is competitive with the sequential destruction method [16] applied to a very sparse network; see table 3. Based on data in [9], MP Monte Carlo considerable outperforms the dagger method and the sequential construction method when applied to the dodecahedron, for Q < 0.05; see table 4. For the same example, the MP is inferior to the Fishman method of bounds [9] by a factor of 1.5-2.5. zyxwvutsrqponmlkjihg zyxwvu zyx zyxwvutsr zyxwvutsrqpo zyxwvutsrqp I ELPERIN ET AL.: ESTIMATION OF NETWORK RELMILITY USING GRAPH EVOLUTION MODELS zyx TABLE 5 Simulation Results for Complete Graphs (all-terminal connectivity). Edge failure probability q = 0.55 Network failure probability The graph K1o 4 5 K20 K25 K30 e.) 0.456*10-’ 0.346-10-2 0.232-10-3 0.147.10-4 0.889*10-5 % wDP <1 <1 <1 <1 <1 581 a P ?UP wCMC wMP 0.54 0.57 0.53 0.50 0.47 70 887 15160 0.27.106 5.08-107 21 169 3280 47220 0.73.107 *) Estimated by MP. The dodecahedron with terminals s = 1 and t =20 (figure la) has a specific feature with respect to the method of bounds. Namely, the bound 1-A in terms of edge-disjoint cuts [9, p 149; 11, p 4631 asymptotically coincides-in this particular case-with the true network unreliability value, since the collection of cuts chosen for A contains all minimum size cuts between s and t (which is far from being so, for general G, s, t). The MP is considerable inferior to the Karp-Luby method of failure sets; see table 4. The KarpLuby method requires extra effort for computing failure sets, in terms of computer time and computer memory; see comments on this issue in [9, p 531. [9] G. S. Fishman, “A comparison of four Monte Carlo methods for estimating the probability of s-t connectedness”, IEEE Truns. Reliubility, vol R-35, 1986, pp 145-154. [lo] G. S. Fishman, “A Monte Carlo sampling plan for estimating reliability parameters and related functions”, Networks, vol17, 1987, pp 169-186. [ 111 G. S. Fishman, “Estimating the S-t reliability function using importance and stratified sampling’’,Operarions Research, ~0137,1989,pp 462-473. [12] F. Harary, Graph Theory, 1969; Addison Wesley. [13] J. E. Hopcroft, J. D. Ullman, “Setmergingalgorithms”, SIAMJ. Computing, vol2, 1973, pp 296-303. [14] R. Karp, M. G. Luby, “A new Monte Carlo method for estimating the failure probability of an n-component system”, 1983; Computer Science Div., Univ. of California. [15] J. B. Kruskal, “Onthe shortest spanning tree of a graph and the travelling salesman problem”, Proc. Amer. Mathematical Soc., vol 7, 1956, pp 48-50. [16] H. Kumamoto, et al, “Dagger sampling Monte Carlo for system unavailability evaluation”, IEEE Trans. Reliability, vol R-29, 1980 Jun, pp 122-125. [17] H. Kumamoto, T. Kazuo, I. Koichi, E. J. Henley, “State transition Monte Carlo for evaluating large, repairable systems”, IEEE Trans. Reliability, VOI R-29, 1980 DE, pp 376-380. [18] M. V. Lomonosov, V. P. Polesskii, “On maximum of probability of connectivity”, Problems of Infonnation Trunsmission, vol 8, n u n 4, 1972. [19] M. V. Lomonosov, “Bernoulli scheme with closure”, Problems of Information Transmission, vol 10, num I, 1974. [20] M. V. Lomonosov, “Tender-spot of a reliable network”, Discrete Applied Mathematics, (to appear). [21] Lucia I. P. Resende, “Implementation of factoring algorithm for reliability evaluation of undirected networks”, IEEE Truns. Reliubility, vol37, 1988 Dw, pp 462-468. [22] A. Satanarayana, R. K. Wood, “A linear time algorithm for k-terminal reliability in series-parallel networks”, SIAM J. Computing, 1985, pp 818-832. zyxwvu ACKNOWLEDGMENT We thank the referees for their constructive and helpful comments. REFERENCES [l] A. Agraval, R. E. Barlow, “A survey of network reliability and domination theory”, Operations Research, vol 32, 1984, pp 478-492. [2] R. E. Barlow, F. Proschan, Statistical Theory of Reliability and Life Testing, 1975; Holt, Rinehart & Winston. [3] C. J. Colboum, % Combinatorics ofNetwork Reliability, 1987; Oxford Univ. Press. [4] C. J. Colbourn, D. D. Harms, “Bounding all-terminal reliability in computer networks”, Networks, 1988, pp 1-12. [5] H. A. David, Order Statistics, 1981, 2d ed; John Wiley & Sons. [6] M. C. Easton, C. K. Wong, “Sequential destruction method for Monte Carlo evaluation of system reliability”, IEEE Truns. Reliability, vol R-29, 1980 Apr, pp 27-32. [7] W. Feller, An Introduction to Probability Theory and Its Applications, vol II, 1966; John Wiley & Sons. [8] G. S. Fishman, “A Monte Carlo sampling plan for estimating network reliability”, Operations Research, vol 34, 1986, pp 581-592. AUTHORS Dr. Tov I. Elperin; The Pearlstone Center for Aeronautid Engineering Studies; Dept. of Mechanical Engineering; Ben Gurion University of the Negev; FQBox 653; Beer-Sheva 84105 ISRAEL. Tov. I. ElpeM: For biography, see IEEE Truns. Reliability, vol 39, 1990 Jun, p 208. Dr.nya B. Gertsbakh; Dept. of Mathematicsand Computer Science; Ben Gurion University of the Negev; FQBox 653; Beer-Sheva 84105 ISRAEL. nya B. Gertsbakh:For biography, see IEEE Trans. Reliability, vol39, 1990 Jun, p 208. Dr. Michael V. Lomonosov; Dept. of Mathematics and Computer Science; Ben Gurion University of the Negev; POBOX653; Beer-Sheva 84105 ISRAEL. Michael V. Lomonosov is an Associate Professor in the Dept. of Mathematics and Computer Science at the Ben Gurion University of the Negev. His research interests include network reliability and network flow theory. zyxwvutsrqpo Manuscript TR89-113 received 1989 December 12; revised 1990 September 22; revised 1991 April 28. IEEE Log Number 01573 4TRF