The capacitated multiple allocation hub location problem: Formulations and algorithms

European Journal of Operational Research 120 (2000) 614±631

www.elsevier.com/locate/orms

Theory and Methodology

The capacitated multiple allocation hub location problem: Formulations and algorithms Jamie Ebery a

a,b,1

, Mohan Krishnamoorthy

a,*

, Andreas Ernst

a,2

, Natashia Boland

b,3

CSIRO Mathematical and Information Sciences, Private Bag 10, Clayton South MDC, Clayton VIC 3169, Australia b Department of Mathematics and Statistics, University of Melbourne, Parkville VIC 3052, Australia Received 1 December 1997; accepted 1 June 1998

Abstract In this paper we consider and present formulations and solution approaches for the capacitated multiple allocation hub location problem. We present a new mixed integer linear programming formulation for the problem. We also construct an ecient heuristic algorithm, using shortest paths. We incorporate the upper bound obtained from this heuristic in a linear-programming-based branch-and-bound solution procedure. We present the results of extensive computational experience with both the heuristic and the exact methods. Ó 2000 Elsevier Science B.V. All rights reserved. Keywords: Location; Linear programming; Heuristics; Branch and bound; Hubs

1. Introduction Hub location problems generally arise in situations where trac (such as mail, telecommunication packets or airline passengers) must be transported from its point of origin to its destination, but where it is expensive or impractical to dedicate exclusive transport links to each origin±

* Corresponding author. Tel.: +61 3 95458042; fax: +61 3 95458080; e-mail: [email protected] 1 E-mail: [email protected] 2 E-mail: [email protected] 3 E-mail: [email protected]

destination (O±D) pair. The problem of locating such interacting hub facilities arises in many applications, some of which include the Civil Aeronautics Board (see [12±14,20,24±29]), airlines (see [2±4,9,23]), postal delivery services (see [12± 15,22]), emergency services (see [10,17]), rejuvenating processes [10], problems on a sphere [5], telecommunications (see [8,18]), designing the architecture of a distributed computing system [16], and applications in engineering [6]. To illustrate the situation, consider a postal service in which we have a set of post code districts that exchange mail (i.e., trac). Each postcode district can be represented by a node. Two postcode districts may be connected by a transportation link, represented as a directed arc. If there is

0377-2217/00/$ - see front matter Ó 2000 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 2 2 1 7 ( 9 8 ) 0 0 3 9 5 - 6

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

an arc joining post code districts i and j, then this could carry mail directly from i to j. The entire set of nodes and arcs is known as a network. Usually in hub location problems, there is a trac stream between each ordered pair of nodes. In other words there exists a stream of trac which must ¯ow from each node, to every other node. Typically it is very expensive to dedicate a transport medium to every trac stream. Hence, hubbing systems are used. Hubs, which are chosen from among the nodes (postcode districts) of the network, act as consolidation centres or sorting centres for trac (mail) and are commonly assumed to be fully interconnected. They facilitate ¯ow of trac and incorporate a volume discount on inter-hub ¯ow. Direct ¯ow between non-hub nodes is not allowed. Instead, ¯ow originating at (or destined for) a non-hub node must be collected at (or distributed by) a hub node. We stipulate that each trac stream must be routed via one or (at most) two hubs. There is an associated cost for routing ¯ow between any two nodes that is based on factors such as time or distance. In some cases there is also a ®xed cost of establishing a hub. There are several variants of hub location problems. For example, we may specify that p of the nodes are to be hubs. These are referred to as p-hub location problems. Hub location problems may include single or multiple allocation of nonhub nodes to hubs. If a node is singly allocated then all trac from that node travels directly to and from only one hub; if it is multiply allocated then its trac may be divided and may travel directly to or from more than one hub. Further, hub location problems may be capacitated or free of any capacities. This means that the quantity of ¯ow that can be redirected or sorted by a hub will be restricted to some maximum volume, known as its capacity and is denoted by Ck for hub k. For an extensive review and classi®cation of hub location problems see Campbell [10]. In this paper we study the capacitated multiple allocation hub location problem (referred to in this paper as the CMAHLP). The choice of problem is motivated by a postal delivery application. In postal delivery, the volume of mail that hubs (sorting centres) can sort is limited by time con-

615

straints. This, in turn, limits the volume of mail that can be sorted at mail sorting centres. Hence we introduce capacity restrictions on the hub nodes. These restrictions are only applied to the trac arriving at the hub directly from non-hub nodes as this is where the sorting takes place. A variant of this problem arises when capacities exist on both, collection as well as distribution (in other words, out-bound ¯ow as well as inbound ¯ow). In the CMAHLP there is also a cost associated with the establishment of a node as a hub. Consequently, the number of hubs for any given problem needs to be derived by the mathematical decision model. We de®ne the CMAHLP more precisely as follows: Given a complete graph G ˆ …N ; A†, where N ˆ f1; . . . ; ng is the node set and A ˆ N  N is the set of directed arcs. There is a trac stream requirement of Wij associated with each ordered pair of nodes …i; j† 2 A; Wij is the volume of the associated trac stream. We do not assume that Wij ˆ Wji or that Wii ˆ 0. Also we are given Cij the cost of sending one unit of trac directly from i to j. Cij is usually proportional to the distance between i and j and is therefore generally Euclidean. Each ¯ow has three separate components: collection (origin node-to-hub), transfer (inter-hub), and distribution (hub-to-destination node). Each component has the respective per unit cost coecient of v, a, d. Thus the cost per unit ¯ow for collection is vCik ; where i is a non-hub node and k is a hub. Similarly inter-hub transfer has cost aCkl where k, l are both hubs and distribution has cost of dClj where l is a hub and j is a non-hub node. Usually, due to the inter-hub transfer eciencies, a < 1 is a transfer discount factor, while v, d ˆ 1. In any case, we assume v, d > a. To the best of our knowledge, there is no solution algorithm or computational study of CMAHLP in the literature. The uncapacitated multiple allocation p-hub median problem (UMApHMP), which is perhaps the closest related problem, has been well researched. The UMApHMP has been studied by Aykin [2,4] who examines a model with capacities on arcs using exact and heuristic approaches. Campbell [10] develops formulations for several variants of UMApHMP. Campbell [11] de®nes two heuristics (ALLFLO and MAXFLO) for

616

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

solving this problem. Klincewicz [21] uses Lagrangean relaxation in conjunction with branch and bound techniques to solve UMApHMP. SkorinKapov et al. [28] give a MILP which tightens the formulation given by Campbell [10] and uses fewer constraints. Ernst and Krishnamoorthy [13] outline a heuristic using shortest paths and obtain exact solutions using two methods, namely explicit enumeration and branch and bound. Furthermore, Ernst and Krishnamoorthy [13] develop an mixed integer linear program (MILP) with O…n3 † variables. Ernst and Krishnamoorthy [14] develop clustering and shortest path schemes to obtain lower bounds which are used in a branch and bound algorithm. A MILP of CMAHLP was given by Campbell [10] which included constraints like Xijkm 6 Yk and Xijkm 6 Ym 8i; j; k; m. We can replace these constraints with their aggregations, as presented in the closely related UMApHMP formulations given by Skorin-Kapov et al. [28]. The result is a more effective formulation. The MILP for CMAHLP given by Campbell [10], restricts the sum of the trac entering any hub via collection and transfer to be less than the capacity of that hub. In our version of the model we only apply a capacity restriction on the volume of trac entering a hub via collection. The reasoning for this is that in postal delivery, once mail has been sorted after collection it does not need to be sorted again for distribution. For the model that we are considering we would simply replace " # X X Wij …Xijkm ‡ Xijmk † ÿ Xijkk 6 Ck Yk 8k 2 N i;j

m

…1† by the following constraints: XX X Wij Xijkm 6 Ck Yk 8k 2 N : i

j

…2†

m

We denote this formulation as CMAHLP-C and present it below. Formulation CMAHLP-C X XXXX Wij Cijkm Xijkm ‡ Fk Yk min i

j

k

m

k

…3†

subject to: XX Xijkm ˆ 1 8i; j 2 N ; k

…4†

m

X Xijkm 6 Yk

8i; j; k 2 N ;

…5†

8i; j; m 2 N ;

…6†

m

X Xijkm 6 Ym k

XX i

