A branch-and-price algorithm for two-echelon electric vehicle routing problem

Motivated by express and e-commerce companies’ distribution practices, we study a two-echelon electric vehicle routing problem. In this problem, fuel-powered vehicles are used to transport goods from a depot to intermediate facilities (satellites) in the first echelon, whereas electric vehicles, which have limited driving ranges and need to be recharged at recharging stations, are used to transfer goods from the satellites to customers in the second echelon. We model the problem as an arc flow model and decompose the model into a master problem and pricing subproblem. We propose a branch-and-price algorithm to solve it. We use column generation to solve the restricted master problem to provide lower bounds. By enumerating all the subsets of the satellites, we generate feasible columns by solving the elementary shortest path problem with resource constraints in the first echelon. Then, we design a bidirectional labeling algorithm to generate feasible routes in the second echelon. Comparing the performance of our proposed algorithm with that of CPLEX in solving a set of small-sized instances, we demonstrate the former’s effectiveness. We further assess our algorithm in solving two sets of larger scale instances. We also examine the impacts of some model parameters on the solution.


Introduction
With the depletion of fossil energy and the urgent need for environmental protection, increasing nations have used electric vehicles (EVs) to replace fuel-powered vehicles for urban transport and logistics. For example, the European Commission's Urban Mobility Package aims to achieve zero emissions of urban logistics by 2030 [20], and its partners agree on a green deal about zero-emission urban logistics in the Netherlands [24]. This means that the use of gasoline-or diesel-powered freight vehicles will be limited in many cities. Consequently, the sales of new energy vehicles have registered phenomenal growth in key markets such as Europe, the Middle East, China, and Australia [31]. The cumulative sales of new energy vehicles in China reached 2.7 million in 2018 [41]. In this environment, many logistics (e.g., SF express) and e-commerce companies (e.g., JD.com, Amazon) have used EVs for logistics and distribution in cities B Zhiguo Wu 17113144@bjtu.edu.cn 1 School of Economics and Management, Beijing Jiaotong University, Beijing, China [51]. Different from gasoline-and diesel-powered vehicles, EVs have limited driving ranges and need en-route recharging, which rely on the volume of batteries and the limited recharging infrastructure because of the immature battery technology. The adoption of EVs gives rise to new challenges for the vehicle scheduling and routing problem, which motivates research on tailoring models to capture the limitations of EVs. Specifically, EVs need to decide when and which recharging stations to recharge when they are running out of battery power. As a result, many studies on the electric vehicle routing problem (EVRP) have recently appeared in the literature [21,26,36,39]. The majority of such works on the EVRP only consider a single echelon where the distribution center delivers goods to the customers directly.
In contrast, the two-echelon (or multi-echelon) distribution system is very efficient for city logistics and e-commerce logistics [13,15]. For a B2C e-commerce firm, customers place orders online on its platform, and then, the orders are picked and collected in the depot, which are generally located on the outskirts of cities. The packages are transported to satellites that are located close to the city centers. After that, the packages are assigned to distributors to deliver to customers. Such a two-echelon distribution system can keep Fig. 1 An instance of the 2E-EVRP large fossil fuel-powered vehicles out of the city centers and use EVs with lower pollution to transfer the packages from the satellites to the customers. Motivated by this real practice of transport logistics, we study the two-echelon electric vehicle routing problem (2E-EVRP) to help express and ecommerce companies address their key transport planning and scheduling problem. In the first echelon, fossil fuelpowered trucks transport goods from the depot to a subset of satellites. The second echelon consists of transferring goods from the satellites to the final customers using EVs. Figure 1 shows an instance with one depot, four satellites, seven customers, and three recharging stations. Our model for the 2E-EVRP is to craft an optimal transport plan to minimize the total cost. We use the arc flow formulation to model the problem and design an exact algorithm to solve it.
We make the following contributions in this study: (1) Motivated by the distribution systems adopted by B2C firms and express firms, we propose the 2E-EVRP and build an arc flow model. (2) We develop an efficient branch-and-price (BP) algorithm to solve the 2E-EVRP. Specifically, we decompose the arc flow model as an integer master problem and apply relaxations of the integer master problem to derive lower bounds for the optimal solution. We design a bidirectional labeling algorithm to generate feasible routes with recharging stations in the second echelon. We adopt the state-space relaxation and the ng-route relaxation as acceleration strategies to reduce the time required to solve the subproblem.
(3) Comparing the performance of the BP algorithm with that of CPLEX in solving a set of small-scale instances, we demonstrate the former's effectiveness. We further assess the performance of the BP algorithm in solving two sets of larger scale instances. We also examine the impacts of some model parameters on the solution and cost.
The rest of the paper is organized as follows: We give a review to identify the research gap in Sect. 2. We formulate the 2E-EVRP as an arc flow model and decompose the model as an integer master problem and apply relaxation to derive the lower bounds in Sect. 3. We propose a BP algorithm with accelerate strategies in Sect. 4. We present the results of numerical studies in Sect. 5. Finally, we conclude the paper and give suggestions for future research in Sect. 6.

