Multi-objective modelling of sustainable closed-loop supply chain network with price-sensitive demand and consumer ’ s incentives

.


Introduction
The world's electronic waste has reached to nearly 60 million tonnes by 2021 (Gill, 2021).Increased electric and electronic products consumption rates, improving economies, shorter product lifecycles and limited repair options are main reasons behind this rising e-waste (Forti, Balde, Kuehr, & Bel, 2020).The global electronic waste in 2019 had a material value of almost £46 billion, which is more than the Gross Domestic Product (GDP) of most countries (World economic forum report 2019).Toxic and hazardous materials in e-waste create several healthrelated problems in people like digestives, neurological, respiratory and bone problems.These materials also affect land and sea animals and damages environment by producing greenhouse gases.In 2019, only 17.4% of global e-waste was collected and adequately recycled (Tiseo, 2021).Clearly, cutting global e-waste is essential to curb the climate change impact and achieve Sustainable Development Goals (SDGs) (D'Adamo, Gastaldi, & Rosa, 2020;Prajapati, Kant, & Shankar, 2019).Several companies have started collecting their end-of-life (EOL) products (to focus on reuse, recycle and repair operations) due to growing public awareness and government regulations (Islam & Huda, 2018;Peng, Shen, Liao, Xue, & Wang, 2020;Raza, 2020 ).The recovery of the used products has received considerable attention due to the value creation following the 'circular economy' concept (Coenen, van der Heijden, & van Riel, 2018;Shekarian, 2020).The closed-loop supply chains (CLSC) integrate reverse logistics with the forward logistics network ( Özceylan, Paksoy, & Bektas ¸, 2014;Khatami, Mahootchi, & Farahani, 2015); and this CLSC approach helps to address rising sustainability pressures by decreasing social and environmental impact without affecting profitability (Shekarian, 2020).A considerable amount of literature has expanded around sustainability and the CLSC due to increasing regulations regarding environmental impact and achievement of SDGs (Kazemi, Modak, & Govindan, 2019;Peng et al., 2020;Raza, 2020;Shekarian, 2020).However, only deterministic demand for return of used products has been considered in previous studies, lacking consideration of uncertain markets and conditions of the returned items and incentives (Jena & Sarmah, 2016;Fathollahi-Fard et al., 2021).
Reverse Logistics (RL) plays vital role in CLSC while meeting several needs (Taleizadeh, Haghighi, & Niaki, 2019).The performance of RL mainly depends on the consumers' willingness to return back their used products (Shaharudin et al., 2015).In any society, very limited numbers of people take the proactive approach and practice social responsibility of returning used products to the collection centres.This results into the low performance of collection centres.The effective collection process makes the CLSC network more efficient (Geyer & Doctori Blass, 2010).
Hence, consumers should get some incentives to return their used products.Scholars proposed some polices to increase the return rates by increasing consumers' willingness (Das & Dutta, 2013;Govindan, Soleimani, & Kannan, 2015).Some of these polices affect the demand of new products as well.For example, Das and Dutta (2013) illustrated how the product exchange policy helped to increase total products' demand as well as rate of return products in CLSC.The quality of returned products in RL is another aspect that needs attention from research community.Guide et al. (2003) shown that the buy-back price of returned products should be dependent on their quality.Dutta, Das, Schultmann, and Fröhling (2016) proposed the optimisation model to determine the optimal buy-back price such that the minimum collection limit set by legislators is achieved as well as the total cost is minimised.
The returned products need to be treated differently based on their qualities (various sizes, different ages and diverse way of uses).However, most of the previous studies considered that the returned products can be used for remanufacturing, recycling or any other processes.The quantity of returned products chosen for recovery process depend on the products' quality.The recent reviews of Peng et al. (2020), Prajapati et al. (2019) and Govindan and Soleimani (2017) on CLSCs emphasised the dire need of considering these issues like decisions for returning products, different qualities of returned products, integration of strategic and tactical decisions over a multi-period, price sensitive demand and consumer's incentives.This research study is motivated by a leading European household appliance manufacturer's supply chain network.The volume of e-waste generated per capita in Europe is considerably higher than the America and Asia.Due to growing population, welfare, and changing lifestyle, household appliances are increasing worldwide (Hischier, Reale, Castellani, & Sala, 2020).Household appliances such as microwaves, ovens and refrigerators have the major share in the global e-waste.It is expected that there will be a 37% growth in the use of these appliances by 2020 compared to 2013 (ProMéxico, 2014).Household appliances bring many benefits, but they equally harm the environment, particularly on SDGs (SDG14: Life below water and SDG15: Life on land); if they are not repaired or safely disposed of (Hischier et al., 2020).Recent European Union law for 'Right to repair' household goods (Harrabin, 2021) is encouraging to support the global efforts towards sustainability.The case company is struggling to efficiently manage the supply chain activities and improve its environmental reputation by balancing the economic growth.The European Union has set the target to become the first climate-neutral region by 2050.Hence this company wants to reduce the environmental impact on supply chain activities by designing a sustainable CLSC network.The design of a supply chain (SC) or CLSC network has been well studied (e.g., Badri, Fatemi Ghomi, & Hejazi, 2017;Prakash, Kumar, Soni, Jain, & Rathore, 2020, Fathollahi-Fard et al. 2018a, 2018b, 2021).However, there is limited research focussing on the design of sustainable CLSC while considering wider complexity of supply chain operations, different qualities of returned products, price sensitive demand and consumers incentives (Farshbaf-Geranmayeh, Taheri-Moghadam, & Torabi, 2020;Govindan & Soleimani, 2017;Peng et al., 2020;Yadegari, Alem-Tabriz, & Zandieh, 2019).Undoubtedly, the CLSC is becoming a topical and ever-growing area in recent times due to increased focus on the United Nation's (UN's) SDGs, circular economy, climate change and other adverse economic impacts (Er Kara, Ghadge, & Bititci, 2020;Fathollahi-Fard & Hajiaghaei-Keshteli, 2018;Fathollahi-Fard, Hajiaghaei-Keshteli, & Mirjalili, 2018a;Peng et al., 2020;Shekarian, 2020).The considered company had well-established procedure to manage the societal impact of their operations hence they were not interested to look at those aspects.Therefore, we have not considered the societal aspects like job opportunities and lost working days in our proposed CLSC model.
The rest of the paper is organised as follows.Section 2 discusses the literature on the design of CLSC networks and sustainability to underline the research gap.Section 3 presents the problem description and the case study.The mathematical formulation is presented in Section 4. Section 5 presents the solution approach consisting of the overall procedure on the multi-objective optimisation approach.The data analysis and findings are presented in Section 6. Research and managerial implications are discussed in Section 7. Finally, Section 8 concludes with a discussion on the key findings, limitations and future research recommendations.

Literature review
Supply Chain Network Design (SCND) is one of the important and well-known problems in SC.Scholars mainly explore facility locationallocation and network flows in supply chain network design considering different aspects of the real-world problem.These aspects include product's features, the number of echelons, stakeholders' requirements and transportation.The problem complexity increases with consideration of forward and reverse logistics in SCND.Recently CLSC has been getting attention from researchers and practitioners due to environmental degradation, climate change, growing waste in SCs and the UN's SDGs.Hence, the relevant previous studies mainly focusing on CLSCs and sustainability are discussed in this section.
The objective of past research studies on CLSC was to minimise total costs by establishing an appropriate reverse logistics network; however, several other aspects need to be addressed while developing a robust CLSC network.Dehghan, Nikabadi, Amiri, and Jabbarzadeh (2018) focused on a multi-product, multi-period CLSC for an edible oil supply chain problem; and Yadegari et al. (2019) optimised a CLSC cost considering the inventory cost and time of the PET bottle industry.Whereas Hajipour, Tavana, Di Caprio, Akhgar, and Jabbari (2019) developed a traceable CLSC network implementing an RFID system to curb product losses and lead time while maximising profit.Furthermore, Atabaki, Khamseh, andMohammadi (2019), Farshbaf-Geranmayeh et al. (2020) and Taleizadeh et al. (2019) considered price-sensitive demand while designing CLSC.Interested readers can refer to the literature review on CLSC and related topics published by Peng et al. (2020), Govindan and Soleimani (2017), Prajapati et al. (2019), Kazemi et al. (2019), Coenen et al. (2018), Raza (2020) and Shekarian (2020).
The research work presented in this paper aims to develop a decision support model for multi-echelon, multi-product, multimodal and multiperiod sustainable CLSC while considering multiple supply chain activities of the company, different quality levels of returned products, price-sensitive demand and consumer's incentives.This paper addresses the research gaps identified in Table 1.The proposed model contributes to CLSC literature by providing a holistic modelling approach by simultaneously integrating multiple supply chain factors such as disposal cost, incentive cost, technology selection cost and processing emissions, which are seldom included in the previous literature.The simultaneous consideration of different quality levels of returned products, price-sensitive demand and consumer's incentives are novel features of the current study.The research contributes to knowledge by accommodating wider SC costs and demand to reduce total emissions from transport, processing, disposal and facilities.Due to the complexity of the model, a hybrid multi-objective algorithm integrating NSGA-II and Co-kriging is implemented along with the epsilon constraint method to solve the problem.Data from a leading European household appliance manufacturer are used to assess and validate the proposed model.

Problem description and case study
This section consists of two subsections.The details about the problem statement including product's features, number of echelons, model type and objective functions is given in subsection 3.1.The second subsection 3.2 presents the supply chain network of the case company and related information.

Problem description
This study develops a multi-echelon, multi-model, multi-product, closed-loop supply chain optimisation model while considering the critical aspects inherent in the strategic planning of a sustainable supply chain for the European household appliance industry.We integrate various concepts (reverse logistics, green supply chain design, pricesensitive demand, consumer's incentives, and facility location allocation) into a multi-objective, mixed-integer, non-linear mathematical model.The multi-objective approach establishes a link between facility location and allocation of customers in a cost-effective manner as a first objective, while keeping in mind the second objective of carbon emission minimisation.This enables the focal company to make strategic decisions based on the evaluation of the Pareto frontier points or consider the trade-offs between two objectivestotal cost and carbon emission.The proposed model minimises the cost and environmental emission of the case company.This model helps the case company to achieve the ISO 14064 certification for quantification and reporting of greenhouse gas emissions and removals.We have used the guidelines from ISO 14064 standard while developing the model.Furthermore, we have applied the guidelines of SA8000 from ISO 26000 category for social responsibility.We looked at the health and safety of workers, working hours, remuneration, forced or compulsory workers, child labour and other factors of the case company while developing this model.

Case study
A European household appliance supply chain has a six-tiered system comprising 21 suppliers, 30 manufacturing sites, 8 distributors, 47 retailers, 8 collection centres, 30 customer zones and 5 disposal sites across Europe (Fig. 1).The various symbols used in Fig. 1  The location of suppliers, disposal centres and customer zones are fixed, and potential location of manufacturers, distributors, collection and remanufacturing centres, and retailers are known.The supply chain network has decided to expand by opening additional production centres, distribution facilities, retail, and collection centres.The exact location of these facilities is determined by solving the optimisation model.The suppliers produce five types of raw materials (modules) and supply them to manufacturers who produce the finished products (three types of household appliances).The distributor facilitates the shipment of products between manufacturers and retailers.The end-customers purchases new products from the retail outlets and return old products in exchange for an incentive.It is assumed that the return products may be of different quality levels, determining the incentive price to be paid to the customer.The returned products are sent to the collection centres, where they are refurbished, disassembled, or sent to the disposal centres depending upon their value or usability.The company uses three types of transportation modes (road, rail, and water) and is interested to see the results of three time periods.
For modelling purposes, it is assumed that end customers do not differentiate between new or refurbished products.The suppliers and disposal centres have infinite capacity, while the manufacturer has limited capacity.When faced with unpredictable demand, the manufacturer and distributor hold inventory for modules and products to avoid any stock-out conditions.Provision of alternative three strategies to reduce emissions is available at the supplier, manufacturer, and disposal centres where extensive fuel consumption occurs.The demand at each retailer is a function of potential customers with the price sensitivity of demand and distance travelled by the customer (described in the formulation section).This implies that a retailer may lose a customer if a competitor offers the commodity at a lower price or is near to the customer.Similarly, the number of used products returned depends on potential customers, incentive sensitivity and distance.
The quality of returned products is of four discrete types{i 1 , i 2 , i 3 , i 4 }.It is assumed that products of quality {i 1 , i 2 } can be recycled to a constituent module, which can be reused, while those of quality {i 3 , i 4 } are non-recyclable (1/3rd of the used returned products) and must be disposed of.The reusable modules are refined and transported to the manufacturers as their quality is nearly at par with that of the modules supplied by the suppliers.The commercially returned goods at the collection centres are refurbished and transferred to distributors, where they enter the product flow.

Mathematical model development
Supply chain sustainability depends on the link between economic, social and ecological/environmental performance (Chaabane et al., 2012).However, our research mainly focuses on effective decision-making concerning the company's economic and ecological aspects of global supply chain design and addresses the following objectives: • To develop a multi-echelon, multi-period, multi-commodity forward/reverse logistics network to minimise the total cost of the closed-loop supply chain.• To minimise the total emissions from the transportation and open facilities operations.
Subject to the constraints to determine the appropriate technology selection, mode of transport selection, best location of manufacturing, distribution, retailing, and collection centres from the available potential locations, allocation of retailers to distributors and collection centres, customer zones to retailers, quantity of modules and products flowing across echelons and incentive price for the used products.
Due to the space constraint, mathematical notations are included in Appendix B.
Objective  1) provides the first objective function of the mathematical formulation, which aims to minimise the total cost comprising of the production cost, disposal cost, processing cost related to disassembly and refurbishing, technology cost, fixed location cost, transportation cost, incentive cost related to the return of used product and inventory cost.
Total Production Cost (TPRC) = Total production cost of modules at the supplier's facilities (TPCS) + Total production cost of the products at the manufacturing plants (TPCM).
Equation ( 2) is the production cost of products considering the production cost at the suppliers and manufacturing centres.Equation (2a) presents the total production cost of the modules at the supplier's facilities, and Equation (2b) provides the total production cost of the products at the manufacturing centres.
Total Disposal Cost: Equation ( 3) presents the total disposal cost pertaining to the products of different qualities.
The total processing cost at the collection and remanufacturing centre is provided as the summation of disassembly cost and refurbishing cost.This is shown by equation (4).Equation (4a) determines the total disassembly cost for different products with their respective qualities.Here, the disassembly is performed only when the product quality isi 1 or i 2 .
Equation ( 4b) provides the total refurbishing cost incurred at the collection centre for different products.
Total Technology Cost: Equation ( 5) presents total technology cost comprising of the cost related to the technologies used by the suppliers to produce modules, cost related to the technologies used by the manufacturing centres to produce the products and the cost associated with the use of technologies for performing the disposal of the products at the disposal sites.
Total Fixed Location Cost: Equation ( 6) depicts the total fixed cost for opening manufacturing plants, distributors, retailers, and collection centres at their respective potential locations.
Total Transportation Cost: TTRC is the total transportation cost which can be computed using equation (7).Equation (7a) shows the total transportation cost for the shipment of modules from the suppliers to the manufacturing plants while considering the fixed and variable cost related to the transportation modes.
Equation ( 7b) indicates the transportation cost related to the movement of products from the manufacturing plants to the distributors while considering the fixed and variable costs for the transportation modes.
Equation (7c) provides the transportation cost for the shipment of products from distributor to retailer, considering the fixed and variable costs for the transportation modes.
Equation ( 7d) illustrates the transportation cost related to the movement of products from retailer to the collection centre considering fixed and variable costs related to the transportation modes.
Equation ( 7e) gives transportation cost related to the shipment of refurbished products from the collection centres to the distributors while considering the fixed and variable costs for the transportation modes.
Equation ( 7f) denotes the transportation cost considering the fixed and variable cost for the transportation modes while shipping the disassembled modules from collection centres to manufacturers.
Equation ( 7g) presents transportation costs comprising the fixed and variable cost for the transportation modes.This equation includes the shipment of products from collection centres with quality i 3 , and i 4 which must be disposed of at the disposal site.
Total Incentive Price: Equation ( 8) computes the overall incentive cost considering the incentive price offered for each product type and the number of used products, with different quality returned to the retailer.
Total Inventory Cost.TINC ¼ Inventory holding cost at the manufacturing plants (INCM) + inventory holding cost at the distribution centres (INCN) (9).
Equation ( 9) provides the total inventory cost comprising the inventory holding cost at the manufacturing plants and the distributors.Equation (9a) presents the inventory holding cost at the manufacturing plants for different modules considering the overall planning horizon.
Equation (9b) gives the inventory holding cost incurred for the distributors while storing different product types for the overall planning horizon.

Objective function II:
Minimize total emissions = Production emissions (EP) + Disposal Emissions (ED) + Processing emissions (EPCC) + Transportation emissions (ET) (10) Equation ( 10) presents the second objective function of the mathematical formulation aiming to minimise the total emission incurred from the production of modules, manufacturing of products, disposal of products and emission for performing the disassembly and refurbishing operation at the collection and remanufacturing centre.This equation also includes the emission incurred from the shipment of new modules, new products, returned products, refurbished products, disassembles modules, and disposed products.
Production Emissions (EP) = Emission produced by suppliers during production (TEML) + Emission produced by manufacturers during production (TEMM) In equation ( 11), EP is the total emissions due to the production of the modules by the suppliers and the manufacturing of the products by the manufacturers.
Equation (11a) computes the total emission incurred for the number of modules of different types produced by the suppliers in the overall planning horizon.
Equation (11b) provides the emission incurred for the number of products of various types manufactured by the manufacturers in the overall planning horizon.
Disposal Emissions: Equation ( 12) calculates the emission incurred for the disposal of products of different types at the disposal sites for the overall planning horizon.
Processing Emissions: EPCC = Emission incurred from the disassembly (EMDAS) + Emission incurred from the refurbishing (EMRFB) (13) Equation ( 13) shows the total emissions incurred at the collection and remanufacturing centres for performing disassembly and refurbishing operations.
Equation ( 13a) finds out emission incurred from the disassembly operation for different types of products with specific qualities at the collection centre for the overall planning horizon.
Equation ( 13b) calculates the emission incurred for performing the refurbishing operation of products at the collection centre for the overall planning horizon.
Transportation Emissions: Equation ( 14) presents the total emissions incurred from the transportation of the new modules, new products, returned products, refurbished products, disassembled modules, and disposed products.
Equation (14a) presents the emission incurred from the transportation of modules from suppliers to manufacturers considering different transportation modes for the overall planning horizon.Equation (14b) presents the relationship for determining the overall emission factor considering the constant and variable emission factor for a specific transportation mode and considering the average weight allowed in a specific transportation mode.
Equation (14c) computes the emission incurred from the shipment of products from the manufacturers to the distributors considering different transportation modes for the overall planning horizon.
Equation (14d) depicts the emission incurred from the transportation of products from the distributors to the retailers considering different transportation modes for the overall planning horizon.
Equation ( 14e) calculates the emission incurred from the transportation of used and returned products from the retailers to the collection centres, considering the different transportation modes for the overall planning horizon.
Equation (14f) estimates the emission incurred for the transportation of refurbished products from the collection centres to the distributors using different transportation modes for the overall planning horizon.
Equation (14g) determines the emission incurred for the shipment of disassembled modules from the collection centres to the manufacturers using the available transportation modes for the overall planning horizon.
Equation ( 14h) provides the emission incurred for the transportation of the disposed products of multiple types from the collection and remanufacturing centres to the disposal sites using multiple transportation modes for the overall planning horizon.Here,I ∈ (i 3 , i 4 ).
Technology Constraints: Constraint ( 15) states that a supplier must use only a certain technology type to produce the modules.
Constraint ( 16) makes sure that a manufacturer must use only a particular type of technology to manufacture the products.
Constraint ( 17) depicts that a disposal site must use a particular type of technology for the disposal of the products.
Transportation mode related Constraints: Constraint ( 18) confirms that a specific transportation mode is used for the shipment of modules from a supplier to a manufacturer, given that the manufacturing facility is open.
Constraint ( 19) guarantees that a particular transportation mode is used for the shipment of products from a certain manufacturer to a specific distributor, given that both the manufacturer and the distributor are open at their respective locations.
Like constraint ( 19), constraints (20) to ( 24) depict the transportation mode selection constraint between distributor and retailer, retailer and collection centre, collection centre and distributor, collection centre and manufacturer and collection centre and disposal facility.
Incentive Price: It is imperative for the organisation to give appropriate importance to the product returns, as it plays a vital role in defining the companycustomer relationship.Customer satisfaction from the perspective of product-return experience can substantially increase sales and thereby enhance the revenue generated for the organisation.Therefore, a satisfactory product-return experience can have a significant impact on the profitability of the company.To quantify a customer's long-term value with the organisation, it is essential to ascertain the role of product return in a customer's purchasing decision.Determining the relationship between costs (related to remanufacturing and collection of used products) and benefits of product returns while considering the drivers and the consequences of products returns can help the supply chain managers to distribute the resources and strategically design their supply chain to maximise its profitability (Petersen and Kumar 2010).
The incentive is given to a customer for returning product v of quality i isG i v .LB i v and UB i v are the lower and upper bounds of the incentive pricing respectively.These lower and upper bounds are a function of the price of the original product pr v and the quality of returned product i.The lower and upper bounds of the incentive pricing are calculated using equations ( 25b) and (25c).
Capacity Constraint for the Manufacturer: Constraint ( 26) ensures that the number of products of each type transported from a specific manufacturer to all distributors should be less than or equal to the manufacturer's storage capacity for a certain product type.
Demand at Retailer: Equation ( 27a) and (27b) are used to determine the demand allocation for each product type of the customer zone to the retailer in each period.Equation (28) depicts the demand constraint in the terms of the total number of products shipped from various distribution centers to the specific retailer should be equal to the demand of the product at the retailer.Equation (29) presents the relationship between the distance of the customer zone and the retailer with the portion of people (expressed by a value within 0 to 1) willing to come to the retailer from the customer zone for buying the product (Kaya & Urek, 2016).
Product Returned to Retailer: Equation ( 30) determines the number of used products of each type with certain quality returned to the retailer from a customer zone in a specific period.Equation (31) helps estimate the portion of people (expressed by a value within 0 to 1) willing to come to the retailer from the customer zone to return the product.Equation (32) computes the total number of used products of specific quality returned to a retailer.
Equation ( 33) estimates the number of commercially returned products of each type to the retailers in each period.
Constraint (34) helps to estimate the number of quality products and disassembles at the collection centres in each period.
Equation ( 35) estimates the number of used and returned products with quality i 3 and i 4 transported from the collection centre to the disposal centre for the disposal of the product.
Equation ( 36) determines the number of products to be refurbished, commercially returned from the retailer to the collection and remanufacturing centre.
Equation ( 36) states that the used product number returned from the retailer to the collection centre comprises the number of commercially returned products at the retailer and the number of used products returned to the retailer with qualityi 3 and i 4 .
Equation ( 38) makes sure that the total number of products moving from a collection center to the distributors in a time period is equal to the total number of products of a particular type refurbished at the collection centre.
Equation ( 39) defines that the total number of modules shipped from a collection centre to manufacturers in a given time depends on the total number of modules obtained from the products disassembled at the collection centre.Equation (40) computes the number of modules obtained at the collection centre after performing the disassembling operation.Mogale et al. Equation (41) shows the relationship between the number of modules required at the manufacture and the total number of products produced at the manufacturing plant.Equation ( 42) ensures that the initial inventory of modules at the manufacturing plant is assumed to be zero.Equation ( 43) helps to determine the inventory of modules at the manufacturing plant for each period.
Equation ( 44) illustrates that the number of products shipped to the retailer using the forward network should be less than or equal to the total number of products at the distributor, which includes new products coming from the manufacturers, refurbished products coming from the collection centre and the inventory of the distributor for the specific product.Equation ( 45) confirms that the initial inventory of the distributor for the specific type of product is assumed to be zero.Equation ( 46) helps to determine the inventory of the distributor for each product type for each period.
In this paper, a well-known epsilon constraint method is applied to solve the developed multi-objective mathematical model.This method can solve the problem instance of the case study effectively.However, to deal with more complex problems in the future, the supply chain case was also interested in metaheuristic solution.Hence, this study overcomes the complexity challenge by developing a two-stage solution approach that hybridises NSGA-II evolutionary algorithm with Surrogate Modelling through Co-Kriging.The proposed hybrid algorithm tests the effectiveness of the results against those obtained from the epsilon constraint method.While the epsilon constraint method obtains the optimal Pareto front, and the NSGA-II is applied to generate initial points of the Pareto front.The Co-Kriging surrogate model interpolates the Pareto frontier using these solution points.

Initial feasible solution generation
The initial feasible solution is obtained by determining the values of the decision variables of the mathematical model.We have taken into consideration the approach adopted by Fathollahi- Fard et al. (2021) in establishing a connection between optimization model and the search space of the algorithm.Initially the binary variables of the mathematical model such as X j l (whether technology j to be used by supplier l or not) and X j s (whether technology j to be used by disposal site s or not) are randomly generated considering the technology constraints given in equations ( 15) and ( 17) respectively.Hence, we obtain the information highlighting the type of technology j to be used by all the suppliers and disposal sites.Initially, the majority of the manufacturing plants,  16), the value of binary variable X j m (whether technology j to be used by manufacturer m or not) is randomly obtained.
If the manufacture is open at location m orY m = 1, in that case using the transportation mode related constraint given in equation ( 18) becomes ∑ w Z wt lm = 1, and accordingly, using the equation we randomly generated the value of the binary variable Z wt lm (whether transportation model w to be used for shipping products from supplier l to manufacturer m in time period t).Furthemore, the value of binary variablesZ wt mn ,Z wt np ,Z wt pr , Z wt rn and Z wt rm are obtained using equations ( 19), ( 20), ( 21), ( 22) and ( 23) respectively, while considering the values of binary variablesY m ,Y n , Y p andY r .Similarly, the value of binary variables Z wt rs is randomly generated using equation ( 24), only when Y r = 1 and the equation ( 24) becomes ∑ w Z wt rs = 1.Now, if a tranportation mode is available during time period t for shipping products from distributor n to retailer p, then the binary variable A t np takes a value 1. Therefore,A t np = 1, if Similarly, binary variables such as A t pr takes a value 1 only when the binary variable Z wt pr is 1.Hence,A t pr = 1, if ∑ w Z wt pr = 1 forp ∈ P, r ∈ R, t ∈ T. The value of the demand of product v from customer zone q allocated to retailer p in time period t, or D t pqv is obtained using equation ( 27b) and the value of parameter dp t pv (which is the demand of product v arising with retailer p in time period t).Now, depending upon the value of γ pq (it is a parameter with a value between 0 and 1, highlighting the willingnes of customer from zone q to come to retailer p for buying products), we generated the value of binary variable A t pq (which highlights whether retailer p serves customer zone p or not), using equation (27a) and the variableD t pqv .The values of decision variablesD t pqv , A t np and Y n are used in equation ( 28) to determine the quantity of product v transported from distributor n to retailer p in time period t or decision variableK t npv .Now, the incentive offered for product v of quality i or value of the decision variable G i v is randomly generated while considered the lower and upper bound as given in equation (25a).Using the value of γ pq and σ (which is the fraction of sales equaling to the number of returned used products), we compute the value of δ pq from equation (31).Now, considering the value ofδ pq ,A t pq ,G i v , and the values of parameters ρ t q andλ, we determine the value of the decision variable K it qpv (it is the quantity of used product v of quality i returned to retailer p from customer zone q in time period t) considering equation ( 30).The quantity of used product v of quality i returned to retailer p in time period t, or decision variable Kr it pv is determined from equation (32) while using the value of varia-bleK it qpv .The value of the decision variable Kda it rv (it is the quantity of product v of quality i disassembled at collection centre r in time period t) is computed using equation (34) while considering the value of decision variables Kr it pv andA t pr .Equation ( 35) helps to estimate the quantity of product v of quality i shipped from collection centre r to disposal centre s in time period t, or K it rsv while using the values of decision variables Kr it pv andA t pr .The quantity of commercially returned products v at retailer p in time period t or the value of decision variable Kc t pv is obtained using equation (33) while considering the values of parameters χ p (it is the fraction determining the number of commercially returned goods at the retailer p) and dp t pv (it is the demand of product v arising at retailer p in time period t).Now the quantity of used and returned product v transported from retailer p to collection centre r in time period t, K t prv is obtained using equation ( 37) which taking into account the values of decision variablesKc t pv , Kr it pv andA t pr .Now, the quantity of product v refurbished at collection centre r in time period t, Krb t rv is computed by using the values of decision variables Kc t pv and A t pr within equation ( 36).Furthermore, the value of variable Krb t rv and equation ( 38) is used to determine the quantity of refurbished product v transported from collection centre r to distributor n in time period t, or decision variableK t rnv .The quantity of module u obtained after performing the disassembling operation at the collection centre r in time period t, Ku t ru is generated using equation ( 40) while considering the values of various parameters such asθ, np vu and decision varia-bleKda it rv .The value of the decision variable K t rmu or the quantity of disassembled module u shipped from collection centre r to manufacturer m in time period t is computed using equation ( 39) and the decision variableKu t ru .The quantity of product v to be transported from manufacturer m to distributor n in time period t, K t mnv is randomly generated while ensuring that value of K t mnv is less than or equal to a mv Y m and thereby satisfying equation ( 26), which is also the capacity constriant of the manufacturer.Now, the value of decision variable H t nv or the inventory of product v at distributor n in time period t is generated using equation ( 44) while considering the values of the decision variablesK t mnv ,K t rnv , K t npv andA t np .Furthermore, the inventory of module u at manufacturing plant m in timer period t, or H t mu is determined by considering equations ( 43) and ( 41).Using equation ( 43), the following can be written, or,H 1 mu = 0. Equation ( 41) can be represented as . Therefore, the value of the decision variable K t lmu (or the quantity of module u shipped from supplier s to manufacturer m in time period t) and H t mu is determined using equation ( 43) and ( 41).Finally, the initial feasible solution for the algorithm is obtained which contains the value of the decision variables of the mathematical model.The initial feasible solution generation takes into account the different constraints of the mathematical model and accordingly determines the feasible solution.

