A Biobjective and Trilevel Programming Model for Hub Location Problem in Design of a Resilient Power Projection Network

Hubs disruptions are taken into account in design of a resilient power projection network. The problem is tackled from a multiple criteria decision-making (MCDM) perspective. Not only the network cost in normal state is considered, but also the cost in the worst-case situation is taken into account. A biobjective and trilevel integer programming model is proposed using game theory. Moreover, we develop a metaheuristic based on tabu search and shortest path algorithm for the resolution of the complex model. Computational example indicates that making tradeoffs between the performances of the network in different situations is helpful for designing a resilient network.


Introduction
Power projection is a term used in military to refer to the capacity of a state to apply national transportation network to rapidly and effectively deploy and sustain forces in and from multiple dispersed locations to respond to crises, to contribute to deterrence, and to enhance regional stability [1].This ability is a crucial element of a state's military power in modern warfare.It can not only guarantee the victory of the war but also help to achieve the military strategic objectives.Generally, in the process of power projection, the troops are firstly consolidated at transportation hubs by means of motorization march and then they are shipped between hubs links using different modes of transportation (highway, railway, air transport, etc.).Finally they are deployed to dispersed locations through motorization march according to the mission demands.The ability to project power relies on safe and efficient transportation network.Hub location-allocation strategies are critical issues in design of the projection network, which have significant influence on the efficiency and survivability of the network.
Traditional studies dealing with the hub location problem in transportation network usually locate hub facilities and allocate spokes to those hubs in order to minimize the total transportation cost in normal state [2][3][4].However, due to the military characteristics of power projection, the network faces intentional attacks from the enemy or interdictor.Some important hub facilities would be disrupted in wartime.In practice, a network would not be completely failed even though one or more nodes have become invalid.The network continues to serve the remaining connected components and would recover through emergency management measures, such as traffic rerouting.The ability of a network to withstand and reduce the impact of disruptions is called resilience [5], which is a key indicator of the network survivability.
In recent years, a number of works have been done emphasizing resilience of transportation systems.Ip and Wang [6] introduced network resilience as a function of the number of reliable paths between all node pairs.Zhao et al. [7] developed a resilience metrics including availability, connectivity, and accessibility in a military supply network.Resilience of a transportation network in Miller-Hooks et al. [8] is defined 2 Mathematical Problems in Engineering as the expected fraction of demand that can be satisfied after disaster with budget constraints.Konak and Bartolacci [9] incorporated network resilience into network designs using a stochastic genetic approach.Janić [10] dealt with estimating the resilience of an air transport network affected by a largescale disruptive event.More relevant works can be found in the review paper of Reggiani et al. [11].
Dealing with the hub location problem in design of a projection network, not only the transportation costs in normal state should be considered, but also the ability of withstanding and reducing the impact of disruptions in wartime should be taken into account.Hence, the hub location problem in power projection network is a multiple criteria decisionmaking (MCDM) problem from the resilience perspective.An MCDM problem refers to a situation in which there are at least two alternative courses of action to choose from and this choice is driven by the desire to meet multiple, often conflicting, objectives [12].There are different classifications of MCDM problems and methods.The classification proposed by Pardalos et al. [13] distinguishes four categories: (1) multiobjective mathematical programming, (2) multiattribute utility theory, (3) outranking relations approach, and (4) preference disaggregation approach.Generally, when the criteria are nonconflicting, an MCDM problem can be solved by combining multiple objectives into one objective using a specific method such as the weighted sum model (WSM) and the analytic hierarchy process (AHP) [14].On the contrary, when the criteria are conflicting, multiobjective optimization based on evolutionary algorithms is often used to identify the Pareto frontier [15].Methods such as multiobjective simulated annealing (MOSA) [16], Pareto iterated local search (PILS) [17], and others have been developed.For a review of MCDM, see de Almeida et al. [18].
Our contribution attempts to analyze the hub location problem in power projection network in view of game theory.We establish a biobjective and trilevel integer programming model, which can be used to make tradeoffs between the performance of the network in normal and interdicted state.A multiobjective tabu search approach based on shortest path algorithm is developed to solve the model.A computational example is presented to illustrate the model and algorithm.The paper ends with some concluding remarks.

Problem Description
In the process of power projection, multiple modes of transportation are used to deliver personnel and equipment, which belongs to intermodal transport.Intermodal transport is characterized by the combination of the advantages of rail and road, rail for long distances and large quantities, and road for collecting and distributing over short or medium distance [19].In this contribution, the motorized march process of power projection can be regarded as the road mode in the intermodal transport.And the backbone transport is usually served by railway or by highway.Commonly, the structure of a power projection network is hub-spoke distribution paradigm [20], in which all traffic moves along spokes connected to the hub at the center.In the power projection network, the starting and end points of the troops can be viewed as nonhub nodes, while the rally points can be regarded as the hub nodes.
It is clear that hubs have higher connectivity than nonhub nodes.In contrast to the failure of nonhubs and links, the failure of hubs has more impact on the resilience of the network [21].Therefore, in this paper, we focus on the issue of hub location strategy dealing with probable hub disruptions in design of a power projection network.When parts of the hub nodes are disrupted and cannot be restored in limited time, common resilience management actions are rerouting the traffic flows through the remaining connected hubs.Traffic flows rerouting can ensure the completion of projection mission.However, the cost and time of power projection will increase.Since the transportation time would greatly exceed the scheduled time if we reroute the traffic flow, we consider the resilience of the power projection network from a cost increasing perspective without time restraint.In this paper, the resilience of the power projection network is defined as follows [5]: where  is the resilience of the power projection network,  before is the transportation costs of the network in normal situation, and  after is the transportation costs of the network when some hub nodes have been disrupted.When dealing with the design of hub location strategy in power projection network, on the one hand, we need to consider minimize to total cost of the network in the normal state, that is Minimize( before ); on the other hand, we hope that the increased cost of the network would not be much higher once some hub nodes are disrupted, that is, Minimize( after ).However, in view of game theory, the enemy may select the hubs for attack which result in the greatest damage to the network [22], that is, Maximize( after ).Therefore, at the beginning of the design we should take the worst-case into account [23], that is, Minimize(Maximize( after )).Besides, after the interdiction, traffic rerouting should follow the leastcost principle; in other words, each origin-destination (O-D) pair is served by its least-cost route via the remaining hubs.Thus, there is  after = Minimize(  ), where   is the transportation cost of the network served by the remaining hubs after attack.
From what is mentioned above, we can conclude that the hub location strategy in power projection networks for resilience perspective is a biobjective and trilevel integer programming problem.

Mathematical Model
3.1.Mathematical Hypothesis.The power projection procedure can be viewed as the intermodal transportation process in hub-spoke network and the basic assumptions are as follows.
(1) The hub network is complete with the presence of a direct hub link between each hub pair.And not all the hubs would be disrupted at the same time.
(2) There is a cost discount factor using hub links due to economies of scale.
(3) A nonhub node can be allocated to more than one hub node and direct links between nonhub nodes are not permitted.
(4) The distance and transportation cost satisfy the triangle inequality, so that a shipment between any O-D pairs may travel at most two hub nodes.
(5) The hubs and links are uncapacitated.
(6) A hub facility loses all its capacity if it was disrupted.However, the traffic demand from or to this node would not be eliminated.That is, the hub node becomes a nonhub node after disruption.(7) The enemy would select the hubs for attack which result in the greatest damage to the network, and we should take the worst-case into account in design.
(8) When some hub nodes are disrupted, traffic rerouting should follow the least-cost rule.  : transportation cost discount factor between hubs using mode .

Decision Variables.
There are five sets of decision variables in the model: 1, if flow from  to  via hub pair (, ) using mode ; 0, otherwise, 1, if flow from  to  via hub pair (, ) using mode  after -disruption; 0, otherwise. (2)

Biobjective and Trilevel Integer Programming Model.
Based on the hypothesis and notations mentioned above, the hub location strategy in power projection networks form resilience perspective is modeled as The mathematical model above is a biobjective and trilevel integer programming.The first level program has two objective functions from the network designer's perspective, one of which seeks to minimize the total network cost in the normal situation presented by objective function (3), while the other seeks to minimize the total network cost in the worst-case (after -disruption) presented by objective function (4).In the objective function (3), the first term in right hand is the network transportation cost; the second term is the fixed costs of hubs.Constraint (5) requires that there are  hubs to be located.Constraints (6) refer to the choice of a specific arc for interhub travel, which ensure that each O-D pair (, ) is assigned to exactly one using one transportation mode between hubs.Constraints ( 7) and ( 8) guarantee that each O-D pair can only be assigned to an existing hub pair.Constraints (9) are binary requirements.
The second level program constitutes a worst-case scenario presented by objective function (10).In the worst case, the enemy (or attacker) would choose the hubs for attack which maximize the total network cost.The first term in the objective function (10) calculates the transportation cost after -disruption; the second term calculates the value of the  hubs which are interdicted.Constraint (11) implies that there are  hubs disrupted.Constraints (12) are integrity constraints.
The third level program is to determine the least-cost route for each O-D pair in the existing network after  hubs failure.The objective function is presented by ( 14) from the users' perspective.If there are  hubs disrupted, traffic flows in the network would be rerouted through the remaining hubs, as presented by constraints (15).Constraints ( 16) and (17) prevent an assignment to a hub that has been disrupted.Constraints (18) are binary requirements.

Multiobjective Tabu Search Algorithm
The model above is a combinatorial problem with multiple objectives and multilayer planning.It is very difficult to solve using quantitative methods.General practices are using the heuristic algorithms [24].In comparison with other heuristic algorithms, tabu search [25,26] is a well-known approach which can efficiently overcome local optimality entrapment through simulating the process of human thinking.It guides the search to explore the feasible region using a short-and a long-term memory of previously visited solutions.This procedure can easily find the best solution and escape from local optima.In this paper, we develop a heuristic approach based on Floyd shortest path algorithm and tabu search to solve the model.Due to the nested relations in the model, feasible solutions can be constructed through three steps.First, a set of  nodes are randomly generated as the hubs in normal situation.Second, all the possibly disrupted cases are enumerated for the current hub set.Then, for each case, the traffic assignment of all O-D pairs to the remaining hubs is determined by shortest path algorithm.

Floyd Shortest Path Algorithm.
During each iteration of the tabu search procedure, assignment problems of all O-D pairs to the given hubs are needed to determine the calculation of the objective functions (3), (10), and (14).The Floyd shortest path algorithm [27] can be employed to get the overall solutions efficiently.
Let  be the set of hubs.The weight matrix of the network is denoted as  = (  ) × .There is Denote  = ( ()  ) × as the route matrix of the network using vertices only from the set {1, 2, . . ., } as intermediate points along the way, where  (0)  =  for  = 0.Then, basic steps for Floyd shortest path algorithm are as follows.
Using the unit transportation cost matrix of the shortest path in the specific network, we can easily calculate the transportation cost terms in the objective functions (3), (10), and ( 14), respectively.

Constitution of Solution and Neighborhood.
For the first level program, denote the set of all nodes as  = [1, 2, . . ., ]  .After randomizing the sequence of the set, the first  nodes in the new set are selected as the set of hubs in the normal situation, denoted by .The neighborhood () is constituted by a swap of a hub node with a nonhub node.For  = 5,  = 2, and randomized set of nodes [2, 4, 1, 5, 3]  , the set of hubs can be selected as [2,4]  For the second level program, given the set of hubs  in normal situation, we can enumerate all the possibly disrupted cases in the worst-case scenario.Suppose that [1 2 3 4 5] is the current set of hubs and there are two hubs disrupted after attack.Then, there are 10 possible combinations of remaining set of hubs, that is, [1 2 3], [1 2 4], [1 2 5], and so forth.The optimization procedure in this level program can be realized by comparing the objective values of all the possible combinations.
The third level program is a facility assignment and path selection problem.The optimal solution can be obtained by the Floyd shortest path algorithm mentioned before.

Nondominated Solutions Set.
The two objectives of the first level program are sometimes conflicting.Generally, the optimization result for the program is not a unique solution but a nondominated solutions set (or Pareto optimal solutions set).Here, a solution,  1 , is said to dominate another solution,  2 , if  1 is not worse than  2 in all objectives and it is strictly better than  2 in at least one objective.A solution is said to be nondominated (or Pareto optimal) if it is not dominated by other solutions in the feasible solution space.All nondominated solutions for the multiobjective problem are called nondominated solutions set or Pareto optimal solution set.4) are used to be the fitness functions for tabu search procedure.

Fitness Function. The objective functions (3) and (
However, in each iteration, only one objective function is randomly selected as the fitness function for the current solution.This idea is motivated by the approach proposed by Kulturel-Konak et al. [28].

Move Rule.
The acceptance of a pairwise exchange between a hub node and a nonhub node is called a move.First, the candidate solutions in the neighborhood of the current solution are sorted by their value of fitness.According to the fitness, the best candidate that is not tabu, or if tabu it satisfies the aspiration criterion, is accepted as the next move for the iteration.
4.6.Tabu List, Tabu Tenure, and Tabu Frequent.During the iterations, the node exchanges are recorded in a short-term memory, called tabu list.The same exchanges are forbidden for a certain number of iterations, called tabu tenure.In this paper, the tabu tenure is set as  = √( − ).The update of the tabu list follows the rule of first in first out and later in later out.In order to avoid the local optimum, the search needs to restart from a different initial solution for several times.Long-term search memory is used to guide the new initial solution away from the nodes with highly visited frequency, called node frequency.

Aspiration Criterion.
The aspiration criteria for the tabu search are as follows: (1) a tabu move can be accepted if it generates a solution which is not dominated by any of the current nondominated solutions; (2) if all the candidate solutions are tabu and there are no candidates that dominate the current nondominated solutions, then the best candidate is selected as the new current solution.

Framework of the Heuristic Algorithm.
Based on what is mentioned above, the framework of the heuristic algorithm is as follows.
Step 1. Initialize Pareto set and tabu list as empty, set the Max Run times of the search and the Max Iteration numbers in each run.The search in each run terminates after "Max Count" numbers of nonimproving moves.The node frequency is set as Node Frequency = 0.
Step 2. Select  nodes randomly as initial hubs set  0 , solve the shortest path problem, and calculate the total network cost  0 1 in normal situation.For the hubs set  0 , enumerate all the possibly disrupted cases in the -disrupted hubs scenario.For each case, the shortest path problem is solved by the third level program (14), and the network cost for the case is obtained.The maximum cost of the cases is substituted into formula (4) and then we can get the total network cost  0 2 in the worst-case scenario.Add { 0 ,  0 1 ,  0 2 } into the Pareto solutions set.Set the initial hubs set as the current solution.
Step 3. Generate a random value  from the [0 1] uniform distribution.If  ≤ 0.5, choose formula (3) as the fitness function; otherwise, select formula (4) as the fitness function.Step 4. Generate the neighborhood set for the current solution.All solutions in the neighborhood are sorted according to their fitness.The first 10 solutions are selected as the candidate solutions.
Step 5.If the best candidate dominates some solutions in the Pareto set, turn to Step 9; otherwise, turn to Step 6.
Step 6.If all candidates are tabu, select the best candidate as the current solution; otherwise, select the best nontabu candidate as the current solution.Check the domination relationship between the current solution and the solutions in the Pareto set.If the current solution is not be dominated by any solution in the Pareto set, turn to Step 7; otherwise, turn to Step 8.
Step 7. Update the current solution, tabu list and tabu frequency.The iteration numbers Iteration + 1. Add the solution into the Pareto set, Count = 0. Turn to Step 10.
Step 8. Update the current solution, tabu list and tabu frequency.Increment the count of iteration numbers, Iteration + 1. Increment the count of consecutive nonimproving moves, Count + 1. Turn to Step 10.
Step 9. Take the candidate solution as the current solution and add it into the Pareto set.Remove the dominated solutions from the Pareto set.Update the tabu list and tabu frequency, Iteration + 1, Count = 0.
Step 10.If Iteration < Max Iteration and Count < Max Count, return to Step 3 for next iteration; otherwise, turn to Step 11.
Step 11.The search restarts from a new initial solution ignoring the nodes with high tabu frequency.Return to Step 2.
Step 12.If the run times reach the maximum limit, terminate the search and output the Pareto set.

Illustrated Example
5.1.Computational Data.Traditional researches on the hub location problem in transportation network usually use the instances from the AP data set [29], American CAB data set [24], or Turkish network data set [4], and so forth.However, these data sets only include the basic data of the distance and demand flow between cities. Besides, all of these data sets only consider the single transportation mode.Since there are no specific cases and data for the research on hub location problem in power projection network, the illustrated computational data in this paper are generated using simulation method.While generating the data, we set some rules in order to make the data conform to the realities of real world logistic networks.For example, the distance between nodes should satisfy the triangle inequality.The unit transportation cost between two nodes through a road link is typically higher compared to the rail link.The rail shipment transit cost typically exceeds that of road shipment.Parameter values in the model are set using the experiences and methods from related literatures [2,3].
Suppose that there are 15 projection nodes; the horizontal and vertical coordinates are sampled randomly from [0, 1000], as shown in Table 1.The spatial distribution of the projection nodes is shown as Figure 1.
The distance matrix is obtained using the Euclidean distance between nodes; that is,   = √(  −   ) 2 + (  −   ) 2 .Suppose that the ground transportation cost for is 0.6 ¥/t⋅km and ground transportation speed is 70 km/h.The railway transportation cost is set as 0.2 ¥/t⋅km and the railway speed is set as 60 km/h.The traffic demands between O-D pairs are generated from the uniform distribution [200,600].The traffic flow matrix is tabulated in Table 2.
The unit transit and operating cost for highway transportation between nodes are generated by sampling from a uniform distribution according to 5%-10% of the ground transportation cost, as shown in Table 3.The unit transit and operating cost of railway transportation are set as twice the cost of highway for the same nodes pair.The fixed cost of opening and operating a hub is generated by sampling from uniform distribution [100, 120]; the unit is ten thousand ¥, as tabulated in Table 4.
The hub-to-hub transportation cost discount factor using highway is set as  1 = 0.8, and the factor using railway is set as  2 = 0.6.

Computational Results and Analysis.
Suppose that there are  = 5 hubs to be located in normal situation and there would be  = 2 hubs disrupted in the worst-case scenario.The test was run on an Intel Core i5 2.2 GHz and 8 GB RAM computer.The tabu tenure was set as 7, and the tabu frequency was set as 5.The relationship between the size of the Pareto set and the number of iterations in a run of tabu search procedure is shown in Figure 2. As Figure 2 shows, the size of the Pareto solutions set no longer changes after 43 iterations, indicating that the algorithm converges.Computational experience herein shows that the search procedure converges after 50 iterations in each run.Thus, we set the maximum iterations in a run of tabu search as 50.In order to avoid the local optimum, the search procedure restarts from a new initial solution for 20 times.The computational problem was solved after 25.4 s CPU time.There are 6 nondominated solutions in the Pareto set, as tabulated in Table 5.
As shown in Table 5, for solution 1 which selects [1 5 8 10 14] as the hubs set, although this instance has the minimum network cost in normal situation, it has the maximum network cost in the worst-case state.In other words, the network has a poor resilience for solution 1.For solution 3 which locates [1 5 8 9 14] as hubs set, it is 10.8% better than solution 1 in the worst-case cost but only 1.7% worse in the normal cost.Solution 3 has the maximum efficiency-cost ratio.Solution 6 which chooses [3 5 8 9 11] as the hubs set has the best network resilience   and the minimum cost in worst-case scenario.However, it imposes on the network an excessive cost in normal situation.Therefore, decision-makers can select different hub location strategies from the Pareto set according to their specific requirements.The tradeoff curve for the hub location strategy is depicted in Figure 3.In general, solution 1 is suitable for the optimum decision-making and the case of lighter threat from the enemy.Solution 6 is suitable for the pessimism decisionmaking and the case of serious threat from the enemy.Since solution 3 is much better than solution 1 in terms of the worstcase cost but not much worse in terms of the normal cost, it is a more reasonable hub location strategy for the design of a resilient power projection network.

Conclusion
The ability of power projection is a crucial element of a state's military power.The projection procedure relies on safe and efficient transportation network.In this paper, we focus the issue on the hub location problems in design of power projection network from a resilience perspective.Not only the cost of the projection network in normal situation is considered, but also the performance in presence of hubs disruptions is taken into account.A biobjective and trilevel integer programming model is proposed for the problem in view of game theory.The model can help to make tradeoffs between the performances of the projection network in different situations.A heuristic approach based on multiobjective tabu search and Floyd shortest path algorithm is developed to solve the model.Computational example shows that the results can provide some useful references for the hub location decision in design of a resilient power projection network.
The intermodal transport in power projection in this paper only considers the modes of road and rail.More transport modes can be taken into account in the future research.In our analysis, it is assumed that the enemy may select the hubs for attack which result in the greatest damage to the network.It is worthwhile to extend the model to cases in which the hubs disruptions following a probability distribution.

Figure 1 :
Figure 1: Spatial distribution of the projection nodes.

Figure 2 :
Figure 2: Convergence curve of the tabu search procedure.

2 Figure 3 :
Figure 3: Tradeoff curve for the hub location strategy.
: traffic demand from origin node  ∈  to destination node  ∈ , where   = 0,   : fixed cost of opening and operating a hub at node ,    : unit transportation cost from node  to node  using transportation mode , where    = 0, ĉ  : drayage and operating cost of one unit transportation via hubs  and , where ĉ  = 0,    =  1  +ĉ   +     + 1  : unit transportation cost from origin node  to destination  via hubs  and  using mode , Parameters.Notations employed in the model are introduced as follows:  = {1, 2, . . ., }: set of all nodes, : set of hubs, : number of hubs, : number of disrupted hubs, where  < ,  = {1, 2}: set of transportation modes between hub pairs, where  = 1 for road transportation and  = 2 for railway transportation,

Table 1 :
Serial number and coordinates of the projection nodes.

Table 3 :
Transit and operating cost for highway transportation between nodes (unit: ¥/t).

Table 4 :
Fix cost of opening and operating hubs (unit: ten thousand ¥).

Table 5 :
Pareto solutions set for the computational example.