The optimal design of the supply chain network plays a vital role in reducing transportation costs and improving system performance (Nozari et al., 2022). Transportation of goods in the supply chain is an important activity that has a great impact on the overall performance of manufacturing companies and the final cost of products by sending and delivering raw materials and finished products on time. For this reason, companies are constantly under pressure to reduce costs as well as increase their service levels (Szmelter-Jarosz et al., 2021). Among other things that companies should continuously monitor is the collection of products distributed by their company. Therefore, in addition to meeting customer demand, companies are also responsible for collecting products. This issue has been announced by governments due to the harmful environmental effects. In order to overcome these pressures, they must regularly redesign their service network and improve logistics processes with better coordination of decisions at strategic and operational levels (Ghahremani-Nahr et al., 2022). Hub location problems have been widely used in recent decades due to their high applications in modern transportation and communication systems. Hubs in the transportation networks of goods and passengers, etc., are places that receive the flow from different sources and direct it to the destination nodes or other hubs; Therefore, the goal of hub location-routing problems is to meet demands by moving people, goods, and information through special facilities called "hubs" that are known as distribution intermediaries between source nodes and destination nodes (Alumur et al., 2021). In a location-routing network, due to the large volume of goods transported from origin to destination, effective design of a network can ensure high-quality services (Li & Wang, 2023). In this network, by reducing the number of links in the source and destination nodes, hubs prevent any direct connection of all nodes, which in turn leads to a reduction in costs, a reduction in $C{O}_2$ gas emissions, and a reduction in the time of goods transfer compared to networks with complete graph connections. (Ghaffarinasab et al., 2022). Also, using the economical factor (0 < α < 1) for communication between hubs also leads to a large reduction in the total network cost (Taherkhani & Alumur, 2019).

In recent years, various types of hub location problems have been studied, one of the most important of which is the hub covering location-routing problem. In these networks, each hub can only transfer goods and services to other nodes within a predefined radius and collect the returned goods. Hub Covering location-routing networks have attracted the attention of many researchers due to various application fields in long-distance communication networks, air transportation, goods transportation, postal services, public and urban transportation, rail transportation, and emergency services (Mokhtarzadeh et al., 2021). As the demand for freight increases with population growth, the importance of the proper design of transportation networks in uncertain conditions becomes apparent (Davari et al., 2013). The use of hubs in transportation networks not only significantly reduces operating costs and network structure preparation costs and delivery time, but also causes better services to customers and reduces the harmful effects of environmental pollution (Hasanzadeh et al., 2018). Therefore, the main concept of hub covering location-routing in the design of a transportation network is that a number of decisions must be made, including the location of hub facilities and non-hub nodes, the allocation of all source and destination nodes (non-hub) directly to the hub, and flow routing. These decisions are aimed at reducing the costs of the total network, reducing the amount of $C{O}_2$ gas emissions and reducing the maximum travel time of vehicles. In addition, considering whether there is a limit to transfer the amount of products through the transportation network or there is no limit, hub location problems are classified into two sub-problems with limited capacity and unlimited capacity. Capacity limitations may be related to network links or hub nodes (Rahmati et al., 2022). These decisions are aimed at reducing the costs of the total network, reducing the amount of $C{O}_2$ gas emissions and reducing the maximum travel time of vehicles. In addition, considering whether there is a limit to transfer the amount of products through the transportation network or there is no limit, hub location problems are classified into two sub-problems with limited capacity and unlimited capacity. Capacity limitations may be related to network links or hub nodes (Rahmati et al., 2022). There is also another type of capacity limitation, which is related to the number of non-hub nodes that can be assigned to a specific hub; Therefore, decision makers can include these considerations during their planning (Zahedi and Nahar, 2020).

According to the stated content, the aim of this research is to design the hub covering location-routing model with delivery and pick-up under uncertainty. To deal with the uncertainty parameters of the model, the two-stage stochastic programming (TSSP) method was used and different scenarios were proposed for each of the uncertainty parameters. Also, the combination of location and routing problems in the hub covering model has led to the use of LP-Metric, MOIWO and MOALO methods to solve the problem. The results of the multi-objective model will be analyzed in order to minimize the total costs, minimize the maximum vehicle travel time and minimize the amount of $C{O}_2$ gas emissions.