Epsilon constraint method
Epsilon constraint method is applied to obtain a set of Pareto optimal solution for the sustainable closed-loop supply chain network problem.This is a one-stage technique, unlike the two-stage NSGA-II and Co-Kriging approach.It is based on the optimisation of a single objective while considering the other objective as a constraint with allowable bounds.Then, the bounds are sequentially modified to obtain other Pareto-optimal solutions.A set of Pareto-optimal solution is obtained by altering the value of epsilon.Finally, based on the decision maker's perspective, the desired compromise solution can be chosen from the Pareto-optimal solutions.

Kriging
Considering the inherent computational complexity in multiobjective problems with many variables and constraints, researchers have focussed largely on approximating Pareto frontiers using initial solutions (Couckuyt, Deschrijver, & Dhaene, 2014).Dixit et al. (2016) used NSGA-II combined with the Co-Kriging approach to solve a performance-based optimisation problem related to supply chain network resilience.On similar lines, this paper focuses on developing a novel application of this approach for solving a CLSC optimisation problem.
Kriging is a Gaussian process interpolation method used to develop approximations from sample data.This surrogate modelling technique is effective in approximating the Pareto frontier in multi-objective opti-misation problems (Couckuyt et al., 2014).It comprises a regression termf(x), which explains the largest variance in sample data, a Gaussian process R(x) with mean zero, variance σ 2 and a correlation matrixρ.The Gaussian process enables the interpolation of residuals and is represented in equation ( 49) and equation ( 50).
where b i (x) are basis functions of x, ∀i ∈ [1, m] and α i , ∀i ∈ [1, m] represent the contribution of each basis function in determiningf(x).Thus, from equation ( 50), f(x) can be defined as the linear combination of the basic functions.For n samples, part of Kriging model is a n × m model matrix, as represented in equation (51).The stochastic process is defined by the n × n correlation matrix ρ as presented in equation ( 52).
Readers are advised to refer the above and the documentation of Oodace toolbox (Couckuyt et al., 2014) for further explanation of the model.

