New Hybrid Algorithm Based on BicriterionAnt for Solving Multiobjective Green Vehicle Routing Problem

Emile Nawej Kayij, Joél Lema Makubikua, Justin Dupar Kampempe Busili

Faculty of Technological and Sciences, University of Kinshasa, Kinshasa, Democratic Republic of the Congo.

DOI: 10.4236/ajor.2023.133003                

Abstract

The main objective of this paper is to propose a new hybrid algorithm for solving the Bi objective green vehicle routing problem (BGVRP) from the BicriterionAnt metaheuristic. The methodology used is subdivided as follows: first, we introduce data from the GVRP or instances from the literature. Second, we use the first cluster route second technique using the k-means algorithm, then we apply the BicriterionAntAPE (BicriterionAnt Adjacent Pairwise Exchange) algorithm to each cluster obtained. And finally, we make a comparative analysis of the results obtained by the case study as well as instances from the literature with some existing metaheuristics NSGA, SPEA, BicriterionAnt in order to see the performance of the new hybrid algorithm. The results show that the routes which minimize the total distance traveled by the vehicles are different from those which minimize the CO₂ pollution, which can be understood by the fact that the objectives are conflicting. In this study, we also find that the optimal route reduces product CO₂ by almost 7.2% compared to the worst route.

Keywords

Metaheuristics, Green Vehicle Routing Problem, Ant Colony Algorithm, Genetic Algorithms, Green Logistics

1. Introduction

Green transportation issues are increasingly showing interest from theoretical, political and social perspectives. Transport activity has a high rate of negative effects on the environment, due to pollutants and the greenhouse effect. The consequences of its emissions are climate change and health complications, reasons why emissions from the transport sector must be reduced, the WHO reports that diesel fumes create several diseases such as lung cancer (International Agency for Cancer Research, 2012). Nevertheless, scientific works that deal with GVRP with the objective of minimizing CO₂ emissions are rare and new in the literature.

Indeed, the combustion of fossil fuels in internal combustion engines produces carbon dioxide CO₂, the propagation of which in the atmosphere contributes to the greenhouse effect and global warming. It should be noted that one liter of fuel produces approximately 2.4 kg of CO₂ during combustion, a value which depends on the type of fuel (petrol, diesel, natural gas ...) and its density.

CO₂ emissions are directly proportional to vehicle fuel consumption. [1] To quantify emissions, several calculation methods exist. Each method takes into account on average 6 to 14 different factors influencing the amount of GHG (Greenhouse Gas) emissions evaporated.

The classic VRP neglects current environmental concerns and deals primarily with economic return. It is therefore important for policy makers to apply restrictive regulations to control greenhouse gas emissions. All of these concerns have drawn the attention of researchers to the Green Vehicle Routing Problem (GVRP) [2] .

Some authors have proposed approaches and algorithms to solve the GVRP. In 2015 A. Moutaoukil et al. [3] worked on the minimization of CO₂ emissions in distribution circuits with heterogeneous fleet, its strength lies in modeling a vehicle routing problem using a heterogeneous fleet, but one of the limits of this proposal lies in the high computation times which only allow small instances to be taken into account. N. Rezaie et al. [4] in 2019 studied a problem of routing green vehicles with time windows considering the heterogeneous fleet of vehicles: two metaheuristic algorithms. The comparison showed that the proposed algorithms provide high-quality solutions with regard to objective functions. The paper has some limitations. First, it can be merely applied to products (or problems) that hard time windows are suitable for them, such as newspaper, meat, vegetable, etc. Second, the practitioners must have the knowledge to use the results of this paper in the real world. Júlio César Ferreira et al. [5] in 2021 proposed for a two-objective green problem a New hybrid optimization algorithm CWNSGA-II and CWTSNSGA-II applied to the distribution of a newspaper, he has the advantage of having obtained superior results for case studies and literature instances, but the difficulty in selecting the best route from a utility function that is often not easy to set up. In 2021 Joydeep Dutta et al. [6] proposed a hybrid multi-objective evolutionary algorithm for open green vehicle routing problem through cluster primary-route secondary approach. This model suggests a solution based on the managers “decision makers” choice from a set of alternative solution.

In this article, the GVRP has been studied, where a central repository provides the demand of customers using a heterogeneous fleet of vehicles. In addition, it was assumed that the rate of fuel consumption depends on the distance traveled and the load of the vehicles as well as the type of vehicle used.

To formulate the problem, a bi-objective model of combinatorial linear programming has been proposed [7] . Thus, the model simultaneously minimizes the total transport costs and the CO₂ emissions of the vehicles during a tour. In the literature, some authors have proven that even simple variants of VRP are NP-hard, like [8] [9] , not to mention the more complex variant of GVRP studied in this article.

In this paper, we make a hybridization of three heuristic approaches to solve the two-objective green vehicle routing problem. This new hybrid algorithm called BicriterionAntAPE is a general, globalizing, high-level hybridization that combines the strengths of each heuristic. The interaction between them is seen by the fact that the heuristics are launched one after the other, each taking the output produced by the previous one. In the literature, most authors use a genetic algorithm to solve the VRP. After finding interesting individual routes, routes are then improved using the 2-opt method followed by the 3-opt [10] . What is new in this article, however, is that the solution found by BicriterionAnt is improved by APE, which is easier to apply instead of the usual r-opt type methods.

Generally the VRP is an NP-hard routing problem defined on a graph $G=\left(V\mathrm{,}E\right)$ where V is a set of n nodes (customers) and E a set of m edges, consisting of a fleet vehicle of capacity Q based in a depot node. Each edge $\left(i\mathrm{,}j\right)$ has a traversal cost $c_{ij}\ge 0$ and a demand $q_i\ge 0$ from customer i as as shown in Figure 1.

The objective function can be either: the minimization of the number of vehicles used, the minimization of the total distance traveled by the vehicles, the

