Three-Phase Methodology Incorporating Scatter Search for Integrated Production , Inventory , and Distribution Routing Problem

This paper proposes the use of scatter search metaheuristic to solve an integrated production, inventory, and distribution routing problem. The problem is based on a single production plant that produces a single product that is delivered to N geographically dispersed customers by a set of homogenous fleet of vehicles.The objective is to construct a production plan and delivery schedule to minimize the total costs and ensuring each customer’s demand ismet over the planning horizon.We assumed that excess production can be stored at the plant or at customer’s sites within some limits, but stockouts due to backordering or backlogging are not allowed. Further testing on a set of benchmark problems to assess the effectiveness of ourmethod is also carried out.We compare our results to the existing metaheuristic algorithms proposed in the literature.


Introduction
Two substantial components to improve the timeliness and consistency of delivery are integrated logistics management and product availability.Integrating production and distribution decisions can have a significant impact on setup, holding, and delivery costs.In general, the problem of coordinating production and transportation is called the productioninventory-distribution routing problem (PIDRP) [1].The research on PIDRP involves some different areas such as vehicle routing, production, and inventory [2].PIDRP is most similar to the inventory routing problem and periodic routing problem, because it requires multiple visits to each customer's sites over the planning horizon.
Chandra and Fisher [3] studied about a single plant, multicustomers in a multiple period by comparing two different methods to investigate the value of coordinating production and distribution.The problem is to minimize the total cost of production, transportation, and inventory.Two different alternative solution approaches are presented to manage this operation, one in which the production scheduling and vehicle routing problems are solved separately and another in which they are coordinated within a single model.The computational results reported a reduction in total operating cost from coordination ranged from 3% to 20%.
Lei et al. [1] presented multifacility, heterogeneous fleet version of the PIDRP that was motivated by a chemical manufacturer with international customers.The problem is to coordinate the production, inventory, and transportation schedules to minimize the total cost over the planning horizon.They proposed a two-phase methodology, by solving a restricted version of the problem by eliminating the routing constraints in phase one and proposed a routing heuristic based on an extended optimal partitioning procedure in phase two to transform the less-than transporter load assignments obtained in phase one into more efficient delivery schedules.
Boudia et al. [4] proposed integer linear model to solve the PIDRP but it failed to solve the large instances, and then a Greedy Randomized Adaptive Search Procedure (GRASP) is developed to tackle the production and distribution decisions simultaneously.Another two improved versions using either a reactive mechanism or a path-relinking process embedded in GRASP are also developed and the results are better than the GRASP alone.
Boudia and Prins [5] presented an alternative method to solve the PIDRP.They solved the problem with memetic algorithm population management (MA|PM) and compared with a two-phase heuristics and GRASP.Computational testing showed that MA|PM can tackle large instances and gave better results compared to two-phase approach and GRASP.
Bard and Nananukul [2,6] developed a mixed integer programming (MIP) model aimed at minimizing production, inventory, and delivery costs.The objective of the problem is to minimize the total cost over the planning horizon without incurring any stockouts at the customer sites.The problem includes a production plant, multiple customers with time varying demand, a finite planning horizon, and a fleet of homogenous vehicles.They developed a threephase methodology centered on tabu search and using allocation model to find a good initial solution.They also combined the features of reactive tabu search algorithm and branch-and-price algorithm by taking efficiency of the tabu search heuristic and the precision of the branch-and-price algorithm.Nananukul [7] extended the idea by improving the clustering of the customers by creating adaptive core clusters in the reactive tabu search algorithms which are used in the clustering process instead of the original data points thus enabling the algorithm to be efficient.Unfortunately the detailed computational results were not given.
Armentano et al. [8] presented two tabu search variants for PIDRP.The first variant involves construction of a shortterm memory and integrating a path relinking procedure, while another one incorporates a longer term memory and integrate the first variant.The algorithms are tested on generated instances and on instances taken from Boudia et al. [4] which involved a single product.Computational results showed that the two variants of tabu search yield good tradeoffs between solution quality and computational time and successfully outperformed Boudia and Prins [5] and Bard and Nananukul [2,9] in all instances.
Adulyasak et al. [10] improves upon the results of Armentano et al. [8] by proposing adaptive large neighborhood search heuristic to take care of binary variables representing the setup and routing variables whilst the continuous variables associated with inventory, production, and quantities delivered are handled by solving a network flow problems.The results outperformed all other know heuristics for PIDRP.
Two inventory replenishment policies, order-up-to level and maximum level were considered in Adulyasak et al. [11] for inventory routing problem and PIDRP.By using adaptive large neighborhood search to obtain the initial solutions and the branch-and-cut algorithm were proposed to solve the different formulations.The authors managed to solve to optimality relatively small instances, 50 customers, three periods, and three vehicle on parallel computers.
Torabi et al. [12] proposed a two-step solution approach to solve an integrated multisites production planning, procurement, and distribution plans.The first step restricts the vehicle routings into direct shipment and solved the full model as mixed integer problem.Scatter search algorithm is employed in the last step by solving the associated consolidation problem in order to improve the solutions.Computational testings showed that this method gives the comparable or improved solution compared to the best solution for original model for 8 out of 10 cases and failed to get a better solution for 2 other problems.However the algorithms were not performed on benchmark instances.
We refer the readers to a review of formulations and solution algorithms in PIDRP by Adulyasak et al. [13] for the state of the art.
In particular, the PIDRP considered in this paper is very similar to that of Bard and Nananukul [2,6] which involves a production plant, multiple customers with time varying demand, a finite planning horizon, and a fleet of homogenous vehicles that delivers the product from the production plant to the customers' sites.Excess production can be stored either at the plant or at the customer sites within some limits, but inventory cannot be transferred between sites and stockouts are not permitted.The objective is to minimize a combination of setup, holding, and routing costs both at the production facility and at customers, without incurring any stockouts at the customer sites.We also propose a three-phase methodology and our algorithm differs from Bard and Nananukul [2,6] in Phase 2 and Phase 3. In Phase 2, we employ the Giant Tour procedure [14], sweep, and savings algorithms which have been shown to be efficient routing procedures and Phase 3 introduces the scatter search algorithm, embedding the new inventory updating mechanism as the improvement methodology.
The remainder of this paper is organized as follows.In Section 2, we present the description of PIDRP and its mathematical formulation.In Section 3, the three-phase scatter search method we employed to solve the PIDRP is discussed in detail and followed by the presentation of computational experiments and results in Section 4. Finally, conclusions are drawn in Section 5.

