# Multidepot Heterogeneous Vehicle Routing Problem for a Variety of Hazardous Materials with Risk Analysis

Authors: Bochen Wang; Qiyuan Qian; Zheyi Tan; Peng Zhang; Aizhi Wu; Yi Zhou
Journal: Scientific Programming
 (2020)

---

## Abstract

This study investigates a multidepot heterogeneous vehicle routing problem for a variety of hazardous materials with risk analysis, which is a practical problem in the actual industrial field. The objective of the problem is to design a series of routes that minimize the total cost composed of transportation cost, risk cost, and overtime work cost. Comprehensive consideration of factors such as transportation costs, multiple depots, heterogeneous vehicles, risks, and multiple accident scenarios is involved in our study. The problem is defined as a mixed integer programming model. A bidirectional tuning heuristic algorithm and particle swarm optimization algorithm are developed to solve the problem of different scales of instances. Computational results are competitive such that our algorithm can obtain effective results in small-scale instances and show great efficiency in large-scale instances with 70 customers, 30 vehicles, and 3 types of hazardous materials.


---

## Body

## 1. Introduction

In the past few years, the rapid development of China’s chemical industry has led to the expansion of the road transportation market for hazardous materials. Thus, the transportation demand of hazardous materials presents a high-speed rising trend. Road transportation is still the main mode of transportation for hazardous materials, with the advantage of its strong mobility, flexibility, continuous transportation capacity, and less restrictive requirements. In contrast, the disadvantages of road transportation are the greatest potential risk, high unit transportation costs, and variability in operating conditions. Due to the characteristics of being flammable, explosive, corrosive, toxic, and radioactive, dangerous goods often cause more serious secondary hazards in traffic accidents, not only the serious loss of personnel and property, but also a huge negative impact on society. Enterprises need to consider the cost of transportation and, more importantly, to ensure the safety of the transportation process, which has brought great challenges. In order to reduce transportation costs, companies will choose shorter transportation routes, and these transportation routes will pass through areas with higher population density, leading to increased potential risks in the transportation process, while choosing a farther transportation route will undoubtedly increase total costs. As reported in this research, the problem of hazardous materials transportation includes two basic objective, minimum cost and risk control [1]. The key to solving such problems is to make a balance between the two points.According to the regulations of the Ministry of Transport of the People’s Republic of China, different types of hazardous materials transportation vehicles have different restrictions during actual operation. For some flammable and explosive products, according to the regulations, transport on specific roads is not allowed without permission. However, limited quantities of dangerous goods with a total mass not exceeding 8000 kg can be transported as ordinary goods. In addition, the hazardous materials manufacturers have the characteristics of limited production and unsuitable long-term storage in the actual situation. During the transportation process, a single visit to a customer cannot take advantage of the vehicle transportation, so as to control the potential risks. It is reasonable for the carrier to complete the transportation task by traversing multiple manufacturers of similar products. This is also the dilemma currently faced by hazardous materials transportation companies. In this paper, we propose a multidepot heterogeneous vehicle routing model in order to find the optimal scheme to complete the transportation task among different manufacturers.The traditional vehicle routing problem (VRP) aims generally to make the transportation cost the lowest under certain restrictions. It is based on general cargo transportation and does not consider the vehicle type, which means that all the vehicles can carry different types of goods. However, in heterogeneous vehicles routing problem, different types of goods must be transported by the specified vehicle type. In hazardous materials transportation, gasoline requires tanker trucks while solid hazardous materials should be carried by special lorry, and they are not interchangeable. In the real business, transportation always takes place between multiple depots and customer points. In addition, risk needs to be considered in the road transportation of hazardous materials. The treatment of risk in this article aims to convert the hazard after the accident into a risk coefficient and introduce uncertain scenarios to predict the probability of risk. The product of the probability of risk and the risk coefficient is the risk cost, which is used to add weighted transportation costs to obtain the total objective. The risks include road accidents, leakage of dangerous goods, and pedestrians on the way. Comprehensive consideration of factors such as transportation costs, multiple depots, heterogeneous vehicles, risks, and multiple accident scenarios can make the study closer to the actual situation. However, in contrast, the more factors we considered, the more difficult it is to solve the problem. Most of the research involves only two or three factors, and few involve more than three factors, which also makes our research valuable.We propose a more complicated mixed integer programming (MIP) model which is close to the actual situation. Existing solvers, such as CPLEX, can only solve small-scale calculation examples. As the scale of the problem becomes larger, the efficiency of CPLEX’s solution continues to decrease. Therefore, we design a bidirectional tuning heuristic and a particle swarm optimization algorithm (PSO) to solve problems at different scales. The results show that the algorithms we designed can effectively solve larger-scale problems and can also find near-accurate solutions to small-scale problems.The reminder of this paper is organized as follows. Section2 reviews related literature. Section 3 formulates a mathematical model. In Section 4, two different algorithms are proposed. Section 5 reports the numerical experiments. Conclusions are summarized in the last section.


## 2. Related Literature

The difference between vehicle routing problem for hazardous materials (VRPHM) and traditional VRP is that VRPHM needs to consider the risk and it is a key issue in research which cannot be ignored. Bula et al. [2] pointed out that the research on VRPHM mainly focuses on the shortest path problem and the routing problem. The former research focuses on controlling risks, and such research has already achieved many results. A population exposure model was proposed by [3], which defines the range of the affected population on each side, and the risk value on each side is summed up to obtain the total risk value on the entire path. Tarantilis and Kiranoudis [4] had also done a similar topic research; the difference was that this research is based on vehicle fleet to meet customer needs, rather than focusing on the shortest path, which is more realistic. Bula et al. [5] defined the risk value of hazardous materials during transportation. The risk value is composed of the probability of accident, the extent of the accident, and the population density of the route. Based on this, a MIP model with the minimum risk value as the objective was established. Cuneo et al. [6] focused on the transportation of fuel oil; they proposed an innovative risk index in the function as far as possible to reduce the risk of accidents during the transportation of fuel oil. The method of most related work is to calculate fixed risk values few papers take stochastic into consideration. However, in other fields, the stochastic is widely used to access uncertainty. Rabbani et al. [7] studied hazardous waste management and controlled the processing costs of industrial wastes by introducing scenarios to defined risks. The risk defined in this paper is much like the method used in [8].Another direction of research is to solve the VRPHM from the perspective of traditional VRP. Pradhananga et al. [9] considered the multiobjective problem of transportation time and risk. Risk in their paper was obtained by multiplying the incidence rate of each arc accident by the population density of the region. Du et al. [10] constructed fuzzy bilevel programming, which is a multidepot problem, and proposed a fuzzy simulation-based heuristic algorithm to solve the problem.To the best of our knowledge, VRPHM with heterogeneous vehicle and uncertainty risks is rarely studied. Homogeneous models are easy to study, but the transportation of hazardous materials usually involves the problem of heterogeneous vehicle. Golden et al. [11] considered the problem of distribution of mixed vehicles with unlimited number, but the article only focused on the capacity of the vehicles. In the general VRP problem, there are many studies on heterogeneous vehicles. A number of effective algorithms were proposed to solve this problem; [12–14] used exact algorithms while [15, 16] designed heuristic algorithms. Due to the fact that VRPHM needs to consider more factors than general VRP, few scholars have studied heterogeneous vehicles.In sum, the contributions of this paper are as follows. First, we propose a complex model to solve multidepot VRPHM, considering risk, heterogeneous vehicles, transportation cost, and time limit. Second, we describe the risk as uncertainty in different scenarios. Third, two effective algorithms are developed.