NSGA II + Co-Kriging
A two-stage approach to solve the multi-objective mixed integer programming model of the CLSC is adopted in this study.NSGA-II is used to determine the Pareto frontier points of the bi-objective model: minimisation of total cost and emissions.The obtained Pareto front is arranged in descending order of crowding distance (Deb et al., 2002).The densely packed points (low crowding distance) are categorised as coarse points and it provides information about a limited space of the Pareto frontier.Data points with high crowding distance are classified as fine points and they provide significant information about further interpolation of the Pareto frontier, hence, decide the shape of the frontier.
The output of fine and coarse data points from the first stage is fed as input to the second stage, determining the interpolated Pareto frontier through Co-Kriging.This approach exploits the correlation between the two types of data points to enrich Pareto frontier and enhance prediction accuracy

μ(x)
The hybrid NSGA-II and Co-Kriging approach is computationally efficient and can handle multi-objective formulations effectively (Dixit et al., 2016).The flowchart depicting the combined NSGA-II and Co-Kriging approach is shown in Fig. 2. It develops enriched Pareto frontier by employing quick interpolation.It also provides variance plots with a degree of confidence associated with it to enable managers to make informed decisions from the interpolated Pareto frontier.

Results and analysis
Results and analysis are discussed in this section.Initially, various type of data and data collection method is explained.Next, we talked about how the mathematical model is solved using an exact and metaheuristic algorithm.The results of facility locations and transportation mode shares are delineated in subsection 6.3 and 6.4 respectively.Subsection 6.5 presents the results of forward and reverse logistics costs and emissions.The optimal incentive prices for returned products are illustrated in subsection 6.6.The objective functions solutions obtained through the NSGA-II algorithm and epsilon constraint method are compared in subsection 6.7.The last subsection 6.8 depicts the Cokriging and variance plots.

