# A New Hybrid Butterfly Optimization Algorithm for Green Vehicle Routing Problem

Authors: Dana Marsetiya Utama; Dian Setiya Widodo; Muhammad Faisal Ibrahim; Shanty Kusuma Dewi
Journal: Journal of Advanced Transportation
 (2020)

---

## Abstract

In the industrial sector, transportation plays an essential role in distribution. This activity impacts climate change and global warming. One of the critical problems in distribution is the green vehicle routing problem (G-VRP). This study focuses on G-VRP for a single distribution center. The objective function is to minimize the distribution costs by considering fuel costs, carbon costs, and vehicle use costs. This research aims to develop the hybrid butterfly optimization algorithm (HBOA) to minimize the distribution costs on G-VRP. It was inspired by the butterfly optimization algorithm (BOA), which was by combining the tabu search (TS) algorithm and local search swap and flip strategies. BOA is a new metaheuristic algorithm that has been successfully applied in various engineering fields. Experiments were carried out to test the parameters of the proposed algorithm and vary the speed of vehicles. The proposed algorithm was also compared with several procedures of prior study. The experimental results proved that the HBOA could minimize the total distribution cost compared to other algorithms. Moreover, the computation time is also included in the analysis.


---

## Body

## 1. Introduction

In distribution, transportation plays an essential role in the industrial sector. It proposes to distribute products from warehouses to customers [1]. Transport and logistics are also significant in the economic development of the world [2]. However, these activities impact climate change and global warming [3, 4]. In China, 30% of the total carbon emissions are caused by the goods transportation sector [5]. Furthermore, in the United States, 28.5% of greenhouse gas is caused by the transportation of goods [6]. Global climate change prevention was declared at the Copenhagen global climate summit in 2009 [7]. Generally, vehicles use fossil fuels as the source of engine-driving energy. Therefore, air pollution is mostly caused by the transportation sector [8, 9]. One of the efforts to solve this problem is determining the right route. The problem of minimizing carbon emissions and fuel energy on this vehicle route is classified as a green vehicle routing problem (G-VRP) [10, 11]. This issue has attracted the attention of many researchers [12, 13]. G-VRP is the development of the classic vehicle routing problem (VRP). The VRP aims to minimize the total cost [14, 15] and the distance of travel [16–18]. However, the G-VRP aims at reducing the environmental impacts, such as reducing fuel consumption and carbon emissions caused by the distribution process [19, 20]. Therefore, several effective procedures have been developed to solve G-VRP. In recent years, there has been an increase in research interest in this problem.The researchers classify G-VRP as NP-hard problem [21–23]. They argue that the search for solutions to these problems is difficult to find at the time of polynomials. Recently, one popular procedure is metaheuristics [24]. Several researchers investigate the issue of G-VRP to reduce fuel consumption and carbon emissions partially. In the fuel consumption problem, some metaheuristic algorithms have been used to solve this problem. Cooray and Rupasinghe [25] proposed genetic algorithms (GA) to reduce energy consumption. Particle swarm optimization (PSO) was proposed by Poonthalir and Nadarajan [26] in a fuel-efficient G-VRP with varying speed. PSO was also used by Norouzi et al. [27] to minimize fuel consumption with time dependency. Zhang et al. [28] offered the ant colony optimization (ACO) algorithm for minimizing fuel consumption in multidepot. Other studies to minimize fuel consumption such as simulated annealing (SA) have been proposed by Kuo [29] and Normasari et al. [30], the revised hybrid intelligent algorithm was developed by Wang et al. [31], and Andelmin and Bartolini [32] offered the heuristic multistart local search procedure. Meanwhile, Macrina et al. [33] proposed a hybrid extensive neighborhood search, and Wang and Lu [10] developed the memetic algorithm with competition.Besides, some researchers have resolved G-VRP to minimize carbon emissions. Some popular algorithms for this problem are GA [34, 35], tabu search (TS) [36–38], the Clarke and Wright algorithm [39], and GA with dynamic programming [40]. The differential evolution algorithm was developed by Kunnapapdeelert and Klinsrisuk [41] to solve G-VRP with pickup and delivery problems. Molina et al. [42] proposed the TS with neighborhood variables to reduce pollutant emissions. On the other hand, several studies on G-VRP to minimize carbon emissions and energy consumption simultaneously have been carried out successfully by researchers. Li et al. [43] proposed a modified PSO to reduce the total costs, including quality loss cost, vehicle operating cost, penalty cost, product freshness cost, emissions cost, and energy cost. TS was offered by Poonthalir and Nadarajan [26] to solve G-VRP by considering heterogeneous fixed fleet. Zhang et al. [7] also proposed the TS algorithm to reduce the total distribution costs, including fuel costs, carbon costs, and vehicle use costs. Hybrid GA was offered by Wang et al. [44] to minimize total cost distribution, which includes carbon emission costs. Shen et al. [45] developed PSO and TS to reduce minimum distribution costs, including penalty costs, the driver salary, fuel costs, and carbon emissions costs. Improved ACO algorithm was proposed by Li et al. [46], and Karagul et al. [47] employed the SA algorithm. Based on the trend of problem-solving methods, advanced metaheuristic algorithms have gained popularity in solving G-VRPs. Even a hybrid algorithm is used to solve this problem, and it has the advantage of solving NP-hard problems. Unfortunately, little research has addressed the use of a hybrid algorithm to solve G-VRP.Recently, one of the advanced algorithms is the butterfly optimization algorithm (BOA). It is a new algorithm that can solve optimization problems proposed by Arora and Singh [48]. BOA has been effectively used to solve problems in various fields. Wen and Cao [49] applied a predicting model for exploring household CO2 emission mitigation strategies. BOA was implemented by Sharma et al. [50] in compression string design, welded beam design, and pressure vessel design. Yıldız et al. [51] used BOA to design automobile suspension components. Several studies have applied the BOA to solve several problems. However, there has not been any research about solving G-VRP using the hybrid butterfly optimization algorithm (HBOA). Those reasons motivate the author to conduct this research. Moreover, although some researchers have investigated G-VRP, minimizing carbon emissions and energy still receives little attention in the research literature. One interesting issue of G-VRP was investigated by Zhang et al. [7]. They solved the G-VRP to minimize the total distribution costs by considering fuel costs, carbon costs, and vehicle use costs with the TS algorithm. Unfortunately, the study of Zhang et al. [7] and previous studies do not consider computation time, and it is an essential aspect of optimization. Therefore, we propose the HBOA to minimize distribution costs that include fuel consumption, carbon emission, and vehicle use costs. There are two main objectives of this research: (1) developing the HBOA to minimize distribution costs of G-VRP and (2) comparing the performance of the proposed algorithm computation time. The HBOA was tested with several experiments to find out the best parameters. It is also compared to several algorithms. This research provides a significant contribution as the HBOA is a new algorithm in the G-VRP.This paper structure is presented as follows: Section2.1 describes assumptions, notations, and problem description; Section 2.2 explains the HBOA algorithm; Section 2.3 presents data and experimental procedures; Section 3 provides results and discussion; and Section 4 presents conclusions and future work.


## 2. Materials and Methods

### 2.1. Assumptions, Notations, and Problem Description

In this section, assumptions, notations, and problem descriptions are based on studies from Zhang et al. [7]. We consider transportation with one distribution center and a set of nodes. Vehicles have equal capacity (homogeneous). The distribution cost considered is fuel consumption cost, carbon emission cost, and a used vehicle cost. Highly total distribution costs require proper distribution management. Therefore, distribution centers need to manage the right transportation routes to minimize the total distribution costs. In this problem, vehicle fuel consumption is based on the distance traveled from node s to node s + 1. We assume that the weight of the additional load M of the vehicle increases fuel consumption p percent. Furthermore, the fuel consumption of unloaded vehicles is also considered in the total distribution costs. Assumptions, notations, and problem descriptions are described in the following section.


#### 2.1.1. Assumptions and Notations

