Multi-Depot Vehicle Routing Optimization Considering Energy Consumption for Hazardous Materials Transportation

Focusing on the multi-depot vehicle routing problem (MDVRP) for hazardous materials transportation, this paper presents a multi-objective optimization model to minimize total transportation energy consumption and transportation risk. A two-stage method (TSM) and hybrid multi-objective genetic algorithm (HMOGA) are then developed to solve the model. The TSM is used to find the set of customer points served by each depot through the global search clustering method considering transportation energy consumption, transportation risk, and depot capacity in the first stage, and to determine the service order of customer points to each depot by using a multi-objective genetic algorithm with the banker method to seek dominant individuals and gather distance to keep evolving the population distribution in the second stage, while with the HMOGA, customer points serviced by the depot and the serviced orders are optimized simultaneously. Finally, by experimenting on two cases with three depots and 20 customer points, the results show that both methods can obtain a Pareto solution set, and the hybrid multi-objective genetic algorithm is able to find better vehicle routes in the whole transportation network. Compared with distance as the optimization objective, when energy consumption is the optimization objective, although distance is slightly increased, the number of vehicles and energy consumption are effectively reduced.


Introduction
Hazardous materials can be toxic, corrosive, inflammable, and explosive, but also are necessary for the development of industry and agriculture.How to choose one or more transportation routes for hazardous materials with low energy consumption and low risk has become a very important topic with practical significance to ensure personal, property, and environmental security as well as economic sustainability.
According to statistics, the daily amount of hazardous materials being transported in China exceeds 1 million tons, and the total amount per year exceeds 400 million tons [1].These large quantities of hazardous materials form a dangerous flowing source on the roads.If an accident occurs, it could cause damage to nearby people and the surrounding environment.Alp et al. [2] studied the effects of meteorologic conditions and wind direction probabilities on the consequences of accidents, and proposed a three-integral mathematical evaluation model for calculating individual and social risks.Leonelli et al. [3] proposed a new personal and social risk assessment model, introducing risk factors such as transportation mode, hazardous materials category, meteorologic conditions, wind direction probability, and seasonal attributes, and considered the impact of population distribution, which greatly improved the accuracy of risk assessment for hazardous materials transportation.Based on a statistical analysis of accidents, Fabiano et al. [4] discussed risk factors from the aspects of road characteristics, weather conditions, and traffic conditions, and proposed a risk assessment model for the hazardous materials transportation process at accident sites.Based on the route segment-specific (location-specific) accident rate, Chakrabarti et al. [5] estimated route segment total risk to measure the average number of persons likely to be exposed to all possible consequence scenarios by computing and comparing the loss of containment and spillage probabilities for different route segments.Taking the length and quality of roads, the population density in different segments of the roads, types of hazardous materials, and types of vehicles into consideration, Faghih-Roohi et al. [6] proposed a dynamic value-at-risk model.Mohammadi et al. [7] designed a reliable hazardous materials transportation network based on four factors: distance from the incident source, number of people exposed to the risk, number of hazardous materials shipments passing that segment, and the characteristics of the road type.With the rapid development of Geographic Information System (GIS) technology, it has also been applied to risk assessment for hazardous materials transportation (Bubbico et al. [8], Sahnoon et al. [9], Chen et al. [10]).Visual map display and data analysis functions provide convenience and improve accuracy for risk assessment during the transportation of hazardous materials.It can be seen that scholars have made many contributions to the problem of transporting hazardous materials.The risk of transporting hazardous materials is related to many uncertain factors, and it is difficult to measure accurately, but it is still a factor that must be considered.
In addition to risk, energy consumption is also an important factor to consider in reducing the transportation costs of enterprises, reducing carbon emissions, and achieving sustainable development of hazardous materials transportation.There are many factors that affect transportation energy consumption, but they can be roughly classified into categories such as vehicle construction, transportation environment, and driver strategies.Factors influencing fuel consumption have been studied by Bigazzi et al. [11], Demir et al. [12], and Suzuki [13].Trucks have to overcome various obstacles while transporting hazardous materials, such as speed, road gradient, and payload.These are key factors that affect energy consumption.Bektaş and Laporte [14] analyzed the relationship between vehicle load, speed, and total cost, and proposed that minimizing the cumulative load alone does not necessarily lead to energy minimization, particularly when there are time window restrictions.Demir et al. [15] derived an optimum driving speed and showed that reductions in emissions could be achieved by varying speed over a network.It should be noted that the optimum driving speed varies to a certain degree between geographic areas due to speed limits and traffic density.Given the fact that the shortest distance may not be the optimal solution for the purpose of lowering fuel consumption, Xiao et al. [16] developed route schedules with lower fuel costs by better managing the trade-off between the total distance and the priority of serving customers with larger demands.Suzuki [17] indicated that significant savings in fuel consumption may be realized by delivering heavy items in the early segments of a tour and lighter items in the later segments.In fact, transportation costs are more related to energy consumption than distance.Therefore, for hazardous materials transportation, considering the factors that affect energy consumption, it is more practical to find one or more low-energy vehicle routes as the optimization objective.
The transportation route is a key factor to guarantee that hazardous materials are transported safely and economically, which many experts and scholars have been studying.Batta et al. classified transportation risk into road risk and node risk [18].Karkazis et al. [19] formulated the routing problem of transporting hazardous materials as an optimization problem, with the objective of minimizing transportation risk and costs.Meng et al. [20] used the dynamic programming method to solve the transportation route optimization problem.Liu et al. [21] applied a multilevel fuzzy comprehensive evaluation method to optimizing hazardous materials transportation routing.Das et al. [22] explored the routing problem in a transportation network with the capacity limit, and developed a multi-objective algorithm to find nondominated solutions.Pradhananga et al. [23] chose minimum transportation time and risk as the optimization objective to search for Pareto optimal solutions.Given the unscientific nature of simply evaluating the total risk of a fleet as a whole without considering individual vehicle risk, Wang et al. [24] developed a two-stage exact algorithm based on the ε-constraint method to solve the proposed problem.Alrukaibi et al. [25] designed a risk/cost algorithm using available data from Kuwait.For their study, incident probability, incident consequence, and risk assessment were used as the algorithm's main criteria to identify available alternative routes.Mohammadi et al. [7] proposed a mathematical model for designing a reliable hazardous materials transportation network on the basis of hub location topology under uncertainties; to cope with the uncertainties, they provided a solution framework based on an integration of the chance-constrained programming approach.Obviously, these studies on the routing problem of hazardous materials transportation have been relatively sufficient, and have gone through a study process from a single objective to multi-objective optimization and from certain to uncertain environments.However, energy consumption is not considered in the model, and to solve the model, most studies turned the multi-objective model into a single objective, and finally got an optimal solution instead of the Pareto optimal solution set.
All of the above studies on the vehicle routing problem (VRP) for hazardous materials transportation assumed that there was only one depot.However, there may, in fact, be more than one vehicle depot in a city, and it is of more practical significance to study the VRP with multiple depots for hazardous materials transportation.At present, only a few studies have focused on the multi-depot vehicle routing problem (MDVRP) for hazardous materials transportation.Zhao et al. [26] considered the return trips between collection centers and recycling centers, and developed a multi-objective model with minimization of cost and risk.This optimization problem with two objectives was then transformed into a single objective problem.Du et al. [27] developed a fuzzy bilevel programming model in which the upper-level formulation allocated customers to depots and the lower level determined the optimal routing for each depot.However, the study only considered transportation risk minimization as the optimization objective.In fact, for hazardous materials transportation, in addition to risk, energy consumption is also a major factor that must be considered.This is mainly due to the fact that energy consumption is a key factor in reducing the transportation costs of enterprises, reducing carbon emissions, and achieving sustainable development.To this end, this study aims to find the optimal vehicle routings for the MDVRP considering both transportation risk and energy consumption.
Many heuristic methods have been developed in the context of the VRP, including the tabu search algorithm [28], the genetic algorithm [29], the ant colony algorithm [30], and so on.Compared to the VRP, the MDVRP is more complicated, because it needs to consider which customer point should be serviced by which depot, in addition to considering the vehicle routing for each depot.It is necessary to coordinate the delivery tasks among multiple depots.Gillett and Johnson presented a clustering procedure and a sweep heuristic for each depot [31].Salhi and Sari proposed a multilevel heuristic method [32].Giosa et al. [33] summarized and proposed a number of heuristics for the two-stage approach to solving the MDVRP.In order to solve the emergency vehicle routing problem with multiple depots, Qin et al. [34] divided it into two steps, and used the nearest assignment and average distance method to transform the MDVRP into multiple VRPs with a single depot in the first step, and then adopted a genetic algorithm to solve the VRPs with a single depot in the second step.Ho et al. [35] developed two hybrid genetic algorithms.The major difference between the two algorithms is that the initial solutions were generated randomly in the first algorithm, and the initialization procedure was incorporated into the Clarke and Wright savings method and the nearest neighbor heuristic in another algorithm.In summary, the methodologies for MDVRP can be categorized into one-stage and two-stage solving methods.
To sum up, a number of studies have explored the traditional VRP with only one depot.These studies move from single objective optimization to multi-objective optimization and from single solution to a nondominated solution set.However, only a few studies have focused on the MDVRP for hazardous materials transportation, and all of them consider a single objective or do not consider transportation energy consumption, due to the complexity of the problem.Furthermore, these studies rarely consider all depots and customer points to find the optimal Pareto solution set to reduce both transportation risk and energy consumption.The existing research has mainly adopted the two-stage method, which cannot make multiple depots coordinate transportation and dynamically solve the problem.In addition, the existing research only considers one objective or converts multiple objectives into one to solve the problem, and they rarely use multi-objective algorithms to seek an optimal Pareto solution set.Different from the current studies on the MDVRP for hazardous materials transportation, this study aims to obtain the Pareto optimal solution set of noncomparable transportation risk and energy consumption.Finally, a two-stage algorithm (TSM) and a one-stage hybrid multi-objective genetic algorithm (HMOGA) are designed, and the solution optimality between these two algorithms is compared and analyzed.
The rest of this paper is organized as follows: Section 2 describes the MDVRP for hazardous materials transportation, Section 3 establishes an optimization model to minimize transportation risk and energy consumption, Section 4 presents a TSM and HMOGA to solve the proposed model, Section 5 gives numerical examples to analyze the effectiveness of the above two algorithms, and Section 6 gives conclusions.