Data collection
Most of the data used for this study came from a leading European household appliance supply chain company.Other relevant data related to environmental parameters were collected from the literature sources.The model was analysed using a commercial software-MATLAB 2018a on a computer with the following specifications: Intel® core™ i7-5500U CPU processor @2.40 GHz, 8.00 GB RAM.Several numerical experiments are implemented to assess the performance of the model and thereby, the effectiveness of this model is established after obtaining the results.
For the transportation mode data, which includes the fixed cost, variable cost, constant and variable emissions factors, the values are presented in Table C1 (Appendix C).The emission-related data was taken from the relevant literature sources.The sensitivity parameters are as follows: price sensitivity = 0.01 and incentive sensitivity = 0.01, ω = 0.5 (Kaya & Urek, 2016).The fraction of commercially returned products at the retailer, FR, is 12%.The fraction of used products to be disposed of FD, is 1/3 as advised by the concerned managers of the case company.
The unit disassembly cost and disposal cost for returned products is presented in Table C2 in Appendix C. The company managers supplied the specific unit production cost of modules and products, and unit refurbishing cost.The approximate unit cost and the weight of the modules at the supplier cities and finished products at the manufacturers were provided by the company managers from the supply chain record books maintained by the firm.The weight of module u 1 , u 2 , u 3 , u 4 and u 5 is 0.25, 0.1, 0.125, 0.05, and 0.07 (kg), respectively.Similarly, the weight of the product v 1 , v 2 and v 3 is 0.83, 0.53 and 0.36 (kg), respectively.They also provided the product-modules matrix to determine the combination of modules required to develop each product (Table C3 in Appendix C).
The emissions due to production of modules at suppliers and finished goods at manufacturers, disposal, disassembly, and refurbishing of returned items were taken from Kim et al. (2008).The selling price of new product V 1 , V 2 and V 3 is £28, £22 and £18 respectively.Similarly, lower and upper bounds of incentive price in our experimental analysis were finalised.Here, e and f are the constants for lower and upper bounds, respectively and for our closed-loop supply chain problem, e and f are assumed as e = 3, f = 1.Using equations ( 56) and ( 57), the lower and upper bound of the incentive price is determined and presented in Table C4 in Appendix C.

