A Multiobjective Route Robust Optimization Model and Algorithm for Hazmat Transportation

. Aiming at route optimization problem of hazardous materials transportation in uncertain environment, this paper presents a multiobjective robust optimization model by taking robust control parameters into consideration. The objective of the model is to minimize not only transportation risk but also transportation time, and a robust counterpartof the model is introduced through applying the Bertsimas-Sim robust optimization theory. Moreover, a fuzzy C-means clustering-particle swarm optimization (FCMC-PSO) algorithm is designed, and the FCMC algorithm is used to cluster the demand points. In addition the PSO algorithm with the adaptive archives grid is used to calculate the robust optimization route of hazmat transportation. Finally, the computational results show the multiobjective route robust optimization model with 3 centers and 20 demand points’ sample studied and FCMC-PSO algorithm for hazmat transportation can obtain different robustness Pareto solution sets. As a result, this study will provide basic theory support for hazmat transportation safeguarding.


Introduction
Hazardous materials (hazmat) refer to products with flammable, poisonous, and corrosive properties that can cause casualties, damage to properties, and environmental pollution and require special protection in the process of transportation, loading, unloading, and storage.In recent years, the demand for hazmat has increased, its freight volume has increased year by year, and the potential transportation risk is also expanding.Practice has proved that the optimization of the transportation route of hazmat can effectively reduce the transportation risk, and it has significant influence to ensure the safety of people along route and protect the surrounding ecological environment.
Many scholars have studied the transportation route optimization problem for hazmat.Rhyne (1994) conducted a statistical analysis of the hazmat transport accident using the diffusion formula [1].Gordon et al. (2000) used the British transport accident to prove the different risk level in different places, which provided a realistic basis for the evacuation of the population and the balance of risks during the accident [2].Power et al. (2000) established a risk-cost analysis model and conducted a systematic analysis to find out the transport plan [3].Wu et al. (2002) studied the vehicle routing optimization problem of hazmat transportation with multidistribution center, and the clustering algorithm was used to solve the complex model [4].Kara et al. (2003) presented a risk model based on the primitive roads to reduce the error of risk estimation, and then they proposed the minimum risk route selection algorithm [5].Fabiano et al. (2005) used the experimental data to calculate the probability and consequence of the transport accident on the specific route of the experimental area and got the optimal transportation scheme [6].Erkut and Ingolfsson (2005) studied the hazmat transport route by considering the shortest distance, the smallest population exposure, and the smallest accident risk [7].Bubbico et al. (2006) analyzed the transportation risk of hazmat and obtained the safety route algorithm by using the experimental data from Italy [8].Wei et al. (2006) analyzed the route selection model under the influence of time-varying conditions without considering the uncertain risk [9].Verma (2009) presented a biobjective optimization model considering the cost and the risk and used the boundary algorithm to solve the model [10] .Wang et al. (2009) established a hazmat transportation path optimization model based on a geographic information system [11].Jassbi et al. (2010) developed a multiobjective optimization framework for the hazmat transportation by minimizing the transport mileage, the number of affected residents, social risks, the probability of accidents, and so on [12].Pradhananga et al. ( 2014) created a dual-objective optimization model with time windows for the hazmat transportation and designed a heuristic algorithm to solve this model [13].Suh-Wen Chiou (2016) proposed a dual-objective and dual-level signal control policy for the hazmat transportation [14].Assadipour et al. (2016) proposed a toll-based dual-level programming method for the hazmat transportation and designed hybrid speed constraints for multiobjective particle swarm optimization algorithm [15].Pamucar et al. (2016) proposed a multiobjective route planning method for hazmat transportation and designed a solution algorithm combining neurofuzzy and artificial bee colony approach [16].Mohammadi et al. (2017) studied the hazmat transportation under uncertain conditions using a mixed integer nonlinear programming model and the metaheuristic algorithm was utilized to solve this model [17].Kheirkhah et al. (2017) established a bilevel optimization model for the hazmat transportation, two heuristic algorithms were designed to solve the dual-level optimization model, and some randomly generated problems were used to verify the applicability and validity of the model [18].Bula et al. (2017) focused on the heterogeneous fleet vehicle routing problem and designed a variant of the variable neighborhood search algorithm to solve the problem [19].Ma et al. (2013) studied the route optimization models and algorithms for hazardous materials transportation under different environments [20].Ma et al. (2018) proposed a road screening algorithm and created a distribution route multiobjective robust optimization for hazardous materials [21].Obviously, although the route optimization of hazmat transportation has made many achievements, the robustness of solution is rarely considered in this field.
Ben-Tal and Neimirovski (1998) proposed robust optimization theory based on ellipsoid uncertainty set [22].Bertsimas and Sim (2003) further put forward adjustable robustness robust optimal theory [23].Based on the adjustable robust optimization theory, this paper will propose a route robust optimization model for hazmat transportation with multiple distribution centers, and the author also designed a kind of improved PSO algorithm.
The rest of this paper is structured as follows: Section 2 introduces the hazmat transportation route problem and establishes a multiobjective route robust optimization model of hazmat transportation.Section 3 presents the FCMC-PSO algorithm.Section 4 is the case study.At the end of the paper, the conclusion is proposed.

