# A Multiobjective Route Robust Optimization Model and Algorithm for Hazmat Transportation

Authors: Changxi Ma; Wei Hao; Ruichun He; Bahman Moghimi
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

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 counterpart of 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.


---

## Body

## 1. 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: Section2 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.


## 2. Multiobjective Route Robust Optimization Model of Hazmat Transportation

### 2.1. 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.


### 2.2. Model

#### 2.2.1. Assumption

There are a few assumptions in this study:(1) Multiple hazmat distribution centers are existent(2) The supply of hazmat distribution centers is adequate(3) Vehicle loading capacity is provided and the demand of each customer is specified(4) Multiple vehicles of the distribution center can service the customer(5) The transportation risk and transportation time are identified among the customer demand points, but they are uncertain number as interval number


#### 2.2.2. Symbol Definition

S 0: set of the Hazmat distribution centers, where S0={i∣i=1,2,…,m} shows that the number distribution center is m and the sequence number of nodes set is 1…m.S 1: set of customer demand points, where S1={i∣i=1,…,n} shows that the number of customer demand points is n and the sequence number of nodes set is 1…n.S: all nodes set in the transportation network, where S=S0∪S1.V d: available transportation vehicle set in the hazmat distribution center, where V={k∣k=1,2,…K},d∈S0.E: road section set among nodes.q i: demand of customer demand point i.L k d: maximum load of transport vehiclekfrom hazmat distribution centerd.r ~ i j: variable transport risk from customer demand pointsi andj, where r~ij∈[rij,rij+r⌢ij](r⌢ij≥0).r i j: transportation risk nominal value from customer demand pointsi andj.r ⌢ i j: deviation of the variable transport risk to its nominal value from customer demand pointsi andj, where r⌢ij≥0.t i j: travel time nominal value from customer demand pointsi andj.t ⌢ i j: deviation of variable travel time to its nominal value from customer demand points i and j, where t⌢ij≥0.t ~ i j: variable transport risk from customer demand pointsi andj, where t~ij∈[tij,tij+t⌢ij](t⌢ij≥0).J i r: the set of columns which all uncertain data r~ij belonging to the ith row of the variable risk matrix are in, here Jir≤n.Γ i r: parameter Γir∈[0,Jir] to adjust robust risk of robust discrete optimization method and control the risk degree of conservatism, where decimal is permitted.Γ i r: maximum integer less than Γir.ψir: the set of column subscriptsjof uncertain data r~ij of line i in the variable risk matrix r~ij.J i t: the set of columns with all uncertain data t~ij belonging to theith row of the variable time matrix, where Jit≤n.Γit: parameter Γit∈[0,Jit] to adjust robust time of robust discrete optimization method and control the time degree of conservatism, where decimal is permitted.Γ i t: the maximum integer less than Γit.ψit: the set of column subscriptj of uncertain data r~ij of linei in variable time matrix t~ij.


#### 2.2.3. Multiobjective Route Robust Optimization Model

(1) m i n  Z 1 = ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ j ∈ S 1 r i j x i j k d + m a x Ψ i r ∪ m r ∣ Ψ i r ⊆ J i r , Ψ i r = Γ i r , m r ∈ J i r ∖ Ψ i r ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ j ∈ Ψ i r r ⌢ i j x i j k d + ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ m r ∈ J i r ∖ Ψ i r Γ i - Γ i r ⌢ i m r x i m r k d - ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S r i 0 x i 0 k d (2) m i n  Z 2 = ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ j ∈ S t i j x i j k d + m a x Ψ i t ∪ m t ∣ Ψ i t ⊆ J i , Ψ i t = Γ i t , m t ∈ J i t \ Ψ i t ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ j ∈ Ψ i t t ⌢ i j x i j k d + ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ m t ∈ J i t \ Ψ i t Γ i - Γ i t ⌢ i m t x i m t k d (3) s.t.  ∑ i ∈ S ∑ j ∈ S 1 ∑ d ∈ S 0 x i j k d q j ≤ L k d ,  ∀ k ∈ V d (4) ∑ j ∈ S 1 ∑ k ∈ V d x i j k d ≤ V d ,  ∀ i ∈ S 0 ,   ∀ d ∈ S 0 (5) m a x i , j ∈ E r i j x i j k d ≤ r m a x ,  ∀ i ∈ S ,   ∀ j ∈ S ,   ∀ d ∈ S 0 ,   ∀ k ∈ V d (6) m a x i , j ∈ E ∑ i ∈ S 0 ∪ S 1 ∑ j ∈ S 0 ∪ S 1 r i j x i j k d - ∑ i ∈ S 1 ∑ j ∈ S 0 r i j x i j k d ≤ R m a x ,  ∀ i ∈ S ,   ∀ j ∈ S ,   ∀ d ∈ S 0 ,   ∀ k ∈ V d (7) ∑ j ∈ S 1 x i j k d - ∑ j ∈ S 1 x j i k d = 0 ,  ∀ i ∈ S 0 ,   ∀ d ∈ S 0 ,   ∀ k ∈ V d (8) ∑ i ∈ S ∑ d ∈ S 0 ∑ k ∈ V d x i j k d = 1 ,  ∀ j ∈ S 1 (9) ∑ j ∈ S ∑ d ∈ S 0 ∑ k ∈ V d x i j k d = 1 ,  ∀ i ∈ S 1 (10) ∑ j ∈ S 0 x i j k d = ∑ j ∈ S 0 x j i k d = 0 ,  ∀ i ∈ S 0 ,   ∀ d ∈ S 0 ,   ∀ k ∈ V d (11) x i j k d = 1,0 ,  ∀ i ∈ S ,   ∀ j ∈ S ,   ∀ d ∈ S 0 ,   ∀ k ∈ V dwhere 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 thresholdrmax set by decision makers. Constraint (6) expresses that the transportation risk of each route must be less than or equal to thresholdRmax set by decision makers. Constraint (7) indicates 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 that 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 xijkd, if the route of Hazmat transport vehiclek from distribution center d containing the segment from nodeito nodej, xijkd =1; else xijkd =0.