Figure 1. Graphical presentation of the VRP.

minimization of the total duration of the routes, the minimization of the total cost of the routes (taking into account the costs of the vehicles, drivers, the minimization of penalties related to violation of constraints, particularly in the case of time windows, the maximization of gains generated by rounds the objectives of minimizing the number of vehicles and the total distance (or duration) of rounds are conflicting because the reduction in the number of vehicles most often leads to an increase in the total distance travelled.

2. Position of the Problem

Solving the multi-objective vehicle routing problem (MOVRP) by so-called exact methods presents many difficulties for medium and large-sized instances. Ant colony algorithms solve this type of problem in a reasonable time and provide good quality approximate solutions [11] .

Ant colony-based algorithms were originally proposed in [12] [13] [14] . They have been successfully applied to several combinatorial problems like the traveling salesman problem, quadratic assignment, vehicle routing, multidimensional knapsack, etc. We then extend this study to the case of multi-objective optimization problems, where the choice of a phenomenal structure is even more difficult. We propose, in this work the resolution of a case study of the GVRP with two objectives to do this, it is essential to choose the appropriate phenomenal strategy, hence the choice of the BicriterionAnt (ACO) algorithm proposed by Iredi et al. [15] .

To calculate CO₂ emissions, we used the method proposed by Zhang et al. [16] . We therefore assumed two fuel consumption rates. The first is for the vehicle when unloaded (empty load) and the second is for the same vehicle when loaded to its maximum capacity. Then, taking into account the load and the distance traveled by the vehicles, the CO₂ emissions per each unit of distance will be determined.

3. Literature Review

The vehicle routing problem was first introduced in the literature by Dantzig and Ramser in 1959 [17] . They work on a study dealing with the truck dispatch problem in a vehicle routing problem (VRP) with capabilities and specific constraints to serve groups of customers in different geographical locations from a central warehouse. These vehicles may or may not return to the central depot.

As pointed out by Reimann [18] , the VRP and all its variants form an important problem in logistics and have been the subject of much research, such as recently by [19] [20] [21] . When the vehicle belongs to the company, the vehicle must return to the central depot, we then speak of VRP closed; otherwise we speak of VRP Open. Sariklis and Powell [22] first introduced the term Open Vehicle Route Problem (OVRP) this was followed by several published papers on OVRP problems using multiple methods/algorithms for their solution. Dutta, Barma, Kar [23] recently applied a modified version of Kruskal’s method and a hybrid algorithm of genetic algorithms to solve OVRP.

Gutierrez et al. [24] studied a new variant of VRP where time is a random variable. Eshtehadi et al. [25] in turn have explored another facet of uncertainty in slotted VRP; they took demand and travel times as uncertain parameters, they discussed three approaches to analyze uncertainty: soft worst case, hard worst case and chance conditions.

Kopper et al. [26] proposed a mathematical method to study heterogeneous GVRPs simply targeting CO₂ emissions, in this paper, different types of vehicles were taken into account as well as the consumption rates of different types of fuels used, payload and mileage.

Yang et al. [27] on the other hand proposed a multi-objective model to deal with GVRP; this mathematical model seeks to minimize operating costs and CO₂ emissions while maximizing the level of customer satisfaction, unlike Masmoudi et al. [28] they considered the fill function in the problem. However, they took into account a homogeneous fleet of vehicles. They first proposed a MILP (Mixed Integer Linear Problem) model of the problem, and then developed a two-step hybrid heuristic algorithm for a home care system where a patient may need the services of several specialists at the same time.

4. Multi Objective Programming

Most real world optimization problems are complex and one of the aspects contributing to this complexity is their multi-objective nature Durillo et al. [29] . In optimization works, including multi-objective optimization, the emphasis is on finding the global optimum or Pareto front yielding the best possible values of the Deb et al. [30] . For others, multi-objective (or multi-criteria) optimization is a process aimed at systematically and simultaneously optimizing a collection of objective functions Marler et al. [31] .

In general, a multi-objective optimization problem can be defined as:

$\min F\left(x\right)=\left(f_1\left(x\right)\mathrm{,}{f}_2\left(x\right)\mathrm{,}\cdots \mathrm{,}{f}_m\left(x\right)\right)$ avec $m\ge 2$

Under the constraints $g\left(x\right)\le 0$ and $h\left(x\right)=0$ where m is the number of functions to optimize, $x=\left(x_1\mathrm{,}\cdots \mathrm{,}{x}_n\right)$ is the vector of decision variables, $F\left(x\right)=\left(f_1\left(x\right)\mathrm{,}{f}_2\left(x\right)\mathrm{,}\cdots \mathrm{,}{f}_m\left(x\right)\right)$ is vector of the criteria to be optimized, $g\left(x\right)$ and $h\left(x\right)$ represent respectively m inequality constraints and p equality constraints. This set of constraints delimits a restricted space of feasible solutions that one will note D.

The image of a solution x in the criteria space is the point $y=\left(y_1\mathrm{,}\cdots \mathrm{,}{y}_m\right)$ with $y_i=f_i\left(x\right),i=1,\cdots ,m$ and $Y=F\left(D\right)$ representing the feasible points in the objective space.

We impose on this set a partial order relation called dominance relation that we will define in the next point.

4.1. Notion of Pareto Dominance

A solution $y=\left(y_1\mathrm{,}\cdots \mathrm{,}{y}_m\right)$ dominates a solution $z=\left(z_1\mathrm{,}\cdots \mathrm{,}{z}_m\right)$ if $\forall i\in \left\{\mathrm{1,}\cdots \mathrm{,}m\right\}$ , $f_i\left(y\right)\le f_i\left(z\right)$ and $\exists i\in \left\{\mathrm{1,}\cdots \mathrm{,}m\right\}$ $f_i\left(y\right)<f_i\left(z\right)$ .

If the solution y dominates the solution z it will be denoted $y\prec z$ . For any pair of solutions y and z one and only one of the following cases can occur: y dominates z or y is dominated by z or aigain y and z are equivalent in the sense of dominance.

Equivalent solutions in the sense of dominance are called in what follows, equivalent solutions in the sense of Pareto or Pareto equivalent solutions or again, non-dominated solutions.

4.2. Optimization by Ant Colonies for Multi-Objective Problems

According to Dorigo et al. [32] ACO algorithms have several advantages in their applications in different, of which several discrete optimization problems are successful. The majority of these problems are NP-hard and the application of ACO algorithms has allowed finding good quality solutions in a reasonable time.

Some variants of the first ACO algorithm named Ant System [12] based on the behavior of artificial ants such as the Max-Min Ant System, Ant Colony System, or Rank-Based Ant System algorithms have emerged, these extensions were implemented due to the fact that despite the encouraging results of Ant System, it was not competitive enough against state of the art algorithms, especially for large-scale problems (such as genetic algorithms).

ACOs are part of combinatorial optimization techniques that have been designed to solve multi-objective problems often having conflicting objectives. An overview of ACO algorithms allowing multi-objective problem solving is available in Garcia-Martinez et al. [33] . These algorithms are gathered under the name of MOACO (Multiple Objective ACO). We can distinguish two families of MOACO algorithms: Pareto and Non-Pareto. MOACO algorithms based on a Pareto evaluation method (Pareto-based MOACO) consist in generating solutions on the Pareto front by applying the principles of Pareto-dominance. The best known algorithms of this family are: CPACO [34] , MOAQ [35] , PACO-MO [36] , BicriterionAnt [37] , BicriterionMC [37] and ACO-bQAP [36] .

4.3. Ant Colony Algorithm

The ACO algorithm was proposed by Dorigo to choose the shortest path in the traveling salesman problem (TSP). It implements the use of ant trail pheromones when searching for the shortest path from their nest to the food source as shown in Figure 2. In this article the BicriterionAnt version particularly caught our attention to solve our problem of touring green vehicles with two objectives.

We first present the initial algorithm Ant Colony System (ACS) which allows finding the shortest Hamiltonian cycle in a graph.

4.4. Data and Rating

● Let $G=\left(V\mathrm{,}E\right)$ be an undirected graph.

● The finite set V of vertices represents customer.

● The finite set $E=\left\{\left(i\mathrm{,}j\right)\mathrm{<mo>\; /\; </mo>}i\mathrm{,}j\in V\right\}$ of edges linking the vertices together.

Figure 2. Overview of ACO.

● $d_{ij}$ represents the edge weights $\left(i\mathrm{,}j\right)\in E$ (the distance between vertex i and vertex j).

● $L^k\left(u\right)$ gives the length of the circuit u representing the total distance traveled during the tour by the ant k.

● ${\tau }_{ij}\left(t\right)$ is the value of the pheromone trail at time t on edge $\left(i\mathrm{,}j\right)$ .

● $n=\left|V\right|$ the total number of customer.

● ${\eta }_{ij}\left(t\right)=\frac{1}{d_{ij}}$ is the constant indicating the visibility of customer j from customer i.

4.5. Transition Choice

Each ant has a memory that will be emptied once the cycle is over, allowing it to avoid returning to a customer already visited. An ant k found in customer i at time t chooses customer j as its destination according to two parameters: the visibility ${\eta }_{ij}$ and the pheromone density ${\tau }_{ij}\left(t\right)$ of the edge connecting it, depending on the value of $p_{ij}^k\left(t\right)$ .

$p_{ij}^k\left(t\right)=\frac{{\left({\tau }_{ij}\left(t\right)\right)}^{\alpha }\cdot {\left({\eta }_{ij}\right)}^{\beta }}{{\displaystyle {\sum }_{j\in N_i}\left({\left({\tau }_{ij}\left(t\right)\right)}^{\alpha }\cdot {\left({\eta }_{ij}\right)}^{\beta }\right)}}\text{if}j\in N_i^k$ (1)

where $p_{ij}^k\left(t\right)$ is the probability that ant k chooses j from i. $N_i^k$ is the set of unvisited neighbors of vertex i by ant k in the current cycle, $\alpha$ and $\beta$ are parameters that control the importance that the ant gives to an edge (intensification/diversification).

4.6. Pheromone Update

At the end of each cycle (each ant has traversed the n vertices that make up the graph), the traces of the pheromones are updated according to the formula:

${\tau }_{ij}\left(t+1\right)=\left(1-\rho \right){\tau }_{ij}\left(t\right)+\Delta {\tau }_{ij}\left(t\right)$ (2)

where $\rho$ is an evaporation coefficient such that $0<\rho <1$ .

$\left(1-\rho \right)$ is the evaporation of the pheromone.

$\Delta {\tau }_{ij}\left(t\right)={\displaystyle \underset{k=1}{\overset{m}{\sum }}}\Delta {\tau }_{ij}^k\left(t\right)$ (3)

represents the quantity of pheromones per unit length deposited on the edge $\left(i\mathrm{,}j\right)$ by the k^th ant between t and $t+n$ , it is given by:

$\Delta {\tau }_{ij}\left(t\right)=\{\begin{array}{ll}\frac{Q}{L^k}\hfill & \text{if}\text{the}\text{ant}k\text{crosses}\text{edge}\left(i\mathrm{,}j\right)\text{in}\text{its}\text{turn}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$

with, $L^k$ : the length of the turn of the k^th ant.

Q: a constant positive number.

5. Presentation of the BicriterionAnt algorithm

The BicriterionAnt algorithm was proposed by Iredi et al. [15] especially to solve the bi-criteria vehicle routing problem. It uses two different pheromone trail structures, $\tau$ and ${\tau }^{\prime }$ , one for each criterion considered.

At each generation, each of the m ants in the colony generates a solution to the problem. During the solution construction step, the ant chooses the next vertex j to visit relative to the probability $p_{ij}^k(\; t\; )$

$p_{ij}^k\left(t\right)=\frac{{\left({\tau }_{ij}\left(t\right)\right)}^{\alpha {\lambda }_k}\cdot {\left({\eta }_{ij}\right)}^{\beta {\lambda }_k}\cdot {\left({{\tau }^{\prime }}_{ij}\left(t\right)\right)}^{\alpha \left(1-{\lambda }_k\right)}\cdot {\left({{\eta }^{\prime }}_{ij}\right)}^{\beta \left(1-{\lambda }_k\right)}}{{\displaystyle {\sum }_{j\in N_i}\left({\left({\tau }_{ij}\left(t\right)\right)}^{\alpha {\lambda }_k}\cdot {\left({\eta }_{ij}\right)}^{\beta {\lambda }_k}\cdot {\left({{\tau }^{\prime }}_{ij}\left(t\right)\right)}^{\alpha \left(1-{\lambda }_k\right)}\cdot {\left({{\eta }^{\prime }}_{ij}\right)}^{\beta \left(1-{\lambda }_k\right)}\right)}}\text{if}j\in N_i^k$ (4)

for each m ants calculate the exchange rule $\lambda$ as follows:

${\lambda }_k=\frac{k-1}{m-1}$ (5)

Once all the ants have built up their solutions, the pheromone trails are evaporated by the usual ACO rule. Since there are two variables for pheromone information and heuristic information, the pheromone evaporation process is performed using:

${\tau }_{ij}\left(t+1\right)=\left(1-\rho \right){\tau }_{ij}\left(t\right)+{\displaystyle \underset{k=1}{\overset{m}{\sum }}}\Delta {\tau }_{ij}^k\left(t\right)$ (6)

${{\tau }^{\prime }}_{ij}\left(t+1\right)=\left(1-\rho \right){{\tau }^{\prime }}_{ij}\left(t\right)+{\displaystyle \underset{k=1}{\overset{m}{\sum }}}\Delta {{\tau }^{\prime }}_{ij}^k\left(t\right)$ (7)

Then, each ant that has generated a solution in the Pareto front in the current cycle is allowed to update the two pheromone structures, $\tau$ and ${\tau }^{\prime }$ , by depositing an amount equal to:

$\Delta {{\tau }^{\prime }}_{ij}^k\left(t\right)=\frac{1}{l}$ (8)

where $l$ is the number of ants that are updating the pheromone trails. The non-dominated solutions generated throughout the execution of the algorithm are put into an external set.

6. Adjacent Pairwise Exchange Algorithm (APE)

APE is a local search heuristic, which allows obtaining achievable results in a reasonable period of time; it is used to bring an improvement in the quality of a given solution. The APE algorithm starts with a given admissible solution, thereafter, it will search in the vicinity of the current solution for any solution to improve the current configuration. At each step of the improvement procedure, the algorithm tests whether the exchange of two adjacent sums (cities) produces a feasible solution to improve the solution and the algorithm continues in this way until we can no longer improve. APE helps the BicriterionAnt algorithm to explore the search area more strongly in order to find other solutions.

7. Algorithm k-Means

Consider a total of n customers. The demand of each customer is known in advance. The n customers are geographically dispersed over a space, and whose position of a customer is represented by its coordinates $\left(x_i\mathrm{,}{y}_i\right)$ with $i=\mathrm{1,}\cdots \mathrm{,}n$ . The coordinate $\left(x_0\mathrm{,}{y}_0\right)$ is gives the position of the central repository. Each customer must be served by one and only one vehicle. Each client belongs to a cluster. Customers who belong to the same cluster are served by the same vehicle, hence the total number of clusters is the total number of vehicles used.

Consider that all n clients are grouped into m clusters and their requests are $d_1\mathrm{,}{d}_2,\cdots \mathrm{,}{d}_n$ .

The number of vehicles used (the number of clusters) can be obtained by:

$m={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\frac{d_i}{Cap}$ (9)

with $Cap$ the load capacity of the vehicle.

Consider that the Euclidean distance $d_{ij}$ between any customers i of a cluster is symmetric, it can be determined by:

$d_{ij}=\sqrt{{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2}$ (10)

where $i=1,\cdots ,n$ and $j=1,\cdots ,n$ .

Let M be the assignment matrix, which represents which client is assigned to which cluster $m_{ij}=1$ if customer i is assigned to cluster j; 0 otherwise.

Thus, the problem of assigning customers to different clusters can be formulated mathematically as

$\min {\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{d}_{ij}{m}_{ij}$ (11)

${\displaystyle \underset{j=1}{\overset{m}{\sum }}}{m}_{ij}=1,i=1,2,\cdots ,n$ (12)

${\displaystyle \underset{i=1}{\overset{n}{\sum }}}{d}_{ij}{m}_{ij}\le Cap,j=1,2,\cdots ,m$ (13)

It is therefore a minimization problem, which seeks to minimize the distance. Constraint (12) represents that each customer is assigned to a single cluster, and constraint (13) is the vehicle or cluster capacity constraint.

The above problem of allocating customers to different clusters can be solved by any clustering method like k-means algorithm. This algorithm finds customer clusters based on customer locations and Euclidean distances to centroids. However, the classical k-means algorithm operates on randomness. The modified k-means algorithm works well in the context of VRP by imposing some specific constraints such as:

● Instead of an arbitrary choice of the value of k, it is appropriate to choose the value of k as the number of vehicles using Equation (9).

● It is better to start by choosing the centroid of the cluster as the depot so that all clients always come to the position closest to the depot.

● It is better to focus first on customers with higher demands than others.

● Instead of choosing customers randomly for different clusters initially, it is better to calculate the Euclidean distance of each centroid and then assign it to the nearest cluster.

Figure 3 shows the structure of our hybrid algorithm BicriterionAntAPE.

8. The Main Formulas for Calculating CO₂ Emissions

The main formulas for calculating CO₂ emissions The Carbon Footprint makes it possible to estimate the impact of an activity on the environment. It was developed in 2004 by the ecological transition agency ADEME, which established that the emission factor for road freight transport is evaluated at 3.07 (kg CO₂/tonnes.km). CO₂ emission factors indicate the quantity of CO₂ emitted during the combustion of a given fuel and for a unit of energy taken.

Four ways to calculate the amount of CO₂ emitted during a tour or service

Figure 3. Flowchart of the hybrid BicriterionAntAPE algorithm.

according to different parameters such as energy consumption rate (TC), emission factor (EF), distance (d), vehicle capacity (CV) and the number of units transported in the vehicle (UT). One of them is the following: The energy source consumption is not known for the service in particular, and the means of transport concerns several beneficiaries. On the one hand, it is necessary to estimate consumption using average consumption and the journey, and on the other hand, to distribute the emissions among the beneficiaries. The calculation formula is

$E_{{\text{CO}}_{\text{2}}}=TC\times d\times EF\left(\frac{UT}{CV}\right)$ (14)

For this purpose, we assumed two fuel consumption rates. The first is for the vehicle of the type k when unloaded ${\mu }_k^0$ , and the second is for the same vehicle when loaded to its maximum capacity ${\mu }_k^{\mathrm{<mo>\; *\; </mo>}}$ then, taking into account the type, load and distance traveled by the vehicles, the CO₂ emissions per unit distance will be determined based on the equation:

$E_{{\text{CO}}_{\text{2}}}=\chi \left({\mu }_k^0+\left(\frac{{\mu }_k^{*}-{\mu }_k^0}{R_k}\right){\beta }_k\right)$ (15)

where $\chi$ : the CO₂ emission rate, $R_k$ : is the maximum load capacity of a vehicle the type k and ${\beta }_k$ : is the vehicle load.

9. Description of the Problem

Our problem consists of two objective functions. The first objective seeks to minimize the overall distance of the transport and the second is that of minimizing the CO₂ emissions of each vehicle during a round. In fact, the most common definition of VRP is the design of optimal routes by a fleet of homogeneous or heterogeneous vehicles located in one or more depots to serve a number of geographically dispersed customers with known demand. This general definition highlights a set of parameters that characterize VRP variants: Transport networks, customers and vehicle fleets. Additional constraints can be added to these three parameters.

10. Formulation of the Selected Problem

Following the brief description of the problem given above, we formulate the problem to be solved as follows:

$f_1=\min {\displaystyle \underset{k=1}{\overset{v}{\sum }}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\displaystyle \underset{j=1}{\overset{n}{\sum }}}{d}_{ij}{x}_{ij}^k$ (16)

$f_2=\min {\displaystyle \underset{k\in K}{\sum }}{\displaystyle \underset{i\in V}{\sum }}{\displaystyle \underset{j\in V}{\sum }}\chi \left({\mu }_k^0{x}_{ij}^k+\left(\frac{{\mu }_k^{*}-{\mu }_k^0}{R_k}\right)q_{ij}^k\right)d_{ij}$ (17)

${\displaystyle \underset{k=1}{\overset{v}{\sum }}}{\displaystyle \underset{j=0}{\overset{n}{\sum }}}{x}_{ij}^k=1i\in V-\left\{0\right\}$ (18)

${\displaystyle \underset{i=0}{\overset{n}{\sum }}}{x}_{ij}^k={\displaystyle \underset{j=0}{\overset{n}{\sum }}}{x}_{ji}^{k}l\in V-\left\{0\right\}k=1,\cdots ,v$ (19)

${\displaystyle \underset{j=1}{\overset{n}{\sum }}}{x}_{0j}^k=1k=1,\cdots ,v$ (20)

${\displaystyle \underset{i=1}{\overset{n}{\sum }}}{x}_{i0}^k=1k=1,\cdots ,v$ (21)

${\displaystyle \underset{i=0}{\overset{n}{\sum }}}{\displaystyle \underset{j=0}{\overset{n}{\sum }}}{s}_i{x}_{ij}^k+{\displaystyle \underset{i=0}{\overset{n}{\sum }}}{\displaystyle \underset{j=0}{\overset{n}{\sum }}}{t}_{ij}{x}_{ij}^k\le Tk=1,\cdots ,v$ (22)

${\displaystyle \underset{i\in U}{\sum }}{\displaystyle \underset{j\in U}{\sum }}{x}_{ij}^k\ge {\displaystyle \underset{j=1}{\overset{n}{\sum }}}{x}_{ij}^k\forall U\subset V-\left\{0\right\};l\in U;k=1,\cdots ,v$ (23)

${\displaystyle \underset{i\in V}{\sum }}{\displaystyle \underset{j\in V_{i\ne j}}{\sum }}D{e}_i{x}_{ij}^k\le R_k{Z}^{k}k=1,\cdots ,v$ (24)

$x_{ij}^k{R}_k\ge q_{ij}^k\forall i\in V-\left\{0\right\}\forall j\in V,k=1,\cdots ,v$ (25)

$x_{ij}^k\mathrm{,}{Z}^k\in \left\{\mathrm{0,1}\right\}\forall i\mathrm{,}j\mathrm{,}k$ (26)

$q_{ij}^k\ge 0\forall i\mathrm{,}j\mathrm{,}k$ (27)

(18) Each client $\forall i\in V-\left\{0\right\}$ is visited once and once. (19) If vehicle k arrives at customer j, it must leave. (20) and (21) Each must return to the depot at the end of their tour. (22) compliance with the maximum duration T of the trips. (23) Elimination of sub-tours to ensure route connectivity. (24) Ensures that each vehicle k could carry a load less than or equal to its maximum capacity $\left(R_k\right)$ . (25) Gets $q_{ij}^k$ for a value greater than zero. (26) $x_{ij}^k,z^k=1$ if the vehicle of type k visits j, and 0 otherwise . (27) $y_{ij}^k$ flux of the arc $\left(i\mathrm{,}j\right)$ loaded in the vehicle of type k ensure that the flux variables are non-negative.

To arrive at routes that produce the minimum of CO₂ emissions, it is essential to build a matrix of CO₂ emissions based on the estimate of CO₂ emitted during the course of each ridge. In this paper, only diesel engines complying with the most recent standards in force are taken into account. For good service performance, the constraint linked to the maximum duration is set at 5 hours for vehicle rounds with an average speed of 55 km/h for urban areas.

11. Analysis of Results

To test the proposed algorithm on the instances in the literature, the Matlab R2021a software allowed us to do the coding, installed on a laptop with an Intel (R) Celeron (R) with a CPU@1.60 GHz processor and 8 GB of RAM, with a 64-bit operating system.

We applied our hybrid BicriterionAnt algorithm to solve some VRP instances from the literature such as Golden_4, Golden_8, Golden_7 proposed by Golden, Wasil, Kelly, and Chao (1998) and X-n167-k10, X-n106-k14, X-n143-k7 proposed by Uchoa, Pecin, Pessoa, Poggi, Subramanian, and Vidal (2013). For NSGAII and SPEA2 the parameters used are: Iterations = 200, Population = 100, Crossover = 0.7, mutation = 0.3.

Case Study

Here we apply our BCAntAPE hybrid method to the problem of a dairy products company has a central depot and a range of uniform delivery vehicles with a maximum capacity of ten tonnes. The demands of its 15 customers are known; we consider the distances between symmetric customers and verify the triangular inequality, as shown in Table 1.

Table 1. Matrix of distance between customers.

Table 2 gives the characteristics of the heavy duty vehicle (HDV) [38] used are:

We tested the algorithms by varying the number of ants from 4, 3, and 2. We considered the following parameters as shown in Table 3.

From the results obtained in Table 4 we see that a minimum CO₂ pollution is not guaranteed for vehicles taking the shortest route. However, several common routes exist (i.e. the visit sequences are the same with different optimization goals).

Table 5 shows the results obtained by NSGAII, SPEA2, BCAnt and BCAntAPE for Golden_4, Golden_8, Golden_7, X-n167-k10, X-n106-k14 instances and X-n143-k7 relative to OF1 and OF2.

12. Performance Measures

Four metrics were used to evaluate the performance of the algorithms, in order to perform a comparative analysis between them. The metrics called M₁, M₂ and M₃, taken from Zitzler et al. [39] , refer respectively to the evaluation of the quality, the distribution and the extension of the Pareto front generated by the algorithm. The fourth metric called Error was taken from Veldhuizen [40] and refers to the percentage of generated solutions that do not belong to the Pareto Front.

In all cases, the Euclidean distance between two points is used as the distance metric, represented by $d\left(p\mathrm{,}q\right)$ , and F is taken as the known Pareto front and

Table 2. The characteristics of HDV.

Table 3. Parameter description table.

Table 4. The statitics of differents algorithms.

Table 5. The summary results of all the instances using BCAnt, BCAntAPE, SPEA2 and NSGAII.

$F^{\prime }$ as the Pareto front generated by the algorithms. The four measures used are defined below.

$M_1\left(F\right)=\frac{1}{\left|F\right|}{\displaystyle \underset{p\in F^{\prime }}{\sum }}\min \left\{d\left(p\mathrm{,}\bar{p}\right)\mathrm{<mo>\; :\; </mo>}\bar{p}\in F\right\}$ (28)

The metric M₂ uses a set $W=\left\{q\in F^{\prime }:d\left(p,q\right)>\delta \right\}$ , given a problem dependent distance value $\delta$ and is calculated as follows:

$M_2\left(F\right)=\frac{1}{\left|F\right|-1}{\displaystyle \underset{p\in F^{\prime }}{\sum }}\left|W\right|$ (29)

The M₃ metric calculates the extension of the generated Pareto front according to:

$M_3\left(F^{\prime }\right)=\sqrt{{\displaystyle \underset{i=1}{\overset{b}{\sum }}}\left\{\max d\left(p_i,q_i\right):p,q\in F^{\prime }\right\}}$ (30)

The error metric is calculated as follows:

$E=\frac{{\displaystyle {\sum }_{i=1}^{F^{\prime }}{n}_i}}{\left|F^{\prime }\right|}$ (31)

where $n_i$ takes the value 0 if the ith solution of $F^{\prime }$ belongs to F, 1 otherwise.

The metrics M₁, M₂ and M₃ have been normalized with respect to their maximum value; in this way the results obtained represent percentages which are used to compare the various algorithms in the set of test problems.

For the metric M₂, the value of $\delta$ was used as ten percent of the distance between the point with the best evaluation in the first objective and the point with the best evaluation in the second objective of the front F, so $\delta$ represents ten percent of the distance between the extreme points of F for each particular problem.

In the used instances of the VRP problems, the BCAntAPE and NSGAII obtain better results than SPEA2 and BCAnt. As shown in Table 6, BCAntAPE for M₂ has the second best average relative distance to the Pareto edge, metric M₃, it gives the best result regarding the average of the Pareto edges closer to the optimal than the other algorithms. We see that BCAntAPE produced the best result for the M₁ metric. Finally, he obtained the second best result in relation to the metric E.

From what precedes the algorithm proposed in this paper, the BCAntAPE has been shown to be as competitive as the other algorithms retained such as NSGAII, SPEA2 and BCAnt.

Table 6. Results of the different metrics.

13. Practical Perspectives

Through the above analyzes, several findings from this research have direct managerial insights into the delivery problem. The results show that the routes which minimize the total distance traveled by the vehicles are different from those which minimize the CO₂ pollution; this can be understood by the fact that the objectives are conflicting.

In this study we also found that the optimal route reduces product CO₂ by almost 7.2% compared to the worst route. Thus, taking road traffic congestion factors into account is an important way to reduce logistics costs and CO₂ emissions.

14. Conclusions

This article studied the multiobjective problem of routing green vehicles. The objective of the proposed model is not only to minimize the total distance traveled by the vehicles of the fleet but also to take into account the carbon emissions of the latter while satisfying the demands of each customer, in addition to that we assumed two fuel consumption rates. The first is for the vehicle of type k when unloaded, and the second is for the same vehicle when loaded to its maximum capacity.

Depending on the characteristics of the selected model, a hybrid metaheuristic has been proposed. This hybrid metaheuristic is based on an extension of the ant colony algorithm BicriterionAnt and an APE improvement heuristic has been associated in order to solve the selected problem.

A case study on a fuel delivery problem was processed to find vehicle routes that pollute less CO₂ and minimize the total distance travelled.

The application of the new BicriterionAntAPE hybrid algorithm is strongly limited to the two-objective problem only and using a single colony for both pheromone structures limits intense exploration of the search area.

Future research might focus on uncertainties for different parameters like travel time, service time, operating costs, etc.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Ressources naturelles Canada. Guide de consommation de carburant 2016. https://www.mcan.gc.ca/energie/efficacite/transports/voirures-camions-legers/achats/7488
[2]	Sheu, J.-B., Chou, Y.-H. and Hu, C.-C. (2005) An Integrated Logistics Operational Model for Green Supply Chain Management. Transportation Research Part E: Logistics and Transportation Review, 41, 287-313. https://doi.org/10.1016/j.tre.2004.07.001
[3]	Moutaoukil, A., et al. (2015) Minimisation des émissions de CO2 dans les circuits de distribution avec flotte hétérogène.
[4]	Rezaei, N., et al. (2019) A Green Vehicle Routing Problem with Time Windows Considering the Heterogeneous Fleet of Vehicles: Two Metaheuristic Algorithms. European Journal of Industrial Engineering, 13, 507-535.
[5]	Ferreira, J.C. and Arns Steiner, M.T. (2021) A Bi-Objective Green Vehicle Routing Problem: A New Hybrid Optimization Algorithm Applied to a Newspaper Distribution. Journal of Geographic Information System, No. 4, 410-433. https://doi.org/10.4236/jgis.2021.134023
[6]	Dutta, J., et al. (2021) A Hybrid Multi-Objective Evolutionary Algorithm for Open Green Vehicle Routing Problem through Cluster Primary-Route Secondary Approach. International Journal of Management Science and Engineering Management, 17, 132-146. https://doi.org/10.1080/17509653.2021.2000901
[7]	Zhalechian, M., Tavakkoli-Moghaddam, R., Zahiri, B. and Mohammadi, M. (2016) Sustainable Design of a Closed-Loop Location-Routing-Inventory Supply Chain Network under Mixed Uncertainty. Transportation Research Part E: Logistics and Transportation Review, 89, 182-214. https://doi.org/10.1016/j.tre.2016.02.011
[8]	Solomon, M.M. (1987) Algorithms for the Vehicle Routing and Scheduling Problems with Time Window Constraints. Operations Research, 35, 254-265. https://doi.org/10.1287/opre.35.2.254
[9]	Bodin, L. and Golden, B. (1981) Classification in Vehicle Routing and Scheduling. Networks, 11, 97-108. https://doi.org/10.1002/net.3230110204
[10]	Baker, B.M. and Ayechew (2003) A Genetic Algorithm for the Vehicle Routing Problem. Computers OR, 30, 787-80. https://doi.org/10.1016/S0305-0548(02)00051-5
[11]	Dorigo, M. and Stuetzle, T. (2009) Ant Colony Optimization: Overview and Recent Advances. IRIDIA Technical Report Series, Technical Report, No. TR/IRIDIA/2009-013, 34 p.
[12]	Dorigo, M. (1992) Optimization, Learning and Natural Algorithms. Ph.D. Thesis, Dipartimento di Elettronica, Politecnico di Milano. (In Italian)
[13]	Dorigo, M., Maniezzo, V. and Colorni, A. (1996) The Ant System: Optimization by a Colony of Cooperating Agents. IEEE Transactions on Systems, Man, and Cybernetics-Part B, 26, 29-41. https://doi.org/10.1109/3477.484436
[14]	Dorigo, M. and Gambardella, L.M. (1997) Ant Colony System: A Cooperative Learning Approach to the Traveling Salesman Problem. IEEE Transactions on Evolutionary Computation, 1, 53-66. https://doi.org/10.1109/4235.585892
[15]	Iredi, S., Merkle, D. and Middendorf, M. (2001) Bi-Criterion Optimization with Multi Colony Ant Algorithms. 1st International Conference on Evolutionary Multi-Criterion Optimization (EMO’01), Vol. 1993, 359-372. https://doi.org/10.1007/3-540-44719-9_25
[16]	Zhang, S., Lee, C.K.M., Choy, K.L., Ho, W. and Ip, W.H. (2014) Design and Development of a Hybrid Artificial Bee Colony Algorithm for the Environmental Vehicle Routing Problem. Transportation Research Part D: Transport and Environment, 31, 85-99. https://doi.org/10.1016/j.trd.2014.05.015
[17]	Dantzig, G.B. and Ramser, J.H. (1959) The Truck Dispatching Problem. Management Science, 6, 80-91. https://doi.org/10.1287/mnsc.6.1.80
[18]	Reimann, M. (2007) Analysing Risk Orientation in a Stochastic VRP. European Journal of Industrial Engineering, 1, 111-130. https://doi.org/10.1504/EJIE.2007.014105
[19]	Gansterer, M. and Hartl, R.F. (2018) Collaborative Vehicle Routing: A Survey. European Journal of Operational Research, 268, 1-12. https://doi.org/10.1016/j.ejor.2017.10.023
[20]	Moons, S., Ramaekers, K., Caris, A. and Arda, Y. (2017) Integrating Production Scheduling and Vehicle Routing Decisions at the Operational Decision Level: A Review and Discussion. Computers Industrial Engineering, 104, 224-245. https://doi.org/10.1016/j.cie.2016.12.010
[21]	Demir, E., Bektaş, T. and Laporte, G. (2014) A Review of Recent Research on Green Road Freight Transportation. European Journal of Operational Research, 237, 775-793. https://doi.org/10.1016/j.ejor.2013.12.033
[22]	Sariklis, D. and Powell, S. (2000) A Heuristic Method for the Open Vehicle Routing Problem. The Journal of the Operational Research Society, 51, 564-573. https://doi.org/10.1057/palgrave.jors.2600924
[23]	Dutta, J., Barma, P., Kar, S. and De, T. (2019) A Modified Kruskal’s Algorithm to Improve Genetic Search for Open Vehicle Routing Problem. International Journal of Business Analytics, 6, 55-57. https://doi.org/10.4018/IJBAN.2019010104
[24]	Gutierrez, A., Dieulle, L., Labadie, N. and Velasco, N. (2018) A Multi-Population Algorithm to Solve the VRP with Stochastic Service and Travel Times. Computers Industrial Engineering, 125, 144-156. https://doi.org/10.1016/j.cie.2018.07.042
[25]	Eshtehadi, R., Fathian, M. and Demir, E. (2017) Robust Solutions to the Pollution-Routing Problem with Demand and Travel Time Uncertainty. Transportation Research Part D: Transport and Environment, 51, 351-363. https://doi.org/10.1016/j.trd.2017.01.003
[26]	Kopfer, H. (2013) Emissions Minimization Vehicle Routing Problem in Dependence of Different Vehicle Classes. Springer, Berlin. https://doi.org/10.1007/978-3-642-35966-8_4
[27]	Yang, B., Hu, Z.-H., Wei, C., Li, S.-Q., Zhao, L. and Jia, S. (2015) Routing with Time-Windows for Multiple Environmental Vehicle Types. Computers Industrial Engineering, 89, 150-161. https://doi.org/10.1016/j.cie.2015.02.001
[28]	Masmoudi, M.A., Hosny, M., Demir, E. and Cheikhrouhou, N. (2018) A Study on the Heterogeneous Fleet of Alternative Fuel Vehicles: Reducing CO2 Emissions by Means of Biodiesel Fuel. Transportation Research Part D: Transport and Environment, 63, 137-155. https://doi.org/10.1016/j.trd.2018.04.025
[29]	Durillo, J.-J., Nebro, A-J., Luna, F., Dorronsoro and Alba, B.E. (2011) jMetal: A Java Framework for Multi-Objective Optimization. Advances in Engineering Software, 42, 760-771. https://doi.org/10.1016/j.advengsoft.2011.05.014
[30]	Deb, K. and Gupta, H. (2005) Searching for Robust Pareto-Optimal Solutions in Multi-Objective Optimization. International Conference on Evolutionary Multi-Criterion Optimization, Vol. 3410, 150-164. https://doi.org/10.1007/978-3-540-31880-4_11
[31]	Marler, R.-T. and Arora, J.-S. (2004) Survey of Multi-Objective Optimization Methods for Engineering. Structural and Multidisciplinary Optimization, 26, 369-395. https://doi.org/10.1007/s00158-003-0368-6
[32]	Dorigo, M., Birattari, M. and Stützle, T. (2006) Ant Colony Optimization: Artificial Ants as a Computational Intelligence Technique. IEEE Computational Intelligence Magazine, 1, 28-39. https://doi.org/10.1109/CI-M.2006.248054
[33]	Garcia-Martinez, C., Cordón, O. and Herrera, F. (2007) A Taxonomy and an Empirical Analysis of Multiple Objective Ant Colony Optimization Algorithms for the Bi-Criteria TSP. European Journal of Operational Research, 180, 116-148. https://doi.org/10.1016/j.ejor.2006.03.041
[34]	Angus, D. and Woodward, C. (2009) Multiple Objective Ant Colony Optimisation. Swarm Intelligence, 3, 69-85. https://doi.org/10.1007/s11721-008-0022-4
[35]	Romero, C.E.M. and Manzanares, E.M. (1999) MOAQ an Ant-Q Algorithm for Multiple Objective Optimization Problems. Genetic and Evolutionary Computing Conference (GECCO), Vol. 1, 894-901.
[36]	Guntsch, M. and Middendorf, M. (2003) Solving Multi-Criteria Optimization Problems with Population-Based ACO. Proceedings of the 2nd International Conference on Evolutionary Multi-Criterion Optimization (EMO), Vol. 2632, 464-478. https://doi.org/10.1007/3-540-36970-8_33
[37]	Lopez-Ibanez, M., Paquete, L. and Stützle, T. (2004) On the Design of ACO for the Biobjective Quadratic Assignment Problem. 4th International Workshop on Ant Algorithms and Swarm Intelligence, Vol. 3172, 214-225. https://doi.org/10.1007/978-3-540-28646-2_19
[38]	Jancovici, M. (2007) Temis-Bilan carbone. Guide des facteurs d’émissions-Calcul des facteurs d’émissions et sources bibliographiques utilisées. Monographie. ADEME France, Paris.
[39]	Zitzler, E., Deb, K. and Thiele, L. (2000) Comparison of Multiobjective Evolutionary Algorithms. Empirical Result. Evolutionary Computation, 8, 173-195. https://doi.org/10.1162/106365600568202
[40]	Van Veldhuizen, D.A. (1999) Multiobjective Evolutionary Algorithms: Classifications, Analyses, and New Innovations. PhD Thesis, Department of Electrical and Computer Engineering, Graduate School of Engineering, Air Force Institute of Technology, Wright-Patterson AFB.