A HYBRID NSGA-II ALGORITHM FOR THE CLOSED-LOOP SUPPLY CHAIN NETWORK DESIGN IN E-COMMERCE

Abstract. Designing the supply chain network is one of the significant areas in e-commerce business management. This concept plays a crucial role in e-commerce systems. For example, location-inventorypricing-routing of an e-commerce supply chain is considered a crucial issue in this field. This field established many severe challenges in the modern world, like maintaining the supply chain for returned items, preserving customers’ trust and satisfaction, and developing an applicable supply chain with cost considerations. The research proposes a multi-objective mixed integer nonlinear programming model to design a closed-loop supply chain network based on the e-commerce context. The proposed model incorporates two objectives that optimize the business’s total profits and the customers’ satisfaction. Then, numerous numerical examples are generated and solved using the epsilon constraint method in GAMS optimization software. The validation of the given model has been tested for the large problems via a hybrid two-level non-dominated sort genetic algorithm. Finally, some sensitivity analysis has been performed to provide some managerial insights.

-How product nature and its elasticity can affect our decision-making in supply chain management?
This research aims to study the closed-loop supply chain structure within E-commerce industries to answer the abovementioned questions. Therefore, a multi-echelon, multi-product closed-loop supply chain for online e-commerce with mixed quality defects returns is designed to optimize the supply chain profit while maximizing customer satisfaction simultaneously. Fuzzy set theory is being used to confront the uncertain nature of the parameters. Additionally, for encountering the high complexity of the proposed model, a metaheuristic algorithm has been employed.
The next section of this paper will examine the literature review on closed-loop supply chain and e-commerce supply chain. The third section is concerned with the methodology and the mathematical model that is used for this study. It also tries to refine the proposed model by eliminating nonlinearity and adding fuzzy set theory to confront uncertainty. The next chapter will concentrate on the algorithms and methods for solving the problem. Chapters five and six give sensitivity analysis and managerial insights regarding the proposed model besides analyzing the proposed algorithm's results. Finally, chapter seven gives the conclusion and future research directions.

