# Solving implementation issues of Permutation Monte Carlo in different network reliability contexts

Héctor Cancela, Leslie Murray, Gerardo Rubino

► To cite this version:

Héctor Cancela, Leslie Murray, Gerardo Rubino. Solving implementation issues of Permutation Monte Carlo in different network reliability contexts. CLEI Electronic Journal., 2025, 28 (3), 10.19153/cleiej.28.3.8. hal-05451049

HAL Id: hal-05451049

https://inria.hal.science/hal-05451049v1

Submitted on 9 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Solving implementation issues
of Permutation Monte Carlo
in different network reliability contexts *

**Héctor Cancela**

Facultad de Ingeniería, Universidad de la República
Montevideo, Uruguay
*cancela@fing.edu.uy

**Leslie Murray**

FCEIA, Universidad Nacional de Rosario
Rosario. Argentina
*lmurray@fceia.unr.edu.ar

**Gerardo Rubino**

INRIA-Rennes, Bretagne Atlantique, Campus de Beaulieu
Rennes, France
*gerardo.rubino@inria.fr

January 2025

<footnote>*This research was partially funded by the PEDECIBA program (Programa de Desarrollo de las Ciencias Básicas - Uruguay), and by Math-AMSUD project 19-MATH-03 RareDep.</footnote>
Abstract

Network reliability computation is an NP-hard problem which has attracted much attention in literature. This problem consists in, given a network where the links may fail or operate with known probabilities, to compute the probability that a given subset of nodes (known as terminals) are connected by the operational links. Given the difficulty to compute the exact value of the network reliability, an alternative which has been much explored in the literature is the use of Monte Carlo estimation methods. In this work we discuss Permutation Monte Carlo, a highly efficient network reliability estimation method. The method is prone to numerical limitations for networks with a high number of links. We discuss this situation and present a simple way to rewrite the algorithm's calculations to make it numerically more stable. We also present a variant of Permutation Monte Carlo implementation for the particular case where all the links in the network under study share the same failure probability distribution (homogeneous network). For this type of network the Permutation Monte Carlo implementation can be redesigned to be much more efficient. We present some experimental test results proving that both proposed implementation variants are extremely efficient.
# Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>2</td><td>Network Static Model</td><td>5</td></tr><tr><td>3</td><td>Dynamic Models</td><td>6</td></tr><tr><td>4</td><td>Permutation Monte Carlo (PMC)</td><td>7</td></tr><tr><td>5</td><td>PMC Implementation Issues</td><td>10</td></tr><tr><td>6</td><td>Alternative Implementation for PMC</td><td>13</td></tr><tr><td>7</td><td>PMC on Network Homogeneous Models</td><td>15</td></tr><tr><td>8</td><td>F-Monte Carlo</td><td>16</td></tr><tr><td>9</td><td>Computational results</td><td>18</td></tr><tr><td>9.1</td><td>PMC Implementation Improvements</td><td>18</td></tr><tr><td>9.2</td><td>PMC and F-Monte Carlo over Homogeneous Models</td><td>23</td></tr><tr><td>10</td><td>Conclusions and future work</td><td>25</td></tr></table>

# 1 Introduction

In recent years, network models have been an important support for the study of a wide variety of problems. The dynamic or static nature, the size, topology and type of connectivity, as well as the different definitions of reliability, make each network model appropriate for a different context.

The present work is based on the so called Network Static Model, a model designed for the analysis of connectivity problems. In this model, the network must guarantee the existence of paths to allow the connection between a set of nodes called terminal nodes. The links of the network are subject to failure, while the nodes are supposed to be perfect (this model can also be extended to include nodes failures). The failure of a link consists of its total interruption and occurs with known probability. Then, there is a probability —called reliability— that the links not failed form paths that guarantee the connection between the terminal nodes.

The computation of the reliability thus defined belongs to the class of NP-hard problem (see [36] and [35]). As all known algorithms have exponential complexity, for
medium to large sized networks it is impossible to compute this measure in reasonable
times [25, 39].

One proposal to avoid this problem is to build an estimate, using Monte Carlo meth-
ods, instead of computing the exact value of network reliability. However, as the links'
reliability grow up to very high values and, therefore, the whole network becomes more
reliable, the most direct Monte Carlo simulation approach, which is Standard Monte
Carlo (Standard MC), becomes less efficient (see for instance [38] for a discussion). This
is because, in a high reliability scenario, almost no replication will find the network in
a state of failure while, in almost all of them the network will be in a perfect working
condition. This causes an increase in the variance of the estimator. We can also justify
this problem more formally as follows.

A parameter that is commonly used to evaluate the efficiency of a Monte Carlo
estimate, a sort of relative error, is the Coefficient of Variation (CV), defined as the
standard deviation divided by the estimated value. Suppose that the unreliability of
the network under study (an extremely low value for a highly reliable network) is Q,
then: $CV = (Q(1 - Q))^{1/2}/Q$. Note that the denominator goes to zero faster than the
numerator, causing the quotient to grow boundless as Q decreases.

Over the years, a large amount of research has aimed to obtain reliability estimators
with lower variance than that of the Standard MC estimator on the Static Network
Model. There is a huge literature on this subject, ranging from methods exploiting
topological or probabilistic properties such as [26], [27], [14], [37], [6], [7] and [28] up
to methods which first transform the original problem into a dynamic model and then
perform the computations in this modified setting, such as [10], [11], [31], [23], [24], [33].
Also many papers discuss other network variants, including aspects such as distance
constraints [4] and dependencies between states [3]. Some state of the art papers can be
found at [13], [8], and [34]. Nevertheless, none of these works have addressed the stability
of the methods from a computational, numerical analysis viewpoint (as opposed as from
a probabilistic/statistical viewpoint).

One of the methods with good statistical properties is Permutation Monte Carlo,
proposed in [11]. Nevertheless, the direct implementation of this method runs into
numerical problems for large networks, which has precluded its application in these
contexts.

In this work, we discuss what is the source of this numerical instability, and we
propose a simple way to solve it by rewriting a key algorithmic component of the method.
At the end of the article we include some experimental tests that demonstrates the
benefits of the proposed modifications.

We also analyze the homogeneous network model, a specific variant in which all the
links fail (or do not fail) with the same probability. The fact that the links share the
same failure probability results in a simplification that some methods exploit to become
more efficient. Although most networks correspond to the most general case, that is,
non-homogeneous models, in many cases they fit the homogeneous model. For this reason, and because there are interesting methods for this latter case, the homogeneous model is worth analyzing. In particular, an important modification can be made to Permutation Monte Carlo in order to make it much more efficient when applied over homogeneous models. We will introduce this variant in Section 7. Within this context we present F-Monte Carlo in Section 8. While most methods in literature are designed for non-homogeneous networks, F-Monte Carlo has been designed to operate exclusively on homogeneous models, which allows it to attain particularly high performance levels in that setting [5]. This high efficiency makes it particularly interesting as a benchmark method. In the last part of the article we present some experimental tests in which we compare our new implementation of Permutation Monte Carlo against F-Monte Carlo.

This article is an extension of a preliminary work presented at CLEI 2024 under the title "Robust Implementation of Permutation Monte Carlo For Network Reliability Estimation" [9], which discussed the numerical problems of Permutation Monte Carlo in the heterogeneous case. In the present work, we add the discussion of the homogeneous model and new experimental tests.

The work is organized as follows. In Section 2 we present the Network Static Model. Sections 3, 4 and 5 discuss the Permutation Monte Carlo method and the numerical instabilities in its usual implementations. In Section 6 we present our proposal for avoiding this stability by rearranging terms in a key computation of the method. Section 7 introduces the homogeneous network model and in Section 8 we present a brief review of F-Monte Carlo. Section 9 presents some numerical results, showing that the method is able to estimate reliability values for networks with hundreds of components. Finally, Section 10 discusses conclusions and future work.

## 2 Network Static Model