## 3. Problem Definition and Formulation

In this section, we describe the multidepot VRPHM, define the parameters, and then formulate a MIP model.


### 3.1. Problem Definition

The objective of this problem is to minimize the transportation costs, the risk costs resulting from the transformation of transportation risks, and the penalty cost of additional working time. We also consider various sudden scenes during transportation and consider the distribution of hazardous materials in multiple depots and heterogeneous vehicles. Different vehicles have different capacities and distribution categories. The vehicle cost in this paper consists of fixed cost and variable cost. Fixed cost refers to the cost unrelated to transportation which includes depreciation cost and insurance cost, while variable cost is related to the transportation distance and chosen route, which includes transportation cost and risk cost. The fixed costs of starting a vehicle and the variable costs per unit of travel distance are also different. The transportation network of hazardous materials can be defined on a directed networkG=O,A, where the set O represents the union of the set D composed of depots and the set N composed of customer points, and A=i,ji,j∈O is an edge set. The network is depicted clearly in Figure 1.Figure 1
A vehicle routing network with three depots and several customers.Let the setK represent the vehicles of fixed cost gk. The maximum load of each vehicle is qk. The vehicle determines whether it can visit the customer based on the type of hazardous materials, described using parameter ci,k. Each vehicle must depart from its own depot 0 to the customer i∈N. After process time ti, the vehicle can travel from customer i∈N to customer j∈N. We introduce a virtual node N+1 to represent the depot that the vehicle returns to after the transportation. This virtual node is the same node as 0. The transportation time is different according to the chosen route which can be expressed as di,j,r,k, where the route set is Ri,j. Besides, the travel time is different from 0 to each customer if the vehicle belongs to different depots. Each vehicle k has the latest working time yk, and an additional overtime costs ctime is needed if the total transportation time is beyond the latest working time. Accidents are very likely to occur in the transportation of hazardous materials. This paper introduces the concept of scenarios with different probability to eliminate such uncertainties. The set of scenarios is represented by S. The risk under different scenarios is related to the route chosen by the vehicle. We use wi,j,r,s to describe the probability of accident that may occur from the customer i∈N to the customer j∈N.Before formulating the model, some assumptions are listed as follows:(1)
All vehicles must return to the depots where they came from after transportation(2)
Each vehicle can only transport one type of hazardous materials(3)
Each customer can be accessed multiple times by different vehicles(4)
Each vehicle serves one customer no more than once(5)
Each vehicle must wait for the preceding vehicle to complete the operation before starting service


### 3.2. Notations

In this section, we define all the parameters and decision variables (Table1).Table 1
Notations.
Indices and setsi,jIndex of customer.NSet of customers.kIndex of vehicle.KSet of vehicles.sIndex of scenario.SSet of scenarios.rIndex of route.Ri,jSet of routes from customeri to customer j.vIndex of type of hazardous materials.VSet of type of hazardous materials.Parametersci,kEquals 1 if vehiclek can visit customer i.di,j,r,kTravel time from customeri to customer j for vehicle k on route r.tiService time at customeri.fi,vTransportation demand of hazardous materialsv at customer i.qkCapacity of vehiclek.ykLatest working time of vehiclek.gkFixed cost of vehiclek.wi,j,r,sRisk probability from customeri to customer j on route r.ctranCoefficient of transportation cost.criskCoefficient of risk cost.ctimeCoefficient of work overtime cost.Decision variablesαi,kBinary, equals 1 if vehiclek serves i.γi,j,k,r,sBinary, equals 1 if vehiclek travels from customer i to customer j by route r on scenario s.εi,k,k′,sBinary, equals 1 if vehiclek′ serves customer i earlier than vehicle k.ηkBinary, equals 1 if vehiclek is started.μk,vBinary, equals 1 if vehiclek transports hazardous material v.βi,k,sFloat, time point when vehiclek arrives at customer i on scenario s.δi,k,sFloat, time point when vehiclek starts serving customer i on scenario s.ζi,k,s,vInteger, quantity of hazardous materialv transported by vehicle k at customer i on scenario s.


### 3.3. Mathematical Model

(1)Minimize∑k∈Kgkηk+1S∑s∈S∑k∈K∑j∈N∪eN∑i∈N∪0∑r∈Ri,jctrandi,j,r,kγi,j,k,r,s+∑s∈S∑k∈K∑j∈N∪eN∑i∈N∪0∑r∈Ri,jcriskwi,j,r,sγi,j,k,r,s+∑s∈S∑k∈Kctimeβe,k,s−yk+,(2)subject to∑k∈Kαi,k≥1,∀i∈N,(3)αi,k≤ci,k,∀i∈N,k∈K,(4)αi,k≤ηk,∀i∈N,k∈K,(5)∑i∈N∪0∑r∈Ri,jγi,j,k,r,s=∑i∈N∪eN∑r∈Ri,jγj,i,k,r,s=αj,k,∀j∈N,k∈K,s∈S,(6)∑j∈N∪eN∑r∈R0,jγ0,j,k,r,s=∑j∈N∪0∑r∈R0,jγj,e,k,r,s=1,∀k∈K,s∈S,(7)∑i∈N∑j∈N∑r∈Ri,jγi,j,k,r,s≤N1−1,N1⊆N,2≤N1≤N,∀k∈K,r∈Ri,j,s∈S,(8)γi,j,k,r,s≤bi,j,r,s,∀i∈N∪0,j∈N∪eN,k∈K,r∈Ri,j,s∈S,(9)∑i∈Nζi,k,s,v≤qk,∀k∈K,s∈S,(10)∑k∈Kζi,k,s,v=fi,v,∀i∈N,s∈S,v∈V,(11)ζi,k,s,v≤fi,vμk,vαi,k,∀i∈N,k∈K,v∈V,s∈S,(12)∑v∈μk,v=ηk,1,∀k∈K,(13)δi,k,s+ti+di,j,r,k≤βj,k,s+M1−γi,j,k,r,s,∀i∈N∪0,j∈N∪eN,k∈K,r∈Ri,j,s∈S,(14)δi,k,s≥βi,k,s,∀i∈N,k∈K,s∈S,(15)δi,k,s+ti≤δi,k′,s+M1−εi,k,k′,s,∀i∈N,k∈K,s∈S,(16)∑k∈K∪0Kεi,k,eK,s=∑k∈K∪eKεi,0K,k,s=1,∀i∈N,k∈K,k∈K,(17)∑k∈K∪0Kεi,k,k′,s=∑k∈K∪eKεi,k′,k,s=1,∀i∈N,k′∈K,s∈S,(18)∑k∈K∑k′∈Kεi,k,k′,s≤K1−1,K1⊆K,2≤K1≤K,∀i∈N,s∈S,(19)αi,k,γi,j,k,r,s,εi,k,k′,s,ηk∈0,1,∀i∈N∪0,j∈N∪e,k∈K,r∈Ri,j,s∈S,(20)βi,k,s,εi,k,k′,s,ζi,k,s≥0,∀i,j∈N,k,k′∈K,s∈S.The objective function (1) minimizes the total costs of vehicle startup, additional overtime costs, and risk conversion costs. Constraints (2) ensure that a customer is visited at least once. Constraints (3) limit whether vehicle k can serve customer i. Constraints (4) show that if vehicle k serves customer i, the vehicle must be started. Constraints (5) make sure that if a customer is visited by a vehicle, other nodes must be visited before and after this node. Constraints (6) mean that, in any scenario, each vehicle must travel from the original depot and return to the same depot. Constraints (7) eliminate subpaths. Constraints (8) impose the condition that whether the vehicles are allowed to choose a route r from customer i to customer j under scenario s is determined by whether the route is passable. For example, some routes may be restricted during the morning and evening rush hours, specific time period, and temporary event. Constraints (9) guarantee that the total loads of a vehicle will not exceed the maximum vehicle capacity. Constraints (10) and (11) ensure that all the demands of hazardous materials must be satisfied. Constraints (12) put the limit that each vehicle can only transport one type of hazardous material. Constraints (13) and (14) describe the time relationship between the time point of arrival and time point of service. Constraints (15) determine the service time of different vehicles at the same customer according to their sequence. Constraints (16) and (17) address the sequence of vehicle access for every customer. Constraints (18) ensure that there are no subpaths in vehicle loop. Constraints (19) and (20) limit the range of decision variables.