This study employed several assumptions in G-VRP, including the following: (1) the route begins and ends at the distribution center; (2) the costs consist of fuel consumption cost, carbon emission cost, and vehicle rental cost; (3) the vehicle has a fixed load capacity for each trip; (4) fuel, emissions, and vehicle usage costs are fixed; (5) vehicle speed is fixed; (6) the demand for each node is fixed; (7) each customer service time is fixed; and (8) this problem considers one distribution center. This study used notations to make it more practical to decipher the problem description. The notations are as follows:TDC: total distribution costRrs: the sth node on the rth route (for example, R23 = 1, the 2 path is 0-3-1-7-0, and node in the 3 is 1)dRrsRrs+1: distance in rth route from node s to node s + 1FCRrsRrs+1: total fuel consumption in rth route from node s to node s + 1KPLRrsRrs+1: the traveled distance per unit fuel in rth route from node s to node s + 1LPHRrsRrs+1: the fuel consumption per unit time of unload vehicle in rth route from node s to node s + 1vRrsRrs+1: the average speed of unloaded vehiclesLRrsRrs+1: load of vehicle in rth route from sth node to s + 1 nodeM: additional load weightp: percentage increase of fuelN: number of vehicles or number of routesVr: number of nodes on route r, r = 1, 2, …, NstRrs: service time at the sth node on the rth routeqRrs: demand at the sth node on the rth routeQ: vehicle capacityCf: fuel consumption cost (fuel prices)Ce: emission carbon cost per unit of fuel consumptionCV: vehicle usage cost per unit time


#### 2.1.2. Problem Description

This study made a mathematical model to describe the problem. The mathematical model is used to minimize distribution costs. The distribution costs considered are fuel cost, carbon emissions cost, and vehicle usage cost. Furthermore, the problem description illustrated is modeled as follows:(1)Min TDC=∑r=1N∑s=1vr−1Cf+Ce×LPHRrsRrs+1×dRrsRrs+1VRrsRrs+1×1+p×LRrsRrs+1M+∑r=1N∑s=1Vr−1CV×dRrsRrs+1VRrsRrs+1+stRrs,subject to(2)∑s=2Vr−1qRrs≤Q,∀r=1,2,…,N,(3)LRrsRrs+1=∑s′=s+1Vr−1qRrs′,∀r=1,2,…,N,(4)Rr1=RrVr=0,∀r=1,2,…,N,(5)N≥0,Vr≥0,Rrs⊂V,∀r=1,2,…,N,∀s=1,2,…,Vr.Equation (1) formulates the objective function in minimizing the total distribution cost, including vehicle use cost, fuel consumption cost, and carbon emissions cost. The cost of fuel consumption and carbon emissions cost considers the increase in fuel consumption (p) for each additional load (M). The fuel consumption per unit time of the unloading vehicle in rth route from node s to node s + 1 also is considered. Constraints (2) and (3) describe formulas to ensure that the total load does not exceed the vehicle capacity. On each route r, the total vehicle load must not exceed the vehicle capacity. The total vehicle load must be ensured that it does not exceed the capacity. It becomes essential in the G-VRP. Constraint (4) shows that the first and last nodes of each vehicle route are the distribution center. As a G-VRP with one distribution center, this constraint ensures that each route starts at the distribution center and also ends at the distribution center. Constraint (5) formulates well-defined decision variables. This constraint defines the number of nodes and routes ≥0, and it describes the decision variable at the sth node on the rth route.


### 2.2. The Proposed Hybrid Butterfly Optimization Algorithm (HBOA)

This study offered HBOA to solve G-VRP. The proposed algorithm was inspired by a BOA metaheuristic algorithm by combining the TS heuristic algorithm and the local search strategy. The main inspiration for the proposed algorithm was from BOA. The BOA was initially proposed by Arora and Singh [48] in 2019. There are two main characteristics in BOA, namely, the fragrance and movement of butterflies. These characteristics distinguish BOA from other algorithms. The basic BOA is shown in Figure 1.Figure 1
Pseudocode of butterfly optimization algorithm (BOA).Although the BOA has been proposed, this algorithm has not yet been satisfactory as it can only solve continuous problems. Meanwhile, the proposed algorithm is expected to solve G-VRP that constitutes sophisticated and discrete characteristics. G-VRP is categorized as an NP-hard combinatorial problem that must be addressed by a discrete search space. Therefore, this study offered a new approach to solve G-VRP. This research proposed five main steps on HBOA, such as (1) convert search agent position to travel order with large rank value (LRV), (2) change the position of 10% search agent based on the TS algorithm, (3) fragrance update, (4) movement of butterflies, and (5) local search. This study proposed an LRV procedure for converting continuous values to discrete values. To improve the BOA performance, this study combined TS and local search algorithms. Swap and flip rules were suggested in the local search strategy. The proposed algorithm is shown in Figure2. The five stages of the proposed algorithm are described in the following section.Figure 2
Pseudocode of hybrid butterfly optimization algorithm (HBOA).


#### 2.2.1. Convert Search Agent Position to Travel Order with Large Rank Value (LRV)

In this section, initializing the search agent position was generated randomly according to the upper and lower limits. The upper and lower limit values were set to determine the position of the BOA agent. At this stage, the search agent’s position was ensured with no repeating numbers on the same search agent (Figure3).Figure 3
Initialization of the search agent position.Furthermore, we proposed the principle of LRV to convert from the position of the search agent (continuous value) to the order of travel (discrete value). LRV is a popular method to convert from continuous value to discrete value in combinatorial problems [52–55]. At this stage, each search agent position value was sorted from the largest value to the smallest one. The LRV representation is shown in Figure 4. However, the illustration of Figure 4(b) could not be applied because, in one position, the search agent had the same value (0.43). In other words, the search agent position could not be applied as the order/route of vehicles visited the same two places.Figure 4
LRV representation. (a) Correct travel sequence and (b) wrong travel sequence.
(a)(b)


#### 2.2.2. Tabu Search Algorithm

In this section, this study proposed that 10% of the initial search agent positions were adjusted to the tabu search (TS) algorithm solution. The TS algorithm is a popular heuristic algorithm widely used to solve G-VRP. This study used the TS procedure developed by Poonthalir and Nadarajan [26]. To solve G-VRP, the five main stages [37] of the TS algorithm comprise (1) representation of solution, (2) initial solution, (3) neighborhood solution, (4) tabu list, and (5) criteria for aspiration and dismissal. Figure 5 represents the stages of the TS algorithm. The TS algorithm used three neighborhood solution rules. These rules were swap (Figure 6), flip (Figure 7), and slide (Figure 8). Swap is a rule in which it swaps two nodes. Flip is the rule of a node exchange by reversing the order of the node.Figure 5
Flowchart of the tabu search algorithm.Figure 6
Swap illustration.Figure 7
Flip illustration.Figure 8
Slide illustration.Meanwhile, a slide is an exchange of nodes by shifting their sequences. For thetth iteration to t, the swap and flip rules were repeated nn−1/2 times. The slide rule was repeated n2 times in each tth iteration. For the solution inspection stage, the TS algorithm checked the tabu test by using the tabu list. It was to avoid repetition in finding a solution. In the aspiration criteria stage, the TS algorithm compared the new solution in the iteration t to the previous solution in the iteration t−1. The new solution would be listed as the best solution if it had a better quality than that of the previous one. Furthermore, the stopping criteria used in the TS algorithm referred to the number of fulfilled iteration.As mentioned earlier, this study proposed that 10% of search agent positions were adjusted to the TS algorithm’s solution. Therefore, the position of the new search agent had to be adjusted to the TS solution. This study also suggested a new position adjustment procedure. The illustration of converting TS solutions to the position of the search agent is exemplified in Figure9.Figure 9
Illustration of tabu search converted solution to the position of search agent.


#### 2.2.3. Update Fragrance

In BOA, each butterfly has a unique fragrance and personality. It is one of the main characteristics that distinguish BOA from any other metaheuristic algorithms. All BOA behaviors are based on the sensory modality (c), stimulus intensity (I), and exponential strength (a). Fragrance (f) is formulated as a function of the physical intensity of the stimulus from BOA. The f formula is presented as follows:(6)f=c.Ia,where f is the value of fragrance that changes in every iteration. This value shows how strong the fragrance is felt by other butterflies. The butterfly stimulus intensity is formulated as I. a is the power exponent that depends on the modality. c formulates the sensory modality. Values of a and c on the used butterfly are in the range [0, 1].


#### 2.2.4. Movement of Butterflies

