Half Open Multi-Depot Heterogeneous Vehicle Routing Problem for Hazardous Materials Transportation

: How to reduce the accidents of hazardous materials has become an important and urgent research topic in the safety management of hazardous materials. In this study, we focus on the half open multi-depot heterogeneous vehicle routing problem for hazardous materials transportation. The goal is to determine the vehicle allocation and the optimal route with minimum risk and cost for hazardous materials transportation. A novel transportation risk model is presented considering the variation of vehicle loading, vehicle types, and hazardous materials category. In order to balance the transportation risk and the transportation cost, we propose a bi-objective mixed integer programming model. A hybrid intelligent algorithm is developed based on the ε -constraint method and genetic algorithm to obtain the Pareto optimal solutions. Numerical experiments are provided to demonstrate the effectiveness of the proposed model. Compared with the close multi-depot heterogeneous vehicle routing problem, the average risk and cost obtained by the proposed bi-objective mixed integer programming model can be reduced by 3.99% and 2.01%, respectively. In addition, compared with the half open multi-depot homogeneous vehicle routing problem, the cost is signiﬁcantly reduced with the acceptable risk.


Introduction
The development of industry and manufacturing increases the transportation of various hazardous materials (hazmat). Hazmat can be transported by road, railway, sea, and pipeline, in which the road transportation is less affected by infrastructures. Thus, hazmat road transportation is still the main transportation mode [1]. In China, it is estimated that approximately 80% of hazmat are transported by road annually and the amount is about 1 billion tons. However, hazmat transportation generates high potential risk due to the low safety of roads and multiple risk factors. Hazmat transportation accidents may lead to leakage, explosion, poisoning, and other accidents, causing casualties, damage to property, and environmental damage. From 2013 to 2017, there were only 356 hazmat transportation accidents in China, but the accidents caused 855 deaths, 2980 injuries, and huge economic loss [2]. For example, on 16 January 2015, a collision accident of a gasoline tanker with three trucks on the Yinmachi Bridge of Rongwu Expressway caused 12 deaths and 11 million RMB economy loss; on 23 May 2017, a gas tanker exploded and caused 15 deaths and 3 injuries on the Zhangjiakou Shijiazhuang Expressway. Therefore, how to effectively reduce the accidents of hazmat transportation has become an important and urgent research topic in the safety management of hazmat.
Hazmat vehicle routing programming and risk evaluation approaches are mainly researched to reduce the hazmat transportation risk [3]. However, there are only a few Sustainability 2021, 13, 1262 2 of 17 works focused on the multi-depot vehicle routing problem (MDVRP) for hazmat transportation, especially the variants of the MDVRP based on real-world applications. The vehicle routing problem (VRP) is a classic combinational optimization problem for logistics transportation that aims to obtain the optimal route with the minimum cost of the constraints on capacities and customer demand. Due to the complexity of real-world transportation problems, there are many variants of the classical VRP in previous literature such as the heterogeneous VRP [4], the dynamic VRP [5], the MDVRP [6], the open VRP [7], etc. Similarly, there is some literature that has studied VRP and its variants for hazmat transportation in recent years. For instance, Du et al. [8] proposed a MDVRP for hazmat transportation to minimize total transportation risk, in which each customer can be serviced by multiple depots of the constraints on vehicles capacities and customer demand. Furthermore, Du et al. [3] formulated a multi-level programming model that considering multiple factors in real-world applications for hazmat MDVRP. Ma et al. [9] focused on the MDVRP for hazmat transportation considering energy consumption and formulated a multi-objective optimization model to minimize total transportation energy consumption and transportation risk. In most VRPs, vehicles are assumed to be homogeneous but this is unrealistic in real-world transportation problems. Moreover, some researchers considered heterogeneous VRP for hazmat transportation, in which vehicles have different characteristics (capacities, fixed cost, etc.). Bula et al. [10] focused on the heterogeneous VRP for hazmat transportation to determine the set of routes with minimizing total transportation risk. For economic reasons, Bula et al. [11] proposed a bi-objective programming model for hazmat heterogeneous VRP that simultaneously minimized total routing risk and total transportation cost. Wang et al. [1] studied the multi-depot heterogeneous VRP for hazmat with considering multiple depots, heterogeneous vehicles, and multiple accident scenarios. To the best of our knowledge, MDVRP with the half open route and heterogeneous vehicles for hazmat transportation is rarely studied. Thus, we simultaneously consider multiple depots, heterogeneous vehicles, half open route with minimizing total transportation risk and total transportation cost. The half open route is different from the close route and the open route. The close route is referred to as a route where the vehicle starts and ends at the same depot. The open route is referred to as a route where the vehicle starts at a depot then services some customers and can end at the customer location. The half open route refers to a route where a vehicle starts at a depot and ends at an arbitrary depot after servicing the last customer (i.e., the same depot or the other depot). For conciseness and comprehensibility, three types of routes including the close, open, and half open are shown in Figure 1. Li et al. [12] proposed the MDVRP with time windows under shared depot resources (similar to the half open route), the goal was to minimize total transportation cost and the results showed that the flexible choice of the end depot can further to reduce total transportation cost. To sum up, the goal of our paper is to discuss the influence on the transportation risk and the transportation cost by considering heterogeneous vehicles and the half open route for hazmat MDVRP. Hazmat vehicle routing programming and risk evaluation approaches are mainly researched to reduce the hazmat transportation risk [3]. However, there are only a few works focused on the multi-depot vehicle routing problem (MDVRP) for hazmat transportation, especially the variants of the MDVRP based on real-world applications. The vehicle routing problem (VRP) is a classic combinational optimization problem for logistics transportation that aims to obtain the optimal route with the minimum cost of the constraints on capacities and customer demand. Due to the complexity of real-world transportation problems, there are many variants of the classical VRP in previous literature such as the heterogeneous VRP [4], the dynamic VRP [5], the MDVRP [6], the open VRP [7], etc. Similarly, there is some literature that has studied VRP and its variants for hazmat transportation in recent years. For instance, Du et al. [8] proposed a MDVRP for hazmat transportation to minimize total transportation risk, in which each customer can be serviced by multiple depots of the constraints on vehicles capacities and customer demand. Furthermore, Du et al. [3] formulated a multi-level programming model that considering multiple factors in real-world applications for hazmat MDVRP. Ma et al. [9] focused on the MDVRP for hazmat transportation considering energy consumption and formulated a multi-objective optimization model to minimize total transportation energy consumption and transportation risk. In most VRPs, vehicles are assumed to be homogeneous but this is unrealistic in real-world transportation problems. Moreover, some researchers considered heterogeneous VRP for hazmat transportation, in which vehicles have different characteristics (capacities, fixed cost, etc.). Bula et al. [10] focused on the heterogeneous VRP for hazmat transportation to determine the set of routes with minimizing total transportation risk. For economic reasons, Bula et al. [11] proposed a bi-objective programming model for hazmat heterogeneous VRP that simultaneously minimized total routing risk and total transportation cost. Wang et al. [1] studied the multi-depot heterogeneous VRP for hazmat with considering multiple depots, heterogeneous vehicles, and multiple accident scenarios. To the best of our knowledge, MDVRP with the half open route and heterogeneous vehicles for hazmat transportation is rarely studied. Thus, we simultaneously consider multiple depots, heterogeneous vehicles, half open route with minimizing total transportation risk and total transportation cost. The half open route is different from the close route and the open route. The close route is referred to as a route where the vehicle starts and ends at the same depot. The open route is referred to as a route where the vehicle starts at a depot then services some customers and can end at the customer location. The half open route refers to a route where a vehicle starts at a depot and ends at an arbitrary depot after servicing the last customer (i.e., the same depot or the other depot). For conciseness and comprehensibility, three types of routes including the close, open, and half open are shown in Figure 1. Li et al. [12] proposed the MDVRP with time windows under shared depot resources (similar to the half open route), the goal was to minimize total transportation cost and the results showed that the flexible choice of the end depot can further to reduce total transportation cost. To sum up, the goal of our paper is to discuss the influence on the transportation risk and the transportation cost by considering heterogeneous vehicles and the half open route for hazmat MDVRP.   For modeling the hazmat half open multi-depot heterogeneous VRP, the most important content is the risk evaluation approach. Generally speaking, hazmat transportation risk is measured as accident probability and accident consequence. Erkut and Ingolfsson [13] reviewed some risk evaluation approaches for hazmat transportation: the traditional risk model [14], the population exposure model [15], the incident probability model [16], the perceived risk model [17], the conditional risk model [18], the mean-variance model, the disutility model, and the minimax model [19], etc. The most frequently used risk evaluation approach is the expected risk model, in which accident probability based on historical data and accident consequence can be measured by the number of affected people, property damage, and environmental damage [20]. One of the common methods to measure accident consequence is the number of affected people in the neighborhood. Garrido and Bronfman [21] proposed that the number of affected people to be calculated by the affected area and surrounding population density, in which the affected area was determined by the exposure radius. Due to the consideration of heterogeneous vehicles, vehicle types are also influential for accident probability and the exposure radius. In addition, most researchers ignored the influence of the variation of vehicle loading on the exposure radius. They assumed the exposure radius was a constant value on each arc. However, vehicle loading will be reduced after servicing a customer, namely, the exposure radius is non-fixed on each arc. In summary, accident probability based on different types of vehicles and the exposure radius based on the variation of vehicle loading, vehicle types, and hazmat category are considered. Furthermore, decision makers generally consider not only reducing the total transportation risk but also the total transportation cost. Thus, the half open multi-depot heterogeneous VRP for hazmat transportation is formulated as a bi-objective mixed integer programming model with minimizing total transportation risk and total transportation cost.
For solution methods, the MDVRP is also an NP-hard problem as well as the VRP. There are many heuristic and meta-heuristic algorithms to solve this problem such as the hybrid genetic algorithm (GA) [22,23], the adaptive large neighborhood search [24,25], the artificial bee colony algorithm [26], the memetic algorithm [27], the two-phase algorithm consisting of the particle swarm optimization algorithm and the tabu search algorithm [28], among others. In real-world transportation problems, transportation risk and transportation cost objectives are generally in conflict [11]. Obviously, when there are multiple objectives, one optimal solution is best on one objective but worst on the others as the conflicts between objectives and incomparable phenomena. For example, the fewer vehicles are used, the lower the transportation cost. However, the reduction of the vehicle number will lead to an increase of unit vehicle loading. Therefore, we can only obtain the Pareto optimal solutions which are also called trade-off solutions for our bi-objective mixed integer programming model to provide more choices for decision makers. There is a large amount of literature about multi-objective optimization algorithms, such as the weighted-sum method, the ε-constraint method, the goal attainment approach, and meta-heuristic algorithms [29]. The ε-constraint method aims to minimize one of the multiple objectives and to limit the others by some limited values which can be used when the solver is exact or heuristic [30]. GA is a classic meta-heuristic algorithm that was introduced by [31]. It can effectively obtain optimal or near-optimal solutions due to its simplicity, easy operation and great flexibility [22]. Thus, the combination between the traditional ε-constraint method and the GA is reasonable and feasible. A hybrid intelligent algorithm based on the ε-constraint method and the GA is designed to eliminate the complexity of the proposed model and avoid significant disadvantages of the weighted-sum method.
In this paper, a corresponding expansion based on the traditional MDVRP for hazmat transportation is studied. The half open route in the relationship between vehicles and depots is considered, and the condition that each vehicle must end at the same depot after servicing the last customer is relaxed. In addition, the hazmat company has a heterogeneous fleet to satisfy customers demand for multiple types of hazmat. Furthermore, comparative results of different problems show that considering the half open route and heterogeneous vehicles for hazmat transportation makes a better balance between risk and cost. The main contributions of this paper are three aspects. Firstly, it presents a new variation of MDVRP for hazmat transportation as the half open multi-depot heterogeneous VRP that considers the influence of the half open route and heterogeneous vehicles. Then, it formulates the proposed problem as a bi-objective mixed integer programming model while minimizing transportation risk and transportation cost. Secondly, it proposes a novel loading-dependent transportation risk model considering the variation of vehicle loading, vehicle types, and hazmat category. Thirdly, it develops a hybrid intelligent algorithm based on the ε-constraint method and an improved GA to obtain Pareto optimal solutions. The effectiveness of the proposed model is demonstrated by numerical experiments.
The remainder of this paper is designed as follows. Section 2 defines a loadingdependent risk model and then formulates a bi-objective mixed integer programming model for minimizing transportation risk and transportation cost. A hybrid intelligent algorithm based on the ε-constraint method and an improved GA is designed in Section 3. Section 4 presents numerical experiments to show the effectiveness of the proposed model and algorithm. Finally, Section 5 provides concluding remarks and briefly discusses future work.

