Minimizing the maximum travel time in a combined model of facility location and network design

Omega 40 (2012) 847–860

Contents lists available at SciVerse ScienceDirect

Omega journal homepage: www.elsevier.com/locate/omega

Minimizing the maximum travel time in a combined model of facility location and network design Ivan Contreras a,n, Elena Ferna´ndez b, Gerhard Reinelt c a

Concordia University and Interuniversity Research Centre on Enterprise Networks, Logistics and Transportation (CIRRELT), Montreal, Canada Statistics and Operations Research Department, Technical University of Catalonia, Barcelona, Spain c Institute of Computer Science, University of Heidelberg, Heidelberg, Germany b

a r t i c l e i n f o

abstract

Article history: Received 27 January 2011 Accepted 23 January 2012 Processed by Associate Editor Pesch. Available online 24 February 2012

This paper presents a combined Facility Location/Network Design Problem which simultaneously considers the location of facilities and the design of its underlying network so as to minimize the maximum customer-facility travel time. The model generalizes the classical p-center problem and has various applications in regional planning, distribution, telecommunications, emergency systems, and other areas. Two mixed integer programming formulations are presented and compared. Several valid inequalities are derived for the most promising of these formulations to strengthen its LP relaxation bound and to reduce the enumeration tree. Numerical results of a series of computational experiments for instances with up to 100 nodes and 500 candidate links are reported. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Facility location Network design Integer programming Valid inequalities

1. Introduction Network location models cover a wide range of applications within the design of transportation and telecommunication systems. These models are frequently used to analyze and determine the location of public and private facilities within a prespecified network. There exist a variety of network location models which can be categorized according to the type of objective that is considered. The min-sum models, such as the classical p-median problem [17] and the uncapacitated facility location problem ([19]), focus on minimizing the overall set-up or operational cost. The min–max models, such as the p-center problem [17], focus on minimizing the largest customer-facility distance. The covering models, such as the set covering location problem [32] and the maximum covering location problem [4], focus on finding the least number of facilities to cover all customers or the maximum number of customers to cover within a prespecified distance, respectively. All these location models assume the underlying network to be given as an input. However, in some contexts the topology of the network plays a crucial role for the optimal location of facilities (see [11]). Economically it may be more attractive to modify the underlying network than to add facilities. For this reason, combined facility location network design models play an important n

Corresponding author. Tel.: þ1 514 8484 2424x3130. E-mail addresses: [email protected], [email protected] (I. Contreras), [email protected] (E. Ferna´ndez), [email protected] (G. Reinelt). 0305-0483/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.omega.2012.01.006

role in situations in which there exists a non-trivial tradeoff between set-up costs for the location of facilities and links of the network as well as operating costs for providing service to customers. These models have a broad range of applications related to the design of distribution systems, air and ground transportation networks, telecommunication networks, hub-andspoke networks, and regional planning efforts (see [23] for details). Facility Location/Network Design Problems (FLNDPs) are a challenging class of combinatorial optimization problems that combine two types of decisions. The location decisions consist of selecting a set of nodes to locate facilities and designing the allocation pattern of nodes to facilities. The network design decisions consist of selecting a set of links to provide the connection of nodes to their allocated facilities. Several costs may affect the structure of optimal solutions such as set-up costs of facilities, set-up costs of links to connect the nodes to facilities, and operating (or service) costs to transport the demand through the network. FLNDPs were initially introduced by Daskin et al. [11] as generalizations of classical facility location problems. The Uncapacitated Facility Location/Network Design Problem (UFLNDP) was proposed by Daskin et al. [11] as a generalization of the uncapacitated facility location problem and further studied by Melkote and Daskin [23]. The UFLNDP considers the minimization of the overall operational costs as well as the set-up costs of facilities and links. As the magnitudes of set-up and service costs in UFLNDP may be non-homogeneous, adding the above three terms into one single objective may be meaningless. One possibility for considering all these costs simultaneously is to minimize

848

I. Contreras et al. / Omega 40 (2012) 847–860

the operating costs while requiring that the overall set-up costs satisfy a given budget constraint. This approach is often applied in network design problems (see, e.g. [15,33]). Adding a budget constraint to FLNDP gives rise to the so-called Facility Location/ Network Design Problem with Budget constraint (FLNDB), which was introduced by Melkote and Daskin [23] and further studied by Cocking et al. [5] and Cocking [6]. Both UFLNDP and FLNDP assume that the amount of demand that can be served by facilities is unbounded. This is valid in situations in which it is known a priori that the facilities will operate significantly below their capacity levels. Nevertheless, there exist situations in which the capacities of the facilities are very constraining. Melkote and Daskin [24] presented the Capacitated Facility Location/Network Design Problem (CFLNDP), which extends the UFLNDP to consider a limited amount of capacity at the facilities. Drezner and Wesolowsky [13] deal with another class of min-sum FLNDPs in which the links can either be one-way or two-way and the service costs consider round trips from the facility to each node and back. Like in many network location models, the min-sum objective is typically useful in cases in which the objective is to minimize the average service cost of the network without taking into consideration if some customers have not been properly allocated in terms of the distance to the facilities. However, there exist contexts in which a cost-based objective is not an appropriate measure of the quality of the solution. This is the case when designing networks for emergency services such as the location of health facilities or emergency units. Service costs represent customer-facility travel times, which are critical for each individual customer and should not be overlooked using a cost-based objective. In this paper we introduce the Center Facility Location/Network Design Problem with Budget constraint (CFLNDB), which consist of selecting a set of nodes to locate facilities and a set of links to allow the connection between customers and facilities. The CFLNDB considers both service costs and design (or set-up) cost. Service costs are taken into account in the objective function which aims at minimizing the maximum travel time between any node to their allocated facilities, rather than minimizing the overall service costs as in the previously proposed FLNDPs. The overall design costs are taken into account by means of a budget constraint which limits the overall set-up costs to a maximum budget. It is assumed that the amount of demand that can be served from the facilities is unbounded. In order to build a more flexible network, we let the model decide the best number of facilities and links to install within a prespecified budget, rather than specifying a priori the number of them. The CFLNDB is a generalization of the classical p-center problem, which has been extensively studied (see [27,28]). However, to the best of the authors’ knowledge the CFLPNDB has not been studied so far. The CFLNDB involves one additional difficulty with respect to classical center location problems, because for tracing the paths that determine the value of the objective function, the location of the facilities and the allocation of nodes to facilities must be known, in addition to the selected arcs. Potential transportation applications of FLNDPs where a min– max approach is required arise in regional planning and land reuse programs. In this case, the government may be simultaneously considering the construction of new roads as well as the location of public facilities such as hospitals, police and fire stations, and schools. Given that these facilities provide essential and emergency services, it is strongly desirable to minimize the maximum travel time from any potential customer to its closest facility. A concrete example of an application of the CFLNDB to a regional planning program appears in the design of a road structure for improving the travel times to health facilities in the Nouna District of Burkina Faso, studied in Cocking et al. [5]

and Cocking [6]. In Nouna, travel times of population to health services become extremely high during the rainy season, since many roads are unusable. Cocking [6] approaches this problem by deciding the number and location of the health facilities over a potential set of villages and the roads to be built or improve to connect the health facilities with a set of nodes representing the population. A min-sum objective is considered and the problem is formulated as a particular case of the FLNDB. However, in this particular application a min–max objective function is more appropriate for the reasons already described. Other papers, such as Syam and Cˆote´ [31] and Murawski and Church [26], have also studied the location of health services, however, they have focused on other types of objectives such as min-sum or maxcovering. Additional applications of FLNDPs where a min–max approach is required arise in the design of less-than-truckload (LTL) distribution systems. In this case, transportation times through the network are critical in achieving delivery deadlines. Therefore, LTL carrier companies are interested in designing a transportation network such that it minimizes the maximum travel time from the breakbulk facilities to the distribution points, subject to budgetary limitations on the design costs for the location of breakbulks and the establishment of a service link between two points. 1.1. Related literature Combined problems of facility location and network design have been studied under different names, according to their particular applications. Contreras and Ferna´ndez [7] presents a unified framework for the general network design problem which encompasses several classical problems involving combined location and network design decisions. Recent surveys on discrete facility location such as ReVelle and Eiselt [27], Klose and Drexl [20], ReVelle et al. [28], and Smith et al. [29], point out some location models in which there exist network design decisions. Melo et al. [25] presents a survey containing more than 120 references of facility location models arising in the design of supply chain networks. The reviewed models include a combination of strategic features such as multilayer facilities, multiple commodities, single/multiple period(s), and deterministic/stochastic parameters. Some of these models also include tactical and operational supply chain decisions such as capacity design, inventory control, procurement, production planning and routing. Srivastava [30] introduces an integrated framework for the location-allocation of facilities to design cost-effective reverse logistic networks. Cordeau et al. [9] propose an integrated model for logistic network design involving strategic decisions regarding the number, location, capacity and technology of manufacturing plants and warehouses. It also considers tactical decisions such as the selection of suppliers, product range assignment, distribution channel, and transportation modes. Some operational decisions, like the routing of raw materials, semi-finished and finished products through the network, are also consider in the model. The literature on problems involving network decisions with min–max objective is scarce. Campbell et al. [3] present a model to minimize the maximum travel time in specific network design problems that involve no location decisions. They consider a series of problems that involve the selection of q arcs in already established networks so as to minimize the total diameter of the upgraded network. In [1], several hub location problems with min–max objectives are described. However, these problems consider the requirement that hub facilities induce a complete graph, making the evaluation of the objective function much easier. We are not aware of any work that considers a min–max objective with joint location and network design decisions