This section explains the phase of the movement of butterflies. There are two main phases in the basic BOA, namely, the initial phase and the movement of the butterfly phase. In the butterfly phase movement, the butterflies move their position as many times as the number of iterations. In this phase, all butterflies in the solution room move to a new position. Then, the fitness value of each butterfly is evaluated. In each iteration, the fitness value of all butterflies is updated. Furthermore, the butterflies produce fragrance in their calculated position based on equation (6). Two movements in BOA are the global search phase and the local search phase. In the global search phase, butterflies take steps towards other butterflies that have the best solution. The global search phase for butterflies is represented in equation (7). Meanwhile, the local search phase is shown in equation (8):(7)Xit+1=Xit+r2g∗−Xitfi,where Xit is the vector solution Xi for the ith butterfly in the iteration t. r shows a random number in the range [0, 1]. g∗ indicates the best solution in the current iteration. The ith butterfly fragrance is represented by fi.(8)Xkt+1=Xit+r2Xjt−Xktfi.Equation (8) indicates the local butterfly search formula. Xjt and Xkt are the j-th and k-th butterflies from the solution room. r is the random number in the range [0, 1]. Movement of butterflies stops until the termination criteria are met. The stopping criterion used is the maximum number of the achieved iteration. After the movement of butterflies, the algorithm produces the best solution based on the fitness values.


#### 2.2.5. Local Search

To improve the BOA performance, this study proposed the local search procedure. Swap and flip were the two local search rules chosen to improve the BOA performance. Figure6 illustrates the proposed swap rules. In this rule, two positions (nodes) were chosen randomly and exchanged. Another local search rule used was flip. In this rule, two nodes were selected randomly and continued to reverse the order of the selected nodes. This rule is illustrated in Figure 7. In the proposed HBOA, for each iteration t, the swap and flip operations were repeated as many as the number of nodes.


### 2.3. Data and Experimental Procedure

#### 2.3.1. Data

In this study, the data of the number of nodes, coordinates, vehicle capacity, and demand were taken from Gaskell [56] and Christofides and Eilon [57]. They used cases with nodes as many as 22 nodes (Table 1) [56], 32 nodes (Table 2) [57], and 50 nodes (Table 3) [57]. Distance (dRrsRrs+1) in rth route from node s to node s+1 is based on formula dRrsRrs+1=Xs−Xs+12+Ys−Ys+12. Meanwhile, the data of the costs and speed data were obtained from Zhang et al. [7]. The fuel cost was 7.3 yuan/liter [7]. The carbon emissions were 0.64 yuan/liter [7]. Furthermore, the vehicle usage fee was 80 yuan/hour [7]. This research employed three categories of vehicle speed (high, medium, and low speed). The high, medium, and low speeds were 107 km/hr, 63 km/hr, and 43 km/hr, respectively. Nine variations of problems (three nodes and three-speed variations) were carried out in this study. Service time for each customer is 0.1 hours. The increase in fuel consumption (p) for each additional load M = 50 is 2%.Table 1
Problem of 22 nodes.
NodeCoordinateDemandXYDepot26623501295272125230125884330926060421727450052182783006282267175724224935082302621509249268110010256267410011265257225122672423001325926525014315233500153292521501631825210017329224250182672131201927519260020303201500212082171752232618175Capacity: 4500.Table 2
Problem of 32 nodes.
NodeCoordinateDemandXYDepot29242501298427700230944540033074644004336475120053204394063214378073224372000832343390093244336001032342975011314435150012311442150133044272501429342116001529641845016261384700172974105501831540765019314406200203213914002132139830022314394130023313378700243043827502529540214002628340640002727939960028271401100029264414500302774392500312904341700323194331100Capacity: 8000.Table 3
Problem of 50 nodes.
NodeCoordinateDemandXYDepot30400137527249493035264164202695403021621471571763198316223952331110512151142411912313229135252314124221153616101652411517272331817334119131392057582821624282242578231657162485210257382826276872730481528436714295848630582719313769113238461233461023346133263562631736636963732229384535153959151440567411017274221101343564114430151645391010463239547253225482555174948281850563710Capacity: 80.


#### 2.3.2. Experimental Procedure

The experiments were designed to determine the effect of HBOA (iteration and population) and speed parameters on the distribution cost and computation time. The experiments were carried out with different parameters. The parameters included the number of populations and iterations. The population parameters used three different levels (10, 50, and 100 populations). The iteration parameters also employed three levels (10, 50, and 100 iterations). This study used the sensory modality of 0.01 and power exponent of 0.1 from BOA parameters. Eighty-one trials were designed in this study. Each result of the trial was recorded for cost and computation time.The HBOA was compared to other algorithms such as BOA [48], TS [7], SA [29], ACO [28], PSO [27], and GA [58]. To compare with several algorithms, this study used one hundred iteration parameters at each vehicle speed in every algorithm. One hundred populations were used in the BOA experiment. Moreover, we used an initial temperature parameter of 1000 and the cooldown factor based on the Kuo [29] formula. One hundred ant populations were adopted for the ACO algorithm. One hundred particles and an inertia weight of 0.5 are used in the PSO algorithm. 100 populations, a crossover probability of 0.8, and a mutation probability of 0.25 were applied in the GA algorithm experiment. The performance was measured using relative error percentage (REP) as presented in equation (9). A positive REP showed that the proposed algorithm is better than the other algorithms. However, a negative REP showed that the proposed algorithm is not competitive as compared to other algorithms.(9)REP=costother algorithms−costproposed algorithmcostproposed algorithm×100%.Besides, this study also compared the computation time in all cases. It was carried out to determine the time efficiency of solving G-VRP. The effect of iteration (t) on distribution costs was also analyzed. This analysis was carried out in 50 nodes, 100 populations, and 100 iterations in the case of medium-speed vehicles. Furthermore, all experiments were conducted with the means of Matlab R2014a software on Windows 8 Intel Celeron with x64-64 2GB RAM processor.


## 2.1. Assumptions, Notations, and Problem Description

In this section, assumptions, notations, and problem descriptions are based on studies from Zhang et al. [7]. We consider transportation with one distribution center and a set of nodes. Vehicles have equal capacity (homogeneous). The distribution cost considered is fuel consumption cost, carbon emission cost, and a used vehicle cost. Highly total distribution costs require proper distribution management. Therefore, distribution centers need to manage the right transportation routes to minimize the total distribution costs. In this problem, vehicle fuel consumption is based on the distance traveled from node s to node s + 1. We assume that the weight of the additional load M of the vehicle increases fuel consumption p percent. Furthermore, the fuel consumption of unloaded vehicles is also considered in the total distribution costs. Assumptions, notations, and problem descriptions are described in the following section.


### 2.1.1. Assumptions and Notations

This study employed several assumptions in G-VRP, including the following: (1) the route begins and ends at the distribution center; (2) the costs consist of fuel consumption cost, carbon emission cost, and vehicle rental cost; (3) the vehicle has a fixed load capacity for each trip; (4) fuel, emissions, and vehicle usage costs are fixed; (5) vehicle speed is fixed; (6) the demand for each node is fixed; (7) each customer service time is fixed; and (8) this problem considers one distribution center. This study used notations to make it more practical to decipher the problem description. The notations are as follows:TDC: total distribution costRrs: the sth node on the rth route (for example, R23 = 1, the 2 path is 0-3-1-7-0, and node in the 3 is 1)dRrsRrs+1: distance in rth route from node s to node s + 1FCRrsRrs+1: total fuel consumption in rth route from node s to node s + 1KPLRrsRrs+1: the traveled distance per unit fuel in rth route from node s to node s + 1LPHRrsRrs+1: the fuel consumption per unit time of unload vehicle in rth route from node s to node s + 1vRrsRrs+1: the average speed of unloaded vehiclesLRrsRrs+1: load of vehicle in rth route from sth node to s + 1 nodeM: additional load weightp: percentage increase of fuelN: number of vehicles or number of routesVr: number of nodes on route r, r = 1, 2, …, NstRrs: service time at the sth node on the rth routeqRrs: demand at the sth node on the rth routeQ: vehicle capacityCf: fuel consumption cost (fuel prices)Ce: emission carbon cost per unit of fuel consumptionCV: vehicle usage cost per unit time


