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.
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. 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.
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 network , where the set represents the union of the set composed of depots and the set composed of customer points, and is an edge set. The network is depicted clearly in Figure 1.
Let the set represent the vehicles of fixed cost . The maximum load of each vehicle is . The vehicle determines whether it can visit the customer based on the type of hazardous materials, described using parameter . Each vehicle must depart from its own depot 0 to the customer . After process time , the vehicle can travel from customer to customer . We introduce a virtual node 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 , where the route set is . Besides, the travel time is different from 0 to each customer if the vehicle belongs to different depots. Each vehicle has the latest working time , and an additional overtime costs 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 . The risk under different scenarios is related to the route chosen by the vehicle. We use to describe the probability of accident that may occur from the customer to the customer .
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).
3.3. Mathematical Model
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 can serve customer . Constraints (4) show that if vehicle serves customer , 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 from customer to customer under scenario 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 and . 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:
: binary, equals 1 if and .
Newly defined constraints are as follows:
4. Algorithm Strategies
The 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.
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 , the variables used to decide whether a vehicle visits a customer. The second type is and , 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 . Then, we move to determining the transportation types and optimize the customers with fixed and . 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 from initialization, determine decision variable and . 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 and to obtain . 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.
|
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 and are obviously related to other decision variables. When is fixed, can be determined by constraints (4) and can also be determined by constraints (5)–(7). Furthermore, the decision variables , related to time, are limited by according to constraints (13). Then the variable can be derived from constraints (13). At this point, we can judge whether the service time 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 , 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 . When we know both and , the load variable can also be determined. We now have clarified the relationship between all the decision variables.
Based on the description above, we set 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 , and its velocity can be defined as . The velocity update formula can be expressed by formula (25) and the position update formula is (26). We use to represent the best personal optimal position of the particle and to represent the global optimal position, where means the number of particles (we take 300) and means the current iteration number:
In the formula for particle velocity update, denotes inertia weight and it is calculated by equation , in which is the number of total iterations (we take 50) and and 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. and are acceleration weights (we take 0.683 for both), while and 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 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.
|
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. 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 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.
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. 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), is a large number based on practical accident losses, and represents additional work cost per hour.
5.2. Performance of Two Heuristics
The comparison between two heuristics and CPLEX is illustrated in Table 2, where means the objective value, means computation time, and is calculated by equation . The subscripts , , represent CPLEX, bidirectional tuning heuristic, and PSO, respectively.
We tested fifteen cases in small scale, shown in Table 2. 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 Table 3, where each case ID means the same as in Table 2. The is calculated by equation . 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.
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.
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.
Data Availability
The 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
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study was supported by Beijing Nova Program of Science and Technology (no. Z191100001119029) and Research on Multi-Depot Heterogeneous Vehicle Routing Problem.