A New Optimization Technique for the Location and Routing Management in Agricultural Logistics

: This paper aims to solve the location and routing problem (LRP) in the agricultural sector with the objective function of fuel cost minimization. Many farmers may have problems when transporting and selling products because of high costs and unfair prices. The proper location of standardized procurement centers and suitable routes will relieve farmers’ problems. This paper includes a realistic constraint that a farm can be visited to collect product more than once. A mathematical model was formulated to be solved by Lingo software, but when the problem size was larger, Lingo was unable to solve the problem within a reasonable processing time. The variable neighborhood strategy adaptive search (VaNSAS) is proposed to solve this LRP. The main contributions of this paper are a real case study problem and the ﬁrst introduction of VaNSAS. Furthermore, the di ﬀ erent combinations of the solution approach are proposed to prove which combination is the best algorithm. The computational results show that VaNSAS can ﬁnd the solutions for all problem sizes in much less processing time compared to Lingo. In medium and large-sized instances, the VaNSAS can reduce processing times by 99.91% and 99.86%, respectively, from solutions obtained by Lingo. Finally, the proposed VaNSAS has been deployed in a case study problem to decide the best locations and transportation routes with the lowest fuel cost.


Introduction
Open innovation plays an important role in driving the achievement of an organization, while integrating inside knowledge with new outside technologies. An organization driving a policy of open innovation will have great opportunities to receive new innovation that can be applied to improve the operating performance of the firm. Currently, we are in the era of the fourth industrial revolution (4IR) and the fourth agricultural revolution (4AR). The 4IR technologies-such as in the areas of robotics, artificial intelligence, and optimization software-will be applied to the field of agriculture. These 4IR technologies will have a significant impact for solving agricultural problems. The 4AR is the period of technological improvement and increase crop productivity. The target is to maximize crop yield and minimize productive factors [1,2]. Due to the fourth agricultural revolution, Road transportation is essential for agricultural sectors in developing countries, but it also has a high transportation cost for farmers due to high fuel prices. When rubber is transported, fuel costs directly impact the total cost, especially when traveling long distances. The total cost for harvesting and transporting natural rubber is 11.55 Baht/Kg. This cost can be divided across two stakeholders, such that a farmer pays 10.53 Baht/Kg and the middlemen pay 1.02 Baht/Kg (assuming passing through two middlemen). The transportation cost for farmer is 1. 29 Baht/Kg from the total cost 10.53 Baht/Kg, and only 0.23 Baht/Kg for the middlemen [8]. Considering transportation in rubber logistics, direct shipment from farmer to middleman is still a weak point that occurs in the collection process, and results in high transportation costs when there are a greater number of routes and vehicles used. Furthermore, there should be an appropriate number of collecting center locations to accommodate the effective transport of rubber. For example, if the number of collecting centers is lower than the demand, farmers have to send rubber a longer distance, which leads to higher fuel costs. Therefore, it is the challenge for the Thai rubber industry to improve logistics operations and reduce transportation costs for a competitive and profitable future.
Thus, this research determines locations for rubber procurement centers and routes of the transportation from rubber farmers to established procurement centers. The problem considers a combination of location planning and routing decisions. Akararungruangkul and Kaewman (2018) [9] have formerly studied the LRP in the agricultural sector. The paper presents a special case of LRP in which the objective function was the amount of fuel used. Road categories were divided into seven types to consider different fuels used in each road type; however, a mathematical model was not formulated for the problem. Subsequently, the adaptive large neighborhood search (ALNS) was applied as the solution approach. Then, the computation performance was compared with the current practice. In this paper, a mathematical model is formulated for the objective of fuel consumption cost minimization, from the perspective of collecting problems. The model was tested by the Lingo program, which is the exact method software. A special constraint was added by allowing for rubber farms that may be served more than once, in cases where they have more rubber than the vehicle capacity. This provides the problem with practical perspective. Moreover, in this research, a new solution method based on a metaheuristic called the variable neighborhood strategy adaptive search (VaNSAS) will be first introduced. This is a new idea of metaheuristics that lets the solution search for the better solution in different areas by applying various local search strategies. The main contributions of this paper are two-fold. First, the formulated problem is from a real case study in the southern region of Thailand. The impact of road types on fuel consumption cost is considered. Alongside this, the paper considers that the rubber quantity may be more than the vehicle capacity, so the rubber farm can be visited more than once. Secondly, this paper presents a new concept of algorithms. The VaNSAS does not appear in any paper published in the literature, which is a method for solving combinatorial problems.
This paper consists of six sections: Section 2 presents the literature review; Section 3 gives an explanation of the problem statement and mathematical model formulation; Section 4 introduces the VaNSAS algorithms; Section 5 reveals the computational studies; and Section 6 is the conclusions of this paper.