### 2.1.2. Problem Description

This study made a mathematical model to describe the problem. The mathematical model is used to minimize distribution costs. The distribution costs considered are fuel cost, carbon emissions cost, and vehicle usage cost. Furthermore, the problem description illustrated is modeled as follows:(1)Min TDC=∑r=1N∑s=1vr−1Cf+Ce×LPHRrsRrs+1×dRrsRrs+1VRrsRrs+1×1+p×LRrsRrs+1M+∑r=1N∑s=1Vr−1CV×dRrsRrs+1VRrsRrs+1+stRrs,subject to(2)∑s=2Vr−1qRrs≤Q,∀r=1,2,…,N,(3)LRrsRrs+1=∑s′=s+1Vr−1qRrs′,∀r=1,2,…,N,(4)Rr1=RrVr=0,∀r=1,2,…,N,(5)N≥0,Vr≥0,Rrs⊂V,∀r=1,2,…,N,∀s=1,2,…,Vr.Equation (1) formulates the objective function in minimizing the total distribution cost, including vehicle use cost, fuel consumption cost, and carbon emissions cost. The cost of fuel consumption and carbon emissions cost considers the increase in fuel consumption (p) for each additional load (M). The fuel consumption per unit time of the unloading vehicle in rth route from node s to node s + 1 also is considered. Constraints (2) and (3) describe formulas to ensure that the total load does not exceed the vehicle capacity. On each route r, the total vehicle load must not exceed the vehicle capacity. The total vehicle load must be ensured that it does not exceed the capacity. It becomes essential in the G-VRP. Constraint (4) shows that the first and last nodes of each vehicle route are the distribution center. As a G-VRP with one distribution center, this constraint ensures that each route starts at the distribution center and also ends at the distribution center. Constraint (5) formulates well-defined decision variables. This constraint defines the number of nodes and routes ≥0, and it describes the decision variable at the sth node on the rth route.


## 2.1.1. Assumptions and Notations

This study employed several assumptions in G-VRP, including the following: (1) the route begins and ends at the distribution center; (2) the costs consist of fuel consumption cost, carbon emission cost, and vehicle rental cost; (3) the vehicle has a fixed load capacity for each trip; (4) fuel, emissions, and vehicle usage costs are fixed; (5) vehicle speed is fixed; (6) the demand for each node is fixed; (7) each customer service time is fixed; and (8) this problem considers one distribution center. This study used notations to make it more practical to decipher the problem description. The notations are as follows:TDC: total distribution costRrs: the sth node on the rth route (for example, R23 = 1, the 2 path is 0-3-1-7-0, and node in the 3 is 1)dRrsRrs+1: distance in rth route from node s to node s + 1FCRrsRrs+1: total fuel consumption in rth route from node s to node s + 1KPLRrsRrs+1: the traveled distance per unit fuel in rth route from node s to node s + 1LPHRrsRrs+1: the fuel consumption per unit time of unload vehicle in rth route from node s to node s + 1vRrsRrs+1: the average speed of unloaded vehiclesLRrsRrs+1: load of vehicle in rth route from sth node to s + 1 nodeM: additional load weightp: percentage increase of fuelN: number of vehicles or number of routesVr: number of nodes on route r, r = 1, 2, …, NstRrs: service time at the sth node on the rth routeqRrs: demand at the sth node on the rth routeQ: vehicle capacityCf: fuel consumption cost (fuel prices)Ce: emission carbon cost per unit of fuel consumptionCV: vehicle usage cost per unit time


## 2.1.2. Problem Description

This study made a mathematical model to describe the problem. The mathematical model is used to minimize distribution costs. The distribution costs considered are fuel cost, carbon emissions cost, and vehicle usage cost. Furthermore, the problem description illustrated is modeled as follows:(1)Min TDC=∑r=1N∑s=1vr−1Cf+Ce×LPHRrsRrs+1×dRrsRrs+1VRrsRrs+1×1+p×LRrsRrs+1M+∑r=1N∑s=1Vr−1CV×dRrsRrs+1VRrsRrs+1+stRrs,subject to(2)∑s=2Vr−1qRrs≤Q,∀r=1,2,…,N,(3)LRrsRrs+1=∑s′=s+1Vr−1qRrs′,∀r=1,2,…,N,(4)Rr1=RrVr=0,∀r=1,2,…,N,(5)N≥0,Vr≥0,Rrs⊂V,∀r=1,2,…,N,∀s=1,2,…,Vr.Equation (1) formulates the objective function in minimizing the total distribution cost, including vehicle use cost, fuel consumption cost, and carbon emissions cost. The cost of fuel consumption and carbon emissions cost considers the increase in fuel consumption (p) for each additional load (M). The fuel consumption per unit time of the unloading vehicle in rth route from node s to node s + 1 also is considered. Constraints (2) and (3) describe formulas to ensure that the total load does not exceed the vehicle capacity. On each route r, the total vehicle load must not exceed the vehicle capacity. The total vehicle load must be ensured that it does not exceed the capacity. It becomes essential in the G-VRP. Constraint (4) shows that the first and last nodes of each vehicle route are the distribution center. As a G-VRP with one distribution center, this constraint ensures that each route starts at the distribution center and also ends at the distribution center. Constraint (5) formulates well-defined decision variables. This constraint defines the number of nodes and routes ≥0, and it describes the decision variable at the sth node on the rth route.


## 2.2. The Proposed Hybrid Butterfly Optimization Algorithm (HBOA)

This study offered HBOA to solve G-VRP. The proposed algorithm was inspired by a BOA metaheuristic algorithm by combining the TS heuristic algorithm and the local search strategy. The main inspiration for the proposed algorithm was from BOA. The BOA was initially proposed by Arora and Singh [48] in 2019. There are two main characteristics in BOA, namely, the fragrance and movement of butterflies. These characteristics distinguish BOA from other algorithms. The basic BOA is shown in Figure 1.Figure 1
Pseudocode of butterfly optimization algorithm (BOA).Although the BOA has been proposed, this algorithm has not yet been satisfactory as it can only solve continuous problems. Meanwhile, the proposed algorithm is expected to solve G-VRP that constitutes sophisticated and discrete characteristics. G-VRP is categorized as an NP-hard combinatorial problem that must be addressed by a discrete search space. Therefore, this study offered a new approach to solve G-VRP. This research proposed five main steps on HBOA, such as (1) convert search agent position to travel order with large rank value (LRV), (2) change the position of 10% search agent based on the TS algorithm, (3) fragrance update, (4) movement of butterflies, and (5) local search. This study proposed an LRV procedure for converting continuous values to discrete values. To improve the BOA performance, this study combined TS and local search algorithms. Swap and flip rules were suggested in the local search strategy. The proposed algorithm is shown in Figure2. The five stages of the proposed algorithm are described in the following section.Figure 2
Pseudocode of hybrid butterfly optimization algorithm (HBOA).


### 2.2.1. Convert Search Agent Position to Travel Order with Large Rank Value (LRV)

In this section, initializing the search agent position was generated randomly according to the upper and lower limits. The upper and lower limit values were set to determine the position of the BOA agent. At this stage, the search agent’s position was ensured with no repeating numbers on the same search agent (Figure3).Figure 3
Initialization of the search agent position.Furthermore, we proposed the principle of LRV to convert from the position of the search agent (continuous value) to the order of travel (discrete value). LRV is a popular method to convert from continuous value to discrete value in combinatorial problems [52–55]. At this stage, each search agent position value was sorted from the largest value to the smallest one. The LRV representation is shown in Figure 4. However, the illustration of Figure 4(b) could not be applied because, in one position, the search agent had the same value (0.43). In other words, the search agent position could not be applied as the order/route of vehicles visited the same two places.Figure 4
LRV representation. (a) Correct travel sequence and (b) wrong travel sequence.
(a)(b)


### 2.2.2. Tabu Search Algorithm

