Solving the Two Echelon Vehicle Routing Problem Using Simulated Annealing Algorithm Considering Drop Box Facilities and Emission Cost: A Case Study of Reverse Logistics Application in Indonesia

: A two echelon distribution system is often used to solve logistics problems. This study considers a two-echelon distribution system in reverse logistics context with the use of drop box facility as an intermediary facility. An optimization model of integer linear programming is proposed, representing a two-echelon vehicle routing problem with a drop box facility (2EVRP-DF). The aim is to ﬁnd the minimum total costs consisting of vehicle transportation costs and the costs to compensate customers who have to travel to access these intermediary facilities. The results are then compared to those of common practice in reverse logistics. In common practice, customers are assumed to go directly to the depot to drop their goods. In addition, this study analyzes the environmental impact by adding a component of carbon emissions emitted by the vehicles. A set of comprehensive computational experiments is conducted. The results indicate that the 2EVRP-DF model can provide optimal costs and lower carbon emissions than the common practice.


Introduction
Indonesia is a country with a large population. Its population reached over 264 million in 2019 [1] and is projected to increase by roughly 12-15% during 2025-2035 become 296-305 million population, respectively [2]. The large number of residents has an impact on the increasing amount of waste. Three cities in Indonesia with the highest volume of waste in 2019 transported per day by type of waste are DKI Jakarta with 4139.88 m 3 , Semarang with 2755.90 m 3 , and Denpasar with 2434.74 m 3 [1]. Besides, a study by Mairizal [3] shows that electronic waste will increase from 7.3 kg/capita in 2021 to 10 kg/capita in 2040. In addition, the island of Java contributes 56% to e-waste generation. There is also an interesting case study conducted at the Galuga landfill, Bogor. This study estimated the The remainder of this paper is structured as follows. Section 2 reviews the literature related to this research. Section 3 shows the formulation of the two-echelon 2EVRP-DF model. Section 4 demonstrates the experimentation. Finally, Section 5 discusses the conclusions and future research.

Literature Review
The basic formula for two-echelon vehicle routing problem (2EVRP) variants, introduced by Perboli [13], is based on the two-echelon capacitated vehicle routing problem (2E-CVRP). The study proposed a mathematical model for the 2E-CVRP. This model was developed based on the classical capacitated vehicle routing problem (CVRP), considering that the delivery is made from the depot to the customer through an intermediary facility in a two-echelon system. Like the CVRP problem, 2E-CVRP is also NP-Hard [14]. Later, Enthoven [15] introduced a two-echelon route problem model with cover options denoted by 2EVRP-CO. The study considered cargo bicycles to deliver goods and locker facilities, becoming intermediary facilities. Furthermore, the results of the 2EVRP-CO model showed that the mileage of the cargo bike can be reduced by 60.4%. Dellaert [15] introduced a two-echelon vehicle route problem with a time window denoted by 2E-VRPTW. The 2E-VRPTW formulation separated urban vehicle routes from urban transport routes as a distinct component of the problem. The study reported that the approach successfully and efficiently solved instances with five satellites and 100 customers. Baldacci [16] proposed a 2EVRP model taking into account work constraints. The goal was to obtain a valid lower bound that breaks 2EVRP into a finite set of multi-depot vehicle routing problems. The study showed promising results with the 2EVRP which can solve the 2EVRP benchmark problem better than in earlier studies. Furthermore, the results showed that the algorithm optimally solves 144 of 153 experimental institutions.
The 2EVRP has been extended to many variants. Liu [17] developed a mathematical model and used a branch-and-cut algorithm to solve a two-echelon capacity vehicle route with grouping constraints denoted by 2E-CVRPGC. This study examined the comparison of the optimal solution between the 2E-CVRP and 2E-CVRPGC models. The results showed that the difference between the optimal solutions was quite significant, with a percentage difference of 38%. In addition, for clustered instances, the difference between 2E-CVRP and 2E-CVRPGC is 55%. Belgin [18] developed two-echelon vehicle routing problem models with simultaneous delivery and retrieval (2E-VRPSPD). In the study, a mathematical model for the problem was proposed and was strengthened using valid and The remainder of this paper is structured as follows. Section 2 reviews the literature related to this research. Section 3 shows the formulation of the two-echelon 2EVRP-DF model. Section 4 demonstrates the experimentation. Finally, Section 5 discusses the conclusions and future research.