Problem Description
As the name implies, the VRP for hazardous materials transportation with one depot is a problem of vehicle scheduling for a specific depot.All customer points are served by a specific depot.However, in reality, it is very common that multiple depots coordinate the transportation and distribution of hazardous materials.The MDVRP for hazardous materials transportation can be defined as the following: (1) There exist several hazardous materials depots, and each depot has enough vehicles to transport the hazardous materials.(2) Multiple customer points exist and will be assigned to different depots.
(3) Each vehicle will service the corresponding customer points and can service several customer points, while each customer point can be serviced only one time by one vehicle.(4) After the delivery task is finished, all transportation vehicles will return to their hazardous materials depot.
Figure 1 is a sketch of the MDVRP for hazardous materials transportation.It can be seen that with three depots, named a, b, and c, there are 15 customer points, and each customer point can be serviced by only one depot.According to the MDVRP, multiple depots can cooperate with the transportation environment according to customer needs to achieve the global optimal objective.
In the MDVRP, because of the strongly corrosive, highly toxic, explosive, and flammable characteristics of hazardous materials, transportation risk is an optimization objective that must be considered.For enterprises, minimizing transportation costs is paramount.However, low risk means high cost to some extent.Therefore, the MDVRP for hazardous materials transportation should be solved according to the idea of multi-objective optimization, and it is very important to find the Pareto solution set.In many reports in the literature, transportation cost is usually measured by distance and time.However, the shortest transportation distance or time in practice does not mean that transportation cost is the lowest.Compared to distance or time, transportation energy consumption is of great significance to reducing business operating costs, lowering carbon emissions, and achieving sustainable development.Therefore, this paper sets two optimization objectives to minimize risk and energy consumption for hazardous materials transportation.

