Solving Hub Location Problems With Profits Using Variable Neighborhood Search

Hub Location Problems

HLPs can be classified by allocation strategies. Three types of allocation strategies: single allocation, multiple allocation, and r-allocation, have been studied in the literature. The single allocation HLP in which each spoke is only allocated to one hub was first studied. The uncapacitated single allocation model was first proposed by O’Kelly ( 13 ). Campbell ( 14 ) and Ernst and Krishnamoorthy ( 15 ) studied different mathematical models. The capacitated single allocation problem has also been studied by Ernst and Krishnamoorthy ( 16 ) and Contreras et al. ( 17 , 18 ). Multiple allocation problems assign each spoke to hubs without limiting the number of hubs. The uncapacitated multiple allocation problem was formulated with integer programming (IP) by Campbell ( 14 ) and Ernst and Krishnamoorthy ( 19 , 20 ). Examples of studies for capacitated multiple allocation problems are Ebery et al. ( 21 ), Rodriguez-Martin and Salazar-Gonzalez ( 22 ) and Contreras et al. ( 23 ). In r-allocation problems, each spoke can be assigned to at most r hubs. Yaman ( 6 ) first introduced a new r-allocation problem to generalize two variants of the multiple and single allocation strategies.

Early studies on HLPs assumed that the hub network is complete with a direct link between each hub pair. However, the complete hub network assumption may be unnecessary in some applications because it increases the total costs of designing hub networks ( 24 ). Therefore, many researchers studied HLPs with incomplete hub networks in which some hub pairs are not linked. There are several examples of designing incomplete hub networks ( 25 – 29 ).

Direct connections (DCs) between spokes are not allowed in traditional HLPs ( 5 ). However, some OD pairs are very close, and flows between them travel through long paths of allocation arcs and hub links. This leads to more costs against the original intention of applying the hub-and-spoke network structure. Gelareh and Nickel ( 28 ), de Camargo et al. ( 29 ) and O’Kelly et al. ( 30 ) took into account the direct arcs when modeling HLPs.

Classical HLPs aim to minimize the total costs of transferring demands. Some studies also consider the additional costs of setting up hubs and hub links. With minimizing total costs in the system, all requests are satisfied between OD pairs. The cost minimization HLPs have been widely applied to many real-life applications. Typical applications can be found in postal delivery services and supply chain management logistics. Lin et al. ( 31 ) proposed a general capacitated p-hub median model and applied it to design a Chinese air cargo network. Çetiner et al. ( 32 ) considered a combined hub location and routing problem in which the first-stage problem is a HLP to minimize total transportation costs. They tested their model on the Turkish postal delivery system. Ishfaq and Sox ( 33 ) extended the p-hub median problem to the domain of intermodal logistics.

Hub Location Problems With Profits

In recent studies, some authors have modeled HLPs to maximize profits, which allow for remaining demands not served. A mixed-integer programming (MIP) formulation and a Lagrangian relaxation (LR) algorithm were presented by Alibeyg et al. ( 34 ) to solve uncapacitated HLPs with profits (HLPPs). Alibeyg et al. ( 35 ) formulated a class of hub network design problems with profit-oriented objectives as IP formulation. An exact algorithm based on the LR and branch-and-bound (BB) was proposed by Alibeyg et al. ( 36 ) to solve the models in Alibeyg et al. ( 35 ). Lin and Lee ( 10 ) compared the differences between cost-oriented objectives and profit-oriented objectives, and formulated the hub network design LTL problem with maximizing total profits under the price elasticity of demands. Taherkhani and Alumur ( 8 ) presented different models for HLPPs with fixed costs under different allocation strategies. The assumptions that the hub network is complete and DCs are not allowed were relaxed in their models. Wu et al. ( 37 ) studied a HLPP with service capacity and flow capacity constraints by formulating a MIP model. They applied a LR-based heuristic to solve the problem which showed good effectiveness. However, they considered the complete hub network in their problem. Taherkhani et al. ( 38 ) modeled the capacitated HLPP with multiple demand classes, which considered the complete hub network. A Benders decomposition method was developed to solve the determined version of the problem, and a Monte Carlo simulation-based algorithm was developed to solve the stochastic one.

There are also some variants of HLPs that consider maximizing profits. Some authors maximized the profits in a competitive scenario. Lüer-Villagra and Marianov ( 3 ) formulated and solved a hub location-pricing problem to maximize profits for the new company in a competitive situation by choosing the optimal hub network design and pricing, and a logit discrete choice model was introduced to capture the customers’ behavior. Zhao et al. ( 39 ), Abbasi and Niknam ( 40 ), and Esmaeili and Sedehzade ( 41 ) also studied the hub network design and pricing problem with competitive situations. Sharma et al. ( 9 ) proposed a two-phase mixed-integer linear programming model for the “coopetitive” profit maximization HLP in the airline industry. Hou et al. ( 42 , 43 ) integrated the revenue management and hub location into a revenue management problem to maximize the profits. Taherkhani et al. ( 44 ) studied profit maximization capacitated multiple allocation HLPs with the complete hub network under two different types of uncertainty: stochastic demands and uncertain revenue. To model the problems in different scenarios, a two-stage stochastic programming model and robust stochastic models were proposed and solved by Benders decomposition algorithms with a sample average approximation scheme. The results show that the proposed algorithms were able to solve optimally instances with up to 75 nodes.

Applications of HLPPs can be found in freight delivery. Lin and Lee ( 10 ) designed a hub network for LTL freight delivery to maximize total profits under the price elasticity of demand. They tested it with the second largest time-definite LTL freight delivery common carrier in Taiwan to show different behaviors between cost minimization and profit maximization.

Algorithms to Solve Hub Location Problems

Most HLPs are NP-hard, which cannot be solved on a large scale in polynomial time ( 11 , 45 ). To solve HLPs to optimality, many solution techniques have been proposed, including BB algorithms ( 8 , 46 , 47 ), LR-based algorithms ( 35 ), and Benders decomposition algorithms ( 29 , 38 , 48 , 49 ). The first work to use a Benders decomposition to solve the uncapacitated multiple allocation HLP was de Camargo et al. ( 48 ). Contreras et al. ( 49 ) developed an enhanced Benders decomposition algorithm for the uncapacitated multiple allocation HLP by a multicut reformulation, Pareto-optimal cuts, reduction tests, and a heuristic procedure. Taherkhani et al. ( 38 ) applied the Benders reformulation to solve a capacitated HLPP with the complete hub network under multiple demand classes.

Since it is difficult to formulate exact algorithms to solve large-scale HLPs, several heuristics have been proposed for this purpose. Tabu search heuristics were developed by Klincewicz ( 50 , 51 ) and Skorin-Kapov and Skorin-Kapov ( 52 ). A simple simulated annealing algorithm was proposed by ( 15 ). LR-based heuristics were studied by Pirkul and Schilling ( 53 ) and Elhedhli and Wu ( 54 ). Kratica et al. ( 55 ) used genetic algorithms to solve the uncapacitated p-hub median problem with single assignments. VNS was proposed by Mladenović and Hansen ( 56 ), and is commonly used to solve HLPs. Ilić et al. ( 57 ) developed a general VNS for solving the uncapacitated single allocation p-hub median problem. They proposed a new nested variable neighborhood descent to reduce the number of local minima. Recently, Mikić et al. ( 58 ) introduced VNS to solve the capacitated modular HLP. Dai et al. ( 59 ) used a VNS method to refine the solution in their implementation of a heuristic-based method called HUBBI. Wandelt et al. ( 60 ) proposed a network contraction approach to speed up the genetic, VNS, and Benders decomposition algorithms in six different HLPs. Applying the contraction technique, they solved the problem with networks of sizes up to 5,000 nodes. Fontes and Goncalves ( 61 ) presented a new hub-and-spoke network problem with multiple allocations of sub-hub nodes for liner shipping operations and applied a variable neighborhood decomposition search algorithm to solve the proposed model.

To the best of the authors’ knowledge, except for Alibeyg et al. ( 34 ) and Wu et al. ( 37 ), HLPPs are mainly solved with exact methods. Alibeyg et al. ( 34 ) presented an uncapacitated HLPP with the complete hub network and proposed a LR-based algorithm, which is able to obtain tight lower and upper bounds, and a primal heuristic was used to find the feasible solution. Wu et al. ( 37 ) also discussed the HLPP under the complete hub network and applied a LR-based heuristic to solve the problem. However, both of those papers considered the complete hub network in their problems. They applied a LR-based heuristic to solve the problem which showed good effectiveness. However, when the problem is the incomplete hub network and relaxes the two-hop constraint, it would not scale up well because LR is often slow and is soon stuck. Therefore, heuristics are needed to solve HLPPs with the incomplete hub network.