### 3.4. Linearizing the Products of Two Variables in the Constraint

In constraint (11), the load quantity contains the form of products of variables μk,v and αi,k. The nonlinear part is not conducive for solver to calculate. To linearize the constraint, some new variables are added as follows. Constraints (11) are then transformed into constraints (21)–(24).A newly defined variable is as follows:θi,k,v: binary, equals 1 if αi,k=1 and μk,v=1.Newly defined constraints are as follows:(21)θi,k,v≤μk,v,∀i∈N,k∈K,v∈V,(22)θi,k,v≤αi,k,∀i∈N,k∈K,v∈V,(23)θi,k,v≥αi,k+μk,v−1,∀i∈N,k∈K,v∈V,(24)ζi,k,s,v≤fi,vθi,k,v,∀i∈N,k∈K,v∈V,s∈S.


## 3.1. Problem Definition

The objective of this problem is to minimize the transportation costs, the risk costs resulting from the transformation of transportation risks, and the penalty cost of additional working time. We also consider various sudden scenes during transportation and consider the distribution of hazardous materials in multiple depots and heterogeneous vehicles. Different vehicles have different capacities and distribution categories. The vehicle cost in this paper consists of fixed cost and variable cost. Fixed cost refers to the cost unrelated to transportation which includes depreciation cost and insurance cost, while variable cost is related to the transportation distance and chosen route, which includes transportation cost and risk cost. The fixed costs of starting a vehicle and the variable costs per unit of travel distance are also different. The transportation network of hazardous materials can be defined on a directed networkG=O,A, where the set O represents the union of the set D composed of depots and the set N composed of customer points, and A=i,ji,j∈O is an edge set. The network is depicted clearly in Figure 1.Figure 1
A vehicle routing network with three depots and several customers.Let the setK represent the vehicles of fixed cost gk. The maximum load of each vehicle is qk. The vehicle determines whether it can visit the customer based on the type of hazardous materials, described using parameter ci,k. Each vehicle must depart from its own depot 0 to the customer i∈N. After process time ti, the vehicle can travel from customer i∈N to customer j∈N. We introduce a virtual node N+1 to represent the depot that the vehicle returns to after the transportation. This virtual node is the same node as 0. The transportation time is different according to the chosen route which can be expressed as di,j,r,k, where the route set is Ri,j. Besides, the travel time is different from 0 to each customer if the vehicle belongs to different depots. Each vehicle k has the latest working time yk, and an additional overtime costs ctime is needed if the total transportation time is beyond the latest working time. Accidents are very likely to occur in the transportation of hazardous materials. This paper introduces the concept of scenarios with different probability to eliminate such uncertainties. The set of scenarios is represented by S. The risk under different scenarios is related to the route chosen by the vehicle. We use wi,j,r,s to describe the probability of accident that may occur from the customer i∈N to the customer j∈N.Before formulating the model, some assumptions are listed as follows:(1)
All vehicles must return to the depots where they came from after transportation(2)
Each vehicle can only transport one type of hazardous materials(3)
Each customer can be accessed multiple times by different vehicles(4)
Each vehicle serves one customer no more than once(5)
Each vehicle must wait for the preceding vehicle to complete the operation before starting service


## 3.2. Notations

In this section, we define all the parameters and decision variables (Table1).Table 1
Notations.
Indices and setsi,jIndex of customer.NSet of customers.kIndex of vehicle.KSet of vehicles.sIndex of scenario.SSet of scenarios.rIndex of route.Ri,jSet of routes from customeri to customer j.vIndex of type of hazardous materials.VSet of type of hazardous materials.Parametersci,kEquals 1 if vehiclek can visit customer i.di,j,r,kTravel time from customeri to customer j for vehicle k on route r.tiService time at customeri.fi,vTransportation demand of hazardous materialsv at customer i.qkCapacity of vehiclek.ykLatest working time of vehiclek.gkFixed cost of vehiclek.wi,j,r,sRisk probability from customeri to customer j on route r.ctranCoefficient of transportation cost.criskCoefficient of risk cost.ctimeCoefficient of work overtime cost.Decision variablesαi,kBinary, equals 1 if vehiclek serves i.γi,j,k,r,sBinary, equals 1 if vehiclek travels from customer i to customer j by route r on scenario s.εi,k,k′,sBinary, equals 1 if vehiclek′ serves customer i earlier than vehicle k.ηkBinary, equals 1 if vehiclek is started.μk,vBinary, equals 1 if vehiclek transports hazardous material v.βi,k,sFloat, time point when vehiclek arrives at customer i on scenario s.δi,k,sFloat, time point when vehiclek starts serving customer i on scenario s.ζi,k,s,vInteger, quantity of hazardous materialv transported by vehicle k at customer i on scenario s.


## 3.3. Mathematical Model

