Modeling and solving the angular constrained minimum spanning tree problem

Computers and Operations Research 112 (2019) 104775

Contents lists available at ScienceDirect

Computers and Operations Research journal homepage: www.elsevier.com/locate/cor

Modeling and solving the angular constrained minimum spanning tree problem Alexandre Salles da Cunha a,1,∗, Abilio Lucena b,2 a b

Universidade Federal de Minas Gerais, Departamento de Ciência da Computação, Belo Horizonte, Brazil Universidade Federal do Rio de Janeiro, Programa de Engenharia de Sistemas e Computação, Rio de Janeiro, Brazil

a r t i c l e

i n f o

Article history: Received 27 December 2018 Revised 17 August 2019 Accepted 20 August 2019 Available online 20 August 2019 Keywords: Combinatorial optimization Angular constrained spanning trees Branch-and-cut algorithms Polyhedral combinatorics

a b s t r a c t Assume one is given an angle α ∈ (0, 2π ] and a complete undirected graph G = (V, E ). The vertices in V represent points in the Euclidean plane. The edges in E represent the line segments between these points, with edge weights equal to segment lengths. Spanning trees of G are called α -spanning trees (α -STs) if, for any i ∈ V, the smallest angle that encloses all line segments corresponding to its i-incident edges does not exceed α . The Angular Constrained Minimum Spanning Tree Problem (α -MSTP) seeks an α -ST with the least weight. The problem arises in the design of wireless communication networks operating under directional antennas. We propose two α -MSTP formulations. One, Fx requiring, in principle, O(2|V| ) inequalities to model the angular constraints (α -AC). For α ∈ (0, π ), however, we show that just O(|V|3 ) of them suffice to attain not only a formulation but also the same Linear Programming relaxation (LPR) bound as the full blown model. The other formulation, Fxy , enlarges the set of Fx variables but manages to model α -AC, compactly. Furthermore, Fxy LPR bounds are proven to dominate their Fx counterparts. That dominance, however, is empirically shown to reduce as α increases. Finally, exact Branch-and-Cut algorithms implemented for the two formulations are shown, empirically, to exchange roles as top performer, throughout the [0, 2π ) range of α values. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction We investigate the Angular Constrained Minimum Spanning Tree Problem (α -MSTP). It is defined over a complete undirected graph G = (V, E ), with n = |V | vertices and m = |E | edges, and takes an angle α ∈ (0, 2π ] as an input parameter. Every vertex of V corresponds to a point in the Euclidean plane and, for simplicity, we simultaneously denote by i a vertex of V and the point it represents in the Euclidean plane. An edge e = {i, j} ∈ E, in turn, corresponds to the line segment defined by points i and j and has a weight equal to the Euclidean distance we , between the two. The weight of a spanning tree equals the sum of the weights for its edges. A spanning tree (V, ET ) of G is an α -spanning tree (α -ST) if, for every vertex i ∈ V, the smallest angle that encloses all line segments

∗

Corresponding author. E-mail addresses: [email protected] (A.S. da Cunha), [email protected] (A. Lucena). 1 Alexandre Salles da Cunha was partially funded by CNPq grant 303928/2018-2 and FAPEMIG grant CEX - PPM-00164/17. 2 Abilio Lucena was partially funded byCNPq Grant 309887/2018-6. https://doi.org/10.1016/j.cor.2019.104775 0305-0548/© 2019 Elsevier Ltd. All rights reserved.

corresponding to its i-incident edges does not exceed α . Among all α −STs, α −MSTP looks for one with minimum weight. The most important defining property for the problem is its angular constraints (α -AC). To better illustrate them, take α = π2 and assume that all i-incident edges for a given graph G are depicted in Fig. 1. As complementary information, the horizontal and vertical coordinates for all corresponding Euclidean plane points are shown in Table 1. From the information provided, edges {i, u} and {i, z} cannot simultaneously belong to a π2 -ST of G. To check that,  and iz  denote the (Euclidean norm) unit vectors diassume that iu rected from point i towards respectively points u and z. For the example in Fig. 1, all unit vectors for i-incident edges are illustrated in Fig. 2. From that figure, one should then note that the circu anti-clockwise lar sector centered at i and obtained by rotating iz  , has a central angle of 7π radians. Likearound i until it meets iu 12  antiwise, the circular sector centered at i, obtained by rotating iu  clockwise around i until it meets iz, has a central angle of 1712π radians. Since the minimum of these two central angles exceeds α = π2 , no circular sector centered at i and angled π2 radians or  and iu  . Notice, howless would simultaneously enclose vectors iz π , a spanever, that if the admissible angle α were enlarged to 712 ning tree of G with {i, z} and {i, u} as its only i-incident edges would not violate α -AC for i. Likewise, if {i, z}, {i, q} and {i, u}

2

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775 Table 1 Input data for points (i.e., vertices) in Figs. 1–3. Coordinates for the Euclidean plane points Vertex of V

Horizontal

Vertical

ࢭ0ip

i z q u t l

0 3 1 −1 −1 2

0 √ √3 3 1 −1 −2

-

π

particular, they are more energy efficient and reduce network congestion, as well as signal interference. Details on the use of directional antennas for wireless networks are found, among others, in Carmin et al. (2011), Ackerman et al. (2013), Yu et al. (2014). 1.1. Our contribution

6

π

3 3π 4 5π 4 7π 4

Fig. 1. Euclidean plane points i, l, q, t, u and z (corresponding coordinates shown in Table 1) and line segments connecting i to every other point.

Fig. 2. Unit vectors corresponding to the line segments (i.e., edges) in Fig. 1.

were the only spanning tree edges incident to i, a similar observation would also apply. The α -MSTP was introduced by Aschner and Katz (2017) and apparently no additional reference exists for the problem. The authors suggested algorithms for some polynomial-time-solvable cases and proved that it is NP-complete, at least for α ∈ { 23π , π }. They also conjectured that an NP-completeness proof is likely to exist for α = π2 . Furthermore, they point out that feasible solutions are only guaranteed to exist for it when α ≥ π3 applies. Finally, they proposed constant factor approximation algorithms for different values of α . In application terms, the problem serves as model for the design of wireless networks that rely on directional antennas. This type of antenna irradiates power in directions that span restricted angles and have several advantages over omni-directional antennas, that irradiate power in all directions (Carmin et al., 2011). In

To the best of our knowledge, no α -MSTP formulation was suggested prior to those we introduce here. Indeed, as one shall see, algebraically modelling α -AC is by no means straightforward. One of our formulations, Fx , only uses variables associated with the edges of G, all of them binary (0-1). It is thus a natural space formulation, under the understanding that the space of the edge variables is the natural one for formulating the problem. To enforce α -AC, the formulation ensures that spanning tree edges must satisfy some α -ST necessary and sufficient conditions, to be introduced later on. These conditions arise from a characterization of admissible and non admissible subsets of α -ST edges and translate into a family of exponentially many cover inequalities. On top of the natural space variables, above, the other formulation, Fxy , also uses an additional family of binary 0-1 variables. These variables, as one shall see, are specifically destined to enforcing α -AC. Accordingly, given that the formulation is not restricted to the natural space variables, it is thus an extended formulation for α -MSTP. Among others, we carried out a polyhedral comparison of our two formulations. It indicates that Fxy is at least as strong as Fx , since the projection of its Linear Programming relaxation (LPR) polyhedron onto the space of the edge variables is contained into that for Fx . We also show that α = π is a critical value for modeling α -MSTP, at least as far as Fx is concerned. When α < π , we prove that only O(n3 ) cover inequalities are indeed necessary to obtain an α -MSTP formulation. Additionally, we also establish that the LPR polyhedron for that formulation is precisely the same one would otherwise obtain if all the exponentially many cover inequalities were used. Accordingly, remaining cover inequalities may therefore be removed from the formulation without any effect on LPR bounds. From an experimental standpoint, we implemented and tested Branch-and-Cut (BC) algorithms (Padberg and Rinaldi, 1991) for solving the two formulations above. Additionally, we also designed and tested a greedy style heuristic akin to the minimum weight spanning tree algorithm of Kruskal (1956). For test instances associated with smaller values of α , LPR bounds for Fxy are significantly stronger than those for Fx . This empirical result nicely complements our theoretical one, which states that Fxy is at least as strong as Fx , for any value of α . Indeed, it may be much stronger, precisely for those instances associated with smaller values of α , which turned out to be the hardest to solve in our computational experiments. Such an LPR pattern carries over to the performance of the corresponding BC algorithms, with the one for Fxy having the edge for smaller values of α . However, as α increases, LPR bounds for the two formulations get progressively closer to one another and the algorithm for Fx becomes more and more competitive. In actual fact, for larger values of α , it produces much better computational results than its counterpart. The two formulations are by then nearly equally strong but Fx involves substantially fewer decision variables and constraints, which gives its corresponding BC algorithm the edge. In summary, none of the BC algorithms introduced here dominates the other (in performance terms) over the whole spectrum of possible values for α . The paper is organized as follows. Section 2 introduces the notation we use. Section 3 presents our two integer programming (IP) formulations for α -MSTP and compares their LPR bounds from a theoretical perspective. Section 4 investigates α related

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775

geometrical properties for Fx . Our greedy α -MSTP heuristic is then described in Section 5. It is followed, in Section 6, with the two BC algorithms. Computational experiments are reported in Section 7. Section 8 offers a summary of our conclusions and presents directions for future investigations. Finally, Appendix A suggests a procedure to separate exactly the exponentially many cover inequalities contained in Fx . Namely, by solving a series of integer programs for every i ∈ V. The procedure remains to be implemented and tested. However, it is bound to be CPU time intensive and unlikely to become an alternative to the heuristic strategy we use. It is presented here, however, for completeness and to eventually attracting interest on this still open α -MSTP related computational complexity issue.

2. Notation For every S⊆V, denote by δ (S) (resp. E(S)) the edges of G with one end vertex (resp. the two end vertices) in S. For simplicity, if S contains a single vertex, say vertex i, we use δ (i) instead of δ ({i}). Extending the notation to a spanning tree T = (V, ET ) of G, we denote by δ T (i) := δ (i) ∩ ET its i-incident edges. Additionally, let D = {di : i ∈ V } define the set of closed unit disks with centers on the points (vertices) of V. Furthermore, let Ri = {ij : j ∈ V \ {i}} denote the unit vectors that are collinear and have the same directions as those vectors having i ∈ V as their initial point and j ∈ Vࢨ{i} as their terminal ones. Furthermore, for any i ∈ V, denote by i the equal vector to (1 0 )T , i.e., a vector with the same magnitude and direction as (1 0 )T but initial point in i ∈ V. Fig. 2 illustrates these unit vectors for the line segments (i.e., edges) in Fig. 1. They are accompanied by the boundary of di , where all terminal points are located. Note that Ri vectors relate directly to the edges of δ (i).

2.1. Angles between i-incident edges We now define ࢭ0ij as the anticlockwise angle that i forms with  i j. For every vector ij ∈ Ri in Fig. 2, corresponding angles ࢭ0ij are shown in Fig. 3. Let us also consider the angle obtained by rotating, say ij, anti , j = k. Denote clockwise around i until it becomes collinear with ik that angle by ࢭjik . Likewise, denote by ࢭkij the similarly obtained  change roles in the scheme above. angle when ij and ik j For a simple procedure to compute ࢭjik and ࢭkij , let 0ik = 0ik ,

if ࢭ0ik ≥ ࢭ0ij , and 0ik = 2π + 0ik , if ࢭ0ik < ࢭ0ij otherwise applies. j

j

It then follows that  jik = 0ik − 0i j . Quite clearly, angle ࢭkij may be similarly computed.

