# An Improved Differential Evolution Algorithm for a Multicommodity Location-Inventory Problem with False Failure Returns

Authors: Congdong Li; Hao Guo; Ying Zhang; Shuai Deng; Yu Wang
Journal: Complexity
 (2018)

---

## Abstract

Customer returns are a common phenomenon in many industries, and they have a significant impact on business organizations and their supply chains. False failure returns are returned products that have no functional or cosmetic defects, and they represent a large body of customer returns in practice. In this paper, we develop a mixed-integer nonlinear programming model to study a multicommodity location-inventory problem in a forward-reverse logistics network. This model minimizes the total cost in this network by considering false failure returns, and it also considers many real-world business scenarios in forward and reverse logistics flows. Moreover, we design a new heuristic approach to solve the model efficiently. Finally, numerical experiments are conducted to validate our solution approach and provide meaningful managerial insights.


---

## Body

## 1. Introduction

Reverse logistics is an emerging topic that has attracted many research attentions because of the great importance of customer returns in today’s business. Nowadays, the percentage of customer returns in total sales varies from 5% to 9% for durable goods, and it is up to 35% for fashion apparel products ‎[1]. Since customer returns are more frequent in e-commerce, return rate in online shopping has exceeded 18% of total sales, and especially, it reaches 74% for fashion products ‎[2]. In the retail industry, the value of customer products has exceeded $100 billion per year in the United States ‎[3]. Although many products are returned because of quality issues, a lot of items that have no defects can still be returned for various reasons such as installation difficulties, bad user experience, or price competition. Those items represent a large body of customer returns, and they are defined as “false failure returns” ‎[4].Forward-reverse logistics network consist of forward flows that are directed from suppliers to customers and reverse flows of returned products from customers to suppliers. Since network design can contribute to reducing 60% the cost of a firm ‎[5], Supply chain network design is an important strategic decision in supply chain management, and it is a big challenge to business managers. For example, many big companies such as Kodak, Xerox, and General Motors have developed their reverse logistics systems and those efforts are very successful ‎[6]. Currently, many business decisions are made jointly to improve the efficiency of business organizations and their supply chains. Location-inventory problem (LIP) ‎[7] is such a problem that integrates strategical facility location decisions and tactical inventory management decisions, and it can be studied by incorporating many real-world business scenarios to improve the performance of a logistics system.In this paper, we study a multicommodity location-inventory problem (MCLIP) in a forward-reverse logistics network in which the reverse flow consists of false failure returns. This problem is formulated as a mixed-integer nonlinear programming model whose objective function is to minimize the total cost in such a network, and it incorporates many business decisions and processes such as facility locations, return processing, and order fulfillment. More specifically, the solution to this problem will determine the number and locations of hybrid distribution-collection centers (HDCCs), how to use those HDCCs to fulfill orders and collect returns, and how to replenish inventories in those HDCCs. Due to the complexity of LIPs, we also design an improved differential evolution algorithm (IDE) to solve the MCLIP efficiently. Computational results are presented to validate the solution approach, and research implications and managerial insights are also discussed.The rest of this paper is organized as follows: In Section2, we review related research works in the literature. In Section 3, we describe the research problem and present the mathematical model. In Section 4, we propose a new heuristic approach to solve the problem efficiently. In Section 5, we show computational results and discuss managerial insights. Section 6 concludes the paper and provides directions for future research.


## 2. Literature Review

The design of logistics networks involves several core components such as facility location and inventory management. LIPs have been extensively studied because of their great importance in practice, and Farahani et al. ‎[8] provide a comprehensive review of such works in the literature. In this section, we review the research works on LIPs from the perspective of forward and reverse logistics, respectively.In the literature, most LIPs are studied under forward logistics networks. For example, Daskin et al. ‎[7] study a LIP with the risk-pooling effect and design a Lagrangian relaxation algorithm as the solution approach. Shen et al. ‎[9] use the same integer programming model as that in Daskin et al.‎[7] and convert it into a set-covering integer programming model and then solve it by column generation. Moreover, LIPs are studied by incorporating many business scenarios and processes. For example, Ozsen et al. ‎[10] consider capacitated warehouses and the impact of multisourcing in addition to Shen et al. ‎[9]. Amiri-Aref et al. ‎[11] study multisourcing location-inventory problem with stochastic demand. Silva et al. ‎[12] and Qu et al. ‎[13] study the stochastic location-inventory problems under different replenishment policies. Mousavi et al. ‎[14] study a seasonal multiple-product location allocation and inventory control problem in a planning horizon consisting of multiple periods. Since business environments are always dynamic, demand uncertainty ‎[15–17] and different risks ‎[18] are also considered.It is emerging to study LIPs in the reverse logistics. Since reverse logistics is an important part of closed-loop supply chains (CLSCs), those works are usually conducted with CLSCs. For example, Diabat et al. ‎[19] study a closed-loop location-inventory problem by considering spare parts, Asl-Najafi et al. ‎[20] study a dynamic closed-loop location-inventory problem under disruption risk, and Kaya et al. ‎[21] and Guo et al. ‎[22] extend such studies by incorporating pricing and vehicle routing decisions in CLSCs, respectively. Since returned products are usually sold to secondary markets in practice, LIPs are also studied by considering CLSCs and secondary markets ‎[23, 24].Although LIPs have been extensively studied from the perspective of reverse logistics, those works are still very limited because the types of customer returns are rarely considered. This paper studies a LIP in a forward-reverse logistics network in which the reverse logistics flow only consists of false failure returns. To the best of our knowledge, this is the first paper that considers false failure returns in the LIP literature, and it makes a significant contribution to the literature because of the great importance of false failure returns in practice.


## 3. Mathematical Model

### 3.1. Problem Description

We consider a three-echelon supply chain network in which a manufacturing plant ships multiple commodities to customer zones via a set of facilities. Moreover, returned products will be collected from customers and shipped to the same facilities. Therefore, those facilities are hybrid distribution-collection centers (HDCCs) [25, 26], and they have many advantages such as cost saving and pollution reduction by improving the efficiency of workflows. In this network, we assume that only false failure returns will be collected. Since those returns do not have any quality issues, they will not be sent to the plant for repair, and instead, they will be resold in the market after some simple processes (e.g., inspection, repackaging, etc.) in HDCCs. Figure 1 shows the forward-reverse logistics network under study. In this network, forward flows will be unidirectional from the plant to HDCCs and then from HDCCs to customer zones, but reverse flows will only be directed from customer zones to HDCCs because false failure returns will not be sent to the plant for repair.Figure 1
A forward-reverse logistics network with false failure returns.This paper studies a multicommodity location-inventory problem (MCLIP) in the network shown in Figure1, and we make the following assumptions to simplify the problem: (i) customer zones are predetermined and fixed; (ii) the returned products will be the same as new items when they are resold to the market; (iii) the distances between HDCCs and customer zones are Euclidean distances. In practice, this problem can be solved to answer the following questions:(1)
What are the best locations for HDCCs in this network?(2)
How to assign those HDCCs to customer zones in the forward and reverse networks, respectively?(3)
How to optimize inventory replenishment decisions at HDCCs?


### 3.2. Notations

Sets W: set of commodities, indexed by w; R: set of candidate HDCC locations, indexed byr; I: set of customer zones, indexed byi.Parameters a r: fixed (yearly) cost of building a HDCC at locationr, for eachr∈R; b r w: fixed administrative and handling cost of placing an order for product w at HDCCr, for eachr∈R and w∈W; c r w: basic transportation cost per shipment of product w from the plant to HDCCr, for eachr∈Rand w∈W; e r w: additional transportation cost per unit of product w from the plant to HDCCr, for eachr∈R and w∈W; d r i w: shipping cost per unit of product w from HDCCr to customer zonei, for eachr∈R,i∈I and w∈W; f r w: repackaging/refurbishing cost per unit of returned product w at HDCCr, for eachr∈R and w∈W; g r w: inspection cost per unit of returned product w at HDCCr, for eachr∈R and w∈W; h r w: (yearly) holding cost per unit of product w at HDCCr, for eachr∈R and w∈W; q w: return rate of product w, for each w∈W; µ i w: mean (daily) demand of product w at customer zonei, for eachi∈Iand w∈W; σ i w 2: variance of (daily) demand of product w at customer zonei, for eachi∈I and w∈W; α: desired percentage of market demand satisfied; z α: standard normal deviation such that P(z≤zα)=α; L: order lead time in days at HDCCs; λ: working days per year; D r w: mean (yearly) demand of product w at HDCCr, for eachr∈R and w∈W.Decision Variables Q r w: order quantity of product w per order at HDCCr; X r = 1 if opening one HDCC at locationr, and 0 otherwise, for eachr∈R; Y r i w = 1 if HDCCr fulfills demand at customer zonei, and 0 otherwise, for eachr∈R,i∈I and w∈W; Z i r w = 1 if HDCCr collects returns from customer zonei, and 0 otherwise, for eachr∈R,i∈I and w∈W.


### 3.3. The Model

In this paper, the research problem is formulated by a mixed-integer nonlinear programming model. The objective of this model is to minimize the total cost in the forward-reverse logistics network, and it comprises (i) location costs, (ii) shipping costs from HDCCs to customer zones, (iii) inventory costs including working inventory and safety stock costs, and (iv) return costs.


#### 3.3.1. Location Cost

Obviously, the total location cost (CL) can be written as(1)CL=∑r∈RarXr.


#### 3.3.2. Shipping Cost from HDCCs to Customer Zones

The total shipping cost (CT) from HDCCs to customer zones in the forward logistics flow can be expressed as follows:(2)CT=λ∑r∈R∑w∈W∑i∈IμiwdriwYriw.


#### 3.3.3. Inventory Cost

Inventory strategy plays an important role in supply chain management. In this study, we assume that HDCCs order new products from the plant by using an (R,Q) policy with type I service, and the total inventory cost includes working inventory and safety stock costs. The (R,Q) model can be approximated by an economic order quantity (EOQ) model [9], and Axsater [27] shows that the maximum relative error introduced by using the EOQ model instead of the (Q,r) policy is 0.118. Therefore, we use an EOQ model to optimize the inventory strategy in this study, and the order frequency and quantity at a HDCC will be calculated by using the total demand of the customer zones served by this HDCC.In general, the working inventory cost includes(1) fixed cost of placing orders, (2) shipping cost from the plant to HDCCs, and (3) holding cost of the working inventory. The individual cost terms are formulated as follows.(1) Order Cost. The annual fixed cost of placing orders for product w from HDCCr to the plant is given by brwDrw/Qrw.(2) Shipping Cost from the Plant to HDCCs. Since false failure returns will be resold to the market after they are processed in HDCCs, shipping quantities from the plant to HDCCs will equal the difference between market demands and return volumes. Therefore, the total shipping cost from the plant to HDCCr per year is given by(3)crwDrwQrw+λerw∑i∈IμiwYriw-∑i∈IqwμiwZirw.(3) Holding Cost of the Working Inventory. In this study, the working inventory comes from two sources which are the new products ordered from the plant and the returned products collected from customer zones. Therefore, the forward and reverse logistics flows will be integrated at HDCCs, and the average inventory of product w at HDCCrduring periodT can be expressed as follows:(4)1T∫0TSrw-Srw-qwSrwTdt=Srw-12θrwT=Qrw1-qw-1-qw2Srw=1+qw21-qwQrw,where Srw is the mean daily demand of product w at HDCC r. Consequently, the holding cost of the working inventory is given by (5)∑i∈I1+qwhrw21-qwQrw.Using Eppen’s risk-pooling result [28], the amount of safety stock to ensure that stockouts occur with a probability of α or less is zαL∑i∈Iσiw2Yriw, and the safety stock cost at HDCCr is given by hrwzαL∑i∈Iσiw2Yriw.Therefore, the total annual inventory cost (CI) of product w at HDCCr is given by (6)CI=brwDrwQrw+crwDrwQrw+λerw∑i∈IμiwYriw-∑i∈IqwμiwZirw+∑i∈I1+qwhrw21-qwQrw+hrwzαL∑i∈Iσiw2Yriw.


