Héctor Cancela
email: cancela@fing.edu.uy

Leslie Murray
email: lmurray@fceia.unr.edu.ar

Gerardo Rubino
email: gerardo.rubino@inria.fr

Reliability Estimation for Stochastic Flow Networks with Dependent Arcs

Keywords: reliability, Creation Process, Destruction Process, Splitting, dependent arcs Static Network Model SFNM Stochastic Flow Network Model SMC Standard Monte Carlo

The Creation and the Destruction Processes (CP and DP) are the basis of many efficient Monte Carlo methods for estimating the unreliability of highly reliable networks on both, the static and the stochastic flow network models, for the case of independent components. Some of these methods are based on the Splitting variance reduction techniques. Due to the Splitting basic mechanism, they operate over CP. Here we propose a Splittingbased Monte Carlo method, using the Marshall-Olkin copula for the case of non independent components. This proposal operates on the DP, because the Marshall-Olkin copula model over networks is quite related -and somehow similar-to DP, which is unusual but here necessary. Before addressing the proposal, the article presents a review of the methods based on CP and DP. At the end, a comparative experimental analysis shows the efficiency of the proposed approach.

NOTATION

G

graph representing a network topology V set of n nodes for any type of network model E set of m links for SNM A set of m arcs for SNFM C set of capacities of the m arcs for SNFM X random vector (X 1 , X 2 , . . . , X m ) indicating the state of links (for SNM) or arcs (for SFNM) x a realization of the random vector X φ(X) structure function for SNM and SFNM F (X) total amount of flow, from the set of nodes that produce flow to the set of nodes that consume flow probability that a trajectory started at threshold i reaches threshold i+1 earlier than t = 1, in ML-DP β i estimator of β i I. BACKGROUND AND PRELIMINARIES In reliability settings, a network is usually modeled as an undirected graph G = (V, E), in which V is a set of n nodes and E, a set of m links connecting these nodes. The links and/or the nodes can fail, and, according to this, the network can be either working or failed. The reliability (unreliability) is the probability that the network is working (failed).

D

A. Static Network Model

In the Static Network Model (SNM), the nodes in V are perfect (they never fail), while the links in E can be in one of two possible states, operational (also called up) or failed (also called down). In the operational state the links provide a perfect connection, whereas in the failed state they perform as if they were removed from the graph.

There is a subset of nodes K ⊂ V , called terminal nodes, such that the network is working if, despite links' failures there is a connected subgraph of G containing all the terminal nodes.

In SNM the links are modeled by the random vector X = (X 1 , . . . , X m ) ∈ {0, 1} m , where X i = 0 means that link i is failed, and X i = 1 means that link i is operational. We say that vector X is the state of the network and we call configuration, denoted x = (x 1 , . . . , x m ), any possible value of X.

According to the configuration, the network assumes one of two possible states, operational or failed. This information is formalized in the so-called structure function, φ, which maps the 2 m possible configurations into {0, 1} as follows: φ(x) = 1 if, when removing from G the links which are failed in x, all the nodes in K belong to the same connected component of the resulting subgraph, and φ(x) = 0 otherwise.

Given the probabilities p i = P{X i = 1} and q i = P{X i = 0} = 1 -p i , i = 1, . . . , m, the unreliability of SNM, ζ, is defined as:

ζ = P{φ(X) = 0} = 1 -E{φ(X)}. (1)
Many algorithms have been designed to compute the exact value of these metrics. However, as their computation is NP-hard [START_REF] Provan | The complexity of counting cuts and of computing the probability that a graph is connected[END_REF], they only apply to small or, at least, medium sized networks, while the case of large networks can only be tackled by simulation.

Standard or crude Monte Carlo (SMC) simulation is the simplest and most straightforward method for estimating network reliability. However, if networks are highly reliable, the efficiency of SMC estimation falls and simulations must be enhanced, for instance, by means of variance reduction methods.

The SMC estimator, ζ, is based on the following indicator variable:

I = 1 if φ(X) = 0, 0 otherwise. (2)
Then, given N independent samples of I, I (i) , i = 1, . . . , N , the SMC unreliability estimator is:

ζ = 1 N N i=1 I (i) . (3)
We immediately have E{ ζ} = ζ and V{ ζ} = ζ(1 -ζ)/N . The ratio V{ ζ} 1/2 /E{ ζ} is interpreted as a statistical bound of the relative error in the estimation. It follows that, as ζ goes to zero, this relative error goes to infinity. In other words, as said before, the SMC unreliability estimator looses precision as the network becomes more reliable (E{ ζ} tends to zero faster than V{ ζ} 1/2 ). It is necessary, therefore, to improve the estimation by means of variance reduction methods [START_REF] Cancela | The recursive variance-reduction simulation algorithm for network reliability evaluation[END_REF], [START_REF] Cancela | Rare events analysis by Monte Carlo techniques in static models[END_REF], [START_REF] Easton | Sequential destruction method for Monte Carlo evaluation of system reliability[END_REF], [START_REF] Fishman | A Monte Carlo Sampling Plan for Estimating Network Reliability[END_REF], [START_REF] Hui | The cross-entropy method for network reliability estimation[END_REF], [START_REF] Lomonosov | On Monte Carlo estimates in network reliability[END_REF], [START_REF] Cancela | On the reliability estimation of stochastic binary systems[END_REF].

B. Stochastic Flow Network

The Stochastic Flow Network Model (SFNM) is based on a directed graph G = (V, A, C), where V is a set of n nodes, A is a set of m arcs connecting these nodes and C the set of the arcs' capacities. The capacity is the maximum amount of flow that can pass through an arc. Nodes are perfect (they never fail) whereas capacities are modeled as a random vector X = (X 1 , . . . , X m ) in which every X i is a multi-valued positive random variable.

Recall that a flow on G is a positive function on the arcs, satisfying Kirchoff's law: for each node u having strictly positives indegree and outdegree, the sum of the flow on the arcs entering u is equal to the sum of the flows on the arcs leaving u. In SFNM, we can more generally define two disjoint subsets of nodes, R ⊂ V and S ⊂ V , such that nodes in R "produce flow", whereas nodes in S "consume flow". If we denote a flow from node x to node y, by f (x, y), then it is easy to prove that:

u∈R v ∈R f (u, v) = v∈S u ∈S f (u, v) =: f tot . (4)
Here, f tot is the total flow circulating in the network, generated in R and consumed in S. Now a configuration of the capacities is x = (x 1 , . . . , x m ), where x i is a capacity value for arc i. A flow f is said to be compatible with the network capacities, written here f ∈ C(x), iff for all arc , we have f ( ) ≤ x . We call F (x) the maximum amount of flow compatible with C(x) that can be sent from R to S, that is,

F (x) = argmax f ∈C(x) f tot .
Thus, F (X) is a random variable. In this context, we are typically given a threshold D > 0 and we want to calculate P{F (X) < D}. The structure function here is φ(X) = 1 when F (X) ≥ D, and φ(X) = 0 otherwise. Therefore, like in SNM, the SFNM unreliability, ζ, is defined as:

ζ = P{φ(X) = 0} = 1 -E{φ(X)}. (5)
The SFNM applies to a wide variety of domains, from data networks, where flows are made up of information packets, to water pipe or electricity distribution networks, to name a few examples. In all these cases reliability is a fundamental parameter, both in performance analysis of existing systems as well as an a target parameter in the design of new ones. Now, the same drawbacks affecting SNM also affect SFNM. This means that, if the network is large, the exact reliability computation is impossible. Estimation by simulation is an appropriate solution to find an approximate value, being SMC the easiest and most straightforward simulation method at hand. However, if the network is highly reliable, SMC loses efficiency and simulations must be improved by means of variance reduction methods.

II. INTRODUCTION

In this section we introduce the two central aspects of the paper, the variance reduction technique that we propose to deal with rarity, and the way we take into account the dependencies between components.

A. Splitting Monte Carlo estimation

One of the most powerful families of techniques to deal with the highly reliable situations is called Splitting, and our paper develops a variant of this class of methods. Splitting is defined on dynamical systems, and the most used models in network reliability, such as our considered SFNM, are static. Our initial step will be transforming our model into a dynamic one, by means of the so-called Creation Process (CP), or the Destruction Process (DP). Both of these stochastic processes have been introduced in [START_REF] Elperin | Estimation of network reliability using graph evolution models[END_REF]. They transform the static SNM into a dynamic process with failures and repairs that occur along an artificial time. In the case of the DP, the guiding idea is as follows: at time t = 0 the network components are operational, and they begin to fail randomly in times proportional to the failure probability of each one. Along with one of the individual failures, the network fails at a time which is, ultimately, proportional to its own failure probability. The CP is similar, but the components and the network start failed, and the events of interest are their "repairs" rather than their "failures".

The Multi-Level CP (ML-CP), proposed in [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF], is an adaptation of the CP to SFNM, for the case of independent components (arcs). However, it does not directly apply to models with dependent components.

Neither CP nor ML-CP, are variance reduction methods. However, they are the basis of many of them, like [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF], [START_REF] Murray | A Splitting algorithm for network reliability estimation[END_REF] and [START_REF] Elperin | Estimation of network reliability using graph evolution models[END_REF], but always on models in which links are independent. In the cases of [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF] and [START_REF] Murray | A Splitting algorithm for network reliability estimation[END_REF], the methods are basically applications of Splitting [START_REF] Glasserman | Splitting for rare event simulation: Analysis of simple cases[END_REF], [START_REF] L'ecuyer | Rare events, splitting and Quasi-Monte Carlo[END_REF], [START_REF] L'ecuyer | Splitting techniques[END_REF], [START_REF] Villén-Altamirano | Analysis of restart simulation: Theoretical basis and sensitivity study[END_REF], [START_REF] Borodina | On splitting scenarios for effective simulation of reliability models[END_REF], [START_REF] Vaisman | Splitting sequential monte carlo for efficient unreliability estimation of highly reliable networks[END_REF], [START_REF] Vaisman | Sequential stratified splitting for efficient monte carlo integration[END_REF] over the corresponding dynamic model (either CP or ML-CP).

B. Dependencies between components

As stated before, SFNM are very useful and used models in different areas. They always work from the knowledge of the structure of a multicomponent system and the behavior of each individual component characterized by their probability of being operational. This last information is enough if we assume independence between components, and this was done until very recently, in almost every dependability analysis. It also means that we try to define as components, subsystems of the given system that can be reasonably assumed to work or fail independently of each other. Now, in real life, components' independence is rare (as an example, see [START_REF] Gonzalez | Analysis of dependencies between failures in the uninett ip backbone networks[END_REF]). Moreover, in many cases, it is difficult to define components this way, and often, the statistical available data comes for more "naturally" defined components that exhibit clear dependencies. Take as an example the case of natural catastrophes (floods, earthquakes, tsunamis, storms), where when components are geographically close they show a high correlation regarding their states (operational or failed). See, for instance, the analysis in our previous works [START_REF] Barrera | Topological optimization of reliable networks under dependent failures[END_REF] or [START_REF] Matus | On the Marshall-Olkin copula model for network reliability under dependent failures[END_REF]. The cost of ignoring dependencies can be high, as shown in previous citations; see also [START_REF] Lapcevic | Impact of dependent failures on the availability of the optical network[END_REF].