In this section, this study proposed that 10% of the initial search agent positions were adjusted to the tabu search (TS) algorithm solution. The TS algorithm is a popular heuristic algorithm widely used to solve G-VRP. This study used the TS procedure developed by Poonthalir and Nadarajan [26]. To solve G-VRP, the five main stages [37] of the TS algorithm comprise (1) representation of solution, (2) initial solution, (3) neighborhood solution, (4) tabu list, and (5) criteria for aspiration and dismissal. Figure 5 represents the stages of the TS algorithm. The TS algorithm used three neighborhood solution rules. These rules were swap (Figure 6), flip (Figure 7), and slide (Figure 8). Swap is a rule in which it swaps two nodes. Flip is the rule of a node exchange by reversing the order of the node.Figure 5
Flowchart of the tabu search algorithm.Figure 6
Swap illustration.Figure 7
Flip illustration.Figure 8
Slide illustration.Meanwhile, a slide is an exchange of nodes by shifting their sequences. For thetth iteration to t, the swap and flip rules were repeated nn−1/2 times. The slide rule was repeated n2 times in each tth iteration. For the solution inspection stage, the TS algorithm checked the tabu test by using the tabu list. It was to avoid repetition in finding a solution. In the aspiration criteria stage, the TS algorithm compared the new solution in the iteration t to the previous solution in the iteration t−1. The new solution would be listed as the best solution if it had a better quality than that of the previous one. Furthermore, the stopping criteria used in the TS algorithm referred to the number of fulfilled iteration.As mentioned earlier, this study proposed that 10% of search agent positions were adjusted to the TS algorithm’s solution. Therefore, the position of the new search agent had to be adjusted to the TS solution. This study also suggested a new position adjustment procedure. The illustration of converting TS solutions to the position of the search agent is exemplified in Figure9.Figure 9
Illustration of tabu search converted solution to the position of search agent.


### 2.2.3. Update Fragrance

In BOA, each butterfly has a unique fragrance and personality. It is one of the main characteristics that distinguish BOA from any other metaheuristic algorithms. All BOA behaviors are based on the sensory modality (c), stimulus intensity (I), and exponential strength (a). Fragrance (f) is formulated as a function of the physical intensity of the stimulus from BOA. The f formula is presented as follows:(6)f=c.Ia,where f is the value of fragrance that changes in every iteration. This value shows how strong the fragrance is felt by other butterflies. The butterfly stimulus intensity is formulated as I. a is the power exponent that depends on the modality. c formulates the sensory modality. Values of a and c on the used butterfly are in the range [0, 1].


### 2.2.4. Movement of Butterflies

This section explains the phase of the movement of butterflies. There are two main phases in the basic BOA, namely, the initial phase and the movement of the butterfly phase. In the butterfly phase movement, the butterflies move their position as many times as the number of iterations. In this phase, all butterflies in the solution room move to a new position. Then, the fitness value of each butterfly is evaluated. In each iteration, the fitness value of all butterflies is updated. Furthermore, the butterflies produce fragrance in their calculated position based on equation (6). Two movements in BOA are the global search phase and the local search phase. In the global search phase, butterflies take steps towards other butterflies that have the best solution. The global search phase for butterflies is represented in equation (7). Meanwhile, the local search phase is shown in equation (8):(7)Xit+1=Xit+r2g∗−Xitfi,where Xit is the vector solution Xi for the ith butterfly in the iteration t. r shows a random number in the range [0, 1]. g∗ indicates the best solution in the current iteration. The ith butterfly fragrance is represented by fi.(8)Xkt+1=Xit+r2Xjt−Xktfi.Equation (8) indicates the local butterfly search formula. Xjt and Xkt are the j-th and k-th butterflies from the solution room. r is the random number in the range [0, 1]. Movement of butterflies stops until the termination criteria are met. The stopping criterion used is the maximum number of the achieved iteration. After the movement of butterflies, the algorithm produces the best solution based on the fitness values.


### 2.2.5. Local Search

To improve the BOA performance, this study proposed the local search procedure. Swap and flip were the two local search rules chosen to improve the BOA performance. Figure6 illustrates the proposed swap rules. In this rule, two positions (nodes) were chosen randomly and exchanged. Another local search rule used was flip. In this rule, two nodes were selected randomly and continued to reverse the order of the selected nodes. This rule is illustrated in Figure 7. In the proposed HBOA, for each iteration t, the swap and flip operations were repeated as many as the number of nodes.


## 2.2.1. Convert Search Agent Position to Travel Order with Large Rank Value (LRV)

In this section, initializing the search agent position was generated randomly according to the upper and lower limits. The upper and lower limit values were set to determine the position of the BOA agent. At this stage, the search agent’s position was ensured with no repeating numbers on the same search agent (Figure3).Figure 3
Initialization of the search agent position.Furthermore, we proposed the principle of LRV to convert from the position of the search agent (continuous value) to the order of travel (discrete value). LRV is a popular method to convert from continuous value to discrete value in combinatorial problems [52–55]. At this stage, each search agent position value was sorted from the largest value to the smallest one. The LRV representation is shown in Figure 4. However, the illustration of Figure 4(b) could not be applied because, in one position, the search agent had the same value (0.43). In other words, the search agent position could not be applied as the order/route of vehicles visited the same two places.Figure 4
LRV representation. (a) Correct travel sequence and (b) wrong travel sequence.
(a)(b)


## 2.2.2. Tabu Search Algorithm

In this section, this study proposed that 10% of the initial search agent positions were adjusted to the tabu search (TS) algorithm solution. The TS algorithm is a popular heuristic algorithm widely used to solve G-VRP. This study used the TS procedure developed by Poonthalir and Nadarajan [26]. To solve G-VRP, the five main stages [37] of the TS algorithm comprise (1) representation of solution, (2) initial solution, (3) neighborhood solution, (4) tabu list, and (5) criteria for aspiration and dismissal. Figure 5 represents the stages of the TS algorithm. The TS algorithm used three neighborhood solution rules. These rules were swap (Figure 6), flip (Figure 7), and slide (Figure 8). Swap is a rule in which it swaps two nodes. Flip is the rule of a node exchange by reversing the order of the node.Figure 5
Flowchart of the tabu search algorithm.Figure 6
Swap illustration.Figure 7
Flip illustration.Figure 8
Slide illustration.Meanwhile, a slide is an exchange of nodes by shifting their sequences. For thetth iteration to t, the swap and flip rules were repeated nn−1/2 times. The slide rule was repeated n2 times in each tth iteration. For the solution inspection stage, the TS algorithm checked the tabu test by using the tabu list. It was to avoid repetition in finding a solution. In the aspiration criteria stage, the TS algorithm compared the new solution in the iteration t to the previous solution in the iteration t−1. The new solution would be listed as the best solution if it had a better quality than that of the previous one. Furthermore, the stopping criteria used in the TS algorithm referred to the number of fulfilled iteration.As mentioned earlier, this study proposed that 10% of search agent positions were adjusted to the TS algorithm’s solution. Therefore, the position of the new search agent had to be adjusted to the TS solution. This study also suggested a new position adjustment procedure. The illustration of converting TS solutions to the position of the search agent is exemplified in Figure9.Figure 9
Illustration of tabu search converted solution to the position of search agent.


## 2.2.3. Update Fragrance

In BOA, each butterfly has a unique fragrance and personality. It is one of the main characteristics that distinguish BOA from any other metaheuristic algorithms. All BOA behaviors are based on the sensory modality (c), stimulus intensity (I), and exponential strength (a). Fragrance (f) is formulated as a function of the physical intensity of the stimulus from BOA. The f formula is presented as follows:(6)f=c.Ia,where f is the value of fragrance that changes in every iteration. This value shows how strong the fragrance is felt by other butterflies. The butterfly stimulus intensity is formulated as I. a is the power exponent that depends on the modality. c formulates the sensory modality. Values of a and c on the used butterfly are in the range [0, 1].


## 2.2.4. Movement of Butterflies

