# Two-Echelon Location-Routing Problem with Time Windows and Transportation Resource Sharing

Authors: Yong Wang; Yaoyao Sun; Xiangyang Guan; Yanyong Guo
Journal: Journal of Advanced Transportation
 (2021)

---

## Abstract

In this work, a two-echelon location-routing problem with time windows and transportation resource sharing (2E-LRPTWTRS) is solved by selecting facility locations and optimizing two-echelon vehicle routes. The optimal solutions improve the efficiency of a logistics network based on the geographical distribution and service time windows of logistics facilities and customers. Furthermore, resource utilization is maximized by enabling resource sharing strategies within and among different logistics facilities simultaneously. The 2E-LRPTWTRS is formulated as a biobjective optimization model, and obtaining the smallest number of required delivery vehicles and the minimum total operating cost are the two objective functions. A two-stage hybrid algorithm composed ofk-means clustering and extended multiobjective particle swarm optimization algorithm is proposed for 2E-LRPTWTRS optimization. A self-adaptive mechanism of flight parameters is introduced and adopted during the iterative process to balance the evolution of particles and improve the efficiency of the two-stage hybrid algorithm. Moreover, 20 small-scale instances are used for an algorithm comparison with multiobjective genetic algorithm and nondominated sorting genetic algorithm-II, and the solutions demonstrate the superiority of the proposed algorithm in optimizing logistics networks. The proposed optimization model and hybrid algorithm are tested by employing a real-world case of 2E-LRPTWTRS in Chongqing, China, and the optimization results verify the positive role of the developed model and algorithm in improving logistics efficiency, reducing operating cost, and saving transportation resources in the operations of two-echelon logistics networks.


---

## Body

## 1. Introduction

For decades, effective modern supply chain management has played an important, positive role in the development of logistics industries [1, 2]. The distribution system in supply chain management involves all operations related to delivering goods from the logistics center to the customers [3]. Therefore, as a key component of the distribution system, solving the vehicle routing problems (VRPs) is essential for achieving effective modern supply chain management [4]. In recent years, the literature in VRPs is evolving toward more complex problems [5]. The complexity stems from various aspects, of which flexibility and integration are the most studied ones [2, 5]. Built upon the pure routing problems, integrated VRPs include a series of decisions, among which the location routing problem (LRP) is an important, for example, [5, 6]. LRPs incorporate network design issues into the traditional VRPs to address the VRPs and the facility location problem simultaneously [7, 8].Traditional two-echelon LRPs (2E-LRPs) involve selecting locations for logistics facilities, assigning customers to appropriate logistics facilities, and optimizing two-echelon vehicle routes [7–9]. However, from a long-term view, the operating costs in pure 2E-LRPs can be increased due to the various requirements of service time within a set up working period [10, 11]. In modern logistics industry, delivery time is increasingly important to guarantee the efficient operations of logistics networks [12, 13]. Thus, several extensions have been conducted on the study of traditional 2E-LRP in the literature, and constraints such as time windows, inventory, and customer demands generally added to traditional 2E-LRPs are not optimized [2, 5]. In a 2E-LRP considering time windows (2E-LRPTW), the effects of violating the service time windows are represented by cost penalties [2, 14, 15], that is, if a customer’s time window is violated, an additional penalty cost is added to the total operating cost of the logistics network. Therefore, 2E-LRPTW can effectively reduce the waiting and delay time in the distribution and improve the efficiency of the logistics network.Moreover, the concept of transportation resource sharing is widely used in urban freight transport industry to reduce the number of required vehicles in the logistics network [16, 17]. In modern society, the growing emphasis on sustainable development has led to the design of resource-saving, environment-friendly logistics networks [18, 19]. As a promising strategy for promoting the sustainability of the logistics network, transportation resource sharing allows each delivery vehicle to run multiple distribution routes of nonoverlapping time windows [17]. Under this sharing strategy, the number of vehicles required for sustaining the operations of a logistics network can be greatly reduced [20]. Therefore, in this paper, the 2E-LRPTW with transportation resource sharing (2E-LRPTWTRS) is proposed to design a sustainable, resource-saving distribution network to achieve the lowest total operating cost as well as the minimum number of required vehicles.The 2E-LRP is defined as a prominent NP-hard problem and usually formulated as an optimization model to find the optimal solution by evolutionary algorithms [2, 21, 22]. Generally, optimization models are proposed with the minimization of total cost as the main objective function in most literature [23–25]. Furthermore, researchers have incorporated clustering into the evolutionary algorithms for simplification of computation owing to factors such as the rapid growth of urban population, the expansion of the market, and the development of logistics industries [4, 26–28]. In the present study, location strategies and vehicle routes are optimized though a biobjective mathematical model, treating the minimization of total cost and number of required vehicles as the two objective functions. Moreover, a two-stage hybrid algorithm including k-means clustering and extended multiobjective particle swarm optimization (EMOPSO) is designed to solve the 2E-LRPTWTRS. The clustering stage seeks to simplify the optimization and provide multiple candidate location strategies, and the second stage optimizes routing with the proposed EMOPSO.The remainder of this paper is organized as follows. Section2 introduces the related studies. Section 3 states the problem and elaborates the setup of 2E-LRPTWTRS. Section 4 articulates the proposed two-stage hybrid algorithm for solving 2E-LRPTWTRS. Section 5 presents the computation results of a real-life case to prove the applicability of the proposed model and algorithm. Section 6 summarizes this study and draws the conclusions.


## 2. Literature Review

A standard two-echelon distribution network involves one level-0 facility and certain number of level-1 facilities [6, 29]. Based on the two-echelon distribution network, the 2E-LRP seeks to find facility locations and optimize vehicle routes simultaneously to satisfy the customers’ demands and has been studied by researchers to design an effective distribution network [7, 8]. Kechmane et al. [30] aimed to address a two-echelon location lot-sizing routing problem, which is addressed through a series of strategies such as location decisions, customers’ assignments, and routing optimization. Ben Mohamed et al. [31] proposed a stochastic two-echelon distribution network design problem including location decisions of depots and transportation scheme, which is applied in addressing the demand uncertainty in the two-echelon network design. Yu et al. [32] focused on waste collection area and proposed a 2E-MOLRP to find the inherent similarities of waste collection applications. Abbassi et al. [33] studied the distribution of nonmedical products in healthcare supply chain logistics and proposed a two-echelon location, distribution problem. None of the studies reviewed above considered the economic effects of delivery time on network design and can result in the increase of the operating cost due to the violation of service time windows.The 2E-LRPTW is studied by researchers to consider the customers’ requirement for service time, which is an extension of the traditional 2E-LRP [2, 6, 34]. In the literature concerning the 2E-LRPTW, researchers commonly convert the time window constraints into a set up penalty cost or a given budget constraint to add to the optimization model [2, 14, 15]. Govindan et al. [23] studied the distribution network in the perishable food industry and integrated time constraints into the sustainable design and optimization of a logistics network due to the increasing concern about the influence of supply chain operations on environment. Bala et al. [35] studied the problem resulting from the distribution of perishable products and proposed a cost-effective delivery plan concerning the consideration of customers’ time windows. Ponboon et al. [34] investigated the LRPTW and discussed the computational results as well as the effect of time windows through a branch-and-price algorithm. Therefore, 2E-LRPTW optimization can contribute to the achievement of an efficient distribution network [15, 36, 37].As another strategy for promoting the achievement of an efficient distribution network, transportation resource sharing is generally introduced and encouraged by researchers in their literature [38–41]. Quintero-Araujo et al. [20] studied resource sharing and considered this concept as an effective, promising strategy to promote the efficient operations of supply chains while addressing the proposed integrated routing problem. Ho and Szeto [42] proposed a multivehicle bike-repositioning problem based on bike-sharing systems in the city and regarded the vehicle-sharing strategy as a green transportation mode. Xu et al. [43] regarded the multiresource allocation scheme as an important factor in the logistics network and integrated the resource-sharing constraint into the optimization model. Transportation resource-sharing strategy is generally encouraged by researchers in the literature concerning the optimization of the two-echelon distribution network. In this study, the transportation resource strategy is introduced, and vehicles can be shared within and among different logistics facilities to reduce transportation resource [17, 44].Moreover, 2E-LRP is a prominent NP-hard problem [45–47]. The 2E-LRP optimization is usually solved by establishing a mathematical model [44, 48, 49]. Döyen et al. [50] developed a two-stage random programming model to solve a humanitarian relief logistics problem and expressed the deterministic equivalent of the model as a mixed-integer linear program. Pichka et al. [51] solved an open 2E-LRP by formulating a three-flow-based mixed-integer linear programming model and evaluated the effectiveness of the model through an extensive experiment. Venkateshan et al. [52] proposed a two-echelon joint continuous-discrete location model, which is applied in the problem of locating a limited number of logistics facilities in a continuous Euclidean space, which can be regarded as the intermediate transshipment points among different stakeholders in supply chain operations. Koç et al. [53] introduced a mixed location-routing problem with time windows and presented a mixed integer programming formulation to minimize vehicle fixed cost and facility cost as well as routing cost. Solving these formulated optimization models requires the design of optimization methodologies and algorithms.Multiple exact or heuristic methodologies have been developed by researchers in LRP to achieve optimization [6, 44, 54]. However, exact approaches are widely applied in solving only small-sized versions of the problem due to its computational limits [55]. In complex, large-scale multiobjective optimization problems and evolutionary methodologies are usually adopted to solve the logistics network optimization, and the clustering method is usually incorporated into the algorithm for a simplification [32, 56–58]. Zamar et al. [59] studied the bale collection optimization and designed a nearest neighbor approach together with a constrained k-means clustering to achieve the minimization of the fuel consumption and travel time simultaneously. Liu et al. [60] proposed a two-step clustering-based hybrid algorithm to solve the two-echelon VRP with mixed vehicles and substantially minimized the transport cost and fuel emissions of the logistics network. Gao et al. [61] developed a hybrid algorithm composed of k-means clustering and ant colony algorithm and introduced k-means clustering for handling the location allocation. Rabbani et al. [62] adopted MOPSO to solve an industrial hazardous waste LRP, seeking to minimize the total cost, site risk, and total transportation risk simultaneously. Nguyen et al. [63] presented four hybrid metaheuristic and constructive heuristics to solve the 2E-LRP arising from transportation applications such as city logistics and employed an additional test to demonstrate the applicability of the algorithm in solving such optimization problems.In the above-reviewed literature, the limitations of 2E-LRP optimization are summarized as follows. (1) Time window constraints in logistics network design are inadequately studied in 2E-LRP optimization. (2) Vehicle sharing is seldom adopted as an effective strategy for promoting sustainable development of logistics networks and improving the utilization of transportation resource in 2E-LRP. (3) A valid mathematical model formulated for the 2E-LRP with time windows and transportation resource sharing is yet to be operationalized. (4) Existing optimization algorithms designed for solving 2E-LRPTWTRS have limited effectiveness and applicability.The contributions of this study are as follows. (1) An effective two-echelon logistics network is designed based on the geographical distribution and time windows of logistics facilities and customers in the proposed 2E-LRPTWTRS. (2) Transportation resource sharing is introduced as an effective strategy for promoting sustainable development. (3) A multiobjective optimization model considering service time windows and resource sharing is proposed to extend the formulation of the 2E-LRP. (4) A two-stage evolutionary algorithm is designed to solve the 2E-LRPTWTRS, and the applicability of the proposed algorithm is proven with a case study.


## 3. Problem Statement and Model Formulation

### 3.1. Problem Statement

The 2E-LRPTWTRS is defined to design an efficient two-echelon logistics network through reasonable location selection of distribution centers and vehicle routing optimization. In the two-echelon distribution network, logistics center (LC) transport cargoes to multiple distribution centers (DCs) by semitrailer trucks in the first echelon, and then, DCs deliver corresponding cargoes to their customers by delivery vehicles in the second echelon. In one working period, each logistics facility including LC and DCs has one operation time window, and each customer has one service time window. By respecting time windows as well as other constraints, the 2E-LRPTWTRS is a more practical problem compared with the traditional 2E-LRP. Figure1 shows the comparison of two-echelon logistics network before and after optimization.Figure 1
Illustration of 2E-LRPTWTRS optimization. (a) Initial logistics network and (b) optimized logistics network.In Figure1(a), the initial logistics network presents an inefficient operation state, although five DCs serve the customers in this distribution area. First, several opened DCs (e.g., DC4 and DC6) serve a few customers, resulting in inefficiency and waste in the utilization of logistics facilities. Second, a series of long-distance deliveries exist due to the unreasonable location selection of DCs and assignment of customers. For instance, for the route (e.g., DC7 ⟶ C1 ⟶ C2 ⟶ C3 ⟶ DC7), the customers are clearly more adjacent to DC1 compared with DC7 while served by DC7 at long distances. Third, vehicles and time window violations are numerous because of unreasonable vehicle scheduling. Therefore, the initial logistics network should be redesigned through relocation selection of DCs and vehicle-routing optimization to obtain improved resource utilization and efficiency. Figure 1(b) shows an optimized logistics network, where DC1, DC2, DC4, and DC6 are selected to serve customers. Through a reasonable location strategy of {DC1, DC2, DC4, DC6} and corresponding vehicle routing optimization, the efficiency of the logistics network is substantially improved. A smaller number of opened DCs enable each opened DC to serve more customers, and the reasonable assignment of customers reduces the distance between each DC and its service customers. The utilization of each opened DC is substantially improved while long-haul deliveries are effectively eliminated. In addition, each vehicle can perform multiple delivery routes with an improved utilization by adopting transportation resource-sharing strategy in the vehicle-routing optimization. For example, three routes include DC1 ⟶ C1 ⟶ C2 ⟶ C3 ⟶ DC1, DC4 ⟶ C13 ⟶ C14 ⟶ C15 ⟶ DC4, and DC6 ⟶ C24 ⟶ C25 ⟶ C26 ⟶ DC6 which are served by the same delivery vehicle. Therefore, the optimized logistics network with transportation resource sharing has a higher efficiency than the initial network.Assuming that the per time unit cost for first-echelon transportation is $35, the per time unit for the second-echelon distribution cost and the penalty cost (earliness and delay penalties) are $20 and $10, respectively. Moreover, the fixed cost of each DC in one working period is $100, and the maintenance costs of each truck and delivery vehicle are $80 and $50, respectively. Table1 presents a result comparison before and after 2E-LRPTRS optimization.Table 1
Comparison before and after 2E-LRPTWTRS optimization.
CaseFixed cost ($)Transportation cost ($)Distribution cost ($)Penalty cost ($)Maintenance cost ($)Total cost ($)Number of trucksNumber of vehiclesBefore optimization7005251320410610356529After optimization4002807600230167013In Table1, a considerable improvement of the logistics network can be achieved through 2E-LRPTWTRS optimization. Total cost decreases from $3565 to $1670, which obtains savings of $1895 including fixed cost savings $300, transportation cost savings $245, distribution cost savings $560, penalty cost savings $410, and maintenance cost savings $380. In addition, the number of trucks is reduced from 2 to 1, and the number of delivery vehicles is reduced from 9 to 3. Therefore, 2E-LRPTWTRS optimization can substantially prompt the efficiency improvement of the two-echelon logistics network and the utilization maximization of transportation resource.


### 3.2. Model Formulation

#### 3.2.1. Definitions

