A Pseudo-Parallel Genetic Algorithm Integrating Simulated Annealing for Stochastic Location-Inventory-Routing Problem with Consideration of Returns in E-Commerce

Facility location, inventory control, and vehicle routes scheduling are three key issues to be settled in the design of logistics system for e-commerce. Due to the online shopping features of e-commerce, customer returns are becoming much more than traditional commerce.This paper studies a three-phase supply chain distribution system consisting of one supplier, a set of retailers, and a single type of product with continuous review (Q, r) inventory policy. We formulate a stochastic location-inventory-routing problem (LIRP) model with no quality defects returns. To solve the NP-hand problem, a pseudo-parallel genetic algorithm integrating simulated annealing (PPGASA) is proposed. The computational results show that PPGASA outperforms GA on optimal solution, computing time, and computing stability.


Introduction
With the booming development of e-commerce industry, especially that B2C online shopping has become a part of people's daily life, returning goods is becoming more and more common.Studies show that there are more customer returns in e-commerce than in the traditional business, and the return rate sometimes can be up to 35% of the initial order in online selling [1,2].Considering the fact that the returned goods in online selling are highly integrated, they usually can reenter market channel after a simple repackaging process instead of depot repair.Here, logistics system as an important support system in e-commerce needs to make some adjustments and improvements.To adapt to the reality of e-commerce market environment, reverse logistics network and highly integrated logistics process should be the necessities.
Facility location, inventory control, and vehicle routes scheduling are three core issues to be settled in strategic/tactical and operational decision levels for logistics system in e-business.Previous work on these three areas is fruitful.In fact, there is a mutually dependent relationship among these problems.Comprehensive optimizing and logistics activities management should be based on this relationship [3].According to this idea, besides location allocation problem and vehicle routing problem, the two-two integration such as inventory-routing problem (IRP), location-routing problem (LRP), and location-inventory problem (LIP) and three-integration problem (location-inventory-routing problem, LIRP) should be further researched.
Many researches concerning the IRP, LRP, and LIP are studied deeply and have made some abundant achievements.However, there are few researches on the integration of LIRP.Some researchers appeal to carry out research on LIRP [4,5].Max Shen and Qi [6] aimed at minimizing the expected total cost resulting in a nonlinear model and proposed an algorithm based on Lagrangian relaxation to solve the model.However, they ignored the routes scheduling when making transportation decision.Liu and Lee [7] are the earliest researcher in the LIRP field; they presented a model for single merchandise and multiretailers LRP with inventory control decisions and established a two-stage heuristic algorithm.But this algorithm is easy to be trapped in local optima.In order to deal with this problem, Liu and Lin [8] established a global optimum heuristics to avoid locally optimal solutions by utilizing an extra mechanism.Ahmadi Javid and Azad [9] considered LIRP with stochastic supply chain system.In the article, they established a heuristic method based on a hybridization of tabu search and simulated annealing to solve the model.Ahmadi-Javid and Seddighi [10] considered LIRP with multisource distribution logistics network.They presented a mixed-integer programming formulation and a three-phase heuristic solves the problem.Guerrero et al. [11] presented the inventory location-routing problem with deterministic demand and proved the capability of the hybrid heuristic algorithm to solve the problem.Results show important savings achieved when compared to decomposed approach and the capability of the algorithm to solve the problem.
Previous researches on reverse logistics mainly focus on independent activities about LIRP.Fleischmann et al. [12] and Jayaraman et al. [13] are interested in determining the location of recycling center with capacity constraints.In recent years, some researches on reverse logistics concerned the integrated system.Lieckens and Vandaele [14] applied a queuing mode in reverse logistics network to solve the facility location problem while considering the impact of inventory costs.Sahyouni et al. [15] developed three generic facility location models that account for both forward and reverse logistics network; Easwaran and Üster [16] proposed a mixed integer linear programming model to optimize the total cost that consists of location, processing, and transportation costs of the multimerchandise closed-loop supply chains.Srivastava [17] established a reverse logistics network optimization model to optimize the location-distribution problem and capacity decisions, and he pointed out that integrated optimization of processing, storage, transportation, and recycling merchandises is one of the directions of future research.The above researches did not consider the influence of returns.
Li et al. [18] proposed a closed-loop LIP model considering returns to minimize the total cost which is produced in both forward and reverse logistics networks.They presented an effective two-stage heuristic algorithm to solve this nonlinear mixed programming model.Li et al. [19] presented a LIRP model considering returns under e-supply chain environment and proposed the HGSAA algorithm as the solving approach.However, these researches did not consider the stochastic demand of customers.
In this paper, we develop a practical stochastic LIRP model considering returns in e-commerce and provide a pseudo-parallel genetic algorithm integrating simulated annealing (PPGASA).To our best knowledge, this work is the first step to introduce returns into the stochastic LIRP in ecommerce, which makes it become more practical.We also provide an effective algorithm named PPGASA to solve the model, which is based on pseudo-parallel genetic algorithm and simulated annealing.Results of numerical examples show that PPGASA outperforms genetic algorithm (GA) on optimal solution, computing time, and computing efficiency and stability.
The remainder of this paper is organized as follows.Section 2 proposes a nonlinear integrated programming model based on forward and reverse logistics networks about stochastic LIRP under e-supply chain environment.Section 3 designs the heuristic algorithm named PPGAISA.Section 4 presents the results of different experiments and corresponding analysis.Section 5 concludes this paper and discusses future research directions.