Literature Review
Supply chain activities are around 10% of gross domestic product [10]. Supply chain activities have a great impact on the costing and profit of a firm. It has been a long time since the location routing problem (LRP) has played a crucial role in strategic planning to successfully compete in the business world, by increasing the performance of interconnections effectively and reducing risk in each area. LRP is combined of two kinds of problems: (1) where to open facilities; (2) how to find a route for products from facilities to customers. LRP is widely applicable to many sectors such as food and drink delivery [11], newspapers distribution [12], hazardous waste collection [13], electric vehicle (EV) charging stations [14], and military operations [15]. There are many factors that influence location selection and routing decisions, including the satisfaction of customers who receive or deliver goods according to the relationship of the business network, and to help achieve more success.
Road cargo transportation is necessary for the local and regional levels, but it has negative impacts on the environment and global climate change. As road transportation activity is a critical source of energy consumption and gas emissions-which are N 2 O, CO 2 , and PM (particulate matter)-green logistics has gained close attention from governments and organizations [16,17]. The classical LRP aims to provide the proper locations and routing of the vehicle to serve all customers. The main goal of classical LRP is total cost minimization, which comprises of the cost of open depots, the vehicles fixed cost, and the transportation cost, which have all been extensively studied [18][19][20][21][22][23]. The GLRP (green LRP) is an extension of classical LRP which aims to minimize environmental impacts by addressing factors such as fuel consumption amounts and carbon emission amounts [16,[24][25][26], and carbon emission costs [27]. Nevertheless, these papers formulated the problem with respect to the distribution problem. From the best of the authors' review, there were few LRP papers that studied the collecting problem and minimized fuel consumption cost in different road categories.
LRP was considered as NP hard as it integrates the capacitated facility location problem with the multi-depot vehicle routing problem [28]. Many heuristic and metaheuristic methods were proposed to solve this problem. In the case of heuristics,   [29] proposed four algorithms of heuristics to solve multi-period multi-echelon inventory transportation problem. The solution approach tested with 15 instances. The results were beneficial for related manufacturing organizations.   [30] proposed a novel approximate algorithm for solving fuel bunker management problem for liner shipping networks. Barreto et al. (2007) [31] proposed various clustering techniques which were contained in a sequential heuristic for solving the LRP. The problem was the total cost minimization and a customer was visited just once. Clustering techniques were also applied by Boudahri et al. (2013) [32] to solve LRP for broiler products. In addition, metaheuristics have been presented for the LRP due to fast and simple features; they provide feasible results within reasonable processing times. Examples of the recognized metaheuristics are simulated annealing (SA) [21,[33][34][35][36], iterated local search (ILS) [37], variable neighborhood search (VNS) [38,39], and adaptive large neighborhood search (ALNS).
VNS was first proposed in 1997 by Mladenovic and Hansen (1997) [40], who introduced the systematic move of the searching area to escape the local optima trap. From a recent survey, VNS has been widely developed both in methodologies and applications. Jarboui et al. (2013) [41] integrated five neighborhood structures as the local search area in the scheme of variable neighborhoods to solve LRP with multiple depots. The results were successfully compared with other methods in the previous study. Schwengerer et al. (2012) [42] suggested VNS to solve the two-echelon location-routing problem (2ELRP). They proposed seven neighborhood structures and included a variable perturbation range. The computational results showed that the solution's quality and runtime were competitive with regard to other methods in the literature. Macedo et al. (2013) [43] studied LRP with capacitated vehicles and depots, and multiple routes (MLRP) is a variant of traditional LRP. They proposed VNS with six neighborhood structures which were used for routing, location, and multiple routes. The results indicated that the method was effective in terms of good results and short processing times. From the brief review above, VNS is shown to be successful for solving the LRP problem.
ALNS was first introduced by Ropke and Pisinger (2006) [44] and developed the large neighborhood search heuristic of Shaw (1997) [45]. The basic idea is to improve the solution at each iteration by deploying the destroy and repair technique to the current solution. Normally, a number of destroy and repair methods are unrestricted, and it is randomly selected at each iteration. Each operator has a weight that is related to the selection probability of an operator, this weight is adaptive and constantly changes based on its past performances. Hemmelmayr et al. (2012) [46] applied ALNS to the two-echelon vehicle routing problem (2E-VRP). This paper was focused on the VRP problem, but the authors also applied their ALNS algorithm to LRP instances. They developed eight new destroy operators and five repair operators. The computational experiments were tested in some instances in the literature; the result revealed that ALNS achieved outstanding results on the LRP. Contardo et al. (2012) [47] proposed ALNS for solving the two-echelon capacitated location-routing problem (2E-CLRP) by applying eight destroy operators and four repair operators. The results showed that ALNS could improve 133 instances out of 147 of the best-known solutions.   [48] presented ALNS as the solution approach for the LRP with intra-route facilities to establish the facilities' locations, whereas the depot's location is predefined. The ALNS can find the new best-known solutions for the benchmark instances from the literature. Moreover, according to Prodhon and Prins (2014) [49], ALNS is the best algorithm when testing on the third best on Barreto and Prodhon sets and Tuzun's instances.
Reviewing literature indicated the effectiveness of metaheuristics methods, especially VNS and ALNS algorithms. Therefore, this paper combines both methods into a new VaNSAS algorithm for solving our LRP problem.

Problem Statement
This model aims to find the best location among all candidate rubber farms, then establish this location as the rubber procurement center. After that, the best transportation routes will be arranged for collecting rubber from farm to procurement center. The rubber quantity from each farm, capacity of the vehicle, and capacity of the procurement center are predefined parameters. The transportation road is different according to the road type, as indicated by different arrow sizes. The rubber farm is allowed to be served more than once by directed shipment (partial shipment), in case rubber quantity is over vehicle capacity, as indicated by the dotted lines. The problem statement is described in Figure 1. to Prodhon and Prins (2014) [49], ALNS is the best algorithm when testing on the third best on Barreto and Prodhon sets and Tuzun's instances.
Reviewing literature indicated the effectiveness of metaheuristics methods, especially VNS and ALNS algorithms. Therefore, this paper combines both methods into a new VaNSAS algorithm for solving our LRP problem.

Problem Statement
This model aims to find the best location among all candidate rubber farms, then establish this location as the rubber procurement center. After that, the best transportation routes will be arranged for collecting rubber from farm to procurement center. The rubber quantity from each farm, capacity of the vehicle, and capacity of the procurement center are predefined parameters. The transportation road is different according to the road type, as indicated by different arrow sizes. The rubber farm is allowed to be served more than once by directed shipment (partial shipment), in case rubber quantity is over vehicle capacity, as indicated by the dotted lines. The problem statement is described in Figure  1. This paper focuses on a case study of rubber logistics in southern Thailand. The case study contains 95 rubber farms. The road categories have been adopted from Akararungruangkul and Kaewman, 2018) [9]; divided into seven types based on average driving speed. The road type and fuel consumption rate are shown in Table 1.  This paper focuses on a case study of rubber logistics in southern Thailand. The case study contains 95 rubber farms. The road categories have been adopted from Akararungruangkul and Kaewman, 2018) [9]; divided into seven types based on average driving speed. The road type and fuel consumption rate are shown in Table 1.  Tables 2 and 3 illustrate the example of road type metric and distance metric for five rubber farms. Then, the fuel consumption metric can be calculated as shown in Table 4. The fuel consumption is calculated by multiplied fuel consumption rate and distance. This metric distinguishes our problem from the classical LRP which aims to optimize the travelling distance, as our model aims to minimize fuel consumption. Let us consider the following routes; 1-4-2-5-3-1 and 1-3-4-2-5-1. These routes have the total distance equal to 177 km and 179 km, but the fuel consumptions are 19.487 L and 18.978 L, respectively. Although route 1-3-4-2-5-1 is longer, our model will select this route for the reason of fuel minimization.

