Robust Optimization for a Bi-Objective Green Closed-Loop Supply Chain with Heterogeneous Transportation System and Presorting Consideration

: In this study, we propose a robust bi-objective optimization model of the green closed-loop supply chain network considering presorting, a heterogeneous transportation system, and carbon emissions. The proposed model is an uncertain bi-objective mixed-integer linear optimization model that maximizes proﬁt and minimizes carbon emissions by considering uncertain costs, selling price, and carbon emissions. The robust optimization approach is implemented using the combined interval and polyhedral, “Interval+ Polyhedral,” uncertainty set to develop the robust counterpart of the proposed model. Robust Pareto optimal solutions are obtained using a lexicographic weighted Tchebycheff optimization approach of the bi-objective model. Intensive computational experiments are conducted and a robust Pareto optimal front is obtained with a probability guarantee that the constraints containing uncertain parameters are not violated (constraint satisfaction).


Introduction
In April of 2018, the concentration of carbon dioxide, CO 2 , in the atmosphere reached 410 parts per million (ppm), which was the highest level in the previous 800,000 years; compared to the previous highest level in 1780 which was just 280 ppm [1,2]. Many crises have resulted from carbon emissions, including rising sea levels and accelerating species extinction [2,3].
The transportation sector was the second-largest source of CO 2 emissions in 2018, around 25% of the total carbon emissions, according to a report by the International Energy Agency [4]. As a result, limiting transportation's excessive CO 2 emission growth has become a critical component of global CO 2 emission reduction [5].
Many governments pay attention to conducting measures and refining legislation to motivate organizations and a broad sector of customers to adopt green and sustainable practices in their production and service activities [6,7]. Approximately 197 countries at the UN climate summit in November 2021 agreed on new measures for carbon emission reduction, including standards for reporting national emissions and the ground rules for trading credits [8].
A green closed-loop supply chain (GCLSC) involves the integration of environmental aspects into all activities of the traditional closed-loop supply chain (CLSC) network [9,10]. Recently, customers' behavior has remarkably changed. They are concerned about the environmental aspects of the products and their price and quality [11,12].
Numerous companies are pursuing greenness and sustainability by employing inspection, recycling, and refurbishing activities to recover the used products collected from end customers to increase their return rates [13]. For instance, Xerox, Canon, Kodak, Dell, and Acer have been pushing green operations. For example, by the year 2012, Xerox had effectively cut their emissions by 42% and their energy consumption by 31% [14].
Transportation fleets that move goods between supply chain hubs are a significant contributor to rising environmental consequences. According to Golicic et al. [15], fleet management has a favorable influence on energy efficiency and reduces fleet-related environmental consequences.
Using a heterogeneous fleet comprised of vehicles of varied sizes, the amount of fuel required to serve a given request portfolio can be drastically reduced [16]. Moreover, it provides a better balance of economic and environmental sustainability simultaneously [17,18].
Government policymakers should incentivize businesses to collect used products from customer zones for reverse logistics activities that in turn consolidate greenness and sustainability [19]. Considering efficient presorting centers at the customer zones substantially separates the low-quality returned products at the early stage from the refurbishment stream, thereby causing a reduction in the induced cost and carbon emissions [20].
Most real-life optimization problems are multi-objective and, usually, are uncertain due to the inherent errors of parameter values [21]. These errors result from inappropriate measurement or estimation due to a lack of knowledge of parameter values [22], a highly dynamic environment, and the physical impossibility to implement a computing solution in real-life problems [23]. Ignoring these uncertainties in optimization problems leads to infeasible and sub-optimal solutions [24].
Stochastic optimization (SO) and robust optimization (RO) are well-known approaches for solving uncertain optimization problems. In the SO approach, the probability of uncertain parameters must be defined accurately, but precisely defining this probability is difficult [23,25]. In addition, this approach devastates the convexity property and alleviates the complexity of the original problem in the chance constraints. However, in the RO approach, the uncertainty of the parameter is predefined in an uncertainty set characterized by its shape and size, and this approach overcomes some drawbacks of the SO approach [26].
As far as we are aware, few studies [27][28][29] have focused on the RO-based multiobjective optimization of CLSC networks. In these investigations, the implemented uncertainty set is the box uncertainty set with a homogeneous transportation system and no presorting centers in client zones. Motivated by the above facts, we developed a bi-objective robust optimization model for the GCLSC by considering presorting and heterogeneous transportation systems. The logical interpretation of the proposed study's framework is depicted in Figure 1. The main contributions of this study are as follows: • Presenting a bi-objective mixed-integer linear programming (MILP) model for the GCLSC while considering presorting and a heterogeneous transportation system as well as uncertain cost, selling price, and carbon emissions uncertainties. • Presenting a robust counterpart model formulation of the bi-objective MILP model of the GCLSC network under uncertainties using the well-known RO approach based on the "Interval + Polyhedral" uncertainty set. • Performing intensive computational experiments using a lexicographic weighted Tchebycheff (LWT) approach to obtain the Pareto optimal solutions of the robust bi-objective model, in which the probability bounds on constraint satisfaction are not violated.
The remainder of this paper is organized as follows. Section 2 presents the relevant literature pertaining to the uncertain GCLSC model from the previous studies. Section 3 discusses the problem description and mathematical modeling of the bi-objective GCLSC model. Section 4 presents the robust counterpart formulation of the bi-objective GCLSC model based on the combined "Interval+Polyhedral" uncertainty set. Section 5 presents the solution approach, LWT approach, and the conducted computational experiments. Finally, the conclusion and possible recommendations for future work are presented in Section 6.

