An efficient matheuristic for the robust multiple allocation p-hub median problem under polyhedral demand uncertainty

Accepted Manuscript

An efficient matheuristic for the robust multiple allocation p-hub median problem under polyhedral demand uncertainty Nader Ghaffarinasab PII: DOI: Reference:

S0305-0548(18)30112-6 10.1016/j.cor.2018.04.021 CAOR 4462

To appear in:

Computers and Operations Research

Received date: Revised date: Accepted date:

12 August 2017 12 March 2018 30 April 2018

Please cite this article as: Nader Ghaffarinasab, An efficient matheuristic for the robust multiple allocation p-hub median problem under polyhedral demand uncertainty, Computers and Operations Research (2018), doi: 10.1016/j.cor.2018.04.021

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • Robust multiple allocation p-hub median problems under polyhedral uncertainty are addressed. • Three different uncertainty models, namely the hose, the hybrid, and the budget models are considered. • Linear MIP formulations are proposed for the problems.

CR IP T

• An efficient tabu search (TS) based matheuristic is proposed for solving the problems.

AC

CE

PT

ED

M

AN US

• Extensive computational experiments are conducted for the considered problems.

1

ACCEPTED MANUSCRIPT

An efficient matheuristic for the robust multiple allocation p-hub median problem under polyhedral demand uncertainty Nader Ghaffarinasaba,∗ of Industrial Engineering, University of Tabriz, Tabriz, Iran

CR IP T

a Department

Abstract

This paper addresses the robust multiple allocation p-hub median problem under polyhedral demand uncertainty. Three variants of polyhedral uncertainty models are used in the problem, namely the hose, the hybrid, and the budget uncertainty models. The problems are formulated as linear mixed integer programming problems and a Tabu Search (TS) based matheuristic approach is proposed to solve the three variants of the

AN US

problem. Extensive computational experiments are conducted based on three well-known data sets in the hub location literature and the results show the capability of the proposed solution algorithm to obtain the optimal solutions in short computational times.

Keywords: Hub location, polyhedral uncertainty, mathematical formulation, matheuristics, tabu search.

M

1. Introduction

Hubs are intermediate facilities serving as transshipment, sorting, and consolidation/break-bulking centers

ED

in many-to-many distribution systems. Instead of traditional point-to-point shipments, in hub networks, each origin-destination (O/D) flow is routed via one or more hub facilities resulting in a network having smaller 5

number of arcs each carrying a substantial amount of flow. As a result, transportation cost in hub networks

PT

tend to exploit the advantage of economies of scale because of high flow volumes, especially on the inter-hub connections. The hub location problem (HLP) deals with locating the hub facilities and allocating the demand

CE

nodes to the hubs in order to route the traffic between O/D pairs while optimizing a criterion of interest (cost, service level, etc.). 10

Hub networks can be classified into different types depending on the way the non-hub nodes are assigned to

AC

the hubs. In this regard, two main types of the allocations are the single allocation and the multiple allocation schemes. In single allocation networks, all the incoming and outgoing traffic to and from any non-hub node are routed through a single hub, whereas in multiple allocation networks, each non-hub node can receive and send flow through more than one hub. Figure 1 illustrates an example of multiple allocation hub network with

15

14 demand points where four nodes (shown as squares) are selected as hub facilities. ∗ Corresponding

author Email address: [email protected], [email protected] (Nader Ghaffarinasab )

Preprint submitted to Computers & Operations Research

May 1, 2018

CR IP T

ACCEPTED MANUSCRIPT

Figure 1: Example of multiple allocation hub network

AN US

From an applicability point of view, both single and multiple allocation networks are used in practice. For example, passenger airline networks typically have multiple allocation because there are flights from some non-hub cities to several or all of an airline’s hubs, whereas less-than-truckload (LTL) trucking networks may have each non-hub node (i.e., end-of-line terminal) assigned to a single break-bulk terminal (i.e., hub). Simi20

larly, some telecommunication networks employ single allocation networks to reduce the cost to construct the

M

network, and others allow or require multiple allocation, as for example to increase reliability and/or provide backups [1].

The HLP belongs to the class of NP hard problems [2]. Furthermore, it involves great amounts of resources

25

ED

and has a major impact on the operational costs and overall efficiency of the service afterwards. Therefore, development of efficient solution methods for solving more realistic, large-scale problem instances in reasonable

PT

time is of utmost importance.

In most of the studies in the literature of the HLP, the input parameters are assumed to be known with certainty. However, in real world applications, deterministic data are seldom available when deciding on the

30

CE

configuration of the hub networks as such decisions are strategic and are made long before the time that real values of these parameters become known. On the other hand, some parameters such as demands, transportation costs, etc. are subject to variation during the service life of the networks as a result of a number of factors

AC

such as changes in population size, seasonal fluctuations, quality and cost of the provided service, economic reforms, political issues, entrance of new competitors, etc. Therefore, the future values of such parameters are poorly predictable in practice and these uncertainties must be taken into account when deciding on configu-

35

ration of an efficient hub network. In order to model the uncertainty in the input parameters, there are a number of approaches depending on the type of the available information on these parameters. If the probability distributions or a set of discrete scenarios for realization of the uncertain parameters are available, techniques such as stochastic programming can be used to optimize the expected values of the considered objective functions. In contrast, there may 3

ACCEPTED MANUSCRIPT

40

be contexts where the probability distributions (or discrete scenarios) are not known and the only available information is the characteristics of a set in which the uncertain parameters are realized (i.e., the uncertainty set). In these cases, robust optimization techniques can be used in order to obtain solutions that have good performance even under worst-case conditions. In this study, we address a multiple allocation hub location problem under polyhedral demand uncertainty.

45

To this end, three different uncertainty sets are considered, namely the hose, the hybrid, and the budget

CR IP T

uncertainty sets. In the hose model, an aggregate upper bound is imposed on the sum of the inbound and outbound traffic adjacent to each node. The hybrid model, on the other hand, generalizes the hose model by considering lower and upper bounds on individual traffic demands. The last model, called the budget model, considers a budget of uncertainty that represents the maximum number of demand parameters outgoing from 50

each node that are allowed to take a value within an interval around their nominal values. A robust optimization approach based on the minimax criterion is utilized to minimize the cost of the network under the

AN US

worst-case scenario.

The problem is formulated as a mixed-integer programming (MIP) model in each case and an efficient matheuristic solution approach, combining Tabu Search (TS) and mathematical programming methods, is 55

developed to solve the problem. Obtained results from the conducted numerical experiments show the efficiency and the effectiveness of the proposed algorithm in solving instances with up to 200 nodes from three well-known benchmarking data sets in short computational times.

M

The remainder of this paper is organized as follows. The next section discusses the relevant literature for the problem at hand. In Section 3, we will present MIP model formulations for the problem under different uncertainty models. The proposed matheuristic solution approach is presented in Section 4. Computational

ED

60

experiments and corresponding results are presented in Section 5. Finally, Section 6 provides conclusions and

PT

some outlooks for future research.

2. Literature review

65

CE

Since 1980s, the HLP has been established as one of the important research areas in the field of facility location theory. The number of published studies in the HLP has grown with an increasing trend in recent

AC

years and different variants of the HLP have been studied so far. The interested readers may refer to [2, 1, 3, 4] as recent surveys on the HLP. The p-hub median problem is an important type of the HLP where a number of p hubs have to be located in

a network aiming at minimizing the total transportation cost (time, distance, etc.) for serving a set of O/D

70

demands. The network connecting the hub facilities is assumed to be a complete graph and the cost of routing on the inter-hub connections is discounted to reflect the economies of scale as a result of more consolidated flows or as a result of using more efficient means of transportation on such connections. The uncapacitated multiple allocation p-hub median problem (UMApHMP) was first formulated by Campbell [5] as a linear MIP model. Later, other formulations for this problem were proposed in [6], [7], and [8]. 4

ACCEPTED MANUSCRIPT

75

Although different variants of the HLP have been widely studied in the literature, the number of works addressing the uncertainty is rather limited. The first work on the uncertain hub location problem is done by Marianov and Serra [9] for location of airline hubs. Hubs are modeled as M/D/c queuing systems. A formula is derived for the probability of the number of customers in the system, which is later used to propose a capacity constraint. This constraint limits the probability of more than a certain number of airplanes in queue

80

to be smaller than or equal to a given value. Yang [10] addresses a problem also in the context of airline

CR IP T

transportation. Stochasticity is associated with the demands and the problem is formulated as a two stage stochastic program. The deterministic equivalent problem is solved where only three scenarios are considered. Direct shipment between non-hub nodes is also allowed. Sim et al. [11] consider a problem with normally distributed stochastic travel times. The probability of the total travel time along a path being not exceeding 85

a given time bound is controlled via chance constraints. Three scenarios with different coefficient of variation are analyzed and a heuristic algorithm is proposed.

AN US

Camargo et al. [12] propose an efficient procedure that simultaneously generates outer-approximation and Benders cuts to tackle the single allocation hub location problem under congestion modeled as a mixed-integer nonlinear program (MINLP). The proposed method is able to optimally solve large instances (up to 200 nodes) 90

in reasonable time. Contreras et al. [13] study the stochastic uncapacitated hub location problems where uncertainty is associated with demands and transportation costs. They show that the stochastic problems with uncertain demands or dependent transportation costs are equivalent to their associated deterministic expected

M

value problem (EVP) where random variables are replaced by their expectations. For the case of uncertain independent transportation costs, they propose a Monte Carlo simulation based algorithm that combines a sample average approximation scheme with a Benders decomposition algorithm to solve problems with stochas-

ED

95

tic independent transportation costs. Zhai et al. [14] consider a new class of two-stage stochastic HLPs with a minimum-risk criterion where uncertain demands are characterized by a random vector.

PT

Alumur et al. [15] address several aspects concerning HLPs under uncertainty. They consider two sources of uncertainty: the set-up costs for the hubs and the flow values to be transported between the nodes. Generic models are presented for single and multiple allocation versions of the problems. Firstly, the two sources of

CE

100

uncertainty are analyzed separately and afterwards a more comprehensive model is proposed, considering all sources of uncertainty. Chaharsooghi et al. [16] consider a reliable multiple allocation HLP where the hubs

AC

are subject to failure during their service life. The problem is formulated as a two-stage stochastic model and a metaheuristic algorithm is proposed for solving it.

105

Shahabi and Unnikrishnan [17] address robust HLPs under single and multiple allocation settings with uncertain demands. The problems are formulated as MINLPs and then transformed into mixed integer conic quadratic programs and an efficient linear relaxation strategy is proposed for solving them. Ghaffari-Nasab et al. [18] study a capacitated robust HLP for both the single and multiple allocation networks under demand uncertainty. The uncertainty is considered in the capacity constraints and the objective function is calculated

110

based of the average (nominal) demand values. Habibzadeh Boukani et al. [19] consider the capacitated

5

ACCEPTED MANUSCRIPT

single and multiple allocation HLPs with uncertain hub setup costs and capacities. The problem is modeled using a minimax regret criterion and test instances from Iranian Aviation Dataset are solved using a standard optimization package. Ghaderi and Rahmaniani [20] address a single allocation HLP under demand and travel time uncertainties. The problem is modeled based on an objective of minimizing the total expected 115

transportation costs, while bounding the relative regret in each scenario and hybrid metaheuristic algorithms are proposed for solving it.

CR IP T

Meraklı and Yaman [21] study a robust multiple allocation HLP under polyhedral demand uncertainty. They model the demand uncertainty using hose and hybrid models. The problem is formulated as linear MIP models using a minimax criterion and two Benders decomposition algorithms are proposed in order to solve 120

the problem on instances with up to 150 nodes. Zetina et al. [22] present robust counterparts for a multiple allocation HLP where the level of conservatism is adjusted by an uncertainty budget. The problem is studied under three different settings: uncertain demands, uncertain transportation costs, and uncertain demands and

AN US

transportation costs. Using a branch-and-cut algorithm implemented on a commercial solver, the authors are able to solve instances with up to 50 nodes from these problems. 125

In another study, Meraklı and Yaman [23] address a capacitated multiple allocation hub location problem under hose demand uncertainty. They consider fixed costs and capacity constraints on hubs and solve instances with up to 50 nodes using two different Benders decomposition algorithms. More recently, Martins de S´ a et al. [24] apply Benders decomposition on the multiple allocation incomplete hub location problem where the

130

M

O/D demands and hub fixed costs are subject to uncertainty. They address the data uncertainty by adopting a robust optimization approach and solve problem instances with up to 100 nodes using the proposed solution

ED

algorithms.

Although integer programming optimization approaches are used to solve various types of the HLP in small sizes, larger instances are usually solved by heuristic, metaheuristic, or matheuristic procedures. In fact, de-

135

PT

velopment of efficient heuristic algorithms has helped researchers to obtain optimal or near-optimal solutions for problems of large sizes in short computational times.

CE

In case of the single allocation HLP, the number of proposed heuristic algorithms is large. O’Kelly [25] proposes two heuristic allocation procedures for solving the uncapacitated single allocation p-hub median problem (USApHMP). A tabu search (TS) heuristic is proposed for the USApHMP in [26]. Ernst and Krishnamoorthy

140

AC

[27] develop a simulated annealing (SA) heuristic for the same problem and show that it is comparable, in both solution quality and computational time, to the TS heuristic in [26]. Ernst and Krishnamoorthy [28] propose heuristic algorithms for solving the capacitated single allocation HLP based on SA and random descent heuristic. Another SA heuristic is proposed for the USApHMP in [29]. Topcuoglu et al. [30] propose a genetic algorithm (GA) for the uncapacitated single allocation HLP. Chen [31] presents a hybrid heuristic based on SA embedded with a tabu list and some improvement procedures for this problem. Silva and Cunha [32] propose

145

three variants of a simple and efficient multi-start TS heuristic as well as a two-stage integrated TS heuristic to solve the USAHLP. Calık et al. [33] develop a TS heuristic for the single allocation hub covering problem

6

ACCEPTED MANUSCRIPT

over incomplete hub networks. Jabalameli et al. [34] propose an SA based heuristic for solving the single allocation p-hub maximal covering problem. Abyazi-Sani and Ghanbari [35] develop an efficient TS based heuristic for solving the uncapacitated single allocation HLP. Yang et al. [36] present an SA based heuristic 150

for solving the single allocation HLP under mixed uncertainty. More recently, Silva and Cunha [37] propose another TS heuristic for the uncapacitated single allocation p-hub maximal covering problem. Some authors have also tackled the multiple allocation HLPs using heuristic algorithms. Campbell [38] pro-

CR IP T

pose a greedy-interchange heuristic to solve the UMApHMP. Ernst and Krishnamoorthy [8] present an LP based branch-and-bound method to obtain exact solutions to the same problem. Marianov et al. [39] propose 155

a three-phase TS heuristic for solving the competitive HLP in a multiple allocation setting. Marianov and Serra. [9] propose another TS heuristic for solving the HLP under congestion. Boland et al. [40] develop preprocessing techniques and fastening constraints for the UMApHMP by exploiting some characteristics of optimal solutions. A GA based solution algorithm is proposed in [41] for solving the multiple allocation com-

160

AN US

petitive hub location and pricing problem in which the inter-hub network was assumed incomplete. Sender et al. [42] present two matheuristics based on a local search and an evolutionary algorithm for the capacitated multiple allocation HLP derived from a practical application in network design of German wagonload traffic. Ghaffarinasab et al. [43] tackle the competitive multiple and single allocation HLPs as Stackelberg games by using SA algorithms. In another study, Ghaffarinasab and Motallebzadeh [44] propose SA based metaheuristic algorithms for solving three variants of the hub interdiction problem under the multiple allocation setting. In this study, we take a more comprehensive and unified approach to model the robust multiple allocation p-

M

165

hub median problem under polyhedral demand uncertainty by considering three different types of uncertainty

ED

sets. Since, to the best of our knowledge, no heuristic solution procedure for solving large-scale instances of these problems has been proposed in the literature, we develop an efficient TS based matheuristic algorithm for solving these problems in short computational times and also study the effect of different input parameters

PT

on characteristics of the final solution by doing an extensive set of computational experiments.

CE

3. Mathematical formulation

Let G = (N, A) be a network, where N is the set of nodes and A is the set of arcs such that A = {(i, j) : i, j ∈ N, i 6= j}. Assume H ⊆ N be a subset of nodes that is available for locating hubs. The number of hubs to be located is fixed and is denoted by p. For all (i, j) ∈ A, let wij denote the amount of flow originated at

AC

170

node i and destined to node j (we assume that wii = 0 for all i ∈ N ), and dij denote the transportation cost

of a unit flow from node i to node j. The transportation cost associated with routing one unit of flow wij via hubs k and m is calculated as: ckm ij = χdik + αdkm + δdmj where χ, α, and δ are cost coefficients associated respectively with collection, transfer, and distribution links on each route. More specifically, α is the volume discount factor to reflect a scale economies on transportation cost via the inter-hub connections (0 ≤ α ≤ 1; χ ≥ α; δ ≥ α). 7

ACCEPTED MANUSCRIPT

Let the variable xkm ij denote the fraction of flow wij that is sent from node i to node j using the link between the hubs k and m. Let also the binary variable yk ∈ {0, 1} be 1 if node k is selected as a hub and 0, otherwise. The problem consists of the selection of p nodes which will act as hubs and determining how the O/D flows will be assigned to the hub nodes so that total transportation cost is minimized. The MIP model for the deterministic version of the UMApHMP proposed by Hamacher et al. [45] can be written as: min

X X X

km wij ckm ij xij

(1)

s.t.:

X

CR IP T

(i,j)∈A k∈H m∈H

yk = p

k∈H

X X

xkm ij = 1

k∈H m∈H

X

X

xkm ij +

m∈H

m∈H|m6=k

xmk ij ≤ yk

∀(i, j) ∈ A

(3)

∀(i, j) ∈ A, k ∈ H

(4)

∀(i, j) ∈ A, k, m ∈ H

(5)

∀k ∈ H

(6)

AN US

xkm ij ≥ 0 yk ∈ {0, 1}

(2)

The objective function (1) minimizes the total transportation cost. 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

designated as hubs. (5) and (6) are the standard domain constraints for the decision variables.

ED

3.1. Case I: Hose demand uncertainty

Under the hose demand uncertainty, it is assumed that the total flow adjacent to each node i ∈ N does defined as [21]:

PT

not exceed a non-negative upper bound value bi . In other words, the hose uncertainty set can be formally

CE

Dhose =

(

w∈

n(n−1) R+

:

X

wij +

j∈N \{i}

X

j∈N \{i}

wji ≤ bi , ∀i ∈ N

)

In the robust multiple allocation p-hub median problem under hose uncertainty, we aim at deciding on the

AC

175

M

is routed via some hub pair. Constraints (4) state that the flows can only be routed via nodes that have been

locations of hubs and the allocation of the O/D pairs in such a manner that the worst case cost over all possible demand realizations within the set Dhose is minimized: min

max

(x,y)∈Ω w∈Dhose

X X X

km wij ckm ij xij

(i,j)∈A k∈H m∈H

where Ω is the set of all solutions (x, y) to the upper level (outer) problem that satisfy constraints (2)-(6). The above model is a nonlinear model that can be linearized using the dual transformation method as in [21]. Based on this linearization scheme, which is used for minimax type robust optimization problems, the dual of 8

ACCEPTED MANUSCRIPT

the lower level (inner) problem is derived and then the problem is converted to a pure (linear) minimization problem. To this end, let (ˆ x, yˆ) ∈ Ω be a given solution for the upper level problem. Then, by fixing the (ˆ x, yˆ) in the lower level problem, we have the following linear programming (LP) problem: X X X

max

w∈Dhose

wij ckm ˆkm ij x ij

(i,j)∈A k∈H m∈H

The above LP is feasible and bounded and therefore, its optimal objective value is the same as that of its dual

CR IP T

problem. Based on this result, we can formulate the robust UMApHMP under the hose demand uncertainty as the following (linear) MIP: min

X

λ i bi

i∈N

λi + λj ≥

X X

km ckm ij xij

k∈H m∈H

λi ≥ 0

∀(i, j) ∈ A

AN US

s.t.: (2) − (6)

∀i ∈ N

where λi is the dual variable associated with the constraint 3.2. Case II: Hybrid demand uncertainty

P

j∈N \{i}

wij +

(7) (8) (9) (10)

P

j∈N \{i}

wji ≤ bi for i ∈ N.

M

The hybrid demand uncertainty generalizes the hose model by restricting every individual demand value wij to an interval of the form [lij , uij ] for all (i, j) ∈ A. In other words, we can formally define the hybrid (

ED

uncertainty set as [21]:

w∈

: lij ≤ wij ≤ uij , ∀(i, j) ∈ A

)

PT

Dhybrid = Dhose ∩

n(n−1) R+

problem.

CE

It should be noted that if we take lij = 0 and uij ≥ min{bi , bj } for all (i, j) ∈ A, then we have Dhybrid = Dhose . P Furthermore, if we take lij = uij and bi ≥ j∈N \{i} (uij +uji ), the robust problem reduces to the deterministic

Based on our earlier discussions over the hose model, the linear MIP model for the robust UMApHMP under

AC

the hybrid demand uncertainty can be derived in a similar manner as follows: min

X

λi bi +

i∈N

X

(uij βij − lij µij )

(11)

(i,j)∈A

s.t.: (2) − (6) λi + λj + βij − µij ≥

X X

km ckm ij xij

∀(i, j) ∈ A

(12)

λi ≥ 0

∀i ∈ N

(13)

βij , µij ≥ 0

∀(i, j) ∈ A

(14)

k∈H m∈H

9

ACCEPTED MANUSCRIPT

where for all (i, j) ∈ A, the variables βij and µij are the dual variables associated with the constraints wij ≤ uij and wij ≥ lij , respectively. 3.3. Case III: Budget demand uncertainty In our last model for polyhedral demand uncertainty, it is assumed that for every (i, j) ∈ A, the flow value wij is realized within an interval of the form [lij , uij ]. To this end, let w ¯ij and rij respectively denote the

CR IP T

center (the nominal value) and the radius of corresponding interval for wij . In other words, we assume that wij ∈ [w ¯ij − rij , w ¯ij + rij ] for every (i, j) ∈ A, where w ¯ij = (lij + uij )/2 and rij = uij − w ¯ij .

For every i ∈ N , we introduce the parameter Γi denoting the maximum number of demand parameters originated from node i whose values are allowed to deviate from the corresponding nominal values. In other words, we define the budget uncertainty set as: ( w∈

n(n−1) R+

: lij ≤ wij ≤ uij , ∀(i, j) ∈ A,

X

j∈N \{i}

|wij − w ¯ij | ≤ Γi , ∀i ∈ N rij

AN US

Dbudget =

)

Note that Γi can take values within the interval [0, |N | − 1] and its role is to adjust the robustness of the proposed method against the level of conservatism of the solution. As a special case, when Γi = 0 for all i ∈ N , we have the deterministic problem. In addition, the case with Γi = |N | − 1 corresponds to the worst

M

case problem where all the outgoing flow values from each node take their worst values from the corresponding intervals of uncertainty.

ED

Similar to the previous cases, our aim is to minimize the worst case cost over all possible demand realizations within the set Dbudget . Therefore, the robust HLP under budget uncertainty can be written as the following minimax model:

X X X

km w ¯ij ckm ij xij

PT

min

(i,j)∈A k∈H m∈H

+

Xn

CE

i∈N

max

Si ⊂N \{i}:|Si |=Γi

X

j∈Si \{i}

rij

X X

k∈H m∈H

km ckm ij xij

o

(15)

s.t.: (2) − (6)

Let the variable 0 ≤ zij ≤ 1, ∀(i, j) ∈ A, denote the absolute deviation from the nominal value of the flow

AC

180

variable wij (in units of rij ). By fixing the (ˆ x, yˆ) in the lower level problem for each i ∈ N , we have the

following LP problem: max

X

j∈N \{i}

s.t.:

X

j∈N \{i}

rij

X X

k∈H m∈H

 ckm ˆkm zij ij x ij

zij ≤ Γi

0 ≤ zij ≤ 1 10

(16) ∀i ∈ N

(17)

∀(i, j) ∈ A

(18)

ACCEPTED MANUSCRIPT

Taking the dual of the above LP problem, the robust HLP under the budget demand uncertainty can be formulated as the following linear MIP model: X X X

km w ¯ij ckm ij xij +

(i,j)∈A k∈H m∈H

s.t.: (2) − (6) θi + ϕij ≥ rij

X

Γi θi +

i∈N

X X

km ckm ij xij

k∈H m∈H

X

ϕij

(19)

(i,j)∈A



∀(i, j) ∈ A

CR IP T

min

θi ≥ 0 ϕij ≥ 0

(20)

∀i ∈ N

(21)

∀(i, j) ∈ A

(22)

where the variables θi and ϕij are the dual variables associated with the constraints (17) and (18), respectively.

AN US

4. Matheuristic solution algorithm

Matheuristics are hybrid algorithms combining mathematical programming techniques with metaheuristic 185

solution approaches [46]. These algorithms are getting more and more popular since they have successfully been implemented to many combinatorial optimization problems (see [47, 48, 49, 50], as some examples). Within the field of facility location and network design, some papers have been published on successful application of

programming techniques is given in [46]. 190

M

matheuristics, e.g., [51, 52, 53, 54]. A general overview on combinations of metaheuristics with mathematical

In our proposed matheuristic algorithm, TS is used to solve the (higher level) master problem which is

ED

comprising the location decisions, whereas a standard commercial solver is used to solve the (lower level) subproblems which are linear programming problems for determining the worst-case O/D flows through the

PT

network nodes.

The tabu search (TS) algorithm proposed by Glover [55], is a local search algorithm which has effectively 195

tackled variety of hard optimization problems. This procedure starts with an initial solution and uses a tabu

CE

list to control the moves in the neighborhood so that trapping in local optima and re-visiting the same solution will not occur. From the current solution, all the non-tabu moves are explored and the best one is selected. This move which might lead to an either better or worse solution than the current solution, is recorded in the

200

AC

tabu list. The future move is among those not listed in the tabu list, unless it fulfills an aspiration criterion, i.e., it has a better objective value than the best found objective. TS would be terminated whenever a termination criterion such as a maximum number of iterations, a fixed number of iterations with no improvement in the best solution, etc. is met. In our implementation, a one-dimensional array is used to represent the solutions. This array of size p includes the numbers associated with the nodes that are selected as hubs. Note that having known the selected hubs

205

(the yk variables), for each O/D pair (i, j), one can easily determine the paths for routing the associated flow wij (the xkm ij variables) by solving a shortest path problem via enumeration. Afterwards, as stated above, the 11

ACCEPTED MANUSCRIPT

worst-case realization of the flow variables (wij ) within the corresponding uncertainty set will be determined by solving the lower level problem using a standard optimization package yielding the value of the objective function for the current solution. 210