#### 2.2.4. 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 setXvrp 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].Objective function is(12)Rr=Γ^r^l+min⁡∑m=1n2rmxm+∑m=1lr^m-r^lxm(13)Tt=Γ^t^l+min⁡∑m=1n2tmxm+∑m=1lt^m-t^lxmConstraint condition is(14)x∈XvrpThen, the optimal objective function value can be obtained asR∗=minl=1,⋯n2+1⁡R(l) and T∗=minl=1,⋯n2+1⁡T(l).


## 2.1. 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.


## 2.2. Model

### 2.2.1. Assumption

There are a few assumptions in this study:(1) Multiple hazmat distribution centers are existent(2) The supply of hazmat distribution centers is adequate(3) Vehicle loading capacity is provided and the demand of each customer is specified(4) Multiple vehicles of the distribution center can service the customer(5) The transportation risk and transportation time are identified among the customer demand points, but they are uncertain number as interval number


### 2.2.2. Symbol Definition

S 0: set of the Hazmat distribution centers, where S0={i∣i=1,2,…,m} shows that the number distribution center is m and the sequence number of nodes set is 1…m.S 1: set of customer demand points, where S1={i∣i=1,…,n} shows that the number of customer demand points is n and the sequence number of nodes set is 1…n.S: all nodes set in the transportation network, where S=S0∪S1.V d: available transportation vehicle set in the hazmat distribution center, where V={k∣k=1,2,…K},d∈S0.E: road section set among nodes.q i: demand of customer demand point i.L k d: maximum load of transport vehiclekfrom hazmat distribution centerd.r ~ i j: variable transport risk from customer demand pointsi andj, where r~ij∈[rij,rij+r⌢ij](r⌢ij≥0).r i j: transportation risk nominal value from customer demand pointsi andj.r ⌢ i j: deviation of the variable transport risk to its nominal value from customer demand pointsi andj, where r⌢ij≥0.t i j: travel time nominal value from customer demand pointsi andj.t ⌢ i j: deviation of variable travel time to its nominal value from customer demand points i and j, where t⌢ij≥0.t ~ i j: variable transport risk from customer demand pointsi andj, where t~ij∈[tij,tij+t⌢ij](t⌢ij≥0).J i r: the set of columns which all uncertain data r~ij belonging to the ith row of the variable risk matrix are in, here Jir≤n.Γ i r: parameter Γir∈[0,Jir] to adjust robust risk of robust discrete optimization method and control the risk degree of conservatism, where decimal is permitted.Γ i r: maximum integer less than Γir.ψir: the set of column subscriptsjof uncertain data r~ij of line i in the variable risk matrix r~ij.J i t: the set of columns with all uncertain data t~ij belonging to theith row of the variable time matrix, where Jit≤n.Γit: parameter Γit∈[0,Jit] to adjust robust time of robust discrete optimization method and control the time degree of conservatism, where decimal is permitted.Γ i t: the maximum integer less than Γit.ψit: the set of column subscriptj of uncertain data r~ij of linei in variable time matrix t~ij.


### 2.2.3. Multiobjective Route Robust Optimization Model