#### 3.3.4. Return Cost

Returns costs are incurred by customer returns, and they include inspection cost, repackaging cost, and shipping cost from customer zones to HDCCs.(1) Inspection Cost.Returned products will be inspected at HDCCs, and the total inspection cost is given by λ∑r∈R∑w∈W∑i∈IqwμiwgrwZriw.(2) Repackaging Cost. The returned products also need to be repackaged before they are resold to the market, and the total repackaging cost is given by λ∑r∈R∑w∈W∑i∈IqwμiwfrwZriw.(3) Shipping Cost from Customer Zones to HDCCs.The total shipping cost from customer zones to HDCCs is given by λ∑r∈R∑w∈W∑i∈IqwμiwdriwZriw.Therefore, the total annual return cost (CR) is given by(7)CR=λ∑r∈R∑w∈W∑i∈Iqwμiwgrw+frw+driwZriw.


#### 3.3.5. Optimization Model

Given the individual cost functions above, the research problem under study can be formulated by a mixed-integer nonlinear program as follows:(8)min⁡Z=CL+CT+CI+CR=∑r∈RarXr+λ∑r∈R∑w∈W∑i∈IμiwdriwYriw+∑r∈R∑w∈WbrwDrwQrw+crwDrwQrw+λerw∑i∈IμiwYriw-∑i∈IqwμiwZirw+∑i∈I1+qwhrw21-qwQrw+hrwzαL∑i∈Iσiw2Yriw+λ∑r∈R∑w∈W∑i∈Iqwμiwgrw+frw+driwZriwsubject to(9)∑r∈RXr≥1;(10)∑r∈RYriw=1,∀i∈I,∀w∈W;(11)∑r∈RZirw=1,∀i∈I,∀w∈W;(12)∑r∈R∑w∈WYriw=W,∀i∈I;(13)∑r∈R∑w∈WZirw=W,∀i∈I;(14)Yriw≤Xr,∀i∈I,∀r∈R,∀w∈W;(15)Zirw≤Xr,∀i∈I,∀r∈R,∀w∈W;(16)∑w∈W∑i∈IYriw≥Xr,∀r∈R;(17)∑i∈IZirw≤M∑i∈IYriw,∀r∈R,∀w∈W;(18)Xr=0,1,∀r∈R;(19)Yriw=0,1,∀i∈I,∀r∈R,∀w∈W;(20)Zirw=0,1,∀i∈I,∀r∈R,∀w∈W;In this model, objective function (8) means that the total cost will be minimized. Constraint (9) means that at least one HDCC will be built. Constraint (10) means that the demand of a product in a customer zone will be fulfilled by only one HDCC. Constraint (11) means that the returns of a product from a customer zone will be collected by only one HDCC. Constraint (12) means that all types of products will be provided to a customer zone by HDCCs in the forward flow. Constraint (13) means that all types of returned products will be collected from a customer zone by HDCCs in the reverse flow. Constraints (14) and (15) mean that a customer zone can be assigned to a HDCC in the forward and reverse logistics flows, respectively, only if this HDCC is built. Constraint (16) means that a HDCC will serve at least one customer zone in the forward logistics network after it is built. Constraint (17) means that a HDCC will collect the returns of a product in the reverse flow only if it fulfills the demand of this product in the forward flow. Constraints (18)–(20) indicate thatXr, Yriw, and Zirw are binary variables.Note that decision variableQrw only appears in the objective function. Since the objective function is convex in Qrw, its optimal value can be obtained by taking the partial derivative of the objective function with respect to it, which gives(21)Qrw∗=∑i∈I21-qwbrw+crwDrw1+qwhrw.Therefore, the objective function can be rewritten as follows:(22)min⁡Z=λ∑r∈R∑w∈W∑i∈IμiwdriwYriw+∑r∈R∑w∈W∑i∈I2λhrw1+qwbrw+crwμiwYriw1-qw+λ∑r∈R∑w∈W∑i∈IerwμiwYriw-qwμiwZirw+∑r∈R∑w∈WhrwzαL∑i∈Iσiw2Yriw+∑r∈RarXr+λ∑r∈R∑w∈W∑i∈Iqwμiwgrw+frw+driwZriw


## 3.1. Problem Description

We consider a three-echelon supply chain network in which a manufacturing plant ships multiple commodities to customer zones via a set of facilities. Moreover, returned products will be collected from customers and shipped to the same facilities. Therefore, those facilities are hybrid distribution-collection centers (HDCCs) [25, 26], and they have many advantages such as cost saving and pollution reduction by improving the efficiency of workflows. In this network, we assume that only false failure returns will be collected. Since those returns do not have any quality issues, they will not be sent to the plant for repair, and instead, they will be resold in the market after some simple processes (e.g., inspection, repackaging, etc.) in HDCCs. Figure 1 shows the forward-reverse logistics network under study. In this network, forward flows will be unidirectional from the plant to HDCCs and then from HDCCs to customer zones, but reverse flows will only be directed from customer zones to HDCCs because false failure returns will not be sent to the plant for repair.Figure 1
A forward-reverse logistics network with false failure returns.This paper studies a multicommodity location-inventory problem (MCLIP) in the network shown in Figure1, and we make the following assumptions to simplify the problem: (i) customer zones are predetermined and fixed; (ii) the returned products will be the same as new items when they are resold to the market; (iii) the distances between HDCCs and customer zones are Euclidean distances. In practice, this problem can be solved to answer the following questions:(1)
What are the best locations for HDCCs in this network?(2)
How to assign those HDCCs to customer zones in the forward and reverse networks, respectively?(3)
How to optimize inventory replenishment decisions at HDCCs?


## 3.2. Notations

Sets W: set of commodities, indexed by w; R: set of candidate HDCC locations, indexed byr; I: set of customer zones, indexed byi.Parameters a r: fixed (yearly) cost of building a HDCC at locationr, for eachr∈R; b r w: fixed administrative and handling cost of placing an order for product w at HDCCr, for eachr∈R and w∈W; c r w: basic transportation cost per shipment of product w from the plant to HDCCr, for eachr∈Rand w∈W; e r w: additional transportation cost per unit of product w from the plant to HDCCr, for eachr∈R and w∈W; d r i w: shipping cost per unit of product w from HDCCr to customer zonei, for eachr∈R,i∈I and w∈W; f r w: repackaging/refurbishing cost per unit of returned product w at HDCCr, for eachr∈R and w∈W; g r w: inspection cost per unit of returned product w at HDCCr, for eachr∈R and w∈W; h r w: (yearly) holding cost per unit of product w at HDCCr, for eachr∈R and w∈W; q w: return rate of product w, for each w∈W; µ i w: mean (daily) demand of product w at customer zonei, for eachi∈Iand w∈W; σ i w 2: variance of (daily) demand of product w at customer zonei, for eachi∈I and w∈W; α: desired percentage of market demand satisfied; z α: standard normal deviation such that P(z≤zα)=α; L: order lead time in days at HDCCs; λ: working days per year; D r w: mean (yearly) demand of product w at HDCCr, for eachr∈R and w∈W.Decision Variables Q r w: order quantity of product w per order at HDCCr; X r = 1 if opening one HDCC at locationr, and 0 otherwise, for eachr∈R; Y r i w = 1 if HDCCr fulfills demand at customer zonei, and 0 otherwise, for eachr∈R,i∈I and w∈W; Z i r w = 1 if HDCCr collects returns from customer zonei, and 0 otherwise, for eachr∈R,i∈I and w∈W.


## 3.3. The Model

In this paper, the research problem is formulated by a mixed-integer nonlinear programming model. The objective of this model is to minimize the total cost in the forward-reverse logistics network, and it comprises (i) location costs, (ii) shipping costs from HDCCs to customer zones, (iii) inventory costs including working inventory and safety stock costs, and (iv) return costs.


### 3.3.1. Location Cost

Obviously, the total location cost (CL) can be written as(1)CL=∑r∈RarXr.


### 3.3.2. Shipping Cost from HDCCs to Customer Zones

The total shipping cost (CT) from HDCCs to customer zones in the forward logistics flow can be expressed as follows:(2)CT=λ∑r∈R∑w∈W∑i∈IμiwdriwYriw.


### 3.3.3. Inventory Cost

Inventory strategy plays an important role in supply chain management. In this study, we assume that HDCCs order new products from the plant by using an (R,Q) policy with type I service, and the total inventory cost includes working inventory and safety stock costs. The (R,Q) model can be approximated by an economic order quantity (EOQ) model [9], and Axsater [27] shows that the maximum relative error introduced by using the EOQ model instead of the (Q,r) policy is 0.118. Therefore, we use an EOQ model to optimize the inventory strategy in this study, and the order frequency and quantity at a HDCC will be calculated by using the total demand of the customer zones served by this HDCC.In general, the working inventory cost includes(1) fixed cost of placing orders, (2) shipping cost from the plant to HDCCs, and (3) holding cost of the working inventory. The individual cost terms are formulated as follows.(1) Order Cost. The annual fixed cost of placing orders for product w from HDCCr to the plant is given by brwDrw/Qrw.(2) Shipping Cost from the Plant to HDCCs. Since false failure returns will be resold to the market after they are processed in HDCCs, shipping quantities from the plant to HDCCs will equal the difference between market demands and return volumes. Therefore, the total shipping cost from the plant to HDCCr per year is given by(3)crwDrwQrw+λerw∑i∈IμiwYriw-∑i∈IqwμiwZirw.(3) Holding Cost of the Working Inventory. In this study, the working inventory comes from two sources which are the new products ordered from the plant and the returned products collected from customer zones. Therefore, the forward and reverse logistics flows will be integrated at HDCCs, and the average inventory of product w at HDCCrduring periodT can be expressed as follows:(4)1T∫0TSrw-Srw-qwSrwTdt=Srw-12θrwT=Qrw1-qw-1-qw2Srw=1+qw21-qwQrw,where Srw is the mean daily demand of product w at HDCC r. Consequently, the holding cost of the working inventory is given by (5)∑i∈I1+qwhrw21-qwQrw.Using Eppen’s risk-pooling result [28], the amount of safety stock to ensure that stockouts occur with a probability of α or less is zαL∑i∈Iσiw2Yriw, and the safety stock cost at HDCCr is given by hrwzαL∑i∈Iσiw2Yriw.Therefore, the total annual inventory cost (CI) of product w at HDCCr is given by (6)CI=brwDrwQrw+crwDrwQrw+λerw∑i∈IμiwYriw-∑i∈IqwμiwZirw+∑i∈I1+qwhrw21-qwQrw+hrwzαL∑i∈Iσiw2Yriw.


### 3.3.4. Return Cost

Returns costs are incurred by customer returns, and they include inspection cost, repackaging cost, and shipping cost from customer zones to HDCCs.(1) Inspection Cost.Returned products will be inspected at HDCCs, and the total inspection cost is given by λ∑r∈R∑w∈W∑i∈IqwμiwgrwZriw.(2) Repackaging Cost. The returned products also need to be repackaged before they are resold to the market, and the total repackaging cost is given by λ∑r∈R∑w∈W∑i∈IqwμiwfrwZriw.(3) Shipping Cost from Customer Zones to HDCCs.The total shipping cost from customer zones to HDCCs is given by λ∑r∈R∑w∈W∑i∈IqwμiwdriwZriw.Therefore, the total annual return cost (CR) is given by(7)CR=λ∑r∈R∑w∈W∑i∈Iqwμiwgrw+frw+driwZriw.