Problem Description
We consider a half open multi-depot heterogeneous VRP composed of multiple depots, multiple customers, and multiple types of vehicles. All vehicles are assigned based on the half open route which means a vehicle starts at a depot and ends at an arbitrary depot (the same depot or others) after servicing customers. An illustrative example with 3 depots and 10 customers is given in Figure 2, where a vehicle starts at depot D1 and ends at depot D1 after servicing customers 1, 2, and 3. In addition, a vehicle starts at depot D2 and ends at depot D1 after servicing customers 9 and 10, etc. The problem is formulated as a bi-objective mixed integer programming model, in which the first objective is minimizing transportation risk and the second objective is minimizing transportation cost. The goal is to determine the vehicle allocation and the optimal route with the constraints of vehicle capacities and customer demand. After obtaining the Pareto optimal solutions by simultaneously minimizing two objectives, decision makers can choose one of them based on their wishes. Some necessary assumptions are listed as follows:

•
Each vehicle starts at a depot and ends at an arbitrary depot after servicing customers. • Each customer can be serviced by multiple depots and multiple types of vehicles. • Each customer can be serviced once and only once by a vehicle. Then, it formulates the proposed problem as a bi-objective mixed integer programming model while minimizing transportation risk and transportation cost. Secondly, it proposes a novel loading-dependent transportation risk model considering the variation of vehicle loading, vehicle types, and hazmat category. Thirdly, it develops a hybrid intelligent algorithm based on the ε-constraint method and an improved GA to obtain Pareto optimal solutions. The effectiveness of the proposed model is demonstrated by numerical experiments. The remainder of this paper is designed as follows. Section 2 defines a loading-dependent risk model and then formulates a bi-objective mixed integer programming model for minimizing transportation risk and transportation cost. A hybrid intelligent algorithm based on the ε-constraint method and an improved GA is designed in Section 3. Section 4 presents numerical experiments to show the effectiveness of the proposed model and algorithm. Finally, Section 5 provides concluding remarks and briefly discusses future work.