The most widespread approach to address dependencies in network models is the Marshall-Olkin copula [START_REF] Marshall | A multivariate exponential distribution[END_REF], which takes the form of a dynamic process in which failures occur sequentially at random times, just like in DP.

The Splitting-like methods published so far only apply to CP. In fact, Splitting perfectly fits CP because, if we formalize the evolution of the system as trajectories in the state space, the rare event is a region of that space that is "hard to reach" and, therefore, the trajectories that get close to it need to be stimulated or "pushed" towards the rare event. If we attempt to apply Splitting over DP, exactly as it is applied over the CP, the event of interest occurs if the trajectories "do not reach" some particular region of the state space. In other words, it will be necessary to prevent DP trajectories from progress beyond some particular bound, which is exactly the opposite of what Splitting does. A Splitting-like method over DP requires, therefore, a change in the framework, which is solved in this article applying Splitting over DP in a way that is significantly different as it it is applied over CP.

C. Contributions and structure of the paper

The main contribution of this paper is a method based on an application of Splitting over the Multi-Level DP (ML-DP) for SFNM, which operates on models with independent and dependent arcs as well.

Sections III and IV are a review of CP and DP, while Sections V and VI look over the ML-CP and the ML-DP. The case of dependent arcs for the ML-DP is approached in Section VII. Some experimental tests are introduced in Section VIII while conclusions and future work are presented in Section IX. The proposed method is introduced in Section VI and adapted to models with dependent arcs in Section VII.

III. THE CREATION PROCESS

The Creation Process (CP) [START_REF] Elperin | Estimation of network reliability using graph evolution models[END_REF] is the underlying mechanism of many variance reduction methods. The CP transforms a SNM into a dynamic model, making the random vector X become a stochastic process X(t) = (X 1 (t), . . . , X m (t)), in which:

X i (t) = 0 if t < τ i , 1 if t ≥ τ i , i = 1, . . . , m. (6)
In short, link i is repaired (switched from failed to operational) at random time τ i . The variables τ 1 , . . . , τ m are independent of each other. The key of CP is that, if τ i is exponentially distributed with rate λ i = -ln(q i ), then P{X i (1) = 1} = p i and P{X i (1) = 0} = q i . Thus, observing the CP at t = 1 is equivalent to observing the original SNM. Therefore, ζ = P{φ(X(1)) = 0}.

Recall that the order statistics of the set {τ 1 , . . . , τ m } is the vector S CP = (τ (1) , . . . , τ (m) ), where τ (k) is the kth-smallest value in the set. In other words, S CP is the set of all the repairing times, arranged in ascending order.

In S CP , there is one particular time, τ (c) , such that:

φ(X(t)) = 0 if t < τ (c) , 1 if t ≥ τ (c) . (7)
Time τ (c) is the repair time of the so called critical link. Then, ∀t, t ≥ τ (c) , there is a connected subgraph of G that includes the subset K of nodes and all the links that at time τ (c) are already repaired. This time evinces the way CP and SNM are related, since the event {τ (c) > 1} in the CP, is the same as {φ(X) = 0} in the corresponding SNM.

Although it is not what it was intended for, the CP can readily be the basis of a SMC estimation. In order to implement this estimation, it is necessary to sample N independent copies of the order statistics

S (i) CP = (τ (i) (1) , . . . , τ (i) (m)
), i = 1, . . . , N , and to compute the following indicator variable, for each one of them:

I (i) CP = 1 if τ (i) (c) > 1 0 otherwise. ( 8
)
Based on this indicator, the SMC unreliability estimator of ζ, subject to CP, is:

ζ CP = 1 N N i=1 I (i) CP . (9)
Again, see that the unreliability estimations introduced so far are based on the following equivalent events: {τ (c) > 1} in this CP-based SMC estimation, and {φ(X) = 0} in SMC over SNM.

Permutation Monte Carlo (PMC) [START_REF] Elperin | Estimation of network reliability using graph evolution models[END_REF] is the first variance reduction method thought of as to fit CP. It is a variant of Conditional Monte Carlo in which the probability of the event {τ (c) > 1} is conditioned on the ordering in which links are repaired (permutation). Calling Ω the set of all possible permutations, then:

ζ = P{τ (c) > 1} = ω∈Ω P{τ (c) > 1 | Ω = ω} γ CP (ω) P{Ω = ω}. ( 10
)
Sampling N independent permutations ω (i) , i = 1, . . . , N , and computing (exactly) N times the function γ CP (ω (i) ) = P{τ (c) > 1 | Ω = ω (i) }, the PMC unreliability estimator, subject to CP, is:

ζ P CP = 1 N N i=1 γ CP (ω (i) ). (11)
A well known result, proved in [START_REF] Elperin | Estimation of network reliability using graph evolution models[END_REF], shows that V{ ζ P CP } ≤ V{ ζ CP }, what makes this PMC estimator more efficient than the corresponding SMC one in [START_REF] Cancela | Rare events analysis by Monte Carlo techniques in static models[END_REF].

Whether CP is the basis of an SMC or a PMC estimator, it is important to see the ordering sequence of repairing times as a trajectory like the one shown in Figure 1. The simplest and most direct way to build such sequence is: (i) sampling all the repairing times and then (ii) sorting them in ascending order. However, based on the properties of the exponential distribution, it is possible to build these sequences much more efficiently, as it will be explained next.

As the target event is {τ (c) > 1}, only repairing times up to τ (c) are of interest. Times between repairs, ∆ i , i = 1, . . . , c, defined as follows, are also exponentially distributed with rates Λ i , i = 1, . . . , c. Calling

T = m i=1 λ i , then Λ 1 = T , Λ 2 = T -λ x
, where link x is the one repaired at time τ (1) , Λ 3 = T -λ x -λ y , where link y is the one repaired at time τ (2) , and so on. On the other hand, the probability that the link repaired at time τ (1) is link x is given by λ x /T , the probability that the link repaired at time τ (2) is link y is given by λ y /(T -λ x ), etc. As a consequence, time τ (c) = c i=1 ∆ i is the sum ofc exponentially distributed random times with rates Λ i , i = 1, . . . , c, and its probability density function is, therefore:

f ∆ 1 +∆ 2 +•••+∆c (t) = c i=1 Λ i c j=1 e -Λ j t c k=1 k =j (Λ k -Λ j ) . ( 12
)
In both CP-based estimations considered so far, i.e. the one in ( 9) and the one in [START_REF] Murray | A Splitting algorithm for network reliability estimation[END_REF], the target value is the probability of the event {τ (c) > 1}. The computation in (9) saves and accumulates '1's for every sampled sequence in which the event {τ (c) > 1} occurs, and '0's otherwise. At the end, all the accumulated values are divided by the number of samples. On the other hand, the mechanism in [START_REF] Murray | A Splitting algorithm for network reliability estimation[END_REF] saves and accumulates the probability of the event {τ (c) > 1} for every sampled sequence -computed by means of ( 12)-and finally all the accumulated values are divided by the number of samples.

Another CP-based variance reduction approach is the one introduced in [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF]. It is a novel application of Splitting [START_REF] Glasserman | Splitting for rare event simulation: Analysis of simple cases[END_REF], [START_REF] L'ecuyer | Rare events, splitting and Quasi-Monte Carlo[END_REF], [START_REF] L'ecuyer | Splitting techniques[END_REF], [START_REF] Villén-Altamirano | Analysis of restart simulation: Theoretical basis and sensitivity study[END_REF], [START_REF] Borodina | On splitting scenarios for effective simulation of reliability models[END_REF], [START_REF] Vaisman | Splitting sequential monte carlo for efficient unreliability estimation of highly reliable networks[END_REF], [START_REF] Vaisman | Sequential stratified splitting for efficient monte carlo integration[END_REF] in which the trajectories of repairing times start at t = 0 and progress towards t = 1. The method uses a partition of the interval [0, 1] formed by a set of N T times, namely 0 < 1 < . . . < N T < 1. These time bounds perform as the thresholds of the splitting mechanism. Thus, whenever a repair time exceeds one of the thresholds, the corresponding trajectory actually crosses this threshold. Then, if the critical link, τ (c) , is not yet part of such trajectory, many (more than one) new trajectories are started at the crossing point in the expectation that at least one of them reaches the next threshold before τ (c) is sampled (see [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF]). Computing the number of trajectories that, at the end of the process have crossed all the thresholds, including the final bound, t = 1, it is possible to obtain an unbiased unreliability estimator, ζ S CP , based on Splitting subject to CP [START_REF] Amrein | A variant of importance splitting for rare event estimation: Fixed number of successes[END_REF], as follows:

ζ S CP = N T i=0 p i , (13)
where p i = n i+1 /n i , i = 1, . . . , N T , being n i the number of crosses through threshold i (n 0 is the number of trajectories started at t = 0 and n N T +1 is the number of trajectories that reach t = 1). The variance of ζ S CP is much smaller than the variance of the crude estimator in [START_REF] Elperin | Estimation of network reliability using graph evolution models[END_REF].

IV. THE DESTRUCTION PROCESS The Destruction Process (DP) [START_REF] Elperin | Estimation of network reliability using graph evolution models[END_REF] is another random process that models networks temporal evolution and -as the CP-it is the basis of variance reduction methods. In DP, links fail (switch from operational to failed) at times τ 1 , . . ., τ m , that is:

X i (t) = 1 if t < τ i , 0 if t ≥ τ i , i = 1, . . . , m. (14)
If τ i is exponentially distributed with rate λ i = -ln(p i ), then P{X i (1) = 1} = p i and P{X i (1) = 0} = q i , meaning that -just like in CP-observing the DP at t = 1 is equivalent than observing the SNM. The unreliability, ζ, is also P{φ(X(1)) = 0}. Now, the order statistics of {τ 1 , . . . , τ m } is the ordered set of failing times S DP = {τ (1) , . . . , τ (m) }. Time τ (c) in S DP is the failing time of the so called critical link. It is the instant since which the remaining operational links -if any-are not enough to guarantee the connection among the nodes in the subset K. See that event {τ (c) < 1} in DP, is the same as event {φ(X) = 0} in SNM.

DP can also be the basis of an SMC estimation. In order to implement this estimation, it is necessary to sample N independent copies of the order statistics, S (i)

DP = {τ (i) (1) , . . . , τ (i) (m) }, i = 1, .
. . , N , and to compute the following indicator variable, for each one of them:

I (i) DP = 1 if τ (i) (c) < 1, 0 otherwise. ( 15
)
Finally, the SMC unreliability estimator, subject to DP, is:

ζ DP = 1 N N i=1 I (i) DP . (16)
PMC also fits DP, in a similar way as it does to CP. Being Ω the set of all possible permutations (accepting now that permutations come from the ordering sequence of failing times), then:

ζ = P{τ (c) < 1} = ω P{τ (c) < 1 | Ω = ω} γ DP (ω) P{Ω = ω}. (17)
Sampling N independent permutations, ω (i) , i = 1, . . . , N , and computing (exactly) N times the function γ DP (ω

(i) ) = P{τ (c) < 1 | Ω = ω (i)
}, the PMC unreliability estimation subject to DP, is:

ζ P DP = 1 N N i=1 γ DP (ω (i) ) (18)
As shown in [START_REF] Elperin | Estimation of network reliability using graph evolution models[END_REF],

V{ ζ P DP } ≤ V{ ζ DP }
, what makes this PMC estimation more efficient than the corresponding SMC estimation in (3).