Relevant Literature
Recently, several studies have focused on the uncertain and multi-objective nature of CLSC optimization models. Uncertainty is inherent in the CLSC models due to the dynamic nature of parameter data and errors that occur when measuring or estimating the input data. While some of these studies have focused on the uncertain CLSC optimization problem with a single objective, which is their primary target of minimizing the total cost incurred in the CLSC network, the rest are concerned with the inherent nature of the uncertainties and conflicting objectives of CLSC optimization networks.
Two streams of literature are reviewed: the first deals with the design of uncertain CLSC with a single objective, whereas the second involves the design of uncertain CLSC with multiple objectives.

Uncertain Single Objective CLSC
The uncertain optimization problem of single objective CLSC models based on the RO technique has recently been developed and expanded; its efficiency and robustness have also been demonstrated. For designing the CLSC, Pishvaee et al. [26] presented an MILP model, which was further extended to tackle uncertain parameters, such as transportation costs, returned products from first-hand customers to collection centers, and the demand of second customers for recovered products using the RO approach; their study objective is a minimization of the total cost.
The RO model for the CLSC network based on the regret value to determine facilities' locations and quantity of flows between facilities in the network were studied by Wang et al. [30].
Almaraj and Trafalis [31] proposed an adjustable RO approach for designing an uncertain CLSC network based on the budget uncertainty set to model the dynamic behavior of the customer demand over the planning horizon to minimize the total cost of the CLSC network. Ramezani et al. [32] proposed a scenario-based optimization model using the min-max regret criterion with uncertain demand and the return rate for designing the CLSC to maximize the total profit using a scenario relaxation algorithm.
Some studies are concerned with real case studies of the uncertain CLSC. For instance, Gholizadeh et al. [33] suggested a robust CLSC optimization model for disposable appliances considering demand, transportation, and operational costs uncertainties using a robust scenario-based stochastic optimization to maximize the profit of the network. Hasani et al. [34] proposed a mixed-integer nonlinear programming (MINLP) model for designing a global CLSC network for a medical device manufacturer, considering uncertain purchasing costs and consumer demand, managed with an RO technique based on the budget uncertainty, "Bertsimas and Sim", set. Jabbarzadeh et al. [35] proposed an SO approach for designing a CLSC network of glass companies by considering disruption risks across different disruption scenarios using the Lagrangian relaxation algorithm.
Gao and Ryan [36] proposed a hybrid model for designing a GCLSC network that combines SO and RO approaches aimed at minimizing the total cost incurred in the network. Samuel et al. [20] presented a scenario-based robust model for designing a GCLSC network of electronic products with considering the uncertain quality of the returned product to minimize the total cost.

