# Multidepot Recycling Vehicle Routing Problem with Resource Sharing and Time Window Assignment

Authors: Yong Wang; Xiuwen Wang; Xiangyang Guan; Jinjun Tang
Journal: Journal of Advanced Transportation
 (2021)

---

## Abstract

This study aims to provide tactical and operational decisions in multidepot recycling logistics networks with consideration of resource sharing (RS) and time window assignment (TWA) strategies. The RS strategy contributes to efficient resource allocation and utilization among recycling centers (RCs). The TWA strategy involves assigning time windows to customers to enhance the operational efficiency of logistics networks. A biobjective mathematical model is established to minimize the total operating cost and number of vehicles for solving the multidepot recycling vehicle routing problem with RS and TWA (MRVRPRSTWA). A hybrid heuristic algorithm including 3Dk-means clustering algorithm and nondominated sorting genetic algorithm- (NSGA-) II (NSGA-II) is designed. The 3D k-means clustering algorithm groups customers into clusters on the basis of their spatial and temporal distances to reduce the computational complexity in optimizing the multidepot logistics networks. In comparison with NSGA algorithm, the NSGA-II algorithm incorporates an elitist strategy, which can improve the computational speed and robustness. In this study, the performance of the NSGA-II algorithm is compared with the other two algorithms. Results show that the proposed algorithm is superior in solving MRVRPRSTWA. The proposed model and algorithm are applied to an empirical case study in Chongqing City, China, to test their applicability in real logistics operations. Four different scenarios regarding whether the RS and TWA strategies are included or not are developed to test the efficacy of the proposed methods. The results indicate that the RS and TWA strategies can optimize the recycling services and resource allocation and utilization and enhance the operational efficiency, thus promoting the sustainable development of the logistics industry.


---

## Body

## 1. Introduction

Sustainability is a constant topic in the logistics industry and research [1, 2]. In reverse logistics, recycling waste products can save resources and facilitate sustainable development [3]. Recycling vehicle routing problem (RVRP), one part of reverse logistics, often appears in logistics operations. For example, empty beer bottles are recycled at convenience stores. An increasing number of enterprises have begun to consider recycling products to improve the utilization of resources due to economic, legal, and social concerns [4, 5]. Therefore, the growing demand for recycling services will stimulate the development of reverse logistics. The design, planning, and operation of sustainable reverse logistics systems pose a challenge to logistics organizations [6].The low transportation resource utilization and irrational transport operations of logistics facilities still appear and need to be addressed [7, 8]. Resource sharing (RS) strategies have been proposed to optimize resource utilization and reduce operating cost in multidepot logistics networks [9]. Customers are reassigned to appropriate logistics facilities through these strategies to avoid unreasonable routing via customer information sharing [10]. Transportation resources are shared within and between logistics facilities to improve transportation resource utilization. Specifically, vehicles can be used to serve customers of a single logistics facility or customers of multiple logistics facilities many times when the time windows of customers allow [11]. In this study, the RS strategy is developed in the multidepot RVRP with time window (MRVRPTW) to optimize the resource utilization.In traditional vehicle routing problems with time window (VRPTW), the service time windows for customers are obtained by the logistics service providers after negotiation with customers or are determined by customers [12, 13]. Logistics service providers often recycle cargoes from customers at the same time every day [14]. In real life, the time windows for most customers are viewed as soft time windows, that is, they can be violated by paying certain penalties [15, 16]. Most logistics service providers try to recycle cargoes within the predetermined time windows in logistics operations to reduce the penalty cost [17]. However, the vehicles cannot provide recycling services for customers in the expected time windows due to long transportation distance and unreasonable service time windows for logistics facilities [18, 19]. Therefore, a set of candidate time windows are assigned to customers in this study, and recycling vehicle routes are designed on the basis of time window assignment (TWA) to optimize the operational efficiency and minimize the operating cost in multidepot RVRP with RS and TWA (MRVRPRSTWA).In this study, the RS and TWA strategies are adopted to improve resource utilization, optimize resource allocation, and enhance the operational efficiency of multidepot recycling logistics networks. A bi-objective mathematical model is presented to obtain the minimum operating cost and the number of vehicles (TNV) in a multidepot logistics network with RS and TWA. A hybrid heuristic algorithm, including a 3Dk-means clustering algorithm and nondominated sorting genetic algorithm- (NSGA-) II (NSGA-II), is devised to optimize the recycling routes. The 3D k-means clustering algorithm groups customers on the basis of their spatial and temporal distances to reduce the computational complexity. NSGA-II is used to assign time windows to customers and find the optimal solutions in the MRVRPRSTWA strategy. A real-world case study is performed to verify the applicability of MRVRPRSTWA in logistics operations. The results confirm the capability of the RS and TWA strategies in optimizing multidepot logistics networks.The rest of this paper is organized as follows. In Section2, relevant literature is reviewed. In Section 3, MRVRPRSTWA is described in detail. In Section 4, a multiobjective mathematical model is formulated to minimize the total operating cost and TNV. In Section 5, a hybrid heuristic algorithm is introduced to solve MRVRPRSTWA. In Section 6, a real-world case is presented to verify the applicability of the proposed method. In Section 7, the conclusions are drawn, and future research directions are proposed.


## 2. Literature Review

With the increasing attention on environmental protection, recycling and related problems have received great concern. The RVRP with time window (RVRPTW) as the main form of reverse logistics has been widely studied to reduce the operating cost in recent years [20]. The multidepot RVRPTW (MRVRPTW) is developed to coordinate the recycling services among logistics facilities [21], and the RS strategy can be considered to optimize the resource utilization in MRVRPTW (MRVRPRSTW) [22]. Candidate time windows are assigned to customers via the TWA strategy to improve operational efficiency compared with the traditional time window in MRVRPRSTW [23]. MRVRPRSTW with TWA (MRVRPRSTWA) investigates the RS and TWA strategies to optimize the resource allocation and operational efficiency in multidepot logistics networks.Increasing interest in recycled products and materials has led to a growing concern in recycling logistics [24, 25]. El korchi and Millet [26] investigated a framework, which used to generate and assess the reverse logistics channel structure, and proposed an environment-friendly structure with high economic benefits to optimize recycling services. Chan et al. [27] proposed a framework of reverse logistics on the basis of the value of reversed products, cost issue, and legislation perspective to optimize reverse logistics. Bai and Sarkis [28] introduced a flexible framework for reverse logistics and a performance measurement of third-party reverse logistics service provider model on the basis of a novel neighborhood rough set method to set the foundation for future research of reverse logistics. Roghanian and Pazhoheshfar [29] presented a probabilistic mixed liner programming model and genetic algorithm to minimize the total cost in reverse logistics. Shaik and Abdul-Kader [30] proposed a multi-criteria performance measurement model considering product lifecycle stages, strategies, capabilities, processes, perspectives, and measures to verify and enhance the enterprise’s performance in reverse logistics.The RS strategy, as a major issue, has drawn increasing attention on the enhanced resource utilization and reduced operating cost [31–33]. An et al. [34] constructed a network DEA model to calculate the cost of logistics network before and after RS and show the performance of the RS strategy. Sun et al. [35] analyzed and evaluated the effect of the RS strategy on alleviating or solving the traffic and environmental problems in megacities from five perspectives of resource, environment, convenience, economy, and governance. Quintero-Araujo et al. [36] investigated different collaborative scenarios on the basis of the RS strategy to obtain a lower logistics cost. The RS strategy is developed to minimize TNV and the total operating cost for improving the efficiency and sustainability of logistics networks [37]. Molina et al. [38] established a mathematical model that aims at maximizing the number of served customers and minimizing the total cost in the case of RS. Wang et al. [39] studied the RS strategy in multiple service periods and multidepot to enhance resource efficiency and refine the resource configuration.Clustering algorithms are always used to group customers before determining vehicle routes in MDVRPTW to simplify the complexity of calculation [40, 41]. Min et al. [42] proposed the maximum-minimum distance clustering method to cluster customers in split-delivery VRP for the better performance of the algorithm. Ge et al. [43] added the service radius and load expansion factors to the clustering algorithm to avoid the vehicle overloading in MDVRPTW. Fan et al. [44] introduced a clustering algorithm based on the temporal-spatial distance to reduce the computational complexity and enhance the quality of initial solution in MDVRPTW. Liu et al. [45] investigated the clustering algorithm that groups customers in terms of the minimum distance and maximum demand to improve the efficiency of solving MDVRPTW.Several algorithms are studied to improve the quality of solution in MRVRPTW [46–48]. Aras et al. [21] established two mixed-integer linear programming models and proposed a Tabu-search-based heuristic algorithm to solve MRVRPTW. Kim et al. [49] developed a two-stage solution procedure to minimize the total distance and cost in MRVRPTW. Liu and He [50] proposed a clustering-based multiple ant colony system algorithm to minimize the total cost for enhancing the route compactness in MRVRPTW. Ramos and Oliveira [51] developed a mathematical model to minimize the variable costs and attain equity, and designed a heuristic algorithm to address MRVRPTW. Liao [4] established the genetic mixed-integer nonlinear programming model to maximize total profit and designed a hybrid genetic algorithm to solve MRVRPTW. Govindan et al. [52] proposed a mathematical model and a fuzzy solution method to minimize the total cost in MRVRPTW.TWA vehicle routing problem (TWAVRP) differs from the traditional VRPTW, where the service time windows are determined by customers. In TWAVRP, candidate time windows are assigned to customers to improve the operational efficiency of logistics networks [53, 54]. Subramanyam et al. [55] introduced a scenario decomposition algorithm to solve TWAVRP for minimizing the expected routing cost of logistics facility. Neves-Moreira et al. [56] presented a two-stage stochastic optimization problem of TWAVRP and addressed it by using a fix-and-optimize-based metaheuristic. Jalilvand et al. [57] developed a stochastic model and a progressive hedging algorithm to minimize the routing cost. The results show that the efficiency of the proposed method has good performance in addressing TWAVRP. A mathematics model is proposed to enhance the operational efficiency and reduce the total cost, and a hybrid heuristic algorithm is used for TWAVRP [58].The aforementioned literature covers many aspects of MRVRPRSTWA but has the following limitations: (1) The RS strategy is rarely considered between and within logistics facilities. (2) The TWA strategy is insufficiently discussed in multidepot logistics network optimization. (3) An efficient mathematical model for optimizing TNV in the presence of RS strategy is lacking. (4) The traditional evolutionary algorithms have limited applicability when solving the MRVRPRSTWA.The contributions of this paper are as follows: (1) The RS strategy is incorporated in the MRVRPRSTWA to maximize resource utilization. (2) The TWA strategy is developed in the MRVRPRSTWA to improve the operational efficiency. (3) A bi-objective integer programming model is proposed to minimize the total operating cost and TNV in multidepot logistics networks with shared resources. (4) A hybrid heuristic algorithm, including 3Dk-means clustering and NSGA-II, is developed to solve MRVRPRSTWA.


## 3. Problem Statement