Problem Description.
Customer returns in e-commerce generally have a high integrity, which do not need to be repaired and can reenter the sales channels after a simple repackaging process [20].Therefore, distribution centers and recycling centers can be merged into merchandise centers (MCs).MC is responsible for distributing normal goods to the retailers of downstream; meanwhile the returned goods are collected to MCs.After repackaging treatment at MCs, returned goods become resalable normal goods.The operations mode of the system is shown in Figure 1.
The supply chain in this study consists of one supplier, multiple MCs, multiple retailers, and a single type of product with continuous review (, ) inventory policy, which is called a three-phase e-commerce logistics system.Firstly, we optimize system construction and select the number of points from retail outlets as MCs.Secondly, we determinate the optimal ordering policy of each delivery path to reduce the inventory cost.Thirdly, we coordinate the arrangements of vehicle routing.
The objective of this problem is to determine the quantity, locations, order times, and order size of MCs and arrange the routes.The final target is to minimize the total cost and improve the efficiency of logistics operations.The involved decisions are as follows: (1) location decisions, the optimal number of MCs, and their locations; (2) inventory decisions, the optimal order times, and order size on a route; (3) routing arrangement, the vehicles delivering goods, and collecting customer returns in the order.

Assumptions
(1) There is a single type of goods.
(2) The total demand on each route is less than or equal to the vehicle capacity.
(3) The vehicle type is homogeneous.
(4) Each route is served by one vehicle.
(5) Each route begins and ends at the same MC.
(6) The capacity of each MC is finite.(7) The forward distribution and reverse collection service could be met at the same time.
(8) The daily demand and returns of each retailer are stochastic and obey the normal distribution.
(9) The returned goods without quality defect.
(10) Returned goods are processed and repackaged at MCs.

Notations
Sets   : optimal order size for MC  on routing .

Model Formulation.
The cost of MC  consists of the following components.
(1) Location Cost.The construction cost of MC  is given by ∑ ∈J (  +  V )  .
(2) Inventory Cost.The annual average total demand is given by   = ∑ ∈I   , each of the order size is given by   =   /  , and the average inventory is given by   /2  .By the central limit theorem, the annual total demand is ∑ ∈J   and standard deviation is √∑ ∈I  2  + 2 ℎ ∑ ,ℎ∈I    ℎ , where  ℎ is the correlation coefficient between retailers  and ℎ.The demands of retailers are independent and identically distributed; therefore  ℎ = 0.The demand variance is ∑ ∈I  2  , the safety stock levels is   √∑ ∈I  2  , the annual cost of placing an order at MC  is given by     , and the annual inventory holding cost at MC  is given by (3) Transportation Cost.The annual cost of the dispatching vehicle at MC  is given by  V   , the annual total distribution costs from every MC to each retailer is given by   ∑ ∈I ∑ ∈J ∑ ∈I ∑ ℎ∈I  ℎ    ℎ , and the annual transportation cost from plant to MC  is given by   ∑ ∈I ∑ ∈K   .
(4) Returns Processing Cost.The annual repackaging cost of returned goods at MC  is given by   ∑ ∈I   .
The objective is to minimize the total cost of the system; we formulate the model  as follows: The objective function ( 1) is to minimize the system's total cost.Constraint (2) ensures that only the retailer who is chosen as MC can provide services.Constraint (3) ensures vehicle cannot be overloaded.Constraint (4) ensures that the total demand of the retailer on the MC path should not exceed the capacity of MC.Constraints ( 5) and ( 6) ensure the continuity of delivery routes.Constraint (7) ensures each retailer is served by the only one vehicle which belongs to a certain MC.Constraint (8) ensures that each route only belongs to one MC.Constraints ( 9)-( 11) ensure the integrality of decision variables.