The structure of the paper is such that in the second part, the research literature is reviewed and the research gap is identified. In the third part, the modeling of the hub covering location-routing problem with delivery and pick-up under uncertainty is discussed. In the fourth part, various solution methods including LP-Metrics, MOIWO and MOALO methods have been investigated. In the fifth section, the analysis of the mathematical model and the sensitivity analysis of the problem are discussed, and in the sixth section, the conclusions and future suggestions of the research are discussed.

In recent years, hub location problems have attracted the attention of many researchers, and various models have been developed, such as covering hub problems, median hubs, P-hubs, etc. According to the focus of this paper on the hub covering location-routing problem, the research conducted has been reviewed in the following.

Lopes et al. (2016) solved the hub location-routing model using GA with random search. They also used the branch and cut algorithm to check the efficiency of their solution method. Computational experiments on several benchmark examples of routing problems showed that the proposed algorithm is effective on medium-to-large-sized examples, while the branch and cut algorithm solves small to medium-sized problems optimally. Karimi (2018) presented a model of a hub covering location-routing with simultaneous delivery and pickup. In this paper, the main goal was to minimize the costs of hub location and vehicle routing considering the limited hubs capacity. The results showed that the Tabu search method is more effective than other methods in finding good-quality solutions. Karimi and Setak (2018) developed a transportation scheduling and hub location problem where the hub network is not fully interconnected. The objective of the proposed model was to maximize the total amount of flow delivered by the last predetermined time and also to minimize the total costs. Wang and Qin (2020) presented a model for the hub covering location problem under uncertainty. The primary goal of this model was to minimize the costs of setting up hubs. To solve this model, a hybrid algorithm named GA-VNS was used and by analyzing several numerical examples, the high efficiency of this method was shown. Rajabi and Darestani (2020) presented a bi-objective optimization model to reduce the time and costs of the hub covering location problem. In this model, the cost, time and capacity of the central hub were considered as fuzzy parameters. NSGA II and NRGA were used to solve the bi-objective model. The results show that there is no significant difference between the CPU. Also, by ranking the solution methods with TOPSIS, the NSGA II has shown its efficiency in solving large numerical examples. Pourmohammadi et al. (2021) presented a multi-objective mathematical model for hub location and routing problem under uncertainty in flow, cost, time and number of job opportunities. The objective of this model is to minimize the total cost of transportation including routing and fixed cost and to maximize employment and regional development as a social responsibility. FIWO algorithm, variable neighborhood search and game theory were used to solve the model. Ghiasi and Vahdani (2021) presented a model for the hub location-routing problem considering balance and capacity decisions. To reduce the total cost of location and routing, GA and ICA have been used. The comparative results of the problem examples show that the GA has a better performance compared to the colonial competition algorithm and considering the capacity and balance considerations can be effective in reducing the costs of the investigated network. Ratli et al. (2022) designed a heuristics method to solve the hub location-routing problem. Their improved algorithm used a sequential variable neighborhood descent with multiple neighborhoods required to match the complex nature of the problem. The calculation results on 912 samples showed that the presented innovative method has the ability to solve 691 problems. Fallah- Tafti et al. (2022) developed a mathematical model for rapid transit networks based on the hub location model, which includes stops (stations) in hub and non-hub levels. In this model, which aimed to minimize the costs of locating the hub and routing the means of transportation, a local search algorithm was used. Hajipour et al. (2022) presented a mathematical model for locating the fire station and facilities assigned to the stations in different periods and emergency situations (wars and natural disasters). This model is designed dynamically according to the demand and covering radius of the facilities, based on traffic and the type of area and phase in different periods. Due to the NP-hardness of the problem, PSO and ABC were used to solve the problem, and the results showed that the provided PSO is better than ABC in terms of the quality of solutions and computing time. Ghaffarinasab and Kara (2022) modeled the multi-allocation hub location problem in which risk is measured using the conditional mean criterion. In this paper, the classic multiple allocation hub location problem, i.e. bladder hub, maximum hub covering and central hub, was investigated. Benders decomposition algorithm was used to solve the mathematical model. Rahmati et al. (2023) modeled a hub location problem with two-stage stochastic profit maximization with uncertain demand. In this model, the advanced sample mean approximation method was used to achieve the appropriate number of scenarios. Also, k-means clustering and self-organizing map clustering algorithms were used for clustering similar samples. Finally, the proposed models were analyzed using the famous Australia Post dataset. Motamedi et al. (2023) used six new techniques in prioritizing customers in the hub optimization problem. In these techniques, each strategy assigns a value to each customer based on distance and demand parameters, where customers are prioritized based on this value. Comparing the no-priority strategy with the best customer prioritization strategy shows that the no-priority strategy has an adverse effect on delivery time. Li et al. (2023) designed a multimodal transportation network based on hub location for emergency relief operations. The nonlinear model presented by them was a bi-objective model with minimization of transportation time and transportation costs. They used the MOGWO to solve this problem. The results showed that the presented model can provide practical management insights for departments and transport companies.

