# Two-Echelon Multidepot Logistics Network Design with Resource Sharing

Authors: Siyu Luo; Yong Wang; Jinjun Tang; Xiangyang Guan; Maozeng Xu
Journal: Journal of Advanced Transportation
 (2021)

---

## Abstract

Resource sharing within a logistics network offers an effective way to solve problems resulting from inefficient and costly operations of individual logistics facilities. However, the existing analysis of resource sharing and profit allocation is still limited. Therefore, this study aims to model resource sharing in two-echelon delivery and pickup logistics networks to improve the overall efficiency and decrease the total network operating cost. A bi-objective integer programming model is first proposed for two-echelon collaborative multidepot pickup and delivery problems with time windows (2E-CMDPDTW) to seek the minimization of operating costs and number of vehicles. Integrating a customer clustering algorithm, a greedy algorithm, and an improved nondominated sorting genetic algorithm-II (Im-NSGA-II), a hybrid method is then designed to handle the 2E-CMDPDTW model. The customer clustering and the greedy algorithms are employed to generate locally optimized initial solutions to accelerate the calculating velocity and guarantee the diversity of feasible solutions. The Im-NSGA-II combines the order crossover operation and the polynomial mutation process to find the optimal solution of the 2E-CMDPDTW. The comparative results show that the proposed hybrid method outperforms the NSGA-II and the multiobjective genetic algorithm. Furthermore, a Shapley value method is used for allocating total profits of established alliances and finding an optimal coalition sequence of the logistics facilities joining alliances based on the strictly monotonic path strategy. Finally, a case study of 2E-CMDPDTW in Chongqing China is conducted to validate the feasibility. Results indicate that this study contributes to long-term partnerships between logistics facilities within multi-echelon logistics networks in practice and contributes to the long-term sustainability of urban logistics pickup and delivery networks’ development.


---

## Body

## 1. Introduction

The multidepot pickup and delivery problem (MDPDP) is a generalization of the multidepot vehicle routing problem (MDVRP) aiming to satisfy customer demands through optimizing vehicle routes among logistics facilities, delivery customers, and pickup customers in a multidepot logistics network. To achieve efficiency improvement under practical constraints, the collaborative multidepot pickup and delivery problem with time windows (CMDPDTW) has been proposed [1, 2]. Demands for more efficient logistics, for example, the urban express delivery, are expanding as the logistics sector rapidly develops and the complexity of logistics networks increases. In 2019, the total expenditure of logistics in China reached 298 trillion RMB, and the express volume of commodities rose by 24% to 63 billion RMB from the year earlier [3]. It is important to improve the performance of the logistics networks simultaneously considering pickup and delivery demands via transportation resource optimization and cost reduction. Logistics pickup and delivery networks (LDPN) are facing a shift from simple structures (e.g., one echelon) to complex ones (e.g., multiple echelons) because of extensive logistics demands and diverse customer requirements. In addition, the optimization of two-echelon logistics pickup and delivery networks (2E-LPDN) and the design of collaborative mechanisms are conducive to build intelligent transportation systems and design sustainable logistics networks.In traditional two-echelon multidepot logistics pickup and delivery networks (2E-MDLPDN), each logistics facility operates independently. In other words, each logistics facility always serves its customers regardless of where a customer is located. The independent operation mode causes operating problems, such as long-haul and cross transportations, long delivery periods, and transportation resource underutilization in the traditional logistics networks [4]. In addition to operating problems within the network, service functionality and capability of logistics facilities are also critical. For a logistics facility, cross-regional demands of customers may cause service unavailable or high service costs, and insufficient or redundant service capabilities are common and result in customer loss or resource underutilization in practice. Thus, reasonable resource configuration and vehicle scheduling among multiple logistics facilities are crucial for designing an optimized logistics network structure.For the underused resources and service capabilities, the mode of a sharing economy provides individuals and groups potential ways to save costs and obtain profits, which is promising and validated in transportation-related issues such as car-sharing services such as Uber and Zipcar [5, 6]. When designing a logistics transportation network, designers usually face the problem of how to simultaneously maximize the profits and minimize the operating cost of the whole network when each logistics facility has idle vehicles and excessive service requirements owing to independent operation. The idea of sharing economy is referred to as a considerable candidate to solve the mentioned problem and further optimize the designed logistics network. Therefore, both the transportation routing optimization and the global resource integration in the 2E-MDLPDN are required and necessary.In the 2E-MDLPDN, products are gathered and then transported among logistics facilities, such as pickup centers and distribution centers, in the first echelon, and pickup and delivery activities between the appropriate logistics facilities and their customers happen in the second echelon. This configuration of the network provides a collaborative basis for logistics facilities in the network, and it is critical and effective to establish collaborative alliances for facilitating resource sharing [7]. In terms of sharing customer information and transportation resources, the scattered customers originally belonging to different logistics facilities can be reallocated to appropriate ones to reduce costs and meet the constraints (e.g., time windows for service). Despite the advantages of collaboration, a critical issue on allocating the cost savings (i.e., profits produced via alliances) is its contributions to the stability of an alliance [8, 9]. With cost savings and reasonable profit allocations, alliances can be more stable and attractive to potential logistics facilities trapped in a traditional network. Therefore, this study aims to provide an effective approach to promote sharing of logistics resources and profit allocation in the two-echelon collaborative multidepot pickup and delivery problem with time windows (2E-CMDPDTW) to optimize the logistics network and minimize logistics costs.In this study, a 2E-CMDPDTW is proposed and solved considering the collaboration among multiple pickup centers and distribution centers within the logistics network. The transportation resources and customer information are shared in the 2E-CMDPDTW based on establishing collaborative alliances, and then, the services of different logistics facilities can be integrated to reduce operating costs and increase revenues of logistics facilities. A bi-objective integer programming model is first proposed to minimize the total operating cost of the whole logistics network and the number of vehicles used by all logistics facilities. Then, a hybrid method is developed to solve the 2E-CMDPDTW optimization, which includes a customer clustering algorithm and a greedy algorithm to generate the initial solution and the improved nondominated sorting genetic algorithm-II (Im-NSGA-II) to find the optimal solution. Finally, a profit allocation method, namely, the Shapley value model, is employed to achieve an efficient allocation of the total profit to the logistics facilities participating in the collaboration, and the strictly monotonic path (SMP) strategy is used for determining the optimal joining sequence that guarantees the alliance stability. Both the developed bi-objective optimization model and the designed intelligent algorithm provide insights into solving two-echelon multidepot logistics network optimization problems with resource sharing in the urban logistics transportation. Besides, sharing resources within the urban logistics transportation system is demonstrated and validated using real-world case data. Consequently, this study is conducive to improving the operation and management of intelligent transportation systems as well as the development of smart cities.The remainder of this study is organized as follows. The related studies are reviewed and summarized in Section2. The problem statement and mathematical model of the 2E-CMDPDTW optimization are presented in Section 3. A three-stage method consisting of customer reassignment and initial solution generation, two-echelon route optimization, and profit allocation and alliance stability analysis is discussed in Section 4. A case study of a 2E-LPDN in Chongqing is conducted to illustrate and demonstrate the proposed model and the hybrid method in Section 5, followed by conclusions and future research directions.


## 2. Literature Review

Traditional multidepot vehicle routing problems with time windows (MDVRPTW) attempt to seek the optimal routes to minimize traveled distances of the vehicles and the total cost and satisfy customers’ demands. Due to the complexity of the network structure and the diversity of customers’ demands, 2E-CMDPDTW referring to an extension of the MDVRPTW concentrates on the collaboration between logistics facilities [7]. Therefore, the literature pertaining to the MDVRPTW, collaboration and resource sharing in logistics networks, solution methods for the 2E-CMDPDTW, and profit allocation methods in collaborative networks are reviewed.


### 2.1. MDVRPTW

In contrast to the traditional VRP, the MDVRP enables vehicles to start from multiple depots and return to the original depot or other depots after providing corresponding pickup or delivery services. In practice, the time windows of pickup and delivery processes can greatly affect customers’ satisfaction degrees with logistics service quality. Thus, related constraints for time windows are commonly introduced in the MDVRP, in which a new problem is produced and named MDVRPTW [10]. The MDVRPTW and its derivatives have received much attention both in the theoretical and practical fields. For example, Rabbouch et al. [11] reviewed the relevant literature concerning the MDVRPTW and proved that the satisfaction of time windows is important and necessary.The MDVRPTW is the same as the MDVRP referring to an NP-hard optimization problem and mainly solved with heuristic algorithms [12]. Many researchers, such as Surekha and Sumathi [13] and the review of Karakatič and Podgorelec [14], employed genetic algorithms or related genetic approaches to solve the MDVRPTW. In addition, Cordeau et al. [15] and Schneider [16] studied a tabu algorithm and an improved tabu algorithm. An ant colony optimization (ACO) is also a prominent method to solving these problems [17–19]. To reduce transportation costs and environmental impacts, a collaborative scheme is integrated into MDVRPTW optimizations [20]. For the collaboration in the multicenter vehicle routing problems (MCVRP), Xu et al. [21] established a mixed-integer linear programming model to study a collaborative multidepot petrol station replenishment problem with multicompartments and the assignment of the customers’ time windows. Considering the horizontal collaboration between distribution centers, Muñoz-Villamizar et al. [22] optimized the urban logistics network to minimize the overall cost and reduce environmental impacts. The existing research of the MDVRPTW provides a research framework for the proposed 2E-CMDPDTW in this study. Moreover, research on collaboration provides a research basis for resource sharing in the MDPDP.


### 2.2. Collaboration and Logistics Resource Sharing

The collaboration and resource sharing of logistics facilities in a logistics network attracts researchers and practitioners as they contribute to cost savings and efficiency improvement [23, 24]. In recent years, many scholars have studied the collaboration and resource sharing of logistics facilities in the MDVRP with simultaneous pickup and delivery problems (PDP) [25, 26]. Quintero-Araujo et al. [27] discussed the horizontal collaboration of logistics facilities in logistics networks and proved that the efficiency of fully collaborative logistics networks is higher than that of semicollaborative or noncollaborative logistics networks. Besides the economic benefits, such as transportation cost reduction, both the society and the environment are beneficial through alleviating traffic congestions and reducing greenhouse gas emissions. To ensure the stability of collaborative logistics networks, the negative impacts of the uncertainties on resource matching are required to be reduced [28]. With the increasing complexity of logistics networks, resource sharing has been modeled into the optimization model of the multi-echelon logistics delivery and pickup network design and can further reduce operating costs of logistics networks [1, 27, 28]. Considering the resource sharing of multiple distribution centers in a two-echelon distribution network, Wang et al. [23] proposed a bi-objective mixed-integer programming model to minimize the operation cost and carbon dioxide emissions of the whole logistics network. Lin [29] studied the coordination and sharing of transportation resources in a PDP with time windows, in which the transportation resources are divided into van trains and foot couriers and proposed a coordination strategy to optimize the service time and cost. For the last-mile PDP, Bhasker et al. [30] presented an optimization model to reduce the delivery and pickup cost of the collaborative logistics network and analyzed the sharing of the extra profits generated by the collaboration of logistics enterprises. Therefore, the resource sharing through the collaboration of multiple pickup centers and distribution centers is integrated to optimize the 2E-MDLPDN.


### 2.3. Relevant Solution Methods for the 2E-CMDPDTW

The previous efforts on solving VRP and their variants provide a solid foundation for the 2E-CMDPDTW, and they can be categorized according to the logistics tasks, configurations of logistics networks, and operating modes. For example, when a task is conducted by a single DC, vehicle routing problems with simultaneous pickups and deliveries (VRPSPD) are usually developed without time window violations [31, 32]; otherwise, simultaneous delivery and pickup problems with time windows (SDPPTW) are formulated [33, 34]. In addition, a multiproduct vehicle routing problem with time windows (MPVRPTW) [35] is formulated when delivery tasks with several products are conducted by a single DC. A basic two-echelon multidepot vehicle routing problem (2E-MDVRP) [36] can be described as an extension of a two-echelon vehicle routing problem with time windows (2E-VRPTW) [37]. Both of them focus on optimizing the vehicle routes in the first- and second-echelon logistics networks to achieve certain objectives (e.g., minimizing logistics operating costs and maximizing customer satisfaction).Due to the diversities of fleets, pickup, and delivery demands, number of depots, and region division, various problems are further developed, such as multivehicle multidepot pickup and delivery problems with time windows (M-MDPDPTW) [38], multidepot heterogeneous-fleet pickup, delivery problems with soft time windows (MDHPDSTW) [39], multidepot vehicle routing problems with multiple pickup and delivery demands (MDVRPMPDD) [40], and two-region multidepot pickup and delivery problems (2R-MDPDP) [41]. Moreover, two-echelon vehicle routing problems with pickup and delivery (2E-VRPPD) [42] and two-echelon vehicle routing problems with simultaneous pickup and delivery (2E-VRPSPD) [43] are formulated based on the number of depots (e.g., single or multiple) and the existence of time window violations. For operating modes, cooperation usually occurs among DCs or multidepots, and thus, horizontal collaboration in multidepot vehicle routing problems (HC-MDVRP) [44] and collaborative multicenter vehicle routing problems (CMCVRP) [45] are further investigated. Similarly, two-echelon logistics pickup and delivery problems (2E-LPDP) [46] are formulated and contribute to performance improvement via cooperation and the satisfaction of customers’ demands. The comparative results of related efforts and this work are presented in Table 1.Table 1
Comparison between the abovementioned research studies and this work.
ReferenceMultiple depotsMultiple echelonsPickup and deliveryTime windowsLogistics collaborationResources sharingMultiple objectivesCustomer clusteringHeuristic algorithmHybrid algorithmProblem studiedSubramanian et al. [31]✓✓✓VRPSPDGoksal et al. [32]✓✓✓VRPSPDWang and Chen [33]✓✓✓SDPPTWMahmoudi and Zhou [34]✓✓SDPPTWGovindan et al. [35]✓✓✓✓✓✓MPVRPTWZhou et al. [36]✓✓✓✓✓2E-MDVRPLi et al. [37]✓✓✓✓2E-VRPTMDridi et al. [38]✓✓✓✓M-MDPDPTWBettinelli et al. [39]✓✓✓MDHPDSTWKachitvichyanukul et al. [40]✓✓✓✓✓MDVRPMPDDSoriano et al. [41]✓✓✓✓✓2R-MDPDPMartins et al. [42]✓✓✓✓✓2E-VRPPDBelgin et al. [43]✓✓✓✓✓2E-VRPSPDPérez-Bernabeu et al. [44]✓✓✓✓✓✓✓HC-MDVRPWang et al. [45]✓✓✓✓✓CMCVRPWang et al. [46]✓✓✓✓✓✓✓2E-LPDPAdelzadeh et al. [47]✓✓✓✓✓✓MDVRPTWThis work✓✓✓✓✓✓✓✓✓✓2E-CMDPDTWNote: ✓ denotes that the content is considered.Research on PDP, MDVRPTW, and their variants established an optimization model with a single objective to optimize the logistics network, especially to minimize the cost, waiting time of customers, and number of transportation tools of the network [31, 37, 45]. However, compared with the optimization model with a single objective, the optimization model with multiobjective can evaluate the efficiency of the logistics network effectively to better optimize the logistics network [47]. For example, Govindan et al. [35] established a bi-objective mixed-integer programming model, and the network cost was considered as the first objective function, and the environmental and social impacts were integrated into the second objective function to optimize the supplier network of multiple commodities. Pérez-Bernabeu et al. [44] developed an iterated local search algorithm to solve the problem of horizontal cooperation between road transportation enterprises and minimize transportation costs and exhaust emissions of vehicles. Besides, the number of vehicles was chosen as one of the objectives. To optimize the MDVRPMPDD, Kachitvichyanukul et al. [40] established a multiobjective metaheuristic model that the objective functions consisted of the total logistics cost, number of vehicles, and number of assigned demands. The PDP is an NP-hard optimization problem, and the 2E-CMDPDTW expands single-echelon logistics networks of the PDP to two-echelon logistics networks. Hence, the 2E-CMDPDTW is also an NP-hard optimization problem and searching for the optimal solution with exact algorithms in large-scale PDPs is cost expensive. Therefore, an efficient heuristic algorithm is required. Goksal et al. [32], Wang et al. [45], and Belgin et al. [43] proposed heuristic solution approaches to study the VRPSPD, CMCVRP, and 2E-VRPSPD, respectively.Despite the powerful computation capacity of heuristic algorithms, more efforts are also made to simplify the design and optimization problems in the two-echelon logistics networks. Clustering algorithms have been combined with heuristic algorithms to solve MDVRP, and the integrated algorithms are capable of effectively reducing computational time [25]. To improve the convergence speed, Crainic et al. [48] developed clustering-based heuristics to solve the 2E-VRPTW. Meihua et al. [49] improved the ACO with multiple neighborhood descent in terms of distance-based clusters. Liu and Liao [50] integrated the k-means clustering with an adaptive large neighborhood search algorithm to simplify the process of solving a collaborative 2E-VRPTW. Moreover, clustering algorithms also play significant roles in improving the performance of heuristic algorithms and simplifying the computing process for large-scale optimization problems. Therefore, in this study, the solution efficiency of the proposed 2E-CMDPDTW can be improved by adopting the clustering algorithm.In recent years, multiobjective heuristic algorithms have become widespread since multiple optimization objectives are inherently conflicting in logistics networks. Paul and Shill [51] proposed an NSGA-II to optimize the two validity indices of the clustering processes, and the results reflected that the proposed NSGA-II has good performance in both gene representation and nongene representation datasets. Vachhani et al. [52] proposed a new agglomerative hierarchical clustering process to accelerate the convergence process of the NSGA-II, and the experiments on the standard test of the multiobjective evolutionary algorithm showed that the improved NSGA-II is better than NSGA-II in terms of the existing diversity method. Due to the high efficiency in solving multiobjective optimization problems, the NSGA-II has been widely used and improved to solve the model for the design and optimization of the logistics network. Bandyopadhyay and Bhattacharya [53] proposed an improved NSGA-II to solve the triobjective problem including minimizing the total delivery cost, the variance of order quantity, and the total inventory cost of the two-echelon logistics network. Li et al. [54] improved a hybrid evolutionary algorithm integrating a local search and the NSGA-II to optimize a reverse logistics network with the minimum total cost and maximum benefits. Wang et al. [46] presented an improved NSGA-II based on the clustering algorithm to solve the MDVRP in a 2E-LPDN.


### 2.4. Profit Allocation

Both the operating efficiency and cost of logistics systems can be improved via resource sharing in collaborative networks. Reasonable profit allocation methods can further facilitate the stable development of collaborative alliances. Considering the fairness preferences of different supply chain members, Wen and Li [55] developed two profit allocation models to analyze the impact of profit allocation on the stability of the supply chain. Focusing on asymmetric information, Liu et al. [56] also presented a profit allocation model for the alliance consisting of different logistics enterprises. Moreover, the game theory is also prominent to solve profit allocation problems. Wang et al. [57] used the game quadratic programming to allocate profits in a cooperative distribution network with time window assignments. In reverse logistics, Liu and Liao [50] proposed an improved Shapley value model to promote the distribution of total profit, and the calculation of transportation costs and carbon emissions is the basis of the model.Despite the abovementioned contributions to 2E-CMDPDTW optimization, its following aspects still need attention. (1) The collaboration among multiple distribution centers or pickup centers in a delivery or pickup network has been discussed, but the resource sharing between two-echelon logistics facilities, especially in the 2E-MDLPDN, is still limited. (2) Multiobjective optimization problems are overlooked. The economic factors such as network operating cost are usually modeled in the 2E-MDLPDN optimization problems. Transportation resources, for instance, the number of vehicles, should be taken into account for the 2E-MDLPDN optimization. (3) Difficulties (e.g., time-consuming convergence and low-quality solution) often occur when using the conventional heuristic algorithms to solve a large-scale 2E-MDLPDN. On one hand, existing customer clustering algorithms developed to improve computational efficiency are not appropriate for the 2E-CMDPDTW proposed in this study. On the other hand, convergence speed acceleration via the integration of heuristic algorithms to improve the quality of the initial feasible solutions needs to be further explored. (4) More efforts should be made on a collaborative 2E-MDLPDN instead of the current focus on horizontal collaboration and profit allocation in a single-echelon logistics network. The profits of a collaborative alliance should be reasonably distributed in a multi-echelon logistics network to stabilize long-term partnerships. In addition, the evaluation of the stability of established alliances and the analysis of the optimal order of each member joining the alliance are necessary.In response to the research needs, the major contributions of this study can be summarized in the following aspects. (1) Constructing a collaborative 2E-MDLPDN, and discussing the influence of resource sharing on the network optimization. (2) Establishing a bi-objective integer programming model for minimizing the total operating cost and the number of vehicles for the 2E-CMDPDTW. (3) Designing a hybrid method integrating the customer clustering process, the greedy algorithm, and the Im-NSGA-II to solve the optimization model efficiently. (4) Proposing an optimal alliance sequence selection method on the basis of the Shapley value method and the SMP method to solve the problems of profit allocation and the order of joining an alliance.


## 2.1. MDVRPTW

In contrast to the traditional VRP, the MDVRP enables vehicles to start from multiple depots and return to the original depot or other depots after providing corresponding pickup or delivery services. In practice, the time windows of pickup and delivery processes can greatly affect customers’ satisfaction degrees with logistics service quality. Thus, related constraints for time windows are commonly introduced in the MDVRP, in which a new problem is produced and named MDVRPTW [10]. The MDVRPTW and its derivatives have received much attention both in the theoretical and practical fields. For example, Rabbouch et al. [11] reviewed the relevant literature concerning the MDVRPTW and proved that the satisfaction of time windows is important and necessary.The MDVRPTW is the same as the MDVRP referring to an NP-hard optimization problem and mainly solved with heuristic algorithms [12]. Many researchers, such as Surekha and Sumathi [13] and the review of Karakatič and Podgorelec [14], employed genetic algorithms or related genetic approaches to solve the MDVRPTW. In addition, Cordeau et al. [15] and Schneider [16] studied a tabu algorithm and an improved tabu algorithm. An ant colony optimization (ACO) is also a prominent method to solving these problems [17–19]. To reduce transportation costs and environmental impacts, a collaborative scheme is integrated into MDVRPTW optimizations [20]. For the collaboration in the multicenter vehicle routing problems (MCVRP), Xu et al. [21] established a mixed-integer linear programming model to study a collaborative multidepot petrol station replenishment problem with multicompartments and the assignment of the customers’ time windows. Considering the horizontal collaboration between distribution centers, Muñoz-Villamizar et al. [22] optimized the urban logistics network to minimize the overall cost and reduce environmental impacts. The existing research of the MDVRPTW provides a research framework for the proposed 2E-CMDPDTW in this study. Moreover, research on collaboration provides a research basis for resource sharing in the MDPDP.