The Network Static Model is a graph $G = (V, E)$ in which a set of nodes, $V$, is connected by a set of links, $E$, with $|V| = n$ and $|E| = m$. Nodes are perfect, only links fail, and they do it independently from each other. Links are modeled by the random vector $\mathbf{X} = (X_1, \dots, X_m) \in \Omega = \{0, 1\}^m$, where $X_i = 0$ means that link $i$ is failed (we say *down*), and is equivalent to removed from the graph, whereas $X_i = 1$ means that link $i$ is operational (we say *up*). We say that vector $\mathbf{X}$ is the state of the network and we call configuration, denoted $\mathbf{x} = (x_1, \dots, x_m)$, to any possible value of $\mathbf{X}$.

The network state space $\Omega$ is the union of two disjoint subspaces, one of which associated to a network *up* condition, and the other one associated to a network *down* condition. Usually, network *up* means that a subset $K \subseteq V$ of nodes is connected by means of the *up* links, and network *down* means that, despite the *up* links, the subset $K \subseteq V$ is disconnected. These two conditions are reflected by the structure function
$\phi(\mathbf{X}) \in \{0, 1\}$ as follows: $\phi(\mathbf{X}) = 0$ means that the network is down, whereas $\phi(\mathbf{X}) = 1$ means that the network is up.

In the most general case, $\mathbb{P}\{X_i = 1\} = p_i$ and $\mathbb{P}\{X_i = 0\} = q_i = 1 - p_i$, $i = 1, \dots, m$, where the $p_i$ value may be different for each different link $i$.

We thus define the network reliability, $R$, as follows:

$$R = \mathbb{P}\{\phi(\mathbf{X}) = 1\}, \qquad (1)$$

and the unreliability, $Q$, as

$$Q = \mathbb{P}\{\phi(\mathbf{X}) = 0\} = 1 - R. \qquad (2)$$

