# Solving a Closed-Loop Location-Inventory-Routing Problem with Mixed Quality Defects Returns in E-Commerce by Hybrid Ant Colony Optimization Algorithm

Authors: Shuai Deng; Yanhui Li; Hao Guo; Bailing Liu
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

This paper presents a closed-loop location-inventory-routing problem model considering both quality defect returns and nondefect returns in e-commerce supply chain system. The objective is to minimize the total cost produced in both forward and reverse logistics networks. We propose a combined optimization algorithm named hybrid ant colony optimization algorithm (HACO) to address this model that is an NP-hard problem. Our experimental results show that the proposed HACO is considerably efficient and effective in solving this model.


---

## Body

## 1. Introduction

According to eMarketer, worldwide business-to-consumer (B2C) e-commerce sales reached $1.471 trillion in 2014, increasing by nearly 20% over 2013 [1]. Customers have grown accustomed to return unwanted products back to the store for any reasons. It is reported that the proportion of customer returns online range from 18% to 74% of original orders under e-commerce environment [2, 3]; thus enterprises take various measures to prevent the appearance of quality defects. However, quality defect is inevitable. So, it is necessary for us to take into account both quality defect returns and nondefect returns; we call it mixed quality defect returns (MQDR), when considering the closed-loop supply chain as a support system in e-commerce environment.As a classic discrete dynamics problem, the customer service level is determined by three important decisions: facility location decision, inventory decision, and transportation decision [4]. Obviously, facility location, inventory control, and transportation optimization are highly related. For example, delivery in small lots and high frequency leads to reducing the in-inventory cost but increases the additional transportation cost. In addition, facility location decision needs to consider inventory decision and distribution decision. Perl and Sirisoponsilp [5] discuss the interdependence between the three key elements. Ballou and Masters [6] provide a schematic representation of the interrelationships among facility location, inventory control, and transportation optimization.In the literature, many papers studied the integration and coordination of any two of the above three decisions: location-inventory problem (LIP), location-routing problem (LRP), and inventory-routing problem (IRP). For reviews on LIP, readers can refer to Erlebacher and Meller [7], Daskin et al. [8, 9], and Liao et al. [10]. For LRP, please refer to Balakrishnan et al. [11], Min et al. [12], and Nagy and Salhi [13]. Refer to Chan et al. [14], Kleywegt et al. [15], and Adelman [16] for IRP.There are few researches about the integration optimization of location-inventory-routing problem (LIRP). Some researchers attempt to carry out research on LIRP [17]. Liu and Lee [18] firstly studied this interesting problem; they proposed a two-phase heuristic method to solve the multidepot location-routing problem (MDLRP) considering inventory optimization. In order to avoid the local optimal solution, Liu and Lin [19] designed a global optimizing heuristic method to find the solutions for LIRP. Shen and Qi [20] presented an algorithm based on Lagrangian relaxation to minimize the inventory and routing costs in strategic location models. They focused on the layout phase and used continuous approximation to get the approximate optimal routing cost, but the vehicle routing was not optimized in their models. Javid and Azad [21] presented a novel LIRP model and proposed heuristic method containing two stages: constructive stage and improvement stage. Ahmadi-Javid and Seddighi [22] presented a mixed-integer programming model and a three-phase heuristic to solve the LIRP with multisource distribution logistics network. Guerrero et al. [23] researched the LIRP with deterministic demand and provided the hybrid algorithm to solve the problem. Zhang et al. [24] proposed a hybrid metaheuristic solution to LIRP considering multiple depots and geographically dispersed customers. Nekooghadirli et al. [25] presented a novel biobjective model of LIRP model considering a multiperiod and multiproduct system. Based on Lagrangian relaxation and a column generation technique, Guerrero et al. [26] developed a relax-and-price heuristic to solve ILRP; they proposed two dependent constraint sets with an exponential nature: Lagrangian relaxation and a column generation technique.However, little research has been conducted on the LIRP considering returns. Li et al. [27] presented the HGSAA algorithm to solve a LIRP model considering returns under e-supply chain environment. To be more consistent with reality, Liu et al. [28] introduced a stochastic demand into LIRP considering returns in e-commerce and proposed a PPGASA algorithm as the solving approach.The above two researches mainly focus on the returns without quality defect but did not consider the MQDR. In this paper, we propose a model of closed-loop LIRP with MQDR. To the best of our knowledge, it is the first time to introduce the MQDR into LIRP in e-commerce. An effective hybrid algorithm named hybrid ant colony optimization (HACO) is provided to solve this model. Results of numerical instances indicate that HACO outperforms ant colony optimization (ACO) on optimal solution, iterations, and computing stability.The remainder of this paper is organized as follows. Section2 presents the mathematical model of LIRP with MQDR. Section 3 proposes the solution approach named HACO. Section 4 analyses the parameters of HACO and shows the results of different experiments. Section 5 gives the conclusion and future research directions.


## 2. Model Formulation