The RS and TWA strategies are proposed in MRVRPRSTWA to optimize the recycling service of the multidepot logistics network. The multidepot logistics network is composed of multiple RCs and customers. Long-distance and crisscross trips often exist in the initial nonoptimal logistics network (Figure1). The transportation resources are shared among logistics facilities, and customers with irrational time windows are reassigned new time windows to adjust optimized recycling routes through the TWA strategy. The comparison of the multidepot logistics network before and after optimization is shown in Figures 1 and 2.Figure 1
Logistics network without the RS and TWA strategies.Figure 2
Logistics network with the RS and TWA strategies.In Figure1, each RC operates independently in the nonoptimal logistics network. Several customers are served by the farther RC rather than the closer one. For example, customer 21 is closer to RC3, while it is served by RC1, resulting in long-distance and crisscross transportation. In addition, RCs can only provide on-time service for two-thirds of customers, and the time windows for the remaining third of customers cannot be met. Relevant measures must be taken to enhance the operational efficiency and optimize the configuration of logistics networks. In Figure 2, centralized transportation is performed among RCs by a fleet of semitrailer trucks and is used to merge and transship goods between RCs, considering changes in customer demand assignment. Resources are shared among RCs, and customers are reassigned to the corresponding RC based on geographical locations and time windows. Irrational transportation phenomenon is reduced, and vehicles are shared among RCs. For example, due to the difference of service time between different routes, V2 is shared between RC1 and RC2, which first serve the customers of RC2 and then the customers of RC1, and thus the vehicle utilization rate is improved and the logistics operating cost is reduced. The violations of time windows are decreased when the TWA strategy is adopted. The RS and TWA strategies call for improving the utilization of transportation resources and operational efficiency and optimizing multidepot logistics networks.This study assumes that the centralized transportation cost (TC) among logistics facilities is $22 per unit time, and the recycling cost from logistics facilities to customers is $20 per unit time. The penalty cost for earliness and delay service is $20 per unit time. The assignment cost (AC) from the initial time window to assigned time (AT) window is $5 per unit time. The maintenance cost (MC) of each semitrailer truck is $300, and that of each vehicle is $200. The AT, TNV, the number of semitrailer trucks (TNS), and the total cost are compared and listed in Table1.Table 1
Comparison of before and after RS and TWA.
CaseTC ($)Recycling cost ($)Penalty cost ($)ATAC ($)TNVTNSMC ($)Total cost ($)Nonoptimal logistics network—17601080——9—18004640Optimal logistics with RS and TWA9901400400562804111004170In Table1, the total cost of reverse logistics network is decreased from $4640 to $4170. In addition, TNV is reduced from nine to five when sharing resources among logistics facilities. When the TWA strategy was applied, the penalty cost decreases by 63%. The RS and TWA strategies can improve the resource utilization and operational efficiency and optimize multidepot logistics networks.


## 4. Model Formulation for MRVRPTWA

### 4.1. Related Definitions and Variables

In this section, the related notations and definitions on sets, parameters, and variables are shown in Table2.Table 2
Notations and definitions in the MRVRPRSTWA.
Set definitionPSet of recycling centers (RCs),P = {p|p = 1, 2, 3…, m} and m is the total number of RCsCSet of all recycling customers,C = {c|c = 1, 2, 3…, h} and h is the total number of customersCfSet of customers with fixed time windowsCaSet of customers with assigned time windowsVSet of vehicles for recycling,V = {v|v = 1, 2, 3…, b} and b is the total number of vehiclesSSet of semitrailer trucks used between RCs,S = {s|s = 1, 2, 3…, w} and w is the total number of vehiclesOvMaximum times for executing recycling routes of vehiclev in a working period, v∈VParametersqiRecycling demand quantity of customeri, i∈CtqmnTransportation quantity from logistics facilitym to n, m, n∈PQsMaximum capacity of semitrailer trucks, s∈SQvMaximum capacity of vehiclev, v∈VQmMaximum capacity of logistics facilitym, m∈PUsUsage cost of semitrailer trucks, s∈S (unit: dollar/km)UvUsage cost of vehiclev, v∈V (unit: dollar/km)DijDistance from logistics facility or customeri to j, i, j∈P∪C, i ≠ j (unit: km)DmnDistance from logistics facilitym to logistics facility n, m, n∈P (unit: km)WNumber of working days in one planning periodBNumber of planning periods in one yearMsAnnual MC of semitrailer trucks, s∈SMvAnnual MC of vehiclev, v∈VFLmFixed cost of logistics facilitym, m∈P[ei,ti]Fixed time window of customeri, i∈Cf[e’ i,l’ i]Expected time window of customeri, i∈Ca[αi,βi]Time window assigned to customeri, i∈Ca[Em,Lm]Service time window of logistics facilitym, m∈PμePenalty cost for early arrival per unit timeμdPenalty cost for delayed arrival per unit timeμaCost coefficient of the customer’s time window from the expected time window to the assigned time window per unit timedtk vmDeparture time of thekth route of vehicle V from logistics facility m, v ∈ V, m ∈ P, k ∈ Ovatk ivArrival time of thekth route of vehicle v at node i, v ∈ V, i ∈ P ∪ C, k ∈ OvtijvTravel time of vehiclev between entities i and j, i, j ∈ P ∪ C, v ∈ V|Ns|Total number of logistic facilities served by semitrailer trucks, s∈SNvTotal number of customers served by vehiclev, v∈VBNBig numberDecision variablesxijvIf vehiclev travels from logistics facility or customer i to j, then xijv = 1; otherwise, xijv = 0, i, j ∈ P∪C, v ∈ Vzk vIf vehiclev has the kth route, then zk v = 1; otherwise, zk v = 0, m ∈ P,v ∈ V, and k ∈ OvωivmIf vehiclev departs from logistics facility m to served customer i, ωivm = 1; otherwise, ωivm = 0, i ∈ C, v ∈ V, and m ∈ PymnsIf semitrailer trucks transports between logistics facility m and n, ymnk = 1; otherwise, ymns = 0, m, n ∈ P, and s ∈ SτimnIf customeri’s logistic facility changes from logistics facility m to n after optimization, τimn = 1; otherwise, τimn = 0, i ∈ C, and m, n ∈ PgmsIf semitrailer trucks departs from logistics facility m, gms = 1; otherwise, gms = 0, s ∈ S, and m ∈ P


### 4.2. Mathematical Model

A bi-objective optimization model for MRVRPRSTWA is formulated to minimize the total operating cost in equation (1) and minimize TNV in equation (2).(1)Min TC=TC1+TC2+TC3,(2)Min MNV=∑v∈Vmin∑i∈P∑j∈Cxijv,1.TC1 contains three components in equation (3): ∑m∈P∑n∈P∑s∈SUs×Dmn×W×ymns represents the centralized TCs by semitrailer trucks among logistics facilities. maxm∈P∑n∈Ptqmn/Qs×Ms/B is the MC of semitrailer trucks. ∑m∈PFCm represents the fixed cost of logistics facilities.(3)TC1=∑m∈P∑n∈P∑s∈SUs×Dmn×W×ymns+maxm∈P∑n∈PtqmnQs×MsB+∑m∈PFLm.TC2 contains two components in equation (4): ∑i∈P∪C∑j∈P∪C∑v∈VUv×Dij×W×xijv is the TC by vehicles from logistics facilities to customers. ∑v∈Vmin∑i∈P∑j∈Cxijv,1×Mv/B represents the MC of vehicles.(4)TC2=∑i∈P∪C∑j∈P∪C∑v∈VUv×Dij×W×xijv+∑v∈Vmin∑i∈P∑j∈Cxijv,1×MvB.TC3 contains five components in equation (5): ∑i∈P∪Cf∑j∈Cf∑v∈V∑k∈Ovxijv×μe×maxei−ativk,0 and ∑i∈P∪Ca∑j∈Ca∑v∈V∑k∈Ovxijv×μe×maxαi−ativk,0 represent the penalty cost caused by vehicles that provide earlier service than time windows of customers. ∑i∈P∪Ca∑j∈Ca∑v∈V∑k∈Ovxijv×μd×maxativk−li,0 and ∑i∈P∪Ca∑j∈Ca∑v∈V∑k∈Ovxijv×μd×maxativk−βi,0 represent the penalty cost caused by vehicles that provide delayed service than time windows of customers. ∑i∈Caμa×minei′−αi,βi−li′ denotes the cost of TWA.(5)TC3=∑i∈P∪Cf∑j∈Cf∑v∈V∑k∈Ovxijv×μe×maxei−ativk,0+∑i∈P∪Cf∑j∈Cf∑v∈V∑k∈Ovxijv×μd×maxativk−li,0+∑i∈P∪Ca∑j∈Ca∑v∈V∑k∈Ovxijv×μe×maxαi−ativk,0+∑i∈P∪Ca∑j∈Ca∑v∈V∑k∈Ovxijv×μd×maxativk−βi,0+∑i∈Caμa×minei′−αi,βi−li′.Subject to(6)∑m∈Pymns=1,∀n∈P,s∈S,m≠n,(7)∑n∈Pymns−∑n∈Pynfs=0,∀m,f∈P,m≠f,s∈S,(8)∑m∈Pynms=1,∀n∈P,s∈S,m≠n,(9)∑m∈P∑n∈Pymns×gps≤Ns−1,∀p∈P,s∈S,p≠m,p≠n,(10)tqmn=∑i∈Cτimn×qi,∀m,n∈P,(11)∑m∈P∑n∈Ptqmn×ymns≤Qs,∀s∈S,(12)∑i∈Pxijv=1,∀j∈C,v∈V,(13)∑i∈C∑j∈Cxijv=1,∀v∈V,(14)∑j∈Pxijv−∑j∈Pxjfv=0,∀i,f∈P∪C,i≠f,v∈V,(15)∑i∈Pxjiv=1,∀j∈C,v∈V,(16)∑i∈C∑j∈Cxijv≤Nv−1,∀v∈V,(17)∑j∈Cqj×xijv×zvk≤Qv,∀i∈P∪C,v∈V,k∈Ov,(18)∑j∈C∑v∈Vqj×ωijvm≤Qm,∀i∈P∪C,m∈P,(19)dtvmk+tmiv−BN1−xmiv≤ativk,∀m∈P,i∈C,v∈V,k∈Ov,(20)dtvmk+tmiv+BN1−xmiv≥ativk,∀m∈P,i∈C,v∈V,k∈Ov,(21)ativk+tijv−BN1−xijv≤atjvk,∀i∈C,j∈C∪P,v∈V,k∈Ov,(22)ativk+tijv+BN1−xijv≥atjvk,∀i∈C,j∈C∪P,v∈V,k∈Ov,(23)Em×xmiv×zvk≤dtvmk≤Lm×xmiv×zvk,∀m∈P,i∈C,v∈V,k∈Ov,(24)Em×xmiv×zvk≤atmvk≤Lm×xmiv×zvk,∀m∈P,i∈C,v∈V,k∈Ov,(25)atmvk+tmnv−M1−zvk+1≤dtvnk+1,∀m,n∈P,v∈V,k∈Ov,(26)xijv=0,1,∀i,j∈P∪C,v∈V,(27)zvk=0,1,∀v∈V,k∈Ov,(28)ωivm=0,1,∀i∈C,v∈V,m∈P,(29)ymns=0,1,∀m,n∈P,∀s∈S,(30)τimn=0,1,∀i∈C,∀m,n∈P,m≠n,(31)gms=0,1,∀m∈P,s∈S.Constraint (6) ensures that each semitrailer truck originally departs from logistics facilities. Constraint (7) is the flow conservation on each logistics facility. Constraint (8) ensures that each semitrailer truck finally returns to logistics facilities. Constraint (9) is used to avoid subtours for semitrailer trucks. Constraint (10) calculates the transportation quantity among logistics facilities. Constraint (11) stipulates that each semitrailer truck has enough capacity to satisfy the transportation quantity of the assigned logistics facilities. Constraint (12) ensures that each vehicle departs from logistics facilities only. Constraint (13) ensures that each customer can be served once by one vehicle. Constraint (14) is the flow conservation on each customer. Constraint (15) ensures that each vehicle finally returns to logistics facilities. Constraint (16) is used to avoid subtours for vehicles. Constraint (17) represents that each vehicle has enough capacity to satisfy the total demand of assigned customers. Constraint (18) ensures that the total demand of customers served by a logistics facility should be beyond its capacity. Constraints (19) and (20) guarantee the continuous departure time of vehicles at logistics facilities. Constraints (21) and (22) guarantee the continuous arrival time of vehicles at customers. Constraints (23) and (24) guarantee the departure and return times of vehicles at logistics facilities. Constraint (25) ensures the continuous departure time of shared routes of each vehicle. Constraints (26)–(31) indicate the binary restrictions on decision variables.


## 4.1. Related Definitions and Variables

In this section, the related notations and definitions on sets, parameters, and variables are shown in Table2.Table 2
Notations and definitions in the MRVRPRSTWA.
Set definitionPSet of recycling centers (RCs),P = {p|p = 1, 2, 3…, m} and m is the total number of RCsCSet of all recycling customers,C = {c|c = 1, 2, 3…, h} and h is the total number of customersCfSet of customers with fixed time windowsCaSet of customers with assigned time windowsVSet of vehicles for recycling,V = {v|v = 1, 2, 3…, b} and b is the total number of vehiclesSSet of semitrailer trucks used between RCs,S = {s|s = 1, 2, 3…, w} and w is the total number of vehiclesOvMaximum times for executing recycling routes of vehiclev in a working period, v∈VParametersqiRecycling demand quantity of customeri, i∈CtqmnTransportation quantity from logistics facilitym to n, m, n∈PQsMaximum capacity of semitrailer trucks, s∈SQvMaximum capacity of vehiclev, v∈VQmMaximum capacity of logistics facilitym, m∈PUsUsage cost of semitrailer trucks, s∈S (unit: dollar/km)UvUsage cost of vehiclev, v∈V (unit: dollar/km)DijDistance from logistics facility or customeri to j, i, j∈P∪C, i ≠ j (unit: km)DmnDistance from logistics facilitym to logistics facility n, m, n∈P (unit: km)WNumber of working days in one planning periodBNumber of planning periods in one yearMsAnnual MC of semitrailer trucks, s∈SMvAnnual MC of vehiclev, v∈VFLmFixed cost of logistics facilitym, m∈P[ei,ti]Fixed time window of customeri, i∈Cf[e’ i,l’ i]Expected time window of customeri, i∈Ca[αi,βi]Time window assigned to customeri, i∈Ca[Em,Lm]Service time window of logistics facilitym, m∈PμePenalty cost for early arrival per unit timeμdPenalty cost for delayed arrival per unit timeμaCost coefficient of the customer’s time window from the expected time window to the assigned time window per unit timedtk vmDeparture time of thekth route of vehicle V from logistics facility m, v ∈ V, m ∈ P, k ∈ Ovatk ivArrival time of thekth route of vehicle v at node i, v ∈ V, i ∈ P ∪ C, k ∈ OvtijvTravel time of vehiclev between entities i and j, i, j ∈ P ∪ C, v ∈ V|Ns|Total number of logistic facilities served by semitrailer trucks, s∈SNvTotal number of customers served by vehiclev, v∈VBNBig numberDecision variablesxijvIf vehiclev travels from logistics facility or customer i to j, then xijv = 1; otherwise, xijv = 0, i, j ∈ P∪C, v ∈ Vzk vIf vehiclev has the kth route, then zk v = 1; otherwise, zk v = 0, m ∈ P,v ∈ V, and k ∈ OvωivmIf vehiclev departs from logistics facility m to served customer i, ωivm = 1; otherwise, ωivm = 0, i ∈ C, v ∈ V, and m ∈ PymnsIf semitrailer trucks transports between logistics facility m and n, ymnk = 1; otherwise, ymns = 0, m, n ∈ P, and s ∈ SτimnIf customeri’s logistic facility changes from logistics facility m to n after optimization, τimn = 1; otherwise, τimn = 0, i ∈ C, and m, n ∈ PgmsIf semitrailer trucks departs from logistics facility m, gms = 1; otherwise, gms = 0, s ∈ S, and m ∈ P


## 4.2. Mathematical Model

A bi-objective optimization model for MRVRPRSTWA is formulated to minimize the total operating cost in equation (1) and minimize TNV in equation (2).(1)Min TC=TC1+TC2+TC3,(2)Min MNV=∑v∈Vmin∑i∈P∑j∈Cxijv,1.TC1 contains three components in equation (3): ∑m∈P∑n∈P∑s∈SUs×Dmn×W×ymns represents the centralized TCs by semitrailer trucks among logistics facilities. maxm∈P∑n∈Ptqmn/Qs×Ms/B is the MC of semitrailer trucks. ∑m∈PFCm represents the fixed cost of logistics facilities.(3)TC1=∑m∈P∑n∈P∑s∈SUs×Dmn×W×ymns+maxm∈P∑n∈PtqmnQs×MsB+∑m∈PFLm.TC2 contains two components in equation (4): ∑i∈P∪C∑j∈P∪C∑v∈VUv×Dij×W×xijv is the TC by vehicles from logistics facilities to customers. ∑v∈Vmin∑i∈P∑j∈Cxijv,1×Mv/B represents the MC of vehicles.(4)TC2=∑i∈P∪C∑j∈P∪C∑v∈VUv×Dij×W×xijv+∑v∈Vmin∑i∈P∑j∈Cxijv,1×MvB.TC3 contains five components in equation (5): ∑i∈P∪Cf∑j∈Cf∑v∈V∑k∈Ovxijv×μe×maxei−ativk,0 and ∑i∈P∪Ca∑j∈Ca∑v∈V∑k∈Ovxijv×μe×maxαi−ativk,0 represent the penalty cost caused by vehicles that provide earlier service than time windows of customers. ∑i∈P∪Ca∑j∈Ca∑v∈V∑k∈Ovxijv×μd×maxativk−li,0 and ∑i∈P∪Ca∑j∈Ca∑v∈V∑k∈Ovxijv×μd×maxativk−βi,0 represent the penalty cost caused by vehicles that provide delayed service than time windows of customers. ∑i∈Caμa×minei′−αi,βi−li′ denotes the cost of TWA.(5)TC3=∑i∈P∪Cf∑j∈Cf∑v∈V∑k∈Ovxijv×μe×maxei−ativk,0+∑i∈P∪Cf∑j∈Cf∑v∈V∑k∈Ovxijv×μd×maxativk−li,0+∑i∈P∪Ca∑j∈Ca∑v∈V∑k∈Ovxijv×μe×maxαi−ativk,0+∑i∈P∪Ca∑j∈Ca∑v∈V∑k∈Ovxijv×μd×maxativk−βi,0+∑i∈Caμa×minei′−αi,βi−li′.Subject to(6)∑m∈Pymns=1,∀n∈P,s∈S,m≠n,(7)∑n∈Pymns−∑n∈Pynfs=0,∀m,f∈P,m≠f,s∈S,(8)∑m∈Pynms=1,∀n∈P,s∈S,m≠n,(9)∑m∈P∑n∈Pymns×gps≤Ns−1,∀p∈P,s∈S,p≠m,p≠n,(10)tqmn=∑i∈Cτimn×qi,∀m,n∈P,(11)∑m∈P∑n∈Ptqmn×ymns≤Qs,∀s∈S,(12)∑i∈Pxijv=1,∀j∈C,v∈V,(13)∑i∈C∑j∈Cxijv=1,∀v∈V,(14)∑j∈Pxijv−∑j∈Pxjfv=0,∀i,f∈P∪C,i≠f,v∈V,(15)∑i∈Pxjiv=1,∀j∈C,v∈V,(16)∑i∈C∑j∈Cxijv≤Nv−1,∀v∈V,(17)∑j∈Cqj×xijv×zvk≤Qv,∀i∈P∪C,v∈V,k∈Ov,(18)∑j∈C∑v∈Vqj×ωijvm≤Qm,∀i∈P∪C,m∈P,(19)dtvmk+tmiv−BN1−xmiv≤ativk,∀m∈P,i∈C,v∈V,k∈Ov,(20)dtvmk+tmiv+BN1−xmiv≥ativk,∀m∈P,i∈C,v∈V,k∈Ov,(21)ativk+tijv−BN1−xijv≤atjvk,∀i∈C,j∈C∪P,v∈V,k∈Ov,(22)ativk+tijv+BN1−xijv≥atjvk,∀i∈C,j∈C∪P,v∈V,k∈Ov,(23)Em×xmiv×zvk≤dtvmk≤Lm×xmiv×zvk,∀m∈P,i∈C,v∈V,k∈Ov,(24)Em×xmiv×zvk≤atmvk≤Lm×xmiv×zvk,∀m∈P,i∈C,v∈V,k∈Ov,(25)atmvk+tmnv−M1−zvk+1≤dtvnk+1,∀m,n∈P,v∈V,k∈Ov,(26)xijv=0,1,∀i,j∈P∪C,v∈V,(27)zvk=0,1,∀v∈V,k∈Ov,(28)ωivm=0,1,∀i∈C,v∈V,m∈P,(29)ymns=0,1,∀m,n∈P,∀s∈S,(30)τimn=0,1,∀i∈C,∀m,n∈P,m≠n,(31)gms=0,1,∀m∈P,s∈S.Constraint (6) ensures that each semitrailer truck originally departs from logistics facilities. Constraint (7) is the flow conservation on each logistics facility. Constraint (8) ensures that each semitrailer truck finally returns to logistics facilities. Constraint (9) is used to avoid subtours for semitrailer trucks. Constraint (10) calculates the transportation quantity among logistics facilities. Constraint (11) stipulates that each semitrailer truck has enough capacity to satisfy the transportation quantity of the assigned logistics facilities. Constraint (12) ensures that each vehicle departs from logistics facilities only. Constraint (13) ensures that each customer can be served once by one vehicle. Constraint (14) is the flow conservation on each customer. Constraint (15) ensures that each vehicle finally returns to logistics facilities. Constraint (16) is used to avoid subtours for vehicles. Constraint (17) represents that each vehicle has enough capacity to satisfy the total demand of assigned customers. Constraint (18) ensures that the total demand of customers served by a logistics facility should be beyond its capacity. Constraints (19) and (20) guarantee the continuous departure time of vehicles at logistics facilities. Constraints (21) and (22) guarantee the continuous arrival time of vehicles at customers. Constraints (23) and (24) guarantee the departure and return times of vehicles at logistics facilities. Constraint (25) ensures the continuous departure time of shared routes of each vehicle. Constraints (26)–(31) indicate the binary restrictions on decision variables.


## 5. Solution Procedure

The hybrid heuristic algorithm is widely applied to solve multiobjective optimization problems [59–61]. This algorithm, including the 3D k-means clustering and NSGA-II algorithms, is proposed to solve the MRVRPRSTWA. The 3D k-means clustering algorithm is often utilized to construct clusters according to the spatial and temporal distances of data [53]. The 3D k-means clustering algorithm can effectively reduce the computational complexity of multidepot logistics networks [37]. The NSGA-II can assign appropriate time windows to customers and obtain Pareto optimal solutions in MRVRPRSTWA. The elitist strategy in NSGA-II can effectively avoid the loss of the best individual and enhance the computational speed and robustness. The elitist strategy is composed of the genetic operations (i.e., selection, crossover, and mutation) and the fast nondominated sorting and crowding distance. The parameters utilized in the flow structure are defined as follows: i is the number of RCs and clustering centers, R is set as the present number of optimization runs, Rmax represents the maximum number of optimization runs, N is set as the present number of generations, and Nmax is set as the maximum number of generations. The flow structure of the hybrid heuristic algorithm is shown in Figure 3.Figure 3
Algorithm flow structure for MRVRPRSTWA.The optimization procedure of MRVRPRSTWA is as follows:Step 1: 3D k-means clustering. Set the number of RCs as the number of clusters and the data from RCs as the initial data of each cluster. Calculate the spatial and temporal distances between cluster centers and customers. Assign each customer to the nearest cluster center.Step 2: the NSGA-II algorithm is used to establish pickup vehicle routes. The initial parameters are set and the initial population is randomly generated.Step 3: genetic operation, including selection, partial-mapped crossover (PMC), and mutation is performed to generate offspring population. The parent and offspring population are combined.Step 4: the new population is selected by performing nondominated sorting and calculating the crowding distance of each individual.Step 5: determine whether the internal termination is satisfied. If the termination condition is unsatisfied, then select customers that accepted the TWA strategy and have time window violations, and assign the appropriate time windows to them; otherwise, return to Step 3.Step 6: calculate the total cost and TNV with the adoption of the TWA strategy by NSGA-II and determine whether the termination condition is satisfied. If the termination condition is unsatisfied, then return to Step 1 to confirm whether the clustering results need to be adjusted; otherwise, terminate the algorithm procedure and find the optimal solution.


### 5.1. 3Dk-Means Clustering

The clustering algorithm is always utilized to construct clusters to each depot in multidepot VRP to reduce computational complexity [62, 63]. The traditional k-means clustering algorithm groups customers on the basis of the spatial distances between customers and clustering centers [64, 65]. The 3D k-means clustering algorithm, a variant of the traditional k-means clustering algorithm, constructs clusters according to the spatial and temporal distances among clustering centers and customers. In this study, the 3D k-means clustering algorithm is applied to optimize MRVRPRSTWA. The spatial and temporal distances based on geographic coordinates and time windows of customers and logistics facilities can be calculated by(32)β1+β2=1,(33)dis=β1xm−xp+ym−yp+β2zm−zp.In equation (32), β1 and β2 represent the proportion of spatial and temporal distances, respectively. In equation (33), (xm, ym, zm) expresses the geographic location and center point of time windows of logistics facility. (xp, yp, zp) expresses the geographic location and center point of time windows of customer. The 3D k-means clustering algorithm procedure is shown in Algorithm 1.Algorithm 1: 3D k-means clustering algorithm.
Input: Number of RCs, geographic coordinates, and time windows of customers and RCsOutput: Clusters, including RCs and their assigned customers(1)
Import the geographic locations and center point of the time windows of RCs and customers as initial data.(2)
Establish 3D vectors on the basis of the initial data and construct 3D data elements.(3)
Seti as the number of clusters on the basis of the number of RCs.(4)
Select the data from RCs as the initial clustering center of each cluster.(5)
Traverse all the data elements.(6)
Calculate the spatial and temporal distances from each element to the clustering center.(7)
Assign each element to the nearest clustering center.(8)
Updatei clusters, and calculate new clustering centers.(9)
If existing customers need to be adjusted among clusters(10)
Then return to step 5;(11)
Else(12)
End if(13)
End for(14)
Output the results ofi clusters.In Algorithm1, the process of clustering can be divided into several components: First, the geographic locations and center point of the time windows of RCs and customers are used to construct the initial 3D data elements. Second, set i on the basis of the number of RCs as the number of clusters and the data from RCs as the initial clustering center of each cluster. Third, traverse data elements and calculate the spatial and temporal distances from each element to every clustering center. Fourth, assign each element to the nearest clustering center, and elements assigned to the same clustering center are placed in one cluster. Fifth, update i clusters, and continue to execute this procedure till no existing customers need to be adjusted among clusters; then, i clusters, including RCs and their served customers, are reported. The clustering results provide the initial population for multidepot vehicle routing optimization.


### 5.2. NSGA-II

Heuristic algorithms have obvious performances in solving multidepot logistics networks [66]. The NSGA-II is an enhanced algorithm based on the NSGA to solve multiobjective evolutionary problems [67, 68]. NSGA-II ranks different individuals via nondominated sorting algorithm and crowding distance comparison operators to force individuals to converge toward and diffuse along the Pareto front. The main procedures of NSGA-II are summarized below.


#### 5.2.1. Selection

The selection operation involves choosing certain individuals every time from the parent population with selection probabilitysp, and the best individual will be inserted into the offspring population. The selection operation will continue until the offspring population size reaches a predetermined size. The specific steps can be described as follows: First, individuals are chosen from the parent population with selection probability sp at random. Then, the individual with the best fitness values is selected and added into the offspring population. The abovementioned procedures are repeated until the offspring population reaches a predetermined size.


#### 5.2.2. PMC

PMC is used as the operator of the crossover operation [69]. The specific procedure of PMC is shown in Figure 4, and the main steps of PMC are as follows: First, randomly select the beginning and ending positions of a gene sequence of a pair of chromosomes. Second, the two sets of genes selected on the two chromosomes are exchanged separately. Third, conduct conflict detection and establish a mapping relationship on the exchanged genes. Fourth, repeat this step until no duplicate genes are present on the same chromosome.Figure 4
Procedure of PMC.The mapping relationship is illustrated in Figure4. For example, a mapping relationship of 3-1and 5-6 genes is presented. Two genes 3 in offspring1 are selected in the second step, and the unselected gene is transformed into 1. Two genes 5 in offspring 1 are selected in the third step, and the unselected gene is transformed into 6. Finally, two children with no conflict genes are generated.


#### 5.2.3. Mutation

Mutation operation aims to change the values of a gene at a certain locus of chromosome to generate a new chromosome [70]. Polynomial mutation is used in this work to prevent local convergence in the course of evolution. The main steps of polynomial mutation are shown in Figure 5.Figure 5
Procedure of mutation.In Figure5, a chromosome is randomly selected from the parent generation in terms of the mutation probability pc. Then, three genes on the parental chromosome are randomly selected, and the exchange and regenerated operation are conducted on the three genes. If duplicate genes appear, then they should be removed and regenerated.


#### 5.2.4. Nondominated Sorting and Crowding Distance Assignment

Several definitions are provided to describe the procedures of nondominated sorting and crowding distance.N represents the set of population, and it can be divided into several subsets (i.e., A1, A2, A3…). The subsets can be indicated as follows: A1 is the set of nondominated population of N, A2 is the set of nondominated population of N–A1, A3 is the set of nondominated population of N–(A1∪A2), and so on. In each individual m, let id(m) be the index of the front to which m belongs and cd(m) be the crowding distance of m. The crowding distance can be calculated in equation (34) as follows:(34)cdm=f1m1−f1m2f1max−f1min+f2m3−f2m4f2max−f2min,where m1 and m2 are the closest populations of m in the same front on the basis of f1; m3 and m4 are the closest populations of m in the same front on the basis of f2; and fmax and fmin are the maximum and minimum values of fi, respectively. In addition, if n is the boundary individual, then the crowding distance of n is infinity, namely, cd(n) = ∞. If the id(m) is smaller than the id(m’) or id(m) is equal to the id(m’), and the cd(m) is larger than the cd(m’), then we can say m is preferred to m’’. The next generation is selected through the abovementioned method, which is nondominated sorting and crowding distance [71].


#### 5.2.5. TWA Strategy

The TWA strategy is designed in the NSGA-II algorithm to assign candidate time windows to customers with improved operational efficiency. The main procedures are presented as follows:Step 1: customers that accepted the TWA strategy with violation of time windows are selected.Step 2: assume that the time window for customer A is [EA, LA], and the actual service time of customer A is i hours earlier (or delayed) than EA. If EA–i (or LA + i) is between the candidate time window [E' A, L' A], then calculate the TWA cost, and compare it with the penalty cost.Step 3: if the TWA cost is lower than the penalty cost, then [E' A, L' A] will assign to customer A; otherwise, keep the initial time window of customer A.Therefore, the appropriate time windows are assigned to the corresponding customers. The process of NSGA-II is described in Algorithm2.Algorithm 2: NSGA-II algorithm operation.
Nondominated Sorting Algorithm-II (NSGA-II)Input: popsize, nodes, R, Rmax, Nmax, sp, cp, and mpOutput: Pareto front optimal solutions(1)
Initialize parameters(2)
# set the population size (popsize), number of customers (nodes), number of generation (R), maximum number of generations (Rmax), maximum number of runs (Nmax), selection probability (sp), crossover probability (cp), and mutation probability (mp)(3)
ForN = 1:Nmax(4)
For R = 1:Rmax(5)
Generate the initial population with size popsize(6)
Objective function evaluation(7)
# compute the objective function to minimize the total operating cost and TNV(8)
Divide popsize into nondominance front and calculate the crowding distance of each individual(9)
For i = 1: R(10)
Implement selection, PMC, and mutation operations to generate offspring population(11)
Combine parent and offspring population, perform nondominated sorting and calculate the crowding distance of each individual(12)
Build the new population on the basis of the parent and offspring populations following partial order(13)
end(14)
end(15)
For n = 1: nodes(16)
Select the customers that accept time window and have time window violation(17)
Select the appropriate time window from candidate time windows to assign to customer(18)
Calculate the TWA cost and compare it with the penalty cost(19)
Assign the time window to the customer or keep the expected time window of customer(20)
end(21)
Find the Pareto front optimal solution(22)
end


## 5.1. 3Dk-Means Clustering

The clustering algorithm is always utilized to construct clusters to each depot in multidepot VRP to reduce computational complexity [62, 63]. The traditional k-means clustering algorithm groups customers on the basis of the spatial distances between customers and clustering centers [64, 65]. The 3D k-means clustering algorithm, a variant of the traditional k-means clustering algorithm, constructs clusters according to the spatial and temporal distances among clustering centers and customers. In this study, the 3D k-means clustering algorithm is applied to optimize MRVRPRSTWA. The spatial and temporal distances based on geographic coordinates and time windows of customers and logistics facilities can be calculated by(32)β1+β2=1,(33)dis=β1xm−xp+ym−yp+β2zm−zp.In equation (32), β1 and β2 represent the proportion of spatial and temporal distances, respectively. In equation (33), (xm, ym, zm) expresses the geographic location and center point of time windows of logistics facility. (xp, yp, zp) expresses the geographic location and center point of time windows of customer. The 3D k-means clustering algorithm procedure is shown in Algorithm 1.Algorithm 1: 3D k-means clustering algorithm.
Input: Number of RCs, geographic coordinates, and time windows of customers and RCsOutput: Clusters, including RCs and their assigned customers(1)
Import the geographic locations and center point of the time windows of RCs and customers as initial data.(2)
Establish 3D vectors on the basis of the initial data and construct 3D data elements.(3)
Seti as the number of clusters on the basis of the number of RCs.(4)
Select the data from RCs as the initial clustering center of each cluster.(5)
Traverse all the data elements.(6)
Calculate the spatial and temporal distances from each element to the clustering center.(7)
Assign each element to the nearest clustering center.(8)
Updatei clusters, and calculate new clustering centers.(9)
If existing customers need to be adjusted among clusters(10)
Then return to step 5;(11)
Else(12)
End if(13)
End for(14)
Output the results ofi clusters.In Algorithm1, the process of clustering can be divided into several components: First, the geographic locations and center point of the time windows of RCs and customers are used to construct the initial 3D data elements. Second, set i on the basis of the number of RCs as the number of clusters and the data from RCs as the initial clustering center of each cluster. Third, traverse data elements and calculate the spatial and temporal distances from each element to every clustering center. Fourth, assign each element to the nearest clustering center, and elements assigned to the same clustering center are placed in one cluster. Fifth, update i clusters, and continue to execute this procedure till no existing customers need to be adjusted among clusters; then, i clusters, including RCs and their served customers, are reported. The clustering results provide the initial population for multidepot vehicle routing optimization.


## 5.2. NSGA-II

Heuristic algorithms have obvious performances in solving multidepot logistics networks [66]. The NSGA-II is an enhanced algorithm based on the NSGA to solve multiobjective evolutionary problems [67, 68]. NSGA-II ranks different individuals via nondominated sorting algorithm and crowding distance comparison operators to force individuals to converge toward and diffuse along the Pareto front. The main procedures of NSGA-II are summarized below.


### 5.2.1. Selection

The selection operation involves choosing certain individuals every time from the parent population with selection probabilitysp, and the best individual will be inserted into the offspring population. The selection operation will continue until the offspring population size reaches a predetermined size. The specific steps can be described as follows: First, individuals are chosen from the parent population with selection probability sp at random. Then, the individual with the best fitness values is selected and added into the offspring population. The abovementioned procedures are repeated until the offspring population reaches a predetermined size.


### 5.2.2. PMC

PMC is used as the operator of the crossover operation [69]. The specific procedure of PMC is shown in Figure 4, and the main steps of PMC are as follows: First, randomly select the beginning and ending positions of a gene sequence of a pair of chromosomes. Second, the two sets of genes selected on the two chromosomes are exchanged separately. Third, conduct conflict detection and establish a mapping relationship on the exchanged genes. Fourth, repeat this step until no duplicate genes are present on the same chromosome.Figure 4
Procedure of PMC.The mapping relationship is illustrated in Figure4. For example, a mapping relationship of 3-1and 5-6 genes is presented. Two genes 3 in offspring1 are selected in the second step, and the unselected gene is transformed into 1. Two genes 5 in offspring 1 are selected in the third step, and the unselected gene is transformed into 6. Finally, two children with no conflict genes are generated.


