Uncertain Multi-Objective Hazardous Materials Transport Route Planning Considering Resilience and Low–Carbon

Hazardous material transport accidents are events with a low probability and high consequence risk. With an increase in the proportion of hazardous materials transported on domestic roads, an increasing number of scholars have begun to study this field. In this study, a multi–objective hazardous materials transport route planning model considering road traffic resilience and low carbon, which considers the uncertainty of demand and time and is under the limit of the time window. It transports many types of hazardous materials from multiple suppliers to multiple retails with three goals (transportation cost, risk, and carbon emission). This model fills the gap in the research on hazardous material transportation in the field of low carbon, and this is the first time that road traffic resilience is considered in the transport of hazardous materials as one of the weight factors of risk calculation. We designed a improved ant colony optimization algorithm (ACO) to obtain the pareto optimal solution set. We compared the improved ACO with genetic algorithm and simulated annealing algorithm. The results show that the improved ACO has better solution quality and space, which verifies the validity and reliability of the improved ACO.


I. INTRODUCTION
The transport of hazardous materials is a type of special transport, that refers to the transport of unconventional materials by special organizations or technical personnel using special vehicles. Hazardous materials are inflammable and explosive materials, hazardous chemicals, radioactive materials, and other objects that can endanger the safety of human life and property. It is a great danger and may cause material loss or casualties if careless. According to the data of National Bureau of Statistics of China [2], the demand for hazardous materials shows a trend of continuous growth. The sales volume of hazardous materials transport vehicles from January to August 2020 increased by 9.7% year-onyear. Therefore, with the increasing demand for hazardous materials and the characteristics of low accident probability The associate editor coordinating the review of this manuscript and approving it for publication was Sotirios Goudos . and high consequence risk, academic circles have conducted in-depth studies on the transport of hazardous materials.
Erkut et al. [3] modeled the risk of hazardous material transportation and selected hazardous material transportation routes according to the preferences of decision makers. Kara et al. [4] built a model based on the two aspects of government and carrier, studied the interaction of transportation cost and risk, and designed a road network for hazardous material transportation. Zhou et al. [28] proposed a combined cut-and-solve and cutting plane method to minimize the transportation risk in hazardous material transportation. Pradhananga et al. [5] considered the transport of hazardous materials with time windows and used a meta-heuristic algorithm for multi-objective optimization to determine the selection path. Wang et al. [11] used the GA-FCM clustering method to study the influence of driver behavior on the transport of hazardous materials. Shen et al. [30] used the eXtreme Gradient Boosting (XGBoost) algorithm to analyze the factors influencing hazardous material transportation VOLUME 11, 2023 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ accidents. Parsafard et al. [26] proposed a shortest path model based on multi-criteria objectives to solve the path problem and determined the optimal path between the bulk terminal and gas station. Sun et al. [25] have designed hazardous materials network, considering the cost fairness between carriers. AlRukaibi et al. [7] conducted a risk analysis of the surrounding environment and crowd safety of the hazardous materials transport route, designed an event simulation scene, and inserted it into GPS for risk analysis. Deng et al. [31] used K-means clustering to classify hazardous material transportation accidents and proposed solutions corresponding to the classification results. Hu et al. [24] considered cost, risk, and emergency response capability and used a compromise weight model and evolutionary algorithm to solve the multi-objective hazardous material transportation problem. Beneventti et al. [13] considered maximizing the minimum weighted sum of distances between hazardous facilities and exposed vulnerable groups to solve the problem of locating hazardous facilities and routing selection in densely populated areas.
In fact, the interaction of road conditions, natural conditions, human factors, and other factors cause complicated transportation conditions, which may lead to uncertainties such as time, demand, and loss of input data, thus affecting the measurement and calculation of risks and costs, and having a great impact on route selection. Pamǔcar et al. [6] input data to solve the problem of hazardous material transport paths on urban roads. Sun [27] studied hazardous material roadrail multi-modal transportation plans in uncertain demand. Zhao et al. [8] regarded distance and population density as fuzzy variables, set accident probability and speed as random, and used fuzzy random theory for modeling and solving. Qu et al. [9] proposed a method for the dynamic optimization of the transport path of hazardous chemicals, minimizing the risk of finding the optimal path for the remaining tasks in the case of uncertainty. Hu et al. [10] set retailer demand as a fuzzy variable, established a credibilistic goal programming model considering vehicle loading changes, and compared the expected value with the given value to select a path. Moghaddam et al. [29] applied the game theory method to find the optimal path with the total distance and total risk as the target in the case of uncertain customer demand and service time. Ma et al. [12] solved the multi-objective problem of vehicle scheduling through robust discrete optimization in the case of uncertain risks. In addition, consider two aspects of different options for the government and enterprise. Ziaei et al. [14] considered various uncertainties such as accident probability and carbon emission factors, and utilized robust counterpart optimization to solve the multiobjective location-routing planning of multi-modal transportation. There have been in-depth studies on the transport of hazardous materials under uncertainty, but there are still some deficiencies in the consideration of low-carbon materials, and there are few studies related to the transport of low carbon hazardous materials.
In terms of traffic resilience, studies have mainly focused on the impact of road or bridge damage on urban road networks, the impact of congestion caused by road disturbances on traffic flow, and supply chain management and operation [15]. Traffic resilience mainly includes three aspects: resistance capacity, absorption capacity, and recovery capacity, which reflect the ability of the traffic system to maintain the original state when a disturbance occurs and the ability of the system to recover the original state after damage [16]. Once traffic resilience is destroyed, there will be a significant impact on the road network, and the risk of hazardous materials transportation will increase accordingly. However, accidents involving hazardous materials during transportation can damage the resilience of road traffic, and the interaction between the two can lead to more serious consequences.
None of the above studies combined resilience and lowcarbon factors in the transport of hazardous materials. Therefore, this study considers road traffic resilience as the weight factor of risk and considers carbon emission as one of the objectives of hazardous materials path planning research. Some contributions of this study are as follows: 1) Carbon emissions is taken as one of the objectives to reduce the impact on the atmosphere and environment in the selection of paths under uncertainty.
2) Road traffic resilience is considered in the selection of hazardous materials transport routes, which is taken as the weight factor of risk parameters to increase research on road traffic resilience in the field of hazardous materials transport.
3) Consider the transport route of many kinds of hazardous materials supplied by multiple suppliers.