(1)Minimize∑k∈Kgkηk+1S∑s∈S∑k∈K∑j∈N∪eN∑i∈N∪0∑r∈Ri,jctrandi,j,r,kγi,j,k,r,s+∑s∈S∑k∈K∑j∈N∪eN∑i∈N∪0∑r∈Ri,jcriskwi,j,r,sγi,j,k,r,s+∑s∈S∑k∈Kctimeβe,k,s−yk+,(2)subject to∑k∈Kαi,k≥1,∀i∈N,(3)αi,k≤ci,k,∀i∈N,k∈K,(4)αi,k≤ηk,∀i∈N,k∈K,(5)∑i∈N∪0∑r∈Ri,jγi,j,k,r,s=∑i∈N∪eN∑r∈Ri,jγj,i,k,r,s=αj,k,∀j∈N,k∈K,s∈S,(6)∑j∈N∪eN∑r∈R0,jγ0,j,k,r,s=∑j∈N∪0∑r∈R0,jγj,e,k,r,s=1,∀k∈K,s∈S,(7)∑i∈N∑j∈N∑r∈Ri,jγi,j,k,r,s≤N1−1,N1⊆N,2≤N1≤N,∀k∈K,r∈Ri,j,s∈S,(8)γi,j,k,r,s≤bi,j,r,s,∀i∈N∪0,j∈N∪eN,k∈K,r∈Ri,j,s∈S,(9)∑i∈Nζi,k,s,v≤qk,∀k∈K,s∈S,(10)∑k∈Kζi,k,s,v=fi,v,∀i∈N,s∈S,v∈V,(11)ζi,k,s,v≤fi,vμk,vαi,k,∀i∈N,k∈K,v∈V,s∈S,(12)∑v∈μk,v=ηk,1,∀k∈K,(13)δi,k,s+ti+di,j,r,k≤βj,k,s+M1−γi,j,k,r,s,∀i∈N∪0,j∈N∪eN,k∈K,r∈Ri,j,s∈S,(14)δi,k,s≥βi,k,s,∀i∈N,k∈K,s∈S,(15)δi,k,s+ti≤δi,k′,s+M1−εi,k,k′,s,∀i∈N,k∈K,s∈S,(16)∑k∈K∪0Kεi,k,eK,s=∑k∈K∪eKεi,0K,k,s=1,∀i∈N,k∈K,k∈K,(17)∑k∈K∪0Kεi,k,k′,s=∑k∈K∪eKεi,k′,k,s=1,∀i∈N,k′∈K,s∈S,(18)∑k∈K∑k′∈Kεi,k,k′,s≤K1−1,K1⊆K,2≤K1≤K,∀i∈N,s∈S,(19)αi,k,γi,j,k,r,s,εi,k,k′,s,ηk∈0,1,∀i∈N∪0,j∈N∪e,k∈K,r∈Ri,j,s∈S,(20)βi,k,s,εi,k,k′,s,ζi,k,s≥0,∀i,j∈N,k,k′∈K,s∈S.The objective function (1) minimizes the total costs of vehicle startup, additional overtime costs, and risk conversion costs. Constraints (2) ensure that a customer is visited at least once. Constraints (3) limit whether vehicle k can serve customer i. Constraints (4) show that if vehicle k serves customer i, the vehicle must be started. Constraints (5) make sure that if a customer is visited by a vehicle, other nodes must be visited before and after this node. Constraints (6) mean that, in any scenario, each vehicle must travel from the original depot and return to the same depot. Constraints (7) eliminate subpaths. Constraints (8) impose the condition that whether the vehicles are allowed to choose a route r from customer i to customer j under scenario s is determined by whether the route is passable. For example, some routes may be restricted during the morning and evening rush hours, specific time period, and temporary event. Constraints (9) guarantee that the total loads of a vehicle will not exceed the maximum vehicle capacity. Constraints (10) and (11) ensure that all the demands of hazardous materials must be satisfied. Constraints (12) put the limit that each vehicle can only transport one type of hazardous material. Constraints (13) and (14) describe the time relationship between the time point of arrival and time point of service. Constraints (15) determine the service time of different vehicles at the same customer according to their sequence. Constraints (16) and (17) address the sequence of vehicle access for every customer. Constraints (18) ensure that there are no subpaths in vehicle loop. Constraints (19) and (20) limit the range of decision variables.


## 3.4. Linearizing the Products of Two Variables in the Constraint

In constraint (11), the load quantity contains the form of products of variables μk,v and αi,k. The nonlinear part is not conducive for solver to calculate. To linearize the constraint, some new variables are added as follows. Constraints (11) are then transformed into constraints (21)–(24).A newly defined variable is as follows:θi,k,v: binary, equals 1 if αi,k=1 and μk,v=1.Newly defined constraints are as follows:(21)θi,k,v≤μk,v,∀i∈N,k∈K,v∈V,(22)θi,k,v≤αi,k,∀i∈N,k∈K,v∈V,(23)θi,k,v≥αi,k+μk,v−1,∀i∈N,k∈K,v∈V,(24)ζi,k,s,v≤fi,vθi,k,v,∀i∈N,k∈K,v∈V,s∈S.


## 4. Algorithm Strategies

The model presented in Section3 can be solved by solvers (e.g., CPLEX) directly in small-scale examples. To deal with large-scale problems, we develop two different heuristics, bidirectional tuning heuristic and particle swarm optimization algorithm (PSO), to solve our proposed model. In Sections 4.1 and 4.2, we describe the framework of the two heuristics, respectively.


### 4.1. Bidirectional Tuning Heuristic

The main idea of bidirectional tuning heuristic is to transform the model into several interrelated subproblems. Solving the subproblems using some fixed decision variables, we can obtain the remaining decision variables and then use the them to gain other decision variables. In the iterative process, if the value of the objective is better than the optimal value, the optimal solution needs to be updated. The exit mechanism of the algorithm can be determined by the number of iterations or by the rules to judge whether the objective value can get a better solution.Before we start solving, the decision variables need to be classified. The decision variables in this paper can be divided into three different types according to their definitions. The first type isαi,k, the variables used to decide whether a vehicle visits a customer. The second type is μk,v and ηk, used to determine the type of hazardous material each vehicle transports. The other decision variables contain the dimensions of random scenarios and are limited by the first two kinds of decision variables, so they are always taken as the variables to be solved. We determine the customer that each vehicle visits and then optimize the transportation types with fixed αi,k. Then, we move to determining the transportation types and optimize the customers with fixed μk,v and ηk. The solving is repeated until the iteration reaches upper limit or no improvement can be obtained. The basic procedure can be defined as follows and the detailed process is shown in Algorithm 1.Step 1: solve the model to find a feasible solution, which can be taken as the initial solution. Record the initial decision variables and the objective value.Step 2: based onαi,k from initialization, determine decision variable ηk and μk,v. Compare the obtained objective value with the historical optimal value, and if the value is better than the historical optimal value, update the optimal value.Step 3: solve the model with fixedηk and μk,v to obtain αi,k. Update the objective value.Step 4: judge the exit condition; if it is satisfied, the optimal solution is obtained and the algorithm stops; otherwise, go to step 2.Algorithm 1: Bidirectional tuning heuristic.
(1)
Initialization(2)
Solve the model, and determineαi,k(3)
Update the optimal solution(4)
While (n < MaxIter)(5)
Solve model with fixedαi,k, determine μk,v and ηk(6)
Update the optimal solution(7)
Solve model with fixedμk,v and ηk, determine αi,k(8)
Update the optimal solution(9)
End while


### 4.2. Particle Swarm Optimization (PSO)

Particle swarm optimization (PSO) has been studied extensively since it was proposed [17]. As one of the swarm intelligence algorithms, PSO has been widely used in vehicle routing problem in the past few years, for example, vehicle routing problem with multiple pickup and delivery [18] and multidepot multitrip vehicle routing problem [19].


#### 4.2.1. Initialization and Velocity Updating Strategy