Literature Review
The basic formula for two-echelon vehicle routing problem (2EVRP) variants, introduced by Perboli [13], is based on the two-echelon capacitated vehicle routing problem (2E-CVRP). The study proposed a mathematical model for the 2E-CVRP. This model was developed based on the classical capacitated vehicle routing problem (CVRP), considering that the delivery is made from the depot to the customer through an intermediary facility in a two-echelon system. Like the CVRP problem, 2E-CVRP is also NP-Hard [14]. Later, Enthoven [15] introduced a two-echelon route problem model with cover options denoted by 2EVRP-CO. The study considered cargo bicycles to deliver goods and locker facilities, becoming intermediary facilities. Furthermore, the results of the 2EVRP-CO model showed that the mileage of the cargo bike can be reduced by 60.4%. Dellaert [15] introduced a two-echelon vehicle route problem with a time window denoted by 2E-VRPTW. The 2E-VRPTW formulation separated urban vehicle routes from urban transport routes as a distinct component of the problem. The study reported that the approach successfully and efficiently solved instances with five satellites and 100 customers. Baldacci [16] proposed a 2EVRP model taking into account work constraints. The goal was to obtain a valid lower bound that breaks 2EVRP into a finite set of multi-depot vehicle routing problems. The study showed promising results with the 2EVRP which can solve the 2EVRP benchmark problem better than in earlier studies. Furthermore, the results showed that the algorithm optimally solves 144 of 153 experimental institutions.
The 2EVRP has been extended to many variants. Liu [17] developed a mathematical model and used a branch-and-cut algorithm to solve a two-echelon capacity vehicle route with grouping constraints denoted by 2E-CVRPGC. This study examined the comparison of the optimal solution between the 2E-CVRP and 2E-CVRPGC models. The results showed that the difference between the optimal solutions was quite significant, with a percentage difference of 38%. In addition, for clustered instances, the difference between 2E-CVRP and 2E-CVRPGC is 55%. Belgin [18] developed two-echelon vehicle routing problem models with simultaneous delivery and retrieval (2E-VRPSPD). In the study, a mathematical model for the problem was proposed and was strengthened using valid and adaptable inequalities from the literature. The study used a hybrid heuristic based environment variable descent The application of the parcel locker in logistics activities has been modelled as a network design or routing optimization problem. Deutsch and Golany [11] developed a packet locker network design problem to solve the last mile distribution problem, selecting the best number, location, and size of intermediate facilities. This study aimed to find the maximum total system profit. The total profit was made up of the revenue from the end customer using the service, deducting the fixed costs and the costs of setting up facilities, discounted shipping costs for the customers, and the potential loss of customers who do not wish to use facilities. Research conducted by Redi [10] developed a two-echelon vehicle route problem model by considering the locker facility denoted by 2E-VRPLF. The study compared the results of the 2E-VRPLF calculation with the home-delivery method, indicating that 2E-VRPLF can reduce transportation costs by 70.4%.
In solving the 2E-CVRP problem, there are various approaches, such as exact approaches, heuristic methods or a combination of both [18]. The column generation technique has been generally used in the exact approach [16]. Examples of this technique are branch and bound or branch and prices [19]. In comparison, the most heuristic method belongs to the category of metaheuristics [20]. Recent developments use an exact approach and a heuristic approach known as math-heuristics [13]. The proposed model in this study is solved using a mathematical programming technique which has embedded an exact approach. Furthermore, a simulated annealing (SA) heuristic is proposed to solve 2EVRP-DF. This is motivated by the success of SA algorithms in solving various VRP variants [10,[21][22][23]. It also addresses the inability of the exact solution approach in solving larger size instances.

The Mathematical Model
The problem is formulated as a mathematical optimization model, namely the twoechelon vehicle routing problem with drop box facilities (2EVRP-DF). This model is divided into two major parts, the first and the second echelon. The first echelon contains a collection of drop box facilities and one depot, which is represented by a set of M = {1 ... |M|}. The distance between nodes i to j is denoted as Cs ij , where i and j are part of |M|. Travel costs from the first echelon will be multiplied by a fare R, so that the total cost for the first echelon is calculated. In addition, the first echelon also adds the emission costs obtained from the total distance Cs ij multiplied by the emission coefficient E v (where v is the index of vehicles) and multiplied by the carbon tax T. In the second echelon, there is a set of customer N = {1 ... |N|}. The customer is associated with demand D n where n ∈ N. first echelon routing decision. Constraint (4) guarantees the continuity of the vehicles in the first echelon. Constraints (5) and (6) make the vehicles start from and return to the central depot after all drop box facilities have been visited. Constraints (7) and (8) limit the capacity of drop box facilities by limiting the total customer requests assigned to each facility. Constraint (9) ensures that only one drop box facility can serve customers. Constraint (10) ensures continuity of flow. Constraint (11) ensures that the vehicle capacity limit is sufficient to serve customer demand.

Transportation Fare
It is assumed that the transportation rates (R) are IDR 3000 per km. This transportation fare is taken into account as a multiplier of the distance in the first and second echelons to obtain the transportation cost per km. This tariff value is based on research conducted by Frans [24] that determined transportation fares based on Ability to Pay and Willingness to Pay.

Carbon Tax & Emission Coefficient
The model considers transportation costs incurred by companies in the first echelon and the emission costs due to pollution generated by vehicles. Ratnawati [25] formulated and proposed an initial carbon tax rate in Indonesia of IDR 80,000 per ton of CO 2 , equivalent to IDR 80 per kg of CO 2 . Thus, we use IDR 80 per kg of CO 2 as the carbon tax amount in this research. Furthermore, Kiris [26] showed that vehicle emission coefficients are related to vehicle type and fuel consumption. For example, a diesel auto 24 mpg vehicle with fuel consumption of 9.8 L/100 km has an emission coefficient of 0.2691 kg CO 2 /km. Thus, this study considers the same type of vehicle with the same emission coefficient of 0.2691 kg CO 2 /km. This study also estimates the emission coefficients of customers' vehicles at the second echelon for experimental purposes. It uses the emission coefficient proposed by Fontaras [27], which is 122.7 CO 2 g/km. Later on, this is converted in kg/km, becoming 0.1227 CO 2 kg/km.

Solving 2EVRP-DF with Simulated Annealing
This study proposed a simulated annealing algorithm to solve the 2EVRP-DF. The main characteristic of SA is its ability to get out of the local optimal by allowing the exploration of a worse solution, even a solution that is unfeasible (with a small probability). The probability of accepting a worse solution is determined by a formula that considers the gap of solution quality with the earlier solution and the parameter denoted as temperature. Several other studies have proven that SA successfully solved the problem of VRP and its variants [28][29][30].
The data used in SA to solve the 2EVRP-DF problem in this study is the same with the ones used in the GUROBI solver, such as the total number of customers, the total number of drop boxes, the total demand of each customer, and the distance between drop boxes and between each customer to each drop box facility. The objective results by the GUROBI solver are then reshaped into the form of a representation solution so that SA can initiate the solution and then perform the optimization process. Parameter settings and procedures of SA will be discussed in detail in the next subsection.

SA Parameters
In the application of the SA algorithm, there are four parameters needed. The first parameter is the initial temperature denoted by T0. The second parameter is the final temperature, which is Tf. The third parameter is Alpha which is the coefficient to control the reduction of temperature. The last parameter is Iiter, representing the number of iterations in finding a new solution at a specific temperature. These parameters need to be fine-tuned. In later discussion, the fine-tuned process for parameters of the SA algorithm shows that it can lead to a better solution quality in an acceptable computational time.

SA Procedure
In the SA algorithm, there are two phases: the initiation phase and the improvement phase. The initiation phase begins with determining the initial solution X, which is carried out randomly. In this phase also, the temperature T is set to the initial temperature T0 and the best solution of X is denoted by XBest. Then, the new solution Y is generated based on the initial solution X in the improvement phase. After obtaining a new solution Y, then Y is compared with X. If the objective value of Y is better than X, then X will be replaced by Y. On the other hand, the SA algorithm will be used to calculate formula (12), denoted by p, to help determine whether Y will replace X or not. Then, it will look for a random value between 0 to 1, symbolized as r. If the value of r is less than p, then X will be replaced with Y. At each iteration, XBest is updated if there is a new best solution found. After Iiter number of iterations, the temperature T is updated using formula (13). The SA algorithm will stop when the temperature T is less than the final temperature Tf. The overall process of SA is shown in Figure 2.

SA Neighborhood
To form a new solution Y, SA uses the neighborhood mechanism consisting of swap, insert, and reverse. Swap is a mechanism for randomly exchanging positions between two nodes in the solution. Insert is a mechanism to randomly reposition one node at a specific position in the solution. The reverse is a mechanism for randomly selecting two nodes at random and to flip the sequence of nodes between the two selected nodes.

SA Solution Representation
The solution of 2EVRP-DF is represented by a series of numbers consisting of the drop box facilities denoted by the set {1, 2, ..., |M|}; −1 dummy zeros represent vehicles returning to the depot to empty capacity, set {|M|, ..., |M| + |N|} represents the customer, and −2 dummy zeros represent the customer assignment separator to the drop box facility. Figure 3 illustrates the solution representation for 2EVRP-DF.
The procedure for finding a solution based on this representation is as follows. First, two lists are created to be introduced into a representation solution, namely a list of drop box facilities and a list of customers. The list of drop box facilities will contain the order of facilities visited by vehicles from the depot. In addition, there is a dummy node with a value of −1, which indicates that the vehicle returns to the depot to empty its vehicle capacity then goes to the next drop box facility. The second list is the customer list. This second list describes the assignment of customers to go to the specified drop box facility. In this list there are also dummy nodes with a value of −2, which are the separators between customers assignments in the second echelon. Figure 3 illustrates the representation of the solution. Take examples of data provided for solution representation in Table 1 and Table 2. In the examples illustrated, each drop box facility has a capacity of 15, demand from customers is 5, and the capacity of vehicles carrying goods is 15. As previously mentioned, the list of solutions

SA Neighborhood
To form a new solution Y, SA uses the neighborhood mechanism consisting of swap, insert, and reverse. Swap is a mechanism for randomly exchanging positions between two nodes in the solution. Insert is a mechanism to randomly reposition one node at a specific position in the solution. The reverse is a mechanism for randomly selecting two nodes at random and to flip the sequence of nodes between the two selected nodes.

SA Solution Representation
The solution of 2EVRP-DF is represented by a series of numbers consisting of the drop box facilities denoted by the set {1, 2, ..., |M|}; −1 dummy zeros represent vehicles returning to the depot to empty capacity, set {|M|, ..., |M| + |N|} represents the customer, and −2 dummy zeros represent the customer assignment separator to the drop box facility. Figure 3 illustrates the solution representation for 2EVRP-DF.
The procedure for finding a solution based on this representation is as follows. First, two lists are created to be introduced into a representation solution, namely a list of drop box facilities and a list of customers. The list of drop box facilities will contain the order of facilities visited by vehicles from the depot. In addition, there is a dummy node with a value of −1, which indicates that the vehicle returns to the depot to empty its vehicle capacity then goes to the next drop box facility. The second list is the customer list. This second list describes the assignment of customers to go to the specified drop box facility. In this list there are also dummy nodes with a value of −2, which are the separators between customers assignments in the second echelon. facility, and then to the drop box facility 1. Then, the vehicle returns to the depot to empty the vehicle and then goes to drop box 4 as the last route. Furthermore, based on the customer list, demand from customer 11 will be handled by the drop box 2 facility, demand from customer 6 will be handled by the drop box 5 facility, demand from customer 7 is handled by the drop box 1 facility, and demands from customers 8, 9, and 10 are handled by the drop box 4 facility. The drop box 3 facility is not visited by vehicles, because all customers have been served.

Results & Discussion
The 2EVRP-DF model in this study was solved using the GUROBI solver (a MILP solver), coded in the AMPL mathematical programming language (AMPL/GUROBI). The experiment was conducted on a PC with an Intel(R) Core(TM) i5 processor at 3.60 GHz, 8 Gb RAM, and running on a Windows 10 64-bit operating system platform. The computational experiments aimed to demonstrate the improvement benefit of considering 2EVRP-DF and the effectiveness of the proposed SA in resolving 2EVRP-DF. The following subsections describe the test instances, parameter settings, and computational results.

Test Instances
Testing data was taken from a drop box dataset located around the Jakarta area [31], and reverse geolocation is carried out using the Google API to obtain its latitude and longitude. For customer locations, uniform random distribution using python programming is carried out around the drop box locations within a radius of 10 km. For small instances, the number of drop box facilities ranges from 5 to 30 and the number of customers ranges from 6 to 25. Meanwhile, for large agencies, the number of drop box facilities is 92, and the number of customers ranges from 70 to 150. After obtaining the drop box facilities and customers' location, the distance between the drop box facility and the customer is calculated using the haversine formula.
Furthermore, as previously mentioned regarding common practice in reverse logistics, we also calculate the distance from the customer to the depot using the haversine formula. We assume that the R rate in the second echelon is also the same as in the first echelon, so to get the cost in common practice we multiply the customer's distance by the rate of Rp3000 per km.

Parameter Setting
In parameter settings, we used the one factor at a time (OFAT) procedure to set one parameter sequentially at a time. The OFAT procedure proved to be effective for selecting parameters for the SA algorithm [32]. In this experiment, each parameter has several options to choose from as follows.  Figure 4b shows that the obtained Tf is equal to 0.5, in order to obtain a minimum objective value, also in an acceptable time. For the Alpha and Iiter parameters the same procedure is used as T0 and Tf. After setting the parameters, from the experimental results, the values of T0 = 90, Tf = 0.5, Alpha = 0.99, and Iiter = 100N. After all these processes, the final parameters are as follows: T0 = 90, Tf = 0.5, Alpha = 0.99, and Iiter = 100N. Figure 4 shows the analysis of the effect of each SA parameter on the objective value and computation time.

Computational Results
The computational experiment was divided into four parts. The first part provides an analysis of the convergence property of the SA algorithm. The second part is the experiment with small size instances of 2EVRPD-DF. The third part is the experiment with large size instances of 2EVRPD-DF. In the second and third part of the experiment, the results obtained by GUROBI, SA, and common practice are compared. Finally, the last part provides a sensitivity analysis of the impact of the increasing number of customers on the emission cost, considering the use of 2EVRPD-DF (result of GUROBI and SA) compared to common practice.
For the first part of the experiment, the small size instance was used, containing 5 to 30 customers and 6 to 25 drop box facilities. The solution comparison of the 2EVRP-DF problem is between AMPL and the proposed SA algorithm. The numerical results are shown in Table 3 Columns 1 & 2 specify the number of customers, and the number of drop box facilities at each instance. Columns 3 to 9 are the results of the calculation of the 2EVRP-DF model, solved using the GUROBI. Columns 3 and 4, respectively, contain the total of the distances in the first and second echelon. Columns 5 and 6, respectively, contain the total emissions generated from the first and second echelon. Column 7 contains the total transportation costs obtained from the sum of the total distances of the first and second echelon, which is then multiplied by the transportation fare. Column 8 contains the total emission costs obtained from the total emissions produced in the first and second echelon, which is then multiplied by the carbon tax. Column 9 is the total cost of transportation and emission costs. The next seven columns are the solutions generated by the SA, each of which has the same definition as the result of the GUROBI solver calculation. The

Computational Results
The computational experiment was divided into four parts. The first part provides an analysis of the convergence property of the SA algorithm. The second part is the experiment with small size instances of 2EVRPD-DF. The third part is the experiment with large size instances of 2EVRPD-DF. In the second and third part of the experiment, the results obtained by GUROBI, SA, and common practice are compared. Finally, the last part provides a sensitivity analysis of the impact of the increasing number of customers on the emission cost, considering the use of 2EVRPD-DF (result of GUROBI and SA) compared to common practice.
For the first part of the experiment, the small size instance was used, containing 5 to 30 customers and 6 to 25 drop box facilities. The solution comparison of the 2EVRP-DF problem is between AMPL and the proposed SA algorithm. The numerical results are shown in Table 3 Columns 1 & 2 specify the number of customers, and the number of drop box facilities at each instance. Columns 3 to 9 are the results of the calculation of the 2EVRP-DF model, solved using the GUROBI. Columns 3 and 4, respectively, contain the total of the distances in the first and second echelon. Columns 5 and 6, respectively, contain the total emissions generated from the first and second echelon. Column 7 contains the total transportation costs obtained from the sum of the total distances of the first and second echelon, which is then multiplied by the transportation fare. Column 8 contains the total emission costs obtained from the total emissions produced in the first and second echelon, which is then multiplied by the carbon tax. Column 9 is the total cost of transportation and emission costs. The next seven columns are the solutions generated by the SA, each of which has the same definition as the result of the GUROBI solver calculation. The last column is the difference between the total cost generated by the GUROBI solver and the solution generated by SA.
The numerical results using the GUROBI solver have shown optimal solution. Furthermore, the numerical results of the proposed SA algorithm have also demonstrated the same results as GUROBI, where the difference in total costs between GUROBI and SA is 0. This indicates that the developed SA algorithm has been validated, because the results given are not smaller or larger than the optimal solution generated by GUROBI.
The second part of the experiment was conducted with the large size instance which has a number of customers of from 70 to 150 and a number of drop box facilities of as many as 92. Each customer has a demand at random value between 1-15 kg, and each drop box facility has a static capacity of 180 kg. Table 4 shows a comparison of the total costs between the 2EVRP-DF model (generated from the GUROBI solver and SA algorithm), and common practice where customers go directly to the depot to drop their items. The total cost generated from the 2EVRP-DF model provides better results than common practice, with an average of 71.95% and 72.09%. This result shows that using the 2EVRP-DF model can provide an average cost advantage of 72.02% compared to the common practice. Furthermore, the SA algorithm gives slightly better results of total costs compared to GUROBI with an average difference of only 0.51%. This is possible because the calculations performed by the GUROBI solver are near optimal. On the other hand, from the computational time point of view, the SA algorithm gives better results than GUROBI with an average computation time of 298.24 s, equivalent to 4.97 min, where the average computational time of GUROBI reaches 23,409.36 s, equivalent to approximately 6.5 h.
In the third part of the experiment, a convergency analysis of the SA algorithm is similar to that of previous studies [33,34]. We compared the convergency of the SA algorithm with other metaheuristics algorithms, in this case a Genetic Algorithm (GA). The instances are selected from the large size dataset which has 70 customers and 92 drop box facilities. Figure 5 shows the convergency of the SA and GA algorithm in finding a solution at a predetermined time, from 0 s to 3000 s. From 0 to 300 s, it is shown that both the SA and GA algorithm is able to search for minimal solutions. For SA, there is a drastic change in the solution at 240, with a total objective of 618,378 and 616,670 in 630 s. Later in the iteration, the SA algorithm made good improvements, but progressing slowly as it approached steady-state. At the final part of the iteration, the SA algorithm shows the same solution until the end of the experiment and shows a convergence property. Meanwhile, GA shows a slower process than SA. The improvement of the solution in GA is shown at 840 and 1890 s with the objective value of 620,123 and 619,556, respectively. At the final part of the iteration, GA also shows the same solution until the end of the experiment and shows a convergence property.
The parameters required by GA are population size, crossover rate, and mutation rate, which are 2000, 0.6, and 0.01, respectively. This means that the probability of population crossover mechanism is 60% and the chance of mutation mechanism is 1%.
The comparison of total objective of SA and GA is shown in Table 5. Both SA and GA demonstrate a good optimization process in finding a minimum solution with the 2EVRP-DF model. SA gives an average total objective of 887,762.69, while GA provides 889,652.29. Although the average difference in the objective values looks small, in solving the 2EVRP-DF model SA outperforms GA by 0.21%, which is in line with several previous studies where SA was able to outperform other metaheuristics algorithms, especially in VRP problems [35][36][37][38]. The comparison of total objective of SA and GA is shown in Table 5. Both SA and GA demonstrate a good optimization process in finding a minimum solution with the 2EVRP-DF model. SA gives an average total objective of 887,762.69, while GA provides 889,652.29. Although the average difference in the objective values looks small, in solving the 2EVRP-DF model SA outperforms GA by 0.21%, which is in line with several previous studies where SA was able to outperform other metaheuristics algorithms, especially in VRP problems [35][36][37][38].
In a further analysis regarding emission costs, the computational results from the 2EVRP-DF model were promising compared to common practice. Table 6 shows the detailed emission calculations for the 2EVRP-DF model and common practice. The 2EVRP-DF model shows that the average total emissions produced by vehicles in the first and second echelon are 47.80 (GUROBI) and 47.56 (GUROBI), where the total emission in common practice averages 133.60. Furthermore, Figure 6 shows that, as more customers are added to the system, the amount of carbon emissions also increases. However, when compared to the 2EVRP-DF model, common practice generates much higher emission costs compared to the 2EVRP-DF. For this reason, it can be concluded that the 2EVRP-DF model can have a positive impact on the environment by reducing the amount of carbon emissions released by vehicles.    added to the system, the amount of carbon emissions also increases. However, when compared to the 2EVRP-DF model, common practice generates much higher emission costs compared to the 2EVRP-DF. For this reason, it can be concluded that the 2EVRP-DF model can have a positive impact on the environment by reducing the amount of carbon emissions released by vehicles. Figure 6. Correlation between total carbon emission and customers.

Conclusions
This study proposes a two-echelon vehicle routing problem with a drop box facility (2EVRP-DF) which is modeled as an integer linear programming problem. The 2EVRP-DF demonstrates the utilization of drop box facilities to carry out reverse logistics activities in urban areas. This study also considers environmental aspects by considering the amount of carbon emissions released by vehicles, and compares the results of 2EVRP-DF with a common practice model where, in reverse logistics activities, customers go directly to the depot rather than intermediate facilities. Numerical examples are carried out to show the effectiveness of using the drop box facility compared to common practice. The use of drop boxes in the 2EVRP-DF model allows an average reduction in total costs (transportation and emissions) of 72.02%. In addition, it is also seen that there is a drastic increase in the number of emissions released by vehicles in common practice, where this does not occur in the 2EVRP-DF model. Furthermore, the experiment using SA to solve the 2EVRP-DF problem model presents less significant results compared to GUROBI. The average difference in solutions between GUROBI and SA is only 0.51%, but in terms of time SA can provide a more acceptable computation time compared to GUROBI which reaches 6.5 h. However, the computational results still need to be tested on more complex characteristics of problems that might occur in practice or real life. In the comparison between SA and GA in solving 2EVRP-DF, SA also outperforms GA in finding a minimum objective solution by an average difference of 0.21%. 35 Colleration of total emission and customers SA GUROBI Common Practice Figure 6. Correlation between total carbon emission and customers.

Conclusions
This study proposes a two-echelon vehicle routing problem with a drop box facility (2EVRP-DF) which is modeled as an integer linear programming problem. The 2EVRP-DF demonstrates the utilization of drop box facilities to carry out reverse logistics activities in urban areas. This study also considers environmental aspects by considering the amount of carbon emissions released by vehicles, and compares the results of 2EVRP-DF with a common practice model where, in reverse logistics activities, customers go directly to the depot rather than intermediate facilities. Numerical examples are carried out to show the effectiveness of using the drop box facility compared to common practice. The use of drop boxes in the 2EVRP-DF model allows an average reduction in total costs (transportation and emissions) of 72.02%. In addition, it is also seen that there is a drastic increase in the number of emissions released by vehicles in common practice, where this does not occur in the 2EVRP-DF model. Furthermore, the experiment using SA to solve the 2EVRP-DF problem model presents less significant results compared to GUROBI. The average difference in solutions between GUROBI and SA is only 0.51%, but in terms of time SA can provide a more acceptable computation time compared to GUROBI which reaches 6.5 h. However, the computational results still need to be tested on more complex characteristics of problems that might occur in practice or real life. In the comparison between SA and GA in solving 2EVRP-DF, SA also outperforms GA in finding a minimum objective solution by an average difference of 0.21%.
For further research, it is recommended to explore various metaheuristic algorithms to solve 2EVRP-DF problems, because there are growing number of developments regarding metaheuristic algorithms available in the literature. It will also be more interesting if a combination of algorithms such as a hybrid algorithm between SA and GA can solve the problem. In term of problem characteristics, it is interesting to consider uncertainty factors such as travel time, demand, and locker availability in the model.