Multiple Allocation Problem

In the multiple allocation HLPP, the hub network is determined with the locations of hubs and the selection of hub links, and then spokes are allocated to hubs with no limit to the number of hubs. The multiple allocation problem with profits is modeled as follows:

max ∑ i ∈ N ∑ j ∈ N ∑ k ∈ N ∑ l ∈ N ( r ij − c ik − c lj ) w ij y ijkl − ( ∑ i ∈ N ∑ k ∈ N ∑ l ∈ N α c kl f ikl + ∑ k ∈ N f k h k + ∑ k ∈ N ∑ l ∈ N g kl z kl )

(1)

s.t.

∑ k ∈ N ∑ l ∈ N y ijkl ≤ 1 i , j ∈ N

(2)

∑ l ∈ N y ijkl + ∑ l ∈ N , l ≠ k y ijlk ≤ h k i , j , k ∈ N

(3)

∑ l ∈ N , l ≠ k f ikl − ∑ l ∈ N , l ≠ k f ilk = ∑ j ∈ N ∑ l ∈ N w ij y ijkl − ∑ j ∈ N ∑ l ∈ N w ij y ijlk i , k ∈ N

(4)

f ikl ≤ O i z kl i , k , l ∈ N , k ≠ l

(5)

z kl ≤ h k k , l ∈ N , k ≠ l

(6)

z kl ≤ h l k , l ∈ N , k ≠ l

(7)

f ikl ≥ 0 i , k , l ∈ N

(8)

h k ∈ { 0 , 1 } k ∈ N

(9)

y ijkl ∈ { 0 , 1 } i , j , k , l ∈ N

(10)

z kl ∈ { 0 , 1 } k , l ∈ N , k ≠ l

(11)

Equation 1 is the objective function denoting total net profits. The first part is the revenue subtracting transportation costs between spokes and hubs for OD pairs. The second part covers all costs related to hubs, including transportation costs through hub links, hub installation costs, and hub link installation costs. Equation 2 ensures that the path through the hubs must be unique for a given OD pair. Equation 3 is a constraint that the path of the demand between OD pairs must be served through the open hubs. Equations 4 and 5 are flow balance constraints. Equations 6 and 7 denote that hub links are set up between hubs. Equations 8 to 11 are domain constraints of decision variables. There are ( | N | 4 + | N | 2 ) binary variables and ( | N | 3 ) real variables. The number of constraints of the model is ( 2 | N | 3 + 3 | N | 2 − 2 | N | ) . Therefore, the complexity of the multiple allocation model without direct links is O ( | N | 4 ) variables and O ( | N | 3 ) constraints.

Single Allocation Problem

For the uncapacitated single allocation HLPP, each spoke can only be allocated to one hub. Thus, it needs a binary variable x ik to denote the specific hub to which the spoke is allocated. Note that when i = k the variable x kk denotes the installation of a hub at the node k . Therefore, the variable h k indicating the same implication is removed in the single allocation model. The model for the single allocation problem is formulated as follows:

max ∑ i ∈ N ∑ j ∈ N ∑ k ∈ N ∑ l ∈ N ( r ij − c ik − c lj ) w ij y ijkl − ( ∑ i ∈ N ∑ k ∈ N ∑ l ∈ N α c kl f ikl + ∑ k ∈ N f k x kk + ∑ k ∈ N ∑ l ∈ N g kl z kl )

(12)

s.t. (2), (4)-(5), (8), (10)-(11)

∑ k ∈ N x ik ≤ 1 i ∈ N

(13)

x ik ≤ x kk i , k ∈ N

(14)

y ijkl ≤ x ik i , j , k , l ∈ N

(15)

y ijkl ≤ x jl i , j , k , l ∈ N

(16)

z kl ≤ x kk k , l ∈ N , k ≠ l

(17)

z kl ≤ x ll k , l ∈ N , k ≠ l

(18)

x ik ∈ { 0 , 1 } i , k ∈ N

(19)

The objective function (Equation 12) is similar to that of the multiple allocation model except for the cost of opening hubs. Equations 2, 4, 5, 8, 10, and 11 in this model have the same explanation as those in the multiple allocation model. Equation 13 denotes that each node is allocated to at most one hub. Equation 14 ensures that nodes can only be allocated to opening hubs. Equations 15 and 16 guarantee that paths of OD pairs pass the allocated hubs. Equations 17 and 18 denote that hub links are set up between open hubs. Equation 19 requires that decision variables x ik are binary. There are ( | N | 4 + 2 | N | 2 − | N | ) binary variables, ( | N | 3 ) real variables, and ( 2 | N | 4 + | N | 3 + 4 | N | 2 − | N | ) constraints in the model of USAI. The complexity of variables and the complexity constraints in USAI are therefore both O ( | N | 4 ) .

R-Allocation Problem

For the r-allocation problem, each spoke can be allocated to r hubs at most, which is the general case for the single allocation problem with r = 1 and the multiple allocation problem with r = | N | . It can be modeled with the same formulation as the single allocation model by replacing Equation 13 with the following constraints:

∑ k ∈ N x ik ≤ r i ∈ N

(20)

Hub Location Problems Allowing for Direct Connections

In the real-world network, direct transportation between OD pairs is sometimes allowed, thus it is necessary to discuss HLPs allowing for DCs. An additional binary variable s ij deciding the installation of a direct link from the node i to the node j is introduced.

The multiple allocation problem allowing for DCs is modeled in the following formulation:

max ∑ i ∈ N ∑ j ∈ N ∑ k ∈ N ∑ l ∈ N ( r ij − c ik − c lj ) w ij y ijkl + ∑ i ∈ N ∑ j ∈ N ( r ij − c ij ) w ij s ij − ( ∑ i ∈ N ∑ k ∈ N ∑ l ∈ N α c kl f ikl + ∑ k ∈ N f k h k + ∑ k ∈ N ∑ l ∈ N g kl z kl + ∑ i ∈ N ∑ j ∈ N q ij s ij )

(21)

s.t. (3)-(11)

s ij + h j ≤ 1 i , j ∈ N

(22)

s ij + h i ≤ 1 i , j ∈ N

(23)

s ij + ∑ k ∈ N ∑ l ∈ N y ijkl ≤ 1 i , j ∈ N

(24)

s ij ∈ { 0 , 1 } i , j ∈ N

(25)

Compared with the version of UMAI, profits from the demand of OD pairs being transported directly and costs operating direct links are considered in the objective function (Equation 21), and four sets of constraints related to DCs are added. Equations 22 and 23 ensure that DCs are between non-hub nodes. Equation 24 ensures that the demand of each OD pair is served by only one path. The model has ( | N | 4 + 2 | N | 2 ) binary variables and ( | N | 3 ) real variables, and ( 2 | N | 3 + 5 | N | 2 − 2 | N | ) constraints. The multiple allocation model allowing for DCs therefore has the same scale as the version without DCs, in the complexity of the number of variables and the number of constraints.

The single allocation problem allowing for DCs is formulated as follows:

max ∑ i ∈ N ∑ j ∈ N ∑ k ∈ N ∑ l ∈ N ( r ij − c ik − c lj ) w ij y ijkl + ∑ i ∈ N ∑ j ∈ N ( r ij − c ij ) w ij s ij − ( ∑ i ∈ N ∑ k ∈ N ∑ l ∈ N α c kl f ikl + ∑ k ∈ N f k x kk + ∑ k ∈ N ∑ l ∈ N g kl z kl + ∑ i ∈ N ∑ j ∈ N q ij s ij )

(26)

s.t. (4), (5), (8), (10), (11), (13)–(19), (24), (25)

s ij + x jj ≤ 1 i , j ∈ N

(27)

s ij + x ii ≤ 1 i , j ∈ N

(28)

Equation 26 is the objective function to maximize total profits. Besides the same constraints as USAI, three sets of constraints (Equations 25, 27, and 28) about DCs are included in USAIDC. Equations 27 and 28 ensure that DCs are only operated between non-hub nodes. The model has ( | N | 4 + 3 | N | 2 − | N | ) binary variables, ( | N | 3 ) real variables, and ( 2 | N | 4 + | N | 3 + 7 | N | 2 − | N | ) constraints, which also has the same complexity as the version not allowing for DCs.

After replacing the constraint (Equation 13) with Equation 20 in USAIDC, the formulation of the r-allocation model allowing for DCs is obtained.

