A Benders decomposition approach for a distribution network design problem with consolidation and capacity considerations

Operations Research Letters 39 (2011) 138–143

Contents lists available at ScienceDirect

Operations Research Letters journal homepage: www.elsevier.com/locate/orl

A Benders decomposition approach for a distribution network design problem with consolidation and capacity considerations Halit Üster a,∗ , Homarjun Agrahari b a

Department of Industrial and Systems Engineering, Texas A&M University, College Station, TX, 77843-3131, United States

b

Division of Service Design and Performance, BNSF Railway, Fort Worth, TX, 76131-1322, United States

article

info

Article history: Received 8 November 2008 Accepted 5 January 2011 Available online 23 February 2011 Keywords: Network design Consolidation Benders decomposition

abstract We develop a model for a strategic freight-forwarding network design problem in which the design decisions involve the locations and capacities of consolidation and deconsolidation centers, and capacities on linehaul linkages as well as the shipment routes from origins to destinations through centers. We devise a solution approach based on Benders decomposition and conduct a computational study that illustrates the efficiency and the effectiveness of the approach. © 2011 Elsevier B.V. All rights reserved.

1. Introduction We consider an integrated distribution network design and site selection problem arising in the context of transportation planning faced by the freight-forwarding industry and large corporations operating over large geographical regions. A freight-forwarder (a third-party logistics provider) acts as an intermediary between a shipper and a carrier (transportation service provider) and is responsible for planning of transportation of the shipper’s freight which, in turn, is handled by the carrier [7, pp. 353]. Clearly, a freight-forwarding company serves many shippers, handling freight with many distinct origins and destinations. By planning for cost effective transportation with the best rates for the shippers via efficient use of cargo space, freightforwarders achieve profitability via transportation economies-ofscale. In the case of large corporations, such as high volume/variety manufacturers or suppliers operating many facilities, the products (e.g. spare parts, replacement parts, returns, and finished products) are frequently shipped among the facilities. For example, assuming that the facilities include four levels: plants, warehouses, stores and customers, the shipments occur between every pair of these locations – back and forth – as well as within each level except inter-customer shipments. Similar to the freight-forwarders, these corporations, as the bearers of the shipment costs, strive for cost effectiveness through transportation economies-of-scale. As a result, efficient load consolidation arises as an important opportunity for profitability in these applications.

∗

Corresponding author. Tel.: +1 979 845 9573. E-mail address: [email protected] (H. Üster).

0167-6377/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.orl.2011.02.003

Motivated by these applications, we consider strategic level multi-commodity network design where each commodity is defined by a unique pair of origin and destination points and a known required flow amount (load). To achieve economies-ofscale, the transportation activities take place in such a way that a commodity, on its trip from origin to destination, is first shipped to a consolidation center where it is combined with several other commodities into a truckload (TL) for a linehaul transfer to a deconsolidation center from which the commodities are shipped to their final destinations. To undertake this type of operation, the locations of the centers with consolidation and/or deconsolidation capabilities are to be determined via selection from a set of candidates. The candidate locations with these capabilities may correspond to existing sites that can be fully or partially dedicated to centers, to warehouse spaces that can be rented or to new sites on which centers are to be established. In either case, specific resources include sorting, loading, and unloading equipment, temporary storage space, and human resources. Thus, along with the location of the centers, it is important to take into account available capacities for specific activities and their associated costs when making strategic decisions on the design of the network. Similarly, the capacity decisions are to be made for linehaul links. For both centers and linehaul links, we consider capacity installments determined by fixed increments of a unit capacity which is conveniently represented by TL capacities. In summary, we explicitly consider the following decisions while designing the distribution network with the above operational characteristics: 1. the locations and capacities of the consolidation and deconsolidation centers;

H. Üster, H. Agrahari / Operations Research Letters 39 (2011) 138–143