The particle swarm optimization algorithm not only considers the personal optimal value, but also records the global optimal value of the entire group in the process of searching for solutions. It has excellent performance in solving optimization problems. The movement of particles depends on the update of their position and velocity; each particle contains the current position and the personal optimal position. The current position of the particles is adjusted by updating the velocity formula, the solution is judged by the specific fitness value, and the optimal position of the entire group is updated in the global range.In the model we proposed, the decision variablesαi,k and μk,v are obviously related to other decision variables. When αi,k is fixed, ηk can be determined by constraints (4) and γi,j,k,r,s can also be determined by constraints (5)–(7). Furthermore, the decision variables βi,k,s, related to time, are limited by γi,j,k,r,s according to constraints (13). Then the variable δi,k,s can be derived from constraints (13). At this point, we can judge whether the service time δi,k,s of each customer is reasonable. If it is not reasonable, we update the time according to the time sequence of starting time and the processing time. It is also necessary to trace back to update variable βi,k,s, which is the subsequent vehicle arrival time of the route. When all decision variables related to time are determined, we can easily confirm the decision variable εi,k,k′,s. When we know both αi,k and μk,v, the load variable ζi,k,s,v can also be determined. We now have clarified the relationship between all the decision variables.Based on the description above, we setμk,v as the outer variable of the algorithm. And the particles that we set up contain information to find the vehicles assigned to customer and to generate visit sequence of each vehicle, which is similar to the approach taken in this work [20]. We define the particle as Pmn, and its velocity can be defined as Velmn. The velocity update formula can be expressed by formula (25) and the position update formula is (26). We use PBestmn to represent the best personal optimal position of the particle and GBestmn to represent the global optimal position, where m means the number of particles (we take 300) and n means the current iteration number:(25)Velmn+1=wnVelmn+c1r1PBestmn−Pmn+c2r2GBestmn−Pmn,(26)Pmn+1=Pmn+Velmn.In the formula for particle velocity update,wn denotes inertia weight and it is calculated by equation wn=wmax−wmax−wminn/N, in which N is the number of total iterations (we take 50) and wmax and wmin are maximum and minimum inertia weight (we take 0.9 and 0.4). This makes the particle swarm have strong global convergence ability at the beginning, but strengthens local convergence ability with the increase of iteration. c1 and c2 are acceleration weights (we take 0.683 for both), while r1 and r2 are random decimals between zero and one.


#### 4.2.2. Procedure of PSO

According to the initialization and velocity updating strategy described above, the basic procedure of the PSO can be depicted as follows and the detailed process is shown in Algorithm2:Step 1: initialize particle swarm related parameters, including swarm size, position information contained in each particle, and particle velocity.Step 2: normalize the information contained in the particles to obtain the vehicle assignment and sequence of customers which can further determine the visit time of each vehicle.Step 3: judge whether the solution generated by each particle meets the capacity constraint. If not, the fitness value of the particle is the maximum value, and if it meets the constraint, the fitness value of the particle can be calculated.Step 4: compare the fitness value with the personal optimal value of the particle. If it is better than the historical personal optimal value, update the personal optimal value and record the position information of the particle.Step 5: update the particle position information and particle velocity according to the formula.Step 6: determine whether the exit condition is satisfied (the gap is small enough or the maximum number of iterations is reached). If the condition is met, exit the loop; otherwise return to step 2.Algorithm 2: Particle swarm optimization (PSO).
(1)
Initialize swarm(2)
For each particle(3)
Parameter conversion (generate customer sequence, vehicle assignment)(4)
Check the constraints(5)
Fitness evaluation(6)
Update best personal solution and global solution(7)
End for(8)
While termination condition not met(9)
Update particle position and velocity(10)
Goto line 2(11)
End while


#### 4.2.3. Coding Rule of PSO

According to line 3 of the particle swarm flow above, parameter transformation refers to the transformation of the position information contained in the particles into the variables needed to solve the model, that is, the vehicle assignment and customer sequence. Figure2 shows the process of assigning vehicles to customers by taking 6 customers of hazardous materials type 1 as example. Figure 2(a) shows the subset of customers for hazardous materials type 1 that need transportation services and the corresponding particle positions. Figure 2(b) shows a subset of vehicles capable of transporting type 1 hazardous materials. Figure 2(c) shows the vehicle that each customer is assigned to based on the location information contained in its particle. If all customers are assigned to one single vehicle and customer demand exceeds the capacity of the vehicle, we then need to adjust the vehicle allocation strategy. We add the requirements in ascending order according to the particle locations of the customers and assign the customers exceeding the capacity to the vehicles with smaller serial numbers in turn until the requirements are fully met.Figure 2
Vehicle assignment mechanism. (a) PSO particle. (b) Vehicle subset. (c) Vehicle ID assigned to each customer.Figure3 depicts the generation of vehicle visiting sequence, and the number of customers is 5. Figure 3(a) represents the location information of each particle, and the customer sequence obtained after sorting it in ascending order is depicted in Figure 3(b). Each iteration of PSO brings about a change in particle position, and a new sequence will be generated after sorting. Better paths can be found by constantly adjusting the sequence.Figure 3
Generation of customer sequence. (a) PSO particle. (b) Customer sequence of each vehicle.


## 4.1. Bidirectional Tuning Heuristic

The main idea of bidirectional tuning heuristic is to transform the model into several interrelated subproblems. Solving the subproblems using some fixed decision variables, we can obtain the remaining decision variables and then use the them to gain other decision variables. In the iterative process, if the value of the objective is better than the optimal value, the optimal solution needs to be updated. The exit mechanism of the algorithm can be determined by the number of iterations or by the rules to judge whether the objective value can get a better solution.Before we start solving, the decision variables need to be classified. The decision variables in this paper can be divided into three different types according to their definitions. The first type isαi,k, the variables used to decide whether a vehicle visits a customer. The second type is μk,v and ηk, used to determine the type of hazardous material each vehicle transports. The other decision variables contain the dimensions of random scenarios and are limited by the first two kinds of decision variables, so they are always taken as the variables to be solved. We determine the customer that each vehicle visits and then optimize the transportation types with fixed αi,k. Then, we move to determining the transportation types and optimize the customers with fixed μk,v and ηk. The solving is repeated until the iteration reaches upper limit or no improvement can be obtained. The basic procedure can be defined as follows and the detailed process is shown in Algorithm 1.Step 1: solve the model to find a feasible solution, which can be taken as the initial solution. Record the initial decision variables and the objective value.Step 2: based onαi,k from initialization, determine decision variable ηk and μk,v. Compare the obtained objective value with the historical optimal value, and if the value is better than the historical optimal value, update the optimal value.Step 3: solve the model with fixedηk and μk,v to obtain αi,k. Update the objective value.Step 4: judge the exit condition; if it is satisfied, the optimal solution is obtained and the algorithm stops; otherwise, go to step 2.Algorithm 1: Bidirectional tuning heuristic.
(1)
Initialization(2)
Solve the model, and determineαi,k(3)
Update the optimal solution(4)
While (n < MaxIter)(5)
Solve model with fixedαi,k, determine μk,v and ηk(6)
Update the optimal solution(7)
Solve model with fixedμk,v and ηk, determine αi,k(8)
Update the optimal solution(9)
End while


## 4.2. Particle Swarm Optimization (PSO)

Particle swarm optimization (PSO) has been studied extensively since it was proposed [17]. As one of the swarm intelligence algorithms, PSO has been widely used in vehicle routing problem in the past few years, for example, vehicle routing problem with multiple pickup and delivery [18] and multidepot multitrip vehicle routing problem [19].


### 4.2.1. Initialization and Velocity Updating Strategy