The literature review shows that there is a research gap in objective functions, constraints and solution methods. Since the minimization of network costs has been discussed in most of the papers. This paper tries to reduce the minimum vehicle traffic time by developing a mathematical model, in addition to reducing costs, and minimizing the amount of $C{O}_2$ gas emissions caused by vehicle traffic. On the other hand, there is no hub covering location model that integrates with the vehicle routing problem. Using the TSSP method can easily control the uncertainty parameters. In addition to delivering goods, the presented model also pick-up returned products at the same time. The literature review shows that the use of meta-heuristic methods is common, but new algorithms have not been used on the hub covering location-routing model. This problem made MOALO and MOIWO to be used in this paper in addition to the LP-Metrics method.

In this section, the hub covering location and routing problem with simultaneous delivery and pick-up uncer uncertainty is discussed. For example, Figure (1) can be considered as an example of the investigated problem.

The main goal in the investigated problem is to determine the optimal number and location of nodes from the network as hubs, assign non-hub nodes to hubs, determine the optimal transportation route to meet uncertainty simultaneous delivery and pickup, and determine the optimal number of vehicles to distribute products. Therefore, to explain the problem better, 7 nodes are considered according to Figure (1). It can be seen that nodes 3 and 6 are selected as hubs in this problem, nodes 1-2-5 are assigned to hub number 6 and nodes 7 are assigned to hub number 3. The way of assigning these nodes to the hubs is based on the covering radius considered for each hub. That is, to distribute the products to the nodes assigned to each hub, vehicles related to that hub are used. In this case, the customers who are within the transportation radius of each vehicle (selected hub) can be allocated to that hub. Also, another goal of this problem is the optimal routing of vehicles, which can be seen according to Figure (1), the optimal route of the 1st vehicle from hub number 6 includes going through nodes 1, 5 and 3 and the optimal route of the 2nd vehicle from hub number 3 as well. It includes going through nodes 7 and 4.

In this paper, parameters such as delivery and pick-up demand, transportation time and transportation cost are taken into consideration in an uncertainty manner, which take on different values in different scenarios, the result of which can be a change in the optimal location of the facility as well as the optimal routing of vehicles. Therefore, the modeling of a hub covering location and routing problem with delivery and pickup under uncertainty can be done according to the following assumptions and scope of the study:

Based on the above assumptions, the mathematical model seeks to minimize the total costs, minimize the maximum vehicles travel time and minimize the amount of $C{O}_2$ gas emissions caused by the traffic of vehicles at the same time. Considering the assumptions of the above model, the number of $i,j,k,m\in N$ nodes and $v\in V$ vehicles are considered. Other parameters and variables for decision making in uncertain conditions can be presented as follows.

Parameters

${\tilde{C}}_{i,j}$	Uncertainty transportation cost between node $i$ and node $j$
${\tilde{T}}_{i,j}$	Uncertainty transport time between node$i$ and node $j$
$D_{i,j}$	The distance between node $i$ and node $j$
${Co}_{i,j,v}$	The emission rate of ${Co}_2$ gas from the traffic of vehicle $v$ between node $i$ and node $j$
${Fk}_k$	The fixed cost of determining node k as a hub based on the quality of service
${FV}_v$	The initial fixed cost of using the vehicle $v$
${\tilde{W}}_{i,j}$	The amount of uncertainty delivery and withdrawal demand between node $i$ and node $j$
$P$	The number of nodes that should be selected as hubs.
$Q$	Number of vehicles in each hub
$R_v$	The maximum covering radius of the vehicle $v$ in the distribution of products based on the quality of service
${Capv}_v$	Vehicle capacity $v$ in product distribution
${Capk}_K$	Maximum hub capacity $k$
$\alpha$	Coefficient of economical factor
$O_i$	The total amount of output flow from node $i$ $O_i=\sum _j{\tilde{W}}_{i,j}\forall i$

Variables

$Tmax$	Maximum travel time by vehicles from hub centers
$F_{i,j,k}$	The amount of flow that starts from node $i$ and passes through hub $j$ and hub $k$
$Y_{i,j,k}$	The amount of flow that starts from node $i$and passes through node $j$ and hub $k$
$X_{i,j,v}$	1; if vehicle v passes through node $i$ and node $j$ 0; otherwise.
$H_k$	1; if hub $k$ is selected 0; otherwise
$Z_{i,v}$	1; if node $i$ is assigned to vehicle $v$ 0; otherwise.