Mathematical Model
This study proposed a model that makes the decision of location routing design while considering environmental impact. The proposed model simultaneously determines the rubber procurement center locations and transportation routes with the objective function to minimize the fuel consumption costs of rubber collection. This combined approach manages location and vehicle routing plans at the same time. It is referred to in the literature surveys on LRPs [28,50]. The details of objective function, constraints, indices, parameters, decision variables, and support variables are shown below Distance from farm i to procurement center j (km) F ij Fuel consumption rate from farm i to procurement center j V Capacity of vehicle (kg) P j Capacity of procurement center (kg) Decision Variables z ij = 1 if farm i is assigned to procurement center j and direct shipment is appeared = 0 otherwise x ij = 1 if travel from farm i to farm j and routing for the rest of direct shipment is appeared = 0 otherwise n i = Number of direct shipments at farm i Support Variables u i Accumulated rubber quantity in vehicle at farm i for sub-tour elimination m j Number of round travels of procurement center j r i Remaining rubber after direct shipment from farm i y j = 1 if procurement center j is chosen = 0 otherwise Objective function The mathematical model shown above can be described as following: (1) Is the objective function which is to minimize the total fuel consumption cost. The first term is direct shipment fuel cost, the second is routing fuel cost. Constraint (2) states that each farm must be assigned to only one procurement center. Constraint (3) determines the number of direct shipments and remaining rubber in case rubber quantity at the farm is over the capacity of vehicle. Constraint (4) ensures that total rubber at the procurement center must not exceed capacity. Constraint (5) determines the transportation route of rubber collection from the farm. Constraint (6) ensures that whenever the vehicle enters farm i, it must leave farm i. Constraint (7) provides the transportation route among the farms which are assigned to the same procurement center. Constraints (8) and (9) are capacity constraints related to vehicle capacity and sub-tour elimination. Finally, constraint (10) is a binary variable.

Solution Approach
This paper introduces the new metaheuristic-based method called variable neighborhood strategy adaptive search (VaNSAS), which is the new concept for solving combinatorial optimization problems. The elemental idea of the VaNSAS is to let the solution search for the better solution in a different area by using several searching methods. The searching method in VaNSAS can be a basic local search, well-known heuristics, modified metaheuristics, and a new designed local search, which is an advantage of VaNSAS that lets the designer put new ideas into the algorithm. Therefore, VaNSAS is very flexible to use for many optimization problems. Generally, the VaNSAS is composed of five steps which are: (1) generate a set of tracks; (2) all tracks select the specified black box; (3) execute the selected black box; (4) improve the track; and (5) repeat steps (1) to (4). The VaNSAS procedure can be explained as follows.