Solution Representation

The goal of HLPPs is to determine hub locations, inter-hub links, assignments of spokes, and travel paths of the served demand between OD pairs. The solution structure includes two parts: the hub network and the assignments of nodes. A solution is denoted as S = ( G , A ) , where G is a directed graph of the hub network and A is a series of assignments of nodes. G = ( H , E ) is the hub network, where H is the hub set and E is the inter-hub link set. In multiple allocation models, nodes can be allocated to any hub. Thus, these problems do not include assignments in the solution. In single allocation problems, A [ i ] represents the hub to which node i is assigned. In r-allocation problems, A [ i ] is a list of hubs with a size less than r .

The path from node i to node j with the first hub k and the last hub l is defined as h ij a = ( k , l ) in multiple allocation problems. The set of nodes that are assigned to the same hub k is defined as a cluster Cl u k in single allocation and r-allocation problems. Because the incomplete network is considered, a matrix C h = [ C kl h ] is introduced to denote unit demand cost to travel between hub pairs. The unit demand cost C kl h from hub k to hub l is obtained from the Dijkstra algorithm.

In problems without DCs, matrix R e = [ r ij e ] is used to record the unit net revenue of each OD pair to update the objective value. For UMAI, the unit net revenue of the OD pair ( i , j ) is r ij e = max { max { r ij − c ik − c jl − α C kl h , k , l ∈ H } , 0 } . For USAI, the unit net revenue of the OD pair ( i , j ) is r ij e = max { r ij − c iA [ i ] − c jA [ j ] − α C A [ i ] A [ j ] h , 0 } . For r-allocation models, the unit net revenue of the OD pair ( i , j ) is r ij e = max { max { r ij − c ik − c jl − α C kl h , k ∈ A [ i ] , l ∈ A [ j ] } , 0 } . The net revenues to serve the demand between OD pair ( i , j ) are r ij e w ij .

In problems allowing for DCs, the costs of transferring flows directly are the sum of transportation costs and setup costs of direct links, that is, C ij dl = c ij w ij + q ij . Thus, the total net revenues of pairs are obtained by subtracting the total costs from the revenues of the demand, expressed as r ij w ij − C ij dl . The total net revenue matrix R t = [ R ij t ] is introduced to record the total net revenues of the demand between each OD pair. Therefore, the total net revenues of the demand between the served OD pair ( i , j ) is R ij t = max { r ij e w ij , r ij w ij − C ij dl } .

Subtracting the operation costs of the hub network from the sum of net revenues of the demand satisfied between OD pairs, the objective value can be obtained. Operation costs C o of the hub network are the sum of the operating costs of the hubs and inter-hub links. They are represented as C o = ∑ h ∈ H f h + ∑ e ∈ E g e , where h is the hub in H , and e is the hub link in E . f and g are fixed costs of hub installation and hub link operation, respectively. Therefore, the objective value Obj of models without DCs is Obj = ∑ i , j ∈ N r ij e w ij − C o . The objective value of models allowing for DCs is Obj = ∑ i , j ∈ N R ij t − C o .

Initial Solution

This section describes how to gain a feasible initial solution. Flows that come in and out of node i are ( O i + D i ) . Nodes are sorted by the decreasing order of ( O i + D i ) . Then, the first | N | 4 nodes are selected as hubs and linked completely to determine the hub network G 0 , and the initial unit demand cost matrix C 0 h between hub pairs is obtained. After the initial hub network is determined, the assignments of the nodes are decided. For two variants of single allocation models, each node is allocated to its closest hub. For r-allocation models, each node is allocated to the first r nearest hubs. This paper does not present assignments of nodes in multiple allocation models. Finally, the objective value and revenue matrix R 0 e (models without DCs) or R 0 t (models allowing for DCs) are computed.

Neighborhood Structure

This section considers two types of local search strategies of removing/adding/replacing hubs and hub links in six problems. In multiple allocation problems, six neighborhood structures are defined for a given solution S . For single allocation problems, eight neighborhood structures are defined for a given solution S . In r-allocation problems, seven neighborhood structures are defined. Besides neighborhood structures of operating hub and hub links, an additional neighborhood structure is also defined to assign spokes for single allocation and r-allocation problems.

The first type of local search strategy is on hubs, including removing/adding/replacing hubs. For a given solution, the hub set is sorted by descending order of total flows through them, which is denoted by H . The spoke set is sorted by increasing order of total flows through them and is denoted by Sp . The objective value Ob j 0 is then calculated. In multiple allocation models, the first hub and the last hub of the path that each OD pair pass are determined. In single allocation and r-allocation models, clusters of hubs are computed.

Removing hubs: This neighborhood structure is denoted by Removing _ hub ( S ) . The hub h is removed and a new hub set H c is obtained. The current hub network is G c = ( H c , E c ) by linking the neighborhood nodes of h with each other. Then, nodes in the cluster of the hub h are reassigned in single allocation problems and r-allocation problems, and the new assignment is A c . For single allocation, nodes are assigned to the closest hub in G c . In r-allocation problems, links between nodes and the first r nearest hubs in H c are activated if they are not linked. Thus, the new solution S c is obtained and the new neighborhood is Removing _ hub ( S c ) .

Adding hubs: This neighborhood structure is denoted by Adding _ hub ( S ) . The node s is selected from the set Sp and added to the hub set, and the current hub set is H c . The new hub s is connected to its nearest hub in H , and the current hub network is obtained, that is, G c = ( H c , E c ) . Nodes are then assigned to the new hub s if c iA [ i ] > c is in single allocation problems. In r-allocation problems, the hub h m is replaced with the hub s in the assignment if min { c iA [ i ] [ j ] , j ∈ H , len ( A [ i ] ) < r } > c is , where h m is arg min { c iA [ i ] [ j ] , j ∈ H , len ( A [ i ] ) < r } . The new assignment is A c .

Replacing hubs: For multiple allocation problems and r-allocation problems, a neighborhood structure is denoted as Replacing _ hubs ( S ) . For single allocation problems, two neighborhood structures are applied to use clusters of hubs denoted by Alternate ( S ) and Locate ( S ) . For Replacing _ hubs ( S ) , the hub h is replaced by the spoke s . Similar to Removing _ hub ( S ) and Adding _ hub ( S ) , the neighborhood nodes of h in the hub set are connected with each other, and s is linked with its closest hub, thus the current hub network is G c = ( H c , E c ) . The assignments of nodes are also changed for single allocation and r-allocation problems, denoted by A c . In r-allocation models, links between nodes in the cluster h and the first r hubs in H c are activated if they are not linked, and the remaining nodes are reassigned in the same way as Adding _ hub ( S ) . In single allocation problems, the same neighborhood structures, called Alternate ( S ) and Locate ( S ) in Ilić et al. ( 57 ), are applied. Alternate ( S ) replaces hub h with a spoke node s in its cluster Cl u h , and the corresponding cluster is Cl u s . The assignments of nodes in Cl u s are changed to s . Locate ( S ) replaces hub h with a spoke node s in clusters of other hubs, and the nodes in the original cluster Cl u h are assigned to the closest hubs and the corresponding clusters and assignments are changed. The details of the neighborhood structures can be found in Ilić et al. ( 57 ).

The second type of local search strategy is on hub links, including removing/adding/replacing hub links. All the connected hub links E l are sorted in decreasing order of lengths and non-existent links E n are sorted in increasing order of lengths. The objective value is then computed for a given solution. In this type of local search strategy, the assignments of spokes are not changed.

Removing hub links: This neighborhood structure is denoted by Removing _ hublink ( S ) . A hub link hl is removed from the set E l and the new hub link set E c is obtained. The current hub network is G c = ( H , E c ) .

Adding hub links: This neighborhood structure is denoted by Adding _ hublink ( S ) . A hub link h l n is selected from the set E n and added to E l . The current hub link set is denoted by E c . The current hub network then is G c = ( H , E c ) .

Replacing hub links: This neighborhood structure is denoted by Replacing _ hublink ( S ) . The hub link hl in E c is removed and replaced with the hub link h l n in E n . The new hub link set and the hub network are then E c and G c = ( H , E c ) , respectively.

In the local search strategies of removing/adding/replacing hubs, spokes are assigned to hubs according to the nearest hub allocation rule in both single and r-allocation problems. To improve the performance of VNS, neighborhood structures of exchanging hubs and spokes are also used. In the single allocation, the operation called Allocate ( S ) is applied to change the assignments of spokes. This operation is proposed by Ilić et al. ( 57 ). For r-allocation problems, an operation of assigning spokes denoted by Assigning _ spokes ( S ) is introduced. In the current assignment of the spoke i , that is, A [ i ] , hubs in A [ i ] are exchanged with hubs in the complement set of A [ i ] . Then, the new assignment of the spoke i is A c [ i ] .