Solving the model
Analysis was conducted utilising the above developed model and data collected from a leading household appliance company based in Europe.The CLSC network of the company with the largest market share and production volume is considered for the proposed mathematical model.
It can be observed that the number of decision variables and constraints related to the problem instance considered for solving the proposed mathematical formulation are 115,014 and 196,292, respectively.Due to this enormous number of decision variables for the CLSC network problem, the exact methods can take a significant amount of computational time to solve such a problem (De et al., 2020;Mogale et al., 2020).For example, the epsilon constraint method took approximately 8397.207 s to solve the proposed problem.Hence, the study employed evolutionary NSGA-II coupled with a Co-Kriging modelling approach which utilises the correlation among input points to interpolate an enriched frontier in large scale problems.The proposed approach for this model developed a Pareto frontier within 3521 s.Moreover, it can be easily understood that using exact methods can tremendously increase the time complexity of the process.Consequently, the vast amount of data available today necessitates the integration of predicting mechanisms like Kriging and stochastic models in traditional evolutionary optimisation techniques to converge to an effective solution space within a shorter time duration.

Facility location results
In this section, the number of open facility location of manufacturer, distributor, retailer, and collection and remanufacturing centre in the CLSC are determined from the Pareto solutions, as presented in Table 2. Since Co-Kriging performs interpolation for inputs of the form (total cost, total emission), The following table provides facility location information obtained from the solutions of epsilon constraint and NSGA-II.The Pareto solutions in Table 2 are arranged in order of increasing cost and decreasing emissions.
It is evident from this table that the number of open facilities are increasing downwards which results in higher cost mainly due to a greater number of open facilities.However, emissions decrease because there are more facilities to cater to the rising demands.This claim is supported by Figs. 3 and 4, highlighting that total facility location cost is the major contributor to total cost while transportation emissions are prime contributors to total carbon emissions.