Literature review
Several systematic reviews of closed-loop supply chains have been undertaken. One of the recent review papers that contain related articles to reverse logistics and closed-loop supply chain from 2007 till 2013 has divided 382 papers into four classes which are reverse logistics, closed-loop supply chain, sustainable supply chain, and green supply chain [14]. Based on this review paper, reverse logistics and CLSC papers mostly focus on recycling, remanufacturing, repairing, and disposing of products to have a sustainable supply chain. However, few papers focus more on the business side of the CLSC and the effect on customer retention and customer acquisition than environmental factors. Both of these groups have been mentioned in this paper, respectively.
Pishvaee et al. [25] integrated the forward and reverse supply chain's design to avoid the sub-optimal answers. Consequently, they studied an integrated forward-reverse logistic network design in a multi-objective, multi-stage model considering multiple capacity levels. After formulating a mixed-integer nonlinear programming model that tries to maximize accountability and minimize costs simultaneously, a dynamic search-based memetic algorithm has been proposed to solve the model in large-size problems. Pishvaee et al. [26] developed mixed-integer linear programming (MILP) for closed-loop supply chain networks that use scenario-based stochastic programming one year later. Since considering parameters as deterministic in strategic planning might cause high costs, robust optimization is being used in this paper despite being computationally expensive. Ramezani et al. [27] took another similar approach two years later. They also used robust optimization to cope with demand and return rate uncertainty in a multi-product, multi-echelon CLSC. For solving the large size problems, they used a scenario relaxation algorithm as only a small subset of scenarios seems applicable and can be used to search for the optimum solution.
Ozceylan et al. [24] integrated strategic and tactical decision making in a CLSC network design problem that intends to minimize the system's total costs using mixed-integer non-linear programming (MINLP) mathematical model. In a gold industry case study, Zohal and Soleimani [39] examined a 7-layer CLSC network (4 echelons in forward and three echelons in reverse) with special attention on CO 2 emissions. They sought to design a network using the ILP model with two objectives of minimizing costs and emissions. In the end, an ant colony-based algorithm was used to solve the problem. Another CLSC model was designed by Forouzanfar et al. [10], which used a multi-objective MINLP model. They used GAMS software to validate and solve the model for small-size problems. The time-consuming process and solving the model with exact algorithms and high complexity of the proposed NP-hard problem led them to develop a metaheuristic algorithm. To meet this end, they proposed a bee colony optimization algorithm and compared its results with the genetic algorithm.
Unlike previous research in reverse logistics, Zhalechian et al. [36] surveyed the design of a sustainable CLSC with its effect on the economy, environment, and society. In this regard, the impact of the designed supply chain on CO 2 emissions, fuel consumption, wasting energy, and job position creation has been considered.
A three-level CLSC with mixed uncertainty, which is modeled by stochastic-possibilistic programming, is being solved using a hybrid metaheuristic algorithm. Iassinovskaia et al. [17] studied closed-loop inventory routing problems for returnable transport items with simultaneous pickup and delivery. A mixed-integer linear program is proposed for the pickup and delivery inventory-routing problem within time windows over one planning horizon in closed-loop supply chains. After testing this MIP programming formulation on small-scale instances, a cluster first-route second metaheuristic is also proposed to reach satisfactory solutions for real problems. Despite the fact that the beneficial aspects that reverse logistics can bring to the environment, designing and maintaining a closed-loop supply chain is costly for any corporation. One study considered the government as an accelerator for encouraging recycling business operations for corporations by offering incentives, fees, and tax exemptions. The difference between this article and other articles was on considering the government's role in the economy. The abovementioned papers discussed the problem in a government-independent context, while here, a government-dependent economy is studied for designing a sustainable supply chain, a Two-echelon reverse supply chain with one manufacturer and one retailer developed, which tries to increase customer willingness to return EOL products by offering discounts and fees. Also, to consider the government's role in the supply chain, some scenarios have been applied on the extended model, investigating the effect of tax exemption and subsidy on both retailer and manufacturer [15].
In a study investigating closed-loop supply chains, Tavakkoli-Moghaddam et al. [33] focused on a locationinventory design of a multiple spare parts' supply chain, which minimizes the total costs of inventory and location-allocation decisions. There exist two types of flow from distribution centers to operational bases (forward demand satisfaction flows) and operational bases to repair centers (backward flow of spare parts). Because of the complicated nature of the proposed MINLP model and consideration of demand as a stochastic parameter with normal distribution, after solving small-size instants of it with GAMS software, they moved to use evolutionary algorithms, specifically particle swarm optimization (PSO) algorithm to cope with the complicated np-hard problem and reach to solutions in reasonable times. Finally, some sensitivity analysis was held for evaluation of the proposed model. A considerable amount of literature has been published on CLSC considering uncertainty since not seeing the unstable and non-deterministic nature of parameters may result in inadequate planning and consequently more costs. The longer the planning horizon is, the more uncertainties we have to tackle to plan and design the supply chain. Hence, ignoring uncertainty will incur dramatic costs in tactical and strategic decisions in supply chain management. Ghahremani-Nahr et al. [11] studied the location and allocation problem with uncertainty in customer demand, return ratio, transport costs, purchase price, and shortage costs in a multiechelon multi-product, and multi-period closed-loop supply chain. They also considered shortage and discount in their proposed MINLP model. They unravel the uncertainty part by using robust fuzzy programming. Moreover, they used a whale optimization algorithm (WOA) to address this problem for large-scale problems.
One particular paper focused on the recycling approach from sustainable development methods such as environmental management, green purchasing, reusing, and remanufacturing. This research aimed to investigate pricing and ordering decisions in a three-echelon (consisting of a collector, recycler, and a manufacturer) dual-channel supply chain for waste recycling processes by using game theory. In the reverse supply chain, a manufacturer capable of recycling waste can receive the EOL products through two channels: collector and recycler. As the names suggest itself, the collector only collects the unrecycled waste, whereas the recycler recycles it. This paper's game theory perspective arises when the manufacturer has to decide between buying from the recycler or buying from the collector and recycling itself [18]. Giri and Dey [13] posed what the manufacturer would do if the collector fails to collect the required amount of waste. Thus, they entered the third party to the Jafari, Hejazi, and Rasti-Barzoki's model by introducing a backup supplier to cover shortage from the collector side. They considered the model as a Stackelberg game in which the manufacturer acts as the leader and other interactions among other parties do not affect the manufacturer's pricing strategies. Mohtashami et al. [23] focused on environmental impacts and energy consumptions within a green closed-loop supply chain. They selected queueing theory to optimize the logistics and waiting time of transportation fleets' network.
Based on Savaskan et al. [29], a CLSC comprises three key components: manufacturer, retailer, and a third party. Each component would behave and make decisions in a cooperative or non-cooperative CLSC environment, usually by game theory. Alamdar et al. [2] considered a multi-level fuzzy CLSC under various cooperative and non-cooperative scenarios. The reason for using fuzzy set theory in their model was considering uncertain demand-related and controllable by selling price and selling effort. The three scenarios focusing on game theory were cooperation between the manufacturer and collector, a coalition between retailer and manufacturer, and cooperation between retailer and collector. A fuzzy multi-objective mixed linear integer programming has been proposed by Fakhrzad and Goodarzian [9] for a production-distribution model in a multi-product, multi-level green closed-loop supply chain network. The objectives of this study were to minimize total cost, minimize gas emissions and maximize reliability.
The last decade has seen numerous studies focusing on the closed-loop supply chain. Accordingly, the uncertainty conditions and the green emissions of facilities are still open issues. In this paper, a new fuzzy multi-objective programming approach presents a production-distribution model to develop a multi-product, multi-period, and multi-level green closed-loop supply chain network problem, which is formulated as multiobjective mixed linear integer programming (MOMILP). In regards to the offered fuzzy multi-objective model, three conflicting goals are exited simultaneously. The objective functions are to minimize the total cost, minimize the gas emissions costs due to vehicle movements between centers, and maximize delivery demand reliability due to the suppliers' reliability. Fuzzy numbers consider the parameters of the model for real-world applications. Another novelty of the proposed model is in the solution methodology. This study provided a new modification of Imperialist Competitive Algorithm (MICA) algorithm in order to address the proposed model. Additionally, this algorithm has been compared with other common metaheuristics such as ant colony optimizer (ACO), genetic algorithm (GA), particle swarm optimizer (PSO) and the well-known Imperialist Competitive Algorithm (ICA) in order to demonstrate the efficiency of the proposed method. Finally, different analyses with a variety of problem complexity in different sizes are performed to assess the performance of algorithms as well as some sensitivity analyses on the efficiency of the model are studied. Finding the best locations for depots, allocation of customers, and the optimal routes to minimize both total cost and the overall violation from the customers' defined time limits was the problem that Mamaghani and Davari [21] focused on. Therefore, they proposed a bi-objective periodic location routing problem with simultaneous pickup and delivery with time windows. Afterward, they used a non-dominated genetic algorithm (NSGA-II) and a non-dominated ranked genetic algorithm (NRGA) to come up with a set of solutions for the proposed problem.
Unlike the abovementioned studies in the closed-loop supply chain, which considered return as a tool for designing a sustainable system, the below articles investigated returns to improve customer satisfaction in online shopping businesses. In 2102, Liu and Xu designed an integrated multiperiod and multicriteria decision-making model for a supply chain network in an e-commerce platform that contains both B2B and B2C transactions. Their proposed model consists of three tiers (manufacturers, retailers, and consumers) who compete and cooperate to reach the network equilibrium. A study by Liu et al. [20] involved a closed-loop location-inventory-routing e-commerce supply chain design. In the literature, usually two of the three sides of the decision-making were studied, which were location inventory problem (LIP), location routing problem (LRP), and inventory routing problem (IRP). Because of the interdependence between facility location, inventory control, and transportation optimization, all these three factors were considered in a three-level optimization LIRP (location inventory routing problem) model.
In this paper, a stochastic LIRP model is solved with a pseudo-parallel genetic algorithm integrating simulated annealing. In 2016, Deng, Li, Guo and Liu published a more sophisticated paper.
Return products under an e-commerce environment can be categorized into non-defect returns (like sizing issues for clothes) and quality defects (quality mismatch with customer requirement). Therefore, the model should take mixed quality defect returns (MQDR) into account. A three-echelon Supply Chain considering decision-making in facility location, inventory, and transportation with both quality defect and non-defect returns has been presented in this article. What is used to address the proposed np-hard problem was the hybrid ant colony optimization algorithm and ABC algorithm. Zhang et al. [37] concentrated mainly on routing and transportation problems in e-commerce logistics trading systems. The rise of online shopping sales from a positive perspective is generating more money and more marginality for e-commerce stakeholders. On a more negative note, the operation departments and logistic services face difficulties in providing massive product delivery services. This nightmare situation, which arises from the booming demand for online shopping, is not the only problem for delivery systems. The increase in fragment transport demand, empty backhauls, idle capacities on lots, and rising fuel prices are more problems that are causing more operational-related costs for an e-commerce business. Thus, online retailers' systematic view of the supply chain needs to be taken so that resource allocation optimized and logistic services quality is enhanced. A carrier collaboration network in the e-commerce logistic system with multiple less-than-truckload carriers and multiple vehicles is considered, aiming to minimize transportation costs. Lastly, a new meta-heuristic approach stochastic plant-pollinated algorithm is developed to solve the MILP model.
Various methods, both exact and metaheuristics, have been proposed to solve multi-objective problems (MOP) on a more solution-oriented note. Some of the most common exact MOP methods are the weighted sum method, lexicographic method, goal programming methods, and Epsilon constraint method [22]. The weighted sum method tries to convert the MOP problem into a single objective problem (SOP) by assigning a weight to each objective function and then calculate the weighted sum of all the objective functions as the new objective function. In the lexicographic method, objective functions are being optimized in order of their importance. In goal programming methods, a goal should be set for each objective function. Next, the new objective function will be defined as the sum of deviation of all objective functions from their goal. Compromise programming is also an exact MOP method similar to goal programming methods. Their difference between these two models is that in goal programming, the goal are being set by the decision-maker while we use the best SOP solution as the goal for the objective function in the compromise programming method [30]. Another useful method for MOP is epsilon constraint in which an objective function would be selected as the main objective function to use for SOP, and the other objective functions will be classified based on an epsilon value. More about this method and the modified version of it, augmented epsilon constraint, will be explained in section four.
Besides exact MOP methods, various multi-objective evolutionary algorithms (MOEA) have been proposed like non-dominated genetic algorithm (NSGA) and multi-objective particle swarm optimizer (MOPSO) [38]. Non-dominated genetic algorithm (NSGA-I), which was introduced by Srinivas and Deb [31] as a MOEA method, ranks the population based on nondomination. Deb et al. [6] proposed the improved version of NSGA-I, named NSGA-II, which the crowding distance concept was added to the previous version [7]. The third version of NSGA had some changes in the selection mechanism and was more useable for MOPs with more than two objective functions [5]. Various improved NSGA-III have been proposed after that to cover many-objective optimization problems [35] and big data optimization problems [34].
A more detailed classification based on seven factors, including the type of uncertainty, single and multiperiod models, single and multi-product models, single and multi-objective models, and the method used to solve the problem is illustrated in Table 1. As shown, most studies only focused on the total supply chain's cost, and few studies have addressed customer satisfaction in their supply chain design. Due to the new era's competitive economy and business, all companies are customer-oriented, and unsatisfying a customer's need will cause a considerable amount of costs in the long run. Although some studies have considered unsatisfied demand as shortage costs besides all other operational costs, only Mamaghani and Davari [21] used a separated objective function to cover customer satisfaction (by defining delivery due date as a customer satisfaction factor). All summed up, it is essential to consider customer satisfaction in our study. We have addressed customer satisfaction by three factors which are 1) Customers' demand coverage 2) Delivery due date and 3) Discount on products as a customer acquisition and retention factor. For the first time, the present research explores the effects of products' discount as a customer satisfaction factor on designing the CLSC network. Additionally, this research addresses the uncertainty within the ecommerce businesses using fuzzy set theory and is student the effect of pricing decision in e-commerce. Finally, few studies used metaheuristics algorithms to solve real-world size problems. This study provides a multi-objective, multi-product, and multi-period closed-loop supply chain design in the e-commerce business. The model tries to answer to the following decisions: (1) What locations should be chosen for fulfillment centers, delivery centers and after-sales departments in order to minimize establishment costs and cover most customers? (2) How much items should be ordered to minimize inventory costs and fulfill maximum customers' demand? (3) How much discount can be used to maximize our profit and keep the customers motivated and satisfied for buying the products? (4) How our transportation and routing should be handled to optimize transportation costs while delivering products within the customer's expected time interval?
Regarding previous works in closed-loop supply chain design, there are several important areas where this study makes original contributions. The present research explores the effects of products' price and discount on customer demands and revenue management. Unlike most of the research focusing more on cost minimization in the supply chain, this study provides new customer satisfaction insights and proposed a bi-objective mathematical model with both cost and customer-oriented objective functions. This paper studies, the impact of uncertain decision-making parameters on the e-commerce supply chain and provides many managerial insights based on analyzing the model. Last but not least, this study proposed a two-stage hybrid metaheuristic method that is applicable for solving big-size problems.