Solution Approach
The abstract idea of solution approach is described as follows.Firstly, we derive optimal order times  *  and optimal order size  *  from Kuhn-Tucker condition (K-T conditions), which also rely on the decision variable   .Secondly, we design a hybrid algorithm based on pseudo-parallel genetic algorithm (PPGA) and simulated annealing (SA) algorithm to obtain the optimal  *  ,  *  , and  * ℎ .
3.1.Inventory Optimization.In order to calculate the minimum of the objective function, make the following operation based on the model.Constraints ( 3)-( 4) are written as follows: In order to obtain the economic order size and the optimal order time, we calculate the gradient of objective function (1) and constraints ( 12) and ( 13): where   =  V +   .We introduce the generalized Lagrange multiplier  *  and  *  into the constraints ( 12) and ( 13), respectively.Let  *  be the K-T point; the K-T conditions are described as follows: where   =   +   .
Let  *  =  *  = 0; then  *  = √ ℎ ∑ ∈I     /2  .Since the original problem is convex programming, the cost of the whole system reaches the minimum value when   =  *  and the optimal order size is into the objective function (1), we can get the model   as follows: (2)∼( 11) , where

Hybrid Algorithm.
As we know, the VRP is an NP-hard problem.The LIRP containing the VRP is more complicated.
It is generally believed that there is no complete, accurate, and efficient analytic algorithm to solve NP-hand problems.Note that bioinspired computation for solving optimization problems is widely applied; we design a hybrid algorithm based on PPGA and SA to solve the proposed model.Although traditional genetic algorithm (GA) has strong global search ability in solving optimization problems, it has defects such as premature and weak local search ability.On the other hand, SA has strong local search ability without premature problem.Therefore, the combination of GA and SA can overcome the defects of each of the two methods, bring into play their respective advantages, and improve the solving efficiency.This hybrid algorithm is named pseudoparallel genetic algorithm integrating simulated annealing (PPGASA).