Shaking

The shaking procedure is used to avoid the local optimal solution, which requires the solution S 0 as the input. The current solution is replaced by a random solution chosen from the operations of removing/adding/replacing hubs. The main steps of the shaking operation are shown in Algorithm 1

OPEN IN VIEWER

Algorithm 1 Shaking procedure
Input: An initial solution S 0 Output: New solution S 1: Generate integer k from 0 to 2; 2: if k==0 then 3:   S = Removing _ hub ( S 0 ) ; 4: else if k==1 then 5:   S = Adding _ hub ( S 0 ) ; 6: else if k==2 then 7: S = Replacing _ hub ( S 0 ) ; 8: end if

Specifically, operations for hubs are selected randomly. The hub and the non-hub node are then randomly chosen to conduct the selected operation. A hub is randomly removed from the current hub set if the operation of removing hubs is selected. A non-hub node is randomly chosen and added to the set of hubs when the selected operation is adding hubs. A hub node is exchanged randomly with a non-hub node when the operation of replacing the hubs is selected. The output is the solution obtained after random operations on hubs.

Overall Heuristic

This section proposes heuristics based on VNS to solve HLPPs, which apply the proposed neighborhoods and the shaking procedure. A sequential strategy is used for exploration of the neighborhoods within deterministic local search, which is defined as Improvement _ procedure ( S ) .

In the improvement procedure, all l max neighborhoods that are used in the search are ordered, that is, l max = 6 for multiple allocation problems, l max = 8 for single allocation problems, and l max = 7 for r-allocation problems. Then, a local search is performed with the neighborhoods by the order. Once a better solution is found in any neighborhood, return to the first neighborhood. If the solution is not improved in the last neighborhood, then the procedure is terminated. The detailed steps of the improvement procedure are shown in Appendix A (Algorithm A1 for multiple allocation models, Algorithm A2 for single allocation models, and Algorithm A3 for r-allocation models).

The steps of VNS-based heuristics are presented in Algorithm 2. For a given initial solution with the hub set H 0 , the hub link set E 0 and assignments A 0 (for single allocation and r-allocation problems), the objective value is Ob j 0 . u is a parameter to record iterations without improvement. Let the output solution S be S 0 and objective value Obj be Ob j 0 . The shaking procedure is applied to the solution S and the current solution S c = ( H c , E c , A c ) ( S c = ( H c , E c ) for multiple allocation models) and the current objective value Ob j c are achieved. The improvement procedure is then performed on S c and the new solution S n and new objective value Ob j n are obtained. If S n is better than S , that is, Ob j n > Obj , then let S = S n , Obj = Ob j n and u = 0 . The algorithm is terminated once the condition is met.

OPEN IN VIEWER

Algorithm 2 Steps of VNS-based heuristics

Input: An initial solution S 0 , the initial objective value Ob j 0 ,   the maximum number of iterations ite r max ,   the maximum number of iterations without improvement ite r nimpr Output: The best solution S , the best net profit Obj 1: S = S 0 , Obj = Ob j 0 , u = 0 ; 2: S c = S 0 , Ob j c = Ob j 0 , k = 1 ; 3: while k ≤ ite r max do 4:  Apply ( S c , Ob j c ) = Shaking ( S ) ; 5:  Apply ( S n , Ob j n ) = Improvement _ procedure ( S c ) 6:  if Ob j n > Obj then 7:    Obj = Ob j n , S = S n , u = 0 8:  else 9:    u = u + 1 10:  end if 11:  if u ≥ ite r nimpr then 12:   return 13:  end if 14:   k = k + 1 15: end while

There are two parameters in the VNS-based heuristics: the maximum number of iterations, ite r max , and the maximum number of iterations without improvement, ite r nimpr . This study set ite r max = 2 | N | and ite r nimpr = | N | . The complete experiments of parameter settings can be found in the next section.

Datasets and Experimental Setups

This study uses datasets that are well known in HLPs: CAB with 25 nodes (CAB25) ( 13 ) and CAB with 100 nodes (CAB100) ( 62 ) to test the experiments. The CAB dataset with 25 nodes is based on airline passenger flows between 25 cities in the U.S.A., accessed in the OR Library ( 63 ). The transportation costs and travel demands between each city pair are provided. The CAB dataset with 100 nodes is the sample of air flows between major U.S. cities, and it is available from O’Kelly ( 62 ). In accordance with Taherkhani and Alumur ( 8 ), demands between each node pair are normalized to make the sum of demands equal to one. The parameter settings of unit revenue, discount of the transportation cost between hub pairs, cost of hub installation, cost of operating hub links, and cost of operating direct links are provided in Taherkhani and Alumur ( 8 ). The parameter settings of the CAB datasets are shown in Table 4.

OPEN IN VIEWER

Table 4.

Parameter Settings for CAB Datasets ( 8 )

Parameters	Descriptions	Values
| N |	Number of nodes	{ 10 , 15 , 20 , 25 , 60 }
w ij	Travel demands	CAB25: provided by OR Library ( 63 ); CAB100: provided by ( 62 )
c ij	Unit transportation cost	CAB25: provided by OR Library ( 63 ); CAB100: provided by O’Kelly ( 62 )
r ij	Revenue of unit demand	Low:1,000, Medium: 1,500, High: 2,000
f k	Cost of the hub installation	Low: 50, Medium: 100, High: 150
g kl	Fixed cost of setting up a hub link	0 . 1 f k
q ij	Fixed cost of setting up a direct link	0 . 2 g kl
α	Transportation cost discount factor	0 . 2 , 0 . 4 , 0 . 6 , 0 . 8

All experiments were carried out with one thread on a laptop with four-core Intel (R) Core (TM) i7-6500 2.50GHz processor and 16GB RAM under Fedora 30 environment. The VNS was implemented with the Python programming language. The exact solutions were obtained using IBM ILOG CPLEX 12.9 and Benders decomposition. The Benders decomposition algorithms are presented in Appendix C and Appendix D. The models were coded in Python with all CPLEX parameters being default. The tolerance gap of the optimal objective value is 10 − 5 gap. The maximum run time was set to 12 h.

Comparison of VNS-Based Heuristics With Exact Methods

This section presents tests of six variants of HLPPs for the CAB25 dataset with | N | ∈ { 10 , 15 , 20 , 25 } using the parameter settings in Table 4. There are tests for 36 instances for each subset with CPLEX, Benders decomposition, and VNS. Each instance was run with VNS 20 times and the median gaps between the objective value of VNS and the optimal solution computed. The initial solution of the VNS-based heuristics is determined by selecting the nodes ranking the first | N | 4 of the total flows as hub nodes and linking them to a complete network. In VNS-based heuristics, the maximum number of iterations was set to ite r max to 2 | N | and the maximum number of iterations without improvement ite r nimpr to | N | . The gap is computed by g a p = o b j o p t − o b j v n s o b j o p t * 100 % , where ob j vns is the objective value of the solution obtained by VNS, and ob j opt is the objective of the optimal solution by CPLEX (if CPLEX can achieve the optimal solution in the time limit) and the Benders decomposition approach.

Results for Multiple Allocation Problems

This section presents the tests of two variants of multiple allocation problems. For UMAI, 139 of 144 instances are solved optimally by VNS for median gaps. Four instances use up the time limit, and three of them do not find the optimal solution with CPLEX in 12 h. For UMAIDC, VNS solves 134 of 144 instances to optimality for the median gaps. In the experiments, CPLEX solves three instances using up the time limit, and one instance is not solved to optimality.

It is more difficult to solve the instances with high revenue than those with medium and low revenues (see Taherkhani and Alumur [ 8 ]). Therefore, Table 5 summarizes the performance of CPLEX, Benders, and VNS for instances with high revenue on multiple allocation models. The results of instances with medium and low revenues are also shown in Appendix B. Two letters are used to represent the cost and revenue parameters. The first letter represents the cost parameter: “L” is the low cost with a value of 50, “M” is the medium cost with a value of 100, and “H” is the high cost with a value of 150. The second letter represents the revenue parameter. The first three columns report parameters of costs and revenues, the problem size, and the discount factor. For each instance in multiple allocation problems, the next two columns show the objective value of CPLEX obtained within the time limit, and the run time of instances. CPLEX uses up the time limit for some instances, and “*” is used to mark them in the column “Time(s).” The next two columns show the objective value and the run time solved by Benders, followed by the median gaps of VNS and the average run time of VNS.