Problem Description and Mathematical Formulation
We consider a production, inventory, and distribution routing problem similar to the one proposed by Bard and Nananukul [2,9].It consists of a single production plant that produces a single product and distributes it to a set of  customers with nonnegative demand   in period  where  = 1, 2, . . .,  and  = 1, . . .,  and limited number of items   can be produced in period  and a limited number of inventories    can be stored by incurring unit holding cost ℎ  at production plant.A fleet of homogeneous vehicles with capacity  delivers the items to the customer's sites, and each vehicle can make at most one trip per period.A limited amount of inventory    can be stored at customers' sites with unit holding cost of ℎ   , and each customer can only be visited at most once per period.
Furthermore, it is assumed that at the end of planning horizon all inventories (both at the production facility and at customer's site) are required to be zero.The objective is to construct a production plan and delivery schedule which minimizes production, inventory, and distribution costs while fulfilling each customer's demand requirement.Let  0 denote the production plant (depot) and { 1 ,  2 , . . .,   } be a set of customers.  is a decision variable denoting the amount to be delivered to customer  in period .The traveling distance (cost) from customer   to customer   (,  = 1, 2, . . ., ) is denoted by   .  equals 1 if there is a route from customer  to customer  (,  = 1, 2, . . ., ) and 0 otherwise.  represents the total quantity on a vehicle before delivering to customer  in period .
The maximum inventory level for production plant is denoted by   max and for customers' sites is denoted by   max, .In this study, the initial inventories at customers' sites are assumed to be zero.
The integrated production, inventory, and distribution routing problem (PIDRP) can be formulated as follows [2]: subject to ≥ 0 and integer, where =   }.The objective function comprises the transportation costs, production setup costs, holding costs at the warehouse and holding costs at the customer sites.Equations ( 2) and ( 3) represent the inventory flow balance equations for production facility and customers, respectively.Equation ( 5) limits production on period  to the capacity of the factory, and (6) allows production in period 0. The total amount available for delivery on period  is limited by the amount of inventories at the factory on period  − 1 as formulated in (4).Equation ( 11) limits the amount delivered to each customer and (7) ensures that if customer  is serviced on period , then it must have a successor on its route, while route continuity is enforced by (8).Equation ( 9) limits the number of vehicles that leaves the factory at period , and (10) keeps track of the load on the vehicles.We note that Adulyasak et al. [10] use slightly different approach in the formulation.
In this study, extending the idea from Bard and Nananukul [6] we propose a three-phase methodology to solve the PIDRP.Our algorithm differs from Bard and Nananukul [6] in Phase 2 and Phase 3. Starting from Phase 1 that solves the allocation model which is the simplified version (relaxed) of the model to determine the amount to be delivered to each customer in each period, Phase 2 routes the customers using the Giant Tour procedure [14], sweep, and savings algorithms to determine the delivery routes for each period.In phase 3, we develop scatter search method by creating composite decision rules and surrogate constraints to improve the initial solutions.We also incorporate the inventory updating based on the forward and backward transfer.
We identify the initial solution in Phase 1 by solving the allocation model as a mixed integer programming to get a set of feasible allocations.The routing variables and routing constraints ( 7)-( 10) are removed and aggregated vehicle capacity constraints are introduced to the allocation model.Since we already deleted the routing constraints in the allocation model, we need an alternative representation for the cost term to determine the approximated cost which is needed to make a delivery to customer  in period .