V. THE MULTI-LEVEL CP CP-based methods are focused on the determination of the probability (extremely low for a highly reliable network) that τ (c) > 1, whereas DP-based methods are focused on the determination of the probability (extremely low for a highly reliable network) that τ (c) < 1. These probabilities are, in both cases, the network unreliability.

Both CP and DP have been designed over networks in which links are modeled as two-valued variables. In [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF], the CP has been adapted -and named Multi-Level CP-to Stochastic Flow Network Models (SFNM) in which the variables that model the arcs can take more than two values. The rest of this section contains a complete review of the ML-CP.

As said in Section II, in SFNM the capacities of the m arcs are represented by a random vector X = (X 1 , . . . , X m ), in which every X i is a multi-valued random variable. Without loss of generality, let us consider that X i , i = 1, . . . , m, are identical discrete random variables distributed as follows:

X i =            M 1 w.p. p 1 , M 2 w.p. p 2 ,
. . .

M n w.p. p n , 0 w.p. p 0 , M 1 > M 2 > . . . > M n > 0, i = 1, . . . , m. (19)
Here, X i = M 1 means that arc i is fully operational, and X i = 0 means that arc i is completely failed. Intermediate values indicate that arc i is partially operational (or partially failed, according to convenience).

When we have X i = M 1 ∀i, F (X) assumes its maximum possible value which, if the network has been properly designed, should be greater than the demand, D. If, due to failures, the capacity of one or more arcs -possibly all of them-is different (lower) than M 1 , then F (X) may be less than D.

Given the distribution of X, the SFNM reliability is P{φ(X) = 1} = E{φ(X)}, and the unreliability, P{φ(X) = 0} = 1 -E{φ(X)}. Many algorithms have been designed to compute the exact value of these terms (see, for example, [START_REF] Soh | An efficient cutset approach for evaluating communication-network reliability with heterogeneous link-capacities[END_REF] and [START_REF] Yeh | An improved sum-of-disjoint-products technique for symbolic multi-state flow network reliability[END_REF]). However, as their computation is NP-hard [START_REF] Provan | The complexity of counting cuts and of computing the probability that a graph is connected[END_REF], they only apply to small or, at least, medium sized networks, whereas the case of large networks can only be estimated by simulation. SMC is the simplest and most straightforward method for obtaining an estimate. However, if networks are highly reliable, the efficiency of SMC estimations fall and simulations must be enhanced by variance reduction methods, like the ones proposed in [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF], [START_REF] Botev | Reliability of stochastic flow networks with continuous link capacities[END_REF], [START_REF] Bulteau | A Monte Carlo simulation of the flow network reliability using importance and stratified sampling[END_REF], [START_REF] Cancela | Highly reliable stochastic flow network reliability estimation[END_REF], [START_REF] Fishman | Evaluating reliability of stochastic flow networks[END_REF], [START_REF] Gertsbakh | Permutational methods for performance analysis of stochastic flow networks[END_REF].

In order to build a variance reduction method, the ML-CP [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF] transforms SFNM into a dynamic model, making the random vector X become the stochastic process X(t) = (X 1 (t), . . . , X m (t)), in which:

X i (t) = 0 if t < τ i M j if t ≥ τ i i = 1, . . . , m and j ∈ {1, . . . , n}. (20)
A particular type of "repair" occurs at time τ i , as the value of X i (t) switches from 0, which means completely failed, to one of the partially operational states, namely, to one of the possible M j > 0. The mechanism proposed for deciding which M j is assigned to X i (t) at time τ i , is explained next.

Let us accept, by now, that there exists a set {λ, t 0 = 0, t 1 , t 2 , . . . , t n-1 , t n = 1} such that λ is the rate of the exponentially distributed random time τ i , and t 1 , t 2 , . . . , t n-1 , a set of fixed times that creates a partition on [0, 1] to let τ i fall into one of the corresponding intervals, with the probability indicated in Figure 2. According to this, if t j-1 < τ i ≤ t j , which happens w.p. p j , j = 1, . . . , n, then X i (t) switches from 0 to M j , and remains in M j forever. If τ i > 1, which happens w.p. p 0 , X i (t) remains in 0 forever. As a result, X i (1) = M j w.p. p j , j = 1, . . . , n, and X i (1) = 0 w.p. p 0 , the same distribution as the one proposed in [START_REF] Gonzalez | Analysis of dependencies between failures in the uninett ip backbone networks[END_REF]. In [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF] it is proved that the parameters that meet all these requirements are:

     λ = -ln(p 0 ), t j = ln(1 -p 1 -p 2 -. . . -p j ) ln(p 0 ) j = 1, . . . , n -1. (21)
Thus, as the arcs are "repaired" as explained, the flow V (X(t)) becomes a constant piece-wise increasing function of time, three replications of which are shown in Figure 3.

Every jump in V (X(t)) corresponds to a repair, whereas not necessarily every arc repair has a jump associated with. Now, the attention must be focused on the evolution of the replications of V (X(t))seen as trajectories-to check whether they are below D at t = 1, or not. This can be the basis of an SMC unreliability estimation, because V (X(1)) < D is an event that indicates the network failure, just as {τ (c) > 1} in CP, or {φ(X) = 0} in SNM.

The value of V (X(t)) can be easily determined at any time, by means of any of the well known algorithms like Ford-Fulkerson [START_REF] Ford | Flows in Networks[END_REF], Edmonds-Karp [START_REF] Edmonds | Theoretical improvements in algorithmic efficiency for network flow problems[END_REF] or Goldberg-Rao [START_REF] Goldberg | Beyond the flow decomposition barrier[END_REF].

If N independent copies of the trajectories V (X(t)) (i) , t ≥ 0, i = 1, . . . , N , are sampled, it is possible to compute the following indicator variable, for each one of them: Finally, the SMC estimation of ζ, subject to the ML-CP, is:

I (i) ML-CP = 1 if V (X(1)) (i) < D, 0 otherwise. ( 22
)
ζ ML-CP = 1 N N i=1 I (i) ML-CP . (23)
As shown in Figure 4, it is not necessary to simulate all the trajectories up to t = 1 because, as V (X(t)) is an increasing function of time, whenever a trajectory reaches D at some time t < 1, it is guaranteed that such trajectory will not be below D at t = 1.

As in the unreliability estimation, the event of interest is the network failure, and considering that here the network failure is modeled by trajectories that at t = 1 are still below D (i.e. trajectories associated with points like (3)), we say that these trajectories -for which the event of interest actually occurs-are successful. Similarly, we call unsuccessful those trajectories for which the event of interest does not occur (for example, those associated with points (1) and ( 2)), even though they model the operational state of the network.

Thus, the SMC unreliability estimation, subject to the ML-CP, ζ ML-CP , is the ratio between the number of successful trajectories and the total number of trajectories initiated at t = 0.

t d t=0 t=1 V Vm (1) (2) (3)
Fig. 4. Successful and unsuccessful trajectories in ML-CP For highly reliable networks, most of the trajectories -almost all of them-will be unsuccessful. A possible interpretation of this fact is as follows: as highly reliable networks are "repaired very soon", trajectories rise up immediately, being extremely rare to find them below D at t = 1. The unreliability estimation in [START_REF] Marshall | A multivariate exponential distribution[END_REF] needs, therefore, to be improved by a variance reduction technique. Otherwise, when the true value of the unreliability is very small, but not zero, the outcome of the estimations will always be zero.

The solution proposed in [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF] consists of the use of Splitting as shown in Figure 5. As the points of interest are the ones associated with trajectories crossing the bound t = 1 below D (the red ones), the candidates are those trajectories that, at any t < 1, are still underneath D. So, trajectories for which V (X(t)) < D at times like 1 , 2 , . . . , are cloned, attempting to finally increase the number of red points. This application of Splitting over the ML-CP [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF], with thresholds 1 , 2 , . . . , achieves an important variance reduction compared to the crude estimation in [START_REF] Marshall | A multivariate exponential distribution[END_REF]. In a setting with N T thresholds (out of which, 0 = 0 and N T +1 = 1 are implicit), if α i is the probability that a trajectory started at threshold i reaches threshold i+1 underneath D, i = 0, . . . , N T , then:

ζ = N T i=0 α i . (24)
An estimator per stage, α i , can be computed as the ratio between the number of trajectories reaching i below D, and the number of trajectories started from i-1 . Once the simulation is complete, and one estimator per stage is computed, the Splitting network unreliability estimation, subject to the ML-CP, is:

ζ S ML-CP = N T i=0 α i . (25)
Some experiments in [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF] show that the variance reduction achieved by Splitting over the ML-CP, produce estimations whose relative errors are less than 1% for unreliabilities in the order of 2.00 × 10 -12 . The duration of the simulations is, approximately, 1.5 million times less than the time that would take to perform SMC simulations, over the same network model, to achieve the same variance and, consequently, the same relative error.

VI. THE MULTI-LEVEL DP

The method to be introduced in this section, named Multi-Level DP, is the core of this article and also one of its main contributions.

Without any loss of generality, let the arcs' capacities be modeled by the following set of discrete random variables:

X i =            M n w.p. p n , M n-1 w.p. p n-1 , . . . M 1 w.p. p 1 , M 0 w.p. p 0 , M n > . . . > M 0 = 0, i = 1, . . . , m. (26)
The ML-DP transforms the random vector X into the stochastic process X(t) = (X 1 (t), . . . , X m (t)), in which:

X i (t) = M n if t < τ i , M j if t ≥ τ i , i = 1, . . . , m and
0 ≤ j ≤ n -1. (27)
A particular type of "failure" occurs at time τ i , as the value of X i (t) switches from M n , which means completely operational, to one of the partially failed states, namely, to one of the possible M j < M n . The mechanism proposed for deciding which M j is assigned to X i (t) at time τ i , is explained next.

Let us accept, by now, that there exists a set {λ, t 0 = 0, t 1 , t 2 , . . . , t n-1 , t n = 1} such that λ is the rate of the exponentially distributed random time τ i , and t 1 , t 2 , . . . , t n-1 , a set of fixed times that creates a partition on [0, 1] to let τ i fall into one of the corresponding intervals, with the probability indicated in Figure 6. According to this, if t j < τ i ≤ t j+1 , which happens w.p. p j , X i (t) switches from M n to M j , and remains in M j forever. If τ i > 1, which happens w.p. p n , X i (t) remains in M n forever. As a result, X i (1) = M j w.p. p j , j = 0, . . . , n, the same distribution as the one proposed in [START_REF] Yeh | An improved sum-of-disjoint-products technique for symbolic multi-state flow network reliability[END_REF]. The parameters that meet all these requirements are:

     λ = -ln(p n ), t j = ln(1 -p 0 -p 1 -. . . -p j-1 ) ln(p n ) , j = 1, . . . , n -1. (28)
Thus, as the arcs "fail" as explained, the flow V (X(t)) becomes a constant piece-wise decreasing function of time, three replications of which are shown in Figure 7.

Every jump in V (X(t)) corresponds to some failure, whereas not necessarily every arc failure has a jump associated with. Now, the attention must be focused on the evolution of the replications of V (X(t)), to see whether they are below D at t = 1, or not. Here, V (X(1)) < D is an event that indicates the network failure, just as {τ (c) < 1} in DP, {φ(X) = 0} in SNM or {F (X) < D} in SFNM. Then, based on the ML-DP model, a SMC unreliability estimation can be obtained as follows.