Multiobjective Route Robust Optimization Model of Hazmat Transportation
. .Problem Definition.Hazmat transportation route robust optimization for multidistribution center is defined as follows: there are several hazmat distribution centers, and each distribution center owns enough hazmat transport vehicles; meanwhile, multiple need points exist which should be assigned to the relevant hazmat distribution center.Vehicles from distribution center will service the corresponding demand points.Each vehicle can service several customer demand points while each customer demand point only can be serviced by one vehicle.After completing the transport mission, the vehicles must return to the distribution center [24].
Uncertainty of hazmat transport refers to the uncertainty of the transportation time and transportation risk, which may be caused by the traffic accident, weather, and traffic density of the road.Compared to ordinary goods transport, hazmat transport is more complicated, and it needs more security demands.Therefore, it is needed to set the goal of minimizing the total hazmat transportation risk.In the process of hazmat transportation, transport time reduction is also necessary.As a consequence, this paper will target minimizing the total transportation time.In conclusion, the scientific transportation routes should be found to guarantee the hazmat transported safely and quickly.
. .Model . . .Assumption.There are a few assumptions in this study: ( = {1, 0} , ∀ ∈ , ∀ ∈ , ∀ ∈  0 , ∀ ∈   (11) where the objective function (1) minimizes the total transportation risk of hazmat [25], and the objective function (2) minimizes the total transportation time of hazmat.Constraint (3) as the load constraint that means any vehicle of any distribution center should satisfy the corresponding load constraint, namely, cannot overload.Constraint (4) means that vehicle number of the distribution center is limited, and vehicles arranged to transport the hazmat should not exceed the distribution center owning.Constraint (5) expresses that the transportation risk of each section must be less than or equal to threshold r max set by decision makers.Constraint (6) expresses that the transportation risk of each route must be less than or equal to threshold R max set by decision makers.Constraint (7)  . . .Robust Counterpart Model.Each objective function of the above multiobjective robust model corresponds to parameter Γ.The purpose is to control the degree of conservatism of the solution.Objective functions ( 1) and ( 2) of the robust optimization model contain "max" extreme value problem.Set feasible solution set X vrp to satisfy all constraints, and robust discrete optimization criterion is used to transform the multiobjective route robust optimization model, and a new robust counterpart of the model is as follows [26,27].

Algorithm
In this section, we propose FCMC-PSO algorithm to solve the multiobjective route optimization problem of hazmat transportation in uncertain environment.The demand points is clustered by the fuzzy C means algorithm, and the transportation route for each demand points is determined based on the adaptive archives grid multiobjective particle swarm optimization [28][29][30].
. .Fuzzy C-Means Clustering.Suppose that n data samples are  = { 1 ,  2 , ⋅ ⋅ ⋅ ,   }, C( ≤C≤n) is the number of types into which the data samples are to be divided, { 1 ,  2 , ⋅ ⋅ ⋅ ,   } indicating the corresponding C categories, U is its similar classification matrix, the clustering centers of all categories are { 1 ,  2 , ⋅ ⋅ ⋅ ,   }, and   (  ) is the membership degree of sample   to category A k (abbreviated as   ).Then the objective function can be expressed as follows: where , it is the synthetic weighted value of the transport risk and time after nondimensionalization between the i-th sample X i and the k-th category center point; m is the characteristics number of the sample; b is the weighting parameter, and the value range is 1 ≤ b ≤ ∞.The fuzzy C-means clustering method is to find an optimal classification, so that the classification can produce the smallest function value J b .It requires a sample for the sum of the membership degree of each cluster is 1, which is satisfied: Formulas ( 17) and ( 18) are used to calculate separately the membership degree   of the sample X i for the category A k and C clustering centers {  }: Let   = { | 2 ≤  < ;   = 0}, for all the i categories,  ∈   ,   = 0.
Using formulas ( 17) and ( 18) to repeatedly modify the cluster center, data membership, and classification, when the algorithm is convergent in theory, we can get the membership degree of the cluster center and each sample for each pattern class; thus the division of the fuzzy clustering is completed.
. .Multiobject PSO Algorithm.Particle swarm algorithm is derived from the study of the predatory behavior of birds.It is used to solve the problem of path optimization.Each particle in the algorithm represents a potential solution, and the fitness value for each particle is determined by the fitness function, and the value of fitness determines the pros and cons of the particle.The particle moves in the N-dimensional solution space and updates the individual position by the individual extremum and the group extremum.In the algorithm, the velocity, position, and fitness value are used to represent the characteristics of the particle.The velocity of the particle determines the direction and distance of the particle movement, and the velocity is dynamically adjusted with the moving experience of its own and other particles.Once the position of the particle is updated, the fitness value will be calculated, and the individual extremum and the population extremum are updated by the fitness values of the new particles, the individual extremum, and the population extremum.Multiobjective particle swarm optimization algorithm is a method based on particle swarm optimization algorithm to solve multiobjective problem.At the same time, the best location of multiple populations exists in the population, and the optimal positions of multiple particles themselves are also found in the iterative process.Therefore, gbest and pbest also need to adopt certain strategies to choose.Aiming at the robust optimization model of hazmat transportation, the key elements of multiobjective particle swarm optimization are as follows.
( ) Individual Coding.In this paper, the method of particle encoding adopts integer encoding, and each particle represents the experienced demand point.For example, when the number of required points is 9, the individual coding is [9 2 1 4 3 6 7 5 8], indicating that the requirement point traversal starting from the distribution center, followed by 9 2 1 4 3 6 7 5 8, and ultimately return to the distribution center.
( ) Fitness Value.In the hybrid particle swarm algorithm, the fitness value is the criterion of judging the quality of the particle.And the fitness function is to facilitate the search and improve the performance of the algorithm.In the paper, the fitness value of the particle is expressed by the objective function of the built model.
( ) Crossover Operation.Crossover operation is the process of replacing the partial structure of the parent individual and reorganizing the new individual.The design of the crossover operation is related to the representation of the coding, the cross-operation design based on the coarranged coding method of the demand point and the distribution center [31,32].The method of integer crossing is adopted.Set the two individuals of the parent as [9 2 1 4 3 6 7 5 8] and [8 3 2 4 1 5 9 7 6], Firstly, two crossover positions are selected, and then the individual is crossed; the operation process can be seen in Figure 1.
The new individuals need to be adjusted if there is a duplicate position, and the adjustment method is to replace the repeated demand points by using the absence of demand points in individuals.For the new individual 1, there are mappings about 2 to 3, 9 to 6, and 5 to 7. The specific adjustments process can be seen in Figure 2.
The strategy of retaining outstanding individuals is used for the owned new individuals, and the particles are updated only when the new particle fitness is better than the old particles.( ) Mutation Operation.The mutation of the particles is to make some changes in some genes of the particle; the mutation can increase the ability of searching particles and increase the diversity of the populations, to avoid falling into the local optimal situation.
The variation is also related to the way the particle is encoded.Based on the lease point and the dispatch center, there are many methods about the coarranged coding and the variation.In this paper, the variation method adopts the individual internal exchange method.For example, for an individual [9 2 3 4 1 5 6 7 8], at first, the mutated positions pos1 and pos2 are selected randomly; and then the positions of two variants are swapped.Assuming that the selected mutation positions are 2 and 6, the mutation operation process can be seen in Figure 3.
The strategy of retaining outstanding individuals is used for the owned new individuals, and the particles are updated only when the new particle fitness is better than the old particles [33].

( ) Multiobject PSO Algorithm Based on Adaptive Archives
Grid.Multiobjective PSO based on adaptive archives grid is a particle swarm optimization algorithm proposed by Coello and Lechuga to solve multiobjective problem [34].Its basic idea is to divide the target space into several hypercubes and to judge the number of noninferiority contained in each hypercube to maintain the external files.In each iteration, if the file does not exceed the given size, then a new nondominated solution will be added to the file.If the file has been filled, the file is maintained according to the density of noninferior solution contained in the hypercube, the noninferior solution is removed from the high-density hypercube, and the noninferior solution with low density is added to ensure the diversity of the population.The algorithm steps are as follows.
Step .Create and initialize a group so that the ex archives external file is empty.
Step .Evaluate all particles and add the noninferior solution to the external file.
Step .Maintain external files according to the adaptive grid method.
Step .Select gbest and pbest for each particle.
Step .Updating the velocity and position of the particles according to the speed formula and the position formula of the particle swarm.Step .Make sure the particles exist in the search space.
Step .If the termination condition is satisfied, the output result algorithm is terminated; if it is not satisfied, Step 2 is executed and the execution is continued.

Case Study
There are 3 distribution centers and 20 demand points, the maximum load for each transport vehicle is 8 ton, and each distribution center has adequate hazmat.The distribution centers are marked as a, b, and c, and the hazmat demand points are marked as 1, 2,..., 20.The demand amount of each demand point is shown in Table 1, and transportation risk and time from each distribution center to demand points and between the demand points are, respectively, shown in Tables 2 and 3.The risk and time nominal values are given in Tables 2  and 3.The transportation risk deviation r (0 ≤ r < 0.5  ) and transportation time deviation t (0 ≤ t < 0.5  ) are provided.
We use FCMC algorithm to calculate the demand points clustering results, and the results can be seen in Table 4. Based on the cluster results, we use the multiobjective PSO to solve the robust optimization problem for each distribution center.The parameters of the algorithm are set as follows: population size is 100, maximum evolution generation is 1000, inertia weight is 0.6, accelerated factor is 1.7, crossover rate is 0.95, and mutation rate is 0.09.The Pareto solution set with different robust control parameters can be obtained by calculation, which are showed in Tables 5-13 and Figures 4-6.
In Table 5, the Pareto optimal solutions sets are obtained by the program running on the condition that uncertain transportation risk and uncertain time are taken the nominal values.It is known from Table 5 that when Γ = 0, the program finds 6 Pareto solutions; in the encoding sequence, the first number a represents the distribution center.While, in the decoding sequence, every alphabet represents one vehicle, the figures behind alphabet show the customer demand points and the corresponding service order.Such as in Table 5, the first decoding sequence indicates that 3 vehicles are needed and each vehicle corresponds to a subroute, they are, respectively, a→9→20→a, a→14→3→a, and a→15→16→a, and it is clearly known that all transport vehicles from distribution center a, after serving the allocated customer demand points, eventually return to distribution center a.When the robustness control parameters Γ = 10 and Γ = 20, Pareto optimal solutions sets obtained by the program running are given in Tables 6 and 7, respectively.Tables 8, 9, and 10 present the optimal Pareto solutions sets for distribution center b when the robustness control parameters are taken 0, 10, and 20.Tables 11, 12, and 13 present the  optimal routes have not changed; when Γ = 20, the total risk optimal route has changed relatively small; however, the total time optimal routes have not changed when the robustness control parameters Γ = 0, Γ = 10, and Γ = 20, which means that uncertainty data are taken nominal value, and the total risk and time optimal routes always have some robustness.
Aiming at distribution center b, when the robustness control parameters Γ = 10 and Γ =20, the total risk optimal routes have not changed; when the robustness control parameters Γ = 0, Γ=10, and Γ =20, the total time optimal vehicle routes do not have any change, which shows the total risk optimal distribution route is relatively stable.Aiming at distribution center c, when the robustness control parameters Γ = 10 and Γ = 20, the total risk optimal vehicle routes do not have any change, but the total time optimal vehicle routes are changed when the robustness control parameters Γ = 0, Γ = 10, and Γ = 20, which shows when Γ = 10 the total risk optimal distribution route is relatively stable, to a certain extent, which can be selected as distribution route, and the total time optimal vehicle route has relatively weak robustness, and if more stable distribution route is needed, it will be expected to increase the robustness control parameter to find the strong robustness distribution route.For solution robustness, as the robustness control parameters get bigger, the corresponding Pareto solution robustness should be enhanced in theory.After a large number of analyses for the actual situation and basic data, we can determine the robustness control parameter value and obtain the corresponding candidate route set.
The strength Pareto genetic algorithm (SPEA) is used to test the efficiency of the FCMC-PSO algorithm.The parameters of the SPEA are set as follows: population size is 50, maximum evolution generation is 300, crossover rate is 0.8, and mutation rate is 0.05.The solution set can be obtained by calculation.The results are shown in Table 17.Compared with SPEA, the convergence iterations and operation time of FCMC-PSO algorithm are reduced.The results show that the FCMC-PSO algorithm designed in this paper not only can obtain a more satisfactory solution but also has faster convergence speed compared with the SPEA.

Conclusion
Hazmat transportation route optimization is an important link to ensure transportation safety of hazmat.In this paper, we take the hazmat transportation route problem with multidistribution center as the research object, considering the transportation risk and transportation time.In addition, an adjustable robustness transportation route multiobjective robust optimization model is established in the end.Speaking of the solution, the FCMC-PSO algorithm is designed in this research.The demand points were assigned through FCMC algorithm in which transportation time and transportation risk are considered.The multiobjective route robust optimization model is solved by multiobject PSO algorithm based on adaptive archives grid.In the end, the example shows that the robust optimization model and FCMC-PSO algorithm can obtain different robustness Pareto solution sets.The robust optimization transportation routes of hazmat will provide basic theory support for safeguarding the transportation safety of hazmat.
In this paper, we only consider the two uncertainties including transportation risk and transportation time.However, there may be some other uncertainties such as customer demand and service time window in the real world.According to this situation, we need to establish the corresponding robust model for future study.Although most of the hazmats are transported by road transportation, the hazmat transportation risks induced by other modes cannot be ignored.The optimization research of multimodal transport modes for the hazmat needs to be further studied.
: deviation of the variable transport risk to its nominal value from customer demand points i and j, where    ≥ 0.   : travel time nominal value from customer demand points i and j.   : deviation of variable travel time to its nominal value from customer demand points  and , where    ≥ 0. t : variable transport risk from customer demand points i and j, where t ∈ [  ,   +    ](    ≥ 0).
(8)icates every hazmat vehicle departing from distribution center should return back to the original distribution center after finishing the transportation task.Constraint(8)and Constraint (9) guarantee (10) each demand point is served once and by one vehicle form a distribution center.Constraint(10)indicates hazmat transport vehicles cannot depart from one distribution center but back to another one.Constraint (11) defines the 0-1 integer variable    , if the route of Hazmat transport vehicle k from distribution center  containing the segment from node i to node j,    =1; else    =0.

Table 2 :
The nominal transportation risk value of hazmat.

Table 3 :
The nominal transportation time value of hazamt.

Table 4 :
Demand points clustering result.