Optimization of LPG Distribution Route Using Variable Neighborhood Tabu Search Algorithm

—PT. Galaxi Energi Pratama (GEP) is one of the biggest distributors of subsidized LPG in Malang Raya area. Currently the route planning is not done very well, which results in a high fuel cost. With the company's main business process being distribution, the planning needs to be improved to maximize the profit. The problem in PT. GEP is classified as the Heterogeneous Vehicle Routing Problem with Multiple Trips (HVRPM). This problem is classified as NP-Hard and requires high computational effort to obtain a good solution so metaheuristic method is preferred. In this research, variable neighborhood tabu search (VNTS) algorithm is developed to solve the HVRPM and implemented to minimize the fuel cost of PT. GEP. The developed algorithm is implemented in the six instances collected from the case study. The generated trips produce a total savings of Rp 150,876 for one operational week, or roughly 18% of the initial cost. The computation time of the algorithm is evaluated by comparing with Simulated Annealing using a problem with the same size. VNTS has a lower average time and is expected to perform competitively when a standardized dataset is used for comparison. The solution quality of the algorithm is then compared with branch-and-bound method. VNTS is able to find one global optimal solution out of the six instances and overall, it performs better than branch-and-bound.


I. INTRODUCTION
ISTRIBUTION is the process of making products or services available for business users or customers. It is an integral part of a company's business process as without it, products cannot reach the customer. A well-planned distribution route can cut fuel expenses and reduce vehicle wear-out, ultimately saving maintenance expenses of the vehicles. A shorter delivery route also means that the product can be received by customer faster, thus improving customer satisfaction.
LPG has become a very common and important product today. Its distribution system needs to be planned well to avoid its absence and keep its price affordable. In Indonesia, LPG is produced by a government-owned company called PT. Pertamina. There are two types of LPG that are sold, which are subsidized and non-subsidized. The subsidized LPG has net weight of 3 kg and is sold only to society from the poor class. Its supply chain consists of PT. Pertamina as the producer, LPG agents, retailers, and the final consumers. Because of the subsidies it receives and its specific market, the distribution process is regulated by the government in terms of its selling price and the quantity that an agent/retailer distributes/sells.
PT. GEP is one of the biggest LPG distributors in Malang Raya area. It distributes 3-kg LPG to retailers that have been registered to PT. Pertamina. Because the product is subsidized, the delivery quantity for each customer varies every day based on a monthly schedule released by Pertamina. Currently, the distribution planning is done by the drivers themselves. The allocation of customers is based on a mutual agreement on the customers that they will regularly visit. This fixed agreement may not be optimal because the regular customers of one driver are sometimes close to the other driver's customers. The planning of the distribution routes itself is done based on each driver's intuition. Nowadays, metaheuristic has been extensively studied for the vehicle routing problem (VRP). This is due to the high complexity of VRP that makes the computation time very long for large problems. The problem in this research is classified as the Heterogeneous VRP with Multiple Trips (HVRPM), following the classification scheme by Caceres-Cruz et al. To the best of the author's knowledge, compared to other variants of VRP, there has not been many studies on metaheuristic to solve the VRPM and very few deals with HVRPM [1]. Several studies that are related to this research are outlined. Paraskevopoulos [3]. The fitness function that is used is a dynamic weighted sum of the total traveling time, penalty to violation of time horizon, and penalty to violation of vehicle capacities. Olivera & Viera develops adaptive memory programming to solve VRPM [4]. The algorithm works by creating a number of initial solutions and inserting attractive sub-routes into an adaptive memory. Then, a new solution is constructed based on that memory and then improved using tabu search. Despaux & Basterrech develops simulated annealing to solve this variant of VRP [5]. The fitness function is the sum of fixed cost, variable cost, and a weighted sum of penalties regarding the constraints. Alonso