Wij

X Xijkm 6 Ck Yk

j

8k 2 N ;

…7†

m

Yk 2 f0; 1g 0 6 Xijkm 6 1

8k 2 N ;

…8†

8i; j; k; m 2 N ;

…9†

where the variable Xijkm denotes the fraction of ¯ow that travels from node i to node j via hubs located at k and m, the variable Yk ; 8k 2 N is de®ned by  1 if node k is a hub; Yk ˆ 0 otherwise: The per unit cost of Xijkm ¯ow is de®ned by Cijkm ˆ vCik ‡ aCkm ‡ dCmj

8i; j; k; m 2 N :

Fk is the ®xed cost associated with the establishment of a hub at node k and Ck is its capacity. Unless otherwise stated, all summation occurs over the set N . In the above formulation, Eq. (4) ensures that the total ¯ow …Wij † from i to j is transferred, while Eq. (5) and Eq. (6) guarantee that transfers only occur via hubs. Eq. (7) ensures that the capacities are being adhered to. With CMAHLP the hubs are capacitated and there is a ®xed cost associated with establishing a given node as a hub. Thus the model not only decides on the optimal location of the hubs and the allocation of non-hub nodes to hubs, but also the number of hubs that are required. Observe also that if the location of the hubs are known a priori, then CMAHLP reduces to a multicommodity ¯ow problem. In this paper we present a new MILP formulation for CMAHLP. This formulation uses fewer variables and constraints than that reported in Campbell [10] or the closely related UMApHMP given by Ernst and Krishnamoorthy [13]. We also

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

describe a new heuristic method for obtaining upper bounds. The optimal solutions to problems are obtained using CPLEX 4.0 [19] in conjunction with the upper bound provided by the heuristic. We test our algorithms on the Australia Post (AP) data set, which was ®rst presented by Ernst and Krishnamoorthy [12]. The ®xed costs and capacities for this data set were given by Ernst and Krishnamoorthy [15]. We also add ®xed costs and capacities to the Civil Aeronautics Board (CAB) data set and test our approaches on this extension to the CAB data set. 2. Mathematical formulations

X X Zik ˆ Wij X Xlji ˆ Wij

l

subject to: X Hk ˆ p; k

#

k

8i; j;

…13†

X X X Ykli ‡ Xkji ÿ Ylki ÿ Zik ˆ 0 j

l

Zik 6

8i; k;

…14†

l

X Wij Hk

8i; k;

…15†

j

Xlji 6 Wij Hl

8i; l; j;

…16† 8i; j; k; l:

…17†

Eq. (11) speci®es the number of hubs whilst Eqs. (12)±(14) represent the divergence equations for a network ¯ow problem. Eq. (15) and Eq. (16) ensures no ¯ow travels directly between non-hub nodes. Formulation UMApHMP-N has …2n3 ‡ n2 ‡ n† variables of which n are binary and …n3 ‡ 3n2 ‡ n ‡ 1† linear constraints. The formulation of CMAHLP that we now develop, implicitly determines the number of hubs that are required for the optimal solution. This is achieved by means of adding a ®xed cost for establishing a hub at node k. Thus the objective function becomes " XX X X v Cik Zik ‡ a Ckl Ykli min i

k

k

l

# XX X i ‡d Clj Xlj ‡ Fk Hk : l

j

…18†

k

Also, Eq. (11) is no longer necessary and can be removed. To take the capacities of the hubs into account we add the following set of constraints: X Zik 6 Ck Hk 8k: …19†

Formulation UMApHMP-N " XX X X v Cik Zik ‡ a Ckl Ykli min k

…12†

l

With minor notational changes and the above de®nitions, Ernst and Krishnamoorthy [13] stated UMApHMP as given below.

XX Clj Xlji ‡d

8i;

j

k

Xlji ; Ykli ; Zik P 0; Hk 2 f0; 1g

CMAHLP may be formulated as a MILP. As a starting point we refer to a formulation for UMApHMP published by Ernst and Krishnamoorthy [13]. We de®ne the parameters of the problem as follows: Cij the cost per unit ¯ow from i to j; Wij the ¯ow from i to j; Fk the cost of establishing a hub at node k; v; a; d are the collection, transfer and distribution coecients respectively; N the number of nodes; Ck the capacity of hub k; p the number of nodes to be selected as hubs. The variables in our formulations are: Zik the ¯ow from node i to hub k; Ykli the ¯ow from node i via hubs k and l; Xlji the ¯ow from node i to node j via hub l;  1 if node k is a hub; Hk ˆ 0 otherwise:

i

617

l

i

…10†

j

…11†

This is enough to specify CMAHLP completely. However, in practice it is found that if we aggregate Eq. (16) in the following manner: X X Xlji 6 Wij Hl 8l; j; …20† i

i

618

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

the optimal solution is found more eciently even though the LP bound in the new formulation is not as tight. Note that this new formulation has …4n2 ‡ 2n† constraints and …2n3 ‡ n2 ‡ n† variables. So a MILP formulation of CMAHLP is

k

‡d

k

XX j

l

subject to: X X Zik ˆ Wij

l

X

Clj Xlji ‡

l

X

i

Fk Hk

…21†

l

k

l

XXX i

j

l

dClj Xlji ‡

X

Fk Hk

…31†

k

subject to: 8i;

…22†

Xlji ˆ Wij

8i; j;

…23†

X Xlji ˆ Wij XX i

Zik 6 Ck Hk

8k;

…24†

Ykli ‡ X

X j

Xkji ÿ

X Ylki ÿ Zik ˆ 0

8i; k;

…25†

l

8i; j;

8i; k;

…26†

Xlji 6

X

Wij Hl

8l; j;

…27†

i

Xlji ; Ykli ; Zik P 0; Hk 2 f0; 1g

8i; j; k; l:

…28†

Note that this formulation has substantially fewer variables than a formulation of CMAHLP that is based on Campbell [10]. Campbell [10] shows usage of O…n4 † variables, whilst CMAHLP-N uses O…n3 † variables. Thus we hope to solve much larger problems with our formulation. We also explored a second formulation in which we reduce the number of variables by making a substitution of the form Zik ˆ

X l

Ykli

8i; k:

…29†

k

Ykli 6 Ck Hk Ykli 6

…32† 8k;

XX i

X X Ykli ˆ Xlji k

Wij Hk

l

XX i

j

i

j

l

Zik 6 X

‡

k

i

X

i

k

Formulation CMAHLP-F XXX …vCik ‡ aCkl †Ykli min

j

k

X

#

i

The new formulation is then given below.

Formulation CMAHLP-N " X X XX min v Cik Zik ‡ a Ckl Ykli i

We also reduce the number of constraints by removing Eqs. (22) and (26) and replacing Eq. (27) with XX XX Ykli 6 Wij Hl 8l: …30†

Wij Hl

…33† 8l;

…34†

j

8i; l;

…35†

j

Xlji ; Ykli P 0; Hk 2 f0; 1g

8i; j; k; l:

…36†

This formulation has …2n3 ‡ n† variables and 2n2 ‡ 2n constraints, which is approximately half the number of constraints of CMAHLP-N. See Section 4 for a comparison of the e€ectiveness of the two formulations.

3. A heuristic algorithm Here we develop a heuristic algorithm to solve CMAHLP. This is motivated by the knowledge that, even with the best formulations, only problems of a relatively small size can be solved optimally. Our heuristic is motivated by the knowledge that, if there are no capacities on the hubs, the all pairs shortest path (APSP) algorithm provides a

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

