Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2018 / Article
Special Issue

Emerging Trends on Optimization and Control under Uncertainty in Transportation and Construction

View this Special Issue

Research Article | Open Access

Volume 2018 |Article ID 4390480 |

Lei Wang, Mark Goh, Ronggui Ding, Vikas Kumar Mishra, "Improved Simulated Annealing Based Network Model for E-Recycling Reverse Logistics Decisions under Uncertainty", Mathematical Problems in Engineering, vol. 2018, Article ID 4390480, 17 pages, 2018.

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

Academic Editor: Frederico R. B. Cruz
Received21 Sep 2018
Revised29 Nov 2018
Accepted03 Dec 2018
Published18 Dec 2018


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.

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) [68], 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 [1417]. 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 [2931].

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. Section 2 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: the identification of RLN; the assessment of RLN; the static analysis of RLN; 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].

In phase , 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 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 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 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.

is defined as a binary and square matrix where if node transports product to node ; otherwise, it is blank. Figure 2 gives an example on the use of to represent a binary RLN [40, 41].

The following notations are used in the formulation of the model:

is total number of nodes in RLN

is total number of functions

is number of simulation periods

is set of nodes in RLN

is set of functions

is set of simulation periods in a given planning horizon, and the value of parameters related to denotes the initial state of the model

is binary and square matrix of RLN,

is matrix of transportation capacity between the nodes,

is matrix of transportation cost between the nodes,

is matrix of function network,

is matrix of node function capacity,

is matrix of node function cost,

is matrix of storage volume of specific product that will be processed by function , product for short, in node

is matrix of unit storage cost of product in node

is matrix of transportation relationship between node and node

is storage quantity of product in node at the end of period

is amount of product shipped from node to node in period

is amount of product shipped from other nodes to node in period

is amount of product shipped from node to other nodes in period

is volume of product increased through predecessor functions of function in node during period

is volume of product decreased through function in node during period

is total function cost in period

is total storage cost in period

is total transportation cost in period

is total opportunity loss in period

is opportunity loss of per unit product

is price of per unite product

is demand for node in period

is benefit in period t,

is total benefit of all the periods

is volume of e-waste collected into node in period

is total volume of e-waste collected in all the periods

is sequence number of distribution function in the set

is sequence number of collection function in the set

is warm-up periods at the start of the simulation process

is number of simulation runs

is warm-up times of simulation

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 matrix related to transportation capacity is put forward, in which refers to the maximal volume of products that can be transported from node to node in a unit period. Analogously, the matrix concerned with the transportation cost within that means the transportation cost of the unit product from node to node . So the matrices and can be obtained by adjusting the parameters that are equal to 1 in the binary matrix of 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 are production, distribution, collection, inspection, recovery, and disposal. Also, the product needs to be processed by the function named product in this study. The relationships between the functions form the function network denoted by the matrix , where is the proportion of product that becomes product .

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 matrix , where notation is the capacity on function of node . Likewise, matrix denotes the function cost of each node, where represents the cost of per unit product produced by function in node .

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 thus is proposed to denote the different storage places for different kinds of products in the nodes, where is the volume of product stored in node . Hence, the matrix is used to denote the storage cost and is the unit storage cost of product in node during one period.

After defining the concepts of and product , a specific matrix about the transportation relationship between two nodes is needed. The binary matrix describes the transportation logic among the nodes and presents the transportation amount, where is the proportion of products in node that delivered to node . As to the type of product delivered, it will depend on the initial function(s) of node , which can be determined by and . The calculation process of transportation volume of product shipped from node to node will be discussed next.

3.2. Mathematical Relationships between Matrices

Hybrid nodes make RLN complex as product in node may be transported to node within function or deal with node itself [10, 42, 43], except the situation that node does not have any function but just works as a storage point. Taking this issue into account, assume that only when node does not have function that the product will be transported to the related node . Maybe node and node both have the function and the volume of product in node exceeds the capacity of node while that of node is not. In this case, product will transport from node to node and the assumption above is updated to “only when node does not have the function or volume of product in node exceeds its capacity of function that product in node will be transported to node ”.