Energy Consumption Evaluation
Rolling resistance, air resistance, and ramp resistance are the main external forces that affect the energy consumption of truck transportation.According to the research in [36], the following formula for calculating power consumption to overcome resistance is used: Equation ( 1) is the power used to overcome rolling resistance during the transportation of hazardous materials, where ( ) Equation ( 2) is the power consumed to overcome air resistance during transportation, where ( ) a P v is the power consumed to overcome air resistance at speed v (kW); a is the air resistance coefficient; a ρ is air density (kg/m 3 ); and z is the windward area of the vehicle, also known as the front area (m 2 ).

( )
, , Equation ( 3) is the power consumed by the truck from node i to node j to overcome ramp resistance during transportation, where ( ) , , g P m v i represents power consumed to overcome ramp resistance at ramp i and transport speed v, and ij i represents the slope of node i to node j.
In addition, the acceleration of trucks while transporting hazardous materials requires a large amount of fuel.According to the literature [37], the energy consumption of trucks accelerating from 0 to v is:

Energy Consumption Evaluation
Rolling resistance, air resistance, and ramp resistance are the main external forces that affect the energy consumption of truck transportation.According to the research in [36], the following formula for calculating power consumption to overcome resistance is used: Equation ( 1) is the power used to overcome rolling resistance during the transportation of hazardous materials, where P r (m, v) represents the power consumed by transport equipment with a total weight of m to overcome rolling resistance at speed v (kW); c r represents the rolling resistance coefficient; g represents gravitational acceleration (m/s 2 ); m represents the total weight of the transport equipment, including the weight of the cargo m c and the vehicle M(t); and v represents the instantaneous speed of the transport equipment (m/s).
Equation ( 2) is the power consumed to overcome air resistance during transportation, where P a (v) is the power consumed to overcome air resistance at speed v (kW); a is the air resistance coefficient; ρ a is air density (kg/m 3 ); and z is the windward area of the vehicle, also known as the front area (m 2 ).
Equation ( 3) is the power consumed by the truck from node i to node j to overcome ramp resistance during transportation, where P g (m, v, i) represents power consumed to overcome ramp resistance at ramp i and transport speed v, and i ij represents the slope of node i to node j.
In addition, the acceleration of trucks while transporting hazardous materials requires a large amount of fuel.According to the literature [37], the energy consumption of trucks accelerating from 0 to v is: The movement of the truck is a real-time change process.Considering that the data are not easy to obtain, this paper uses the average speed v and the average slope i to calculate the average energy consumption during transportation.Therefore, the average energy consumption of hazardous materials transportation between 2 nodes is: where d ij represents the distance of transportation from node i to node j (km); v represents the average speed during transportation (m/s); i represents the average slope; and n num represents the average number of accelerations per kilometer (times/km).

MDVRP Multi-Objective Model for Hazardous Materials Transportation
Hazardous materials transportation risk is usually defined as the product of the accident rate and the potential loss during transportation.The probability of accidents occurring is often affected by time, weather conditions, vehicle type, loading form, and road conditions, and the severity of the consequences are related to the number of people affected, property losses, and weather conditions.In this paper, we use the risk measurement model proposed by Erkut et al. [38] to calculate hazardous materials transportation risk, and regard the number of people affected on both sides of the road as the main indicator for measuring risk: where r ij is the risk of road segment (persons), a ij is the accident rate, p ij is the population density (persons/km 2 ), and r d is the impact radius of the accident (m).
The MDVRP multi-objective model for hazardous materials transportation is: ∑ In the above model, as the objective functions, Equations ( 7) and ( 8), respectively, are used to minimize total transportation risk and energy consumption, and it is worth mentioning that the risk of the delivery vehicle on its return to the depot from the final customer point is not calculated in total risk, while energy is calculated in total transportation energy.In the objective function, S 0 is the depot set, here S 0 = {i|i = 1, 2, • • • , m}; m is the total number of hazardous materials depots; S 1 is the customer point set, S 1 = {i|i = 1, 2, • • • , n}, where n is the number of customer points; S is the transportation network nodes set, here S = S 0 ∪ S 1 ; V d is the vehicle set of hazardous materials depot d and ijk defines the 0-1 integer variable.Equation ( 9) is the load constraint, which means any vehicle of any depot should satisfy the corresponding load constraint, namely, it cannot overload.In the equation, q j is the quantity demanded at customer point i and L d k is the load capacity of vehicle k from depot d.Equation ( 10) means the number of vehicles in the depot is limited, and the number of vehicles arranged to transport hazardous materials should not exceed the number owned by the depot.Equation (11) indicates that every vehicle departing from the depot should return to the original depot after finishing the transportation task.Equations ( 12) and ( 13) guarantee that each customer point is serviced once and by one vehicle from a depot.Equation ( 14) means the vehicles cannot depart from one depot and return to another depot.Equation ( 15) constrains the maximum risk on the segment.Equation ( 16) denotes that if the route of vehicle k from depot d for transporting hazardous materials contains the road segment from node i to node j, x d ijk is equal to 1, else x d ijk equals 0.