2. the selection of linehaul transfer links (between consolidation and deconsolidation centers) and their capacities in terms of the number of truckload trips, and 3. the routing of commodities from their origins to destinations through centers. As a result, our objective pertains to minimizing the overall design cost, which consists of LTL and TL type transportation costs and center location/selection costs. Since the origin-toconsolidation-center and deconsolidation-center-to-destination shipments involve smaller loads (individual commodities), it is reasonable to assume that these shipments are in LTL (LessThan-Truckload) mode in which, typically, a constant dollar value per weight per mile of shipment is used in calculating transportation costs. On the other hand, since the linehaul transfers (consolidation-to-deconsolidation) involve full truckload shipments, we assume a TL mode of transportation for this component and calculate the associated costs by employing per mile shipment costs for dispatching a TL. We assume that the commodities with large loads and/or the ones whose origin and destination nodes are in close proximity can be shipped directly without consolidation efforts. To the best of our knowledge, the problem we pose in this study is not addressed in the previous literature, however, the areas of hub location (HL) as well as fixed charge network design (FCND) are closely related. HL problems are concerned with locating hubs and assigning the nodes of a physical network to these hubs so that the total cost of fixed hub locations and transportation is minimized over the network. The economies-of-scale for inter-hub transfers provide the motivation for locating hubs. Alternative formulations and algorithmic studies are found for the capacitated case (where a hub has a capacity limitation for the flow it can handle) in [8,12, 18,20] and, for the uncapacitated case, in [1,2,10,11]. An excellent review of the HL problems can be found in [9]. Several differences, that also underline the generalizations we address, between HL problems and our study can be outlined as follows. First, we consider explicit commodity-based routing decisions, which is not the case in hub location problems. More specifically, we are interested in assigning commodities to consolidation and deconsolidation centers as opposed to HL problems in which the assignment of a node to a hub implicitly determines the assignment of the commodities destined-to and originating-from that particular node to that hub. Second, in HL, the cost structures for from-node-to-hub, from-hub-to-node, and inter-hub transfer are typically assumed to be of LTL. To reflect the economies-of-scale, an estimate of the per-unit per-mile cost that is significantly lower (discounted) is employed on the inter-hub links (compared to other connections). In our case, we consider truckload shipments on linehaul links and adopt an associated cost structure that explicitly captures consolidation practice of interest. Third, we explicitly consider capacity decisions under fixed capacity installment costs on transfer links as well as at the centers. Therefore, despite operational similarities, the problem of interest has fundamental differences from the HL problem. In FCND problems, given an underlying network, a set of commodities with specific origins and destinations and fixed and variable arc costs, the focus is on determining the routing of commodities that minimizes the total fixed and variable costs on the arcs. Depending on the application, the arcs can be uncapacitated [13,16] or capacitated [14,15,4]. Studies considering piecewise linear concave costs include [3], which considers a Lagrangian approach, and [19], which provides a linearprogramming based heuristic approach. The problem at hand distinguishes itself from the general FCND problems in several ways. First, we consider the consolidation and deconsolidation activities as well as locational decisions explicitly in our model. Second, linehaul transfer costs are not linear in

139

flow, as is typically considered in FCND, with some exceptions such as [14,15], but rather, are represented by a step function based on consolidated TL shipments. Third, our problem has a structure that requires the flow of commodities to use a path with at most three hops (i.e., three arcs including origin-to-center, linehaul, center-to-destination). Clearly, if center locations and physical nodes corresponding to commodity origins coincide in the optimal solution to our problem, a commodity’s path from originto-destination may include a lower number of hops. In the rest of the paper, we first provide a compact model formulation (Section 2). In doing so, we take the perspective of a freight-forwarder or a corporation’s own unit for logistics planning where both utilize a common carrier for service. Our formulation facilitates the use of a Benders decomposition framework for its solution. For this purpose, we also generate strong Benders cuts that greatly improve the performance of the Benders algorithm (Section 3) as we illustrate in our computational studies (Section 4). We summarize our results and conclusions as well as future research directions in Section 5. 2. Problem formulation In our formulation, we utilize five sets with the following formations. The set P = {1, . . . , N } represents the set of N commodities. The origin and destination nodes of the commodities are given by F ={f1 , . . . , fN } and T ={t1 , . . . , tN }, respectively. In addition, the sets J and K represent the candidate center locations where the consolidation and deconsolidation activities can take place, respectively. The rest of the notation is as follows: Parameters:

wi the demand for commodity i ∈ P U S Cj Dk

β α

dfi j dkti djk