### 3.3.5. Optimization Model

Given the individual cost functions above, the research problem under study can be formulated by a mixed-integer nonlinear program as follows:(8)min⁡Z=CL+CT+CI+CR=∑r∈RarXr+λ∑r∈R∑w∈W∑i∈IμiwdriwYriw+∑r∈R∑w∈WbrwDrwQrw+crwDrwQrw+λerw∑i∈IμiwYriw-∑i∈IqwμiwZirw+∑i∈I1+qwhrw21-qwQrw+hrwzαL∑i∈Iσiw2Yriw+λ∑r∈R∑w∈W∑i∈Iqwμiwgrw+frw+driwZriwsubject to(9)∑r∈RXr≥1;(10)∑r∈RYriw=1,∀i∈I,∀w∈W;(11)∑r∈RZirw=1,∀i∈I,∀w∈W;(12)∑r∈R∑w∈WYriw=W,∀i∈I;(13)∑r∈R∑w∈WZirw=W,∀i∈I;(14)Yriw≤Xr,∀i∈I,∀r∈R,∀w∈W;(15)Zirw≤Xr,∀i∈I,∀r∈R,∀w∈W;(16)∑w∈W∑i∈IYriw≥Xr,∀r∈R;(17)∑i∈IZirw≤M∑i∈IYriw,∀r∈R,∀w∈W;(18)Xr=0,1,∀r∈R;(19)Yriw=0,1,∀i∈I,∀r∈R,∀w∈W;(20)Zirw=0,1,∀i∈I,∀r∈R,∀w∈W;In this model, objective function (8) means that the total cost will be minimized. Constraint (9) means that at least one HDCC will be built. Constraint (10) means that the demand of a product in a customer zone will be fulfilled by only one HDCC. Constraint (11) means that the returns of a product from a customer zone will be collected by only one HDCC. Constraint (12) means that all types of products will be provided to a customer zone by HDCCs in the forward flow. Constraint (13) means that all types of returned products will be collected from a customer zone by HDCCs in the reverse flow. Constraints (14) and (15) mean that a customer zone can be assigned to a HDCC in the forward and reverse logistics flows, respectively, only if this HDCC is built. Constraint (16) means that a HDCC will serve at least one customer zone in the forward logistics network after it is built. Constraint (17) means that a HDCC will collect the returns of a product in the reverse flow only if it fulfills the demand of this product in the forward flow. Constraints (18)–(20) indicate thatXr, Yriw, and Zirw are binary variables.Note that decision variableQrw only appears in the objective function. Since the objective function is convex in Qrw, its optimal value can be obtained by taking the partial derivative of the objective function with respect to it, which gives(21)Qrw∗=∑i∈I21-qwbrw+crwDrw1+qwhrw.Therefore, the objective function can be rewritten as follows:(22)min⁡Z=λ∑r∈R∑w∈W∑i∈IμiwdriwYriw+∑r∈R∑w∈W∑i∈I2λhrw1+qwbrw+crwμiwYriw1-qw+λ∑r∈R∑w∈W∑i∈IerwμiwYriw-qwμiwZirw+∑r∈R∑w∈WhrwzαL∑i∈Iσiw2Yriw+∑r∈RarXr+λ∑r∈R∑w∈W∑i∈Iqwμiwgrw+frw+driwZriw


## 3.3.1. Location Cost

Obviously, the total location cost (CL) can be written as(1)CL=∑r∈RarXr.


## 3.3.2. Shipping Cost from HDCCs to Customer Zones

The total shipping cost (CT) from HDCCs to customer zones in the forward logistics flow can be expressed as follows:(2)CT=λ∑r∈R∑w∈W∑i∈IμiwdriwYriw.


## 3.3.3. Inventory Cost

Inventory strategy plays an important role in supply chain management. In this study, we assume that HDCCs order new products from the plant by using an (R,Q) policy with type I service, and the total inventory cost includes working inventory and safety stock costs. The (R,Q) model can be approximated by an economic order quantity (EOQ) model [9], and Axsater [27] shows that the maximum relative error introduced by using the EOQ model instead of the (Q,r) policy is 0.118. Therefore, we use an EOQ model to optimize the inventory strategy in this study, and the order frequency and quantity at a HDCC will be calculated by using the total demand of the customer zones served by this HDCC.In general, the working inventory cost includes(1) fixed cost of placing orders, (2) shipping cost from the plant to HDCCs, and (3) holding cost of the working inventory. The individual cost terms are formulated as follows.(1) Order Cost. The annual fixed cost of placing orders for product w from HDCCr to the plant is given by brwDrw/Qrw.(2) Shipping Cost from the Plant to HDCCs. Since false failure returns will be resold to the market after they are processed in HDCCs, shipping quantities from the plant to HDCCs will equal the difference between market demands and return volumes. Therefore, the total shipping cost from the plant to HDCCr per year is given by(3)crwDrwQrw+λerw∑i∈IμiwYriw-∑i∈IqwμiwZirw.(3) Holding Cost of the Working Inventory. In this study, the working inventory comes from two sources which are the new products ordered from the plant and the returned products collected from customer zones. Therefore, the forward and reverse logistics flows will be integrated at HDCCs, and the average inventory of product w at HDCCrduring periodT can be expressed as follows:(4)1T∫0TSrw-Srw-qwSrwTdt=Srw-12θrwT=Qrw1-qw-1-qw2Srw=1+qw21-qwQrw,where Srw is the mean daily demand of product w at HDCC r. Consequently, the holding cost of the working inventory is given by (5)∑i∈I1+qwhrw21-qwQrw.Using Eppen’s risk-pooling result [28], the amount of safety stock to ensure that stockouts occur with a probability of α or less is zαL∑i∈Iσiw2Yriw, and the safety stock cost at HDCCr is given by hrwzαL∑i∈Iσiw2Yriw.Therefore, the total annual inventory cost (CI) of product w at HDCCr is given by (6)CI=brwDrwQrw+crwDrwQrw+λerw∑i∈IμiwYriw-∑i∈IqwμiwZirw+∑i∈I1+qwhrw21-qwQrw+hrwzαL∑i∈Iσiw2Yriw.


## 3.3.4. Return Cost

Returns costs are incurred by customer returns, and they include inspection cost, repackaging cost, and shipping cost from customer zones to HDCCs.(1) Inspection Cost.Returned products will be inspected at HDCCs, and the total inspection cost is given by λ∑r∈R∑w∈W∑i∈IqwμiwgrwZriw.(2) Repackaging Cost. The returned products also need to be repackaged before they are resold to the market, and the total repackaging cost is given by λ∑r∈R∑w∈W∑i∈IqwμiwfrwZriw.(3) Shipping Cost from Customer Zones to HDCCs.The total shipping cost from customer zones to HDCCs is given by λ∑r∈R∑w∈W∑i∈IqwμiwdriwZriw.Therefore, the total annual return cost (CR) is given by(7)CR=λ∑r∈R∑w∈W∑i∈Iqwμiwgrw+frw+driwZriw.


## 3.3.5. Optimization Model

Given the individual cost functions above, the research problem under study can be formulated by a mixed-integer nonlinear program as follows:(8)min⁡Z=CL+CT+CI+CR=∑r∈RarXr+λ∑r∈R∑w∈W∑i∈IμiwdriwYriw+∑r∈R∑w∈WbrwDrwQrw+crwDrwQrw+λerw∑i∈IμiwYriw-∑i∈IqwμiwZirw+∑i∈I1+qwhrw21-qwQrw+hrwzαL∑i∈Iσiw2Yriw+λ∑r∈R∑w∈W∑i∈Iqwμiwgrw+frw+driwZriwsubject to(9)∑r∈RXr≥1;(10)∑r∈RYriw=1,∀i∈I,∀w∈W;(11)∑r∈RZirw=1,∀i∈I,∀w∈W;(12)∑r∈R∑w∈WYriw=W,∀i∈I;(13)∑r∈R∑w∈WZirw=W,∀i∈I;(14)Yriw≤Xr,∀i∈I,∀r∈R,∀w∈W;(15)Zirw≤Xr,∀i∈I,∀r∈R,∀w∈W;(16)∑w∈W∑i∈IYriw≥Xr,∀r∈R;(17)∑i∈IZirw≤M∑i∈IYriw,∀r∈R,∀w∈W;(18)Xr=0,1,∀r∈R;(19)Yriw=0,1,∀i∈I,∀r∈R,∀w∈W;(20)Zirw=0,1,∀i∈I,∀r∈R,∀w∈W;In this model, objective function (8) means that the total cost will be minimized. Constraint (9) means that at least one HDCC will be built. Constraint (10) means that the demand of a product in a customer zone will be fulfilled by only one HDCC. Constraint (11) means that the returns of a product from a customer zone will be collected by only one HDCC. Constraint (12) means that all types of products will be provided to a customer zone by HDCCs in the forward flow. Constraint (13) means that all types of returned products will be collected from a customer zone by HDCCs in the reverse flow. Constraints (14) and (15) mean that a customer zone can be assigned to a HDCC in the forward and reverse logistics flows, respectively, only if this HDCC is built. Constraint (16) means that a HDCC will serve at least one customer zone in the forward logistics network after it is built. Constraint (17) means that a HDCC will collect the returns of a product in the reverse flow only if it fulfills the demand of this product in the forward flow. Constraints (18)–(20) indicate thatXr, Yriw, and Zirw are binary variables.Note that decision variableQrw only appears in the objective function. Since the objective function is convex in Qrw, its optimal value can be obtained by taking the partial derivative of the objective function with respect to it, which gives(21)Qrw∗=∑i∈I21-qwbrw+crwDrw1+qwhrw.Therefore, the objective function can be rewritten as follows:(22)min⁡Z=λ∑r∈R∑w∈W∑i∈IμiwdriwYriw+∑r∈R∑w∈W∑i∈I2λhrw1+qwbrw+crwμiwYriw1-qw+λ∑r∈R∑w∈W∑i∈IerwμiwYriw-qwμiwZirw+∑r∈R∑w∈WhrwzαL∑i∈Iσiw2Yriw+∑r∈RarXr+λ∑r∈R∑w∈W∑i∈Iqwμiwgrw+frw+driwZriw


## 4. Solution Approach

Location problems are NP-hard in general [29], and LIPs can be more complicated. Evolutionary algorithms are widely used to solve LIPs, and differential evolution (DE), which was first introduced by Storn and Price [30], is popular method to solve nonlinear and nonconvex functions [31]. Although DE has a strong global search ability to solve such problems, its performance is not always guaranteed due to its weaknesses such as the poor local search ability. To obtain a more stable and effective approach, we design an improved differential evolution algorithm (IDE) to solve MCLIP. More specifically, IDE enhances the performance of DE by two ways: First, IDE introduces opposition-based learning (OBL) [32] to generate better initial populations. Second, it uses a novel mutation factorF based on the Gaussian distributionN(0,1) to enhance the diversity of populations.


### 4.1. Individuals and Operators

Generally, a DE algorithm has four main operations: initialization, mutation, crossover, and selection. After an initial population is created in the initialization step, the algorithm will execute mutation, crossover, and selection operations iteratively to improve the population until a stop criterion is satisfied. IDE follows the same procedure, and it uses notations shown in Table1.Table 1
IDE notations.
Notation Explanation Np population size N dimension (i.e., number of customer zones in this study) G maximum number of generations for evolution K number of consecutive iterations if no better individual Zi,g solution vector of individualsi in generation g (i=1, 2, …, Np) V Z i , g mutant vector of individualsi in generation g (i=1, 2, …, Np) U Z i , g trial vector of individualsi in generation g (i=1, 2, …, Np) Fi,g mutation factor of individualsi in generation g (i=1, 2, …, Np) CR crossover factor


#### 4.1.1. Initialization

When DE is applied to solve optimization problems, an encoding-decoding scheme will be needed to convert individuals in DE to the solutions to the optimization model, and vice versa. In IDE, an individual is represented by the matrix shown in (23).(23)Zi,g=x1,1i,g,x1,2i,g,…,x1,Ni,g,y1,N+1i,g,y1,N+2i,g,…,y1,2Ni,g…xW,1i,g,xW,2i,g,…,xW,Ni,g,yW,N+1i,g,yW,N+2i,g,…,yW,2Ni,g,i=1,2,…,Np;g=1,2,…,G;This matrix hasW rows and 2N columns. For the research problem under study, W andN denote the numbers of commodities and customers, respectively, and the first and lastN columns are related to the forward and reverse flows, respectively. More specifically, if xw,ji,g=r, then the demand of product w in customer zonej will be fulfilled by HDCCr in the forward flow, and if yw,ji,g=r, then the returns of product w from customer zonej will be collected by HDCCr in the reverse flow.Figure2 shows an example matrix for a MCLIP that has two commodities, ten customer zones, and five candidate HDCC locations. For product 1, HDCC 1 fulfills demands in customer zones 1, 6, and 10 and collects returns from customer zones 1, 7, and 9; HDCC 2 fulfills demands in customer zones 4 and 7 and collects returns from customer zones 2 and 6; HDCC 3 fulfills demands in customer zones 2, 3, 5, and 8 and collects returns from customer zones 4, 5, and 8; HDCC 5 fulfills the demand in customer zone 9 and collects returns from customer zones 3 and 10. For product 2, location decisions can be also explained in the same manner.Figure 2
An example individual.The aim of encoding is transforming the initial individual to practical solution of LIP. The formulation of decoding is as follows.To decode an individual to a solution to the optimization problem, we can use a mechanism shown in (24) and (25).(24)xw,ji,0=roundxL+rand·xU-xL;i=1,2,…,Np;j=1,2,…,N;w=1,2,…,W;(25)yw,ji,0=UroundyL+rand·lengthU-1;i=1,2,…,Np;j=N+1,N+2,…,2N;w=1,2,…,W;where U=unique(xw,1i,g,xw,2i,g,…,xw,Ni,g) are distinct HDCCs that provide productw in the forward flow, rand is a random variable that is uniformly distributed on [0,1],xL andxU are the lower and upper bounds of xw,ji,0, andround is the rounding function. Obviously, we havexL =yL = 1, andxU is decided by the number of candidate HDCC locations. The expression of yw,ji,0 ensures that product w can be returned to a HDCC in the reverse flow only if it provides this product in the forward flow, which is reinforced by Constraint (17) in the optimization model shown in Section 3.3.5.IDE introduce a new mechanism based on OBL to generate better initial populations. More specifically, IDE will improve an initial population that is generated randomly by checking its opposite population by the following steps.Step 1.
Generate an initial populationZi,0 randomly: (26)Zi,0=x1,1i,0,x1,2i,0,…,x1,Ni,0,y1,N+1i,0,y1,N+2i,0,…,y1,2Ni,0…xW,1i,0,xW,2i,0,…,xW,Ni,0,yW,N+1i,0,yW,N+2i,0,…,yW,2Ni,0,i=1,2,…,Np;Step 2.
Obtain its opposite populationZi,0¯:(27)Zi,0¯=x1,1i,0¯,x1,2i,0¯,…,x1,Ni,0¯,y1,N+1i,0¯,y1,N+2i,0¯,…,y1,2Ni,0¯…xW,1i,0¯,xW,2i,0¯,…,xW,Ni,0¯,yW,N+1i,0¯,yW,N+2i,0¯,…,yW,2Ni,0¯xw,ji,0¯=xL+xU-xw,ji,0;i=1,2,…,Np;j=1,2,…,N;w=1,2,…,W;where xw,ji,0¯ is the entry in row w and columnj of the ith individual in the opposite population and yw,ji,0¯ (i=1, 2, …,Np;j=N+1,N+2, …, 2N; w=1, 2, …, W) is calculated by (25).Step 3.
Select theNp best solutions from Zi,0∪Zi,0¯ as initial population.


#### 4.1.2. Mutation

The mutation strategy plays an important role DE to enhance the global optimization ability and convergence speed. Usually, DE will use DE/rand/1/bin as the mutation operation that will produce a new vector by adding the weighted difference of two randomly selected vectors to the third vector in generationg[30]. In this study, we design a new mutation strategy to improve the diversity of populations, and new mutant vectors vxw,ji,g+1 and vyw,ji,g+1 are generated by (28) and (29).(28)vxw,ji,g+1=roundxw,jr1,g+Fi,g·xw,jr2,g-xw,jr3,g;i=1,2,…,Np;j=1,2,…N;w=1,2,…,W;(29)vyw,ji,g+1=roundyw,jr1,g+Fi,g·yw,jr2,g-yw,jr3,g;i=1,2,…,Np;j=N+1,N+2,…,2N;w=1,2,…,W;wherei≠r1≠r2≠r3,r1,r2,r3 are randomly selected from {1,2,…,Np}, Fi,g =F·ξi,g and it amplifies the difference between two vectors,Fis an initial mutation factor, and ξi,g is independently generated from the normal distribution N(0,1). In the mutation step, the mutation factor will be updated first and then used to generate new mutant vectors, which will improve local and global search abilities and generate better mutant vectors. If vxw,ji,g+1and vyw,ji,g+1 in a new mutant vector exceed the upper and lower bounds, then the vector will be generated again within its range.


#### 4.1.3. Crossover

The crossover operation mixes mutate and target vectors to generate trial vectors and increase the diversity of populations. In IDE, the trial vector is generated by (30) and (31).(30)uxw,ji,g+1=vxw,ji,grand⁡i<CRxw,ji,gotherwisei=1,2,…,Np;j=1,2,…,N;w=1,2,…,W;(31)uyw,ji,g+1=vyw,ji,grand⁡i<CRyw,ji,gotherwisei=1,2,…,Np;j=N+1,N+2,…,2N;w=1,2,…,W;where rand(i) is random variable that is uniformly distributed on [0,1], and CR∈[0,1] is a predefined crossover rate.


#### 4.1.4. Feasibility Correction

A new individual that is generated by mutation or crossover operations may not be feasible by violating Constraint (17). Therefore, a feasibility correction procedure is developed to exclude infeasible individuals. In this study, a new individual will be checked after it is generated. If Constraint (17) is violated, then it will be replaced by a new individual that is randomly generated by (25).


#### 4.1.5. Selection

After a new population is generated, the selection operation will be executed to evaluate the objective values of all trial vectors. In IDE, a trial vectorf(UXi,g+1) will be compared with the corresponding target vectorf(Xi,g) on their objective values by using a greedy criterion[30]. More specifically, the selection operation can be expressed as follows:(32)Zi,g+1=UZi,g+1fUZi,g+1<fZi,gZi,gotherwise.


#### 4.1.6. Stop Criterion

In DE, the mutation, crossover, and selection operations will be executed iteratively until a stop criterion is satisfied. In this study, IDE will stop until any of the criteria below is satisfied.(1) No better solution is found in consecutiveK iterations.(2) The maximum number of iterations, i.e.,G, is reached.


### 4.2. Algorithm Flow

In summary, IDE consists of the following steps.Step 1.
Setg = 0 andk = 0.Step 2.
Create an initial population that consists ofNp individuals by applying OBL, and calculate the objective values of all individuals in the population.Step 3.
Execute mutation operation.Step 4.
Execute crossover operation.Step 5.
Perform feasibility correction.Step 6.
Execute selection operation.Step 7.
Setg =g + 1. Calculate the objective values of all individuals in the new population. If the minimal objective value is less than that in the old population, setk = 1. Otherwise, setk =k + 1.Step 8.
Ifk =K org =G, stop the algorithm. Otherwise, go to Step 4.


## 4.1. Individuals and Operators

Generally, a DE algorithm has four main operations: initialization, mutation, crossover, and selection. After an initial population is created in the initialization step, the algorithm will execute mutation, crossover, and selection operations iteratively to improve the population until a stop criterion is satisfied. IDE follows the same procedure, and it uses notations shown in Table1.Table 1
IDE notations.
Notation Explanation Np population size N dimension (i.e., number of customer zones in this study) G maximum number of generations for evolution K number of consecutive iterations if no better individual Zi,g solution vector of individualsi in generation g (i=1, 2, …, Np) V Z i , g mutant vector of individualsi in generation g (i=1, 2, …, Np) U Z i , g trial vector of individualsi in generation g (i=1, 2, …, Np) Fi,g mutation factor of individualsi in generation g (i=1, 2, …, Np) CR crossover factor


### 4.1.1. Initialization

When DE is applied to solve optimization problems, an encoding-decoding scheme will be needed to convert individuals in DE to the solutions to the optimization model, and vice versa. In IDE, an individual is represented by the matrix shown in (23).(23)Zi,g=x1,1i,g,x1,2i,g,…,x1,Ni,g,y1,N+1i,g,y1,N+2i,g,…,y1,2Ni,g…xW,1i,g,xW,2i,g,…,xW,Ni,g,yW,N+1i,g,yW,N+2i,g,…,yW,2Ni,g,i=1,2,…,Np;g=1,2,…,G;This matrix hasW rows and 2N columns. For the research problem under study, W andN denote the numbers of commodities and customers, respectively, and the first and lastN columns are related to the forward and reverse flows, respectively. More specifically, if xw,ji,g=r, then the demand of product w in customer zonej will be fulfilled by HDCCr in the forward flow, and if yw,ji,g=r, then the returns of product w from customer zonej will be collected by HDCCr in the reverse flow.Figure2 shows an example matrix for a MCLIP that has two commodities, ten customer zones, and five candidate HDCC locations. For product 1, HDCC 1 fulfills demands in customer zones 1, 6, and 10 and collects returns from customer zones 1, 7, and 9; HDCC 2 fulfills demands in customer zones 4 and 7 and collects returns from customer zones 2 and 6; HDCC 3 fulfills demands in customer zones 2, 3, 5, and 8 and collects returns from customer zones 4, 5, and 8; HDCC 5 fulfills the demand in customer zone 9 and collects returns from customer zones 3 and 10. For product 2, location decisions can be also explained in the same manner.Figure 2
An example individual.The aim of encoding is transforming the initial individual to practical solution of LIP. The formulation of decoding is as follows.To decode an individual to a solution to the optimization problem, we can use a mechanism shown in (24) and (25).(24)xw,ji,0=roundxL+rand·xU-xL;i=1,2,…,Np;j=1,2,…,N;w=1,2,…,W;(25)yw,ji,0=UroundyL+rand·lengthU-1;i=1,2,…,Np;j=N+1,N+2,…,2N;w=1,2,…,W;where U=unique(xw,1i,g,xw,2i,g,…,xw,Ni,g) are distinct HDCCs that provide productw in the forward flow, rand is a random variable that is uniformly distributed on [0,1],xL andxU are the lower and upper bounds of xw,ji,0, andround is the rounding function. Obviously, we havexL =yL = 1, andxU is decided by the number of candidate HDCC locations. The expression of yw,ji,0 ensures that product w can be returned to a HDCC in the reverse flow only if it provides this product in the forward flow, which is reinforced by Constraint (17) in the optimization model shown in Section 3.3.5.IDE introduce a new mechanism based on OBL to generate better initial populations. More specifically, IDE will improve an initial population that is generated randomly by checking its opposite population by the following steps.Step 1.
Generate an initial populationZi,0 randomly: (26)Zi,0=x1,1i,0,x1,2i,0,…,x1,Ni,0,y1,N+1i,0,y1,N+2i,0,…,y1,2Ni,0…xW,1i,0,xW,2i,0,…,xW,Ni,0,yW,N+1i,0,yW,N+2i,0,…,yW,2Ni,0,i=1,2,…,Np;Step 2.
Obtain its opposite populationZi,0¯:(27)Zi,0¯=x1,1i,0¯,x1,2i,0¯,…,x1,Ni,0¯,y1,N+1i,0¯,y1,N+2i,0¯,…,y1,2Ni,0¯…xW,1i,0¯,xW,2i,0¯,…,xW,Ni,0¯,yW,N+1i,0¯,yW,N+2i,0¯,…,yW,2Ni,0¯xw,ji,0¯=xL+xU-xw,ji,0;i=1,2,…,Np;j=1,2,…,N;w=1,2,…,W;where xw,ji,0¯ is the entry in row w and columnj of the ith individual in the opposite population and yw,ji,0¯ (i=1, 2, …,Np;j=N+1,N+2, …, 2N; w=1, 2, …, W) is calculated by (25).Step 3.
Select theNp best solutions from Zi,0∪Zi,0¯ as initial population.