OPEN IN VIEWER

Table 5.

Multiple Allocation Problems Solutions With High Revenue

| N |	Pra.	α	Multiple allocation model						Multiple allocation model allowing for direct connections
			CPLEX		Benders		VNS		CPLEX		Benders		VNS
			Obj	Time (s)	Obj	Time (s)	Median gap (%)	Average time (s)	Obj	Time (s)	Obj	Time (s)	Median gap (%)	Average time (s)
10	L	0.2	1,397.27	17.43	1,397.27	12.63	0.00	19.01	1,398.76	35.08	1,398.78	12.40	0.00	12.12
		0.4	1,300.20	21.52	1,300.20	5.99	0.00	6.18	1,306.68	46.56	1,306.68	3.95	0.00	3.26
		0.6	1,240.33	3.67	1,240.33	5.00	0.00	3.52	1,290.93	0.63	1,290.93	2.10	0.00	0.25
		0.8	1,200.58	14.70	1,200.58	4.89	0.00	2.29	1,290.93	0.64	1,290.93	1.18	0.00	0.25
	M	0.2	1,170.69	3.97	1,170.69	6.85	0.00	2.29	1,205.22	8.24	1,205.22	3.60	0.00	0.36
		0.4	1,117.89	4.30	1,117.89	4.86	0.00	1.23	1,205.22	0.70	1,205.22	2.25	0.00	0.28
		0.6	1,081.31	8.96	1,081.31	4.87	0.00	1.68	1,205.22	0.60	1,205.22	1.93	0.00	0.26
		0.8	1,078.80	1.01	1,078.80	3.54	0.00	1.57	1,205.22	0.43	1,205.22	1.18	0.00	0.27
	H	0.2	1,057.88	0.57	1,057.88	3.97	0.00	1.37	1,127.95	0.63	1,127.95	3.25	0.00	0.29
		0.4	1,007.89	1.08	1,007.89	4.52	0.00	1.34	1,127.95	0.53	1,127.95	2.95	0.00	0.24
		0.6	978.80	1.44	978.80	3.94	0.00	1.22	1,127.95	0.57	1,127.95	2.23	0.00	0.25
		0.8	978.80	0.48	978.80	3.87	0.00	1.53	1,127.95	0.48	1,127.95	1.14	0.00	0.23
15	L	0.2	1,177.66	2,045.15	1,177.67	146.59	0.20	40.08	1,188.42	1,221.73	1,188.42	64.72	0.00	32.48
		0.4	1,037.80	696.73	1,037.81	48.22	0.00	42.86	1,051.26	1,063.91	1,051.27	41.50	0.01	11.68
		0.6	916.13	1,963.52	916.14	53.41	0.00	38.94	966.48	90.26	966.48	27.29	0.00	5.59
		0.8	843.21	1,174.31	843.21	38.38	0.00	18.15	942.21	5.76	942.21	21.52	0.00	2.65
	M	0.2	900.21	290.77	900.21	52.78	0.00	24.20	925.39	77.38	925.39	47.04	0.00	5.72
		0.4	780.58	955.92	780.58	43.77	0.00	11.48	862.17	36.41	862.17	34.55	0.00	4.81
		0.6	706.42	636.12	706.42	36.36	0.00	6.63	832.79	5.02	832.79	20.42	0.00	1.99
		0.8	706.42	6.09	706.42	37.05	0.00	5.68	832.79	3.10	832.79	18.20	0.00	2.36
	H	0.2	701.40	52.44	701.40	42.84	0.00	5.34	780.80	5.62	780.80	29.08	0.00	4.32
		0.4	638.18	20.03	638.18	39.34	0.00	5.30	738.20	3.40	738.20	23.47	0.00	2.92
		0.6	617.26	3.67	617.26	33.85	0.00	2.80	738.20	2.72	738.20	15.94	0.00	2.86
		0.8	617.26	4.34	617.26	25.83	0.00	3.05	738.20	2.59	738.20	17.12	0.00	2.90
20	L	0.2	1,235.80	18,005.96	1,235.80	241.64	0.00	85.89	1,237.34	16,536.04	1,237.35	313.82	0.00	68.85
		0.4	1,089.63	26,798.68	1,089.63	235.10	0.00	52.50	1,097.03	31,019.81	1,097.04	203.17	0.00	40.21
		0.6	967.81	19,195.13	967.81	221.14	0.00	60.21	993.72	7,496.58	993.73	214.90	0.00	22.28
		0.8	908.67	4,573.99	908.67	183.94	0.00	26.88	965.95	333.65	965.95	111.54	0.00	20.28
	M	0.2	958.90	3,443.96	958.91	177.21	0.00	65.74	979.82	1,285.79	979.83	174.23	0.00	37.08
		0.4	856.59	523.58	856.60	97.01	0.00	32.49	884.20	371.37	884.21	135.95	0.00	18.28
		0.6	787.47	186.75	787.47	150.56	0.00	19.26	850.74	64.41	850.74	164.65	0.01	14.21
		0.8	773.38	57.64	773.38	106.42	0.00	26.11	850.74	41.76	850.74	81.30	0.00	9.32
	H	0.2	778.33	174.64	778.33	120.81	0.00	36.34	797.80	273.68	797.81	143.99	0.00	20.39
		0.4	703.12	166.79	703.12	115.20	0.00	18.65	765.95	26.69	765.95	154.26	0.00	9.42
		0.6	677.47	370.55	677.47	75.82	0.00	13.90	765.95	36.22	765.95	93.90	0.00	7.77
		0.8	673.38	284.47	673.38	100.11	0.00	24.63	765.95	35.78	765.95	69.14	0.00	9.72
25	L	0.2	1,162.92	*	1,162.92	858.99	0.00	250.24	1,166.32	*	1,166.32	893.92	0.00	155.19
		0.4	1,008.43	*	1,008.46	1,087.39	0.00	179.21	1,024.43	*	1,024.43	842.47	0.00	89.75
		0.6	898.06	*	898.24	900.95	0.00	115.65	926.48	23,270.90	926.49	592.33	0.00	57.00
		0.8	839.40	26,435.50	839.40	582.50	0.00	80.50	898.47	*	901.10	395.40	0.00	61.50
	M	0.2	911.26	29,504.19	911.27	851.93	0.00	66.07	927.43	2,317.84	927.44	628.26	0.00	52.23
		0.4	803.72	1,551.94	803.73	435.42	0.00	32.35	823.75	1,608.95	823.75	422.80	0.00	51.00
		0.6	717.73	7,584.62	717.73	561.17	0.00	30.15	777.74	723.23	777.74	298.47	0.00	43.95
		0.8	683.21	*	690.90	439.87	0.00	28.91	777.74	143.40	777.74	370.60	0.00	23.68
	H	0.2	738.08	1,294.94	738.08	650.32	0.00	28.35	748.03	3,756.20	748.04	436.89	0.00	40.05
		0.4	633.72	4,853.83	633.73	391.88	0.00	26.53	694.52	149.31	694.52	471.41	0.00	18.49
		0.6	599.18	728.05	599.18	349.44	0.00	36.64	694.52	177.14	694.52	241.64	0.00	14.92
		0.8	599.18	218.24	599.18	236.08	0.00	34.56	694.52	145.86	694.52	141.94	0.00	14.06

Note: |N| = the node number of instances; Pra. = the cost parameter; L = the low cost with a value of 50; M = the medium cost with a value of 100, H = the high cost with a value of 150. α = the discount factor with {0.2, 0.4, 0.6, 0.8}; Obj =the objective value; * = CPLEX uses up the time limit.

From Table 5, it is found that 43 of 48 instances can obtain the optimal solution with VNS for both multiple allocation problems. The run time of VNS is much faster than CPLEX and Benders. The model with DCs results in higher net profits compared with the one without DCs. When hub installation costs increase, values of net profit decrease, and the run time of CPLEX, Benders, and VNS tends to decrease. When discount factor α increases with the same values of installation cost and revenue, the net profit value decreases, and the run time of CPLEX, Benders, and VNS also tends to decrease.

Results for Single Allocation Problems

The CAB25 dataset was also tested with node number | N | ∈ { 10 , 15 , 20 , 25 } for single allocation problems. In the experiments, 134 of 144 instances of USAI obtain the optimal solutions by VNS. Six instances solved by CPLEX take more than 12 h, and two of them are solved with gaps. In the variant with DCs, 142 of 144 instances are solved to optimality by VNS. CPLEX solves three instances using up the time limit, but all of them are solved optimally.

