A two-echelon location routing problem considering sustainability and hybrid open and closed routes under uncertainty

Location-routing is an extremely important problem in supply chain management. In the location-routing problem, decisions are made about the location of facilities such as distribution centers as well as the set of vehicle routes. Today, organizations seek to reduce the transportation cost by outsourcing leading to a particular kind of transportation problems known as open routing. However, the increasing attention to environment have led to paying more attention to environmental issues and reducing the environmental impacts of logistics activities. To this end, in this paper, both open and closed routes were simultaneously addressed by developing a multi-objective mixed integer linear programming model that included three economic, environmental, and social responsibility aspects. The three objective functions of the proposed model encompass the minimization of total costs and greenhouse gas emissions, and the maximization of employment rate and economic development. Also, in this study, a different type of routing was considered in each echelon. A small-sized problem instance was solved using the Augmented Epsilon Constraint (AEC) method with the CPLEX Optimizer Solver for the validation of the proposed model. Moreover, the sensitivity analysis was performed to investigate the effect of changing main parameters on the values of the objective function. Due to the NP-Hardness of the problem, two efficient metaheuristic algorithms of Non-dominated Sorting Genetic Algorithm (NSGA-II) and Multi-Objective Stochastic Fractal Search (MOSFS) were exploited to solve the medium and large size problems. The performance of the algorithms was compared on the basis of six different well-known indexes of Time, MID, RAS, Diversity, Spacing, and SNS. According to the obtained results, the performance of the MOSFS algorithm was %20, %9, %11.22, %10.03, and %19.06 higher than the performance of the NSGA-II on the basis of SNS, RAS, MID, Diversity, and Time indexes, respectively. On the other hand, the NSGA-II performance was %6.3 higher than the MOSFS performance in terms of Spacing index.


Introduction
To take advantage of competitive and economic advantages, it is essential for the organization to focus on minimizing costs or maximizing profits [1]. On the other hand, the production activities have detrimental environmental impacts due to the overconsumption of natural resources [2]. Road transportation also faces remarkable environmental problems and issues such as the reduction of fossil fuels usage and greenhouse gases emission. To resolve these issues, governments and organizations must take practical initiatives and change transportation policies to provide sustainable logistics operations [3]. The environmental impacts have been considered as one of the primary aspects of sustainable development due to raising awareness about the pollutants and their negative impacts on humans' lives [4].
Sustainable development also addresses social responsibility issues [5]. Sustainable development is an initial strategic goal for the supply chain [6]. Most of the previous studies on sustainability problems have investigated the economic and environmental aspects. However, limited studies have addressed the social responsibility aspect. This lack of research is due to the complexity of social responsibility modeling [7].
Sustainable supply chain problems are classified into two main categories: green and sustainable supply chain [8] so that the green supply chain (the first category) is the subset of sustainable supply chain (the second category) [9]. A supply chain consists of suppliers, manufacturing plants, and customers as well as materials flow, information flow and cash flow between the components. The management of a supply chain system can be divided into three levels: strategic, technical, and operational levels. An effective cooperation among the entire supply chain system can lead to a significant reduction in total costs. Supply chain network design usually initiates with selecting potential locations and determining the capacity required for factories [10]. Leng et al. proposed a bi-objective model for the sustainable location-routing problem and solved it using well-known multi-objective evolutionary algorithms [11]. In addition, Galindres et al. [12] considered a multi-objective sustainable capacitated location routing problem and solved it using exact and approximate solution methods.
A considerable amount of investment is required for setting up new facilities (factories, depots, or warehouses). However, these facilities can be utilized for a lengthy duration [13]. The performance of the entire supply chain network is considerably influenced by the locations of factories. Location-routing problems (LRP) integrate two basic problems in logistics. In LRP, the decisions on the locations of various facilities are integrated with the decisions on the vehicles' routes. Independent decisions on either of these two categories can significantly affect the global optimal solution [14].
Delivery of goods from origin to destination is often done through one or more intermediate facilities such as hubs and warehouses. These types of distribution systems are commonly known as multi-echelon systems. Each echelon can be defined as the connector between two adjacent levels. Here, an echelon may be any kind of facility. Multi-echelon distribution systems are often used by public and private sector organizations in distribution networks to implement transportation systems and traffic planning strategies. The common examples of multi-echelon distribution systems are logistics systems and multi-mode urban transportation systems. Twoechelon distribution systems are a special type of multi-echelon distribution systems in which the distribution network consists of three levels. Two-echelon distribution systems have received more attention in recent years because of their wide application in daily jobs.
A distribution network that consists of three distinct sets of vertices corresponding to potential factory locations (origins), potential secondary (mid-level) facility locations and destinations (customers) is known as two-echelon distribution network, in which the locations of customers are constant and predetermined, however, the locations of the required plants and/or the peripheral and secondary provisions and facilities are not predefined. Fig. 1 depicts an overview of a two-echelon location-routing problem (2E-LRP).
Nowadays, transportation companies often outsource transportation in one or more echelons to decrease the transportation cost and consequently the cost of the final product. This matter has led to the emerge of a new type of vehicle routing problems called open routing (The vehicles do not return to the origin after serving the last customer and an open route is created). In practice, the open vehicle routing problem (OVRP) relates to a situation in which a company does not possess its own fleet of vehicles or its current fleet of vehicles is not able to meet the demands of all customers. In this case, all or a part of the distribution activities are contractually outsourced to a third-party logistics company [15].
Nevertheless, the reduction of the adverse and negative effects of transporting operations and activities on the environment is a critical issue. Consequently, green transportation systems have been developed in manufacturing and distribution industry. Different approaches such as Pollution Routing Problem (PRP) may be employed in the vehicle routing problems for reducing the greenhouse gases emissions. The goal of PRP is to select the vehicle routes with the lowest greenhouse gas emissions, especially carbon dioxide emissions [7].
According to the aforementioned issues, this research seeks to present a sustainable three-objective model for 2E-LRP under uncertainty. The problem consists of different three levels including factories, depots, and customers. In addition to the minimization of the total costs of the whole distribution system, the amount of CO 2 emissions is minimized, which is one of the major criteria and factors for sustainable development. Also, the social responsibility aspect (employment rate and the community development) is considered in the proposed model. In this problem, the location of depots is determined and the routing of the vehicles is distinguished on both two echelons. In addition, both open and closed vehicle routing problems are taken into account. The proposed model is called the hybrid open and closed sustainable two-echelon location-routing problem (COM-S2ELRP). The results of solving this multiobjective mathematical programming model determine the optimal location of the depots, the routing of the vehicles, and the optimal amount of transportation in both two echelons, in addition to the allocation of customers and depots.