A. Development of VNTS Algorithm for HVRPM
The VNS framework is adopted from the study by Cheikh et al. (2015) and the construction of initial solution is based on the study by Paraskevopoulos et al. (2008). The construction procedure is called semi-parallel construction heuristic and was designed for Heterogeneous VRP with Time Windows (HVRPTW).

B. Algorithm Validation
The algorithm is validated by running the algorithm using a small dataset and comparing with the result of exact method (branch-and-bound) using LINGO 18.0. The dataset is from a study by Setiawan et al. (2019) which consists of six customer nodes and three vehicles of two different types. Minor modifications are made to remove the multiple products and add more time components in the time horizon constraint.

C. Parameter Experimentation
There are 5 parameters in the VNTS. The first two are the construction parameters are the weight of components inside the greedy function, which are α 1 and α 2 . The values of α 1 and α 2 are ranged between 0 and 100, with increments of 10%. Next, the parameters of TS will be experimented: MaxTabuIter ∈ {10, 20} and tabu tenure ∈ {10, 20, 30}. Lastly the parameter of VNS is the maximum iterations without any improvement. They are set as MaxVNTSIter ∈ {5, 10, 15}.

D. Algorithm Implementation in Case Study
After the optimal values of the parameters are determined, the VNTS algorithm that has been coded in MATLAB is run using the six instances from the case study, corresponding to the six days of observation. Lastly, the algorithm's performance during its implementation is evaluated.

E. Evaluation of Algorithm Performance
The performance of the developed VNTS algorithm is evaluated based on its average computation time and the solution quality. The dataset used for the evaluation is derived from the instances collected from the case study. The VNTS time is compared with the average time of another metaheuristic, which is Simulated Annealing from Despaux & Basterrech (2014). Next, for the evaluation of the solution quality, the algorithm is compared with exact method because there is lack of standardized dataset for HVRPM. VNTS is compared with branch-and-bound in LINGO.

F. Output Analysis
The output of algorithm implementation is analyzed. Furthermore, the result of evaluation of algorithm performance is discussed.

A. HVRPM Description
The VRP with Multiple Trips (VRPM) is a variant of VRP where each vehicle can serve multiple trips or routes. A trip is defined as a sequence of customer services started at and followed by a visit to the depot and without intermediate stop at the depot. Its characteristics and constraints are similar to the classical CVRP but there is one additional constraint which is concerned with the time horizon. This constraint states that for each vehicle, the total time of all its trips must not exceed a certain time limit. In pure VRPM, the total time only consists of total traveling time [9]. The extension of this variant considers service time at customer and/or loading time at the depot. Meanwhile, the heterogeneous VRP (HVRP) is an extension of VRP where the vehicles are heterogeneous and exist in limited number. The vehicles may differ from one another in terms of capacity, variable cost, and fixed cost. In the case study of this research, the time horizon of each vehicle is also heterogeneous. From these definitions, HVRPM can be defined as a problem where a fleet of heterogeneous vehicles can perform multiple trips within their time horizon and the model seeks to minimize the total cost of these trips.

B. Mathematical Formulation
The mathematical model for HVRPM is developed by modifying the 4-index model of VRPM by Cattaruzza et al. (2016) to suit the problem of the case study.

C. LPG Distribution Case Study (PT. GEP)
PT. GEP is one of the biggest LPG distributors in Malang Raya area that focuses on subsidized LPG. However, it has a subsidiary company that distributes non-subsidized LPG. PT. GEP owns two vehicles that are heterogeneous in capacity, variable cost, and working hour (time horizon). The data for the case study is collected from six days of observation and will be solved individually. The number of customers of the six instances ranges from 13 to 18. The total fuel cost of six days is Rp 823,439.