capacity per TL for linehaul transfers base capacity (in TLs) for consolidation and deconsolidation centers cost of base capacity for consolidation at center j ∈ J cost of base capacity for deconsolidation at center k ∈ K TL transportation cost per mile (for linehaul transfers) LTL transportation cost per unit per mile (for collection and distribution) distance between fi , i ∈ P and consolidation center j ∈ J distance between deconsolidation center k ∈ K and ti , i∈P distance between centers j ∈ J and k ∈ K .

Decision Variables: zijk fraction of commodity i’s demand assigned to link (j, k), j ∈ J, k ∈ K yjk capacity in number of TLs installed on the link (j, k), j ∈ J, k ∈ K vjc number of base capacity units for consolidation at center j∈J vkd number of base capacity units for deconsolidation at center k ∈ K. For the sake of brevity in formulation, we also let Tjk and Wijk represent β djk and wi α (dfi j + dkti ), respectively. Then, the problem of interest can be formulated as follows:

(MnP) Min

−−−

Wijk zijk +

i∈P j∈J k∈K

+

−

Cj vjc +

j∈J

−

Dk vkd

−−

Tjk yjk

j∈J k∈K

(1)

k∈K

subject to

−− j∈J k∈K

zijk = 1 ∀ i ∈ P ,

(2)

140

H. Üster, H. Agrahari / Operations Research Letters 39 (2011) 138–143

−

wi zijk ≤ U yjk ∀ j ∈ J , k ∈ K ,

(3)

Then, for fixed values of y¯ , v¯ c , and v¯ d , the Benders subproblem BSP is given as the following linear program.

yjk ≤ S vjc

(4)

(BSP) Min

i∈P

−

∀ j ∈ J,

k∈K

−

−−−

Wijk zijk

(7)

i∈P j∈J k∈K

yjk ≤ S vkd

∀ k ∈ K,

(5)

zijk ≥ 0, yjk , vjc , vkd ∈ Z+

subject to

−−

j∈J

∀ i ∈ P , j ∈ J, k ∈ K .

(6)

The first term in the objective function (1) is the sum of the collection and distribution costs for the commodities. The last three terms represent the sum of costs of installing capacity on transfer links and for consolidation and deconsolidation at the centers, respectively. The first constraint set (2) guarantees that each required commodity flow is achieved. The constraint set (3) dictates the capacity requirements/restrictions on transfer links. The constraint sets (4) and (5) provide capacity requirements/restrictions for consolidation and deconsolidation activities at the centers, respectively. Lastly, the constraint set (6) controls the nonnegativity and integrality of the variables. 3. Benders decomposition based solution approach In general terms, the Benders decomposition technique involves decomposing the overall formulation (after a reformulation) into a master problem and a subproblem, and then solving them iteratively by utilizing the solution of one in the other [6]. The master problem involves a mixed integer program with only one continuous variable that is introduced as an auxiliary variable to facilitate the interaction between the master and the subproblem. The subproblem, on the other hand, is a linear program incorporating the integer variables as parameters whose values are determined by solving the master problem. In an iteration of the overall solution procedure, we solve the current master problem and determine a lower bound for the overall problem along with the corresponding values of the integer variables. We then solve the dual of the subproblem and obtain an upper bound by utilizing its objective value along with the cost components implied by the master problem solution. The bounds are updated and if the stopping criterion is not met, we generate a Benders cut using the dual subproblem solution. We add this cut to the master problem, and continue iterations until a stopping criterion is met. As stopping criteria, we employ a small percentage gap between the best upper and lower bounds and a maximum number of iterations reached, and whichever is achieved first terminates the procedure. To see the special structure of our formulation that facilitates a decomposition approach, notice that the integer variables y, vc , and vd model capacity related decisions and the continuous variable z models routing (shipment quantities and paths) decisions. We observe that for known capacities, i.e., for fixed y, vc , and vd values, the routing problem is a linear program that can be solved efficiently. Utilizing this structure, in our Benders decomposition based solution approach, the master problem includes the integer capacity decisions and the subproblem determines the routing of commodities under given capacity levels. 3.1. Benders subproblem and its dual In order to generate a Benders reformulation of our problem MnP, we first present the subproblem formulation. For this purpose, let y¯ , v¯ c , and v¯ d represent given values of capacity installations on links and centers that guarantee a feasible routing.

