Zdravko I Botev

Pierre L'ecuyer

Reliability Estimation for Networks with Minimal Flow Demand and Random Link Capacities

Keywords: network reliability, stochastic flow network, Conditional Monte Carlo, permutation Monte Carlo, generalized splitting

We consider a network whose links have random capacities and in which a certain target amount of flow must be carried from some source nodes to some destination nodes. Each destination node has a fixed demand that must be satisfied and each source node has a given supply. We want to estimate the unreliability of the network, defined as the probability that the network cannot carry the required amount of flow to meet the demand at all destination nodes. When this unreliability is very small, which is our main interest in this paper, standard Monte Carlo estimators become useless because failure to meet the demand is a rare event. We propose and compare two different methods to handle this situation, one based on a conditional Monte Carlo approach and the other based on generalized splitting. We find that the first is more effective when the network is highly reliable and not too large, whereas for a larger network and/or moderate reliability, the second is more effective.

Introduction

Network reliability estimation problems are commonplace in various application areas such as transportation, communication, and power distribution systems; see for example [START_REF] Gertsbakh | Models of network reliability[END_REF]. In many of those problems, the states of certain network components are subject to uncertainty and there is a set of conditions under which the network is operational, and one wishes to estimate the network unreliability, defined as the probability u that the network is in a failed state (i.e., is not operational). When u is very small, a standard (crude) Monte Carlo (MC) approach that merely generates the component states, computes the indicator function that the network is operational or not, and averages over n independent runs to estimate u, is unsatisfactory because the relative error (defined as the standard deviation of the estimator divided by the expected value u) of the MC estimator goes to infinity when u → 0.

One reliability problem that has received a lot of attention is the static network reliability estimation problem, in which each link of the network is failed with a given probability and the network is operational when a given (specific) subset of the nodes are all connected. Effective estimation methods have been developed for this problem when u is small; see [START_REF] Botev | Static network reliability estimation via generalized splitting[END_REF][START_REF] Botev | Static network reliability estimation under the Marshall-Olkin copula[END_REF][START_REF] Elperin | Estimation of network reliability using graph evolution models[END_REF][START_REF] Fishman | A Monte Carlo sampling plan for estimating network reliability[END_REF][START_REF] Gertsbakh | Models of network reliability[END_REF][START_REF] Lomonosov | Combinatorics and reliability Monte Carlo[END_REF] and the references therein. The model considered in this paper is more general. Instead of having only a binary state (up or down), each link has a random capacity that can take many possible values, there is a fixed demand that must be satisfied at certain nodes (called the destination nodes), a fixed supply is available at some other nodes (the source nodes), and the network is operational when it can carry the flow to satisfy all the demands. As a special case, there can be a single source node and a single destination node, with a fixed demand, and the network is operational when the maximum flow that can be sent from the source to the destination reaches the demand. We will describe our methods in this particular setting to simplify the notation, but the methods apply to the general setting as well. The case of links with binary states is a special case. The several methods developed for this special case do not readily apply to the network flow setting considered here, but we show how two of the best available methods for the binary case, permutation Monte Carlo (PMC) and generalized splitting (GS), can be adapted to this problem. The adaptation is not straightforward.

The PMC method [START_REF] Elperin | Estimation of network reliability using graph evolution models[END_REF][START_REF] Gertsbakh | Models of network reliability[END_REF][START_REF] Lomonosov | Combinatorics and reliability Monte Carlo[END_REF]] constructs an artificial continuous-time Markov chain (CTMC) defined as follows. Each capacity is assumed to have a discrete distribution over a finite set of possible values. This can approximate a continuous distribution if needed. We assume that all links start at their minimal capacity, and the capacity of one link may increase each time the CTMC has a jump. The CTMC is constructed so that the probability that the network is failed at time 1 is equal to u. PMC generates the discrete-time Markov chain (DTMC) underlying the CTMC, i.e., only the sequence of states that are visited until the network is operational, and conditional on that sequence it computes the probability that the network is failed at time 1, as an estimator of u. This conditional probability can be computed by exploiting the property that the failure time has a phase-type conditional distribution, whose cumulative distribution function (cdf) and density can be expressed in terms of matrix exponentials. We show how to adapt and apply the PMC principle to our problem. The CTMC construction is quite different than for the binary case. We also prove, under certain conditions, that the resulting PMC estimator has bounded relative error (BRE) when u → 0 for a given network.

GS [START_REF] Botev | Efficient Monte Carlo simulation via the generalized splitting method[END_REF] is a rare-event estimation method where the rare event is the intersection of a nested sequence of events and its probability is the product of conditional probabilities. Each conditional probability is estimated thanks to resampling strategies, making the overall estimation more accurate than a direct estimation of the rare-event probability itself. The application of GS to this problem was discussed in [START_REF] Botev | Reliability of stochastic flow networks with continuous link capacities[END_REF] for the situation in which the capacities have a continuous distribution, and experimental results were reported for a small example. But the GS algorithm proposed there does not work in general when the capacities have a discrete distribution. We show however that GS can be applied in the discrete case if we combine it with the same CTMC construction as for PMC. The GS algorithm does not have BRE in the asymptotic regime when u → 0, but it becomes more efficient that PMC when the size of the network increases. The relative error typically increases (empirically) as O(log u).

