A highly efficient exact algorithm for the uncapacitated multiple allocation p-hub center problem

Article history: Received October 10, 2019 Received in revised format: October 28, 2019 Accepted December 25, 2019 Available online December 26, 2019 Globalization and increasing competition in global markets have forced businesses to provide a high level of service to their customers. Time-sensitive transportation systems which are used in transportation of perishable goods, express mail delivery, and emergency services are playing a very important role in this regard. This paper addresses the problem of uncapacitated multiple allocation p-hub center problem (UMApHCP) which is fundamental in proper functioning of time-sensitive transportation systems. A mixed-integer programming formulation is proposed for the problem and a highly efficient Benders decomposition algorithm is developed for solving it. The proposed algorithm is capable of solving large-scale instances of the problem to optimality in order of seconds. . by the authors; licensee Growing Science, Canada 20 20 ©


Introduction
Hubs are key facilities that act as intermediate centers for switching, sorting, and consolidating the commodities in many-to-many transportation systems and telecommunication networks. Generally, direct shipping of these commodities is not an economical decision as establishing direct linkages between each origin-destination (O/D) pair is extremely costly. Therefore, smaller number of connections are built in hub networks and each connection carry large traffic volume which in turn makes it possible to exploit economies of scale in transportation costs, especially on the inter-hub connections.
The hub location problem (HLP) deals with locating the hub facilities in the network and determining the way the non-hub nodes are allocated to the established hubs and the pattern based on which the O/D flow are routed via the hub facilities in such a way that a specific objective function is optimized. Regarding way the non-hub nodes are allocated to the hubs, we have two main types of hub networks. If each non-hub node send and receive the associated traffic through exactly one hub, then we deal with a single allocation hub network. On the other hand, if there is no restriction on the number of hubs with which each non-hub node can communicate, then the resulting network is called a multiple allocation hub network. A typical multiple allocation hub network is depicted in Fig. 1. In this figure, circles show the non-hub nodes, whereas the triangles represent the hub facilities. Most of the studies on the HLP deal with cost minimization objectives including fixed cost for establishing hubs and variable transportation costs. Although reducing cost is the primary concern in many practical problems in transportation, postal services, and telecommunication industries, it sometimes leads to unsatisfactory solutions in terms of service equitability as the distance between the O/D pairs might be excessively large in the worst case. This is an undesirable outcome particularly when we deal with the delivery of perishable or time-sensitive items or when the aim is to provide a socially equitable service for different O/D traffic which is essential for a sustainable transportation system. To remedy this drawback, one option is to use hub center problem in order to minimize the maximum distance (or cost) between the O/D pairs. In many situations in practice, the transported commodities are too sensitive to transportation time and need fast delivery. A possible real life application of this problem could be the location of hubs for transportation of perishable items (fruits, vegetables, seafood, etc.) where the hubs need to be located in a way that the clients could be served as soon as possible to avoid the damage of the products. Another application could be location of hubs in postal services where the express services are required by some customers for different reasons. In this work we address the uncapacitated multiple allocation p-hub center (UMApHCP). The goal is to locate a fixed number of hub facilities in such a way that the largest distances between the O/D pairs is minimized. The considered problem is modeled as mixed-integer linear program (MILP) and a highly efficient exact solution approach based on Benders decomposition is proposed to solve it. Extensive computational experiments are conducted to study the efficiency of the proposed solution algorithm and the effect of different input parameters on the optimal solutions of the problem.
The remainder of this paper is organized as follows. The literature review is presented in Section 2. A mathematical formulation for the problem is presented in Section 3. Section 4 describes the proposed Benders decomposition algorithm. Numerical results are presented in Section 5. Finally, Section 6 concludes the paper and provides some directions for future research.