### 4.1.2. Mutation

The mutation strategy plays an important role DE to enhance the global optimization ability and convergence speed. Usually, DE will use DE/rand/1/bin as the mutation operation that will produce a new vector by adding the weighted difference of two randomly selected vectors to the third vector in generationg[30]. In this study, we design a new mutation strategy to improve the diversity of populations, and new mutant vectors vxw,ji,g+1 and vyw,ji,g+1 are generated by (28) and (29).(28)vxw,ji,g+1=roundxw,jr1,g+Fi,g·xw,jr2,g-xw,jr3,g;i=1,2,…,Np;j=1,2,…N;w=1,2,…,W;(29)vyw,ji,g+1=roundyw,jr1,g+Fi,g·yw,jr2,g-yw,jr3,g;i=1,2,…,Np;j=N+1,N+2,…,2N;w=1,2,…,W;wherei≠r1≠r2≠r3,r1,r2,r3 are randomly selected from {1,2,…,Np}, Fi,g =F·ξi,g and it amplifies the difference between two vectors,Fis an initial mutation factor, and ξi,g is independently generated from the normal distribution N(0,1). In the mutation step, the mutation factor will be updated first and then used to generate new mutant vectors, which will improve local and global search abilities and generate better mutant vectors. If vxw,ji,g+1and vyw,ji,g+1 in a new mutant vector exceed the upper and lower bounds, then the vector will be generated again within its range.


### 4.1.3. Crossover

The crossover operation mixes mutate and target vectors to generate trial vectors and increase the diversity of populations. In IDE, the trial vector is generated by (30) and (31).(30)uxw,ji,g+1=vxw,ji,grand⁡i<CRxw,ji,gotherwisei=1,2,…,Np;j=1,2,…,N;w=1,2,…,W;(31)uyw,ji,g+1=vyw,ji,grand⁡i<CRyw,ji,gotherwisei=1,2,…,Np;j=N+1,N+2,…,2N;w=1,2,…,W;where rand(i) is random variable that is uniformly distributed on [0,1], and CR∈[0,1] is a predefined crossover rate.


### 4.1.4. Feasibility Correction

A new individual that is generated by mutation or crossover operations may not be feasible by violating Constraint (17). Therefore, a feasibility correction procedure is developed to exclude infeasible individuals. In this study, a new individual will be checked after it is generated. If Constraint (17) is violated, then it will be replaced by a new individual that is randomly generated by (25).


### 4.1.5. Selection

After a new population is generated, the selection operation will be executed to evaluate the objective values of all trial vectors. In IDE, a trial vectorf(UXi,g+1) will be compared with the corresponding target vectorf(Xi,g) on their objective values by using a greedy criterion[30]. More specifically, the selection operation can be expressed as follows:(32)Zi,g+1=UZi,g+1fUZi,g+1<fZi,gZi,gotherwise.


### 4.1.6. Stop Criterion

In DE, the mutation, crossover, and selection operations will be executed iteratively until a stop criterion is satisfied. In this study, IDE will stop until any of the criteria below is satisfied.(1) No better solution is found in consecutiveK iterations.(2) The maximum number of iterations, i.e.,G, is reached.


## 4.1.1. Initialization

When DE is applied to solve optimization problems, an encoding-decoding scheme will be needed to convert individuals in DE to the solutions to the optimization model, and vice versa. In IDE, an individual is represented by the matrix shown in (23).(23)Zi,g=x1,1i,g,x1,2i,g,…,x1,Ni,g,y1,N+1i,g,y1,N+2i,g,…,y1,2Ni,g…xW,1i,g,xW,2i,g,…,xW,Ni,g,yW,N+1i,g,yW,N+2i,g,…,yW,2Ni,g,i=1,2,…,Np;g=1,2,…,G;This matrix hasW rows and 2N columns. For the research problem under study, W andN denote the numbers of commodities and customers, respectively, and the first and lastN columns are related to the forward and reverse flows, respectively. More specifically, if xw,ji,g=r, then the demand of product w in customer zonej will be fulfilled by HDCCr in the forward flow, and if yw,ji,g=r, then the returns of product w from customer zonej will be collected by HDCCr in the reverse flow.Figure2 shows an example matrix for a MCLIP that has two commodities, ten customer zones, and five candidate HDCC locations. For product 1, HDCC 1 fulfills demands in customer zones 1, 6, and 10 and collects returns from customer zones 1, 7, and 9; HDCC 2 fulfills demands in customer zones 4 and 7 and collects returns from customer zones 2 and 6; HDCC 3 fulfills demands in customer zones 2, 3, 5, and 8 and collects returns from customer zones 4, 5, and 8; HDCC 5 fulfills the demand in customer zone 9 and collects returns from customer zones 3 and 10. For product 2, location decisions can be also explained in the same manner.Figure 2
An example individual.The aim of encoding is transforming the initial individual to practical solution of LIP. The formulation of decoding is as follows.To decode an individual to a solution to the optimization problem, we can use a mechanism shown in (24) and (25).(24)xw,ji,0=roundxL+rand·xU-xL;i=1,2,…,Np;j=1,2,…,N;w=1,2,…,W;(25)yw,ji,0=UroundyL+rand·lengthU-1;i=1,2,…,Np;j=N+1,N+2,…,2N;w=1,2,…,W;where U=unique(xw,1i,g,xw,2i,g,…,xw,Ni,g) are distinct HDCCs that provide productw in the forward flow, rand is a random variable that is uniformly distributed on [0,1],xL andxU are the lower and upper bounds of xw,ji,0, andround is the rounding function. Obviously, we havexL =yL = 1, andxU is decided by the number of candidate HDCC locations. The expression of yw,ji,0 ensures that product w can be returned to a HDCC in the reverse flow only if it provides this product in the forward flow, which is reinforced by Constraint (17) in the optimization model shown in Section 3.3.5.IDE introduce a new mechanism based on OBL to generate better initial populations. More specifically, IDE will improve an initial population that is generated randomly by checking its opposite population by the following steps.Step 1.
Generate an initial populationZi,0 randomly: (26)Zi,0=x1,1i,0,x1,2i,0,…,x1,Ni,0,y1,N+1i,0,y1,N+2i,0,…,y1,2Ni,0…xW,1i,0,xW,2i,0,…,xW,Ni,0,yW,N+1i,0,yW,N+2i,0,…,yW,2Ni,0,i=1,2,…,Np;Step 2.
Obtain its opposite populationZi,0¯:(27)Zi,0¯=x1,1i,0¯,x1,2i,0¯,…,x1,Ni,0¯,y1,N+1i,0¯,y1,N+2i,0¯,…,y1,2Ni,0¯…xW,1i,0¯,xW,2i,0¯,…,xW,Ni,0¯,yW,N+1i,0¯,yW,N+2i,0¯,…,yW,2Ni,0¯xw,ji,0¯=xL+xU-xw,ji,0;i=1,2,…,Np;j=1,2,…,N;w=1,2,…,W;where xw,ji,0¯ is the entry in row w and columnj of the ith individual in the opposite population and yw,ji,0¯ (i=1, 2, …,Np;j=N+1,N+2, …, 2N; w=1, 2, …, W) is calculated by (25).Step 3.
Select theNp best solutions from Zi,0∪Zi,0¯ as initial population.


## 4.1.2. Mutation

The mutation strategy plays an important role DE to enhance the global optimization ability and convergence speed. Usually, DE will use DE/rand/1/bin as the mutation operation that will produce a new vector by adding the weighted difference of two randomly selected vectors to the third vector in generationg[30]. In this study, we design a new mutation strategy to improve the diversity of populations, and new mutant vectors vxw,ji,g+1 and vyw,ji,g+1 are generated by (28) and (29).(28)vxw,ji,g+1=roundxw,jr1,g+Fi,g·xw,jr2,g-xw,jr3,g;i=1,2,…,Np;j=1,2,…N;w=1,2,…,W;(29)vyw,ji,g+1=roundyw,jr1,g+Fi,g·yw,jr2,g-yw,jr3,g;i=1,2,…,Np;j=N+1,N+2,…,2N;w=1,2,…,W;wherei≠r1≠r2≠r3,r1,r2,r3 are randomly selected from {1,2,…,Np}, Fi,g =F·ξi,g and it amplifies the difference between two vectors,Fis an initial mutation factor, and ξi,g is independently generated from the normal distribution N(0,1). In the mutation step, the mutation factor will be updated first and then used to generate new mutant vectors, which will improve local and global search abilities and generate better mutant vectors. If vxw,ji,g+1and vyw,ji,g+1 in a new mutant vector exceed the upper and lower bounds, then the vector will be generated again within its range.


## 4.1.3. Crossover

The crossover operation mixes mutate and target vectors to generate trial vectors and increase the diversity of populations. In IDE, the trial vector is generated by (30) and (31).(30)uxw,ji,g+1=vxw,ji,grand⁡i<CRxw,ji,gotherwisei=1,2,…,Np;j=1,2,…,N;w=1,2,…,W;(31)uyw,ji,g+1=vyw,ji,grand⁡i<CRyw,ji,gotherwisei=1,2,…,Np;j=N+1,N+2,…,2N;w=1,2,…,W;where rand(i) is random variable that is uniformly distributed on [0,1], and CR∈[0,1] is a predefined crossover rate.


## 4.1.4. Feasibility Correction

A new individual that is generated by mutation or crossover operations may not be feasible by violating Constraint (17). Therefore, a feasibility correction procedure is developed to exclude infeasible individuals. In this study, a new individual will be checked after it is generated. If Constraint (17) is violated, then it will be replaced by a new individual that is randomly generated by (25).


## 4.1.5. Selection

After a new population is generated, the selection operation will be executed to evaluate the objective values of all trial vectors. In IDE, a trial vectorf(UXi,g+1) will be compared with the corresponding target vectorf(Xi,g) on their objective values by using a greedy criterion[30]. More specifically, the selection operation can be expressed as follows:(32)Zi,g+1=UZi,g+1fUZi,g+1<fZi,gZi,gotherwise.


## 4.1.6. Stop Criterion

In DE, the mutation, crossover, and selection operations will be executed iteratively until a stop criterion is satisfied. In this study, IDE will stop until any of the criteria below is satisfied.(1) No better solution is found in consecutiveK iterations.(2) The maximum number of iterations, i.e.,G, is reached.