The related notations and variables adopted to formulate the 2E-LRPTWTRS model are defined in Table2. Owing to the reality and generality, the developed model is subject to three assumptions.Table 2
Notations and description used in the 2E-LRPTWTRS model.
DescriptionSetZSet of LC and DCs,Z=0,1,2,3,...,z, and 0 denotes the LCISet of DCs andI=1,2,3,...,iJSet of customers andJ=1,2,3,...,jSSet of semitrailer trucks andS=1,2,3,...,sVSet of delivery vehicles andV=1,2,3,...,vRsSet of sequence for transportation routes performed by trucks and Rs=1,2,3,…,m,s∈SRvSet of sequence for distribution routes performed by vehiclev and Rv=1,2,3,…,n,v∈VUsSet of trucks serving DCsUvSet of delivery vehicles serving customersNsmSet of DCs served by trucks in the mth transportation route, s∈S,m∈RsNvnSet of customers served by vehiclev in the nth distribution route, v∈V,n∈RvParameterCsCapacity of trucksCvCapacity of vehiclevCiCapacity of DCiQiDemand from DCi, i∈IQjDemand from customerj, j∈JDziDistance from LC or DCi to DC d, z,i∈ZDijDistance from DC or customerd to customer c, i,j∈I∪JfsFuel consumption of trucksfvFuel consumption of vehiclevPsDiesel pricePvGasoline priceWNumber of working periods in one yearλePenalty cost per time unit for earlinessλlPenalty cost per time unit for delayαz,βzTime window for logistics facilityz, z∈Zαj,βjTime window for customerj, j∈JDT0smTime trucks departs from LC in the mth route, s∈S,m∈RsDTivnTime vehiclev departs from DC i in the nth route, i∈I,v∈V,n∈RvATzsmTime trucks arrives at LC or DC z in the mth route, z∈Z,s∈S,m∈RsATjvnTime vehiclev arrives at DC or customer j in the nth route, j∈I∪J,v∈V,n∈RvMMA very large numberMsMaintenance cost of trucksMvMaintenance cost of delivery vehiclevFiFixed cost of DC I,i∈IRsNumber of transportation routes served by trucks, s∈SRvNumber of distribution routes served by vehiclev, v∈VNsmNumber of DCs served by trucks in the mth route, s∈S,m∈RsJTotal number of customers to be served in the given distribution areaNvnNumber of customers served by vehiclev in the nth route, v∈V,n∈RvUsNumber of trucks serving DCsUvNumber of vehicles serving customersDecision variableyzismIf trucks travels from logistics facility z to i in the mth route, yzism=1, otherwise, yzism=0, z,i∈Z,s∈S,m∈RsxijvnIf vehiclev travels from DC or customer d to customer c in the nth route, xijvn=1, otherwise, xijvn=0, i,j∈I∪J,v∈V,n∈Rvy0isIf trucks travels from LC 0 to DC i,y0is=1, otherwise, y0is=0, i∈I,s∈SxijvIf vehiclev travels from DC i to customer j, xijv=1, otherwise, xijv=0,i∈I,j∈J,v∈VssIf trucks is used to serve customers, ss=1, otherwise, ss=0,s∈SvvIf delivery vehiclev is used to serve customers, vv=1, otherwise, vv=0,v∈VOiIf DCi is selected to serve customers, Oi=1, otherwise, Oi=0,i∈IOijIf customerj is served by DC i,Oij=1, otherwise, Oij=0,i∈I,j∈JAssumption 1.
Within one working period, each customer is served exactly once, and LC cannot serve customers directly.Assumption 2.
Each candidate DC’s location is given and known, and not all existing DCs are selected.Assumption 3.
Every vehicle can only depart from one DC and must return to the same DC after finishing delivery.


#### 3.2.2. Model Formulation

The 2E-LRPTWTRS is formulated as a biobjective optimization model in this section. The two objectives are minimizing the total cost and the number of delivery vehicles. The objective functions are presented in equations (1) and (2). In Equation (1), the minimization of total cost F1 contains four components, namely, C1, C2, C3, and C4. Equation (2) expresses the minimization of required delivery vehicles:(1)minF1=C1+C2+C3+C4,(2)F2=min∑v∈Vvv⋅min∑i∈I∑j∈J∑n∈Rvxijvn,1,In equation (3), C1 denotes the fixed cost within a working period. Fi⋅Oi is the fixed cost of the opened DCs:(3)C1=∑i∈IFi×Oi.In equation (4), C2 denotes the total transportation cost and maintenance cost for semitrailer trucks within one working period. yzism⋅Lzi⋅fs⋅Ps and Ms/W×Us are the transportation cost and maintenance cost in the first echelon, respectively:(4)C2=∑z,i∈Z∑s∈S∑m∈Rsyzism⋅Lzi⋅fs⋅Ps+MsW×Us.In equation (5), C3 denotes the total distribution cost and maintenance cost for delivery vehicles within one working period. yijvn⋅Lij⋅fv⋅Pv and Mv/W×Uv are the distribution cost and maintenance cost in the second echelon, respectively:(5)C3=∑i,j∈I∪J∑v∈V∑n∈Rvxijvn⋅Lij⋅fv⋅Pv+MvW×Uv.In equation (6), C4 denotes the total penalty cost of trucks and vehicles for earliness or delay. λe⋅maxαi−∑s∈S∑m∈RsATismw,0 and λl⋅max∑s∈S∑m∈RsATism−βi,0 are the penalty cost of semitrailer trucks. λe⋅maxαj−∑v∈V∑n∈RvATjvn,0andλl⋅max∑v∈V∑n∈RvATjvn−βj,0 are the penalty cost of delivery vehicles:(6)C4=∑i∈Iλe⋅maxαi−∑s∈S∑m∈RsATism,0+∑i∈D∑w∈Wλl⋅max∑s∈S∑m∈RsATism−βi,0+∑j∈Jλe⋅maxαj−∑v∈V∑n∈RvATjvn,0+∑c∈Cλl⋅max∑v∈V∑n∈RvATjvn−βj,0,which subject to(7)∑z∈Z∑s∈S∑m∈Rsyzism=1,i∈I,(8)∑i∈I∪J∑v∈V∑n∈Rvxijvn=1,j∈J,(9)∑z∈Zyzism−∑h∈Zyihsm=0,i∈I,s∈S,m∈Rs,(10)∑i∈I∪Jxijvn−∑l∈I∪Jxjlvn=0,j∈J,v∈V,n∈Rv,(11)∑s∈Sss⋅min∑i∈I∑m∈Rsy0ism,1=Us,(12)∑v∈Vvv⋅min∑i∈I∑j∈J∑n∈Rvxijvn,1=Uv,(13)∑i∈Iy0is=Rs,s∈S,(14)∑i∈I∑j∈Jxijv=Rv,4v∈V,(15)∑i∈Iy0ism≥∑i∈Iy0ism+1,s∈S,m∈Rs,m≠Ks,(16)∑i∈I∑j∈Jxijvn≥∑i∈I∑j∈Jxijvn+1,v∈V,n∈Rv,n≠Rv,(17)∑i∈IQi×∑z∈Zyzism≤Cs,s∈S,m∈Rs,(18)∑j∈JQj×∑i∈I∪Jxijvn≤Cv,v∈V,n∈Rv,(19)∑j∈JQj×Oij≤Ci,i∈I,(20)∑z,i∈Nsmyzism≤Nsm−1,s∈S,m∈Rs,(21)∑i,j∈Nvnxijvn≤Nvn−1,v∈V,n∈Rv,(22)α0≤DT0sm≤β0,s∈S,m∈Rs,(23)α0≤AT0sm≤β0,s∈S,m∈Rs,(24)αi⋅∑j∈Jxijvn≤DTivn≤βi⋅∑j∈Jxijvn,i∈I,v∈V,n∈Rv,(25)αi⋅∑j∈Jxjivn≤ATivn≤βi⋅∑j∈Jxjivn,i∈I,v∈V,n∈Rv,(26)AT0sm≤DT0sm+1,s∈S,m∈Rs,m≠Rs,(27)∑i∈IATivn≤∑i∈IDTivn+1,v∈V,n∈Rv,n≠Rv,(28)xijvn≤Oij,i∈I,j∈J,v∈V,n∈Rv,(29)xjhvn≤Oij×Oih,j,h∈J,i∈I,v∈V,n∈Rv,(30)yzism=0,1,z,i∈Z,s∈S,m∈Rs,(31)xijvn=0,1,i,j∈I∪J,v∈V,n∈Rv,(32)y0is=0,1,i∈I,s∈S,(33)xijv=0,1,i∈I,j∈J,v∈V,(34)ss=0,1,s∈S,(35)vv=0,1,v∈V,(36)Oi=0,1,i∈I,(37)Oij=0,1,i∈I,j∈J.Constraints (7) and (8) ensure that each node (i.e., DC or customer) is served exactly once within one service period. Constraints (9) and (10) are flow conservation constraints. Constraints (11) and (12) separately formulate the total number of trucks and delivery vehicles within one service period. Constraints (13) and (14) separately formulate the total number of routes executed by each truck and delivery vehicle within one service period. Constraints (15) and (16) separately ensure the sequence of routes performed by each truck and delivery vehicle. Constraints (17)–(19) are capacity constraints. Constraints (20) and (21) eliminate the subtours. Constraints (22)–(25) ensure that the time window constraints are respected. Constraints (26) and (27) ensure the continuous departure time of two sequential routes performed by each truck and delivery vehicle, respectively. Constraints (28) and (29) ensure that only nodes within the same cluster can be connected. Constraints (30)–(37) define the binary decision variables.


## 3.1. Problem Statement

The 2E-LRPTWTRS is defined to design an efficient two-echelon logistics network through reasonable location selection of distribution centers and vehicle routing optimization. In the two-echelon distribution network, logistics center (LC) transport cargoes to multiple distribution centers (DCs) by semitrailer trucks in the first echelon, and then, DCs deliver corresponding cargoes to their customers by delivery vehicles in the second echelon. In one working period, each logistics facility including LC and DCs has one operation time window, and each customer has one service time window. By respecting time windows as well as other constraints, the 2E-LRPTWTRS is a more practical problem compared with the traditional 2E-LRP. Figure1 shows the comparison of two-echelon logistics network before and after optimization.Figure 1
Illustration of 2E-LRPTWTRS optimization. (a) Initial logistics network and (b) optimized logistics network.In Figure1(a), the initial logistics network presents an inefficient operation state, although five DCs serve the customers in this distribution area. First, several opened DCs (e.g., DC4 and DC6) serve a few customers, resulting in inefficiency and waste in the utilization of logistics facilities. Second, a series of long-distance deliveries exist due to the unreasonable location selection of DCs and assignment of customers. For instance, for the route (e.g., DC7 ⟶ C1 ⟶ C2 ⟶ C3 ⟶ DC7), the customers are clearly more adjacent to DC1 compared with DC7 while served by DC7 at long distances. Third, vehicles and time window violations are numerous because of unreasonable vehicle scheduling. Therefore, the initial logistics network should be redesigned through relocation selection of DCs and vehicle-routing optimization to obtain improved resource utilization and efficiency. Figure 1(b) shows an optimized logistics network, where DC1, DC2, DC4, and DC6 are selected to serve customers. Through a reasonable location strategy of {DC1, DC2, DC4, DC6} and corresponding vehicle routing optimization, the efficiency of the logistics network is substantially improved. A smaller number of opened DCs enable each opened DC to serve more customers, and the reasonable assignment of customers reduces the distance between each DC and its service customers. The utilization of each opened DC is substantially improved while long-haul deliveries are effectively eliminated. In addition, each vehicle can perform multiple delivery routes with an improved utilization by adopting transportation resource-sharing strategy in the vehicle-routing optimization. For example, three routes include DC1 ⟶ C1 ⟶ C2 ⟶ C3 ⟶ DC1, DC4 ⟶ C13 ⟶ C14 ⟶ C15 ⟶ DC4, and DC6 ⟶ C24 ⟶ C25 ⟶ C26 ⟶ DC6 which are served by the same delivery vehicle. Therefore, the optimized logistics network with transportation resource sharing has a higher efficiency than the initial network.Assuming that the per time unit cost for first-echelon transportation is $35, the per time unit for the second-echelon distribution cost and the penalty cost (earliness and delay penalties) are $20 and $10, respectively. Moreover, the fixed cost of each DC in one working period is $100, and the maintenance costs of each truck and delivery vehicle are $80 and $50, respectively. Table1 presents a result comparison before and after 2E-LRPTRS optimization.Table 1
Comparison before and after 2E-LRPTWTRS optimization.
CaseFixed cost ($)Transportation cost ($)Distribution cost ($)Penalty cost ($)Maintenance cost ($)Total cost ($)Number of trucksNumber of vehiclesBefore optimization7005251320410610356529After optimization4002807600230167013In Table1, a considerable improvement of the logistics network can be achieved through 2E-LRPTWTRS optimization. Total cost decreases from $3565 to $1670, which obtains savings of $1895 including fixed cost savings $300, transportation cost savings $245, distribution cost savings $560, penalty cost savings $410, and maintenance cost savings $380. In addition, the number of trucks is reduced from 2 to 1, and the number of delivery vehicles is reduced from 9 to 3. Therefore, 2E-LRPTWTRS optimization can substantially prompt the efficiency improvement of the two-echelon logistics network and the utilization maximization of transportation resource.


## 3.2. Model Formulation

### 3.2.1. Definitions

The related notations and variables adopted to formulate the 2E-LRPTWTRS model are defined in Table2. Owing to the reality and generality, the developed model is subject to three assumptions.Table 2
Notations and description used in the 2E-LRPTWTRS model.
DescriptionSetZSet of LC and DCs,Z=0,1,2,3,...,z, and 0 denotes the LCISet of DCs andI=1,2,3,...,iJSet of customers andJ=1,2,3,...,jSSet of semitrailer trucks andS=1,2,3,...,sVSet of delivery vehicles andV=1,2,3,...,vRsSet of sequence for transportation routes performed by trucks and Rs=1,2,3,…,m,s∈SRvSet of sequence for distribution routes performed by vehiclev and Rv=1,2,3,…,n,v∈VUsSet of trucks serving DCsUvSet of delivery vehicles serving customersNsmSet of DCs served by trucks in the mth transportation route, s∈S,m∈RsNvnSet of customers served by vehiclev in the nth distribution route, v∈V,n∈RvParameterCsCapacity of trucksCvCapacity of vehiclevCiCapacity of DCiQiDemand from DCi, i∈IQjDemand from customerj, j∈JDziDistance from LC or DCi to DC d, z,i∈ZDijDistance from DC or customerd to customer c, i,j∈I∪JfsFuel consumption of trucksfvFuel consumption of vehiclevPsDiesel pricePvGasoline priceWNumber of working periods in one yearλePenalty cost per time unit for earlinessλlPenalty cost per time unit for delayαz,βzTime window for logistics facilityz, z∈Zαj,βjTime window for customerj, j∈JDT0smTime trucks departs from LC in the mth route, s∈S,m∈RsDTivnTime vehiclev departs from DC i in the nth route, i∈I,v∈V,n∈RvATzsmTime trucks arrives at LC or DC z in the mth route, z∈Z,s∈S,m∈RsATjvnTime vehiclev arrives at DC or customer j in the nth route, j∈I∪J,v∈V,n∈RvMMA very large numberMsMaintenance cost of trucksMvMaintenance cost of delivery vehiclevFiFixed cost of DC I,i∈IRsNumber of transportation routes served by trucks, s∈SRvNumber of distribution routes served by vehiclev, v∈VNsmNumber of DCs served by trucks in the mth route, s∈S,m∈RsJTotal number of customers to be served in the given distribution areaNvnNumber of customers served by vehiclev in the nth route, v∈V,n∈RvUsNumber of trucks serving DCsUvNumber of vehicles serving customersDecision variableyzismIf trucks travels from logistics facility z to i in the mth route, yzism=1, otherwise, yzism=0, z,i∈Z,s∈S,m∈RsxijvnIf vehiclev travels from DC or customer d to customer c in the nth route, xijvn=1, otherwise, xijvn=0, i,j∈I∪J,v∈V,n∈Rvy0isIf trucks travels from LC 0 to DC i,y0is=1, otherwise, y0is=0, i∈I,s∈SxijvIf vehiclev travels from DC i to customer j, xijv=1, otherwise, xijv=0,i∈I,j∈J,v∈VssIf trucks is used to serve customers, ss=1, otherwise, ss=0,s∈SvvIf delivery vehiclev is used to serve customers, vv=1, otherwise, vv=0,v∈VOiIf DCi is selected to serve customers, Oi=1, otherwise, Oi=0,i∈IOijIf customerj is served by DC i,Oij=1, otherwise, Oij=0,i∈I,j∈JAssumption 1.
Within one working period, each customer is served exactly once, and LC cannot serve customers directly.Assumption 2.
Each candidate DC’s location is given and known, and not all existing DCs are selected.Assumption 3.
Every vehicle can only depart from one DC and must return to the same DC after finishing delivery.


### 3.2.2. Model Formulation

