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

Abstract

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 find 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.

1. 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³, Semarang with 2755.90 m³, and Denpasar with 2434.74 m³ [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 average daily removal of microplastics from the leachate channel at 80,640 ± 604.80 particles. The daily release of mesoplastic from leachate channels into the aquatic environment was estimated at 618,240 ± 1905.45 particles [4].

With regard to the increasing level of waste, awareness of the importance of waste management to create a sustainable environment is increasing as well, particularly in the context of recycling process or waste disposal for products that cannot be recycled [5]. Reverse Logistics (RL) plays an essential role in the success of waste management [6]. Reverse logistics is a part of the supply chain, where the activities are reversed from common logistics. RL consists of reverse product flow activities, starting with the final customer until the goods are returned to the manufacturer. RL goals can vary from proper dumping at the end of a product’s life cycle, renewals, product recollections, to being part of any eco-friendly business program. Besides maintaining good environmental sustainability, RL is also beneficial for the company’s image. It is the trade-off between the costs incurred to support the process and the benefit generated by the actions.

For RL in urban areas, end customers usually returns goods to sellers or retailers directly. This common practice often encounters problems that are not profitable for customers, such as traffic congestion and environmental issues such as air pollution. In addition, if the distance between the customer and the seller is far, then the shipping cost of the goods to be returned is charged to the customer who returns the goods.

An alternative to the common practice for RL activities for customer goods is the pick-up method. This can save time and is more cost-effective than the common practice given correct planning. In the pick-up method, the customer can pick up or send the package through an intermediary facility near the customer’s residence or office. The customer’s pick-up location is usually in a communal area that is easily accessible. A more conventional form of parcel locker as the intermediary facility is seen in the application of the receiving box, which is a fixed locker for receiving packages. It is commonly installed outside the consumer’s house. In a study conducted by Lemke [7], almost 15% of customers tend to be more likely to use lockers if the locker location is close to their home. Many recent studies on logistics activity in general are concerned not only with the economic perspective but also with the sustainable and green aspects of the operations [8,9]. Thus, consideration of emission costs is often included when discussing sustainable and green logistics activity. Furthermore, studies regarding the utilization of parcel lockers to enhance reverse logistics are also gaining increasing interest from researchers [10,11,12]. However, among the research available on parcel lockers applied in reverse logistics, to the best of our knowledge no available study proposes the analysis of using parcel lockers or drop box facilities for reverse logistics considering emission cost.

This study discusses reverse logistics activities such as product returns or refurbishment by utilizing a network of drop box facilities. The problem is formulated as an optimization model denoted as the two-echelon vehicle routing problem with drop box facilities (2EVRP-DF). The 2EVRP-DF problem model is used to find the minimum value of transportation costs that need to be incurred by vehicles that pick-up goods at intermediary facilities, and travel costs that must be incurred by customers who travel to intermediary facilities. In addition to calculating the cost of transportation, we also consider environmental factors, namely by adding the carbon tax component and the emission coefficient to the problem model, which is the main contribution of this study. Figure 1 illustrates the reverse logistics considering the drop box facilities.

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.

2. 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 and a local search algorithm called VND-LS. The result showed effectiveness in finding a good solution for 2EVRPSPD.

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.

3. Methods

3.1. The Mathematical Model

The problem is formulated as a mathematical optimization model, namely the two-echelon 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\in N$. The assignment of customers is limited by the capacity of the drop box facility Qs_m, where $m\in \left\{2\dots \left|M\right|\right\}$. The vehicle capacity Qv_v also limits route decisions to the first echelon where v is one of the vehicles, which is also represented by the set V = {1 ... |V|}. The following are the decision variables contained in the 2EVRP-DF model.

3.1.1. Decision Variables

• X_i,j,v: equals 1 if the route from i to j is made, otherwise 0
• Z_i: equals 1 if the drop box facility i is open, otherwise 0
• S_i,j: equals 1 if the drop box facility i serving the customer j, otherwise 0
• Q_i: The total amount of demand served in drop box facility i
• a_i,v: The total amount demand in drop box facility i carried by the vehicle v

3.1.2. Problem Formulation

Formula (1) defines the objective function, which aims to find a solution with the minimum total cost, consists of the first echelon’s transportation costs, the emission costs, and the costs that the customer must incur regarding the drop box facility in the second echelon. Constraints (2) and (3) ensure that the drop box facility is serviced once in the 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.

$$\begin{gathered} \min \mathrm{R}{\displaystyle \sum }_{i=1}^{\left|M\right|}{\displaystyle \sum }_{j=i}^{\left|M\right|}{\displaystyle \sum }_{v=1}^{\left|V\right|}C{s}_{i,j}·X_{i,j,v}+\mathrm{T}{\displaystyle \sum }_{i=1}^{\left|M\right|}{\displaystyle \sum }_{j=i}^{\left|M\right|}{\displaystyle \sum }_{v=1}^{\left|V\right|}C{s}_{i,j}·X_{i,j,v}·E_v+\mathrm{R}{\displaystyle \sum }_{i=1}^{\left|N\right|}{\displaystyle \sum }_{j=2}^{\left|M\right|}C{c}_{i,j}·s_{i,j} \\ +T {\displaystyle \sum }_{i=1}^{\left|N\right|}{\displaystyle \sum }_{j=2}^{\left|M\right|}C{c}_{i,j}·s_{i,j}· E2 \end{gathered}$$

(1)

$${\displaystyle {\displaystyle \sum }_{i=1}^{\left|M\right|}{\displaystyle \sum }_{v=1}^{\left|V\right|}}{X}_{i,j,v}=z_j\forall j\in \left\{2\dots \left|M\right|\right\}$$

(2)

$${{\displaystyle \sum }}_{j=1}^{\left|M\right|}{{\displaystyle \sum }}_{v=1}^{\left|V\right|}{X}_{i,j,v}=z_i\forall i\in \left\{2\dots \left|M\right|\right\}$$

(3)

$${{\displaystyle \sum }}_{i=1}^{\left|M\right|}{X}_{i,l,v}-{{\displaystyle \sum }}_{j=1}^{\left|M\right|}{X}_{l,j,v}=0\forall l\in \left\{2\dots \left|M\right|\right\},\forall v\in \left\{1\dots \left|V\right|\right\}$$

(4)

$${{\displaystyle \sum }}_{j=2}^{\left|M\right|}{X}_{1,j,v}\le 1\forall v\in \left\{1\dots \left|V\right|\right\}$$

(5)

$${{\displaystyle \sum }}_{i=2}^{\left|M\right|}{X}_{i,1,v}\le 1\forall v\in \left\{1\dots \left|V\right|\right\}$$

(6)

$${{\displaystyle \sum }}_{i=1}^{\left|N\right|}{D}_i·s_{i,j}\le Q_j\forall j\in \left\{2\dots \left|M\right|\right\}$$

(7)

$$Q_j\le Q{s}_{j }·z_j\forall j\in \left\{2\dots \left|M\right|\right\}$$

(8)

$${{\displaystyle \sum }}_{j=2}^{\left|M\right|}{s}_{ij}=1\forall i\in \left\{1\dots \left|N\right|\right\}$$

(9)

$$a_{i,v}+Q_j-L.\left(1-x_{i,j,v}\right)\le a_{j,v}\forall i\in \left\{1\dots \left|M\right|\right\},\forall j\in \left\{2\dots \left|M\right|\right\},\forall v\in \left\{1\dots \left|V\right|\right\}$$

(10)

$${{\displaystyle \sum }}_{i=1}^{\left|M\right|}{{\displaystyle \sum }}_{j=2}^{\left|M\right|}{x}_{i,j,v}·Q_j\le Q{v}_v\forall v\in \left\{1\dots \left|V\right|\right\}$$

(11)

3.1.3. 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.

3.1.4. 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₂, equivalent to IDR 80 per kg of CO₂. Thus, we use IDR 80 per kg of CO₂ 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₂/km. Thus, this study considers the same type of vehicle with the same emission coefficient of 0.2691 kg CO₂/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₂ g/km. Later on, this is converted in kg/km, becoming 0.1227 CO₂ kg/km.

3.2. 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.

3.2.1. 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.

3.2.2. 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.

$$p=exp\left(\frac{f\left(y\right)-f\left(x\right)}{T}\right)$$

(12)

$$T=T\ast Alpha$$

(13)

3.2.3. 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.

3.2.4. 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. Example of location, capacity, and demand for 2EVRP-DF.

Node Facilities & Customers	Latitude	Longitude	Capacity	Demand
0	−6.2701	106.8376	-	-
1	−6.25631	106.8126	15	-
2	−6.25693	106.8518	15	-
3	−6.25465	106.8138	15	-
4	−6.28424	106.8118	15	-
5	−6.25153	106.8247	15	-
6	−6.24495	106.8287	-	5
7	−6.26012	106.8155	-	5
8	−6.2746	106.8212	-	5
9	−6.28718	106.8012	-	5
10	−6.27253	106.8207	-	5
11	−6.24853	106.8441	-	5

and

Table 2. Illustration of cost matrix between nodes in 2EVRP-DF.

Nodes	0	1	2	3	4	5	6	7	8	9	10	11
0	0	3.16	2.15	3.14	3.26	2.51	2.96	2.68	1.88	4.45	1.89	2.5
1	3.16	0	4.33	0.23	3.11	1.43	2.18	0.53	2.24	3.66	2.01	3.58
2	2.15	4.33	0	4.21	5.37	3.06	2.88	4.03	3.92	6.53	3.86	1.27
3	3.14	0.23	4.21	0	3.3	1.24	1.96	0.64	2.36	3.88	2.13	3.41
4	3.26	3.11	5.37	3.3	0	3.91	4.75	2.71	1.49	1.22	1.63	5.34
5	2.51	1.43	3.06	1.24	3.91	0	0.86	1.39	2.59	4.74	2.38	2.17
6	2.96	2.18	2.88	1.96	4.75	0.86	0	2.23	3.4	5.6	3.19	1.75
7	2.68	0.53	4.03	0.64	2.71	1.39	2.23	0	1.73	3.4	1.49	3.41
8	1.88	2.24	3.92	2.36	1.49	2.59	3.4	1.73	0	2.62	0.24	3.85
9	4.45	3.66	6.53	3.88	1.22	4.74	5.6	3.4	2.62	0	2.7	6.4
10	1.89	2.01	3.86	2.13	1.63	2.38	3.19	1.49	0.24	2.7	0	3.72
11	2.5	3.58	1.27	3.41	5.34	2.17	1.75	3.41	3.85	6.4	3.72	0

. 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 {2,5,1,−1,4,11,−2,6,−2,7,−2,8,9,10} is divided into two, namely the list of drop box facilities and the list of customers. The list of drop box facilities will contain {2,5,1,−1,4}. On the other hand, the customer list is {11,−2,6,−2,7,−2,8,9,10}. Based on an example list of drop box facilities and customers, the following solutions can be explained. The first facility visited by the vehicle is the drop box facility 2. Then the vehicle goes to the drop box 5 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.

4. 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.

4.1. 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.

4.2. 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.

• Initial Temp. (T0): 30, 60, 90
• Final Temp. (Tf): 5, 0.5, 0.05
• Alpha: 0.25, 0.5, 0.9, 0.99
• Iiter: 25N, 50N, 100N, 250N, 300N (N = number of customers)

Figure 4a shows the experimental determination of the T0 parameter, while the other parameters are fixed. The T0 determination experiment shows the number 90 in the search for the minimum objective value in an acceptable time. Then the search for the Tf parameter was carried out with the value T0 = 90, and the other parameters were fixed. 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.

4.3. 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. Comparison of results between GUROBI & SA on small size instance.

Customer (N)	Total Demand (kg)	Drop Box (M)	2EVRP-DF—GUROBI							2EVRP-DF—SA							Difference
			E1 (km)	E2 (km)	Emission E1 (kg CO2)	Emission E2 (kg CO2)	Transport Cost (IDR)	Total Emission Cost (IDR)	Total Cost (IDR)	E1 (km)	E2 (km)	Emission E1 (kg CO2)	Emission E2 (kg CO2)	Transport Cost (IDR)	Total Emission Cost (IDR)	Total Cost (IDR)
5	40	25	3.80	5.23	9.03	1.02	27,090	82	27,172	3.80	5.23	1.02	0.64	27,090	133	27,223	0
10	73	25	7.70	7.49	15.19	2.07	45,570	166	45,736	7.70	7.49	2.07	0.92	45,570	240	45,810	0
15	95	25	7.88	10.22	18.10	2.12	54,300	170	54,470	7.88	10.22	2.12	1.25	54,300	270	54,570	0
20	124	25	10.62	13.16	23.78	2.86	71,340	229	71,569	10.62	13.16	2.86	1.61	71,340	358	71,698	0
25	143	25	10.82	15.29	26.11	2.91	78,330	233	78,563	10.82	15.29	2.91	1.88	78,330	383	78,713	0
30	174	6	12.80	41.04	30.28	3.44	90,840	275	91,115	13.05	41.04	3.51	5.04	162,270	684	162,954	0
30	174	11	13.05	31.46	54.09	3.51	162,270	281	162,551	13.42	31.46	3.61	3.86	134,640	598	135,238	0
30	174	16	13.42	29.29	44.88	3.61	134,640	289	134,929	13.56	29.29	3.65	3.59	128,550	580	129,130	0
30	174	21	13.56	18.27	42.85	3.65	128,550	292	128,842	13.65	18.27	3.67	2.24	95,760	473	96,233	0
30	174	25	13.65	17.48	31.92	3.67	95,760	294	96,054	12.80	17.48	3.44	2.14	90,840	448	91,288	0

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. Comparison of results between 2EVRP-DF and common practice on large size instance.

Customer (N)	Total Demand (kg)	Drop Box (M)	2EVRP-DF—GUROBI				2EVRP-DF—SA				Common Practice (CP)			Cost Diff. GUROBI—CP (%)	Cost Diff. SA—CP (%)	Cost Diff. GUROBI—SA (%)	Time Diff. GUROBI—SA (s)
			Transport Cost (Rp)	Emission Cost (IDR)	Total Cost (IDR)	Solved Time (s)	Transport Cost (Rp)	Emission Cost (IDR)	Total Cost (IDR)	Solved Time (s)	Transport Cost (IDR)	Total Emission Cost (IDR)	Total Cost (IDR)
70	527	92	618,894	2833	621,727	23,822.00	613,822	2825	616,647	135.20	1,817,917	5948	1,823,866	65.91	66.19	0.82	23,686.80
80	630	92	712,940	3219	716,159	27,824.40	712,044	3209	715,253	163.00	2,240,801	7332	2,248,133	68.14	68.18	0.13	27,661.40
90	704	92	763,551	3391	766,942	26,585.00	762,265	3386	765,651	229.60	2,551,443	8348	2,559,792	70.04	70.09	0.17	26,355.40
100	786	92	814,160	3549	817,709	27,099.00	812,360	3529	815,889	246.90	2,929,518	9585	2,939,103	72.18	72.24	0.22	26,852.10
110	873	92	874,412	3734	878,146	24,057.20	871,158	3727	874,885	287.70	3,242,748	10,610	3,253,358	73.01	73.11	0.37	23,769.50
120	964	92	951,347	4050	955,397	23,137.00	951,347	4050	955,397	314.70	3,647,314	11,934	3,659,248	73.89	73.89	-	22,822.30
130	1035	92	1,010,491	4262	1,014,753	13,455.20	1,008,732	4239	1,012,971	393.90	3,951,341	12,929	3,964,270	74.40	74.45	0.18	13,061.30
140	1101	92	1,095,068	4568	1,099,636	18,812.40	1,078,093	4502	1,082,595	427.30	4,344,548	14,215	4,358,763	74.77	75.16	1.55	18,385.10
150	1194	92	1,159,221	4812	1,164,033	25,892.00	1,145,798	4776	1,150,574	485.90	4,673,080	15,290	4,688,371	75.17	75.46	1.16	25,406.10
Average			888,898	3824	892,722	23,409.36	883,958	3805	887,763	298.24	3,266,523	10,688	3,277,212	71.95	72.09	0.51	23,111.11

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. Comparison of total objective between SA and GA.

Customer (N)	Total Demand (kg)	Drop Box (M)	2EVRP-DF—SA							2EVRP-DF—GA							Difference
			E1 (km)	E2 (km)	Emission E1 (kg CO2)	Emission E2 (kg CO2)	Transport Cost (IDR)	Total Emission Cost (IDR)	Total Cost (IDR)	E1 (km)	E2 (km)	Emission E1 (kg CO2)	Emission E2 (kg CO2)	Transport Cost (IDR)	Total Emission Cost (IDR)	Total Cost (IDR)
70	527	92	69.78	134.82	18.78	16.54	613,822	1502	616,647	69.06	136.52	18.58	16.75	616,729	2827	619,556	0.47
80	630	92	75.13	162.22	20.22	19.90	712,044	1617	715,253	75.74	161.91	20.38	19.87	712,940	3220	716,159	0.13
90	704	92	76.20	177.89	20.50	21.83	762,265	1640	765,651	76.27	177.86	20.52	21.82	762,393	3388	765,781	0.02
100	786	92	74.38	196.41	20.01	24.10	812,360	1601	815,889	75.52	195.87	20.32	24.03	814,159	3548	817,708	0.22
110	873	92	74.88	215.51	20.15	26.44	871,158	1612	874,885	74.56	216.52	20.06	26.57	873,254	3731	876,985	0.24
120	964	92	80.02	237.09	21.53	29.09	951,347	1723	955,397	80.02	237.09	21.53	29.09	951,346	4050	955,396	0.00
130	1035	92	80.15	256.10	21.57	31.42	1,008,732	1725	1,012,971	80.22	256.48	21.59	31.47	1,010,112	4245	1,014,357	0.14
140	1101	92	83.25	276.11	22.40	33.88	1,078,093	1792	1,082,595	84.14	276.13	22.64	33.88	1,080,807	4522	1,085,329	0.25
150	1194	92	87.68	294.25	23.60	36.10	1,145,798	1888	1,150,574	87.17	296.43	23.46	36.37	1,150,814	4786	1,155,601	0.43
Average			77.94	216.71	20.97	26.59	883,957.67	3805.02	887,762.69	78.08	217.20	21.01	26.65	885,839.39	3812.90	889,652.29	0.21

. 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. Emission cost comparison between 2EVRP-DF and common practice.

Customer (N)	Total Demand (kg)	Drop Box (M)	2EVRP-DF—GUROBI						2EVRP-DF—SA						Common Practice (CP)
			E1 (km)	E2 (km)	Emmision E1 (kg CO2)	Emmision E2 (kg CO2)	Total Emission (kg CO2)	Emission Cost (IDR)	E1 (km)	E2 (km)	Emmision E1 (kg CO2)	Emmision E2 (kg CO2)	Total Emission (kg CO2)	Emission Cost (IDR)	Transport Cost (IDR)	Total Emission Cost (IDR)	Total Cost (IDR)
70	527	92	69.03	137.27	18.58	16.84	35.42	2833	69.78	134.82	18.78	16.54	35.32	2825	605.97	74.35	5948
80	630	92	75.74	161.91	20.38	19.87	40.25	3219	75.13	162.22	20.22	19.90	40.12	3209	746.93	91.65	7332
90	704	92	76.24	178.28	20.52	21.87	42.39	3391	76.20	177.89	20.50	21.83	42.33	3386	850.48	104.35	8348
100	786	92	75.52	195.87	20.32	24.03	44.35	3549	74.38	196.41	20.01	24.10	44.11	3529	976.51	119.82	9585
110	873	92	74.53	216.94	20.06	26.62	46.68	3734	74.88	215.51	20.15	26.44	46.59	3727	1080.92	132.63	10,610
120	964	92	80.02	237.09	21.53	29.09	50.63	4050	80.02	237.09	21.53	29.09	50.63	4050	1215.77	149.18	11,934
130	1035	92	81.56	255.27	21.95	31.32	53.27	4262	80.15	256.10	21.57	31.42	52.99	4239	1317.11	161.61	12,929
140	1101	92	84.09	280.93	22.63	34.47	57.10	4568	83.25	276.11	22.40	33.88	56.28	4502	1448.18	177.69	14,215
150	1194	92	86.95	299.46	23.40	36.74	60.14	4812	87.68	294.25	23.60	36.10	59.70	4776	1557.69	191.13	15,290
Average			78.19	218.11	21.04	26.76	47.80	3824	77.94	216.71	20.97	26.59	47.56	3805	1088.84	133.60	10,688

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.

5. 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.