A. UNCERTAINTY THEORY
To solve the uncertainty of time and demand in the transport of hazardous materials, this study introduces the uncertainty theory to deal with the above uncertain factors. Li and Liu [18] provide the following four axioms, where is a nonempty set representing the sample space, L is the power set of .
∈ L, is an element of the power set, representing events. M represents the uncertain measure (credibility measure) of an event.
Axiom 1 (Normality Axiom): The set function that satisfies these four axioms is an uncertain measure. The following are some definitions and theorems of the uncertainty theory: Definition 1: An uncertain variable ξ is defined on the uncertain measure space (M, L, ), (x) is an uncertainty distribution function, where x is a value of the set of real numbers R; then, Definition 2: ξ is an uncertain variable defined in the uncertain measure space (M, L, ); then, the expected value of the uncertain variable is Definition 3: (x) is an uncertainty distribution function, if (x) is continuous and strictly increasing, and ∀x ∈ R, is a regular uncertainty distribution. Its inverse function is the inverse uncertainty distribution −1 (α), where α is an uncertainty measure.

B. ROAD TRAFFIC RESILIENCE
Resilience mainly refers to the four aspects of system performance: robustness, redundancy, resources, and rapidity [16]. Road traffic resilience is the ability of a road traffic system to absorb and recover from disturbances [21]. It is expressed as the weighted sum of the resilience of all nodes, which is determined by the weighted average resilience of the roads [22]. To study the system's ability to maintain normal operation under disturbance and to recover to its original level after it is damaged, resilience has been used to evaluate the performance of the traffic system under disturbance [21]. The state changes in the road system under the disturbance factors are shown in Fig 1. where t 0 , t 1 , t 2 represent the time when the disturbance starts, reaches the lowest state, and recovers to the original level, respectively. G(t 0 ), G(t 1 ), G(t 2 ) are the corresponding system states at the aforementioned moments. In the process of disturbance, the resilience index (19) of the system's average accumulation performance is shown in Eq. (7).
The resources and time consumed in the process from recession to recovery should be considered when calculating the traffic resilience. In this study, the recovery days were simply taken as the index to calculate road traffic resilience. It is added to the risk as one of the weight factors of the risk calculation; that is, the constraint of road traffic toughness is considered in the path selection.

III. MODEL BUILDING A. PROBLEM DESCRIPTION
Vehicle routing problem refers to the optimal collection or distribution route in a central warehouse for a fleet serving a group of customers in a transportation network [1]. In this study, low-carbon and road traffic resilience are added to the traditional hazardous materials transport route problem, and road traffic resilience is considered as one of the weights of the transportation risk calculation for the first time. Owing to the complex actual situation, the travel time between two customers may change in real time during transportation, and the volume of hazardous materials will change with the fluctuation of market demand. In the process of selecting the path, the resilience values of the nodes are integrated into the resilience values of the road section, which is one of the conditions of the optimal path. Simultaneously, considering the constraints of the time window, vehicle capacity, and other aspects, this study builds a hazardous materials transport route model with transport cost, risk and carbon emission as the target, which transports a variety of hazardous materials from multiple starting points. The road network structure designed to address this problem is illustrated in Fig. 2.
The network structure is composed of three main parts: supplier nodes, retail and common nodes. The supplier and retail nodes have time windows. In the regulated road network considering various factors such as time, demands, and resilience limits, a model aiming at transport costs, risks and carbon emissions is built for path planning from multiple starting points for many types of hazardous materials.

B. ASSUMPTION AND PARAMETERS
For be convenience in building the model and obtaining the final solution, some assumptions are made for the model as follows: 1) Compared with the traditional transport of hazardous materials, this problem involves the transport of many types of hazardous materials, so that a retailer can be visited for many times. However, when transporting the same types of hazardous materials, a retailer can only be visited once.
2) Each vehicle has a capacity limit, and the total amount of materials transported each time should not be exceed the capacity limit.
3) Each path must start and end with the same supplier node.
4) The travel time between two points connected in the path is uncertain. 5) A vehicle can only visit retail nodes within the capacity limit, and the number of retails visited is uncertain. 6) A vehicle can only pass the node once during transportation, and the number of vehicles is unlimited.
7) The demand of each retail is uncertain, and each retail has its own time window.
To facilitate a description of the transport of hazardous materials in this study, the required parameter symbols are defined in Table1.

C. MODEL 1) TRANSPORT COST
In many studies, transport costs are mainly considered as fixed and variable costs of vehicles. Fixed costs mainly refer to related personal expenses, vehicle management, insurance costs, etc. Variable costs include fuel consumption cost, tire consumption costs, etc. The transportation costs calculation formula is as follows: Owing to many uncertain factors, it is difficult to obtain accurate results in risk assessment. This study adopted the two-factor method [20] for risk modeling.
where R represents risk, P represents the probability of risk occurrence, and C represents the consequence. In this study, consequence refers to the number of exposed populations. Many studies [5] and [12] do not consider the change in the quantity of hazardous materials in transit, but in reality, the greater the quantity, the more serious the consequences. In this study, the risk is in a state of dynamic change with the choice of the route. If shipment is completed at the last retail node, the risk is set to zero when the vehicle returns to the initial node. Then, it is assumed that the path set before the shipment is completed is path, the probability of accident is Pro, the influence radius is r, and W is the weight matrix of road traffic resilience. The calculation formula is as follows: s , s ∈path ∩ L where the original value of D is the original quantity of vehicle m and s is the retailer set belonging to the set path. In transportation, D becomes smaller after unloading through each retail node, which reduces the risk of subsequent transportation.