Generate a Set of Tracks
A track is the representative for the solution. A set of tracks is initially generated in terms of the real number. For example, if five tracks are generated and each track contains 15 elements, a set of tracks will be obtained as shown in Table 5. The number of tracks (NT) in each iteration is equal to five in this case-it is a predefined parameter that depends on designer decisions and problem size. Number of elements (NE) can be determined by equation 11, where NF is number of farm and NP is number of procurement centers. If the problem has five farms and two procurement center, then NE = 2(5)+2 = 12 elements.

Track Interpretation Method
Track interpretation is the procedure that changes the code in each track to the solution of the problem. In order to illustrate the interpretation method, the related information is provided in Table 6. Let us consider five rubber farms with aggregate rubber of 38,270 kg, with the procurement centers chosen holding a capacity of 25,000 kg each. The number of procurement centers (NP) is set by aggregate rubber divided by 80% of their capacity, so in this problem, NP is equal to 2. The track interpretation method can be explained by the following steps.
Step 1: Divide elements into three groups: A (depot or procurement center) = 5 elements, B (farm) = 5 elements and C (sequence) = 2 elements, as shown in Table 7. A sequence is used for setting the priority of a depot when assigning a farm to a depot. Step 2: Sort the element value in each group by ascending order, as shown in Table 8. Step 3: Group A: select the first two numbers (NP = 2) to be the procurement center, so farm 5 and farm 2 will be depots.
Group C: The sequence is 2 and 1; this means depot 2 is the first priority to assign, so start assigning a farm to depot 2, then to depot 5.
Step 4: Group B: start assigning farm 3 and farm 4 to depot 2 until it satisfies the capacity constraint. The total rubber is farm 2 (depot) + farm 3 + farm 4 = 19,160 kg, but the vehicle load is farm 3 + farm 4 = 12,640 kg. Then, the assignment will be stopped and depot 5 will be started. The next number is farm 1, so assign it to depot 5. The partial shipment occurs in farm 1 because it has more rubber than vehicle capacity. The next number is farm 2, but it is skipped because it has been selected to be a depot already.
This step will allow the initial solution to always obtain feasible state because the assignment method will satisfy the constraint of vehicle capacity and depot capacity at all iterations.
From Figure 2, farm 1 will be visited twice according to our special constraint. The objective function will be calculated by using fuel consumption data given in Table 4. For route 2-3-4-2, the fuel consumption cost is (3.136 + 4.508 + 3.920) × 28 = 323.792 baht. For route 5-1-5-1-5, (2.430 × 2 × 2) × 28 = 272.160 baht. Thus, the total fuel consumption of this track is 323.792 + 272.160 = 595.952 baht. Group C: The sequence is 2 and 1; this means depot 2 is the first priority to assign, so start assigning a farm to depot 2, then to depot 5.
Step 4: Group B: start assigning farm 3 and farm 4 to depot 2 until it satisfies the capacity constraint. The total rubber is farm 2 (depot) + farm 3 + farm 4 = 19,160 kg, but the vehicle load is farm 3 + farm 4 = 12,640 kg. Then, the assignment will be stopped and depot 5 will be started. The next number is farm 1, so assign it to depot 5. The partial shipment occurs in farm 1 because it has more rubber than vehicle capacity. The next number is farm 2, but it is skipped because it has been selected to be a depot already.
This step will allow the initial solution to always obtain feasible state because the assignment method will satisfy the constraint of vehicle capacity and depot capacity at all iterations.  Figure 2, farm 1 will be visited twice according to our special constraint. The objective function will be calculated by using fuel consumption data given in Table 4. For route 2-3-4-2, the fuel consumption cost is (3.136 + 4.508 + 3.920) × 28 = 323.792 baht. For route 5-1-5-1-5, (2.430 × 2 × 2) × 28 = 272.160 baht. Thus, the total fuel consumption of this track is 323.792 + 272.160 = 595.952 baht.

All Tracks Select the Specified Black Box
Black box is the algorithm that is designed to improve the solution. This step allows the tracks obtained from previous step to improve themselves by selecting to enter in a black box. The black box methods and the black box numbers are unrestricted. It can be a simple local search or wellknown heuristics or any metaheuristics, such as genetic algorithm (GA), simulated annealing (SA), and particle swarm optimization (PSO). It is based on the designer's decision and the complexity of the problem. In each iteration, a track will randomly select the preferred black box. This step means that a single solution can change the neighborhood at any the time. The selection of a black box will be performed using the probability in Equations (12)

All Tracks Select the Specified Black Box
Black box is the algorithm that is designed to improve the solution. This step allows the tracks obtained from previous step to improve themselves by selecting to enter in a black box. The black box methods and the black box numbers are unrestricted. It can be a simple local search or well-known heuristics or any metaheuristics, such as genetic algorithm (GA), simulated annealing (SA), and particle swarm optimization (PSO). It is based on the designer's decision and the complexity of the problem. In each iteration, a track will randomly select the preferred black box. This step means that a single solution can change the neighborhood at any the time. The selection of a black box will be performed using the probability in Equations (12)- (13).
where P bt is probability of selecting black box b in iteration t. I bt-1 is a bonus point; it is given 1 point if the black box of the last iteration contains the best solution, but given zero if this is not the case. A bt-1 is the average objective value of every track which selects black box b in the last iteration. N bt-1 is the number of tracks that select black box b in the last t iteration. Z bt-1 is the best objective obtained from black box b in the last iteration. F, K, and M are predefined parameters.