### 5.2.3. Mutation

Mutation operation aims to change the values of a gene at a certain locus of chromosome to generate a new chromosome [70]. Polynomial mutation is used in this work to prevent local convergence in the course of evolution. The main steps of polynomial mutation are shown in Figure 5.Figure 5
Procedure of mutation.In Figure5, a chromosome is randomly selected from the parent generation in terms of the mutation probability pc. Then, three genes on the parental chromosome are randomly selected, and the exchange and regenerated operation are conducted on the three genes. If duplicate genes appear, then they should be removed and regenerated.


### 5.2.4. Nondominated Sorting and Crowding Distance Assignment

Several definitions are provided to describe the procedures of nondominated sorting and crowding distance.N represents the set of population, and it can be divided into several subsets (i.e., A1, A2, A3…). The subsets can be indicated as follows: A1 is the set of nondominated population of N, A2 is the set of nondominated population of N–A1, A3 is the set of nondominated population of N–(A1∪A2), and so on. In each individual m, let id(m) be the index of the front to which m belongs and cd(m) be the crowding distance of m. The crowding distance can be calculated in equation (34) as follows:(34)cdm=f1m1−f1m2f1max−f1min+f2m3−f2m4f2max−f2min,where m1 and m2 are the closest populations of m in the same front on the basis of f1; m3 and m4 are the closest populations of m in the same front on the basis of f2; and fmax and fmin are the maximum and minimum values of fi, respectively. In addition, if n is the boundary individual, then the crowding distance of n is infinity, namely, cd(n) = ∞. If the id(m) is smaller than the id(m’) or id(m) is equal to the id(m’), and the cd(m) is larger than the cd(m’), then we can say m is preferred to m’’. The next generation is selected through the abovementioned method, which is nondominated sorting and crowding distance [71].


### 5.2.5. TWA Strategy

The TWA strategy is designed in the NSGA-II algorithm to assign candidate time windows to customers with improved operational efficiency. The main procedures are presented as follows:Step 1: customers that accepted the TWA strategy with violation of time windows are selected.Step 2: assume that the time window for customer A is [EA, LA], and the actual service time of customer A is i hours earlier (or delayed) than EA. If EA–i (or LA + i) is between the candidate time window [E' A, L' A], then calculate the TWA cost, and compare it with the penalty cost.Step 3: if the TWA cost is lower than the penalty cost, then [E' A, L' A] will assign to customer A; otherwise, keep the initial time window of customer A.Therefore, the appropriate time windows are assigned to the corresponding customers. The process of NSGA-II is described in Algorithm2.Algorithm 2: NSGA-II algorithm operation.
Nondominated Sorting Algorithm-II (NSGA-II)Input: popsize, nodes, R, Rmax, Nmax, sp, cp, and mpOutput: Pareto front optimal solutions(1)
Initialize parameters(2)
# set the population size (popsize), number of customers (nodes), number of generation (R), maximum number of generations (Rmax), maximum number of runs (Nmax), selection probability (sp), crossover probability (cp), and mutation probability (mp)(3)
ForN = 1:Nmax(4)
For R = 1:Rmax(5)
Generate the initial population with size popsize(6)
Objective function evaluation(7)
# compute the objective function to minimize the total operating cost and TNV(8)
Divide popsize into nondominance front and calculate the crowding distance of each individual(9)
For i = 1: R(10)
Implement selection, PMC, and mutation operations to generate offspring population(11)
Combine parent and offspring population, perform nondominated sorting and calculate the crowding distance of each individual(12)
Build the new population on the basis of the parent and offspring populations following partial order(13)
end(14)
end(15)
For n = 1: nodes(16)
Select the customers that accept time window and have time window violation(17)
Select the appropriate time window from candidate time windows to assign to customer(18)
Calculate the TWA cost and compare it with the penalty cost(19)
Assign the time window to the customer or keep the expected time window of customer(20)
end(21)
Find the Pareto front optimal solution(22)
end


## 5.2.1. Selection

The selection operation involves choosing certain individuals every time from the parent population with selection probabilitysp, and the best individual will be inserted into the offspring population. The selection operation will continue until the offspring population size reaches a predetermined size. The specific steps can be described as follows: First, individuals are chosen from the parent population with selection probability sp at random. Then, the individual with the best fitness values is selected and added into the offspring population. The abovementioned procedures are repeated until the offspring population reaches a predetermined size.


## 5.2.2. PMC

PMC is used as the operator of the crossover operation [69]. The specific procedure of PMC is shown in Figure 4, and the main steps of PMC are as follows: First, randomly select the beginning and ending positions of a gene sequence of a pair of chromosomes. Second, the two sets of genes selected on the two chromosomes are exchanged separately. Third, conduct conflict detection and establish a mapping relationship on the exchanged genes. Fourth, repeat this step until no duplicate genes are present on the same chromosome.Figure 4
Procedure of PMC.The mapping relationship is illustrated in Figure4. For example, a mapping relationship of 3-1and 5-6 genes is presented. Two genes 3 in offspring1 are selected in the second step, and the unselected gene is transformed into 1. Two genes 5 in offspring 1 are selected in the third step, and the unselected gene is transformed into 6. Finally, two children with no conflict genes are generated.


## 5.2.3. Mutation

Mutation operation aims to change the values of a gene at a certain locus of chromosome to generate a new chromosome [70]. Polynomial mutation is used in this work to prevent local convergence in the course of evolution. The main steps of polynomial mutation are shown in Figure 5.Figure 5
Procedure of mutation.In Figure5, a chromosome is randomly selected from the parent generation in terms of the mutation probability pc. Then, three genes on the parental chromosome are randomly selected, and the exchange and regenerated operation are conducted on the three genes. If duplicate genes appear, then they should be removed and regenerated.


## 5.2.4. Nondominated Sorting and Crowding Distance Assignment

Several definitions are provided to describe the procedures of nondominated sorting and crowding distance.N represents the set of population, and it can be divided into several subsets (i.e., A1, A2, A3…). The subsets can be indicated as follows: A1 is the set of nondominated population of N, A2 is the set of nondominated population of N–A1, A3 is the set of nondominated population of N–(A1∪A2), and so on. In each individual m, let id(m) be the index of the front to which m belongs and cd(m) be the crowding distance of m. The crowding distance can be calculated in equation (34) as follows:(34)cdm=f1m1−f1m2f1max−f1min+f2m3−f2m4f2max−f2min,where m1 and m2 are the closest populations of m in the same front on the basis of f1; m3 and m4 are the closest populations of m in the same front on the basis of f2; and fmax and fmin are the maximum and minimum values of fi, respectively. In addition, if n is the boundary individual, then the crowding distance of n is infinity, namely, cd(n) = ∞. If the id(m) is smaller than the id(m’) or id(m) is equal to the id(m’), and the cd(m) is larger than the cd(m’), then we can say m is preferred to m’’. The next generation is selected through the abovementioned method, which is nondominated sorting and crowding distance [71].


## 5.2.5. TWA Strategy

The TWA strategy is designed in the NSGA-II algorithm to assign candidate time windows to customers with improved operational efficiency. The main procedures are presented as follows:Step 1: customers that accepted the TWA strategy with violation of time windows are selected.Step 2: assume that the time window for customer A is [EA, LA], and the actual service time of customer A is i hours earlier (or delayed) than EA. If EA–i (or LA + i) is between the candidate time window [E' A, L' A], then calculate the TWA cost, and compare it with the penalty cost.Step 3: if the TWA cost is lower than the penalty cost, then [E' A, L' A] will assign to customer A; otherwise, keep the initial time window of customer A.Therefore, the appropriate time windows are assigned to the corresponding customers. The process of NSGA-II is described in Algorithm2.Algorithm 2: NSGA-II algorithm operation.
Nondominated Sorting Algorithm-II (NSGA-II)Input: popsize, nodes, R, Rmax, Nmax, sp, cp, and mpOutput: Pareto front optimal solutions(1)
Initialize parameters(2)
# set the population size (popsize), number of customers (nodes), number of generation (R), maximum number of generations (Rmax), maximum number of runs (Nmax), selection probability (sp), crossover probability (cp), and mutation probability (mp)(3)
ForN = 1:Nmax(4)
For R = 1:Rmax(5)
Generate the initial population with size popsize(6)
Objective function evaluation(7)
# compute the objective function to minimize the total operating cost and TNV(8)
Divide popsize into nondominance front and calculate the crowding distance of each individual(9)
For i = 1: R(10)
Implement selection, PMC, and mutation operations to generate offspring population(11)
Combine parent and offspring population, perform nondominated sorting and calculate the crowding distance of each individual(12)
Build the new population on the basis of the parent and offspring populations following partial order(13)
end(14)
end(15)
For n = 1: nodes(16)
Select the customers that accept time window and have time window violation(17)
Select the appropriate time window from candidate time windows to assign to customer(18)
Calculate the TWA cost and compare it with the penalty cost(19)
Assign the time window to the customer or keep the expected time window of customer(20)
end(21)
Find the Pareto front optimal solution(22)
end


## 6. Implementation and Analysis

### 6.1. Algorithm Comparison

