1.INTRODUCTION
In recent decades, various industrial sectors have been affected by the growth of environmental concerns, introduction of new regulations and the increase in the rate of return as a result of shortened product life cycle. Accordingly, with the intention of passing the obligatory regulations, the manufacturers have been compelled to retrieve and recover products that are at the end of their life cycle (see Ilgin and Gupta, 2010).
Planning for the collection of used products is related to the reverse logistics network design, which includes not only selecting appropriate paths and vehicles for transporting goods but also determining the optimal locations and capacities for opening the collection and recovery centers. In this area, there are some challenging decisions to be made when various products and varying levels of quality are concerned. This is why a wide variety of mathematical optimization models have been published to deal with different situations in network design, which usually considers different stages for reverse logistics networks.
The literature addresses the collection of returned products in three ways: thirdparty logistics (CruzRivera and Ertel, 2009; de Figueiredo and Mayerle, 2008; Govindan et al., 2012; Krikke et al., 2008; Min and Ko, 2008), opening collection centers by the remanufacturer/ manufacturer (Aras and Aksen, 2008; Östlin et al., 2008; Pokharel and Liang, 2012; Tagaras and Zikopoulos, 2008), and the use of retailers (Choi et al., 2013; Lee et al., 2011; Wojanowski et al., 2007).
The majority of the used products do not have any value in terms of functionality (Schultmann et al., 2006); only a reasonable amount of profit is made when some of them are collected and recycled. Imposing penalties for noncollected, endoflife (EOL) products or pollutant limitations generated during the production process are samples of new approaches that have actively been pursued in the past decade. For example, in the Waste Electrical and Electronic Equipment (WEEE) directive, governments oblige manufacturers to take the responsibility for the entire life cycle of their products.
Generally, customers have no motivation to return their EOL products after use (Guide et al., 2003). Unfortunately, there are few papers that have examined the effect of incentive price on the amount of returned products. The rate of return by the customers is related to the incentive price provided by the collectors, which is determined by various factors, such as distance from the collection or recycling centers and scraps quality. A significant contribution that the present research intends to make and sets it apart from the past researches is that it considers the price effect on the probability amount of returns, taking this fact into account that customers do not have an exact price in mind for selling their products. Therefore, in this paper, the maximum offered price is expressed as a triangular fuzzy number.
This study applies such concepts as distance effect, capacity and penalties simultaneously so that in addition to determining the volume of the products that must be purchased during different periods, it can calculate the optimal number, location, and capacity of the centers to reach the intended goals.
In short, this study attempts to illustrate incentivedependent returns in dynamic network design, which provides a more realistic expression of the offered incentives effect on customers' reaction, and includes different stages. Furthermore, we will evaluate a new procedure for solving this issue, which benefits from the linearization and an exact solution approach.
The remainder of this paper is organized as follows. After a brief review of the relevant literature in Section 2, in Section 3, we define the problem and introduce the objective functions in detail. Next, the mathematical formulation for the model is developed. In Section 4, a tabu search based heuristic is described for solving the problem. In Section 5, a numerical example of its occurrence is used to show the validation of the heuristic approach and investigate model applicability, following which the computational results are presented for the algorithm. Finally, the concluding remarks are presented in the last section.
2.LITERATURE REVIEW
Logistics literature is remarkably rich in papers in the context of product recovery and recycling. Fleischmann et al. (2000) reviewed some case studies on logistics network design for product recovery in different industries, after which he identified general characteristics for such logistics networks. Besides, they denote five groups of activities that appear to be recurrent in EOL product recovery networks: collection, inspection/ separation, reprocessing, redistribution, and disposal.
Additionally, the majority of studies on reverse flows formulate discrete facility locationallocation models. For example, Jayaraman et al. (2003) modeled the reverse logistics of hazardous products with a multilevel warehouse location model. The objective of the model is to find the optimal number and location of collection and refurbishing facilities with the corresponding flow of the hazardous products. In their mixedinteger linear programming model, Jayaraman et al. (2003) present the number of returned products at each originating site.
Pishvaee et al. (2009) classified the literature about the logistics network into three sections: forward logistics, reverse logistics, and integrated forward and reverse logistics. Based on their classified literature, there are not many papers that consider the logistics network design in multiperiods (dynamic) environment. In dynamic environment, Min et al. (2006) proposed a dynamic nonlinear mixedinteger programming model for the deterministic logistics network involving both spatial and temporal consolidation of returned products. Furthermore, Ko and Evans (2007) considered a network operated by a thirdparty logistics service provider. They presented a dynamic mixedinteger nonlinear programming model for the simultaneous design of the forward and reverse network.
In the literature, few papers have attempted to show the relationship between incentives and the product return. Some recent articles are de Figueiredo and Mayerle (2008), Aras and Aksen (2008) and Aksen et al. (2009). However, all of them have just considered the location of collection centers and allocation of customer zones to them without appraising cost and revenue issues in multistage nature of reverse logistics network design.
For instance, Aras and Aksen (2008) formulated a mixedinteger nonlinear facility locationallocation model to find both the optimal locations of a predetermined number of collection centers and the optimal incentive values for different return types. In another study, Aksen et al. (2009) provided a bilevel mixedinteger nonlinear programming to determine the location of collection centers and the optimal price offered for product return. Clearly, the price offered by the company influences the quality level and the rate of returned products. Therefore, they assume that reservation price follows the right triangular distribution (RTD).
To structure the literature review on logistics network design and to identify future research avenues, we have classified design problems of logistics network according to six general specifications: general structure of the network, modeling type, objective functions, problem definition, the network stages included and the solution method. We developed a coding system as shown in Table 1, and the recent mathematical models available in the literature have been coded based on this system in Table 2. Categories based on closed loop or opened loop in Table 2 indicate refer to Fleischmann et al. (2000) who expressed the differences between the two networks in their article.
As shown in Table 2, there are few papers that determine the incentive price for product return in reverse logistics network design. These papers have just considered the simplest case of network design, which is comprised of single period, single product, and single stage. Moreover, only the location of collection and inspection centers (CICs) has been determined while other stages of the network have been neglected.
This paper extends the reverse logistics literature by considering the fuzzy relationship between the proposed incentive price and product acquisition. It also considers a dynamic design of reverse logistics with multiproduct, multiquality levels, multistage and limited budget in each period. Furthermore, in this research, a tabu search based heuristic has been provided so as to find a solution for the model. The coding of the model under question is provided in Table 2.
3.PROBLEM DEFINITION
The reverse logistics network stages that are considered in this research have been illustrated in Figure 1. Based on this figure, the used products are collected from customer zones and sent to the collection/inspection centers. Afterwards, the products are delivered to different recovery centers for recycling. The following are the assumptions considered in the design of such a network:

The model is multiperiod.

The model considers multiple products, which have a different but known quality level.

The potential locations of collection/inspection and recovery facilities are known and fixed.

Costs parameters (setup, fixed cost, variables, nonutilized capacity, noncollected returns, transportation, and holding costs) are known for each location, product, and time period.

Holding cost depends on the residual inventory at the end of each period.

Recovery options operate independently.

Customers themselves dispose of used products that are not collected by the collectors.

The budget is only intended for the construction of facilities. It should also be noted that this amount is limited and that the remainder at the end of each period can be used in other periods.
According to Aksen et al. (2009) the product holders' decision to return the products can be modeled by using the notion of consumer surplus. The product holder i which has the used product type p with quality level q at time t will make a return if the collection center offers a unit incentive B_{pqt} that is at least as large as a reservation price B^{0}_{pqit} . It is assumed that all customers in different locations have the same mental model in responding to the similar offered price. We assume that B^{0}_{pqit} follows the RTD, whose density function is given in (1) and Figure 2. The proportion P_{pqit} of customers in zone i who are willing to return their used products of type p with quality level of q at time t when the collectors offer incentive B_{pqt} per product is calculated from Eq. (2).
It is noted that B^{0}_{pqit} takes values in an interval between 0 and b_{pqit}, which represents the maximum contemplated incentive by a customer for product type p with quality level q. However, considering the deterministic value for the upper bound of this interval may not be suitable since customers do not consider the exact price for the sale of their products. Therefore, we make use of a triangular fuzzy form, defined by b_{pqit,} to show the greatest amount possible.
In the next section, based on the aforementioned characteristics of the network, a fuzzy mixedinteger nonlinear programming model is presented.
3.1.Model Formulation
The following notations are used in the formulation of the mentioned design of reverse logistics network:
Sets
P set of used product types
T set of time periods
I set of fixed customer locations
J set of potential CICs
H set of potential recovery locations
L set of transportation mode
M set of recovery facilities (RFs) outputs
N set of capacity levels available for facilities
Q set of used products' quality levels
Parameters
oc_{jt}^{n} cost of establishing CIC j with capacity level n in period t
fc_{jt} operating fixed cost for CIC j in period t
vc_{pqjt} variable operating cost of product p with quality level q for CIC j in period t
w_{jt}^{n} capacity with level n for CIC j in period t
hc_{jt} inventory holding cost at CIC j in period t
bc_{t} available budget for CICs in period t
or_{ht}^{n} cost of establishing RF h with capacity level n in period t
fr_{ht} operating fixed cost for RF h in period t
vr_{pqht} variable operating cost of product p with quality level q for RF h in period t
pc_{ht} penalty per unit of nonutilized capacity at RF h in period t
car_{ht}^{n} capacity with level n for RF h in period t
s_{mht} sale price of final product m generated at RF h in period t
α_{mpqht} generation fraction of product m from used product p with quality level q at RF h in period t
br_{t} available budget for RFs in period t
d1_{i} distance between customer location i and CIC j
d2_{jh} distance between CIC j and RF h
ct_{lt} handling cost of per unit product traveled in unit distance by transportation mode l in period t
cp_{lt} handling cost of product p using transportation mode l in period t
ca1_{ijlt} flow capacity from customer location i to CIC j using transportation mode l in period t
ca2_{jhlt} flow capacity from CIC j to RF h using transportation mode l in period t
r_{pqit} quantity of used product p with quality level q in customer zone i in period t
uc_{pt} penalty per unit of noncollected product p in period t
b_{pqit} contemplated price by customer i for product p with quality level q in period t
Decision variables
Q1_{pqijlt} quantity of product p with quality level q shipped from customer zone i to CIC j using transportation mode l in period t
Q2_{pqjhlt} quantity of product p with quality level q shipped from CIC j to RF h using transportation mode l in period t
B_{pqt} offered price for used product p with quality level q in period t
X_{jt}^{n} binary variable, it takes 1 if CIC j with capacity level n is installed in period t, and 0 if otherwise
Y_{ht}^{n} binary variable, it takes 1 if RF h with capacity level n is installed in period t, and 0 if otherwise
Dependent variable
I_{pjt} inventory level of used product p collected at CIC j in period t
U_{pit} quantity of uncollected product p from customer i in period t
S_{C}_{t} remaining budget for CICs at the end of period t
S_{r}_{t} remaining budget for RFs at the end of period t
p_{pqit} percent of customers in zone i who are willing to return their product of type p with quality level q in period t
ϕ_{jt}^{n} binary variable, it takes 1 if CIC j with capacity level n is installed in period t or before t, and 0 if otherwise
ξ_{ht}^{n} binary variable, it take 1 if RF h with capacity level n is installed in period t or before of t and 0 otherwise
Based on the aforementioned parameters and indices, the fuzzy mixedinteger nonlinear programming (FMIN LP) model is developed. The objective function is to minimize the total costs minus potential revenues from recovering the used product. The objective function includes the following issues:
Opening and operating cost of collection and inspection centers:
it needs to be mentioned that if the facility opens during a period, it will remain open in the subsequent periods. Establishing the cost of CIC and RF facilities occurs once they are established for the first time. Operating fixed cost is charged in each period after the establishment of the facility.
Opening and operating cost of recovery centers:
Transportation cost:
Revenues from recovering used product plus the remaining budget in the final period:
these revenues are subtracted from costs at objective function.
The objective function is to minimize the summation of Eqs. (3)–(9) minus Eq. (10). The constraints for the model are as follows.
Constraints (11) and (12) ensure that a facility can be established in each location at most in one capacity level. Eqs. (13) and (14) represent capacity limitations on product flow between different nodes. Constraint (15) illustrates the capacity restrictions on recovery facilities in terms of the number of total products entering them. Eq. (16) assures the inventory balance of used products at CICs during periods. Eq. (17) determines the initial inventory of products collected at CICs. Constraint (18) indicates the inventory capacity limitation at CICs. In Eq. (19), it is shown that the quantity of returned products is related to the potential quantity in customer zones. On the other hand, this equation bounds the amount of returned products. Constraint (20) represents possible uttermost quantity of various products that can be collected from customer zones. Eq. (21) shows the relationship between the maximum proportion of returns and the offered incentive. Eqs. (22) to (25) assure budget counterbalance among different periods. Fuzzy Eq. (26) guarantees that the offered prices do not exceed the RTD upper bound. Eqs. (27) and (28) specify facilities that have been opened in the previous periods. Finally, constraints in set (29) enforce the nonnegativity restrictions on the corresponding decision variables and constraints in set (30) enforce the integrality restrictions on the binary variables.
In this paper, a linear ranking function is applied to convert the fuzzy parameter into the crisp equivalent number through the use of the first index of Yager (1978) and Yager (1981). Thus, by applying the index of Ronald and considering the triangular fuzzy number of b_{pqit} = m^{L}_{pqit}, m^{C}_{pqit}, m^{R}_{pqit} , the aforementioned FMINLP problem is transformed into the crisp equivalent problem by replacing constraints (21) and (26) with the following equations.
where d^{L}_{pqit} and d^{R}_{pqit} are the lateral margins (right and left, respectively) of the triangular fuzzy number with central point of m^{C}_{pqit} (Figure 3).
4.SOLUTION APPROACH
Since our model has been known as NPhard, it is difficult to obtain high quality solutions for the relatively largesized examples in a reasonable timeframe using commercial software. Therefore, a heuristic algorithm is developed to solve it.
In this section, we propose an iterative algorithm based on tabu search procedure. Tabu search belongs to a class of local search techniques and is based on avoiding local optima by using memory structures called tabu lists. These lists temporarily record visited solutions and prevent them from being cycled around.
Tabu search performs a search for the solution space by moving from one identified solution to the best solution in a subset of the neighborhood of the current solution. Since there is no necessity to improve the solution in all iterations, a tabu mechanism is provided to prevent the turning process in a sequential series of solutions. This way, the exploration process is not allowed to go back to the previous encountered solutions.
Another feature of the algorithm is called ‘diversification' and ensures the search process is not limited to a restricted part of the solution space. In addition, intensification is identified as yet another characteristic of the approach that consists of a greedy local search around the bestrecognized solutions.
4.1.Proposed Tabu Search Approach
Due to the type of the model and the dependence of variables, the complexity of this problem is extremely high; this is why it is necessary to use a proper method for providing solutions. In this paper, a heuristic approach is applied, which utilizes two procedures, one being approximate and the other one exact, in each iteration.
We use tabu search for providing solutions to the problem and we take advantage of CPLEX libraries to achieve an exact solution at each stage. First, the model is converted to a convex linear form through determining the neighbor solutions for some variables, including φ, ξ (or X and Y) and B by tabu search algorithm (the full explanation is provided in the initial solution). Afterwards, the model is solved through the use of CPLEX software libraries. Finally, the obtained values by CPLEX for objective functions, Q1 and Q2, return to the main algorithm and this process would be repeated to satisfy the stop condition. We terminate the algorithm as soon as the maximum number of iterations limit is met. The algorithm structure can be seen in Figure 4, where the used notations are defined as follows:
Num_Iter number of performed iterations
Max_Iter maximum number of iterations
Num_notchange number of iterations throughout which the incumbent does not improve
Max_notchange maximum amount of iterations throughout which the incumbent does not improve
Obj objective value of newly generated neighboring solution
Best_Neigh objective value of the best neighboring solution
Obj* objective value of the incumbent
The main parts of the algorithm are the following:
Initial solution generation:
As mentioned above, the initial solutions are generated only for the offered price and location of CICs and RFs. These centers are randomly assigned to the available places after reviewing and examining allocation requirements. When the assigning condition, the setup cost is lower than the remaining budget, is satis fied, these places are removed from the list of the available places and the associated variables X and Y, as well as the dependent variables φ and ξ, take 1. Besides, the capacity of centers is randomly selected from among the intended cases for sites (n) in the same way.
Now the amount of B_{pqt} is estimated by a random value between zero and the highest price calculated among customers for the used product p with quality level q in period t.
Generate the neighbor solutions:
At this stage, the entire possible neighborhoods must be investigated. In our algorithm, initially, the location is assigned to the CICs or the RFs without taking the feasibility into account (the procedure is illustrated in Algorithm 1 in Appendix). In Figure 5, a visual representation of generating neighbor solutions is depicted, where rows and columns show potential centers and periods, respectively and black dots indicate construction of the centers.
Afterward their feasibility is considered by functions given in Algorithm 2, and according to what was said about generating initial solutions, the value of offered prices (B) would be determined for every feasible answer. Then, a simplex search will be performed on the model using CPLEX library.
Tabu list:
In any iteration that has achieved a better answer than the existing answers, the current answer is transferred to tabu list to avoid the loop and duplication of answers. This way, the tabu list is specifically designed for objective functions, X and Y. After generating a neighboring solution and determining the objective value, the answer is primarily compared with the current solution. Afterwards, its nonattendance in the tabu list is checked (Algorithm 3). The comparisons are made using the objective values and if they are equal, they will be conducted through the comparison of X and Y.
Releasing condition:
In one of the iterations, if the best solution obtained in the neighborhood is better than the current solution, but has been in the tabu list, it is released from the list and instead, the last solution is listed.
According to the components explained in detail above, the structure of tabu search algorithm is given as Algorithm 4.
Algorithm 4. Main structure of the tabu search algorithm
Begin
//TL: Tabu List, contains Collection and Recovery Centers
and Fitness;
CapLevCol, CapLevRec: indicate Capacity Level of each
established Centers;
Cur, New, N1 and B: indicate Current, New, Best of
Neighbors and Best Solution, respectively//
Initialize Solution:
Cur_φ; Cur_ξ;
Cur_X; Cur_Y;
Cur_Q1, Cur_Q2, Cur_Fitness ← Simplex(Cur_ φ,
Cur_ ξ);
While iteration < Max_iteration do
N1Fitness = NB (Cur_ φ, Cur_ ξ);
If N1_Fitness < Cur_Fitness then
Tabu_List ← Cur_ φ, Cur_ ξ, Cur_Fitness;
Cur_Solution ← N1_Solution;
If TL_Checker (Cur_X, Cur_Y, Cur_Fitness) then
//Aspiration
Remove Cur_Solution from TL
End If.
If N1_Fitness < Best_Fitness then
Best_Solution ← N1_Solution;
End If.
Else
If TL_Checker (N1_X, N1_Y, N1_Fitness) then
Cur_Solution ← N1_Solution;
End if.
End If.
iteration = iteration+1;
End While.
End.
5.NUMERICAL EXAMPLE
To evaluate the performance of the proposed approach on small and mediumsized problems, some illustrative examples have been developed. The size of the investigated examples and their parameters value are specified in Tables 3 and 4, respectively. In Table 4, parameters have been randomly generated using uniform distributions.
Table 5 gives a report of the results obtained by the heuristic method coded in C++ as a subalgorithm, and numerous reruns at some iteration in the first instance. The resultant behavior of the proposed model is illustrated by a group of titles defined in the following:

Is calculated as the sum of the total quantity of uncollected returns from customers,

Is the remaining budget for making CICs and RFs at the final date of planned periods,

Total income that has been earned by the sale of new products,

Total cost is the sum of all costs that are spent in every period of the considered planning horizon,

Average time taken by the algorithm for solving the model,
According to the analysis conducted on different iterations, as shown in Table 5, it is observed that there was considerable amelioration among the obtained results although there was not a significant improvement at iterations higher than 200. Figure 6 depicts the objective values that are calculated from the difference between the total costs and incomes.
After this, the approach accuracy assessment results have been summarized using two solving methods, Lingo 13.0 solver and presented tabu search algorithm, in Table 6. All experiments were run in an Intel Core i3 CPU, at 2.13 GHz and with 4.00 GB of RAM. Moreover, the solution of Lingo is a local optimal solution, which can be known from the results in Table 6. To make a comparison, we have used the initial dimensions expressed in Table 3 and have randomly generated ten examples with different values of parameters in the range given in Table 4. According to this table and what is shown in Figure 7, it can be observed that in addition to reducing the time spent on solving the problems through tabu search significantly, the improving results obtained by tabu search for different examples is approximately 10% to 70% better than the Lingo results.
The results of the other instances are presented in Table 7 to show the effectiveness of the proposed tabu search approach in solving higherdimensional examples. All the runs were done at 200 iterations in the same manner as before. As it can be seen, the suggested approach is able to solve the problems that could not be solved by the commercial software in a finite time, though the solving time has dramatically surged with the increasing size of the problems.
As the results in Table 7 show, it can be understood which heuristic approach causes a significant improvement not only in the average time spent by the algorithm for solving the model but also in the objective values on medium to highsized problems.
6.CONCLUSION
To sum up, since there are not many articles that discuss financial incentive issues in reverse logistics, this paper investigated a dynamic mathematical model to appraise the incentive effect on return quantity of used products in a reverse logistics network design, which included different stages. In the optimization model, the RTD applying the fuzzy technique was used to interpret this relation. Furthermore, other effective parameters have been applied to configure a mixed integer nonlinear programming model for locating facilities and allocating customers to them, determining the optimal suggested price and planning the recovery strategy in a product return network.
Moreover, a heuristic method based on the tabu search procedure has been proposed for problem solving, which has benefited from the linearization and the exact solution approach. This was followed by comparing the results with those of Lingo commercial software for giving illustrative examples to analyze and validate the method. The computational results show the efficiency and effectiveness of the developed heuristic solution method when time complexity is addressed.
This study has developed several points that can be further investigated in future researches. For instance, the accuracy and efficiency of the proposed method could be improved. Furthermore, a number of verification and validation methods could be helpful in testing the accuracy and consistency of the process. Another potential extension to the setting investigated in this paper may consider the inclusion of different patterns for the behavior responding to the incentives offered. Finally, considering other stages of the reverse logistics network and improving the efficiency of the proposed method may be beneficial in testing the process.