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

Authors: Bailing Liu; Hui Chen; Yanhui Li; Xiang Liu
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

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.


---

## Body

## 1. 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 e-commerce, 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. Section2 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.


## 2. Mathematical Model

### 2.1. 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.Figure 1
Diagram of the LIRP network.The supply chain in this study consists of one supplier, multiple MCs, multiple retailers, and a single type of product with continuous review(
Q
,
r
) 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.


### 2.2. 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.


### 2.3. Notations

Sets
I: set of retailers;
J: set of potential MCs;
K: set of routes.Constants
a: the order lead time;
b: the maximum capacity of the vehicle;
d
i: the mean (daily) demand at retailers;
e
j
k: the annual average demand in the service of route k at MC
j;
f
j: the fixed cost at MC
j, j
∈
J;
f
v: the management cost at MC
j;
w
j: maximum capacity at MC
j;
s
g
h: the distance from every MC to each retailer;
c
j: the unit of product transportation cost from plant to MC
j;
c
h: the inventory holding cost per unit of merchandise per year;
c
b: repackaging cost of unit returned merchandise at MC
j;
c
o: fixed cost of placing an order at MC
j;
c
v: fixed cost of dispatching vehicles per time at MC
j;
c
m: the delivering cost per unit distance;
f
l
(
x
): retailers demand probability density function in the lead time a;
q
i: the quantity of merchandises returned by the retailer i per day;
α: the inventory level of service;
λ: working days per year.Decision Variables
X
j: 1, if candidate MC
j is selected as a MC location, and 0 otherwise;
Y
i
j
k: 1, if the retailer i is assigned to MC
j on routing k, and 0 otherwise;
Z
g
h
j
k: 1, if node g is from node h served by a MC
j on routing k, and 0 otherwise;
N
j
k: optimal order times for MC
j on routing k;
Q
j
k: optimal order size for MC
j on routing k.


### 2.4. Model Formulation

The cost ofMC
j consists of the following components.(1) Location Cost. The construction cost of MC
j is given by ∑
j
∈
J
(
f
j
+
f
v
)
X
j.(2) Inventory Cost. The annual average total demand is given by e
j
k
=
∑
i
∈
I
λ
d
i, each of the order size is given by Q
j
k
=
e
j
k
/
N
j
k, and the average inventory is given by e
j
k
/
2
N
j
k. By the central limit theorem, the annual total demand is ∑
j
∈
J
e
j
k and standard deviation is ∑
i
∈
I
σ
i
2
+
2
ρ
g
h
∑
g
,
h
∈
I
σ
g
σ
h, where ρ
g
h is the correlation coefficient between retailers g and h. The demands of retailers are independent and identically distributed; therefore ρ
g
h
=
0. The demand variance is ∑
i
∈
I
a
σ
i
2, the safety stock levels is z
α
∑
i
∈
I
a
σ
i
2, the annual cost of placing an order at MC
j is given by c
o
N
j
k, and the annual inventory holding cost at MC
j is given by ∑
i
∈
I
∑
k
∈
K
c
h
(
λ
(
d
i
+
q
i
)
Y
i
j
k
/
2
N
j
k
+
z
α
a
∑
i
∈
I
Y
i
j
k
σ
i
2
).(3) Transportation Cost. The annual cost of the dispatching vehicle at MC
j is given by c
v
N
j
k, the annual total distribution costs from every MC to each retailer is given by λ
c
m
∑
i
∈
I
∑
j
∈
J
∑
g
∈
I
∑
h
∈
I
s
g
h
d
i
Z
g
h
j
k, and the annual transportation cost from plant to MC
j is given by c
j
∑
i
∈
I
∑
k
∈
K
e
j
k.(4) Returns Processing Cost. The annual repackaging cost of returned goods at MC
j is given by λ
c
b
∑
i
∈
I
q
i.The objective is to minimize the total cost of the system; we formulate the modelP as follows:(1)
P
:
min
⁡
Z
=
∑
j
∈
J
f
j
+
f
v
X
j
+
∑
j
∈
J
∑
i
∈
I
∑
k
∈
K
c
j
e
j
k
Y
i
j
k
+
∑
j
∈
J
∑
i
∈
I
∑
k
∈
K
λ
c
b
q
i
Y
i
j
k
+
λ
∑
j
∈
J
∑
k
∈
K
∑
g
∈
I
∑
h
∈
I
c
m
s
g
h
d
i
-
q
i
Z
g
h
j
k
+
∑
j
∈
J
∑
k
∈
K
c
o
+
c
v
N
j
k
+
∑
i
∈
I
∑
j
∈
J
∑
k
∈
K
c
h
+
z
α
a
∑
i
∈
I
Y
i
j
k
σ
i
2
λ
d
i
+
q
i
Y
i
j
k
2
N
j
k
+
z
α
a
∑
i
∈
I
Y
i
j
k
σ
i
2
(2)
s
.
t
.
Y
i
j
k
≤
X
j
,
∀
j
∈
J
,
∀
k
∈
K
,
∀
i
∈
I
,
(3)
hhhhh
λ
∑
i
∈
I
(
d
i
+
q
i
)
Y
i
j
k
N
j
k
≤
b
,
∀
k
∈
K
,
∀
j
∈
J
,
(4)
hhhhh
∑
k
∈
K
∑
j
∈
J
e
j
k
Y
i
j
k
N
j
k
≤
w
j
,
∀
j
∈
J
,
(5)
hhhhh
∑
h
∈
I
Z
g
h
j
k
≤
Y
g
j
k
,
∀
g
∈
Ι
,
∀
j
∈
J
,
∀
k
∈
K
,
(6)
hhhhh
∑
g
∈
I
Z
g
h
j
k
≤
Y
h
j
k
,
∀
h
∈
I
,
∀
j
∈
J
,
∀
k
∈
K
,
(7)
hhhhh
∑
j
∈
I
∑
k
∈
K
Y
i
j
k
=
1
,
∀
i
∈
I
,
∀
j
∈
J
,
(8)
hh
∑
j
∈
X
j
=
1
,
j
∈
J
,
Y
i
j
k
≤
1
,
∀
i
∈
I
,
∀
k
∈
K
,
(9)
hhhhh
X
j
∈
0,1
,
j
∈
J
,
(10)
hhhhh
Y
i
j
k
∈
0,1
,
∀
i
∈
I
,
j
∈
J
,
∀
k
∈
K
,
(11)
hhhhh
Z
g
h
j
k
∈
0,1
,
∀
g
∈
Ι
,
∀
h
∈
I
,
j
∈
J
,
∀
k
∈
K
.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.