Execute the Selected Black Box
This paper proposes that three black boxes contain three algorithms which are (1) SWAP heuristic; (2) adaptive large neighborhood search (ALNS); and (3) variable neighborhood search (VNS). The procedure of these algorithms will be revealed as follows. . The procedure of SWAP shows as the following steps.
Step 1: Initial solution obtained from the track that chooses to operate in this black box.
Step 2: Randomly select two indices (a, b) = (farm 2, farm 5), then exchange their position. Example of SWAP is illustrated in Figure 3. . The procedure of SWAP shows as the following steps.
Step 1: Initial solution obtained from the track that chooses to operate in this black box.
Step 2: Randomly select two indices (a, b) = (farm 2, farm 5), then exchange their position. Example of SWAP is illustrated in Figure 3.
Step 3: Perform interpretation process of a track.
Step 4: Repeat step 2 to 4 until all positions have been exchanged.
Step 5: Return the best solution. ALNS consists of a destroying and a repairing operation; each operation contains a set of destroying and repairing methods called operators. The operator will be assigned a weight which is adjusted in every iteration based on its success. The weight reflects the operator's performance and the selection of an operator in each iteration is based on its weight. The operator which has a high performance will have more chance to be selected repeatedly. The procedure of ALNS shows as the following steps.
Step 1: Initial solution obtained from the track that chooses to operate in this black box; let it be s.
Step 2: Setting the initial solution s is the best solution s*.
Step 3: Setting the weight and selected probability of each operator.
Step 4: Randomly select destroying and repairing operators according to their weight.
Step 5: Apply the operators to the current solution, then a new solution (s') will be obtained.
Step 6: Perform interpretation process of a track.
Step 7: Perform acceptance criterion. If the acceptance criterion is met, then s  s'. If s is better than s*, then s*  s'.
Step 9: Repeat step 3 to 8 until stopping criterion is met.
Step 10: Return the best solution. The acceptance criterion in Step 7 is performed for accepting or rejecting a new solution (s'). This paper uses a greedy acceptance principle which only accept s', if it is better than the current solution Step 3: Perform interpretation process of a track.
Step 4: Repeat step 2 to 4 until all positions have been exchanged.
Step 5: Return the best solution.

Adaptive Large Neighborhood Search (ALNS)
ALNS consists of a destroying and a repairing operation; each operation contains a set of destroying and repairing methods called operators. The operator will be assigned a weight which is adjusted in every iteration based on its success. The weight reflects the operator's performance and the selection of an operator in each iteration is based on its weight. The operator which has a high performance will have more chance to be selected repeatedly. The procedure of ALNS shows as the following steps.
Step 1: Initial solution obtained from the track that chooses to operate in this black box; let it be s.
Step 2: Setting the initial solution s is the best solution s*.
Step 3: Setting the weight and selected probability of each operator.
Step 4: Randomly select destroying and repairing operators according to their weight.
Step 5: Apply the operators to the current solution, then a new solution (s') will be obtained.
Step 6: Perform interpretation process of a track.
Step 7: Perform acceptance criterion. If the acceptance criterion is met, then s ← s'.
If s is better than s*, then s* ← s'.
Step 9: Repeat step 3 to 8 until stopping criterion is met.
Step 10: Return the best solution. The acceptance criterion in Step 7 is performed for accepting or rejecting a new solution (s'). This paper uses a greedy acceptance principle which only accept s', if it is better than the current solution (s).
This paper designs three destroying operators and two repairing operators which will be explained in the below section.

Destroying Operators
Destroying operators are used to remove a part of the current solution. Degree of removal (d) is randomly selected in each iteration by varying amounts, from 10% to 40% of farms, representing light removal to strong removal. After selection of d, the q farms will then be removed from the current solution. An example of the destroying operation is shown in Figure 4. Destroying operators are used to remove a part of the current solution. Degree of removal (d) is randomly selected in each iteration by varying amounts, from 10% to 40% of farms, representing light removal to strong removal. After selection of d, the q farms will then be removed from the current solution. An example of the destroying operation is shown in Figure 4. •

Random Removal
This operator randomly removes q farms from the initial solution.
Step 1: Randomly choose degree of removal, let d = 30%. If the initial solution contains seven farms, two farms will be removed.
Step 2: Remove those farms to the farm pool. •

Worst Removal
This operator removes the worst farm that generates the highest cost in the initial solution.
Step 1: Randomly choose degree of removal, let d = 30%. If the initial solution contains seven farms, two farms will be removed.
Step 2: Calculate the total cost of the initial solution.
Step 3: Remove farm one by one from the initial solution, then calculate the different value by comparing it to the total cost of the initial solution.
Step 4: Ranking farm by decreasing order based on different values in Step 3.
Step 5: Remove the first two farms which have large values to the farm pool. •

Related Removal
Step 1: Randomly choose degree of removal, let d = 30%. If the initial solution contains seven farms, two farms will be removed.
Step 2: Randomly choose a farm, called a seed farm.
Step 3: Calculate the value of total cost when traveling from the seed farm to each farm.
Step 3: Ranking farms by descending order based on values in Step 3.
Step 4: Remove q-1 farms which have the least value to the farm pool.

Repairing Operators
When a destroying operator is executed, then a repairing operator is applied to re-insert the farms in different positions by expecting the better solution. An example of repairing operation is shown in Figure 5.  •

Random Removal
This operator randomly removes q farms from the initial solution.
Step 1: Randomly choose degree of removal, let d = 30%. If the initial solution contains seven farms, two farms will be removed.
Step 2: Remove those farms to the farm pool. •

Worst Removal
This operator removes the worst farm that generates the highest cost in the initial solution.
Step 1: Randomly choose degree of removal, let d = 30%. If the initial solution contains seven farms, two farms will be removed.
Step 2: Calculate the total cost of the initial solution.
Step 3: Remove farm one by one from the initial solution, then calculate the different value by comparing it to the total cost of the initial solution.
Step 4: Ranking farm by decreasing order based on different values in Step 3.
Step 5: Remove the first two farms which have large values to the farm pool. •

Related Removal
Step 1: Randomly choose degree of removal, let d = 30%. If the initial solution contains seven farms, two farms will be removed.
Step 2: Randomly choose a farm, called a seed farm.
Step 3: Calculate the value of total cost when traveling from the seed farm to each farm.
Step 3: Ranking farms by descending order based on values in Step 3.
Step 4: Remove q-1 farms which have the least value to the farm pool.

Repairing Operators
When a destroying operator is executed, then a repairing operator is applied to re-insert the farms in different positions by expecting the better solution. An example of repairing operation is shown in Figure 5.
Step 1: Randomly choose degree of removal, let d = 30%. If the initial solution contains seven farms, two farms will be removed.
Step 2: Randomly choose a farm, called a seed farm.
Step 3: Calculate the value of total cost when traveling from the seed farm to each farm.
Step 3: Ranking farms by descending order based on values in Step 3.
Step 4: Remove q-1 farms which have the least value to the farm pool.

Repairing Operators
When a destroying operator is executed, then a repairing operator is applied to re-insert the farms in different positions by expecting the better solution. An example of repairing operation is shown in Figure 5. •

Random Insertion
The farms from the farm pool will be randomly inserted back to the solution. •

Best Insertion
The algorithm will search the best position by considering the minimal total cost.
Step 1: Determine a position to insert a farm, then calculate the total cost.
Step 2: Insert a farm into the lowest total cost position.
Step 3: Re-do step 1 and step 2 until all farms in the farm pool are inserted.

• Farms Exchange (N 1 )
A neighbor of a solution S of N 1 is acquired by exchanging two random farms from the same route or same depot, from the different route or a different depot. Example of a farm's exchange structure is shown in Figure 6. The farms from the farm pool will be randomly inserted back to the solution. •

Best Insertion
The algorithm will search the best position by considering the minimal total cost.
Step 1: Determine a position to insert a farm, then calculate the total cost.
Step 2: Insert a farm into the lowest total cost position.
Step 3: Re-do step 1 and step 2 until all farms in the farm pool are inserted.

Farms Exchange (N1)
A neighbor of a solution S of N1 is acquired by exchanging two random farms from the same route or same depot, from the different route or a different depot. Example of a farm's exchange structure is shown in Figure 6. •

Depots Exchange (N2)
A neighbor of a solution S of N2 is acquired by exchanging the farms between two random depots. All farms in a depot will be re-assigned to another depot by exchanging each other. An example of a depot exchange structure is shown in Figure 7.

• Depots Exchange (N 2 )
A neighbor of a solution S of N 2 is acquired by exchanging the farms between two random depots. All farms in a depot will be re-assigned to another depot by exchanging each other. An example of a depot exchange structure is shown in Figure 7. The farms from the farm pool will be randomly inserted back to the solution. •

Best Insertion
The algorithm will search the best position by considering the minimal total cost.
Step 1: Determine a position to insert a farm, then calculate the total cost.
Step 2: Insert a farm into the lowest total cost position.
Step 3: Re-do step 1 and step 2 until all farms in the farm pool are inserted.

Farms Exchange (N1)
A neighbor of a solution S of N1 is acquired by exchanging two random farms from the same route or same depot, from the different route or a different depot. Example of a farm's exchange structure is shown in Figure 6. •

Depots Exchange (N2)
A neighbor of a solution S of N2 is acquired by exchanging the farms between two random depots. All farms in a depot will be re-assigned to another depot by exchanging each other. An example of a depot exchange structure is shown in Figure 7.

• Depot Status Change (N 3 )
A neighbor of a solution S of N 3 is acquired by closing an open depot and opening another one, then re-assigning all farms to the new depot. An example of a depot status change structure is shown in Figure 8.
where is track i, element j at iteration t + 1. is track i, element j at the last iteration. and are predefined parameters.
is the first randomly selected track and is the second randomly selected track.
is the personal best track and is the global best track.

Repeat
Step 1 to 4 Repeat steps 2 to 4 iteratively. VaNSAS is the new type of metaheuristics for which the target is to let algorithms search in many different areas to obtain the best solution possible. The searching can move to have more diversification and intensification at all times depending on the black box methods that are designed to use. The VaNSAS algorithm is shown in Algorithm 2.
where Z ijt+1 is track i, element j at iteration t + 1. Z ijt is track i, element j at the last iteration. α and β are predefined parameters. Z 2 jt is the first randomly selected track and Z 3 jt is the second randomly selected track. Z pb ijt is the personal best track and Z gb ijt is the global best track.

Repeat
Step 1 to 4 Repeat steps 2 to 4 iteratively. VaNSAS is the new type of metaheuristics for which the target is to let algorithms search in many different areas to obtain the best solution possible. The searching can move to have more diversification and intensification at all times depending on the black box methods that are designed to use. The VaNSAS algorithm is shown in Algorithm 2.

Algorithm 2 VaNSAS
Input Number of farms, fuel consumption rate, distance, rubber quantity, related capacity Output Fuel cost Begin While I less than predefined number of iterations.

Computational Framework and Result
The proposed heuristic has been tested with three groups of testing instances. Group 1: small-sized cases, contains five to nine rubber farms. Group 2: medium-sized cases, contains 10-30 rubber farms. Group 3: medium-sized cases, contains 35-50 rubber farms and the case study problem. The computational framework is shown in Table 9. The complexity of the problem instances is explained as follows. The number of constraints and decision variables of Group 1 are 255-1,107 and 55-171 respectively. The number of constraints and decision variables of Group 2 are 9,720-30,780 and 820-1,830 respectively. The number of constraints and decision variables of Group 3 are 47,985-894,045 and 2,485-18,145 respectively. All proposed algorithms have been compared with the solutions obtained from Lingo v.11 software. The lower bound was compared in case the optimal solution was not obtainable within limited time. The heuristics were tested for five runs, then the best solution reported among the five runs. Each method was set a number of iterations with 1000 as the stopping criterion. The mathematical model and the proposed algorithms were coded by Visual studio 2019 on PC Intel Core i5 CPU 2.70 GHz Ram 4 GB operated on Window 8.1 Pro. All statistical tests were conducted at the 95% confidence level.
The algorithm consists of two methods of black box selection and track improvement. The combinations of the proposed algorithms are defined, as detailed in Table 10. Table 10. Name of the proposed algorithms

Algorithms Name Definition
VaNSAS-1 Using black box selection equation (12) + track improvement equation (14) VaNSAS-2 Using black box selection equation (12) + track improvement equation (15) VaNSAS-3 Using black box selection equation (13) + track improvement equation (14) VaNSAS-4 Using black box selection equation (13) + track improvement equation (15) In order to assess the value of parameters, the DOE method of factorial design was used via Minitab software (Minitab Inc.). Figure 9 shows   (15) In order to assess the value of parameters, the DOE method of factorial design was used via Minitab software (Minitab Inc.). Figure 9 shows the optimization plot result, the value of parameters which are used in the proposed algorithms are shown as follow: F = 0.90, K = 0.10, M = 0.90, α = 0.90 and β = 0.90.  Table 11 shows the testing results of all sized instances. For small-sized instances, both Lingo and proposed algorithms could find the optimal solutions within short computational times. Because the problems were small, the difference was not obvious.
For medium-sized instances, Lingo could not find the optimal solutions within a computational time of 48 h, the solutions were only the best it could find during the time restriction. On the other hand, proposed algorithms could find the better solution within shorter computational time. Table  12 illustrates percentage gap between two methods that calculated solutions using equation 16. The proposed VaNSAS obtained the percentage gap of 3.05%, 3.02%, 2.45%, and 2.82%, respectively. The statistical test in Table 13 shows that all proposed VaNSAS perform significantly differently from Lingo's solution and all VaNSAS are non-significantly different from each other. Besides, the computational time of VaNSAS outperformed Lingo, consuming an average time of only 2.5 min, 2.6 min, 2.2 min, and 2.27 min, respectively.  Table 11 shows the testing results of all sized instances. For small-sized instances, both Lingo and proposed algorithms could find the optimal solutions within short computational times. Because the problems were small, the difference was not obvious.
For medium-sized instances, Lingo could not find the optimal solutions within a computational time of 48 h, the solutions were only the best it could find during the time restriction. On the other hand, proposed algorithms could find the better solution within shorter computational time. Table 12 illustrates percentage gap between two methods that calculated solutions using equation 16. The proposed VaNSAS obtained the percentage gap of 3.05%, 3.02%, 2.45%, and 2.82%, respectively. The statistical test in Table 13 shows that all proposed VaNSAS perform significantly differently from Lingo's solution and all VaNSAS are non-significantly different from each other. Besides, the computational time of VaNSAS outperformed Lingo, consuming an average time of only 2.5 min, 2.6 min, 2.2 min, and 2.27 min, respectively.
where %gap is the percentage difference between Lingo and VaNSAS, S P is the solution obtained from the proposed VaNSAS, and S L is the solution obtained from Lingo. For large-sized instances, Lingo used 120 h to find the solutions, but it could find only lower bound. However, the VaNSAS used only 9-10 min to find the reasonable solutions. Although these solutions may not be optimal, the statistical test in Table 14 reveals that there is a non-significant difference between Lingo's solutions and VaNSAS's solutions. Moreover, all VaNSAS are non-significantly different from each other in solving the large problem.  The percentage gap results in Table 12 were analyzed by using the ANOVA method to validate whether VaNSAS executed differently when the problem size was larger. The testing results are shown in Table 15. From Table 15, when the problem size was larger, all proposed VaNSAS performed significantly differently from Lingo's solutions. Nevertheless, the statistical results indicated that the solutions from VaNSAS were non-significantly different from Lingo's solutions in small-sized and large-sized instances. Furthermore, the VaNSAS obtained the better results than Lingo in medium-sized instances within limited time.
VaNSAS-1 showed the best performance when considered percentage gap in all sized instances. Therefore, this approach has been implemented to the case study problem and Table 16 shows the results of the case study. There are nine procurement centers which buy and collect the rubber from all farms across 22 transportation routes. The total fuel consumption cost is minimized to 13,280.45 baht. The special constraint that a farm can be visited more than once was applied in route nos. 10 and 20; farm nos. 4 and 47 were visited twice because they had more rubber than the vehicle's capacity.
The last experiment was conducted to show VaNSAS-1 performance by verifying the solution plot in Figure 10. The main idea of VaNSAS was to let the algorithm constantly change the searching area by applying diversification and intensification concepts. It can escape from the trap of local optima, then move to other areas, as implied by the turning point of the graph line. Therefore, VaNSAS can find often the new optimal solution, it can find new best solution during the all simulation time.

Conclusions and Future Research
This research has presented a new solution approach for solving location routing problems in agricultural logistics. The problem contained two main stages: (1) the selection of rubber procurement locations, and (2) the routing management of rubber collection. In this problem, the farm could generate more rubber than the vehicle's capacity, so this paper added the special constraint of letting the farm be visited more than once. The mathematical model has been formulated to represent the case study problem with the objective function to minimize the total fuel consumption cost.
The heuristic proposed in this research was called the variable neighborhood strategy adaptive

Conclusions and Future Research
This research has presented a new solution approach for solving location routing problems in agricultural logistics. The problem contained two main stages: (1) the selection of rubber procurement locations, and (2) the routing management of rubber collection. In this problem, the farm could generate more rubber than the vehicle's capacity, so this paper added the special constraint of letting the farm be visited more than once. The mathematical model has been formulated to represent the case study problem with the objective function to minimize the total fuel consumption cost.
The heuristic proposed in this research was called the variable neighborhood strategy adaptive search (VaNSAS). The VaNSAS comprised of three searching strategies which were: (1) swap, (2) adaptive large neighborhood search (ALNS), and (3) variable neighborhood search (VNS). These strategies were selected randomly to improve solutions in each searing iteration. Furthermore, we split the proposed heuristic into four sub-heuristics; VaNSAS-1, VaNSAS-2, VaNSAS-3, and VaNSAS-4, to evaluate the performance of the black box selection formula and track improvement formula.
The computational result showed that VaNSAS-1 outperformed the other proposed heuristics because it obtained the smallest gap when compared to Lingo's results. This implied the mixing of black box selection Formula 12 and track improvement Formula 14. Formula 12 selected a black box based on the intensification concept because it included a weighting of the best objective obtained from the black box. Thus, the best black box in the last iteration had more chance to be selected repeatedly. Formula 14 improved track information based on the diversification concept since it combined with the two random track's information to create the new track. Therefore, VaNSAS-1 has been selected to solve the case study problem. The results shown that nine farms were selected to be rubber procurement centers and 22 transportation routes were designed with the lowest total fuel cost. The computational experiment indicates that the new optimization technique can solve this case study problem effectively. The solution shows several farms were served twice according to the special constraint. These results describe the contributions of this research.
Akararungruangkul and Kaewman [9] applied modified differential evolution (MDE) for special case of LRP with fuel used minimization. For case study problem which contain 110 farms, the proposed MDE can save 15.05% of fuel amount comparing to current practice. The extended experiment of VaNSAS-1 reveals that 5.48% of fuel cost can be saved from traditional LRP. However, the problem size and the case problem are different. The larger the problem is, the greater savings it can obtain. The different distancesin different cases effect the comparison of both studies as well.
This paper based on static supply assumption; the rubber quantity of each farm is predefined number acquired from historical data. This is the limitation because the location and routing plan may not be suitable when the rubber quantity is changed. Besides, the data collection needs to be prompt and update in real-time. For future research, the authors suggest that another real-world constraint should be considered, such as stochastic supply according to daily rubber tapping of farmers. This will affect the selection of procurement center locations and routing paths as a daily dynamic move. A mobile application could be developed to support the gathering of daily rubber amounts from each farm, and the data will be sent for processing via a proposed algorithm to obtain the suitable procurement center locations and routes. Another suggestion is to apply VaNSAS to other optimization problems, allowing this algorithm to be opened for innovation in other industries and organizations.