# Improved Simulated Annealing Based Network Model for E-Recycling Reverse Logistics Decisions under Uncertainty

Authors: Lei Wang; Mark Goh; Ronggui Ding; Vikas Kumar Mishra
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Electronic waste recycle (e-recycling) is gaining increasing importance due to greater environmental concerns, legislation, and corporate social responsibility. A novel approach is explored for designing the e-recycling reverse logistics network (RLN) under uncertainty. The goal is to obtain a solution, i.e., increasing the storage capacity of the logistics node, to achieve optimal or near-optimal profit under the collection requirement set by the government and the investment from the enterprise. The approach comprises two parts: a matrix-based simulation model of RLN formed for the uncertainty of demand and reverse logistics collection which calculates the profit under a given candidate solution and simulated annealing (SA) algorithm that is tailored to generating solution using the output of RLN model. To increase the efficiency of the SA algorithm, network static analysis is proposed for getting the quantitative importance of each node in RLN, including the static network generation process and index design. Accordingly, the quantitative importance is applied to increase the likelihood of generating a better candidate solution in the neighborhood search of SA. Numerical experimentation is conducted to validate the RLN model as well as the efficiency of the improved SA.


---

## Body

## 1. Introduction

Electronic waste recycling (e-recycling) is gaining notice in all organizations due to a combination of environmental, legal, and social factors [1]. Typically, each person disposes of 11 kg of electronic waste (e-waste), or the equivalent of 73 mobile phones each year, while almost 60% of the population do not know how to recycle the e-waste in Singapore [2]. The National Environmental Agency (NEA) has embarked on the Extended Producer Responsibility program to build an e-waste recycling management system by 2021 [3], which is a reverse logistics collection initiative [4]. In the framework, the producers of the covered electrical and electronic equipment are responsible for collecting and properly treating the e-waste under the standard raised by NEA [3]. Therefore, the e-recycling problem for firms is on how to design the RL network, e.g., increasing the storage capacity of each node and maximizing the firm’s profit under a specified deposit–refund threshold set by the government [5].The design of an RLN must follow from the forward logistics network (FLN) [6–8], which is a “system whose constituent parts include material suppliers, production facilities, distribution services and customers linked together by feed forward flow of materials and feedback flow of information” [9]. By convention in logistics, the nodes that perform the roles of collection and inspection in RLN may also undertake the function of distribution in FLN [4, 10]. Hence, a firm’s revenue and sales volume are influenced by the e-recycling as the collection and the subsequent activities may impact the FLN’s flexibility of meeting the customer demand [10, 11]. For the purpose of maximizing the enterprise profit, the demand should be considered when designing the RLN. It is clear that both the customer demand of the electronic product and the reverse logistics collection of e-waste are uncertain [1, 10], which make the RLN design problem be a classical stochastic optimization problem.It is reported in the literature that typically mixed integer linear programming (MILP) is used for designing and modeling the RLN under the deterministic environment [7, 12, 13]. Stochastic MILP model is developed considering the uncertainty in demand and in reverse logistics collection volume [14]. For solving the stochastic MILP model, two approaches are commonly used. (a) The first approach is two-stage stochastic programming, which is adopted by most of the researchers in this area, where the decision variables are categorized into two groups: the first and second stages [14–17]. For instance, Khatami et al. [14] designed an RLN and integrated it into an existing multiproduct forward network, which was redesigned as well. In their study, the first stage decision variables are strategic and independent from scenarios; on the contrary, second stage decision variables are tactical and scenario-based. (b) The second approach is converting the stochastic MILP model into an equivalent deterministic model under a given confidence level, which can be easily solved by classical technology [18]. All the approaches mentioned above have two basic parts: the evaluation and optimization of the RLN design solution [1]. The evaluation part describes reverse logistics situation and computes objective function under a given solution [19, 20], such as the collection cost [21], logistics cost [8, 22], and enterprise profit [23]. The optimization part selects a suitable solution based on the output of the evaluation part and generates a new candidate solution. The optimization part is usually according to the solution techniques, which include the branch-and-cut [24, 25], branches-and-bound [25], and Benders’ decomposition [14, 26] in the standard mathematical optimization techniques of MILP and simulated annealing (SA) algorithm [20, 27], genetic algorithm [18, 28], and swarm optimization algorithm [29] from the heuristic methods. As the RLN design problem belongs to a class of combinatorial optimization problems that require computational resources at an exponentially growing rate when the number of decision variables increases, hence, heuristic methods are usually adopted when the number of decision variables is large [29–31].Simulation is an effective tool to address the stochastic problems, which has been tested for having better “goodness-of-fit” to the real-life circumstances [1, 32]. This study thus builds a simulation model of capturing the e-cycling RLN given the uncertainty in demand and e-waste collection, which serves as the evaluation part [10, 33]. The RLN design solution changes the capacity of the nodes in RLN, under which the value of the objective function is calculated by the simulation model. So for the optimization part, as the combination of the decision variables is tremendous in this study, simulated annealing (SA) algorithm, which is widely used in the field of logistics network design [1, 34], is selected to function the optimization role. SA generates candidate solutions and tests them based on the evaluation results from the RLN simulation model. For obtaining the evaluation result, sufficient runs of the simulation model need to be accomplished until the output become stable enough. Therefore, integrating the simulation model and the SA algorithm is costly on the respect of computing time [35], which leads to the requirement of improving the SA. Nevertheless, even though heuristic methods are commonly used to obtain a good feasible solution in the design of RLN, researchers care for the application of heuristics to a specific problem and pay less attention to the improvement of them [36, 37]. To address this problem, following the social network analysis, network static analysis is applied to enhance the performance of the SA algorithm, which is the main academic contribution of this paper. Specifically, on the basis of the proposed network generation process, the RLN is converted to a static network whose nodes are the same in terms of the attribute, which is the foundation of conducting the network static analysis. Meanwhile, a new network index is put forward by referencing the concept of betweenness in the theory of social network analysis. Then the quantitative importance of the nodes in RLN is calculated based on the network index value of the nodes in the static network. Accordingly, when generating a new candidate solution in the neighborhood search part of SA, the storage capacity of an RLN node with a higher importance is more likely to be increased.The rest of this paper is structured as follows. Section2 provides the optimization framework of designing the RLN. The simulation model of RLN is presented in Section 3. Section 4 discusses the network static analysis formed by the network generation and index design of the network, and the improved SA algorithm is described in Section 5. Computational results are reported in Section 6, and, finally, a summary of this study and some possible future works are given in Section 7.


## 2. The Optimization Framework for RLN Design

The optimization framework for RLN design consists of four phases:(1) the identification of RLN; (2) the assessment of RLN; (3) the static analysis of RLN; (4) the dynamic analysis of RLN. Figure 1 illustrates this framework. It is noted that, in the general structure of reverse logistics, products from the producers are delivered to a quantity of geographically dispersed customers through forward activities and the returned e-waste are taken from customers and transported to the nodes with specific function, e.g., recovery and disposal, depending on the decision taken to either recapture value or dispose of it [1].Figure 1
The optimization framework for e-recycling RLN design.In phase(1), the nodes are identified based on the logistics process and geographic position. Specifically, the nodes those pertaining to the reverse process as well as the forward process are taken into account [10]. The transportation between different geographic positions causes transport time and cost. As the hybrid nodes within more than one function are considered in this study, so after identifying the functions of a node, the function interactions within a node need to be identified. This phenomenon can be regarded as the inner interaction of a node. Correspondingly, the external interaction is the transportation relationship between one pair of nodes.In phase(2) of the network assessment, the capacity and cost of specific function of each node as well as those of transportation between the nodes are evaluated, which are termed as the node function evaluation and node interaction evaluation, respectively. Meanwhile, the function interaction is evaluated as quantitative sequence between two functions. Therefore, RLN is described by the matrices, which is underlying the simulation model of RLN.Phase(3) conducted the network static analysis to gain the nodes quantitative importance. This analysis has two steps: one describes the static network (SN) by the adjacent matrix and another is the design of the network index. Thus, a matrix generation process is implied to obtain an SN from the matrices of RLN and index is put forward based on the e-recycling characteristic and social network analysis theory. Accordingly, the importance of each node is calculated eventually.The network dynamic analysis of phase(4) consists of three activities: (a) the construction of RLN simulation model; (b) the improvement of SA algorithm using the output of network static analysis; (c) the dynamic analysis of RLN after combining the improved SA and RLN model, and accordingly the network design solution as well as sensitivity analysis result is obtained.


## 3. The Matrix-Based Simulation Model of RLN