3

Table 2 Angles ࢭjik computed for every distinct pair of ordered vectors  in Fig. 2, with a highlight for the largest corresponding ij and ik angle for every individual vector, ij. j

z q u t l

k z

q

u

t

l

Largest ࢭjik

-

π

7π 12 5π 12

13π 12 11π 12

19π 12 17π 12

19π 12 11π 6 19π 12 3π 2 3π 2

11π 6 17π 12 11π 12 5π 12

6

-

19π 12 13π 12 7π 12

-

3π 2

π

π 2

-

3π 2

π π 2

-

 happen to be the same vector, i.e., if  =  When ij and ik 0i j 0ik  are difapplies,  jik = ki j = 0 then results. Conversely, if ij and ik ferent vectors,  jik + ki j = 2π then holds.  in Fig. 2, For every ordered pair of distinct vectors ij and ik Table 2 indicates their corresponding angle ࢭjik . Furthermore, for every individual vector ij in that figure, the largest of these values is also highlighted in the table. 2.2. Circular sectors and the α -MSTP angular constraints We will now extend the procedure above to compute a particularly important circular sector for disk di , i ∈ V. Namely, the smallest angled circular sector that simultaneously encompasses all vectors in a given Zi , Zi ⊆Ri . The resulting procedure, as one shall see, will allow us to model α -AC, as required by α -MSTP. In order to describe it, let s(i)⊆δ (i) be the set of edges of G corresponding to the vectors in Zi . Note that by rotating any given ij ∈ Zi anticlockwise around i, it is straightforward to compute the angles ࢭjil that ij forms with every il ∈ Zi \ {ij}. Out of these, assume that ࢭjim is the one with the largest value and consider the circular sector it implies in disk di . Such a sector clearly encloses all vectors in Zi and should be kept for the moment. Note, however, that it may have a central angle larger or smaller than other similarly computed circular sectors, initiated at Zi vectors other than ij. Accordingly, apply the procedure above to every individual ij ∈ Zi . Then, amongst all circular sectors with central angles ࢭjim , keep the one with the smallest angle. Denote it by ࢮi and its central angle by θ (ࢮi ) and note that this is the circular sector we seek. From the discussion above, identifying ࢮi thus follows from the identification of one circular sector for every ij ∈ Zi , i.e., the one with the largest central angle ࢭjim . Among these, the one with the smallest central angle then defines ࢮi . For the vertices indicated in Fig. 1 and Zi = Ri , note that θ (i ) = 32π then results (see Table 2 for details). 2.3. Locally feasible subsets of i-incident edges For every i ∈ V and every ij ∈ Ri , define

Li j := { j} ∪ {k ∈ V

\ {i, j} : ki j ≤ α}

and note that k ∈ Lij if and only if one reaches or goes past ij by  anticlockwise around i by an angle of α . Therefore, all rotating ik vertices k ∈ Lij are such that k0i j , which is equal to either ࢭ0ij or

Fig. 3. Angles ࢭ0ij defined by every different vector ij in Fig. 2, with horizontal axis support for i depicted for reference.

0i j + 2π , by definition, satisfies k0i j ∈ [0ik , 0ik + α ]. For the unit vectors Ri in Fig. 2, Table 3 shows the different sets Lij that apply to every different α ∈ {π , 23π , 32π }. Feasibility of a spanning tree T = (V, ET ) may thus be checked one vertex i ∈ V at a time, as follows. Define iT as the smallest angled circular sector that simultaneously encompasses all vectors in Zi = {ij ∈ Ri : {i, j} ∈ δT (i )}. If i is a leaf of T, we define θ (iT ) =

4

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775

spanning tree edge e = {i, j} ∈ δ (i ) for each i ∈ V, and (ii) if yi j = 1 holds, the vertices that may communicate directly with i are constrained to those k ∈ Vࢨ{i}: j ∈ Lik . Based on the information above, an optimal Integer Programming (IP) solution to α -MSTP is thus obtained by solving

Table 3 Subsets of vertices Lij for every vector ij ∈ Ri in Fig. 2 and every angle α in {π , 23π , 32π }.

α j z q u t l

π

2π 3

3π 2

{t, l, z} {l, z, q} {l, z, q, u} {q, u, t} {u, t, l}

{l, z} {l, z, q} {z, q, u} {u, t} {t, l}

{u, t, l, z} {t, l, z, q} {t, l, z, q, u} {l, z, q, u, t} {q, u, t, l}



min

 

we xe : (x, y ) ∈ Fxy ∩ B

,

(5)

where Fxy is the polyhedral region defined as

0. Otherwise, if |δ T (i)| ≥ 2 applies, θ (iT ) is then given by





3m

e∈E







xe = n − 1

(6)

e∈E



xe ≤ |S | − 1

S ⊂ V, S = ∅

(7)

i∈V

(8)

yi j ≤ xe

e = {i, j} ∈ E, i < j

(9)

y ji ≤ xe

e = {i, j} ∈ E, i > j

(10)

yik

e = {i, j} ∈ E, i < j

(11)

y jk

e = {i, j} ∈ E, i > j

(12)

xe ≥ 0

e = {i, j} ∈ E

(13)

3. Integer programming formulations

yi j ≥ 0

i ∈ V, ij ∈ Ri ,

(14)

Let G = (V, E ) be a graph corresponding to Euclidean plane points, as previously described. Additionally, denote by T the set of all spanning trees for it. Problem α -MSTP is then generically described as the following combinatorial optimization problem:

and an LPR corresponding to (5) is given by

θ (iT ) = min

max

{i, j}∈δT (i )

{i,k}∈δT (i )\{i, j}

 jik

.

(1)

If θ (iT ) ≤ α holds for every i ∈ V, T is an α -ST and is therefore feasible for the problem. Computing iT could be made simpler, as we describe next. Assume that the edges of T incident to i are δT (i ) = {{i, v1 }, {i, v2 }, . . . , {i, v p }} and that they are sorted according to 0iv1 ≤ 0iv2 ≤ · · · ≤ 0iv p , p = |δT (i )|. Given this sorting, note that iv1 forms a largest possible angle with iv p . Likewise, iv2 does that with iv1 and so on until, finally, iv p−1 is the corresponding pair for iv p . Assuming that v0 ≡ vp , for simplicity, θ (i ) may then be

e∈E ( S )



yi j = 1

ij∈Ri

xe ≤

k∈Li j

T

efficiently computed in O(n) time by direct comparison of just p angles as follows:

θ (iT ) = min

a=1,...,p

min







va iva−1 .

(2)

we

e∈ET

(V, ET ) ∈ T

θ (iT ) ≤ α

(3) i ∈ V.

(4)

Two α -MTSP formulations are introduced here, Fx and Fxy , and they basically differ on the way they characterize α -AC (4), in terms of linear inequalities. To simplify the presentation, Fxy will be described first. 3.1. Formulation Fxy Formulation Fxy involves two sets of variables. One, x = {xe ∈ B : e ∈ E }, for selecting spanning tree edges, with xe = 1 holding if e ∈ E is chosen for the spanning tree and xe = 0 applying otherwise. The other set, y = {yi j ∈ B : i, j ∈ V, ij ∈ Ri }, is used to align the directional antennas, so as to allow applicable pairs of vertices to communicate directly. Accordingly, yi j = 1 implies that the antenna placed at i has one of its extreme transmission rays collinear with ij. The other extreme ray is then positioned, anticlockwise, α degrees away from ij. In that case, one of the two conditions for feasible transmission between i and j is naturally enforced, i.e., i can see j. However, antenna positioning at j should also allow j to see i. Note as well that one should also simultaneously guarantee feasible transmission from i to every other vertex sharing a spanning tree edge with i, and back from it to i. As we shall see, Fxy provides such a guarantee by imposing, without loss of generality, the two conditions that follow: (i) yi j = xe = 1 holds for a single



xe ≤



k∈L ji



min





we xe : (x, y ) ∈ Fxy ,

e∈E

with value w(F xy ). On their own, constraints (6),(7) and (13) define an integral polytope whose vertices are the incidence vectors for spanning trees of G Edmonds (1971). Constraint (6), in particular, enforces that the required number of (spanning tree) edges is selected. In turn, inequalities (7), alongside (13), guarantee that the selected edges form a cycle free subgraph. Note as well that inequalities

xe ≤ 1, e = {i, j} ∈ E,

(15)

are naturally implied by (7), written for sets S = {i, j}. Remaining constraints in (6)-(14) are in charge of enforcing α AC, (4). In particular, constraints (8) impose that exactly one of the vectors of Ri is chosen for establishing the positioning of the i-located antenna. Together with (14), they also ensure that y ≤ 1 necessarily holds. As for (9) and (10), they respectively ensure, provided e = {i, j} is a spanning tree edge, that ij is a candidate for determining the positioning of the i-located antenna and that ji plays a similar role for j. These inequalities are complemented with (11) and (12), which guarantee that an edge e = {i, j} is only taken into a spanning tree if the antennas positioned at both i and j allow direct communication between the two. Finally, all these observations put together ensure that Fxy is a valid formulation for α -MSTP. 3.2. Formulation Fx Differently from formulation Fxy , Fx only uses the previously defined variables x. Furthermore, as we shall see, Fx redefines the generic α -AC, (4), as a family of valid linear inequalities, solely

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775

based on variables x. Namely, a family of inequalities that prevents the use of α -ST illegal subsets of edges. Proposition 1, that follows, characterizes these subsets. Proposition 1. Consider a subset of edges {{i, v1 }, . . . , {i, v p }}, p = |s(i )|, indexed so that

s(i)⊆δ (i),

0iv1 ≤ · · · ≤ 0iv p

vk ivk−1 > α



min



we xe : x ∈ Fx ∩ R

m

.

(22)

e∈E



(17)



Complementing the description above, note that an LPR for (22) is given by

min





we xe : x ∈ Fx ∩ R

m

,

(23)

e∈E

k=1

with value w(F x ).

holds,

|δT (i ) ∩ s(i )| ≤ |s(i )| − 1

(18)