Table 6 summarizes solutions of single allocation problems with high revenue. The results of instances with medium and low revenues for single allocation problems are reported in Appendix B. The same columns are used in Table 6 as those in Table 5. “*” is also used to mark the instances which are solved by CPLEX using up the time limit in the column “Time(s)” of CPLEX. Observing Table 6, it is found that the model with DCs results in higher net profits compared with the variant without DCs. The VNS performs better for the model allowing for DCs than the one without DCs for high revenue. In the model without DCs, the solutions obtained by VNS have gaps in seven instances, while 48 instances with the high revenue of the model allowing for DCs all obtain the optimal solution.

OPEN IN VIEWER

Table 6.

Single Allocation Problems Solutions With High Revenue

| N |	Pra.	α	Single allocation model						Single allocation model allowing for direct connections
			CPLEX		Benders		VNS		CPLEX		Benders		VNS
			Obj	Time (s)	Obj	Time (s)	Median gap (%)	Average time (s)	Obj	Time (s)	Obj	Time (s)	Median gap (%)	Average time (s)
10	L	0.2	1,394.05	11.68	1,394.05	17.62	0.00	2.21	1,394.05	14.19	1,394.05	18.12	0.00	2.12
		0.4	1,274.03	9.95	1,274.03	9.56	0.45	1.13	1,290.93	17.53	1,290.93	14.16	0.00	0.44
		0.6	1,176.89	16.98	1,176.89	6.03	0.09	0.59	1,290.93	4.52	1,290.93	3.71	0.00	0.28
		0.8	1,110.33	13.88	1,110.33	4.61	0.00	0.64	1,290.93	4.30	1,290.93	1.70	0.00	0.33
	M	0.2	1,164.99	7.88	1,164.99	5.27	0.00	0.48	1,205.22	11.10	1,205.22	4.67	0.00	0.31
		0.4	1,105.69	6.18	1,105.69	3.91	0.00	0.50	1,205.22	3.90	1,205.22	2.79	0.00	0.35
		0.6	1,047.37	7.71	1,047.37	4.56	0.00	0.46	1,205.22	3.26	1,205.22	2.08	0.00	0.36
		0.8	990.12	9.16	990.12	5.45	0.00	0.43	1,205.22	4.20	1,205.22	1.79	0.00	0.33
	H	0.2	1,054.01	4.52	1,054.01	3.05	0.00	0.46	1,127.95	3.00	1,127.95	3.37	0.00	0.35
		0.4	995.69	3.21	995.69	2.85	0.00	0.34	1,127.95	4.11	1,127.95	2.62	0.00	0.34
		0.6	937.37	4.89	937.37	3.87	0.00	0.39	1,127.95	2.86	1,127.95	2.61	0.00	0.34
		0.8	918.95	3.15	918.95	2.82	0.00	0.21	1,127.95	2.87	1,127.95	1.75	0.00	0.31
15	L	0.2	1,165.51	3,383.75	1,165.51	159.87	0.43	4.41	1,175.39	4,655.54	1,175.39	125.49	0.00	4.24
		0.4	992.30	2,827.97	992.30	75.41	0.05	5.26	1,011.83	2,904.09	1,011.83	82.46	0.00	2.54
		0.6	840.94	2,877.65	840.94	78.24	0.00	3.97	942.21	142.57	942.21	30.61	0.00	1.46
		0.8	742.03	445.89	742.03	28.87	0.97	1.48	942.21	24.29	942.21	14.11	0.00	1.33
	M	0.2	885.40	107.95	885.40	44.86	0.00	2.90	914.76	144.84	914.76	31.20	0.00	1.77
		0.4	736.29	242.36	736.29	52.55	0.33	2.20	833.41	101.28	833.41	25.52	0.00	1.41
		0.6	667.26	130.86	667.26	24.08	0.00	0.80	832.79	35.49	832.79	14.27	0.00	1.16
		0.8	667.26	38.83	667.26	16.37	0.00	0.70	832.79	33.66	832.79	11.17	0.00	1.12
	H	0.2	690.76	96.68	690.76	17.77	0.00	1.43	770.17	45.90	770.17	22.65	0.00	1.88
		0.4	617.26	67.30	617.26	16.29	0.00	0.78	738.20	37.40	738.20	15.25	0.00	1.28
		0.6	617.26	26.75	617.26	12.09	0.00	0.93	738.20	21.51	738.20	11.16	0.00	1.18
		0.8	617.26	30.22	617.26	10.90	0.00	0.77	738.20	21.15	738.20	14.86	0.00	1.25
20	L	0.2	1,221.92	8,446.00	1,221.92	245.80	0.00	10.81	1,224.47	17,868.13	1,224.47	352.60	0.00	14.36
		0.4	1,050.44	*	1,050.44	296.88	0.00	8.07	1,058.96	*	1,058.96	355.25	0.00	10.22
		0.6	883.76	*	886.63	353.89	0.00	6.39	965.95	2,598.66	965.95	150.03	0.00	3.17
		0.8	795.66	8,543.51	795.66	217.96	0.00	5.16	965.95	265.00	965.95	95.71	0.00	3.14
	M	0.2	944.05	1,033.23	944.05	213.15	0.00	6.77	964.42	880.61	964.42	147.96	0.00	5.54
		0.4	815.10	1,185.41	815.10	102.82	0.00	4.65	850.74	1,641.09	850.74	163.33	0.00	3.46
		0.6	736.23	897.22	736.23	134.72	0.00	3.47	850.74	216.18	850.74	85.45	0.00	3.13
		0.8	712.83	822.65	712.83	69.27	0.00	2.66	850.74	205.00	850.74	41.48	0.00	3.16
	H	0.2	96.55	36.08	766.37	88.39	0.00	4.88	786.29	896.34	786.29	164.23	0.00	6.36
		0.4	677.86	767.96	677.86	66.78	0.00	3.15	765.95	187.96	765.95	88.51	0.00	2.83
		0.6	662.83	265.07	662.83	55.06	0.00	2.15	765.95	218.44	765.95	50.35	0.00	3.48
		0.8	662.83	208.50	662.83	35.68	0.00	2.67	765.95	147.40	765.95	46.63	0.00	3.87
25	L	0.2	1,149.28	*	1,149.28	594.54	0.00	20.51	1,153.19	*	1,153.19	996.82	0.00	17.25
		0.4	975.52	*	975.52	965.45	0.00	19.43	991.14	*	991.14	903.83	0.00	19.15
		0.6	830.82	*	830.82	606.46	1.58	10.48	888.19	15,826.47	888.19	494.38	0.00	6.56
		0.8	736.42	*	736.42	611.83	0.90	8.98	888.19	1,002.34	888.19	275.12	0.00	6.91
	M	0.2	895.55	10,885.93	895.55	692.85	0.00	12.27	915.20	7,030.10	915.20	477.96	0.00	12.88
		0.4	765.00	7,003.91	765.00	276.16	0.00	7.80	787.35	12,709.38	787.35	392.63	0.00	11.75
		0.6	649.18	13,204.39	649.18	400.75	0.00	5.47	777.74	1,607.73	777.74	283.72	0.00	5.43
		0.8	649.18	1,500.68	649.18	186.76	0.00	5.50	777.74	1,402.66	777.74	151.92	0.00	5.99
	H	0.2	725.55	4,074.29	725.55	281.47	0.00	7.12	735.80	5,344.98	735.80	319.45	0.00	10.53
		0.4	599.18	5,891.26	599.18	311.78	0.00	6.55	694.52	1,682.25	694.52	243.43	0.00	5.47
		0.6	599.18	1,776.76	599.18	128.62	0.00	5.06	694.52	1,165.25	694.52	118.09	0.00	6.52
		0.8	599.18	2,624.90	599.18	120.10	0.00	5.20	694.52	872.34	694.52	126.59	0.00	4.56

Note: |N| = the node number of instances; Pra. = the cost parameter; L = the low cost with a value of 50; M = the medium cost with a value of 100, H = the high cost with a value of 150. α = the discount factor with {0.2, 0.4, 0.6, 0.8}; Obj =the objective value; * = CPLEX uses up the time limit.

Compared with the solutions in Table 5, the net profit values of single allocation models are lower than those of multiple allocation models for instances with the same parameter settings. CPLEX solves instances with multiple allocation models more easily than those with the single allocation strategy, while VNS solves single allocation models faster than multiple allocation models because of the larger space to search in multiple allocation problems.