Literature review
There are many variants of VRP [6,19,25]. However, there are three streams of literature related to our work, namely the two-echelon vehicle routing problem (2E-VRP), the electric vehicle routing problem (EVRP), and the two-echelon electric vehicle routing problem (2E-EVRP).

Two-echelon vehicle routing problem
Crainic et al. [12] considered the 2E-VRP in early time, which is formally introduced by Perboli et al. and Crainic et al. [11,45]. Cuda et al. [14] gave an extensive review of the problem. In the following, we only review the works that are most closely related to our study.
There are some studies that provide exact algorithms for the problem. Perboli and Tadei [43] proposed two metaheuristics and some new families of valid inequalities derived from the traveling salesman problem (TSP) and the vehicle routing problem (VRP). Jepsen et al. [33] presented an edge flow model and applied a relaxation, which is similar to the formulation introduced by Perboli et al. [45]. Based on this, they proposed a branch-and-cut algorithm to solve the problem. Baldacci et al. [5] built a new model of the 2E-VRP and derived some new lower bounds. Based on this, they proposed an exact method that transforms the 2E-VRP into multi-depot capacitated vehicle routing problems (MDCVRPs) with a limited set. Dellaert [15] proposed an arc-based formulation and developed a branch-and-price algorithm that can solve instances of up to five satellites and 100 customers. Perboli et al. [44] further gave two groups of valid inequalities: one was an extension of the existing VRP valid inequalities and the other consisted of explicit valid inequalities for the 2E-VRP.
Because of the hard-solving of the problem, researchers have developed heuristics algorithms for it. Hemmelmayr et al. [29] proposed an adaptive large neighborhood search (ALNS) heuristic for the 2E-VRP. Breunig et al. [8] presented a meta-heuristic, in which some tailored operators were designed to optimize the selection of satellites. Grangier et al. [28] considered integrated constraints (e.g., time windows, synchronization, multiple trips) in the second echelon and proposed an ALNS to solve the problem. Li et al. [37] adopted Clarke and Wright's savings heuristic algorithm, and further improved through a local search phase.
The above works assume that the transport in both echelons performed by gasoline-or diesel-powered vehicles, and do not consider the vehicle en-route recharging. In contrast, we consider the 2E-EVRP, where EVs are used in the second echelon, so we need to consider when and where to recharge the EVs.

Electric vehicle routing problem
Since 2000, green supply chain [1,17,18,52] and city logistics [27,38,48] have attracted extensive interest, and the use of EVs is thought of as an effective means for realizing green distribution. There has been a growing interest in the EVRP, since Gonçalves et al. [26] first studied it. Pelletier et al. [42] gave a comprehensive review and suggested some research perspectives for the EVRP. In the following, we only review the works that are most closely related to our study.
Liao et al. [39] considered the battery charging or swapping in the EVRP, proposed a dynamic programming to solve the EVRP with the shortest travel time, and developed an algorithm for the fixed tour EV touring problem. Hiermann et al. [30] studied a mixed vehicle routing problem that involves various types of EVs with different acquisition costs, driving ranges, and payloads. The results showed that the type and number of EVs to be used depend on the typical customer distribution in the urban area. Desaulniers et al. [16] considered four variants of the EVRP and designed customized branch-and-price algorithms for each variant. Montoya et al. [40] studied the EVRP with piecewise linear charging functions and proposed a hybrid meta-heuristic to solve it. Zhang et al. [53] introduced a new model considering the energy consumption caused by coal-fired generation. Keskin and Çatay [35] compared the total recharging costs under three recharging configurations, i.e., normal, fast, and super-fast, where partial recharging was permitted. Froger et al. [22] further studied the EVRP with a general nonlinear charging function, assuming that the EVs can be fully or partially recharged. Cortés-Murcia et al. [10] proposed a new variant that allows at most one customer to visit when an EV is recharging en-route. Hua and Xu [23] studied the location of charging stations to serve all EVs in a given area. Tahami et al. [49] provided a compact formulation of the EVPR, and three exact methods were introduced to solve the formulation and its augmented variants. The above works focus on the single-echelon EVRP. In contrast, we consider the 2E-EVRP in which EVs are used in the second echelon, while fossil fuel-powered trucks are used in the first echelon.

Two-echelon electric vehicle routing problem
To our best knowledge, there are only three studies on the 2E-EVRP up to now [7,34,50]. Jie et al. [34] studied the 2E-EVRP in which the battery of EV can be replaced at the stations. They proposed a hybrid algorithm combining the column generation (CG) with an ALNS to deal with the problem. Wang et al. [50] studied the 2E-EVPR with time windows of satellites and customers. The randomly generated instances were solved by GUROBI. In our paper, we focus on the exact algorithm to find the optimal solution of 2E-EVRP.
Breunig et al. [7] studied the 2E-EVRP in which EVs were only used in the second echelon, and proposed an exact algorithm for it based on the algorithm proposed by Baldacci et al. [5] for 2EVRP. The algorithm consists of two main steps: (i) generate the set of all first echelon routes and calculate an initial lower bound; (ii)fix the selected set of first echelon routes and transform the resulting problem into MDCVRP, which is further solved by the procedure proposed in Baldacci and Mingozzi [2]. Different from them, we propose an arc flow model and decompose the model into an integer master problem and a pricing subproblem. We then design a branch-and-price algorithm to solve the problem.

