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

(is 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. (e 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. (e 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.


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.
us, 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]. e 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. is 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.
e 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. e 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. e 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. e 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. erefore, we design a bidirectional tuning heuristic and a particle swarm optimization algorithm (PSO) to solve problems at different scales. e results show that the algorithms we designed can effectively solve larger-scale problems and can also find near-accurate solutions to small-scale problems. e reminder of this paper is organized as follows. Section 2 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.

Related Literature
e 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. e 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. e 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. e 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. e 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][13][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. ird, two effective algorithms are developed.

Problem Definition and Formulation
In this section, we describe the multidepot VRPHM, define the parameters, and then formulate a MIP model. e 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. e 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. e fixed costs of starting a vehicle and the variable costs per unit of travel distance are also different. e transportation network of hazardous materials can be defined on a directed network G � (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, j) | i, j ∈ O is an edge set. e network is depicted clearly in Figure 1.
Let the set K represent the vehicles of fixed cost g k . e maximum load of each vehicle is q k . e vehicle determines whether it can visit the customer based on the type of hazardous materials, described using parameter c i,k . Each vehicle must depart from its own depot 0 to the customer i ∈ N. After process time t i , 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. is virtual node is the same node as 0.
e transportation time is different according to the chosen route which can be expressed as d i,j,r,k , where the route set is R i,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 y k , and an additional overtime costs c time 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. is paper introduces the concept of scenarios with different probability to eliminate such uncertainties. e set of scenarios is represented by S. e risk under different scenarios is related to the route chosen by the vehicle. We use w i,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

Notations.
In this section, we define all the parameters and decision variables (Table1).

Mathematical Model.
Minimize k∈K g k η k + 1 |S| s∈S k∈K j∈N∪ e(N) subject to i∈N j∈N r∈R i,j  Scientific Programming k∈K k′∈K e 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.

Scientific Programming 5
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)

Algorithm Strategies
e model presented in Section 3 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.

Bidirectional Tuning Heuristic.
e 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. e 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. e decision variables in this paper can be divided into three different types according to their definitions. e first type is α i,k , the variables used to decide whether a vehicle visits a customer. e second type is μ k,v and η k , used to determine the type of hazardous material each vehicle transports. e 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 . en, we move to determining the transportation types and optimize the customers with fixed μ k,v and η k . e solving is repeated until the iteration reaches upper limit or no improvement can be obtained. e 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.

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].

Initialization and Velocity Updating Strategy.
e 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. e movement of particles depends on the update of their position and velocity; each particle contains the current position and the personal optimal position. e 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 c 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 c i,j,k,r,s according to constraints (13). en 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 P n m , and its velocity can be defined as Vel n m . e velocity update formula can be expressed by formula (25) and the position update formula is (26). We use PBest n m to represent the best personal optimal position of the particle and GBest n m to represent the global optimal position, where m means the number of particles (we take 300) and n means the current iteration number: Vel n+1 m � w n Vel n m + c 1 r 1 PBest n m − P n m + c 2 r 2 GBest n m − P n m , In the formula for particle velocity update, w n denotes inertia weight and it is calculated by equation w n � w max − (w max − w min )n/N, in which N is the number of total iterations (we take 50) and w max and w min are maximum and minimum inertia weight (we take 0.9 and 0.4).
is makes the particle swarm have strong global convergence ability at the beginning, but strengthens local convergence ability with the increase of iteration. c 1 and c 2 are acceleration weights (we take 0.683 for both), while r 1 and r 2 are random decimals between zero and one.

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 Algorithm 2: 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.

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. Figure 2 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 3 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 (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 ALGORITHM 1: Bidirectional tuning heuristic. Scientific Programming 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.

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.

Parameter Setting.
We have designed multiple groups of randomly generated numerical examples for both smallscale 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. e 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. e three  (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 ALGORITHM 2: Particle swarm optimization (PSO). 8 Scientific Programming cost coefficients are estimated under China's hazardous material transportation market. c tran � 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), c risk � 10000 is a large number based on practical accident losses, and c time � 50 represents additional work cost per hour.

Performance of Two
Heuristics. e comparison between two heuristics and CPLEX is illustrated in Table 2, where Z means the objective value, T means computation time, and GAP is calculated by equation (Z − Z C )/Z C . e subscripts C, B, P represent CPLEX, bidirectional tuning heuristic, and PSO, respectively.
We tested fifteen cases in small scale, shown in Table 2. e 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%. e 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 nearoptimal solution.
We list a group of large-scale instances of customer number from 20 to 70 to further verify the effectiveness of the algorithm. e result is shown in Table 3, where each case ID means the same as in Table 2. e GAP is calculated by equation (Z C − 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  erefore, we had to use the results computed by solving the problem after relaxation to judge the effectiveness of algorithm. e 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. e computation time increases linearly as the problem size increases, which is much better than CPLEX and bidirectional tuning heuristic.
To 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.

Conclusion
e 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. e 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.
e 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. e 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.
Data Availability e raw/processed data required to reproduce the findings cannot be shared at this time as the data also form a part of an ongoing study.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.