A Two-Echelon Electric Vehicle Routing Problem with Time Windows and Battery Swapping Stations

Abstract

Driven by the new laws and regulations concerning the emission of greenhouse gases, it is becoming more and more popular for enterprises to adopt cleaner energy. This research proposes a novel two-echelon vehicle routing problem consisting of mixed vehicles considering battery swapping stations, which includes one depot, multiple satellites with unilateral time windows, and customers with given demands. The fossil fuel-based internal combustion vehicles are employed in the first echelon, while the electric vehicles are used in the second echelon. A mixed integer programming model for this proposed problem is established in which the total cost, including transportation cost, handling cost, fixed cost of two kinds of vehicles, and recharging cost, is minimized. Moreover, based on the variable neighborhood search, a metaheuristic procedure is developed to solve the problem. To validate its effectiveness, extensive numerical experiments are conducted over the randomly generated instances of different sizes. The computational results show that the proposed metaheuristic can produce a good logistics scheme with high efficiency.

1. Introduction

With the development of industrialization and urbanization, the significant increase in greenhouse gas (GHS) emissions caused by overexploitation of a large number of fossil fuel resources has led to increasingly serious air pollution in cities recently. For example, in January 2013, a large-scale haze continued to cover Central and Eastern China [1]. Since the emission was owed to large fossil fuel-based internal combustion vehicles (ICVs), this is identified as the main cause for fog haze weather, and large ICVs with heavy loads are starting to be prevented from entering urban areas by more and more governments in China. Thus, it is quite common now that large ICVs are only permitted to transport demands from the depots located in the suburbs of cities to the satellites near customers, while small vehicles with lower emissions (e.g., electric vehicles) are employed to transport demands from satellites to customers, which arouses the interests of the investigators in the two-echelon vehicle routing problems (2E-VRP) [2].

Due to the environmental advantages, electric vehicles (EVs) have gradually played an important role in transportation. The environmentally friendly features of EVs, such as free of GHG, lower noise pollution, and higher energy efficiency, can help transportation companies to gain a green image or even subsidies for environmental protection. The governments of China, the US, Japan, the EU, and other countries have strongly promoted the use of EVs and encouraged supply chains to apply EVs in constructing a sustainable logistics distribution network [3]. As early in 2008, EVs were introduced in Beijing to serve the transportation of the Olympic Games. Similarly, in Vienna, fully electric vehicles have been operating for several years. In 2017, Shenzhen, China replaced more than 16,300 buses in the city with EVs, making it the city with the largest fleet of EVs in the world [4]. As one of the largest e-commerce companies in China, JD.com (accessed on 1 September 2019) is gradually replacing its existing vehicles with EVs [5]. In 2018, the amount of EV ownership in China was ranked first in the world [6]. Furthermore, EV sharing systems are being considered recently [7]. The routing problem considering EVs has become a research hotspot [8].

In this work, a new two-echelon vehicle routing problem with time windows, including one depot, multiple satellites, and customers, as well as several battery swapping stations (BSSs), noted as 2E-EVRPTW-BSS, was considered. Compared with existing investigations, two new features of the problem can be highlighted. First, considering the limited capacity of EV batteries, insufficient charging installations for long-distance transportation, as well as transportation efficiency, the mixed-fleet scheme is considered. That is, conventional heavy-loaded trucks are used from the depot, which is usually far away from the urban area, to the satellites, which are supposed to be located around the urban fringe, for the first echelon, while EVs are employed from the satellites to the customers within the city for the second echelon due to the environmental considerations or regulations. Secondly, the optimization of the transportation strategies of the two echelons is effectively linked via the time constraint of the satellites as an intermediate, while they are usually optimized separately according to previous research. This is essential to reduce or even avoid the waiting time of the demands at satellites in real life. A hybrid metaheuristic based on the variable neighborhood search (VNS) is developed to solve this new problem, which has been proven to be both effective and efficient by extensive numerical experiments.

The rest of the research is organized as follows. In Section 2, a review of the related literature is introduced. In Section 3, a mixed integer linear programming model is proposed to formulate the problem under consideration. Section 4 details the proposed algorithm for solving this 2E-EVRPTW-BSS, while numerical experiments are carried out and analyzed in Section 5. Finally, Section 6 summarizes the conclusions of this paper.

2. Literature Review

The problem of 2E-EVRPTW-BSS mainly includes two streams of concerns: the electric vehicle routing problem (EVRP) and the two-echelon vehicle routing problem (2E-VRP). Hence, the survey is organized from these two aspects.

2.1. Related Research on the EVRP

Considering the threat of generated pollution in the transportation system, research on the Green Vehicle Routing Problem (GVRP) has become a hot issue [9]. As an extension of the GVRP, the EVRP has been studied for many years, such as by Cru-Chávez et al. [10] and Cattaruzza et al. [11]. In the past ten years, with the popularity of EVs around the world, there has been a myriad of research carried out on the use of EVs in the field of the transportation system [12].

The first to consider optimization algorithms for EVs were Conrad and Figliozzi [13]. In their paper, the batteries of EVs can be recharged at customer locations with extra recharging costs. The research regarding the technical and market background of EV demand distribution, as well as a review of EVRP’s research work were presented by Pelletier et al. [14]. The main studies on EVs in the field of transportation research were summarized in their work. Many advantages provided by EVs were given by them, such as frequent stop-and-go, low driving speed, and reduction in air pollution. Time window constraints for customer delivery and restrictions on recharging stations (RSs) were studied by Schneider et al. [15]. Four variants of the EVRP considering time windows categorized by the allowable charging times of each route and the battery charging scheme (full recharge or partial recharge) were studied by Desaulniers et al. [16]. The BSS location routing problem with capacitated EVs, which aims to determine the locations of BSSs and their routing plan, was proposed by Yang and Sun [17]. An online routing problem regarding EVs was presented by Adler and Mirchandani [18]. Their aim was to make an occasional detour for vehicles by booking a battery to benefit future vehicles, thereby minimizing the average delay of all vehicles. For a given EV transportation network, how to achieve the solution with the shortest transportation time was stated in Liao et al. [19]. From an environmental perspective, Zhang et al. [20] studied EVRP considering RSs to find the minimized energy consumption routing plan. Karakatič [21] considered the nonlinear recharging times at stations and optimized the EVRPTW with multiple depots. Raeesi et al. [22] researched the Electric Vehicle Routing Problem with Time Windows and Synchronized Mobile Battery Swapping (EVRPTW-SMBS), which considered Battery Swapping Vans to recharge EVs in a designated time and space.