Uncertain Multi-Objective CLSC Network
Due to the realistic multi-objective and uncertain nature of the CLSC design and optimization problems, researchers have paid attention to designing robust models that hedge against uncertain realization and provide optimal solutions for the conflicting nature of objectives.
The uncertainty in a multi-objective CLSC network using the SO approach is usually characterized by a probability distribution on the parameters. Zhalechian et al. [37] proposed an MINLP model to minimize the total cost and environmental impacts and maximize the positive social impacts of designing a supply chain network using a hybrid meta-heuristic algorithm; they applied SO and modified the game theory to handle uncertainty.
Ghasemzadeh et al. [38] suggested an MILP model to develop a stochastic CLSC network aimed to minimize Eco-indicator 99 and maximize profit. The model implemented in a real-world case study of the tire manufacturing industry, a two-stage stochastic optimization approach was implemented to cope with uncertainties in their study.
Moheb-Alizadeh et al. [39] developed a stochastic integrated multi-objective MINLP model, in which the design of a sustainable CLSC network considers sustainability outcomes and the efficiency of facility resource utilization, using a Lagrangian relaxation algorithm. A chance-constrained programming model has been proposed for handling uncertainties in a GCLSC network design aimed to minimize the expected value and variance of the total of cost, CO 2 emissions are controlled by providing a novel chance constraint approach using the general algebraic modeling system (GAMS) software; four multi-objective decision-making methods were applied by Abad and Pasan [40].
The fuzzy robust optimization (FRO) approach was implemented for dealing with uncertain parameters when the available information is vague.
Soleimani et al. [41] developed a multi-objective model for designing a CLSC network using an FRO approach, and their objective was to maximize the total profit and the responsiveness of customer demand while minimizing lost working days due to occupational accidents. In addition, a genetic algorithm was implemented for solving their model. Nayeri et al. [42] proposed an MILP model for designing a sustainable CLSC network for a water tank industry using an FRO approach to handle uncertainties in order to minimize the cost and environmental impact and maximize the social impacts using the goal programming approach.
Hybridization between SO and FRO approaches was implemented by Yu and Solvang [43] for handling uncertainties in designing the CLSC network to minimize the total cost and carbon footprint induced by the network. Scenario-based robust optimization (SRO) was also implemented to handle the uncertain parameters of the design of CLSC networks.
Ruimin et al. [44] proposed a robust MINLP model to deal with GCLSC by considering two conflicting objectives simultaneously, the economic cost and environmental impact, using scenario-based optimization.
An uncertain multi-objective MILP model of the CLSC network design was studied to minimize the total costs, maximize the on-time delivery of the products purchased from suppliers, and maximize the quality of the produced products on the forward chain that can be recovered in the reverse supply chain using an SRO method [45]. Gholizadeh et al. [46] proposed a bi-objective model to minimize the environmental impact and maximize profit for a dairy CLSC using a heuristic approach and robust optimization scenario-based. Some studies have implemented an interval-based RO approach where the uncertain parameters are predefined within an uncertainty set which characterized by its shape and size [22]. Darestani and Hemmati [28] used the RO approach for designing a CLSC network for perishable products to minimize the environmental and cost objectives. Yavari and Geraeli [47] proposed an MILP model for a multi-period and multi-product GCLSC network of dairy products with limited life shelf under uncertainties of demands, rate of return, and the quality of returned products based on the RO approach using a heuristic method to minimize the total cost and the environmental pollutants.
Jiao et al. [48] proposed a data-driven model in CLSC to mitigate recovery uncertainty and greenhouse emissions (GHE) based on uncertain customer demand to minimize the total cost using the chance constraint method to control the GHE and define the perturbation of uncertain demand using the RO approach based on the "Bertsimas and Sim" uncertainty set.