IV. METHODOLOGY
The metaheuristic method used in this research is based on variable neighborhood search (VNS). Variable neighborhood search (VNS) is a metaheuristic method proposed by Mladenović & Hansen in 1997. It applies a systematic change of neighborhood in two phases: the descent phase to find the local optimum and the perturbation phase to escape from the corresponding valley. The descent phase is performed in a deterministic way while the perturbation is done stochastically. The method of this research is called variable neighborhood tabu search (VNTS) because the descent phase uses tabu search (TS) rather than the commonly used variable neighborhood descent (VND). VNTS was first proposed by Moreno Pérez et al. (2003). The idea behind VNTS is to avoid getting trapped in a local optimum by allowing nonimproving moves. However, this presents the risk of cycling back to previously visited solutions. The use of tabu list in TS helps to minimize this risk.
The VNS framework and neighborhood structures are adopted from the study by Cheikh et al. (2015) while the construction of initial solution and the experimented values of tabu search parameters are based on the study by Paraskevopoulos et al. (2008). The stopping condition of the VNTS is the maximum number of iterations without any improvement, dubbed MaxVNTSIter. The algorithm is then executed using MATLAB software. The procedure of the VNTS algorithm is given in Pseudocode 1 on Figure 1.

A. Construction of Initial Solution
The construction method is called semi-parallel construction heuristic by Paraskevopoulos et al. (2008) and was proposed for HVRP with Time Windows (HVRPTW). The modifications in this research are removing the time windows, adding the checking for the time horizon constraint, and adjusting the mechanism of vehicle removal. The last is done because in HVRPM, vehicles are removed from available vehicles only if the time horizon has been reached.
First, two lists are defined: Cs as the list of unassigned customers and Vk as the list of available vehicles. At each iteration and for each vehicle, the algorithm finds a seed customer based on the furthest feasible customer from the depot. Then the algorithm simultaneously constructs a route for each vehicle using an insertion-based method until either the maximum capacity or time horizon is reached. At this point, the trips are not permanent yet so customers can be assigned to the trips of multiple vehicles. The insertion is done based on a weighted greedy function Φ ijuk that consists of two components (originally six). The first indicates the traveling time increase caused by the insertion while the second forces the algorithm to prioritize customers with large demands in order to maximize the utilization of the vehicle's capacity.
During the above process, the feasibility of the insertion regarding the capacity and time horizon is also checked. After a trip has been constructed for each vehicle, they are compared based on a measure called the Average Cost per Unit Transferred (ACUTk) that is calculated by (13). It indicates the total cost that is incurred in order to carry one unit of customer's demand using vehicle k. The trip with the smallest value of ACUTk is added to the partially constructed solution and the customers are removed from Cs. In the start of the next iteration, if no seed customer can be inserted into a vehicle, then the vehicle is removed from Vk.