3) CARBON EMISSION
For carbon emissions, the estimation method described in the literature [19] was adopted. Where T m represents the turnover of vehicle m,q represents the average loading capacity of the vehicles, and B represents the carbon emission factor.
A hazardous material transportation path planning model was constructed, which had three objectives: transportation cost, risk, and carbon emissions. Constraints such as the time, load, and flow are as follow: Constraints (8)-(10) refer to the three objectives of transportation cost, risk, and carbon emissions, respectively. Constraints (11) and (12) indicate that when transporting the same type of hazardous materials, one retailer can only be served by one vehicle. Constraint (13) represents the flow constraint. Constraint (14) indicates that the demand for retailers should match the quantity of materials used in transportation. Constraint (15) indicates that the load of each vehicle should be consistent with the demand of retailers. Constraint (16) indicates that the load on each vehicle should not exceed its loading capacity. Constraint (17) indicates that each vehicle can carry only one type of hazardous material. Constraints (18)- (24) are time constraints, in which the travel time is uncertain. Constraints (18) and (19) indicate that the vehicle should arrive at the retail before the end of the time window, where T is a large positive number. Constraint (20) represents the retail service time, and Constraints (21) and (22) are constraints on the arrival and departure times of the supplier, where T is a large positive number. Constraints (23) and (24) refer to the time at which the vehicle leaves the supplier or returns to the supplier during its time window.

IV. ALGORITHM DESIGN
The TSP is a typical NP-hard problem. This study considers the transportation of a variety of materials from multiple suppliers under the uncertainty of demand and time. The complexity is high, and it is difficult to use the traditional operational research method to solve this problem. Therefore, an ACO is designed to obtain the Pareto solution set of transportation cost, risk, and carbon emissions under constraint conditions (11)- (24). ACO is a probabilistic algorithm that simulates ant foraging behavior and can be used to determine the optimal path. In the process of finding the food source, ants use pheromones as an index, and pheromones change dynamically. New pheromones are added to the environment, while old pheromones are lost [17]. By continuously stacking pheromones, more ants are attracted to release them, thereby forming a positive feedback process that moves the result in the direction of the optimal solution. The combination of the local pheromone update and global pheromone update in the algorithm can accelerate the convergence speed, enhance the randomness of the algorithm, and prevent premature algorithm phenomena.
ACO is widely used in practical complex industrial problems, particularly TSP problems. Compared with other optimization algorithms, ACO has fewer control parameters, and the parameters are conveniently adjusted and optimized; thus, it has attracted increasing attention [32], [33], [34], [35], [36]. Therefore, we used the ant colony optimization algorithm to solve the multi-objective dangerous goods transportation optimization problem. In this study, we improve the ant colony algorithm, which can perform local and global pheromone update to find better solutions faster.
In this study, the road network structure is fixed. Compared with other path-planning articles, the selection of nodes should not only rely on the selection probability, but also consider the constraints of the road network structure. In the transportation process, because vehicle load changes dynamically and risk changes dynamically. The pheromone of the previous iteration affects the solution of the next iteration, but multiple factors (such as cost, risk, and carbon emissions) are considered in the selection probability and pheromone calculation. These factors restrict each other, and appropriate ACO parameters are obtained through multiple experiments, which can effectively prevent the algorithm from falling into local optima. In addition, after each ant completes the task, the local information is updated. In each iteration, the optimal solutions of all the ants are updated globally, which can accelerate the convergence speed of the algorithm. The main process of the algorithm is as follows Step 1: Set the iteration termination condition and the number of ant colony k.
Step 2: Initialize. First, the pheromone matrix is initialized and τ ij =0 (the pheromone from i to j is set to zero). Second, the starting supplier node is randomly selected. In addition, the type of materials transported was selected randomly because the simultaneous transportation of different hazardous materials may produce significant consequences. Therefore, only one type of material can be transported. If the demand of all retails for the type of materials is not zero, the type of material can be selected.
Step 3: Select the node. The probability of selecting the next node is determined according to the transfer probability, and the transfer probability formula is shown in Equation (25): s ∈ the next set of optional points (25) τ ij is the pheromone from node i to node j, θ and β are the parameters of the algorithm, η ij is usually the reciprocal of distance in traditional ACO. However, in this study, the value is the reciprocal of the product of distance and risk in Eq. (26), which can consider many factors in selecting nodes to find high-quality solution sets quickly.
After obtaining the probability of all the next optional nodes, the next node is selected using the roulette. Choose the cumulative probability of node i, J i = i j=1 p k mj , randomly generate a number r ∈ (0, 1), if the J i−1 < r < J i , select node i.
Step 4: Tabu the points selected in Step 3 that meet the constraint conditions in the tabu table.
Step 5: Determine whether a full path has been completed. In all constraint conditions whether there are alternative retail nodes, if yes, then jump to Step 3, until all retails meet the conditions are selected. After the last retail, the materials have been finished, and the shortest path algorithm returns to the original supplier. The full path was considered complete.
Step 6: Local pheromone update. After all materials are delivered, the ant completes the task. Then, the pheromone matrix was updated. The formula is shown in (27): For pheromone increment τ ij , cost, risk and carbon emissions are all added to it for pheromone update, and the three are normalized with a weight of 1. This is an improvement that is different from other ACO papers [5] and [12], where cost, risk and carbon emissions are added to the pheromone iteration to accelerate convergence and obtain Pareto optimal solution. After updating the local pheromone, the number of ants is updated to k = k+1.
Step 7: Non-dominated sorting. The results of all ants are sorted in a non-dominated manner to obtain different hierarchical results. The results of the first layer are the Pareto optimal solution sets for this iteration.
Step 8: Global pheromone update. The first layer results of the non-dominant sort in step 7 are used as the optimal solution to update the pheromone, and the formula is the same in Eq. (27).
Step 10: Check whether the iteration termination conditions are satisfied. If yes, the program is exited and the results are saved. Otherwise, return to the Step 2 to continue the update iteration. A flow chart of the ACO is shown in Fig. 4.
The pseudocode of the improved ACO is shown in Fig. 5.