## 4.2. Algorithm Flow

In summary, IDE consists of the following steps.Step 1.
Setg = 0 andk = 0.Step 2.
Create an initial population that consists ofNp individuals by applying OBL, and calculate the objective values of all individuals in the population.Step 3.
Execute mutation operation.Step 4.
Execute crossover operation.Step 5.
Perform feasibility correction.Step 6.
Execute selection operation.Step 7.
Setg =g + 1. Calculate the objective values of all individuals in the new population. If the minimal objective value is less than that in the old population, setk = 1. Otherwise, setk =k + 1.Step 8.
Ifk =K org =G, stop the algorithm. Otherwise, go to Step 4.


## 5. Computational Results

In this section, numerical analysis is presented to validate IDE and provide managerial insights. First, we show how to obtain the optimal setting in IDE for the best possible performance. Second, the solution to an example MCLIP is presented. Third, sensitivity analysis is conducted to study the influences of the parameters in MCLIP. Last, IDE is compared with Lingo 11 and DE to validate its performance. All the experiments are implemented by Java JDK 1.7 on a Windows PC (AMD A10-9600P RADEON R5, 10 COMPUTE CORES 4C+6G 2.40GHz; RAM: 4.00 GB DDR; OS: Windows 10).To get started, we use an example MCLIP that consists of eight candidate HDCC locations, eighty customer zones, and three commodities. The candidate HDCC locations and customer zones are uniformly distributed on a grid of[0,50]×[0,50], and the parameters are shown in Table 2.Table 2
MCLIP parameters.
Parameters Value Parameters Value Parameters Value a r U [ 1000,1500 ] b r w U [ 5,10 ] c r w 10 e r w 5 f r w U [ 1,3 ] g r w U [ 1,3 ] h r w U [ 1,3 ] q w U [ 0.3,0.6 ] µ i w U [ 20,30 ] σ i w ​ 2 σ i ​ 2=µi z α 1.96 L 1 Remark:U[a,b] denotes the uniform distribution over [a,b].


### 5.1. Parameter Analysis

The performance of DE can be significantly affected by its parameters or the combinations of those parameters. In this subsection, we test the performance of IDE under differentFandCR values. In this experiment, we setNp= 4S(S=N∗W) and the algorithm was executed 30 times under each setting, and the numerical results are shown in Table 3.Table 3
IDE parameter test results.
CR F 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 OS Ratio 0.1 96.67% 56.67% 36.67% 3.33% 0.00% 0.00% 0.00% 0.00% 0.00% CPU Time 28.42 22.74 20.13 18.90 17.79 17.16 17.15 16.78 17.33  OS Ratio 0.2 96.67% 83.33% 23.33% 10.00% 3.33% 0.00% 0.00% 0.00% 0.00% CPU Time 42.01 35.09 31.17 29.26 27.49 26.76 26.41 26.22 27.18  OS Ratio 0.3 0.00% 56.67% 50.00% 20.00% 3.33% 0.00% 0.00% 0.00% 0.00% CPU Time 47.40 46.46 43.12 40.61 38.78 37.51 36.65 36.72 37.19  OS Ratio 0.4 0.00% 0.00% 0.00% 3.33% 6.67% 0.00% 0.00% 0.00% 0.00% CPU Time 46.89 47.66 46.78 46.04 45.43 44.88 44.63 44.37 44.32  OS Ratio 0.5 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% CPU Time 47.01 47.00 45.49 46.61 46.18 46.28 46.08 44.13 45.53  OS Ratio 0.6 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% CPU Time 45.89 47.56 47.17 46.26 44.70 46.78 47.01 45.96 44.89  OS Ratio 0.7 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% CPU Time 47.45 47.03 46.94 47.53 44.71 45.39 46.70 45.14 45.11  OS Ratio 0.8 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% CPU Time 47.23 47.82 46.28 45.75 45.43 45.27 43.85 46.25 46.59  OS Ratio 0.9 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% CPU Time 48.08 44.82 44.95 46.98 46.45 47.66 45.06 45.77 43.33 Remark: in Table3, the optimal solutions are obtained by Lingo 11, and “OS Ratio” means the percentage of finding optimal solutions by IDE.From Table3, we can see thatF andCR have a significant impact on the searching ability and convergence speed of IDE, which are indicated by OS ratio and CPU time, respectively. Obviously, IDE has the best performance whenF = 0.1 andCR = 0.1, and this setting will be used in the subsequent experiments.


### 5.2. Illustrative Example

This subsection shows the optimal to the example MCLIP. Table4 shows the optimal solution which indicates HDCC locations and the corresponding customer zones in the forward and reverse flows, and Table 5 shows the optimal value and individual costs in the optimal solution.Table 4
Optimal solution to the example MCLIP.
HDCC Location Commodity Customer Zone Customer Zone (Forward Flow) (Reverse Flow) 1 1 48,64,75 48,64,75 2 48,64,75 48,64,75 3 48,64,75 48,64,75  2 1 4,9.10,11,26,29,30,39,44,45,47,57,68,72 4,9.10,11,26,29,30,39,44,45,47,57,68,72 2 4,9.10,11,26,29,30,39,44,45,47,57,68,72 4,9.10,11,26,29,30,39,44,45,47,57,68,72 3 4,9.10,11,26,29,30,39,44,45,47,57,68,72 4,9.10,11,12,26,29,30,39,44,45,47,57,68,71,72  3 1 12,14,18,22,32,34,43,60,70,71 12,14,18,22,32,34,43,60,70,71 2 12,14,18,22,32,34,43,60,70,71 12,14,18,22,32,34,38,43,60,70,71,79 3 12,14,18,22,32,34,43,60,70,71 14,18,22,32,34,43,60,70  4 1 2,7,13,35,52,56,58,63,76,80 2,7,35,52,56,58,80 2 2,7,13,35,52,56,58,63,76,80 2,7,13,35,52,56,58,63,76,80 3 2,7,13,35,52,56,58,63,76,80 2,7,8,13,21,35,51,52,56,58,61,63,76,80  5 1 5,15,16,19,31,37,38,42,49,62,77,79 5,15,16,19,31,37,38,42,49,62,77,79 2 5,15,16,19,31,37,38,42,49,62,77,79 5,15,16,19,31,37,42,49,62,77, 3 5,15,16,19,31,37,38,42,49,62,77,79 5,15,16,19,31,37,38,42,49,62,77,79  6 1 3,6,17,20,23,24,36,67,69 3,17,20,23,24,36,67,69 2 3,6,17,20,23,24,36,67,69 3,6,17,20,23,24,36,67,69 3 3,6,17,20,23,24,36,67,69 3,6,17,20,23,24,36,67,69  7 1 8,21,25,28,46,51,53,54,61,65,73,74 6,8,13,21,25,27,2846,51,53,54,61,63,65,73,74,76 2 8,21,25,28,46,51,53,54,61,65,73,74 8,21,25,28,46,51,53,54,61,65,73,74 3 8,21,25,28,46,51,53,54,61,65,73,74 25,28,46,53,54,65,73,74  8 1 1,27,33,40,41,50,55,59,66,78 1,33,40,41,50,55,59,66,78 2 1,27,33,40,41,50,55,59,66,78 1,27,33,40,41,50,55,59,66,78 3 1,27,33,40,41,50,55,59,66,78 1,27,33,40,50,55,59,66,78Table 5
Optimal value and individual costs.
C L C T C I C R Total 10262 16651085.12 4665316.68 11617106.42 32943770.22


### 5.3. Sensitivity Analysis

The optimal solution of a MCLIP can be affected significantly by important business parameters. In this subsection, sensitivity analysis is conducted onqw, driw, brw, crw, erw, hrw, grw, and frw to study their influences. More specifically, a parameter will be tested by changing its value by (-30%, 30%) and fixing the other parameters and the numerical results are shown in Table 6.Table 6
Sensitivity analysis results.
Δ -30.00% -20.00% -10.00% 0.00% 10.00% 20.00% 30.00% q w C I 5967659.67 5532973.87 5098809.85 4665316.678 4232715.54 3801356.65 3371846.74 Gap 27.92% 18.60% 9.29% 0.00% -9.27% -18.52% -27.73% C R 8131974.50 9293685.14 10455395.78 11617106.42 12778817.06 13940527.71 15102238.35 Gap -30.00% -20.00% -10.00% 0.00% 10.00% 20.00% 30.00% OV 30760981.29 31488006.13 32215552.76 32943770.22 33672879.73 34403231.48 35135432.21 Gap -6.63% -4.42% -2.21% 0.00% 2.21% 4.43% 6.65%  d r i w C T 11655759.59 13320868.10 14985976.61 16651085.12 18316193.63 19981302.15 21646410.66 Gap -30.00% -20.00% -10.00% 0.00% 10.00% 20.00% 30.00% C R 9141586.34 9969691.59 10793951.68 11617106.42 12437398.32 13255554.35 14073565.94 Gap -21.31% -14.18% -7.09% 0.00% 7.06% 14.10% 21.15% OV 25472924.60 27966138.37 30455506.97 32943770.22 35429170.64 37912435.18 40395555.28 Gap -22.68% -15.11% -7.55% 0.00% 7.54% 15.08% 22.62%  b r w C I 4659545.89 4661513.59 4663436.18 4665316.678 4667157.74 4668961.77 4670730.94 Gap -0.12% -0.08% -0.04% 0.00% 0.04% 0.08% 0.12% OV 32937999.43 32939967.13 32941889.73 32943770.22 32945611.28 32947415.32 32949184.48 Gap -0.02% -0.01% -0.01% 0.00% 0.01% 0.01% 0.02%  c r w C I 4657224.85 4660009.99 4662704.48 4665316.678 4667853.73 4670321.79 4672726.23 Gap -0.17% -0.11% -0.06% 0.00% 0.05% 0.11% 0.16% OV 32935678.39 32938463.54 32941158.03 32943770.22 32946307.27 32948775.33 32951179.77 Gap -0.02% -0.02% -0.01% 0.00% 0.01% 0.02% 0.02%  e r w C I 3292753.68 3750274.68 4207795.68 4665316.678 5122837.68 5580358.68 6037879.68 Gap -29.42% -19.61% -9.81% 0.00% 9.81% 19.61% 29.42% OV 31571207.22 32028728.22 32486249.22 32943770.22 33401291.22 33858812.22 34316333.22 Gap -4.17% -2.78% -1.39% 0.00% 1.39% 2.78% 4.17%  h r w C I 4650403.35 4655668.91 4660623.12 4665316.678 4669787.84 4674066.35 4678175.84 Gap -0.32% -0.21% -0.10% 0.00% 0.10% 0.19% 0.28% OV 32928856.89 32934122.46 32939076.67 32943770.22 32948241.39 32952519.89 32956629.38 Gap -0.05% -0.03% -0.01% 0.00% 0.01% 0.03% 0.04%  g r w C R 11058612.56 11245420.25 11431721.41 11617106.42 11802239.42 11987372.42 12172505.42 Gap -4.81% -3.20% -1.60% 0.00% 1.59% 3.19% 4.78% OV 32385276.36 32572084.05 32758385.21 32943770.22 33128903.22 33314036.22 33499169.22 Gap -1.70% -1.13% -0.56% 0.00% 0.56% 1.12% 1.69%  f r w C R 11146187.90 11303872.70 11461479.65 11617106.42 11770266.78 11922189.48 12073560.78 Gap -4.05% -2.70% -1.34% 0.00% 1.32% 2.63% 3.93% OV 32472851.69 32630536.49 32788143.45 32943770.22 33096930.58 33248853.28 33400224.58 Gap -1.43% -0.95% -0.47% 0.00% 0.46% 0.93% 1.39% Remark: in Table6, “OV” means the optimal value, and “Gap” means the change of the optimal value. Note that the gap will be zero if no parameter is changed.From Table6, we can see that optimal values can change significantly due to the change of those parameters. For example, the optimal value will decrease by 6.63% if qw decreases by 30%, and the optimal value will increase by 6.65% if qw increases by 30%. The influences of the parameters on individual costs are summarized as follows: When qw increases,CR will increase,CI and will decrease, butCL andCT will not change. When driw increases,CT andCR will increase butCL andCI will not change. When brw, crw, erw, hrw increase, onlyCI will increase and the other costs will not change. When grw or frw increases, onlyCI will increase. We can see that, in general, driw and qw have the most significant influences. Moreover, the optimal value is more sensitive when driw, qw, and erw are relatively small. In practice, business managers may focus more on those parameters to stabilize the supply chain performance.