The 2E-LRPTWTRS is formulated as a biobjective optimization model in this section. The two objectives are minimizing the total cost and the number of delivery vehicles. The objective functions are presented in equations (1) and (2). In Equation (1), the minimization of total cost F1 contains four components, namely, C1, C2, C3, and C4. Equation (2) expresses the minimization of required delivery vehicles:(1)minF1=C1+C2+C3+C4,(2)F2=min∑v∈Vvv⋅min∑i∈I∑j∈J∑n∈Rvxijvn,1,In equation (3), C1 denotes the fixed cost within a working period. Fi⋅Oi is the fixed cost of the opened DCs:(3)C1=∑i∈IFi×Oi.In equation (4), C2 denotes the total transportation cost and maintenance cost for semitrailer trucks within one working period. yzism⋅Lzi⋅fs⋅Ps and Ms/W×Us are the transportation cost and maintenance cost in the first echelon, respectively:(4)C2=∑z,i∈Z∑s∈S∑m∈Rsyzism⋅Lzi⋅fs⋅Ps+MsW×Us.In equation (5), C3 denotes the total distribution cost and maintenance cost for delivery vehicles within one working period. yijvn⋅Lij⋅fv⋅Pv and Mv/W×Uv are the distribution cost and maintenance cost in the second echelon, respectively:(5)C3=∑i,j∈I∪J∑v∈V∑n∈Rvxijvn⋅Lij⋅fv⋅Pv+MvW×Uv.In equation (6), C4 denotes the total penalty cost of trucks and vehicles for earliness or delay. λe⋅maxαi−∑s∈S∑m∈RsATismw,0 and λl⋅max∑s∈S∑m∈RsATism−βi,0 are the penalty cost of semitrailer trucks. λe⋅maxαj−∑v∈V∑n∈RvATjvn,0andλl⋅max∑v∈V∑n∈RvATjvn−βj,0 are the penalty cost of delivery vehicles:(6)C4=∑i∈Iλe⋅maxαi−∑s∈S∑m∈RsATism,0+∑i∈D∑w∈Wλl⋅max∑s∈S∑m∈RsATism−βi,0+∑j∈Jλe⋅maxαj−∑v∈V∑n∈RvATjvn,0+∑c∈Cλl⋅max∑v∈V∑n∈RvATjvn−βj,0,which subject to(7)∑z∈Z∑s∈S∑m∈Rsyzism=1,i∈I,(8)∑i∈I∪J∑v∈V∑n∈Rvxijvn=1,j∈J,(9)∑z∈Zyzism−∑h∈Zyihsm=0,i∈I,s∈S,m∈Rs,(10)∑i∈I∪Jxijvn−∑l∈I∪Jxjlvn=0,j∈J,v∈V,n∈Rv,(11)∑s∈Sss⋅min∑i∈I∑m∈Rsy0ism,1=Us,(12)∑v∈Vvv⋅min∑i∈I∑j∈J∑n∈Rvxijvn,1=Uv,(13)∑i∈Iy0is=Rs,s∈S,(14)∑i∈I∑j∈Jxijv=Rv,4v∈V,(15)∑i∈Iy0ism≥∑i∈Iy0ism+1,s∈S,m∈Rs,m≠Ks,(16)∑i∈I∑j∈Jxijvn≥∑i∈I∑j∈Jxijvn+1,v∈V,n∈Rv,n≠Rv,(17)∑i∈IQi×∑z∈Zyzism≤Cs,s∈S,m∈Rs,(18)∑j∈JQj×∑i∈I∪Jxijvn≤Cv,v∈V,n∈Rv,(19)∑j∈JQj×Oij≤Ci,i∈I,(20)∑z,i∈Nsmyzism≤Nsm−1,s∈S,m∈Rs,(21)∑i,j∈Nvnxijvn≤Nvn−1,v∈V,n∈Rv,(22)α0≤DT0sm≤β0,s∈S,m∈Rs,(23)α0≤AT0sm≤β0,s∈S,m∈Rs,(24)αi⋅∑j∈Jxijvn≤DTivn≤βi⋅∑j∈Jxijvn,i∈I,v∈V,n∈Rv,(25)αi⋅∑j∈Jxjivn≤ATivn≤βi⋅∑j∈Jxjivn,i∈I,v∈V,n∈Rv,(26)AT0sm≤DT0sm+1,s∈S,m∈Rs,m≠Rs,(27)∑i∈IATivn≤∑i∈IDTivn+1,v∈V,n∈Rv,n≠Rv,(28)xijvn≤Oij,i∈I,j∈J,v∈V,n∈Rv,(29)xjhvn≤Oij×Oih,j,h∈J,i∈I,v∈V,n∈Rv,(30)yzism=0,1,z,i∈Z,s∈S,m∈Rs,(31)xijvn=0,1,i,j∈I∪J,v∈V,n∈Rv,(32)y0is=0,1,i∈I,s∈S,(33)xijv=0,1,i∈I,j∈J,v∈V,(34)ss=0,1,s∈S,(35)vv=0,1,v∈V,(36)Oi=0,1,i∈I,(37)Oij=0,1,i∈I,j∈J.Constraints (7) and (8) ensure that each node (i.e., DC or customer) is served exactly once within one service period. Constraints (9) and (10) are flow conservation constraints. Constraints (11) and (12) separately formulate the total number of trucks and delivery vehicles within one service period. Constraints (13) and (14) separately formulate the total number of routes executed by each truck and delivery vehicle within one service period. Constraints (15) and (16) separately ensure the sequence of routes performed by each truck and delivery vehicle. Constraints (17)–(19) are capacity constraints. Constraints (20) and (21) eliminate the subtours. Constraints (22)–(25) ensure that the time window constraints are respected. Constraints (26) and (27) ensure the continuous departure time of two sequential routes performed by each truck and delivery vehicle, respectively. Constraints (28) and (29) ensure that only nodes within the same cluster can be connected. Constraints (30)–(37) define the binary decision variables.


## 3.2.1. Definitions

The related notations and variables adopted to formulate the 2E-LRPTWTRS model are defined in Table2. Owing to the reality and generality, the developed model is subject to three assumptions.Table 2
Notations and description used in the 2E-LRPTWTRS model.
DescriptionSetZSet of LC and DCs,Z=0,1,2,3,...,z, and 0 denotes the LCISet of DCs andI=1,2,3,...,iJSet of customers andJ=1,2,3,...,jSSet of semitrailer trucks andS=1,2,3,...,sVSet of delivery vehicles andV=1,2,3,...,vRsSet of sequence for transportation routes performed by trucks and Rs=1,2,3,…,m,s∈SRvSet of sequence for distribution routes performed by vehiclev and Rv=1,2,3,…,n,v∈VUsSet of trucks serving DCsUvSet of delivery vehicles serving customersNsmSet of DCs served by trucks in the mth transportation route, s∈S,m∈RsNvnSet of customers served by vehiclev in the nth distribution route, v∈V,n∈RvParameterCsCapacity of trucksCvCapacity of vehiclevCiCapacity of DCiQiDemand from DCi, i∈IQjDemand from customerj, j∈JDziDistance from LC or DCi to DC d, z,i∈ZDijDistance from DC or customerd to customer c, i,j∈I∪JfsFuel consumption of trucksfvFuel consumption of vehiclevPsDiesel pricePvGasoline priceWNumber of working periods in one yearλePenalty cost per time unit for earlinessλlPenalty cost per time unit for delayαz,βzTime window for logistics facilityz, z∈Zαj,βjTime window for customerj, j∈JDT0smTime trucks departs from LC in the mth route, s∈S,m∈RsDTivnTime vehiclev departs from DC i in the nth route, i∈I,v∈V,n∈RvATzsmTime trucks arrives at LC or DC z in the mth route, z∈Z,s∈S,m∈RsATjvnTime vehiclev arrives at DC or customer j in the nth route, j∈I∪J,v∈V,n∈RvMMA very large numberMsMaintenance cost of trucksMvMaintenance cost of delivery vehiclevFiFixed cost of DC I,i∈IRsNumber of transportation routes served by trucks, s∈SRvNumber of distribution routes served by vehiclev, v∈VNsmNumber of DCs served by trucks in the mth route, s∈S,m∈RsJTotal number of customers to be served in the given distribution areaNvnNumber of customers served by vehiclev in the nth route, v∈V,n∈RvUsNumber of trucks serving DCsUvNumber of vehicles serving customersDecision variableyzismIf trucks travels from logistics facility z to i in the mth route, yzism=1, otherwise, yzism=0, z,i∈Z,s∈S,m∈RsxijvnIf vehiclev travels from DC or customer d to customer c in the nth route, xijvn=1, otherwise, xijvn=0, i,j∈I∪J,v∈V,n∈Rvy0isIf trucks travels from LC 0 to DC i,y0is=1, otherwise, y0is=0, i∈I,s∈SxijvIf vehiclev travels from DC i to customer j, xijv=1, otherwise, xijv=0,i∈I,j∈J,v∈VssIf trucks is used to serve customers, ss=1, otherwise, ss=0,s∈SvvIf delivery vehiclev is used to serve customers, vv=1, otherwise, vv=0,v∈VOiIf DCi is selected to serve customers, Oi=1, otherwise, Oi=0,i∈IOijIf customerj is served by DC i,Oij=1, otherwise, Oij=0,i∈I,j∈JAssumption 1.
Within one working period, each customer is served exactly once, and LC cannot serve customers directly.Assumption 2.
Each candidate DC’s location is given and known, and not all existing DCs are selected.Assumption 3.
Every vehicle can only depart from one DC and must return to the same DC after finishing delivery.


## 3.2.2. Model Formulation

The 2E-LRPTWTRS is formulated as a biobjective optimization model in this section. The two objectives are minimizing the total cost and the number of delivery vehicles. The objective functions are presented in equations (1) and (2). In Equation (1), the minimization of total cost F1 contains four components, namely, C1, C2, C3, and C4. Equation (2) expresses the minimization of required delivery vehicles:(1)minF1=C1+C2+C3+C4,(2)F2=min∑v∈Vvv⋅min∑i∈I∑j∈J∑n∈Rvxijvn,1,In equation (3), C1 denotes the fixed cost within a working period. Fi⋅Oi is the fixed cost of the opened DCs:(3)C1=∑i∈IFi×Oi.In equation (4), C2 denotes the total transportation cost and maintenance cost for semitrailer trucks within one working period. yzism⋅Lzi⋅fs⋅Ps and Ms/W×Us are the transportation cost and maintenance cost in the first echelon, respectively:(4)C2=∑z,i∈Z∑s∈S∑m∈Rsyzism⋅Lzi⋅fs⋅Ps+MsW×Us.In equation (5), C3 denotes the total distribution cost and maintenance cost for delivery vehicles within one working period. yijvn⋅Lij⋅fv⋅Pv and Mv/W×Uv are the distribution cost and maintenance cost in the second echelon, respectively:(5)C3=∑i,j∈I∪J∑v∈V∑n∈Rvxijvn⋅Lij⋅fv⋅Pv+MvW×Uv.In equation (6), C4 denotes the total penalty cost of trucks and vehicles for earliness or delay. λe⋅maxαi−∑s∈S∑m∈RsATismw,0 and λl⋅max∑s∈S∑m∈RsATism−βi,0 are the penalty cost of semitrailer trucks. λe⋅maxαj−∑v∈V∑n∈RvATjvn,0andλl⋅max∑v∈V∑n∈RvATjvn−βj,0 are the penalty cost of delivery vehicles:(6)C4=∑i∈Iλe⋅maxαi−∑s∈S∑m∈RsATism,0+∑i∈D∑w∈Wλl⋅max∑s∈S∑m∈RsATism−βi,0+∑j∈Jλe⋅maxαj−∑v∈V∑n∈RvATjvn,0+∑c∈Cλl⋅max∑v∈V∑n∈RvATjvn−βj,0,which subject to(7)∑z∈Z∑s∈S∑m∈Rsyzism=1,i∈I,(8)∑i∈I∪J∑v∈V∑n∈Rvxijvn=1,j∈J,(9)∑z∈Zyzism−∑h∈Zyihsm=0,i∈I,s∈S,m∈Rs,(10)∑i∈I∪Jxijvn−∑l∈I∪Jxjlvn=0,j∈J,v∈V,n∈Rv,(11)∑s∈Sss⋅min∑i∈I∑m∈Rsy0ism,1=Us,(12)∑v∈Vvv⋅min∑i∈I∑j∈J∑n∈Rvxijvn,1=Uv,(13)∑i∈Iy0is=Rs,s∈S,(14)∑i∈I∑j∈Jxijv=Rv,4v∈V,(15)∑i∈Iy0ism≥∑i∈Iy0ism+1,s∈S,m∈Rs,m≠Ks,(16)∑i∈I∑j∈Jxijvn≥∑i∈I∑j∈Jxijvn+1,v∈V,n∈Rv,n≠Rv,(17)∑i∈IQi×∑z∈Zyzism≤Cs,s∈S,m∈Rs,(18)∑j∈JQj×∑i∈I∪Jxijvn≤Cv,v∈V,n∈Rv,(19)∑j∈JQj×Oij≤Ci,i∈I,(20)∑z,i∈Nsmyzism≤Nsm−1,s∈S,m∈Rs,(21)∑i,j∈Nvnxijvn≤Nvn−1,v∈V,n∈Rv,(22)α0≤DT0sm≤β0,s∈S,m∈Rs,(23)α0≤AT0sm≤β0,s∈S,m∈Rs,(24)αi⋅∑j∈Jxijvn≤DTivn≤βi⋅∑j∈Jxijvn,i∈I,v∈V,n∈Rv,(25)αi⋅∑j∈Jxjivn≤ATivn≤βi⋅∑j∈Jxjivn,i∈I,v∈V,n∈Rv,(26)AT0sm≤DT0sm+1,s∈S,m∈Rs,m≠Rs,(27)∑i∈IATivn≤∑i∈IDTivn+1,v∈V,n∈Rv,n≠Rv,(28)xijvn≤Oij,i∈I,j∈J,v∈V,n∈Rv,(29)xjhvn≤Oij×Oih,j,h∈J,i∈I,v∈V,n∈Rv,(30)yzism=0,1,z,i∈Z,s∈S,m∈Rs,(31)xijvn=0,1,i,j∈I∪J,v∈V,n∈Rv,(32)y0is=0,1,i∈I,s∈S,(33)xijv=0,1,i∈I,j∈J,v∈V,(34)ss=0,1,s∈S,(35)vv=0,1,v∈V,(36)Oi=0,1,i∈I,(37)Oij=0,1,i∈I,j∈J.Constraints (7) and (8) ensure that each node (i.e., DC or customer) is served exactly once within one service period. Constraints (9) and (10) are flow conservation constraints. Constraints (11) and (12) separately formulate the total number of trucks and delivery vehicles within one service period. Constraints (13) and (14) separately formulate the total number of routes executed by each truck and delivery vehicle within one service period. Constraints (15) and (16) separately ensure the sequence of routes performed by each truck and delivery vehicle. Constraints (17)–(19) are capacity constraints. Constraints (20) and (21) eliminate the subtours. Constraints (22)–(25) ensure that the time window constraints are respected. Constraints (26) and (27) ensure the continuous departure time of two sequential routes performed by each truck and delivery vehicle, respectively. Constraints (28) and (29) ensure that only nodes within the same cluster can be connected. Constraints (30)–(37) define the binary decision variables.


## 4. Solution Methodology

This section introduces how PSO is adapted to address multiobjective optimization problems, that is, EMOPSO. As a two-stage hybrid algorithm, the proposed EMOPSO first introducesk-means clustering to simplify the 2E-LRPTWTRS and provide multiple candidate location strategies, and then, the MOPSO is conducted to solve the two-echelon routing optimization in terms of the different provided location strategies. Figure 2 demonstrates the solving procedure of EMOPSO applied in 2E-LRPTWTRS optimization. The involved parameters adopted in Figure 2 are defined as follows: Iter is the count of optimization runs, max_Iter is the maximum optimization runs, t denotes the number of iterations, max_t denotes the maximum iterations, and gbest and pbest,l represent the global best position and the lth particle’s personal best position, respectively.Figure 2
Flowchart of the hybrid algorithm for 2E-LRPTWTRS optimization.Figure2 shows that k-means clustering and the proposed EMOPSO constitute the two-stage hybrid algorithm. In Stage 1, based on customers’ geographical coordinates, k-means clustering enables customers with geographically close distance to be set in the same cluster, remarkably simplifying the following vehicle routing optimization. In Stage 2, the proposed EMOPSO algorithm provides the optimized vehicle routes in terms of the different location strategies obtained in the first stage. First, the particles in the swarm and the external repository are initialized. Second, each particle is evaluated to determine the dominance relations. Third, the external repository is updated, and the adaptive grids are constructed. Fourth, gbest and pbest,l are determined and then the particles are updated. Fifth, if the count of optimization runs has reached the set, the maximum value is checked: if yes, the optimal results are reported; otherwise, the process returns to the second step, and the following steps are repeated.