(1) m i n  Z 1 = ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ j ∈ S 1 r i j x i j k d + m a x Ψ i r ∪ m r ∣ Ψ i r ⊆ J i r , Ψ i r = Γ i r , m r ∈ J i r ∖ Ψ i r ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ j ∈ Ψ i r r ⌢ i j x i j k d + ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ m r ∈ J i r ∖ Ψ i r Γ i - Γ i r ⌢ i m r x i m r k d - ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S r i 0 x i 0 k d (2) m i n  Z 2 = ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ j ∈ S t i j x i j k d + m a x Ψ i t ∪ m t ∣ Ψ i t ⊆ J i , Ψ i t = Γ i t , m t ∈ J i t \ Ψ i t ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ j ∈ Ψ i t t ⌢ i j x i j k d + ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ m t ∈ J i t \ Ψ i t Γ i - Γ i t ⌢ i m t x i m t k d (3) s.t.  ∑ i ∈ S ∑ j ∈ S 1 ∑ d ∈ S 0 x i j k d q j ≤ L k d ,  ∀ k ∈ V d (4) ∑ j ∈ S 1 ∑ k ∈ V d x i j k d ≤ V d ,  ∀ i ∈ S 0 ,   ∀ d ∈ S 0 (5) m a x i , j ∈ E r i j x i j k d ≤ r m a x ,  ∀ i ∈ S ,   ∀ j ∈ S ,   ∀ d ∈ S 0 ,   ∀ k ∈ V d (6) m a x i , j ∈ E ∑ i ∈ S 0 ∪ S 1 ∑ j ∈ S 0 ∪ S 1 r i j x i j k d - ∑ i ∈ S 1 ∑ j ∈ S 0 r i j x i j k d ≤ R m a x ,  ∀ i ∈ S ,   ∀ j ∈ S ,   ∀ d ∈ S 0 ,   ∀ k ∈ V d (7) ∑ j ∈ S 1 x i j k d - ∑ j ∈ S 1 x j i k d = 0 ,  ∀ i ∈ S 0 ,   ∀ d ∈ S 0 ,   ∀ k ∈ V d (8) ∑ i ∈ S ∑ d ∈ S 0 ∑ k ∈ V d x i j k d = 1 ,  ∀ j ∈ S 1 (9) ∑ j ∈ S ∑ d ∈ S 0 ∑ k ∈ V d x i j k d = 1 ,  ∀ i ∈ S 1 (10) ∑ j ∈ S 0 x i j k d = ∑ j ∈ S 0 x j i k d = 0 ,  ∀ i ∈ S 0 ,   ∀ d ∈ S 0 ,   ∀ k ∈ V d (11) x i j k d = 1,0 ,  ∀ i ∈ S ,   ∀ j ∈ S ,   ∀ d ∈ S 0 ,   ∀ k ∈ V dwhere 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 thresholdrmax set by decision makers. Constraint (6) expresses that the transportation risk of each route must be less than or equal to thresholdRmax set by decision makers. Constraint (7) indicates 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 that 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 xijkd, if the route of Hazmat transport vehiclek from distribution center d containing the segment from nodeito nodej, xijkd =1; else xijkd =0.


### 2.2.4. 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 setXvrp 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].Objective function is(12)Rr=Γ^r^l+min⁡∑m=1n2rmxm+∑m=1lr^m-r^lxm(13)Tt=Γ^t^l+min⁡∑m=1n2tmxm+∑m=1lt^m-t^lxmConstraint condition is(14)x∈XvrpThen, the optimal objective function value can be obtained asR∗=minl=1,⋯n2+1⁡R(l) and T∗=minl=1,⋯n2+1⁡T(l).


## 2.2.1. Assumption

There are a few assumptions in this study:(1) Multiple hazmat distribution centers are existent(2) The supply of hazmat distribution centers is adequate(3) Vehicle loading capacity is provided and the demand of each customer is specified(4) Multiple vehicles of the distribution center can service the customer(5) The transportation risk and transportation time are identified among the customer demand points, but they are uncertain number as interval number


## 2.2.2. Symbol Definition

S 0: set of the Hazmat distribution centers, where S0={i∣i=1,2,…,m} shows that the number distribution center is m and the sequence number of nodes set is 1…m.S 1: set of customer demand points, where S1={i∣i=1,…,n} shows that the number of customer demand points is n and the sequence number of nodes set is 1…n.S: all nodes set in the transportation network, where S=S0∪S1.V d: available transportation vehicle set in the hazmat distribution center, where V={k∣k=1,2,…K},d∈S0.E: road section set among nodes.q i: demand of customer demand point i.L k d: maximum load of transport vehiclekfrom hazmat distribution centerd.r ~ i j: variable transport risk from customer demand pointsi andj, where r~ij∈[rij,rij+r⌢ij](r⌢ij≥0).r i j: transportation risk nominal value from customer demand pointsi andj.r ⌢ i j: deviation of the variable transport risk to its nominal value from customer demand pointsi andj, where r⌢ij≥0.t i j: travel time nominal value from customer demand pointsi andj.t ⌢ i j: deviation of variable travel time to its nominal value from customer demand points i and j, where t⌢ij≥0.t ~ i j: variable transport risk from customer demand pointsi andj, where t~ij∈[tij,tij+t⌢ij](t⌢ij≥0).J i r: the set of columns which all uncertain data r~ij belonging to the ith row of the variable risk matrix are in, here Jir≤n.Γ i r: parameter Γir∈[0,Jir] to adjust robust risk of robust discrete optimization method and control the risk degree of conservatism, where decimal is permitted.Γ i r: maximum integer less than Γir.ψir: the set of column subscriptsjof uncertain data r~ij of line i in the variable risk matrix r~ij.J i t: the set of columns with all uncertain data t~ij belonging to theith row of the variable time matrix, where Jit≤n.Γit: parameter Γit∈[0,Jit] to adjust robust time of robust discrete optimization method and control the time degree of conservatism, where decimal is permitted.Γ i t: the maximum integer less than Γit.ψit: the set of column subscriptj of uncertain data r~ij of linei in variable time matrix t~ij.