Problem Statement and Mathematical Modeling
The GCLSC proposed in this study compromises a set of customer zones where the returned products are collected, a set of presorting centers located at these zones, a recycling center, and an inspection and refurbishment (IR) center with different capacity sizes: small, medium, and big. The structure of the proposed GCLSC network (Figure 2) is the same structure presented by Samuel et al. [20]. Transportation between the facilities is a heterogeneous transportation system comprising different vehicle fleets that are varied in capacity, cost, and carbon emissions. The returned products at customer zones are either subject to presorting or nonpresorting, depending on the quality of the returned products. For high or very low-quality products, it is obvious that they can be transferred directly to the IR and recycling centers, respectively. Otherwise, the returned products are presorted in the presorting center depending on the quality of these products and the efficiency of the presorting center. A fraction of these products are sent to the recycling center, and the rest are sent to the IR center.
To maximize the profit that results from the difference between the selling of refurbished products and the total cost incurred in the network and minimize the total carbon emissions due to processes and transportation activities in the network, we consider further assumptions as follows: • The refurbished products are in a condition as good as new ones. • Selling price, costs, and carbon emissions are subject to uncertainty with an unknown probability distribution. • The capacities, costs, and carbon emissions due to transportation modes are predefined.
Costs and emissions are based on the round trip distances between network facilities. • The locations of customers zones, the recycling center, and the IR center are predefined. • The capacities of the facilities are predefined and deterministic. • Only one transportation mode can be used between any two consecutive facilities in the GCLSC network. Based on the above assumptions and description, an MILP model of the proposed bi-objective GCLSC network and its nomenclature presented at the end of this paper.
On the one hand, the first objective function,g pro , of the proposed model is the maximization of the profit; it comprises four parts: the first part,g 1 , represents the selling price of the refurbished products at the IR center; the second part,g 2 , represents the fixed costs of opening facilities in the entire network; the third part,g 3 , represents the transportation costs due to transportation activities in the GCLSC network using heterogeneous modes of transportation of vehicles v; the last part,g 4 , represents the cost of processing activities in the network, including presorting, inspection, refurbishment, and recycling.
On the other hand, the second objective function,g emi , of the proposed model is the minimization of the carbon emissions due to the various activities in the GCLSC network; it comprises two parts: the first part, em tr , represents the emissions due to transportation activities using different modes of transportation vehicles in the network part of the objective function. The second part, em pr , is the emissions due to different processing activities in the network, such as presorting, inspection, refurbishing, and recycling.
Constraint (8) ensure the equilibrium at customer zone k; that is, the amount of returned product p equals the quantities sent to the recycling center for scrapping and IR center without presorting in addition to the amount sent for presorting centers.
∑ v∈V X ins pkv + X pre pk + ∑ v∈V X rec pkv = r pk ∀p, k Constraint (9) guarantees that the amount sent for presorting cannot exceed the amount of returned products.
∑ p∈P X pre pk ≤ U k ∑ p∈P r pk ∀k (9) Constraint (10) governs the quantity sent from the presorting center to the IR center, which comprises good-quality products and a fraction of bad-quality products due to the inefficiency of the presorting process.
∑ v∈V Y ins pkv = X pre pk q pk + 1 − q pk β ∀p, k Constraint (11) guarantees that the amount sent from the customer zone k to the IR center of size s without presorting cannot exceed the maximum capacity of the biggest size of the IR center.
Constraint (12) governs the quantity sent from the presorting center to the recycling center for scrapping.
Constraint (13) ensures that the presorted quantity is presorted is directed to the IR center or the recycling center for each product p and customer zone k.
Constraint (14) governs the quantity of product p sent to the IR center from customer zone k, considering both presorted and non-presorted products. Constraint (15) guarantees the equilibrium at the IR center that the quantity of inspected product p is directed either for refurbishing or scrapping.
Constraint (16) determines the quantity of a product p directed from the IR center to the recycling center; this quantity comprises a fraction due to the presorting center inefficiency, and the rest due to the bad quality of non-presorted products.
Constraint (17) stipulates that the quantity sent from the customer zones in the presence and absence of presorting in addition to the IR center to the recycling center cannot exceed its capacity.
Constraint (18) Set of constraints with respect to multimode transportation vehicles v, constraints (20)- (25) ensure that only one type of vehicle is used to transport the products between every two facilities in the proposed GCLSC network.
The set of constraints (26)- (31) guarantee that the quantity of transported products between two facilities cannot exceed the capacity of the total number of the vehicles of selected transportation mode v.
∑ p∈P w p X rec pkv ≤ ca vNkv ∀k, v The set of constraints (32)-(37) that ensure the total number of vehicles of transportation mode v if this mode already has been selected.
The set of constraint (38) ensures that the total number of transported quantities between the facilities using the selected vehicle modes v equals the total number of returned products p.