𝑓 𝐶
represents the fixed cost of making delivery to customer  on period ,    denotes the variable cost of delivering one item to customer  in period , and    is valued to 1 if delivery is made to customer  in period  and 0 otherwise.
As in Bard and Nananukul [2], we divide the problem into two cases, for problem instances with  2  ≤ 500 and for instances with  2  > 500.Since all the instances we considered are  2  > 500, we introduce the variable cost term ∑       , where    is approximated by the cost of making a delivery to customer  directly from the depot divided by the total demand of customer  in period .
The allocation model of the PIDRP is formulated as follows: Additional new constraints are as follows: Allocation model modifies the objective function in the full model and replaces the routing variables with additional parameter to represent the costs (the second term in (1a)).The term ∑       is the approximated transportation cost replacing the term ∑ ∈ ∑ ∈ 0 ∑ ∈ 0     , the actual distance (transportation) cost.Allocation model also uses the same constraints but eliminating the routing constraints (constraints ( 7)-( 9)) and adding additional two constraints.Equation ( 14) limits the total amount that can be delivered on period  to a fixed percentage of the total transportation capacity.The parameter testing by Bard and Nananukul [2] showed that the percentage value of 80% always yielded feasible solutions.Additional variable    is included in constraint (15) to keep track of whether customer  receives a delivery in period .Constraint (15) was a modification of constraint (11) by replacing the routing variables   with    .

Scatter Search
The scatter search metaheuristic has successfully been applied to a widespread variety of vehicle routing problems.Corberán et al. [15] proposed a scatter search to solve a real-life problem with multiple objectives.Two different heuristics were used to construct the initial trial solutions in the scatter search.Two simple exchange procedures, insertion and swap, are used to improve the solutions.The combination method is based on a voting scheme.The algorithms tested on real data show that scatter search can solve the practical problem efficiently.A more recent application of scatter search is by Mota et al. [16], who presented a scatter search for vehicle routing problem with split demands.Local search is adopted as the improvement method.Four kinds of critical clients are defined to produce new solutions.The algorithm was tested on a set of benchmark instances and it was found that scatter search algorithm always produces the least number of vehicles compared the tabu search developed by the authors.Scatter search is an evolutionary metaheuristic that operates on a set of solutions, which the scatter search literature refers to as the reference set (Refset).The evolution of the Refset is achieved by way of combining reference solutions to yield trial solutions with combination of attributes not present in the previous set of solutions.The Refset is a collection of "good" solutions found during the search, where the meaning of "good" is not limited to quality as measured by the objective function value.For instance, a solution may be good because it provides diversity with respect to other solutions in the reference set.In fact, some implementations of scatter search divide the Refset into two subsets, consisting, respectively, of solution quality and diversity.Scatter search was first introduced by Glover [17].Glover made a template in a version customized for nonlinear optimization problems with continuous variables.Laguna and Marti [18] published the first book on scatter search, containing introductory tutorials and advanced techniques such as the use of memory and path relinking.
The scatter search terminology that is used in this paper is similar to Laguna and Marti [18].The algorithm is made up of several distinct steps: A diversification generator, an improvement method, a reference set update, a subset generation method which operates on the reference set in order to produce a subset of its solutions as a basis for creating combined solutions, and a solution combination method which transforms a given subset of solutions produced by the subset generation method into one or more combined solutions vectors.
The scatter search procedure stops when a termination criterion-either the maximum number of iterations, Max-Iter, is reached, or the reference set does not change, or improvement does not warrant further iterations.
The scatter search algorithm can be formally stated as shown in Algorithm 1.