Given N independent copies of the trajectories V (X(t)) (i) , t ≥ 0, i = 1, . . . , N , it is possible to compute the following indicator variable, for each one of them:

I (i) ML-DP = 1 if V (X(1)) (i) < D 0 otherwise. (29)
Finally, the SMC estimation of ζ, subject to the ML-DP, is:

ζ ML-DP = 1 N N i=1 I (i) ML-DP . (30)
As shown in Figure 8, it is not necessary to simulate all the trajectories until t = 1 because, as V (X(t)) is a decreasing function of time, whenever V (X(t)) < D for some t < 1, it is guaranteed that such trajectory will be below D at t = 1.

Here, trajectories are classified exactly as in section V. According to this, the ones associated with points like (3) are successful, whereas the ones associated with points like (1) and (2) are unsuccessful.

Thus, the SMC unreliability estimation, ζ ML-DP , is the ratio between the number of successful trajectories and the total number of trajectories initiated at t = 0.

(3) Fig. 8. Successful and unsuccessful trajectories in ML-DP For highly reliable networks, most of the trajectories -almost all of them-will be unsuccessful. A possible interpretation of this fact is as follows: as highly reliable networks "almost never fail", trajectories are most of the time close to V max , being extremely rare to find them below D, at t = 1.

Thus, if the networks under analysis are highly reliable, the estimation in [START_REF] Fishman | Evaluating reliability of stochastic flow networks[END_REF] needs to be improved by a variance reduction technique, otherwise, when the true value of the unreliability is very small, but not zero, the outcome of the estimations will always be zero.

The solution proposed in this article consists of an application of Splitting over the ML-DP, as shown in Figure 9.

Keep in mind that for highly reliable networks, the arcs are most of the time at their maximum capacity, M n , or at least at capacity values that are close to M n . Then, if all -or almost all-the arcs are most of the time close to their maximum capacity, trajectories tend to be all the time close to V max , and they rarely fall down.

As we are interested in the probability that the trajectories reach the bound D earlier than t = 1, it is necessary to stimulate trajectories that tend to deviate from V max , "pushing" them down to D. Splitting can do this, and it does it by cloning recursively the trajectories that crosses the thresholds located between V max and D. Formally speaking, the candidates are those trajectories for which V (X(t)) < V max , at times t < 1. Then, trajectories down-crossing the thresholds 1 , 2 , . . ., earlier than t = 1 are cloned, finally increasing the number of red points. In a general setting we have N T thresholds (out of which, 0 = V max and N T +1 = D are implicit). Calling β i the probability that a trajectory started at threshold i reaches threshold i+1 earlier than t = 1, and accepting that ζ is the probability that a trajectory that starts from V max at t = 0, reaches D earlier than t = 1, then:

ζ = N T i=0 β i . ( 31
)
As we are interested in an estimator for ζ,we need a set of estimators for β i , i = 0, . . . , N T . The estimators β i can be computed as the ratio between the number of trajectories reaching i earlier than t = 1, and the number of trajectories started from i-1 .

Finally, the Splitting network unreliability estimation, subject to the ML-DP, is:

ζ S ML-DP = N T i=0 β i . (32)
In the example shown in Figure 9 six trajectories are launched at t = 0, all of them with capacity V (X(0)) = V max . They evolve and two of them reach 1 earlier that t = 1, therefore, β 0 = 2/6 = 0.33. As a consequence, four new trajectories are launched from 1 , two from every crossing point (we consider a splitting factor of 2). These four trajectories evolve and two of them reach 2 producing β 1 = 2/4 = 0.50. Finally, β 2 = 3/4 = 0.75, what makes ζ S ML-DP = 0.33 × 0.50 × 0.75 = 0.12.

VII. SFNM WITH NON-INDEPENDENT ARC FAILURES

In all the models used so far in this article, a failure affecting a single component (either a link in SNM or an arc in SFNM) is independent from failures affecting the other components. This hypothesis, though an important mathematical simplification, is often quite far from real networks behavior.

Systems with dependencies have been studied in recent times, with a particular emphasis in the last few years. Works like [START_REF] Barrera | Topological optimization of reliable networks under dependent failures[END_REF], [START_REF] Matus | On the Marshall-Olkin copula model for network reliability under dependent failures[END_REF], [START_REF] Botev | Static network reliability estimation under the Marshall-Olkin copula[END_REF] and [START_REF] Lin | Quantifying the impact of correlated failures on stochastic flow network reliability[END_REF] address models with dependent failures and show that the dependency has a high impact in the reliability value.

Think for example in the very simple SFNM shown in Figure 10. According to the capacities, X 1 and X 2 , there is a maximum amount of flow generated in node s that can reach node t. This maximum is simply X 1 + X 2 . Suppose that the service requirement is such that D = 0, meaning that any amount of flow from s to t above zero is enough and, therefore, the only failure condition corresponds to the event X 1 + X 2 = 0. The unreliability of this network is then, ζ = P{X 1 + X 2 = 0} or, considering that

s t X1 X2 Fig. 10. Simple Stochastic Flow Network X 1 , X 2 ≥ 0: ζ = P{X 1 = 0 ∧ X 2 = 0} (33)
Then, when both arcs are independent, the unreliability, called now ζ i , is:

ζ i = 1 = 0} × P{X 2 = 0} (34)
This setting, in which the behavior of every component is independent from those of the rest, is the background idea of all the models considered so far in this article.

Suppose now that, due to any reason, the behavior of both arcs is positive correlated in the following sense: if the event X 2 = 0 has actually occurred, the probability of event X 1 = 0 is increased with respect to the corresponding marginal, that is P{X 1 = 0 | X 2 = 0} > P{X 1 = 0}. Suppose also that failures in arc 1 influence arc 2 in the same way, that is P{X 2 = 0 | X 1 = 0} > P{X 2 = 0}. Then, the expressions above can be easily adapted to show that now the network unreliability, called ζ d , and ζ i are related as follows:

ζ d = P{X 1 = 0 ∧ X 2 = 0} > P{X 1 = 0} × P{X 2 = 0} = ζ i . ( 35
)
The opposite (ζ d < ζ i ) occurs if the behavior of the arcs is negative correlated. Thus, to get closer to a realistic value of the unreliability (either as and estimate or an exact computation) it is important to find the way to consider dependencies among the network components. Think now of an SFNM, in which the arcs capacity distribution is not necessarily the same for every arc. In the ML-DP introduced earlier in this article, the arcs fail independently from each other at times determined by the exponentially distributed random variables τ i , i = 1, . . . , m. No matter which M j < M n is the capacity after τ i , the failure actually occurs at time τ i , whose probability density function is λ i e -λ i x over [0, ∞). Thus:

P{τ i > x} = e -λ i x , i = 1, . . . , m. (36)
This distribution satisfies the lack of memory property so, given the times x, y ≥ 0,

P{τ i > x + y | X > y} = P{τ i > x}, i = 1, . . . , m, (37)
indicating that the residual lifetime of each arc is independent of its age. For single continuous random variables, the only probability distribution that satisfies the lack of memory property is the exponential distribution. The lack of memory property can be extended to the random vector (τ 1 , . . . , τ m ), over [0, ∞) m if, given the times x 1 , . . . , x m , and y ≥ 0:

P{τ 1 > x 1 + y, . . . , τ m > x m + y | τ 1 > y, . . . , τ m > y} = P{τ 1 > x 1 , . . . , τ m > x m }. ( 38
)
This implies that the residual lifetime of the arcs is independent of the age of the whole set of arcs. If m = 1, as stated above, τ 1 is exponentially distributed. For the case of m = 2, Marshall and Olkin have shown in [START_REF] Marshall | A multivariate exponential distribution[END_REF] that, if τ 1 and τ 2 are both exponentially distributed with parameters λ 1 , λ 2 > 0, there exists a third parameter λ 1,2 ≥ 0, such that the following bivariate distribution:

P{τ 1 > x 1 , τ 2 > x 2 } = e (-λ 1 x 1 -λ 2 x 2 -λ 1,2 max {x 1 ,x 2 }) ( 39
)
satisfies the lack of memory property.

The distributions referred to in ( 36) and ( 39) verify the following two properties: they are exponentials, so that every τ i is useful to model the lifetime of an arc, and they satisfy the lack of memory property. It is worth noting that, in the case of the bivariate distribution, the parameter λ 1,2 is responsible for a dependence between τ 1 and τ 2 . As said in [START_REF] Mai | Sampling exchangeable and hierarchical Marshall-Olkin distributions[END_REF], expression [START_REF] Rubino | Rare event simulation using Monte Carlo methods[END_REF] makes the marginal distributions of τ 1 and τ 2 be coupled correctly in order to let the dependence structure guarantee the lack of memory property. See also that, if λ 1,2 = 0, the dependency is extinguished and τ 1 and τ 2 become independent.

Marshall and Olkin also proved in [START_REF] Marshall | A multivariate exponential distribution[END_REF] that for m > 2, it is possible to find a distribution with the lack of memory property in which the marginals are exponentials, if there is a distribution with these two properties for all the (m -1)-subvectors of (τ 1 , . . . , τ m ). As these distributions exist for m = 1 and m = 2, it follows inductively that they exist ∀ m. The general case of this distribution (m ≥ 1), called Marshall-Olkin distribution, is going to be introduced next.

Let I ⊂ V = {1, . . . , m} be a non empty subset of the set of arcs, and let C be the collection of all the possible sets I, except for the empty set. A random vector (τ 1 , . . . , τ m ) has the Marshall-Olkin distribution if:

P{τ 1 > x 1 , . . . , τ m > x m } = e - I∈C λ I max i∈I {x i } , (40)
where λ I ≥ 0. As there are 2 m different sets in V , but there is no need to consider the empty set,

| C |= 2 m -1. If, for example, m = 3, V = {1, 2, 3} and C = {{1}, {2}, {3}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}}, then (τ 1 , τ 2 , τ 3 ) has the Marshall-Olkin distribution if P{τ 1 > x 1 , τ 2 > x 2 , τ 3 > x 3 } = e ψ , where ψ = -λ 1 x 1 -λ 2 x 2 -λ 3 x 3 -λ 1,2 max{x 1 , x 2 } -λ 1,3 max {x 1 , x 3 } -λ 2,3 max {x 2 , x 3 } -λ 1,2,3 max {x 1 , x 2 , x 3 }.
When the arcs are independent, their expected lifetimes are 1/λ i , i = 1, . . . , m, and this is all what a simulation needs for sampling the failing times. When the hypothesis of independence does not hold, and one attempts to model the network in terms of the Marshall-Olkin distribution, some of the parameters -precisely those who measure the degree of dependency-may not be easy to determine. To overcome this difficulty, the model can be expressed in terms of a different set of variables, which is more explicitly related to the behavior of the components. These new variables, called here E I , I ∈ C, are independent and exponentially distributed with parameter λ I . These variables are the components of the so-called Marshall-Olkin copula and they are related to the vector (τ 1 , . . . , τ m ), as follows:

τ k = min{E I | I ∈ C, k ∈ I}, k = 1, . . . , m (41)
where