As we all know, customers’ return in e-commerce is higher than traditional commerce. Because of personal dissatisfaction, or a mistaken purchase of the wrong product, some of the returns are without quality defects. These returns can reenter into the market after a simple repackaging process without being recovered [29]. While the other returns result from quality defects, which need to be sent back to the plant and be recovered.In order to meet the needs of MQDR, the merchandise center (MC) is necessary to deliver normal merchandises to the demand points (DPs) of downstream and collect the returned merchandises. MC integrates the functions of distribution center and recycling center and provides quality inspection and repackaging services. Meanwhile the returned merchandises are collected to MCs. Returned merchandises without quality defects become resalable normal items after repackaging treatment at MCs. The plant will recover the returns with quality defects and bring them to the market again.The operation mode of the system is shown in Figure1. The closed-loop supply chain in this paper consists of one plant, multiple MCs, multiple DPs, and a single type of product with continuous inventory policy under the e-commerce environment.Figure 1
Closed-loop supply chain for a single product.The goal of this study is to decide the quantity and location of MCs and arrange the vehicle routes and determine the ordering times on each route. To minimize the total cost of logistics operations, this problem involves the following three decisions: (1) location decisions: obtain the optimal number of MCs and their locations; (2) inventory management: determine the ordering times on each route; (3) routing optimization: arrange the vehicles to delivery merchandises and collect returns.To benefit from the risk of MQDR, we take assumptions (1)–(8) from Li et al. [27] into consideration: since the single-product system is researched in this paper, assumption (1) is necessary; in the capacitated vehicle routing problem, assumption (2) should be satisfied [30]; assumption (3) eliminates the indeterminacy from the different type of vehicle; assumption (4) means that each DP is well served by the only vehicle route [31]; assumption (5) ensures that each route will return to the same MC after traversing; assumption (6) follows the early published papers considering uncapacitated MCs [32]; assumption (7) takes MCs as the distribution center and recycling center; assumption (8) is a simplification of the reality [33].The returned merchandises without quality defect are processed and repackaged at MCs, while others will be shipped back to the plant for reprocessing after a predetermined quantity at the MCs. Assume that the demand at each retailer is known and letR be the set of candidate MCs. Let S be the set of DPs and let M be the number of DPs. Let V be set of vehicles from the MCs to DPs. Let U
=
R
∪
S. Let the following notation denote the decisions of the firm:N
r
v: ordering times of M
C
r on routing v.
X
i
j
r
v: =1, if node j is served by M
C
r on routing v from node i, or 0 otherwise.
Y
i
r
v: =1, if node i is assigned to M
C
r on routing v, or 0 otherwise.
Z
r: =1, if M
C
r is selected as an MC location, or 0 otherwise.
U
j
k: auxiliary variable avoiding the subtour constraints in route k.According to the aforementioned assumptions, the inventory levels depend on both demand and the quantity of MQDR. So, during each replenishment cycle, the holding cost of MCs isλ
h
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
d
i
+
q
i
+
w
i
/
2
N
r
v
Y
i
r
v, where h is annual inventory holding per unit merchandises, d
i is mean (daily) demand for D
P
i, and q
i and w
i are quantity of merchandises without and with quality defect returned by D
P
i per day.In order to exactly describe the logistic distribution costs. Letc
r be the transportation costs per unit product from plan to M
C
r. Let l be the delivering cost per unit distance. Let s
i
j be the distance from node i to node j. And let λ be the working days per year. The total transportation costs from plant to DPs through MCs can be expressed, respectively, as λ
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
c
r
(
d
i
-
q
i
+
w
i
)
Y
i
r
v and λ
l
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
∑
j
∈
U
s
i
j
d
i
X
i
j
r
v.So the cost of forward distribution is(1)
∑
r
∈
R
f
r
Z
r
+
∑
v
∈
V
∑
r
∈
R
b
r
N
r
v
+
∑
v
∈
V
∑
r
∈
R
e
r
N
r
v
+
λ
h
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
d
i
+
q
i
+
w
i
2
N
r
v
Y
i
r
v
+
λ
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
c
r
d
i
-
q
i
+
w
i
Y
i
r
v
+
λ
l
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
∑
j
∈
U
s
i
j
d
i
X
i
j
r
v
,where e
r denote the fixed cost of dispatching vehicles per time at M
C
r, f
r denote the fixed (annual) administrative and construction cost of M
C
r, and b
r denote the ordering cost per unit product from plant to M
C
r.We letk to be the returning cost per unit of merchandise from DPs to MCs, so the total reverse transportation costs from DPs back to MCs are(2)
λ
k
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
q
i
+
w
i
Y
i
r
v
.The cost of deal with mixed quality defects is(3)
λ
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
p
r
q
i
Y
i
r
v
+
λ
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
a
r
q
i
+
w
i
Y
i
r
v
.We adopta
r as the inspecting cost per unit for the returned product and p
r as the repackaging cost of unit returned merchandise without quality problem at MCs.In summary, the model is formulated as follows:(4)
min
⁡
Z
=
∑
r
∈
R
f
r
Z
r
+
∑
v
∈
V
∑
r
∈
R
b
r
N
r
v
+
∑
v
∈
V
∑
r
∈
R
e
r
N
r
v
+
λ
h
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
d
i
+
q
i
+
w
i
2
N
r
v
Y
i
r
v
+
λ
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
c
r
d
i
-
q
i
+
w
i
Y
i
r
v
+
λ
l
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
∑
j
∈
U
s
i
j
d
i
X
i
j
r
v
+
λ
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
p
r
q
i
Y
i
r
v
+
λ
k
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
q
i
+
w
i
Y
i
r
v
+
λ
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
a
r
q
i
+
w
i
s
Y
i
r
v
.It is easy to find that the objective function (4) is convex in N
r
v. We can simplify the objective function by solving N
r
v. Consequently, the optimization solution of N
r
v can be obtained by taking the derivative of the function with respect to N
r
v; the result is as follows:(5)
N
r
v
∗
=
∑
i
∈
S
λ
h
d
i
+
q
i
+
w
i
Y
i
r
v
2
e
r
+
f
r
.The optimization problem (4) given a known N
r
v
∗ can now be written as(6)
min
⁡
Z
=
2
λ
h
∑
r
∈
R
∑
v
∈
V
∑
i
∈
S
f
r
+
e
r
d
i
+
q
i
+
w
i
Y
i
r
v
+
λ
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
c
r
d
i
-
q
i
+
w
i
Y
i
r
v
+
λ
l
∑
v
∈
V
∑
r
∈
R
∑
i
∈
U
∑
j
∈
U
s
i
j
d
i
X
i
j
r
v
+
λ
k
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
q
i
+
w
i
Y
i
r
v
+
λ
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
p
r
q
i
Y
i
r
v
+
λ
∑
v
∈
V
∑
r
∈
R
∑
i
∈
S
a
r
q
i
+
w
i
Y
i
r
v
+
∑
r
∈
R
f
r
Z
r
(7)
s.t.
∑
r
∈
R
Z
r
≥
1
,
r
∈
R
(8)
∑
v
∈
V
∑
i
∈
S
Y
i
r
v
≤
Z
r
,
r
∈
R
(9)
∑
i
∈
S
d
i
Y
i
r
v
≤
g
,
r
∈
R
,
v
∈
V
(10)
∑
v
∈
V
∑
r
∈
R
Y
i
r
v
=
1
,
i
∈
S
(11)
∑
v
∈
V
∑
r
∈
R
∑
i
∈
U
X
i
j
r
v
=
1
,
j
∈
S
(12)
∑
i
∈
U
X
i
k
r
v
-
∑
j
∈
U
X
k
j
r
v
=
0
,
k
∈
S
,
r
∈
R
,
v
∈
V
(13)
U
l
k
-
U
j
k
+
M
x
l
j
k
≤
M
-
1
,
l
,
j
∈
S
,
k
∈
R
(14)
-
Y
i
j
v
+
∑
k
∈
U
X
i
k
r
v
+
X
k
j
r
v
≤
1
,
i
∈
S
,
j
,
r
∈
R
,
v
∈
V
(15)
Z
r
=
0,1
,
r
∈
R
(16)
Y
i
r
v
=
0,1
,
i
∈
S
,
r
∈
R
,
v
∈
V
(17)
X
i
j
r
v
=
0,1
,
i
∈
U
,
j
∈
S
,
r
∈
R
,
v
∈
V
.The objective function (6) is to minimize the total cost; (7) ensure the selected MC is not empty; (8) ensure each DP is traversed by a unique vehicle which belongs to a certain MC; (9) ensure the amount of each delivery from MC on each route must be within vehicle capacity; (10) ensure that each route has only one vehicle; (11) ensure each DP must be followed by exactly one note; (12) ensure every DP node of the system will be serviced before it gives services to the others; (13) ensure the subtour is eliminated. Equation (14) ensures that each DP is assigned to an MC when there exists a route that starts from the MC passing through the DP. Equations (15)–(17) ensure the nonnegativity and integrality of decision variables.


## 3. Solution Approach

Like the VRP, the closed-loop LIRP is also an NP-hard problem, since it includes the VRP and is more complex than VRP. Generally speaking, there does not exist a complete, efficient, and accurate analytic algorithm to address NP-hard problems; ant colony optimization (ACO) has been proved very successful and widely applied to solve the static and dynamic problems as an EC algorithm [34]. However, ACO does not distinguish ant behavior results. The pheromone concentration will distribute in every direction for the iteration. Therefore, this leads to low searching efficiency. The algorithm may get caught in local optimization if we do not take preventive measures. On the other hand, ABC provides an effective institution to find the global optimal solution from the trapping of local optimal solution [35]. So, in this study, we present a hybrid ant colony optimization algorithm based on the combination of ACO and ABC to solve the above LIRP model.


### 3.1. Initialize Solution