## 2.2. Collaboration and Logistics Resource Sharing

The collaboration and resource sharing of logistics facilities in a logistics network attracts researchers and practitioners as they contribute to cost savings and efficiency improvement [23, 24]. In recent years, many scholars have studied the collaboration and resource sharing of logistics facilities in the MDVRP with simultaneous pickup and delivery problems (PDP) [25, 26]. Quintero-Araujo et al. [27] discussed the horizontal collaboration of logistics facilities in logistics networks and proved that the efficiency of fully collaborative logistics networks is higher than that of semicollaborative or noncollaborative logistics networks. Besides the economic benefits, such as transportation cost reduction, both the society and the environment are beneficial through alleviating traffic congestions and reducing greenhouse gas emissions. To ensure the stability of collaborative logistics networks, the negative impacts of the uncertainties on resource matching are required to be reduced [28]. With the increasing complexity of logistics networks, resource sharing has been modeled into the optimization model of the multi-echelon logistics delivery and pickup network design and can further reduce operating costs of logistics networks [1, 27, 28]. Considering the resource sharing of multiple distribution centers in a two-echelon distribution network, Wang et al. [23] proposed a bi-objective mixed-integer programming model to minimize the operation cost and carbon dioxide emissions of the whole logistics network. Lin [29] studied the coordination and sharing of transportation resources in a PDP with time windows, in which the transportation resources are divided into van trains and foot couriers and proposed a coordination strategy to optimize the service time and cost. For the last-mile PDP, Bhasker et al. [30] presented an optimization model to reduce the delivery and pickup cost of the collaborative logistics network and analyzed the sharing of the extra profits generated by the collaboration of logistics enterprises. Therefore, the resource sharing through the collaboration of multiple pickup centers and distribution centers is integrated to optimize the 2E-MDLPDN.


## 2.3. Relevant Solution Methods for the 2E-CMDPDTW

The previous efforts on solving VRP and their variants provide a solid foundation for the 2E-CMDPDTW, and they can be categorized according to the logistics tasks, configurations of logistics networks, and operating modes. For example, when a task is conducted by a single DC, vehicle routing problems with simultaneous pickups and deliveries (VRPSPD) are usually developed without time window violations [31, 32]; otherwise, simultaneous delivery and pickup problems with time windows (SDPPTW) are formulated [33, 34]. In addition, a multiproduct vehicle routing problem with time windows (MPVRPTW) [35] is formulated when delivery tasks with several products are conducted by a single DC. A basic two-echelon multidepot vehicle routing problem (2E-MDVRP) [36] can be described as an extension of a two-echelon vehicle routing problem with time windows (2E-VRPTW) [37]. Both of them focus on optimizing the vehicle routes in the first- and second-echelon logistics networks to achieve certain objectives (e.g., minimizing logistics operating costs and maximizing customer satisfaction).Due to the diversities of fleets, pickup, and delivery demands, number of depots, and region division, various problems are further developed, such as multivehicle multidepot pickup and delivery problems with time windows (M-MDPDPTW) [38], multidepot heterogeneous-fleet pickup, delivery problems with soft time windows (MDHPDSTW) [39], multidepot vehicle routing problems with multiple pickup and delivery demands (MDVRPMPDD) [40], and two-region multidepot pickup and delivery problems (2R-MDPDP) [41]. Moreover, two-echelon vehicle routing problems with pickup and delivery (2E-VRPPD) [42] and two-echelon vehicle routing problems with simultaneous pickup and delivery (2E-VRPSPD) [43] are formulated based on the number of depots (e.g., single or multiple) and the existence of time window violations. For operating modes, cooperation usually occurs among DCs or multidepots, and thus, horizontal collaboration in multidepot vehicle routing problems (HC-MDVRP) [44] and collaborative multicenter vehicle routing problems (CMCVRP) [45] are further investigated. Similarly, two-echelon logistics pickup and delivery problems (2E-LPDP) [46] are formulated and contribute to performance improvement via cooperation and the satisfaction of customers’ demands. The comparative results of related efforts and this work are presented in Table 1.Table 1
Comparison between the abovementioned research studies and this work.
ReferenceMultiple depotsMultiple echelonsPickup and deliveryTime windowsLogistics collaborationResources sharingMultiple objectivesCustomer clusteringHeuristic algorithmHybrid algorithmProblem studiedSubramanian et al. [31]✓✓✓VRPSPDGoksal et al. [32]✓✓✓VRPSPDWang and Chen [33]✓✓✓SDPPTWMahmoudi and Zhou [34]✓✓SDPPTWGovindan et al. [35]✓✓✓✓✓✓MPVRPTWZhou et al. [36]✓✓✓✓✓2E-MDVRPLi et al. [37]✓✓✓✓2E-VRPTMDridi et al. [38]✓✓✓✓M-MDPDPTWBettinelli et al. [39]✓✓✓MDHPDSTWKachitvichyanukul et al. [40]✓✓✓✓✓MDVRPMPDDSoriano et al. [41]✓✓✓✓✓2R-MDPDPMartins et al. [42]✓✓✓✓✓2E-VRPPDBelgin et al. [43]✓✓✓✓✓2E-VRPSPDPérez-Bernabeu et al. [44]✓✓✓✓✓✓✓HC-MDVRPWang et al. [45]✓✓✓✓✓CMCVRPWang et al. [46]✓✓✓✓✓✓✓2E-LPDPAdelzadeh et al. [47]✓✓✓✓✓✓MDVRPTWThis work✓✓✓✓✓✓✓✓✓✓2E-CMDPDTWNote: ✓ denotes that the content is considered.Research on PDP, MDVRPTW, and their variants established an optimization model with a single objective to optimize the logistics network, especially to minimize the cost, waiting time of customers, and number of transportation tools of the network [31, 37, 45]. However, compared with the optimization model with a single objective, the optimization model with multiobjective can evaluate the efficiency of the logistics network effectively to better optimize the logistics network [47]. For example, Govindan et al. [35] established a bi-objective mixed-integer programming model, and the network cost was considered as the first objective function, and the environmental and social impacts were integrated into the second objective function to optimize the supplier network of multiple commodities. Pérez-Bernabeu et al. [44] developed an iterated local search algorithm to solve the problem of horizontal cooperation between road transportation enterprises and minimize transportation costs and exhaust emissions of vehicles. Besides, the number of vehicles was chosen as one of the objectives. To optimize the MDVRPMPDD, Kachitvichyanukul et al. [40] established a multiobjective metaheuristic model that the objective functions consisted of the total logistics cost, number of vehicles, and number of assigned demands. The PDP is an NP-hard optimization problem, and the 2E-CMDPDTW expands single-echelon logistics networks of the PDP to two-echelon logistics networks. Hence, the 2E-CMDPDTW is also an NP-hard optimization problem and searching for the optimal solution with exact algorithms in large-scale PDPs is cost expensive. Therefore, an efficient heuristic algorithm is required. Goksal et al. [32], Wang et al. [45], and Belgin et al. [43] proposed heuristic solution approaches to study the VRPSPD, CMCVRP, and 2E-VRPSPD, respectively.Despite the powerful computation capacity of heuristic algorithms, more efforts are also made to simplify the design and optimization problems in the two-echelon logistics networks. Clustering algorithms have been combined with heuristic algorithms to solve MDVRP, and the integrated algorithms are capable of effectively reducing computational time [25]. To improve the convergence speed, Crainic et al. [48] developed clustering-based heuristics to solve the 2E-VRPTW. Meihua et al. [49] improved the ACO with multiple neighborhood descent in terms of distance-based clusters. Liu and Liao [50] integrated the k-means clustering with an adaptive large neighborhood search algorithm to simplify the process of solving a collaborative 2E-VRPTW. Moreover, clustering algorithms also play significant roles in improving the performance of heuristic algorithms and simplifying the computing process for large-scale optimization problems. Therefore, in this study, the solution efficiency of the proposed 2E-CMDPDTW can be improved by adopting the clustering algorithm.In recent years, multiobjective heuristic algorithms have become widespread since multiple optimization objectives are inherently conflicting in logistics networks. Paul and Shill [51] proposed an NSGA-II to optimize the two validity indices of the clustering processes, and the results reflected that the proposed NSGA-II has good performance in both gene representation and nongene representation datasets. Vachhani et al. [52] proposed a new agglomerative hierarchical clustering process to accelerate the convergence process of the NSGA-II, and the experiments on the standard test of the multiobjective evolutionary algorithm showed that the improved NSGA-II is better than NSGA-II in terms of the existing diversity method. Due to the high efficiency in solving multiobjective optimization problems, the NSGA-II has been widely used and improved to solve the model for the design and optimization of the logistics network. Bandyopadhyay and Bhattacharya [53] proposed an improved NSGA-II to solve the triobjective problem including minimizing the total delivery cost, the variance of order quantity, and the total inventory cost of the two-echelon logistics network. Li et al. [54] improved a hybrid evolutionary algorithm integrating a local search and the NSGA-II to optimize a reverse logistics network with the minimum total cost and maximum benefits. Wang et al. [46] presented an improved NSGA-II based on the clustering algorithm to solve the MDVRP in a 2E-LPDN.


## 2.4. Profit Allocation

Both the operating efficiency and cost of logistics systems can be improved via resource sharing in collaborative networks. Reasonable profit allocation methods can further facilitate the stable development of collaborative alliances. Considering the fairness preferences of different supply chain members, Wen and Li [55] developed two profit allocation models to analyze the impact of profit allocation on the stability of the supply chain. Focusing on asymmetric information, Liu et al. [56] also presented a profit allocation model for the alliance consisting of different logistics enterprises. Moreover, the game theory is also prominent to solve profit allocation problems. Wang et al. [57] used the game quadratic programming to allocate profits in a cooperative distribution network with time window assignments. In reverse logistics, Liu and Liao [50] proposed an improved Shapley value model to promote the distribution of total profit, and the calculation of transportation costs and carbon emissions is the basis of the model.Despite the abovementioned contributions to 2E-CMDPDTW optimization, its following aspects still need attention. (1) The collaboration among multiple distribution centers or pickup centers in a delivery or pickup network has been discussed, but the resource sharing between two-echelon logistics facilities, especially in the 2E-MDLPDN, is still limited. (2) Multiobjective optimization problems are overlooked. The economic factors such as network operating cost are usually modeled in the 2E-MDLPDN optimization problems. Transportation resources, for instance, the number of vehicles, should be taken into account for the 2E-MDLPDN optimization. (3) Difficulties (e.g., time-consuming convergence and low-quality solution) often occur when using the conventional heuristic algorithms to solve a large-scale 2E-MDLPDN. On one hand, existing customer clustering algorithms developed to improve computational efficiency are not appropriate for the 2E-CMDPDTW proposed in this study. On the other hand, convergence speed acceleration via the integration of heuristic algorithms to improve the quality of the initial feasible solutions needs to be further explored. (4) More efforts should be made on a collaborative 2E-MDLPDN instead of the current focus on horizontal collaboration and profit allocation in a single-echelon logistics network. The profits of a collaborative alliance should be reasonably distributed in a multi-echelon logistics network to stabilize long-term partnerships. In addition, the evaluation of the stability of established alliances and the analysis of the optimal order of each member joining the alliance are necessary.In response to the research needs, the major contributions of this study can be summarized in the following aspects. (1) Constructing a collaborative 2E-MDLPDN, and discussing the influence of resource sharing on the network optimization. (2) Establishing a bi-objective integer programming model for minimizing the total operating cost and the number of vehicles for the 2E-CMDPDTW. (3) Designing a hybrid method integrating the customer clustering process, the greedy algorithm, and the Im-NSGA-II to solve the optimization model efficiently. (4) Proposing an optimal alliance sequence selection method on the basis of the Shapley value method and the SMP method to solve the problems of profit allocation and the order of joining an alliance.


## 3. Problem Statement and Mathematical Model for 2E-CMDPDTW

### 3.1. Problem Statement

Collaboration between logistics facilities in the network contributes to resource integration and optimization so as to minimize the overall operating cost and maximize the revenue of individual logistics facilities. In this study, the 2E-CMDPDTW consists of different logistics facilities containing a logistics center (LC), multiple pickup centers (PCs), and multiple distribution centers (DCs) as the first echelon and a certain number of delivery and pickup customers (Cs) as the second echelon in the 2E-LPDN. Of these facilities, the LC is responsible for collecting and distributing products; and, PCs and DCs provide pickup and delivery service for their customers, respectively. As for the transportation resources, transportation tasks between different logistics facilities in the first-echelon are carried out by semitrailer trucks. For the second-echelon transportation service between a logistics facility and its customers, vehicles are utilized to provide door-to-door services. Fixed time windows are requested by all facilities and customers in the network. If customers were not served within their time windows, the network operating cost would increase along with the decrease in service efficiency.The structures of the 2E-CMDPDTW optimization before the collaboration and after the collaboration are illustrated in Figure1. Without collaboration, logistics facilities operate independently. In comparison, collaborating logistics facilities in the alliance can share customers and transportation resources. In Figure 1, the numbers attached to the links represent the travel time of semitrailer trucks or vehicles on the routes.Figure 1
The 2E-LPDN before and after collaboration. (a) Two-echelon noncollaborative pickup and delivery network. (b) Two-echelon collaborative pickup and delivery network.
(a)(b)In the noncollaborative 2E-LPDN before optimization (Figure1(a)), the LC serves DCs, PCs, and its customers with various demands simultaneously; however, DCs and PCs only serve their customers with delivery and pickup demands, respectively. Three problems arise consequently. (1) There are unreasonable transportation routes (e.g., long-haul and cross transportation). For example, C2 is closer to the LC but served by DC1. This situation can be observed in practice in terms of customer loyalty and long-term contractual partnership. (2) Customers’ time windows are frequently violated. The vehicle, for instance, starting from the LC must wait for open time windows when arriving at C13 and C14. The demand of C12 cannot be satisfied on time since the arrival of the vehicle from DC2 is beyond the given time window. (3) Transportation resources are underutilized. For example, the demand for DC1 served by the LC can be satisfied with a vehicle rather than a semitrailer truck. In other words, other facilities can also be served with the truck. In response to the aforementioned problems, collaboration with resource sharing is effective for efficiency improvement and reducing resource underutilization.Figure1(b) presents the collaborative 2E-LPDN after optimization. The LC plays an important role in organizing delivery tasks and pickup tasks and coordinating logistics facilities within the network. The advantages of collaboration are concluded as follows. First of all, with customer information sharing, the customers in the alliance can be reassigned to new service providers closer to them. In addition, the transportation routes both in the first and second echelon can be optimized to guarantee the given service time windows. Finally, sharing transportation resources, such as vehicles and semitrailer trucks, can effectively avoid idle transportation resources in the collaborative logistics network. For example, the customers re-arranged in the route LC ⟶ C21 ⟶ C22 ⟶ C23 ⟶ C24 ⟶ LC, and the other route PC2 ⟶ C32 ⟶ C33 ⟶ PC2 can be served orderly by the two vehicles on a single transportation task. As a consequence, the aim of reducing the total cost and the number of vehicles is achieved in the optimized 2E-CMDPDTW.Table2 shows the comparative results of the costs, number of trucks, and number of vehicles before and after the 2E-CMDPDTW optimization. Assume that the transportation cost per unit time of the semitrailer trucks equals $30, and all of the delivery cost and pickup cost of vehicles equal $15. The unit time penalty cost of time windows violation is $10, and each of DCs and PCs pays for a fixed collaborative cost of $100 when deciding to collaborate.Table 2
Comparison before and after collaboration in the two-echelon pickup and delivery network.
CaseTransportation cost ($)Pickup and delivery cost ($)Penalty cost ($)Collaborative cost ($)Total operating cost ($)Time windows violationNumber of trucksNumber of vehiclesBefore collaboration660132070020505412After collaboration450825050017750111In Table2, the total operating cost of the optimized 2E-CMDPDTW is decreased by $275, and 3 semitrailer trucks and 1 vehicle are saved. In the network before collaboration, the number of violating time windows equals 5, while all the customers’ time windows are satisfied in the collaborative one. The results reflect that the total operating cost, the number of trucks, the number of vehicles, and the number of violating time windows in the 2E-CMDPDTW are reduced. Therefore, resource sharing can reduce operating costs and save transportation resources in the entire network.


### 3.2. Definitions

In this section, the definition, notations, and decision variables used in this study for the 2E-CMDPDTW optimization are introduced. The definitions and the decision variables are listed in Tables3 and 4, respectively.Table 3
Definitions and notations related to the mathematical model.
SetsDefinitionsLSet of LC,L=L0,l∈LISet of DCs,I=i|i=1,2,3,…,hJSet of PCs,J=j|j=1,2,3,…,h′UdSet of all delivery customers in the collaborative 2E-LPDN,Ud=ud|ud=1,2,3,…,uUpSet of all pickup customers in the collaborative 2E-LPDN,Up=up|up=1,2,3,…,u′MSet of all logistics facilities, containing the LC, DCs, and PCs in the collaborative 2E-LPDN,M=m|m=1,2,3,…,h,h+1,…,h+h′,h+h′+1Ud1Set of delivery customers served by the LC,Ud1=ud1|ud1=1,2,3,…,u1Ud2Set of delivery customers served by the DCs,Ud2=ud2|ud2=1,2,3,…,u2Up1Set of pickup customers served by the LC,Up1=up1|up1=1,2,3,…,u1′Up2Set of pickup customers served by the PCs,Up2=up2|up2=1,2,3,…,u2′DSet of the DCs and their delivery customers,D=I∪Ud2=d|d=1,2,3,…,h,h+1,…,h+u2PSet of the PCs and their pickup customers,P=J∪Up2=p|p=1,2,3,…,h′,h′+1,…,h′+u2′LdSet of the LC and its delivery customers,Ld=L∪Ud1=ld|ld=1,2,3,…,u1+1LpSet of the LC and its pickup customers,Lp=L∪Up1=lp|lp=1,2,3,…,u2+1SSet of semitrailer trucks used for product transportation among LC, DCs, and PCs,S=s|s=1,2,3,…,s′VSet of vehicles used for product transportation between logistics facilities and customers,V=v|v=1,2,3,…,v′OsSet of the DCs and the PCs served by semitrailer trucks, s∈SSivdSet of delivery customers whose delivery service are provided by vehiclev from DC i, Sivd⊆Ud2, i∈I and v∈VSjvpSet of pickup customers whose pickup service are provided by vehiclev from PC j, Sjvp⊆Up2, j∈J and v∈VSlvdSet of delivery customers whose delivery service are provided by vehiclev from LC l, Slvd⊆Ud1, l∈L and v∈VSlvpSet of pickup customers whose pickup service are provided by vehiclev from LC l, Slvp⊆Up1, l∈L and v∈VParametersLmm′sDistance of semitrailer trucks starts from logistics facility m to m′, s∈S, mandm′∈M,andm≠m′Ldud2vDistance of vehiclev starts from the DC or delivery customer d to ud2, v∈V, d∈D, and ud2∈Ud2Lpup2vDistance of vehiclev starts from the PC or pickup customer p to up2, v∈V, p∈P, and up2∈Up2Lldud1vDistance of vehiclev starts from the LC or delivery customer ld to ud1, v∈V, ld∈Ld, and ud1∈Ud1Llpup1vDistance of vehiclev starts from the LC or pickup customer lp to up1, v∈V, lp∈Lp, and up1∈Up1fsDiesel fuel consumption rate per mile for semitrailer trucks (gallon/mile), s∈SfvPetrol fuel consumption rate per mile for vehiclev (gallon/mile), v∈VpsUnit price of diesel used in semitrailer trucks (dollar/gallon)pvUnit price of petrol used in vehicles (dollar/gallon)MsAverage annual cost for the maintenance of semitrailer trucksMvAverage annual cost for the maintenance of vehiclesBNumber of working periods per yearQiDemand quantity of DCi for delivery, i∈IQjDemand quantity of PCj for pickup, j∈IRud1Demand quantity of delivery customerud1 served by the LC, ud1∈Ud1Rup1Demand quantity of pickup customerup1 served by the LC, up1∈Up1Rud2Demand quantity of delivery customerud2 served by the DCs, ud2∈Ud2Rup2Demand quantity of pickup customerup2 served by the PCs, up2∈Up2KsApproved load capacity of semitrailer trucks, s∈SKvApproved load capacity of vehiclev, v∈VFmFixed cost of logistics facilitym, m∈MθVariable cost coefficient for products’ transportationcmDiscount on the costs of logistics facilitym, if collaboration occurs (dollar), m∈MOsNumber of DCs and PCs that provides transportation services by semitrailer trucks, s∈SSivdNumber of delivery customers whose delivery service are provided by vehiclev from DC i, i∈I and v∈VSjvpNumber of pickup customers whose pickup service are provided by vehiclev from PC j, j∈J and v∈VSlvdNumber of delivery customers whose delivery service are provided by vehiclev from LC l, l∈L and v∈VSlvpNumber of pickup customers whose pickup service are provided by vehiclev from LC l, l∈L and v∈Vαi,βiTime window of the DCi, i∈Iαj,βjTime window of the PCj, j∈Jαud,βudTime window of delivery customerud, ud∈Udαup,βupTime window of pick customerup, up∈UpTsMaximum travel time allowed for semitrailer trucks on every route, s∈STvMaximum travel time allowed for vehiclev on every route, v∈VweThe unit penalty cost for semitrailer trucks and vehicles arriving earlywdThe unit penalty cost for semitrailer trucks and vehicles arriving delayatisTime of semitrailer trucks arrival at the DC i, s∈S and i∈IatjsTime of semitrailer trucks arrival at the PC j, s∈S and j∈JativTime of vehiclev arrival at DC i, v∈V and i∈IatjvTime of vehiclev arrival at PC j, v∈V and j∈JatudvTime of vehiclev arrival at delivery customer ud, v∈V and ud∈UdatupvTime of vehiclev arrival at pickup customer up, v∈V and up∈Upatud2vTime of vehiclev arrival at delivery customer ud2, v∈V and ud2∈Ud2atup2vTime of vehiclev arrival at pickup customer up2, v∈V and up2∈Up2atud1vTime of vehiclev arrival at delivery customer ud1, v∈V and ud1∈Ud1atup1vTime of vehiclev arrival at pickup customer up1, v∈V and up1∈Up1dtmsTime of semitrailer trucks depart from logistics facility m, m∈M and s∈SdtdvTime of vehiclev depart from DC i or delivery customer ud2, v∈V, i∈I, and ud2∈Ud2dtpvTime of vehiclev depart from PC j or pickup customer up2, v∈V, j∈J, and up2∈Up2dtldvTime of vehiclev depart from LC l or delivery customer ud1, v∈V, ld∈Ld, and ud1∈Ud1dtlpvTime of vehiclev depart from LC l or pickup customer up1, v∈V, lp∈Lp, and up1∈Up1dtud2vTime of vehiclev depart from delivery customer ud2, v∈V and ud2∈Ud2dtup2vTime of vehiclev depart from pickup customer up2, v∈V and up2∈Up2drlsReturn time of semitrailer trucks to LC l, l∈L and s∈SdrlvReturn time of vehiclev to LC l, l∈L and v∈Vtmm′sTime of semitrailer trucks travel from logistics facility m to m′, m∈M and s∈StmisTime of semitrailer trucks travel from logistics facility m to DC i, m∈M, i∈I, and s∈StmjsTime of semitrailer trucks travel from logistics facility m to PC j, m∈M, i∈I, and s∈Stdud2vTime of vehiclev travel from DC i or ud2 to ud2, v∈V, d∈D, i∈I, and ud2∈Ud2tpup2vTime of vehiclev travel from PC j or up2 to up2, v∈V, p∈P, j∈J, and up2∈Up2tldud1vTime of vehiclev travel from LC l or ud1 to ud1, v∈V, ld∈Ld, and ud1∈Ud1tlpup1vTime of vehiclev travel from LC l or up1 to up1, v∈V, lp∈Lp, and up1∈Up1tud2ivTime of vehiclev travel from ud2 to DC i, v∈V, i∈I, and ud2∈Ud2tup2jvTime of vehiclev travel from up2 to PC j, v∈V, j∈J, and up2∈Up2Table 4
Decision variables related to the mathematical model.
VariablesDefinitionsxmm′sxmm′s=1, if semitrailer truck s starts from logistics facility m to m′; otherwise, xmm′s=0, mandm′∈M and s∈Sxmisxmis=1, if semitrailer truck s starts from logistics facility m to DC i; otherwise, xmis=0, m∈M and s∈Sxmjsxmjs=1, if semitrailer truck s starts from logistics facility m to PC j; otherwise, xmjs=0, m∈M and s∈Sxdud2vxdud2v=1, if vehicle v starts from d to ud2; otherwise, xdud2v=0, d∈D, v∈V, and ud2∈Ud2xpup2vxpup2v=1, if vehicle v starts from p to up2; otherwise, xpup2v=0, p∈P, v∈V, and up2∈Up2xldud1vxldud1v=1, if vehicle v starts from ld to ud1; otherwise, xldud1v=0, ld∈Ld, ud1∈Ud1, and v∈Vxlpup1vxlpup1v=1, if vehicle v starts from lp to up1; otherwise, xlpup1v=0, lp∈Lp, up1∈Up1, and v∈Vxmuvxmuv=1, if vehicle v starts from logistics facility m to customer u; otherwise, xmuv=0, m∈M, and u∈Ud∪UpZmZm=1, if logistics facility m agrees to join the alliance to cooperate; otherwise, Zm=0, m∈MqmisQuantity of products transported by semitrailer trucks from logistics facility m to DC i, m∈M, s∈S, and i∈IqmjsQuantity of products transported by semitrailer trucks from logistics facility m to PC j, m∈M, s∈S, and j∈JqmudvQuantity of products delivered by vehiclev from logistics facility m to ud, v∈V, m∈M, and ud∈UdqmupvQuantity of products collected by vehiclev from logistics facility m to up, v∈V, m∈M, and up∈Upqiud2vQuantity of products delivered by vehiclev from DC i to ud2, v∈V, i∈I, and ud2∈Ud2qjup2vQuantity of products collected by vehiclev from PC j to up2, v∈V, j∈J, and up2∈Up2qdud2vQuantity of products delivered by vehiclev from DC i or ud2 to ud2, v∈V, d∈D, i∈I, and ud2∈Ud2qpup2vQuantity of products collected by vehiclev from PC j or up2 to up2, v∈V, p∈P, j∈J, and up2∈Up2qldud1vQuantity of products delivered by vehiclev from LC l or ud1 to ud1, v∈V, ld∈Ld, and ud1∈Ud1qlpup1vQuantity of products collected by vehiclev from LC l or up1 to up1, v∈V, lp∈Lp, and up1∈Up1