## 2.1. 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.Figure 1
Diagram of the LIRP network.The supply chain in this study consists of one supplier, multiple MCs, multiple retailers, and a single type of product with continuous review(
Q
,
r
) 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.


## 2.2. 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.


## 2.3. Notations

Sets
I: set of retailers;
J: set of potential MCs;
K: set of routes.Constants
a: the order lead time;
b: the maximum capacity of the vehicle;
d
i: the mean (daily) demand at retailers;
e
j
k: the annual average demand in the service of route k at MC
j;
f
j: the fixed cost at MC
j, j
∈
J;
f
v: the management cost at MC
j;
w
j: maximum capacity at MC
j;
s
g
h: the distance from every MC to each retailer;
c
j: the unit of product transportation cost from plant to MC
j;
c
h: the inventory holding cost per unit of merchandise per year;
c
b: repackaging cost of unit returned merchandise at MC
j;
c
o: fixed cost of placing an order at MC
j;
c
v: fixed cost of dispatching vehicles per time at MC
j;
c
m: the delivering cost per unit distance;
f
l
(
x
): retailers demand probability density function in the lead time a;
q
i: the quantity of merchandises returned by the retailer i per day;
α: the inventory level of service;
λ: working days per year.Decision Variables
X
j: 1, if candidate MC
j is selected as a MC location, and 0 otherwise;
Y
i
j
k: 1, if the retailer i is assigned to MC
j on routing k, and 0 otherwise;
Z
g
h
j
k: 1, if node g is from node h served by a MC
j on routing k, and 0 otherwise;
N
j
k: optimal order times for MC
j on routing k;
Q
j
k: optimal order size for MC
j on routing k.


## 2.4. Model Formulation

The cost ofMC
j consists of the following components.(1) Location Cost. The construction cost of MC
j is given by ∑
j
∈
J
(
f
j
+
f
v
)
X
j.(2) Inventory Cost. The annual average total demand is given by e
j
k
=
∑
i
∈
I
λ
d
i, each of the order size is given by Q
j
k
=
e
j
k
/
N
j
k, and the average inventory is given by e
j
k
/
2
N
j
k. By the central limit theorem, the annual total demand is ∑
j
∈
J
e
j
k and standard deviation is ∑
i
∈
I
σ
i
2
+
2
ρ
g
h
∑
g
,
h
∈
I
σ
g
σ
h, where ρ
g
h is the correlation coefficient between retailers g and h. The demands of retailers are independent and identically distributed; therefore ρ
g
h
=
0. The demand variance is ∑
i
∈
I
a
σ
i
2, the safety stock levels is z
α
∑
i
∈
I
a
σ
i
2, the annual cost of placing an order at MC
j is given by c
o
N
j
k, and the annual inventory holding cost at MC
j is given by ∑
i
∈
I
∑
k
∈
K
c
h
(
λ
(
d
i
+
q
i
)
Y
i
j
k
/
2
N
j
k
+
z
α
a
∑
i
∈
I
Y
i
j
k
σ
i
2
).(3) Transportation Cost. The annual cost of the dispatching vehicle at MC
j is given by c
v
N
j
k, the annual total distribution costs from every MC to each retailer is given by λ
c
m
∑
i
∈
I
∑
j
∈
J
∑
g
∈
I
∑
h
∈
I
s
g
h
d
i
Z
g
h
j
k, and the annual transportation cost from plant to MC
j is given by c
j
∑
i
∈
I
∑
k
∈
K
e
j
k.(4) Returns Processing Cost. The annual repackaging cost of returned goods at MC
j is given by λ
c
b
∑
i
∈
I
q
i.The objective is to minimize the total cost of the system; we formulate the modelP as follows:(1)
P
:
min
⁡
Z
=
∑
j
∈
J
f
j
+
f
v
X
j
+
∑
j
∈
J
∑
i
∈
I
∑
k
∈
K
c
j
e
j
k
Y
i
j
k
+
∑
j
∈
J
∑
i
∈
I
∑
k
∈
K
λ
c
b
q
i
Y
i
j
k
+
λ
∑
j
∈
J
∑
k
∈
K
∑
g
∈
I
∑
h
∈
I
c
m
s
g
h
d
i
-
q
i
Z
g
h
j
k
+
∑
j
∈
J
∑
k
∈
K
c
o
+
c
v
N
j
k
+
∑
i
∈
I
∑
j
∈
J
∑
k
∈
K
c
h
+
z
α
a
∑
i
∈
I
Y
i
j
k
σ
i
2
λ
d
i
+
q
i
Y
i
j
k
2
N
j
k
+
z
α
a
∑
i
∈
I
Y
i
j
k
σ
i
2
(2)
s
.
t
.
Y
i
j
k
≤
X
j
,
∀
j
∈
J
,
∀
k
∈
K
,
∀
i
∈
I
,
(3)
hhhhh
λ
∑
i
∈
I
(
d
i
+
q
i
)
Y
i
j
k
N
j
k
≤
b
,
∀
k
∈
K
,
∀
j
∈
J
,
(4)
hhhhh
∑
k
∈
K
∑
j
∈
J
e
j
k
Y
i
j
k
N
j
k
≤
w
j
,
∀
j
∈
J
,
(5)
hhhhh
∑
h
∈
I
Z
g
h
j
k
≤
Y
g
j
k
,
∀
g
∈
Ι
,
∀
j
∈
J
,
∀
k
∈
K
,
(6)
hhhhh
∑
g
∈
I
Z
g
h
j
k
≤
Y
h
j
k
,
∀
h
∈
I
,
∀
j
∈
J
,
∀
k
∈
K
,
(7)
hhhhh
∑
j
∈
I
∑
k
∈
K
Y
i
j
k
=
1
,
∀
i
∈
I
,
∀
j
∈
J
,
(8)
hh
∑
j
∈
X
j
=
1
,
j
∈
J
,
Y
i
j
k
≤
1
,
∀
i
∈
I
,
∀
k
∈
K
,
(9)
hhhhh
X
j
∈
0,1
,
j
∈
J
,
(10)
hhhhh
Y
i
j
k
∈
0,1
,
∀
i
∈
I
,
j
∈
J
,
∀
k
∈
K
,
(11)
hhhhh
Z
g
h
j
k
∈
0,1
,
∀
g
∈
Ι
,
∀
h
∈
I
,
j
∈
J
,
∀
k
∈
K
.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.


