# Exact reliability optimization for series-parallel graphs using convex
envelopes

Javiera Barrera  Faculty of Engineering and Sciences, Universidad Adolfo
Ibáñez, Santiago, Chile<EMAIL_ADDRESS>,  Eduardo Moreno  Faculty of
Engineering and Sciences, Universidad Adolfo Ibáñez, Santiago, Chile
<EMAIL_ADDRESS>,  Gonzalo Muñoz  Institute of Engineering Sciences,
Universidad de O’Higgins, Rancagua, Chile<EMAIL_ADDRESS>and  Pablo
Romero  Faculty of Engineering, Universidad de la República, Montevideo,
Uruguay  Faculty of Natural and Exact Sciences, Universidad de Buenos Aires,
Argentina<EMAIL_ADDRESS>
###### Abstract.

Given its wide spectrum of applications, the classical problem of all-terminal
network reliability evaluation remains a highly relevant problem in network
design.

The associated optimization problem—to find a network with the best possible
reliability under multiple constraints—presents an even more complex
challenge, which has been addressed in the scientific literature but usually
under strong assumptions over failures probabilities and/or the network
topology.

In this work, we propose a novel reliability optimization framework for
network design with failures probabilities that are independent but not
necessarily identical. We leverage the linear-time evaluation procedure for
network reliability in the series-parallel graphs of Satyanarayana and Wood
[29] to formulate the reliability optimization problem as a mixed-integer
nonlinear optimization problem. To solve this nonconvex problem, we use
classical convex envelopes of bilinear functions, introduce custom cutting
planes, and propose a new family of convex envelopes for expressions that
appear in the evaluation of network reliability. Furthermore, we exploit the
refinements produced by spatial branch-and-bound to locally strengthen our
convex relaxations. Our experiments show that, using our framework, one can
efficiently obtain optimal solutions in challenging instances of this problem.

###### Key words and phrases:

network reliability, reliability optimization, series-parallel graphs, convex
envelopes, nonlinear optimization

The research leading to these results received funding from grants ANID
Fondecyt Regular 1200809 (J.B., E.M.); ANID Fondecyt Iniciación 11190515
(G.M.); ANID ANILLO ACT192094 (E.M., J.B.); Program Math Amsud 19-MATH-03
(J.B.,P.R.)

##  1\. Introduction

Network design to optimize the reliability of a network is one of the most
classic problems in the optimization literature, with applications in many of
the aspects of our lives and certainly in many more in the future. In the
recent special issue to celebrate the 50-year anniversary of the _Networks_
journal, a beautiful and exhaustive review by Brown et al. [9] was published,
which covers the different notions, theory and applications of network
reliability, along with future directions. In a similar way, Pérez-Rosés [25]
revisited the 60 years of this problem, initially started by E.F. Moore and
C.E Shannon [22] in 1956. This paper proposes a new methodology to advance in
one of the most classical and defying problems in this area, the exact
reliability optimization of a network.

Consider a simple undirected graph $G=(V,E)$ representing a network whose
edges fail according to a random probability—these failures are assumed to be
independent. We consider one of the most classical notions of reliability: the
_all-terminal reliability_ of $G$, which is defined as the probability that
$G$ remains connected.

The exact reliability evaluation belongs to the class of $\mathcal{NP}$-hard
problems [3], even when only two terminals are required to be connected [31]
or under identical link failure probabilities [27]. However, for certain
specific families of graphs, the problem can be solved in polynomial time. In
the particular case of series-parallel graphs, Satyanarayana and Wood [29]
provide a set of reliability-preserving reductions to simplify the graph and
its structure that make it possible to compute the reliability of a series-
parallel network in linear time. This work is our starting point. This type of
reduction can also be applied to general graphs to reduce the size and
complexity of the problem. For example, a recent work [17] extended these
results to obtain a parametrized algorithm for computing the reliability of
graphs with small treewidth in linear time. Nonetheless, this approach has not
been thoroughly exploited in the literature.

As Brown et al. [9] note in their recent review, one of the most relevant and
less-traveled paths is the optimization of network reliability. That is, given
a limited set of resources (e.g., number of edges), how should one select a
subgraph of $G$ that maximizes its reliability. Given the hardness of
computing the exact reliability of a network, using this metric as the
objective function for an optimization problem can be highly impractical in
general. However, since the computation of the reliability in certain families
of graphs is tractable and due to the sharp progress we have witnessed in
mixed-integer nonlinear optimization technology, an efficient method for
reliability optimization for particular classes of graphs is plausible.

A survey on reliability optimization is provided in [7]. Most previous works
have assumed either identical link costs or independent link failures with
identical probability; in these cases, the optimal topologies are highly
symmetrical [34]. A foundational work on this problem for the all-terminal
reliability problem was carried out by Boesch [6]. In this work, a fixed
budget is considered under identical costs (i.e., the network has a precise
number of links), and identical independent link failures with probability
$\rho$ are assumed. The goal is to design a graph whose all-terminal
reliability is maximized in a uniform sense, i.e., over the whole compact set
$\rho\in[0,1]$, among all the graphs with the same number of nodes and links.
These graphs are known as _uniformly most-reliable graphs_ or UMRGs. It is
known that UMRGs must have the maximum tree number and maximum connectivity,
providing evidence of the symmetry of the optimal graphs under the strong
assumptions of independence and identical costs. Nevertheless, several
historical conjectures regarding the construction and existence of UMRGs
remain open [2]. In the restricted class of series-parallel graphs, their
UMRGs are characterized in [23].

One way of tackling the difficulties associated with exact reliability
computation is to use meta-heuristics to provide a _good_ solution for the
problem; this approach was highly studied during the 1990s [1, 15, 12]. Later
work on this subject includes hybrid ACO [14] and self-tuning heuristics [13].
More recently, [24] combines heuristic techniques with branch-and-bound
methods. However, the main drawback of these techniques is that heuristics
cannot provide guarantees on the optimality of their solutions. A different
approach to overcome this problem is to use simulations of the failure process
and embed these sampled scenarios into an optimization model, known as the
sample average approximation method [19]. This approach allows one to consider
failure models without assuming independent or identical failures
probabilities but generates an approximated formulation. For example, [30]
provides a powerful methodology to minimize the cost of a network meeting two-
terminal reliability constraints, which considers probabilistic cuts in the
branch-and-bound tree. Similarly, [4] presents a reliability optimization
model with dependent link failures ruled by the Marshall-Olkin copula [8, 20].
In [26], a stochastic network flow model using scenarios for the two-terminal
reliability case is provided.

Since the nonlinear functions representing the reliability of the problem are
neither convex nor concave functions, straightforward mathematical
optimization formulations seem to be of little use in reliability optimization
problems. Nonetheless, one can rely on convexification techniques to
approximate the problem with a tractable alternative and exploit spatial
branch-and-bound to further refine such approximations. This is a key
component of state-of-the-art mixed-integer nonlinear programming (MINLP)
technology and, to the best of our knowledge, it has not been exploited for
the general network reliability problem. For example, convex envelopes of
functions can be efficiently computed under special cases, which yield strong
convex relaxations. This idea is used in [33] to optimize the design of a
reliable chemical plant, which can be represented as a series graph. Similar
techniques have been used for other network-related problems, for example, on
AC optimal power flow problems [10].

Our paper follows this previous idea. We consider the optimization reliability
problem with independent failure but with different failure probability among
edges. In this setting, the reduction techniques remain valid [29]. Based on
these reliability-preserving serial-parallel reductions, a convex MINLP
formulation for the reliability optimization problem is obtained for series-
parallel graphs. In this formulation, each reduction generates new constraints
(a linear number). We provide tight convex envelopes for the functions
appearing from the reduction process, which, combined with the refinements
carried out in spatial branch-and-bound trees, allows us to obtain the (exact)
optimal solution efficiently. This idea could also be extended to other
families of graphs: for instance, one can rely on series-parallel reductions
to decrease the size of a problem.

This article is organized as follows. Section 2 presents series and parallel
reductions that preserve all-terminal reliability, following the work of
Satyanarayana and Wood [29]. Section 3 presents the main contributions of this
work. Specifically, a nonlinear and nonconvex formulation of the reliability
optimization problem is introduced in Subsection 3.1, considering series-
parallel reductions. Subsection 3.2 introduces convex envelopes associated
with the series-parallel reductions, using classical McCormick envelopes [21]
and a novel envelope for series-type reductions (Theorem 1). Further
improvements to the resulting convex optimization problem are addressed in
Subsection 3.5. The computational effectiveness of our proposal is studied in
Section 3.4. Finally, Section 4 presents concluding remarks and directions for
future work.