zijk = 1 ∀ i ∈ P ,

(8)

j∈J k∈K

−

wi zijk ≤ U y¯ jk ∀ j ∈ J , k ∈ K ,

(9)

i∈P

zijk ≥ 0 ∀ i ∈ P , j ∈ J , k ∈ K .

(10)

To obtain the dual of this subproblem (DBSP), we define the dual variables µi and λjk for constraints (8) and (9), respectively. Then, we have (DBSP)

−

Max

µi −

i∈P

−−

U y¯ jk λjk

(11)

j∈J k∈K

subject to

µi − wi λjk ≤ Wijk ∀ i ∈ P , j ∈ J , k ∈ K ,

(12)

µi -free, λjk ≥ 0 ∀ i ∈ P , j ∈ J , k ∈ K .

(13)

Let B denote the set of all extreme points of the DBSP polyhedron b given by (12) and (13), and µbi , δijk , λbjk and ηb , b ∈ B , denote the values of the associated dual variables and objective function, respectively. Also letting η∗ be the the optimal objective value, we must have η∗ ≥ ηb , ∀ b ∈ B , and, thus, the DBSP can be restated as minη≥0 {η : η ≥ ηb , ∀ b ∈ B }, where

ηb =

−

µbi −

i∈P

−−

U y¯ jk λbjk

∀ b ∈ B.

j∈J k∈K

3.2. Reformulation of MnP and the Benders master problem Utilizing the above representation of the DBSP that is based on the extreme points of its polyhedron, we can reformulate the MnP as Min

η+

−−

Tjk yjk +

− j∈J

j∈J k∈K

Cj vjc +

−

Dk vkd

(14)

k∈K

subject to

−

yjk ≤ S vjc

∀ j ∈ J,

(15)

yjk ≤ S vkd

∀ k ∈ K,

(16)

k∈K

− j∈J

η≥

− i∈P

µbi −

−−

U yjk λbjk

∀ b ∈ B,

(17)

j∈J k∈K

η ≥ 0, yjk , vjc , vkd ∈ Z+ ∀ i ∈ P , j ∈ J , k ∈ K .

(18)

This reformulation presents a difficulty due to a very large number of constraints of type (17). Moreover, at optimality, not all of the constraints in (17) will be binding. Therefore, in the iterative Benders decomposition approach [6], one works with a relaxed version of MnP by considering only a subset of these constraints in each iteration. We denote this subset by (BendersCutSet) which includes the constraints (17) generated via solving the DBSPs in the previous iterations. The relaxed formulation thus obtained is termed the Benders master problem, BMP as given below, and its optimum solution clearly provides a lower bound on the optimum solution to MnP. In order to obtain stronger lower bounds with improved solution times while solving the BMP we also include the surrogate

H. Üster, H. Agrahari / Operations Research Letters 39 (2011) 138–143

constraint (23) into the BMP. Although not needed for correct formulation of MnP, (23) ensures enough total capacity installment on the transfer links. Specifically, the constraint (23) is essential for the boundedness of the DBSP (i.e., for the feasibility of the BSP).

−−

(BMP) Min η +

Tjk yjk +

j∈J k∈K

− j∈J

Cj vjc +

−

Dk vkd

(19)

k∈K

yjk ≤ S vjc

∀ j ∈ J,

(20)

yjk ≤ S vkd

∀ k ∈ K,

(21)

k∈K

− j∈J

(22)

(BendersCutSet)

−− j∈J k∈K

U yjk ≥

−

wi

(23)

i∈I

η ≥ 0, yjk , vjc , vkd ∈ Z+ ∀ i ∈ P , j ∈ J , k ∈ K .

determine the optimum µ ¯ i , ∀i ∈ P and λ¯ jk , ∀(j, k) ∈ O values by solving the following linear program. Reduced (DBSP)

−

Max

i∈P

subject to

−

141

µi −

−

U y¯ jk λjk

(25)

(j,k)∈O

subject to

µi − wi λjk ≤ Wijk ∀ i ∈ P , (j, k) ∈ O ,

(26)

µi -free, λjk ≥ 0 ∀ i ∈ P , (j, k) ∈ O .

(27)