Literature review
Routes are expanded from the initial facilities to the depots and from the depots to the customers to hand over goods to the customers. In this type of LRP, routes from the initial facilities to the depots are called the first echelon routes, and routes from the depots to the customers are called second the echelon routes [16]. Jacobsen & Madsen investigated 2E-LRP for the first time considering the newspaper distribution system in Denmark [17]. Lin and Lei [18] developed a 2E-LRP model that included a number of depots as well as two kinds of customers. Crainic et al. [19] solved 2E-LRP using a heuristic algorithm. Nguyen et al. [20] suggested a model for 2E-LRP presuming a focal depot and various potential locations for other depots. Also, the set-up costs of depots were different and their capacities were considered limited.
Martínez-Salazar et al. [21] presented a bi-objective 2E-LRP model including cost and traveled distance minimization and solved it using two metaheuristic algorithms. Also, Rahmani et al. [22] developed a mixed integer programming (MIP) model with several assumptions for 2E-LRP and solved it employing two heuristic approaches. Moreover, Vidović et al. [23] provided a MIP model for 2E-LRP to design a collecting and recycling system for the non-hazardous recyclable waste. Ouhader and El kyal [24] suggested a sustainable 2E-LRP for shipping products from suppliers to customers. The model included three objective functions: minimization of emissions and costs as well as maximization of job opportunities. Finally, the Epsilon Constraint method (ECM) was exploited as a solution method. Zhao et al. [23] considered a heterogeneous transportation fleet for 2E-LRP and exploited a heuristic approach to tackle it.
Pichka et al. [25] introduced a new type of 2E-LRP known as the two-echelon open location-routing problem (2E-OLRP) presuming the vehicles do not turn back to the depots. Darvish et al. [26] proposed a flexible 2E-LRP considering flexible delivery time and distribution network design. Amiri et al. [27] considered time window for 2E-LRP and utilized the Lagrangian Relaxation method as the solution approach.
Cao et al. [28] reviewed the studies on 2E-LRP which shows there is no research on open and closed routing problems. Lu et al. [29] proposed a 2E-LRP model considering closed routing and applied two metaheuristic algorithms for solving the model. Liu and Jiang [30] developed an MIP model for the combination of open and closed vehicle routing problem (COMVRP) and employed a metaheuristic algorithm to deal with it. Yu and Lin [31] applied the Simulated Annealing (SA) algorithm to solve the open vehicle routing problem.
Several studies conducted on the sustainable vehicle routing problems have addressed economic, environmental, and social responsibility aspects. For example, Navazi et al. [32] incorporated sustainability into location-routing problem considering economic (cost minimization), environmental (emission minimization) and social responsibility (customer satisfaction) aspects for collecting expired products. They solved the problem using NSGA-II and MOPSO. Also, Navazi et al. [33] proposed a sustainable location-routing model for product delivery to customers. In this study, four goals of reducing costs, environmental pollutions, driving accidents and increasing customer satisfaction through timely delivery were considered. Tayebi Araghi et al. [34] developed a model for the green multi-facilities open location-routing problem with planar facility locations and uncertain customer. In this model, the random location of facilities was included as an effective factor in supply chain costs. Also, the location of the depots, the allocation of vehicles and the selection of the routes were taken into account to minimize the emission of carbon dioxide in the entire supply chain. Nucamendi-Guillén et al. [35] suggested a mixed integer linear programming (MILP) model for the multi-depot open location-routing problem with a heterogeneous fixed fleet (MD-OLRP) to solve the problem of a company in collecting raw consumable materials. In addition, Momeni et al. [36] used vehicle routing in environmental protection problem for the first time.
Liu and Liu [37] suggested a sustainable two-stage stochastic model. Zhang et al. [38] provided a sustainable location-routing model considering multiple depots in emergency conditions to minimize relief cost, traveling time, and CO 2 emission. Nekooghadirli et al. [39] presented a bi-objective model for the location-routing and inventory problem considering uncertain traveling time and customer demand. They also utilized several metaheuristic algorithms to deal with the problem. Tang et al. [40] introduced a model for the sustainable location-routing and inventory problem. In this study, environmental issues based on consumer behavior were investigated. Ebrahimi [41] proposed a model with multiple objective functions for random allocation routing problem taking the concepts of discount and sustainability into account. The objectives comprise economic (minimizing costs), environmental (minimizing pollutions), and social responsibility (maximizing responsiveness to customers) aspects. Finally, the model was solved using ECM.
Wang et al. [42] modelled 2E-LRP considering multiple periods and shared transportation resource to find optimal facility locations as well as optimal routes in different periods during the decision horizon. Gandra et al. [43] investigated the impact of loading restrictions on the two-echelon location-routing problem. Fallahtafti et al. [44] presented a multi-objective two-echelon location-routing problem for cash logistics to reduce the risk of theft while transporting cash and solved the problem using meta-heuristic approaches. Bassey and Zelibe [45] suggested a mixed integer nonlinear programming model for the problem of two-echelon inventory location model with response time requirement and lateral transshipment. Cheng et al. [46] proposed a multi-period two-echelon location-routing model (MP-2ELRP) to minimize the cost and duration of cleaning up the waste caused by accidents using Temporary Disaster Waste Management Sites (TDWMSs). In these sites, waste was stored and processed before being sent to disposal sites. Liu et al. [47] examined the sustainable location, routing, and inventory problem. Moreover, Babaee et al. [48] addressed the sustainable location routing problem. Furthermore, Masoudipour et al. [49] presented a model for the sustainable closed-loop supply chain problem. Du et al. [50] presented a multi-objective optimization model for two-echelon joint delivery location-routing problem considering carbon emissions and operational costs under online shopping. Mahmoodirad et al. [51] investigated the sustainable multi-objective multi-product location-routing-inventory problem considering open routing and direct transportation. They also employed several metaheuristic algorithms such as NSGA-II and MOGWO. Zandkarimkhani et al. [52] developed a model for designing a sustainable open loop supply chain network considering the location-routing problem with a combined approach based on Fuzzy AHP and Fuzzy TOPSIS methods. Raeisi and Jafarzadeh Ghoushchi [53] developed a robust fuzzy multi-objective location-routing model for hazardous waste problem under uncertain conditions to overcome waste disposal and prevent the spread of COVID-19. They   Fig. 2. An example of a feasible solution to an instance problem. solved the model using meta-heuristic algorithms and compared the solutions. In addition, Ben Mohamed et al. [54] suggested a model for the two-echelon stochastic multi-period capacitated location-routing problem (2E-SM-CLRP) considering stochastic and time-varying demand and varying costs.
The review of studies on the two-echelon location-routing problems (2E-LRP) shows that open and closed routes have not been considered simultaneously in previous studies. In this paper, a sustainable mixed integer linear programming (MILP) model including economic aspect (cost minimization), environmental aspect (minimization of transportation emissions) and social responsibility aspect (job creation and community development) is proposed for the problem of two-echelon location-routing considering the hybrid open and closed routes under uncertainty conditions. In addition, in this research, a distinct kind of routing is examined in each echelon.