good basis for a heuristic approach for multiple allocation hub location problems (see [13,14]). Given a set of uncapacitated hubs, the APSP will optimally allocate nodes (multiply) to hubs. However, since we have hub capacities for the CMAHLP, we can develop a heuristic based on a method for hub selection, the APSP procedure and a rerouting procedure for ¯ow that violates capacity restrictions. We start o€ with an initial nomination of hubs. To this, we apply the APSP. If the capacities are satis®ed at any stage, we have a feasible solution for the CMAHLP and we note its value. If not, we reroute excess ¯ow using a heuristic procedure. Since the above approach is biased by our initial nomination of hubs, we attempt a greedy single-interchange procedure to investigate if other potential hubs (in the neighbourhood of the nominated hub set) yield better solutions. We repeat this process until some termination criteria are satis®ed. Furthermore, it is possible to use an exact multicommodity ¯ow approach to perform the allocation of non-hub nodes to the nominated hubs in an ecient manner. However, since this is a computationally expensive procedure, we only apply it to the ®ve best hub solutions found during the heuristic approach. We start by creating an ordered list L, which is a permutation of all the nodes. We refer to the nodes in L as prioritised hub candidates. Initially the candidates are ranked in descending order of Ci =Fi for each node i 2 N . We refer to the ith element of L as L…i†; for each i 2 N P . De®ne P W the total ¯ow of the system by: W ˆ i j Wij . De®ne H …L; h† ˆ fL…i†: i 6 hg where h 2 N . Initially the ®rst h P nodes are taken as trial hubs, where h ˆ minft: tiˆ1 CL…i† P Wg. Trial hubs are nodes that are chosen to be hubs in a given iteration of the heuristic. Suppose that at some iteration of the heuristic L ˆ …; ; ; . . . ; n† and h ˆ 3, then H …L; h† ˆ f1; 2; 3g, i.e. 1; 2 and 3 are the trial hubs. We also de®ne z…H † to be the objective value calculated by the heuristic (this might not necessarily be the true objective value) corresponding to H ˆ H …L; h†. The algorithm below gives an outline of the heuristic. We keep track of the ®ve best solutions found by the algorithm, denoted as H1 ; . . . ; H5

619

(Note that Hk ˆ …Lk ; hk †, where Lk is the list corresponding to the kth best solution and hk is the number of hubs in the set). The main for loop performs a descent in which hubs are swapped one at a time. In each iteration z is the best solution found so far. The ®nal two procedures of the heuristic perform a search in the vicinity of the three best solutions and evaluates the exact objective value for the ®ve best sets of hubs found. Main Heuristic begin initialise L to be nodes in descending P order of Ci =Fi t h Minft: iˆ1 CL…i† P Wg L Lo for (i ˆ 0 to 2) do 1 z if(i ˆ 2) P h Minft: tiˆ1 CL1 …i† P Wg ‡ 1 L L1 for (t ˆ 1 to h) do # swap every node with the tth hub Procedure Scenarios…t; h; H1 ; H2 ; H3 ;  H4 ; H5 ; L; z † repeat # add a hub until no improvements z z h h‡1 Procedure Scenarios…h; h; H1 ; H2 ; H3 ; H4 ; H5 ; L; z † while (z < z ) if(i ˆ 0) L Lo Procedure Cycle…0; L† Procedure The Boil…; H1 ; H2 ; H3 † Procedure Multicommodity…H1 ; H2 ; H3 ; H4 ; H5 † end. The following paragraphs describe the procedures that are referred to above. In any given procedure the input parameters (I) are given to the left of the semicolon, and the input/output parameters (I/O) are given to the right. For example procedure APSP(I,I;I/O,I/O).

620

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

3.1. Procedure Scenarios …t; h; H1 ; H2 ; H3 ; H4 ; H5 ; L; z † This procedure takes nodes in H …L; h†, uses them as trial hubs and evaluates the hub selection using a heuristic. Every other node is then swapped with the trial hub L…t† and the resulting selection of hubs is evaluated using the same heuristic. The hub receiving the ¯ow from i to j is stored in sij . The shortest path distance from i to j is given by dij . All ¯ow is initially routed along the shortest paths irrespective of hubs capacities. The ¯ows are then adjusted by re-routing some of the ¯ow away from hubs exceeding their capacity, to obtain a feasible solution. A feasible solution can always be found since only those scenarios with total hub capacity greater than total ¯ow are evaluated. It is important to note that we always ensure that there are no duplicate solutions in the set of ®ve best solutions. This is important for the e€ectiveness of ThefBoil and the Multicommodity procedure. 

Procedure Scenarios …t; h; H1 ; H2 ; H3 ; H4 ; H5 ; L; z † begin repeat n ‡ 1 ÿ h times Procedure P P APSP…h;L; d; s† z i j dij Wij if(a hub has exceeded it's capacity) Procedure ReRouting…h;L; s; d; z† z if(z > z) then z if(z < z…H5 †) then update H1 ; . . . ; H5 Procedure Swap…t; h ‡ 1;L† Procedure Cycle…h;L† return 3.2. Procedure APSP…h; L; d; s† Given a set of hubs, this procedure creates a network with three layers: All nodes in the bottom layer, hub nodes in the centre layer and all nodes in the top layer. Nodes in layers 1 and 3 are connected to those in layer 2 and nodes in layer 2 are fully inter-connected. The arc costs in this network are given by the collection (layer 1 to layer 2), transfer (layer 2 to layer 2) and distribution (layer 2 to layer

3) costs. It then uses the Floyd±Warshall ``All Pairs Shortest Path'' algorithm (see Ref. [1]) to assign nodes to hubs irrespective of capacities. For details on the use of APSP to hub location problems, refer to Ref. [13]. sij represents the successor of node i, that is the hub which sorts the trac from i to j. The shortest distances from i to j are returned in dij . Also note that this algorithm has complexity O…hn2 †, where h is the number of trial hubs. 3.3. Procedure ReRouting (h; L; s; d; z) Rerouting is performed by removing hubs from the network that have reached their capacity and rerouting the excess ¯ow via hubs that have not yet  to denote the set reached their capacity. We use H of hubs that have not yet reached their capacity. The excess ¯ow is rerouted via the cheapest alter (Note that it is possible for this alnative in H ternative path to include two new hubs that were not used before). However this procedure does not guarantee optimality, since the ¯ow to be rerouted is chosen arbitrarily. Optimal solutions could be obtained by solving a linear program which is too computationally expensive to perform for each iteration of the heuristic. Note however that we do solve the ¯ow problem exactly for the ®ve best sets of hubs. It should also be noted that Pn if a node's outgoing ¯ow exceeds its capacity … jˆ1 Wij > Ci † then it is never selected as a potential hub since such a node is unlikely to be a hub in the optimal solution. We use bk to denote the volume of ¯ow that is collected by hub k. Procedure ReRouting…h; L; s; d; z† begin P P bk i j:sij ˆk Wij  H fk: bk < Ck g while(9k s.t bk > Ck ) do for (i,j) 2 N  N if(bsij > Csij ) argmin k  cil ‡ dlj l2H / minfWij ; bsij ÿ Csij ; Ck ÿ bk g z z ‡ /…cik ‡ dkj ÿ dij † bsij ÿ / bsij bk ‡ / bk   ÿ fkg H if(bk P Ck ) then H return

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

3.4. Procedure Swap (i; j; L) This process simply swaps the node L…i† with node L…j†. So if L ˆ …1; 2; 3; 4; 5; 6; 7; 8† after one application of Swap…1; 5;L† we would now have L ˆ …5; 2; 3; 4; 1; 6; 7; 8†.

3.5. Procedure Cycle (h; L) A Cycle is the following sequence of Swaps: Procedure Cycle…h;L† begin for (k ˆ 1 to n ÿ h ÿ 1) do Procedure Swap…h ‡ k; h ‡ k ‡ 1;L† return; So if we have list L ˆ …1; 2; 3; 4; 5; 6† with h ˆ 3, after one application of Cycle, L ˆ …1; 2; 3; 5; 6; 4†.

Procedure

H2 ; H3 ; H4 ; H5 ;L; z †

621 

Scenarios…t; h; H1 ;

until no improvement in H1 ; . . . ; H3 return 3.7.Procedure Multicommodity …H1 ; H2 ; H3 ; H4 ; H5 † This procedure takes the ®ve best solutions found and solves the Linear Programming formulation given by CMAHLP-F with the Hk ®xed accordingly for each solution. That is we solve the LP for each of the ®ve solutions using the hubs that were found for each corresponding solution. We selected the CMAHLP-F in preference to the CMAHLP-N formulation since in practice CMAHLP-F is considerably more ecient than CMAHLP-N for the case where the hubs are known a priori. We solve the linear programs using CPLEX 4.0 [19]. 4. Computational results