The remainder of the paper is organized as follows. In Section 2, we formulate the network flow model considered in this paper. In Section 3, we construct a CTMC which permits one to apply PMC to this model, for the case where each capacity is distributed over a finite set. In Section 4, we explain how to apply GS to this model. We report numerical experiments in Section 5. Our experimental results agree with the fact that PMC has BRE when u → 0, under appropriate conditions. It can accurately estimate extremely small values of u when the network is not too large. When the network gets larger and u is not too small, on the other hand, GS becomes more effective than PMC.

The model

Let G = (V , E ) be a graph with a set of nodes V and a set of links E with cardinality m = |E |. For i = 1, . . . , m, link i has a random integer-valued flow capacity X i with discrete marginal distribution p i (x) = P[X i = x] over the set

X i = {c i,0 , . . . , c i,b i }, 0 ≤ c i,0 < c i,1 < • • • < c i,b i < ∞.
This is a standard assumption; see [START_REF] Botev | Reliability of stochastic flow networks with continuous link capacities[END_REF] and the references therein. Thus, the random network state X = (X 1 , . . . , X m ) belongs to the space X = ∏ m i=1 X i and has joint pdf p(x) = P(X = x), for x ∈ X . We also make the standard independence assumption (see [START_REF] Alexopoulos | Capacity expansion in stochastic flow networks[END_REF][START_REF] Bulteau | A new importance sampling Monte Carlo method for a flow network reliability problem[END_REF][START_REF] Chou | An efficient and robust design optimisation of multi-state flow network for multiple commodities using generalised reliability evaluation algorithm and edge reduction method[END_REF]) that p(x) = ∏ m i=1 P[X i = x i ] and that the nodes do not fail.

To keep the notation and the exposition simple, in the remainder of the paper we describe the model and the methods under the assumption that there is a single source and a single destination. The generalization to multiple sources and destinations is straightforward, as explained below. The fixed demand level at the destination is d net > 0 and the maximum flow that can be carried from the source to the destination is a random variable Ψ(X), which is a function of the link capacities. The well-known max-flow min-cut theorem says that the maximum value of a flow from a source to a destination is equal to the minimum capacity of a cut in the network. Efficient algorithms are available to compute Ψ(X); for example the Ford-Fulkerson algorithm.

We are interested in estimating the unreliability of the flow network, defined here as

u = P[Ψ(X) < d net ] = ∑ {x∈X :Ψ(x)<d net } p(x);
that is, the probability that the maximum flow Ψ(X) fails to meet the demand. This problem was considered in [START_REF] Fishman | Monte Carlo: Concepts, algorithms, and applications[END_REF], for example. In the particular case where X i = {0, 1} for each i and d net = 1, we have an instance of the static network reliability problem mentioned in the introduction, with the source and destination as the selected set of nodes to be connected. To generalize to multiple sources and destinations, we would assume a fixed demand d i at each destination node i, a fixed supply s i at each source node i, and the event {Ψ(X) < d net } would be replaced by the event that the network does not have sufficient capacity to send flow to satisfy all the demands from the available supplies.

For small networks, it is possible to compute and store most of the minimal cutsets or pathsets and use them to obtain exact or approximate values for u; see [START_REF] Jane | Computing multi-state two-terminal reliability through critical arc states that interrupt demand[END_REF][START_REF] Zuo | An efficient method for reliability evaluation of multistate networks given all minimal path vectors[END_REF] for example. But for large networks, no polynomial-time algorithm is known for computing u exactly [START_REF] Colbourn | The combinatorics of network reliability[END_REF], and one must rely on approximations or on estimation via Monte Carlo. Of particular interest is the situation in which the network is highly reliable, i.e., u is a very small rare-event probability, because crude Monte Carlo then becomes ineffective.

Several Monte Carlo variance-reduction methods have been proposed for network reliability estimation in rare-event situations; see, e.g., [START_REF] Botev | Static network reliability estimation via generalized splitting[END_REF][START_REF] Botev | Static network reliability estimation under the Marshall-Olkin copula[END_REF][START_REF] Cancela | On the RVR simulation algorithm for network reliability evaluation[END_REF][START_REF] Elperin | Estimation of network reliability using graph evolution models[END_REF][START_REF] Gertsbakh | Models of network reliability[END_REF][START_REF] Ramirez-Marquez | A Monte-Carlo simulation approach for approximating multi-state two-terminal reliability[END_REF][START_REF] L'ecuyer | Approximate zero-variance importance sampling for static network reliability estimation[END_REF][START_REF] Tuffin | An adaptive zero-variance importance sampling approximation for static network dependability evaluation[END_REF] and the references given there. Most of these methods are for the special case of independent links with binary states and nodes that never fail. Some have been extended to links with three possible states [START_REF] Gertsbakh | Permutational methods for performance analysis of stochastic flow networks[END_REF][START_REF] Gertsbakh | Network reliability Monte Carlo with nodes subject to failure[END_REF][START_REF] Gertsbakh | Ternary networks: Reliability and Monte Carlo[END_REF], but this remains restrictive. We now describe how two of the most efficient methods, PMC and GS, can be adapted to our model.

Reformulating the model as a CTMC and applying PMC

We now show how to construct an artificial CTMC for this static model, which will permit us to apply PMC as described in the introduction. This CTMC construction differs from that used in [START_REF] Botev | Static network reliability estimation under the Marshall-Olkin copula[END_REF][START_REF] Gertsbakh | Models of network reliability[END_REF].

Constructing the CTMC