While there are different closed formulas to express $R$ and $Q$, many of which seem simple to compute, all of them involve an exponential number of terms, so that their direct computation is out of reach, except for very small networks. As mentioned in Section 1, well known theoretical results classify exact reliability computation as an NP-hard problem (actually a \#P-hard problem, see [36] and [35]), so, unless $P = NP$, there is no hope of finding a polynomial method for performing this computation.

# 3 Dynamic Models

The difficulty in the exact computation of the reliability leads to the search for estimators, the best being those that allow a precise result with low calculation effort. Hence, in this context, one of the most important requirements for an estimator —ultimately, a random variable— is its low variance.

One of the ways to address the problem of obtaining a low variance estimator is to transform the static model into a dynamic one, provided that both models have the same reliability, and making the estimation on the latter one. Such transformation is the basis of different methods, like Splitting and Permutation Monte Carlo, originally intended to operate over dynamic models. An extensive review of these ideas can be found in [11], [22], [15], [29], [30], [1] and [17].

Permutation Monte Carlo is based on considering time $t$, along which the links are repaired (they go from down to up), provided that they are down at time $t = 0$. This is the guiding idea of many important research lines, some of which have been also adapted to the case of multiple states per component, like the work in [17]. An original and efficient approach, one of whose applications can be found in [21], is based on a distribution called D-spectrum, proposed to make a particular description of the sequences of repairs, for the case of homogeneous systems. For a complete survey of these methods, readers are advised to look over [19] and [20].

Permutation Monte Carlo consists of assigning a probability distribution to the links' repair times, and analyzing the distribution of the instant at which the subset of nodes in $K$ becomes connected.
# 4 Permutation Monte Carlo (PMC)

In PMC the Network Static Model is transformed into a dynamic model by the so called *Creation Process* (CP), proposed in [11]. At this transformation, the network state $\mathbf{X} = (X_1, \dots, X_m)$, becomes the stochastic process $\mathbf{X}(t) = (X_1(t), \dots, X_m(t)) \in \Omega = \{0, 1\}^m$, in which:

$$X_i(t) = \begin{cases} 0 & \text{if } t < \tau_i \\ 1 & \text{otherwise,} \end{cases} \quad i = 1, \dots, m. \qquad (3)$$

In words, at $t = 0$ all the links are *down*, or simply do not exist, while at times $\tau_i$ link $i$ is repaired or created (it goes *up*). Every $\tau_i$ is exponentially distributed: $\tau_i \sim \text{Exp}(\lambda_i)$, with $\lambda_i = -\ln(q_i)$, then:

$$\mathbb{P}\{X_i(1) = 1\} = \mathbb{P}\{\tau_i \le 1\} = 1 - e^{\ln(q_i)} = 1 - q_i = r_i \qquad (4)$$

$$\mathbb{P}\{X_i(1) = 0\} = \mathbb{P}\{\tau_i > 1\} = e^{\ln(q_i)} = q_i \qquad (5)$$

This means that, to observe CP at $t = 1$ is the same as observing the Network Static Model.

If at $t = 1$ all the links are *up* with probability equal to its own individual reliability, then the probability that the network is *up* at $t = 1$ is $R$. Similarly, the probability that the network is *down* at $t = 1$ is $Q$. Finally, $R = \mathbb{E}\{\phi(\mathbf{X}(1))\}$ and $Q = \mathbb{E}\{1 - \phi(\mathbf{X}(1))\}$.

There are $m!$ different ways of ordering the links' repairs, each of which is a *permutation* of the set $T = \{\tau_1, \tau_2, \dots, \tau_m\}$, where $\tau_i$ is the repairing time of link $i$, $i = 1, \dots, m$. Calling $\Pi$ to the set of all the *permutations* in $T$, the general form of some $\pi \in \Pi$ is:

$$\pi = \{\tau_{(1)}, \tau_{(2)}, \dots, \tau_{(m)}\}, \qquad (6)$$

Figure 1 shows a possible sequence of repairs. It is important to distinguish the reference to the order in the sequence (the index in brackets) and the number of link repaired, because time $\tau_{(i)}$ is the repair time of, say, the j-th link, but not necessarily $i = j$. Actually, the terms $\tau_{(i)}$, $i = 1, \dots, m$, are the order statistics of $T$. The most important value in this set is $\tau_{(c)}$, the time at which $\phi(\mathbf{X}(t))$ goes to 1, that is:

$$\phi(\mathbf{X}(t)) = \begin{cases} 0 & \text{if } t < \tau_{(c)} \\ 1 & \text{otherwise.} \end{cases} \qquad (7)$$

Thus, the event $\phi(\mathbf{X}) = 0$ in the Network Static Model is the same as event $\tau_{(c)} \ge 1$ in CP which, in the end, means that $Q = \mathbb{P}\{\tau_{(c)} \ge 1\}$. Then, given the following indicator variable:
<figure><img src="image_2.png" /><figcaption>Figure 1: The Creation Process (CP)</figcaption></figure>

$$I_P = \begin{cases} 1 & \text{if } \tau_{(c)} \ge 1 \\ 0 & \text{otherwise,} \end{cases} \qquad (8)$$

we can build a crude or standard unreliability estimation, $\hat{Q}_S$, based on the CP, as follows:

$$\hat{Q}_S = \frac{1}{N} \sum_{i=1}^{N} I_P^{(i)}, \qquad (9)$$

where $I_P^{(i)}$, $i = 1, \dots, N$, is a set of $N$ independent copies of variable $I_P$.

In order to improve the estimation in (9), it is important to find a mathematical description of the time $\tau_{(c)}$. Considering the times between repairs, as shown in Figure 1, time $\tau_{(c)}$ is composed as follows:

$$\tau_{(c)} = \sum_{i=1}^{c} \Delta_i \qquad (10)$$

where the times $\Delta_i$ can be computed, due to the properties of the exponential distribution.

Time $\Delta_1$ is the smallest (earliest) out of a set of $m$ exponentially distributed times, each one of them with rate $\lambda_i$, $i = 1, \dots, m$. So, the rate of $\Delta_1$ is $\Lambda_1 = \sum_{i=1}^{m} \lambda_i$, whereas the link repaired at time $\tau_{(1)}$ can be sampled from a discrete distribution in which the probability of each link is $\lambda_i/\Lambda_1$, $i = 1, \dots, m$. Assuming that link $z$ is sampled at time $\tau_{(1)}$, $\Delta_2$ is exponentially distributed with rate $\Lambda_2 = \Lambda_1 - \lambda_z$, and the link repaired at time $\tau_{(2)}$ follows a discrete distribution in which the probability of each link is $\lambda_i/\Lambda_2$ $i = 1, \dots, m$, $i \neq z$.

Calling $\lambda_{(i)}$ to the rate of the exponential corresponding to the link repaired at time $\tau_{(i)}$, we can present the definition of the variables $\Delta_i$ as follows:

$$\Delta_1 \sim \text{Exp}(\lambda_{(1)} + \lambda_{(2)} + \dots + \lambda_{(c)} + \dots + \lambda_{(m)}) \qquad (11)$$

$$\Delta_2 \sim \text{Exp}(\lambda_{(2)} + \dots + \lambda_{(c)} + \dots + \lambda_{(m)}) \qquad (12)$$

$$\vdots \qquad (13)$$

$$\Delta_c \sim \text{Exp}(\lambda_{(c)} + \dots + \lambda_{(m)}) \qquad (14)$$
The random variables $\Delta_1, \dots, \Delta_c$, and, therefore the variable $\tau_{(c)} = \sum_{i=1}^c \Delta_i$, are determined by the permutation $\pi \in \Pi$. Thus, according to the Total Probability Theorem, the network unreliability can be expressed as:

$$Q = \mathbb{P}\{\tau_{(c)} \ge 1\} = \sum_{\pi} \underbrace{\mathbb{P}\{\tau_{(c)} \ge 1 \mid \Pi = \pi\}}_{\gamma(\pi)} \mathbb{P}\{\Pi = \pi\} \quad (15)$$

This expression is the basis of PMC which, in the end, is an application of a more general method known as Conditional Monte Carlo. The term $\gamma(\pi)$ is the probability that the network is up at $t = 1$ when the links are repaired according to the permutation $\pi$. Thus, if for any sampled permutation, the probability $\gamma(\pi)$ can be computed, the following Monte Carlo estimation is possible:

$$\hat{Q}_P = \frac{1}{N} \sum_{i=1}^{N} \gamma(\pi^{(i)}), \quad (16)$$

where $\pi^{(i)}$, $i = 1, \dots, N$, is a set of $N$ independent replications of the permutation $\pi$. It is important to notice that, in order to compute $\gamma(\pi^{(i)})$, it is not necessary to sample the times $\Delta_1, \dots, \Delta_c$, for every replication, but only the order of repairs (permutation).

Next we show an analysis of the random variables underlying each of the two recently introduced estimators. In the case of PMC we have:

$$Q_P = \sum_{\pi} \gamma(\pi) \mathbb{P}\{\Pi = \pi\}, \quad (17)$$

from which it is simple to infer that:

$$\mathbb{V}\{Q_P\} = \sum_{\pi} \gamma(\pi)^2 \mathbb{P}\{\Pi = \pi\} - \left(\sum_{\pi} \gamma(\pi) \mathbb{P}\{\Pi = \pi\}\right)^2. \quad (18)$$

Now we address the variance of $Q_S$, the random variable underlying the standard
estimator:

$$ \mathbb{V}\{Q_S\} = Q(1 - Q) \tag{19} $$

$$ = \sum_{\pi} \gamma(\pi) \mathbb{P}\{\Pi = \pi\} \left( 1 - \sum_{\pi} \gamma(\pi) \mathbb{P}\{\Pi = \pi\} \right) \tag{20} $$

$$ = \sum_{\pi} \gamma(\pi) \mathbb{P}\{\Pi = \pi\} - \left( \sum_{\pi} \gamma(\pi) \mathbb{P}\{\Pi = \pi\} \right)^2 \pm \sum_{\pi} \gamma(\pi)^2 \mathbb{P}\{\Pi = \pi\} \tag{21} $$

$$ = \sum_{\pi} \gamma(\pi)^2 \mathbb{P}\{\Pi = \pi\} - \left( \sum_{\pi} \gamma(\pi) \mathbb{P}\{\Pi = \pi\} \right)^2 + \sum_{\pi} \gamma(\pi) \mathbb{P}\{\Pi = \pi\} - \sum_{\pi} \gamma(\pi)^2 \mathbb{P}\{\Pi = \pi\} \tag{22} $$

$$ = \mathbb{V}\{Q_P\} + \underbrace{E(\gamma(\pi)) - E(\gamma(\pi)^2)}_{\ge 0} \tag{23} $$

what proves that $\mathbb{V}\{Q_P\} \le \mathbb{V}\{Q_S\}$ and, ultimately, shows that the PMC estimation is more efficient than standard estimation.

# 5 PMC Implementation Issues

Once a permutation, $\pi$, is sampled, $\gamma(\pi)$ can be computed according to the distribution determined in [2], that in this context reads:

$$ \mathbb{P}\{\Delta_1 + \dots + \Delta_c > t \mid \pi\} = \tag{24} $$

$$ \left( \prod_{i=1}^{c} \Lambda_i \right) \sum_{j=1}^{c} \frac{e^{-\Lambda_j t}}{\Lambda_j \prod_{\substack{k=1 \\ k \neq j}}^{c} (\Lambda_k - \Lambda_j)} \tag{25} $$

For $t = 1$ this expression becomes:

$$ \mathbb{P}\{\Delta_1 + \dots + \Delta_c > 1 \mid \pi\} = \tag{26} $$

$$ \left( \prod_{i=1}^{c} \Lambda_i \right) \sum_{j=1}^{c} \frac{e^{-\Lambda_j}}{\Lambda_j \prod_{\substack{k=1 \\ k \neq j}}^{c} (\Lambda_k - \Lambda_j)} \tag{27} $$

The use of this formula is affected by both, the number of links and their reliability. As a numerical example, if $q_i = 1 \times 10^{-6}$, a typical value in high-reliability benchmark
<table><tr><td>Λ<sub>1</sub> = 170.34</td><td>S<sub>1</sub> = 1.22E-70</td></tr><tr><td>Λ<sub>2</sub> = 161.83</td><td>S<sub>2</sub> = -8.98E-66</td></tr><tr><td>Λ<sub>3</sub> = 153.31</td><td>S<sub>3</sub> = 3.08E-61</td></tr><tr><td>Λ<sub>4</sub> = 144.79</td><td>S<sub>4</sub> = -6.53E-57</td></tr><tr><td>Λ<sub>5</sub> = 136.28</td><td>S<sub>5</sub> = 9.54E-53</td></tr><tr><td>Λ<sub>6</sub> = 127.76</td><td>S<sub>6</sub> = -1.02E-48</td></tr><tr><td>Λ<sub>7</sub> = 119.24</td><td>S<sub>7</sub> = 8.17E-45</td></tr><tr><td>Λ<sub>8</sub> = 110.72</td><td>S<sub>8</sub> = -5.03E-41</td></tr><tr><td>Λ<sub>9</sub> = 102.21</td><td>S<sub>9</sub> = 2.38E-37</td></tr><tr><td>Λ<sub>10</sub> = 93.69</td><td>S<sub>10</sub> = -8.67E-34</td></tr><tr><td>Λ<sub>11</sub> = 85.17</td><td>S<sub>11</sub> = 2.38E-30</td></tr><tr><td>Λ<sub>12</sub> = 76.65</td><td>S<sub>12</sub> = -4.82E-27</td></tr><tr><td>Λ<sub>13</sub> = 68.14</td><td>S<sub>13</sub> = 6.77E-24</td></tr><tr><td>Λ<sub>14</sub> = 59.62</td><td>S<sub>14</sub> = -5.95E-21</td></tr><tr><td>Λ<sub>15</sub> = 51.10</td><td>S<sub>15</sub> = 2.48E-18</td></tr></table>

$$ \mathbb{P}\{\Delta_1 + \dots + \Delta_c > 1 | \pi\} = 1.26414\text{E-16} \qquad (29) $$

Even when the network under analysis is rather small (for larger networks, the problem is much worse), some of the terms are extremely small and their representation, and also the operations in which they are involved —by means of a programming language— can result in precision losses and numerical problems. As a reference, the single-precision floating point representation usually employed in programming languages can only represent with full precision values greater or equal to $1.18 \times 10^{-38}$ (and even allowing for precision losses, so-called sub-normal representation, only allows for numbers greater or equal to $1.4 \times 10^{-45}$). We can see that in this example there are many terms that are smaller than these values, and this then affects the Monte Carlo estimation.

An alternative would be to employ double precision floating numbers and arithmetic, which is available in most modern programming languages. Double precision representation can only represent with full precision values greater or equal to $2.23 \times 10^{-308}$. While this seems a very large range (and would indeed cover all the numbers discussed in the previous example), if we consider a similar example in which the network is composed of 90 links with $p_i = 0.9998$, $i = 1, \dots, 90$, and we still consider the replication for which $c = 15$ we find, for instance, that $S_1 = 9.42 \times 10^{-318}$.

Think now of a network model composed of 130 links with $p_i = 0.9998$, $i = 1, \dots, 90$. Here, for $c = 125$, we have $\prod_{i=1}^{c} \Lambda_i = 8.34 \times 10^{333}$ (in this case, exceeding the largest number that can be represented in double precision floating numbers, $1.80 \times 10^{308}$).

This shows that for medium-sized networks, with more than 50 links, even using double-precision arithmetic is not sufficient to correctly implement the PMC method.
# 6 Alternative Implementation for PMC

As we mentioned in the previous section, the main numerical instabilities for PMC appear when trying to compute $\gamma(\pi) = \mathbb{P}\{\Delta_1 + \dots + \Delta_c > 1 | \pi\}$. We observe that it is possible to re-arrange the formula, to ensure that neither very small nor very large terms appear:

$$ \mathbb{P}\{\Delta_1 + \dots + \Delta_c > 1 | \pi\} = \sum_{j=1}^{c} \frac{\left(\prod_{i=1}^{c} \Lambda_i\right) e^{-\Lambda_j}}{\Lambda_j \prod_{\substack{k=1 \\ k \neq j}}^{c} (\Lambda_k - \Lambda_j)} \quad (30) $$

$$ = \sum_{j=1}^{c} \frac{\Lambda_j \prod_{\substack{i=1 \\ i \neq j}}^{c} \Lambda_i}{\Lambda_j \prod_{\substack{k=1 \\ k \neq j}}^{c} (\Lambda_k - \Lambda_j)} e^{-\Lambda_j} \quad (31) $$

$$ = \sum_{j=1}^{c} \underbrace{\left( \prod_{\substack{i=1 \\ i \neq j}}^{c} \frac{\Lambda_i}{\Lambda_i - \Lambda_j} \right)}_{\alpha_j} e^{-\Lambda_j} \quad (32) $$

Even when the size of the terms $\Lambda_i$ and $\Lambda_i - \Lambda_j$ may grow excessively, it is quite easy to accept that the size of the quotients $\alpha_j$ won't grow disproportionately and they will
be perfectly manageable for a programming language:

$$
\begin{align*}
\alpha_1 &= 11628 \\
\alpha_2 &= -171360 \\
\alpha_3 &= 1175720 \\
\alpha_4 &= -4979520 \\
\alpha_5 &= 14549535 \\
\alpha_6 &= -31039008 \\
\alpha_7 &= 49884120 \\
\alpha_8 &= -61395840 \\
\alpha_9 &= 58198140 \\
\alpha_{10} &= -42325920 \\
\alpha_{11} &= 23279256 \\
\alpha_{12} &= -9405760 \\
\alpha_{13} &= 2645370 \\
\alpha_{14} &= -465120 \\
\alpha_{15} &= 38760
\end{align*}
$$

Let us now transform the terms $e^{-\Lambda_j}$ to allow proper operations with the terms $\alpha_j$:

$$
e^{-\Lambda_j} = e^{-(\lambda_{(j)} + \dots + \lambda_{(c)} + \dots + \lambda_{(m)})} \qquad (33)
$$

$$
= e^{-\lambda_{(j)}} \dots e^{-\lambda_{(c)}} \dots e^{-\lambda_{(m)}} \qquad (34)
$$

$$
= e^{-(\ln(q_{(j)}))} \dots e^{-(\ln(q_{(c)}))} \dots e^{-(\ln(q_{(m)}))} \quad (35)
$$

$$
= e^{\ln(q_{(j)})} \dots e^{\ln(q_{(c)})} \dots e^{\ln(q_{(m)})} \qquad (36)
$$

$$
= q_{(j)} \dots q_{(c)} \dots q_{(m)} \tag{37}
$$

Finally, the computation of $\mathbb{P}\{\Delta_1 + \dots + \Delta_c > 1 \mid \pi\}$ can be solved by the following
recurrence, up to $i = c - 1$.

$$
\chi_{i+1} = q_{(i+1)}(\chi_i + \alpha_{i+1}) \quad (38)
$$

starting by $\chi_1 = \alpha_1 q_{(1)}$, and setting to zero $\alpha_j$, $j = c + 1, \dots, m$. See this in a small
example in which $c = 3$ for a network composed of 7 links:

$$ \sum_{j=1}^{3} \alpha_j e^{-\Lambda_j} = \alpha_1 e^{-\Lambda_1} + \alpha_2 e^{-\Lambda_2} + \alpha_3 e^{-\Lambda_3} \quad (39) $$

$$ = \alpha_1 q_{(1)} q_{(2)} q_{(3)} q_{(4)} q_{(5)} q_{(6)} q_{(7)} \quad (40) $$

$$ + \alpha_2 q_{(2)} q_{(3)} q_{(4)} q_{(5)} q_{(6)} q_{(7)} \quad (41) $$

$$ + \alpha_3 q_{(3)} q_{(4)} q_{(5)} q_{(6)} q_{(7)} \quad (42) $$

$$ = q_{(3)} q_{(4)} q_{(5)} q_{(6)} q_{(7)} (\alpha_1 q_{(1)} q_{(2)} + \alpha_2 q_{(2)} + \alpha_3) \quad (43) $$

$$ = q_{(3)} q_{(4)} q_{(5)} q_{(6)} q_{(7)} (q_{(2)}(\alpha_1 q_{(1)} + \alpha_2) + \alpha_3) \quad (44) $$

Now, the values involved in the operations —that should be performed in the order determined by the recurrence—are the $\alpha_i$, which are not too large, and the $q_i$ which are small, but not extremely small. Thus, we can approach the target value through operations none of which involve numbers in a range that may cause a numerical problem.

# 7 PMC on Network Homogeneous Models

A particular case of the Network Static Model is that in which $P\{X_i = 1\} = p$ and $P\{X_i = 0\} = q = 1 - p, i = 1, \dots, m$. This type of model, where all the links share the same distribution, is called homogeneous model. There are algorithms designed to operate exclusively on homogeneous models, for example F-Monte Carlo [5] and D-Spectrum [17, 19]. Permutation Monte Carlo allows to significantly simplify its implementation when applied over homogeneous network models. Let us briefly recall the basics of Permutation Monte Carlo before we introduce this simplification.

In Figure 1 we see a possible sequence of repairs. In order to compute $\gamma(\pi) = P\{\tau_{(c)} \ge 1 | \Pi = \pi\}$ by the formula (27), whether we use the calculation variants proposed in Section 5 or not, we need to identify the links repaired at every $\tau_{(i)}$, that is, the permutation $\pi$.

As explained in Section 4, $\pi$ can be sampled as follows: the link repaired at time $\tau_{(1)}$ comes from a discrete distribution in which the probability of each link is $\lambda_i/\Lambda_1, i = 1, \dots, m$, where $\Lambda_1 = \sum_{i=1}^m \lambda_i$. If the link repaired at time $\tau_{(1)}$ is $z$, the link repaired at time $\tau_{(2)}$ follows a discrete distribution in which the probability of each link is $\lambda_i/\Lambda_2, i = 1, \dots, m, i \neq z$, where $\Lambda_2 = \Lambda_1 - \lambda_z$. And this continues until some link is repaired at time $\tau_{(c)}$.

Let us now apply these ideas to the an homogeneous model in which $\lambda_i = \lambda, i = 1, \dots, m$. Now, $\Lambda_1 = \sum_{i=1}^m \lambda_i = m\lambda$, therefore $\tau_{(1)}$ is sampled from a discrete distribution in which the probability of each link is $\lambda/\Lambda_1 = 1/m$. Then, no matter which link would
have been sampled first, the next link to be repaired comes from a discrete distribution in which the probability of each of the $m-1$ remaining links has probability $1/(m-1)$. This means that all the links can be successively sampled from uniform distributions, out of the set of links that have not yet been sampled.

See also that $\Lambda_i = (m-i-1)\lambda$, $i=1, \dots, c$, what shows that $\Lambda_i$ only depends on $i$, but not on which is the $i$-th repaired link. Many permutations share the same $c$, which is now the only necessary value to compute the probability $\mathbb{P}\{\tau_{(c)} \ge 1\}$. This significant simplification let us transform $\gamma(\pi)$ into $\gamma(c)$ as follows:

$$ \gamma(c) = \mathbb{P}\{\Delta_1 + \dots + \Delta_c > t \mid c\} \qquad (46) $$

Thus, once a replication is sampled, the outcome of the formula (25) only depends on the number of links repaired, but not on which those links are. There are, therefore, $m$ possible values for expression (46) —a many as possible values of $c$— which can be computed only once, at the beginning of the simulation. If we save these values into an array $P[c]$, $c = 1, \dots, m$, the mechanism of PMC can be summarized as follows: create two sets, $L_0$ containing all the links with their variables in 0 and $L_1$ empty, then (i) sample uniformly one link from $L_0$, set its variable to 1 and move it into $L_1$, (ii) compute $\phi(\mathbf{X}_1)$, where $\mathbf{X}_1$ is the set of links in $L_1$ and if $\phi(\mathbf{X}_1) = 0$ go back to (i); otherwise make $c = |L_1|$ and return $P[c]$ as the outcome of the replication.

After performing $N$ independent replications, we can make the unreliability estimation indicated in (16), replacing each $\gamma(\pi^{(i)})$ by the outcome of the algorithm introduced in the prior paragraph.

In short, when it operates over homogeneous models, the Permutation Monte Carlo algorithm can be simplified because of two reasons: firstly because in every new permutation each link is sampled from a uniform distribution, secondly —and this the main reason— because we can avoid the computation of Formula (27) in each replication.

# 8 F-Monte Carlo

The F-Monte Carlo method [5] is a simulation technique for estimating network unreliability that works exclusively over homogeneous models. It is based on an exact unreliability formula, in which some unknown, and difficult-to-calculate, coefficients are replaced by unbiased estimators.

In order to introduce the formula in which the method is based, it is necessary to partition the space $\Omega = \{0, 1\}^m$ as follows:

$$ \Omega_j = \left\{ \mathbf{x} \in \Omega : \sum_{i=1}^{m} x_i = m - j \right\} \quad j = 0, \dots, m, $$
All the configurations in $\Omega_j$ are composed of $j$ zeros and $m-j$ ones, therefore $\mathbb{P}\{\mathbf{X} = \mathbf{x} | \mathbf{x} \in \Omega_j\} = q^j p^{m-j}$. Since $|\Omega_j| = \binom{m}{j}$, we have:

$$ \sum_{j=0}^{m} \binom{m}{j} q^j p^{m-j} = 1. $$

Let $\alpha_j \le 1$, $j = 0, \dots, m$, be the fraction of configurations $\mathbf{x} \in \Omega_j$ for which $\phi(\mathbf{x}) = 0$, then:

$$ \sum_{j=0}^{m} \alpha_j \binom{m}{j} q^j p^{m-j} = Q. \qquad (47) $$

There is only one configuration $\mathbf{x}$ in $\Omega_0$, and for this configuration $\phi(\mathbf{x}) = 1$ (none of the links are down), therefore, $\alpha_0 = 0$. There is also only one configuration $\mathbf{x} \in \Omega_m$, and for this configuration $\phi(\mathbf{x}) = 0$ (all the links are down), therefore $\alpha_m = 1$. Based on this we can write:

$$ Q = \sum_{j=1}^{m} \alpha_j \binom{m}{j} q^j p^{m-j} = q^m + \sum_{j=1}^{m-1} \alpha_j \binom{m}{j} q^j p^{m-j} \qquad (48) $$

Note that only the coefficients $\alpha_j$, $j = 1, \dots, m-1$, are unknown. F-Monte Carlo is a simulation method aimed at estimating the network unreliability $Q$, based on expression (48). It consists in replacing the coefficients $\alpha_j$ by the unbiased estimators $\hat{\alpha}_j$, $j = 1, \dots, m-1$. Replacements made, the network unreliability estimator reads:

$$ \hat{Q}_{\mathrm{F}} = \sum_{j=1}^{m} \hat{\alpha}_j \binom{m}{j} q^j p^{m-j} = q^m + \sum_{j=1}^{m-1} \hat{\alpha}_j \binom{m}{j} q^j p^{m-j}. \qquad (49) $$

The idea for obtaining the estimators $\hat{\alpha}_j$ is to make Crude Monte Carlo estimations in the subspaces $\Omega_j$, $j = 1, \dots, m-1$.

Since the parameters $\alpha_j$ in (47) are just the fraction of configurations $\mathbf{x} \in \Omega_j$ for which $\phi(\mathbf{x}) = 0$, no matter the probability distribution of the links, we can make the estimations in every $\Omega_j$ using different probability values, other than $p$ and $q$. In fact, sampling under $p$ and $q$ will produce a high number of configurations $\mathbf{x}$ for which $\phi(\mathbf{x}) = 1$, and a few or none for which $\phi(\mathbf{x}) = 0$, causing an increase in the variance of the estimators $\hat{\alpha}_j$ and, therefore, in the variance of $\hat{Q}_{\mathrm{F}}$.

To make efficient estimations, we may use the random variable $\mathbf{Y} = (Y_1, \dots, Y_m) \in \Omega = \{0, 1\}^m$, where $Y_i = 0$ means that link $i$ is down, $Y_i = 1$ means that link $i$ is up
and $\mathbb{P}\{Y_i = 0\} = \mathbb{P}\{Y_i = 1\} = 0.5 \ \forall i$. Based on $\mathbf{Y}$ we may partition the state space $\Omega = \{0, 1\}^m$ as follows:

$$ \Omega_j = \{ \mathbf{y} \in \Omega : \sum_{i=1}^{m} y_i = m - j \}, \quad j = 0, \dots, m. $$

In every $\Omega_j$ we may define $\mathbf{Y}_j$, a random vector of size $m$, composed of $j$ zeros and $m-j$ ones. Every $\mathbf{Y}_j$ is uniformly distributed with $\mathbb{P}\{\mathbf{Y}_j = \mathbf{y}_j \mid \mathbf{y}_j \in \Omega_j\} = 1/(m_j!)$. As $\alpha_j$ is the fraction of configurations $\mathbf{x} \in \Omega_j$ for which $\phi(\mathbf{x}) = 0$, $\alpha_j$ is also the fraction of values of $\mathbf{y}_j$ for which $\phi(\mathbf{y}_j) = 0$.

To build the estimators $\hat{\alpha}_j$ we may generate $N_j$ samples of $\mathbf{Y}_j$. The process for obtaining these samples can be implemented as follows: (i) fill with 1s all the components of a vector $\mathbf{y}_j$, (ii) choose uniformly, and fill with 0s, $j$ of the $m$ components of $\mathbf{y}_j$. There is an algorithm called RPERM, proposed by Fishman in [12], that implements this process very efficiently. Thus, if, for every $j$, we execute this algorithm $N_j$ times, we can produce $N_j$ samples, $\mathbf{Y}_j^{(i)}$, $i = 1, \dots, N_j$ and then build the estimator of $\alpha_j$ as follows:

$$ \hat{\alpha}_j = \frac{1}{N_j} \sum_{i=1}^{N_j} 1 - \phi(\mathbf{Y}_j^{(i)}), \quad j = 1, \dots, m-1, \qquad (50) $$

Once we have built these estimators, we have all we need to compute the network reliability estimator, $\hat{Q}_F$, introduced in (49).

# 9 Computational results

## 9.1 PMC Implementation Improvements

The purpose of the experiments presented in this section is twofold, first to quantify the performance improvement of PMC against the Standard MC method, and then to show the benefit of the methodology proposed in this article.

The Standard MC unreliability estimation, $\hat{Q}_S$, is based on $N$ independent copies of the variable $\mathbf{X}$, sampled from the original distribution: $\mathbb{P}\{X_i = 1\} = p$ and $\mathbb{P}\{X_i = 0\} = q = 1 - p$, $i = 1, \dots, m$. According to this, the estimator is:

$$ \hat{Q}_S = 1 - \frac{1}{N} \sum_{i=1}^{N} \phi(\mathbf{X}^{(i)}) \qquad (51) $$

The usual measure to compare a simulation method like PMC against the standard method is the speedup value. We call it here $W_{S/P}$, and we compute it by the following
quotient:

$$W_{S/P} = \frac{t_S \times V_S}{t_P \times V_P} \qquad (52)$$

where $t_S$ and $t_P$ are the simulation times of Standard MC and PMC, and $V_S$ and $V_P$ are the respective variances. This measure allows to take into consideration both, the improvement in precision obtained by the method under consideration (in this case PMC), and the eventual difference in computational effort between the methods.

First of all we present a case much used in literature, the Dodecahedron network (see [7], [33], [32], [18] and [40]), shown in Figure 2. With the same topology, we make five different experiments, corresponding to different individual link reliability values: $r_i = 0.9; 0.99; 0.999; 0.9999; 0.99999$; corresponding to stricter reliability values and more demanding and difficult numerical estimations. Table 1 rows correspond to the individual link reliabilities; the columns show the unreliability estimation computed by PMC, and the speedups over the Standard MC method. It is possible to observe that in the low-reliability case, Standard MC is a bit more efficient than PMC; but when reliability values grow, PMC has a striking advantage, achieving a speedup of many orders of magnitude. This is also shown in graphical format in Figure 3.

<figure><img src="image_3.png" /><figcaption>Figure 2: Dodecahedron Network</figcaption></figure>

As a second test case we present the results obtained over the G-20 network. Figure 4 shows the general topology of G-n graphs. We have chosen the case where $n = 20$, corresponding to a network with 42 nodes and 82 links. Table 2 includes the following columns: individual link reliabilities, unreliability estimation computed by PMC, and speedup over the Standard MC method. Again we report five different experiments, corresponding to the link reliability values $r_i = 0.9; 0.99; 0.999; 0.9999; 0.99999$.

In third place we present results obtained over the Grid-n × n network (shown in Figure 5) for the case of $n = 5$. The Grid-5×5 network has 25 nodes and 40 links. The results are shown in Table 3, including again the same columns: individual link
<table><caption>Table 1: PMC vs. Standard MC over Dodecahedron Network</caption><thead><tr><th>r<sub>i</sub></th><th>Q&#x0302;<sub>P</sub></th><th>W<sub>S/P</sub></th></tr></thead><tbody><tr><td>0.9</td><td>2.87E-03</td><td>6.15E-01</td></tr><tr><td>0.99</td><td>2.03E-06</td><td>3.32E+01</td></tr><tr><td>0.999</td><td>1.97E-09</td><td>2.21E+04</td></tr><tr><td>0.9999</td><td>1.96E-12</td><td>2.14E+07</td></tr><tr><td>0.99999</td><td>1.96E-15</td><td>2.12E+10</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 3: Bar chart of speedup order of magnitudes for the Dodecahedron network</figcaption></figure>

<table><caption>Table 2: PMC vs. Standard MC over G-20 Network</caption><thead><tr><th>r<sub>i</sub></th><th>Q&#x0302;<sub>P</sub></th><th>W<sub>S/P</sub></th></tr></thead><tbody><tr><td>0.9</td><td>5.75E-02</td><td>1.17E-01</td></tr><tr><td>0.99</td><td>2.43E-04</td><td>5.03E-01</td></tr><tr><td>0.999</td><td>2.07E-06</td><td>2.07E+01</td></tr><tr><td>0.9999</td><td>2.03E-08</td><td>1.85E+03</td></tr><tr><td>0.99999</td><td>2.03E-10</td><td>1.82E+05</td></tr></tbody></table>

reliabilities, unreliability estimation computed by PMC, and speedup over the Standard MC method. Five different experiments are reported for each topology, corresponding to the link reliability values $r_i = 0.9; 0.99; 0.999; 0.9999; 0.99999$.
<figure><img src="image_5.png" /><figcaption>Figure 4: G-n Networks</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: PMC vs. Standard MC over Grid-n × n Network</figcaption></figure>

From the results reported in the three tables, we can see that our implementation of the PMC estimator is able to obtain good estimations for the cases of very reliable networks, which are a challenge for the Standard MC; without being hampered by the small numbers/potential underflow effect of the original PMC proposal. The results are consistently good for networks of different sizes.

Now we show, in Table 4, the last set of results from experiments designed to try
out the main proposal in this article. On the left hand side of the table, shaded in gray,
we have results that correspond to simulations carried out through the classic PMC
implementation, that is, before the proposed improvements. On the right hand side of
the table, we have simulations performed on the same network models, with a PMC
implementation in which we include the modifications proposed to avoid numerical in-
stabilities. Every experiment whose result is shown as correspond to a simulation
affected by extremely large or extremely small numbers and, as a consequence, an error
in the final unreliability estimation (the implementation stops with an exception error
message, without giving a numerical estimate). The case of extremely large numbers
comes up with the growth of the network in terms of the number of links, the extremely
small numbers correspond to the increment of the individual link reliability. The Grid-
n × n networks, for n = 8, n = 9 and n = 10 have, respectively, 112, 144 and 180
links.
<table><caption>Table 3: PMC vs. Standard MC over Grid-5×5 Network</caption><thead><tr><th>r<sub>i</sub></th><th>Q&#x0302;<sub>P</sub></th><th>W<sub>S/P</sub></th></tr></thead><tbody><tr><td>0.9</td><td>2.45E-02</td><td>2.01E-01</td></tr><tr><td>0.99</td><td>2.06E-04</td><td>1.55E+00</td></tr><tr><td>0.999</td><td>2.02E-06</td><td>9.97E+01</td></tr><tr><td>0.9999</td><td>2.01E-08</td><td>9.48E+03</td></tr><tr><td>0.99999</td><td>2.01E-10</td><td>6.95E+05</td></tr></tbody></table>

<table><caption>Table 4: $\hat{Q}_P$ over Grid-$n \times n$ Networks</caption><thead><tr><th rowspan="2">r<sub>i</sub></th><th colspan="3">Standard implementation</th><th colspan="3">Proposed implementation</th></tr><tr><th>n = 8</th><th>n = 9</th><th>n = 10</th><th>n = 8</th><th>n = 9</th><th>n = 10</th></tr></thead><tbody><tr><td>0.9</td><td>2.45E-02</td><td>2.37E-02</td><td>*****</td><td>2.45E-02</td><td>2.37E-02</td><td>2.32E-02</td></tr><tr><td>0.99</td><td>2.34E-04</td><td>*****</td><td>*****</td><td>2.34E-04</td><td>1.80E-04</td><td>1.97E-04</td></tr><tr><td>0.999</td><td>2.52E-06</td><td>*****</td><td>*****</td><td>2.52E-06</td><td>1.57E-06</td><td>2.05E-06</td></tr><tr><td>0.9999</td><td>2.55E-08</td><td>*****</td><td>*****</td><td>2.55E-08</td><td>1.55E-08</td><td>2.09E-08</td></tr><tr><td>0.99999</td><td>2.55E-10</td><td>*****</td><td>*****</td><td>2.55E-10</td><td>1.54E-10</td><td>2.09E-10</td></tr><tr><td>0.999999</td><td>*****</td><td>*****</td><td>*****</td><td>2.55E-12</td><td>1.54E-12</td><td>2.09E-12</td></tr></tbody></table>

The experiments were conducted on a DELL Inspiron 15 computer, processor Intel(R) Core(TM) i5-8250U, operating system Linux, kernel version 5.4.0-150. All the implementations were programmed in C language, compiler gcc, version 7.5.0.

The number of replications for the PMC executions in Tables 1, 2 and 3 were, respectively, $10 \times 10^6$, $1 \times 10^6$ and $4 \times 10^6$, whereas, the number of replications for the Standard MC simulations were $100 \times 10^6$, $20 \times 10^6$ and $50 \times 10^6$. These values were chosen so that all the simulations would run for approximately 100 seconds.

In the case of Table 4, in which the only purpose is to distinguish those implementations that work correctly from those that do not, a number of $10^5$ replications per PMC run was enough.

From the results reported in the four tables, we can see that our implementation of the PMC estimator is able to obtain good estimations for the cases of very reliable networks, which are a challenge for the Standard MC; without being hampered by the small numbers/potential underflow effect of the original PMC proposal. The results are consistently good for networks of different sizes, up to 180 links, which is quite demanding computationally and numerically.

While the Standard MC has decreasing precision levels for more reliable graphs, we can see that PMC has a better performance, leading to increasing gains when com-
pared against Standard MC. The algorithmic ideas behind PMC are stable even for very high link reliability values, and our implementation also removes numerical unstabilities, leading to these very good results.

## 9.2 PMC and F-Monte Carlo over Homogeneous Models

The experiments shown in this section focus on the PMC variant proposed for the homogeneous case called here H-PMC. Again we base the experimental comparison on the speedup, now defined as follows:

$$ \mathbb{W}_{m/H} = \frac{t_m \times \mathbb{V}_m}{t_H \times \mathbb{V}_H} \qquad (53) $$

where *H* refers to H-PMC and *m* to the method to which H-PMC is compared.

We run the simulations on the benchmark networks shown in Figures (6), (7), (8), and (9) which are composed as follows:

<table><thead><tr><td>Network</td><td>Nodes</td><td>Links</td></tr></thead><tbody><tr><td>Ring-1×16</td><td>16</td><td>16</td></tr><tr><td>Ring-2×16</td><td>32</td><td>48</td></tr><tr><td>Ring-3×16</td><td>48</td><td>80</td></tr><tr><td>Ring-4×16</td><td>64</td><td>112</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 6: Ring-1×16 Network</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Ring-2×16 Network</figcaption></figure>

In Tables 5 and 6 we evaluate the efficiency increase of the proposed implementation for homogeneous systems. To do so, we estimate the unreliability twice, firstly using
<figure><img src="image_9.png" /><figcaption>Figure 8: Ring-3×16 Network</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Ring-4×16 Network</figcaption></figure>

the standard PMC implementation with the improvements proposed in Section 6 ($\hat{Q}_P$), secondly using the implementation proposed in Section 7 for homogeneous systems ($\hat{Q}_H$). Since the proposed modification does not modify the accuracy, the improvement in speedup is only due to a reduction in the execution times. Specifically, the execution times of the improved variant are in the order of two thirds of the times of the standard variant. And this improvement grows, slightly, along with the growth of the reliability and the network size.

<table><caption>Table 5: PMC vs. H-PMC over Ring-2×16</caption><thead><tr><th>r<sub>i</sub></th><th>Q&#x0302;<sub>P</sub></th><th>Q&#x0302;<sub>H</sub></th><th>W<sub>P/H</sub></th></tr></thead><tbody><tr><td>0.9</td><td>1.04E-02</td><td>1.04E-02</td><td>1.46</td></tr><tr><td>0.99</td><td>2.65E-06</td><td>2.63E-06</td><td>1.52</td></tr><tr><td>0.999</td><td>2.02E-09</td><td>1.98E-09</td><td>1.53</td></tr><tr><td>0.9999</td><td>1.95E-12</td><td>1.92E-12</td><td>1.53</td></tr></tbody></table>

<table><caption>Table 6: PMC vs. H-PMC over networks with r<sub>i</sub> = 0.999</caption><thead><tr><th>Network</th><th>Q&#x0302;<sub>P</sub></th><th>Q&#x0302;<sub>H</sub></th><th>W<sub>P/H</sub></th></tr></thead><tbody><tr><td>Ring-1×16</td><td>6.35E-05</td><td>6.36E-05</td><td>1.51</td></tr><tr><td>Ring-2×16</td><td>2.02E-09</td><td>1.98E-09</td><td>1.53</td></tr><tr><td>Ring-3×16</td><td>1.81E-09</td><td>1.69E-09</td><td>1.58</td></tr><tr><td>Ring-4×16</td><td>2.11E-09</td><td>1.81E-09</td><td>1.76</td></tr></tbody></table>
In Tables 7 and 8 we show an experimental comparison between the variant of PMC suitable for homogeneous systems and F-Monte Carlo. The unreliability is estimated by both methods, $\hat{Q}_H$ and $\hat{Q}_F$, and the corresponding speedup $W_{F/H}$ is evaluated. An improvement is observed which, as in the two previous experiments, grows along with the growth of the network reliability and the network size.

<table><caption>Table 7: F-MC vs. H-PMC over Ring-3×16</caption><thead><tr><th>r<sub>i</sub></th><th>Q&#x0302;<sub>F</sub></th><th>Q&#x0302;<sub>H</sub></th><th>W<sub>F/H</sub></th></tr></thead><tbody><tr><td>0.9</td><td>1.05E-02</td><td>1.04E-02</td><td>1.17</td></tr><tr><td>0.99</td><td>2.67E-06</td><td>2.63E-06</td><td>3.27</td></tr><tr><td>0.999</td><td>1.96E-09</td><td>1.98E-09</td><td>3.95</td></tr><tr><td>0.9999</td><td>1.88E-12</td><td>1.92E-12</td><td>4.01</td></tr></tbody></table>

<table><caption>Table 8: F-MC vs. H-PMC over networks with r<sub>i</sub> = 0.999</caption><thead><tr><th>Network</th><th>Q&#x0302;<sub>F</sub></th><th>Q&#x0302;<sub>H</sub></th><th>W<sub>F/H</sub></th></tr></thead><tbody><tr><td>Ring-1×16</td><td>6.35E-05</td><td>6.36E-05</td><td>2.01</td></tr><tr><td>Ring-1×16</td><td>1.96E-09</td><td>1.98E-09</td><td>4.04</td></tr><tr><td>Ring-1×16</td><td>1.91E-09</td><td>1.69E-09</td><td>6.22</td></tr><tr><td>Ring-1×16</td><td>2.00E-09</td><td>1.81E-09</td><td>7.46</td></tr></tbody></table>

# 10 Conclusions and future work

In this work, we discuss the PMC simulation method for estimating network reliability values, and we show that the original formulation involves the computation of both extremely small and extremely large values, which for large networks exceed usual floating point precision capacities. We present an alternative way to apply the method, by rewriting terms of a key formula, avoiding this issue and resulting in a more robust implementation. Computational experiments show that the improved PMC implementation proposal is able to estimate the reliability in medium and large sized networks where the standard implementation did not work, as numerical exceptions resulted in aborted executions. It is interesting to note that these problems arise even in cases where the link reliabilities and the overall network reliability are not extremely high, i.e, the problem is with the intermediate steps of the computations, and not with the expected final results. The numerical results, presented over a variety of network topologies of
different sizes and structures, show the efficiency and robustness of the PMC method when implemented following our proposal.

We also show how to take advantage of some issues in the implementation of PMC that allow it to be adapted very efficiently to the homogeneous case, in which all links share the same probability distribution. We shown, experimentally, that the variant adjusted to this case is particularly efficient, outperforming, when tested on homogeneous networks, both the standard PMC and F-Monte Carlo, the latter being a method specifically designed for homogeneous networks.

While both implementations improve both the robustness and the efficiency of the PMC method, we do not have a theoretical result which allows to know in advance the boundaries in network size and link reliability values which guarantee robust computations. Future work can include studying both this issue, and also apply similar analysis to other simulation methods from literature. This aspect of the numerical implementation of the methods is not much studied in previous literature, and for practical applications it can be of relevance.

## References

[1] Michael Amrein and Hans R. Künsch. A variant of importance splitting for rare event estimation: Fixed number of successes. *ACM Trans. Model. Comput. Simul.*, 21(2):13:1–13:20, February 2011.

[2] Márton Balázs. Sum of independent exponentials. Online notes - https://people.maths.bris.ac.uk/mb13434/sumexp.pdf, 2005.

[3] Z. I. Botev, P. L'Ecuyer, R. Simard, and B. Tuffin. Static network reliability estimation under the Marshall-Olkin copula. *ACM Transactions on Modeling and Computer Simulation*, 26(2):Article 14, 2016.

[4] Eduardo Canale, Franco Robledo, Pablo Romero, and Pablo Sartor. Monte Carlo methods in diameter-constrained reliability. *Optical Switching and Networking*, 14, Part 2(0):134–148, 2014. Special Issue on RNDM 2013.

[5] Eduardo Canale, Franco Robledo, Pablo Romero, and Pablo Sartor. Monte Carlo methods in diameter-constrained reliability. *Optical Switching and Networking*, 14:134–148, 2014. Special Issue on RNDM 2013.

[6] H. Cancela and M. El Khadiri. A recursive variance-reduction algorithm for estimating communication-network reliability. *IEEE Transactions on Reliability*, 44(4):595–602, December 1995.
[7] H. Cancela and M. El Khadiri. The recursive variance-reduction simulation algorithm for network reliability evaluation. *IEEE Transactions on Reliability*, 52(2):207-212, June 2003.

[8] H. Cancela, M. El Khadiri, and G. Rubino. Rare events analysis by Monte Carlo techniques in static models. In G. Rubino and B. Tuffin, editors, *Rare Event Simulation using Monte Carlo Methods*, chapter 7, pages 145-170. John Wiley & Sons, 2009.

[9] Hector Cancela, Leslie Murray, and Gerardo Rubino. Robust implementation of Permutation Monte Carlo for network reliability estimation. In *Proceedings of the L Latin American Computer Conference (CLEI 2024)*, Bahia Blanca, Argentina, pages 1-8, 2024.

[10] M.C. Easton and C.K. Wong. Sequential destruction method for Monte Carlo evaluation of system reliability. *IEEE Transactions on Reliability*, R-29(1), April 1980.

[11] T. Elperin, I. B. Gertsbakh, and M. Lomonosov. Estimation of network reliability using graph evolution models. *IEEE Transactions on Reliability*, 40(5):572-581, Dec 1991.

[12] George S. Fishman. *Monte Carlo: concepts, algorithms, and applications*. Springer series in operations research. Springer-Verlag, New York, 1996.

[13] G.S. Fishman. A comparison of four Monte-Carlo methods for estimating the probability of s-t connectedness. *IEEE Transactions on Reliability*, R-35(2):145-155, June 1986.

[14] G.S. Fishman. A Monte Carlo sampling plan for estimating network reliability. *Operational Research*, 34, 1986.

[15] M. J. J. Garvels. *The Splitting Method in Rare Event Simulation*. PhD thesis, Faculty of mathematical Science, University of Twente, The Netherlands, 2000.

[16] Ilya Gertsbakh, Eyal Neuman, and Radislav Vaisman. Monte Carlo for estimating exponential convolution. *Communications in Statistics - Simulation and Computation*, 44(10):2696-2704, 2015.

[17] Ilya Gertsbakh, Reuven Rubinstein, Yoseph Shpungin, and Radislav Vaisman. Permutational methods for performance analysis of stochastic flow networks. *Probability in the Engineering and Informational Sciences*, 28(1):21-38, 2014.
[18] Ilya Gertsbakh, Yoseph Shpungin, and Radislav Vaisman. *Ternary Networks: Reliability and Monte Carlo*. Springer Publishing Company, Incorporated, 2014.

[19] Ilya B. Gertsbakh and Yoseph Shpungin. *Models of Network Reliability: Analysis, Combinatorics, and Monte Carlo*. CRC Press, Inc., USA, 1st edition, 2009.

[20] Ilya B. Gertsbakh, Yoseph Shpungin, and Radislav Vaisman. *Ternary Networks - Reliability and Monte Carlo*. Springer Briefs in Electrical and Computer Engineering. Springer, 2014.

[21] Ilya B. Gertsbakh, Yoseph Shpungin, and Radislav Vaisman. D-spectrum and reliability of a binary system with ternary components. *Probability in the Engineering and Informational Sciences*, 30(1):25–39, 2016.

[22] Paul Glasserman, Philip Heidelberger, Perwez Shahabuddin, and Tim Zajic. Splitting for rare event simulation: Analysis of simple cases. In *Proceedings of the 28th conference on Winter simulation*, WSC 1996, Coronado, CA, USA, December 8-11, 1996, pages 302–308, 1996.

[23] K-P Hui, N. Bean, M. Kraetzl, and D. Kroese. The tree cut and merge algorithm for estimation of network reliability. *Probability in the Engineering and Informational Sciences*, 17(1):23–45, 2003.

[24] K-P Hui, N. Bean, M. Kraetzl, and D. Kroese. The cross-entropy method for network reliability estimation. *Annals of Operations Research*, 134:101–118(18), February 2005.

[25] Kazuhiro Kobayashi, Hozumi Morohosi, and Tatsuo Oyama. Applying path-counting methods for measuring the robustness of the network-structured system. *International Transactions in Operational Research*, 16(3):371–389, 2009.

[26] H. Kumamoto, K. Tanaka, and K. Inoue. Efficient evaluation of system reliability by Monte Carlo method. *IEEE Transactions on Reliability*, R-26(5):311–315, Dec. 1977.

[27] H. Kumamoto, K. Tanaka, K. Inoue, and E.J. Henley. Dagger-sampling Monte Carlo for system unavailability evaluation. *IEEE Transactions on Reliability*, R-29(2):122–125, June 1980.

[28] P. L'Ecuyer, G. Rubino, S. Saggadi, and B. Tuffin. Approximate Zero-Variance Importance Sampling for Statistic Network Reliability Estimation. *IEEE Transactions on Reliability - TR*, 60(3):590–604, 2011.
[29] Pierre L'Ecuyer, Valérie Demers, and Bruno Tuffin. Rare events, splitting, and Quasi-Monte Carlo. *ACM Trans. Model. Comput. Simul.*, 17(2), April 2007.

[30] Pierre L'Ecuyer, François Le Gland, Pascal Lezaud, and Bruno Tuffin. Splitting techniques. In Gerardo Rubino and Bruno Tuffin, editors, *Rare Event Simulation using Monte Carlo Methods*, chapter 3, pages 39–61. John Wiley & Sons, 2009.

[31] M. Lomonosov. On Monte Carlo estimates in network reliability. *Probability in the Engineering and Informational Sciences*, 8:245–264, 1994.

[32] Aleksandr Moshnikov. Evaluation of network reliability and element importance metrics. In *Majorov International Conference on Software Engineering and Computer Systems*, 2020.

[33] L. Murray, H. Cancela, and G. Rubino. A splitting algorithm for network reliability. *IIE Transactions*, 45(2):177–189, 2013.

[34] H. Pérez-Rosés. Sixty years of network reliability. *Mathematics in Computer Science*, 12(3):275–293, 2018.

[35] Scott J. Provan and Michael O. Ball. The Complexity of Counting Cuts and of Computing the Probability that a Graph is Connected. *SIAM Journal on Computing*, 12(4):777–788, 1983.

[36] A. Rosenthal. Computing the reliability of complex networks. *SIAM Journal on Applied Mathematics*, 32(2):384–393, 1977.

[37] S. M. Ross. A new simulation estimator of system reliability. *Journal of Applied Mathematics and Stochastic Analysis*, 7(3), 1994.

[38] G. Rubino. Network reliability evaluation. In K. Bagchi and J. Walrand, editors, *State-of-the art in performance modeling and simulation*, chapter 11, pages 275–302. Gordon and Breach Books, 1998.

[39] L. Schäfer, S. García, and V. Srithammavanh. Simplification of inclusion-exclusion on intersections of unions with application to network systems reliability. *Reliability Engineering and System Safety*, 173:23–33, 2018.

[40] Stefano Sebastio, Kishor S. Trivedi, Dazhi Wang, and Xiaoyan Yin. Fast computation of bounds for two-terminal network reliability. *European Journal of Operational Research*, 238(3):810–823, 2014.