Problem Description
We consider a half open multi-depot heterogeneous VRP composed of multiple depots, multiple customers, and multiple types of vehicles. All vehicles are assigned based on the half open route which means a vehicle starts at a depot and ends at an arbitrary depot (the same depot or others) after servicing customers. An illustrative example with 3 depots and 10 customers is given in Figure 2, where a vehicle starts at depot D1 and ends at depot D1 after servicing customers 1, 2, and 3. In addition, a vehicle starts at depot D2 and ends at depot D1 after servicing customers 9 and 10, etc. The problem is formulated as a bi-objective mixed integer programming model, in which the first objective is minimizing transportation risk and the second objective is minimizing transportation cost. The goal is to determine the vehicle allocation and the optimal route with the constraints of vehicle capacities and customer demand. After obtaining the Pareto optimal solutions by simultaneously minimizing two objectives, decision makers can choose one of them based on their wishes. Some necessary assumptions are listed as follows:  Each vehicle starts at a depot and ends at an arbitrary depot after servicing customers.  Each customer can be serviced by multiple depots and multiple types of vehicles.  Each customer can be serviced once and only once by a vehicle. For conciseness and comprehensibility, the notations we used are listed in Table 1.

Model Formulation
The transportation risk is a measure of accident probability and accident consequence, in which the accident consequence is usually measured as the number of exposed people. Moreover, the number of affected people around an accident is measured by the affected area and population density [32]. The population density is always assumed as constant values, and the affected area based on the exposure radius of an accident is formulated as follows: where p s ij and r sk ij are historical data that depend on the vehicle types. Since part of the hazmat will be unloaded when the vehicle arrives at each customer, the transportation risk is dynamic. Thus, we should consider the loading-dependent risk model and assume the function relationship between the vehicle loading and the exposure radius as follows: where α and β are constant values based on vehicle types and hazmat category, in which α is the ratio of the maximum exposure radius to maximum vehicle loading, the unit of α is km/ton. Apparently, the maximum vehicle loading and maximum exposure radius are different for different types of vehicles. In addition, β ≥ 1, and the exposure radius is linearly dependent on the vehicle loading when β = 1. For example, the sum of exposure radius of two vehicles for loading 1 ton is equal to the exposure radius of a vehicle for loading 2 tons when these vehicles have accidents. However, the exposure radius of a vehicle for loading 2 tons is bigger than two vehicles for loading 1 ton in more cases, thus β > 1 is more realistic. Similar to the literature [10,11,33], the transportation risk on arc (i, j) can be formulated as: The transportation risk on arc (i, j) as shown in yellow rectangle of Figure 3 is measured by accident probability which is based on vehicle type, the exposed people of the rectangular affected area which is depending on the product of the length of the arc and the exposure radius and the population density.
is km/ton. Apparently, the maximum vehicle loading and maximum exposure radius are different for different types of vehicles. In addition, 1   , and the exposure radius is linearly dependent on the vehicle loading when 1   . For example, the sum of exposure radius of two vehicles for loading 1 ton is equal to the exposure radius of a vehicle for loading 2 tons when these vehicles have accidents. However, the exposure radius of a vehicle for loading 2 tons is bigger than two vehicles for loading 1 ton in more cases, thus 1   is more realistic. Similar to the literature [10,11,33], the transportation risk on arc   , i j can be formulated as: The transportation risk on arc   , i j as shown in yellow rectangle of Figure 3 is measured by accident probability which is based on vehicle type, the exposed people of the rectangular affected area which is depending on the product of the length of the arc and the exposure radius and the population density. In the half open multi-depot heterogeneous VRP for hazmat transportation, transportation cost is usually considered by decision makers due to economic reasons. The transportation cost consists of fixed cost and routing cost when the vehicle k of type s travels on arc   , i j , which is defined by:

A Bi-Objective Mixed Integer Programming Model
Based on the aforementioned descriptions of assumptions, notations, risk, and cost models, we formulated the half open multi-depot heterogeneous VRP for hazmat transportation as the following bi-objective mixed integer programming model: In the half open multi-depot heterogeneous VRP for hazmat transportation, transportation cost is usually considered by decision makers due to economic reasons. The transportation cost consists of fixed cost and routing cost when the vehicle k of type s travels on arc (i, j), which is defined by:

A Bi-Objective Mixed Integer Programming Model
Based on the aforementioned descriptions of assumptions, notations, risk, and cost models, we formulated the half open multi-depot heterogeneous VRP for hazmat transportation as the following bi-objective mixed integer programming model: ∑ i∈V x sk ∑ i∈V ∑ j∈V y sk ij ≤ q s z sk , ∀s ∈ S, k ∈ K x sk ij = z sk = {0, 1}, ∀i, j ∈ V, s ∈ S, k ∈ K There are two objective functions of the mathematical model. The first objective function (5) minimizes the transportation risk. The second one (6) minimizes the transportation cost. Constraints (7)-(8) guarantee that each customer must be visited once and only once by a vehicle. Constraint (9) is that a vehicle starts at a depot and ends at an arbitrary depot. Constraint (10) is a flow balance for each node. Constraint (11) ensures sub-tour elimination. Constraint (12) states that customer demand should be satisfied. Constraint (13) denotes that vehicle loading must not exceed the maximum vehicle capacity. Constraint (14) indicates the relationship between two decision variables. Constraints (15)-(17) define three decision variables.

Algorithm
As aforementioned, since VRP is an NP-hard optimization problem, the above MDVRP is certainly an NP-hard problem. Due to the complexity of the proposed model and objectives attributes, there are many heuristic and meta-heuristic algorithms to solve this problem. One of the most frequently used meta-heuristic algorithms in the literature is the GA. GA is a widely used method of networks, as well as VRP. However, the bi-objective mixed integer programming model cannot be easily solved by the classic GA; thus we use the ε-constraint method by transforming the bi-objective problem into minimizing one of the two objectives and limiting another objective. Then, a hybrid intelligent algorithm composed of the ε-constraint method and an improved GA is designed. In the following, the principles of the ε-constraint method and the GA are developed.