Since the natural number is an efficient coding method for these problems, the sequence of solutions is composed of candidate MCs (1,2
,
…
,
R) and DPs, which are indicated by (R
+
1
,
…
,
R
+
S). The candidate solution of our proposed model will be described well by those natural number sequences. As an example, Figure 1 fully interpreted the perceptions of our method, which refers to the individual feasible solution: {
1
8
13
5
3
9
14
15
10
11
6
4
7
16
17
12
}.In the HACO, the moving strategy of the ant in nodei is depending on the pseudorandom proportional rule. The rule indicates that the ant has both exploiting and exploration ability, which means the ant is guided by the pheromone trails as well as the heuristic information. In this case, the ant has a higher degree of exploring unknown knowledge. The connected function of the pheromone values τ
i
j and the heuristic values η
i
j is shown as(18)
P
i
j
k
=
τ
i
j
α
η
i
j
β
∑
k
∈
a
l
l
o
w
k
τ
i
k
α
η
i
k
β
j
∈
a
l
l
o
w
k
0
o
t
h
e
r
w
i
s
e
,where τ
i
j is the density of pheromone remaining on the edge (
i
,
j
), η
i
j is the inverse value of distance between node i and node j, α and β are user-defined parameters for corresponding pheromone concentration and heuristic information, and a
l
l
o
w
(
k
) is the remaining nodes to be visited by ant k.


### 3.2. The ABC Phase

In order to improve the performance of global searching of our algorithm, the paper applied the scout bee searching phase into the ACO. Scout bees are free bees used for finding a new better solution from the neighbor known solution. As soon as a scout bee finds a new solution, she turns into an employed bee. If there is no improvement in the quality of solution, the bee will abandon that source and continue to search for another new solution.The searching function of scout bees is as(19)
x
^
i
j
=
x
m
a
x
-
rand
⁡
0,1
∗
x
m
a
x
-
x
i
j
.To meet the requirements for coding sequence type, we described two operations to complete scouts searching process, namely, random array reverse (RAR) and random swap (RS).Step 1.
Set the initial number of scout been and probability p
0.Step 2.
Generate two positions randomly nameda and b, for each a
<
b.Step 3.
Get a random probability0
<
p
<
1; if p
>
p
0 turn to Step 4; otherwise, turn to Step 5.Step 4.
Reverse the array between positionsa and b as a new solution.Step 5.
Swap the position ofa and b as a new solution.Step 6.
Calculate the cost of new solution.Step 7.
Keep the best solution to the next iteration and return to Step2.The pseudocodes of ABC are shown in Pseudocode1.Pseudocode 1:Pseudocode of an ABC framework.
Procedure: ABC
Input: the initial sequence, the number of scout been and p
0
Output: the better sequence
Begin
Taken and p
0
whilei
<
n
p
=
rand
⁡
[
0,1
]
takea
<
b
ifp
>
p
0
Reverse the array between positiona and b
else
Swap the position ofa and b
Output: the better sequence
End


### 3.3. Global Pheromone Trail Update

The global pheromone updating rule is triggered at the end of iteration to reward tours that are in line with the objective of impedance minimization. This strategy is applied to reinforce the pheromone density on the sets of edges belonging to the inspect tour and to increase the likelihood that this tour will also be selected by other ant agents. The rule of global pheromone updating is given by(20)
τ
i
j
=
1
-
ρ
τ
i
j
+
ρ
Δ
τ
i
j
b
s
i
,
j
∈
T
b
s
,where (21)
Δ
τ
i
j
b
s
=
L
L
b
s
i
,
j
∈
g
l
o
b
a
l
b
e
s
t
t
o
u
r
0
o
t
h
e
r
w
i
s
e
.L is a constant initial pheromone. L
b
s is the cost of the best of all the tours produced by all m agents from the beginning of the iteration. ρ
∈
(
0,1
] is the pheromone evaporation coefficient. To improve the pheromone trail quality, a part of the worst result is removed.


### 3.4. Local Pheromone Trail Update

In addition to the global pheromone trail updating rule, the selected ants will update the local pheromone trail in the process of passing ana
r
c
(
i
,
j
). It is opposite to the normal pheromone trail updating rule that increases the pheromone density while ants cross over arcs. The purpose of using the local pheromone trail update rule is to prevent stagnation behavior because the arc becomes less desirable for the following ants. The rule of local pheromone updating is as follows:(22)
τ
i
j
=
1
-
ξ
τ
i
j
+
ξ
τ
0
,where τ
0 is a constant at the beginning of pheromone trails and ξ is a user-defined coefficient that lowers the pheromone density of arcs traversed by the intelligent ants.


### 3.5. Algorithm Flow

Step 1.
Get the formulas for solvingN
r
v.Step 2.
Set the initial parameters for the model: set of candidate MCsR, set of DPs S, set of vehicles V, inspecting cost a
r, ordering cost b
r, transportation costs c
r, daily demand d
i, dispatching vehicles cost e
r, fixed (annual) administrative and construction cost f
r, vehicle capacity g, holding cost h, and returning cost k.Step 3.
Parameter setting for HACO is as follows: ants numberm, evolution terminate iteration M, pheromone concentration impact factor α, heuristic information pheromones impact factor β, evaporation rate of the pheromone ρ, constant initial pheromone L, and mutation probability array reverse p
0.Step 4.
Using unit matrixτ
i
j, calculate the probability(23)
P
i
j
k
=
τ
i
j
α
η
i
j
β
∑
k
∈
a
l
l
o
w
k
τ
i
k
α
η
i
k
β
j
∈
a
l
l
o
w
k
0
o
t
h
e
r
w
i
s
e
.Step 5.
Ant solutions generation module: each ant will generate a feasible solution after traversing the DPs.Step 6.
Best ant solution module: after calculating each ant’s solution, select the best solution which is known as the iteration best to compare with the global best. Keep the next best solution as the next global best.Step 7.
Scout bee module: random selection probability0
<
p
<
1; if p
>
p
0, turn to random array reverse operation. Otherwise, turn to random swap operation.Step 8.
Pheromone updating module: update the information pheromones as follows:(24)
τ
i
j
=
1
-
ρ
τ
i
j
+
ρ
Δ
τ
i
j
b
s
i
,
j
∈
T
b
s
,
τ
i
j
=
1
-
ξ
τ
i
j
+
ξ
τ
0
.Step 9.
Termination module: if the parent optimal solution and offspring optimal solution are equal during continuousM generations, stop the algorithm. Otherwise, return to Step 3 after M increments.Step 10.
Output.The pseudocodes of HACO are shown in Pseudocode2.Pseudocode 2:Pseudocode of the proposed HACO.
Procedure: HACO for LIRP
Input: coordinates of nodes, demands and returns of DPs, MC parameters, vehicle capacity, HACO parameters
Output: the best solution (include routes, MCs locations, order times and order size)
Begin
Takemax
⁡
t
whilemax
⁡
t
<
M
for 1 tom
Foraging Behavior of Ants
τ
i
j
=
(
1
-
ξ
)
τ
i
j
+
ξ
τ
0
Calculate individual total cost Tcost(m)
end
Tcost_best =min
⁡
(
T
c
o
s
t
(
m
)
)
for 1 tom
/
2 Neighbor range
Scout bee searching the neighbor range
a
=
r
o
u
n
d
(
r
a
n
d
)
b
=
r
o
u
n
d
(
r
a
n
d
)
makea
<
b
ifp
<
0.5
v
e
c
h
=
[
v
e
c
h
(
1
:
a
-
1
)
v
e
c
h
(
b
:
-
1
:
a
)
v
e
c
h
(
(
b
+
1
)
:
e
n
d
)
];
else
exchange(a
,
b)
end
if Tcost(m
/
2) < Tcost_best
Tcost_best = Tcost(m
/
2)
end
end
for 1 tom
if Tcost(m) > Tcost_ave
for 1 ton
τ
i
j
=
(
1
-
ρ
)
τ
i
j
+
ρ
Δ
τ
i
j
b
s
end
end
end
end
Output: the best solution
End