I. Contreras et al. / Omega 40 (2012) 847–860

without making additional requirements on the topology induced by the selected arcs. Recently, some papers have appeared dealing with combined models of facility location and network design. Murawski and Church [26] study a FLNDP in which the objective is to maximize the demand that is captured subject to a budget constraint on the total installation cost of the selected arcs. Costa et al. [10] present a network design model with location ingredients arising in electrical distribution systems. The problem consists of selecting links to locate two types of edge facilities as well as locating node facilities to allow the connection of edge facilities. The objective is to minimize the sum of design costs of edge and node facilities and the operational costs. Gollowitzer and Ljubic´ [16] studies connected facility location problems, which combine facility location and Steiner trees. These problems consider the location of facilities, the assignment of customers and the connection of open facilities via a Steiner tree. The objective is to minimize the design cost for the facilities and the links of the Steiner tree as well as the service costs. The rest of the paper is organized as follows. In the following section, the CFLNDB is formally defined and its complexity established. In Section 3, we present two different mixed integer programming formulations. One of them is based on multicommodity type variables to trace the paths in the network and the other is based on the fact that optimal solutions to CFLNDB are associated with rooted forests. Given that the former formulation turns out to be the most promising when used in a general purpose solver, in Section 4 we show how it can be successively improved by means of different types of valid inequalities. The description of the computational results and the analysis of the obtained results are given in Section 5. Before presenting the obtained numerical results, we analyze the effect of employing tight upper bounds to reinforce some constraints of the formulation, and for simple but effective reduction tests. Section 6 draws some conclusions and some comments on possible lines for future research.

849

node. As a consequence, any arc not belonging to any such path can be removed from any feasible solution without affecting its feasibility or its objective function value. Moreover, given that it is assumed that nodes where facilities are located are self-assigned and there are no capacity constraints, no path Pi will contain more than one facility. Therefore, arcs between two facilities will never be part of an optimal solution. This implies that there exists an optimal solution to CFLNDB that does not contain any arc connecting two facilities or a non-facility node and a facility different from the one it is allocated to. Thus, each component of such a solution contains one single facility which is its root node. In addition, if one such component contains a cycle, say C, then there exists one arc a A C not belonging to any of the paths Pi, i A V. Therefore, arc a can also be removed without affecting the feasibility or deteriorating the objective function value of the resulting solution. As a consequence there exists an optimal solution which is a rooted forest. Let R denote the set of rooted forests of D, and RB the subset of rooted forests that satisfy the budget constraint, i.e., 8 9 < = X X X RB ¼ R A R : f rk þ cij r B : : ; T AR T AR k

k

ði,jÞ A Ak

For a rooted forest R we write R ¼ fT 1 ,T 2 , . . . ,T r g, where – T k ¼ ðV k ,Ak ,r k Þ is a rooted tree with V k D V, Ak D A, and root node r k A V k , for k ¼ 1, . . . ,r, – fV 1 ,V 2 , . . . ,V k g defines a partition of V. For RA R and iA V given, we call the root of i the root node of the unique Vk such that i A V k , and we denote by tði,RÞ the travel time from node i A V k to its root node rk, which is the total travel time of the unique path Pi from i to rk in R. Then, the objective function value of R A R is given by the maximum travel time from any node iA V to its root node max tði,RÞ: iAV

Therefore, using the above notation, CFLNDB can be stated as, 2. The Center Facility Location/Network Design Problem Let D ¼ ðV,AÞ be a directed graph with node set V ¼ f1; 2, . . . ,ng and arc set A. Let f i Z0 denote the set-up (or design) cost of selecting node i as a facility. For each arc (i,j), let cij Z 0 and t ij Z 0 denote its construction cost and travel time, respectively. Furthermore, B Z0 gives the budget on the overall installation costs (facilities plus arcs). Feasible solutions to the Center Facility Location/Network Design Problem with Budget constraint (CFLNDB) consist of: – a set of nodes of V to locate facilities, – an allocation of non-facility nodes to facilities, – a subset of arcs that connect each node to its allocated facility, such that the total set-up cost of the facilities plus the construction cost of the selected arcs does not exceed the budget B. The value of a feasible solution to CFLNDB is given by the maximum travel time from any node to its allocated facility, in the graph induced by the selected arcs. Then the problem is to find a feasible solution of minimum cost. In the case of FLND and FLNDB, when all set-up and transportation costs are non-negative, it is known that there exists an optimal solution which is a rooted forest (see [23]). As we will see next, this property also holds for CFLNDB. Observe first that the only arcs that are required in a feasible solution are the ones that define the paths Pi from each node iA V to its allocated facility. For i A V the length of Pi determines the travel time from i to its root

ðCFLNDBÞ

min max tði,RÞ:

R A RB i A V

Observe that the well known p-center problem defined on a complete graph, with assignment costs dij, i,j A V satisfying the triangle inequality, is a particular case of CFLNBD with design costs fi ¼1, i A V, construction costs cij ¼ 0, i,j A V, travel times t ij ¼ dij , i,j A V, and budget B ¼p. Since the p-center problem is NP-hard (see, e.g. [18,22]), we can also establish the NP-hardness of CFLNDB. Note that without the budget constraint the node optimality property can be proved for CFLNBD using similar arguments as for the center location problem [17]. That is, even if facilities can be located at the edges of the network, an optimal solution exists in which all the facilities are located at the nodes of the network. However, like in other center location problems with additional knapsack type constraints, after introducing the budget constraint, the node optimality property no longer holds. In addition, in order to allow facilities to be located at the edges of the network, the facility set-up cost function should be defined for all the points of the edges. In any case, like it is usual in most location problems on networks, we assume in the following that facilities can only be located at the nodes of the network. 3. Mathematical programming formulations of the CFLNBD In classical location problems service costs, representing the assignment costs of customers to facilities, are known in advance.

850

I. Contreras et al. / Omega 40 (2012) 847–860