### 3.3. Mathematical Model

In this study, a bi-objective mathematical model is developed to optimize the total operating cost of the network and the number of vehicles in the 2E-CMDPDTW. The objective functions in equations (1) and (2) are as follows:(1)MinF1=TC1+TC2+PC1+PC2+TC3,(2)MinF2=∑m∈M∑v∈V∑u∈Ud∪Upxmuv.TC1 expresses the total operating cost of semitrailer trucks, and it includes the total transportation cost and the total maintenance cost in alliances:(3)TC1=∑m,m′∈M,m≠m′∑s∈Sxmm′s×Lmm′s×fs×ps+maxs∈S∑m∈L∑i∈Ixmis,∑m∈L∑j∈Jxmjs×MsB,where ∑m,m′∈M,m≠m′∑s∈Sxmm′s×Lmm′s×fs×ps represents the transportation cost of semitrailer trucks and maxs∈S∑m∈L∑i∈Ixmis,∑m∈L∑j∈Jxmjs×Ms/B depicts the maintenance cost of semitrailer trucks.TC2 expresses the total delivery cost, pickup cost, and maintenance cost of vehicles in the second echelon in a working period:(4)TC2=∑d,ud2∈Dd≠ud2∑v∈Vxdud2v×Ldud2v×fv×pv+∑m∈I∑u∈Ud2∑v∈Vxmuv×MvB+∑p,up2∈Pp≠up2∑v∈Vxpup2v×Lpup2v×fv×p+∑m∈J∑u∈Up2∑v∈Vxmuv×MvB+∑ld,ud1∈Ldld≠ud1∑v∈Vxldud1v×Lldud1v×fv×pv+∑m∈L∑u∈Ud1∑v∈Vxmuv×MvB+∑lp,up1∈Lplp≠up1∑v∈Vxlpup1v×Llpup1v×fv×pv+∑m∈L∑u∈Up1∑v∈Vxmuv×MvB,where ∑d,ud2∈Dd≠ud2∑v∈Vxdud2v×Ldud2v×fv×pv+∑m∈I∑u∈Ud2∑v∈Vxmuv×Mv/B is the delivery cost and maintenance cost of delivery vehicles starting from DCs, ∑p,up2∈Pp≠up2∑v∈Vxpup2v×Lpup2v×fv×p+∑m∈J∑u∈Up2∑v∈Vxmuv×Mv/B represents the pickup cost and maintenance cost of pickup vehicles starting from PCs, ∑ld,ud1∈Ldld≠ud1∑v∈Vxldud1v×Lldud1v×fv×pv+∑m∈L∑u∈Ud1∑v∈Vxmuv×Mv/B depicts the delivery cost and maintenance cost of delivery vehicles from the LC, and ∑lp,up1∈Lplp≠up1∑v∈Vxlpup1v×Llpup1v×fv×pv+∑m∈L∑u∈Up1∑v∈Vxmuv×Mv/B denotes the pickup cost and maintenance cost of pickup vehicles from the LC.PC1 expresses the penalty cost of the semitrailer trucks for early and delayed arrival at the DCs and PCs:(5)PC1=∑m∈M∑i∈I∑s∈Sqmis×we×maxαi−atis,0+∑m∈M∑i∈I∑s∈Sqmis×wd×maxatis−βi,0+∑m∈M∑j∈J∑s∈Sqmjs×we×maxαj−atjs,0+∑m∈M∑j∈J∑s∈Sqmjs×wd×maxatjs−βj,0,where ∑m∈M∑i∈I∑s∈Sqmis×we×maxαi−atis,0 and ∑m∈M∑i∈I∑s∈Sqmis×wd×maxatis−βi,0 indicate the penalty cost of earliness and delay of the semitrailer trucks for DCs and ∑m∈M∑j∈J∑s∈Sqmjs×we×maxαj−atjs,0 and ∑m∈M∑j∈J∑s∈Sqmjs×wd×maxatjs−βj,0 denote the penalty cost of earliness and delay of the semitrailer trucks for PCs.PC2 represents the penalty cost of the vehicles for early and delayed arrival at the delivery customers and pickup customers:(6)PC2=∑m∈M∪Ud∑ud∈Ud∑v∈Vqmudv×we×maxaud−atudv,0+∑m∈M∪Ud∑ud∈Ud∑v∈Vqmudv×wd×maxatudv−βud,0+∑m∈M∪Up∑up∈Up∑v∈Vqmupv×we×maxaup−atupv,0+∑m∈M∪Up∑up∈Up∑v∈Vqmupv×wd×maxatupv−βup,0,where ∑m∈M∪Ud∑ud∈Ud∑v∈Vqmudv×we×maxaud−atudv,0 and ∑m∈M∪Ud∑ud∈Ud∑v∈Vqmudv×wd×maxatudv−βud,0 depict the penalty cost of earliness and delay of the vehicles for delivery customers and ∑m∈M∪Up∑up∈Up∑v∈Vqmupv×we×maxaup−atupv,0 and ∑m∈M∪Up∑up∈Up∑v∈Vqmupv×wd×maxatupv−βup,0 represent the penalty cost of earliness and delay of the vehicles for pickup customers.TC3 expresses the total cost of logistics facilities including the fixed and variable costs and the collaborative discount for establishing alliances between logistics facilities:(7)TC3=∑m∈MFm+∑i∈I∑m∈Mθqmis+∑j∈J∑m∈Mθqmjs−∑m∈MZmcm,where ∑m∈MFm represents the total fixed cost of the LC, DCs, and PCs, ∑i∈I∑m∈Mθqmis and ∑j∈J∑m∈Mθqmjs depict the total variable cost of DCs and PCs, respectively, and ∑m∈MZmcm denotes the collaboration discount for establishing alliances in terms of government subsidies.The constraints of the model are divided into three parts such as constraints for the first echelon, constraints for the second echelon, and constraints for the first and the second echelons.


#### 3.3.1. Constraints for the First Echelon

Constraint (8) ensures that, if m′=L0, a semitrailer truck starts from the LC and returns to the LC after the task is completed. Otherwise, if m′ is a DC or PC, this constraint imposes the balance of semitrailer trucks entering and leaving the logistics facility. Constraints (9) and (10) guarantee that each DC and PC only can be served by one semitrailer truck, respectively. Constraint (11) eliminates subtours of every route in the first echelon. Constraints (12) and (13) express that the demands of each DC and PC can be satisfied, respectively. Constraints (14) and (15) ensure that the maximum load of semitrailer trucks is able to satisfy the total demand of the logistics facilities which are served by the semitrailer trucks. Constraints (16) and (17) imply that semitrailer trucks serve DCs and PCs within their time windows, respectively. Constraint (18) represents that the semitrailer truck arrives at the DC. Constraint (19) defines that the semitrailer truck arrives at the PC. Constraints (20) and (21) denote that the time of each semitrailer truck departure from and return to the logistics facilities meet the time window of the LC. Constraint (22) limits that the total time a semitrailer truck travels on each route must be less than or equal to the allowed time of each route:(8)∑m∈M,m≠m′xmm′s−∑m∈M,m≠m′xm′ms=0,∀m′∈Mand∀s∈S,(9)∑m∈M,m≠i∑s∈Sxmis≤1,∀i∈I,(10)∑m∈M,m≠j∑s∈Sxmjs≤1,∀j∈J,(11)∑m∈Os∑m′∈Os,m′≠mxmm′s=Os−1,∀s∈SandOs≥2,(12)∑s∈S∑m∈Mqmis×xmis≥Qi,∀i∈Iandi≠m,(13)∑s∈S∑m∈Mqmjs×xmjs≥Qi,∀j∈Jandj≠m,(14)∑i∈I∑m∈Mqmis×xmis≤Ks,∀s∈S,(15)∑j∈J∑m∈Mqmjs×xmjs≤Ks,∀s∈S,(16)αi≤atis≤βi,∀i∈Iand∀s∈S,(17)αj≤atjs≤βj,∀j∈Jand∀s∈S,(18)dtms+tmis≤atis,∀m∈M,∀i∈I,and∀s∈S,(19)dtms+tmjs≤atjs,∀m∈M,∀j∈J,and∀s∈S,(20)αl≤dtms≤βl,∀l∈L,∀s∈S,andm=l0,(21)αl≤drls≤βl,∀l∈L,∀s∈S,andm=l0,(22)∑m,m′∈M,m≠m′tmm′s×xmm′s≤Ts,s∈S.


#### 3.3.2. Constraints for the Second Echelon

Constraints (23) and (24) ensure that a delivery vehicle serving DCs’ customers originates from a DC to serve customers and returns to the DC when it completed a task. Constraint (25) imposes the balance of delivery vehicles entering and leaving the delivery customers served by DCs. Constraints (26) and (27) guarantee that a pickup vehicle serving PCs’ customers originates from a PC and returns to the PC after completing a task. Constraint (28) imposes the balance of pickup vehicles entering and leaving the pickup customers served by PCs. Constraints (29) and (30) express that a delivery vehicle serving the LC originates from the LC and returns to the LC. Constraint (31) imposes the balance of delivery vehicles entering and leaving the delivery customers served by the LC. Constraints (32) and (33) imply that a pickup vehicle serving the LC originates from the LC and returns to the LC. Constraint (34) imposes the balance of pickup vehicles entering and leaving the pickup customers served by the LC. Constraints (35)–(38) are used to eliminate subtours of every route in the second echelon. Constraints (39)–(42) are used to guarantee the demands of each delivery or pickup customer can be satisfied. Constraints (43)–(46) stipulate that the maximum load of each vehicle can satisfy the total demand of customers served by this vehicle. Constraints (47) and (48) represent that all vehicles serve delivery and pickup customers within their time windows. Constraints (49)–(52) define that the vehicle arrives at the delivery and pickup customers. Constraints (53)–(56) limit that the total travel time for a vehicle must be less than or equal to the allowed time of each route:(23)∑ud2∈Ud2xdud2v=1,∀d∈Iand∀v∈V,(24)∑ud2∈Ud2xud2dv=1,∀d∈Iand∀v∈V,(25)∑ud2∈Dxdud2v−∑ud2∈Dxud2dv=0,∀d∈D,d≠ud2,and∀v∈V,(26)∑up2∈Up2xpup2v=1,∀p∈Jand∀v∈V,(27)∑up2∈Up2xup2pv=1,∀p∈Jand∀v∈V,(28)∑up2∈Pxpup2v−∑up2∈Pxup2pv=0,∀p∈P,p≠up2,and∀v∈V,(29)∑ud1∈Ud1xldud1v=1,∀ld∈Ldand∀v∈V,(30)∑ud1∈Ud1xud1ldv=1,∀ld∈Ldand∀v∈V,(31)∑ud1∈Ldxldud1v−∑ud1∈Ldxud1ldv=0,∀ld∈Ld,ld≠ud1,and∀v∈V,(32)∑up1∈Up1xlpup1v=1,∀lp∈Lpand∀v∈V,(33)∑up1∈Up1xup1lpv=1,∀lp∈Lpand∀v∈V,(34)∑up1∈Lpxlpup1v−∑up1∈Lpxup1lpv=0,∀lp∈Lp,lp≠up1,and∀v∈V,(35)∑d∈Sivd∑ud2∈Sivdxdud2v=Sivd−1,∀i∈I,d≠ud2,∀v∈V,Sivd⊆Ud2,andSivd≥2,(36)∑p∈Sjvp∑up2∈Sjvpxpup2v=Sjvp−1,∀j∈J,p≠up2,∀v∈V,Sjvp⊆Up2,andSjvp≥2,(37)∑ld∈Slvd∑ud1∈Slvdxldud1v=Slvd−1,∀l∈L,ld≠ud1,∀v∈V,Slvd⊆Ud1,andSlvd≥2,(38)∑lp∈Slvp∑up1∈Slvpxlpup1v=Slvp−1,∀l∈L,lp≠up1,∀v∈V,Slvp⊆Up1,andSlvp≥2,(39)∑d∈D∑v∈Vqdud2v×xdud2v≥Rud2,∀ud2∈Ud2andd≠ud2,(40)∑p∈P∑v∈Vqpup2v×xpup2v≥Rup2,∀up2∈Up2andp≠up2,(41)∑ld∈Ld∑v∈Vqldud1v×xldud1v≥Rud1,∀ud1∈Ud1andld≠ud1,(42)∑lp∈Lp∑v∈Vqlpup1v×xlpup1v≥Rup1,∀up1∈Up1andlp≠up1,(43)∑d∈D∑ud2∈Ud2,d≠ud2qdud2v×xdud2v≤Kv,∀v∈V,(44)∑p∈P∑up2∈Up2,p≠up2qpup2v×xpup2v≤Kv,∀v∈V,(45)∑ud1∈Ud1∑ld∈Ld,ld≠ud1qldud1v×xldud1v≤Kv,∀v∈V,(46)∑up1∈Up1∑lp∈Lp,lp≠up1qlpup1v×xlpup1v≤Kv,∀v∈V,(47)αud≤atudv≤βud,∀ud∈Udand∀v∈V,(48)αup≤atupv≤βup,∀up∈Upand∀v∈V,(49)dtdv+tdud2v≤atud2v,∀ud2∈Ud2,ud2≠d,∀d∈D,and∀v∈V,(50)dtpv+tpup2v≤atup2v,∀up2∈Up2,up2≠p,∀p∈P,and∀v∈V,(51)dtldv+tldud1v≤atud1v,∀ud1∈Ud1,ud1≠ld,∀ld∈Ld,and∀v∈V,(52)dtlpv+tlpup1v≤atup1v,∀up1∈Up1,up1≠lp,∀lp∈Lp,and∀v∈V,(53)∑d,ud2∈D,d≠ud2tdud2v×xdud2v≤Tv,∀v∈V,(54)∑p,up2∈P,p≠up2tpup2v×xpup2v≤Tv,∀v∈V,(55)∑ld,ud1∈Ld,ld≠ud1tldud1v×xldud1v≤Tv,∀v∈V,(56)∑lp,up1∈Lp,lp≠up1tlpup1v×xlpup1v≤Tv,∀v∈V.


#### 3.3.3. Constraints for the First and the Second Echelon

Constraints (57) and (58) are used to connect the transportation of products in the first echelon and the product delivery and pickup process in the second echelon. Constraint (57) specifies the product quantity transported from facility m to DC i, and the delivery quantity is equal to the total demands of the customers served by DC i. Constraint (58) defines the product quantity transported from logistics facility m to PC j, and the pickup quantity is equal to the total demands of the customers served by PC j. Constraints (59) and (60) define that all vehicles must serve DCs and PCs within their time windows, respectively. Constraints (61) and (62) represent that the vehicle arrives at the DCs and PCs, respectively. Constraints (63) and (64) denote that each vehicle departure from and return to the LC meets the time windows of the LC:(57)∑m∈M∑s∈Sqmis=∑v∈V∑ud2∈Ud2qiud2v×xiud2v,∀i∈Iandi≠m,(58)∑m∈M∑s∈Sqmjs=∑v∈V∑up2∈Up2qjup2v×xjup2v,∀j∈Jandj≠m,(59)αi≤ativ≤βi,∀v∈V,(60)αj≤atjv≤βj,∀v∈V,(61)dtud2v+tud2iv≤ativ,∀ud2∈Ud2,∀i∈I,and∀v∈V,(62)dtup2v+tup2jv≤atjv,∀up2∈Up2,∀j∈J,and∀v∈V,(63)αl≤dtldv≤βl,∀l∈L,ld=l0,and∀v∈V,(64)αl≤drlv≤βl,∀l∈Land∀v∈V.


## 3.1. Problem Statement

Collaboration between logistics facilities in the network contributes to resource integration and optimization so as to minimize the overall operating cost and maximize the revenue of individual logistics facilities. In this study, the 2E-CMDPDTW consists of different logistics facilities containing a logistics center (LC), multiple pickup centers (PCs), and multiple distribution centers (DCs) as the first echelon and a certain number of delivery and pickup customers (Cs) as the second echelon in the 2E-LPDN. Of these facilities, the LC is responsible for collecting and distributing products; and, PCs and DCs provide pickup and delivery service for their customers, respectively. As for the transportation resources, transportation tasks between different logistics facilities in the first-echelon are carried out by semitrailer trucks. For the second-echelon transportation service between a logistics facility and its customers, vehicles are utilized to provide door-to-door services. Fixed time windows are requested by all facilities and customers in the network. If customers were not served within their time windows, the network operating cost would increase along with the decrease in service efficiency.The structures of the 2E-CMDPDTW optimization before the collaboration and after the collaboration are illustrated in Figure1. Without collaboration, logistics facilities operate independently. In comparison, collaborating logistics facilities in the alliance can share customers and transportation resources. In Figure 1, the numbers attached to the links represent the travel time of semitrailer trucks or vehicles on the routes.Figure 1
The 2E-LPDN before and after collaboration. (a) Two-echelon noncollaborative pickup and delivery network. (b) Two-echelon collaborative pickup and delivery network.
(a)(b)In the noncollaborative 2E-LPDN before optimization (Figure1(a)), the LC serves DCs, PCs, and its customers with various demands simultaneously; however, DCs and PCs only serve their customers with delivery and pickup demands, respectively. Three problems arise consequently. (1) There are unreasonable transportation routes (e.g., long-haul and cross transportation). For example, C2 is closer to the LC but served by DC1. This situation can be observed in practice in terms of customer loyalty and long-term contractual partnership. (2) Customers’ time windows are frequently violated. The vehicle, for instance, starting from the LC must wait for open time windows when arriving at C13 and C14. The demand of C12 cannot be satisfied on time since the arrival of the vehicle from DC2 is beyond the given time window. (3) Transportation resources are underutilized. For example, the demand for DC1 served by the LC can be satisfied with a vehicle rather than a semitrailer truck. In other words, other facilities can also be served with the truck. In response to the aforementioned problems, collaboration with resource sharing is effective for efficiency improvement and reducing resource underutilization.Figure1(b) presents the collaborative 2E-LPDN after optimization. The LC plays an important role in organizing delivery tasks and pickup tasks and coordinating logistics facilities within the network. The advantages of collaboration are concluded as follows. First of all, with customer information sharing, the customers in the alliance can be reassigned to new service providers closer to them. In addition, the transportation routes both in the first and second echelon can be optimized to guarantee the given service time windows. Finally, sharing transportation resources, such as vehicles and semitrailer trucks, can effectively avoid idle transportation resources in the collaborative logistics network. For example, the customers re-arranged in the route LC ⟶ C21 ⟶ C22 ⟶ C23 ⟶ C24 ⟶ LC, and the other route PC2 ⟶ C32 ⟶ C33 ⟶ PC2 can be served orderly by the two vehicles on a single transportation task. As a consequence, the aim of reducing the total cost and the number of vehicles is achieved in the optimized 2E-CMDPDTW.Table2 shows the comparative results of the costs, number of trucks, and number of vehicles before and after the 2E-CMDPDTW optimization. Assume that the transportation cost per unit time of the semitrailer trucks equals $30, and all of the delivery cost and pickup cost of vehicles equal $15. The unit time penalty cost of time windows violation is $10, and each of DCs and PCs pays for a fixed collaborative cost of $100 when deciding to collaborate.Table 2
Comparison before and after collaboration in the two-echelon pickup and delivery network.
CaseTransportation cost ($)Pickup and delivery cost ($)Penalty cost ($)Collaborative cost ($)Total operating cost ($)Time windows violationNumber of trucksNumber of vehiclesBefore collaboration660132070020505412After collaboration450825050017750111In Table2, the total operating cost of the optimized 2E-CMDPDTW is decreased by $275, and 3 semitrailer trucks and 1 vehicle are saved. In the network before collaboration, the number of violating time windows equals 5, while all the customers’ time windows are satisfied in the collaborative one. The results reflect that the total operating cost, the number of trucks, the number of vehicles, and the number of violating time windows in the 2E-CMDPDTW are reduced. Therefore, resource sharing can reduce operating costs and save transportation resources in the entire network.


## 3.2. Definitions

In this section, the definition, notations, and decision variables used in this study for the 2E-CMDPDTW optimization are introduced. The definitions and the decision variables are listed in Tables3 and 4, respectively.Table 3
Definitions and notations related to the mathematical model.
SetsDefinitionsLSet of LC,L=L0,l∈LISet of DCs,I=i|i=1,2,3,…,hJSet of PCs,J=j|j=1,2,3,…,h′UdSet of all delivery customers in the collaborative 2E-LPDN,Ud=ud|ud=1,2,3,…,uUpSet of all pickup customers in the collaborative 2E-LPDN,Up=up|up=1,2,3,…,u′MSet of all logistics facilities, containing the LC, DCs, and PCs in the collaborative 2E-LPDN,M=m|m=1,2,3,…,h,h+1,…,h+h′,h+h′+1Ud1Set of delivery customers served by the LC,Ud1=ud1|ud1=1,2,3,…,u1Ud2Set of delivery customers served by the DCs,Ud2=ud2|ud2=1,2,3,…,u2Up1Set of pickup customers served by the LC,Up1=up1|up1=1,2,3,…,u1′Up2Set of pickup customers served by the PCs,Up2=up2|up2=1,2,3,…,u2′DSet of the DCs and their delivery customers,D=I∪Ud2=d|d=1,2,3,…,h,h+1,…,h+u2PSet of the PCs and their pickup customers,P=J∪Up2=p|p=1,2,3,…,h′,h′+1,…,h′+u2′LdSet of the LC and its delivery customers,Ld=L∪Ud1=ld|ld=1,2,3,…,u1+1LpSet of the LC and its pickup customers,Lp=L∪Up1=lp|lp=1,2,3,…,u2+1SSet of semitrailer trucks used for product transportation among LC, DCs, and PCs,S=s|s=1,2,3,…,s′VSet of vehicles used for product transportation between logistics facilities and customers,V=v|v=1,2,3,…,v′OsSet of the DCs and the PCs served by semitrailer trucks, s∈SSivdSet of delivery customers whose delivery service are provided by vehiclev from DC i, Sivd⊆Ud2, i∈I and v∈VSjvpSet of pickup customers whose pickup service are provided by vehiclev from PC j, Sjvp⊆Up2, j∈J and v∈VSlvdSet of delivery customers whose delivery service are provided by vehiclev from LC l, Slvd⊆Ud1, l∈L and v∈VSlvpSet of pickup customers whose pickup service are provided by vehiclev from LC l, Slvp⊆Up1, l∈L and v∈VParametersLmm′sDistance of semitrailer trucks starts from logistics facility m to m′, s∈S, mandm′∈M,andm≠m′Ldud2vDistance of vehiclev starts from the DC or delivery customer d to ud2, v∈V, d∈D, and ud2∈Ud2Lpup2vDistance of vehiclev starts from the PC or pickup customer p to up2, v∈V, p∈P, and up2∈Up2Lldud1vDistance of vehiclev starts from the LC or delivery customer ld to ud1, v∈V, ld∈Ld, and ud1∈Ud1Llpup1vDistance of vehiclev starts from the LC or pickup customer lp to up1, v∈V, lp∈Lp, and up1∈Up1fsDiesel fuel consumption rate per mile for semitrailer trucks (gallon/mile), s∈SfvPetrol fuel consumption rate per mile for vehiclev (gallon/mile), v∈VpsUnit price of diesel used in semitrailer trucks (dollar/gallon)pvUnit price of petrol used in vehicles (dollar/gallon)MsAverage annual cost for the maintenance of semitrailer trucksMvAverage annual cost for the maintenance of vehiclesBNumber of working periods per yearQiDemand quantity of DCi for delivery, i∈IQjDemand quantity of PCj for pickup, j∈IRud1Demand quantity of delivery customerud1 served by the LC, ud1∈Ud1Rup1Demand quantity of pickup customerup1 served by the LC, up1∈Up1Rud2Demand quantity of delivery customerud2 served by the DCs, ud2∈Ud2Rup2Demand quantity of pickup customerup2 served by the PCs, up2∈Up2KsApproved load capacity of semitrailer trucks, s∈SKvApproved load capacity of vehiclev, v∈VFmFixed cost of logistics facilitym, m∈MθVariable cost coefficient for products’ transportationcmDiscount on the costs of logistics facilitym, if collaboration occurs (dollar), m∈MOsNumber of DCs and PCs that provides transportation services by semitrailer trucks, s∈SSivdNumber of delivery customers whose delivery service are provided by vehiclev from DC i, i∈I and v∈VSjvpNumber of pickup customers whose pickup service are provided by vehiclev from PC j, j∈J and v∈VSlvdNumber of delivery customers whose delivery service are provided by vehiclev from LC l, l∈L and v∈VSlvpNumber of pickup customers whose pickup service are provided by vehiclev from LC l, l∈L and v∈Vαi,βiTime window of the DCi, i∈Iαj,βjTime window of the PCj, j∈Jαud,βudTime window of delivery customerud, ud∈Udαup,βupTime window of pick customerup, up∈UpTsMaximum travel time allowed for semitrailer trucks on every route, s∈STvMaximum travel time allowed for vehiclev on every route, v∈VweThe unit penalty cost for semitrailer trucks and vehicles arriving earlywdThe unit penalty cost for semitrailer trucks and vehicles arriving delayatisTime of semitrailer trucks arrival at the DC i, s∈S and i∈IatjsTime of semitrailer trucks arrival at the PC j, s∈S and j∈JativTime of vehiclev arrival at DC i, v∈V and i∈IatjvTime of vehiclev arrival at PC j, v∈V and j∈JatudvTime of vehiclev arrival at delivery customer ud, v∈V and ud∈UdatupvTime of vehiclev arrival at pickup customer up, v∈V and up∈Upatud2vTime of vehiclev arrival at delivery customer ud2, v∈V and ud2∈Ud2atup2vTime of vehiclev arrival at pickup customer up2, v∈V and up2∈Up2atud1vTime of vehiclev arrival at delivery customer ud1, v∈V and ud1∈Ud1atup1vTime of vehiclev arrival at pickup customer up1, v∈V and up1∈Up1dtmsTime of semitrailer trucks depart from logistics facility m, m∈M and s∈SdtdvTime of vehiclev depart from DC i or delivery customer ud2, v∈V, i∈I, and ud2∈Ud2dtpvTime of vehiclev depart from PC j or pickup customer up2, v∈V, j∈J, and up2∈Up2dtldvTime of vehiclev depart from LC l or delivery customer ud1, v∈V, ld∈Ld, and ud1∈Ud1dtlpvTime of vehiclev depart from LC l or pickup customer up1, v∈V, lp∈Lp, and up1∈Up1dtud2vTime of vehiclev depart from delivery customer ud2, v∈V and ud2∈Ud2dtup2vTime of vehiclev depart from pickup customer up2, v∈V and up2∈Up2drlsReturn time of semitrailer trucks to LC l, l∈L and s∈SdrlvReturn time of vehiclev to LC l, l∈L and v∈Vtmm′sTime of semitrailer trucks travel from logistics facility m to m′, m∈M and s∈StmisTime of semitrailer trucks travel from logistics facility m to DC i, m∈M, i∈I, and s∈StmjsTime of semitrailer trucks travel from logistics facility m to PC j, m∈M, i∈I, and s∈Stdud2vTime of vehiclev travel from DC i or ud2 to ud2, v∈V, d∈D, i∈I, and ud2∈Ud2tpup2vTime of vehiclev travel from PC j or up2 to up2, v∈V, p∈P, j∈J, and up2∈Up2tldud1vTime of vehiclev travel from LC l or ud1 to ud1, v∈V, ld∈Ld, and ud1∈Ud1tlpup1vTime of vehiclev travel from LC l or up1 to up1, v∈V, lp∈Lp, and up1∈Up1tud2ivTime of vehiclev travel from ud2 to DC i, v∈V, i∈I, and ud2∈Ud2tup2jvTime of vehiclev travel from up2 to PC j, v∈V, j∈J, and up2∈Up2Table 4
Decision variables related to the mathematical model.
VariablesDefinitionsxmm′sxmm′s=1, if semitrailer truck s starts from logistics facility m to m′; otherwise, xmm′s=0, mandm′∈M and s∈Sxmisxmis=1, if semitrailer truck s starts from logistics facility m to DC i; otherwise, xmis=0, m∈M and s∈Sxmjsxmjs=1, if semitrailer truck s starts from logistics facility m to PC j; otherwise, xmjs=0, m∈M and s∈Sxdud2vxdud2v=1, if vehicle v starts from d to ud2; otherwise, xdud2v=0, d∈D, v∈V, and ud2∈Ud2xpup2vxpup2v=1, if vehicle v starts from p to up2; otherwise, xpup2v=0, p∈P, v∈V, and up2∈Up2xldud1vxldud1v=1, if vehicle v starts from ld to ud1; otherwise, xldud1v=0, ld∈Ld, ud1∈Ud1, and v∈Vxlpup1vxlpup1v=1, if vehicle v starts from lp to up1; otherwise, xlpup1v=0, lp∈Lp, up1∈Up1, and v∈Vxmuvxmuv=1, if vehicle v starts from logistics facility m to customer u; otherwise, xmuv=0, m∈M, and u∈Ud∪UpZmZm=1, if logistics facility m agrees to join the alliance to cooperate; otherwise, Zm=0, m∈MqmisQuantity of products transported by semitrailer trucks from logistics facility m to DC i, m∈M, s∈S, and i∈IqmjsQuantity of products transported by semitrailer trucks from logistics facility m to PC j, m∈M, s∈S, and j∈JqmudvQuantity of products delivered by vehiclev from logistics facility m to ud, v∈V, m∈M, and ud∈UdqmupvQuantity of products collected by vehiclev from logistics facility m to up, v∈V, m∈M, and up∈Upqiud2vQuantity of products delivered by vehiclev from DC i to ud2, v∈V, i∈I, and ud2∈Ud2qjup2vQuantity of products collected by vehiclev from PC j to up2, v∈V, j∈J, and up2∈Up2qdud2vQuantity of products delivered by vehiclev from DC i or ud2 to ud2, v∈V, d∈D, i∈I, and ud2∈Ud2qpup2vQuantity of products collected by vehiclev from PC j or up2 to up2, v∈V, p∈P, j∈J, and up2∈Up2qldud1vQuantity of products delivered by vehiclev from LC l or ud1 to ud1, v∈V, ld∈Ld, and ud1∈Ud1qlpup1vQuantity of products collected by vehiclev from LC l or up1 to up1, v∈V, lp∈Lp, and up1∈Up1


## 3.3. Mathematical Model

In this study, a bi-objective mathematical model is developed to optimize the total operating cost of the network and the number of vehicles in the 2E-CMDPDTW. The objective functions in equations (1) and (2) are as follows:(1)MinF1=TC1+TC2+PC1+PC2+TC3,(2)MinF2=∑m∈M∑v∈V∑u∈Ud∪Upxmuv.TC1 expresses the total operating cost of semitrailer trucks, and it includes the total transportation cost and the total maintenance cost in alliances:(3)TC1=∑m,m′∈M,m≠m′∑s∈Sxmm′s×Lmm′s×fs×ps+maxs∈S∑m∈L∑i∈Ixmis,∑m∈L∑j∈Jxmjs×MsB,where ∑m,m′∈M,m≠m′∑s∈Sxmm′s×Lmm′s×fs×ps represents the transportation cost of semitrailer trucks and maxs∈S∑m∈L∑i∈Ixmis,∑m∈L∑j∈Jxmjs×Ms/B depicts the maintenance cost of semitrailer trucks.TC2 expresses the total delivery cost, pickup cost, and maintenance cost of vehicles in the second echelon in a working period:(4)TC2=∑d,ud2∈Dd≠ud2∑v∈Vxdud2v×Ldud2v×fv×pv+∑m∈I∑u∈Ud2∑v∈Vxmuv×MvB+∑p,up2∈Pp≠up2∑v∈Vxpup2v×Lpup2v×fv×p+∑m∈J∑u∈Up2∑v∈Vxmuv×MvB+∑ld,ud1∈Ldld≠ud1∑v∈Vxldud1v×Lldud1v×fv×pv+∑m∈L∑u∈Ud1∑v∈Vxmuv×MvB+∑lp,up1∈Lplp≠up1∑v∈Vxlpup1v×Llpup1v×fv×pv+∑m∈L∑u∈Up1∑v∈Vxmuv×MvB,where ∑d,ud2∈Dd≠ud2∑v∈Vxdud2v×Ldud2v×fv×pv+∑m∈I∑u∈Ud2∑v∈Vxmuv×Mv/B is the delivery cost and maintenance cost of delivery vehicles starting from DCs, ∑p,up2∈Pp≠up2∑v∈Vxpup2v×Lpup2v×fv×p+∑m∈J∑u∈Up2∑v∈Vxmuv×Mv/B represents the pickup cost and maintenance cost of pickup vehicles starting from PCs, ∑ld,ud1∈Ldld≠ud1∑v∈Vxldud1v×Lldud1v×fv×pv+∑m∈L∑u∈Ud1∑v∈Vxmuv×Mv/B depicts the delivery cost and maintenance cost of delivery vehicles from the LC, and ∑lp,up1∈Lplp≠up1∑v∈Vxlpup1v×Llpup1v×fv×pv+∑m∈L∑u∈Up1∑v∈Vxmuv×Mv/B denotes the pickup cost and maintenance cost of pickup vehicles from the LC.PC1 expresses the penalty cost of the semitrailer trucks for early and delayed arrival at the DCs and PCs:(5)PC1=∑m∈M∑i∈I∑s∈Sqmis×we×maxαi−atis,0+∑m∈M∑i∈I∑s∈Sqmis×wd×maxatis−βi,0+∑m∈M∑j∈J∑s∈Sqmjs×we×maxαj−atjs,0+∑m∈M∑j∈J∑s∈Sqmjs×wd×maxatjs−βj,0,where ∑m∈M∑i∈I∑s∈Sqmis×we×maxαi−atis,0 and ∑m∈M∑i∈I∑s∈Sqmis×wd×maxatis−βi,0 indicate the penalty cost of earliness and delay of the semitrailer trucks for DCs and ∑m∈M∑j∈J∑s∈Sqmjs×we×maxαj−atjs,0 and ∑m∈M∑j∈J∑s∈Sqmjs×wd×maxatjs−βj,0 denote the penalty cost of earliness and delay of the semitrailer trucks for PCs.PC2 represents the penalty cost of the vehicles for early and delayed arrival at the delivery customers and pickup customers:(6)PC2=∑m∈M∪Ud∑ud∈Ud∑v∈Vqmudv×we×maxaud−atudv,0+∑m∈M∪Ud∑ud∈Ud∑v∈Vqmudv×wd×maxatudv−βud,0+∑m∈M∪Up∑up∈Up∑v∈Vqmupv×we×maxaup−atupv,0+∑m∈M∪Up∑up∈Up∑v∈Vqmupv×wd×maxatupv−βup,0,where ∑m∈M∪Ud∑ud∈Ud∑v∈Vqmudv×we×maxaud−atudv,0 and ∑m∈M∪Ud∑ud∈Ud∑v∈Vqmudv×wd×maxatudv−βud,0 depict the penalty cost of earliness and delay of the vehicles for delivery customers and ∑m∈M∪Up∑up∈Up∑v∈Vqmupv×we×maxaup−atupv,0 and ∑m∈M∪Up∑up∈Up∑v∈Vqmupv×wd×maxatupv−βup,0 represent the penalty cost of earliness and delay of the vehicles for pickup customers.TC3 expresses the total cost of logistics facilities including the fixed and variable costs and the collaborative discount for establishing alliances between logistics facilities:(7)TC3=∑m∈MFm+∑i∈I∑m∈Mθqmis+∑j∈J∑m∈Mθqmjs−∑m∈MZmcm,where ∑m∈MFm represents the total fixed cost of the LC, DCs, and PCs, ∑i∈I∑m∈Mθqmis and ∑j∈J∑m∈Mθqmjs depict the total variable cost of DCs and PCs, respectively, and ∑m∈MZmcm denotes the collaboration discount for establishing alliances in terms of government subsidies.The constraints of the model are divided into three parts such as constraints for the first echelon, constraints for the second echelon, and constraints for the first and the second echelons.


### 3.3.1. Constraints for the First Echelon

Constraint (8) ensures that, if m′=L0, a semitrailer truck starts from the LC and returns to the LC after the task is completed. Otherwise, if m′ is a DC or PC, this constraint imposes the balance of semitrailer trucks entering and leaving the logistics facility. Constraints (9) and (10) guarantee that each DC and PC only can be served by one semitrailer truck, respectively. Constraint (11) eliminates subtours of every route in the first echelon. Constraints (12) and (13) express that the demands of each DC and PC can be satisfied, respectively. Constraints (14) and (15) ensure that the maximum load of semitrailer trucks is able to satisfy the total demand of the logistics facilities which are served by the semitrailer trucks. Constraints (16) and (17) imply that semitrailer trucks serve DCs and PCs within their time windows, respectively. Constraint (18) represents that the semitrailer truck arrives at the DC. Constraint (19) defines that the semitrailer truck arrives at the PC. Constraints (20) and (21) denote that the time of each semitrailer truck departure from and return to the logistics facilities meet the time window of the LC. Constraint (22) limits that the total time a semitrailer truck travels on each route must be less than or equal to the allowed time of each route:(8)∑m∈M,m≠m′xmm′s−∑m∈M,m≠m′xm′ms=0,∀m′∈Mand∀s∈S,(9)∑m∈M,m≠i∑s∈Sxmis≤1,∀i∈I,(10)∑m∈M,m≠j∑s∈Sxmjs≤1,∀j∈J,(11)∑m∈Os∑m′∈Os,m′≠mxmm′s=Os−1,∀s∈SandOs≥2,(12)∑s∈S∑m∈Mqmis×xmis≥Qi,∀i∈Iandi≠m,(13)∑s∈S∑m∈Mqmjs×xmjs≥Qi,∀j∈Jandj≠m,(14)∑i∈I∑m∈Mqmis×xmis≤Ks,∀s∈S,(15)∑j∈J∑m∈Mqmjs×xmjs≤Ks,∀s∈S,(16)αi≤atis≤βi,∀i∈Iand∀s∈S,(17)αj≤atjs≤βj,∀j∈Jand∀s∈S,(18)dtms+tmis≤atis,∀m∈M,∀i∈I,and∀s∈S,(19)dtms+tmjs≤atjs,∀m∈M,∀j∈J,and∀s∈S,(20)αl≤dtms≤βl,∀l∈L,∀s∈S,andm=l0,(21)αl≤drls≤βl,∀l∈L,∀s∈S,andm=l0,(22)∑m,m′∈M,m≠m′tmm′s×xmm′s≤Ts,s∈S.


### 3.3.2. Constraints for the Second Echelon

Constraints (23) and (24) ensure that a delivery vehicle serving DCs’ customers originates from a DC to serve customers and returns to the DC when it completed a task. Constraint (25) imposes the balance of delivery vehicles entering and leaving the delivery customers served by DCs. Constraints (26) and (27) guarantee that a pickup vehicle serving PCs’ customers originates from a PC and returns to the PC after completing a task. Constraint (28) imposes the balance of pickup vehicles entering and leaving the pickup customers served by PCs. Constraints (29) and (30) express that a delivery vehicle serving the LC originates from the LC and returns to the LC. Constraint (31) imposes the balance of delivery vehicles entering and leaving the delivery customers served by the LC. Constraints (32) and (33) imply that a pickup vehicle serving the LC originates from the LC and returns to the LC. Constraint (34) imposes the balance of pickup vehicles entering and leaving the pickup customers served by the LC. Constraints (35)–(38) are used to eliminate subtours of every route in the second echelon. Constraints (39)–(42) are used to guarantee the demands of each delivery or pickup customer can be satisfied. Constraints (43)–(46) stipulate that the maximum load of each vehicle can satisfy the total demand of customers served by this vehicle. Constraints (47) and (48) represent that all vehicles serve delivery and pickup customers within their time windows. Constraints (49)–(52) define that the vehicle arrives at the delivery and pickup customers. Constraints (53)–(56) limit that the total travel time for a vehicle must be less than or equal to the allowed time of each route:(23)∑ud2∈Ud2xdud2v=1,∀d∈Iand∀v∈V,(24)∑ud2∈Ud2xud2dv=1,∀d∈Iand∀v∈V,(25)∑ud2∈Dxdud2v−∑ud2∈Dxud2dv=0,∀d∈D,d≠ud2,and∀v∈V,(26)∑up2∈Up2xpup2v=1,∀p∈Jand∀v∈V,(27)∑up2∈Up2xup2pv=1,∀p∈Jand∀v∈V,(28)∑up2∈Pxpup2v−∑up2∈Pxup2pv=0,∀p∈P,p≠up2,and∀v∈V,(29)∑ud1∈Ud1xldud1v=1,∀ld∈Ldand∀v∈V,(30)∑ud1∈Ud1xud1ldv=1,∀ld∈Ldand∀v∈V,(31)∑ud1∈Ldxldud1v−∑ud1∈Ldxud1ldv=0,∀ld∈Ld,ld≠ud1,and∀v∈V,(32)∑up1∈Up1xlpup1v=1,∀lp∈Lpand∀v∈V,(33)∑up1∈Up1xup1lpv=1,∀lp∈Lpand∀v∈V,(34)∑up1∈Lpxlpup1v−∑up1∈Lpxup1lpv=0,∀lp∈Lp,lp≠up1,and∀v∈V,(35)∑d∈Sivd∑ud2∈Sivdxdud2v=Sivd−1,∀i∈I,d≠ud2,∀v∈V,Sivd⊆Ud2,andSivd≥2,(36)∑p∈Sjvp∑up2∈Sjvpxpup2v=Sjvp−1,∀j∈J,p≠up2,∀v∈V,Sjvp⊆Up2,andSjvp≥2,(37)∑ld∈Slvd∑ud1∈Slvdxldud1v=Slvd−1,∀l∈L,ld≠ud1,∀v∈V,Slvd⊆Ud1,andSlvd≥2,(38)∑lp∈Slvp∑up1∈Slvpxlpup1v=Slvp−1,∀l∈L,lp≠up1,∀v∈V,Slvp⊆Up1,andSlvp≥2,(39)∑d∈D∑v∈Vqdud2v×xdud2v≥Rud2,∀ud2∈Ud2andd≠ud2,(40)∑p∈P∑v∈Vqpup2v×xpup2v≥Rup2,∀up2∈Up2andp≠up2,(41)∑ld∈Ld∑v∈Vqldud1v×xldud1v≥Rud1,∀ud1∈Ud1andld≠ud1,(42)∑lp∈Lp∑v∈Vqlpup1v×xlpup1v≥Rup1,∀up1∈Up1andlp≠up1,(43)∑d∈D∑ud2∈Ud2,d≠ud2qdud2v×xdud2v≤Kv,∀v∈V,(44)∑p∈P∑up2∈Up2,p≠up2qpup2v×xpup2v≤Kv,∀v∈V,(45)∑ud1∈Ud1∑ld∈Ld,ld≠ud1qldud1v×xldud1v≤Kv,∀v∈V,(46)∑up1∈Up1∑lp∈Lp,lp≠up1qlpup1v×xlpup1v≤Kv,∀v∈V,(47)αud≤atudv≤βud,∀ud∈Udand∀v∈V,(48)αup≤atupv≤βup,∀up∈Upand∀v∈V,(49)dtdv+tdud2v≤atud2v,∀ud2∈Ud2,ud2≠d,∀d∈D,and∀v∈V,(50)dtpv+tpup2v≤atup2v,∀up2∈Up2,up2≠p,∀p∈P,and∀v∈V,(51)dtldv+tldud1v≤atud1v,∀ud1∈Ud1,ud1≠ld,∀ld∈Ld,and∀v∈V,(52)dtlpv+tlpup1v≤atup1v,∀up1∈Up1,up1≠lp,∀lp∈Lp,and∀v∈V,(53)∑d,ud2∈D,d≠ud2tdud2v×xdud2v≤Tv,∀v∈V,(54)∑p,up2∈P,p≠up2tpup2v×xpup2v≤Tv,∀v∈V,(55)∑ld,ud1∈Ld,ld≠ud1tldud1v×xldud1v≤Tv,∀v∈V,(56)∑lp,up1∈Lp,lp≠up1tlpup1v×xlpup1v≤Tv,∀v∈V.


### 3.3.3. Constraints for the First and the Second Echelon

Constraints (57) and (58) are used to connect the transportation of products in the first echelon and the product delivery and pickup process in the second echelon. Constraint (57) specifies the product quantity transported from facility m to DC i, and the delivery quantity is equal to the total demands of the customers served by DC i. Constraint (58) defines the product quantity transported from logistics facility m to PC j, and the pickup quantity is equal to the total demands of the customers served by PC j. Constraints (59) and (60) define that all vehicles must serve DCs and PCs within their time windows, respectively. Constraints (61) and (62) represent that the vehicle arrives at the DCs and PCs, respectively. Constraints (63) and (64) denote that each vehicle departure from and return to the LC meets the time windows of the LC:(57)∑m∈M∑s∈Sqmis=∑v∈V∑ud2∈Ud2qiud2v×xiud2v,∀i∈Iandi≠m,(58)∑m∈M∑s∈Sqmjs=∑v∈V∑up2∈Up2qjup2v×xjup2v,∀j∈Jandj≠m,(59)αi≤ativ≤βi,∀v∈V,(60)αj≤atjv≤βj,∀v∈V,(61)dtud2v+tud2iv≤ativ,∀ud2∈Ud2,∀i∈I,and∀v∈V,(62)dtup2v+tup2jv≤atjv,∀up2∈Up2,∀j∈J,and∀v∈V,(63)αl≤dtldv≤βl,∀l∈L,ld=l0,and∀v∈V,(64)αl≤drlv≤βl,∀l∈Land∀v∈V.


## 3.3.1. Constraints for the First Echelon

Constraint (8) ensures that, if m′=L0, a semitrailer truck starts from the LC and returns to the LC after the task is completed. Otherwise, if m′ is a DC or PC, this constraint imposes the balance of semitrailer trucks entering and leaving the logistics facility. Constraints (9) and (10) guarantee that each DC and PC only can be served by one semitrailer truck, respectively. Constraint (11) eliminates subtours of every route in the first echelon. Constraints (12) and (13) express that the demands of each DC and PC can be satisfied, respectively. Constraints (14) and (15) ensure that the maximum load of semitrailer trucks is able to satisfy the total demand of the logistics facilities which are served by the semitrailer trucks. Constraints (16) and (17) imply that semitrailer trucks serve DCs and PCs within their time windows, respectively. Constraint (18) represents that the semitrailer truck arrives at the DC. Constraint (19) defines that the semitrailer truck arrives at the PC. Constraints (20) and (21) denote that the time of each semitrailer truck departure from and return to the logistics facilities meet the time window of the LC. Constraint (22) limits that the total time a semitrailer truck travels on each route must be less than or equal to the allowed time of each route:(8)∑m∈M,m≠m′xmm′s−∑m∈M,m≠m′xm′ms=0,∀m′∈Mand∀s∈S,(9)∑m∈M,m≠i∑s∈Sxmis≤1,∀i∈I,(10)∑m∈M,m≠j∑s∈Sxmjs≤1,∀j∈J,(11)∑m∈Os∑m′∈Os,m′≠mxmm′s=Os−1,∀s∈SandOs≥2,(12)∑s∈S∑m∈Mqmis×xmis≥Qi,∀i∈Iandi≠m,(13)∑s∈S∑m∈Mqmjs×xmjs≥Qi,∀j∈Jandj≠m,(14)∑i∈I∑m∈Mqmis×xmis≤Ks,∀s∈S,(15)∑j∈J∑m∈Mqmjs×xmjs≤Ks,∀s∈S,(16)αi≤atis≤βi,∀i∈Iand∀s∈S,(17)αj≤atjs≤βj,∀j∈Jand∀s∈S,(18)dtms+tmis≤atis,∀m∈M,∀i∈I,and∀s∈S,(19)dtms+tmjs≤atjs,∀m∈M,∀j∈J,and∀s∈S,(20)αl≤dtms≤βl,∀l∈L,∀s∈S,andm=l0,(21)αl≤drls≤βl,∀l∈L,∀s∈S,andm=l0,(22)∑m,m′∈M,m≠m′tmm′s×xmm′s≤Ts,s∈S.


## 3.3.2. Constraints for the Second Echelon

Constraints (23) and (24) ensure that a delivery vehicle serving DCs’ customers originates from a DC to serve customers and returns to the DC when it completed a task. Constraint (25) imposes the balance of delivery vehicles entering and leaving the delivery customers served by DCs. Constraints (26) and (27) guarantee that a pickup vehicle serving PCs’ customers originates from a PC and returns to the PC after completing a task. Constraint (28) imposes the balance of pickup vehicles entering and leaving the pickup customers served by PCs. Constraints (29) and (30) express that a delivery vehicle serving the LC originates from the LC and returns to the LC. Constraint (31) imposes the balance of delivery vehicles entering and leaving the delivery customers served by the LC. Constraints (32) and (33) imply that a pickup vehicle serving the LC originates from the LC and returns to the LC. Constraint (34) imposes the balance of pickup vehicles entering and leaving the pickup customers served by the LC. Constraints (35)–(38) are used to eliminate subtours of every route in the second echelon. Constraints (39)–(42) are used to guarantee the demands of each delivery or pickup customer can be satisfied. Constraints (43)–(46) stipulate that the maximum load of each vehicle can satisfy the total demand of customers served by this vehicle. Constraints (47) and (48) represent that all vehicles serve delivery and pickup customers within their time windows. Constraints (49)–(52) define that the vehicle arrives at the delivery and pickup customers. Constraints (53)–(56) limit that the total travel time for a vehicle must be less than or equal to the allowed time of each route:(23)∑ud2∈Ud2xdud2v=1,∀d∈Iand∀v∈V,(24)∑ud2∈Ud2xud2dv=1,∀d∈Iand∀v∈V,(25)∑ud2∈Dxdud2v−∑ud2∈Dxud2dv=0,∀d∈D,d≠ud2,and∀v∈V,(26)∑up2∈Up2xpup2v=1,∀p∈Jand∀v∈V,(27)∑up2∈Up2xup2pv=1,∀p∈Jand∀v∈V,(28)∑up2∈Pxpup2v−∑up2∈Pxup2pv=0,∀p∈P,p≠up2,and∀v∈V,(29)∑ud1∈Ud1xldud1v=1,∀ld∈Ldand∀v∈V,(30)∑ud1∈Ud1xud1ldv=1,∀ld∈Ldand∀v∈V,(31)∑ud1∈Ldxldud1v−∑ud1∈Ldxud1ldv=0,∀ld∈Ld,ld≠ud1,and∀v∈V,(32)∑up1∈Up1xlpup1v=1,∀lp∈Lpand∀v∈V,(33)∑up1∈Up1xup1lpv=1,∀lp∈Lpand∀v∈V,(34)∑up1∈Lpxlpup1v−∑up1∈Lpxup1lpv=0,∀lp∈Lp,lp≠up1,and∀v∈V,(35)∑d∈Sivd∑ud2∈Sivdxdud2v=Sivd−1,∀i∈I,d≠ud2,∀v∈V,Sivd⊆Ud2,andSivd≥2,(36)∑p∈Sjvp∑up2∈Sjvpxpup2v=Sjvp−1,∀j∈J,p≠up2,∀v∈V,Sjvp⊆Up2,andSjvp≥2,(37)∑ld∈Slvd∑ud1∈Slvdxldud1v=Slvd−1,∀l∈L,ld≠ud1,∀v∈V,Slvd⊆Ud1,andSlvd≥2,(38)∑lp∈Slvp∑up1∈Slvpxlpup1v=Slvp−1,∀l∈L,lp≠up1,∀v∈V,Slvp⊆Up1,andSlvp≥2,(39)∑d∈D∑v∈Vqdud2v×xdud2v≥Rud2,∀ud2∈Ud2andd≠ud2,(40)∑p∈P∑v∈Vqpup2v×xpup2v≥Rup2,∀up2∈Up2andp≠up2,(41)∑ld∈Ld∑v∈Vqldud1v×xldud1v≥Rud1,∀ud1∈Ud1andld≠ud1,(42)∑lp∈Lp∑v∈Vqlpup1v×xlpup1v≥Rup1,∀up1∈Up1andlp≠up1,(43)∑d∈D∑ud2∈Ud2,d≠ud2qdud2v×xdud2v≤Kv,∀v∈V,(44)∑p∈P∑up2∈Up2,p≠up2qpup2v×xpup2v≤Kv,∀v∈V,(45)∑ud1∈Ud1∑ld∈Ld,ld≠ud1qldud1v×xldud1v≤Kv,∀v∈V,(46)∑up1∈Up1∑lp∈Lp,lp≠up1qlpup1v×xlpup1v≤Kv,∀v∈V,(47)αud≤atudv≤βud,∀ud∈Udand∀v∈V,(48)αup≤atupv≤βup,∀up∈Upand∀v∈V,(49)dtdv+tdud2v≤atud2v,∀ud2∈Ud2,ud2≠d,∀d∈D,and∀v∈V,(50)dtpv+tpup2v≤atup2v,∀up2∈Up2,up2≠p,∀p∈P,and∀v∈V,(51)dtldv+tldud1v≤atud1v,∀ud1∈Ud1,ud1≠ld,∀ld∈Ld,and∀v∈V,(52)dtlpv+tlpup1v≤atup1v,∀up1∈Up1,up1≠lp,∀lp∈Lp,and∀v∈V,(53)∑d,ud2∈D,d≠ud2tdud2v×xdud2v≤Tv,∀v∈V,(54)∑p,up2∈P,p≠up2tpup2v×xpup2v≤Tv,∀v∈V,(55)∑ld,ud1∈Ld,ld≠ud1tldud1v×xldud1v≤Tv,∀v∈V,(56)∑lp,up1∈Lp,lp≠up1tlpup1v×xlpup1v≤Tv,∀v∈V.


## 3.3.3. Constraints for the First and the Second Echelon

Constraints (57) and (58) are used to connect the transportation of products in the first echelon and the product delivery and pickup process in the second echelon. Constraint (57) specifies the product quantity transported from facility m to DC i, and the delivery quantity is equal to the total demands of the customers served by DC i. Constraint (58) defines the product quantity transported from logistics facility m to PC j, and the pickup quantity is equal to the total demands of the customers served by PC j. Constraints (59) and (60) define that all vehicles must serve DCs and PCs within their time windows, respectively. Constraints (61) and (62) represent that the vehicle arrives at the DCs and PCs, respectively. Constraints (63) and (64) denote that each vehicle departure from and return to the LC meets the time windows of the LC:(57)∑m∈M∑s∈Sqmis=∑v∈V∑ud2∈Ud2qiud2v×xiud2v,∀i∈Iandi≠m,(58)∑m∈M∑s∈Sqmjs=∑v∈V∑up2∈Up2qjup2v×xjup2v,∀j∈Jandj≠m,(59)αi≤ativ≤βi,∀v∈V,(60)αj≤atjv≤βj,∀v∈V,(61)dtud2v+tud2iv≤ativ,∀ud2∈Ud2,∀i∈I,and∀v∈V,(62)dtup2v+tup2jv≤atjv,∀up2∈Up2,∀j∈J,and∀v∈V,(63)αl≤dtldv≤βl,∀l∈L,ld=l0,and∀v∈V,(64)αl≤drlv≤βl,∀l∈Land∀v∈V.


## 4. Hybrid Solution Approach

To solve the developed mathematical model, a hybrid method including the customers clustering, the greedy algorithm, and the Im-NSGA-II is proposed, and a three-stage method is designed for the 2E-CMDPDTW optimization. In the first stage, customer clustering is used to reassign the customers in the 2E-LPDN. Then, based on the clustering results and the greedy algorithm, initial feasible solutions can be generated. In the second stage, the NSGA-II is employed to search for the optimal route of the two-echelon network. Especially, the second-echelon vehicle routes are first optimized, and then, the first-echelon semitrailer truck routes are generated based on the optimal route of vehicles. The first stage and the second stage are conducted using the proposed hybrid heuristic algorithm. In the third stage, the Shapley value method is used for distributing profits to alliance members. Finally, according to the principle of the SMP, the sequence of the logistics facilities joining the alliance is determined, and the optimal scheme for alliance participants to join the alliance is obtained.


### 4.1. Hybrid Heuristics for the 2E-CMDPDTW Optimization

The NSGA-II proposed by Deb et al. [58] has been prominent in academic and practice domains in recent years. It is widely accepted that the NSGA-II can solve problems with multiple objectives effectively and efficiently. In the 2E-CMDPDTW, logistics facilities with different functionalities and amounts of customers with delivery and pickup demands are considered, which increases the computational complexity of the network. Therefore, the customer clustering algorithm and the greedy algorithm are introduced into the hybrid heuristic approach. The two-echelon logistics network is partitioned off as small regions using customer clustering, and the greedy algorithm aims to preliminarily generate the local optimal initial solution of each region [59]. Finally, the improved genetic operators are employed in the NSGA-II to find the global solution. The hybrid method is designed and implemented, as shown in Figure 2.Figure 2
Flowchart of the hybrid heuristic algorithm.The parameters in Figure2 can be explained as follows: k = {k1, k2} represents the number of clusters, in which k1 and k2 are the numbers of cluster centers, pop is the initial population, t and tmax represent the current and the maximum generations, respectively, and run and rmax denote the current and the maximum optimization runs, respectively.


#### 4.1.1. Customer Clustering Algorithm

Customer clustering has been proved to be an effective method to simplify the large-scale and multidepot logistics network [7, 60]. Its basic process refers to reassign customers to logistics facilities closer to them. On one hand, to improve convergence speed and reduce computing efforts, the customer clustering algorithm is employed to generate initial feasible solutions. On the other hand, the customer clustering algorithm plays an important role in avoiding long-haul and cross transportation and reducing transportation costs. Algorithm 1 shows the procedure of the designed clustering algorithm.Algorithm 1: Procedure of the designed customer clustering algorithm.
Customer clustering algorithmInputs: Geographical coordinates of logistics facilities and geographical coordinates of customersOutputs: Clustering results(1)
Step 1: Set the number of LC and DCs k1 occupied in delivery tasks, and k2 indicates the number of LC and PCs handling product pickup demands.(2)
Step 2: Select elements in k1 and k2 as centers of corresponding clusters.(3)
For customers with delivery demands(4)
Step 3: Evaluate the distance between each delivery customer and each cluster center belonging to k1.(5)
For customers with pickup demands(6)
Step 4: Evaluate the distance between each pickup customer and each cluster center belonging to k2.(7)
If the same distances between a certain customer and several cluster centers existStep 5: Assign the customer to one of the cluster centers randomly.ElseStep 6: Assign the customer to the closest cluster center.End(8)
Step 7: Reassess the similarity of customers in each center, and update the clusters.(9)
Step 8: Repeat Steps 3-6 until every customer is assigned to the closest cluster center.(10)
Print the clustering results


#### 4.1.2. Greedy Algorithm

The greedy algorithm referring to an approximate algorithm for optimization problems is usually employed to generate initial solutions [61]. The basic principle of the greedy algorithm is to guarantee the best choice in each search. With the greedy algorithm, this study can generate the local optimal solutions instead of the random ones, thereby improving the convergence speed. The major steps for implementing the greedy algorithm are introduced as follows.Step 1: initialize inputs according to the outputs of customer clustering.Cl1,Cdi2, and Cpj2 represent the sets of customers served by LC l, DC i, and PC j, respectively. Both customers with delivery and pickup demands can be served by the LC; thus, Rd1 and Rp1 represent the delivery and pickup route sets of vehicles from LC l. Similarly, Rdi2 and Rpj2 represent the route sets of vehicles from DC i and PC j.Step 2: define conditions for adding subroutes. A new subroute can be added into the vehicle route if the two customers in the subroute appeared less than twice in the vehicle route, either of the two customers does not appear in the vehicle route, and their time windows will not be violated.Step 3: generate the distance matrix of the customers and sort alternative subroutes. Calculate the distances between the customers with the same kind of requirements, and generate the distance matrix. Then, the lengths of alternative subroutes in the matrix are sorted from the shortest to the longest.Step 4: add the alternative subroutes according to the sorted results. If the alternative subroute fits the conditions above, it will be added to one of the route sets (e.g.,Rd1,Rp1, Rdi2, or Rpj2) according to the customer sets (i.e., Cl1,Cdi2, and Cpj2). In addition, if the customers in the subroutes have appeared in the vehicle route, the subroute will be directly connected to the existing route with the same customer. Repeat this step until the transportation quantity exceeds the maximum capacity of a vehicle.Step 5: generate closed vehicle routes. According to the clustering results, the customers in a certain vehicle route are served by the same facility so that a vehicle route is closed by adding the facility its customers belong to.


#### 4.1.3. The Improved Nondominated Sorting Genetic Algorithm-II