## 3. Solution Approach

The abstract idea of solution approach is described as follows. Firstly, we derive optimal order timesN
j
k
* and optimal order size Q
j
k
* from Kuhn-Tucker condition (K-T conditions), which also rely on the decision variable Y
i
j
k. Secondly, we design a hybrid algorithm based on pseudo-parallel genetic algorithm (PPGA) and simulated annealing (SA) algorithm to obtain the optimal X
j
*, Y
i
j
k
*, and Z
g
h
j
k
*.


### 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:(12)
g
N
j
k
=
b
-
λ
∑
i
∈
I
d
i
+
q
i
Y
i
j
k
N
j
k
≥
0
,
∀
k
∈
K
,
∀
j
∈
J
,
(13)
h
N
j
k
=
w
j
-
∑
k
∈
K
∑
i
∈
I
e
j
k
Y
i
j
k
N
j
k
≥
0
,
∀
j
∈
J
.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):(14)
∇
z
N
j
k
=
∑
j
∈
I
∑
k
∈
K
d
-
j
k
-
λ
c
h
∑
i
∈
I
d
i
+
q
i
Y
i
j
k
2
N
j
k
2
,
∇
g
N
j
k
=
λ
∑
i
∈
I
d
i
+
q
i
Y
i
j
k
N
j
k
2
,
∀
k
∈
K
,
∀
j
∈
J
,
∇
h
N
j
k
=
∑
k
∈
K
∑
i
∈
I
e
j
k
Y
i
j
k
N
j
k
2
,
∀
j
∈
J
,where d
-
j
k
=
c
v
+
c
o.We introduce the generalized Lagrange multiplierγ
j
k
* and λ
j
k
* into the constraints (12) and (13), respectively. Let N
j
k
* be the K-T point; the K-T conditions are described as follows:(15)
∑
j
∈
I
∑
k
∈
K
d
-
j
k
-
c
h
∑
i
∈
I
λ
u
i
Y
i
j
k
2
N
j
k
*
2
-
γ
j
k
*
∑
i
∈
I
λ
u
i
Y
i
j
k
N
j
k
*
2
-
λ
j
k
*
∑
k
∈
K
∑
i
∈
I
e
j
k
Y
i
j
k
N
j
k
*
2
=
0
,
γ
j
k
*
b
-
∑
i
∈
I
λ
u
i
Y
i
j
k
N
j
k
*
=
0
,
λ
j
k
*
w
j
-
∑
k
∈
K
∑
i
∈
I
e
j
k
Y
i
j
k
N
j
k
*
=
0
,
γ
j
k
*
,
λ
j
k
*
≥
0
,where u
i
=
d
i
+
q
i.Letγ
j
k
*
=
λ
j
k
*
=
0; then N
j
k
*
=
c
h
∑
i
∈
I
λ
u
i
Y
i
j
k
/
2
d
-
j
k. Since the original problem is convex programming, the cost of the whole system reaches the minimum value when N
j
k
=
N
j
k
* and the optimal order size is Q
j
k
*
=
∑
i
∈
I
e
j
k
Y
i
j
k
/
N
j
k
*
=
2
d
-
j
k
∑
i
∈
I
e
i
j
2
Y
i
j
k
/
λ
c
h
u
i.PuttingN
j
k
*
=
c
h
∑
i
∈
I
λ
u
i
Y
i
j
k
/
2
d
-
j
k into the objective function (1), we can get the model P
N as follows:(16)
P
N
:
min
⁡
Z
¯
=
∑
j
∈
J
+
c
b
∑
i
∈
I
q
i
Y
i
j
k
+
c
-
h
∑
i
∈
I
σ
-
i
2
Y
i
j
k
f
-
X
j
hhhhhhhhhhhh
+
∑
k
∈
K
∑
i
∈
I
σ
-
i
2
Y
i
j
k
d
-
-
j
k
∑
i
∈
I
e
j
k
Y
i
j
k
+
c
-
j
∑
i
∈
I
u
i
Y
i
j
k
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
c
b
∑
i
∈
I
q
i
Y
i
j
k
+
c
-
h
∑
i
∈
I
σ
-
i
2
Y
i
j
k
s
.
t
.
2


~


11
,where d
-
-
j
k
=
2
c
h
d
-
j
k
=
2
c
h
c
v
+
∑
g
∈
I
∑
h
∈
I
λ
c
m
(
d
i
-
q
i
)
+
c
r
q
i
s
g
h
Z
g
h
j
k
+
c
o, f
-
=
f
j
+
f
v, c
-
j
=
c
j, c
-
h
=
λ
c
h
z
α, and σ
-
i
2
=
a
σ
i
2.


### 3.2. 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 pseudo-parallel genetic algorithm integrating simulated annealing (PPGASA).