Background
The HLP has become an important and well-studied area of research since 1980s within the field of facility location theory. The number of published studies on the HLP has had an increasing trend in recent years and different variants of the HLP have been studied so far. Recent surveys can be found in Alumur and Kara (2008), Campbell and O'Kelly (2012), and Farahani et al. (2013) for the interested readers. One of the important types of the HLP is the p-hub center problem which aims to serve the O/D demand pairs by establishing a fixed number of p hubs in such a way that the maximum transportation cost (time, distance, etc.) is minimized. This problem has many applications in transportation of time-sensitive commodities and hence it has been studied by numerous researches. Mixed integer programming (MIP) formulations for the single and multiple allocation p-hub center problem was proposed for the first time by Campbell (1994). Kara and Tansel (2000) developed three different integer linear programming (ILP) formulations of the single-allocation p-hub center problem, and found computationally that one of the formulations to be consistently superior. Kratica and Stanimirović (2006) devised a genetic algorithm (GA) for the UMApHCP. Ernst et al. (2009) proposed a new ILP formulation for the uncapacitated single allocation p-hub center problem (USApHCP) using a new decision variable for the maximum collection and distribution cost/time between a hub and the associated O/D nodes. Sim et al. (2009) considered a p-hub center problem with normally distributed stochastic travel times and employed chance constraints to control the probability of the total travel time along a path to not exceed a given time bound. Meyer et al. (2009) presented a 2-phase algorithm for the USApHCP where a set of potential optimal hub combinations was computed using a shortest path based branch and bound. The allocation phase was done using a reduced sized formulation giving the optimal solution. Furthermore, they developed a heuristic based on an ant colony optimization (ACO) approach to provide good upper bound for their algorithm. Bashiri et al. (2013) presented a GA based heuristic to solve the fuzzy capacitated p-hub center problem. Yang et al. (2013) developed a hybrid particle swarm optimization (PSO) algorithm for fuzzy p-hub center problem where the travel times were modeled using normal fuzzy vectors. Brimberg et al. (2017a) proposed a basic variable neighborhood search (VNS) heuristic for the UMApHCP. In another work, Brimberg et al. (2017b) proposed a general variable neighborhood search heuristic for the USApHCP.
Benders decomposition (Benders, 1962) is a method based on row generation which is suitable for solving MIP problems. In this procedure, the problem is reformulated by partitioning it into two simpler problems and a cutting plane approach is used to solve the reformulated problem. Benders decomposition algorithm has been successfully applied to the HLP in numerous research works. de  (2011) developed an accelerated Benders decomposition algorithm for an uncapacitated multiple allocation HLP tailored for urban transport and liner shipping network design problem. The authors extended the same Benders decomposition algorithm for the multi-period uncapacitated multiple assignment HLP under budget constraint in (Gelareh et al., 2015). Contreras et al. (2012) proposed a Benders algorithm for an extension of the classical capacitated multiple allocation HLP in which the amount of capacity installed at the hubs is a decision variable. Another Benders decomposition algorithm is devised for the manyto-many hub location routing problem by de Camargo et al.

Mathematical formulations
Assume that G = (N, A) is a graph with N as the set of nodes and A as the set of arcs such that A={(i,j): i,j ∈ N, i < j }. Let H ⊆ N represent the subset of candidate nodes for locating hubs. The number of hubs to be located is fixed and is denoted by p. For all i, j ∈ N, let dij represent the distance between nodes i and j. The total routing distance for the O/D flow associated with pair (i,j) ∈ A via hubs k and m in that order is calculated as: where α is the discount factor to reflect economies of scale on the inter-hub connections (0 < α < 1), which can be interpreted as a speed-up (or cost decrease) factor incurred by usage of a faster (or more efficient) means of transport on such connections.
To develop an MIP formulation for uncapacitated multiple allocation p-hub center problem (UMApHCP), we use the binary variable zk ∈ {0,1} to be 1 if node k is selected as a hub and 0, otherwise. Also, the non-negative variable xijkm denotes the fraction of traffic associated with the O/D pair (i,j) ∈ A, that is routed via hubs k and m in that order. Further, let β represent the maximum distance between the O/D pairs in the network. The problems consist of selecting p nodes as hubs and determining how the O/D flows will be assigned to these hubs in such a way that the maximum distances between O/D pairs is minimized. The MIP model for the UMApHCP can be written as: (1) min  s.t: The objective function (1) minimizes the maximum distances between O/D pairs. Constraint (2) determines the number of hubs to be located in the network. Constraints (3) assure that the whole flow associated with each O/D pair is routed via some hub pair. Constraints (4) state that the flows can only be routed via nodes that have been designated as hubs. Constraints (5) calculate the maximum distance between the O/D pairs. (6), (7), and (8) are the standard domain constraints for the decision variables.