Transport mode shares
This section shows that an intermodal transport network is not always conducive to the company's profitability (Agamez-Arias & Moyano-Fuentes, 2017).In other words, it is not the best, most costeffective strategy to control emissions.Table 3 shows that though rail and water transport modes reduce emissions to a huge extent, the cost incurred in employing them is not feasible relative to road transport which offers minimum logistics cost.Another important insight derived from this case study is that the shift from road to rail transport modes offers the biggest emissions reduction.However, this fact cannot be generalised to all types of supply chain problems because of location restrictions on railway terminal availability.For example, if the rail terminal is far from the facility point, then road transport must be employed to transport the products.Nevertheless, we can conclude that, for this model, a shift from road to rail and water to rail transportation reduces emissions significantly.

Forward and reverse logistics cost and emissions
This section presents the various types of cost incurred by the firm while operating a CLSC.It is evident from Figs. 3 and 4 that facility location and transportation costs are major contributors to the total cost and, together, they constitute 59% of total cost.Similarly, transportation and manufacturing emissions contribute 64% to total emissions.From Fig. 5 it can be stated that the forward supply chain incurs more cost and emissions than the reverse supply chain.This implies that managers should try to promote reverse logistics as it is aligned with the concept of a green supply chain.
Accurate recycling, refurbishing and disposal of used items effectively lowers emissions by reducing reliance on the production of new products to an appreciable degree.In the model, ERM ≫>≫ ERN > ERS implies that most of the products returned by customers are recyclable and, hence, incur less emissions due to the transportation of commercially returned and disposable items.Here, ERM refers to the carbon emission incurred for the shipment of disassembled modules from the collection centres to the manufacturers.ERN refers to the carbon emission incurred for the transportation of refurbished products from the collection centres to the distributors.ERS refers to the carbon emissions incurred for the transportation of the disposed products from the collection and remanufacturing centres to the disposal sites.
Emissions due to processing are minimal; therefore, investing in technology at the collection centre is not a feasible option as it unnecessarily adds to the technology/investment cost.However, investment in technology at the disposal centre is justified because even though the total number of products to be disposed of is less, the carbon emissions due to the disposable products are considerably higher.Investment in technology at the disposal centre helps to lower the carbon emissions incurred from the disposable products.Transportation emissions in forward logistics (ELM + EMN + ENP) are considerably higher than in reverse logistics (EPR + ERM + ERN + ERS) because the total number of products transferred in forward logistics is higher.Here, ELM refers to the carbon emissions incurred from the transportation of modules from suppliers to manufacturers, EMN and ENP refer to carbon emission incurred from the shipment of products from the manufacturers to the distributors and the transportation of products from the distributors to the retailers, respectively.EPR refers to the carbon emissions incurred from the transportation of used and returned products from the retailers to the collection and remanufacturing centres.

Optimal incentive price for returned products
In this section, incentive prices for the three types of returned products with four qualities each are assessed.It is assumed that the retailer has the discretion to decide the quality of the returned products in accordance with the specified standards and provide appropriate incentives to customers.We can infer from Fig. 6 that for quality i 1 and i 2 of each product, the optimal incentive price is more towards the upper bound and for quality i 3 and i 4 it is towards the lower bound.This shows the organisation's inclination to lure customers to return used products that are in good condition.
Moreover, it can recover costly goods of acceptable quality at prices below 50% of the price of new products.Hence, they can be brought back into the market circulation at a lower cost to the company.Although the firm pays slightly more for these products than others, as cost of refurbishing is lower, it helping to generate higher profit.On the other hand, the low-quality returned products are a liability to the firm as incineration and disposal of these products is very costly, considering the government's environmental regulations.However, even for lowquality used products, the incentive provided throws light upon the policymaker's perspective of promoting the facilities to recover a certain proportion of their products and ensure proper disposal.