V. EXPERIMENTAL EXAMPLE AND RESULTS
In this study, an experimental example was designed to verify the effectiveness of the proposed model and algorithm.
To ensure the reliability of the example, a road network structure with 30 nodes was randomly generated. See Fig 3. For other parameters, see Table 2-4, Fig 3 has 30 nodes, including three supplier nodes (green), ten retail nodes(yellow), and the rest are common nodes(blue).  Tables 2 and 3 show the time windows of supplier nodes and retail nodes and retail demand, respectively. The retail demand satisfies the linear uncertainty distribution. Table 4 lists the distance and average population density of the nodes.
The proposed problem was executed on an AMD Ryzen 7 4700U (Radeon Graphics 2.00 GHz) PC with a Windows 10 operating system.
The travel time in transportation was calculated by the distance and speed, of which the speed was set to 40. Consider-   ing the disturbance of the actual road conditions, an uncertain variable is added in the time calculation, which conforms to the normal uncertainty distribution with an expectation of 0.05 and a standard deviation of 0.1. Therefore, the travel time is uncertain. Each vehicle had a loading capacity of 13, a service time of 1 at the retail nodes and the influence radius is 0.3.
In addition, the parameters of ACO (θ, β, ρ, Q) are very important, including the information heuristic factor, expectation heuristic factor, information volatile factor, and pheromone strength. If θ and ρ are too large, the algorithm is not easy to converge. On the contrary, the algorithm will fall into local optimum. If β and Q are too small, the algorithm cannot converge easily. By contrast, the algorithm falls into a local optimum. Many experiments have shown that when θ = 0.7, β = 5.0, ρ = 0.5, Q = 100, the performance of the algorithm is better.
After one thousand iterations, the Pareto optimal solution set was obtained, as shown in Table 6. Fig 6 shows the results of the normalization of the Pareto optimal solution set.  Table 6 shows, from left to right, the number of vehicles completing the transportation task, the total cost, total risk, total carbon emission, total time, total distance in the transportation process and the node resilience values of all paths.
In Table 6 and Fig 6, except for the total risk, the change trends of the other results are approximately the same. As the number of vehicles increased, the total cost, total carbon emissions, total time, and total distance also increased, but the total risk decreased. This is because there is a cargo weight factor in total risk. The greater the number of vehicles, the less cargo is shared by each vehicle on average. The smaller the cargo weight factor, the lower the risk of each road section, and the total risk decreases accordingly. Solution 1 has the minimum total risk, and all other values are maximum. Solution 9 has the maximum total risk, and all other values are minimum. Solution 1 is compared to solution 9, the number of vehicles, total cost, total carbon emission, total time and total distance increased by 44.44%, 45.80%, 75.21%, 72.56% and 77.30% respectively, and the total risk decreased by 203.77%. Although all values except risk increased, the reduction in risk was significantly higher than the increase in the other values.
When the number of vehicles was the same, the total risk was inversely correlated with the total distance, whereas the remaining values were positively correlated with the total distance. This is because with an increase in distance, the in-transit time of vehicles increases, and the time period for hazardous materials to be exposed to the road network becomes longer, which may increase the possibility of risk occurrence and the total risk value increases.
An increase in the resilience value in the selected path indicates that the path has a strong ability to resist risks and recover to the original state after the occurrence of risks; thus, the corresponding risk value becomes smaller. In practice, decision-makers can choose a solution according to the degree of importance they can accept.   Taking solution 9 as an example, the path is presented in Table 5 and Fig 8. The Pareto optimal solution set obtained after several ACO calculations is shown in Fig 7. The same Pareto   solution set can be obtained after each operation, indicating that the algorithm exhibits strong stability and validity.
We ran the proposed improved ACO, genetic algorithm (GA) and simulated annealing algorithm (SA) three times and compared the best experimental results, as shown in Fig. 9. In addition, we calculated and compared the convergence indices [23] for the three algorithms. The calculation formula is shown in (28).
where is the calculated Pareto solution set, and * is the reference point set that is obtained uniformly from the Pareto fronts. | * | denotes the number of Pareto fronts. The smaller the index IGD, the better the algorithm convergence. The convergence indices of the improved ACO, genetic algorithm, and simulated annealing algorithm are 0.183, 0.211, and 0.205, respectively. It can be concluded that the ACO designed in this study has a preferable convergence and validity.
Many points with different colors and shapes represent the Pareto optimal solution set obtained by the three algorithms in Fig 9. In combination with Table 7, the solution set of the improved ACO has 12 solutions, whereas the GA and SA have 8 and 9 solutions, respectively. Thus, the improved ACO has a wider solution space. In terms of the quality of the solutions, the improved ACO usually has better solutions. The highlighted bold parts in Table 7 are optimal solutions of GA and SA, but they are superior to ACO only in terms of cost, carbon emissions, and higher risk. Therefore, the improved ACO has better solution quality and space, and better feasibility and stability. This can provide decision-makers with more feasible and high-quality solutions for practical application problems, which has great practical value.