The Robust Model of Bi-Objective Green Closed-Loop Supply Chain
RO is designed to deal with a lack of information. The uncertain parameters in robust optimization are taken at their worst case values; therefore, the robust optimization approach results in a solution that is immunized against uncertainty [24].
A general uncertain MILP model is given by P : min where x j and y k represent continuous and integer decision variables, respectively,c j and e k are the uncertain coefficients of the objective function, andã ij ,b ik , andd i are the true values of the ith constraint coefficients that are subjected to uncertainty; the P model can be reformulated as follows: Without loss of generality, we consider the general ith constraint, where the true value of uncertain constraint parametersã ij ,b ik , andd i are a ij +â ij ξ ij , b ik +b ik ξ ik , and d i +d i ξ i0 , respectively; a ij , b ik , and d i are the nominal value of the parameters,â ij ,b ik , andd i are the perturbations of the parameters around their nominal values respectively, and ξ ij , ξ ik , and ξ i0 are random independent variables that take values of [−1, 1]. Hence, the ith constraint can be reformulated as follows: J i and K i are the sets that contain the uncertain parameters in the ith constraint. For simplicity, we consider ξ = {ξ ij , ξ ik , ξ i0 } and U as the predefined set of random variables.
To become immune against infeasibility that would result from any realization of uncertainty, the ith constraint can be rewritten as follows: The robust counterpart optimization formulation of the ith constraint based on the combined interval and polyhedral, "Interval + polyhedral" uncertainty set is as follows: where θ ij , η ik , ϕ i0 , and ω i are positive dual variables, and Γ represents the adjustable size parameter of the "Interval + Polyhedral" uncertainty set which reflects the degree of conservativeness. Figure 3 illustrates different geometrical representations of the combined "Interval + Polyhedral" uncertainty set based on the value of the adjustable parameter Γ. The proof of robust counterpart formulation of the ith constraint, which is shown in Equation (47) is available in Li et al. [51] and Bertsimas and Sim [52]. Following the practical guide of the robust optimization approach proposed by Gorissen et al. [23] when the coefficients of the objective functions have uncertainty, constraints have been introduced to model the coefficients in the objective functions and approximate their values using auxiliary variables.
The robust counter reformulation of the P MOMILP based on the "Interval + Polyhedral" uncertainty set induced an MILP model, RC I−P − P MOMILP , as follow: ω 1 + θ pk ≥ŝ pk Z re f pkv ∀p ∈ P, k ∈ K, v ∈ V (50) ω 1 + θ pk pre ≥p pre pk X pre pk ∀p ∈ P, k ∈ K (60) ω 2 + η1 col p ≥ em col p r pk ∀p ∈ P, k ∈ K (73) ∀p ∈ P, k ∈ K (74) and set of constraints that have been depicted in Equations (8)- (42). The auxiliary variables, δ pro and δ emi , are introduced to approximate the profit and carbon emissions, respectively. A set of constraints are presented in Equations (49)-(65) to model the robust counterpart of the first uncertain objective function (Maxg pro ).
A set of constraints are presented in Equations (66)-(79) to model the robust counterpart of the second uncertain objective function (Ming emi ).

Probability Upper Bound of Constraint Violation in the Robust Optimization
In contrast to stochastic programming, RO does not require a known probability distribution to handle uncertain parameters. Probabilistic guarantees can be used to determine the lower bound on constraint satisfaction based on the desired constraint violation.
The probability upper bound of the constraint violation of the robust counterpart formulation based on the "Interval + Polyhedral" uncertainty set is proposed by Bertsimas and Sim [52] and Li et al. [53]. It is bounded B(|Q i |, Γ), where |Q i | is the cardinality of the uncertain parameters in the ith constraint. Assuming the probability distribution of uncertainty is symmetric, independent, and bounded, it is given as follows: The proof of upper bounds on the probability of constraint violation is provided by Bertsimas and Sim [52] and Li et al. [53].