Comparison of objective function solutions of NSGA-II and epsilon constraint
In Table 4 we can see that NSGA-II + Co-Kriging converges to an effective solution compared to the solution obtained from epsilon constraint with a significant reduction in computational time.The difference between the average total cost obtained using NSGA-II + Cokriging algorithm and the epsilon constraint method is reported as a solution gap for the first objective.Similarly, the solution gap for the second objective is also determined.The solution gap for the first objective is 1.47% and 1.72% for the second objective, proving the near optimality of the integrated approach.Also, the Pareto front approximated by NSGA-II + Co-Kriging is quite close to that of epsilon constraint.The red numerical values denote the points predicted by Co-Kriging using NSGA-2 inputs.The Pareto solution points obtained from the integrated approach closely mimic the epsilon constraint solution.This aspect advocates the application of Co-Kriging for multi-objective problems with many variables.The efficiency of Co-Kriging is solely based on the fact that it employs the correlation between the input points to quickly predict and interpolate the Pareto frontier in the solution space.Thus, it enhances the decision-making capabilities of managers.

Co-Kriging and variance plots
The computationally challenging bi-objective problem formulation is solved using NSGA-II for 50 iterations considering 200 initial solutions to generate two fine and five coarse points.The red squares in Fig. 7a denote the input points obtained by NSGA-II.In this figure, the red dotted line is the Pareto frontier interpolated by Co-Kriging to join the NSGA-II input points.At the same time, the thick black curve represents the interpolation of the Pareto frontier into farther space.The Co-Kriging plot was obtained instantaneously, which establishes that the integrated approach is computationally less expensive and efficaciously    deals with multi-objective formulations.Besides, the Co-Kriging provides a variance plot for predicting the variance in total emissions with respect to the total cost incurred, as shown in Fig. 7b.The NSGA-II + Cokriging approach provides an eclectic combination of two conflicting objectives through an enriched Pareto front and enables decisionmakers to select the best trade-off solution with an appropriate degree of prediction accuracy.
In the variance plot (Fig. 7b), the variance for total emissions is minimum in the range [3.88 × 10 6 £, 4.15 × 10 6 £] which implies that the total emissions can be predicted with maximum accuracy for this range of total cost.The variance plot poses a very interesting conclusion.We note that variance is minimum in the region where NSGA-II points are defined.However, moving away from the points, Co-Kriging tries to predict the Pareto frontier using the correlation among the input points.Moving farther, the correlation between the input and predicted points gradually decreases; hence, variance increases gradually on either side, as seen in Fig. 7b.
The integrated Co-Kriging approach reduces our reliance on computationally time-consuming NSGA-II to give an enriched solution.Since it considers only initial coarse and fine data points, it significantly reduces the total number of simulations (Dixit et al., 2016).Managers    can make informed decisions by evaluating the trade-off between the objective functions based on the efficient interpolated frontier, because Co-Kriging interpolates the Pareto front even in the regions where NSGA-II points are not present.Decision-makers can employ the variance plots to discern additional information regarding the prediction accuracy of each Pareto frontier point.

Research and managerial implications
This section initially presents our major contributions and findings in terms of research implications in subsection 7.1.Then, we have discussed how the proposed model and findings are beneficial to the managers to make informed decisions and choices.

Research implications
This paper makes several implications in terms of mathematical model development and application of hybrid metaheuristic algorithm.Lack of simultaneous consideration of multi-product, multi-mode and multi-period characteristics in the mathematical model is evident from existing studies (Dehghan et al., 2018;Nayeri et al., 2020).Most of the authors focused on a single economic objective of the sustainable development of industry sectors; thus, the environmental or green perspective is still limited (Hajipour et al., 2019;Shekarian, 2020;Yadegari et al., 2019).Most of the real-life problems are multi-objective in nature; hence researchers should pay more attention to the development of multi-objective models.Scholars should provide several Pareto optimal solutions to decision-makers, which will help them to make informed decisions (Mogale et al., 2020).The novel non-linear programming models and methods -particularly heuristics, hybrid algorithms and metaheuristics need to be developed to deal with the real-life complicated and complex CLSC problems (Peng et al., 2020;Govindan et al., 2015).The heuristic, metaheuristic and simulation studies are more realistic and applicable in practice than exact methods (Kim, Chung, Kang, & Jeong, 2018;Sahebjamnia et al., 2018).

Managerial implications
The proposed decision support model enables managers to make informed choices and determine the trade-off between cost and emissions.Large scale companies can implement this decision support model to reduce the burden on their raw material resources since used products can be disassembled and brought into the supply chain circulation of products.Such steps can reduce up to 1.5 million tonnes of electrical waste generated in the UK each year, thus contributing significantly to reduce carbon emissions and landfill costs (Harrabin, 2020).This study also concludes that a shift from road to rail transportation and water to rail transportation reduces emissions significantly.Since several international governmental bodies (e.g., COP 26) are taking stringent measures to urge large organisations to adopt appropriate steps to curb carbon emissions, there is burgeoning scope for the application of efficient technology and intermodal transportation to develop a sustainable business environment.
The development of enriched Pareto fronts, in less computational time, justifies the use of predictive modelling along with evolutionary algorithms.The variance plots enable managers to ascertain the range in which they can make decisions with high prediction accuracy.This paper sheds light on the effects of incentive pricing of returned goods in reverse logistics network.This helps to increase the rate of returned products from consumers.The decision-makers can do the planning of returned products based on the analysis of the incentive pricing results and different quality levels of product.The analytical approach to solving complex problem provides us with a trade-off solution, which can help managers make strategic decisions under dynamic environments.Furthermore, the developed multi-objective approach offers the possibility to respond to governmental legislation while simultaneously optimising the total cost components of their CLSC.

Conclusion and future scope of research
Keeping sustainability principles in mind, this paper presented a multi-product, multi-period, mixed-integer, mathematical optimisation model for sustainable CLSC.Initially, the model was solved using the epsilon constraint method.Considering a large number of decision variables and constraints in the model, a hybrid meta-heuristic approach (e.g., NSGA-II + Co-Kriging a surrogate modelling technique) was used to determine the optimal trade-off between costs and emissions in the CLSC network in minimum computational time.The NSGA-II was run for a small number of simulations to obtain coarse and fine points, which are the inputs to Co-Kriging resulting in Pareto solutions.The solution gap of 1.47% from optimal solution of total cost and 1.72% from optimal solution of total emission clearly shows that this integrated approach closely approximates the results of the epsilon constraint approach and is computationally less expensive.The results show that costs and emissions in forward logistics supply chains are significantly greater than those in a reverse supply chain.This is not a unique finding, but the implications can play a major role in convincing the supply chain policymakers of firms to incorporate reverse logistics into their operations and, hence, promote remanufacturing of goods.The major contributors to total cost include transportation cost, facility location technology, production, and manufacturing cost, while the major contributors to emission consist of production and manufacturing-based emissions.Analysis of the incentive pricing results shows that policymakers should fix the incentives of Quality 1 and Quality 2 of each product towards the upper bound and the incentives of Quality 3 and Quality 4 more towards the lower bound.
Similar to other studies, the current study has some limitations, which provides avenues for future research.Discussed model was assessed based on the data collected from a single company and thus its difficult to generalise the findings.Further improvements to the model can be made by incorporating uncertainty of parameters to develop a multi-objective stochastic optimisation problem.Our future scope will focus on integrating bounded rational decisions and risk pooling effects in the present model.Another research avenue could be the development or modification of existing algorithms to make the solution computationally inexpensive.The application of recently developed algorithms like red deer algorithm (Fathollahi-Fard et al., 2020c), whale optimisation algorithm (Fathollahi-Fard et al., 2021) and social engineering optimiser (Fathollahi-Fard et al., 2020a) can be tested on the current model.This research can be extended to the multi-greenhouse gas emissions case by normalising other greenhouse gases using climate change effects.The impact of logistics and supply chain networks on the social factors of sustainability, such as number of job opportunities created, lost days due to various reasons, public health, accidents, and noise pollution could be explored (Mogale et al., 2020).The implications of various standards like ISO 14000 for environmental performance and ISO 26000 for social responsibility on CLSC network remains unexplored and provides another avenue for future research.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence D.G. Mogale et al. the work reported in this paper.Average weight allowed in transportation mode w ρ t q Potential customers in customer zone q in time period t Ω Coefficient of price sensitivity of demand γ pq Parameter between 0 and 1 based on distances between retailer p and customer zone q to estimate the portion of people going from q to p to buy products dp t pv Demand of product v arising at retailer p in time period t λ Incentive sensitivity of collected amount δpq Parameter between 0 and 1 based on distances between retailer p and customer zone q to determine the portion of people going from q to p to return products Equals to 1 if technology j is used by supplier l, otherwise 0 Equals to 1 if technology j is used by manufacturer m, otherwise 0

X j s
Equals to 1 if technology j is used by disposal site s, otherwise 0

Ym
Equals to 1 if manufacturing plant is opened at location m, otherwise 0 Yn Equals to 1 if distribution centre is opened at location n, otherwise 0 Yp Equals to 1 if a retail facility is opened at location p, otherwise 0 Yr Equals to 1 if a collection centre is opened at location r, otherwise 0 Z wt lm Equals to 1 if transportation mode w is used for shipping modules from supplier l to manufacturer m in time period t, otherwise 0. Z wt mn Equals to 1 if transportation mode w is used for shipping products from manufacturer m to distributor n in time period t, otherwise 0.

Z wt np
Equals to 1 if transportation mode w is used for shipping products from distributor n to retailer p in time period t, otherwise 0.

Z wt pr
Equals to 1 if transportation mode w is used for shipping return and used products from retailer p to collection centre r in time period t, otherwise 0.

Z wt rn
Equals to 1 if transportation mode w is used for shipping refurbished products from collection centre r to distributor n in time period t, otherwise 0.

Z wt rm
Equals to 1 if transportation mode w is used for shipping disassembled products from collection centre r to manufacturer m in time period t, otherwise 0.

Z wt rs
Equals to 1 if transportation mode w is used for shipping disposed products from collection centre r to disposal centre s in time period t, otherwise 0.

A t np
Equals to 1 if in the forward network, retailer p receives shipment from distributor n in time period t, otherwise 0

A t pr
Equals to 1 if in the reverse network, retailer p returns product to collection centre r in time period t, otherwise 0 A t pq Equals to 1 if customer zone q is served by retailer p in time period t, otherwise 0 Continuous variables:

Table C1
Data related to the parameters of the mathematical model.

Table C3
Modules requirement for each Product.

Table C4
Bounds on incentive price to customers.
are described below.1: Transportation of raw materials (modules), 2: Transportation of manufactured products, 3: Transportation of products, X: Purchase and return of goods by customers, A: Transfer of commercially returned and used products, B: Transfer of refurbished goods, C: Transfer of disassembled modules from old products and D: Scrapped products.
. The function values of coarse and fine data points is represented as y c = ( of a model Y c (x) from coarse points (X c , y c ), followed by developing a second Kriging model Y d (x) on the residuals of coarse and fine points (X e , y d ), where y d = y e − ω.μ c (X e ).The parameter ω is estimated as part of the MLE of the second Kriging model.The choice of the correlation and regression function of both Kriging models can be adjusted separately for the coarse data and the residuals, respectively.The Co-Kriging model is defined by equations (53)-(55).The block matrices M(M c , M d ), F(F c , F d ), cor(x) and ρ(ρ c , ρ d ) are used in two separate kriging models μ c (x) and μ d (x), respectively.In this model, σ 2 c and σ 2 d are the process variances of the two Kriging models Y c (x) and Y d (x) respectively.

Fig. 5 .
Fig. 5. Forward and Reverse Supply Chain cost and emission.

uv
Fraction determining the number of commercially returned goods at the retailer p θ Fraction of used products returned which are apt for disposal σ Fraction of sales equalling to the number of returned used products npvu Number of modules of u used for making one unit of product v ch t Inventory holding cost for per unit module u per time period t ch t Inventory holding cost for per unit product v per time period t prv Price of the original product v

v
u shipped from supplier s to manufacturer m in time period t K t mnv Quantity of product v transported from manufacturer m to distributor n in time period t K t npv Quantity of product v transported from distributor n to retailer p in time period t K t prv Quantity of used and returned product v transported from retailer p to collection centre r in time period t K t rnv Quantity of refurbished product v transported from collection centre r to distributor n in time period t K t rmu Quantity of disassembled module u shipped from collection centre r to manufacturer m in time period t K it rsv Quantity of product v of quality i shipped from collection centre r to disposal centre s in time period t K it qpv Quantity of used product v of quality i returned to retailer p from customer zone q in time period t Kda it rv Quantity of product v of quality i disassembled at collection centre r in time period t Krb t rv Quantity of product v refurbished at collection centre r in time period t Ku t ru Quantity of module u obtained after performing the disassembling operation at the collection centre r in time period t Kc t pv Quantity of commercially returned products v at retailer p in time period t Kr it pv Quantity of used product v of quality i returned to retailer p in time period t D t pqv Demand of product v from customer zone q allocated to retailer p in time period t H t mu Inventory of module u at manufacturing plant m in timer period t H t nv Inventory of product v at distributor n in timer period t G i Incentive price offered for product v of quality i Appendix C. Data See Table C1-C4.

Table 1
Overview of relevant CLSC models in literature.
centres, retail facilities and collection centres are kept open and hence the value of respective binary variablesY m ,Y n , Y p and Y r are obtained accordingly.Using the values ofY m , and technology constraint given in equation ( D.G.Mogale et al.distribution

Table 2
Open facility location corresponding to Pareto solution.

Table 3
Shares of cost and emission for each transportation mode.

Table 4
Comparison of solutions of epsilon constraint and NSGA-II.