For each i, let P(X i = c i,k ) = p i (c i,k ) = r i,k > 0 for k = 0, . . . , b i . Define independent exponential random variables Y i,1 , . . . ,Y i,b i with rates λ i,1 , . . . , λ i,b i , respectively, where the λ i,k still have to be chosen. Suppose that the capacity of link i is c i,0 from time T i,0 = 0 to time T i,1 = min(Y i,1 , . . . ,Y i,b i ) (exclusive), after that it is c i,1 from time T i,1 to time T i,2 = min(Y i,2 , . . . ,Y i,b i ), it is c i,2 from time T i,2 to time T i,3 = min(Y i,3 , . . . ,Y i,b i ), and so on, and finally it is

c i,b i from time T i,b i to T i,b i +1 = ∞.
Under this process, the capacity of link i at time γ ≥ 0 is given by

X i (γ) = c i,k for T i,k ≤ γ < T i,k+1 and 0 ≤ k ≤ b i (1) = max k {c i,k : T i,k ≤ γ}. (2)
The times T i,1 , . . . , T i,b i are not necessarily all distinct; often, many of them are equal, so that the number of jumps at which the capacity changes can be much smaller than b i . For example, if

T i,1 = Y i,b i , then we have T i,1 = T i,2 = • • • = T i,b i . As another example, if b i = 3 and Y i,2 < Y i,1 < Y i,3 , then 0 < T i,1 = T i,2 < T i,3
and the capacity of link i jumps from c i,0 to c i,2 at time T i,2 = Y i,2 and jumps again from c i,2 to c i,3 at time T i,3 = Y i,3 . In general, the process {X i (γ), γ ≥ 0} has an upward jump at each of the distinct jump times T i,k .

To show that this process is a CTMC, suppose that we are at time γ ≥ 0 and X i (γ) = c i,k . Then we know that Y i,k ≤ γ and that Y i, > γ for all > k. The Y i, for < k can be anything, but they have no influence on the process trajectory after time γ. This means that the current state X i (γ) contains all the relevant information that needs to be known at time γ to generate the future of the process.

The capacity X i (γ) of link i at time γ ≥ 0 satisfies

P[X i (γ) ≤ c i,k ] = P[min(Y i,k+1 , . . . ,Y i,b i ) > γ] = exp[-γ(λ i,k+1 + • • • + λ i,b i )].
If we select the λ i,k 's so that the last expression equals r i,0 + • • • + r i,k for each k when γ = 1, then X i (1) has the exact same distribution as X i , the capacity of link i in the original static model. This is equivalent to having

λ i,k+1 + • • • + λ i,b i = -ln(r i,0 + • • • + r i,k ).
To achieve this, it suffices to put

λ i,b i = -ln(r i,0 + • • • + r i,b i -1 ) = -ln(1 -r i,b i ) (3)
and then for k = b i -1, . . . , 1 (in descending succession):

λ i,k = -ln(r i,0 + • • • + r i,k-1 ) -λ i,k+1 -• • • -λ i,b i . (4)
Note that (4) can be rewritten as

λ i,k = -ln(r i,0 + • • • + r i,k-1 ) + ln(r i,0 + • • • + r i,k
), which can never be negative. We have proved the following.

Proposition 1. If we select λ i,b i , λ i,b i -1 , .
. . , λ i,1 according to ( 3) and ( 4) and the process X i (•) as in ( 2), then λ i,k ≥ 0 for each k, {X i (γ), γ ≥ 0} is a CTMC process, and X i (1) has exactly the same discrete distribution as the capacity of link i in the original model:

P[X i (1) = c i,k ] = r i,k for k = 0, . . . , b i . As a result, X(1) = (X 1 (1), . . . , X m (1)
) has the same distribution as X and one has u

= P[Ψ(X(1)) < d net ].

Applying PMC

Under the assumption that all links are independent, a simple way of applying PMC to this model is as follows. Generate all the Y i,k 's independently with their rates λ i,k , put them in a large vector

Y = (Y 1,1 , . . . ,Y 1,b 1 , . . . ,Y m,1 , . . . ,Y m,b m ) of size κ = b 1 + • • • + b m
, and sort this vector in increasing order to obtain

Y π(1) ≤ Y π(2) ≤ • • • ≤ Y π(κ) ,
where

π( j) = (i, k) if Y i,k
is in position j in the sorted vector, so that π = (π(1), . . . , π(κ)) can be seen as the permutation of the pairs (i, k) that corresponds to the sort. This permutation gives an ordering of the κ pairs (i, k). When scanning those pairs in the given order, each pair (i, k) corresponds to a potential capacity increase for link i. The capacity increases if and only if no pair (i, k ) for k > k has occurred before. Conditional on π, one can add those pairs in the given order and update the capacities accordingly, until the maximum flow in the network reaches d net . Suppose this occurs when adding the pair

(i, k) = π(C) for some integer C > 0. Let T C = Y π(C) .
The (unbiased) conditional (PMC) estimator of u is then

P[Ψ(X(1)) < d net | π] = P[T C > 1 | π] = P[T C > 1 | π(1), . . . , π(C)] = P[A 1 + • • • + A C > 1],
where

A 1 = Y π(1) is an exponential random variable with rate Λ 1 = ∑ m i=1 ∑ b i k=1 λ i,k , each A j = Y π( j) -Y π( j-1
) is an exponential random variable with rate Λ j = Λ j-1λ π( j-1) for j = 2, . . . ,C, and these A j 's are independent. Given π and C, T C = A 1 + • • • + A C is the sum of C independent exponential random variables with rates Λ 1 , . . . , Λ C , which has a phase-type distribution, whose complementary cdf is given by