A branch-price-and-cut algorithm proposed by Desaulniers et al. [16] aimed to optimally solve four variants of the problem in their work. They also found the optimal solution for instances with up to 100 customers and 21 RSs. However, the EVRP was more commonly solved by heuristics. The electric vehicle routing problem considering time windows and recharging stations (E-VRPTW) was stated by Schneider et al. [15] and the hybrid heuristic combining VNS and Tabu search was considered to solve the E-VRPTW. To solve the BSS location routing problem involving capacitated vehicles, methods on how to extend algorithms to the vehicle routing problems considering intermediate stops were shown in Hof et al. [23]. Their work improved the best solutions previously found for the instances of Yang and Sun [17] by using the proposed extended adaptive variable neighborhood search (AVNS) algorithm. An adaptive large neighborhood search (ALNS) algorithm was introduced by Keskin and Çatay [24] for the EVRP considering partial recharging. Breunig et al. [25] applied EVs in the second echelon of the two-echelon transportation system, which was named as the electric two-echelon vehicle routing problem (E2EVRP).

2.2. Related Research on the 2E-VRP

The two-echelon vehicle routing problem (2E-VRP) can be regarded as a special case of the two-echelon location routing problem (2E-LRP), where the locations of the depots and satellites are usually given. Cuda et al. [26] gave the review of the 2E-VRP and the 2E-LRP. The routing problem was first applied to the two-echelon transportation system by Jacobsen and Madsen [27]. After that, the definition of a 2E-VRP was given by Crainic et al. [28]. The 2E-VRP transportation system was first proposed in Perboli et al. [29]. According to their research, it can be concluded that a mathematical model of the 2E-VRP based on the flow of the demand of each arc was first introduced. Constraints such as time windows of customers and synchronization were considered in the 2E-VRP in Grangier et al. [30]. Li et al. [31] considered the real-time transshipment capacity variation of the 2E-VRP, and in their work, a mathematical model was given. Some exact algorithms were adopted to solve the 2E-VRP, and they can be found in Jepsen et al. [32], Baldacci et al. [33], and Santos et al. [34]. However, most of the 2E-VRPs were solved by heuristics. Crainic et al. [35] introduced a multi-start heuristic that solved the two-echelon transportation subproblems iteratively. Breunig et al. [36] proposed the large neighborhood search (LNS) combined with local search to solve the 2E-VRP and only applying ICVs both in two echelons. An adaptive large neighborhood search (ALNS) algorithm with kinds of destroy and repair operations was proposed by Hemmelmayr et al. [37] to solve the problem. Wang et al. [38] combined the VNS algorithm and the integer programming algorithm to solve the two-echelon capacitated vehicle routing problem with environmental considerations (2E-CVRP-E). Kitjacharoenchai et al. [39] presented a new two-echelon routing problem, which considered a synchronized truck-drone operation by allowing several drones to fly from a truck. Boysen et al. [40] presented and discussed research on the last-mile delivery.

The EVRP and the two-echelon network issues have drawn more and more public attention in academic research and practical application in the literature. In order to determine the urban distribution strategy under the limit of battery range, Jie et al. [41] first studied the 2E-VRP in the context of EV application. Breunig et al. [25] used EVs in the second-echelon transportation, and they observed that the detour miles due to recharging decrease as a function of the charging stations’ density; however, they did not consider the time windows of customers and the unilateral time windows of satellites. Considering that most distribution companies follow the customer-oriented service standard nowadays, how to satisfy the requests of customers, such as time windows, is more important. For the two-echelon transportation system, studying unilateral time windows of satellites can avoid extra wait at satellites for vehicles. Wang et al. [42] proposed the mathematical model of the 2E-VRP in the context of ICV and EV application, which considered the time windows.

Based on the survey, it is not difficult to find that only a limited number of papers with a focus on 2E-VRP involving EVs, where the research on the mixed fleet and integrated optimization for two echelons are rare. Motivated by the state of the art of the investigations and the managerial requirements from the real applications, this paper aims at providing a new solution to this problem.

3. Problem Formulation

3.1. Problem Description

The 2E-EVRPTW-BSS system considers a depot, multiple satellites, a number of customers with deterministic demands, and a group of BSSs that have no restrictions on the capacity of EVs. In addition, the positions of all these nodes have been given and fixed.

Heavy loaded ICVs are employed in the first echelon to transport demands from the depot to each satellite, while EVs are considered in the second echelon to transport demands from satellites to each customer, as depicted in Figure 1. EVs are supposed to have a limited driving range and have to visit BSSs to swap their batteries before they run out of battery power. In addition, to promise the timeliness of delivery, this work also considers the time windows of customers, which in turn require the latest arrival time of ICVs at the satellites, i.e., there is also a time window for the arrival time of ICVs at each satellite. The 2E-EVRPTW-BSS aims to minimize the total costs, including the transportation cost of two echelons, fixed cost of the vehicles, and battery swapping cost.

3.2. Model

The 2E-EVRPTW-BSS is defined on a complete digraph G = (V, A), where the set of nodes V includes customers, satellites, and BBSs. Demands are transported from the depot to the set of customers. Denote q_i as the demand of the ith customer. The number of available ICVs in the first echelon is K₁, while that of EVs in the second echelon is K₂. The capacity of the ICVs and EVs are Q₁ and Q₂, respectively.