##  2\. Definitions and reliability-preserving reductions

Consider an undirected graph $G=(V,E)$. Nodes are perfectly reliable, but
links may fail with independent probabilities $q_{e}$ for $e=1,\ldots|E|$. Let
us denote by $p_{e}=1-q_{e}$ the elementary reliability of the link $e$. We
denote by $\mathcal{R}_{G}(p)$ the _all-terminal reliability_ of graph $G$:
the probability that $G$ is connected.

Given two graphs $G_{1}$ and $G_{2}$ with two distinguished vertices
$s(G_{i})$ and $t(G_{i})$, a _series composition_ of $G_{1}$ and $G_{2}$,
denoted by $G_{1}+_{S}G_{2}$, is the disjoint union of both graphs, merging
$t(G_{1})$ with $s(G_{2})$. In this case, $s(G_{1}+_{S}G_{2})=s(G_{1})$ and
$t(G_{1}+_{S}G_{2})=t(G_{2})$. Similarly, a _parallel composition_
$G_{1}+_{P}G_{2}$ is the disjoint union of both graphs, merging $s(G_{1})$
with $s(G_{2})$ (thus $s(G_{1})=s(G_{2})=s(G_{1}+_{S}G_{2})$ and merging
$t(G_{1})$ with $t(G_{2})$ (thus $t(G_{1})=t(G_{2})=t(G_{1}+_{P}G_{2})$).

A graph $G$ is a _series-parallel graph_ if it can be obtained from a sequence
of series-parallel compositions starting from the single edges $G_{e}=\\{e\\}$
for $e=1\ldots|E|$. Formally, let $\mathcal{G}_{0}=\\{G_{e}:e=1\ldots|E|\\}$
be the set of single edges. Iteratively, we construct
$\mathcal{G}_{i}=\mathcal{G}_{i-1}\bigcup(G_{|E|+i}\setminus(G_{j}\cup
G_{k}))$ where $G_{|E|+i}=G_{j}+_{\odot}G_{k}$ is either a series or parallel
composition of graphs $G_{j}$ and $G_{k}$ in $\mathcal{G}_{i-1}$. Note that
the number of connected components in $\mathcal{G}_{i}$ is $|E|-i$ because
this number decreases by one in each iteration. Therefore,
$\mathcal{G}_{|E|-1}$ contains only one element, which is a connected series-
parallel graph. We denote by
$\mathcal{S}_{G}=\left[G_{|E|+i}=G_{j}+_{\odot}G_{k}\right]_{i=1}^{|E|-1}$
this sequence of series-parallel compositions to construct $G$. Note that this
sequence is not unique for a given graph $G$.

Satyanarayana and Wood [29] presented a set of reliability-preserving
transformations for computing the reliability of series-parallel graphs based
on a sequence of series-parallel compositions.These transformations based on
series-parallel compositions can also be used in general graphs to reduce
their size in their reliability computation.

When two edges $e_{j},e_{k}$ are _in parallel_ , with reliabilities
$p_{j},p_{k}$, these edges can be replaced by a new single edge $e_{i}$ with
reliability $p_{i}=1-(1-p_{j})(1-p_{k})$, which is the probability that at
most one of these two links fails.

If two edges $e_{j},e_{k}$ are _in series_ , with reliabilities $p_{j},p_{k}$,
at least one of them must remain operational to keep the graph connected.
Thus, if we replace these two edges with a new edge $e_{i}$, the reliability
of this edge must consider this event. Let $\mathcal{A}$ be the event that
$e_{j}$ and $e_{j}$ do not fail simultaneously. Hence, the reliability of the
graph satisfies

$\displaystyle\mathcal{R}_{G}$ $\displaystyle=\mathbb{P}[G\text{ is
connected}|\mathcal{A}]\cdot\mathbb{P}[\mathcal{A}]+\underbrace{\mathbb{P}[G\text{
is connected}|\mathcal{A}^{C}]}_{0}\cdot\mathbb{P}[\mathcal{A}^{C}]$
$\displaystyle=\mathbb{P}[G\text{ is
connected}|\mathcal{A}]\cdot\mathbb{P}[\mathcal{A}]$

Therefore, in the case of a series reduction replacing edges $e_{j}$ and
$e_{k}$ with a new edge $e_{i}$, the reliability of the resulting graph must
be multiplied by $\mathbb{P}[\mathcal{A}]=1-(1-p_{j})(1-p_{k})$, and the
reliability of the new edge $e_{i}$ is the probability that both edges are
operational conditional to the event that at least one of them remains
operational, that is $p_{j}\cdot p_{k}$ normalized by the probability of
$\mathcal{A}$:

$p_{i}=\frac{p_{j}p_{k}}{1-(1-p_{j})(1-p_{k})}.$

Finally, note that the sequence $\mathcal{S}_{G}$ that _constructs_ the graph
$G$ can be used to compute the reliability of $G$. That is, if
$G_{|E|+i}=G_{j}+_{\odot}G_{k}$ and $G_{j}$ and $G_{k}$ corresponds to an edge
of $E$, then $G_{|E|+i}$ is a new edge that replaces the original two edges
with a new edge representing its series/parallel composition. Therefore,
applying this sequence iteratively, at each step of the sequence
$\mathcal{S}_{G}$, the subgraphs $G_{j},G_{k}$ in a composition are edges.

We formalize the reliability computation of a series-parallel graph $G$ in
Algorithm 1.

Algorithm 1 Compute the reliability of a series-parallel graph $G$ from its
composition sequence $\mathcal{S}_{G}$.

Composition sequence
$\mathcal{S}_{G}=\left[G_{|E|+i}=G_{j}+_{\odot}G_{k}\right]_{i=1}^{|E|-1}$

$p_{e}$ reliabilities of edges $e\in 1\ldots|E|$

for $i=1\ldots|E|$ do

$G_{i}\leftarrow\\{i\\}$

$Y_{i}\leftarrow p_{i}$

$\Omega_{i}\leftarrow 1$

end for

for $i=1\ldots|E|-1$ do

if $G_{|E|+i}=G_{j}+_{P}G_{k}$ then

$Y_{|E|+i}\leftarrow 1-(1-Y_{j})(1-Y_{k})$

$\Omega_{|E|+i}\leftarrow 1$

else if $G_{|E|+i}=G_{j}+_{S}G_{k}$ then

$Y_{|E|+i}\leftarrow\frac{Y_{j}Y_{k}}{1-(1-Y_{j})(1-Y_{k})}$

$\Omega_{|E|+i}\leftarrow 1-(1-Y_{j})(1-Y_{k})$

end if

end for

return $\mathcal{R}_{G}:=Y_{2|E|-1}\cdot\prod_{i=1}^{2|E|-1}\Omega_{i}$

In other words, $Y_{i}$ represents the reliability of the edge $i$ for
$i\leq|E|$, or the reliability of the edge resulting from the reduction
$G_{|E|+i}=G_{j}+_{\odot}G_{k}$ for $i>|E|$. Similarly, $\Omega_{|E|+i}$
represents the reliability factor from the reduction $G_{|E|+i}$, which is
either $\Omega_{|E|+i}=1$ (parallel composition) or
$\Omega_{|E|+i}=1-(1-p_{j})(1-p_{k})$ (series composition). When the graph has
been reduced to a single edge, the reliability of $G$ is equal to the
reliability of this edge ($Y_{2|E|-1}$) multiplied by all the factors
$\Omega_{i}$. Algorithm 1 allows one to compute the all-terminal reliability
of $G$ in linear time.

##  3\. Optimizing the reliability of a series-parallel graph

###  3.1. A nonlinear optimization model for a series-parallel graph

The aforementioned results provide a procedure to compute the resulting
reliability $\mathcal{R}_{G}$ of a series-parallel graph $G$. Our main
interest is in studying the network design problem of selecting the subgraph
that maximizes reliability given a set of constraints. Specifically, given a
graph $G=(V,E)$, we are interested in the selection of a subset of edges
$F\subseteq E$ satisfying the given constraints such that the reliability of
the graph $(V,F)$ is maximized.

We now proceed to formulate this problem as an MINLP. Let $X_{e}\in\\{0,1\\}$
for $e\in E$ be binary variables indicating whether $e\in F$, and let $AX\leq
b$ be a given set of the arbitrary linear constraints that any valid $X$ must
satisfy. These can be, for example, an upper bound on the number of edges to
be considered. Let $\mathcal{R}(p)$ be the reliability of network $G$ given
the probability vector $p\in[0,1]^{E}$—each component $p_{e}$ is the
reliability of $e\in E$. We note that if a link $e\in E$ is not considered in
$F$, this is equivalent to assuming that its elementary reliability is 0.
Then, we can formulate our problem as:

$\displaystyle\max$ $\displaystyle\
\mathcal{R}(p_{1}X_{1},p_{2}X_{2},\ldots,p_{|E|}X_{|E|})$ $\displaystyle
s.t.\,\,$ $\displaystyle AX\leq b$ $\displaystyle
X_{e}\in\\{0,1\\}\quad\forall e\in 1\ldots|E|$

If $G$ is a series-parallel graph, we can apply the reliability-preserving
reductions over the composition sequence
$\mathcal{S}_{G}=\left[G_{|E|+i}=G_{j}+_{\odot}G_{k}\right]_{i=1}^{|E|-1}$
that constructs $G$. Following the idea and notation behind Algorithm 1, we
define the continuous variables $Y_{i}\in[0,1]$ and $\Omega_{i}\in[0,1]$ for
each $i=1\ldots 2|E|-1$ to represent the reliability and the correction factor
of each step of the sequence. Additionally, we define the continuous variables
$\bar{\Omega}_{i}\in[0,1]$ to represent the product of the correction factors.
Using these variables, we can formulate the problem as follows:

(1a) $\displaystyle\max\ $ $\displaystyle R$ (1b) $\displaystyle AX$
$\displaystyle\leq b$ (1c) $\displaystyle Y_{i}$ $\displaystyle=p_{i}X_{i}$
$\displaystyle i=1\ldots|E|$ (1d) $\displaystyle\Omega_{i}$ $\displaystyle=1$
$\displaystyle i=1\ldots|E|$ (1e) $\displaystyle Y_{|E|+i}$
$\displaystyle=1-(1-Y_{j})(1-Y_{k})$ $\displaystyle
i=1\ldots|E|-1:G_{|E|+i}=G_{j}+_{P}G_{k}$ (1f) $\displaystyle\Omega_{|E|+i}$
$\displaystyle=1$ $\displaystyle i=1\ldots|E|-1:G_{|E|+i}=G_{j}+_{P}G_{k}$
(1g) $\displaystyle Y_{|E|+i}$
$\displaystyle=\frac{Y_{j}Y_{k}}{1-(1-Y_{j})(1-Y_{k})}$ $\displaystyle
i=1\ldots|E|-1:G_{|E|+i}=G_{j}+_{S}G_{k}$ (1h) $\displaystyle\Omega_{|E|+i}$
$\displaystyle=1-(1-Y_{j})(1-Y_{k})$ $\displaystyle
i=1\ldots|E|-1:G_{|E|+i}=G_{j}+_{S}G_{k}$ (1i) $\displaystyle\bar{\Omega}_{1}$
$\displaystyle=\Omega_{1}$ $\displaystyle i=2\ldots 2|E|-1$ (1j)
$\displaystyle\bar{\Omega}_{i}$
$\displaystyle=\bar{\Omega}_{i-1}\cdot\Omega_{i}$ $\displaystyle i=2\ldots
2|E|-1$ (1k) $\displaystyle R$
$\displaystyle=Y_{2|E|-1}\cdot\bar{\Omega}_{2|E|-1}$ (1l) $\displaystyle
Y_{i},\Omega_{i},\bar{\Omega}_{i}$ $\displaystyle\in[0,1]$ $\displaystyle
i=1\ldots 2|E|-1$ (1m) $\displaystyle X_{i}$ $\displaystyle\in\\{0,1\\}$
$\displaystyle i=1\ldots|E|$

Constraints (1c)-(1d) correspond to the main decision variables, indicating
whether an edge $e$ is considered in the subgraph $F$. In the latter case,
edge $e$ has elementary reliability equal to zero. Constraints (1e)-(1f) model
a parallel reduction $G_{|E|+i}=G_{j}+_{P}G_{k}$, in which case the new edge
$i$ has reliability $1-(1-p_{i})(1-p_{k})$ and there is no reliability
correction factor ($\Omega_{|E|+i}=1$). Constraints (1g)-(1h) model a series
reduction, where the new edge has reliability
$\frac{p_{j}p_{k}}{1-(1-p_{j})(1-p_{k})}$ and the reliability correction
factor is $\Omega_{|E|+i}=1-(1-p_{j})(1-p_{k})$. Constraints (1i)-(1j) ensure
that $\bar{\Omega}_{i}$ is equal to the cumulative product of the factors
$\Omega_{i}$, that is, $\bar{\Omega}_{i}=\prod_{j<i}\Omega_{j}$. Finally,
constraint (1k) provides the reliability of the graph $R$, which is the
operational probability after the last reduction, i.e., when $G$ has been
reduced to a single edge.

Note that this problem can be seen as a _mixed-integer quadratically
constrained program_ (MIQCP); the left side of constraint (1g) can also be
written as $Y_{i}\cdot\Omega_{i}=Y_{j}\cdot Y_{k}$. However, all quadratic
constraints are nonconvex, and thus this model can be challenging for most of
the current nonlinear optimization solvers.

###  3.2. Convex envelopes for the problem

The main issue with model (1) is that the resulting constraints involve
nonconvex and nonconcave functions. In fact, the three bivariate functions
$f_{1}(x,y)=xy$, $f_{2}(x,y)=1-(1-x)(1-y)$ and
$f_{3}(x,y)=\tfrac{xy}{1-(1-x)(1-y)}$ appearing in (1) are neither convex nor
concave functions111Without loss of generality, we assume that
$f_{3}(0,0)=\lim_{(x,y)\rightarrow(0,0)}f_{3}(x,y)=0$..

One common approach for generating a (possibly strong) convex relaxation is to
use the _envelopes_ of these functions. The concave envelope of $f(x,y)$ over
a given domain $D$ is the smallest concave overestimator ${f}^{ave}(x,y)\geq
f(x,y)$ for all $(x,y)\in D$ and can be used to relax a constraint of type
$f(x,y)\geq z$ with a convex constraint ${f}^{ave}(x,y)\geq z$. Similarly, the
convex envelope of $f(x,y)$ is the largest convex underestimator
${f}^{vex}(x,y)\leq f(x,y)$ for all $(x,y)\in D$ and can be used to relax a
constraint $f(x,y)\leq z$. In our setting, this implies that we can relax an
equality constraint $z=f_{i}(x,y)$ with two convex constraints
${f}^{vex}_{i}(x,y)\leq z\leq{f}^{ave}_{i}(x,y)$.

Since our optimization problem only considers equality constraints, in
principle, we should aim at computing both convex and concave envelopes.
However, due to the structure of our problem, constraints
${f}^{vex}_{i}(x,y)\leq z$ are not necessary. We show below that
${f}^{ave}_{i}(x,y)$ are all increasing functions in both variables, and
considering that we are maximizing reliability, along with the simple
structure of our constraints, $z\leq{f}^{ave}_{i}(x,y)$ will always be active
in an optimal solution of the resulting convex relaxation. This is expected,
as $f_{1}(x,y),f_{2}(x,y)$ and $f_{3}(x,y)$ are increasing functions in both
variables in the square $[0,1]\times[0,1]$.

For the case of $f_{1}(x,y)=x\cdot y$, its envelopes are well-known and can be
obtained by the McCormick envelopes [21]. Let us assume that $x,y\in[0,1]$,
and let $L_{x},L_{y}$ and $U_{x},U_{y}$ be lower and upper bounds for $x$ and
$y$; then, the concave envelope of $f_{1}(x,y)$ is given by:

(2) $f_{1}(x,y)=x\cdot y\leq\begin{cases}U_{x}\cdot y+x\cdot L_{y}-U_{x}\cdot
L_{y}&\text{if }y-L_{y}\geq\frac{U_{y}-L_{y}}{U_{x}-L_{x}}\cdot(x-L_{x})\\\
x\cdot U_{y}+L_{x}\cdot y-L_{x}\cdot U_{y}&\text{if not}\end{cases}$

This envelope is increasing in both variables when the variable bounds are
nonnegative.

For the case of $f_{2}(x,y)=1-(1-x)(1-y)=x+y-x\cdot y$, its concave envelope
can be obtained using the convex envelope of $f_{1}(x,y)$, resulting in the
following piecewise linear function:

(3)
$f_{2}(x,y)=1-(1-x)(1-y)\leq\begin{cases}x\cdot(1-L_{y})+y\cdot(1-L_{x})+L_{x}L_{y}&\text{if
}U_{y}-y\geq\frac{U_{y}-L_{y}}{U_{x}-L_{x}}\cdot(x-L_{x})\\\
x\cdot(1-U_{y})+y\cdot(1-U_{x})+U_{x}U_{y}&\text{if not}\end{cases}$

Since in our case both variable bounds are less than 1, this envelope is also
increasing in both variables.

For the case of $f_{3}(x,y)=\tfrac{x\cdot y}{x+y-xy}$, an explicit formula for
its concave envelopes is not known. Here, we provide its concave envelope for
the case of $L_{x}=L_{y}=0$.

######  Theorem 1.

The concave envelope of $f_{3}(x,y)=\tfrac{x\cdot y}{x+y-xy}$ for
$(x,y)\in[0,U_{x}]\times[0,U_{y}]$ for $U_{x},U_{y}\leq 1$ is given by:

(4) ${f}^{ave}_{3}(x,y)=\begin{cases}\frac{x\cdot y}{x+y-U_{x}\cdot
y}&\text{if }x/U_{x}\geq y/U_{y}\\\ \frac{x\cdot y}{x+y-x\cdot U_{y}}&\text{if
not}\end{cases}$

###### Proof.

It is easy to see that ${f}^{ave}_{3}(x,y)\geq f_{3}(x,y)$. On the other hand,
note that for $(x,y)\in[0,U_{x}]\times[0,U_{y}]$,

${f}^{ave}_{3}(x,y)=\min\left\\{\frac{x\cdot y}{x+y-U_{x}y},\frac{x\cdot
y}{x+y-xU_{y}}\right\\}$

The Hessian matrix of $\frac{x\cdot y}{x+y-U_{x}y}$ is:

$\mathbf{H}\left(\frac{x\cdot y}{x+y-U_{x}\cdot
y}\right)=\frac{2(1-U_{x})}{(x+y-U_{x}\cdot y)^{3}}\begin{bmatrix}-y^{2}&xy\\\
xy&-x^{2}\end{bmatrix}$

which is a negative semidefinite-matrix, and thus it is a concave function in
$[0,U_{x}]\times[0,U_{y}]$. The same result can be obtained for $\frac{x\cdot
y}{x+y-x\cdot U_{y}}$ by exchanging $x$ and $y$. Finally, since
${f}^{ave}_{3}(x,y)$ is the minimum of these two concave functions, we
conclude that ${f}^{ave}_{3}(x,y)$ is concave over $[0,U_{x}]\times[0,U_{y}]$.

Finally, we need to show that ${f}^{ave}_{3}(x,y)$ is the _smallest_ concave
overestimator of $f_{3}(x,y)$. We note that if $y=\lambda x$, then the
function $f_{3}(x,\lambda x)=\tfrac{\lambda x}{1+\lambda-\lambda x}$ is a
convex function on $x$. In fact, $\tfrac{\partial f_{3}(x,\lambda x)}{\partial
x}=2\lambda^{2}(1+\lambda)/(1+\lambda-\lambda x)^{3}$, which is positive for
any $\lambda>0$. Therefore, the best possible overestimator over the line
$y=\lambda x$ is given by the linear function interpolating the origin and the
intersection of $(x,\lambda x)$ with either $x=U_{x}$ or $y=U_{y}$. We show
that ${f}^{ave}_{3}(x,y)$ is a function satisfying this condition. In fact, if
$U_{y}/U_{x}\geq\lambda$, then $y=\lambda x$ intersects first with $x=U_{x}$
and ${f}^{ave}_{3}(x,\lambda x)=\tfrac{\lambda x}{1+\lambda-\lambda U_{x}}$,
which is a linear function, and $f_{3}(U_{x},\lambda
U_{x})={f}^{ave}_{3}(U_{x},\lambda U_{x})$. Otherwise, if
$U_{y}/U_{x}\leq\lambda$, then $y=\lambda x$ intersects first with $y=U_{y}$
and ${f}^{ave}_{3}(x,\lambda x)=\tfrac{\lambda x}{1+\lambda-U_{y}}$, which is
a linear function such that
$f_{3}(U_{y}/\lambda,U_{y})={f}^{ave}_{3}(U_{y}/\lambda,U_{y})$. ∎∎

This idea of exploiting the convexity of $f_{3}$ over the rays
$f_{3}(x,\lambda x)$ can also be extended to find concave envelopes for other
functions satisfying this property; this is further elaborated in parallel
work [5]. Note that, as anticipated, ${f}^{ave}_{3}(x,y)$ is an increasing
function in both variables in $[0,U_{x}]\times[0,U_{y}]$; it can be easily
verified that

$\nabla\left(\frac{x\cdot y}{x+y-U_{x}\cdot
y}\right)=\frac{1}{(-U_{x}y+x+y)^{2}}\left((1-U_{x})y^{2},x^{2}\right)$

which is a nonnegative vector whenever $U_{x}\leq 1$. The other part of the
definition of ${f}^{ave}_{3}(x,y)$ can be verified similarly.

In Figure 1, we show the three concave envelopes we have discussed in this
section. Using these envelopes, we can formulate a mixed-integer convex
nonlinear problem that can be solved more efficiently than the original model.
This provides a tractable overestimation of the reliability of the resulting
graph.

Figure 1. Concave envelopes for $x\cdot y$, $x+y-x\cdot y$ and
$\tfrac{xy}{x+y-x\cdot y}$ for $0\leq x\leq 1$ and $0\leq y\leq 1$.

###  3.3. A mixed integer convex approximation

Using the concave envelopes of the bivariate functions resulting from the
series and parallel reductions, we can replace the corresponding constraint
from model (1a)-(1m) and obtain the following mixed-integer convex
optimization approximation model:

(5a) $\displaystyle\max\ R$ (5b) $\displaystyle AX\leq b$ (5c) $\displaystyle
Y_{i}$ $\displaystyle=p_{i}X_{i},\quad\Omega_{i}=1,\quad\bar{\Omega}_{i}=1$
$\displaystyle i=1\ldots|E|$ (5d) $\displaystyle Y_{|E|+i}$ $\displaystyle\leq
Y_{j}(1-U_{k})+(1-U_{j})Y_{k}+U_{j}U_{k}$ $\displaystyle
i=1\ldots|E|-1:G_{i}=G_{j}+_{P}G_{k}$ (5e) $\displaystyle Y_{|E|+i}$
$\displaystyle\leq Y_{j}(1-L_{k})+(1-L_{j})Y_{k}+L_{j}L_{k}$ $\displaystyle
i=1\ldots|E|-1:G_{i}=G_{j}+_{P}G_{k}$ (5f) $\displaystyle\Omega_{|E|+i}$
$\displaystyle=1,\quad\bar{\Omega}_{|E|+i}=\bar{\Omega}_{|E|+i-1}$
$\displaystyle i=1\ldots|E|-1:G_{i}=G_{j}+_{P}G_{k}$ (5g) $\displaystyle
Y_{|E|+i}$
$\displaystyle\leq\begin{cases}\frac{Y_{j}Y_{k}}{Y_{j}+Y_{k}-U_{j}Y_{k}}&\text{if
}x/U_{x}\geq y/U_{y}\\\ \frac{Y_{j}Y_{k}}{Y_{j}+Y_{k}-Y_{j}U_{k}}&\text{if
not}\end{cases}$ $\displaystyle i=1\ldots|E|-1:G_{i}=G_{j}+_{S}G_{k}$ (5h)
$\displaystyle\Omega_{|E|+i}$ $\displaystyle\leq
Y_{j}(1-U_{k})+(1-U_{j})Y_{k}+U_{j}U_{k}$ $\displaystyle
i=1\ldots|E|-1:G_{i}=G_{j}+_{S}G_{k}$ (5i) $\displaystyle\Omega_{|E|+i}$
$\displaystyle\leq Y_{j}(1-L_{k})+(1-L_{j})Y_{k}+L_{j}L_{k}$ $\displaystyle
i=1\ldots|E|-1:G_{i}=G_{j}+_{S}G_{k}$ (5j) $\displaystyle\bar{\Omega}_{|E|+i}$
$\displaystyle\leq
U^{\bar{\Omega}}_{|E|+i-1}\cdot\Omega_{|E|+i}+\bar{\Omega}_{|E|+i-1}\cdot
L^{\Omega}_{|E|+i}-U^{\bar{\Omega}}_{i-1}\cdot L^{\Omega}_{i}$ $\displaystyle
i=1\ldots|E|-1:G_{i}=G_{j}+_{S}G_{k}$ (5k) $\displaystyle\bar{\Omega}_{|E|+i}$
$\displaystyle\leq
L^{\bar{\Omega}}_{i-1}\cdot\Omega_{i}+\bar{\Omega}_{i-1}\cdot
U^{\Omega}_{i}-L^{\bar{\Omega}}_{i-1}\cdot U^{\Omega}_{i}$ $\displaystyle
i=1\ldots|E|-1:G_{i}=G_{j}+_{S}G_{k}$ (5l) $\displaystyle R$
$\displaystyle\leq U_{2|E|-1}\cdot\bar{\Omega}_{2|E|-1}+Y_{z}\cdot
L^{\bar{\Omega}}_{2|E|-1}-U_{2|E|-1}\cdot L^{\bar{\Omega}}_{2|E|-1}$ (5m)
$\displaystyle R$ $\displaystyle\leq
L_{2|E|-1}\cdot\bar{\Omega}_{2|E|-1}+Y_{z}\cdot
U^{\bar{\Omega}}_{2|E|-1}-L_{2|E|-1}\cdot U^{\bar{\Omega}}_{2|E|-1}$ (5n)
$\displaystyle Y_{i}$ $\displaystyle,\Omega_{i},\bar{\Omega}_{i}\in[0,1]$
$\displaystyle i=1\ldots 2|E|-1$ (5o) $\displaystyle X_{i}$
$\displaystyle\in\\{0,1\\}$ $\displaystyle i=1\ldots|E|$

where the constants $L_{i}$ and $U_{i}$ are valid lower and upper bounds for
$Y_{i}$ and $L^{\Omega}_{i}$, $U^{\Omega}_{i}$, $L^{\bar{\Omega}}_{i}$ and
$U^{\bar{\Omega}}_{i}$ are valid lower and upper bounds for variables
$\Omega_{i}$ and $\bar{\Omega}_{i}$, respectively. These upper bounds can be
precomputed by assigning $U_{i}=p_{i}$ and $L_{i}=0$ for all $i\in 1\ldots|E|$
and then applying the corresponding functions $f_{1}$, $f_{2}$ or $f_{3}$ to
these bounds for each series or parallel composition in $\mathcal{S}_{G}$.

All constraints in previous model are linear, except for inequality (5g).
However, since ${f}^{ave}_{3}(x,y)$ is concave, we can enforce this nonlinear
constraint with linear constraints given by its tangent hyperplane. Given a
point $(x^{*},y^{*})$, we upper bound ${f}^{ave}_{3}(x,y)$ by the linear
constraint ${f}^{ave}_{3}(x^{*},y^{*})+\frac{\partial{f}^{ave}_{3}}{\partial
x}(x^{*},y^{*})\cdot(x-x^{*})+\frac{\partial{f}^{ave}_{3}}{\partial
y}(x^{*},y^{*})\cdot(y-y^{*})$, which is

(6) $\displaystyle{f}^{ave}_{3}(x,y)$
$\displaystyle\leq\begin{cases}\left(\frac{y^{*}}{x^{*}+y^{*}-U_{x}\cdot
y^{*}}\right)^{2}\cdot(1-U_{x})\cdot
x+\left(\frac{x^{*}}{x^{*}+y^{*}-U_{x}y^{*}}\right)^{2}\cdot y&\text{if
}x/U_{x}\geq y/U_{y}\\\
\left(\frac{y^{*}}{x^{*}+y^{*}-x^{*}U_{y}}\right)^{2}\cdot
x+\left(\frac{x^{*}}{x^{*}+y^{*}-x^{*}U_{y}}\right)^{2}\cdot(1-U_{y})\cdot
y&\text{if not}\end{cases}$

These linear constraints can be added dynamically to optimization solvers
during their optimization procedures.

###  3.4. Computational experiments

####  3.4.1. Instances and implementation

To test the effectiveness of our methodology, we generated random series-
parallel graphs with $|E|=50$ and $|E|=100$ edges in the following way. The
elementary reliabilities $p_{e}$ are generated uniformly at random between
$0.9$ and $1.0$ for each edge $e$. To generate each instance, we start with a
perfect bipartite matching. Clearly, the number of connected components in $G$
is initially $|E|$. Then, we iteratively select two components uniformly at
random and connect them either in series or parallel (with equal probability),
which diminishes the number of connected components in $G$ by one. The process
is repeated until the resulting graph is connected. This procedure is detailed
in Algorithm 2.

For the additional constraints (Eq. 1b), we impose a cardinality constraint
that only a $\alpha$ fraction of the edges can be selected: $\sum_{e\in
E}X_{e}\leq\alpha|E|$. Low values of $\alpha$ ($\alpha<0.5$) lead to either
infeasible problems or a reduced set of feasible solutions, which makes the
problem easy to solve. Similarly, high values of $\alpha$ ($\alpha>0.8$)
encourage the solution to include most of the edges, also making the problem
easy to solve. For these reasons, in our experiments, we use $\alpha=0.8$,
which is a high value where our models behave very well, and $\alpha=0.6$,
which yields the most challenging instances of our problem.

In our computational experiments, we compare two different configurations:

  1. Convex envelope cuts:

The model described in Section 3.3, including the cuts from (6) to approximate
our concave envelope from Theorem 1.

  2. Without cuts:

The model from Section 3.3, considering a simpler constraint to bound the
nonlinear concave constraint (5g), which we describe next.

The second configuration is constructed to understand the effectiveness of the
concave approximation provided in Theorem 1 for series compositions. This
simpler general approximation for $f_{3}(x,y)$ is given by constraints
$f_{3}(x,y)\leq x$ and $f_{3}(x,y)\leq y$; these constraints are obtained when
replacing either $x$ or $y$ in the denominator of $f_{3}(x,y)$ with 1 or by
replacing $U_{x}$ or $U_{y}$ with 1 in ${f}^{ave}_{3}(x,y)$. This corresponds
to the tangent hyperplanes of ${f}^{ave}_{3}(x,y)$ on $(x,y)=(0,U_{y})$ and
$(x,y)=(U_{x},0)$.

Algorithm 2 Generation of random instances of size $k$

$E\leftarrow\\{e_{1},\ldots e_{k}\\}$ (disjoint edges)

$p_{e}\leftarrow U[0.9,1]$ for all $e\in E$ $\triangleright$ Random initial
reliabilities

$\mathcal{G}\leftarrow E$

while $|ConnectedComponents(\mathcal{G})|>1$ do

$g_{1},g_{2}\leftarrow$ random subgraphs from $\mathcal{G}$ (chosen
equiprobable)

$u\leftarrow U[0,1]$

if $u<0.5$ then

$g=g_{1}+_{P}g_{2}$ $\triangleright$ Parallel connection

else

$g=g_{1}+_{S}g_{2}$ $\triangleright$ Series connection

end if

$\mathcal{G}\leftarrow\mathcal{G}\setminus\\{g_{1}\cup g_{2}\\}\cup g$
$\triangleright$ Replace $g_{1}$ and $g_{2}$ in $\mathcal{G}$ by $g$

end while

These models were implemented with Python 3.7 using the IBM® Decision
Optimization CPLEX® Modeling for Python (DOcplex.MP) v2.11 of CPLEX Studio
v12.10 [18]. All CPLEX parameters have their default values, and no cut
manager was implemented for the additional cuts (6). All computations were
made on machines running Linux under x86_64 architecture running in a single
thread.

####  3.4.2. Computational results

Figure 2. Comparison between the reported objective value and the true
reliability and the optimal reliability of the resulting network

Figure 2 (left) compares the objective value reported by our model (x-axis)
for each instance versus the true reliability of the resulting solution for
the two configurations described. Similarly, Figure 2 (right) compares the
objective value reported versus the true optimal solution of the problem, for
the cases where the latter is known.222True optimal solutions are obtained
using the results of the following sections but presented here for
illustrative purposes. Each subfigure is divided in four, depending on the
number of edges (horizontal) and the budget (vertical) of the problem.

In these experiments, we can see that both configurations behave extremely
well for a budget of $\alpha=0.8$: the resulting objective value of the model
is not only very close to the real reliability of the solution but also very
close to the true optimal solution of the problem. These relative errors are
less than $1\%$ for all instances and configurations evaluated.

For $\alpha=0.6$, this behavior changes drastically. The reported objective
value of the problem (an upper bound on the reliability) differs considerably
for many instances in both configurations. Moreover, it seems that this effect
is more pronounced for lower reliabilities. A similar effect is observed when
we compare the reported reliability with the optimal reliability of the
problem.

These experiments show that, in challenging instances, the overestimation of
the reliability provided by the concave envelopes alone is not good enough for
estimating the true reliability of the resulting graphs, leading to suboptimal
solutions for our problem. Moreover, they indicate that little is gained from
using cutting planes to iteratively approximate ${f}^{ave}_{3}$, in comparison
to its simple approximation based on two tangent hyperplanes, at least in this
setting. In the next section, we study different improvements to avoid these
issues and where the cutting planes associated with ${f}^{ave}_{3}$ have
considerably more impact. These improvements lead the model to the true
optimal solution in most instances.

###  3.5. Further improvements to the model

####  3.5.1. Providing the true reliability of the solution

The previous model considers overestimators of the nonlinear functions, which
can lead to overestimations of the true reliability of the resulting network.
However, based on the selected edges of the graph (the values of the $X$
variables), we can compute the “true” reliability of the resulting graph, and
we can improve our model by introducing _combinatorial Benders cuts_ [11] to
avoid this problem.

Specifically, given a feasible solution $\hat{X}$ we can compute the resulting
reliability
$\mathcal{R}_{\hat{X}}=\mathcal{R}(p_{1}\hat{X}_{1},\ldots,p_{|E|}\hat{X}_{|E|})$
in linear time using Algorithm 1 and obtain the following valid inequality—a
combinatorial Benders cut associated with $\hat{X}$—that bounds the value of
$R$

(7) $R\leq\mathcal{R}_{\hat{X}}+\sum_{e\in
E:\hat{X}_{e}=1}(1-X_{e})+\sum_{e\in E:\hat{X}_{e}=0}X_{e}$

This constraint implies that if $\hat{X}$ is the optimal solution to the
problem, then $R=\mathcal{R}_{\hat{X}}$. For other feasible solutions
$X^{*}\neq\hat{X}$, this constraint is trivially satisfied because $R\leq 1$.
We can strengthen this inequality using the following observation. Note that
the all-terminal reliability of a graph
$\mathcal{R}(p_{1}X_{1},p_{2}X_{2},\ldots,p_{|E|}X_{|E|})$ only increases when
a new edge is added. In other words, the reliability of any solution $X$ such
that $\hat{X}_{e}=0\Rightarrow X_{e}=0$ must be smaller than that of
$\hat{X}$. Therefore, the combinatorial cuts can be strengthened to

(8) $R\leq\mathcal{R}_{\hat{X}}+\sum_{e\in E:\hat{X}_{e}=0}X_{e}.$

By the same reasoning, if a feasible solution of the problem contains all
edges selected in $\hat{X}$, then its reliability will be at least
$\mathcal{R}_{\hat{X}}$. Hence, we can enforce this lower bound for the
reliability by adding the following constraint each time that an incumbent
solution $\hat{X}$ has been found:

(9) $R\geq\mathcal{R}_{\hat{X}}-\sum_{e\in E:\hat{X}_{e}=1}(1-X_{e}).$

This reasoning can also be expanded to other variables. The monotonicity
exhibited by the reliability function also holds for $Y_{i}$, $\Omega_{i}$ and
$\bar{\Omega}_{i}$ as functions of $X$, because all three functions
$f_{1}(x,y)$, $f_{2}(x,y)$ and $f_{3}(x,y)$ are nondecreasing in both
dimensions. Therefore, given a solution $\hat{X}$, we can apply Algorithm 1 to
compute the values of variables $Y$, $\Omega$ and $\bar{\Omega}$ associated
with this solution and derive similar cuts for all variables of the problem
each time that a new incumbent solution $\hat{X}$ is found during the branch-
and-bound process.

######  Remark 1.

Including these inequalities during a branch-and-bound procedure ensures that
we obtain an optimal solution to the original problem (1) at the end of the
optimization routine, unless a time limit is reached. Whenever a feasible
solution $\hat{X}$ to the relaxation is found, the combinatorial Benders cut
ensures that $R=\mathcal{R}_{\hat{X}}$ if $\hat{X}$ is reported optimal. The
resulting optimization routine might be impractical though, as it might resort
to a costly enumeration if the relaxation (5) is not tight enough. The next
improvements aim at better approximating the problem and making the tree-
search more efficient.

####  3.5.2. Improving inequalities on the branch-and-bound tree

Our proposed model considers the best possible concave overestimators for each
function. These envelopes depend on the lower and upper bounds for each
variable, and even if these bounds are tight, the envelopes may not provide a
tight approximation of the functions in the whole feasible region.
Nevertheless, the branch-and-bound procedure of MILP solvers is based on
imposing new bounds on the variables while branching and thus improving the
relaxations in each node. These bounds are only valid locally, but we can use
them to obtain the concave envelopes based on these new bounds, which yields
better local approximations of each function.

If the branch-and-bound process fixes a variable $X_{e}$ to 0 or 1, we can
propagate this decision to improve the lower and upper bounds for all
variables corresponding to reductions that include this edge. For instance, if
the fixed variable is $X_{j}=1$, then $Y_{j}=p_{j}$. Thus, for the case of a
parallel reduction, we can improve (5d) and (5e) by adding the local linear
constraint $Y_{i}=1-(1-p_{j})(1-Y_{k})$. In the case of a series reduction, we
can improve (5g) and its associated cuts (6) by adding the inequality

(10) $Y_{i}\leq\left(\frac{p_{j}}{Y_{k}^{*}+p_{j}-Y^{*}_{k}\cdot
p_{j}}\right)^{2}\cdot(Y_{k}-Y_{k}^{*})+\frac{p_{j}Y_{k}^{*}}{Y_{k}^{*}+p_{j}-Y^{*}_{k}\cdot
p_{j}}$

where $Y_{k}^{*}$ is the current solution of variable $Y_{k}$ at the node of
the branch-and-bound tree. This is possible because the function $f_{3}(p,y)$
is concave on $y$ for any fixed value of $p$, so (10) corresponds to the
gradient of this function on $y=Y_{k}^{*}$. Similar improvements can be
included for the remaining equations involving $\Omega_{i}$ and
$\bar{\Omega}_{i}$.

Note that our concave envelopes for series reduction (Theorem 1) only apply
when the lower bounds are equal to $0$, so this approximation cannot be
improved if the lower bounds are improved. However, we can still add linear
constraints in this case: since $Y_{i}=f_{3}(Y_{j},Y_{k})$ is increasing and
concave for a fixed $Y_{j}$ or $Y_{k}$, we can overestimate this function by
two hyperplanes tangent to the point $(L_{j},L_{k})$.

$\displaystyle Y_{i}$ $\displaystyle\leq\left.\frac{\partial f_{3}}{\partial
x}\right|_{(L_{j},L_{k})}\cdot(Y_{j}-L_{j})+\left.\frac{\partial
f_{3}}{\partial
y}\right|_{(U_{j},L_{k})}\cdot(Y_{k}-L_{k})+f_{3}(L_{j},L_{k})$
$\displaystyle=\frac{L_{k}^{2}}{(L_{j}+L_{k}-L_{j}L_{k})^{2}}\cdot(Y_{j}-L_{j})+\frac{U_{j}^{2}}{(U_{j}+L_{k}-U_{j}L_{k})^{2}}\cdot(Y_{k}-L_{k})+\frac{L_{j}L_{k}}{L_{j}+L_{k}-L_{j}L_{k}}$

To see that these cuts are valid, note that function $f_{3}(x,y)$ is concave
on $x$ for a fixed $y$, so in particular, this is a valid upper bound for
$y=L_{k}$. On the other hand, the partial derivative $\frac{\partial
f_{3}}{\partial y}$ is increasing with respect to $x$, attaining its maximum
value on $x=U_{j}$. Since $f_{3}$ is also concave for a fixed $x$, then it is
a valid bound for $x=U_{j}$ and then for all
$(x,y)\in[L_{j},U_{j}]\times[L_{k},U_{k}]$.

Interchanging the roles of $x$ and $y$ in $f_{3}(x,y)$, we can also bound
$Y_{i}$ for a series composition by

$\displaystyle Y_{i}$ $\displaystyle\leq\left.\frac{\partial f_{3}}{\partial
x}\right|_{(L_{j},U_{k})}\cdot(Y_{j}-L_{j})+\left.\frac{\partial
f_{3}}{\partial
y}\right|_{(L_{j},L_{k})}\cdot(Y_{k}-L_{k})+f_{3}(L_{j},L_{k})$
$\displaystyle=\frac{U_{k}^{2}}{(L_{j}+U_{k}-L_{j}U_{k})^{2}}\cdot(Y_{j}-L_{j})+\frac{L_{j}^{2}}{(L_{j}+L_{k}-L_{j}L_{k})^{2}}\cdot(Y_{k}-L_{k})+\frac{L_{j}L_{k}}{L_{j}+L_{k}-L_{j}L_{k}}$

See Figure 3 for an example on how these hyperplanes improve the
overestimation of $f_{3}$ on $(L_{j},L_{k})$.

Figure 3. Function $\tfrac{xy}{x+y-x\cdot y}$ for $0.3\leq x\leq 0.8$ and
$0.4\leq y\leq 0.9$ (in orange), its concave envelope over
$[0,0.8]\times[0,0.9]$ (in blue) and the additional tangent hyperplanes at
$(0.3,0,4)$ (in green)

###  3.6. Computational experiments for the true optimal solution for the
problem

In this second set of experiments, we now include the combinatorial Benders
cut (8) to ensure that the optimal solution provides the true reliability for
the problem. We include these cuts for the previous two configurations convex
envelope cuts and without cuts, and we add a third configuration:

Improved envelope cuts::

The model described in Section 3.3, including the improvements based on the
local bounds provided by the branch-and-bound tree (see Subsection 3.5.2).

Combinatorial cuts for computing the exact reliability (§3.5.1) are
implemented as LazyContraintCallback, and improved cuts in branch-and-bound
(§3.5.2) are implemented as UserCutCallback. For this set of experiments, we
set a time limit of 3 hours for each problem.

Additionally, to benchmark our models with other solvers, we solve the
original model (1a)-(1m) using the MINLP solver SCIP v7.02 [16] compiled with
the parameters for better performance for this kind of nonlinear nonconvex
problem.

Figure 4. Performance profiles of the different configurations ($|E|=50$)

Figure 4 shows the performance profiles of the different configurations for
instances with $|E|=50$ edges. The figure on the left shows the percentage of
instances solved up to optimality before a given time ($x$-axis). For the
instances that are not solved within the time limit, the figure on the right
shows the percentage of instances attaining a given optimality gap (in log
scale).

Let us first analyze the instances for $\alpha=0.8$. As shown in previous
experiments, for these cases, all methods provide a good approximation of the
true and optimal reliabilities. Therefore, the configurations without the
improved envelope cuts behave well, solving 90% the instances within the time
limit of 3 hours, and the unsolved instances present a small optimality gap—a
2% gap in the worst case. However, note that there is a nonnegligible portion
of instances that require more than an hour to solve; this indicates that the
model is able to find a good solution but that it cannot quickly prove its
optimality because it needs to visit a large branch-and-bound tree to discard
all other potential solutions. On the other hand, the improved envelope cuts
behave drastically differently, solving all instances in less than one minute.
This can be explained by this configuration’s ability to locally adapt concave
envelopes during the branch-and-bound tree, providing better estimations and
thus better bounds, which yield a smaller branch-and-bound tree.

This better approximation of the nonconcave functions is even more relevant
for $\alpha=0.6$. While the combinatorial Benders cuts help in fixing the
mismatch between the reported reliability and the true reliability discussed
in Section 3.4.2, the weak estimations provided by the global approximation of
the functions prevent the solver from finding better solutions and/or proving
optimality efficiently: only 14% of the instances are solved within the time
limit when the improved envelope cuts are not included. This number increases
to 88% when these cuts are included, with most of these instances being solved
in just a few minutes.

Figure 5. Performance profiles of the different configurations ($|E|=100$)

The dominance of the improved envelope cuts also occurs for the instances with
$|E|=100$ edges. We present these results in Figure 5. We first note that the
problem is considerably harder to solve. For instance, the configuration
_without cuts_ can solve only 31% of the instances for $\alpha=0.8$ within the
time limit of three hours. Adding the convex envelopes of Theorem 1 improves
this metric, but only marginally. Nevertheless, the optimality gap obtained by
these configurations is good, with more than 95% of the instances finishing
with a gap of less than 1%. In the improved envelope cuts setting, 96% of the
instances are solved for $\alpha=0.8$, and the remaining instances finish with
an optimality gap less than $0.1\%$.

The problem becomes much more challenging for $\alpha=0.6$ and 100 edges. This
is expected, because due to the sequential construction process of the graph,
the differences between the nonlinear functions and their concave envelopes
are propagated into the overall approximation quality and become more relevant
when the number of steps in the construction sequence (i.e., the number of
edges) is large. In fact, without considering the improved envelope cuts, the
solver is not able to solve any instance, and the optimality gaps are
substantial for most of the cases. The performance improves when including the
improved envelope cuts, resulting in 17% of the instances being solved and
obtaining better optimality gaps for the unsolved instances.

####  3.6.1. Comparison with MINLP solver

To benchmark the proposed model against current state-of-art solvers for
nonlinear optimization models, we also solve the problem using the SCIP
solver. SCIP is among the best general-purpose solvers that are able to deal
with nonconvex constraints. It implements multiple bounding techniques, some
of which are similar to those studied in this paper, along with spatial
branch-and-bound based on linear outer-approximations of the problem. For more
details, see [32].

Figures 4 and 5 show performance profiles of SCIP in comparison with our
approach. We can see that for $|E|=50$, SCIP behaves in a similar way to the
improved envelope cuts, being slightly slower for $\alpha=0.8$. However, for
larger problems with $|E|=100$ edges, SCIP’s performance decreases
considerably, and it is outperformed by our proposed improved envelope cuts.
SCIP can solve only half of the instances for $\alpha=0.8$ and only 6% of the
instances for $\alpha=0.6$, reaching the time limit with optimality gaps that
are worse than the basic configuration without cuts, in most cases.

Figure 6. Final optimality gaps and number of nodes traversed in the branch-
and-bound tree ($|E|=100$, $\alpha=0.6$)

Figure 6 shows the optimality gaps (log scale) versus the best-bound objective
value obtained by the improved envelope cuts setting and SCIP on the unsolved
instances for the case of $100$ edges and $\alpha=0.6$. This figure indicates
that the problems become harder when the reliability of the problem is lower:
the optimality gaps are large when the objective bound is low. This can be
explained because this is the region where the difference between the
nonconvex functions and their concave envelopes differ the most (recall Figure
2), so the approximation is not sufficiently tight to lead the solver to prove
the optimality of the solutions. This is also correlated with the number of
nodes traversed in the branch-and-bound: the number of branch-and-bound nodes
are smaller in harder instances, indicating that the subproblems at each node
are harder to solve (probably because they include a larger number of cuts).
Interestingly, similar behaviors also occur in SCIP, even if the latter is
able to visit 10 times more nodes of the branch-and-bound tree.

##  4\. Conclusions and further extensions

We provide an optimization framework to solve network design problems for
maximizing the all-terminal reliability problem on series-parallel graphs when
failure probabilities are independent but not identical. Our approach exploits
the use of concave envelopes of the nonconcave functions that can be
implemented successfully using current optimization solvers, something that
has not been explored thus far in this context.

The special properties of the functions that appear in reliability
optimization allow us to derive envelopes that can be refined and exploited in
the solution process. Computational experiments show that it is highly
beneficial to perform such refinements of the concave envelopes along the
branch-and-bound process and thus provide better local approximations for the
nonlinear functions. If this is not done, the solver faces difficulties in
obtaining good solutions or proving optimality.

These techniques can be extended to more general contexts of network
reliability optimization. For example, similar ideas can be used for
$K$-terminal reliabilities, where the functions associated to other
reliability-preserving reductions (see [28]) could also be approximated by
their concave envelopes in a similar way. However, these concave envelopes are
not known and can be difficult to find in closed form. Therefore, developing
new techniques, such as those presented in [5], that can handle these
functions is a promising future direction that can considerably widen the
applicability of our proposed framework. Additionally, the use of
convex/concave envelopes for reliability optimization can also be applied to
more general families of graphs. In fact, the reductions discussed in this
paper apply to any graph and allow us to reduce the size of the problem. The
smaller problem can be solved using other optimization techniques such as
sample average approximation. This appears very promising, in particular for
graphs with small treewidth, as recently discussed in [17].

## References

  * [1] Hosam MF AboElFotoh and Loulwa S Al-Sumait, _A neural approach to topological optimization of communication networks, with reliability constraints_ , IEEE Transactions on Reliability 50 (2001), no. 4, 397–408.
  * [2] Kassie Archer, Christina Graves, and David Milan, _Classes of uniformly most reliable graphs for all-terminal reliability_ , Discret. Appl. Math. 267 (2019), 12–29.
  * [3] Michael O. Ball and J. Scott Provan, _The complexity of counting cuts and of computing the probability that a graph is connected_ , SIAM J. Computing 12 (1983), 777–788.
  * [4] Javiera Barrera, Héctor Cancela, and Eduardo Moreno, _Topological optimization of reliable networks under dependent failures_ , Operations Research Letters 43 (2015), no. 2, 132 – 136.
  * [5] Javiera Barrera, Eduardo Moreno, and Gonzalo Muñoz, _Convex envelopes for ray-concave functions_ , arXiV 2105.03532 (2021).
  * [6] Francis T. Boesch, _On unreliability polynomials and graph connectivity in reliable network synthesis_ , Journal of Graph Theory 10 (1986), no. 3, 339–352.
  * [7] F.T. Boesch, A. Satyanarayana, and C.L. Suffel, _A survey of some network reliability analysis and synthesis results_ , Networks 54 (2009), no. 2, 99–107.
  * [8] Zdravko I. Botev, Pierre L’Ecuyer, Richard J. Simard, and Bruno Tuffin, _Static network reliability estimation under the Marshall-Olkin copula_ , ACM Transactions on Modeling and Computer Simulation 26 (2016), no. 2, 14:1–14:28.
  * [9] Jason I Brown, Charles J Colbourn, Danielle Cox, Christina Graves, and Lucas Mol, _Network reliability: Heading out on the highway_ , Networks 77 (2021), no. 1, 146–160.
  * [10] Michael Bynum, Anya Castillo, Jean-Paul Watson, and Carl D Laird, _Strengthened SOCP relaxations for ACOPF with McCormick envelopes and bounds tightening_ , Computer Aided Chemical Engineering 44 (2018), 1555–1560.
  * [11] Gianni Codato and Matteo Fischetti, _Combinatorial benders’ cuts for mixed-integer linear programming_ , Operations Research 54 (2006), no. 4, 756–766.
  * [12] Darren L Deeter and Alice E Smith, _Economic design of reliable networks_ , IIE Transactions 30 (1998), no. 12, 1161–1174.
  * [13] Berna Dengiz and Cigdem Alabas-Uslu, _A self-tuning heuristic for the design of communication networks_ , Journal of the Operational Research Society 66 (2015), no. 7, 1101–1114.
  * [14] Berna Dengiz, Fulya Altiparmak, and Onder Belgin, _Design of reliable communication networks: A hybrid ant colony optimization algorithm_ , IIE Transactions 42 (2010), no. 4, 273–287.
  * [15] Berna Dengiz, Fulya Altiparmak, and Alice E Smith, _Efficient optimization of all-terminal reliable networks, using an evolutionary approach_ , IEEE Transactions on Reliability 46 (1997), no. 1, 18–26.
  * [16] Gerald Gamrath, Daniel Anderson, Ksenia Bestuzheva, Wei-Kun Chen, Leon Eifler, Maxime Gasse, Patrick Gemander, Ambros Gleixner, Leona Gottwald, Katrin Halbig, Gregor Hendel, Christopher Hojny, Thorsten Koch, Pierre Le Bodic, Stephen J. Maher, Frederic Matter, Matthias Miltenberger, Erik Mühmer, Benjamin Müller, Marc E. Pfetsch, Franziska Schlösser, Felipe Serrano, Yuji Shinano, Christine Tawfik, Stefan Vigerske, Fabian Wegscheider, Dieter Weninger, and Jakob Witzig, _The SCIP Optimization Suite 7.0_ , Technical report, Optimization Online, March 2020.
  * [17] Amir Kafshdar Goharshady and Fatemeh Mohammadi, _An efficient algorithm for computing network reliability in small treewidth_ , Reliability Engineering & System Safety 193 (2020), 106665.
  * [18] IBM, _Ibm decision optimization cplex modeling for python_ , 2019.
  * [19] Anton J Kleywegt, Alexander Shapiro, and Tito Homem-de Mello, _The sample average approximation method for stochastic discrete optimization_ , SIAM Journal on Optimization 12 (2002), no. 2, 479–502.
  * [20] Omar Matus, Javiera Barrera, Eduardo Moreno, and Gerardo Rubino, _On the Marshall-Olkin copula model for network reliability under dependent failures_ , IEEE Transactions on Reliability 68 (2019), no. 2, 451–461.
  * [21] Garth P McCormick, _Computability of global solutions to factorable nonconvex programs: Part I—Convex underestimating problems_ , Mathematical Programming 10 (1976), no. 1, 147–175.
  * [22] Edward F Moore and Claude E Shannon, _Reliable circuits using less reliable relays_ , Journal of the Franklin Institute 262 (1956), no. 3, 191–208.
  * [23] Eric M Neufeld and Charles J Colbourn, _The most reliable series-parallel networks_ , Networks 15 (1985), no. 1, 27–32.
  * [24] Omer Ozkan, Murat Ermis, and Ilker Bekmezci, _Reliable communication network design: The hybridisation of metaheuristics with the branch and bound method_ , Journal of the Operational Research Society 71 (2020), no. 5, 784–799.
  * [25] Hebert Pérez-Rosés, _Sixty years of network reliability_ , Mathematics in Computer Science 12 (2018), no. 3, 275–293.
  * [26] J.L. Pulsipher and V.M. Zavala, _Measuring and optimizing system reliability: a stochastic programming approach_ , TOP 28 (2020), 626–645.
  * [27] A. Rosenthal, _Computing the reliability of complex networks_ , SIAM Journal on Applied Mathematics 32 (1977), no. 2, 384–393.
  * [28] A. Satyanarayana and Mark K. Chang, _Network reliability and the factoring theorem_ , Networks 13 (1983), no. 1, 107–120.
  * [29] A. Satyanarayana and R. Wood, _A linear-time algorithm for computing k-terminal reliability in series-parallel networks_ , SIAM Journal on Computing 14 (1985), no. 4, 818–832.
  * [30] Yongjia Song and James R. Luedtke, _Branch-and-cut approaches for chance-constrained formulations of reliable network design problems_ , Mathematical Programming Computation 5 (2013), no. 4, 397–432.
  * [31] L. Valiant, _The complexity of enumeration and reliability problems_ , SIAM Journal on Computing 8 (1979), no. 3, 410–421.
  * [32] Stefan Vigerske and Ambros Gleixner, _SCIP: Global optimization of mixed-integer nonlinear programs in a branch-and-cut framework_ , Optimization Methods and Software 33 (2018), no. 3, 563–593.
  * [33] Yixin Ye, Ignacio E Grossmann, and Jose M Pinto, _Mixed-integer nonlinear programming models for optimal design of reliable chemical plants_ , Computers & Chemical Engineering 116 (2018), 3–16.
  * [34] Leila Zia and David W Coit, _Redundancy allocation for series-parallel systems using a column generation approach_ , IEEE Transactions on Reliability 59 (2010), no. 4, 706–717.