Open-access FASTER MIXED-INTEGER QUADRATIC CONIC PROGRAMS FOR TWO COMPETITIVE MULTIPLE ALLOCATION P-HUB LOCATION PROBLEMS

ABSTRACT

We present new mixed-integer quadratic conic formulations for two variants of the multiple-allocation p-hub location problems in a competitive environment. The problems consist in locating p hubs so that an entrant company can establish its hub-and-spoke network to provide transportation services for pairs of origin-destination that exchange flows in a competitive market. The objective is to maximize the entrant’s market share when compared to its competitors. Both problems assume that the paths used to route the flows have one or at most two hubs. However, whereas the first problem allows an origin-destination to be serviced by multiple routes, the second problem requires that a single path be used. Here, we show that instead of maximizing the entrant’s market share, it is computationally more interesting to minimize the market lost so that equivalent, but more suitable programs to conic solvers can be obtained. When solved by a commercial conic programming solver, our proposed formulations achieve average speedups of 89 times for the multi-path variant and 37 times for the single-path variant, as shown in our extensive computational experiments on solving well-known datasets. Therefore, this work not only reformulates the problem but also substantially outperforms all prior works, demonstrating the practical applicability and effectiveness of our approaches.

Keywords:
hub location; competitive facility location; mixed-integer quadratic conic programming; mixed-integer non-linear optimization

1 INTRODUCTION

We address two variants of the multiple-allocation p-hub location problem (O’Kelly, 1986, 1987) in a competitive market. Generally speaking, the problems consist of locating p hub facilities to design a hub-and-spoke network for an entrant company that wishes to offer transportation services in a competitive market where pairs of origin-destination points that exchange flows can choose the services of different providers. This problem encapsulates a strategic decision-making process, where the objective is not merely to optimize network efficiency, but to do so in a competitive environment, in which customer choice is influenced by factors such as price, travel time, and service attractiveness. Such scenarios are common-place in industries such as air transportation, freight logistics, and telecommunications, where the design of an efficient hub-and-spoke network becomes a decisive factor for competitiveness and market share acquisition.

Hubs are facilities responsible for consolidating flows from different origins before routing them either to their final destinations or to another hub that will then distribute these flows to their final destinations. Bundling flows in interhub connections allows the entrant company to realize scale economies, which can lead to lower transportation costs than directly connecting the origin-destination points (O’Kelly, 1986; O’Kelly & Bryan, 1998). Traditionally, it is assumed that any path connecting a pair of origin-destination points has one or at most two hubs, although there are other variants of hub location problems that completely disregard this assumption (Contreras & O’Kelly, 2019).

Here, different companies compete non-cooperatively for the market by offering transportation services with varying relative utilities, e.g., price, travel time, service and transportation convenience, or company attractiveness. Customers usually choose the most appealing transportation service based on these relative utilities (Huff, 1964, 1966). The entrant company must then take these characteristics into account during the design of its hub-and-spoke network if it wants to maximize the captured market share (Eiselt & Marianov, 2009; Tiwari et al., 2021a,b). Note that, generally, companies focus on their individual gains, which are greatly affected by the location of their hubs. The success of an entrant company in capturing a larger market share is substantially influenced by the attractiveness of its hubs and their ability to sway customer behavior (Drezner & Drezner, 2001).

In the literature, this problem has been modeled with two different assumptions regarding how origin-destination points are connected to the hub network. Eiselt & Marianov (2009) and Tiwari et al. (2021a) have assumed that every pair of origin-destination points can interact with every installed hub of the network, thus allowing a multitude of paths to serve any pair of origin-destination points. We have dubbed this problem the multi-path multiple allocation p-hub location problem (MPMAPHLP). Tiwari et al. (2021b) have criticized this assumption, arguing that it prevents the aggregation of flows at the installed hubs and, therefore, hinders the exploitation of scale economies in inter-hub connections.

Tiwari et al. (2021b) have then proposed the assumption that only a single path should be used to route the flows of a given pair of origin-destination points to make it resemble the traditional multiple allocation variants of the hub location literature (Campbell & O’Kelly, 2012). We have named this variant the single path multiple allocation p-hub location problem (SPMAPHLP). Regardless of the adopted assumption, the problem gives rise to a difficult mixed-integer non-linear program that challenges existing off-the-shelf solvers. Thus, to solve such a problem, one needs to rely on different solution strategies.

For the problem involving the MPMAPHLP, Eiselt & Marianov (2009) have proposed a multi-start two-phase metaheuristic. In the first phase, for a given number of iterations, a random solution is generated and then improved by means of a local search that swaps an installed hub for a closed one. Whenever an improved solution is found, it is saved in a fixed-size list containing the overall best solutions from all iterations. The hubs of the solutions in this list form a concentration set that is used in the second phase of their algorithm. Again, for a given number of iterations, p hubs are randomly selected from the concentration set to form a solution to be improved using a local search procedure that now swaps two installed hubs for two closed ones belonging to the concentration set. The overall best solution is saved and returned at the end of the method.

For their turn, Tiwari et al. (2021a) have proposed four different approaches to solve the MPMAPHLP. Tiwari et al. have devised a Kelley cutting-plane algorithm, a reformulation based on a mixed-integer quadratic conic program (MIQCP), and Lagrangian relaxation-based algorithms using their first two methods as solvers for the Lagrangian subproblems. Tiwari et al. have made the source code for their implementation of their methods available online. Unfortunately, during the coding of their methods, Tiwari et al. typed the variable representing the scale economies in the inter-hub connections as an integer instead of setting it to float or double, thereby rendering their results biased or less insightful. For a given instance, the induced hub-and-spoke network was the same for the different tested values for the scale economy, as reported in their computational experiment tables, preventing a proper analysis of their reported results.

Tiwari et al. (2021b) have also proposed the same four previous solution strategies to address the SPMAPHLP. Once again, Tiwari et al. have typed their variable to represent the scale economy in their computer source code as an integer variable instead of a float or double, thus leading to biased or questionable results. Moreover, though they have assumed that only one path is used to route the flows for a given pair of origin-destination points, thus requiring a binary variable to represent whether a path is used or not, Tiwari et al. have relaxed this requirement by setting these variables to be non-negative within zero and one, which allowed their procedures to return dual bounds for the problem.

In their computational experiments, Tiwari et al. (2021a,b) have shown that their MIQCPs were competitive compared to their other devised approaches for solving instances with different arrangements for the number of hubs and scale economies in the Civil Aeronautics Board (CAB) dataset, using CPLEX as their conic solver set to run for one or two hours, respectively. Their MIQCPs have proven to be interesting tools for designing hub-and-spoke networks.

Here, we show a better way to reformulate both the MPMAPHLP and SPMAPHLP problems to obtain equivalent ones, but with the objective function of minimizing the market loss instead of maximizing the captured market share. This simple yet crucial approach leads to more lean and mean MIQCPs that have attained speed-ups of 89 and 37 times on average for the MPMAPHLP and SPMAPHLP problems, respectively, substantially outperforming the MIQCPs of Tiwari et al. (2021a,b) when solving the same set of computational experiments.

The remainder of this work is organized as follows. Section 2 presents the adopted notation and definitions and formally describes both problems. Section 3 underscores the differences between the proposed MIQCPs and those in the literature; while Section 4 discusses the attained numerical experiments. Finally, Section 5 concludes with our final remarks and suggestions for future work.

2 NOTATION, DEFINITIONS AND FORMULATIONS

Both MPMAPHLP and SPMAPHLP problems use the same underlying notation and definitions. Let N be the set of points scattered in a geographic region, in which each point iN is both a hub candidate and an origin or a destination for some demands. Let W={(i,j)N×N:ij} be the set of pairs of origin-destination points , with each pair (i, j) ∈ W exchanging w ij > 0 units of flow. Let I={(k,l)N×N} also be the set of all possible interhub connections or hub-pairs. There is also a set C of competitors that are already operating with a hub-and-spoke structure and competing for their respective market share, and an entrant company e that wishes to enter this market.

Although demand information is common knowledge for all competitors, no other information or infrastructure is shared between them; i.e., the companies do not cooperate. Moreover, each company aC ∪ {e} must operate its own sets H aN and Ia={(k,l)Ha×Ha} of installed hubs and active interhub connections, respectively, that may be used to form paths to route the captured demands. Recall that each path has one or two hubs, and that whenever two different hubs are used to transport the bundle of demand flows, scale economies can be exploited. Here, scale economies are represented by the scalar parameter α ∈ [0, 1].

Any entrant company seeks to capture the maximum market share possible by installing p hubs to form paths that can attract most of the demand from the origin-destination points (i, j) ∈ W . Note that each path iklj, with (i, j) ∈ W and (k, l) ∈ I a , operated by a company aC ∪ {e}, has a perceived utility uiklja can be written as for the demand w ij , which follows a probabilistic choice model that can be written as

u i k l j a = A k l a γ ( T i k l j a ) β + ( 1 - γ ) ( B i k l j a ) δ .

Akla represents the attractiveness of the interhub connection (k, l) ∈ I a and its respective hubs, while Biklja and Tiklja are the cost and travel time, respectively, of the path. Parameters β and δ are scalar values that estimate the impact of travel time and cost, respectively, on the selection behavior of the demands. Higher power values indicate a dislike for longer travel times or higher transportation costs. Parameter γ represents the weight of the relative importance of each factor (cost and time) in the utility of the path. Moreover, parameters Biklja and Tiklja are equal to c ik + αc kl + c lj and t ik +t kl +t lj , respectively, for all (i, j) ∈ W , (k, l) ∈ I, and aC ∪{e}, where c ij and t ij , i, jN are the unitary transportation cost and the travel time between points i, j, respectively.

The utility parameters are then used in an attraction function fije that returns the amount of demand w ij for the origin-destination pair (i, j) ∈ W that the entrant company e will capture with its available paths formed by its set I e of hub-pairs or interhub connections or:

f i j e = w i j ( k , l ) I e u i k l j e a C { e } ( k , l ) I a u i k l j a = w i j ( k , l ) I e u i k l j e ( k , l ) I e u i k l j e + U i j c

in which Uijc=aC(k,l)Iauiklja is introduced to simplify notation hereafter.

Note that while sets H a and I a are known for competitors aC, sets H e and I e that maximize the total captured demand must be found for the entrant company e. Note also that the cardinality of set H e is equal to p, since exactly p hubs will be installed for the entrant company e. Furthermore, to find these sets or the underlined configuration of hubs and active interhub connections, function fije needs to be properly adjusted to comply with the respective assumptions of the MPMAPHLP and SAMAPHLP problems, which, given the nature of the function fije, are initially modeled as mixed-integer non-linear programs (MINLP). Table 1 summarizes the main notation. However, additional definitions for auxiliary parameters and decision variables might be required when recasting the current mixed-integer nonlinear programs into conic programs in Section 3.

Table 1
Summary of the main notation.

2.1 The maximization form for the MPMAPHLP problem

In the MPMAPHLP problem, an entrant company e wants to enter a market for the transportation of flows between pairs of origin-destination points that seek the most attractive service among competing companies. The entrant company e plans to maximize its captured market share by installing p hubs to form its hub-and-spoke network. It is assumed that each competitor can use multiple paths that go through different hubs or hub-pairs to serve any pair of points (i, j) ∈ W.

To formulate the MPMAPHLP, Eiselt & Marianov (2009) have defined the decision variables of location y k ∈ {0, 1}, which assume a value of one if a hub is installed at a point kN, and zero otherwise; and the activation of hub-pairs or interhub connections ω kl ∈ {0, 1}, which assume a value of one if the interhub connection or hub-pair (k, l) ∈ I is set, and zero otherwise. Note that the decision variables ω kl assemble the set I e of the entrant company. With the aid of these decision variables, Eiselt & Marianov have adjusted the attraction function for the pair of origin-destination points (i, j) ∈ W as

f i j ( ω ) = w i j ( k , l ) I u i k l j e ω k l ( k , l ) I u i k l j e ω k l + U i j c

and then modeled the problem as follows:

MAXMPMAPHLP max ( i , j ) W w i j ( k , l ) I u i k l j e ω k l ( k , l ) I u i k l j e ω k l + U i j c 1 s.t.: ω k l y k ( k , l ) I 2 ω k l y l ( k , l ) I 3 k N y k = p 4 y k { 0 , 1 } k N 5 ω k l { 0 , 1 } ( k , l ) I 6

The objective function (1) maximizes the captured demand of the entrant company. Constraints (2) and (3) are the activation constraints. The hub-pair or interhub connection (k, l) ∈ I can only be active if the respective hubs are installed. Constraint (4) guarantees that p hubs will be installed. The domain of the decision variables are shown in constraints (5) and (6).

Note that the MAXMPMAPHLP formulation is an integer nonlinear program (INLP) that can be transformed into a mixed-integer nonlinear program (MINLP) by relaxing the decision variables ω kl ∈ [0, 1], (k, l) ∈ I. As the attractiveness function is a concave, smooth, and differentiable function, and since we wish to maximize the captured market share, whenever both y k and y l assume the value of one, then w kl will also assume the value of one. See constraints (2) and (3). Nevertheless, whenever either y k or y l is equal to zero, w kl assumes the value of zero. Hence, we shall replace constraints (6) hereafter with their relaxations.

2.2 The maximization form for the SPMAPHLP problem

The formulation MAXMPMAPHLP for the MPMAPHLP problem assumes that every point is connected to every hub, thereby enabling a large number of paths to serve any pair of origin-destination points that exchange flows. Tiwari et al. (2021b) have then proposed a variant that restricts the number of paths to a single one, giving rise to the single-path multiple allocation p-hub location (SPMAPHLP). Tiwari et al. have then adapted the attractiveness function using decision variables x iklj ∈ {0, 1} that assume the value of one if the path (i, k, l, j), (i, j) ∈ W and (k, l) ∈ I, is used to route the flow demand w ij and zero otherwise. The attractiveness function is written as

f i j ( x ) = ( k , l ) I u i k l j e x i k l j ( k , l ) I u i k l j e x i k l j + U i j c

for (i, j) ∈ W . Tiwari et al. have then embedded the attractiveness function in the formulation of Hamacher et al. (2004) proposed for the multiple allocation hub location problem, resulting in the following program.

MAXSPMAPHLP max ( i , j ) W w i j ( k , l ) I u i k l j e x i k l j ( k , l ) I u i k l j e x i k l j + U i j c 7 s.t.: ( k , l ) I x i k l j = 1 ( i , j ) W 8 l N x i k l j + l N { k } x i l k j y k ( i , j ) W , k N 9 k N y k = p 10 y k { 0 , 1 } k N 11 x i k l j { 0 , 1 } ( i , j ) W , ( k , l ) I 12

The objective function (7) maximizes the total captured demand of the entrant company. Constraints (8) impose that only one path will be available for the demand of the pair of origin-destination points (i, j) ∈ W . Constraints (9) ensure that only paths with installed hubs are available. As before, constraint (10) guarantees that p hubs will be installed. Constraints (11) and (12) show the domain of the decision variables.

To avoid branching on the x variables, Tiwari et al. (2021b) have relaxed them in their computation experiments, thereby attaining dual bounds for the problem. Even though the objective function is a concave, smooth, and differentiable function, and the problem has a maximization objective function, there is no guarantee that integer-valued x will be obtained when their domain is relaxed.

Though constraints (8) are assignment constraints, they actually act as knapsack constraints in the formulation, given the maximization sense of the objective function. Therefore, relaxing the x variables may lead to a solution that has them fractional at the end of a solution procedure. To be consistent with the computational experiments ofTiwari et al. (2021b), we will replace constraints (8) with their relaxed counterparts hereafter.

Moreover, Marín et al. (2006) have shown that the formulation of Hamacher et al. (2004) can produce optimal solutions with routes having more than two hubs when the triangular inequalities do not hold. Instead of using the constraints proposed by Marín et al. to avoid these potentially unfeasible paths, we favor the following constraints.

l N : l j x i i l j y i - y j ( i , j ) W (13)

k N : k i x i k j j y j - y i ( i , j ) W (14)

x i i j j y i + y j - 1 ( i , j ) W (15)

Constraints (13) forbid the origin-destination pair (i, j) from flowing through two other hubs once the origin i is a hub, while the destination j is not. Constraints (14) forbid the origin-destination pair (i, j) from flowing through two other hubs once the destination j is a hub, but the origin i is not. Constraints (15) enforce the flow through arc (i, j) once both the origin i and destination j are hubs.

Note that:

  • yi = 1 and y j = 0 yield constraints (13) active while constraints (14) and (15) are rendered inactive.

  • yi = 0 and y j = 1 yield constraints (14) active while constraints (13) and (15) are rendered inactive.

  • yi = 1 and y j = 1 yield constraints (15) active while constraints (13) and (14) are rendered inactive.

All three of these constraint families are rendered ineffective if none of the above occurs.

Herein, formulations based only on the formulation of Hamacher et al. (2004) are called original, whereas those using the above constraints are labeled fixed. We have analyzed both variants in our computational experiments to evaluate the impact of these additional constraints. Constraints (13)-(15) have been successfully applied to render the necessary corrections in de Camargo et al. (2009), and we prefer using them since the (i, j) (origin-destination) pair decomposition is preserved. This decomposition is a clue for any technique exploiting problem structure and might be the subject of future research.

2.3 Our minimization form for the MPMAPHLP and SPMAPHLP problems

Both MPMAPHLP and SPMAPHLP problems have been modeled as maximization problems. Nevertheless, it is possible to model them as equivalent problems that minimize the loss of market share. This simple sense modification results in a mathematical structure that has many advantages, as it leads to lean and mean second-order quadratic conic programs, as we will see in the next section.

We know that the sum of the captured market and the loss of market share is equal to the total demand. So we can model both MPMAPHLP and SPMAPHLP problems as the minimization of the loss of market share. Note that, for any pair of origin-destination points (i, j) ∈ W in the MPMAPHLP problem, we have the following relation:

f i j ( ω ) + ϕ i j ( ω ) = w i j ϕ i j ( ω ) = w i j - w i j ( k , l ) I u i k l j e ω k l ( k , l ) I u i k l j e ω k l + U i j c = w i j 1 - ( k , l ) I u i k l j e ω k l ( k , l ) I u i k l j e ω k l + U i j c = w i j U i j c ( k , l ) I u i k l j e ω k l + U i j c

in which φ i j (ω) is the loss of market share. After dividing the right-hand side by Uijc, we get

ϕ i j ( ω ) = w i j 1 ( k , l ) I π i k l j ω k l + 1

in which πiklj=uikljeUijc (i, j) ∈ W , and (k, l) ∈ I. We can then model the MPMAPHLP problem as follows:

MINMPMAPHLP min ( i , j ) W w i j 1 ( k , l ) I π i k l j e ω k l + 1 s.t.: 2 - 5 ω k l [ 0 , 1 ] ( k , l ) I

Applying the same rationality to the SPMAPHLP problem, we obtain the following formulation:

MINSPMAPHLP min ( i , j ) W w i j 1 ( k , l ) I π i k l j x i k l j + 1 s.t.: 8 - 11 x i k l j [ 0 , 1 ] ( i , j ) W , ( k , l ) I

All formulations, MAXMPMAPHLP, MAXSPMAPHLP, MINMPMAPHLP, and MINSPMAPHLP, are mixed-integer nonlinear programs that cannot be directly solved by most off-the-shelf solvers. Nevertheless, they are amenable to being recast into conic programs, which can be solved by solvers like Gurobi and CPLEX.

3 MIXED-INTEGER QUADRATIC CONIC PROGRAMS

Although conic modeling still yields a nonlinear program, the resulting formulation has a special type of mathematical structure that can be solved by off-the-shelf solvers. Conic optimization optimizes a linear function over conic inequalities (Ben-Tal & Nemirovski, 2001), having such an approach an easy computational implementation and many successful applications (Benson & Ü mit Sag˘lam, 2014). Here, we will streamline how Tiwari et al. (2021a,b) have recast their MAXMPMAPHLP and MAXSPMAPHLP formulations, and then we will show how we have recast our MINMPMAPHLP and MINSPMAPHLP programs.

3.1 Recasting MAXMPMAPHLP formulation into a conic program

Using Atamtürk et al. (2012)’s guidelines, Tiwari et al. (2021a) have recast the MAXMPMAPHLP problem into a mixed-integer quadratic conic program after introducing auxiliary variables Q ij ≥ 0, V ij ≥ 0, and (i, j) ∈ W, in which Qij=1(k,l)Iuiklje ωkl+Uijc and Vij=Qij (k,l)Iuiklje ωkl. They have reworked the former relation as Qij(k,l)Iuiklje ωkl+QijUijc=1 to obtain a new relation Vij+QijUijc=1, (i, j) ∈ W . Returning yet to the first relation, Tiwari et al. have introduced other decision variables r ij ≥ 0 to replace its denominator or rij=(k,l)Iuiklje ωkl+Uijc, (i, j) ∈ W , which allowed them to write the inequality Qij1rij or Q ij r ij ≥ 1. Observing the second relation, they have rewritten it as Vij(k,l)Iuiklje ωklrij or Vijrij(k,l)Iuiklje ωkl.

As the decision variables ω kl were originally defined as binary variables, the last two relations can be rewritten as second-order rotated conic constraints and then transformed into standard second-order constraints or (Q ij + r ij )2 ≥ 2 + Q ij 2 + r ij 2 and (Vij+rij)22(k,l)Iuiklje ωkl2+Vij2+rij2, for all (i, j) ∈ W . Gathering all these relations, Tiwari et al. (2021a) have proposed the following mixed-integer second-order conic program for the MAXMPMAPHLP formulation:

CPMAXMPMAPHLP max ( i , j ) W w i j V i j s.t.: 2 - 5 V i j + U i j c Q i j = 1 ( i , j ) W r i j = ( k , l ) I u i k l j e ω k l + U i j c ( i , j ) W ( V i j + r i j ) 2 2 ( k , l ) H e u i k l j e ω k l 2 + V i j 2 + r i j 2 ( i , j ) W ( Q i j + r i j ) 2 2 + Q i j 2 + r i j 2 ( i , j ) W Q i j , V i j , r i j 0 ( i , j ) W ω k l [ 0 , 1 ] ( i , j ) I

3.2 Recasting MAXSPMAPHLP formulation into a conic program

Using the same rationale in Tiwari et al. (2021a), Tiwari et al. (2021b) have performed the same steps for the MAXSPMAPHLP formulation to obtain the following conic program:

CPMAXSPMAPHLP max ( i , j ) W w i j V i j s.t.: 8 - 11 V i j + U i j c Q i j = 1 ( i , j ) W r i j = ( k , l ) I u i k l j e x i k l j + U i j c ( i , j ) W ( V i j + r i j ) 2 2 ( k , l ) H e u i k l j e x i k l j 2 + V i j 2 + r i j 2 ( i , j ) W ( Q i j + r i j ) 2 2 + Q i j 2 + r i j 2 ( i , j ) W Q i j , V i j , r i j 0 ( i , j ) W x i k l j [ 0 , 1 ] ( i , j ) W , ( k , l ) I

3.3 Recasting MINMPMAPHLP formulation into a conic program

We make use of the rotated cone in our reformulation, which takes the following form xyz ⇐⇒ ||(2x, yz)||2y + z, where || · ||2 is the L 2 norm and x, y, z are decision variables of suitable sizes. Let the denominator of the attractiveness function be equal to r ij = ∑(k,l)∈I π iklj ω kl + 1; then we define Vij=1rij. Note that r ij ≥ 1. Now, we define Vij=1rij or the rotated cone inequality V ij r ij ≥ 1, for (i, j) ∈ W. Therefore, the MINMPMAPHLP problem is equivalent to the following conic program:

CPMINMPMAPHLP min ( i , j ) W w i j V i j s.t.: 2 - 5 r i j = ( k , l ) I π i k l j ω k l + 1 ( i , j ) W V i j r i j 1 ( i , j ) W r i j 1 ( i , j ) W 0 V i j 1 ( i , j ) W ω k l [ 0 , 1 ] ( i , j ) I

The objective function is linear, and the constraints are now either linear or conic quadratic, making it a general mixed-integer conic quadratic program suitable for solving with solvers like Gurobi and CPLEX.

3.4 Recasting MINSPMAPHLP formulation into a conic program

Repeating the same steps, but now for the MINSPMAPHLP problem, we obtain the following conic program:

CPMINSPMAPHLP min ( i , j ) W w i j V i j s.t.: 8 - 11 and 13 - 15 r i j = ( k , l ) I π i k l j x i k l j + 1 ( i , j ) W V i j r i j 1 ( i , j ) W r i j 1 ( i , j ) W 0 V i j 1 ( i , j ) W x i k l j [ 0 , 1 ] ( i , j ) W , ( k , l ) I

Note that our conic programs, i.e., the CPMINMPMAPHLP and CPMINSPMAPHLP formulations, have fewer decision variables and constraints, yet possess an interesting mathematical structure suitable for the solvers. Their leaner and meaner structure explains the superior performance in the computational experiments carried out and presented in the next section.

However, it is important to stress at this point that all the presented models in this article are mixed-integer nonlinear or mixed-integer conic programs. While continuous conic programs may be efficiently solved by means of interior-point methods, the presence of integer decision variables means that the worst-case performance of all of them is non-polynomial. Hence, from a strict theoretical point of view, all the formulations deployed here as well as those already available in the literature are NP-hard.

4 COMPUTATIONAL EXPERIMENTS AND PERFORMANCE

We have conducted a comparative analysis between the results of our conic formulations and those of Tiwari et al. (2021a) and Tiwari et al. (2021b). The computational experiments were conducted using the standard datasets commonly utilized in hub location problems. The Civil Aeronautics Board (CAB) dataset encompasses data on passenger flow volumes and intercity distances for 495 US cities. Our experiments have included instances with 10, 15, 20, 25, 30, 40, and 50 nodes. The customers’ utility was calculated according to the procedure detailed in Tiwari et al. (2021b)’s study.

All computational tests were conducted on a Dell Workstation Precision 5860, featuring an Intel(R) Xeon(R) W3-2423 CPU @ 4.2GHz and 256 GB of RAM. The workstation operated in a 64-bit Linux environment. The methods were coded in C++ and utilized CPLEX Optimization Studio 22.1.1. CPLEX was configured to run sequentially using a single thread. For all runs, we have established the default termination criterion of an optimality gap of 0.01% or a total running time limit of 3600s (one hour). All source codes and instances are available at https://github.com/risargo/MIQCPMAPHLP. Note that, in all the following tables, the entries marked with ‘-’ mean unavailable/impossible-to-determine/out-of-memory, while empty entries mean values equal to zero.

The results for the multiple path problem variant are displayed in Table 2, which reports the computational times in seconds (cpu (s)) and the number of explored branch and bound nodes (#bbn) taken by the maximization formulation of Tiwari et al. (2021a) and our minimization approach to solve the instances indicated by the columns labeled n, p, and α or the size of the instances in terms of the total number of nodes, the number of hubs to be installed, and the economy of scale set on the inter-hub connections, respectively.

Table 2
This table presents a comparison of the performance of the algorithms in solving instances with sizes n = {10, 15, 20, 25, 30, 40} while varying parameters p and α accordingly. We have included the average running times in seconds (cpu), the total count of branch-and-bound nodes (#bbn), and the final gaps.

The winning formulations with respect to the total computational time have their entries in bold. Table 2 also reports the captured market share (%) and the final optimality gap (gap (%)) when not null, and the ratio between the computational times of the maximization formulation and our minimization one. After each batch of instances with the same number of nodes (n) and hubs to be installed (p), a row displays the averages of the main results for that batch.

Reviewing Table 2, it becomes evident that, for problem sizes n between 10 and 30, using the minimization model is significantly more efficient than the maximization formulation. Note that the former demonstrates a remarkable performance of being 355 times faster than the latter, while being, on average, 89 times faster.

Although it can be argued that this performance decreases with the increase in instance sizes, it is important to note that, under the given time restrictions, the maximization formulation fails to provide primal and dual bounds for all cases involving 40 nodes with p = 5 and p = 6, while the minimization formulation provides an optimality gap for all of these instances within one hour of computational time.

Finally, the benchmarking methodology proposed by Dolan & Moré (2002) is used to render Figure 1, in which the dominance of our minimization approach is further confirmed for the multiple path variant. Our minimization formulation is unequivocally superior in performance to its maximization counterpart.

Figure 1
Performance profile for the maximization and minimization formulations for the multiple path variant.

The benchmark profile graph by Dolan & Moré (2002), formally called performance profile, is a crucial methodological tool in the comparative analysis of numerical optimization solvers. It aims to provide an objective and robust assessment of the relative performance of a set of algorithms (S) on a battery of test problems (P), overcoming the limitations of simple aggregate metrics such as averages or totals.

The construction of the profile is based on the relative performance ratio τ p,s , which compares the cost of a solver s on problem p (typically CPU time or number of iterations) with the best cost obtained by any solver for that same problem. The result is a cumulative distribution function, ρ s (τ), which represents the probability that a solver s will solve a problem with a cost at most τ times higher than the best cost.

The interpretation of the performance profile reveals the efficiency and robustness characteristics of each solver. The intersection of the curve with the vertical axis at τ = 1 indicates absolute efficiency, denoting the fraction of problems for which the solver is the fastest (i.e., τ p,s = 1). As the factor τ increases on the horizontal axis, the curve rises, indicating the robustness of the solver, which is the probability that it will solve a problem at a cost reasonably close to the minimum.

The curve that visually dominates the graph-being consistently higher and further to the left-is, by definition, representative of the superior solver, as it combines the highest probability of being the best solver (efficiency) with the highest likelihood of solving problems in general (robustness). Note that, when a solver fails to solve an instance, the benchmark profile methodology assigns an infinite performance ratio to that solver for that specific problem, implying thus that the solver is unable to solve that instance. For instance, a solver that fails on half of the problems will have its performance profile curve stopped at a maximum height of 0.5 on the y-axis. This is ensured by assigning an infinite performance ratio to the failed cases.

Table 3 presents the computational results for the single-path problem introduced by Tiwari et al. (2021b). The table has the same display structure as the previous one. Note that the results also show a similar trend. For problem sizes ranging from 10 to 30, our minimization formulation performs better than its maximization counterpart, presenting slightly inferior performance in only six cases out of 150.

Table 3
This table presents a comparison of the performance of the algorithms in solving instances with sizes n = {10, 15, 20, 25, 30, 40} while varying parameters p and α accordingly. We have included the average running times in seconds (cpu), the total count of branch-and-bound nodes (#bbn), and the final gaps.

For instances with 40 nodes, our minimization model was able to optimally solve or return small optimality gaps for 17 of the 30 instances, while the maximization model failed to solve 24 instances within the given time limit and returned small optimality gaps for the others. Nonetheless, the benchmark performance profile shown in Figure 2 leaves no doubt about the superiority of our minimization model. It is important to note that both the original single-path maximization and minimization models did not take into account the concerns raised by Marín et al. (2006) who pointed out that the Hamacher et al. (2004)’s formulation can redeem incorrect hub-and-spoke networks when constraints like (13)-(15) are not added to the formulation.

Figure 2
Performance profile for the maximization and minimization formulations for the single path original variant.

Hence, we have assessed the impact of these constraints on the single-path maximization and minimization models, reporting their computational results in Table 4. Recall that constraints (13)-(15) prevent the formation of pathological paths with more than two hubs on them. As can be seen in Table 4, these fixing constraints severely impact both formulations, substantially increasing the computational effort required to solve the same set of instances.

Table 4
This table presents a comparison of the performance of the algorithms in solving instances with sizes n = {10, 15, 20, 25, 30, 40} while varying parameters p and α accordingly. We have included the average running times in seconds (cpu), the total count of branch-and-bound nodes (#bbn), and the final gaps.

On the one hand, Tiwari et al. (2021b)’s formulation performance seems to be more affected for problem sizes ranging from 10 to 30, with the impact being more pronounced for the latter size, as the established time limit is reached more often before finding an optimal solution. On the other hand, analyzing the 40 node instances, it is hard to define a winning model, as both seem to be tied. Unfortunately, this effect is not observed in the performance profile shown in Figure 3, as it clearly favors our minimization model.

Figure 3
Performance profile for the maximization and minimization formulations for the single path fixed variant.

One may argue that the given time limit may interfere with the results, especially for p = 5, p = 6, and n = 40, since these results do not clearly differ from one model to the other, making it hard to declare a winning model. Therefore, we have designed an additional computational experiment in which we have removed the imposition of a time limit while solving an instance with n = 50 and p = 5. The idea is to break the ties and further investigate both behavior of both formulations’ without the interference of a time limit while addressing larger problems.

Table 5 reports the attained results. While clearly demonstrating the dominance of the minimization model once again, since the maximization model of Tiwari et al. (2021b) fails to render bounds for all instances due to a lack of memory, meaning that it would take more than 256 GB of RAM to run the model properly, our minimization formulation was able to reach proven optimality for all instances in under 6.5 hours.

Table 5
This table presents a comparison of the performance of the algorithms in solving instances with size n = {50} with parameter p = 5 and parameter varying α accordingly. All time limits removed. We have included the average running times in seconds (cpu), the total count of branch-and-bound nodes (#bbn), and the final gaps.

Please recall the highlighted difference between empty entries versus entries marked with ’-’, when not even the continuous bounds are available. Nevertheless, it could be argued that the new formulations are still very resource-intensive, considering the memory usage and the solution times achieved. This is a clear indication that the proposed formulations would require a decomposition effort in order to improve their scalability. This might be a direction for future research, provided that for a fixed hub structure, the origin-destination pairwise decomposition is preserved by constraints (13) - (15), easily allowing for a Benders’ decomposition scheme, for instance.

5 CONCLUSION AND FINAL REMARKS

In conclusion, our proposed models offer a simpler, more robust, and efficient approach to both problem variants with single and multiple paths, outperforming the previous models proposed by Tiwari et al. (2021a) and Tiwari et al. (2021b). Instead of maximizing the market share of an entrant company, we have adopted the strategy of minimizing the loss of market share. This has represented a significant advancement over traditional assumptions.

Future research could explore more advanced methodologies for handling large-scale problems, such as the generalized Benders Decomposition (Geoffrion, 1972; de Camargo et al., 2009; Miranda Jr et al., 2011; Miranda Jr. et al., 2011; de Camargo & Miranda Jr, 2012; de Camargo & Miranda, 2012) or the Outer Approximation technique combined with the Benders’ Decomposition method (de Camargo et al., 2011) to further improve solution efficiency. In addition, investigating sources of uncertainty in problem data, such as variability in transportation flows, could make the model more realistic and adaptable to practical scenarios.

Finally, despite the growing body of articles in this field, there is a lack of studies that include the topic of environmental sustainability as a factor in the attractiveness of the choice of transport. Furthermore, the present study have focused on a single-objective market-share formulation, but real-world hub-network design is inherently multi-objective. Investigating extensions that simultaneously maximize the captured flow, minimize transportation costs, carbon emissions, and embed richer user-preference models as a way to optimize service attractiveness represents a promising avenue for future research.

Acknowledgements

The authors thank CNPq for its financial support, grants 305478/2022-2 and 423065/2021-1.

Data Availability

All the source codes and instances used to produce the results of this article are available in GitHub at https://github.com/risargo/MIQCPMAPHLP or in Open Science Framework (OSF) at https://osf.io/wq6jt/?view only=78700465d980479d897b6f8082740d21. The data might also be available from the authors upon request.

References

  • ATAMTÜRK A, BERENGUER G & SHEN ZJ. 2012. A conic integer programming approach to stochastic joint location-inventory problems. Operations Research, 60(2): 366-381.
  • BENTAL A & NEMIROVSKI A. 2001. Lectures on modern convex optimization: analysis, algorithms, and engineering applications. SIAM.
  • BENSON HY & ÜMIT SAĞLAM. 2014. Mixed-Integer Second-Order Cone Programming: A Survey. chap. Chapter 2, pp. 13-36. Informs. Available at: https://pubsonline.informs.org/doi/abs/10.1287/educ.2013.0115
    » https://pubsonline.informs.org/doi/abs/10.1287/educ.2013.0115
  • CAMPBELL JF & O’KELLY ME. 2012. Twenty-five years of hub location research. Transportation Science, 46(2): 153-169.
  • CONTRERAS I & O’KELLY ME. 2019. Hub Location Problems. pp. 327-363. Cham: Springer International Publishing. Available at: https://doi.org/10.1007/978-3-030-32177-2_12
    » https://doi.org/10.1007/978-3-030-32177-2_12
  • DE CAMARGO RS, DE MIRANDA JR G & FERREIRA RP. 2011. A hybrid outer-approximation/benders decomposition algorithm for the single allocation hub location problem under congestion. Operations Research Letters, 39(5): 329-337.
  • DE CAMARGO RS & MIRANDA G. 2012. Single allocation hub location problem under congestion: Network owner and user perspectives. Expert Systems with Applications, 39(3): 3385-3391.
  • DE CAMARGO RS & MIRANDA JR G. 2012. Addressing congestion on single allocation hub-and-spoke networks. Pesquisa Operacional, 32: 465-496.
  • DE CAMARGO RS, MIRANDA JR G, FERREIRA RPM & LUNA H. 2009. Multiple allocation hub-and-spoke network design under hub congestion. Computers & Operations Research , 36(12): 3097-3106.
  • DOLAN ED & MORÉ JJ. 2002. Benchmarking optimization software with performance profiles. Mathematical programming, 91: 201-213.
  • DREZNER T & DREZNER Z. 2001. A Note on Applying the Gravity Rule to the Airline Hub Problem. Journal of Regional Science, 41(1): 67-72. Available at: https://onlinelibrary.wiley.com/doi/abs/10.1111/0022-4146.00207
    » https://onlinelibrary.wiley.com/doi/abs/10.1111/0022-4146.00207
  • EISELT HA & MARIANOV V. 2009. A conditional p-hub location problem with attraction functions. Computers & Operations Research , 36(12): 3128-3135.
  • GEOFFRION A. 1972. Generalized Benders decomposition. Journal of Optimization Theory and Applications, 10: 237-260.
  • HAMACHER HW, LABBÉ M, NICKEL S & SONNEBORN T. 2004. Adapting polyhedral properties from facility to hub location problems. Discrete Applied Mathematics, 145(1): 104-116. Available at: https://www.sciencedirect.com/science/article/pii/S0166218X04000691 graph Optimization IV.
    » https://www.sciencedirect.com/science/article/pii/S0166218X04000691
  • HUFF DL. 1964. Defining and estimating a trading area. Journal of marketing, 28(3): 34-38.
  • HUFF DL. 1966. A Programmed Solution for Approximating an Optimum Retail Location. Land Economics, 42(3): 293-303.
  • 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(1): 274-292. Available at: https://www.sciencedirect.com/science/article/pii/S0377221704007994
    » https://www.sciencedirect.com/science/article/pii/S0377221704007994
  • MIRANDA JR G, DE CAMARGO RS, PINTO LR, CONCEIÇÃO SV & FERREIRA RPM. 2011. Hub location under hub congestion and demand uncertainty: the Brazilian case study. Pesquisa Operacional , 31: 319-349.
  • MIRANDA JR G, LUNA HPL, DE CAMARGO RS & PINTO LR. 2011. Tree network design avoiding congestion. Applied mathematical modelling, 35(9): 4175-4188.
  • O’KELLY ME. 1986. The location of interacting hub facilities. Transportation science, 20(2): 92-106.
  • O’KELLY ME. 1987. A quadratic integer program for the location of interacting hub facilities. European journal of operational research, 32(3): 393-404.
  • O’KELLY ME & BRYAN D. 1998. Hub location with flow economies of scale. Transportation research part B: Methodological, 32(8): 605-616.
  • TIWARI R, JAYASWAL S & SINHA A. 2021a. Alternate solution approaches for competitive hub location problems. European Journal of Operational Research, 290: 68-80.
  • TIWARI R, JAYASWAL S & SINHA A. 2021b. Competitive hub location problem: Model and solution approaches. Transportation research part B: Methodological , 146: 237-261.

Edited by

  • Editor responsible for the review
    Editor-in-Chief: Annibal Parracho Sant’Anna.

Publication Dates

  • Publication in this collection
    06 Feb 2026
  • Date of issue
    2026

History

  • Received
    02 July 2025
  • Accepted
    02 Nov 2025
location_on
Sociedade Brasileira de Pesquisa Operacional Rua Mayrink Veiga, 32 - sala 601 - Centro, 20090-050 , Tel.: +55 21 2263-0499 - Rio de Janeiro - RJ - Brazil
E-mail: sobrapo@sobrapo.org.br
rss_feed Acompanhe os números deste periódico no seu leitor de RSS
Ir para o topo Reportar erro