#### 3.2.1. 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 unrepeatedV
+
S natural numbers constitute a sequence, which represents an individual. While 1,2
,
…
,
V indicate candidate MCs and V
+
1
,
…
,
V
+
S 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(17)
f
k
=
1
Z
¯
.(3) Selection. 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 N; the fitness value of the individual k is f
k and a number ξ
∈
[
0,1
] is generated randomly.If∑
k
=
1
i
-
1
f
k
/
∑
k
=
1
N
f
k
<
ξ
≤
∑
k
=
1
i
f
k
/
∑
k
=
1
N
f
k, then individual i is selected to be replicated.(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.


#### 3.2.2. 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. AfterN-generation’s independent evolution, the dual-channel exchange mechanism is invoked as follows.Step 1.
Select the individualA with the minimal fitness function value from the unsuccessful subgroup.Step 2.
Select randomly individualB from the successful subgroup.Step 3.
Iff
B
<
f
A, A enter successful subgroup and B enter unsuccessful subgroup; otherwise, remain unchanged.Step 4.
Select the individualC with the maximal fitness function value from the successful subgroup.Step 5.
Select randomly individualD from the unsuccessful subgroup.Step 6.
Iff
D
>
f
C, D enter unsuccessful subgroup and C enter successful subgroup; otherwise, remain unchanged.


#### 3.2.3. 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 p
=
exp
⁡
-
1
/
f
new
-
1
/
f
old
/
T > random digit, where T is annealing temperature, f
new is the objective function value of the new solution, and f
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 t
k
+
1
=
α
t
k, k
≥
0, 0
<
α
<
1; the nearer α is to 1, the slower the temperature decreases.


#### 3.2.4. Termination

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


### 3.3. Algorithm Flow

Step 0. Get the formulas for solving N
j
k
* and Q
j
k
*.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 g; the population size N; evolution terminate generation M; crossover probability p
c; mutation probability p
m; temperature of the cooling coefficient α; the initial annealing temperature T; and so on.Step 2. Calculate the fitness value f
k of an individual. If the parent optimal solution and offspring optimal solution are equal during continuous M 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 N generations, the dual channel exchange mechanism is executed.Step 4. If f
i
<
f
j
(
j
>
i
), accept the new individual; otherwise, accept the new individuals with the probability p
=
exp
⁡
(
-
1
/
f
new
-
1
/
f
old
/
T
).Step 5. Update the annealing temperature, and return to Step 2.The pseudocodes of HGSAA are shown in Pseudocode1.Pseudocode 1:Pseudocode of the proposed PPGASA.
Procedure: PPGASA for LIRP
Input: coordinates of nodes, demands and returns of retailers, MC parameters, vehicle capacity, PPGASA parameters
Output: the best solution (include routes, MCs locations, order times and order size)
Begin
Take pop
for i = 1 to pop_size
Calculate individual fitness value  f
k
end
f
best
=
m
a
x
{
f
k
}
n = 1
k = 1
T = 100
while n
≤ 300
select operator
if random ≤ p
c  then
crossover operator
end
if random ≤ p
m  then
mutation operator
end
newpop
for m = 1to newpop_size
if objvalue(m) < objvalue1(m)
Sel(m,:) = Z(m,:)
else
Sel(m,:) = Cos(m,:)
end
end
pops1 = zhuanh(Sel)
[objvalue1] = calobjvalue(pops1, d)
for i = 1: popsize
if objvalue(i) > objvalue1(i)
objvalue(i) = objvalue1(i)
Z(i,:) = Sel(i,:)
else
p = rand
if p ≤ exp((objvalue1(i) − objvalue(i))/T)
objvalue(i) = objvalue1(i)
Z(i,:) = Sel(i,:)
end
end
end
gen = gen + 1
pops2 = zhuanh(Z)
[objvalue2] = calobjvalue(pops2, d)
[hh, kk] = min(objvalue2)
GY(gen, 1) = hh
GX(gen,:) = pops2(kk,:)
end
T = α * T
end
output the best solution
end


## 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:(12)
g
N
j
k
=
b
-
λ
∑
i
∈
I
d
i
+
q
i
Y
i
j
k
N
j
k
≥
0
,
∀
k
∈
K
,
∀
j
∈
J
,
(13)
h
N
j
k
=
w
j
-
∑
k
∈
K
∑
i
∈
I
e
j
k
Y
i
j
k
N
j
k
≥
0
,
∀
j
∈
J
.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):(14)
∇
z
N
j
k
=
∑
j
∈
I
∑
k
∈
K
d
-
j
k
-
λ
c
h
∑
i
∈
I
d
i
+
q
i
Y
i
j
k
2
N
j
k
2
,
∇
g
N
j
k
=
λ
∑
i
∈
I
d
i
+
q
i
Y
i
j
k
N
j
k
2
,
∀
k
∈
K
,
∀
j
∈
J
,
∇
h
N
j
k
=
∑
k
∈
K
∑
i
∈
I
e
j
k
Y
i
j
k
N
j
k
2
,
∀
j
∈
J
,where d
-
j
k
=
c
v
+
c
o.We introduce the generalized Lagrange multiplierγ
j
k
* and λ
j
k
* into the constraints (12) and (13), respectively. Let N
j
k
* be the K-T point; the K-T conditions are described as follows:(15)
∑
j
∈
I
∑
k
∈
K
d
-
j
k
-
c
h
∑
i
∈
I
λ
u
i
Y
i
j
k
2
N
j
k
*
2
-
γ
j
k
*
∑
i
∈
I
λ
u
i
Y
i
j
k
N
j
k
*
2
-
λ
j
k
*
∑
k
∈
K
∑
i
∈
I
e
j
k
Y
i
j
k
N
j
k
*
2
=
0
,
γ
j
k
*
b
-
∑
i
∈
I
λ
u
i
Y
i
j
k
N
j
k
*
=
0
,
λ
j
k
*
w
j
-
∑
k
∈
K
∑
i
∈
I
e
j
k
Y
i
j
k
N
j
k
*
=
0
,
γ
j
k
*
,
λ
j
k
*
≥
0
,where u
i
=
d
i
+
q
i.Letγ
j
k
*
=
λ
j
k
*
=
0; then N
j
k
*
=
c
h
∑
i
∈
I
λ
u
i
Y
i
j
k
/
2
d
-
j
k. Since the original problem is convex programming, the cost of the whole system reaches the minimum value when N
j
k
=
N
j
k
* and the optimal order size is Q
j
k
*
=
∑
i
∈
I
e
j
k
Y
i
j
k
/
N
j
k
*
=
2
d
-
j
k
∑
i
∈
I
e
i
j
2
Y
i
j
k
/
λ
c
h
u
i.PuttingN
j
k
*
=
c
h
∑
i
∈
I
λ
u
i
Y
i
j
k
/
2
d
-
j
k into the objective function (1), we can get the model P
N as follows:(16)
P
N
:
min
⁡
Z
¯
=
∑
j
∈
J
+
c
b
∑
i
∈
I
q
i
Y
i
j
k
+
c
-
h
∑
i
∈
I
σ
-
i
2
Y
i
j
k
f
-
X
j
hhhhhhhhhhhh
+
∑
k
∈
K
∑
i
∈
I
σ
-
i
2
Y
i
j
k
d
-
-
j
k
∑
i
∈
I
e
j
k
Y
i
j
k
+
c
-
j
∑
i
∈
I
u
i
Y
i
j
k
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
c
b
∑
i
∈
I
q
i
Y
i
j
k
+
c
-
h
∑
i
∈
I
σ
-
i
2
Y
i
j
k
s
.
t
.
2


~