## 3.1. Initialize Solution

Since the natural number is an efficient coding method for these problems, the sequence of solutions is composed of candidate MCs (1,2
,
…
,
R) and DPs, which are indicated by (R
+
1
,
…
,
R
+
S). The candidate solution of our proposed model will be described well by those natural number sequences. As an example, Figure 1 fully interpreted the perceptions of our method, which refers to the individual feasible solution: {
1
8
13
5
3
9
14
15
10
11
6
4
7
16
17
12
}.In the HACO, the moving strategy of the ant in nodei is depending on the pseudorandom proportional rule. The rule indicates that the ant has both exploiting and exploration ability, which means the ant is guided by the pheromone trails as well as the heuristic information. In this case, the ant has a higher degree of exploring unknown knowledge. The connected function of the pheromone values τ
i
j and the heuristic values η
i
j is shown as(18)
P
i
j
k
=
τ
i
j
α
η
i
j
β
∑
k
∈
a
l
l
o
w
k
τ
i
k
α
η
i
k
β
j
∈
a
l
l
o
w
k
0
o
t
h
e
r
w
i
s
e
,where τ
i
j is the density of pheromone remaining on the edge (
i
,
j
), η
i
j is the inverse value of distance between node i and node j, α and β are user-defined parameters for corresponding pheromone concentration and heuristic information, and a
l
l
o
w
(
k
) is the remaining nodes to be visited by ant k.


## 3.2. The ABC Phase

In order to improve the performance of global searching of our algorithm, the paper applied the scout bee searching phase into the ACO. Scout bees are free bees used for finding a new better solution from the neighbor known solution. As soon as a scout bee finds a new solution, she turns into an employed bee. If there is no improvement in the quality of solution, the bee will abandon that source and continue to search for another new solution.The searching function of scout bees is as(19)
x
^
i
j
=
x
m
a
x
-
rand
⁡
0,1
∗
x
m
a
x
-
x
i
j
.To meet the requirements for coding sequence type, we described two operations to complete scouts searching process, namely, random array reverse (RAR) and random swap (RS).Step 1.
Set the initial number of scout been and probability p
0.Step 2.
Generate two positions randomly nameda and b, for each a
<
b.Step 3.
Get a random probability0
<
p
<
1; if p
>
p
0 turn to Step 4; otherwise, turn to Step 5.Step 4.
Reverse the array between positionsa and b as a new solution.Step 5.
Swap the position ofa and b as a new solution.Step 6.
Calculate the cost of new solution.Step 7.
Keep the best solution to the next iteration and return to Step2.The pseudocodes of ABC are shown in Pseudocode1.Pseudocode 1:Pseudocode of an ABC framework.
Procedure: ABC
Input: the initial sequence, the number of scout been and p
0
Output: the better sequence
Begin
Taken and p
0
whilei
<
n
p
=
rand
⁡
[
0,1
]
takea
<
b
ifp
>
p
0
Reverse the array between positiona and b
else
Swap the position ofa and b
Output: the better sequence
End


## 3.3. Global Pheromone Trail Update

The global pheromone updating rule is triggered at the end of iteration to reward tours that are in line with the objective of impedance minimization. This strategy is applied to reinforce the pheromone density on the sets of edges belonging to the inspect tour and to increase the likelihood that this tour will also be selected by other ant agents. The rule of global pheromone updating is given by(20)
τ
i
j
=
1
-
ρ
τ
i
j
+
ρ
Δ
τ
i
j
b
s
i
,
j
∈
T
b
s
,where (21)
Δ
τ
i
j
b
s
=
L
L
b
s
i
,
j
∈
g
l
o
b
a
l
b
e
s
t
t
o
u
r
0
o
t
h
e
r
w
i
s
e
.L is a constant initial pheromone. L
b
s is the cost of the best of all the tours produced by all m agents from the beginning of the iteration. ρ
∈
(
0,1
] is the pheromone evaporation coefficient. To improve the pheromone trail quality, a part of the worst result is removed.


## 3.4. Local Pheromone Trail Update

In addition to the global pheromone trail updating rule, the selected ants will update the local pheromone trail in the process of passing ana
r
c
(
i
,
j
). It is opposite to the normal pheromone trail updating rule that increases the pheromone density while ants cross over arcs. The purpose of using the local pheromone trail update rule is to prevent stagnation behavior because the arc becomes less desirable for the following ants. The rule of local pheromone updating is as follows:(22)
τ
i
j
=
1
-
ξ
τ
i
j
+
ξ
τ
0
,where τ
0 is a constant at the beginning of pheromone trails and ξ is a user-defined coefficient that lowers the pheromone density of arcs traversed by the intelligent ants.


## 3.5. Algorithm Flow

Step 1.
Get the formulas for solvingN
r
v.Step 2.
Set the initial parameters for the model: set of candidate MCsR, set of DPs S, set of vehicles V, inspecting cost a
r, ordering cost b
r, transportation costs c
r, daily demand d
i, dispatching vehicles cost e
r, fixed (annual) administrative and construction cost f
r, vehicle capacity g, holding cost h, and returning cost k.Step 3.
Parameter setting for HACO is as follows: ants numberm, evolution terminate iteration M, pheromone concentration impact factor α, heuristic information pheromones impact factor β, evaporation rate of the pheromone ρ, constant initial pheromone L, and mutation probability array reverse p
0.Step 4.
Using unit matrixτ
i
j, calculate the probability(23)
P
i
j
k
=
τ
i
j
α
η
i
j
β
∑
k
∈
a
l
l
o
w
k
τ
i
k
α
η
i
k
β
j
∈
a
l
l
o
w
k
0
o
t
h
e
r
w
i
s
e
.Step 5.
Ant solutions generation module: each ant will generate a feasible solution after traversing the DPs.Step 6.
Best ant solution module: after calculating each ant’s solution, select the best solution which is known as the iteration best to compare with the global best. Keep the next best solution as the next global best.Step 7.
Scout bee module: random selection probability0
<
p
<
1; if p
>
p
0, turn to random array reverse operation. Otherwise, turn to random swap operation.Step 8.
Pheromone updating module: update the information pheromones as follows:(24)
τ
i
j
=
1
-
ρ
τ
i
j
+
ρ
Δ
τ
i
j
b
s
i
,
j
∈
T
b
s
,
τ
i
j
=
1
-
ξ
τ
i
j
+
ξ
τ
0
.Step 9.
Termination module: if the parent optimal solution and offspring optimal solution are equal during continuousM generations, stop the algorithm. Otherwise, return to Step 3 after M increments.Step 10.
Output.The pseudocodes of HACO are shown in Pseudocode2.Pseudocode 2:Pseudocode of the proposed HACO.
Procedure: HACO for LIRP
Input: coordinates of nodes, demands and returns of DPs, MC parameters, vehicle capacity, HACO parameters
Output: the best solution (include routes, MCs locations, order times and order size)
Begin
Takemax
⁡
t
whilemax
⁡
t
<
M
for 1 tom
Foraging Behavior of Ants
τ
i
j
=
(
1
-
ξ
)
τ
i
j
+
ξ
τ
0
Calculate individual total cost Tcost(m)
end
Tcost_best =min
⁡
(
T
c
o
s
t
(
m
)
)
for 1 tom
/
2 Neighbor range
Scout bee searching the neighbor range
a
=
r
o
u
n
d
(
r
a
n
d
)
b
=
r
o
u
n
d
(
r
a
n
d
)
makea
<
b
ifp
<
0.5
v
e
c
h
=
[
v
e
c
h
(
1
:
a
-
1
)
v
e
c
h
(
b
:
-
1
:
a
)
v
e
c
h
(
(
b
+
1
)
:
e
n
d
)
];
else
exchange(a
,
b)
end
if Tcost(m
/
2) < Tcost_best
Tcost_best = Tcost(m
/
2)
end
end
for 1 tom
if Tcost(m) > Tcost_ave
for 1 ton
τ
i
j
=
(
1
-
ρ
)
τ
i
j
+
ρ
Δ
τ
i
j
b
s
end
end
end
end
Output: the best solution
End