Relevant Operations of PPGA.
Parallel genetic algorithm (PGA) can improve the operation speed of GA and maintain population diversity, inhibiting the occurrence of "premature" phenomenon.Generally, the PGA is executed on a parallel machine or a LAN.But a single processor machine is applicable if real-time requirements are not needed, which is so called pseudo-parallel genetic algorithm (PPGA) [21].
(1) Encoding.In the GA, a population of candidate solutions (called individuals) to an optimization problem is evolved toward better solutions.GA cannot directly handle the settings of the problem space, and it must transform them into chromosome or individual.Binary coding is the most commonly used method of coding GA.That is, a binary character set {0, 1} represents the problem space of candidate solutions.
This study adopts the natural number coding method; using unrepeated  +  natural numbers constitute a sequence, which represents an individual.While 1, 2, . . .,  indicate candidate MCs and +1, . . ., + indicate retailers, the encoding can describe a candidate solution of the above optimization problem.
(2) Fitness Function.The fitness function of GA is called the evaluation function and is used to determine the degree of individuals in the population of the index, which is based on the objective function for the problem to evaluate.In this study, the fitness function is defined as ( (4) Crossover.Crossover refers to the part of the structure of two parent individuals to replace the restructuring and generate new individual operation.Through the crossover, search ability of GA can be leap increase.Partially matched crossover (PMX) [23] is used in this paper.
Step 1. Select two parent individuals randomly from the population.
Step 2. Generate two random cut points to represent the mapped segments.
Step 3. Exchange the segments of the two parent individuals to produce two new individuals.
Step 4. Determine the mapping relations between two segments.
Step 5. Legalize two new individuals with mapping relationship through repair strategy.
(5) Mutation.The basic content of mutation operators is to change the value for some loci of individuals in the population on the genes.A simple and efficient mutation operation, that is, swap mutation [24], is used.The details are as follows.
Step 1. Select one parent individual randomly from the population.
Step 2. Generate two random numbers to represent the mutation points.
Step 3. Swap the positions of these two mutation points to produce a new individual.
Compared with other mutations, studies show that convergence rate of this method has a greater advantage in control population.It can effectively prevent premature convergence of GA and avoid the occurrence of local optimal solution.

Dual-Channel Exchange Mechanism.
We calculate the fitness function values using the initializing population, sort the individuals by descending of the fitness function values, and then divide the population equally into two groups.The one with larger values is called successful subgroup, and the other one is called unsuccessful subgroup.In order to maintain the independent evolution of the two subgroups and balance the overall exploring ability and the convergence speed, information exchange should occur only in the appropriate stages.After -generation's independent evolution, the dual-channel exchange mechanism is invoked as follows.
Step 1. Select the individual  with the minimal fitness function value from the unsuccessful subgroup.
Step 2. Select randomly individual  from the successful subgroup.
Step 3. If   <   ,  enter successful subgroup and  enter unsuccessful subgroup; otherwise, remain unchanged.
Step 4. Select the individual  with the maximal fitness function value from the successful subgroup.
Step 5. Select randomly individual  from the unsuccessful subgroup.
Step 6.If   >   ,  enter unsuccessful subgroup and  enter successful subgroup; otherwise, remain unchanged.

Relevant Operations of SA (1) The Annealing Process of Accepting the New Individual.
In order to prevent the population from local optimization, the Metropolis acceptance criteria in SA are applied into the PPGA in this paper.If the new solution is greater than the old solution, then accept new individual and hold it to the next generation.Otherwise, the new individual is accepted with the probability  = exp(−(1/ new −1/ old )/) > random digit, where  is annealing temperature,  new is the objective function value of the new solution, and  old is the objective function value of the old solution.
(2) Temperature Amended Criterion.One of the key steps in the process of SA is to determine the update function of temperature.The function is used to continuously reduce the temperature value.When its temperature reduces to approximately zero, the final solution is considered as the global optimal solution.The updated function is  +1 =   ,  ≥ 0, 0 <  < 1; the nearer  is to 1, the slower the temperature decreases.

Termination.
In this paper, the termination condition is that the fitness has reached a plateau such that successive  iterations no longer produce better results.

Algorithm Flow
Step 0. Get the formulas for solving  *  and  *  .
Step 1. Set the initial parameters: coordinates of the retailers and the candidate MCs; demands and returns of the retailers; the maximum capacity of the vehicle ; the population size ; evolution terminate generation ; crossover probability   ; mutation probability   ; temperature of the cooling coefficient ; the initial annealing temperature ; and so on.
Step 2. Calculate the fitness value   of an individual.If the parent optimal solution and offspring optimal solution are equal during continuous  generations, the algorithm stops and outputs the current optimal solution; otherwise, go to the next step.
Step 3. Divide the population into successful and unsuccessful subgroups by individual fitness value.Each subgroup independently performs selection, crossover, and mutation operations.After each independent evolution of  generations, the dual channel exchange mechanism is executed.
Step 5. Update the annealing temperature, and return to Step 2.
The pseudocodes of HGSAA are shown in Pseudocode 1.

Computational Experiments and Results Analysis
An example is used to illustrate the proposed heuristic method.The data of Gaskell 67 − 29 × 5 comes from the LRP database at University of Aveiro [25].Gaskell 67 is the name of this instance; 29 × 5 means there are 29 retailers and 5 candidate MCs.The coordinate of all nodes and the demands of retailers are given by the database.Daily demands of customers are uniformly generated from  (20,1).The other parameters are as follows: the inventory holding cost per unit of goods per year ℎ = 6, the vehicle capacity  = 500, the  , d)    In this subsection, we will discuss five parameters' impact on the PPGASA.The performances of PPGASA vary with the different values of these parameters, which are shown in Tables 1-5.
Figures 2-6 represent the parameters' effect on the objective function values.The data was normalized through two dimensions, that is, cost and time, and three indicators, that is, mean, standard deviation, and coefficient of variation.Actually, in order to facilitate the discovery of regular pattern, diagram (b) in each figure is a different view of diagram (a), which is resorted by the best objective function value.From these figures, we found that the performance of PPGASA reaches the best when   = 0.7∼0.9,  = 0.0005∼0.001,q = 0.85∼0.9,T = 100, and  = 20.