### 4.1.K-Means Clustering

K-means clustering is a commonly method adopted for partitioning the data to simplify the optimization problems [2, 27, 56]. Given a set of customer points J = {1, 2, 3, …, |J|}, where each customer point is a two-dimensional vector including latitude and longitude. With the minimization of dissimilarity as an objective function, k -means clustering algorithm is used to partition the customers into k clusters A = {A1, A2,A3, …, Ak}. Based on the Euclidean distance from each customer point j to the corresponding center, the objective function B is defined as(38)B=min∑a=1k∑j=1Jzaj⋅d2j,Aa,where if j ∈ Aa, zaj = 1; otherwise, zaj = 0. d2(j, Aa) is the distance between customer j and the center. Therefore, the objective function B is to minimize the total distance between the customers and the corresponding center. Algorithm 1 shows the procedure of k-means clustering algorithm.Algorithm 1:K-means clustering.
Input: the number of clusters k and set of |J| customer points J = {1, 2, 3, …, |J|}Output: clustering results including k centers, k clusters, and membership matrixSteps:(1)
Initialize thek centers by selecting k customer points from J randomly(2)
Repeat:(i)
Assign each customer point to the closest cluster by calculating and comparing the distance between the customer point and each cluster center(ii)
Calculate the objective functionB using equation (2)(iii)
Update each cluster center(3)
Until each center stays unchangeable(4)
Export the resultsIn Algorithm1, the k-means clustering algorithm is conducted with the following procedures. First, the k centers are initialized by selecting k customer points from J randomly, and the initial membership matrix is constructed. Second, all the customer points are traversed, and the distances between the customer points and each cluster center are computed. Third, the distance between the customer point and each cluster center is compared, and then, the customer to the closest cluster is assigned. Fourth, the above steps are repeated until each center stays unchangeable. The candidate location strategies can be provided through the clustering results, and then, the vehicle routes are initialized and optimized in terms of different location strategies with the EMOPSO algorithm.


### 4.2. EMOPSO Algorithm

In the PSO first developed by Kennedy and Eberhart [64, 65], the potential solutions are represented with a swarm of particles. For the PSO algorithm, two parts, namely, position and velocity, identify each particle, and two leaders, namely, personal best (pbest,l) and global best (gbest), update each particle. Several modifications are made in EMOPSO to adapt PSO in multiobjective optimization problems: (1) the selection of gbest in particle updating, (2) the determination of optimal Pareto front, and (3) the introduction of self-adaptive flight parameter mechanism. The procedure of the proposed EMOPSO is shown in Algorithm 2. The involved operations including the external repository and the self-adaptive flight parameter mechanism are detailed in the following sections.Algorithm 2: EMOPSO.
Input:(1)
max_Iter: the maximum iteration(2)
nP: the particle size(3)
nR: the size of external repository(4)
wl: the lth particle inertia weight(5)
c1l: the lth particle personal learning coefficient(6)
c2l: the lth particle global learning coefficient(7)
nG: the count of grids per dimensionOutput: the optimal solutions(1)
Initialize particles in the swarm (position and velocity)(2)
Evaluate each particle by calculating the corresponding objective function(3)
Initialize the external repository using nondominated solution maintaining(4)
Iter = 0(5)
ForIter = 1 : max_Iter(6)
For each particle(7)
Update each particle’s position and velocity(8)
Evaluate the newly generated particle(9)
Update thepbest,l(10)
End for(11)
Update the external repository and control the repository size(12)
Selectgbest from the external repository by adaptive grids(13)
Adjust the parameters by self-adaptive flight parameter mechanism(14)
Iter = Iter + 1(15)
End for(16)
Report the optimal results


#### 4.2.1. External Repository and Selection ofgbest

The main objective of constructing an external repository is to store the nondominated particles during the updating. The repository has two key parts, namely, repository controller and grid. The repository controller is introduced to estimate whether the particle should be added to the repository. At each iteration, the nondominated particles are compared with the existing particles stored in the external repository to decide whether the existing particle should be replaced. Figure3 illustrates the possible scenarios when conducting repository control. The rectangle is the external repository, the circle denotes the particle, the blue circle represents the newly generated particles, and the red circle denotes the particle dominated by the new particle.Figure 3
Possible scenarios for repository controlling: (a) scenario 1, (b) scenario 2, (c) scenario 3, and (d) scenario 4.
(a)(b)(c)(d)In Figure3(a), the repository is empty, and the current new nondominated particle can be accepted. In Figure 3(b), none of the existing particles within the repository is dominated by the new particle, and the new particle cannot be accepted to add to the repository. In Figure 3(c), particles are dominated by the new particle, and the dominated particles should be replaced by the new particle. Finally, in Figure 3(d), the stored particles reach the allowed size of the repository, and the adaptive grids can be constructed. The adaptive grids are established to obtain a well-distributed Pareto front. The core idea is to divide the objective function space into multiple regions, and Figure 4 presents the idea of the adaptive grids.Figure 4
Illustration of adaptive grid idea. (a) The insertion of a new particle and (b) the insertion of a new particle outside the current grids.
(a)(b)In Figure4, UP represents the set of newly generated particles. In Figure 4(a), the objective function space is divided into 7 × 7 grids, and the new particle replaces its dominated particle and is added to the grids. In Figure 4(b), the new particle is outside the current grids; then, the grids should be reconstructed, and each particle within it should be relocated.


#### 4.2.2. Self-adaptive Flight Parameter Mechanism

In the iteration,gbest and pbest,l, the involved parameters, are key components for balancing the evolution status such as stagnation, convergence, and diversity. Therefore, to adjust and determine these flight parameters, a self-adaptive flight parameter mechanism is introduced to balance the global search and local exploitation by collecting the dominating relation and diversity information. In the evolutionary process of the external repository, after a dominance test, if the dominated particles are replaced and discarded, then parameters wl and c1l should be smaller, whereas parameter c2l should be larger. Moreover, if the dominated particles are retained in the repository, then parameters wl and c1l should be larger, whereas parameter c2l should be smaller. Therefore, the self-adaptive flight parameter mechanism is designed based on dominating relation as(39)APlt=Lmint+LmaxtLmaxt+Llt,where APl(t) is the lth particle’s adaptive parameter, Lmin(t) is the minimum distance from gbest among all the particles, Lmax(t) is the maximum distance from gbest among all the particles, Ll(t) is the distance between gbest and the lth particle, wl(t) is the lth particle’s inertia weight at tth iteration, and c1l(t) and c2l(t) represent the lth particle’s flight parameters at tth iteration. The adaptive parameters can be expressed as equations (40)–(42):(40)wlt=wlt−1,pbest,lt−1=pbest,lt,wlt−1×1−APlt,pbest,lt−1≺pbest,lt,wlt−1×1+APlt,pbest,lt−1≻pbest,lt,(41)c1lt=c1lt−1,pbest,lt−1=pbest,lt,c1lt−1×1−APlt,pbest,lt−1≺pbest,lt,c1lt−1×1+APlt,pbest,lt−1≻pbest,lt,(42)c2lt=c2lt−1,pbest,lt−1=pbest,lt,c2lt−1×1−APlt,pbest,lt−1≺pbest,lt,c2lt−1×1+APlt,pbest,lt−1≻pbest,lt.During the search, such a self-adaptive flight parameter mechanism can help modify the flight parameters, which enables the proposed EMOPSO algorithm to generate a better optimal solution, that is, the convergence during the search of EMOPSO can be effectively reduced and avoided by the self-adaptive flight parameter mechanism, which is beneficial for pushing the Pareto front forward.


## 4.1.K-Means Clustering

K-means clustering is a commonly method adopted for partitioning the data to simplify the optimization problems [2, 27, 56]. Given a set of customer points J = {1, 2, 3, …, |J|}, where each customer point is a two-dimensional vector including latitude and longitude. With the minimization of dissimilarity as an objective function, k -means clustering algorithm is used to partition the customers into k clusters A = {A1, A2,A3, …, Ak}. Based on the Euclidean distance from each customer point j to the corresponding center, the objective function B is defined as(38)B=min∑a=1k∑j=1Jzaj⋅d2j,Aa,where if j ∈ Aa, zaj = 1; otherwise, zaj = 0. d2(j, Aa) is the distance between customer j and the center. Therefore, the objective function B is to minimize the total distance between the customers and the corresponding center. Algorithm 1 shows the procedure of k-means clustering algorithm.Algorithm 1:K-means clustering.
Input: the number of clusters k and set of |J| customer points J = {1, 2, 3, …, |J|}Output: clustering results including k centers, k clusters, and membership matrixSteps:(1)
Initialize thek centers by selecting k customer points from J randomly(2)
Repeat:(i)
Assign each customer point to the closest cluster by calculating and comparing the distance between the customer point and each cluster center(ii)
Calculate the objective functionB using equation (2)(iii)
Update each cluster center(3)
Until each center stays unchangeable(4)
Export the resultsIn Algorithm1, the k-means clustering algorithm is conducted with the following procedures. First, the k centers are initialized by selecting k customer points from J randomly, and the initial membership matrix is constructed. Second, all the customer points are traversed, and the distances between the customer points and each cluster center are computed. Third, the distance between the customer point and each cluster center is compared, and then, the customer to the closest cluster is assigned. Fourth, the above steps are repeated until each center stays unchangeable. The candidate location strategies can be provided through the clustering results, and then, the vehicle routes are initialized and optimized in terms of different location strategies with the EMOPSO algorithm.


## 4.2. EMOPSO Algorithm

In the PSO first developed by Kennedy and Eberhart [64, 65], the potential solutions are represented with a swarm of particles. For the PSO algorithm, two parts, namely, position and velocity, identify each particle, and two leaders, namely, personal best (pbest,l) and global best (gbest), update each particle. Several modifications are made in EMOPSO to adapt PSO in multiobjective optimization problems: (1) the selection of gbest in particle updating, (2) the determination of optimal Pareto front, and (3) the introduction of self-adaptive flight parameter mechanism. The procedure of the proposed EMOPSO is shown in Algorithm 2. The involved operations including the external repository and the self-adaptive flight parameter mechanism are detailed in the following sections.Algorithm 2: EMOPSO.
Input:(1)
max_Iter: the maximum iteration(2)
nP: the particle size(3)
nR: the size of external repository(4)
wl: the lth particle inertia weight(5)
c1l: the lth particle personal learning coefficient(6)
c2l: the lth particle global learning coefficient(7)
nG: the count of grids per dimensionOutput: the optimal solutions(1)
Initialize particles in the swarm (position and velocity)(2)
Evaluate each particle by calculating the corresponding objective function(3)
Initialize the external repository using nondominated solution maintaining(4)
Iter = 0(5)
ForIter = 1 : max_Iter(6)
For each particle(7)
Update each particle’s position and velocity(8)
Evaluate the newly generated particle(9)
Update thepbest,l(10)
End for(11)
Update the external repository and control the repository size(12)
Selectgbest from the external repository by adaptive grids(13)
Adjust the parameters by self-adaptive flight parameter mechanism(14)
Iter = Iter + 1(15)
End for(16)
Report the optimal results


### 4.2.1. External Repository and Selection ofgbest

The main objective of constructing an external repository is to store the nondominated particles during the updating. The repository has two key parts, namely, repository controller and grid. The repository controller is introduced to estimate whether the particle should be added to the repository. At each iteration, the nondominated particles are compared with the existing particles stored in the external repository to decide whether the existing particle should be replaced. Figure3 illustrates the possible scenarios when conducting repository control. The rectangle is the external repository, the circle denotes the particle, the blue circle represents the newly generated particles, and the red circle denotes the particle dominated by the new particle.Figure 3
Possible scenarios for repository controlling: (a) scenario 1, (b) scenario 2, (c) scenario 3, and (d) scenario 4.
(a)(b)(c)(d)In Figure3(a), the repository is empty, and the current new nondominated particle can be accepted. In Figure 3(b), none of the existing particles within the repository is dominated by the new particle, and the new particle cannot be accepted to add to the repository. In Figure 3(c), particles are dominated by the new particle, and the dominated particles should be replaced by the new particle. Finally, in Figure 3(d), the stored particles reach the allowed size of the repository, and the adaptive grids can be constructed. The adaptive grids are established to obtain a well-distributed Pareto front. The core idea is to divide the objective function space into multiple regions, and Figure 4 presents the idea of the adaptive grids.Figure 4
Illustration of adaptive grid idea. (a) The insertion of a new particle and (b) the insertion of a new particle outside the current grids.
(a)(b)In Figure4, UP represents the set of newly generated particles. In Figure 4(a), the objective function space is divided into 7 × 7 grids, and the new particle replaces its dominated particle and is added to the grids. In Figure 4(b), the new particle is outside the current grids; then, the grids should be reconstructed, and each particle within it should be relocated.


### 4.2.2. Self-adaptive Flight Parameter Mechanism

In the iteration,gbest and pbest,l, the involved parameters, are key components for balancing the evolution status such as stagnation, convergence, and diversity. Therefore, to adjust and determine these flight parameters, a self-adaptive flight parameter mechanism is introduced to balance the global search and local exploitation by collecting the dominating relation and diversity information. In the evolutionary process of the external repository, after a dominance test, if the dominated particles are replaced and discarded, then parameters wl and c1l should be smaller, whereas parameter c2l should be larger. Moreover, if the dominated particles are retained in the repository, then parameters wl and c1l should be larger, whereas parameter c2l should be smaller. Therefore, the self-adaptive flight parameter mechanism is designed based on dominating relation as(39)APlt=Lmint+LmaxtLmaxt+Llt,where APl(t) is the lth particle’s adaptive parameter, Lmin(t) is the minimum distance from gbest among all the particles, Lmax(t) is the maximum distance from gbest among all the particles, Ll(t) is the distance between gbest and the lth particle, wl(t) is the lth particle’s inertia weight at tth iteration, and c1l(t) and c2l(t) represent the lth particle’s flight parameters at tth iteration. The adaptive parameters can be expressed as equations (40)–(42):(40)wlt=wlt−1,pbest,lt−1=pbest,lt,wlt−1×1−APlt,pbest,lt−1≺pbest,lt,wlt−1×1+APlt,pbest,lt−1≻pbest,lt,(41)c1lt=c1lt−1,pbest,lt−1=pbest,lt,c1lt−1×1−APlt,pbest,lt−1≺pbest,lt,c1lt−1×1+APlt,pbest,lt−1≻pbest,lt,(42)c2lt=c2lt−1,pbest,lt−1=pbest,lt,c2lt−1×1−APlt,pbest,lt−1≺pbest,lt,c2lt−1×1+APlt,pbest,lt−1≻pbest,lt.During the search, such a self-adaptive flight parameter mechanism can help modify the flight parameters, which enables the proposed EMOPSO algorithm to generate a better optimal solution, that is, the convergence during the search of EMOPSO can be effectively reduced and avoided by the self-adaptive flight parameter mechanism, which is beneficial for pushing the Pareto front forward.


## 4.2.1. External Repository and Selection ofgbest

