# Design of a Multiobjective Reverse Logistics Network Considering the Cost and Service Level

Authors: Shuang Li; Nengmin Wang; Zhengwen He; Ada Che; Yungao Ma
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

Reverse logistics, which is induced by various forms of used products and materials, has received growing attention throughout this decade. In a highly competitive environment, the service level is an important criterion for reverse logistics network design. However, most previous studies about product returns only focused on the total cost of the reverse logistics and neglected the service level. To help a manufacturer of electronic products provide quality postsale repair service for their consumer, this paper proposes a multiobjective reverse logistics network optimisation model that considers the objectives of the cost, the total tardiness of the cycle time, and the coverage of customer zones. The Nondominated Sorting Genetic Algorithm II (NSGA-II) is employed for solving this multiobjective optimisation model. To evaluate the performance of NSGA-II, a genetic algorithm based on weighted sum approach and Multiobjective Simulated Annealing (MOSA) are also applied. The performance of these three heuristic algorithms is compared using numerical examples. The computational results show that NSGA-II outperforms MOSA and the genetic algorithm based on weighted sum approach. Furthermore, the key parameters of the model are tested, and some conclusions are drawn.


---

## Body

## 1. Introduction

Reverse logistics is the process of planning, implementing, and controlling the flow of raw materials, in-process inventory, finished goods, and related information from the point of consumption to the point of recovery or the point of proper disposal [1]. Reverse logistics involves activities such as the return, reconditioning, refurbishment and recycling of products, and packaging. In recent years, for a variety of economic, environmental, or legislative reasons, reverse logistics has received increasing attention from industry and academia.Reverse logistics network design is a major strategic issue that determines the number, location and capacity of the collection points and centralised return centres. Effective reverse logistics is believed to benefit a company substantially. First, it can reduce the cost and improve the utilisation rate of materials [2]. Second, it can increase the profit of an enterprise and build a good enterprise reputation [3, 4]. Moreover, it is a good way to improve customer satisfaction and loyalty and then to maintain a sustainable competitive advantage [5, 6].The establishment of appropriate performance measures is a key factor in implementing a successful reverse logistics system. The common performance measures include cost minimisation, customer satisfaction maximisation, cycle time minimisation, flexibility, and the overall efficiency of the reverse logistics system. These standards are usually divided into two parts, the cost and the service [7]. Many researchers have studied reverse logistics network design for different industries, but most of them focused only on the overall cost, and few researchers have considered the service level to be another criterion in their model. The minimisation of the cost is commonly a major concern to be considered when building a reverse logistics network system, but the service level is also a key factor when determining the survival and development of a company under the current economic environment, which is driven by customer values. To service providers, both the service level and the total service cost are major concerns [8]. A well-managed reverse logistics network cannot only provide important cost savings in procurement, recovery, disposal, inventory holding, and transportation but can also help in customer retention [9]. Amini et al. [5] argued that the management of service activities such as a repair service, product upgrades, and product disposal can form an important part of a corporate strategy. For the manufacturers that produce electronic products such as computers, mobile phones, and cameras, it is of paramount importance to improve customer satisfaction as well as to provide products of high quality.This paper proposes a multiobjective reverse logistics network optimisation model for a manufacturer’s postsale repair service. Two additional objectives that are associated with the service level are considered besides the minimisation of the overall cost. Thus, three objectives are considered in this paper. These objectives are the minimisation of the total reverse logistics cost, the minimisation of the total tardiness of the cycle time, and the maximisation of the rate of customer zones covered within the acceptable service coverage of collection points (coverage of customer zones). Our purpose is to find a set of nondominated solutions that determine the number and location of the collection points and the repair centres among the potential facility locations as well as the associated transportation flows between the customer zones and the service facilities. To deal with multiobjective and obtain a set of Pareto-optimal solutions, NSGA-II is implemented for the proposed model. To evaluate the performance of NSGA-II, a genetic algorithm based on weighted sum approach (called GA_WS hereafter) and MOSA are also applied. Computational experiments are conducted to compare the performance of the three algorithms, and the key parameters of the model are tested.This paper is organised as follows. In the next section, a literature review on a reverse logistics network design is presented. In Section3, a multiobjective reverse logistics network design problem that involves a postsale service level is formulated. Section 4 develops a multiobjective optimisation model for the proposed problem. Section 5 describes the three multiobjective evolutionary algorithms (MOEAs), namely, NSGA-II, GA_WS, and MOSA. Section 6 gives the computational results for comparing the performance of the three algorithms and testing the key parameters of the model. Finally, the conclusions are presented and future research directions are highlighted in Section 7.


## 2. Literature Review

Although reverse logistics is a new field that has obtained attention only over the past decade, many scholars have performed research on reverse logistics. Fleischmann et al. [10] divided reverse logistics into three main areas, namely, distribution planning, inventory control, and production planning. They reviewed quantitative mathematical models proposed in the literature for each of these areas. Some researchers have proposed decision conceptual frameworks on reverse logistics such as De Brito and Dekker [3] and Lambert et al. [4]. They both proposed a decision conceptual framework for reverse logistics in terms of strategic, tactic, and operational decisions. Furthermore, they noted that reverse logistics network design belongs to strategic decisions.The design of a product recovery network is one of the important and challenging problems in the field of reverse logistics [11]. Many researchers have conducted quantitative analyses of product recovery networks and have proposed mathematical models. However, the majority of studies have concentrated only on the overall cost or profit, including Alumur et al. [12], Barros et al. [13], Cruz-Rivera and Ertel [14], Das and Chowdhury [15], Hu et al. [16], Jayaraman et al. [17], Lee and Dong [18], Lieckens and Vandaele [19], Min et al. [20], and Salema et al. [21]. They each proposed a mathematical model and considered the total cost or profit to be a single objective function in their studies. However, in the real world, there are no design tasks that are single objective problems [22]. The design problems usually involve trade-offs among multiple and conflicting objectives, such as cost, resource utilisation, and service level. Some researchers have studied the multiobjective optimisation of reverse logistic networks.Ioannis [23] proposed a multiobjective model for locating disposal or treatment facilities and for transporting hazardous waste along the links of a transportation network. The objectives considered the minimisation of the total operating cost, the minimisation of the total perceived risk, the equitable distribution of risk among population centres and the equitable distribution of disutility caused by the operation of treatment facilities, and a goal programming approach was proposed to solve the problem. Pati et al. [24] formulated a mixed integer goal programming model to assist in the proper management of the paper recycling logistics system. The objectives considered were a reduction in the reverse logistics cost, product quality improvement through increased segregation at the source, and environmental benefits through increased wastepaper recovery. Ahluwalia and Nema [25] proposed a multiobjective reverse logistics model for integrated computer waste management. This model was based on an integer linear programming approach with the objective of minimising the environmental risk as well as the cost. The studies reviewed above about product recovery or recycling network design all considered the environmental aspects in addition to the overall cost and found a Pareto-optimal solution or a restrictive set of Pareto-optimal solutions based on their solution approaches for the problem.The substantial challenge of globalisation and the fierce competition of markets prompt more and more manufacturing firms to implement effective reverse logistics networks for product returns and to provide better consumer postsale service. Amini et al. [5] indicated that an important means for companies to differentiate themselves as well as to increase profitability in highly competitive environments is through the use of service management. They further argued that one of the most important service management activities is repair services, which represent important opportunities to create profit streams and to strengthen customer loyalty. Du and Evans [8] proposed a biobjective reverse logistics network optimisation model for a manufacturer’s postsale repair service, considering the service level in addition to the overall cost. The objectives were the minimisation of the overall costs and the minimisation of the total tardiness of the cycle time. The solution approach consisted of a combination of three algorithms: a scatter search, the dual simplex method, and the constraint method. Zarandi et al. [26] addressed a closed-loop supply chain distribution network design problem in which reverse flows were imported into forward model. They proposed three multiobjective models considering the covering objectives as the measure of service level besides the total cost. A fuzzy goal programming approach was developed for the problem.To help a computer manufacturer implement quality postsale repair service, this paper proposes a multiobjective reverse logistics network optimisation model considering both cost and service level. In this model, both the total tardiness of the cycle time and the coverage of customer zones are considered as measures of the service level. As computers are daily necessities, if the products can be repaired and returned to the customers within a satisfactory period, customer satisfaction can be definitely improved [5]. Because customers who return products do not prefer long distances, the collection points must be located within a certain maximum distance from them [9]. Thus, a multiobjective optimisation model is proposed for the reverse logistic network optimisation problem.