The NSGA-II, a variant of the NSGA, contributes to balancing multiple objectives and generating the optimal decision on optimization problems [62, 63]. Considering the proposed model, this paper improves the method of generating the initial population and using the more appropriate genetic operators based on the original NSGA-II operating process. To improve the performance of global search and maintain the diversity of a population, the genetic operators including order crossover (OX) [45] and polynomial mutation [57] are exploited. The major procedures of implementing the Im-NSGA-II are presented in Figure 3.Figure 3
Major procedures of the Im-NSGA-II.In Figure3, the individuals in the initial population are first generated using the greedy algorithm and the scanning algorithm enabling the diversity of the population. Second, according to the nondomination sorting method [64], the individuals in population Pt are sorted. Each individual in Pt is assigned with a rank value representing the nondominated level. Selection, cross-over, and mutation operators are then employed to generate offspring in Qt. A new population Rt can be generated by combining the individuals in Pt and Qt. Third, a new population Pt+1 is generated by performing nondominated sorting and crowding calculation in Rt. The above process is repeated until the optimal Pareto solution is selected from all the Pareto solutions.Generally, it is hard to pursue all optimal objectives simultaneously in a bi-objective problem since they can be conflicting in nature [65]. According to these objectives, the solutions referring to Pareto solutions cannot be compared simply. Figure 4 shows the optimal Pareto solution of the bi-objective problem. Both individuals A and B are in the optimal Pareto front, and the others dominated by the solutions in the optimal Pareto front are in the feasible region. For the bi-objective problem, when objective 1 (e.g., number of vehicles) is deemed as the selection criterion, individual B is selected as the optimal Pareto solution because of fewer vehicles.Figure 4
Pareto solutions of the multiobjective optimization problem.The Im-NSGA-II is proposed to seek bi-objective optimal Pareto solutions in the 2E-CMDPDTW, as shown in Algorithm2. The processes of nondominated sorting, crowding distance calculating, and genetic operations are presented as follows.Algorithm 2: Procedure of the proposed Im-NSGA-II.
Im-NSGA-II algorithmInputs: Population sizepopsize, maximum number of generations tmax, maximum number of optimization runs rmax, crossover probability pcro, mutation probability pmutOutputs: Optimal Pareto solutions(1)
Initialize parameters.(2)
Set system parameters:popsize, tmax, rmax(3)
Set genetic parameters:pcro, pmutImplement(4)
for run = 1: rmax(5)
Generate the initial population by greedy algorithm and scanning algorithm(6)
fori = 1: tmax(7)
fori = 1: popsize(8)
Evaluate the objective function(9)
End(10)
Carry outorder crossover and polynomial mutation for the population(11)
Generate the offspring populationQi(12)
Combine the parent populationPi and the offspring population Qi to generate a new population Ri(13)
Perform thenon-dominated sorting and crowding distance calculating for the Ri(14)
Select individuals from theRi to generate a new population Pi+1(15)
end(16)
end(17)
Obtain the optimal solutions(1)Nondominated Sorting and Crowding Distance Calculating Processes. The bi-objective values of each individual in the population can be first calculated. Then, the Pareto fronts are generated in accordance to the dominant relationships between the individuals [66]. Finally, the crowding distance is calculated for the individuals, which are selected to generate a new population for iteration based on the rank of Pareto fronts and crowding distances. The calculation process of the crowding distances of individuals is depicted in Figure 5.Figure 5
Calculation of crowding distances of individuals.In Figure5, pr and pl are two adjacent individuals of individual p with two objective values of O1 and O2. The crowding distance of individual p is expressed by Dp, where Dp=pr×O1−pl×O1+pr×O2−pl×O2. When selecting individuals to generate a new population, the individuals with the smaller rank in the nondominated solution set will be preferred. If the total number of the selected individuals is bigger than the size of the population, the individuals with the longer crowding distances in the last selected rank will be selected and added to the population [66]. After the maximum number of iterative operations, the solution in the optimal Pareto front (i.e., Rank 1) is the optimal Pareto solution of the bi-objective optimization problem.(2)Genetic Operation. Genetic operation mainly includes crossover operation and mutation operation. The application of the crossover operator can improve the search capability of the Im-NSGA-II. The convergence speed, stability, and global search capability can be improved by employing proper mutation operators [67]. The main processes of crossover operation and mutation operation are as follows.(1)
The order crossover operator is employed in the designed Im-NSGA-II as it outperforms many other operators [68]. First, two chromosomes such as Parent 1 and Parent 2 in the parent generation are selected to eliminate the values representing the logistics facility, and the customers encoding in the two chromosomes are both served by the same facility. Then, randomly select the starting and ending positions of the two preprocessed parent chromosomes. A section of the chromosome is retrieved in terms of two random and disparate points to generate proto-children. Finally, the rest positions in the proto-children will be filled with the genes not retrieved in Parent 1.(2)
The polynomial mutation plays an important role in multiobjective optimization problems [69]. Considering its advantages, the polynomial mutation is used in our proposed Im-NSGA-II. The offspring vi1,t+1 are obtained from the parent individuals vi1,t in equation (65), and the relative parameters are calculated in equation (66):(65)vi1,t+1=vi1,t+βi,(66)βi=2ui1/ηm+1−1,ui<0.5,1−21−ui1/1+ηm,ui≥0.5,whereηm represents the distribution index and ui is a random value, ui∈0,1.


### 4.2. Profit Allocation and Alliance Stability

#### 4.2.1. Shapley Value Method

Alliances consisting of logistics facilities contribute to cost savings and efficiency improvement. The cost savings referring to the profit of the alliance should be distributed to alliance members to motivate them to take profitable actions such as resource sharing [70]. Accordingly, the methods for profit allocation are important and necessary. The Shapley value method is utilized in this study. Here, A and As are defined as the grand alliance consisting of all logistics facilities in the 2E-CMDPDTW and a suballiance including some of them, respectively. The profit of each logistic facility m in the grand alliance A is represented by φmA,v and calculated in equation (67). In general, an alliance organizer responsible to send messages and coordinate members is deserved to gain extra profits. Thus, a given proportion σ∈0,1 of the total profit of the alliance is rewarded to the organizer. As for each suballiance, the total profit can be calculated in equation (68):(67)φmA,v=∑m∈As,As⊆AAs−1!A−As!A!×vAs−vAs−m,(68)vAs=1−σ×max∑m∈Ascm−cAs,0,where c({m}) represents the operating cost of a logistics facility m before taking part in the alliance As and c(As) depicts the total operating cost of suballiance As.


#### 4.2.2. Strictly Monotonic Path Method

The total profit of the alliance is significantly affected by the adding sequence of logistics facilities since the process of establishing the alliance is dynamic. The logistics facility with higher cost savings can be prioritized to join the alliance. The SMP is used to evaluate and select schemes of cooperative organization establishment [69]. In this study, the optimal sequence for logistics facilities to join the alliance is generated based on the SMP method. Let π represent a set of logistics facilities adding in the alliance with sequences and ηm,π,kdenote the decreased percentage of the operating cost of the logistics facility m when joining in the alliance. The decreased percentage can be calculated in equation (69):(69)ηm,π,k=φm∪πμ≤k,μ∈Aμ,vcm,πm≤k.As a consequence, the optimal sequence of adding in the alliance can be generated by conducting the SMP method. The selection processes are introduced as follows. According to the diagonal principle [57], the minimum decreased percentages on the diagonal values of the cost decreased percentage matrix are selected. Among these selected decreased percentages, the sequence with the highest decreased percentage is selected as the optimal sequence. If the lowest decreased percentages on the diagonal values are equal, comparisons of the second-lowest decreased percentages of different alliance sequences are then conducted. Among them, the optimal sequence is also determined by the highest decreased percentage. Finally, repeat the above processes until the optimal alliance sequence is obtained.


#### 4.2.3. Alliance Stability

Alliances are the basis for resource sharing and profit allocation in the 2E-CMDPDTW. The stability of the established alliance is the major focus. Despite the implemented profit allocation, the alliance stability is required to be evaluated. Based on the proximity between the schemes of different profit allocation methods and core center values, alliance stability can be guaranteed [68]. The core center value for each alliance is calculated in equation (70):(70)vA−vA−mvA×α+∑c∈A,c≠myc=vA−m,where v(A) represents the total profit of the grand alliance A, α is a variable that controls the profit changes before and after the logistics facility m joining the alliance A, and yc is the core center value of the logistics facility m joining the grand alliance A.


## 4.1. Hybrid Heuristics for the 2E-CMDPDTW Optimization

The NSGA-II proposed by Deb et al. [58] has been prominent in academic and practice domains in recent years. It is widely accepted that the NSGA-II can solve problems with multiple objectives effectively and efficiently. In the 2E-CMDPDTW, logistics facilities with different functionalities and amounts of customers with delivery and pickup demands are considered, which increases the computational complexity of the network. Therefore, the customer clustering algorithm and the greedy algorithm are introduced into the hybrid heuristic approach. The two-echelon logistics network is partitioned off as small regions using customer clustering, and the greedy algorithm aims to preliminarily generate the local optimal initial solution of each region [59]. Finally, the improved genetic operators are employed in the NSGA-II to find the global solution. The hybrid method is designed and implemented, as shown in Figure 2.Figure 2
Flowchart of the hybrid heuristic algorithm.The parameters in Figure2 can be explained as follows: k = {k1, k2} represents the number of clusters, in which k1 and k2 are the numbers of cluster centers, pop is the initial population, t and tmax represent the current and the maximum generations, respectively, and run and rmax denote the current and the maximum optimization runs, respectively.


### 4.1.1. Customer Clustering Algorithm

Customer clustering has been proved to be an effective method to simplify the large-scale and multidepot logistics network [7, 60]. Its basic process refers to reassign customers to logistics facilities closer to them. On one hand, to improve convergence speed and reduce computing efforts, the customer clustering algorithm is employed to generate initial feasible solutions. On the other hand, the customer clustering algorithm plays an important role in avoiding long-haul and cross transportation and reducing transportation costs. Algorithm 1 shows the procedure of the designed clustering algorithm.Algorithm 1: Procedure of the designed customer clustering algorithm.
Customer clustering algorithmInputs: Geographical coordinates of logistics facilities and geographical coordinates of customersOutputs: Clustering results(1)
Step 1: Set the number of LC and DCs k1 occupied in delivery tasks, and k2 indicates the number of LC and PCs handling product pickup demands.(2)
Step 2: Select elements in k1 and k2 as centers of corresponding clusters.(3)
For customers with delivery demands(4)
Step 3: Evaluate the distance between each delivery customer and each cluster center belonging to k1.(5)
For customers with pickup demands(6)
Step 4: Evaluate the distance between each pickup customer and each cluster center belonging to k2.(7)
If the same distances between a certain customer and several cluster centers existStep 5: Assign the customer to one of the cluster centers randomly.ElseStep 6: Assign the customer to the closest cluster center.End(8)
Step 7: Reassess the similarity of customers in each center, and update the clusters.(9)
Step 8: Repeat Steps 3-6 until every customer is assigned to the closest cluster center.(10)
Print the clustering results


### 4.1.2. Greedy Algorithm

The greedy algorithm referring to an approximate algorithm for optimization problems is usually employed to generate initial solutions [61]. The basic principle of the greedy algorithm is to guarantee the best choice in each search. With the greedy algorithm, this study can generate the local optimal solutions instead of the random ones, thereby improving the convergence speed. The major steps for implementing the greedy algorithm are introduced as follows.Step 1: initialize inputs according to the outputs of customer clustering.Cl1,Cdi2, and Cpj2 represent the sets of customers served by LC l, DC i, and PC j, respectively. Both customers with delivery and pickup demands can be served by the LC; thus, Rd1 and Rp1 represent the delivery and pickup route sets of vehicles from LC l. Similarly, Rdi2 and Rpj2 represent the route sets of vehicles from DC i and PC j.Step 2: define conditions for adding subroutes. A new subroute can be added into the vehicle route if the two customers in the subroute appeared less than twice in the vehicle route, either of the two customers does not appear in the vehicle route, and their time windows will not be violated.Step 3: generate the distance matrix of the customers and sort alternative subroutes. Calculate the distances between the customers with the same kind of requirements, and generate the distance matrix. Then, the lengths of alternative subroutes in the matrix are sorted from the shortest to the longest.Step 4: add the alternative subroutes according to the sorted results. If the alternative subroute fits the conditions above, it will be added to one of the route sets (e.g.,Rd1,Rp1, Rdi2, or Rpj2) according to the customer sets (i.e., Cl1,Cdi2, and Cpj2). In addition, if the customers in the subroutes have appeared in the vehicle route, the subroute will be directly connected to the existing route with the same customer. Repeat this step until the transportation quantity exceeds the maximum capacity of a vehicle.Step 5: generate closed vehicle routes. According to the clustering results, the customers in a certain vehicle route are served by the same facility so that a vehicle route is closed by adding the facility its customers belong to.


### 4.1.3. The Improved Nondominated Sorting Genetic Algorithm-II

The NSGA-II, a variant of the NSGA, contributes to balancing multiple objectives and generating the optimal decision on optimization problems [62, 63]. Considering the proposed model, this paper improves the method of generating the initial population and using the more appropriate genetic operators based on the original NSGA-II operating process. To improve the performance of global search and maintain the diversity of a population, the genetic operators including order crossover (OX) [45] and polynomial mutation [57] are exploited. The major procedures of implementing the Im-NSGA-II are presented in Figure 3.Figure 3
Major procedures of the Im-NSGA-II.In Figure3, the individuals in the initial population are first generated using the greedy algorithm and the scanning algorithm enabling the diversity of the population. Second, according to the nondomination sorting method [64], the individuals in population Pt are sorted. Each individual in Pt is assigned with a rank value representing the nondominated level. Selection, cross-over, and mutation operators are then employed to generate offspring in Qt. A new population Rt can be generated by combining the individuals in Pt and Qt. Third, a new population Pt+1 is generated by performing nondominated sorting and crowding calculation in Rt. The above process is repeated until the optimal Pareto solution is selected from all the Pareto solutions.Generally, it is hard to pursue all optimal objectives simultaneously in a bi-objective problem since they can be conflicting in nature [65]. According to these objectives, the solutions referring to Pareto solutions cannot be compared simply. Figure 4 shows the optimal Pareto solution of the bi-objective problem. Both individuals A and B are in the optimal Pareto front, and the others dominated by the solutions in the optimal Pareto front are in the feasible region. For the bi-objective problem, when objective 1 (e.g., number of vehicles) is deemed as the selection criterion, individual B is selected as the optimal Pareto solution because of fewer vehicles.Figure 4
Pareto solutions of the multiobjective optimization problem.The Im-NSGA-II is proposed to seek bi-objective optimal Pareto solutions in the 2E-CMDPDTW, as shown in Algorithm2. The processes of nondominated sorting, crowding distance calculating, and genetic operations are presented as follows.Algorithm 2: Procedure of the proposed Im-NSGA-II.
Im-NSGA-II algorithmInputs: Population sizepopsize, maximum number of generations tmax, maximum number of optimization runs rmax, crossover probability pcro, mutation probability pmutOutputs: Optimal Pareto solutions(1)
Initialize parameters.(2)
Set system parameters:popsize, tmax, rmax(3)
Set genetic parameters:pcro, pmutImplement(4)
for run = 1: rmax(5)
Generate the initial population by greedy algorithm and scanning algorithm(6)
fori = 1: tmax(7)
fori = 1: popsize(8)
Evaluate the objective function(9)
End(10)
Carry outorder crossover and polynomial mutation for the population(11)
Generate the offspring populationQi(12)
Combine the parent populationPi and the offspring population Qi to generate a new population Ri(13)
Perform thenon-dominated sorting and crowding distance calculating for the Ri(14)
Select individuals from theRi to generate a new population Pi+1(15)
end(16)
end(17)
Obtain the optimal solutions(1)Nondominated Sorting and Crowding Distance Calculating Processes. The bi-objective values of each individual in the population can be first calculated. Then, the Pareto fronts are generated in accordance to the dominant relationships between the individuals [66]. Finally, the crowding distance is calculated for the individuals, which are selected to generate a new population for iteration based on the rank of Pareto fronts and crowding distances. The calculation process of the crowding distances of individuals is depicted in Figure 5.Figure 5
Calculation of crowding distances of individuals.In Figure5, pr and pl are two adjacent individuals of individual p with two objective values of O1 and O2. The crowding distance of individual p is expressed by Dp, where Dp=pr×O1−pl×O1+pr×O2−pl×O2. When selecting individuals to generate a new population, the individuals with the smaller rank in the nondominated solution set will be preferred. If the total number of the selected individuals is bigger than the size of the population, the individuals with the longer crowding distances in the last selected rank will be selected and added to the population [66]. After the maximum number of iterative operations, the solution in the optimal Pareto front (i.e., Rank 1) is the optimal Pareto solution of the bi-objective optimization problem.(2)Genetic Operation. Genetic operation mainly includes crossover operation and mutation operation. The application of the crossover operator can improve the search capability of the Im-NSGA-II. The convergence speed, stability, and global search capability can be improved by employing proper mutation operators [67]. The main processes of crossover operation and mutation operation are as follows.(1)
The order crossover operator is employed in the designed Im-NSGA-II as it outperforms many other operators [68]. First, two chromosomes such as Parent 1 and Parent 2 in the parent generation are selected to eliminate the values representing the logistics facility, and the customers encoding in the two chromosomes are both served by the same facility. Then, randomly select the starting and ending positions of the two preprocessed parent chromosomes. A section of the chromosome is retrieved in terms of two random and disparate points to generate proto-children. Finally, the rest positions in the proto-children will be filled with the genes not retrieved in Parent 1.(2)
The polynomial mutation plays an important role in multiobjective optimization problems [69]. Considering its advantages, the polynomial mutation is used in our proposed Im-NSGA-II. The offspring vi1,t+1 are obtained from the parent individuals vi1,t in equation (65), and the relative parameters are calculated in equation (66):(65)vi1,t+1=vi1,t+βi,(66)βi=2ui1/ηm+1−1,ui<0.5,1−21−ui1/1+ηm,ui≥0.5,whereηm represents the distribution index and ui is a random value, ui∈0,1.


## 4.1.1. Customer Clustering Algorithm

Customer clustering has been proved to be an effective method to simplify the large-scale and multidepot logistics network [7, 60]. Its basic process refers to reassign customers to logistics facilities closer to them. On one hand, to improve convergence speed and reduce computing efforts, the customer clustering algorithm is employed to generate initial feasible solutions. On the other hand, the customer clustering algorithm plays an important role in avoiding long-haul and cross transportation and reducing transportation costs. Algorithm 1 shows the procedure of the designed clustering algorithm.Algorithm 1: Procedure of the designed customer clustering algorithm.
Customer clustering algorithmInputs: Geographical coordinates of logistics facilities and geographical coordinates of customersOutputs: Clustering results(1)
Step 1: Set the number of LC and DCs k1 occupied in delivery tasks, and k2 indicates the number of LC and PCs handling product pickup demands.(2)
Step 2: Select elements in k1 and k2 as centers of corresponding clusters.(3)
For customers with delivery demands(4)
Step 3: Evaluate the distance between each delivery customer and each cluster center belonging to k1.(5)
For customers with pickup demands(6)
Step 4: Evaluate the distance between each pickup customer and each cluster center belonging to k2.(7)
If the same distances between a certain customer and several cluster centers existStep 5: Assign the customer to one of the cluster centers randomly.ElseStep 6: Assign the customer to the closest cluster center.End(8)
Step 7: Reassess the similarity of customers in each center, and update the clusters.(9)
Step 8: Repeat Steps 3-6 until every customer is assigned to the closest cluster center.(10)
Print the clustering results


## 4.1.2. Greedy Algorithm

The greedy algorithm referring to an approximate algorithm for optimization problems is usually employed to generate initial solutions [61]. The basic principle of the greedy algorithm is to guarantee the best choice in each search. With the greedy algorithm, this study can generate the local optimal solutions instead of the random ones, thereby improving the convergence speed. The major steps for implementing the greedy algorithm are introduced as follows.Step 1: initialize inputs according to the outputs of customer clustering.Cl1,Cdi2, and Cpj2 represent the sets of customers served by LC l, DC i, and PC j, respectively. Both customers with delivery and pickup demands can be served by the LC; thus, Rd1 and Rp1 represent the delivery and pickup route sets of vehicles from LC l. Similarly, Rdi2 and Rpj2 represent the route sets of vehicles from DC i and PC j.Step 2: define conditions for adding subroutes. A new subroute can be added into the vehicle route if the two customers in the subroute appeared less than twice in the vehicle route, either of the two customers does not appear in the vehicle route, and their time windows will not be violated.Step 3: generate the distance matrix of the customers and sort alternative subroutes. Calculate the distances between the customers with the same kind of requirements, and generate the distance matrix. Then, the lengths of alternative subroutes in the matrix are sorted from the shortest to the longest.Step 4: add the alternative subroutes according to the sorted results. If the alternative subroute fits the conditions above, it will be added to one of the route sets (e.g.,Rd1,Rp1, Rdi2, or Rpj2) according to the customer sets (i.e., Cl1,Cdi2, and Cpj2). In addition, if the customers in the subroutes have appeared in the vehicle route, the subroute will be directly connected to the existing route with the same customer. Repeat this step until the transportation quantity exceeds the maximum capacity of a vehicle.Step 5: generate closed vehicle routes. According to the clustering results, the customers in a certain vehicle route are served by the same facility so that a vehicle route is closed by adding the facility its customers belong to.


## 4.1.3. The Improved Nondominated Sorting Genetic Algorithm-II