The main objective of constructing an external repository is to store the nondominated particles during the updating. The repository has two key parts, namely, repository controller and grid. The repository controller is introduced to estimate whether the particle should be added to the repository. At each iteration, the nondominated particles are compared with the existing particles stored in the external repository to decide whether the existing particle should be replaced. Figure3 illustrates the possible scenarios when conducting repository control. The rectangle is the external repository, the circle denotes the particle, the blue circle represents the newly generated particles, and the red circle denotes the particle dominated by the new particle.Figure 3
Possible scenarios for repository controlling: (a) scenario 1, (b) scenario 2, (c) scenario 3, and (d) scenario 4.
(a)(b)(c)(d)In Figure3(a), the repository is empty, and the current new nondominated particle can be accepted. In Figure 3(b), none of the existing particles within the repository is dominated by the new particle, and the new particle cannot be accepted to add to the repository. In Figure 3(c), particles are dominated by the new particle, and the dominated particles should be replaced by the new particle. Finally, in Figure 3(d), the stored particles reach the allowed size of the repository, and the adaptive grids can be constructed. The adaptive grids are established to obtain a well-distributed Pareto front. The core idea is to divide the objective function space into multiple regions, and Figure 4 presents the idea of the adaptive grids.Figure 4
Illustration of adaptive grid idea. (a) The insertion of a new particle and (b) the insertion of a new particle outside the current grids.
(a)(b)In Figure4, UP represents the set of newly generated particles. In Figure 4(a), the objective function space is divided into 7 × 7 grids, and the new particle replaces its dominated particle and is added to the grids. In Figure 4(b), the new particle is outside the current grids; then, the grids should be reconstructed, and each particle within it should be relocated.


## 4.2.2. Self-adaptive Flight Parameter Mechanism

In the iteration,gbest and pbest,l, the involved parameters, are key components for balancing the evolution status such as stagnation, convergence, and diversity. Therefore, to adjust and determine these flight parameters, a self-adaptive flight parameter mechanism is introduced to balance the global search and local exploitation by collecting the dominating relation and diversity information. In the evolutionary process of the external repository, after a dominance test, if the dominated particles are replaced and discarded, then parameters wl and c1l should be smaller, whereas parameter c2l should be larger. Moreover, if the dominated particles are retained in the repository, then parameters wl and c1l should be larger, whereas parameter c2l should be smaller. Therefore, the self-adaptive flight parameter mechanism is designed based on dominating relation as(39)APlt=Lmint+LmaxtLmaxt+Llt,where APl(t) is the lth particle’s adaptive parameter, Lmin(t) is the minimum distance from gbest among all the particles, Lmax(t) is the maximum distance from gbest among all the particles, Ll(t) is the distance between gbest and the lth particle, wl(t) is the lth particle’s inertia weight at tth iteration, and c1l(t) and c2l(t) represent the lth particle’s flight parameters at tth iteration. The adaptive parameters can be expressed as equations (40)–(42):(40)wlt=wlt−1,pbest,lt−1=pbest,lt,wlt−1×1−APlt,pbest,lt−1≺pbest,lt,wlt−1×1+APlt,pbest,lt−1≻pbest,lt,(41)c1lt=c1lt−1,pbest,lt−1=pbest,lt,c1lt−1×1−APlt,pbest,lt−1≺pbest,lt,c1lt−1×1+APlt,pbest,lt−1≻pbest,lt,(42)c2lt=c2lt−1,pbest,lt−1=pbest,lt,c2lt−1×1−APlt,pbest,lt−1≺pbest,lt,c2lt−1×1+APlt,pbest,lt−1≻pbest,lt.During the search, such a self-adaptive flight parameter mechanism can help modify the flight parameters, which enables the proposed EMOPSO algorithm to generate a better optimal solution, that is, the convergence during the search of EMOPSO can be effectively reduced and avoided by the self-adaptive flight parameter mechanism, which is beneficial for pushing the Pareto front forward.


## 5. Computational Experiments

### 5.1. Case Study

A real case of 2E-LRPTWTRS in Chongqing, China, is employed to present the effectiveness of the proposed optimization model and hybrid algorithm. The logistics network composed of one LC, six candidate DCs (DC1, DC2, …, DC6), and 150 customers (C1, C2, …, C150) in this selected distribution area is shown in Figure5. Here, LC is represented as pentagon, stars, squares, crosses, triangles, and circles, and diamonds refer to the DCs and their corresponding service customers.Figure 5
Spatial distribution of logistics facilities and customers.


#### 5.1.1. Data Description

According to actual surveys and related references [2, 17, 37], the involved parameter values used in the proposed optimization model and EMOPSO can be summarized as follows: capacity of the truck Cs = 1000, capacity of the vehicle Cv = 200, fuel consumption of the truck fs = 0.25, fuel consumption of the vehicle fv = 0.2, diesel price Ps = 20, gasoline price Pv = 15, maintenance cost of the truck Ms = 1500, maintenance cost of the vehicle Mv = 500, penalty cost per time unit for earliness λe = 10, penalty cost per time unit for delay λl = 25, maximum iteration max_Iter = 600, particle size nP = 150, the size of external repository nR = 75, the count of grids per dimension nG = 7, initial particle inertia weight w = 1, initial particle personal learning coefficient c1 = 1, and initial particle global learning coefficient c2 = 2. Moreover, Tables 3 and 4 show the initial characteristics of logistics facilities and initial distribution routes, respectively.Table 3
Initial characteristics of logistics facilities.
FacilitiesSymbolFixed cost ($)Number of service customersTime windowsLC——[0, 22]DC120026[1, 19]DC219024[0.5, 17.5]DC322025[0, 18]DC425026[2.5, 21]DC517025[2, 18]DC624024[1.5, 18]Total1270150Table 4
Initial distribution routes of logistics network.
Distribution centersOptimal distribution routesNumber of vehiclesDC1DC1 ⟶ C128 ⟶ C32 ⟶ C113 ⟶ C123 ⟶ C106 ⟶ DC1DC1 ⟶ C108 ⟶ C146 ⟶ C114 ⟶ C143 ⟶ C23 ⟶ DC1DC1 ⟶ C88 ⟶ C145 ⟶ C59 ⟶ C9 ⟶ C7 ⟶ DC1DC1 ⟶ C71 ⟶ C73 ⟶ C82 ⟶ C17 ⟶ DC1DC1 ⟶ C4 ⟶ C69 ⟶ C116 ⟶ C3 ⟶ C86 ⟶ DC131DC2DC2 ⟶ C38 ⟶ C12 ⟶ C34 ⟶ C97 ⟶ DC2DC2 ⟶ C62 ⟶ C142 ⟶ C115 ⟶ C8 ⟶ C14 ⟶ DC2DC2 ⟶ C148 ⟶ C57 ⟶ C91 ⟶ C84 ⟶ DC2DC2 ⟶ C21 ⟶ C45 ⟶ C37 ⟶ C129 ⟶ DC2DC2 ⟶ C42 ⟶ C122 ⟶ C56 ⟶ C29 ⟶ C44 ⟶ DC2DC2 ⟶ C51 ⟶ C121 ⟶ C112 ⟶ C26 ⟶ C139 ⟶ DC2DC3DC3 ⟶ C50 ⟶ C104 ⟶ C79 ⟶ C70 ⟶ DC3DC3 ⟶ C149 ⟶ C80 ⟶ C43 ⟶ C63 ⟶ C118 ⟶ C132 ⟶ DC3DC3 ⟶ C133 ⟶ C101 ⟶ C81 ⟶ C92 ⟶ DC3DC3 ⟶ C95 ⟶ C109 ⟶ C120 ⟶ C2 ⟶ DC3DC3 ⟶ C20 ⟶ C5 ⟶ C72 ⟶ DC3DC3 ⟶ C24 ⟶ C100 ⟶ C22 ⟶ C105 ⟶ DC3DC4DC4 ⟶ C39 ⟶ C64 ⟶ C16 ⟶ C150 ⟶ C78 ⟶ DC4DC4 ⟶ C58 ⟶ C138 ⟶ C135 ⟶ C130 ⟶ DC4DC4 ⟶ C80⟶ C52 ⟶ C127 ⟶ C96 ⟶ C94 ⟶ DC4DC4 ⟶ C34 ⟶ C49 ⟶ C119 ⟶ C33 ⟶ C15 ⟶ DC4DC4 ⟶ C25 ⟶ C65 ⟶ C53 ⟶ C89 ⟶ DC4DC5DC5 ⟶ C48 ⟶ C75 ⟶ C99 ⟶ C55 ⟶ C107 ⟶ DC5DC5 ⟶ C41 ⟶ C141 ⟶ C19 ⟶ C30 ⟶ C28 ⟶ DC5DC5 ⟶ C6 ⟶ C61 ⟶ C67 ⟶ C66 ⟶ DC5DC5 ⟶ C74 ⟶ C134 ⟶ C137 ⟶ C90 ⟶ CDC ⟶ DC5DC5 ⟶ C140 ⟶ C40 ⟶ C11 ⟶ C46 ⟶ C18 ⟶ DC5DC6DC6 ⟶ C10 ⟶ C124 ⟶ C27 ⟶ C13 ⟶ C111 ⟶ DC6DC6 ⟶ C144 ⟶ C147 ⟶ C1 ⟶ C98 ⟶ C68 ⟶ C110 ⟶ DC6DC6 ⟶ C117 ⟶ C47 ⟶ C126 ⟶ C131 ⟶ C54 ⟶ DC6DC6 ⟶ C125 ⟶ C60 ⟶ C103 ⟶ C35 ⟶ C76 ⟶ DC6Tables3 and 4 show that, with all the existing DCs selected for serving customers, each DC serves a few customers and results in a low utilization of logistics facilities. Moreover, the delivery routes without vehicle sharing can lead to a large number of delivery vehicles for the two-echelon logistics network. Therefore, determining the reasonable location strategy and routing scheme through 2E-LRPTWTRS optimization to reduce the total cost is important for optimizing the logistics network.


#### 5.1.2. Optimization Results

The optimization of 2E-LRPTWTRS can achieve a substantial simplification throughk-means clustering. Before the operation of the proposed EMOPSO to optimize the vehicle routes, an effective customer clustering can avoid searching several unnecessary solutions and increasing the computation time, which seriously improves the efficiency of the algorithm. Therefore, to avoid such unreasonable phenomenon, customers are clustered through k-means clustering. In customer clustering, the main challenge is the determination of the k value. Silhouette analysis is introduced to evaluate the performance of customer clustering for each possible k value and address this challenge [2, 4]. The data description shows that 150 customers are in this studied case, and the set of possible k values for the case is reasonably determined as {3, 4, 5}. Figure 6 shows the customer clustering results and silhouette values in terms of possible k values.Figure 6
Customer (a) clustering results and (b) silhouette values.
(a)(b)Figure6(a) shows that three different clustering results are obtained because of the three different possible k values, and the circles with the same color are assigned to the same cluster. Moreover, one center is in each cluster, and the customers of the same cluster are with shorter distance to the corresponding center compared with the other centers. In Figure 6(b), the silhouette values are provided for evaluating the clustering performance, and the majority of the elements have a positive value while few negative values are observed for each clustering scenario, which denotes that a good clustering performance can be obtained through the introduced k-means lustering method. Based on the simplification through k-means clustering, three different candidate location strategies can be obtained. Figure 7 presents the location strategies, and Table 5 shows the optimization results of each strategy.Figure 7
Location strategies obtained by clustering results. (a) Location strategy of {DC3, DC4, DC5}, (b) location strategy of {DC3, DC4, DC5, DC6}, and (c) location strategy of {DC2, DC3, DC4, DC5, DC6}.
(a)(b)(c)Table 5
Result comparison among different location strategies.
Location strategyFixed cost ($)Transportation cost ($)Distribution cost ($)Penalty cost ($)Total cost ($)Number of vehicles{DC3, DC4, DC5}64017631924040487{DC3, DC4, DC5, DC6}88020124841635815{DC2, DC3, DC4, DC5, DC6}10702352600840135Figure7 shows that the selected location of DCs for each location strategy is marked with red circles, and the customers of each selected DC are marked with a black polygon. Table 6 compares the optimization results of these three location strategies to determine the best location strategy, that is, the strategy with the minimum operating cost and maximum utilization of transportation resource. In Table 5, the location strategy of {DC3, DC4, DC5, DC6} requires $3581 in operating cost and achieves five delivery vehicles for distribution. Compared with the two other strategies, the operating cost decreases by 467 and 432, while the number of vehicles increases by 2. Therefore, the location strategy of {DC3, DC4, DC5, DC6} can be determined as the best strategy for the 2E-LRPTWTRS optimization. Then, Table 7 presents the two-echelon vehicle routes optimized by the proposed EMOPSO algorithm.Table 6
Optimal distribution routes and shared vehicles.
Shared vehiclesOptimal distribution routesTime windowV1DC4 ⟶ C140 ⟶ C102 ⟶ C96 ⟶ C90 ⟶ C135 ⟶ DC4[14.4, 18.4]DC4 ⟶ C46 ⟶ C34 ⟶ C97 ⟶ C100 ⟶ DC4[7, 11.2]]DC4 ⟶ C18 ⟶ C12 ⟶ C24 ⟶ C25 ⟶ DC4[4.9, 6.4]DC4 ⟶ C65 ⟶ C103 ⟶ C77 ⟶ C87 ⟶ DC4[11.4, 13.9]DC4 ⟶ C15 ⟶ C49 ⟶ C11 ⟶ C40 ⟶ C58 ⟶ DC4[1.7, 3.2]DC3 ⟶ C45 ⟶ C37 ⟶ C5⟶ C20 ⟶ C3 ⟶ DC3[3.2, 4.7]DC3 ⟶ C71 ⟶ C148 ⟶ C79 ⟶ C113 ⟶ DC3[10.7, 14.6]DC3 ⟶ C7 ⟶ C44 ⟶ C23 ⟶ C51 ⟶ DC3[0.9, 4.2]V2DC4 ⟶ C80 ⟶ C89 ⟶ C105 ⟶ C127 ⟶ DC4[11.7, 13.7]DC4 ⟶ C8 ⟶ C52 ⟶ C60 ⟶ C2 ⟶ C61 ⟶ DC4[2.6, 5.1]DC4 ⟶ C119 ⟶ C116 ⟶ C67 ⟶ C125 ⟶ DC4[15.4, 19]DC3 ⟶ C132 ⟶ C81 ⟶ C101 ⟶ C108 ⟶ C146 ⟶ DC3[13, 17.5]DC3 ⟶ C4 ⟶ C21 ⟶ C42 ⟶ C29 ⟶ C59 ⟶ DC3[0.7, 5.6]DC3 ⟶ C9 ⟶ C32 ⟶ C43 ⟶ C73 ⟶ C86 ⟶ DC3[5.8, 12.6]V3DC4 ⟶ C38 ⟶ C33 ⟶ C14 ⟶ C31 ⟶ DC4[1.2, 3.3]DC4 ⟶ C53 ⟶ C6 ⟶ C92 ⟶ C114 ⟶ C143 ⟶ C138 ⟶ DC4[4.4, 17.5]V4DC3 ⟶ C26 ⟶ C16 ⟶ C57 ⟶ C56 ⟶ C50 ⟶ DC3[1.5, 4.2]DC3 ⟶ C84 ⟶ C78 ⟶ C70 ⟶ C149 ⟶ DC3[11.3, 13]DC3 ⟶ C88⟶ C150 ⟶ C109 ⟶ C122 ⟶ C82 ⟶ DC3[14.2, 18.2]DC6 ⟶ C10 ⟶ C13 ⟶ C39 ⟶ C124 ⟶ C147 ⟶ DC6[5.6, 12.4]DC6 ⟶ C27 ⟶ C1 ⟶ C48 ⟶ C55 ⟶ C62 ⟶ DC6[1.3, 5.1]DC6 ⟶ C93 ⟶ C72 ⟶ C99 ⟶ C120 ⟶ DC6[15.8, 18.5]DC5 ⟶ C111 ⟶ C110 ⟶ C117 ⟶ C76 ⟶ C68 ⟶ C134 ⟶ DC5[10.8, 18.6]DC5 ⟶ C41 ⟶ C28 ⟶ C54 ⟶ C36 ⟶ C19 ⟶ DC51[2.4, 4.5]DC5 ⟶ C47 ⟶ C8 ⟶ DC5[1.2, 2.1]V5DC3 ⟶ C129 ⟶ C145 ⟶ C91 ⟶ C69 ⟶ DC3[11.6, 16.1]DC3 ⟶ C85 ⟶ C112 ⟶ C118 ⟶ C121 ⟶ C123 ⟶ DC3[11.4, 13.3]DC3 ⟶ C95 ⟶ C128 ⟶ C106 ⟶ C139 ⟶ C104 ⟶ C100 ⟶ DC3[11.7, 15.8]DC6 ⟶ C17 ⟶ C63 ⟶ C83 ⟶ C107 ⟶ C64 ⟶ DC6[2, 15.8]DC6 ⟶ C115 ⟶ C142 ⟶ C75 ⟶ C144 ⟶ C130 ⟶ C94 ⟶ DC6[10.9, 13.4]DC5 ⟶ C30 ⟶ C35 ⟶ C137 ⟶ C131 ⟶ C66 ⟶ DC5[3.8, 15.2]DC5 ⟶ C98 ⟶ C141 ⟶ C126 ⟶ C136 ⟶ C74 ⟶ DC5[13.2, 14.1]Table 7
Comparison before and after 2E-LRPTWTRS optimization.
Location strategiesFixed cost ($)Transportation cost ($)Distribution cost ($)Penalty cost ($)Total cost ($)Number of vehiclesBefore optimization{DC1, DC2, DC3, DC4, DC5, DC6}1270387402449573031After optimization{DC3,DC4, DC5, DC6}88020124841635815Table7 shows that 32 distribution routes after 2E-LRPTWTRS optimization while the optimized number of delivery vehicles is 5 because of transportation resource sharing. Transportation resource sharing denotes each vehicle can execute multiple routes of nonoverlapping time windows. For example, seven distribution routes are served by delivery V5 after optimization in Table 8. According by the described sharing, a substantial reduction in transportation resource can be achieved in 2E-LRPTWTRS optimization. Table 6 and Figure 8 show the result comparison before and after 2E-LRPTWTRS optimization.Table 8
Comparison among different vehicle-sharing scenarios.
CaseFixed cost ($)Transportation cost ($)Distribution cost ($)Penalty cost ($)Total cost ($)Number of vehicles1880301309719429711288025630011641539388020124841635815Figure 8
Result comparison before and after 2E-LRPTWTRS optimization.With all existing DCs selected for serving customers in the initial logistics network, the inefficiency of the logistics network is reflected by the high operating cost and low utilization of transportation resource. In the logistics network after 2E-LRPTWTRS optimization, DC3, DC4, DC5, and DC6 are selected to serve customers, and a substantial improvement in the efficiency of the network is obtained. First, the optimized total operating cost is $3581, which saves by $2149 compared with the initial operating cost. Second, the required number of delivery vehicles is 5, which saves 26 compared with the initial required vehicles. Therefore, the proposed 2E-LRPTWTRS optimization methodologies are effective for two-echelon logistics network optimization.


#### 5.1.3. Related Analysis and Discussion

Three cases are considered to evaluate and determine the optimal sharing scheme, prove the contribution of transportation resource sharing, and guarantee a good performance of the sharing scheme simultaneously in the 2E-LRPTWTRS optimization. (1) Transportation resource-sharing strategy is not adopted in the optimization. (2) Transportation resource-sharing strategy is proposed while vehicles can be shared just within the same logistics facilities. (3) Vehicles can be shared within and among different logistics facilities simultaneously in the optimization process. Table8 and Figure 9 show the comparison of the optimization results.Figure 9
Result comparison among different vehicle-sharing scenarios.Table8 and Figure 9 show the optimization results in terms of different vehicle-sharing scenarios. Among the three scenarios, the transportation cost, distribution cost, penalty cost, and numbers of vehicles can obtain the minimum values in the third scenario, which can contribute to the lowest total operating cost in the 2E-LRPTWTRS optimization. Compared with the two other scenarios, the total cost of the third scenario separately decreases by $716 and $572, and the required number of vehicles separately reduces by 6 and 4. Therefore, sharing vehicles within and among different facilities can promote the effectiveness of 2E-LRPTWTRS optimization.


### 5.2. Algorithm Comparison

Following the detailed procedure of the proposed EMOPSO algorithm elaborated in Section5, an algorithm comparison is conducted to demonstrate the superiority of the proposed algorithm in this section. Multiobjective genetic algorithm (MOGA) [2] and nondominated sorting genetic algorithm-II (NSGA-II) [65, 66] are thus selected for the algorithm comparison owing to the good performance and widespread introduction in optimizing logistics network. The comparison is conducted with 20 MDVRPTW instances designed by Cordeau (https://neo.lcc.uma.es/vrp/vrp-instances/multiple-depot-vrp-with-time-windows-instances/), and Table 9 shows the involved characteristics of the instances.Table 9
Description of instances.
InstanceDatasetsNumber of DCsNumber of customersVehicle capacity1Pr-014482002Pr-024961953Pr-0341441904Pr-0441921855Pr-0542401806Pr-0642881757Pr-076722008Pr-0861441909Pr-09621618010Pr-10628817011Pr-1144820012Pr-1249619513Pr-13414419014Pr-14419218515Pr-15424018016Pr-16428817517Pr-1767220018Pr-18614419019Pr-19621618020Pr-206288170Table9 shows that the 20 datasets are different from one another in the number of DCs, customers, and vehicle capacity. The three algorithms, namely, MOGA, NSGA-II, and the proposed EMOPSO, are used to calculate the optimal total operating cost (TC), the required vehicles (NV), and the waiting time (T) of each instance. The involved parameters of the algorithm comparison are set as follows: population size popsize = 150, maximum number of generations genmax = 300, crossover probability crosp = 0.9, and mutation probability mutp = 0.1 in MOGA and NSGA-II; maximum iteration max_Iter = 300, particle size nP = 150, the size of external repository nR = 40, the count of grids per dimension nG = 5, initial particle inertia weight w = 1, and initial particle personal learning coefficient c1 = 1, and the initial particle global learning coefficient c2 = 2 in EMOPSO. The calculation results provided by the three algorithms are presented in Table 10.Table 10
Optimization result comparison by each algorithm.
InstanceEMOPSONSGA-IIMOGATC ($)NVT (min)TC ($)NVT (min)TC ($)NVT (min)1157661617648241656723218846172182727199772932099822216110322228113042357131825021421267315225265815192777162327151721630391518326117313200163372854132328321334290114308334416263469174034251639940011735434618374169173910436719264601192845382130111695627196172618748241219907222167734200593313229482125419222479102514248911182629122726531426152896141331311429310116281631191714321619263354202417284915192969162929981735183468162134851727357419241941061621445818344210193620465819284919203847962139Average288713213069142930271530t-test9.3E − 087.9E − 072.0E − 071.4E − 081.8E − 084.8E − 07In Table10, the calculation of t-test shows a substantial difference among the optimal results calculated by the algorithms and proves the reasonability of the comparison among the three algorithms. By comparing the average values of TC, NV, and WT, the results provided by EMOPSO are all superior to those provided by MOGA and NSGA-II. First, the average operating cost in EMOPSO is $2887, which saves by $182 compared with MOGA and $140 compared with NSGA-II. Second, compared with 14 required vehicles in MOGA and 15 required vehicles in NSGA-II, the average required number of vehicles in EMOPSO is 13. Third, EMOPSO obtains the minimum average value of waiting time as 21, compared with 29 in MOGA and 30 in NSGA-II. Therefore, the proposed EMOPSO is more effective than the two other algorithms.


### 5.3. Management Insights

The design of a two-echelon logistics network is crucial for modern supply chain management because the logistics operating cost constitutes the majority of the expenses of companies. In this study, the proposed 2E-LRPTWTRS optimization provides an effective method with respect to tackling the facility and routing decision simultaneously, and the introduction of transportation resource sharing considerably improves the utilization of delivery vehicles. The management insights concluded are as follows:(1)
An effective location strategy consists of determining the reasonable number and geographical location of opened logistics facilities. Compared with the initial logistics network, the location strategy provided by the proposed 2E-LRPTWTRS optimization enables the customers of a distribution area to be served by fewer DCs with a higher distribution efficiency. In addition, the vehicle-routing optimization based on the determined location strategy can design a more well-organized distribution routes for the two-echelon logistics network, which greatly reduces a series of unreasonable distribution phenomenon in the initial logistics network. Therefore, the proposed 2E-LRPTWTRS optimization can provide logistics enterprises an effective reference for tackling the location decisions and vehicle-routing scheme simultaneously to guarantee the efficient operations of the two-echelon logistics network.(2)
The introduction of transportation resource-sharing strategy in the proposed 2E-LRPTWTRS optimization plays an important role in reducing logistics cost for a two-echelon logistics network. Different distribution routes have varied starting and ending times because of the time window constraints. Different from the phenomenon that each distribution route requires one delivery vehicle in the traditional distribution, a delivery vehicle can be shared within and among different logistics facilities and perform multiple distribution routes only if the time window allows, which can remarkably reduce the number of required vehicles in the two-echelon logistics network. This strategy not only satisfies the requirement of logistics enterprises in reducing operating cost but also responds to the call of our government to promote sustainable development, which is demonstrated as a promising strategy to be encouraged.


## 5.1. Case Study

A real case of 2E-LRPTWTRS in Chongqing, China, is employed to present the effectiveness of the proposed optimization model and hybrid algorithm. The logistics network composed of one LC, six candidate DCs (DC1, DC2, …, DC6), and 150 customers (C1, C2, …, C150) in this selected distribution area is shown in Figure5. Here, LC is represented as pentagon, stars, squares, crosses, triangles, and circles, and diamonds refer to the DCs and their corresponding service customers.Figure 5
Spatial distribution of logistics facilities and customers.


### 5.1.1. Data Description

According to actual surveys and related references [2, 17, 37], the involved parameter values used in the proposed optimization model and EMOPSO can be summarized as follows: capacity of the truck Cs = 1000, capacity of the vehicle Cv = 200, fuel consumption of the truck fs = 0.25, fuel consumption of the vehicle fv = 0.2, diesel price Ps = 20, gasoline price Pv = 15, maintenance cost of the truck Ms = 1500, maintenance cost of the vehicle Mv = 500, penalty cost per time unit for earliness λe = 10, penalty cost per time unit for delay λl = 25, maximum iteration max_Iter = 600, particle size nP = 150, the size of external repository nR = 75, the count of grids per dimension nG = 7, initial particle inertia weight w = 1, initial particle personal learning coefficient c1 = 1, and initial particle global learning coefficient c2 = 2. Moreover, Tables 3 and 4 show the initial characteristics of logistics facilities and initial distribution routes, respectively.Table 3
Initial characteristics of logistics facilities.
FacilitiesSymbolFixed cost ($)Number of service customersTime windowsLC——[0, 22]DC120026[1, 19]DC219024[0.5, 17.5]DC322025[0, 18]DC425026[2.5, 21]DC517025[2, 18]DC624024[1.5, 18]Total1270150Table 4
Initial distribution routes of logistics network.
Distribution centersOptimal distribution routesNumber of vehiclesDC1DC1 ⟶ C128 ⟶ C32 ⟶ C113 ⟶ C123 ⟶ C106 ⟶ DC1DC1 ⟶ C108 ⟶ C146 ⟶ C114 ⟶ C143 ⟶ C23 ⟶ DC1DC1 ⟶ C88 ⟶ C145 ⟶ C59 ⟶ C9 ⟶ C7 ⟶ DC1DC1 ⟶ C71 ⟶ C73 ⟶ C82 ⟶ C17 ⟶ DC1DC1 ⟶ C4 ⟶ C69 ⟶ C116 ⟶ C3 ⟶ C86 ⟶ DC131DC2DC2 ⟶ C38 ⟶ C12 ⟶ C34 ⟶ C97 ⟶ DC2DC2 ⟶ C62 ⟶ C142 ⟶ C115 ⟶ C8 ⟶ C14 ⟶ DC2DC2 ⟶ C148 ⟶ C57 ⟶ C91 ⟶ C84 ⟶ DC2DC2 ⟶ C21 ⟶ C45 ⟶ C37 ⟶ C129 ⟶ DC2DC2 ⟶ C42 ⟶ C122 ⟶ C56 ⟶ C29 ⟶ C44 ⟶ DC2DC2 ⟶ C51 ⟶ C121 ⟶ C112 ⟶ C26 ⟶ C139 ⟶ DC2DC3DC3 ⟶ C50 ⟶ C104 ⟶ C79 ⟶ C70 ⟶ DC3DC3 ⟶ C149 ⟶ C80 ⟶ C43 ⟶ C63 ⟶ C118 ⟶ C132 ⟶ DC3DC3 ⟶ C133 ⟶ C101 ⟶ C81 ⟶ C92 ⟶ DC3DC3 ⟶ C95 ⟶ C109 ⟶ C120 ⟶ C2 ⟶ DC3DC3 ⟶ C20 ⟶ C5 ⟶ C72 ⟶ DC3DC3 ⟶ C24 ⟶ C100 ⟶ C22 ⟶ C105 ⟶ DC3DC4DC4 ⟶ C39 ⟶ C64 ⟶ C16 ⟶ C150 ⟶ C78 ⟶ DC4DC4 ⟶ C58 ⟶ C138 ⟶ C135 ⟶ C130 ⟶ DC4DC4 ⟶ C80⟶ C52 ⟶ C127 ⟶ C96 ⟶ C94 ⟶ DC4DC4 ⟶ C34 ⟶ C49 ⟶ C119 ⟶ C33 ⟶ C15 ⟶ DC4DC4 ⟶ C25 ⟶ C65 ⟶ C53 ⟶ C89 ⟶ DC4DC5DC5 ⟶ C48 ⟶ C75 ⟶ C99 ⟶ C55 ⟶ C107 ⟶ DC5DC5 ⟶ C41 ⟶ C141 ⟶ C19 ⟶ C30 ⟶ C28 ⟶ DC5DC5 ⟶ C6 ⟶ C61 ⟶ C67 ⟶ C66 ⟶ DC5DC5 ⟶ C74 ⟶ C134 ⟶ C137 ⟶ C90 ⟶ CDC ⟶ DC5DC5 ⟶ C140 ⟶ C40 ⟶ C11 ⟶ C46 ⟶ C18 ⟶ DC5DC6DC6 ⟶ C10 ⟶ C124 ⟶ C27 ⟶ C13 ⟶ C111 ⟶ DC6DC6 ⟶ C144 ⟶ C147 ⟶ C1 ⟶ C98 ⟶ C68 ⟶ C110 ⟶ DC6DC6 ⟶ C117 ⟶ C47 ⟶ C126 ⟶ C131 ⟶ C54 ⟶ DC6DC6 ⟶ C125 ⟶ C60 ⟶ C103 ⟶ C35 ⟶ C76 ⟶ DC6Tables3 and 4 show that, with all the existing DCs selected for serving customers, each DC serves a few customers and results in a low utilization of logistics facilities. Moreover, the delivery routes without vehicle sharing can lead to a large number of delivery vehicles for the two-echelon logistics network. Therefore, determining the reasonable location strategy and routing scheme through 2E-LRPTWTRS optimization to reduce the total cost is important for optimizing the logistics network.


### 5.1.2. Optimization Results

The optimization of 2E-LRPTWTRS can achieve a substantial simplification throughk-means clustering. Before the operation of the proposed EMOPSO to optimize the vehicle routes, an effective customer clustering can avoid searching several unnecessary solutions and increasing the computation time, which seriously improves the efficiency of the algorithm. Therefore, to avoid such unreasonable phenomenon, customers are clustered through k-means clustering. In customer clustering, the main challenge is the determination of the k value. Silhouette analysis is introduced to evaluate the performance of customer clustering for each possible k value and address this challenge [2, 4]. The data description shows that 150 customers are in this studied case, and the set of possible k values for the case is reasonably determined as {3, 4, 5}. Figure 6 shows the customer clustering results and silhouette values in terms of possible k values.Figure 6
Customer (a) clustering results and (b) silhouette values.
(a)(b)Figure6(a) shows that three different clustering results are obtained because of the three different possible k values, and the circles with the same color are assigned to the same cluster. Moreover, one center is in each cluster, and the customers of the same cluster are with shorter distance to the corresponding center compared with the other centers. In Figure 6(b), the silhouette values are provided for evaluating the clustering performance, and the majority of the elements have a positive value while few negative values are observed for each clustering scenario, which denotes that a good clustering performance can be obtained through the introduced k-means lustering method. Based on the simplification through k-means clustering, three different candidate location strategies can be obtained. Figure 7 presents the location strategies, and Table 5 shows the optimization results of each strategy.Figure 7
Location strategies obtained by clustering results. (a) Location strategy of {DC3, DC4, DC5}, (b) location strategy of {DC3, DC4, DC5, DC6}, and (c) location strategy of {DC2, DC3, DC4, DC5, DC6}.
(a)(b)(c)Table 5
Result comparison among different location strategies.
Location strategyFixed cost ($)Transportation cost ($)Distribution cost ($)Penalty cost ($)Total cost ($)Number of vehicles{DC3, DC4, DC5}64017631924040487{DC3, DC4, DC5, DC6}88020124841635815{DC2, DC3, DC4, DC5, DC6}10702352600840135Figure7 shows that the selected location of DCs for each location strategy is marked with red circles, and the customers of each selected DC are marked with a black polygon. Table 6 compares the optimization results of these three location strategies to determine the best location strategy, that is, the strategy with the minimum operating cost and maximum utilization of transportation resource. In Table 5, the location strategy of {DC3, DC4, DC5, DC6} requires $3581 in operating cost and achieves five delivery vehicles for distribution. Compared with the two other strategies, the operating cost decreases by 467 and 432, while the number of vehicles increases by 2. Therefore, the location strategy of {DC3, DC4, DC5, DC6} can be determined as the best strategy for the 2E-LRPTWTRS optimization. Then, Table 7 presents the two-echelon vehicle routes optimized by the proposed EMOPSO algorithm.Table 6
Optimal distribution routes and shared vehicles.
Shared vehiclesOptimal distribution routesTime windowV1DC4 ⟶ C140 ⟶ C102 ⟶ C96 ⟶ C90 ⟶ C135 ⟶ DC4[14.4, 18.4]DC4 ⟶ C46 ⟶ C34 ⟶ C97 ⟶ C100 ⟶ DC4[7, 11.2]]DC4 ⟶ C18 ⟶ C12 ⟶ C24 ⟶ C25 ⟶ DC4[4.9, 6.4]DC4 ⟶ C65 ⟶ C103 ⟶ C77 ⟶ C87 ⟶ DC4[11.4, 13.9]DC4 ⟶ C15 ⟶ C49 ⟶ C11 ⟶ C40 ⟶ C58 ⟶ DC4[1.7, 3.2]DC3 ⟶ C45 ⟶ C37 ⟶ C5⟶ C20 ⟶ C3 ⟶ DC3[3.2, 4.7]DC3 ⟶ C71 ⟶ C148 ⟶ C79 ⟶ C113 ⟶ DC3[10.7, 14.6]DC3 ⟶ C7 ⟶ C44 ⟶ C23 ⟶ C51 ⟶ DC3[0.9, 4.2]V2DC4 ⟶ C80 ⟶ C89 ⟶ C105 ⟶ C127 ⟶ DC4[11.7, 13.7]DC4 ⟶ C8 ⟶ C52 ⟶ C60 ⟶ C2 ⟶ C61 ⟶ DC4[2.6, 5.1]DC4 ⟶ C119 ⟶ C116 ⟶ C67 ⟶ C125 ⟶ DC4[15.4, 19]DC3 ⟶ C132 ⟶ C81 ⟶ C101 ⟶ C108 ⟶ C146 ⟶ DC3[13, 17.5]DC3 ⟶ C4 ⟶ C21 ⟶ C42 ⟶ C29 ⟶ C59 ⟶ DC3[0.7, 5.6]DC3 ⟶ C9 ⟶ C32 ⟶ C43 ⟶ C73 ⟶ C86 ⟶ DC3[5.8, 12.6]V3DC4 ⟶ C38 ⟶ C33 ⟶ C14 ⟶ C31 ⟶ DC4[1.2, 3.3]DC4 ⟶ C53 ⟶ C6 ⟶ C92 ⟶ C114 ⟶ C143 ⟶ C138 ⟶ DC4[4.4, 17.5]V4DC3 ⟶ C26 ⟶ C16 ⟶ C57 ⟶ C56 ⟶ C50 ⟶ DC3[1.5, 4.2]DC3 ⟶ C84 ⟶ C78 ⟶ C70 ⟶ C149 ⟶ DC3[11.3, 13]DC3 ⟶ C88⟶ C150 ⟶ C109 ⟶ C122 ⟶ C82 ⟶ DC3[14.2, 18.2]DC6 ⟶ C10 ⟶ C13 ⟶ C39 ⟶ C124 ⟶ C147 ⟶ DC6[5.6, 12.4]DC6 ⟶ C27 ⟶ C1 ⟶ C48 ⟶ C55 ⟶ C62 ⟶ DC6[1.3, 5.1]DC6 ⟶ C93 ⟶ C72 ⟶ C99 ⟶ C120 ⟶ DC6[15.8, 18.5]DC5 ⟶ C111 ⟶ C110 ⟶ C117 ⟶ C76 ⟶ C68 ⟶ C134 ⟶ DC5[10.8, 18.6]DC5 ⟶ C41 ⟶ C28 ⟶ C54 ⟶ C36 ⟶ C19 ⟶ DC51[2.4, 4.5]DC5 ⟶ C47 ⟶ C8 ⟶ DC5[1.2, 2.1]V5DC3 ⟶ C129 ⟶ C145 ⟶ C91 ⟶ C69 ⟶ DC3[11.6, 16.1]DC3 ⟶ C85 ⟶ C112 ⟶ C118 ⟶ C121 ⟶ C123 ⟶ DC3[11.4, 13.3]DC3 ⟶ C95 ⟶ C128 ⟶ C106 ⟶ C139 ⟶ C104 ⟶ C100 ⟶ DC3[11.7, 15.8]DC6 ⟶ C17 ⟶ C63 ⟶ C83 ⟶ C107 ⟶ C64 ⟶ DC6[2, 15.8]DC6 ⟶ C115 ⟶ C142 ⟶ C75 ⟶ C144 ⟶ C130 ⟶ C94 ⟶ DC6[10.9, 13.4]DC5 ⟶ C30 ⟶ C35 ⟶ C137 ⟶ C131 ⟶ C66 ⟶ DC5[3.8, 15.2]DC5 ⟶ C98 ⟶ C141 ⟶ C126 ⟶ C136 ⟶ C74 ⟶ DC5[13.2, 14.1]Table 7
Comparison before and after 2E-LRPTWTRS optimization.
Location strategiesFixed cost ($)Transportation cost ($)Distribution cost ($)Penalty cost ($)Total cost ($)Number of vehiclesBefore optimization{DC1, DC2, DC3, DC4, DC5, DC6}1270387402449573031After optimization{DC3,DC4, DC5, DC6}88020124841635815Table7 shows that 32 distribution routes after 2E-LRPTWTRS optimization while the optimized number of delivery vehicles is 5 because of transportation resource sharing. Transportation resource sharing denotes each vehicle can execute multiple routes of nonoverlapping time windows. For example, seven distribution routes are served by delivery V5 after optimization in Table 8. According by the described sharing, a substantial reduction in transportation resource can be achieved in 2E-LRPTWTRS optimization. Table 6 and Figure 8 show the result comparison before and after 2E-LRPTWTRS optimization.Table 8
Comparison among different vehicle-sharing scenarios.
CaseFixed cost ($)Transportation cost ($)Distribution cost ($)Penalty cost ($)Total cost ($)Number of vehicles1880301309719429711288025630011641539388020124841635815Figure 8
Result comparison before and after 2E-LRPTWTRS optimization.With all existing DCs selected for serving customers in the initial logistics network, the inefficiency of the logistics network is reflected by the high operating cost and low utilization of transportation resource. In the logistics network after 2E-LRPTWTRS optimization, DC3, DC4, DC5, and DC6 are selected to serve customers, and a substantial improvement in the efficiency of the network is obtained. First, the optimized total operating cost is $3581, which saves by $2149 compared with the initial operating cost. Second, the required number of delivery vehicles is 5, which saves 26 compared with the initial required vehicles. Therefore, the proposed 2E-LRPTWTRS optimization methodologies are effective for two-echelon logistics network optimization.