This makes it possible to model the center objective as a telescopic sum [14] or by means of covering type constraints, which limit the maximum assignment cost of a customer to its allocated facility [12]. In contrast, in bottleneck network design problems, the above techniques become extremely involved, since service costs are not known in advance. These now represent communication costs between pairs of customers, which depend on the selected arcs, and the objective is to minimize the maximum communication cost given by the diameter of the network induced by the selected arcs (see, e.g. [3]). Moreover, in order to evaluate the objective it is not enough to know which are the selected arcs, as the paths connecting the different pairs of nodes must be traced. In CFLNDP service costs represent access times of nodes to their allocated facilities in the network induced by the selected arcs, and the objective is to minimize the maximum service cost. Therefore, these problems involve one additional difficulty, because for tracing the paths that determine the value of the objective function, the location of the facilities and the allocation of nodes to facilities must be known, in addition to the selected arcs. In this section we present two different formulations for CFLNBD. The first one makes use of multicommodity type decision variables, like it is common in other network design problems (see for instance [7,21]). In the second formulation, which exploits the structure induced by optimal solutions, decision variables are defined to represent feasible rooted trees R A RB . 3.1. Formulation with multicommodity type decision variables We define several sets of decision variables to represent feasible solutions to the CFLNBD. The location of the facilities is represented by the binary variables ( 1 if a facility is located at node i, zi ¼ 0 otherwise, for each iA V. The assignment of nodes to facilities is modeled with binary variables x, where for each i,j A V with j a i, ( 1 if node i is assigned to facility j, xij ¼ 0 otherwise: The arcs selected to connect the nodes with their allocated facilities are modeled with binary variables y, where for each ðk,mÞ A A, ( 1 if arc ðk,mÞ is constructed, ykm ¼ 0 otherwise: The paths from nodes to their allocated facilities, which need to be traced in order to evaluate the objective function, are modeled by means of multicommodity type binary variables X, where for each i,j A V with j a i, and ðk,mÞ A A, we set ( 1 if arc ðk,mÞ is on the path from i to j, ij X km ¼ 0 otherwise: Finally, we define a continuous variable T that gives the maximal accumulated travel time for all paths. Then, CFLNDB can now be written as the (mixed) integer program min T, X ðk,mÞ A A

t km X ijkm r T

ð1Þ 8i,j A V,

ð2Þ

X

xij þzi ¼ 1

8i A V,

ð3Þ

j A V\fig

X

X ijik ¼ xij

8i,j A V,j a i,

ð4Þ

ði,kÞ A A

X

X ijkm 

ðk,mÞ A A

X

X ðm,kÞ A A,k a i

X ijkj ¼ xij

X ijmk ¼ 0

8i,j,mA V,j a i,m,

8i,j A V,j a i,

ð5Þ

ð6Þ

ðk,jÞ A A

xij r zj

8i,j A V,j a i,

X ijkm r ykm

8i,j, A V,j a i,ðk,mÞ A A,

X X f i zi þ ckm ykm r B, iAV

ð7Þ ð8Þ ð9Þ

ðk,mÞ A A

X ijkm A f0; 1g

8i,j, A V,j ai 8ðk,mÞ A A,

ð10Þ

xij ,zi A f0; 1g

8i,j, A V,

ð11Þ

ykm A f0; 1g

8ðk,mÞ A A:

ð12Þ

Constraints (2) keep track of the longest accumulated travel time among all possible paths. The set of equalities (3) states that, if node i is not a facility, then it has to be assigned to one. Constraints (4)–(6) are a particular case of flow conservation constraints, by which a path connecting nodes i and j is constructed, only if node i is assigned to a facility located at node j. As usual, in this type of constraints, for each pair i,j A V one equality is redundant and can be eliminated. Inequalities (7) ensure that nodes are assigned to open facilities, whereas constraints (8) guarantee that all the arcs of the paths are constructed. The budget limitation on the total set-up costs for the location of facilities and the construction of links is imposed in (9). Finally, ((10)–(12) are non-negativity and integrality conditions on the decision variables. In the following formulation (1)–(12) above is referred to as formulation F.

3.2. Alternative formulation Formulation F uses multicommodity type decision variables, which are commonly used in network design models (see, e.g. [2]). This type of variables is indeed very intuitive but, on the other hand, requires four indices, which most often limits the size of the instances that can be solved efficiently. For this reason we have developed an alternative formulation for CFLNDB which uses fewer variables and exploits the rooted tree structure of optimal solutions. We use the same variables z and y as in F and define some new decision variables:

 Non-facility nodes are partitioned in two categories: leaves and intermediate nodes. Leaves are modeled using two-index binary variables h which also indicate the root node (facility) to which the leaf is allocated. For each i,j A V with j ai, we set ( 1 if node i is a leaf assigned to facility j, hij ¼ 0 otherwise: Intermediate nodes are modeled with binary variables q. For these nodes we do not explicitly indicate the facility they are

I. Contreras et al. / Omega 40 (2012) 847–860

assigned to. For each iA V we define ( 1 if node i is an intermediate node of some path; qi ¼ 0 otherwise:

 Finally, we define continuous variables g for the travel times. For every j A V, gj gives the maximal accumulated travel time of any path of R from a leaf to node j. In any rooted forest, for each j A V, at least one path starting at a leaf and containing j exists. If a facility is located at j, then j is a root of the forest, and the path terminates at j. Otherwise, if j is not a facility, the path continues towards a root node. In any case, gj represents the maximal accumulated travel time of any path that starts at a leaf and enters node j. As travel times are assumed to be non-negative, in any feasible solution, the maximum of all the gj values corresponds to the total travel time from some leaf node to the facility it is allocated to, and indicates the objective function value of the solution. CFLNDB can now be written as the (mixed) integer program min max g j , jAV

ðg i þ t ij Þyij r g j X

8ði,jÞ A A,

hij þqi þzi ¼ 1

8i A V,

ð13Þ ð14Þ

j A V\fig

X

X

hij þqi ¼

j A V\fig

yij

8i A V,

ð15Þ

ði,jÞ A A

hij þ qj r

X

ykj

8i,j A V,

ð16Þ

ðk,jÞ A A

yij r zj þ qj hij r zj

8ði,jÞ A A,

ð17Þ

8i,j A V, j a i,

ð18Þ

X X f i zi þ ckm ykm rB, iAV

g i Z0

ð19Þ

ðk,mÞ A A

8iA V,

hij ,zi ,qi A f0; 1g ykm A f0; 1g

ð20Þ 8i,j A V,

ð21Þ

8ðk,mÞ A A:

ð22Þ

Constraints (13) keep track of the longest accumulated travel time among all paths that terminate at node j A V. The set of constraints (14) establishes a partition of the set of nodes of the network: every node i A V is either the beginning of some path, an intermediate node, or it is used to locate a facility. The paths from nodes to facilities are defined by constraints (15)–(17). When iA V is a leaf or an intermediate node, constraints (15) ensure that exactly one arc leaving node i exists. However, when i A V is a facility, constraints (15) ensure, together with constraints (14), that there is no arc leaving node i. Inequalities (16) ensure that, when j A V is an intermediate node or a facility, then some arc entering node j exists. The set of complementary constraints (17) guarantees that the end node of an arc is either an intermediate node or a facility. Alternatively, taking constraints (14) into P account, inequalities (17) can be expressed as ðj,kÞ A A hjk þ yij r 1, indicating that when j A V is a leaf, no arc may enter node j. Constraints (18) ensure that no leaf is assigned to a node which is not a facility. Together with constraints (15), they also guarantee that no path exists which ends at a node which is not a facility. Constraint (19) is the budget constraint that guarantees that the total set-up costs for the location of facilities and the

851

construction of links does not exceed a pre-specified value. Finally, constraints (20)–(22) are the non-negativity and integrality constraints, respectively. Formulation (14)–(22) is a min–max non-linear mixed integer program that cannot be used directly by any commercial solver to obtain optimal solutions to CFLNDB. Thus, we linearize constraints (13) with the classical big-M technique, where M is an upper bound on the value of any feasible solution (e.g. the sum of all the travel costs). Then, we obtain the following mixed integer linear programming formulation: ðCP0 Þ

min T,

T Z gi

8i A V,

g j Z g i þt ij yij Mð1yij Þ

ð23Þ 8ði,jÞ A A,

ð14Þ2ð22Þ:

ð24Þ

Note that, by constraints (23), variable T is the maximum accumulated travel time among all nodes j A V and, thus, it represents the optimal value to the problem. Let CP denote the formulation that results from CP0 when the integrality conditions hij A f0; 1g are substituted by hij Z 0, i,j A V. An optimal solution to CP0 can be obtained by solving CP. If in CP the variables z, q and y remain binary, then constraints (14) and P (15) imply that in any feasible solution to CP, j A V\fig hij takes the value 0 or 1, for all iA V. Moreover, in such a solution the values of the variables h do not determine the objective function value. Thus, if in a feasible solution hij1 ¼ D1 40 and hij2 ¼ D2 40 for some i A V, then the solution with hij1 ¼ D1 þ D2 4 0 and hij2 ¼ 0 and all other variables as before, is also feasible and has the same objective function value. Therefore, in the following we will focus on the solution to CP. 3.3. Comparison of formulations We next perform a comparison, between the multicommodity flow-based formulation F and the alternative formulation CP, to observe their effectiveness and limitations for solving CFLNDB instances. Table 1 summarizes the numerical results obtained with both formulations on two sets of randomly generated instances (see Section 5 for a detailed description of the instances generation process), each of them with 16 instances of 10 and 20 nodes, respectively. All instances were solved with the MIP Xpress optimizer, under a Windows environment on a VAIO VGNSR19XN, with a processor Intel(R) Core(TM)2 P8400 at 2.26 GHz and 3 GB of RAM. We enabled the generation of cuts in order to observe the potential improvement of the LP relaxation bound by using them. Default values were used for all other parameters. A CPU time limit of 1 h was set. Table 1 has two blocks of columns, one for each dimension n A f10; 20g. In turn, each block has two columns, one for each of the two formulations F and CP. The first three rows give the average gaps in percent between the optimal value (or bestknown value, when the optimal value was not known at termination) and the value of (1) the LP bound of the formulation, (2) the LP bound of the formulation enhanced with the cuts generated at the root node of the search tree, and (3) the best lower bound at termination. The row labeled ‘‘max % g’’ gives the maximum gaps at termination among all the instances of the corresponding group. Row labeled ‘‘Avg. time’’ depicts the average computing times (in s), while rows labeled ‘‘mint’’ and ‘‘maxt’’ give the minimum and maximum computing times among all the instances of each subset, respectively. Entries in row labeled ‘‘Avg. nodes’’ are the average number of nodes explored in the search tree. Finally, the row labeled ‘‘Opt. solutions’’ gives the number of

852

I. Contreras et al. / Omega 40 (2012) 847–860

The travel times of the arcs used in feasible solutions yield trivial lower bounds on the accumulated travel times to reach the end-nodes of such arcs.

Table 1 Comparison of formulations CP and F. n ¼10

n¼ 20

CP

F

CP

F

100.00 41.45 0.00 0.00

91.96 45.72 21.51 85.19

100.00 43.51 4.65 37.35

94.48 84.79 78.43 100.00

6.82 1.06 52.63

1783.33 46.31 3600.00

1065.93 70.48 3600.00

3530.17 2834.66 3600.00

Avg. nodes

3274.23

3379.42

157,118.87

67.02

Opt. solutions

16/16

10/16

13/16

2/16

Avg. Avg. Avg. max

% % % %

gap LP gap LP þcuts gap 1 h g

Avg. time mint maxt

Proposition 1. For every ðj,kÞ A A, a valid inequality for Q is g k Zt jk yjk :

ð25Þ

These inequalities can be strengthened to take the travel time associated with chains defined by two consecutive arcs into account. In particular, in any feasible solution to CP, for given (i,j) and (j,k), the value t ij ðyij þ yjk 1Þ is positive only when the two consecutive arcs ði,jÞ and ðj,kÞ are used. Clearly, in this case, the accumulated travel time to reach node k will be at least t ij þ t jk . Thus, we have Proposition 2. For every (i,j) and ðj,kÞ A A, a valid inequality for Q is

instances of each group that could be solved to optimality within the time limit. The results of Table 1 provide clear evidence of the inherent difficulty of CFLNDP, even when considering small size instances. First, the LP bounds obtained with both formulations are extremely weak. The value of the initial LP bound of formulation CP is 0 in all of the 32 considered instances. The initial gaps of formulation F are over 90% in all of the 32 instances. Second, the addition of cuts generated by Xpress reduces the LP gap up to 40%, in both formulations. Third, none of the formulations allows us to solve optimally all the considered instances within 1 h of computing time. The results of Table 1 also provide clear evidence of the superiority of formulation CP over formulation F, in terms of their effectiveness for solving CFLNDP instances to optimality. With CP we obtain the optimal solution of all the 10 nodes instances, whereas with F we only obtain the optimal solution of 10 out of the 16 instances. The average gap at termination of the 6 unsolved instances is close to 22%, with a maximum gap at termination slightly over 85%. Also, observe that the average time to solve the 10 nodes instances is less than 7 s for CP while this average is above 30 min for F. Furthermore, these CPU times range between 1 and 52 s for CP, and between 46 and 3600 s for F. With the larger instances, the superiority of CP is even more evident. CP was able to solve to optimality 13 out of the 16 instances of 20 nodes, whereas F was only able to solve to optimality 2 of them. Moreover, at termination the average gap is less than 5% for CP and a little more than 78% for F. In view of these results we concentrate on improving CP in order to solve larger CFLNDP instances. In the next sections, we present some complementary techniques to obtain a stronger CP formulation. First, we derive several types of valid inequalities. Second, we indicate how tight upper bounds can be used both to further tighten some of the valid inequalities, and to fix the optimal values of some variables.

g k Zt ij ðyij þyjk 1Þ þt jk yjk :

ð26Þ

It is clear that inequalities (26) extend inequalities (25) to paths with two arcs. These inequalities can be extended further to longer paths. However, the number of inequalities increases with the number of arcs of the path and their contribution to the improvement of the lower bound does not compensate for the additional computational burden due to the increase of the size of the formulation. Proposition 3. For every i A V, a valid inequality for Q is X T Zg i þ t ij yij :

ð27Þ

ði,jÞ A A

Proof. Consider i A V. By constraints (14) and (15), in any feasible solution to CP, the out-degree of any node i A V is at most 1, so P P that ði,jÞ A A yij r 1. If ði,jÞ A A yij ¼ 0, then inequality (27) reduces to constraint (23). Otherwise, if yij ¼ 1 for some ði,jÞ A A, then inequality (27) reduces to T Z g i þ t ij , which is implied by constraints (24). & Since constraints (23) are dominated by inequalities (27), they can be directly replaced by them. P In any feasible solution to CP we have j A V hij r 1 for any i A V. Therefore, the following family of inequalities is also valid. Proposition 4. Let aij denote the length of the shortest path, relative to travel times t, from i to j in D. For every i A V, a valid inequality for Q is X aij hij : ð28Þ TZ jAV

For j A V, let bj ¼ minft rj þt jk 9ðr,jÞ,ðj,kÞ A Ag. If j is an intermediate node in a feasible solution to CP, there is (at least) one arc entering node j, and (exactly) one arc leaving node j. Therefore, the maximum travel time T must be at least bj . Therefore, we have: Proposition 5. For j A V, let bj ¼ minft rj þ t jk 9ðr,jÞ,ðj,kÞ A Ag. Then, a valid inequality for Q is

4. Valid Inequalities for CP In this section we derive some valid inequalities for CP. The addition of these inequalities does not change the feasible domain of CP, since each of them must be satisfied by all feasible integer solutions. However, they yield a reinforced formulation of the LP relaxation of CP, as fractional solutions to the linear programming relaxation of CP might fail to satisfy them, thus, producing improved LP bounds. As we will see, the addition of these inequalities allows to solve the CFLNDP more efficiently. Throughout the following, Q denotes the convex hull of feasible solutions to CP associated with rooted forests.

T Z bj qj :

ð29Þ

Taking into account that no leaf can be at the same time an intermediate node, inequalities (28) and (29) can be combined into one single family. Corollary 1. Let bj ¼ minft rj þt jk 9ðr,jÞ,ðj,kÞ A Ag and let aij denote the length of the shortest path from i to j in D. For every i A V, a valid inequality for Q is X TZ aij hij þ bi qi : ð30Þ jAV

I. Contreras et al. / Omega 40 (2012) 847–860

Inequalities (29) can be somewhat strengthened to take those arcs into account which actually define a given rooted forest as stated in the following result. Let r1 and k1, respectively, denote the indices of the predecessor and the successor of a given node j, in the incoming and outgoing arcs which determine the value of bj (with ties broken arbitrarily). Proposition 6. For every j A V, a valid inequality for Q is X ðt ji t jk1 Þyji : T Z bj qj þ

ð31Þ

Proof. Consider again a rooted forest associated with a feasible solution to CP and a node j A V. If j is an intermediate node of the P P rooted forest, ði,jÞ A A yij Z 1 and ðj,iÞ A A yji ¼ 1. If i0 A V is any of the predecessors of j in the rooted forest, then the actual travel time from i0 to the (unique) successor of j can be expressed as X X t ji yji ¼ t i0 j þ t jk1 þ ðt ji t jk1 Þyji t i0 j þ ðj,iÞ A A

X

ðj,iÞ A A

ðt ji t jk1 Þyji ¼ bj þ

ðj,iÞ A A

X

ðt ji t jk1 Þyji :

ðj,iÞ A A

Therefore, when j A V is an intermediate node, inequality (31) holds since in this case qj ¼1. Note also that, when j A V is not an intermediate node, the inequality also holds. If j is a facility, there exists no arc leaving j in the rooted forest, so that the right hand side of the inequality takes the value zero. Otherwise, j is a leaf, and the right hand side of the inequality is a trivial lower bound on the value of the arc that leaves node j in the solution. & Another reinforcement of constraints (29) gives rise to the family of inequalities T Z bj qj þ ðt ij t r1 j Þyij ,

8i,j A V:

ð32Þ

The rationale for inequalities (32) is quite similar to that of inequalities (31). Now we focus on the arcs entering node j in the current rooted forest. Given that several such arcs may exist simultaneously, it is no longer correct to accumulate the travel time of all of them to the travel time from node j to its successor. Therefore, each potential predecessor of node j gives rise to one valid inequality. Proposition 7. Define gij ¼ t ij þ minft ri þ t jk 9ðr,iÞ,ðj,kÞ A Ag, for every ði,jÞ A A. Then, a valid inequality for Q is T Z gij ðqi þ qj þ yij 2Þ:

ð33Þ

The validity of inequalities (33) follows as the right hand side of the inequality is positive only if both i and j are intermediate nodes and, in addition, arc ði,jÞ is used. In this case, it is clear that the accumulated travel time must be at least gij . Proposition 8. For every (i,j), (j,k) and ðk,mÞ A A a valid inequality for Q is T Z t ij yij þ t jk yjk þt km ðykm þyjk 1Þ:

Proposition 9. Let 0 ¼ H1i o H2i o    o HW i denote the different possible values of the travel times from i to other nodes in the original graph D, ordered by increasing values and let Sri ¼ fj A V,j ai9aij o Hri g. For iA V and r ¼ 2, . . . ,W i a valid inequality for Q is 0 1 X X r@ ð35Þ hij  z j A: T Z g i þ Hi qi þ jAV

ðj,iÞ A A

Z t r1 j þ t jk1 þ

853

ð34Þ

Inequality (34) is valid since its right hand side accumulates the travel times of all possible subpaths of the path ði,j,k,mÞ when some of its arcs, ði,jÞ, ðj,kÞ, or (k,m) are used consecutively in a feasible solution to CP. Inequalities (33) and (34) could also be extended to consider longer paths. However, similarly to inequalities (25), the computational burden caused by the increase in the size of the formulation does not compensate for the improvement on the lower bound associated with the linear programming relaxation.

j A Sri

Proof. By constraints (14), for a given i A I we have X X X qi þ hij  zj ¼ 1zi  zj : jAV

j A Sri

j A Sri

As Sri is the set of nodes that, in the original digraph, can be i reached from i in a total travel time strictly smaller than HW i , the above expression takes the value 1 in a feasible solution to CP only if, in the associated rooted tree, the travel time from node i to its allocated root node is at least Hri. Therefore, inequality (35) is valid, because gi represents the maximal accumulated time to node i, from any leave that is connected with its root node with a path going through to node i. & 4.1. Effectiveness of the valid inequalities For evaluating the effect of the valid inequalities we have incorporated them into formulation CP and run new series of computational experiments with the same test instances as before, plus 16 additional instances with 30 nodes. Since the rationale of the cover type inequalities (35) is quite different from all the other ones, we have run the new experiments in two separated phases, in order to evaluate the effect of inequalities (35) independently. Some preliminary testing indicated that the best combination of inequalities for the first phase was (25), (30) and (31), since for all other combinations the added improvement did not compensate for the increase in the CPU time. Therefore, for the first phase of the new experiments, formulation CP has been reinforced with inequalities (25), (30) and (31), while the second phase uses inequalities (35) in addition. Table 2 gives the results of the first phase. The meaning of the columns is as follows. Column ‘‘Opt’’ gives the number of instances (out of the 16 instances in the corresponding group) that could be solved optimally within 1 h. Column ‘‘% gap’’ depicts the average gaps between the LP and optimal values of CP in percent, i.e., 100ðOptLPÞ=Opt. (The value of the best-known solution was taken instead of the optimal value, for instances that could not be solved to optimality within the maximum time limit.) Columns ‘‘nodes’’ and ‘‘maxnod ’’ indicate, respectively, the average and maximum number of nodes explored in the branch and bound tree, while columns ‘‘time’’ and ‘‘maxt ’’ give the average and maximum running time in seconds over the 16 instances, respectively. Finally, column ‘‘maxg ’’ depicts the maximum gap between the optimal and the LP value, among all the instances of the group. As can be seen, a considerable improvement can be obtained by just including these families of inequalities. Even if the gap between the LP bound and the optimal value is, on the average,

Table 2 Results of formulation CP reinforced with inequalities (25), (30) and (31). n

Opt

% gap

maxg

Nodes

maxnod

Time

maxt

10 20 30

16 15 12

45.11 45.41 45.60

66.14 72.48 76.93

490.06 35,215.38 86,045.88

2975 319,671 332,584

1.64 376.42 1418.95

7.45 3600.00 3600.00

854

I. Contreras et al. / Omega 40 (2012) 847–860

still very big (over 45%), the computational cost for exactly solving the instances has decreased considerably. In particular, for the instances with n ¼10, the average CPU time required to solve the instances has been reduced from 7 s to less than 2 s. The improvement can be further appreciated for n ¼20. Now all but one instance could be solved within the maximal time limit. This improvement is also reflected in the average number of explored nodes, which has been reduced by one order of magnitude, as well as in the average required CPU time, which has decreased from above 15 min to a little more than 6 min. Despite the remarkable effect of the inequalities added in the first phase of the new experiments, in general, the improvement that was obtained was not good enough to solve all instances with n¼ 30 and nearly no instance with n ¼40 in the maximum time limit of 1 h. On the contrary, as illustrated in Table 3, when the cover type inequalities (35) are also used, the additional improvement allows not only to solve the smaller instances much more efficiently, but also to solve instances on up to n¼50 (see again Section 5 for the description of the instances) within the time limit. It is interesting to observe that inequalities (35) have a very small influence on the quality of the LP gap which, on average, is reduced by less than 5%, but an important impact on the overall effort needed to exactly solve the instances. Now the average number of explored nodes has decreased to less than 10% of the number of nodes for the instances of Table 2, resulting in average reductions on the required CPU times over 85%. Observe also that the computational effort required to solve exactly the new sets of 40 and 50 nodes instances is quite small, since the maximum number of explored nodes is less than 16,000 resulting in a maximum CPU time slightly over 15 min. Unfortunately, the reinforcement of formulation CP obtained with the different types of inequalities was still not good enough to allow us to handle instances of larger sizes. In the next section we describe alternative tools that we have used to further reinforce the formulation.

5. Computational study For our experiments all programs have been coded in C and run under Windows on a VAIO VGN-SR19XN, with a processor Intel(R) Core(TM)2 P8400 at 2.26 GHz and 3 GB of RAM. The solver Xpress v2.2.0 has been used to solve formulation CP. Default parameters were used, except for the branching rules where we first branched on the z variables, then on the q variables

Table 3 Results of the additional reinforcement with inequalities (35). n

Opt

% gap

maxg

Nodes

maxnod

Time

maxt

10 20 30 40 50

16 16 16 16 16

43.33 41.32 40.99 42.52 42.65

65.01 65.74 53.62 50.40 53.63

3.19 148.88 210.00 502.50 2693.31

14 1409 2143 4947 13,541

1.35 13.55 74.02 264.08 1018.13

2.95 33.84 276.98 1230.34 2955.41

and, finally, on the y variables. The time limit was set to 10 h of CPU time (36,000 s). 5.1. Test instances Our generation procedure is largely based on that of Melkote and Daskin [23,24] which, in turn, is based on that of Balakrishnan et al. [2]. We have generated networks with a number of nodes n A f10; 20,30; 40,50; 60,70; 80,90; 100g and, for each fixed number of nodes, all possible combinations with a number of links m A f2n,3n,4n,5ng and a budget B A f300n,400ng. For fixed values of n and m, first we randomly generate the (x,y)-coordinates of the nodes from a uniform distribution in ð½0; 100,½0; 100Þ, and the facility set-up costs from an integer uniform distribution in ½1200; 1700. Then, the desired number of links is randomly selected, and for each such link we calculate the Euclidean distance between its end nodes and round it to the nearest integer dij. The construction costs of the links are symmetric and are taken to be cij ¼ cji ¼ kdij ð1 þrÞ, where the proportionality factor varies in f2; 10g, and r is randomly generated from a uniform distribution in ½0:2,0:2. This procedure generates 16 instances for each fixed number of nodes, giving a total of 160 data sets. From each data set we have generated two different instances, one with symmetric travel times and another one with asymmetric travel times. In the first case, the travel time of each arc is given by t ij ¼ t ji ¼ dij ð1 þ rÞ where again r is randomly generated from a uniform distribution in ½0:2,0:2. In the second case, the travel times in each of the two possible directions are t ij ¼ dij ð1 þ r 1 Þ and t ji ¼ dij ð1 þr 2 Þ where, as before, r 1 ,r 2 are randomly generated from a uniform distribution in ½0:2,0:2. In total, we have generated a set of 320 instances. Each instance has been labeled as ‘‘InMkB’’, where ‘‘I’’¼‘‘S’’ for symmetric instances and ‘‘I’’¼‘‘A’’ for asymmetric ones; n is the number of nodes, ‘‘M’’ takes the values 2, 3, 4 or 5 depending on the relation m=n between the number of arcs and the number of nodes; ‘‘k’’ takes the value of the parameter k; and ‘‘B’’¼‘‘T’’ when the budget is 300n, whereas ‘‘B’’¼‘‘L’’ when the budget is 400n. Thus, for example, the instance labeled as S050410L has been generated with symmetric travel times, n ¼50 nodes, m ¼ 4  50 arcs, a k¼10, and a budget B ¼ 400  50. A summary of the design for the symmetric case is shown in Table 4. 5.2. Preliminary tests: effect of the upper bound As it is usual in min–max formulations, the quality of a known upper bound U is crucial to improve formulation CP in two complementary aspects. On the one hand, such a bound allows to obtain a tighter expression of constraints (24), as we can substitute the big M value on the right hand side of constraints (24) with the value U. On the other hand, we can apply simple elimination tests, derived from the validity of inequalities (25), (28) and (29), to reduce the size of the instances. Let U denote the value of the solution obtained with a heuristic. (i) (ii) (iii) (iv)

If If If If

U ot ij , then yij ¼ 0. U o aij , then hij ¼ 0. U o bi , then qi ¼0. there is no j A V, j ai such that t ij rM, then zi ¼1.

Table 4 Instance generation for the symmetric case. k A f2; 10g; r  U½0:2,0:2. n

m

B

f

c

t

10r,r ¼ 1, . . . ,10

2n, 3n, 4n, 5n

300n, 400n

[1200, 1700]

kdij

t ij ¼ t ji ¼ dij þ r

I. Contreras et al. / Omega 40 (2012) 847–860

Table 5 summarizes the effect of using a good upper bound U, both to reinforce constraints (24) and in the elimination tests. Columns denoted by ‘‘fix z’’, ‘‘fix q’’, ‘‘fix h’’, and ‘‘fix y’’, indicate the average number of z, q, h and y variables that have been fixed, respectively. In our experiments with instances on up to 50 nodes, the upper bound U is obtained with the GRASP heuristic of Contreras et al. [8].

855

The obtained results indicate that taking the value of the upper bound into account considerably reduces the overall computational burden required to solve the instances to optimality. In general, the improvement in the LP bound is moderate and the LP gaps remain above 30%. In contrast, the required CPU times have decreased by at least one order of magnitude and the maximum CPU time for any of the considered instances has decreased from

Table 5 Effectiveness of the upper bound. n

Opt

% gap

maxg

Nodes

maxnod

Time

maxt

Fix z

Fix q

Fix h

Fix y

10 20 30 40 50

16 16 16 16 16

36.87 27.84 33.48 30.81 37.99

65.00 56.36 52.51 44.75 42.63

1.25 44.71 31.36 413.25 332.38

5 589 280 3237 2561

0.15 0.80 2.58 10.68 31.66

0.52 4.80 9.12 66.43 152.06

0.13 0.29 0.29 0.50 0.63

3.50 7.75 10.62 14.93 20.38

45.75 286.00 709.65 1360.50 2204.63

9.88 22.29 30.86 47.43 53.25

Table 6 Computational results for medium-size symmetric instances. Opt

U

LP

% dev heur

% gap LP

Time

S050202L S050210L S050302L S050310L S050402L S050410L S050502L S050510L S050202T S050210T S050302T S050310T S050402T S050410T S050502T S050510T

16 19 17 17 14 16 17 18 20 21 20 21 19 19 22 21

19 19 17 21 16 17 18 20 21 21 21 21 20 19 23 23

9.21 12.25 10.00 10.09 10.16 9.18 11.22 11.67 11.71 14.37 11.73 12.95 11.45 11.02 13.25 13.73

15.79 0.00 0.00 19.05 12.50 5.88 5.56 10.00 4.76 0.00 4.76 0.00 5.00 0.00 4.35 8.70

42.44 35.53 41.18 40.65 27.43 42.63 34.00 35.17 41.45 31.57 41.35 38.33 39.74 42.00 39.77 34.62

23.32 5.29 20.47 70.83 15.21 21.37 12.37 25.68 29.09 7.16 9.38 46.04 41.23 9.92 25.72 106.87

S060202L S060210L S060302L S060310L S060402L S060410L S060502L S060510L S060202T S060210T S060302T S060310T S060402T S060410T S060502T S060510T

18 17 18 19 17 18 22 23 20 20 24 25 17 23 28 29

19 17 19 19 26 23 25 26 22 20 27 26 26 23 29 30

10.38 11.11 10.74 10.46 11.75 10.39 11.59 12.92 11.97 12.07 12.32 12.67 11.75 12.57 14.04 15.76

5.26 0.00 5.26 0.00 34.62 21.74 12.00 11.54 9.09 0.00 11.11 3.85 34.62 0.00 3.45 3.33

42.33 34.65 40.33 44.95 30.88 42.28 47.32 43.83 40.15 39.65 48.67 49.32 30.88 45.35 49.86 45.66

S070202L S070210L S070302L S070310L S070402L S070410L S070502L S070510L S070202T S070210T S070302T S070310T S070402T S070410T S070502T S070510T

18 16 15 20 17 16 19 19 21 19 20 23 20 22 26 26

19 18 18 21 18 19 21 26 23 21 21 27 24 22 29 31

10.36 10.13 8.74 10.96 9.93 9.78 10.50 11.32 11.82 11.61 10.76 13.23 11.46 11.66 13.14 13.93

5.26 11.11 16.67 4.76 5.56 15.79 9.52 26.92 8.70 9.52 4.76 14.81 16.67 0.00 10.34 16.13

42.44 36.69 41.73 45.20 41.59 38.88 44.37 40.42 43.71 38.89 46.20 42.48 42.70 47.00 49.46 46.42

Nodes

Fix z

Fix q

Fix h

Fix y

235 1 56 4453 1 98 1 334 552 1 1 2023 641 1 138 1490

1 1 1 0 1 1 0 2 1 0 0 0 0 0 0 0

21 27 24 10 29 19 27 24 16 22 14 5 20 15 17 15

2218 2258 2256 2200 2260 2242 2256 2228 2176 2212 2200 2074 2186 2190 2162 2166

50 54 54 16 98 102 94 80 32 26 24 6 46 70 36 38

979.57 34.75 47.49 42.06 12.40 33.77 7074.11 72.56 71.90 32.11 949.28 2003.80 13.28 24.27 569.72 90.22

18,976 229 286 418 1 65 249,202 1403 1160 26 48,930 53,104 1 1 16,810 2372

0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0

21 31 16 18 9 16 5 11 13 17 5 5 9 16 2 7

3214 3264 3200 3160 3138 3206 3192 3156 3116 3180 2982 2946 3138 3206 3096 3068

68 94 72 80 6 12 0 2 38 52 14 28 6 12 0 0

85.08 29.36 157.40 2606.60 46.89 39.95 63.15 50.56 175.70 238.32 138.47 6532.77 54.18 212.00 133.57 227.42

540 1 4502 54,717 116 1 314 74 2387 3512 1414 1,152,500 134 3777 2279 4764

1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0

16 20 14 18 21 17 10 4 10 10 9 11 9 10 2 1

4370 4424 4392 4334 4496 4488 4432 4292 4222 4298 4302 4088 4358 4402 4218 4132

52 64 76 54 20 8 2 2 18 32 38 18 8 4 0 0

856

I. Contreras et al. / Omega 40 (2012) 847–860

nearly 50 min to 152 s. These results suggest that the value of the LP gap is not a meaningful indicator of the actual difficulty of the problem, since instances with large LP gaps are solved efficiently in small computing times. Additional computational tests indicate that, when an upper bound is known, tightening inequalities (24) is equally important as applying the elimination tests since the improvement is considerably smaller when only inequalities (24) are reinforced or when only the elimination tests are applied. Next we will analyze in detail the numerical results that we have obtained in the computational experiments with instances on up to 100 nodes. 5.3. Computational experiments The algorithm that we apply in our computational experiments is the following. First, we apply the heuristic of Contreras and Ferna´ndez [8] to obtain an upper bound U, and we apply the variable elimination tests for fixing as many variables as possible. Then, we consider the formulation CP associated with

the resulting instance, reinforced with inequalities (25), (30), (31) and (35), where the value M in right hand side of constraints (24) is substituted by the upper bound U. Finally the solver of Xpress v2.2.0 is called for solving to optimality the reinforced formulation. It is clear that a valid upper bound can also be obtained by Xpress. However, we have observed that usually, it is not until many nodes have been explored in the search tree that the solver of Xpress produces relatively good bounds. Since for the overall optimization process it is important to obtain as soon as possible the benefits of having a reasonable upper bound, we use the heuristic of Contreras and Ferna´ndez [8] and reinforce the formulation CP before resorting to Xpress. The results of our experiments with the symmetric instances are summarized in Tables 6 and 7, whereas the results of the asymmetric instances are summarized in Tables 8 and 9. The meaning of the columns is as follows: columns ‘‘Opt’’ show the values of the optimal/best-known solutions. The following two columns, ‘‘U’’ and ‘‘LP’’, give the values of the feasible solution

Table 7 Computational results of large symmetric instances. Opt

U

LP

% dev heur

% gap LP

Time

Nodes

Fix z

Fix q

Fix h

Fix y

S080202L S080210L S080302L S080310L S080402L S080410L S080502L S080510L S080202T S080210T S080302T S080310T S080402T S080410T S080502T S080510T

16 15 17 17 15 16 20 21 19 19 22 21 19 20 23 28

18 16 21 19 17 17 22 25 20 21 25 24 23 23 27 32

9.06 8.99 10.45 9.96 9.83 9.43 11.00 11.92 10.52 10.17 12.28 11.91 10.90 10.74 13.23 14.35

11.11 6.25 19.05 10.53 11.76 5.88 9.09 16.00 5.00 9.52 12.00 12.50 17.39 13.04 14.81 12.50

43.38 40.07 38.53 41.41 34.47 40.06 45.00 43.24 44.63 46.47 41.52 43.29 42.63 46.30 42.48 48.75

240.02 86.77 204.41 437.64 27.85 84.33 861.06 147.87 3574.00 381.62 568.89 7438.72 91.59 113.04 6777.90 36,000.00

2087 235 1338 4990 1 165 18,706 1596 18,328 5785 8446 149,912 368 446 152,500 719,447

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

23 30 21 12 25 18 11 10 15 13 8 4 10 9 8 4

5766 5872 5700 5652 5918 5882 5728 5724 5668 5738 5460 5502 5734 5730 5444 5496

80 72 44 30 12 8 8 4 46 24 18 12 4 2 2 0

S090202L S090210L S090302L S090310L S090402L S090410L S090502L S090510L S090202T S090210T S090302T S090310T S090402T S090410T S090502T S090510T

16 14 18 17 15 15 19 17 19 19 18 23 18 18 23 24

17 18 18 19 19 18 21 21 21 21 23 23 21 19 30 25

8.87 8.96 9.45 9.23 8.73 9.18 10.09 10.38 10.04 10.23 10.80 10.72 10.03 10.54 11.96 12.14

5.88 22.22 0.00 10.53 21.05 16.67 9.52 19.05 9.52 9.52 21.74 0.00 14.29 5.26 23.33 4.00

44.56 36.00 47.50 42.13 41.80 38.80 46.89 38.94 47.16 39.82 43.16 53.39 44.28 41.44 48.00 47.22

1428.05 5752.97 36,000.00 3756.47 155.16 58.23 416.24 219.02 3797.82 13,439.35 36,000.00 36,000.00 110.01 87.53 1658.25 18,581.75

11,212 72,881 789,410 44,109 747 1 4468 1655 54,227 124,458 689,189 269,909 262 185 16,911 314,288

0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0

21 21 14 13 20 11 12 8 14 5 6 3 14 2 5 22

7404 7370 7198 7420 7510 7442 7396 7328 7198 7224 7002 7076 7484 7090 7214 9188

62 50 10 8 6 4 2 4 10 28 14 14 4 0 2 56

S100202L S100210L S100302L S100310L S100402L S100410L S100502L S100510L S100202T S100210T S100302T S100310T S100402T S100410T S100502T S100510T

16 16 18 17 14 16 20 17 17 19 19 20 17 18 24 22

17 20 18 18 17 18 22 22 19 22 21 22 18 20 25 24

8.63 8.86 8.81 9.14 8.28 8.28 10.05 9.44 9.58 10.22 10.20 10.54 9.34 9.39 10.06 11.06

5.88 20.00 0.00 5.56 17.65 11.11 9.09 22.73 10.53 13.64 9.52 9.09 5.56 10.00 4.00 8.33

46.06 44.63 51.06 46.24 40.86 48.25 49.75 44.47 43.65 46.21 46.32 47.30 45.06 47.83 50.71 49.73

2374.91 36,000.00 36,000.00 36,000.00 208.25 898.67 36,000.00 8238.82 4520.20 13,944.22 36,000.00 36,000.00 958.89 36,000.00 29,324.69 535.74

18,424 124,563 110,178 325,003 1111 4968 418,168 100,946 15,343 121,565 384,265 630,934 9268 425,134 404,316 3264

2 0 0 1 0 1 0 0 0 1 0 0 1 0 0 0

15 19 15 13 16 10 6 6 11 11 12 8 10 11 2 3

8918 9018 9112 9264 9300 9212 9172 9138 9040 8774 8766 8834 9212 9220 9044 9028

44 38 26 8 8 8 2 2 22 34 14 8 8 2 2 0

I. Contreras et al. / Omega 40 (2012) 847–860

857

Table 8 Computational results of medium-size asymmetric instances. Opt

U

LP

% dev heur

% gap LP

A050202L A050210L A050302L A050310L A050402L A050410L A050502L A050510L A050202T A050210T A050302T A050310T A050402T A050410T A050502T A050510T

16 21 18 23 18 19 15 23 17 26 22 28 21 23 19 27

18 24 21 27 21 22 17 24 18 30 27 29 23 27 19 31

9.95 12.25 11.65 12.81 11.58 12.01 10.33 13.12 12.13 14.23 12.49 15.64 12.80 13.79 10.91 15.08

11.11 12.50 14.29 14.81 14.29 13.64 11.76 4.17 5.56 13.33 18.52 3.45 8.70 14.81 0.00 12.90

37.81 41.67 35.28 44.30 35.67 36.79 31.13 42.96 28.65 45.27 43.23 44.14 39.05 40.04 42.58 44.15

16.75 36.44 19.89 23.32 343.31 22.14 14.71 74.07 7.21 71.23 34.55 47.92 36.15 35.60 26.60 1077.96

A060202L A060210L A060302L A060310L A060402L A060410L A060502L A060510L A060202T A060210T A060302T A060310T A060402T A060410T A060502T A060510T

17 23 18 22 15 18 17 17 21 27 20 27 19 22 18 21

20 24 21 22 19 20 19 20 22 30 23 29 22 26 21 25

10.74 12.63 10.80 12.90 9.57 10.92 9.95 9.33 12.83 15.34 12.40 15.93 10.94 12.74 11.65 11.05

15.00 4.17 14.29 0.00 21.05 10.00 10.53 15.00 4.55 10.00 13.04 6.90 13.64 15.38 14.29 16.00

36.82 45.09 40.00 41.36 36.20 39.33 41.47 45.12 38.90 43.19 38.00 41.00 42.42 42.09 35.28 47.38

A070202L A070210L A070302L A070310L A070402L A070410L A070502L A070510L A070202T A070210T A070302T A070310T A070402T A070410T A070502T A070510T

15 19 17 17 15 18 15 15 17 24 20 24 17 22 18 21

19 23 19 22 17 19 17 18 21 30 22 26 22 25 21 23

9.26 10.89 9.93 9.97 8.84 9.79 10.11 8.63 11.05 13.55 11.43 12.59 10.35 11.82 11.30 10.57

21.05 17.39 10.53 22.73 11.76 5.26 11.76 16.67 19.05 20.00 9.09 7.69 22.73 12.00 14.29 8.70

38.27 42.68 41.59 41.35 41.07 45.61 32.60 42.47 35.00 43.54 42.85 47.54 39.12 46.27 37.22 47.15

obtained with the heuristic and the values of the LP relaxation, respectively. Columns ‘‘% dev heur’’ give the values of the percent deviation between the best-known solution and the optimal value, i.e., 100ðUOptÞ=U, whereas ‘‘% gap LP’’ provide the values of the percent deviation of the optimal/best-known solution and the value of the LP relaxation, i.e., 100ðOptLPÞ=Opt. The meaning of the remaining columns is the same as in the tables described in the previous sections. For the instances that did not reach the time limit of 36,000 s, the entries in columns ‘‘Opt’’, are their optimal values, whereas for the other instances, these entries are the values of the best solution found, although optimality of such solutions could not be proven. As can be seen, we solved all medium-size instances, both symmetric and asymmetric (except for instance A070510T), although as the sizes of the instances increased the time limit was not enough, especially for the larger 100 nodes instances. For all the instances that were not solved optimally after 1 h of CPU time, we recorded the value of the bestknown solution and the lower bound at that point. Thus, in the cases that the instances were still not solved to optimality at

Time

Nodes

Fix z

Fix q

Fix h

Fix y

147 520 68 51 7053 1 1 2430 1 2299 545 2003 310 141 174 36,957

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

17 9 19 8 18 12 24 11 17 4 5 6 14 5 17 5

2213 2160 2243 2124 2209 2193 2243 2116 2213 2061 2155 2092 2176 2099 2208 1928

22 4 18 0 27 37 106 67 22 1 4 0 12 13 79 25

44.16 125.36 35.40 29.63 288.95 68.87 33.14 40.37 26.33 256.07 16.72 83.41 37.14 310.64 38.03 300.67

409 4156 134 31 3186 894 90 205 1 12,189 1 2487 36 19,301 101 12,726

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

19 11 16 12 14 15 20 10 18 6 11 7 8 5 15 0

3247 3094 3244 3214 3190 3152 3184 3144 3212 2915 3198 3066 3099 2957 3131 2971

11 9 9 8 49 24 90 86 9 2 2 2 33 3 70 30

36.27 51.98 68.05 61.82 46.91 386.32 57.67 576.29 30.07 179.32 358.32 502.74 131.27 2660.00 117.70 36,000.00

1 152 414 150 106 9635 161 10,043 1 4856 7811 15,982 1318 81,932 925 1,690,759

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

13 10 21 8 16 11 28 12 11 0 12 6 9 3 13 4

4507 4377 4465 4375 4463 4393 4468 4382 4451 4184 4382 4271 4302 4189 4332 4195

6 1 18 3 70 35 125 86 0 0 8 2 29 10 60 34

termination after 10 h, we could compare the above bounds with the final ones. We observed that, in all cases, the improvement of the lower bound was nearly negligible (only in one case the increase was more than one unit), whereas in most of the cases there was some improvement in the value of the upper bound. Therefore, we assume that, probably, the best-known solutions are not optimal ones. However, in all cases, the gap between the best-known solution and the lower bound never exceeded two units, so best-known solutions are very close to the optimal ones. Next we discuss briefly how the characteristics of the instances may affect their difficulty for being solved. In our computational experiments we have observed that the increase in the difficulty for solving the instances as their dimensions grow depends strongly on the type of instance. In particular, some types of instances can still be solved, in general, within small CPU times even when their dimensions increase, whereas other types of instances become very time consuming already for not so large dimensions. Roughly speaking, it seems that asymmetric instances are somewhat easier to solve than symmetric ones.

858

I. Contreras et al. / Omega 40 (2012) 847–860

Table 9 Computational results of large asymmetric instances. Opt

U

LP

% dev heur

% gap LP

A080202L A080210L A080302L A080310L A080402L A080410L A080502L A080510L A080202T A080210T A080302T A080310T A080402T A080410T A080502T A080510T

15 17 15 16 15 17 14 20 18 22 19 21 18 21 15 20

19 19 19 22 20 22 17 23 21 27 22 25 21 25 18 23

8.32 10.24 9.22 9.52 9.14 9.73 8.77 11.44 9.75 12.24 10.51 11.36 10.64 11.91 9.97 11.44

21.05 10.53 21.05 27.27 25.00 22.73 17.65 13.04 14.29 18.52 13.64 16.00 14.29 16.00 16.67 13.04

44.53 39.76 38.53 40.50 39.07 42.76 37.36 42.80 45.83 44.36 44.68 45.90 40.89 43.29 33.53 42.80

A090202L A090210L A090302L A090310L A090402L A090410L A090502L A090510L A090202T A090210T A090302T A090310T A090402T A090410T A090502T A090510T

14 18 16 17 13 16 13 17 17 21 18 22 16 20 15 20

19 20 18 22 17 19 15 18 20 25 22 24 18 24 18 21

8.53 9.47 8.78 10.23 7.80 9.38 8.47 9.09 9.81 11.22 10.02 12.08 8.98 10.76 9.33 10.48

26.32 10.00 11.11 22.73 23.53 15.79 13.33 5.56 15.00 16.00 18.18 8.33 11.11 16.67 16.67 4.76

A100202L A100210L A100302L A100310L A100402L A100410L A100502L A100510L A100202T A100210T A100302T A100310T A100402T A100410T A100502T A100510T

15 15 13 16 13 17 12 15 17 19 15 21 15 18 14 19

19 18 17 20 17 19 15 19 21 25 20 26 19 22 17 20

8.18 8.66 8.23 8.96 7.79 8.61 7.75 8.34 9.29 9.96 9.27 10.55 8.73 9.77 8.60 9.64

21.05 16.67 23.53 20.00 23.53 10.53 20.00 21.05 19.05 24.00 25.00 19.23 21.05 18.18 17.65 5.00

Nodes

Fix z

Fix q

Fix h

Fix y

72.62 102.98 167.56 116.67 58.31 1443.24 147.18 1251.97 104.83 86.67 102.16 109.16 90.54 2313.24 894.19 1236.80

2242 6658 13,766 7597 11 237,274 10,248 175,137 659 50 439 584 274 64,834 4405 17,513

1 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0

7 14 14 7 7 7 23 6 6 3 8 3 4 4 18 6

5787 5833 5849 5741 5668 5511 5787 5537 5704 5570 5738 5620 5625 5371 5724 5537

8 9 2 3 17 12 102 25 5 0 1 1 14 7 78 25

39.07 44.29 45.13 39.82 40.00 41.38 34.85 46.53 42.29 46.57 44.33 45.09 43.88 46.20 37.80 47.60

103.93 8590.71 103.41 194.25 125.67 4856.88 183.21 36,000.00 125.21 529.04 2653.62 1014.78 266.60 36,000.00 747.63 36,000.00

151 139,324 273 1412 435 36,610 648 358,280 478 5353 24,709 12,390 2353 997,708 9558 394,389

0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0

10 6 11 4 12 11 25 13 9 0 6 2 10 4 14 8

7415 7383 7490 7329 7362 7388 7403 7295 7371 7134 7322 7235 7294 7102 7212 7081

7 6 12 4 30 30 122 72 7 3 3 3 22 3 75 30

45.47 42.27 36.69 44.00 40.08 49.35 35.42 44.40 45.35 47.58 38.20 49.76 41.80 45.72 38.57 49.26

241.32 128.31 973.27 1121.65 189.66 36,000.00 150.51 2612.67 299.64 28,087.39 217.79 527.26 985.44 33,978.30 412.39 36,000.00

958 130 5622 11,124 732 422,411 211 23,059 1407 329,204 583 5219 7466 308,327 1698 391,615

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

7 10 17 9 15 7 20 10 5 3 8 1 8 3 13 10

9147 9218 9218 9182 9159 8961 9260 8909 9037 8850 9056 8855 9019 8726 9140 8814

9 14 5 0 34 19 117 45 7 1 0 0 16 11 75 35

Symmetric Instances

3000.00

L T

2500.00

Time

L T

2500.00

2000.00

2000.00

1500.00

1500.00

1000.00

1000.00

500.00

500.00

0.00

Asymmetric instances

3000.00

0.00 m=2n, k=2

m=2n, k=10

m=3n, k=2

m=3n, k=10

m=4n, k=2

m=4n, k=10

m=5n, k=2

m=5n, k=10

m=2n, k=2

m=2n, k=10

m=3n, k=2

m=3n, k=10

m=4n, k=2

m=4n, k=10

m=5n, k=2

m=5n, k=10

Fig. 1. Comparison of average CPU times between instances ‘‘L’’ and ‘‘T’’.

Figs. 1–3 allow for a better comparison among the different classes of instances, symmetric (‘‘S’’) and asymmetric (‘‘A’’), and for a deeper analysis of the influence of the choice of parameters on the difficulty of the instances. In all cases, the left chart

corresponds to ‘‘S’’ instances and the right one to instances of class ‘‘A’’, so it is possible to observe whether or not the influence of the parameter’s choice is similar on both classes. The figures represent average CPU times of instances generated with fixed

I. Contreras et al. / Omega 40 (2012) 847–860

859

Symmetric Instances

Asymmetric instances

3000.00

3000.00

k=2 k=10

2500.00 2000.00

2000.00

1500.00

1500.00

1000.00

1000.00

500.00

500.00

0.00

k=2 k=10

2500.00

0.00

m=2n, L

m=2n, T

m=3n, L

m=3n, T

m=4n, L

m=4n, T

m=5n, L

m=2n, L

m=5n, T

m=2n, T

m=3n, L

m=3n, T

m=4n, L

m=4n, T

m=5n, L

m=5n, T

Fig. 2. Comparison of average CPU times between instances with k¼ 2 and k¼10.

Symmetric Instances

3000.00 2500.00 2000.00

Asymmetric instances

3000.00

m=2n

2500.00

m=3n m=4n

2000.00

m=5n

1500.00

1500.00

1000.00

1000.00

500.00

500.00

m=2n m=3n m=4n m=5n

0.00

0.00 k=2, L

k=10, L

k=2, T

k=2, L

k=10, T

k=10, L

k=2, T

k=10, T

Fig. 3. Comparison of average CPU times between instances with k¼ 2 and k¼10.

sets of parameters over all dimensions from 50 to 100. Furthermore, all instances requiring more than 3600 s of CPU have been considered as ‘‘difficult’’, so that for calculating these averages, CPU times over 1 h have been ‘‘truncated’’ to 3600 s. In this way, all values lie in a rather small range and relatively small differences can be appreciated at the same scale. Fig. 1 illustrates the influence of the value of the budget B. For ‘‘S’’ instances the average CPU times of the ‘‘T’’ instances are always bigger than those of the ‘‘L’’ instances, for each considered set of parameters, whereas for the ‘‘A’’ instances the higher difficulty of the ‘‘T’’ instances is not so evident, although, in general also ‘‘T’’ instances appear to require more computational time. Fig. 2 illustrates the influence of the parameter k, that relates the length of an arc with its travel time. In the class ‘‘A’’, instances with k¼10 are, in general, more difficult, although this is not so evident in the class ‘‘S’’. Finally, Fig. 3 illustrates the influence of the number m of arcs density on the difficulty of an instance. Now, instances with m ¼ 4n seem to be easier, and those with m ¼ 3n seem to be harder for the ‘‘S’’ class, whereas instances with m ¼ 3n seem to be easier and instances with m ¼ 4n or m ¼ 5n seem to be harder for the ‘‘A’’ class.

6. Conclusions and remarks In this paper we have studied the center Facility Location/ Network Design Problem. We have presented two integer programming formulations for the problem, which have been compared empirically. The multicommodity-based formulation was clearly outperformed by the formulation that exploits the structure of the problem, which has been successively reinforced with different types of valid inequalities. We have generated a set

of benchmark instances with both symmetric and asymmetric travel costs with different characteristics and varying number of nodes up to 100. The numerical results of the computational experiments indicate that, in general, symmetric instances are somewhat more difficult to solve than the asymmetric ones. However, in both cases, the vast majority of instances on up to 80 nodes were solved optimally within 1 h of CPU time. For instances with 90 and 100 nodes in many cases more CPU time was required, although most of them could be solved to optimality within 10 h of CPU time. Future research will focus on different formulations for the problem and especially those suitable for being addressed with column generation methodology. Another interesting topic for further research is to consider capacities both in the arcs or the facilities.

Acknowledgments This research has been partially supported by grant MTM200914039-C06-05 of the Spanish Ministry of Education and Science and by ERDF funds. The research of the first author has been partially supported through grant 197243/217997 from the National Council of Science and Technology, Me´xico and through contract DFG Re776/10-1 of the University of Heidelberg. The research of the second author has been partially supported through grant PR2008-0116 of the Spanish Ministry of Science and Education. This support is gratefully acknowledged. Thanks are due to two referees for their valuable comments. References [1] Alumur S, Kara BY. Network hub location problems: the state of the art. European Journal of Operational Research 2008;190:1–21.

860

I. Contreras et al. / Omega 40 (2012) 847–860

[2] Balakrishnan A, Magnanti TL, Wong RT. A dual ascent procedure for largescale uncapacitated network design. Operations Research 1989;37:716–40. [3] Campbell AM, Lowe T, Zhang L. Upgrading arcs to minimize the maximum travel time in a network. Networks 2006;47:72–80. [4] Church R, ReVelle C. The maximal covering location problem. Regional Science Association Papers 1974;32:101–18. [5] Cocking C, Fleßa S, Reinelt G. Locating health facilities in Nouna District, Burkina Faso. In: Haasis H-D, Kopfer H, Schnberger J, editors. Operations research proceedings; 2005. [6] Cocking C. Solutions to facility location–network design problems. PhD thesis. Germany: University of Heidelberg; 2008. [7] Contreras I, Ferna´ndez E. General network design: a unified view of combined location and network design problems. European Journal of Operational Research, doi:10.1016/j.ejor.2011.11.009, in press. [8] Contreras I, Ferna´ndez E, Reinelt G. Heuristic solutions for the center facility location network design problem with budget constraint, in preparation. [9] Cordeau J-F, Pasin F, Solomon MM. An integrated model for logistics network design. Annals of Operations Research 2006;144:59–82. [10] Costa A, Franc-a PM, Lyra Filho C. Two-level network design with intermediate facilities: an application to electrical distribution systems. Omega 2011;39:3–13. [11] Daskin MS, Hurter AP, Van Buer MG. Toward an integrated model of facility location and transportation network design. Working Paper. Evanston, IL, USA: The Transportation Center, Northwestern University; 1993. [12] Daskin MS. Network and discrete location: models, algorithms and applications. New York: John Wiley and Sons, Inc.; 1995. [13] Drezner Z, Wesolowsky G. Network design: selection and design of links and facility location. Transportation Research Part A 2003;37:241–56. [14] Elloumi S, Labbe´ M, Pochet Y. New formulation and resolution method for the p-center problem. INFORMS Journal on Computing 2004;16:84–94. [15] Ghiani G, Laporte G, Musmanno R. Fixed-charge, network design models. In: Introduction to logistics systems planning and control. Wiley; 2004. [16] Gollowitzer S, Ljubic´ I. MIP models for connected facility location: a theoretical and computational study. Computers & Operations Research 2011;38:435–49. [17] Hakimi SL. Optimum locations of switching centers and the absolute centers and medians of a graph. Operations Research 1964;12:450–9.

[18] Kariv O, Hakimi SL. An algorithmic approach to network location problems, I: the p-centers. SIAM Journal of Applied Mathematics 1979;37:513–37. [19] Kuehn AA, Hamburger MJ. A heuristic program for locating warehouses. Management Science 1963;9:643–66. [20] Klose A, Drexl A. Facility location models for distribution system design. European Journal of Operational Research 2005;162:4–29. [21] Li X, Aneja YP, Huo J. Using branch-and-price approach to solve the directed network design problem with relays. Omega 2012;40:672–9. [22] Megiddo N, Supowit KJ. On the complexity of some common geometric location problems. SIAM Journal on Computing 1984;13:182–96. [23] Melkote S, Daskin MS. An integrated model of facility location and transportation network design. Transportation Research Part A 2001;35: 515–38. [24] Melkote S, Daskin MS. Capacitated facility location/network design problems. European Journal of Operational Research 2001;129:481–95. [25] Melo MT, Nickel S, Saldanha-da-Gama F. Facility location and supply chain management—a review. European Journal of Operational Research 2009;196: 401–12. [26] Murawski L, Church RL. Improving accessibility to rural health services: the maximal covering network improvement problem. Socio-Economic Planning Sciences 2009;43:102–10. [27] ReVelle CS, Eiselt HA. Location analysis: a synthesis and survey. European Journal of Operational Research 2005;165:1–19. [28] ReVelle CS, Eiselt HA, Daskin MS. A bibliography for some fundamental problem categories in discrete location science. European Journal of Operational Research 2008;184:817–48. [29] Smith HK, Laporte G, Harper PR. Location analysis: highlights of growth to maturity. Journal of the Operational Research Society 2009;60:140–8. [30] Srivastava SK. Network design for reverse logistics. Omega 2008;36:535–48. [31] Syam S, Cˆote´ MJ. A location-allocation model for service providers with application to not-for-profit health care organizations. Omega 2010;38: 157–66. [32] Toregas C, Swain R, ReVelle C, Bergmann L. The location of emergency service facilities. Operations Research 1971;19:1363–73. [33] Wong RT. Probabilistic analysis of a network design problem heuristic. Networks 1985;15:347–63.