Modeling the problem

$Min{OBF}_1={\displaystyle \sum _{i,j,k}}{\tilde{C}}_{j,k}{Y}_{i,j,k}+{\displaystyle \sum _{i,j,k}}\alpha {\tilde{C}}_{j,k}{F}_{i,j,k}+{\displaystyle \sum _k}{Fk}_k{H}_k+{\displaystyle \sum _{i,v}}{FV}_v{Z}_{iv}$	(1)
$Min{OBF}_2=Tmax$	(2)
$Min{OBF}_3={\displaystyle \sum _{i,j,v}}{CO}_{i,j,v}{X}_{i,j,v}$	(3)
$s.t.:$
${\displaystyle \sum _{\begin{array}{c}k\\ k\ne i\end{array}}}\left(Y_{i,i,k}+F_{i,i,k}\right)=O_i,\forall i$	(4)
${\displaystyle \sum _k}\left(Y_{i,k,j}+F_{i,k,j}\right)={\displaystyle \sum _k}\left(Y_{i,j,k}+F_{i,j,k}\right)+{\tilde{W}}_{i,j},\forall i,ji\ne j$	(5)
$Y_{i,k,j}\le O_i{\displaystyle \sum _v}{X}_{k,j,v},\forall i,j,kk\ne j$	(6)
${\displaystyle \sum _{\begin{array}{c}j\\ j\ne i\end{array}}}{X}_{i,j,v}-{\displaystyle \sum _{\begin{array}{c}j\\ j\ne i\end{array}}}{X}_{j,i,v}=0,\forall i,v$	(7)
${\displaystyle \sum _k}{H}_k=P$	(8)
${\displaystyle \sum _j}\left(F_{i,k,j}+F_{i,j,k}\right)=O_i{H}_k,\forall i,k$	(9)
$Z_{i,v}\le {\displaystyle \sum _{\begin{array}{c}j\\ j\ne i\end{array}}}{X}_{j,i,v},\forall i,v$	(10)
${\displaystyle \sum _v}{Z}_{i,v}=Q{H}_i,\forall i$	(11)
${\displaystyle \sum _i}{Z}_{i,v}=1,\forall v$	(12)
${\displaystyle \sum _{\begin{array}{c}j,v\\ j\ne i\end{array}}}{X}_{i,j,v}=1+\left(Q-1\right)H_i,\forall i$	(13)
$Tmax\ge {\displaystyle \sum _{i,j}}{\tilde{T}}_{i,j}{X}_{i,j,v},\forall v$	(14)
${\displaystyle \sum _{i,j}}{D}_{i,j}{X}_{i,j,v}\le R_v,\forall v$	(15)
${\displaystyle \sum _{i,j}}{\tilde{W}}_{i,j}{X}_{i,j,v}\le {Capv}_v,\forall v$	(16)
${\displaystyle \sum _{i,j}}{F}_{i,j,k}\le {Capk}_k{H}_k,\forall k$	(17)
$X_{i,j,v},H_k,Z_{i,v}\in \left\{\mathrm{0,1}\right\},\forall i,j,k,v$	(18)
$Tmax,F_{i,j,k},Y_{i,j,k}\ge 0,\forall i,j,k$	(19)

Equation (1) reduces the total costs. These costs include the cost of transporting products between nodes, the cost of locating nodes as hubs, and the fixed costs of using vehicles. Equation (2) shows the second objective function of the problem and includes minimizing the maximum travel time between nodes. Equation (3) deals with the minimization of CO2 gas emissions caused by vehicle traffic. Equation (4) calculates the total output current from each node. Equation (5) provides flow balance between nodes and ensures that the amount of pickup and delivery occurs simultaneously at each node. Equation (6) guarantees that the total output flow from each node can be carried by only one of the available vehicles. Equation (7) guarantees that each vehicle visits each node (customer) only once. Equation (11) guarantees that the vehicles must pass through the selected hubs. Equations (12) and (13) guarantee that each route should be assigned to only one vehicle. Equation (14) shows the maximum travel time between nodes. Equation (15) calculates the covering radius resulting from the traffic of each vehicle. The relation (16) guarantees that each vehicle cannot move more than its capacity. The relation (17) ensures that only the hubs that have the ability to exchange products between their nodes as much as less than their capacity are selected. The relation (18) and (19) shows the type of the decision variables.