The design structure matrix method introduced by Steward [38] is useful for representing and analyzing the relations among system components, which is similar to the adjacent matrix in the social network analysis [39] and complex network analysis [35, 40]. In this paper, the matrix is used to describe RLN and extended to presenting the attributes of the node and arc as well. Mathematical relationships among the matrices and ultimately a matrix-based simulation model of RLN will be presented in this section.M = m i j is defined as a binary and square matrix where mij=1 if node ni transports product to node nj; otherwise, it is blank. Figure 2 gives an example on the use of M to represent a binary RLN [40, 41].Figure 2
Illustration of binary RLN.The following notations are used in the formulation of the model:N is total number of nodes in RLNK is total number of functionsT is number of simulation periodsN = 1,2 , … , i , … , N is set of nodes in RLNK = 1,2 , … , k , … , K is set of functionsT = 0,1 , … , t , … , T is set of simulation periods in a given planning horizon, and the value of parameters related to t=0 denotes the initial state of the modelM = m i j is binary and square matrix of RLN, ∀i,j∈N,j≠iM T A = m i j t a is matrix of transportation capacity between the nodes, ∀i,j∈N,j≠iM T C = m i j t c is matrix of transportation cost between the nodes, ∀i,j∈N,j≠iM F N = m k l f n is matrix of function network, ∀k,l∈K,l≠kM F A = m i k f a is matrix of node function capacity, ∀i∈N,k∈KM F C = m i k f c is matrix of node function cost, ∀i∈N,k∈KM F S = m i k f s is matrix of storage volume of specific product that will be processed by function k, product k for short, in node i,∀i∈N,k∈KM s c = m i k s c is matrix of unit storage cost of product k in node i,∀i∈N,k∈KM T N = m i j t n is matrix of transportation relationship between node i and node j,∀i,j∈N,j≠iF S i k t is storage quantity of product k in node i at the end of period t,∀i∈N,k∈K,t∈TT r i j t k is amount of product k shipped from node j to node i in period t,∀i,j∈N,t∈T,j≠iI i t k is amount of product k shipped from other nodes to node i in period t,∀i∈N,k∈K,t∈TO i t k is amount of product k shipped from node i to other nodes in period t,∀i∈N,k∈K,t∈TU i t ( k ) is volume of product k increased through predecessor functions of function k in node i during period t,∀i∈N,k∈K,t∈TD i k t is volume of product k decreased through function k in node i during period t,∀i∈N,k∈K,t∈TC f t is total function cost in period t,∀t∈TC s t is total storage cost in period t,∀t∈TC p t is total transportation cost in period t,∀t∈TC o t is total opportunity loss in period t,∀t∈Tp 1 is opportunity loss of per unit productp 2 is price of per unite productD i t is demand for node i in period t,∀i∈N,t∈TF t is benefit in period t, ∀t∈TF is total benefit of all the periodsC o l i t is volume of e-waste collected into node i in period t,∀i∈N,t∈TC o l is total volume of e-waste collected in all the periodsk 1 is sequence number of distribution function in the set Kk 2 is sequence number of collection function in the set KW D is warm-up periods at the start of the simulation processN s i m is number of simulation runsW T S is warm-up times of simulationT h r e s h o l d is accuracy degree of stable value


### 3.1. Matrices about the Attributes of Node and Arc in RLN

The binary RLN shows the basic transportation relationships between each pair of nodes. From this, parameters regarding the attributes of the arc are taken into account, such as transport capacity and cost. Therefore, the matrixMTA=mijta related to transportation capacity is put forward, in which mijta refers to the maximal volume of products that can be transported from node i to node j in a unit period. Analogously, the matrix MTC=mijtc concerned with the transportation cost within mijtc that means the transportation cost of the unit product from node i to node j. So the matrices MTA and MTC can be obtained by adjusting the parameters that are equal to 1 in the binary matrix of M to the specific capacity and cost, respectively.After defining the attribute of the arc, the attribute of another key component of RLN, i.e., node, should be considered, which mainly depends on the logistics functions. Agrawal et al. make a summary of the reverse logistics functions as reuse, repair, remanufacturing, recycling, and disposal based on the previous study [1]. The functions of the forward process are mainly related to production and distribution [10]. As for e-recycling, the products returning from customs, i.e., e-waste, will deal with the function of disposal largely [3]. Accordingly, the functions that repair and remanufacture are regarded as recovery. So the functions included in the set K are production, distribution, collection, inspection, recovery, and disposal. Also, the product needs to be processed by the function k named product k in this study. The relationships between the functions form the function network denoted by the matrix MFN=mklfn, where mklfn is the proportion of product k that becomes product l.The node attributes are defined as function capacity, function cost, and function storage. Specifically, function capacity is the highest rate at which goods are produced. The function capacity of all the nodes in RLN is represented by a matrixMFA=mikfa, where notation mikfa is the capacity on function k of node i. Likewise, matrix MFC=mikfc denotes the function cost of each node, where mikfc represents the cost of per unit product produced by function k in node i.The storage at the node should also be categorized based on the function as the node sends different products to responding nodes. Except this, there are hybrid nodes in the network generating various kinds of products in the node consequently, e.g., collection and inspection [10, 42]. The storage matrix MFS=mikfs thus is proposed to denote the different storage places for different kinds of products in the nodes, where mikfs is the volume of product k stored in node i. Hence, the matrix MSC=miksc is used to denote the storage cost and miksc is the unit storage cost of product k in node i during one period.After defining the concepts ofMFS and product k, a specific matrix MTN=mijtn about the transportation relationship between two nodes is needed. The binary matrix M describes the transportation logic among the nodes and MTN presents the transportation amount, where mijtn is the proportion of products in node i that delivered to node j. As to the type of product delivered, it will depend on the initial function(s) of node j, which can be determined by MFA and MFN. The calculation process of transportation volume of product k shipped from node i to node j will be discussed next.


### 3.2. Mathematical Relationships between Matrices

Hybrid nodes make RLN complex as productk in node i may be transported to node j within function k or deal with node i itself [10, 42, 43], except the situation that node i does not have any function but just works as a storage point. Taking this issue into account, assume that only when node i does not have function k that the product k will be transported to the related node j. Maybe node i and node j both have the function k and the volume of product k in node i exceeds the capacity of node i while that of node j is not. In this case, product k will transport from node i to node j and the assumption above is updated to “only when node i does not have the function k or volume of product k in node i exceeds its capacity of function k that product k in node i will be transported to node j”.Hence, using the notations and the matrices, the RLN mode can be formulated as follows.When nodei does not have function k (mikfc=0) or its volume of product k exceeds its capacity (FSikt>mikfa), node i will transport product k to directed node j based on the value of mijtn, which possesses function k or serves as pure storage nodes without any functions and its downstream nodes have function k. The amount of the product transported for node i to node j, i.e., Trijtk, is stated as follows:(1)Trijtk=FSikt×mijtn.And the output amount of product k of node i during period t is found from(2)Oitk=∑j=1,j≠iNTrijtk=∑j=1,j≠iNFSikt×mijtn.So for the input amount of product k from other nodes to node i, the equation is(3)Iitk=∑j=1,j≠iNTrjitk=∑j=1,j≠iNFSjkt×mjitn.What is more, the volume of product k, i.e., the product that needs to deal with function k, in node i will decrease when node i has function k based on the given capacity mikfa. Therefore,(4)Ditk=min⁡FSikt,mikfa.Meanwhile, if node i has the predecessor functions of function k, the volume of product k will increase as(5)Uitk=∑l=1,l≠kKmin⁡FSilt,milcf⁡.The value of FSikt will change as time going on depending on the four processes: the input from the other nodes, the output to the other nodes, the increase through predecessor functions of function k in node i, and the decrease through function k in node i. Mathematically, we have(6)FSikt+1=FSikt-Oitk+Iitk-Ditk+Uitk.Until now, the basic running mechanism of RLN model is developed using the matrices and relationships mentioned above. Subsequently, parameters and corresponding mathematical relationships will be put forward to assess the model and end in the evaluation formula of the model, which is related to the objective function and works as the key interface between RLN model and SA algorithm.Two factors, cost and benefit, will be considered to assess the model and the other factors can be subsumed under these two factors [35]. Besides the cost incurred through function, storage, and transportation, opportunity loss is considered as the e-recycling disturbs the product flow of the forward logistics [4, 10]. What is more, there are lots of hybrid nodes that distribute products to the customer and collect e-waste from the customer; the big amount of e-waste may lead to more opportunity loss as e-waste occupies the storage place of products [4], which make the opportunity loss significant in the design of e-recycling RLN.So the function cost in periodt is(7)Cft=∑n=1N∑k=1Kmikfc×min⁡FSikt-1,mikfa.The storage cost in period t is(8)Cst=∑i=1N∑k=1Kmiksc×FSikt-1+FSikt2.The transportation cost in period t is(9)Cpt=∑i=1N∑j≠iN∑k=1Kmijtc×Trijtk.To determine the opportunity loss, suppose the demand for node i follows a known distribution in the period t, denoted by Dit. The opportunity loss of per unit of a product and the number of distribution function in function set K are denoted by p1 and k1, respectively, and opportunity loss in period t can be calculated as follows:(10)Cot=∑i=1Np1×Dit-min⁡FSik1t,mik1fa,Dit⁡.Likewise, the benefit during period t depends on the unit price p2 which is given by (11)Ft=∑i=1Np2×min⁡FSik1t,mik1fa,Dit-Cft-Cst-Cpt-Cot.And the total benefit F is the sum of Ft, which works as the evolution formula of RLN model:(12)F=∑t=1TFt.Suppose the volume of e-waste collected into node i also follows a known distribution in period t, denoted by Colit. Let the number of collection in set K be k2, so the total collection volume Col is(13)Col=∑t=1T∑i=1Nmin⁡mik2fa,Colit.The model of RLN is programmed on MATLAB with due considerations given to the following:(a) Warm-up duration (WD) is the volume of periods at the start of one simulation process, after which the duration of simulation will start: WD is concerned with the complexity and size of RLN.(b) RLN simulation model is stochastic because of the uncertainty of demand and collection, so the arithmetic mean ofF of simulation iterations is used for getting the stable value and the criterion principle is formulated as [40](14)∑q=1NsimFqNsim-∑q=1Nsim-1FqNsim-1<Threshold,whenNsim≥2.Threshold is set as 0.01% of the values of F [35], which is calculated by the warm-up process in section (c).(c) Warm-up times of simulation (WTS) is used to avoid the accidental event that the criterion principle is satisfied when Nsim is small and the stable value of F is not obtained. When the simulation model of RLN runs several times whose number is greater than WTS, a temp arithmetic mean of F, which informs the order magnitude of the stable value, is obtained to calculate the value of Threshold [35].The stable value ofF will be used in the neighborhood search of SA to find a candidate solution on capacity settings which will be discussed in the section regarding SA algorithm thereinafter.


## 3.1. Matrices about the Attributes of Node and Arc in RLN

The binary RLN shows the basic transportation relationships between each pair of nodes. From this, parameters regarding the attributes of the arc are taken into account, such as transport capacity and cost. Therefore, the matrixMTA=mijta related to transportation capacity is put forward, in which mijta refers to the maximal volume of products that can be transported from node i to node j in a unit period. Analogously, the matrix MTC=mijtc concerned with the transportation cost within mijtc that means the transportation cost of the unit product from node i to node j. So the matrices MTA and MTC can be obtained by adjusting the parameters that are equal to 1 in the binary matrix of M to the specific capacity and cost, respectively.After defining the attribute of the arc, the attribute of another key component of RLN, i.e., node, should be considered, which mainly depends on the logistics functions. Agrawal et al. make a summary of the reverse logistics functions as reuse, repair, remanufacturing, recycling, and disposal based on the previous study [1]. The functions of the forward process are mainly related to production and distribution [10]. As for e-recycling, the products returning from customs, i.e., e-waste, will deal with the function of disposal largely [3]. Accordingly, the functions that repair and remanufacture are regarded as recovery. So the functions included in the set K are production, distribution, collection, inspection, recovery, and disposal. Also, the product needs to be processed by the function k named product k in this study. The relationships between the functions form the function network denoted by the matrix MFN=mklfn, where mklfn is the proportion of product k that becomes product l.The node attributes are defined as function capacity, function cost, and function storage. Specifically, function capacity is the highest rate at which goods are produced. The function capacity of all the nodes in RLN is represented by a matrixMFA=mikfa, where notation mikfa is the capacity on function k of node i. Likewise, matrix MFC=mikfc denotes the function cost of each node, where mikfc represents the cost of per unit product produced by function k in node i.The storage at the node should also be categorized based on the function as the node sends different products to responding nodes. Except this, there are hybrid nodes in the network generating various kinds of products in the node consequently, e.g., collection and inspection [10, 42]. The storage matrix MFS=mikfs thus is proposed to denote the different storage places for different kinds of products in the nodes, where mikfs is the volume of product k stored in node i. Hence, the matrix MSC=miksc is used to denote the storage cost and miksc is the unit storage cost of product k in node i during one period.After defining the concepts ofMFS and product k, a specific matrix MTN=mijtn about the transportation relationship between two nodes is needed. The binary matrix M describes the transportation logic among the nodes and MTN presents the transportation amount, where mijtn is the proportion of products in node i that delivered to node j. As to the type of product delivered, it will depend on the initial function(s) of node j, which can be determined by MFA and MFN. The calculation process of transportation volume of product k shipped from node i to node j will be discussed next.


## 3.2. Mathematical Relationships between Matrices

Hybrid nodes make RLN complex as productk in node i may be transported to node j within function k or deal with node i itself [10, 42, 43], except the situation that node i does not have any function but just works as a storage point. Taking this issue into account, assume that only when node i does not have function k that the product k will be transported to the related node j. Maybe node i and node j both have the function k and the volume of product k in node i exceeds the capacity of node i while that of node j is not. In this case, product k will transport from node i to node j and the assumption above is updated to “only when node i does not have the function k or volume of product k in node i exceeds its capacity of function k that product k in node i will be transported to node j”.Hence, using the notations and the matrices, the RLN mode can be formulated as follows.When nodei does not have function k (mikfc=0) or its volume of product k exceeds its capacity (FSikt>mikfa), node i will transport product k to directed node j based on the value of mijtn, which possesses function k or serves as pure storage nodes without any functions and its downstream nodes have function k. The amount of the product transported for node i to node j, i.e., Trijtk, is stated as follows:(1)Trijtk=FSikt×mijtn.And the output amount of product k of node i during period t is found from(2)Oitk=∑j=1,j≠iNTrijtk=∑j=1,j≠iNFSikt×mijtn.So for the input amount of product k from other nodes to node i, the equation is(3)Iitk=∑j=1,j≠iNTrjitk=∑j=1,j≠iNFSjkt×mjitn.What is more, the volume of product k, i.e., the product that needs to deal with function k, in node i will decrease when node i has function k based on the given capacity mikfa. Therefore,(4)Ditk=min⁡FSikt,mikfa.Meanwhile, if node i has the predecessor functions of function k, the volume of product k will increase as(5)Uitk=∑l=1,l≠kKmin⁡FSilt,milcf⁡.The value of FSikt will change as time going on depending on the four processes: the input from the other nodes, the output to the other nodes, the increase through predecessor functions of function k in node i, and the decrease through function k in node i. Mathematically, we have(6)FSikt+1=FSikt-Oitk+Iitk-Ditk+Uitk.Until now, the basic running mechanism of RLN model is developed using the matrices and relationships mentioned above. Subsequently, parameters and corresponding mathematical relationships will be put forward to assess the model and end in the evaluation formula of the model, which is related to the objective function and works as the key interface between RLN model and SA algorithm.Two factors, cost and benefit, will be considered to assess the model and the other factors can be subsumed under these two factors [35]. Besides the cost incurred through function, storage, and transportation, opportunity loss is considered as the e-recycling disturbs the product flow of the forward logistics [4, 10]. What is more, there are lots of hybrid nodes that distribute products to the customer and collect e-waste from the customer; the big amount of e-waste may lead to more opportunity loss as e-waste occupies the storage place of products [4], which make the opportunity loss significant in the design of e-recycling RLN.So the function cost in periodt is(7)Cft=∑n=1N∑k=1Kmikfc×min⁡FSikt-1,mikfa.The storage cost in period t is(8)Cst=∑i=1N∑k=1Kmiksc×FSikt-1+FSikt2.The transportation cost in period t is(9)Cpt=∑i=1N∑j≠iN∑k=1Kmijtc×Trijtk.To determine the opportunity loss, suppose the demand for node i follows a known distribution in the period t, denoted by Dit. The opportunity loss of per unit of a product and the number of distribution function in function set K are denoted by p1 and k1, respectively, and opportunity loss in period t can be calculated as follows:(10)Cot=∑i=1Np1×Dit-min⁡FSik1t,mik1fa,Dit⁡.Likewise, the benefit during period t depends on the unit price p2 which is given by (11)Ft=∑i=1Np2×min⁡FSik1t,mik1fa,Dit-Cft-Cst-Cpt-Cot.And the total benefit F is the sum of Ft, which works as the evolution formula of RLN model:(12)F=∑t=1TFt.Suppose the volume of e-waste collected into node i also follows a known distribution in period t, denoted by Colit. Let the number of collection in set K be k2, so the total collection volume Col is(13)Col=∑t=1T∑i=1Nmin⁡mik2fa,Colit.The model of RLN is programmed on MATLAB with due considerations given to the following:(a) Warm-up duration (WD) is the volume of periods at the start of one simulation process, after which the duration of simulation will start: WD is concerned with the complexity and size of RLN.(b) RLN simulation model is stochastic because of the uncertainty of demand and collection, so the arithmetic mean ofF of simulation iterations is used for getting the stable value and the criterion principle is formulated as [40](14)∑q=1NsimFqNsim-∑q=1Nsim-1FqNsim-1<Threshold,whenNsim≥2.Threshold is set as 0.01% of the values of F [35], which is calculated by the warm-up process in section (c).(c) Warm-up times of simulation (WTS) is used to avoid the accidental event that the criterion principle is satisfied when Nsim is small and the stable value of F is not obtained. When the simulation model of RLN runs several times whose number is greater than WTS, a temp arithmetic mean of F, which informs the order magnitude of the stable value, is obtained to calculate the value of Threshold [35].The stable value ofF will be used in the neighborhood search of SA to find a candidate solution on capacity settings which will be discussed in the section regarding SA algorithm thereinafter.