1 -F(γ | π) = P[T C > γ | π] = e t 1 exp(Q γ)1 (5)
where e t 1 = (1, 0, . . . , 0), 1 = (1, . . . , 1) t (the t means "transposed"), and

Q =         -Λ 1 Λ 1 0 • • • 0 0 -Λ 2 Λ 2 0 . . . 0 0 . . . . . . 0 . . . 0 0 -Λ C-1 Λ C-1 0 • • • 0 0 -Λ C         .
Reliable and fast computation of ( 5) is discussed in [START_REF] Botev | Static network reliability estimation via generalized splitting[END_REF][START_REF] Botev | Static network reliability estimation under the Marshall-Olkin copula[END_REF].

To compute the critical number C at which the flow reaches the demand, we must be able to update efficiently the maximum flow in the network each time we increase the capacity of one link. We do this as explained in Section 4 of [START_REF] Botev | Reliability of stochastic flow networks with continuous link capacities[END_REF]. We refer to this algorithm as the incremental maximum flow algorithm.

To estimate u by PMC, for a fixed threshold d net , we simulate n independent realizations

W 1 , . . . ,W n of W = W (π) = P[T C > 1 | π] (6)
and take the average Wn = (1/n) ∑ n i=1 W i . Compared with the crude Monte Carlo estimator that would take the indicator

I = I[Ψ(X(1)) < d net ] in place of W , it is always true that Var[W ] < Var[I], because W = E[I | π].
The estimators discussed so far are for a single (fixed) demand d net . With PMC, it is also possible to estimate u = u(d net ) as a function of the demand d net , over some interval, using the same simulations for all demands. To do this, for any given permutation π, we can compute C = C(d net ) as a function of the demand over the interval of interest. This would be a step function, often with just a few jumps. Then we compute W for all values of C that are visited over this interval. This provides an estimator W (d net ) of u(d net ) as a function of d net . By averaging the n realizations of this estimator, we obtain a functional estimator of u(d net ) over the interval of interest.

Improved PMC

The PMC strategy described earlier can be improved by removing some useless jumps. First, whenever c i,k ≥ d net for k < b i , we can immediately remove all the jumps (i, k + 1), . . . , (i, b i ), because when the capacity of a link has reached d net , it is useless to increase it further. Capacity levels larger than d net can in fact be all reset to d net right away in the model, the probability of d net in the new model being taken as the the probability of values larger than d net in the initial model. For simplicity, when the demand is fixed, we assume in our algorithm that this has been done already, so that c i,b i ≤ d net for all i, and then there is no need to remove those useless capacity levels.

Second, the jump times Y i,k that do not change the capacity of link i can also be removed. That is, whenever π( j) = (i, k) and the capacity of link i has already reached a value c i,k > c i,k , i.e., (i, k ) = π( j ) for some j < j, then there is no need to consider the pair (i, k) when it is encountered in the permutation, so we can remove the corresponding jump.

Let π be the permutation obtained after removing all those pairs (i, k) from the sorted vector, and C the corresponding value of C in this reduced permutation. As soon as the max flow reaches d net , we have found C. When we encounter π( j) = (i, k) and the previous capacity of link i was c i,k < c i,k , the capacity of link i jumps to c i,k and we must decrease

Λ j by λi,k = λ i,k +1 + • • • + λ i,k ,
because the jumps that correspond to (i, k + 1), . . . , (i, k) can now be removed from consideration. draw Y i,1 , . . . ,Y i,b i with the appropriate rates λ i,k

3:

let L i = {(i, 1), . . . , (i, b i )} 4: min ← (i, b i ) and λi,b i ← λ i,b i 5: for k = b i -1 to 1 do 6: if Y i,k > Y min then 7:
remove (i, k) from the list L i 8:

λmin ← λmin + λ i,k

9:

S i,k ← 0 // this jump is deactivated S i,k ← 1 // this jump is activated 14: merge the sorted lists L 1 , . . . , L m into a single list sorted by increasing order, Y π(1) , . . . ,Y π( κ) 15:

Λ 1 ← λ 1,1 + • • • + λ 1,b 1 + • • • + λ m,b m 16: j ← 0 17: X ← (c 1,0 , . . . , c m,0 ) 18: while maximum flow Ψ(X) < d net do 19: j ← j + 1 20: if S π( j) = 1 then 21: (i, k) ← π( j)
// this jump has not been removed or executed 22:

Λ j+1 ← Λ j -λi,k 23: S i,k ← 0 24: X i ← c i,k // increase capacity of i-th link 25:
Filter() // do nothing (default), or FilterSingle, FilterAll, etc. 26: C ← j // the critical jump number 27: return

W ← P A 1 + • • • + A k-1 > 1 | π, C .
Algorithm 1 describes this reduced version of PMC in a more formal way. It returns one realization of W . Indentation delimits the scope of the loops. In the first for loop, for each link i, the algorithm generates the exponential random variables Y i,k and then immediately eliminates those that correspond to (useless) jump times at which the capacity of the link does not change. This preliminary filtering is very easy and efficient to apply and may eliminate a significant fraction of the jumps, especially for links that have many capacity levels. The remaining jumps are sorted in a single list (for all links) and each one receives a Boolean tag S π( j) , initialized to 1, which means that this jump is currently scheduled to occur.