The particle swarm optimization algorithm not only considers the personal optimal value, but also records the global optimal value of the entire group in the process of searching for solutions. It has excellent performance in solving optimization problems. The movement of particles depends on the update of their position and velocity; each particle contains the current position and the personal optimal position. The current position of the particles is adjusted by updating the velocity formula, the solution is judged by the specific fitness value, and the optimal position of the entire group is updated in the global range.In the model we proposed, the decision variablesαi,k and μk,v are obviously related to other decision variables. When αi,k is fixed, ηk can be determined by constraints (4) and γi,j,k,r,s can also be determined by constraints (5)–(7). Furthermore, the decision variables βi,k,s, related to time, are limited by γi,j,k,r,s according to constraints (13). Then the variable δi,k,s can be derived from constraints (13). At this point, we can judge whether the service time δi,k,s of each customer is reasonable. If it is not reasonable, we update the time according to the time sequence of starting time and the processing time. It is also necessary to trace back to update variable βi,k,s, which is the subsequent vehicle arrival time of the route. When all decision variables related to time are determined, we can easily confirm the decision variable εi,k,k′,s. When we know both αi,k and μk,v, the load variable ζi,k,s,v can also be determined. We now have clarified the relationship between all the decision variables.Based on the description above, we setμk,v as the outer variable of the algorithm. And the particles that we set up contain information to find the vehicles assigned to customer and to generate visit sequence of each vehicle, which is similar to the approach taken in this work [20]. We define the particle as Pmn, and its velocity can be defined as Velmn. The velocity update formula can be expressed by formula (25) and the position update formula is (26). We use PBestmn to represent the best personal optimal position of the particle and GBestmn to represent the global optimal position, where m means the number of particles (we take 300) and n means the current iteration number:(25)Velmn+1=wnVelmn+c1r1PBestmn−Pmn+c2r2GBestmn−Pmn,(26)Pmn+1=Pmn+Velmn.In the formula for particle velocity update,wn denotes inertia weight and it is calculated by equation wn=wmax−wmax−wminn/N, in which N is the number of total iterations (we take 50) and wmax and wmin are maximum and minimum inertia weight (we take 0.9 and 0.4). This makes the particle swarm have strong global convergence ability at the beginning, but strengthens local convergence ability with the increase of iteration. c1 and c2 are acceleration weights (we take 0.683 for both), while r1 and r2 are random decimals between zero and one.


### 4.2.2. Procedure of PSO

According to the initialization and velocity updating strategy described above, the basic procedure of the PSO can be depicted as follows and the detailed process is shown in Algorithm2:Step 1: initialize particle swarm related parameters, including swarm size, position information contained in each particle, and particle velocity.Step 2: normalize the information contained in the particles to obtain the vehicle assignment and sequence of customers which can further determine the visit time of each vehicle.Step 3: judge whether the solution generated by each particle meets the capacity constraint. If not, the fitness value of the particle is the maximum value, and if it meets the constraint, the fitness value of the particle can be calculated.Step 4: compare the fitness value with the personal optimal value of the particle. If it is better than the historical personal optimal value, update the personal optimal value and record the position information of the particle.Step 5: update the particle position information and particle velocity according to the formula.Step 6: determine whether the exit condition is satisfied (the gap is small enough or the maximum number of iterations is reached). If the condition is met, exit the loop; otherwise return to step 2.Algorithm 2: Particle swarm optimization (PSO).
(1)
Initialize swarm(2)
For each particle(3)
Parameter conversion (generate customer sequence, vehicle assignment)(4)
Check the constraints(5)
Fitness evaluation(6)
Update best personal solution and global solution(7)
End for(8)
While termination condition not met(9)
Update particle position and velocity(10)
Goto line 2(11)
End while


### 4.2.3. Coding Rule of PSO

According to line 3 of the particle swarm flow above, parameter transformation refers to the transformation of the position information contained in the particles into the variables needed to solve the model, that is, the vehicle assignment and customer sequence. Figure2 shows the process of assigning vehicles to customers by taking 6 customers of hazardous materials type 1 as example. Figure 2(a) shows the subset of customers for hazardous materials type 1 that need transportation services and the corresponding particle positions. Figure 2(b) shows a subset of vehicles capable of transporting type 1 hazardous materials. Figure 2(c) shows the vehicle that each customer is assigned to based on the location information contained in its particle. If all customers are assigned to one single vehicle and customer demand exceeds the capacity of the vehicle, we then need to adjust the vehicle allocation strategy. We add the requirements in ascending order according to the particle locations of the customers and assign the customers exceeding the capacity to the vehicles with smaller serial numbers in turn until the requirements are fully met.Figure 2
Vehicle assignment mechanism. (a) PSO particle. (b) Vehicle subset. (c) Vehicle ID assigned to each customer.Figure3 depicts the generation of vehicle visiting sequence, and the number of customers is 5. Figure 3(a) represents the location information of each particle, and the customer sequence obtained after sorting it in ascending order is depicted in Figure 3(b). Each iteration of PSO brings about a change in particle position, and a new sequence will be generated after sorting. Better paths can be found by constantly adjusting the sequence.Figure 3
Generation of customer sequence. (a) PSO particle. (b) Customer sequence of each vehicle.


## 4.2.1. Initialization and Velocity Updating Strategy

The particle swarm optimization algorithm not only considers the personal optimal value, but also records the global optimal value of the entire group in the process of searching for solutions. It has excellent performance in solving optimization problems. The movement of particles depends on the update of their position and velocity; each particle contains the current position and the personal optimal position. The current position of the particles is adjusted by updating the velocity formula, the solution is judged by the specific fitness value, and the optimal position of the entire group is updated in the global range.In the model we proposed, the decision variablesαi,k and μk,v are obviously related to other decision variables. When αi,k is fixed, ηk can be determined by constraints (4) and γi,j,k,r,s can also be determined by constraints (5)–(7). Furthermore, the decision variables βi,k,s, related to time, are limited by γi,j,k,r,s according to constraints (13). Then the variable δi,k,s can be derived from constraints (13). At this point, we can judge whether the service time δi,k,s of each customer is reasonable. If it is not reasonable, we update the time according to the time sequence of starting time and the processing time. It is also necessary to trace back to update variable βi,k,s, which is the subsequent vehicle arrival time of the route. When all decision variables related to time are determined, we can easily confirm the decision variable εi,k,k′,s. When we know both αi,k and μk,v, the load variable ζi,k,s,v can also be determined. We now have clarified the relationship between all the decision variables.Based on the description above, we setμk,v as the outer variable of the algorithm. And the particles that we set up contain information to find the vehicles assigned to customer and to generate visit sequence of each vehicle, which is similar to the approach taken in this work [20]. We define the particle as Pmn, and its velocity can be defined as Velmn. The velocity update formula can be expressed by formula (25) and the position update formula is (26). We use PBestmn to represent the best personal optimal position of the particle and GBestmn to represent the global optimal position, where m means the number of particles (we take 300) and n means the current iteration number:(25)Velmn+1=wnVelmn+c1r1PBestmn−Pmn+c2r2GBestmn−Pmn,(26)Pmn+1=Pmn+Velmn.In the formula for particle velocity update,wn denotes inertia weight and it is calculated by equation wn=wmax−wmax−wminn/N, in which N is the number of total iterations (we take 50) and wmax and wmin are maximum and minimum inertia weight (we take 0.9 and 0.4). This makes the particle swarm have strong global convergence ability at the beginning, but strengthens local convergence ability with the increase of iteration. c1 and c2 are acceleration weights (we take 0.683 for both), while r1 and r2 are random decimals between zero and one.


## 4.2.2. Procedure of PSO

According to the initialization and velocity updating strategy described above, the basic procedure of the PSO can be depicted as follows and the detailed process is shown in Algorithm2:Step 1: initialize particle swarm related parameters, including swarm size, position information contained in each particle, and particle velocity.Step 2: normalize the information contained in the particles to obtain the vehicle assignment and sequence of customers which can further determine the visit time of each vehicle.Step 3: judge whether the solution generated by each particle meets the capacity constraint. If not, the fitness value of the particle is the maximum value, and if it meets the constraint, the fitness value of the particle can be calculated.Step 4: compare the fitness value with the personal optimal value of the particle. If it is better than the historical personal optimal value, update the personal optimal value and record the position information of the particle.Step 5: update the particle position information and particle velocity according to the formula.Step 6: determine whether the exit condition is satisfied (the gap is small enough or the maximum number of iterations is reached). If the condition is met, exit the loop; otherwise return to step 2.Algorithm 2: Particle swarm optimization (PSO).
(1)
Initialize swarm(2)
For each particle(3)
Parameter conversion (generate customer sequence, vehicle assignment)(4)
Check the constraints(5)
Fitness evaluation(6)
Update best personal solution and global solution(7)
End for(8)
While termination condition not met(9)
Update particle position and velocity(10)
Goto line 2(11)
End while