Principles of ε-Constraint Method
A multi-objective optimization problem can be formulated as follows: where F 1 (x), F 2 (x), F 3 (x), . . . , F n (x) are n objective functions of a multi-objective optimization problem, and x is the decision variable belonging to the set of feasible solutions B.
For two feasible solutions x ∈ B and x * ∈ B, x dominates x * , if F i (x) ≤ F i (x * ) for all i ∈ {1, 2, . . . , n} with at least one strict inequality. A feasible solution x * ∈ B is said to be Pareto optimal solution (non-dominated solution), if and only if no other solution x ∈ B exists such that x dominates x * . The non-dominated set of the entire feasible decision space is called the Pareto optimal set and the boundary defined by the set of all points mapped from the Pareto optimal set is called the Pareto front.
The proposed model has two objective functions subject to a set of constraints. According to the precedence relationship of decision variables, hazmat can be transported on an arc only when one type of vehicle is used. Therefore, the risk objective function is regarded as the constraint on the cost objective function. Then, the ε-constraint method for solving our model is developed as follows: (1) Transform a bi-objective mixed integer programming model into a single objective model with limiting another single objective model, then we determine problem P which is minimizing transportation cost with transportation risk and other constraints, Problem P : constraints (7) where ε 2 is an upper limit of the value of risk function. Take problem P as an instance, the range of ε 2 is determined by the ideal point and the nadir point [34], the ideal point (Risk I , Cost I ) of our problem can be respectively obtained by solving every single objective model, the problem P 1 , i.e., minimize the transportation risk model with constraints (7)-(16) and the problem P 2 , i.e., minimize transportation cost model with constraints (7)- (16). Furthermore, the nadir point (Risk N , Cost N ) can be respectively obtained by using the optimal solution of cost function as an additional constraint of risk function with constraints (7)-(16) namely problem P 3 , similarly, using the optimal solution of risk function as an additional constraint of cost function with constraints (7)-(16) namely problem P 4 .
where ε 2 is changed with different m. (3) Solve problem P with an initial ε 2 to compute an optimal solution, then add the optimal solution to problem P to obtain optimal risk and cost vector. Then, repeat to solve problem P with different ε 2 and above steps to obtain the Pareto front.

Principles of GA
An improved GA is designed to solve every single objective function P 1 -P 4 based on the process of the ε-constraint method, taking problem P 1 as an instance to develop the process of the GA and the hybrid intelligent algorithm.

Representation Structure
An optimal solution for problem P 1 consists of determining optimal routes and vehicle types. Furthermore, the half open multi-depot heterogeneous VRP for hazmat transportation exists some problems to be determined: (1) which type of vehicle is assigned to service customers; (2) at which depot should vehicle start and end; (3) which route of vehicle to service customers; and (4) whether satisfying vehicles capacities and customer demand. Notably, a solution represents a chromosome. For the convenience of consideration, we assume that there are 3 depots, 10 customers, and 3 types of vehicles.
For instance, as shown in Figure 4, the chromosome structure is (D1, 1, 10, 9, D2; D3, 4, 5, 8, D2; D2, 6, 7, 3, 2, D3), in which red squares represent depots and yellow squares represent customers. One type of vehicle starts at depot D1 and services customers 1, 10, and 9 in order, then ends at depot D2, etc. Simply speaking, a vehicle starts at a depot to service several customers with satisfying vehicle capacity and customer demand, then can end at an arbitrary depot. If all conditions are satisfied, the optimal routes consist of customers, depots, and vehicle types are determined.
where 2  is changed with different m .
(3) Solve problem ' P with an initial 2  to compute an optimal solution, then add the optimal solution to problem P to obtain optimal risk and cost vector. Then, repeat to solve problem ' P with different 2  and above steps to obtain the Pareto front.

Principles of GA
An improved GA is designed to solve every single objective function 1 4 P P  based on the process of the  -constraint method, taking problem 1 P as an instance to develop the process of the GA and the hybrid intelligent algorithm.

Representation Structure
An optimal solution for problem 1 P consists of determining optimal routes and vehicle types. Furthermore, the half open multi-depot heterogeneous VRP for hazmat transportation exists some problems to be determined: (1) which type of vehicle is assigned to service customers; (2) at which depot should vehicle start and end; (3) which route of vehicle to service customers; and (4) whether satisfying vehicles capacities and customer demand. Notably, a solution represents a chromosome. For the convenience of consideration, we assume that there are 3 depots, 10 customers, and 3 types of vehicles. For instance, as shown in Figure 4, the chromosome structure is (D1, 1, 10, 9, D2; D3, 4, 5, 8, D2; D2, 6, 7, 3, 2, D3), in which red squares represent depots and yellow squares represent customers. One type of vehicle starts at depot D1 and services customers 1, 10, and 9 in order, then ends at depot D2 , etc. Simply speaking, a vehicle starts at a depot to service several customers with satisfying vehicle capacity and customer demand, then can end at an arbitrary depot. If all conditions are satisfied, the optimal routes consist of customers, depots, and vehicle types are determined.

Initialization Process
First, randomly sort arrays 1, 2, 10 that represent customers. Then, randomly select a number from 1 to 10 as the number of assigned vehicles. Here, we assume that 3 vehicles are assigned. Then, generate 3 empty cells and place one of 10 numbers randomly sorted customers in each empty cell in order. The remaining 7 customers are randomly placed in 1 of 3 cells. Moreover, according to the assigned 3 vehicles, 6 depots are randomly generated. Insert 3 cells with the customer representation between every 2 depots. Finally, 3 cells are rearranged to form a complete chromosome. Repeat the above procedures _ pop size times to obtain _ pop size chromosomes.

Fitness Function
Problem 1 P is minimizing the transportation cost, the fitness function equals to the