## 2.2.3. Multiobjective Route Robust Optimization Model

(1) m i n  Z 1 = ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ j ∈ S 1 r i j x i j k d + m a x Ψ i r ∪ m r ∣ Ψ i r ⊆ J i r , Ψ i r = Γ i r , m r ∈ J i r ∖ Ψ i r ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ j ∈ Ψ i r r ⌢ i j x i j k d + ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ m r ∈ J i r ∖ Ψ i r Γ i - Γ i r ⌢ i m r x i m r k d - ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S r i 0 x i 0 k d (2) m i n  Z 2 = ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ j ∈ S t i j x i j k d + m a x Ψ i t ∪ m t ∣ Ψ i t ⊆ J i , Ψ i t = Γ i t , m t ∈ J i t \ Ψ i t ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ j ∈ Ψ i t t ⌢ i j x i j k d + ∑ d ∈ S 0 ∑ k ∈ V d ∑ i ∈ S ∑ m t ∈ J i t \ Ψ i t Γ i - Γ i t ⌢ i m t x i m t k d (3) s.t.  ∑ i ∈ S ∑ j ∈ S 1 ∑ d ∈ S 0 x i j k d q j ≤ L k d ,  ∀ k ∈ V d (4) ∑ j ∈ S 1 ∑ k ∈ V d x i j k d ≤ V d ,  ∀ i ∈ S 0 ,   ∀ d ∈ S 0 (5) m a x i , j ∈ E r i j x i j k d ≤ r m a x ,  ∀ i ∈ S ,   ∀ j ∈ S ,   ∀ d ∈ S 0 ,   ∀ k ∈ V d (6) m a x i , j ∈ E ∑ i ∈ S 0 ∪ S 1 ∑ j ∈ S 0 ∪ S 1 r i j x i j k d - ∑ i ∈ S 1 ∑ j ∈ S 0 r i j x i j k d ≤ R m a x ,  ∀ i ∈ S ,   ∀ j ∈ S ,   ∀ d ∈ S 0 ,   ∀ k ∈ V d (7) ∑ j ∈ S 1 x i j k d - ∑ j ∈ S 1 x j i k d = 0 ,  ∀ i ∈ S 0 ,   ∀ d ∈ S 0 ,   ∀ k ∈ V d (8) ∑ i ∈ S ∑ d ∈ S 0 ∑ k ∈ V d x i j k d = 1 ,  ∀ j ∈ S 1 (9) ∑ j ∈ S ∑ d ∈ S 0 ∑ k ∈ V d x i j k d = 1 ,  ∀ i ∈ S 1 (10) ∑ j ∈ S 0 x i j k d = ∑ j ∈ S 0 x j i k d = 0 ,  ∀ i ∈ S 0 ,   ∀ d ∈ S 0 ,   ∀ k ∈ V d (11) x i j k d = 1,0 ,  ∀ i ∈ S ,   ∀ j ∈ S ,   ∀ d ∈ S 0 ,   ∀ k ∈ V dwhere 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 thresholdrmax set by decision makers. Constraint (6) expresses that the transportation risk of each route must be less than or equal to thresholdRmax set by decision makers. Constraint (7) indicates 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 that 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 xijkd, if the route of Hazmat transport vehiclek from distribution center d containing the segment from nodeito nodej, xijkd =1; else xijkd =0.


## 2.2.4. 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 setXvrp 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].Objective function is(12)Rr=Γ^r^l+min⁡∑m=1n2rmxm+∑m=1lr^m-r^lxm(13)Tt=Γ^t^l+min⁡∑m=1n2tmxm+∑m=1lt^m-t^lxmConstraint condition is(14)x∈XvrpThen, the optimal objective function value can be obtained asR∗=minl=1,⋯n2+1⁡R(l) and T∗=minl=1,⋯n2+1⁡T(l).


## 3. 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–30].


### 3.1. Fuzzy C-Means Clustering