Solving Methods
Different from single objective optimization, possible conflicts exist between multiple objectives in the multi-objective optimization, and an improved sub-objective might cause the performance of the other one to decrease.Multi-objective optimization generally gets a set of Pareto solutions; elements of the set are called Pareto optimal solutions.In reality, decision-makers select one or more solutions from a Pareto optimal solution set as the optimal solution of a multi-objective optimization problem based on personal preference.For a multi-objective optimization problem, the typical algorithms include niched Pareto genetic algorithm (NPGA) [39], nondominated sorting genetic algorithm-II (NSGA-II) [40], and strong Pareto evolutionary algorithm (SPEA2) [41].These algorithms have better solution efficiency for particular problems.However, for a specific problem, these algorithms cannot be used directly, but often need to be modified according to the problem properties.
This paper mainly uses genetic algorithm theory to design a method of solving MDVRP for hazardous materials transportation.Based on biological characteristics, the genetic algorithm realizes diversity and global search by operating the population composed of potential solutions.It focuses on the set of individuals, which is consistent with the Pareto solution set of the multi-objective optimization problem.Moreover, the genetic algorithm does not need many mathematical prerequisites, but can deal with all types of objective functions and constraints, and has good performance for combinatorial optimization.

Two-Stage Method
Solving the MDVRP for hazardous materials transportation has a certain complexity, and the algorithm needs to solve two problems: one is to choose appropriate customer points for hazardous materials depots, and the other is to design a reasonable distribution service order for each depot to service the selected customer points.The method is designed in two stages: In stage one, we use the global search cluster to convert the problem into several multi-objective optimization problems with single depots, which decreases the solving complexity of the problem.In stage two, we design a multi-objective genetic algorithm to solve the routing problem transformed in stage one, then we can obtain the corresponding routing for each depot to service its customer points.