This section explains the phase of the movement of butterflies. There are two main phases in the basic BOA, namely, the initial phase and the movement of the butterfly phase. In the butterfly phase movement, the butterflies move their position as many times as the number of iterations. In this phase, all butterflies in the solution room move to a new position. Then, the fitness value of each butterfly is evaluated. In each iteration, the fitness value of all butterflies is updated. Furthermore, the butterflies produce fragrance in their calculated position based on equation (6). Two movements in BOA are the global search phase and the local search phase. In the global search phase, butterflies take steps towards other butterflies that have the best solution. The global search phase for butterflies is represented in equation (7). Meanwhile, the local search phase is shown in equation (8):(7)Xit+1=Xit+r2g∗−Xitfi,where Xit is the vector solution Xi for the ith butterfly in the iteration t. r shows a random number in the range [0, 1]. g∗ indicates the best solution in the current iteration. The ith butterfly fragrance is represented by fi.(8)Xkt+1=Xit+r2Xjt−Xktfi.Equation (8) indicates the local butterfly search formula. Xjt and Xkt are the j-th and k-th butterflies from the solution room. r is the random number in the range [0, 1]. Movement of butterflies stops until the termination criteria are met. The stopping criterion used is the maximum number of the achieved iteration. After the movement of butterflies, the algorithm produces the best solution based on the fitness values.


## 2.2.5. Local Search

To improve the BOA performance, this study proposed the local search procedure. Swap and flip were the two local search rules chosen to improve the BOA performance. Figure6 illustrates the proposed swap rules. In this rule, two positions (nodes) were chosen randomly and exchanged. Another local search rule used was flip. In this rule, two nodes were selected randomly and continued to reverse the order of the selected nodes. This rule is illustrated in Figure 7. In the proposed HBOA, for each iteration t, the swap and flip operations were repeated as many as the number of nodes.


## 2.3. Data and Experimental Procedure

### 2.3.1. Data

In this study, the data of the number of nodes, coordinates, vehicle capacity, and demand were taken from Gaskell [56] and Christofides and Eilon [57]. They used cases with nodes as many as 22 nodes (Table 1) [56], 32 nodes (Table 2) [57], and 50 nodes (Table 3) [57]. Distance (dRrsRrs+1) in rth route from node s to node s+1 is based on formula dRrsRrs+1=Xs−Xs+12+Ys−Ys+12. Meanwhile, the data of the costs and speed data were obtained from Zhang et al. [7]. The fuel cost was 7.3 yuan/liter [7]. The carbon emissions were 0.64 yuan/liter [7]. Furthermore, the vehicle usage fee was 80 yuan/hour [7]. This research employed three categories of vehicle speed (high, medium, and low speed). The high, medium, and low speeds were 107 km/hr, 63 km/hr, and 43 km/hr, respectively. Nine variations of problems (three nodes and three-speed variations) were carried out in this study. Service time for each customer is 0.1 hours. The increase in fuel consumption (p) for each additional load M = 50 is 2%.Table 1
Problem of 22 nodes.
NodeCoordinateDemandXYDepot26623501295272125230125884330926060421727450052182783006282267175724224935082302621509249268110010256267410011265257225122672423001325926525014315233500153292521501631825210017329224250182672131201927519260020303201500212082171752232618175Capacity: 4500.Table 2
Problem of 32 nodes.
NodeCoordinateDemandXYDepot29242501298427700230944540033074644004336475120053204394063214378073224372000832343390093244336001032342975011314435150012311442150133044272501429342116001529641845016261384700172974105501831540765019314406200203213914002132139830022314394130023313378700243043827502529540214002628340640002727939960028271401100029264414500302774392500312904341700323194331100Capacity: 8000.Table 3
Problem of 50 nodes.
NodeCoordinateDemandXYDepot30400137527249493035264164202695403021621471571763198316223952331110512151142411912313229135252314124221153616101652411517272331817334119131392057582821624282242578231657162485210257382826276872730481528436714295848630582719313769113238461233461023346133263562631736636963732229384535153959151440567411017274221101343564114430151645391010463239547253225482555174948281850563710Capacity: 80.


### 2.3.2. Experimental Procedure

The experiments were designed to determine the effect of HBOA (iteration and population) and speed parameters on the distribution cost and computation time. The experiments were carried out with different parameters. The parameters included the number of populations and iterations. The population parameters used three different levels (10, 50, and 100 populations). The iteration parameters also employed three levels (10, 50, and 100 iterations). This study used the sensory modality of 0.01 and power exponent of 0.1 from BOA parameters. Eighty-one trials were designed in this study. Each result of the trial was recorded for cost and computation time.The HBOA was compared to other algorithms such as BOA [48], TS [7], SA [29], ACO [28], PSO [27], and GA [58]. To compare with several algorithms, this study used one hundred iteration parameters at each vehicle speed in every algorithm. One hundred populations were used in the BOA experiment. Moreover, we used an initial temperature parameter of 1000 and the cooldown factor based on the Kuo [29] formula. One hundred ant populations were adopted for the ACO algorithm. One hundred particles and an inertia weight of 0.5 are used in the PSO algorithm. 100 populations, a crossover probability of 0.8, and a mutation probability of 0.25 were applied in the GA algorithm experiment. The performance was measured using relative error percentage (REP) as presented in equation (9). A positive REP showed that the proposed algorithm is better than the other algorithms. However, a negative REP showed that the proposed algorithm is not competitive as compared to other algorithms.(9)REP=costother algorithms−costproposed algorithmcostproposed algorithm×100%.Besides, this study also compared the computation time in all cases. It was carried out to determine the time efficiency of solving G-VRP. The effect of iteration (t) on distribution costs was also analyzed. This analysis was carried out in 50 nodes, 100 populations, and 100 iterations in the case of medium-speed vehicles. Furthermore, all experiments were conducted with the means of Matlab R2014a software on Windows 8 Intel Celeron with x64-64 2GB RAM processor.


## 2.3.1. Data

In this study, the data of the number of nodes, coordinates, vehicle capacity, and demand were taken from Gaskell [56] and Christofides and Eilon [57]. They used cases with nodes as many as 22 nodes (Table 1) [56], 32 nodes (Table 2) [57], and 50 nodes (Table 3) [57]. Distance (dRrsRrs+1) in rth route from node s to node s+1 is based on formula dRrsRrs+1=Xs−Xs+12+Ys−Ys+12. Meanwhile, the data of the costs and speed data were obtained from Zhang et al. [7]. The fuel cost was 7.3 yuan/liter [7]. The carbon emissions were 0.64 yuan/liter [7]. Furthermore, the vehicle usage fee was 80 yuan/hour [7]. This research employed three categories of vehicle speed (high, medium, and low speed). The high, medium, and low speeds were 107 km/hr, 63 km/hr, and 43 km/hr, respectively. Nine variations of problems (three nodes and three-speed variations) were carried out in this study. Service time for each customer is 0.1 hours. The increase in fuel consumption (p) for each additional load M = 50 is 2%.Table 1
Problem of 22 nodes.
NodeCoordinateDemandXYDepot26623501295272125230125884330926060421727450052182783006282267175724224935082302621509249268110010256267410011265257225122672423001325926525014315233500153292521501631825210017329224250182672131201927519260020303201500212082171752232618175Capacity: 4500.Table 2
Problem of 32 nodes.
NodeCoordinateDemandXYDepot29242501298427700230944540033074644004336475120053204394063214378073224372000832343390093244336001032342975011314435150012311442150133044272501429342116001529641845016261384700172974105501831540765019314406200203213914002132139830022314394130023313378700243043827502529540214002628340640002727939960028271401100029264414500302774392500312904341700323194331100Capacity: 8000.Table 3
Problem of 50 nodes.
NodeCoordinateDemandXYDepot30400137527249493035264164202695403021621471571763198316223952331110512151142411912313229135252314124221153616101652411517272331817334119131392057582821624282242578231657162485210257382826276872730481528436714295848630582719313769113238461233461023346133263562631736636963732229384535153959151440567411017274221101343564114430151645391010463239547253225482555174948281850563710Capacity: 80.


## 2.3.2. Experimental Procedure