Problem description
In this study, a multi-product, multi-period, multi-echelon closed-loop supply chain network in E-commerce is designed under various capacity levels and uncertainty. This research was inspired by the case of the retail business model in e-commences. Despite the marketplace business model, which uses the seller's resources and items for servicing customers, customers' demand should be forecasted in the retail business model. Based on the company's budget and warehouse capacity, all purchases are taking place. Figure 1 illustrates the forward and reverses logistics in the proposed closed-loop supply chain.
In the proposed location-inventory-pricing-routing problem, four questions should be answered after solving the mathematical model.
(1) Location: what locations will be chosen for establishing fulfillment centers, distribution centers, and aftersale departments? Binary decision variables are defined to answer whether we are constructing these facilities or not, which will detail the mathematical model formulation. (2) Inventory: How many of each product is our demand? How many products should we order from each supplier for the replenishment part, considering our demand, inventory capacity, and ordering costs? (3) Routing: From which distribution center customers get their order, and in what order products are delivered to customers are in which the demand is covered, the desired time period of the customer is satisfied, and transportation costs are minimized? (4) Pricing: what would be the desired and optimum price for each product, which can increase the margin while keeping customers satisfied by proposing considerable discounts?
The proposed mathematical formula aims the answer the abovementioned questions. In the forward flow network, ordered products ship from suppliers to fulfillment centers. Fulfillment centers receive inventory and store products for current or future customer's needs. Then, demanded products will be sent to distribution centers, where all products deliver to customer zones. In the reverse flow (which is noticeable with red arrows in Fig. 1), a proportion of the product will be sent back to the distribution centers and then to the fulfillment centers. If the fulfillment center has an after-sales department or not, two things can happen. If the fulfillment center cannot handle returned products (due to the lack of an after-sale department), all returned products would be destroyed. On a more positive note, if the fulfillment center is equipped with an after-sales department, it can use a portion of the returned products as a new consignment for future demands. Therefore, this model will show the tradeoff between the costs of adding after-sales departments and money savings for using returned products. Throughout this paper, the term elasticity coefficient will refer to the correlation between a product's demand and the amount of discount on that product. The proposed model will maximize the total profit of the supply chain and improve customer satisfaction by minimizing unmet demand, minimizing delayed orders, and maximizing returned and acquired customers caused by proposed discounts.