Phase 2. Therefore, fixing the µ as determined in the first phase to ¯ , secondly, we solve the following models – one for each (j, k) µ for which the y¯ jk is equal to zero – to obtain the values of such dual variables so that the resulting set of dual variables provide a strong Benders cut for inclusion into the (BendersCutSet). Max

−

−

U λjk

(28)

(j,k)∈O ′

(24)

In each iteration of the Benders approach, first a master problem is solved to obtain the values of the integer variables y, vc , and vd . Then, these values are used to solve the DBSP to obtain a set of values for dual variables µb , δb , and λb (i.e., an extreme point b ∈ B of the dual polyhedron) and a new cut in the form of (17) to include into (BendersCutSet). Note that, when the DBSP is solved for fixed y¯ , v¯ c , and v¯ d values, an upper bound for MnP can easily be calculated by adding DBSP’s objective value and the total of fixed cost components from the corresponding BMP solution (i.e. the objective value of the corresponding BMP excluding the value of η). 3.3. Solving subproblem DBSP and strengthened Benders cuts Observe that the structure of DBSP resembles the dual of the transportation problem. Specifically, treating each linehaul transfer link (j, k) with a non-zero y¯ jk value as a supply point with available supply of U y¯ jk and each commodity as a demand point with a demand of wi , BSP is a transportation problem with shipment costs Wijk . Thus, due to well-known inherent degeneracy in transportation problems, DBSP may have alternative optimal solutions leading to the possibility of multiple Benders cuts that can be generated. In such situations, it has been observed that (e.g. in contexts of capacitated facility location and network design) obtaining strengthened Benders cuts in an iteration improves the approach significantly [17,21,22]. Magnanti and Wong [17] define the concept of strongness of a cut as follows. Given an optimization problem miny∈Y , z ∈ℜ {z : f (u) + y g (u) ≤ z ∀ u ∈ U }, the cut f (u1 ) + y g (u1 ) ≤ z dominates (is stronger than) the cut f (u2 ) + y g (u2 ) ≤ z, if f (u1 )+y g (u1 ) ≥ f (u2 )+y g (u2 ) ∀ y ∈ Y with a strict inequality for at least one y ∈ Y . Then, due to the existence of multiple solutions in DBSP, a solution to DBSP should be determined so that it provides a strengthened Benders cut. For this purpose, we proceed similarly to the two-phase approach given by Van Roy [21] for the capacitated facility location problem. For our specific problem, observe that the dual values associated with y¯ jk parameters with a value of zero do not have any impact on the optimum objective value of the DBSP. That is, when y¯ jk = 0, we can modify its coefficient −U λjk without changing the objective function value, provided that the feasibility is maintained, i.e. it satisfies constraint (12). Then, to solve DBSP efficiently, we proceed as follows: Phase 1. We solve a reduced DBSP and obtain a set of dual variable values λjk associated with links (j, k) with a nonzero y¯ jk value as well as optimum µi values for all i ∈ P . To obtain the reduced DBSP, let O be the set of links (j, k) where y¯ jk ≥ 1. Then, we

subject to

µ ¯ i − wi λjk ≤ Wijk ∀ i ∈ P , (j, k) ∈ O ′

(29)

λjk ≥ 0 ∀ (j, k) ∈ O

(30)

′

where O ′ = (J × K ) \ O . The optimum solution to (28)–(30) can be found as follows. The constraint (29) can be written as λjk ≥ rijk , where rijk = (µ ¯ i − Wijk )/wi , ∀ i ∈ P , (j, k) ∈ O ′ . Letting rˆjk = maxi∈P {ri }, the optimum λjk solution to (28)–(30) is given by

 rˆ λ¯ jk = jk 0

if rˆjk > 0 otherwise

∀ (j, k) ∈ O ′ .