I:k∈I λ I > 0 ∀ k = 1, . . . , m. (42)
Expression (42) shows that every arc must belong to, at least, one I with λ I > 0. Now, we can analyze the model in terms -only-of the set of variables E I , with the following interpretation. All the arcs in I start to operate at time t = 0, and they fail simultaneously, as the consequence of a shock occurring at time E I . Any arc of the network fails when the first shock including it occurs, just as pointed out in (41).

It remains to be proved that the failing times determined by the variables E I , I ∈ C, in the Marshall-Olkin copula model, are such that the vector (τ 1 , . . . , τ m ) follows the Marshall-Olkin distribution.

For simplicity, let us first analyze the case of networks composed of only two arcs, for which the life times are τ 1 and τ 2 . Here, the surviving condition is:

For arc 1:

τ 1 > x 1 , For arc 2: τ 2 > x 2 , → (τ 1 > x 1 ∧ τ 2 > x 2 ),
In terms of the Marshall-Olkin copula, this condition becomes:

For arc 1:

E 1 > x 1 ∧ E 1,2 > x 1 , For arc 2: E 2 > x 2 ∧ E 1,2 > x 2 , → (E 1 > x 1 ∧ E 2 > x 2 ∧ E 1,2 > max{x 1 , x 2 }). (43)
Let us now go back to the general case of a network composed of m arcs. Then, expression (43) takes the form:

(E I > max i∈I {x i }, ∀I ∈ C)
Finally:

P{τ 1 > x 1 , . . . , τ m > x m } = P{E I > max i∈I {x i }, ∀I ∈ C} (44)
= I∈C e -λ I max i∈I {x i } (45) = e - I∈C λ I max i∈I {x i } . (46)
This proves that, when the failing times are determined by the variables E I , I ∈ C, vector (τ 1 , . . . , τ m ) follows the Marshall-Olkin distribution. As a result, in order to build a network model under the Marshall-Olkin distribution, one has to focus on the failing times, and proceed as follows:

1) Select the arcs for which one knows the mean time of failures that depend only on the arc and that are not related with the rest of the network. With this information, set properly the parameters λ I of the variables E I for which I contains only the selected arc. 2) Select sets of arcs for which one knows the mean time of the simultaneous failure of the set (shock), for failures that depend only on the selected set and that are not related with the rest of the network. With this information, set the parameters λ I of the variables E I for which I contains the elements of the selected set. Every I ∈ C not considered in any of these two items is implicitly associated with variables E I with λ I = 0.

VIII. EXPERIMENTAL TESTS

This section shows the results of some experiments in which the SFNM unreliability is estimated by means of Splitting over the ML-DP, for the case of non-independent arc failures.

Let us first point out how Splitting over the ML-DP fits the model introduced in Section VII. The time to failure in sets composed of only one arc are modeled by exponentially distributed random variables E I , I ∈ {{1}, {2}, . . . , {m}}, and it is reasonable to consider all the sets, because each one of these variables Fig. 12. Lattice 5 × 5 Network models the individual behavior of every arc, i.e. the behavior without consideration of the dependencies, which is always the first description at hand in the analysis of systems' components. Besides, the sets composed of more than one arc are modeled in terms of variables E I , where I is a reference to a set of arcs. For example, if I = {1, 3, 12}, the variable E I models the simultaneous failure of arcs 1, 3 and 12. It is usual, therefore, to consider only a few number of sets, out of all the possible ones, because to observe the joint behavior of a number of arcs, that may not even be close to each other, is not an easy task or, at least, is not as easy as observing only one arc.

All the experiments are based on the networks shown in Figures 11 and12. These benchmark topologies are some of those usually employed in the literature; and we experiment with a wide range of arc failure probability values, as this parameter has a larger effect in the performance of the methods than the topology itself. In both of them all the sets composed of one arc have the same probability distribution, thus, we will refer to them as I s (the s indicates single). The sets composed of more than one arc also share the same probability distribution; we will refer to them as I m (the m indicates multiple). In other words, I s is a reference to any of the sets composed of one arc, whereas I m a reference to any of the sets composed of more than one arc. In the case of the Fishman Network, the sets I m are used to model failures in nodes, because they consider the simultaneous failure of the arcs that touches the same node (for both, incoming and outgoing arcs). There is one set I m per node, except for s = 0 and t = 9. In the case of the Lattice 5 × 5 Network there are four sets I m clearly indicated (shaded) in Figure 12.

The dynamics of the sets I s is simple and clear: the arcs start at their maximum capacity value M n and, at a time determined by the variable E Is they fail, what in terms of the mechanism introduced in Section VI, means that their capacities drop to some M j < M n . The dynamics of the sets I m is quite similar, with the only difference of a simultaneous change in the capacity value of all the arcs that belong to the set. Thus, at a time determined by the variable E Im , the capacity of all the arcs in the set drop, but not necessarily all of them to the same value. The ML-DP is flexible enough to permit the definition of simultaneous failures in which the arcs drop to capacity values that are not necessarily equal. See also that the capacity value after a failure of an arc as a single set, does not necessarily need to be the same as the capacity after the failure of the same arc, as part of a set. And the same holds for one arc as part of different sets.

For both of the networks used in the experiments, we propose the capacity distributions shown in Table I, Type A and Type B. These sets of capacities were selected with the purpose of testing Splitting over the ML-DP over two different network models. See that, for the sets I s , in the Type B the capacity values are a bit more scattered than in Type A, thus, the variance for Type B is larger that the variance for Type A. For the sets I m , the distribution, and therefore the variance, is the same for both types. Assuming that, in real models, a failure affecting any arc as part of a multiple set should be rarer than a failure affecting the same arc itself, we proposed that the probability that the capacity of a single arc drops to zero is ten times larger than the probability that the capacity of all the arcs of a multiple set drop to zero.

All the probabilities in Table I are based on p 0 -the probability that any arc is completely failed, without consideration of dependencies-. This value is used as a parameter to identify the experiments, as it is a reference to the different instances of the network models. According to the explanations given in Sections VI and VII, the rates λ Is of the variables E Is are -ln(p 8 ) for Type A and -ln(p 10 ) for Type B, whereas the rates λ Im of the variables E Im are -ln(p 8 ) for Type A and -ln(p 10 ) for Type B. Tables II and III show some results in which the efficiency of Splitting over the ML-DP is evaluated in terms of the Speedup, which is defined as follows:

Speedup = t ML-DP × V ML-DP t S ML-DP × V S ML-DP , (47)
where the numerator refers to the execution time and the variance of SMC estimation of ζ, subject to the ML-DP (see [START_REF] Fishman | Evaluating reliability of stochastic flow networks[END_REF]), and the denominator refers to the same parameters for the method under test, i.e. Splitting over the ML-DP. Assuming that the most efficient method is the one referred to in the denominator, the Speedup is an indication of the improvement in the execution time to achieve the same variance as the worst method, or, the variance reduction achieved by the best method when the execution times of both methods are equal. The splitting implementation in all the experiments is Fixed-Splitting, with a splitting factor of 2. The number N of trajectories launched was selected in order to obtain a relative error (RE) of, approximately, 3%. As expected, the Speedup increases as the network becomes more reliable. The method under test seems to be more efficient for networks of Type B, in which the variance of the capacity distribution is higher. In the second set of experiments the objective was to observe the influence of the dependency on the unreliability estimations. To do this, the distribution of the sets I s was fixed to the one corresponding to p 0 = 0.00001. Table IV reports the results. The values of the probability p 0 were proposed, and the corresponding values of p 8 , for Type A, and p 10 , for Type B, were set accordingly. For extremely low values of p 0 , the failures of the sets I m have no influence and the unreliability estimated is the same as the corresponding network model in which the arcs are all independent from each other. But when the values of p 0 become larger -and extremely larger-the unreliability is mostly determined by the shocks affecting sets of multiple arcs. See that when the difference in the probabilities is large in favor of the sets I m , making that the arcs almost never fail individually, the behavior of Lattice 5 × 5 Network is only determined by the four shaded sets of arcs, each of which fails with probability p 0 . In this case it is easy to see that the unreliability becomes p 0 4 + 4p 0 3 (1 -p 0 ) + 2p 0 2 (1 -p 0 ) 2 . For p 0 = 0.001, p 0 = 0.01 and p 0 = 0.1 this last expression returns, respectively, the values 2.00 × 10 -06 , 2.00 × 10 -04 and 1.99 × 10 -02 which are almost identical to the values reported in the last three lines of Table IV. It remains to explain the way trajectories are cloned at the cross of a threshold. Let us suppose that, as shown in Figure 13 (a), arc a fails at time t k making V (X(t)) drop from v 1 to v 2 . If the splitting factor is 2, two new trajectories should start at point (1), both of them with the state of the original trajectory at this point. The state of the original trajectory is composed of two elements: (i) the set of arcs that have already failed at time t k -and (ii) the arc a failing and moving V (X(t)) from v 1 to v 2 . We must consider, then, that two trajectories that verify (i) and (ii) must emerge from point (1), as shown in Figure 13 (b).

In order to continue these two new trajectories, for each one of them, a new arc will be sampled at time t > t k and, eventually V (X(t)) will drop down to some v < v 2 , as shown in Figure 13 (c). In summary, at the crossing in point (1), restart two new trajectories at t k , both with V (X(t)) = v 2 and accept that, at point (1), there is one incoming and two outgoing trajectories. It is possible that, as shown in Figure 14 (a), due to the failure of some arc a, a trajectory crosses more than one threshold simultaneously. The mechanism is essentially the same, the new trajectories emerging at every crossing point must retain the same state. Trajectories are sequentially cloned, as shown in Figure 14 (b) and, at the end, it is only necessary to restart many new trajectories at t k (eight, in this example), all of them with V (X(t)) = v 2 , and to make account of the incoming and outgoing trajectories at all the crossing points, as shown in Figure 14 (c). The optimal number of thresholds -which is referred to as N T in Tables II, III and IV-is selected after a few pilot runs. It is interesting to point out that, in many cases, this optimal number is such that simulations are plenty of crosses of multiple thresholds simultaneously. Let us analyze, for example, the Lattice 5 × 5 Network, Type A. In this case, due to the possible state of the arcs, V (X(t)) ∈ {0, 2, 4, 6, 8, 10, 12, 16}. The thresholds will, therefore, partition the space [V max , D], which, for the selected model, is [START_REF] Borodina | On splitting scenarios for effective simulation of reliability models[END_REF][START_REF] Lomonosov | On Monte Carlo estimates in network reliability[END_REF]. The minimum size of the step to traverse this space, for trajectories that start at V max and end at D, is 2. This suggests that the only condition under which there will not be crosses of multiple thresholds simultaneously is N T ≤ 4 and only in the extremely fortuitous event in which all the steps are of size 2. Looking at the values of N T in all the tables, it is clear that in most simulations multiple crosses occur.

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

The SFNM is a particular case of a Multi Component System Model (MCSM). In a typical MCSM there is a set of individual components that cooperate with each other, so that the whole system provides some type of service. The individual behavior of the components is known, and the way they interact with each other is also known. But the global response of a MCSM, as a function of the set of components, is usually unknown and sometimes very difficult to determine. The reliability of an MCSM is defined exactly as we have defined the reliability of an SFNM, that is, the probability that the system can deliver service as expected.