## 4.2.3. Coding Rule of PSO

According to line 3 of the particle swarm flow above, parameter transformation refers to the transformation of the position information contained in the particles into the variables needed to solve the model, that is, the vehicle assignment and customer sequence. Figure2 shows the process of assigning vehicles to customers by taking 6 customers of hazardous materials type 1 as example. Figure 2(a) shows the subset of customers for hazardous materials type 1 that need transportation services and the corresponding particle positions. Figure 2(b) shows a subset of vehicles capable of transporting type 1 hazardous materials. Figure 2(c) shows the vehicle that each customer is assigned to based on the location information contained in its particle. If all customers are assigned to one single vehicle and customer demand exceeds the capacity of the vehicle, we then need to adjust the vehicle allocation strategy. We add the requirements in ascending order according to the particle locations of the customers and assign the customers exceeding the capacity to the vehicles with smaller serial numbers in turn until the requirements are fully met.Figure 2
Vehicle assignment mechanism. (a) PSO particle. (b) Vehicle subset. (c) Vehicle ID assigned to each customer.Figure3 depicts the generation of vehicle visiting sequence, and the number of customers is 5. Figure 3(a) represents the location information of each particle, and the customer sequence obtained after sorting it in ascending order is depicted in Figure 3(b). Each iteration of PSO brings about a change in particle position, and a new sequence will be generated after sorting. Better paths can be found by constantly adjusting the sequence.Figure 3
Generation of customer sequence. (a) PSO particle. (b) Customer sequence of each vehicle.


## 5. Numerical Experiments

In this part, we use extensive numerical experiments to verify whether the heuristic proposed in this article is efficient as well as the effectiveness of the model. All the experiments were conducted on a workstation with two Xeon E5-2643 CPUs (24 cores) of 3.4 GHz and 128 G of memory. All programs were compiled with C#, and the mixed integer programming model was solved by CPLEX.


### 5.1. Parameter Setting

We have designed multiple groups of randomly generated numerical examples for both small-scale and large-scale problems. In the small-scale example, we calculated the situation of 5, 10, and 15 customers with 2 or 3 types of hazardous materials, respectively, while in larger scale, we calculated the problem with customer number from 20 to 70 of three different hazardous materials. Since a customer could have multiple transportation demands with heterogeneous vehicles, each customer can be viewed as multiple customers, which makes solving on a larger scale more difficult. The remaining parameters were generated under some specific rules; for example, each customer had at least 1 type of hazardous materials transportation demand and the quantity was limited to a certain range. Each vehicle had a fixed latest working time. The three cost coefficients are estimated under China’s hazardous material transportation market.ctran=200 means the transportation cost in an hour which is calculated by speed, fuel consumption, and fuel cost (60 km/h × 0.66 L/km × 5.05 yuan/L), crisk=10000 is a large number based on practical accident losses, and ctime=50 represents additional work cost per hour.


### 5.2. Performance of Two Heuristics

The comparison between two heuristics and CPLEX is illustrated in Table2, where Z means the objective value, T means computation time, and GAP is calculated by equation Z−ZC/ZC. The subscripts C, B, P represent CPLEX, bidirectional tuning heuristic, and PSO, respectively.Table 2
Comparison of two heuristics and CPLEX in small-scale problems.
Cases IDCPLEXBituningComparisonPSOComparisonZCTCZBTBTB/TCGAPB (%)ZPTPTP/TCGAPP (%)A5-3-2-1102370.2102370.21.00.00102659.145.50.27A5-3-2-2101980.2101980.21.00.00102334.321.50.34A5-3-2-394210.294210.10.50.0094614.221.00.42A5-3-3-1147392.0147391.30.70.00149217.33.71.23A5-3-3-2142491.4142490.20.10.001428724.817.70.27A5-3-3-3137891.7137890.60.40.001394657.333.71.14A10-3-2-112041144.1120416.30.00.001224824.90.21.72A10-3-2-212217335.0122171.10.00.001245731.30.11.96A10-3-2-312527606.3125274.30.00.001263629.70.00.87A10-3-3-115570358.141557022.20.10.001570342.70.10.85A10-3-3-216655602.91669263.80.10.221688932.70.11.40A10-3-3-319774610.31977498.90.20.002000339.10.11.16A15-6-2-1235511576.223551215.20.10.002407961.10.02.24A15-6-2-2211151352.121115144.50.10.002157760.20.02.19A15-6-2-3247902161.124790218.70.10.002516970.40.01.53Avg.—516.8—51.80.30.01—33.39.61.17We tested fifteen cases in small scale, shown in Table2. The case ID A5-3-2-1 means 5 customers, 3 vehicles, and 2 hazardous material types of instance 1. When the customer turns to 20, the proposed model cannot obtain a feasible solution in 7200 seconds. It is obvious that CPLEX can only maintain high calculation efficiency in extremely small-scale examples which contain less than 10 customers. As the customer number and types of hazardous materials increase, the calculating time spent by CPLEX grows exponentially. Meanwhile, the bidirectional tuning heuristic can also obtain the optimal solution without increasing time significantly. In terms of small-scale problem, bidirectional tuning performed the best among the three methods with an average of 51.8 seconds of computation time and average gap of about 0.01%. The PSO spent the shortest time with an average of 33.3 seconds and average gap of about 1.17%, which could be considered as a near-optimal solution.We list a group of large-scale instances of customer number from 20 to 70 to further verify the effectiveness of the algorithm. The result is shown in Table3, where each case ID means the same as in Table 2. The GAP is calculated by equation ZC−Z/Z. Due to the complexity of the model, when the customer number is over 30 with vehicle number of 10, the bidirectional tuning heuristic cannot find optimal solution in 7200 seconds. We find that even if it takes more time to search, the solution is not improved. Therefore, we had to use the results computed by solving the problem after relaxation to judge the effectiveness of algorithm. The method proved to be feasible in other research [21]. However, even after relaxation, CPLEX still takes a long time to find a feasible solution. A long calculation time means no practical value, so it is acceptable that we set the computation time as 7200 seconds to get the solution within the effective time. From the data of Table 3, we can see an evident improvement of the PSO’s gap value. Moreover, the improvement increases with the scale of problem. The computation time increases linearly as the problem size increases, which is much better than CPLEX and bidirectional tuning heuristic.Table 3
Comparison of two heuristics and CPLEX in large-scale problems.
Case IDCPLEX (7200)BituningPSOComparisonZCZBTBZPTPTP/TBGAPB (%)GAPP (%)A20-6-3-13650735915809.136404124.40.151.650.28A20-6-3-23541234148815.334972119.70.153.701.26A30-10-3-149132—>720046878341.3——4.81A30-10-3-253486—>720050540438.5——5.83A40-15-3-169194—>720063715923.4——8.60A40-15-3-271866—>720066598732.2——7.91A50-20-3-186210—>7200782381229.2——10.19A50-20-3-2100626—>7200880681521.2——14.26A60-25-3-1115044—>7200990652671.9——16.13A60-25-3-2114854—>7200979653016.5——17.24A70-30-3-1172731—>72001347307407.3——28.21A70-30-3-2139333—>72001134366175.8——22.83To summarize, in small-scale problem, bidirectional tuning heuristic can obtain the same optimal solution as CPLEX, while PSO can obtain near-accurate solution. On the contrary, PSO spent less time to find the solution, followed by the bidirectional tuning heuristic. In large-scale problem, PSO shows great capability with high-quality solution in a short time.