Some assumptions are made for the model.

(1)

Direct delivery from the depot to customers is prohibited.

(2)

Each satellite can be visited by ICVs multiple times in the first echelon.

(3)

Only the first echelon allows split deliveries.

(4)

There is no direct travel between satellites in the second echelon.

(5)

The supply at the depot is sufficient.

(6)

The loading and service time at customers are assumed as being included into the travel time, and thus can be ignored.

(7)

The violation of time windows is not allowed.

The variables and parameters for the model are summarized in

Table 1. Variables and parameters of the 2E-EVRPTW-BSS.

Symbol	Description
Sets
$V_0$	$\mathrm{Depot}, V_0=\left\{0\right\}$.
$V_S$	$\mathrm{Set} \mathrm{of} \mathrm{satellites}, V_S=\{1,2,...,m\}$.
$V_C$	$\mathrm{Set} \mathrm{of} \mathrm{customers}, V_C=\{1,2,...,n\}$.
$V_B$	$\mathrm{Set} \mathrm{of} \mathrm{BSSs}, V_B=\{1,2,...,l\}$.
$V_1$	$V_1=V_0\cup V_S$.
$V_2$	$V_2=V_S\cup V_B\cup V_C$.
$V$	$V=V_0\cup V_S\cup V_B\cup V_C$.
K1, K2	Set of the ICVs in the 1st-echelon and set of EVs in the 2nd-echelon.
K	$K=K_1\cup K_2$.
Parmeters
M	A large positive number.
Tmax	The max running time of the proposed VNS algorithm.
$Q_1,Q_2$	The capacity of the ICV and EV, respectively.
$B$	The battery capacity of the EVs.
c1, c2	Cost of unit time from i to j for ICV and EV, respectively.
c3	Battery swapping cost.
$c_f^1,c_f^2$	The fixed cost of ICV and EV.
$\left[e_i,l_i\right]$	The time window for the deriving at node i.
$a$	The electric power consumption of each EV for unit time.
qi	The demand for node i.
tij	Travel time from node i to node j.
Decision and variables
$t_i^1,t_i^2$	Arrival time at the node i in the first echelon and the second echelon, respectively.
ts	Latest arrival time of ICV at the satellite s.
$b_{ik}^{+}$	The remaining battery power of EV k when it arrives at node i in the 2nd-echelon.
$b_{ik}^{-}$	The remaining battery power of EV k when it leaves node i in the 2nd-echelon.
uik	The remaining demand in vehicle k at node i.
xijk	1 if vehicle k travels arc (i,j) in the 1st-echelon; 0 otherwise.
zijsk	1 if EV k travels arc (i,j) from the satellite s; 0 otherwise.

.

In our model, the total cost to be minimized includes three parts:

(1)

Travelling costs of ICVs and EVs:

$$f_1=c_1{\displaystyle \sum _{k\in K_1}{\displaystyle \sum _{i\in V_1\backslash \left\{j\right\}}{\displaystyle \sum _{j\in V_1}{t}_{ij}{x}_{ijk}}}}$$

(1)

$$f_2=c_2{\displaystyle \sum _{k\in K_2}{\displaystyle \sum _{s\in V_S}{\displaystyle \sum _{i\in V_2\backslash \left\{j\right\}}{\displaystyle \sum _{j\in V_2}{t}_{ij}{z}_{ijsk}}}}}$$

(2)

(2)

Battery swapping costs:

$$f_3=c_3{\displaystyle \sum _{k\in K_2}{\displaystyle \sum _{s\in V_S}{\displaystyle \sum _{i\in V_B\backslash \left\{j\right\}}{\displaystyle \sum _{j\in V_2\backslash i}{z}_{ijsk}}}}}$$

(3)

(3)

Fixed cost of ICVs and EVs:

$$f_4=c_f^1({\displaystyle \sum _{k\in K_1}{\displaystyle \sum _{j\in V_S}{x}_{0jk}}})$$

(4)

$$f_5=c_f^2({\displaystyle \sum _{s\in V_S}{\displaystyle \sum _{k\in K_2}{\displaystyle \sum _{j\in V_C\cup V_B}{z}_{sjsk}}}})$$

(5)

The model is presented as follows.

Objective:

$$\mathrm{Min}f=f_1+f_2+f_3+f_4+f_5$$

(6)

Subject to:

$${\displaystyle \sum _{i\in V_1}{x}_{ijk}}={\displaystyle \sum _{i\in V_1}{x}_{jik}} \forall j\in V_1,i\ne j,k\in K_1$$

(7)

$${\displaystyle \sum _{k\in K_1}{\displaystyle \sum _{i\in V_1}{x}_{ijk}}\ge 1} \forall j\in V_S,i\ne j$$

(8)

$${\displaystyle \sum _{i\in V_S}{u}_{ik}}\le Q_1 \forall k\in K_1$$

(9)

$$t_j^1\ge t_i^1+t_{ij}-M(1-x_{ijk}) \forall i\in V_1,j\in V_S,i\ne j,k\in K_1$$

(10)

$$t_j^2\le t_i^2+t_{ij}+M(1-x_{ijk}) \forall i\in V_1,j\in V_S,i\ne j,k\in K_1$$

(11)

$$t_s^1\le t_s \forall s\in V_S$$

(12)

$${\displaystyle \sum _{j\in V_2}{z}_{ijsk}}={\displaystyle \sum _{j\in V_2}{z}_{jisk}} \forall i\in V_2,i\ne j,k\in K_2,s\in V_S$$

(13)

$${\displaystyle \sum _{k\in K_2}{\displaystyle \sum _{s\in V_S}{\displaystyle \sum _{i\in V_2}{z}_{ijsk}}}}=1 \forall j\in V_C,i\ne j$$