Then these jumps are "executed" in chronological order, by increasing the corresponding capacities, until the critical jump number C is found. After that, W can be computed. The Boolean variables S π( j) are used in the optional Filter() subroutine, which can be used to try to eliminate further useless jumps after a jump is executed and the corresponding capacity is increased (this is discussed in Section 3.4). Algorithm 1 would be invoked n times, independently, and u would be estimated by the average of the n realizations of W .

Other variants of the algorithm can be considered and some might be more efficient, but this is not completely clear. For example, instead of generating all the variables Y i,k at the beginning, one may think of generating the permutation π directly without generating those Y i,k , as was done in [START_REF] Gertsbakh | Models of network reliability[END_REF] for the binary case. This appears complicated and we did not implement it.

Removing jumps having no impact on maximum flow

In Algorithm 1, in the case where Filter() does nothing, all pairs (i, k) for which the capacity of link i increases are retained in π and the corresponding jumps are executed. But it sometimes occurs that increasing the capacity of link i to c i,k (or more) is useless because it can no longer have an impact on the event that the maximum flow exceeds the demand or not. In this case, one can cancel (deactivate) all the future jumps related to the capacity of link i. In our implementation, these future jumps are canceled by setting their Boolean variables S i,k to 0. Increasing the capacity of link i is useless in particular if it is already possible to send d net units of flow between the two nodes connected by link i. This obviously happens if the capacity of link i is already of d net , which is trivial to verify, but under our assumption (made at the beginning of Section 3.3) that a link has no capacity level above d net , this cannot happen, and our algorithm ignores this possibility.

Increasing the capacity of link i is also useless when d net units of flow can be sent in total, either directly on link i or indirectly via other links. This is generally harder (more costly) to verify. To detect it, one can run a max-flow algorithm to compute how much flow can be sent between these two nodes. This can be done each time the capacity of a link is increased. Algorithm 2 does this only for the link i whose capacity has just been increased, at each step j. Since the link i generally changes at every step j, we have a different max-flow problem (for a different pair of nodes) at each step. For this reason, in our implementation we recompute the max-flow from scratch at each step j. Of course, this brings significant overhead. Algorithm 3 is even more ambitious: it computes the max-flow between nodes for all pairs of nodes. Then for each link i for which the current max-flow between the corresponding two nodes meets the demand, it cancels all future jumps associated with that link. Doing this at each step j might be too costly, so in our implementation the user selects a positive integer ν and does it only at every ν steps, i.e., when j is a multiple of ν. We compute the max-flow for all pairs of nodes using the algorithm of [START_REF] Gusfield | Very simple methods for all pairs network flow analysis[END_REF], which is a simplified variant of the Gomory-Hu method [START_REF] Gomory | Multi-terminal network flows[END_REF]. This algorithm computes the all-pairs max-flow by applying |V | -1 times a max-flow algorithm for one pair of nodes, which is generally more efficient than applying a max-flow algorithm for each link. Algorithm 3 also recomputes the maximum flows from scratch each time it is called, rather than reusing computations from the previous time and just updating the max flows. (In fact, the |V | -1 pairs of nodes for which the max-flow is computed in the algorithm change from one call to the next, and we are not aware of an effective incremental algorithm that would reuse and just update the previous computations.)

Algorithm 2 : FilterSingle f ← compute maximum flow between terminal nodes of link i if f ≥ d net then for k = 1 to b i do if S i,k = 1 then S i,k ← 0 Λ j+1 ← Λ j+1 -λi,k
Algorithm 3 : FilterAll if j mod ν = 0 then {F v,w } ← max flow between all pairs of nodes (v, w), computed via Gusfield's algorithm

for all i = (v, w) ∈ E do if F v,w ≥ d net then for k = 1 to b i do if S i,k = 1 then S i,k ← 0 Λ j+1 ← Λ j+1 -λi,k

Bounded relative error for PMC

For a single run, the crude MC estimator

I = I[Ψ(X(1)) < d net ] of u, which is a Bernoulli random variable with mean u, has variance u(1 -u), so its relative error (RE) is RE[I] = u(1 -u)/u ≈ u -1/2 → ∞ when u → 0.
With n runs, the variance is divided by n and the RE by √ n. When u is very small, we may need an excessively large n to obtain a sufficiently small RE. With PMC, the RE is sometimes much better behaved than with MC. In this section, we obtain conditions under which the PMC estimator W has bounded relative error (BRE), i.e., RE[W ] remains bounded when u → 0. The proofs have some similarity with those in [START_REF] Botev | Static network reliability estimation under the Marshall-Olkin copula[END_REF].

Suppose the probabilities r i,k = r i,k (ε) in our model depend on some parameter ε in a way that u = u(ε) → 0 when ε → 0. In what follows, the quantities in the model are assumed implicitly to depend on ε. A non-negative quantity that may depend on ε is O(1) if it remains bounded when ε → 0. It is Θ(1) if it is bounded and also bounded away from 0, when ε → 0.

In our setting, the vector Y and the permutation π have finite length κ and C is bounded by κ. The number of possible permutations is therefore finite. Let p(π) be the probability of permutation π. Proposition 2. (i) If p(π) = Θ(1) for all π, then the PMC estimator has BRE.

(ii) This holds in particular if λ i,k i /λ j,k j = Θ(1) for all i, j, k i , k j (we then say that the rates are balanced).

Proof. 1), which implies BRE.