Results for r-Allocation Problems

To analyze the performance of r-allocation models, the same CAB25 instances are tested with node number | N | ∈ { 10 , 15 , 20 , 25 } , and r = 2 . In the experiments, 132 of 144 instances of the r-allocation problem without DCs obtain the optimal solutions by VNS. CPLEX solves 14 instances using up the time limit, and 12 of them are solved with gaps. In the model with DCs, 130 of 144 instances are solved to optimality by VNS. When solved by CPLEX, five instances use up the time limit, and two of them are solved with gaps.

Table 7 reports the results obtained with r-allocation models with high revenue. The results of instances with medium and low revenues for r-allocation problems are shown in Appendix B. In the model without DCs, it is found that eight of 48 instances are solved by VNS with non-zero gaps, and only one instance has a gap of over 1%. CPLEX uses up the time limit for nine of 48 instances and obtains a solution with non-zero gaps for seven of them. In the model allowing for DCs, all of 48 instances are solved to optimality by VNS and four of them use up the time limit when solved by CPLEX.

OPEN IN VIEWER

Table 7.

r-Allocation Problems Solutions With High Revenue

| N |	Pra.	α	r-Allocation model						r-Allocation model allowing for direct connections
			CPLEX		Benders		VNS		CPLEX		Benders		VNS
			Obj	Time (s)	Obj	Time (s)	Median gap (%)	Average time (s)	Obj	Time (s)	Obj	Time (s)	Median gap (%)	Average time (s)
10	L	0.2	1,397.26	27.92	1,397.26	10.74	0.00	3.97	1,397.69	37.28	1,397.69	13.02	0.00	2.12
		0.4	1,296.72	31.28	1,296.72	7.02	0.00	2.35	1,303.46	74.79	1,303.46	8.64	0.00	0.44
		0.6	1,234.36	14.52	1,234.36	6.61	0.00	1.80	1,290.93	0.91	1,290.93	4.91	0.00	0.28
		0.8	1,191.18	27.77	1,191.18	8.16	1.04	1.27	1,290.93	0.93	1,290.93	2.63	0.00	0.33
	M	0.2	1,169.71	6.91	1,169.71	6.82	0.00	1.45	1,205.22	9.02	1,205.22	4.95	0.00	0.31
		0.4	1,117.89	4.62	1,117.89	7.44	0.00	0.66	1,205.22	0.97	1,205.22	4.13	0.00	0.35
		0.6	1,081.31	9.61	1,081.31	6.89	0.00	0.87	1,205.22	0.84	1,205.22	3.74	0.00	0.36
		0.8	1,078.80	1.39	1,078.80	4.44	0.00	1.06	1,205.22	0.88	1,205.22	1.92	0.00	0.33
	H	0.2	1,057.88	0.67	1,057.88	4.63	0.00	1.08	1,127.95	0.76	1,127.95	3.95	0.00	0.35
		0.4	1,007.89	1.03	1,007.89	4.43	0.00	0.87	1,127.95	0.76	1,127.95	3.86	0.00	0.34
		0.6	978.80	3.19	978.80	5.00	0.00	0.98	1,127.95	0.69	1,127.95	3.53	0.00	0.34
		0.8	978.80	1.18	978.80	4.22	0.00	1.25	1,127.95	0.82	1,127.95	2.93	0.00	0.31
15	L	0.2	1,174.89	5,287.78	1,174.89	137.33	0.07	8.19	1,187.22	2,467.73	1,187.23	69.22	0.00	4.24
		0.4	1,025.77	3,625.79	1,025.78	66.17	0.31	11.80	1,042.36	3,376.71	1,042.36	75.90	0.00	2.54
		0.6	892.78	13,621.05	892.79	83.95	0.03	7.32	966.48	114.52	966.48	43.36	0.00	1.46
		0.8	817.62	18,436.48	817.63	114.25	0.24	3.78	942.21	10.90	942.21	13.48	0.00	1.33
	M	0.2	899.00	436.67	899.01	47.46	0.00	6.58	925.39	96.37	925.39	52.95	0.00	1.77
		0.4	777.15	1,206.99	777.16	47.91	0.00	5.44	862.17	41.96	862.17	30.34	0.00	1.41
		0.6	706.42	2,896.10	706.42	53.20	0.00	6.09	832.79	7.82	832.79	16.34	0.00	1.16
		0.8	706.42	32.35	706.42	19.83	0.00	4.45	832.79	7.56	832.79	14.35	0.00	1.12
	H	0.2	701.40	88.75	701.40	53.51	0.00	2.83	780.80	9.34	780.80	25.81	0.00	1.88
		0.4	638.18	34.63	638.18	31.66	0.00	2.77	738.20	10.57	738.20	17.25	0.00	1.28
		0.6	617.26	8.14	617.26	16.14	0.00	1.84	738.20	7.34	738.20	8.83	0.00	1.18
		0.8	617.26	8.66	617.26	24.61	0.00	1.52	738.20	5.87	738.20	8.70	0.00	1.25
20	L	0.2	1,235.30	7,327.84	1,235.30	260.18	0.00	18.95	1,237.00	21,000.78	1,237.00	304.21	0.00	14.36
		0.4	1,081.56	*	1,081.48	259.32	0.00	22.62	1,091.17	*	1,091.17	323.56	0.00	10.22
		0.6	952.91	*	954.37	367.82	0.84	13.95	992.98	9,033.65	992.99	210.21	0.00	3.17
		0.8	889.45	*	890.86	154.09	0.77	16.36	965.95	736.21	965.95	135.54	0.00	3.14
	M	0.2	958.90	4,512.37	958.91	288.55	0.00	15.95	979.83	2,343.45	979.83	182.69	0.00	5.54
		0.4	856.59	1,711.02	856.60	110.45	0.00	11.96	884.20	2,110.64	884.21	151.51	0.00	3.46
		0.6	787.47	642.42	787.47	153.58	0.00	6.12	850.74	379.99	850.74	144.81	0.00	3.13
		0.8	773.38	939.92	773.38	102.39	0.00	11.30	850.74	131.84	850.74	70.73	0.00	3.16
	H	0.2	778.33	379.11	778.33	137.53	0.00	14.66	797.80	642.89	797.81	137.36	0.00	6.36
		0.4	703.12	693.83	703.12	67.47	0.00	5.76	765.95	247.58	765.95	107.71	0.00	2.83
		0.6	677.47	679.01	677.47	119.97	0.00	5.61	765.95	167.47	765.95	57.03	0.00	3.48
		0.8	673.38	31,646.98	673.38	65.26	0.00	11.81	765.95	115.06	765.95	59.87	0.00	3.87
25	L	0.2	1,162.57	*	1,162.58	809.31	0.00	50.46	1,165.98	*	1,165.98	1,226.61	0.00	17.25
		0.4	1,002.23	*	1,002.24	1,389.84	0.00	51.87	1,019.23	*	1,019.24	1,370.07	0.00	19.15
		0.6	880.16	*	887.23	1,432.92	0.00	27.26	925.99	25,747.41	925.99	506.06	0.00	6.56
		0.8	829.03	*	832.10	487.68	0.20	39.04	898.25	*	901.10	667.65	0.00	6.91
	M	0.2	911.26	24,059.11	911.27	756.85	0.00	37.39	927.43	4,052.59	927.44	787.85	0.00	12.88
		0.4	803.72	12,341.73	803.73	488.48	0.00	20.44	823.74	4,185.33	823.75	384.57	0.00	11.75
		0.6	717.23	*	717.23	430.82	0.00	26.57	777.74	2,725.07	777.74	363.97	0.00	5.43
		0.8	683.21	*	690.90	380.68	0.00	21.79	777.74	983.35	777.74	233.53	0.00	5.99
	H	0.2	738.08	3,774.87	738.08	505.12	0.00	20.85	748.03	11,142.71	748.04	421.78	0.00	10.53
		0.4	633.72	15,855.35	633.73	271.21	0.00	19.33	694.52	1,528.71	694.52	270.92	0.00	5.47
		0.6	599.18	3,073.49	599.18	347.41	0.00	30.27	694.52	1,090.49	694.52	130.25	0.00	6.52
		0.8	599.18	1,447.71	599.18	272.94	0.00	26.81	694.52	973.00	694.52	83.53	0.00	4.56

Note: r = 2; |N| = the node number of instances; Pra. = the cost parameter; L = the low cost with a value of 50; M = the medium cost with a value of 100, H = the high cost with a value of 150. α = the discount factor with {0.2, 0.4, 0.6, 0.8}; Obj =the objective value; * = CPLEX uses up the time limit.