Hence, using the notations and the matrices, the RLN mode can be formulated as follows.

When node does not have function ( or its volume of product exceeds its capacity (), node will transport product to directed node based on the value of , which possesses function or serves as pure storage nodes without any functions and its downstream nodes have function . The amount of the product transported for node to node , i.e., , is stated as follows:And the output amount of product of node during period is found fromSo for the input amount of product from other nodes to node , the equation isWhat is more, the volume of product , i.e., the product that needs to deal with function k, in node will decrease when node has function based on the given capacity . Therefore,Meanwhile, if node has the predecessor functions of function , the volume of product will increase asThe value of 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 in node , and the decrease through function in node . Mathematically, we haveUntil 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 period isThe storage cost in period isThe transportation cost in period isTo determine the opportunity loss, suppose the demand for node follows a known distribution in the period , denoted by . The opportunity loss of per unit of a product and the number of distribution function in function set are denoted by and , respectively, and opportunity loss in period can be calculated as follows:Likewise, the benefit during period depends on the unit price which is given by And the total benefit is the sum of , which works as the evolution formula of RLN model:Suppose the volume of e-waste collected into node also follows a known distribution in period , denoted by . Let the number of collection in set be , so the total collection volume isThe 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: 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 of of simulation iterations is used for getting the stable value and the criterion principle is formulated as [40] is set as 0.01% of the values of [35], which is calculated by the warm-up process in section (c).

(c) Warm-up times of simulation () is used to avoid the accidental event that the criterion principle is satisfied when is small and the stable value of is not obtained. When the simulation model of RLN runs several times whose number is greater than , a temp arithmetic mean of , which informs the order magnitude of the stable value, is obtained to calculate the value of [35].

The stable value of 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., . 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 in node . So the node in SN is the subnode generated from node in RLN, whose sequence in the matrix of SN is equal to . 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.

The matrix of 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 and node function capacity . The generation process of SN also can be regarded as the fusion of matrices of , , and .

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 works as a link between node (FT node) and node , so product from node will pass through node and goes into node , which can deal with the function in node . So the subnodes of node are identified on the basis of the functions of nodes and . Meanwhile, if either node or node is the PS node as well, it is necessary to search the upstream nodes of node and downstream nodes of node , 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 of 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 .

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 node 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 has the potential to give rise to the perturbation to the downstream nodes’ operations and eventually supply disruptions [41]. In this situation, node 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 node for a particular pair of nodes and is defined as the proportion of geodesics, the shortest paths, connecting that pair that passes through [45]. This will be adjusted as follows: nodes and are regarded as start nodes and end nodes, respectively, in SN instead of the arbitrary two nodes. Meanwhile, all the paths that connect and are taken into account instead of the shortest paths. Then, the dependency of node on node denotes the congestion probability of caused by which can be calculated by all the start nodes and corresponding path. And also, the influence of the node while it is stockout is concerned with the end node and its corresponding path from to , which can be calculated by the sum of its proportions of all the end nodes. So the betweenness of node is defined as the sum of the betweenness proportions of for all pairs including start and end nodes. The calculating process of the betweenness of node is thus as follows: first, finding out all the start-end paths whose number is and then calculating the times of node ( passing by the paths. The index of betweenness of node g ( is defined as To account for the quantitative importance of node 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 is normalized as And then quantitative importance of node in RLN is given by

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 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 of . 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 be the changed storage capacity of node that belongs to a finite set whose range is from zero to the value that is nearest to but less than ) where is the unit facility of storage, is the corresponding cost, and is the total investment for increasing storage capacity.

The solution representation comprises a set denoted by [] where is the current capacity of node , composed by the initial capacity and the changing capacity , as shown in Figure 4.

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 () and collection 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 important 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 and , which is a measure of the relative amount calculated by the current storage capacity of each node ():Let be the adjustment profitability of node , stated asAnd will be used to define the solution given aswhere generates a single uniformly distributed random number in the interval .

The capacity of the node responding to maximal in the set of will increase by the value of unit change, , 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 maximum cannot be increased further and those of other nodes are zero. Then, the capacity of the node with maximum will be decreased and the capacity of the node with the second maximum will be increased. What is more, the capacity of the node with minimum may not be reduced as its current capacity is equal to the initial capacity; therefore the node with the second minimum 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 variable , denoted by 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, represents the process of calling the algorithm of RLN simulation model, whose value is under the specific ; is the set of quantitative importance of each node in RLN; represents the total number of iterations that the perturbation (neighborhood search) should repeat at a particular temperature; and are the initial and final temperature, respectively; is the Boltzmann constant used in computing the probability of accepting a worse solution; and is the coefficient of the cooling schedule [29].

Generate an initial solution X by the equal division method;

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.

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 Figures 6 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.

The initial storage capacity is shown in Table 1, respectively [10, 35].


Capacity/1000 Unit44201010491516151613

A 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 temperature , Boltzmann constant , 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:




All the 27 combinations of the parameters are tested and using 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.

For the setting of and , as these two parameters have the same function that increasing the likelihood of avoiding getting stuck at a local optima, the value of is set by experience and the value of is selected modestly. In detail, when testing the combinations of , we find that the best solution is obtained before the temperature reducing to 0.2. For statistical and computational convenience, 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 Table 2, along with the RIN model parameter settings, which are adopted from an existing study; detailed discussions can be found in [51].


360 periods1000(unit)
20 periods100(°C)
100 times0.01(°C)

6.2. Static Analysis

Using the matrices , and , 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 , 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.



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 (×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 Figure 10, 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.

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 Figure 11.

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 Figure 12.

From Figure 12, 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 , 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 Table 4.

Best value ($)NSR(times)

Improved SA10.1710.0610.120.32884413601.55117.94
Standard SA10.159.9710.070.48907469687.90120.63

Table 4 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 , 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 [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. Table 5 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 . 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.

Number of nodes in the problemsAve percentage gap (%)NSR percentage gap (%)


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 () in RLN is calculated by summing up the betweenness values of the nodes ( corresponding to node 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].

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare no conflicts of interest.


This work was supported by AStar grant (R-385-000-049-305), National Natural Science Foundation of China (71572094), and Natural Science Foundation of Shandong Province (ZR2015GM015). The authors acknowledge the anonymous referees for their valuable and constructive criticism on this work, which improved the content substantially.


  1. S. Agrawal, R. K. Singh, and Q. Murtaza, “A literature review and perspectives in reverse logistics,” Resources, Conservation & Recycling, vol. 97, pp. 76–92, 2015. View at: Publisher Site | Google Scholar
  2. P. Bhunia, NEA study: Singapore produces 60,000 tonnes of e-waste annually, only 6% recycled, 2018.
  3. N. E. Agency, NEA To Implement E-waste Management System For Singapore By 2021, 2018.
  4. H. S. Kilic, U. Cebeci, and M. B. Ayhan, “Reverse logistics system design for the waste of electrical and electronic equipment (WEEE) in Turkey,” Resources, Conservation & Recycling, vol. 95, no. 1, pp. 120–132, 2015. View at: Publisher Site | Google Scholar
  5. R. Wojanowski, V. Verter, and T. Boyaci, “Retail-collection network design under deposit-refund,” Computers & Operations Research, vol. 34, no. 2, pp. 324–345, 2007. View at: Publisher Site | Google Scholar
  6. E. Behmanesh and J. Pannek, “The Effect of Various Parameters of Solution Methodology on a Flexible Integrated Supply Chain Model,” Mathematical Problems in Engineering, vol. 2018, Article ID 5935268, 14 pages, 2018. View at: Publisher Site | Google Scholar | MathSciNet
  7. K. Govindan, M. Fattahi, and E. Keyvanshokooh, “Supply chain network design under uncertainty: a comprehensive review and future research directions,” European Journal of Operational Research, vol. 263, no. 1, pp. 108–141, 2017. View at: Publisher Site | Google Scholar | MathSciNet
  8. E. Keyvanshokooh, M. Fattahi, S. M. Seyed-Hosseini, and R. Tavakkoli-Moghaddam, “A dynamic pricing approach for returned products in integrated forward/reverse logistics network design,” Applied Mathematical Modelling: Simulation and Computation for Engineering and Environmental Systems, vol. 37, no. 24, pp. 10182–10202, 2013. View at: Publisher Site | Google Scholar | MathSciNet
  9. G. C. Stevens, “Integrating the supply chain,” International Journal of Physical Distribution & Materials Managemen, vol. 19, no. 8, pp. 3–8, 1989. View at: Google Scholar
  10. M. Fattahi and K. Govindan, “Integrated forward/reverse logistics network design under uncertainty with pricing for collection of used products,” Annals of Operations Research, vol. 253, no. 1, pp. 193–225, 2017. View at: Publisher Site | Google Scholar | MathSciNet
  11. A. Alshamsi and A. Diabat, “A reverse logistics network design,” Journal of Manufacturing Systems, vol. 37, no. 3, pp. 589–598, 2015. View at: Publisher Site | Google Scholar
  12. X. Bing, J. M. Bloemhof-Ruwaard, and J. G. A. J. Van Der Vorst, “Sustainable reverse logistics network design for household plastic waste,” Flexible Services and Manufacturing Journal, vol. 26, no. 1-2, pp. 119–142, 2014. View at: Publisher Site | Google Scholar
  13. M. I. Gomes, A. P. Barbosa-Povoa, and A. Q. Novais, “Modelling a recovery network for WEEE: A case study in Portugal,” Waste Management, vol. 31, no. 7, pp. 1645–1660, 2011. View at: Publisher Site | Google Scholar
  14. M. Khatami, M. Mahootchi, and R. Z. Farahani, “Benders' decomposition for concurrent redesign of forward and closed-loop supply chain network with demand and return uncertainties,” Transportation Research Part E: Logistics and Transportation Review, vol. 79, pp. 1–21, 2015. View at: Publisher Site | Google Scholar
  15. V. de Rosa, M. Gebhard, E. Hartmann, and J. Wollenweber, “Robust sustainable bi-directional logistics network design under uncertainty,” International Journal of Production Economics, vol. 145, no. 1, pp. 184–198, 2013. View at: Publisher Site | Google Scholar
  16. L. J. Zeballos, C. A. Méndez, A. P. Barbosa-Povoa, and A. Q. Novais, “Multi-period design and planning of closed-loop supply chains with uncertain supply and demand,” Computers & Chemical Engineering, vol. 66, pp. 151–164, 2014. View at: Publisher Site | Google Scholar
  17. C. Cheng, M. Qi, Y. Zhang, and L.-M. Rousseau, “A two-stage robust approach for the reliable logistics network design problem,” Transportation Research Part B: Methodological, vol. 111, pp. 185–202, 2018. View at: Publisher Site | Google Scholar
  18. E. Roghanian and P. Pazhoheshfar, “An optimization model for reverse logistics network under stochastic environment by using genetic algorithm,” Journal of Manufacturing Systems, vol. 33, no. 3, pp. 348–356, 2014. View at: Publisher Site | Google Scholar
  19. M. Kiliç, G. Ulusoy, and F. S. Şerifoǧlu, “A bi-objective genetic algorithm approach to risk mitigation in project scheduling,” International Journal of Production Economics, vol. 112, no. 1, pp. 202–216, 2008. View at: Publisher Site | Google Scholar
  20. M. S. Pishvaee, K. Kianfar, and B. Karimi, “Reverse logistics network design using simulated annealing,” The International Journal of Advanced Manufacturing Technology, vol. 47, no. 1-4, pp. 269–281, 2010. View at: Publisher Site | Google Scholar
  21. T. Thadewald and H. Büning, “Jarque-Bera test and its competitors for testing normality—a power comparison,” Journal of Applied Statistics, vol. 34, no. 1-2, pp. 87–105, 2007. View at: Publisher Site | Google Scholar | MathSciNet
  22. A. Ç. Suyabatmaz, F. T. Altekin, and G. Şahin, “Hybrid simulation-analytical modeling approaches for the reverse logistics network design of a third-party logistics provider,” Computers & Industrial Engineering, vol. 70, no. 1, pp. 74–89, 2014. View at: Publisher Site | Google Scholar
  23. B. Ayvaz, B. Bolat, and N. Aydin, “Stochastic reverse logistics network design for waste of electrical and electronic equipment,” Resources, Conservation & Recycling, vol. 104, pp. 391–404, 2015. View at: Publisher Site | Google Scholar
  24. O. Listeş, “A generic stochastic model for supply-and-return network design,” Computers & Operations Research, vol. 34, no. 2, pp. 417–442, 2007. View at: Publisher Site | Google Scholar
  25. X. Qin, X. Liu, and L. Tang, “A two-stage stochastic mixed-integer program for the capacitated logistics fortification planning under accidental disruptions,” Computers & Industrial Engineering, vol. 65, no. 4, pp. 614–623, 2013. View at: Publisher Site | Google Scholar
  26. E. Keyvanshokooh, S. M. Ryan, and E. Kabir, “Hybrid robust and stochastic optimization for closed-loop supply chain network design using accelerated Benders decomposition,” European Journal of Operational Research, vol. 249, no. 1, pp. 76–92, 2016. View at: Publisher Site | Google Scholar | MathSciNet
  27. M. Fattahi, M. Mahootchi, K. Govindan, and S. M. Moattar Husseini, “Dynamic supply chain network design with capacity planning and multi-period pricing,” Transportation Research Part E: Logistics and Transportation Review, vol. 81, pp. 169–202, 2015. View at: Publisher Site | Google Scholar
  28. H. J. Ko and G. W. Evans, “A genetic algorithm-based heuristic for the dynamic integrated forward/reverse logistics network for 3PLs,” Computers & Operations Research, vol. 34, no. 2, pp. 346–366, 2007. View at: Publisher Site | Google Scholar
  29. P. Fu, H. Li, X. Wang, J. Luo, S.-L. Zhan, and C. Zuo, “Multiobjective Location Model Design Based on Government Subsidy in the Recycling of CDW,” Mathematical Problems in Engineering, vol. 2017, Article ID 9081628, 9 pages, 2017. View at: Publisher Site | Google Scholar
  30. J. Bautista and J. Pereira, “Modeling the problem of locating collection areas for urban waste management. An application to the metropolitan area of Barcelona,” Omega , vol. 34, no. 6, pp. 617–629, 2006. View at: Publisher Site | Google Scholar
  31. K. Schweiger and R. Sahamie, “A hybrid Tabu Search approach for the design of a paper recycling network,” Transportation Research Part E: Logistics and Transportation Review, vol. 50, no. 1, p. 98, 2013. View at: Publisher Site | Google Scholar
  32. A. Jayant, P. Gupta, and S. Garg, “Simulation modelling and analysis of network design for closed-loop supply chain: a case study of battery industry,” Procedia Engineering, vol. 97, pp. 2213–2221, 2014. View at: Google Scholar
  33. K. S. Moghaddam, “Fuzzy multi-objective model for supplier selection and order allocation in reverse logistics systems under supply and demand uncertainty,” Expert Systems with Applications, vol. 42, no. 15-16, pp. 6237–6254, 2015. View at: Publisher Site | Google Scholar
  34. K. Govindan, H. Soleimani, and D. Kannan, “Reverse logistics and closed-loop supply chain: a comprehensive review to explore the future,” European Journal of Operational Research, vol. 240, no. 3, pp. 603–626, 2015. View at: Publisher Site | Google Scholar | MathSciNet
  35. L. Wang, Network dynamic analysis for project governance risk, Management Science and Engineering, Shandong University, 2017.
  36. S. Kara, F. Rugrungruang, and H. Kaebernick, “Simulation modelling of reverse logistics networks,” International Journal of Production Economics, vol. 106, no. 1, pp. 61–69, 2007. View at: Publisher Site | Google Scholar
  37. L. F. Dai and X. F. Wang, “Research on mixed intelligent arithmetic of reuse reverse logistics centers' location model,” Advanced Materials Research, vol. 945-949, pp. 3246–3251, 2014. View at: Publisher Site | Google Scholar
  38. D. Steward, “The design structure matrix: A method for managing the deisgn of complex systems,” IEEE Transactions on Engineering Management, vol. 28, no. 1981, p. pp, 1981. View at: Google Scholar
  39. L. C. Freeman, “Centrality in social networks conceptual clarification,” Social Networks, vol. 1, no. 3, pp. 215–239, 1978-1979. View at: Publisher Site | Google Scholar
  40. C. Fang and F. Marle, “A simulation-based risk network model for decision support in project risk management,” Decision Support Systems, vol. 52, no. 3, pp. 635–644, 2012. View at: Publisher Site | Google Scholar
  41. Y. Kim, T. Y. Choi, T. Yan, and K. Dooley, “Structural investigation of supply networks: A social network analysis approach,” Journal of Operations Management, vol. 29, no. 3, pp. 194–211, 2011. View at: Publisher Site | Google Scholar
  42. S. M. Hatefi and F. Jolai, “Robust and reliable forward-reverse logistics network design under demand uncertainty and facility disruptions,” Applied Mathematical Modelling, vol. 38, no. 9-10, pp. 2630–2647, 2014. View at: Publisher Site | Google Scholar | MathSciNet
  43. J. Guo, X. Liu, and J. Jo, “Dynamic joint construction and optimal operation strategy of multi-period reverse logistics network: a case study of Shanghai apparel E-commerce enterprises,” Journal of Intelligent Manufacturing, vol. 28, no. 3, pp. 819–831, 2017. View at: Publisher Site | Google Scholar
  44. Q. Luo and D. Zhong, “Using social network analysis to explain communication characteristics of travel-related electronic word-of-mouth on social networking sites,” Tourism Management, vol. 46, pp. 274–282, 2015. View at: Publisher Site | Google Scholar
  45. J. Scott, “Popularity, Mediation and Exclusion,” in Social Network Analysis, M. Steele, Ed., pp. 95–112, SAGE, London, UK, 2017. View at: Google Scholar
  46. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” The Journal of Chemical Physics, vol. 21, no. 6, pp. 1087–1092, 1953. View at: Publisher Site | Google Scholar
  47. S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” Science, vol. 220, no. 4598, pp. 671–680, 1983. View at: Publisher Site | Google Scholar | MathSciNet
  48. H. R. Krikke, A. Van Harten, and P. C. Schuur, “Business case Océ: reverse logistic network re-design for copiers,” OR Spectrum, vol. 21, no. 3, pp. 381–409, 1999. View at: Publisher Site | Google Scholar
  49. P. J. M. van Laarhoven and E. H. L. Aarts, Simulated Annealing: Theory and Applications, Springer, 1987. View at: Publisher Site
  50. S. Bahrami, F. Doulati Ardejani, and E. Baafi, “Application of artificial neural network coupled with genetic algorithm and simulated annealing to solve groundwater inflow problem to an advancing open pit mine,” Journal of Hydrology, vol. 536, pp. 471–484, 2016. View at: Publisher Site | Google Scholar
  51. D.-H. Lee and M. Dong, “Dynamic network design for reverse logistics operations under uncertainty,” Transportation Research Part E: Logistics and Transportation Review, vol. 45, no. 1, pp. 61–71, 2009. View at: Publisher Site | Google Scholar
  52. W. M. Mendenhall, T. L. Sincich, and N. S. Boudreau, Statistics for Engineering and the Sciences, Chapman and Hall/CRC, Boca Raton, FL, USA, 6th edition, 2016.
  53. N. Balakrishnan, T. Colton, B. Everitt, W. Piegorsch, F. Ruggeri, and J. L. Teugels, Wiley StatsRef: Statistics Reference Online, John Wiley & Sons, Ltd, Chichester, UK, 2014.
  54. G. V. Glass, “Testing Homogeneity of Variances,” American Educational Research Journal, vol. 3, no. 3, pp. 187–190, 1966. View at: Publisher Site | Google Scholar
  55. X.-B. Ma, F.-C. Lin, and Y. Zhao, “An adjustment to the Bartlett's test for small sample size,” Communications in Statistics—Simulation and Computation, vol. 44, no. 1, pp. 257–269, 2015. View at: Publisher Site | Google Scholar | MathSciNet

Copyright © 2018 Lei Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

More related articles

 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles