Tight bounds from a path based formulation for the tree of hub location problem

Computers & Operations Research 36 (2009) 3117 -- 3127

Contents lists available at ScienceDirect

Computers & Operations Research journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / c o r

Tight bounds from a path based formulation for the tree of hub location problem Ivan Contreras,a ∗ , Elena Fernández a , Alfredo Marín b a b

Dpt. d' Estadística i Investigació Operativa, Universitat Politècnica de Catalunya, Barcelona, Spain Departamento de Estadística e Investigación Operativa, Universidad de Murcia, Murcia, Spain

A R T I C L E

I N F O

Available online 25 December 2008 Keywords: Hub location Spanning trees Lagrangean relaxation

A B S T R A C T

This paper considers the tree of hub location problem. We propose a four index formulation which yields much tighter LP bounds than previously proposed formulations, although at a considerable increase of the computational burden when obtained with a commercial solver. For this reason we propose a Lagrangean relaxation, based on the four index formulation, that exploits the structure of the problem by decomposing it into independent subproblems which can be solved quite efficiently. We also obtain upper bounds by means of a simple heuristic that is applied at the inner iterations of the method that solves the Lagrangean dual. As a consequence, the proposed Lagrangean relaxation produces tight upper and lower bounds and enable us to address instances up to 100 nodes, which are notably larger than the ones previously considered in the literature, with sizes up to 20 nodes. Computational experiments have been performed with benchmark instances from the literature. The obtained results are remarkable. For most of the tested instances we obtain or improve the best known solution and for all tested instances the deviation between our upper and lower bounds, never exceeds 10%. © 2009 Published by Elsevier Ltd.

1. Introduction Hub location is an area that is receiving increasing attention in discrete location due both to the variety and economic impact of potential applications, and to the theoretical interest of the different formulations for these problems. Hub location problems have their origin in transportation systems where several sites send and receive some product and where it is convenient to use transshipment points (hubs), where the product is collected and distributed to reduce transportation costs. A common feature of hub location problems is that a discount factor is applied to the usage of the links that connect any pair of hubs. Traditionally, many hub location models assume that the hubs are fully interconnected, that is to say, that there exists a link connecting any pair of hubs, which can be used by applying the corresponding discount factor. However, it is widely accepted that there exist many applications in which the backbone network (i.e. the network connecting the facilities) is not fully interconnected (see, for instance the works of O'Kelly and Miller [1] or Klincewicz [2]). Moreover, some hub location models have been proposed recently where two hubs are not necessarily connected by some link (see, for

∗ Corresponding author. Tel.: +34 934 502 997. E-mail addresses: [email protected] (I. Contreras), [email protected] (E. Fernández), [email protected] (A. Marín). 0305-0548/$ - see front matter © 2009 Published by Elsevier Ltd. doi:10.1016/j.cor.2008.12.009

instance [3,4]). For example, the so-called hub-arc location models [5,6] impose no condition on the structure of the arcs that connect hubs, and they do not even require these arcs to define one single connected component. Another example of a hub location model not requiring full interconnection between hubs is the tree of hub location problem (THLP) that was recently proposed in Contreras et al. [7], which is the problem that we address in this work. It is a single-allocation hub location problem where a fixed number of hubs have to be located on a network, with the particularity that it is required that the hubs are connected by means of a (non-directed) tree. The THLP is defined on a directed graph where it is assumed that for each pair of nodes there exists a known amount of flow that must be sent through the network. To this end p hubs must be located on the network and connected by means of a tree. Each node must be allocated to one single hub and all the flow from/to the node must leave/enter it through its allocated hub. Excepting the arcs that connect each node with its allocated hub, the only arcs that can be used for routing the flows must be links connecting hubs. There is a unit transportation cost associated with each arc. As usual, if hubs are located at both end-nodes of the arc, a discount factor is applied to the unit cost. The objective is to minimize the operation costs of the system, which depend on the amount of flow that circulates through each arc, and on the links to which the discount factor is applied. Potential applications of models where facilities are connected by means of a tree arise when the cost of the links between

3118

I. Contreras et al. / Computers & Operations Research 36 (2009) 3117 -- 3127

facilities is very high, and as a consequence full interconnection is prohibitive. For sending all the flows in the network a path must exist between each pair of facilities, i.e., a connected graph must be built. But due to the high cost of the connections, connectivity must be achieved using the minimum number of links. Specific applications of such problems arise mostly in telecommunications [8,9] and in transportation (see, for instance, the recent work of Nguyen and Knippel [10] for an excellent description of the practical relevance of tree-backbone problems in small package delivery). One concrete example of an application of the THLP is the design of the high-speed train network in Spain, which is currently under construction and which is planned to be completed by 2020. This train network has been designed with a tree structure and it is intended that, when finished, each city with more than 10,000 inhabitants will be within 50 km of some high-speed train station. Another application is the design of rapid transit systems for urban areas where there are flows of users (citizens) who must travel between origin/destination pairs. More often than not, such systems are designed with a tree structure, where the users are allocated to the closest station (hub). The TLHP is an NP-hard problem which involves location, design, and routing decisions. The location decisions refer to where to locate the hubs on the network. The design decisions refer to how to interconnect the selected hubs so as to define a tree, and to the (single) allocation of the nodes to the selected hubs. The routing decisions refer to how to route the flows between vertices of the network through the tree of hubs. Note that in the THLP the location and the design decisions determine the routing decisions, since for each pair of nodes of the network there exists only one path connecting them in the tree of hubs. In terms of the classification of O'Kelly and Miller [1], the THLP falls within the so-called Protocol B which contains hub location problems with single allocation of nodes to hubs, and hub-to-hub links. In terms of the classification of Klincewicz [2], the THLP is a tree-star hub location problem, in which the backbone network is a tree and the single allocation of non-hub nodes to hubs defines a set of stars. Different problems where the backbone network consists of a tree have been studied in the literature. Among them we can mention the optimum communication spanning tree problem (OCSTP) [8,11] and the minimum sum violation tree problem (MSVTP) [10]. While the objective in the MSVTP is to minimize the sum of the violations of the delivery due dates in transportation systems, in the OCSTP the objective is to minimize the operation costs of the system which coincides with that of the THLP. In fact, the OCSTP is a subproblem of the THLP, since the latter can be easily transformed into an OCSTP if the location of the hubs and the allocation of non-hubs to hubs are known. In addition to the above problems and to other hub-location problems, the THLP is also related to other types of problems where location, design, or routing decisions are considered. To mention just a few: the location of tree-shaped facilities (also called extensive facilities) on a network [12–14]; rapid transit network design problems that locate lines on rail transit systems [15]; and some locationrouting problems, where a route is located and customers are assigned to the “closest” node in the route, like the ring-star problem [16,17]. The similarities between the THLP and the well-known single allocation p-hub location problem (p-HLPSA) are the following: (i) both problems are hub location ones where a discount factor is applied to the links that connect a pair of hubs; (ii) both problems locate exactly p hubs; and (iii) in both problems single allocation of non-hub nodes to hubs is required. The main dissimilarity is, however, that in the p-HLPSA the backbone is a complete network, whereas in the THLP the backbone network consists of a tree as explained above. In turn, this basic distinction implies other important differences in the

structure of the feasible solutions that affect to all relevant design and solution methods. In particular, it is important to recall that in optimal solutions to the p-HLPSA the path that sends the flow between two nodes contains at most three links connecting at most two consecutive hubs. However, in the THLP the subgraph associated with any feasible solution contains one single path between each pair of nodes, which may consist of up to p − 1 links connecting p consecutive hubs. In Contreras et al. [7] the THLP was introduced, and a three indices formulation was presented and computationally tested in terms of its LP bound and of the effort required by a standard optimization solver to solve it. The proposed formulation was successful in solving small instances of the THLP up to 40 nodes. The lower bounds that it produced were not tight enough so as to address efficiently the optimal solution to larger instances. Observing the quality of the LP bounds produced by the different formulations for classical single allocation hub location problems, it is clear that the best LP bounds are usually obtained with four indices formulations (see, for instance, the works by Campbell [18], by Skorin-Kapov et al. [19], or by Ernst and Krishnamoorthy [20]). However, in practice, three indices formulations have been preferred to four indices formulations, since the latter are not always manageable due to the enormous number of variables that they involve. One possibility for obtaining the tight LP bounds associated with four indices formulations avoiding the computational burden implied by the solution of the LP relaxation is to resort to Lagrangean relaxation. It is well known that Lagrangean relaxation produces lower bounds at least as tight as the LP ones (see, for instance [21]). The potential advantage of Lagrangean relaxation is that depending on the relaxation it may be possible to decompose the Lagrangean function into easier smaller subproblems which can be easily and efficiently solved, thus allowing to obtain tight lower bounds with a manageable computational effort for larger instances. Recently, this methodology has been successfully used by Contreras et al. [22,23] who, by means of a Lagrangean relaxation of a four indices formulation, have optimally solved instances up to 200 nodes for the capacitated hub location problem with single assignment. In this work we present a formulation of the THLP that uses two sets of two indices binary variables, and one set of four indices binary variables. The two indices binary variables in the first set, z, define the allocation of nodes to hubs, whereas the two indices binary variables of the second set, y, indicate whether or not a pair of nodes corresponds to a pair of adjacent hubs in the tree of hubs. The set ij of four indices binary variables xkm indicate whether or not the flow emanating from node i with destination in node j is routed via two adjacent hubs m and k. We propose a Lagrangean relaxation of the formulation that allows decomposing the Lagrangean function into two independent subproblems, one in terms of the z variables, and the second one in terms of the y and x variables. Moreover, we show that each of these independent subproblems can be expressed in terms of only the two indices variables and we indicate how to obtain their corresponding solutions. We also present a fast heuristic to obtain good quality feasible solutions to the problem, that is applied at every iteration of the subgradient optimization that is used for solving the Lagrangean dual. The numerical results obtained in a series of computational experiment with instances up to 100 nodes indicate that the four index formulation produces tighter LP bounds than those obtained with the previous three index formulations of Contreras et al. [7]. The rest of this work is structured as follows: Section 2 describes the THLP and the four index formulation. Section 3 presents the proposed Lagrangean relaxation, and analyzes the structure of the subproblems and their solutions. Section 4 presents the proposed primal heuristic to obtain feasible solutions. The description of the computational experiments and the analysis of the obtained results are

I. Contreras et al. / Computers & Operations Research 36 (2009) 3117 -- 3127

given in Section 5. Finally, in Section 6 we present some conclusions and remarks.

Path for sending the flow from i to j (i≠ j)

2. The THLP

i non-hub

Consider a complete digraph without loops G = (N, A) whose set of nodes, N = {1, . . . , n}, represents the set of origins and destinations of a certain product that is routed through G via some transshipment nodes that are called hubs. When some of the end-nodes of the arc (r, s) ∈ A is not a hub, each unit of product that traverses (r, s) incurs a cost crs  0, whereas when both r and s are hubs a discount factor 0    1 is applied, and the unit cost associated with arc (r, s) is crs . Triangle inequality of costs cij is not assumed. We do not assume symmetry either, i.e., costs cij and cji can be different. For each pair i, j ∈ N, Wij  0 denotes the flow from origin i to  destination j. For each vertex i, Oi = j Wij denotes the total flow  emanating from i, whereas Di = j Wji denotes the total flow with destination in node i. Without loss of generality we assume that the graph of flows between nodes is connected (see [7]). Any node of N can be chosen to become a hub, and there is a fixed number p of nodes of N, 2  p  n − 1, that must be chosen to be hubs. For each pair i, j ∈ N, if neither i nor j are hubs the flow Wij must go from i to j through one or more intermediate hubs. When i or j are hubs Wij can go directly from i to j, although it is also possible that Wij uses one or more intermediate hubs. No node which is not a hub can be used to transship the product. We also require single allocation and, thus, every node i which is not a hub must be assigned (allocated) to one single hub, so that all the flow from/to i to/from any other node must pass through this hub. We require the hubs to be connected by means of a (non-directed) tree. Observe that since the links that connect hubs are required to define a tree (that we denote small tree), the structure that results after adding to the small tree the links connecting non-hub nodes to their allocated hub nodes defines a tree in G (that we denote large tree). Thus, in the large tree the path between each pair i and j will be unique. In this path all the transshipment nodes will be hubs so all the arcs, excepting (possibly) the first and last ones, link hubs. The cost per unit of flow of the path between a pair i and j is the sum of the costs of the arcs in the path, where the discount factor is applied to all the inter-hub links. We want (i) to locate p hubs; (ii) to link them so as to define a tree; and (iii) to allocate every non-hub node to one single hub, in such a way that the overall cost be minimized. The cost is the sum of the routing costs associated with the flows between all the elements in N × N. Given (i, j) ∈ N × N, the routing cost of the flow Wij is given by the cost per unit of flow of the unique path from i to j in the large tree which traverses at least one hub times the amount of flow Wij . We next present the formulation for THLP that we will use throughout this work. All the variables in the formulation are binary. The location of hubs as well as the allocations of non-hubs to hubs, are modeled by means of the following z-variables:  zik =

1 0

if the non-hub i is allocated to the hub k otherwise

∀i, k, i  k.

When i=k, variable zkk represents whether or not a hub is located at node k. The z variables will be referred to as location/allocation variables. The links that connect hubs in the small tree are modeled by means of the following y variables  ykm =

1 0

if the edge {k, m} links hubs m and > k otherwise

∀k, ∀mk.

3119

zim

j non-hub

i non-hub

m1

i 1

i

j hub

i

i hub

i

i hub

r-1

2

1

mr xi jm

xi j m m

mr

mr xi j m

xi j m m

2

1

r-1

1

2

j zj m r

mr

mr

m2

xi j m m

j xi j mr j

m2

m1 x i jim1

j hub

m2

m1 xijim1

j non-hub

1 2

j z j mr

xi jmr-1m r

xi j m m

m1 zim1

mr

m2

xi j m

r-1

mr

j xi j mr j

Fig. 1. Variables at value 1 in the path for sending flow between i and j.

Finally, we define the following four indices x variables. For all i, j, m, k ∈ N, k  m let ij xkm

=

⎧ ⎨1 ⎩

0

if the flow from i to j traverses arc (k, m) connecting hubs k and m, otherwise. ij

That is, xkm = 1 means that the link connecting hubs m and k, m  k, belongs to the small tree, and also that, when sending the flow from i to j, this link is used in the direction from k to m. Note that in any optimal solution, for any pair i, j, the path that is used for sending the flow from i to j will not use any arc entering node i, nor any arc leaving node j, since such a path would not be a simple one. ij Therefore, for any i, j, all variables xkm with k = j or m = i will take value zero in any optimal solution. We only define them for ease of presentation, since, otherwise we should exclude the indices k = j or m = i throughout the paper. The x variables will be used, together with the z variables, to formulate the unique path that exists in the large tree between each pair of vertices. Before giving the formulation of the problem we workout the constraints that are necessary to model such paths. 2.1. Formulating the paths between pairs of vertices Recall that in the THLP for each pair of vertices i, j ∈ N, the flow Wij must be sent through the (unique) path of the large tree that connects i and j. Recall also that this path must contain at least one hub and that the only possible non-hub nodes in the path are i and/or j, so that the only arcs in the paths that may not belong to the small tree are the first and the last ones. Fig. 1 illustrates the four possible types of paths and the relation between the z and the x variables in each case. In particular, let i − m1 − m2 · · · − mr − j denote the path of the large tree from i to j. The x variables are used to represent the arcs of these paths that belong to the small tree, whereas some z variables are used to represent the arcs of the paths that do not belong to the small tree. When i is a hub (i.e. zii = 1), then the first arc of the path belongs to the small tree and it is represented by ij variable xim =1. When i is not a hub, the first arc in the path is (i, m1 ), 1 i  m1 which does not correspond to a link of the small tree. Such arc is represented by means of variable zim1 = 1. When the destination node j is a hub (i.e. zjj = 1), arc (mr , j) belongs to the small tree and it

3120

I. Contreras et al. / Computers & Operations Research 36 (2009) 3117 -- 3127 ij

is represented by variable xm j = 1. When j is not a hub, then the last r arc in the path is (mr , j), mr  j, which does not correspond to a link of the small tree and is represented by means of variable zjmr = 1. ij

Therefore, for each pair of nodes i, j, we use variables x·· together with variables zik , i  k, and variables zjk , j  k for formulating a path between i and j. This is done with the flow equilibrium constraints associated with the network Nij = (V, Aij ), where some nodes are connected by more than one arc (although each copy of the arc is associated with a different variable). In particular the set of arcs Aij and their associated variables are: • One copy of arc (i, m) associated with variable zim , for all m  i. • One copy of arc (m, j) associated with variable zjm , for all m  j. ij

• One copy of arc (i, m) associated with variable xim , for all m  i. ij

• One copy of arc (m, j) associated with variable xmj , for all m, m  i and m  j. ij • One copy of arc (k, m) associated with variable xkm , for all k, m, k  m, k  i, k  j, m  i, m  j. Therefore, using the node-arc incidence matrix of the network Nij , with a right hand side vector having a “+1” in the row corresponding to vertex i and a “−1” in the row of vertex j (the rest of the components being zero) the constraints that model paths between the pair i, j, i  j in the network Nij are  mi

ij xim

 mk

−

+ zji +



ij

 mj

zim = 1,

ij

 mk

xmj − zij −

mj

Therefore, when i  j the set of constraints that define the path that connect i and j in the large tree can be expressed as  mk

mi

ij

mk

xmk − zik = 0

∀k.

(1)

mi

zik − zik = 0 if k  i and k  j,  zim = −1, − which, in fact, coincide with the single assignment constraints.

ij

xmk − zik = 0



if k  i, j,

2.2. Formulation for the THLP Taking into account the above considerations, the formulation that we propose for the THLP is the following:

zjm = −1.

mj

The left hand side of the first constraint (k = i) represents the arcs outgoing node i. One such arc must be selected in the path from i to j. When i is a hub this arc will connect i with another hub m, unless j is assigned to i, in which case zji = 1. Otherwise, if i is not a hub, the first arc in the path from i to j will connect i with the hub m it is allocated to. As explained above, no optimal solution will use any arc entering node i in the path that sends the flow from i to j. Thus, in  ij this equation we can subtract the term m  i xmi or omit such term. When k  i and k  j, there will be an arc leaving node k in the path from i to j if and only if there is one arc entering node k. The candidate arcs for leaving node k are those that assume that k is a hub and that the flow from i to j goes through k and, in particular, from k to another hub m (possibly j), or from k to non-hub j, in the case that j is not a hub and it is allocated to hub k (variable zjk = 1). The candidate arcs for entering node k also assume that node k is a hub and the path from i to j goes through k. In particular, they are either arcs from another hub m (possibly i) to k or, in the case that i is not a hub but it is allocated to hub k, the arc (i, k) which is represented by the variable zik = 1. The flow equilibrium constraint associated with k = j has a symmetric interpretation to that associated with k = i. Since each node has a single assignment (later we will prove that  the node it is assigned to node is a hub), we have that m∈N zim = 1. Therefore, the flow equilibrium constraint associated with k = i can be rewritten as ij



ij

xkm + zjk −

In addition, it is well known that in flow equilibrium systems associated with node-arc incidence matrices there is one redundant constraint. Thus, in our formulation, in the block corresponding to the pair i, j, we arbitrarily eliminate the constraint associated with row k = j. When i = j, the large tree contains a (non-empty) path from i to i only if i is not a hub (zii = 0). Now, the corresponding (trivial) set of constraints is given by  zim = 1,

(F)



mj

mi

mi

xkm + zjk −

Similarly, the flow equilibrium constraint associated with k = j can be rewritten as  ij  ij xmj − zij − [1 − zjj ] = −1 ⇐⇒ − xmj + zjj − zij = 0. −

xim + zji + [1 − zii ] = 1 ⇐⇒

 mi

ij

xim + zji − zii = 0.

min

 i

s.t.



(cik Oi + cki Di )zik +

  i

k

zik = 1

j

k mk

Wij ckm xijkm

∀i,

(2)

k



zkk = p,

(3)

k

 mk ij

ij

xkm + zjk −

 mk

ij

xmk − zik = 0

ij

xkm + xmk  ykm ∀i, j, k, ∀m > k,  ykm = p − 1, k m>k ij xkm  0

∀i, j, k, i  j, k  j,

(4) (5) (6)

∀i, j, k, m, k  m,

zik ∈ {0, 1} ykm ∈ {0, 1}

∀i, k, ∀k, ∀m > k.

The second term of the objective function represents the costs of the arcs (k, m) that are used with discount (i.e. arcs connecting two hubs). Observe that the total flow through one such arc is the sum of the flows between pairs (i, j) that use arc (k, m) for sending the flow Wij . The first term of the objective function accounts for the costs of the links of the large tree that do not belong to the small tree. Each such link represents the allocation of a non-hub node i to a hub node k. Thus, the total (outgoing plus incoming) flow through the link is Oi + Di . Constraint (2) ensure that each node is allocated to one node (in Proposition 1 we prove that this node must be a hub), whereas constraint (3) set the number of hubs to p. As explained throughout Section 2.1, constraint (4) model the condition that the variables x and z define paths in the large tree between each pair of vertices.

I. Contreras et al. / Computers & Operations Research 36 (2009) 3117 -- 3127

Constraint (5) relate the x and the y variables. By definition, the only x variables at value one are associated with arcs that connect two hubs in the paths that are used for sending the flow between nodes. These arcs are the links of the small tree which are represented by the y variables. Constraint (6) is a necessary condition to guarantee that y defines a tree which on its own is not enough. As we next see, together with constraints (5) and (4) we can proof that the y variables define a spanning tree. Proposition 1. Let zˆ , yˆ , xˆ vectors of appropriate dimensions, zˆ i,k , yˆ k,m ∈ ij {0, 1} (∀i, k, m, m > k), xˆ kn  0 (∀i, j, k, m, k  m) that satisfy constraints (2)–(6). Then, 1. The set of edges such that yˆ km = 1 defines a spanning tree on the subgraph induced by the nodes such that zˆ ii = 1. 2. Non-hub nodes are allocated to hub nodes, i.e. zˆ ik  zˆ kk , ∀i, k. Proof. 1. Observe first that when a hub is located at node i (zˆ ii = 1) then, since p  2, there exists j  i such that zˆ ji = 0. Therefore, for the pair i, j the constraint (4) for k = i guarantees that there exists m ij such that xˆ im = 1. In turn, constraint (5) ensure that either yˆ im = 1 or yˆ mi = 1 (depending on whether or not i < m). That is, for each hub i there exists at least one node m such that yˆ im = 1 or yˆ mi = 1. Observe also that the yˆ variables at value one define a connected structure, and the only yˆ variables at value one are the ones associated with edges connecting two hubs. This follows since the flow equilibrium constraint (4) guarantee that xˆ variables at value one define a structure that connects all hub nodes. Since the number of hub nodes is fixed to p, the number of xˆ variables at value one repij resenting arcs between hubs (i.e. xˆ km = 1 with k  m) must be at least p − 1. Since constraint (6) sets the number of positive yˆ 's to p − 1, the number of the above xˆ 's at value one is exactly p − 1. In turn, this implies that the only yˆ 's at value one are precisely those assoij ij ciated with edges (k, m), m > k such that xˆ km + xˆ mk = 1 for some i, j. Therefore, the edges associated with yˆ 's at value one define a connected structure of p − 1 edges that connect all p hubs, thus defining a spanning tree on the graph induced by the hubs. 2. To see that zˆ ik  zˆ kk ,∀i, k, observe that when i is a non-hub node which is allocated to hub k (i.e. zˆ ik = 1) then, since p  2, there exists j s.t. zˆ jk = 0. Therefore, constraint (4) for the pair i, j and the index k,

Let  F denote the formulation F reinforced with constraints y defines a forest, zkm + ykm  zmm

zmk + ykm  zkk

(7) ∀k, ∀m > k,

(8)

∀k, ∀m > k.

(9)

Since a tree is a forest with one single component, Proposition 1 guarantees that any feasible solution to F also satisfies constraint (7). Also by Proposition 1, the flow constraint (4) guarantee that feasible solutions to formulation F satisfy conditions (8) and (9). Thus, the set of feasible solutions to F and  F coincides. However, it is easy to check that optimal solutions to the LP relaxation of formulation F need not satisfy constraints (7)–(9). Hence, the LP bound of formulation  F is at least as good as the LP bound of formulation F. Therefore, taking into account the above mentioned properties of the Lagrangean bounds, in order to reinforce the quality of our Lagrangean bound, throughout this section we work with formulation  F. Now, we formulate the Lagrangean relaxation that results when relaxing in a Lagrangean fashion the flow constraint (4), together with the new sets of constraints (8) and (9). While it is clear that the Lagrangean bound associated with this relaxation is dominated by that of the relaxation of only the flow constraint (4), it is also clear that (i) we obtain a Lagrangean function that can be decomposed in two independent subproblems, which would not be possible if we do not relax constraints (8) and (9) and (ii) the obtained Lagrangean lower bound still dominates the LP bound of formulation  F. This idea of reinforcing the Lagrangean function with valid inequalities to the original problem which are, in turn, relaxed in a Lagrangean fashion, with the only objective of improving the resulting Lagrangean relaxation, is not new (see, for instance [25]). The Lagrangean function that we obtain when relaxing constraints (4), (8), and (9) in a Lagrangean fashion with a multipliers vector (w, u, v) of appropriate dimension is L(w, u, v)= min

 i

+ +

(cik Oi +cki Di )zik +

k

 i

ji kj



⎛ wijk ⎝

  i

 mk

ij xkm

j

k mk

+ zjk −

Wij ckm xijkm

 mk

⎞ ij xmk

− zik ⎠

ukm (zkm + ykm − zmm )

k m>k

ij xˆ kh

guarantees that there exists h such that = 1. In turn, constraint (5) implies that yˆ hk = 1 or yˆ kh = 1. As we have seen, the only yˆ 's at value one are those associated with edges connecting two hubs. Thus, zˆ kk = 1. 

3121

+ s.t.





vkm (zmk + ykm − zkk )

k m>k

zik = 1

∀i,

(10)

k

 Remark. Since the z variables are required to be binary, and node-arc incidence matrices are totally unimodular, there exists some optimal solution to model F where the x variables will be integer valued. However, integrality conditions on the z and the y variables are needed to guarantee that their values are Boolean.

zkk = p,

(11)

k



ykm = p − 1, k m>k ij ij xkm + xmk  ykm ∀i, j, k, ∀m > k, y defines a forest, ij

xkm  0

3. Lagrangean relaxation for THLP Lagrangean relaxation (LR) is a well known method to solve largescale combinatorial optimization problems. It exploits the inherent structure of the problems to obtain lower bounds on the value of the optimal solution. For a recent review of LR techniques and applications, see [24]. Before describing the relaxation let us recall (see, for instance [21]) that (i) any Lagrangean relaxation of a given formulation will give a lower bound at least as good as the LP bound of the formulation and (ii) when the considered Lagrangean relaxation has the integrality property, the bound of the Lagrangean dual will coincide with the LP bound.

∀i, j, k, m, m  k,

(14) (15)

zik ∈ {0, 1} ∀i, k,

(16) (17)



cik zik +



hkm ykm k m>k     ij ij + F km xkm i j k mk i

s.t.

(13)

ykm ∈ {0, 1} ∀k, ∀m > k. After some algebra L(w, u, v) can be expressed as L(w, u, v) = min

(12)

k

(10)–(17)

3122

I. Contreras et al. / Computers & Operations Research 36 (2009) 3117 -- 3127

where the coefficients of the objective function are:   ⎧     ⎪ ⎪ uji + vij − wiji if i = k, ⎪ cii (Oi + Di ) − ⎪ ⎪ i j< i ⎨  j> i j  i • cik = (cik Oi + cki Di ) + uik + wijk − wijk if i < k, ⎪ ⎪ ji j  i,j  k ⎪   ⎪ ⎪ ⎩ (cik Oi + cki Di ) + vki + wijk − wijk if i > k, ji

• hkm = ukm + vkm

•

j  i,j  k

∀k, ∀m > k, and

⎧ Wii ckm ⎪ ⎪ ⎨ ij Wij ckm + wijk − wijm F km = ⎪ ⎪ Wij ckj + wijk ⎩ Wij cjm − wijm

if if if if

i = j, i  j, i  j, i  j,

∀k, m, k  j, m  j, m  k, k  j, m = j, m  k, k = j, m  j, m  k.

Note that L(w, u, v) can be decomposed into two independent subproblems: (1) a problem Lz (w, u, v) in the space of z variables and (2) a problem Lyx (w, u, v) in the space of y and x variables. The subproblem in the space of the z variables is Lz (w, u, v) = min

 i

cik zik

k

s.t. (10), (11), (16). The subproblem in the space of the y and x variables is Lyx (w, u, v) = min

 k m>k

hkm ykm +

  i

j

k mk

ij

ij

∀i, j. However, when ykm = 1, variable xkm will take the value one ij F km

< 0. Similarly, when ymk = 1, variable ij ij ij will take the value one whenever F mk < F km and F mk < 0. There-

whenever ij xmk

ij ij F km  F mk

and

fore, the x variables can be eliminated from subproblem Lyx (w, u, v), since the contribution to the objective function when ykm = 1 can  ij ij be expressed as hkm = hkm + i j min{0, F km , F mk }.An equivalent formulation for subproblem Lyx (w, u, v) is

Lyx (w, u, v) = min



hkm ykm

k m>k

s.t.

(12), (14), (17).

Observe that the optimal solution to subproblem Lyx (w, u, v) can be obtained by applying Kruskal's algorithm and stopping when p−1 edges have been selected. While solutions to Lyx (w, u, v) do not contain cycles, they may contain more than one connected component and, thus, it is possible that they do not define a tree. The reader will indeed notice that if we would further restrict the y variables to define a tree of cardinality p − 1, subproblem Lyx (w, u, v) would be a minimum weighted (p − 1)-tree problem which is NP-hard (see [26]). This remark justifies our choice in the reinforcement of the y variables in the formulation  F.

ij

F km xkm

3.3. The solution to the Lagrangean dual

s.t. (12), (14), (15), (17). 3.1. The solution to subproblem Lz (w, u, v) Without constraint (11), the optimal solution to Lz (w, u, v) would assign each index i to the index k with minimum value of the coefficient cik . But constraint (11) requires that exactly p indices are “self-assigned”. For selecting such indices, ∀i ∈ N, let k(i) = k if k ∈ argmin{cik : k ∈ N} and let k ∈ argmin{cik : k ∈ N} be arbitrarily selected, otherwise. Let also K = {i ∈ N : k(i) = i}, and p1 = |K|: • When p1 = p the optimal solution to Lz (w, u, v) is given by zkk = 1, k ∈ K; zik(i) = 1, i ∈ N\K; and zik = 0, otherwise. • When p1 < p, we identify the p2 = p − p1 indices in N\K with smallest increment from “best assignment” to “self-assignment”. That is ∀i ∈ N\K, i = cii − cik(i) > 0. Let i1 , i2 , . . . , ip2 denote the p2 indices in N\K with the smallest 's. Then, the optimal solution to Lz (w, u, v) is given by zkk = 1, k ∈ K ∪ {i1 , i2 , . . . , ip2 }; zik(i) = 1, i ∈/ K ∪ {i1 , i2 , . . . , ip2 }; and zik = 0, otherwise. • When p1 > p, we identify the p2 = p1 − p indices in K with smallest increment from “self-assignment” to “second-best assignment”. That is ∀i ∈ K, let m(i) ∈ argmin{cim : m ∈ N, m  i} be arbitrarily selected. Now, ∀i ∈ K, i = cim(i) − cii > 0. Let i1 , i2 , . . . , ip2 denote the p2 indices in K with the smallest 's. Then, the optimal solution to Lz (w, u, v) is given by zkk = 1, k ∈ K\{i1 , i2 , . . . , ip2 }; zim(i) = 1, k ∈ {i1 , i2 , . . . , ip2 }; zik(i) = 1, i ∈ N\K; and zik = 0, otherwise. Note that subproblem Lz (w, u, v) has the integrality property, since its optimal solution does not change if integrality constraints on the z variables are substituted by the conditions 0  zik  1, ∀i, k. 3.2. The solution to subproblem Lyx (w, u, v) The optimal solution to subproblem Lyx (w, u, v) requires to set at value one exactly p−1 variables y. If one such variable, say ykm , m > k, ij ij takes the value zero, then constraint (13) indicate that xkm = xmk = 0,

Proposition 2. For a given vector of multipliers (w, u, v) the Lagrangean function L(w, u, v) can be solved in O(n4 ). Proof. The time required to solve subproblem Lz (w, u, v) is dominated by that needed for obtaining the coefficients cik , ∀i, k. For a given pair i, k, cik can be computed in O(n). Thus all the coefficients can be obtained in O(n3 ). Finding all the coefficients k(i) and m(i) can be done in O(n2 log n), and selecting the elements to change from the best assignment requires O(p). The time required to solve subproblem  Lyx (w, u, v) is also dominated by that needed for obtaining the coefficients hkm , ∀k, ∀m > k. ij

Each F km can be obtained in O(1), so that for a given pair k, m, m > k,  ij ij the computation of hkm =hkm + i j min{0, F km , F mk } requires O(n2 ). Therefore, the calculation of all the coefficients hkm , requires O(n4 ). Once the coefficients are known, the solution to Lyx (w, u, v) can be obtained in O(n + p log p).  In order to obtain the best lower bound we solve the Lagrangean dual that is given by

(D)

zD = max w

L(w, u, v) u  0, v  0.

Given that both Lz (w, u, v) and  Lyx (w, u, v) have the integrality property, L(w, u, v) also has the integrality property and, thus zD is at least the LP bound to formulation F (zD would coincide with the LP bound if constraints (8) and (9) would not have been considered). As usual, D can be solved by subgradient optimization. A scheme of the subgradient algorithm is depicted in Algorithm 1, where the subgradient vector is partitioned in three blocks (, , ) corresponding to the relaxed blocks of constraints (4), (8), and (9), respectively. For a ˆ u, ˆ v), ˆ let z(w, ˆ u, ˆ v), ˆ y(w, ˆ u, ˆ v) ˆ and x(w, ˆ u, ˆ v) ˆ denote an given vector (w,

I. Contreras et al. / Computers & Operations Research 36 (2009) 3117 -- 3127

ˆ u, ˆ v),the ˆ optimal solution to L(w, components of a subgradient vector ˆ u, ˆ v) ˆ are given by of L(w, u, v) at (w,  ij  ij ˆ u, ˆ v) ˆ = ˆ u, ˆ v) ˆ + zjk (w, ˆ u, ˆ v) ˆ − ˆ u, ˆ v) ˆ ijk (w, x (w, x (w, mk

km

mk

mk

ˆ u, ˆ v), ˆ − zik (w, ˆ ˆ ˆ ˆ ˆ ˆ ˆ u, ˆ v) ˆ − zmm (w, ˆ u, ˆ v), ˆ mk (w, u, v) = zkm (w, u, v) + ykm (w, ˆ u, ˆ v) ˆ = zmk (w, ˆ u, ˆ v) ˆ + ykm (w, ˆ u, ˆ v) ˆ − zkk (w, ˆ u, ˆ v). ˆ mk (w, The output of the algorithm is a lower bound zD and z∗ denotes a known upper bound on the optimal value of the original problem. k The parameter  is halved after 25 consecutive iterations without improving the lower bound and is reset to the value 2 every 100 iterations. The termination criteria are: (i) when all the components of the subgradient vector are zero (optimality has been achieved); (ii) when the lower bound has not improved in the last 50 iterations; and (iii) when a limit of 2000 on the number of iterations has been reached. Algorithm 1. Subgradient method. Iteration 0 k

Initialize zD = −∞; w0 = 0;  = 2. Let z∗ be a known upper bound on the optimal value. Iteration k Step 1. Solve the Lagrangean function L(wk , uk , vk ) Step 2. if (L(wk , uk , vk ) > zD )zD ← L(wk , uk , vk ) Step 3. Evaluate the subgradient ((wk , uk , vk ), (wk , uk , vk ), (wk , uk , vk )) Step 4. Calculate the step length ∗ k ,uk ,vk )| ←  ((wk ,uk ,vk ),|z−L(w (wk ,uk ,vk ),(wk ,uk ,vk )) 2 Step5. (wk+1 , uk+1 , vk+1 ) ← (wk , uk , vk )

tk

by the coefficients hkm . For the assignment of non-hubs to hubs we proceed as follows. Let i ∈/ H, and let k(i) be such that zi,k(i) (w, u, v)=1. Then if k(i) ∈ H, we assign i to k(i). Otherwise, we assign i to  k(i) ∈ H such that ci,k(i) = min{cik : k ∈ H}. Finally, the local search explores two different neighborhoods in which the set of open hubs does not change. The node-shift neighborhood considers all solutions that can be reached from the current one by changing the allocation of exactly one non-hub node, whereas the arc-shift neighborhood considers all solutions that can be reached from the current one by changing exactly one arc in the tree of hubs. Let s = (H, a, T(H)) be the current solution, where a : N −→ H denotes the assignment of nodes to hubs, then Nnode-shift (s) = {s = (H, a , T(H)) : ∃!i ∈ N\H, a (i)  a(i)} and Narc-shift (s) = {s = (H, a, T (H)) : ∃!(i, j) ∈ T(H), ∃!(i , j ) ∈ T (H), (i, j) ∈/ T (H), (i , j ) ∈/ T(H)}.

For exploring Nnode-shift we consider all pairs of the form (i, j) with i ∈ H, j ∈ N\H and a(j)  i. In the case of Narc-shift , we only consider the subsets of arcs (i, j) ∈ T(H) and (i, j ) ∈ T (H) where i is a hub node of degree one (in the tree). In both cases we perform the first improving move. Each time the heuristic yields a feasible solution we try to update the upper bound that we use in step length calculation of the subgradient method for subsequent iterations. 5. Computational experience

k

+ tk ((wk , uk , vk ), (wk , uk , vk ), (wk , uk , vk )) Step 6. k ← k + 1

4. The heuristic for obtaining feasible solutions In order to obtain feasible solutions, and thus upper bounds that can be compared to the lower bounds generated by the Lagrangean dual, we propose a fast primal heuristic that is applied at each iteration of the subgradient optimization, and uses the information of the optimal solution to the Lagrangean function of the current iteration. Let z(w, u, v) denote the optimal solution to subproblem Lz (w, u, v), and y(w, u, v), x(w, u, v), the optimal solution to subproblem  Lyx (w, u, v). The heuristic has four main steps: • • • •

3123

Select the set of hubs. Build the tree that connects the hubs. Assign non-hub nodes to hubs. Perform a simple local search.

Two different strategies have been tested for selecting the set of hubs. The first one takes as hubs the nodes such that zkk (w, u, v) = 1. The second strategy considers the pairs k, m with m > k such that ykm (w, u, v) = 1 by increasing values of the coefficients hkm . Then, edges are considered in turn and their end-nodes k and m are selected as hubs (each of them independently) if zkk (w, u, v) = 1 or zmm (w, u, v) = 1, until p hubs are selected or until all edges such that ykm (w, u, v) = 1 have been considered. In the latter case the set of edges such that ykm (w, u, v) = 1 is again considered (again in the order given by increasing values of the coefficients hkm ) and their endnodes are selected as hubs (if they were not selected previously), until p hubs have been selected. Let H denote the set of indices of the p hubs that have been selected. Next, we build the minimum spanning tree T(H) connecting all the hubs in H, relative to the costs given

A computational study was carried out in order to test the performance of the formulation F and to compare it with formulation of Contreras et al. [7], when data for comparison are available. The formulation of Contreras et al. [7] uses three indices decision variables, Xikm , that represent the amount of flow with origin in node i that traverses the arc (k, m), where k and m are two hubs. This formulation, which we denote F3, also uses location/assignment variables (z) and design variables (y), which are exactly the same as in formulation F. In formulation F3 variables X and z are related by means of flow type constraints that have been used to formulate other discrete hub-location problems like, for instance, Ernst and Krishnamoorthy [20] and Marín [27]. We have used the AP (Australian Post) set of instances which can be downloaded from mscmga.ms.ic.ac.uk/jeb/orlib/phubinfo.html. This data set, which is commonly used in the literature, consists of the Euclidean distances between 200 cities in Australia, a code to reduce the size of the set by grouping cities, and the values of Wij (postal flow between cities). The subgradient optimization algorithm for solving the Lagrangean dual and the embedded heuristic of Section 4, have been coded in C. For comparative purposes, we have used the solver Xpress-Mosel v2.2.0, in a Windows environment, for solving model F3 of [7] and the LP relaxation of model F. All programs have been run under Windows environment on a HP mobile workstation with a processor Intel(R) Core(TM)2 T7600 with 2.33 GHz and 4 GB of RAM. We have considered instances with a number of nodes, n, ranging from 10 to 100. For ease of presentation instances have been classified in three groups: small (n = 10, 20, 25), medium (n = 40, 50, 60), and large (n = 70, 75, 90, 100). For each instance size, we combined three different values of p, namely 3, 5, and 8, and three discount factors: 0.2, 0.5, and 0.8. The results of the small instances are given in Table 1. For these instances we did the following experiments: • We used the solver Xpress-Mosel with formulation F3 of [7]. The results of this experiment are referred to as LP-F3 when they

3124

I. Contreras et al. / Computers & Operations Research 36 (2009) 3117 -- 3127

Table 1 Results with small size instances. n

10

p

3

5

8

20

3

5

8

25

3

5

8



Optimal

Upper bound

% gap

% deviation

BB-F3

LR

LP-F3

LP-F

BB-F3

Time (s) LR

LP-F3

BB-F3

LP-F

LR

0.2 0.5 0.8 0.2 0.5 0.8 0.2 0.5 0.8

52541.0 63166.8 72640.8 34340.0 49418.7 64013.2 20513.4 39288.1 57953.4

52541.0 63166.8 72640.8 34340.0 49418.7 64013.2 20513.4 39288.1 57953.4

52541.0 63166.8 72640.8 34340.0 49418.7 64013.2 20513.4 39288.1 57953.4

2.59 6.74 8.54 5.90 11.04 14.05 13.46 18.19 20.14

0.00 0.00 0.00 0.03 0.00 0.00 1.08 2.91 3.55

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 0.01 0.19 0.07 0.14 0.21 1.27 3.05 3.78

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

9.3 9.6 1.3 0.9 4.0 9.9 23.5 82.4 125.9

0.3 0.4 0.6 0.3 0.4 0.6 0.4 0.4 0.4

0.2 0.4 0.7 0.3 0.7 1.0 1.0 0.8 0.7

0.5 0.5 0.8 0.2 0.5 0.8 0.2 0.5 0.8

58761.2 69515.9 78177.6 46480.4 61061.4 73592.9 35295.1 52294.3 68272.7

58761.2 69515.9 78177.6 46480.4 61061.4 73592.9 35295.1 52294.3 68272.7

58761.2 69515.9 78177.6 46480.4 61061.4 73592.9 35383.1 52494.5 71870.8

1.91 3.57 5.32 3.60 6.86 9.08 6.46 10.26 12.57

0.00 0.00 0.00 0.00 0.00 0.29 0.05 0.13 0.43

0.00 0.00 0.00 0.00 0.00 1.70 2.89 8.23 10.96

0.00 0.13 0.21 0.03 0.38 1.56 1.08 0.81 1.55

0.5 2.1 2.1 0.7 1.4 2.1 0.4 1.2 1.9

15.6 103.5 137.4 198.1 2190.9 3600.0 3600.0 3600.0 3600.0

55.3 132.2 117.2 106.4 176.2 252.1 234.4 221.0 222.7

3.3 5.4 10.9 4.6 12.8 13.7 17.1 18.3 17.2

0.2 0.5 0.8 0.2 0.5 0.8 0.2 0.5 0.8

60602.3 70130.9 79442.4 47432.7 61046.7 73569.9 37295.6 54043.7 *70072.5

60602.3 70130.9 79442.4 47432.7 61046.7 73569.9 37295.6 54043.7 70072.5

60602.3 70130.9 79442.4 47432.7 61046.7 73710.1 37707.7 55340.2 72076.2

1.02 2.72 5.19 3.27 6.19 8.43 5.61 10.35 13.18

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.73 2.42

0.00 0.00 0.00 0.00 2.04 5.45 4.63 10.57 13.94

0.11 0.00 0.37 0.00 0.04 0.40 0.22 2.30 3.51

2.1 5.3 9.0 2.0 4.9 7.4 1.8 4.2 5.5

67.6 194.1 2550.7 661.1 3600.0 3600.0 3600.0 3600.0 3600.0

408.7 507.1 911.2 566.0 1006.6 1558.0 1276.1 2421.9 2353.0

20.3 9.1 23.8 9.7 13.8 29.8 36.8 37.8 36.5

correspond to the LP relaxation of F3 and as BB-F3, when they correspond to the branch and bound algorithm of Xpress applied to the integer formulation of the problem. • We used the solver Xpress-Mosel to solve the LP relaxation of formulation F. The results of these experiments are referred to as LP-F. • We solved the Lagrangean dual of Section 3 with subgradient optimization and performed the primal heuristic at some inner iterations. The results of these experiments are referred to as LR. For all the experiments a time limit of one hour of cpu was set. The optimal values for these instances are given in column under optimal. Instance with n = 25, p = 8 and  = 0.8, which is marked with an asterisk, could not be optimally solved within this time limit. For this instance we give the value of the best known solution. The two columns under upper bound give the value of the best solution found with the exact algorithm of Xpress applied to formulation F3 (column under BB-F3), and the best solution found with the heuristic of Section 4 applied at the inner iterations of the subgradient optimization algorithm (column under LR). The next two columns under % gap allow us to compare the values of the LP relaxations of models F3 and F. In particular, in column under LP-F3 we give the percent deviation with respect to the optimal solution of the LP relaxation of model F3 (100(Optimal − zLP-F3 )/Optimal), whereas in column under LP-F we give the percent deviation with respect to the optimal solution of the bound obtained with the LP relaxation of formulation F (100(Optimal−zLP-F )/Optimal). The next two columns combine somewhat information of the previous columns, since the entries of these two columns give, for each formulation F3 and F, the percent deviation between their best upper and lower bounds. In the case of formulation F3 when Xpress was able to optimally solve the instance, these two bounds coincide so that the percent deviation between the bounds is zero. However, for the instances that could not be optimally solved within the time limit of one hour of cpu, this

entry gives the percent deviation 100(zupper − zlower )/zlower between the best upper bound and the best lower bound at termination. In the case of formulation F, the entries of this column give the deviation 100(zheur − zLR )/zlower between the best solution found by the heuristic and the lower bound of the Lagrangean relaxation. Finally, the columns under time give the cpu times, in seconds, required by our different experiments. In particular, in column under LP − F3 we give the time needed by Xpress to solve the LP relaxation of formulation F3; in column under BB − F3 we give the time needed by Xpress to optimally solve formulation F3 (or 3600 when the time limit was reached); in column under LP − F we give the time needed by Xpress to solve the LP relaxation of formulation F; and in column LR we give the time required by the Lagrangeanrelaxation to obtain both the upper and lower bounds. The obtained results allow us to make several observations. On the one hand, we can appreciate the quality of the lower bounds obtained with formulation F. Observe that the percent gaps between the lower bound and the optimal solution are always very small (only in four out of the 27 instances of the set this gap is over 1%), and this percent gap is indeed much smaller than that of formulation F3 (although the cpu times required to solve the Lagrangean dual are always larger than those needed to solve the LP relaxation of formulation F3). On the other hand, the results under column % deviation indicate that, despite its simplicity, the heuristic gives solutions that have small deviations from the lower bounds (and, thus, from the optimal solutions). However, it is also clear that for the smaller instances of 10 nodes and for the instances with 20 nodes and small values of p, formulation F3 can be preferred, since it allows to optimally solve the instances which, in general, is not possible with the Lagrangean dual. However, as the sizes of the instances increase (or, especially, as p increases) the advantage of formulation F3 disappears, since it is no longer possible to optimally solve the instances in one hour of cpu time, and then the percent deviation between the upper and lower bounds is

I. Contreras et al. / Computers & Operations Research 36 (2009) 3117 -- 3127

3125

Table 2 Results with medium size instances. n

40

p

3

5

8

50

3

5

8

60

3

5

8



Best

Upper bound

% deviation

Time (s)

BB-F3

LR

LP-F3

BB-F3

LR

LP-F3

BB-F3

LR

0.2 0.5 0.8 0.2 0.5 0.8 0.2 0.5 0.8

62543.7 72383.2 80724.8 52603.0 65289.3 76965.1 43766.7 59455.8 74488.5

62629.6 73512.7 82624.3 53855.8 67804.6 83239.4 48693.0 61927.3 77736.0

62543.7 72383.2 80724.8 52603.0 65289.3 76965.1 43766.7 59455.8 74488.5

1.80 3.78 4.95 4.50 6.86 9.20 7.27 10.89 13.90

1.46 4.31 6.49 6.16 9.87 15.55 16.12 13.64 16.83

1.44 0.00 0.16 1.07 0.47 2.02 2.45 3.06 5.71

35.9 166.0 248.3 91.1 207.7 370.9 99.0 198.1 278.2

3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0

86.1 55.5 81.9 125.4 132.7 131.9 179.0 145.0 153.7

0.2 0.5 0.8 0.2 0.5 0.8 0.2 0.5 0.8

62504.2 72891.2 80719.8 51799.1 65423.6 76714.2 43805.6 59794.8 74393.2

62540.6 77979.4 85022.5 52358.4 75013.5 89688.7 51470.2 70948.2 n.a.

62504.2 72891.2 80719.8 51799.1 65423.6 76714.2 43805.6 59794.8 74393.2

1.95 4.31 4.95 3.69 7.10 8.92 6.58 11.14 13.38

1.80 10.44 9.55 4.39 18.62 21.78 20.07 24.74 n.a.

1.88 0.29 0.31 1.13 0.85 1.76 2.49 4.05 6.04

142.0 1390.8 1820.5 245.8 1512.8 2971.6 360.0 861.2 2166.1

3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0

211.2 221.5 266.2 253.9 327.1 293.8 334.5 343.8 326.9

0.2 0.5 0.8 0.2 0.5 0.8 0.2 0.5 0.8

62934.4 73411.3 81528.8 51438.2 65508.0 77805.8 43994.2 59746.2 74087.8

n.a. n.a. n.a. n.a. n.a. n.a. n.a. n.a. n.a.

62934.4 73411.3 81528.8 51438.2 65508.0 77805.8 43994.2 59746.2 74087.8

1.75 n.a. n.a. 3.37 n.a. n.a. 7.43 n.a. n.a.

n.a. n.a. n.a. n.a. n.a. n.a. n.a. n.a. n.a.

2.42 0.49 0.39 1.45 0.88 2.61 4.20 4.01 5.72

1644.4 3600.0 3600.0 2242.1 3600.0 3600.0 2099.5 3600.0 3600.0

3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0 3600.0

420.8 351.0 571.2 522.3 580.3 582.1 571.5 642.8 612.6

considerably larger than that of the Lagrangean relaxation. Note that the cpu time needed to obtain the lower bound with the Lagrangean dual is quite small (it never reaches 40 seconds of cpu time) and, excepting for the very small instances of 10 nodes, it is considerably smaller than that required to solve the LP relaxation of formulation F. This advantage of the Lagrangean dual relative to the LP relaxation of formulation F becomes larger as the sizes of the instances increase. The results of medium size instances are given in Table 2. For these instances it was no longer possible to solve the LP relaxation of formulation F within the time limit of one hour. Therefore, we used the solver Xpress-Mosel with formulation F3 of [7], both to solve the LP relaxation and to obtain the best upper and lower bounds within one hour of cpu time. Now, column under best gives the value of the best-known solution, and the next two columns under upper bound give the values of the best solutions found with the exact algorithm of Xpress applied to formulation F3 (column under BB-F3), and the best solution found with the heuristic of Section 4 applied at the inner iterations of the subgradient optimization algorithm (column under LR). The next three columns under % deviation give the percent deviation between the best upper and (i) the lower bound of the LP relaxation of formulation F3 (column under LP − F3), (ii) the branch and bound applied to formulation F3 (column under BB−F3), and (iii) the Lagrangean relaxation of Section 3 (column under LR). Finally, the three columns under time give the cpu times in seconds required by each experiment. As before, column under LP-F3 gives the time needed by Xpress to solve the LP relaxation of formulation F3; column under BB-F3 gives the time needed by Xpress with formulation F3; and in column LR we give the time required for solving the Lagrangean dual, including the time needed by the heuristic at the inner iterations of the subgradient optimization. As can be seen, in all cases it was not possible to solve optimally formulation F3 and the time limit of 3600 seconds of cpu time was reached. Moreover, for most of the 60 nodes instances (six instances) it was not possible to solve the LP relaxation of F3 within the time

limit, so no valid lower bound (nor upper bound) was obtained. For three more instances (one with n = 50, p = 8,  = 0.8; another one with n = 60, p = 3,  = 0.2; and another one with n = 60, p = 5,  = 0.2) Xpress could solve the LP relaxation in less than one hour but no upper bound is available, since the remaining time was not enough so as to find a feasible solution. For the remaining 17 instances, the percent deviation between the upper bound obtained with the simple heuristic was always better than that obtained with formulation F3, and the lower bound obtained with the Lagrangean relaxation was also always better than the best lower bound obtained with F3. As a consequence, the percent deviations of the Lagrangean relaxation of formulation F are always considerably smaller than those of formulation F3. It is worth noting that the percent deviation between the upper and lower bounds of the Lagrangean relaxation never exceeded 6.04%. Letting aside the more difficult instances with p = 8, the largest percent deviation reduces to 2.61%. It is also worth noting that the times required by the Lagrangean relaxation are still very small for instances of these sizes. This becomes even more evident if one takes into account the times required by Xpress with formulation F3. We conclude this section with the results of our computational experiments with the large instances. Now we report on the results obtained with the Lagrangean relaxation of Section 3, which are depicted in Table 3. Column under best gives the value of the best solution found with the heuristic of Section 4 applied at the inner iterations of the subgradient optimization algorithm. Column under lower bound gives the value of the lower bound obtained with the Lagrangean relaxation. Column under % deviation gives the percent deviation between the upper bound best and the lower bound lower bound. Finally, column under time gives the cpu times in seconds required by the Lagrangean relaxation to obtain both the lower and upper bounds. For these experiments we have released the time limit of one hour of cpu and we have increased the maximum number of iterations in the termination criterion of the subgradient optimization algorithm to 3000 (instead of 2000).

3126

I. Contreras et al. / Computers & Operations Research 36 (2009) 3117 -- 3127

Table 3 Results with large instances.

70

3

5

8

75

3

5

8

90

3

5

8

100

3

5

8

Best LR

Lower bound LR

% deviation LR

Time (s) LR

0.2 0.5 0.8 0.2 0.5 0.8 0.2 0.5 0.8

63341.1 73497.7 81681.8 52978.1 66407.8 78083.3 44922.6 61134.6 75723.6

61239.2 73281.2 81372.6 51289.5 65525.6 76074.4 42782.3 57907.2 69947.8

3.31 0.29 0.37 3.18 1.32 2.57 4.76 5.27 7.62

518.9 692.7 518.9 1070.6 1075.2 938.0 1044.4 1076.2 1117.4

0.2 0.5 0.8 0.2 0.5 0.8 0.2 0.5 0.8

63412.2 73549.6 81726.1 57725.1 66152.8 78271.4 45521.0 61885.4 76641.3

61330.2 72941.5 80697.1 56904.9 65520.6 75945.7 42802.7 58191.5 69890.6

3.28 0.46 0.65 1.42 0.96 2.97 5.97 5.95 8.80

884.0 1562.5 1400.4 1667.3 1668.4 1475.6 1722.1 1637.2 1616.6

0.2 0.5 0.8 0.2 0.5 0.8 0.2 0.5 0.8

63270.8 73259.8 81404.0 53021.8 66170.4 78146.9 45514.0 62164.7 76000.0

60795.6 72941.5 80697.1 50994.3 65223.1 75473.2 41972.2 57542.7 68859.0

3.91 0.43 0.86 3.82 1.43 3.42 7.78 7.43 9.39

1805.5 2266.2 2752.4 3245.5 2835.7 2776.8 3304.5 3289.1 3138.4

0.2 0.5 0.8 0.2 0.5 0.8 0.2 0.5 0.8

63442.5 73415.9 81473.2 53548.7 67064.8 78376.2 45816.0 62942.8 76412.4

61087.6 72789.0 80669.7 51005.3 65302.5 75289.0 41875.8 57669.9 68882.9

3.71 0.85 0.98 4.75 2.62 3.93 8.60 8.37 9.85

5072.6 5142.5 5157.0 5215.6 5173.6 5007.3 5231.9 4968.4 5230.5

In our opinion the obtained results are very good. Despite the difficulty and the large sizes of the instances, we are able to obtain very good lower bounds and very good feasible solutions with reasonable computation times. Again this becomes more evident if we take into account the difficulties of Xpress for solving the LP relaxation of Formulation F3 which, as we have observed in Table 2, gives much weaker lower bounds than the Lagrangean relaxation of formulation F and requires higher computation times. The percent deviations that we obtain are still very small, since they never exceed 10%. In fact, if we exclude the difficult instances with p = 8 the percent deviation is always below 4%. For most of the instances, the limit of 3000 iterations of the subgradient algorithm has been reached now, meaning that the subgradient optimization algorithm stopped before it converged to the optimal solution to the Lagrangean dual. Therefore, it is most likely that if more iterations were allowed, the subgradient optimization algorithm would yield a better lower bound and, thus, smaller percent deviations. 6. Conclusions and further research In this paper we have proposed a four index formulation for the THLP that empirically yields tighter LP bounds than those of previously proposed formulations. However, the computational burden required to obtain these bounds with a commercial solver becomes unaffordable as the sizes of the instances increase. For reducing the time needed for obtaining these tight bounds, we have proposed a Lagrangean relaxation that can be decomposed into independent subproblems which can be solved quite efficiently. We also obtain

tight upper bounds with a simple heuristic that is applied at the inner iterations of the subgradient optimization that is used for solving the Lagrangean dual. While the instances previously considered in the literature for the THLP did not exceed 20 nodes, the proposed algorithm has enabled us to handle efficiently instances up to 100 nodes. We have run a series of computational experiments on a set of 90 benchmark instances of sizes ranging from 10 to 100. In most cases the percent deviations between our upper and lower bounds do not exceed 4%, and they never exceed 10%. In our opinion these results are remarkable taking into account the required computation times and the difficulty of the considered problem. Although the obtained results are very good, in general, the optimality of the solutions could not be guaranteed. Thus, a natural extension of this work is to consider exact algorithms for solving THLP. Acknowledgments The authors thank Plan Nacional de Investigación Científica, Desarrollo e Innovación Tecnológica (I + D + I)—projects MTM200614961-C05-(01,04)—, Fundación Séneca—project 02911/PI/05—, and ERDF funds for the support given to the research which led to this article. The research of the first author has been partially supported through Grant 197243/217997 from CONACYT, México. References [1] O'Kelly ME, Miller HJ. The hub network design problem: a review and synthesis. Journal of Transport Geography 1994;2(1):31–40. [2] Klincewicz JG. Hub location in backbone/tributary network design: a review. Location Science 1998;6:307–35. ¨ [3] Nickel S, Schobel A, Sonneborn T. Hub location problems in urban traffic networks. In: Niittymaki J, Pursula M, editors. Mathematics methods and optimization in transportation systems. Dordrecht: Kluwer Academic Publishers; 2001. [4] Gelareh S. Hub location models in public transport planning. PhD thesis, University of Saarlandes, Germany; 2008. [5] Campbell JF, Ernst AT, Krishnamoorthy M. Hub arc location problems: part I—introduction and results. Management Science 2005;51(10):1540–55. [6] Campbell JF, Ernst AT, Krishnamoorthy M. Hub arc location problems: part II—formulations and optimal algorithms. Management Science 2005;51(10):1556–71. [7] Contreras I, Fernández E, Marín A. The tree of hubs location problem. European Journal of Operational Research, 2007, submitted for publication. Available at http://www-eio.upc.es/∼elena/. [8] Hu TC. Optimum communication spanning trees. SIAM Journal on Computing 1974;3(3):188–95. [9] Nguyen VH, Knippel A. On tree star network design. In: INOC 2007: international network optimization conference, Spa, Belgium; 2007. [10] Chen H, Campbell AM, Thomas BW. Network design for time-constrained delivery. Naval Research Logistics 2008;55:493–515. [11] Wu BY. Approximation algorithms for the optimal p-source communication spanning tree. Discrete Applied Mathematics 2004;143:31–42. [12] Hakimi SL, Schmeichel EF, Labbé M. On locating path or tree shaped facilities on networks. Networks 1993;23:543–55. [13] Kim TU, Lowe TJ, Tamir A, Ward JE. On the location of a tree-shaped facility. Networks 1996;28:167–75. [14] Puerto J, Tamir A. Locating tree-shaped facilities using the ordered median objective. Mathematical Programming 2005;102:313–38. [15] Marín A. An extension to rapid transit network design problems. TOP 2007;15(2):231–241. [16] Labbé M, Laporte G, Rodríguez-Martín I, Salazar-González JJ. The ring star problem: polyhedral analysis and exact algorithm. Networks 2004;43(3):177– 89. [17] Laporte G, Rodríguez-Martín I. Locating a cycle in a transportation or a telecommunications network. Networks 2007;50(1):92–108. [18] Campbell JF. Integer programming formulations of discrete hub location problems. European Journal of Operational Research 1994;72:387–405. [19] Skorin-Kapov D, Skorin-Kapov J, O'Kelly M. Tight linear programming relaxations of uncapacitated p-hub median problems. European Journal of Operational Research 1996;94:582–93. [20] Ernst AT, Krishnamoorthy M. Solution algorithms for the capacitated single allocation hub location problem. Annals of Operations Research 1999;86:141– 59. [21] Geoffrion AM. Lagrangean relaxation for integer programming. Mathematical Programming Study 1974;2:82–114.

I. Contreras et al. / Computers & Operations Research 36 (2009) 3117 -- 3127

[22] Contreras I, Díaz J, Fernández E. Lagrangean relaxation for the capacitated hub location problem with single assignment. Operations Research Spectrum, 2009, doi:10.1007/500291-008-0159-y. [23] Contreras I, Díaz J, Fernández E. Branch and Price for Large-Scale capacitated hub location problems with single assignment. INFORMS Journal on Computing, 2008, submitted for publication. [24] Guignard M. Lagrangean relaxation. TOP 2003;11:151–228. ¨ Fernández E, Jornsten ¨ [25] Barceló J, Hallefjord A, K. Lagrangean relaxation and constraint generation procedures for capacitated plant location problems. Operations Research Spectrum 1990;12:79–88.

3127

[26] Garey MR, Johnson DS. Computers and intractability: a guide to the theory of Np-completeness. New York: Freeman; 1979. [27] Marín Á. Formulating and solving splittable capacitated multiple allocation hub location problems. Computers & Operations Research 2005;32: 3093–109.