Suppose that n data samples areX={X1,X2,⋯,Xn},C(2≤C≤n) is the number of types into which the data samples are to be divided, {A1,A2,⋯,Ac} indicating the corresponding C categories,U is its similar classification matrix, the clustering centers of all categories are {V1,V2,⋯,Vc}, and μk(Xi) is the membership degree of sample Xi to categoryAk (abbreviated as μik). Then the objective function can be expressed as follows:(15)JbU,V=∑i=1n∑k=1cμikbdik2where dik=d(Xi-Vk)=∑j=1mXij-Vkj2, it is the synthetic weighted value of the transport risk and time after nondimensionalization between thei-th sampleXi and thek-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 Jb. It requires a sample for the sum of the membership degree of each cluster is 1, which is satisfied:(16)∑j=1cμjXi=1,i=1,2,…,nFormulas (17) and (18) are used to calculate separately the membership degree μik of the sampleXi for the categoryAk and C clustering centers {Vi}:(17)μik=1∑j=1cdik/djk2/b-1LetIk={i∣2≤C<n;dik=0}, for all the i categories, i∈Ik, μik=0.(18)Vij=∑k=1nμikbXkj∑k=1nμikbUsing 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.


### 3.2. 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.(  1) 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.(  2) 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.(  3) 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.Figure 1
Crossover operation.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 Figure2.Figure 2
Adjustment operation.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.(  4) 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 Figure3.Figure 3
Mutation operation.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].(  5) 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 1.
Create and initialize a group so that theex_archives external file is empty.Step 2.
Evaluate all particles and add the noninferior solution to the external file.Step 3.
Maintain external files according to the adaptive grid method.Step 4.
Selectgbest andpbest for each particle.Step 5.
Updating the velocity and position of the particles according to the speed formula and the position formula of the particle swarm.Step 6.
Make sure the particles exist in the search space.Step 7.
If the termination condition is satisfied, the output result algorithm is terminated; if it is not satisfied, Step2 is executed and the execution is continued.


## 3.1. Fuzzy C-Means Clustering

Suppose that n data samples areX={X1,X2,⋯,Xn},C(2≤C≤n) is the number of types into which the data samples are to be divided, {A1,A2,⋯,Ac} indicating the corresponding C categories,U is its similar classification matrix, the clustering centers of all categories are {V1,V2,⋯,Vc}, and μk(Xi) is the membership degree of sample Xi to categoryAk (abbreviated as μik). Then the objective function can be expressed as follows:(15)JbU,V=∑i=1n∑k=1cμikbdik2where dik=d(Xi-Vk)=∑j=1mXij-Vkj2, it is the synthetic weighted value of the transport risk and time after nondimensionalization between thei-th sampleXi and thek-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 Jb. It requires a sample for the sum of the membership degree of each cluster is 1, which is satisfied:(16)∑j=1cμjXi=1,i=1,2,…,nFormulas (17) and (18) are used to calculate separately the membership degree μik of the sampleXi for the categoryAk and C clustering centers {Vi}:(17)μik=1∑j=1cdik/djk2/b-1LetIk={i∣2≤C<n;dik=0}, for all the i categories, i∈Ik, μik=0.(18)Vij=∑k=1nμikbXkj∑k=1nμikbUsing 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.


## 3.2. 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.(  1) 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.(  2) 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.(  3) 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.Figure 1
Crossover operation.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 Figure2.Figure 2
Adjustment operation.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.(  4) 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 Figure3.Figure 3
Mutation operation.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].(  5) 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 1.
Create and initialize a group so that theex_archives external file is empty.Step 2.
Evaluate all particles and add the noninferior solution to the external file.Step 3.
Maintain external files according to the adaptive grid method.Step 4.
Selectgbest andpbest for each particle.Step 5.
Updating the velocity and position of the particles according to the speed formula and the position formula of the particle swarm.Step 6.
Make sure the particles exist in the search space.Step 7.
If the termination condition is satisfied, the output result algorithm is terminated; if it is not satisfied, Step2 is executed and the execution is continued.