The proposed TS algorithm uses two input parameters, namely Ltabu and Imax . Ltabu represents the size of the tabu list and Imax denotes the number of iterations with no improvement in the best solution. In other words, the algorithm would be terminated when it cannot improve the best solution in the last Imax iterations.

CR IP T

To produce an initial solution to our algorithm, we randomly select p nodes from the set H as hub nodes. A simple neighborhood structure, called the “Swap Hub” operator, is employed to explore the solution space. 215

Based on this operator, for a given pair of a hub node and a non-hub node, the two nodes are swapped. In other words, the hub node becomes non-hub and the non-hub node becomes hub. Having altered the set of open hubs using the above mentioned operator, the allocations of flows are determined by solving shortest path problems via enumeration and the flow volumes as well as the corresponding value of the objective function

220

AN US

are determined by solving the lower level problem using the standard solver.

To elaborate more precisely, we start our algorithm by generating an initial solution Xp consisting of the locations of the selected p hubs. The tabu list is initially empty (i.e., TL = {}). We set the iteration counter I to 0 and set the Xp and its objective function f (Xp ) as the best solution Xp∗ and best objective value f ∗

found so far, respectively (i.e., Xp∗ ← Xp , f ∗ ← f (Xp )). At each iteration, we select a random hub node r ∈ Xp and replace it with all the non-hub nodes s ∈ H \ Xp one at a time. Specifically, for each s ∈ H \ Xp

we generate a new solution Xp00 based on the current solution (Xp ) by using the “Swap Hub” operator as

M

225

Xp00 ← Xp \ {r} ∪ {s}. Subsequently, the flow assignments are determined based on the new set of hubs using

ED

enumeration and the corresponding objective function, f (Xp00 ), is calculated. Then the best non-hub node

which has resulted in the least objective function after applying the “Swap Hub” is determined and labeled / TL, the tabu list is as s∗ and the corresponding solution is labeled as Xp0 (i.e., Xp0 ← Xp \ {r} ∪ {s∗ }). If s∗ ∈ updated based on the first-in-first-out (FIFO) rule. In other words, the node that has just been removed from

PT

230

the set of hubs (i.e., the node r) is not allowed to become hub in the next Ltabu iterations. We update the

CE

current solution as Xp ← Xp0 . In addition, if the objective function f (Xp0 ) of the new solution Xp0 is smaller

than the best objective so far, f ∗ , we also set Xp∗ ← Xp0 , f ∗ ← f (Xp0 ), and I ← 0. Otherwise, if s∗ ∈ TL,

we accept the new solution only if it satisfies the aspiration criterion, i.e., f (Xp0 ) < f ∗ . In this case, we set Xp∗ ← Xp0 and f ∗ ← f (Xp0 ). We also update the current solution as Xp ← Xp0 . If the best known solution

AC

235

is not improved in the last iteration, we set I ← I + 1. The algorithm is terminated when the number of iteration with no improvement in the best solution (I) exceeds the Imax . The pseudo-codes of the proposed matheuristic algorithm is illustrated in Algorithm 1.

5. Computational experiments 240

In order to test the efficiency of the proposed matheuristic algorithms, we use three well-known data sets from the literature of the HLP, namely the CAB, the TR, and the AP data sets. The CAB data set introduced 12

ACCEPTED MANUSCRIPT

Algorithm 1 : The proposed TS based matheuristic Generate a random initial set of hubs Xp ;

2:

Initialize the tabu list TL={};

3:

Determine the flow assignments based on Xp using enumeration;

4:

Calculate f (Xp ) by solving the lower level problem using standard solver;

5:

I ← 0, f ∗ ← f (Xp ), Xp∗ ← Xp ;

6:

while I < Imax do

7:

I ← I + 1;

8:

ftest ← ∞;

9:

Select a random hub node r ∈ Xp ;

10: 11:

CR IP T

1:

for each s ∈ H \ Xp do

Generate a new solution Xp00 based on Xp as Xp00 ← Xp \ {r} ∪ {s}; Determine the flow assignments based on Xp00 using enumeration;

13:

Calculate f (Xp00 ) by solving the lower level problem using standard solver;

14:

if f (Xp00 ) < ftest then

15:

ftest ← f (Xp00 );

s∗ ← s;

16:

end if end for

19:

Xp0 ← Xp \ {r} ∪ {s∗ };

20:

if s∗ ∈ / TL then

ED

18:

21:

Update TL;

22:

if f (Xp0 ) < f ∗ then

Xp∗ ← Xp0 , f ∗ ← f (Xp0 );

24:

I ← 0;

27:

Xp ← Xp0 ;

else

if f (Xp0 ) < f ∗ then

AC

28:

CE

26:

end if

PT

23:

25:

M

17:

AN US

12:

29:

Xp∗ ← Xp0 , f ∗ ← f (Xp0 );

30:

I ← 0;

Xp ← Xp0 ;

31:

32: 33:

end if end if

34:

end while

35:

return Xp∗ , f ∗ .

13

ACCEPTED MANUSCRIPT

by O’Kelly [25] is based on the airline passenger interactions between 25 US cities in 1970 evaluated by the Civil Aeronautics Board (CAB). As in most of the published works, the original flow values in the CAB data set are scaled so that their sum equals to 1. The second data set used in our computational experiments is 245

the TR data set [56] which is based on the cargo flows between 81 cities of Turkey. The flow values in the TR data set are scaled in such a way that their sum equals 1000. For the CAB and TR data sets, the values of the parameters χ and δ are set to 1, whereas the parameter α is considered at four levels as: α ∈ {0.2, 0.4, 0.6, 0.8}.

CR IP T

The third data set is the Australia Post (AP) data set introduced by Ernst and Krishnamoorthy [27]. The AP data set is based on a postal delivery system in Sydney, Australia and consists of 200 nodes representing 250

postal districts. In our experiments, we use instances with |N | ∈ {25, 40, 50, 75, 100, 125, 150, 200} from the AP data set. The values of the parameters χ, α, and δ in the AP data set are selected as 3, 0.75, and 2, respectively. Contrary to the case of CAB and TR data sets, the AP data set includes demands originating and ending at the same node, i.e., wii 6= 0, ∀i ∈ N . However, for the uniformity of computation, we set wii = 0

255

AN US

for all i ∈ N in the AP data set. For all the three data sets, it is assumed that all the demand nodes are candidate locations for installing hub facilities (i.e., H = N ).

As in [21], the value of the upper bounds adjacent to each node is defined as bi =

P

j∈N \{i} (wij

+ wji )

for all i ∈ N . For the hybrid model, on the other hand, the lower and upper bounds on the individual flow values are respectively set as lij = max{0, (1 − ψ)wij } and uij = (1 + ψ)wij for all (i, j) ∈ A, with ψ ∈ {0.2, 0.4, 0.6, 0.8, 1, 2}. For the budget uncertainty model, the radius of the uncertainty intervals is selected as rij = ρw ¯ij with ρ ∈ {0.2, 0.5}. Furthermore, the value of the budget of uncertainty (Γi ) is selected

M

260

from among the set {0, 1, 2, 5, 10, 15, 20} for different problem instances.

ED

The proposed matheuristic algorithm is implemented in JAVA and the mathematical models for the subproblems are solved by CPLEX version 12.6 using the dual simplex algorithm. We set the time limit of three hours. All the experiments have been run on a computer with Intel(R) Core(TM) i3-3220 CPU of 3.30 GHz

PT

and 16 GB of RAM, using the Microsoft Windows 7 operating system. In our computational experiments, each problem instance is solved by the TS algorithm for five times and the best solutions obtained are reported.

CE

5.1. Parameters tuning

Parameter tuning is important for proper functioning of metaheuristic (or matheuristic) algorithms as it has a great influence on the efficiency and effectiveness of these algorithms. As mentioned earlier, the

AC

265

proposed TS based matheuristic has two parameters that need to be tuned, namely the Ltabu and Imax . In order to obtain a good value for these parameters, we apply the full factorial method, widely used for optimizing process parameters in design of experiments (DOE) [57]. We use three different levels for each of these parameters as Ltabu ∈ {4, 7, 10} and Imax ∈ {5, 10, 15}. Therefore, we have nine different combinations of the design parameters that need to be considered in the experiments. Further, in order to study the effect of the type and size of the data sets on the optimal values of the design parameters, we conduct the full factorial experiments separately on the CAB, TR, and AP data sets (respectively as small, medium, and large data sets). To this end, we select 30 instances of the problem with hybrid demand uncertainty 14

ACCEPTED MANUSCRIPT

from each data set for which the value of optimal solution is already found by the Benders algorithm in [21]. For the CAB and TR data sets, the 30 instances are obtained by combining the problem parameters as p ∈ {4, 5}, α ∈ {0.4, 0.6, 0.8}, and ψ ∈ {0.2, 0.4, 0.6, 0.8, 1}, whereas for the AP data set, the 30 instances are selected as |N | ∈ {125, 150}, p ∈ {3, 4, 5}, and ψ ∈ {0.2, 0.4, 0.6, 0.8, 1}. The selected instances are then solved by the proposed TS under all nine different combinations of the design parameters (Ltabu and Imax ). In order to be able to compare the objective values of different instances, we use the percentage gap (%GAP) criterion

%GAP =

CR IP T

in our experimentations as follows:

OFT S − OFOP T × 100 OFOP T

(23)

where OFOP T and OFT S represent the optimal objective value (obtained by the Benders algorithm [21]) and the objective value found by the proposed TS, respectively. The designed full factorial experiments are 270

conducted by means of analysis of variance (ANOVA) in R programming language [58] at significance level of

AN US

0.05. As a followup to ANOVA test, we perform the Fisher’s least significant difference (LSD) test [57] for pairwise comparisons of the parameter level in case they are shown to be significantly different by ANOVA. The results obtained by the conducted experiments are summarized in Figures 2 - 4 for the CAB, TR, and AP instances, respectively. These figures show the means plots as well as the corresponding 95% significance 275

intervals for the %GAP values.

Based on ANOVA results, the mean %GAPs under the three levels of the parameter Ltabu do not differ

M

significantly for all the three data sets. Hence, the three vertical lines representing the confidence intervals in the plots for Ltabu are shown in blue color in Figures 2 - 4. However, the difference between the mean %GAPs

level value (the value 5) results in solutions that are significantly different from those of the two other levels (the values 10 and 15). However, the solutions obtained by the two levels 10 and 15 do not differ significantly.

2.0

●

%GAP

●

1.5

●

1.4 1.2

%GAP

1.6

2.5

1.8

CE

PT

For this reason, the vertical lines corresponding to the values of 10 and 15 are shown in red color.

1.0

1.0

●

●

0.5

●

0.8

AC

280

ED

under the three levels of the parameter Imax are shown to be significant. The LSD test indicates that the first

4

7

10

5

Ltabu

10

15

Imax

Figure 2: The means and 95% confidence interval plots for %GAP with the CAB data set

15

1.5

2.5

2.0

3.0

ACCEPTED MANUSCRIPT

●

1.0

1.0

●

1.5

%GAP

%GAP

2.0

●

CR IP T

0.5

0.5

●

●

●

4

7

10

5

Ltabu

10

15

Imax

2.5

1.4

3.0

AN US

1.6

Figure 3: The means and 95% confidence interval plots for %GAP with the TR data set

%GAP

1.0

●

0.8 0.4

7

● ●

10

5

ED

4

0.5

M

0.6

1.0

●

0.0

%GAP

●

1.5

2.0

1.2

●

Ltabu

10

15

Imax

PT

Figure 4: The means and 95% confidence interval plots for %GAP with the AP data set

Based on the above results, we select the value of the parameter Ltabu as 7 for the CAB and TR data set

set to 10 for the TR and AP data sets. The reason is that the results under the two values 10 and 15 are not significantly different but setting the value to 10 results in smaller computational times. Table 1 shows the selected values for the algorithm parameters for different data sets.

AC

285

CE

and as 10 for the AP data set. The parameter Imax is set as 15 for the CAB data set. However, its value is

Table 1: Selected parameter values for different test data

data set

Ltabu

Imax

CAB

7

15

TR

7

10

AP

10

10

16

ACCEPTED MANUSCRIPT

5.2. Results for the hose and hybrid uncertainty cases Table 2 shows the results obtained by solving the problem under the hose and hybrid uncertainty models 290

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 column shows the objective function (OF) value of best solution obtained by the TS (in five runs of the algorithm). Since the problem instances are solved to optimality using Benders decomposition algorithm in [21], we compare the quality of our solutions to the corresponding optimal

295

CR IP T

solutions. To this end, we have reported the percentage gap (%GAP) between the best objective value found by the TS and the optimal objective value. These values are shown below the objective values in the tables. The average CPU time (in seconds) for the five runs of the TS algorithm as well as the number of successful runs (SR) in which the TS has found its best solution are reported in the next column. Finally, the last row (CP U (s)) shows the average computational times for solving instances using the proposed TS algorithm. The results reported in Table 2 indicate that the proposed TS algorithm is able to obtain the optimal solution for all the instances of the CAB data set in very short computational times. All the optimal solutions have

AN US

300

been obtained in less than half of a second. The small solution times as well as large values of successful runs (SR) show the efficiency of the proposed solution algorithm. Note that as the value of the parameter ψ increases, the optimal objective function value for the hybrid model rises. Furthermore, the optimal objective values of the hybrid model are in general smaller than that of the hose model. However, as the value of ψ increases, the difference between the objectives of the two models diminish. This is due to the fact that the

M

305

uncertainty set for the hose model is larger than the corresponding set for the hybrid model. Nevertheless, as the value of parameter ψ gets larger, the two uncertainty sets get closer to each other. For this reason, the

ED

optimal objective values for the hybrid model with ψ = 2 are very close to those of the corresponding instances under the hose model. Observe also that for a given number of the opened hubs p, the optimal objective value increases as the value of the discount factor (α) increases. In contrast, for a fixed value of the discount factor,

CE

PT

the value of the total transportation cost decreases as the number of opened hubs (p) increases.

AC

310

17

ACCEPTED MANUSCRIPT

Table 2: Results for the problem under hose and hybrid uncertainty with the CAB data set Hybrid (ψ = 0.2)

5

Hybrid (ψ = 0.8)

Hybrid (ψ = 1)

Hybrid (ψ = 2)

Hose

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

0.2

1007.72

0.38

1019.41

0.27

1031.10

0.25

1042.80

0.24

1054.49

0.25

1054.99

0.19

1054.99

0.19

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

0.4

1095.61

0.32

1118.73

0.30

1141.84

0.29

1164.96

0.32

1188.08

0.30

1190.79

0.22

1190.79

0.14

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

0.6

1172.03

0.3

1206.98

0.37

1241.87

0.36

1269.64

0.33

1297.42

0.32

1319.78

0.25

1319.78

0.15

0.00

5

0.00

5

0.00

5

0.00

5

0.00

0.8

1222.71

0.38

1256.55

0.38

1290.39

0.36

1318.08

0.41

1342.32

0.00

4

0.00

5

0.00

5

0.00

5

0.00

0.2

770.59

0.34

788.02

0.34

805.36

0.43

822.70

0.37

839.44

0.00

5

0.00

5

0.00

5

0.00

5

0.00

0.4

893.41

0.39

927.19

0.37

960.96

0.37

994.66

0.4

1024.40

0.00

5

0.00

5

0.00

5

0.00

5

0.00

0.6

996.94

0.37

1044.22

0.36

1091.50

0.38

1136.48

0.33

1180.64

0.00

5

0.00

5

0.00

5

0.8

1079.13

0.35

1136.03

0.35

1190.64

0.32

0.00

5

0.00

5

0.00

5

0.2

635.69

0.31

652.91

0.36

670.12

0.37

0.00

5

0.00

5

0.00

5

0.4

788.62

0.31

821.96

0.35

854.22

0.39

0.00

5

0.00

5

0.00

4

0.6

914.26

0.40

962.07

0.30

1009.88

0.30

0.00

2

0.00

5

0.00

0.8

1013.03

0.31

1074.31

0.34

1135.59

0.00

5

0.00

2

0.00

3

0.00

1

0.00

2

0.00

5

0.00

5

0.2

547.75

0.35

565.50

0.29

583.25

0.27

601.00

0.31

618.74

0.28

639.79

0.26

646.72

0.18

0.00

3

0.00

0.4

711.42

0.39

746.51

0.00

1

0.00

0.6

855.24

0.28

905.78

0.00

1

0.8

974.35

0.33

1037.38

5

0.00

0.00

0.34

0.00

5

0.00

1244.66

0.36

1293.22

CPU(s)

5

0.00

5

0.00

5

0.38

1417.49

0.26

1418.84

0.18

5

0.00

5

0.00

5

0.36

845.12

0.32

845.26

0.16

5

0.00

5

0.00

5

0.42

1036.58

0.33

1037.64

0.16

5

0.00

5

0.00

5

0.34

1209.00

0.32

1213.09

0.15

5

0.00

5

0.00

5

0.35

1359.06

0.33

1367.93

0.18

0.00

5

0.00

5

0.00

5

0.00

5

687.33

0.32

704.54

0.33

722.29

0.28

726.44

0.15

0.00

5

0.00

5

0.00

5

0.00

5

886.47

0.30

918.73

0.34

954.92

0.29

967.16

0.15

0.00

5

0.00

4

0.00

3

0.00

5

1057.70

0.33

1105.51

0.29

1156.82

0.27

1170.07

0.16

0.00

1196.86

2

0.00

4

0.00

2

0.00

5

0.39

1251.39

0.31

1326.78

0.27

1343.21

0.16

2

0.00

5

0.00

3

0.00

1

0.00

4

0.00

5

0.29

781.60

0.29

816.68

0.31

851.77

0.31

899.59

0.23

914.10

0.17

2

0.00

2

0.00

4

0.00

2

0.00

2

0.00

3

0.28

956.32

0.32

1005.79

0.24

1055.19

0.24

1112.80

0.23

1129.91

0.16

2

0.00

5

0.00

4

0.00

2

0.00

3

0.00

5

0.29

1098.67

0.37

1158.20

0.28

1215.17

0.29

1298.23

0.22

1322.23

0.17

2

0.00

4

0.00

4

0.00

2

0.00

2

0.00

0.33

0.34

0.33

0.32

0.27

2 0.16

CE

CP U (s)

0.00

2 0.34

CR IP T

OF

AN US

4

Hybrid (ψ = 0.6)

CPU(s)

M

3

Hybrid (ψ = 0.4)

OF

ED

2

α

PT

p

The results for solving the hose and hybrid models with the TR data set are presented in Table 3. As

AC

can be seen from the presented results, the proposed TS algorithm has found the optimal solution for all the instances in quite short computational times. Note that the solution times for the hose model are considerably

315

shorter than the solution time for the hybrid model. This can be the attributed to the smaller time needed for solving the lower level problem of the hose model that has smaller number of variables in comparison with the lower level problem of the hybrid model. As in the case of the CAB data set, one can observe that the optimal value of the objective function for the hybrid model when the parameter ψ has larger values, is very

close (identical for some instances) to those of the hose model. Furthermore, the changes in the optimal value 320

of the objective function as a result of varying values for the input parameters p and α are similar to the

18

ACCEPTED MANUSCRIPT

corresponding changes with the CAB data set. Table 3: Results for the problem under hose and hybrid uncertainty with the TR data set Hybrid (ψ = 0.2)

5

Hybrid (ψ = 0.8)

Hybrid (ψ = 1)

Hybrid (ψ = 2)

Hose

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

0.2

781669.72

14.92

791979.72

15.16

802289.71

13.59

812599.71

13.82

822909.71

11.32

823485.27

10.92

826877.58

4.08

0.00

5

0.00

1

0.00

5

0.00

4

0.00

5

0.00

5

0.00

5

0.4

840112.66

15.99

859638.82

15.78

879164.99

15.47

892040.38

15.31

902575.19

12.94

902575.19

11.51

902575.19

3.71

0.00

5

0.00

5

0.00

5

0.00

4

0.00

0.6

883983.61

16.23

910717.95

17.37

926290.29

17.42

940622.96

15.44

954955.62

0.00

5

0.00

4

0.00

3

0.00

5

0.00

0.8

909256.41

17.16

938955.05

17.81

959486.50

20.00

978472.84

18.24

996333.67

0.00

5

0.00

5

0.00

5

0.00

5

0.00

0.2

669320.24

13.38

678208.41

14.29

687096.58

13.66

695984.75

14.38

704872.91

0.00

2

0.00

2

0.00

2

0.00

5

0.00

0.4

743571.26

15.46

760945.74

12.80

778320.22

16.29

795694.70

15.99

812263.87

0.00

5

0.00

5

0.00

3

0.00

2

0.00

0.6

802850.50

14.59

827328.13

14.45

851179.46

13.33

874652.79

15.32

896670.92

0.00

2

0.00

5

0.00

5

0.8

845601.96

14.26

876890.63

16.18

908179.30

16.05

0.00

5

0.00

1

0.00

1

0.2

580050.10

15.75

589882.64

14.74

598397.47

17.14

0.00

2

0.00

5

0.00

5

0.4

676223.28

14.24

694784.44

18.02

713345.61

18.31

0.00

2

0.00

4

0.00

5

0.6

755449.94

14.87

780891.70

16.40

804676.28

12.38

0.00

4

0.00

5

0.00

0.8

811709.80

14.4

843479.75

13.16

875182.54

0.00

5

0.00

2

0.00

0.2

501839.91

17.86

511185.49

18.84

520391.93

0.00

5

0.00

0.4

613491.23

14.21

631820.52

0.00

2

0.00

0.6

705452.52

11.72

732038.97

0.00

5

0.00

0.8

779668.60

12.83

812942.30

1

0.00

0.00

12.41

14.87

0.00

650149.82

0.00

5

0.00

4

954955.62

12.97

954955.62

3.81

CR IP T

5

5

0.00

5

0.00

5

16.93

996504.70

14.62

996504.70

4.13

5

0.00

5

0.00

5

13.22

704872.91

12.78

704872.91

4.01

5

0.00

5

0.00

5

14.15

812263.87

14.08

812263.87

5.09

5

0.00

5

0.00

5

14.8

896670.92

12.97

896670.92

4.08

0.00

5

0.00

5

0.00

5

0.00

5

939110.71

14.14

963781.26

15.37

968747.12

13.89

968747.41

4.55

0.00

5

0.00

5

0.00

5

0.00

5

606822.73

14.15

615247.99

16.43

618170.78

18.61

619704.92

5.68

0.00

5

0.00

3

0.00

1

0.00

5

731857.44

14.01

749377.80

12.56

751689.75

20.31

751736.11

6.04

0.00

2

0.00

2

0.00

4

0.00

5

828223.76

12.04

849488.45

12.26

856918.94

17.79

856956.89

5.71

0.00

5

0.00

3

0.00

2

0.00

2

906333.49

12.74

933544.51

14.19

947749.84

17.27

950994.70

3.99

2

0.00

2

0.00

4

0.00

4

0.00

5

16.47

529385.67

17.34

538379.41

19.09

540666.63

19.36

541609.30

6.87

2

0.00

2

0.00

3

0.00

5

0.00

2

14.94

668479.11

14.90

685959.90

18.59

691650.49

16.69

693039.20

4.87

3

0.00

5

0.00

2

0.00

4

0.00

2

0.00

2

11.73

757265.98

13.54

782009.47

17.23

806752.95

13.70

816929.99

14.81

821577.20

5.53

3

0.00

3

0.00

5

0.00

4

0.00

5

0.00

3

16.18

846138.12

11.78

879333.95

13.76

912157.71

12.17

928125.96

16.46

935014.05

4.31

5

0.00

2

0.00

2

0.00

3

0.00

4

0.00

15.33

15.24

14.93

14.67

15.32

2 4.78

CE

CP U (s)

4

5

13.45

CPU(s)

16.94

AN US

4

Hybrid (ψ = 0.6)

CPU(s)

M

3

Hybrid (ψ = 0.4)

OF

ED

2

α

PT

p

Tables 4 and 5 present the results for solving the hose and hybrid models with the AP data set for small and medium (|N | ≤ 75) and large (|N | ≥ 100) instances, respectively. All the instances with up to 50 nodes, 325

AC

have already been solved to optimality using the Benders algorithm in [21]. However, for the instances with more than 50 nodes, the problem is solved for only three values of p as 3, 4, and 5 using the Benders algorithm.

Furthermore, in case of the Hybrid model, the instances with 200 are not solved because of long computational times. Therefore, the %GAP values are only reported for the instances that have been solved in [21] in Tables 4 and 5. The presented results show that the proposed TS algorithm is able to find the optimal solutions for the hose

330

and hybrid models with the AP data set. From a solution time perspective, one can observe that the large instances of the hybrid model from the AP data set are solved in less than four minutes, on average. For the

19

ACCEPTED MANUSCRIPT

hose model, the solution times are even smaller and the average CPU time for this model is less than one minute. Table 4: Results for the problem under hose and hybrid uncertainty with the AP data set (|N | ≤ 75) Hybrid (ψ = 0.2)

75

Hybrid (ψ = 0.8)

Hybrid (ψ = 1)

Hybrid (ψ = 2)

Hose

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

2

165060.80

0.31

168819.02

0.27

172577.24

0.30

176335.46

0.30

180093.68

0.28

187247.20

0.26

203814.57

0.16

0.00

5

0.00

5

0.00

5

0.00

5

0.00

3

147422.11

0.26

151519.33

0.27

155317.65

0.31

158889.51

0.22

162461.37

0.00

5

0.00

5

0.00

3

0.00

1

0.00

4

133170.64

0.30

137000.41

0.20

140830.19

0.24

144659.96

0.20

148172.11

0.00

5

0.00

2

0.00

5

0.00

2

0.00

5

119026.66

0.19

122661.84

0.19

126292.66

0.17

129914.99

0.25

133537.32

0.00

2

0.00

3

0.00

5

0.00

3

0.00

2

171404.41

1.26

175697.36

1.22

179990.31

1.47

184283.25

1.22

188576.20

0.00

5

0.00

5

0.00

5

0.00

5

0.00

3

153747.16

0.97

157672.41

1.00

161597.66

1.22

165522.92

1.20

169448.17

0.00

5

0.00

5

0.00

4

4

139463.84

0.82

143129.51

0.98

146795.19

0.95

0.00

5

0.00

2

0.00

5

5

129982.82

0.90

133609.26

0.83

137235.70

0.89

0.00

2

0.00

1

0.00

3

2

173131.97

2.75

177272.91

2.63

181413.84

2.47

0.00

5

0.00

5

0.00

5

3

155228.15

1.83

159126.30

2.78

163024.46

1.82

0.00

3

0.00

5

0.00

3

4

140720.60

2.12

144354.06

1.90

147987.53

1.87

0.00

2

0.00

3

0.00

5

130029.85

2.04

133816.84

1.55

137577.93

0.00

5

0.00

2

0.00

2

176713.81

13.38

181152.69

14.40

185591.56

-

4

-

3

158616.03

10.22

162688.30

0.00

4

0.00

4

143378.36

10.64

147155.22

0.00

5

0.00

2

0.00

5

0.00

4

0.00

4

0.00

2

0.00

2

5

133621.16

10.90

137394.13

9.84

141167.10

11.87

144940.08

10.64

148713.05

10.22

155507.19

10.67

170306.35

5.60

5

0.00

4

0.00

5

0.00

5

0.00

2

0.00

5

0.00

0.00

3.68

0.00

5

0.00

CR IP T

0.00

5

0.00

150460.86

0.82

154126.54

168723.49

0.33

182598.15

CPU(s)

5 0.12

5

0.00

5

0.00

5

0.19

154566.97

0.22

166162.91

0.11

4

0.00

3

0.00

3

0.21

139672.40

0.21

152274.07

0.13

5

0.00

5

0.00

5

1.26

196166.30

1.28

209111.18

0.48

5

0.00

5

0.00

5

1.20

176035.95

1.26

189952.43

0.46

5

0.00

5

0.00

5

0.88

160622.60

0.81

176189.84

0.57

0.00

5

0.00

5

0.00

3

0.00

3

140862.13

0.95

144488.57

0.76

150883.31

0.89

165649.37

0.55

0.00

5

0.00

2

0.00

2

0.00

2

185554.78

2.41

189695.72

2.44

197309.05

2.29

211318.98

1.14

0.00

5

0.00

5

0.00

5

0.00

5

166922.61

2.17

170820.77

2.25

177595.47

2.08

191842.19

1.15

0.00

5

0.00

5

0.00

5

0.00

5

151620.99

1.68

155254.45

2.13

161910.24

2.05

177383.68

0.94

3

0.00

2

0.00

2

0.00

1

0.00

6

2.18

141339.02

2.17

145100.10

1.88

151722.01

1.71

166131.78

1.43

1

0.00

3

0.00

2

0.00

4

0.00

3

13.44

189777.48

15.32

193811.57

13.20

201768.16

13.07

215849.01

3.30

5

-

5

-

5

-

5

-

5

-

5

10.69

166760.57

11.94

170832.83

13.90

174905.10

14.91

181884.48

11.76

196368.51

3.73

5

0.00

2

0.00

5

0.00

5

0.00

4

0.00

5

9.18

150932.07

10.85

154708.93

10.24

158485.78

10.49

165109.53

9.77

181077.10

3.61

3.62

3.87

3.98

AC

CE

CP U (s)

5 0.21

AN US

50

Hybrid (ψ = 0.6)

CPU(s)

M

40

Hybrid (ψ = 0.4)

OF

ED

25

α

PT

p

20

3.91

3.67

3 1.47

ACCEPTED MANUSCRIPT

Table 5: Results for the problem under hose and hybrid uncertainty with the AP data set (|N | ≥ 100) Hybrid (ψ = 0.2)

150

200

Hybrid (ψ = 0.8)

Hybrid (ψ = 1)

Hybrid (ψ = 2)

Hose

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

2

177186.48

33.38

181474.44

35.70

185762.40

35.80

190050.37

40.45

194338.33

37.64

202290.52

39.30

217300.63

9.83

0.00

4

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

3

158994.51

29.80

163006.54

34.01

167018.57

31.86

171030.60

35.82

175042.62

28.90

181994.20

32.41

196754.67

9.47

0.00

5

0.00

4

0.00

4

0.00

5

0.00

5

0.00

5

0.00

1

4

144217.26

32.00

147958.05

32.45

151698.84

29.72

155439.64

23.76

159180.43

27.31

166020.92

26.84

181884.09

11.52

0.00

5

0.00

2

0.00

1

0.00

2

0.00

5

135171.84

37.50

138973.88

26.22

142764.70

32.45

146518.15

27.70

150271.61

0.00

2

0.00

2

0.00

4

0.00

3

0.00

2

177656.21

87.85

181911.00

87.84

186165.79

80.79

190420.58

86.05

194675.38

0.00

5

0.00

4

0.00

5

0.00

5

0.00

3

159417.26

78.29

163452.67

72.08

167488.07

71.48

171523.48

85.24

175538.13

0.00

2

0.00

5

0.00

5

0.00

5

0.00

4

144632.27

66.57

148396.25

67.02

152160.22

55.44

155924.20

56.79

159681.40

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

135542.63

70.42

139343.37

73.11

143144.12

96.67

146944.86

70.52

150745.61

0.00

5

0.00

2

0.00

4

2

178241.48

174.82

182556.63

175.33

186871.77

186.59

-

5

-

5

-

5

3

159983.60

195.08

164032.32

141.64

168081.04

163.53

0.00

5

0.00

5

0.00

3

4

145198.98

131.85

148975.27

138.93

152751.56

159.88

0.00

5

0.00

3

0.00

5

5

136037.51

137.41

139833.66

140.78

143629.82

125.98

0.00

2

0.00

4

0.00

3

2

180445.86

563.48

184728.93

643.99

189012.00

589.69

-

5

-

5

-

3

162018.57

509.06

166148.88

521.50

170279.19

-

5

-

5

-

4

146753.59

522.53

150587.62

513.85

154421.65

-

1

-

5

139058.88

543.00

143209.66

5

-

-

200.82

5

0.00

3

0.00

2

157031.30

31.94

172098.88

11.90

2

0.00

5

0.00

5

83.10

202645.51

77.56

217967.72

22.66

5

0.00

5

0.00

5

75.36

182590.08

69.54

197275.77

25.71

5

0.00

4

0.00

5

78.60

166538.70

69.83

182518.12

23.01

3

0.00

5

0.00

5

67.68

157452.82

69.84

172420.17

24.12

0.00

4

0.00

2

0.00

3

0.00

2

191186.91

183.34

195487.73

190.34

203638.50

176.59

219010.32

46.73

-

5

-

5

-

5

-

5

172129.76

147.66

176178.48

150.15

183335.89

170.64

198361.42

36.96

0.00

5

0.00

5

0.00

5

0.00

4

156527.85

141.71

160304.14

140.25

167163.01

156.64

183373.34

43.27

0.00

5

0.00

5

0.00

5

0.00

3

147425.97

164.30

151222.12

123.90

158001.70

138.72

173381.56

39.80

0.00

2

0.00

3

0.00

5

0.00

2

193295.08

612.34

197526.10

567.85

205085.38

539.28

219688.55

132.19

5

-

5

-

5

-

5

-

5

519.21

177035.89

477.37

178538.25

516.99

185873.16

592.44

199944.64

125.84

185433.91

171.75

176171.21

148.43

5

-

2

-

5

-

5

443.99

158255.68

624.84

162089.71

611.38

169026.31

665.45

5

-

5

-

5

-

3

-

2

459.45

146939.90

518.02

150413.01

661.14

154335.38

516.90

160943.16

471.69

2

-

3

-

2

-

1

-

197.74

196.32

214.94

202.85

4

4

3

2

208.04

55.20

PT

CP U (s)

CPU(s)

29.17

CR IP T

125

Hybrid (ψ = 0.6)

CPU(s)

AN US

100

Hybrid (ψ = 0.4)

OF

M

α

ED

p

5.3. Results for the budget uncertainty case Tables 6 and 7 present the results obtained by solving the problem under the budget uncertainty set with

CE

335

the CAB data set for ρ=0.2 and 0.5, respectively. Since these instances are not previously solved in the

AC

literature, we have solved the problem based on the proposed MIP model using CPLEX. Therefore, in this case, the %GAP values are calculated based on the optimal solutions obtained by CPLEX. Note that in case of Γi = 0, the problem is equivalent to the deterministic UMApHMP with the nominal values for the O/D

340

demands where these parameters are not allowed to deviate from their nominal values. The results reported in Tables 6 and 7 reveal that the proposed TS algorithm is able to obtain the optimal solutions for all the instances of the CAB data set. From a solution time perspective, it is shown that the TS solves all the problem instances in very short CPU times. Note that for any combination of the parameters ρ, p, and α, as the value of the uncertainty budget (Γi ) rises, the optimal objective value increases. This

345

is due to the fact that the cost of immunizing the system against a higher uncertainty level (larger number 21

ACCEPTED MANUSCRIPT

of uncertain parameters) is greater than the corresponding cost for a lower uncertainty level. Furthermore, it can be seen that when the radius of uncertainty intervals (ρ) increase, the optimal objective values get larger as the uncertain parameters can take their values from wider uncertainty intervals. In addition, as the number of opened hubs (p) increases or as the value of the discount factor reduces, the total transportation cost decreases.

4

Γi = 2

Γi = 5

Γi = 10

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

0.2

996.02

0.31

1035.21

0.29

1057.08

0.27

1100.67

0.00

5

0.00

5

0.00

5

0.00

0.4

1072.49

0.26

1113.61

0.27

1136.98

0.28

1184.97

0.00

5

0.00

5

0.00

5

0.00

0.6

1137.08

0.25

1180.11

0.25

1205.68

0.26

1258.24

0.00

5

0.00

5

0.00

4

0.00

0.8

1180.02

0.28

1225.09

0.25

1253.03

0.26

1308.25

0.00

3

0.00

0.2

752.91

0.24

779.01

0.00

5

0.00

0.4

859.64

0.24

888.49

0.00

5

0.00

0.6

949.23

0.29

981.43

0.00

5

0.00

0.8

1020.04

0.25

1056.32

0.00

5

0.00

0.2

618.48

0.24

636.60

0.00

5

0.4

754.49

0.25

0.00

5

0.6

866.45

0.25

0.00

1

0.00

2

0.00

2

0.00

5

0.00

5

0.8

951.76

0.27

986.56

0.25

1010.02

0.27

1052.96

0.28

1091.28

0.31

0.00

2

0.00

5

0.00

3

0.00

3

0.00

3

0.2

530.00

0.26

545.02

0.23

555.35

0.24

578.76

0.30

603.67

0.35

0.00

3

0.00

2

0.00

3

0.00

5

0.00

3

0.4

676.34

0.26

697.40

0.28

711.22

0.28

740.76

0.27

771.64

0.31

0.00

3

0.00

3

0.00

5

0.00

2

0.00

2

0.6

804.70

0.28

832.06

0.26

849.98

0.35

885.12

0.28

919.70

0.28

0.00

2

0.00

3

0.00

2

0.00

4

0.00

5

0.8

910.35

0.30

944.14

0.24

966.13

0.27

1006.39

0.29

1042.99

0.27

5

0.00

3

0.00

4

0.00

3

0.00

CE

PT

5

Γi = 1

OF

CP U (s)

0.00

0.00

0.26

5

0.00

5

0.00

0.26

794.78

0.26

826.22

778.06

CPU(s)

OF

SR

%GAP

SR

0.29

1142.45

0.34

CPU(s)

5

0.00

5

0.28

1230.00

0.31

5

0.00

5

0.29

1305.13

0.32

4

0.00

5

0.28

1355.80

0.30

5

0.00

5

0.30

860.15

0.31

AN US

3

Γi = 0 (deterministic)

5

0.00

5

0.00

5

0.00

5

0.24

906.46

0.28

943.85

0.27

981.90

0.29

5

0.00

5

0.00

5

0.00

5

0.25

1002.82

0.26

1044.87

0.29

1085.55

0.29

5

0.00

4

0.00

5

0.00

5

0.25

1080.83

0.25

1126.81

0.29

1168.66

0.28

5

0.00

5

0.00

4

0.00

5

0.28

648.35

0.29

675.27

0.27

704.44

0.31

1

0.00

5

0.00

5

0.00

5

0.24

793.29

0.25

826.39

0.29

860.83

0.32

M

2

α

0.00

5

0.00

5

0.00

5

0.00

5

896.46

0.28

916.13

0.30

955.03

0.30

991.63

0.27

ED

p

CR IP T

Table 6: Results for the problem under budget uncertainty with the CAB data set (ρ = 0.2)

0.26

0.27

AC

350

22

0.28

5 0.30

ACCEPTED MANUSCRIPT

Table 7: Results for the problem under budget uncertainty with the CAB data set (ρ = 0.5)

5

Γi = 2

Γi = 5

Γi = 10

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

0.2

996.02

0.33

1093.99

0.28

1148.67

0.28

1257.63

0.36

1362.10

0.44

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

0.4

1072.49

0.27

1175.28

0.27

1233.71

0.30

1353.70

0.34

1466.26

0.43

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

0.6

1137.08

0.30

1244.66

0.29

1308.57

0.29

1439.97

0.34

1557.21

0.36

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

0.8

1180.02

0.31

1292.68

0.37

1362.55

0.29

1500.60

0.33

1619.46

0.36

0.00

5

0.00

5

0.00

4

0.00

0.2

752.91

0.27

818.16

0.26

856.95

0.29

935.55

0.00

5

0.00

5

0.00

5

0.00

0.4

859.64

0.28

931.77

0.26

976.70

0.30

1070.16

0.00

5

0.00

5

0.00

5

0.00

0.6

949.23

0.30

1029.74

0.26

1083.21

0.29

1188.32

0.00

5

0.00

5

0.00

3

0.00

0.8

1020.04

0.25

1110.74

0.28

1172.02

0.30

1286.96

0.00

5

0.00

5

0.00

5

0.00

0.2

618.48

0.37

663.79

0.35

693.15

0.30

760.45

0.00

5

0.00

0.4

754.49

0.31

813.41

0.00

2

0.00

0.6

866.45

0.26

941.48

0.00

4

0.00

0.8

951.76

0.32

1038.77

0.00

2

0.00

0.2

530.00

0.28

567.54

0.00

4

0.00

0.4

676.34

0.26

729.00

0.00

5

0.6

804.70

0.27

0.00

2

0.8

910.35

0.28

0.00 CP U (s)

4

0.29

CPU(s)

CR IP T

4

Γi = 1

OF

5

0.00

5

0.36

1020.36

0.36

5

0.00

4

0.28

1165.29

0.32

5

0.00

5

0.33

1290.03

0.30

2

0.00

5

0.33

1391.59

0.36

1

0.00

5

0.31

833.38

0.32

AN US

3

Γi = 0 (deterministic)

5

0.00

2

0.00

5

0.00

2

0.27

851.50

0.34

934.20

0.26

1020.34

0.31

1

0.00

5

0.00

5

0.00

5

0.30

990.65

0.31

1086.60

0.27

1179.41

0.29

3

0.00

5

0.00

3

0.00

5

0.27

1097.42

0.30

1204.77

0.29

1300.56

0.27

4

0.00

5

0.00

4

0.00

5

0.30

593.37

0.31

651.91

0.27

714.18

0.31

5

0.00

2

0.00

3

0.00

5

0.24

763.55

0.29

837.39

0.26

914.58

0.30

M

2

α

0.00

4

0.00

2

0.00

3

0.00

5

873.10

0.33

917.89

0.27

1005.74

0.27

1092.19

0.29

0.00

5

0.00

3

0.00

5

0.00

2

994.48

0.30

1049.45

0.30

1150.12

0.27

1241.60

0.27

2

0.00

3

0.00

4

0.00

ED

p

0.00

0.29

0.30

0.30

2 0.33

PT

The results for solving the problem under the budget uncertainty with the TR data set are presented in Tables 8 and 9. We note that as the instances from the TR data set are large (|N | = 81), CPLEX is not

CE

able to solve these instances because of excessive memory requirements that makes it impossible to solve such instances using the available resources. For this reason, we only report the solutions obtained by the proposed 355

TS algorithm and no %GAP values are reported.

AC

Note from Tables 8 and 9 that the proposed matheuristic algorithm is able to solve all the instances in short computational times. The average time required for solving these instances is less than 10 seconds and it can be verified that the algorithm has found the optimal solutions for all the instances in deterministic case (Γi = 0) as the optimal objective values for the deterministic problem are reported in [21]. As in case of the

360

CAB data set, here also it can be observed that the total cost associated with the larger values of the budget of uncertainty is greater than the total cost for smaller values of the budget of uncertainty. Furthermore, the change in the value of the objective function value as a result of altering the values of other input parameters, such as p, α, and ρ, is similar to the case of the CAB data set.

23

ACCEPTED MANUSCRIPT

Table 8: TS results for the problem under budget uncertainty with the TR data set (ρ = 0.2)

4

5

CPU(s)

%GAP

SR

0.2

770985.04

6.81

0.00

4

0.4

820586.50

6.26

0.00

5

0.6

857219.51

6.11

0.00

5

0.8

878672.80

5.66

0.00

5

0.2

660218.05

6.34

0.00

1

0.4

726196.77

5.82

0.00

5

0.6

778077.05

5.86

0.00

5

0.8

814313.29

5.65

0.00

5

0.2

570217.55

6.62

0.00

2

0.4

657662.12

6.97

0.00

3

0.6

729447.41

5.95

0.00

2

0.8

777778.51

6.88

0.00

2

0.2

492494.33

8.82

0.00

3

0.4

595161.93

6.85

0.00

5

0.6

678419.46

6.13

0.00

2

0.8

744056.84

9.14

0.00

Γi = 2

CPU(s)

OF

SR 787277.11

6.92

838260.81

5.92

876368.28

6.49

898983.00

6.69

674697.71

6.59

742672.01

6.28

796461.29

7.21

834217.37

6.17

580921.12

6.01

671458.95

7.76

746070.28

7.92

797590.35

6.00

504405.71

8.81

607937.72

7.45

694351.14

6.51

763294.56

7.55

6.03

847808.51

6.58

886412.07

6.34

909044.62

6.74

681599.52

6.75

751332.83

6.32

805465.13

7.91

843545.21

6.36

585728.99

7.35

677317.26

7.37

753064.35

7.40

805976.38

6.38

506073.61

11.49

613296.11

8.46

700893.26

7.18

771038.04

8.12

865850.29

7.45

905177.23

6.60

928443.82

6.71

696572.46

6.75

766237.54

6.64

821598.87

8.10

860861.41

7.46

5

5

3

5

3

5

5

5

1

5

3

4

2

2

3

5

2

3

1

5

598235.78

8.01

691577.66

7.32

771293.30

9.09

823365.76

7.34

516845.64

10.29

626045.86

7.64

833408.09

8.35

887160.01

7.49

927455.36

7.50

715223.56

8.63

786930.38

9.05

849457.25

7.96

904176.01

7.95

OF

CPU(s)

862762.04

8.43

918134.53

9.01

945140.56

8.69

959553.78

8.77

SR

5

5

5

5

4

951230.59

7.59

713723.54

7.32

5

4

785610.31

8.37

842670.87

7.35

882842.86

7.19

969353.28

8.36

727300.48

9.32

800705.16

7.61

859045.50

7.61

899809.80

8.13

626259.00

9.53

723588.93

7.82

803397.40

7.80

858673.27

8.10

539647.70

10.80

653318.99

10.50

746803.92

7.53

821746.42

7.45

5

709572.27

8.54

788172.83

7.84

842852.92

8.51

529708.86

8.18

641081.72

8.12

732568.56

9.79

806319.15

9.35

8.63

913117.52

9.06

636014.47

9.22

734803.48

10.44

818307.79

9.38

872326.84

8.67

548047.43

9.39

663698.69

9.63

758323.54

10.77

834166.84

7.58

4

2

3

2

5

5

5

3

3

5

1

3

2

3

5

5

2

3

5

871881.91

3

3

3

8.45

3

2

2

813018.30

5

5

1

9.56

5

5

8.97

738872.30

5

5

614133.29

8.03

2

4

5

984096.35

5

2

2

2

Γi = 20

CPU(s) SR

5

3

3

OF

SR

2

3

5

Γi = 15

CPU(s)

4

5

5

OF

4

5

5

1

7.81

5

5

6.62

813538.54

5

5

Γi = 10

CPU(s) SR

796233.68

5

2

OF

SR

5

2

Γi = 5

CPU(s)

2

3

2

3

3

3

4

3

2

6.89

7.30

7.81

8.15

8.45

9.06

PT

CP U (s)

Γi = 1 OF

CR IP T

OF

AN US

3

Γi = 0 (deterministic)

M

2

α

ED

p

Tables 10 - 13 present the results for solving the problem under the budget uncertainty for the AP data set. The instances of up to 50 nodes (|N | ≤ 50) are solved using both CPLEX and the proposed TS algorithm.

CE

365

It is shown that the TS has solved all the problem instances with up to 50 nodes to optimality. The larger

AC

instances (|N | ≥ 75) are not solved by CPLEX due to excessive memory requirements. However, as can be seen from Tables 10 and 11, the proposed TS algorithm is able to solve all the small and medium instances in average CPU time of around 2 seconds. On the other hand, the results reported in Tables 12 and 13 show that

370

the large instances of the AP data set are solved within very short times as the average solution time for these instances is less one minute. It can also be verified (based on the results from [21]) that the TS has found the optimal solutions for all the instances in deterministic case (Γi = 0). Therefore, it can be concluded that

the proposed matheuristic algorithm is able to solve all the instances of up to 200 nodes under three different types of polyhedral uncertainty sets very efficiently. 375

24

ACCEPTED MANUSCRIPT

Table 9: TS results for the problem under budget uncertainty with the TR data set (ρ = 0.5)

4

5

CPU(s)

%GAP

SR

0.2

770985.04

6.93

0.00

5

0.4

820586.50

6.19

0.00

5

0.6

857219.51

5.96

0.00

5

0.8

878672.80

6.70

0.00

5

0.2

660218.05

7.11

0.00

5

0.4

726196.77

6.65

0.00

5

0.6

778077.05

7.25

0.00

4

0.8

814313.29

6.43

0.00

5

0.2

570217.55

9.21

0.00

2

0.4

657662.12

8.89

0.00

4

0.6

729447.41

5.75

0.00

2

0.8

777778.51

5.78

0.00

5

0.2

492494.33

6.41

0.00

5

0.4

595161.93

9.88

0.00

5

0.6

678419.46

8.09

0.00

3

0.8

744056.84

7.64

0.00

Γi = 2

CPU(s)

OF

Γi = 5

CPU(s)

SR 811153.20

6.20

864772.27

6.04

905091.43

6.22

929448.31

6.22

692959.74

6.50

767384.85

6.67

824037.66

6.30

864073.48

7.13

596976.47

7.64

692154.20

6.20

771004.58

6.83

827292.95

7.20

516029.14

7.79

627101.41

7.52

718248.65

8.83

SR 833544.62

6.53

888641.52

6.57

930200.91

6.45

954602.34

6.33

711720.19

6.42

789036.92

6.65

846547.26

6.59

887393.08

7.36

608996.15

8.23

706799.97

8.04

788489.75

7.33

848957.81

6.63

526442.51

7.89

5

6.61

1003100.35

6.71

750783.04

7.65

826298.68

7.22

886881.59

6.65

930683.58

7.23

4

3

5

2

2

5

5

5

1

4

5

5 640263.13

10.26

742450.97

8.15

827340.88

7.93

891077.73

7.89

553372.60

12.02

640497.38

8.59

734603.97

7.75

810765.05

987020.28

8.39

1032809.14

7.35

8.56

850449.15

8.44

9.19

1029560.28

8.08

OF

CPU(s)

999865.52

9.15

1064456.58

8.99

1077022.13

8.37

1113055.19

8.81

SR

5

5

5

5 1060067.27

7.20

793562.59

8.00

5

4

874730.62

8.73

938760.56

7.45

985637.21

7.88

1105374.01

8.25

827504.59

9.22

912467.75

8.42

978638.28

8.22

1028054.57

7.95

710299.60

10.39

822479.14

9.77

913481.22

10.23

980015.41

7.61

610240.81

9.68

740554.58

8.64

849380.59

9.97

937837.48

8.19

5

5

9.77

875419.79

6.89

940464.52

8.63

585530.64

8.49

709961.41

8.44

813792.22

7.95

899074.21

9.37

8.04

1010845.67

7.43

1061323.86

9.02

733007.05

8.09

850515.51

8.00

944209.10

10.33

1011915.95

8.75

631377.07

11.96

766503.84

8.57

878179.65

8.17

969081.79

8.44

5

3

5

2

4

3

5

5

3

2

5

4

3

2

5

2

2

2

5

5

1

5

3

943250.60

5

3

3

8.79

2

4

787437.49

856280.17

5

5

7.84

8.58

4

5

680006.90

1142231.67

5

5

4

770429.70

966603.53

5

4

8.75

Γi = 20

CPU(s) SR

5

4

2

9.23

7.79

3

3

2

926480.64

2

672371.76

OF

SR

5

5

Γi = 15

CPU(s)

5

5

4

5

977113.80

4

5

OF

5

2

2

7.18

7.38

5

3

6.75

933745.98

5

5

791469.70

7.70

5

5

5

876806.76

5

5

Γi = 10

CPU(s) SR

5

5

5

OF

2

5

3

5

3

4

5

2

3

6.88

7.29

8.07

8.13

8.89

8.82

PT

CP U (s)

Γi = 1 OF

CR IP T

OF

AN US

3

Γi = 0 (deterministic)

M

2

α

ED

p

To have a better visualization on the effect of the budget of uncertainty (Γi ) and the radius of flow

CE

parameter intervals (ρ) on the total transportation cost of the network, we plot the total cost for different values of Γi and ρ with an instance of 100 nodes (|N |=100) from the AP data set. Two values are used for the radius of interval as ρ ∈ {0.2, 0.5}. However, since for each node i ∈ N , we have |N | − 1 different outgoing flow values that can change simultaneously, we selected 11 different values for the budget of uncertainty (Γi )

AC

380

varying from 0 to |N | − 1. The resulting plots are depicted in Figure 5.

25

ACCEPTED MANUSCRIPT

Table 10: TS results for the problem under budget uncertainty with the AP data set (|N | ≤ 75, ρ = 0.2)

75

Γi = 2

Γi = 5

Γi = 10

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

2

161302.58

0.22

165006.38

0.19

167986.04

0.19

174420.06

0.20

181781.09

0.21

0.00

5

0.00

5

0.00

5

0.00

5

0.00

5

3

143324.89

0.19

146558.97

0.21

149172.43

0.22

154718.02

0.20

161252.60

0.21

0.00

3

0.00

3

0.00

2

0.00

3

0.00

4

4

129326.76

0.19

132146.31

0.19

134397.61

0.18

139262.74

0.19

145213.41

0.25

0.00

2

0.00

5

0.00

2

0.00

3

0.00

5

5

115391.48

0.20

117947.73

0.16

119868.69

0.21

124161.48

0.18

129423.24

0.20

0.00

4

0.00

3

0.00

3

0.00

2

167111.47

0.72

170387.27

0.68

172577.36

0.73

177401.75

0.00

4

0.00

5

0.00

5

0.00

3

149821.91

0.72

152765.47

0.82

154709.88

0.82

159070.47

0.00

5

0.00

5

0.00

5

0.00

4

135798.16

0.83

138498.53

0.94

140231.85

1.11

144157.38

0.00

5

0.00

4

0.00

2

0.00

5

126356.39

0.83

128844.10

0.81

130438.67

0.81

134004.64

0.00

2

0.00

3

0.00

2

0.00

2

168991.03

1.32

172107.84

1.45

174104.78

1.40

178520.54

0.00

5

0.00

5

3

151329.99

1.32

154123.77

1.64

0.00

5

0.00

5

4

137087.13

1.67

139639.30

1.60

0.00

5

0.00

4

5

126236.27

1.71

128542.71

1.71

0.00

4

0.00

2

2

172274.94

4.78

175114.15

5.14

-

5

-

5

3

154543.76

5.57

157118.85

5.01

-

3

-

5

-

2

-

3

-

5

4

139601.50

5.86

141924.47

5.88

143357.05

5.49

146278.05

7.23

149758.47

6.29

-

2

-

5

-

5

-

4

-

3

5

129848.19

4.93

132049.34

6.05

133279.33

6.49

135901.69

5.67

139105.71

6.15

2

-

2

-

3

-

5

-

-

1.94

2

0.00

2

0.75

183325.74

0.82

5

0.00

4

0.96

164344.63

0.99

5

0.00

3

1.15

148887.39

0.85

3

0.00

2

0.83

138446.16

1.05

3

0.00

5

1.52

183760.81

1.48

0.00

5

0.00

5

0.00

5

155903.56

1.58

159786.77

1.43

164484.15

1.56

0.00

5

0.00

5

0.00

5

141293.67

1.53

144793.97

1.58

148997.51

1.68

0.00

5

0.00

3

0.00

2

130004.53

1.74

133127.82

1.75

136897.00

1.73

0.00

3

0.00

5

0.00

3

176836.01

5.37

180410.32

5.41

184701.40

5.21

-

5

-

5

-

5

158677.39

5.11

161922.06

5.32

165832.35

5.51

2.03

2.06

2.15

3 2.14

AC

CE

PT

CP U (s)

CPU(s)

CR IP T

50

Γi = 1

OF

AN US

40

Γi = 0 (deterministic)

M

25

p

ED

|N |

(a) ρ = 0.2

(b) ρ = 0.5

Figure 5: Effect of the budget of uncertainty (Γi ) and the radius of intervals (ρ) on the total cost (AP data set, |N |=100)

26

ACCEPTED MANUSCRIPT

Table 11: TS results for the problem under budget uncertainty with the AP data set (|N | ≤ 75, ρ = 0.5)

75

Γi = 2

Γi = 5

Γi = 10

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

CPU(s)

OF

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

%GAP

SR

2

161302.58

0.23

170562.07

0.19

178011.25

0.19

194096.28

0.20

212498.87

0.23

0.00

5

0.00

5

0.00

5

0.00

4

0.00

3

3

143324.89

0.19

151410.09

0.19

157943.76

0.17

171807.71

0.22

188144.16

0.23

0.00

4

0.00

5

0.00

5

0.00

4

0.00

4

4

129326.76

0.16

136375.63

0.17

142003.90

0.20

154166.71

0.19

169043.40

0.21

0.00

4

0.00

2

0.00

3

0.00

3

0.00

2

5

115391.48

0.19

121782.11

0.16

126584.51

0.28

137316.47

0.18

150470.88

0.19

0.00

5

0.00

5

0.00

5

0.00

2

167111.47

0.72

175300.97

0.75

180776.20

0.70

192837.16

0.00

5

0.00

5

0.00

4

0.00

3

149821.91

0.69

157180.80

0.72

162041.83

0.77

172943.32

0.00

4

0.00

3

0.00

5

0.00

4

135798.16

0.83

142549.09

0.78

146882.38

0.86

156696.21

0.00

2

0.00

1

0.00

5

0.00

5

126356.39

0.92

132575.66

0.86

136562.09

0.79

145477.02

0.00

5

0.00

3

0.00

2

0.00

2

168991.03

1.42

176783.06

1.52

181775.42

1.51

192660.89

0.00

5

0.00

5

3

151329.99

1.45

158314.43

1.55

0.00

5

0.00

5

4

137087.13

1.75

143467.56

1.99

0.00

2

0.00

5

5

126236.27

1.50

132002.37

2.07

0.00

3

0.00

4

2

172274.94

4.89

179372.98

5.00

-

5

-

5

3

154543.76

5.15

160981.49

5.42

-

5

-

4

-

1

-

2

-

3

4

139601.50

5.53

145408.91

5.33

148990.37

6.14

156292.87

6.00

164993.91

6.32

-

5

-

2

-

5

-

2

-

3

5

129848.19

6.45

135351.07

7.33

138426.03

6.77

144981.92

6.72

152991.98

6.79

2

-

2

-

3

-

2

-

CP U (s)

2.00

CPU(s)

CR IP T

50

Γi = 1

OF

4

0.00

4

0.79

207647.15

0.84

5

0.00

5

0.76

186128.70

0.98

3

0.00

2

0.93

168521.24

1.03

5

0.00

1

1.07

156550.64

0.95

2

0.00

5

1.52

205771.55

1.60

AN US

40

Γi = 0 (deterministic)

M

25

p

ED

|N |

0.00

5

0.00

5

0.00

5

162763.92

1.59

172471.93

1.73

184215.37

1.85

0.00

2

0.00

5

0.00

5

147603.48

1.53

156354.23

1.86

166863.06

1.97

0.00

4

0.00

4

0.00

5

135544.55

1.61

143352.79

1.94

152775.74

1.74

0.00

4

0.00

2

0.00

4

183677.62

5.03

192613.39

5.28

203341.10

5.55

-

4

-

5

-

5

164877.84

5.24

172989.52

6.21

182765.23

5.84

2.13

2.09

2.22

3 2.27

PT

Note that the case with Γi =0 corresponds to the deterministic problem in which no flow parameter is allowed to deviate from its nominal value. In contrast, the case with Γi =99 corresponds to the worst case problem where all the outgoing flow values from each node take their worst (the largest, in our case) values from the corresponding intervals of uncertainty. As noted earlier, the total transportation cost rises as the

CE

385

value of the budget of uncertainty increases. However, it can be observed from Fig. 5 that the grow in the

AC

total cost is much faster when the budget of uncertainty is low compared to higher values of the budget of uncertainty. Furthermore, one can see from this figure that the grow in the total transportation cost for the case of ρ=0.5 is substantially larger that the corresponding increase for the case of ρ=0.2.

390

The results reported in Tables 2-13 show that the proposed TS based matheuristic algorithm solves all the instances of the three problem variants very efficiently in terms of solution time and quality. It should, however, be noted that as the number of installed hubs (p) increases, the number of successful runs (SR) tend to decrease in general. This shows that as p grows, it becomes harder for the TS algorithm to find the optimal solution. Considering the optimality gap, one can observe that the percentage gap (%GAP) between

27

ACCEPTED MANUSCRIPT

Table 12: TS results for the problem under budget uncertainty with the AP data set (|N | ≥ 100, ρ=0.2)

OF

CPU(s)

Γi = 1 OF

SR 100

2

172898.52

11.67

175760.70

11.77

3

154982.48

12.66

157538.17

13.70

4

140476.46

14.54

142797.13

15.62

5

131369.79

11.57

133589.73

16.74

2

173401.42

26.76

176219.89

28.87

3

155381.86

23.50

157908.94

24.97

4

5

5

4

140868.29

26.59

143160.35

25.97

5

131741.88

26.87

133887.69

29.56

2

173926.34

39.59

176715.62

38.84

3

155934.88

39.14

158466.36

41.20

4

141422.69

40.96

143719.13

44.82

5

132241.36

54.76

134386.98

45.58

2

176016.07

101.41

178757.65

95.29

3

157888.26

96.88

160354.00

115.16

4

142919.56

115.64

145141.85

110.63

5

135274.17

111.57

137405.14

121.46

5

144157.00

15.78

134758.95

18.14

177753.62

26.18

159256.87

28.82

144412.80

25.10

135007.04

35.16

178247.52

39.65

159814.93

47.29

144972.18

51.02

135502.96

51.15

180039.09

103.01

5

161467.32

110.52

4

146187.64

114.33

1

138738.13

181.02

5

5

3

3

2

137151.85

14.26

180724.82

27.72

161977.39

29.51

146865.11

29.34

137285.15

40.19

5

2

5

5

3

2

Γi = 15

CPU(s)

181204.08

45.71

162515.68

48.05

147404.00

49.42

137758.48

66.50

182681.34

105.55

163843.10

110.89

184565.37

12.87

165339.15

12.73

149918.14

17.18

148351.99

121.90

140731.67

139.09

187654.74

13.05

168095.85

14.13

OF

CPU(s)

190176.53

14.93

170386.10

14.08

152385.64

14.22

154425.11

13.69

SR

5

5

5

5

5

140009.43

19.26

184420.81

25.38

5

5

165327.33

26.44

149859.60

25.85

139994.84

31.89

142284.99

23.60

187328.22

27.12

167958.34

29.09

152195.28

27.50

142149.39

30.40

187690.45

45.84

168366.50

48.51

152597.03

45.74

142522.76

56.21

188370.20

115.90

168881.25

125.05

152915.60

125.80

144729.00

150.44

4

165779.01

42.98

150308.36

45.56

140400.73

46.20

185873.64

111.70

166650.78

114.01

150910.97

127.06

142825.90

146.54

143961.13

29.57

190054.49

49.60

170504.16

46.05

154485.32

51.91

144284.85

50.66

190485.34

116.86

170791.49

134.74

154634.15

145.87

146527.09

133.51

5

3

2

2

5

5

5

3

5

5

5

1

2

5

5

5

3

27.43

5

2

4

154134.72

4

3

2

26.23

5

5

3

170140.33

4

5

4

27.28

4

3

46.15

189737.44

5

5

184821.86

16.80

4

3

5

144241.15

4

5

5

3

Γi = 20

CPU(s) SR

5

5

5

OF

SR

5

5

5

5

2

5

4

4

4

3

3

2

48.76

54.61

54.50

53.24

55.79

56.20

PT

47.13

17.18

5

ED

5

146747.09

5

2

OF

5

5

5

5

13.93

3

1

3

161844.87

2

5

3

12.84

5

5

4

180610.63

5

3

5

395

13.41

2

4

CP U (s)

158992.77

5

5

200

13.22

3

Γi = 10

CPU(s) SR

177383.97

5

2

OF

SR

5

2

Γi = 5

CPU(s)

5

5

150

OF

SR

5

125

Γi = 2

CPU(s)

CR IP T

Γi = 0 (deterministic)

AN US

p

M

|N |

the best found solution by the TS algorithm and the optimal objective is zero for all the instances for which

CE

the optimal solution is known. Moreover, based on the results obtained by our experiments, the average gap (over five runs of the algorithm) never exceeded three percent which is an indication of the robustness and

AC

efficiency of the proposed algorithm.

6. Conclusions

400

In this research the robust multiple allocation p-hub median problem under polyhedral demand uncertainty

was considered. It was assumed that the flow volumes between O/D pairs were of uncertain nature and the underlying uncertainty was modeled using three types of uncertainty set called the hose, the hybrid, and the budget uncertainty sets. The problems were formulated as linear mixed integer programming problems and an efficient matheuristic solution approach combining the Tabu Search (TS) metaheuristic and mathematical

28

ACCEPTED MANUSCRIPT

Table 13: TS results for the problem under budget uncertainty with the AP data set (|N | ≥ 100, ρ=0.5)

OF

CPU(s)

Γi = 1 OF

SR 100

2

172898.52

11.64

3

154982.48

12.60

4

140476.46

12.67

5

131369.79

19.27

2

173401.42

23.47

3

155381.86

26.41

4

140868.29

28.91

5

131741.88

33.44

2

173926.34

42.43

3

155934.88

42.51

4

141422.69

54.76

5

132241.36

50.69

2

176016.07

97.55

3

157888.26

111.52

4

142919.56

114.96

5

135573.93

135.22

180029.56

12.27

161371.70

12.09

146278.13

12.86

136851.60

14.73

180426.10

25.17

161699.55

24.40

146598.44

25.84

137106.41

28.05

180899.54

40.82

162263.59

44.58

147163.80

46.71

137605.40

57.20

182870.03

97.77

164052.61

121.44

148475.28

123.37

140703.66

136.40

5

149677.82

14.84

139842.70

16.63

184253.98

27.42

165069.40

28.19

149729.56

30.94

139904.77

36.26

184729.29

42.94

165635.00

43.02

150296.42

46.19

140395.36

46.06

186073.61

105.07

5

166835.89

103.56

151089.78

151.37

143057.75

161.25

5

5

5

2

2

3

4

17.12

191686.22

25.38

171870.69

27.95

155860.34

27.08

145572.59

42.22

3

3

5

5

5

2

192120.68

45.40

172386.88

50.63

156375.96

48.56

146034.16

60.23

192679.25

105.28

172775.34

106.43

156500.65

148.05

202064.02

13.77

180874.16

14.17

164016.68

14.20

148468.28

148.40

Γi = 20

CPU(s)

OF

CPU(s)

216093.55

14.35

193491.52

15.45

175306.59

14.01

SR 209789.09

14.44

187765.91

14.87

170118.21

16.58

5

SR

4

5

5

5

5

152968.88

16.15

200927.94

26.10

5

5

180245.53

27.82

163298.84

29.22

152326.66

33.88

158657.78

19.34

208189.29

28.51

186823.07

30.89

169090.88

29.84

157701.51

49.01

208284.21

51.68

187013.93

51.11

169223.84

58.55

157867.17

58.29

206901.41

118.53

185370.72

122.08

167866.92

145.56

158921.90

127.19

3

180545.20

53.40

163546.91

54.79

152621.57

53.11

200660.00

113.41

179794.56

140.90

162879.31

135.83

154157.79

132.41

162230.87

31.18

214158.48

47.55

192358.08

56.68

173992.07

60.84

162272.37

66.71

212189.25

113.19

190146.34

128.85

172206.03

117.29

162956.83

128.24

2

5

3

2

5

5

5

1

5

4

5

2

5

5

2

5

4

35.00

2

5

3

173983.46

2

3

5

34.07

5

5

3

192278.03

3

5

4

30.34

3

5

45.70

214190.13

5

2

201132.72

17.40

5

2

4

163418.38

2

5

5

1

OF

SR

5

5

5

5

1

5

3

3

2

5

3

4

51.48

54.94

55.87

56.55

58.53

56.95

PT

51.13

145860.50

Γi = 15

CPU(s)

5

5

ED

5

15.78

5

3

5

156153.04

4

3

OF

5

5

5

2

13.03

3

5

5

172138.45

5

5

1

12.41

5

3

5

192160.59

5

3

2

405

12.76

5

5

CP U (s)

165008.19

5

5

200

12.49

2

Γi = 10

CPU(s) SR

184085.71

5

3

OF

SR

5

2

Γi = 5

CPU(s)

5

3

150

OF

SR

5

125

Γi = 2

CPU(s)

CR IP T

Γi = 0 (deterministic)

AN US

p

M

|N |

programming techniques was proposed to solve the three variants of the problem.

CE

Extensive computational experiments based on three well-known data sets from the literature were conducted to analyze the efficiency and the effectiveness of the proposed matheuristic algorithm as well as to study the

AC

effect of different input parameters on the final solutions of all the three models. Computational results showed that the proposed algorithm was able to solve large-scale instances of the problem in very short computational

410

times.

This research can be extended in several directions. To this end, one may focus on applying the proposed solution algorithm on other versions of the robust HLP. Moreover, instead of using classical HLP assumptions such as flow-independent economies of scale on only inter-hub links, more realistic assumptions like flowdependent scale economies on all links can be used to better reflect the real-world situations.

29

ACCEPTED MANUSCRIPT

415

Acknowledgments The author sincerely thanks the anonymous referees for their valuable comments and suggestions that have helped improve the paper considerably.

References

420

CR IP T

References

[1] J. F. Campbell, M. E. O’Kelly, Twenty-five years of hub location research, Transportation Science 46 (2) (2012) 153–169.

[2] S. Alumur, B. Y. Kara, Network hub location problems: The state of the art, European journal of operational research 190 (1) (2008) 1–21.

425

AN US

[3] R. Z. Farahani, M. Hekmatfar, A. B. Arabani, E. Nikbakhsh, Hub location problems: A review of models, classification, solution techniques, and applications, Computers & Industrial Engineering 64 (4) (2013) 1096 – 1109.

[4] I. Contreras, Hub location problems, in: G. Laporte, S. Nickel, F. Saldanha da Gama (Eds.), Location

M

science, Springer, 2015, pp. 311–344.

[5] J. F. Campbell, Location and allocation for distribution systems with transshipments and transportion economies of scale, Annals of Operations Research 40 (1) (1992) 77–99.

ED

430

[6] J. F. Campbell, Integer programming formulations of discrete hub location problems, European Journal of Operational Research 72 (2) (1994) 387–405.

PT

[7] D. Skorin-Kapov, J. Skorin-Kapov, M. O’Kelly, Tight linear programming relaxations of uncapacitated p-hub median problems, European Journal of Operational Research 94 (3) (1996) 582 – 593. [8] A. T. Ernst, M. Krishnamoorthy, Exact and heuristic algorithms for the uncapacitated multiple allocation

CE

435

p-hub median problem, European Journal of Operational Research 104 (1) (1998) 100–112.

AC

[9] V. Marianov, D. Serra, Location models for airline hubs behaving as m/d/c queues, Computers & Operations Research 30 (7) (2003) 983–1003.

[10] T.-H. Yang, Stochastic air freight hub location and flight routes planning, Applied Mathematical Mod-

440

elling 33 (12) (2009) 4424–4430. [11] T. Sim, T. J. Lowe, B. W. Thomas, The stochastic p-hub center problem with service-level constraints, Computers & Operations Research 36 (12) (2009) 3166 – 3177.

30

ACCEPTED MANUSCRIPT

[12] R. S. de Camargo, G. de Miranda Jr., R. P. Ferreira, A hybrid outer-approximation/benders decomposition algorithm for the single allocation hub location problem under congestion, Operations Research 445

Letters 39 (5) (2011) 329–337. [13] I. Contreras, J.-F. Cordeau, G. Laporte, Stochastic uncapacitated hub location problem, European Journal of Operational Research 212 (3) (2011) 518–528.

CR IP T

[14] H. Zhai, Y. Liu, W. Chen, Applying minimum-risk criterion to stochastic hub location problems, Procedia Engineering 29 (2012) 2313–2321, 2012 International Workshop on Information and Electronics 450

Engineering.

[15] S. A. Alumur, S. Nickel, F. S. da Gama, Hub location under uncertainty, Transportation Research Part B: Methodological 46 (4) (2012) 529–543.

AN US

[16] S. Chaharsooghi, F. Momayezi, N. Ghaffarinasab, An adaptive large neighborhood search heuristic for solving the reliable multiple allocation hub location problem under hub disruptions, International Journal 455

of Industrial Engineering Computations 8 (2) (2017) 191–202.

[17] M. Shahabi, A. Unnikrishnan, Robust hub network design problem, Transportation Research Part E: Logistics and Transportation Review 70 (2014) 356–373.

M

[18] N. Ghaffari-Nasab, M. Ghazanfari, E. Teimoury, Robust optimization approach to the design of huband-spoke networks, The International Journal of Advanced Manufacturing Technology 76 (5-8) (2015) 1091–1110.

ED

460

[19] F. H. Boukani, B. F. Moghaddam, M. S. Pishvaee, Robust optimization approach to capacitated single and multiple allocation hub location problems, Computational and Applied Mathematics 35 (1) (2016)

PT

45–60.

[20] A. Ghaderi, R. Rahmaniani, Meta-heuristic solution approaches for robust single allocation p-hub median problem with stochastic demands and travel times, The International Journal of Advanced Manufacturing

CE

465

Technology 82 (9-12) (2016) 1627–1647. [21] M. Meraklı, H. Yaman, Robust intermodal hub location under polyhedral demand uncertainty, Trans-

AC

portation Research Part B: Methodological 86 (2016) 66 – 85.

[22] C. A. Zetina, I. Contreras, J.-F. Cordeau, E. Nikbakhsh, Robust uncapacitated hub location, Transporta-

470

tion Research Part B: Methodological 106 (Supplement C) (2017) 393–410.

[23] M. Meraklı, H. Yaman, A capacitated hub location problem under hose demand uncertainty, Computers & Operations Research 88 (2017) 58 – 70. [24] E. M. de S´ a, R. Morabito, R. S. de Camargo, Benders decomposition applied to a robust multiple allocation incomplete hub location problem, Computers & Operations Research 89 (2018) 31–50. 31

ACCEPTED MANUSCRIPT

475

[25] M. E. O’Kelly, A quadratic integer program for the location of interacting hub facilities, European Journal of Operational Research 32 (3) (1987) 393–404. [26] D. Skorin-Kapov, J. Skorin-Kapov, On tabu search for the location of interacting hub facilities, European Journal of Operational Research 73 (3) (1994) 502–509. [27] A. T. Ernst, M. Krishnamoorthy, Efficient algorithms for the uncapacitated single allocation p-hub median problem, Location science 4 (3) (1996) 139–154.

CR IP T

480

[28] A. T. Ernst, M. Krishnamoorthy, Solution algorithms for the capacitated single allocation hub location problem, Annals of Operations Research 86 (1999) 141–159.

[29] S. Abdinnour-Helm, Using simulated annealing to solve the p-hub median problem, International Journal of Physical Distribution & Logistics Management 31 (3) (2001) 203–220.

[30] H. Topcuoglu, F. Corut, M. Ermis, G. Yilmaz, Solving the uncapacitated hub location problem using

AN US

485

genetic algorithms, Computers & Operations Research 32 (4) (2005) 967–984. [31] J.-F. Chen, A hybrid heuristic for the uncapacitated single allocation hub location problem, Omega 35 (2) (2007) 211–220.

[32] M. R. Silva, C. B. Cunha, New simple and efficient heuristics for the uncapacitated single allocation hub location problem, Computers & Operations Research 36 (12) (2009) 3152–3165.

M

490

ED

[33] H. Calık, S. A. Alumur, B. Y. Kara, O. E. Karasan, A tabu-search based heuristic for the hub covering problem over incomplete hub networks, Computers & Operations Research 36 (12) (2009) 3088 – 3096. [34] M. S. Jabalameli, F. Barzinpour, A. Saboury, N. Ghaffari-Nasab, A simulated annealing-based heuristic

495

PT

for the single allocation maximal covering hub location problem, International Journal of Metaheuristics 2 (1) (2012) 15–37.

CE

[35] R. Abyazi-Sani, R. Ghanbari, An efficient tabu search for solving the uncapacitated single allocation hub location problem, Computers & Industrial Engineering 93 (2016) 99 – 109.

AC

[36] K. Yang, L. Yang, Z. Gao, Planning and optimization of intermodal hub-and-spoke network under mixed uncertainty, Transportation Research Part E: Logistics and Transportation Review 95 (2016) 248 – 266.

500

[37] M. R. Silva, C. B. Cunha, A tabu search heuristic for the uncapacitated single allocation p-hub maximal covering problem, European Journal of Operational Research 262 (3) (2017) 954 – 965.

[38] J. F. Campbell, Hub location and the p-hub median problem, Operations Research 44 (6) (1996) 923–935. [39] V. Marianov, D. Serra, C. ReVelle, Location of hubs in a competitive environment, European Journal of Operational Research 114 (2) (1999) 363–371.

32

ACCEPTED MANUSCRIPT

505

[40] N. Boland, M. Krishnamoorthy, A. T. Ernst, J. Ebery, Preprocessing and cutting for multiple allocation hub location problems, European Journal of Operational Research 155 (3) (2004) 638–653. [41] A. L¨ uer-Villagra, V. Marianov, A competitive hub location and pricing problem, European Journal of Operational Research 231 (3) (2013) 734–744. [42] J. Sender, T. Siwczyk, P. Mutzel, U. Clausen, Matheuristics for optimizing the network in German wagonload traffic, EURO Journal on Computational Optimization (2016) 1–26.

CR IP T

510

[43] N. Ghaffarinasab, A. Motallebzadeh, Y. Jabarzadeh, B. Y. Kara, Efficient simulated annealing based solution approaches to the competitive single and multiple allocation hub location problems, Computers & Operations Research 90 (Supplement C) (2018) 173 – 192.

[44] N. Ghaffarinasab, A. Motallebzadeh, Hub interdiction problem variants: Models and metaheuristic solution algorithms, European Journal of Operational Research 267 (2) (2018) 496 – 512.

AN US

515

[45] H. W. Hamacher, M. Labb´e, S. Nickel, T. Sonneborn, Adapting polyhedral properties from facility to hub location problems, Discrete Applied Mathematics 145 (1) (2004) 104 – 116. [46] V. Maniezzo, T. St¨ utzle, S. Voß, Matheuristics: Hybridizing Metaheuristics and Mathematical Programming, Annals of Information Systems, Springer US, 2009.

[47] G. H. Fonseca, H. G. Santos, E. G. Carrano, Integrating matheuristics and metaheuristics for timetabling,

M

520

Computers & Operations Research 74 (2016) 108 – 117.

ED

[48] S.-W. Lin, K.-C. Ying, Optimization of makespan for no-wait flowshop scheduling problems using efficient matheuristics, Omega 64 (2016) 115 – 125.

525

PT

[49] P. Grangier, M. Gendreau, F. Lehud, L.-M. Rousseau, A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking, Computers & Operations Research 84 (2017) 116 –

CE

126.

[50] L. Fanjul-Peyro, F. Perea, R. Ruiz, Models and matheuristics for the unrelated parallel machine scheduling

AC

problem with additional resources, European Journal of Operational Research 260 (2) (2017) 482 – 493. [51] Y. Kochetov, N. Kochetova, A. Plyasunov, A matheuristic for the leader-follower facility location and

530

design problem, in: Proceedings of the 10th Metaheuristics International Conference (MIC 2013), Vol. 32, Citeseer, 2013.

[52] D. Aksen, N. Aras, A matheuristic for leader-follower games involving facility location-protectioninterdiction decisions, in: Metaheuristics for Bi-level Optimization, Springer, 2013, pp. 115–151. [53] F. Stefanello, O. C. B. de Ara´ ujo, F. M. M¨ uller, Matheuristics for the capacitated p-median problem, 535

International Transactions in Operational Research 22 (1) (2015) 149–167. 33

ACCEPTED MANUSCRIPT

[54] X. Li, K. Wei, Y. Aneja, P. Tian, Y. Cui, Matheuristics for the single-path design-balanced service network design problem, Computers & Operations Research 77 (2017) 141 – 153. [55] F. Glover, Future paths for integer programming and links to artificial intelligence, Computers & Operations Research 13 (5) (1986) 533–549. [56] P. Z. Tan, B. Y. Kara, A hub covering model for cargo delivery systems, Networks 49 (1) (2007) 28–39.

CR IP T

[57] D. C. Montgomery, Design and analysis of experiments, 8th Edition, John Wiley & Sons, 2013.

CE

PT

ED

M

AN US

[58] R: A Language and Environment for Statistical Computing, http://www.r-project.org/.

AC

540

34