Problem formulation
In this section, we formulate the 2E-EVRP as an arc flow model and decompose the model as an integer master problem. The relaxation of the integer master problem is used to derive lower bounds for the optimal solution.

Arc flow model
We model the 2E-EVRP on graph G = {V , A}. The vertex set V = {0} ∪ S ∪ C ∪ F, where 0 represents the depot, S and C represents the sets of satellites and customers, respectively, and F is the set of recharging stations at second echelon. The arc set A = A 1 ∪ A 2 consists of two subsets, namely the set of arcs at first echelon A 1 = {(i, j)|i, j ∈ {0} ∪ S, i = j} and the set of arcs at second echelon Suppose the trucks and the EVs are located at the depot and satellites, respectively. The trucks are homogeneous and their load capacity is Q 1 . The EVs are also homogeneous, and their load capacity and battery capacity are Q 2 and Q 3 , respectively. When an EV travels on arc (i, j), it incurs transport cost c i j and electricity consumption e i j . We assume that c i j and e i j are linear functions of the distance d i j , i.e., c i j = αd i j and e i j = βd i j , α, β > 0. Let E be the maximum recharging time that is required to recharge Q 3 units of electricity, we have E = γ Q 3 , where γ > 0 is the charging rate. Here, we assume that α = β = γ = 1 to simplify the presentation and analysis. In fact, the values of these parameters do not affect the effectiveness of our proposed algorithm.
Let the fixed costs of a truck and an EV be u 1 and u 2 , respectively. Suppose the battery is fully recharged on a visit to a recharging station and an EV can visit a recharging station many times. However, consecutive visits to two recharging stations are forbidden.
Each customer i ∈ C has a demand q i . We assume that the depot and satellites are large enough to handle the demand of all customers. The satellites are intermediate stations that are used to transfer goods from the depot to the customers. Each route of the first echelon starts from and ends with the depot after visiting some satellites. The goods delivered by a truck cannot exceed the load capacity Q 1 . Each satellite can be visited by more than one truck. However, a satellite can be visited by a truck no more than once. Any route of the second echelon starts from a satellite, visits a subset of customers and recharging stations, and then returns to the same satellite. The total quantity of goods delivered by an EV cannot exceed the load capacity Q 2 and the remaining battery power of EVs cannot be negative. Each customer can only be visited once. Table 1 summarizes the notation used throughout the paper. Then, the 2E-EVRP problem can be modeled as the following mixed integer programming (MIP): The remaining battery power when EV k arrives and leaves vertex i from satellite s, respectively;

b ik
The amount of goods delivered to satellite i by truck k; x i jk Binary decision variable demonstrating whether truck k passes through arc (i, j); y i jsk Binary decision variable demonstrating whether EV k from satellite s passes through arc (i, j). Min The objective function (1) is to minimize the total cost including the fixed cost and the transport cost in both echelons. Constraints (2) guarantee that a truck visits a satellite no more than once. Constraints (3) and (8) show the flow balance of vehicles. Constraints (4) and (9) ensure the capacity of vehicles cannot be exceeded. Constraints (5) indicate the balance of goods receiving and delivery, which means that all goods received by a satellite from the depot must be delivered to customers. Constraint (6) prevents the situation that an EV returns to a different satellite from the starting satellite. Constraints (7) specify that each customer can only be visited once. Constraints (10) reveals the relationship of battery power when an EV transfers from vertex i to vertex j, which also eliminates the sub-tours in the route. Constraints (11) impose that the battery power is full when an EV departs from a satellite or a recharging station. Constraints (12) keep the battery power of an EV the same when it arrives at or leaves a customer.