The NSGA-II, a variant of the NSGA, contributes to balancing multiple objectives and generating the optimal decision on optimization problems [62, 63]. Considering the proposed model, this paper improves the method of generating the initial population and using the more appropriate genetic operators based on the original NSGA-II operating process. To improve the performance of global search and maintain the diversity of a population, the genetic operators including order crossover (OX) [45] and polynomial mutation [57] are exploited. The major procedures of implementing the Im-NSGA-II are presented in Figure 3.Figure 3
Major procedures of the Im-NSGA-II.In Figure3, the individuals in the initial population are first generated using the greedy algorithm and the scanning algorithm enabling the diversity of the population. Second, according to the nondomination sorting method [64], the individuals in population Pt are sorted. Each individual in Pt is assigned with a rank value representing the nondominated level. Selection, cross-over, and mutation operators are then employed to generate offspring in Qt. A new population Rt can be generated by combining the individuals in Pt and Qt. Third, a new population Pt+1 is generated by performing nondominated sorting and crowding calculation in Rt. The above process is repeated until the optimal Pareto solution is selected from all the Pareto solutions.Generally, it is hard to pursue all optimal objectives simultaneously in a bi-objective problem since they can be conflicting in nature [65]. According to these objectives, the solutions referring to Pareto solutions cannot be compared simply. Figure 4 shows the optimal Pareto solution of the bi-objective problem. Both individuals A and B are in the optimal Pareto front, and the others dominated by the solutions in the optimal Pareto front are in the feasible region. For the bi-objective problem, when objective 1 (e.g., number of vehicles) is deemed as the selection criterion, individual B is selected as the optimal Pareto solution because of fewer vehicles.Figure 4
Pareto solutions of the multiobjective optimization problem.The Im-NSGA-II is proposed to seek bi-objective optimal Pareto solutions in the 2E-CMDPDTW, as shown in Algorithm2. The processes of nondominated sorting, crowding distance calculating, and genetic operations are presented as follows.Algorithm 2: Procedure of the proposed Im-NSGA-II.
Im-NSGA-II algorithmInputs: Population sizepopsize, maximum number of generations tmax, maximum number of optimization runs rmax, crossover probability pcro, mutation probability pmutOutputs: Optimal Pareto solutions(1)
Initialize parameters.(2)
Set system parameters:popsize, tmax, rmax(3)
Set genetic parameters:pcro, pmutImplement(4)
for run = 1: rmax(5)
Generate the initial population by greedy algorithm and scanning algorithm(6)
fori = 1: tmax(7)
fori = 1: popsize(8)
Evaluate the objective function(9)
End(10)
Carry outorder crossover and polynomial mutation for the population(11)
Generate the offspring populationQi(12)
Combine the parent populationPi and the offspring population Qi to generate a new population Ri(13)
Perform thenon-dominated sorting and crowding distance calculating for the Ri(14)
Select individuals from theRi to generate a new population Pi+1(15)
end(16)
end(17)
Obtain the optimal solutions(1)Nondominated Sorting and Crowding Distance Calculating Processes. The bi-objective values of each individual in the population can be first calculated. Then, the Pareto fronts are generated in accordance to the dominant relationships between the individuals [66]. Finally, the crowding distance is calculated for the individuals, which are selected to generate a new population for iteration based on the rank of Pareto fronts and crowding distances. The calculation process of the crowding distances of individuals is depicted in Figure 5.Figure 5
Calculation of crowding distances of individuals.In Figure5, pr and pl are two adjacent individuals of individual p with two objective values of O1 and O2. The crowding distance of individual p is expressed by Dp, where Dp=pr×O1−pl×O1+pr×O2−pl×O2. When selecting individuals to generate a new population, the individuals with the smaller rank in the nondominated solution set will be preferred. If the total number of the selected individuals is bigger than the size of the population, the individuals with the longer crowding distances in the last selected rank will be selected and added to the population [66]. After the maximum number of iterative operations, the solution in the optimal Pareto front (i.e., Rank 1) is the optimal Pareto solution of the bi-objective optimization problem.(2)Genetic Operation. Genetic operation mainly includes crossover operation and mutation operation. The application of the crossover operator can improve the search capability of the Im-NSGA-II. The convergence speed, stability, and global search capability can be improved by employing proper mutation operators [67]. The main processes of crossover operation and mutation operation are as follows.(1)
The order crossover operator is employed in the designed Im-NSGA-II as it outperforms many other operators [68]. First, two chromosomes such as Parent 1 and Parent 2 in the parent generation are selected to eliminate the values representing the logistics facility, and the customers encoding in the two chromosomes are both served by the same facility. Then, randomly select the starting and ending positions of the two preprocessed parent chromosomes. A section of the chromosome is retrieved in terms of two random and disparate points to generate proto-children. Finally, the rest positions in the proto-children will be filled with the genes not retrieved in Parent 1.(2)
The polynomial mutation plays an important role in multiobjective optimization problems [69]. Considering its advantages, the polynomial mutation is used in our proposed Im-NSGA-II. The offspring vi1,t+1 are obtained from the parent individuals vi1,t in equation (65), and the relative parameters are calculated in equation (66):(65)vi1,t+1=vi1,t+βi,(66)βi=2ui1/ηm+1−1,ui<0.5,1−21−ui1/1+ηm,ui≥0.5,whereηm represents the distribution index and ui is a random value, ui∈0,1.


## 4.2. Profit Allocation and Alliance Stability

### 4.2.1. Shapley Value Method

Alliances consisting of logistics facilities contribute to cost savings and efficiency improvement. The cost savings referring to the profit of the alliance should be distributed to alliance members to motivate them to take profitable actions such as resource sharing [70]. Accordingly, the methods for profit allocation are important and necessary. The Shapley value method is utilized in this study. Here, A and As are defined as the grand alliance consisting of all logistics facilities in the 2E-CMDPDTW and a suballiance including some of them, respectively. The profit of each logistic facility m in the grand alliance A is represented by φmA,v and calculated in equation (67). In general, an alliance organizer responsible to send messages and coordinate members is deserved to gain extra profits. Thus, a given proportion σ∈0,1 of the total profit of the alliance is rewarded to the organizer. As for each suballiance, the total profit can be calculated in equation (68):(67)φmA,v=∑m∈As,As⊆AAs−1!A−As!A!×vAs−vAs−m,(68)vAs=1−σ×max∑m∈Ascm−cAs,0,where c({m}) represents the operating cost of a logistics facility m before taking part in the alliance As and c(As) depicts the total operating cost of suballiance As.


### 4.2.2. Strictly Monotonic Path Method

The total profit of the alliance is significantly affected by the adding sequence of logistics facilities since the process of establishing the alliance is dynamic. The logistics facility with higher cost savings can be prioritized to join the alliance. The SMP is used to evaluate and select schemes of cooperative organization establishment [69]. In this study, the optimal sequence for logistics facilities to join the alliance is generated based on the SMP method. Let π represent a set of logistics facilities adding in the alliance with sequences and ηm,π,kdenote the decreased percentage of the operating cost of the logistics facility m when joining in the alliance. The decreased percentage can be calculated in equation (69):(69)ηm,π,k=φm∪πμ≤k,μ∈Aμ,vcm,πm≤k.As a consequence, the optimal sequence of adding in the alliance can be generated by conducting the SMP method. The selection processes are introduced as follows. According to the diagonal principle [57], the minimum decreased percentages on the diagonal values of the cost decreased percentage matrix are selected. Among these selected decreased percentages, the sequence with the highest decreased percentage is selected as the optimal sequence. If the lowest decreased percentages on the diagonal values are equal, comparisons of the second-lowest decreased percentages of different alliance sequences are then conducted. Among them, the optimal sequence is also determined by the highest decreased percentage. Finally, repeat the above processes until the optimal alliance sequence is obtained.


### 4.2.3. Alliance Stability

Alliances are the basis for resource sharing and profit allocation in the 2E-CMDPDTW. The stability of the established alliance is the major focus. Despite the implemented profit allocation, the alliance stability is required to be evaluated. Based on the proximity between the schemes of different profit allocation methods and core center values, alliance stability can be guaranteed [68]. The core center value for each alliance is calculated in equation (70):(70)vA−vA−mvA×α+∑c∈A,c≠myc=vA−m,where v(A) represents the total profit of the grand alliance A, α is a variable that controls the profit changes before and after the logistics facility m joining the alliance A, and yc is the core center value of the logistics facility m joining the grand alliance A.


## 4.2.1. Shapley Value Method

Alliances consisting of logistics facilities contribute to cost savings and efficiency improvement. The cost savings referring to the profit of the alliance should be distributed to alliance members to motivate them to take profitable actions such as resource sharing [70]. Accordingly, the methods for profit allocation are important and necessary. The Shapley value method is utilized in this study. Here, A and As are defined as the grand alliance consisting of all logistics facilities in the 2E-CMDPDTW and a suballiance including some of them, respectively. The profit of each logistic facility m in the grand alliance A is represented by φmA,v and calculated in equation (67). In general, an alliance organizer responsible to send messages and coordinate members is deserved to gain extra profits. Thus, a given proportion σ∈0,1 of the total profit of the alliance is rewarded to the organizer. As for each suballiance, the total profit can be calculated in equation (68):(67)φmA,v=∑m∈As,As⊆AAs−1!A−As!A!×vAs−vAs−m,(68)vAs=1−σ×max∑m∈Ascm−cAs,0,where c({m}) represents the operating cost of a logistics facility m before taking part in the alliance As and c(As) depicts the total operating cost of suballiance As.


## 4.2.2. Strictly Monotonic Path Method

The total profit of the alliance is significantly affected by the adding sequence of logistics facilities since the process of establishing the alliance is dynamic. The logistics facility with higher cost savings can be prioritized to join the alliance. The SMP is used to evaluate and select schemes of cooperative organization establishment [69]. In this study, the optimal sequence for logistics facilities to join the alliance is generated based on the SMP method. Let π represent a set of logistics facilities adding in the alliance with sequences and ηm,π,kdenote the decreased percentage of the operating cost of the logistics facility m when joining in the alliance. The decreased percentage can be calculated in equation (69):(69)ηm,π,k=φm∪πμ≤k,μ∈Aμ,vcm,πm≤k.As a consequence, the optimal sequence of adding in the alliance can be generated by conducting the SMP method. The selection processes are introduced as follows. According to the diagonal principle [57], the minimum decreased percentages on the diagonal values of the cost decreased percentage matrix are selected. Among these selected decreased percentages, the sequence with the highest decreased percentage is selected as the optimal sequence. If the lowest decreased percentages on the diagonal values are equal, comparisons of the second-lowest decreased percentages of different alliance sequences are then conducted. Among them, the optimal sequence is also determined by the highest decreased percentage. Finally, repeat the above processes until the optimal alliance sequence is obtained.


## 4.2.3. Alliance Stability

Alliances are the basis for resource sharing and profit allocation in the 2E-CMDPDTW. The stability of the established alliance is the major focus. Despite the implemented profit allocation, the alliance stability is required to be evaluated. Based on the proximity between the schemes of different profit allocation methods and core center values, alliance stability can be guaranteed [68]. The core center value for each alliance is calculated in equation (70):(70)vA−vA−mvA×α+∑c∈A,c≠myc=vA−m,where v(A) represents the total profit of the grand alliance A, α is a variable that controls the profit changes before and after the logistics facility m joining the alliance A, and yc is the core center value of the logistics facility m joining the grand alliance A.


## 5. Empirical Analyses

### 5.1. Algorithm Comparison and Analysis

The comparison between the Im-NSGA-II and the other two famous multiobjective optimization algorithms is conducted to evaluate the performance of the designed algorithm. In this section, the NSGA-II and the multiobjective genetic algorithm (MOGA) are employed since they are of interest to researchers in multiobjective optimization fields [71, 72]. The logistics networks in the modified “Solomon datasets” [73] and the “MDVRPTW instances2” [69] are selected as benchmarks. The parameters in the datasets are shown in Table 5.Table 5
Parameters in the datasets.
InstanceSolomon datasetsMDVRPTW datasetsNumber of customersNumber of delivery customersNumber of delivery customersNumber of depots1—pr_0148242442C101—50252533C202—50252534R101—50252535R204—50252536RC104—50252537—pr_0772363668—pr_0296484849C102—1005050410C201—1005050411R102—1005050412RC105—1005050413RC203—1005050414—pr_031447272415C1_2_1—200100100516R1_2_2—200100100517RC1_2_4—200100100618—pr_09216108108619—pr_05240120120420—pr_062881441444The parameters for optimization algorithm implementation are listed as follows,popsize = 100, gmax = 250, pcro = 0.9, pmut = 0.1, and runmax = 10. The mentioned three algorithms are implemented ten times for each instance in the datasets, and the optimal results are compared in Table 6.Table 6
Comparative results of three algorithms.
InstanceIm-NSGA-IINSGA-IIMOGACostVehicleTimeCostVehicleTimeCostVehicleTime11,3838601,4069671,4428672685664712769719670368156676566375766541,37212651,48212741,47612795836466874477901476686966390767598177671,62910681,72410791,689118282,06412772,19512862,136128491,36811831,46510961,325109810970585992599989792111,6747841,7988971,776795121,93815882,01516982,0331895131,3806861,4429971,4367106142,687171082,793191122,822181171511,9542514212,1032715112,262271531610,1392514510,6572615310,263291671714,2042314914,1092415314,77924157183,053241643,266271723,44725171193,335242083,472262213,59126222204,057292934,285323154,29033312Average3,314141083,423151183,45615119t test−4.283E + 00−8.942E + 00−4.574E + 00−9.401E + 00p value3.626E − 042.004E − 081.840E − 048.849E − 09In Table6, the designed Im-NSGA-II outperforms the NSGA-II and the MOGA in most instances. In Instances 1 to 7, the Im-NSGA-II can obtain the minimum operating costs. Despite that the computation time by the Im-NSGA-II is longer than that by the NSGA-II and the MOGA, the minimum operating cost and number of vehicles are guaranteed in Instance 3. In Instances 8 to 14, excepting for Instances 9, the Im-NSGA-II can obtain the minimum values of the operating cost and the computation time, and the number of vehicles is not larger than that by the other algorithms. In Instance 9, both the operating cost and the number of vehicles by the Im-NSGA-II are higher than those in MOGA, while the computation time is minimum. In Instance 17, the Im-NSGA-II obtains a higher operating cost than NSGA-II, but both the number of vehicles and the computation time are the minimum.


### 5.2. Data Source

A real-world 2E-LPDN in Chongqing China is utilized to illustrate and evaluate the proposed model and solution method in the 2E-CMDPDTW. In Figure6, a 2E-LPDN consists of one LC, two DCs (i.e., DC1 and DC2), two PCs (i.e., PC1 and PC2), 90 delivery customers, and 92 pickup customers. The logistics center labeled LC is represented by a red node with the letter L. The customers served by the LC are divided into two parts: the delivery customers are represented by red nodes with a dot symbol, and the pickup customers are represented by red nodes with a star symbol. DC1 and DC2 along with their customers are represented by blue and green nodes with a dot symbol, respectively. Similarly, nodes in yellow and white with a star symbol indicate PC1 and PC2 along with their customers, respectively. In addition, the four facilities are labeled DC1, DC2, PC1, and PC2, respectively, and distinguished from their customers. The characteristics of the LC, DCs, PCs, and their customers are presented in Table 7.Figure 6
Geographic layouts of the LC, DCs, PCs, and customers.Table 7
The characteristics of the LC, DCs, PCs, and customers.
FacilitiesSymbolNumber of customersCustomer typeServed customersCustomer symbolLC28Delivery customersD1∼D18Pickup customersP1∼P10DC132Delivery customersD19∼D50DC240Delivery customersD51∼D90PC138Pickup customersP11∼P48PC244Pickup customersP49∼P90


### 5.3. Parameter Setting and Optimization Results

Using the proposed model and solution method, the logistics network is redesigned with the optimal route in terms of operating cost savings and transportation resource reduction. First of all, the customer clustering algorithm is implemented to reassign the 182 customers, and sequentially, the greed algorithm is conducted to generate initial feasible solutions. The results of customer clustering are shown in Figure7. Finally, the Im-NSGA-II is run to pursue optimal Pareto solutions. The parameters in this case study are based on literature [7, 25, 46, 69] and practical investigations. The parameters referring to the operation of a semitrailer truck and a vehicle are set as Ks = 1700, Kv = 200, fs = 6.6, fv = 3.4, Ms = 15000, Mv = 1500, ps = 3.4, and pv = 3.2. The parameters related to the operation and collaboration of logistics facilities are set as F1 = 1200, F2 = 800, F3 = 750, F4 = 850, F5 = 900, c1 = 400, c2 = 550, c3 = 450, c4 = 600, c5 = 600, we = 0.2, wd = 0.3, and B = 52. Besides, the parameters involved in the Im-NSGA-II are popsize = 100, tmax = 500, pcro = 0.9, pmut = 0.1, and runmax = 2. The established alliances are organized by the LC in the logistics network. Therefore, at most 20 types of alliances are available. According to the proposed model and hybrid heuristic algorithm, we calculate the comparison results of optimization schemes for all forms of alliances shown in Table 8. Figure 8 illustrates the comparison of the optimization schemes.Figure 7
Results of customer clustering.Table 8
Comparative results of the costs, the number of trucks, and the number of vehicles.
AllianceInitial networkOptimized networkSavings and reductionCostVehicleTruckCostVehicleTruckCostVehicleTruck{LC}6,4925—5,7935—6990—{DC1}5,4145—4,9514—4631—{DC2}7,0576—6,5795—4781—{PC1}5,7235—5,0564—6671—{PC2}7,2695—6,6905—5790—{LC, DC1}11,9051018,827913,07810{LC, DC2}13,54811110,3591013,18910{LC, PC1}12,2141018,648913,56710{LC, PC2}13,76110110,070913,69110{LC, DC1, DC2}18,96216213,1111425,85120{LC, DC1, PC1}17,62815211,6821315,94721{LC, DC1, PC2}19,17415214,1041315,07021{LC, PC1, DC2}19,27116213,2141416,05721{LC, PC1, PC2}19,48315213,9441415,54011{LC, PC2, DC2}20,81716214,6361416,18121{LC, DC1, DC2, PC1}26,23121317,6041828,62731{LC, DC1, DC2, PC2}24,89721316,6621828,23531{LC, PC1, PC2, DC1}26,54021317,7671828,77331{LC, PC1, PC2, DC2}25,46321316,5381928,92521{LC, DC1, DC2, PC1, PC2}31,95426420,16123211,79332Figure 8
Comparison of the costs, number of trucks, and number of vehicles of alliances.In Table8 and Figure 8, the total operating cost, the number of trucks, and the number of vehicles are reduced compared with the initial 2E-LPDN. For example, the decreased percentages of the total cost, the number of vehicles, and the number of semitrailer trucks can reach 36.9%, 11.5%, and 50%, respectively. In addition, an increasing trend of the reduction of the cost savings and the transportation resource is observed as the alliance members increase. For example, the decreased percentage of the total cost increases from 30.86% to 32.89% as well as the number of vehicles from 12.5% to 14.29%, when PC1 has been a member of the alliance {LC, DC1, DC2}. Therefore, establishing the grand alliance must be most effective to redesign and optimize the logistics network in this case. Table 9 shows the comparative results of the initial and the optimized 2E-CMDPDTW. Figure 9 illustrates the comparative results of the network before and after establishing the grand alliance.Table 9
Comparative results of the initial and the optimized 2E-CMDPDTW.
CaseTransportation cost ($)Pickup and delivery cost ($)Penalty cost ($)Maintenance cost ($)Collaborative cost ($)Total operating cost ($)Number of trucksNumber of vehiclesInitial 2E-CMDPDTW79201821039201904031954426Optimized 2E-CMDPDTW5859938110811240260020161223Figure 9
Comparative results of the network before and after establishing the grand alliance.In Table9 and Figure 9, the total cost of the grand alliance includes the transportation cost, pickup and delivery cost, penalty cost, maintenance cost, and collaborative cost. Through establishing a grand alliance, the resource sharing in the 2E-CMDPDTW can be realized to effectively reduce the total cost and the number of vehicles in the 2E-LPDN. For example, the total cost of the initial 2E-CMDPDTW is 31,954, while the total cost of the optimized 2E-CMDPDTW is 20,161 with a decrease value of 11,793. Moreover, the number of trucks and vehicles is 4 and 26 in the initial network, respectively. By contrast, with a value of 2 and 23, the numbers of trucks and numbers of vehicles are reduced by 2 and 3 in the optimized network, respectively.


### 5.4. Profit Allocation and Sequential Coalition Selection

It is widely accepted that economic benefits are the focus of the majority of business strategies. Similarly, potential profits earned in alliances contribute to decisions making on participation and alliance member sequence, and vice versa. According to the results in the case study, the logistics network can earn an overall profit in terms of the grand alliance. Therefore, the establishment and the stability of alliances are depending on effective profit allocations to alliance members. According to practitioners, the reward percentage is given asσ=0.12, which means that 12% of the total profit is rewarded to the alliance provider. The results of the profit allocation in the case are presented, as shown in Table 10.Table 10
Profit of each logistics facility in the 2E-CMDPDTW.
AllianceInitial costOptimized costvAsφmA,v{LC}6,4925,793699(699, 0, 0, 0, 0){DC1}5,4144,951463(0, 463, 0, 0, 0){DC2}7,0576,579478(0, 0, 478, 0, 0){PC1}5,7235,056667(0, 0, 0, 667, 0){PC2}7,2696,690579(0, 0, 0, 0, 579){LC, DC1}11,9058,8273,078(1657, 1421, 0, 0, 0){LC, DC2}13,54810,3593,189(1705, 0, 1484, 0, 0){LC, PC1}12,2148,6483,567(1799, 0, 0, 1768, 0){LC, PC2}13,76110,0703,691(1706, 0, 0, 0, 1985){LC, DC1, DC2}18,96213,1115,851(2389, 1699, 1762, 0, 0){LC, DC1, PC1}17,62811,6825,947(2091, 1755, 0, 2101, 0){LC, DC1, PC2}19,17414,1045,070(2030, 1088, 0, 0, 1952){LC, DC2, PC1}19,27113,2146,057(2139, 0, 1817, 2101, 0){LC, PC1, PC2}19,48313,9445,540(2166, 0, 0, 1428, 1946){LC, DC2 PC2}20,81714,6366,181(2210, 0, 1886, 0, 2085){LC, DC1, DC2, PC1}26,23117,6048,627(2442, 1922, 1998, 2265, 0){LC, DC1, DC2, PC2}24,89716,6628,235(2394, 1853, 1902, 0, 2086){LC, DC1, PC1, PC2}26,54017,7678,773(2465, 1899, 0, 2276, 2133){LC, DC2, PC1, PC2}25,46316,5388,925(2471, 0, 2002, 2306, 2146){LC, DC1, DC2, PC1, PC2}31,95420,16111,793(2684, 2069, 2165, 2326, 2549)Compared to the network without optimization, the LC, DC1, DC2, PC1, and PC2 can save costs with values of 699, 463, 478, 677, and 579 in the optimized network. After optimizing the network and establishing the grand alliance, the total cost can be reduced with a value of 11,793 referring to the profit. Accordingly, these logistics facilities can acquire allocated profits with values of 2,684, 2,069, 2165, 2326, and 2,549, respectively. Figure10 depicts the processes along with the decreased percentages of cost savings of the logistics facilities participating in the alliance. For transparency, two examples are presented: (1) the processes complying with the SMP are represented by solid lines. (2) The others are depicted by dash lines. Starting from each logistic facility, the best coalition sequences satisfying SMP are presented, as shown in Table 11. The optimal sequence {PC1, LC, DC1, PC2, DC2} is obtained based on the diagonal rule, and the decreased proportions of costs for each facility in it are illustrated in Figure 11.Figure 10
Cost reduction percentages when new members join the alliance.Table 11
Best coalition sequences satisfying SMP from each logistics facility.
π1 = {DC1, LC, DC2, PC1, PC2}π2 = {DC2, LC, PC1, DC1, PC2}PlayeriDC1LCDC2PC1PC2PlayeriDC2LCPC1DC1PC2η (i, π, 1)8.6%η (i, π, 1)6.8%η (i, π, 2)26.3%25.5%η (i, π, 2)21.0%26.3%η (i, π, 3)31.4%36.8%25.0%η (i, π, 3)25.8%33.0%36.7%η (i, π, 4)35.5%37.6%28.3%39.6%η (i, π, 4)28.3%37.6%39.6%35.5%η (i, π, 5)41.3%38.2%30.7%40.6%35.1%η (i, π, 5)30.7%41.3%40.7%38.2%35.1%π3 = {PC1, LC, DC1, PC2, DC2}π4 = {PC2, LC, DC2, PC1, DC1}PlayeriPC1LCDC1PC2DC2PlayeriPC2LCDC2PC1DC1η (i, π, 1)11.7%η (i, π, 1)8.0%η (i, π, 2)30.9%27.7%η (i, π, 2)27.3%26.3%η (i, π, 3)36.7%32.2%32.4%η (i, π, 3)28.7%34.0%26.7%η (i, π, 4)39.8%38.0%35.1%29.3%η (i, π, 4)29.5%38.1%28.4%40.3%η (i, π, 5)40.6%41.3%38.2%36.0%30.7%η (i, π, 5)35.1%41.3%30.7%40.6%38.2%π5 = {LC, PC1, DC1, PC2, DC2}PlayeriLCPC1DC1PC2DC2η (i, π, 1)10.8%η (i, π, 2)27.7%30.9%η (i, π, 3)32.2%36.7%32.4%η (i, π, 4)38.0%39.8%35.1%29.3%η (i, π, 5)41.3%40.6%38.2%35.1%30.7%Figure 11
Cost reduction percentages of the grand alliance.In Figure11, as the first member of the alliance, PC1 acquires 11.7% cost savings. The costs of PC1 and LC are decreased by 30.9% and 27.7% when the second member LC participates. Accordingly, when the grand alliance is generated, 40.6%, 41.3%, 38.2%, 36.0%, and 30.7% costs are reduced for PC1, LC, DC1, PC2, and DC2, respectively.