B. Neighborhood Structures
There are four neighborhood structures that are used in the following sequence [3]: 1. Customer insertion: This operator removes a customer from a position i and inserts into another position j of different trip (inter-route), of the same trip (intra-route), or of a newly created trip. In other studies, this neighborhood is also referred to as "Relocate" [2], [

C. Shaking Phase
The shaking phase makes use of the first, third, and fourth neighborhood. In every shaking, h neighborhoods are chosen based on a probability distribution and then applied consecutively. Both of these parameters are referenced from Cheikh et al. (2015) with the number of moves, h, being three moves and the probability distribution Pq being P{q 1 , q 3 , q 4 } = {0.6, 0.2, 0.2}.

D. Descent Phase by Tabu Search
The descent phase makes use of the four neighborhood structures which are applied sequentially based on the VNS scheme. The fitness function to evaluate a solution is the objective function of the mathematical model, as shown in (1). The stopping condition of the tabu search is the maximum iterations without any improvement, dubbed MaxTabuIter.

A. Algorithm Validation
Algorithm validation is done using the dataset by Setiawan et al. (2019) for Heterogeneous VRP with Multiple Trips and Multiple Products. Two modifications are made. First, the multiple products into a single product by multiplying the products' quantity by their respective weight and summing the results. Second, service time and loading time are added to the time horizon constraint. Both of these time components are arbitrarily set to 0.5. Finally after the modifications, the dataset is solved using LINGO 18.0 and used as a benchmark for the VNTS algorithm.

B. Parameter Experimentation
Parameter experimentation is done using one of the six instances from the case study. The first two parameters are the construction parameters, α 1 and α 2 . The optimal pair of values for these are 0.6 and 0.4. Here, customer proximity is highly prioritized over capacity utilization. This is because in the case study, the customer demand is relatively high compared to the vehicles' capacity, causing the customer insertion to become inflexible. In fact, a high value of α 2 does not affect the average remaining capacity at all and while failing to do so, the algorithm sacrifices proximity.
The two next parameters are the tabu tenure and MaxTabuIter. The optimal values are 20 and 20. Finally, the selected value of MaxVNTSIter is 15.

C. Algorithm Implementation
The algorithm is run for 10 trials for each instance in the case study. It produces a total cost of Rp 672,562 from the six instances. This saves Rp 150,876 for one operational week, or roughly 18% of the initial cost as Table 1 shows.

D. Evaluation of Algorithm Performance
The algorithm performance during its implementation is evaluated in terms of the computation time and solution quality. Regarding the time, VNTS solves the instances of the case study with an average time of 8 to 23.7 seconds. In addition, the average time is also compared with another metaheuristic, which is Simulated Annealing (SA) to solve HVRPM with Time Windows (HVRPMTW) by Despaux & Basterrech (2014). A dummy dataset is made from the instances of the case study to match the problem size of the study. With problems of the same size, SA solves HVRPMTW with an average time of 35 seconds. On the other hand, the developed VNTS algorithm solves HVRPM within 28.8 seconds in average.
Next, the solution quality of the VNTS algorithm is evaluated by comparing it with B&B method in LINGO. This is done twice for each dataset. First, the data from the case study is solved with specific time limits. To make the evaluation fair, the computation time limit for each dataset is set as the total computation time of 10 trials of VNTS algorithm, following Table 2. Second, the data is solved again using LINGO but with the maximum computation time being 2 hours. The purpose of second run is to obtain the best possible result within reasonable time, also in hopes of getting global optimal results.

A. Analysis of Algorithm Implementation
The algorithm is implemented using the best parameter settings that have been identified, which are α 1 = 0.6, α 2 = 0.4, tabu tenure = 20, MaxTabuIter = 20, and MaxVNTSIter = 15. The algorithm is run for 10 trials and the generated trips produce a total cost of Rp 672,562. If the new trips are implemented, the company will be able to save Rp 150,876 or 18% of its initial cost as Table 3 and Table 4 shows.
There are two factors that are improved in the new generated trips. The first factor is the re-assignment customers to the vehicles. In the new trips, it can be seen that the customers assigned to each vehicle differs from the initial trips. This indicates that the fixed mutual agreement between the drivers regarding their regular customers is not optimal. This improvement is facilitated by the inter-route moves from the first and third neighborhood as well as the second and fourth neighborhood structure. The first improvement causes more customers to be assigned to vehicle 2. This leads to a reduction in the total cost because vehicle 2 has a lower fuel cost than vehicle 1. However, there is also a tradeoff where its capacity is lower, meaning that it cannot bring as much load in one trip. With this tradeoff, distant customers with small demand should be assigned to vehicle 2 so that the vehicle is not overloaded and the high distance is negated through the vehicle's low fuel cost.
Besides reducing the total cost, the re-assignment to vehicle 2 also causes the work load between the two vehicles to become more balanced. Initially, the ratio of working time (the sum of loading/unloading time, traveling time, and service time in one day) in one week is roughly 27:10 between vehicles 1 and 2. In the new trips, this ratio decreases to 65% to 35% or approximately 20:10. It is to be noted, however, that the scope of this ratio only includes the distribution of subsidized LPG. For vehicle 2, half of its working hour is allocated to distributing non-subsidized LPG by customer order as well. Therefore the workload of vehicle 2 will increase again to a certain extent as Table 5 and Table  6 shows.
The second improvement in the new trips is the routing of customers within a vehicle. With the same customers as the initial trip, the algorithm produces a better result by rearranging which customers are served in which trip. For instance, on day 6, the customers of vehicle 1 are the same with that from the initial trips but the arrangement of the customers differs greatly.
In the initial trips, the capacity is not utilized efficiently so the remaining capacity ranges from 30 to 50 units. Moreover, the vehicle needs to perform one last trip only to serve one customer, leaving 220 units of remaining capacity. In contrast, the new trips are able to maximize the capacity utilization, leaving only 10 to 45 units of remaining capacity. As a result, the travel time of vehicle 1 decreases and the number of trips performed by the vehicle also decreases. With this improvement alone, a savings of Rp 41,090 is obtained on day 6.

B. Analysis of Algorithm Performance
The first aspect of the algorithm performance is the average computation time. For the six instances from the case study, the average time ranges from 8 to 23.7 seconds which the author considers reasonable. Next, VNTS is compared with Simulated Annealing for HVRPMTW by Despaux & Basterrech (2014). Admittedly, the referenced study is more complex with the additional time window constraint. Nevertheless, it is used as a benchmark because it is the study with most similar variant of VRP and it tests the algorithm in a relatively small dataset. Based on the comparison, VNTS solves a problem of the same size in a fairly lower computation time.
Next, the solution quality of the VNTS algorithm is analyzed by comparing it with the result of exact method. In the first comparison, the case study is solved in LINGO and the maximum computation time is set equal to the total   Table 3, the GAP values are calculated. In this comparison, VNTS produces results that are better or equal to B&B. On day 1 and day 3 to 6, the VNTS generates better solutions than B&B, indicated by the GAP value being negative. For these instances, the GAP ranges from -0.0001 to -0.0395 as Table 7 shows. In the form of cost, this difference ranges from Rp 12 to Rp 4,843. On day 2, the best cost generated by VNTS is equal to that by B&B because the GAP is zero. Therefore in this comparison, the VNTS performs better than B&B.
The next comparison is done by running the LINGO software for two hours. This is done in hopes of getting a global or near-global optimal solution using exact method so that the solution quality of VNTS can further be evaluated. From the six instances, only the costs of day 2, 3, 4, and 6 are improved so the GAP values increase only for these days.
For the dataset of day 2, the best cost by VNTS is higher than the exact method so the GAP value increases from 0 to 0.0007. In this case, the cost difference is negligible since it is only Rp 72. On day 3 and day 4, the VNTS solutions still have lower cost than B&B, with a difference of Rp 2,643 and Rp 3,755 respectively as Table 8 shows. The GAP value of day 2 is far smaller than the absolute value of the GAP on day 3-4. So at this point, VNTS outperforms B&B. Lastly for day 6, the solution of the exact method is globally optimal and the corresponding GAP value is zero. The developed VNTS algorithm succeeds to obtain a global optimal solution within a notably lower computation time, approximately 67 times faster.

VII. CONCLUSION
In this research, a Variable Neighborhood Tabu Search (VNTS) algorithm to solve the Heterogeneous VRP with Multiple Trips (HVRPM) has been developed. The developed algorithm is validated by comparing the result with the global optimal result from LINGO software. Based on the results of three trials, the algorithm produces the same solution as the global optimal solution and the algorithm is considered valid.
Then it is run for 10 trials using the six instances from the case study of PT. GEP. The trips generated by the algorithm is able to reduce the total fuel cost of one operational week by Rp 150,876 or 18% of the initial cost. The performance of the algorithm is evaluated using the instances from the case study. The computation time is reasonable and is competitive when compared with Simulated Annealing. The solution quality is compared with branch-and-bound method. VNTS achieves one global optimal solution out of six instances and overall, the quality of the solutions are better.