Initialization Process
First, randomly sort arrays 1, 2, 10 that represent customers. Then, randomly select a number from 1 to 10 as the number of assigned vehicles. Here, we assume that 3 vehicles are assigned. Then, generate 3 empty cells and place one of 10 numbers randomly sorted customers in each empty cell in order. The remaining 7 customers are randomly placed in 1 of 3 cells. Moreover, according to the assigned 3 vehicles, 6 depots are randomly generated. Insert 3 cells with the customer representation between every 2 depots. Finally, 3 cells are rearranged to form a complete chromosome. Repeat the above procedures pop_size times to obtain pop_size chromosomes.

Fitness Function
Problem P 1 is minimizing the transportation cost, the fitness function equals to the reciprocal of the cost objective. The chromosome with a higher fitness value has a bigger chance of being selected. Then, normalize the fitness value to calculate the cumulative probability. Repeat the above procedures pop_size times.

Selection Process
We employed the method of spinning the roulette wheel to select chromosomes. First, according to the cumulative probability, select two individuals from the population. Then, treat these two chromosomes as parents to crossover and mutation. The bigger the fitness value, the bigger the selection probability. Repeat pop_size times to obtain new chromosomes.

Crossover Process
The partial matching crossover method is used. Denote p c as the crossover probability. Divide chromosomes into pairs. Select two chromosomes from parents, then randomly generate two integers between customers 1 and 10 to determine two positions. Then, crossover the number in the middle of the two positions. Eliminate conflicts by using a partial mapping method for repeated numbers in an individual, without repeating digital reservations. Finally, reinsert customers into their original locations and obtain two children. Repeat the above process pop_size times to get pop_size chromosomes. The crossover process is shown in Figure 5.

Selection Process
We employed the method of spinning the roulette wheel to select chromosomes. First, according to the cumulative probability, select two individuals from the population. Then, treat these two chromosomes as parents to crossover and mutation. The bigger the fitness value, the bigger the selection probability. Repeat _ pop size times to obtain new chromosomes.

Crossover Process
The partial matching crossover method is used. Denote c p as the crossover probability. Divide chromosomes into pairs. Select two chromosomes from parents, then randomly generate two integers between customers 1 and 10 to determine two positions. Then, crossover the number in the middle of the two positions. Eliminate conflicts by using a partial mapping method for repeated numbers in an individual, without repeating digital reservations. Finally, reinsert customers into their original locations and obtain two children. Repeat the above process _ pop size times to get _ pop size chromosomes. The crossover process is shown in Figure 5.

Mutation Process
Denote m p as the mutation probability and repeat the following process _ pop size times. Select one chromosome from the parents, randomly generate two integers between customers 1 and 10 from one child chromosome obtained in the previous step to determine two positions, then swap these two numbers in each chromosome. The mutation process is shown in Figure 6.

Mutation Process
Denote p m as the mutation probability and repeat the following process pop_size times. Select one chromosome from the parents, randomly generate two integers between customers 1 and 10 from one child chromosome obtained in the previous step to determine Sustainability 2021, 13, 1262 10 of 17 two positions, then swap these two numbers in each chromosome. The mutation process is shown in Figure 6.

Mutation Process
Denote m p as the mutation probability and repeat the following process _ pop size times. Select one chromosome from the parents, randomly generate two integers between customers 1 and 10 from one child chromosome obtained in the previous step to determine two positions, then swap these two numbers in each chromosome. The mutation process is shown in Figure 6.

The Hybrid Intelligent Algorithm
In Section 3.1, the ε-constraint method was used to solve the proposed model and principles of the ε-constraint method were presented. Moreove, for every single objective function, the GA was designed to obtain optimal solutions. Thus, a hybrid intelligent algorithm is proposed to solve the bi-objective mixed integer model, and the complete algorithm steps are developed as follows: Step 1. Compute problem P 1 , P 2 , P 3 , and P 4 by the GA.
Step 2. Determine the range (Risk I , Risk N ) of the risk function. Step 4. Solve the problem P with ε 2 , and compute an optimal solution x * by the GA.
Step 5. Solve the risk function and the cost function with x * , obtain the set of values {Risk(x * ), Cost(x * )}.
Step 6. Repeat to solve problem P with different ε 2 by the GA. If ε 2 < Risk I , go to step 4.
Step 7. Obtain the Pareto front.