## 4. Network Static Analysis

NSA is used for improving SA algorithm through increasing the probability of generating a better candidate solution in the neighborhood search, which outputs the quantitative importance of the nodes in the respect of maximizing the value of evaluation formula, i.e.,F. Accordingly, this section continues by discussing static network (SN) generation process based on the matrices of RLN and then put forward the index of SN which will be used for calculating the importance of nodes in RLN.


### 4.1. Network Generation Process

Previous studies have tried to find the important nodes in the logistics network, closed-loop supply chain network, etc. by NSA [41], e.g., social network analysis. However, node attributes vary with the number of functions in RLN, which make the NSA short of premises of application as it requires the nodes at one on the attribute [35]. We will divide the nodes in RLN to subnodes until all the nodes are unified on the respect of attribute and the subnodes are regarded as the stored place for specific product k in node i. So the node in SN is the subnode generated from node i in RLN, whose sequence in the matrix of SN is equal to i-1∗N+k.The arc in SN is thus regarded as the directed relationship between the nodes that obey a product flow logic, which can be obtained from the matrices of RLN. Figure 3 illustrates SN by presenting its generation process.Figure 3
The generation process of SN.Notes.→ Product flow in nodes and related data from the matrix MTN. ⇢ Product flow among nodes and related data from the matrix MFN. ▭ Functions and related data from the matrix MFA.The matrix ofMTN reflects the transportation relationships among the nodes denoting the distribution centers, factory, etc. Relatively, in the same node, the link of different functions denotes the sequence of the functions which can be obtained from the function network MFN and node function capacity MFA. The generation process of SN also can be regarded as the fusion of matrices of MTN, MFN, and MFC.There are pure storage nodes in the real word and typical logistics models [1, 10, 44]. Such kind of nodes serve as pure storage (PS) nodes but without any functions in RLN, and the identification of their subnodes in SN is dependent on their connected nodes in RLN. To be more specific, PS node has different places, which are corresponding to subnodes in SN, for the product received from nodes pointing to the PS nodes in RLN, which can be named as the first-tier (FT) node. But the FT nodes may also have the products that will not be delivered to the PS node. The nodes directed by PS node thus are used to judge whether the product of FT node will be delivered to PS node or not. For example, a PS node i works as a link between node h (FT node) and node j, so product from node h will pass through node i and goes into node j, which can deal with the function in node j. So the subnodes of node i are identified on the basis of the functions of nodes h and j. Meanwhile, if either node h or node j is the PS node as well, it is necessary to search the upstream nodes of node h and downstream nodes of node j, respectively.


### 4.2. Network Analysis Index

The node which takes up the important position in the network is more likely to be improved in respect of storage capacity. Meanwhile, individual connections are aggregated to define a measure of connectivity between each node pair, named the network index and reflecting the importance of the node in the network with the specific structure. Hence, as to the SN, the index should be designed based on the capacity setting of the node that impacts the value ofF fundamentally. Specifically, the importance of one node is concerned with the output of the network related to two aspects: one is the extent of the node causing product backlog to the upstream nodes and another is whether the node is likely to lead to the stockout to the downstream nodes or not. If the node with high importance has less capacity, the normal product flow will be severely affected and lead to the high storage cost and opportunity loss and ultimately less F.It is clear that the network is the set of nodes and arcs, but it can also be regarded as the gathering of subsets that denote the paths involving several nodes and corresponding arcs. If nodei is passed through lots of paths, it will have a high probability to be the bottleneck of product flow which makes upstream nodes cannot deliver product to it as usual but store the product by themselves and end in the product backlog. Meanwhile, node i has the potential to give rise to the perturbation to the downstream nodes’ operations and eventually supply disruptions [41]. In this situation, node i works as a significant intermediate node in the network; this is related to the classical issue “betweenness” in the theory of social network analysis.In the social network analysis, the index regarding with betweenness is a measure of the centrality of a node in the network and is normally calculated based on the shortest paths between each pair nodes [45]. The reason for taking the shortest path into account is that the nodes mainly interact through the shortest path in the social network [44, 45], e.g., communication network. What is more, every node may have the chance to spread information to other nodes and the spread stops randomly. But in the SN, the product must start from the start node, such as the production node and collection node, and end in the distribution node or disposal node. Each path from the start node to end node plays a part in SN no matter it is the shortest or not.So based on the idea of betweenness of social network theory and the difference between social network and static network, the calculating index of node quantitative importance is put forward by adjusting the index of betweenness in social network analysis. In the social network analysis, betweenness proportion of nodeg for a particular pair of nodes y and z is defined as the proportion of geodesics, the shortest paths, connecting that pair that passes through g [45]. This will be adjusted as follows: nodes y and z are regarded as start nodes and end nodes, respectively, in SN instead of the arbitrary two nodes. Meanwhile, all the paths that connect y and z are taken into account instead of the shortest paths. Then, the dependency of node y on node g denotes the congestion probability of g caused by y which can be calculated by all the start nodes and corresponding path. And also, the influence of the node g while it is stockout is concerned with the end node z and its corresponding path from g to z, which can be calculated by the sum of its proportions of all the end nodes. So the betweenness of node g is defined as the sum of the betweenness proportions of g for all pairs including start and end nodes. The calculating process of the betweenness of node g is thus as follows: first, finding out all the start-end paths whose number is W and then calculating the times of node g (wg) passing by the paths. The index of betweenness of node g (Bg) is defined as (15)Bg=wgW.To account for the quantitative importance of node i in RLN which is from zero to one and calculated by the sum of the betweenness of its corresponding nodes in SN, betweenness of node g is normalized as (16)Bg′=Bg∑g=1N×KBg.And then quantitative importance Si of node i in RLN is given by(17)Si=∑i-1×K+1≤g≤i×KBg′=∑i-1×K+1≤g≤i×KBg∑g=1N×KBg.


## 4.1. Network Generation Process

Previous studies have tried to find the important nodes in the logistics network, closed-loop supply chain network, etc. by NSA [41], e.g., social network analysis. However, node attributes vary with the number of functions in RLN, which make the NSA short of premises of application as it requires the nodes at one on the attribute [35]. We will divide the nodes in RLN to subnodes until all the nodes are unified on the respect of attribute and the subnodes are regarded as the stored place for specific product k in node i. So the node in SN is the subnode generated from node i in RLN, whose sequence in the matrix of SN is equal to i-1∗N+k.The arc in SN is thus regarded as the directed relationship between the nodes that obey a product flow logic, which can be obtained from the matrices of RLN. Figure 3 illustrates SN by presenting its generation process.Figure 3
The generation process of SN.Notes.→ Product flow in nodes and related data from the matrix MTN. ⇢ Product flow among nodes and related data from the matrix MFN. ▭ Functions and related data from the matrix MFA.The matrix ofMTN reflects the transportation relationships among the nodes denoting the distribution centers, factory, etc. Relatively, in the same node, the link of different functions denotes the sequence of the functions which can be obtained from the function network MFN and node function capacity MFA. The generation process of SN also can be regarded as the fusion of matrices of MTN, MFN, and MFC.There are pure storage nodes in the real word and typical logistics models [1, 10, 44]. Such kind of nodes serve as pure storage (PS) nodes but without any functions in RLN, and the identification of their subnodes in SN is dependent on their connected nodes in RLN. To be more specific, PS node has different places, which are corresponding to subnodes in SN, for the product received from nodes pointing to the PS nodes in RLN, which can be named as the first-tier (FT) node. But the FT nodes may also have the products that will not be delivered to the PS node. The nodes directed by PS node thus are used to judge whether the product of FT node will be delivered to PS node or not. For example, a PS node i works as a link between node h (FT node) and node j, so product from node h will pass through node i and goes into node j, which can deal with the function in node j. So the subnodes of node i are identified on the basis of the functions of nodes h and j. Meanwhile, if either node h or node j is the PS node as well, it is necessary to search the upstream nodes of node h and downstream nodes of node j, respectively.


## 4.2. Network Analysis Index