(14)

$${\displaystyle \sum _{i\in V_2\backslash \left\{j\right\}}{\displaystyle \sum _{j\in V_2}{q}_j{z}_{ijsk}}}\le Q_2 \forall s\in V_S,k\in K_2$$

(15)

$${\displaystyle \sum _{k\in K_1}{u}_{sk}}\ge {\displaystyle \sum _{k\in K_2}{\displaystyle \sum _{i\in V_2\backslash \left\{j\right\}}{\displaystyle \sum _{j\in V_2}{q}_j{z}_{ijsk}}}} \forall s\in V_S$$

(16)

$$e_i\le t_i^2\le l_i \forall i\in V_C$$

(17)

$$t_j^2\ge t_i^2+t_{ij}-M(1-z_{ijsk}) \forall i\in V_2,j\in V_C\cup V_B,i\ne j,k\in K_2$$

(18)

$$t_j^2\le t_i^2+t_{ij}+M(1-z_{ijsk}) \forall i\in V_2,j\in V_C\cup V_B,i\ne j,k\in K_2$$

(19)

$$b_{ik}^{-}=b_{ik}^{+} \forall i\in V_C,\forall k\in K_2$$

(20)

$$b_{ik}^{-}=B \forall i\in V_S\cup V_B,k\in K_2$$

(21)

$$b_{jk}^{+}\le b_{ik}^{-}-a{t}_{ij}{z}_{ijsk}+B(1-z_{ijsk}) \forall i,j\in V_2,i\ne j,s\in V_S,k\in K_2,$$

(22)

$$b_{jk}^{+}\ge b_{ik}^{-}-a{t}_{ij}{z}_{ijsk}-B(1-z_{ijsk}) \forall i,j\in V_2,i\ne j,s\in V_S,k\in K_2,$$

(23)

$$b_{jk}^{-},b_{jk}^{+}\ge 0 \forall j\in V_2,k\in K_2$$

(24)

$$x_{ijk}\in \{0,1\} \forall i,j\in V_1,i\ne j,k\in K_1$$

(25)

$$z_{ijsk}\in \{0,1\} \forall i,j\in V_2,i\ne j,s\in V_S,k\in K_2$$

(26)

$$u_{ik}\ge 0 \forall i\in V_1, \forall k\in K_1$$

(27)

$$t_i^1,t_j^2\ge 0 \forall i\in V_S,j\in V_2$$

(28)

The main constraints of the above model can be divided into three parts. The constraints of the first echelon and the second echelon are the Formulas (7)–(11) and (13)–(23), respectively, while the constraint (12) links the two echelons through the latest arrival time at each satellite.

More specifically, constraints (7) and (13) ensure the flow conservation for ICVs and EVs, respectively. Constraint (8) implies that each satellite can be visited more than once by ICVs. Constraints (9) and (15) promise that the loading cannot exceed the capacities of ICV and EV, respectively. Constraints (10) and (18) give the upper bound of the arrival time of two echelons. Constraints (11) and (19) specify lower bounds of the arrival time of two echelons. Constraint (12) requires that the arrival time of ICV at each satellite s should be earlier than t_s. Constraint (14) guarantees that each customer can be serviced by an EV only once. Constraint (16) ensures fully meeting the demands of each customer. Constraint (17) promises that each EV should arrive at each customer within the time window. Constraint (20) means that the EV does not consume battery power during the operation of a customer. Constraint (21) guarantees that the EV is fully charged when it starts from the satellite or a BSS. Constraints (22) and (23) keep track of the state of battery power of EVs starting from node i to node j.

4. Algorithms for the 2E-EVRPTW-BSS

Considering the classical VRP is NP hard [36], its variant 2E-EVRPTW-BSS is undoubtedly an NP-hard problem as well. Hence, we adopt a metaheuristic to achieve the optimal solution within a reasonable time.

The proposed metaheuristic is developed based on a VNS algorithm. Mladenović and Hansen [43] first proposed the VNS algorithm. In their work, the systematic changes of the neighborhood were employed within the local search to get rid of the local optimal solutions, which can be used to successfully solve myriad optimization problems in reality [44,45].

In this study, the proposed algorithm begins with an initial solution S₀ = ($S_1^0$,$S_2^0$), where $S_1^0$ and $S_2^0$ are the initial solutions of the first and second echelons, respectively. The shaking neighborhood structures (denoted as N_k, k = 1,2,…, k_max and k_max is the number of shaking neighborhood structures) are first used to improve the efficiency in searching for solution space, followed by the local search neighborhood structures (denoted as the set N_l, l = 1,2,…, l_max, and l_max is the number of local search neighborhood structures). A solution $S\prime =(S_1^{\prime },S_2^{\prime })$ is first generated in a shaking neighborhood structure N_k in each iteration, then a local search procedure is applied to produce a better solution S’’ for improvement.

The iterative process stops until the shaking neighborhoods have been exhausted or the given termination criterion is met—the maximum running time (T_max) in this paper. The overall iterative structure of the proposed VNS is presented in Algorithm 1, and the details will be discussed in the following sub-sections.

Algorithm 1. VNS.

Input: initial solution S0 = ($S_1^0$,$S_2^0$).

A set of neighborhood structures Nk, k = 1,2,…, kmax for shaking phase.

S*←S0

k, l = 1

Output: best feasible solution achieved S*.

1. Whilek ≤ kmax

2.  Shaking: pick a random solution $S_2^{\prime }$ from the kth neighborhood of S*.

Construct $S_1^{\prime }$ based on $S_2^{\prime }$. Update the solution as $S\prime =(S_1^{\prime },S_2^{\prime })$

3.  While l ≤ lmax

4.   Local search $S_2^{\prime \prime }$← the best solution from the lth neighborhood of $S_2^{\prime }$.

5.    Construct$S_1^{\prime \prime }$ based on $S_2^{\prime \prime }$ Update the current solution as $S\prime \prime =(S_1^{\prime \prime },S_2^{\prime \prime })$.

6.   If$f_{eva}(S\prime \prime )\le f_{eva}(S\prime )$then

7.    $S\prime \leftarrow S\prime \prime$

8.    L = 1

9.   Else

10.     L = l + 1

11.   If$f_{eva}(S\prime \prime )\le f_{eva}(S$*) then

12.   S*←S’’

13.   K = 1

14.   Else:

15.   K = k + 1

16.   If the termination condition is met, then return S*.

17.   End while

18. End while

4.1. Solution Representation and Encoding

A full solution is composed of two parts; the first part includes a sequence of satellite nodes that are visited by ICVs that start and end at the depot for the first echelon, while the second part contains a sequence of customer nodes (including BBSs) which start from and end at a given satellite node for the second echelon.

As mentioned before, a full solution to 2E-EVRPTW-BSS is denoted as the set S = {S₁,S₂}, where $S_1=\{r_1^1,r_2^1,...,r_{K_1}^1\}$ is the set of routes for the first echelon and $S_2=\{r_1^2,r_2^2,...,r_{K_2}^2\}$ is the set of routes for the second echelon. For each route in S₁, $r_k^1=\left(0,s_1^k,s_2^k,...,0\right),k=1,2,...,K_1$ represents a sequence of satellites traveled by the ICV k which starts and ends at the depot 0. Similarly, for each route in S₂, $r_k^2=\left(s_k,v_1^k,v_2^k,...,s_k\right),k=1,2,...,K_2$ means a sequence of customers (including BSSs) traveled by the EV k which starts and ends at the satellite s_k. Let S($r_k^1$) denote the set of satellites in route $r_k^1$. In route $r_k^2$, C($r_k^2$) and B($r_k^2$) are denoted as the sets of customers and BSSs, respectively. An example is depicted in Figure 2. The full solution S = (S₁, S₂) for 2E-EVRPTW-BSS includes two parts, with S₁ = {(0,1,3,0), (0,2,0), (0,1,2,0)} and S₂ = {(1,5,16,6,7,1), (1,8,9,10,1), (2,11,12,18,13,2), (3,14,19,15,3)}.

In our study, the optimization for the routes of the two echelons are interrelated. Specifically, the routes of the first echelon are generated based on that of the second echelon. Hence, the solution can be encoded with the routes of only the second echelon, i.e., the encoded solution is composed of only $S_2=\{r_1^2,r_2^2,...,r_{K_2}^2\}$. The full solution will be generated through a decoding process.

4.2. Initial Solution

The initial solution (only the second echelon) is constructed based on Schneider et al. [15]. Firstly, for each satellite, the nearest customer node i₀ is chosen, and all of the customers are sorted in ascending order of the angle between the customer itself, the satellite, and the i₀. Then, the customers are inserted into one route iteratively. BSSs are inserted into the route when a violation of electric power occurs. If a violation of capacity occurs, a new route is created. To consider time window requirements, in each route, a customer u can be inserted between successive vertices i and j only when $e_i\le e_u\le e_j$.

4.3. Solution Decoding

The main task for decoding is to determine the full solution including two echelons. Considering the latest arrival time at each satellite, the first echelon routes can be constructed based on those of the second echelon in this paper. The demand delivered by the EVs from each satellite and the latest arriving time of the ICVs at each satellite can also be derived according to all the second-echelon routes departing from and returning to it. Based on the simplified least-cost insertion heuristic, the generation of the first-echelon routes can be described as follows [25].

Step 1. According to the second-echelon routes, determine the latest arrival time of ICVs (t_s) and the required demand for each satellite;

Step 2. Back-and-forth routes are created for each satellite whose demand is larger than capacity Q₁ (delivery split). That is, if the demand of satellite s is larger than Q₁, construct a back-and-forth route (0, i, 0), and this process will continue until the remaining demand of the satellite is smaller than Q₁;

Step 3. After Step 2, all satellites visited by the second-echelon routes are sorted in ascending order of the latest arriving time t_s;

Step 4. Select an ICV, insert sequentially the visited satellites to generate its route until the arrival time of the ICV is later than t_s of the last inserted satellite or the load of the ICV exceeds Q₁. Then, select the next ICV to repeat the inserting process until all the satellites are visited. The first-echelon routes are thus constructed, and a complete solution is derived (as indicated in Figure 2).

4.4. Solution Evaluation

A solution is evaluated by Formula (29):

$$f_{eva}\left(S\right)=f\left(S\right)+\alpha p_{cap}\left(S\right)+\beta p_{tw}\left(S\right)+\gamma p_{bat}\left(S\right)$$

(29)

where f(S) is the total cost as represented by Formula (6), p_cap(S) shows the violation of total capacity defined by (30), p_tw(S) shows the time window violation at customers as indicated in (31), p_bat(S) measures the battery capacity violation defined by (33), and $\alpha ,\beta ,\gamma$ are penalty factors of the violations. The penalty factors are large enough values, which can effectively eliminate the possibility of the solution containing time window violation being the optimal solution.

The total capacity violation of a solution S can be calculated by

$$p_{cap}(S)={\displaystyle \sum _{k=1}^{K_2}\max ({\displaystyle \sum _{i\in C(r_k^2)}{q}_i-Q_2},0})$$

(30)

There are three cases leading to the violation of time windows, i.e., arrival time $t_s^1$ > t_s, $t_i^2$ < e_i and $t_i^2$ > l_i (I ∈ V_C). Hence, the total time window violation of a solution can be calculated by

$$p_{wt}(S)={\displaystyle \sum _{k=1}^{K_1}{\displaystyle \sum _{s\in S(r_k^1)}\max (t_s^1-t_s,0)}}+{\displaystyle \sum _{k=1}^{K_2}{\displaystyle \sum _{i\in C(r_k^2)}\max (t_i^2-l_i},0})+{\displaystyle \sum _{k=1}^{K_2}{\displaystyle \sum _{i\in C(r_k^2)}\max (e_i-t_i^2},0})$$

(31)

To calculate battery capacity violations, for each route $r_k^2\in S_2,k=1,2,...,K_2$, we define the variable ${\mathsf{\Phi }}_{v_i}^k$, k = 1,2,…, K₂ to represent the remaining battery power for the EV k when leaving node v_i.

$${\mathsf{\Phi }}_{v_i}^k=\{\begin{array}{c}B\\ {\mathsf{\Phi }}_{v_{i-1}}^k-a\cdot t_{v_{i-1}{v}_i}\end{array}\begin{array}{c}{v}_{i-1}\in B(r_k^2)\cup \left\{v_0^2\right\}\\ otherwise\end{array}\begin{array}{c}i=1,\dots ,m\end{array};k=1,...,K_2$$

(32)

Obviously, ${\mathsf{\Phi }}_{v_i}^k$ should be non-negative. The battery capacity violation of a route $r_k^2$ can thus be defined as

$$p_{bat}(S)={\displaystyle \sum _{k=0}^{K_2}{\displaystyle \sum _{i\in C(r_k^2)}\max (-{\mathsf{\Phi }}_i^k},0)}$$

(33)

A varying penalty scheme is adopted to accelerate the search process, which is similar to the work of Vidal et al. [46]. When an infeasible solution still violates the constraints after the penalty, each penalty factor will be multiplied by 10 and the local search will be conducted again. If the infeasible solutions still exist after two rounds of multiplication by 10, no further attempts will be made. Therefore, this procedure cannot ensure the elimination of all the infeasible solutions, but it can effectively reduce the number of them at a low cost.

4.5. Neighborhoods and Local Search

The proposed VNS of this paper adopted several frequently used inter- and intra-route neighborhood structures.

2-Opt*. As a modification of the 2-opt operator originally proposed by Lin [47], the 2-Opt* operator was applied to VRPTW by Potvin and Rousseau [48]. For the two routes that depart from the same satellite, this operator replaces the arc (i, I + 1) of one route and arc (j, j + 1) of another route with arcs (i, j + 1) and (j, I + 1) (Figure 3a). With respect to 2E-EVRPTW-BSS, there is at least one BBS among the nodes i, i + 1, j, j + 1 and no two successive BBSs in the generated routes, compared with the 2-Opt* operator in traditional problems in which only customers are swapped.

Relocate. This operator deletes and replaces a node from one route and then inserts it into another route (Figure 3b) or another location in the same route (Figure 3c). It is performed only on the sequence of customers.

Exchange. The exchange operator swaps the positions of two nodes, which is applied to customers. This operator can be used for both intra- (Figure 3d) and inter-route (Figure 3e) moves.

Swap (${\lambda }_1,{\lambda }_2$). This operator swaps ${\lambda }_1$ consecutive customers in one route with ${\lambda }_2$ consecutive customers in another route (Figure 3f), or just swaps in the same route (Figure 3g).

Shift $\left({\lambda }_1,0\right)$. This operator removes ${\lambda }_1$ consecutive customers and inserts them into another position in the same route (Figure 3h) or another route (Figure 3i).

Satellite-change. This operator was proposed by Wang et al. [38]. Subject to the constraints on the second-echelon EVs per satellite, this operator exchanges the current satellite with another satellite and inserts the old satellite into a position of the route (Figure 3j).

StationInRe. This operator, proposed by Schneider et al. [15], performs insertions and removals of BSS. This operator is applied to the arcs starting from a BBS (k, j), $k\in V_B$. If the arc (k, j) does not belong to the current solution, the stationInRe will insert the BBS k between j and its predecessor, as depicted in Figure 3k. If the arc already exists, the BSS k will be removed, as shown in Figure 3l.

The proposed VNS procedure in Algorithm 1 contains two phases to generate a neighborhood, i.e., shaking and local search. In the shaking phase, five inter-route neighborhood structures (2-Opt*, inter-Relocate, inter-Exchange, inter-Swap, and inter-Shift) and one intra-route neighborhood structure (Satellite-Change) are sequentially used to generate the neighborhood.

In the local search phase, five intra-route neighborhood structures, i.e., intra-Relocate, intra-Exchange, intra-Swap, intra-Shift, and Station-InRe, are sequentially applied. Local search neighborhoods are performed in a cyclic manner [49]. Each neighborhood is explored until no further improvement can be found, then the next neighborhood structure is selected. When all the five neighborhood structures have been performed, the search process repeats again. The procedure terminates when the local optimum is reached in every single neighborhood.

5. Numerical Experiments

5.1. Experiment Description and Parameter Settings

Considering the 2E-EVRPTW-BSS is a new model, there are currently no benchmark instances available to evaluate the algorithm performance of this problem; therefore, the test instances are proposed based on the instances of Breunig et al. [25], which are designed for the 2E-VRP. The locations of the depot, satellites and customers are directly used from Breunig et al. [25], and for each instance, several customer nodes are chosen to be the BSSs [50]. Customer demands are randomly generated between 10 and 150 kg.

This paper generates the locations of the BSSs randomly, and the number of BSSs is set to be 1/5 of the nodes in the network according to Schneider et al. [15]. The adjacency matrix of transportation time between any pair of nodes is calculated by dividing Euclidean distance by the ICV speed in the first echelon or the EV speed in the second echelon. The relative parameters are set as shown in

Table 2. Parameter settings for 2E-EVRPTW-BSS model.

Parameter	Value (Unit)	Parameter	Value (Unit)
c1	3	Q1	400 (kg)
c2	5	Q2	280 (kg)
c3	5	v1	50 (km/h)
$c_f^1$	50	v2	40 (km/h)
$c_f^2$	80	$\alpha$	500
a	3 (kWh)	$\beta$	500
B	100 (kWh)	$\gamma$	500

.

5.2. Computational Results and Analysis

The proposed VNS is coded in Python 3.6 and is implemented on an Intel Xeon E3-123v3 (8 GB RAM) computer. The optimal solutions to some small-sized instances can be obtained by using the GUROBI 8.0.1 solver.

5.2.1. Comparison between the Proposed VNS and GUROBI for the Small-Sized Instances

The results obtained by the GUROBI solver and the proposed VNS are presented in

Table 3. Results for the small-sized instances.

Inst	Instance Size			GUROBI		VNS			Gap (VNS/G)
	m	n	l	Obj	Time (s)	Best	Avg	Time (s)
1-1	1	5	1	239.38	9.04	239.38	239.38	0.01	0.00%
1-2	1	5	1	236.29	10.07	236.29	236.29	0.01	0.00%
1-3	1	5	1	233.68	13.01	233.68	233.68	0.01	0.00%
1-4	1	5	1	234.66	10.67	234.66	234.66	0.01	0.00%
1-5	1	5	1	229.40	12.44	229.40	229.40	0.01	0.00%
1-6	2	10	2	469.55	2736.22	469.55	469.55	3.06	0.00%
1-7	2	10	2	476.69	987.44	476.69	476.69	3.04	0.00%
1-8	2	10	2	457.64	3827.14	457.64	457.64	3.07	0.00%
1-9	2	10	2	409.42	2736.12	409.42	409.42	4.04	0.00%
1-10	2	10	2	472.30	3412.34	472.30	472.30	4.04	0.00%
1-11	2	15	3	735.66	4500	733.07	733.07	17.15	−0.35%
1-12	2	15	3	750.79	4500	748.79	748.79	18.36	−0.27%
1-13	2	15	3	822.44	4500	821.67	821.67	19.77	−0.09%
1-14	2	15	3	620.76	4500	620.76	620.76	18.34	0.00%
1-15	2	15	3	841.86	4500	841.86	841.86	26.06	0.00%

Note: the optimal values are marked in bold. Gap (x/y) = 100% × (Obj (x) − Obj (y))/Obj (y).

. The instance size is based on the number of satellites (m), customers (n), and BSSs (l), as well as the number of ICVs and EVs, denoted as n (ICV) and n (EV), respectively. Each instance is operated 10 times, and “Best” is the best-found solution and “Avg” is the average optimal solution for the 10 instances. The objective value and CPU time for GUROBI and VNS are listed in the column of Obj and Time, respectively.

It is easy to find in Table 3 that the proposed VNS achieves the optimal solution for all the 10 instances whose optimum can be reached by GUROBI. However, the running time of the VNS is incredibly short compared with that of GUROBI. With the size of the instance growing, the GUROBI solver fails to produce the exact solutions for the last five instances even when the running time is up to 4500 s, while the proposed VNS algorithm can still find solutions with good performance efficiently. In fact, the proposed VNS can arrive at the solution within 5 s for most small-sized instances and consumes only 26 s even for the instances that GUROBI is hard to handle in a reasonable time.

5.2.2. Comparison between Two Metaheuristics for the Middle- and Large-Sized Instances

In order to evaluate the performance of the proposed VNS for problems with a larger size, this paper compares the proposed VNS with another metaheuristic presented by Wang et al. [38]. Their algorithm does not consider the BBS and time windows; hence, some adaptions have to be made to fit the requirement for solving the 2E-EVRPTW-BSS. The algorithm is denoted as Wang (2017) in

Table 4. Results for the middle-sized instances.

Inst	Instance Size			Wang (2017)			VNS			Gap (VNS/Wang (2017))
	m	n	l	Best	Avg	Time (s)	Best	Avg	Time (s)
2-1	2	20	4	857.05	857.05	64.77	857.05	857.05	73.45	0.00%
2-2	2	20	4	986.73	986.73	78.05	986.73	986.73	80.19	0.00%
2-3	2	20	4	853.67	853.67	100.40	853.67	853.67	123.35	0.00%
2-4	2	30	6	1310.07	1310.07	987.17	1310.07	1310.07	1056.69	0.00%
2-5	2	30	6	1415.36	1415.36	1209.68	1415.36	1415.36	1233.37	0.00%
2-6	2	30	6	1440.31	1440.31	1444.67	1440.31	1440.31	1425.45	0.00%
2-7	2	40	8	1915.31	1915.31	1433.08	1915.31	1915.31	1433.29	0.00%
2-8	2	40	8	1915.59	1915.59	2004.07	1919.59	1915.59	1904.75	0.21%
2-9	2	40	8	1918.16	1918.16	2012.39	1918.16	1918.16	1983.08	0.00%
2-10	2	50	10	2309.42	2309.42	2100.89	2307.70	2307.70	2142.66	−0.07%
2-11	2	50	10	2376.50	2376.50	2387.79	2377.84	2322.84	2223.04	0.06%
2-12	2	50	10	67,586.62	67,586.62	2366.45	60,314.84	60,314.84	2245.89	−10.76%
2-13	3	60	12	81,753.85	81,753.85	3104.56	77,195.14	77,195.14	3004.35	−5.58%
2-14	3	60	12	2818.74	2818.74	3399.97	2802.56	2802.56	3179.97	−0.57%
2-15	3	60	12	2783.41	2783.41	3374.67	2761.43	2761.43	3113.37	−0.79%

and

Table 5. Results for the large-sized instances.

Inst	Instance Size			Wang (2017)		VNS		Gap (VNS/Wang (2017))
	m	n	l	Best	Avg	Best	Avg
3-1	5	100	20	4816.65	4816.65	4788.05	4788.05	−0.60%
3-2	5	100	20	4529.14	4529.14	4522.59	4522.59	−0.14%
3-3	5	100	20	4592.02	4592.02	4457.66	4457.66	−3.01%
3-4	10	100	20	4640.24	4640.24	4513.64	4513.64	−2.80%
3-5	10	100	20	4452.50	4452.50	4352.32	4352.32	−2.30%
3-6	10	100	20	4616.22	4616.22	4555.90	4555.90	−1.32%
3-7	10	100	20	4441.32	4441.32	4382.85	4382.85	−1.33%
3-8	10	100	20	4497.03	4497.03	4432.90	4432.90	−1.45%
3-9	10	100	20	4710.86	4710.86	4589.38	4589.38	−2.65%
3-10	10	100	20	4769.99	4769.99	4695.48	4695.48	−1.59%
3-11	10	100	20	4313.99	4313.99	4247.62	4247.62	−1.56%
3-12	10	100	20	4260.96	4260.96	4251.71	4251.71	−0.22%
3-13	10	100	20	4630.31	4630.31	4562.67	4562.67	−1.48%
3-14	10	100	20	4411.85	4411.85	4348.67	4348.67	−1.45%
3-15	10	200	20	9407.12	9407.12	9277.62	9277.62	−1.40%
3-16	10	200	40	9216.41	9216.41	9142.42	9142.42	−0.81%
3-17	10	200	40	9984.68	9984.68	9851.93	9851.93	−1.35%

.

The results of the middle-sized instances and the large-sized instances are reported in Table 4 and Table 5, respectively. For the middle-sized instances, the algorithm terminates when all the shaking neighborhoods have been exhausted. The objective value and CPU time for the two algorithms, as well as the gaps between the two algorithms, are listed in Table 4. For the large-sized instances, however, we set the max running time T_max of 2000 s as the termination criterion due to the over-long time for exploration. Hence, only the objective values and gaps are listed in Table 5.

We can find from Table 4 that the performance of the two algorithms is quite similar. Although the objective values obtained by the proposed VNS is slightly worse than that of Wang (2017), it exhibits the advantage in running time when the instance size is larger. Namely, the proposed VNS is more efficient in the computation for larger instances. This can be explained that the proposed VNS tries to optimize the whole solution including two echelons, while Wang (2017) optimizes only the sequence of the second-echelon and generates a feasible sequence for the first echelon based on the second part. Hence, the proposed VNS can find a better solution more easily with fewer explorations in local research.

The advantage of the proposed VNS in computation efficiency is even more significant for the large-sized instances. It is fairly apparent in Table 5 that the proposed VNS can achieve a better solution than Wang (2017) with the same running time for almost all the instances (with only one exception). Therefore, the proposed VNS shows a good performance in solution quality and significant advantages in computation efficiency for the 2E-EVRPTW-BSS with a larger size, which has special value in practical applications.

5.2.3. Comparison of Carbon Emission between ICVs and EVs in the Second-Echelon

One of the important advantages of EVs is the reduction in carbon emissions. To illustrate the effects of emission reduction, we compare the emissions of our solution with that of fossil fuel-based internal combustion vehicles (ICVs) for the second echelon. The unit emissions of EVs are calculated based on the formulation of Zhang et al. (2018) [20], and the unit emissions of ICVs are calculated based on the formulation of Wang et al. (2017) [38]. The parameters of an EV are adopted from Pelletier et al. (2019) [51], and the parameters of an ICV are adopted from Wang et al. (2017) [38]. We assume that the road grade angle is 0, the vehicle acceleration is 0, and the EV and the ICV have the same travel speed.

The results are shown in

Table 6. Results of emissions of ICVs and EVs in the second echelon.

Inst	Instance Size		EV	ICV	Gap (EV/ICV)
	m	n
1–6	2	10	86.21	107.12	−19.52%
1–10	2	10	58.47	82.66	−29.26%
2–10	2	50	559.22	694.88	−19.52%
2–11	2	50	513.77	638.41	−19.52%
3–1	5	100	752.01	934.45	−19.52%

, where “EV” means the emissions when only EVs are adopted and “ICV” means the emissions when only ICVs are adopted.

From last column in Table 6, it can be concluded that using EVs in transportation systems rather than ICVs is more environmental. The emissions of EVs can be reduced by about 20%. Thus, compared to ICVs, applying EVs is more environmentally friendly.

6. Conclusions

In this research, a novel two-echelon vehicle routing problem including mixed vehicles considering battery swapping stations (2E-EVRPTW-BSS) is considered, in which the integration of optimization for the two echelons is implemented by restricting the latest arrival time of ICVs at satellites in the first echelon to meet the requirement of EVs in the second echelon. The objective of the 2E-VERPTW-BSS aims to minimize the total cost. Various constraints related to vehicle capacity, battery capacity, time windows, and battery swapping stations are considered based on the requirements from real operations.

As a variant of the VRP, this problem is obviously NP hard and untractable by using exact algorithms. In addition, the integration optimization for a two-echelon transportation system inevitably increases the difficulty in solving this problem. Hence, a kind of VNS algorithm is proposed for solving the 2E-EVRPTW-BSS problem in this paper. To evaluate the performance of the proposed VNS, the algorithm developed by Wang et al. [37] is adapted according to the requirements of 2E-VERPTW-BSS for comparison. The extensive numerical experiments on different sizes of instances are conducted to examine the effectiveness and efficiency of the proposed VNS in this paper. It is found that the proposed VNS can easily achieve the optimal solutions for the small-sized instances. For larger instances, the proposed VNS demonstrates a higher efficiency in searching for the solutions of high performance compared with the algorithm of Wang et al. [37] and can find better solutions in the same running time.

The research of this paper is a kind of extension to the existing investigations. Of course, there are some other factors from the real demands worthy of being considered. For example, uncertainties are almost everywhere in the real world. Hence, in future studies, uncertainties in elements such as transportation time, electric power consumption, etc., can be included in the models, and the optimization under uncertainty can be considered.