### 5.1.3. Related Analysis and Discussion

Three cases are considered to evaluate and determine the optimal sharing scheme, prove the contribution of transportation resource sharing, and guarantee a good performance of the sharing scheme simultaneously in the 2E-LRPTWTRS optimization. (1) Transportation resource-sharing strategy is not adopted in the optimization. (2) Transportation resource-sharing strategy is proposed while vehicles can be shared just within the same logistics facilities. (3) Vehicles can be shared within and among different logistics facilities simultaneously in the optimization process. Table8 and Figure 9 show the comparison of the optimization results.Figure 9
Result comparison among different vehicle-sharing scenarios.Table8 and Figure 9 show the optimization results in terms of different vehicle-sharing scenarios. Among the three scenarios, the transportation cost, distribution cost, penalty cost, and numbers of vehicles can obtain the minimum values in the third scenario, which can contribute to the lowest total operating cost in the 2E-LRPTWTRS optimization. Compared with the two other scenarios, the total cost of the third scenario separately decreases by $716 and $572, and the required number of vehicles separately reduces by 6 and 4. Therefore, sharing vehicles within and among different facilities can promote the effectiveness of 2E-LRPTWTRS optimization.


## 5.1.1. Data Description

According to actual surveys and related references [2, 17, 37], the involved parameter values used in the proposed optimization model and EMOPSO can be summarized as follows: capacity of the truck Cs = 1000, capacity of the vehicle Cv = 200, fuel consumption of the truck fs = 0.25, fuel consumption of the vehicle fv = 0.2, diesel price Ps = 20, gasoline price Pv = 15, maintenance cost of the truck Ms = 1500, maintenance cost of the vehicle Mv = 500, penalty cost per time unit for earliness λe = 10, penalty cost per time unit for delay λl = 25, maximum iteration max_Iter = 600, particle size nP = 150, the size of external repository nR = 75, the count of grids per dimension nG = 7, initial particle inertia weight w = 1, initial particle personal learning coefficient c1 = 1, and initial particle global learning coefficient c2 = 2. Moreover, Tables 3 and 4 show the initial characteristics of logistics facilities and initial distribution routes, respectively.Table 3
Initial characteristics of logistics facilities.
FacilitiesSymbolFixed cost ($)Number of service customersTime windowsLC——[0, 22]DC120026[1, 19]DC219024[0.5, 17.5]DC322025[0, 18]DC425026[2.5, 21]DC517025[2, 18]DC624024[1.5, 18]Total1270150Table 4
Initial distribution routes of logistics network.
Distribution centersOptimal distribution routesNumber of vehiclesDC1DC1 ⟶ C128 ⟶ C32 ⟶ C113 ⟶ C123 ⟶ C106 ⟶ DC1DC1 ⟶ C108 ⟶ C146 ⟶ C114 ⟶ C143 ⟶ C23 ⟶ DC1DC1 ⟶ C88 ⟶ C145 ⟶ C59 ⟶ C9 ⟶ C7 ⟶ DC1DC1 ⟶ C71 ⟶ C73 ⟶ C82 ⟶ C17 ⟶ DC1DC1 ⟶ C4 ⟶ C69 ⟶ C116 ⟶ C3 ⟶ C86 ⟶ DC131DC2DC2 ⟶ C38 ⟶ C12 ⟶ C34 ⟶ C97 ⟶ DC2DC2 ⟶ C62 ⟶ C142 ⟶ C115 ⟶ C8 ⟶ C14 ⟶ DC2DC2 ⟶ C148 ⟶ C57 ⟶ C91 ⟶ C84 ⟶ DC2DC2 ⟶ C21 ⟶ C45 ⟶ C37 ⟶ C129 ⟶ DC2DC2 ⟶ C42 ⟶ C122 ⟶ C56 ⟶ C29 ⟶ C44 ⟶ DC2DC2 ⟶ C51 ⟶ C121 ⟶ C112 ⟶ C26 ⟶ C139 ⟶ DC2DC3DC3 ⟶ C50 ⟶ C104 ⟶ C79 ⟶ C70 ⟶ DC3DC3 ⟶ C149 ⟶ C80 ⟶ C43 ⟶ C63 ⟶ C118 ⟶ C132 ⟶ DC3DC3 ⟶ C133 ⟶ C101 ⟶ C81 ⟶ C92 ⟶ DC3DC3 ⟶ C95 ⟶ C109 ⟶ C120 ⟶ C2 ⟶ DC3DC3 ⟶ C20 ⟶ C5 ⟶ C72 ⟶ DC3DC3 ⟶ C24 ⟶ C100 ⟶ C22 ⟶ C105 ⟶ DC3DC4DC4 ⟶ C39 ⟶ C64 ⟶ C16 ⟶ C150 ⟶ C78 ⟶ DC4DC4 ⟶ C58 ⟶ C138 ⟶ C135 ⟶ C130 ⟶ DC4DC4 ⟶ C80⟶ C52 ⟶ C127 ⟶ C96 ⟶ C94 ⟶ DC4DC4 ⟶ C34 ⟶ C49 ⟶ C119 ⟶ C33 ⟶ C15 ⟶ DC4DC4 ⟶ C25 ⟶ C65 ⟶ C53 ⟶ C89 ⟶ DC4DC5DC5 ⟶ C48 ⟶ C75 ⟶ C99 ⟶ C55 ⟶ C107 ⟶ DC5DC5 ⟶ C41 ⟶ C141 ⟶ C19 ⟶ C30 ⟶ C28 ⟶ DC5DC5 ⟶ C6 ⟶ C61 ⟶ C67 ⟶ C66 ⟶ DC5DC5 ⟶ C74 ⟶ C134 ⟶ C137 ⟶ C90 ⟶ CDC ⟶ DC5DC5 ⟶ C140 ⟶ C40 ⟶ C11 ⟶ C46 ⟶ C18 ⟶ DC5DC6DC6 ⟶ C10 ⟶ C124 ⟶ C27 ⟶ C13 ⟶ C111 ⟶ DC6DC6 ⟶ C144 ⟶ C147 ⟶ C1 ⟶ C98 ⟶ C68 ⟶ C110 ⟶ DC6DC6 ⟶ C117 ⟶ C47 ⟶ C126 ⟶ C131 ⟶ C54 ⟶ DC6DC6 ⟶ C125 ⟶ C60 ⟶ C103 ⟶ C35 ⟶ C76 ⟶ DC6Tables3 and 4 show that, with all the existing DCs selected for serving customers, each DC serves a few customers and results in a low utilization of logistics facilities. Moreover, the delivery routes without vehicle sharing can lead to a large number of delivery vehicles for the two-echelon logistics network. Therefore, determining the reasonable location strategy and routing scheme through 2E-LRPTWTRS optimization to reduce the total cost is important for optimizing the logistics network.


## 5.1.2. Optimization Results