The experiments were designed to determine the effect of HBOA (iteration and population) and speed parameters on the distribution cost and computation time. The experiments were carried out with different parameters. The parameters included the number of populations and iterations. The population parameters used three different levels (10, 50, and 100 populations). The iteration parameters also employed three levels (10, 50, and 100 iterations). This study used the sensory modality of 0.01 and power exponent of 0.1 from BOA parameters. Eighty-one trials were designed in this study. Each result of the trial was recorded for cost and computation time.The HBOA was compared to other algorithms such as BOA [48], TS [7], SA [29], ACO [28], PSO [27], and GA [58]. To compare with several algorithms, this study used one hundred iteration parameters at each vehicle speed in every algorithm. One hundred populations were used in the BOA experiment. Moreover, we used an initial temperature parameter of 1000 and the cooldown factor based on the Kuo [29] formula. One hundred ant populations were adopted for the ACO algorithm. One hundred particles and an inertia weight of 0.5 are used in the PSO algorithm. 100 populations, a crossover probability of 0.8, and a mutation probability of 0.25 were applied in the GA algorithm experiment. The performance was measured using relative error percentage (REP) as presented in equation (9). A positive REP showed that the proposed algorithm is better than the other algorithms. However, a negative REP showed that the proposed algorithm is not competitive as compared to other algorithms.(9)REP=costother algorithms−costproposed algorithmcostproposed algorithm×100%.Besides, this study also compared the computation time in all cases. It was carried out to determine the time efficiency of solving G-VRP. The effect of iteration (t) on distribution costs was also analyzed. This analysis was carried out in 50 nodes, 100 populations, and 100 iterations in the case of medium-speed vehicles. Furthermore, all experiments were conducted with the means of Matlab R2014a software on Windows 8 Intel Celeron with x64-64 2GB RAM processor.


## 3. Results and Discussion

### 3.1. The Comparison of Varied Parameters and Speed towards Costs

Table4 shows the results obtained from eighty-one experiments with variations of nodes, speed, iteration, and population. It shows that the minimum distribution cost solution is produced in the population parameters and high iterations. Therefore, overall, these results suggest that population parameters and significant iterations effectively minimize distribution costs for G-VRPs . It is interesting to note that with successive increases in both iteration and population, the distribution cost declined. It shows that the number of iterations and large population minimized the total costs. Besides, the speed of the vehicle affects the total distribution cost. Low speed requires high distribution costs. Average distribution costs are produced from medium-speed vehicles. However, high speed results in small distribution costs. Therefore, this result shows that the high speed reduced the total costs in the case of G-VRP. Overall, these findings are consistent with the findings reported by Zhang et al. [7]. Overall, these findings are consistent with Zhang et al. [7], which indicate variation, speed, iteration, and population effect in the distribution cost.Table 4
Results of the comparison between varied parameter and speed towards costs (yuan).
VrIterationHigh speedMedium speedLow speedPopulationPopulationPopulation10501001050100105010022102814.42813.12385.92971.52967.726543949.63621.73442502786.32666.92403.52909.32587.52518.53739.335773180.51002211.82161.02020.72857.22423.62300.93511.13065.7280832103392.62806.52687.4395537773261.54487.93889.53753.8502955.92732.62412.93627.33295.93068.54171.33657.53049.61002922.12647.32399.63116.62836.32614.53442.32922.22998.350107629.07398.36950.191278528.58381.7113911107810095507090.37002.66289.590908321.3798710991980690661006876.96614.36097.2888381937327.41042393338836The results of the iteration (t) effect on the distribution costs are shown in Figure 10. It illustrates the algorithm effectiveness that can be seen from the impact of iteration on the distribution costs. From the data in Figure 10, the cost of distribution decreased as the iteration was added. Besides, the effect of iteration on the distribution costs shows that the convergence curve on HBOA is better than other algorithms. The HBOA produces better total distribution in each iteration compared to BOA [48], TS [7], SA [29], ACO [28], PSO [27], and GA [58]. The results of this study indicate that the proposed algorithm is effectively used to solve the G-VRP.Figure 10
Effects of iteration (t) of each algorithm on distribution cost in the medium speed.


### 3.2. The Comparison of Varied Parameters and Speed towards Computation Time

Table5 illustrates the experimental comparison between the varied parameters and speed towards the computation time. Small populations and iterations result in less computation time. However, large populations and iterations require considerable computation time. Therefore, the most apparent finding that emerges from the analysis is that the number of iteration and large population increased the computational time. Another important finding is that different vehicle speed did not appear to affect the computation time. Low speed, medium speed, and high speed produce relatively the same computation time.Table 5
Results of the comparison between varied parameters and speed towards computation time (seconds).
VrIterationHigh speedMedium speedLow speedPopulationPopulationPopulation105010010501001050100221021.68110.05305.1921.76110.17305.5121.76110.41305.9650111.79512.59707.73113.37513.59707.91114.71515.14709.35100205.341139.501334.64206.361140.881336.20206.451141.161337.93321026.98116.47312.6727.05116.58312.9827.06116.83313.4450143.18546.20747.43144.76547.20747.61146.09548.74749.05100267.981212.401413.01269.001213.781414.57269.091214.061416.30501032.28122.88320.1432.35123.00320.4632.36123.25320.9150174.57579.81787.13176.15580.81787.31177.48582.35788.75100330.621285.301491.37331.641286.681492.93331.731286.961494.66Furthermore, the number of nodes affected the computational time. Cases with 22 nodes need little computation time. However, in the case of 50 nodes, the computation time required is considerable. Therefore, the experimental results show that the number of nodes increased the computing time. The comparison between the proposed algorithm’s computation time and several other algorithms in the medium speed is presented in Figure11. It can be seen that the proposed algorithm provided a relatively higher computation time as compared to several other algorithms, such as BOA [48], TS [7], SA [29], ACO [28], PSO [27], and GA [58]. Besides, the addition of nodes also increased the time significantly. Therefore, it can be concluded that the number of nodes has a significant effect on the computational time. It confirms the findings of Oesterle and Bauernhansl [21] and Braekers et al. [23], stating that VRP is an NP-hard problem. Based on these results, further research is expected to be carried out to reduce the computation time so that the algorithm may become more efficient.Figure 11
Comparison between the computation time of the proposed algorithm and several other algorithms in the medium speed.Although HBOA produces considerable computation time, the resulting total distribution costs are minimal. The small total cost of distribution is one of the most critical decisions in operations management. Decision-makers prefer to choose decisions with minimal total distribution costs because they provide benefits. Conversely, decision-makers pay less attention to computation time because short computing time does not guarantee a minimal total distribution cost.


### 3.3. The Comparison of Algorithms

Table6 shows the comparison of the REP values between the proposed algorithm and other algorithms. As shown in Figure 3, the REP values (based on equation (9)) of BOA [48], TS [7], SA [29], ACO [28], PSO [27], and GA [58] were 9%, 46%, 27%, 31%, 28%, and 23%, respectively. The positive values from REP indicate that the proposed algorithm is more effective in solving G-VRP. Data processing results show that there is no REP average, which results in a negative REP value in the comparison algorithm. The order of the algorithm that has the smallest to largest positive REP is BOA [48], GA [58], SA [29], PSO [27], ACO [28], and TS [7]. Thus, the findings confirm that HBOA is more competitive as compared to other algorithms. In other words, HBOA can significantly improve the quality of the G-VRP solution.Table 6
Comparison of distribution cost and the relative error percentage (REP) values between the proposed algorithm and other algorithms.
SpeedVrDistribution cost (yuan)REP (%)HBOABOATSSAACOPSOGABOATSSAACOPSOGAHigh222020.72161.02990.82355.12572.424162347.674817272016322399.62647.33567.631173107.13100.43009.3104930292925506097.26614.38694.67868.57713.777077515.984329272623Medium222300.92423.63444.72870.83240.72988.1267555025413016322614.52836.34021.53720.83775.43768.63686.285442444441507327.47987.09148.58472.383828375.28292.892516141413Low2228083065.74299.13650.43791.63757.33585.295330353428322998.33657.54875.94279.44326.34289.24193.12263434443405088369333.0117799880.610667.59972.195536331221138Average94627312823The experimental results show that the HBOA can produce a minimal total distribution cost. This result is an essential strength of the HBO algorithm. Unfortunately, there is a contradiction in the resulting computation time. The HBOA requires a relatively high computation time compared to other algorithms. However, in large nodes (50 nodes), the resulting computation time can compete with the TS algorithm. The high computation time of HBOA is caused by the large computation time TS algorithm that used to replace 10% of search agents BOA. Furthermore, the swap and flip procedure require repetition in each iteration, requiring a large computation time. In addition, this study used the sensory modality of 0.01 and power exponent of 0.1 from BOA parameters. In future investigations, it may be possible to use different sensory modality parameters and power exponent to test the quality of the solution (total distribution cost and computation time).