Updating the Inventory Level (Step 7).
In the inventory updating we propose two types of moves: a forward and backward transfer.The aim of the forward transfer is to reduce the inventory holding cost without increasing drastically the transportation cost.In the backward transfer the preference is given to the suppliers with the lower holding cost in order to determine whether the transportation and the inventory holding cost can be further consolidated.Examples of the forward and backward transfers are illustrated in Figures 1  and 2, respectively.Backward Transfer.Figure 1 shows an example of backward transfer where the routings are separated by zeros and   and   are the pick-up quantity and the inventory, respectively.Assuming the coordinates of the 5 customers are  1 (−2, 0),  2 (1, 1),  3 (4, 2),  4 (3, 5), and  5 (1, −2) and the depot is located at (0, 0) the holding costs per unit for each customer are ℎ 1 = 12, ℎ 2 = 9, ℎ 3 = 6, ℎ 4 = 3, and ℎ 5 = 6 and the vehicle capacity is 10.
Step  The selection of period and customer to be transferred is favorable to the lower holding cost.In this example, we select customer 4 in period 5.The saving is found by increasing the inventory cost and decrease in routing.
Initial routings are for periods 4 and 5, with the route cost 18.771591 and 24.36395 and 0 holding cost.According to the inventory updating mechanism, we transfer customer 4 in the period 5 to period 4. As the same customer is visited in period 4, we aggregate the amount to be transferred with the existing delivery quantity and we note that the resulting aggregated amount does not violate the capacity constraint.After the transfer of delivery amount by 2 units, we have a holding cost 6 with ℎ 4 = 3.This will reduce the distance in period 5.The overall savings after the transfer is 64.2431 − 58.5812 = 5.6619.Forward Transfer.Figure 2 shows a forward transfer and the objective is to reduce the holding cost whilst achieving a balance between the inventory and the transportation costs.
The selection of period and supplier to be transferred is biased towards customers with high holding cost.In this example, we select customer in period 1.Note that we limit the transfer to at most 2 periods only.This is to ensure that the increase in the routing cost is not exceedingly high.
The demands for customer 1 in periods  the resultant holding cost for periods 1, 2, and 3 are 81, 24, and 36, respectively, and the total cost, including the routing cost for all 3 periods, is 240.3398.Customer 1 is not visited in periods 2 and 3, so we apply forward transfer by inserting customer 1 to period 3 according to the best insertion.Note that inserting customer 1 in period 2 results in the violation of vehicle capacity constraint.The saving after the transfer is 240.3398− 153.4129 = 86.9269.

Diversification Generation Method.
Diversification generation method generates a set of diverse solutions which is denoted as , and the size of the population of solutions is represented by  (i.e., || = ).Few different algorithms are applied within the diversification generation method to obtain initial solutions.This method produces feasible solutions that can be used as trial solutions for the scatter search procedure.The solutions are created using sweep algorithm [19], savings algorithm [20], and Giant Tour Procedure [14].The Giant Tour Procedure starts by constructing cost network, which considers delivery amounts and vehicle capacity constraints.The construction starts by calculating the cost from the depot 0 to the customer   and from customer   going back to the depot as the cost of the arc (0,   ).If the combined deliveries of customers   and   are less than the vehicle's capacity, then the arc of (  ,   ) is added into the existing cost network.The construction is continued until the vehicle could not accommodate the new customer, and we start with a new vehicle.The process is continued until there are no more arcs connecting to the last customer in the giant tour.The giant tour is then partitioned using Djikstra's algorithm to form routes.In order to improve the solutions obtained from sweep and saving algorithms, we applied Giant Tour Procedure following the order of the routes generated by sweep and savings algorithms.The procedure is able to improve on the solution by collapsing a few vehicles in each period.

Improvement Method.
The solutions generated by the diversification generation method and combination method (see Section 3.5) are subjected to the improvement method.We apply 1-0 exchange, 1-1 interchange, 2- * as interroute procedure, and 2- as intraroute in the improvement method to improve the solutions.
(1) 1-0 Exchange (insertion).This exchange method removes a customer form one route and inserts it into another route.Select a customer   ,  ∈ {1, . . ., ), from the route   in the period ,  ∈ {1, . . ., }, and insert it into another route   within the same period.The customer is inserted into another route if it does not violate any constraints and reduces the routing cost.
(2) 1-1 Interchange (swap).This method starts with removing two customers from their initial routes within the same period.The customer   ,  ∈ {1, . . ., ), from route   is exchanged with   ,  ∈ {1, . . ., ), from route   .The exchange is applied if it could produce a shorter length of the routing cost or reduce the number of vehicles.
(3) 2- * .2- * algorithm is similar to the 2-, but, instead of deleting edges within the same route, it deletes two different edges on the different routes and then reconnects them by considering the lower routing cost.
(4) 2-.This method selects two different edges within the same route and then deletes and reconnects them to other edges.The move is accepted if the resulting routing cost is lower than the previous cost.The process is continued until no further improvement is found.