## 3. Problem Definition

The problem considered in this paper is from a company that is a manufacturer of a variety of computer equipment. With the purchase of a computer, the customers are promised a two-year warranty. Because of improvements in people’s living standards and fast developments in our society, the demand for computers increases every year. Increased purchases and higher quality standards have dramatically increased the volume of returned products. It is time for the company to implement an effective reverse logistics network that satisfies capacity limitations and demand requirements for the postsale repair service of the increased volume of product returns.The company considers taking advantage of the established facilities of the third-party logistics provider as collection points and establishing centralised return centres to guarantee product safety and technology privacy. In this way, customers return products that need repair to the collection point. The products are then transported to the repair centres. After repair, the products are quickly delivered back to the collection points. Then, customers are called to take the products back from the collection points.The company considers the service level to be another important criterion aside from the total cost of the reverse logistics network and focuses on the convenience of customers when returning their products and the efficiency of the repair service. To make it convenient for customers to return products, the collection points will be set within a certain maximised distance to the customer zones, specifically to maximise the rate of the customer zones covered within the acceptable service distance of collection points. To improve customer satisfaction, the returned products should be repaired and sent back to customers within their expected cycle time as often as possible, to minimise the total tardiness of the cycle time. Thus, the coverage of customer zones and the total tardiness of the cycle time are used to measure the service level of the repair service.To address the problems that face the company, we formulate the reverse logistics network problem as a multiobjective integer nonlinear programming model. Three objectives are considered:(1) minimisation of the total reverse logistics cost, (2) minimisation of the total tardiness of the cycle time, and (3) maximisation of the coverage of customer zones.The main issues to be addressed by this model are the following.(1)
Which locations are to be chosen for the collection points and the repair centres?(2)
How many collection points and repair centres are needed?(3)
How to best arrange the transportation flows that start from the customer zones, go through the collection points and repair centres and then go back to the consumers?To summarise, this paper proposes a multiobjective optimisation model for a three-echelon reverse logistics network design problem, which determines the optimal location and the number of both the collection points and repair centres and the transportation flows between the customer zones and the facility sites.


## 4. Problem Formulation

Prior to developing the multiobjective optimisation model, we make the following assumptions.(1)
The possibility of direct shipment from customers to a repair centre is ruled out.(2)
Given the small volume of individual returns from customers, a collection point has a sufficient capacity to hold the returned products.(3)
The transportation costs between the customers and their nearest collection points are negligible because of the short distances between the customers and their nearest collection point.(4)
The location/allocation plan covers a planning horizon within which no substantial changes are incurred in the customer demands and in the transportation infrastructure.


### 4.1. Indices

i:
Index for the customer zones(i∈I)j:
Index for the collection points(j∈J)k:
Index for the repair centres(k∈K).


### 4.2. Model Parameters