VI. DISCUSSION AND CONCLUSION
In this study, a multi-objective transport model of hazardous materials was developed considering road traffic resilience and low-carbon in the case of uncertain demand and time. The load changes during transportation; therefore, the risk is also dynamically changed. The model considers the transportation cost, risk, and carbon emissions as the objective. Applying the parallel computing power of the improved ACO to obtain the Pareto optimal solution, compared with the GA and SA, the improved ACO has better solution quality and solution space, which reflects the superiority and reliability of the proposed algorithm. This provides decision-makers with better and more diverse solutions in practical situations. The main contributions of this study are as follows: 1) For the first time, road traffic resilience is considered as one of the risk weight factors in the transport of hazardous materials.
2) The integration of low-carbon factors into the transport of hazardous materials and increasing research on the transport of hazardous materials in the low-carbon field.
3) Consider the transport of multi-starting point and multikind hazardous materials.
With an increase in the number of vehicles, transportation cost and carbon emissions increase and risk value decreases, but the reduction in risk value is obviously higher than the increase in transportation cost and carbon emissions. By comparing some solutions, the risk reduction range is 3-5 times the increase range of other values, indicating that risk can be effectively reduced by considering resilience in the transport of hazardous material.
Although the effectiveness and reliability of the proposed algorithm have been verified, there are still some limitations, and much work will be carried out in the future research. As follows: 1) In the calculation of road traffic resilience, only the consumed days of recovery were considered, and the resources consumed by recovery were not considered. In future studies, the days consumed and the resources should be considered together.
2) The supply quantity was unlimited. Future research can consider the limitations of supply and possible impact and corresponding solutions when demand cannot be met.
3) The model considered only two uncertain factors (demand and time). A variety of uncertain factors can be added in the future, such as the environment and driver behavior, which can enrich the model to make it more realistic. 4) We can cite the latest studies about path planning to go deep into the study.