Reference Set Update
Method. denotes the reference set and its size is denoted by  (i.e., || = ).Solutions are included in the reference set by a measure of quality (objective value) or diversity.Solutions are denoted by   , ( ∈ 1, . . .), and it is assumed that  1 is the best solution and   is the worst solution according to the objective value.The reference set consists of two subsets, in which the subset of high quality solutions denoted by  1 contains the  1 best solutions, and the subset of diverse solutions denoted by  2 contains the  2 diverse solutions; hence  =  1 +  2 .
The initial reference set consists of  1 best solutions that belong to  and a number of  2 solutions from  that maximize the minimum distance to .The distance between the two solutions is calculated by adding the number of noncommon arcs of each solution before the combination.We adopt Russell and Chiang [21] for updating the  by composing the new reference set of ( 1 +  2 )/2 best solutions from the original reference set, and the remaining solutions are created from newly generated (combined) solutions.The solutions in the reference set are not changed until all combinations of solutions, generated by the combination method, are performed as prescribed by the subset generation method.This indicates the use of a static update of the reference set, where the set of solutions generated by the combination method is denoted as  and the subset of solutions chosen from  is denoted by .

Subset Generation Method.
Subset generation method is the foundation for constructing new solutions in scatter search.The subsets are built based on the reference set.This method generates subsets of the solutions in the reference set which are combined in the solution combination method.In our implementation, we restrict the method to select representative subsets by using Subset Type-1 [22], which consists of all 2-element subsets comprising of two different solutions which is explained further in the next section.

Solution Combination Method.
We adapt the solution combination method in Torabi et al. [12] that use Subset Type-1 to create new solutions and applied to all pairs of subsets generated by subset generation method.Each subset comprises of two different solutions that is combined to generate a now solution.This method is divided into two steps.
Step 1 aims to combine only the common elements of the combined routes and Step 2 assigns the demand of the remaining customers.
Step 1.Let  be a solution with  routes and  a solution with  routes, where   is the th route for solution ,  = {1, . . ., }, and   is the th route for solution .The solutions  and  are combined as follows: (1.1) Build a matrix  ×  where its components (  ,   ) have the number of common arcs between the route  of solution  and route  of solution .
(1.2) Choose the components (  ,   ) which have the greater number of common arcs; if there is a tie, then select randomly.

Mathematical Problems in Engineering
(1. 3) The combined route is performed by the routes' common elements (arcs).
(1.4) Each combined route is excluded from the list (delete the line and column referring to components of (  ,   )).
Step 2. This step aims to fulfill the customers who have not been served by Step 1 and done through the insertion procedure.This procedure is done as follows: (2.1) Pick customer ,  ∈ (  ,   ), which is the farthest one from the depot and is going to be inserted.
(2.2) Based on the initial solutions  and , all the routes in which customer  is inserted are verified and also all the edges in which   = 1 (customer  is served before customer ) or   = 1 (customer  is served after customer ).
(2.3)For each  belonging to one of the combined routes in Step 1, we calculate the cost for inserting customer  before customer  or the cost for inserting customer  after customer .
(2.4) Choose each possible insertion with the minimum cost.If there is no feasible position to be inserted, construct a new route.
The procedure stops when it reaches the maximum iteration.

Results and Discussion
All the algorithms are written in Matlab 7.7 and performed on a 3.1 GHz processor with 8 GB of RAM.For the algorithm testing, we generate randomly 14 data sets to test the algorithm that comprises of 12, 20, 50, and 100 customers with 5, 10, 14, and 21 periods.The locations of the customers are generated randomly in a square of 100 × 100.The locations of the customers for the 20 customers are extended from the 12 customers instance by adding 8 new randomly generated customers.Similarly, the 50 customers instance is extended from the 20 customers by generating randomly 30 additional customers.The same procedure was applied for the 100 customers instance by generating randomly an additional 50 customers' locations.
All of the data sets have demands in every period, with exception for case 50 customers.The holding cost for each customer is generated within the range [1,10] and the demands are generated randomly within the range [0, 50].The vehicle capacity is fixed to 100 and the depot is located at (0, 0) for all data sets.
Table 1 shows the results, the number of vehicles, and the CPU time for scatter search metaheuristic.CPLEX 12.4 is allowed to run 10800 seconds (3 hours) but, in most cases, CPLEX terminates prematurely because of being out of memory.It is observed that CPLEX failed to get lower bounds for large sized problems.We note that the formulation presented in Section 2 is weak because of the imposed number of vehicles and is able to solve to optimality for a very small problem [9].However the CPLEX solutions are used as guideline in order to ensure that our solutions are correct.Table 1 illustrates the CPLEX solutions and the solutions obtained by our algorithm.Furthermore, we test the performance of scatter search algorithm on the set of instances generated by Boudia et al. [4] comprising of three subsets of 30 instances each with 50, 100, and 200 customers over a planning horizon of 20 periods.Every subset of instances has a limited number of vehicles, with 5, 9, and 13 vehicles, respectively.
All of the data sets have demands in every period.The holding costs for each customer and the production plant are assigned at 1 for all instances, and the vehicle's capacity is fixed at 8000 for instances with 50 and 100 customers and 12000 for instances with 200 customers.The production capacities are fixed at 50000, 120000, and 240000 for instances with 50, 100, and 200 customers, respectively, and the storage capacity at the plant varies between one and a half and twice of the production capacity.The production setup cost is proportional to production capacity.(2007) and it also performs well for large instances when compared to the reactive tabu search algorithms of Bard and Nananukul [2].However when compared to Armentano et al. [8] the difference between the results obtained by our algorithm is between 4% and 6% as worst off.Our algorithm outperforms only in one instance, instance 1 in 50-customer problem.When compared to the algorithms of Adulyasak et al. [10] our algorithm performs poorly and they are in between 10% and 12%.

Conclusion
In this paper, we present scatter search metaheuristic algorithms for a finite horizon, multiperiod, and multiproduct PIDRP with no-split deliveries and no-backlogging.We propose three-phase methodology that implements the allocation model in the first phase and construct the routes in Phase 2 for the vehicle routing problem.In Phase 3, we develop a scatter search algorithm by creating composite . 2 Boudia and Prins [5]. 3 Bard and Nananukul [2]. 4 Armentano et al. [8]. 5 Adulyasak et al. [10,11].
decision rules and surrogate constraints to improve the initial solutions.Phase 1 is similar to Bard and Nananukul [2] but second phase uses the giant tour, savings, and sweep algorithms to construct the routes.We observe that our algorithm is superior when compared to the memetic/population management of Boudia and Prins [5] for the 50 customers instances but performs quite poorly when compared to the reactive tabu search of Bard and Nananukul [2].For large instances, the 200customer problem, our algorithm is very much superior when compared to Bard and Nananukul [2], producing better results in 28 out of 30 instances and comparable to Boudia and Prins [5].However our algorithm performs poorly when compared to Armentano et al. [8] and Adulyasak et al. [10].Adulyasak et al. produced the best known results so far.For future research we need to further fine-tune our algorithm in order to be competitive with the current best known results.

Figure 1 :
Figure 1: Example of a backward transfer.

Figure 2 :
Figure 2: Example of a forward transfer.
1. Generate Solutions -Generate trial solutions using Giant Tour Procedure, sweep algorithm and savings algorithm Step 2. Improve the Solutions -Apply 2- to improve solutions generated in Step 1 Step 3. Build the reference set -put |Refset| =  1 +  2 solutions in the reference set Step 4. Initialize best solutions -make  as the best solution in the current Refset, Generate subset -generate 2 pairs of solution (2 element subset) of the Refset.Step 6. Combine solutions -generate new combined solutions from pairs of 2 element subsets in Step 5.Step 7. Improve the solutions -update inventory level, and for the affected periods (transfer from and transfer to periods) apply within the same period 1-0 and 1-1 exchange, and 2-.Step 8. Update reference set -update the reference set by maintaining the number of solutions inside the Refset by replacing the existing solutions with the better combined solutions.

Table 1 :
Results from randomly generated data sets.

Table 2 :
Results for 50 customers instances.

Table 4 :
Results for 200 customers instances.
comparable to MA-MP algorithms of Boudia and Prins

Table 5 :
Average total costs obtained by different heuristics.