The optimization of 2E-LRPTWTRS can achieve a substantial simplification throughk-means clustering. Before the operation of the proposed EMOPSO to optimize the vehicle routes, an effective customer clustering can avoid searching several unnecessary solutions and increasing the computation time, which seriously improves the efficiency of the algorithm. Therefore, to avoid such unreasonable phenomenon, customers are clustered through k-means clustering. In customer clustering, the main challenge is the determination of the k value. Silhouette analysis is introduced to evaluate the performance of customer clustering for each possible k value and address this challenge [2, 4]. The data description shows that 150 customers are in this studied case, and the set of possible k values for the case is reasonably determined as {3, 4, 5}. Figure 6 shows the customer clustering results and silhouette values in terms of possible k values.Figure 6
Customer (a) clustering results and (b) silhouette values.
(a)(b)Figure6(a) shows that three different clustering results are obtained because of the three different possible k values, and the circles with the same color are assigned to the same cluster. Moreover, one center is in each cluster, and the customers of the same cluster are with shorter distance to the corresponding center compared with the other centers. In Figure 6(b), the silhouette values are provided for evaluating the clustering performance, and the majority of the elements have a positive value while few negative values are observed for each clustering scenario, which denotes that a good clustering performance can be obtained through the introduced k-means lustering method. Based on the simplification through k-means clustering, three different candidate location strategies can be obtained. Figure 7 presents the location strategies, and Table 5 shows the optimization results of each strategy.Figure 7
Location strategies obtained by clustering results. (a) Location strategy of {DC3, DC4, DC5}, (b) location strategy of {DC3, DC4, DC5, DC6}, and (c) location strategy of {DC2, DC3, DC4, DC5, DC6}.
(a)(b)(c)Table 5
Result comparison among different location strategies.
Location strategyFixed cost ($)Transportation cost ($)Distribution cost ($)Penalty cost ($)Total cost ($)Number of vehicles{DC3, DC4, DC5}64017631924040487{DC3, DC4, DC5, DC6}88020124841635815{DC2, DC3, DC4, DC5, DC6}10702352600840135Figure7 shows that the selected location of DCs for each location strategy is marked with red circles, and the customers of each selected DC are marked with a black polygon. Table 6 compares the optimization results of these three location strategies to determine the best location strategy, that is, the strategy with the minimum operating cost and maximum utilization of transportation resource. In Table 5, the location strategy of {DC3, DC4, DC5, DC6} requires $3581 in operating cost and achieves five delivery vehicles for distribution. Compared with the two other strategies, the operating cost decreases by 467 and 432, while the number of vehicles increases by 2. Therefore, the location strategy of {DC3, DC4, DC5, DC6} can be determined as the best strategy for the 2E-LRPTWTRS optimization. Then, Table 7 presents the two-echelon vehicle routes optimized by the proposed EMOPSO algorithm.Table 6
Optimal distribution routes and shared vehicles.
Shared vehiclesOptimal distribution routesTime windowV1DC4 ⟶ C140 ⟶ C102 ⟶ C96 ⟶ C90 ⟶ C135 ⟶ DC4[14.4, 18.4]DC4 ⟶ C46 ⟶ C34 ⟶ C97 ⟶ C100 ⟶ DC4[7, 11.2]]DC4 ⟶ C18 ⟶ C12 ⟶ C24 ⟶ C25 ⟶ DC4[4.9, 6.4]DC4 ⟶ C65 ⟶ C103 ⟶ C77 ⟶ C87 ⟶ DC4[11.4, 13.9]DC4 ⟶ C15 ⟶ C49 ⟶ C11 ⟶ C40 ⟶ C58 ⟶ DC4[1.7, 3.2]DC3 ⟶ C45 ⟶ C37 ⟶ C5⟶ C20 ⟶ C3 ⟶ DC3[3.2, 4.7]DC3 ⟶ C71 ⟶ C148 ⟶ C79 ⟶ C113 ⟶ DC3[10.7, 14.6]DC3 ⟶ C7 ⟶ C44 ⟶ C23 ⟶ C51 ⟶ DC3[0.9, 4.2]V2DC4 ⟶ C80 ⟶ C89 ⟶ C105 ⟶ C127 ⟶ DC4[11.7, 13.7]DC4 ⟶ C8 ⟶ C52 ⟶ C60 ⟶ C2 ⟶ C61 ⟶ DC4[2.6, 5.1]DC4 ⟶ C119 ⟶ C116 ⟶ C67 ⟶ C125 ⟶ DC4[15.4, 19]DC3 ⟶ C132 ⟶ C81 ⟶ C101 ⟶ C108 ⟶ C146 ⟶ DC3[13, 17.5]DC3 ⟶ C4 ⟶ C21 ⟶ C42 ⟶ C29 ⟶ C59 ⟶ DC3[0.7, 5.6]DC3 ⟶ C9 ⟶ C32 ⟶ C43 ⟶ C73 ⟶ C86 ⟶ DC3[5.8, 12.6]V3DC4 ⟶ C38 ⟶ C33 ⟶ C14 ⟶ C31 ⟶ DC4[1.2, 3.3]DC4 ⟶ C53 ⟶ C6 ⟶ C92 ⟶ C114 ⟶ C143 ⟶ C138 ⟶ DC4[4.4, 17.5]V4DC3 ⟶ C26 ⟶ C16 ⟶ C57 ⟶ C56 ⟶ C50 ⟶ DC3[1.5, 4.2]DC3 ⟶ C84 ⟶ C78 ⟶ C70 ⟶ C149 ⟶ DC3[11.3, 13]DC3 ⟶ C88⟶ C150 ⟶ C109 ⟶ C122 ⟶ C82 ⟶ DC3[14.2, 18.2]DC6 ⟶ C10 ⟶ C13 ⟶ C39 ⟶ C124 ⟶ C147 ⟶ DC6[5.6, 12.4]DC6 ⟶ C27 ⟶ C1 ⟶ C48 ⟶ C55 ⟶ C62 ⟶ DC6[1.3, 5.1]DC6 ⟶ C93 ⟶ C72 ⟶ C99 ⟶ C120 ⟶ DC6[15.8, 18.5]DC5 ⟶ C111 ⟶ C110 ⟶ C117 ⟶ C76 ⟶ C68 ⟶ C134 ⟶ DC5[10.8, 18.6]DC5 ⟶ C41 ⟶ C28 ⟶ C54 ⟶ C36 ⟶ C19 ⟶ DC51[2.4, 4.5]DC5 ⟶ C47 ⟶ C8 ⟶ DC5[1.2, 2.1]V5DC3 ⟶ C129 ⟶ C145 ⟶ C91 ⟶ C69 ⟶ DC3[11.6, 16.1]DC3 ⟶ C85 ⟶ C112 ⟶ C118 ⟶ C121 ⟶ C123 ⟶ DC3[11.4, 13.3]DC3 ⟶ C95 ⟶ C128 ⟶ C106 ⟶ C139 ⟶ C104 ⟶ C100 ⟶ DC3[11.7, 15.8]DC6 ⟶ C17 ⟶ C63 ⟶ C83 ⟶ C107 ⟶ C64 ⟶ DC6[2, 15.8]DC6 ⟶ C115 ⟶ C142 ⟶ C75 ⟶ C144 ⟶ C130 ⟶ C94 ⟶ DC6[10.9, 13.4]DC5 ⟶ C30 ⟶ C35 ⟶ C137 ⟶ C131 ⟶ C66 ⟶ DC5[3.8, 15.2]DC5 ⟶ C98 ⟶ C141 ⟶ C126 ⟶ C136 ⟶ C74 ⟶ DC5[13.2, 14.1]Table 7
Comparison before and after 2E-LRPTWTRS optimization.
Location strategiesFixed cost ($)Transportation cost ($)Distribution cost ($)Penalty cost ($)Total cost ($)Number of vehiclesBefore optimization{DC1, DC2, DC3, DC4, DC5, DC6}1270387402449573031After optimization{DC3,DC4, DC5, DC6}88020124841635815Table7 shows that 32 distribution routes after 2E-LRPTWTRS optimization while the optimized number of delivery vehicles is 5 because of transportation resource sharing. Transportation resource sharing denotes each vehicle can execute multiple routes of nonoverlapping time windows. For example, seven distribution routes are served by delivery V5 after optimization in Table 8. According by the described sharing, a substantial reduction in transportation resource can be achieved in 2E-LRPTWTRS optimization. Table 6 and Figure 8 show the result comparison before and after 2E-LRPTWTRS optimization.Table 8
Comparison among different vehicle-sharing scenarios.
CaseFixed cost ($)Transportation cost ($)Distribution cost ($)Penalty cost ($)Total cost ($)Number of vehicles1880301309719429711288025630011641539388020124841635815Figure 8
Result comparison before and after 2E-LRPTWTRS optimization.With all existing DCs selected for serving customers in the initial logistics network, the inefficiency of the logistics network is reflected by the high operating cost and low utilization of transportation resource. In the logistics network after 2E-LRPTWTRS optimization, DC3, DC4, DC5, and DC6 are selected to serve customers, and a substantial improvement in the efficiency of the network is obtained. First, the optimized total operating cost is $3581, which saves by $2149 compared with the initial operating cost. Second, the required number of delivery vehicles is 5, which saves 26 compared with the initial required vehicles. Therefore, the proposed 2E-LRPTWTRS optimization methodologies are effective for two-echelon logistics network optimization.


## 5.1.3. Related Analysis and Discussion

Three cases are considered to evaluate and determine the optimal sharing scheme, prove the contribution of transportation resource sharing, and guarantee a good performance of the sharing scheme simultaneously in the 2E-LRPTWTRS optimization. (1) Transportation resource-sharing strategy is not adopted in the optimization. (2) Transportation resource-sharing strategy is proposed while vehicles can be shared just within the same logistics facilities. (3) Vehicles can be shared within and among different logistics facilities simultaneously in the optimization process. Table8 and Figure 9 show the comparison of the optimization results.Figure 9
Result comparison among different vehicle-sharing scenarios.Table8 and Figure 9 show the optimization results in terms of different vehicle-sharing scenarios. Among the three scenarios, the transportation cost, distribution cost, penalty cost, and numbers of vehicles can obtain the minimum values in the third scenario, which can contribute to the lowest total operating cost in the 2E-LRPTWTRS optimization. Compared with the two other scenarios, the total cost of the third scenario separately decreases by $716 and $572, and the required number of vehicles separately reduces by 6 and 4. Therefore, sharing vehicles within and among different facilities can promote the effectiveness of 2E-LRPTWTRS optimization.


## 5.2. Algorithm Comparison

Following the detailed procedure of the proposed EMOPSO algorithm elaborated in Section5, an algorithm comparison is conducted to demonstrate the superiority of the proposed algorithm in this section. Multiobjective genetic algorithm (MOGA) [2] and nondominated sorting genetic algorithm-II (NSGA-II) [65, 66] are thus selected for the algorithm comparison owing to the good performance and widespread introduction in optimizing logistics network. The comparison is conducted with 20 MDVRPTW instances designed by Cordeau (https://neo.lcc.uma.es/vrp/vrp-instances/multiple-depot-vrp-with-time-windows-instances/), and Table 9 shows the involved characteristics of the instances.Table 9
Description of instances.
InstanceDatasetsNumber of DCsNumber of customersVehicle capacity1Pr-014482002Pr-024961953Pr-0341441904Pr-0441921855Pr-0542401806Pr-0642881757Pr-076722008Pr-0861441909Pr-09621618010Pr-10628817011Pr-1144820012Pr-1249619513Pr-13414419014Pr-14419218515Pr-15424018016Pr-16428817517Pr-1767220018Pr-18614419019Pr-19621618020Pr-206288170Table9 shows that the 20 datasets are different from one another in the number of DCs, customers, and vehicle capacity. The three algorithms, namely, MOGA, NSGA-II, and the proposed EMOPSO, are used to calculate the optimal total operating cost (TC), the required vehicles (NV), and the waiting time (T) of each instance. The involved parameters of the algorithm comparison are set as follows: population size popsize = 150, maximum number of generations genmax = 300, crossover probability crosp = 0.9, and mutation probability mutp = 0.1 in MOGA and NSGA-II; maximum iteration max_Iter = 300, particle size nP = 150, the size of external repository nR = 40, the count of grids per dimension nG = 5, initial particle inertia weight w = 1, and initial particle personal learning coefficient c1 = 1, and the initial particle global learning coefficient c2 = 2 in EMOPSO. The calculation results provided by the three algorithms are presented in Table 10.Table 10
Optimization result comparison by each algorithm.
InstanceEMOPSONSGA-IIMOGATC ($)NVT (min)TC ($)NVT (min)TC ($)NVT (min)1157661617648241656723218846172182727199772932099822216110322228113042357131825021421267315225265815192777162327151721630391518326117313200163372854132328321334290114308334416263469174034251639940011735434618374169173910436719264601192845382130111695627196172618748241219907222167734200593313229482125419222479102514248911182629122726531426152896141331311429310116281631191714321619263354202417284915192969162929981735183468162134851727357419241941061621445818344210193620465819284919203847962139Average288713213069142930271530t-test9.3E − 087.9E − 072.0E − 071.4E − 081.8E − 084.8E − 07In Table10, the calculation of t-test shows a substantial difference among the optimal results calculated by the algorithms and proves the reasonability of the comparison among the three algorithms. By comparing the average values of TC, NV, and WT, the results provided by EMOPSO are all superior to those provided by MOGA and NSGA-II. First, the average operating cost in EMOPSO is $2887, which saves by $182 compared with MOGA and $140 compared with NSGA-II. Second, compared with 14 required vehicles in MOGA and 15 required vehicles in NSGA-II, the average required number of vehicles in EMOPSO is 13. Third, EMOPSO obtains the minimum average value of waiting time as 21, compared with 29 in MOGA and 30 in NSGA-II. Therefore, the proposed EMOPSO is more effective than the two other algorithms.


## 5.3. Management Insights

The design of a two-echelon logistics network is crucial for modern supply chain management because the logistics operating cost constitutes the majority of the expenses of companies. In this study, the proposed 2E-LRPTWTRS optimization provides an effective method with respect to tackling the facility and routing decision simultaneously, and the introduction of transportation resource sharing considerably improves the utilization of delivery vehicles. The management insights concluded are as follows:(1)
An effective location strategy consists of determining the reasonable number and geographical location of opened logistics facilities. Compared with the initial logistics network, the location strategy provided by the proposed 2E-LRPTWTRS optimization enables the customers of a distribution area to be served by fewer DCs with a higher distribution efficiency. In addition, the vehicle-routing optimization based on the determined location strategy can design a more well-organized distribution routes for the two-echelon logistics network, which greatly reduces a series of unreasonable distribution phenomenon in the initial logistics network. Therefore, the proposed 2E-LRPTWTRS optimization can provide logistics enterprises an effective reference for tackling the location decisions and vehicle-routing scheme simultaneously to guarantee the efficient operations of the two-echelon logistics network.(2)
The introduction of transportation resource-sharing strategy in the proposed 2E-LRPTWTRS optimization plays an important role in reducing logistics cost for a two-echelon logistics network. Different distribution routes have varied starting and ending times because of the time window constraints. Different from the phenomenon that each distribution route requires one delivery vehicle in the traditional distribution, a delivery vehicle can be shared within and among different logistics facilities and perform multiple distribution routes only if the time window allows, which can remarkably reduce the number of required vehicles in the two-echelon logistics network. This strategy not only satisfies the requirement of logistics enterprises in reducing operating cost but also responds to the call of our government to promote sustainable development, which is demonstrated as a promising strategy to be encouraged.


## 6. Conclusions

In this study, vehicle-sharing and time window constraints are introduced to optimize the 2E-LRPTWTRS. A 2E-LRPTWTRS optimization model is established to address the facility location problem and VRP with time windows and resource sharing to reduce the total operating cost as well as required number of delivery vehicles simultaneously. The contributions of this study mainly consist of two aspects. First, establishing a biobjective optimization model seeks to achieve the minimization of the total cost and the number of required vehicles. Second, a two-stage hybrid algorithm, which is composed ofk-means clustering in the first stage and EMOPSO algorithm in the second stage, is designed. Then, 20 small-scale instances are used to conduct an algorithm comparison. The superiority of EMOPSO is demonstrated by comparing the optimization results with MOGA and NSGA-II.The practical importance of the abovementioned model and algorithm are further proven through a real case of 2E-LRPTWTRS in Chongqing. First, the computation results ofk-means clustering provide three candidate location strategies, and the final optimal location strategy is determined through optimization result comparison. By comparing the optimization results, the location strategy of {DC3, DC4, DC5, DC6} is determined as the best strategy to proceed with the subsequent vehicle routing optimization. Second, the vehicle-routing optimization with transportation resource sharing substantially eliminates the possible unreasonable transportation phenomenon and greatly saves the transportation resource. Compared with the initial logistics network, the number of vehicles decreases from 31 to 5 and the total operating cost decreases from $5730 to $3581. Therefore, the required vehicles and total operating cost save by 26 and $2149 separately through the proposed model and algorithm.The analysis results of transportation resource sharing demonstrate the contribution of such strategy to the optimization results as well as the rationality of adapting vehicle sharing within and among different logistics facilities. In selecting the scheme of vehicle sharing, three scenarios are considered, and a result optimization among these three scenarios is conducted to determine the optimal sharing scheme for the proposed 2E-LRPTWTRS optimization. The comparison results show the scheme of sharing vehicles within and among different logistics facilities simultaneously is better than the other scenarios because of the minimum total operating cost and required vehicles. Therefore, transportation resource-sharing strategy can contribute to addressing the proposed 2E-LRPTWTRS optimization with good performance.This study introduces effective, applicable methodologies to solve 2E-LRPTWTRS and presents several related practical references for logistics managers. The limitations of the present study and the directions of future work are as follows: (1) considering multiechelon logistics network in 2E-LRPTWTRS optimization is an interesting direction, (2) considering determining multiple service periods in the two-echelon logistics network is worth studying, (3) the integration of exact algorithms and hybrid heuristics approaches can be considered to improve the accuracy and performance of optimal solutions, and (4) dynamic customer demands should be considered as another research direction in two-echelon location-routing problem.



---

*Source: 8280686-2021-05-10.xml*