The node which takes up the important position in the network is more likely to be improved in respect of storage capacity. Meanwhile, individual connections are aggregated to define a measure of connectivity between each node pair, named the network index and reflecting the importance of the node in the network with the specific structure. Hence, as to the SN, the index should be designed based on the capacity setting of the node that impacts the value ofF fundamentally. Specifically, the importance of one node is concerned with the output of the network related to two aspects: one is the extent of the node causing product backlog to the upstream nodes and another is whether the node is likely to lead to the stockout to the downstream nodes or not. If the node with high importance has less capacity, the normal product flow will be severely affected and lead to the high storage cost and opportunity loss and ultimately less F.It is clear that the network is the set of nodes and arcs, but it can also be regarded as the gathering of subsets that denote the paths involving several nodes and corresponding arcs. If nodei is passed through lots of paths, it will have a high probability to be the bottleneck of product flow which makes upstream nodes cannot deliver product to it as usual but store the product by themselves and end in the product backlog. Meanwhile, node i has the potential to give rise to the perturbation to the downstream nodes’ operations and eventually supply disruptions [41]. In this situation, node i works as a significant intermediate node in the network; this is related to the classical issue “betweenness” in the theory of social network analysis.In the social network analysis, the index regarding with betweenness is a measure of the centrality of a node in the network and is normally calculated based on the shortest paths between each pair nodes [45]. The reason for taking the shortest path into account is that the nodes mainly interact through the shortest path in the social network [44, 45], e.g., communication network. What is more, every node may have the chance to spread information to other nodes and the spread stops randomly. But in the SN, the product must start from the start node, such as the production node and collection node, and end in the distribution node or disposal node. Each path from the start node to end node plays a part in SN no matter it is the shortest or not.So based on the idea of betweenness of social network theory and the difference between social network and static network, the calculating index of node quantitative importance is put forward by adjusting the index of betweenness in social network analysis. In the social network analysis, betweenness proportion of nodeg for a particular pair of nodes y and z is defined as the proportion of geodesics, the shortest paths, connecting that pair that passes through g [45]. This will be adjusted as follows: nodes y and z are regarded as start nodes and end nodes, respectively, in SN instead of the arbitrary two nodes. Meanwhile, all the paths that connect y and z are taken into account instead of the shortest paths. Then, the dependency of node y on node g denotes the congestion probability of g caused by y which can be calculated by all the start nodes and corresponding path. And also, the influence of the node g while it is stockout is concerned with the end node z and its corresponding path from g to z, which can be calculated by the sum of its proportions of all the end nodes. So the betweenness of node g is defined as the sum of the betweenness proportions of g for all pairs including start and end nodes. The calculating process of the betweenness of node g is thus as follows: first, finding out all the start-end paths whose number is W and then calculating the times of node g (wg) passing by the paths. The index of betweenness of node g (Bg) is defined as (15)Bg=wgW.To account for the quantitative importance of node i in RLN which is from zero to one and calculated by the sum of the betweenness of its corresponding nodes in SN, betweenness of node g is normalized as (16)Bg′=Bg∑g=1N×KBg.And then quantitative importance Si of node i in RLN is given by(17)Si=∑i-1×K+1≤g≤i×KBg′=∑i-1×K+1≤g≤i×KBg∑g=1N×KBg.


## 5. Improved Simulated Annealing Algorithm

Metropolis et al. first introduced SA in 1953 [46] and it was later popularized by Kirkpatrick et al. [47]. SA has been applied to many hard combinatorial optimization problems, including the design of reverse logistics network and related closed-loop supply chain network [29, 48]. The SA’s neighborhood search method is a random fashion that lets every factor in the candidate solution have the same probability to be changed [49], which increases the probability of local optima and computing time [50]. After the static analysis of SN, the quantitative importance related to the value of F in RLN dynamic model is obtained, which can help to find a better candidate solution efficiently. Therefore, in this section, different parts of the improved SA are introduced, and the relations between them are explained.


### 5.1. Solution Representation and Initial Solution

The capacity improvement of the specific node will increase the profitability of RLN and lead to a better value ofF. The improvement of the node capacity, which is mainly with regard to the storage capacity, is discrete as the unit facility of inventory related to the specific volume of the product. Let ∆Ci be the changed storage capacity of node i that belongs to a finite set ∆C whose range is from zero to the value that is nearest to but less than Invest/(S×Cper) where Cper is the unit facility of storage, S is the corresponding cost, and Invest is the total investment for increasing storage capacity.The solution representationX comprises a set Ci denoted by [C1,C2,…,Ci,…,CN] where Ci is the current capacity of node i, composed by the initial capacity ICi and the changing capacity ∆Ci, as shown in Figure 4.Figure 4
Illustration of solution representation.It is assumed that the enterprise should achieve the standard raised by the government, presented by the volume of collection product [3, 4]. Therefore, the enterprise should put out specific money to increase the capacity of the network. Accordingly, the standard raised by the government and the investment concerned with capacity improvement are the constraint conditions. So the objective function is to maximize the value of total benefit by increasing the capacity of nodes under the constraint of total investment (Invest) and collection standard (Standard):(18)maxF=max∑t=1T∑i=1Np2⁡×min⁡FSik1t,mik1fa,Dit-∑n=1N∑k=1Kmikfc×min⁡FSikt-1,mikfa-∑i=1N∑k=1Kmiksc,F=max×FSikt-1+FSikt2-∑i=1N∑j≠iN∑k=1Kmijtc×Trijtk-∑i=1Np1×Dit-min⁡FSik1t,mik1fa,Dit⁡(19)S.T.∑i=1N∆Ci∗S≤Invest,(20)∑t=1T∑i=1Nmin⁡mik2fa,Colit≥Standard.All the nodes in RLN have the probability of increasing storage capacity, so the initial solution is obtained by increasing the capacity of all the nodes averagely, named as equal division method. Note that the increasing quantity is discrete, so the average may not be an integer, and then the integer less and most near to the average will be taken as the increasing capacity of the nodes and the remainder will be assigned randomly [29].


### 5.2. Definition of the Neighborhood Search

In the standard SA algorithm, the candidate solutions are always defined randomly. NSA generates the quantitative importantSi of each node. So when trying to find a candidate solution in the improved SA, the adjustment of each node’s capacity depends on the difference between Si and Ci′, which is a measure of the relative amount calculated by the current storage capacity of each node (Ci):(21)Ci′=Ci∑i=1NCi.Let Pi be the adjustment profitability of node ni, stated as(22)Pi=1+Si-Ci′2.And Pi′ will be used to define the solution given as(23)Pi′=Pi×rand0,1,where rand(a,b) generates a single uniformly distributed random number in the interval (a,b).The capacity of the node responding to maximal in the set ofP′ will increase by the value of unit change, Sper, and reversely that is related to the minimum will be decreased by the same value.There are two conditions that need to be considered as the increased capacity of the node with the maximumPi′ cannot be increased further and those of other nodes are zero. Then, the capacity of the node with maximum Pi′ will be decreased and the capacity of the node with the second maximum Pi′ will be increased. What is more, the capacity of the node with minimum Pi′ may not be reduced as its current capacity is equal to the initial capacity; therefore the node with the second minimum Pi′ will be searched until the capacity of this node is higher than its initial value.As a consequence, in the proposed solution approach, the improved SA algorithm is applied to swiftly search the feasible domain of the variableICi, denoted by X in the pseudocode of the improved SA in Algorithm 1. This is an extension of the pseudocode presented in [29], which is related to the standard SA. In Algorithm 1, obj(X) represents the process of calling the algorithm of RLN simulation model, whose value is F under the specific X; S=Si is the set of quantitative importance of each node in RLN; Iiter represents the total number of iterations that the perturbation (neighborhood search) should repeat at a particular temperature; T0 and Tf are the initial and final temperature, respectively; K is the Boltzmann constant used in computing the probability of accepting a worse solution; and α is the coefficient of the cooling schedule [29].Algorithm 1:Pseudocode of proposed SA improved by the output of NSA.
Begin Generate an initial solution X by the equal division method; X b e s t = X ; F X = o b j X ; F b e s t = F X; T = T 0; w h i l e   T > T f I = 0; w h i l e   I < I i t e r r 1 = r a n d o m 0,1; G e n e r a t e   a   n e w   s o l u t i o n   Y   f r o m   X   t h r o u g h   t h e   i m p r o v e d   n e i g h b o r h o o d   s e a r c h   m e t h o d   w i t h i n   S; F Y = o b j Y; Δ = F Y - F X; i f   Δ ≥ 0 X = Y ; F X = F Y ; e l s e r = r a n d o m 0,1; i f   r < exp ⁡ Δ / ( K × T ) X = Y ; F X = F Y ;  i f   F X ≥ F b e s t X b e s t = X ; F b e s t = F X ; I = I + 1;  T = T × α;  e n d


## 5.1. Solution Representation and Initial Solution