Dataset and Experimental Setup
We conduct numerical experiments on the bi-objective mixed integer programming model of the half open multi-depot heterogeneous VRP for hazmat transportation, which is associated with 3 depots, 10 customers, and 3 types of unlimited vehicles. As shown in Figure 7, red squares and yellow circles respectively represent depots and customers, i.e., their locations are predetermined and the distance (km) between customers and depots are calculated by the Euclidean distance. The quantities of customer demand are randomly given in customer order, namely, (3, 2, 1, 4, 5, 6, 5, 3, 2, 4) (ton). Moreover, population density is also randomly generated. Vehicle attributes are presented in Table 2, in which only one category of hazmat is considered, and the relevant parameters of the GA are as follows: p c = 0.9, p m = 0.09, pop_size = 100 and iterations = 1000. Furthermore, to illustrate the effectiveness of the proposed model, we compare the half open multi-depot heterogeneous VRP with the close multi-depot heterogeneous VRP and the half open multi-depot homogeneous VRP. The above algorithm was run with the personal computer environment: Inter (R) Core (TM) i5-5350U CPU, 1.80GHz, 8.00 GB, and Windows 10 operating system. Sustainability 2021, 13, x FOR PEER REVIEW 12 of 18  Since our algorithm ends up with Pareto optimal solutions, how to select a Pareto optimal vector to compare is the next significant step. Here, we adopt a method to select a Pareto optimal vector of all Pareto optimal solutions which is: To demonstrate whether it is beneficial to adopt the half open route for hazmat transportation, we compare the optimal solutions obtained from the half open multi-depot heterogeneous VRP and the close multi-depot heterogeneous VRP. The close multi-depot heterogeneous VRP is considered a special situation differing from the half open multi-depot heterogeneous VRP. To obtain Pareto optimal solutions, the bi-objective mixed integer programming model for the close multi-depot heterogeneous VRP will have a little change compared to the proposed model. The constraint (9)

 
. In other words, a vehicle must start and end at  Since our algorithm ends up with Pareto optimal solutions, how to select a Pareto optimal vector to compare is the next significant step. Here, we adopt a method to select a Pareto optimal vector of all Pareto optimal solutions which is: (1) Calculate new values of risk R and cost C with different Pareto optimal vectors, in which R = R−R min R max −R min , C = C−C min C max −C min , where R and C are different for each Pareto optimal vector; (2) Predefine weight coefficients λ 1 and λ 2 , set λ 1 = λ 2 = 0.5 for convenience; (3) Calculate the values of minλ 1 R + λ 2 C and choose a Pareto optimal vector.
To demonstrate whether it is beneficial to adopt the half open route for hazmat transportation, we compare the optimal solutions obtained from the half open multi-depot heterogeneous VRP and the close multi-depot heterogeneous VRP. The close multi-depot heterogeneous VRP is considered a special situation differing from the half open multidepot heterogeneous VRP. To obtain Pareto optimal solutions, the bi-objective mixed integer programming model for the close multi-depot heterogeneous VRP will have a little change compared to the proposed model. The constraint (9) in the proposed model ∑ i∈I ∑ j∈J x sk ij = ∑ j∈J ∑ i∈I x sk ji ≤ 1, ∀s ∈ S, k ∈ K, which means that a vehicle starts at a depot and ends at an arbitrary depot will be translated to ∑ i∈I ∑ j∈J x sk ij = ∑ j∈J ∑ i∈I x sk ji = 1, ∀s ∈ S, k ∈ K. In other words, a vehicle must start and end at the same depot in this situation. Then, we compare the results of the half open route with the close route in two ways: (1) Calculate the difference values of all Pareto optimal solutions for the half open route and the close route when M = 20, and 15 different Pareto optimal vectors are obtained. As shown in Figure 8, the Pareto fronts of two problems are obtained. It is clear that all Pareto optimal solutions of the half open route are better than the close route. The results show that the average risk can be reduced by 3.99% and the average cost can be reduced by 2.01%, namely, it is better to adopt the half open route. (2) Choose a Pareto optimal vector using the aforementioned method. As shown in Table 3, in this situation, risk can be reduced by 3.55% and cost can be reduced by 2.36%. For the close route, decision makers may wish a vehicle can carry more hazmat but will cause more risk and cost. For convenience, the half open optimal route is shown in Figure 9.  (2) Choose a Pareto optimal vector using the aforementioned method. As shown in Table 3, in this situation, risk can be reduced by 3.55% and cost can be reduced by 2.36%. For the close route, decision makers may wish a vehicle can carry more hazmat but will cause more risk and cost. For convenience, the half open optimal route is shown in Figure 9.  For the homogeneous VRP, only one type of vehicle can be used, including light vehicle, medium vehicle, and heavy vehicle. The bi-objective mixed integer programming model will have a big change which is variables x sk ij , y sk ij , and z sk change to x k ij , y k ij , and z k . In addition, the transportation risk model (5), cost model (6), and constraints (7) z . In addition, the transportation risk model (5), cost model (6), and constraints (7)-(16) have the same change. Then, use the same method to compare the Pareto optimal vectors of two problems.
As shown in Table 4, firstly, compare the results of using heterogeneous vehicles with the light vehicle, the risk can be increased by 1.12% and the cost can be reduced by 3.99%. The results indicate that when the hazmat company only has one type of the light vehicle, more vehicles are assigned to service customers and will create more cost. Though the risk is reduced, the increase of cost is unacceptable. Secondly, compare the results of using heterogeneous vehicles with the medium vehicle, the risk can be increased by 1.06% and the cost can be reduced by 4.56%. The results indicate that only using one type of the medium vehicle causes more cost than heterogeneous vehicles, due to the medium vehicle having a bigger fixed cost and unit distance transportation cost than the light vehicle. Thirdly, compare the results of using heterogeneous vehicles with the heavy vehicle, the risk can be reduced by 2.48% and the cost can be reduced by 5.29%. The results indicate that only using one type of the heavy vehicle, both risk and cost are increased because the heavy vehicle has the biggest fixed cost, unit distance transportation cost, and carries more hazmat. To sum up, the use of heterogeneous vehicles reduced the cost by 3.99% and 4.56% in terms of acceptable risk compared to light vehicle and medium vehicle. Moreover, using heterogeneous vehicles can reduce risk and cost by 2.48% and 5.29%, respectively, compared to heavy vehicle. Thus, considering heterogeneous vehicles is more feasible to be implemented in real-world applications. As shown in Table 4, firstly, compare the results of using heterogeneous vehicles with the light vehicle, the risk can be increased by 1.12% and the cost can be reduced by 3.99%.
The results indicate that when the hazmat company only has one type of the light vehicle, more vehicles are assigned to service customers and will create more cost. Though the risk is reduced, the increase of cost is unacceptable. Secondly, compare the results of using heterogeneous vehicles with the medium vehicle, the risk can be increased by 1.06% and the cost can be reduced by 4.56%. The results indicate that only using one type of the medium vehicle causes more cost than heterogeneous vehicles, due to the medium vehicle having a bigger fixed cost and unit distance transportation cost than the light vehicle. Thirdly, compare the results of using heterogeneous vehicles with the heavy vehicle, the risk can be reduced by 2.48% and the cost can be reduced by 5.29%. The results indicate that only using one type of the heavy vehicle, both risk and cost are increased because the heavy vehicle has the biggest fixed cost, unit distance transportation cost, and carries more hazmat. To sum up, the use of heterogeneous vehicles reduced the cost by 3.99% and 4.56% in terms of acceptable risk compared to light vehicle and medium vehicle. Moreover, using heterogeneous vehicles can reduce risk and cost by 2.48% and 5.29%, respectively, compared to heavy vehicle. Thus, considering heterogeneous vehicles is more feasible to be implemented in real-world applications.