## 4. 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 Table1, 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^m(0≤r^m<0.5rm) and transportation time deviation t^m(0≤t^m<0.5rm) are provided.Table 1
Customer demands.
Demandpoint 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20  Demand(ton) 2 1.5 4.5 3 1.5 4 2.5 3 3 4.5 2 1.5 2.5 3.5 2 2 2.5 3.5 2.5 3Table 2
The nominal transportation risk value of hazmat.
r i j a b c 1 2 3 …… 18 19 20 a 0   39 71 59 …… 38 80 52 b  0  30 30 49 …… 74 67 78 c   0 77 57 37 …… 49 34 64 1 39 30 77 0 32 30 …… 35 60 49 2 71 30 57 32 0 67 …… 33 65 55 3 59 49 37 30 67 0 …… 39 39 71 …… …… …… …… …… …… …… …… …… …… …… 18 38 74 49 35 33 39 …… 0 55 34 19 80 67 34 60 65 39 …… 55 0 67 20 52 78 64 49 55 71 …… 34 67 0Table 3
The nominal transportation time value of hazamt.
r i j a b c 1 2 3 …… 18 19 20 a 0   92 103 43 …… 107 97 44 b  0  37 36 88 …… 95 57 115 c   0 41 35 101 …… 40 74 64 1 92 37 41 0 75 102 …… 64 101 80 2 103 36 35 75 0 62 …… 108 84 95 3 43 88 101 102 62 0 …… 103 77 85 …… …… …… …… …… …… …… …… …… …… …… 18 107 95 40 64 108 103 …… 0 52 43 19 97 57 74 101 84 77 …… 52 0 86 20 44 115 64 80 95 85 …… 43 86 0We use FCMC algorithm to calculate the demand points clustering results, and the results can be seen in Table4.Table 4
Demand points clustering result.
Distribution center Demand points a 3,9,14,15,16,20 b 1,2,6,8,10,11,12 c 4,5,7,13,17,18,19Based 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 Tables5–13 and Figures 4–6.Table 5
Pareto solution set of robust control parametersΓ=0 for distribution center a.
Encoding Decoding Total risk Total time 9-20-14-3-15-16 a-9-20-a-14-3-a-15-16 237 536 9-16-14-3-15-20 a-9-16-a-14-3-a-15-20 261 448 15-20-9-14-3-16 a-15-20-9-a-14-3-a-16 247 520 14-3-16-9-20-15 a-14-3-a-16-9-20-a-15 253 469 15-14-16-9-3-20 a-15-14-16-a-9-3-a-20 291 404 14-3-15-9-16-20 a-14-3-a-15-9-16-a-20 285 425Table 6
Pareto solution set of robust control parametersΓ=10 for distribution center a.
Encoding Decoding Total risk Total time 3-14-20-9-16-15 a-3-14-a-20-9-16-a-15 317.817 637.81 9-20-16-14-3-15 a-9-20-16-a-14-3-a-15 319.893 621.749 9-20-14-16-3-15 a-9-20-a-14-16-a-3-15 377.505 520.018 20-9-16-14-3-15 a-20-9-16-a-14-3-a-15 337.903 532.141 9-3-15-14-16-20 a-9-3-a-15-14-16-a-20 385.542 463.803 9-3-20-15-14-16 a-9-3-a-20-15-a-14-16 354.151 527.306 9-20-14-3-15-16 a-9-20-a-14-3-a-15-16 305.49 657.577 20-15-14-3-9-16 a-20-15-a-14-3-a-9-16 359.107 525.349Table 7
Pareto solution set of robust control parametersΓ=20 for distribution center a.
Encoding Decoding Total risk Total time 3-14-20-9-16-15 a-3-14-a-20-9-16-a-15 323.453 761.423 9-3-15-20-14-16 a-9-3-a-15-20-a-14-16 381.844 560.087 9-20-16-14-3-15 a-9-20-16-a-14-3-a-15 357.98 644 9-16-20-14-3-15 a-9-16-20-a-14-3-a-15 391.378 558.043 9-20-14-16-3-15 a-9-20-a-14-16-a-3-15 415.592 542.269 9-3-15-14-16-20 a-9-3-a-15-14-16-a-20 423.629 486.054 9-20-15-14-3-16 a-9-20-15-a-14-3-a-16 378.69 606.48 3-14-9-20-16-15 a-3-14-a-9-20-16-a-15 316.995 866.785 9-3-20-15-14-16 a-9-3-a-20-15-a-14-16 392.239 549.557 3-14-16-9-20-15 a-3-14-a-16-9-20-a-15 336.135 750.943 9-20-14-3-15-16 a-9-20-a-14-3-a-15-16 343.577 679.828Table 8
Pareto solution set of robust control parametersΓ=0 for distribution center b.
Encoding Decoding Total risk Total time 12-10-11-1-2-6-8 b-12-10-11-b-1-2-6-b-8 289 582 11-10-1-12-6-2-8 b-11-10-b-1-12-6-b-2-8 303 471 12-6-10-11-1-2-8 b-12-6-b-10-11-b-1-2-8 302 530 12-10-11-1-6-8-2 b-12-10-11-b-1-6-b-8-2 293 573 1-12-10-2-6-8-11 b-1-12-10-b-2-6-b-8-11 296 572 8-6-12-10-11-1-2 b-8-6-b-12-10-11-b-1-2 299 549 2-1-6-12-10-11-8 b-2-1-6-b-12-10-11-b-8 284 627 8-6-1-12-2-10-11 b-8-6-b-1-12-2-b-10-11 363 458Table 9
Pareto solution set of robust control parametersΓ=10 for distribution center b.
Encoding Decoding Total risk Total time 2-6-8-12-1-11-10 b-2-6-b-8-12-1-b-11-10 415.472 552.167 1-6-2-11-10-12-8 b-1-6-2-b-11-10-12-b-8 348.014 799.695 2-8-6-1-12-11-10 b-2-8-b-6-1-12-b-11-10 399.364 631.219 1-12-6-11-10-8-2 b-1-12-6-b-11-10-b-8-2 375.264 660.062 2-12-1-6-8-11-10 b-2-12-1-b-6-8-b-11-10 485.04 533.589 2-8-6-12-1-11-10 b-2-8-b-6-12-1-b-11-10 402.203 557.881 1-6-2-12-10-8-11 b-1-6-2-b-12-10-b-8-11 346.521 872.417 1-6-12-11-10-8-2 b-1-6-12-b-11-10-b-8-2 356.381 724.274 2-6-1-11-10-12-8 b-2-6-1-b-11-10-12-b-8 381.268 645.953Table 10
Pareto solution set of robust control parametersΓ=20 for distribution center b.
Encoding Decoding Total risk Total time 2-6-1-12-10-8-11 b-2-6-1-b-12-10-b-8-11 417.863 778.297 1-6-2-12-10-8-11 b-1-6-2-b-12-10-b-8-11 352.157 1028.237 2-6-1-11-10-12-8 b-2-6-1-b-11-10-12-b-8 419.356 677.732 1-6-2-11-10-12-8 b-1-6-2-b-11-10-12-b-8 353.65 955.516 2-8-6-1-12-11-10 b-2-8-b-6-1-12-b-11-10 437.452 662.997 2-8-6-12-1-11-10 b-2-8-b-6-12-1-b-11-10 440.291 589.659 2-6-8-12-1-11-10 b-2-6-b-8-12-1-b-11-10 453.559 583.946 1-6-12-11-10-8-2 b-1-6-12-b-11-10-b-8-2 362.017 880.094 1-12-6-11-10-8-2 b-1-12-6-b-11-10-b-8-2 380.9 815.883 2-12-1-6-8-10-11 b-2-12-1-b-6-8-b-10-11 532.501 576.115Table 11
Pareto solution set of robust control parametersΓ=0 for distribution center c.
Encoding Decoding Total risk Total time 17-4-19-7-13-18-5 c-17-4-19-c-7-13-c-18-5 278 678 7-5-18-17-4-19-13 c-7-5-18-c-17-4-19-c-13 294 569 19-4-13-7-5-18-17 c-19-4-13-c-7-5-18-c-17 298 566 7-4-13-17-19-5-18 c-7-4-13-c-17-19-5-c-18 382 524 7-5-18-17-19-4-13 c-7-5-18-c-17-19-4-c-13 305 545 17-4-13-19-5-7-18 c-17-4-13-c-19-5-7-c-18 361 539 17-18-19-4-13-7-5 c-17-18-c-19-4-13-c-7-5 320 544 18-5-19-17-4-13-7 c-18-5-19-c-17-4-13-c-7 359 543 19-5-7-17-18-13-4 c-19-5-7-c-17-18-c-13-4 379 530Table 12
Pareto solution set of robust control parametersΓ=10 for distribution center c.
Encoding Decoding Total risk Total time 18-17-13-4-19-7-5 c-18-17-c-13-4-19-c-7-5 433.965 640.855 18-5-7-17-4-19-13 c-18-5-7-c-17-4-19-c-13 374.262 714.495 4-19-13-7-5-18-17 c-4-19-13-c-7-5-18-c-17 328.774 1068.644 18-17-19-4-13-7-5 c-18-17-c-19-4-13-c-7-5 432.616 649.645 18-17-13-4-19-5-7 c-18-17-c-13-4-19-c-5-7 442.263 632.41 5-18-7-17-4-19-13 c-5-18-7-c-17-4-19-c-13 361.711 813.376 18-17-7-5-19-13-4 c-18-17-c-7-5-19-c-13-4 565.935 629.522 18-5-13-17-4-19-7 c-18-5-13-c-17-4-19-c-7 369.796 741.803 13-4-18-17-19-5-7 c-13-4-c-18-17-c-19-5-7 538.139 632.006 18-17-13-19-4-7-5 c-18-17-c-13-19-4-c-7-5 415.204 661.122 4-19-13-17-18-5-7 c-4-19-13-c-17-18-5-c-7 345.113 937.323 18-17-13-19-4-5-7 c-18-17-c-13-19-4-c-5-7 421.737 652.677 4-19-13-17-5-18-7 c-4-19-13-c-17-5-18-c-7 333.451 975.738 4-19-13-17-18-7-5 c-4-19-13-c-17-18-c-7-5 351.435 928.3 18-5-7-19-4-17-13 c-18-5-7-c-19-4-17-c-13 388.454 682.357Table 13
Pareto solution set of robust control parametersΓ=20 for distribution center c.
Encoding Decoding Total risk Total time 18-5-13-17-4-19-7 c-18-5-13-c-17-4-19-c-7 418.386 797.49 13-4-19-17-18-7-5 c-13-4-19-c-17-18-c-7-5 565.673 677.462 13-4-19-17-18-5-7 c-13-4-19-c-17-18-5-c-7 558.373 686.485 13-4-19-18-17-7-5 c-13-4-19-c-18-17-c-7-5 588.331 674.181 18-17-13-19-4-7-5 c-18-17-c-13-19-4-c-7-5 463.794 716.809 4-19-13-17-5-18-7 c-4-19-13-c-17-5-18-c-7 339.087 1301.559 18-17-19-4-13-7-5 c-18-17-c-19-4-13-c-7-5 481.205 705.332 4-19-13-17-18-5-7 c-4-19-13-c-17-18-5-c-7 350.749 1263.143 18-17-13-4-19-7-5 c-18-17-c-13-4-19-c-7-5 482.554 696.543 18-17-13-19-4-5-7 c-18-17-c-13-19-4-c-5-7 470.326 708.364 4-19-13-18-17-7-5 c-4-19-13-c-18-17-c-7-5 379.729 1250.84 4-19-13-17-18-7-5 c-4-19-13-c-17-18-c-7-5 357.071 1254.121 4-19-13-7-5-18-17 c-4-19-13-c-7-5-18-c-17 334.41 1394.464 4-17-19-18-5-7-13 c-4-17-19-c-18-5-7-c-13 405.101 1245.591 18-5-7-17-19-4-13 c-18-5-7-c-17-19-4-c-13 449.149 738.192 18-17-13-4-19-5-7 c-18-17-c-13-4-19-c-5-7 490.853 688.097 5-18-7-17-4-19-13 c-5-18-7-c-17-4-19-c-13 409.798 895.155 18-5-7-17-4-19-13 c-18-5-7-c-17-4-19-c-13 422.852 770.182Figure 4
Pareto optimal solution distribution of robust control parametersΓ=0, Γ=10, and Γ=20 for distribution center a.Figure 5
Pareto optimal solution distribution of robust control parametersΓ=0, Γ=10, and Γ=20 for distribution center b.Figure 6
Pareto optimal solution distribution of robust control parametersΓ=0, Γ=10, and Γ=20 for distribution center c.In Table5, 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 Pareto solutions sets for distribution center c when the robustness control parameters are taken 0, 10, and 20. The above tables are analyzed similar to Table 5. Figures 4, 5, and 6, respectively, show the Pareto optimal solutions for distribution centers a, b, and c when the robustness control parameters Γ = 0, Γ = 10, and Γ = 20.When the robustness control parametersΓ = 0, Γ = 10, and Γ = 20, each objective optimal distribution route for all vehicles from distribution centers a, b, c is shown in Tables 14, 15, and 16. Aiming at distribution center a, when the robustness control parameters Γ = 0 and Γ = 10, the total risk 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.Table 14
Pareto optimal solutions for distribution center a whenΓ = 0, Γ = 10 and Γ = 20.
Γ Distribution center a Risk optimal route Time optimal route 0 a-9-20-a a-15-14-16-a a-14-3-a a-9-3-a a-15-16-a a-20-a  10 a-9-20-a a-15-14-16-a a-14-3-a a-9-3-a a-15-16-a a-20-a  20 a-9-20-16-a a-15-14-16-a a-14-3-a a-9-3-a a-15-a a-20-aTable 15
Pareto optimal solutions for distribution center b whenΓ = 0, Γ = 10 and Γ = 20.
Γ Distribution center b Risk optimal route Time optimal route 0 b-2-1-6-b b-8-6-b b-12-10-11-b b-1-12-2-b b-8-b b-10-11-b  10 b-1-6-2-b b-8-6-b b-12-10-b b-2-12-1-b b-8-11-b b-10-11-b  20 b-1-6-2-b b-8-6-b b-12-10-b b-2-12-1-b b-8-11-b b-10-11-bTable 16
Pareto optimal solutions for distribution center c whenΓ = 0, Γ = 10 and Γ = 20.
Γ Distribution center c Risk optimal route Time optimal route 0 c-17-4-19-c c-7-4-13-c c-7-13-c c-17-19-5-c c-18-5-c c-18-c  10 c-4-19-13-c c-13-4-c c-7-5-18-c c-7-5-19-c c-17-c c-18-17-c  20 c-4-19-13-c c-13-4-19-c c-7-5-18-c c-7-5-c c-17-c c-18-17-cFor 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 Table17. 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.Table 17
Performance comparison between FCMC-PSO algorithm and SPEA.
Index FCMC-PSO SPEA Γ value 0 20 40 0 20 40 Convergence Iterations 33 39 51 58 72 97 Run time/(s) 21 33 49 40 52 88


## 5. 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.



---

*Source: 2916391-2018-10-09.xml*