Hub location under hub congestion and demand uncertainty: the Brazilian case study

Hub location under hub congestion and demand uncertainty: the Brazilian case study

Gilberto de Miranda Junior^I,*; Ricardo Saraiva de Camargo^II; Luiz Ricardo Pinto^III; Samuel Vieira Conceição^IV; Ricardo Poley Martins Ferreira^V

^IDepartamento de Engenharia de Produção, Universidade Federal de Minas Gerais, Belo Horizonte, MG, Brasil. E-mail: [e-mail]

^IIDepartamento de Engenharia de Produção, Universidade Federal de Minas Gerais, Belo Horizonte, MG, Brasil. E-mail: [e-mail]

^IIIDepartamento de Engenharia de Produção, Universidade Federal de Minas Gerais, Belo Horizonte, MG, Brasil. E-mail: [e-mail]

^IVDepartamento de Engenharia de Produção, Universidade Federal de Minas Gerais, Belo Horizonte, MG, Brasil. E-mail: [e-mail]

^VDepartamento de Engenharia Mecânica, Universidade Federal de Minas Gerais, Belo Horizonte, MG, Brasil. E-mail: [e-mail]

ABSTRACT

In this work, a mixed integer nonlinear programming model combining direct service links, demand uncertainty and congestion effects is proposed. This model is efficiently solved by Generalized Benders Decomposition, for instances of moderate sizes and reasonable number of scenarios. The deployed algorithms are further used for re-designing the Brazilian air transportation network, enabling the analysis of future demand scenarios and providing decision support about the optimal investment policy for Brazil.

Keywords: hub-and-spoke networks, Benders Decomposition, large scale optimization, mixed integernon-linear programming, demand uncertainty.

1 INTRODUCTION

In the last few years, the Brazilian air transportation network has gone through several operational breakdowns. The apparent reason is the lack of capacity or infra-structure on the main airports of the system. But a closer look of the network reveals that the total installed capacity is much larger than the total traffic volume to be handled, indicating that the network seems to operate in a sub-optimal way. As a consequence, the network performance has been severely compromised, reducing the opportunity of the involved logistic operators for attaining the economies of scale that lead to sustainable revenues.

In order to provide the necessary operational efficiency to this system, as well as to define optimal investment policies considering the system expansion and future demand sceneries, a mixed integer non-linear programming formulation is proposed. This model aims to redesign the current network, optimizing the use of the installed capacity and balancing the flows through the network. Here, the underlying structure adopted for the Brazilian air transportation network is a customized hub-and-spoke topology. The desired effect is eventually to deactivate some of the current network hubs while activating several others to handle the additional third party traffic.

Hub-and-spoke systems have been largely employed in the telecommunication and transportation areas (see [18, 43, 4, 42, 14]). This class of systems arises when commodities (passenger, cargo parcels, telecommunication packets) of a set of origins must be sent to a set of destinations. Instead of establishing a direct connection for each origin-destination pair, facilities that serve as switching, transshipment, sorting and distribution nodes, designated as hubs, are used as the only valid intermediate points in a path from an origin to a destination. Flows from different origins are gathered at these hub facilities prior to be routed to an intermediate hub or to be delivered to their final destinations.

The employment of these hub facilities and the routing of consolidated flows through inter-hub links allow the centralization of commodity handling and the transportation cost per unit of flow to be less expensive than directly shipping via a non-hub network structure. In other words, the hub networks take advantage of scale economies on inter-hub connections [58, 59].

Usually, the main concerns of hub network design are the location of hub facilities and the allocation of non-hub nodes to hubs at the least expensive network structure, accounting for that the installation and the transportation costs. This is done by selecting node candidates as hubs and assigning origin and destination nodes to these hubs, given the flow volumes between each origin-destination node pair, the installation costs of hubs and the transportation costs between nodes of the network. While hubs are in general fully connected among them, non-hub nodes can be single allocated, meaning that a non-hub can be assigned to exactly one hub only; and multiple allocated, meaning that a non-hub can be connected to more than one hub of the network. Alternative configurations can be considered in order to ensure the application of a specific protocol [62]. For further basic notions, different considerations and formulations of hub networks see the comprehensive surveys of Campbell [12, 13] and of Campbell et al. [14] and of Alumur & Kara [2].

The minimization of the installation and transportation costs approach generates solutions that tend to overload a small number of hubs, forcing the congestion effects into consideration. The hub literature has been dealing with the congestion cost effects on hub networks implicitly, when these effects are represented by constraints, and explicitly, when the congestion costs are expressed on the objective function.

Grove & O'Kelly [36] are among the first authors to study the effects of congestion on hub networks. They have demonstrated how the schedule delays of airline systems are influencedby the amount of flow at the hubs by simulating a single assignment hub network with fixed hub locations.

A great variety of researchers has tackled the congestion effects restricting the amount of flow transiting through a hub by means of capacity constraints. Aykin [3] has devised a lagrangian relaxation approach, while Ernst & Krishnamoorthy [27] have used simulated annealing and random descent methods. Ebery et al. [23] have imposed capacities only on incoming flows. They have developed a branch-and-bound based on a shortest path heuristic. Solutions obtained by their method may have flows from a hub to itself routed via another hub as explained by Campbell et al. [14]. Campbell et al. [17, 15, 16] have proposed the hub arc location problem, in which the hub network is seen as a three-layered network. The bottom one has the origin nodes of the demand, the middle layer has the hub nodes, and the top one has the destination nodes of the demand.

Yaman & Carello [72] have addressed the problem of solving single allocation hub location problems with modular link capacities. They have proposed a non-linear formulation where the size of inter-hub connections are treated as stepwise, instead of the linear approach proposed by Labbé et al. [44] to the quadratic capacitated hub location problem with single assignment. In order to accomplish that, capacity restrictions have been imposed on the amount of traffic flowing through the hubs rather than only considering the incoming traffic. The problem has been solved by means of an branch-and-cut algorithm based on a linearization with an exponential number of constraints and a two-level local search method. Both methods have been then combined in a heuristic concentration (see [68, 67]) scheme. Marín [54] also has addressed the capacitated multiple allocation hub location problems, but assumed that the flow between a given origin-destination pair can be split into several routes.

Another interesting work is the one presented by Kara & Tansel [41]. They have modeled the transient times at hubs in addition to the travel times, having as objective function the minimization of the latest-arrival times, and consequently, the excessive delays at the hubs. Yaman et al. [73] have presented a similar formulation for a Turkish cargo delivery system where transient times have been taken into account as well as journey times.

A work that stands out is the one of Marianov & Serra [53]. They have modeled the hub network as an M/D/c queuing network, proposing capacity constraints based on the probability of waiting customers in the system. Due to the computational complexity of these constraints, they have linearized them and then solved the resulting model by a tabu search algorithm. A similar work has been proposed by Rodriguez et al. [65], but a simulated annealing algorithm has been used.

Costa et al. [22] have presented a multi-criteria formulation to the capacitated single allocation hub location problem. Besides the traditional cost minimizing function, they have considered alternatively the minimization of the time the hubs take to process the flow or the minimization of the maximum service time of the hubs. For instance, Rodriguez-Martin and Salazar-Gonzalez [66] have proposed a Benders decomposition algorithm [6] to tackle a capacitated hub location problem based on a multi-commodity formulation where the arcs connecting the hubs are not assumed to create a complete graph. Their model is similar to the one presented by Sridhar and Park [71] for the fixed-charge capacitated network design problem.

Modeling congestion through capacity constraint on the flows does not mimic the exponential nature of the congestion effects. Elhedhli & Hu [24] have been the first ones to consider the costs of the congestion effects explicitly on the objective function. Using a convex cost function that increases exponentially as more flows go through the hubs, they have proposed a non-linear model to a single assignment hub location problem. The congestion convex cost function is a power-law function of airport usage relative to its capacity and it has been widely used to estimate delay costs [35]. They have linearized their model and then solved it using lagrangian relaxation.

On the other hand, stochastic location problems where demand uncertainties are taken into account have been addressed by G. Albareda-Sambola et al. [1], Hochreiter & Pflug [40], Snyder et al. [70], Gabor & Van Ommeren [31], Snyder & Daskin [69], Laporte et al. [46], Louveaux & Peeters [48], Laporte et al. [45] and França & Luna [28]. In all these works, it is established that stochastic modeling provides a superior design by considering uncertainty in important parameters. When dealing with HS network design, it is interesting to consider the demand uncertainty, since, according to Camargo et al. [10], the impact of demand components in the optimal HS network under congestion is considerable. Another interesting conclusion of Camargo et al. [10] is that under huge congestion, the marginal utility of installing a new hub in the HS network is reduced. This counter-intuitive effect is due to the HS classical protocol - if a given demand is considerably high, it is not interesting to have this link on the HS network, since it will always imply in high congestion. It is clear that if a given node is overloaded with his own demand, it is unable to handle third party traffic. It is also uninteresting to route such high flows trough a hub, under penalty of severely impact the hub performance. A further step could be the allowance of direct routes between non hub nodes, where no economy of scale is observed but no congestion either.

In this paper, we explore further the congestion effects written as a convex cost function similarly to [24], but addressing the multiple allocation hub location problem, considering also effects of demand uncertainty combined with direct service links. We propose a non-linear mixed integer programming problem based on a classical formulation [37] due to its linear programming bound quality when compared to others [12, 25, 26, 61]. In our formulation the number of hubs on a route is limited to two, even if there is a route with a lower cost using more than two hubs.

Although, one of the main strategies to handle non-linear problems is to linearize them, we have developed a Generalized Benders Decomposition (GBD) [33] to cope with our non-linear mixed integer program. The problem is decomposed in two smaller problems: at a higher level, named as master problem (MP), the location decisions are made; while at an inferior level, known as subproblem (SP), the flow balance and the congestion are handled. The MP is an integer programming problem, while the SP is a non linear convex transportation problem. Using our approach, we have been able to solve some large problems of the CAB and AP standard data sets to optimality. We have also used the deployed algorithm to deal with a real hub location case, studying the Brazilian air transportation network.

This paper is organized as follows: in Section , the model formulation and the used notation are introduced. In Section , the Generalized Benders Decomposition is developed to solve the problem and it is also demonstrated how the subproblems can be optimally solved. Computational experiments and conclusion remarks are presented in Section and , respectively.

2 MODEL FORMULATION

In this section the our model formulation efforts are presented. The basic components of the model are the following sets: let N be the set of node locations that exchange traffic and let H be the set of node candidates to become hubs, H ⊆ N. For any pair of nodes i and j (i, j ∈ N), we have w_ij, the flow from origin i to destination j that is routed through either one or two installed hubs. Usually, we have w_ij ≠ w_ji.

Further, let f_k be the fixed cost of installing a hub at node k ∈ H and let c_ijkm be the transportation cost per unit of flow from node i to j routed via hubs k and m (i, j ∈ N and k, m ∈ H).This transportation cost is the composition of three cost segments: c_ijkm = c_ik + αc_km + c_mj, where c_ik and c_mj are the standard transportation cost per unit from location i (j) to hub k (m), and αc_km is the discounted transportation cost between hubs k and m. The discount factor 0 < α < 1 represents the scale economies on the inter-hub connection. If only one hub is used in a given route, we have k = m and no discount factor is applied.

The following decision variables are defined:

y_k =

x_ijkm> 0 is the flow from origin i to destination j (i, j ∈ N) that is routed through hubs k and m (k, m ∈ H) in that order.

In order to start the development of the congested version of the uncapacitated multiple allocation hub location problem (UMAHLP), the formulation due to Hamacher et al. [37] is stated as:

s. t.:

The model (1)-(5) is reputed to have a tight linear programming bound [37], being one of the strongest formulations for the UMAHLP available. One of its major drawbacks is its size. However, as demonstrated by the pioneering work of Geoffrion & Graves [32], structure is more important than size when dealing with large scale mathematical programs.

In Camargo et al. [11] the above formulation has been used to solve instances up to 200 nodes with the aid of Benders decomposition. In this sense, formulation - is adopted as a basis to develop a congested approach to hub location due to its nice theoretical properties and good structure. Hamacher et al. [37] have demonstrated that constraints (3) are facet defining. In addition, for a fixed hub structure, the remaining formulation is fully decomposable for each i - j pair. In spite of this, O'Kelly & Brian [60] have shown that having the same discount factor α for all the inter-hub connections may be unrealistic. Optimal solutions where large discounts are given on links having small flows can be generated in this way. Overcoming this drawback, a hub location formulation with flow-dependant economies of scale can be addressed, as proposed by Camargo et al. [9].

In order to allow the use of direct service links, is necessary to add a new set of variables:

θ_ij> 0 is the amount of flow transported through a direct connection for pair i - j.

It is also necessary to modify constraints (2) to include the direct paths, producing:

Implying to rewrite the objective function as:

Moreover, Marín et al. [55] have shown how formulation (1)-(5) can allow optimal solutions with routes having more than two hubs when the triangular inequalities are not satisfied. Instead of using the constraints derived by Marín to deal with those possible unfeasible paths, the following additional constraints are preferable in a decomposition framework, since they preserve the i - j pair decomposition:

Constraints (8) make sure that the flow w_ij goes through route i - j - i - j when i and j are hubs. Constraints (9) and (10) limit the number of hubs to two when i or j are hubs. Furthermore, the already deployed subproblem solution algorithm of Camargo et al. [11] can be partially reused. Remark that the above cited origin-destination special structure is a consequence of the interchangeability of indexes i, j, k and m, i.e. H ≡ N. It is also straight-forward to modify the model and the deployed algorithms if H ⊆ N is required. Since we are considering the congestion effects as a function of the total hub flow, we need an extra decision variable for each hub to account it:

g_k> 0 is the total flow which goes through hub k ∈ H.

Further, a set of additional coupling constraints have to be added in order to amount the total flow in a given hub:

A congestion cost function τ_k(g_k) is assumed as given for each installed hub k ∈ H. This function is considered increasing on [0, +∞), proper convex and smooth. It is responsible for distributing the flow load along each installed hub, thus reducing the overall congestion expressed in monetary values. The congestion cost function used is a power law where τ_k= e, and e and b are positive constants, b > 1. This class of functions has been already used by [24] to the single allocation case and by [10] in his recent work. The congested multiple allocation hub location problem with direct service can now be stated as follows:

s. t.:

In the above model the objective function (12) minimizes the total cost, e.g., the sum of theinstallation, congestion and linear transportation costs. Constraints (13) determine the total amount of flow which goes through hub k. Constraints (14) assure that the flow for every pair i - j is routed via some hub pair or use a direct service link. Once again, if only one hub is used, we have k = m. Constraints (15) guarantee that the flow going through a hub only happens if that hub is installed. Constraints (19), (20) and (21) avoid problematic routes as discussed above. Constraints (16), (17) and (22) are the non-negativity constraints of the continuous variables x_ijkm, r_ij and g_k, while the constraints restrict the integer variables y_k to be 0 or 1.

In order to consider effects of demand uncertainty it is possible to adopt a discrete set of states of the world S, associating each one to a probability p^s, s ∈ S, rewriting the demand parameter as , the flow from origin i to destination j in the state of the world s ∈ S. Clearly, the formulation (12)-(22) must be generalized to deal with the new stochastic parameters. We have chosen to replicate the sets of decision variables x_ijkm, θ_ij, and g_k and associated constraints for each state of the world, in order to minimize the sum of fixed (installation) costs plus the expected variable (linear transportation and congestion) cost components. The first step is to redefine the above cited decision variables:

> 0 is the flow from origin i to destination j (i, j ∈ N) that is routed through hubs k and m (k, m ∈ H) in that order, in the state of the world s ∈ S.

> 0 is the amount of flow transported through a direct connection for pair i - j in the state of the world s ∈ S.

> 0 is the total flow which goes through hub k ∈ H in the state of the world s ∈ S.

Rewriting formulation (12)-(22) by introducing the new sets of decision variables leads to the following stochastic programming problem, that we are calling hub location problem under hub congestion and demand uncertainty:

s.t.:

As the above model has not capacity constraints, there is always a feasible solution for (23)-(33). Furthermore, for any given fixed structure of hubs, i.e. a feasible vector y, the model transforms into a non-linear convex transportation model. This non-linear transportation problem can then be solved by a feasible directions search method such as [29] or by the proximal decomposition algorithm of Mahey et al. [52]. Gathering all of these features, we have been motivated in deploying the Generalized Benders decomposition method (GBD) [6, 33] to tackle the problem.

Concerning GBD, the algorithm performance is very sensitive to the overall contribution of the non-linear terms, since it is a first order approach. This is particularly true if the associated mixed integer linear program has a poor linear programming bound (see [51, 28, 63]). For this reason, formulations having less variables but weaker lower bounds, as those proposed by Ernst & Krishnamoorthy [26], may be disregarded for this kind of application.

3 APPLYING THE GENERALIZED BENDERS DECOMPOSITION (GBD) METHOD

In 1962, Benders [6] proposed a partitioning method for solving mixed linear and nonlinear integer programming problems. Benders has defined a relaxation algorithm which partitions the original problem into two simpler problems: the relaxed master problem (RMP) and the subproblem (SP). The first one is a relaxed version of the original problem with the set of integer variables and its associated constraints; while the last one is the original problem with the values of the integer variables fixed by the RMP.

In general terms, Benders decomposition relies on a projection problem manipulation, followed by solution strategies of dualization, outer linearization and relaxation. The complicating variables (integer variables) of the original problem are projected out, resulting into an equivalent model with fewer variables, but many more constraints. However, when attaining optimality, a large number of these constraints will not be binding. Hence a strategy of relaxation can be devised to ignore all but a few of these constraints.

Benders [6] has suggested a procedure to add these constraints on demand, adopting two levels of coordination. At a higher level, the RMP is responsible for fixing the values of the integer variables and for providing a lower bound for the problem. At a lower level, the subproblem is solved in order to obtain an upper bound and a violated cut on the master problem space. This valid inequality is known as Benders cut and will enrich the RMP in the next iteration. The procedure stops when the upper bound and lower bound converge towards the optimal value.

The Benders decomposition has been successfully applied on solving large scale problems since the pioneering paper of Geoffrion and Graves [32], the uncapacitated network design problem with undirected arcs of Magnanti et al. [49], the locomotive and car assignment problem of Cordeau et al. [20, 21], the cellular manufacturing system design of Heragu & Chen [38], the multicasting network design problems of Miranda [57] and Randazzo & Luna [64] and the hub location problem with economies of scale of Camargo et al. [9]. The original procedure has been specialized to deal with stochastic programming problems, receiving sometimes the name of L-Shaped Decomposition Method on the literature of stochastic programming, due to the form of the constraint matrices generated, amenable to dual decomposition techniques.

Several authors have addressed other extensions and improvement of the Benders decomposition method. Magnanti & Wong [50] have proposed an approach for selecting the dual variables and thus improving the overall performance of the algorithm. McDaniel & Devine [56] have developed an algorithm for accelerating the decomposition by removing temporarily the integrality constraints of the RMP. Balas & Bergthaller [5] have devised a nice approach for generating strong cuts prior of the beginning of the Benders iterations. Geoffrion [33] has extended the basic algorithm to non-linear problems, formulating the Generalized Benders Decomposition (GBD). The GBD has been employed on different problems. Cai et al. [8] have used it on a nonconvex water resource management problem, Hoang [39] has applied it on a nonlinear model for network design, while Mahey et al. [51] have solved capacity and flow assignments in telecommunication networks by means of GBD.

When GBD is used to tackle the congestioned hub design problem, the RMP chooses a hub structure configuration represented by the decision variables y, and the SP solves the non-linear transportation problem with price sensitive demands on the resulting network. The dual solution of the transportation problem is then used to generate a Benders cut that will redefine the hub structure configuration in the next iteration. Adjustment of the non-linear objective function is necessary at the subproblem. This is accomplished by dualizing the coupling constraints associated to the non-linearity. After dualization, the resulting subproblems are relatively easy to solve.

In the next Sections 3.1 and 3.2, we present the relaxed master problem and the subproblem, respectively, of the congestioned hub-and-spoke network design problem.

3.1 Relaxed Master Problem (RMP)

From the viewpoint of mathematical programming, a projection of the problem (24)-(33) onto the space of the structural variables y can be done, resulting on the following problem to be solved:

where Y = {y ∈ {0, 1}| for y fixed, there are feasible flows satisfying (24)-(28) and (30)-(33) and z(y)^s is the following subproblem to be solved for each state of the world s:

subject to (26), (30)-(32) for y fixed and the appropriate s, and where

G = {(x, g, θ)| x > 0 and g > 0 and θ > 0, satisfying (24), (25), (27), (28), (33)}.

For any hub structure configuration there is a feasible path connecting the pair of demands i - j, hence, there is no need for further feasibility constraints on the domain of the projected problem (34). Further, since the subproblem has a convex differentiable objective function and linear constraints, its Karush-Kuhn-Tucker conditions are necessary and sufficient for optimality. This makes the problem consequently amenable to dualization techniques (see [33, 47]).

After associating the vectors of dual variables u > 0, v > 0, t > 0 and r > 0 to the coupling constraints (260, (30)-(32) and since there is no duality gap, the optimal value of the subproblem (35) for any y ∈ Y can be given by:

Recalling the interchangeability of indexes i, j, k and m, z(y)^s can be rewritten in a more compact form as:

Defining as:

for all k ∈ H, can now be stated as:

Hence, the whole problem (34) is equivalent to:

As the supremum is the least upper bound, the mathematical program (23)-(32) can be represented as the following master problem (MP) with the help of variable η > 0:

s. t.:

Problem (38)-(40) is solved by the GBD through a strategy of relaxation, i.e. ignoring all but a few of the constraints (39). This strategy requires a procedure which adds iteratively the constraints to the MP as needed. So, at any iteration h, the optimal value z(y^h)^s is found from (37) after solving the non-linear transportation subproblem for a given installed hub structure y^h and recovering the vector of optimal multipliers (u^s)^h, (v^s)^h, (r^s)^h and (t^s)^h (u^s = (u^s)^h, v^s = (v^s)^h, r^s = (r^s)^h and t^s = (t^s)^h). Thus the optimal value of z(y^h)^s is given by:

From (39), we have the following constraint associated to u^h, v^h, r^h and t^h for the state of the world s:

Isolating the minimum given by (41), we get the following Benders cut based on the y^h and (u^s)^h, (v^s)^h, (r^s)^h and (t^s)^h for the state of the world s:

Resulting then on the following relaxed master problem:

s. t:

where L is the maximum number of Benders cuts to be taken in order to attain the optimal solution. Note that, at each iteration h, new Benders cuts of the form (43) can be generated provided that the optimal values of z(y^h)^s and the vector of multipliers (u^s)^h, v^s)^h, (r^s)^h and (t^s)^h have been recovered. These values are attained by solving the subproblems, the subject of the next Section 3.2.

3.2 Subproblems

For a hub structure y^h fixed by the RMP at iteration h, the subproblem z(y^h)^s, see equation (35), to be solved at the state of the world s is given by:

s. t.:

Since the function τ_k () is increasing, the coupling constraints (46) can be dualized by associating the Lagrangian multipliers β^s > 0. Thus the resulting problem is separable into two subproblems: a linear transportation problem, d_L(β^s), having only the and the variables; and another, a convex flow assignment transportation problem, d_NL(β^s), having only the variables;

subject to constraints (47)-(52).

With the correspondent dual variables b^s ∈ , the optimal value z(y^h)^s can then be computed as:

Recalling that d_NL(β^s) has a convex differentiable objective function and linear constraints, its Karush-Kuhn-Tucker conditions are necessary and sufficient for optimality. So, for each k ∈ H and a fixed hub structure , the optimal solution of ()^h should minimize its correspondent parcel of the d_NL(β) objective function, what implies, for iteration h:

The optimal values of ()^h can be obtained by employing an algorithm based on the Flow Deviation algorithm (FD) [34], described in the Section 3.2.1.

3.2.1 The Non-linear Transportation Problem

The subproblem z(y^h)^s, given by (45)-(53), is a non-linear transportation problem that can be solved by means of a feasible directions search method such as [29] or the proximal decomposition algorithm of Mahey et al. [52]. In the present paper, a specialized Flow Deviation [34] algorithm has been deployed.

The Flow Deviation algorithm is an optimization method designed to find the minimum of a convex function subject to linear constraints. The algorithm has been shown to be efficient since the early works of Fratta et al. [30] and Cantor & Gerla [19]. In order to apply the method to the subproblem (45)-(53) and thus compute the optimal values of , a more adequate formulation has to be applied. An arc-node formulation with side constraints has been used. The role of the side constraints is to restrain the set of feasible paths from an origin to a destination.

The specialized Flow Deviation algorithm takes advantage of the fact that each commodity may only follow a reduced set of routes. Because of this, when evaluating the shortest path for a given i - j pair, it is not necessary to deploy an algorithm like Dijkstra, but only enumerate this reduced set. Furthermore, as it is not necessary to assess the flow of each commodity in each arc during the running of the algorithm, but just the total flow in each arc, specially for each hub, it is possible to save computer memory and processing time.

Once computed the values of (g^s)^h and hence the values of (β^s)^h (see (54)), the subproblem d_L(β^s) can now be solved, subject of the next Section 3.2.2.

3.2.2 The Linear Transportation Problem

After fixing the optimal values of (β^s)^h and (g^s)^h, the subproblem d_L(β^s) can now be solvedand thus the values of (u^s)^h, (v^s)^h, (r^s)^h and (t^s)^h are now available. Instead of solving the d_L(β^s) by means of simplex or interior point methods, it can be solved by specializing an algorithm which relies on the complementary slackness condition and on the linear programming duality theory.

As proved in Camargo et al. [11] the primal problem d_L(β^s) is always feasible and bounded for any fixed feasible y_k satisfying constraints (47)-(52) as well as its dual. The linear programming dual is also decomposable for each i - j pair and state of the world s:

s. t.:

Proposition 3.1 The optimal solution of the dual problem (55)-(67) can be found by inspection.

Proof. The solution of (55)-(68) is started by setting as the shortest path for a given i - j pair considering the installed structure y^h at the state of the world s (see Proposition 2 of [11]).

Starting values for are already available, when k ≠ i and k ≠ j, from the algorithm proposed in [11] by adding the respective values of β^s to the linear transportation costs. In order to find the optimal values of the other dual variables, please refer to the four steps procedure described in Proposition 1 of [10], repeating it to each state of the world s.

□

4 COMPUTATIONAL EXPERIMENTS

This section is composed by three distinct subsections, each one aiming to explore a different feature of the proposed model and deployed algorithms. In the first subsection, the superior quality of the stochastic design is illustrated by comparing the outcomes of deterministic and stochastic modeling for a toy instance. The second subsection investigates the scalability of the deployed algorithms, testing them for benchmark data sets. In the third section, the case study is described in detail, and the obtained results are properly discussed.

4.1 Deterministic versus stochastic hub-and-spoke network design

The goal here is to show the superior quality of stochastic hub location if compared to deterministic hub network design. The adopted parameters for the hub congestion functions were e = 0.0001 and b = 2. These experiences have been made using XPRESS-MP QMIP solver. As working problem instance, an adapted version of CAB10 has been produced. The generated demand matrix is showed below, having only three distinct sceneries:

The adopted probability vector was p = [0.333 0.333 0.333]^T. In order to establish a fair comparison, a deterministic equivalent demand matrix has been produced by taking w_ij = Σ_{s ∈ S}p^s. This deterministic equivalent instance was solved by formulation (12)-(22). The resulting optimal hub structures (y vectors) are depicted ahead:

Deterministic: y* = [0 1 0 0 1 0 0 1 0 1]^T

Stochastic: y* = [0 0 0 0 1 0 0 1 0 1]^T

In spite of the detailed report of a single experience - for the sake of simplicity - several other tests have been carried out. They were done by tuning the probability vector and the hub installation costs. In general, the stochastic hub network design is a bit more conservative, saving fixed costs, when compared to the deterministic solution. It is also true that the same hub structure has been observed many times, meaning that a deterministic solution is almost always optimal or near-optimal. However, considering the costs involved in such networks and the long term impact on its performance, the risk of designing a non-optimal network by disregarding the uncertainty is definitely not affordable.

4.2 The GBD Algorithm Scalability

In order to test the performance of the deployed generalized Benders decomposition algorithm, the CAB data set of the US Civil Aeronautics Board and the AP data set introduced by Ernst and Krishnamoorthy [25, 27] have been used. The former includes problems with 10, 15, 20 and 25 nodes, while from the later problems with 10 to 200 nodes can be obtained.

The hub installation costs have been generated for all instances using a gaussian distribution with average equals to f_o and variance set to 40% to assess how the different fixed costs vary in real problems. Where f_o represents the scaled difference in objective value between a scenario in which there is a virtual hub located in the center of mass and a scenario in which all nodes are hubs, as introduced by Ebery et al. [23]. For alternative works considering these test beds see [60, 24, 23, 7, 11, 10]. Furthermore, as done by Ebery et al. [23], the nodes with higher flows have been selected to have the higher fixed costs, hardening, in general, the selection of potential hubs.

The generalized Benders decomposition algorithm was implemented with the aid of the general purpose mixed integer programming package XPRESS-MP from Dash optimization. The package routines have been used on the solution of the mixed integer master program. The experiments have been carried out on a regular desktop micro-computer equipped with a Intel-Pentium IV 3.0 GHz processor and 2Gb of RAM memory. The value of ε has been set to 0.0001%.

The attained results are displayed in Table 1. In this table, the columns Variables and Iterations account for the total number of decision variables of the solved instance (integer, nonlinear and continuous) and the total number of global iterations of GBD algorithm. The column SP time fraction (%) accounts the ratio of time spent to solve all the subproblems and the total solution time. The optimization of hub location problems under uncertainty is expensive, specially concerning the subproblem solution where a large scale non-linear program may be solved for each scenery of each iteration. In spite of this fact, some large problems could be solved in reasonable time, for a long term strategic decision problem.

The large number of iterations is explained by the aggressive congestion cost functions that were adopted. For all these tests the parameters e and b were taken as 0.5 and 2 respectively, implying a congestion cost ratio [10] scaling from 20% to 40%. Since GBD is an outer linearization strategy, it is particularly sensitive as the nonlinearities become dominant in a given MINLP.

4.3 Designing the Brazilian Air Transportation Network

Brazil is the fifth largest country of the world, and its area is 8.5 millions of squared kilometers, occupying about 47% of the area of South American continent. The population of the country is 184 million of inhabitants. Due to its great territorial extension and to the size its population, the country needs to have a great number of airports to allow the transportation of passengers and cargo. Nowadays, the country has 1759 private airports and 739 public airports, totaling 2498 airports.

Brazil is the second country of the world that more it possesses airports, behind just of the United States of America. The boards of INFRAERO and ANAC, two government agencies, are in charge of the management of 67 of those airports, that concentrate 97% of the commercial air transportation of the country. Those airports move about 2 million landings and takeoffs, and almost 50 million passengers and 1.2 million of tons per year (basis: year 2007).

In the last few years, Brazilian GDP is growing more than the average global GDP increase, around 5% year. This means an expected high demand increase for air transportation, since each 1% of increase on GDP means 1.15% of demand increase in the air transportation matrix. The passenger traffic in Brazil is strongly concentrated in few airports located in the cities of São Paulo, Rio de Janeiro and Brasília. The existent airports in those cities are the current passenger hubs of the Brazilian network.

Since the governmental investment profile has not matched the increasing traffic of the last few years, these hubs are severely overloaded, and the whole network is experiencing huge delays and several flight cancelations. On the other hand, several other airports appear to have idle capacity, such as the airports in the cities of Belo Horizonte, Salvador and Curitiba.

In order to use the deployed algorithms, several modeling assumptions have been made. Because of the large computational costs involved, only the 30 larger Brazilian airports were selected to compose the experience. Since INFRAERO and ANAC (Brazilian Civil Aeronautics Board) have historical series of the observed demands for each origin-destination pair in the network, as well as the distance matrix, the data mining phase was not a trouble. Another assumption was the use of the most economic available aircrafts, since the airlines always seek for the lowest operational costs.

The hub installation costs were estimated on the minimal infra-structure investment necessary to promote a node to the hub status. In some cases, this investment means only to improve air traffic control capabilities or to enlarge the passengers terminal or the cargo manipulation area. Sometimes this may imply to enlarge the runway length or to duplicate the existent runway or even to provide a better transportation system to reach the airport. These parameters were established according to the information provided by ANAC and INFRAERO.

In addition, a slightly different congestion function was used to account the congestion costs. For this set of experiences a function were the congestion effects are perceived only if a flow above 75% of the hub nominal capacity is reached was adopted:

τ_k(g_k) = max{0, e(g_k - a_k)^b}

This assumption is consistent with the most works on the literature about congestion cost functions. The capacity array for the Brazilian airports has been provided by ANAC and INFRAERO and is given in Table 2. All the capacities are expressed in passengers per year.

A standard discount factor α of 0.5 for the inter-hub connections was adopted, since almost all the airlines establish a minimum of 50% of tickets sold by flight leg to be sustainable. The values of the parameters e and b were taken as 1.0 and 2 producing very aggressive congestion effects (see [10]). To mimic the trend of growth of the traffic matrix, a triangular probability distribution was adopted, from 2011 to 2023. The central trend was taken forecasting a mean growth of 5.75% per year, corresponding to 5% of mean growth in GDP. The optimistic and pessimistic forecasts - the extremes of the triangular distribution - were based on the historical data from ANAC. Discrete distributions with ten sceneries were taken as approximations of the adopted triangular distributions. This small number of sceneries is necessary to ensure a reasonable computational cost, dealing with a little loss of accuracy.

The obtained optimal hub structures - y vectors - and the expected total traffic for each hub are shown in Table 3. All the proposed instances have run in computing times under 15 hours of CPU. We must remark that only the flows on the hubs were displayed in Table 3, meaning that if a node is not a hub his traffic is not shown in Table 3. The resulting network designs are plotted graphically in Figures 1, 2, 3 and 4. In these figures, the relative thickness of a given link is a measure of its traffic intensity.

It is important to recall that the capacity array informs just the nominal capacities, since the true capacity of an airport is not known with accuracy. For instance, in 2007 SBSP has registeredmore than 18 millions of passengers, but its nominal capacity is only 15 millions. So, it is not unexpected to observe flows slightly above that nominal value, if it is preferable to tolerate some congestion to avoid the installation of new hubs. Other effect that is remarkable is the transfer of part of the increasing demand to the direct service links, as a way to save fixed or congestion costs. This flow transfer explains why the total hub traffic does not grows in the exactexpected rate.

It is possible to see that the infrastructure investment may be directed to the airports of Belo Horizonte (SBCF), Recife (SBRF), Curitiba (SBCT) and Salvador (SBSV). These potential new hubs have idle capacity and may be better exploited in order to reduce the operational costs and flight delays on the Brazilian air transportation network. It is possible to see that it is also necessary to develop a strategy of expansion of all the network hubs, under penalty of a new traffic saturation in the near future. This expansion is particularly easy for SBCF and SBSV, complaining with data of INFRAERO and ANAC.

It was also observed an increasing transfer of traffic to the direct service links in order to avoid the hub congestion costs in 2019 and 2023, what recalls the urgency of a system expansion strategy. Another interesting effect is that the very huge traffic volume of the pair SBSP-SBRJ was not enough to make SBRJ a hub, since this particular airport has no capacity to deal withthird party traffic.

Another specially problematic airport is SBSP, having a preferential location and huge demands, but also having an insufficient infra-structure and being hard to expand, since it is inside the city of São Paulo. As the demand is increased from 2007 to 2023, no hub is removed but several hubs are added to the network. This result is quite interesting and could not be anticipated a priori. Additional research may investigate how these network designs are affected by allowing capacity expansions. However, capacity expansion costs are hard to estimate for several of the airports under study, except for SBCF and SBSV (where these costs are assumed small).

5 CONCLUDING REMARKS

In this article a model for accounting demand uncertainty and congestion costs for designing multiple allocation hub and spoke networks has been proposed. In addition, the proposed model allows the establishment of direct service links if they are economically interesting.

Generalized Benders decomposition was the selected technique to tackle the proposed problem. This technique has shown enough computing power to overcome large instances in reasonable time, although the computing times observed were expressively superior to those in [10], due to the most stressing congestion cost functions chosen.

The deployed algorithms were also used to redesign the Brazilian air transportation network, and the attained designs were properly discussed. The resulting network designs indicate that is necessary to redistribute flow on several hubs of the network in order to reduce the flight delays and other congestion effects. It has also pointed out that a system expansion is necessary in the near future under penalty of a new saturation of the installed infra-structure.

Future work must be devoted on evaluating how the capacity expansion of the major hubs may impact the proposed optimal hub-and-spoke networks.

ACKNOWLEDGEMENTS

We would like to thank the government of the state of Minas Gerais, Brazil, for partially supporting this work. We also thank the traffic engineers Erivelton Guedes and Aurélio Braga for the data and the insights they have provided from ANAC and INFRAERO. The authors' research has been partially funded by CNPq, project 480757/2008-9, and FAPEMIG, project PPM-00399-09.

Received February 1, 2010 / Accepted March 17, 2011

• [1] ALBAREDA-SAMBOLA M, FERNANDEZ E & LAPORTE G. 2007. Heuristic and lower bound for a stochastic location-routing problem. European Journal of Operational Research, 179(3): 940-955.
• [2] ALUMUR S & KARA BY. 2008. Network hub location problems: The state of the art. European Journal of Operational Research, 190: 1-21.
• [3] AYKIN T. 1994. Lagragian relaxation based approaches to capacitated hub-and-spoke network design problem. European Journal of Operational Research, 79: 501-523.
• [4] AYKIN T. 1995. Network policies for hub-and-spoke systems with applications to the air transportation system. Transportation Science, 29: 201-221.
• [5] BALAS E & BERGTHALLER C. 1983. Benders method revisited. Journal of Computational and Applied Mathematics, 9: 3-12.
• [6] BENDERS JF. 1962. Partitioning procedures for solving mixed-variables programming problems. Numerisch Mathematik, 4: 238-252.
• [7] BOLAND N, KRISHNAMOORTHY M, ERNST A & EBERY J. 2004. Preprocessing and cutting for multiple allocation hub location problem. European Journal of Operational Research, 155: 638-653.
• [8] CAI X, MCKINNEY DC, LASDON LS & WATKINS JR DW. 2001. Solving large nonconvex water resources management models using generalized Benders decompositon. Operations Research, 49(2): 235-245.
• [9] CAMARGO RS, MIRANDA G & LUNA HPL. 2009. Benders decomposition for hub location problems with economies of scale. Transportation Science, 43: 86-97.
• [10] CAMARGO RS DE, MIRANDA JR G, FERREIRA RPM & LUNA HP. 2009. Multiple allocationhub-and-spoke network design under hub congestion. Computers and Operations Research, 36(12): 3097-3106.
• [11] CAMARGO RS DE, MIRANDA JR G & LUNA HP. 2008. Benders decomposition for the uncapacitated multiple allocation hub location problem. Computers and Operations Research, 35(4): 1047-1064.
• [12] CAMPBELL JF. 1994. Integer programming formulations of discrete hub location problems. European Journal of Operational Research, 72: 387-405.
• [13] CAMPBELL JF. 1994. A survey of network hub location. Studies in Locational Analysis, 6: 31-49.
• [14] CAMPBELL JF, ERNST AT & KRISHNAMOORTHY M. 2002. Hub location problems. In: DREZNER Z & HAMACHER HW (Eds.). Facility Location: Application and Theory, pages 373-407. Springer.
• [15] CAMPBELL JF, ERNST AT & KRISHNAMOORTHY M. 2005. Hub arc location problems: Part I - Introduction on results. Management Science, 51(10): 1540-1555.
• [16] CAMPBELL JF, ERNST AT & KRISHNAMOORTHY M. 2005. Hub arc location problems: Part II - Formulations and optimal algorithms. Management Science, 51(10): 1556-1571.
• [17] CAMPBELL JF, STIEHR G, ERNST AT & KRISHNAMOORTHY M. 2003. Solving hub arc location problems on a cluster of workstations. Parallel Computing, 29: 555-574.
• [18] CAMPBELL JF. 1992. Location and allocation for distribution systems with transshipments and transportation economies of scale. Annals of Operations Research, 40: 77-99.
• [19] CANTOR DG & GERLA M. 1974. Optimal routing in a packet-swiched computer network. IEEE Transations on Computers, c-23(10): 1062-1069, October.
• [20] CORDEAU JF, SOUMIS F & DESROSIERS J. 2000. A Benders decomposition approach for the locomotive and car assignment problem. Transportation Science, 34: 133-149.
• [21] CORDEAU JF, SOUMIS F & DESROSIERS J. 2001. Simultaneous assignment of locomotives and cars to passenger trains. Operations Research, 49(4): 531-548.
• [22] COSTA M DA G, CAPTIVO ME & CLÍMACO J. 2008. Capacitated single allocation hub location problem - a bi-criteria approach. Computers and Operations Research, 35: 3671-3695.
• [23] EBERY J, KRISHNAMOORTHY M, ERNST A & BOLAND N. 2000. The capacitated multiple allocation hub location problema: Formulations and algorithms. European Journal of Operational Research, 120: 614-631.
• [24] ELHEDHLI S & HU FX. 2005. Hub-and-spoke network design with congestion. Computers & Operations Research, 32: 1615-1632.
• [25] ERNST AT & KRISHNAMOORTHY M. 1996. Efficient algorithms for the uncapacitated single allocation p-hub median problem. Location Science, 4: 139-154.
• [26] ERNST AT & KRISHNAMOORTHY M. 1998. Exact and heuristic algorithms for the uncapacitated multiple allocation p-hub median problem. European Journal of Operational Research, 104: 100-112.
• [27] ERNST AT & KRISHNAMOORTHY M. 1999. Solution algorithms for the capacitated single allocation hub location problem. Annals of Operations Research, 86: 141-159.
• [28] FRANÇA PM & LUNA HPL. 1982. Solving stochastic transportation-location problem by generalized benders decomposition. Transportation Science, 16(2): 113-126.
• [29] FRANK M & WOLFE P. 1956. An algorithm for quadratic programming. Naval Research Logistics Quarterly, 3: 95-110.
• [30] FRATTA M, GERLA M & KLEINROCK L. 1973. The flow deviation method: An approach to store-and-forward communication network design. Networks, 3: 97-133.
• [31] GABOR AF & VAN OMMEREN JCW. 2006. An approximation algorithm for a facility location problem with stochastic demands and inventories. Operations Research Letters, 34(3): 257-263.
• [32] GEOFFRION AM & GRAVES GW. 1974. Multicomodity distribution system design by Benders decomposition. Management Science, 20: 822-844.
• [33] GEOFFRION AM. 1972. Generalized Benders decomposition. Journal of optimization Theory and Applications, 10(4): 237-260.
• [34] GERLA M. 1973. The Design of Store-and-Forward (S/F) Networks for Computer Communications PhD thesis, University of California, Los Angeles, EUA.
• [35] GILLEN D & LEVINSON D. 1999. The full cost of air travel in the california corridor. In: PRESENTED IN THE 78^th ANNUAL MEETING OF TRANSPORTATION RESEARCH BOARD, Washington, DC, January, 10-14.
• [36] GROVE GP & O'KELLY ME. 1986. Hub networks and simultated schedule delay. Papers of the Regional Science Association, 59: 103-119.
• [37] HAMACHER HW, LABBÉ M, NICKEL S & SONNEBORN T. 2004. Adapting polyhedral properties from facility to hub location problems. Discrete Applied Mathematics, 145: 104-116.
• [38] HERAGU SS & CHEN J-S. 1998. Optimal solution of cellular manufacturing system design: Benders decomposition approach. European Journal of Operational Research, 107: 175-192.
• [39] HOANG HH. 2005. Topological optimization of networks: A nonlinear mixed integer model employing generalized benders decomposition. IEEE Transations on Automatic Control, AC-27: 164-169.
• [40] HOCHREITER R & PFLUG GC. 2007. Financial scenario generation for stochastic multistage decision processes as facility location problems. Annals of Operations Research, 152: 257-272.
• [41] KARA BY & TANSEL BC. 2003. The latest arrival hub location problem. Management Science, 47: 1408-1420.
• [42] KLINCEWICZ JG. 1998. Hub location in backbone/tributary network design: a review. Location Science, 6: 307-335.
• [43] KUBY MJ & ROBERT G. 1993. The hub network design problem with stopovers and feeders: The case of federal express. Transportation Research A, 27: 01-12.
• [44] LABBÉ M, YAMAN H & GOURDIN E. 2005. A branch and cut algorithm for hub location problems with single assignment. Mathematical Programming, 102: 371-405.
• [45] LAPORTE G, LOUVEAUX FV & MERCURE H. 1989. Models and exact-solutions for a class of stochastic location-routing problems. European Journal OF Operational Research, 39(1): 71-78.
• [46] LAPORTE G, LOUVEAUX FV & VANHAMME L. 1994. Exact solution to a location problem with stochastic demands. Transportation Science, 28(2): 95-103.
• [47] LASDON LS. 1970. Optimization Theory for Large Systems Macmillan, London.
• [48] LOUVEAUX FV & PEETERS D. 1992. A dual-based procedure for stochastic facility location. Operations Research, 40(3): 564-573.
• [49] MAGNANTI TL, MIRCHANDANI P & WONG RT. 1986. Tailoring Benders decomposition for uncapacitated network design. Mathematical Programming Study, 26: 112-154.
• [50] MAGNANTI TL & WONG RT. 1981. Accelerating benders decomposition: Algorithmic enhancement and model selection criteria. Operations Research, 29(3): 464-483.
• [51] MAHEY P, BENCHAKROUN A & BOYER F. 2001. Capacity and flow assignment of data networks by generalized benders decomposition. Journal of Global Optimization, 20(2): 169-189.
• [52] MAHEY P, OUOROU A, LEBLANC L & CHIFFLET J. 1998. A new proximal decomposition algorithm for routing in telecommunication networks. Networks, 31: 227-238.
• [53] MARIANOV V & SERRA D. 2003. Location models for airline hubs behaving as m/d/c queues. Computers and Operations Research, 30: 983-1003.
• [54] MARÍN A. 2005. Formulating and solving splittable capacitated multiple allocation hub location problems. Computer and Operations Research, 32: 3093-3109.
• [55] MARÍN A, CÁNOVAS L & LANDETE M. 2006. New formulations for the uncapacitated multiple allocation hub location problem. European Journal of Operational Research, 172: 274-292.
• [56] MCDANIEL D & DEVINE M. 1977. A modified Benders partitioning algorithm for mixed integer programming. Management Science, 24(3): 312-319, November.
• [57] MIRANDA G JR. 2004. Facility Location and Network Design with Congestion Costs and Interdependency PhD thesis, Universidade Federal de Minas Gerais, Departamento de Ciência da Computação.
• [58] O'KELLY ME. 1986. The location of interacting hub facilities. Transportation Science, 20: 92-106.
• [59] O'KELLY ME. 1987. A quadratic integer program for the location of interacting hub facilities.European Journal of Operational Research, 32: 393-404.
• [60] O'KELLY ME & BRYAN DL. 1998. Hub location with flow economies of scale. Transportation Research Part B, 32(8): 605-616.
• [61] O'KELLY ME, BRYAN DL, SKORIN-KAPOV D & SKORIN-KAPOV J. 1996. Hub network design with single and multiple allocation: A computational study. Location Science, 4(3): 125-138.
• [62] O'KELLY ME & MILLER HJ. 1994. The hub network design problem. Journal of Transport Geography, 2: 31-40.
• [63] OUOROU O, LUNA HPL & MAHEY P. 2001. Multicommodity network expansion under elastic demands. Optimization and Engineering, 2: 277-292.
• [64] RANDAZZO CD & LUNA HP. 2001. A comparison of optimal methods for local access uncapacitated network design. The Annals of Operations Research, 106: 263-286.
• [65] RODRIGUEZ V, ALVAREZ MJ & BARCOS L. 2007. Hub location under capacity constraints. Transportation Research Part E, 43: 495-505.
• [66] RODRIGUEZ-MARTIN I & SALAZAR-GONZALEZ JJ. 2008. Solving a capacitated hub location problem. European Journal of Operational Research, 184: 468-479.
• [67] ROSING KE & HODGSON MJ. 2002. Heuristic concentration for the p-median: an example demonstrating how and why it works. Computers and Operations Research, 29: 1317-1330.
• [68] ROSING KE & REVELLE CS. 1997. Heuristic concentration: two stage solution construction. European Journal of Operational Research, 97: 75-86.
• [69] SNYDER LV & DASKIN MS. 2006. Stochastic p-robust location problems. IIE Transactions, 38(11): 971-985.
• [70] SNYDER LV, DASKIN MS & TEO CP. 2007. The stochastic location model with risk pooling. European Journal of Operational Research, 179(3): 1221-1238.
• [71] SRIDHAR VV & PARK JS. 2000. Benders-and-cut algorithm for fixed-charge capacitated network design problem. European Journal of Operational Research, 125: 622-632.
• [72] YAMAN H & CARELLO G. 2005. Solving the hub location problem with modular link capacities. Computers and Operations Research, 32: 3227-3245.
• [73] YAMAN H, KARA BY & TANSEL BC. 2007. The latest arrival hub location problem for cargo delivery systems with stopovers. Transportation Research Part B, 41: 906-919.

*

Corresponding author

Publication Dates

• Publication in this collection
  05 Aug 2011
• Date of issue
  Aug 2011

History

• Accepted
  17 Mar 2011
• Received
  01 Feb 2010