### 5.4. Performance Analysis

In this subsection, IDE is compared with Lingo 11 and DE to validate its performance in terms of solution accuracy and CPU time efficiency, and the three approaches are tested on a set of small- (i.e., 20×5×2/3/4, 40×6×2/3/4), medium- (i.e., 50×5/6/7×3, 70×6/7/8×3), and large-sized problems (i.e., 80×8×2/3/4, 100×8/9/10×3). Since population sizeNp is an important parameter in IDE and it may vary according to the problem sizes, its value needs to be set properly to get the best possible performance. According to Figures 3–5 (in which “red line” represents solution accuracy in terms of the percentage of finding the optimal solutions, and “blue line” represents the mean CPU time), IDE will have the best performance whenNp= 3S, 4S, 5S for small-, medium-, and large-sized problems, respectively. Therefore, we will use these settings in the subsequent experiments.Figure 3
Analysis ofNpfor small-sized problems (40×6×4).Figure 4
Analysis ofNpfor medium-sized problems (70×8×3).Figure 5
Analysis ofNpfor large-sized problems (100×10×3).Table7 shows the comparison results between IDE, DE, and Lingo 11. To avoid biased results, IDE and DE were executed 30 times for each test instance. Therefore, in this table, “O.V.” means the optimal value, “O.V. (Best)” and “O.V. (Mean)” represent the best and mean optimal values in 30 runs, respectively, “OS Ratio” means the percentage of finding the optimal solutions, and “S.D.” means the standard deviation of optimal values in 30 runs.Table 7
Lingo 11 versus DE versus IDE.
Instance Lingo 11 DE IDE (T1-T2) / T1 I R | W | O.V. CPU Time O.V. O.V. OS Ratio CPU Time S.D. O.V. O.V. OS Ratio CPU Time S.D. (S) (Best) (Mean) (S) - T1 (Best) (Mean) (S) - T2 20 5 2 6987580 6 6987580.30 6997722.08 21/30 0.14 0.01 6987580.30 6987580.30 30/30 0.14 0.01 -3.95% 3 10355510 15 10355509.01 10358913.98 26/30 0.34 0.05 10355509.01 10355509.01 30/30 0.35 0.03 -2.63% 4 13487440 27 13487441.55 13488007.01 29/30 0.65 0.05 13487441.55 13487441.55 30/30 0.67 0.05 -3.08%  40 6 2 12598710 53 12598709.04 12608013.39 16/30 1.01 0.06 12598709.04 12598750.74 29/30 0.89 0.04 11.68% 3 18642000 135 18642000.91 18643688.34 26/30 2.71 0.08 18642000.91 18642000.91 30/30 2.50 0.08 7.71% 4 24367940 263 24367942.66 24368298.67 26/30 5.71 0.14 24367942.66 24367942.66 30/30 5.19 0.10 9.02%  50 5 3 23313920 123 23313920.93 23316057.84 26/30 4.85 0.32 23313920.93 23313920.93 30/30 5.03 0.21 -3.68% 6 23017930 229 23017934.37 23018947.56 27/30 6.78 0.16 23017934.37 23017934.37 30/30 6.05 0.12 10.76% 7 19471180 326 19471179.29 19471323.45 28/30 9.57 0.17 19471179.29 19471179.29 30/30 7.45 0.16 22.14%  70 6 3 31652380 472 31652384.45 31652401.08 29/30 16.17 0.50 31652384.45 31652384.45 30/30 16.26 0.60 -0.56% 7 30150360 1394 30150355.84 30150535.41 27/30 21.51 0.47 30150355.84 30150355.84 30/30 18.28 0.42 14.99% 8 28295930 1993 28295927.33 28296407.75 25/30 22.89 0.43 28295927.33 28295939.60 29/30 20.85 0.53 8.92%  80 8 2 21939770 721 21939767.66 21939991.65 28/30 13.28 0.24 21939767.66 21939767.66 30/30 12.05 0.31 9.31% 3 32943770 2974 32943770.22 32943973.24 29/30 41.51 1.15 32943770.22 32943770.22 30/30 35.76 0.72 13.83% 4 43234620 >3600 43234618.55 43234779.15 29/30 91.60 1.84 43234618.55 43234618.55 30/30 74.05 1.35 19.15%  100 8 3 39081930 >3600 39081579.93 39081760.66 28/30 115.95 2.07 39081579.93 39081579.93 30/30 69.93 1.45 39.69% 9 38614670 >3600 39767341.71 40015199.70 0/30 102.18 0.90 38614396.07 38614592.75 29/30 83.84 1.72 17.95% 10 36450520 >3600 37621661.96 37827792.61 0/30 124.76 0.74 36450173.57 36450173.57 30/30 88.87 2.74 28.76%From Table7, we can see that IDE is the most effective and efficient method to solve MCLIPs, and its advantages over Lingo 11 and DE are summarized as follows:(1)
The optimal values obtained by IDE and Lingo 11 were almost identical, but IDE has better runtimes than Lingo 11. This indicates that IDE is much more efficient than Lingo 11 with the same solution accuracy.(2)
IDE and DE can achieve similar solution accuracies on small- and medium-sized problems. However, the optimal values from DE are significantly larger than those from IDE on large-sized problems. Particularly, DE cannot find optimal solutions for the problems whose sizes are over (100×9×3). This indicates that IDE has a great global search ability and is more stable and consistent than DE.(3)
IDE has better runtimes than DE in most test instances, and it has a greater advantage on large-sized problems. This indicates that IDE is more efficient than DE.According to the comparison results, we can see that IDE is the best approach to solve MCLIPs compared with DE and Lingo 11.


## 5.1. Parameter Analysis

The performance of DE can be significantly affected by its parameters or the combinations of those parameters. In this subsection, we test the performance of IDE under differentFandCR values. In this experiment, we setNp= 4S(S=N∗W) and the algorithm was executed 30 times under each setting, and the numerical results are shown in Table 3.Table 3
IDE parameter test results.
CR F 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 OS Ratio 0.1 96.67% 56.67% 36.67% 3.33% 0.00% 0.00% 0.00% 0.00% 0.00% CPU Time 28.42 22.74 20.13 18.90 17.79 17.16 17.15 16.78 17.33  OS Ratio 0.2 96.67% 83.33% 23.33% 10.00% 3.33% 0.00% 0.00% 0.00% 0.00% CPU Time 42.01 35.09 31.17 29.26 27.49 26.76 26.41 26.22 27.18  OS Ratio 0.3 0.00% 56.67% 50.00% 20.00% 3.33% 0.00% 0.00% 0.00% 0.00% CPU Time 47.40 46.46 43.12 40.61 38.78 37.51 36.65 36.72 37.19  OS Ratio 0.4 0.00% 0.00% 0.00% 3.33% 6.67% 0.00% 0.00% 0.00% 0.00% CPU Time 46.89 47.66 46.78 46.04 45.43 44.88 44.63 44.37 44.32  OS Ratio 0.5 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% CPU Time 47.01 47.00 45.49 46.61 46.18 46.28 46.08 44.13 45.53  OS Ratio 0.6 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% CPU Time 45.89 47.56 47.17 46.26 44.70 46.78 47.01 45.96 44.89  OS Ratio 0.7 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% CPU Time 47.45 47.03 46.94 47.53 44.71 45.39 46.70 45.14 45.11  OS Ratio 0.8 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% CPU Time 47.23 47.82 46.28 45.75 45.43 45.27 43.85 46.25 46.59  OS Ratio 0.9 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% CPU Time 48.08 44.82 44.95 46.98 46.45 47.66 45.06 45.77 43.33 Remark: in Table3, the optimal solutions are obtained by Lingo 11, and “OS Ratio” means the percentage of finding optimal solutions by IDE.From Table3, we can see thatF andCR have a significant impact on the searching ability and convergence speed of IDE, which are indicated by OS ratio and CPU time, respectively. Obviously, IDE has the best performance whenF = 0.1 andCR = 0.1, and this setting will be used in the subsequent experiments.


## 5.2. Illustrative Example

This subsection shows the optimal to the example MCLIP. Table4 shows the optimal solution which indicates HDCC locations and the corresponding customer zones in the forward and reverse flows, and Table 5 shows the optimal value and individual costs in the optimal solution.Table 4
Optimal solution to the example MCLIP.
HDCC Location Commodity Customer Zone Customer Zone (Forward Flow) (Reverse Flow) 1 1 48,64,75 48,64,75 2 48,64,75 48,64,75 3 48,64,75 48,64,75  2 1 4,9.10,11,26,29,30,39,44,45,47,57,68,72 4,9.10,11,26,29,30,39,44,45,47,57,68,72 2 4,9.10,11,26,29,30,39,44,45,47,57,68,72 4,9.10,11,26,29,30,39,44,45,47,57,68,72 3 4,9.10,11,26,29,30,39,44,45,47,57,68,72 4,9.10,11,12,26,29,30,39,44,45,47,57,68,71,72  3 1 12,14,18,22,32,34,43,60,70,71 12,14,18,22,32,34,43,60,70,71 2 12,14,18,22,32,34,43,60,70,71 12,14,18,22,32,34,38,43,60,70,71,79 3 12,14,18,22,32,34,43,60,70,71 14,18,22,32,34,43,60,70  4 1 2,7,13,35,52,56,58,63,76,80 2,7,35,52,56,58,80 2 2,7,13,35,52,56,58,63,76,80 2,7,13,35,52,56,58,63,76,80 3 2,7,13,35,52,56,58,63,76,80 2,7,8,13,21,35,51,52,56,58,61,63,76,80  5 1 5,15,16,19,31,37,38,42,49,62,77,79 5,15,16,19,31,37,38,42,49,62,77,79 2 5,15,16,19,31,37,38,42,49,62,77,79 5,15,16,19,31,37,42,49,62,77, 3 5,15,16,19,31,37,38,42,49,62,77,79 5,15,16,19,31,37,38,42,49,62,77,79  6 1 3,6,17,20,23,24,36,67,69 3,17,20,23,24,36,67,69 2 3,6,17,20,23,24,36,67,69 3,6,17,20,23,24,36,67,69 3 3,6,17,20,23,24,36,67,69 3,6,17,20,23,24,36,67,69  7 1 8,21,25,28,46,51,53,54,61,65,73,74 6,8,13,21,25,27,2846,51,53,54,61,63,65,73,74,76 2 8,21,25,28,46,51,53,54,61,65,73,74 8,21,25,28,46,51,53,54,61,65,73,74 3 8,21,25,28,46,51,53,54,61,65,73,74 25,28,46,53,54,65,73,74  8 1 1,27,33,40,41,50,55,59,66,78 1,33,40,41,50,55,59,66,78 2 1,27,33,40,41,50,55,59,66,78 1,27,33,40,41,50,55,59,66,78 3 1,27,33,40,41,50,55,59,66,78 1,27,33,40,50,55,59,66,78Table 5
Optimal value and individual costs.
C L C T C I C R Total 10262 16651085.12 4665316.68 11617106.42 32943770.22