Employing the above specific approach – reduced DBSP in Phase 1, followed by the efficient solution approach (i.e., not relying on an LP solver) for Phase 2 – to address subproblem and generation of Benders cuts greatly enhances the solution times in the BD framework. It also provides the ability to solve large instances without having memory problems, since only the Phase 1 problem, whose size is significantly reduced, relies on an LP solver. 3.4. Benders decomposition framework Having defined the Benders dual subproblem DBSP and its solution so as to generate strengthened cuts and the master problem BMP which we solve using CPLEX, we next outline the overall algorithm given in Display 1. However, we first note that, through our computational studies, the solution time to obtain ∗ an optimal solution to BMP, ZBMP , which is a lower bound on the original problem, increases quickly as the number of cuts is increased due to growth in the problem size. This clearly affects the overall (solution time) performance of the algorithm in an adverse manner. Thus, whenever we solve the BMP, we allow CPLEX to continue with branching until an optimality gap (CPLEX parameter EpGap) of g% is reached for the Benders iterations more than mitr. In our computational study, (Section 4), we employ g values of 5, 4, 3, 2, and 0.1 for mitr values of 0, 40, 60, 80, and 200, respectively. Upon termination, we record the lower and upper bounds to the BMP as Z LBMP and Z UBMP . Since Z LBMP is a lower bound on Z ∗BMP , it is also a lower bound on the original problem, i.e., it is the lower bound obtained at an iteration itr, denoted by LBitr . On the other hand, to calculate a valid upper bound UBitr at an iteration itr, we need to utilize the Z UBMP , which is a feasible solution to the BMP. Specifically, we add the fixed costs obtained in this feasible solution (Z UBMP − η¯ ) and the optimum routing cost in the corresponding routing problem (Z DBSP ), thus obtaining the objective value of a feasible solution to the original problem.

142

H. Üster, H. Agrahari / Operations Research Letters 39 (2011) 138–143

Display 1 Benders Decomposition Based Algorithm 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:

12: 13: 14: 15: 16: 17: 18:

initialize LB = −∞, UB = ∞, ϵ = 0.02, gap = 2ϵ , itr=0, Maxitr=400, BendersCutSet= ∅ Solve BMP with early termination criteria to obtain ZL and Z U with its associated (¯y, v¯ c , v¯ d , η) ¯ BMP BMP while (itr ≤ Maxitr and gap > ϵ ) do itr = itr + 1 Substituting y¯ , solve DBSP to obtain Z DBSP and (µ, λ) Obtain UBitr as Z U − η¯ + Z DBSP BMP itr if UB < UB then UB = UBitr end if Add the cut generated to (BendersCutSet) Solve BMP with early termination criteria to obtain ZL and Z U with its associated (¯y, v¯ c , v¯ d , η) ¯ BMP BMP itr L Set LB = Z BMP if LBitr > LB then LB = LBitr end if gap = (UB-LB)/UB end while return UB and the corresponding solution (z, y, vc , vd )

We first initialize the algorithm with best upper bound (UB) and lower bound (LB) values and the algorithmic parameters ϵ and Maxitr with specific values as shown in line 1. Next, we solve an initial BMP without any Benders cuts (BendersCutSet = ∅). Then, the main iteration (while) loop is started in which we first solve the DBSP, update the UB if necessary, generate Benders cuts using its solution, and solve the BMP for new values of the lower bound and the integer variables. Upon termination, the upper bound and the corresponding solution are returned. 4. Computational study To test the efficiency in terms of solution time and quality of our solution method, we solve randomly generated problems as well as data having demand distribution characteristics that reflect real applications. All of the computational studies were performed on a machine with Pentium Core 2 quad 3.0 GHz CPU with 8.0 GB RAM. The algorithms were implemented using C++ utilizing STL (Standard Template Library) and Concert Technology 2.0 with CPLEX 9.0 which we employ to solve the BMP and reduced DBSP. 4.1. Experimental setup Linehaul consolidated shipments mainly occur between large geographical regions; for example, shipments between eastern and western regions of the continental US or between northern and southern parts of the eastern (or western or central) US. To represent such underlying geography in our randomly generated problem instances, we create two identical squares of size 50 separated by a certain fixed distance A horizontally. The left square consists of only origin nodes and the right square consists of only destination nodes. We generate an equal number 500 of uniformly distributed point coordinates in each of the squares representing the potential physical origin and destination nodes. We also select M distinct nodes out of this set of points in each square to be the potential center locations and calculate all distances using the randomly generated point coordinates and the Euclidean norm. Then, we randomly select N distinct pairs from potential physical origin and destination nodes to determine N commodities. We note that this way of network generation approach is similar to