11
,where d
-
-
j
k
=
2
c
h
d
-
j
k
=
2
c
h
c
v
+
∑
g
∈
I
∑
h
∈
I
λ
c
m
(
d
i
-
q
i
)
+
c
r
q
i
s
g
h
Z
g
h
j
k
+
c
o, f
-
=
f
j
+
f
v, c
-
j
=
c
j, c
-
h
=
λ
c
h
z
α, and σ
-
i
2
=
a
σ
i
2.


## 3.2. 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 pseudo-parallel genetic algorithm integrating simulated annealing (PPGASA).


### 3.2.1. 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 unrepeatedV
+
S natural numbers constitute a sequence, which represents an individual. While 1,2
,
…
,
V indicate candidate MCs and V
+
1
,
…
,
V
+
S 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(17)
f
k
=
1
Z
¯
.(3) Selection. 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 N; the fitness value of the individual k is f
k and a number ξ
∈
[
0,1
] is generated randomly.If∑
k
=
1
i
-
1
f
k
/
∑
k
=
1
N
f
k
<
ξ
≤
∑
k
=
1
i
f
k
/
∑
k
=
1
N
f
k, then individual i is selected to be replicated.(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.


### 3.2.2. 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. AfterN-generation’s independent evolution, the dual-channel exchange mechanism is invoked as follows.Step 1.
Select the individualA with the minimal fitness function value from the unsuccessful subgroup.Step 2.
Select randomly individualB from the successful subgroup.Step 3.
Iff
B
<
f
A, A enter successful subgroup and B enter unsuccessful subgroup; otherwise, remain unchanged.Step 4.
Select the individualC with the maximal fitness function value from the successful subgroup.Step 5.
Select randomly individualD from the unsuccessful subgroup.Step 6.
Iff
D
>
f
C, D enter unsuccessful subgroup and C enter successful subgroup; otherwise, remain unchanged.


### 3.2.3. 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 p
=
exp
⁡
-
1
/
f
new
-
1
/
f
old
/
T > random digit, where T is annealing temperature, f
new is the objective function value of the new solution, and f
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 t
k
+
1
=
α
t
k, k
≥
0, 0
<
α
<
1; the nearer α is to 1, the slower the temperature decreases.


### 3.2.4. Termination

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


## 3.2.1. 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 unrepeatedV
+
S natural numbers constitute a sequence, which represents an individual. While 1,2
,
…
,
V indicate candidate MCs and V
+
1
,
…
,
V
+
S 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(17)
f
k
=
1
Z
¯
.(3) Selection. 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 N; the fitness value of the individual k is f
k and a number ξ
∈
[
0,1
] is generated randomly.If∑
k
=
1
i
-
1
f
k
/
∑
k
=
1
N
f
k
<
ξ
≤
∑
k
=
1
i
f
k
/
∑
k
=
1
N
f
k, then individual i is selected to be replicated.(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.


## 3.2.2. 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. AfterN-generation’s independent evolution, the dual-channel exchange mechanism is invoked as follows.Step 1.
Select the individualA with the minimal fitness function value from the unsuccessful subgroup.Step 2.
Select randomly individualB from the successful subgroup.Step 3.
Iff
B
<
f
A, A enter successful subgroup and B enter unsuccessful subgroup; otherwise, remain unchanged.Step 4.
Select the individualC with the maximal fitness function value from the successful subgroup.Step 5.
Select randomly individualD from the unsuccessful subgroup.Step 6.
Iff
D
>
f
C, D enter unsuccessful subgroup and C enter successful subgroup; otherwise, remain unchanged.


## 3.2.3. 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 p
=
exp
⁡
-
1
/
f
new
-
1
/
f
old
/
T > random digit, where T is annealing temperature, f
new is the objective function value of the new solution, and f
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 t
k
+
1
=
α
t
k, k
≥
0, 0
<
α
<
1; the nearer α is to 1, the slower the temperature decreases.


## 3.2.4. Termination

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


## 3.3. Algorithm Flow

Step 0. Get the formulas for solving N
j
k
* and Q
j
k
*.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 g; the population size N; evolution terminate generation M; crossover probability p
c; mutation probability p
m; temperature of the cooling coefficient α; the initial annealing temperature T; and so on.Step 2. Calculate the fitness value f
k of an individual. If the parent optimal solution and offspring optimal solution are equal during continuous M 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 N generations, the dual channel exchange mechanism is executed.Step 4. If f
i
<
f
j
(
j
>
i
), accept the new individual; otherwise, accept the new individuals with the probability p
=
exp
⁡
(
-
1
/
f
new
-
1
/
f
old
/
T
).Step 5. Update the annealing temperature, and return to Step 2.The pseudocodes of HGSAA are shown in Pseudocode1.Pseudocode 1:Pseudocode of the proposed PPGASA.
Procedure: PPGASA for LIRP
Input: coordinates of nodes, demands and returns of retailers, MC parameters, vehicle capacity, PPGASA parameters
Output: the best solution (include routes, MCs locations, order times and order size)
Begin
Take pop
for i = 1 to pop_size
Calculate individual fitness value  f
k
end
f
best
=
m
a
x
{
f
k
}
n = 1
k = 1
T = 100
while n
≤ 300
select operator
if random ≤ p
c  then
crossover operator
end
if random ≤ p
m  then
mutation operator
end
newpop
for m = 1to newpop_size
if objvalue(m) < objvalue1(m)
Sel(m,:) = Z(m,:)
else
Sel(m,:) = Cos(m,:)
end
end
pops1 = zhuanh(Sel)
[objvalue1] = calobjvalue(pops1, d)
for i = 1: popsize
if objvalue(i) > objvalue1(i)
objvalue(i) = objvalue1(i)
Z(i,:) = Sel(i,:)
else
p = rand
if p ≤ exp((objvalue1(i) − objvalue(i))/T)
objvalue(i) = objvalue1(i)
Z(i,:) = Sel(i,:)
end
end
end
gen = gen + 1
pops2 = zhuanh(Z)
[objvalue2] = calobjvalue(pops2, d)
[hh, kk] = min(objvalue2)
GY(gen, 1) = hh
GX(gen,:) = pops2(kk,:)
end
T = α * T
end
output the best solution
end


## 4. Computational Experiments and Results Analysis

An example is used to illustrate the proposed heuristic method. The data ofGaskell
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 N
(
20,1
). The other parameters are as follows: the inventory holding cost per unit of goods per year h
=
6, the vehicle capacity b
=
500, the delivering cost per unit distance c
m
=
1, the handling cost per unit product c
b
=
3, fixed cost of dispatching vehicles per time c
v
=
18, and the conversion constant λ
=
300. The customers’ daily returns are uniformly generated from N
(
3,1
). c
o, f
v are uniformly generated from U
[
6,10
]
,
U
[
30,70
].


### 4.1. Parameters Discussion

Based on Matlab 7.0 platform, we programmed the PPGASA and then run it 30 times on a computer (CPU: Intel Core i3-2130 @ 3.4 GHz, RAM: 2.85 GB DDR, OS: Windows 7).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 Tables1–5.Table 1
The results with differentp
m (cost: million yuan, CPU time: second).
p
m
Cost
CPU time
Range
Mean
Standard deviation
Coefficient of variation
Range(second)
Mean(second)
Standard deviation
Coefficient of variation
0.0001
1.28
2.17
0.34
0.27
0.76
1.46
0.17
0.28
0.0005
0.99
1.65
0.27
0.16
0.86
1.60
0.25
0.15
0.001
1.60
2.15
0.52
0.33
1.34
1.54
0.21
0.30
0.005
1.30
2.1
0.40
0.30
0.86
1.53
0.17
0.27
n = 20, p
c = 0.9, q = 0.9, and T = 100.Table 2
The results with differentp
c (cost: million yuan, CPU time: second).
p
c
Cost
CPU time
Range
Mean
Standard deviation
Coefficient of variation
Range
Mean
Standard deviation
Coefficient of variation
0.1
2.70
1.97
0.65
0.33
4.26
2.92
1.13
0.39
0.4
2.26
2.18
0.65
0.30
2.01
1.42
0.57
0.40
0.7
2.91
2.17
0.84
0.39
1.34
1.42
0.42
0.30
0.9
0.99
1.65
0.27
0.16
0.86
1.60
0.25
0.15
n = 20, p
m = 0.0005, q = 0.9, and T = 100.Table 3
The results with differentn (cost: million yuan, CPU time: second).
n
Cost
CPU time
Range
Mean
Standarddeviation
Coefficient of variation
Range
Mean
Standarddeviation
Coefficient of variation
10
1.23
1.70
0.52
0.31
0.43
0.75
0.12
0.16
15
1.21
1.72
0.46
0.27
0.39
0.39
0.16
0.14
20
0.99
1.65
0.27
0.16
0.86
1.60
0.25
0.15
30
1.18
1.86
0.37
0.20
2.53
2.55
0.67
0.26
p
m = 0.0005, p
c = 0.9, q = 0.9, and T = 100.Table 4
The results with differentT (cost: million yuan, CPU time: second).
T
Cost
CPU time
Range
Mean
Standarddeviation
Coefficient of variation
Range
Mean
Standarddeviation
Coefficient of variation
10
1.96
2.12
0.72
0.34
0.72
1.50
0.42
0.28
40
2.05
2.09
0.61
0.29
0.78
1.59
0.47
0.29
70
1.10
1.98
0.58
0.29
0.92
1.54
0.39
0.27
100
0.99
1.65
0.27
0.16
0.86
1.60
0.25
0.15
p
m = 0.0005, p
c = 0.9, q = 0.9, and n = 20.Table 5
The results with differentq (cost: million yuan, CPU time: second).
q
Cost
CPU time
Range
Mean
Standarddeviation
Coefficient of variation
Range
Mean
Standarddeviation
Coefficient of variation
0.1
1.56
1.98
0.68
0.35
0.77
1.57
0.43
0.27
0.4
1.68
1.93
0.61
0.32
1.53
1.58
0.38
0.24
0.7
1.60
2.13
0.73
0.34
0.59
1.51
0.44
0.29
0.9
0.99
1.65
0.27
0.16
0.86
1.60
0.25
0.15
p
m = 0.0005, p
c = 0.9, T = 100, and n = 20.Figures2–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 p
c = 0.7~0.9, p
m = 0.0005~0.001,q = 0.85~0.9,T = 100, and n
=
20.Figure 2
Optimization results on different valuesp
m.
(a)
(b)Figure 3
Optimization results on different valuesp
c.
(a)
(b)Figure 4
Optimization results on different valuesn.
(a)
(b)Figure 5
Optimization results on different valuesT.
(a)
(b)Figure 6
Optimization results on different valuesq.
(a)
(b)


### 4.2. 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. Table6 shows the solution. Figure 7 shows topological structure of the supply chain. MCs were established at j
1, j
2, j
3, j
4, j
5 points and there are eight vehicles distribution routes.Table 6
Solution of instanceGaskell 67 − 29 × 5.
MC
Routing number
Routing
Order times
Order interval(day)
Order quantity
Safety stock
j
1
k
1
i
6 → i
21 → i
1
0
32
9
656
630
j
2
k
2
i
17 → i
9 → i
1
2 → i
18 → i
23 → i
19 → i
14
31
10
590
610
j
3
k
3
i
2
53
6
702
744
k
4
i
25 → i
5
31
10
387
400
j
4
k
5
i
3 → i
4 → i
1 → i
24 → i
29 → i
26
46
7
613
658
k
6
i
28 → i
15
35
9
531
558
j
5
k
7
i
16 → i
13 → i
7 → i
11 → i
8
40
8
555
592
k
8
i
20 → i
22
35
9
253
144Figure 7
Topological structure of the network.For comparison, GA is programmed by Matlab 7.0 as well, and the instanceGaskel
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.Figure 8
Trends of optimal objective function value by GA.Figure 9
Trends of optimal objective function value by PPGASA.The fluctuation cure of optimal objective function values obtained from the 100 times experiments is shown in Figures10 and 11, respectively.Figure 10
The fluctuation curve of optimal objective function value by PPGASA.Figure 11
The fluctuation curve of optimal objective function value by GA.The optimal objective function value and the CPU time of these two algorithms are shown in Table7. As shown in this table, the cost of PPGASA is significantly lower than GA (
p
≪
0.05
), and the efficiency of algorithms measured by CPU time in seconds is significant difference between the two algorithms (
p
≪
0.05
).Table 7
Statistical results of optimal objective function value of two algorithms (cost: yuan, CPU time: second).
Algorithm
Mean
Standard deviation
Coefficient of variation
Significance test
t value
p value (Sig. 2-tailed)
Cost
PPGASA
1650568
271320
0.16
−9.144
0.000
GA
2134686
454600
0.21
CPU time
PPGASA
1.60
0.25
0.15
6.826
0.000
GA
1.38
0.25
0.18


### 4.3. Extended Experiments

In this section, a series of experiments are given to show that PPGASA is more efficient and stable than GA. Similar with Section4.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 ofGaskell
67
-
22
×
5.Each instance was calculated 100 times by PPGASA and GA, respectively; the results are shown in Tables8 and 9.Table 8
Optimal objective function values of two algorithms (yuan).
Instance name
Algorithm
Mean
Standard deviation
Coefficient of variation
Significance test
t value
p value (Sig. 2-tailed)
Perl 183 – 12 × 2
PPGASA
27274.15
3072.70
0.1127
−15.71
0.000
GA
59644.00
20373.95
0.3416
Gaskell 67 – 22 × 5
PPGASA
176597.12
51842.11
0.2936
−2.069
0.040
GA
193606.40
63820.20
0.3296
Gaskell 67 – 36 × 5
PPGASA
1532353.00
240454.91
0.1569
−24.268
0.000
GA
3266300.00
672809.81
0.2060
Perl 183 – 55 × 15
PPGASA
4368430.00
336908.31
0.0771
−3.877
0.000
GA
4571132.00
399811.55
0.0874
Christofides 69 − 75 × 10
PPGASA
5653803.00
481169.53
0.0851
−2.113
0.036
GA
58256800.00
655839.90
0.0113
Perl 183 – 85 × 7
PPGASA
6095438.00
219212.14
0.0360
−3.388
0.001
GA
6221775.00
301680.88
0.0485
Christofides 69 – 100 × 10
PPGASA
6848350.00
617689.91
0.0902
−3.391
0.001
GA
7179759.00
757338.00
0.1055Table 9
CPU time of two algorithms (seconds).
Instance name
Algorithm
Mean
Standard deviation
Coefficient of variation
significance test
t value
p value (Sig. 2-tailed)
Perl 183 − 12 × 2
PPGASA
0.5325
0.0931
0.1748
−2.488
0.014
GA
0.5712
0.1243
0.2176
Gaskell 67 − 22 × 5
PPGASA
1.3804
0.2670
0.1934
−3.158
0.002
GA
1.5322
0.3996
0.2608
Gaskell 67 − 36 × 5
PPGASA
1.8469
0.2226
0.1205
−4.704
0.000
GA
2.0334
0.3282
0.1614
Perl 183 − 55 × 15
PPGASA
3.8288
0.2372
0.0620
−5.766
0.000
GA
4.1037
0.4135
0.1008
Christofides 69 − 75 × 10
PPGASA
8.2635
1.6029
0.1940
−10.607
0.000
GA
11.5331
2.6330
0.2283
Perl 183 − 85 × 7
PPGASA
9.3730
0.9723
0.1037
−5.326
0.000
GA
10.4171
1.7021
0.1634
Christofides 69 − 100 × 10
PPGASA
10.3791
1.5966
0.1538
−13.805
0.000
GA
14.4457
2.4756
0.1714


### 4.4. Result Analysis

Tables8 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 Tables8 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 than GA since the coefficient of PPGASA is less than GA. At the same time, all the p values (Sig. 2-tailed) are less than 0.05; it means that it is difference with remarkably significance between PPGASA and GA. According to all of the above, PPGASA is more effective than GA in resolving such problems.


## 4.1. Parameters Discussion

Based on Matlab 7.0 platform, we programmed the PPGASA and then run it 30 times on a computer (CPU: Intel Core i3-2130 @ 3.4 GHz, RAM: 2.85 GB DDR, OS: Windows 7).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 Tables1–5.Table 1
The results with differentp
m (cost: million yuan, CPU time: second).
p
m
Cost
CPU time
Range
Mean
Standard deviation
Coefficient of variation
Range(second)
Mean(second)
Standard deviation
Coefficient of variation
0.0001
1.28
2.17
0.34
0.27
0.76
1.46
0.17
0.28
0.0005
0.99
1.65
0.27
0.16
0.86
1.60
0.25
0.15
0.001
1.60
2.15
0.52
0.33
1.34
1.54
0.21
0.30
0.005
1.30
2.1
0.40
0.30
0.86
1.53
0.17
0.27
n = 20, p
c = 0.9, q = 0.9, and T = 100.Table 2
The results with differentp
c (cost: million yuan, CPU time: second).
p
c
Cost
CPU time
Range
Mean
Standard deviation
Coefficient of variation
Range
Mean
Standard deviation
Coefficient of variation
0.1
2.70
1.97
0.65
0.33
4.26
2.92
1.13
0.39
0.4
2.26
2.18
0.65
0.30
2.01
1.42
0.57
0.40
0.7
2.91
2.17
0.84
0.39
1.34
1.42
0.42
0.30
0.9
0.99
1.65
0.27
0.16
0.86
1.60
0.25
0.15
n = 20, p
m = 0.0005, q = 0.9, and T = 100.Table 3
The results with differentn (cost: million yuan, CPU time: second).
n
Cost
CPU time
Range
Mean
Standarddeviation
Coefficient of variation
Range
Mean
Standarddeviation
Coefficient of variation
10
1.23
1.70
0.52
0.31
0.43
0.75
0.12
0.16
15
1.21
1.72
0.46
0.27
0.39
0.39
0.16
0.14
20
0.99
1.65
0.27
0.16
0.86
1.60
0.25
0.15
30
1.18
1.86
0.37
0.20
2.53
2.55
0.67
0.26
p
m = 0.0005, p
c = 0.9, q = 0.9, and T = 100.Table 4
The results with differentT (cost: million yuan, CPU time: second).
T
Cost
CPU time
Range
Mean
Standarddeviation
Coefficient of variation
Range
Mean
Standarddeviation
Coefficient of variation
10
1.96
2.12
0.72
0.34
0.72
1.50
0.42
0.28
40
2.05
2.09
0.61
0.29
0.78
1.59
0.47
0.29
70
1.10
1.98
0.58
0.29
0.92
1.54
0.39
0.27
100
0.99
1.65
0.27
0.16
0.86
1.60
0.25
0.15
p
m = 0.0005, p
c = 0.9, q = 0.9, and n = 20.Table 5
The results with differentq (cost: million yuan, CPU time: second).
q
Cost
CPU time
Range
Mean
Standarddeviation
Coefficient of variation
Range
Mean
Standarddeviation
Coefficient of variation
0.1
1.56
1.98
0.68
0.35
0.77
1.57
0.43
0.27
0.4
1.68
1.93
0.61
0.32
1.53
1.58
0.38
0.24
0.7
1.60
2.13
0.73
0.34
0.59
1.51
0.44
0.29
0.9
0.99
1.65
0.27
0.16
0.86
1.60
0.25
0.15
p
m = 0.0005, p
c = 0.9, T = 100, and n = 20.Figures2–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 p
c = 0.7~0.9, p
m = 0.0005~0.001,q = 0.85~0.9,T = 100, and n
=
20.Figure 2
Optimization results on different valuesp
m.
(a)
(b)Figure 3
Optimization results on different valuesp
c.
(a)
(b)Figure 4
Optimization results on different valuesn.
(a)
(b)Figure 5
Optimization results on different valuesT.
(a)
(b)Figure 6
Optimization results on different valuesq.
(a)
(b)


## 4.2. 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. Table6 shows the solution. Figure 7 shows topological structure of the supply chain. MCs were established at j
1, j
2, j
3, j
4, j
5 points and there are eight vehicles distribution routes.Table 6
Solution of instanceGaskell 67 − 29 × 5.
MC
Routing number
Routing
Order times
Order interval(day)
Order quantity
Safety stock
j
1
k
1
i
6 → i
21 → i
1
0
32
9
656
630
j
2
k
2
i
17 → i
9 → i
1
2 → i
18 → i
23 → i
19 → i
14
31
10
590
610
j
3
k
3
i
2
53
6
702
744
k
4
i
25 → i
5
31
10
387
400
j
4
k
5
i
3 → i
4 → i
1 → i
24 → i
29 → i
26
46
7
613
658
k
6
i
28 → i
15
35
9
531
558
j
5
k
7
i
16 → i
13 → i
7 → i
11 → i
8
40
8
555
592
k
8
i
20 → i
22
35
9
253
144Figure 7
Topological structure of the network.For comparison, GA is programmed by Matlab 7.0 as well, and the instanceGaskel
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.Figure 8
Trends of optimal objective function value by GA.Figure 9
Trends of optimal objective function value by PPGASA.The fluctuation cure of optimal objective function values obtained from the 100 times experiments is shown in Figures10 and 11, respectively.Figure 10
The fluctuation curve of optimal objective function value by PPGASA.Figure 11
The fluctuation curve of optimal objective function value by GA.The optimal objective function value and the CPU time of these two algorithms are shown in Table7. As shown in this table, the cost of PPGASA is significantly lower than GA (
p
≪
0.05
), and the efficiency of algorithms measured by CPU time in seconds is significant difference between the two algorithms (
p
≪
0.05
).Table 7
Statistical results of optimal objective function value of two algorithms (cost: yuan, CPU time: second).
Algorithm
Mean
Standard deviation
Coefficient of variation
Significance test
t value
p value (Sig. 2-tailed)
Cost
PPGASA
1650568
271320
0.16
−9.144
0.000
GA
2134686
454600
0.21
CPU time
PPGASA
1.60
0.25
0.15
6.826
0.000
GA
1.38
0.25
0.18


## 4.3. Extended Experiments

In this section, a series of experiments are given to show that PPGASA is more efficient and stable than GA. Similar with Section4.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 ofGaskell
67
-
22
×
5.Each instance was calculated 100 times by PPGASA and GA, respectively; the results are shown in Tables8 and 9.Table 8
Optimal objective function values of two algorithms (yuan).
Instance name
Algorithm
Mean
Standard deviation
Coefficient of variation
Significance test
t value
p value (Sig. 2-tailed)
Perl 183 – 12 × 2
PPGASA
27274.15
3072.70
0.1127
−15.71
0.000
GA
59644.00
20373.95
0.3416
Gaskell 67 – 22 × 5
PPGASA
176597.12
51842.11
0.2936
−2.069
0.040
GA
193606.40
63820.20
0.3296
Gaskell 67 – 36 × 5
PPGASA
1532353.00
240454.91
0.1569
−24.268
0.000
GA
3266300.00
672809.81
0.2060
Perl 183 – 55 × 15
PPGASA
4368430.00
336908.31
0.0771
−3.877
0.000
GA
4571132.00
399811.55
0.0874
Christofides 69 − 75 × 10
PPGASA
5653803.00
481169.53
0.0851
−2.113
0.036
GA
58256800.00
655839.90
0.0113
Perl 183 – 85 × 7
PPGASA
6095438.00
219212.14
0.0360
−3.388
0.001
GA
6221775.00
301680.88
0.0485
Christofides 69 – 100 × 10
PPGASA
6848350.00
617689.91
0.0902
−3.391
0.001
GA
7179759.00
757338.00
0.1055Table 9
CPU time of two algorithms (seconds).
Instance name
Algorithm
Mean
Standard deviation
Coefficient of variation
significance test
t value
p value (Sig. 2-tailed)
Perl 183 − 12 × 2
PPGASA
0.5325
0.0931
0.1748
−2.488
0.014
GA
0.5712
0.1243
0.2176
Gaskell 67 − 22 × 5
PPGASA
1.3804
0.2670
0.1934
−3.158
0.002
GA
1.5322
0.3996
0.2608
Gaskell 67 − 36 × 5
PPGASA
1.8469
0.2226
0.1205
−4.704
0.000
GA
2.0334
0.3282
0.1614
Perl 183 − 55 × 15
PPGASA
3.8288
0.2372
0.0620
−5.766
0.000
GA
4.1037
0.4135
0.1008
Christofides 69 − 75 × 10
PPGASA
8.2635
1.6029
0.1940
−10.607
0.000
GA
11.5331
2.6330
0.2283
Perl 183 − 85 × 7
PPGASA
9.3730
0.9723
0.1037
−5.326
0.000
GA
10.4171
1.7021
0.1634
Christofides 69 − 100 × 10
PPGASA
10.3791
1.5966
0.1538
−13.805
0.000
GA
14.4457
2.4756
0.1714


## 4.4. Result Analysis

Tables8 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 Tables8 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 than GA since the coefficient of PPGASA is less than GA. At the same time, all the p values (Sig. 2-tailed) are less than 0.05; it means that it is difference with remarkably significance between PPGASA and GA. According to all of the above, PPGASA is more effective than GA in resolving such problems.


## 5. 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 methods. A heuristic algorithm named PPGASA based on PPGA and SA is designed to solve our model.
(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.



---

*Source: 586581-2015-03-02.xml*