Experiment Results.
For getting the results of the instance, we run another 100 times on the same computer.One of the minimum values of objective function in the 100 experiments is 963850.Table 6 shows the solution.Figure 7 shows topological structure of the supply chain.MCs were established at  1 ,  2 ,  3 ,  4 ,  5 points and there are eight vehicles distribution routes.
For comparison, GA is programmed by Matlab 7.0 as well, and the instance Gaskel 67 − 29 × 5 was run 100 times on the same computer.Figure 8 shows the trend of optimal objective function values along with the evolution generations using PPGASA.Figure 9 shows the trend of optimal objective function value along with the evolution generations using GA.
The fluctuation cure of optimal objective function values obtained from the 100 times experiments is shown in Figures 10 and 11, respectively.
The optimal objective function value and the CPU time of these two algorithms are shown in Table 7.As shown in this table, the cost of PPGASA is significantly lower than GA ( ≪ 0.05), and the efficiency of algorithms measured by CPU time in seconds is significant difference between the two algorithms ( ≪ 0.05).

Extended Experiments.
In this section, a series of experiments are given to show that PPGASA is more efficient and stable than GA.Similar with Section 4.2, all the data in our experiments comes from LRP database of the University of Aveiro [25].In order to ensure the demands of retailers are not more than the vehicle capacity, we need to enumerate some instances.In this study, the daily demands are set as 1/15 of corresponding demands of Gaskell 67 − 22 × 5.
Each instance was calculated 100 times by PPGASA and GA, respectively; the results are shown in Tables 8 and 9.
4.4.Result Analysis.Tables 8 and 9 show the mean, standard deviation, coefficient of variation, and result of significance test of 7 instances in different scales on cost optimization and CPU time.
Observe from Tables 8 and 9 that the mean cost and CPU time of PPGASA are less than GA, which means our approach obtains better approximate optimal solution with less computational cost.Further, our approach is more stable

Conclusion and Future Research
There exists a higher return rate in e-commerce than in traditional business.Generally, the returned goods have no quality defect and have great integrity.Just after a simple repackaging process, the returned goods can reenter the sales channels, which put forward high requirements to the logistics system that support the operation of e-commerce.In this paper, we formulate the stochastic LIRP model with returns in e-commerce and provide an effective heuristic algorithm named PPGASA to solve it.The main contributions are summarized as follows.
(1) We establish a stochastic LIRP model considering no quality defects returns in e-commerce to minimize the total cost of both forward and reverse logistics networks.It is very useful to help managers make the right decisions for e-commerce operations.
(2) The stochastic LIRP model with returns is a NP-hand problem and very hard to be solved by traditional   (3) Results of experimental data show that PPGASA outperforms GA on optimal solution, computing time, and computing stability.PPGASA is an excellent choice to solve the proposed LIRP model effectively.
While our model and approach increase the threshold of existing literature in this topic, we recognize a number of ways our research could be embellished by.Some extensions can   be done in further work.Considering the variety of the types of goods, the multiple products model should be proposed.In reality, decision makers are always in front of imprecise and vague operational conditions [26].Considering the fuzzy demand of customs or related fuzzy costs, more practical LIRP model should be developed.Moreover, differential evolution algorithms (DEs) have turned out to be one of the best evolutionary algorithms in a variety of fields [27,28].In the future, we may use an improved DE to find better solutions for the LIRPs.

Figure 6 :
Figure 6: Optimization results on different values .

Figure 7 :
Figure 7: Topological structure of the network.

Figure 8 :
Figure 8: Trends of optimal objective function value by GA.

Figure 9 :
Figure 9: Trends of optimal objective function value by PPGASA.

Figure 10 :
Figure 10: The fluctuation curve of optimal objective function value by PPGASA.

6 Figure 11 :
Figure 11: The fluctuation curve of optimal objective function value by GA.
[22]election.Selection of superior individual from the population and elimination of inferior individual operation is called selection.The objective of selection is to optimize the individual (or solution) directly inherited to the next generation or by paired crossover generate new individual inherited to the next generation.Selection is based on individual fitness evaluation.Wheel selection operator[22](also known as proportional selection operator) is used in this study.Suppose the population size is ; the fitness value of the individual  is   and a number  ∈ [0, 1] is generated randomly.If ∑ −1 =1   / ∑

Table 2 :
The results with different   (cost: million yuan, CPU time: second).

Table 3 :
The results with different  (cost: million yuan, CPU time: second).

Table 5 :
The results with different  (cost: million yuan, CPU time: second).

Table 7 :
Statistical results of optimal objective function value of two algorithms (cost: yuan, CPU time: second).

Table 8 :
Optimal objective function values of two algorithms (yuan).

Table 9 :
CPU time of two algorithms (seconds).