Dantzig-Wolfe decomposition
We shall use the idea of Dantzig-Wolfe decomposition to propose an algorithm for the problem. The MIP model (1)-(16) will be decomposed into an integer master problem and a pricing subproblem by the Dantzig-Wolfe decomposition. Let R 1 and R 2 be the feasible routes sets of the first and second echelon, respectively. Let R s 2 be the set of feasible routes in the second echelon starting from and ending at the satellite s, where R 2 = s∈S R s 2 . Let λ r be an integer variable that indicates the times route r ∈ R 1 is implemented. For each route r ∈ R 1 ∪ R 2 , c r denotes its transport cost and a i r (i ∈ V \{0}) equals to 1 if vertex i is visited by route r and 0 otherwise. Let p s r be the quantity of the goods delivered to satellite s using route r ∈ R 1 . For a route r ∈ R s 2 , θ s r is a decision variable that is equal to 1 if route r is a component of the optimal distribution strategy and 0 otherwise. Table 2 summarizes the notation of the integer master problem used throughout the paper. The integer master problem (IMP) is shown as follows: The objection function (17) is the sum of the fixed cost and transport cost in both echelons. Constraints (18) ensure the flow balance at each satellite, which means the total amount of goods received by a satellite from the depot equals that shipped from the satellite to customers. Constraints (19) ensure that each customer is visited only once in the second echelon. Constraints (20) guarantee that the capacity of the trucks is not exceeded. Constraint (21) imposes the relationship between a s r and p s r that a route of the first echelon cannot deliver goods to a satellite if the route does not visit it. Constraints (22)-(24) restrict the domains of the associated variables.
We relax constraints (22) and (23) as λ r ≥ 0 and θ s r ≥ 0 to generate lower bounds for the model. The linear relaxation of the IMP (17)-(24) is called the master problem. The lower bound of the master problem can be strengthened by imposing the capacity constraint which imposes a minimum number of vehicles used at the first echelon to transport the goods from the depot to satellites.

Branch-and-price algorithm
For model (17)- (25), there are too many variables to solve by a standard MIP solver or the branch-and-bound algorithm.
Therefore, we propose a BP algorithm for it based on its structure. Noting that each route r ∈ R 1 ∪ R 2 corresponds to a variable (λ r or θ s r ) in the master problem, we keep the set of routes R 1 unchanged and replace the sets R s 2 (s ∈ S) by the initial subsets of routes R s We call the master problem with the variables (columns) λ r (r ∈ R 1 ) and θ s r (s ∈ S, r ∈ R s 2 ) the restricted master problem (RMP). Then, we adopt the column generation (CG) algorithm to solve the RMP. We branch the RMP into two subproblems based on a branching strategy if the solution of the RMP is fractional. Therefore, our algorithm comprises three main components: (i) an insertion heuristic (IH) to generate the initial subsets of routes R s 2 ⊂ R s 2 (s ∈ S), (ii) a CG algorithm that evaluates the lower bounds of IMP (17)-(25) by solving the RMP, and (iii) a branching strategy that forces the solution of the RMP to be integers.

The branch-and-price algorithm framework
In this subsection, we discuss the BP algorithm, which is called Algorithm 1. At first, the initial set of variables (columns) ζ is generated. Based on the initial set of columns ζ , we call CG to calculate the new lower bound. Then, we explore the branch nodes of search tree T and update the lower bound and upper bound. The optimal solution r * is obtained when the lower bound L B equals the upper bound U B.

Algorithm 1 The BP algorithm
Input: Data on the depot, satellites, customers, and vehicles, and the maximum running time T max ; Output: The optimal solution r * ; 1: Initialize the columns set ζ =∅ and the search tree T =∅; 2: Call IH to generate a feasible solution ψ of the second echelon and get the initial columns set ζ . Calculate U B 0 , U B ← U B 0 (see Sect. 4.2); 3: Call CG to solve P 0 , and get the solution S 0 and its cost C Choose a branch node P such that C(P) = min P ∈T {C(P )}. Obtain its solution S; if S 0 is an integer solution and C(P 0 ) < U B then 9: U B ← C(P 0 ) and prune the branch nodes P such that C(P) ≥ U B; 10: else if S 0 is a non-integer solution then 11: Branch P 0 into P 1 and P 2 , call CG to get the optimal solutions S 1 and S 2 and their costs C(P 1 ) and C(P 2 ); 12: T ← T ∪ P 1 , P 2 , L B ← min P ∈T