## 3.1. The Comparison of Varied Parameters and Speed towards Costs

Table4 shows the results obtained from eighty-one experiments with variations of nodes, speed, iteration, and population. It shows that the minimum distribution cost solution is produced in the population parameters and high iterations. Therefore, overall, these results suggest that population parameters and significant iterations effectively minimize distribution costs for G-VRPs . It is interesting to note that with successive increases in both iteration and population, the distribution cost declined. It shows that the number of iterations and large population minimized the total costs. Besides, the speed of the vehicle affects the total distribution cost. Low speed requires high distribution costs. Average distribution costs are produced from medium-speed vehicles. However, high speed results in small distribution costs. Therefore, this result shows that the high speed reduced the total costs in the case of G-VRP. Overall, these findings are consistent with the findings reported by Zhang et al. [7]. Overall, these findings are consistent with Zhang et al. [7], which indicate variation, speed, iteration, and population effect in the distribution cost.Table 4
Results of the comparison between varied parameter and speed towards costs (yuan).
VrIterationHigh speedMedium speedLow speedPopulationPopulationPopulation10501001050100105010022102814.42813.12385.92971.52967.726543949.63621.73442502786.32666.92403.52909.32587.52518.53739.335773180.51002211.82161.02020.72857.22423.62300.93511.13065.7280832103392.62806.52687.4395537773261.54487.93889.53753.8502955.92732.62412.93627.33295.93068.54171.33657.53049.61002922.12647.32399.63116.62836.32614.53442.32922.22998.350107629.07398.36950.191278528.58381.7113911107810095507090.37002.66289.590908321.3798710991980690661006876.96614.36097.2888381937327.41042393338836The results of the iteration (t) effect on the distribution costs are shown in Figure 10. It illustrates the algorithm effectiveness that can be seen from the impact of iteration on the distribution costs. From the data in Figure 10, the cost of distribution decreased as the iteration was added. Besides, the effect of iteration on the distribution costs shows that the convergence curve on HBOA is better than other algorithms. The HBOA produces better total distribution in each iteration compared to BOA [48], TS [7], SA [29], ACO [28], PSO [27], and GA [58]. The results of this study indicate that the proposed algorithm is effectively used to solve the G-VRP.Figure 10
Effects of iteration (t) of each algorithm on distribution cost in the medium speed.


## 3.2. The Comparison of Varied Parameters and Speed towards Computation Time

Table5 illustrates the experimental comparison between the varied parameters and speed towards the computation time. Small populations and iterations result in less computation time. However, large populations and iterations require considerable computation time. Therefore, the most apparent finding that emerges from the analysis is that the number of iteration and large population increased the computational time. Another important finding is that different vehicle speed did not appear to affect the computation time. Low speed, medium speed, and high speed produce relatively the same computation time.Table 5
Results of the comparison between varied parameters and speed towards computation time (seconds).
VrIterationHigh speedMedium speedLow speedPopulationPopulationPopulation105010010501001050100221021.68110.05305.1921.76110.17305.5121.76110.41305.9650111.79512.59707.73113.37513.59707.91114.71515.14709.35100205.341139.501334.64206.361140.881336.20206.451141.161337.93321026.98116.47312.6727.05116.58312.9827.06116.83313.4450143.18546.20747.43144.76547.20747.61146.09548.74749.05100267.981212.401413.01269.001213.781414.57269.091214.061416.30501032.28122.88320.1432.35123.00320.4632.36123.25320.9150174.57579.81787.13176.15580.81787.31177.48582.35788.75100330.621285.301491.37331.641286.681492.93331.731286.961494.66Furthermore, the number of nodes affected the computational time. Cases with 22 nodes need little computation time. However, in the case of 50 nodes, the computation time required is considerable. Therefore, the experimental results show that the number of nodes increased the computing time. The comparison between the proposed algorithm’s computation time and several other algorithms in the medium speed is presented in Figure11. It can be seen that the proposed algorithm provided a relatively higher computation time as compared to several other algorithms, such as BOA [48], TS [7], SA [29], ACO [28], PSO [27], and GA [58]. Besides, the addition of nodes also increased the time significantly. Therefore, it can be concluded that the number of nodes has a significant effect on the computational time. It confirms the findings of Oesterle and Bauernhansl [21] and Braekers et al. [23], stating that VRP is an NP-hard problem. Based on these results, further research is expected to be carried out to reduce the computation time so that the algorithm may become more efficient.Figure 11
Comparison between the computation time of the proposed algorithm and several other algorithms in the medium speed.Although HBOA produces considerable computation time, the resulting total distribution costs are minimal. The small total cost of distribution is one of the most critical decisions in operations management. Decision-makers prefer to choose decisions with minimal total distribution costs because they provide benefits. Conversely, decision-makers pay less attention to computation time because short computing time does not guarantee a minimal total distribution cost.


## 3.3. The Comparison of Algorithms

Table6 shows the comparison of the REP values between the proposed algorithm and other algorithms. As shown in Figure 3, the REP values (based on equation (9)) of BOA [48], TS [7], SA [29], ACO [28], PSO [27], and GA [58] were 9%, 46%, 27%, 31%, 28%, and 23%, respectively. The positive values from REP indicate that the proposed algorithm is more effective in solving G-VRP. Data processing results show that there is no REP average, which results in a negative REP value in the comparison algorithm. The order of the algorithm that has the smallest to largest positive REP is BOA [48], GA [58], SA [29], PSO [27], ACO [28], and TS [7]. Thus, the findings confirm that HBOA is more competitive as compared to other algorithms. In other words, HBOA can significantly improve the quality of the G-VRP solution.Table 6
Comparison of distribution cost and the relative error percentage (REP) values between the proposed algorithm and other algorithms.
SpeedVrDistribution cost (yuan)REP (%)HBOABOATSSAACOPSOGABOATSSAACOPSOGAHigh222020.72161.02990.82355.12572.424162347.674817272016322399.62647.33567.631173107.13100.43009.3104930292925506097.26614.38694.67868.57713.777077515.984329272623Medium222300.92423.63444.72870.83240.72988.1267555025413016322614.52836.34021.53720.83775.43768.63686.285442444441507327.47987.09148.58472.383828375.28292.892516141413Low2228083065.74299.13650.43791.63757.33585.295330353428322998.33657.54875.94279.44326.34289.24193.12263434443405088369333.0117799880.610667.59972.195536331221138Average94627312823The experimental results show that the HBOA can produce a minimal total distribution cost. This result is an essential strength of the HBO algorithm. Unfortunately, there is a contradiction in the resulting computation time. The HBOA requires a relatively high computation time compared to other algorithms. However, in large nodes (50 nodes), the resulting computation time can compete with the TS algorithm. The high computation time of HBOA is caused by the large computation time TS algorithm that used to replace 10% of search agents BOA. Furthermore, the swap and flip procedure require repetition in each iteration, requiring a large computation time. In addition, this study used the sensory modality of 0.01 and power exponent of 0.1 from BOA parameters. In future investigations, it may be possible to use different sensory modality parameters and power exponent to test the quality of the solution (total distribution cost and computation time).


## 4. Conclusion

This study discussed the green vehicle routing problem (G-VRP). The main objective of this research was to develop HBOA to minimize the distribution costs on G-VRP. This research successfully developed HBOA to solve G-VRP. The HBOA is proposed based on the BOA, which is improved with TS and local search procedures such as swap and flip. The experimental results show that the increase in population parameters and the HBOA iteration can minimize the total distribution costs. To test the algorithm performance, this algorithm was compared with several procedures. The experimental results proved that the HBOA produced a minimum total distribution cost than other algorithms. Therefore, the proposed algorithm is more competitive than the comparison algorithm. In the computation time, the results showed that the number of nodes significantly affects the computational time in HBOA. However, the proposed algorithm provides a relatively higher computation time compared to several other algorithms. Therefore, further research needs to be done to reduce the computation time so that the algorithm may become more efficient. Moreover, sensory modality and power exponent parameters need to be tested at various values. Future research should also aim at developing algorithms and problems with dynamic vehicle speeds, multidepot (distribution centers), and perishable products.



---

*Source: 8834502-2020-12-22.xml*