Fig. 1. Progression of UB and LB values with (solid) and without (dotted) strong cuts.

the approach used in other distribution network design studies including [3,5,19]. We consider a truck capacity U as 8 units, the base capacity S as 2 units, an A value of 0 (i.e., a single square representing the large region) or 50, and the values of α and β as 1 and 5, respectively. As for the size of the instances, we consider N values of 3000, 4000, 5000, 6000, 8000, and 10000; and M values of 20, 30, and 40. A triplet of (A, N , M ) values in a data set represents a problem class for which we solve 10 instances generated randomly. We generate the demand data by employing a demand structure that may more closely reflect real applications. Specifically, it can safely be assumed that the volume between some relatively small number of node pairs (e.g. warehouse-to-warehouse shipments) may be much larger than others. Based on this observation, we can categorize the demand distribution for the commodities. For this purpose, we generate the demands (wi ) in three categories representing 70%, 20%, and 10% of the commodities with demands randomly generated using U[0.1, 0.3], U(0.3, 0.6], and U(0.6, 0.9], respectively. 4.2. Computational results To evaluate the impact of employing strong cuts, we solve each instance with our Benders algorithm, given in Display 1, with and without strong cuts. The progression of UB and LB values over the iterations around the optimum is plotted in Fig. 1. Clearly, the use of strong cuts is very effective in speeding up convergence and they affect the quality of both the upper and lower bounds. Although the Fig. 1 is for one instance, in our numerical studies we observe similar high performance in convergence with strong cuts in other instances. In all of the following computational experiments, we employ strong cuts generated as described in Section 3.3. Our computational results are summarized in Table 1, in which we show average, minimum, and maximum solution times over 10 instances in each row corresponding to a class (A, N , M ). With only one exception, all of the instances are solved within 2% optimality in less than 400 iterations. An instance in class (A, N , M ) = (0, 10,000, 40) is solved to 2.07% optimality in 400 iterations; the maximum runtime of 4435 s belongs to this instance. Another outlier instance occurs in class (A, N , M ) = (50, 10,000, 40) with a runtime of 4834, as shown in Table 1, however, for this instance, the optimality gap is less than 2.0%.

H. Üster, H. Agrahari / Operations Research Letters 39 (2011) 138–143

143

Table 1 Summary of runtimes (in seconds). N

M = 20

A=0

Min

Ave

Max

M = 30 Min

Ave

Max

M = 40 Min

Ave

Max

3000 4000 5000 6000 8000 10,000

25 42 77 57 16 168

98 81 126 121 204 456

325 127 211 246 315 924

40 70 75 117 414 175

167 252 311 444 756 741

322 857 718 1215 1075 1137

37 173 166 196 613 823

326 567 609 1014 919 1939

814 1134 1139 2391 1461 4435

14 26 32 24 119 55

49 82 148 129 218 329

116 364 305 320 283 840

27 48 75 131 75 164

204 109 280 475 504 717

763 218 564 964 941 1537

140 255 190 332 149 483

374 449 625 513 1162 1772

739 655 1328 712 2150 4834

A = 50 3000 4000 5000 6000 8000 10,000

In Table 1, we observe that in general the runtimes are well within acceptable levels for the optimality gap of 2.0%. Even for quite large instances with 10,000 commodities (N value) and high number of potential consolidation and deconsolidation centers (M value), average runtimes are mostly well below 30 min and even the maximum runtimes are low. Clearly, increasing the problem size has an impact on average solution times, however, this does not appear to be too restrictive even when large instances are considered. For comparison purposes, using CPLEX to solve MnP with the addition of (23), we were only able to solve the instances (A, N , M ) classes including (0, 3000, 20), (0, 4000, 20), (50, 3000, 20), and (50, 4000, 20) with average solution times 168, 154, 139, and 159 s, and maximum solution times of 293, 267, 266, and 475 s, respectively. For the rest of the classes, the high memory requirement, largely due to the increased size of the model, is the main problem. In the Benders Decomposition approach, this issue is very efficiently alleviated since the problem is decomposed into smaller size main and subproblems, and more importantly, we solve for a large number of variables in DBSP using the second phase algorithm without relying on CPLEX. To emphasize the importance of this latter point, we also note that when we attempt to solve DBSP without reduction, runtimes become excessive and, for larger instances, memory again becomes an issue. 5. Concluding remarks This study can be extended to incorporate transportation cost functions that reflect economies-of-distance along with economies-of-scale on the linehaul links. Furthermore, our model can be extended to incorporate vehicle routing considerations in the origin-to-consolidation and deconsolidation-to-destination portions of the operations. Acknowledgements The authors thank the associate editor and referees for their helpful comments on an earlier version of this paper. References [1] S. Abdinnour-Helm, A hybrid heuristic for the uncapacitated hub location problem, European Journal of Operational Research 106 (2) (1998) 489–499.