In this section, benchmark instances are performed to test the quality of NSGA-II in MRVRPRSTWA for comparing its results with those received by multiobjective evolutionary algorithm (MOEA) [72] and multiobjective particle swarm optimization (MOPSO) [73]. The multidepot VRPTW (MDVRPTW) datasets obtained from the database of the NEO research group (https://neo.lcc.uma.es/vrp/vrp-instances/multiple-depot-vrp-with-time-windows-instances/) are illustrated in Table 3.Table 3
Relevant settings of data instances.
InstancesDatasetsNumber of depotsNumber of customersNumber of candidate time windowsVehicle capacity1MDPR144822002MDPR249621953MDPR3414431904MDPR4419241855MDPR5424041806MDPR6428851757MDPR767222008MDPR8614431909MDPR96216418010MDPR106288517011MDPR11448220012MDPR12496219513MDPR134144319014MDPR144192418515MDPR154240418016MDPR164288517517MDPR17672220018MDPR186144319019MDPR196216418020MDPR2062885170Table3 presents the relevant characteristics of 20 benchmarks. The parameters for NSGA-II and MOEA are as follows: population size popsize = 150, selection probability sp = 0.6, crossover probability cp = 0.9, and mutation probability mp = 0.1. The parameters for MOPSO are as follows: maximum iteration number max_it = 200, inertia weight ω = 0.9, personal confidences pc = 2, and social learning confidences gc = 3. The results of the three algorithms, which contain total cost (cost), TNV (vehicle), and computation time (Ct), are compared in Table 4.Table 4
Comparison of the results of the three algorithms on different instances.
InstancesNSGA-IIMOEAMOPSOCost ($)VehicleCt (s)Cost ($)VehicleCt (s)Cost ($)VehicleCt (s)1845031841274271891068951792961741951490692211283772043137168210161371021115067921541576410225189401324319923142315158291023121070152381804112244617032112512141215267193381324771695292502037312255173259261820048122542360115248185501026991796810266255261726122484142681019119112782350015287194541128811833532071472092108673321712105985216158401023114772923113127437224191651324217094112331415789102302011014257180531224615159291023721171152411711911251161807412241202601424619231132651715932825121234132792017312271181815810257194111126419367112791919998122622448216298214611327820195371126825779173172074212288Average15479923720019132501752011248t-test−12.955−5.529p-value3.51502E-111.24E-05Table4 exhibits an adequate difference among the optimal results of those three algorithms on the basis of the values of t-test and p-value. The NSGA-II algorithm is superior to MOEA and MOPSO in obtaining the optimal cost in MRVRPRSTWA. The costs of MOEA and MOPSO are higher than that of NSGA-II in 20 benchmark instances. For instance, the average cost of NSGA-II is $15479. However, the costs of MOEA and MOPSO are $20019 and $17520, respectively. The minimum number of vehicles calculated by NSGA-II is nine compared to TNV of MOEA and MOPSO. The NSGA-II outperforms MOEA and MOPSO in terms of the minimum computation time. The average computation time of NSGA-II is 237s, which is lower than the computation times of MOEA and MOPSO. Therefore, the NSGA-II algorithm can achieve better performance in solving MRVRPRSTWA compared to MOEA and MOPSO.


### 6.2. Data Source

A real reverse logistics network in Chongqing, China, is used to verify the applicability and feasibility of the proposed solution in MRVRPRSTWA. The logistics network has five RCs (i.e., RC1, RC2, RC3, RC4, and RC5) and 183 customers. The spatial distribution of logistics network is shown in Figure6. The RCs and customers are marked as stars and rhombuses, respectively. Rhombuses in various colors represent customers served by different RCs. The initial vehicle routes are shown in Table 5.Figure 6
Spatial distribution of RCs and customers.Table 5
Initial vehicle routes before optimization.
FacilitiesRoutesRC1RC1⟶P175⟶P174⟶P167⟶P157⟶P163⟶P172⟶P169⟶P176⟶P162⟶P179⟶P155⟶RC1RC1⟶P180⟶P177⟶P171⟶RC1RC1⟶P150⟶P159⟶P152⟶P165⟶P168⟶P161⟶P170⟶P154⟶P181⟶P158⟶RC1RC1⟶P182⟶P164⟶P156⟶P183⟶P160⟶P153⟶P166⟶P151⟶P178⟶P173⟶RC1RC2RC2⟶P98⟶P107⟶P95⟶P97⟶P99⟶P93⟶P94⟶P84⟶P85⟶P104⟶P103⟶RC2RC2⟶P82⟶P100⟶P89⟶P90⟶P101⟶P96⟶P88⟶P105⟶P87⟶P102⟶P92⟶RC2RC2⟶P91⟶P86⟶P106⟶P83⟶RC2RC3RC3⟶P122⟶121⟶P113⟶P114⟶P115⟶P108⟶P131⟶RC3RC3⟶P124⟶P149⟶P112⟶P147⟶P140⟶P144⟶P141⟶RC3RC3⟶P126⟶P143⟶P120⟶P145⟶P123⟶P125⟶P128⟶P139⟶RC3RC3⟶P148⟶P137⟶P110⟶P138⟶P118⟶P134⟶P136⟶P142⟶P132⟶P129⟶RC3RC3⟶P135⟶P146⟶P133⟶P109⟶P111⟶P116⟶P117⟶P119⟶P130⟶P127⟶RC3RC4RC4⟶P22⟶P16⟶P40⟶P47⟶R42RC4⟶P14⟶P5⟶P9⟶P3⟶P36⟶P33⟶P18⟶P12⟶P30⟶P11⟶P4⟶P35⟶P6⟶RC4RC4⟶P15⟶P7⟶P39⟶P28⟶P37⟶P38⟶P23⟶RC4RC4⟶P2⟶P26⟶P10⟶P24⟶P17⟶P8⟶P45⟶P20⟶P29⟶P19⟶P21⟶RC4RC4⟶P13⟶P1⟶P32⟶P25⟶P34⟶P41⟶P44⟶P43⟶P27⟶P46⟶P31⟶RC4RC5RC5⟶P52⟶P49⟶P58⟶P66⟶P79⟶P57⟶P68⟶P50⟶P53⟶P51⟶P55⟶P71⟶RC5RC5⟶P59⟶P74⟶P76⟶P62⟶P61⟶P70⟶P78⟶P60⟶P80⟶P54⟶P67⟶P56⟶RC5RC5⟶P64⟶P81⟶P77⟶P48⟶P69⟶P73⟶P72⟶P63⟶RC5RC5⟶P75⟶P65⟶RC5In Table5, customers are served by 21 vehicle routes among five RCs. Each RC is operated independently, which results in the waste of resources. For example, three customers (i.e., P180, P177, and P171) belonging to RC1 and two customers (i.e., P75 and P65) belonging to RC5 are served by a vehicle, respectively. The irrational vehicle routes should be decreased to enhance the operational efficiency of the multidepot logistics network.


### 6.3. Relevant Parameter Setting

The objective of MRVRPRSTWA is to obtain the lowest operating cost and TNV in the multidepot reverse logistics network. The related parameters used in the optimization model and the hybrid heuristic algorithm initialization are shown in Table6. The time window of RC1 is [6, 18]. RC2 and RC3 have the same time window [8, 20]. The time window of RC4 and RC5 is the same as [10, 22]. The candidate time windows are constructed on the basis of the real situation of customers that accepted TWA. Hence, four candidate time windows are constructed: [8, 10], [10, 11], [12, 15], and [16, 18].Table 6
Parameter values utilized in the computational experiments.
NotationDefinitionβ1Coefficient of spatial distanceβ2Coefficient of temporal distanceQsMaximum capacity of the semitrailer truckQvMaximum capacity of vehicleUsUsage cost of the semitrailer truck (unit: dollar/km)UvUsage cost of vehiclev (unit: dollar/km)MsAnnual MC of the semitrailer truckMvAnnual MC of the vehicleμePenalty cost for early arrival per unit timeμdPenalty cost for delayed arrival per unit timeμaAC per unit timeWNumber of working days in one planning periodBNumber of planning periods in one yearFL1Fixed cost of RC1FL2Fixed cost of RC2FL3Fixed cost of RC3FL4Fixed cost of RC4FL5Fixed cost of RC5PopsizePopulation sizeNmaxMaximum number of generations runsRmaxMaximum number of runsspSelection probabilitycpCrossover probabilitympMutation probability


### 6.4. Optimization Results

Customers are clustered on the basis of the 3Dk-means clustering algorithm to reduce the computational complexity in the multidepot reverse logistics network. Customers are reassigned to the RC with the nearest spatial and temporal distances. The results of the 3D k-means clustering algorithm are shown in Table 7 and Figure 7.Table 7
3Dk-means clustering results in MRVRPRSTWA.
FacilitiesRC1RC2RC3RC4RC5CustomersP7 P15 P39 P42 P48 P69 P72 P73 P78 P93 P99 P110 P118 P138 P130 P131 P177 P180 P182 P33 P36 P40 P43 P44 P49 P126 P128 P132P11 P17 P24 P25 P32 P34 P53 P91 P92 P97 P102 P133 P135 P146 P148 P160 P183 P115 P116 P117 P119 P120 P61 P62 P70 P66 P67 P71P1 P2 P4 P6 P10 P13 P26 P28 P29 P35 P41 P50 P51 P60 P74 P76 P82 P83 P86 P87 P95 P98 P103 P104 P106 P107 P109 P112 P113 P121 P122 P124 P129 P136 P137 P149 P153 P154 P156 P158 P164 P166P8 P14 P16 P19 P20 P21 P22 P23 P31 P37 P38 P45 P46 P47 P55 P59 P63 P68 P80 P84 P85 P89 P94 P114 P123 P134 P139 P141 P142 P144 P145 P147 P151 P155 P161 P165 P168 P170 P171 P173 P179 P181P3 P5 P9 P12 P18 P27 P30 P52 P54 P56 P57 P58 P64 P65 P75 P77 P79 P81 P88 P90 P96 P100 P101 P105 P108 P111 P125 P127 P140 P143 P150 P152 P157 P159 P162 P163 P167 P169 P172 P174 P175 P176 P178Figure 7
Customer clustering results in MRVRPRSTWA.Table7 illustrates that 183 customers are reassigned to the corresponding logistics facilities and grouped into five clusters. RC1 is allocated 28 customers, RC2 is allocated 28 customers, RC3 is allocated 42 customers, and RC4 and RC5 are allocated 42 and 43 customers, respectively. The spatial and temporal distances of the customers are shown in Figure 7, wherein x represents the longitude of coordinate, y indicates the latitude of coordinate, and z denotes the mid value of the time window. The optimized vehicle routes with the RS and TWA strategies are shown in Table 8.Table 8
Vehicle routes with the RS and TWA strategies.
VehiclesFacilitiesRoutesVehicle 1RC1RC1⟶P128⟶P138⟶P7⟶P69⟶P78⟶P15⟶P36⟶P42⟶RC1RC3RC3⟶P1⟶P122⟶P50⟶P95⟶P109⟶P51⟶P41⟶P129⟶P13⟶P166⟶P158⟶RC3RC3RC3⟶P98⟶P124⟶P86⟶P112⟶P82⟶P136⟶P28⟶P137⟶P113⟶P104⟶RC3Vehicle 2RC2RC2⟶P133⟶P11⟶P160⟶P17⟶P71⟶P115⟶P34⟶P148⟶P120⟶P32⟶P91⟶RC2RC2RC2⟶P102⟶P24⟶P25⟶P146⟶P119⟶P183⟶P70⟶P9⟶P97⟶RC2RC2RC2⟶P61⟶P62⟶P135⟶P116⟶P53⟶P117⟶P66⟶P67⟶RC2Vehicle 3RC3RC3⟶P106⟶P83⟶P87⟶P29⟶P103⟶P10⟶P76⟶P74⟶P4⟶P60⟶P2⟶P107⟶RC3RC4RC4⟶P31⟶P141⟶P59⟶P68⟶P147⟶P45⟶P151⟶P55⟶P38⟶P123⟶P94⟶RC4RC5RC5⟶P154⟶P26⟶P121⟶P35⟶P176⟶P175⟶P90⟶P179⟶P172⟶RC5Vehicle 4RC4RC4⟶P14⟶P114⟶P37⟶P84⟶P85⟶P134⟶P173⟶P156⟶P181⟶P178⟶P168⟶RC4RC4RC4⟶P63⟶P21⟶P20⟶P145⟶P19⟶P6⟶P89⟶P22⟶P47⟶P161⟶RC4RC4RC4⟶P144⟶P165⟶P8⟶P46⟶P139⟶P142⟶P16⟶P80⟶P155⟶P23⟶RC4Vehicle 5RC5RC5⟶P150⟶P88⟶P152⟶P57⟶P100⟶P125⟶P162⟶P56⟶P5⟶P96⟶P159⟶RC5RC5RC5⟶P108⟶P65⟶P54⟶P81⟶P30⟶P77⟶P75⟶P101⟶P111⟶RC5Vehicle 6RC1RC1⟶P39⟶P48⟶P93⟶P180⟶P131⟶P44⟶P118⟶P177⟶P110⟶RC1RC1RC1⟶P72⟶P99⟶P73⟶P182⟶P132 ⟶P126⟶P49⟶P33⟶P130⟶P43⟶P40⟶RC1Vehicle 7RC5RC5⟶P3⟶P169⟶P140⟶P9⟶P153⟶P167⟶P164⟶P52⟶P171⟶P170⟶P174⟶RC5RC5RC5⟶P105⟶P79⟶P18⟶P12⟶P58⟶P157⟶P163⟶P64⟶P143⟶P27⟶P127⟶RC5Table8 exhibits that seven vehicles are shared on different routes to serve customers in the multidepot reverse logistics network. For example, vehicle 1 is shared in three routes. One route (e.g., RC1 ⟶ P128 ⟶ P138 ⟶ P7 ⟶ P69 ⟶ P78 ⟶ P15 ⟶ P36 ⟶ P42 ⟶ RC1) departs from RC1. Meanwhile, two routes (e.g., RC3 ⟶ P1 ⟶ P122 ⟶ P50 ⟶ P95 ⟶ P109 ⟶ P51 ⟶ P41 ⟶ P129 ⟶ P13 ⟶ P166 ⟶ P158 ⟶ RC3, RC3 ⟶ P98 ⟶ P124 ⟶ P86 ⟶ P112 ⟶ P82 ⟶ P136 ⟶ P28 ⟶ P137 ⟶ P113 ⟶ P104 ⟶ RC3) depart from RC3. Vehicle sharing can improve the resource utilization and operational efficiency of logistics facilities and promote sustainable development among logistics networks. The comparison of results before and after optimization is shown in Table 9 and Figure 8.Table 9
Comparison of the results before and after optimization.
ScenariosTransportation cost ($)Penalty cost ($)Assignment cost ($)Fixed cost ($)Number of semitrailer trucksNumber of vehiclesTotal cost ($)Before optimization46756000—6350—2138025After optimization13071566215063502721373Figure 8
Comparison of the results before and after optimization.In Table9 and Figure 8, the total cost is reduced from $38025 to $21373 in the initial and optimized reverse logistics network on the basis of the RS and TWA strategies in the MRVRPRSTWA. In addition, the TC is decreased from $4675 to $1307. The sum of the penalty cost and AC after optimization is lower than that before optimization through the TWA strategy. The total number of vehicles is decreased to seven when the vehicles are shared among five RCs.


### 6.5. Analysis and Discussion

Four cases are considered to verify the applicability of the proposed method in MRVRPRSTWA. In Case 1, each logistics facility operates independently. In Case 2, the RS strategy is considered among logistics facilities. In Case 3, the TWA strategy is adopted among logistics facilities. In Case 4, the RS and TWA strategies are adopted among logistics facilities. The calculation results of the four scenarios are shown in Table10 and Figure 9.Table 10
Comparison of the optimization results of the four scenarios.
ScenariosThe total cost ($)Violated time (min)Assigned time (min)The number of semitrailer trucksThe number of vehiclesCase 138025333——21Case 227939188—212Case 332616107247217Case 4213738721527Figure 9
Comparison of the total cost, violated time, and TNV.In Table10 and Figure 9, the total cost, violated time, and TNV are significantly reduced in the optimized logistics network with the RS and TWA strategies. TNV in Case 2 is lower than that in Case 1, which shares vehicles among reverse logistics facilities. The violated time is significantly reduced when the TWA strategy is adopted among logistics facilities in Case 3 compared to that of Case 1. The total cost, violated time, and TNV of Case 4 have the minimum value of $21373, 87 min, and 7, respectively. The total cost, violated time, and TNV are compared in Figure 9. Therefore, the proposed method in MRVRPRSTWA can effectively reallocate resources and improve the operational efficiency.


### 6.6. Implications

The RS and TWA strategies are considered to optimize recycling service in the multidepot reverse logistics network. The implications obtained from this work are described as follows:(1)
The RS strategy contributes to the resource utilization of reverse logistics network and promotes the sustainable development of the intelligent logistics system. Vehicles often serve customers on one route in a working period, thereby resulting in low vehicle utilization. Vehicles serve customers on several routes in a working period based on satisfying customer service time window with the adoption of the RS strategy. This strategy shares the resources within and between logistics facilities to obtain the maximum resource utilization in MRVRPRSTWA.(2)
The TWA strategy can enhance the operational efficiency of logistics networks. The candidate time windows of customers with high penalty costs on account of the violation of time windows are assigned to corresponding customers through the TWA strategy to achieve the maximum operational efficiency in MRVRPRSTWA. The TWA strategy is of great theoretical and practical significance to the sustainable development of reverse logistics networks.(3)
The RS and TWA strategies in reverse logistics network can reduce the negative impact of logistics enterprises on the environment, optimize resource allocation, maximize benefits, promote the development of green logistics, and provide related development decisions for logistics enterprises and government departments. The established mathematical model and the hybrid heuristic algorithm can also provide a reliable reference for the applicability and effectiveness of the RS and TWA strategies in addressing MRVRPRSTWA.


## 6.1. Algorithm Comparison

In this section, benchmark instances are performed to test the quality of NSGA-II in MRVRPRSTWA for comparing its results with those received by multiobjective evolutionary algorithm (MOEA) [72] and multiobjective particle swarm optimization (MOPSO) [73]. The multidepot VRPTW (MDVRPTW) datasets obtained from the database of the NEO research group (https://neo.lcc.uma.es/vrp/vrp-instances/multiple-depot-vrp-with-time-windows-instances/) are illustrated in Table 3.Table 3
Relevant settings of data instances.
InstancesDatasetsNumber of depotsNumber of customersNumber of candidate time windowsVehicle capacity1MDPR144822002MDPR249621953MDPR3414431904MDPR4419241855MDPR5424041806MDPR6428851757MDPR767222008MDPR8614431909MDPR96216418010MDPR106288517011MDPR11448220012MDPR12496219513MDPR134144319014MDPR144192418515MDPR154240418016MDPR164288517517MDPR17672220018MDPR186144319019MDPR196216418020MDPR2062885170Table3 presents the relevant characteristics of 20 benchmarks. The parameters for NSGA-II and MOEA are as follows: population size popsize = 150, selection probability sp = 0.6, crossover probability cp = 0.9, and mutation probability mp = 0.1. The parameters for MOPSO are as follows: maximum iteration number max_it = 200, inertia weight ω = 0.9, personal confidences pc = 2, and social learning confidences gc = 3. The results of the three algorithms, which contain total cost (cost), TNV (vehicle), and computation time (Ct), are compared in Table 4.Table 4
Comparison of the results of the three algorithms on different instances.
InstancesNSGA-IIMOEAMOPSOCost ($)VehicleCt (s)Cost ($)VehicleCt (s)Cost ($)VehicleCt (s)1845031841274271891068951792961741951490692211283772043137168210161371021115067921541576410225189401324319923142315158291023121070152381804112244617032112512141215267193381324771695292502037312255173259261820048122542360115248185501026991796810266255261726122484142681019119112782350015287194541128811833532071472092108673321712105985216158401023114772923113127437224191651324217094112331415789102302011014257180531224615159291023721171152411711911251161807412241202601424619231132651715932825121234132792017312271181815810257194111126419367112791919998122622448216298214611327820195371126825779173172074212288Average15479923720019132501752011248t-test−12.955−5.529p-value3.51502E-111.24E-05Table4 exhibits an adequate difference among the optimal results of those three algorithms on the basis of the values of t-test and p-value. The NSGA-II algorithm is superior to MOEA and MOPSO in obtaining the optimal cost in MRVRPRSTWA. The costs of MOEA and MOPSO are higher than that of NSGA-II in 20 benchmark instances. For instance, the average cost of NSGA-II is $15479. However, the costs of MOEA and MOPSO are $20019 and $17520, respectively. The minimum number of vehicles calculated by NSGA-II is nine compared to TNV of MOEA and MOPSO. The NSGA-II outperforms MOEA and MOPSO in terms of the minimum computation time. The average computation time of NSGA-II is 237s, which is lower than the computation times of MOEA and MOPSO. Therefore, the NSGA-II algorithm can achieve better performance in solving MRVRPRSTWA compared to MOEA and MOPSO.


## 6.2. Data Source

A real reverse logistics network in Chongqing, China, is used to verify the applicability and feasibility of the proposed solution in MRVRPRSTWA. The logistics network has five RCs (i.e., RC1, RC2, RC3, RC4, and RC5) and 183 customers. The spatial distribution of logistics network is shown in Figure6. The RCs and customers are marked as stars and rhombuses, respectively. Rhombuses in various colors represent customers served by different RCs. The initial vehicle routes are shown in Table 5.Figure 6
Spatial distribution of RCs and customers.Table 5
Initial vehicle routes before optimization.
FacilitiesRoutesRC1RC1⟶P175⟶P174⟶P167⟶P157⟶P163⟶P172⟶P169⟶P176⟶P162⟶P179⟶P155⟶RC1RC1⟶P180⟶P177⟶P171⟶RC1RC1⟶P150⟶P159⟶P152⟶P165⟶P168⟶P161⟶P170⟶P154⟶P181⟶P158⟶RC1RC1⟶P182⟶P164⟶P156⟶P183⟶P160⟶P153⟶P166⟶P151⟶P178⟶P173⟶RC1RC2RC2⟶P98⟶P107⟶P95⟶P97⟶P99⟶P93⟶P94⟶P84⟶P85⟶P104⟶P103⟶RC2RC2⟶P82⟶P100⟶P89⟶P90⟶P101⟶P96⟶P88⟶P105⟶P87⟶P102⟶P92⟶RC2RC2⟶P91⟶P86⟶P106⟶P83⟶RC2RC3RC3⟶P122⟶121⟶P113⟶P114⟶P115⟶P108⟶P131⟶RC3RC3⟶P124⟶P149⟶P112⟶P147⟶P140⟶P144⟶P141⟶RC3RC3⟶P126⟶P143⟶P120⟶P145⟶P123⟶P125⟶P128⟶P139⟶RC3RC3⟶P148⟶P137⟶P110⟶P138⟶P118⟶P134⟶P136⟶P142⟶P132⟶P129⟶RC3RC3⟶P135⟶P146⟶P133⟶P109⟶P111⟶P116⟶P117⟶P119⟶P130⟶P127⟶RC3RC4RC4⟶P22⟶P16⟶P40⟶P47⟶R42RC4⟶P14⟶P5⟶P9⟶P3⟶P36⟶P33⟶P18⟶P12⟶P30⟶P11⟶P4⟶P35⟶P6⟶RC4RC4⟶P15⟶P7⟶P39⟶P28⟶P37⟶P38⟶P23⟶RC4RC4⟶P2⟶P26⟶P10⟶P24⟶P17⟶P8⟶P45⟶P20⟶P29⟶P19⟶P21⟶RC4RC4⟶P13⟶P1⟶P32⟶P25⟶P34⟶P41⟶P44⟶P43⟶P27⟶P46⟶P31⟶RC4RC5RC5⟶P52⟶P49⟶P58⟶P66⟶P79⟶P57⟶P68⟶P50⟶P53⟶P51⟶P55⟶P71⟶RC5RC5⟶P59⟶P74⟶P76⟶P62⟶P61⟶P70⟶P78⟶P60⟶P80⟶P54⟶P67⟶P56⟶RC5RC5⟶P64⟶P81⟶P77⟶P48⟶P69⟶P73⟶P72⟶P63⟶RC5RC5⟶P75⟶P65⟶RC5In Table5, customers are served by 21 vehicle routes among five RCs. Each RC is operated independently, which results in the waste of resources. For example, three customers (i.e., P180, P177, and P171) belonging to RC1 and two customers (i.e., P75 and P65) belonging to RC5 are served by a vehicle, respectively. The irrational vehicle routes should be decreased to enhance the operational efficiency of the multidepot logistics network.


## 6.3. Relevant Parameter Setting

The objective of MRVRPRSTWA is to obtain the lowest operating cost and TNV in the multidepot reverse logistics network. The related parameters used in the optimization model and the hybrid heuristic algorithm initialization are shown in Table6. The time window of RC1 is [6, 18]. RC2 and RC3 have the same time window [8, 20]. The time window of RC4 and RC5 is the same as [10, 22]. The candidate time windows are constructed on the basis of the real situation of customers that accepted TWA. Hence, four candidate time windows are constructed: [8, 10], [10, 11], [12, 15], and [16, 18].Table 6
Parameter values utilized in the computational experiments.
NotationDefinitionβ1Coefficient of spatial distanceβ2Coefficient of temporal distanceQsMaximum capacity of the semitrailer truckQvMaximum capacity of vehicleUsUsage cost of the semitrailer truck (unit: dollar/km)UvUsage cost of vehiclev (unit: dollar/km)MsAnnual MC of the semitrailer truckMvAnnual MC of the vehicleμePenalty cost for early arrival per unit timeμdPenalty cost for delayed arrival per unit timeμaAC per unit timeWNumber of working days in one planning periodBNumber of planning periods in one yearFL1Fixed cost of RC1FL2Fixed cost of RC2FL3Fixed cost of RC3FL4Fixed cost of RC4FL5Fixed cost of RC5PopsizePopulation sizeNmaxMaximum number of generations runsRmaxMaximum number of runsspSelection probabilitycpCrossover probabilitympMutation probability


## 6.4. Optimization Results

Customers are clustered on the basis of the 3Dk-means clustering algorithm to reduce the computational complexity in the multidepot reverse logistics network. Customers are reassigned to the RC with the nearest spatial and temporal distances. The results of the 3D k-means clustering algorithm are shown in Table 7 and Figure 7.Table 7
3Dk-means clustering results in MRVRPRSTWA.
FacilitiesRC1RC2RC3RC4RC5CustomersP7 P15 P39 P42 P48 P69 P72 P73 P78 P93 P99 P110 P118 P138 P130 P131 P177 P180 P182 P33 P36 P40 P43 P44 P49 P126 P128 P132P11 P17 P24 P25 P32 P34 P53 P91 P92 P97 P102 P133 P135 P146 P148 P160 P183 P115 P116 P117 P119 P120 P61 P62 P70 P66 P67 P71P1 P2 P4 P6 P10 P13 P26 P28 P29 P35 P41 P50 P51 P60 P74 P76 P82 P83 P86 P87 P95 P98 P103 P104 P106 P107 P109 P112 P113 P121 P122 P124 P129 P136 P137 P149 P153 P154 P156 P158 P164 P166P8 P14 P16 P19 P20 P21 P22 P23 P31 P37 P38 P45 P46 P47 P55 P59 P63 P68 P80 P84 P85 P89 P94 P114 P123 P134 P139 P141 P142 P144 P145 P147 P151 P155 P161 P165 P168 P170 P171 P173 P179 P181P3 P5 P9 P12 P18 P27 P30 P52 P54 P56 P57 P58 P64 P65 P75 P77 P79 P81 P88 P90 P96 P100 P101 P105 P108 P111 P125 P127 P140 P143 P150 P152 P157 P159 P162 P163 P167 P169 P172 P174 P175 P176 P178Figure 7
Customer clustering results in MRVRPRSTWA.Table7 illustrates that 183 customers are reassigned to the corresponding logistics facilities and grouped into five clusters. RC1 is allocated 28 customers, RC2 is allocated 28 customers, RC3 is allocated 42 customers, and RC4 and RC5 are allocated 42 and 43 customers, respectively. The spatial and temporal distances of the customers are shown in Figure 7, wherein x represents the longitude of coordinate, y indicates the latitude of coordinate, and z denotes the mid value of the time window. The optimized vehicle routes with the RS and TWA strategies are shown in Table 8.Table 8
Vehicle routes with the RS and TWA strategies.
VehiclesFacilitiesRoutesVehicle 1RC1RC1⟶P128⟶P138⟶P7⟶P69⟶P78⟶P15⟶P36⟶P42⟶RC1RC3RC3⟶P1⟶P122⟶P50⟶P95⟶P109⟶P51⟶P41⟶P129⟶P13⟶P166⟶P158⟶RC3RC3RC3⟶P98⟶P124⟶P86⟶P112⟶P82⟶P136⟶P28⟶P137⟶P113⟶P104⟶RC3Vehicle 2RC2RC2⟶P133⟶P11⟶P160⟶P17⟶P71⟶P115⟶P34⟶P148⟶P120⟶P32⟶P91⟶RC2RC2RC2⟶P102⟶P24⟶P25⟶P146⟶P119⟶P183⟶P70⟶P9⟶P97⟶RC2RC2RC2⟶P61⟶P62⟶P135⟶P116⟶P53⟶P117⟶P66⟶P67⟶RC2Vehicle 3RC3RC3⟶P106⟶P83⟶P87⟶P29⟶P103⟶P10⟶P76⟶P74⟶P4⟶P60⟶P2⟶P107⟶RC3RC4RC4⟶P31⟶P141⟶P59⟶P68⟶P147⟶P45⟶P151⟶P55⟶P38⟶P123⟶P94⟶RC4RC5RC5⟶P154⟶P26⟶P121⟶P35⟶P176⟶P175⟶P90⟶P179⟶P172⟶RC5Vehicle 4RC4RC4⟶P14⟶P114⟶P37⟶P84⟶P85⟶P134⟶P173⟶P156⟶P181⟶P178⟶P168⟶RC4RC4RC4⟶P63⟶P21⟶P20⟶P145⟶P19⟶P6⟶P89⟶P22⟶P47⟶P161⟶RC4RC4RC4⟶P144⟶P165⟶P8⟶P46⟶P139⟶P142⟶P16⟶P80⟶P155⟶P23⟶RC4Vehicle 5RC5RC5⟶P150⟶P88⟶P152⟶P57⟶P100⟶P125⟶P162⟶P56⟶P5⟶P96⟶P159⟶RC5RC5RC5⟶P108⟶P65⟶P54⟶P81⟶P30⟶P77⟶P75⟶P101⟶P111⟶RC5Vehicle 6RC1RC1⟶P39⟶P48⟶P93⟶P180⟶P131⟶P44⟶P118⟶P177⟶P110⟶RC1RC1RC1⟶P72⟶P99⟶P73⟶P182⟶P132 ⟶P126⟶P49⟶P33⟶P130⟶P43⟶P40⟶RC1Vehicle 7RC5RC5⟶P3⟶P169⟶P140⟶P9⟶P153⟶P167⟶P164⟶P52⟶P171⟶P170⟶P174⟶RC5RC5RC5⟶P105⟶P79⟶P18⟶P12⟶P58⟶P157⟶P163⟶P64⟶P143⟶P27⟶P127⟶RC5Table8 exhibits that seven vehicles are shared on different routes to serve customers in the multidepot reverse logistics network. For example, vehicle 1 is shared in three routes. One route (e.g., RC1 ⟶ P128 ⟶ P138 ⟶ P7 ⟶ P69 ⟶ P78 ⟶ P15 ⟶ P36 ⟶ P42 ⟶ RC1) departs from RC1. Meanwhile, two routes (e.g., RC3 ⟶ P1 ⟶ P122 ⟶ P50 ⟶ P95 ⟶ P109 ⟶ P51 ⟶ P41 ⟶ P129 ⟶ P13 ⟶ P166 ⟶ P158 ⟶ RC3, RC3 ⟶ P98 ⟶ P124 ⟶ P86 ⟶ P112 ⟶ P82 ⟶ P136 ⟶ P28 ⟶ P137 ⟶ P113 ⟶ P104 ⟶ RC3) depart from RC3. Vehicle sharing can improve the resource utilization and operational efficiency of logistics facilities and promote sustainable development among logistics networks. The comparison of results before and after optimization is shown in Table 9 and Figure 8.Table 9
Comparison of the results before and after optimization.
ScenariosTransportation cost ($)Penalty cost ($)Assignment cost ($)Fixed cost ($)Number of semitrailer trucksNumber of vehiclesTotal cost ($)Before optimization46756000—6350—2138025After optimization13071566215063502721373Figure 8
Comparison of the results before and after optimization.In Table9 and Figure 8, the total cost is reduced from $38025 to $21373 in the initial and optimized reverse logistics network on the basis of the RS and TWA strategies in the MRVRPRSTWA. In addition, the TC is decreased from $4675 to $1307. The sum of the penalty cost and AC after optimization is lower than that before optimization through the TWA strategy. The total number of vehicles is decreased to seven when the vehicles are shared among five RCs.


## 6.5. Analysis and Discussion

Four cases are considered to verify the applicability of the proposed method in MRVRPRSTWA. In Case 1, each logistics facility operates independently. In Case 2, the RS strategy is considered among logistics facilities. In Case 3, the TWA strategy is adopted among logistics facilities. In Case 4, the RS and TWA strategies are adopted among logistics facilities. The calculation results of the four scenarios are shown in Table10 and Figure 9.Table 10
Comparison of the optimization results of the four scenarios.
ScenariosThe total cost ($)Violated time (min)Assigned time (min)The number of semitrailer trucksThe number of vehiclesCase 138025333——21Case 227939188—212Case 332616107247217Case 4213738721527Figure 9
Comparison of the total cost, violated time, and TNV.In Table10 and Figure 9, the total cost, violated time, and TNV are significantly reduced in the optimized logistics network with the RS and TWA strategies. TNV in Case 2 is lower than that in Case 1, which shares vehicles among reverse logistics facilities. The violated time is significantly reduced when the TWA strategy is adopted among logistics facilities in Case 3 compared to that of Case 1. The total cost, violated time, and TNV of Case 4 have the minimum value of $21373, 87 min, and 7, respectively. The total cost, violated time, and TNV are compared in Figure 9. Therefore, the proposed method in MRVRPRSTWA can effectively reallocate resources and improve the operational efficiency.


## 6.6. Implications

The RS and TWA strategies are considered to optimize recycling service in the multidepot reverse logistics network. The implications obtained from this work are described as follows:(1)
The RS strategy contributes to the resource utilization of reverse logistics network and promotes the sustainable development of the intelligent logistics system. Vehicles often serve customers on one route in a working period, thereby resulting in low vehicle utilization. Vehicles serve customers on several routes in a working period based on satisfying customer service time window with the adoption of the RS strategy. This strategy shares the resources within and between logistics facilities to obtain the maximum resource utilization in MRVRPRSTWA.(2)
The TWA strategy can enhance the operational efficiency of logistics networks. The candidate time windows of customers with high penalty costs on account of the violation of time windows are assigned to corresponding customers through the TWA strategy to achieve the maximum operational efficiency in MRVRPRSTWA. The TWA strategy is of great theoretical and practical significance to the sustainable development of reverse logistics networks.(3)
The RS and TWA strategies in reverse logistics network can reduce the negative impact of logistics enterprises on the environment, optimize resource allocation, maximize benefits, promote the development of green logistics, and provide related development decisions for logistics enterprises and government departments. The established mathematical model and the hybrid heuristic algorithm can also provide a reliable reference for the applicability and effectiveness of the RS and TWA strategies in addressing MRVRPRSTWA.


## 7. Conclusions

This work studies the RS and TWA strategies in multidepot recycling logistics networks. The RS strategy can effectively improve resource utilization and optimize resource configuration in logistics networks. A bi-objective optimization model is proposed for MRVRPRSTWA with the minimum operating cost and TNV. A hybrid heuristic algorithm, which includes 3Dk-means clustering algorithm and NSGA-II algorithm, is also designed to address MRVRPRSTWA. The 3D k-means clustering algorithm is used to cluster customers in the first step of the hybrid heuristic algorithm for reducing the computational complexity of multidepot vehicle routing problem. Next, the NSGA-II algorithm is used to find the optimal vehicle routes and assign candidate time windows to customers.A real-world case study in Chongqing City, China, is analyzed to verify the performance of the proposed mathematical model and the hybrid heuristic algorithm in MRVRPRSTWA. The computational results demonstrate that the total operating cost and TNV can be decreased through the proposed method. The total operating cost is reduced from $38025 to $21373, and TNV is minimized from 21 to 7. The comparison of the performance of three algorithms shows that the NSGA-II has a better performance in solving MRVRPRSTWA than MOEA and MOPSO. Four different scenarios with and without RS and TWA strategies are presented to demonstrate the applicability of the proposed method. The results of four scenarios indicate that the operational efficiency of the logistics network can be improved via the RS and TWA strategies.In this study, the RS and TWA strategies are encouraged to promote the sustainable development of multidepot logistics network. On the one hand, the RS strategy among logistics facilities can optimize the resource configuration and enhance the resource utilization. On the other hand, the TWA strategy can improve the operational efficiency of logistics network by assigning candidate time windows to customers. Although only the RS and TWA strategies can improve the operational efficiency of logistics network to a certain extent, the greatest benefits can be obtained by simultaneously adopting both, which is demonstrated in section6.5. The RS and TWA strategies can be used as references when policy-makers make decisions for logistics development.Many research questions can be studied in the future on the basis of this work. (1) Delivery activity can be considered on the basis of the MRVRPRSTWA to construct pickup and delivery logistics network. (2) The proposed algorithm can be combined with exact algorithms to enhance the accuracy and efficiency of computation. (3) Collaboration among logistics facilities and profit allocation mechanism can be investigated on the MRVRPRSTWA. (4) The environmental impact in the MRVRPRSTWA can be explored in the future.



---

*Source: 2327504-2021-05-17.xml*