For any of these two models, the exact computation of the reliability is impossible or almost impossible when the size of the system, measured by the number of components (arcs in the SFNM), is large. The individual features and properties of the components is not a crucial issue as long as the way they interact with each other is clearly defined and known. The parameter that usually makes the reliability determination extremely difficult (the one that gives the problem its NP-hard character) is the number of components.

In the problem addressed in the paper, the components whose individual behavior is known are the arcs, and the way they cooperate with each other, so that the network can transport a certain amount of flow, is given by its topology. Here, the parameter whose value makes it necessary to improve simulation by variance reduction methods is the number of arcs. It is possible -and relatively simple-to transform the problem and to make the same evaluation considering failures in the nodes, or even considering different sets of nodes generating and consuming flow. But the difficulties in developing variance reduction methods, whatever form the network takes, will be about the same, and the resulting methods will be quite similar, as long as the number of components remain high or extremely high.

These are the reasons why, sometimes, important simplifications are made on the benchmark models used in the development of variance reduction methods, like for example, the use of only two terminal nodes or the fact that nodes are perfect and never fail. However, the change in difficulty resulting from disregarding these simplifications is minimal compared to the intrinsic complexity of the variance reduction methods.

For these types of simplified models, a complexity analysis made in [START_REF] Garvels | The splitting method in rare event simulation[END_REF] claims that, as a function of a rarity parameter, the complexity of the standard simulation methods is exponential in the number of components of the system, whereas the complexity of the Splitting-like methods is polynomial (quadratic). Similar results can be found, for instance, in Chapter 3 of [START_REF] Rubino | Rare event simulation using Monte Carlo methods[END_REF].

IX. CONCLUSIONS

The main proposal of this article is the introduction of a Monte Carlo method designed to estimate the unreliability of highly reliable stochastic flow networks. It is based on an application of Splitting over the well-known Destruction Process, a model in which the components fail at exponentially distributed times. This makes the new procedure, called Splitting over the ML-DP, suitable to work on networks in which the components are non-independent, in terms of the Marshall-Olkin copula.

The paper also proposes a review of techniques based on CP and DP [START_REF] Elperin | Estimation of network reliability using graph evolution models[END_REF], which are graph evolution models originally conceived to work on SNM. The CP was adapted in [START_REF] Cancela | Efficient estimation of stochastic flow network reliability[END_REF], and the DP was adapted here in this article, to SFNM. Due to the fundamentals of the Marshall-Olkin copula model, that refers to failures that take place in exponentially distributed times, it was inevitable to prefer DP over CP to build the new method. And the only way Splitting can be applied to enhance the DP is by partitioning -by means of thresholds-the domain of V (X(t)) and not the temporal domain.

The Marshall-Olkin copula model, so far the most accepted model to work with networks in which the components are not independent, is firstly mentioned in this article over a general setting, and then it is adapted to Splitting over the ML-DP with a high level of detail and providing many implementation guidelines.

The experimental tests show a high level of efficiency, reflected in the order of the estimated relative error bounds reached in very low simulation times. The performance improvement is also reflected in the values of the Speedup, indicating, as expected for any efficient variance reduction method, that the improvements are more significant as the networks become more reliable. While the experiments were performed on two benchmark topologies, and other network configurations could result in other values, the tests employed a wide variety of probabilities for the different capacity distribution functions. This is a most important factor for the method performance, and in this sense we think the experimental results give a clear indication of the performance of the method.

It is interesting and somehow unexpected that the most efficient parameter settings are those that cause crosses of multiple thresholds simultaneously.

Given that the article introduces a new proposal, it is expected that many of its features will require further analysis. One of the topics to be considered as part of our future work is a theoretical analysis of robustness, for example, in terms of the measures discussed in [START_REF] Cancela | New measures of robustness in rare event simulation[END_REF].

Although the variance of the Splitting estimator is extremely difficult to compute exactly in these contexts, another interesting subject to consider, at least experimentally, is to check if there are conditions under which the model verifies the Bounded Relative Error property. The design of techniques to determine -or, at least, to approximate-the optimal number of thresholds and the multiplying Splitting factor are also of interest.

The capacity distributions that we worked with in this article are, in all cases, discrete. However, it is possible, and this is another interesting subject for future work, to adapt Splitting over the ML-DP to the case of continuous capacity distributions.

Fig. 1 .

 1 Fig. 1. Ordering sequence of repairing times.

Fig. 2 .

 2 Fig. 2. The ML-CP mechanism for assigning values to Xi(t)

Fig. 3 .

 3 Fig. 3. Increasing of V (X(t)) as the arcs are repaired.

Fig. 5 .

 5 Fig. 5. Splitting over the ML-CP

Fig. 6 .

 6 Fig. 6. The ML-DP mechanism for assigning values to Xi(t)

Fig. 7 .

 7 Fig. 7. Decreasing of V (X(t)) as the arcs fail.

Fig. 9 .

 9 Fig. 9. Splitting over the ML-DP

Fig. 11 .

 11 Fig. 11. Fishman Network

Fig. 13 .

 13 Fig. 13. Cloning a trajectory after crossing a single threshold

Fig. 14 .

 14 Fig. 14. Cloning a trajectory after crossing multiple thresholds

  unbiased estimator of ζ where: A=empty means SMC, A=P means PMC and A=S means Splitting. B is either CP, DP, ML-CP or ML-DP.

		the flow demand in SFNM
	ζ	network unreliability
	ζ	unbiased estimator of ζ
	ζ A B
	N	sample size of a Monte Carlo experiment
	N T	number of thresholds in Splitting

i i-th threshold in both, ML-CP and ML-DP α i probability that a trajectory started at threshold i reaches threshold i+1 underneath D, in ML-CP α i estimator of α i β i

  = 1 -p 0 -p 2 -p 4 10 w.p. p 10 = (1 -p 1 -p 0 ) × 2/3 I s 4 w.p. p 4 = p 0 /100 9 w.p. p 9 = (1 -p 1 -p 0 ) × 1/3 2 w.p. p 2 = p 0 /10 1 w.p. p 1 = p 0 /2 0 w.p. p 0 0 w.p. p 0 I m 8 w.p. p 8 = 1 -p 0 10 w.p. p 10 = 1 -p 0 0 w.p. p 0 = p 0 /10 0 w.p. p 0 = p 0 /10

	TABLE I
	CAPACITY DISTRIBUTION OF THE NETWORKS USED IN THE EXPERIMENTAL TESTS
	Type A	Type B
	8 w.p. p 8

TABLE II EFFICIENCY

 II EVALUATION OF THE UNRELIABILITY OVER NETWORKS TYPE A

	Network	p 0	N	N T	ζ S ML-DP	RE	t S ML-DP	Speedup
		0.001	8.00E+03	4	1.52E-02 2.98%	0.20 sec.	1
	Fishman	0.0001	4.50E+05	7	1.44E-04 3.03%	1.68 sec.	9
	D = 15	0.00001	6.00E+06 13 1.37E-06 3.00%	15.87 sec.	98
		0.000001 1.60E+08 18 1.40E-08 2.88%	6.46 min.	400
		0.0000001 2.10E+09 24 1.38E-10 2.82% 123.29 min.	2,186
		0.001	1.00E+05	6	4.58E-03 1.17%	22.93 sec.	1
	Lattice 5 × 5	0.0001	1.00E+06	9	4.70E-05 2.55%	17.01 sec.	14
	D = 7	0.00001	1.00E+07 15 4.78E-07 2.81%	87.82 sec.	193
		0.000001 1.00E+08 22 4.71E-09 3.14% 16.97 min.	1,411
		0.0000001 1.00E+09 27 4.74E-11 4.96% 186.67 min.	5,097

TABLE III EFFICIENCY

 III EVALUATION OF THE UNRELIABILITY OVER NETWORKS TYPE B

	Network	p 0	N	N T	ζ S ML-DP	RE	t S ML-DP	Speedup
		0.01	1.00E+04	8	5.80E-03 2.48%	0.64 sec.	2
		0.001	5.00E+05 14 5.82E-05 1.13%	25.93 sec.	17
	Fishman	0.0001	1.00E+06 18 5.41E-07 2.42%	76.28 sec.	156
	D = 15	0.00001 5.00E+06 26 5.98E-09 2.56% 11.58 min.	1,223
		0.000001 2.50E+07 30 5.91E-11 4.03% 93.35 min.	5,800
		0.01	1.50E+04	8	2.00E-03 2.66%	6.37 sec.	2
		0.001	1.00E+05 15 1.85E-05 2.62%	58.48 sec.	20
	Lattice 5 × 5	0.0001	1.00E+06 19 1.87E-07 2.90% 10.11 min.	118
	D = 7	0.00001 6.00E+06 26 1.91E-09 2.96% 116.06 min.	1,083
		0.000001 2.50E+07 27 1.90E-11 5.65% 386.32 min.	9,321

ACKNOWLEDGMENT

This work was partially supported by the Math AmSud project "RareDep" (Rare events analysis in multi-component systems with dependent components), 2019-2020-2021.


I Ajayi

F Kaddour
email: fatima.kaddour@anfr.fr

L Mroueh
email: lina.mroueh@isep.fr

Impact of Imperfect Channel State Information on Physical Layer Security by Precoding and Diversity

Keywords: Imperfect CSI, Adaptive Matched-Filter Precoding, Diversity, PLS, OFDM, FDD, TDD

Physical layer security (PLS) is an emerging paradigm that makes use of wireless channel characteristics to provide security. Many PLS schemes require knowledge of the channel state information (CSI). However, CSI is usually imperfect due to factors such as noisy feedback channels, channel estimation errors and outdated CSI. In this paper, we investigate the impact of imperfect CSI on the secrecy and error rate performances of a PLS scheme that combines adaptive matched-filter (MF) precoding and diversity in orthogonal frequency division multiplexing (OFDM). Particularly, we derive the secrecy capacity and error rate expressions for the legitimate and eavesdropper's channels under imperfect CSI assumption. The impact of the imperfect CSI is studied via theoretical as well as numerical techniques in frequency-selective Rayleigh fading wiretap channel. The analysis is done in both frequency division duplex (FDD) and time division duplex (TDD) modes.

I. INTRODUCTION

Physical layer security (PLS) techniques take advantage of the characteristics of wireless channels such as noise, fading, interference, diversity, dispersion, etc to ensure an intended receiver successfully decodes a transmitted data while preventing an eavesdropper from doing so [START_REF] Hamamreh | Classifications and Applications of Physical Layer Security Techniques for Confidentiality: A Comprehensive Survey[END_REF]. The use of artificial noise injection to enhance security by degrading the eavesdropper's channel can be seen in [START_REF] Hong | Artificial-Noise-Aided Secure MIMO Wireless Communications via Intelligent Reflecting Surface[END_REF]- [START_REF] Goel | Guaranteeing secrecy using artificial noise[END_REF]. Beamforming, a technique that involves adjusting the signal amplitudes and phases to form a strong beam towards a direction of interest, has also been used for PLS [START_REF] Liang | A Physical Layer Security Approach Based on Optical Beamforming for Indoor Visible Light Communication[END_REF]- [START_REF] Song | Beamformer Design for Physical Layer Security in Dual-Polarized Millimeter Wave Channels[END_REF]. The use of channel coding techniques for security is addressed in [START_REF] Nooraiepour | Randomized serially concatenated ldgm codes for the gaussian wiretap channel[END_REF], [START_REF] Masuda | Low Complexity Decoding of Downlink Chaos NOMA Scheme with Physical Layer Security[END_REF]. Optimal power allocation as a PLS technique was worked on in [START_REF] Chen | Power-efficient joint resource allocation for multiuser wiretap OFDM channels[END_REF], [START_REF] Karachontzitis | Security-aware maxmin resource allocation in multiuser OFDMA downlink[END_REF]. In wireless communications, diversity is used to mitigate the effects of multipath fading [START_REF] Alireza | Impact of a Time-Varying Rician Fading Channel on the Performance of Alamouti Transmit Diversity Technique[END_REF]- [START_REF] Damen | Diagonal Algebraic Space Time Block Codes[END_REF]. Some works in the literature address the use of diversity in providing PLS [START_REF] Allen | Secure Space-Time Block Coding without Transmitter CSI[END_REF]- [START_REF] Yusuf | Enhancing Physical-Layer Security In Wireless Communications Using Signal Space Diversity[END_REF].