Benders reformulation
Benders decomposition is a row generation based exact solution method that can be applied to solve large-scale mixed integer programming problems (Benders, 1962). In this technique, the problem is reformulated using a smaller number of variables and a large number of constraints. The problem is then solved using a cutting plane approach in which the relaxed problem, called as the master problem, is solved at each iteration by adding cutting planes found by the subproblem. Benders reformulation exploits the fact that computational burden of an MIP problem substantially increases with the problem size. Therefore, a divide-and-conquer scheme is employed to decompose a single large problem into smaller problems which can be solved more efficiently in terms of the computational time and resources. Motivated by this fact, in this work we apply Benders decomposition to the UMApHCP. In classical implementation of Benders algorithm, we need to solve the master problem at each iteration. However, we use a modern implementation within branch-and-cut framework where the master problem is solved in a single attempt and the cuts are added on the fly whenever required by utilizing recent developments in off-the-shelf solvers. We separate Benders cuts whenever a candidate integer solution is found in the branch-and-cut tree of the master problem. By doing so, the computational effort needed to solve an integer problem at each iteration is considerably reduced.

The subproblem and the master problem
By fixing the binary location variables vector as ẑ = z , the subproblem (SP) for the UMApHCP can be written as: Let λij, θijk, and μij be the dual variable associated with constraints (10), (11), and (12) in the SP, respectively. Hence, the dual subproblem (DSP) for UMApHCP can be written as follows: We can now write the master problem (MP) for UMApHCP as follows: in which (λ t ,θ t ,μ t ) is the t th extreme point of the feasible solution space of the DSP. Observe that constraint (24) assures the installation of p hubs in the network which in turn guarantees the feasibility of the SP at any iteration. Therefore, there is no need for adding feasibility cuts to the MP in our problems.

Solving the dual subproblem
In this section we devise algorithms for solving the DSP by inspection (without using standard solver) which resulting in substantial savings is the solution time of the problems. To this end, we first determine the optimal values of μij variables. Due to the maximization sense in the objective function of the DSP (15)-(21) and according to constraints (18), we can conclude that the value of variable μij associated with the O/D pair with longest distance will take 1 and the remaining variables will have 0 value in an optimal solution. We present a procedure for determining the optimal values of μij variables in Algorithm 1. In the above algorithm H1 is the set of opened hubs at the current iteration. Having obtained the optimal values for the μij variables, the DSP for the UMApHCP can now be reduced to the following linear programming problem: ' ' i j    This problem can be solved by inspection and the corresponding optimal dual variables (λi'j' and θi'j'k) can be determined using the complementary slackness conditions as explained by Contreras et al. (2011a). Furthermore, since for only one O/D pair (i.e., (i',j')) the values of the dual variables are nonzero, our master problem reduces to the following much smaller MIP which can be solved more efficiently.

Test problems
We conduct an extensive set of computational experiments in order to demonstrate the efficiency of the proposed Benders decomposition algorithms. To this end, we use three well-known data sets from the hub location literature, namely the CAB, the TR, and the AP data sets. The CAB data set constitute the airline passenger interactions data between 25 US cities in 1970 and is introduced by O'Kelly (1987). Our second data set is the TR data set (Tan and Kara, 2007) which is based on the cargo flows between 81 cities of Turkey. The third data set is the Australia Post (AP) data set which is introduced by Ernst and Krishnamoorthy (1996)  We adopt a modern implementation of the Benders algorithm within a branch-and-cut framework, where the master problem is solved on a single tree and the optimality cuts are added one at a time using the lazy constraint callback function available in CPLEX. The time limit of two hours is used in our experiments. Table 1 shows the results obtained by solving the UMApHCP with the CAB data set. The columns entitled p and α denote the number of opened hubs and the value of the discount factor, respectively. The next two columns show the optimal objective function value and the corresponding set of opened hubs. The next columns show the solution times (in seconds) for the MIP model and the Benders decomposition algorithm.