The r-allocation models are more difficult to solve than the single and multiple allocation models when CPLEX is used. Compared with multiple allocation models, additional sets of binary variables and corresponding sets of constraints need to be defined. They are more difficult than single allocation models with regard to solving efficiency, because allowing for r hubs instead of one hub leads to larger solution space. When solved by VNS, it requires more CPU time than the single allocation models, but less than the multiple allocation models because the search spaces of r-allocation problems are smaller than multiple allocation versions but larger than single allocation variants.

To show the performance of the heuristics for solving different variants of HLPPs, Figure 1 shows the average time of CAB25 instances solved by VNS, CPLEX, and enhanced Benders decomposition for different network sizes. The average run time of VNS, CPLEX, and enhanced Benders decomposition for different network sizes is visualized in each sub-figure. The y-axis of each figure is on a log scale. As shown in Figure 1, the VNS of multiple allocation problems is two orders of magnitude faster than CPLEX, while it is one order of magnitude faster than that of Benders. The VNS of single allocation problems is two to three orders of magnitude faster than CPLEX, while it is one to two orders of magnitude faster than that of Benders. The VNS of r-allocation problems is two to three orders of magnitude faster than CPLEX, and it is one order of magnitude faster than that of Benders. Figure 2 shows the average of gaps of CAB25 instances with the node number | N | ∈ { 10 , 15 , 20 , 25 } . Minimum gaps, maximum gaps, and median gaps are shown. For each size network, the average of median gaps are less than 0.5%. Figures 3 and 4 report the number of instances solved by three methods to optimality in 1 h and in all tests. Except for the instances that cannot be solved to optimality, VNS solves over 90% of instances to optimality in 1 h.

OPEN IN VIEWER

Figure 1.

Average run time of instances solved by VNS, CPLEX, and enhanced Benders decomposition with the same number of nodes. The number of nodes is | N | ∈ { 10 , 15 , 20 , 25 } : (a) UMAI model, (b) UMAIDC model, (c) USAI model, (d) USAIDC model, (e) UrAI model, and (f) UrAIDC model; U = uncapacitated problem; MA = multiple allocation; SA = single allocation; rA = r-allocation; I = the problem with incomplete hub network.

OPEN IN VIEWER

Figure 2.

The average gaps of variable neighborhood search solutions compared with optimal solutions for instances with the same number of nodes. The number of nodes is | N | ∈ { 10 , 15 , 20 , 25 } : (a) UMAI model, (b) UMAIDC model, (c) USAI model, (d) USAIDC model, (e) UrAI model, and (f) UrAIDC model; U = uncapacitated problem; MA = multiple allocation; SA = single allocation; rA = r-allocation; I = the problem with incomplete hub network.

OPEN IN VIEWER

Figure 3.

Numbers of instances solved by three methods to optimality in one hour. The number of nodes is | N | ∈ { 10 , 15 , 20 , 25 } : (a) UMAI model, (b) UMAIDC model, (c) USAI model, (d) USAIDC model, (e) UrAI model, and (f) UrAIDC model; U = uncapacitated problem; MA = multiple allocation; SA = single allocation; rA = r-allocation; I = the problem with incomplete hub network.

OPEN IN VIEWER

Figure 4.

Numbers of instances solved by three methods to optimality in all tests. The number of nodes is | N | ∈ { 10 , 15 , 20 , 25 } : (a) UMAI model, (b) UMAIDC model, (c) USAI model, (d) USAIDC model, (e) UrAI model, and (f) UrAIDC model; U = uncapacitated problem; MA = multiple allocation; SA = single allocation; rA = r-allocation; I = the problem with incomplete hub network.

Figures 5 and 6 show the run time of CPLEX and VNS for instances with different parameter settings and node number | N | ∈ { 20 , 25 } . Because it is more difficult to solve the instances with high revenue than those with medium and low revenues (see Taherkhani and Alumur [8]), the results with high revenue are reported first, in Figure 5. It is found that the run times of CPLEX and VNS have a positive correlation for most instances. The instance which is difficult for CPLEX to solve also uses more time when solved by VNS. Instances with lower values of the discount factor tend to take a longer solution time given the fixed network scale and hub installation cost. When the network scale and the value of the discount factor are the same, solving the instances with the lower cost value of the hub installation takes a longer time. Figure 6 shows the run time of CPLEX and VNS for the instances with low hub installation costs. Figure 6 permits similar conclusions to those from Figure 5.

OPEN IN VIEWER

Figure 5.

Run time of instances with high revenue solved by CPLEX and variable neighborhood search (VNS). The revenue is 2,000; the number of nodes is | N | ∈ { 20 , 25 } . The x-axis are instances with different costs of hub installation and discount factors, where the letter represents the value of hub installation costs: L = 50, M = 100, H = 150, and discount factor = { 0 . 2 , 0 . 4 , 0 . 6 , 0 . 8 } . The y-axis is run time in log scale. “CPLEX_20” is the solution obtained by CPLEX for the instances with 20 nodes; “VNS_20” is the solution obtained by VNS for instances with 20 nodes; “CPLEX_25” is the solution obtained by CPLEX for instances with 25 nodes; “VNS_25” is the solution obtained by VNS for instances with 25 nodes: (a) UMAI model, (b) UMAIDC model, (c) USAI model, (d) USAIDIodel, (e) UrAI model, and (f) UrAIDC model; U = uncapacitated problem; MA = multiple allocation; SA = single allocation; rA = r-allocation; I = the problem with incomplete hub network.

OPEN IN VIEWER

Figure 6.

Run time of instances with low hub installation costs solved by CPLEX and variable neighborhood search. The number of nodes is | N | ∈ { 20 , 25 } . The x-axis are the instances with different values of revenues and discount factors; the letter represents the value of revenues: L = 1,000, M = 1,500; H = 2,000, and discount factor = { 0 . 2 , 0 . 4 , 0 . 6 , 0 . 8 } . The y-axis is run time in log scale. “CPLEX_20” is the solution obtained by CPLEX for the instances with 20 nodes; “VNS_20” is the solution obtained by VNS for instances with 20 nodes; “CPLEX_25” is the solution obtained by CPLEX for instances with 25 nodes; “VNS_25” is the solution obtained by VNS for instances with 25 nodes: (a) UMAI model, (b) UMAIDC model, (c) USAI model, (d) USIC model, (e) UrAI model, and (f) UrAIDC model; U = uncapacitated problem; MA = multiple allocation; SA = single allocation; rA = r-allocation; I = the problem with incomplete hub network.

Discussion of VNS-Based Heuristics Setups

Initial Solution

This section compares the determined initial solution with random initial solutions. The determined initial solution is generated by picking up nodes ranking the first | N | 4 of total flows as hub nodes and linking them completely. For the random solution, the hubs are select randomly from the node set, and then linked completely. The CAB25 instance is tested with the node number | N | = 25 , the high revenue (2,000), the low hub installation cost (50), and the discount factor α = 0 . 2 . Twenty random initial solutions are tested for each HLP and the median objective value of the solution is obtained to compare it with the solution obtained by the VNS starting with Ihe determined initial solution. The reIults are shown in Table 8. The column “Time (s)” is the average time of 20 tests. The tested instance can be solved to optimality with both the determined initial solution and the random initial solution for six models. The average time of the case with the determined initial solution is longer than that of the case with a random initial solution for models without DCs, while it is the opposite for models allowing for DCs. Therefore, the heuristics are not sensitive to the initial solutions.

OPEN IN VIEWER

Table 8.

Comparison of Determined Initial Solution and Random Initial Solution

Problems	Determined			Random
	Maximum gaps (%)	Median gaps (%)	Time (s)	Maximum gaps (%)	Median gaps (%)	Time (s)
UMAI	0.00	0.00	250.24	0.00	0.00	210.59
UMAIDC	0.00	0.00	155.19	0.00	0.00	170.47
USAI	0.00	0.00	20.51	0.00	0.00	12.63
USAIDC	0.00	0.00	17.25	0.00	0.00	20.00
UrAI	0.00	0.00	50.46	0.00	0.00	49.95
UrAIDC	0.00	0.00	36.58	0.00	0.00	45.68

Note: Revenue = 2,000; hub installation cost = 50; α = 0 . 2 ; number of nodes = 25; U = uncapacitated problem; MA = multiple allocation; SA = single allocation; rA = r-allocation; I = the problem with incomplete hub network.