Precoding is an interesting way to enhance PLS and we take a look at some of the research works in this direction. In [START_REF] Cruz | Wireless Physicallayer Security Using Precoding and an Active Eavesdropper[END_REF], the authors proposed the use of zero-forcing (ZF) and the minimum means squared error (MMSE) precoders for PLS in a wiretap system where there is spatial decorrelation between the main channel and the wiretap channel. The performance of the eavesdropper is completely degraded because it carries out blind equalization as there is no channel state information (CSI) leakage to it. A precoded orthogonal space-time block coding (POSTBC) in a multiple-input single-output (MISO) setup that minimizes the error rate only at the legitimate receiver was investigated in [START_REF] Hamamreh | A practical physical-layer security method for precoded OSTBC-based systems[END_REF]. In the paper, the authors also proposed a new technique called "precoding along with partial pre-equalizing" (PCPPE) that further improves the security of the system. This is achieved by using a new precoder design that is composed of both the original precoder and a newly designed unitary matrix. In [START_REF] Maeng | Precoder Design for mmWave UAV Communications with Physical Layer Security[END_REF], a precoder that incorporates the CSI of the eavesdropper in a millimetre-wave (mmWave) unmanned aerial vehicles (UAV) system was proposed. Symbol-level precoding to counteract learningassisted eavesdropping in downlink multiuser-MISO system was worked on in [START_REF] Mayouche | Learning-Assisted Eavesdropping and Symbol-Level Precoding Countermeasures for Downlink MU-MISO Systems[END_REF]. The authors in [START_REF] Nguyen | Performance Analysis of Matched-Filter Precoded MISO-OFDM Systems in the Presence of Imperfect CSI[END_REF] studied the impact of imperfect CSI on the mean square error (MSE) performance of a downlink MISO orthogonal frequency division multiplexing (OFDM) systems using a matched-filter (MF) precoder. Some of the works in the literature that employ MF precoding in multi-antenna scenarios can be seen in [START_REF] Al-Nahari | Physical layer security using massive multiple-input and multiple-output: passive and active eavesdroppers[END_REF]- [START_REF] Nguyen | Secure Massive MIMO With the Artificial Noise-Aided Downlink Training[END_REF]. These are in contrast to our work that considers a single antenna multi-carrier scenario.

Studies on the impact of imperfect CSI in different PLS schemes can be seen in [START_REF] Yang | Secure Massive MIMO Under Imperfect CSI: Performance Analysis and Channel Prediction[END_REF]- [START_REF] Dung | Secrecy Performance Analysis for MIMO Relay System with Transmit/Receive Antenna Selection under Imperfect CSI[END_REF]. The authors in [START_REF] Yang | Performance Analysis of Secure Communication in Massive MIMO with Imperfect Channel State Information[END_REF] studied secure communications in a multiuser massive multiple input multiple output (MU-mMIMO) system with imperfect CSI due to outdated CSI and channel estimation errors. The obtained results showed a significant reduction in secrecy capacity due to imperfect CSI. In [START_REF] Kumar | Physical layer security for multiuser multi-eavesdroppermulti-input multi-output (MIMO) system in the presence ofimperfect feedback[END_REF], the impact of imperfect CSI in an MU-MIMO system that combines selection transmission at the transmitter and maximum ratio combining (MRC) at the receivers was studied. It was observed that PLS performance in terms of probability of non zero secrecy capacity and secrecy outage probability degrades with a rise in imperfect CSI. In [START_REF] Dung | Secrecy Performance Analysis for MIMO Relay System with Transmit/Receive Antenna Selection under Imperfect CSI[END_REF], the authors analysed the secrecy performance of a MIMO relay system under imperfect CSI. They concluded that the saturated minimum secrecy outage probability and maximum secrecy capacity depends on the severity of the imperfect CSI.

In this paper, we study the impact of imperfect CSI on a PLS scheme that combines repetition scheme for diversity and adaptive MF precoding in OFDM. The study investigates the use of MF precoding in a single antenna frequency selective channel scenario. Precoding transfers most of the signal processing to the transmitter and the receiver only carries out symbol detection. This means that the accuracy of the CSI at the transmitter will have a greater impact on the system's performance compared to the accuracy of the CSI at the receiver. We assume that imperfect CSI with equal estimation error variance is available at the three nodes: sender (Alice), legitimate receiver (Bob) and eavesdropper (Eve). Repetition scheme is adopted to provide diversity gain to Bob but the MF precoder ensures that Eve loses the diversity gain. This leads to a security gap between Bob and Eve measured in terms of bit error rate (BER) and secrecy capacity. We summarize our major contributions as follows:

• We study the use of adaptive MF precoding in a single-antenna multi-carrier transmission mode, combined with repetition diversity.

• We derive the expressions for the received signals, conditional secrecy capacity and conditional QPSK BER for Bob and Eve under imperfect CSI conditions.

• These derivations are verified numerically and theoretically.

• We provide an analysis of this imperfect CSI on the security performance of the system.

• We consider frequency division duplex (FDD) and time division duplex (TDD) modes. In FDD mode, there is CSI feedback between Alice and Bob and this CSI leaks to Eve. In this mode, we assume that Eve is passive but fully aware of the CSI of the wiretap and main channels. This is the worst case for security. Channel reciprocity is adopted in TDD mode and Eve only carries out blind equalization.

This paper is organized as follows. In Section II, we describe the system model. Section III is dedicated to the impact of imperfect CSI. Performance evaluation is discussed in Section IV and conclusions drawn in Section V.

Notations: Vectors are denoted by bold letters, whereas individual vector elements are denoted by normal letters. Norm-2 is defined by . . Conjugate and absolute value are symbolized by (.) * and |.| respectively. II. SYSTEM MODEL Fig. 1 shows the block diagram for the OFDM communication system. Alice transmits information signals to Bob in the presence of Eve. At Alice, the instantaneous CSI is used to design an adaptive MF precoder unique to all subcarriers. After this, diversity encoding takes place on two uncorrelated subcarriers before OFDM modulation and transmission. An advantage of this scheme is a reduced receiver complexity compared to the legacy repetition scheme since channel equalization is done at Alice. The receiver simply adds the two received signals for every diversity pair and carries out symbol detection. It is important to note that other transmitter diversity approaches such as diagonal algebraic space-time (DAST) coding [START_REF] Damen | Diagonal Algebraic Space Time Block Codes[END_REF] and Alamouti scheme [START_REF] Alamouti | A Simple Transmit Diversity Technique for Wireless Communications[END_REF] can also be adopted with minimal modifications to match the specifics of the schemes.

We adopt a scenario in which Bob and Eve are located at separate locations. Hence, their channels exhibit uncorrelated propagation in a rich scattering environment and the CSI of both channels are uncorrelated. To simplify the expressions and without loss of generality, the input-output relations are specified in frequency domains only. We model the imperfect CSI on the i th subcarrier as given in [START_REF] Hyadi | An Overview of Physical Layer Security in Wireless Communication Systems With CSIT Uncertainty[END_REF]:

H i = √ 1 -ε Ĥi + √ ε Hi (1)
where H is the actual channel gain without errors and Ĥ is the imperfect channel gain with errors.

The estimation error H is a zero-mean unit variance complex Gaussian random variable CN (0, 1), and it is independent of H. The variance of the estimation error is denoted by ε ∈ [0, 1].

By employing the legacy repetition scheme for diversity, Alice transmits two precoded data symbols on two uncorrelated subcarriers in the same OFDM block. For every diversity pair (i, j), the received signals after OFDM demodulation are

r i = H i x i + n i , r j = H j x j + n j (2)
where x i and x j are the precoded data symbols while n i and n j are the complex additive white Gaussian noise (AWGN) terms for the subcarrier pair (H i , H j ).

x i = P i s, x j = P j s

and P is the adaptive MF precoder that is designed according to the instantaneous CSI of the legitimate receiver. It is expressed as

P i = √ 2H * i H i, j , P j = √ 2H * j H i, j (4)
where H i, j is the vector representing the subcarrier pair, (H i ,H j ), of the main channel under perfect CSI assumption. For a total of N available subcarriers, there are N/2 subcarrier pairs. By substituting (3) and ( 4) in ( 2) and simply summing the received signals on

From ( 5), the instantaneous signal to noise ratio (SNR) at Bob is given as

γ ib1 = H i, j 2 γ a ( 6
)
where γ a is the average SNR at the receiver.

The calculation of the BER for QPSK constellation is readily available in the literature when the decision variables are Gaussian random variables [START_REF] Cho | On the general BER expression of one and two dimensional amplitude modulation[END_REF]

BER(SNR) = 1 2 erfc SNR 2 (7)
Therefore, by conditioning on the set of variables H i and H j , we can obtain the conditional QPSK error probability corresponding the subcarrier pair (i,j) at Bob

BER b1 (i, j) H i ,H j = 1 2 erfc γ ib1 2 (8)
The final BER is obtained by averaging the conditional BER on the variables H i H j for all subcarrier pairs (i,j).

Similarly, when the CSI is perfect, the received signal at Eve, the instantaneous SNR and conditional QPSK BER can be respectively written as

se1 = √ 2 H ie H * i + H je H * j H i, j s + n ie + n je (9) γ ie1 = αα * H i, j 2 γ a ( 10
) BER e1 (ie, je) H ie ,H je = 1 2 erfc γ ie1 2 (11)
where

α = H ie H * i + H je H * j (12)
H ie, je is a vector representing a wiretap channel subpair (H ie , H je ) when the CSI is perfect. n ie and n je are the complex AWGN components for the subcarrier pair.

Secrecy capacity is the positive difference between the main channel capacity and eavesdropper's channel capacity [START_REF] Barros | Secrecy Capacity of Wireless Channels[END_REF]. A positive value means secrecy is achievable and a zero implies there is no secrecy guarantee. We measure secrecy capacity in bps/Hz (or bits/channel use). Similar to the error rate performance, the secrecy capacity also is obtained by averaging the conditional secrecy capacity on the channel gain variables H i and H j for all subcarrier pairs (i,j). The conditional channel capacities of Bob and Eve and the conditional secrecy capacity are respectively expressed as

C b1 (i, j) H i ,H j = 1 2 log 2 (1 + γ ib1 ) (13)
C e1 (ie, je) H ie ,H je = 1 2 log 2 (1 + γ ie1 ) (14)
C s1 H i ,H j ,H ie ,H je = C b1 -C e1 , if γ ib1 > γ ie1 0, otherwise. (15)
The factor of half in ( 13) and ( 14) is because only half of the available bandwidth is used for the transmission. For N available subcarriers, N/2 unique symbols are transmitted.