{C(P };
13: end if 14: end while 15: return r * Since the number of satellites is generally relatively small, we enumerate all the subsets of the satellites with the depot in the first echelon. Then, we find the shortest route of each subset that starts from and ends at the depot by solving the ESPPRC. These shortest routes constitute the set of the first echelon routes R 1 . In the second echelon, we apply an insertion heuristic (Algorithm 2) to generate a feasible solution ψ of the second echelon, where ψ = s∈S R s 2 . Then, we get the initial set of variables as ζ = {λ r |r ∈ R 1 } ∪ {θ s r |s ∈ S, r ∈ R s 2 }.
Once the feasible solution ψ of the second echelon is fixed, the 2E-EVRP can be transformed into solving the split delivery vehicle routing problem at the first echelon. We use the CPLEX to solve this problem and derive the initial upper bound U B 0 .
We denote T as the search tree to represent the original problem and its subproblems. The RMP with the initial set of columns is regarded as the original problem and is thought of as the root node of T . In each iteration, we select an unfathomed node from T which is denoted as P 0 . The optimal solution and its corresponding value of P 0 are denoted as S 0 and C(P 0 ), respectively. The upper bound in the current iteration is U B. If it is an integer solution, the U B will be updated and the unnecessary branch nodes will be pruned. If S 0 is not an integer solution and C(P 0 ) < U B, we branch P 0 into two subproblems P 1 and P 2 (child nodes) and add P 1 and P 2 to the search tree T . The optimal solutions of P 1 and P 2 can be obtained by invoking the CG algorithm. Since P 0 is fathomed, we remove P 0 from the search tree T . Then, the lower bound will be updated.

Insertion heuristic algorithm
In this subsection, we design an insertion heuristic (IH) algorithm to obtain a feasible solution ψ of the second echelon, which consists of the routes that satisfy the resource constraints without en-route recharging. The algorithm consists of two steps: (i) assign customers to their nearest satellites; (ii) route the assigned customers of each satellite. We summarize IH algorithm as Algorithm 2.

Algorithm 2 Insertion heuristic
Input: Distance d i j , cost c i j , and electricity consumption e i j for each arc (i, j) ∈ A 2 , demand q i for each customer i ∈ C, EVs' capacity Q 2 , and the maximum recharging time E; Output:A feasible solution ψ of the second echelon; 1: For each s ∈ S, initialize set s = ∅, D s = ∅; ****** assign customers to its nearest satellite ****** 2: for i ∈ C do 3: for s ∈ S do 5: Calculate the distance d is between customer i and satellite s; Insert customer u ∈ s with the minimum inrsertion cost into route P n s , such that the load capacity and battery capacity constraints are satisfied, and remove u from s ; 17: until no such customer in s exists 18: n ← n + 1; 19: end while 20: end for 21: return ψ = s∈S P n s

Column generation algorithm
In this subsection, we describe the column generation (CG) algorithm that generates lower bounds for the IMP (17)- (25). To state the algorithm, we first introduce the pricing subproblem, which aims at finding a column with the minimum reduced cost.
For a given satellite s ∈ S, let α s and π i (i ∈ C) be the dual variables associated with constraints (18) and (19), respectively. Suppose there is a route r ∈ R 2 starting from satellite s. Let y i j be a binary variable to indicate whether this route passes through arc (i, j). h i j means the remaining battery power when EV leaves vertex i. a i denotes whether customer i is visited by this route. Then, the pricing subproblem is shown as follows: The objective function (26) seeks for the column with minimum reduced cost. Constraints (27)-(28) imply the relationship between y i j and a i . Constraints (29)-(30) restrict the route starting from and ending at the same satellite. Constraints (31)- (35) are the same as constraints (8)- (12). Constraints (36)- (38) are the variable domains.
The pricing subproblem is subject to cost, capacity, and electricity constraints, which are called resource constraints. We can re-formulate the subproblem as an ESPPRC in G =(V , A 2 ) where V = S ∪ F ∪ C. Since the ESPPRC is strongly NP-hard, we adopt dynamic programming to solve it.
To compute the reduced costs of the feasible routes in G , we replace the arc cost c i j by a modified cost c i j as follows: The variables π i and q i are set to 0 if i ∈ S ∪ F. For any feasible route r , the reduced cost c r equals the sum of the modified costs of the arcs in the route.
The CG algorithm is an iterative algorithm that alternates between solving the RMP and the pricing subproblem. We denote the set of feasible routes as R 2 , which is generated by a labeling algorithm.

Algorithm 3 CG Algorithm
Input:The RMP of a branch node and its column set ζ ; Output:The optimal solution S * ; 1: do { 2: Solve the RMP to get a solution S, and dual variables α s and π i ; 3: Use the dual variables α s and π i to construct the pricing subproblem; 4: Call the labeling algorithm to generate a feasible route set R 2 (see Sect. 4.3.1); 5: Calculate the reduced cost c r for each r ∈ R 2 and c k = min r ∈R 2 {c r };

Labeling algorithm
In this subsection, we discuss the labeling algorithm that generates a set of feasible routes R 2 of the second echelon. Each feasible route in R 2 starts from a satellite and ends at the same satellite, and satisfies a set of resource constraints. To present conveniently, we call the starting satellite s and the ending satellite t. Therefore, each feasible route corresponds to an s − t path. In the labeling algorithm, each label of vertex i represents a partial path from a satellite s to vertex i. Labels extend forwardly using the resource extension functions (REFs) and those that violate the resource constraints will be discarded [32]. Furthermore, a dominance rule is introduced to eliminate partial paths from which the s − t path with the minimum reduced cost cannot be yielded by extending them.
Because of the symmetry of our model, we use a bidirectional search to generate the set of feasible routes of the second echelon. The bidirectional search consists of two mono-directional searches, namely the forward search from s to t and the backward search from t to s. To present clearly the bidirectional search, we introduce the mono-directional search first.

Mono-directional search
A mono-directional forward search extends the labels from s to t. In the search, a partial path r that starts from s to vertex i can be represented by a label where X i is the reduced cost along path r , Y i is the cumulative load transferred along path r , Z i is the cumulative recharging time required, since last recharge along path r , and H cust k i is the number of visiting times for customer k ∈ C along path r .
The components of the initial label at satellite s are set to 0. Suppose that (X i , Y i , Z i , (H cust k i ) k∈C , s) is a label of vertex i ∈ s ∪ C ∪ F and j ∈ C ∪ F is a successor of i, i.e., arc (i, j) ∈ A 2 . When the partial path r extends from i to j, a new label (X j , Y j , Z j , (H cust k j ) k∈C , s) of vertex j ∈ C ∪ F is generated according to the following REFs: To avoid enumerating all the feasible partial paths, each new generated label is examined by the dominance rule. Given that REFs (40)- (43) are nondecreasing, we define the dominance rule as follows: Definition 1 Let L n = (X n , Y n , Z n , (H cust k n ) k∈C , s), n ∈ {1, 2}, be two labels associated with different paths starting from the same satellite s and ending at the same vertex i. Label L 2 is dominated by L 1 if and at least one of the inequalities is strict.
Bidirectional search In the bidirectional search, labels are extended both forward from satellite s to its successors and backward from t to its predecessors. Labels, REFs, and the dominance rules of the backward search can be defined analogously as those of the forward search.
In the bidirectional search, we select a resource as the critical resource for which its value is nondecreasing as the path extends. We choose the capacity of the vehicles as the critical resource and set the midpoint as Q 2 /2. The restriction condition for the forward search is A forward label is no allowed to be extended if its component Y exceeds Q 2 /2. Simultaneously, backward labels are extended backward from t except Y > Q 2 /2. Let L f w i and L bw i be the sets of forward and backward labels associated with vertex i, respectively. The forward and backward labels should be merged to form an s − t path. Let be the labels of a forward path and a backward path starting from the same satellite s and ending at the same vertex i, respectively. An s − t path can be yielded by merging these two labels only if the labels meet the following conditions: otherwise (50) and the additional conditions, which correspond to two cases of the merged vertex.

Acceleration strategies
We adopt two strategies to reduce the computing time to solve the pricing subproblem. One strategy is the statespace relaxation, which is proposed by Christofides et al. [9]. This strategy is achieved by removing the label component (H cust k ) k∈C that is used to record the number of visiting times. The other strategy is the ng-route relaxation introduced by Baldacci et al. [3]. The state-space relaxation has been used to solve the TSP [4,46]. In their paper, the component (H cust k ) k∈C is replaced by the sum of all customers' visits δ = k∈C H cust k . Here, we remove the component (H cust k ) k∈C from the label instead of replacing it by δ to obtain the label with relaxation. Let L f w and L bw be the sets of forward and backward labels with relaxation, respectively. The labeling algorithm is performed to extend the labels of L f w and L bw in each round of the CG algorithm. If the labeling algorithm fails to generate columns with negative reduced costs, it is implemented with the set of labels L f w and L bw .
For each vertex i ∈ C ∪ F, let N i be its neighborhood that contains the |N i | − 1 closest vertices (customers and recharging stations) and i itself. Let P = {s, i 1 , i 2 , ..., i p } be a feasible path associated with the subset (P) ⊆ V (P), where (P) is A ng-route is a non-necessarily elementary path P = {s, i 1 , i 2 , ..., i p =i} starting from a satellite s, visiting some vertices, and ending at vertex i, such that i / ∈ (P ), where P = {s, i 1 , i 2 , ..., i p−1 }is a ng-path. The label L i p is allowed to extend to a vertex i p+1 , such that i p+1 / ∈ (P) and all of the label components at vertex i p+1 satisfy the conditions. When N i contains all of the vertices, such that N i ⊇ C, it is obvious that label L i p generates the elementary paths. If |N i | is small enough, it is difficult to produce a near-elementary route, which means that the subproblem can be solved easily. To trade off the computing time and solution quality, we set |N i | to be 10 in the numerical studies.

Branching strategy
To derive an integer solution, we resort to branching. The branching strategy consists of three types of branching decisions [47]. Let {λ r |r ∈ R 1 } and { θ s r |s ∈ S, r ∈ R s 2 } be the solutions of the latest RMP. The variables {λ r |r ∈ R 1 } are preferred to branch first. Suppose λ w is the most fractional variable, i.e., its decimal part is closest to 0.5. Two child nodes are created, namely one adds the constraint λ r ≥ λ w and the other adds the constraint λ r ≤ λ w .
If all of the variables {λ r |r ∈ R 1 } are integers, we evaluate the number of vehicles r ∈ R s 2 θ s r used by each satellite s ∈ S.
Assume that satellite k uses the most fractional number of vehicles. Similarly, two child nodes are produced, namely one adds the constraint The third branching decision is used if {λ r |r ∈ R 1 } and { r ∈ R s 2 θ s r |s ∈ S} are integers. Let t i jr be the number of times that a vehicle travels the arc (i, j) ∈ A 2 . The total number of times all vehicles passed through arc (i, j) is r ∈ R s 2 t i jr θ s r . Assuming that (u, v) ∈ A 2 is the most fractional arc, we create two child nodes, namely one adds the constraint r ∈R s 2 t uvr θ s r ≤ 0 and the other adds the constraint r ∈R s 2 t uvr θ s r ≥ 1.
The second and third branching decisions are imposed by adding inequalities to the RMP. Each inequality has a corresponding dual variable. For a given satellite s ∈ S, let I s B R be the set of the inequalities that have been added to the RMP at the current branch node and β s be the set of the corresponding dual variables of the second decision. Let I s B A be the set of inequalities and γ s be the set of dual variables of the third decision. Then, the subproblem is If I s B R or I s B A is empty, the corresponding part in (54) is 0.
To reduce the memory requirement, we adopt the depthfirst search strategy to explore the search tree.

Numerical studies
In this section, we present numerical studies to evaluate the effectiveness of the BP algorithm. We present the test instances used in the studies in Sect. 5.1 and the performance of the algorithm in Sect. 5.2. We present the sensitivity analysis of recharging station density, battery capacity, and fixed costs of EVs in Sect. 5.3.
We coded the algorithm in Python 3.6.3 and used CPLEX 12.8.0 to solve the RMP. We conducted the numerical studies on a computer with an Inter(R) Core(TM) i7-9700K CPU (3.6GHz) and 32GB of RAM memory.

The test instances
As the 2E-EVRP problem in our paper is different from those considered in existing studies, no benchmarks' instances are available to test the performance of our proposed algorithm. Therefore, we generate new instances for the numerical studies.
First, we generated a set of small-sized instances, labeled as Set 1, on the basis of instances in Breunig et al. [7]. We generated second and third sets of instances, labeled as Set 2 and Set 3, respectively, which are based on the instances proposed by Breunig et al. [7] and Desaulniers et al. [16], respectively. Set 1 has 12 instances, each of which has ten customers and no more than two recharging stations. Set 2 has 12 instances, each of which has 20 customers and two recharging stations. Set 3 also has 12 instances, each of which has 20 customers and two recharging stations. To apply the instances to our problem, we need to preprocess these instances. First of all, we randomly removed some customers if the number of customers is great than our instances. Afterward, for the instances from Breunig et al. [7], we removed the recharging function as well as the capacity limit of the satellites. For the instances from Desaulniers et al. [16], we randomly inserted the satellites into the graph. All the instances can be downloaded by accessing URL: https:// github.com/Novigram/Insatances-for-2E-EVRP.git.

Computational results
In this subsection, we report the results of applying our BP algorithm to solve the generated instances. First, we compare the performance of the BP algorithm with that of CPLEX (arc flow model) on Set 1 (the small-sized instances) and report the results in Table 3. We give the names of the instances in the first column, e.g., "E-c10-s3-r1-1" means that the first instance, which has ten customers, three satellites, and one recharging station. This naming rule also applies to Set 2 and Set 3. The column "Opt" reports the optimal solution of CPLEX. If CPLEX cannot find the optimal solution within 3,600 s, the best upper bound is considered as the best approximate solution. The columns "UB" and "LB" denote the upper bound and lower bound founded on terminating the algorithm. The column "T(s)" shows the computation time (seconds) of CPLEX and the BP algorithm, respectively. The column "Gap" reports the percentage gap between the upper bound and lower bound.
defined as 100(UB-LB)/UB. From Table 3, we see that the BP algorithm can solve all of the 12 instances and find the optimal solutions. In addition, the BP algorithm solves the instances (except E-c10-s3-r1-5) faster than CPLEX. On average, the BP algorithm takes 13.9727 s, whereas CPLEX takes 47.2367 s to solve an instance in Set 1. The BP algorithm is 70.42% faster than CPLEX. Table 3 shows the correctness of the model and the effectiveness of our algorithm.
Next, we evaluate the performance of the BP algorithm in solving some larger scale problems (Set 2 and Set 3). We set the maximum CPU running time as 100,000 s. We report the computational results in Table 4.
In Table 4, the upper bound in column "UB" marked with an asterisk means that the BP algorithm finds the optimal solution. The column "RLB" denotes the root lower bound, which is the optimal value of the RMP at the root node.
From Table 4, we see that the BP algorithm performs well in solving the 2E-EVRP in moderate computing time. Overall, the BP algorithm solves 15 out of the 24 instances within the maximum running time. The BP algorithm finds the optimal solutions for seven instances in Set 2. For the other instances, the average percentage gap between the upper bound and lower bound is small, which is approximately 1.1 %. For the instances in Set 3, the BP algorithm solves eight instances within the maximum running time. For the other instances, the percentage gaps between the upper and lower bounds are very small and acceptable.
We can also observe that some instances are very hard to solve (E-c20-s2-r2-1, E-c20-s2-r2-2, etc.), while other are very easy to solve (E-c20-s2-r2-5, E-c20-s2-r2-9, etc.), despite they have the same dimensions. The computing time is mainly affected by the number of branch nodes. Each branch process means that the CG has to be executed many times, which is very time-consuming. If an integer solution of a branch node is produced, this node is no longer branched. Therefore, the earlier the integer solution is generated, the fewer of branch nodes are produced. If the search tree has fewer branches, the running time will be smaller.
In Sets 5 and 6, we fixed the numbers of the depot, satellites, and recharging stations for each instance. However, we varied the number of recharging stations n r in the different groups of Set 4. For each instance of Set 4, we added the additional recharging stations in the graph when n r changes. Suppose the x-coordinates of all the vertices (depot, satellites, and recharging stations) are in the interval [x min , x max ] and the y-coordinates are in the interval [y min , y max ]. We randomly located the additional recharging stations in the rectangular area [x min , x max ]×[y min , y max ]. Figure 2 depicts the relationship between the number of recharging stations and the total cost. Figure 3 shows the changes in the total cost concerning the battery capacity. Figure 4 shows the impacts of the fixed cost of the EVs on the transport cost (transport dis-  Impact of fixed cost on the transport cost tance). The "Gap" in Figs. 2 and 3 means the percentage gaps between the total costs with and without the battery capacity restriction. Since the total cost consists of the transport and fixed costs, changing the fixed cost directly affects the total cost, which cannot show the impact on the transport cost. Therefore, the "Gap" in Fig. 4 is the percentage gap between the transport costs with and without the vehicle fixed cost. In Fig. 2, the points on each line n ∈ {1,2,3,4,5} represent the five instances with n r = {2,4,6,8,10}, which are derived from instance n. Overall, the percentage gap decreases with the recharging station density, but the decline is slower as the recharging station density becomes larger. In other words, the total cost will decrease with increasing recharging station density. However, when the recharging station density reaches a certain level, adding new recharging stations will result in a slower decline in the total cost.
In addition, there is no significant change in the percentage gaps on the line 4 and line 5, i.e., the recharging stations density has a slight impact on the total cost. This is because the optimal values of instance 4 and instance 5 with n r = 2 are very close to their corresponding optimal values without the battery capacity restriction. When new recharging stations are added, there is no better option that can significantly narrow the gap. Thus, the location of the recharging stations is very important. If the recharging stations are located at appropriate positions, the BP algorithm can generate promising solutions even if the recharging station density is very low.
In Fig. 3, the points on each line n ∈ {1,2,3,4,5} represent the five instances with different battery capacity, which are derived from instance n. When the battery capacity is small, a slight increase in the battery capacity can considerably reduce the transport cost. When the battery capacity is large, increasing the battery capacity has a small impact on the total cost. Moreover, when the recharging stations density is high, the effect of increasing the battery capacity to save cost is not obvious. In other words, when the recharging infrastructure is adequate, the battery capacity has a small impact on the transport cost.
In Fig. 4, the points on each line n ∈ {1,2,3,4,5} represent the five instances with different fixed costs, which are derived from instance n. Figure 4 shows the impact of the fixed costs of the EVs on the transport cost in the line 1 and line 3. When the fixed cost increases to a certain extent, the vehicle will change from the original distribution route to visiting more customers, so the number of vehicles used is reduced. However, the fixed cost has no impact on the transport cost in the other three lines. This is mainly caused by one or more of the following reasons: (i) The high loading rates of the EVs make it impossible for the EVs to access more customers.
(ii) The distribution of customers in each region is relatively centralized, so the distances between the different regions are long.
(iii) The battery capacity of the EVs is small and there is no available recharging station nearby, which makes the EVs unable to undertake more distribution tasks.

Conclusions
In this paper, motivated by the express and e-commerce companies' distribution practices, we study the 2E-EVRP, which is an extension of the 2E-VRP considering battery capacity and en-route recharging. We propose a BP algorithm to solve the problem, which uses enumeration to generate feasible routes for the first echelon and the labeling algorithm for the second echelon. Conducting computational experiments to evaluate the effectiveness of the BP algorithm, we show that the algorithm can find the optimal solutions efficiently. We demonstrate the effectiveness of the BP algorithm by comparing it with CPLEX for solving small-scale test instances. We further assess the performance of the BP algorithm in solving middle-scale test instances. We also explore the impacts of recharging station density, battery capacity, and fixed cost of the EVs on the optimal solution.
Based on our study, we suggest some topics worthy of future research. In practice, a host of other factors, e.g., time windows, the limited number of EVs, nonlinear recharging, and discharging characteristics, need to be considered when manager vehicle schedules. Extending our model and algorithm to consider other factors is a topic worth studying in the future. In this paper, we propose a BP algorithm to solve the problem. Other solution algorithms, be they exact or approximate, could be devised to solve the problem.
Author Contributions Zhiguo Wu designed the proposed algorithm and conducted the numerical experiments. Juliang Zhang adjusted the structure and checked the grammar of this paper. Zhiguo Wu is the major contributor in writing the manuscript. All authors read and approved the final manuscript. Availability of data and materials All data generated or analyzed during this study are included in this published article [and its supplementary information files].

Conflicts of interest
We declare that we have no conflict of interest.
Ethics approval Not applicable.

Consent to participate Not applicable.
Consent for publication Written informed consent for publication was obtained from all participants.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.