- Email: [email protected]
Hub interdiction problem variants: Models and metaheuristic solution algorithms
Accepted Manuscript
Hub interdiction problem variants: Models and metaheuristic solution algorithms Nader Ghaffarinasab, Alireza Motallebzadeh PII: DOI: Reference:
S0377-2217(17)31069-X 10.1016/j.ejor.2017.11.058 EOR 14846
To appear in:
European Journal of Operational Research
Received date: Revised date: Accepted date:
2 September 2016 23 November 2017 23 November 2017
Please cite this article as: Nader Ghaffarinasab, Alireza Motallebzadeh, Hub interdiction problem variants: Models and metaheuristic solution algorithms, European Journal of Operational Research (2017), doi: 10.1016/j.ejor.2017.11.058
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 • New hub interdiction problems are introduced to the hub location literature. • New bilevel and single level mathematical formulations are proposed for the problems. • Different multi-stage CA constraints are presented and analyzed.
AC
CE
PT
ED
M
AN US
• Extensive computational study are carried out on the problems.
CR IP T
• Efficient simulated annealing based heuristics are proposed to solve the problems.
1
ACCEPTED MANUSCRIPT
Hub interdiction problem variants: Models and metaheuristic solution algorithms Nader Ghaffarinasaba,∗, Alireza Motallebzadehb a Department b Faculty
of Industrial Engineering, University of Tabriz, Tabriz, Iran of Economics, Management and Business, University of Tabriz, Tabriz, Iran
CR IP T
Abstract Hub location problem (HLP) is one of the strategic problems encountered in designing transportation and telecommunication networks. Regardless of the considered objective in design of hub networks, such as cost minimization or service level maximization, the located hubs can be subject to natural or intentional disruptions after installation. In this paper, we address the multiple allocation p-hub median, p-hub maximal covering, and p-hub center problems under intentional disruptions. In each case, the problem is considered as a Stackelberg game where the leader locates
AN US
p hubs to optimize his/her objective function, whereas the follower tries to identify and interdict r hubs that their loss would diminish the network performance the most. Bilevel and single level mathematical formulations are presented to model the problem from the leader’s and the follower’s perspectives. Furthermore, efficient Simulated Annealing (SA) heuristics are proposed for solving the problems. Extensive computational experiments show the capability of the proposed SA algorithms to obtain the optimal solutions in short computational times. Some managerial insights are also derived based on the obtained numerical results.
M
Keywords: Location, Hub interdiction problems, Competitive models, Bilevel programming, Simulated annealing.
ED
1. Introduction
Hub networks play a major role in reducing cost and enhancing service level in many transportation, telecommunications and computer networks. These networks provide efficient service between many origins-destination (O/D)
PT
pairs of nodes via a set of hubs that serve as switching and flow consolidation points. The hub location problem (HLP) deals with locating hub facilities and allocating non-hub nodes (spokes) to the installed hubs in order to route
CE
the traffic between O/D pairs [4]. Regarding the way the non-hub nodes are allocated to the hubs, there are two main types of allocation in the hub networks called single and multiple allocation. In a single allocation network, all the incoming and outgoing traffic to/from every demand node is routed through a single hub, whereas in a multiple
AC
allocation network, which is the underlying network in this paper, each demand node can receive and send flow through more than one hub. Most of the papers published in the literature of the HLP focus on location and network design decisions in a deterministic environment where the reliability and robustness issues in cases of service disruption are not addressed. However, in real world applications, hub facilities can be lost or stop operating for a variety of reasons, such as natural disasters, inclement weather, equipment malfunctioning, labor strikes, and intentional attacks [42]. In other ∗ Corresponding
author Email addresses: [email protected], [email protected] (Nader Ghaffarinasab ), [email protected] (Alireza Motallebzadeh)
Preprint submitted to European Journal of Operational Research
December 1, 2017
ACCEPTED MANUSCRIPT
words, disasters whether intentional or accidental, can sometimes render a system inoperable or inefficient. In this paper, we consider the multiple allocation p-hub median, p-hub maximal covering, and p-hub center problems under intentional disruptions. The problems are tackled as leader-follower Stackelberg games where the decision maker who makes the initial location decisions is called leader and the interdictor is called follower. Accordingly, for each of the three HLP variants, two models are developed: a single level mixed integer programming (MIP) model on behalf of the follower and a bilevel model on behalf of the leader. To the best of our knowledge, the hub interdiction problem under covering and center objectives are modeled and solved for the first time in the literature in this paper.
CR IP T
We also introduce some multi-stage closest assignment (CA) constraints and compare their performance with those of other existing constraints in terms of solution time by a commercial solver. Furthermore, to solve the above mentioned problems, efficient solution algorithms based on simulated annealing (SA) are proposed. In case of the follower’s problem, the single level models are also solved using a standard optimization package that enables us to evaluate the performance of the proposed SA algorithms as well as the proposed mathematical models.
The remainder of this paper is organized as follows. The next section reviews the relevant literature to the problem.
AN US
In Section 3, we will present the mathematical formulations for the problems. The proposed SA algorithms are presented in Section 4. Computational experiments and corresponding results are presented in Section 5. Finally, some concluding remarks and directions for future works are given in Section 6.
2. Literature Review
M
2.1. Hub location problem
Hub location models can be categorized according to the type of objective that is considered. The most frequently addressed models are minisum hub location models, such as the p-hub median problem and the uncapacitated hub
ED
location problem (with fixed costs), that focus on minimizing the overall setup and operational costs. O’Kelly [50] developed the first quadratic mathematical formulation for uncapacitated single allocation p-hub median problem
PT
(USApHMP). In case of uncapacitated multiple allocation p-hub median problem (UMApHMP), the first mathematical model was introduced by Campbell [10]. Linear integer programming formulations for single and multiple allocation p-hub median problem were presented in [11]. Skorin-Kapov et al. [60] proposed new models for unca-
CE
pacitated single and multiple allocation hub location problems having tight linear programming relaxations. Ernst and Krishnamoorthy [27] introduced an efficient flow based mathematical formulation for the single allocation p-hub
AC
median problem that was later extended to model the multiple allocation p-hub median problem [28]. Kara proved that the p-hub median problem belongs to the class of NP-hard problems [36]. Another important variant of the HLP is the p-hub maximal covering problem which aims at locating a fixed number of hub facilities in such a way that the total covered flow associated with the O/D pairs is maximized. The p-hub maximal covering problem is closely related to the hub set covering problem in which the objective is to cover all the O/D pairs using a minimum number of opened hubs. Campbell [11] formulated the p-hub maximal covering and the hub set covering problems as linear MIPs. Kara and Tansel [37] proposed new models for the single allocation hub set covering problem and proved that it is NP-hard. Later, Wagner [64] introduced novel formulations for both single and multiple allocation hub set covering problems. Peker and Kara [53] addressed the p-hub maximal covering
3
ACCEPTED MANUSCRIPT
problem considering partial coverage, an extension of the traditional binary coverage. The third important variant of the HLP is the p-hub center problem which is defined based on the minimax criterion in which the maximum service cost (distance, time, etc.) between all the O/D pairs is minimized. Campbell [11] was first to introduce the p-hub center problem in the hub location literature. He presented formulations for both single and multiple allocation p-hub center problems. Kara and Tansel [38] addressed different linear formulations for the single allocation p-hub center problem. Later, Ernst et al. [26] developed new mathematical formulations for single and multiple allocation p-hub center problems and showed that their single allocation model is superior to the
CR IP T
one in [38]. They also proved that the problem is NP-hard. Yaman and Elloumi [66] introduced new models for star p-hub center problem to minimize the length of the longest O/D pairs in a star hub network. Contreras et al. [20] presented a generalized model for the classical p-center problem that combines facility location and network design problems by deciding simultaneously on the location of facilities and the design of its underlying network. Although integer programming optimization approaches are used to solve small hub problems, larger instances of the HLP are usually solved by heuristic or metaheuristic procedures. Campbell [12] presented a greedy-interchange
AN US
heuristic for the UMApHMP. Ernst and Krishnamoorthy [28] proposed a heuristic and an exact solution algorithm based on shortest paths for the same problem. Ebery et al. [25] developed an efficient shortest path based heuristic for the capacitated multiple allocation p-hub median problem. Jabalameli et al. [35] proposed an SA heuristic for the single allocation p-hub maximal covering problem. Hamacher and Meyer [33] proposed a binary search algorithm for solving the single allocation p-hub center problem which exploits the relationship between p-hub center and hub covering problems. An adaptive large neighborhood search (ALNS) algorithm was proposed by Chaharsooghi et al.
M
[15] for a reliable multiple allocation HLP with hub failure possibility. Silva and Cunha [59] presented an efficient tabu search (TS) algorithm for solving the uncapacitated single allocation p-hub maximal covering problem. More
ED
recently, Ghaffarinasab et al. [32] developed efficient simulated annealing based heuristics for solving the competitive single and multiple allocation hub location problems. For more details on the HLP and recent advances in this field,
PT
the interested readers are referred to surveys [4, 13, 30], and [19]. 2.2. Facility location with disruption and interdiction In the most of the classical facility location models, it is assumed that systems will operate as designed during its
CE
whole life. However, the network elements in real world environments are subject to various disruptions which can substantially affect the system’s performance. Disruptions, either caused by natural events or intentional attacks,
AC
may lead to unavailability of some elements in the service networks. Drezner [24] was the first researcher who studied location problems with unreliable facilities in p-median problems. Later, Drezner’s problem was tackled with different heuristic method by Lee [40]. Berman et al. [6] studied the location of unreliable facilities on a network where a customer may not be able to receive satisfactory service depending on the distance from the assigned facility. Synder and Daskin [61] considered the reliable p-median problem and the uncapacitated location problem assuming that the facilities are subject to random disruption with equal probabilities. Snyder et al. [62] conducted a review of different reliability considerations in facility location models. Berman et al. [7] introduced a reliable model based on the different failure probabilities for various facilities. In another work, Berman et al. [8] considered location of unreliable facilities where a customer have incomplete information about the state of operating facilities. Cui et
4
ACCEPTED MANUSCRIPT
al. [21] and Li and Ouyang [44] developed continuum approximation models to analyze the reliable facility location problems. Lim et al. [46] considered design of a reliable supply chain network in the presence of random facility disruptions with the option of hardening selected facilities. Peng et al. [54] addressed another reliable networks design problem considering a robust optimization measure in which the relative regret in each scenario is required to be no more than a constant value. Hong et al. [34] proposed two robust mathematical models to solve the emergency response facility location and transportation problems which are responsible for distributing relief goods after disaster. Lei and Tong [43] extended the reliable facility location model developed in [21] to introduce a level
CR IP T
assignment model for the p-median problem under non-uniform failure probabilities. Berman et al. [9] studied a class of location models with unreliable facilities and correlated failures with median and center objectives. More recently, Akg¨ un et al. [1] addressed a facility location problem in the context of disaster management in the pre-disaster phase of operations. They developed an optimization model that minimizes the risk that a demand point may be exposed to because it is not supported by the located facilities.
Unlike the above studies which focus on reliability under random disruption, there are a number of studies addressing
AN US
deliberate disruptions. Deliberate and rational attacks to the installed facilities was first addressed through the use of interdiction models for military applications [58], but later they were also used in business environments such as facility location and transportation network design. Church et al. [18] introduced the r-interdiction median problem (RIMP) and the r-interdiction covering problem (RICP) which were used to identify the set of facilities that, if lost, would harm service delivery the most. Later, Church et al. [57] considered the r-interdiction median problem with fortification (RIMPF) where the objective was to identify the most effective way of protecting a subset of existing
M
facilities so that the impact of the most disruptive attack to r unprotected facilities is minimized. Lei [41] introduced the r-interdiction center (RICent) problem in which the objective is to identify and interdict a
ED
fixed number of installed facilities in such a manner that the worst-case performance of the system (largest distance between demand points and the serving facilities) is further diminished. Aksen et al. [2] considered a budgetconstrained version of the RIMPF in which the facilities are assumed to have different protection costs and flexible
PT
capacities which can be expanded to accommodate the demand of customers previously assigned to interdicted facilities. Church and Scaparra [16] introduced a probabilistic version of RIMP where an attempted interdiction is
CE
successful only with a given probability. Losada et al. [47] further extended the model in [16] by assuming that the probability of a successful disrupting attack depends on its intensity and on the amount of disrupting resources used in the attack. Liberatore et al. [45] addressed a stochastic variant of RIMPF with uncertain number of interdictions.
AC
O’Hanley and Church [49] proposed a new version of RICP that maximizes a combination of initial coverage by p facilities and the minimum coverage level following the loss of the most critical r facilities. A new variation of the RIMP considering capacitated facilities and outsourcing option with partial interdiction was introduced by Aksen et al. [3] in which an attacked facility may preserve part of its capacity and the unmet demand after interdiction can be outsourced at some cost. According to the literature, only three papers have been published so far on hub location under intentional disruptions. Lei [42] introduced the hub interdiction median (HIM) problem which can be used to identify the set of critical facilities in a hub network. The author also proposed two bilevel formulations called the hub median protection problem (HMPP) and the HMPP with fortification (HMPP-F) where the former aims at locating hubs in anticipation 5
ACCEPTED MANUSCRIPT
of a worst-case facility loss and the latter tries to fortify a subset of existing hubs in a network to minimize postinterdiction system cost. Parvaresh et al. [51] formulated the multiple allocation p-hub median problem under intentional disruptions as a bilevel model in which the follower’s objective is to identify and interdict those hubs which would deteriorate the system’s efficiency in the worst possible way. For solving the problem, they proposed two algorithms based on SA and solved instances of up to 60 nodes using their heuristics. In a follow up work, Parvaresh et al. [52] extended their former model with two objective functions for the upper level problem where the leader aims at locating its hubs to minimize normal and worst-case transportation cost. They developed two
CR IP T
multi-objective metaheuristic algorithms based on SA and Tabu Search (TS) to solve this problem with instances of up to 50 nodes.
3. Mathematical Formulation
In this section, single level and bilevel mathematical models are presented for the three variants of the hub interdiction problem. In our proposed models, it is assumed that both the leader and the follower have complete
AN US
information about the game and will act rationally. Our models are based on the uncapacitated multiple allocation hub location problem under the assumption that the installed hubs are fully interconnected. Furthermore, we assume that transportation cost on the inter-hub connections are discounted by a fixed discount factor and that the direct connections between non-hub nodes are not allowed. The following notations are used hereafter: Sets and parameters
set of all demand nodes (i, j, k, m ∈ N )
H
set of nodes that are available for locating hubs H ⊆ N
wij
amount of flow originated at node i and destined to node j
dij
transportation cost (distance, time, etc.) of a unit flow from node i to node j
α
transfer cost per unit flow and per unit distance on inter-hub connections (0 ≤ α ≤ 1)
dijkm
cost of routing one unit of flow wij via hubs k and m, calculated as dijkm = dik + αdkm + dmj
p
number of hubs to be located by the leader
r
number of hubs to be interdicted by the follower
ED
PT
CE
β
M
N
covering radius
AC
Decision variables Xijkm Zk
Yk
fraction 1, = 0, 1, = 0,
of flow originated at node i and destined to node j routed via hubs k and m if node k is selected as a hub, otherwise if the hub located at node k is interdicted, otherwise
6
ACCEPTED MANUSCRIPT
Vij
=
Uij
=
1, if flow originated at node i and destined to node j is covered after interdiction, 0, otherwise 1, if the maximum O/D cost in the network is associated with node pair (i, j), 0, otherwise
CR IP T
3.1. r-hub interdiction median problem Assume that the leader is already operating in a market with hubs located at set Z ⊆ H of p nodes. The follower then attempts to attack r of the leader’s hubs (r ≤ p) in such a way that the leader’s operational cost is maximized. Define Fijkm as the set of existing hub pairs (q, s) for which the cost of routing one unit of flow between nodes i and j via hubs q and s is higher than the cost of routing the same flow via hubs k and m. In other words:
∀i, j ∈ N , k, m ∈ Z
AN US
Fijkm = {(q, s) ∈ Z × Z|diq + αdqs + dsj > dik + αdkm + dmj }
The follower’s problem then consists of the selection of a set of r installed hubs which will be disrupted so that total cost incurred by the leader is maximized. The MIP model for the multiple allocation r-hub interdiction median problem (r-HIMP) as an improved variant of the original HIM formulation presented by Lei [41, 42] can now be written as: XXX X
wij dijkm Xijkm
i∈N j∈N k∈Z m∈Z
s.t.:
X
Yk = r
X X
(1) (2)
ED
k∈Z
M
max
∀i, j ∈ N
(3)
∀i, j ∈ N , k ∈ Z
(4)
∀i, j ∈ N , k, m ∈ Z
(5)
Xijkm ≥ 0
∀i, j ∈ N , k, m ∈ Z
(6)
Yk ∈ {0, 1}
∀k ∈ Z
(7)
Xijkm = 1
k∈Z m∈Z
Xijkm +
X
PT
X
m∈Z
X
Xijmk ≤ 1 − Yk
Xijqs ≤ Yk + Ym
AC
CE
(q,s)∈Fijkm
m∈Z|m6=k
The objective function (1) maximizes the total transportation cost of the leader. Constraint (2) determines the number of hubs to be interdicted by the follower. Constraints (3) state that for each pair of nodes i and j, the whole flow wij must be routed via some hub pair. Constraints (4) ensure that flows can be routed via an installed hub if that hub is not interdicted. Constraints (5) are called the (multi-stage) closest assignment (CA) constraints which are generalized from the single-stage CA constraints originally introduced by Wagner and Falkson (W-F) [65]. They ensure that each O/D flow must be routed via a least-cost pair of hubs that are not disrupted. To be more precise, these CA constraints imply that for a pair of hubs (k, m), if both k and m are kept open (i.e., Yk + Ym = 0), the O/D flow wij cannot be routed via a pair of hubs which are more expensive than routing via k and m. Finally, constraints 7
ACCEPTED MANUSCRIPT
(6) and (7) are non-negativity and binary constraints, respectively. Now, let us assume that the leader is aware of the subsequent attack by the follower, and thus tries to locate its p hubs in such a way that the total transportation cost after the upcoming attack by the follower (including r interdictions) be minimum. The bilevel model for the hub interdiction median problem (BHIMP) can be written as: min H1 (Z) X s.t.: Zk = p
(8) (9)
k∈H
∀k ∈ H
where H1 (Z) = max
XXX X
wij dijkm Xijkm
i∈N j∈N k∈Z m∈Z
Yk ≤ Zk
AN US
s.t.: (2) − (7)
CR IP T
Zk ∈ {0, 1}
∀k ∈ Z
(10)
(11)
(12)
The objective function (8) minimizes the highest possible level of transportation cost after disruption of r hubs which is calculated by (11) as the objective of the lower level follower’s problem. Constraint (9) enforces the number of installed hubs by the leader to be p. Constraints (12) state that the follower can only interdict the hubs that have already been opened by the leader.
M
3.2. r-hub interdiction maximal covering problem
In the maximal covering variant of the problem, it is assumed that the leader’s objective is to maximize the total
ED
covered flows in the network. The O/D flow wij is said to be covered if it is routed via a pair of hubs in such a way that corresponding unit routing cost is no more than a prespecified covering radius (β). The follower aims to hit r of the installed hubs in such a manner that the total covered flows by the leader after interdiction is minimized. To
PT
model the problem, for each node pair (i, j), we define Mij as the set of hub pairs (k, m) through which the cost of routing one unit of flow wij is less than β. In other words:
CE
Mij = {(k, m) ∈ Z × Z|dik + αdkm + dmj ≤ β}
∀i, j ∈ N
Based on the above discussion, the multiple allocation r-hub interdiction maximal covering problem (r-HIMCP)
AC
solves the following MIP model: min
XX
wij Vij
(13)
i∈N j∈N
s.t.:(2), (7) Vij ≥ 1 − Yk − Ym
∀i, j ∈ N , (k, m) ∈ Mij
(14)
Vij ∈ {0, 1}
∀i, j ∈ N
(15)
The objective function (13) minimizes the total flows covered in the network. Constraints (14) stipulate that an O/D flow wij must be covered unless all the hub pairs (k, m) that currently cover it (i.e., (k, m) ∈ Mij ) are interdicted. 8
ACCEPTED MANUSCRIPT
We should note that the maximal covering version of the hub interdiction problem implicitly enforces CA conditions. In other words, in any optimal solution of this problem, all O/D flows are routed through a minimum cost path. Hence, there is no need to explicitly add CA constraints to this formulation. The bilevel model for the hub interdiction maximal covering problem (BHIMCP) can be proposed as: max H2 (Z)
(16)
where H2 (Z) = min
XX
wij Vij
i∈N j∈N
s.t.: (2), (7), (12), (14), (15)
CR IP T
s.t.:(9), (10)
(17)
The objective function (16) maximizes the lowest possible level of total covered flows after disruption of r hubs which
AN US
is calculated by (17) as the objective of the lower level follower’s problem. 3.3. r-hub interdiction center problem
Hub center problems are important for the delivery of perishable or time sensitive items. While in the median and maximal covering versions of the problem, the average-case performance of the system is diminished by the attacker, in the hub center interdiction problem, the follower plans to worsen the worst-case performance. In other
M
words, assuming a set of p existing hubs for the leader, the follower aims to select and disrupt r hubs so that highest transportation between all O/D pairs is maximized. The MIP model for the multiple allocation r-hub interdiction
ED
center problem (r-HICP) can then be written as:
PT
max
max
(i,j)∈N ×N
{
X X
k∈Z m∈Z
dijkm Xijkm }
(18)
s.t.:(2) − (7)
The objective function (18) maximizes the maximum transportation cost between all O/D pairs. Note that the above
CE
model is a nonlinear problem of maximax type which can be linearized by defining new set of variables. Let γ be a nonnegative variable acting as a surrogate for the maximum of O/D costs in the network. The MIP model for the
AC
multiple allocation r-hub interdiction center problem (r-HICP) can now be written as: max γ
(19)
s.t.:(2) − (7) X X γ≤ dijkm Xijkm + (1 − Uij )M k∈Z m∈Z
XX
∀i, j ∈ N
Uij = 1
(20) (21)
i∈N j∈N
γ≥0
(22)
Uij ∈ {0, 1}
∀i, j ∈ N 9
(23)
ACCEPTED MANUSCRIPT
where M is a sufficiently large number. An appropriate value for this number which results in a tight formulation is M = (2 + α)dmax where dmax = maxi,j∈N {dij }. The bilevel model for the hub interdiction center problem (BHICP) can be written as: min H3 (Z)
(24)
s.t.:(9), (10)
H3 (Z) = max γ s.t.: (2) − (7), (12), (20) − (23)
CR IP T
where (25)
The objective function (24) minimizes the highest transportation cost between all O/D pairs after disruption of r hubs which is calculated by (25) as the objective of the lower level follower’s problem.
Bilevel models are hard to solve even for a small number of decision variables [5], [22]. Therefore, to solve the
AN US
proposed models, we have developed efficient metaheuristic solution algorithms that are going to be discussed in detail in Section 4. 3.4. Alternative CA constraints
As stated earlier, CA constraints ensure that each O/D flow is routed through a least-cost pair of hubs that are not
M
disrupted. So far, different types of (single-stage) CA constraints have been introduced for discrete facility location problems (see [31, 41, 29] for comparison between various CA constraints). These single-stage CA constraints can be generalized as multi-stage CA constraints to be used in HLP models. For instance, Lei [41, 42] extended the
ED
single-stage CA constraints formulated by Church and Cohon (C-C) [17] and used it in his hub interdiction median (HIM) model as follows:
Xijqs + Xijkm ≥ 1 − Yk − Ym
PT
X
(q,s)∈Cijkm
∀i, j ∈ N , k, m ∈ Z
(26)
CE
where for all i, j ∈ N and k, m ∈ Z, the set Cijkm is defined as follows: Cijkm = {(q, s) ∈ Z × Z|diq + αdqs + dsj < dik + αdkm + dmj or
AC
diq + αdqs + dsj = dik + αdkm + dmj and (q < k or q = k and s < m)}
Constraints (26) maintain that if both hubs k and m are not disrupted (i.e., 1 − Yk − Ym = 1) and if no route has a P lower cost than route i-k-m-j (i.e., (q,s)∈Cijkm Xijqs = 0 ), then the route i-k-m-j must be assigned or chosen (i.e., Xijkm = 1).
Another type of CA constraints is proposed by Rojeski and ReVelle (R-R) [55] which has been extended to its multi-stage counterpart by Lei [41] as: X
(q,s)∈Cijkm
Qqs + Xijkm ≥ 1 − Yk − Ym
10
∀i, j ∈ N , k, m ∈ Z
(27)
ACCEPTED MANUSCRIPT
where Qkm = 1 if both hubs k and m are kept open. In other words, this variable can be defined as the product Qkm = (1 − Yk )(1 − Ym ). However, in order to linearize this equality, the following constraints are added along with (27): ∀k, m ∈ Z
(28)
Qkm ≤ 1 − Yk
∀k, m ∈ Z
(29)
Qkm ≤ 1 − Ym
∀k, m ∈ Z
(30)
Qkm ≥ 0
∀k, m ∈ Z
(31)
CR IP T
Qkm ≥ 1 − Yk − Ym
The logic behind the R-R constraints is as follows: If none of the hubs k and m are interdicted (i.e., 1 − Yk − Ym = 1), and there is no other open hub pair providing a lower routing cost than that of the hub pair (k, m) to the O/D flow P wij (i.e., (q,s)∈Cijkm Qqs = 0), then the flow wij must be routed via the hubs k and m (i.e., Xijkm = 1). As a further example, we extend the single-stage CA constraints originally introduced by Dobson and Karmarkar
AN US
(D-K) [23] to its multi-stage equivalent as: Qqs + Xijkm ≤ 1
∀i, j ∈ N , k, m ∈ Z, q, s ∈ Cijkm
(32)
along with the constraints (28)-(31). The idea of the D-K constraints is that the O/D flow wij is not routed via the hubs k and m (i.e., Xijkm = 0) if there is another pair of open hubs which provide a shorter path than that of the hubs k and m (i.e., 1 − Qqs = 0).
M
In Section 5, we will conduct a set of computational experiments in order to compare the performance of the four presented CA constraints, namely the W-F, C-C, R-R and D-K constraints in terms of solution time by using a
ED
commercial solver.
4. Metaheuristic Solution Algorithm
PT
In this section, we describe in detail the proposed simulated annealing (SA) metaheuristic algorithm for solving the six problems presented in Section 3. Inspired by heat treatment of metals in metallurgical industries, SA is
CE
a metaheuristic optimization algorithm which is effective in solving combinational optimization problems. SA was ˇ developed in 1953 by Metropolis et al. [48] and was independently described by Kirkpatrick et al. [39] and Cern` y [14]. The algorithm adopts the general concepts in gradual cooling process of metals for solving optimization problems.
AC
To this end, the SA algorithm starts with an initial solution at a high temperature and consecutively moves to the neighboring solutions via algorithm loops letting the temperature be reduced through the algorithm iterations. At each step, if the newly generated neighboring solution is better, in terms of the objective value, than the current solution, the current solution is replaced by the new solution. Otherwise, the algorithm accepts the new solution with a probability e−|∆E|/T where ∆E is the difference between the objective function values of the current and the new solutions and T is the current temperature. In early stages where the temperature is too high, there is a high probability to accept poor solutions. However, with a gradual decrease in temperature, there will be less probability to accept a bad solution in final stages. The algorithm stops if either the final temperature limit is reached or there is no improvements in a number of consecutive algorithm iterations. 11
ACCEPTED MANUSCRIPT
In the remainder of this section, we discuss several features of our proposed SA algorithms such as solution representation scheme, initial solution generation procedure, neighbor generation mechanisms, parameters used in the algorithms, and pseudo-codes for the proposed algorithms. 4.1. Solution Representation One of the important components of any metaheuristic algorithm is the solution representation. It has to be designed in such a way that facilitates neighbor generation and fast calculation of the objective functions. It must also
CR IP T
guarantee accessibility to the entire solution space [56]. In this paper, solutions are represented by one-dimensional arrays that contain numbers associated with the opened and operational hubs. For the leader’s problem, the array’s size is p comprising the set of hubs Z ⊆ H opened by the leader, whereas for the follower’s problem the representation array is of size r containing a subset Y of opened hubs that are interdicted by the follower (Y ⊆ Z). Note that having known the set of opened and operational (not interdicted) hubs, for each O/D pair (i, j), one can easily determine the paths for routing the associated flow wij by solving a shortest path problem or using enumeration.
AN US
4.2. Initial Solution Generation
The initial solution can affect the performance of a metaheuristic algorithm. Indeed, the quality of the final solution obtained by the metaheuristic algorithms relies to some extent on the quality of the initial solution. To produce an initial solution for the leader’s problem, we randomly select a subset Z ⊆ H of size p to serve as hub nodes. In case of the follower’s problem, an initial solution is generated by randomly selecting a subset Y ⊆ Z
M
showing the r hubs that are disrupted by the follower. The routing decisions are then determined based on the set
4.3. Neighborhood Structures
ED
of opened and operational hubs.
We define and apply a single operator for generating neighboring solutions in our proposed SA algorithms. This operator is called “Swap” and is used to alter one of the hubs in the solution representation array. In case of the
PT
leader’s problem, we first select a hub node and a non-hub node, randomly. Afterwards, the selected hub node becomes non-hub and the selected non-hub node becomes hub. For the follower’s problem, on the other hand, a
CE
disrupted hub and a non-disrupted hub are selected randomly and are changed afterwards. 4.4. Parameters used in the SA procedure
AC
Five parameters, namely T0 , Tf , δ1 , δ2 , and Imax are used in the proposed SA algorithms. T0 denotes the initial
temperature, whereas Tf is the final temperature below which the SA procedure is terminated. δ1 and δ2 are used as cooling rate parameters in the follower’s and the leader’s problems, respectively. Finally, Imax is the number of generated solutions at each temperature of the algorithm. 4.5. The overall SA Algorithms To solve the follower’s problem, it is assumed that the leader has already located his/her hubs on a set of candidate nodes Z ⊆ H. We start our algorithms with generating an initial solution Y ⊆ Z, as described previously. Then the follower’s objective value, f (Z, Y), can be calculated based on either (1), (13), or (18), depending on the variant of 12
ACCEPTED MANUSCRIPT
the HLP under consideration. We set the initial temperature as T0 which is an input parameter of the algorithm. Ybest denotes the best solution found so far and fbest denotes the corresponding objective function value. At each
temperature, we generate a new solution Y 0 by using the current solution based on the “Swap” operator presented
above. Subsequently, the objective function of the new solution f (Z, Y 0 ) is calculated. We define ∆E as the difference between the objectives of the new and current solutions, i.e., ∆E = f (Z, Y 0 ) − f (Z, Y). If ∆E > 0 for the r-HIMP
and r-HICP, or ∆E < 0 for the r-HIMCP, we update the current solution as Y ← Y 0 and if the objective f (Z, Y 0 )
of the new solution Y 0 is even better than best objective fbest so far, we also set fbest ← f (Z, Y 0 ) and Ybest ← Y 0 .
CR IP T
Otherwise, if ∆E ≤ 0 for the r-HIMP and r-HICP, or ∆E ≥ 0 for the r-HIMCP, we generate a random number ρ
from the interval [0,1]. If ρ is larger than e−|∆E|/T , we update the current solution as Y ← Y 0 . In other words, we
accept the solutions of worse quality with probability e−|∆E|/T to help the algorithm not get trapped in local optima. We repeat this cycle for Imax times at each temperature with the best solution obtained so far. Subsequently, we reduce the temperature at each iteration as T ← δ × T . The algorithm is terminated when the current temperature T drops below the prespecified final temperature Tf .
AN US
The SA algorithm for the leader’s problem is more complex than that of the follower’s problem as we aim at optimizing a bilevel model. In this case, we have a nested algorithm comprising an upper-level SA for optimizing the leader’s decisions and a lower-level SA that optimizes the follower’s decisions whenever a new solution is generated for the leader by the upper-level SA. We start the algorithm by generating an initial solution for the leader Z and improve this solution through the subsequent algorithm iterations. At each iteration of the upper-level SA where a new solution Z 0 is generated for the leader, the lower-level SA is run to obtain the corresponding best solution for the
M
follower Y ∗ (Z 0 ). In this case we define ∆E = f (Z 0 , Y ∗ (Z 0 )) − f (Z, Y ∗ (Z)). The new solution Z 0 is accepted if ∆E < 0 for the BHIMP and BHICP, or ∆E > 0 for the BHIMCP. Furthermore, if the objective f (Z 0 , Y ∗ (Z 0 )) of
ED
the new solution is even better than the best objective so far (fbest ), we set fbest ← f (Z 0 , Y ∗ (Z 0 )) and Zbest ← Z 0 .
The pseudo-codes of the proposed SA algorithms for the follower’s and the leader’s problems are illustrated in
PT
Algorithms 1 and 2, respectively.
5. Computational Experiments
CE
In this section, we present the results obtained from our computational experiments. Extensive experimentations are conducted in order to test the efficiency of the proposed mathematical models and the CA constraints as well as the SA algorithms. To this end, we use two famous data sets from the literature of the HLP: the CAB and the
AC
TR data sets. The CAB data set introduced by O’Kelly [50] is based on the airline passenger interactions between 25 US cities in 1970 evaluated by the Civil Aeronautics Board (CAB). This data set has been used by most of the hub location researchers in the literature in which all the 25 nodes are candidates for being hubs (|H| = 25). The second data set that is used in our computational experiments is the TR data set [63] which is based on the cargo flows between 81 cities of Turkey where only 22 of these cities are candidate nodes for locating hubs (|H| = 22). In P both the data sets, the O/D flows are scaled by the sum of flows (i.e., i,j wij ). The proposed SA algorithms are
implemented in Microsoft Visual C# 2013 (version 5.0). Moreover, the mathematical models r-HIMP, r-HIMCP, and r-HICP are solved independently using CPLEX version 12.6. All the experiments have been run on a computer
13
ACCEPTED MANUSCRIPT
Algorithm 1 : SA for the follower’s problem (T0 , Tf , δ1 , Imax , α, Z, r) 1: Generate a random initial solution Y 2: Calculate f (Z, Y) 3: T ← T0 , fbest ← f (Z, Y), Ybest ← Y, I ← 0 4: while T > Tf do
7:
Generate a new solution Y 0 based on Y using “Swap” operator Calculate f (Z, Y 0 )
8:
∆E ← f (Z, Y 0 ) − f (Z, Y)
9:
if ∆E > 0 (or ∆E < 0) then Y ← Y 0 , f (Z, Y) ← f (Z, Y 0 )
10: 11:
CR IP T
6:
for I < Imax do
else ρ ← rand(0, 1)
12:
if ρ > e−|∆E|/T then
13:
Y ← Y 0 , f (Z, Y) ← f (Z, Y 0 );
14:
end if
15: 16:
end if
17:
if f (Z, Y) > fbest (or f (Z, Y) < fbest ) then end if
20:
Y ← Ybest
21:
I ←I +1
22:
end for
23:
I←0
24:
T ← δ1 × T
ED
19:
M
Ybest ← Y, fbest ← f (Z, Y)
18:
PT
25: end while 26: return Ybest , fbest
AN US
5:
equipped with Intel(R) Core(TM) i3-3220 CPU of 3.30 GHz and 16 GB of RAM, using the Microsoft Windows 7
CE
operating system.
5.1. Test instances and parameters setting
AC
We have generated and used a large number of instances from the CAB and TR data sets. The number of hubs
to be opened by the leader (p) for the CAB data set is taken as 3, 4, or 5, whereas in the TR data set, p is chosen as 6, 10, or 14. The number of disrupted hubs by the follower are selected as r ∈ {0, 1, ..., (p − 1)} for the CAB instances and as r ∈ {0, 1, ..., 5} for the TR instances. The parameter α is considered at three levels: 0.4, 0.6, and 0.8 for both data sets. For solving the maximal covering variant, different levels for the covering radius parameter (β) are considered. For the CAB data set, β takes a value as 500, 1000, or 1500, whereas for the TR data set the covering radius is chosen as 400, 700, or 1000. The parameters of the proposed SA algorithms are tuned by setting a good trade-off between time and quality of
14
ACCEPTED MANUSCRIPT
Algorithm 2 : SA for the leader’s (bilevel) problem (T0 , Tf , δ2 , α, p, r) 1: Generate a random initial solution Z
2: Get Y ∗ (Z) and f (Z, Y ∗ (Z)) by solving the follower’s problem using Algorithm 1 3: T ← T0 , fbest ← f (Z, Y ∗ (Z)), Zbest ← Z 4: while T > Tf do
Get Y ∗ (Z 0 ) and f (Z 0 , Y ∗ (Z 0 )) by solving the follower’s problem using Algorithm 1
7:
∆E ← f (Z 0 , Y ∗ (Z 0 )) − f (Z, Y ∗ (Z))
8:
if ∆E < 0 (or ∆E > 0) then
9: 10: 11: 12:
Z ← Z 0 , f (Z, Y ∗ (Z)) ← f (Z 0 , Y ∗ (Z 0 )) else ρ ← rand(0, 1)
if ρ > e−|∆E|/T then Z ← Z 0 , f (Z, Y ∗ (Z)) ← f (Z 0 , Y ∗ (Z 0 ))
13: 14:
CR IP T
6:
Generate a new solution Z 0 based on Z using “Swap” operator
end if
AN US
5:
15:
end if
16:
if f (Z, Y ∗ (Z)) < fbest (or f (Z, Y ∗ (Z)) > fbest ) then
17:
Zbest ← Z, fbest ← f (Z, Y ∗ (Z))
18:
end if
19:
T ← δ2 × T
20: end while
M
21: return Zbest , fbest
ED
the final solutions. In a set of preliminary experiments, different combinations of parameters were tested on a large number of test instances and the values reported in Table 1 have been selected as the best values which lead to
PT
high-quality solutions in short CPU times both for the leader’s and the follower’s problems under the two data sets.
AC
CE
Table 1: SA parameters used for different data sets Median
Max. Cov.
Center
CAB
TR
CAB
TR
CAB
TR
T0
100
200
100
200
100
200
Tf
1
1
1
1
1
1
δ1
0.985
0.985
0.995
0.995
0.995
0.999
δ2
0.99
0.99
0.992
0.992
0.992
0.997
Imax
4
4
4
4
4
4
5.2. Comparison of the performance of different CA constraints In the first set of our computational experiments, we try to test the performance of the four different multi-stage CA constraints introduced in Section 3, namely the W-F, C-C, R-R, and D-K constraints. To this end, 12 test instances with different p and r values have been selected from the CAB and TR data sets. For each instance, the r-hub interdiction median problem (r-HIMP) and the r-hub interdiction center problem (r-HICP) are solved using
15
ACCEPTED MANUSCRIPT
different CA constraints under three distinct values of α by CPLEX. The existing p hubs are assumed to be located based on the solutions of the p-hub median and the p-hub center problems, respectively. Since the model proposed for the r-hub interdiction maximal covering problem (r-HIMCP) does not need CA constraints, we have not included this problem in the current part of our experiments. The results obtained by solving the r-HIMP are reported in Table 2. The first column in this table represent the name of the data set, whereas the second and third columns show the number of the leader’s (existing) hubs and the number of hubs which are disrupted by the follower, respectively. The next five columns report the solution results for the case of α = 0.4. More specifically, the column labeled as
CR IP T
“OF” gives the optimal value of the objective function and the four columns under the heading “CPU(s)” show the solution times, in seconds, that are spent to obtain the optimal solutions by CPLEX using different CA constraints. Similarly, the results for the cases of α = 0.6 and 0.8 are reported in the next columns. We used a time limit of 10800 second (3 hours) for CPLEX in our experiments.
Table 2: CPLEX results with different CA constraints for the r-HIMP α = 0.4
5
TR
6
10
OF
α = 0.6
CPU(s)
OF
C-C
R-R
D-K
W-F
C-C
R-R
D-K
OF
CPU(s) R-R
D-K
1
1088.55
0.27
0.50
0.53
0.92
1140.78
0.13
0.30
0.30
0.55
1190.91
0.15
0.22
0.17
0.54
2
1718.25
0.72
0.84
0.42
0.64
1413.48
0.59
0.70
0.72
0.81
1422.42
0.80
0.90
0.42
0.80
3
3517.63
0.14
0.18
0.24
0.61
3517.63
0.14
0.20
0.13
0.59
3517.63
0.13
0.14
0.13
0.55
1
1014.74
0.51
0.43
0.74
1.44
1095.45
0.40
0.42
0.60
1.44
1159.49
0.20
0.21
0.39
1.09
2
1606.23
1.79
2.05
0.85
1.43
1656.07
1.82
1.96
0.80
1.11
1623.38
1.10
1.14
0.53
1.07
3
1903.06
2.15
2.09
1.10
2.08
1946.07
1.68
1.88
1.37
1.53
1930.37
1.89
1.62
0.75
2.70
1
666.11
8.91
15.18
28.32
48.42
734.69
5.77
14.23
33.53
49.62
786.67
5.30
12.06
34.89
53.11
2
860.75
17.83
169.44
72.85
1869.90
867.44
14.76
125.27
88.12
502.91
899.02
10.65
111.13
70.24
565.70
3
1147.94
15.52
345.21
128.52
1656.50
1080.74
23.57
203.95
105.91
890.22
1096.33
13.33
194.80
91.79
500.18
1
538.78
33.11
78.90
278.13
2441.04
636.98
30.79
78.10
303.77
1301.61
710.74
28.89
64.61
328.04
1042.51
2
605.50
557.97
2132.21
530.51
> 10800
675.49
454.60
2340.90
466.20
> 10800
742.41
289.61
1813.45
424.89
> 10800
3
701.95
709.52
2941.45
917.88
> 10800
725.17
524.10
3896.84
814.43
> 10800
779.72
397.02
2628.02
582.44
> 10800
112.37
474.04
163.34
> 2301.92
88.20
555.40
151.32
> 2029.20
62.42
402.36
127.89
> 1980.69
Avg.
W-F
α = 0.8
CPU(s)
AN US
4
r
M
CAB
p
ED
Data set
W-F
C-C
PT
Note from Table 2 that the solution time for the model with W-F constraints is the smallest among the four CA constraints. The average times for solving the instances using the W-F constraints are less than two minutes for the
CE
three values of α. The second best performance among the four CA constraints belongs to the R-R constraints with the average solution times of less than three minutes. The third best performance in terms of solution times belongs to the C-C constraints, whereas the weakest performance is associated with the D-K constraints. Observe that the
AC
model with the D-K constraints was not able to solve the two largest instances within the time limit of three hours. Note also that, for the instances of the CAB data set, the four CA conditions have more or less similar solution times. However, for larger instances, as in the TR data set, the differences between the solution times associated with different CA constraints become more recognizable. Table 3 shows the results obtained by CPLEX for solving the r-HICP with different CA constraints.
16
ACCEPTED MANUSCRIPT
Table 3: CPLEX results with different CA constraints for the r-HICP α = 0.4
CAB
p
4
5
TR
6
10
r
OF
α = 0.6
CPU(s)
OF
C-C
R-R
D-K
OF
W-F
C-C
R-R
D-K
CPU(s) R-R
D-K
1
2736.53
0.47
0.78
0.77
1.10
2752.40
0.26
0.62
0.76
0.92
4915.25
0.25
0.53
0.68
0.81
2
5007.65
1.28
1.19
0.59
1.27
5007.65
0.90
0.92
0.57
1.06
5053.12
0.83
0.77
0.58
0.82
3
5451.58
0.53
0.54
0.86
1.23
5451.58
0.77
0.36
0.34
0.67
5406.80
0.59
0.28
0.32
0.74
1
2390.81
0.58
0.64
1.24
1.49
2583.17
0.33
0.48
0.60
1.34
2824.84
0.29
0.36
0.54
1.41
2
3012.90
1.55
1.99
1.29
1.97
4598.85
1.07
2.13
1.29
1.81
4598.85
0.78
1.78
2.61
2.77
3
5007.65
1.01
8.51
1.62
2.81
5229.61
1.64
7.76
3.76
2.96
5229.61
1.24
3.10
2.01
2.34
1
1522.20
20.03
43.15
35.53
60.13
1728.20
14.19
26.28
33.97
57.25
1876.60
12.39
15.69
34.83
61.12
2
1834.20
205.84
543.74
140.39
10065.01
2282.00
65.10
352.25
86.07
2757.83
2282.00
32.07
493.28
87.77
5674.85
3
3014.00
126.58
1205.43
513.81
1438.85
3336.00
67.39
519.78
639.19
1
1317.00
71.19
127.82
305.90
888.34
1610.40
69.47
119.12
414.17
2
1514.80
1063.41
4731.96
793.51
> 10800
1719.00
854.86
4961.58
1318.31
3
2282.00
1286.67
> 10800
> 10800
> 10800
1856.00
909.92
> 10800
> 10800
231.60
> 1455.48
> 1049.63
> 2838.52
165.49
> 1399.27
> 1108.25
Avg.
W-F
α = 0.8
CPU(s)
W-F
C-C
CR IP T
Data set
1201.07
3032.00
46.80
320.13
112.56
1126.72
1870.60
45.06
182.41
365.71
950.48
> 10800
1915.20
260.55
6410.62
1221.65
> 10800
> 10800
1934.20
1184.03
> 10800
> 10800
> 10800
132.07
> 1519.08
> 1052.44
> 2391.83
> 2229.30
406.57
It can be seen from Table 3 that the solution times for the r-HICP are, on average, larger than the corresponding
AN US
times for the r-HIMP. Furthermore, one can observe that, similar to the case of median problem, the W-F constraints perform better than the other three types of multi-stage CA constraints for the center problem. The R-R, C-C, and D-K constraints are placed in the second, the third, and the fourth ranks in terms of the average solution time. The results also show that, it is possible for CPLEX to solve all the instances within the three hour time limit only by using W-F constraints. For all the other three CA constraints, at least one of the instances are not solved within this time limit. Based on the reported results in Tables 2 and 3, we use the W-F constraints for solving problems by
M
the commercial solver in the rest of our computational experiments. 5.3. Results for hub interdiction median problems
ED
Results obtained by solving the r-HIMP on the CAB data set are presented in Table 4. In case of the follower’s problem, it is assumed that the leader has already located its hubs on a subset Z ⊆ H of nodes. To this end, we
PT
assume that the leader is unaware of the upcoming attack by the follower, and hence selects its p hubs based on the solution of the uncapacitated multiple allocation p-hub median problem (UMApHMP) aiming to minimize its
AC
CE
transportation cost.
17
ACCEPTED MANUSCRIPT
Table 4: Results for the r-HIMP with the CAB data set α = 0.4 p
r
CPLEX
α = 0.6 SA
CPLEX
α = 0.8 SA
CPLEX
SA
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
0
859.63
0.21
859.63
0.02
949.23
0.11
949.23
0.02
1020.03
0.09
1020.03
0.02
1
1206.33
0.23
1206.33
0.02
1263.61
0.12
1263.61
0.02
1309.29
0.12
1309.29
0.02
2
3517.63
0.08
3517.63
0.01
3517.63
0.07
3517.63
0.01
3517.63
0.06
3517.63
0.01
0
754.48
0.27
754.48
0.01
866.44
0.16
866.44
0.02
951.75
0.07
951.75
0.02
1
1088.55
0.27
1088.55
0.02
1140.78
0.13
1140.78
0.02
1190.91
0.15
1190.91
0.02
2
1718.25
0.72
1718.25
0.02
1413.48
0.59
1413.48
0.02
1422.42
0.80
1422.42
0.02
3
3517.63
0.14
3517.63
0.01
3517.63
0.14
3517.63
0.01
3517.63
0.13
3517.63
0.01
0
676.33
0.37
676.33
0.02
804.70
0.28
804.70
0.02
910.35
0.13
910.35
0.02
1
1014.74
0.51
1014.74
0.03
1095.45
0.40
1095.45
0.04
1159.49
0.20
1159.49
0.04
2
1606.23
1.79
1606.23
0.02
1656.07
1.82
1656.07
0.02
1623.38
1.10
1623.38
0.02
3
1903.06
2.15
1903.06
0.02
1946.07
1.68
1946.07
0.02
1930.37
1.89
1930.37
0.02
4
3517.63
0.39
3517.63
0.02
3517.63
0.24
3517.63
0.02
3517.63
0.23
3517.63
0.02
Avg.
1781.71
0.59
1781.71
0.01
1807.39
0.47
1807.39
0.02
1839.24
0.41
1839.24
0.02
3
4
5
CR IP T
OF
AN US
Note that for each value of p, when the value of r equals 0, the leader’s cost is the same as the optimal cost of UMApHMP, as the follower has not interdicted any hubs. Observe also that, as the value of r gets larger, the leader’s cost increases substantially. This can be interpreted as a result of the fact that the leader has decided on the location of its hubs without being aware of upcoming disruptions. Furthermore, the results show that higher values for the discount factor (α) result in slightly elevated costs for the leader as the inter-hub transportation cost
M
grow with the value of the discount factor. It can be seen from the above table that the proposed SA algorithm has obtained the optimal solution for all instances within a negligible CPU time, indicating the high efficiency of the SA
model (1)-(7) is also efficient.
ED
algorithm. Observe also that, the average solution times for CPLEX are less than a second that means the proposed
Table 5 shows the results for solving the r-HIMP on the TR data set. Here, it is also assumed that, the leader has
AC
CE
PT
already selected its p hubs based on the solution of the UMApHMP.
18
ACCEPTED MANUSCRIPT
Table 5: Results for the r-HIMP with the TR data set α = 0.4
10
14
α = 0.6
CPLEX
SA
α = 0.8
CPLEX
SA
CPLEX
SA
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
0
558.70
2.76
558.70
0.10
650.42
1.56
650.42
0.10
725.25
1.70
725.25
0.10
1
666.11
8.91
666.11
0.65
734.69
5.77
734.69
0.68
786.67
5.30
786.67
0.66
2
860.75
17.83
860.75
0.46
867.44
14.76
867.44
0.47
899.02
10.65
899.02
0.47
3
1147.94
15.52
1147.94
0.31
1080.74
23.57
1080.74
0.32
1096.33
13.33
1096.33
0.32
4
1557.13
13.62
1557.13
0.20
1235.56
23.39
1235.56
0.19
1244.05
18.21
1244.05
0.20
5
1783.13
14.94
1783.13
0.12
1783.13
13.71
1783.13
0.12
1783.13
11.86
1783.13
0.13
0
475.45
7.81
475.45
0.10
581.77
7.91
581.77
0.10
675.96
8.87
675.96
0.10
1
538.78
33.11
538.78
1.88
636.98
30.79
636.98
2.08
710.74
28.89
710.74
1.93
2
605.50
557.97
605.50
1.57
675.49
454.60
675.49
1.70
742.41
289.61
742.41
1.60
3
701.95
709.52
701.95
1.13
725.17
524.10
725.17
1.19
779.72
397.02
779.72
1.14
4
806.50
871.29
806.50
0.95
807.77
629.98
807.77
0.99
847.52
533.12
847.52
0.95
5
1019.68
1109.26
1019.68
0.64
992.61
969.82
992.61
0.67
1015.33
555.11
1015.33
0.66
0
427.26
23.06
427.26
0.10
545.55
22.32
545.55
0.10
653.22
21.53
653.22
0.10
1
490.97
124.32
490.97
7.15
600.84
123.79
600.84
7.26
673.14
127.34
673.14
7.27
2
549.96
1605.60
549.96
6.14
639.51
1216.38
639.51
6.32
692.83
1060.78
692.83
6.51
3
600.79
2339.64
600.69
5.35
686.87
1650.53
686.87
5.39
722.71
1460.90
722.71
5.56
4
659.22
2919.40
659.22
4.42
732.95
1848.92
732.95
4.59
760.15
2121.35
760.15
4.91
5
749.46
3569.30
749.46
4.10
810.10
2880.21
810.10
3.88
796.59
2526.00
796.59
4.08
819.40
774.65
819.40
1.96
821.53
580.11
821.53
2.00
866.93
510.64
866.93
2.03
Avg.
CR IP T
6
r
AN US
p
Note from the above table that the SA algorithm has obtained the optimal solution for all instances of the TR
M
data set. From a solution time perspective, it is shown that CPLEX solves all the instances in reasonable times. However, the solution times for the proposed SA are too small as the SA solves the instances in around 2 seconds, on average, which shows how efficient the proposed SA algorithm is even for solving large instances.
ED
Tables 6 and 7 show the results obtained by solving the BHIMP using the proposed SA algorithm for the CAB and TR data sets, respectively. Since it is not possible to directly solve the bilevel models by using CPLEX, we have
PT
only reported the results from the SA in these tables. Table 6: SA results for the BHIMP with the CAB data set
AC
CE
p
r
α = 0.4
α = 0.6
α = 0.8
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
1
1206.33
17.91
1263.61
17.88
1301.60
16.59
2
1511.81
7.95
1511.81
7.84
1511.81
8.40
1
1048.01
25.75
1087.43
22.08
1131.99
20.86
2
1315.41
18.35
1401.93
18.09
1384.68
18.45
3
1536.93
9.11
1596.49
9.41
1523.65
7.97
1
891.74
37.50
976.68
38.35
1033.03
40.45
2
1358.86
22.74
1242.10
22.91
1385.17
23.50
3
1457.88
16.56
1427.66
16.10
1654.53
17.09
4
1878.22
8.12
1536.93
7.98
1718.52
8.05
Avg.
1356.13
18.22
1338.29
17.85
1405.00
17.93
3
4
5
Results presented in Tables 6 and 7 reveal that the leader can decrease its transportation cost by acting rationally
19
ACCEPTED MANUSCRIPT
Table 7: SA results for the BHIMP with the TR data set r
10
14
α = 0.6
α = 0.8
CPU(s)
OF
CPU(s)
OF
CPU(s)
1
654.13
715.45
725.29
673.13
779.73
642.17
2
753.66
703.63
809.98
690.95
852.80
687.23
3
1024.90
410.38
957.23
409.62
976.98
409.25
4
1029.26
240.87
1036.56
241.33
1083.37
239.09
5
1140.65
107.46
1140.65
112.28
1140.65
109.56
1
532.82
1957.38
628.90
1940.12
707.41
1917.71
2
583.93
1708.12
671.73
1607.70
737.82
1626.26
3
671.15
1309.55
710.76
1309.87
769.82
1298.30
4
741.80
960.43
786.07
958.97
832.30
961.11
5
947.13
701.81
902.83
695.09
945.90
695.99
1
474.90
3962.61
578.61
3828.11
673.14
3928.73
2
509.41
3347.16
609.43
3723.12
692.83
3288.31
3
570.18
2918.01
661.48
2879.43
717.09
2717.38
4
592.15
2209.65
697.69
2274.19
754.83
2281.49
5
682.30
1992.16
727.22
1549.64
Avg.
CR IP T
6
α = 0.4 OF
AN US
p
772.64
1985.50
765.38
1823.65
779.32
1555.29
828.67
1508.41
in anticipation of the follower’s attack. For example, as can be seen in Table 7, in case p = 6 and r = 5 with α = 0.8, if the leader locates its hubs using the bilevel BHIMP model, then the transportation cost will reduce considerably
M
from 1783.13 (when the leader locates its hubs based on the UMApHMP) to 1140.65. Note also that, the CPU times for solving the BHIMP are much longer than the CPU times needed to solve the r-HIMP. This is due to the bilevel nature of the BHIMP that, as mentioned in the previous section, requires the proposed SA algorithm to solve the
ED
follower’s problem whenever it obtains a new solution for the leader (nested SA algorithm). However, the times are still very good for a strategic planning problem such as locating facilities in a competitive environment.
PT
5.4. Results for hub interdiction maximal covering problems Tables 8 and 9 present the results for solving the r-HIMCP with the CAB and TR data sets, respectively. In
CE
this case, it is assumed that the leader has already selected its p hubs based on the solution of the uncapacitated multiple allocation p-hub maximal covering problem (UMApHMCP) with the aim of maximizing the total covered flows in the network. As stated earlier, in this problem, the follower tries to minimize the leader’s total covered flows
AC
by interdicting r hubs from among the installed p hubs. To show the covering radius considered for each instance, a new column labeled as “β” is added in the corresponding tables. The columns labeled as “OF (%)” show the P total covered flows by the leader as percentage of the total flows within the network ( i,j wij ). As can be seen, the SA algorithm has again obtained the optimal solutions for all instances in less than one and three seconds, on
average, for the CAB and TR data sets, respectively. It should also be noted that the CPU times for CPLEX are also very short even for the large TR instances which shows the efficiency of the proposed mathematical model for the r-HIMCP. Observe from Tables 8 and 9 that, in case the leader decides on its hubs based on UMApHMCP, the follower can dramatically shrink the total value of covered flows by interdicting some r hubs. Moreover, by increasing the value 20
ACCEPTED MANUSCRIPT
Table 8: Results for the r-HIMCP with the CAB data set α = 0.4
5
3
4
5
3
4
5
SA
OF (%)
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
0
37.05
0.20
37.05
0.05
34.87
0.12
34.87
0.05
28.41
0.02
28.41
0.05
1
14.71
0.05
14.71
0.05
16.03
0.02
16.03
0.05
14.64
0.01
14.64
0.05
2
6.70
0.03
6.70
0.05
6.70
0.07
6.70
0.04
6.70
0.01
6.70
0.05
0
42.90
0.24
42.90
0.05
38.17
0.06
38.17
0.05
31.31
0.03
31.31
0.05
1
16.59
0.11
16.59
0.06
15.55
0.06
15.55
0.07
17.54
0.03
17.54
0.07
2
8.21
0.15
8.21
0.06
8.17
0.07
8.17
0.06
9.60
0.06
9.60
0.06
3
0.48
0.02
0.48
0.05
1.00
0.02
1.00
0.05
2.90
0.02
2.90
0.05
0
47.08
0.33
47.08
0.05
41.07
0.05
41.07
0.05
32.91
0.05
32.91
0.05
1
19.08
0.10
19.08
0.08
18.45
0.09
18.45
0.08
19.14
0.19
19.14
0.09
2
10.43
0.08
10.43
0.08
11.08
0.18
11.08
0.12
11.08
0.05
11.08
0.11
3
2.02
0.04
2.02
0.09
3.90
0.19
3.90
0.07
3.90
0.22
3.90
0.07
4
0.48
0.05
0.48
0.05
1.00
0.05
1.00
0.05
1.00
0.33
1.00
0.05
17.14
0.11
17.14
0.06
16.33
0.08
16.33
0.06
14.92
0.08
14.92
0.06
0
65.90
0.07
65.90
0.05
61.90
0.02
61.90
0.05
57.49
0.02
57.49
0.05
1
16.81
0.07
16.81
0.05
22.51
0.02
22.51
0.05
20.10
0.18
20.10
0.06
2
3.56
0.02
3.56
0.05
3.56
0.01
3.56
0.05
2.71
0.01
2.71
0.05
0
74.44
0.12
74.44
0.05
67.01
0.03
67.01
0.05
62.14
0.02
62.14
0.05
1
30.82
0.04
30.82
0.06
49.28
0.05
49.28
0.07
24.75
0.04
24.75
0.06
2
8.30
0.03
8.30
0.07
9.11
0.27
9.11
0.06
7.37
0.41
7.37
0.06
3
3.56
0.02
3.56
0.05
3.56
0.02
3.56
0.05
2.71
0.02
2.71
0.05
0
82.04
0.26
82.04
0.05
72.71
0.08
72.71
0.05
66.57
0.04
66.57
0.05
1
48.34
0.04
48.34
0.08
54.98
0.05
54.98
0.08
40.54
0.09
40.54
0.09
2
15.95
0.08
15.95
0.08
14.81
0.05
14.81
0.12
13.58
0.14
13.58
0.12
3
8.30
0.07
8.30
0.09
8.21
0.06
8.21
0.07
7.37
0.06
7.37
0.07
4
3.56
0.08
3.56
0.05
3.56
0.08
3.56
0.05
2.71
0.08
2.71
0.05
30.13
0.07
30.13
0.06
30.93
0.06
30.93
0.06
25.67
0.09
25.67
0.06
0
94.47
0.05
94.47
0.05
84.55
0.03
84.55
0.05
81.31
0.02
81.31
0.05
1
57.60
0.16
57.60
0.05
59.76
0.02
59.76
0.05
38.08
0.02
38.08
0.06
2
7.57
0.02
7.57
0.05
7.57
0.02
7.57
0.04
7.57
0.01
7.57
0.05
0
96.96
0.08
96.96
0.05
88.51
0.03
88.51
0.05
83.86
0.03
83.86
0.05
1
69.64
0.04
69.64
0.06
71.65
0.07
71.65
0.07
52.75
0.04
52.75
0.07
2
29.43
0.03
29.43
0.06
26.04
0.05
26.04
0.06
26.77
0.05
26.77
0.06
3
7.57
0.03
7.57
0.05
7.57
0.03
7.57
0.05
7.57
0.06
7.57
0.04
0 1 2
99.08
0.06
99.08
0.05
90.18
0.05
90.18
0.05
84.64
0.06
84.64
0.05
71.21
0.11
71.21
0.08
73.32
0.05
73.32
0.08
54.31
0.08
54.31
0.08
30.18
0.05
30.18
0.10
26.63
0.08
26.63
0.12
38.50
0.07
38.50
0.12
7.90
0.25
7.90
0.08
7.90
0.33
7.90
0.07
23.02
0.16
23.02
0.07
2.00
0.05
2.00
0.05
6.30
0.05
6.30
0.05
7.57
0.04
7.57
0.05
47.80
0.07
47.80
0.06
45.83
0.06
45.83
0.06
42.16
0.05
42.16
0.06
CE
3 4
Avg.
CPLEX
CPU(s)
Avg. 1500
α = 0.8 SA
OF (%)
Avg. 1000
CPLEX
CR IP T
4
α = 0.6 SA
AN US
3
CPLEX
M
500
r
ED
p
PT
β
AC
of r, the reduction in the coverage percentage gets more considerable, as expected. For example, in the CAB data set with p = 5, α = 0.4, β = 1500, and r = 0, if the leader chooses its hubs according to optimal solution of UMApHMCP, the flow coverage percentage will be 99.08%. However, as the follower decides to disrupt 1, 2, 3, or 4 hubs, the total covered flows by the leader will fall to 71.21%, 30.18%, 7.90%, or 2.00%, respectively. Another interesting observation is that, for smaller covering radius (β) values, the follower’s attack can be more effective with the same number of interdictions. For instance, consider the case p = 14, r = 5, and α = 0.4 in the TR data set. In case the covering radius is 700, the five interdictions drop the leader’s coverage from 92.48% to 35.24%, whereas for β = 1000, the leader’s coverage before and after the interdiction are 99.96% and 79.10%, respectively. Furthermore,
21
ACCEPTED MANUSCRIPT
the results show that as the discount factor value (α) grows, the coverage percentage falls down as a result of elevated inter-hub transportation cost. Table 9: Results for the r-HIMCP with the TR data set α = 0.4
10
14
6
10
10
14
Avg.
SA
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
0
30.04
0.91
30.04
1.00
23.79
0.61
23.79
1.00
17.89
0.79
17.89
1.00
1
15.35
0.95
15.35
1.32
15.90
0.50
15.90
1.32
12.93
0.43
12.93
1.31
2
10.21
0.49
10.21
0.95
10.78
0.59
10.78
0.90
9.03
0.50
9.03
0.93
3
6.27
0.46
6.27
0.63
6.20
0.65
6.20
0.63
5.69
0.49
5.69
0.71
4
3.79
0.73
3.79
0.41
3.54
0.48
3.54
0.39
3.66
0.71
3.66
0.39
5
1.71
0.57
1.71
0.28
1.71
0.60
1.71
0.26
1.79
0.48
1.79
0.26
0
40.93
1.83
40.93
1.00
29.76
1.36
29.76
1.00
21.74
1.98
21.74
1.00
1
30.87
3.61
30.87
3.72
20.92
1.62
20.92
3.68
17.01
1.59
17.01
3.71
2
20.48
2.08
20.48
3.15
15.98
1.71
15.98
3.17
12.67
1.62
12.67
3.11
3
14.47
1.84
14.47
2.26
11.62
1.85
11.62
2.28
9.68
6.16
9.68
2.25
4
10.36
1.63
10.36
1.96
9.41
1.59
9.41
1.92
7.64
1.63
7.64
1.92
5
6.77
7.20
6.77
1.35
7.53
2.03
7.53
1.39
6.17
1.56
6.17
1.30
0
48.42
3.33
48.42
1.00
33.85
2.99
33.85
1.00
24.03
2.93
24.03
1.00
1
39.33
3.13
39.33
7.09
27.92
2.43
27.92
7.14
20.05
2.38
20.05
7.22
2
32.08
2.49
32.08
6.09
22.10
8.43
22.10
6.11
17.17
2.44
17.17
6.20
3
25.87
2.69
25.87
5.33
4
20.83
2.89
20.83
4.46
5
15.99
4.11
15.99
3.81
20.76
2.27
20.76
2.54
0
73.25
0.66
73.25
1.00
1
52.17
0.70
52.17
1.38
2
31.48
0.73
31.48
0.98
3
19.20
0.74
19.20
0.63
4
9.74
0.87
9.74
0.43
5
3.70
0.71
3.70
0.27
0
86.89
1.70
86.89
1
66.73
2.40
66.73
2
48.18
6.71
48.18
3
33.98
2.47
33.98
4
26.42
5
19.57
0
92.48
1
83.62
2
70.11
18.46
2.46
18.46
5.47
14.39
3.61
14.39
5.22
14.52
2.42
14.52
4.46
11.44
4.93
11.44
4.51
12.20
4.32
12.20
3.76
9.97
4.56
9.97
3.77
15.89
2.03
15.89
2.54
12.53
2.15
12.38
2.54
60.98
0.54
60.98
1.00
52.19
0.47
52.19
1.10
47.40
0.61
47.40
1.37
42.16
0.83
42.16
1.31
34.97
0.70
34.97
0.95
31.74
2.81
31.74
0.95
19.34
0.52
19.34
0.64
18.20
0.53
18.20
0.62
11.05
0.50
11.05
0.42
10.73
0.56
10.73
0.42
4.89
0.45
4.89
0.28
4.89
4.77
4.89
0.27
1.00
69.03
2.88
69.03
0.98
57.12
1.49
57.12
0.90
3.77
57.74
2.27
57.74
3.74
49.64
1.97
49.64
3.74
3.18
49.58
2.05
49.58
3.14
43.18
1.90
43.18
3.12
2.28
37.12
2.57
37.12
2.31
32.82
1.86
32.82
2.29
2.34
26.42
1.97
25.06
1.59
25.06
1.92
22.63
1.62
22.63
1.89
1.83
19.57
1.37
20.50
1.91
20.50
1.35
19.75
1.78
19.75
1.34
2.86
92.48
1.00
73.83
4.56
73.83
1.00
59.77
2.46
59.77
1.00
4.77
83.62
7.00
68.21
3.93
68.21
7.11
57.84
3.72
57.84
7.03
4.48
70.11
6.03
65.02
3.35
65.02
6.14
55.37
2.80
55.37
6.21
3
52.14
4.70
52.14
5.23
55.73
4.30
55.73
5.54
49.06
2.96
49.06
5.24
4
42.82
4.39
42.82
4.50
42.21
4.40
42.21
4.49
39.61
3.76
39.61
4.49
5
35.24
4.18
35.24
3.77
31.77
4.94
31.77
3.82
31.39
4.47
31.39
3.79
47.09
2.62
47.09
2.54
43.02
2.33
43.02
2.56
37.67
2.26
37.67
2.53
CE AC
6
CPLEX
OF (%)
Avg. 1000
α = 0.8 SA
CPU(s)
PT
14
CPLEX
OF (%)
Avg. 700
α = 0.6 SA
CR IP T
6
CPLEX
AN US
400
r
M
p
ED
β
0
98.04
0.55
98.04
1.00
89.93
0.51
89.93
1.00
81.51
0.46
81.51
1.00
1
82.19
4.26
82.19
1.41
72.99
2.45
72.99
1.36
71.29
1.92
71.29
1.35
2
59.88
2.00
59.88
0.99
39.46
1.60
39.46
1.00
43.45
1.35
43.45
0.96
3
30.20
1.65
30.20
0.66
23.94
1.65
23.94
0.65
31.73
1.49
31.73
0.64
4
14.98
1.55
14.98
0.44
16.84
1.62
16.84
0.44
21.17
1.48
21.17
0.44
5
8.70
0.53
8.70
0.28
8.70
0.81
8.70
0.30
8.70
1.03
8.70
0.29
0
99.86
1.62
99.86
1.00
94.53
1.46
94.53
1.00
84.51
1.41
84.51
1.00
1
98.48
9.91
98.48
3.84
88.07
5.61
88.07
3.85
81.13
6.76
81.13
3.81
2
90.87
5.99
90.87
3.21
76.21
6.56
76.21
3.21
78.83
2.48
78.83
3.20
3
78.19
8.14
78.19
2.31
58.84
3.77
58.84
2.31
73.37
5.24
73.37
2.34
4
67.12
12.01
67.12
2.00
44.49
5.44
44.49
1.98
53.86
6.40
53.86
2.06
5
44.04
12.18
44.04
1.39
26.08
2.93
26.08
1.39
40.72
5.33
40.72
1.36
0
99.96
4.27
99.96
1.04
96.41
2.69
96.41
1.05
85.81
5.14
85.81
1.09
1
99.08
26.09
99.08
7.16
93.96
21.70
93.96
7.14
83.96
14.01
83.96
6.95
2
96.80
8.86
96.80
6.25
89.19
5.17
89.19
6.21
81.92
4.41
81.92
6.17
3
93.79
14.06
93.79
5.38
83.06
6.27
83.06
5.34
79.65
5.40
79.65
5.31
4
88.69
26.74
88.69
4.54
75.92
8.99
75.92
4.53
74.79
7.57
74.79
4.53
5
79.10
38.90
79.10
3.81
64.50
38.38
64.50
4.18
66.63
23.13
66.63
3.82
73.88
9.96
73.88
2.59
63.50
6.53
63.50
2.60
63.50
5.27
63.50
2.57
22
ACCEPTED MANUSCRIPT
Tables 10 and 11 summarize the results of solving the BHIMCP using the proposed SA algorithm for the CAB and TR data sets, respectively. It can be seen from these tables that considering the upcoming disruptions in leader’s decisions helps him/her prevent huge losses in the total covered flows after the follower’s interdictions. For instance, in the CAB data set with p = 5, α = 0.4, and β = 1500, if the leader chooses its hubs according to optimal solution of UMApHMCP, its total covered flows will plummet from 99.08% to 2.00% when the follower interdicts 4 hubs, whereas in case the leader decides on the location of its hubs based on the BHIMCP, its total covered flow will stay above 52% after the follower’s attack of the same intensity. From a CPU time perspective, it can be observed that
OF (%)
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
1
19.97
54.91
19.66
52.98
CR IP T
the proposed SA solves the BHIMCP with the CAB instances in less than 75 seconds, on average. The solution
19.66
53.00
2
13.38
25.14
7.72
24.93
11.16
25.45
1
16.59
96.39
22.19
97.74
19.66
101.18
2
8.71
57.24
8.85
56.09
11.03
56.90
times for the TR instances, though much longer than those of the CAB instances, are quite acceptable for such a large scale bilevel MIP problem.
Table 10: SA results for the BHIMCP with the CAB data set
500
3
4
5
r
7.69
27.05
7.69
25.53
7.72
25.64
33.39
161.72
26.45
169.25
19.14
158.41
2
14.32
95.30
20.66
93.62
13.76
96.69
3
8.38
62.41
4
2.90
11.79
55.52
8.85
52.73
25.19
4.04
25.66
1.57
26.85
13.93
67.26
14.34
66.81
12.51
66.32
1
44.21
57.08
39.97
56.10
44.84
53.35
2
39.97
27.11
36.83
27.03
8.75
27.30
1
39.20
110.34
54.09
104.11
52.57
136.71
2
17.29
56.37
39.97
56.89
15.40
56.38
3
9.47
31.08
3.66
28.27
2.71
28.71
1
64.48
172.69
56.32
168.69
54.06
163.03
2
16.23
98.54
19.55
97.54
23.07
98.36
3
12.05
55.32
16.56
58.90
13.13
61.31
PT
4
AC
CE
5
4
3.68
30.10
3.56
31.09
2.71
29.35
27.40
70.96
30.06
69.85
24.14
72.72
1
71.12
69.19
68.50
89.08
68.26
69.82
2
58.72
29.80
57.50
29.31
62.87
29.58
1
71.55
98.13
72.84
99.59
72.84
101.96
2
55.71
61.21
62.77
56.12
64.89
55.40
3
13.67
29.51
52.30
29.31
58.72
32.11
1
75.80
164.75
74.47
168.55
72.12
159.61
2
59.17
103.21
66.23
106.09
63.36
107.09
3
58.80
65.11
63.79
63.19
59.45
63.10
4
52.30
30.15
58.66
35.91
28.21
29.33
57.43
72.34
64.12
75.24
61.19
72.00
Avg.
1500
3
4
5
Avg.
α = 0.8
1
ED
3
α = 0.6
3
Avg. 1000
α = 0.4
AN US
p
M
β
23
ACCEPTED MANUSCRIPT
Table 11: SA results for the BHIMCP with the TR data set
6
10
14
6
CPU(s)
OF (%)
CPU(s)
1
20.98
702.29
15.96
669.57
12.93
678.65
2
14.64
427.46
11.97
418.87
10.08
417.69
3
8.62
320.42
8.15
326.11
7.47
310.66
4
5.61
175.24
5.09
194.12
4.07
224.90
5
1.87
97.00
2.20
78.29
2.20
77.90
1
31.53
2464.91
23.69
2402.08
17.52
2458.94
2
25.10
2060.72
17.98
2072.12
14.92
1966.24
3
17.41
1823.73
14.26
1755.72
11.75
1722.50
4
13.12
1495.01
11.42
1485.62
9.64
1453.09
5
11.40
1010.98
8.69
1008.55
7.41
1008.69
1
39.39
5428.22
28.00
5960.60
20.08
5899.21
2
32.28
5221.50
22.10
4146.33
17.96
5715.84
3
25.87
4948.01
18.76
4652.70
15.46
4581.11
4
22.92
4444.45
15.66
4109.87
12.56
3645.88
5
18.57
2977.73
12.54
3018.91
10.02
3074.12
19.29
2239.84
14.43
2153.30
11.60
2215.69
1
57.20
646.62
49.40
733.08
44.13
624.23
2
43.71
443.03
40.19
429.37
38.39
426.94
3
34.62
366.46
28.25
341.09
29.16
331.07
4
25.20
185.72
17.27
178.69
25.09
187.98
5
16.93
90.87
16.93
121.68
11.50
84.80
1
77.75
2487.87
62.80
2401.69
54.05
2531.53
2
68.79
2026.81
56.06
2046.96
48.84
2013.65
3
56.45
1749.89
47.54
1852.62
44.03
1791.32
4
41.69
1546.60
38.98
1494.46
37.29
1541.08
5
25.98
970.76
34.85
1063.05
33.32
981.34
1
88.15
5028.92
69.96
5856.44
57.97
5393.62
2
83.01
5004.26
66.16
5128.45
55.69
4537.91
3
68.82
4590.43
56.44
4402.63
51.78
4327.11
4
63.22
4052.19
53.58
3789.18
47.54
3811.55
48.63
3190.87
49.46
3235.19
44.18
3285.95
53.34
2158.75
45.86
2204.97
41.53
2124.67
1
90.79
642.43
81.20
643.09
76.90
639.54
2
66.77
454.90
74.48
440.53
68.75
434.67
3
55.36
350.25
49.45
349.93
57.99
325.81
4
41.09
210.33
48.66
214.66
46.97
223.90
5
37.64
106.99
37.64
96.61
37.18
91.45
1
98.55
3813.87
91.79
3751.44
82.53
3890.29
2
95.25
3155.20
87.03
3113.09
79.72
3207.09
3
86.64
3005.69
82.39
2814.50
75.18
2778.65
4
75.96
2307.07
75.93
2289.30
57.24
2391.02
5
52.19
1485.20
45.89
1543.53
56.07
1410.71
1
99.51
5330.66
94.67
5732.49
84.43
5653.48
2
98.75
4367.81
90.47
4941.17
82.80
5038.05
3
95.16
4149.00
88.98
4379.11
80.58
4241.21
4
89.63
3565.24
86.39
3844.61
75.28
3901.87
5
83.52
3054.74
79.26
3192.73
72.71
3013.17
77.79
2399.96
74.28
2489.79
68.96
2482.73
PT
14
5
Avg. 6
AC
CE
1000
10
14
Avg.
α = 0.8
OF (%)
ED
10
α = 0.6
CPU(s)
Avg. 700
α = 0.4 OF (%)
CR IP T
400
r
AN US
p
M
β
24
ACCEPTED MANUSCRIPT
5.5. Results for hub interdiction center problems Results obtained by solving the r-HICP by using CPLEX as well as the proposed SA algorithm on the CAB and TR data sets are presented in Tables 12 and 13, respectively. To solve the follower’s problem, it is assumed that the leader has already located its hubs based on the solution of the uncapacitated multiple allocation p-hub center problem (UMApHCP). In this case, the follower tries to maximize the maximum distance between all O/D pairs by attacking some of the leader’s hubs.
α = 0.4 p
r
α = 0.6
CPLEX
SA
α = 0.8
CPLEX
SA
CPLEX
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
0
2064.66
0.39
2064.66
0.05
2243.76
0.38
2243.76
0.05
1
5007.65
0.44
5007.65
0.06
4915.25
0.23
4915.25
0.07
2
5451.58
0.60
5451.58
0.06
5056.95
0.16
5056.95
0.06
0
1774.44
0.47
1774.44
0.05
2127.12
0.33
2127.12
0.05
1
2736.53
0.47
2736.53
0.07
2752.40
0.26
2752.40
0.07
2
5007.65
1.28
5007.65
0.06
5007.65
0.90
5007.65
0.07
3
5451.58
0.53
5451.58
0.06
5451.58
0.77
5451.58
0.06
0
1599.73
0.58
1599.73
0.05
1916.15
0.35
1916.15
1
2390.81
0.58
2390.81
0.10
2583.17
0.33
2
3012.90
1.55
3012.90
0.07
4598.85
3
5007.65
1.01
5007.65
0.06
4
5451.58
0.70
5451.58
Avg.
3746.39
0.72
3746.39
4
CPU(s)
OF
CPU(s)
2515.57
0.23
2515.57
0.05
4933.31
0.21
4933.31
0.06
5148.16
0.43
5148.16
0.06
2437.70
0.24
2437.70
0.05
4915.25
0.25
4915.25
0.08
5053.12
0.83
5053.12
0.07
5406.80
0.59
5406.80
0.06
0.05
2288.78
0.38
2288.78
0.05
2583.17
0.09
2824.84
0.29
2824.84
0.09
1.07
4598.85
0.07
4598.85
0.78
4598.85
0.07
5229.61
1.64
5229.61
0.08
5229.61
1.24
5229.61
0.08
0.07
5451.58
0.90
5451.58
0.06
5451.58
0.85
5451.58
0.05
0.06
3944.50
0.61
3944.50
0.06
4233.63
0.53
4233.63
0.06
M
5
SA
OF
AN US
OF 3
CR IP T
Table 12: Results for the r-HICP with the CAB data set
Table 13: Results for the r-HICP with the TR data set
CPLEX CPU(s)
0
1216.00
2.98
1
1522.20
20.03
2
1834.20
205.84
3
3014.00
4
14
CPLEX
α = 0.8 SA
CPLEX
SA
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
1216.00
0.10
1467.00
2.01
1467.00
0.10
1755.40
2.14
1755.40
0.10
1522.20
0.68
1728.20
14.19
1728.20
0.64
1876.60
12.39
1876.60
0.68
1834.20
0.50
2282.00
65.10
2282.00
0.47
2282.00
32.07
2282.00
0.45
126.58
3014.00
0.32
3336.00
67.39
3336.00
0.31
3032.00
46.80
3032.00
0.30
3600.00
109.33
3600.00
0.23
3474.00
90.72
3474.00
0.20
3474.00
96.39
3474.00
0.19
5
3748.00
18.81
3748.00
0.12
3748.00
14.92
3748.00
0.13
3748.00
14.37
3748.00
0.12
0
1114.20
7.86
1114.20
0.50
1434.80
7.72
1434.80
0.50
1755.40
8.85
1755.40
0.50
1
1317.00
71.19
1317.00
1.89
1610.40
69.47
1610.40
1.89
1870.60
45.06
1870.60
1.87
2
1514.80
1063.41
1514.80
1.59
1719.00
854.86
1719.00
1.60
1915.20
260.55
1915.20
1.57
3
2282.00
1286.67
2282.00
1.11
1856.00
909.92
1856.00
1.13
1934.20
1184.03
1934.20
1.13
4
2330.00
1271.66
2330.00
0.97
2096.00
825.82
2096.00
0.94
2096.00
467.43
2096.00
0.94
5
2352.00
2328.94
2352.00
0.68
2282.00
1866.93
2282.00
0.67
2352.00
1052.48
2352.00
0.66
0
1114.20
25.89
1114.20
1.00
1434.80
25.19
1434.80
1.10
1755.40
23.63
1755.40
1.00
1
1277.20
232.11
1277.20
4.00
1518.80
201.80
1518.80
3.51
1876.60
174.75
1876.60
3.52
2
1382.80
1111.73
1382.80
3.09
1608.80
1039.72
1608.80
3.10
1915.80
638.95
1915.80
3.03
3
1690.00
2442.43
1690.00
2.63
1750.40
1536.54
1750.40
2.63
1966.00
2831.04
1966.00
2.63
4
1964.00
3423.79
1964.00
2.24
1827.20
2649.48
1827.20
2.23
2096.00
2508.86
2096.00
2.23
5
2282.00
8728.92
2282.00
1.92
1964.00
8058.31
1964.00
1.90
2282.00
3419.66
2282.00
1.86
1975.25
1248.78
1975.25
1.30
2063.18
1016.67
2063.18
1.28
2221.28
712.19
2221.28
1.26
AC
10
α = 0.6 SA
PT
OF
CE
6
r
ED
α = 0.4 p
Avg.
Note that the proposed SA algorithm is able to obtain the optimal solutions of the r-HICP for all instances in the 25
ACCEPTED MANUSCRIPT
CAB and TR data sets in very short computational times. By comparing the CPU times between the SA heuristic and CPLEX, the efficiency of the SA algorithm becomes more evident. It can be seen that by increasing the discount factor (α), the average service level provided by the leader is diminished as the inter-hub transportation cost grow with the discount factor value. Observe also that, similar to previous cases, larger values for r result in considerably worsened objective function values for the leader. For example, for p = 14, r = 0, and α = 0.4 in Table 13, the service level provided by the leader as the maximum cost between the O/D pairs is 1114.20 based on the solution of the UMApHCP. However, after the follower attacks the network by disrupting 4 hubs, the maximum distance amounts
CR IP T
to 1964.00, showing an increment of 849.80 units, which can be considered as a serious damage to the service level provided by the leader.
Tables 14 and 15 show the results for solving the BHICP problem for the CAB and TR data sets, respectively. The results reported in these tables indicate that the proposed SA algorithm is able to solve the problem instances in very reasonable times.
3
4
5
r
α = 0.4 OF
CPU(s)
1
3012.90
85.10
2
3373.35
50.13
1
2454.39
143.80
2
3179.67
81.02
3
3529.58
41.18
1
2200.26
2
3012.90
3
4003.38
α = 0.6
α = 0.8
OF
CPU(s)
OF
CPU(s)
3012.90
77.32
3179.67
84.74
3471.87
50.06
3373.35
49.72
2645.59
141.91
2720.59
150.77
4072.25
79.10
3082.54
83.19
5451.58
49.73
3707.23
39.98
M
p
AN US
Table 14: SA results for the BHICP with the CAB data set
155.62
2427.42
165.09
2606.41
162.91
134.76
3529.58
131.45
3471.87
134.18
92.03
89.10
3745.39
92.83
4180.17
59.12
4598.85
44.58
4072.25
46.90
Avg.
3216.29
93.64
3726.64
92.04
3328.81
93.91
ED
4329.71
4
PT
Similar to the previous cases, it can be observed that the leader can benefit by solving the BHICP to reduce the adverse effects of the follower’s disruptive action. For instance, as stated above, for p = 14 and α = 0.4 in the TR
CE
data set, if the leader locates its hubs based on the UMApHCP, the follower will diminish the maximum O/D cost by 849.80 units using 4 interdictions, whereas acting based on the BHICP, the leader’s service level will be worsened at most by 409.40 units. Therefore, it is highly beneficial for the leader to solve the bilevel model and locate its hubs
AC
accordingly.
5.6. Optimal hub locations on map In this section, we graphically illustrate how the optimal hub locations for the leader and optimal interdictions
by the follower change as the leader makes its decisions with and without being aware of the follower’s upcoming attack. The analysis is conducted for all the three variants of the hub interdiction problem. Figure 1 depicts the optimal solutions for the leader and the follower for the hub interdiction median problem on the CAB data set with p = 5, r = 2, and α = 0.8. The leader’s installed hubs are shown as squares, whereas the follower’s interdicted hubs are shown as cross marks. In part (a), it is assumed that the leader is unaware of the follower’s subsequent attack 26
ACCEPTED MANUSCRIPT
Table 15: SA results for the BHICP with the TR data set r
10
14
α = 0.6
α = 0.8
CPU(s)
OF
CPU(s)
OF
CPU(s)
1
1468.20
792.38
1663.80
788.61
1838.40
809.84
2
1590.00
543.34
1772.00
533.11
1912.00
517.02
3
2338.00
320.73
2352.00
321.09
2282.00
311.31
4
2684.00
200.99
2352.00
207.41
2732.00
210.34
5
2684.00
93.11
2512.00
92.16
3190.00
93.79
1
1254.40
2410.84
1518.80
2406.59
1791.40
2388.43
2
1420.40
2483.43
1608.80
2238.96
1870.60
2063.05
3
1488.40
1876.38
1715.00
1733.65
1882.60
2109.98
4
2330.00
1449.11
1906.60
1678.11
1964.20
1820.11
5
2352.00
1022.32
2282.00
1201.85
2338.00
972.10
1
1254.40
4624.10
1518.80
4677.14
1791.40
4600.28
2
1376.40
3984.46
1605.60
4011.35
1870.60
4126.63
3
1468.20
4261.13
1676.60
4283.20
1882.00
4310.68
4
1523.60
4997.60
1737.80
3590.14
1964.00
4097.01
5
1836.00
3107.19
1804.53
2144.47
Avg.
CR IP T
6
α = 0.4 OF
AN US
p
1964.00
3244.09
2282.00
3951.50
1879.05
2067.16
2106.08
2158.80
and locates its hubs based on the solution of the UMApHMP, whereas part (b) depicts the optimal solutions in case
CE
PT
ED
M
the leader locates its hubs based on the BHIMP. As can be seen from the Figure 1, the optimal hub locations as well
(a) leader acts based on the UMApHMP
(b) leader acts based on the BHIMP
AC
Figure 1: Optimal hub locations and interdictions for hub median problem with p = 5, r = 2 and α = 0.8
as interdictions are different for the cases in which the leader acts based on the r-HIMP and BHIMP. This shows the importance of solving the BHIMP, and locating the hubs accordingly, for the leader. Figure 2 depicts the optimal solutions for the leader and the follower for the hub interdiction maximal covering problem with β = 1500. Here, it can also be seen that the optimal network configurations are quite different as the leader acts whether based on the r-HIMCP or BHIMCP. Finally, the optimal solutions for the leader and the follower for the hub interdiction center problem are depicted in Figure 3. Similar to the two previous cases, it can also be observed that the optimal set of hubs for the leader as well as the optimal interdictions for the follower are different 27
CR IP T
ACCEPTED MANUSCRIPT
(a) leader acts based on the UMApHMCP
(b) leader acts based on the BHIMCP
Figure 2: Optimal hub locations and interdictions for hub maximal covering problem with p = 5, r = 2, α = 0.8 and β = 1500
ED
M
AN US
based on whether the leader acts in anticipation of the follower’s attack or not.
(a) leader acts based on the UMApHCP
(b) leader acts based on the BHICP
PT
Figure 3: Optimal hub locations and interdictions for hub center problem with p = 5, r = 2 and α = 0.8
CE
5.7. Evaluating the leader’s objective value improvements Results presented in the previous sections reveal that, when the leader acts in anticipation of follower’s attack, i.e., he/she acts based on the solution of the BHIMP, BHIMCP, or BHICP, the effect of follower’s disruptive action
AC
is less severe than that of the case where the leader is unaware of the upcoming attack and acts based on the solution of the UMApHMP, UMApHMCP, or UMApHCP, respectively. In other words, acting based on the bilevel models, creates value for the leader. As a last part of our numerical analyses, in this section, we evaluate the changes (in percent) in the leader’s objective function value when he/she acts according to the bilevel models based on the presented results in preceding sections. To this end, we define and use a metric called the Objective Function Change Percentage (OFCP), as follows: OFCP =
OFB − OFS × 100% OFS
28
ACCEPTED MANUSCRIPT
in which OFB is the average value of leader’s objective function when he/she acts based on the bilevel models (BHIMP, NBHIMCP, or BHICP), whereas OFS is the average value of leader’s objective function when the locations are selected according to the single level models (UMApHMP, UMApHMCP, or UMApHCP). The averages are taken with respect to the three discount factor (α) values. The OFCP values for different instances of the CAB data set with the three variants of the problem are presented in Table 16. The first two columns show the number of hubs opened by the leader and the number of hubs interdicted by the follower, respectively. The third column give the OFCP values for the hub median variant of the problem. Similarly, the next three columns show the OFCP values
the OFCP values for the hub center variant of the problem.
CR IP T
for the hub maximal covering variant under the three covering radius (β) values. Finally, the last column present
Table 16: Average OFCP values for the CAB data set
p
Median
Maximal covering β=500
β=1000
β=1500
Center
-0.20%
30.65%
117.13%
33.74%
-38.04%
2
-57.02%
60.50%
770.30%
688.60%
-34.73%
1
-4.47%
17.63%
39.11%
11.95%
-19.24%
2
-9.93%
10.05%
193.22%
122.97%
-31.42%
3
-55.87%
427.40%
61.14%
449.05%
-22.20%
1
-11.26%
39.37%
21.55%
11.84%
-7.24%
2
-18.41%
49.56%
32.72%
98.05%
-17.99%
3
-21.45%
195.52%
74.79%
368.93%
-21.91%
4
-51.35%
243.15%
1.22%
776.94%
-21.42%
Avg.
-25.55%
119.31%
145.69%
284.67%
-23.80%
4
ED
5
AN US
1
M
3
r
Observe that for the CAB data set, the average improvement in the leader’s objective value is 25.55% for hub median problem and 23.80% for the hub center variant. For the hub maximal covering variant, on the other hand, the
PT
improvements are substantially higher than the other two variants. The leader’s objective has improved respectively 119.31%, 145.69%, and 284.67%, on average, for β = 500, 1000, and 1500 by acting based on the bilevel models.
CE
An interesting observation is that, the average OFCP values get larger as the follower disrupts more hubs. In other words, the benefit of using bilevel models are more considerable for the leader when facing a powerful interdictor. Table 17 shows the OFCP values for different instances of the TR data set with the three variants of the problem.
AC
As can be seen from the Table 17, the average improvements in the value of leader’s objective with the TR data
set are 7.65% and 10.51% for hub median and center variants, respectively. However, for the hub maximal covering variant, the leader’s objective has improved 16.57%, 49.18%, and 49.85%, on average, for β = 400, 700, and 1000, respectively. Considering the improvement percentages for the CAB and TR data sets, one can observe that the improvements are larger for the CAB data set (especially for the maximal covering variant) compared to those of the TR data set. This difference can partially be explained by the existence of some very high traffic nodes in the CAB data set which if interdicted, inflict a high damage to the network’s efficiency. Therefore, proper location of hubs by the leader as a result of acting based on the bilevel model, can prevent such harmful damages. Another
29
ACCEPTED MANUSCRIPT
Table 17: Average OFCP values for the TR data set
r
14
Center
β=400
β=700
β=1000
-1.29%
12.88%
6.35%
9.90%
-3.05%
2
-8.02%
22.22%
24.54%
47.07%
-17.57%
3
-11.00%
33.48%
62.20%
89.59%
-25.69%
4
-21.99%
34.39%
114.34%
158.01%
-26.36%
5
-36.03%
20.35%
236.50%
330.88%
-25.42%
1
-0.92%
5.73%
11.77%
1.94%
-4.86%
2
-1.48%
18.05%
23.24%
6.54%
3
-2.50%
21.39%
42.44%
16.07%
-16.24%
4
-4.13%
24.70%
59.17%
26.39%
-4.92%
5
-7.65%
34.34%
57.39%
39.07%
-0.20%
1
-2.17%
0.19%
3.06%
0.58%
-2.31%
2
-3.75%
1.39%
7.54%
1.53%
3
-3.07%
2.33%
12.81%
3.20%
-7.02%
4
-5.00%
9.30%
31.85%
4.97%
-11.24%
5
-5.76%
7.78%
44.58%
12.02%
-6.83%
-7.65%
16.57%
49.18%
49.85%
-10.51%
1
10
Maximal covering
Avg.
CR IP T
6
Median
-4.84%
-1.12%
AN US
p
M
important factor which may contribute to this difference is the limited number of candidate nodes for locating the hubs in the TR data set (|H| = 22). This limited set of candidate hub nodes make the action space of the decision makers limited, which in turn results in less difference between the set of leader’s hubs when he/she decides based
ED
on either of single level or bilevel models. Nevertheless, based on the results presented in the above tables, it can be concluded that the leader will substantially benefit via acting based on the bilevel models by reducing the adverse
6. Conclusions
PT
effects of the follower’s attack.
CE
In this paper, we addressed three different variants of the hub interdiction problem which has attracted the attention of researchers in recent years. To this end, we considered the multiple allocation p-hub median, the p-
AC
hub maximal covering, and the p-hub center problems under intentional disruptions by assuming a leader-follower Stackelberg game where the decision maker who makes the initial location decisions is called the leader, whereas the one who decides to interdict the installed hubs is called the follower. New formulations were proposed in each case which model the problems from both the leader’s and the follower’s perspectives. Furthermore, different types of CA constraints were examined in order to identify the one with the best performance in terms of solution time by a commercial solver. To solve the problems, we proposed efficient simulated annealing (SA) heuristics. Extensive computational experiments based on the CAB and TR data sets were conducted to analyze different properties of the problems and to evaluate the performance of the proposed SA heuristics as well as the mathematical formulations. The results demonstrated the high efficiency of the presented models and the SA algorithms both in terms of solution
30
ACCEPTED MANUSCRIPT
quality and computational time. Our research can be extended in several ways. An interesting research direction is to consider capacity restrictions on hubs or a limited budget for installing and/or disrupting the hubs. In addition, the assumption that an attack on a given hub is always successful causing complete disruption can be relaxed resulting in a problem in which attacks are successful only with a given probability or only cause partial disruptions. Moreover, developing exact methods such as branch-and-bound or decomposition algorithms might be of interest for the researchers.
CR IP T
Acknowledgments The authors sincerely thank the three anonymous referees for their valuable comments and suggestions that have helped improve the paper considerably.
References
AN US
˙ G¨ [1] Akg¨ un, I., um¨ u¸sbuga, F., Tansel, B., 2015. Risk based facility location by using fault tree analysis in disaster management. Omega 52, 168 – 179.
[2] Aksen, D., Piyade, N., Aras, N., 2010. The budget constrained r-interdiction median problem with capacity expansion. Central European Journal of Operations Research 18 (3), 269–291. [3] Aksen, D., engl Akca, S., Aras, N., 2014. A bilevel partial interdiction problem with capacitated facilities and
M
demand outsourcing. Computers & Operations Research 41, 346 – 358.
[4] Alumur, S., Kara, B. Y., 2008. Network hub location problems: The state of the art. European Journal of
ED
Operational Research 190 (1), 1–21.
[5] Bard, J. F., 1998. Practical bilevel optimization: algorithms and applications. Kluwer Academic Publishers,
PT
Dordrecht, The Netherlands.
[6] Berman, O., Drezner, Z., Wesolowsky, G. O., 2003. Locating service facilities whose reliability is distance
CE
dependent. Computers & Operations Research 30 (11), 1683–1695. [7] Berman, O., Krass, D., Menezes, M. B., 2007. Facility reliability issues in network p-median problems: strategic
AC
centralization and co-location effects. Operations Research 55 (2), 332–350. [8] Berman, O., Krass, D., Menezes, M. B., 2009. Locating facilities in the presence of disruptions and incomplete information. Decision Sciences 40 (4), 845–868.
[9] Berman, O., Krass, D., Menezes, M. B., 2013. Location and reliability problems on a line: Impact of objectives and correlated failures on optimal location patterns. Omega 41 (4), 766 – 779. [10] Campbell, J. F., 1992. Location and allocation for distribution systems with transshipments and transportion economies of scale. Annals of operations research 40 (1), 77–99.
31
ACCEPTED MANUSCRIPT
[11] Campbell, J. F., 1994. Integer programming formulations of discrete hub location problems. European Journal of Operational Research 72 (2), 387–405. [12] Campbell, J. F., 1996. Hub location and the p-hub median problem. Operations Research 44 (6), 923–935. [13] Campbell, J. F., O’Kelly, M. E., 2012. Twenty-five years of hub location research. Transportation Science 46 (2), 153–169.
rithm. Journal of optimization theory and applications 45 (1), 41–51.
CR IP T
ˇ [14] Cern` y, V., 1985. Thermodynamical approach to the traveling salesman problem: An efficient simulation algo-
[15] Chaharsooghi, S., Momayezi, F., Ghaffarinasab, N., 2017. An adaptive large neighborhood search heuristic for solving the reliable multiple allocation hub location problem under hub disruptions. International Journal of Industrial Engineering Computations 8 (2), 191–202.
[16] Church, R., Scaparra, M. P., 2007. Analysis of Facility Systems’ Reliability When Subject to Attack or a Natural
AN US
Disaster. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 221–241.
[17] Church, R. L., Cohon, J. L., 1976. Multiobjective location analysis of regional energy facility siting problems. Tech. rep., Brookhaven National Lab., Upton, NY (USA).
[18] Church, R. L., Scaparra, M. P., Middleton, R. S., 2004. Identifying critical infrastructure: the median and
M
covering facility interdiction problems. Annals of the Association of American Geographers 94 (3), 491–502. [19] Contreras, I., 2015. Hub location problems. In: Location Science. Springer International Publishing, pp. 311–344.
ED
[20] Contreras, I., Fernndez, E., Reinelt, G., 2012. Minimizing the maximum travel time in a combined model of facility location and network design. Omega 40 (6), 847 – 860. [21] Cui, T., Ouyang, Y., Shen, Z.-J. M., 2010. Reliable facility location design under the risk of disruptions.
PT
Operations Research 58 (4-part-1), 998–1011.
CE
[22] Dempe, S., 2002. Foundations of bilevel programming. Springer Science & Business Media. [23] Dobson, G., Karmarkar, U. S., 1987. Competitive location on a network. Operations Research 35 (4), 565–574.
AC
[24] Drezner, Z., 1987. Heuristic solution methods for two location problems with unreliable facilities. Journal of the Operational Research Society 38 (6), 509–514.
[25] Ebery, J., Krishnamoorthy, M., Ernst, A., Boland, N., 2000. The capacitated multiple allocation hub location problem: Formulations and algorithms. European Journal of Operational Research 120 (3), 614–631.
[26] Ernst, A. T., Hamacher, H., Jiang, H., Krishnamoorthy, M., Woeginger, G., 2009. Uncapacitated single and multiple allocation p-hub center problems. Computers & Operations Research 36 (7), 2230–2241. [27] Ernst, A. T., Krishnamoorthy, M., 1996. Efficient algorithms for the uncapacitated single allocation p-hub median problem. Location science 4 (3), 139–154. 32
ACCEPTED MANUSCRIPT
[28] Ernst, A. T., Krishnamoorthy, M., 1998. Exact and heuristic algorithms for the uncapacitated multiple allocation p-hub median problem. European Journal of Operational Research 104 (1), 100–112. [29] Espejo, I., Mar´ın, A., Rodr´ıguez-Ch´ıa, A. M., 2012. Closest assignment constraints in discrete location problems. European Journal of Operational Research 219 (1), 49–58. [30] Farahani, R. Z., Hekmatfar, M., Arabani, A. B., Nikbakhsh, E., 2013. Hub location problems: A review of models, classification, solution techniques, and applications. Computers & Industrial Engineering 64 (4), 1096–
CR IP T
1109. [31] Gerrard, R. A., Church, R. L., 1996. Closest assignment constraints and location models: properties and structure. Location Science 4 (4), 251–270.
[32] Ghaffarinasab, N., Motallebzadeh, A., Jabarzadeh, Y., Kara, B. Y., 2018. Efficient simulated annealing based solution approaches to the competitive single and multiple allocation hub location problems. Computers &
AN US
Operations Research 90 (Supplement C), 173 – 192.
[33] Hamacher, H. W., Meyer, T., 2006. Hub cover and hub center problems. Technische Universit¨ at Kaiserslautern, Fachbereich Mathematik.
[34] Hong, J.-D., Xie, Y., Jeong, K.-Y., 2012. Development and evaluation of an integrated emergency response facility location model. Journal of Industrial Engineering and Management 5 (1), 4.
M
[35] Jabalameli, M. S., Barzinpour, F., Saboury, A., Ghaffari-Nasab, N., 2012. A simulated annealing-based heuristic for the single allocation maximal covering hub location problem. International Journal of Metaheuristics 2 (1),
ED
15–37.
[36] Kara, B., 1999. Modeling and analysis of issues in hub location problems. Doctor of Philosophy Thesis, Bilkent
PT
University, Ankara, Turkey.
[37] Kara, B., Tansel, B., 2003. The single-assignment hub covering problem: Models and linearizations. Journal of
CE
the Operational Research Society 54 (1), 59–64. [38] Kara, B. Y., Tansel, B. C., 2000. On the single-assignment p-hub center problem. European Journal of Opera-
AC
tional Research 125 (3), 648–655. [39] Kirkpatrick, S., Vecchi, M. P., et al., 1983. Optimization by simmulated annealing. science 220 (4598), 671–680. [40] Lee, S.-D., 2001. On solving unreliable planar location problems. Computers & Operations Research 28 (4), 329–344.
[41] Lei, T. L., 2010. Location modeling utilizing closest and generalized closest transport/interaction assignment constructs. University of California, Santa Barbara. [42] Lei, T. L., 2013. Identifying critical facilities in hub-and-spoke networks: A hub interdiction median problem. Geographical Analysis 45 (2), 105–122. 33
ACCEPTED MANUSCRIPT
[43] Lei, T. L., Tong, D., 2013. Hedging against service disruptions: an expected median location problem with site-dependent failure probabilities. Journal of Geographical Systems 15 (4), 491–512. [44] Li, X., Ouyang, Y., 2010. A continuum approximation approach to reliable facility location design under correlated probabilistic disruptions. Transportation research part B: methodological 44 (4), 535–548. [45] Liberatore, F., Scaparra, M. P., Daskin, M. S., 2011. Analysis of facility protection strategies against an uncertain number of attacks: the stochastic r-interdiction median problem with fortification. Computers & Operations
CR IP T
Research 38 (1), 357–366. [46] Lim, M., Daskin, M. S., Bassamboo, A., Chopra, S., 2010. A facility reliability problem: formulation, properties, and algorithm. Naval Research Logistics (NRL) 57 (1), 58–70.
[47] Losada, C., Scaparra, M. P., Church, R. L., Daskin, M. S., 2012. The stochastic interdiction median problem with disruption intensity levels. Annals of Operations Research 201 (1), 345–365.
AN US
[48] Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., Teller, E., 1953. Equation of state calculations by fast computing machines. The journal of chemical physics 21 (6), 1087–1092. [49] O’Hanley, J. R., Church, R. L., 2011. Designing robust coverage networks to hedge against worst-case facility losses. European Journal of Operational Research 209 (1), 23–36.
[50] O’Kelly, M. E., 1987. A quadratic integer program for the location of interacting hub facilities. European Journal
M
of Operational Research 32 (3), 393–404.
[51] Parvaresh, F., Golpayegany, S. H., Husseini, S. M., Karimi, B., 2013. Solving the p-hub median problem under
ED
intentional disruptions using simulated annealing. Networks and Spatial Economics 13 (4), 445–470. [52] Parvaresh, F., Husseini, S. M., Golpayegany, S. H., Karimi, B., 2014. Hub network design problem in the
PT
presence of disruptions. Journal of Intelligent Manufacturing 25 (4), 755–774. [53] Peker, M., Kara, B. Y., 2015. The p-hub maximal covering problem and extensions for gradual decay functions.
CE
Omega 54, 158–172.
[54] Peng, P., Snyder, L. V., Lim, A., Liu, Z., 2011. Reliable logistics networks design with facility disruptions.
AC
Transportation Research Part B: Methodological 45 (8), 1190–1211. [55] Rojeski, P., ReVelle, C., 1970. Central facilities location under an investment constraint. Geographical Analysis 2 (4), 343–360.
[56] Saboury, A., Ghaffari-Nasab, N., Barzinpour, F., Jabalameli, M. S., 2013. Applying two efficient hybrid heuristics for hub location problem with fully interconnected backbone and access networks. Computers & Operations Research 40 (10), 2493–2507. [57] Scaparra, M. P., Church, R. L., 2008. A bilevel mixed-integer program for critical infrastructure protection planning. Computers & Operations Research 35 (6), 1905–1923. 34
ACCEPTED MANUSCRIPT
[58] Scaparra, M. P., Church, R. L., 2015. Location problems under disaster events. In: Laporte, G., Nickel, S., Saldanha da Gama, F. (Eds.), Location Science. Springer International Publishing, pp. 623–642. [59] Silva, M. R., Cunha, C. B., 2017. A tabu search heuristic for the uncapacitated single allocation p-hub maximal covering problem. European Journal of Operational Research 262 (3), 954 – 965. [60] Skorin-Kapov, D., Skorin-Kapov, J., O’Kelly, M., 1996. Tight linear programming relaxations of uncapacitated p-hub median problems. European Journal of Operational Research 94 (3), 582–593.
Transportation Science 39 (3), 400–416.
CR IP T
[61] Snyder, L. V., Daskin, M. S., 2005. Reliability models for facility location: the expected failure cost case.
[62] Snyder, L. V., Scaparra, M. P., Daskin, M. S., Church, R. L., 2006. Planning for disruptions in supply chain networks. Tutorials in operations research 2, 234–257.
AN US
[63] Tan, P. Z., Kara, B. Y., 2007. A hub covering model for cargo delivery systems. Networks 49 (1), 28–39. [64] Wagner, B., 2008. Model formulations for hub covering problems. Journal of the Operational Research Society 59 (7), 932–938.
[65] Wagner, J., Falkson, L., 1975. The optimal nodal location of public facilities with price-sensitive demand. Geographical Analysis 7 (1), 69–83.
M
[66] Yaman, H., Elloumi, S., 2012. Star p-hub center problem and star p-hub median problem with bounded path
AC
CE
PT
ED
lengths. Computers & Operations Research 39 (11), 2725–2732.
35
Hub interdiction problem variants: Models and metaheuristic solution algorithms Nader Ghaffarinasab, Alireza Motallebzadeh PII: DOI: Reference:
S0377-2217(17)31069-X 10.1016/j.ejor.2017.11.058 EOR 14846
To appear in:
European Journal of Operational Research
Received date: Revised date: Accepted date:
2 September 2016 23 November 2017 23 November 2017
Please cite this article as: Nader Ghaffarinasab, Alireza Motallebzadeh, Hub interdiction problem variants: Models and metaheuristic solution algorithms, European Journal of Operational Research (2017), doi: 10.1016/j.ejor.2017.11.058
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 • New hub interdiction problems are introduced to the hub location literature. • New bilevel and single level mathematical formulations are proposed for the problems. • Different multi-stage CA constraints are presented and analyzed.
AC
CE
PT
ED
M
AN US
• Extensive computational study are carried out on the problems.
CR IP T
• Efficient simulated annealing based heuristics are proposed to solve the problems.
1
ACCEPTED MANUSCRIPT
Hub interdiction problem variants: Models and metaheuristic solution algorithms Nader Ghaffarinasaba,∗, Alireza Motallebzadehb a Department b Faculty
of Industrial Engineering, University of Tabriz, Tabriz, Iran of Economics, Management and Business, University of Tabriz, Tabriz, Iran
CR IP T
Abstract Hub location problem (HLP) is one of the strategic problems encountered in designing transportation and telecommunication networks. Regardless of the considered objective in design of hub networks, such as cost minimization or service level maximization, the located hubs can be subject to natural or intentional disruptions after installation. In this paper, we address the multiple allocation p-hub median, p-hub maximal covering, and p-hub center problems under intentional disruptions. In each case, the problem is considered as a Stackelberg game where the leader locates
AN US
p hubs to optimize his/her objective function, whereas the follower tries to identify and interdict r hubs that their loss would diminish the network performance the most. Bilevel and single level mathematical formulations are presented to model the problem from the leader’s and the follower’s perspectives. Furthermore, efficient Simulated Annealing (SA) heuristics are proposed for solving the problems. Extensive computational experiments show the capability of the proposed SA algorithms to obtain the optimal solutions in short computational times. Some managerial insights are also derived based on the obtained numerical results.
M
Keywords: Location, Hub interdiction problems, Competitive models, Bilevel programming, Simulated annealing.
ED
1. Introduction
Hub networks play a major role in reducing cost and enhancing service level in many transportation, telecommunications and computer networks. These networks provide efficient service between many origins-destination (O/D)
PT
pairs of nodes via a set of hubs that serve as switching and flow consolidation points. The hub location problem (HLP) deals with locating hub facilities and allocating non-hub nodes (spokes) to the installed hubs in order to route
CE
the traffic between O/D pairs [4]. Regarding the way the non-hub nodes are allocated to the hubs, there are two main types of allocation in the hub networks called single and multiple allocation. In a single allocation network, all the incoming and outgoing traffic to/from every demand node is routed through a single hub, whereas in a multiple
AC
allocation network, which is the underlying network in this paper, each demand node can receive and send flow through more than one hub. Most of the papers published in the literature of the HLP focus on location and network design decisions in a deterministic environment where the reliability and robustness issues in cases of service disruption are not addressed. However, in real world applications, hub facilities can be lost or stop operating for a variety of reasons, such as natural disasters, inclement weather, equipment malfunctioning, labor strikes, and intentional attacks [42]. In other ∗ Corresponding
author Email addresses: [email protected], [email protected] (Nader Ghaffarinasab ), [email protected] (Alireza Motallebzadeh)
Preprint submitted to European Journal of Operational Research
December 1, 2017
ACCEPTED MANUSCRIPT
words, disasters whether intentional or accidental, can sometimes render a system inoperable or inefficient. In this paper, we consider the multiple allocation p-hub median, p-hub maximal covering, and p-hub center problems under intentional disruptions. The problems are tackled as leader-follower Stackelberg games where the decision maker who makes the initial location decisions is called leader and the interdictor is called follower. Accordingly, for each of the three HLP variants, two models are developed: a single level mixed integer programming (MIP) model on behalf of the follower and a bilevel model on behalf of the leader. To the best of our knowledge, the hub interdiction problem under covering and center objectives are modeled and solved for the first time in the literature in this paper.
CR IP T
We also introduce some multi-stage closest assignment (CA) constraints and compare their performance with those of other existing constraints in terms of solution time by a commercial solver. Furthermore, to solve the above mentioned problems, efficient solution algorithms based on simulated annealing (SA) are proposed. In case of the follower’s problem, the single level models are also solved using a standard optimization package that enables us to evaluate the performance of the proposed SA algorithms as well as the proposed mathematical models.
The remainder of this paper is organized as follows. The next section reviews the relevant literature to the problem.
AN US
In Section 3, we will present the mathematical formulations for the problems. The proposed SA algorithms are presented in Section 4. Computational experiments and corresponding results are presented in Section 5. Finally, some concluding remarks and directions for future works are given in Section 6.
2. Literature Review
M
2.1. Hub location problem
Hub location models can be categorized according to the type of objective that is considered. The most frequently addressed models are minisum hub location models, such as the p-hub median problem and the uncapacitated hub
ED
location problem (with fixed costs), that focus on minimizing the overall setup and operational costs. O’Kelly [50] developed the first quadratic mathematical formulation for uncapacitated single allocation p-hub median problem
PT
(USApHMP). In case of uncapacitated multiple allocation p-hub median problem (UMApHMP), the first mathematical model was introduced by Campbell [10]. Linear integer programming formulations for single and multiple allocation p-hub median problem were presented in [11]. Skorin-Kapov et al. [60] proposed new models for unca-
CE
pacitated single and multiple allocation hub location problems having tight linear programming relaxations. Ernst and Krishnamoorthy [27] introduced an efficient flow based mathematical formulation for the single allocation p-hub
AC
median problem that was later extended to model the multiple allocation p-hub median problem [28]. Kara proved that the p-hub median problem belongs to the class of NP-hard problems [36]. Another important variant of the HLP is the p-hub maximal covering problem which aims at locating a fixed number of hub facilities in such a way that the total covered flow associated with the O/D pairs is maximized. The p-hub maximal covering problem is closely related to the hub set covering problem in which the objective is to cover all the O/D pairs using a minimum number of opened hubs. Campbell [11] formulated the p-hub maximal covering and the hub set covering problems as linear MIPs. Kara and Tansel [37] proposed new models for the single allocation hub set covering problem and proved that it is NP-hard. Later, Wagner [64] introduced novel formulations for both single and multiple allocation hub set covering problems. Peker and Kara [53] addressed the p-hub maximal covering
3
ACCEPTED MANUSCRIPT
problem considering partial coverage, an extension of the traditional binary coverage. The third important variant of the HLP is the p-hub center problem which is defined based on the minimax criterion in which the maximum service cost (distance, time, etc.) between all the O/D pairs is minimized. Campbell [11] was first to introduce the p-hub center problem in the hub location literature. He presented formulations for both single and multiple allocation p-hub center problems. Kara and Tansel [38] addressed different linear formulations for the single allocation p-hub center problem. Later, Ernst et al. [26] developed new mathematical formulations for single and multiple allocation p-hub center problems and showed that their single allocation model is superior to the
CR IP T
one in [38]. They also proved that the problem is NP-hard. Yaman and Elloumi [66] introduced new models for star p-hub center problem to minimize the length of the longest O/D pairs in a star hub network. Contreras et al. [20] presented a generalized model for the classical p-center problem that combines facility location and network design problems by deciding simultaneously on the location of facilities and the design of its underlying network. Although integer programming optimization approaches are used to solve small hub problems, larger instances of the HLP are usually solved by heuristic or metaheuristic procedures. Campbell [12] presented a greedy-interchange
AN US
heuristic for the UMApHMP. Ernst and Krishnamoorthy [28] proposed a heuristic and an exact solution algorithm based on shortest paths for the same problem. Ebery et al. [25] developed an efficient shortest path based heuristic for the capacitated multiple allocation p-hub median problem. Jabalameli et al. [35] proposed an SA heuristic for the single allocation p-hub maximal covering problem. Hamacher and Meyer [33] proposed a binary search algorithm for solving the single allocation p-hub center problem which exploits the relationship between p-hub center and hub covering problems. An adaptive large neighborhood search (ALNS) algorithm was proposed by Chaharsooghi et al.
M
[15] for a reliable multiple allocation HLP with hub failure possibility. Silva and Cunha [59] presented an efficient tabu search (TS) algorithm for solving the uncapacitated single allocation p-hub maximal covering problem. More
ED
recently, Ghaffarinasab et al. [32] developed efficient simulated annealing based heuristics for solving the competitive single and multiple allocation hub location problems. For more details on the HLP and recent advances in this field,
PT
the interested readers are referred to surveys [4, 13, 30], and [19]. 2.2. Facility location with disruption and interdiction In the most of the classical facility location models, it is assumed that systems will operate as designed during its
CE
whole life. However, the network elements in real world environments are subject to various disruptions which can substantially affect the system’s performance. Disruptions, either caused by natural events or intentional attacks,
AC
may lead to unavailability of some elements in the service networks. Drezner [24] was the first researcher who studied location problems with unreliable facilities in p-median problems. Later, Drezner’s problem was tackled with different heuristic method by Lee [40]. Berman et al. [6] studied the location of unreliable facilities on a network where a customer may not be able to receive satisfactory service depending on the distance from the assigned facility. Synder and Daskin [61] considered the reliable p-median problem and the uncapacitated location problem assuming that the facilities are subject to random disruption with equal probabilities. Snyder et al. [62] conducted a review of different reliability considerations in facility location models. Berman et al. [7] introduced a reliable model based on the different failure probabilities for various facilities. In another work, Berman et al. [8] considered location of unreliable facilities where a customer have incomplete information about the state of operating facilities. Cui et
4
ACCEPTED MANUSCRIPT
al. [21] and Li and Ouyang [44] developed continuum approximation models to analyze the reliable facility location problems. Lim et al. [46] considered design of a reliable supply chain network in the presence of random facility disruptions with the option of hardening selected facilities. Peng et al. [54] addressed another reliable networks design problem considering a robust optimization measure in which the relative regret in each scenario is required to be no more than a constant value. Hong et al. [34] proposed two robust mathematical models to solve the emergency response facility location and transportation problems which are responsible for distributing relief goods after disaster. Lei and Tong [43] extended the reliable facility location model developed in [21] to introduce a level
CR IP T
assignment model for the p-median problem under non-uniform failure probabilities. Berman et al. [9] studied a class of location models with unreliable facilities and correlated failures with median and center objectives. More recently, Akg¨ un et al. [1] addressed a facility location problem in the context of disaster management in the pre-disaster phase of operations. They developed an optimization model that minimizes the risk that a demand point may be exposed to because it is not supported by the located facilities.
Unlike the above studies which focus on reliability under random disruption, there are a number of studies addressing
AN US
deliberate disruptions. Deliberate and rational attacks to the installed facilities was first addressed through the use of interdiction models for military applications [58], but later they were also used in business environments such as facility location and transportation network design. Church et al. [18] introduced the r-interdiction median problem (RIMP) and the r-interdiction covering problem (RICP) which were used to identify the set of facilities that, if lost, would harm service delivery the most. Later, Church et al. [57] considered the r-interdiction median problem with fortification (RIMPF) where the objective was to identify the most effective way of protecting a subset of existing
M
facilities so that the impact of the most disruptive attack to r unprotected facilities is minimized. Lei [41] introduced the r-interdiction center (RICent) problem in which the objective is to identify and interdict a
ED
fixed number of installed facilities in such a manner that the worst-case performance of the system (largest distance between demand points and the serving facilities) is further diminished. Aksen et al. [2] considered a budgetconstrained version of the RIMPF in which the facilities are assumed to have different protection costs and flexible
PT
capacities which can be expanded to accommodate the demand of customers previously assigned to interdicted facilities. Church and Scaparra [16] introduced a probabilistic version of RIMP where an attempted interdiction is
CE
successful only with a given probability. Losada et al. [47] further extended the model in [16] by assuming that the probability of a successful disrupting attack depends on its intensity and on the amount of disrupting resources used in the attack. Liberatore et al. [45] addressed a stochastic variant of RIMPF with uncertain number of interdictions.
AC
O’Hanley and Church [49] proposed a new version of RICP that maximizes a combination of initial coverage by p facilities and the minimum coverage level following the loss of the most critical r facilities. A new variation of the RIMP considering capacitated facilities and outsourcing option with partial interdiction was introduced by Aksen et al. [3] in which an attacked facility may preserve part of its capacity and the unmet demand after interdiction can be outsourced at some cost. According to the literature, only three papers have been published so far on hub location under intentional disruptions. Lei [42] introduced the hub interdiction median (HIM) problem which can be used to identify the set of critical facilities in a hub network. The author also proposed two bilevel formulations called the hub median protection problem (HMPP) and the HMPP with fortification (HMPP-F) where the former aims at locating hubs in anticipation 5
ACCEPTED MANUSCRIPT
of a worst-case facility loss and the latter tries to fortify a subset of existing hubs in a network to minimize postinterdiction system cost. Parvaresh et al. [51] formulated the multiple allocation p-hub median problem under intentional disruptions as a bilevel model in which the follower’s objective is to identify and interdict those hubs which would deteriorate the system’s efficiency in the worst possible way. For solving the problem, they proposed two algorithms based on SA and solved instances of up to 60 nodes using their heuristics. In a follow up work, Parvaresh et al. [52] extended their former model with two objective functions for the upper level problem where the leader aims at locating its hubs to minimize normal and worst-case transportation cost. They developed two
CR IP T
multi-objective metaheuristic algorithms based on SA and Tabu Search (TS) to solve this problem with instances of up to 50 nodes.
3. Mathematical Formulation
In this section, single level and bilevel mathematical models are presented for the three variants of the hub interdiction problem. In our proposed models, it is assumed that both the leader and the follower have complete
AN US
information about the game and will act rationally. Our models are based on the uncapacitated multiple allocation hub location problem under the assumption that the installed hubs are fully interconnected. Furthermore, we assume that transportation cost on the inter-hub connections are discounted by a fixed discount factor and that the direct connections between non-hub nodes are not allowed. The following notations are used hereafter: Sets and parameters
set of all demand nodes (i, j, k, m ∈ N )
H
set of nodes that are available for locating hubs H ⊆ N
wij
amount of flow originated at node i and destined to node j
dij
transportation cost (distance, time, etc.) of a unit flow from node i to node j
α
transfer cost per unit flow and per unit distance on inter-hub connections (0 ≤ α ≤ 1)
dijkm
cost of routing one unit of flow wij via hubs k and m, calculated as dijkm = dik + αdkm + dmj
p
number of hubs to be located by the leader
r
number of hubs to be interdicted by the follower
ED
PT
CE
β
M
N
covering radius
AC
Decision variables Xijkm Zk
Yk
fraction 1, = 0, 1, = 0,
of flow originated at node i and destined to node j routed via hubs k and m if node k is selected as a hub, otherwise if the hub located at node k is interdicted, otherwise
6
ACCEPTED MANUSCRIPT
Vij
=
Uij
=
1, if flow originated at node i and destined to node j is covered after interdiction, 0, otherwise 1, if the maximum O/D cost in the network is associated with node pair (i, j), 0, otherwise
CR IP T
3.1. r-hub interdiction median problem Assume that the leader is already operating in a market with hubs located at set Z ⊆ H of p nodes. The follower then attempts to attack r of the leader’s hubs (r ≤ p) in such a way that the leader’s operational cost is maximized. Define Fijkm as the set of existing hub pairs (q, s) for which the cost of routing one unit of flow between nodes i and j via hubs q and s is higher than the cost of routing the same flow via hubs k and m. In other words:
∀i, j ∈ N , k, m ∈ Z
AN US
Fijkm = {(q, s) ∈ Z × Z|diq + αdqs + dsj > dik + αdkm + dmj }
The follower’s problem then consists of the selection of a set of r installed hubs which will be disrupted so that total cost incurred by the leader is maximized. The MIP model for the multiple allocation r-hub interdiction median problem (r-HIMP) as an improved variant of the original HIM formulation presented by Lei [41, 42] can now be written as: XXX X
wij dijkm Xijkm
i∈N j∈N k∈Z m∈Z
s.t.:
X
Yk = r
X X
(1) (2)
ED
k∈Z
M
max
∀i, j ∈ N
(3)
∀i, j ∈ N , k ∈ Z
(4)
∀i, j ∈ N , k, m ∈ Z
(5)
Xijkm ≥ 0
∀i, j ∈ N , k, m ∈ Z
(6)
Yk ∈ {0, 1}
∀k ∈ Z
(7)
Xijkm = 1
k∈Z m∈Z
Xijkm +
X
PT
X
m∈Z
X
Xijmk ≤ 1 − Yk
Xijqs ≤ Yk + Ym
AC
CE
(q,s)∈Fijkm
m∈Z|m6=k
The objective function (1) maximizes the total transportation cost of the leader. Constraint (2) determines the number of hubs to be interdicted by the follower. Constraints (3) state that for each pair of nodes i and j, the whole flow wij must be routed via some hub pair. Constraints (4) ensure that flows can be routed via an installed hub if that hub is not interdicted. Constraints (5) are called the (multi-stage) closest assignment (CA) constraints which are generalized from the single-stage CA constraints originally introduced by Wagner and Falkson (W-F) [65]. They ensure that each O/D flow must be routed via a least-cost pair of hubs that are not disrupted. To be more precise, these CA constraints imply that for a pair of hubs (k, m), if both k and m are kept open (i.e., Yk + Ym = 0), the O/D flow wij cannot be routed via a pair of hubs which are more expensive than routing via k and m. Finally, constraints 7
ACCEPTED MANUSCRIPT
(6) and (7) are non-negativity and binary constraints, respectively. Now, let us assume that the leader is aware of the subsequent attack by the follower, and thus tries to locate its p hubs in such a way that the total transportation cost after the upcoming attack by the follower (including r interdictions) be minimum. The bilevel model for the hub interdiction median problem (BHIMP) can be written as: min H1 (Z) X s.t.: Zk = p
(8) (9)
k∈H
∀k ∈ H
where H1 (Z) = max
XXX X
wij dijkm Xijkm
i∈N j∈N k∈Z m∈Z
Yk ≤ Zk
AN US
s.t.: (2) − (7)
CR IP T
Zk ∈ {0, 1}
∀k ∈ Z
(10)
(11)
(12)
The objective function (8) minimizes the highest possible level of transportation cost after disruption of r hubs which is calculated by (11) as the objective of the lower level follower’s problem. Constraint (9) enforces the number of installed hubs by the leader to be p. Constraints (12) state that the follower can only interdict the hubs that have already been opened by the leader.
M
3.2. r-hub interdiction maximal covering problem
In the maximal covering variant of the problem, it is assumed that the leader’s objective is to maximize the total
ED
covered flows in the network. The O/D flow wij is said to be covered if it is routed via a pair of hubs in such a way that corresponding unit routing cost is no more than a prespecified covering radius (β). The follower aims to hit r of the installed hubs in such a manner that the total covered flows by the leader after interdiction is minimized. To
PT
model the problem, for each node pair (i, j), we define Mij as the set of hub pairs (k, m) through which the cost of routing one unit of flow wij is less than β. In other words:
CE
Mij = {(k, m) ∈ Z × Z|dik + αdkm + dmj ≤ β}
∀i, j ∈ N
Based on the above discussion, the multiple allocation r-hub interdiction maximal covering problem (r-HIMCP)
AC
solves the following MIP model: min
XX
wij Vij
(13)
i∈N j∈N
s.t.:(2), (7) Vij ≥ 1 − Yk − Ym
∀i, j ∈ N , (k, m) ∈ Mij
(14)
Vij ∈ {0, 1}
∀i, j ∈ N
(15)
The objective function (13) minimizes the total flows covered in the network. Constraints (14) stipulate that an O/D flow wij must be covered unless all the hub pairs (k, m) that currently cover it (i.e., (k, m) ∈ Mij ) are interdicted. 8
ACCEPTED MANUSCRIPT
We should note that the maximal covering version of the hub interdiction problem implicitly enforces CA conditions. In other words, in any optimal solution of this problem, all O/D flows are routed through a minimum cost path. Hence, there is no need to explicitly add CA constraints to this formulation. The bilevel model for the hub interdiction maximal covering problem (BHIMCP) can be proposed as: max H2 (Z)
(16)
where H2 (Z) = min
XX
wij Vij
i∈N j∈N
s.t.: (2), (7), (12), (14), (15)
CR IP T
s.t.:(9), (10)
(17)
The objective function (16) maximizes the lowest possible level of total covered flows after disruption of r hubs which
AN US
is calculated by (17) as the objective of the lower level follower’s problem. 3.3. r-hub interdiction center problem
Hub center problems are important for the delivery of perishable or time sensitive items. While in the median and maximal covering versions of the problem, the average-case performance of the system is diminished by the attacker, in the hub center interdiction problem, the follower plans to worsen the worst-case performance. In other
M
words, assuming a set of p existing hubs for the leader, the follower aims to select and disrupt r hubs so that highest transportation between all O/D pairs is maximized. The MIP model for the multiple allocation r-hub interdiction
ED
center problem (r-HICP) can then be written as:
PT
max
max
(i,j)∈N ×N
{
X X
k∈Z m∈Z
dijkm Xijkm }
(18)
s.t.:(2) − (7)
The objective function (18) maximizes the maximum transportation cost between all O/D pairs. Note that the above
CE
model is a nonlinear problem of maximax type which can be linearized by defining new set of variables. Let γ be a nonnegative variable acting as a surrogate for the maximum of O/D costs in the network. The MIP model for the
AC
multiple allocation r-hub interdiction center problem (r-HICP) can now be written as: max γ
(19)
s.t.:(2) − (7) X X γ≤ dijkm Xijkm + (1 − Uij )M k∈Z m∈Z
XX
∀i, j ∈ N
Uij = 1
(20) (21)
i∈N j∈N
γ≥0
(22)
Uij ∈ {0, 1}
∀i, j ∈ N 9
(23)
ACCEPTED MANUSCRIPT
where M is a sufficiently large number. An appropriate value for this number which results in a tight formulation is M = (2 + α)dmax where dmax = maxi,j∈N {dij }. The bilevel model for the hub interdiction center problem (BHICP) can be written as: min H3 (Z)
(24)
s.t.:(9), (10)
H3 (Z) = max γ s.t.: (2) − (7), (12), (20) − (23)
CR IP T
where (25)
The objective function (24) minimizes the highest transportation cost between all O/D pairs after disruption of r hubs which is calculated by (25) as the objective of the lower level follower’s problem.
Bilevel models are hard to solve even for a small number of decision variables [5], [22]. Therefore, to solve the
AN US
proposed models, we have developed efficient metaheuristic solution algorithms that are going to be discussed in detail in Section 4. 3.4. Alternative CA constraints
As stated earlier, CA constraints ensure that each O/D flow is routed through a least-cost pair of hubs that are not
M
disrupted. So far, different types of (single-stage) CA constraints have been introduced for discrete facility location problems (see [31, 41, 29] for comparison between various CA constraints). These single-stage CA constraints can be generalized as multi-stage CA constraints to be used in HLP models. For instance, Lei [41, 42] extended the
ED
single-stage CA constraints formulated by Church and Cohon (C-C) [17] and used it in his hub interdiction median (HIM) model as follows:
Xijqs + Xijkm ≥ 1 − Yk − Ym
PT
X
(q,s)∈Cijkm
∀i, j ∈ N , k, m ∈ Z
(26)
CE
where for all i, j ∈ N and k, m ∈ Z, the set Cijkm is defined as follows: Cijkm = {(q, s) ∈ Z × Z|diq + αdqs + dsj < dik + αdkm + dmj or
AC
diq + αdqs + dsj = dik + αdkm + dmj and (q < k or q = k and s < m)}
Constraints (26) maintain that if both hubs k and m are not disrupted (i.e., 1 − Yk − Ym = 1) and if no route has a P lower cost than route i-k-m-j (i.e., (q,s)∈Cijkm Xijqs = 0 ), then the route i-k-m-j must be assigned or chosen (i.e., Xijkm = 1).
Another type of CA constraints is proposed by Rojeski and ReVelle (R-R) [55] which has been extended to its multi-stage counterpart by Lei [41] as: X
(q,s)∈Cijkm
Qqs + Xijkm ≥ 1 − Yk − Ym
10
∀i, j ∈ N , k, m ∈ Z
(27)
ACCEPTED MANUSCRIPT
where Qkm = 1 if both hubs k and m are kept open. In other words, this variable can be defined as the product Qkm = (1 − Yk )(1 − Ym ). However, in order to linearize this equality, the following constraints are added along with (27): ∀k, m ∈ Z
(28)
Qkm ≤ 1 − Yk
∀k, m ∈ Z
(29)
Qkm ≤ 1 − Ym
∀k, m ∈ Z
(30)
Qkm ≥ 0
∀k, m ∈ Z
(31)
CR IP T
Qkm ≥ 1 − Yk − Ym
The logic behind the R-R constraints is as follows: If none of the hubs k and m are interdicted (i.e., 1 − Yk − Ym = 1), and there is no other open hub pair providing a lower routing cost than that of the hub pair (k, m) to the O/D flow P wij (i.e., (q,s)∈Cijkm Qqs = 0), then the flow wij must be routed via the hubs k and m (i.e., Xijkm = 1). As a further example, we extend the single-stage CA constraints originally introduced by Dobson and Karmarkar
AN US
(D-K) [23] to its multi-stage equivalent as: Qqs + Xijkm ≤ 1
∀i, j ∈ N , k, m ∈ Z, q, s ∈ Cijkm
(32)
along with the constraints (28)-(31). The idea of the D-K constraints is that the O/D flow wij is not routed via the hubs k and m (i.e., Xijkm = 0) if there is another pair of open hubs which provide a shorter path than that of the hubs k and m (i.e., 1 − Qqs = 0).
M
In Section 5, we will conduct a set of computational experiments in order to compare the performance of the four presented CA constraints, namely the W-F, C-C, R-R and D-K constraints in terms of solution time by using a
ED
commercial solver.
4. Metaheuristic Solution Algorithm
PT
In this section, we describe in detail the proposed simulated annealing (SA) metaheuristic algorithm for solving the six problems presented in Section 3. Inspired by heat treatment of metals in metallurgical industries, SA is
CE
a metaheuristic optimization algorithm which is effective in solving combinational optimization problems. SA was ˇ developed in 1953 by Metropolis et al. [48] and was independently described by Kirkpatrick et al. [39] and Cern` y [14]. The algorithm adopts the general concepts in gradual cooling process of metals for solving optimization problems.
AC
To this end, the SA algorithm starts with an initial solution at a high temperature and consecutively moves to the neighboring solutions via algorithm loops letting the temperature be reduced through the algorithm iterations. At each step, if the newly generated neighboring solution is better, in terms of the objective value, than the current solution, the current solution is replaced by the new solution. Otherwise, the algorithm accepts the new solution with a probability e−|∆E|/T where ∆E is the difference between the objective function values of the current and the new solutions and T is the current temperature. In early stages where the temperature is too high, there is a high probability to accept poor solutions. However, with a gradual decrease in temperature, there will be less probability to accept a bad solution in final stages. The algorithm stops if either the final temperature limit is reached or there is no improvements in a number of consecutive algorithm iterations. 11
ACCEPTED MANUSCRIPT
In the remainder of this section, we discuss several features of our proposed SA algorithms such as solution representation scheme, initial solution generation procedure, neighbor generation mechanisms, parameters used in the algorithms, and pseudo-codes for the proposed algorithms. 4.1. Solution Representation One of the important components of any metaheuristic algorithm is the solution representation. It has to be designed in such a way that facilitates neighbor generation and fast calculation of the objective functions. It must also
CR IP T
guarantee accessibility to the entire solution space [56]. In this paper, solutions are represented by one-dimensional arrays that contain numbers associated with the opened and operational hubs. For the leader’s problem, the array’s size is p comprising the set of hubs Z ⊆ H opened by the leader, whereas for the follower’s problem the representation array is of size r containing a subset Y of opened hubs that are interdicted by the follower (Y ⊆ Z). Note that having known the set of opened and operational (not interdicted) hubs, for each O/D pair (i, j), one can easily determine the paths for routing the associated flow wij by solving a shortest path problem or using enumeration.
AN US
4.2. Initial Solution Generation
The initial solution can affect the performance of a metaheuristic algorithm. Indeed, the quality of the final solution obtained by the metaheuristic algorithms relies to some extent on the quality of the initial solution. To produce an initial solution for the leader’s problem, we randomly select a subset Z ⊆ H of size p to serve as hub nodes. In case of the follower’s problem, an initial solution is generated by randomly selecting a subset Y ⊆ Z
M
showing the r hubs that are disrupted by the follower. The routing decisions are then determined based on the set
4.3. Neighborhood Structures
ED
of opened and operational hubs.
We define and apply a single operator for generating neighboring solutions in our proposed SA algorithms. This operator is called “Swap” and is used to alter one of the hubs in the solution representation array. In case of the
PT
leader’s problem, we first select a hub node and a non-hub node, randomly. Afterwards, the selected hub node becomes non-hub and the selected non-hub node becomes hub. For the follower’s problem, on the other hand, a
CE
disrupted hub and a non-disrupted hub are selected randomly and are changed afterwards. 4.4. Parameters used in the SA procedure
AC
Five parameters, namely T0 , Tf , δ1 , δ2 , and Imax are used in the proposed SA algorithms. T0 denotes the initial
temperature, whereas Tf is the final temperature below which the SA procedure is terminated. δ1 and δ2 are used as cooling rate parameters in the follower’s and the leader’s problems, respectively. Finally, Imax is the number of generated solutions at each temperature of the algorithm. 4.5. The overall SA Algorithms To solve the follower’s problem, it is assumed that the leader has already located his/her hubs on a set of candidate nodes Z ⊆ H. We start our algorithms with generating an initial solution Y ⊆ Z, as described previously. Then the follower’s objective value, f (Z, Y), can be calculated based on either (1), (13), or (18), depending on the variant of 12
ACCEPTED MANUSCRIPT
the HLP under consideration. We set the initial temperature as T0 which is an input parameter of the algorithm. Ybest denotes the best solution found so far and fbest denotes the corresponding objective function value. At each
temperature, we generate a new solution Y 0 by using the current solution based on the “Swap” operator presented
above. Subsequently, the objective function of the new solution f (Z, Y 0 ) is calculated. We define ∆E as the difference between the objectives of the new and current solutions, i.e., ∆E = f (Z, Y 0 ) − f (Z, Y). If ∆E > 0 for the r-HIMP
and r-HICP, or ∆E < 0 for the r-HIMCP, we update the current solution as Y ← Y 0 and if the objective f (Z, Y 0 )
of the new solution Y 0 is even better than best objective fbest so far, we also set fbest ← f (Z, Y 0 ) and Ybest ← Y 0 .
CR IP T
Otherwise, if ∆E ≤ 0 for the r-HIMP and r-HICP, or ∆E ≥ 0 for the r-HIMCP, we generate a random number ρ
from the interval [0,1]. If ρ is larger than e−|∆E|/T , we update the current solution as Y ← Y 0 . In other words, we
accept the solutions of worse quality with probability e−|∆E|/T to help the algorithm not get trapped in local optima. We repeat this cycle for Imax times at each temperature with the best solution obtained so far. Subsequently, we reduce the temperature at each iteration as T ← δ × T . The algorithm is terminated when the current temperature T drops below the prespecified final temperature Tf .
AN US
The SA algorithm for the leader’s problem is more complex than that of the follower’s problem as we aim at optimizing a bilevel model. In this case, we have a nested algorithm comprising an upper-level SA for optimizing the leader’s decisions and a lower-level SA that optimizes the follower’s decisions whenever a new solution is generated for the leader by the upper-level SA. We start the algorithm by generating an initial solution for the leader Z and improve this solution through the subsequent algorithm iterations. At each iteration of the upper-level SA where a new solution Z 0 is generated for the leader, the lower-level SA is run to obtain the corresponding best solution for the
M
follower Y ∗ (Z 0 ). In this case we define ∆E = f (Z 0 , Y ∗ (Z 0 )) − f (Z, Y ∗ (Z)). The new solution Z 0 is accepted if ∆E < 0 for the BHIMP and BHICP, or ∆E > 0 for the BHIMCP. Furthermore, if the objective f (Z 0 , Y ∗ (Z 0 )) of
ED
the new solution is even better than the best objective so far (fbest ), we set fbest ← f (Z 0 , Y ∗ (Z 0 )) and Zbest ← Z 0 .
The pseudo-codes of the proposed SA algorithms for the follower’s and the leader’s problems are illustrated in
PT
Algorithms 1 and 2, respectively.
5. Computational Experiments
CE
In this section, we present the results obtained from our computational experiments. Extensive experimentations are conducted in order to test the efficiency of the proposed mathematical models and the CA constraints as well as the SA algorithms. To this end, we use two famous data sets from the literature of the HLP: the CAB and the
AC
TR data sets. The CAB data set introduced by O’Kelly [50] is based on the airline passenger interactions between 25 US cities in 1970 evaluated by the Civil Aeronautics Board (CAB). This data set has been used by most of the hub location researchers in the literature in which all the 25 nodes are candidates for being hubs (|H| = 25). The second data set that is used in our computational experiments is the TR data set [63] which is based on the cargo flows between 81 cities of Turkey where only 22 of these cities are candidate nodes for locating hubs (|H| = 22). In P both the data sets, the O/D flows are scaled by the sum of flows (i.e., i,j wij ). The proposed SA algorithms are
implemented in Microsoft Visual C# 2013 (version 5.0). Moreover, the mathematical models r-HIMP, r-HIMCP, and r-HICP are solved independently using CPLEX version 12.6. All the experiments have been run on a computer
13
ACCEPTED MANUSCRIPT
Algorithm 1 : SA for the follower’s problem (T0 , Tf , δ1 , Imax , α, Z, r) 1: Generate a random initial solution Y 2: Calculate f (Z, Y) 3: T ← T0 , fbest ← f (Z, Y), Ybest ← Y, I ← 0 4: while T > Tf do
7:
Generate a new solution Y 0 based on Y using “Swap” operator Calculate f (Z, Y 0 )
8:
∆E ← f (Z, Y 0 ) − f (Z, Y)
9:
if ∆E > 0 (or ∆E < 0) then Y ← Y 0 , f (Z, Y) ← f (Z, Y 0 )
10: 11:
CR IP T
6:
for I < Imax do
else ρ ← rand(0, 1)
12:
if ρ > e−|∆E|/T then
13:
Y ← Y 0 , f (Z, Y) ← f (Z, Y 0 );
14:
end if
15: 16:
end if
17:
if f (Z, Y) > fbest (or f (Z, Y) < fbest ) then end if
20:
Y ← Ybest
21:
I ←I +1
22:
end for
23:
I←0
24:
T ← δ1 × T
ED
19:
M
Ybest ← Y, fbest ← f (Z, Y)
18:
PT
25: end while 26: return Ybest , fbest
AN US
5:
equipped with Intel(R) Core(TM) i3-3220 CPU of 3.30 GHz and 16 GB of RAM, using the Microsoft Windows 7
CE
operating system.
5.1. Test instances and parameters setting
AC
We have generated and used a large number of instances from the CAB and TR data sets. The number of hubs
to be opened by the leader (p) for the CAB data set is taken as 3, 4, or 5, whereas in the TR data set, p is chosen as 6, 10, or 14. The number of disrupted hubs by the follower are selected as r ∈ {0, 1, ..., (p − 1)} for the CAB instances and as r ∈ {0, 1, ..., 5} for the TR instances. The parameter α is considered at three levels: 0.4, 0.6, and 0.8 for both data sets. For solving the maximal covering variant, different levels for the covering radius parameter (β) are considered. For the CAB data set, β takes a value as 500, 1000, or 1500, whereas for the TR data set the covering radius is chosen as 400, 700, or 1000. The parameters of the proposed SA algorithms are tuned by setting a good trade-off between time and quality of
14
ACCEPTED MANUSCRIPT
Algorithm 2 : SA for the leader’s (bilevel) problem (T0 , Tf , δ2 , α, p, r) 1: Generate a random initial solution Z
2: Get Y ∗ (Z) and f (Z, Y ∗ (Z)) by solving the follower’s problem using Algorithm 1 3: T ← T0 , fbest ← f (Z, Y ∗ (Z)), Zbest ← Z 4: while T > Tf do
Get Y ∗ (Z 0 ) and f (Z 0 , Y ∗ (Z 0 )) by solving the follower’s problem using Algorithm 1
7:
∆E ← f (Z 0 , Y ∗ (Z 0 )) − f (Z, Y ∗ (Z))
8:
if ∆E < 0 (or ∆E > 0) then
9: 10: 11: 12:
Z ← Z 0 , f (Z, Y ∗ (Z)) ← f (Z 0 , Y ∗ (Z 0 )) else ρ ← rand(0, 1)
if ρ > e−|∆E|/T then Z ← Z 0 , f (Z, Y ∗ (Z)) ← f (Z 0 , Y ∗ (Z 0 ))
13: 14:
CR IP T
6:
Generate a new solution Z 0 based on Z using “Swap” operator
end if
AN US
5:
15:
end if
16:
if f (Z, Y ∗ (Z)) < fbest (or f (Z, Y ∗ (Z)) > fbest ) then
17:
Zbest ← Z, fbest ← f (Z, Y ∗ (Z))
18:
end if
19:
T ← δ2 × T
20: end while
M
21: return Zbest , fbest
ED
the final solutions. In a set of preliminary experiments, different combinations of parameters were tested on a large number of test instances and the values reported in Table 1 have been selected as the best values which lead to
PT
high-quality solutions in short CPU times both for the leader’s and the follower’s problems under the two data sets.
AC
CE
Table 1: SA parameters used for different data sets Median
Max. Cov.
Center
CAB
TR
CAB
TR
CAB
TR
T0
100
200
100
200
100
200
Tf
1
1
1
1
1
1
δ1
0.985
0.985
0.995
0.995
0.995
0.999
δ2
0.99
0.99
0.992
0.992
0.992
0.997
Imax
4
4
4
4
4
4
5.2. Comparison of the performance of different CA constraints In the first set of our computational experiments, we try to test the performance of the four different multi-stage CA constraints introduced in Section 3, namely the W-F, C-C, R-R, and D-K constraints. To this end, 12 test instances with different p and r values have been selected from the CAB and TR data sets. For each instance, the r-hub interdiction median problem (r-HIMP) and the r-hub interdiction center problem (r-HICP) are solved using
15
ACCEPTED MANUSCRIPT
different CA constraints under three distinct values of α by CPLEX. The existing p hubs are assumed to be located based on the solutions of the p-hub median and the p-hub center problems, respectively. Since the model proposed for the r-hub interdiction maximal covering problem (r-HIMCP) does not need CA constraints, we have not included this problem in the current part of our experiments. The results obtained by solving the r-HIMP are reported in Table 2. The first column in this table represent the name of the data set, whereas the second and third columns show the number of the leader’s (existing) hubs and the number of hubs which are disrupted by the follower, respectively. The next five columns report the solution results for the case of α = 0.4. More specifically, the column labeled as
CR IP T
“OF” gives the optimal value of the objective function and the four columns under the heading “CPU(s)” show the solution times, in seconds, that are spent to obtain the optimal solutions by CPLEX using different CA constraints. Similarly, the results for the cases of α = 0.6 and 0.8 are reported in the next columns. We used a time limit of 10800 second (3 hours) for CPLEX in our experiments.
Table 2: CPLEX results with different CA constraints for the r-HIMP α = 0.4
5
TR
6
10
OF
α = 0.6
CPU(s)
OF
C-C
R-R
D-K
W-F
C-C
R-R
D-K
OF
CPU(s) R-R
D-K
1
1088.55
0.27
0.50
0.53
0.92
1140.78
0.13
0.30
0.30
0.55
1190.91
0.15
0.22
0.17
0.54
2
1718.25
0.72
0.84
0.42
0.64
1413.48
0.59
0.70
0.72
0.81
1422.42
0.80
0.90
0.42
0.80
3
3517.63
0.14
0.18
0.24
0.61
3517.63
0.14
0.20
0.13
0.59
3517.63
0.13
0.14
0.13
0.55
1
1014.74
0.51
0.43
0.74
1.44
1095.45
0.40
0.42
0.60
1.44
1159.49
0.20
0.21
0.39
1.09
2
1606.23
1.79
2.05
0.85
1.43
1656.07
1.82
1.96
0.80
1.11
1623.38
1.10
1.14
0.53
1.07
3
1903.06
2.15
2.09
1.10
2.08
1946.07
1.68
1.88
1.37
1.53
1930.37
1.89
1.62
0.75
2.70
1
666.11
8.91
15.18
28.32
48.42
734.69
5.77
14.23
33.53
49.62
786.67
5.30
12.06
34.89
53.11
2
860.75
17.83
169.44
72.85
1869.90
867.44
14.76
125.27
88.12
502.91
899.02
10.65
111.13
70.24
565.70
3
1147.94
15.52
345.21
128.52
1656.50
1080.74
23.57
203.95
105.91
890.22
1096.33
13.33
194.80
91.79
500.18
1
538.78
33.11
78.90
278.13
2441.04
636.98
30.79
78.10
303.77
1301.61
710.74
28.89
64.61
328.04
1042.51
2
605.50
557.97
2132.21
530.51
> 10800
675.49
454.60
2340.90
466.20
> 10800
742.41
289.61
1813.45
424.89
> 10800
3
701.95
709.52
2941.45
917.88
> 10800
725.17
524.10
3896.84
814.43
> 10800
779.72
397.02
2628.02
582.44
> 10800
112.37
474.04
163.34
> 2301.92
88.20
555.40
151.32
> 2029.20
62.42
402.36
127.89
> 1980.69
Avg.
W-F
α = 0.8
CPU(s)
AN US
4
r
M
CAB
p
ED
Data set
W-F
C-C
PT
Note from Table 2 that the solution time for the model with W-F constraints is the smallest among the four CA constraints. The average times for solving the instances using the W-F constraints are less than two minutes for the
CE
three values of α. The second best performance among the four CA constraints belongs to the R-R constraints with the average solution times of less than three minutes. The third best performance in terms of solution times belongs to the C-C constraints, whereas the weakest performance is associated with the D-K constraints. Observe that the
AC
model with the D-K constraints was not able to solve the two largest instances within the time limit of three hours. Note also that, for the instances of the CAB data set, the four CA conditions have more or less similar solution times. However, for larger instances, as in the TR data set, the differences between the solution times associated with different CA constraints become more recognizable. Table 3 shows the results obtained by CPLEX for solving the r-HICP with different CA constraints.
16
ACCEPTED MANUSCRIPT
Table 3: CPLEX results with different CA constraints for the r-HICP α = 0.4
CAB
p
4
5
TR
6
10
r
OF
α = 0.6
CPU(s)
OF
C-C
R-R
D-K
OF
W-F
C-C
R-R
D-K
CPU(s) R-R
D-K
1
2736.53
0.47
0.78
0.77
1.10
2752.40
0.26
0.62
0.76
0.92
4915.25
0.25
0.53
0.68
0.81
2
5007.65
1.28
1.19
0.59
1.27
5007.65
0.90
0.92
0.57
1.06
5053.12
0.83
0.77
0.58
0.82
3
5451.58
0.53
0.54
0.86
1.23
5451.58
0.77
0.36
0.34
0.67
5406.80
0.59
0.28
0.32
0.74
1
2390.81
0.58
0.64
1.24
1.49
2583.17
0.33
0.48
0.60
1.34
2824.84
0.29
0.36
0.54
1.41
2
3012.90
1.55
1.99
1.29
1.97
4598.85
1.07
2.13
1.29
1.81
4598.85
0.78
1.78
2.61
2.77
3
5007.65
1.01
8.51
1.62
2.81
5229.61
1.64
7.76
3.76
2.96
5229.61
1.24
3.10
2.01
2.34
1
1522.20
20.03
43.15
35.53
60.13
1728.20
14.19
26.28
33.97
57.25
1876.60
12.39
15.69
34.83
61.12
2
1834.20
205.84
543.74
140.39
10065.01
2282.00
65.10
352.25
86.07
2757.83
2282.00
32.07
493.28
87.77
5674.85
3
3014.00
126.58
1205.43
513.81
1438.85
3336.00
67.39
519.78
639.19
1
1317.00
71.19
127.82
305.90
888.34
1610.40
69.47
119.12
414.17
2
1514.80
1063.41
4731.96
793.51
> 10800
1719.00
854.86
4961.58
1318.31
3
2282.00
1286.67
> 10800
> 10800
> 10800
1856.00
909.92
> 10800
> 10800
231.60
> 1455.48
> 1049.63
> 2838.52
165.49
> 1399.27
> 1108.25
Avg.
W-F
α = 0.8
CPU(s)
W-F
C-C
CR IP T
Data set
1201.07
3032.00
46.80
320.13
112.56
1126.72
1870.60
45.06
182.41
365.71
950.48
> 10800
1915.20
260.55
6410.62
1221.65
> 10800
> 10800
1934.20
1184.03
> 10800
> 10800
> 10800
132.07
> 1519.08
> 1052.44
> 2391.83
> 2229.30
406.57
It can be seen from Table 3 that the solution times for the r-HICP are, on average, larger than the corresponding
AN US
times for the r-HIMP. Furthermore, one can observe that, similar to the case of median problem, the W-F constraints perform better than the other three types of multi-stage CA constraints for the center problem. The R-R, C-C, and D-K constraints are placed in the second, the third, and the fourth ranks in terms of the average solution time. The results also show that, it is possible for CPLEX to solve all the instances within the three hour time limit only by using W-F constraints. For all the other three CA constraints, at least one of the instances are not solved within this time limit. Based on the reported results in Tables 2 and 3, we use the W-F constraints for solving problems by
M
the commercial solver in the rest of our computational experiments. 5.3. Results for hub interdiction median problems
ED
Results obtained by solving the r-HIMP on the CAB data set are presented in Table 4. In case of the follower’s problem, it is assumed that the leader has already located its hubs on a subset Z ⊆ H of nodes. To this end, we
PT
assume that the leader is unaware of the upcoming attack by the follower, and hence selects its p hubs based on the solution of the uncapacitated multiple allocation p-hub median problem (UMApHMP) aiming to minimize its
AC
CE
transportation cost.
17
ACCEPTED MANUSCRIPT
Table 4: Results for the r-HIMP with the CAB data set α = 0.4 p
r
CPLEX
α = 0.6 SA
CPLEX
α = 0.8 SA
CPLEX
SA
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
0
859.63
0.21
859.63
0.02
949.23
0.11
949.23
0.02
1020.03
0.09
1020.03
0.02
1
1206.33
0.23
1206.33
0.02
1263.61
0.12
1263.61
0.02
1309.29
0.12
1309.29
0.02
2
3517.63
0.08
3517.63
0.01
3517.63
0.07
3517.63
0.01
3517.63
0.06
3517.63
0.01
0
754.48
0.27
754.48
0.01
866.44
0.16
866.44
0.02
951.75
0.07
951.75
0.02
1
1088.55
0.27
1088.55
0.02
1140.78
0.13
1140.78
0.02
1190.91
0.15
1190.91
0.02
2
1718.25
0.72
1718.25
0.02
1413.48
0.59
1413.48
0.02
1422.42
0.80
1422.42
0.02
3
3517.63
0.14
3517.63
0.01
3517.63
0.14
3517.63
0.01
3517.63
0.13
3517.63
0.01
0
676.33
0.37
676.33
0.02
804.70
0.28
804.70
0.02
910.35
0.13
910.35
0.02
1
1014.74
0.51
1014.74
0.03
1095.45
0.40
1095.45
0.04
1159.49
0.20
1159.49
0.04
2
1606.23
1.79
1606.23
0.02
1656.07
1.82
1656.07
0.02
1623.38
1.10
1623.38
0.02
3
1903.06
2.15
1903.06
0.02
1946.07
1.68
1946.07
0.02
1930.37
1.89
1930.37
0.02
4
3517.63
0.39
3517.63
0.02
3517.63
0.24
3517.63
0.02
3517.63
0.23
3517.63
0.02
Avg.
1781.71
0.59
1781.71
0.01
1807.39
0.47
1807.39
0.02
1839.24
0.41
1839.24
0.02
3
4
5
CR IP T
OF
AN US
Note that for each value of p, when the value of r equals 0, the leader’s cost is the same as the optimal cost of UMApHMP, as the follower has not interdicted any hubs. Observe also that, as the value of r gets larger, the leader’s cost increases substantially. This can be interpreted as a result of the fact that the leader has decided on the location of its hubs without being aware of upcoming disruptions. Furthermore, the results show that higher values for the discount factor (α) result in slightly elevated costs for the leader as the inter-hub transportation cost
M
grow with the value of the discount factor. It can be seen from the above table that the proposed SA algorithm has obtained the optimal solution for all instances within a negligible CPU time, indicating the high efficiency of the SA
model (1)-(7) is also efficient.
ED
algorithm. Observe also that, the average solution times for CPLEX are less than a second that means the proposed
Table 5 shows the results for solving the r-HIMP on the TR data set. Here, it is also assumed that, the leader has
AC
CE
PT
already selected its p hubs based on the solution of the UMApHMP.
18
ACCEPTED MANUSCRIPT
Table 5: Results for the r-HIMP with the TR data set α = 0.4
10
14
α = 0.6
CPLEX
SA
α = 0.8
CPLEX
SA
CPLEX
SA
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
0
558.70
2.76
558.70
0.10
650.42
1.56
650.42
0.10
725.25
1.70
725.25
0.10
1
666.11
8.91
666.11
0.65
734.69
5.77
734.69
0.68
786.67
5.30
786.67
0.66
2
860.75
17.83
860.75
0.46
867.44
14.76
867.44
0.47
899.02
10.65
899.02
0.47
3
1147.94
15.52
1147.94
0.31
1080.74
23.57
1080.74
0.32
1096.33
13.33
1096.33
0.32
4
1557.13
13.62
1557.13
0.20
1235.56
23.39
1235.56
0.19
1244.05
18.21
1244.05
0.20
5
1783.13
14.94
1783.13
0.12
1783.13
13.71
1783.13
0.12
1783.13
11.86
1783.13
0.13
0
475.45
7.81
475.45
0.10
581.77
7.91
581.77
0.10
675.96
8.87
675.96
0.10
1
538.78
33.11
538.78
1.88
636.98
30.79
636.98
2.08
710.74
28.89
710.74
1.93
2
605.50
557.97
605.50
1.57
675.49
454.60
675.49
1.70
742.41
289.61
742.41
1.60
3
701.95
709.52
701.95
1.13
725.17
524.10
725.17
1.19
779.72
397.02
779.72
1.14
4
806.50
871.29
806.50
0.95
807.77
629.98
807.77
0.99
847.52
533.12
847.52
0.95
5
1019.68
1109.26
1019.68
0.64
992.61
969.82
992.61
0.67
1015.33
555.11
1015.33
0.66
0
427.26
23.06
427.26
0.10
545.55
22.32
545.55
0.10
653.22
21.53
653.22
0.10
1
490.97
124.32
490.97
7.15
600.84
123.79
600.84
7.26
673.14
127.34
673.14
7.27
2
549.96
1605.60
549.96
6.14
639.51
1216.38
639.51
6.32
692.83
1060.78
692.83
6.51
3
600.79
2339.64
600.69
5.35
686.87
1650.53
686.87
5.39
722.71
1460.90
722.71
5.56
4
659.22
2919.40
659.22
4.42
732.95
1848.92
732.95
4.59
760.15
2121.35
760.15
4.91
5
749.46
3569.30
749.46
4.10
810.10
2880.21
810.10
3.88
796.59
2526.00
796.59
4.08
819.40
774.65
819.40
1.96
821.53
580.11
821.53
2.00
866.93
510.64
866.93
2.03
Avg.
CR IP T
6
r
AN US
p
Note from the above table that the SA algorithm has obtained the optimal solution for all instances of the TR
M
data set. From a solution time perspective, it is shown that CPLEX solves all the instances in reasonable times. However, the solution times for the proposed SA are too small as the SA solves the instances in around 2 seconds, on average, which shows how efficient the proposed SA algorithm is even for solving large instances.
ED
Tables 6 and 7 show the results obtained by solving the BHIMP using the proposed SA algorithm for the CAB and TR data sets, respectively. Since it is not possible to directly solve the bilevel models by using CPLEX, we have
PT
only reported the results from the SA in these tables. Table 6: SA results for the BHIMP with the CAB data set
AC
CE
p
r
α = 0.4
α = 0.6
α = 0.8
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
1
1206.33
17.91
1263.61
17.88
1301.60
16.59
2
1511.81
7.95
1511.81
7.84
1511.81
8.40
1
1048.01
25.75
1087.43
22.08
1131.99
20.86
2
1315.41
18.35
1401.93
18.09
1384.68
18.45
3
1536.93
9.11
1596.49
9.41
1523.65
7.97
1
891.74
37.50
976.68
38.35
1033.03
40.45
2
1358.86
22.74
1242.10
22.91
1385.17
23.50
3
1457.88
16.56
1427.66
16.10
1654.53
17.09
4
1878.22
8.12
1536.93
7.98
1718.52
8.05
Avg.
1356.13
18.22
1338.29
17.85
1405.00
17.93
3
4
5
Results presented in Tables 6 and 7 reveal that the leader can decrease its transportation cost by acting rationally
19
ACCEPTED MANUSCRIPT
Table 7: SA results for the BHIMP with the TR data set r
10
14
α = 0.6
α = 0.8
CPU(s)
OF
CPU(s)
OF
CPU(s)
1
654.13
715.45
725.29
673.13
779.73
642.17
2
753.66
703.63
809.98
690.95
852.80
687.23
3
1024.90
410.38
957.23
409.62
976.98
409.25
4
1029.26
240.87
1036.56
241.33
1083.37
239.09
5
1140.65
107.46
1140.65
112.28
1140.65
109.56
1
532.82
1957.38
628.90
1940.12
707.41
1917.71
2
583.93
1708.12
671.73
1607.70
737.82
1626.26
3
671.15
1309.55
710.76
1309.87
769.82
1298.30
4
741.80
960.43
786.07
958.97
832.30
961.11
5
947.13
701.81
902.83
695.09
945.90
695.99
1
474.90
3962.61
578.61
3828.11
673.14
3928.73
2
509.41
3347.16
609.43
3723.12
692.83
3288.31
3
570.18
2918.01
661.48
2879.43
717.09
2717.38
4
592.15
2209.65
697.69
2274.19
754.83
2281.49
5
682.30
1992.16
727.22
1549.64
Avg.
CR IP T
6
α = 0.4 OF
AN US
p
772.64
1985.50
765.38
1823.65
779.32
1555.29
828.67
1508.41
in anticipation of the follower’s attack. For example, as can be seen in Table 7, in case p = 6 and r = 5 with α = 0.8, if the leader locates its hubs using the bilevel BHIMP model, then the transportation cost will reduce considerably
M
from 1783.13 (when the leader locates its hubs based on the UMApHMP) to 1140.65. Note also that, the CPU times for solving the BHIMP are much longer than the CPU times needed to solve the r-HIMP. This is due to the bilevel nature of the BHIMP that, as mentioned in the previous section, requires the proposed SA algorithm to solve the
ED
follower’s problem whenever it obtains a new solution for the leader (nested SA algorithm). However, the times are still very good for a strategic planning problem such as locating facilities in a competitive environment.
PT
5.4. Results for hub interdiction maximal covering problems Tables 8 and 9 present the results for solving the r-HIMCP with the CAB and TR data sets, respectively. In
CE
this case, it is assumed that the leader has already selected its p hubs based on the solution of the uncapacitated multiple allocation p-hub maximal covering problem (UMApHMCP) with the aim of maximizing the total covered flows in the network. As stated earlier, in this problem, the follower tries to minimize the leader’s total covered flows
AC
by interdicting r hubs from among the installed p hubs. To show the covering radius considered for each instance, a new column labeled as “β” is added in the corresponding tables. The columns labeled as “OF (%)” show the P total covered flows by the leader as percentage of the total flows within the network ( i,j wij ). As can be seen, the SA algorithm has again obtained the optimal solutions for all instances in less than one and three seconds, on
average, for the CAB and TR data sets, respectively. It should also be noted that the CPU times for CPLEX are also very short even for the large TR instances which shows the efficiency of the proposed mathematical model for the r-HIMCP. Observe from Tables 8 and 9 that, in case the leader decides on its hubs based on UMApHMCP, the follower can dramatically shrink the total value of covered flows by interdicting some r hubs. Moreover, by increasing the value 20
ACCEPTED MANUSCRIPT
Table 8: Results for the r-HIMCP with the CAB data set α = 0.4
5
3
4
5
3
4
5
SA
OF (%)
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
0
37.05
0.20
37.05
0.05
34.87
0.12
34.87
0.05
28.41
0.02
28.41
0.05
1
14.71
0.05
14.71
0.05
16.03
0.02
16.03
0.05
14.64
0.01
14.64
0.05
2
6.70
0.03
6.70
0.05
6.70
0.07
6.70
0.04
6.70
0.01
6.70
0.05
0
42.90
0.24
42.90
0.05
38.17
0.06
38.17
0.05
31.31
0.03
31.31
0.05
1
16.59
0.11
16.59
0.06
15.55
0.06
15.55
0.07
17.54
0.03
17.54
0.07
2
8.21
0.15
8.21
0.06
8.17
0.07
8.17
0.06
9.60
0.06
9.60
0.06
3
0.48
0.02
0.48
0.05
1.00
0.02
1.00
0.05
2.90
0.02
2.90
0.05
0
47.08
0.33
47.08
0.05
41.07
0.05
41.07
0.05
32.91
0.05
32.91
0.05
1
19.08
0.10
19.08
0.08
18.45
0.09
18.45
0.08
19.14
0.19
19.14
0.09
2
10.43
0.08
10.43
0.08
11.08
0.18
11.08
0.12
11.08
0.05
11.08
0.11
3
2.02
0.04
2.02
0.09
3.90
0.19
3.90
0.07
3.90
0.22
3.90
0.07
4
0.48
0.05
0.48
0.05
1.00
0.05
1.00
0.05
1.00
0.33
1.00
0.05
17.14
0.11
17.14
0.06
16.33
0.08
16.33
0.06
14.92
0.08
14.92
0.06
0
65.90
0.07
65.90
0.05
61.90
0.02
61.90
0.05
57.49
0.02
57.49
0.05
1
16.81
0.07
16.81
0.05
22.51
0.02
22.51
0.05
20.10
0.18
20.10
0.06
2
3.56
0.02
3.56
0.05
3.56
0.01
3.56
0.05
2.71
0.01
2.71
0.05
0
74.44
0.12
74.44
0.05
67.01
0.03
67.01
0.05
62.14
0.02
62.14
0.05
1
30.82
0.04
30.82
0.06
49.28
0.05
49.28
0.07
24.75
0.04
24.75
0.06
2
8.30
0.03
8.30
0.07
9.11
0.27
9.11
0.06
7.37
0.41
7.37
0.06
3
3.56
0.02
3.56
0.05
3.56
0.02
3.56
0.05
2.71
0.02
2.71
0.05
0
82.04
0.26
82.04
0.05
72.71
0.08
72.71
0.05
66.57
0.04
66.57
0.05
1
48.34
0.04
48.34
0.08
54.98
0.05
54.98
0.08
40.54
0.09
40.54
0.09
2
15.95
0.08
15.95
0.08
14.81
0.05
14.81
0.12
13.58
0.14
13.58
0.12
3
8.30
0.07
8.30
0.09
8.21
0.06
8.21
0.07
7.37
0.06
7.37
0.07
4
3.56
0.08
3.56
0.05
3.56
0.08
3.56
0.05
2.71
0.08
2.71
0.05
30.13
0.07
30.13
0.06
30.93
0.06
30.93
0.06
25.67
0.09
25.67
0.06
0
94.47
0.05
94.47
0.05
84.55
0.03
84.55
0.05
81.31
0.02
81.31
0.05
1
57.60
0.16
57.60
0.05
59.76
0.02
59.76
0.05
38.08
0.02
38.08
0.06
2
7.57
0.02
7.57
0.05
7.57
0.02
7.57
0.04
7.57
0.01
7.57
0.05
0
96.96
0.08
96.96
0.05
88.51
0.03
88.51
0.05
83.86
0.03
83.86
0.05
1
69.64
0.04
69.64
0.06
71.65
0.07
71.65
0.07
52.75
0.04
52.75
0.07
2
29.43
0.03
29.43
0.06
26.04
0.05
26.04
0.06
26.77
0.05
26.77
0.06
3
7.57
0.03
7.57
0.05
7.57
0.03
7.57
0.05
7.57
0.06
7.57
0.04
0 1 2
99.08
0.06
99.08
0.05
90.18
0.05
90.18
0.05
84.64
0.06
84.64
0.05
71.21
0.11
71.21
0.08
73.32
0.05
73.32
0.08
54.31
0.08
54.31
0.08
30.18
0.05
30.18
0.10
26.63
0.08
26.63
0.12
38.50
0.07
38.50
0.12
7.90
0.25
7.90
0.08
7.90
0.33
7.90
0.07
23.02
0.16
23.02
0.07
2.00
0.05
2.00
0.05
6.30
0.05
6.30
0.05
7.57
0.04
7.57
0.05
47.80
0.07
47.80
0.06
45.83
0.06
45.83
0.06
42.16
0.05
42.16
0.06
CE
3 4
Avg.
CPLEX
CPU(s)
Avg. 1500
α = 0.8 SA
OF (%)
Avg. 1000
CPLEX
CR IP T
4
α = 0.6 SA
AN US
3
CPLEX
M
500
r
ED
p
PT
β
AC
of r, the reduction in the coverage percentage gets more considerable, as expected. For example, in the CAB data set with p = 5, α = 0.4, β = 1500, and r = 0, if the leader chooses its hubs according to optimal solution of UMApHMCP, the flow coverage percentage will be 99.08%. However, as the follower decides to disrupt 1, 2, 3, or 4 hubs, the total covered flows by the leader will fall to 71.21%, 30.18%, 7.90%, or 2.00%, respectively. Another interesting observation is that, for smaller covering radius (β) values, the follower’s attack can be more effective with the same number of interdictions. For instance, consider the case p = 14, r = 5, and α = 0.4 in the TR data set. In case the covering radius is 700, the five interdictions drop the leader’s coverage from 92.48% to 35.24%, whereas for β = 1000, the leader’s coverage before and after the interdiction are 99.96% and 79.10%, respectively. Furthermore,
21
ACCEPTED MANUSCRIPT
the results show that as the discount factor value (α) grows, the coverage percentage falls down as a result of elevated inter-hub transportation cost. Table 9: Results for the r-HIMCP with the TR data set α = 0.4
10
14
6
10
10
14
Avg.
SA
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
0
30.04
0.91
30.04
1.00
23.79
0.61
23.79
1.00
17.89
0.79
17.89
1.00
1
15.35
0.95
15.35
1.32
15.90
0.50
15.90
1.32
12.93
0.43
12.93
1.31
2
10.21
0.49
10.21
0.95
10.78
0.59
10.78
0.90
9.03
0.50
9.03
0.93
3
6.27
0.46
6.27
0.63
6.20
0.65
6.20
0.63
5.69
0.49
5.69
0.71
4
3.79
0.73
3.79
0.41
3.54
0.48
3.54
0.39
3.66
0.71
3.66
0.39
5
1.71
0.57
1.71
0.28
1.71
0.60
1.71
0.26
1.79
0.48
1.79
0.26
0
40.93
1.83
40.93
1.00
29.76
1.36
29.76
1.00
21.74
1.98
21.74
1.00
1
30.87
3.61
30.87
3.72
20.92
1.62
20.92
3.68
17.01
1.59
17.01
3.71
2
20.48
2.08
20.48
3.15
15.98
1.71
15.98
3.17
12.67
1.62
12.67
3.11
3
14.47
1.84
14.47
2.26
11.62
1.85
11.62
2.28
9.68
6.16
9.68
2.25
4
10.36
1.63
10.36
1.96
9.41
1.59
9.41
1.92
7.64
1.63
7.64
1.92
5
6.77
7.20
6.77
1.35
7.53
2.03
7.53
1.39
6.17
1.56
6.17
1.30
0
48.42
3.33
48.42
1.00
33.85
2.99
33.85
1.00
24.03
2.93
24.03
1.00
1
39.33
3.13
39.33
7.09
27.92
2.43
27.92
7.14
20.05
2.38
20.05
7.22
2
32.08
2.49
32.08
6.09
22.10
8.43
22.10
6.11
17.17
2.44
17.17
6.20
3
25.87
2.69
25.87
5.33
4
20.83
2.89
20.83
4.46
5
15.99
4.11
15.99
3.81
20.76
2.27
20.76
2.54
0
73.25
0.66
73.25
1.00
1
52.17
0.70
52.17
1.38
2
31.48
0.73
31.48
0.98
3
19.20
0.74
19.20
0.63
4
9.74
0.87
9.74
0.43
5
3.70
0.71
3.70
0.27
0
86.89
1.70
86.89
1
66.73
2.40
66.73
2
48.18
6.71
48.18
3
33.98
2.47
33.98
4
26.42
5
19.57
0
92.48
1
83.62
2
70.11
18.46
2.46
18.46
5.47
14.39
3.61
14.39
5.22
14.52
2.42
14.52
4.46
11.44
4.93
11.44
4.51
12.20
4.32
12.20
3.76
9.97
4.56
9.97
3.77
15.89
2.03
15.89
2.54
12.53
2.15
12.38
2.54
60.98
0.54
60.98
1.00
52.19
0.47
52.19
1.10
47.40
0.61
47.40
1.37
42.16
0.83
42.16
1.31
34.97
0.70
34.97
0.95
31.74
2.81
31.74
0.95
19.34
0.52
19.34
0.64
18.20
0.53
18.20
0.62
11.05
0.50
11.05
0.42
10.73
0.56
10.73
0.42
4.89
0.45
4.89
0.28
4.89
4.77
4.89
0.27
1.00
69.03
2.88
69.03
0.98
57.12
1.49
57.12
0.90
3.77
57.74
2.27
57.74
3.74
49.64
1.97
49.64
3.74
3.18
49.58
2.05
49.58
3.14
43.18
1.90
43.18
3.12
2.28
37.12
2.57
37.12
2.31
32.82
1.86
32.82
2.29
2.34
26.42
1.97
25.06
1.59
25.06
1.92
22.63
1.62
22.63
1.89
1.83
19.57
1.37
20.50
1.91
20.50
1.35
19.75
1.78
19.75
1.34
2.86
92.48
1.00
73.83
4.56
73.83
1.00
59.77
2.46
59.77
1.00
4.77
83.62
7.00
68.21
3.93
68.21
7.11
57.84
3.72
57.84
7.03
4.48
70.11
6.03
65.02
3.35
65.02
6.14
55.37
2.80
55.37
6.21
3
52.14
4.70
52.14
5.23
55.73
4.30
55.73
5.54
49.06
2.96
49.06
5.24
4
42.82
4.39
42.82
4.50
42.21
4.40
42.21
4.49
39.61
3.76
39.61
4.49
5
35.24
4.18
35.24
3.77
31.77
4.94
31.77
3.82
31.39
4.47
31.39
3.79
47.09
2.62
47.09
2.54
43.02
2.33
43.02
2.56
37.67
2.26
37.67
2.53
CE AC
6
CPLEX
OF (%)
Avg. 1000
α = 0.8 SA
CPU(s)
PT
14
CPLEX
OF (%)
Avg. 700
α = 0.6 SA
CR IP T
6
CPLEX
AN US
400
r
M
p
ED
β
0
98.04
0.55
98.04
1.00
89.93
0.51
89.93
1.00
81.51
0.46
81.51
1.00
1
82.19
4.26
82.19
1.41
72.99
2.45
72.99
1.36
71.29
1.92
71.29
1.35
2
59.88
2.00
59.88
0.99
39.46
1.60
39.46
1.00
43.45
1.35
43.45
0.96
3
30.20
1.65
30.20
0.66
23.94
1.65
23.94
0.65
31.73
1.49
31.73
0.64
4
14.98
1.55
14.98
0.44
16.84
1.62
16.84
0.44
21.17
1.48
21.17
0.44
5
8.70
0.53
8.70
0.28
8.70
0.81
8.70
0.30
8.70
1.03
8.70
0.29
0
99.86
1.62
99.86
1.00
94.53
1.46
94.53
1.00
84.51
1.41
84.51
1.00
1
98.48
9.91
98.48
3.84
88.07
5.61
88.07
3.85
81.13
6.76
81.13
3.81
2
90.87
5.99
90.87
3.21
76.21
6.56
76.21
3.21
78.83
2.48
78.83
3.20
3
78.19
8.14
78.19
2.31
58.84
3.77
58.84
2.31
73.37
5.24
73.37
2.34
4
67.12
12.01
67.12
2.00
44.49
5.44
44.49
1.98
53.86
6.40
53.86
2.06
5
44.04
12.18
44.04
1.39
26.08
2.93
26.08
1.39
40.72
5.33
40.72
1.36
0
99.96
4.27
99.96
1.04
96.41
2.69
96.41
1.05
85.81
5.14
85.81
1.09
1
99.08
26.09
99.08
7.16
93.96
21.70
93.96
7.14
83.96
14.01
83.96
6.95
2
96.80
8.86
96.80
6.25
89.19
5.17
89.19
6.21
81.92
4.41
81.92
6.17
3
93.79
14.06
93.79
5.38
83.06
6.27
83.06
5.34
79.65
5.40
79.65
5.31
4
88.69
26.74
88.69
4.54
75.92
8.99
75.92
4.53
74.79
7.57
74.79
4.53
5
79.10
38.90
79.10
3.81
64.50
38.38
64.50
4.18
66.63
23.13
66.63
3.82
73.88
9.96
73.88
2.59
63.50
6.53
63.50
2.60
63.50
5.27
63.50
2.57
22
ACCEPTED MANUSCRIPT
Tables 10 and 11 summarize the results of solving the BHIMCP using the proposed SA algorithm for the CAB and TR data sets, respectively. It can be seen from these tables that considering the upcoming disruptions in leader’s decisions helps him/her prevent huge losses in the total covered flows after the follower’s interdictions. For instance, in the CAB data set with p = 5, α = 0.4, and β = 1500, if the leader chooses its hubs according to optimal solution of UMApHMCP, its total covered flows will plummet from 99.08% to 2.00% when the follower interdicts 4 hubs, whereas in case the leader decides on the location of its hubs based on the BHIMCP, its total covered flow will stay above 52% after the follower’s attack of the same intensity. From a CPU time perspective, it can be observed that
OF (%)
CPU(s)
OF (%)
CPU(s)
OF (%)
CPU(s)
1
19.97
54.91
19.66
52.98
CR IP T
the proposed SA solves the BHIMCP with the CAB instances in less than 75 seconds, on average. The solution
19.66
53.00
2
13.38
25.14
7.72
24.93
11.16
25.45
1
16.59
96.39
22.19
97.74
19.66
101.18
2
8.71
57.24
8.85
56.09
11.03
56.90
times for the TR instances, though much longer than those of the CAB instances, are quite acceptable for such a large scale bilevel MIP problem.
Table 10: SA results for the BHIMCP with the CAB data set
500
3
4
5
r
7.69
27.05
7.69
25.53
7.72
25.64
33.39
161.72
26.45
169.25
19.14
158.41
2
14.32
95.30
20.66
93.62
13.76
96.69
3
8.38
62.41
4
2.90
11.79
55.52
8.85
52.73
25.19
4.04
25.66
1.57
26.85
13.93
67.26
14.34
66.81
12.51
66.32
1
44.21
57.08
39.97
56.10
44.84
53.35
2
39.97
27.11
36.83
27.03
8.75
27.30
1
39.20
110.34
54.09
104.11
52.57
136.71
2
17.29
56.37
39.97
56.89
15.40
56.38
3
9.47
31.08
3.66
28.27
2.71
28.71
1
64.48
172.69
56.32
168.69
54.06
163.03
2
16.23
98.54
19.55
97.54
23.07
98.36
3
12.05
55.32
16.56
58.90
13.13
61.31
PT
4
AC
CE
5
4
3.68
30.10
3.56
31.09
2.71
29.35
27.40
70.96
30.06
69.85
24.14
72.72
1
71.12
69.19
68.50
89.08
68.26
69.82
2
58.72
29.80
57.50
29.31
62.87
29.58
1
71.55
98.13
72.84
99.59
72.84
101.96
2
55.71
61.21
62.77
56.12
64.89
55.40
3
13.67
29.51
52.30
29.31
58.72
32.11
1
75.80
164.75
74.47
168.55
72.12
159.61
2
59.17
103.21
66.23
106.09
63.36
107.09
3
58.80
65.11
63.79
63.19
59.45
63.10
4
52.30
30.15
58.66
35.91
28.21
29.33
57.43
72.34
64.12
75.24
61.19
72.00
Avg.
1500
3
4
5
Avg.
α = 0.8
1
ED
3
α = 0.6
3
Avg. 1000
α = 0.4
AN US
p
M
β
23
ACCEPTED MANUSCRIPT
Table 11: SA results for the BHIMCP with the TR data set
6
10
14
6
CPU(s)
OF (%)
CPU(s)
1
20.98
702.29
15.96
669.57
12.93
678.65
2
14.64
427.46
11.97
418.87
10.08
417.69
3
8.62
320.42
8.15
326.11
7.47
310.66
4
5.61
175.24
5.09
194.12
4.07
224.90
5
1.87
97.00
2.20
78.29
2.20
77.90
1
31.53
2464.91
23.69
2402.08
17.52
2458.94
2
25.10
2060.72
17.98
2072.12
14.92
1966.24
3
17.41
1823.73
14.26
1755.72
11.75
1722.50
4
13.12
1495.01
11.42
1485.62
9.64
1453.09
5
11.40
1010.98
8.69
1008.55
7.41
1008.69
1
39.39
5428.22
28.00
5960.60
20.08
5899.21
2
32.28
5221.50
22.10
4146.33
17.96
5715.84
3
25.87
4948.01
18.76
4652.70
15.46
4581.11
4
22.92
4444.45
15.66
4109.87
12.56
3645.88
5
18.57
2977.73
12.54
3018.91
10.02
3074.12
19.29
2239.84
14.43
2153.30
11.60
2215.69
1
57.20
646.62
49.40
733.08
44.13
624.23
2
43.71
443.03
40.19
429.37
38.39
426.94
3
34.62
366.46
28.25
341.09
29.16
331.07
4
25.20
185.72
17.27
178.69
25.09
187.98
5
16.93
90.87
16.93
121.68
11.50
84.80
1
77.75
2487.87
62.80
2401.69
54.05
2531.53
2
68.79
2026.81
56.06
2046.96
48.84
2013.65
3
56.45
1749.89
47.54
1852.62
44.03
1791.32
4
41.69
1546.60
38.98
1494.46
37.29
1541.08
5
25.98
970.76
34.85
1063.05
33.32
981.34
1
88.15
5028.92
69.96
5856.44
57.97
5393.62
2
83.01
5004.26
66.16
5128.45
55.69
4537.91
3
68.82
4590.43
56.44
4402.63
51.78
4327.11
4
63.22
4052.19
53.58
3789.18
47.54
3811.55
48.63
3190.87
49.46
3235.19
44.18
3285.95
53.34
2158.75
45.86
2204.97
41.53
2124.67
1
90.79
642.43
81.20
643.09
76.90
639.54
2
66.77
454.90
74.48
440.53
68.75
434.67
3
55.36
350.25
49.45
349.93
57.99
325.81
4
41.09
210.33
48.66
214.66
46.97
223.90
5
37.64
106.99
37.64
96.61
37.18
91.45
1
98.55
3813.87
91.79
3751.44
82.53
3890.29
2
95.25
3155.20
87.03
3113.09
79.72
3207.09
3
86.64
3005.69
82.39
2814.50
75.18
2778.65
4
75.96
2307.07
75.93
2289.30
57.24
2391.02
5
52.19
1485.20
45.89
1543.53
56.07
1410.71
1
99.51
5330.66
94.67
5732.49
84.43
5653.48
2
98.75
4367.81
90.47
4941.17
82.80
5038.05
3
95.16
4149.00
88.98
4379.11
80.58
4241.21
4
89.63
3565.24
86.39
3844.61
75.28
3901.87
5
83.52
3054.74
79.26
3192.73
72.71
3013.17
77.79
2399.96
74.28
2489.79
68.96
2482.73
PT
14
5
Avg. 6
AC
CE
1000
10
14
Avg.
α = 0.8
OF (%)
ED
10
α = 0.6
CPU(s)
Avg. 700
α = 0.4 OF (%)
CR IP T
400
r
AN US
p
M
β
24
ACCEPTED MANUSCRIPT
5.5. Results for hub interdiction center problems Results obtained by solving the r-HICP by using CPLEX as well as the proposed SA algorithm on the CAB and TR data sets are presented in Tables 12 and 13, respectively. To solve the follower’s problem, it is assumed that the leader has already located its hubs based on the solution of the uncapacitated multiple allocation p-hub center problem (UMApHCP). In this case, the follower tries to maximize the maximum distance between all O/D pairs by attacking some of the leader’s hubs.
α = 0.4 p
r
α = 0.6
CPLEX
SA
α = 0.8
CPLEX
SA
CPLEX
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
0
2064.66
0.39
2064.66
0.05
2243.76
0.38
2243.76
0.05
1
5007.65
0.44
5007.65
0.06
4915.25
0.23
4915.25
0.07
2
5451.58
0.60
5451.58
0.06
5056.95
0.16
5056.95
0.06
0
1774.44
0.47
1774.44
0.05
2127.12
0.33
2127.12
0.05
1
2736.53
0.47
2736.53
0.07
2752.40
0.26
2752.40
0.07
2
5007.65
1.28
5007.65
0.06
5007.65
0.90
5007.65
0.07
3
5451.58
0.53
5451.58
0.06
5451.58
0.77
5451.58
0.06
0
1599.73
0.58
1599.73
0.05
1916.15
0.35
1916.15
1
2390.81
0.58
2390.81
0.10
2583.17
0.33
2
3012.90
1.55
3012.90
0.07
4598.85
3
5007.65
1.01
5007.65
0.06
4
5451.58
0.70
5451.58
Avg.
3746.39
0.72
3746.39
4
CPU(s)
OF
CPU(s)
2515.57
0.23
2515.57
0.05
4933.31
0.21
4933.31
0.06
5148.16
0.43
5148.16
0.06
2437.70
0.24
2437.70
0.05
4915.25
0.25
4915.25
0.08
5053.12
0.83
5053.12
0.07
5406.80
0.59
5406.80
0.06
0.05
2288.78
0.38
2288.78
0.05
2583.17
0.09
2824.84
0.29
2824.84
0.09
1.07
4598.85
0.07
4598.85
0.78
4598.85
0.07
5229.61
1.64
5229.61
0.08
5229.61
1.24
5229.61
0.08
0.07
5451.58
0.90
5451.58
0.06
5451.58
0.85
5451.58
0.05
0.06
3944.50
0.61
3944.50
0.06
4233.63
0.53
4233.63
0.06
M
5
SA
OF
AN US
OF 3
CR IP T
Table 12: Results for the r-HICP with the CAB data set
Table 13: Results for the r-HICP with the TR data set
CPLEX CPU(s)
0
1216.00
2.98
1
1522.20
20.03
2
1834.20
205.84
3
3014.00
4
14
CPLEX
α = 0.8 SA
CPLEX
SA
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
OF
CPU(s)
1216.00
0.10
1467.00
2.01
1467.00
0.10
1755.40
2.14
1755.40
0.10
1522.20
0.68
1728.20
14.19
1728.20
0.64
1876.60
12.39
1876.60
0.68
1834.20
0.50
2282.00
65.10
2282.00
0.47
2282.00
32.07
2282.00
0.45
126.58
3014.00
0.32
3336.00
67.39
3336.00
0.31
3032.00
46.80
3032.00
0.30
3600.00
109.33
3600.00
0.23
3474.00
90.72
3474.00
0.20
3474.00
96.39
3474.00
0.19
5
3748.00
18.81
3748.00
0.12
3748.00
14.92
3748.00
0.13
3748.00
14.37
3748.00
0.12
0
1114.20
7.86
1114.20
0.50
1434.80
7.72
1434.80
0.50
1755.40
8.85
1755.40
0.50
1
1317.00
71.19
1317.00
1.89
1610.40
69.47
1610.40
1.89
1870.60
45.06
1870.60
1.87
2
1514.80
1063.41
1514.80
1.59
1719.00
854.86
1719.00
1.60
1915.20
260.55
1915.20
1.57
3
2282.00
1286.67
2282.00
1.11
1856.00
909.92
1856.00
1.13
1934.20
1184.03
1934.20
1.13
4
2330.00
1271.66
2330.00
0.97
2096.00
825.82
2096.00
0.94
2096.00
467.43
2096.00
0.94
5
2352.00
2328.94
2352.00
0.68
2282.00
1866.93
2282.00
0.67
2352.00
1052.48
2352.00
0.66
0
1114.20
25.89
1114.20
1.00
1434.80
25.19
1434.80
1.10
1755.40
23.63
1755.40
1.00
1
1277.20
232.11
1277.20
4.00
1518.80
201.80
1518.80
3.51
1876.60
174.75
1876.60
3.52
2
1382.80
1111.73
1382.80
3.09
1608.80
1039.72
1608.80
3.10
1915.80
638.95
1915.80
3.03
3
1690.00
2442.43
1690.00
2.63
1750.40
1536.54
1750.40
2.63
1966.00
2831.04
1966.00
2.63
4
1964.00
3423.79
1964.00
2.24
1827.20
2649.48
1827.20
2.23
2096.00
2508.86
2096.00
2.23
5
2282.00
8728.92
2282.00
1.92
1964.00
8058.31
1964.00
1.90
2282.00
3419.66
2282.00
1.86
1975.25
1248.78
1975.25
1.30
2063.18
1016.67
2063.18
1.28
2221.28
712.19
2221.28
1.26
AC
10
α = 0.6 SA
PT
OF
CE
6
r
ED
α = 0.4 p
Avg.
Note that the proposed SA algorithm is able to obtain the optimal solutions of the r-HICP for all instances in the 25
ACCEPTED MANUSCRIPT
CAB and TR data sets in very short computational times. By comparing the CPU times between the SA heuristic and CPLEX, the efficiency of the SA algorithm becomes more evident. It can be seen that by increasing the discount factor (α), the average service level provided by the leader is diminished as the inter-hub transportation cost grow with the discount factor value. Observe also that, similar to previous cases, larger values for r result in considerably worsened objective function values for the leader. For example, for p = 14, r = 0, and α = 0.4 in Table 13, the service level provided by the leader as the maximum cost between the O/D pairs is 1114.20 based on the solution of the UMApHCP. However, after the follower attacks the network by disrupting 4 hubs, the maximum distance amounts
CR IP T
to 1964.00, showing an increment of 849.80 units, which can be considered as a serious damage to the service level provided by the leader.
Tables 14 and 15 show the results for solving the BHICP problem for the CAB and TR data sets, respectively. The results reported in these tables indicate that the proposed SA algorithm is able to solve the problem instances in very reasonable times.
3
4
5
r
α = 0.4 OF
CPU(s)
1
3012.90
85.10
2
3373.35
50.13
1
2454.39
143.80
2
3179.67
81.02
3
3529.58
41.18
1
2200.26
2
3012.90
3
4003.38
α = 0.6
α = 0.8
OF
CPU(s)
OF
CPU(s)
3012.90
77.32
3179.67
84.74
3471.87
50.06
3373.35
49.72
2645.59
141.91
2720.59
150.77
4072.25
79.10
3082.54
83.19
5451.58
49.73
3707.23
39.98
M
p
AN US
Table 14: SA results for the BHICP with the CAB data set
155.62
2427.42
165.09
2606.41
162.91
134.76
3529.58
131.45
3471.87
134.18
92.03
89.10
3745.39
92.83
4180.17
59.12
4598.85
44.58
4072.25
46.90
Avg.
3216.29
93.64
3726.64
92.04
3328.81
93.91
ED
4329.71
4
PT
Similar to the previous cases, it can be observed that the leader can benefit by solving the BHICP to reduce the adverse effects of the follower’s disruptive action. For instance, as stated above, for p = 14 and α = 0.4 in the TR
CE
data set, if the leader locates its hubs based on the UMApHCP, the follower will diminish the maximum O/D cost by 849.80 units using 4 interdictions, whereas acting based on the BHICP, the leader’s service level will be worsened at most by 409.40 units. Therefore, it is highly beneficial for the leader to solve the bilevel model and locate its hubs
AC
accordingly.
5.6. Optimal hub locations on map In this section, we graphically illustrate how the optimal hub locations for the leader and optimal interdictions
by the follower change as the leader makes its decisions with and without being aware of the follower’s upcoming attack. The analysis is conducted for all the three variants of the hub interdiction problem. Figure 1 depicts the optimal solutions for the leader and the follower for the hub interdiction median problem on the CAB data set with p = 5, r = 2, and α = 0.8. The leader’s installed hubs are shown as squares, whereas the follower’s interdicted hubs are shown as cross marks. In part (a), it is assumed that the leader is unaware of the follower’s subsequent attack 26
ACCEPTED MANUSCRIPT
Table 15: SA results for the BHICP with the TR data set r
10
14
α = 0.6
α = 0.8
CPU(s)
OF
CPU(s)
OF
CPU(s)
1
1468.20
792.38
1663.80
788.61
1838.40
809.84
2
1590.00
543.34
1772.00
533.11
1912.00
517.02
3
2338.00
320.73
2352.00
321.09
2282.00
311.31
4
2684.00
200.99
2352.00
207.41
2732.00
210.34
5
2684.00
93.11
2512.00
92.16
3190.00
93.79
1
1254.40
2410.84
1518.80
2406.59
1791.40
2388.43
2
1420.40
2483.43
1608.80
2238.96
1870.60
2063.05
3
1488.40
1876.38
1715.00
1733.65
1882.60
2109.98
4
2330.00
1449.11
1906.60
1678.11
1964.20
1820.11
5
2352.00
1022.32
2282.00
1201.85
2338.00
972.10
1
1254.40
4624.10
1518.80
4677.14
1791.40
4600.28
2
1376.40
3984.46
1605.60
4011.35
1870.60
4126.63
3
1468.20
4261.13
1676.60
4283.20
1882.00
4310.68
4
1523.60
4997.60
1737.80
3590.14
1964.00
4097.01
5
1836.00
3107.19
1804.53
2144.47
Avg.
CR IP T
6
α = 0.4 OF
AN US
p
1964.00
3244.09
2282.00
3951.50
1879.05
2067.16
2106.08
2158.80
and locates its hubs based on the solution of the UMApHMP, whereas part (b) depicts the optimal solutions in case
CE
PT
ED
M
the leader locates its hubs based on the BHIMP. As can be seen from the Figure 1, the optimal hub locations as well
(a) leader acts based on the UMApHMP
(b) leader acts based on the BHIMP
AC
Figure 1: Optimal hub locations and interdictions for hub median problem with p = 5, r = 2 and α = 0.8
as interdictions are different for the cases in which the leader acts based on the r-HIMP and BHIMP. This shows the importance of solving the BHIMP, and locating the hubs accordingly, for the leader. Figure 2 depicts the optimal solutions for the leader and the follower for the hub interdiction maximal covering problem with β = 1500. Here, it can also be seen that the optimal network configurations are quite different as the leader acts whether based on the r-HIMCP or BHIMCP. Finally, the optimal solutions for the leader and the follower for the hub interdiction center problem are depicted in Figure 3. Similar to the two previous cases, it can also be observed that the optimal set of hubs for the leader as well as the optimal interdictions for the follower are different 27
CR IP T
ACCEPTED MANUSCRIPT
(a) leader acts based on the UMApHMCP
(b) leader acts based on the BHIMCP
Figure 2: Optimal hub locations and interdictions for hub maximal covering problem with p = 5, r = 2, α = 0.8 and β = 1500
ED
M
AN US
based on whether the leader acts in anticipation of the follower’s attack or not.
(a) leader acts based on the UMApHCP
(b) leader acts based on the BHICP
PT
Figure 3: Optimal hub locations and interdictions for hub center problem with p = 5, r = 2 and α = 0.8
CE
5.7. Evaluating the leader’s objective value improvements Results presented in the previous sections reveal that, when the leader acts in anticipation of follower’s attack, i.e., he/she acts based on the solution of the BHIMP, BHIMCP, or BHICP, the effect of follower’s disruptive action
AC
is less severe than that of the case where the leader is unaware of the upcoming attack and acts based on the solution of the UMApHMP, UMApHMCP, or UMApHCP, respectively. In other words, acting based on the bilevel models, creates value for the leader. As a last part of our numerical analyses, in this section, we evaluate the changes (in percent) in the leader’s objective function value when he/she acts according to the bilevel models based on the presented results in preceding sections. To this end, we define and use a metric called the Objective Function Change Percentage (OFCP), as follows: OFCP =
OFB − OFS × 100% OFS
28
ACCEPTED MANUSCRIPT
in which OFB is the average value of leader’s objective function when he/she acts based on the bilevel models (BHIMP, NBHIMCP, or BHICP), whereas OFS is the average value of leader’s objective function when the locations are selected according to the single level models (UMApHMP, UMApHMCP, or UMApHCP). The averages are taken with respect to the three discount factor (α) values. The OFCP values for different instances of the CAB data set with the three variants of the problem are presented in Table 16. The first two columns show the number of hubs opened by the leader and the number of hubs interdicted by the follower, respectively. The third column give the OFCP values for the hub median variant of the problem. Similarly, the next three columns show the OFCP values
the OFCP values for the hub center variant of the problem.
CR IP T
for the hub maximal covering variant under the three covering radius (β) values. Finally, the last column present
Table 16: Average OFCP values for the CAB data set
p
Median
Maximal covering β=500
β=1000
β=1500
Center
-0.20%
30.65%
117.13%
33.74%
-38.04%
2
-57.02%
60.50%
770.30%
688.60%
-34.73%
1
-4.47%
17.63%
39.11%
11.95%
-19.24%
2
-9.93%
10.05%
193.22%
122.97%
-31.42%
3
-55.87%
427.40%
61.14%
449.05%
-22.20%
1
-11.26%
39.37%
21.55%
11.84%
-7.24%
2
-18.41%
49.56%
32.72%
98.05%
-17.99%
3
-21.45%
195.52%
74.79%
368.93%
-21.91%
4
-51.35%
243.15%
1.22%
776.94%
-21.42%
Avg.
-25.55%
119.31%
145.69%
284.67%
-23.80%
4
ED
5
AN US
1
M
3
r
Observe that for the CAB data set, the average improvement in the leader’s objective value is 25.55% for hub median problem and 23.80% for the hub center variant. For the hub maximal covering variant, on the other hand, the
PT
improvements are substantially higher than the other two variants. The leader’s objective has improved respectively 119.31%, 145.69%, and 284.67%, on average, for β = 500, 1000, and 1500 by acting based on the bilevel models.
CE
An interesting observation is that, the average OFCP values get larger as the follower disrupts more hubs. In other words, the benefit of using bilevel models are more considerable for the leader when facing a powerful interdictor. Table 17 shows the OFCP values for different instances of the TR data set with the three variants of the problem.
AC
As can be seen from the Table 17, the average improvements in the value of leader’s objective with the TR data
set are 7.65% and 10.51% for hub median and center variants, respectively. However, for the hub maximal covering variant, the leader’s objective has improved 16.57%, 49.18%, and 49.85%, on average, for β = 400, 700, and 1000, respectively. Considering the improvement percentages for the CAB and TR data sets, one can observe that the improvements are larger for the CAB data set (especially for the maximal covering variant) compared to those of the TR data set. This difference can partially be explained by the existence of some very high traffic nodes in the CAB data set which if interdicted, inflict a high damage to the network’s efficiency. Therefore, proper location of hubs by the leader as a result of acting based on the bilevel model, can prevent such harmful damages. Another
29
ACCEPTED MANUSCRIPT
Table 17: Average OFCP values for the TR data set
r
14
Center
β=400
β=700
β=1000
-1.29%
12.88%
6.35%
9.90%
-3.05%
2
-8.02%
22.22%
24.54%
47.07%
-17.57%
3
-11.00%
33.48%
62.20%
89.59%
-25.69%
4
-21.99%
34.39%
114.34%
158.01%
-26.36%
5
-36.03%
20.35%
236.50%
330.88%
-25.42%
1
-0.92%
5.73%
11.77%
1.94%
-4.86%
2
-1.48%
18.05%
23.24%
6.54%
3
-2.50%
21.39%
42.44%
16.07%
-16.24%
4
-4.13%
24.70%
59.17%
26.39%
-4.92%
5
-7.65%
34.34%
57.39%
39.07%
-0.20%
1
-2.17%
0.19%
3.06%
0.58%
-2.31%
2
-3.75%
1.39%
7.54%
1.53%
3
-3.07%
2.33%
12.81%
3.20%
-7.02%
4
-5.00%
9.30%
31.85%
4.97%
-11.24%
5
-5.76%
7.78%
44.58%
12.02%
-6.83%
-7.65%
16.57%
49.18%
49.85%
-10.51%
1
10
Maximal covering
Avg.
CR IP T
6
Median
-4.84%
-1.12%
AN US
p
M
important factor which may contribute to this difference is the limited number of candidate nodes for locating the hubs in the TR data set (|H| = 22). This limited set of candidate hub nodes make the action space of the decision makers limited, which in turn results in less difference between the set of leader’s hubs when he/she decides based
ED
on either of single level or bilevel models. Nevertheless, based on the results presented in the above tables, it can be concluded that the leader will substantially benefit via acting based on the bilevel models by reducing the adverse
6. Conclusions
PT
effects of the follower’s attack.
CE
In this paper, we addressed three different variants of the hub interdiction problem which has attracted the attention of researchers in recent years. To this end, we considered the multiple allocation p-hub median, the p-
AC
hub maximal covering, and the p-hub center problems under intentional disruptions by assuming a leader-follower Stackelberg game where the decision maker who makes the initial location decisions is called the leader, whereas the one who decides to interdict the installed hubs is called the follower. New formulations were proposed in each case which model the problems from both the leader’s and the follower’s perspectives. Furthermore, different types of CA constraints were examined in order to identify the one with the best performance in terms of solution time by a commercial solver. To solve the problems, we proposed efficient simulated annealing (SA) heuristics. Extensive computational experiments based on the CAB and TR data sets were conducted to analyze different properties of the problems and to evaluate the performance of the proposed SA heuristics as well as the mathematical formulations. The results demonstrated the high efficiency of the presented models and the SA algorithms both in terms of solution
30
ACCEPTED MANUSCRIPT
quality and computational time. Our research can be extended in several ways. An interesting research direction is to consider capacity restrictions on hubs or a limited budget for installing and/or disrupting the hubs. In addition, the assumption that an attack on a given hub is always successful causing complete disruption can be relaxed resulting in a problem in which attacks are successful only with a given probability or only cause partial disruptions. Moreover, developing exact methods such as branch-and-bound or decomposition algorithms might be of interest for the researchers.
CR IP T
Acknowledgments The authors sincerely thank the three anonymous referees for their valuable comments and suggestions that have helped improve the paper considerably.
References
AN US
˙ G¨ [1] Akg¨ un, I., um¨ u¸sbuga, F., Tansel, B., 2015. Risk based facility location by using fault tree analysis in disaster management. Omega 52, 168 – 179.
[2] Aksen, D., Piyade, N., Aras, N., 2010. The budget constrained r-interdiction median problem with capacity expansion. Central European Journal of Operations Research 18 (3), 269–291. [3] Aksen, D., engl Akca, S., Aras, N., 2014. A bilevel partial interdiction problem with capacitated facilities and
M
demand outsourcing. Computers & Operations Research 41, 346 – 358.
[4] Alumur, S., Kara, B. Y., 2008. Network hub location problems: The state of the art. European Journal of
ED
Operational Research 190 (1), 1–21.
[5] Bard, J. F., 1998. Practical bilevel optimization: algorithms and applications. Kluwer Academic Publishers,
PT
Dordrecht, The Netherlands.
[6] Berman, O., Drezner, Z., Wesolowsky, G. O., 2003. Locating service facilities whose reliability is distance
CE
dependent. Computers & Operations Research 30 (11), 1683–1695. [7] Berman, O., Krass, D., Menezes, M. B., 2007. Facility reliability issues in network p-median problems: strategic
AC
centralization and co-location effects. Operations Research 55 (2), 332–350. [8] Berman, O., Krass, D., Menezes, M. B., 2009. Locating facilities in the presence of disruptions and incomplete information. Decision Sciences 40 (4), 845–868.
[9] Berman, O., Krass, D., Menezes, M. B., 2013. Location and reliability problems on a line: Impact of objectives and correlated failures on optimal location patterns. Omega 41 (4), 766 – 779. [10] Campbell, J. F., 1992. Location and allocation for distribution systems with transshipments and transportion economies of scale. Annals of operations research 40 (1), 77–99.
31
ACCEPTED MANUSCRIPT
[11] Campbell, J. F., 1994. Integer programming formulations of discrete hub location problems. European Journal of Operational Research 72 (2), 387–405. [12] Campbell, J. F., 1996. Hub location and the p-hub median problem. Operations Research 44 (6), 923–935. [13] Campbell, J. F., O’Kelly, M. E., 2012. Twenty-five years of hub location research. Transportation Science 46 (2), 153–169.
rithm. Journal of optimization theory and applications 45 (1), 41–51.
CR IP T
ˇ [14] Cern` y, V., 1985. Thermodynamical approach to the traveling salesman problem: An efficient simulation algo-
[15] Chaharsooghi, S., Momayezi, F., Ghaffarinasab, N., 2017. An adaptive large neighborhood search heuristic for solving the reliable multiple allocation hub location problem under hub disruptions. International Journal of Industrial Engineering Computations 8 (2), 191–202.
[16] Church, R., Scaparra, M. P., 2007. Analysis of Facility Systems’ Reliability When Subject to Attack or a Natural
AN US
Disaster. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 221–241.
[17] Church, R. L., Cohon, J. L., 1976. Multiobjective location analysis of regional energy facility siting problems. Tech. rep., Brookhaven National Lab., Upton, NY (USA).
[18] Church, R. L., Scaparra, M. P., Middleton, R. S., 2004. Identifying critical infrastructure: the median and
M
covering facility interdiction problems. Annals of the Association of American Geographers 94 (3), 491–502. [19] Contreras, I., 2015. Hub location problems. In: Location Science. Springer International Publishing, pp. 311–344.
ED
[20] Contreras, I., Fernndez, E., Reinelt, G., 2012. Minimizing the maximum travel time in a combined model of facility location and network design. Omega 40 (6), 847 – 860. [21] Cui, T., Ouyang, Y., Shen, Z.-J. M., 2010. Reliable facility location design under the risk of disruptions.
PT
Operations Research 58 (4-part-1), 998–1011.
CE
[22] Dempe, S., 2002. Foundations of bilevel programming. Springer Science & Business Media. [23] Dobson, G., Karmarkar, U. S., 1987. Competitive location on a network. Operations Research 35 (4), 565–574.
AC
[24] Drezner, Z., 1987. Heuristic solution methods for two location problems with unreliable facilities. Journal of the Operational Research Society 38 (6), 509–514.
[25] Ebery, J., Krishnamoorthy, M., Ernst, A., Boland, N., 2000. The capacitated multiple allocation hub location problem: Formulations and algorithms. European Journal of Operational Research 120 (3), 614–631.
[26] Ernst, A. T., Hamacher, H., Jiang, H., Krishnamoorthy, M., Woeginger, G., 2009. Uncapacitated single and multiple allocation p-hub center problems. Computers & Operations Research 36 (7), 2230–2241. [27] Ernst, A. T., Krishnamoorthy, M., 1996. Efficient algorithms for the uncapacitated single allocation p-hub median problem. Location science 4 (3), 139–154. 32
ACCEPTED MANUSCRIPT
[28] Ernst, A. T., Krishnamoorthy, M., 1998. Exact and heuristic algorithms for the uncapacitated multiple allocation p-hub median problem. European Journal of Operational Research 104 (1), 100–112. [29] Espejo, I., Mar´ın, A., Rodr´ıguez-Ch´ıa, A. M., 2012. Closest assignment constraints in discrete location problems. European Journal of Operational Research 219 (1), 49–58. [30] Farahani, R. Z., Hekmatfar, M., Arabani, A. B., Nikbakhsh, E., 2013. Hub location problems: A review of models, classification, solution techniques, and applications. Computers & Industrial Engineering 64 (4), 1096–
CR IP T
1109. [31] Gerrard, R. A., Church, R. L., 1996. Closest assignment constraints and location models: properties and structure. Location Science 4 (4), 251–270.
[32] Ghaffarinasab, N., Motallebzadeh, A., Jabarzadeh, Y., Kara, B. Y., 2018. Efficient simulated annealing based solution approaches to the competitive single and multiple allocation hub location problems. Computers &
AN US
Operations Research 90 (Supplement C), 173 – 192.
[33] Hamacher, H. W., Meyer, T., 2006. Hub cover and hub center problems. Technische Universit¨ at Kaiserslautern, Fachbereich Mathematik.
[34] Hong, J.-D., Xie, Y., Jeong, K.-Y., 2012. Development and evaluation of an integrated emergency response facility location model. Journal of Industrial Engineering and Management 5 (1), 4.
M
[35] Jabalameli, M. S., Barzinpour, F., Saboury, A., Ghaffari-Nasab, N., 2012. A simulated annealing-based heuristic for the single allocation maximal covering hub location problem. International Journal of Metaheuristics 2 (1),
ED
15–37.
[36] Kara, B., 1999. Modeling and analysis of issues in hub location problems. Doctor of Philosophy Thesis, Bilkent
PT
University, Ankara, Turkey.
[37] Kara, B., Tansel, B., 2003. The single-assignment hub covering problem: Models and linearizations. Journal of
CE
the Operational Research Society 54 (1), 59–64. [38] Kara, B. Y., Tansel, B. C., 2000. On the single-assignment p-hub center problem. European Journal of Opera-
AC
tional Research 125 (3), 648–655. [39] Kirkpatrick, S., Vecchi, M. P., et al., 1983. Optimization by simmulated annealing. science 220 (4598), 671–680. [40] Lee, S.-D., 2001. On solving unreliable planar location problems. Computers & Operations Research 28 (4), 329–344.
[41] Lei, T. L., 2010. Location modeling utilizing closest and generalized closest transport/interaction assignment constructs. University of California, Santa Barbara. [42] Lei, T. L., 2013. Identifying critical facilities in hub-and-spoke networks: A hub interdiction median problem. Geographical Analysis 45 (2), 105–122. 33
ACCEPTED MANUSCRIPT
[43] Lei, T. L., Tong, D., 2013. Hedging against service disruptions: an expected median location problem with site-dependent failure probabilities. Journal of Geographical Systems 15 (4), 491–512. [44] Li, X., Ouyang, Y., 2010. A continuum approximation approach to reliable facility location design under correlated probabilistic disruptions. Transportation research part B: methodological 44 (4), 535–548. [45] Liberatore, F., Scaparra, M. P., Daskin, M. S., 2011. Analysis of facility protection strategies against an uncertain number of attacks: the stochastic r-interdiction median problem with fortification. Computers & Operations
CR IP T
Research 38 (1), 357–366. [46] Lim, M., Daskin, M. S., Bassamboo, A., Chopra, S., 2010. A facility reliability problem: formulation, properties, and algorithm. Naval Research Logistics (NRL) 57 (1), 58–70.
[47] Losada, C., Scaparra, M. P., Church, R. L., Daskin, M. S., 2012. The stochastic interdiction median problem with disruption intensity levels. Annals of Operations Research 201 (1), 345–365.
AN US
[48] Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., Teller, E., 1953. Equation of state calculations by fast computing machines. The journal of chemical physics 21 (6), 1087–1092. [49] O’Hanley, J. R., Church, R. L., 2011. Designing robust coverage networks to hedge against worst-case facility losses. European Journal of Operational Research 209 (1), 23–36.
[50] O’Kelly, M. E., 1987. A quadratic integer program for the location of interacting hub facilities. European Journal
M
of Operational Research 32 (3), 393–404.
[51] Parvaresh, F., Golpayegany, S. H., Husseini, S. M., Karimi, B., 2013. Solving the p-hub median problem under
ED
intentional disruptions using simulated annealing. Networks and Spatial Economics 13 (4), 445–470. [52] Parvaresh, F., Husseini, S. M., Golpayegany, S. H., Karimi, B., 2014. Hub network design problem in the
PT
presence of disruptions. Journal of Intelligent Manufacturing 25 (4), 755–774. [53] Peker, M., Kara, B. Y., 2015. The p-hub maximal covering problem and extensions for gradual decay functions.
CE
Omega 54, 158–172.
[54] Peng, P., Snyder, L. V., Lim, A., Liu, Z., 2011. Reliable logistics networks design with facility disruptions.
AC
Transportation Research Part B: Methodological 45 (8), 1190–1211. [55] Rojeski, P., ReVelle, C., 1970. Central facilities location under an investment constraint. Geographical Analysis 2 (4), 343–360.
[56] Saboury, A., Ghaffari-Nasab, N., Barzinpour, F., Jabalameli, M. S., 2013. Applying two efficient hybrid heuristics for hub location problem with fully interconnected backbone and access networks. Computers & Operations Research 40 (10), 2493–2507. [57] Scaparra, M. P., Church, R. L., 2008. A bilevel mixed-integer program for critical infrastructure protection planning. Computers & Operations Research 35 (6), 1905–1923. 34
ACCEPTED MANUSCRIPT
[58] Scaparra, M. P., Church, R. L., 2015. Location problems under disaster events. In: Laporte, G., Nickel, S., Saldanha da Gama, F. (Eds.), Location Science. Springer International Publishing, pp. 623–642. [59] Silva, M. R., Cunha, C. B., 2017. A tabu search heuristic for the uncapacitated single allocation p-hub maximal covering problem. European Journal of Operational Research 262 (3), 954 – 965. [60] Skorin-Kapov, D., Skorin-Kapov, J., O’Kelly, M., 1996. Tight linear programming relaxations of uncapacitated p-hub median problems. European Journal of Operational Research 94 (3), 582–593.
Transportation Science 39 (3), 400–416.
CR IP T
[61] Snyder, L. V., Daskin, M. S., 2005. Reliability models for facility location: the expected failure cost case.
[62] Snyder, L. V., Scaparra, M. P., Daskin, M. S., Church, R. L., 2006. Planning for disruptions in supply chain networks. Tutorials in operations research 2, 234–257.
AN US
[63] Tan, P. Z., Kara, B. Y., 2007. A hub covering model for cargo delivery systems. Networks 49 (1), 28–39. [64] Wagner, B., 2008. Model formulations for hub covering problems. Journal of the Operational Research Society 59 (7), 932–938.
[65] Wagner, J., Falkson, L., 1975. The optimal nodal location of public facilities with price-sensitive demand. Geographical Analysis 7 (1), 69–83.
M
[66] Yaman, H., Elloumi, S., 2012. Star p-hub center problem and star p-hub median problem with bounded path
AC
CE
PT
ED
lengths. Computers & Operations Research 39 (11), 2725–2732.
35