Managing facility disruption in hub-and-spoke networks: formulations and efficient solution methods

Hub disruption result in substantially higher transportation cost and customer dissatisfaction. In this study, first a mathematical model to design hub-and-spoke networks under hub failure is presented. For a fast and inexpensive recovery, the proposed model constructs networks in which every single demand point will have a backup hub to be served from in case of disruption. The problem is formulated as a mixed integer quadratic program in a way that could be linearized without significantly increasing the number of variables. To further ease the model’ computational burden, indicator constraints are employed in the linearized model. The resulting formulation produced optimal solutions for small and some medium size instances. To tackle large problems, three efficient particle swarm optimisation-based metaheuristics which incorporate efficient solution representation, short-term memory and special crossover operator are proposed. We present the results for two scenarios relating to high and low probabilities of hub failures and provide managerial insight. The computational results, using problem instances with various sizes taken from CAB and TR datasets, confirm the effectiveness and efficiency of the proposed problem formulation and our new solution techniques.


Introduction
The classical hub location problem aims to find suitable locations for hubs, allocate the existing demand to these facilities, and direct the flow between origin-destination pairs at minimum cost. In this class of locational problems it is common to assume there is a link between every hub pair, there is no direct path between non hub nodes, and there is economies of scale for using the inter-hub connections (Alumur and Kara 2008). Hub location problems are categorised into two distinctive groups namely single and multiple allocation problems. In the former, all incoming and outgoing traffic to and from every node is transferred via a single hub, while in the latter each node is allowed to receive and send flow through more than one hub. The focus of this paper is on the Uncapacitated Single Allocation p-Hub Location problem where the number of hub facilities to open is predetermined. This problem is also known as the Single Allocation p-Hub Median problem.
Early research on hub-and-spoke systems focused on developing efficient mathematical formulations and solution techniques for this paradigm. O'Kelly (1986O'Kelly ( , 1987 presented the first mathematical model for the single allocation p-hub median problem. Another pioneer of the hub location research, Campbell (1994), developed a linear integer formulation for the problem. Based on the idea of multicommodity flow, Ernst and Krishnamoorthy (1996) proposed a new set of formulations for both single and multiple allocation cases. A widely used formulation in the literature for single and multiple allocation p-hub median problem is the work of Skorin-Kapov et al. (1996). The MIP models proposed by Skorin-Kapov et al. (1996) yield to tight linear relaxation.
Over the years, a number of solution techniques e.g., approximation and exact methods have been developed and successfully applied to solve instances of the hub locationallocation problem. Examples of such methods include Simulated Annealing (Ernst and Krishnamoorthy 1999), Genetic Algorithm (Kratica et al. 2007), Hybrid GA and Tabu Search (Abdinnour-Helm 1998), Lagrangian Relaxation (Elhedhli and Wu 2010;Contreras et al. 2009), Benders Decomposition (Contreras et al. 2012), Branch and Bound (Ernst and Krishnamoorthy 1998), and Branch and Cut (Yaman and Carello 2005)

among others.
More recent studies focus on the extension of the classical version and aim to develop more realistic hub location-allocation problems such as models with congestion (e.g., Elhedhli and Wu 2010;Azizi et al. 2017), flow dependent economies of scale (Camargo et al. 2009;Campbell et al. 2005), and competitive hub location models (Eiselt and Marianov 2009). For an overview of hub location problem the reader may refer to Campbell and O'Kelly (2012) and Contreras (2015).
In traditional approaches to network design, supply/facility disruption is often presumed to be a rare event. These approaches are mainly concerned with the location of the facilities and the allocation of demand to these facilities in such a way to minimise the total network cost. In practice, however, it is likely that over the course of time one or even more than one facility become disrupted. The cause of disruption may vary from severe weather condition and natural disasters to labour dispute and man-made sabotage. While a set of disruption causes (e.g., earth quake) may occur less frequently, a large number of other causes are likely to strike a supply chain network. For instance, to deal with frequent supply disruption Wall-Mart uses an emergency operations centre. The centre is responsible to prepare for and mitigate the effects of man-made and/or natural disasters (Snyder et al. 2016).
Research have shown that even small disruption in supply may have devastating impacts. In 1998, strikes at two General Motors parts plants resulted in closure of 100 other parts plants and then 26 assembly plants. In the context of locational problems, facility disruption result in excessive transportation cost as customers initially served by these facilities now must be served by other operational facilities in the network (Snyder and Daskin 2005). Azad et al. (2012) proposed a capacitated supply chain network design model under random disruptions. The formulation aims to determine the locations and types of distribution centres and allocate customers to each opened facility. Azad et al. (2012) considered partially disrupted facilities which can operate with a fraction of their initial capacities. The authors considered a "general supply chain network" which is different from the topology of the huband-spoke system. In a hub-and-spoke network the flow is processed at least in one facility, is bidirectional and there are inter-hub flow interactions.
Hub location problems with reliability consideration have received very limited attention in the literature (Kim and O'Kelly 2009). Kim and O'Kelly (2009) formulated two p-hub location problems in telecommunication network with reliability consideration. The first model, p-hub maximum reliability, maximises the expected network flow. The second model, p-hub mandatory dispersion, specifies the optimal hub locations to increase network reliability and prevents excessive concentration of interacting flows in particular hub facilities. Kim and O'Kelly (2009) investigated both single and multiple allocation cases without considering backup hubs and rerouting the flows. An et al. (2015) proposed a mixed integer nonlinear formulation for reliable hub-andspoke problem. The authors successfully solved the formulation for problem instances up to 25 nodes using branch and bound and Lagrangian methods. An et al. (2015) compared the number of passengers served in a classical and that in a reliable network. Using a numerical example, they show that a reliable network can transport more passengers by its regular routes than a network with classical configuration. They also showed that the configuration of classical and reliable networks may not be identical. An et al. (2015) assumed (a) the flow transported in the network is symmetrical and (b) in designing alternative routes a multiple allocation is used (regardless of single allocation or multiple allocation structure used to determine regular routes). Furthermore, in the proposed model by An et al. (2015) flow transported between a pair of nodes which are allocated to the same hub facility will be rerouted through a single backup hub i.e., restricting a pair of demand points initially assigned to the same hub to be rerouted via a single backup hub. Tran et al. (2017) proposed a mixed integer nonlinear model to find the optimal location of predetermined number of hubs. They linearize the model using a specialized flow network called a probability lattice and used a tabu search algorithm to solve problem instances up to 25 nodes with symmetrical flows. Tran et al. (2017) used a level-set approach to allocate spokes to hubs in a specific order. This method of allocation though effective may prevent the model to be extended by e.g., incorporating capacity constraints and/or adapting multiple allocation scheme. Azizi et al. (2016) also proposed a mixed integer nonlinear formulation to design hub-and-spoke network taking into account the probability of hub failure and re-routing cost. They solved a number of relatively small problem instances to optimality and proposed an evolutionary algorithm to solve large instances. Similar to An et al. (2015), the authors showed that configuration of networks with and without hub failure consideration may not be the same. Azizi et al. (2016) consider a backup facility for each hub and assumed that once a facility becomes unavailable all demand points initially served by the disrupted facility is served by one of the operating facilities in the network. Reallocating the affected nodes to only one of the operating facilities in the network though effective in cases where demand points need to communicate via a single hub, may not lead to optimal solution in other cases.
In this paper, we present models that simultaneously minimises the day-to-day network operations cost and the expected transportation cost that incurs when one of the hub facilities experience short term disruption. In doing so, our proposed models seek to find the optimal locations for the hub facilities, allocate the demand to these hubs and find the least expensive backup facility for each demand point in the network. The resulting network(s) is expected to satisfy the demand with a minimum additional cost in the event of supply disruption. Our models intend to provide robust solutions (i.e., networks) that perform well even when part of the network fails. Handling instances with and without symmetrical flow, reallocating : Affected hub : Operational hub : Affected Inter-hub routes : Affected node-hub route : Re-assignments To further elaborate our approach, consider the hub-and-spoke system with 10 nodes and 3 hubs presented in Fig. 1. In such a network if any of the three hub facilities is disrupted, either one or the two other operating hubs in the network could be used as backup facility (ies) to satisfy the demand initially served by the disrupted facility. For instance, if hub 5 is disrupted then to maintain network operations, demand point 1 could be reallocated to either hub 9 or 7. Similarly, demand point 10 could also be reassigned to one of the two operating facilities as depicted in Fig. 1.
In short, the contribution of the paper include: (i) An investigation into a reliable hub-and-spoke system incorporating heterogeneous probability of hub failure (ii) A non-linear model purposely designed to simultaneously construct networks with and without symmetrical flow and select backup facilities for every demand possibly affected by a disruption (iii) An improved linear model with the use of indicator constraints and (iv) Efficient PSO-based meta-heuristics for large instances.
The remainder of this paper is organised as follow. In Sect. 2, we describe the problem and present non-linear, linear and improved linear formulations. The details of the three proposed PSO based algorithms are presented in Sect. 3 followed by computational results and discussions in Sect. 5. Concluding remarks are presented in Sect. 6.

Problem statement and formulation
In the single allocation p-hub median problem, each node is assigned to a particular hub and all incoming and outgoing flow are routed through that hub. The total number of nodes in the network is assumed to be N , each node is considered a potential site and the number of hubs to open is p. A path from origin spoke i to destination spoke j includes three parts: collection from spoke i to the first hub k, transfer between first hub k and the second hub m, and distribution from hub m to destination j. The cost per unit flow along this route, i → k → m → j, is calculated as χ ×c ik +α×c km +δ×c m j where χ and δ are coefficients of collection and distribution respectively and α is the inter-hub discount factor. Let λ i j be the amount of flow to be routed from origin i to destination j, the transportation cost from i to j routed via hubs k and m, C ikm j is then calculated as C ikm j = λ i j χ × c ik + α × c km + δ × c m j .

Notation
Our model incorporates the following decision variables: hub location and allocation variable z, the route selection x, and the backup selection variable u. These variables are defined as follows: z ik = 1 if node i is assigned to hub k and =0 otherwise; x ikm j = 1 if flow from i to j passes through hub k and m and = 0 otherwise; u ikn = 1 if n is the backup for node i when hub k in the path from origin i to any destination j is disrupted and = 0 otherwise.
We consider location-specific probability of hub failure (heterogeneous). Each hub k (or m) has a probability of failure q k (q m ) and we assume that only one hub will be disrupted at any given time.

Model formulation
To formulate the problem, in the event of any hub disruption we recognise four different types of flows namely: (1) flow that is not affected by the disruption (2) flow that uses the disrupted facility as its first hub (3) flow that uses the disrupted facility as the second hub (4) flow that originated from or destined to a hub facility.
The cost of transporting the flow in part of the network which is not affected by disruption could be calculated by subtracting the transportation cost of those routes affected by hub failure from the total network transportation cost as follows (Azizi et al. 2016): Alternatively, this cost could be calculated by adding up the transportation cost in all routes for which the disrupted hub is neither their first nor their second hub. In other words, the total cost of interactions between nodes in the right hand side of the network in Fig. 1 namely nodes 2, 3, 7, 9, 8, 6 and 4 gives the transportation cost in part of the network that is not affected by disruption. This is presented as where Index l represent the assumed disrupted facility. We conducted a computational experiment and found that formulating this cost as a function of R l and Q l may result in a weak LP relaxation which is likely to lead to a lengthy computational time. Therefore, in our formulation Eq.
(1) is used to calculate the cost of transporting the flow in part of the network which is not affected by disruption.
In case of hub disruption, we assume the affected demand points are served by other facilities in the network. We denote ikn as the rerouting cost of the flow through backup facility n when hub k in the path from origin i to destination j fails. ikn is expressed as ikn = m =k j C inm j x ikm j ∀i, k, n; i = k (3) Figure 2 illustrates how a flow is rerouted via backup facility n when the first hub k is unavailable. The new route to transfer the flow from origin i to destination j is highlighted by solid lines. Similarly, β jmn is the cost of rerouting the flow through backup n if hub m fails (the first hub for the route initiate at node j), see Fig. 3.
We formulated the problem as a bi-objective optimisation problem and the two components (F1 & F2) of the objective function of t he model are presented as follows: The costs F1 and F2 refer to the transportation cost in a regular situation and the expected transportation cost resulted from hub failures respectively. Objective F2 consists of five terms. The first term computes the transportation cost of the flow that is not affected by disruption. The second term in F2 calculates the rerouting cost of the flow originated from node i through backup facilities when hub k in the path from origin i to destination j is disrupted and the path includes two hubs; the third term calculates the rerouting cost when the only hub in the path from origin i to destination j is disrupted. In a path with a single hub, origin (i) and destination ( j) are initially assigned to the same node (i.e., k). When hub k is disrupted, the two demand points (i and j) would be either allocated to the same or to different backup facilities. The forth term in F2 calculates the rerouting cost of the flow through backup facilities when hub m is disrupted in the path from i to j. In this study it is assumed that when a hub is disrupted, it cannot send or receive any flow to or from other nodes in the network. Therefore, the last term in the objective function F2 (i.e., the fifth term) penalizes the loss of flow/demand in disrupted situations where the source or destination of the flow is a hub. In our study, the penalty cost of losing a unit flow, ϕ i j , is considered twice as much as the transportation cost of a unit flow between origin i and destination j (Azizi et al. 2016). The proposed model minimizes the weighted objective function w F1+(1 − w) F2 where 0 ≤ w ≤ 1. The first formulation of the Reliable Single Allocation p-Hub Location problem which is denoted by (RSApHL-I) for short is presented as follows Subject to: Constraints (8)-(12) are classical constraints for the single allocation p-hub location problem (Skorin-Kapov et al. 1996). Constraint (8) ensures every node is assigned to only one hub whereas Constraint (9) limits the number of hubs to be opened to " p". Constraint (10) is to ensure node i is assigned to an open hub. Constraints (11) and (12) guarantee that all the traffic between an origin-destination pair has been routed via the hub sub-network. Constraint (13) guarantees each node i has only one backup (n) and it differs from the hub initially assigned to (i.e., k). Constraint (14) ensures the candidate backup n for the affected node i is an open "hub" facility. Equality constraint (15) calculate the transportation cost of the flow that is not routed via hub l. Constraints (16) compute the rerouting cost of flow through backup facility n that is originally transported via the first hub k. Constraints (17) calculate the rerouting cost through backup n that is originally transported through hub m. Constrains (18) and (19) are the standard integrality constraints. The proposed formulation for RSApHL-I is a Mixed Integer Quadratic Programing (MIQP) and includes n 2 + 4n 3 + n 4 variables. In the next section we will show that because of our novel modelling approach linearizing the model using conventional techniques will not significantly increase the number of variables. Our approach is similar to the technique proposed by Chaovalitwongse et al. (2004) to linearize multi-quadratic 0-1 programming problems.

Linearization
The non-linear terms in the above objective function resulted from the multiplication of binary and non-binary variables of x, u, z, and γ i jl , ikn , β jmn . To linearize the model, these terms are substituted by three continuous and binary variables y i jl , ikn , ψ ikmnj and ξ jmn as follows.
To enforce the above relationships and using the Big-M idea the following sets of constraints are added to the formulations. Here, we refer to Big-Ms as M 1 , M 2 and M 3 .
ikn ≤ ikn ∀i, k, n The resulting linear model for Reliable Single Allocation p-Hub Location problem which we call (RSApHL-II) is presented as follow.
Subject to: The number of variables in the above MIP model increases to n 2 +7n 3 + n 4 + n 5 . A summary of the number of variables and constraints before and after the linearization is given in Table 1. RSApHL-II (MIP) n 2 + 7n 3 + n 4 + n 5 1 + n + 2n 2 + 16n 3 + 4n 5

Linear model with indicator constraints
In MILP formulations, it is common to model constraints that either hold or are relaxed depending on the value of a binary variable. These constraints are often modelled using the so called Big-M formulations to enforce the relationship between a new variable and its corresponding nonlinear term. These Big-M formulations though useful are often blamed for unstable behaviour in MILP model. Indicator Constraints have the advantage of preventing the problems associated with Big-M formulations. Similar to Big-M formulations, indicator constraints use a binary variable to turn on or turn off the enforcement of a constraint but without using a big M. Further information about the indicator constraints could be found in the paper by Bonami et al. (2015). In our linear formulation i.e., RSApHL-II we replaced all the constraints formulated using Big Ms by indicator constraints and set 0 ≤ x ikm j ≤ 1. For instance, using indicator constraints the following can be written logically as The new formulation with indicator constraints which we refer to as (RSApHL-III) is given below for the sake of simplicity and convenience.

Preliminary study
To evaluate and compare the performance of the proposed (linear) formulations with and without indicator constraints, we conducted a limited computational experiment as described in the following. A relatively small problem instance with 10 nodes and 3 hubs is selected from the U.S. Civil Aeronautics Board which is known as CAB dataset (O'Kelly 1987). We generated three instances of the original problem by setting the objective function weight w to 0.3, 0.5 and 0.7. The discount parameter α is set to 0.2; the coefficients of collection χ and distribution cost δ are set to one per unit for all three instances. The problems are solved using CPLEX 12.6 with default options values. The algorithm is run for 5000 s on Intel Core i5 PC with 3.2 GHz processor with 4 GB of RAM. The computational results are presented in Table 2. The results pertaining to RSApHL-II model show that CPLEX is unable to produce the optimal solution for the problem instance with w = 0.3 in 5000 s but manages to locate the optimal solution for the other two instances with w = 0.5 and 0.7.
For the RSApHL-III model, results presented in Table 2 clearly show that the same problem instances could be solved more effectively and with significantly less computational efforts when all Big-M formulations are replaced with indicator constraints. All optimal solutions are now discovered while requiring between 20 and 25% of the (computational) times reported for instances of RSApHL-II. This is an interesting result which is achieved by a combination of an effective formulation and the use of the recently made available modelling tool i.e., indicator constraints.
This improved formulation will be used to provide when possible lower and upper bounds for benchmarking purposes. As will be shown in our computational results, the above improved model formulation could only be used to solve relatively small to medium size problem instances to optimality.
Due to the limitation in solving the above model formulation to optimality for realistically sized problem instances with CPLEX, one way forward is to design an efficient metaheuristic. In this study, we developed three metaheuristics based on the well-known evolutionary algorithm of Particle Swarm Optimisation (PSO). These algorithms are discussed in the following section.

Particle swarm optimization
Particle swarm optimization (PSO) is one of the most efficient nature-inspired optimization algorithms. It was proposed by Kennedy and Eberhart (1995) about two decades ago. The algorithm was developed based on the motion of a flock of birds searching for foods and aimed to mainly optimise continuous non-linear optimisation problems.
In PSO based algorithms, at the beginning of the evolutionary process a set of particles which is called a swarm is generated randomly. Each member of a swarm is called a particle. During the search process, particles are flown through a hyperspace. A particle's status on the space is characterized by two elements namely position and velocity. Particles may change their position in the search space similar to flying birds searching for food in the sky. While searching the hyperspace, a particle adjusts its (new) velocity according to its own best experience, the best experience of all particles in the swarm and its previous velocity. The particle's new velocity and current position are then used to determine its new position.
A large portion of the published research works on PSO is dedicated to continuous optimization problems whereas there is a lack of research dealing with its discrete counterpart namely combinatorial optimization problems (Liu et al. 2008). To the best of our knowledge, our proposed algorithms is the first PSO-based optimisation technique to be developed for solving this interesting though challenging class of hub location problem.
In this paper two variations of a basic PSO and a hybrid PSO algorithm that incorporate crossover and memory are proposed. The two versions which we refer to as PSO-v1 and PSO-v2, differ in terms of solutions representation and the mechanism by which a continuous particle position is transformed into one or more discrete solution(s). The hybrid PSO algorithm on the other hand, is an extension of PSO-v2 which incorporates a short-term memory and a crossover operator, both of which are in our knowledge, new ingredients that are successfully embedded for the first time into the PSO search methodology.

PSO-v1 algorithm
In PSO-v1 algorithm, each particle is represented by a (1 × n) array in which n represent the number of nodes. We define the position of a particle by With a (i) kl being its ith position value. We represent a group of m particles (i.e., a population of particles) at each iteration l by Particle position PSO-v1 algorithm starts by generating random continuous position values for the m particles in the population in the range where p is the number of hubs and 0 < ε < 1; in this study, we assume ε = 0.01. Figure 4 illustrates a typical particle position for a problem with 10 nodes and 3 hubs which is populated with continuous random numbers in the range of [1, 3.99].
Particle velocity In PSO-based algorithms, particles fly to new areas in the search space using a velocity function. In this study, the initial velocities are generated randomly according to the following expression: v k0 is the initial velocity of particle p k0 with v min and v max being the pre-defined range of velocity values. We assume v min is zero and we define v max as a function of the number of nodes (n) and hubs ( p) in the problem where θ is a small value between 0 < θ < 1; here we assume θ = 0. Updating particle velocity and position At iteration l the kth particle velocity is updated using the following formula v k, l+1 = ωv kl + C 1 r 1 (Pbest kl − p kl ) + C 2 r 2 Pbest gbest,l − p kl ; k = 1, . . . , n (54) where C 1 and C 2 are acceleration coefficients; ω the inertia factor; and r 1 and r 2 are two independent random numbers uniformly distributed in the range [0 , 1]. Once all the n particles new velocities are determined, the corresponding particles' positions are then updated by adding the new particles velocities to their current positions as follow Also it is common practice in PSO algorithms to control the search by limiting the magnitude of particle's elements. In this study we set Transforming a continuous particle position into one or more discrete solutions Due to the continuous nature of the particle position in PSO, standard encoding schemes cannot be directly adopted for discrete location problems such as hub and spoke systems. One procedure to transform a continuous particle position into a discrete solution in PSO-v1 is described via the following example. Consider the particle' position presented in Fig. 5. As mentioned earlier, the number of columns/elements in a particle position array corresponds to the number of nodes in the problem (in this example the number of elements are 10). The integer part of the largest element of the array refers to the number of hubs which in the above example is 3. The one dimensional array with continuous values presented in Fig. 4 could be encoded as an array with integer values by taking its integer part [a n ] where a n represent the value of each element in the position array. The corresponding particle position array with integer values for the example in Fig. 4 is presented in Fig. 5. All elements in the particle position that have the same value constitute a cluster. For example, the following three clusters {3, 7}, {2, 5, 10} and {1, 4, 6, 8, 9} are presented in Fig. 6a.
As each member of a cluster (i.e., a node) is a potential site for a hub facility, a discrete feasible solution can be generated by connecting the selected hubs. For example, a possible configuration is given in Fig. 6b where nodes 3, 2, and 4 corresponding to clusters 1, 2 and 3 are chosen. Considering the combination of potential hub facilities in each cluster node, a quite large number of feasible solutions could be constructed from one single particle position. In the example presented in Fig. 6a, a total of 30 (3 × 2 × 5) solutions could be generated. In general, if n i refers to the number of nodes in cluster i, the total number of configurations become p i=1 n i . Clearly, as the number of hubs and/nodes increases the number of solutions/possibilities could increase significantly. Examining all possible solutions produced by changing hub facilities in each cluster could be, particularly, for relatively large instances computationally expensive. Therefore, it is more practical to either randomly select one solution among all  Fig. 6 a Three cluster nodes associated with particle array in Fig. 6, b a typical solution generated from the three cluster nodes Start i=1, initialise K0 (say K0=0.5) Do while i < SwarmSize create a one-dimensional particle position array with continuous values Take the integer part of the particle elements (i.e., ) Identify the p cluster nodes in the particle position array Generate K solutions from potential solutions; Evaluate the cost for each configuration and select the solution with the least cost. Apply a perturbation to increase diversification Fig. 7 Algorithmic steps to convert a continuous particle position into a discrete solution  possible solutions or generate a fraction of potential solutions and select the best one with the lowest total cost.

Loop End
We set up a limited number of experiments using instances with 10, 15, 20, and 25 nodes and concluded that examining a fraction of the potential solutions (e.g., k 0 ) and selecting the best one among the subgroup could provide good overall best results.
In our proposed PSO-v1 algorithm, once a discrete solution is found, a small perturbation (e.g., reassigning a node to another hub) is also applied to improve the diversification in particles population. The main steps to transfer a continuous particle position into a discrete solution are given in Fig. 7 with Swarmsize representing the population size.
Backup selection PSO-v1 encoding scheme transfers a continuous particle position to a discrete solution to a single allocation median hub location problem. Although the selection process of backup facilities could also be a part of the solution representation, to retain the simplicity of our solution representation we opted to treat this decision separately. Backup facilities for the existing demand points could be chosen (1) randomly (2) by examining all possibilities or (3) based on proximity to a potential backup facility. Examining all combinations (i.e., complete enumeration) will provide the best solution as it guarantees the selection of the most economically suitable backup but it is computationally expensive even for small problem instances. We therefore compared the two other options using a number of test problems and concluded that in large instances selecting nearby backups to demand nodes outperforms the random backup selection strategy. For the small to medium size problem instances, random backup selection may result in better quality solutions.

PSO-v2 algorithm
In this variant PSO-v2, a continuous particle position is directly encoded/mapped to only one discrete solution. The main difference(s) between this version and the earlier one are the solution representations and the mechanisms through which a continuous particle position is encoded into a solution.
Particle position and representation Here, a particle is represented by a matrix of size 2 × n where n stands for the number of nodes. Entities in the first row of the matrix represent hub facilities. Therefore, in a problem with p hubs, we deal with only the first p elements in the first row of the matrix and the remainders are used to rank the entire elements later. For instance, see Fig. 8 for a problem with 10 nodes and 3 hubs.
Entities in the second row of the matrix denote the allocation of each node to one of the hubs specified in the first row.
To generate an initial particle position, the first row of the matrix is populated with continuous numbers generated using the following equation we assume p min = 0 and p max = p × n. Similar to PSO-v1 position array, the second row of the matrix is filled with randomly generated continuous numbers in the range of [1, p + 1 − ε] where 0 < ε < 1.
Particle velocity Similar to that in PSO-v1, the initial values of particles' velocities, v k0 , are generated using Eq. (52) with v max and v min defined as in PSO-V1 (Eq. 53). Particles velocities and positions are updated in each iteration using Eqs. (54) and (44) respectively.
Transforming a continuous particle position into a single discrete solution To transfer a continuous particle position into a discrete solution, first all entities in the first row of the particle position matrix are ranked. Then p elements of the first row are considered as potential hub locations. For instance, in the particle position presented in Fig. 8 the smallest number is the fifth element "0.7" and is ranked "1"; the third element "1.1" is ranked second (i.e., "2"); and the first element "1.6" is ranked third (i.e., "3") etc. After ranking the numbers in the first row of the particle position matrix, the three potential hub facilities ( p assume to be 3) are located as node 3 (first hub), 5(second hub) and 2(third hub).
To determine the allocation of non-hub nodes to the hubs in the first row, we consider the integer parts of the elements in the second row in Fig. 8.   1.6 2.5 1.1 3.6 0.7 5.6 2.1 7.1 3.4 4.2 3.1 2.95 1.3 3.2 2.7 3.6 1.6 3.3 3.4 2.1

Start
i=1 Do while i < SwarmSize Take a two-dimensional particle position matrix with continuous values Rank all elements in the first row Consider the first p elements in the first row as potential hub locations Take the integer parts of the numbers in the second row Evaluate the solution If infeasible, repair the solution Fig. 11 Steps to convert a continuous particle position into a discrete solution:  For instance, in Fig. 9 the third and the seventh elements of the second row are identical and have the value of "1" which means nodes 3 and 7 are allocated to the first hub which is hub "3".

Loop End
Repair algorithm During the encoding process of PSO-v2 occasionally an infeasible solution may be generated as a given hub may not be necessary allocated to itself. In such situation, a simple repair algorithm to amend the solution and change it into a feasible one is used. For instance, in Fig. 9, the second node (node "2") has been assigned to the second hub (hub "5") instead of being allocated to itself (i.e., node "2") which is the third hub in the first row. To repair this solution, the value of the second element (of the second row) is simply replaced by value "3". The resulting repaired solution is given in Fig. 10.
The transfer of a continuous particle position into a discrete solution is presented in Fig. 11.

Hybrid PSO algorithm
Our preliminary results for relatively small problem instances show that both variants (i.e., PSO-v1 and PSO-v2) perform well in finding good quality solutions to the RSApHL-III.
Nevertheless in the majority of cases the same solutions are obtained by PSO-v2 requiring relatively shorter computational times.
Research have shown that PSO-based optimization algorithms are highly randomised and, in general, are susceptible to early convergence (i.e., could be trapped into a local optima too early in the search). To overcome this drawback and hence improve the performance of the PSO-2 algorithm, two new features are imbedded into the algorithm to make our hybrid PSO which we refer to H-PSO for short. (a) The first feature is a memory which is setup by collecting a number of global best particle positions at the end of each iteration. The size of the memory is assumed to be the same as that of the PSO swarm size e.g., 20. (b) The second feature is a crossover operator aiming to improve the quality of the solutions in the PSO population as well as diversifying the search. Once a memory is completed, solutions stored in the memory are combined using a special type of crossover operator. To perform the crossover operation, two parent particles are selected each time randomly to produce an offspring particle. The fitness of the newly generated solution is then compared with those in the PSO population; the offspring replaces the first inferior particle found in the population. The offspring particle is ignored if it cannot replace any of the particles in the population. Upon completion of the new population the memory is reset before the next iteration of the algorithm begins. The crossover operator is described next.
The proposed crossover operator The solution representation used in H-PSO is the same as that in PSO-v2 (2 dimensional matrix). Classical crossover operators are often designed for one-dimensional arrays and could not be directly applied here. Therefore, we propose a special type of crossover that combines the classical "single point" with what we name a "constructive" crossover operator.
To produce a new offspring, a template 2-dimentional matrix for the offspring solution is first constructed. Next, two parents are marked randomly from the current population. Then, one of the two parents is selected randomly and the first element on the "first row" of its matrix is transferred to the same place on the offspring particle template matrix. The second element of the offspring template is taken from the other parent. These steps are repeated until the first row of the offspring array is populated with elements transferred from the two parents.
Upon completion of the first row of the offspring particle matrix, a classical single-point crossover is applied to the second row of the parents' particles to populate the second row of the offspring particle matrix. In a classical single-point crossover operation, one point is selected on the second row of the parent particle matrix. One of the parents is then selected randomly and all elements up to that point on its array are transferred into the offspring matrix. The rest of the locations on the offspring matrix (i.e., second row) are filled by the elements beyond the crossover point on the other parent chromosome. The steps of the proposed crossover operator are presented in Fig. 12.

Computational results and analysis
We tested our proposed PSO-based heuristics on 126 benchmark problems and solved our formulation for 54 of these instances using CPLEX. These instances are derived from U.S. Civil Aeronautics Board (CAB) (O'Kelly 1987) and Turkish Postal System (TR) datasets (Çetiner 2003). The problem instances are generated by setting the number of nodes N to 10, 15, 20, 25 and 55; the number of open hubs " p" to 3 and 5; the discount factor α to 0.2, 0.4 Single point-constructive crossover operator i=1 While +1 Do create a template matrix for an offspring particle select two particles from the current swarm j=1 select one of the two particles while j<p+1 Do Transfer item j from the first row in the particle matrix to the same place in the offspring particle matrix select the other particle End while Apply the classical single-point crossover operation to the second rows of the two particles and complete the offspring matrix End while Fig. 12 The proposed crossover operator for H-PSO and 0.8; and the objective weight w to 0.3, 0.5, and 0.7. The coefficients of collection χ and distribution cost δ are set to 1 per unit for all test problems (χ = δ =1). The PSO parameters are as follows. The swarm size is set to 25 for test problems with 10 nodes and 20 for the rest of instances. C1 and C2 are set to 2 and the inertia factor, ω, to 1 for all test problems. The computational times are presented in Table 3. Each test problem is run 10 times and the best results are reported. The MIP version of RSApHL-III model is coded in AIMMS and solved using CPLEX 12.6 with CPLEX options set to their default values. All algorithms are run on an Intel Core i5 PC with 3.2 GHz processor with 4 GB of RAM. In this section we analyse two scenarios reflecting real life applications. In the first one we consider networks that are quite highly susceptible to hub disruptions. Whereas in the second, we attempt to address other networks in which hub disruption is less likely to occur. The former scenario includes a network in which hubs are assumed to be, for instance, distribution centres. In these types of networks the probabilities of hub failures are expected to be relatively high. In the second scenario, such as hub-and-spoke networks in airline industry, the chance of hub failure is believed to be relatively low though extremely costly if that happens.
For the case with high probability of hub failure, the probabilities, q i , for all nodes in the network i = 1 . . . n are generated randomly from U [0.1, 0.3]. For the other case, we use the low probabilities of hub failures as given in the paper by An et al. (2015).

CAB dataset
In CAB dataset the flow is transported equally in both directions (i.e., symmetrical flow). In such networks, the size of the problem can be reduced by calculating the transportation cost in just one direction and then doubling the cost to present the total cost of transporting the flow. This setting has been adopted by a number of studies and, to some extent, reduces the computational time required to solve problems with symmetrical flow. We begin our analysis by evaluating the computational performance of the proposed PSO algorithms and the CPLEX using 18 small benchmark problems with 10 nodes and 3 and 5 hubs. A summary of the (average) results is presented in the first row of Table 4. CPLEX solved all 18 benchmark problems in a reasonable computational times with 17 within 1-11 min and the other requiring just 18 min with a large portion of this time being spent in closing the last 5% gap. The average computational time to solve all 18 problems is 372.5 s. Nevertheless these are still significantly high when compared against the time required by any of the PSO-based algorithms.
The results in Table 4 further show that the three PSO-based algorithms performed well though H-PSO found more optimal solutions than the other two algorithms (i.e., . For test problems with 10 nodes, the optimal locations of hubs and allocations of non-hubs to the selected hub facilities as well as backup facilities for each node are presented in Table 5. For example, in this table the optimal value of the objective function for the problem with 3 hubs, discount factor of 0.2, and the objective weight of 0.3 is 367.24; the locations of the hubs are 4 (Chicago), 2 (Baltimore) and 7 (Dallas-FW); the backup facility for nodes 1, 5, 6, 9 is 2 (Baltimore) and for nodes 3, 8 and 10 is 4 (Chicago) (Fig. 13).
For problem instances with 15 nodes and 3 and 5 hubs, CPLEX provides lower and upper bounds for 15 instances and finds only lower bounds for the remaining 3. A comparison of the average costs and the average computational times in Table 4 indicate that both H-PSO and PSO-v2 perform better than PSO-v1. Similar to the results for test problems with 10 nodes, no significant difference between the results provided by H-PSO and PSO-v2 is observed.
In test problems with 15 nodes, our proposed PSO based algorithms improved the upper bounds for the 18 instances with H-PSO algorithm achieving improvement up to 20%.
CPLEX did not provide any solution for larger instances with 20 and 25 nodes due to the lack of sufficient memory. A summary of the computational results for these 36 large problems are also presented in Table 4. Here, the three PSO algorithms could be clearly ranked with H-PSO on top in terms of both solutions quality and the computational times followed by PSO-v2 and PSO-v1.
Further observations and analysis Incorporating the probability of hub failure in decision making process of the hub location and allocation problem is a proactive strategy aimed at reducing the excessive financial and related non-financial costs when one or more facilities become unavailable. However, managers are less likely to use the resulting network with reliability consideration if a significant increase in the regular transportation cost incurs. To investigate this we consider a problem with 10 nodes and 3 hubs and empirically analyse the effect of hub failure consideration on the regular network transportation cost. The benchmark problem is taken from the CAB dataset and generated by setting the objective function weight w and the discount parameter of α to 0.3 and 0.8 respectively. The optimal network configuration for this problem without hub failure consideration is depicted in Fig. 14a with its corresponding transportation cost being 716.98 units. The resulting network when the possibility of hub failure (between 10 and 30%) is taken into consideration is depicted in Fig. 14b. The network presented in Fig. 14b is based on our approach and realistically reflects the real world. The (unweighted) regular transportation cost of the network in Fig. 14b is 752.65 unit which is 4.7% higher than that of the network provided by the classical model. This finding suggests that using our proposed model that takes to account facility unavailability may not result in a significant increase in regular transportation cost while providing a more robust configuration in case of hub failure. Furthermore, as illustrated in Fig. 14   suggests that solving a classical hub location problem to optimality and then deciding on the backup facilities is a simple but an ill-informed strategy that will lead to an inferior solution. This is an important strategic decision that ought not to be taken lightly.

TR dataset
The proposed formulation and algorithms have been also tested on 54 benchmark problems derived from TR dataset with non-symmetrical flow ranging from 10 to 55 nodes. The average costs and computational times over 18 problem instances with 10 nodes and 3 and 5 hubs are presented in the first row of Tables 6. Similar to the CAB dataset, CPLEX solved all benchmark problems with 10 nodes in reasonable computational times. It appears that, in general, CPLEX has less difficulty in solving instances from TR dataset as the average computational time required (280 sec) is about 25% lower than that used for CAB dataset (373 sec). This could be due to the dataset's structure and properties. Further examination of the results in Table 6 reveals that all three proposed PSO algorithms performed well in solving small to medium size instances and for the larger problems H-PSO outperforms the other two. This is the same performance pattern as the one observed earlier from the CAB dataset.

Scenario II: the case of low hub failure probability
In some applications of the hub-and-spoke systems it is common to assume a low probability of hub failure. To examine how the proposed formulation behaves under these circumstances and to assess what impact this change might have on the final solutions, we tested the same problem instances from the CAB dataset but with low probability of hub failures instead. As  Fig. 15 a Optimal network to classical single allocation hub location problem, b optimal network to proposed reliable single allocation hub location problem with low probabilities of hub failure mentioned earlier, in this study we use the low probabilities of hub failures as given in the paper by An et al. (2015).
To be consistent with the analysis conducted in the first scenario, all other parameters used in the formulation and the PSO algorithm are the same.
Furthermore, as H-PSO proved to be more powerful than the other two proposed algorithms from now on we only discuss the results provided by this algorithm. We begin with the small instances with 10 nodes and 3 and 5 hubs from the CAB dataset.
For convenience, the average cost of the optimal/best solutions for all problem instances tested from CAB dataset with low and high probabilities of hub failures are also presented in Table 7. From the computational point of view, it is clear that instances with low probabilities of hub failures could be solved with the proposed formulation (i.e., RSApHLP-III) much faster than that with high probabilities. For instance, the average computational times to solve all 18 benchmark problems with 10 nodes and 3 and 5 hubs are 127 and 373 s respectively.
H-PSO also solved two instances to optimality and improved the upper bounds for a number of those problem instances. The quality of the solutions provided by the H-PSO are clearly comparable with those provided by CPLEX in extended computational times. This observation demonstrates the efficiency and stability of our proposed hybrid PSO algorithm.
A total of 36 larger instances with 20 and 25 nodes are also solved using our hybrid PSO (H-PSO) algorithm. A summary of the results and detail description of each solution for these classes of benchmark problems are presented in Tables 7, 8 respectively. In Table 8 we report the total network cost, the hub locations, demand allocation and the selected backup facilities. For instance, for the 20 nodes instance with 3 hubs, discount factor 0.2, and the objective function weight of 0.3, the cost of the best solution found is 306.76 units; the selected hub locations are 20 (Pittsburgh), 7 (Dallas-FW), and 8 (Denver); the backup facility for nodes 1-4, 6, 12-14, and 17-19 is 7 (Dallas-FW); for nodes 5,9-10, 15 and 16 is hub 8 (Denver); and for node 11 it is hub 20 (Pittsburgh).

Conclusion
In this study, a mixed integer quadratic programming (MIQP) formulation is proposed for the reliable single allocation hub location problem with heterogeneous probability of hub failure. The model assigned a backup facility to every demand point in the network for a fast and a low-cost recovery after hub failure. The objective function of the model minimises the weighted transportation cost in regular situation and the expected cost resulted from hub failures.
The proposed formulation is developed in such way to maintain the size of the problem after linearization. To improve the computational efficiency of the resulting MILP model, all constraints initially formulated using Big Ms are substituted with logical constraints known as Indicator Constraints. Using the proposed formulation, we solved all small instances and some medium size problem instances to optimality and provided lower and/or upper bounds for other problems. To tackle large instances, three novel algorithms based on PSO namely PSO-v1, PSO-v2 and H-PSO are proposed. In PSO-v1 and PSO-v2, we investigated the effect of two different solution representations and the mechanisms through which a continuous particle is transferred to a discrete solution. It was found that PSO-v2 that maps a continuous particle to a single discrete solution performs better than PSO-v1. The performance of the PSO-v2 is then further improved by introducing two new interesting features namely a crossover operator and a memory to keep track of good solutions that are explored during the search. To our knowledge, this is the first time the hybridisation of PSO with those attributes is examined.
We presented and discussed the computational results for two scenarios. In the first scenario it is assumed that hub disruption is more likely to occur and therefore the selected probabilities of hub failures are set relatively high. Whereas in the second a hub disruption is assumed to be an event with a low probability of occurrence. Our numerical results show that even with low levels of probabilities of hub failures the regular transportation cost of the reliable network will not increase significantly and the amount of cost saving remains considerable. This finding, however, is based on experiments involving relatively small instances. Further research required to examine optimal solutions to large problem instances.
In this study the single allocation version of the hub location problem is considered, however the proposed model and our metaheuristic could be easily adopted to address the multiple allocation cases. Other venues of future research including related problems that incorporate capacity constraints and fixed costs are also worth examining.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.