is then valid for every α -ST. Proof. The proof is by contradiction and assumes the existence of an α -angled i-centered circular sector that encompasses all unit  : {i, v } ∈ s(i )}. Note, in this case, that one may vectors in Zi = {iv k k restrict attention to circular sectors with their initial (extreme) rays  ∈ Z is collinear with a vector in Zi . Accordingly, assume that iv i k such a vector. Namely, one that when rotated anticlockwise around i by an angle of α coincides or goes past every other vector in  ∈ Z from qualifying Zi . Note, however, that (17) prevents any iv i k for the role. A contradiction is thus reached and the result then follows.  Let us now restate condition (17) in terms of

v1 ∈ Liv p ∧ v2 ∈ Liv1 ∧ · · · ∧ v p ∈ Liv p−1 ,

(19)

Proposition 2. Let s(i ) = {{i, v1 }, {i, v2 }, . . . , {i, v p }} be a subset of δ (i) which satisfies (16) and is such that

∨v1 ∈ Liv p ∨ v2 ∈ Liv1 ∨ · · · ∨ v p ∈ Liv p−1

(20)

holds. Then, there must exist a circular sector centered at i and with a central angle at most α which encloses all vectors in Zi = {ij : {i, j} ∈ s ( i )}. Proof. Assume, without loss of generality, that the first clause in (19) is not satisfied, i.e., that v1 ∈ Liv p holds. Then, a circular sector enclosing all vectors in Z is obtained by rotating iv1 anticlockwise i

around i by angle of α . Accordingly, the result thus follows.



Let us call α -ST non-admissible, a nonempty subset of edges s(i)⊆δ (i) satisfying conditions (16) and (17) (or equivalently, condition (19)). Conversely, we will call α -ST locally admissible, a set satisfying (16) and (20). The modeling and algorithmic importance of Propositions 1 and 2 should not be undervalued for the easiness in stating and proving them. Just by itself, Proposition 2 provides a quite effective and easily-verifiable procedure for building an α -MSTP constructive heuristic. Indeed, it is the core to our greedy heuristic to the problem, to be described in Section 5. Furthermore, as we shall see next, by combining the two propositions one is capable of modeling α -AC, (4), exclusively in terms of the canonical variables x. Lemma 1. Let Fx be the polyhedral region defined by the intersection of constraints (6),(7), and (13) with the exponentially many cover inequalities

xe ≤ |s ( i )| − 1, s ( i ) ⊆ δ ( i ),

e∈s ( i )

|s(i )| ≥ 2, s(i ) satisfies (16 ) − (17 ).

3.3. A polyhedral comparison of Fx and Fxy For the remainder of this section we will be comparing Fx and Fxy , in terms of the strength of their corresponding LPR bounds. To that aim, we use Proposition 3, which is presented next and that follows from Proposition 1. Proposition 3. Assume one is given a proper subset of edges s(i) ⊂ δ (i), i ∈ V. Additionally, assume that s(i ) = {{i, v1 }, . . . , {i, v p }}

satisfies (16) and (17) and let Lsi j(i ) = Li j ∩ {k ∈ V \ {i} : {i, k} ∈ s(i )} be the subset of vertices of Lij contained in {v1 , . . . , v p }. Then

Lsiv(i ) ∩ · · · ∩ Lsiv(pi ) = ∅

(24)

1

and

so as to characterize feasible sets of α -ST edges, as defined by the following proposition.





Proof. The proof is given by Propositions 1 and 2.

(16)



Then, α −MSTP can be modeled by the following integer program

s (i ) =

applies and v0 ≡ vp is used, for convenience. If p 

5

(21)

Liv1 ∩ · · · ∩ Liv p = ∅

(25)

necessarily hold for s(i). Proof. Let us assume, once again, that v0 ≡ vp holds, for convenience. Validity of (24) is established by contradiction, after

p assuming the existence of vu ∈ k=1 Lsiv(i ) . Accordingly, take vu ∈ k

{v1 , . . . , v p } and note that vu ∈ Lsiv(i) should then hold. However, u−1 by hypothesis, vu ivu−1 > α is assumed and vu ∈ Livu−1 must then result. Additionally, given that Lsiv(i ) ⊆ Livk for every k ∈ {1, . . . , p}, k

one cannot have vu ∈ Lsiv(i ) . A contradiction is thus reached and u−1 (24) then follows. Validity of (25) follows from the validity of (24) and, in order to prove it, initially set s (i ) = s(i ). Then pick any edge {i, j} ∈ δ (i)ࢨs (i) and, without loss of generality, let {i, vu }, {i, vu−1 } ∈ s(i ) be such that 0i j ∈ [0ivu−1 , 0ivu ]. Now note, by hypothesis, that vu ivu−1 > α (or equivalently, vu ∈ Livu−1 ). Furthermore, given that 0i j ∈ [0ivu−1 , 0ivu ], note that vu ∈ Li j also applies. Accords (i )∪{i, j}

ingly, vu ∈ Liv s (i )∪{i, j}

Liv

1

u−1

s (i )∪{i, j}

∩ Liv

2

s (i )∪{i, j}

, since Liv u−1

s (i )∪{i, j}

∩ · · · ∩ Liv p

⊆ Livu−1 , and therefore

= ∅.

Now, update s (i) ← s (i) ∪ {i, j} and recursively apply the same reasoning until s (i ) = δ (i ) is attained, and the proof is thus obtained.



One of the main reasons for introducing Proposition 3 is to uncover the relation





e={i, j}∈s(i ) k∈Li j

yik ≤ (|s(i )| − 1 )



yik

(26)

k∈∪e={i, j}∈s(i ) Li j

for s(i) ⊂ δ (i), i ∈ V. As one shall see, (26) is used in our proof that formulation Fxy is at least as strong as formulation Fx , as stated in the theorem that follows. Theorem 1. (x, y ) ∈ Fxy ⇒ x ∈ Fx and, therefore, w(F xy ) ≥ w(Fx ).

6

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775

Proof. Pick any (x, y ) ∈ Fxy , i ∈ V and s(i)⊆δ (i). One then has:  e∈s ( i )

xe ≤





e={i, j}∈s (i ) k∈Li j

≤ ( |s ( i )| − 1 )

summing up (11) for{i, j} ∈ s(i )

yik 

 k∈ e={i, j}∈s (i ) Li j

yik

due to (25), every k appears at most

|s(i )| − 1times in



Li j ,

e={i, j}∈s (i )

and y ≤ 1 = |s ( i )| − 1

due to (8)

Inequalities (21) are thus satisfied by x in (x, y). Accordingly, it then follows that x ∈ Fx and that w(Fxy ) ≥ w(Fx ) therefore holds.  4. Geometrical properties of F x as a function of α We now concentrate on a particular relaxation of Fx , Fx23 , that restricts the use of cover inequalities (21) to subsets s(i)⊆δ (i) with |s(i)| ∈ {2, 3}. Fx23 is of special interest since, as we shall see, for α < π it defines not only an α -MSTP formulation but also implies a polytope that is exactly the same as Fx . The result holds despite the fact that Fx contains O(2n ) cover inequalities while Fx23 contains only O(n3 ) of them. We motivate the investigation, with the example depicted in Fig. 4. It shows that, in general, Fx23 does not define a formulation for α −MSTP. In particular, for α ∈ [π , 32π ), the set s(i ) = δ (i ) of four edges used in the example is α -ST non-admissible, as their corresponding unit vectors in the picture confirm. Note as well that all proper subsets of that s(i) turn out to be α -ST locally admissible. Accordingly, quite clearly Fx23 is not an α -MSTP formulation for the example described. However, the situation changes when α < π is considered and Fx23 then becomes a formulation for the problem. Theorem 2. For a graph G = (V, E ) as previously defined and α < π , assume that a set of edges s(i ) = {{i, v1 }, . . . , {i, v p }}, p ≥ 4, is given. Assume as well that s(i) satisfies condition (19) and is therefore α ST non-admissible for G. Then, there must exist a subset s (i) ⊂ s(i), with |s (i)| ∈ {2, 3}, that also satisfies (19) and, therefore, is α -ST nonadmissible for G as well. Consequently, Fx23 is an α -MSTP formulation, if α < π .

Fig. 5. Figures for the proof of Theorem 2. The figure was constructed for α = 56π .

and vq ivt > α . If both hold together, define s (i ) = {{i, vt }, {i, vq }}. Since vq ∈ Livt and vt ∈ Livq , edges in s (i) are α -ST non-admissible and xivt + xivq ≤ 1 is thus a valid inequality for α −MSTP. Therefore, from now on, we assume that all subsets of s(i) with exactly two edges are α -ST locally admissible. Now assume, yet again, that v0 ≡ vp holds, for convenience.   Analogously to (2), define θ (is(i ) ) = min vz ivz−1 : z = 1, . . . , p . By hypothesis, s(i) is non-admissible, so that θ (is(i ) ) > α must then hold. Without loss of generality, consider two consecutive vectors in Zi , say iv1 , iv2 , and their corresponding edges in s(i), i.e., {i, v1 }, {i, v2 }. Since the edges in s(i) satisfy condition (19), v2 ∈ Liv1 and consequently v1 ∈ Liv2 (or, equivalently, v1 iv2 ≤ α ) must necessarily hold since, otherwise, {{i, v1 }, {i, v2 }} would be α -ST nonadmissible.  ∈ Z ,  Now define ivt = arg min{0ivk : iv i v1 ivk > α}. Such an ivt k must necessarily exist since, otherwise, θ (is(i ) ) would not exceed α . Furthermore, vt iv1 ≤ α (or, equivalently, vt ∈ Liv1 ) must then apply since, otherwise, {{i, v1 }, {i, vt }} would be α -ST non-admissible. Accordingly, one thus have

v2 ∈ Liv1 and v1 ∈ Livt .

(27)

Then, consider the following two possibilities, respectively illustrated in Fig. 5(a) and (b). •

 : {i, v } ∈ s(i )} and recall that the edges in s(i) Proof. Let Zi = {iv k k are defined so that 0iv1 ≤ 0iv2 ≤ · · · 0iv p holds. Then pick any two distinct s(i) edges, for instance, {i, vt }, {i, vq }, q, t ∈ {1, . . . , p}, with q < t. Since α < π , we may simultaneously have vt ivq > α •

vt iv2 > α (or, equivalently, vt ∈ Liv2 ). Combined with (27), we then have 0iv1 ≤ 0iv2 ≤ 0ivt , v2 ∈ Liv1 , vt ∈ Liv2 , v1 ∈ Livt . Set s (i ) = {{i, v1 }, {i, v2 }, {i, vt }} and note that s (i) satisfies condition (19). Therefore, s (i) is α -ST nonadmissible and cover inequality xiv1 + xiv2 + xivt ≤ 2 must be satisfied. vt iv2 ≤ α (or, equivalently, vt ∈ Liv2 ). Then, by the definition of s(i) and ivt , there must exist ivt−1 such that v1 ivt−1 ≤ α and vt ivt−1 > α . Since α < π , we have that vt−1 iv1 > α , i.e., vt−1 ∈ Liv1 . Now, define s (i ) = {{i, v1 }, {i, vt−1 }, {i, vt }} and note that: 0iv1 ≤ 0ivt−1 ≤ 0ivt , vt ∈ Livt−1 , v1 ∈ Livt and vt−1 ∈ Liv1 . Therefore, s (i) is α -ST nonadmissible and the proof is complete. 

Fig. 4. Example showing that, in general, Fx23 does not define a formulation for α MSTP. Unit vectors depicted in the figure are such that ziq = qiu = uit = tiz = π . Taking any α ∈ [π , 3π ) and setting s (i ) = {{i, z}, {i, q}, {i, u}, {i, t }}, s(i) is α -ST 2 2 non-admissible and thus x(s(i)) ≤ 3 defines a valid inequality for α −MSTP. However, it can be easily checked that all subsets s (i) ⊂ s(i), s (i ) = ∅ are α -ST locally admissible.

From a geometrical standpoint, the latter results can be interpreted as follows. No matter how close to π parameter α gets from below, Fx23 always define a formulation for α −MSTP. At the critical value α = π , and for values of α larger than π , Fx23 may well contain integer vectors violating α -AC, (4). Therefore, in order to have an x-space α -MSTP formulation for α ≥ π , additional cover inequalities (21) are thus required. Namely those corresponding to s(i)⊆δ (i), |s(i)| ≥ 4. As it turns out, polyhedral region Fx has another quite interesting property. Namely, Fx = Fx23 applies when α < π holds. In

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775

other words, inequalities (21) for s(i)⊆δ (i), |s(i)| ≥ 4, then become redundant in the polyhedral description of Fx . Accordingly, polynomial time separation of cover inequalities (21) thus results when α < π , given that only O(n3 ) of these inequalities need be checked. Theorem 3, that follows, formalizes this result. Theorem 3. If α < π then Fx =

Fx23 .

Proof. Since Fx23 is a relaxation for Fx , Fx ⊆ Fx23 obviously holds. Accordingly, all one needs to shown is that Fx23 ⊆ Fx also applies. Pick any x ∈ Fx23 and any s (i)⊆δ (i) that simultaneously satisfies condition (19) and |s (i)| ∈ {2, 3}. Clearly,

 applies. Now pick any s(i)⊆δ (i) with e∈s (i ) xe ≤ |s (i )| − 1 then

 s (i) ⊂ s(i). Since x ≤ 1, e∈ (s (i )\s (i )) xe ≤ |s (i ) \ s (i )| thus results.

 Hence, e ∈ s ( i ) xe = e∈s (i ) xe + e∈ (s (i )\s (i )) xe ≤ |s (i )| − 1 + |s (i ) \  s (i )| = |s(i )| − 1, x ∈ Fx and the results follows.  5. A greedy α-MSTP heuristic The heuristic introduced here is based on Proposition 2, that provides an efficient procedure for identifying whether or not any subset of edges s(i)⊆δ (i), i ∈ V, is α −ST locally admissible for G. As for any Kruskal-like minimum spanning tree heuristic, ours select edges of G in nondecreasing weight order, provided they do no form a cycle with previously selected edges. Additionally, for our specific case, one must also check whether or not edge inclusion maintains α -AC feasibility. Assume E = {e1 , . . . , em }, we1 ≤ we2 ≤ · · · ≤ wem , is given. Accordingly, edge e1 is necessarily selected at the first iteration of the procedure and T = (V, ET ), ET = {e1 } thus results. At iteration k ≥ 2, an edge e = {i, j} ∈ E only qualifies as a candidate for eventually entering T if and only if i and j do not belong to a same component of T. Provided that is the case, e = {i, j} is only finally accepted in T if δ T (i) ∪ {i, j} and δ T (j) ∪ {i, j} are both α -ST locally admissible for G. The procedure is carried out until |ET | = n − 1 holds, and an α -ST is thus obtained for G, or, else, no edge of E remains to be checked and the algorithm fails to find such a tree. The procedure, that we shall call Greedy α -ST (G-α -ST), is formalized as follows: (Initialization) ET ← ∅, k ← 1 While |ET | = n − 1 ∧ k ≤ m do - Pick the next candidate edge ek = {i, j}. - If (V, ET ∪ {i, j}) is cycle free then: ∗ Set s(i ) = δT (i ) ∪ {i, j}. Order edges in s(i) according to (16). ∗ Set s( j ) = δT ( j ) ∪ {i, j}. Order edges in s(j) according to (16). ∗ If condition (20) holds for s(i) and s(j), update ET ← ET ∪ {i, j}. - k←k+1 G-α -ST clearly runs in polynomial time. Edges of G can be sorted in O(mlog (n)) time and checking whether or not (V, ET ∪ {i, j}) is cycle free may be done in O(1) time, by using union-find data structures for disjoint sets (see Cap. 21 in Cormen et al., 2009). On the other hand, the edges of δ T (i) can be kept ordered in O(log (n)) time per edge inclusion. Accordingly, checking if s(i) satisfies (20) may be carried out in O(n) time and the overall running time for G-α -ST is therefore no worse than O(mn). 6. Branch-and-cut algorithms In this section, we introduce two BC algorithms for α −MSTP. The first one, BCFXY is based on formulation Fxy . The second, BCFX, relies on Fx . Implementation details for BCFXY are presented first.

7

6.1. A BC algorithm for Fxy

BCFXY starts off by solving the linear program   min



we xe : (x, y ) ∈ F xy ,

(28)

e∈E

where F xy is a relaxation of Fxy given by the intersection of (6), (8)–(14) and (15). Let (x, y ) ∈ [0, 1]3m be an optimal solution for (28). BCFXY then atempts to identify Subtour Elimination Constraints (SECs) (7) that violate x, i.e., looks for subsets of vertices S ⊂ V for which

e∈E (S ) xe > |S| − 1 holds. Separation of inequalities (7) is initially carried out heuristically, complemented with an exact approach, if necessary (see Section 6.1.1, for details). Violated SECs are then appended to F xy , thus resulting in a strengthened version of that program, which is then re-optimized. The procedure is repeated until x violates no SECs and, therefore, LP relaxation bound w(Fxy ) is thus obtained. At this stage, if (x, y ) ∈ B3m holds, it defines an optimal solution for (5). Conversely, if it is fractional, BCFXY then branches on fractional variables. The cutting plane algorithm just described is applied at every every node of the BC tree. BCFXY is implemented under the XPRESS mixed integer optimization package, release 8.4 (2017). XPRESS is thus responsible for solving LPRs, (28), and managing the BC tree. It uses default options to choose a variable to branch and implements a best-first search strategy. Additional features offered by XPRESS, such as automatic cut generation and primal heuristics, are kept switched off. Likewise, multi-threading is not used as well.

6.1.1. Separation of subtour elimination constraints The cutting plane engine in BCFXY implements both exact and heuristic separation of SECs. At every BC node, exact separation is only carried out if heuristic separation fails. For the descriptions that follow, assume that x ∈ [0, 1]m denotes the x component of the point (x, y ), to be separated, and define E := {{i, j} ∈ E : xi j > 0}. Accordingly, G = (V, E ) defines the support

graph for x. A SEC is assumed to be violated if e∈E (S ) xe − |S| + 1 > 10−2 holds. All violated SECs, identified either by the heuristic or the exact algorithm, are appended to the relaxation. Our SEC separation heuristic generates a Maximum Weight Spanning Tree of G, under edge costs {ce = xe : e ∈ E }. It uses Kruskal’s algorithm (Kruskal, 1956) for that purpose and whenever an edge is selected for entering the solution, the connected component that thus results in that iteration of the algorithm is then checked for SEC violation at x. For the exact separation of SECs, we rely on the algorithm in Padberg and Wolsey (1983). It solves (n − 2 ) max-flow/min-cut problems. These are conveniently defined over directed networks originating from G. Every problem identifying a subset of vertices S ⊂ V ensured to contain some specific, pre-defined, vertices of G. Under the conditions imposed, these subsets S are the ones most prone for SEC violation at x. Each of them is then checked for that condition. If no violation occurs, x then belongs to the convex hull of incidence vectors for spanning trees of G and BCFXY therefore resorts to branching. 6.1.2. α -MSTP BCFXY upper bounds At the very start, BCFXY is initialized with the α -MSTP upper bound given by heuristic G − α -ST (see Section 5). Aiming at improving that bound, the heuristic is also applied right before branching from every applicable BC node. In that case, it is then run under complementary weights, {(1 − xe )we : e ∈ E }, generated for the LP relaxation at hand.

8

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775 Table 4 Computational results - α = π3 . Instance

w( · )

n

BLB

BUB

t(s)

opt ? y/n

eil51 eil51 eil51 berlin52 berlin52 berlin52 st70 st70 st70 pr76 pr76 pr76 eil76 eil76 eil76 rat99 rat99 rat99 kroA100 kroA100 kroA100 kroB100 kroB100 kroB100 kroC100 kroC100 kroC100 kroD100 kroD100 kroD100 kroE100 kroE100 kroE100 eil101 eil101 eil101 lin105 lin105 lin105

15 30 51 15 30 52 20 40 70 25 50 76 25 50 76 35 70 99 35 70 100 35 70 100 35 70 100 35 70 100 35 70 100 35 70 101 35 70 105

w(Fxy )

w ( Fx )

BCFXY

BCFX

BCFXY

BCFX

BCFXY

BCFX

BCFXY

BCFX

243.1 363.0 517.6 4292.5 6167.3 8037.7 371.9 563.9 803.0 38892.3 74164.6 120116.4 357.5 475.1 650.1 556.5 1104.3 1563.2 16251.5 22465.2 25666.0 15583.2 23387.9 26911.3 16537.5 21324.2 25415.4 16498.0 21811.6 25426.3 17114.4 23127.6 26538.1 447.9 690.1 750.2 5261.5 11825.5 17007.4

203.5 316.2 442.6 4066.6 5481.9 7360.4 345.2 503.4 699.1 35231.9 68523.7 111340.1 306.1 412.2 562.1 493.7 964.6 1362.6 14706.6 20199.8 23195.1 13853.3 20750.5 23762.0 14568.3 19368.7 22857.5 14707.2 19592.9 23049.1 15452.0 20776.5 23802.4 382.8 589.5 652.3 4987.2 11026.1 15723.8

258.1 409.9 558.7 5262.7 7090.3 9105.0 420.5 628.6 857.8 41809.0 83009.5 132689.1 395.0 511.8 686.2 633.2 1177.8 1639.5 19305.4 24450.0 26845.5 17263.2 25368.7 28148.3 18639.7 23182.0 26708.4 19250.4 23304.3 26760.9 19786.2 24777.7 27575.1 499.1 729.7 777.2 5751.0 12568.5 17708.4

258.1 393.0 507.0 5262.7 6845.3 7985.2 420.5 604.6 782.8 41809.0 78401.5 120327.6 395.0 472.4 609.9 605.0 1071.1 1466.1 17790.5 22278.9 24907.6 17080.0 23078.8 25568.1 17803.4 21493.5 24525.9 18081.5 21865.0 24760.0 18841.6 22728.3 25437.1 466.5 645.0 686.1 5751.0 12058.8 16598.5

258.1 409.9 587.2 5262.7 7090.3 9896.8 420.5 628.6 909.9 41809.0 83009.5 139876.8 395.0 551.3 756.1 633.2 1268.8 1845.5 19707.6 26615.5 31309.7 17263.2 27334.7 30250.1 18639.7 25111.2 31366.9 19250.4 25309.5 30751.6 20012.2 26553.4 30774.2 499.6 791.7 931.1 5751.0 13300.2 19388.1

258.1 413.4 605.7 5262.7 7142.9 10075.5 420.5 635.2 927.1 41809.0 84784.7 142218.9 395.0 560.6 786.3 633.2 1278.0 1846.5 19866.1 27211.7 31775.2 17263.2 28529.1 32465.0 18657.5 25312.0 31877.3 19283.4 25620.1 31132.3 20117.1 27161.0 31873.5 503.9 824.6 908.2 5751.0 13253.0 19795.7

0.2 200.6 0.5 90.9 1.5 305.4 0.8 563.3 14.2 250.3 30.5 125.0 1711.2 88.3 -

0.5 0.8 3.0 4.9 2242.5 2314.6 -

y y n y y n y y n y y n y n n y n n n n n y n n y n n y n n n n n n n n y n n

y n n y n n y n n y n n y n n n n n n n n n n n n n n n n n n n n n n n y n n

6.2. A BC algorithm for Fx

in an updated reinforced relaxation F x to Fx , and the algorithm then iterates.

BCFX and BCFXY share several common features. For instance, as before, XPRESS is also in charge of solving LPRs and managing the search tree. Likewise, a best-first search strategy is also implemented by the solver. Furthermore, BCFX separates SECs precisely as BCFXY does and calls to heuristic G − α -ST are carried out exactly as previously described in Section 6.1.2. Finally, BCFX also branches on variables, according to XPRESS default settings. The main difference between BCFX and BCFXY is that, in addition to SECs, BCFX also separates cover inequalities (21). In order to describe how that separation is carried out, consider the following relaxation of (23):

 min



 we xe : x ∈ F x ,

(29)

e∈E

where polyhedral region F x ⊃ Fx denotes the intersection of constraints (6), (13), and (15). Note that neither SECs (7) nor cover inequalities (21) are initially included in F x . The cutting plane algorithm we use for reinforcing F x always separates SEC and cover inequalities at every solution round to (29). Violated inequalities found are appended to F x , thus resulting

6.2.1. Separation of cover inequalities The separation of cover inequalities, (21), naturally decomposes into n independent problems, one for every i ∈ V. However, apart from our previously described procedure for the α < π case, no other polynomial time algorithm appears to exist for solving these problems. Indeed, for α ≥ π , determining the complexity class for this separation problem apparently remains an open question. Denote by x ∈ [0, 1]m an optimal solution to (29) and define δ (i ) := {e = {i, j} ∈ δ (i ) : xe > 0}. In Appendix A we suggest an exact separation scheme for inequalities (21), for any value of α . It consists in solving a sequence of |δ (i )| − 1 integer programs for every i ∈ V. Differently from Appendix A, BCFX implements another separation strategy for cover inequalities. Namely, one which distinguishes the cases α < π and α ≥ π , in accordance with Theorem 2. Accordingly, when α < π applies, BCFX only separates (21) for sets s(i ) ⊆ δ (i ), |s(i)| ∈ {2, 3}, since Fx = Fx23 then holds. Separation, in this case, carried out in polynomial time by simple enumeration of these sets. For α ≥ π , however, we partition cover inequalities (21) into two subsets, with a distinct separation strategy for each of them.

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775

9

Table 5 Computational results - α = 23π . Instance

n

w( · )

BLB

BUB

t(s)

opt ? y/n

eil51 eil51 eil51 berlin52 berlin52 berlin52 st70 st70 st70 pr76 pr76 pr76 eil76 eil76 eil76 rat99 rat99 rat99 kroA100 kroA100 kroA100 kroB100 kroB100 kroB100 kroC100 kroC100 kroC100 kroD100 kroD100 kroD100 kroE100 kroE100 kroE100 eil101 eil101 eil101 lin105 lin105 lin105

15 30 51 15 30 52 20 40 70 25 50 76 25 50 76 35 70 99 35 70 100 35 70 100 35 70 100 35 70 100 35 70 100 35 70 101 35 70 105

w(Fxy )

w ( Fx )

BCFXY

BCFX

BCFXY

BCFX

BCFXY

BCFX

BCFXY

BCFX

201.2 307.5 425.4 3765.0 5145.5 7032.0 334.6 471.7 652.7 33415.1 66015.4 106023.3 297.3 406.1 546.3 485.0 937.7 1323.0 13805.2 19180.0 22108.6 13474.6 19671.2 22346.3 14180.1 18288.8 21806.1 13857.4 18859.7 21753.6 14672.0 19490.6 22524.2 372.5 580.8 639.1 4748.7 10443.0 14777.5

200.9 307.3 424.6 3765.0 5138.7 7031.8 334.1 471.7 652.5 33415.1 66015.2 106023.3 297.1 405.7 545.6 485.0 937.0 1322.0 13791.3 19154.6 22060.1 13474.6 19667.1 22329.8 14180.1 18287.5 21803.9 13857.4 18829.2 21753.6 14672.0 19490.6 22524.2 372.3 579.4 638.8 4741.9 10443.0 14777.5

209.3 315.3 436.5 4307.2 5347.2 7488.7 340.3 505.1 687.6 35435.9 71920.1 110485.9 307.9 423.6 561.3 504.1 970.2 1362.5 15214.5 19985.2 23024.9 14192.2 20710.3 23231.0 14709.1 19383.9 23558.5 14475.2 19748.4 22444.3 14937.0 20180.5 23317.3 384.3 601.1 654.3 4888.5 10834.0 15219.4

209.3 315.3 436.5 4307.2 5347.2 7470.2 340.3 505.1 687.6 35435.9 71920.1 112197.6 307.9 423.6 565.1 504.1 970.2 1367.1 15214.5 20352.8 23117.0 14192.2 20817.2 23266.6 14709.1 19383.9 22793.4 14475.2 20003.6 22617.4 14937.0 20180.5 23630.3 384.3 602.4 656.1 4888.5 10834.0 15243.1

209.3 315.3 436.5 4307.2 5347.2 7488.7 340.3 505.1 687.6 35435.9 71920.1 115777.8 307.9 423.6 567.8 504.1 970.2 1367.3 15214.5 20374.3 23976.5 14192.2 20844.5 23596.1 14709.1 19383.9 23558.5 14475.2 20104.8 22668.8 14937.0 20180.5 23914.4 384.3 605.3 670.5 4888.5 10834.0 15392.8

209.3 315.3 436.5 4307.2 5347.2 7491.3 340.3 505.1 687.6 35435.9 71920.1 115765.1 307.9 423.6 565.1 504.1 970.2 1367.1 15214.5 20352.8 23981.0 14192.2 20817.2 23512.9 14709.1 19383.9 23710.6 14475.2 20003.6 22617.4 14937.0 20180.5 23834.6 384.3 604.9 670.4 4888.5 10834.0 15423.0

0.1 0.6 12.0 0.2 1.3 1448.4 0.1 9.1 365.5 0.8 755.4 1.2 876.7 3.5 598.4 169.1 3.7 1.5 2889.2 7013.3 3.4 0.6 252.3 2.6 2.0 250.9 -

0.0 0.2 1.6 0.1 0.2 0.0 3.9 146.1 0.1 1148.1 0.2 43.7 724.9 0.4 28.5 4609.6 16.3 836.0 0.5 369.6 0.3 221.7 0.5 687.4 3745.7 0.1 12.4 1.2 0.5 97.5 -

y y y y y y y y y y y n y y n y y n y n n y n n y y y y n n y y n y n n y y n

y y y y y n y y y y y n y y y y y y y y n y y n y y n y y y y y n y n n y y n

For α ≥ π and sets s(i ) ⊆ δ (i ), i ∈ V, |s(i)| ∈ {2, 3}, we implement the same polynomial time separation scheme indicated above. However, when |s(i)| ≥ 4 applies and x is fractional, we simply skip separation for the corresponding cover inequalities. Indeed, they are only separated when x is integral valued. Accordingly, a lazyconstraints separation strategy is implemented in that case. Such a separation, one should note, is carried out in polynomial time, by simple inspection. The cover inequalities separation scheme implemented by BCFX prevent us from incurring the cost of solving the integer programs in Appendix A. Furthermore, as our empirical results indicate, it does not bring a significant loss in lower bound quality. Finally, if x is integer, Appendix A shows that every violated inequality is violated by exactly the same amount. In particular, assume violation for a cover inequality defined by a set s (i ) ⊂ δ (i ). Then, every cover inequality defined by a set s(i), s (i ) ⊂ s(i ) ⊆ δ (i ), is violated as well. Accordingly, if x is integer, in order to attain the same separation benefits offered by the procedure in Appendix A, one just needs to check whether or not s(i ) = δ (i ) is α -ST locally admissible. 7. Computational experiments A description follows of the computational experiments we carried out for BCFX and BCFXY. All algorithms involved were imple-

mented in C and compiled with gcc, with optimization flags -O3 turned on, under the Linux operating system (release LTS 14.04). A computer equipped with an Intel XEON E5645 processor, running at 2.4GHz and having 32Gb of RAM memory (12Mb of cache memory), was used for the experiments.

7.1. Test instances No test instance specifically designed for α -MSTP appears to exist in the literature. Accordingly, we then opted for adapting to it graphs corresponding to Euclidean Traveling Salesman Problem (ETSP) instances, found in the TSPLIB (Reinelt, 1991). More specifically, we used their underlying set of Euclidean plane points for that purpose. Accordingly, out of every set of points indicated above, three distinct graphs G = (V, E ) were generated. One, corresponding exactly to the TSPLIB graph. The other two, involving a lesser number of vertices, i.e., points. Take, for example, TSPLIB instance eil51, that involves 51 vertices. Graphs with n ∈ {15, 30, 51} vertices were then generated out of the original set of 51 Euclidean plane points. The largest of them, the TSPLIB graph itself. The other two, involving respectively the first 15 and the first 30 points. Edge sets E were always complete, no matter the value of n. Furthermore, for every applicable pair of vertices, their corresponding Euclidean distances, {we : e = {i, j} ∈ E }, were taken as edge weights. Full

10

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775 Table 6 Computational results - α = π . Instance

n

w( · )

BLB

BUB

t(s)

opt ? y/n

eil51 eil51 eil51 berlin52 berlin52 berlin52 st70 st70 st70 pr76 pr76 pr76 eil76 eil76 eil76 rat99 rat99 rat99 kroA100 kroA100 kroA100 kroB100 kroB100 kroB100 kroC100 kroC100 kroC100 kroD100 kroD100 kroD100 kroE100 kroE100 kroE100 eil101 eil101 eil101 lin105 lin105 lin105

15 30 51 15 30 52 20 40 70 25 50 76 25 50 76 35 70 99 35 70 100 35 70 100 35 70 100 35 70 100 35 70 100 35 70 101 35 70 105

w(Fxy )

w ( Fx )

BCFXY

BCFX

BCFXY

BCFX

BCFXY

BCFX

BCFXY

BCFX

178.2 271.2 387.6 3264.2 4699.7 6305.7 289.3 412.5 582.7 28001.7 54386.6 88940.5 265.2 355.8 489.0 403.7 808.8 1153.7 11923.2 16689.6 19220.8 11465.3 17271.9 19855.5 12354.4 15893.0 18896.1 12424.0 16247.9 19057.2 12841.0 17157.7 19933.4 329.9 516.8 578.5 4148.3 9093.0 13183.1

178.2 271.2 387.6 3264.2 4699.7 6305.7 288.6 412.5 581.1 28001.7 54386.6 88940.5 265.2 355.7 488.9 403.7 808.8 1153.7 11923.2 16689.6 19202.4 11465.3 17271.9 19854.4 12354.4 15893.0 18896.1 12417.4 16247.9 19057.2 12841.0 17157.7 19933.4 329.9 516.8 578.5 4146.7 9086.5 13183.1

178.2 271.3 388.4 3264.2 4739.1 6402.7 289.9 415.4 582.9 28179.9 54870.8 89932.5 265.2 357.2 489.2 403.7 808.8 1154.3 12047.7 16729.0 19261.9 11516.4 17294.9 19934.9 12499.3 15986.6 18966.1 12503.6 16375.9 19160.5 13062.4 17235.3 20074.5 329.9 516.8 579.5 4178.8 9118.6 13188.8

178.2 271.3 388.4 3264.2 4739.1 6402.7 289.9 415.4 582.9 28179.9 54870.8 89932.5 265.2 357.2 489.2 403.7 808.8 1154.3 12047.7 16729.0 19261.9 11516.4 17294.9 19934.9 12499.3 15986.6 18966.1 12503.6 16375.9 19160.5 13062.4 17235.3 20074.5 329.9 516.8 579.5 4178.8 9118.6 13188.8

178.2 271.3 388.4 3264.2 4739.1 6402.7 289.9 415.4 582.9 28179.9 54870.8 89932.5 265.2 357.2 489.2 403.7 808.8 1154.3 12047.7 16729.0 19261.9 11516.4 17294.9 19934.9 12499.3 15986.6 18966.1 12503.6 16375.9 19160.5 13062.4 17235.3 20074.5 329.9 516.8 579.5 4178.8 9118.6 13188.8

178.2 271.3 388.4 3264.2 4739.1 6402.7 289.9 415.4 582.9 28179.9 54870.8 89932.5 265.2 357.2 489.2 403.7 808.8 1154.3 12047.7 16729.0 19261.9 11516.4 17294.9 19934.9 12499.3 15986.6 18966.1 12503.6 16375.9 19160.5 13062.4 17235.3 20074.5 329.9 516.8 579.5 4178.8 9118.6 13188.8

0.0 0.1 1.1 0.0 0.3 1.8 0.1 0.4 2.5 0.1 1.1 13.8 0.1 1.0 3.9 0.2 2.2 9.3 0.3 2.8 8.9 0.3 2.4 15.9 0.4 3.3 10.4 0.3 3.3 15.3 0.6 3.9 26.0 0.2 2.9 20.2 0.3 2.6 10.7

0.0 0.0 0.1 0.0 0.1 1.4 0.0 0.0 0.1 0.0 0.7 51.1 0.0 0.1 0.2 0.0 0.1 0.2 0.1 0.2 0.5 0.0 0.1 1.2 0.1 0.4 0.9 0.0 0.5 7.4 0.1 0.5 13.2 0.0 0.2 2.8 0.0 0.1 0.2

y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y

y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y

double precision was used for computing these distances. In total, 39 distinct graphs were thus generated. Values α ∈ { π3 , 23π , π , 76π } were then associated with each of the 39 graphs above, thus resulting in 156 distinct α -MSTP instances. Note, however, that α values smaller than π3 may eventually imply infeasible instances (Ackerman et al., 2013).

7.2. Numerical results Tables 4–7 present BCFX and BCFXY computational results for the instances indicated above. Each table corresponding to a different value of α . The first column in every table identifies a TSPLIB instance. The second, an α -MSTP instance originating from it, parameter n indicating its corresponding number of vertices. Every subsequent group of two columns displays a different type of performance information, respectively for BCFXY and BCFX, in that order. Namely, root node LP relaxation bound w( · ); lower an upper bounds, respectively BLB and BUB, at the end of the BC run; total CPU time t(s), in seconds; and an indication (opt ?), “y” (yes) or “n” (no), as to whether or not the algorithm solved the instance to proven optimality, under the imposed CPU time limit of two hours. An instance is considered solved (to certified optimality) if the absolute BC duality gap for it, (BUB-BLB), falls below 10−5 . Note, in this case, that corresponding rounded up values for BLB and

rounded down values for BUB, do match. If, instead, CPU time limit is hit before an optimality certificate is attained, symbol “-” replaces CPU time in the corresponding column entry. Fig. 6 displays performance profiles (Dolan and Moré, 2002) for BCFX and BCFXY. Two profile plots appear in every frame, one plot for each BC algorithm, every frame corresponding to a different value of α . Plots were generated by taking M = {BCFX, BCFXY} as the set of algorithms under comparison and I as the set of 39 corresponding α -MSTP instances, one set of instances for every different value of α ∈ { π3 , 23π , π , 76π }. The CPU time algorithm a ∈ M takes to solve instance g ∈ I is denoted by tg,a . For every algorithm a ∈ M, Fig. 6(a)–(d) plot their corresponding functions t fa (τ ) = |I1| |{g ∈ I : min{t g,a:s∈M} ≤ τ }|. g,s

Profiles were generated from the tg,a values reported in Tables 4–7. Note, however, that tg,a = 7200 is always taken whenever algorithm a fails to obtain an optimality certificate for instance g, within our imposed CPU time limit of 7200 seconds. Such a choice of tg,a value, we stress, has no impact on the conclusions reached on the relative performance of the algorithms investigated (see Dolan and Moré, 2002, for a detailed discussion on the subject). It should be noted that fa (τ ) corresponds to a certain fraction of the test instances solved by algorithm a. Namely, those that are solved in CPU times within a factor τ (τ ≥ 1) of the CPU times taken by the best performing algorithm. As an example, take τ = 10 and consider Fig. 6(c), that corresponds to the α = π case.

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775

11

Table 7 Computational results - α = 76π . Instance

n

w( · )

BLB

BUB

t(s)

opt ? y/n

eil51 eil51 eil51 berlin52 berlin52 berlin52 st70 st70 st70 pr76 pr76 pr76 eil76 eil76 eil76 rat99 rat99 rat99 kroA100 kroA100 kroA100 kroB100 kroB100 kroB100 kroC100 kroC100 kroC100 kroD100 kroD100 kroD100 kroE100 kroE100 kroE100 eil101 eil101 eil101 lin105 lin105 lin105

15 30 51 15 30 52 20 40 70 25 50 76 25 50 76 35 70 99 35 70 100 35 70 100 35 70 100 35 70 100 35 70 100 35 70 101 35 70 105

w(Fxy )

w ( Fx )

BCFXY

BCFX

BCFXY

BCFX

BCFXY

BCFX

BCFXY

BCFX

178.2 267.9 381.3 3259.2 4536.1 6225.9 278.1 404.3 573.9 27429.3 53421.1 87451.7 260.3 347.0 475.9 402.1 794.9 1136.4 11725.0 16404.7 18893.4 11353.9 16978.5 19583.3 12057.6 15516.8 18595.1 12177.0 16062.0 18845.6 12585.7 17026.8 19567.5 327.3 509.3 567.0 4075.1 9032.7 13095.6

178.2 267.9 381.3 3259.2 4536.0 6225.9 278.1 404.3 573.9 27429.3 53421.1 87451.7 260.3 347.0 475.8 402.0 794.9 1136.3 11725.0 16404.7 18893.4 11353.9 16978.5 19583.3 12057.6 15516.8 18595.1 12177.0 16062.0 18843.1 12585.7 17026.8 19567.5 327.3 509.3 567.0 4075.1 9032.7 13095.6

178.2 267.9 381.3 3259.2 4540.8 6229.3 278.1 404.3 573.9 27457.9 53449.7 87480.3 260.3 347.0 475.9 402.4 795.2 1136.7 11754.7 16410.6 18907.6 11353.9 16978.5 19591.4 12191.3 15523.8 18595.1 12177.0 16067.7 18845.6 12638.5 17062.0 19572.1 327.3 509.3 567.0 4075.1 9032.7 13095.6

178.2 267.9 381.3 3259.2 4540.8 6229.3 278.1 404.3 573.9 27457.9 53449.7 87480.3 260.3 347.0 475.9 402.4 795.2 1136.7 11754.7 16410.6 18907.6 11353.9 16978.5 19591.4 12191.3 15523.8 18595.1 12177.0 16067.7 18845.6 12638.5 17062.0 19572.1 327.3 509.3 567.0 4075.1 9032.7 13095.6

178.2 267.9 381.3 3259.2 4540.8 6229.3 278.1 404.3 573.9 27457.9 53449.7 87480.3 260.3 347.0 475.9 402.4 795.2 1136.7 11754.7 16410.6 18907.6 11353.9 16978.5 19591.4 12191.3 15523.8 18595.1 12177.0 16067.7 18845.6 12638.5 17062.0 19572.1 327.3 509.3 567.0 4075.1 9032.7 13095.6

178.2 267.9 381.3 3259.2 4540.8 6229.3 278.1 404.3 573.9 27457.9 53449.7 87480.3 260.3 347.0 475.9 402.4 795.2 1136.7 11754.7 16410.6 18907.6 11353.9 16978.5 19591.4 12191.3 15523.8 18595.1 12177.0 16067.7 18845.6 12638.5 17062.0 19572.1 327.3 509.3 567.0 4075.1 9032.7 13095.6

0.0 0.1 0.3 0.0 0.2 0.7 0.0 0.1 0.6 0.1 0.3 1.1 0.0 0.3 1.3 0.1 0.7 2.2 0.2 1.0 2.5 0.1 0.6 2.0 0.2 0.8 2.6 0.1 1.1 2.9 0.2 1.5 3.8 0.1 0.6 4.4 0.1 0.7 3.0

0.0 0.0 0.1 0.0 0.0 0.1 0.0 0.0 0.1 0.0 0.1 0.1 0.0 0.0 0.1 0.0 0.1 0.1 0.0 0.1 0.1 0.0 0.0 0.1 0.0 0.1 0.1 0.0 0.1 0.1 0.0 0.1 0.1 0.0 0.1 0.4 0.0 0.1 0.1

y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y

y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y y

It shows, for about 60% of the 39 instances, that CPU times spent by BCFXY are bounded from above by 10 times the CPU times attained by the best performing algorithm (i.e, BCFX or BCFXY, whatever applies). Algorithms with performance profile plots above the others are likely to be faster, for the class of instances considered. Let us now focus on the results in Table 4, obtained for instances with the smallest value of α we consider, i.e., α = π3 . For these instances, BCFXY LPRs, w(Fxy ), are typically much stronger than their BCFX counterparts, w(Fx ). Measured through the ra100(w(Fxy )−w(Fx )) tio , w(Fxy ) turns out to be 12.1% stronger, on the w ( Fx ) average. Accordingly, BCFXY clearly outperforms BCFX in that respect. Furthermore, BCFXY solves all 6 instances solved by BCFX, plus 8 additional ones. To solve that former group of 6 instances, the total CPU time required by BCFXY (105.5 s) is about 44 times smaller than that required by BCFX (4566.3 s). Finally, performance profile plots in Fig. 6(a) also confirm the dominance of BCFXY over BCFX for the α = π3 instances. Indeed, BCFXY also outperforms BCFX for the 25 instances left unsolved by both algorithms. At the CPU time limit, lower bounds BLB for each of these 25 instances are then stronger for BCFXY. Likewise, for 23 of them, feasible solutions returned by BCFXY are also stronger. Contrary to the results above, performance profile plots in Fig. 6(b)–(d), for instances with α ∈ { 23π , π , 76π }, indicate that BCFX is much faster than BCFXY when α > π3 . Even for the small-

est of the these three values, i.e., α = 23π , computational results already lean in favor of BCFX. Among others, in terms of the number of optimality certificates attained under the imposed CPU time limit. While BCFX solved to proven optimality 30 out of the 39 instances involved, only 26 of them were solved by BCFXY. Furthermore, for 5 out of the 7 instances left unsolved by both algorithms, BCFX found better feasible solutions, BUB, within the CPU time provided. Additionally, for all of the previous 7 instances, BLB lower bounds are also stronger for BCFX. Both algorithms solved to proven optimality all instances with α ∈ {π , 76π }, with BCFX typically running faster than BCFXY. Let us now analyse the results over the whole spectrum of α ∈ { π3 , 23π , π , 76π } values we considered. While w(Fxy ) was significantly stronger than w(Fx ) for α = π3 (i.e., 12.1% stronger, on the average), the two lower bounds are nearly equally as strong for the remaining cases. More precisely, the average ratios 100(w(Fxy )−w(Fx )) attained for α = 23π and α = π are quite small, of w (F ) x

0.07% and 0.02%, respectively. For the α = 76π case, the average ratio is even smaller, i.e., less than 0.01%. Accordingly, from a lower bound perspective, Fxy has no significant advantage over Fx for the α ≥ 23π values we consider. Furthermore, the number of binary 0-1 variables in formulation Fx is just one third of that for Fxy . Additionally, α -AC enforcing cover inequalities are always treated as cutting planes. By contrast, constraints (11) and (12), that play a similar role for Fxy , become

12

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775

Fig. 6. Performance profiles for algorithms

BCFXY and BCFX for different values of α ∈ { π3 , π , π3 , 76π }.

denser as α grows. As a result, average CPU time spent per BC enumeration node becomes smaller for BCFX. Consequently, within a same CPU time span, significantly more BC nodes are investigated by BCFX and heuristic G-α -ST is therefore more frequently called by it. Precisely for these two reasons, BCFX also dominates BCFXY for those instances left unsolved by both algorithms. Take the α = 23π case, for example, and note that although w(Fxy ) is at least as strong as w(Fx ) for these instances, better lower bounds BLB are attained by BCFX at the CPU time limit. Note that bounds w(Fx ) and w(Fxy ) come quite close to one another for our α ≥ 23π instances. Accordingly, the computational impact of BCFX not carrying out fractional point separation of cover inequalities for sets s(i), |s(i)| ≥ 4, for our α ≥ π instances, is apparently very marginal. As our computational experiments indicate, α -MSTP becomes easier to solve as α values grow. Indeed, more optimality certificates are then attained by both algorithms. One of the reasons −w(· )) is that LP duality gaps, ( 100(BUB ), are substantially reduced BUB with the increase of α . Considering only those instances solved to proven optimality by BCFX, average duality gaps are 11.0%, 4.7%, 0.4% and 0.1%, respectively for α equal to π3 , 23π , π , and 76π . Instances with α = π3 are, by far, the hardest for our algorithms to solve. To highlight this point, let us analyse the 24 of them that remained unsolved by both algorithms. In particular, consider their −BLB ) CPU time limit duality gaps, ( 100(BUB ). They are, on the averBUB

age, 8.8% for BCFXY and 17.5% for BCFX, what is considerably high in both cases. Still, to attain these gaps, heuristic G-α -ST proved particularly important. Indeed, it was then our only source for feasible α -MSTP solutions. No LP relaxation at any BC tree node resulting naturally integral for these instances. The best-first search strategy implemented by our BC algorithms being to blame, in that case. Accordingly, it is fair to say that the α -MSTP upper bounds displayed in Table 4 originate, predominately, from G-α -ST. As we indicated before, G-α -ST is not guaranteed to always finding a feasible α -MSTP solution. Therefore, we address how effective it was, in practice, in attaining that goal, right before the BC algorithms are called. In Table 8, we display some G-α -ST computational results. For each instance and value of α , the table gives the initial upper bound (IUB) found by G-α -ST, the best upper bound (BUB) known, obtained by either BCFX or BCFXY, and the  −BU B corresponding gap (in percentage values), defined as 100 IU BBUB . IUB is precisely the upper bound used to initialize the BC algorithms BCFX and BCFXY. If the heuristic fails to find any IUB, ∞ appears in the corresponding entry. G-α -ST computational times are not presented in the table because they are negligible, always less than 0.01 s. Note that for just two instances the heuristic failed to find an α -ST. As one would expect, G-α -ST solution quality increases with the increase of α . Average values attained were approximately 7.7%, 3.4%, 1.6% and 0.8% away from BUB averages, for α = π3 , 23π , π and 76π , respectively. In summary, we believe that

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775

13

Table 8 Computational results obtained by the α −MSTP heuristic G-α -ST. IUB is the cost of the heuristic solution and BUB is the best known upper bound. Instance

eil51 eil51 eil51 berlin52 berlin52 berlin52 st70 st70 st70 pr76 pr76 pr76 eil76 eil76 eil76 rat99 rat99 rat99 kroA100 kroA100 kroA100 kroB100 kroB100 kroB100 kroC100 kroC100 kroC100 kroD100 kroD100 kroD100 kroE100 kroE100 kroE100 eil101 eil101 eil101 lin105 lin105 lin105

n

15 30 51 15 30 52 20 40 70 25 50 76 25 50 76 35 70 99 35 70 100 35 70 100 35 70 100 35 70 100 35 70 100 35 70 101 35 70 105

α=

π

α=

3

2π 3

α=π

α=

7π 6

IUB

BUB

gap (%)

IUB

BUB

gap (%)

IUB

BUB

gap (%)

IUB

BUB

gap (%)

299.9 437.2 636.4 6056.5 ∞ 10586.3 436.8 695.1 1007.3 ∞ 90469.3 147087.1 437.3 596.0 815.7 653.4 1289.6 1961.6 21852.8 28525.7 32216.9 20150.5 29213.8 33094.5 19131.8 26088.1 32074.2 21479.9 27765.7 31897.8 21097.1 28773.2 33128.7 562.2 855.3 964.2 6164.2 14016.8 20316.8

258.1 409.9 587.2 5262.7 7090.3 9896.8 420.5 628.6 909.9 41809.0 83009.5 139876.8 395.0 551.3 756.1 633.2 1268.8 1845.5 19707.6 26615.5 31309.7 17263.2 27334.7 30250.1 18639.7 25111.2 31366.9 19250.4 25309.5 30751.6 20012.2 26553.4 30774.2 499.6 791.7 906.0 5751.0 13209.7 19388.1

16.2 6.7 8.4 15.1 7.0 3.9 10.6 10.7 9.0 5.2 10.7 8.1 7.9 3.2 1.6 6.3 10.9 7.2 2.9 16.7 6.9 9.4 2.6 3.9 2.3 11.6 9.7 3.7 5.4 8.4 7.7 12.5 8.0 6.4 7.2 6.1 4.8

219.8 340.0 453.0 4361.8 5561.2 7988.8 341.9 511.3 696.3 36450.0 74763.4 119636.5 321.5 450.9 588.7 509.0 991.0 1397.7 16017.1 21130.4 24123.7 15088.5 21948.3 24366.2 14776.5 19706.7 23860.2 15317.9 20746.5 23698.3 14991.1 21024.0 24522.9 409.1 628.8 695.4 4957.9 11123.5 15673.8

209.3 315.3 436.5 4307.2 5347.2 7488.7 340.3 505.1 687.6 35435.9 71920.1 115589.2 307.9 423.6 565.1 504.1 970.2 1367.1 15214.5 20352.8 23976.5 14192.2 20817.2 23565.3 14709.1 19383.9 23558.5 14475.2 20003.6 22617.4 14937.0 20180.5 23825.3 384.3 605.3 669.2 4888.5 10834.0 15376.3

5.0 7.8 3.8 1.3 4.0 6.7 0.5 1.2 1.3 2.9 4.0 3.5 4.4 6.4 4.2 1.0 2.1 2.2 5.3 3.8 0.6 6.3 5.4 3.4 0.5 1.7 1.3 5.8 3.7 4.8 0.4 4.2 2.9 6.4 3.9 3.9 1.4 2.7 1.9

179.0 272.2 393.9 3267.0 4749.5 6443.5 298.2 420.5 590.9 29428.2 56137.8 91685.8 268.3 359.2 497.5 410.4 824.1 1176.0 12144.9 16948.5 19676.0 11869.2 17472.4 20165.2 13033.4 16294.2 19311.8 12687.8 16531.6 19432.4 13081.4 17496.5 20533.0 338.3 526.0 591.9 4283.8 9180.8 13300.5

178.2 271.3 388.4 3264.2 4739.1 6402.7 289.9 415.4 582.9 28179.9 54870.8 89932.5 265.2 357.2 489.2 403.7 808.8 1154.3 12047.7 16729.0 19261.9 11516.4 17294.9 19934.9 12499.3 15986.6 18966.1 12503.6 16375.9 19160.5 13062.4 17235.3 20074.5 329.9 516.8 579.5 4178.8 9118.6 13188.8

0.5 0.4 1.4 0.1 0.2 0.6 2.9 1.2 1.4 4.4 2.3 1.9 1.2 0.6 1.7 1.7 1.9 1.9 0.8 1.3 2.1 3.1 1.0 1.2 4.3 1.9 1.8 1.5 1.0 1.4 0.1 1.5 2.3 2.5 1.8 2.1 2.5 0.7 0.8

179.0 267.9 381.3 3259.2 4565.0 6248.4 278.1 414.0 579.7 27579.3 53694.6 87848.6 260.6 350.0 482.0 409.7 804.0 1151.8 12031.2 16660.5 18922.7 11353.9 17166.4 19828.2 12195.2 15554.3 18645.9 12412.2 16197.6 19197.3 12776.6 17185.7 19736.1 328.9 510.9 573.7 4080.2 9059.3 13121.6

177.7 267.1 380.1 3259.2 4505.3 6229.3 278.1 404.3 573.8 27457.0 53448.8 87479.4 259.5 346.9 475.6 402.3 795.2 1136.6 11754.7 16402.6 18907.6 11353.9 16978.5 19563.4 12191.3 15523.8 18595.1 12177.0 16067.7 18845.6 12638.5 17062.0 19572.1 327.3 508.2 566.0 4075.1 9032.7 13095.6

0.7 0.3 0.3 0.0 1.3 0.3 0.0 2.4 1.0 0.4 0.5 0.4 0.4 0.9 1.3 1.8 1.1 1.3 2.4 1.6 0.1 0.0 1.1 1.4 0.0 0.2 0.3 1.9 0.8 1.9 1.1 0.7 0.8 0.5 0.5 1.4 0.1 0.3 0.2

G-α -ST achieved the goals it was designed for, i.e., quickly providing good quality α -MSTP upper bounds, either at initialization or throughout the execution of our BC algorithms. 7.2.1. Computational results for sparse instances Wireless communication networks may well be sparse, for instance, due to the existence of obstacles in between two given problem locations. In that case, the edge between the corresponding pair of vertices may eventually be absent from G. In accordance with that and despite the fact that α -MSTP was originally posed over complete Euclidean graphs, we also conducted computational experiments for sparse instances of the problem. In what follows, we discuss the main conclusions we reached out of them. We do not present results in detail, in order to keep the paper at a manageable size. We generated sparse α -MSTP instances from the same TSPLIB instances we generated our dense α -MSTP ones. More than that, sparse instances were generated directly from the complete graphs corresponding to our dense α -MSTP instances. Namely, out of each different complete graph involved, two sparse graphs were generated. One, with an acceptance probability of 0.15 applied to every complete graph edge, thus aiming at obtaining a resulting sparse graph Gd with density d = 15%. The other, with an edge acceptance probability of 0.25, thus aiming at a resulting sparse graph Gd of density d = 25%. These instances can be obtained from the authors, on request.

In total, 39 sparse instances were tested, for each value of d ∈ {0.15, 0.25} and α ∈ { π3 , 23π , π , 76π }. All d = 25% instances turned out to be feasible. For the d = 15% instances, however, some of them are infeasible. Namely, 6, 4, 3 and 3 instances, respectively for α values of π3 , 23π , π , and 76π . While BCFX always detected infeasibility right at the root node, BCFX had to resort to branching in one out of the 16 infeasible instances, in order to do so. For just 2 out of the 33 feasible instances corresponding to d = 15% and α = π3 , BCFXY failed to provide optimality certificates, within the imposed CPU time limit of two hours. For all other 31 instances, BCFXY succeeded in that task. Success rates for BCFX were similar for α ∈ { 23π , π , 23π } and d ∈ {0.15, 0.25}, but much worse than BCFXY’s for α = π3 . More specifically, BCFX left 13 (resp. 17) unsolved α = π3 instances for d = 0.15 (resp. d = 0.25). Overall, considering the performance of our BC algorithms, sparse instances do not appear to bring any relevant new information over what we already acquired from dense instances experiments. The previous relative performance of the BC algorithms is not changed for sparse instances and the practical difficulty in solving them complies with the previously observed pattern for dense instances. Namely, instances with α < π remain the hardest to solve to proven optimality. Likewise, BCFXY performs better than BCFX for α < π instances while the reverse is true for the α ≥ π ones. The only noticeable difference is that the performance advantage BCFX attains over BCFXY, for α ≥ π instances,

14

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775

is relatively smaller here than it was before, for dense instances. The reason for that is the positive impact sparse instances bring to Fxy density. In particular, its α -AC imposing constraints become a lot sparser. That fact, in turn, allows its corresponding LPRs to be solved much faster than its dense instance counterparts. As for heuristic G-α -ST, the main points to highlight are the following. As one would expect, G-α -ST performance is worse for low density instances and smaller values of α . Under these settings, whenever G-α -ST fails to find a feasible solution, the corresponding instance is most of the time feasible, as attested by our BC algorithms. In that case, it is likely that G-α -ST found feasible solutions for it, throughout the BC enumeration process, at nodes other than the root. In any case, for instances with intended d = 15% densities, feasible solutions were obtained, for our single trial G-α -ST run, for 8, 29, 32 and 34 instances, respectively for α values of π3 , 23π , π , and 76π . Corresponding figures for the intended d = 25% instances are 19, 37, 38 and 38. In relation with the conclusions above, it is important to remark that our BC algorithms implement a best-first-search strategy, which is known to be very prone on not generating Linear Programming relaxations that are naturally integer feasible. Accordingly, it is thus very likely that G-α -ST indeed succeeded in finding feasible α -MSTP solutions at BC nodes other than the root. Therefore, the fact that it might eventually fail to find a root node feasible solution is not, just by itself, a serious impairment, given that G-α -ST operates at every BC node. Finally, when G-α -ST succeeds in finding a feasible α -MSTP solution at the root node of the BC tree, the resulting gaps are higher (sometimes twice as high) than those previously attained for corresponding dense instances. For d = 15% and α values of π3 , 23π , π , and 76π , they are respectively 14.4%, 4.6%, 2.1%, and 1.8%. Corresponding figures for d = 25% are respectively 13.4%, 5.7%, 1.9%, and 1.6%.

Computational experiments indicate that the BC algorithm for Fxy , i.e., BCFXY, is our option of choice for instances with small values of α . The LP bounds attained by Fxy are, in that case, substantially stronger than those for Fx . Furthermore, that advantage also translates into BCFXY running faster than its counterpart for Fx , i.e., BCFX. For larger values of α , however, that choice is reversed. LP relaxation bounds for Fxy and Fx then tend to be, virtually, equally strong. Furthermore, α -AC enforcing constraints tend to be denser in the BCFXY case, with the increase of α . That, in turn, positively impacts on the CPU time required for computing LP relaxation bounds w(Fx ) and BCFX is then pushed into the lead. Various interesting investigation topics are open for future research. For the α ≥ π case, computational results suggest that only very marginal benefits should be expected from exact separation of larger cover inequalities, i.e., those involving four or more edges. However, establishing the computational complexity for that separation is certainly an interesting open problem. Another attractive proposition is to attempt to characterize additional valid inequalities for the problem. In that respect, we plan to further investigate the projection of Fxy onto the space of the edge variables, x. More specifically, we aim at hopefully characterizing the entire set of families of valid inequalities one may generate from the projection cone rays. The idea appears promising since we not only have theoretical guarantees that Fxy is at least as strong as Fx , but our computational experiments also show that it may actually be significantly stronger, at least for instances with small (and therefore more challenging) values of α . Accordingly, an enhanced BC algorithm might arise from this proposed investigation. Finally, the IP formulations proposed here may be adapted for modeling and solving several other interesting and closely related network design problems. To name just one, we refer the reader to the Orientation and Power Assignment Problem (Aschner and Katz, 2017).

8. Conclusions

Appendix A. Separating inequalities (21) by solving integer programs

Our investigation appears to be the very first devoted to IP formulations and exact algorithms for α -MSTP. This is a combinatorial optimization problem with a strong computational geometry flavour, that finds applications in the design of wireless networks operating under directional antennas. The key issue addressed here is how to represent α -MSTP angular constraints, i.e, the fundamental defining property for the problem. Two alternatives are suggested for defining these α -AC, each of them leading to a distinct linear integer programming formulation of the problem. BC algorithms for solving these formulations are then proposed and tested. The algorithms are accompanied by a heuristic, G-α -ST, that aims at finding feasible α -MSTP solutions. One of our formulations, Fxy , involves two sets of decision variables. One, x ∈ Bm , to select the spanning tree edges. The other, y ∈ B2m , to enforce the angular constraints. The other formulation, Fx , only uses the natural space variables, x. Additionally, it relies on combinatorial arguments to characterize subsets of input graph edges as α -MSTP admissible or not. That characterization translates into exponentially many cover inequalities, that Fx then uses to enforce α -AC. Formulation Fxy was analytically shown to be at least as strong as Fx , in LP relaxation bound terms. As an offspring to our main investigation, we disclosed the fact that α = π is a critical value for formulation Fx . In particular, for α < π , we identified a key O(n3 ) subset of the exponentially many cover inequalities. Namely, one that suffices to attain not only an α −MSTP formulation but also a polyhedron that is exactly the same as that for the full blown formulation Fx , containing all cover inequalities. Remaining cover inequalities, in this case, being redundant and therefore having no polyhedral description effect.

So far we failed to obtain any polynomial time algorithm for the exact separation of inequalities (21), for the α ≥ π case. Accordingly, we propose, next, a computationally intensive procedure to do so. More efficient separation procedures are likely to exist for (21). Ours, however, serves as an initial contribution towards formalizing and solving the problem. As one shall see, our procedure requires the solution of O(n) sequences of integer programs, one sequence for every applicable vertex i ∈ V. In order to describe it, let x be the point to be separated and assume that our corresponding support graph for it, G = (V, E ), is such that δ (i ) := {e ∈ δ (i ) : xe > 0} applies. Additionally, define Z i := {ij ∈ Ri : e = {i, j} ∈ δ (i )}. For every i ∈ V and K ∈ {2, . . . , |δ (i )|}, we formulate an integer program denoted by (SEP)i,K . If (SEP)i,K is feasible, it returns a subset s(i ) ⊆ δ (i ) of cardinality K which is most likely to violate an inequality (21), for the conditions imposed. Otherwise, if it is infeasible, it serves as a proof that no such violated inequality exists. Accordingly, solving (SEP)i,K for every i ∈ V and all values of

K ∈ {2, . . . , |δ (i )|}, suffices to separate cover inequalities (21), exactly. Every problem (SEP)i,K is formulated with binary 0-1 deci-

sion variables z = {zek ∈ B : e = {i, j} ∈ δ (i ), k = 1, . . . , K }. In particular, zek = 1 applies if and only if • •

edge e is included in the set s(i) one is looking for, and angle ࢭ0ij , defined by the unit vector ij associated with e = {i, j} ∈ s(i ), is the kth smallest angle amongst all similarly defined s(i)-edge angles.

A.S. da Cunha and A. Lucena / Computers and Operations Research 112 (2019) 104775

In order to illustrate the explanations above, consider the previously defined variables z and assume that δ (i ) is the set of 5 edges depicted in Fig. 1. Assume as well that K is fixed to 3 and that one wishes to identify, through z, the subset s(i ) = {{i, z}, {i, u}, {i, l }}, all of them assumed to belong to δ (i ). Accordingly, given that ࢭ0iz < ࢭ0iu < ࢭ0il and K = |s(i )| = 3 hold, the only non-null vari1 = z 2 = z 3 = 1. ables corresponding to s(i) must necessarily be ziz iu il Problem (SEP)i,K is thus defined as:

(SEP)i,K

max

K 



(1 − xe )zek

(A.1)

k=1 e={i, j}∈δ (i ) K 

zek ≤ 1

e = {i, j} ∈ δ (i )

zek = 1

k = 1, . . . , K



zkf +1

(A.2) (A.3)

e = {i, j} ∈ δ (i )

(A.4)

k = 1, . . . , K − 1 e = {i, j} ∈ δ (i )

(A.5)

f ={i,u}∈δ (i )\{e}:ui j >α



ze1 ≤

zKf

f ={i,u}∈δ (i )\{e}: jiu >α

zek ∈ B

e = {i, j} ∈ δ (i ) k = 1, . . . , K

that there exists no subset of K non-admissible edges in δ (i ). Let us assume otherwise. Namely, that (SEP)i,K is feasible and that an optimal solution for it, z ∈ BK |δ (i )| , is such that K 



(1 − xe )zke ≥ 1

k=1 e={i, j}∈δ (i )

e={i, j}∈δ (i )

zek ≤

enforce that e = {i, j} cannot be assigned an index k = 1, with ze1 = 1 thus holding, if the following additional condition is not enforced. Namely, that the circular sector obtained by rotating ij anticlockwise around i by α radians does not contain the Kth indexed edge, f, for which zKf = 1 would then hold. Program (SEP)i,K may eventually be infeasible, thus implying

z(x )i,K :=

k=1



15

holds. In that case all inequalities (21), for sets s(i ) ⊆ δ (i ) with exactly K edges, are satisfied by x. However, if z(x )i,K < 1 holds,

set s(i ) := {e = {i, j} ∈ δ (i ) : Kk=1 ze = 1} would then imply an inequality (21) violated by x. Finally note that if xe = 1 holds for every e ∈ δ (i ), any feasible solution for (SEP)i,K will have a same objective function value z(x )i,K = 0. Accordingly, any of them would imply a violated inequality (21). In that case, for such a particular i ∈ V, the separation of inequalities (21) would not require solving the whole sequence of problems (SEP)i,K : K = 2, . . . , |δ (i )|. It would be enough, in this case, to order the edges in δ (i ) according to (16) and check if condition (17) applies to the entire set δ (i ).

(A.6)

Constraints (A.2) impose that every edge in δ (i ) must be assigned to at most one index in {1, . . . , K }. In turn, constraints (A.3) enforce that every index k ∈ {1, . . . , K } must be assigned to exactly one δ (i ) edge. Accordingly, exactly K edges of δ (i ) must therefore be included in the set s(i) one is aiming for. Constraints (A.4) and (A.5) enforce that the edges of s(i), sorted according to (16), must satisfy conditions (17), that characterize sets of i-incident non-admissible edges. More precisely, assume that ࢭ0ij is the kth (k < K) smallest angle in {0il : l ∈ {1, . . . , K }, {i, l } ∈ s(i )} and that zek = 1, e = {i, j}, thus hold. Accordingly, e must not be within the circular sector of α radians and ini , where f = {i, u} ∈ s(i ) is tial extreme ray collinear with a vector iu further defined as follows. Namely, f must be such that ࢭ0iu is the (k + 1 )th smallest angle amongst those corresponding to s(i) edges and zkf +1 = 1 thus holds. Such a condition is enforced by (A.4). Constraints (A.5), in turn, model similar conditions applying to a pair of edges e = {i, j} and f = {i, u}, e and f belonging to s(i) and such that ze1 = zKf = 1 holds. Note, in this case, that constraints (A.5)

References Ackerman, E., Gelander, T., Pinchasi, R., 2013. Ice-creams and wedge graphs. Comput. Geom. 46, 213–218. Aschner, R., Katz, M.J., 2017. Bounded-angle spanning tree: modeling networks with angular constraints. Algorithmica 77, 349–373. Carmin, P., Katz, M.J., Lotker, Z., Rosén, A., 2011. Connectivity guarantees for wireless networks with directional antennas. Comput. Geom. 44, 477–485. Cormen, T., Leiserson, C., Rivest, R., 2009. Introduction to Algorithms, 3rd ed. MIT Press. Dolan, E.D., Moré, J.J., 2002. Benchmark optimization software with performance profiles. Math. Program. 91 (2), 201–213. Edmonds, J., 1971. Matroids and the greedy algorithm. Math. Program. 1 (1), 127–136. Kruskal, J., 1956. On the shortest spanning subtree of a graph and the traveling salesman problem. Proc. Am. Math. Soc. 7, 48–50. Padberg, M., Rinaldi, G., 1991. A branch-and-cut algorithm for the resolution of large-scale symmetric traveling salesman problems. SIAM Rev. 33 (1), 60–100. Padberg, M.W., Wolsey, L.A., 1983. Trees and cuts. Ann. Discrete Math. 17, 511–517. Reinelt, G., 1991. TSPLIB - a traveling salesman problem library. ORSA J. Comput. 3 (4), 376–384. XPRESS mixed integer optimization package, release 8.4.2017. FICO XPRESS. Yu, Z., Teng, J., Bai, X., Xuan, D., Jia, W., 2014. Connected coverage in wireless networks with directional antennas. ACM Trans. Sens. Netw. 10 (3), 51:1–51:28.