Solution Approach
. . , f m (x)}. The general formulation of MOM is given below: Definition 1. A decision vector x * ∈ S is called Pareto optimal, or efficient, if there is no other decision vector x ∈ S such that x dominates x * . If x * ∈ S is a Pareto optimal solution, the vector f (x * ) is said to be a non-dominated point in the objective space [54]. For instance, suppose that S ⊂ R 3 and its feasible region f ⊂ R 2 , the continuous thick line contains all Pareto optimal solutions and f (x * ) as one of these solutions (Figure 4).

Definition 2.
A decision vector x * ∈ S is a weakly Pareto optimal, or weakly efficient, if there is no feasible solution x ∈ S such that f (x) < f (x * ) is a weakly efficient solution, and the vector f (x * ) is said to be a weakly non-dominated point in the objective space.

Definition 3.
A Pareto optimum solution is said to be supported if there are positive weights λ 1 , λ 2 , . . . , λ k that make the solution optimal in terms of the linear combination (weighted sums x ∈ S with coefficients λ 1 , λ 2 , . . . , λ k . If this is not the case, the solution is referred to be non-supported [55].
The solution approach implemented in this study is the LWT method which was developed by Steuer and Choo [56] to obtain Pareto optimal solutions by minimizing the maximum deviation of the objective functions from their values where they are optimized individually. The LWT approach outperforms existing techniques because of the following advantages [57,58]: • Achieves adequate optimal solutions of MOMs. • Provides both supported and non-supported Pareto solutions, whereas the weighted sum method, which is widely utilized in the literature, can only achieve supported Pareto solutions. • Obtains Pareto optimal solutions in both the convex and concave space hulls.
The mathematical formulation of the LWT approach for MOM is given below: The mathematical formulation of LWT approach for the robust counterpart model, RC I−P − P MOMILP , of the bi-objective GCLSC based on the "Interval + Polyhedral" uncertainty set is as follows: LWT − (RC I−P − P MOMILP ) : Min α and set of constraints that have been depicted in Equations (8)- (42) and Equations (49)-(79).

Computational Experiments
In this section, to evaluate the performance of the proposed robust model of the GCLSC network by considering a presorting and heterogeneous transportation system using the LWT approach discussed in Section 4.2, a set of computational experiments was conducted using GAMS/CPLEX optimization solver. Most of the model input parameters were adopted from the literature [20], which are based on realistic data. Vehicle capacity and CO 2 emissions were calculated for road transport operations as in Mohammed et al. [49], and the transportation cost based on vehicle type was calculated as in Wangsa and Wee [59] as shown in Table 2. The number of returned products, selling price of refurbished products at customer zones, operations costs, and weight of returned products are presented in Table 3. Capacities and fixed opening costs of different sizes of the IR center are presented in Table 4. Distances between customer zones, the IR center, and the recycling center are presented in Table 5, and the rest of the model parameters are presented in Table 6.   For model verification and to show the robustness of the RO approach, we consider that the perturbation of uncertain parameters are 5% and 10% around their nominal values assuming that the upper bounds of constraint violation are 10%, 15%, and 20%. In other words, the lower bounds of constraint satisfaction are 90%, 85%, and 80%, respectively. Five different weight sets (λ pro and λ emi ) were considered for both the objective functions.
Increasing the lower bound of the probability of constraint satisfaction results in an increase in the conservative of the obtained solutions; that is, increasing the size of the uncertainty set Γ increases the costs that result in a decrease in the profit from the economic aspect. However, from the environmental perspective, increasing the uncertainty set leads to more carbon emission trade-offs between the conservatives and the values of profit and carbon emissions.
Robust objective function values were achieved at lower levels than the deterministic model. The values of the objective functions worsened as the perturbation around the nominal values increased (Tables 7-12). On the one hand, when the perturbation was increased from 5% to 10%, for example, at a lower bound of the probability of constraint satisfaction of 80% and weights of λ pro = λ emi = 0.5, the profit declined by 6.62% and 13.42%, respectively, from their deterministic values. The value of carbon emissions, on the other hand, decreased by 1.10% and 2.29% from their deterministic levels, as shown in Tables 7 and 8, respectively. For the lower bound of constraint satisfaction at 80% with 5% data perturbations, the deviation of the robust profit decreased from its deterministic value by 5.24% to 7.29%. Meanwhile, the robust values of the carbon emission increase ranging from 1.09% to 1.52% (Table 7). By increasing the data perturbations to 10% the robust values for the profit and carbon emissions deviated by 10.47% to 14.59% and 2.12% to 3.07%, respectively (Table 8).
Increasing the lower bound of constraints satisfaction to 85% with 5% perturbations resulted in an increase in the deviation, ranging from 6.03% to 7.86% and 1.09% to 1.62% for the profit and carbon emissions, respectively, (Table 9). However, considering 10% perturbation, the deviation ranged from 12.07% to 15.42% and 2.12% to 3.26% for the profit and carbon emissions, respectively (Table 10).
For the lower bound of the probability of constraints satisfaction at 90% with 5% perturbation of the uncertain parameters, the deviation ranged from 6.57% to 8.40% and 1.21% to 1.74% for the profit and carbon emissions respectively (Table 11). In addition, considering 10% data perturbations, the deviation ranged from 13.15% to 16.25% and 2.61% to 3.52% for the profit and carbon emissions, respectively (Table 12).
Moreover, the robust optimal Pareto front of the bi-objective GCLSC model with 90% probability of constraints satisfaction with 5% and 10% is presented in Figure 5a,b, respectively. By increasing the deviation from the nominal values, the optimal values of the objective functions become worse as shown in the Pareto front representation in Figure 5.   On one hand, the obtained robust Pareto front of the proposed bi-objective GCLSC model is achieved with the levels of 80%, 85%, and 90% probability of constraints satisfaction and 5% and 10% of data perturbations, revealing that our model yields more conservative solutions that are immune against the uncertainty and comply with the nature of conflicting bi-objective functions.
On the other hand, increasing the probability of constraints satisfaction increases the uncertain space beside the integer variables of the model, hence increasing the model's complexity.

Conclusions and Outlook
This study provided a robust bi-objective optimization model for the GCLSC network by considering a presorting and heterogeneous transportation system. The proposed biobjective MILP model aims to maximize the profit and minimize the total carbon emissions incurred in the entire network by considering the selling price, costs, and carbon emissions uncertainties. A robust optimization approach based on the "Interval + Polyhedral" uncertainty set was adopted to tackle these uncertainties. An LWT approach was implemented to obtain the Pareto optimal solutions, which were hedges against perturbations of the uncertain parameters. A set of computational experiments was conducted to validate and show the performance of the proposed model.
Pareto optimal fronts of the robust bi-objective GCLSC model were examined considering the maximum probability of constraint violation. In other words, the lower bound of constraint satisfaction was considered at three levels: 80%, 85%, and 90% combined with 5% and 10% perturbations of the uncertain parameters around their nominal values. The lower the bound of constraint satisfaction, the more conservative as well as the less profit and more carbon emissions. Comparing the obtained robust solution with the deterministic solution in the experiments showed that our proposed model is more conservative and robust. Finally, the limitations of the proposed study, which can provide possible recommendations for future research, are as follows: the amount of time needed for computing is directly proportional to the number of integer variables; in other words, the amount of time needed for computation increases as the number of integer variables increases [60]. Small size instances are used in this paper for model validation. For large-scale problems, the model becomes computationally hard to obtain the optimal solution using the GAMS solver. Therefore, it is one of our future extension directions to propose heuristic or metaheuristic approaches to handle the large-size instances of the developed GCLSC.
However, future studies should extend the proposed robust bi-objective model to a multi-objective model by considering a social objective that contributes to the welfare and prosperity of society, such as maximizing the number of job opportunities and fair salaries for employees. Enhancement can be made using the quality of the robust solutions, the iterative solution framework proposed by Li and Floudas [61], and different uncertainty sets shapes such as ellipsoidal, polyhedral, combined "ellipsoidal and polyhedral", and combined "box + ellipsoidal + polyhedral".