Due to the uncertainty of the parameters of delivery and pick-up demand, transportation time and transportation cost, the TSSP method has been used to control the model. Therefore, by defining $s\in S$ as the number of scenarios, parameters and decision variables need to be redefined and modified as follows.

Sets

$S$	The number of scenarios $s\in S$

parameters

$C_{i,j,s}$	Transportation cost between node $i$ and node $j$ in scenario $s$
$T_{i,j,s}$	Transportation time between node $i$ and node$j$ in scenario $s$
$W_{i,j,s}$	The amount of delivery and pick-up demand between node $i$ and node $j$ in scenario $s$
$P_s$	The probability of occurrence of scenario $s$

Also, the amount of the total output flow from node $i$ can be defined as follows.

$O_{is}={\displaystyle \sum _j}{W}_{i,j,s},\forall i,s$

(20)

Decision making variables

$F_{i,j,k,s}$	The amount of flow that starts from node $i$ in scenario $s$ and passes through hub $j$ and hub $k$
$Y_{i,j,k,s}$	The amount of flow that starts from node $i$ and passes through node $j$ and hub $k$ in scenario $s$
$X_{i,j,v,s}$	1; if, in scenario $s$, vehicle $v$ passes through node $i$ and node $j$0; otherwise.
$Z_{i,v,s}$	1; if node $i$ is assigned to vehicle $v$ in scenario s0; otherwise.

Three different methods have been used to solve the mathematical model. These methods include LP-Metrics, MOALO and MOIWO. For this purpose, a brief description of the solution methods has been given below.

In this section, the solution of the multi-objective model of the hub covering location and routing problem with simultaneous delivery and pickup has been discussed. For this purpose and to show the output of the above model, a small numerical example has been designed and the multi-objective problem has been solved using the LP-Metrics method. In the following, in order to solve the problem in large scale, MOALO and MOIWO are used and statistical and non-statistical comparisons are shown on the results of this problem. Finally, the most efficient algorithm has been shown in terms of computational indicators and its performance compared to other algorithms.

In this paper, the modeling and solving of a new problem of covering hub location-routing with simultaneous delivery and pick-up under uncertainty was discussed. The non-determinism of the problem data led to the use of the TSSP method to control the uncertainty parameters. The designed model took different strategic decisions such as hub location and tactical decisions such as vehicle routing at the same time. Examining the results of solving the model with the LP-Metrics method showed that the costs related to the network have increased with the reduction of $C{O}_2$ gas emissions and the maximum time of vehicle traffic. Also, by examining the economical factor, it was observed that the more this value decreases, the more vehicles tend to move through the hubs and it leads to the reduction of the total network costs. This change in the route of vehicles leads to a reduction in the amount of $C{O}_2$ gas emissions, as well as a reduction in the amount of $C{O}_2$ gas emissions.

MOIWO and MOALO were used to analyze numerical examples in large scale and evaluate different solution methods. Investigations of these solution methods on different numerical examples also showed that there was a conflict between the objective functions, and with the decrease of the second and third objective functions, the first objective function increased. Also, the solution time of these solution methods has been much less compared to the LP-Metrics method. The LP-Metrics method has not been able to solve numerical examples of very large scale. Based on the results, it was observed that the LP-Metrics method was only able to solve sample problems number 1 to 3. While meta-heuristic algorithms with results close to the LP-Metrics method have solved all the numerical examples. In these analyses, the maximum percentage of relative difference between the results of meta-heuristic algorithms and the LP-Metrics method was less than 2.83%. Also, on average, MOALO has the highest NPF and the lowest SM and CPU-time. While MOIWO has won the most MSI. On the other hand, by examining the significant relationship between the averages of the indicators, it was also observed that there is a significant difference between the averages of the computing time.

The results obtained from this model help the managers in the design of the fleet network to manage the costs of the total network, the delivery time of the goods and reduce the harmful environmental effects. Also, on the other hand, the existence of the two-stage stochastic programming method also helps managers in analyzing different scenarios and their impact on strategic and tactical decisions. Due to the high efficiency of the solution methods proposed in this paper, managers can use MOIWO and MOALO to achieve the desired results in case studies with larger sizes. The limitations in this research, such as not considering the queuing system in the distribution and collection of products, can be suggested as one of the future researches. Also, the use of uncertainty control methods such as fuzzy robust can also increase the richness of the mathematical model.