III. IMPACT OF IMPERFECT CSI

To achieve secure communication in the presence of an eavesdropper, the adaptive MF precoder ensures that only Bob maintains the diversity gain provided by the repetition scheme. Since the eavesdropped signals are precoded using the legitimate user's CSI, Eve will lose the diversity gain and experience a higher error rate and lower channel capacity than Bob. This degradation of Eve compared to Bob guarantees secure communication in the presence of Eve as Bob can decode transmitted symbols at a higher SNR compared to Eve.

To investigate the impact of imperfect CSI on the system, we assume that imperfect CSI with equal estimation noise variance is available at Alice, Bob and Eve. The system will no longer be optimal for Bob because the precoded symbols do not maximize the SNR anymore as is expected in MF precoding. As this estimation error variance increases, the intersymbol interference increases and Bob loses the diversity gain and thus the security gap over Eve decreases.

The received signal at Bob under imperfect CSI is expressed as ri = H i xi + n i , r j = H j x j + n j [START_REF] Akuon | Secure Signal and Space Alamouti Scheme[END_REF] where xi = Pi s, x j = Pj s [START_REF] Yusuf | Enhancing Physical-Layer Security In Wireless Communications Using Signal Space Diversity[END_REF] and P, the adaptive MF precoder that is designed according to the imperfect instantaneous CSI

Pi = √ 2 Ĥ * i Ĥi, j , Pj = √ 2 Ĥ * j Ĥi, j (18)
where Ĥi, j is the vector representing a particular subcarrier pair, ( Ĥi , Ĥj ), of the main channel under imperfect CSI assumption. By substituting ( 1), ( 17) and ( 18) in ( 16), the received signal at Bob is shown in [START_REF] Hamamreh | A practical physical-layer security method for precoded OSTBC-based systems[END_REF] below

From [START_REF] Hamamreh | A practical physical-layer security method for precoded OSTBC-based systems[END_REF], the instantaneous SNR at Bob can be written as

γ ib2 = (1 -ε) Ĥi, j 2 γ a εγ a + 1 (20)
Similar to earlier derivations, the conditional QPSK BER at Bob is given as

BER b2 (i, j) H i ,H j = 1 2 erfc γ ib2 2 (21)
Following the same convention, the received signal (shown in ( 22) below), instantaneous SNR and conditional QPSK BER at Eve are respectively expressed as

γ ie2 = (1 -ε)(β β * )γ a (εγ a + 1) Ĥi, j 2 (23)
BER e2 (ie, je) H ie ,H je = 1 2 erfc γ ie2 2 (24)
where

β = Ĥie Ĥ * i + Ĥje Ĥ * j (25)
Ĥie and Ĥje are the channel frequency response for the wiretap channel (channel between Alice and Eve) when the CSI is imperfect.

In this case, the conditional channel capacities for Bob and Eve and the conditional secrecy capacity can be respectively expressed as IV. PERFORMANCE EVALUATION Table I is the summary of the parameters and their values for the simulations carried out. We adopt an OFDM system with a total of N = 128 subcarriers. The symbols are QPSK modulated. Firstly, we assume a perfect CSI (ε = 0), then imperfect CSI is assumed with error variances of ε = 0.02 and 0.10 respectively.

C b2 (i, j) Ĥi , Ĥj = 1 2 log 2 (1 + γ ib2 ) (26)
C e2 (ie, je) Ĥie , Ĥje = 1 2 log 2 (1 + γ ie2 ) (27)
C s2 Ĥi , Ĥj , Ĥie , Ĥje = C b2 -C e2 , if γ ib2 > γ ie2 0, otherwise. (28)
For a fair comparison, we compare the error rate and secrecy performances of Bob and Eve under the assumptions of the same average SNR and estimation error variances at both of them. Due to spatial decorrelation and rich scattering, the main channel and wiretap channel are uncorrelated.

A. Bit Error Rate (BER)

is a Quality of Service (QoS) related metric that can also be used to analyse the security of a system. The channel with higher BER is less secure and more degraded compared to the channel with lower BER under similar conditions and the same information signal. The difference in the BER is a measure of the security gap in the system [START_REF] Hamamreh | Classifications and Applications of Physical Layer Security Techniques for Confidentiality: A Comprehensive Survey[END_REF]. Fig. 2 shows the BER performance of the system. We present the numerical and theoretical results.

When the CSI is perfect (ε = 0), it can be observed that Bob outperforms Eve. Bob maintains a diversity of 2 but Eve loses the diversity gain in FDD mode and is fully degraded in TDD mode. In FDD mode, Eve is aware of the main channel instantaneous CSI and is able to use this instantaneous CSI for equalization. However, since the eavesdropped symbols have been precoded with an uncorrelated CSI, it leads to loss of diversity gain compared to Bob. The higher BER for Eve compared to Bob is an indication of the security gap in the system. In TDD mode, the main channel CSI is not available to Eve, Eve carries out blind equalization and is completely degraded.

Next, we look at the impact of imperfect CSI on the BER performance. When ε = 0.02, we observe in Fig. 2 a slight increase in the error rates of Bob and Eve. However, the system is still secure as the gap between Bob and Eve remains significant. As the error variance increases to ε = 0.1, the error rate significantly increases at both Bob and Eve. Bob completely loses the diversity gain. They both begin to exhibit an error floor around an average SNR of 20dB. The system is highly sensitive to channel estimation errors.

B. Secrecy Capacity

In Fig. 3, we plot the secrecy capacity of the system against the average SNR. As expected, we sb2 Average SNR(dB) see that the secrecy capacity is highest when there is no estimation error (ε = 0), an ideal assumption. However, when the CSI is imperfect, the secrecy capacity reduces as the variance of the estimation error increases. The higher the magnitude of the estimation error variance, the lower the overall secrecy capacity in the system. This is in agreement with the BER performance. At all estimation error variances studied, the secrecy capacity remains positive. This means that although the imperfect CSI causes performance degradation at Bob and Eve, the channel capacity of Bob remains higher than the channel capacity of Eve.

= 2(1 -ε) Ĥi, j s + √ 2ε( Hi Ĥ * i + Hj Ĥ * j ) Ĥi, j s + n i + n j (19) se2 = √ 1 -ε( Ĥie Ĥ * i + Ĥje Ĥ * j ) + √ ε( Hie Ĥ * i + Hje Ĥ * j ) Ĥi, j / √ 2 s + n ie + n je (22

V. CONCLUSION

We studied the impact of imperfect CSI on a system that combines the popular MF precoding and repetition scheme to provide PLS. The adaptive MF precoder is designed based on the instantaneous CSI of the main channel and ensures that only the legitimate receiver experiences the diversity gain from the repetition scheme. The precoder deteriorates the error rate and secrecy performance of the eavesdropper. Since this scheme depends on CSI estimation, imperfect CSI has a significant impact on the secrecy and error rate performance of the system. We derived the expressions for the received signals, secrecy capacity and BER for the legitimate receiver and the eavesdropper under perfect and imperfect CSI conditions. Using numerical and theoretical simulations, we observe that as the CSI estimation error variance increases, the secrecy capacity reduces and BER significantly increases leading to system degradation. To mitigate the effect of imperfect CSI and improve the security in the system, a future work on this paper will be the proposal of low resource consuming real time CSI prediction algorithms. Alice will then use the predicted CSI for adaptive precoding. When the CSI prediction is highly accurate, Bob's performance will be optimal and near the perfect CSI performance while the effects of Imperfect CSI will remain at Eve. This work has been supported by ISEP and ANFR.

Figure 1 :

 1 Figure 1: OFDM system model employing diversity and precoding for security

Figure 2 :Figure 3 :

 23 Figure 2: Bob and Eve error rate performances with imperfect CSI at Alice, Bob and Eve, ε = 0, 0.02, 0.1

Table I :

 I Simulation Parameters

	Parameter	Values
	Modulation Scheme	QPSK
	Error Variance(ε)	0, 0.02, 0.1
	Number of subcarriers	128

  )

	0	2	4	6	8	10	12	14	16	18	20

8th International Conference on Electrical Engineering, Computer Science and Informatics (EECSI) 20-21 October 2021


What is metabarcoding?

• Some DNA fragments are highly conserved within a species and variable between species. These are the genetic markers or barcodes.

• Metabarcoding, by identifying barcodes through sequencing, allows blind identification of all species present in a sample at once. • Implemented an optimized process to performed long read metabarcoding on 16S and 16S-23S with PacBio HiFi technology :

• 16S-23S region is more informative than full length 16S

• Problem remains with the database exhaustivity that require time to be completed

Conclusions

••

  Gene encoding the 16S subunit of ribosomal RNA • Highly conserved gene in all bacteria • Highly conserved regions (for primers) interleaved with variable regions (bacteria identification) • Choice of the regions depends on the community to be analyzed Ribosomal RNA: bacterial markers Comparison of different hypervariable regions of 16S rRNA Number of taxa identified at each taxonomical using individual and concatenated hypervariable regions of the 16S rRNA gene. • Analysis of a combination of sequences from V2 to V9 regions identified more taxa. Pig faeces (hard to extract, inhibitors) 32 faeces samples from 16 individuals (GenPhySE) exposed or not to mycotoxine Fumonisin B1 • Short read (V3-V4) data available validation on short read (V3-V4 region)  Slight biases when sequencing mock DNA  Zymo DNA miniprep kit gave the best results, choosen for metabarcoding (Mock bact) https://www.pacb.com/wp-content/uploads/Multiplexed-Amplicon-Library-Preparation-Using-SMRTbell-Express-Template-Prep-Kit-2.0-%E2%80%93-Customer-Training.pdf  Protocols have been tested with 16S and 16S-23S barcodes on mock and pig faeces samples https://www.pacb.com/wp-content/uploads/Multiplexed-Amplicon-Library-Preparation-Using-SMRTbell-Express-Template-Prep-Kit-2.0-%E2%80%93-Customer-Training.pdf  Protocols have been tested with 16S and 16S-23S barcodes on mock and pig faeces samples www.pacb.com/wp-content/uploads/Multiplexed-Amplicon-Library-Preparation-Using-SMRTbell-Express-Template-Prep-Kit-2.0-%E2%80%93-Customer-Training.pdf  Protocols have been tested with 16S and 16S-23S barcodes on mock and pig faeces samples A longer barcode allows to describe more precisely a simple community.



  This pipeline allows to perform assembly, taxonomic affiliation of contigs, functional annotation and binning of contigs. /forgemia.inra.fr/genotoul-bioinfo/metagwgs •Genomic signature •Differential contig depth  The Genotoul bio-info platform has developed a metagenomic shotgun pipeline for the analysis of Illumina and PacBio HiFi sequencing data. The PacBio HiFi part was developed during SeqOccIn project.

http://bioinfo.genotoul.fr @BioinfoGenotoul

Bioinformatics process for metabarcoding

Sequencing Demultiplexing

Each index is assignes to its original sample

Denoising and clustering

Sequences are groupes by similarity. Amplification and sequencing errores are masked.

Taxonomic affiliation

Each group is affiliated with a taxon through a reference database.

Abundance table

Detect and count the taxa present in each sample.