Numerical results
The results reported in Table 1 show that the proposed Benders algorithm is able to obtain the optimal solution for all the instances of the CAB data set in very short computational time. The time needed to obtain the optimal solution by using the Benders algorithm is in order of hundredth of a second, whereas the average time taken by CPLEX to solve the MIP model is around 200 seconds. Note that for a given number of the opened hubs p, the optimal set of hub facilities change substantially as the value of discount factor varies. Observe also that for a fixed value of p, the optimal objective value increases as the value of the discount factor (α) increases. In contrast, for a fixed value of the discount factor, the value of the total transportation cost decreases as the number of opened hubs (p) increases.  Table 2 shows the results obtained by solving the UMApHCP with the TR data set. As can be seen from the results, the proposed Benders algorithm has solved all the instances in quite short computational times. Note that the maximum solution times for different instances of the real-sized TR data set are around half a minute. Nevertheless, the MIP model could not be solved by CPLEX because of excessive memory requirements. Moreover, the changes in the optimal value of the objective function as a result of varying values for the input parameters p and α are similar to the corresponding changes with the CAB data set.  Table 3 presents the results of solving the problem for the AP data set with small instancs (|N|=25 and 50). It can be seen from this table that the Benders algorithm all the instances in less than one second of CPU time. The MIP model could solve th instances with |N|=25 using CPLEX in average time of around 220 seconds. Note that for the problem with |N|=50 only one of the instances could be solved by CPLEX within the allowed time of two hours. For the remaining instances which was not solved within the two hours limit, the gap percentage between the corresponding upper and lower bounds are reported inside the parentheses. The results obtained by solving the medium size AP instances (|N|=75 and 100) are presented in Table  4. While none of the instances could be solved as MIP model by CPLEX, the Benders algorithm could solve all the instances of medium sizes in less than three seconds. These results clearly demonstrate the high efficiency of the proposed solution algorithm.
Finally, Table 5 shows the results obtained by solving the large-scale AP instances (|N|=150 and 200). All the instances of size 150 are solved in less than 8 seconds. For the instances with 200 nodes, the maximum solution time is around 70 seconds. To the best of our knowledge, this is the first time that the instances of up to 200 nodes for the UMApHCP are solved to optimality in the literature. The solution times are extremely small for solving such large problem instance. This shows that even larger instances can be solved using the proposed exact solution algorithm in short computational times. Since the UMApHCP belongs to the class of NP-hard problems, proposing an exact solution algorithm which is able to solve large-scale instances of the problem in order of seconds is of utmost theoretical and practical importance.

Conclusions
We considered the uncapacitated multiple allocation p-hub center problem (UMApHCP) and formulated the problem as linear mixed integer programming model. In order to solve large-scale instances of the problem, a highly efficient Benders decomposition algorithm was proposed. An extensive set of computational experiments were done in order to analyze the efficiency of the proposed solution algorithm and to study the effect of different input parameters on the final solutions of the model. Obtained results demonstrate the capability of the proposed exact algorithm to solve large-scale instances of the problem in short computational times for the first time in the literature. Since the UMApHCP belongs to the class of NP-hard problems, proposing an exact solution algorithm which is able to solve large-scale instances of the problem in order of seconds is of utmost theoretical and practical importance. An interesting line for further extension of this research is to consider some sources of uncertainty (such as transportation costs, distances, etc.) in the proposed model. Furthermore, one may relax the classical HLP assumptions such as complete inter-hub network or fixed flow-independent discount factor for representing the economies of scale in order to deal with more realistic real-world situations. Finally, the proposed solution algorithm can be applied to other variants of the hub location problem with service level considerations such as hub covering problems.