3.6. Procedure The Boil …; H1 ; H2 ; H3 † The The Boil ®nds three local minima by using a descent procedure starting from the three best solutions. The ®nal solutions returned by the procedure are minimal with respect to the neighbourhood operation of swapping any hub node with any non-hub node. For solution H1 we accept any improvement, while for solutions H2 and H3 we use a steepest descent method. Note that the algorithm never searches the same neighbourhood twice as H1 ; . . . ; H3 are always di€erent and because our implementation includes a check to skip any of these if they were already searched in a previous iteration of the repeat loop. Procedure The Boil…; H1 ; H2 ; H3 † begin repeat for (j ˆ 1 to 3) do Lj Lj for (h ˆ hj ÿ 1 to hj ‡ 1) do if(j ˆ 1) then L Lj for (t ˆ h to 1) do if(j > 1) then L Lj

In this section we present results for the heuristic that is described in this paper. We also discuss the e€ectiveness of the three formulations of CMAHLP discussed in this paper: CMAHLP-C, CMAHLP-N and CMAHLP-F. All computational experiments were performed on a DEC 3000/700 (200MHz alpha chip). We solved the three formulations using CPLEX 4.0 [19] with the preprocessor and heuristic turned o€. This was because the pre-processor seldom aided the problem and the CPLEX 4.0 [19] heuristic was inferior to the one that we described in this paper. When solving any problem we used the solution obtained from our heuristic (which was coded in C) as an upper bound. We do not use the con®gurations obtained by the heuristic to aid any of the exact approaches. Our computational experiments were carried out using the CAB and the AP data sets. The CAB data set comes from The Civil Aeronautics Board, USA and has been the most popular data set used to benchmark algorithms for p-hub problems. It includes problems of size n ˆ 10; 15, 20 and 25 each with a 2 f0:2; 0:4; 0:6; 0:8; 1:0g and v; d ˆ 1. For a detailed description of the CAB data set see

622

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

[25]. This data set however, does not contain ®xed costs or capacities for nodes, so we generated them using a model that is described in Section 4.1. The AP data set is derived from a real-world hub location problem for Australia Post and is described in Ernst and Krishnamoorthy [12]. The AP data set, which includes ®xed costs and capacities for the nodes, can be obtained from OR-Library [7]. This data set consists of 200 nodes, which represent postcode districts. Problems of size n ˆ 10; 20; 25; 40; 50; 100; 200 can be obtained from the data set. Other signi®cant di€erences between the AP and CAB data set are that with AP in general we have the following: v 6ˆ d, Wij 6ˆ Wji , Wii 6ˆ 0. 4.1. Generation of ®xed cost and capacities Given an arbitrary data set which contains the …x; y† co-ordinates and the volumes of ¯ows Wij for nodes, we generate capacities for each node according to   n 3di Oi Oi ; ‡ …37† Ci ˆ p 5dm Om where di is the distance of node i from the centre of mass of the system, Oi is the total out¯ow from node i, dm and Om are the maximum values of di and Oi respectively, and p is an integer, which can be used to control the number of hubs in the optimal solution. When p is increased, we could expect more hubs to be present in the optimal solution. The ®xed costs are calculated via two methods. Firstly it is evaluated as a function of the distance from the centre of mass according to   3di : …38† Fid ˆ fo 1 ÿ 5dm Secondly we evaluate ®xed costs as a function of capacity according to     Ci ‡ Oi 1 ‡ : …39† Fic ˆ fo 5 Cm ‡ Om 2 Cm is the maximum value of Ci . Here, fo represents the scaled di€erence in objective value between a scenario in which a central node (the node closest

to the centre of mass) is a hub and a scenario in which all nodes are hubs. The function is scaled based on p. That is, XX

fo ˆ

i

ÿ

…vCih ‡ dChj †Wij

j

XX i

!, aCij Wij

p;

…40†

j

where h is the node that is closest to the centre of mass. In the ®rst part of the above function, we need only to collect and distribute ¯ows as there is only one hub. In the latter part, since we assume that all nodes are hubs, we only need to transfer ¯ows. We use these functions to generate ®xed costs and capacities for the CAB data set. We obtain problems of the form nFp. In this notation, n indicates the number of nodes in the problem and p is an integer. F is a label that indicates the type of ®xed cost model that is used. It takes the label d if the ®xed costs are derived based on distance from the centre of mass (as given by Eq. (38)). It takes the label c is the ®xed costs are dependent on capacities (as given by Eq. (39)). For each problem type, the ®xed costs and capacities have been derived in such a way that the number of hubs in the optimal solution will generally increase as p increases. For the AP data set, the ®xed costs and capacities have already been provided by Ernst and Krishnamoorthy [15]. Two types of ®xed costs, loose (L) and tight (T) are considered. In the data set with tight ®xed costs, nodes with larger ¯ows have higher ®xed costs. So in general it is harder to choose which nodes should be hubs in such cases. Similarly two types of capacities, loose (L) and tight (T) are considered. So for each problem with n nodes, there are four cases to be considered: LL, LT, TL and TT. A problem in the table is denoted by nFC, where n is the number of nodes. F and C are the types of ®xed costs and capacities respectively. In Tables 1 and 2 we give the problem number, problem name, the optimal solution value along with the corresponding set of optimal hubs for each problem that we have considered. In all, Tables 1 and 2 present the problem data for 100

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

623

Table 1 CMAHLP solutions for the AP data set Problem Opt Sol Opt Hubs

10LL 221032.70 147

10TL 257558.06 4 5 10

10LT 246495.03 1 4 5 10

10TT 257558.06 4 5 10

Problem Opt Sol Opt Hubs

20LL 230385.47 7 14

20TL 266877.50 7 19

20LT 246433.77 6 8 14

20TT 287574.78 1 10 19

Problem Opt Sol Opt Hubs

25LL 233636.69 8 18

25TL 305982.08 9 23

25LT 263760.56 9 12 19

25TT 333221.75 12 14 24

Problem Opt Sol Opt Hubs

40LL 238313.45 14 29

40TL 297204.10 14 19

40LT 261547.11 14 26 30

40TT 339965.84 14 25 38

Problem Opt Sol Opt Hubs

50LL 234399.93 15 35

50TL 314845.56 3 24

50LT 264999.73 6 26 32 48

50TT 391287.33 6 26 41 48

test problems. This includes 20 AP problems and 80 CAB problems. Note that for the CAB problems, there are 16 for each value of a ˆ 0:2; 0:4; 0:6; 0:8:1:0. The optimal solutions and the optimal locations of hubs that are presented in Table 1 and Table 2 are presented here for the sake of completeness. We obtain these using our MILP approach. The time required to obtain these solutions and other relevant information are presented in subsequent discussions. From Table 2 we observe that, the value of p that we chose in our ®xed cost models (38, 39) are generally equal to (and in some cases, less than) the number of hubs present in the optimal solution. One exception is Problem 20c2. We also generated problems with larger values of p in the ®xed cost models (for example 15d9, 25c14 and so on). In all of these cases, the actual number of hubs that were obtained was less that the value of p chosen in the corresponding ®xed cost model. 4.2. Performance of heuristic In this section, we test the e€ectiveness of our heuristic. We present the results of computational experiments on 20 problems from the AP data set and 16 problems from the CAB data set (with

a ˆ 0:2). While we use these 36 problems to analyse the e€ectiveness of our heuristic, we also present the performance of our heuristic over all of the other CAB problems in Tables 4, 6 and 7. In Table 3 `Prob' indicates the problem name according to nFC or nFp depending on whether the problem is from the AP or CAB data set respectively. In the application of our heuristic, it is possible (at any stage of the algorithm) for a set of allocations to violate the capacity constraints on the chosen hubs (although the total ¯ow across the network will not exceed the sum of all the hub capacities). In such cases, we apply the rerouting procedure. `Gap_r' is the upper gap obtained if only the rerouting procedure is applied to all trial hubs found by the heuristic. Additionally, it is possible for us to retain the ®ve best sets of hubs obtained during the algorithm. We can then apply an exact multicommodity ¯ow algorithm to these sets of hubs to resolve the allocations while satisfying all the capacity constraints. In general, we could apply the multicommodity ¯ow procedure to every set of hubs obtained by the heuristic. In practice this would be far too computationally expensive. Hence, we only apply it to the ®ve best sets of hubs obtained by the heuristic. `Gap_u' is the upper gap obtained if the multicommodity ¯ow procedure is applied at the end.

624

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