(i) Note that u ≥ P[T C > 1 | π]p(π) = W (π)p(π) for any π. If p(π) = Θ(1), then W (π)/u = O(1) and max π W (π)/u = O(1). Therefore E[W 2 /u 2 ] = O(
(ii) Note that p(π) = ∏ κ j=1 λ π( j) /Λ π( j) . Under the given assumption, λ π( j) /Λ π( j) = Θ(1) for all j, which implies that p(π) = Θ(1) for all π.

As a concrete illustration of an asymptotic regime in which the r i,k depend on ε, we define a regime similar to one that has been widely used for highly reliable Markovian systems [START_REF] Nakayama | General conditions for bounded relative error in simulations of highly reliable Markovian systems[END_REF][START_REF] Rubino | Markovian models for dependability analysis[END_REF][START_REF] Shahabuddin | Importance sampling for the simulation of highly reliable Markovian systems[END_REF]. Suppose that link i is operating at capacity c i,k . with probability

r i,k = a i,k ε d i,k
for some constants a i,k > 0 and d i,k > 0 independent of ε, for all k ∈ {0, . . . , b i -1} (that is, not at full capacity). This implies that r i,b i = 1 -∑ b i -1 j=0 r i,k = Θ(1) for all i. That is, the event that any link is not at full capacity is a rare event. This implies that failure to meet the demand is a rare event, and therefore RE[I] → ∞ when ε → 0. More specifically, any state vector x = (c 1,k 1 , . . . , c m,k m ) for which Ψ(x) < d net has probability P(X

= x) = ∏ m i=1 r i,k i = ε d(x) (1 + o(1)) for some d(x) > 0. Let d min = min{d(x) : Ψ(x) < d net }. Then u = ∑ {x:Ψ(x)<d net } P(X = x) = Θ(ε d min ).
On the other hand, we have Proposition 3. In the setting just defined, the PMC estimator has BRE.

Proof. Recall that for 1 ≤ i ≤ m, λ i,b i = -ln(∑ b i -1 k=0 r i,k ) = Θ(ln(ε)). Moreover, for all k < b i , λ i,k = ln(∑ k =0 r i, )-ln(∑ k-1 =0 r i, ) = Θ(ln(ε)).
The conditions of Proposition 2 (ii) are then verified, hence the result.

The results of this section apply to the improved PMC variants as well; the proofs are easily adapted.

A generalized splitting algorithm

Botev et al. [START_REF] Botev | Reliability of stochastic flow networks with continuous link capacities[END_REF] have explained how to adapt the GS algorithm proposed and studied in [START_REF] Botev | Efficient Monte Carlo simulation via the generalized splitting method[END_REF][START_REF] Botev | Static network reliability estimation via generalized splitting[END_REF] to the stochastic flow problem considered here, but for the situation where the capacities have a continuous distribution. The aim of the algorithm is to obtain a sample of realizations of X which is approximately a sample from the distribution of X conditional on Ψ(X) < d net . The estimator is then given by the realized sample size (which is random) divided by its largest possible value. The algorithm uses intermediate demand levels

d net = d τ < • • • < d 1 < d 0 ,
where d 0 is the maximal possible flow, achieved when each link i is at its maximal capacity c i,b i . These levels and their number τ are fixed a priori and chosen so that P[Ψ(X) < d t | Ψ(X) < d t-1 ] ≈ 1/s for t = 1, . . . , τ -1, and at most 1/s for t = τ, where s is a small integer also fixed (usually and in all our experiments in this paper, s = 2). The levels are estimated by pilot runs, as explained in [START_REF] Botev | Efficient Monte Carlo simulation via the generalized splitting method[END_REF][START_REF] Botev | Static network reliability estimation via generalized splitting[END_REF]. The algorithm starts by sampling X from its original distribution. If Ψ(X) < d 1 , it resamples each coordinate of X conditional on Ψ(X) < d 1 , via Gibbs sampling, repeats this s times, and keeps the states X for which Ψ(X) < d 2 (their number is in {0, . . . , s}). At each level t = 3, . . . , τ, this type of resampling is applied to each state that has been retained at the previous step (for which Ψ(X) < d t-1 ), by resampling that state twice from its distribution conditional on Ψ(X) < d t-1 , and retaining the states for which Ψ(X) < d t . At the last level, we count the number N of chains for which Ψ(X) < d τ , and return W = N/s τ-1 as an estimator of u. This is repeated n times independently, to produce n independent realizations of W , say W 1 , . . . ,W n , whose average Wn is an unbiased estimator of u. This estimator does not have BRE, because the RE increases with the number of levels; the RE is typically (roughly) proportional tolog u (see [START_REF] Botev | Dependent failures in highly-reliable static networks[END_REF] for a proof of this result in an idealized setting). It can also handle large networks.

In general, this GS algorithm is not directly applicable when the capacities have discrete distributions, because then Ψ(X) also has a discrete distribution and it may happen that this distribution is too coarse (e.g., all the probability mass is on just a few possible values). Then it may be impossible to select levels d t for which P[Ψ(X)

< d t | Ψ(X) < d t-1 ] ≈ 1/s for t = 1, . . . , τ -1.
It is nevertheless possible to apply GS in that case by constructing the vector Y as for PMC in the previous section, and resampling this vector instead of X. Recall that Y is a vector of κ independent exponential random variables. The GS algorithm will operate similarly as the one described above, except that now the levels 0 = γ 0 < γ 1 < • • • < γ τ = 1 are on T C , and the resampling at each step t is for Y and is conditional on T C > γ t-1 . This is valid because {Ψ(X(γ t )) < d net } = {T C > γ t }. The corresponding GS procedure operates in the same way as the GS algorithm with anti-shocks in [START_REF] Botev | Static network reliability estimation under the Marshall-Olkin copula[END_REF]. We first generate Y from its original distribution. Then at each level t, we take each state (realization or modification of Y) that has been retained at the previous step (for which T C > γ t-1 ), we resample all its coordinates s times (i.e., for s Gibbs sampling steps, where each step starts from the result of the previous step) from its distribution conditional on T C > γ t-1 , to obtain two new states, and we retain the states for which T C > γ t . At the last level, we count the number N of chains for which T C > γ τ = 1, and return W = N/s τ-1 as an estimator of u. The resampling of Y conditional on T C > γ t-1 via Gibbs sampling can be done in a similar way as in [START_REF] Botev | Static network reliability estimation via generalized splitting[END_REF][START_REF] Botev | Static network reliability estimation under the Marshall-Olkin copula[END_REF]. We first select a permutation π of the κ coordinates of the vector Y. Then for j = 1, . . . , κ, we resample Y π( j) as follows: If π( j) = (i, k), the current capacity X i (γ t-1 ) of link i is less than c i,k (or equivalently min(Y i,k , . . . ,Y i,b i ) > γ t-1 ), and by changing the current capacity of link i to c i,k (or equivalently changing Y i,k to 0) we would have T C < γ t-1 (the maximum flow would meet the demand), then we resample Y i,k from its exponential density truncated to (γ t-1 , ∞). Otherwise we resample Y i,k from its original exponential density. To sample from the truncated density, it suffices to generate Y i,k from the original density and add γ t-1 .

Numerical Examples

In this section, we provide some numerical examples that compare the PMC and GS algorithms, and show how they behave when u → 0. In these examples, we parameterize the models by ε in a way that u = u(ε) → 0 when ε → 0, exactly as in Section 3.5, in the asymptotic regime when the probability that links are not operating at full capacity is getting close to zero. For all variants of PMC, we used formula (2) in [START_REF] Botev | Static network reliability estimation under the Marshall-Olkin copula[END_REF] with high precision arithmetic to compute (5).

Experimental setting

We used the same experimental protocol as in [START_REF] Botev | Static network reliability estimation under the Marshall-Olkin copula[END_REF], comparing four methods. Method PMC refers to Algorithm 1 without any filtering step. PMC-Single and PMC-All refer to PMC combined with the filtering as in Algorithms 2 and 3 respectively. Method GS refers to generalized splitting, implemented as described in Section 4. The splitting levels were determined via the adaptive Algorithm 3 of [START_REF] Botev | Static network reliability estimation via generalized splitting[END_REF], with n 0 = 500 and s = 2. The levels were estimated using a single run of the adapative algorithm, and these same levels were used for every independent replication of the GS algorithm.

For each example and method, we report the unreliability estimate Wn , its empirical relative error RE[ Wn ] = S n /( √ n Wn ) where S 2 n is the empirical variance, and the work-normalized relative variance (WNRV) of Wn , defined as WNRV[ Wn ] = T × RE 2 [ Wn ], where T is the total CPU time (in seconds) for the n runs of the algorithm. One must keep in mind that T and the WNRV depend on the software and hardware used for the computations. The experiments were run on Intel Xeon E5-2680 CPUs, on a linux cluster. The sample size for every algorithm was n = 5 × 10 4 .

For each example we use the following model. Each link i has the capacity levels {0, 1, . . . , b i }, i.e., c i,k = k for k = 0, . . . , b i . We take r

i,k = P(X i = c i,k ) = ρ b i -k-1 ε for k < b i and r i,b i = 1 - ∑ b i -1 k=0 ρ b i -k-1 ε
, where ρ, ε and {b i } are model parameters.

A 4 × 4 lattice graph

Our first example uses the 4 × 4 lattice graph, which has 16 nodes and 24 links. The flow has to be sent from one corner to the opposite corner. We take b i = 8, ρ = 0.6 and d net = 10, and let ε range from 10 -4 to 10 -13 . Table 1 reports the values of Wn , RE[ Wn ] and WNRV[ Wn ], for methods GS and PMC, for some values of ε. Figure 1 shows plots of RE[ Wn ] and WNRV[ Wn ] for all four methods. We see that for this small example, PMC-All always has the smallest RE, followed by PMC-Single. These REs increase very slowly when ε decreases, for the values we have tried. They should eventually stabilize when ε → 0. In terms of work-normalized relative variance, i.e., when taking the computing time into account, GS is the most effective method when ε is not too small, but when ε gets smaller, GS requires a larger simulation effort while the effort required by PMC variants remains approximately stable, so these PMC methods eventually catch up, in agreement with our asymptotic results. Among them, PMC-Single has the smallest WNRV. ε considered, PMC-All has the smallest RE while GS wins in terms of WNRV.

ε = 10 -4 ε = 10 -5 ε = 10 -6 ε = 10 -7 ε = 10 -

A dodecahedron network

In this example we use the well-known dodecahedron network (Figure 3), with 20 nodes and 30 links, often used as a standard benchmark in network reliability estimation [START_REF] Botev | Static network reliability estimation via generalized splitting[END_REF][START_REF] Cancela | A recursive variance-reduction algorithm for estimating communication-network reliability[END_REF][START_REF] Cancela | Rare event analysis by Monte Carlo techniques in static models[END_REF][START_REF] Cancela | Analysis and improvements of path-based methods for Monte Carlo reliability evaluation of static models[END_REF][START_REF] Tuffin | An adaptive zero-variance importance sampling approximation for static network dependability evaluation[END_REF]. Here we took ρ = 0.7, b i = 4 and d net = 5. Note that when ε is very small, most of the failures will occur because there is not enough capacity in the three links connected to node 1 (links 1, 2, 3), or not enough capacity in the three links connected to node 20 (links [START_REF] Rubino | Markovian models for dependability analysis[END_REF][START_REF] Shahabuddin | Importance sampling for the simulation of highly reliable Markovian systems[END_REF][START_REF] Tuffin | An adaptive zero-variance importance sampling approximation for static network dependability evaluation[END_REF]. These are the two bottleneck cuts.

Table 2 reports the values of Wn , RE[ Wn ], and WNRV[ Wn ], for the GS and PMC methods, for different values of ε. We see that the estimates Wn agree very well across the two methods. Figure 4 shows RE[ Wn ] and WNRV[ Wn ] as functions of ε, for all four methods. We see that PMC-All has by far the smallest RE for all ε, and it also wins in terms of WNRV, except for ε > 10 -5 where GS wins. The latter case is approximately when u ≥ 7 × 10 -11 , which is already pretty small. When ε decreases, the WNRV increases for GS in part because the RE increases, but also because the computing time increases. The figure shows what happens when ε gets very small.

Algorithm 1 :

 1 PMC algorithm for multi-state flow network 1: for i = 1 to m do 2:

Figure 1 :

 1 Figure 1: RE (left) and WNRV (right) for four methods, for the 4 × 4 lattice graph.

Figure 2 :

 2 Figure 2: RE (left) and WNRV (right) for four methods, for the 6 × 6 lattice graph.

Figure 3 :

 3 Figure 3: A dodecahedron graph with 20 nodes and 30 links (figure taken from [3]).

10 - 3

 3 10 -5 10 -7 10 -9 10 -11 10 -13 10 -15 10 -17 10 -5 10 -7 10 -9 10 -11 10 -13 10 -15 10 -

Figure 4 :

 4 Figure 4: RE (left) and WNRV (right) for four methods, for the dodecahedron example

8

  Wn for PMC 2.99 × 10 -5 2.98 × 10 -6 2.99 × 10 -7 2.99 × 10 -8 2.99 × 10 -9 RE[ Wn ] for PMC 3.16 × 10 -2 3.34 × 10 -2 3.41 × 10 -2 3.69 × 10 -2 3.74 × 10 -2 WNRV[ Wn ] for PMC 3.17 × 10 -2 3.20 × 10 -2 3.17 × 10 -2 3.62 × 10 -2 3.54 × 10 -2 Wn for GS 2.98 × 10 -5 2.99 × 10 -6 2.99 × 10 -7 2.98 × 10 -8 2.99 × 10 -9 RE[ Wn ] for GS 3.43 × 10 -2 3.32 × 10 -2 3.15 × 10 -2 3.30 × 10 -2 4.33 × 10 -2 WNRV[ Wn ] for GS 1.07 × 10 -2 1.43 × 10 -2 1.68 × 10 -2 2.35 × 10 -2 2.86 × 10 -2

Table 1 :

 1 Estimation of u, RE, and WNRV for some values of ε, for the 4 × 4 lattice example 5.3 6 × 6 lattice graph Figure 2 shows plots of RE[ Wn ] and WNRV[ Wn ] for a 6 × 6 lattice graph, with 36 nodes and 60 links. Again, the flow has to be sent from one corner to the opposite corner. Here, for all values of

		10 -1.2	Method	GS
				PMC PMC-All				10 -0.4
	Relative Error	10 -1.6 10 -1.4		PMC-Single			Work Normalized Relative Variance	10 -1 10 -0.8 10 -0.6	Method	GS PMC
								10 -1.2	PMC-All PMC-Single
			10 -6	10 -8	10 -10	10 -12	10 -6	10 -8	10 -10	10 -12

  7.06 × 10 -13 7.06 × 10 -15 7.05 × 10 -17 RE[ Wn ] for PMC 8.63 × 10 -2 7.21 × 10 -2 6.68 × 10 -2 5.97 × 10 -2 5.86 × 10 -2 WNRV[ Wn ] for PMC 1.85 × 10 -1 1.37 × 10 -1 1.14 × 10 -1 8.95 × 10 -2 8.57 × 10 -2 Wn for GS 7.07 × 10 -9 7.06 × 10 -11 7.07 × 10 -13 7.07 × 10 -15 7.05 × 10 -17 RE[ Wn ] for GS 3.95 × 10 -2 4.30 × 10 -2 4.58 × 10 -2 5.17 × 10 -2 4.97 × 10 -2 WNRV[ Wn ] for GS 3.13 × 10 -2 4.52 × 10 -2 6.00 × 10 -2 8.81 × 10 -2 1.06 × 10 -1

Table 2 :

 2 Estimation of u, RE, and WNRV for some values of ε, for the dodecahedron example

Acknowledgments

This work has been supported by the Australian Research Council under DE140100993 Grant to Z. I. Botev, an NSERC-Canada Discovery Grant, a Canada Research Chair, and an Inria International Chair to P. L'Ecuyer. P. L'Ecuyer acknowledges the support of the Faculty of Science Visiting Researcher Award at UNSW. We are grateful to Rohan Shah, who performed the numerical experiments.