### 5.5. Alliance Stability

To demonstrate the effectiveness of the employed profit allocation method, the comparisons are conducted between the Shapley value method and other methods such as the equal profit method (EPM) [74], the game quadratic programming (GQP) [57], and the minimum costs-remaining savings (MCRS) [75]. As mentioned in Section 4.2.3, the comparative results of the proximity between profit allocation schemes generated by different methods and the center values are presented, as shown in Table 12 and demonstrated in Figure 12.Table 12
Comparative results between different schemes and center value.
MethodProfit allocation schemesCenterDistanceShapley(2684, 2069, 2165, 2326, 2549)(2702, 2316, 2029, 2312, 2434)305EPM(2913, 2132, 2277, 2118, 2353)429MCRS(3113, 2032, 2077, 2218, 2353)517GQP(2706, 2059, 2077, 2428, 2547)308Figure 12
Comparison of distances between different methods and center value.In Table12 and Figure 12, the value of the core center is (2702, 2316, 2029, 2312, 2434). Comparative results between the core center value and schemes using different methods suggest that the profit allocation scheme via the MCRS method is the farthest from the core center. The profit of DC1 in MCRS is less than that in the other profit allocation schemes, and thus, DC1 is more likely to refuse to join the alliance to further improve its profit. Similarly, both the profit of LC in GQP and the profit of PC1 in EPM can be improved. The distances between profit allocation schemes using the EPM, MCRS, and GQP methods are 429, 517, and 308, which are not the closest. Therefore, the stability and development of the alliance will be threatened. It is obvious that the closest scheme is generated by the Shapley value method. Therefore, the logistics facilities will participate in the grand alliance and construct a long-term partnership when the profit allocation scheme is set as (2,684, 2,069, 2,165, 2,326, 2,549).


## 5.1. Algorithm Comparison and Analysis

The comparison between the Im-NSGA-II and the other two famous multiobjective optimization algorithms is conducted to evaluate the performance of the designed algorithm. In this section, the NSGA-II and the multiobjective genetic algorithm (MOGA) are employed since they are of interest to researchers in multiobjective optimization fields [71, 72]. The logistics networks in the modified “Solomon datasets” [73] and the “MDVRPTW instances2” [69] are selected as benchmarks. The parameters in the datasets are shown in Table 5.Table 5
Parameters in the datasets.
InstanceSolomon datasetsMDVRPTW datasetsNumber of customersNumber of delivery customersNumber of delivery customersNumber of depots1—pr_0148242442C101—50252533C202—50252534R101—50252535R204—50252536RC104—50252537—pr_0772363668—pr_0296484849C102—1005050410C201—1005050411R102—1005050412RC105—1005050413RC203—1005050414—pr_031447272415C1_2_1—200100100516R1_2_2—200100100517RC1_2_4—200100100618—pr_09216108108619—pr_05240120120420—pr_062881441444The parameters for optimization algorithm implementation are listed as follows,popsize = 100, gmax = 250, pcro = 0.9, pmut = 0.1, and runmax = 10. The mentioned three algorithms are implemented ten times for each instance in the datasets, and the optimal results are compared in Table 6.Table 6
Comparative results of three algorithms.
InstanceIm-NSGA-IINSGA-IIMOGACostVehicleTimeCostVehicleTimeCostVehicleTime11,3838601,4069671,4428672685664712769719670368156676566375766541,37212651,48212741,47612795836466874477901476686966390767598177671,62910681,72410791,689118282,06412772,19512862,136128491,36811831,46510961,325109810970585992599989792111,6747841,7988971,776795121,93815882,01516982,0331895131,3806861,4429971,4367106142,687171082,793191122,822181171511,9542514212,1032715112,262271531610,1392514510,6572615310,263291671714,2042314914,1092415314,77924157183,053241643,266271723,44725171193,335242083,472262213,59126222204,057292934,285323154,29033312Average3,314141083,423151183,45615119t test−4.283E + 00−8.942E + 00−4.574E + 00−9.401E + 00p value3.626E − 042.004E − 081.840E − 048.849E − 09In Table6, the designed Im-NSGA-II outperforms the NSGA-II and the MOGA in most instances. In Instances 1 to 7, the Im-NSGA-II can obtain the minimum operating costs. Despite that the computation time by the Im-NSGA-II is longer than that by the NSGA-II and the MOGA, the minimum operating cost and number of vehicles are guaranteed in Instance 3. In Instances 8 to 14, excepting for Instances 9, the Im-NSGA-II can obtain the minimum values of the operating cost and the computation time, and the number of vehicles is not larger than that by the other algorithms. In Instance 9, both the operating cost and the number of vehicles by the Im-NSGA-II are higher than those in MOGA, while the computation time is minimum. In Instance 17, the Im-NSGA-II obtains a higher operating cost than NSGA-II, but both the number of vehicles and the computation time are the minimum.


## 5.2. Data Source

A real-world 2E-LPDN in Chongqing China is utilized to illustrate and evaluate the proposed model and solution method in the 2E-CMDPDTW. In Figure6, a 2E-LPDN consists of one LC, two DCs (i.e., DC1 and DC2), two PCs (i.e., PC1 and PC2), 90 delivery customers, and 92 pickup customers. The logistics center labeled LC is represented by a red node with the letter L. The customers served by the LC are divided into two parts: the delivery customers are represented by red nodes with a dot symbol, and the pickup customers are represented by red nodes with a star symbol. DC1 and DC2 along with their customers are represented by blue and green nodes with a dot symbol, respectively. Similarly, nodes in yellow and white with a star symbol indicate PC1 and PC2 along with their customers, respectively. In addition, the four facilities are labeled DC1, DC2, PC1, and PC2, respectively, and distinguished from their customers. The characteristics of the LC, DCs, PCs, and their customers are presented in Table 7.Figure 6
Geographic layouts of the LC, DCs, PCs, and customers.Table 7
The characteristics of the LC, DCs, PCs, and customers.
FacilitiesSymbolNumber of customersCustomer typeServed customersCustomer symbolLC28Delivery customersD1∼D18Pickup customersP1∼P10DC132Delivery customersD19∼D50DC240Delivery customersD51∼D90PC138Pickup customersP11∼P48PC244Pickup customersP49∼P90


## 5.3. Parameter Setting and Optimization Results

Using the proposed model and solution method, the logistics network is redesigned with the optimal route in terms of operating cost savings and transportation resource reduction. First of all, the customer clustering algorithm is implemented to reassign the 182 customers, and sequentially, the greed algorithm is conducted to generate initial feasible solutions. The results of customer clustering are shown in Figure7. Finally, the Im-NSGA-II is run to pursue optimal Pareto solutions. The parameters in this case study are based on literature [7, 25, 46, 69] and practical investigations. The parameters referring to the operation of a semitrailer truck and a vehicle are set as Ks = 1700, Kv = 200, fs = 6.6, fv = 3.4, Ms = 15000, Mv = 1500, ps = 3.4, and pv = 3.2. The parameters related to the operation and collaboration of logistics facilities are set as F1 = 1200, F2 = 800, F3 = 750, F4 = 850, F5 = 900, c1 = 400, c2 = 550, c3 = 450, c4 = 600, c5 = 600, we = 0.2, wd = 0.3, and B = 52. Besides, the parameters involved in the Im-NSGA-II are popsize = 100, tmax = 500, pcro = 0.9, pmut = 0.1, and runmax = 2. The established alliances are organized by the LC in the logistics network. Therefore, at most 20 types of alliances are available. According to the proposed model and hybrid heuristic algorithm, we calculate the comparison results of optimization schemes for all forms of alliances shown in Table 8. Figure 8 illustrates the comparison of the optimization schemes.Figure 7
Results of customer clustering.Table 8
Comparative results of the costs, the number of trucks, and the number of vehicles.
AllianceInitial networkOptimized networkSavings and reductionCostVehicleTruckCostVehicleTruckCostVehicleTruck{LC}6,4925—5,7935—6990—{DC1}5,4145—4,9514—4631—{DC2}7,0576—6,5795—4781—{PC1}5,7235—5,0564—6671—{PC2}7,2695—6,6905—5790—{LC, DC1}11,9051018,827913,07810{LC, DC2}13,54811110,3591013,18910{LC, PC1}12,2141018,648913,56710{LC, PC2}13,76110110,070913,69110{LC, DC1, DC2}18,96216213,1111425,85120{LC, DC1, PC1}17,62815211,6821315,94721{LC, DC1, PC2}19,17415214,1041315,07021{LC, PC1, DC2}19,27116213,2141416,05721{LC, PC1, PC2}19,48315213,9441415,54011{LC, PC2, DC2}20,81716214,6361416,18121{LC, DC1, DC2, PC1}26,23121317,6041828,62731{LC, DC1, DC2, PC2}24,89721316,6621828,23531{LC, PC1, PC2, DC1}26,54021317,7671828,77331{LC, PC1, PC2, DC2}25,46321316,5381928,92521{LC, DC1, DC2, PC1, PC2}31,95426420,16123211,79332Figure 8
Comparison of the costs, number of trucks, and number of vehicles of alliances.In Table8 and Figure 8, the total operating cost, the number of trucks, and the number of vehicles are reduced compared with the initial 2E-LPDN. For example, the decreased percentages of the total cost, the number of vehicles, and the number of semitrailer trucks can reach 36.9%, 11.5%, and 50%, respectively. In addition, an increasing trend of the reduction of the cost savings and the transportation resource is observed as the alliance members increase. For example, the decreased percentage of the total cost increases from 30.86% to 32.89% as well as the number of vehicles from 12.5% to 14.29%, when PC1 has been a member of the alliance {LC, DC1, DC2}. Therefore, establishing the grand alliance must be most effective to redesign and optimize the logistics network in this case. Table 9 shows the comparative results of the initial and the optimized 2E-CMDPDTW. Figure 9 illustrates the comparative results of the network before and after establishing the grand alliance.Table 9
Comparative results of the initial and the optimized 2E-CMDPDTW.
CaseTransportation cost ($)Pickup and delivery cost ($)Penalty cost ($)Maintenance cost ($)Collaborative cost ($)Total operating cost ($)Number of trucksNumber of vehiclesInitial 2E-CMDPDTW79201821039201904031954426Optimized 2E-CMDPDTW5859938110811240260020161223Figure 9
Comparative results of the network before and after establishing the grand alliance.In Table9 and Figure 9, the total cost of the grand alliance includes the transportation cost, pickup and delivery cost, penalty cost, maintenance cost, and collaborative cost. Through establishing a grand alliance, the resource sharing in the 2E-CMDPDTW can be realized to effectively reduce the total cost and the number of vehicles in the 2E-LPDN. For example, the total cost of the initial 2E-CMDPDTW is 31,954, while the total cost of the optimized 2E-CMDPDTW is 20,161 with a decrease value of 11,793. Moreover, the number of trucks and vehicles is 4 and 26 in the initial network, respectively. By contrast, with a value of 2 and 23, the numbers of trucks and numbers of vehicles are reduced by 2 and 3 in the optimized network, respectively.


## 5.4. Profit Allocation and Sequential Coalition Selection

It is widely accepted that economic benefits are the focus of the majority of business strategies. Similarly, potential profits earned in alliances contribute to decisions making on participation and alliance member sequence, and vice versa. According to the results in the case study, the logistics network can earn an overall profit in terms of the grand alliance. Therefore, the establishment and the stability of alliances are depending on effective profit allocations to alliance members. According to practitioners, the reward percentage is given asσ=0.12, which means that 12% of the total profit is rewarded to the alliance provider. The results of the profit allocation in the case are presented, as shown in Table 10.Table 10
Profit of each logistics facility in the 2E-CMDPDTW.
AllianceInitial costOptimized costvAsφmA,v{LC}6,4925,793699(699, 0, 0, 0, 0){DC1}5,4144,951463(0, 463, 0, 0, 0){DC2}7,0576,579478(0, 0, 478, 0, 0){PC1}5,7235,056667(0, 0, 0, 667, 0){PC2}7,2696,690579(0, 0, 0, 0, 579){LC, DC1}11,9058,8273,078(1657, 1421, 0, 0, 0){LC, DC2}13,54810,3593,189(1705, 0, 1484, 0, 0){LC, PC1}12,2148,6483,567(1799, 0, 0, 1768, 0){LC, PC2}13,76110,0703,691(1706, 0, 0, 0, 1985){LC, DC1, DC2}18,96213,1115,851(2389, 1699, 1762, 0, 0){LC, DC1, PC1}17,62811,6825,947(2091, 1755, 0, 2101, 0){LC, DC1, PC2}19,17414,1045,070(2030, 1088, 0, 0, 1952){LC, DC2, PC1}19,27113,2146,057(2139, 0, 1817, 2101, 0){LC, PC1, PC2}19,48313,9445,540(2166, 0, 0, 1428, 1946){LC, DC2 PC2}20,81714,6366,181(2210, 0, 1886, 0, 2085){LC, DC1, DC2, PC1}26,23117,6048,627(2442, 1922, 1998, 2265, 0){LC, DC1, DC2, PC2}24,89716,6628,235(2394, 1853, 1902, 0, 2086){LC, DC1, PC1, PC2}26,54017,7678,773(2465, 1899, 0, 2276, 2133){LC, DC2, PC1, PC2}25,46316,5388,925(2471, 0, 2002, 2306, 2146){LC, DC1, DC2, PC1, PC2}31,95420,16111,793(2684, 2069, 2165, 2326, 2549)Compared to the network without optimization, the LC, DC1, DC2, PC1, and PC2 can save costs with values of 699, 463, 478, 677, and 579 in the optimized network. After optimizing the network and establishing the grand alliance, the total cost can be reduced with a value of 11,793 referring to the profit. Accordingly, these logistics facilities can acquire allocated profits with values of 2,684, 2,069, 2165, 2326, and 2,549, respectively. Figure10 depicts the processes along with the decreased percentages of cost savings of the logistics facilities participating in the alliance. For transparency, two examples are presented: (1) the processes complying with the SMP are represented by solid lines. (2) The others are depicted by dash lines. Starting from each logistic facility, the best coalition sequences satisfying SMP are presented, as shown in Table 11. The optimal sequence {PC1, LC, DC1, PC2, DC2} is obtained based on the diagonal rule, and the decreased proportions of costs for each facility in it are illustrated in Figure 11.Figure 10
Cost reduction percentages when new members join the alliance.Table 11
Best coalition sequences satisfying SMP from each logistics facility.
π1 = {DC1, LC, DC2, PC1, PC2}π2 = {DC2, LC, PC1, DC1, PC2}PlayeriDC1LCDC2PC1PC2PlayeriDC2LCPC1DC1PC2η (i, π, 1)8.6%η (i, π, 1)6.8%η (i, π, 2)26.3%25.5%η (i, π, 2)21.0%26.3%η (i, π, 3)31.4%36.8%25.0%η (i, π, 3)25.8%33.0%36.7%η (i, π, 4)35.5%37.6%28.3%39.6%η (i, π, 4)28.3%37.6%39.6%35.5%η (i, π, 5)41.3%38.2%30.7%40.6%35.1%η (i, π, 5)30.7%41.3%40.7%38.2%35.1%π3 = {PC1, LC, DC1, PC2, DC2}π4 = {PC2, LC, DC2, PC1, DC1}PlayeriPC1LCDC1PC2DC2PlayeriPC2LCDC2PC1DC1η (i, π, 1)11.7%η (i, π, 1)8.0%η (i, π, 2)30.9%27.7%η (i, π, 2)27.3%26.3%η (i, π, 3)36.7%32.2%32.4%η (i, π, 3)28.7%34.0%26.7%η (i, π, 4)39.8%38.0%35.1%29.3%η (i, π, 4)29.5%38.1%28.4%40.3%η (i, π, 5)40.6%41.3%38.2%36.0%30.7%η (i, π, 5)35.1%41.3%30.7%40.6%38.2%π5 = {LC, PC1, DC1, PC2, DC2}PlayeriLCPC1DC1PC2DC2η (i, π, 1)10.8%η (i, π, 2)27.7%30.9%η (i, π, 3)32.2%36.7%32.4%η (i, π, 4)38.0%39.8%35.1%29.3%η (i, π, 5)41.3%40.6%38.2%35.1%30.7%Figure 11
Cost reduction percentages of the grand alliance.In Figure11, as the first member of the alliance, PC1 acquires 11.7% cost savings. The costs of PC1 and LC are decreased by 30.9% and 27.7% when the second member LC participates. Accordingly, when the grand alliance is generated, 40.6%, 41.3%, 38.2%, 36.0%, and 30.7% costs are reduced for PC1, LC, DC1, PC2, and DC2, respectively.


## 5.5. Alliance Stability

To demonstrate the effectiveness of the employed profit allocation method, the comparisons are conducted between the Shapley value method and other methods such as the equal profit method (EPM) [74], the game quadratic programming (GQP) [57], and the minimum costs-remaining savings (MCRS) [75]. As mentioned in Section 4.2.3, the comparative results of the proximity between profit allocation schemes generated by different methods and the center values are presented, as shown in Table 12 and demonstrated in Figure 12.Table 12
Comparative results between different schemes and center value.
MethodProfit allocation schemesCenterDistanceShapley(2684, 2069, 2165, 2326, 2549)(2702, 2316, 2029, 2312, 2434)305EPM(2913, 2132, 2277, 2118, 2353)429MCRS(3113, 2032, 2077, 2218, 2353)517GQP(2706, 2059, 2077, 2428, 2547)308Figure 12
Comparison of distances between different methods and center value.In Table12 and Figure 12, the value of the core center is (2702, 2316, 2029, 2312, 2434). Comparative results between the core center value and schemes using different methods suggest that the profit allocation scheme via the MCRS method is the farthest from the core center. The profit of DC1 in MCRS is less than that in the other profit allocation schemes, and thus, DC1 is more likely to refuse to join the alliance to further improve its profit. Similarly, both the profit of LC in GQP and the profit of PC1 in EPM can be improved. The distances between profit allocation schemes using the EPM, MCRS, and GQP methods are 429, 517, and 308, which are not the closest. Therefore, the stability and development of the alliance will be threatened. It is obvious that the closest scheme is generated by the Shapley value method. Therefore, the logistics facilities will participate in the grand alliance and construct a long-term partnership when the profit allocation scheme is set as (2,684, 2,069, 2,165, 2,326, 2,549).


## 6. Conclusions

The collaboration between logistics facilities is an effective management strategy to promote resource sharing, provide considerable profits for logistics enterprises, and improve the performance of urban pickup and delivery networks. To improve efficiency and reduce the total operating cost, the sharing of customer information and transportation resources is considered. First, a bi-objective integer programming model is developed to minimize the total operating cost and the number of vehicles in the 2E-CMDPDTW. Second, a hybrid method including the customer clustering algorithm, the greedy algorithm, and the Im-NSGA-II is developed to solve the mathematical model. Then, the Shapley value method is used for profit allocation in the alliance established in the optimized network, and the optimal sequence of joining the alliance can be generated using the SMP. Finally, the alliance stability is discussed by comparing different profit allocation methods {Fernie, 2007 #3217}.An empirical case in Chongqing, China, is used to test the effectiveness of the proposed model and solution method for solving the 2E-CMDPDTW. Compared to the noncollaborative network, the operating cost and the number of vehicles are reduced by $11,793 and 3, respectively. The solution speed and quality of the proposed Im-NSGA-II are better than the existing multiobjective optimization algorithms, such as NSGA-II and MOGA. A reasonable profit allocation scheme for the alliance is provided by the Shapley value method. Taking into account the profit-seeking nature of the logistics facilities, the alliance organizer is more likely to select the profit allocation scheme provided by the Shapley value method rather than the GQP, MCRS, and EPM methods. The profits of $2326, $2684, $2069, $2549, and $2165 are allocated to the five logistics facilities in the optimal scheme {PC1, LC, DC1, PC2, DC2} generated by the SMP method. Establishing a grand alliance consisting of all logistics facilities in the network maximizes cost savings and transportation resource utilization. Besides, the reasonable profit distribution schemes are conducive to the stable development of the alliance, and they are determined by the appropriate distribution methods.From the operational perspective, resource sharing has the potential to share operating risks and contributes to cost savings for logistics enterprises. Another key reason for promoting resource sharing is to achieve on-time deliveries via customer information sharing. As for urban transportation governance, resource sharing plays a significant role in optimizing logistics delivery and pickup networks and mitigating urban traffic congestions. Therefore, resource sharing among logistics facilities provides a promising way for logistics enterprises to reduce costs and improve efficiency and for the government to improve the regional logistics service capacity and build more efficient urban transportation networks.Despite an effective approach to solving 2E-CMDPDTW is proposed in this study, the following concerns need to be addressed in future research. (1) In terms of logistics network design, customers with dynamic demands or multiple service time windows can be modeled in a two-echelon collaborative multideport logistics network. (2) For the 2E-MDLPDN optimization, more noneconomic factors (e.g., environmental factors) can be considered. (3) In the optimization solution approach, exact algorithms can be combined with the heuristic algorithm to obtain high-quality optimal solutions and a faster convergence speed. (4) In the aspect of logistics collaboration, the contributions of different alliance members to the alliance stability can be further explored on the basis of different profit allocation methods.



---

*Source: 6619539-2021-10-11.xml*