The capacity improvement of the specific node will increase the profitability of RLN and lead to a better value ofF. The improvement of the node capacity, which is mainly with regard to the storage capacity, is discrete as the unit facility of inventory related to the specific volume of the product. Let ∆Ci be the changed storage capacity of node i that belongs to a finite set ∆C whose range is from zero to the value that is nearest to but less than Invest/(S×Cper) where Cper is the unit facility of storage, S is the corresponding cost, and Invest is the total investment for increasing storage capacity.The solution representationX comprises a set Ci denoted by [C1,C2,…,Ci,…,CN] where Ci is the current capacity of node i, composed by the initial capacity ICi and the changing capacity ∆Ci, as shown in Figure 4.Figure 4
Illustration of solution representation.It is assumed that the enterprise should achieve the standard raised by the government, presented by the volume of collection product [3, 4]. Therefore, the enterprise should put out specific money to increase the capacity of the network. Accordingly, the standard raised by the government and the investment concerned with capacity improvement are the constraint conditions. So the objective function is to maximize the value of total benefit by increasing the capacity of nodes under the constraint of total investment (Invest) and collection standard (Standard):(18)maxF=max∑t=1T∑i=1Np2⁡×min⁡FSik1t,mik1fa,Dit-∑n=1N∑k=1Kmikfc×min⁡FSikt-1,mikfa-∑i=1N∑k=1Kmiksc,F=max×FSikt-1+FSikt2-∑i=1N∑j≠iN∑k=1Kmijtc×Trijtk-∑i=1Np1×Dit-min⁡FSik1t,mik1fa,Dit⁡(19)S.T.∑i=1N∆Ci∗S≤Invest,(20)∑t=1T∑i=1Nmin⁡mik2fa,Colit≥Standard.All the nodes in RLN have the probability of increasing storage capacity, so the initial solution is obtained by increasing the capacity of all the nodes averagely, named as equal division method. Note that the increasing quantity is discrete, so the average may not be an integer, and then the integer less and most near to the average will be taken as the increasing capacity of the nodes and the remainder will be assigned randomly [29].


## 5.2. Definition of the Neighborhood Search

In the standard SA algorithm, the candidate solutions are always defined randomly. NSA generates the quantitative importantSi of each node. So when trying to find a candidate solution in the improved SA, the adjustment of each node’s capacity depends on the difference between Si and Ci′, which is a measure of the relative amount calculated by the current storage capacity of each node (Ci):(21)Ci′=Ci∑i=1NCi.Let Pi be the adjustment profitability of node ni, stated as(22)Pi=1+Si-Ci′2.And Pi′ will be used to define the solution given as(23)Pi′=Pi×rand0,1,where rand(a,b) generates a single uniformly distributed random number in the interval (a,b).The capacity of the node responding to maximal in the set ofP′ will increase by the value of unit change, Sper, and reversely that is related to the minimum will be decreased by the same value.There are two conditions that need to be considered as the increased capacity of the node with the maximumPi′ cannot be increased further and those of other nodes are zero. Then, the capacity of the node with maximum Pi′ will be decreased and the capacity of the node with the second maximum Pi′ will be increased. What is more, the capacity of the node with minimum Pi′ may not be reduced as its current capacity is equal to the initial capacity; therefore the node with the second minimum Pi′ will be searched until the capacity of this node is higher than its initial value.As a consequence, in the proposed solution approach, the improved SA algorithm is applied to swiftly search the feasible domain of the variableICi, denoted by X in the pseudocode of the improved SA in Algorithm 1. This is an extension of the pseudocode presented in [29], which is related to the standard SA. In Algorithm 1, obj(X) represents the process of calling the algorithm of RLN simulation model, whose value is F under the specific X; S=Si is the set of quantitative importance of each node in RLN; Iiter represents the total number of iterations that the perturbation (neighborhood search) should repeat at a particular temperature; T0 and Tf are the initial and final temperature, respectively; K is the Boltzmann constant used in computing the probability of accepting a worse solution; and α is the coefficient of the cooling schedule [29].Algorithm 1:Pseudocode of proposed SA improved by the output of NSA.
Begin Generate an initial solution X by the equal division method; X b e s t = X ; F X = o b j X ; F b e s t = F X; T = T 0; w h i l e   T > T f I = 0; w h i l e   I < I i t e r r 1 = r a n d o m 0,1; G e n e r a t e   a   n e w   s o l u t i o n   Y   f r o m   X   t h r o u g h   t h e   i m p r o v e d   n e i g h b o r h o o d   s e a r c h   m e t h o d   w i t h i n   S; F Y = o b j Y; Δ = F Y - F X; i f   Δ ≥ 0 X = Y ; F X = F Y ; e l s e r = r a n d o m 0,1; i f   r < exp ⁡ Δ / ( K × T ) X = Y ; F X = F Y ;  i f   F X ≥ F b e s t X b e s t = X ; F b e s t = F X ; I = I + 1;  T = T × α;  e n d


## 6. Illustrative Example

In this section, we highlight the application of the proposed approach to a computing example adopted from the classical reverse logistics model and the key issues concerned with the e-recycling problem. The following analysis shows the implementation of the proposed approach, and the effectiveness is indicated by comparing the improved SA and the standard one.The SA program is implemented in MATLAB R2017b on a Windows 10 PC with Intel Core i7 6700 CPU at 3.40 GHz and 16 GB of RAM.


### 6.1. RLN Modeling

The RLN in this study is adapted from [10] and then adjusted based on the characteristic of e-recycling and investigation about Extended Producer Responsibility program in Singapore, such as the high propagation of scrap product in the collected products. Figure 5 shows the binary RLN.Figure 5
The binary RLN.Based on the binary RLN, the matrices related to transportation capacity, cost, and propagation can be obtained as well as the matrices regarding node function capacity, node function cost, and the relationship between the functions. The matrices are shown in Figures6 and 7, respectively [10, 23]. Figure 7 shows that node 09 ships products to node 05 and node 10, respectively, under the proportion of 100%, which means that products to node 05 and to node 10 are different. In other words, node 09 delivers all the products of a specific kind to node 05 and another kinds of products are transported to node 10 totally.Figure 6
The matrix of the function network.Figure 7
The matrix of the transportation network.The initial storage capacityIC is shown in Table 1, respectively [10, 35].Table 1
Initial storage capacities of nodes.
Node 01 02 03 04 05 06 07 08 09 10 Capacity/1000 Unit 44 20 10 10 49 15 16 15 16 13A significant feature of SA is a procedure where the worse solution is accepted and replaces a better solution with a probability and eventually reduces the chance of getting stuck in the local optima. The probability of accepting a worse solution should be high when the temperature is high and reduces in the annealing process. The value of the probability is related to the initial temperatureT0, Boltzmann constant K, and the cooling schedule coefficient α. Therefore, the SA’s parameters tuning is mainly about the joint configuration of these parameters, which is determined by a preliminary experiment. The parameter values tested are as follows:T o = 50,100,150;α = 0.01,0.02,0.03;K = 100,200,300;All the 27 combinations of the parameters are tested and usingTo=100,α=0.01,K=300 seems to satisfy the requirement to the change of the probability of accepting the worse solution, as shown in Figure 8. In Figure 8, the blue full line represents the value of the probability of accepting a worse solution, and the red dot line denotes the trend of the probability in the annealing process.Figure 8
The trend of the probability of accepting a worse solution usingTo=100,α=0.01,K=300.For the setting ofIiter and Tf, as these two parameters have the same function that increasing the likelihood of avoiding getting stuck at a local optima, the value of Iiter is set by experience and the value of Tf is selected modestly. In detail, when testing the combinations of To,α,K, we find that the best solution is obtained before the temperature reducing to 0.2. For statistical and computational convenience, Tf is set as 0.01. Even though this increases the SA computing time, it is not a limiting factor since the redundancy computing time is less than 25 minutes using the software MATLAB on a normal PC.The settings of SA parameters are shown in Table2, along with the RIN model parameter settings, which are adopted from an existing study; detailed discussions can be found in [51].Table 2
Parameter values.
Parameter Value Parameter Value T 360 periods C p e r 1000(unit) W D 20 periods T 0 100(°C) WTS 100 times T f 0.01(°C) Standard 1775(unit) K 300 I n v e s t 600,000$ I i t e r 5 S 12,000$ α 0.01


### 6.2. Static Analysis

Using the matricesMFN,MFA, and MTN, the SN is generated, as shown in Figure 9. In Figure 9, the last two digits denote the function and others denote the node. For example, in “502”, 02 refers to the second function in K, i.e., distribution, and “5” represents the fifth node in RLN. Based on the SN, the quantitative importance of each node is calculated, as shown in Table 3.Table 3
Node quantitative importance calculation.
SN RLN SN RLN Node B g B g ′ Node S i Node B g B g ′ Node S i 101 0.29 0.05 01 0.1 802 0.18 0.03 08 0.03 102 0.29 0.05 902 0.59 0.1 09 0.40 202 0.35 0.06 02 0.06 903 0.35 0.06 302 0.18 0.03 03 0.03 904 0.71 0.12 402 0.18 0.03 04 0.03 905 0.59 0.1 502 0.88 0.15 05 0.15 906 0.12 0.02 602 0.18 0.03 06 0.03 1006 0.12 0.02 10 0.02 702 0.18 0.03 07 703 0.35 0.06 0.15 704 0.35 0.06Figure 9
Static network.


### 6.3. Computational Results and Sensitivity Analysis