## 5.1. Parameter Setting

We have designed multiple groups of randomly generated numerical examples for both small-scale and large-scale problems. In the small-scale example, we calculated the situation of 5, 10, and 15 customers with 2 or 3 types of hazardous materials, respectively, while in larger scale, we calculated the problem with customer number from 20 to 70 of three different hazardous materials. Since a customer could have multiple transportation demands with heterogeneous vehicles, each customer can be viewed as multiple customers, which makes solving on a larger scale more difficult. The remaining parameters were generated under some specific rules; for example, each customer had at least 1 type of hazardous materials transportation demand and the quantity was limited to a certain range. Each vehicle had a fixed latest working time. The three cost coefficients are estimated under China’s hazardous material transportation market.ctran=200 means the transportation cost in an hour which is calculated by speed, fuel consumption, and fuel cost (60 km/h × 0.66 L/km × 5.05 yuan/L), crisk=10000 is a large number based on practical accident losses, and ctime=50 represents additional work cost per hour.


## 5.2. Performance of Two Heuristics

The comparison between two heuristics and CPLEX is illustrated in Table2, where Z means the objective value, T means computation time, and GAP is calculated by equation Z−ZC/ZC. The subscripts C, B, P represent CPLEX, bidirectional tuning heuristic, and PSO, respectively.Table 2
Comparison of two heuristics and CPLEX in small-scale problems.
Cases IDCPLEXBituningComparisonPSOComparisonZCTCZBTBTB/TCGAPB (%)ZPTPTP/TCGAPP (%)A5-3-2-1102370.2102370.21.00.00102659.145.50.27A5-3-2-2101980.2101980.21.00.00102334.321.50.34A5-3-2-394210.294210.10.50.0094614.221.00.42A5-3-3-1147392.0147391.30.70.00149217.33.71.23A5-3-3-2142491.4142490.20.10.001428724.817.70.27A5-3-3-3137891.7137890.60.40.001394657.333.71.14A10-3-2-112041144.1120416.30.00.001224824.90.21.72A10-3-2-212217335.0122171.10.00.001245731.30.11.96A10-3-2-312527606.3125274.30.00.001263629.70.00.87A10-3-3-115570358.141557022.20.10.001570342.70.10.85A10-3-3-216655602.91669263.80.10.221688932.70.11.40A10-3-3-319774610.31977498.90.20.002000339.10.11.16A15-6-2-1235511576.223551215.20.10.002407961.10.02.24A15-6-2-2211151352.121115144.50.10.002157760.20.02.19A15-6-2-3247902161.124790218.70.10.002516970.40.01.53Avg.—516.8—51.80.30.01—33.39.61.17We tested fifteen cases in small scale, shown in Table2. The case ID A5-3-2-1 means 5 customers, 3 vehicles, and 2 hazardous material types of instance 1. When the customer turns to 20, the proposed model cannot obtain a feasible solution in 7200 seconds. It is obvious that CPLEX can only maintain high calculation efficiency in extremely small-scale examples which contain less than 10 customers. As the customer number and types of hazardous materials increase, the calculating time spent by CPLEX grows exponentially. Meanwhile, the bidirectional tuning heuristic can also obtain the optimal solution without increasing time significantly. In terms of small-scale problem, bidirectional tuning performed the best among the three methods with an average of 51.8 seconds of computation time and average gap of about 0.01%. The PSO spent the shortest time with an average of 33.3 seconds and average gap of about 1.17%, which could be considered as a near-optimal solution.We list a group of large-scale instances of customer number from 20 to 70 to further verify the effectiveness of the algorithm. The result is shown in Table3, where each case ID means the same as in Table 2. The GAP is calculated by equation ZC−Z/Z. Due to the complexity of the model, when the customer number is over 30 with vehicle number of 10, the bidirectional tuning heuristic cannot find optimal solution in 7200 seconds. We find that even if it takes more time to search, the solution is not improved. Therefore, we had to use the results computed by solving the problem after relaxation to judge the effectiveness of algorithm. The method proved to be feasible in other research [21]. However, even after relaxation, CPLEX still takes a long time to find a feasible solution. A long calculation time means no practical value, so it is acceptable that we set the computation time as 7200 seconds to get the solution within the effective time. From the data of Table 3, we can see an evident improvement of the PSO’s gap value. Moreover, the improvement increases with the scale of problem. The computation time increases linearly as the problem size increases, which is much better than CPLEX and bidirectional tuning heuristic.Table 3
Comparison of two heuristics and CPLEX in large-scale problems.
Case IDCPLEX (7200)BituningPSOComparisonZCZBTBZPTPTP/TBGAPB (%)GAPP (%)A20-6-3-13650735915809.136404124.40.151.650.28A20-6-3-23541234148815.334972119.70.153.701.26A30-10-3-149132—>720046878341.3——4.81A30-10-3-253486—>720050540438.5——5.83A40-15-3-169194—>720063715923.4——8.60A40-15-3-271866—>720066598732.2——7.91A50-20-3-186210—>7200782381229.2——10.19A50-20-3-2100626—>7200880681521.2——14.26A60-25-3-1115044—>7200990652671.9——16.13A60-25-3-2114854—>7200979653016.5——17.24A70-30-3-1172731—>72001347307407.3——28.21A70-30-3-2139333—>72001134366175.8——22.83To summarize, in small-scale problem, bidirectional tuning heuristic can obtain the same optimal solution as CPLEX, while PSO can obtain near-accurate solution. On the contrary, PSO spent less time to find the solution, followed by the bidirectional tuning heuristic. In large-scale problem, PSO shows great capability with high-quality solution in a short time.


## 6. Conclusion

The heterogeneous vehicle routing problem for hazardous materials is a practical problem that deserves to be studied. It is difficult to balance the total cost under the premise of considering many influencing factors. In this paper, we study the multidepot VRPHM with risk analysis and propose mixed integer programming to transform the actual problem into a model, which can be solved by mathematical methods. The main contributions we made in this research are listed as follows: First, we consider comprehensive factors including transportation costs, multiple depots, heterogeneous vehicles, risks, and multiple accident scenarios which make the study closer to the actual situation. Furthermore, the risk is defined as uncertainty in different scenarios, and we consider heterogeneous vehicles which are rarely studied. For solving the problem in an efficient way, we design a bidirectional tuning heuristic and particle swarm optimization (PSO) to be applied to different scales of problem.The numerical experiments show that our proposed algorithm can be used in small-scale problem with faster speed than CPLEX. And in large-scale problem (70 customers, 30 vehicles, and 3 types), PSO can find high-quality feasible solution in acceptable time.However, this model can still be improved in terms of, for example, the definition of risk with uncertainty. The data we tested in numerical experiments are randomly generated in a reasonable range, which may be a little difference from the actual situation. And more algorithms can be tried to solve this complex problem.



---

*Source: 8839628-2020-08-28.xml*