Table 2 CMAHLP solutions for the CAB data set (a ˆ 0:2; 0:4; 0:6; 0:8; 1:0) Problem (a ˆ 0:2) Opt Sol Opt Hubs Problem (a ˆ 0:2† Opt Sol Opt Hubs

10c2 2843.50 9 10 20c2 3149.80 6 18 19

10d2 1224.66 79 20d2 1770.20 12 17

10c4 2500.00 2 4 10 20c4 2759.83 4 13 18 19

10d4 888.28 3478 20d4 1319.33 4 12 17

15c2 4073.77 9 15 25c2 3470.20 8 25

15d2 1775.98 4 12 25d2 1944.89 17 22

15c4 3601.53 4 8 13 25c4 3124.22 4 8 18

15d4 1401.80 4 12 14 25d4 1509.55 4 12 17

Problem (a ˆ 0:4) Opt Sol Opt Hubs Problem (a ˆ 0:4) Opt Sol Opt Hubs

10c2 2569.56 9 10 20c2 2880.29 6 18 19

10d2 1184.20 79 20d2 1735.38 12 17

10c4 2260.15 2 4 10 20c4 2534.39 4 10 18

10d4 899.44 347 20d4 1319.37 4 12 17

15c2 3673.83 5 15 25c2 3207.10 8 25

15d2 1720.46 4 12 25d2 1918.98 12 17

15c4 3267.62 4 8 13 25c4 2910.25 4 8 18

15d4 3263.10 4 12 14 25d4 1517.24 4 12 17

Problem (a ˆ 0:6) Opt Sol Opt Hubs Problem (a ˆ 0:6) Opt Sol Opt Hubs

10c2 2280.26 9 10 20c2 2592.73 6 18 19

10d2 1123.01 79 20d2 1675.89 12 17

10c4 2010.84 2 4 10 20c4 2308.97 4 10 18

10d4 885.75 347 20d4 1308.51 4 12 17

15c2 3270.23 9 15 25c2 2932.14 8 25

15d2 1659.78 4 12 25d2 1865.37 4 12 17

15c4 2918.41 4 8 13 25c4 2686.76 4 8 18

15d4 1393.28 4 12 14 25d4 1508.23 4 12 17

Problem (a ˆ 0:8) Opt Sol Opt Hubs Problem (a ˆ 0:8) Opt Sol Opt Hubs

10c2 1972.85 9 10 20c2 2291.47 6 18 19

10d2 1048.15 79 20d2 1595.66 7 17

10c4 1744.35 2 4 10 20c4 2063.40 4 10 18

10d4 867.11 347 20d4 1281.88 4 12 17

15c2 2855.16 79 25c2 2650.70 8 25

15d2 1587.89 4 12 25d2 1794.66 12 25

15c4 2568.23 4 8 13 25c4 2788.48 4 8 18

15d4 1373.55 4 12 14 25d4 1480.42 4 12 17

Problem (a ˆ 1:0) Opt Sol Opt Hubs Problem (a ˆ 1:0) Opt Sol Opt Hubs

10c2 1657.65 9 19 20c2 1961.60 5 7 18

10d2 969.31 79 20d2 1473.07 7 17

10c4 1461.73 2 4 10 20c4 1792.43 4 10 18

10d4 839.31 347 20d4 1242.88 4 12 17

15c2 2406.35 79 25c2 2361.57 8 25

15d2 1509.07 4 12 25d2 1669.93 12 25

15c4 2200.27 4 8 13 25c4 2188.46 4 8 18

15d4 1340.68 4 7 12 25d4 1427.53 4 12 17

Here `Gap' is de®ned as the relative deviation of the heuristic solution from the optimal solution, expressed as a percentage. We also present `CPU_r' and `CPU_u' which re¯ect the computational e€ort required to obtain `Gap_r' and `Gap_u' respectively. From Table 3 we observe that the rerouting procedure on its own is not sucient to yield good upper bounds. The heuristic with only the rerouting procedure is able to produce an average upper gap of 0.78% over all the problems in the table. In contrast, the heuristic with the multicommodity ¯ow algorithm ®nds the optimal solution in all but one of the test problems and has an average upper gap of 0.04%. Of course, the improvement in the upper gap (obtained by solving the multicommodity ¯ow)

requires more computational e€ort. In some problems CPU_u is as much as 20 times CPU_r. However, from our experiments, we found that even a 0.25% decrease in the upper gap often results in signi®cant savings in the overall time to solve the problem exactly (using any of the formulations that we described). Thus, the increased computational e€ort expended in getting a tighter upper bound at the root node pays o€ in reduced overall times for solving CMAHLP exactly. This is because the branch and bound procedure is able to prune substantial segments of the search tree. So, despite the computational overhead, we can conclude that it is useful to perform the multicommodity ¯ow at the end of the heuristic. In general, we found that problems in the AP data set with tight capacities and ®xed costs are

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631 Table 3 Heuristic results for AP and CAB …a ˆ 0:2† problems Prob

Gap_r (%)

CPU_r (s)

Gap_u (%)

CPU_u (s)

10LL 10TL 10LT 10TT 20LL 20TL 20LT 20TT 25LL 25TL 25LT 25TT 40LL 40TL 40LT 40TT 50LL 50TL 50LT 50TT 10c2 10d2 10c4 10d4 15c2 15d2 15c4 15d4 20c2 20d2 20c4 20d4 25c2 25d2 25c4 25d4

0.00 0.00 0.00 0.00 0.00 0.00 1.72 1.13 0.18 0.00 0.71 4.37 0.00 0.18 0.66 1.27 0.00 0.22 1.45 5.84 1.29 3.34 0.45 0.00 0.31 0.00 0.52 0.56 1.29 0.00 1.37 0.00 0.00 0.00 1.14 0.00

0.05 0.04 0.06 0.04 0.30 0.22 0.37 0.45 0.55 0.46 1.47 1.23 2.01 2.21 4.65 9.66 6.29 5.76 16.92 20.15 0.03 0.03 0.13 0.07 0.10 0.10 0.30 0.20 0.79 0.26 0.51 0.56 0.74 0.50 1.15 0.84

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.36 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.61 0.62 0.59 0.53 4.24 3.77 4.57 4.73 8.51 8.26 10.85 9.98 37.71 36.46 47.73 52.83 92.83 88.72 119.93 120.59 0.42 0.49 0.59 0.63 1.32 1.37 1.73 1.70 4.80 3.99 4.76 4.89 8.39 9.99 7.40 7.52

hard problems that test the eciency of the heuristic. Additionally, it would appear that CAB problems of type ncp are also dicult to solve. During its search for solutions, our heuristic encounters and tries various combinations of hubs. Since we know the optimal sets of hubs for all of the test problems, we tried to see if our heuristic found this optimal set during its search. Indeed, we observed that our heuristic obtained the optimal set of hubs for all problems in Table 3. However, in some cases, the heuristic is not able to

625

retain the optimal hubs because the rerouting procedure often produces poor quality solutions. For example, for problems 50TT the optimal set of hubs is encountered, but rejected since the rerouting procedure (operated on the optimal set of hubs) yields a solution that is much worse than the ®ve best solutions that at that point in time. The rerouting procedure, when applied to the optimal set of hubs for 50TT produces a solution that is 15.91% away from the optimal solution. Taking all of the above issues into consideration, we conclude that it is useful and necessary to adopt the multicommodity ¯ow procedure at the end of the heuristic. It is, however, possible for us to consider the application of the multicommodity ¯ow approach on more than (merely) the ®ve best sets of hubs. From Table 6 we see that the heuristic produces a gap of 1.13% for problem 10d2 with a ˆ 0:6. If, however, we were to consider the application of the multicommodity ¯ow approach to ®fteen sets of hubs for this problem, we would ®nd the optimal solution. This would, of course, come at an additional computational overhead. After some empirical testing, we found that, by considering ®ve sets of hubs, we are able to trade o€ solution quality and computational eciency. 4.3. E€ectiveness of the formulations The e€ectiveness of the three formulations CMAHLP-C, CMAHLP-F and CMAHLP-N are compared in Table 4 using only a subset of the full range of test data that is available. In Table 4, Gap_l is the lower gap which is de®ned as the ratio of the relaxed LP solution and the optimal solution, expressed as a percentage. TN represents the number of nodes that were evaluated in the branch and bound tree. CPU is the number of cpu seconds required to solve the problem. I_R is the number of simplex iterations that were performed at the root node and I/TN is the average number of iterations per tree node that were performed during the branch and bound search. From Table 4, we observe that CMAHLP-C, the formulation of CMAHLP provided by Campbell [10], has the smallest lower gap (or, in some cases,

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2

0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4

0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8

ÿ ÿ ÿ ÿ ÿ ÿ ÿ ÿ

10c2 10d2 10c4 10d4 15c2 15d2 15c4 15d4

10c2 10d2 10c4 10d4 15c2 15d2 15c4 15d4

10c2 10d2 10c4 10d4 15c2 15d2 15c4 15d4

10LL 10TL 10LT 10TT 20LL 20TL 20LT 20TT

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 2.07 0.00 0.00 0.00 0.00 0.00 0.00

0.00 1.32 0.00 0.00 0.00 0.00 0.00 0.00

Gap_u

Problem a I_R

I/TN

0.00 0.00 0.00 0.00 0.00 0.00 0.80 1.32

0.08 0.22 0.62 0.00 3.67 0.00 2.58 0.17

0.86 0.26 0.76 0.00 3.55 0.00 2.37 0.00

0.85 0.00 0.88 0.00 3.98 0.00 2.27 0.00

0. 0 0 0 0 0 8 16

2 8 8 0 99 0 73 4

4 4 9 0 79 0 61 2

6 3 10 0 94 0 56 0

9.94 16.15 7.94 10.82 2994.58 3177.31 5680.47 7993.72

39.68 26.24 33.77 18.18 4483.33 746.05 2969.90 521.61

77.44 24.30 36.05 18.05 5093.08 889.29 3007.00 607.57

70.82 44.82 42.89 14.47 4955.70 825.81 3066.85 621.72

34 71.6 129.4 ÿ 659.1 ÿ 601.9 302.3

506.0 54.3 154.6 ÿ 955.0 ÿ 684.8 76.5

328.7 67.7 153.8 ÿ 829.4 ÿ 765.8 ÿ

4293 ÿ 73.18 ÿ 3575 ÿ 4794 ÿ 83408 ÿ 87835 ÿ 110480 494.6 108061 765.4

12148 6925 6265 5688 47909 45921 35990 26761

8270 8244 6909 6137 42497 43756 34298 40051

8998 13027 8319 5448 42287 49213 31628 36272

2.72 3.47 4.05 2.97 1.14 2.60 4.10 2.89

1.80 5.33 2.62 2.82 5.82 2.84 5.10 4.86

1.50 2.44 1.60 1.01 4.72 0.00 3.11 2.39

1.11 0.73 1.28 0.61 4.60 0.00 2.66 1.05

13 11 21 7 4 10 28 27

6 28 28 20 130 6 93 32

6 10 18 10 104 0 79 7

8 3 12 4 95 0 71 4

TN

1.62 2.17 3.30 1.40 70.48 146.63 97.38 107.92

5.57 4.90 5.08 2.08 225.33 23.54 68.30 20.24

5.42 2.83 3.47 0.85 179.99 14.51 79.38 13.72

3.77 2.13 2.06 0.58 134.25 16.40 64.08 8.85

CPU (s)

Formulation CMAHLP-N

CPU (s)

Gap_l (%)

TN

Formulation CMAHLP-C

Gap_l (%)

Table 4 Comparison of formulations

829 1548 1122 1129 17847 29199 16961 18903

3234 1183 1388 889 9577 11200 6154 4936

3106 1919 1451 743 12541 10471 5918 6185

2555 2482 1316 711 11707 12044 6390 5093

I_R

55.3 50.7 68.2 49.7 563.3 633.6 233.1 257.3

213.0 74.9 65.2 34.0 234.8 235.8 116.2 82.8

204.0 67.6 73.2 17.9 238.1 ÿ 159.2 130.3

122.0 95.0 65.3 16.5 257.8 ÿ 184.3 137.5

I/TN

66.80 62.87 45.38 41.76 60.95 54.10 42.62 32.14

18.24 28.60 9.07 14.66 21.33 29.50 13.39 19.59

22.19 44.09 13.44 28.04 26.49 42.14 18.62 31.46

23.0 49.74 14.5 33.1 28.29 47.98 20.46 37.23

368 266 385 178 3694 1775 1829 647

39 75 41 56 190 78 141 121

41 100 57 98 210 66 137 118

50 91 64 85 218 71 155 152

TN

6.37 5.22 7.15 3.70 886.00 449.60 558.37 217.97

2.25 3.30 1.82 1.66 59.07 13.99 28.17 12.49

2.27 3.65 2.48 2.56 74.45 12.89 31.60 17.32

2.64 2.94 2.77 2.08 69.73 16.81 34.31 20.90

CPU (s)

Formulation CMAHLP-F Gap_l (%)

121 134 130 155 211 330 249 379

130 150 125 132 182 301 166 205

129 136 120 121 170 244 168 205

127 140 120 126 176 249 174 202

I_R

13.8 15.5 14.2 16.0 40.4 41.9 43.8 45.2

31.5 24.6 22.8 17.6 56.9 40.9 41.2 24.6

30.0 23.2 23.3 15.2 61.4 39.7 46.3 32.4

31.6 22.1 26.0 16.6 65.2 56.0 48.3 33.8

I/TN

626 J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

no gap at all) over all the test problems. This is because it is a stronger formulation when compared to CMAHLP-N and CMAHLP-F. CMAHLP-C has O…n4 † variables, while the other two formulations have O…n3 † variables. This increased tightness comes at the cost of much higher computational times for solving CMAHLP-C. Using this formulation, we are only able to solve problems of size n ˆ 10 and 20 for the AP data set and problems of size n ˆ 10 and 15 with a ˆ 0:2; 0:8 for the CAB data set were solved in a reasonable amount of computing time. In general, it takes more than an order of magnitude of computational e€ort to solve CMAHLP-C than it does for solving either CMAHLP-N or CMAHLP-F, even though, in general, it requires more (or, in some cases, equal number of) tree nodes to be searched to solve the latter two formulations. So, even with the limited amount of testing that we performed using CMAHLP-C, we can easily conclude that it is not as e€ective as our two new formulations. As a result we will no longer consider CMAHLP-C in the subsequent analysis and discussions. When we examine the results for the CAB problems in Table 4 above it is much harder to determine whether formulation CMAHLP-N or CMAHLP-F is better with a ˆ 0:2. For example problem 15c2 with a ˆ 0:2 is solved in approximately half the time using CMAHLP-F as opposed

627

to CMAHLP-N. However the converse is true for problem 15d4. When a ˆ 0:8 it is clear that CMAHLP-F is superior to CMAHLP-N in terms of the computational e€ort required. So, we need to solve the remaining CAB problems using both CMAHLP-N and CMAHLP-F before arriving at any conclusions on which of the two is more suited for solving these types of problems. The computational e€ort required for solving AP problems using the formulation CMAHLP-N is considerably less than that required for solving CMAHLP-F although the latter requires far fewer iterations to solve the LP relaxation at the root node (denoted by I_R). Due to lack of memory arising from the large size of the tree search, it was not possible to solve some of the larger AP problems using CMAHLP-F. Thus, in Table 5, which presents results for the larger AP problems, we only test formulation CMAHLP-N. 4.4. Solving the AP and CAB data sets In Table 5 we present results for the large AP problems. As discussed in Section 4.3, we use only CMAHLP-N for solving the larger AP problems. Tables 6 and 7 present the results for the large number of problems present in the CAB data set. We note that the total CPU seconds reported in Tables 5±7 do not include the CPU time consumed

Table 5 Results for larger AP problems Problem

25LL 25TL 25LT 25TT 40LL 40TL 40LT 40TT 50LL 50TL 50LT 50TT

Formulation CMAHLP-N Gap_ l (%)

TN

CPU (s)

I_R

I/TN

1.31 2.89 3.29 4.33 2.00 2.82 1.84 2.43 1.43 4.23 3.57 2.68

14 6 50 25 24 8 8 16 10 + 54 62

457.96 805.75 798.70 866.99 11687.75 39722.20 8966.49 5278.76 60640.03 + 79762.58 68887.91

64226 93135 46890 59595 544687 720614 467265 190662 1022345 1322307 1000988 525840

475.6 2411.8 466.9 904.3 2143.9 19168.6 2342.8 2306.4 14580.2 + 5672.2 4831.0

628

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

by the heuristic (as presented in Table 3). The upper gaps for all the AP problems and a few of the CAB problems have already been provided in Table 3. So, we only provide the heuristic upper bounds for the CAB problems. It appears that the AP problems of type TL are the most dicult problems in this data set. We are unable to solve problem 50TL in under 300000 CPU seconds (hence a `+' sign in the CPU secs entry against this problem). Although we do not obtain a solution for 50TL using formulation CMAHLP-N, we have, however, managed to solve this problem using

the formulation CMAHLP-F in 38885 CPU seconds. Indeed, this is how we are able to quote an optimal solution for Problem 50TL in Table 1. The formulation CMAHLP-F performs better than formulation CMAHLP-N for only problems 40TL and 50TL in the AP data set. We note that the time for solving Problem 40TL using CMAHLP-F is 10174 CPU secs. However, for some of the other AP problems CMAHLP-F is unable to obtain solutions because the branch and bound tree grows to an unmanagable size, thus exceeding the speci®ed memory allocation for storing the tree.

Table 6 Results for CAB problems with a ˆ 0:2; 0:4; 0:6 Problem

a

Gap_u Formulation CMAHLP-N (%)

Gap_l (%)

TN

Formulation CMAHLP-F

CPU (s)

I_R

I/TN

Gap_l (%)

TN

CPU (s)

I_R

I/TN

20c2 20d2 20c4 20d4 25c2 25d2 25c4 25d4

0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

3.57 0.38 3.07 0.31 1.82 0.88 1.77 1.92

70 3 78 2 20 8 37 4

836.25 114.31 194.20 44.56 1511.60 674.78 595.10 266.54

34322 37649 19944 20535 117263 77388 60284 58079

884.0 853.7 271.4 217.5 1963 2021.3 474.9 986.5

34.11 53.72 27.08 42.48 31.47 52.55 26.39 43.78

325 152 551 188 362 280 642 418

404.84 98.47 442.68 75.40 1322.88 493.77 1658.91 502.31

461 613 378 436 3887 4265 2984 3764

106.0 82.3 71.7 48.3 141.7 98.9 103.3 80.61

20c2 20d2 20c4 20d4 25c2 25d2 25c4 25d4

0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4

5.89 0.00 0.00 0.00 0.00 0.00 0.00 0.00

4.70 2.27 3.36 2.03 2.36 1.91 2.18 3.00

103 20 84 7 21 34 34 5

1138.70 442.57 236.34 44.75 2013.97 3348.62 703.88 556.12

32237 59031 18878 15289 125527 141472 66818 68068

659.0 761.8 249.5 134.4 1936.9 1896.9 479.3 1445.6

31.96 48.12 24.98 35.99 29.15 47.41 24.09 37.56

378 198 426 160 266 302 485 398

598.75 148.28 366.72 71.62 913.71 460.35 1032.83 389.95

421 536 369 420 3876 3967 2647 2983

127.4 73.2 71.1 48.1 130.0 92.3 96.0 66.5

10c2 10d2 10c4 10d4 15c2 15d2 15c4 15d4 20c2 20d2 20c4 20d4 25c2 25d2 25c4 25d4

0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6

0.00 1.13 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

1.66 3.89 2.50 0.97 5.17 0.94 3.92 3.48 5.76 5.69 3.68 3.30 3.43 3.56 2.85 4.30

6 12 22 9 125 4 74 12 115 30 81 10 22 40 45 7

5.53 2.82 4.30 1.15 207.17 16.65 72.27 10.40 979.66 425.02 230.57 42.70 1699.90 2431.03 672.42 519.96

3017 1761 1611 1045 11862 11823 6026 4536 41191 55643 16629 15088 113156 146956 43316 64313

195.3 73.2 71.0 25.7 248.0 111.0 157.8 91.2 623.8 778.7 210.4 147.6 1671.2 1209.3 373.3 1324.7

20.70 37.02 11.82 21.45 24.26 36.00 16.45 25.50 29.02 42.55 22.01 29.40 26.31 41.56 21.40 31.12

42 91 49 73 209 74 124 124 291 210 278 128 234 289 366 316

2.15 3.10 2.15 1.84 54.52 11.55 29.42 13.77 258.03 119.38 202.47 48.65 677.13 613.21 748.40 370.06

125 154 117 122 180 262 161 201 447 549 396 489 3755 4425 2534 3205

27.4 21.4 24.4 15.6 54.5 36.2 45.3 26.4 92.3 68.8 63.3 46.2 126.7 89.7 95.1 64.7

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

629

Table 7 Results for CAB problems with a ˆ 0:8; 1:0 Problem

a

Gap_u

Formulation CMAHLP-N

(%)

Gap_l (%)

TN

CPU (s)

I_R

I/TN

Gap_l (%)

Formulation CMAHLP-F TN

CPU (s)

I_R

I/TN

20c2 20d2 20c4 20d4 25c2 25d2 25c4 25d4

0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8

0.00 2.07 0.00 0.00 0.00 1.42 13.85 0.00

7.00 8.71 4.58 5.36 4.42 4.88 3.97 6.20

134 63 74 17 36 32 67 16

1162.57 565.88 199.72 76.14 2237.18 2837.60 936.94 457.97

36369 33963 16951 13397 104283 122515 35636 48472

533.5 501.1 176.1 228.0 1287.6 1807.6 415.6 571.4

25.19 36.01 17.97 22.20 23.10 35.11 28.20 24.50

290 307 169 122 189 256 308 224

285.48 196.47 108.27 54.07 578.29 437.55 1156.65 214.28

505 568 404 526 3853 4324 2846 3491

90.4 67.4 58.2 48.8 139.0 102.6 155.8 64.3

10c2 10d2 10c4 10d4 15c2 15d2 15c4 15d4 20c2 20d2 20c4 20d4 25c2 25d2 25c4 25d4

1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0

0.00 2.55 0.00 0.00 0.00 0.00 0.00 0.00 0.76 0.00 0.00 0.00 0.00 2.28 0.00 0.00

2.88 7.02 4.02 5.90 5.74 5.45 6.41 8.46 7.63 10.13 6.01 8.59 5.99 8.88 5.64 7.92

14 47 26 52 122 16 175 129 131 78 107 62 66 68 153 39

4.37 5.68 3.67 3.25 156.25 31.10 106.96 34.83 1098.63 526.97 281.27 100.21 2652.05 3093.47 1250.66 259.75

2371 1023 965 682 8949 7584 3794 2939 29156 24325 12243 9456 88812 104198 28357 26917

70.9 56.0 50.3 27.1 194.4 205.4 94.5 51.5 558.4 453.4 177.5 138.0 1058.1 1090.6 240.9 207.2

15.67 20.35 5.66 8.49 17.38 23.90 10.13 14.63 20.38 28.36 13.71 16.01 20.12 29.06 14.77 18.34

25 75 29 56 151 108 171 151 219 305 164 130 220 249 240 121

1.52 2.30 1.19 1.60 31.64 14.60 28.45 28.14 146.35 134.17 107.32 42.87 438.20 386.50 436.55 113.91

244 229 198 191 365 453 307 315 1123 1011 656 719 4642 4805 3180 3314

38.8 22.3 23.5 19.8 57.7 39.5 42.1 25.9 93.6 63.6 76.1 50.1 130.8 110.4 106.1 78.3

We note that the heuristic either ®nds the optimal solution or tight upper gaps (represented by Gap_u) for most of the problems. There are some anomalies in which the heuristic yields bad results. One of them is Problem 20c2 (a ˆ 0:4) with a upper gap of 5.89%. Another one is Problem 25c4 (a ˆ 0:8) in which the upper gap is as high as 13.85%. The reason for this is that the heuristic gets caught in a locally optimal neighbourhood. In our heuristic, we only perform a 1-interchange swap procedure. In other words, one potential hub which is not in the trial set is swapped for a hub that is in the trial set. In these two problems, each of the individual hubs in the optimal solution has limited capacity. However, together, they are able to satisfy the capacity restrictions. Thus, in performing a 1-interchange procedure, the optimal hubs are often rejected from consideration in the trial set.

We found that, if we had implemented a 2-interchange procedure in our heuristic, we are able to get the optimal hubs in these problems at a much larger computational e€ort (since 2-interchange is computationally more expensive than the 1-interchange procedure that is found to be adequate, except in these two cases). In general, the lower gaps (denoted by Gap_l) are reasonably small over all the AP problems even as n increases. Because of the relative tightness of the lower gap, we need only a few tree nodes (in most cases) to solve the problem. In Table 6, we present results for the CAB problems with a ˆ f0:2; 0:4; 0:6g. Table 7 provides results for the CAB problems with a ˆ f0:8; 1:0g. In Table 6, we do not repeat the results that have already been provided to the CAB problems (with n ˆ 10; 15 and a ˆ 0:2; 0:4; 0:8). We compare formulations CMAHLP-N and CMAHLP-F by examining Tables 6 and 7. The lower gap at the root

630

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

node for CMAHLP-F is signi®cantly large in comparison to the relatively tight lower bounds produced at the root node when using CMAHLP-N. However, since the computational e€ort required by CMAHLP-F at each tree node is relatively small (as indicated by I/TN), we are able to bridge the large lower gap quickly. We observe that as a increases, CMAHLP-F becomes considerably more e€ective. This can probably be attributed to the fact as a increases the lower gap for CMAHLP-F decreases while the gap for CMAHLP-N increases. With a ˆ 0:2 the average gap of CMAHLP-F is approximately 22 times that of CMAHLP-N, where on the other hand when a ˆ 1:0 the average gap of CMAHLP-F is only about 3 times that of CMAHLP-N. In the case of a ˆ 0:8; 1:0, every problem type is solved more quickly when using CMAHLP-F rather than CMAHLP-N. Even for smaller values of a, it appears as though it is better to use CMAHLP-F for the CAB problems in which the value of p is small. For example, regardless of the value of a, CMAHLP-F can solve problems 20c2, 20d2, 25c2, 25d2 (in which p ˆ 2) quicker than CMAHLP-N. Thus, we can conclude that, given our test examples, we would be better o€ using CMAHLP-F for problems in which p is small and for problems with large values of a.

quickly using our formulations than using existing formulations. The existing formulation of CMAHLP in the literature is unable to solve larger hub location problems in a reasonable amount of computational time. However, due to the smaller size of our new formulations, we are able to solve moderately-sized capacitated multiple allocation hub location problems. We also tested the e€ectiveness of the two new formulations that we presented and provided pointers to the types of problem sets that each of them will be suitable for. Real world problems are, however, much larger than those presented and solved in this paper. Further research must concentrate on increasing the size of hub location problems that can be solved optimally. This could be achieved either by adopting better approaches or by developing a new formulation that incorporates the best features of the two successful formulations that are described in this paper.

5. Conclusions

References

In this paper we described a new MILP for the capacitated multiple allocation hub location problem (CMAHLP) which had fewer variables and constraints than those previously presented in the literature. We also develop a variation of this new formulation of CMAHLP. We also described an ef®cient constructive heuristic algorithm for solving the CMAHLP. This uses a shortest path approach and yields tight upper bounds for most of the problems that we tested. We incorporate the upper bound obtained from this heuristic in a linearprogramming-based branch-and-bound solution procedure. Our new formulations are weaker than the strong MILP formulations presented in the literature. Yet, we are able to solve problems more

Acknowledgements We acknowledge the helpful and constructive comments provided by the two anonymous referees.

[1] R.K. Ahuja, T.L. Magnanti, J.B. Orlin, Network Flows: Theory Algorithms and Applications, Prentice Hall, Englewood Cli€s, NJ, 1993. [2] T. Aykin, Networking policies for hub-and-spoke systems with applications to the air transportation system, Working paper, Rutgers Graduate School of Management, Newark, NJ, 1993. [3] T. Aykin, Lagrangian relaxation based approaches to capacitated hub-and-spoke network design problem, European Journal of Operational Research 79 (1994) 501± 523. [4] T. Aykin, Networking policies for hub-and-spoke systems with applications to the air transportation system, Transportation Science 29 (3) (1995) 201±221. [5] T. Aykin, G. Brown, Interacting new facilities and location-allocation problems, Transportation Science 26 (1992) 212±222. [6] N. Barton, P. Cleary, J. Ha, M. Krishnamoorthy, J. Mooney, N. Stokes, Industrial case studies at the

J. Ebery et al. / European Journal of Operational Research 120 (2000) 614±631

[7] [8] [9] [10] [11] [12] [13]

[14]

[15] [16] [17]

interface between engineering and applied mathematics, in: Proceedings of the 1st Biennial Australian Engineering Mathematics Conference, Melbourne, 1994. J.E. Beasley, Or-Library: Distributing test problems by electronic mail, Journal of the Operational Research Society 41 (11) (1990) 1069±1072. R. Boorstyn, H. Frank, Large-scale network topological optimisation, IEEE Transactions on Communications 25 (1977) 29±47. J.F. Campbell, Location-allocation for distribution systems with transshipments and transportation economies of scale, Annals of Operations Research 40 (1992) 77±99. J.F. Campbell, Integer programming formulations of discrete hub location problems, European Journal of Operational Research 72 (1994) 387±405. J.F. Campbell, Hub location and the p-hub median problem, Operations Research 44 (6) (1996) 923±935. A.T. Ernst, M. Krishnamoorthy, Ecient algorithms for the uncapacitated single allocation p-hub median problem, Location Science 4 (1996) 139±154. A.T. Ernst, M. Krishnamoorthy, Exact and heuristic algorithms for the uncapacitated multiple allocation phub median problem, European Journal of Operational Research 104 (1) (1998) 100±112. A.T. Ernst, M. Krishnamoorthy, An exact solution approach based on shortest-paths for p-hub median problems, INFORMS Journal on Computing 10 (2) (1998) 149±162. A.T. Ernst, M. Krishnamoorthy, Solution algorithms for the capacitated single allocation hub location problem, Annals of Operations Research, forthcoming. B. Gavish, Models for con®guring large-scale distributed computing systems, AT & T Technical Journal 64 (2) (1985) 491±531. S.L. Hakimi, Optimum locations of switching centres and the absolute centres and medians of a graph, Operations Research 12 (1964) 450±459.

631

[18] S.L. Hakimi, Optimum distribution of switching centres in a communication network and some related graph theoretic problems, Operations Research 13 (1965) 462±475. [19] CPLEX Optimization, Using the CPLEX callable library (version 4.0), CPLEX Optimization, Suite 279, 930 Tahoe Blvd., Bldg 802, Incline Valley, NV 89451±9436, USA, 1996. [20] J.G. Klincewicz, Heuristics for the p-hub location problem, European Journal of Operational Research 53 (1) (1991) 25±37. [21] J.G. Klincewicz, Dual Algorithms for the uncapacitaed hub location problem, Location Science 4 (1996) 173±184. [22] M. Krishnamoorthy, G. Mills, D. Sier, Strategic con®guration of the mail processing network: Location-allocation modelling ± Stage-1, CSIRO Technical Report DMS-C94/ 9 Division of Mathematics and Statistics CSIRO Australia, 1994. [23] M. O'Kelly, Activity levels at hub facilities in interacting networks, Geographical Analysis 18 (1986) 343±356. [24] M. O'Kelly, The location of interacting hub facilities, Transportation Science 20 (1986) 92±106. [25] M. O'Kelly, A quadratic integer program for the location of interacting hub facilities, European Journal of Operational Research 32 (1987) 393±404. [26] M. O'Kelly, Hub facility location with ®xed costs, Papers in Regional Science: The Journal of the RSAI 71 (1992) 293±306. [27] D. Skorin-Kapov, J. Skorin-Kapov, On tabu search for the location of interacting hub facilities, European Journal of Operational Research 73 (1994) 502±509. [28] D. Skorin-Kapov, J. Skorin-Kapov, M. O'Kelly, Tight linear programming relaxations of uncapacitated p-hub median problems, European Journal of Operational Research 94 (3) (1996) 582±593. [29] K. Smith, M. Krishnamoorthy, M. Palaniswami, Neural versus traditional approaches to the location of interacting hub facilities, Location Science 4 (1996) 155±171.