The simulation model of RLN and the improved SA are both programmed on MATLAB, and the latter serves as the main program that calls the former. Through running the programs, the optimized storage capacity is obtained, which is[49251515491532164913] (×1000 unit).Based on the optimized node capacity, more values of parameters are computed by running the RIL simulation model. The result is in period-unite, which is helpful and convenient for the decisions regarding different durations.In Figure10, the size of the node reflects the volume of products stored in this node and the color presents the volume of product handled with the function(s) performed by the node.Figure 10
The volume of products stored in and dealt with each node.Node 03, node 07, and node 09 are selected for sensitive analysis. The increased capacity of node 09 is more than that of any other node; the increased capacity of node 07 is only less than node 09; and the increased capacity of node 03 is equal to the average increased capacity. The range of increased capacity of each node is from 1 to 20 (1000 Unit), and the result is shown in Figure11.Figure 11
Profits for different improved capacities of nodes.The profit increases obviously when node 09 capacity improves, but there is no significant distinction on profit when node 03 capacity changes. What is more, the profit presents a decreasing trend when node 07 capacity increases. It is because that node 09 is a hybrid node within the function of collection, inspection, and recovery and thus has a great impact on the e-cycling directly and indirectly. The RLN is facing the situation that lots of e-waste need to be collected and dealt with properly. Increasing the capacity of node 09 contributes to increasing the collection volume as well as the efficiency of dealing with the e-waste. The indirect impact of node 09 is related to Node 07 performing the collection function. The collected products by node 07 will be dealt with by node 09; therefore, if node 07 capacity is just increasing, there will be more collection products stored in node 07 leading to a higher store cost. This is the reason why the profit decreases when node 07 capacity increases. Node 03 does not serve any function related to the reverse logistics; it receives product from its upstream node (node 02) and distributes them to the customers. Increasing the capacity of node 03 is useless when it cannot receive more products from node 02.


### 6.4. Contrastive Analysis of Improved and Standard SAs

For the purpose of indicating the efficiency of improved SA, it is compared with the standard SA which just differs on the neighborhood search. The change in the objective value of each SA as the number of simulation runs increases, which is based on the temperature decreasing, is shown in Figure12.Figure 12
Objective value changing process of each SA.From Figure12, the dotted line regarding standard SA fluctuates at the start of the simulation, which is caused by the pure random neighborhood search style. On the contrary, the full line related to improved SA quickly rises, which means that it generates a better candidate solution at each step, with the help of the node quantitative importance. What is more, the final value of profit, i.e., the best value of F, obtained by the improved SA is higher than that of standard SA, which demonstrates the effectiveness of the improved SA. Meanwhile, the improved SA searches the final profit with less number of running the simulation than the standard SA, which signifies the efficiency of the improved SA.To further verify the above conclusions, both of the standard SA and improved SA run 20 times. For each kind of SA, a sample of the best value and a sample of the number of simulation runs (NSR) are got. The maximum (Max), minimum (Min), average (Ave), and standard deviation (SD) of the samples are described in Table4.Table 4
Descriptive analysis of the samples.
Best value (×104$) NSR(times) Max Min Ave SD Max Min Ave SD Improved SA 10.17 10.06 10.12 0.32 884 413 601.55 117.94 Standard SA 10.15 9.97 10.07 0.48 907 469 687.90 120.63Table4 shows that the improved SA performs better than the standard SA with respect to the effectiveness (best value) as well as the efficiency (NSR). But it is still necessary to judge whether the difference is significant or not in a statistical sense and the significance level is set as 5% by convention.With respect to the samples of best value, we should judge whether the samples obey the normal distribution and have equal variance, which are the foundations of comparing the two samples by parametric test [52]. A two-sided goodness-of-fit test, named Lilliefors test, is adapted to judge whether the sample obeys the normal distribution [53], the p value of the improved SA sample is 0.5000 and that of the standard SA is 0.1182, which means that both of them obey the normal distribution. And then, Bartlett’s test is used to judge whether the two samples have equal variance [54]; the p value is 0.0514, which indicates that the test accepts the null hypothesis that the variances are equal across the two samples. Accordingly, the parametric test is conducted as the two basic assumptions are satisfied. The returned p value is 0.0014, which indicates the two samples are different from the lens of statistics theory [55]. Therefore, it is reasonable to make the conclusion that the improved SA is better than the standard SA on the respect of effectiveness from the lance of statistics theory.By the same token, when comparing the NSR samples, the Lilliefors test shows that both of the samples obey the normal distribution. The p value calculated by Bartlett’s test is 0.9237, which is in favor of accepting the null hypotheses that the variances are equal across the two samples [54]. Therefore, the two basic foundations of conducting the parametric test are satisfied. The parametric test returns the p value as 0.0392 [52], which is in favor of the conclusion that these two samples are different with respect to the average. Therefore, the NSR of the improved SA is less than that of the standard SA based on the statistics theory, which indicates the efficiency of the former.In order to verify the conclusion that the improve SA performs better than standard SA with respect to effectiveness and efficiency, two new problems are designed based on the problem in the illustrative example, which has 10 nodes and named as 10-node problem. One of the new problems is designed within 5 nodes by selecting nodes from the 10-node problem following the principle that all the functions should be involved. Similarly, a 15-node problem is designed by adding five new nodes into the RIN of the 10-node problems. In order to compare the performance of the two SAs on the three problems, two percentage gaps (%) are proposed: Ave percentage gap calculated as (Ave of standard SA – Ave of improved SA)/Ave of standard SA and NSR percentage gap calculated as (NSR of improved SA – NSR of standard SA)/NSR of improved SA. Table5 shows the result, which demonstrates that the improved SA performs better than the standard SA in dealing with all the three problems. The Ave percentage gaps are small; it is because that both of the two SAs running efficient times until the temperature decrease to Tf. Even though the values of Ave percentage gaps is small, the differences between the two SAs is significant, which demonstrate the effectiveness of the improved SA. The NSR percentage gaps of the three problems are all relatively big which demonstrates the efficiency of the improved SA.Table 5
Percentage gaps of Ave and NSR between the two SAs.
Number of nodes in the problems Ave percentage gap (%) NSR percentage gap (%) 5 1.30 8.72 10 0.50 14.36 15 0.70 19.66


## 6.1. RLN Modeling

The RLN in this study is adapted from [10] and then adjusted based on the characteristic of e-recycling and investigation about Extended Producer Responsibility program in Singapore, such as the high propagation of scrap product in the collected products. Figure 5 shows the binary RLN.Figure 5
The binary RLN.Based on the binary RLN, the matrices related to transportation capacity, cost, and propagation can be obtained as well as the matrices regarding node function capacity, node function cost, and the relationship between the functions. The matrices are shown in Figures6 and 7, respectively [10, 23]. Figure 7 shows that node 09 ships products to node 05 and node 10, respectively, under the proportion of 100%, which means that products to node 05 and to node 10 are different. In other words, node 09 delivers all the products of a specific kind to node 05 and another kinds of products are transported to node 10 totally.Figure 6
The matrix of the function network.Figure 7
The matrix of the transportation network.The initial storage capacityIC is shown in Table 1, respectively [10, 35].Table 1
Initial storage capacities of nodes.
Node 01 02 03 04 05 06 07 08 09 10 Capacity/1000 Unit 44 20 10 10 49 15 16 15 16 13A significant feature of SA is a procedure where the worse solution is accepted and replaces a better solution with a probability and eventually reduces the chance of getting stuck in the local optima. The probability of accepting a worse solution should be high when the temperature is high and reduces in the annealing process. The value of the probability is related to the initial temperatureT0, Boltzmann constant K, and the cooling schedule coefficient α. Therefore, the SA’s parameters tuning is mainly about the joint configuration of these parameters, which is determined by a preliminary experiment. The parameter values tested are as follows:T o = 50,100,150;α = 0.01,0.02,0.03;K = 100,200,300;All the 27 combinations of the parameters are tested and usingTo=100,α=0.01,K=300 seems to satisfy the requirement to the change of the probability of accepting the worse solution, as shown in Figure 8. In Figure 8, the blue full line represents the value of the probability of accepting a worse solution, and the red dot line denotes the trend of the probability in the annealing process.Figure 8
The trend of the probability of accepting a worse solution usingTo=100,α=0.01,K=300.For the setting ofIiter and Tf, as these two parameters have the same function that increasing the likelihood of avoiding getting stuck at a local optima, the value of Iiter is set by experience and the value of Tf is selected modestly. In detail, when testing the combinations of To,α,K, we find that the best solution is obtained before the temperature reducing to 0.2. For statistical and computational convenience, Tf is set as 0.01. Even though this increases the SA computing time, it is not a limiting factor since the redundancy computing time is less than 25 minutes using the software MATLAB on a normal PC.The settings of SA parameters are shown in Table2, along with the RIN model parameter settings, which are adopted from an existing study; detailed discussions can be found in [51].Table 2
Parameter values.
Parameter Value Parameter Value T 360 periods C p e r 1000(unit) W D 20 periods T 0 100(°C) WTS 100 times T f 0.01(°C) Standard 1775(unit) K 300 I n v e s t 600,000$ I i t e r 5 S 12,000$ α 0.01


## 6.2. Static Analysis