Global Search Cluster
Deciding which depot services which customer points is related to the customer point itself, the current depot, and other customer points that are serviced by the current depot.It also means that there is an inverse relationship between the risk and energy of the customer point, the current depot, and other customer points that are serviced by the current depot to determine whether the depot serves the customer point.By defining intimacy in order to determine which customer points the depot will service, customer points will be serviced by the depot with the highest intimacy.
Intimacy is defined as follows: (rW Here, on the hypothesis that customer point j is serviced by depot d, f (i, d) is defined as the affinity between depot d and customer point i; (rW) ij refers to dimensionless risk and energy-weighted average on network nodes i to j, r ij is the dimensionless nominal risk, and W ij is the dimensionless nominal energy.In the above equations, a, b, α, and β are the weighting coefficients used to adjust the factor magnitude; L d is the current capacity of the depot; S q (d) refers to the customer points assigned to depot d; the initial state for S q (d) is the depot itself; and |S q (d)| is the assigned customer point number to depot d.In this way, all customer points are divided into |S 0 | groups, where |S 0 | is the number of the depot.

Multi-Objective Genetic Algorithm
The designed multi-objective genetic algorithm constructs a Pareto optimal solution set by the banker method, using the gather density method to keep evolving the population distribution.The basic process of the algorithm is shown in Figure 2, in which Pop refers to population, |Pop| is the population size and its consistent size is N, Paretos is the Pareto optimal solution set, |Paretos| is the number of Pareto optimal solutions in the set and its value takes 2N as the maximum, and the end condition of the algorithm is the generation limit.Here, the evolutionary group is a population with a fixed number of chromosomes.The chromosomes of the offspring population are composed of all the Pareto individuals in the Pareto pool, all the individuals in the Pareto set of the current population, and the individuals in the current population selected by the selection operator.
(1) Encoding and Decoding Chromosomes We adopted a natural number coding method.For example, a network has one hazardous materials depot and eight customer points.If the depot is marked 0 and the customer points are respectively marked 1-8, the randomly generated sequence 0 3 5 1 4 2 8 7 6 becomes one of the chromosomes.Because in the model a few vehicles are needed to service the customer points, chromosomes obtained by this coding method also have to be decoded.The chromosomes are decoded based on a greedy strategy, by inserting the customer points into the route according to the order of the genes in the chromosome without violating load capacity constraint.If another customer point is added, it will violate the load constraint, which means that we should assign another vehicle to service this customer point.Assuming that there are eight customer points and their quantity demand, successively, is 3 tons, 2 tons, 4 tons, 1 ton, 2 tons, 2 tons, 2 tons, and 1 ton, and vehicle load capacity is 8 tons, the chromosome 0 3 5 1 4 2 8 7 6 decodes as follows: (

1) Encoding and Decoding Chromosomes
We adopted a natural number coding method.For example, a network has one hazardous materials depot and eight customer points.If the depot is marked 0 and the customer points are respectively marked 1-8, the randomly generated sequence 0 3 5 1 4 2 8 7 6 becomes one of the chromosomes.Because in the model a few vehicles are needed to service the customer points, chromosomes obtained by this coding method also have to be decoded.The chromosomes are decoded based on a greedy strategy, by inserting the customer points into the route according to the order of the genes in the chromosome without violating load capacity constraint.If another customer point is added, it will violate the load constraint, which means that we should assign another vehicle to service this customer point.Assuming that there are eight customer points and their quantity demand, successively, is 3 tons, 2 tons, 4 tons, 1 ton, 2 tons, 2 tons, 2 tons, and 1 ton, and vehicle load capacity is 8 tons, the chromosome 0 3   To chromosome A and B respectively, delete the genes behind the mating area that appeared in the mating area, then get their offspring chromosomes A" = 0 3 4 5 6 7 1 8 2, B" = 0 6 5 1 8 2 3 4 7.
(2) Mutation operator: reverse transcription method In the chromosome two different genes are first randomly selected except the genes representing hazardous materials depot, and then a reverse operation is performed on the selected genes.
(3) Constructing Pareto Optimal Solution Set The banker method is not a kind of backtracking method, and constructing new Pareto optimal solutions does not need to compare with the existing Pareto optimal solutions.Set Paretos as the Pareto optimal solutions set and Pop as evolution group set, then adopt the banker method to construct Pareto optimal solution set for Pop as follows: Initialize the Pareto optimal solutions set as Paretos.
Take out an individual (usually the first one) from Pop, as a banker, compare the banker with other individuals, and delete the individual dominated by the banker.
After the banker is compared with all individuals in Pop, if any individual dominates the banker, delete the banker, otherwise add the banker into Paretos.
Repeat steps and until Pop becomes null.
(4) Calculating Individual Gathering Density Evolving population distribution is an important part of the multi-objective optimization.Here, the gathering density method is used to keep the evolving population distribution, specifically adopting gather distance to represent the individual gather density, and the greater gather distance shows the smaller gather density.Before calculating each individual distance, we need to order the individual according to its objective function value.Set distance I i as the gather distance for individual i, and I j i as the value of objective j for individual i.If individual i contains m objectives, the gather distance for individual i will be calculated by Equation (20):

Hybrid Multi-Objective Genetic Algorithm
When solving the MDVRP for hazardous materials transportation, we need to consider the customer points serviced by each depot as well as the service order.Therefore, we designed a hybrid multi-objective genetic algorithm (HMOGA) to solve the MDVRP for hazardous materials transportation.The algorithm adopts the same method shown in Section 4.1.2in the aspects of Pareto optimal solution set construction and keeping the evolving population distribution.Other specific operations are as follows: (1) Chromosome Coding and Decoding The hybrid multi-objective genetic algorithm adopts a hybrid coding mode combining the customer point and distribution service.The coding mode divides the chromosome into two parts: part 1 is to select the customer points for each depot to service, with a length equal to the number of all customer points and gene values are depot numbers; part 2 is to determine the service order for each customer point, with a length also equal to the number of all customer points, while genes are customer points.Combining part 1 with part 2 and according to the loading capacity constraint of hazardous materials transportation vehicles, we adopt a greedy selection strategy to decode the chromosome.For example, if there are three depots and 10 customer points, the hybrid coding chromosomes are as follows: Chromosome part 1: a→c→b→c→a→b→a→b→c→b Chromosome part 2: 3→9→2→10→1→8→6→7→5→4 (2) Genetic Operator Design Similar to coding, genetic operators are also designed in two parts.Chromosome part 1 adopts two-point crossover and simple mutation operators; chromosome part 2 uses partially mapped crossover and reverse transcription operators.
(1) Crossover operator Two-point crossover operator: First, according to certain probability, freely choose two individual chromosomes as the parent generation, and select two genes in part 1 of the selected chromosomes as the intersections; second, swap the genes between the intersections of the selected chromosomes, generating the offspring chromosomes.Partially mapped crossover operator: Use the partially mapped crossover operator in Section 4.1.2to complete the crossover operation.
(2) Mutation operator Simple mutation operator: With a certain probability, randomly select a chromosome, freely choose a gene on part 1, and change the gene to its allele; if the selected gene is b, change it to a or c.Reverse transcription operator: same as reverse transcription operations in Section 4.1.2.

Cases Design
An enterprise has three hazardous materials depots and 20 customer points with hazardous materials.The depots are designated a, b, and c, and the customer points are numbered 1-20.Each depot uses trucks to deliver hazardous materials, and the dead weight, payload, capacity, engine efficiency, and other vehicle parameters of the trucks are identical.Assume that each depot has enough hazardous materials and transportation vehicles.This paper uses the affected population model in the traditional risk model to measure the risk of each road segment.In addition to the length of the road segment, the traditional affected population model also needs some relevant parameters such as affected radius, accident rate, and population density.For energy consumption, a middle-level method is used.Trucks perform delivery tasks by overcoming rolling resistance, air resistance, and ramp resistance during the transportation of hazardous materials.To calculate rolling resistance and ramp resistance, in addition to the total weight of the vehicle, the rolling resistance coefficient, the ramp resistance coefficient, the road ramp, and other parameters are required.For air resistance, the air resistance coefficient, the vehicle windward area, and the air density are required parameters.In order to better verify the solution effect of the designed methods, we designed two cases, case 1 and case 2.
In case 1, the risk triangle fuzzy number of the first three hazardous materials depots and 20 customer points in the literature [27] are used as the source data of road segment length, population density, accident rate, and road ramp from each depot to and between customer points.The road segment length is obtained by subtracting the lower limit from the upper limit of the fuzzy number; the population density around the road segment is obtained by subtracting the lower limit from the value with maximum possibilities of the fuzzy number and multiplying by 150; the accident rate is obtained by multiplying the value with maximum possibilities of the fuzzy number by 10 −5 ; the road ramp is obtained by subtracting the value with maximum possibilities from the upper limit of the fuzzy number.
In case 2, the road ramp is a random number between −5% and 5%, and the accident rate is a random number between 10 −6 and 5 × 10 −5 .The road segment length, the population density on both sides of the road, and the number of customer points are the same as the corresponding data in case 1.The road ramp and the accident rate are shown in Tables A1-A3 in Appendix A, where the letters indicate the depot number and the figures indicate the number of the customer point.
The parameters related to vehicles transporting hazardous materials in the two cases are shown in Table 1.According to the setting of the affected radius of hazardous materials in the literature [42], the accident affected radius in both cases is 1000 m.The demand amount of each customer point in the two cases is measured by weight (t), as shown in Table 2.

Results Analysis
The parameters for the global search clustering method in the first stage are set as a = 0.6, b = 0.4, α = 0.9, β = 0.1.Table 3 shows the customer points serviced by each hazardous materials depot in case 1 and case 2 by calculating the clustering algorithm proposed in the first stage.In the second stage, the parameters of the multi-objective genetic algorithm are set as follows: population size is 100, maximum evolution generation is 100, partially mapped crossover rate is 0.95, and reverse transcription mutation rate is 0.2.Based on the clustering results for depots a, b, and c in case 1 and case 2, we obtain their Pareto solution set of customer points serviced by each depot, as shown in Tables 4 and 5.
The parameters of the HMOGA are set as follows: population size is 200, maximum evolution generation is 200, two-point crossover rate is 0.6, partially mapped crossover rate is 0.9, simple mutation rate is 0.3, and reverse transcription mutation rate is 0.1.Tables 6 and 7 show the Pareto optimal solution sets obtained by solving the HMOGA in case 1 and case 2, respectively.In the "Routes" column of Tables 4-7, each letter and the figures behind it indicate the transport route of a delivery vehicle.For example, "a-1-4-19-20-a-14-17" in Table 4 indicates the need for two delivery vehicles.Its transportation routes are a-1-4-19-20-a and a-14-17-a.From the above tables, it can be seen that the methods designed here can all obtain the Pareto optimal solution set, and the solutions in the set are not comparable.Figures 3 and 4 are the Pareto optimal distribution curves by solving the HMOGA in case 1 and case 2, respectively.It can be seen from Figures 3 and 4 that the HMOGA can obtain the Pareto solution distributed in the whole solution space instead of focusing on a certain area, where the first point and the last point in each figure correspond to the optimal solution of energy consumption and risk, respectively.Decision-makers can select the appropriate routes for hazardous materials transport vehicles according to the actual situation.If the transportation task is currently carried out according to the fourth solution in Table 4, the decision-maker can choose the first solution in the table to carry out the transportation task when an accident occurs between customer points 2 and 3 and traffic is prohibited.
Figures 4 and 5 show optimal risk value and energy value, respectively, of each depot obtained by the two methods in case 1.  Figures 4 and 5 show optimal risk value and energy value, respectively, of each depot obtained by the two methods in case 1.
It can be seen from Figures 5 and 6 in case 1 that compared to the single objective risk and energy consumption optimal solutions, the HMOGA model performed nearly six times and 1.1% better than the TSM.The huge difference in Figure 5 is mainly due to the characteristics of data used in case 1.In the optimal risk solution, depot an optimized by HMOGA only serves three customer points (1, 4, and 6), and the accident rate and population density of the road segment between depot a and the three customer points are relatively small.However, there are six customer points (1,4,14,17,19,20) serviced by depot an obtained by the TSM in case 1, and the accident rate and population density between depot a and customer points 17 and 20 is relatively large.Therefore, the risk value of depot a obtained by the HMOGA is much smaller than the optimal risk value obtained by the TSM.Figures 4 and 5 show optimal risk value and energy value, respectively, of each depot obtained by the two methods in case 1.
It can be seen from Figures 5 and 6 in case 1 that compared to the single objective risk and energy consumption optimal solutions, the HMOGA model performed nearly six times and 1.1% better than the TSM.The huge difference in Figure 5 is mainly due to the characteristics of data used in case 1.In the optimal risk solution, depot an optimized by HMOGA only serves three customer points (1, 4, and 6), and the accident rate and population density of the road segment between depot a and the three customer points are relatively small.However, there are six customer points (1,4,14,17,19,20) serviced by depot an obtained by the TSM in case 1, and the accident rate and population density between depot a and customer points 17 and 20 is relatively large.Therefore, the risk value of depot a obtained by the HMOGA is much smaller than the optimal risk value obtained by the TSM.Figures 4 and 5 show optimal risk value and energy value, respectively, of each depot obtained by the two methods in case 1.
It can be seen from Figures 5 and 6 in case 1 that compared to the single objective risk and energy consumption optimal solutions, the HMOGA model performed nearly six times and 1.1% better than the TSM.The huge difference in Figure 5 is mainly due to the characteristics of data used in case 1.In the optimal risk solution, depot an optimized by HMOGA only serves three customer points (1, 4, and 6), and the accident rate and population density of the road segment between depot a and the three customer points are relatively small.However, there are six customer points (1,4,14,17,19,20) serviced by depot an obtained by the TSM in case 1, and the accident rate and population density between depot a and customer points 17 and 20 is relatively large.Therefore, the risk value of depot a obtained by the HMOGA is much smaller than the optimal risk value obtained by the TSM.It can be seen from Figures 5 and 6 in case 1 that compared to the single objective risk and energy consumption optimal solutions, the HMOGA model performed nearly six times and 1.1% better than the TSM.The huge difference in Figure 5 is mainly due to the characteristics of data used in case 1.In the optimal risk solution, depot an optimized by HMOGA only serves three customer points (1, 4, and 6), and the accident rate and population density of the road segment between depot a and the three customer points are relatively small.However, there are six customer points (1,4,14,17,19,20) serviced by depot an obtained by the TSM in case 1, and the accident rate and population density between depot a and customer points 17 and 20 is relatively large.Therefore, the risk value of depot a obtained by the HMOGA is much smaller than the optimal risk value obtained by the TSM.For optimal energy consumption, the optimization result of the HMOGA only saves 1.1% compared with the optimization result of the TSM.This is mainly due to the strong dependence of energy consumption on the length of the road segment.In case 1, the length and slope of the road segment have little changes, and the slope is only uphill, so the optimization effect of the HMOGA is not obvious.
Similarly, Figures 7 and 8 show the optimal risk value and energy value of each depot in case 2, respectively.In case 2, since the slope and accident rate of the road segment are randomly generated, for the optimal solutions of risk and energy consumption, the HMOGA has a nice optimization effect compared with the TSM.As can be seen from Figures 7 and 8, the optimal risk and energy consumption values obtained by the HMOGA are reduced by 12.4% and 54.8%, respectively, compared with the values obtained by the TSM in case 2. For optimal energy consumption, the optimization result of the HMOGA only saves 1.1% compared with the optimization result of the TSM.This is mainly due to the strong dependence of energy consumption on the length of the road segment.In case 1, the length and slope of the road segment have little changes, and the slope is only uphill, so the optimization effect of the HMOGA is not obvious.
Similarly, Figures 7 and 8 show the optimal risk value and energy value of each depot in case 2, respectively.For optimal energy consumption, the optimization result of the HMOGA only saves 1.1% compared with the optimization result of the TSM.This is mainly due to the strong dependence of energy consumption on the length of the road segment.In case 1, the length and slope of the road segment have little changes, and the slope is only uphill, so the optimization effect of the HMOGA is not obvious.
Similarly, Figures 7 and 8 show the optimal risk value and energy value of each depot in case 2, respectively.In case 2, since the slope and accident rate of the road segment are randomly generated, for the optimal solutions of risk and energy consumption, the HMOGA has a nice optimization effect compared with the TSM.As can be seen from Figures 7 and 8, the optimal risk and energy consumption values obtained by the HMOGA are reduced by 12.4% and 54.8%, respectively, compared with the values obtained by the TSM in case 2.  For optimal energy consumption, the optimization result of the HMOGA only saves 1.1% compared with the optimization result of the TSM.This is mainly due to the strong dependence of energy consumption on the length of the road segment.In case 1, the length and slope of the road segment have little changes, and the slope is only uphill, so the optimization effect of the HMOGA is not obvious.
Similarly, Figures 7 and 8 show the optimal risk value and energy value of each depot in case 2, respectively.In case 2, since the slope and accident rate of the road segment are randomly generated, for the optimal solutions of risk and energy consumption, the HMOGA has a nice optimization effect compared with the TSM.As can be seen from Figures 7 and 8, the optimal risk and energy consumption values obtained by the HMOGA are reduced by 12.4% and 54.8%, respectively, compared with the values obtained by the TSM in case 2. In case 2, since the slope and accident rate of the road segment are randomly generated, for the optimal solutions of risk and energy consumption, the HMOGA has a nice optimization effect compared with the TSM.As can be seen from Figures 7 and 8, the optimal risk and energy consumption values obtained by the HMOGA are reduced by 12.4% and 54.8%, respectively, compared with the values obtained by the TSM in case 2.
As to single depots, the optimal solution obtained by the TSM may be better than the HMOGA, but for the entire transportation network, the HMOGA can find better transportation routes than the TSM.Accordingly, to solve the MDVRP for hazardous materials transportation, the HMOGA has certain advantages in practical application.
For the MDVRP for hazardous materials transportation, most studies aim to reduce transportation costs, with distance as the main optimization objective.In order to compare the optimization effect when transportation distance and risk are used as the optimization objective, in addition to considering transportation risk, we test case 2. Table 8 shows the optimal vehicle route obtained by taking energy consumption and distance as the optimization objective in case 2. It can be seen from Table 8 that for the MDVRP for hazardous materials transportation, compared with distance as the optimization objective, when energy consumption is optimized, although total distance increases by 2.3%, total energy consumption reduces by 47.9% and the number of vehicles decreases by one.Therefore, it is more practical to consider the energy consumption of truck transportation in the MDVRP.

Conclusions
Hazardous materials transportation routing is one of the basic elements to assure safe transportation.In reality, it is very common that some depots coordinate the transportation and distribution of hazardous materials.
Designing scientific and reasonable transportation routes considering factors such as transportation risk and energy consumption can help deliver hazardous materials to each customer point safely, quickly, and economically.We took the MDVRP for hazardous materials transportation as the research object, considering transportation risk and energy consumption, and established a hazardous materials transportation routing multi-objective optimization model.In the model, the risk of the road segment is measured by the number of people who could be affected if an accident occurred, and energy consumption is measured by the trucks overcoming resistance during transportation.
As for solving the model, we designed a TSM and a HMOGA.In the TSM, the purpose of the first stage is to find the set of customer points serviced by each depot through the global search clustering method considering transportation energy consumption, transportation risk, and depot capacity, and to determine the service order for customer points to each depot by using a multi-objective genetic algorithm with the banker method to seek dominant individuals and gather distance to keep evolving population distribution in the second stage.The HMOGA combines the solution method of the two stages into one stage.In the design of the algorithm, customer points serviced by the depot and serviced orders are optimized simultaneously.In the end, we compared the above two methods through two cases, and the results show that the methods can obtain different Pareto solution sets.The TSM and HMOGA have their own advantages, but the HMOGA is able to find better transportation routes for the whole transportation network.
In addition, when the length and slope of the road segment do not change much, energy consumption does not change much, and decision-makers can use transport risk as the main basis for decision-making.When the difference between the slope and length of the road segment is large, it is important to optimize the energy consumption to find the Pareto optimal solution set for hazardous materials transportation.The decision-maker can select the appropriate solution to perform the transportation task in the Pareto solution according to the actual situation.Compared with distance, when energy consumption is the optimization objective, although the distance is slightly increased, the number of vehicles and amount of energy consumption are effectively reduced.
Due to the lack of actual data, such as accident probability and slope information of road segments, we designed the example by simply improving the data in the existing literature.Although the rationality of the model and algorithm is verified, a large-scale actual road network study is still necessary for further research.In addition, we found the Pareto optimal solution set, and for each solution in the set, it is difficult to measure which one is better.Hence, in combination with other conditions, choosing one or more Pareto solutions as hazardous materials transportation routes is the focus of the next phase of research.
represents the power consumed by transport equipment with a total weight of m to overcome rolling resistance at speed v (kW); r c represents the rolling resistance coefficient; g represents gravitational acceleration (m/s 2 ); m represents the total weight of the transport equipment, including the weight of the cargo c m and the vehicle M(t); and v represents the instantaneous speed of the transport equipment (m/s).

Transportation route 1 21 Figure 2 .
Figure 2. Basic flow of the multi-objective genetic algorithm.

Figure 2 .
Figure 2. Basic flow of the multi-objective genetic algorithm.(2)Designing a Genetic Operator The selection operator, crossover operator, and mutation operator play important roles in deciding solution efficiency of the genetic algorithm.The basic flow of the selection operator is shown in Figure 1.

( 1 )
Crossover operator: partially mapped crossover Randomly select a mating area in the selected two chromosomes, for example, A = 0 To chromosomes A and B, add their own mating area to the other chromosome, for instance, to the two chromosomes in , after operation to obtain A = 0

Figure 3 .
Figure 3. Pareto optimal solution distribution solved by HMOGA in case 1.

Figure 4 .
Figure 4. Pareto optimal solution distribution solved by HMOGA in case 2.

Figure 5 .
Figure 5. Minimum risk obtained by the two algorithms for each depot in case 1.

Figure 3 .
Figure 3. Pareto optimal solution distribution solved by HMOGA in case 1.

Figure 3 .
Figure 3. Pareto optimal solution distribution solved by HMOGA in case 1.

Figure 4 .
Figure 4. Pareto optimal solution distribution solved by HMOGA in case 2.

Figure 5 .
Figure 5. Minimum risk obtained by the two algorithms for each depot in case 1.

Figure 4 .
Figure 4. Pareto optimal solution distribution solved by HMOGA in case 2.

Figure 3 .
Figure 3. Pareto optimal solution distribution solved by HMOGA in case 1.

Figure 4 .
Figure 4. Pareto optimal solution distribution solved by HMOGA in case 2.

Figure 5 .Figure 5 .
Figure 5. Minimum risk obtained by the two algorithms for each depot in case 1.

Figure 6 .
Figure 6.Minimum energy obtained by the two algorithms for each depot in case 1.

Figure 7 .
Figure 7. Minimum risk obtained by the two algorithms for each depot in case 2.

Figure 8 .
Figure 8. Minimum energy consumption obtained by the two algorithms for each depot in case 2.

Figure 6 .
Figure 6.Minimum energy obtained by the two algorithms for each depot in case 1.

Figure 6 .
Figure 6.Minimum energy obtained by the two algorithms for each depot in case 1.

Figure 7 .
Figure 7. Minimum risk obtained by the two algorithms for each depot in case 2.

Figure 8 .
Figure 8. Minimum energy consumption obtained by the two algorithms for each depot in case 2.

Figure 7 .
Figure 7. Minimum risk obtained by the two algorithms for each depot in case 2.

Figure 6 .
Figure 6.Minimum energy obtained by the two algorithms for each depot in case 1.

Figure 7 .
Figure 7. Minimum risk obtained by the two algorithms for each depot in case 2.

Figure 8 .
Figure 8. Minimum energy consumption obtained by the two algorithms for each depot in case 2.

Figure 8 .
Figure 8. Minimum energy consumption obtained by the two algorithms for each depot in case 2.

Table 2 .
Demand amount of each customer point in case 1 and case 2.

Table 3 .
Customer points for each depot in case 1 and case 2.

Table 4 .
Pareto solution sets for customer points serviced by each depot in case 1.

Table 5 .
Pareto solution sets for customer points serviced by each depot in case 2.

Table 7 .
Pareto optimal solution set solved by the HMOGA in case 2.

Table 8 .
Optimal vehicle routes obtained by using energy consumption and distance as the optimization objectives in case 2.