ri:
The daily volume of products returned by the customer zoneiw:
The annual working daysa:
The annual cost of rent for a collection pointh:
The handling cost of a product unitek:
The annual average construction cost of a repair centredij:
The distance from customer zonei to collection point jdjk:
The distance from the collection pointj to the repair centre ktm:
The average number of hours for repairing a piece of returned producttjk:
The round trip transportation time between the collection pointj and the repair centre kte:
The cycle time expected by customersmk:
The maximum capacity of a repair centrel:
The maximum allowable distance from a given customer zone to a collection point (service coverage)(4.1)pij={1ifdij≤l,0otherwise.cp:
The minimum number of  open collection pointsrc:
The minimum number of established repair centresM:
An arbitrarily chosen large number.We have thef(Xjk,djk)=Eαβ function for the freight rate, where α is a discount rate according to the volume of shipment between collection point j and repair centre k; β is a penalty rate applied for the distance between collection point j and repair centre k; E is a unit freight rate.Xjk=∑iriYijWjk: volume of products returned from collection point j to repair centre k(4.2)α={1Xjk≤p1,α1p1<Xjk≤p2,α2Xjk>p2,(4.3)β={1djk≤q1,β1q1<djk≤q2,β2djk>q2.p1,p2: Volume of returned products for a discountq1,q2: Distance between collection point j and repair centre k for penalties.


### 4.3. Decision Variables

(4.4)
Y
i
j
=
{
1
if
consumer
zone
i
is
allocated
to
collection
point
j
,
0
otherwise
,
W
j
k
=
{
1
if
collection
point
j
is
allocated
to
repair
center
k
,
0
otherwise
,
G
k
=
{
1
if
a
repair
center
is
established
at
site
k
,
0
otherwise
,
Z
j
=
{
1
if
a
collection
point
is
established
at
site
j
,
0
otherwise
.


### 4.4. Mathematical Formulation

(4.5)
min
f
1
=
a
∑
j
Z
j
+
h
w
∑
i
r
i
+
∑
k
e
k
G
k
+
∑
k
∑
j
(
∑
i
r
i
w
Y
i
j
W
j
k
)
×
f
(
X
j
k
,
d
j
k
)
,
(4.6)
min
f
2
=
∑
k
∑
j
∑
i
r
i
Y
i
j
W
j
k
max
{
t
m
+
t
j
k
-
t
e
,
0
}
,
(4.7)
max
f
3
=
∑
i
∑
j
r
i
p
i
j
Y
i
j
∑
i
r
i
,
subject to
(4.8)∑jYij=1,∀i∈I,(4.9)∑kWjk=Zj,∀j∈J,(4.10)Zj≤∑iYij≤M·Zj,∀j∈J,(4.11)Gk≤∑jWjk≤M·Gk,∀k∈K,(4.12)∑j∑iriYijWjk≤mkGk,∀k∈K,(4.13)cp≤∑jZj,(4.14)rc≤∑kGk,(4.15)Yij,Wjk,Zj,Gk∈(0,1)∀i∈I,∀j∈J,∀k∈K.The objective function (4.5) minimises the total cost of the reverse logistics network, including the annual rent cost of the collection points, the annual average construction cost of the repair centres, the materials handling cost, and the transportation cost. The objective function (4.6) minimises the total tardiness of the cycle time. The objective function (4.7) maximises the coverage of consumer zones. Constraint (4.8) assures that a customer zone is assigned to a single collection point. Constraint (4.9) ensures that an established collection point is allocated to a single repair centre. Constraint (4.10) prevents any customer zone from being assigned to the unopened collection points and assures that there must be some customer zones that are assigned to an opened collection point. Constraint (4.11) prevents any return flows from the collection points to the unopened repair centres and ensures that there must be some collection points that are assigned to an opened repair centre. Constraint (4.12) ensures that the total volume of the products returned from the collection points does not exceed the maximum capacity of a repair centre. Constraint (4.13) and Constraint (4.14) maintain a minimum number of collection points and repair centres for product return. Constraint (4.15) assures binary integer values for the decision variables Yij,Wjk,Zj,andGk.


### 4.5. Complexity of the Model

The maximum covering problem is usually defined as follows: given a connected network with demand at nodes, locate one or more facility sites at nodes in such a way as to maximise the coverage of demand nodes [27]. In the proposed model of this paper, some collection points are located to maximise the coverage of customer zones (returned products). The covered customer zones are within the acceptable service distance of collection points. Without any loss of generality, let ∑jZj=cp. If the proposed model is simplified to the problem with a single objective of maximizing the coverage of customer zones (objective function (4.7)) subject to constraint (4.8), constraint (4.10), constraint (4.13) and constraint (4.15), the simplified model can be regarded as a special case of the maximum covering model. Because the maximum covering problem has been proved to be NP-hard [28], our model which considers two additional objectives besides the maximisation of coverage of customer zones must be NP-hard as well.


## 4.1. Indices

i:
Index for the customer zones(i∈I)j:
Index for the collection points(j∈J)k:
Index for the repair centres(k∈K).


## 4.2. Model Parameters

ri:
The daily volume of products returned by the customer zoneiw:
The annual working daysa:
The annual cost of rent for a collection pointh:
The handling cost of a product unitek:
The annual average construction cost of a repair centredij:
The distance from customer zonei to collection point jdjk:
The distance from the collection pointj to the repair centre ktm:
The average number of hours for repairing a piece of returned producttjk:
The round trip transportation time between the collection pointj and the repair centre kte:
The cycle time expected by customersmk:
The maximum capacity of a repair centrel:
The maximum allowable distance from a given customer zone to a collection point (service coverage)(4.1)pij={1ifdij≤l,0otherwise.cp:
The minimum number of  open collection pointsrc:
The minimum number of established repair centresM:
An arbitrarily chosen large number.We have thef(Xjk,djk)=Eαβ function for the freight rate, where α is a discount rate according to the volume of shipment between collection point j and repair centre k; β is a penalty rate applied for the distance between collection point j and repair centre k; E is a unit freight rate.Xjk=∑iriYijWjk: volume of products returned from collection point j to repair centre k(4.2)α={1Xjk≤p1,α1p1<Xjk≤p2,α2Xjk>p2,(4.3)β={1djk≤q1,β1q1<djk≤q2,β2djk>q2.p1,p2: Volume of returned products for a discountq1,q2: Distance between collection point j and repair centre k for penalties.


## 4.3. Decision Variables

(4.4)
Y
i
j
=
{
1
if
consumer
zone
i
is
allocated
to
collection
point
j
,
0
otherwise
,
W
j
k
=
{
1
if
collection
point
j
is
allocated
to
repair
center
k
,
0
otherwise
,
G
k
=
{
1
if
a
repair
center
is
established
at
site
k
,
0
otherwise
,
Z
j
=
{
1
if
a
collection
point
is
established
at
site
j
,
0
otherwise
.


## 4.4. Mathematical Formulation

(4.5)
min
f
1
=
a
∑
j
Z
j
+
h
w
∑
i
r
i
+
∑
k
e
k
G
k
+
∑
k
∑
j
(
∑
i
r
i
w
Y
i
j
W
j
k
)
×
f
(
X
j
k
,
d
j
k
)
,
(4.6)
min
f
2
=
∑
k
∑
j
∑
i
r
i
Y
i
j
W
j
k
max
{
t
m
+
t
j
k
-
t
e
,
0
}
,
(4.7)
max
f
3
=
∑
i
∑
j
r
i
p
i
j
Y
i
j
∑
i
r
i
,
subject to
(4.8)∑jYij=1,∀i∈I,(4.9)∑kWjk=Zj,∀j∈J,(4.10)Zj≤∑iYij≤M·Zj,∀j∈J,(4.11)Gk≤∑jWjk≤M·Gk,∀k∈K,(4.12)∑j∑iriYijWjk≤mkGk,∀k∈K,(4.13)cp≤∑jZj,(4.14)rc≤∑kGk,(4.15)Yij,Wjk,Zj,Gk∈(0,1)∀i∈I,∀j∈J,∀k∈K.The objective function (4.5) minimises the total cost of the reverse logistics network, including the annual rent cost of the collection points, the annual average construction cost of the repair centres, the materials handling cost, and the transportation cost. The objective function (4.6) minimises the total tardiness of the cycle time. The objective function (4.7) maximises the coverage of consumer zones. Constraint (4.8) assures that a customer zone is assigned to a single collection point. Constraint (4.9) ensures that an established collection point is allocated to a single repair centre. Constraint (4.10) prevents any customer zone from being assigned to the unopened collection points and assures that there must be some customer zones that are assigned to an opened collection point. Constraint (4.11) prevents any return flows from the collection points to the unopened repair centres and ensures that there must be some collection points that are assigned to an opened repair centre. Constraint (4.12) ensures that the total volume of the products returned from the collection points does not exceed the maximum capacity of a repair centre. Constraint (4.13) and Constraint (4.14) maintain a minimum number of collection points and repair centres for product return. Constraint (4.15) assures binary integer values for the decision variables Yij,Wjk,Zj,andGk.


## 4.5. Complexity of the Model

The maximum covering problem is usually defined as follows: given a connected network with demand at nodes, locate one or more facility sites at nodes in such a way as to maximise the coverage of demand nodes [27]. In the proposed model of this paper, some collection points are located to maximise the coverage of customer zones (returned products). The covered customer zones are within the acceptable service distance of collection points. Without any loss of generality, let ∑jZj=cp. If the proposed model is simplified to the problem with a single objective of maximizing the coverage of customer zones (objective function (4.7)) subject to constraint (4.8), constraint (4.10), constraint (4.13) and constraint (4.15), the simplified model can be regarded as a special case of the maximum covering model. Because the maximum covering problem has been proved to be NP-hard [28], our model which considers two additional objectives besides the maximisation of coverage of customer zones must be NP-hard as well.


## 5. The Solution Approach

Traditionally, there are several algorithms for solving the multiobjective optimisation problems, including theε-constraint method, objective programming, and weighting approach. However, these approaches can be used to find only a Pareto-optimal solution or a restrictive set of Pareto-optimal solutions. To obtain a diverse set of Pareto-optimal solutions and enable the decision maker for evaluating a greater number of solutions, MOEAs are widely applied. They can simultaneously address a set of possible solutions in one single run. Among all of the MOEAs, NSGA-II proposed by Deb et al. [29] is one of the most popular algorithms with more accuracy and a higher convergence speed. It has been widely used by many researchers, such as Fallah-Mehdipour et al. [30], Gutjahr et al. [31], Kannan et al. [32], Lin and Yeh [33], and Saadatseresht et al. [34]. With these facts in mind, NSGA-II is selected to solve the NP-hard problem proposed in this paper.To validate the results obtained using NSGA-II, GA_WS, which has been widely applied since the last century [22, 35, 36], is adopted as well. Since both NSGA-II and GA_WS are the extended states of Genetic algorithm (GA), MOSA proposed by Suppapitnarm et al. [37] is also employed to investigate the effectiveness of NSGA-II further.


### 5.1. Common Features

Representation is an important issue for a successful implementation of GA and SA. This paper adopts the binary coding method. Each chromosome is based on a single-dimensional array, which consists of binary values and represents decision variables in terms of the opening or closing of collection points and repair centres. For example, the representation of a chromosome is illustrated in Figure1. Each collection point or repair centre has one gene that represents an opening or closing decision. The chromosome has 10 collection points and 5 repair centres. As shown in Figure 1, collection point 1, collection point 6, collection point 7, repair centre 1, and repair centre 4 are open.Figure 1
A representation scheme: cp stands for a collection point and rc stands for a repair centre.The decision variables for the opening or closing of the collection points and repair centres can be obtained when generating the initial population of the GA or the new solutions in MOSA. Other decision variables, which involve the transport flow from customer zones to the repair centres through the collection points, can be obtained by two assignment algorithms. The first assignment algorithm is used for obtaining the total daily demand of the opened collection points. In other words, each customer zone should be assigned to the nearest collection points because we assumed that there is sufficient capacity at each collection point as a result of the small volume of returns. The second assignment algorithm is used for assigning opened collection points to an appropriate repair centre according to a capacity limitation. To solve this problem, we applied the Vogel method for a transportation problem.


### 5.2. NSGA-II

NSGA-II derives a new generation from the current generation by a mechanism that includes three different modules: (i) fast nondominated sort, (ii) crowding distance assignment, and (iii) the crowded comparison operator.A solutioni is said to constrained-dominate a solution j if any of the following conditions are true.(1)
Solutioni is feasible, and solution j is not feasible.(2)
Solutionsi and j are both infeasible, but solution i has a smaller overall constraint violation.(3)
Solutionsi and j are feasible, and solution i dominates solution j [29].First, the fast nondominated sorting approach is used to sort the initial population into different ranks according to their nondomination level. Second, the density estimation is employed for individuals of each nondomination level to obtain an estimate of the density of the solutions surrounding them. And then the genetic operators are used to create a child population. Thereafter, the created child population is combined with the parent population to form a combined population of size2N. Then the entire population is sorted according to nondominated. The new population (size N) is formed by adding solutions from the first front until the size exceeds N. To choose exactly the population members, the crowded-comparison operator is used to sort the solutions of the last accepted front, and the best solutions that are needed to fill all of the population slots are chosen. The new population of size N is now used for selection, crossover, and mutation to create a new population.The genetic operators are described below.(
1
) Select
As a selection mechanism, a binary tournament selection strategy was adopted by forming two teams of chromosomes. Each team consists of two chromosomes that are randomly selected from the current population. The solution with the lower (better) rank is selected if the two solutions are from different front. The solution with the higher crowding distance is selected if both of the solutions belong to the same front.(
2
) Crossover
The crossover operator generates new children by exchanging parts of the strings of a pair of selected parents. Here, we employed the two-point crossover in which one point is used for locating the collection points and the other point is used for locating repair centres. The two crossover points are randomly selected. Figure2 is an illustration of the crossover operator.Figure 2
An illustration of the crossover operator.(
3
) Mutation
A mutation is usually performed by modifying a gene within a chromosome. Here, we applied the multipoint mutation by randomly selecting five bit values of opening/closing decision variables of collection points and three bit values of opening/closing decision variables of repair centres. Then, the bit value is changed from 0 to 1 or from 1 to 0.


### 5.3. GA_WS

To determine the weight values, we adopted the adaptive weighting approach proposed by Gen and Cheng [38]. This method determines the weights based on the ideal point that was generated in each evolutionary process. Figure 3 illustrates the weight strategy in the objective space. The adaptive weight of each objective for an individual in the current generation is determined using (5.1)
(5.1)wi=1fimax-fimin,
where fimax and fimin are the maximum and minimum values of the ith objective in the current generation, respectively.Figure 3
Illustration of Gen and Cheng’s weight strategy.In terms of the problem consisting ofm objective functions, the weighted-sum objective function of a given solution x is decided by the following
(5.2)f(x)=∑i=1mwi(fi-fimin)=∑i=1mfi-fiminfimax-fimin.The smaller the weighted-sum objective value is, the better the individual. The fitness function is formed by adding a penalty to the total objective function. The penaltyfunction=M*px, where M is the penalty value, which is larger than any possible objective function value. px is calculated by the following
(5.3)px=∑k(∑j∑iriYijWjk-mkGk)if∑j∑iriYijWjk>mkGk,∀k∈K,otherwise0.The initial population is generated randomly. Here, we apply the same genetic operators as those used in NSGA-II. Different from NSGA-II, the selection operator is based on the fitness value of the individual. To find multiple optimal solutions in one single simulation run, the Pareto-optimal set is created by the nondominated solutions in the initial population and is updated by new individuals that are obtained with genetic operators at every generation.


### 5.4. MOSA

The description of the MOSA procedure is as follows.(1)
An initial solutionx is randomly generated. If it is infeasible, regenerate it until it is feasible. Calculate the objective function values of the initial solution x and put it into the Pareto-optimal set.(2)
The mutation operator in GA_WS is chosen as a moving strategy to obtain a new solutiony. If it is infeasible, regenerate it until it is feasible. Then the new solution y is compared with every solution in the Pareto-optimal set and the Pareto-optimal set is updated.(3)
If the new solutiony is accepted into the Pareto-optimal set, substitute the current solution x with y. Otherwise, the new solution  y  is accepted based on the following probability:  P=min(1,∏i=1Mexp{(fi(x)-fi(y))/Ti}), where Ti is the current temperature. If it is still not accepted, the current solution remains.(4)
Randomly select a solution from the Pareto-optimal set every certain number (iter_num) of generations as the current solution and go on searching.(5)
In every iteration, when the number of evaluated solutions reaches the set value (iter_max) or the number of times of sequential rejection of a new solution exceeds the set value (un_max), the temperature is reduced at a certain rate (b). The search process terminates when the current temperature drops to the final temperature or the total number of the evaluated solutions reaches the set value.


## 5.1. Common Features

Representation is an important issue for a successful implementation of GA and SA. This paper adopts the binary coding method. Each chromosome is based on a single-dimensional array, which consists of binary values and represents decision variables in terms of the opening or closing of collection points and repair centres. For example, the representation of a chromosome is illustrated in Figure1. Each collection point or repair centre has one gene that represents an opening or closing decision. The chromosome has 10 collection points and 5 repair centres. As shown in Figure 1, collection point 1, collection point 6, collection point 7, repair centre 1, and repair centre 4 are open.Figure 1
A representation scheme: cp stands for a collection point and rc stands for a repair centre.The decision variables for the opening or closing of the collection points and repair centres can be obtained when generating the initial population of the GA or the new solutions in MOSA. Other decision variables, which involve the transport flow from customer zones to the repair centres through the collection points, can be obtained by two assignment algorithms. The first assignment algorithm is used for obtaining the total daily demand of the opened collection points. In other words, each customer zone should be assigned to the nearest collection points because we assumed that there is sufficient capacity at each collection point as a result of the small volume of returns. The second assignment algorithm is used for assigning opened collection points to an appropriate repair centre according to a capacity limitation. To solve this problem, we applied the Vogel method for a transportation problem.


## 5.2. NSGA-II

NSGA-II derives a new generation from the current generation by a mechanism that includes three different modules: (i) fast nondominated sort, (ii) crowding distance assignment, and (iii) the crowded comparison operator.A solutioni is said to constrained-dominate a solution j if any of the following conditions are true.(1)
Solutioni is feasible, and solution j is not feasible.(2)
Solutionsi and j are both infeasible, but solution i has a smaller overall constraint violation.(3)
Solutionsi and j are feasible, and solution i dominates solution j [29].First, the fast nondominated sorting approach is used to sort the initial population into different ranks according to their nondomination level. Second, the density estimation is employed for individuals of each nondomination level to obtain an estimate of the density of the solutions surrounding them. And then the genetic operators are used to create a child population. Thereafter, the created child population is combined with the parent population to form a combined population of size2N. Then the entire population is sorted according to nondominated. The new population (size N) is formed by adding solutions from the first front until the size exceeds N. To choose exactly the population members, the crowded-comparison operator is used to sort the solutions of the last accepted front, and the best solutions that are needed to fill all of the population slots are chosen. The new population of size N is now used for selection, crossover, and mutation to create a new population.The genetic operators are described below.(
1
) Select
As a selection mechanism, a binary tournament selection strategy was adopted by forming two teams of chromosomes. Each team consists of two chromosomes that are randomly selected from the current population. The solution with the lower (better) rank is selected if the two solutions are from different front. The solution with the higher crowding distance is selected if both of the solutions belong to the same front.(
2
) Crossover
The crossover operator generates new children by exchanging parts of the strings of a pair of selected parents. Here, we employed the two-point crossover in which one point is used for locating the collection points and the other point is used for locating repair centres. The two crossover points are randomly selected. Figure2 is an illustration of the crossover operator.Figure 2
An illustration of the crossover operator.(
3
) Mutation
A mutation is usually performed by modifying a gene within a chromosome. Here, we applied the multipoint mutation by randomly selecting five bit values of opening/closing decision variables of collection points and three bit values of opening/closing decision variables of repair centres. Then, the bit value is changed from 0 to 1 or from 1 to 0.


## 5.3. GA_WS

To determine the weight values, we adopted the adaptive weighting approach proposed by Gen and Cheng [38]. This method determines the weights based on the ideal point that was generated in each evolutionary process. Figure 3 illustrates the weight strategy in the objective space. The adaptive weight of each objective for an individual in the current generation is determined using (5.1)
(5.1)wi=1fimax-fimin,
where fimax and fimin are the maximum and minimum values of the ith objective in the current generation, respectively.Figure 3
Illustration of Gen and Cheng’s weight strategy.In terms of the problem consisting ofm objective functions, the weighted-sum objective function of a given solution x is decided by the following
(5.2)f(x)=∑i=1mwi(fi-fimin)=∑i=1mfi-fiminfimax-fimin.The smaller the weighted-sum objective value is, the better the individual. The fitness function is formed by adding a penalty to the total objective function. The penaltyfunction=M*px, where M is the penalty value, which is larger than any possible objective function value. px is calculated by the following
(5.3)px=∑k(∑j∑iriYijWjk-mkGk)if∑j∑iriYijWjk>mkGk,∀k∈K,otherwise0.The initial population is generated randomly. Here, we apply the same genetic operators as those used in NSGA-II. Different from NSGA-II, the selection operator is based on the fitness value of the individual. To find multiple optimal solutions in one single simulation run, the Pareto-optimal set is created by the nondominated solutions in the initial population and is updated by new individuals that are obtained with genetic operators at every generation.


## 5.4. MOSA

The description of the MOSA procedure is as follows.(1)
An initial solutionx is randomly generated. If it is infeasible, regenerate it until it is feasible. Calculate the objective function values of the initial solution x and put it into the Pareto-optimal set.(2)
The mutation operator in GA_WS is chosen as a moving strategy to obtain a new solutiony. If it is infeasible, regenerate it until it is feasible. Then the new solution y is compared with every solution in the Pareto-optimal set and the Pareto-optimal set is updated.(3)
If the new solutiony is accepted into the Pareto-optimal set, substitute the current solution x with y. Otherwise, the new solution  y  is accepted based on the following probability:  P=min(1,∏i=1Mexp{(fi(x)-fi(y))/Ti}), where Ti is the current temperature. If it is still not accepted, the current solution remains.(4)
Randomly select a solution from the Pareto-optimal set every certain number (iter_num) of generations as the current solution and go on searching.(5)
In every iteration, when the number of evaluated solutions reaches the set value (iter_max) or the number of times of sequential rejection of a new solution exceeds the set value (un_max), the temperature is reduced at a certain rate (b). The search process terminates when the current temperature drops to the final temperature or the total number of the evaluated solutions reaches the set value.


## 6. Computational Results

In this section, the performance of the three algorithms was compared on five sets of fifteen randomly generated instances, and the key parameters of the model were tested.


### 6.1. Date Generation

To compare the performance of the three algorithms, five sets of fifteen randomly generated instances were created in Table1. The number of daily returned products by each customer zone was randomly generated to be between 10 and 50. The round trip transportation time between the collection point j  and the repair centre k was computed as follows: tjk=djk×0.6. Other parameters of the model are shown in Table 2.Table 1
Parameters of five problem sets.
Problem sets
1
2
3
4
5
Number of customer zones
50
60
80
100
120
Number of potential collection points
20
30
35
40
50
Number of potential repair centres
8
10
10
12
12
Location scale of all points
0~100
0~100
0~150
0~150
0~200Table 2
Parameters of the model.
Parameter
Index
Value
Annual renting cost of the collection points
a
200
Handling cost per unit product
h
0.1
Capacity of a repair centre
m
k
2000
Minimum number of open collection points
cp
5
Minimum number of established repair centres
rc
1
Average working hours to repair a returned product
t
m
10
α
1
0.8
Discount rate with respect to the shipping volume
α
2
0.6
p
1
200 units
p
2
400 units
β
1
1.1
Penalty rate with respect to the shipping distance
β
2
1.2
q
1
25
q
2
60
Working days per year
w
250
Cost of establishing a repair centre
e
k
3000
Service coverage
l
25
Unit standard transportation cost
E
1
Expected cycle time
t
e
30


### 6.2. Comparison of the Three Algorithms

To compare the Pareto-optimal set that is obtained by three algorithms on the same basis, the number of solutions searched is taken to be the stopping criterion. According to the different sizes of the five Problem sets, the number of solutions searched is set to be 120,000, 122,000, 128,000, 130,000, and 132,000. Because both GA_WS and NSGA-II are based on a GA, the same parameters are set for both. Based on extensive experiments, the parameters are as follows: the population size = 40, the crossover rate = 0.96, and the mutationrate=0.07. Because the population size of GA_WS or NSGA-II is 40, the maximum numbers of generations for GA_WS and NSGA-II on the five Problem sets are 3,000, 3,050, 3,200, 3,250, and 3,300, respectively. For MOSA, the initial temperature Ts is calculated from the following
(6.1)Ts=ΔSUMmaxlnProb,
where ΔSUMmax is the maximum value of the sum of differences between the maximum and the minimum of the three objective values, which are chosen from the objective values of 30 randomly generated neighbours of the initial solution. The value of Prob is set to 0.95 to ensure that a feasible inferior solution is accepted with a probability of 0.95. Other parameters are as follows: iter_max=1000,un_max=50, iter_num=40, and b=0.95. The search process terminates when the total number of evaluated solutions reaches the set value. These three algorithms are all coded in the C++ programming language in the VC++6.0 environment and are executed on a Dell Intel Core 2 Duo computer with a speed of 2.10 GHz and with 2.00 GB of memory.The Pareto-optimal set is the set of Pareto-optimal solutions that consists of all decision vectors for which the corresponding objective vectors cannot be improved in a given dimension without worsening another solution [39]. To evaluate the performance of the three algorithms, we adopted the following standards. Standard (1), the average number of Pareto-optimal solutions; Standard (2), the average number of nondominated solutions; Standard (3), the average ratio of the Pareto-optimal solutions. Due to the stochastic nature of the suggested algorithms, these standards were obtained by the three algorithms over 10 runs. Standard (2) and Standard (3) were calculated in the following manner. Let P1, P2, and P3 be the sets of Pareto-optimal solutions that are obtained from one run of GA_WS, NSGA-II, and MOSA, respectively, and let P be the union of the sets of Pareto-optimal solutions (i.e., P=P1∪P2∪P3), with the result that P consists of only nondominated solutions. The number of Pareto-optimal solutions in Pi that are not dominated by any other solutions in P is calculated by (6.2). The ratio of Pareto-optimal solutions in Pi that are not dominated by any other solutions in P is calculated by (6.3):
(6.2)Numbernondominatedsolutions=Pi-{X∈Pi∣∃Y∈P:Y≺X},(6.3)Rpos(Pi)=|Pi-{X∈Pi∣∃Y∈P:Y≺X}||Pi|,
where Y≺X means that the solution X is dominated by the solution Y. The higher the ratio Rpos(Pi) is, the better the solution set Pi [22].We calculated the three standards by running the three algorithms for all of the fifteen instances over 10 runs. The computational results are shown in Table3. As can be observed from the table, the average numbers of Pareto-optimal solutions are approximately equal with NSGA-II and MOSA, and GA_WS is inferior to NSGA-II or MOSA for all of the instances expect for three instances. The comparison of the three algorithms with respect to the average number of nondominated solutions shows that NSGA-II performs the best among the three algorithms for all of the instances expect for instance 10. The average ratio of Pareto-optimal solutions on NSGA-II changes between 81% and 100%. This ratio is between 29% and 100% on MOSA, and the ratio changes between 9% and 100% on GA_WS. These results suggest that NSGA-II tends to find more solutions with higher quality than the other two algorithms. This advantage grows when the problem size becomes larger.Table 3
Comparison of the performance of the three algorithms.
Standard(1)
Standard(2)
Standard(3)
NSGA-II
MOSA
GA_WS
NSGA-II
MOSA
GA_WS
NSGA-II
MOSA
GA_WS
Problem set 1
Instance 1
10
10
10
10
10
10
1.00
1.00
1.00
Instance 2
8.7
8.7
8.6
8.7
8.3
7.5
1.00
0.96
0.87
Instance 3
6
6
6
6
6
6
1.00
1.00
1.00
Problem set 2
Instance 4
13
13
12.6
13
12.8
10.3
1.00
0.98
0.82
Instance 5
12
11.7
10.3
12
11.7
10.3
1.00
1.00
1.00
Instance 6
25.2
23
18.5
21.6
19.2
17.1
0.86
0.83
0.93
Problem set 3
Instance 7
25.4
23.3
19.4
25.4
21.5
11.8
1.00
0.93
0.62
Instance 8
27.3
29.7
31.7
26.3
11.8
10.6
0.96
0.41
0.34
Instance 9
19.4
16
15.7
19.4
15.6
11.3
1.00
0.98
0.72
Problem set 4
Instance 10
8
7.7
6.7
6.4
6.6
6
0.81
0.87
0.89
Instance 11
16
13.9
9.5
15.2
9.1
6.5
0.95
0.66
0.69
Instance 12
31
31.2
22.6
30.7
21.3
2
0.99
0.68
0.09
Problem set 5
Instance 13
30.5
48.7
47.6
29.3
13.6
4.4
0.96
0.29
0.09
Instance 14
18.9
21.1
17.6
18.4
9.6
13.7
0.97
0.45
0.78
Instance 15
14.6
16.3
12.6
14.6
11.7
5.4
1.00
0.74
0.45The computation times on NSGA-II, MOSA, and GA_WS for all of the fifteen instances are shown in Table4. It can be observed from this table that the computational time for GA_WS is the shortest time among all of the three algorithms. This result occurs because that once MOSA generates a new solution, it is compared with each solution in the Pareto-optimal set to determine whether it can be accepted. Once it is accepted, the Pareto-optimal set is updated. These operators increase the computational time of MOSA. GA_WS updates the Pareto-optimal set every generation. Although NSGA-II obtains the Pareto-optimal set only from the population of the last generation, it applies a fast nondominated sorting approach and a crowded-comparison approach to evaluate the individuals of every generation, which makes its computational time longer than GA_WS. The computational time in MOSA is longer than in NSGA-II for all of the instances except for the instances in Problem set 5. When the size of the instances becomes larger in Problem set 5, the computation times on NSGA-II become longer than for MOSA.Table 4
Comparison of the three algorithms in terms of the computational time.
CPU times (s)
NSGA-II
MOSA
GA_WS
Problem set 1
Instance 1
51.34
65.05
43.28
Instance 2
59.79
68.87
44.78
Instance 3
27.66
32.80
25.78
Problem set 2
Instance 4
161.37
214.64
114.91
Instance 5
63.92
83.50
62.11
Instance 6
57.93
59.37
51.25
Problem set 3
Instance 7
189.52
238.05
136.74
Instance 8
89.89
92.35
81.01
Instance 9
104.11
107.26
86.83
Problem set 4
Instance 10
134.88
144.69
121.89
Instance 11
138.22
139.36
124.94
Instance 12
154.35
168.12
135.13
Problem set 5
Instance 13
220.63
207.66
200.92
Instance 14
227.20
213.55
202.63
Instance 15
244.02
233.79
212.03It can also be observed from Table4 that, except for the instance 4 and instance 7, the computational time for each algorithm increases as the size of the problem set increases. Instance of problem set 1 contains 1188 binary variables and 106 constraints. Instance of problem set 2 contains 2140 binary variables and 140 constraints. Instance of problem set 3 contains 3195 binary variables and 170 constraints. Instance of problem set 4 contains 4532 binary variables and 204 constraints. Instance of problem set 5 contains 6662 binary variables and 244 constraints. As the number of variables and constraints increases, the computational time for each of the algorithms increases accordingly.


### 6.3. Model Experiments with Sensitivity Analysis

Sensitivity experiments were conducted on the maximum capacity of the repair centresmk, the expected cycle time te, and the service coverage of the collection points l to see how these parameters affect the objective function values and the Pareto-optimal set. We set mk to be 1000, 1500, and 2000; te to be 25, 30, and 35; l to be 20, 25, and 30, respectively. For each instance, the Pareto-optimal sets obtained by the three algorithms under different parameter values were compared. The results of sensitivity experiments conducted on the three parameters are similar for the three algorithms. For the limited space, we only give some illustrations of the changes in the three parameters for the Pareto-optimal set obtained by NSGA-II.The experiments show that the maximum capacity of the repair centres affects only the outcome of the Pareto-optimal set of large problem sets, including Problem set 3, Problem set 4, and Problem set 5. For instances of these problem sets, as the maximum capacity of the repair centres increases, the total reverse logistics cost decreases. Table5 illustrates the change in the maximum capacity of the repair centres on the Pareto-optimal set obtained by NSGA-II. Here,  f1, f2, and f3 represent the three objective function values, respectively. For Problem set 1 and Problem set 2, the Pareto-optimal set appears to be insensitive to changes in the maximum capacity of the repair centres. This scenario can be explained as follows. For Problem set 1 and Problem set 2, the volume of the daily returned products is relatively small, so that a small capacity in the repair centres can meet the demand. Therefore, the increase in the maximum capacity of the repair centres does not affect the cost of the Pareto-optimal solutions. However, when the problem size becomes larger, a larger capacity of the repair centres is required to reduce the cost. If the capacity of the repair centre is still small, then more repair centres are needed and the cost will be higher.Table 5
Sensitivity analysis with respect to the maximum capacity of a repair centre.
m
k
=
1000
m
k
=
1500
m
k
=
2000
f
1
f
2
f
3
f
1
f
2
f
3
f
1
f
2
f
3
548970
432.12
100
535120
432.12
100
532120
432.12
100
571825
168.58
100
555475
168.58
100
555475
168.58
100
575420
131.85
100
561570
131.85
100
558570
131.85
100
555475
235.79
99.07
552475
235.79
99.07
552475
235.79
99.07
542095
0
98.21
541845
0
98.21
538845
0
98.21
539170
189.14
93.46
541445
132.10
98.21
538445
132.10
98.21
542050
163.43
93.43
539095
0
97.21
530695
6.68
93.19
539050
1274.64
93.43
533695
6.68
93.19
532375
0
91.64
535375
0
91.64
533170
7531.29
92.47
529450
500.09
86.90
529450
7531.29
86.90
535375
0
91.64
526450
500.09
85.90
534055
228.44
83.38
532375
0
90.64
528445
0
78.81
531445
0
77.30
529450
500.09
86.89
521725
0
72.24
528655
6999.83
72.24
526450
500.09
85.90
528445
0
78.81
521725
0
72.24The experiments show that the expected cycle time and the service coverage of collection points affect the outcome of the Pareto-optimal set of almost all of the instances. As the expected cycle time increases, the total tardiness of the cycle time decreases, and the cost of Pareto-optimal solutions decreases when the expected cycle time increases, to some extent. Table6 illustrates the change in the expected cycle time for the Pareto-optimal set obtained by NSGA-II. The coverage rate ascends obviously with an increase in the service coverage, and the cost also decreases when the expected cycle time increases, to some extent. Table 7 illustrates the change in the service coverage of the Pareto-optimal set obtained by NSGA-II. This phenomenon illustrates the trade-off among the objective functions, as we expected. When the consumers have lower requirements for the service coverage and the expected cycle time, some of the cost can be saved. In return, the company must pay more money to meet the higher requirements of consumers.Table 6
Sensitivity analysis with respect to the expected cycle time.
t
e
=
25
t
e
=
30
t
e
=
35
f
1
f
2
f
3
f
1
f
2
f
3
f
1
f
2
f
3
439355
2421.61
100
439355
777.64
100
407345
0
100
411335
3266.69
100
411335
1246.69
100
407135
0
97.77
407345
4350.20
100
407345
1660.2
100
402235
78.21
97.77
428755
2421.61
97.77
428755
777.64
97.77
402945
0
89.63
407135
4236.99
97.77
407135
1616.99
97.77
398750
0
77.59
402235
4755.20
97.77
402235
1915.20
97.77
411590
696.78
96.54
411590
0
96.54
406935
3266.69
89.63
406935
1246.69
89.63
402945
4350.20
89.63
402945
1660.20
89.63
419225
0
86.83
407190
0
86.17
407190
696.78
86.17
398750
43.91
77.59
398750
198.91
77.59
404875
0
51.51
413475
0
77.07
400375
0
41.14
409075
0
66.70
395925
963.65
41.14
404875
0
51.51
400375
0
41.14
395925
1118.65
41.14Table 7
Sensitivity analysis with respect to the service coverage.
l
=
20
l
=
25
l
=
30
f
1
f
2
f
3
f
1
f
2
f
3
f
1
f
2
f
3
346220
0
100
342995
0
100
332325
0
100
342230
76.45
100
331395
75.13
100
331395
75.13
100
339535
1652.45
100
332325
0
97.93
328925
7254.31
100
342995
0
97.93
330375
457.11
97.93
330375
457.11
100
341695
457.11
97.52
328470
535.17
85.59
329570
1548.32
97.52
331395
75.13
94.06
326020
6224.36
85.59
327815
9259.53
95.62
336125
0
91.99
326645
507.98
58.04
328470
535.17
94.47
332325
0
91.59
328595
50.87
58.04
326020
6224.36
94.47
330375
457.11
91.59
325520
1277.29
51.93
326645
507.98
68.88
328470
535.17
75.85
324325
0
33.26
328595
50.87
68.88
326020
6224.36
75.85
324325
0
39.77
328595
50.87
55.56
326645
507.98
55.56


## 6.1. Date Generation

To compare the performance of the three algorithms, five sets of fifteen randomly generated instances were created in Table1. The number of daily returned products by each customer zone was randomly generated to be between 10 and 50. The round trip transportation time between the collection point j  and the repair centre k was computed as follows: tjk=djk×0.6. Other parameters of the model are shown in Table 2.Table 1
Parameters of five problem sets.
Problem sets
1
2
3
4
5
Number of customer zones
50
60
80
100
120
Number of potential collection points
20
30
35
40
50
Number of potential repair centres
8
10
10
12
12
Location scale of all points
0~100
0~100
0~150
0~150
0~200Table 2
Parameters of the model.
Parameter
Index
Value
Annual renting cost of the collection points
a
200
Handling cost per unit product
h
0.1
Capacity of a repair centre
m
k
2000
Minimum number of open collection points
cp
5
Minimum number of established repair centres
rc
1
Average working hours to repair a returned product
t
m
10
α
1
0.8
Discount rate with respect to the shipping volume
α
2
0.6
p
1
200 units
p
2
400 units
β
1
1.1
Penalty rate with respect to the shipping distance
β
2
1.2
q
1
25
q
2
60
Working days per year
w
250
Cost of establishing a repair centre
e
k
3000
Service coverage
l
25
Unit standard transportation cost
E
1
Expected cycle time
t
e
30


## 6.2. Comparison of the Three Algorithms

To compare the Pareto-optimal set that is obtained by three algorithms on the same basis, the number of solutions searched is taken to be the stopping criterion. According to the different sizes of the five Problem sets, the number of solutions searched is set to be 120,000, 122,000, 128,000, 130,000, and 132,000. Because both GA_WS and NSGA-II are based on a GA, the same parameters are set for both. Based on extensive experiments, the parameters are as follows: the population size = 40, the crossover rate = 0.96, and the mutationrate=0.07. Because the population size of GA_WS or NSGA-II is 40, the maximum numbers of generations for GA_WS and NSGA-II on the five Problem sets are 3,000, 3,050, 3,200, 3,250, and 3,300, respectively. For MOSA, the initial temperature Ts is calculated from the following
(6.1)Ts=ΔSUMmaxlnProb,
where ΔSUMmax is the maximum value of the sum of differences between the maximum and the minimum of the three objective values, which are chosen from the objective values of 30 randomly generated neighbours of the initial solution. The value of Prob is set to 0.95 to ensure that a feasible inferior solution is accepted with a probability of 0.95. Other parameters are as follows: iter_max=1000,un_max=50, iter_num=40, and b=0.95. The search process terminates when the total number of evaluated solutions reaches the set value. These three algorithms are all coded in the C++ programming language in the VC++6.0 environment and are executed on a Dell Intel Core 2 Duo computer with a speed of 2.10 GHz and with 2.00 GB of memory.The Pareto-optimal set is the set of Pareto-optimal solutions that consists of all decision vectors for which the corresponding objective vectors cannot be improved in a given dimension without worsening another solution [39]. To evaluate the performance of the three algorithms, we adopted the following standards. Standard (1), the average number of Pareto-optimal solutions; Standard (2), the average number of nondominated solutions; Standard (3), the average ratio of the Pareto-optimal solutions. Due to the stochastic nature of the suggested algorithms, these standards were obtained by the three algorithms over 10 runs. Standard (2) and Standard (3) were calculated in the following manner. Let P1, P2, and P3 be the sets of Pareto-optimal solutions that are obtained from one run of GA_WS, NSGA-II, and MOSA, respectively, and let P be the union of the sets of Pareto-optimal solutions (i.e., P=P1∪P2∪P3), with the result that P consists of only nondominated solutions. The number of Pareto-optimal solutions in Pi that are not dominated by any other solutions in P is calculated by (6.2). The ratio of Pareto-optimal solutions in Pi that are not dominated by any other solutions in P is calculated by (6.3):
(6.2)Numbernondominatedsolutions=Pi-{X∈Pi∣∃Y∈P:Y≺X},(6.3)Rpos(Pi)=|Pi-{X∈Pi∣∃Y∈P:Y≺X}||Pi|,
where Y≺X means that the solution X is dominated by the solution Y. The higher the ratio Rpos(Pi) is, the better the solution set Pi [22].We calculated the three standards by running the three algorithms for all of the fifteen instances over 10 runs. The computational results are shown in Table3. As can be observed from the table, the average numbers of Pareto-optimal solutions are approximately equal with NSGA-II and MOSA, and GA_WS is inferior to NSGA-II or MOSA for all of the instances expect for three instances. The comparison of the three algorithms with respect to the average number of nondominated solutions shows that NSGA-II performs the best among the three algorithms for all of the instances expect for instance 10. The average ratio of Pareto-optimal solutions on NSGA-II changes between 81% and 100%. This ratio is between 29% and 100% on MOSA, and the ratio changes between 9% and 100% on GA_WS. These results suggest that NSGA-II tends to find more solutions with higher quality than the other two algorithms. This advantage grows when the problem size becomes larger.Table 3
Comparison of the performance of the three algorithms.
Standard(1)
Standard(2)
Standard(3)
NSGA-II
MOSA
GA_WS
NSGA-II
MOSA
GA_WS
NSGA-II
MOSA
GA_WS
Problem set 1
Instance 1
10
10
10
10
10
10
1.00
1.00
1.00
Instance 2
8.7
8.7
8.6
8.7
8.3
7.5
1.00
0.96
0.87
Instance 3
6
6
6
6
6
6
1.00
1.00
1.00
Problem set 2
Instance 4
13
13
12.6
13
12.8
10.3
1.00
0.98
0.82
Instance 5
12
11.7
10.3
12
11.7
10.3
1.00
1.00
1.00
Instance 6
25.2
23
18.5
21.6
19.2
17.1
0.86
0.83
0.93
Problem set 3
Instance 7
25.4
23.3
19.4
25.4
21.5
11.8
1.00
0.93
0.62
Instance 8
27.3
29.7
31.7
26.3
11.8
10.6
0.96
0.41
0.34
Instance 9
19.4
16
15.7
19.4
15.6
11.3
1.00
0.98
0.72
Problem set 4
Instance 10
8
7.7
6.7
6.4
6.6
6
0.81
0.87
0.89
Instance 11
16
13.9
9.5
15.2
9.1
6.5
0.95
0.66
0.69
Instance 12
31
31.2
22.6
30.7
21.3
2
0.99
0.68
0.09
Problem set 5
Instance 13
30.5
48.7
47.6
29.3
13.6
4.4
0.96
0.29
0.09
Instance 14
18.9
21.1
17.6
18.4
9.6
13.7
0.97
0.45
0.78
Instance 15
14.6
16.3
12.6
14.6
11.7
5.4
1.00
0.74
0.45The computation times on NSGA-II, MOSA, and GA_WS for all of the fifteen instances are shown in Table4. It can be observed from this table that the computational time for GA_WS is the shortest time among all of the three algorithms. This result occurs because that once MOSA generates a new solution, it is compared with each solution in the Pareto-optimal set to determine whether it can be accepted. Once it is accepted, the Pareto-optimal set is updated. These operators increase the computational time of MOSA. GA_WS updates the Pareto-optimal set every generation. Although NSGA-II obtains the Pareto-optimal set only from the population of the last generation, it applies a fast nondominated sorting approach and a crowded-comparison approach to evaluate the individuals of every generation, which makes its computational time longer than GA_WS. The computational time in MOSA is longer than in NSGA-II for all of the instances except for the instances in Problem set 5. When the size of the instances becomes larger in Problem set 5, the computation times on NSGA-II become longer than for MOSA.Table 4
Comparison of the three algorithms in terms of the computational time.
CPU times (s)
NSGA-II
MOSA
GA_WS
Problem set 1
Instance 1
51.34
65.05
43.28
Instance 2
59.79
68.87
44.78
Instance 3
27.66
32.80
25.78
Problem set 2
Instance 4
161.37
214.64
114.91
Instance 5
63.92
83.50
62.11
Instance 6
57.93
59.37
51.25
Problem set 3
Instance 7
189.52
238.05
136.74
Instance 8
89.89
92.35
81.01
Instance 9
104.11
107.26
86.83
Problem set 4
Instance 10
134.88
144.69
121.89
Instance 11
138.22
139.36
124.94
Instance 12
154.35
168.12
135.13
Problem set 5
Instance 13
220.63
207.66
200.92
Instance 14
227.20
213.55
202.63
Instance 15
244.02
233.79
212.03It can also be observed from Table4 that, except for the instance 4 and instance 7, the computational time for each algorithm increases as the size of the problem set increases. Instance of problem set 1 contains 1188 binary variables and 106 constraints. Instance of problem set 2 contains 2140 binary variables and 140 constraints. Instance of problem set 3 contains 3195 binary variables and 170 constraints. Instance of problem set 4 contains 4532 binary variables and 204 constraints. Instance of problem set 5 contains 6662 binary variables and 244 constraints. As the number of variables and constraints increases, the computational time for each of the algorithms increases accordingly.


## 6.3. Model Experiments with Sensitivity Analysis

Sensitivity experiments were conducted on the maximum capacity of the repair centresmk, the expected cycle time te, and the service coverage of the collection points l to see how these parameters affect the objective function values and the Pareto-optimal set. We set mk to be 1000, 1500, and 2000; te to be 25, 30, and 35; l to be 20, 25, and 30, respectively. For each instance, the Pareto-optimal sets obtained by the three algorithms under different parameter values were compared. The results of sensitivity experiments conducted on the three parameters are similar for the three algorithms. For the limited space, we only give some illustrations of the changes in the three parameters for the Pareto-optimal set obtained by NSGA-II.The experiments show that the maximum capacity of the repair centres affects only the outcome of the Pareto-optimal set of large problem sets, including Problem set 3, Problem set 4, and Problem set 5. For instances of these problem sets, as the maximum capacity of the repair centres increases, the total reverse logistics cost decreases. Table5 illustrates the change in the maximum capacity of the repair centres on the Pareto-optimal set obtained by NSGA-II. Here,  f1, f2, and f3 represent the three objective function values, respectively. For Problem set 1 and Problem set 2, the Pareto-optimal set appears to be insensitive to changes in the maximum capacity of the repair centres. This scenario can be explained as follows. For Problem set 1 and Problem set 2, the volume of the daily returned products is relatively small, so that a small capacity in the repair centres can meet the demand. Therefore, the increase in the maximum capacity of the repair centres does not affect the cost of the Pareto-optimal solutions. However, when the problem size becomes larger, a larger capacity of the repair centres is required to reduce the cost. If the capacity of the repair centre is still small, then more repair centres are needed and the cost will be higher.Table 5
Sensitivity analysis with respect to the maximum capacity of a repair centre.
m
k
=
1000
m
k
=
1500
m
k
=
2000
f
1
f
2
f
3
f
1
f
2
f
3
f
1
f
2
f
3
548970
432.12
100
535120
432.12
100
532120
432.12
100
571825
168.58
100
555475
168.58
100
555475
168.58
100
575420
131.85
100
561570
131.85
100
558570
131.85
100
555475
235.79
99.07
552475
235.79
99.07
552475
235.79
99.07
542095
0
98.21
541845
0
98.21
538845
0
98.21
539170
189.14
93.46
541445
132.10
98.21
538445
132.10
98.21
542050
163.43
93.43
539095
0
97.21
530695
6.68
93.19
539050
1274.64
93.43
533695
6.68
93.19
532375
0
91.64
535375
0
91.64
533170
7531.29
92.47
529450
500.09
86.90
529450
7531.29
86.90
535375
0
91.64
526450
500.09
85.90
534055
228.44
83.38
532375
0
90.64
528445
0
78.81
531445
0
77.30
529450
500.09
86.89
521725
0
72.24
528655
6999.83
72.24
526450
500.09
85.90
528445
0
78.81
521725
0
72.24The experiments show that the expected cycle time and the service coverage of collection points affect the outcome of the Pareto-optimal set of almost all of the instances. As the expected cycle time increases, the total tardiness of the cycle time decreases, and the cost of Pareto-optimal solutions decreases when the expected cycle time increases, to some extent. Table6 illustrates the change in the expected cycle time for the Pareto-optimal set obtained by NSGA-II. The coverage rate ascends obviously with an increase in the service coverage, and the cost also decreases when the expected cycle time increases, to some extent. Table 7 illustrates the change in the service coverage of the Pareto-optimal set obtained by NSGA-II. This phenomenon illustrates the trade-off among the objective functions, as we expected. When the consumers have lower requirements for the service coverage and the expected cycle time, some of the cost can be saved. In return, the company must pay more money to meet the higher requirements of consumers.Table 6
Sensitivity analysis with respect to the expected cycle time.
t
e
=
25
t
e
=
30
t
e
=
35
f
1
f
2
f
3
f
1
f
2
f
3
f
1
f
2
f
3
439355
2421.61
100
439355
777.64
100
407345
0
100
411335
3266.69
100
411335
1246.69
100
407135
0
97.77
407345
4350.20
100
407345
1660.2
100
402235
78.21
97.77
428755
2421.61
97.77
428755
777.64
97.77
402945
0
89.63
407135
4236.99
97.77
407135
1616.99
97.77
398750
0
77.59
402235
4755.20
97.77
402235
1915.20
97.77
411590
696.78
96.54
411590
0
96.54
406935
3266.69
89.63
406935
1246.69
89.63
402945
4350.20
89.63
402945
1660.20
89.63
419225
0
86.83
407190
0
86.17
407190
696.78
86.17
398750
43.91
77.59
398750
198.91
77.59
404875
0
51.51
413475
0
77.07
400375
0
41.14
409075
0
66.70
395925
963.65
41.14
404875
0
51.51
400375
0
41.14
395925
1118.65
41.14Table 7
Sensitivity analysis with respect to the service coverage.
l
=
20
l
=
25
l
=
30
f
1
f
2
f
3
f
1
f
2
f
3
f
1
f
2
f
3
346220
0
100
342995
0
100
332325
0
100
342230
76.45
100
331395
75.13
100
331395
75.13
100
339535
1652.45
100
332325
0
97.93
328925
7254.31
100
342995
0
97.93
330375
457.11
97.93
330375
457.11
100
341695
457.11
97.52
328470
535.17
85.59
329570
1548.32
97.52
331395
75.13
94.06
326020
6224.36
85.59
327815
9259.53
95.62
336125
0
91.99
326645
507.98
58.04
328470
535.17
94.47
332325
0
91.59
328595
50.87
58.04
326020
6224.36
94.47
330375
457.11
91.59
325520
1277.29
51.93
326645
507.98
68.88
328470
535.17
75.85
324325
0
33.26
328595
50.87
68.88
326020
6224.36
75.85
324325
0
39.77
328595
50.87
55.56
326645
507.98
55.56


## 7. Conclusions

In this paper, we presented a multiobjective integer nonlinear programming model for a three-echelon reverse logistic network design problem. This model considered not only the traditional cost factor but also the service level, which was represented by the total tardiness of the cycle time and the coverage of consumer zones. The model can help a computer manufacturer decide the optimal number and location of collection points and repair centres and the transportation arrangement of the returned products from the customer zones to the repair centres through the collection points after a trade-off of the total cost of the reverse logistics and the service quality level.Because it is an NP-hard problem with the property of multiobjective, NSGA-II was adopted for the proposed model. To evaluate the performance of NSGA-II, GA_WS and MOSA were also applied. The performance of the three algorithms was compared on five sets of fifteen randomly generated instances. The comparative analysis showed that NSGA-II and MOSA outperformed GA_WS in terms of the average numbers of Pareto-optimal solutions, and NSGA-II tended to find solutions with the highest quality among the three algorithms. Finally, three key parameters of the models were tested, including the maximum capacity of the repair centres, the expected cycle time, and the service coverage of the collection points. The analysis results showed that the maximum capacity of the repair centres only affected the outcome of the Pareto-optimal set of large problem sets. The expected cycle time and the service coverage of the collection points affected the outcome of the Pareto-optimal set of almost all of the instances.Note that the study in this paper is based on the assumption that the location/allocation plan covers a planning horizon within which no substantial changes are incurred in terms of customer demands and the transportation infrastructure. However, this scenario is often not the case in many practical circumstances. Thus, the dynamic multiperiod reverse logistics network design problem must be addressed in further studies.



---

*Source: 928620-2012-10-16.xml*