## 4. Computational Experiments and Results Analysis

In this section, numerical simulations are given to illustrate the performance of HACO compared with the traditional ACO. Both algorithms in this paper are compiled byMatlab R2014a and run on a computer with 8 GB main memory and 3.6 GHZ CPU. All instances come from the LRP database in University of Aveiro [36].


### 4.1. Parameters Discussion

Parameter values selection is crucial to the efficiency of algorithms. An example named Gaskell 67-22 × 5 from the database, which contains the nodes coordinate and the DPs demand, is used to determine the optimal parameter. Gaskell 67 is the instance’s name and 22 × 5 means 5 candidate MCs for 22 DPs. The inventory holding costh
=
2, the vehicle capacity g
=
500, transportation costs c
r
=
2, returning cost k
=
2, working days λ
=
300, and the delivering cost per unit distance l
=
0.7. The other parameters of the instance are as follows: a
r
~
U
(
16,20
); b
r
~
U
(
16,20
); c
r
~
U
(
6,10
); e
r
~
U
(
21,25
); q
i
~
U
(
12,25
); and w
i
~
U
(
2,5
).The parameters of algorithm are initialized as follows: ant’s numberm
=
50, evolution terminate iteration M
=
100, pheromone concentration impact factor α
=
1, heuristic information pheromones impact factor β
=
5, evaporation rate of the pheromone ρ
=
0.1, constant L
=
150, mutation probability array reverse p
0
=
0.5.We run the program 50 times on the same computer. The performance of ACO and HACO varies with the different values of the parameters, which are shown in Tables1–6. In these tables, the symbol C.V. means the coefficient of variation.Table 1
Results with differentm (cost: million CNY; iterations: times).
m
Cost
Iterations
Mean
Std. dev.
C.V.
Mean
Std. dev.
C.V.
30
ACO
32.4
0.96
0.0295
168.6
30.64
0.1817
HACO
33.46
1.12
0.0336
125.72
29.02
0.2308
40
ACO
32.13
0.96
0.03
173.48
32.95
0.1899
HACO
33.37
0.93
0.0277
117.82
15.4
0.1307
50
ACO
32.22
0.74
0.0229
157.19
40.58
0.2581
HACO
32.96
0.96
0.0292
119.26
17.47
0.1465
60
ACO
32.21
0.67
0.0208
158.12
39.39
0.2491
HACO
32.72
0.86
0.0263
121.86
17.37
0.1425
70
ACO
31.96
0.71
0.0223
162.88
32.48
0.1994
HACO
32.77
1.13
0.0344
121.06
23.7
0.1958
α
=
1, β
=
5, ρ
=
0.1, L
=
150, and p
0
=
0.5.Table 2
Results with differentα (cost: million CNY; iterations: times).
α
Cost
Iterations
Mean
Std. dev.
C.V.
Mean
Std. dev.
C.V.
0.25
ACO
31.92
0.54
0.0168
189.3
73.28
0.3871
HACO
31.63
0.79
0.025
189.38
69.64
0.3677
0.5
ACO
31.85
0.78
0.0246
181.86
57.34
0.3153
HACO
31.52
0.91
0.029
182.24
63.94
0.3509
1
ACO
32.22
0.74
0.0229
157.19
40.58
0.2581
HACO
33.06
1.04
0.0315
123.28
14.02
0.1137
1.25
ACO
32.13
0.79
0.0247
154.24
30.12
0.1953
HACO
33.11
1.09
0.0331
115.78
11.61
0.1002
1.5
ACO
32.14
0.9
0.0279
154.96
25.68
0.1658
HACO
32.97
1.11
0.0335
116.2
10.23
0.0881
m
=
50, β
=
5, ρ
=
0.1, L
=
150, and p
0
=
0.5.Table 3
Results with differentβ (cost: million CNY; iterations: times).
β
Cost
Iterations
Mean
Std. dev.
C.V.
Mean
Std. dev.
C.V.
3
ACO
31.91
0.91
0.0286
192.78
42.08
0.2183
HACO
33.06
1.19
0.0359
148.8
54.84
0.3686
4
ACO
32.19
0.8
0.0249
170.2
41.89
0.2461
HACO
33.05
1.14
0.0344
126
42.17
0.3347
5
ACO
32.22
0.74
0.0229
157.19
40.58
0.2581
HACO
32.96
0.96
0.0292
119.26
17.47
0.1465
6
ACO
32.1
0.77
0.0241
159.38
31.99
0.2007
HACO
33.12
1.19
0.0358
119.86
22.4
0.1869
7
ACO
32.42
0.84
0.0258
154.69
23.92
0.1546
HACO
33.25
1.07
0.0322
111.68
8.62
0.0772
α
=
1, m
=
50, ρ
=
0.1, L
=
150, and p
0
=
0.5.Table 4
Results with differentρ (cost: million CNY; iterations: number).
ρ
Cost
Iterations
Mean
Std. dev.
C.V.
Mean
Std. dev.
C.V.
0.1
ACO
32.22
0.74
0.0229
157.19
40.58
0.2581
HACO
32.12
0.93
0.0289
138.4
21.02
0.1518
0.2
ACO
32.38
0.87
0.0269
133.27
20.57
0.1543
HACO
32.7
0.97
0.0296
124.12
13.01
0.1048
0.3
ACO
32.65
0.94
0.0289
132.8
31.09
0.2341
HACO
32.75
1.45
0.0443
122.1
13.09
0.1072
0.4
ACO
33.03
0.94
0.0284
131.86
28.76
0.2181
HACO
33.31
1.1
0.033
122.98
20.2
0.1643
0.5
ACO
33.41
0.87
0.0261
157.27
67.72
0.4306
HACO
33.55
1.08
0.0323
116.24
15.44
0.1328
α
=
1, m
=
50, β
=
5, L
=
150, and p
0
=
0.5.Table 5
Results with differentL (cost: million CNY; iterations: times).
L
Cost
Iterations
Mean
Std. dev.
C.V.
Mean
Std. dev.
C.V.
50
ACO
31.99
0.72
0.0224
159.72
33.39
0.21
HACO
33.07
1.22
0.0369
123.7
25.96
0.2098
100
ACO
31.9
0.67
0.0209
159.8
30.73
0.19
HACO
32.88
1.09
0.033
121.63
18.91
0.1555
150
ACO
32.22
0.74
0.0229
157.19
40.58
0.2581
HACO
33.02
1.05
0.0319
121.14
23.2
0.1915
200
ACO
32.11
0.8
0.0251
167.66
38.61
0.23
HACO
32.9
1.12
0.0339
120.8
22.05
0.1825
250
ACO
32.08
0.87
0.0271
164.08
34.16
0.21
HACO
33.1
1.09
0.033
121.87
28.95
0.2376
α
=
1, m
=
50, β
=
5, ρ
=
0.1, and p
0
=
0.5.Table 6
Results with differentp
0 (cost: million CNY; iterations: times).
p
0
Cost
Iterations
Mean
Std. dev.
C.V.
Mean
Std. dev.
C.V.
0.3
33.07
1.22
0.0369
123.70
25.96
0.2098
0.4
32.88
1.09
0.0330
121.63
18.91
0.1555
0.5
33.02
1.05
0.0319
121.14
23.20
0.1915
0.6
32.90
1.12
0.0339
120.80
22.05
0.1825
0.7
33.10
1.09
0.0330
121.87
28.95
0.2376
α
=
1, m
=
50, β
=
5, ρ
=
0.1, and L
=
150.Tables1–6 represent the parameters’ effect on the objective function values. The data was normalized through two dimensions, that is, cost and iterations, and three indicators, that is, mean, standard deviation, and coefficient of variation. Actually, in order to find the minimal cost, we usually take the parameter values, where the cost is lower and more stable. From the discussion, we found that HACO reaches the best performance when m
=
60, α
=
0.5, β
=
5, ρ
=
0.1, L
=
100, and p
0
=
0.4, while ACO reaches the best performance when m
=
70, α
=
0.5, β
=
3, ρ
=
0.1, and L
=
100.


### 4.2. Computational Experiment

To get a reliable conclusion, we run another 50 times on the same computer with the best parameter values in Gaskell 67-22 × 5. One of the best solutions of objective function in the 50 experiments of HACO is 30.2 million CNY. Table7 shows the solution. MCs were established at MC1, MC2, and MC5 with five vehicles distribution routes. Figure 2 shows topological structure of the closed-loop supply chain.Table 7
The solution of Gaskell 67-22 × 5.
MC
Routing number
Routing
Order times
MC1
V1
1
→
12
→
13
→
9
→
10
→
26
→
1
17
V2
1
→
16
→
15
→
18
→
14
→
1
18
MC2
V3
2
→
23
→
24
→
27
→
25
→
17
→
2
37
MC5
V4
5
→
22
→
20
→
21
→
8
→
19
→
5
24
V5
5
→
7
→
6
→
11
→
5
18Figure 2
Topological structure of the network.Figure3 shows their trend of optimal objective function values along with iterations. The fluctuation curves of optimal objective function value are varied by different algorithm, which are shown in Figures 4(a) and 4(b).Figure 3
Trends of objective function value.Figure 4
The fluctuation curve of optimal objective function value.
(a)
ACO
(b)
HACOAs shown in Figure3, the cost and iterations of HACO are lower than ACO; and in Figure 4, the range and mean value of the minimum cost of HACO are also lower than ACO, which both imply that HACO is more efficient than ACO in solving the LIRP.


### 4.3. Extended Experiments

In this section, a series of instances are given to show that HACO is more efficient and stable than classical software and ACO. In order to ensure the demands of DPs are not more than the vehicle capacity, we need to enumerate some instances. In this paper, the daily demands are set as 1/10 of corresponding demands of the database.As we know, Lingo is a representative classical optimization software tool. Thus we usedLingo 11.0 to solve the problem by using a small-sized instance named Peal 183-12 × 2 and two medium-sized instance named Gaskell 67-22 × 5 and Gaskell 67-36 × 5; the results are shown in Table 8.Table 8
Comparisons between HACO and Lingo.
Instance name
Perl 183-12
×
2
Gaskell 67-22
×
5
Gaskell 67-36
×
5
CPU time
Cost
CPU time
Cost
CPU time
Cost
Lingo
523 s
715392
>1 hour
∖
>1 hour
∖
HACO
24.5 s
709152
44 s
32691202
99 s
31550446Each instance was run 50 times by HACO and ACO with their optimized parameters values, respectively; the results are shown in Tables9 and 10.Table 9
Optimal objective function values of two algorithms (CNY).
Instance name
Algorithm
Mean
Std. dev.
C.V.
t value
Significance test
p value (sig. 1-tailed)
Perl 183-12
×
2
ACO
711347.1
6979.11
0.0098
1.6838
0.048
HACO
709237.9
5454.48
0.0077
Gaskell 67-22
×
5
ACO
31929818.3
627210.66
0.0196
1.7118
0.045
HACO
31718084.1
609606.04
0.0192
Gaskell 67-36
×
5
ACO
3288921.3
43585.25
0.0133
1.6602
0.050
HACO
3275052
39868.81
0.0122
Perl 183-55
×
15
ACO
285145.3
3779.15
0.0133
1.8391
0.034
HACO
283831.6
3351.00
0.0118
Christofides 69-75
×
10
ACO
419695.8
4948.15
0.0118
2.5020
0.007
HACO
417298.6
4627.32
0.0111
Perl 183-85
×
7
ACO
486011.0
4503.39
0.0093
1.8080
0.037
HACO
484424.1
4270.44
0.0088
Christofides 69-100
×
10
ACO
458235.6
6387.00
0.0139
1.6796
0.048
HACO
456441.6
5317.70
0.0117Table 10
Iterations of two algorithms (times).
Instance name
Algorithm
Mean
Std. dev.
C.V.
t value
Significance test
p value (sig. 1-tailed)
Perl 183-12
×
2
ACO
193.32
78.70
0.4071
0.6323
0.264
HACO
183.90
70.04
0.3808
Gaskell 67-22
×
5
ACO
191.34
86.36
0.4513
0.7831
0.218
HACO
180.70
43.32
0.2397
Gaskell 67-36
×
5
ACO
254.44
76.84
0.3020
0.2012
0.421
HACO
251.84
49.42
0.1962
Perl 183-55
×
15
ACO
241.74
77.29
0.3197
0.1000
0.460
HACO
240.46
47.33
0.1968
Christofides 69-75
×
10
ACO
245.24
79.10
0.3226
1.2340
0.110
HACO
229.92
62.30
0.2710
Perl 183-85
×
7
ACO
247.18
94.89
0.3839
0.4355
0.332
HACO
240.80
65.75
0.2730
Christofides 69-100
×
10
ACO
248.06
97.18
0.3918
0.2627
0.397
HACO
244.08
65.48
0.2681


### 4.4. Result Analysis

According to Table8, we found that (1) for the small-sized instance, HACO can obtain better result than Lingo within less time and (2) for medium-sized instances, Lingo cannot get the global optimization within 1 hour, while HACO can solve the problem in a short time.Observe, from Tables9 and 10, that HACO is more efficient than ACO for the following reasons. (1) The cost of HACO is significantly lower than ACO (p
<
0.05); (2) the difference of the number of iteration between HACO and ACO is not significant (p
>
0.05); (3) HACO is more stable than ACO as the coefficient of variation (C.V.) is lower. To sum up, our algorithm reduces the cost with the same number of iterations compared with ACO.By improving pheromone updates and bee colony searching, we improve the solution quality of the algorithm and make it useful as a guide for the ant searching process. Observed from the results of numerical simulations, HACO can get better result with a fewer number of iterations. Hence, comparing with ACO, HACO is adopted as a better approach in solving this LIRP with MQDR.


## 4.1. Parameters Discussion

Parameter values selection is crucial to the efficiency of algorithms. An example named Gaskell 67-22 × 5 from the database, which contains the nodes coordinate and the DPs demand, is used to determine the optimal parameter. Gaskell 67 is the instance’s name and 22 × 5 means 5 candidate MCs for 22 DPs. The inventory holding costh
=
2, the vehicle capacity g
=
500, transportation costs c
r
=
2, returning cost k
=
2, working days λ
=
300, and the delivering cost per unit distance l
=
0.7. The other parameters of the instance are as follows: a
r
~
U
(
16,20
); b
r
~
U
(
16,20
); c
r
~
U
(
6,10
); e
r
~
U
(
21,25
); q
i
~
U
(
12,25
); and w
i
~
U
(
2,5
).The parameters of algorithm are initialized as follows: ant’s numberm
=
50, evolution terminate iteration M
=
100, pheromone concentration impact factor α
=
1, heuristic information pheromones impact factor β
=
5, evaporation rate of the pheromone ρ
=
0.1, constant L
=
150, mutation probability array reverse p
0
=
0.5.We run the program 50 times on the same computer. The performance of ACO and HACO varies with the different values of the parameters, which are shown in Tables1–6. In these tables, the symbol C.V. means the coefficient of variation.Table 1
Results with differentm (cost: million CNY; iterations: times).
m
Cost
Iterations
Mean
Std. dev.
C.V.
Mean
Std. dev.
C.V.
30
ACO
32.4
0.96
0.0295
168.6
30.64
0.1817
HACO
33.46
1.12
0.0336
125.72
29.02
0.2308
40
ACO
32.13
0.96
0.03
173.48
32.95
0.1899
HACO
33.37
0.93
0.0277
117.82
15.4
0.1307
50
ACO
32.22
0.74
0.0229
157.19
40.58
0.2581
HACO
32.96
0.96
0.0292
119.26
17.47
0.1465
60
ACO
32.21
0.67
0.0208
158.12
39.39
0.2491
HACO
32.72
0.86
0.0263
121.86
17.37
0.1425
70
ACO
31.96
0.71
0.0223
162.88
32.48
0.1994
HACO
32.77
1.13
0.0344
121.06
23.7
0.1958
α
=
1, β
=
5, ρ
=
0.1, L
=
150, and p
0
=
0.5.Table 2
Results with differentα (cost: million CNY; iterations: times).
α
Cost
Iterations
Mean
Std. dev.
C.V.
Mean
Std. dev.
C.V.
0.25
ACO
31.92
0.54
0.0168
189.3
73.28
0.3871
HACO
31.63
0.79
0.025
189.38
69.64
0.3677
0.5
ACO
31.85
0.78
0.0246
181.86
57.34
0.3153
HACO
31.52
0.91
0.029
182.24
63.94
0.3509
1
ACO
32.22
0.74
0.0229
157.19
40.58
0.2581
HACO
33.06
1.04
0.0315
123.28
14.02
0.1137
1.25
ACO
32.13
0.79
0.0247
154.24
30.12
0.1953
HACO
33.11
1.09
0.0331
115.78
11.61
0.1002
1.5
ACO
32.14
0.9
0.0279
154.96
25.68
0.1658
HACO
32.97
1.11
0.0335
116.2
10.23
0.0881
m
=
50, β
=
5, ρ
=
0.1, L
=
150, and p
0
=
0.5.Table 3
Results with differentβ (cost: million CNY; iterations: times).
β
Cost
Iterations
Mean
Std. dev.
C.V.
Mean
Std. dev.
C.V.
3
ACO
31.91
0.91
0.0286
192.78
42.08
0.2183
HACO
33.06
1.19
0.0359
148.8
54.84
0.3686
4
ACO
32.19
0.8
0.0249
170.2
41.89
0.2461
HACO
33.05
1.14
0.0344
126
42.17
0.3347
5
ACO
32.22
0.74
0.0229
157.19
40.58
0.2581
HACO
32.96
0.96
0.0292
119.26
17.47
0.1465
6
ACO
32.1
0.77
0.0241
159.38
31.99
0.2007
HACO
33.12
1.19
0.0358
119.86
22.4
0.1869
7
ACO
32.42
0.84
0.0258
154.69
23.92
0.1546
HACO
33.25
1.07
0.0322
111.68
8.62
0.0772
α
=
1, m
=
50, ρ
=
0.1, L
=
150, and p
0
=
0.5.Table 4
Results with differentρ (cost: million CNY; iterations: number).
ρ
Cost
Iterations
Mean
Std. dev.
C.V.
Mean
Std. dev.
C.V.
0.1
ACO
32.22
0.74
0.0229
157.19
40.58
0.2581
HACO
32.12
0.93
0.0289
138.4
21.02
0.1518
0.2
ACO
32.38
0.87
0.0269
133.27
20.57
0.1543
HACO
32.7
0.97
0.0296
124.12
13.01
0.1048
0.3
ACO
32.65
0.94
0.0289
132.8
31.09
0.2341
HACO
32.75
1.45
0.0443
122.1
13.09
0.1072
0.4
ACO
33.03
0.94
0.0284
131.86
28.76
0.2181
HACO
33.31
1.1
0.033
122.98
20.2
0.1643
0.5
ACO
33.41
0.87
0.0261
157.27
67.72
0.4306
HACO
33.55
1.08
0.0323
116.24
15.44
0.1328
α
=
1, m
=
50, β
=
5, L
=
150, and p
0
=
0.5.Table 5
Results with differentL (cost: million CNY; iterations: times).
L
Cost
Iterations
Mean
Std. dev.
C.V.
Mean
Std. dev.
C.V.
50
ACO
31.99
0.72
0.0224
159.72
33.39
0.21
HACO
33.07
1.22
0.0369
123.7
25.96
0.2098
100
ACO
31.9
0.67
0.0209
159.8
30.73
0.19
HACO
32.88
1.09
0.033
121.63
18.91
0.1555
150
ACO
32.22
0.74
0.0229
157.19
40.58
0.2581
HACO
33.02
1.05
0.0319
121.14
23.2
0.1915
200
ACO
32.11
0.8
0.0251
167.66
38.61
0.23
HACO
32.9
1.12
0.0339
120.8
22.05
0.1825
250
ACO
32.08
0.87
0.0271
164.08
34.16
0.21
HACO
33.1
1.09
0.033
121.87
28.95
0.2376
α
=
1, m
=
50, β
=
5, ρ
=
0.1, and p
0
=
0.5.Table 6
Results with differentp
0 (cost: million CNY; iterations: times).
p
0
Cost
Iterations
Mean
Std. dev.
C.V.
Mean
Std. dev.
C.V.
0.3
33.07
1.22
0.0369
123.70
25.96
0.2098
0.4
32.88
1.09
0.0330
121.63
18.91
0.1555
0.5
33.02
1.05
0.0319
121.14
23.20
0.1915
0.6
32.90
1.12
0.0339
120.80
22.05
0.1825
0.7
33.10
1.09
0.0330
121.87
28.95
0.2376
α
=
1, m
=
50, β
=
5, ρ
=
0.1, and L
=
150.Tables1–6 represent the parameters’ effect on the objective function values. The data was normalized through two dimensions, that is, cost and iterations, and three indicators, that is, mean, standard deviation, and coefficient of variation. Actually, in order to find the minimal cost, we usually take the parameter values, where the cost is lower and more stable. From the discussion, we found that HACO reaches the best performance when m
=
60, α
=
0.5, β
=
5, ρ
=
0.1, L
=
100, and p
0
=
0.4, while ACO reaches the best performance when m
=
70, α
=
0.5, β
=
3, ρ
=
0.1, and L
=
100.


## 4.2. Computational Experiment

To get a reliable conclusion, we run another 50 times on the same computer with the best parameter values in Gaskell 67-22 × 5. One of the best solutions of objective function in the 50 experiments of HACO is 30.2 million CNY. Table7 shows the solution. MCs were established at MC1, MC2, and MC5 with five vehicles distribution routes. Figure 2 shows topological structure of the closed-loop supply chain.Table 7
The solution of Gaskell 67-22 × 5.
MC
Routing number
Routing
Order times
MC1
V1
1
→
12
→
13
→
9
→
10
→
26
→
1
17
V2
1
→
16
→
15
→
18
→
14
→
1
18
MC2
V3
2
→
23
→
24
→
27
→
25
→
17
→
2
37
MC5
V4
5
→
22
→
20
→
21
→
8
→
19
→
5
24
V5
5
→
7
→
6
→
11
→
5
18Figure 2
Topological structure of the network.Figure3 shows their trend of optimal objective function values along with iterations. The fluctuation curves of optimal objective function value are varied by different algorithm, which are shown in Figures 4(a) and 4(b).Figure 3
Trends of objective function value.Figure 4
The fluctuation curve of optimal objective function value.
(a)
ACO
(b)
HACOAs shown in Figure3, the cost and iterations of HACO are lower than ACO; and in Figure 4, the range and mean value of the minimum cost of HACO are also lower than ACO, which both imply that HACO is more efficient than ACO in solving the LIRP.


## 4.3. Extended Experiments

In this section, a series of instances are given to show that HACO is more efficient and stable than classical software and ACO. In order to ensure the demands of DPs are not more than the vehicle capacity, we need to enumerate some instances. In this paper, the daily demands are set as 1/10 of corresponding demands of the database.As we know, Lingo is a representative classical optimization software tool. Thus we usedLingo 11.0 to solve the problem by using a small-sized instance named Peal 183-12 × 2 and two medium-sized instance named Gaskell 67-22 × 5 and Gaskell 67-36 × 5; the results are shown in Table 8.Table 8
Comparisons between HACO and Lingo.
Instance name
Perl 183-12
×
2
Gaskell 67-22
×
5
Gaskell 67-36
×
5
CPU time
Cost
CPU time
Cost
CPU time
Cost
Lingo
523 s
715392
>1 hour
∖
>1 hour
∖
HACO
24.5 s
709152
44 s
32691202
99 s
31550446Each instance was run 50 times by HACO and ACO with their optimized parameters values, respectively; the results are shown in Tables9 and 10.Table 9
Optimal objective function values of two algorithms (CNY).
Instance name
Algorithm
Mean
Std. dev.
C.V.
t value
Significance test
p value (sig. 1-tailed)
Perl 183-12
×
2
ACO
711347.1
6979.11
0.0098
1.6838
0.048
HACO
709237.9
5454.48
0.0077
Gaskell 67-22
×
5
ACO
31929818.3
627210.66
0.0196
1.7118
0.045
HACO
31718084.1
609606.04
0.0192
Gaskell 67-36
×
5
ACO
3288921.3
43585.25
0.0133
1.6602
0.050
HACO
3275052
39868.81
0.0122
Perl 183-55
×
15
ACO
285145.3
3779.15
0.0133
1.8391
0.034
HACO
283831.6
3351.00
0.0118
Christofides 69-75
×
10
ACO
419695.8
4948.15
0.0118
2.5020
0.007
HACO
417298.6
4627.32
0.0111
Perl 183-85
×
7
ACO
486011.0
4503.39
0.0093
1.8080
0.037
HACO
484424.1
4270.44
0.0088
Christofides 69-100
×
10
ACO
458235.6
6387.00
0.0139
1.6796
0.048
HACO
456441.6
5317.70
0.0117Table 10
Iterations of two algorithms (times).
Instance name
Algorithm
Mean
Std. dev.
C.V.
t value
Significance test
p value (sig. 1-tailed)
Perl 183-12
×
2
ACO
193.32
78.70
0.4071
0.6323
0.264
HACO
183.90
70.04
0.3808
Gaskell 67-22
×
5
ACO
191.34
86.36
0.4513
0.7831
0.218
HACO
180.70
43.32
0.2397
Gaskell 67-36
×
5
ACO
254.44
76.84
0.3020
0.2012
0.421
HACO
251.84
49.42
0.1962
Perl 183-55
×
15
ACO
241.74
77.29
0.3197
0.1000
0.460
HACO
240.46
47.33
0.1968
Christofides 69-75
×
10
ACO
245.24
79.10
0.3226
1.2340
0.110
HACO
229.92
62.30
0.2710
Perl 183-85
×
7
ACO
247.18
94.89
0.3839
0.4355
0.332
HACO
240.80
65.75
0.2730
Christofides 69-100
×
10
ACO
248.06
97.18
0.3918
0.2627
0.397
HACO
244.08
65.48
0.2681


## 4.4. Result Analysis

According to Table8, we found that (1) for the small-sized instance, HACO can obtain better result than Lingo within less time and (2) for medium-sized instances, Lingo cannot get the global optimization within 1 hour, while HACO can solve the problem in a short time.Observe, from Tables9 and 10, that HACO is more efficient than ACO for the following reasons. (1) The cost of HACO is significantly lower than ACO (p
<
0.05); (2) the difference of the number of iteration between HACO and ACO is not significant (p
>
0.05); (3) HACO is more stable than ACO as the coefficient of variation (C.V.) is lower. To sum up, our algorithm reduces the cost with the same number of iterations compared with ACO.By improving pheromone updates and bee colony searching, we improve the solution quality of the algorithm and make it useful as a guide for the ant searching process. Observed from the results of numerical simulations, HACO can get better result with a fewer number of iterations. Hence, comparing with ACO, HACO is adopted as a better approach in solving this LIRP with MQDR.


## 5. Conclusion and Future Research

With the development of e-commerce, customers’ return keeps a high rate with MQDR, which can be reentered into markets after being repackaged or recovered. In this research, we built a closed-loop LIRP model considering both quality defect returns and nondefect returns; we call it MQDR in this paper. We perform an extensive computational study and observe the following interesting results.(1)
Considering MQDR are computationally beneficial for the formulation presented, the MQDR and closed-loop pattern with returns are features of the proposed problem in e-commerce, which is never considered in previous work.
(2)
Since the evolutionary computation algorithm has been proved successfully in tackling NP-hard problem, a hybrid algorithm is proposed by combining ACO algorithm and ABC algorithm to solve the LIRP. HACO integrated the scout bee searching phase into the ACO to improve the global searching ability.
(3)
The performance of HACO is evaluated by using the instances in the LRP database, and HACO outperforms ACO on convergence, optimal solution, and computing stability. This numerical study shows the efficiency and effectiveness of the solution method.However, developing other elements for the LIRP will lead to further research directions. And analyzing the model under the dynamic demand of customs and a time-varying demand can be a valuable subject. The design of experiments and verification by discrete dynamics simulation should be established. Fruit fly optimization algorithm (FOA) as one of the best EC algorithms has attracted the attention of various researchers [37]. It is important to apply these models and algorithms to the operation and management of enterprises to improve the decision-making efficiency of e-commerce logistics system.



---

*Source: 6467812-2016-06-30.xml*