Assumptions
Following assumptions are made for formulating this problem: -Suppliers and fulfillment centers have limited capacities.
-The location of suppliers and customer zones are known.
-There are a set of potential sites for locating fulfillment centers with different capacity levels.
-The after-sale department can only be opened in one location where a fulfillment center exists.
-There are a set of potential sites for locating distribution centers.
-The location of all nodes has to be determined.
-Different types of products are considered in the CLSP.
-The unsatisfied demand of customers is back-ordered.
-Customer demand and the fraction of returned products are considered as uncertain parameters.
-Fulfillment centers with the after-sales department can handle mixed quality defect returns.

Notations
Sets Index Definition S Set of suppliers (∀s ∈ S) F Set of fulfillment centers (∀f ∈ F ) N Set of capacity levels available for fulfillment centers (∀n ∈ N ) D Set of distribution center (∀d ∈ D) C Set of customer zone (∀c ∈ C) V Aggregate index of customer zones and distribution centers ( Quantity of returned product i shipped from distribution center d to fulfillment center f in period t R c idct Quantity of returned product i shipped from customer zone c to distribution center d in period t UDict Quantity of non-satisfied demand of product i of customer zone c in period t 1 if vehicle r is assigned to distribution center d in period t yruvt 1 if there is a route of vehicle r from node u to node v in period t atvt Arrival time vehicle for delivering orders to node v in period t utct 1 if orders does not deliver to customer zone c on time in period t

Mathematical formulation
The first objective function (3.1) tends to maximize the total supply chain revenue, which can be classified under blow terms: -Total income (first term).
-Costs of constructing fulfillment centers and distribution centers (second and third terms) and equipping fulfillment centers with after sale department (4th term). -Costs of opening distribution centers (5th term).
-Logistic costs from supplier to fulfillment center (8th term) and from the fulfillment center to distribution center (9th and 10th term). -Collection and distribution costs in distribution centers (11th and 12th term). -Vehicles' transportation costs for delivering products to customers (13th term).
The second objective function (3.2) tries to improve customer satisfaction by minimizing unsatisfied demands and the number of delayed or early-shipped orders while maximizing the customer acquisition and retention gained by discounting. Constraint (3.3) sends the total orders from suppliers to the fulfillment centers. Constraint (3.4) is the capacity constraint associated with suppliers. Constraint (3.5) ensures that the selected fulfillment center is established before it receives orders from suppliers. Constraint (3.6) states that one capacity level must be selected for a fulfillment center at most. Constraint (3.7) implies the minimum required inventory threshold for establishing fulfillment center f with capacity level n. Constraint (3.8) guarantees that the inventory level of a fulfillment center does not exceed its capacity. Constraint (3.9) controls the input and output volumes of the fulfillment center. The constraint in inequality (3.10) ensures that products flow to a distribution center only if it is open at that period. Constraint (3.11) checks that a distribution center can be opened if it is established already. Equation (3.12) is a balance constraint on the distribution center and ensures that the product's input flow equals the output flow. Constraint (3.13) guarantees that customer demand must be met until the last moment. It is worth mentioning that the demand of product i in customer zone c in time t is influenced by three factors, which are 1-base demand of the customer zone for product i, 2-discount level for the product i and 3-elasticity coefficient of that product for each discount level. Constraint (3.14) presents the number of returned products a period after being bought by customers. Constraint (3.15), like (3.12), is a balance constraint and assures that input and output flow of return products remain equal. Constraint (3.16) implies that an after-sale department can be opened when a fulfillment center has been opened at the same site. Constraint (3.17) shows that returned non-defect products will come back to the supply chain again and are sent to the fulfillment center as a consignment for future sales if the mentioned fulfillment center has an after-sales department. Constraint (3.18) states that no vehicle can start their trip from a distribution center unless established and has flows from or to customer centers. Constraint (3.19) makes sure that only one distribution center is in each route. Constraint (3.20) refers to the capacity of distribution centers on using the number of vehicles. Constraint (3.21) and (3.22) indicate that a route can be made from the distribution center to the customer zone and vice versa if there is a product flow between them. Constraint (3.23) also says that a vehicle can travel between two customer zones when they are on the same route and receive their demand from the same distribution center. Equation (3.24) limits the customer zones to get their demand from at least one distribution center. Constraint (3.25) ensures that each customer zone should be a route if a product flow toward them. Constraint (3.26) is related to the sub-tour elimination. Equation (3.27) is related to flow conservation and guarantees that the route remains circular and that when a vehicle visits a customer center, it should leave the visited node. Constraint (3.28) makes sure that only one distribution center is in each route. Constraint (3.29) and (3.30) calculates the arrival time of vehicles for each node. Constraint (3.31) and (3.32) also examine whether a node receives its orders on time or not. Constraint (3.33) states that each product can have only one discount level in each period. Constraint (3.34) states that the beginning of the route is the starting point of time. Constraint (3.35)-(3.37) set the initial value for inventory, unsatisfied demand, and product flow toward customers, respectively. Finally, constraints (3.38) and (3.39) define the type of decision variables.

Linearization and fuzzy programming
The abovementioned mathematical model has two drawbacks, which are: (1) It is an MINLP model that is pretty hard to solve for large instances.
(2) All parameters are considered deterministic, even parameters like the ratio of return and demand, which are uncertain.
Therefore, the proposed model will be modified using linearization methods and fuzzy set theory in this section [1].
Considering constraint (3.13), regarding the fact that demand is by nature an uncertain parameter; we use fuzzy set theory and chance constraint programming (CCP) to modify the constraints. Suppose that dem ict is a trapezoidal fuzzy number ( dem ict1 , dem ict2 , dem ict3 , dem ict4 ) and we want the necessity of this constraint to be more than ϑ d .
Consequently, the equivalent crisp parametric constraint can be written as follows: The same process should work on constraint (3.14) too.
Knowing that uncertainty parameter of the constraint has to be formed with a satisfaction level of at least ϑ α the crisp format of the above fuzzy constraint would be as follows: Not only constraint (3.17) suffers from uncertain parameters, but also it is a nonlinear equation. Therefore, first, we linearize it and then transform it to a crisp constraint. The linearized version of the equation (3.17) is as follows: After applying CCP we would have: For transferring constraints (3.21) and (3.22) into a linear format, we introduce an auxiliary decision variable named W dct to cope with nonlinearity problems. Therefore, instead of equation (3.21)  Lastly, there are nonlinear expressions and uncertain parameters in the first objective function (3.1) in the 1st and 5th term. For the first term, which is concerned with total revenue, we introduce another auxiliary variable named G iqct to linearize our equations.
Also, for the 5th term of the objective function (3.1), we add three below constraints to the model.
In addition, we use the expected value of the objective function instead for fuzzy programming. Therefore, the first objective function would be like this:

Solution methods
The proposed model in this study is a mixed-integer nonlinear mathematical model with uncertain parameters. While we took care of nonlinearity and uncertainty in the previous section and turned the MINLP model into a MILP model, solving the transformed model with exact methods is still a computationally expensive task. Because the remained problem is an Np-hard problem, using heuristic and metaheuristic algorithms seems to be a reasonable alternative for large instances in this sense. Therefore, after defining three numerical examples with small, medium, and large problem size, we used augmented epsilon constraint as an exact method to solve the small and medium-sized problems and adopted a hybrid NSGA II algorithm as a metaheuristics method to resolve the large-sized instance. GAMS software is used for the augmented epsilon algorithm, and for the hybrid NSGA II, both GAMS and MATLAB have been used. The following will explain both methods in more detail.

Augmented epsilon constraint method
The selected multi-objective decision making (MODM) method for solving the problem and finding exact solutions was augmented epsilon constraint [4]. One of the posterior methods is applied to a multi-objective integer programming problem and can generate the exact Pareto set if the step size is appropriately chosen. Initially, we solve the problem for each objective function to extract the ideal and anti-ideal points. If one of the single objective problems has multiple answers with the same objective function, the solution that optimizes the other objective function will be chosen (lexicographic optimization method). After generating the Payoff table, we divided the interval of ideal and anti-ideal points for each objective function into several equal segments and assigned segment edges as epsilon for the below SOPs. It should be noted that the decision-maker selects one of the objective functions based on their priorities.
Where e is substantially small and r k is the range of kth objective function.

Two-level non-dominated sorting genetic algorithm II
In this section, we explain the metaheuristic method used for solving the model in large-sized instances. The used metaheuristic for solving this problem is the non-dominated sorting genetic algorithm 2 (NSGA II) developed by Deb et al. [7]. However, the used chromosomes in this problem do not contain all the mathematical model variables (due to colossal solution space). We divided the main problem into two sub-problems so that the integer and binary variables fall into the first sub-problem and the rest of the variables, which are positive continuous variables, remain in the second sub-problem. All variables in the first sub-problem will be presented in the chromosome and chosen randomly in the algorithm. After choosing one set of random variables for the first sub-problem, the second sub-problems will be solved using linear programming. In the second sub-problem, all the integer and binary variables will be considered as model parameters. As all the remaining variables are continuous and there is no nonlinear expression in the model, the second sub-problem can be solved using linear optimization methods. If we consider n as the size of all integer variables in our problem, then the first sub-problem's complexity would be O (nlgn) [12]. The complexity of the second sub-problem, which is being solved using linear optimization methods, is not related to the input variables' size but is related to the number of constraints. Consequently, the proposed procedure can solve the problem with polynomial complexity. In the following section, we will describe how the algorithm is being used to solve the proposed problem. The process of the abovementioned two-level NSGA II algorithm can be seen in the following pseudocode: Procedure Two level NSGA-II(n,g,f(x)) -Initialize Population P -Generate Random Population -size n -Extract decision variables from chromosomes -Call GAMS to solve the LP model

End Procedure
Used chromosomes for representing binary variables have six parts: (1) A 2 × f shape matrix which displays fulfillment related variables ( Figure 2).
The index of the matrix is related to the index of fulfillment center (f ). The number in the first row represents what capacity level is being used (U f n ). The second row also depicts if the fulfillment center is equipped with an after-sales department (U f ).
(2) A t × d shape matrix which includes distribution center variables (Figure 3).
The columns' index represents the distribution center (d), and the index of rows illustrates periods (t). The numbers of the matrix determines whether the distribution center is open or not in the period (Z dt ). If one column have only zero values, then the distribution center will not select for constructions at all (Z d ).
(3) A t × r shape matrix which demonstrates assignment of vehicles to distribution centers (Figure 4).
The columns' index represents the vehicle index (r), and the index of rows illustrates periods (t). Numbers in the matrix show that vehicle r is assigned to which distribution center (d) at time t. If the number is 0 then its means that that vehicle is not assigned to any distribution center. This matrix is a representation for variable Y rdt (4) A 2 × c × t shape matrix involves information about the routing problem in the main problem ( Figure 5).  The indexes of this matrix represent customer zones (c) and periods (t). Numbers in the first row indicate what vehicle is assigned to a customer zone c at time t. numbers in the second row also represent the sequence of customer zones that vehicles visit. If we sort the second-row numbers in the ascending order for each vehicle, it shows the sequence of customer zones that each vehicle visits (y ruvt ).
(5) A i × t shape matrix which depicts the discount level of product i in period t ( Figure 6).
(6) A single continuous number between 0 and 1 which is the weight of the first objective function and takes care of solving the second sub-problem using the weighted sum method.
After generating the mentioned chromosome (the first population, cross over, or mutation), we would solve the second sub-problem with GMAS software and linear programming solvers. Note that all binary and integer variables are now model parameters in the second sub-problem. Additionally, the second sub-problem will be solved using the weighted sum method with the weight assigned from the first sub-problem in chromosome generation.
Genetic operators are used in genetic algorithms to generate diversity and to combine existing solutions into others. The main difference among them is that the former operates on one chromosome, that is, they are  unary, while the latter are binary operators. To solve the described problem, we have used uniform crossover for continuous numbers between 0 and 1 (second row of 4th part of the chromosome) and a k-point crossover for integer numbers.
Additionally, Scramble and flip bit mutation has been used for mutation of continuous and integer values.

Hypervolume metric (HM)
Finding optimum and diverse Pareto solutions are two objectives of every MODM method. Various metrics have been defined to evaluate these two objectives [3]. In this study, we have used a hypervolume indicator, also known as S metric, which measures the covered space of found Pareto frontier by the algorithm compared to a reference point to validate our proposed algorithm. For each Pareto frontier solution found by the algorithm, v i represents the v-dimensional Lebesgue measure by that solution.
Larger HV values mean better methods and better Pareto solutions.

Numerical examples
This section defines multiple instances of the problem with different sizes and tries to solve them with the last part's proposed algorithm. Before turning into the metaheuristic algorithm for solving the problem, we evaluate our proposed model with exact methods (augmented epsilon constraint) only for small and mediumsize instances. After ensuring the model's validation, we run the proposed algorithm for all instances, including No. of discount levels (Q)   1  3  1  1  2  5  3  3  2  3  2  6  2  3  3  10  7  4  5  6  3  8  3  5  5  20  13  7  6  10   Table 3. Source of randomly generated parameters.
Parameter Range Parameter Range a large-sized problem, which its size is close to standardized real-world problems. Table 2 shows the size and characteristics of each test problem. The algorithm was executed on a laptop, with an Intel(R) Core (TM) i7-8550U CPU @ 1.80 GHz 1.99 GHz processor and 8 GB RAM of memory. Table 3 depicts the values of defined parameters in test problems, which are generated randomly based on uniform distribution.
Parameter tuning in meta-heuristics is one of the pivotal factors that affect the performance of the algorithms. Taguchi experimental design is one of the best methods for determining the parameters of the algorithm. The most critical parameters in the NSGA II algorithm are the population size, iteration number, and crossover and mutation percentages. The results of Taguchi methods are presented in Table 4 and Figure 9.

Results
Initially, we solve each objective function's problem to extract the ideal and anti-ideal points using the lexicographic optimization method, as shown in Table 5.   Table 6. Decision variables for the test problem 1 (discount levels). The optimum values for each decision variable solving the first test problem with the first objective function can be seen in the tables below and Figure 10. Table 6 shows the optimum values for discount levels. Table 7 depicts the inventory side variables such as storage and ordering variables. Tables 8-10 demonstrate the routing and allocation side of the forward and backward logistics problem. Figure 10. The structure CLSC for the test problem 1. Table 8. Decision variables for the test problem 1 (routing and allocation variables in forward logistics). Table 9. Decision variables for the test problem 1 (routing and allocation variables in forward logistics). After generating the Payoff table, we solve the small and medium-size problems on both exact (augmented epsilon constraint) and metaheuristic (two-level NSGA II) algorithms. The results for each test problem are presented in Table 12. Figure 11 and Table 11 also compare the results of the exact methods and the proposed hybrid metaheuristic method for the first two test problems (Table 12).
As shown in Table 11 and Figure 11, the gap of the metaheuristic algorithm is reasonably small. Thus, it would be a suitable replacement for the exact method, especially for large-size problems. Therefore, down to the fact that using the exact method for a fuzzy nonlinear model would take a considerable amount of time, we will use the proposed metaheuristic to overcome this complexity. The obtained solutions from the proposed metaheuristic algorithm can be seen in Table 13. What is worth paying attention to is the direction of objective functions in Figure 12. Considering the first question of the research in the introduction part, it seems that Table 10. Decision variables for the test problem 1 (routing and allocation variables in backward logistics).   Figure 11. Comparison of exact method and the proposed method for small and medium-size instances.  these two conflicting objective functions cannot be optimized at the same time. Therefore, the decision-maker is to decide which solution suits better for the business and its preferences (Table 13). Table 14 depicts the HV indicator of both exact and metaheuristic algorithm and Figure 13 shows the HV metric over iterations in the metaheuristic algorithm. As it can be seen, the metaheuristic method provides an HV metric almost as good as the exact method.

Sensitivity analysis
Sensitivity analysis is a beneficial way to understand the model more promptly and grasp valuable insights from it. The behavior of two objective functions in various inventory costs and distribution centers' opening  costs can be seen in Figures 14 and 15. As it will be expected, by increasing both unit inventory holding cost of the product (h i ) and opening cost of distribution centers (od d ), the first objective function (total profit) will shrink down. However, the reason for the growth in the second objective function while boosting costs should be addressed. The higher unit inventory cost, the lower preference to store products, the less demand coverage, more uncovered demands, and lower customer satisfaction. The same thing can be said for the opening costs of distribution centers. With higher opening costs, the model tends to select fewer distribution centers for operating in the supply chain, and consequently, lower demands will be covered. Another mentionable fact is that changing the opening costs of distribution centers, which are closer to customer zones, has more impact on customer satisfaction objective function. It is to be noted that the below results were obtained using the weighted sum method of normalized objective functions.
Another set of interesting parameters that can be studied in sensitivity analysis would be the threshold and capacity of fulfillment centers (thr n and Cap f n ), illustrated in Figures 16 and 17. Having more fulfillment capacity leads to cover more demand, which brings more customer satisfaction (lower values for second objective function) and more sales and more profit (higher amount of first objective function). One interesting point is the severe changes of the first objective function on some capacity ranges and increases in the second objective   function (around 4,000 and 8,000). This behavior's motive is that in these values, the optimum solution tends to decrease the number of fulfillment centers as it can cover the demand with lower fulfillment centers, which causes lower construction costs. In these points, the model decides to cover less demand (higher values of the second objective) but have lower construction costs and higher profit. Turning into the minimum threshold for constructing fulfillment centers, the higher threshold leads to limitations on constructing fulfillment centers and limitations on covering demand, making a lower profit and lower customer satisfaction. It also should be pointed out that the below figures acquired optimizing only the first objective function. One pivotal note that is worth mentioning here is that by having a detailed look at Figure 14 throw Figure 17, it can be seen that both objective functions are moving in the same direction simultaneously. When the first objective function increases, the second function decreases, and both objective functions are moving toward a better solution. When the first objective function decreases, the second function increases, and both objective functions move to a worse point. Therefore, this same direction change of objective functions is being obtained from managing the supply chain's limitations, such as the capacity of fulfillment centers and their threshold, opening, and inventory costs.
In the introduction section, we introduced three central questions of this study. In the previous section, we witnessed the answer to the first question. Two objective functions of this model do not move towards optimism in one direction and have conflicting objectives. Therefore, instead of one single optimum answer, we have a set of non-dominated solutions. The decision-maker is to decide which answer should be chosen from the Second objective function First objective function Average of fulfillment centers' capacity  Pareto-optimal solutions. The next two questions will be answered in this section using sensitivity analysis. We investigate the correlation between uncertain parameters and objective functions. In Equation (3.41), we introduced the ϑ d parameter for the necessity of constraint (3.13). The left side of the constraint (3.13) should remain within an interval, which was calculated in Equations (3.42) and (3.43). The lower necessity we set for the demand parameter, the wider the interval range we have. This state can be seen in the following figure.
The below charts provide the objective function values of SOP models per each necessity amount for the constraint (3.13). What stands out from the right-side chart (Fig. 19a) is that when we give lower necessity to the model, the objective function value improves. The reason for this behavior is that while running the SOP model for the second objective function, constraints (3.42) and (3.43) are active and inactive, respectively, which means that the SOP model of the 2nd objective functions only consider the lower amounts of demand and that is the reason that the number of unsatisfied demands will shrink by reducing the necessity of constraint (3.13). On the other hand, the SOP model of the first objective function takes the higher amount of demand into consideration as by having more demand, the profit and revenue will be more. The reason for some fluctuations and revenue reductions (in 0.4 and 0.7 necessity) is that covering more demand needs more facilities (more fulfillment centers and more delivery centers) which will generate more construction costs. Figure 19 shows the behavior of both objective functions in the sensitivity analysis of constraint (3.13). The third question was related to the effect of product nature and its elasticity on supply chain management. The Figure 20 were obtained after changing the elasticity coefficient for all products and applying and running two SOP models for each objective function. The graph below shows that the more discount influences demand, the more profit and revenue we can have. However, customer satisfaction seems to act differently with higher elasticities. There is insightful information within the below charts, proposing discount as a powerful tool for controlling demand and consequently managing our supply chain. A more in-depth analysis of the usage of discounts for managing demand and customers will be discussed in the managerial insights section.

Discussions and managerial insights
The present study adds to the growing research of e-commerce supply chain management. This research, compared to the other studies, involves decision-making in uncertain environments. It also involves strategic, tactical, and operational decision-making in one model considering having location problem as long-term decision-making, inventory management as mid-term decision-making, and routing problem as short-term decision-making. It covered supply chain internal decision variables and contained some decision variables as leverages that could have been used to control the external factors (such as discount variables that are able to control demand, the external parameter of the supply chain). We witnessed from the right-side chart of figure 18 that when we give lower necessity to the model, the objective function value improves. Also, the profit and revenue (first objective function) will be more by having more demand. It is apparent from these two charts that while having more demand can be favorable due to having more revenue, if the supply chain cannot cover and manage this demand, it may have a negative effect on customer satisfaction. Because the customers' needs will not be fulfilled, it can decrease the demand amount from customers and consequently revenue in the future.
Additionally, from figure 19, we saw that the more elasticity, the more profit. Discount can play a significant role in controlling our demand and consequently managing our supply chain. If the supply chain cannot cover customers' total demand, price discount can act as leverage and decrease and control the demand. If we are overstock or have much deadstock with one product, the discount leverage can play a role and increase the demand to eliminate our surplus stock. Based on the following sensitivity analysis, the more elasticity we have, the more influential the discount leverage can perform. While we can play with demand using discount leverage in higher elasticity values, unsatisfied demands can also get out of hand more quickly in high elasticities ( Figure 20).

Conclusions
The present research aimed to examine the design of a closed-loop supply chain within an e-commerce platform. We developed a multi-objective, multi-product, multi-period, and multi-layer closed-loop supply chain. We formulated the problem as a mixed-integer non liner programming (MINLP) model to maximize two different objective functions. Since the problem was computationally expensive and np-hard, we used a hybrid NSGA-II algorithm beside our exact method, augmented epsilon constraint. This research's primary purpose was to answer three questions were about decision making with having cost-oriented objectives and customer-centric objectives simultaneously, decision-making in an uncertain environment of e-commerce, and finally, the impact of product's demand uncertain nature on supply chain management.
One of the significant findings of this study is that two conflicting objective functions can improve when the supply chain's limitations are considered. As we saw in the Discussions part, the amount of customers' demand coverage should be planned considering the supply chain's resources and capacities. In this sense, both objective functions can be in the same direction, and we can increase the business's profit while having satisfied customers. However, if the limitations of operational resources or logistics were ignored, increasing the product's supply and demand coverage would negatively impact customer-related objectives. In addition, these findings have significant implications for the understanding of how products' demand nature and their elasticity can play a role in controlling the demand and consequently managing our business.
The generalizability of these results is subject to certain limitations. For instance, various factors have not been seen in this research, such as environmental effects, impacts on employment, and traditional retail. The lack of information on demand also limits the study. In the last decade, big data technologies and machine learning algorithms have witnessed a dramatic shift that predicts and estimates customers' needs and demand more accurately. Despite these limitations, the study certainly adds some information about the dynamics of e-commerce and its decision-makings to our understanding.