- Email: [email protected]
Minimum cost dominating tree sensor networks under probabilistic constraints
Accepted Manuscript
Minimum cost dominating tree sensor networks under probabilistic constraints Pablo Adasme, Rafael Andrade, Abdel Lisser PII: DOI: Reference:
S1389-1286(16)30387-5 10.1016/j.comnet.2016.11.011 COMPNW 6053
To appear in:
Computer Networks
Received date: Revised date: Accepted date:
6 May 2016 9 November 2016 12 November 2016
Please cite this article as: Pablo Adasme, Rafael Andrade, Abdel Lisser, Minimum cost dominating tree sensor networks under probabilistic constraints, Computer Networks (2016), doi: 10.1016/j.comnet.2016.11.011
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Minimum cost dominating tree sensor networks under probabilistic constraints
CR IP T
Pablo Adasmea , Rafael Andradeb , Abdel Lisserc a Departamento
AN US
de Ingenier´ıa El´ ectrica, Universidad de Santiago de Chile, Avenida Ecuador 3519, Santiago, Chile. [email protected] b Departamento de Estat´ ıstica e Matem´ atica Aplicada, Universidade Federal do Cear´ a, Campus do Pici, BL 910, CEP 60.455-760 - Fortaleza, Cear´ a, Brasil. [email protected] c Laboratoire de Recherche en Informatique, Universit´ e de Paris-Sud 11, Ada Lovelace, 91405 Paris, France. [email protected]
Abstract
There is an increasing interest in planning sensor networks by considering both the impact of distances among sensors and the risk that power consumption leads to a very small network lifetime. A sensor failure can affect sensors in its neighborhood and compromise the network data communication. Weather
M
conditions may cause the power consumption of data communication to vary with uncertainty. This work introduces a compact probabilistic optimization
ED
approach to handle this problem while considering jointly or separately dependence among power consumption of the links of the network in a unified framework. We explore the concept of copulas in a dominating arborescence
PT
(DA) model for directed graphs, extended accordingly to handle the uncertain parameters. We give a proof of the DA model correctness and show that it can
CE
solve to optimality some benchmark instances of the deterministic dominating tree problem. Numerical results for the probabilistic approach show that our model tackles randomly generated instances with up to 120 nodes.
AC
Keywords: Dominating tree, probabilistic constraints, compact formulation, mixed integer linear programming, sensor network design
Preprint submitted to Journal of Computer Networks
November 14, 2016
ACCEPTED MANUSCRIPT
1. Introduction A classic version of the dominating tree (DT) sensor network problem [24] concerns determining a minimum cost tree T in an undirected graph G = (V, E)
5
CR IP T
of set of nodes V and set of edges E, with ce ≥ 0 being the cost (weight-distance) of edge e ∈ E, where any node v ∈ V is in T or is adjacent to a node in this tree.
Link cost distances are related as usual [24] with the power required to establish
the topology of the involved dominating tree connections. The topology cost of
a tree T is the sum of the cost distances of its edges. The DT problem is known to be NP-Hard [24].
AN US
In some sensor networks, with each edge {u, v} ∈ E, of extremities u, v ∈ V ,
10
there is associated an uncertain power consumption dξuv necessary to maintain the data transmission between u and v operational, where ξ denotes a random variable of known probability distribution function. Since power consumption in each sensor may vary due to weather or radio interference, sometimes the available power (battery autonomy of sensors) to keep the network operational
M
15
becomes insufficient. This may cause loss of communication among the sensors. To overcome this problem, in this paper we consider probabilistic constraints
ED
related to power consumption in the classic DT problem. The resulting probabilistically constrained dominating tree (PCDT) problem consists in finding a minimum cost distance dominating tree T = (V (T ), E(T )) of G, with set of
PT
20
nodes V (T ) and set of edges E(T ) subject to a probabilistic constraint related to the risk that the total power consumption of each sensor u ∈ V (T ) exceeds a
CE
given limit dm corresponding to its power autonomy. The power consumption of P a sensor u is defined as the sum {u,v}∈E(T ) dξuv of all links in E(T ) connected to u in T .
AC
25
Another classic optimization approach for planning sensor networks consid-
ers the minimum number of connected dominating sensors [11]. In this case, the problem is equivalent to the one in [4, 24] when the distance between sensors equals 1. The novelty here is to provide new ideas for the development of
30
communication protocols for sensor networks by considering both the impact of
ACCEPTED MANUSCRIPT
distances among sensors and the risk that power consumption of sensors exceeds a power threshold available to keep the network operational. In this sense, our work is a first attempt to answer this question by considering jointly or sepa-
35
CR IP T
rately dependence among power consumption of the links of the network in a unified framework.
Both the stochastic mathematical optimization model and the theory to
make it tractable by mixed integer linear programming are contributions to future algorithm developments for practical applications. In this sense, this
step of our research focuses mainly on the theory rather than implementing and
evaluating the overhead of the proposed models in practice. Notice however,
AN US
40
that any dominating tree topology proposed as solution by our model can be used in real implementations as well as in [4, 11, 24]. The reader is referred to [8] for further details on the involved problems in sensor networks. 1.1. An application of a dominating tree sensor network
In Figure 1, we show an example of an undirected connected sensor net-
M
45
work. The links of the graph in Figure 1 represent the possible data exchange
ED
among the 8 sensors {A, B, . . . , H}. Assume that two sensors connected by a link are at a unit distance of each other and that values near each communication link represent, deterministically the power consumption required to keep each link operational. The objective is to determine the minimum distance
PT
50
dominating tree topology for the communication among sensors respecting the power autonomy of each sensor. In fact, if the power autonomy of each sensor is
CE
unlimited, then the optimal dominating tree topology should be the unique link T = ({A, B}) of unit distance, since minimizing the number of connected dominating nodes solves this problem [11]. However, if with each sensor we associate
AC
55
a battery autonomy limit of 5 units to be used with all its communication links, then T = ({A, B}) is no more feasible because the power consumption required by this link is equal to 5.9 units and this value exceeds the available power for sensors A and B. By inspection, the optimal solution for this example is the
60
dominating tree topology T 0 = ({A, C}, {C, D}, {B, D}), of distance equal to 3,
ACCEPTED MANUSCRIPT
with total link power consumption of sensors A, B, C, and D being equal to 1, 1, 5, and 5, respectively. Of course, we assume that the remaining nodes not in T 0 will communicate with their nearest neighbors in T 0 and that the required
E
4
F
1 2
1
H
2 5.9
A C
4
B 1 D
1
G
CR IP T
power to accomplish this task is negligible for the problem.
AN US
Figure 1: Example of an undirected connected sensor network.
Thus, in order to deal with the probabilistic constraints when power con-
65
sumption of the existing links vary according to, e.g., extreme and uncertain weather conditions, we propose a polynomial compact formulation for the DT
M
problem which is based on dominating arborescence of a directed graph and give a proof of its correctness. The compact model allows us to obtain determin70
istic equivalent mixed integer linear programming formulations for the PCDT
ED
problem. Numerical results show that the new approach allows us to solve to optimality random and benchmark instances for the DT and PCDT problems
PT
with up to 120 nodes. To the best of our knowledge, this is the first work reporting proven optimal solutions for challenging benchmark instances of the classic 75
DT problem [4, 26]. Moreover, the probabilistic modelling approach allows us
CE
to handle joint and separate probabilistic constraints in a unified framework.
AC
1.2. Probabilistic constraints Stochastic programming is a powerful technique to deal with the uncertainty
of some input parameters of a mathematical programming problem [23]. These
80
parameters are assumed to behave as random variables which are distributed according to a given probability distribution function. In particular, the probabilistically constrained approach imposes a probability of occurrence for some
ACCEPTED MANUSCRIPT
(or all) of the constraints of a mathematical model. This means that some of the constraints will be satisfied, at least for a given percentage [1, 5, 10, 15, 21, 22]. 85
As far as we know, no attempt has been made to deal with probabilistic con-
CR IP T
straints in the problem of designing sensor networks with a dominating tree topology under uncertain power consumption. Probabilistic constraints can be considered either separately [1, 3, 5, 10, 15, 16, 19, 20, 21, 22] or jointly [5, 6, 7, 17].
In this paper, we consider both joint and separate probabilistic constraints
90
for the PCDT problem. For this purpose, we assume that the row vectors of the
AN US
matrix (dξij ) are joint dependent multivariate normally distributed vectors with
known means and covariance matrices. The dependence of the random vectors is handled by means of copulas [6, 13]. A copula is a distribution function 95
β : [0, 1]|V | → [0, 1] of some |V |-dimensional random vector whose marginals are uniformly distributed on [0, 1]. We apply useful properties of Gumbel-Hougaard copulas to describe the dependence between the rows of the matrix (dξuv ). The
M
latter allows us to handle joint and separate probabilistic constraints in a unified framework. The theory of copula distributions was developed in probability theory and mathematical statistics [13, 18]. Copula distributions allow us to find
ED
100
deterministic convex reformulations for probabilistic constraints while simultaneously handling dependent random variables. In the stochastic programming
PT
field, this has not been well developed yet. In fact, this is the most important reason for using copula distributions so far. Several copula distributions have been proposed such as Clayton, Gumbel and Frank copulas, among others. In
CE
105
this paper, we mainly focus on the use of Gumbel-Hougaard copula as it is the one that leads to more tractable deterministic reformulations as shown in [6],
AC
and also because it allows a unified framework to handle joint and separate probabilistic constraints simultaneously. In Section 4, we give a brief formal
110
description of this copula. In [13], the authors use copulas to come up with convexity results for chance
constrained problems with dependent random right hand side whereas in [6] the authors apply Gumbel-Hougaard copulas to a quadratic optimization prob-
ACCEPTED MANUSCRIPT
lem subject to joint probabilistic constraints. In this paper, we transform the 115
probabilistic constraints into equivalent deterministic second order conic constraints that we further linearize. Consequently, we obtain deterministic equiv-
CR IP T
alent mixed integer linear programming (MILP) formulations for the PCDT problem. Additionally, we investigate the impact of using valid inequalities referred to as generalized sub-tour elimination constraints in all our proposed 120
models [11].
The paper is organized as follows. In Section 2, we introduce both separate and joint probabilistic constraints for the problem. Subsequently, we introduce
AN US
a new compact polynomial formulation for the DT problem and give a proof of its correctness. Then, in Sections 3 and 4 we present deterministic equiva125
lent MILP formulations using separate and joint probabilistic constraints, respectively. Then, in Section 5 we conduct numerical experiments in order to compare all the proposed models for random and benchmark [26] instances for the DT and PCDT problems. Finally, in Section 6 we conclude the paper and
2. Problem formulation
ED
130
M
provide some insights for future research.
In this section, we present a new variant of the classic DT problem while
PT
considering both separate and joint probabilistic constraints that we incorporate in the proposed DT model leading to the PCDT problem formulation. Then, we introduce a compact polynomial formulation for the problem which is based on arborescence of directed graphs and give a proof of its correctness.
CE
135
AC
2.1. Probabilistic dominating tree problem Let G = (V, E) be an undirected sensor network with set of nodes V and
set of connections E. We propose working with an equivalent directed graph H = (V, A), where the set of arcs A is obtained from E by replacing each edge
140
{u, v} ∈ E by the two arcs (u, v) and (v, u) ∈ A. We associate with each arc
(u, v) ∈ A a cost cuv > 0 and an uncertain power consumption dξuv > 0 given
ACCEPTED MANUSCRIPT
by a random variable ξ of known probability distribution function. We assume that the matrices c = (cuv ) and dξ = (dξuv ) are symmetric. Given a directed graph H = (V, A), an arborescence T of H is a rooted directed acyclic subgraph of H in which, given a root vertex r ∈ V , for any other vertex v ∈ V \ {r}, there
CR IP T
145
is exactly one directed path from r to v in T . An arborescence T of a directed
graph H is dominating if every node of H or is in T or is adjacent to a node
in T . The problem is to find a minimum cost dominating arborescence (DA) T of H where the probability that the sum of power consumption of the arcs 150
adjacent to each node u of T exceeds a positive maximum value dm is limited
AN US
above by a certain risk parameter α ∈ [0, 0.5). Notice that the parameter α must
lie in this interval in order to keep the convexity of the deterministic equivalent probabilistic constraints throughout the paper [2, 6, 15].
We represent a DA T of H by a vector x ∈ {0, 1}|A| where xuv = 1 if arc 155
(u, v) belongs to T and xuv = 0, otherwise. We use binary variables yv , for all v ∈ V , to indicate if a node v belongs to T . If it is the case, we set yv = 1;
M
− + and set yv = 0, otherwise. Define N(v) = {u|(u, v) ∈ A}, N(v) = {u|(v, u) ∈ A},
− + N(v) = N(v) ∪N(v) , and N[v] = N(v) ∪{v}. The following probabilistic constraints
ED
model the involved uncertainties X P dξuv xuv ≤ dm , ∀u ∈ V ≥ (1 − α) v∈N +
(1)
160
PT
(u)
where P {·} denotes a probability measure. We assume that the row vectors of
CE
the input matrix d = (duv (ξ)) are dependent multivariate normally distributed vectors with known means and covariance matrices. The probabilistic constraint (1) imposes the condition that each node power constraint must be satisfied at
AC
least for (1 − α)% of the cases where α represents the risk of not satisfying
165
some of these constraints. We remark that the probabilistic constraints (1) are formulated as joint probabilistic constraints under the assumption of row dependence. Next, we further consider the case of separate probabilistic constraints where the row vectors in matrix (dξuv ) are independent multivariate normally
ACCEPTED MANUSCRIPT
distributed vectors. In this case, constraints (1) can be replaced by X P dξuv xuv ≤ dm ≥ (1 − α), ∀u ∈ V v∈N +
(2)
170
CR IP T
(u)
Notice that the statement “∀u ∈ V ” in (1) is inside the probability notation.
This means that the probabilities should be jointly satisfied ∀u ∈ V whereas
in (2), they should be satisfied in a separate way. We handle the dependence between these random vectors by means of copulas [6, 13]. Now, we propose a new compact polynomial formulation for the DT problem that allows to handle efficiently the constraints (1) and (2) as it is based on directed graphs.
AN US
175
2.2. A compact polynomial formulation for the DT problem
Let associate with each node v ∈ V a non-negative continuous variable πv and let fv represent the degree of node v in the undirected graph G. A polynomial compact formulation for the dominating arborescence (DA) problem in the corresponding directed graph H is DTc :
M
180
X
min
{x,y,π}
(3)
X
(4)
ED
X
cuv xuv ,
(u,v)∈A
xuv =
v∈V
(u,v)∈A
yv − 1,
AC
CE
PT
πv − πu + |V |(1 − xuv ) ≥ 1,
∀(u, v) ∈ A,
(5)
πv − πu − |V |(1 − xuv ) ≤ 1, ∀(u, v) ∈ A, X X xvu + xuv ≥ yv , ∀v ∈ V,
(6)
X
+ u∈N(v)
X
u∈N(v)
X
− u∈N(v)
(7)
− u∈N(v)
+ u∈N(v)
xvu +
X
− u∈N(v)
yu ≥ 1,
xuv ≤ fv yv ,
∀v ∈ V,
xuv ≤ 1,
1 ≤ πv ≤ |V | − 2,
∀v ∈ V, ∀v ∈ V,
x ∈ {0, 1}|A| , y ∈ {0, 1}|V | .
∀v ∈ V,
(8) (9) (10) (11) (12)
ACCEPTED MANUSCRIPT
where cuv represents the cost of arc (u, v) ∈ A;
CR IP T
xuv is a binary variable denoting if arc (u, v) belongs (xuv = 1) or not (xuv = 0) to the DA; 185
πv is a continuous variable corresponding to the position in the sequence of arcs starting from the vertex of the DA that has no incoming arc to v; and
yv is a binary variable indicating if vertex v belongs (yv = 1) or not (yv = 0)
190
AN US
to the DA;
The total cost of a feasible DA T is minimized in the objective function (3). P Constraint (4) establishes that T has v∈V yv − 1 arcs. Constraints (5)-(6) ensure that if an arc (u, v) belongs to the solution, then πv − πu = 1; otherwise
both constraints become redundant. These constraints together with (10) avoid cycles in any feasible solution. The idea behind these constraints is to explore
195
M
simultaneously the space of primal and dual variables of the problem structure, based on a distance property. Indeed, every two end nodes of any arc in the
ED
solution are at a unit distance from each other. In fact, for this problem, we can show that (6) are redundant in the presence of (10). Concerning constraints (7), they establish that if node v is in the solution,
200
PT
then at least one arc is adjacent to v. Constraints (8) are valid inequalities for DTc , indicating that if node v is not in the solution, then no arc is adjacent to
CE
v; otherwise, at most fv arcs can be adjacent to this node. Constraints (9) are domination constraints and ensure that any node v in V has a neighbor in T . Constraints (10) limit the number of arcs arriving at a node v ∈ V to at
AC
most one. Constraints (11) bound the possible values for the potential variables
205
π. This is due to the fact that any directed path in the solution arborescence spans at most |V | − 2 nodes. In fact, since DTc is a minimization problem and
the π variables have one as lower bound, the upper bound |V | − 2 on these
variables is redundant. Finally, (12) are the domain of the decision variables.
ACCEPTED MANUSCRIPT
Theorem 1. An optimal solution of DTc referred to H gives a dominating tree 210
of G of minimum cost. Proof. Initially, we show that any feasible solution (¯ x, y¯, π ¯ ) of DTc is an
CR IP T
acyclic connected arborescence of H. Suppose that x ¯ induces a cycle of arcs
(v1 , v2 ), (v2 , v3 ), · · · , (vp , v1 ), with vi ∈ V for i = 1, . . . , p, for some p ≤ |V | − 3.
By constraints (5) and (6), πv2 = πv1 + 1, πv3 = πv2 + 1, . . ., πvp = πvp−1 + 1 and πv1 = πvp + 1. This implies that πv1 < πv2 < · · · < πvp < πv1 cannot hold.
Thus, the directed graph induced by x ¯ is acyclic. Now suppose that x ¯ is acyclic
and w.l.o.g. induces two connected arborescences whose sets of nodes induced
AN US
by the non null components of y¯ are denoted by D1 ⊂ V and D2 ⊂ V , with
D1 ∪ D2 being a non-connected dominating set of G. By constraints (7), both extremities u and v of an arc (u, v) such that x ¯uv = 1 belong to the solution, i.e., y¯u = y¯v = 1. Note that arcs corresponding to non null components of x ¯ have both extremities or in D1 , denoted by A(D1 ), or in D2 , denoted by A(D2 ).
X
x ¯uv = |D1 | − 1,
X
x ¯uv = |D2 | − 1,
ED
connected, we have that
M
As the directed graphs induced by D1 and D2 are assumed to be individually
(u,v)∈A(D1 )
PT
and
CE
Moreover,
X
(u,v)∈A(D2 )
x ¯uv =
(u,v)∈A
X
(u,v)∈A(D1 )
and
AC
X
v∈D1
Consequently,
(u,v)∈A
x ¯uv =
X
v∈D1
y¯v − 1 +
x ¯uv ,
y¯v = |D1 |, y¯v = |D2 |.
X
y¯v − 1 =
v∈D2
X
(u,v)∈A(D2 )
X
v∈D2
X
x ¯uv +
X
v∈D1 ∪D2
y¯v − 2 =
X
v∈V
y¯v − 2,
ACCEPTED MANUSCRIPT
which violates constraint (4). Therefore, any optimal solution of DTc must be a connected arborescence of G. Finally, by dropping the orientation of its arcs we obtain a dominating tree of minimum cost of G.
215
CR IP T
In Figure 2 we show some ingredients involved in the proof of Theorem 1.
Initially, we determine the corresponding directed graph H from the undirected graph in Figure 1. A bi-directed arrow represents two arcs that replace the undirected link between its corresponding sensors.
For the sake of comparison with PCDT, the power required by each arc is
220
AN US
equal to that of the original link as depicted in the directed graph H. Assume
that the power available at each sensor is of 5 units. The cost between two sensors is equal to 1. We show five situations avoided or allowed by this theorem to show how an optimal DA leads to an optimal DT while guaranteeing acyclic and connected solutions.
If we were interested in monitoring the sensor network area by using communication protocols as suggested in [8], where the main purpose is to determine
M
225
the minimum connected dominating set of sensors [11], the solution should be the arborescence {(B, A)} in (i), of unitary cost distance, constituted of only one
ED
arc and two dominating sensors B and A. The remaining sensors are dominated and communicate with the dominating ones by using the dotted links representing their nearest neighbor (sensor) in this arborescence. However, if we intend
PT
230
to consider the interference of natural events of the sensor network on the power consumption between sensors B and A (because the required power to keep this
CE
communication link exceeds the available power at these sensors), this solution is no longer feasible. On the other hand, in (ii) when we relax the connectivity requirement, we obtain two arborescences {(A, C)} and {(D, B)} that together
AC
235
constitute a (non-connected) dominating set presenting minimum cost distance equal to two. Their arcs power consumption respect the battery autonomy of the involved sensors, but the number of arcs is two while the number of dominating sensors is four, thus violating constraint (4). Unfortunately, as we can
240
see in this situation, the use of cut inequalities [11] is not helpful to force ob-
ACCEPTED MANUSCRIPT
taining a connected solution because the minimum number of non-dominating sensors needed to connect these two sets is zero. Recall that the term “cut inequality” [11] is a valid inequality defined on
245
CR IP T
the decision variables y to cut off an integer solution associated with a proper subset S = {i|yi = 1} of vertices inducing a disconnected graph. If a subset of
vertices S is not connected, then at least one vertex in S¯ = V \ S must belong to P the solution and the corresponding cut inequality is i∈S¯ yi ≥ 1. These valid inequalities are usually used within exact algorithms [11].
Subsequently, notice that individually the sets {A, C} and {D, B} are not dominating. Although, their corresponding arborescences can be connected by
AN US
250
adding, for instance, the arc (C, D). We cannot impose that at least an arc must be connected to the set {A, C} or to {D, B} because their complement with respect to the remaining sensors are dominating sets [11]. When a set S and ¯ are both non-dominating, at least an arc must connect S its complement, say S, 255
to S¯ [11]. In (iii), we have a fictitious cycle induced by the arcs (C, A), (A, B),
M
(C, D) and (D, B). We show why it cannot belong to a feasible DA. Initially, if we assume that πC = π ¯ , constraints (5) and (6) impose that πA = πD = π ¯+1
ED
and πB = π ¯ + 2. Thus, these constraints are not sufficient to avoid cycles of this type in the absence of constraints (10) that are required for the model to 260
be correct. The optimal DA for this example is {(C, A), (C, D), (D, B)}, rooted
PT
at sensor C as depicted in (iv), and has cost distance equal to three. The sum of the arcs power consumption at sensors A, B, C, and D do not exceed the
CE
available power at each sensor. The values for the π variables associated with this solution can be πC = 1, πA = πD = 2, πB = 3. The remaining values of
265
the π variables are irrelevant for this solution. By dropping the orientation of
AC
the arcs in this DA, we obtain the DT in (v) that is the optimal solution for the undirected graph in Figure 1. The fact that DTc is represented as an arborescence of a directed graph
allows us to handle efficiently the probabilistic constraints (1) and (2).
270
Hereafter, we refer to the PCDT problem as the problem obtained by the constraints defining the DTc model with the additional probabilistic constraints
ACCEPTED MANUSCRIPT
Digraph H
4 E
F 5.9
A
B
2
H
G C
E
F
A
B
1 C
F
5.9
A
B 1
H
D
G
D
4 E
F 5.9
A
C
G
H
D
1
C
4
E
E
F
A 1
1 C
D
4
F
A
G
H
B
1
1
C
M
H
B
G
D
(iii)
(ii)
(i)
B
1
AN US
E
1
1
H
1
CR IP T
2
1
4
G
D
(v)
(iv)
Figure 2: Example of how Theorem 1 works to determine a dominating tree from an optimal
ED
dominating arborescence while avoiding infeasible solutions.
PT
(2) or (1), for the separate or joint probabilistically constrained approaches, respectively. We further consider the resulting MILP model obtained from DTc
AC
CE
while replacing constraints (8) with the following set of constraints
275
X
xuv
≤
yu ,
∀(u, v) ∈ A,
(13)
xuv
≤
(14)
xuv
≤
yv , ∀(u, v) ∈ A, X yj − 1, ∀i ∈ V ,
(u,v)∈A(N[i] )
(15)
j∈N[i]
and constraints (10) with the following ones X
− u∈N(v)
xuv ≤ yv ,
∀v ∈ V.
(16)
ACCEPTED MANUSCRIPT
We denote this resulting MILP model by DTcV i . In particular, constraints (15) are referred to as a strengthened version of the generalized sub-tour elimination constraints [11]. Finally, we denote hereafter the corresponding linear
280
CR IP T
programming (LP) relaxations of DTc and DTcV i by LPc and LPcV i , respectively. 3. DT problem with separate probabilistic constraints
In this section, we propose a deterministic equivalent MILP formulation for
the PCDT problem using separate probabilistic constraints. For this purpose,
we transform the probabilistic constraints (2) into deterministic equivalent sec-
285
AN US
ond order conic constraints [15]. Subsequently, we linearize the conic constraints to obtain an MILP formulation for the PCDT problem.
In order to obtain a deterministic MILP formulation, we assume that the row vectors in matrix (dξuv ) in (2) are independent multivariate random variables normally distributed with known means d¯uv , ∀(u, v) ∈ A. Also, let (σ u ) = 290
M
u (σiv ), ∀ (i, v) ∈ A, u ∈ V , be the root square matrix of the corresponding covari-
ance matrices of each row vector (dξu• ), ∀ u ∈ V , that we assume to be positive
semidefinite. In this case, notice that the constraints (2) can be equivalently
ED
written as [15, 22]
PT
v u u uX X d¯uv xuv + F −1 (1 − α)u t
+ v∈N(u)
i∈V
X
+ v∈N(u) \{i}
2
ux ≤d , σiv uv m
∀ u ∈ V, (17)
CE
where F −1 (·) denotes the inverse of the standard normal cumulative distribution function. The condition for the parameter α ∈ [0, 0.5) is mandatory for the constraints (17), otherwise these constraints would not be convex. Proposition 1. A deterministic equivalent MILP formulation for the PCDT
AC
295
problem with separate probabilistic constraints can be written as M IP1
min
{x,y,π,ϕ}
X
cuv xuv ,
(u,v)∈A
subject to (4) − (12) and
(18)
ACCEPTED MANUSCRIPT
(F −1 (1 − α))2 d2m +
X
i∈V
X
u 2 (σiv ) xuv +
+ v∈N(u) \{i}
(−2dm d¯uv + d¯2uv )xuv +
+ v∈N(u)
X
+ + v∈N(u) k∈N(u) \{i,v}
X
+ v∈N(u)
d¯uv d¯uk ϕkuv ,
(19)
∀u ∈ V,
xuv + xuk − 1 ≤ ϕkuv ≤ xuv
(20)
∀u, v, k ∈ V, (u 6= v 6= k),
∀u, v, k ∈ V, (u 6= v 6= k),
Proof. We write constraint (17) equivalently as [15, 20] v 2 u u uX X X ux ≤d − F −1 (1 − α)u σiv d¯uv xuv , uv m t + v∈N(u) \{i}
+ v∈N(u)
M
Notice that in order for (24) to make sense, we must have
(23)
∀u ∈ V. (24) P
+ v∈N(u)
d¯uv xuv ≤
dm , ∀u ∈ V . Thus, taking the square of both sides in (24) we obtain 2 2 X X X u (F −1 (1 − α))2 σiv xuv ≤ dm − d¯uv xuv ,
ED
300
(21)
(22)
AN US
ϕ ∈ {0, 1}|V ×V ×V | .
i∈V
u u k σiv σik ϕuv ≤
+ + v∈N(u) k∈N(u) \{v}
∀u ∈ V, X d¯uv xuv ≤ dm , ϕkuv ≤ xuk
X
X
CR IP T
X
i∈V
+ v∈N(u) \{i}
+ v∈N(u)
CE
PT
that after solving each quadratic term, is equivalent to X X X X u 2 (F −1 (1 − α))2 (σiv ) xuv + d2m +
X
AC
+ v∈N(u)
i∈V
+ v∈N(u) \{i}
(−2dm d¯uv + d¯2uv )xuv +
X
+ + v∈N(u) k∈N(u) \{i,v}
X
∀u ∈ V,
u u σiv σik xuv xuk ≤
d¯uv d¯uk xuv xuk ,
+ + v∈N(u) k∈N(u) \{v}
∀u ∈ V.
(25)
Finally, if we define ϕkuv = xuv xuk , ∀u, v, k ∈ V , with u 6= v 6= k, and use
standard methods [9, 12] to linearize ϕkuv as shown in (21) and (22), then (25) can be written as in (19). Finally, constraints (4)-(12) ensure that feasible
305
solutions are dominating arborescences, thus concluding the proof.
ACCEPTED MANUSCRIPT
We further consider the resulting MILP model obtained from M IP1 while replacing constraints (8) with the set of constraints (13)-(15) and constraints (10) with the constraints (16). Hereafter, we denote this MILP model by M IP1V i .
310
CR IP T
Finally, we denote the corresponding LP relaxations of M IP1 and M IP1V i by LP1 and LP1V i , respectively.
4. DT problem with joint probabilistic constraints
In this section, we obtain a deterministic equivalent MILP formulation for the PCDT problem using joint probabilistic constraints. We transform the prob-
315
AN US
abilistic constraints into deterministic equivalent second order conic constraints and linearize them in order to obtain an MILP formulation. In particular, our proposed model is based on the family of Gumbel-Hougaard copula. For a
deeper comprehension related to the theory of copulas we refer the reader to [25] and to the book [18]. The existence and uniqueness of a copula, for any
320
M
distribution function and all its marginals are guaranteed by Sklar’s Theorem [18]. In this paper, we consider the Gumbel-Hougaard copula which can be 1/θ |V | X θ (− ln tk ) βθ (t) = exp − k=1
ED
written as
(26)
β(t) =
|V | Y
tk
k=1
CE
PT
for given θ ≥ 1. In particular when θ = 1, this copula takes the form
and represents the joint distribution of independent random variables. Next,
AC
consider the following transformations P
+ v∈N(u)
ωu (x) = v u u uP t i∈V
(dξuv − d¯uv )xuv P
+ v∈N(u) \{i}
2
ux σiv uv
ACCEPTED MANUSCRIPT
dm − gu (x) = v u u uP t i∈V
where
P
i∈V
P
+ v∈N(u) \{i}
P
d¯uv xuv
+ v∈N(u)
P
+ v∈N(u) \{i}
2
2
ux σiv uv
CR IP T
and
u σiv xuv is strictly positive. The random variable ωu (x)
has one-dimensional standard normal distribution which is independent of the
variable x = (xuv ). Note that the probabilistic constraint (1) can be equivalently written as
AN US
325
P {ωu (x) ≤ gu (x), ∀u ∈ V } ≥ (1 − α)
(27)
Lemma 1. [6] If the random vector (ω1 (x), . . . , ω|V | (x))T , where ωu (x) has one-dimensional standard normal distribution, has a joint distribution driven
330
M
by the Gumbel-Hougaard copula βθ with some θ ≥ 1, then the constraint (1) is equivalent to the set of constraints
+ v∈N(u)
ED
v u u u X X 1/θ u d¯uv xuv + F −1 (1 − α)zu t
PT
∀u ∈ V X zu = 1, u∈V
zu ≥ 0,
∀u ∈ V
i∈V
X
+ v∈N(u) \{i}
2
ux ≤d , σiv uv m
(28) (29)
CE
where F −1 (·) is the inverse of the standard normal cumulative distribution function.
AC
Notice that when the term
P
+ v∈N(u)
xuv equals zero for a particular u ∈ V ,
constraints (28) are trivially satisfied. Moreover, the value of zu = 0 in this case
335
since the root square term in (28) equals zero as well. The following convexity lemma is also required Lemma 2. [6] If (1 − α) ≥ [0, 1].
1 2
1/θ and θ ≥ 1, then F −1 (1 − α)zu is convex on
ACCEPTED MANUSCRIPT
340
Hence, in this paper we assume that 12 ≤ (1 − α) < 1. Notice that when 1/θ θ → ∞, the function F −1 (1 − α)zu → F −1 (1 − α) which is equivalent to a separate probabilistic constraint [15, 22]. This shows that the Gumbel-
CR IP T
Hougaard copula allows us to handle joint and separate probabilistic constraints in a unified framework. Also notice that for values of θ ∈ [1, ∞) each zu , ∀u ∈ V
in constraint (29) will always be positive, otherwise the constraint (28) would 345
never be satisfied since F −1 (1) → ∞, unless the root square term equals zero.
In this case, we assume that the row vectors in matrix (dξuv ) in (1) are de-
pendent multivariate random variables normally distributed with known means
AN US
u d¯uv , ∀(u, v) ∈ A. Also, let (σ u ) = (σiv ), ∀(i, v) ∈ A, u ∈ V be defined as in the
previous section. Thus, a deterministic equivalent MILP model can be obtained 350
by means of the following proposition.
Proposition 2. Assume that there exists a sufficiently large number of line
M
segments l = {1, . . . , L} for fixed α ∈ [0, 0.5) and θ ≥ 1 in order to approximate 2 1/θ from above the univariate convex function F −1 (1 − α)zu on [0, 1], ∀u ∈
V , then the PCDT problem obtained by replacing the probabilistic constraint (1) with the deterministic joint probabilistic constraints (28)-(29) can be equivalently written as
min
{x,y,π,λ,µ,ϕ,φ,z}
CE AC
X
cuv xuv ,
(30)
(u,v)∈A
subject to (4) − (12), (20) − (23), (29) and X X X X u 2 (σiv ) λuv +
PT
M IP2
ED
355
i∈V
d2m +
+ v∈N(u) \{i}
X
+ + v∈N(u) k∈N(u) \{i,v}
(−2dm d¯uv + d¯2uv )xuv +
+ v∈N(u)
X
+ v∈N(u)
u u k σiv σik µuv ≤
X
d¯uv d¯uk ϕkuv ,
+ k∈N(u) \{v}
∀u ∈ V,
(31)
φu ≥ al + bl zu ,
∀u ∈ V, l = {1, . . . , L},
Mxuv + φu − M ≤ λuv ≤ Mxuv , λuv ≤ φu , λuv ≥ 0
∀u, v ∈ V, (u 6= v), ∀u, v ∈ V, (u 6= v),
∀u, v ∈ V, (u 6= v),
(32) (33) (34) (35)
ACCEPTED MANUSCRIPT
Mϕkuv + φu − M ≤ µkuv ≤ Mϕkuv , ∀u, v, k ∈ V, (u 6= v; u, v 6= k),
(36)
µkuv ≤ φu ,
(37)
∀u, v, k ∈ V, (u 6= v; u, v 6= k).
(38)
CR IP T
µkuv ≥ 0,
∀u, v, k ∈ V, (u 6= v; u, v 6= k),
Proof. From Lemma 1 and considering (20), we can write (28) as
i∈V
+ v∈N(u) \{i}
or equivalently
1/θ
(F −1 ((1 − α)zu ))2 d2m +
X
X
i∈V
X
(−2dm d¯uv + d¯2uv )xuv +
X
+ v∈N(u)
X
X
X
d¯uv xuv ,
+ + v∈N(u) k∈N(u) \{i,v}
X
2
u u σiv σik xuv xuk ≤
d¯uv d¯uk xuv xuk ,
+ + v∈N(u) k∈N(u) \{v}
(39)
As before, we define ϕkuv = xuv xuk , ∀u, v, k ∈ V , with u 6= v 6= k, and linearize 1/θ
these variables according to (21) and (22). Now, let φu ≥ (F −1 ((1 − α)zu ))2 =
ED
360
al + bl zu , ∀u ∈ V, l = {1, . . . , L} be a tight linear approximation of (F −1 ((1 − 1/θ
α)zu ))2 from above for L sufficiently large where each (al , bl ) denotes the coeffi-
PT
cients of each linear segment uniformly selected on the interval (0,1]. Notice that this interval is implied by constraints (29). Also, define nonnegative variables λuv = φu xuv , ∀ u, v ∈ V, (u 6= v) and µkuv = φu ϕkuv , ∀ u, v, k ∈ V, (u 6= v 6= k).
CE
365
The λ and µ variables can be linearized according to (33)-(35) and to (36)-(38), respectively, with M being a bigM positive constant. Then, after developing 1/θ
AC
(39) and approximating (F −1 ((1−α)zu ))2 by φu , ∀u ∈ V ; and replacing φu xuv
by λuv , ∀ u, v ∈ V, (u 6= v); and φu ϕkuv by µkuv , ∀ u, v, k ∈ V, (u 6= v 6= k), we
370
∀u ∈ V,
M
∀u ∈ V.
u σiv xuv ≤ dm −
u 2 (σiv ) xuv +
+ v∈N(u) \{i}
+ v∈N(u)
2
X
AN US
X 1/θ (F −1 ((1 − α)zu ))2
obtain the resulting constraints (31) under the condition that (32) and (29) must be verified in order to guaranty that the φu variables represent a very close ap1/θ
proximation of (F −1 ((1 − α)zu ))2 , ∀u ∈ V . Finally, the remaining constraints
ACCEPTED MANUSCRIPT
of this model ensure that feasible solutions are dominating arborescences, thus concluding the proof.
CR IP T
As a direct consequence of Proposition 2, we have the following corollary.
375
Corollary 1. The optimal objective function value of M IP2 is an upper bound for the optimal objective function value of M IP1 .
Proof. The proof is immediate due to the convexity of the function (F −1 ((1 − 1/θ
α)zu ))2 for values of zu ∈ [0, 1], α ∈ [0, 0.5) and θ ≥ 1 since (F −1 ((1 − 1/θ
α)zu ))2 ≥ (F −1 (1 − α))2 . This clearly implies that the feasible set of M IP2 is
AN US
380
contained in the feasible set of M IP1 .
Similarly as for M IP1 , we further consider the resulting MILP model obtained from M IP2 by replacing the constraints (8) with the set of constraints (13)-(15) and the constraints (10) with the constraints (16). Finally, we denote this new 385
model by M IP2V i and the corresponding LP relaxations of M IP2 and M IP2V i
M
by LP2 and LP2V i , respectively.
ED
5. Numerical results
In this section, we present numerical results for all the proposed models:
390
PT
DTc , DTcV i , LPc , LPcV i , M IP1 , M IP1V i , LP1 , LP1V i , M IP2 , M IP2V i , LP2 and LP2V i . We implement a Matlab program using CPLEX 12.6 for solving the MILP and LP models [14]. The numerical experiments have been carried
CE
out on an Intel(R) 64 bits core (TM) with 3.4 GH and 8 GB of RAM under Windows 7. We consider only connected disk graphs G = (V, E) where each
AC
disk represents the transmission range of a node v ∈ V . All instances considered
395
in this work are generated according to the procedure reported in [4, 26] which 2 is as follows: the cost of each edge {u, v} ∈ E is determined as cuv = wuv where
wuv represents the Euclidean distance between nodes u and v. An edge {u, v} between nodes u and v of V belongs to E if cuv is smaller than or equal to the considered node transmission range. In each instance, we distribute randomly
ACCEPTED MANUSCRIPT
400
|V | nodes in an area of 500 × 500 m2 and consider node transmission ranges equal to 100 m, 150 m and 350 m. The same transmission range applies for all nodes of a given instance. We further consider 18 benchmark DT instances
CR IP T
from [4, 26] with |V | = 50 and |V | = 100 with distinct transmission ranges of 100 m, 125 m and 150 m.
Concerning the uncertain parameters, each entry in matrix d¯ = (d¯uv ),
405
∀{u, v} ∈ E, is uniformly and randomly generated from the interval (0, 20].
u Similarly, we generate auxiliary matrices σau of entries (σiv )a , ∀u ∈ V , whose
values of their elements are uniformly and randomly drawn from the interval
410
AN US
[0, 1]. Each matrix σau , ∀u ∈ V , is multiplied by its transpose, thus originating
a new symmetric positive semidefinite matrix σ ˆ u := (σau ) × (σau )T . Finally, we
determine each matrix σ u from the one σ ˆ u by setting σ u := (ˆ σ u )1/2 , ∀u ∈ V
where (ˆ σ u )1/2 represents the square root of σ ˆ u . The values of α and dm are set to α = 0.1 and dm =
max d¯uv × γ, respectively, where γ ∈ [0.5, 3.2] is
{u,v}∈E
arbitrarily chosen. For the models M IP2 , M IP2V i , LP2 and LP2V i , we set θ = 1 1/θ
in the function (F −1 ((1 − α)zu ))2 . Notice that θ = 1 is the smallest value for
M
415
this parameter. It allows to obtain more conservative optimal solutions in the 1/θ
ED
corresponding models. When θ → ∞, then F −1 ((1 − α)zu ) → F −1 (1 − α) which corresponds to the case of separate probabilistic constraints that we use in our models M IP1 , M IP1V i , LP1 , and LP1V i . Considering values of θ in [1, ∞), we can obtain optimal solutions for both joint to separate probabilistic
PT
420
constraints in a unified framework. With respect to the line segments used in 1/θ
CE
the linear approximation of (F −1 ((1 − α)zu ))2 , we generate values of al and bl , ∀l ∈ {1, . . . , L}, by considering the interval (10−5 , 0.99995), for each zu ∈ V ,
with equally spaced subintervals of length 10−2 .
AC
425
In Figure 3, we plot the function (F −1 ((1 − α)z
1/θ
))2 for different values of
θ ∈ {1, 2, 10, 200}. From these curves, we confirm that this function is convex
and that for large values of θ, the curves become constant. In our numerical experiments, we observed that the length 10−2 is enough for these models to guarantee that they can find the optimal solution of each instance. Also, as the
ACCEPTED MANUSCRIPT
16 θ =1 θ =2 θ = 10 θ = 200
12
8 6 4 2 0
0
0.1
0.2
0.3
0.4
0.5
Values of z
0.8
0.9
1
Curves for the normal inverse cumulative distribution function for θ ∈
{1, 2, 10, 200}.
430
0.7
AN US
Figure 3:
0.6
CR IP T
10
1
1 22 2 1/θ F −1 (1 − α)z
14
range variation (see the vertical axis) of the curves in Figure 3 are not larger than 15 for values of z ∈ (0, 1), we set M = 20 for the joint models M IP2 , M IP2V i , LP2 and LP2V i .
M
We adopt the following organization for the tables presented hereafter. Each table reports results obtained with two models specified accordingly in the set of columns referred to M odelI and M odelII . In the last three rows of each table we
ED
435
present minimum, maximum and average (when all column data are available and, if not, we indicate as ‘na’) values for each column. Column (#) gives the
PT
instance identifier. The number of nodes and the density of edges in percentage (considering the maximum possible number of edges) of each instance are given by columns (|V |) and (Den %), respectively. Columns (Opt), (B&Bn),
CE
440
(T ime(s)), (LR) and (T imeLR (s)) of each ‘M odel∗ ’ give, respectively, the optimal solution value or the best incumbent solution value obtained by CPLEX
AC
in a time limit of 10 hours, the number of CPLEX branch and bound nodes, the CPU time in seconds to obtain the optimal or incumbent solution, the value
445
of the linear relaxation of the corresponding MILP formulation and the CPU time in seconds to obtain the linear relaxed solution. Finally, integrality gaps (in percentage) are reported in columns (GapI %) and (GapII %) referred to
ACCEPTED MANUSCRIPT
#
|V |
Den %
1
10
2
M odelII = DTcV i
M odelI = DTc
Integrality gaps
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
80.0
5954.64
0
0.18
994.06
0.20
5954.64
0
0.17
5954.64
0.16
83.31
0
12
80.3
2178.31
0
0.27
2178.31
0.25
2178.31
0
0.28
2178.31
0.25
0
0
3
14
81.3
11547.03
34
0.20
1056.38
0.16
11547.03
0
0.20
11547.03
0.19
90.85
0
4
15
58.1
29432.11
4
0.26
11550.21
0.17
29432.11
0
0.18
29432.11
0.19
60.76
0
5
16
84.2
1361.40
0
0.19
1361.40
0.19
1361.40
6
18
64.1
1665.67
0
0.20
838.55
0.19
1665.67
7
20
71.6
2041.29
0
0.19
1088.98
0.18
2041.29
8
22
87.9
62.34
0
0.20
62.34
0.17
62.34
9
24
65.9
18423.85
1194
0.53
3537.69
0.20
18423.85
10
25
76.3
2650.94
0
0.20
123.33
0.17
2650.94
11
30
76.8
3423.67
0
0.25
1124.68
0.18
3423.67
12
35
76.0
1755.49
0
0.23
916.62
0.21
1755.49
13
40
70.5
105.81
0
0.23
105.81
0.19
105.81
14
45
62.0
1292.91
0
0.29
443.64
0.21
1292.91
15
50
73.1
899.27
0
0.27
481.87
0.21
899.27
0
16
60
80.3
307.19
0
0.36
64.70
0.24
307.19
17
70
72.8
1473.30
508
1.61
231.18
0.26
18
80
78.4
482.77
3
0.66
156.02
19
90
76.3
5.63
0
0.59
5.63
20
100
69.1
2159.86
32990
14.95
21
110
71.5
192.58
0
1.03
22
120
79.0
330.56
0
1.20
Min.
58.10
5.63
0
0.18
Max.
87.90
29432.11
32990
Ave.
74.34
3988.48
1578.77
CR IP T
Opt
0.17
1361.40
0.20
0
0
0
0.19
1665.67
0.19
49.66
0
0
0.18
2041.29
0.17
46.65
0
0
0.20
62.34
0.19
0
0
0
0.20
13968.51
0.17
80.80
24.18
0
0.21
2650.94
0.20
95.35
0
0
0.25
3423.67
0.21
67.15
0
0
0.26
1755.49
0.22
47.79
0
0
0.31
105.81
0.23
0
0
0
0.30
1292.91
0.26
65.69
0
0.38
899.27
0.28
46.42
0
0
0.59
307.19
0.45
78.94
0
1473.30
0
0.81
876.69
0.48
84.31
40.49
0.34
482.77
0
1.10
482.77
0.64
67.68
0
0.42
5.63
0
1.45
5.63
0.79
0
0
49.50
0.40
2159.86
0
1.61
2159.86
0.90
97.71
0
5.50
0.54
192.58
0
2.10
192.58
1.21
97.14
0
27.90
0.65
330.56
0
3.12
330.56
1.71
91.55
0
5.50
0.16
5.63
0
0.17
5.63
0.16
0
0
14.95
11550.21
0.65
29432.11
0
3.12
29432.11
1.71
97.71
40.49
1.09
1200.19
0.26
3988.48
0
0.64
3758.84
0.42
56.89
2.93
M
AN US
0
ED
Table 1: Numerical results for models DTc and DTcV i using a transmission range of 350 m.
M odelI and M odelII , respectively. These gaps, for each model, are obtained according to the expression 100 Opt−LR Opt . In Tables 1-3, 4-6 and 7-9, we present numerical results for DTc , LPc , DTcV i , LPcV i , M IP1 , LP1 , M IP1V i , LP1V i , and
PT
450
M IP2 , LP2 , M IP2V i , LP2V i , respectively, for graph instances with transmission
CE
ranges of 350 m, 150 m and 100 m, respectively. The set of instances is the same in all these tables. This allows to compare the objective function values for all the proposed models with and without probabilistic constraints. In particular, it allows to compare all the models using separate and joint probabilistic con-
AC 455
straints.
We observe in Tables 1-3 that models DTc , DTcV i , M IP1 , M IP1V i
were able to solve to optimality instances with up to 120 nodes, while with model M IP2V i we were able to solve to optimality only instances with up to 60 nodes. This shows that solving the problem with joint probabilistic constraints
ACCEPTED MANUSCRIPT
M odelII = M IP1V i
M odelI = M IP1
Integrality gaps
|V |
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
1
10
11190.53
0
0.21
994.06
0.24
11190.53
0
0.27
5954.64
0.19
91.12
46.79
2
12
16008.16
0
0.36
2178.31
0.30
16008.16
0
0.30
2178.31
0.28
86.39
86.39
3
14
68487.36
0
0.25
1056.38
0.23
68487.36
0
0.28
11547.03
0.23
98.46
83.14
4
15
111723.24
0
0.37
11550.21
0.24
111723.24
0
0.36
29432.11
0.24
89.66
73.66
5
16
5850.25
0
0.33
1361.40
0.27
5850.25
0
6
18
29967.21
0
0.51
838.55
0.31
29967.21
0
7
20
101438.94
0
0.52
1088.98
0.37
101438.94
0
8
22
4061.08
0
0.73
62.34
0.47
4061.08
0
9
24
91649.72
37
1.28
3537.69
0.53
91649.72
0
10
25
5097.69
0
0.96
123.33
0.62
5097.69
0
11
30
81299.10
0
1.36
1124.68
1.09
81299.10
0
12
35
41773.18
0
2.28
916.62
1.59
41773.18
0
13
40
113464.73
0
3.44
124.56
2.38
113464.73
0
14
45
75417.12
0
6.63
443.64
3.17
75417.12
0
15
50
42249.56
0
7.33
481.87
16
60
47108.72
0
8.89
64.70
17
70
91211.93
0
13.39
231.18
18
80
57358.00
0
23.12
156.02
19
90
59042.83
0
36.55
5.63
20
100
37954.74
79749
135.12
49.50
21
110
55372.76
65
75.59
5.50
22
120
48777.20
20
113.13
27.90
Min.
4061.08
0
0.21
5.50
Max.
113464.73
79749
135.12
11550.21
Ave.
54386.54
3630.50
19.65
1361.40
0.27
76.73
76.73
0.42
1665.67
0.31
97.20
94.44
0.51
2041.29
0.38
98.93
97.99
0.78
62.34
0.47
98.46
98.46
1.22
14350.58
0.56
96.14
84.34
1.00
2650.94
0.65
97.58
48.00
1.62
3423.67
1.04
98.62
95.79
2.43
1760.55
1.65
97.81
95.79
3.62
730.10
2.40
99.89
99.36
6.50
1420.45
3.27
99.41
98.12
AN US
0.34
5.01
42249.56
0
8.13
899.27
5.01
98.86
97.87
9.98
47108.72
0
9.59
307.19
10.01
99.86
99.35
15.55
91211.93
0
14.07
1177.32
16.34
99.75
98.71
26.52
57358.00
0
24.28
531.59
26.52
99.73
99.07
40.54
59042.83
0
37.01
5.63
40.37
99.99
99.99
53.81
37954.74
0
85.95
2159.86
53.16
99.87
94.31
77.67
55372.76
0
78.90
192.58
77.15
99.99
99.65
121.83
48777.20
0
120.40
330.56
122.79
99.94
99.32
0.23
4061.08
0
0.27
5.63
0.19
76.73
46.79
121.83
113464.73
0
120.40
29432.11
122.79
99.99
99.99
16.48
54386.54
0
18.09
3826.50
16.51
96.56
89.42
M 1201.04
CR IP T
#
Table 2: Numerical results for models M IP1 and M IP1V i using a transmission range of 350
460
ED
m.
is significantly harder than solving it with separate probabilistic constraints.
PT
Note that, in general, optimal solution values for these instances follow the relations Opt(DTc ) = Opt(DTcV i ) ≤ Opt(M IP1 ) = Opt(M IP1V i ) ≤ Opt(M IP2 ) =
CE
Opt(M IP2V i ). Moreover, we observe that optimal solution values of M IP1 and
M IP2 are significantly larger than those of the classic problem DTc . In fact, in
465
these tables, optimal solution values of M IP1 and M IP2 are the same with the
AC
exception of the instance 4 which is larger for M IP2 . We also see, in terms of CPU times and number of CPLEX branch and bound nodes, that models DTcV i ,
M IP1V i and M IP2V i are more effective, presenting significantly smaller values than those obtained with models DTc , M IP1 and M IP2 , respectively. This can
470
be explained by the fact that, for the models DTcV i , M IP1V i and M IP2V i , we
ACCEPTED MANUSCRIPT
M odelII = M IP2V i
M odelI = M IP2
Integrality gaps
|V |
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
1
10
11190.53
0
0.33
994.06
0.21
11190.53
0
0.37
5954.64
0.25
91.12
46.79
2
12
16008.16
0
0.62
2178.31
0.27
16008.16
0
0.50
2178.31
0.28
86.39
86.39
3
14
68487.36
0
0.64
1056.38
0.28
68487.36
0
0.47
11547.03
0.27
98.46
83.14
4
15
165378.94
0
1.75
11550.21
0.30
165378.94
0
1.30
29432.11
0.29
93.02
82.20
5
16
5850.25
0
1.31
1361.40
0.34
5850.25
0
6
18
29967.21
0
2.28
838.55
0.41
29967.21
0
7
20
101438.94
0
3.04
1088.98
0.51
101438.94
0
8
22
4061.08
0
5.66
62.34
0.72
4061.08
0
9
24
91649.72
0
9.11
3537.69
0.77
91649.72
0
10
25
5097.69
0
9.48
123.33
0.92
5097.69
0
11
30
81299.10
0
22.71
1124.68
1.58
81299.10
0
12
35
41773.18
0
60.05
916.62
2.49
41773.18
0
13
40
113464.73
0
88.56
105.81
3.78
113464.73
0
14
45
75417.12
0
177.05
443.64
5.28
75417.12
275
15
50
42249.56
17816
23009.00
481.87
16
60
*
*
*
64.70
Min.
4061.08
0
0.33
62.34
Max.
165378.94
17816
23009
11550.21
Ave.
na
na
na
1620.53
1.29
1361.40
0.37
76.73
76.73
1.95
1665.67
0.42
97.20
94.44
2.76
2041.29
0.52
98.93
97.99
5.51
62.34
0.67
98.46
98.46
8.35
13968.51
0.81
96.14
84.76
9.02
2650.94
0.92
97.58
48.00
21.35
3423.67
1.55
98.62
95.79
59.32
1755.49
2.63
97.81
95.80
89.20
105.81
3.69
99.91
99.91
192.60
1292.91
5.46
99.41
98.29
AN US
‘*’: Instance not solved.
CR IP T
#
7.88
42249.56
2555
9071.78
899.27
11.04
98.85
97.87
14.71
47108.72
27539
13640.87
307.19
15.15
*
99.35
0.21
4061.08
0
0.37
62.34
0.25
76.73
46.79
14.71
165378.94
27539
13640.87
29432.11
15.15
99.91
99.91
2.52
56277.64
1898.06
1444.16
4915.41
2.77
na
86.61
Table 3: Numerical results for models M IP2 and M IP2V i using a transmission range of 350
M
m.
use generalized sub-tour elimination constraints as suggested in [11]. With this,
ED
the resulting linear relaxed solution values are considerably tighter than those of the models DTc , M IP1 and M IP2 . This is confirmed by the corresponding integrality gaps of these models. In Tables 4-6, we observe similar trends as for the Tables 1-3. We solve
PT
475
instances to optimality with up to 120 nodes with DTcV i and M IP1V i , while
CE
M IP2V i was able to solve instances only with up to 40 nodes. The average density of the instances with transmission range of 150 m in Tables 4-6 is 23.69%. In Tables 1-3, the average density is 74.34% for a transmission range of 350 m. In general, we observe that for sparser graphs in Tables 4-6, all the proposed
AC
480
models are harder to solve to optimality than for the instances in Tables 13. In particular, we do not solve the instances 13-22 with DTc and M IP1 in Tables 4-5 since the performance of these models deteriorate rapidly in terms of CPU time and number of CPLEX branch and bound nodes with respect to
ACCEPTED MANUSCRIPT
#
|V |
Den %
1
10
2
M odelII = DTcV i
M odelI = DTc
Integrality gaps
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
40.0
23092.52
0
0.16
18034.58
0.17
23092.52
0
0.18
20072.22
0.16
21.90
13.08
12
27.3
48404.68
1608
0.31
15473.25
0.17
48404.68
39
0.17
38899.51
0.16
68.03
19.64
3
14
31.9
49313.37
2996
0.81
12285.70
0.17
49313.37
0
0.19
44603.10
0.19
75.09
9.55
4
15
31.4
30541.11
1161
0.33
6222.58
0.18
30541.11
0
0.17
28662.80
0.16
79.63
6.15
5
16
23.3
55752.72
3322
1.01
10735.29
0.17
55752.72
6
18
21.6
81463.64
7774
2.17
18346.67
0.17
81463.64
7
20
32.1
24496.14
9014
1.93
3435.11
0.17
24496.14
8
22
18.2
84785.20
1640
0.47
29765.84
0.20
84785.20
9
24
19.2
80526.73
80956
19.13
27545.97
0.19
80526.73
10
25
21.0
81788.71
71415
18.38
11574.38
0.17
81788.71
11
30
23.4
68319.99
3609789
1497.94
14513.16
0.19
68319.99
12
35
20.8
62317.15
5407920
2673.93
4679.39
0.18
62317.15
13
40
23.2
*
*
*
8068.38
0.22
49488.68
14
45
17.7
*
*
*
4691.96
0.26
75328.01
15
50
21.6
*
*
*
2501.11
0.21
63502.41
16
60
21.4
*
*
*
2057.92
0.23
17
70
20.8
*
*
*
1104.35
18
80
21.3
*
*
*
19
90
20.3
*
*
20
100
22.3
*
21
110
21.2
22
120
CR IP T
Opt
0.19
51464.70
0.20
80.74
7.69
92
0.20
69614.28
0.17
77.48
14.55
0
0.18
20395.01
0.18
85.98
16.74
0
0.20
74023.19
0.19
64.89
12.69
43
0.25
64755.40
0.20
65.79
19.59
95
0.27
74403.74
0.18
85.85
9.03
673
1.28
52705.93
0.23
78.75
22.85
1237
0.84
50575.31
0.20
92.49
18.84
2145
1.32
42907.63
0.29
*
13.30
3315
4.59
57792.93
0.21
*
23.28
381682
571.11
41574.57
0.34
*
34.53
54085.65
490056
1272.71
37004.99
0.37
*
31.58
0.25
52190.98
147292
441.37
40132.70
0.45
*
23.10
1464.58
0.30
54560.08
20742
100.69
43144.40
0.52
*
20.92
*
1234.25
0.35
52221.99
1535848
12811.46
33703.59
0.64
*
35.46
*
*
472.70
0.34
39500.38
742084
13318.08
21983.80
4.58
*
44.35
*
*
*
838.56
0.44
46160.84
1068092
26345.97
29937.91
8.32
*
35.14
21.3
*
*
*
384.73
0.62
35746.04
97228
789.64
27000.90
0.87
*
24.46
Min.
17.70
23092.52
0
0.16
384.73
0.17
23092.52
0
0.17
20072.22
0.16
21.90
6.15
Max.
40
84785.20
5407920
2673.93
29765.84
0.62
84785.20
1535848
26345.97
74403.74
8.32
92.49
44.35
Ave.
23.69
na
na
na
8883.20
0.24
55163.04
204121.04
2530.04
43879.93
0.85
na
20.75
M
‘*’: Instance not solved.
AN US
0
485
ED
Table 4: Numerical results for models DTc and DTcV i using a transmission range of 150 m.
models DTcV i and M IP1V i . Similarly, we do not solve the instances 11-14 in
PT
Table 6 with M IP2 . Notice that the instance 14 in Table 6 cannot be solved to optimality in 10 hours. If we compare the average CPU time and number
CE
of CPLEX branch and bound (B&B) nodes required by DTcV i for the instances in Tables 1 and 4, we see an increment of 2530 seconds and of 204121 B&B
490
nodes, respectively. Similarly, if we compare the average CPU time and number
AC
of B&B nodes required by M IP1V i for the instances in Table 2 and 5, we see an increment of 117 seconds and of 313 B&B nodes, respectively. For M IP2V i ,
these values are also larger in Table 6 than in Table 3. In particular, we observe in Tables 4 and 5 that solving M IP1V i requires less CPU time and B&B nodes
495
than for DTcV i . With respect to the optimal solution values for these instances,
ACCEPTED MANUSCRIPT
M odelII = M IP1V i
M odelI = M IP1
Integrality gaps
|V |
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
1
10
23092.52
0
0.25
18034.58
0.18
23092.52
0
0.23
20072.22
0.18
21.90
13.08
2
12
49644.58
228
0.81
15473.25
0.19
49644.58
0
0.45
38899.51
0.19
68.83
21.64
3
14
49313.37
0
0.58
12285.70
0.20
49313.37
0
0.44
44603.10
0.23
75.09
9.55
4
15
37742.02
609
1.64
6222.58
0.23
37742.02
0
0.45
30479.76
0.24
83.51
19.24
5
16
65018.51
2670
3.76
10735.29
0.27
65018.51
0
6
18
105718.12
1846
4.41
18409.93
0.31
105718.12
0
7
20
28618.15
4354
13.29
3435.11
0.34
28618.15
0
8
22
90372.38
0
2.31
29804.10
0.41
90372.38
0
9
24
93457.94
8954
68.55
27565.03
0.45
93457.94
0
10
25
109801.97
12591
91.68
11574.38
0.59
109801.97
99
11
30
72915.68
298589
3098.10
14513.16
0.81
72915.68
211
12
35
65624.59
1954271
36000
4679.39
1.17
65624.59
0
13
40
*
*
*
8068.38
2.52
54640.40
51
14
45
*
*
*
4741.48
2.37
86892.79
1887
15
50
*
*
*
2501.11
16
60
*
*
*
2082.34
17
70
*
*
*
1104.35
18
80
*
*
*
1467.46
19
90
*
*
*
1247.48
20
100
*
*
*
472.70
21
110
*
*
*
844.37
22
120
*
*
*
386.83
Min.
23092.52
0
0.25
386.83
Max.
109801.97
1954271
36000
29804.10
Ave.
na
na
na
8893.13
51723.09
0.28
83.49
20.45
0.70
69614.28
0.30
82.59
34.15
1.10
20395.01
0.35
88.00
28.73
1.25
74023.19
0.40
67.02
18.09
2.23
65052.04
0.45
70.51
30.39
7.68
74414.69
0.53
89.46
32.23
13.31
53238.70
0.84
80.09
26.99
12.51
51020.42
1.33
92.86
22.25
25.73
43270.73
1.85
*
20.81
181.61
58056.61
2.80
*
33.19
AN US
0.55
3.27
74332.33
144
71.71
42282.05
4.14
*
43.12
5.92
79566.37
1424
246.78
37115.07
10.13
*
53.35
9.12
74225.29
1026
370.96
40653.33
24.82
*
45.23
13.75
96169.10
0
238.23
44211.81
51.48
*
54.03
19.78
119041.25
0
279.12
35418.41
67.32
*
70.25
28.92
105560.31
0
326.00
23133.77
176.51
*
78.08
40.62
110563.48
2044
779.15
30730.80
203.25
*
72.21
50.47
143280.95
21
415.24
27496.97
497.72
*
80.81
0.18
23092.52
0
0.23
20072.22
0.18
21.90
9.55
50.47
143280.95
2044
779.15
74414.69
497.72
92.86
80.81
8.26
78890.55
313.95
135.24
44359.34
47.51
na
37.63
M
‘*’: Instance not solved.
CR IP T
#
ED
Table 5: Numerical results for models M IP1 and M IP1V i using a transmission range of 150 m.
PT
we observe again the relations Opt(DTcV i ) ≤ Opt(M IP1V i ) ≤ Opt(M IP2V i ). However, in this case, we observe that the optimal solution values obtained with M IP2V i are larger than those obtained with M IP1V i for the most part
CE
of the instances in Tables 5 and 6. As before, all the proposed models which
500
use generalized sub-tour elimination constraints present smaller CPU times and
AC
B&B nodes than the other models, as well as tighter linear relaxed solution values. In Tables 7-9, we confirm the trends observed in Tables 1-6. In this case,
we solve instances with transmission range equal to 100 m and with up to
505
90, 120 and 30 nodes to optimality with models DTcV i , M IP1V i and M IP2V i , respectively. The average density for these instances is 15.23%. In this case, we
ACCEPTED MANUSCRIPT
M odelII = M IP2V i
M odelI = M IP2
Integrality gaps
|V |
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
1
10
31209.13
19
0.59
18034.58
0.22
31209.13
0
0.72
20072.22
0.21
42.21
35.68
2
12
49644.58
292
1.20
15473.25
0.23
49644.58
0
0.81
38899.51
0.22
68.83
21.64
3
14
58875.71
10
2.56
12285.70
0.28
58875.71
0
1.51
44603.10
0.27
79.13
24.24
4
15
50636.68
1637
7.57
6222.58
0.30
50636.68
17
2.20
30479.76
0.29
87.71
39.81
5
16
66653.52
1086
10.37
10735.29
0.31
66653.52
0
6
18
109854.11
1128
19.88
18350.30
0.37
109854.11
0
7
20
28618.15
2124
40.34
3435.11
0.48
28618.15
0
8
22
119731.92
4415
103.47
29765.84
0.58
119731.92
965
9
24
93457.94
37747
946.02
27545.97
0.80
93457.94
0
10
25
120619.93
3951
250.34
11574.38
1.00
120619.93
141
11
30
*
*
*
14513.16
1.32
74851.40
1037
12
35
*
*
*
4679.39
2.05
69707.63
1147
13
40
*
*
*
8068.38
3.06
68765.65
40751
14
45
*
*
*
4701.45
5.96
-
16585
Min.
28618.15
10
0.59
3435.11
Max.
120619.93
37747
946.02
29765.84
Ave.
na
na
na
13241.81
1.61
51493.66
0.34
83.89
22.74
2.71
69614.28
0.39
83.30
36.63
5.48
20395.01
0.49
88.00
28.73
27.71
74023.19
0.64
75.14
38.18
19.72
64899.85
0.75
70.53
30.56
28.54
74404.22
0.89
90.40
38.32
204.77
53035.09
1.68
*
29.15
324.69
50591.14
2.74
*
27.42
31901.30
43031.95
4.46
*
37.42
36000
57792.92
5.18
-
-
AN US
‘*’: Instance not solved.
CR IP T
#
0.22
28618.15
5.96
120619.93
1.21
na
0
0.72
20072.22
0.21
42.21
21.64
40751
36000
74404.22
5.18
90.40
39.81
4331.64
4894.41
49523.99
1.32
na
na
‘-’: No solution found by CPLEX in 10 hours of execution time.
Table 6: Numerical results for models M IP2 and M IP2V i using a transmission range of 150
M
m.
do not solve the instances 12-22 with models DTc and M IP1 in Tables 7-8 and do not solve the instances 12-14 in Table 9 neither with model M IP2 nor with model
510
ED
M IP2V i . For the latter, these instances proved to be infeasible which evidences the more restrictive nature of the joint probabilistic constraints. In particular,
PT
notice that for the instance 11 in Table 9, model M IP2 cannot obtain a feasible solution in 10 hours while M IP2V i allows to obtain the optimal solution of this instance in less than 32 seconds. This shows the remarkable effectiveness of using
CE
generalized sub-tour elimination constraints in M IP2 . The remaining instances
515
in Table 9 are solved to optimality in less than 13 seconds. The maximum CPU time and number of B&B nodes required by DTcV i for the instances in Table
AC
4 are 26345.97 seconds and 1535848 B&B nodes, respectively. While these values in Table 7 are 10 hours and 3283374 B&B nodes, respectively. Similarly, the maximum CPU time and number of B&B nodes required by M IP1V i for
520
the instances in Tables 5 and 8 are 779.15 seconds and 2044 B&B nodes, and 6842.07 seconds and 7316 B&B nodes, respectively. Once again, we see that
ACCEPTED MANUSCRIPT
#
|V |
Den %
1
10
2
M odelII = DTcV i
M odelI = DTc
Integrality gaps
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
24.4
37059.44
109
0.18
12988.63
0.17
37059.44
0
0.16
32670.26
0.16
64.95
11.84
12
24.2
29210.83
1743
0.34
6144.18
0.17
29210.83
0
0.19
25446.84
0.17
78.97
12.89
3
14
18.7
43698.31
859
0.28
22948.40
0.19
43698.31
0
0.17
36559.41
0.17
47.48
16.34
4
15
30.5
15058.40
132
0.26
6452.20
0.18
15058.40
0
0.19
14399.43
0.16
57.15
5
16
28.3
37183.35
176606
23.74
7988.78
0.17
37183.35
6
18
19.0
51756.80
37094
5.55
8945.88
0.16
51756.80
7
20
19.5
37129.01
19444
3.78
14016.64
0.18
37129.01
8
22
12.6
53246.64
1907
0.44
24492.05
0.17
53246.64
CR IP T
Opt
0.76
26066.52
0.17
78.52
29.90 8.77
0
0.17
47217.18
0.17
82.72
53
0.20
33568.69
0.18
62.25
9.59
5
0.20
42900.75
0.19
54.00
19.43
140
0.27
68481.14
0.19
69.78
13.12
125
0.23
51847.69
0.18
68.20
17.13
0
0.34
56431.89
0.18
79.51
14.63
3506
1.75
69276.38
0.24
*
26.10
4
0.62
58193.85
0.42
*
14.81
1722
0.97
72773.23
0.20
*
13.46
24
13.4
78826.94
228457
43.12
23817.61
0.17
78826.94
25
13.3
62563.97
12894
3.02
19895.33
0.17
62563.97
11
30
11.7
66102.25
141751
40.52
13543.69
0.20
66102.25
12
35
11.4
*
*
*
14587.62
0.18
93746.39
13
40
13.1
*
*
*
11592.88
0.33
68309.31
14
45
9.9
*
*
*
17672.45
0.19
84094.27
15
50
13.5
*
*
*
9096.80
0.20
61561.21
15843
14.48
46175.67
0.22
*
24.99
16
60
10.1
*
*
*
8056.37
0.26
81849.38
1842122
1529.41
59355.98
0.36
*
27.48
17
70
9.9
*
*
*
9542.68
0.25
90175.30
804550
1065.83
65187.84
0.39
*
27.71
18
80
10.7
*
*
*
7484.63
0.25
69333.96
3283374
36000
50158.65
0.70
*
27.66
19
90
10.3
*
*
*
5790.72
0.30
65868.76
865667
4886.80
53093.35
0.55
*
19.40
20
100
10.1
*
*
*
5111.89
0.33
77769.10
2565329
36000
61545.34
1.26
*
20.86
21
110
10.3
*
*
*
3552.84
0.41
60757.72
2193242
36000
45797.55
1.33
*
24.62
22
120
10.2
*
*
*
3136.19
1.01
-
-
-
41738.95
2.87
*
-
Min.
9.90
15058.40
109.00
0.18
3136.19
0.16
15058.40
0
0.16
14399.43
0.16
47.48
4.38
Max.
30.50
78826.94
228457
43.12
24492.05
1.01
93746.39
3283374
36000
72773.23
2.87
82.72
29.90
Ave.
15.23
na
na
na
11675.38
0.25
na
na
na
48131.20
0.47
na
na
M
AN US
9 10
‘*’: Instance not solved.
ED
‘-’: No solution found by CPLEX due to a shortage of memory in at most 10 hours of execution time.
Table 7: Numerical results for models DTc and DTcV i using a transmission range of 100 m.
PT
solving M IP1V i is easier than solving DTcV i . Observe that the optimal solution values obtained with M IP2V i are larger than those obtained with M IP1V i for many instances in Tables 8 and 9. In order to give more insight with respect to the parameter values of θ ∈ [1, ∞), in Figure 4 we plot curves for the objective
CE 525
function values of M IP2V i and LP2V i for values of θ ∈ [1, 30]. We choose this
AC
interval as we noticed in our numerical results that for large values of θ, these curves remain constant. Consequently, it would suffice to solve M IP1V i instead of M IP2V i for larger values of θ as they are equivalent models. By doing so, we
530
4.38
3137
cover all range of optimal solution values for θ ∈ [1, ∞). In other words, we cover joint and separate probabilistic constraints in a unified framework. In Figure 4, we randomly generate 20 instances with dimensions of |V | ∈ {20, 30} while using
ACCEPTED MANUSCRIPT
M odelII = M IP1V i
M odelI = M IP1
Integrality gaps
|V |
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
1
10
37135.50
0
0.24
12988.63
0.19
37135.50
0
0.26
32670.26
0.19
65.02
12.02
2
12
29210.83
311
0.62
6144.18
0.20
29210.83
0
0.31
25446.84
0.20
78.97
12.89
3
14
46506.67
0
0.51
23066.16
0.22
46506.67
0
0.39
37640.09
0.23
50.40
19.07
4
15
17534.73
1366
2.89
6452.20
0.24
17534.73
0
0.47
14399.43
0.24
63.20
17.88
5
16
37371.33
6155
11.37
7988.78
0.27
37371.33
206
6
18
53950.88
350
3.35
8945.88
0.30
53950.88
0
7
20
39473.32
6800
22.24
14016.64
0.42
39473.32
0
8
22
54899.55
21
2.56
24531.80
0.38
54899.55
0
9
24
81104.56
28899
135.35
23817.61
0.51
81104.56
11
10
25
62563.97
3374
22.88
19895.33
0.50
62563.97
0
11
30
84006.04
2935324
36000
13543.69
0.85
84006.04
0
12
35
*
*
*
14587.62
1.20
102660.28
7316
13
40
*
*
*
11592.92
1.86
83699.46
808
14
45
*
*
*
17787.13
2.59
100856.63
0
15
50
*
*
*
9096.80
16
60
*
*
*
8077.72
17
70
*
*
*
9589.57
18
80
*
*
*
7511.77
19
90
*
*
*
5805.59
20
100
*
*
*
5177.71
21
110
*
*
*
3594.44
22
120
*
*
*
3148.52
Min.
17534.73
0
0.24
3148.52
Max.
84006.04
2935324
36000
24531.80
Ave.
na
na
na
26106.59
0.25
78.62
30.14
0.62
47469.20
0.28
83.42
12.01
0.99
33783.35
0.34
64.49
14.41
1.05
42900.75
0.41
55.32
21.86
3.96
68568.85
0.45
70.63
15.46
2.50
51849.05
0.52
68.20
17.13
3.83
56750.11
0.88
83.87
32.45
186.99
70581.10
1.35
*
31.25
49.44
59177.87
2.03
*
29.30
25.85
72947.29
2.54
*
27.67
AN US
1.12
3.18
82085.96
1089
211.40
47330.98
3.70
*
42.34
5.44
94646.26
4341
1155.25
60609.40
8.11
*
35.96
8.41
113202.67
2628
1308.36
68833.86
11.92
*
39.19
12.59
115524.99
63
513.74
51738.12
40.11
*
55.21
18.02
97746.86
4828
4277.56
54613.25
59.67
*
44.13
24.73
135766.52
1243
1274.63
62456.70
82.06
*
54.00
34.09
84259.69
3550
5918.51
46655.89
157.72
*
44.63
44.66
119676.87
5766
6842.07
45150.42
189.31
*
62.27
0.19
17534.73
0
0.26
14399.43
0.19
50.40
12.01
44.66
135766.52
7316
6842.07
72947.29
189.31
83.87
62.27
7.31
76085.61
1447.68
989.96
48985.42
25.56
na
30.51
M
‘*’: Instance not solved.
11698.21
CR IP T
#
ED
Table 8: Numerical results for models M IP1 and M IP1V i using a transmission range of 100 m.
PT
transmission ranges equal to 100 m and 150 m. Thus, we compute the average optimal solution values of the MILP model and its linear relaxation. We do not 535
consider large size instances because the trends of the curves remain the same
CE
and also because the CPU times become rapidly prohibitive. From Figure 4, we observe that the objective function values of M IP2V i are considerably larger
AC
for values of θ ∈ [1, 10]. Finally, we observe that the linear relaxation remains constant for all these instances.
540
In Table 10, we present numerical results for 18 benchmark instances [4, 26]
with dimensions of |V | ∈ {50, 100} and transmission ranges equal to 150 m, 125
m and 100 m. We report results only for DTcV i and M IP1V i . We do not present
numerical results for M IP2V i for these instances because CPLEX runs out of
ACCEPTED MANUSCRIPT
M odelII = M IP2V i
M odelI = M IP2
Integrality gaps
|V |
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
1
10
37135.50
0
0.42
12988.63
0.20
37135.50
0
0.39
32670.26
0.21
65.02
12.02
2
12
30341.18
237
1.28
6144.18
0.23
30341.18
0
0.62
25446.84
0.23
79.75
16.13
3
14
46506.67
0
1.01
23066.16
0.28
46506.67
0
0.83
37640.09
0.25
50.40
19.07
4
15
17534.73
96
2.80
6452.20
0.29
17534.73
0
1.51
14399.43
0.30
63.20
17.88
5
16
37371.33
4393
27.08
7988.78
0.34
37371.33
157
5.90
26071.30
0.33
78.62
30.24
6
18
56179.84
1749
16.85
8945.88
0.41
56179.84
3
4.34
47469.20
0.39
84.08
15.50
7
20
41537.16
7481
110.38
14016.64
0.47
41537.16
3
7.91
33777.43
0.65
66.26
18.68
8
22
54899.55
111
10.58
24531.80
0.62
54899.55
0
2.62
42900.75
0.56
55.32
21.86
9
24
81104.56
31519
450.44
23817.61
0.75
81104.56
43
10.78
68499.04
0.73
70.63
15.54
10
25
66078.65
2969
98.75
19895.33
0.97
66078.65
38
12.15
52077.55
0.90
69.89
21.19
11
30
-
447732
36000
13543.69
1.31
84090.22
5
31.40
56538.15
1.43
-
32.76
12
35
Inf
Inf
Inf
12988.62
0.33
Inf
Inf
Inf
32670.26
0.23
Inf
Inf
13
40
Inf
Inf
Inf
11630.40
3.03
Inf
Inf
Inf
58843.38
3.77
Inf
Inf
14
45
Inf
Inf
Inf
17709.88
4.13
Inf
Inf
Inf
72884.68
5.32
Inf
Inf
Min.
17534.73
0
0.42
6144.18
Max.
81104.56
447732
36000
24531.80
Ave.
na
na
na
14551.41
‘Inf’: Infeasible solution.
AN US
‘-’: No solution found with CPLEX in 10 hours.
CR IP T
#
0.20
17534.73
0
0.39
14399.43
0.21
50.40
12.02
4.13
84090.22
157
31.40
72884.68
5.32
84.08
32.76
0.95
na
na
na
42992.02
1.09
na
na
Table 9: Numerical results for models M IP2 and M IP2V i using a transmission range of 100
4
M
x 10
6
MIP2V i
ED
Optimal solutions
7
LP2V i
5
4
10 15 20 25 Values of θ with |V|=20 and using a transmission range of 100m .
9
LP2V i 7
Optimal solutions
Optimal solutions
CE 5.5
AC
6
7
6
5
5
10 15 20 25 Values of θ with |V|=20 and using a transmission range of 150m .
30
4
x 10
6.5
MIP2V i
8
4
7.5
x 10
6
30
PT
5
4
10
Optimal solutions
m.
MIP2V i LP2V i
10 15 20 25 Values of θ with |V|=30 and using a transmission range of 100m .
30
x 10
5.5 MIP2V i LP2V i
5
4.5
5
10 15 20 25 Values of θ with |V|=30 and using a transmission range of 150m .
Figure 4: Objective function values for M IP2V i and LP2V i with θ ∈ [1, 30].
30
ACCEPTED MANUSCRIPT
#
Den. %
M odelI = DTcV i Opt
B&Bn
T ime(s)
LR
M odelII = M IP1V i T imeLR (s)
Opt
B&Bn
T ime(s)
LR
Integrality gaps T imeLR (s)
GapI %
GapII %
ins-050-1
22.5
647.75
930
2.36
569.33
0.29
694.59
0
57.03
569.68
4.72
12.11
17.98
ins-050-2
21.4
863.69
2147
7.32
751.17
0.29
907.24
43
66.68
751.17
4.94
13.03
17.20
ins-050-3
22.0
743.94†
819
2.84
659.58
0.29
793.46
9
60.44
659.58
4.21
11.34
16.87
ins-100-1
21.0
876.69
148438
2085.21
755.25
0.78
1020.90
17
610.92
755.25
140.19
13.85
26.02
ins-100-2
21.8
657.35
36441
681.26
576.02
0.79
745.14
0
545.82
576.19
124.98
12.37
22.67
ins-100-3
20.5
722.87
3057
98.61
636.80
0.68
884.43
CR IP T
Benchmark instances with |V | = 50 and |V | = 100 and transmission range of 150 m.
82
689.12
636.83
135.62
11.91
27.99
Benchmark instances with |V | = 50 and |V | = 100 and transmission range of 125 m. 15.8
802.95
2677
5.41
696.11
0.26
835.45
135
69.36
696.13
4.00
13.31
16.68
ins-050-2
15.7
1055.10
101943
136.01
899.01
0.33
1122.71
3356
601.42
899.16
4.53
14.79
19.91
ins-050-3
15.3
877.77
13144
23.36
764.52
0.34
914.92
411
99.06
764.65
4.75
12.90
16.42
ins-100-1
16.1
943.01
92870
1087.90
839.79
0.82
1018.00
61
1214.40
839.79
129.58
10.95
17.51
ins-100-2
15.8
917.00
1229927
22802.25
754.12
1.06
963.19
198
1261.36
754.12
69.93
17.76
21.71
ins-100-3
14.7
998.18
542634
9003.54
839.74
0.66
1077.58
413
1760.61
839.74
111.58
15.87
22.07
AN US
ins-050-1
Benchmark instances with |V | = 50 and |V | = 100 and transmission range of 100 m. ins-050-1
10.0
1204.41
2985
3.93
1101.23
0.36
1300.09
6181
591.84
1101.23
3.75
8.57
15.30
ins-050-2
9.6
1340.44
249177
136.92
1099.81
0.40
1386.78
18091
1464.24
1101.46
3.82
17.95
20.57
ins-050-3
10.3
1316.39
374158
257.42
1131.34
0.38
1426.21
8610
834.99
1131.62
4.89
14.06
20.66
ins-100-1
10.8
1217.47
2020217
25155.92
1008.94
2.28
1256.97
3152
11689.91
1008.94
82.98
17.13
19.73
ins-100-2
10.6
1128.40
795701
6789.46
945.70
0.72
1204.59
11402
36000
948.02
50.17
16.19
21.29
ins-100-3
9.7
1252.99
1594116
13866.27
1055.37
0.72
1358.58
3355
36000
1055.64
78.13
15.77
22.30
‘†’: In [26] the reported optimal solution value is 743.74. The correct is really 743.94.
M
Table 10: Numerical results for Benchmark instances [26].
memory rapidly. From Table 10, we confirm that our proposed deterministic compact model DTcV i can solve all these instances to optimality in less than
ED
545
7 hours of CPU time. Consequently, we further confirm that the solutions
PT
obtained with the meta-heuristic approaches proposed by the authors in [4, 26] are the optimal ones (with the exception of instance ins-050-3 with transmission range of 150 m). These values were not proven optimal in [26]. In particular, we observe that the CPU time required by CPLEX to solve DTcV i increases when
CE 550
the density of the graphs decreases. Regarding the optimal solution values
AC
of M IP1V i , we observe that they are larger than those obtained with DTcV i . Conversely, we observe that the CPU times are lower except for the instances ins100-2 and ins-100-2 with transmission range of 100 m that CPLEX cannot solve
555
with M IP1V i in 10 hours of execution time. Again, we observe that the CPU times increase when the density of the graphs is smaller. Finally, we observe that the CPU time to solve the linear relaxations is significantly larger for M IP1V i
ACCEPTED MANUSCRIPT
than for DTcV i , whereas the number of B&B nodes is significantly larger for the latter model. The integrality gaps confirm that the linear relaxation of DTcV i is tighter than the one of M IP1V i .
CR IP T
560
6. Conclusion
In this paper, we introduce a new probabilistically constrained dominating
tree problem which consists in finding a minimum cost dominating tree T of
an undirected sensor network G subject to a probabilistic constraint related 565
to the risk that the total power consumption of each sensor in T exceeds a
AN US
given limit. To handle uncertain sensor power consumption of the links of the network, we consider both separate and joint probabilistic constraints in an im-
proved compact polynomial formulation for the classic dominating tree problem. The dependence among power sensor consumption of the links of the network is 570
handled by means of copulas [6, 13]. The compact formulation is based on ar-
M
borescences of directed graphs for which we give a proof of its correctness. The compact model allows to obtain deterministic equivalent mixed integer linear programming formulations for the probabilistic counterpart. Consequently, we
575
ED
solve random and challenging benchmark instances for the DT and PCDT problems with up to 120 nodes to optimality. Solving these problems to optimality
PT
is extremely important as it allows to compare the efficiency of any other algorithmic procedure to be proposed as part of future research. Additionally, new modeling approaches while considering the overhead and other metrics of sensor
CE
networks should be considered as variants of the DT and PCDT problems in
580
order to generate new protocols. Finally, it is important to note that the probabilistic modelling approach allows to handle joint and separate probabilistic
AC
constraints in a unified framework.
Acknowledgments The authors are grateful to the five referees and to the Editor Professor
585
Stavros Toumpis for their valuable comments and suggestions that helped to
ACCEPTED MANUSCRIPT
improve significantly our paper. The first author also acknowledges the financial support of the USACH/DICYT Project 061513VC DAS.
References 590
CR IP T
References
[1] P. Adasme and A. Lisser, Stochastic and semidefinite optimization for
scheduling in orthogonal frequency division multiple access networks, Journal of Scheduling 17(5), 445-469, (2014).
AN US
[2] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004. 595
[3] A. Charnes, W. Cooper and G. Symonds, Cost horizons and certainty equivalents: an approach to stochastic programming of heating oil, Management Science, vol. 4, 235-263, (1958).
M
[4] S. N. Chaurasia and A. Singh, A hybrid heuristic for dominating tree problem, Soft Computing, 20 (1), 377-397, (2016). [5] J. Cheng, E. Delage and A. Lisser, Distributionally Robust Stochas-
ED
600
tic Knapsack Problem, SIAM Journal on Optimization 24(3), 1485-1506
PT
(2014).
[6] J. Cheng, M. Houda and A. Lisser, Chance constrained 0–1 quadratic pro-
CE
grams using copulas, Optimization Letters 9(7), 1283-1295, (2015). 605
[7] J. Cheng and A. Lisser, A second-order cone programming approach for
AC
linear programs with joint probabilistic constraints, Operations Research Letters, vol. 40, 325-328, (2012).
[8] J. Carle and D. Simplot-Ryl, Energy-efficient area monitoring for sensor
610
networks, Computer, 37(2), 40-46, (2004).
[9] R. Fortet, Applications de l’alg`ebre de boole en recherche operationelle, Revue Francaise de Recherche Operationelle, 4, 17-26, (1960).
ACCEPTED MANUSCRIPT
[10] A. Gaivoronski, A. Lisser, R. Lopez and X. Hu, Knapsack problem with probability constraints, Journal of Global Optimization, 49, 397-413, (2011). [11] B. Gendron, A. Lucena, A. Salles da Cunha and L. Simonetti, Benders
CR IP T
615
Decomposition, Branch-and-Cut, and Hybrid Algorithms for the Minimum
Connected Dominating Set Problem, INFORMS Journal on Computing 26(4), 645-657, (2014).
[12] A. Gupte, S. Ahmed, M. Cheon and S. Dey, Solving Mixed Integer Bilinear
AN US
Problems Using MILP Formulations, SIAM Journal on Optimization 23(2),
620
721-744, (2013).
[13] R. Henrion and C. Strugarek, Convexity of chance constraints with dependent random variables: the use of copulae. In: Bertocchi M., Consigli G., Dempster M.A.H. (eds.), Chapter: Stochastic Optimization Methods in Finance and Energy Vol. 163 of the series International Series in Operations
625
M
Research & Management Science, 427-439, Springer, New York (2011).
ED
[14] IBM ILOG CPLEX Optimization Studio Information Center. [15] M. Lobo, L. Vandenberghe, S. Boyd and H. Lebret, Linear Algebra and its Applications, 284, 193-228, (1998). [16] J. Luedtke, S. Ahmed and G. Nemhauser, An integer programming ap-
PT
630
proach for linear programs with probabilistic constraints, Mathematical
CE
Programming, 122, 247-272, (2010). [17] L. Miller and H. Wagner, Chance-constrained programming with joint con-
AC
straints, Operations Research, 13, 930-945, (1965).
635
[18] R. Nelsen, An Introduction to Copulas, 2nd edn. Springer, New York, (2006).
[19] A. Pr´ekopa, On probabilistic constrained programming, Proceedings of the Princeton Symposium on Mathematical Programming, Princeton University Press, Princeton, (1970).
ACCEPTED MANUSCRIPT
640
[20] A. Pr´ekopa, Stochastic Programming, 600. Kluwer Academic Publishers, Dordrecht Boston, (1995). [21] N. Sahinidis, Optimization under uncertainty: State-of-the-art and oppor-
CR IP T
tunities, Computers and Chemical Engineering, 28, 971-983, (2004).
[22] R. Schultz, L. Stougie and M. van der Vlerk, Two-stage stochastic integer programming: A survey, Statistica Neerlandica, 50, 404-416, (1996).
645
[23] A. Shapiro, D. Dentcheva and A. Ruszczynski, Lectures on stochastic programming: Modeling and theory, 436. SIAM Philadelphia, Series on Opti-
AN US
mization, Vol. 9 of MPS/SIAM, Philadelphia, (2009).
[24] I. Shin, Y. Shen and M. Thai, On approximation of dominating tree in wireless sensor networks, Optimization Letters 4, 393-403, (2010).
650
[25] V. Singh, S. Jain and A. Tyagi, Risk and reliability analysis: a handbook
M
for civil and environmental engineers, vol. 1. ASCE, (2007). [26] S. Sundar and A. Singh, New heuristic approaches for the dominating tree
AC
CE
PT
ED
problem, Applied Soft Computing 13, 4695-4703, (2013).
ACCEPTED MANUSCRIPT
Pablo Adasme is an associate researcher and full professor in computer sci-
655
ence at the Electrical Engineering Department of the Universidad de Santiago de Chile. He received the PhD degrees from the Universidad de Santiago de
CR IP T
Chile and from the University of Paris Sud 11, in 2010 in Chile and France, respectively. Currently, his main research interests are stochastic and semidefi660
nite programming, combinatorial optimization problems, and metaheuristic al-
PT
ED
M
AN US
gorithms applied to wireless communication networks and energy problems.
Rafael Andrade received the B.Sc. degree in Computer Science from State University of Cear (UECE, Fortaleza, Brazil) in 1997 and in 1999, the M.Sc. degree in Computer Science and System Engineering from Federal University of
CE 665
Rio de Janeiro (UFRJ, Rio de Janeiro, Brazil). From 1999 to 2002 he realized his
AC
Ph.D. thesis at France Tlcom R&D laboratory (FTR&D, Issy-les-Moulineaux, France). In 2002, he obtained his Ph.D. degree in Informatics from Paris-Nord University (Paris 13, Villetaneuse, France). He was assistant professor at Paris-
670
Sud University and at Paris-Nord University, in the Laboratory for Computer Science (LRI) from 2002 to 2003, and in the Computer Science Laboratory of Paris-Nord (LIPN) from 2003 to 2004, respectively. In the occasion of a sabbat-
ACCEPTED MANUSCRIPT
ical year in 2012, he was invited professor at LRI (Paris-Sud). Currently, he is Associate Professor at Federal University of Cear (UFC, Fortaleza, Brazil) in the 675
department of statistics and applied mathematics where he has been working
CR IP T
since 2004. His main research interests are stochastic, robust and deterministic programming, operations research, combinatorial optimization, telecommunica-
ED
M
AN US
tion network design and network optimization problems.
Abdel Lisser is a full Professor in Computer Science laboratory of University
PT
680
of Paris Sud since 2001. He was heading research group at France Telecom from Research Center 1996 to 2001 and research engineer at France Telecom
CE
Research Center from 1988 to 1996. He got the Ph.D. at the University of Paris Dauphine in Computer Science in 1987 and the Habilitation thesis at the University of Paris Nord in 2000. His main research area is combinatorial
AC
685
and stochastic optimization with application to telecommunication and recently energy problems.
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Minimum cost dominating tree sensor networks under probabilistic constraints Pablo Adasme, Rafael Andrade, Abdel Lisser PII: DOI: Reference:
S1389-1286(16)30387-5 10.1016/j.comnet.2016.11.011 COMPNW 6053
To appear in:
Computer Networks
Received date: Revised date: Accepted date:
6 May 2016 9 November 2016 12 November 2016
Please cite this article as: Pablo Adasme, Rafael Andrade, Abdel Lisser, Minimum cost dominating tree sensor networks under probabilistic constraints, Computer Networks (2016), doi: 10.1016/j.comnet.2016.11.011
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Minimum cost dominating tree sensor networks under probabilistic constraints
CR IP T
Pablo Adasmea , Rafael Andradeb , Abdel Lisserc a Departamento
AN US
de Ingenier´ıa El´ ectrica, Universidad de Santiago de Chile, Avenida Ecuador 3519, Santiago, Chile. [email protected] b Departamento de Estat´ ıstica e Matem´ atica Aplicada, Universidade Federal do Cear´ a, Campus do Pici, BL 910, CEP 60.455-760 - Fortaleza, Cear´ a, Brasil. [email protected] c Laboratoire de Recherche en Informatique, Universit´ e de Paris-Sud 11, Ada Lovelace, 91405 Paris, France. [email protected]
Abstract
There is an increasing interest in planning sensor networks by considering both the impact of distances among sensors and the risk that power consumption leads to a very small network lifetime. A sensor failure can affect sensors in its neighborhood and compromise the network data communication. Weather
M
conditions may cause the power consumption of data communication to vary with uncertainty. This work introduces a compact probabilistic optimization
ED
approach to handle this problem while considering jointly or separately dependence among power consumption of the links of the network in a unified framework. We explore the concept of copulas in a dominating arborescence
PT
(DA) model for directed graphs, extended accordingly to handle the uncertain parameters. We give a proof of the DA model correctness and show that it can
CE
solve to optimality some benchmark instances of the deterministic dominating tree problem. Numerical results for the probabilistic approach show that our model tackles randomly generated instances with up to 120 nodes.
AC
Keywords: Dominating tree, probabilistic constraints, compact formulation, mixed integer linear programming, sensor network design
Preprint submitted to Journal of Computer Networks
November 14, 2016
ACCEPTED MANUSCRIPT
1. Introduction A classic version of the dominating tree (DT) sensor network problem [24] concerns determining a minimum cost tree T in an undirected graph G = (V, E)
5
CR IP T
of set of nodes V and set of edges E, with ce ≥ 0 being the cost (weight-distance) of edge e ∈ E, where any node v ∈ V is in T or is adjacent to a node in this tree.
Link cost distances are related as usual [24] with the power required to establish
the topology of the involved dominating tree connections. The topology cost of
a tree T is the sum of the cost distances of its edges. The DT problem is known to be NP-Hard [24].
AN US
In some sensor networks, with each edge {u, v} ∈ E, of extremities u, v ∈ V ,
10
there is associated an uncertain power consumption dξuv necessary to maintain the data transmission between u and v operational, where ξ denotes a random variable of known probability distribution function. Since power consumption in each sensor may vary due to weather or radio interference, sometimes the available power (battery autonomy of sensors) to keep the network operational
M
15
becomes insufficient. This may cause loss of communication among the sensors. To overcome this problem, in this paper we consider probabilistic constraints
ED
related to power consumption in the classic DT problem. The resulting probabilistically constrained dominating tree (PCDT) problem consists in finding a minimum cost distance dominating tree T = (V (T ), E(T )) of G, with set of
PT
20
nodes V (T ) and set of edges E(T ) subject to a probabilistic constraint related to the risk that the total power consumption of each sensor u ∈ V (T ) exceeds a
CE
given limit dm corresponding to its power autonomy. The power consumption of P a sensor u is defined as the sum {u,v}∈E(T ) dξuv of all links in E(T ) connected to u in T .
AC
25
Another classic optimization approach for planning sensor networks consid-
ers the minimum number of connected dominating sensors [11]. In this case, the problem is equivalent to the one in [4, 24] when the distance between sensors equals 1. The novelty here is to provide new ideas for the development of
30
communication protocols for sensor networks by considering both the impact of
ACCEPTED MANUSCRIPT
distances among sensors and the risk that power consumption of sensors exceeds a power threshold available to keep the network operational. In this sense, our work is a first attempt to answer this question by considering jointly or sepa-
35
CR IP T
rately dependence among power consumption of the links of the network in a unified framework.
Both the stochastic mathematical optimization model and the theory to
make it tractable by mixed integer linear programming are contributions to future algorithm developments for practical applications. In this sense, this
step of our research focuses mainly on the theory rather than implementing and
evaluating the overhead of the proposed models in practice. Notice however,
AN US
40
that any dominating tree topology proposed as solution by our model can be used in real implementations as well as in [4, 11, 24]. The reader is referred to [8] for further details on the involved problems in sensor networks. 1.1. An application of a dominating tree sensor network
In Figure 1, we show an example of an undirected connected sensor net-
M
45
work. The links of the graph in Figure 1 represent the possible data exchange
ED
among the 8 sensors {A, B, . . . , H}. Assume that two sensors connected by a link are at a unit distance of each other and that values near each communication link represent, deterministically the power consumption required to keep each link operational. The objective is to determine the minimum distance
PT
50
dominating tree topology for the communication among sensors respecting the power autonomy of each sensor. In fact, if the power autonomy of each sensor is
CE
unlimited, then the optimal dominating tree topology should be the unique link T = ({A, B}) of unit distance, since minimizing the number of connected dominating nodes solves this problem [11]. However, if with each sensor we associate
AC
55
a battery autonomy limit of 5 units to be used with all its communication links, then T = ({A, B}) is no more feasible because the power consumption required by this link is equal to 5.9 units and this value exceeds the available power for sensors A and B. By inspection, the optimal solution for this example is the
60
dominating tree topology T 0 = ({A, C}, {C, D}, {B, D}), of distance equal to 3,
ACCEPTED MANUSCRIPT
with total link power consumption of sensors A, B, C, and D being equal to 1, 1, 5, and 5, respectively. Of course, we assume that the remaining nodes not in T 0 will communicate with their nearest neighbors in T 0 and that the required
E
4
F
1 2
1
H
2 5.9
A C
4
B 1 D
1
G
CR IP T
power to accomplish this task is negligible for the problem.
AN US
Figure 1: Example of an undirected connected sensor network.
Thus, in order to deal with the probabilistic constraints when power con-
65
sumption of the existing links vary according to, e.g., extreme and uncertain weather conditions, we propose a polynomial compact formulation for the DT
M
problem which is based on dominating arborescence of a directed graph and give a proof of its correctness. The compact model allows us to obtain determin70
istic equivalent mixed integer linear programming formulations for the PCDT
ED
problem. Numerical results show that the new approach allows us to solve to optimality random and benchmark instances for the DT and PCDT problems
PT
with up to 120 nodes. To the best of our knowledge, this is the first work reporting proven optimal solutions for challenging benchmark instances of the classic 75
DT problem [4, 26]. Moreover, the probabilistic modelling approach allows us
CE
to handle joint and separate probabilistic constraints in a unified framework.
AC
1.2. Probabilistic constraints Stochastic programming is a powerful technique to deal with the uncertainty
of some input parameters of a mathematical programming problem [23]. These
80
parameters are assumed to behave as random variables which are distributed according to a given probability distribution function. In particular, the probabilistically constrained approach imposes a probability of occurrence for some
ACCEPTED MANUSCRIPT
(or all) of the constraints of a mathematical model. This means that some of the constraints will be satisfied, at least for a given percentage [1, 5, 10, 15, 21, 22]. 85
As far as we know, no attempt has been made to deal with probabilistic con-
CR IP T
straints in the problem of designing sensor networks with a dominating tree topology under uncertain power consumption. Probabilistic constraints can be considered either separately [1, 3, 5, 10, 15, 16, 19, 20, 21, 22] or jointly [5, 6, 7, 17].
In this paper, we consider both joint and separate probabilistic constraints
90
for the PCDT problem. For this purpose, we assume that the row vectors of the
AN US
matrix (dξij ) are joint dependent multivariate normally distributed vectors with
known means and covariance matrices. The dependence of the random vectors is handled by means of copulas [6, 13]. A copula is a distribution function 95
β : [0, 1]|V | → [0, 1] of some |V |-dimensional random vector whose marginals are uniformly distributed on [0, 1]. We apply useful properties of Gumbel-Hougaard copulas to describe the dependence between the rows of the matrix (dξuv ). The
M
latter allows us to handle joint and separate probabilistic constraints in a unified framework. The theory of copula distributions was developed in probability theory and mathematical statistics [13, 18]. Copula distributions allow us to find
ED
100
deterministic convex reformulations for probabilistic constraints while simultaneously handling dependent random variables. In the stochastic programming
PT
field, this has not been well developed yet. In fact, this is the most important reason for using copula distributions so far. Several copula distributions have been proposed such as Clayton, Gumbel and Frank copulas, among others. In
CE
105
this paper, we mainly focus on the use of Gumbel-Hougaard copula as it is the one that leads to more tractable deterministic reformulations as shown in [6],
AC
and also because it allows a unified framework to handle joint and separate probabilistic constraints simultaneously. In Section 4, we give a brief formal
110
description of this copula. In [13], the authors use copulas to come up with convexity results for chance
constrained problems with dependent random right hand side whereas in [6] the authors apply Gumbel-Hougaard copulas to a quadratic optimization prob-
ACCEPTED MANUSCRIPT
lem subject to joint probabilistic constraints. In this paper, we transform the 115
probabilistic constraints into equivalent deterministic second order conic constraints that we further linearize. Consequently, we obtain deterministic equiv-
CR IP T
alent mixed integer linear programming (MILP) formulations for the PCDT problem. Additionally, we investigate the impact of using valid inequalities referred to as generalized sub-tour elimination constraints in all our proposed 120
models [11].
The paper is organized as follows. In Section 2, we introduce both separate and joint probabilistic constraints for the problem. Subsequently, we introduce
AN US
a new compact polynomial formulation for the DT problem and give a proof of its correctness. Then, in Sections 3 and 4 we present deterministic equiva125
lent MILP formulations using separate and joint probabilistic constraints, respectively. Then, in Section 5 we conduct numerical experiments in order to compare all the proposed models for random and benchmark [26] instances for the DT and PCDT problems. Finally, in Section 6 we conclude the paper and
2. Problem formulation
ED
130
M
provide some insights for future research.
In this section, we present a new variant of the classic DT problem while
PT
considering both separate and joint probabilistic constraints that we incorporate in the proposed DT model leading to the PCDT problem formulation. Then, we introduce a compact polynomial formulation for the problem which is based on arborescence of directed graphs and give a proof of its correctness.
CE
135
AC
2.1. Probabilistic dominating tree problem Let G = (V, E) be an undirected sensor network with set of nodes V and
set of connections E. We propose working with an equivalent directed graph H = (V, A), where the set of arcs A is obtained from E by replacing each edge
140
{u, v} ∈ E by the two arcs (u, v) and (v, u) ∈ A. We associate with each arc
(u, v) ∈ A a cost cuv > 0 and an uncertain power consumption dξuv > 0 given
ACCEPTED MANUSCRIPT
by a random variable ξ of known probability distribution function. We assume that the matrices c = (cuv ) and dξ = (dξuv ) are symmetric. Given a directed graph H = (V, A), an arborescence T of H is a rooted directed acyclic subgraph of H in which, given a root vertex r ∈ V , for any other vertex v ∈ V \ {r}, there
CR IP T
145
is exactly one directed path from r to v in T . An arborescence T of a directed
graph H is dominating if every node of H or is in T or is adjacent to a node
in T . The problem is to find a minimum cost dominating arborescence (DA) T of H where the probability that the sum of power consumption of the arcs 150
adjacent to each node u of T exceeds a positive maximum value dm is limited
AN US
above by a certain risk parameter α ∈ [0, 0.5). Notice that the parameter α must
lie in this interval in order to keep the convexity of the deterministic equivalent probabilistic constraints throughout the paper [2, 6, 15].
We represent a DA T of H by a vector x ∈ {0, 1}|A| where xuv = 1 if arc 155
(u, v) belongs to T and xuv = 0, otherwise. We use binary variables yv , for all v ∈ V , to indicate if a node v belongs to T . If it is the case, we set yv = 1;
M
− + and set yv = 0, otherwise. Define N(v) = {u|(u, v) ∈ A}, N(v) = {u|(v, u) ∈ A},
− + N(v) = N(v) ∪N(v) , and N[v] = N(v) ∪{v}. The following probabilistic constraints
ED
model the involved uncertainties X P dξuv xuv ≤ dm , ∀u ∈ V ≥ (1 − α) v∈N +
(1)
160
PT
(u)
where P {·} denotes a probability measure. We assume that the row vectors of
CE
the input matrix d = (duv (ξ)) are dependent multivariate normally distributed vectors with known means and covariance matrices. The probabilistic constraint (1) imposes the condition that each node power constraint must be satisfied at
AC
least for (1 − α)% of the cases where α represents the risk of not satisfying
165
some of these constraints. We remark that the probabilistic constraints (1) are formulated as joint probabilistic constraints under the assumption of row dependence. Next, we further consider the case of separate probabilistic constraints where the row vectors in matrix (dξuv ) are independent multivariate normally
ACCEPTED MANUSCRIPT
distributed vectors. In this case, constraints (1) can be replaced by X P dξuv xuv ≤ dm ≥ (1 − α), ∀u ∈ V v∈N +
(2)
170
CR IP T
(u)
Notice that the statement “∀u ∈ V ” in (1) is inside the probability notation.
This means that the probabilities should be jointly satisfied ∀u ∈ V whereas
in (2), they should be satisfied in a separate way. We handle the dependence between these random vectors by means of copulas [6, 13]. Now, we propose a new compact polynomial formulation for the DT problem that allows to handle efficiently the constraints (1) and (2) as it is based on directed graphs.
AN US
175
2.2. A compact polynomial formulation for the DT problem
Let associate with each node v ∈ V a non-negative continuous variable πv and let fv represent the degree of node v in the undirected graph G. A polynomial compact formulation for the dominating arborescence (DA) problem in the corresponding directed graph H is DTc :
M
180
X
min
{x,y,π}
(3)
X
(4)
ED
X
cuv xuv ,
(u,v)∈A
xuv =
v∈V
(u,v)∈A
yv − 1,
AC
CE
PT
πv − πu + |V |(1 − xuv ) ≥ 1,
∀(u, v) ∈ A,
(5)
πv − πu − |V |(1 − xuv ) ≤ 1, ∀(u, v) ∈ A, X X xvu + xuv ≥ yv , ∀v ∈ V,
(6)
X
+ u∈N(v)
X
u∈N(v)
X
− u∈N(v)
(7)
− u∈N(v)
+ u∈N(v)
xvu +
X
− u∈N(v)
yu ≥ 1,
xuv ≤ fv yv ,
∀v ∈ V,
xuv ≤ 1,
1 ≤ πv ≤ |V | − 2,
∀v ∈ V, ∀v ∈ V,
x ∈ {0, 1}|A| , y ∈ {0, 1}|V | .
∀v ∈ V,
(8) (9) (10) (11) (12)
ACCEPTED MANUSCRIPT
where cuv represents the cost of arc (u, v) ∈ A;
CR IP T
xuv is a binary variable denoting if arc (u, v) belongs (xuv = 1) or not (xuv = 0) to the DA; 185
πv is a continuous variable corresponding to the position in the sequence of arcs starting from the vertex of the DA that has no incoming arc to v; and
yv is a binary variable indicating if vertex v belongs (yv = 1) or not (yv = 0)
190
AN US
to the DA;
The total cost of a feasible DA T is minimized in the objective function (3). P Constraint (4) establishes that T has v∈V yv − 1 arcs. Constraints (5)-(6) ensure that if an arc (u, v) belongs to the solution, then πv − πu = 1; otherwise
both constraints become redundant. These constraints together with (10) avoid cycles in any feasible solution. The idea behind these constraints is to explore
195
M
simultaneously the space of primal and dual variables of the problem structure, based on a distance property. Indeed, every two end nodes of any arc in the
ED
solution are at a unit distance from each other. In fact, for this problem, we can show that (6) are redundant in the presence of (10). Concerning constraints (7), they establish that if node v is in the solution,
200
PT
then at least one arc is adjacent to v. Constraints (8) are valid inequalities for DTc , indicating that if node v is not in the solution, then no arc is adjacent to
CE
v; otherwise, at most fv arcs can be adjacent to this node. Constraints (9) are domination constraints and ensure that any node v in V has a neighbor in T . Constraints (10) limit the number of arcs arriving at a node v ∈ V to at
AC
most one. Constraints (11) bound the possible values for the potential variables
205
π. This is due to the fact that any directed path in the solution arborescence spans at most |V | − 2 nodes. In fact, since DTc is a minimization problem and
the π variables have one as lower bound, the upper bound |V | − 2 on these
variables is redundant. Finally, (12) are the domain of the decision variables.
ACCEPTED MANUSCRIPT
Theorem 1. An optimal solution of DTc referred to H gives a dominating tree 210
of G of minimum cost. Proof. Initially, we show that any feasible solution (¯ x, y¯, π ¯ ) of DTc is an
CR IP T
acyclic connected arborescence of H. Suppose that x ¯ induces a cycle of arcs
(v1 , v2 ), (v2 , v3 ), · · · , (vp , v1 ), with vi ∈ V for i = 1, . . . , p, for some p ≤ |V | − 3.
By constraints (5) and (6), πv2 = πv1 + 1, πv3 = πv2 + 1, . . ., πvp = πvp−1 + 1 and πv1 = πvp + 1. This implies that πv1 < πv2 < · · · < πvp < πv1 cannot hold.
Thus, the directed graph induced by x ¯ is acyclic. Now suppose that x ¯ is acyclic
and w.l.o.g. induces two connected arborescences whose sets of nodes induced
AN US
by the non null components of y¯ are denoted by D1 ⊂ V and D2 ⊂ V , with
D1 ∪ D2 being a non-connected dominating set of G. By constraints (7), both extremities u and v of an arc (u, v) such that x ¯uv = 1 belong to the solution, i.e., y¯u = y¯v = 1. Note that arcs corresponding to non null components of x ¯ have both extremities or in D1 , denoted by A(D1 ), or in D2 , denoted by A(D2 ).
X
x ¯uv = |D1 | − 1,
X
x ¯uv = |D2 | − 1,
ED
connected, we have that
M
As the directed graphs induced by D1 and D2 are assumed to be individually
(u,v)∈A(D1 )
PT
and
CE
Moreover,
X
(u,v)∈A(D2 )
x ¯uv =
(u,v)∈A
X
(u,v)∈A(D1 )
and
AC
X
v∈D1
Consequently,
(u,v)∈A
x ¯uv =
X
v∈D1
y¯v − 1 +
x ¯uv ,
y¯v = |D1 |, y¯v = |D2 |.
X
y¯v − 1 =
v∈D2
X
(u,v)∈A(D2 )
X
v∈D2
X
x ¯uv +
X
v∈D1 ∪D2
y¯v − 2 =
X
v∈V
y¯v − 2,
ACCEPTED MANUSCRIPT
which violates constraint (4). Therefore, any optimal solution of DTc must be a connected arborescence of G. Finally, by dropping the orientation of its arcs we obtain a dominating tree of minimum cost of G.
215
CR IP T
In Figure 2 we show some ingredients involved in the proof of Theorem 1.
Initially, we determine the corresponding directed graph H from the undirected graph in Figure 1. A bi-directed arrow represents two arcs that replace the undirected link between its corresponding sensors.
For the sake of comparison with PCDT, the power required by each arc is
220
AN US
equal to that of the original link as depicted in the directed graph H. Assume
that the power available at each sensor is of 5 units. The cost between two sensors is equal to 1. We show five situations avoided or allowed by this theorem to show how an optimal DA leads to an optimal DT while guaranteeing acyclic and connected solutions.
If we were interested in monitoring the sensor network area by using communication protocols as suggested in [8], where the main purpose is to determine
M
225
the minimum connected dominating set of sensors [11], the solution should be the arborescence {(B, A)} in (i), of unitary cost distance, constituted of only one
ED
arc and two dominating sensors B and A. The remaining sensors are dominated and communicate with the dominating ones by using the dotted links representing their nearest neighbor (sensor) in this arborescence. However, if we intend
PT
230
to consider the interference of natural events of the sensor network on the power consumption between sensors B and A (because the required power to keep this
CE
communication link exceeds the available power at these sensors), this solution is no longer feasible. On the other hand, in (ii) when we relax the connectivity requirement, we obtain two arborescences {(A, C)} and {(D, B)} that together
AC
235
constitute a (non-connected) dominating set presenting minimum cost distance equal to two. Their arcs power consumption respect the battery autonomy of the involved sensors, but the number of arcs is two while the number of dominating sensors is four, thus violating constraint (4). Unfortunately, as we can
240
see in this situation, the use of cut inequalities [11] is not helpful to force ob-
ACCEPTED MANUSCRIPT
taining a connected solution because the minimum number of non-dominating sensors needed to connect these two sets is zero. Recall that the term “cut inequality” [11] is a valid inequality defined on
245
CR IP T
the decision variables y to cut off an integer solution associated with a proper subset S = {i|yi = 1} of vertices inducing a disconnected graph. If a subset of
vertices S is not connected, then at least one vertex in S¯ = V \ S must belong to P the solution and the corresponding cut inequality is i∈S¯ yi ≥ 1. These valid inequalities are usually used within exact algorithms [11].
Subsequently, notice that individually the sets {A, C} and {D, B} are not dominating. Although, their corresponding arborescences can be connected by
AN US
250
adding, for instance, the arc (C, D). We cannot impose that at least an arc must be connected to the set {A, C} or to {D, B} because their complement with respect to the remaining sensors are dominating sets [11]. When a set S and ¯ are both non-dominating, at least an arc must connect S its complement, say S, 255
to S¯ [11]. In (iii), we have a fictitious cycle induced by the arcs (C, A), (A, B),
M
(C, D) and (D, B). We show why it cannot belong to a feasible DA. Initially, if we assume that πC = π ¯ , constraints (5) and (6) impose that πA = πD = π ¯+1
ED
and πB = π ¯ + 2. Thus, these constraints are not sufficient to avoid cycles of this type in the absence of constraints (10) that are required for the model to 260
be correct. The optimal DA for this example is {(C, A), (C, D), (D, B)}, rooted
PT
at sensor C as depicted in (iv), and has cost distance equal to three. The sum of the arcs power consumption at sensors A, B, C, and D do not exceed the
CE
available power at each sensor. The values for the π variables associated with this solution can be πC = 1, πA = πD = 2, πB = 3. The remaining values of
265
the π variables are irrelevant for this solution. By dropping the orientation of
AC
the arcs in this DA, we obtain the DT in (v) that is the optimal solution for the undirected graph in Figure 1. The fact that DTc is represented as an arborescence of a directed graph
allows us to handle efficiently the probabilistic constraints (1) and (2).
270
Hereafter, we refer to the PCDT problem as the problem obtained by the constraints defining the DTc model with the additional probabilistic constraints
ACCEPTED MANUSCRIPT
Digraph H
4 E
F 5.9
A
B
2
H
G C
E
F
A
B
1 C
F
5.9
A
B 1
H
D
G
D
4 E
F 5.9
A
C
G
H
D
1
C
4
E
E
F
A 1
1 C
D
4
F
A
G
H
B
1
1
C
M
H
B
G
D
(iii)
(ii)
(i)
B
1
AN US
E
1
1
H
1
CR IP T
2
1
4
G
D
(v)
(iv)
Figure 2: Example of how Theorem 1 works to determine a dominating tree from an optimal
ED
dominating arborescence while avoiding infeasible solutions.
PT
(2) or (1), for the separate or joint probabilistically constrained approaches, respectively. We further consider the resulting MILP model obtained from DTc
AC
CE
while replacing constraints (8) with the following set of constraints
275
X
xuv
≤
yu ,
∀(u, v) ∈ A,
(13)
xuv
≤
(14)
xuv
≤
yv , ∀(u, v) ∈ A, X yj − 1, ∀i ∈ V ,
(u,v)∈A(N[i] )
(15)
j∈N[i]
and constraints (10) with the following ones X
− u∈N(v)
xuv ≤ yv ,
∀v ∈ V.
(16)
ACCEPTED MANUSCRIPT
We denote this resulting MILP model by DTcV i . In particular, constraints (15) are referred to as a strengthened version of the generalized sub-tour elimination constraints [11]. Finally, we denote hereafter the corresponding linear
280
CR IP T
programming (LP) relaxations of DTc and DTcV i by LPc and LPcV i , respectively. 3. DT problem with separate probabilistic constraints
In this section, we propose a deterministic equivalent MILP formulation for
the PCDT problem using separate probabilistic constraints. For this purpose,
we transform the probabilistic constraints (2) into deterministic equivalent sec-
285
AN US
ond order conic constraints [15]. Subsequently, we linearize the conic constraints to obtain an MILP formulation for the PCDT problem.
In order to obtain a deterministic MILP formulation, we assume that the row vectors in matrix (dξuv ) in (2) are independent multivariate random variables normally distributed with known means d¯uv , ∀(u, v) ∈ A. Also, let (σ u ) = 290
M
u (σiv ), ∀ (i, v) ∈ A, u ∈ V , be the root square matrix of the corresponding covari-
ance matrices of each row vector (dξu• ), ∀ u ∈ V , that we assume to be positive
semidefinite. In this case, notice that the constraints (2) can be equivalently
ED
written as [15, 22]
PT
v u u uX X d¯uv xuv + F −1 (1 − α)u t
+ v∈N(u)
i∈V
X
+ v∈N(u) \{i}
2
ux ≤d , σiv uv m
∀ u ∈ V, (17)
CE
where F −1 (·) denotes the inverse of the standard normal cumulative distribution function. The condition for the parameter α ∈ [0, 0.5) is mandatory for the constraints (17), otherwise these constraints would not be convex. Proposition 1. A deterministic equivalent MILP formulation for the PCDT
AC
295
problem with separate probabilistic constraints can be written as M IP1
min
{x,y,π,ϕ}
X
cuv xuv ,
(u,v)∈A
subject to (4) − (12) and
(18)
ACCEPTED MANUSCRIPT
(F −1 (1 − α))2 d2m +
X
i∈V
X
u 2 (σiv ) xuv +
+ v∈N(u) \{i}
(−2dm d¯uv + d¯2uv )xuv +
+ v∈N(u)
X
+ + v∈N(u) k∈N(u) \{i,v}
X
+ v∈N(u)
d¯uv d¯uk ϕkuv ,
(19)
∀u ∈ V,
xuv + xuk − 1 ≤ ϕkuv ≤ xuv
(20)
∀u, v, k ∈ V, (u 6= v 6= k),
∀u, v, k ∈ V, (u 6= v 6= k),
Proof. We write constraint (17) equivalently as [15, 20] v 2 u u uX X X ux ≤d − F −1 (1 − α)u σiv d¯uv xuv , uv m t + v∈N(u) \{i}
+ v∈N(u)
M
Notice that in order for (24) to make sense, we must have
(23)
∀u ∈ V. (24) P
+ v∈N(u)
d¯uv xuv ≤
dm , ∀u ∈ V . Thus, taking the square of both sides in (24) we obtain 2 2 X X X u (F −1 (1 − α))2 σiv xuv ≤ dm − d¯uv xuv ,
ED
300
(21)
(22)
AN US
ϕ ∈ {0, 1}|V ×V ×V | .
i∈V
u u k σiv σik ϕuv ≤
+ + v∈N(u) k∈N(u) \{v}
∀u ∈ V, X d¯uv xuv ≤ dm , ϕkuv ≤ xuk
X
X
CR IP T
X
i∈V
+ v∈N(u) \{i}
+ v∈N(u)
CE
PT
that after solving each quadratic term, is equivalent to X X X X u 2 (F −1 (1 − α))2 (σiv ) xuv + d2m +
X
AC
+ v∈N(u)
i∈V
+ v∈N(u) \{i}
(−2dm d¯uv + d¯2uv )xuv +
X
+ + v∈N(u) k∈N(u) \{i,v}
X
∀u ∈ V,
u u σiv σik xuv xuk ≤
d¯uv d¯uk xuv xuk ,
+ + v∈N(u) k∈N(u) \{v}
∀u ∈ V.
(25)
Finally, if we define ϕkuv = xuv xuk , ∀u, v, k ∈ V , with u 6= v 6= k, and use
standard methods [9, 12] to linearize ϕkuv as shown in (21) and (22), then (25) can be written as in (19). Finally, constraints (4)-(12) ensure that feasible
305
solutions are dominating arborescences, thus concluding the proof.
ACCEPTED MANUSCRIPT
We further consider the resulting MILP model obtained from M IP1 while replacing constraints (8) with the set of constraints (13)-(15) and constraints (10) with the constraints (16). Hereafter, we denote this MILP model by M IP1V i .
310
CR IP T
Finally, we denote the corresponding LP relaxations of M IP1 and M IP1V i by LP1 and LP1V i , respectively.
4. DT problem with joint probabilistic constraints
In this section, we obtain a deterministic equivalent MILP formulation for the PCDT problem using joint probabilistic constraints. We transform the prob-
315
AN US
abilistic constraints into deterministic equivalent second order conic constraints and linearize them in order to obtain an MILP formulation. In particular, our proposed model is based on the family of Gumbel-Hougaard copula. For a
deeper comprehension related to the theory of copulas we refer the reader to [25] and to the book [18]. The existence and uniqueness of a copula, for any
320
M
distribution function and all its marginals are guaranteed by Sklar’s Theorem [18]. In this paper, we consider the Gumbel-Hougaard copula which can be 1/θ |V | X θ (− ln tk ) βθ (t) = exp − k=1
ED
written as
(26)
β(t) =
|V | Y
tk
k=1
CE
PT
for given θ ≥ 1. In particular when θ = 1, this copula takes the form
and represents the joint distribution of independent random variables. Next,
AC
consider the following transformations P
+ v∈N(u)
ωu (x) = v u u uP t i∈V
(dξuv − d¯uv )xuv P
+ v∈N(u) \{i}
2
ux σiv uv
ACCEPTED MANUSCRIPT
dm − gu (x) = v u u uP t i∈V
where
P
i∈V
P
+ v∈N(u) \{i}
P
d¯uv xuv
+ v∈N(u)
P
+ v∈N(u) \{i}
2
2
ux σiv uv
CR IP T
and
u σiv xuv is strictly positive. The random variable ωu (x)
has one-dimensional standard normal distribution which is independent of the
variable x = (xuv ). Note that the probabilistic constraint (1) can be equivalently written as
AN US
325
P {ωu (x) ≤ gu (x), ∀u ∈ V } ≥ (1 − α)
(27)
Lemma 1. [6] If the random vector (ω1 (x), . . . , ω|V | (x))T , where ωu (x) has one-dimensional standard normal distribution, has a joint distribution driven
330
M
by the Gumbel-Hougaard copula βθ with some θ ≥ 1, then the constraint (1) is equivalent to the set of constraints
+ v∈N(u)
ED
v u u u X X 1/θ u d¯uv xuv + F −1 (1 − α)zu t
PT
∀u ∈ V X zu = 1, u∈V
zu ≥ 0,
∀u ∈ V
i∈V
X
+ v∈N(u) \{i}
2
ux ≤d , σiv uv m
(28) (29)
CE
where F −1 (·) is the inverse of the standard normal cumulative distribution function.
AC
Notice that when the term
P
+ v∈N(u)
xuv equals zero for a particular u ∈ V ,
constraints (28) are trivially satisfied. Moreover, the value of zu = 0 in this case
335
since the root square term in (28) equals zero as well. The following convexity lemma is also required Lemma 2. [6] If (1 − α) ≥ [0, 1].
1 2
1/θ and θ ≥ 1, then F −1 (1 − α)zu is convex on
ACCEPTED MANUSCRIPT
340
Hence, in this paper we assume that 12 ≤ (1 − α) < 1. Notice that when 1/θ θ → ∞, the function F −1 (1 − α)zu → F −1 (1 − α) which is equivalent to a separate probabilistic constraint [15, 22]. This shows that the Gumbel-
CR IP T
Hougaard copula allows us to handle joint and separate probabilistic constraints in a unified framework. Also notice that for values of θ ∈ [1, ∞) each zu , ∀u ∈ V
in constraint (29) will always be positive, otherwise the constraint (28) would 345
never be satisfied since F −1 (1) → ∞, unless the root square term equals zero.
In this case, we assume that the row vectors in matrix (dξuv ) in (1) are de-
pendent multivariate random variables normally distributed with known means
AN US
u d¯uv , ∀(u, v) ∈ A. Also, let (σ u ) = (σiv ), ∀(i, v) ∈ A, u ∈ V be defined as in the
previous section. Thus, a deterministic equivalent MILP model can be obtained 350
by means of the following proposition.
Proposition 2. Assume that there exists a sufficiently large number of line
M
segments l = {1, . . . , L} for fixed α ∈ [0, 0.5) and θ ≥ 1 in order to approximate 2 1/θ from above the univariate convex function F −1 (1 − α)zu on [0, 1], ∀u ∈
V , then the PCDT problem obtained by replacing the probabilistic constraint (1) with the deterministic joint probabilistic constraints (28)-(29) can be equivalently written as
min
{x,y,π,λ,µ,ϕ,φ,z}
CE AC
X
cuv xuv ,
(30)
(u,v)∈A
subject to (4) − (12), (20) − (23), (29) and X X X X u 2 (σiv ) λuv +
PT
M IP2
ED
355
i∈V
d2m +
+ v∈N(u) \{i}
X
+ + v∈N(u) k∈N(u) \{i,v}
(−2dm d¯uv + d¯2uv )xuv +
+ v∈N(u)
X
+ v∈N(u)
u u k σiv σik µuv ≤
X
d¯uv d¯uk ϕkuv ,
+ k∈N(u) \{v}
∀u ∈ V,
(31)
φu ≥ al + bl zu ,
∀u ∈ V, l = {1, . . . , L},
Mxuv + φu − M ≤ λuv ≤ Mxuv , λuv ≤ φu , λuv ≥ 0
∀u, v ∈ V, (u 6= v), ∀u, v ∈ V, (u 6= v),
∀u, v ∈ V, (u 6= v),
(32) (33) (34) (35)
ACCEPTED MANUSCRIPT
Mϕkuv + φu − M ≤ µkuv ≤ Mϕkuv , ∀u, v, k ∈ V, (u 6= v; u, v 6= k),
(36)
µkuv ≤ φu ,
(37)
∀u, v, k ∈ V, (u 6= v; u, v 6= k).
(38)
CR IP T
µkuv ≥ 0,
∀u, v, k ∈ V, (u 6= v; u, v 6= k),
Proof. From Lemma 1 and considering (20), we can write (28) as
i∈V
+ v∈N(u) \{i}
or equivalently
1/θ
(F −1 ((1 − α)zu ))2 d2m +
X
X
i∈V
X
(−2dm d¯uv + d¯2uv )xuv +
X
+ v∈N(u)
X
X
X
d¯uv xuv ,
+ + v∈N(u) k∈N(u) \{i,v}
X
2
u u σiv σik xuv xuk ≤
d¯uv d¯uk xuv xuk ,
+ + v∈N(u) k∈N(u) \{v}
(39)
As before, we define ϕkuv = xuv xuk , ∀u, v, k ∈ V , with u 6= v 6= k, and linearize 1/θ
these variables according to (21) and (22). Now, let φu ≥ (F −1 ((1 − α)zu ))2 =
ED
360
al + bl zu , ∀u ∈ V, l = {1, . . . , L} be a tight linear approximation of (F −1 ((1 − 1/θ
α)zu ))2 from above for L sufficiently large where each (al , bl ) denotes the coeffi-
PT
cients of each linear segment uniformly selected on the interval (0,1]. Notice that this interval is implied by constraints (29). Also, define nonnegative variables λuv = φu xuv , ∀ u, v ∈ V, (u 6= v) and µkuv = φu ϕkuv , ∀ u, v, k ∈ V, (u 6= v 6= k).
CE
365
The λ and µ variables can be linearized according to (33)-(35) and to (36)-(38), respectively, with M being a bigM positive constant. Then, after developing 1/θ
AC
(39) and approximating (F −1 ((1−α)zu ))2 by φu , ∀u ∈ V ; and replacing φu xuv
by λuv , ∀ u, v ∈ V, (u 6= v); and φu ϕkuv by µkuv , ∀ u, v, k ∈ V, (u 6= v 6= k), we
370
∀u ∈ V,
M
∀u ∈ V.
u σiv xuv ≤ dm −
u 2 (σiv ) xuv +
+ v∈N(u) \{i}
+ v∈N(u)
2
X
AN US
X 1/θ (F −1 ((1 − α)zu ))2
obtain the resulting constraints (31) under the condition that (32) and (29) must be verified in order to guaranty that the φu variables represent a very close ap1/θ
proximation of (F −1 ((1 − α)zu ))2 , ∀u ∈ V . Finally, the remaining constraints
ACCEPTED MANUSCRIPT
of this model ensure that feasible solutions are dominating arborescences, thus concluding the proof.
CR IP T
As a direct consequence of Proposition 2, we have the following corollary.
375
Corollary 1. The optimal objective function value of M IP2 is an upper bound for the optimal objective function value of M IP1 .
Proof. The proof is immediate due to the convexity of the function (F −1 ((1 − 1/θ
α)zu ))2 for values of zu ∈ [0, 1], α ∈ [0, 0.5) and θ ≥ 1 since (F −1 ((1 − 1/θ
α)zu ))2 ≥ (F −1 (1 − α))2 . This clearly implies that the feasible set of M IP2 is
AN US
380
contained in the feasible set of M IP1 .
Similarly as for M IP1 , we further consider the resulting MILP model obtained from M IP2 by replacing the constraints (8) with the set of constraints (13)-(15) and the constraints (10) with the constraints (16). Finally, we denote this new 385
model by M IP2V i and the corresponding LP relaxations of M IP2 and M IP2V i
M
by LP2 and LP2V i , respectively.
ED
5. Numerical results
In this section, we present numerical results for all the proposed models:
390
PT
DTc , DTcV i , LPc , LPcV i , M IP1 , M IP1V i , LP1 , LP1V i , M IP2 , M IP2V i , LP2 and LP2V i . We implement a Matlab program using CPLEX 12.6 for solving the MILP and LP models [14]. The numerical experiments have been carried
CE
out on an Intel(R) 64 bits core (TM) with 3.4 GH and 8 GB of RAM under Windows 7. We consider only connected disk graphs G = (V, E) where each
AC
disk represents the transmission range of a node v ∈ V . All instances considered
395
in this work are generated according to the procedure reported in [4, 26] which 2 is as follows: the cost of each edge {u, v} ∈ E is determined as cuv = wuv where
wuv represents the Euclidean distance between nodes u and v. An edge {u, v} between nodes u and v of V belongs to E if cuv is smaller than or equal to the considered node transmission range. In each instance, we distribute randomly
ACCEPTED MANUSCRIPT
400
|V | nodes in an area of 500 × 500 m2 and consider node transmission ranges equal to 100 m, 150 m and 350 m. The same transmission range applies for all nodes of a given instance. We further consider 18 benchmark DT instances
CR IP T
from [4, 26] with |V | = 50 and |V | = 100 with distinct transmission ranges of 100 m, 125 m and 150 m.
Concerning the uncertain parameters, each entry in matrix d¯ = (d¯uv ),
405
∀{u, v} ∈ E, is uniformly and randomly generated from the interval (0, 20].
u Similarly, we generate auxiliary matrices σau of entries (σiv )a , ∀u ∈ V , whose
values of their elements are uniformly and randomly drawn from the interval
410
AN US
[0, 1]. Each matrix σau , ∀u ∈ V , is multiplied by its transpose, thus originating
a new symmetric positive semidefinite matrix σ ˆ u := (σau ) × (σau )T . Finally, we
determine each matrix σ u from the one σ ˆ u by setting σ u := (ˆ σ u )1/2 , ∀u ∈ V
where (ˆ σ u )1/2 represents the square root of σ ˆ u . The values of α and dm are set to α = 0.1 and dm =
max d¯uv × γ, respectively, where γ ∈ [0.5, 3.2] is
{u,v}∈E
arbitrarily chosen. For the models M IP2 , M IP2V i , LP2 and LP2V i , we set θ = 1 1/θ
in the function (F −1 ((1 − α)zu ))2 . Notice that θ = 1 is the smallest value for
M
415
this parameter. It allows to obtain more conservative optimal solutions in the 1/θ
ED
corresponding models. When θ → ∞, then F −1 ((1 − α)zu ) → F −1 (1 − α) which corresponds to the case of separate probabilistic constraints that we use in our models M IP1 , M IP1V i , LP1 , and LP1V i . Considering values of θ in [1, ∞), we can obtain optimal solutions for both joint to separate probabilistic
PT
420
constraints in a unified framework. With respect to the line segments used in 1/θ
CE
the linear approximation of (F −1 ((1 − α)zu ))2 , we generate values of al and bl , ∀l ∈ {1, . . . , L}, by considering the interval (10−5 , 0.99995), for each zu ∈ V ,
with equally spaced subintervals of length 10−2 .
AC
425
In Figure 3, we plot the function (F −1 ((1 − α)z
1/θ
))2 for different values of
θ ∈ {1, 2, 10, 200}. From these curves, we confirm that this function is convex
and that for large values of θ, the curves become constant. In our numerical experiments, we observed that the length 10−2 is enough for these models to guarantee that they can find the optimal solution of each instance. Also, as the
ACCEPTED MANUSCRIPT
16 θ =1 θ =2 θ = 10 θ = 200
12
8 6 4 2 0
0
0.1
0.2
0.3
0.4
0.5
Values of z
0.8
0.9
1
Curves for the normal inverse cumulative distribution function for θ ∈
{1, 2, 10, 200}.
430
0.7
AN US
Figure 3:
0.6
CR IP T
10
1
1 22 2 1/θ F −1 (1 − α)z
14
range variation (see the vertical axis) of the curves in Figure 3 are not larger than 15 for values of z ∈ (0, 1), we set M = 20 for the joint models M IP2 , M IP2V i , LP2 and LP2V i .
M
We adopt the following organization for the tables presented hereafter. Each table reports results obtained with two models specified accordingly in the set of columns referred to M odelI and M odelII . In the last three rows of each table we
ED
435
present minimum, maximum and average (when all column data are available and, if not, we indicate as ‘na’) values for each column. Column (#) gives the
PT
instance identifier. The number of nodes and the density of edges in percentage (considering the maximum possible number of edges) of each instance are given by columns (|V |) and (Den %), respectively. Columns (Opt), (B&Bn),
CE
440
(T ime(s)), (LR) and (T imeLR (s)) of each ‘M odel∗ ’ give, respectively, the optimal solution value or the best incumbent solution value obtained by CPLEX
AC
in a time limit of 10 hours, the number of CPLEX branch and bound nodes, the CPU time in seconds to obtain the optimal or incumbent solution, the value
445
of the linear relaxation of the corresponding MILP formulation and the CPU time in seconds to obtain the linear relaxed solution. Finally, integrality gaps (in percentage) are reported in columns (GapI %) and (GapII %) referred to
ACCEPTED MANUSCRIPT
#
|V |
Den %
1
10
2
M odelII = DTcV i
M odelI = DTc
Integrality gaps
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
80.0
5954.64
0
0.18
994.06
0.20
5954.64
0
0.17
5954.64
0.16
83.31
0
12
80.3
2178.31
0
0.27
2178.31
0.25
2178.31
0
0.28
2178.31
0.25
0
0
3
14
81.3
11547.03
34
0.20
1056.38
0.16
11547.03
0
0.20
11547.03
0.19
90.85
0
4
15
58.1
29432.11
4
0.26
11550.21
0.17
29432.11
0
0.18
29432.11
0.19
60.76
0
5
16
84.2
1361.40
0
0.19
1361.40
0.19
1361.40
6
18
64.1
1665.67
0
0.20
838.55
0.19
1665.67
7
20
71.6
2041.29
0
0.19
1088.98
0.18
2041.29
8
22
87.9
62.34
0
0.20
62.34
0.17
62.34
9
24
65.9
18423.85
1194
0.53
3537.69
0.20
18423.85
10
25
76.3
2650.94
0
0.20
123.33
0.17
2650.94
11
30
76.8
3423.67
0
0.25
1124.68
0.18
3423.67
12
35
76.0
1755.49
0
0.23
916.62
0.21
1755.49
13
40
70.5
105.81
0
0.23
105.81
0.19
105.81
14
45
62.0
1292.91
0
0.29
443.64
0.21
1292.91
15
50
73.1
899.27
0
0.27
481.87
0.21
899.27
0
16
60
80.3
307.19
0
0.36
64.70
0.24
307.19
17
70
72.8
1473.30
508
1.61
231.18
0.26
18
80
78.4
482.77
3
0.66
156.02
19
90
76.3
5.63
0
0.59
5.63
20
100
69.1
2159.86
32990
14.95
21
110
71.5
192.58
0
1.03
22
120
79.0
330.56
0
1.20
Min.
58.10
5.63
0
0.18
Max.
87.90
29432.11
32990
Ave.
74.34
3988.48
1578.77
CR IP T
Opt
0.17
1361.40
0.20
0
0
0
0.19
1665.67
0.19
49.66
0
0
0.18
2041.29
0.17
46.65
0
0
0.20
62.34
0.19
0
0
0
0.20
13968.51
0.17
80.80
24.18
0
0.21
2650.94
0.20
95.35
0
0
0.25
3423.67
0.21
67.15
0
0
0.26
1755.49
0.22
47.79
0
0
0.31
105.81
0.23
0
0
0
0.30
1292.91
0.26
65.69
0
0.38
899.27
0.28
46.42
0
0
0.59
307.19
0.45
78.94
0
1473.30
0
0.81
876.69
0.48
84.31
40.49
0.34
482.77
0
1.10
482.77
0.64
67.68
0
0.42
5.63
0
1.45
5.63
0.79
0
0
49.50
0.40
2159.86
0
1.61
2159.86
0.90
97.71
0
5.50
0.54
192.58
0
2.10
192.58
1.21
97.14
0
27.90
0.65
330.56
0
3.12
330.56
1.71
91.55
0
5.50
0.16
5.63
0
0.17
5.63
0.16
0
0
14.95
11550.21
0.65
29432.11
0
3.12
29432.11
1.71
97.71
40.49
1.09
1200.19
0.26
3988.48
0
0.64
3758.84
0.42
56.89
2.93
M
AN US
0
ED
Table 1: Numerical results for models DTc and DTcV i using a transmission range of 350 m.
M odelI and M odelII , respectively. These gaps, for each model, are obtained according to the expression 100 Opt−LR Opt . In Tables 1-3, 4-6 and 7-9, we present numerical results for DTc , LPc , DTcV i , LPcV i , M IP1 , LP1 , M IP1V i , LP1V i , and
PT
450
M IP2 , LP2 , M IP2V i , LP2V i , respectively, for graph instances with transmission
CE
ranges of 350 m, 150 m and 100 m, respectively. The set of instances is the same in all these tables. This allows to compare the objective function values for all the proposed models with and without probabilistic constraints. In particular, it allows to compare all the models using separate and joint probabilistic con-
AC 455
straints.
We observe in Tables 1-3 that models DTc , DTcV i , M IP1 , M IP1V i
were able to solve to optimality instances with up to 120 nodes, while with model M IP2V i we were able to solve to optimality only instances with up to 60 nodes. This shows that solving the problem with joint probabilistic constraints
ACCEPTED MANUSCRIPT
M odelII = M IP1V i
M odelI = M IP1
Integrality gaps
|V |
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
1
10
11190.53
0
0.21
994.06
0.24
11190.53
0
0.27
5954.64
0.19
91.12
46.79
2
12
16008.16
0
0.36
2178.31
0.30
16008.16
0
0.30
2178.31
0.28
86.39
86.39
3
14
68487.36
0
0.25
1056.38
0.23
68487.36
0
0.28
11547.03
0.23
98.46
83.14
4
15
111723.24
0
0.37
11550.21
0.24
111723.24
0
0.36
29432.11
0.24
89.66
73.66
5
16
5850.25
0
0.33
1361.40
0.27
5850.25
0
6
18
29967.21
0
0.51
838.55
0.31
29967.21
0
7
20
101438.94
0
0.52
1088.98
0.37
101438.94
0
8
22
4061.08
0
0.73
62.34
0.47
4061.08
0
9
24
91649.72
37
1.28
3537.69
0.53
91649.72
0
10
25
5097.69
0
0.96
123.33
0.62
5097.69
0
11
30
81299.10
0
1.36
1124.68
1.09
81299.10
0
12
35
41773.18
0
2.28
916.62
1.59
41773.18
0
13
40
113464.73
0
3.44
124.56
2.38
113464.73
0
14
45
75417.12
0
6.63
443.64
3.17
75417.12
0
15
50
42249.56
0
7.33
481.87
16
60
47108.72
0
8.89
64.70
17
70
91211.93
0
13.39
231.18
18
80
57358.00
0
23.12
156.02
19
90
59042.83
0
36.55
5.63
20
100
37954.74
79749
135.12
49.50
21
110
55372.76
65
75.59
5.50
22
120
48777.20
20
113.13
27.90
Min.
4061.08
0
0.21
5.50
Max.
113464.73
79749
135.12
11550.21
Ave.
54386.54
3630.50
19.65
1361.40
0.27
76.73
76.73
0.42
1665.67
0.31
97.20
94.44
0.51
2041.29
0.38
98.93
97.99
0.78
62.34
0.47
98.46
98.46
1.22
14350.58
0.56
96.14
84.34
1.00
2650.94
0.65
97.58
48.00
1.62
3423.67
1.04
98.62
95.79
2.43
1760.55
1.65
97.81
95.79
3.62
730.10
2.40
99.89
99.36
6.50
1420.45
3.27
99.41
98.12
AN US
0.34
5.01
42249.56
0
8.13
899.27
5.01
98.86
97.87
9.98
47108.72
0
9.59
307.19
10.01
99.86
99.35
15.55
91211.93
0
14.07
1177.32
16.34
99.75
98.71
26.52
57358.00
0
24.28
531.59
26.52
99.73
99.07
40.54
59042.83
0
37.01
5.63
40.37
99.99
99.99
53.81
37954.74
0
85.95
2159.86
53.16
99.87
94.31
77.67
55372.76
0
78.90
192.58
77.15
99.99
99.65
121.83
48777.20
0
120.40
330.56
122.79
99.94
99.32
0.23
4061.08
0
0.27
5.63
0.19
76.73
46.79
121.83
113464.73
0
120.40
29432.11
122.79
99.99
99.99
16.48
54386.54
0
18.09
3826.50
16.51
96.56
89.42
M 1201.04
CR IP T
#
Table 2: Numerical results for models M IP1 and M IP1V i using a transmission range of 350
460
ED
m.
is significantly harder than solving it with separate probabilistic constraints.
PT
Note that, in general, optimal solution values for these instances follow the relations Opt(DTc ) = Opt(DTcV i ) ≤ Opt(M IP1 ) = Opt(M IP1V i ) ≤ Opt(M IP2 ) =
CE
Opt(M IP2V i ). Moreover, we observe that optimal solution values of M IP1 and
M IP2 are significantly larger than those of the classic problem DTc . In fact, in
465
these tables, optimal solution values of M IP1 and M IP2 are the same with the
AC
exception of the instance 4 which is larger for M IP2 . We also see, in terms of CPU times and number of CPLEX branch and bound nodes, that models DTcV i ,
M IP1V i and M IP2V i are more effective, presenting significantly smaller values than those obtained with models DTc , M IP1 and M IP2 , respectively. This can
470
be explained by the fact that, for the models DTcV i , M IP1V i and M IP2V i , we
ACCEPTED MANUSCRIPT
M odelII = M IP2V i
M odelI = M IP2
Integrality gaps
|V |
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
1
10
11190.53
0
0.33
994.06
0.21
11190.53
0
0.37
5954.64
0.25
91.12
46.79
2
12
16008.16
0
0.62
2178.31
0.27
16008.16
0
0.50
2178.31
0.28
86.39
86.39
3
14
68487.36
0
0.64
1056.38
0.28
68487.36
0
0.47
11547.03
0.27
98.46
83.14
4
15
165378.94
0
1.75
11550.21
0.30
165378.94
0
1.30
29432.11
0.29
93.02
82.20
5
16
5850.25
0
1.31
1361.40
0.34
5850.25
0
6
18
29967.21
0
2.28
838.55
0.41
29967.21
0
7
20
101438.94
0
3.04
1088.98
0.51
101438.94
0
8
22
4061.08
0
5.66
62.34
0.72
4061.08
0
9
24
91649.72
0
9.11
3537.69
0.77
91649.72
0
10
25
5097.69
0
9.48
123.33
0.92
5097.69
0
11
30
81299.10
0
22.71
1124.68
1.58
81299.10
0
12
35
41773.18
0
60.05
916.62
2.49
41773.18
0
13
40
113464.73
0
88.56
105.81
3.78
113464.73
0
14
45
75417.12
0
177.05
443.64
5.28
75417.12
275
15
50
42249.56
17816
23009.00
481.87
16
60
*
*
*
64.70
Min.
4061.08
0
0.33
62.34
Max.
165378.94
17816
23009
11550.21
Ave.
na
na
na
1620.53
1.29
1361.40
0.37
76.73
76.73
1.95
1665.67
0.42
97.20
94.44
2.76
2041.29
0.52
98.93
97.99
5.51
62.34
0.67
98.46
98.46
8.35
13968.51
0.81
96.14
84.76
9.02
2650.94
0.92
97.58
48.00
21.35
3423.67
1.55
98.62
95.79
59.32
1755.49
2.63
97.81
95.80
89.20
105.81
3.69
99.91
99.91
192.60
1292.91
5.46
99.41
98.29
AN US
‘*’: Instance not solved.
CR IP T
#
7.88
42249.56
2555
9071.78
899.27
11.04
98.85
97.87
14.71
47108.72
27539
13640.87
307.19
15.15
*
99.35
0.21
4061.08
0
0.37
62.34
0.25
76.73
46.79
14.71
165378.94
27539
13640.87
29432.11
15.15
99.91
99.91
2.52
56277.64
1898.06
1444.16
4915.41
2.77
na
86.61
Table 3: Numerical results for models M IP2 and M IP2V i using a transmission range of 350
M
m.
use generalized sub-tour elimination constraints as suggested in [11]. With this,
ED
the resulting linear relaxed solution values are considerably tighter than those of the models DTc , M IP1 and M IP2 . This is confirmed by the corresponding integrality gaps of these models. In Tables 4-6, we observe similar trends as for the Tables 1-3. We solve
PT
475
instances to optimality with up to 120 nodes with DTcV i and M IP1V i , while
CE
M IP2V i was able to solve instances only with up to 40 nodes. The average density of the instances with transmission range of 150 m in Tables 4-6 is 23.69%. In Tables 1-3, the average density is 74.34% for a transmission range of 350 m. In general, we observe that for sparser graphs in Tables 4-6, all the proposed
AC
480
models are harder to solve to optimality than for the instances in Tables 13. In particular, we do not solve the instances 13-22 with DTc and M IP1 in Tables 4-5 since the performance of these models deteriorate rapidly in terms of CPU time and number of CPLEX branch and bound nodes with respect to
ACCEPTED MANUSCRIPT
#
|V |
Den %
1
10
2
M odelII = DTcV i
M odelI = DTc
Integrality gaps
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
40.0
23092.52
0
0.16
18034.58
0.17
23092.52
0
0.18
20072.22
0.16
21.90
13.08
12
27.3
48404.68
1608
0.31
15473.25
0.17
48404.68
39
0.17
38899.51
0.16
68.03
19.64
3
14
31.9
49313.37
2996
0.81
12285.70
0.17
49313.37
0
0.19
44603.10
0.19
75.09
9.55
4
15
31.4
30541.11
1161
0.33
6222.58
0.18
30541.11
0
0.17
28662.80
0.16
79.63
6.15
5
16
23.3
55752.72
3322
1.01
10735.29
0.17
55752.72
6
18
21.6
81463.64
7774
2.17
18346.67
0.17
81463.64
7
20
32.1
24496.14
9014
1.93
3435.11
0.17
24496.14
8
22
18.2
84785.20
1640
0.47
29765.84
0.20
84785.20
9
24
19.2
80526.73
80956
19.13
27545.97
0.19
80526.73
10
25
21.0
81788.71
71415
18.38
11574.38
0.17
81788.71
11
30
23.4
68319.99
3609789
1497.94
14513.16
0.19
68319.99
12
35
20.8
62317.15
5407920
2673.93
4679.39
0.18
62317.15
13
40
23.2
*
*
*
8068.38
0.22
49488.68
14
45
17.7
*
*
*
4691.96
0.26
75328.01
15
50
21.6
*
*
*
2501.11
0.21
63502.41
16
60
21.4
*
*
*
2057.92
0.23
17
70
20.8
*
*
*
1104.35
18
80
21.3
*
*
*
19
90
20.3
*
*
20
100
22.3
*
21
110
21.2
22
120
CR IP T
Opt
0.19
51464.70
0.20
80.74
7.69
92
0.20
69614.28
0.17
77.48
14.55
0
0.18
20395.01
0.18
85.98
16.74
0
0.20
74023.19
0.19
64.89
12.69
43
0.25
64755.40
0.20
65.79
19.59
95
0.27
74403.74
0.18
85.85
9.03
673
1.28
52705.93
0.23
78.75
22.85
1237
0.84
50575.31
0.20
92.49
18.84
2145
1.32
42907.63
0.29
*
13.30
3315
4.59
57792.93
0.21
*
23.28
381682
571.11
41574.57
0.34
*
34.53
54085.65
490056
1272.71
37004.99
0.37
*
31.58
0.25
52190.98
147292
441.37
40132.70
0.45
*
23.10
1464.58
0.30
54560.08
20742
100.69
43144.40
0.52
*
20.92
*
1234.25
0.35
52221.99
1535848
12811.46
33703.59
0.64
*
35.46
*
*
472.70
0.34
39500.38
742084
13318.08
21983.80
4.58
*
44.35
*
*
*
838.56
0.44
46160.84
1068092
26345.97
29937.91
8.32
*
35.14
21.3
*
*
*
384.73
0.62
35746.04
97228
789.64
27000.90
0.87
*
24.46
Min.
17.70
23092.52
0
0.16
384.73
0.17
23092.52
0
0.17
20072.22
0.16
21.90
6.15
Max.
40
84785.20
5407920
2673.93
29765.84
0.62
84785.20
1535848
26345.97
74403.74
8.32
92.49
44.35
Ave.
23.69
na
na
na
8883.20
0.24
55163.04
204121.04
2530.04
43879.93
0.85
na
20.75
M
‘*’: Instance not solved.
AN US
0
485
ED
Table 4: Numerical results for models DTc and DTcV i using a transmission range of 150 m.
models DTcV i and M IP1V i . Similarly, we do not solve the instances 11-14 in
PT
Table 6 with M IP2 . Notice that the instance 14 in Table 6 cannot be solved to optimality in 10 hours. If we compare the average CPU time and number
CE
of CPLEX branch and bound (B&B) nodes required by DTcV i for the instances in Tables 1 and 4, we see an increment of 2530 seconds and of 204121 B&B
490
nodes, respectively. Similarly, if we compare the average CPU time and number
AC
of B&B nodes required by M IP1V i for the instances in Table 2 and 5, we see an increment of 117 seconds and of 313 B&B nodes, respectively. For M IP2V i ,
these values are also larger in Table 6 than in Table 3. In particular, we observe in Tables 4 and 5 that solving M IP1V i requires less CPU time and B&B nodes
495
than for DTcV i . With respect to the optimal solution values for these instances,
ACCEPTED MANUSCRIPT
M odelII = M IP1V i
M odelI = M IP1
Integrality gaps
|V |
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
1
10
23092.52
0
0.25
18034.58
0.18
23092.52
0
0.23
20072.22
0.18
21.90
13.08
2
12
49644.58
228
0.81
15473.25
0.19
49644.58
0
0.45
38899.51
0.19
68.83
21.64
3
14
49313.37
0
0.58
12285.70
0.20
49313.37
0
0.44
44603.10
0.23
75.09
9.55
4
15
37742.02
609
1.64
6222.58
0.23
37742.02
0
0.45
30479.76
0.24
83.51
19.24
5
16
65018.51
2670
3.76
10735.29
0.27
65018.51
0
6
18
105718.12
1846
4.41
18409.93
0.31
105718.12
0
7
20
28618.15
4354
13.29
3435.11
0.34
28618.15
0
8
22
90372.38
0
2.31
29804.10
0.41
90372.38
0
9
24
93457.94
8954
68.55
27565.03
0.45
93457.94
0
10
25
109801.97
12591
91.68
11574.38
0.59
109801.97
99
11
30
72915.68
298589
3098.10
14513.16
0.81
72915.68
211
12
35
65624.59
1954271
36000
4679.39
1.17
65624.59
0
13
40
*
*
*
8068.38
2.52
54640.40
51
14
45
*
*
*
4741.48
2.37
86892.79
1887
15
50
*
*
*
2501.11
16
60
*
*
*
2082.34
17
70
*
*
*
1104.35
18
80
*
*
*
1467.46
19
90
*
*
*
1247.48
20
100
*
*
*
472.70
21
110
*
*
*
844.37
22
120
*
*
*
386.83
Min.
23092.52
0
0.25
386.83
Max.
109801.97
1954271
36000
29804.10
Ave.
na
na
na
8893.13
51723.09
0.28
83.49
20.45
0.70
69614.28
0.30
82.59
34.15
1.10
20395.01
0.35
88.00
28.73
1.25
74023.19
0.40
67.02
18.09
2.23
65052.04
0.45
70.51
30.39
7.68
74414.69
0.53
89.46
32.23
13.31
53238.70
0.84
80.09
26.99
12.51
51020.42
1.33
92.86
22.25
25.73
43270.73
1.85
*
20.81
181.61
58056.61
2.80
*
33.19
AN US
0.55
3.27
74332.33
144
71.71
42282.05
4.14
*
43.12
5.92
79566.37
1424
246.78
37115.07
10.13
*
53.35
9.12
74225.29
1026
370.96
40653.33
24.82
*
45.23
13.75
96169.10
0
238.23
44211.81
51.48
*
54.03
19.78
119041.25
0
279.12
35418.41
67.32
*
70.25
28.92
105560.31
0
326.00
23133.77
176.51
*
78.08
40.62
110563.48
2044
779.15
30730.80
203.25
*
72.21
50.47
143280.95
21
415.24
27496.97
497.72
*
80.81
0.18
23092.52
0
0.23
20072.22
0.18
21.90
9.55
50.47
143280.95
2044
779.15
74414.69
497.72
92.86
80.81
8.26
78890.55
313.95
135.24
44359.34
47.51
na
37.63
M
‘*’: Instance not solved.
CR IP T
#
ED
Table 5: Numerical results for models M IP1 and M IP1V i using a transmission range of 150 m.
PT
we observe again the relations Opt(DTcV i ) ≤ Opt(M IP1V i ) ≤ Opt(M IP2V i ). However, in this case, we observe that the optimal solution values obtained with M IP2V i are larger than those obtained with M IP1V i for the most part
CE
of the instances in Tables 5 and 6. As before, all the proposed models which
500
use generalized sub-tour elimination constraints present smaller CPU times and
AC
B&B nodes than the other models, as well as tighter linear relaxed solution values. In Tables 7-9, we confirm the trends observed in Tables 1-6. In this case,
we solve instances with transmission range equal to 100 m and with up to
505
90, 120 and 30 nodes to optimality with models DTcV i , M IP1V i and M IP2V i , respectively. The average density for these instances is 15.23%. In this case, we
ACCEPTED MANUSCRIPT
M odelII = M IP2V i
M odelI = M IP2
Integrality gaps
|V |
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
1
10
31209.13
19
0.59
18034.58
0.22
31209.13
0
0.72
20072.22
0.21
42.21
35.68
2
12
49644.58
292
1.20
15473.25
0.23
49644.58
0
0.81
38899.51
0.22
68.83
21.64
3
14
58875.71
10
2.56
12285.70
0.28
58875.71
0
1.51
44603.10
0.27
79.13
24.24
4
15
50636.68
1637
7.57
6222.58
0.30
50636.68
17
2.20
30479.76
0.29
87.71
39.81
5
16
66653.52
1086
10.37
10735.29
0.31
66653.52
0
6
18
109854.11
1128
19.88
18350.30
0.37
109854.11
0
7
20
28618.15
2124
40.34
3435.11
0.48
28618.15
0
8
22
119731.92
4415
103.47
29765.84
0.58
119731.92
965
9
24
93457.94
37747
946.02
27545.97
0.80
93457.94
0
10
25
120619.93
3951
250.34
11574.38
1.00
120619.93
141
11
30
*
*
*
14513.16
1.32
74851.40
1037
12
35
*
*
*
4679.39
2.05
69707.63
1147
13
40
*
*
*
8068.38
3.06
68765.65
40751
14
45
*
*
*
4701.45
5.96
-
16585
Min.
28618.15
10
0.59
3435.11
Max.
120619.93
37747
946.02
29765.84
Ave.
na
na
na
13241.81
1.61
51493.66
0.34
83.89
22.74
2.71
69614.28
0.39
83.30
36.63
5.48
20395.01
0.49
88.00
28.73
27.71
74023.19
0.64
75.14
38.18
19.72
64899.85
0.75
70.53
30.56
28.54
74404.22
0.89
90.40
38.32
204.77
53035.09
1.68
*
29.15
324.69
50591.14
2.74
*
27.42
31901.30
43031.95
4.46
*
37.42
36000
57792.92
5.18
-
-
AN US
‘*’: Instance not solved.
CR IP T
#
0.22
28618.15
5.96
120619.93
1.21
na
0
0.72
20072.22
0.21
42.21
21.64
40751
36000
74404.22
5.18
90.40
39.81
4331.64
4894.41
49523.99
1.32
na
na
‘-’: No solution found by CPLEX in 10 hours of execution time.
Table 6: Numerical results for models M IP2 and M IP2V i using a transmission range of 150
M
m.
do not solve the instances 12-22 with models DTc and M IP1 in Tables 7-8 and do not solve the instances 12-14 in Table 9 neither with model M IP2 nor with model
510
ED
M IP2V i . For the latter, these instances proved to be infeasible which evidences the more restrictive nature of the joint probabilistic constraints. In particular,
PT
notice that for the instance 11 in Table 9, model M IP2 cannot obtain a feasible solution in 10 hours while M IP2V i allows to obtain the optimal solution of this instance in less than 32 seconds. This shows the remarkable effectiveness of using
CE
generalized sub-tour elimination constraints in M IP2 . The remaining instances
515
in Table 9 are solved to optimality in less than 13 seconds. The maximum CPU time and number of B&B nodes required by DTcV i for the instances in Table
AC
4 are 26345.97 seconds and 1535848 B&B nodes, respectively. While these values in Table 7 are 10 hours and 3283374 B&B nodes, respectively. Similarly, the maximum CPU time and number of B&B nodes required by M IP1V i for
520
the instances in Tables 5 and 8 are 779.15 seconds and 2044 B&B nodes, and 6842.07 seconds and 7316 B&B nodes, respectively. Once again, we see that
ACCEPTED MANUSCRIPT
#
|V |
Den %
1
10
2
M odelII = DTcV i
M odelI = DTc
Integrality gaps
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
24.4
37059.44
109
0.18
12988.63
0.17
37059.44
0
0.16
32670.26
0.16
64.95
11.84
12
24.2
29210.83
1743
0.34
6144.18
0.17
29210.83
0
0.19
25446.84
0.17
78.97
12.89
3
14
18.7
43698.31
859
0.28
22948.40
0.19
43698.31
0
0.17
36559.41
0.17
47.48
16.34
4
15
30.5
15058.40
132
0.26
6452.20
0.18
15058.40
0
0.19
14399.43
0.16
57.15
5
16
28.3
37183.35
176606
23.74
7988.78
0.17
37183.35
6
18
19.0
51756.80
37094
5.55
8945.88
0.16
51756.80
7
20
19.5
37129.01
19444
3.78
14016.64
0.18
37129.01
8
22
12.6
53246.64
1907
0.44
24492.05
0.17
53246.64
CR IP T
Opt
0.76
26066.52
0.17
78.52
29.90 8.77
0
0.17
47217.18
0.17
82.72
53
0.20
33568.69
0.18
62.25
9.59
5
0.20
42900.75
0.19
54.00
19.43
140
0.27
68481.14
0.19
69.78
13.12
125
0.23
51847.69
0.18
68.20
17.13
0
0.34
56431.89
0.18
79.51
14.63
3506
1.75
69276.38
0.24
*
26.10
4
0.62
58193.85
0.42
*
14.81
1722
0.97
72773.23
0.20
*
13.46
24
13.4
78826.94
228457
43.12
23817.61
0.17
78826.94
25
13.3
62563.97
12894
3.02
19895.33
0.17
62563.97
11
30
11.7
66102.25
141751
40.52
13543.69
0.20
66102.25
12
35
11.4
*
*
*
14587.62
0.18
93746.39
13
40
13.1
*
*
*
11592.88
0.33
68309.31
14
45
9.9
*
*
*
17672.45
0.19
84094.27
15
50
13.5
*
*
*
9096.80
0.20
61561.21
15843
14.48
46175.67
0.22
*
24.99
16
60
10.1
*
*
*
8056.37
0.26
81849.38
1842122
1529.41
59355.98
0.36
*
27.48
17
70
9.9
*
*
*
9542.68
0.25
90175.30
804550
1065.83
65187.84
0.39
*
27.71
18
80
10.7
*
*
*
7484.63
0.25
69333.96
3283374
36000
50158.65
0.70
*
27.66
19
90
10.3
*
*
*
5790.72
0.30
65868.76
865667
4886.80
53093.35
0.55
*
19.40
20
100
10.1
*
*
*
5111.89
0.33
77769.10
2565329
36000
61545.34
1.26
*
20.86
21
110
10.3
*
*
*
3552.84
0.41
60757.72
2193242
36000
45797.55
1.33
*
24.62
22
120
10.2
*
*
*
3136.19
1.01
-
-
-
41738.95
2.87
*
-
Min.
9.90
15058.40
109.00
0.18
3136.19
0.16
15058.40
0
0.16
14399.43
0.16
47.48
4.38
Max.
30.50
78826.94
228457
43.12
24492.05
1.01
93746.39
3283374
36000
72773.23
2.87
82.72
29.90
Ave.
15.23
na
na
na
11675.38
0.25
na
na
na
48131.20
0.47
na
na
M
AN US
9 10
‘*’: Instance not solved.
ED
‘-’: No solution found by CPLEX due to a shortage of memory in at most 10 hours of execution time.
Table 7: Numerical results for models DTc and DTcV i using a transmission range of 100 m.
PT
solving M IP1V i is easier than solving DTcV i . Observe that the optimal solution values obtained with M IP2V i are larger than those obtained with M IP1V i for many instances in Tables 8 and 9. In order to give more insight with respect to the parameter values of θ ∈ [1, ∞), in Figure 4 we plot curves for the objective
CE 525
function values of M IP2V i and LP2V i for values of θ ∈ [1, 30]. We choose this
AC
interval as we noticed in our numerical results that for large values of θ, these curves remain constant. Consequently, it would suffice to solve M IP1V i instead of M IP2V i for larger values of θ as they are equivalent models. By doing so, we
530
4.38
3137
cover all range of optimal solution values for θ ∈ [1, ∞). In other words, we cover joint and separate probabilistic constraints in a unified framework. In Figure 4, we randomly generate 20 instances with dimensions of |V | ∈ {20, 30} while using
ACCEPTED MANUSCRIPT
M odelII = M IP1V i
M odelI = M IP1
Integrality gaps
|V |
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
1
10
37135.50
0
0.24
12988.63
0.19
37135.50
0
0.26
32670.26
0.19
65.02
12.02
2
12
29210.83
311
0.62
6144.18
0.20
29210.83
0
0.31
25446.84
0.20
78.97
12.89
3
14
46506.67
0
0.51
23066.16
0.22
46506.67
0
0.39
37640.09
0.23
50.40
19.07
4
15
17534.73
1366
2.89
6452.20
0.24
17534.73
0
0.47
14399.43
0.24
63.20
17.88
5
16
37371.33
6155
11.37
7988.78
0.27
37371.33
206
6
18
53950.88
350
3.35
8945.88
0.30
53950.88
0
7
20
39473.32
6800
22.24
14016.64
0.42
39473.32
0
8
22
54899.55
21
2.56
24531.80
0.38
54899.55
0
9
24
81104.56
28899
135.35
23817.61
0.51
81104.56
11
10
25
62563.97
3374
22.88
19895.33
0.50
62563.97
0
11
30
84006.04
2935324
36000
13543.69
0.85
84006.04
0
12
35
*
*
*
14587.62
1.20
102660.28
7316
13
40
*
*
*
11592.92
1.86
83699.46
808
14
45
*
*
*
17787.13
2.59
100856.63
0
15
50
*
*
*
9096.80
16
60
*
*
*
8077.72
17
70
*
*
*
9589.57
18
80
*
*
*
7511.77
19
90
*
*
*
5805.59
20
100
*
*
*
5177.71
21
110
*
*
*
3594.44
22
120
*
*
*
3148.52
Min.
17534.73
0
0.24
3148.52
Max.
84006.04
2935324
36000
24531.80
Ave.
na
na
na
26106.59
0.25
78.62
30.14
0.62
47469.20
0.28
83.42
12.01
0.99
33783.35
0.34
64.49
14.41
1.05
42900.75
0.41
55.32
21.86
3.96
68568.85
0.45
70.63
15.46
2.50
51849.05
0.52
68.20
17.13
3.83
56750.11
0.88
83.87
32.45
186.99
70581.10
1.35
*
31.25
49.44
59177.87
2.03
*
29.30
25.85
72947.29
2.54
*
27.67
AN US
1.12
3.18
82085.96
1089
211.40
47330.98
3.70
*
42.34
5.44
94646.26
4341
1155.25
60609.40
8.11
*
35.96
8.41
113202.67
2628
1308.36
68833.86
11.92
*
39.19
12.59
115524.99
63
513.74
51738.12
40.11
*
55.21
18.02
97746.86
4828
4277.56
54613.25
59.67
*
44.13
24.73
135766.52
1243
1274.63
62456.70
82.06
*
54.00
34.09
84259.69
3550
5918.51
46655.89
157.72
*
44.63
44.66
119676.87
5766
6842.07
45150.42
189.31
*
62.27
0.19
17534.73
0
0.26
14399.43
0.19
50.40
12.01
44.66
135766.52
7316
6842.07
72947.29
189.31
83.87
62.27
7.31
76085.61
1447.68
989.96
48985.42
25.56
na
30.51
M
‘*’: Instance not solved.
11698.21
CR IP T
#
ED
Table 8: Numerical results for models M IP1 and M IP1V i using a transmission range of 100 m.
PT
transmission ranges equal to 100 m and 150 m. Thus, we compute the average optimal solution values of the MILP model and its linear relaxation. We do not 535
consider large size instances because the trends of the curves remain the same
CE
and also because the CPU times become rapidly prohibitive. From Figure 4, we observe that the objective function values of M IP2V i are considerably larger
AC
for values of θ ∈ [1, 10]. Finally, we observe that the linear relaxation remains constant for all these instances.
540
In Table 10, we present numerical results for 18 benchmark instances [4, 26]
with dimensions of |V | ∈ {50, 100} and transmission ranges equal to 150 m, 125
m and 100 m. We report results only for DTcV i and M IP1V i . We do not present
numerical results for M IP2V i for these instances because CPLEX runs out of
ACCEPTED MANUSCRIPT
M odelII = M IP2V i
M odelI = M IP2
Integrality gaps
|V |
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
Opt
B&Bn
T ime(s)
LR
T imeLR (s)
GapI %
GapII %
1
10
37135.50
0
0.42
12988.63
0.20
37135.50
0
0.39
32670.26
0.21
65.02
12.02
2
12
30341.18
237
1.28
6144.18
0.23
30341.18
0
0.62
25446.84
0.23
79.75
16.13
3
14
46506.67
0
1.01
23066.16
0.28
46506.67
0
0.83
37640.09
0.25
50.40
19.07
4
15
17534.73
96
2.80
6452.20
0.29
17534.73
0
1.51
14399.43
0.30
63.20
17.88
5
16
37371.33
4393
27.08
7988.78
0.34
37371.33
157
5.90
26071.30
0.33
78.62
30.24
6
18
56179.84
1749
16.85
8945.88
0.41
56179.84
3
4.34
47469.20
0.39
84.08
15.50
7
20
41537.16
7481
110.38
14016.64
0.47
41537.16
3
7.91
33777.43
0.65
66.26
18.68
8
22
54899.55
111
10.58
24531.80
0.62
54899.55
0
2.62
42900.75
0.56
55.32
21.86
9
24
81104.56
31519
450.44
23817.61
0.75
81104.56
43
10.78
68499.04
0.73
70.63
15.54
10
25
66078.65
2969
98.75
19895.33
0.97
66078.65
38
12.15
52077.55
0.90
69.89
21.19
11
30
-
447732
36000
13543.69
1.31
84090.22
5
31.40
56538.15
1.43
-
32.76
12
35
Inf
Inf
Inf
12988.62
0.33
Inf
Inf
Inf
32670.26
0.23
Inf
Inf
13
40
Inf
Inf
Inf
11630.40
3.03
Inf
Inf
Inf
58843.38
3.77
Inf
Inf
14
45
Inf
Inf
Inf
17709.88
4.13
Inf
Inf
Inf
72884.68
5.32
Inf
Inf
Min.
17534.73
0
0.42
6144.18
Max.
81104.56
447732
36000
24531.80
Ave.
na
na
na
14551.41
‘Inf’: Infeasible solution.
AN US
‘-’: No solution found with CPLEX in 10 hours.
CR IP T
#
0.20
17534.73
0
0.39
14399.43
0.21
50.40
12.02
4.13
84090.22
157
31.40
72884.68
5.32
84.08
32.76
0.95
na
na
na
42992.02
1.09
na
na
Table 9: Numerical results for models M IP2 and M IP2V i using a transmission range of 100
4
M
x 10
6
MIP2V i
ED
Optimal solutions
7
LP2V i
5
4
10 15 20 25 Values of θ with |V|=20 and using a transmission range of 100m .
9
LP2V i 7
Optimal solutions
Optimal solutions
CE 5.5
AC
6
7
6
5
5
10 15 20 25 Values of θ with |V|=20 and using a transmission range of 150m .
30
4
x 10
6.5
MIP2V i
8
4
7.5
x 10
6
30
PT
5
4
10
Optimal solutions
m.
MIP2V i LP2V i
10 15 20 25 Values of θ with |V|=30 and using a transmission range of 100m .
30
x 10
5.5 MIP2V i LP2V i
5
4.5
5
10 15 20 25 Values of θ with |V|=30 and using a transmission range of 150m .
Figure 4: Objective function values for M IP2V i and LP2V i with θ ∈ [1, 30].
30
ACCEPTED MANUSCRIPT
#
Den. %
M odelI = DTcV i Opt
B&Bn
T ime(s)
LR
M odelII = M IP1V i T imeLR (s)
Opt
B&Bn
T ime(s)
LR
Integrality gaps T imeLR (s)
GapI %
GapII %
ins-050-1
22.5
647.75
930
2.36
569.33
0.29
694.59
0
57.03
569.68
4.72
12.11
17.98
ins-050-2
21.4
863.69
2147
7.32
751.17
0.29
907.24
43
66.68
751.17
4.94
13.03
17.20
ins-050-3
22.0
743.94†
819
2.84
659.58
0.29
793.46
9
60.44
659.58
4.21
11.34
16.87
ins-100-1
21.0
876.69
148438
2085.21
755.25
0.78
1020.90
17
610.92
755.25
140.19
13.85
26.02
ins-100-2
21.8
657.35
36441
681.26
576.02
0.79
745.14
0
545.82
576.19
124.98
12.37
22.67
ins-100-3
20.5
722.87
3057
98.61
636.80
0.68
884.43
CR IP T
Benchmark instances with |V | = 50 and |V | = 100 and transmission range of 150 m.
82
689.12
636.83
135.62
11.91
27.99
Benchmark instances with |V | = 50 and |V | = 100 and transmission range of 125 m. 15.8
802.95
2677
5.41
696.11
0.26
835.45
135
69.36
696.13
4.00
13.31
16.68
ins-050-2
15.7
1055.10
101943
136.01
899.01
0.33
1122.71
3356
601.42
899.16
4.53
14.79
19.91
ins-050-3
15.3
877.77
13144
23.36
764.52
0.34
914.92
411
99.06
764.65
4.75
12.90
16.42
ins-100-1
16.1
943.01
92870
1087.90
839.79
0.82
1018.00
61
1214.40
839.79
129.58
10.95
17.51
ins-100-2
15.8
917.00
1229927
22802.25
754.12
1.06
963.19
198
1261.36
754.12
69.93
17.76
21.71
ins-100-3
14.7
998.18
542634
9003.54
839.74
0.66
1077.58
413
1760.61
839.74
111.58
15.87
22.07
AN US
ins-050-1
Benchmark instances with |V | = 50 and |V | = 100 and transmission range of 100 m. ins-050-1
10.0
1204.41
2985
3.93
1101.23
0.36
1300.09
6181
591.84
1101.23
3.75
8.57
15.30
ins-050-2
9.6
1340.44
249177
136.92
1099.81
0.40
1386.78
18091
1464.24
1101.46
3.82
17.95
20.57
ins-050-3
10.3
1316.39
374158
257.42
1131.34
0.38
1426.21
8610
834.99
1131.62
4.89
14.06
20.66
ins-100-1
10.8
1217.47
2020217
25155.92
1008.94
2.28
1256.97
3152
11689.91
1008.94
82.98
17.13
19.73
ins-100-2
10.6
1128.40
795701
6789.46
945.70
0.72
1204.59
11402
36000
948.02
50.17
16.19
21.29
ins-100-3
9.7
1252.99
1594116
13866.27
1055.37
0.72
1358.58
3355
36000
1055.64
78.13
15.77
22.30
‘†’: In [26] the reported optimal solution value is 743.74. The correct is really 743.94.
M
Table 10: Numerical results for Benchmark instances [26].
memory rapidly. From Table 10, we confirm that our proposed deterministic compact model DTcV i can solve all these instances to optimality in less than
ED
545
7 hours of CPU time. Consequently, we further confirm that the solutions
PT
obtained with the meta-heuristic approaches proposed by the authors in [4, 26] are the optimal ones (with the exception of instance ins-050-3 with transmission range of 150 m). These values were not proven optimal in [26]. In particular, we observe that the CPU time required by CPLEX to solve DTcV i increases when
CE 550
the density of the graphs decreases. Regarding the optimal solution values
AC
of M IP1V i , we observe that they are larger than those obtained with DTcV i . Conversely, we observe that the CPU times are lower except for the instances ins100-2 and ins-100-2 with transmission range of 100 m that CPLEX cannot solve
555
with M IP1V i in 10 hours of execution time. Again, we observe that the CPU times increase when the density of the graphs is smaller. Finally, we observe that the CPU time to solve the linear relaxations is significantly larger for M IP1V i
ACCEPTED MANUSCRIPT
than for DTcV i , whereas the number of B&B nodes is significantly larger for the latter model. The integrality gaps confirm that the linear relaxation of DTcV i is tighter than the one of M IP1V i .
CR IP T
560
6. Conclusion
In this paper, we introduce a new probabilistically constrained dominating
tree problem which consists in finding a minimum cost dominating tree T of
an undirected sensor network G subject to a probabilistic constraint related 565
to the risk that the total power consumption of each sensor in T exceeds a
AN US
given limit. To handle uncertain sensor power consumption of the links of the network, we consider both separate and joint probabilistic constraints in an im-
proved compact polynomial formulation for the classic dominating tree problem. The dependence among power sensor consumption of the links of the network is 570
handled by means of copulas [6, 13]. The compact formulation is based on ar-
M
borescences of directed graphs for which we give a proof of its correctness. The compact model allows to obtain deterministic equivalent mixed integer linear programming formulations for the probabilistic counterpart. Consequently, we
575
ED
solve random and challenging benchmark instances for the DT and PCDT problems with up to 120 nodes to optimality. Solving these problems to optimality
PT
is extremely important as it allows to compare the efficiency of any other algorithmic procedure to be proposed as part of future research. Additionally, new modeling approaches while considering the overhead and other metrics of sensor
CE
networks should be considered as variants of the DT and PCDT problems in
580
order to generate new protocols. Finally, it is important to note that the probabilistic modelling approach allows to handle joint and separate probabilistic
AC
constraints in a unified framework.
Acknowledgments The authors are grateful to the five referees and to the Editor Professor
585
Stavros Toumpis for their valuable comments and suggestions that helped to
ACCEPTED MANUSCRIPT
improve significantly our paper. The first author also acknowledges the financial support of the USACH/DICYT Project 061513VC DAS.
References 590
CR IP T
References
[1] P. Adasme and A. Lisser, Stochastic and semidefinite optimization for
scheduling in orthogonal frequency division multiple access networks, Journal of Scheduling 17(5), 445-469, (2014).
AN US
[2] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004. 595
[3] A. Charnes, W. Cooper and G. Symonds, Cost horizons and certainty equivalents: an approach to stochastic programming of heating oil, Management Science, vol. 4, 235-263, (1958).
M
[4] S. N. Chaurasia and A. Singh, A hybrid heuristic for dominating tree problem, Soft Computing, 20 (1), 377-397, (2016). [5] J. Cheng, E. Delage and A. Lisser, Distributionally Robust Stochas-
ED
600
tic Knapsack Problem, SIAM Journal on Optimization 24(3), 1485-1506
PT
(2014).
[6] J. Cheng, M. Houda and A. Lisser, Chance constrained 0–1 quadratic pro-
CE
grams using copulas, Optimization Letters 9(7), 1283-1295, (2015). 605
[7] J. Cheng and A. Lisser, A second-order cone programming approach for
AC
linear programs with joint probabilistic constraints, Operations Research Letters, vol. 40, 325-328, (2012).
[8] J. Carle and D. Simplot-Ryl, Energy-efficient area monitoring for sensor
610
networks, Computer, 37(2), 40-46, (2004).
[9] R. Fortet, Applications de l’alg`ebre de boole en recherche operationelle, Revue Francaise de Recherche Operationelle, 4, 17-26, (1960).
ACCEPTED MANUSCRIPT
[10] A. Gaivoronski, A. Lisser, R. Lopez and X. Hu, Knapsack problem with probability constraints, Journal of Global Optimization, 49, 397-413, (2011). [11] B. Gendron, A. Lucena, A. Salles da Cunha and L. Simonetti, Benders
CR IP T
615
Decomposition, Branch-and-Cut, and Hybrid Algorithms for the Minimum
Connected Dominating Set Problem, INFORMS Journal on Computing 26(4), 645-657, (2014).
[12] A. Gupte, S. Ahmed, M. Cheon and S. Dey, Solving Mixed Integer Bilinear
AN US
Problems Using MILP Formulations, SIAM Journal on Optimization 23(2),
620
721-744, (2013).
[13] R. Henrion and C. Strugarek, Convexity of chance constraints with dependent random variables: the use of copulae. In: Bertocchi M., Consigli G., Dempster M.A.H. (eds.), Chapter: Stochastic Optimization Methods in Finance and Energy Vol. 163 of the series International Series in Operations
625
M
Research & Management Science, 427-439, Springer, New York (2011).
ED
[14] IBM ILOG CPLEX Optimization Studio Information Center. [15] M. Lobo, L. Vandenberghe, S. Boyd and H. Lebret, Linear Algebra and its Applications, 284, 193-228, (1998). [16] J. Luedtke, S. Ahmed and G. Nemhauser, An integer programming ap-
PT
630
proach for linear programs with probabilistic constraints, Mathematical
CE
Programming, 122, 247-272, (2010). [17] L. Miller and H. Wagner, Chance-constrained programming with joint con-
AC
straints, Operations Research, 13, 930-945, (1965).
635
[18] R. Nelsen, An Introduction to Copulas, 2nd edn. Springer, New York, (2006).
[19] A. Pr´ekopa, On probabilistic constrained programming, Proceedings of the Princeton Symposium on Mathematical Programming, Princeton University Press, Princeton, (1970).
ACCEPTED MANUSCRIPT
640
[20] A. Pr´ekopa, Stochastic Programming, 600. Kluwer Academic Publishers, Dordrecht Boston, (1995). [21] N. Sahinidis, Optimization under uncertainty: State-of-the-art and oppor-
CR IP T
tunities, Computers and Chemical Engineering, 28, 971-983, (2004).
[22] R. Schultz, L. Stougie and M. van der Vlerk, Two-stage stochastic integer programming: A survey, Statistica Neerlandica, 50, 404-416, (1996).
645
[23] A. Shapiro, D. Dentcheva and A. Ruszczynski, Lectures on stochastic programming: Modeling and theory, 436. SIAM Philadelphia, Series on Opti-
AN US
mization, Vol. 9 of MPS/SIAM, Philadelphia, (2009).
[24] I. Shin, Y. Shen and M. Thai, On approximation of dominating tree in wireless sensor networks, Optimization Letters 4, 393-403, (2010).
650
[25] V. Singh, S. Jain and A. Tyagi, Risk and reliability analysis: a handbook
M
for civil and environmental engineers, vol. 1. ASCE, (2007). [26] S. Sundar and A. Singh, New heuristic approaches for the dominating tree
AC
CE
PT
ED
problem, Applied Soft Computing 13, 4695-4703, (2013).
ACCEPTED MANUSCRIPT
Pablo Adasme is an associate researcher and full professor in computer sci-
655
ence at the Electrical Engineering Department of the Universidad de Santiago de Chile. He received the PhD degrees from the Universidad de Santiago de
CR IP T
Chile and from the University of Paris Sud 11, in 2010 in Chile and France, respectively. Currently, his main research interests are stochastic and semidefi660
nite programming, combinatorial optimization problems, and metaheuristic al-
PT
ED
M
AN US
gorithms applied to wireless communication networks and energy problems.
Rafael Andrade received the B.Sc. degree in Computer Science from State University of Cear (UECE, Fortaleza, Brazil) in 1997 and in 1999, the M.Sc. degree in Computer Science and System Engineering from Federal University of
CE 665
Rio de Janeiro (UFRJ, Rio de Janeiro, Brazil). From 1999 to 2002 he realized his
AC
Ph.D. thesis at France Tlcom R&D laboratory (FTR&D, Issy-les-Moulineaux, France). In 2002, he obtained his Ph.D. degree in Informatics from Paris-Nord University (Paris 13, Villetaneuse, France). He was assistant professor at Paris-
670
Sud University and at Paris-Nord University, in the Laboratory for Computer Science (LRI) from 2002 to 2003, and in the Computer Science Laboratory of Paris-Nord (LIPN) from 2003 to 2004, respectively. In the occasion of a sabbat-
ACCEPTED MANUSCRIPT
ical year in 2012, he was invited professor at LRI (Paris-Sud). Currently, he is Associate Professor at Federal University of Cear (UFC, Fortaleza, Brazil) in the 675
department of statistics and applied mathematics where he has been working
CR IP T
since 2004. His main research interests are stochastic, robust and deterministic programming, operations research, combinatorial optimization, telecommunica-
ED
M
AN US
tion network design and network optimization problems.
Abdel Lisser is a full Professor in Computer Science laboratory of University
PT
680
of Paris Sud since 2001. He was heading research group at France Telecom from Research Center 1996 to 2001 and research engineer at France Telecom
CE
Research Center from 1988 to 1996. He got the Ph.D. at the University of Paris Dauphine in Computer Science in 1987 and the Habilitation thesis at the University of Paris Nord in 2000. His main research area is combinatorial
AC
685
and stochastic optimization with application to telecommunication and recently energy problems.
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT