https://orcid.org/0000-0002-2939-4811
Figures
Abstract
Stochastic and robust optimization approaches often result in sub-optimal solutions for the uncertain p-hub median problem when continuous design parameters are discretized to form different environmental scenarios. To solve this problem, this paper proposes a triangular fuzzy number model for the Non-Strict Uncapacitated Multi-Allocation p-hub Median Problem. To enhance the quality and the speed of optimization, a novel optimization approach, combining the triangular fuzzy number evaluation index with the Genetic-Tabu Search algorithm, is proposed. During the iterations of the Genetic-Tabu Search algorithm for finding the optimal solution, the fitness of fuzzy hub schemes is calculated by considering the relative positional relationships of triangular fuzzy number membership functions. This approach directly addresses the triangular fuzzy number model and ensures the integrity of information in the p-hub problem as much as possible. It is verified by the classic Civil Aeronautics Board and several self-constructed data sets. The results indicate that, compared to the traditional Genetic Algorithm and Tabu Search algorithm, the Genetic-Tabu Search algorithm reduces average computation time by 49.05% and 40.93%, respectively. Compared to traditional random, robust, and real-number-based optimization approaches, the proposed optimization approach reduces the total cost in uncertain environments by 1.47%, 2.80%, and 8.85%, respectively.
Citation: Wang Y, Zhu T, Yuan K, Li X (2024) A fuzzy interval optimization approach for p-hub median problem under uncertain information. PLoS ONE 19(3): e0297295. https://doi.org/10.1371/journal.pone.0297295
Editor: Mazyar Ghadiri Nejad, Cyprus International University Faculty of Engineering: Uluslararasi Kibris Universitesi Muhendislik Fakultesi, TURKEY
Received: September 19, 2023; Accepted: January 2, 2024; Published: March 15, 2024
Copyright: © 2024 Wang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: the National Science Foundation of China (Grant U2033213), the Open Fund of Key Laboratory of Flight Techniques and Flight Safety of the CAAC (Grant FZ2021KF12), and the Construction Plan of Scientific Research and Innovation Team for Civil Avi-ation Flight University of China (Grant JG2022-21).
Competing interests: This study was supported by the National Science Foundation of China (Grant U2033213), the Open Fund of Key Laboratory of Flight Techniques and Flight Safety of the CAAC (Grant FZ2021KF12), the Construction Plan of Scientific Research and Innovation Team for Civil Aviation Flight University of China (Grant JG2022-24), and the Fundamental Research Funds for the Central Universities (Grant PHD2023-030, ZHMH2022-014).
1. Introduction
The p-hub median problem is one of the hub location problems that involves determining the locations of p hubs, the itineraries of all Origin-to-Destination (OD) pairs, and the link relations between points, such that the total system cost is minimized. Such kind of a problem can be further divided into Single Allocation p-Hub Median Problem (SApHMP) and Multiple Allocation p-Hub Median Problem (MApHMP) based on how hubs and non-hubs are connected [1]. Early interests in such kind of a problem were mainly focused on modeling skills. O’Kelly [2] proposed a nonlinear programming model with a quadratic objective function. Building upon the former research, Campbell [3] proposed a Non-Strict Uncapacitated Multi-Allocation p-Hub Median Problem (NSUMApHMP), which specifies that non-hubs can be interconnected and multiple paths can be selected for transportation. To efficiently address the classical p-hub center problem, Skorin-Kapov [4] developed new mixed linear formulations with a tight linear programming relaxation which effectively simplifies the model’s complexity. The aforementioned research developed several classical models of the p-hub median problem.
However, the assumption of deterministic design parameters in hub design models above is unrealistic and problematic due to inherent uncertainty, and will inevitably lead to sub-optimal selection of hub locations. To improve the anti-risk capability in uncertain environments, many scholars attempt to treat these design parameters (such as unit flow cost and demand, etc.) as random variables to generate environment scenarios. Stochastic optimization was one of the main effective approaches. It selected the hub location with the best expectation value for all scenarios. Typical representative findings are included as follows. Contreras et al. [5] solved the uncapacitated multiple allocation p-hub median problem (UMApHMP) by turning the random expectation model into a deterministic expectation model with known probability distributions. Zhai et al. [6] developed a two-stage model that serves the demand to solve the UMApHMP by minimizing the expected cost in each scenario. Ghaderi and Rahmaniani [7] modeled the uncapacitated single allocation p-hub median problem (USApHMP) with stochastic demand and travel time to obtain the optimal expectation value of the objective function in each scenario. Peiró et al. [8] developed a stochastic programming modeling framework for the r-allocation p-hub median problem with nonstop services. Based on these stochastic optimization models, heuristic algorithms are proposed to improve their ability to solve. Azizi et al. [9] proposed an exact algorithm and a Genetic Algorithm to solve the model that regards the total expected congestion cost as the objective function of USApHMP. Shang et al. [10] developed a memetic algorithm to solve MApHMP with stochastic flow and demand. Alvarez et al. [11] designed a meta-heuristic framework with Monte Carlo simulation to solve SApHMP, in which flows and demands are regarded as stochastic variables. These above researches provide feasible approaches to solve the uncertainty in the p-hub median problem.
Stochastic optimization methods typically necessitate specifying the exact probabilities for each environmental scenario, which is challenging in reality due to the inability to predefine all possible scenarios. Hence, scholars have introduced robust optimization methods, which optimize hub locations by searching for the best solutions under the worst-case scenarios, accommodating partial situations. Typical representative findings are included as follows. Ghaffarinasab [12] optimized the UMApHMP by minimizing the worst-case cost with the hybrid uncertainty set. Wang et al. [13] presented a novel two-stage distributed robust optimization model by minimizing the worst-case cost for the two-allocation p-hub median problem under the case that the demand has partially distributed information. Meraklč and Yaman [14] used the hose uncertainty set to minimize the worst-case cost within all possible demands for UMApHMP. These approaches obtain the optimal solution in the worst-case scenario by robust optimization models. To further solve the robust optimization model, Ghaffarinasab et al. [15] utilized the hybrid demand uncertainty sets to transform the robust USApHMP model into a linear MIP model, and the Tabu-based mate-heuristic solution algorithms were proposed to solve it. To construct a mean-robustness programming model to solve USApHMP, Ge [16] minimized the average of total cost and the weighted sum of risks in all scenarios. Ghaffarinasab [17] devised a precise solving algorithm based on Benders decomposition to address the UMApHMP. These stochastic optimization researches above provide feasible approaches to the p-hub median problem under uncertain cases with unknown scenarios.
The above two aforementioned approaches partially address the uncertainty in the p-hub median problem. However, due to their limited utilization of continuous parameter information and lack of specification of distribution probabilities, the generated results are likely to remain sub-optimal [18]. Many scholars attempted to use interval-valued data and fuzzy mathematical methods to solve the hub location problem. Therefore, the interval-valued numbers with upper and lower limits are introduced into the study. Ghaffari-Nasab et al. [19] represented the uncertain parameters by the midpoint of the two interval-valued numbers to solve the capacitated single and multiple allocation hub location problems. Habibzadeh Boukani et al. [20] transformed a kind of two interval-valued parameter model into an exact real number type-based one using the expectations of the interval cost and capacity. Zetina et al. [21] developed an exact real number type-based model to represent the uncertain parameters with the interval width values of the two interval-valued parameters. Talbi and Todosijevic [22] quantified the solution uncertainty by the widths of two interval parameters. Similarly, some studies generated the interval values of parameters such as the travel time to overcome the shortcomings of the p-hub model. Amin-Naseri et al. [23] established a robust bi-objective uncapacitated single allocation p-hub median problem by setting the minimization of travel time uncertainty as the objective function. De sá [24] introduced uncertain budget factors and set demand and travel time as interval numbers to overcome the defects of exact models. These above researches on interval numbers successfully leverage the continuous information of parameters in the p-hub median problem, rather than relying solely on a finite set of discrete scenarios.
To raise the accuracy of the decision results, decision-makers should provide not only the upper and lower bounds but also the most probable values of related parameters while selecting the location of hubs. Triangular fuzzy numbers can more realistically represent the parameters of hub location problems. Therefore, many Scholars construct a triangular fuzzy number and convert it into a precise model for a solution. Maharjan [25] transformed the temporary logistics hub location-allocation model with triangular fuzzy numbers into an exact real number type-based one based on a fuzzy chance constraint of the credibility distribution theory. Kayvanfar [26] proposed a method to solve the multi-objective and multi-echelon supply-distribution model based on the fuzzy triangular reference set weight. TootoonI et al. [27] developed triangular fuzzy type I and II programming approaches with interval models to solve the single allocation ordered median hub location problem. These above researches use triangular fuzzy numbers to solve the of hub-location allocation problem. Furthermore, in terms of path allocation, Tirkolaee et al. [28] described the demand as a triangular fuzzy number and established a fuzzy chance-constrained programming model to solve the fuzzy multi-trip location-routing problem. Sun [29] employed a defuzzification, linearization, and weighted sum method to present a crisp linear model to solve the routing problem in the road-rail intermodal transportation network. Compared with interval numbers, these methods with triangular fuzzy numbers fully utilize the information of the parameters in the p-hub median problem.
Through a comprehensive review of studies on the p-hub median problem under uncertain conditions, the following conclusions can be drawn: (1) Current research often transforms interval models into various precise real-number models for solution purposes. (2) Existing algorithms have not directly addressed mathematical programming models with triangular fuzzy number parameters. These research frameworks only leverage partial characteristics of triangular fuzzy number parameters, rather than considering their complete information, which may compromise solution accuracy. To further enhance the quality and precision of decision solutions, this study proposes a novel optimization approach building upon the models and algorithms employed by predecessors in addressing the p-hub median problem. This approach integrates fuzzy mathematics theory with heuristic algorithms, incorporating the membership function evaluation index of triangular fuzzy numbers into heuristic algorithms. Consequently, it establishes an algorithm capable of fully utilizing information from uncertain parameters and directly solving the p-hub median problem with triangular fuzzy numbers. In comparison to prior work, the contributions of the present study can be summarized as follows.
(1) A fuzzy NSUMApHMP model based on the triangular fuzzy number is proposed, capable of determining the optimal hub network under uncertain conditions and revealing paths with the lowest cost at different scales of the p-hub median problem.
(2) A customized heuristic algorithm is devised with a triangular fuzzy number evaluation index as the fitness operator, allowing for the direct solution of the triangular fuzzy number p-hub median problem. The evaluation index for triangular fuzzy numbers is based on the relative positional relationships derived from membership function deductions, employed to assess the cost-effectiveness of hub location schemes. This customized heuristic algorithm is an improved Genetic-Tabu Search algorithm utilizing an enhanced Floyd-Warshall algorithm to explore the shortest paths for all OD pairs in the hub network.
The remaining sections of this paper are organized as follows. Section 2 introduces an NSUMApHMP model based on triangular fuzzy numbers. A fuzzy evaluation index is constructed as a fitness operator to assess the merits of various fuzzy hub schemes. Section 3 devises a Genetic-Tabu Search algorithm, considering the proposed evaluation index as a fitness operator. In Section 4, the effectiveness of the proposed algorithm and optimization approach is validated based on the classical Civil Aviation Bureau (CAB) dataset and several self-generated datasets. The sensitivity analysis verifies the impact of triangular fuzzy number type characteristics on the hub scheme. Finally, Section 5 concludes the entire paper.
2. Mathematical formulation
This research aims to generate high-quality decision-making solutions to the p-hub media problem by fully utilizing the information in triangular fuzzy number parameters. Therefore, a unique optimization approach based on a domain search heuristic algorithm is proposed. During the domain search process, the triangular fuzzy number evaluation index is used as a fitness operator to sort the cost of hub solutions. This approach can solve the interval model of the p-hub medium problem directly, avoiding the information loss caused by converting the interval model into an accurate real number model. This section introduces the research problem and mathematical model first, and then the establishment process of the evaluation index for the hub scheme of triangular fuzzy number type as a fitness operator in the algorithm.
2.1. Model description and assumptions
The advantage of Non-Strict Uncapacitated Multi Allocation p-hub Median Problem (NSUMApHMP) in p-hub median problem is that it allows for the allocation of multiple hubs, which results in better adaption to complex transportation networks and demand distribution in reality [30]. Based on this model, this study proposed a novel fuzzy NSUMApHMP model considering passenger flow and cost as uncertain parameters. In this model, triangular fuzzy numbers are used to represent uncertain passenger flow and cost, in order to better fit the real air transportation situation. By incorporating triangular fuzzy numbers into the model, the utilization of continuous parameters in the optimization-solving process is ensured, thereby avoiding the information loss caused by parameter discretization. The aim of the fuzzy NSUMApHMP model is to determine the locations of p hub cities, the allocation of non-hub cities to hub cities, and the transfer flow ratio of each OD pair (i,j) via hubs k and m, with the objective of minimizing the total system cost. Following classic assumptions listed in Campbell [3] and O’Kelly [2], the assumptions of our paper can be depicted as follows.
(1) Assume that all hubs are fully connected to each other. Each non-hub node can be connected to any multiple hub nodes, and the non-hub nodes can also be directly connected to each other.
(2) Suppose that any one OD passenger flow can be transferred between two hub nodes at most.
(3) The hub’s capacity and the hub network’s links are infinite.
(4) The number of required hubs to locate is fixed in advance.
(5) Assume that there are economies of scale between hub nodes such that the discount factor can be applied to the unit flow cost of routing between a pair of hubs.
(6) The design parameters including the unit flow cost and demand are α (0≤α≤1) considered to be the triangular fuzzy number types.
2.2. Notations
(1) D: the total number of nodes in a network.
(2) p: the number of selected hubs.
(3) : the passenger demand volume between node i and j, in the form of triangular fuzzy number type, generally has
, but for simplicity, it assumes
.
(4) : the unit flow cost from node i directly to node j, in the form of triangular fuzzy number type.
(5) : the unit flow cost of transporting a passenger from node i to node j via hub k and m, where
.
(6) Xijkm: the proportion of an OD passenger demand volume transported via path i→k→m→j (where k, m are hubs) to the total OD passenger demand volume.
(7) Zij: the proportion of an OD passenger demand volume directly transported to the total OD passenger demand volume.
(8) Yk: 1 if node k is selected as the hub, 0 otherwise.
2.3. Mathematical model
The triangular fuzzy number type-based NSUMApHMP model based on Campbell’s 4-index model [3] can be written as follows.
(1)
s.t.
(2)
(3)
(4)
(5)
(6)
(7)
(8)
The objective function is to minimize the sum of the transportation cost of the hub scheme of the triangular fuzzy number type. Constraint (2) ensures that the total number of hubs selected equals p. Constraints (3) ensure that each OD passenger demand volume can be totally transported. Constraints (4)-(5) ensure that each non-hub node can be connected to multiple hub nodes. Constraints (6)-(8) indicate the type and the range of values of decision variables.
The objective function of this model contains parameters of the triangular fuzzy number type, and it is apparent that the proposed model preserves the full information of these uncertain parameters. However, this makes the solution challenging as the objective function of the proposed model is the total transportation cost of triangular fuzzy number type, which cannot be directly solved using linear programming methods. To assess the superiority of the fuzzy hub scheme and obtain the global optimal solution, an evaluation index and a customized algorithm are required.
2.4. Triangular fuzzy number evaluation index based on membership function
To design a heuristic algorithm based on triangular fuzzy numbers, it is necessary to compare the fitness of each hub scheme of the triangular fuzzy number type during the algorithm iteration process. The following subsection provides an evaluation index that can be used to compare the cost of hub schemes in any scenario and rank multiple hub schemes. By comparing the relative positions of membership functions for transport cost based on triangular fuzzy numbers, a comparison of hub schemes is conducted.
2.4.1. Basic definition.
Definition 1 [31] Let be said as a triangular fuzzy number where
. The lower bound of
is defined as AL, the upper bound of
is represented as AN and the most probable value of
is expressed as AM. Characteristic membership function of
can be expressed as Eq (9).
(9)
Note that when
, the triangular fuzzy number
is indicated as an ordinary positive real number. Several definitions related to operators of triangular fuzzy numbers are shown as follows.
(1) Addition: .
(2) Multiplication: ,
(λ>0).
Definition 2 [32] Let X be a domain, and the set of all fuzzy sets on X is defined as F(X). If and
, the definition of
can be written as Eq (10).
(10)
It calls
as the α-cut of the fuzzy set
. Similarly, for ∀α∈[0,1],
.
2.4.2. Comparison of hub schemes of triangular fuzzy number type.
The positions of numbers besides two hub schemes of the triangular fuzzy number type ( and
) can express the position relationship between them in the coordinate system. As shown in Fig 1, if there are no overlapping parts of two hub schemes of the triangular fuzzy number type, it’s easy to be denoted by Scenario 1 that
, all parts of
are to the right of
. At this point,
is strictly superior to
. As shown in Fig 2, in Scenario 2, if BN<AL, all parts of
are to the left of
. Similarly,
is strictly superior to
.
If two hub schemes of the triangular fuzzy number type overlap, all eight possible scenarios can be depicted from Figs 3–10. These eight specific scenarios can be expressed using the following formulas.
Scenario 3: . Scenario 4:
.
Scenario 5: . Scenario 6:
.
Scenario 7: . Scenario 8:
.
Scenario 9: . Scenario 10:
.
As shown above, all overlapping relationships between two hub schemes of the triangular fuzzy number type can be divided into at most six regions, namely ①, ②, ③, ④, ⑤, and ⑥.
2.4.3. Construction of the evaluation index based on position relationship of membership function.
Based on the overlapping area of membership function, the fuzzy preference relations can be expressed as
(11)
Geometrically, S1 shows the area of the left part of
superior to the left part of
. S2 shows the area of the left part of
superior to the left part of
. S3 shows the area of the right part of
superior to the right part of
. S4 shows the area of the right part of
superior to the right part of
. S5 shows the area of the overlapping part of
and
.
For direct calculation, Eq (11) can be transformed into several mathematical expressions containing AL, AM, and AN for through Ri ({Ri|i∈1,2,…6} corresponds to ①, ②, ③, ④, ⑤, and ⑥ in Figs 1–10, respectively), as shown in Table 1.
It can be further summarized as
(12)
and
(13)
Note that there’s no area overlapping in Scenario 1 and Scenario 2 as shown in Figs 1 and 2, respectively. The corresponding values of
are obtained by Eq (11) equaling 0 and 1 respectively. We find the superiority of
and
in Scenario 1 is only related to area ③ as well as in Scenario 2. Therefore, two new formulas can be added to represent Scenario 1 and Scenario 2 (as shown in Table 2).
Then a unique numerical representation of the position relation can be formed and is applicable to all scenarios. When is reduced, the position of
is more inclined to the left of
. Therefore, the inferiority of
to
is increased. Conversely, the position of
is more inclined to the right of
when
is increased. Similarly, the superiority of
to
is increased. The range of values of the
for each scenario tested is shown as Eq (14).
(14)
Comparing the quality of multiple hub schemes of triangular fuzzy number type requires a reference point. The fuzzy mean is proposed as a reference point
, which is calculated using
. A larger
value indicates that the fuzzy number is closer to the right than the fuzzy mean, making it better than the other fuzzy numbers.
Definition 3 Let the two hub schemes of the triangular fuzzy number type be and
respectively. Then the fuzzy average value can be defined as Eq (15).
(15)
Where
is defined as the fuzzy average of
and
. Then
can be used to evaluate
and
as follow.
(1) Calculate and
respectively.
(2) Compare the magnitude of both values.
1) If , then
.
2) If , then
.
3) If , then
.
Using the following steps, the index can evaluate the quality of the multiple triangular fuzzy numbers.
Step 1: Let .
Step 2: Find the .
Step 3: Calculate .
Step 4: Evaluate based on the magnitude of indicator values.
3. Triangular fuzzy number-oriented Genetic-Tabu Search algorithm
Section 2 finished the establishment of the model and the design of the fitness operator in the algorithm structure. In this section, to obtain high-quality hub solutions in an efficient solution time, we further design a complete heuristic algorithm based on the proposed fuzzy NSUMApHMP. A Genetic-Tabu Search algorithm is proposed. It generates superior initial domain space using the Genetic Algorithm, and then the optimal hub scheme is obtained through iterative optimization by the global search capability of the Tabu Search algorithm.
3.1. Generating initial solution
The Tabu Search (TS) algorithm is well-known for its local search capabilities and memory functions. The quality of the initial solution can significantly affect its performance, particularly in large-scale problems. The convergence speed of TS algorithm is slow when the initial solution is sub-optimal, which limits its practical applicability to real-world problems.
To improve the efficiency of TS algorithm in solving large-scale problems, a hybrid strategy of Genetic-Tabu Search algorithm, which combines the advantages of Genetic Algorithm and Tabu Search algorithm, was proposed. It overcomes the shortcomings of TS in initial solutions of large-scale problems. The Genetic Algorithm is used to generate a high-quality initial solution, which are served as the input for the Tabu Search algorithm. The following pseudocode example illustrates the process of generating initial solutions.
Algorithm 1: Initialization based on Genetic algorithm
Input:
Scale (n), Number of hubs(p)
output:
Initial hub allocation
while reaching stopping conditions
/* Initialization of the population (Hubs)
for 1 to n
randomly generate a chromosome that includes the network structure (hub selection and connection)
end for
/* Chromosome crossover
while the offspring population < n
select parents using the roulette wheel selection method, swap hub locations, and then reassign connections based on the updated hub locations
/* Chromosome mutation
mutate one hub and reassign connection methods
end while
/* Calculate fitness value
for 1 to n
calculate the total network cost as the fitness of each chromosome
end for
end while
/* Return output (Hubs)
3.2. Defining the neighborhood space and its candidate solution set
Given that there are p(n−p) feasible solutions in the neighborhood space, all obtained neighborhood solutions are served as the candidate solution set to be generated for the next step.
3.3. Designing fitness operator
As mentioned above, the fitness operator is basis of the evaluation index. It can be depicted as the following pseudocode example. In addition, the fitness operator will be also used in the Floyd-Warshall algorithm to solve the shortest path problem for node connectivity.
Algorithm 2: Computation of the fitness operator
Input:
size of the candidate solution set (m),cost parameters of solution i in the candidate solution set , fuzzy average of the objective function values
,
output:
evaluating results of the solution (ranki)
/* Conduct evaluation
for i in m
if
if
if
if , then
else
else
else
if , then
else
else
if
if , then
else
else
if , then
else
if , then
else
end for
/* Return output ranki
3.4. Determining Tabu list and length
To minimize the roundabout search, the Tabu list follows the First-In-First-Out (FIFO) principle and uses the fixed setting method to determine the Tabu length. Each exchanged hub node is given a Tabu length and then is added to the Tabu list. It can be released only if its Tabu length equals zero or if it satisfies the aspiration criterion.
3.5. Setting aspiration criterion
The aspiration criterion adopted in this paper includes two aspects.
(1) If a Tabu candidate solution is better than the current optimal solution, it is updated with the current solution and the optimal solution.
(2) If all candidate solutions are Tabu and are inferior to the current optimal solution, then the optimal Tabu candidate solution is updated to the current solution.
3.6. Procedure of Genetic-Tabu Search algorithm
By inputting the passenger demand volume parameter and unit flow cost parameter of the triangular fuzzy number type into the Genetic-Tabu Search algorithm, the optimal hub solution is obtained through iteration of the algorithm process. The whole procedure of the Genetic-Tabu Search algorithm can be depicted in Fig 11, where Sbest is the optimal solution, l is the number of iterations, lMax is the maximum number of iterations, S0 is the initial solution, Sl is the current solution, is the best solution among the current candidate solutions, and CountBest is the maximum number of times to keep the optimal solution.
As shown in Fig 11, the algorithm begins by generating an initial solution using a Genetic Algorithm. In each iteration, for each non-hub node and each hub node in the set, individual exchanges are performed. Among these exchanges, the algorithm selects the optimal candidate solution. If this candidate solution is superior to the current best solution, the solution updates. Otherwise, the solution is checked for Tabu status. If all candidate solutions are Tabu, the best Tabu solution is chosen as the current solution. After an iteration, the algorithm updates the current best solution, the current solution, and the Tabu list. Subsequently, it checks whether the stopping criteria are met. If they are, the results are output, and the algorithm terminates. Otherwise, the process continues with further exchanges and updates until the stopping criteria are satisfied.
4. Numerical example
To verify the effectiveness of the customized algorithm and the proposed approach in an uncertain environment in this research, we conducted numerical example calculation experiments. In this section, the Civil Aeronautics Board data set (CAB) introduced by O’Kelly [2] and several self-constructed data sets are used to verify the proposed optimization approach and customized algorithm. The experiment is implemented under a computer environment with Windows 10, 16G of computing memory, and 2.6GHz CPU, using Java JDK version 1.8. The Tabu length, the stopping condition, and the maximum number of iterations are set to 10, 50, and 1000, respectively. To test the algorithm’s effectiveness in solving large-scale problems, the CAB is used to randomly generate a network with 200 nodes, and the original 25 nodes with real data (called for self-constructed data set for simplicity) were selected to verify the reliability of the model.
4.1. Data generation and its illustration
Self-constructed data set and triangular fuzzy number-based parameters for unit flow cost and demand volume are generated using the CAB data set. The dataset’s real number data are set as the most probable values AM in numerical experiments. A method is developed to generate all related parameter values in the proposed model, ensuring randomness and representing all overlapping scenarios. Additionally, a method is developed to generate probability distributions conforming to triangular fuzzy numbers for stochastic simulation experiments.
(1) Generation of self-constructed data set
To evaluate the performance of the algorithm on a large scale, the 25 city nodes from the CAB dataset were taken as the basic nodes and the CAB dataset expanded to a self-constructed data set with 200 nodes. Subsequently, 175 unique virtual nodes were randomly generated within a 2000×2000 grid to ensure that the triangular inequality holds between all nodes. The unit flow cost of the generated virtual nodes is proportional to their distance from other nodes, while the corresponding passenger flow is randomly generated within the range of the minimum to maximum demand volume from the original CAB dataset.
(2) Generation of triangular fuzzy number parameters
The demand volume and unit flow cost
are assumed to yield normal distribution and are independent of each other. The data in the CAB data set are taken as the mean μ and the variance is set to βμ(β<1). To investigate the variation of the hub solution when the demand volume and unit flow cost are changed, data sets with narrow intervals and wide intervals are designed. The specific randomly generated data format is shown in Table 3.
(3) Data generation method for simulation experiments
The real number data for the simulation experiments should to satisfy the triangular fuzzy number probability distribution. Here, the probability density function of the triangular fuzzy number by 0.5(AL,AN) is calculated with the membership function in Eq (9) in Section 2.
(16)
Transforming Eq (16) into a probability density function, and then we use the inverse transformation method to obtain the stochastic simulation formula Eq (17) for variable x. r is a random number that yields [0, 1] uniform distribution.
4.2. The effectiveness of the algorithm
To evaluate the effectiveness of the proposed customized algorithm in solving the NSUMApHMP, the CPU time and accuracy of the proposed algorithm are compared with other algorithms on the self-constructed data set, including the Genetic Algorithm (GA) from Qin Z [33], the Tabu Search algorithm (TS), and the Branch-and-Bound algorithm (B&B) in the GUROBI solver. GTS, TS, and GA are compared with symmetric triangular fuzzy numbers and real numbers generated from the self-constructed data set. Since the GUROBI solver cannot directly utilize fuzzy numbers, it can only perform calculations using the real numbers generated from the self-constructed data set. The experimental instances were set with a range of 10 to 200 nodes, encompassing 11 different scales, and a discount factor α of 0.8 was applied. In addition, for the parameter settings of GA, the experiments are set with a population size of 400, a crossover probability of 0.2, probabilities of 0.05 and 0.3 for the hub and linkage condition group variants respectively, a discount factor α = 0.8, and the stopping conditions equal to the values in GTS.
The accuracy of all algorithms was compared using GAP as shown in Eq (18).
(18)
For GTS, TS, and GA in real number data, their UBs are the upper bounds obtained by the GUROBI solver. Their LBs is the optimal objective function value. For GTS, TS, and GA in triangular fuzzy number data, their UBs are the upper bounds obtained by the GTS. The results for real number data and triangular fuzzy number data are displayed in Tables 4 and 5, respectively.
Based on the performance of various algorithms using the real-number data in Table 4, it can be observed that the B&B algorithm can only obtain effective solutions for up to 40 nodes within 3 hours. From the perspective of CPU time, for the same problem size, GA requires the longest CPU time, followed by TS, while the average CPU time of GTS is 49.05% and 40.93% shorter than that of GA and TS, respectively. The improvement of the GTS algorithm lies in its enhanced utilization of the advantages of the Genetic Algorithm during the search process. In terms of the solution accuracy measure, TS can maintain the lowest GAP for small-scale instances with nodes less than 80, attributed to its characteristics of neighborhood search. With the increase in scale, the average GAP for TS and GA became 2.79% and 2.77%, respectively. In comparison, GTS, utilizing GA to generate high-quality initial solutions and using them as the initial neighborhood for search, converges to superior solutions during problem-solving. It results in lower GAP of 2.09% in large-scale instances. This insight suggests that, in practical applications, particularly for problems involving a larger number of nodes, the integration of various algorithmic strengths can be employed to enhance computational efficiency.
Table 5 presents the results of calculations using triangular fuzzy numbers, and it is evident that B&B algorithm cannot solve them due to its nature as an exact algorithm. By comparing the average CPU time of the other three algorithms in Table 4, the average CPU time in Table 5 of GTS, TS, and GA increased by 11.74%, 14.32%, and 15.35% respectively. This indicates that although introducing a triangular fuzzy number evaluation index can better express uncertainty, it also comes with a certain computational cost. Nevertheless, as shown in Table 5, GTS continues to maintain faster computation time and higher computational accuracy, with a CPU time of 872.85s and the lowest GAP.
4.3. The effectiveness of the optimization approach under uncertainty
To validate the effectiveness of the proposed optimization approach (Approach 1) in an uncertain environment, we compared it with three other classical approaches: a stochastic optimization approach (Approach 2) [14], a robust optimization approach (Approach 3) [34], and an approach based on real-point data types [4] (Approach 4). First, we generated passenger demand and unit flow cost parameters based on the triangular fuzzy number data in Table 3. Next, the parameters were used to solve the hub schemes for all four approaches. Approach 1 employed a customized algorithm, while the other approaches used the GUROBI solver. Subsequently, as described in Section 5.1 (Method (2)), the generated real-type parameters for passenger demand and unit flow cost were applied in a simulation experiment. The specific steps of the simulation experiment involved randomly generated real values for passenger demand volume and the unit flow cost, which are within the specified range of the generated triangular fuzzy number parameters. These real parameters were then substituted into the hub schemes of the four approaches to calculate the corresponding total cost. Through simulation experiments, the adaptability of the four approaches in uncertain environments can be compared.
In this experimental section, the other parameters were set as the narrow interval and the wide interval condition: the scale was set to 10, 15, 20, and 25, the number of hubs was set to 3, 4, and 5, and the discount factor was set to 0.4 and 0.8. The obtained hub schemes by the four approaches under different instances with given uncertain parameters are presented in Tables 6 and 7.
(1) Under the same conditions, the hub solutions of the four approaches vary with different scales. As the scale increases, the differences in hub solutions among the four methods become more significant. This implies that in real-world aviation operations, as the scale of the problem expands and demand within the network increases, it leads to different layouts for real aviation hubs.
(2) When the discount factor changes, the hub solutions of the four approaches change. This indicates that the discount factor significantly influences the decision outcomes. The reason is that changes in the discount factor led to different trade-offs in transportation costs, affecting the selection of real aviation hub locations and the distribution of network traffic.
(3) When the width of the interval changes, the hub solutions of the first three approaches, which consider uncertainty, change. Since Approach 4 does not consider uncertainty, the hub solutions remain unchanged. This suggests that approaches considering uncertainty can more flexibly adapt to changes in the real aviation operational environment.
To further compare the adaptability of these approaches in uncertain environments, Figs 12–15 display the average cost of their generated solutions over 1000 simulation experiments.
The results in Figs 12–14 indicate that Approach 4, which does not consider any uncertainty, exhibits the highest average total cost in each scenario, indicating its inferior adaptability in uncertain environments. Approaches 2 and 3, employing random optimization and robust optimization, respectively, show a reduced average total cost compared with Approach 4. In most cases, Approach 1, incorporating the triangular fuzzy number evaluation index, provides the lowest average total cost solution. This is because Approaches 2 and 3 can only utilize information limited to several discrete and finite scenarios. In contrast, Approach 1, employed the proposed evaluation index and utilized continuous uncertain parameters represented by triangular fuzzy numbers.
Furthermore, variations in the discount factor also impact the hub selection of the four approaches. A detailed analysis of the fluctuations in hub schemes (changes in hub nodes) under different discount factor conditions for the four approaches is presented in Table 8, comparing the specific changes in hub schemes between Tables 6 and 7. The hub configuration of Approaches 1 to 4 changed 48,39,40 and 8 times, respectively. This indicates that the evaluation index of our approach in this study can more accurately balance changes in transportation cost within the aviation hub network. It also suggests that our approach in this study is better suited to cope with external environmental fluctuations.
The above results imply that, in addressing the p-hub median problem in the aviation transportation industry, our optimization approach, which considers continuous uncertainty factors, can offer superior and robust solutions. It provides a more comprehensive modeling of uncertainty, enabling better adaptation to changes in uncertainty levels and discount factors. Consequently, it aids in optimizing the location of aviation hubs and the allocation of traffic, enhancing operational efficiency, and reducing cost.
4.4. Sensitivity analysis
To further investigate how the uncertainty level of triangular fuzzy number parameters affects the total transportation cost of hub solutions, this section examines the proposed optimization approach in obtaining solutions for NSUMApHMP in an uncertain environment. The uncertainty levels of passenger demand volume and unit flow cost are reflected by variations in the interval width and the most probable values, with a discount factor of 0.8. The specific results of the hub scheme are provided in Tables A1 to A5 in S1 Appendix.
4.4.1. Impact of interval width.
(1) When the most probable value remains unchanged, increase the interval width
As shown in Fig 16, the lower bound of the parameters gradually moves towards zero, while the upper bound increases up to 1.2 times the original upper bound, with a step size of one-tenth of 0.2AR. The lower bound step size is one-eleventh of AL to ensure the effectiveness of the Floyd-Warshall algorithm. The upper and lower bound of the interval are simultaneously expanded outward by 10 steps. Hub scheme differences between each step and the original scheme are compared to obtain the hub fluctuation.
Fig 16 illustrates that, with the increase in interval width, the disparity between the initial hub configuration and the outcome in NSUMApHMP grows, indicating an escalation in problem uncertainty. This phenomenon is attributed to alterations in the relative positioning of the membership functions of triangular fuzzy number parameters, leading to variations in hub schemes. Additionally, as the scale increases, the extent of change in the hub schemes for a network with 25 nodes (depicted by the deep segment in Fig 16) is noticeably greater than that of networks with 10, 15, and 20 nodes (depicted by the light segments in Fig 16). This indicates a positive correlation between the impact of interval width on the hub scheme and the size of the network. In computations involving large-scale p-hub median problems, changes in interval width are more likely to impact hub scheme choices.
(2) When the most probable value remains unchanged, narrows the interval width
To narrow the interval, the lower bound is adjusted to approach and the upper bound is adjusted to approach
, while keeping both bounds from the most probable value. The step size for the upper bound is set to one-tenth of
and the step size for the lower bound is set to one-eleventh of
. The results of the variation in hub schemes are depicted in Fig 17.
According to Fig 17, it can be observed that when the interval width narrows, the hub scheme remains unchanged. This result is attributed to the fact that when the width of the parameters approaches the most probable value, the triangular fuzzy number parameters gradually transform into real numbers, leading to the reduction in the degree of uncertainty. This implies that as the range of critical parameters becomes clearer and more focused, the selection of hub schemes becomes more stable, no longer influenced by the broad range of uncertainty.
4.4.2. Impact of the most probable value.
The lower and upper bounds for the passenger flow and cost parameters remain fixed, while the most probable value within the corresponding interval is allowed to change. To achieve this, we designed one-tenth interval width steps that move from the lower to the upper bound. The most probable value moves from the lower limit to the upper limit in increments of the specified step size. The results of the variation in hub schemes are depicted in Fig 18.
Fig 18 demonstrates that variations in the most probable values significantly impact hub schemes across the majority of scales. Sensitivity to changes in the most probable values is higher across various scales compared with the sensitivity to interval width depicted in Fig 16. This implies that, for solutions to the p-hub median problem, accurate estimation of the most probable values is crucial for the rational planning and operation of hubs.
4.4.3. Impact of both interval width and the most probable value.
(1) Change the most probable value as the interval widens.
Similar to the approach in Section 4.4.1, we gradually increase the interval width here. However, in this case, the most probable value is set to move within a fixed range as the interval widens. The results of the variation in hub schemes are depicted in Fig 19. In comparison to Fig 16, the hub schemes obtained in Fig 19 exhibit greater variability, especially at the scale of 25 nodes. This increased variability is attributed to the combined intensification of uncertainty in the aviation hub network represented by the most probable value and interval width. This indicates that in practical applications of addressing the p-hub median problem, as the uncertainty in the most probable value and interval width increases, hub schemes become more flexible and dynamic, which means greater adaptability and robustness.
(2) Change the most probable value as the interval narrows.
To ensure that the most probable value stays within the bound of the new interval during movement, it is set to transition from to
. The range of variation for the upper bound is
, and the range of variation for the lower bound is
. The resulting hub schemes generated during parameter variation are documented in Fig 20. By comparing Figs 17 and 20, we observe significant fluctuations in hub schemes when both the interval narrows and the most probable value changes simultaneously. This indicates that as data uncertainty decreases, triangular fuzzy numbers transform into real numbers, converting the fuzzy problem into a deterministic problem represented by the most probable value. Therefore, the variation in hub schemes at this point is primarily influenced by the most probable value.
5. Conclusions
In this study, we proposed a linear programming model based on triangular fuzzy number parameters to investigate the NSUMApHMP under uncertain conditions. To address the challenges of directly solving linear programming models with triangular fuzzy number parameters, we customized a heuristic algorithm that employs a triangular fuzzy number evaluation index as the fitness operator. To accurately evaluate the hub cost of triangular fuzzy numbers during the algorithmic iterations, we designed a ranking formula based on the membership function. To ensure both speed and precision in solving the problem, we developed a Genetic-Tabu Search algorithm, which utilizes the memorization principles of Tabu Search to find the optimal solution on top of the initial search space generated by the Genetic Algorithm. The computational results, obtained using both the classical CAB dataset and a self-constructed dataset, demonstrate the effectiveness of our proposed approach.
The simulation experiment results demonstrated that the proposed optimization approach consistently provides excellent solutions for NSUMApHMP under uncertain conditions. To validate the effectiveness of the customized algorithm, a detailed analysis was conducted. By comparing the customized algorithm with traditional TS and GA, we revealed that GTS has advantages in both solution quality and speed for resolving NSUMApHMP, as it reduced the runtime by 49.05% and 40.93%, respectively, and delivered high-quality solutions for large-scale instances. In terms of adaptability to uncertain environments, this optimization approach outperformed classical random optimization, robust optimization, and precise real-number optimization approaches, offering cost-effective and more flexible hub schemes. Sensitivity analysis indicated that both the interval width and the most probable value influence hub scheme selection. Specifically, the impact of the interval width on hub schemes was positively correlated with the scale size, while the influence of the most probable value was generally more significant across various scales than the interval width. Therefore, decision-makers should focus on anticipating and predicting changes in the most probable value and exercise effective control over the range of uncertainties, particularly as the scale gradually increases.
Despite the findings of this study on NSUMApHMP in uncertain environments, there are still some limitations that warrant further investigation. Firstly, our analysis primarily focused on hub schemes within the route network, without delving into the allocation issues of aircraft in the aviation system. However, in practical aviation operations, the unit flow cost is contingent upon the aircraft types and frequencies allocated to the routes. Simultaneously, the allocation of aircraft types and frequencies on routes is dependent on passenger demand along the origin-destination (OD) paths. This interdependency creates a mutual correlation between fleet allocation and route network design. Future research efforts could involve joint optimization of hub networks and fleet planning in uncertain environments to more comprehensively address real-world aviation transport problems. Secondly, this study relied solely on triangular fuzzy numbers to represent uncertainty. In future research, exploration of other symbolic data, such as gradient fuzzy numbers, could be considered to provide a more refined representation of uncertainty, thereby yielding a more practical decision plan.
Supporting information
S1 Appendix. The specific results of the hub scheme of Section 4.4.
https://doi.org/10.1371/journal.pone.0297295.s001
(DOCX)
References
- 1. Campbell J. F. Hub location and the p-hub median problem. Oper. Res. 1996, 44, 923–935.
- 2. O’kelly M. E. A quadratic integer program for the location of interacting hub facilities. Eur. J. Oper. Res. 1987, 32, 393–404.
- 3. Campbell J. F. Integer programming formulations of discrete hub location problems. Eur. J. Oper. Res. 1994, 72, 387–405.
- 4. Skorin-Kapov D.; Skorin-Kapov J.; O’Kelly M. Tight linear programming relaxations of uncapacitated p-hub median problems. Eur. J. Oper. Res. 1996, 94, 582–593.
- 5. Contreras I.; Cordeau J. F.; Laporte G. Stochastic uncapacitated hub location. Eur. J. Oper. Res. 2011, 212, 518–528.
- 6. Zhai H.; Liu Y. K.; Yang K. Modeling two-stage UHL problem with uncertain demands. Appl. Math. Model. 2016, 40, 3029–3048.
- 7. Ghaderi A.; Rahmaniani R. Meta-heuristic solution approaches for robust single allocation p-hub median problem with stochastic demands and travel times. Int. J. Adv. Manuf. Technol. 2016, 82, 1627–1647.
- 8. Peiró J.; Corberán Á.; Martí R.; Saldanha-da-Gama F. Heuristic solutions for a class of stochastic uncapacitated p-hub median problems. Transport. Sci. 2019, 53, 1126–1149.
- 9. Azizi N.; Vidyarthi N.; Chauhan S. S. Modelling and analysis of hub-and-spoke networks under stochastic demand and congestion. Ann. Oper. Res. 2018, 264, 1–40.
- 10. Shang X.; Yang K.; Jia B.; Gao Z. The stochastic multi-modal hub location problem with direct link strategy and multiple capacity levels for cargo delivery systems. Transportmetrica. A. 2021, 17, 380–410.
- 11. Alvarez Fernandez S.; Ferone D.; Juan A.; Tarchi D. A simheuristic algorithm for video streaming flows optimisation with QoS threshold modelled as a stochastic single-allocation p-hub median problem. J. Simul. 2022, 16, 480–493.
- 12. Ghaffarinasab N. An efficient matheuristic for the robust multiple allocation p-hub median problem under polyhedral demand uncertainty. Comput. Oper. Res. 2018, 97, 31–47.
- 13. Wang C.; Liu Y.; Zhang P.; Yang G. Two-Stage Distributionally Robust Optimization for a Two-Allocation p-Hub Median Problem. Journal of Uncertain Systems 2021, 14, 2150004.
- 14. Meraklı M.; Yaman H. Robust intermodal hub location under polyhedral demand uncertainty. Transport. Res. B-Meth. 2016, 86, 66–85.
- 15. Ghaffarinasab N.; Zare Andaryan A.; Ebadi Torkayesh A. Robust single allocation p-hub median problem under hose and hybrid demand uncertainties: models and algorithms. Int. J. Manag. Sci. Eng. 2020, 15, 184–195.
- 16. Ge W. A mean-robustness stochastic programming model for p-hub median problem. Int. J. Elec. Eng. Educ. 2019, 0020720919879371.
- 17. Ghaffarinasab N. Exact algorithms for the robust uncapacitated multiple allocation p-hub median problem. Optim. Lett. 2022, 16, 1745–1772.
- 18. Pourghader Chobar A., Adibi M. A., & Kazemi A. A novel multi-objective model for hub location problem considering dynamic demand and environmental issues. Journal of industrial engineering and management studies. 2021, 8, 1–31.
- 19. Ghaffari-Nasab N.; Ghazanfari M.; Teimoury E. Robust optimization approach to the design of hub-and-spoke networks. Int. J. Adv. Manuf. Technol. 2015, 76, 1091–1110.
- 20. Habibzadeh Boukani F.; Farhang Moghaddam B.; Pishvaee M. S. Robust optimization approach to capacitated single and multiple allocation hub location problems. Comput. Appl. Math. 2016, 35, 45–60.
- 21. Zetina C. A.; Contreras I.; Cordeau J. F.; Nikbakhsh E. Robust uncapacitated hub location. Transport. Res. B-Meth. 2017, 106, 393–410.
- 22. Talbi E. G.; Todosijević R. The robust uncapacitated multiple allocation p-hub median problem. Comput. Ind. Eng. 2017, 110, 322–332.
- 23. Amin-Naseri M. R.; Yazdekhasti A.; Salmasnia A. Robust bi-objective optimization of uncapacitated single allocation p-hub median problem using a hybrid heuristic algorithm. Neural Comput. Appl. 2018, 29, 511–532.
- 24. de Sá E. M.; Morabito R.; de Camargo R. S. Benders decomposition applied to a robust multiple allocation incomplete hub location problem. Comput. Oper. Res. 2018, 89, 31–50.
- 25. Maharjan R.; Hanaoka S. A credibility-based multi-objective temporary logistics hub location-allocation model for relief supply and distribution under uncertainty. Socio-Econ. Plan. Sci. 2020, 70, 100727.
- 26. Kayvanfar V.; S. Sajadieh M.; Moattar Husseini S. M.; Karimi B. Analysis of a multi-echelon supply chain problem using revised multi-choice goal programming approach. Kybernetes 2018, 47, 118–141.
- 27. Tootooni B.; Sadegheih A.; Khademi Zare H.; Vahdatzad M. A. A novel type I and II fuzzy approach for solving single allocation ordered median hub location problem. Iranian Journal of Operations Research 2020, 11, 65–79.
- 28. Tirkolaee E. B.; Abbasian P.; Weber G. W. Sustainable fuzzy multi-trip location-routing problem for medical waste management during the COVID-19 outbreak. Sci. Total Environ. 2021, 756, 143607. pmid:33220997
- 29. Sun Y. Green and reliable freight routing problem in the road-rail intermodal transportation network with uncertain parameters: a fuzzy goal programming approach. J. Adv. Transport. 2020, 2020, 1–21.
- 30.
Zhou, Y., Yao, D., Sun, Q., & Wu, Q. Application of Genetic Algorithm in p-hub Airline Network Design Problem. In 2nd International Conference on Electronics, Network and Computer Engineering (ICENCE 2016). 2016, September, 298–303. Atlantis Press.
- 31. Van Laarhoven P. J.; Pedrycz W. A fuzzy extension of Saaty’s priority theory. Fuzzy Set Syst. 1983, 11, 229–241.
- 32. Dubois D.; Hüllermeier E.; Prade H. On the representation of fuzzy rules in terms of crisp rules. Inform. Sciences 2003, 151, 301–326.
- 33. Qin Z.; Gao Y. Uncapacitated p-hub location problem with fixed costs and uncertain flows. J. Intell. Manuf. 2017, 28, 705–716.
- 34. Bai M. G.; Jiang T.; Zhu J. F. The Robust optimization of the hub and spoke airline network design based on Tabu search. Chinese journal of Mathematics Asthmatics Practice and Theory 2008, 38, 60–69.