Case 3: Comparison of Parameters
Different initial parameters are generated from the perspective of fixed cost and accident probability of vehicles to demonstrate the effectiveness of the proposed model. As shown in Table 5, instances 2, 3, and 4 have different accident probability from instance 1. For instance 2, the risk can be reduced by 1.90% and the cost can be increased by 1.43%. The results show that the accident probability of different types of vehicles will influence the number of vehicles used, but the consideration of the fixed cost of vehicles is not important. Meanwhile, for instances 3 and 4, the risk can be respectively reduced by 2.55% and 3.89%. In addition, the cost can be increased by 5.71% and 4.92%. The results indicate that when the accident probability changes, the risk is the first objective to minimize and cost may increase. For instances 5, 6, and 7, the fixed cost of different types of vehicles is changed. For instance 5, the cost can be reduced by 1.10% and the risk can be reduced by 1.53%. The results show that the heavy vehicle may be more used. For instance 6, the risk can be increased by 0.95% which is acceptable, and the cost can be increased by 2.17% because of the increase in the fixed cost. For instance 7, the risk can be reduced by 1.02% and the cost can be increased by 2.53%. To sum up, the fixed cost also influences the vehicles used and the total transportation cost. Moreover, the change of risk is acceptable. Comparing instances 8, 9, and 10, accident probability and fixed cost will influence risk and cost. In instances 11, 12, and 13, when the interval gets bigger of accident probability and fixed cost we obtain the same results. In a word, the fixed cost and accident probability of vehicles are two important factors affecting the risk and cost of hazmat transportation. Therefore, how to choose the vehicle to use and how to arrange the vehicle route according to the actual customers demand are the main issues for decision makers to consider. Furthermore, we conclude that risk and cost are two conflict objectives. Once decision makers want less risk, it will cause more cost, thus decision makers should choose one of Pareto optimal solutions with respective wishes.

Conclusions
In this study, a half open multi-depot heterogeneous VRP for hazmat transportation was studied that considered heterogeneous vehicles located at multiple depots, and each vehicle starting at a depot and ending at an arbitrary depot after servicing the last customer. The main contributions are listed as follows: first, we considered a real-world application that a hazmat company has multiple depots and heterogeneous vehicles to service customers, in which all vehicles were assigned based on the half open route. Furthermore, the transportation risk model was defined as the loading-dependent risk considering the variation of vehicle loading, vehicle types, and hazmat category. Then, we proposed a bi-objective mixed integer programming model to achieve a good balance between transportation risk and transportation cost. To solve the proposed model in an efficient way, we designed a hybrid intelligent algorithm composed of the ε-constraint method and an improved GA to obtain Pareto optimal solutions. By comparing the proposed problem with the close multi-depot heterogeneous VRP, the computational results showed that the average transportation risk and cost of the former can be reduced by 3.99% and 2.01%, respectively. Thus, the consideration of the half open route that applying on hazmat transportation is reasonable. Furthermore, the contrastive results of the proposed problem and the half open multi-depot homogeneous VRP indicated that only using light vehicles and medium vehicles increased the transportation cost by 3.99% and 4.56% with acceptable risk reduction, and only using heavy vehicles increased the transportation risk and cost by 2.48% and 5.29%. These numerical experiments showed that considering the half open route and heterogeneous vehicles for hazmat transportation by decision makers creates a better balance between risk and cost. This work is considered from the perspective of a hazmat company, but it can be extended to the perspective of multiple hazmat companies. These companies can reduce transportation risk and save transportation cost by sharing their own depots.
Further research is required to consider some realistic factors, such as traffic restriction and customer satisfaction, and improve the running efficiency of the algorithms. Moreover, we plan to focus on the MDVRP for hazmat transportation and its variants. The transportation risk also needs to be measured combined with data prediction method [35] and realistic traffic conditions [36]. Furthermore, more realistic situations should be taken into account, such as the multi-depot location routing problem, the multi-period MDVRP for hazmat transportation.