The proposed mathematical model
As stated in the previous section, both open and closed routes have not been simultaneously investigated so far. In the present research, a sustainable MILP model consisting all three aspects of economic (cost minimization), environmental (minimization of transportation emissions) and social responsibility (job creation and community development) is rendered for 2E-LRP taking the hybrid open and closed routes under uncertainty conditions into account. Also, a different type of routing is considered in each echelon.
This section proposes a three-objective MILP model for COM-S2ELRP under uncertainty conditions. In the proposed model, the minimization of total costs and greenhouse gas emissions together with the maximization of employment rate and community development through establishing the facilities in deprived and undeveloped regions is considered.
The real-world assumptions are incorporated into this model in such a way that the closed routes are assumed for the first echelon (each trip should finish at the beginning point that is a factory) and the open routes are assumed for the second echelon (succeeding the satisfaction of the last customer's demand, the vehicle do not turn back to the depot, in other words, the fleet is outsourced based on a contract for transporting goods). Fig. 2 shows an overview of a feasible solution to an instance single-period single-product problem. As shown in the above figure, decision on the facility location is made only at the second level, and at the first level decision on the location is made, while decision on routing is made in both echelons.
The sets, indices, parameters, and decision variables of the proposed three-objective MILP model are presented as follows:

Sets and indices
The set of the whole network notes Np The set of factories at the first level of the distribution network N d The second level facility set (depots, warehouses or distribution centers) Nc The set of customers at the third level A The set of all connecting arcs between network nodes V The set of available vehicles in each facility (factory) of the first level of the distribution network i, j, k, m The index of each network node v The vehicle index in the first echelon The fixed cost of using each vehicle in the first echelon FV2

Parameters related to the economic aspect
The fixed cost of using each vehicle in the second echelon di The demand of customer i cij Direct (Euclidean) distance (cost corresponding to transport) between two nodes i and j ∈ N NV Maximum number of available vehicles in each depot CP k The capacity of the first level facility (factory) CD k The capacity of the second level facility (depot) CV v 1 The capacity limit of each vehicle in the first echelon CV2 The capacity limit of each vehicle in the second echelon

Parameters related to the environmental aspect
The amount of CO 2 emissions of the empty vehicles of the first echelon (kg CO 2 /km) Eve2 The amount of CO 2 emissions of the empty vehicles of the second echelon (kg CO 2 /km) The amount of CO 2 emissions of the full vehicles of the first echelon (kg CO 2 /km)

E vf2
The amount of CO 2 emissions of the full vehicles of the second echelon (kg CO 2 /km)

Parameters related to the social responsibility aspect ωem
The importance factor of employment The importance factor of community development JO k Number of permanent jobs created corresponding to each depot in node (region) k EV k The rate of increase in the economic value of node k corresponding to a depot set-up in node (region) k ur k The unemployment rate of node (region) k rd k The level of economic development in node (region) k

Decision variables
Binary variable corresponding to depot set-up (1: if depot k is set up, otherwise: 0) s v ij Binary variable corresponding to routing in the first echelon (1: if vehicle v is traveled between nodes i and j in the first echelon, otherwise: 0) xij Binary variable corresponding to routing in the second echelon (1: if a route is created between nodes i and j in the second echelon, otherwise: 0) w ik Binary variable corresponding to the assignment of customers to depots (1: if customer i is assigned to depot k, otherwise: 0) The integer variable corresponding to the amount of remaining load in vehicle v in the first echelon traveling from node i to node j Uij The integer variable corresponding to the amount of remaining load in the vehicles of the second echelon traveling from node i to node j Consider a complete directional graph network G = (N,A). In which, N = N p ∪ N d ∪ N c represents the set of all network nodes. N p , N d , and N c denote factories (the first level), potential locations for depots (the second level), and customers (the third level), respectively. Also, A = {(i, j) : i, jεN, i ∕ = j} denotes the set of all network arcs. In addition, the following relation exists for all i, j, and d ∈ N: c id + c dj ≥ c ij , in which, c ij represents the cost of traveling between nodes i and j. Every customer has a demand equal to the amount of

Three-objective mixed integer linear programming (MILP) model
Now, the three-objective MILP model is presented for the COM-S2ELRP problem in hand: Equation (1) minimizes the total cost of the entire system and comprises five components: (1) The transportation costs of the first echelon vehicles, (2) The transportation costs of the second echelon vehicles, transportation costs of second-class travel by secondary vehicles, the fixed costs of setting up depots in the second level, the fixed costs of employing vehicles in the first echelon, and the fixed costs of utilizing vehicles in the second echelon.
Equation (2) minimizes the total amount of CO 2 emissions of the whole system and consists of two components: the amount of CO 2 emissions of the first echelon, and the amount of CO 2 emissions of the second echelon. It is noteworthy that the amount of CO 2 emissions in both two echelons depends on the type of vehicle and its load and distance traveled, also, the vehicles are heterogeneous in the first echelon and homogeneous in the second echelon.
The third objective function deals with social responsibility. According to ISO 26000 guidance on social responsibility, community participation and development is one of the main aspects of social responsibility. There are two main methods to calculate the social responsibility aspect: (1) Employment, (2) Regional development In the proposed model, setting up a depot as a facilitator of intermediate transfer can lead to the participation and development of the community. In fact, employment and economic development are the two main reasons for establishing facilities. In other words, job opportunities resulting from a depot set-up are considered. Social value is calculated using the number of created job opportunities and the unemployment rate in that region. In this regard, creating job opportunities in a region with a higher unemployment rate leads to greater social value. Regional development rate and economic value of established depots are addressed to assess social value. Hence, higher social value means more importance given to the less developed regions. The economic development and employment criteria are presented as follows: Both above criteria are taken into account in the third objective function of the model shown in Equation (3): Subject to: Constraints related to the first echelon: Equation (4) ensures that the N p number of factories must be established at the first level.
Constraint (5) states that each of the available vehicles in each factory cannot travel more than one route. ∑ Constraints (6), (7), and (8) ensure that for each node in the first and second levels corresponding to the first echelon routes, exactly one arc with one vehicle enters and one arc with the same vehicle exits.
Constraint (9) ensures that the first echelon routes originating from the existing factories go to a depot that already has been set up.
Constraint (10) shows that the total demand met by the factories in the first echelon cannot exceed their capacity.
Constraint (11) hinders the formation of a route between factories in the first echelon.
Constraint (12) refers to two basic principles in the first echelon: first, that it guarantees the balance of flow in the routes created in the first echelon, and second, that this constraint hinders the formation of sub-loops in the first echelon.
Constraint (13) ensures that the total load of a vehicle in the first echelon cannot exceed the capacity of that vehicle.
Constraint (14) mandates that the first echelon vehicle returning to the origin factory at the end of the tour (closed route) cannot have any load and must be empty.
Constraints related to the second echelon: Constraint (15) indicates that the second echelon routes can originate from a depot that already has been set up.
Constraint (16) ensures that each customer receives service exactly once.
∑ jεN/Np Constraint (17) ensures that for each node in the second echelon the number of incoming arcs (routes) is equal to the number of out coming arcs (routes). ∑ jεN/Np Constraint (18) denotes the flow balance in the second echelon.
Constraint (19) guarantees that the remaining demand of the customers (the remaining load of the vehicle in the second echelon) cannot exceed the capacity of the vehicle. ∑ Constraint (20) mandates that the total demand of customers assigned to a particular depot must be met by vehicles sent from the same depot.
Constraint (21) ensures that in the second echelon the amount of remaining load in the vehicle after serving the last customer is equal to zero (open route).
Constraints (22) and (23) Constraint (24) mandates that each customer is assigned to only one depot. ∑ i∈Nc d i w ik ≤ CD k y k ∀k ∈ N d (25) Constraint (25) ensures that the total demand of customers assigned to a depot cannot exceed the capacity of that depot.

The proposed robust possibilistic programming model
Because of the insufficiency of information and uncertainty associated with the model parameters in the previous section, determining the exact values of the parameters will be very difficult and impossible in some case. The model presented in the previous section for the sustainable location-routing problem is a multi-objective mixed integer programming model in which the fuzzy numbers are utilized for the non-deterministic cost and capacity parameters. To deal with the uncertainty of the problem, a robust possibilistic programing model is developed. Possibilistic programming is utilized to tackle the uncertainty associated with the objective function coefficients and the model constraints. In possibilistic programming, the historical data are considered together with the decision-maker's opinion. In this study, a robust possibilistic programming model based on the chance-based possibilistic programming model presented by Pishvaee et al. [55].
A solution to an optimization problem is a robust solution if it encompasses feasibility robustness and optimality robustness. Feasibility robustness means that the solution must remain feasible for almost all possible scenarios of uncertain parameters. Optimality robustness means that the value of the objective function for a robust solution must have a minimum deviation from its optimal value for almost all amounts of uncertain parameters.

The proposed robust programming model
First, the necessity measure, that is applied to defuzzify the fuzzy numbers, is explained. If the uncertain model parameters are presumed based on the chance and trapezoid fuzzy numbers denoted by ω = (ω 1 , ω 2 , ω 3 , ω 4 ), the equivalent crisp number is obtained using Formula (36): The necessity of the constraint ω ≤ r is defined in Formula (37): Nec denotes the necessity measure in the chance-based constraint programming.
Based on the upper formula, if α ≥ 0.5: Relation (38) can be directly used to convert the fuzzy constraints into their equivalent crisp constraints. In this study, the trapezoidal fuzzy numbers are employed for the fixed and variable costs, shown in Fig. 3: According to Formula (36) and Relation (38), the uncertain mathematical programming model presented in the previous section can be defuzzified shown in Equation (39) and Constraints (40)(41)(42) :   Fig. 3. Fuzzy parameter with trapezoidal distribution.
Here, the converted parts corresponding to the first objective function and the constraints of the primary model that have fuzzy numbers are presented. The robust model considering feasibility robustness and optimality robustness is expressed using Equation (43) and Constraints (44)(45)(46): In the above model, the minimum confidence level of the fuzzy constraints (α, β, and γ) must be determined by the decision-maker in such a way that like the sensitivity analysis method, the amounts of the parameters must be changed. As the number of fuzzy constraints increases, the number of tests required to determine the appropriate values of confidence levels increases significantly. Also, the model is not sensitive to the deviation of the objective function value from its optimal value, and this can impose a high risk on the decision-maker in real-world problems. Therefore, the application of the robust possibilistic approach is effective to minimize the imposed risk [55].
The first expression E(Z 1 ) represents the expected value of the first objective function. The second expression, σ× (Z 1 Max − Z 1 Min ), represents the difference between the two boundary values of Z 1 Max and Z 1 Min , which are obtained by placing the upper bound amounts and the lower bound amounts of the parameters in the objective function Z 1 , respectively. σ indicates the importance of this expression compared to other expressions of the objective function. Therefore, the second expression results in minimizing the maximum negative and positive deviations from the expected optimal value of Z 1 . In addition, this expression controls the optimality robustness of the solution vector. The subsequent expressions specify the confidence level associated with every fuzzy constraint in which δ i is the possible deviation penalty of constraint i that contains the uncertain parameter. The last three expressions of the first objective function determine the difference between the value used in the chance-based constraint with the best value of the uncertain parameter. In addition, these last three expressions control the feasibility robustness.

Solving methodology
In this study, the Augmented Epsilon Constraint (AEC) method is used to solve the multi-objective optimization model and obtain the efficient Pareto frontier. Due to the NP-Hardness of the location-routing problem, two metaheuristic algorithms named NSGA-II and MOSFS are employed to deal with this problem.

The Augmented Epsilon Constraint (AEC) method
In the AEC method, first the appropriate range for changing the objective function values (e i ) must be determined and then the Pareto front is obtained for different values of (e i ). For better implementation of this method, the appropriate ranges of epsilons (e i ) should be determined using the Lex method. The two main steps in this method contain: determining the range of the e i values, and solving the AEC model. The AEC model of this research is presented as Model (47) [56].

Initial solution representation method
The continuous solution representation is employed for this problem. The solution representation contains a string of decimal numbers ranging from 0 to 1 with the length of Nc + Nd − 1 + V + Nd + V − 1. V represents the set of available vehicles in each facility, Nd denotes the second-echelon facility, and Nc denotes the total number of customers at the second echelon. For instance, a solution representation for a numerical example with V = 2, Nd = 3, Nc = 5 can be shown as follows: The permutation of the numbers of this part is obtained by arranging in a descending order. In this permutation, the numbers that are greater than Nc are considered as the separators, so that any series of numbers that are less than Vc , which are placed in sequence, corresponding with a route starting from a facility and ending at the last point of that group. In this example, the formed route initiates from the first facility and finishes after passing the route of 3 − 4 − 5 − 1.
The second part of the solution representation with the length V is corresponding to the allocation of the vehicles to the factories. Same as the first part, the routes linking factories and facilities are found after sorting in descending order using separators. 4 2 3 1 In this example, the route that can be formed starts from the second factory and passes through 2-3-1. It should be noted that facilities that are not used will be removed from this route. In this example, as only the first facility is utilized, the route originates from the second factory to the first facility and vice versa. So far, the factories and routes have been determined. Then, according to the route, the number of commodities and goods loaded at the starting of the route, the number of goods unloaded at each node and the traveled length can be calculated.

Non-dominated Sorting Genetic Algorithm (NSGA-II)
Genetic algorithm (GA) is one of the well-known metaheuristic algorithms that has emerged from biological models of living organisms. In this algorithm, the natural selection process is simulated so that the fittest and most desirable solutions are chosen to produce offspring of the next generation. The NSGA-II algorithm is a multi-objective GA that contains the following steps [57,58]: (1) Creating an initial population (2) Calculating the fitness criteria (3) Non-dominated sorting the population and calculating crowding distance Performing crossover and mutation operations to produce new offspring (4) Combining the initial population and the population created by the crossover and mutation operations (5) Replacing the initial (parent) population with the fittest (the most desirable) members of the combined population created in the previous step (6) All steps are repeated until the desirable generation (or optimal solutions) are achieved.

Crossover operator (One-point crossover)
In the first step, it randomly selects a pair of chromosomes (solutions). In the second step, it randomly selects a location along the chromosome string, and finally in the third step, it swaps the amounts of the two strings regarding the location selected in the previous step. It should be noted that the crossover rate must be adjusted appropriately to determine what percentage of the current population is selected for crossover operation. An example of the crossover operation is represented in Fig. 4.

Mutation operator
This operation is usually applied to a single chromosome (solution), which causes some of the chromosome genes to be randomly changed to produce a new offspring. This operator increases the scatter of solutions and reduces the probability of getting stuck in a local optimum. In this research, a gene is randomly selected and a new amount is randomly assigned to that given gene. An example of the mutation operation is shown in Fig. 5.

Multi-Objective Stochastic Fractal Search (MOSFS)
The Stochastic Fractal Search (SFS) algorithm is an efficient metaheuristic algorithm introduced by Salimi [59]. This algorithm has been used in various research topics. The Multi-Objective version of this algorithm (MOSFS) is able to solve complex multi-objective optimization problems by obtaining an efficacious set of the best non-dominated solutions with increased diversity. This algorithm is also capable of efficiently searching the solution space. Moreover, exploration and exploitation features are ensured by systematic random walks together with adaptive jump distance. Furthermore, the MOSFS algorithms contains fewer parameters to be tuned. For more information about the characteristics and advantages of the MOSFS algorithms, please refer to Ref. [60].

Problem design
A set of Sterle standard problems are used to test and evaluate the proposed model [61]. These problems are known according to the Fig. 4. A crossover representation for the chromosome of the sample problem.
M. Hajghani et al. location of facilities and customers in three groups including I1, I2 and I3, and each group is known by 31 problems, represented in Table 1.
In this paper, the group of I2 problems with the following specifications is considered, shown in Table 2.

Adjusting the parameters of the algorithms using the Taguchi method
The performance of any metaheuristic algorithm is strongly affected by the amount set for its parameters.

Parameter setting for the NSGA-II and MOSFS algorithms
Four NSGA-II parameters including MaxIt (maximum iteration), NPOP (number of initial population), PC (crossover rate), and PM (mutation rate) must be set at their optimal levels. Three MOSFS parameters including MaxIt (maximum iteration), Diff (number of initial population), and Walk (harmonic memory coefficient) must be set at their optimal levels. In this section, three levels are considered for each parameter. For this purpose, problem number 15, which encompasses 3 factories, 10 depots, and 20 customers, is chosen. For each parameter, three levels of low (1), medium (2) and high (3) are defined separately to solve the problem, shown in Table 3.
After determining the levels for the algorithm parameters using the Taguchi method and the MINITAB software, the required experiments are designed, shown in Table 4. It should be noted that each experiment is performed 10 times and their average is recorded to reduce the experimental errors

The evaluation metrics
Six indices have been proposed in the literature so far for the performance analysis of the multi-objective optimization algorithms [62] These indices are explained as follows.

Time
The less run time the better performance of an algorithm.

MID
The mean of deviations of the Pareto solutions from the ideal point I sol = min(z 1 , z 2 ) (in which the values of both objective functions are optimal or the origin of the coordinates can be considered as the ideal point for the minimization problems) is measured by MID (Mean Ideal Distance) index [63]. For computing the MID value of maximizing objective function, the values of Pareto solutions are become inverse and the ideal point is considered (0). The value of MID is calculated using Equation (48).
where ‖ pa − I sol ‖ 2 denotes the Euclidean distance between pa ℇ F(A) and the ideal point. The lower MID value the better performance  of an algorithm.

Diversity
The Diversity index, which is computed using Equation (49), measures the diameter of a space cube used by the end values of objective functions for a set of non-dominated solutions. The more Diversity value the better [64].

RAS
The value of the RAS (Rate of Achievement to two objectives Simultaneously) is calculated using Equation (51) [65].
in which, F i = min{f 1i , f 2i }. The smaller RAS value the better.

SNS
The value of SNS (Spread of Non-dominance Solutions) is calculated by Equation (52) [65].
The higher SNS values the more diversity of solutions, in other words, the better solution quality. After ten runs, the average is reported for each test result, presented in Table 5.
In order to obtain an output for each experiment, all the indices are normalized using Equations (53) and (54), represented in Table 6.
In this normalization method, the indices with a negative nature are converted into the positive indices. The indices are prioritized using Goal Programming based on their importance and a weight is considered for each index accordingly [64,65]. Subsequently, the total weights of the indices for each experiment are calculated using Equation (55) according to the importance weights.
According to the Response values, S N is computed and the parameters' levels are determined. The MINITAB outputs are depicted in Figs. 6 and 7.
The levels of the NSGA-II and MOSFS parameters are represented in Table 7.

Computational results
The proposed model was coded in a notebook system with Intel Core™ i5 processor, 4 Gb RAM, and Microsoft Windows 10 Ultimate operating system. First, the AEC method as an exact method was used through the GAMS optimization software to validate the proposed model. Then, the sensitivity analysis was performed to investigate the effects of changing the main parameters of the model on the optimal values of the objective functions. Subsequently, the NSGA-II and MOSFS metaheuristic algorithms were employed through MATLAB version 2018a software for solving the large size problems. Finally, the performance of NSGA-II and MOSFS were compared in terms of six different indexes.

Model validation
An instance problem with 20 nodes, in which spatial coordinates were specified and Euclidean distances were considered, was designed and solved using the AEC method together with two NSGAII and MOSFS metaheuristic algorithms. At the first level of this problem, there are three factories with capacities of 15,000, 15,000 and 30,000 units. At the second level, there exist different five potential locations for setting up depots with the same capacity of 8000 units. Also, the set-up cost of each depot is 150. In addition, there are 12 customers with deterministic and known demands. The transport fleet of the first echelon is limited and heterogeneous, and each factory has three vehicles with different types. The capacities of these vehicles are 7000, 14000, and 20,000 units, respectively, and the usage fixed costs of vehicles are 30, 20, and 10 monetary units, respectively. The transport fleet of the second echelon is homogenous and unlimited. The capacity of these vehicles is the same and equal to 5000 units, also, the usage fixed cost of

Table 7
The NSGA-II and MOSFS parameters. The Pareto frontiers obtained by using the AEC, NSGA-II, and MOSFS are depicted in Fig. 8. The Pareto frontiers obtained by AEC, NSGA-II, and MOSFS are presented in Fig. 9. As shown in Fig. 9, the established depots and the created routes in the first and second echelons are plotted. depots are located in four locations out of the five potential locations, and routes are originated from depots to deliver goods to customers (6 trips are started from 4 depots; two trips are started from each of depots 5 and 7 and one trip is started from each of depots 6 and 4). Taking into account the fixed cost of set-up depots as well as the fixed cost of using vehicles in the first and second echelon, the total cost is equal to 1650 monetary units, which is the highest total cost among all of the obtained Pareto solutions. Regarding the second objective function, the model chooses more vehicles with less environmental pollution so that the number of trips in both echelons is increased and the distance traveled per trip is reduced which results in reducing CO 2 emissions to the lowest level of 194.548 among all of the obtained Pareto solutions. With regard to the third objective function, the employment and economic development rates rise as a result of setting up four depots which lead to rising the third objective function value to the maximum value of 17.861 among all of the obtained Pareto solutions.

Sensitivity analysis
The effect of the changes in the fixed cost of set-up depots (FD k ) on the values of three objective functions is depicted in Fig. 10. In this scenario, the first objective function value (Z 1 ) is considered optimal and the other model parameters are fixed.
As can be seen in Fig. 10, increasing the amount of depots set-up cost leads to increasing the value of the first objective function (Z 1 ) and the value of the second objective function (Z 2 ). However, the value of the third objective function (Z 3 ) is decreased. Because, if we set up a number of depots, the total costs of the entire system will be increased. Also, the regional development and employment rates Fig. 8. Pareto frontiers obtained by using the AEC, NSGA-II, and MOSFS. will be increased. However, the total amount of CO 2 emissions will be decreased.
The effect of changes in the capacity of depots (CD k ) on the values of three objective functions is illustrated in Fig. 11. In this scenario, the first objective function value (Z 1 ) is also optimal and the other model parameters are fixed.
As shown in Fig. 11, increasing the capacity of the depots leads to decreasing the values of the three objective functions (Z 1 , Z 2 ,Z 3 ). Because, increasing the capacity of the depots results in decreasing the number of established depots that means decreasing the depots   set-up costs and transportation costs. On the other hand, increasing the capacity of the depots leads to decreasing the total amount of CO 2 emissions, because the distance traveled by vehicles is reduced. Eventually, increasing the capacity of the depots results in decreasing the value of the third objective function, which is related to the social responsibility, as increasing the capacity of depots means decreasing the number of established depots and reducing the regional development and employment rates.

Comparison of the performance of the NSGA-II and MOSFS algorithms
Subsequently, the group of I2 problems 31 problems (the specifications of the problems are presented in Table 1) introduced by Sterle [61] were used to compare the performance of two metaheuristic algorithms. For each metaheuristic algorithm, every problem was run ten times and the average was reported as the final solution, the values of all the indices are shown in Table 8.
The comparison of the performance of two NSGA-II and MOSFS metaheuristic algorithms is illustrated in Fig. 10 in terms of Time, MID, Diversity, Spacing, RAS, SNS. As shown in Fig. 12 (a), the solution time of MOSFS is less than the solution time of NSGAII. As a result, MOSFS algorithm performs better according to Time. The less MID value, the better performance of the algorithm. Fig. 12(b) indicates that MOSFS outperforms NSGA-II based on the MID index. The more Diversity value the better. Therefore, the MOSFS algorithm performs better in terms of this measure, displayed in Fig. 12(c). As can be seen in Fig. 12(d), in some problems, MOSFS performs better in terms of Spacing, but, NSGA-II outperforms in some other problems. The less RAS value the better. According to Fig. 12(e), the MOSFS algorithm has a better performance based on the RAS index. The more SNS value the better. As can be seen in Fig. 12(f), there is not considerable difference between the two algorithms for small and medium size problems, however the performance of the MOSFS algorithm is better for the larger problems.
It can be concluded from the above figures that there is no considerable difference in the performance of two algorithms for the small size problems. However, as the problem size grows, this difference rises significantly and the MOSFS algorithm performed better. For example, based on the Time index, as the problem size increases, the computational time for both algorithms increases, but the run time of the MOSFS algorithm increases with a slighter slope than the NSGAII algorithm. As can be seen in Table 8, for the small-sized problems, the NSGA-II algorithm outperforms the MOSFS algorithm in terms of Time, Diversity, RAS, and SNS indices. For the medium-sized problems, NSGA-II outperforms based on Spacing, and MOSFS performs better in terms of the other five indices. Finally, for the large-sized problems, the MOSFS algorithm outperforms in terms of all six indices. In general, NSGAII outperforms only according to Spacing and MOSFS outperforms based on Time, RAS, Diversity, MID, and SNS.

Conclusion
During the decades, governments have paid a lot of attention the increase of economic growth and employment rate. On the other hand, economic and population growth, and technological advancement have had significant detrimental impacts on the environment and natural resources. Therefore, sustainable development based on three main aspects including economy, environment, and social responsibility plays a crucial role in the performance of supply chains. Nowadays, many transportation operations of the supply chains are performed by logistic companies. The concept of open routing has been introduced due to the limited number of available vehicles in a logistic company as a percentage of transportation activities are performed by a logistic firm.
For this purpose, in this paper, both open and closed routes were taken into account in the sustainable 2E-LRP. The hybrid closed and open routes in a sustainable 2E-LRP can be mentioned as the main contribution of the proposed model. A multi-objective mixed integer linear programming (MILP) model including minimization of costs and CO 2 emissions and maximization of social responsibility (creating job opportunities and community development) was developed for the problem in hand. It should be noted that optimizing the three objectives of the problem is not possible as these objectives are in conflict with each other. The proposed model seeks to minimize the total costs by reducing the use of vehicles through reducing the number of transportation routes in the two-echelon distribution network (due to the incurred fixed cost of transportation per vehicle) and increasing the allowable load of vehicles. In addition, this model attempts to minimize the CO 2 emissions by increasing the number of vehicles and lowering their loads (the amount of CO 2 emissions is directly related to the amount of loading). Also, the proposed model tries to increase job opportunities and economic growth by establishing more facilities. If the decision-makers concentrate on increasing social responsibility and reducing pollution, the first (cost) objective function increases significantly increases. Therefore, a trade-off solution should be provided for all of the three aspects of sustainability.
First, a small-sized problem was designed to validate the proposed multi-objective model. The model was solved using the exact AEC method together with two metaheuristic NSGA-II and MOSFS algorithms. Due to the NP-Hardness of the problem, the metaheuristic algorithms named NSGA-II and MOSFS were exploited through the MATLAB software to solve the large size problems and obtain the Pareto solutions. The results showed that the MOSFS algorithm obtained better Pareto frontier in terms of six indices. As some suggestions for future research, it is recommended to consider time dependencies and time windows for customers. Also, the reverse routes may be considered in the model. Moreover, customers can be categorized based on their consumption attitudes.

Author contribution statement
Masoud Hajghani: Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper. Mohammad Ali Forghani: Conceived and designed the experiments. Ali Heidari: Performed the experiments; Analyzed and interpreted the data. Mohammad Khalilzadeh: Analyzed and interpreted the data; Wrote the paper.

Funding statement
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Data availability statement
Data will be made available on request.

Declaration of interest's statement
The authors declare no competing interests.