Using the matricesMFN,MFA, and MTN, the SN is generated, as shown in Figure 9. In Figure 9, the last two digits denote the function and others denote the node. For example, in “502”, 02 refers to the second function in K, i.e., distribution, and “5” represents the fifth node in RLN. Based on the SN, the quantitative importance of each node is calculated, as shown in Table 3.Table 3
Node quantitative importance calculation.
SN RLN SN RLN Node B g B g ′ Node S i Node B g B g ′ Node S i 101 0.29 0.05 01 0.1 802 0.18 0.03 08 0.03 102 0.29 0.05 902 0.59 0.1 09 0.40 202 0.35 0.06 02 0.06 903 0.35 0.06 302 0.18 0.03 03 0.03 904 0.71 0.12 402 0.18 0.03 04 0.03 905 0.59 0.1 502 0.88 0.15 05 0.15 906 0.12 0.02 602 0.18 0.03 06 0.03 1006 0.12 0.02 10 0.02 702 0.18 0.03 07 703 0.35 0.06 0.15 704 0.35 0.06Figure 9
Static network.


## 6.3. Computational Results and Sensitivity Analysis

The simulation model of RLN and the improved SA are both programmed on MATLAB, and the latter serves as the main program that calls the former. Through running the programs, the optimized storage capacity is obtained, which is[49251515491532164913] (×1000 unit).Based on the optimized node capacity, more values of parameters are computed by running the RIL simulation model. The result is in period-unite, which is helpful and convenient for the decisions regarding different durations.In Figure10, the size of the node reflects the volume of products stored in this node and the color presents the volume of product handled with the function(s) performed by the node.Figure 10
The volume of products stored in and dealt with each node.Node 03, node 07, and node 09 are selected for sensitive analysis. The increased capacity of node 09 is more than that of any other node; the increased capacity of node 07 is only less than node 09; and the increased capacity of node 03 is equal to the average increased capacity. The range of increased capacity of each node is from 1 to 20 (1000 Unit), and the result is shown in Figure11.Figure 11
Profits for different improved capacities of nodes.The profit increases obviously when node 09 capacity improves, but there is no significant distinction on profit when node 03 capacity changes. What is more, the profit presents a decreasing trend when node 07 capacity increases. It is because that node 09 is a hybrid node within the function of collection, inspection, and recovery and thus has a great impact on the e-cycling directly and indirectly. The RLN is facing the situation that lots of e-waste need to be collected and dealt with properly. Increasing the capacity of node 09 contributes to increasing the collection volume as well as the efficiency of dealing with the e-waste. The indirect impact of node 09 is related to Node 07 performing the collection function. The collected products by node 07 will be dealt with by node 09; therefore, if node 07 capacity is just increasing, there will be more collection products stored in node 07 leading to a higher store cost. This is the reason why the profit decreases when node 07 capacity increases. Node 03 does not serve any function related to the reverse logistics; it receives product from its upstream node (node 02) and distributes them to the customers. Increasing the capacity of node 03 is useless when it cannot receive more products from node 02.


## 6.4. Contrastive Analysis of Improved and Standard SAs

For the purpose of indicating the efficiency of improved SA, it is compared with the standard SA which just differs on the neighborhood search. The change in the objective value of each SA as the number of simulation runs increases, which is based on the temperature decreasing, is shown in Figure12.Figure 12
Objective value changing process of each SA.From Figure12, the dotted line regarding standard SA fluctuates at the start of the simulation, which is caused by the pure random neighborhood search style. On the contrary, the full line related to improved SA quickly rises, which means that it generates a better candidate solution at each step, with the help of the node quantitative importance. What is more, the final value of profit, i.e., the best value of F, obtained by the improved SA is higher than that of standard SA, which demonstrates the effectiveness of the improved SA. Meanwhile, the improved SA searches the final profit with less number of running the simulation than the standard SA, which signifies the efficiency of the improved SA.To further verify the above conclusions, both of the standard SA and improved SA run 20 times. For each kind of SA, a sample of the best value and a sample of the number of simulation runs (NSR) are got. The maximum (Max), minimum (Min), average (Ave), and standard deviation (SD) of the samples are described in Table4.Table 4
Descriptive analysis of the samples.
Best value (×104$) NSR(times) Max Min Ave SD Max Min Ave SD Improved SA 10.17 10.06 10.12 0.32 884 413 601.55 117.94 Standard SA 10.15 9.97 10.07 0.48 907 469 687.90 120.63Table4 shows that the improved SA performs better than the standard SA with respect to the effectiveness (best value) as well as the efficiency (NSR). But it is still necessary to judge whether the difference is significant or not in a statistical sense and the significance level is set as 5% by convention.With respect to the samples of best value, we should judge whether the samples obey the normal distribution and have equal variance, which are the foundations of comparing the two samples by parametric test [52]. A two-sided goodness-of-fit test, named Lilliefors test, is adapted to judge whether the sample obeys the normal distribution [53], the p value of the improved SA sample is 0.5000 and that of the standard SA is 0.1182, which means that both of them obey the normal distribution. And then, Bartlett’s test is used to judge whether the two samples have equal variance [54]; the p value is 0.0514, which indicates that the test accepts the null hypothesis that the variances are equal across the two samples. Accordingly, the parametric test is conducted as the two basic assumptions are satisfied. The returned p value is 0.0014, which indicates the two samples are different from the lens of statistics theory [55]. Therefore, it is reasonable to make the conclusion that the improved SA is better than the standard SA on the respect of effectiveness from the lance of statistics theory.By the same token, when comparing the NSR samples, the Lilliefors test shows that both of the samples obey the normal distribution. The p value calculated by Bartlett’s test is 0.9237, which is in favor of accepting the null hypotheses that the variances are equal across the two samples [54]. Therefore, the two basic foundations of conducting the parametric test are satisfied. The parametric test returns the p value as 0.0392 [52], which is in favor of the conclusion that these two samples are different with respect to the average. Therefore, the NSR of the improved SA is less than that of the standard SA based on the statistics theory, which indicates the efficiency of the former.In order to verify the conclusion that the improve SA performs better than standard SA with respect to effectiveness and efficiency, two new problems are designed based on the problem in the illustrative example, which has 10 nodes and named as 10-node problem. One of the new problems is designed within 5 nodes by selecting nodes from the 10-node problem following the principle that all the functions should be involved. Similarly, a 15-node problem is designed by adding five new nodes into the RIN of the 10-node problems. In order to compare the performance of the two SAs on the three problems, two percentage gaps (%) are proposed: Ave percentage gap calculated as (Ave of standard SA – Ave of improved SA)/Ave of standard SA and NSR percentage gap calculated as (NSR of improved SA – NSR of standard SA)/NSR of improved SA. Table5 shows the result, which demonstrates that the improved SA performs better than the standard SA in dealing with all the three problems. The Ave percentage gaps are small; it is because that both of the two SAs running efficient times until the temperature decrease to Tf. Even though the values of Ave percentage gaps is small, the differences between the two SAs is significant, which demonstrate the effectiveness of the improved SA. The NSR percentage gaps of the three problems are all relatively big which demonstrates the efficiency of the improved SA.Table 5
Percentage gaps of Ave and NSR between the two SAs.
Number of nodes in the problems Ave percentage gap (%) NSR percentage gap (%) 5 1.30 8.72 10 0.50 14.36 15 0.70 19.66


## 7. Conclusion

This study has explored an approach within a simulation model and improved SA algorithm for addressing the RLN design issue, i.e., increasing the storage capacity of each node of RLN to maximize the firm’s profit. The simulation model abstracts the e-cycling reverse logistics with consideration of the forward logistics. In view of uncertain demand and e-waste collection volume, the simulation model is capable of analyzing the enterprise profit under a given RLN design solution. The SA algorithm is used to select a suitable solution based on the output of the simulation model and generate a new candidate solution. Integrating the simulation model and SA, an optimal solution for designing RLN can be obtained.Further, SA is improved using the network static analysis. In particular, a network generation process is proposed for changing the RLN to a static network, which satisfies the requirement of conducting the network static analysis. A new network index regarding betweenness is designed with the consideration of the characteristics of reverse logistics. The quantitative importance of a node (A) in RLN is calculated by summing up the betweenness values of the nodes (A1,A2,…) corresponding to node A in the static network. Then the quantitative importance is applied in the neighborhood search part of SA, which lets the storage capacity of a node with a higher importance be more likely to be increased. Therefore, the chance of generating a better candidate solution is increased. Our approach of integrating the RLN simulation model and improved SA is well illustrated with an example. Two more examples are generated and discussed as well. The contrastive analysis within three different examples shows that the improved SA is beyond the standard SA on the respects of effectiveness and efficiency.We can refine this research by (a) taking into account other static networks related to transport cost, time and put forward corresponding indexes to make the NSA more accurate and (b) considering control module regarding overstock and stockout problems into the RLN simulation model for the purpose of designing RLN under different operation strategies, which bring the method closer to reality. Our approach can be applied to other network optimization problems including both dynamic network analysis and static network analysis, such as risk interaction network [40] and project governance network [35].



---

*Source: 4390480-2018-12-18.xml*