[2] S. Abdinnour-Helm, M.A. Venkataramanan, Solution approaches to hub location problems, Annals of Operations Research 78 (1998) 31–50. [3] A. Amiri, H. Pirkul, New formulation and relaxation to solve a concave-cost network flow problem, Journal of Operational Research Society 48 (1997) 278–287. [4] A. Atamtürk, On capacitated network design cut-set polyhedra, Mathematical Programming, Series B 92 (2002) 425–437. [5] A. Balakrishnan, S.C. Graves, A composite algorithm for a concave-cost network flow problem, Networks 19 (1989) 175–202. [6] J.F. Benders, Partitioning procedures for solving mixed-variables programming problems, Numerische Mathematik 4 (1962) 238–252. [7] D.J. Bowersox, D.J. Closs, M.B. Cooper, Supply Chain Logistics Management, The McGraw-Hill Companies, Inc., New York, NY, 2002. [8] J.F. Campbell, Integer programming formulations of discrete hub location problems, European Journal of Operational Research 72 (1994) 387–405. [9] J.F. Campbell, A.T. Ernst, M. Krishnamoorthy, Hub location problems, in: Z. Drezner, H.W. Hamacher (Eds.), Facility Location: Applications and Theory, Springer-Verlag, New York, NY, 2002, pp. 373–407. [10] J.F. Chen, A hybrid heuristic for the uncapacitated single allocation hub location problem, Omega 35 (2) (2007) 211–220. [11] C.B. Cunha, M.R. Silva, A genetic algorithm for the problem of configuring a hub-and-spoke network for a LTL trucking company in Brazil, European Journal of Operational Research 127 (3) (2007) 747–758. [12] J. Ebery, M. Krishnamoorthy, A.T. Ernst, N. Boland, The capacitated multiple allocation hub location problem: Formulation and algorithms, European Journal of Operational Research 120 (2000) 614–631. [13] K. Holmberg, J. Hellstrand, Solving the uncapacitated network design problem by a Lagrangean heuristic and branch-and-bound, Operations Research 46 (2) (1998) 247–259. [14] T.L. Magnanti, P. Mirchandani, R. Vachani, The convex hull of two core capacitated network design problems, Mathematical Programming 60 (1993) 233–250. [15] T.L. Magnanti, P. Mirchandani, R. Vachani, Modeling and solving the twofacility capacitated network loading problem, Operations Research 43 (1) (1995) 142–157. [16] T.L. Magnanti, P. Mireault, R. Wong, Tailoring benders decomposition for uncapacitated network design problem, Mathematical Programming 26 (1986) 112–154. [17] T.L. Magnanti, R.T. Wong, Accelerating Benders decomposition: Algorithmic enhancement and model selection criteria, Operations Research 29 (3) (1981) 464–484. [18] A. Marin, Formulating and solving splittable capacitated multiple allocation hub location problems, Computers and Operations Research 32 (12) (2005) 3093–3109. [19] A. Muriel, F.N. Munshi, Capacitated multicommodity network flow problems with piecewise linear concave costs, IIE Transactions 36 (2003) 683–696. [20] M. Randall, Solution approaches for the capacitated single allocation hub location problem using ant colony optimisation, Journal Computational Optimization and Applications 39 (2) (2008) 239–261. [21] T.J. Van Roy, A cross decomposition algorithm for capacitated facility location, Operations Research 34 (1986) 145–163. [22] P. Wentges, Accelerating Benders decomposition for the capacitated facility location problem, Mathematical Methods of Operations Research 44 (2) (1996) 267–290.