## 5.3. Sensitivity Analysis

The optimal solution of a MCLIP can be affected significantly by important business parameters. In this subsection, sensitivity analysis is conducted onqw, driw, brw, crw, erw, hrw, grw, and frw to study their influences. More specifically, a parameter will be tested by changing its value by (-30%, 30%) and fixing the other parameters and the numerical results are shown in Table 6.Table 6
Sensitivity analysis results.
Δ -30.00% -20.00% -10.00% 0.00% 10.00% 20.00% 30.00% q w C I 5967659.67 5532973.87 5098809.85 4665316.678 4232715.54 3801356.65 3371846.74 Gap 27.92% 18.60% 9.29% 0.00% -9.27% -18.52% -27.73% C R 8131974.50 9293685.14 10455395.78 11617106.42 12778817.06 13940527.71 15102238.35 Gap -30.00% -20.00% -10.00% 0.00% 10.00% 20.00% 30.00% OV 30760981.29 31488006.13 32215552.76 32943770.22 33672879.73 34403231.48 35135432.21 Gap -6.63% -4.42% -2.21% 0.00% 2.21% 4.43% 6.65%  d r i w C T 11655759.59 13320868.10 14985976.61 16651085.12 18316193.63 19981302.15 21646410.66 Gap -30.00% -20.00% -10.00% 0.00% 10.00% 20.00% 30.00% C R 9141586.34 9969691.59 10793951.68 11617106.42 12437398.32 13255554.35 14073565.94 Gap -21.31% -14.18% -7.09% 0.00% 7.06% 14.10% 21.15% OV 25472924.60 27966138.37 30455506.97 32943770.22 35429170.64 37912435.18 40395555.28 Gap -22.68% -15.11% -7.55% 0.00% 7.54% 15.08% 22.62%  b r w C I 4659545.89 4661513.59 4663436.18 4665316.678 4667157.74 4668961.77 4670730.94 Gap -0.12% -0.08% -0.04% 0.00% 0.04% 0.08% 0.12% OV 32937999.43 32939967.13 32941889.73 32943770.22 32945611.28 32947415.32 32949184.48 Gap -0.02% -0.01% -0.01% 0.00% 0.01% 0.01% 0.02%  c r w C I 4657224.85 4660009.99 4662704.48 4665316.678 4667853.73 4670321.79 4672726.23 Gap -0.17% -0.11% -0.06% 0.00% 0.05% 0.11% 0.16% OV 32935678.39 32938463.54 32941158.03 32943770.22 32946307.27 32948775.33 32951179.77 Gap -0.02% -0.02% -0.01% 0.00% 0.01% 0.02% 0.02%  e r w C I 3292753.68 3750274.68 4207795.68 4665316.678 5122837.68 5580358.68 6037879.68 Gap -29.42% -19.61% -9.81% 0.00% 9.81% 19.61% 29.42% OV 31571207.22 32028728.22 32486249.22 32943770.22 33401291.22 33858812.22 34316333.22 Gap -4.17% -2.78% -1.39% 0.00% 1.39% 2.78% 4.17%  h r w C I 4650403.35 4655668.91 4660623.12 4665316.678 4669787.84 4674066.35 4678175.84 Gap -0.32% -0.21% -0.10% 0.00% 0.10% 0.19% 0.28% OV 32928856.89 32934122.46 32939076.67 32943770.22 32948241.39 32952519.89 32956629.38 Gap -0.05% -0.03% -0.01% 0.00% 0.01% 0.03% 0.04%  g r w C R 11058612.56 11245420.25 11431721.41 11617106.42 11802239.42 11987372.42 12172505.42 Gap -4.81% -3.20% -1.60% 0.00% 1.59% 3.19% 4.78% OV 32385276.36 32572084.05 32758385.21 32943770.22 33128903.22 33314036.22 33499169.22 Gap -1.70% -1.13% -0.56% 0.00% 0.56% 1.12% 1.69%  f r w C R 11146187.90 11303872.70 11461479.65 11617106.42 11770266.78 11922189.48 12073560.78 Gap -4.05% -2.70% -1.34% 0.00% 1.32% 2.63% 3.93% OV 32472851.69 32630536.49 32788143.45 32943770.22 33096930.58 33248853.28 33400224.58 Gap -1.43% -0.95% -0.47% 0.00% 0.46% 0.93% 1.39% Remark: in Table6, “OV” means the optimal value, and “Gap” means the change of the optimal value. Note that the gap will be zero if no parameter is changed.From Table6, we can see that optimal values can change significantly due to the change of those parameters. For example, the optimal value will decrease by 6.63% if qw decreases by 30%, and the optimal value will increase by 6.65% if qw increases by 30%. The influences of the parameters on individual costs are summarized as follows: When qw increases,CR will increase,CI and will decrease, butCL andCT will not change. When driw increases,CT andCR will increase butCL andCI will not change. When brw, crw, erw, hrw increase, onlyCI will increase and the other costs will not change. When grw or frw increases, onlyCI will increase. We can see that, in general, driw and qw have the most significant influences. Moreover, the optimal value is more sensitive when driw, qw, and erw are relatively small. In practice, business managers may focus more on those parameters to stabilize the supply chain performance.


## 5.4. Performance Analysis

In this subsection, IDE is compared with Lingo 11 and DE to validate its performance in terms of solution accuracy and CPU time efficiency, and the three approaches are tested on a set of small- (i.e., 20×5×2/3/4, 40×6×2/3/4), medium- (i.e., 50×5/6/7×3, 70×6/7/8×3), and large-sized problems (i.e., 80×8×2/3/4, 100×8/9/10×3). Since population sizeNp is an important parameter in IDE and it may vary according to the problem sizes, its value needs to be set properly to get the best possible performance. According to Figures 3–5 (in which “red line” represents solution accuracy in terms of the percentage of finding the optimal solutions, and “blue line” represents the mean CPU time), IDE will have the best performance whenNp= 3S, 4S, 5S for small-, medium-, and large-sized problems, respectively. Therefore, we will use these settings in the subsequent experiments.Figure 3
Analysis ofNpfor small-sized problems (40×6×4).Figure 4
Analysis ofNpfor medium-sized problems (70×8×3).Figure 5
Analysis ofNpfor large-sized problems (100×10×3).Table7 shows the comparison results between IDE, DE, and Lingo 11. To avoid biased results, IDE and DE were executed 30 times for each test instance. Therefore, in this table, “O.V.” means the optimal value, “O.V. (Best)” and “O.V. (Mean)” represent the best and mean optimal values in 30 runs, respectively, “OS Ratio” means the percentage of finding the optimal solutions, and “S.D.” means the standard deviation of optimal values in 30 runs.Table 7
Lingo 11 versus DE versus IDE.
Instance Lingo 11 DE IDE (T1-T2) / T1 I R | W | O.V. CPU Time O.V. O.V. OS Ratio CPU Time S.D. O.V. O.V. OS Ratio CPU Time S.D. (S) (Best) (Mean) (S) - T1 (Best) (Mean) (S) - T2 20 5 2 6987580 6 6987580.30 6997722.08 21/30 0.14 0.01 6987580.30 6987580.30 30/30 0.14 0.01 -3.95% 3 10355510 15 10355509.01 10358913.98 26/30 0.34 0.05 10355509.01 10355509.01 30/30 0.35 0.03 -2.63% 4 13487440 27 13487441.55 13488007.01 29/30 0.65 0.05 13487441.55 13487441.55 30/30 0.67 0.05 -3.08%  40 6 2 12598710 53 12598709.04 12608013.39 16/30 1.01 0.06 12598709.04 12598750.74 29/30 0.89 0.04 11.68% 3 18642000 135 18642000.91 18643688.34 26/30 2.71 0.08 18642000.91 18642000.91 30/30 2.50 0.08 7.71% 4 24367940 263 24367942.66 24368298.67 26/30 5.71 0.14 24367942.66 24367942.66 30/30 5.19 0.10 9.02%  50 5 3 23313920 123 23313920.93 23316057.84 26/30 4.85 0.32 23313920.93 23313920.93 30/30 5.03 0.21 -3.68% 6 23017930 229 23017934.37 23018947.56 27/30 6.78 0.16 23017934.37 23017934.37 30/30 6.05 0.12 10.76% 7 19471180 326 19471179.29 19471323.45 28/30 9.57 0.17 19471179.29 19471179.29 30/30 7.45 0.16 22.14%  70 6 3 31652380 472 31652384.45 31652401.08 29/30 16.17 0.50 31652384.45 31652384.45 30/30 16.26 0.60 -0.56% 7 30150360 1394 30150355.84 30150535.41 27/30 21.51 0.47 30150355.84 30150355.84 30/30 18.28 0.42 14.99% 8 28295930 1993 28295927.33 28296407.75 25/30 22.89 0.43 28295927.33 28295939.60 29/30 20.85 0.53 8.92%  80 8 2 21939770 721 21939767.66 21939991.65 28/30 13.28 0.24 21939767.66 21939767.66 30/30 12.05 0.31 9.31% 3 32943770 2974 32943770.22 32943973.24 29/30 41.51 1.15 32943770.22 32943770.22 30/30 35.76 0.72 13.83% 4 43234620 >3600 43234618.55 43234779.15 29/30 91.60 1.84 43234618.55 43234618.55 30/30 74.05 1.35 19.15%  100 8 3 39081930 >3600 39081579.93 39081760.66 28/30 115.95 2.07 39081579.93 39081579.93 30/30 69.93 1.45 39.69% 9 38614670 >3600 39767341.71 40015199.70 0/30 102.18 0.90 38614396.07 38614592.75 29/30 83.84 1.72 17.95% 10 36450520 >3600 37621661.96 37827792.61 0/30 124.76 0.74 36450173.57 36450173.57 30/30 88.87 2.74 28.76%From Table7, we can see that IDE is the most effective and efficient method to solve MCLIPs, and its advantages over Lingo 11 and DE are summarized as follows:(1)
The optimal values obtained by IDE and Lingo 11 were almost identical, but IDE has better runtimes than Lingo 11. This indicates that IDE is much more efficient than Lingo 11 with the same solution accuracy.(2)
IDE and DE can achieve similar solution accuracies on small- and medium-sized problems. However, the optimal values from DE are significantly larger than those from IDE on large-sized problems. Particularly, DE cannot find optimal solutions for the problems whose sizes are over (100×9×3). This indicates that IDE has a great global search ability and is more stable and consistent than DE.(3)
IDE has better runtimes than DE in most test instances, and it has a greater advantage on large-sized problems. This indicates that IDE is more efficient than DE.According to the comparison results, we can see that IDE is the best approach to solve MCLIPs compared with DE and Lingo 11.


## 6. Conclusion and Future Research

In this paper, we study a multicommodity location-inventory problem by considering false failure returns in a forward-reverse logistics network. To the best of our knowledge, this is the first work that considers false failure returns to optimize joint location-inventory decisions, and it contributes to the literature because of the great importance of false failure returns in practice.This research can be extended in several directions: First, other types of returns such as defective items can be considered to make this study more practical. Second, the (R,Q) model is adopted to replenish inventories in this work. However, there are other popular inventory control policies in practice, and it will be interesting to explore them and compare their performances in forward-reverse logistics networks. Third, it will be valuable to study location-inventory-routing problems by considering false failure returns.



---

*Source: 1967398-2018-10-24.xml*