Electric Vehicle Routing Problem with Battery Swapping Considering Energy Consumption and Carbon Emissions

In this paper, we study an electric vehicle routing problem while considering the constraints on battery life and battery swapping stations. We first introduce a comprehensive model consisting of speed, load and distance to measure the energy consumption and carbon emissions of electric vehicles. Second, we propose a mixed integer programming model to minimize the total costs related to electric vehicle energy consumption and travel time. To solve this model efficiently, we develop an adaptive genetic algorithm based on hill climbing optimization and neighborhood search. The crossover and mutation probabilities are designed to adaptively adjust with the change of population fitness. The hill climbing search is used to enhance the local search ability of the algorithm. In order to satisfy the constraints of battery life and battery swapping stations, the neighborhood search strategy is applied to obtain the final optimal feasible solution. Finally, we conduct numerical experiments to test the performance of the algorithm. Computational results illustrate that a routing arrangement that accounts for power consumption and travel time can reduce carbon emissions and total logistics delivery costs. Moreover, we demonstrate the effect of adaptive crossover and mutation probabilities on the optimal solution.


Introduction
As cleaner alternatives to fossil fuel-based vehicles (FFVs), electric freight vehicles (EFVs) have been widely used in logistics and transportation to deal with environmental challenges. Compared with conventional FFVs, EFVs are cleaner and quieter, which means less air pollution and noise. According to the US Department of Energy (DOE), EFVs can reduce greenhouse gas emissions by 48% compared to the traditional FFVs running on gasoline and diesel [1]. In order to reduce greenhouse gas emissions and improve fuel efficiency, governments have taken corresponding measures to support EFVs. The U.S. National Petroleum Council has confirmed that electric vehicles are a very promising tool for reducing urban air pollution [2]. Europe is committed to reducing its greenhouse gas emissions by 40% until 2030 [3]. In the USA, the maximum tax credit for electric vehicles at the national level is $7500, and states can adopt further incentives to facilitate increased purchases [4]. The U.S. Environmental Protection Agency (EPA) together with the National Highway Traffic Safety Administration (NHTSA) has adopted measures to promote clean-energy vehicles [5]. In 2013, the cabinet of the State Council of China approved the Energy Saving and New Energy Vehicle Industry Development Plan, aiming to In Section 3, we describe the BVRP-EC problem, introduce a comprehensive measurement method of electric vehicle energy consumption and propose a mathematical model for BVRP-EC. In Section 4, the algorithm of AGE-HN is presented in detail to solve BVRP-EC. We use the computational experiments to test and analyze the algorithm's performance in Section 5. Finally, Section 6 concludes this work and provides several potential research avenues for the future.

Literature Review
There has been a large amount of research conducted on the electric vehicle issues related to this paper. We divided the relevant literature into the following three streams: energy consumption evaluation of electric vehicles, the electric vehicle routing problem and location optimization of battery swapping stations.

Energy Consumption Evaluation of Electric Vehicles
Due to the fact that the carbon emissions of electric vehicles are directly proportional to their energy consumption, the evaluation of energy consumption is of great importance to reduce the impact of carbon emissions.  introduce the Electric Vehicle Routing Problem (EVRP) and formulate a corresponding mathematical model [23]. They illustrate the importance of using the energy-minimizing objective rather than the classical distance-minimizing objective. In the context of an electric vehicle routing problem, Pelletier et al. (2019) propose a robust optimization framework to take these energy consumption uncertainties into account [24]. By comparing the nominal solution with the robust solution, they also prove the importance of considering the uncertainty of energy consumption. Zhang et al. (2015) study a vehicle routing problem with fuel consumption and carbon emissions and design a tabu search algorithm named RS-TS for solving this problem [25]. Goeke and Schneider (2015) present a battery energy consumption model based on mileage, road slope, vehicle speed and load capacity [26]. Based on the efficiency and load of an energy storage system (ESS), Hwang et al. (2020) define the application standard to maximize the fuel efficiency of a series hybrid electric bus (SHEB) [27]. Zeng et al. (2018) propose a control strategy based on an adaptive simplified equivalent consumption minimum strategy (ECMS) [28], and the simulation results show that the fuel consumption of this strategy is reduced by 16.43% in the test.
Although the above-mentioned papers analyze the effect of energy consumption in electric vehicle transport operations, they do not consider the case of battery swapping. In this paper, we will incorporate the costs of energy consumption into the electric vehicle routing problem under the battery swapping mode.

Electric Vehicle Routing Problem
Unlike a traditional vehicle routing problem involving FFVs, electric vehicle routing problems incorporate EFVs into the delivery routes. The research on the electric vehicle routing problem has been relatively limited until now. Under the condition of a sharing economy,  investigate the urban logistics distribution routing problem of electric vehicles while considering a carbon tax and time of use prices [29]. Moreover, they put forward an optimization model of the electric vehicle routing problem. Conrad and Figliozzi (2011) take the charging demand of electric vehicles into account [30]. They propose a routing optimization model for rechargeable electric vehicles with the goal of minimizing the vehicle use cost and vehicle distribution cost (including the vehicle travel distance cost, service time cost and charging cost). Gerhard et al. (2016) introduce electric vehicles of various models and formulate a vehicle routing model with time window constraints to minimize the number of vehicles and distribution distance [31]. Bahrami et al. (2020) study the plug-in hybrid vehicle routing problem, establish an energy management model with the goal of minimizing total energy consumption, and propose an accurate branch pricing method and a heuristic algorithm [32]. Alejandro et al. (2017) propose a piecewise linear method to fit the nonlinear process of battery charging [33]. Compared with the traditional charging function, the linearization method has a higher power control accuracy and reduces the calculation error for the charging time. Erdogan and Miller (2012) explore the electric vehicle routing problem and put forward two methods: an improved Clarke Wright saving algorithm and a density based clustering algorithm [34].
Most of the research reviewed above focuses on the electric vehicle routing problem in the context of the charging mode. However, they do not consider the effect of battery swapping and corresponding stations. In this paper, we consider the total costs related to driving time and energy consumption. Moreover, we consider the fact that even though a vehicle fleet is fixed, a battery station may be visited several times by that fleet's EFVs to accomplish a delivery task.

Location Optimization of Battery Swapping Stations
To ensure that electric vehicles never get stranded in their delivery routes, they must go to battery swapping stations to replace their batteries before they run out. Obviously, the location of the battery swapping stations has a direct effect on the routing problem. Jie et al. (2019) study a two-echelon capacitated electric vehicle routing problem with a swapping station [11]. They demonstrate the applicability of the proposed model by conducting extensive computational experiments. They show the efficiency of a column generation and an adaptive large neighborhood search (CG-ALNS). Wang et al. (2019) establish a new multi-period equilibrium planning model of electric vehicle charging stations [35]. They focus on the problem of mobile charging station (MCS) positioning and relocation and obtain an overall optimal solution. Yang and Sun (2015) first examine the location vehicle routing problem for an electric vehicle charging and swapping station [36]. Taking the battery swapping station and vehicle routing as the joint optimization objectives, an integer programming model for electric vehicles with capacity constraints is then established. Worley et al. (2012) take the capacity and mileage of vehicles as constraints and regard the location of charging stations and the vehicle routing problem of electric vehicles as joint optimization objectives [37]. They also establish an integer programming model. However, due to the complexity of the model, they provide no solution algorithm. Due to the limited travel range of e-buses, the deployment of a switching station network is very important. Liu et al. (2016) present a new method to determine the capacity and location of a battery swapping station (BSS), respectively [38]. By minimizing the net losses of a distribution system, they propose a model for optimal locations of BSSs. Moon et al. (2020) propose a mathematical programming model for the deployment of e-bus network switching stations [39]. Their purpose is to determine the location of battery swapping stations and the scheduling problem.
The extant studies exert their efforts on the strategic location optimization of battery swapping stations. However, as a type of public infrastructure, the locations of battery swapping stations are often fixed. As a result, the electric vehicle routing problem is less commonly studied with joint considerations of the fixed fleet, a limited driving range and a limited loading capacity.
This paper studies the BVRP-EC problem and builds a mixed integer programming model to minimize the costs related to power consumption and travel time. An adaptive genetic algorithm called AGE-HN is proposed based on hill climbing optimization and a neighborhood search. In this paper, we draw several conclusions as follows. First, the adaptive genetic algorithm is fast and effective in searching for satisfactory solutions. Second, a routing arrangement that takes power consumption and travel time into account can reduce the carbon emissions and total logistics distribution costs. Third, the adaptive crossover and mutation probabilities can avoid the algorithm of AGE-HN falling into a local optimum and improve the global search ability of the algorithm. Fourth, AGE-HN can converge towards the optimal solution faster than a traditional genetic algorithm.

Problem Description
For convenience, we first define some mathematical symbols: a complete directed graph G = (V, E) represents the freight distribution network for BVRP-EC, where V = C ∪ F is the set of nodes; C = {1, . . . , n} is the set of customer nodes; F = 0, n + 1, . . . , n + f is the set of battery swapping stations; Node 0 represents the depot-it is assumed that the depot also has a battery swapping station; f + 1 is the maximum number of available battery swapping stations; E = (i, j) i j ∈ V is the set of arcs; K = 1, 2, . . . , k represents a set of EFVs with the load capacity Q, and k is the maximum number of available EFVs; the distance of each arc (i, j) is set to be d ij ; v ij is the travel speed of EFVs in arc (i, j); p ij represents the energy consumption of EFVs in arc (i, j); l ij indicates the load carried by the electric vehicle in arc (i, j); the service time of EFVs at the battery swapping station f is s f , f ∈ F; q i represents the demand for goods at customer i, i ∈ C; T is the battery capacity of EFVs; y ijk represents the remaining electricity when the vehicle k completes the visit on arc(i, j). The decision variables are required in the following: Binary variable x ijk is 1 when the vehicle k visits arc (i, j), otherwise it is 0.
For the problem of BVRP-EC, there is a vehicle fleet of EFVs with a limited number of vehicles, load capacity and battery capacity. There are a set of customers with a deterministic level of demand, a depot and a set of battery swapping stations. Their locations are given and fixed. Each EFV starts from the depot to deliver goods for customers, and then returns to the depot after completing the task. Each customer is visited only once. The battery swapping station can replace depleted batteries with fully charged ones, which takes a much shorter time than the charging mode. We assume that the EFVs have full battery power after they reach the depot and the battery swapping stations. Unlike the traditional vehicle routing problem with FFVs, this assumption usually allows vehicles to visit the customer only once while BVRP-EC takes the battery swapping electric-vehicle as a carrier, which may need to be charged in the battery swapping station. In order to accomplish the delivery task, EFVs are allowed to visit a battery swapping station one or more times.
Meanwhile, the constraints of battery swapping stations and the battery driving range are considered. We define the usage cost of EFVs per hour as α, which is positively proportional to driving time, consisting of a driver's salary, the tolls for crossing the roads and bridges, etc. The battery swapping station collects the fee of the battery swapping service based on the power consumption of EFVs. In this context, we assume that the cost of power consumption is β per kilowatt hour. BVRP-EC will optimize the EFVs routing to minimize the total cost of vehicle usage and power consumption. In addition, as discussed above, the following conditions should be satisfied: A battery swapping station can be visited several times by EFVs; the load capacity of each vehicle must exceed the total customer demand in its route; each vehicle only has a one-way trip; and the driving range is limited because of battery capacity.

Calculation of Energy Consumption and Carbon Emissions of EFVs
There are different perspectives for vehicle energy consumption, such as Well to Wheel (WTW) and Tank to Wheel (TTW) [40]. The majority of the transportation fuel consumption and emissions use TTW mode, where the fuel loaded by the freight is converted into automobile power. In this paper, we take energy consumption and carbon emissions into account from the perspective of TTW. In logistics operations, the energy consumption of electric vehicles is affected by many factors, such as the speed, acceleration, weight, frontal area, distance and road conditions [41]. Among them, the dynamic change of vehicle speeds is affected by air resistance and road congestion, which then affects the consumption of electric power. The weight of electric vehicles includes empty weight and load, and the load decreases when the number of customers to be delivered is reduced. Following [15], the energy consumption of electric vehicles on arc (i, j) can be calculated as follows: where a ij = a + g sin θ ij + gC r cos θ ij is a constant related to the road section; parameter a represents vehicle acceleration; parameter g is the gravity constant; the road angle is θ ij and the rolling resistance coefficient is C r ; w is the weight of the empty electric vehicle; γ = 0.5C d Aρ is the coefficient related to EFVs; C d is the traction coefficient; A is the windward area of the vehicle; and ρ is the air density. The energy consumption calculation method used in this paper can comprehensively reflect the influence of the key factors, such as speed, load and driving distance, which makes the calculation of electric vehicle energy consumption more accurate and consistent with its practical applications. The carbon emissions of electric vehicles are related to the energy consumption on arc (i, j) and the energy emissions factor e of the electric vehicles. Thus, the measurement of carbon emissions is given as follows: where the electric energy emission factor e is the carbon dioxide emission per unit of electric energy. This factor is related to the type of electric vehicle and power generation and is a constant in a specific logistics distribution. Here, we set it to e = 0.69 kg/kwh [42,43]. In reality, due to economic and scientific uncertainties, there are a lot of uncertain factors in the calculation of electric vehicle carbon emissions, such as emission levels, prediction methods, related emission concentrations, temperature levels and so on. All these factors play important roles in the assessment of carbon emission costs [44]. Furthermore, there are also many uncertain elements that may indirectly affect vehicle carbon emissions by affecting future transportation activities, such as household motorization levels and structures, building costs and so on [45]. In order to facilitate calculations, we mainly consider some closely relevant factors such as vehicle speed, weight and travel distance, etc. In the future research, taking the uncertainty of carbon emissions into consideration will become one of our critical research directions.

Models
According to the calculation of EFVs energy consumption and carbon emissions, the mathematical model of BVRP-EC is described as follows: k∈K i∈V k∈K j∈V The objective function (3) minimizes the total costs consisting of the electric vehicle use cost and power consumption cost. Constraints (4) and (5) ensure that each customer node is visited only once.
Constraints (6) are the limited depot and the number of vehicles, which means that all EFVs start from the depot and return to it after completing the task. Constraints (7) are the flow conservation equations. Constraints (8) and (9) are the loading capacity of the electric vehicles. Constraints (10) mean that the electric vehicle is fully charged after charging in the battery swapping station. Constraints (11) state that the remaining power consumption after the electric vehicle passes through an arc (i, j) does not exceed the remaining power consumption on the node i minus the power consumption on this arc. Constraints (12) ensure that the electric vehicle has enough residual power to return to the depot. Constraints (13) reflect the binary decision variables.
It can be seen from the objective function of BE-LRP that there are two special cases in this problem. First, let α = 0, and this model will minimize the power consumption of electric vehicles. As explained earlier, since carbon emissions are directly proportional to the power consumption of electric vehicles, they are also equivalent to minimizing carbon emissions in transportation. Second, let β = 0, and the objective of the model is to minimize the travel time cost in transportation. Therefore, by considering the travel time cost and low-carbon cost of this model, it can reflect comprehensive evaluation of the economic and environmental performance in the distribution of electric vehicles.
The objective function of BVRP-EC considers the total costs of vehicle use and power consumption, while a traditional VRP usually takes the minimum total travel distance as the goal. Meanwhile, BVRP-EC also imposes constraints related to electric vehicle charging, power consumption and a fixed fleet. Moreover, it allows electric vehicles to visit a battery swapping station multiple times so as to satisfy delivery requirements. This also consists with the reality of logistics distribution.

Solution Algorithm for BVRP-EC
The classical VRP is a NP-hard problem [13] and the BVRP-EC is one of its extensions. In order to improve the solution efficiency, it is necessary to design a heuristic algorithm. Therefore, this paper proposes an adaptive genetic algorithm based on hill climbing optimization and neighborhood search (AGE-HN). An adaptive genetic algorithm is an intelligent algorithm whose parameters can be adjusted flexibly according to the search process [46]. Compared with the classic genetic algorithm, the crossover and mutation probabilities are adaptively adjusted with the population fitness. In order to strengthen the local search ability of the algorithm and ensure the algorithm converges to the optimal solution quickly, this paper also introduces the hill climbing search algorithm. The hill climbing algorithm uses the iteration process to improve the current solution until the solution reaches a "peak". In addition, in order to meet the constraints of battery range and battery swapping stations, a neighborhood search of the battery swapping stations is carried out for each route in the optimal solution generated by the adaptive genetic algorithm to further improve the optimal solution and finally obtain the optimal feasible solution. The basic flow chart of AGE-HN is shown in Figure 1.

Coding, Population Initialization and Fitness Function
Based on the characteristics of the BVRP-EC model, we use a natural number coding strategy. Firstly, we determine the individual chromosome length based on the number of depots, customer nodes and battery swapping stations. Then, individual coding element is added to the EFV routing. The natural number code is expressed as follows: The depot is 0; customer nodes are (1, 2, 3, . . . , n); and battery swapping stations are denoted by (n + 1, n + 2, n + 3, . . . , n + f ). We can decide the minimum number of electric vehicles for good delivery in terms of the total demand i∈C q i of all the customer nodes and the EFV's loading capacity. For example, we can assume that there are 10 customers, and at least 2 vehicles are needed, which is expressed by 2 chromosomes. The number of chromosomes is equal to the number of electric vehicles. The customer service node set of the first vehicle is (3,4,5,6). The position of genes (3,4,5,6) in the chromosome can be arranged randomly. Similarly, the service customer node set of the second vehicle is (7,8,9,1,2), and the position of genes (7,8,9,1,2) in the chromosome can be arranged randomly as well. As a result, a complete vehicle routing (0, 3, 4, 5, 6, 0, 7, 8, 9, 1, 2, 0) is created by these two chromosomes via the depot 0.
In the population, each individual corresponds to its own delivery route. The judgment of individual fitness mainly depends on two aspects: (a) Individual objective function value; (b) whether the individual meets the constraints. This paper adopts the natural number coding method, and the encoding process takes the constraints of the problem into account, such as the number of electric vehicles, vehicle load and driving range. In line with minimum value of objective function for BVRP-EC, the fitness function is represented by a maximum value, δ = 1/Z.

Selection Strategy
This paper uses a selection strategy combining the best individual and the Roulette method. First, the fitness of each individual is calculated. Second, individuals are ranked in descending order according to their fitness value. Individuals with the highest fitness value are retained to enter the subsequent genetic operation. Finally, the roulette method is used to determine whether the remaining individuals can enter the follow-up operation. The Roulette method calculates the ratio of individual fitness to total fitness in the population, δ i / i∈Ω δ i , where Ω is the set of all individuals in the population.
The larger the ratio is, the greater is the probability of individual selection.

Crossover Strategy
The crossover strategy adopted in this paper is to deal with chromosome genes properly and improve the solution space, so as to produce more excellent individuals. It is assumed that δ, δ min and δ max represent the average fitness, minimum fitness and maximum fitness of the population, respectively. We also set the parameters µ = (0.5, 1) and ν = (0, 0.5). δ/δ max is the ratio of the average fitness to the maximum fitness within the population, which reflects the fitness distribution within the population. The closer the value is to 0.5, the higher the population concentration becomes. δ min /δ max is defined as the ratio of minimum fitness to maximum fitness, reflecting the fitness concentration of the whole population. When it gets close to 0, it means that the algorithm may fall into the local optimal solution. On the other hand, when δ/δ max ≥ µ, δ min /δ max ≥ ν, this represents a group of individuals with centralized fitness, otherwise the individual fitness of the population is relatively scattered.
The crossover operation can affect the searching ability of the solution space through the operations of the chromosome gene, retaining excellent individuals and producing new ones. In order to avoid premature convergence, inbreeding and inferior individuals harming the population, the following adaptive crossover probability is given: where p c represents the crossover probability, U is the set of individual components of the parent, x 1 i and x 2 i are the ith component of the parent, x min and x max are the minimum and maximum components of an individual, respectively, while δ 1 and δ 2 represent the two fitness levels of the parent.
In the early stage of the algorithm, when the individuals in the population are relatively scattered, there is a large genetic difference among paternal individuals. The higher the fitness of the parent, the greater the probability of crossover operation. In the later stage of the algorithm, the individuals in the population tend to concentrate and there are two situations: (1) The genetic difference between the parents is small and the algorithm may converge to the optimal solution. It is necessary to improve the probability of the crossover operation so the average fitness can further approach the maximum value.
(2) Due to the large difference of gene performance among parents, it is easy to lead to a local optimal solution. This means a larger cross probability is needed.
In this paper, we decide whether to use the crossover strategy by judging the relationship between the crossover probability p c and the generated random real number in (0, 1). The crossover operation is performed only if p c is greater than this random value, as shown below. (1) Two chromosomes are randomly selected as the parent chromosomes of the crossover operation. (2) For two vehicle routes in the given parent chromosome, the crossover position is randomly selected and the sequential crossover operation is performed. (3) Two new chromosomes are generated by merging the daughter chromosomes of the two vehicle routes and removing the elements with continuous 0 values.

Mutation Strategy
Mutation operation imitates chromosome mutation in gene inheritance, thus increasing population and chromosome diversity in a solution space. In this paper, the adaptive mutation probability is used to enhance the diversity of individuals and make the algorithm jump out of the local optimal search process as soon as possible. The adaptive mutation probability is calculated as follows: In the early stage of the algorithm, when the individuals in the population are scattered, the population diversity in the solution space is high and a fixed mutation probability is used. However, in the later stage of the algorithm, the individuals in the population are concentrated and the population diversity in the solution space is low. In this case, the mutation probability is adaptively adjusted according to the fitness. In other words, if the fitness of the individual is greater than the average fitness, the probability of mutation operation decreases with the fitness value, otherwise the probability of mutation operation increases.
In this paper, we first decide whether a mutation operation is carried out by judging whether the adaptive crossover probability p m is greater than the real number randomly generated in (0, 1). If it is less than the real number, the mutation operation is canceled; otherwise, the following operations are conducted. (1) Following the aforementioned rules, the selected chromosome is regarded as the parent chromosome. (2) In the parent chromosome, the mutation positions of all vehicle routes are determined randomly, and the crossover mutation operation is carried out. (3) The successive 0 elements are deleted and the daughter chromosomes of vehicle routes are merged to generate a new chromosome.

Hill Climbing Optimization and Termination Criteria
The hill climbing optimization of the algorithm is mainly used in the following two stages. In the first stage, the initial population is optimized and the individual fitness of the population is improved. In the second stage, after the genetic process of the algorithm is executed, the optimal fitness individuals are optimized to realize the local optimization of the optimal solution of the genetic operation and the solution is further improved. The following steps are performed for the hill climbing optimization: Step 1. For the sub vehicle routes represented by chromosomes, the position of the gene is randomly selected and the crossover operation is performed.
Step 2. The fitness values of each individual before and after the crossover operation are compared. If it becomes larger, then the original chromosome is replaced; otherwise, the original chromosome remains unchanged.
Step 3. Determine whether the termination condition has been reached. If not, repeat steps 1 and 2 until the specified number of iterations is reached.
There are three common types of termination rules. First, when a given objective functions value is reached, the algorithm is terminated. This method usually needs to know the optimal objective value in advance. Second, if the current solution is not improved after several continuous iterations, the algorithm will terminate. Because the number of continuous iterations is difficult to estimate in advance, the algorithm may not be able to jump out of the local optimal solution. Third, based on the required calculation time and accuracy, the maximum number of iterations can be set to ensure that the algorithm converges to the optimal solution. Therefore, this paper selects the third termination rule.

Neighborhood Search for Battery Swapping Stations
According to the location of the customer nodes and the battery swapping stations, neighborhood search is used to ensure that each route in the optimal solution meets the constraints of battery driving range and battery swapping stations. As the location of the battery swapping stations is given, when the electric vehicle is driven to a customer node in the order of delivery, whether to carry out the battery swapping operation is decided according to whether the remaining power can reach the next customer node or not. If the remaining power is not enough to go to the next customer node, neighborhood search will find the nearest batter-swapping station based on the driving radius between the two customer nodes. The neighborhood search procedure for reaching a possible battery swapping station is shown in Figure 2. The detailed steps of the neighborhood search are as follows: Step 1. Calculate the remaining power once an electric vehicle arrives at a customer node.
Step 2. Judge whether the remaining power can reach the location of the next destination (customer node or depot). If yes, go to step 3; otherwise, go to step 4.
Step 3. Go to the next destination node and terminate the battery swapping search.
Step 4. Find the nearest battery swapping station locations and record the number of battery swapping stations that can be reached.
Step 5. Establish the nearest route scheme. Calculate the distance among the customer node, each battery swapping station and the next destination node.
Step 6. Compare the distance to each battery swapping station, and finally choose the shortest one.
Step 7. Stop the neighborhood search, go to Step 3.

AGE-HN Algorithm
According to the previous algorithm design and the key steps, the main steps of the AGE-HN algorithm are provided as follows: Step 1. Set the initial population size as Inn and the maximum number of algorithm iterations as Gnmax. Calculate the distance between each two vertices. Code the individual and generate the initial population which meets the constraints.
Step 2. Calculate the fitness of the initial population and record the individual with lowest fitness.
Step 3. Optimize the individual with the lowest fitness in the population using hill climbing, replace the original individual and set the initial iteration number of hill climbing, t = 0.
Step 4. Determine whether the maximum number of climbing optimization Tmax has been reached. If so, record the optimal solution of the current population; otherwise, go to Step 3, set t = t + 1.
Step 5. Rank the fitness level of the current population and calculate the selection probability of each individual. If an individual has the highest fitness in the current population, then the selection probability is set to be 1; for others, the selection probability is determined by the roulette rule. Set the initial iteration number of the algorithm as s = 0.
Step 6. Perform a crossover operation: According to the adaptive crossover probability, calculate the crossover probability of each individual. Randomly generate a real value between (0, 1). If the crossover probability is larger than the random value, perform the crossover operation; if the crossover probability is less than the random value, the crossover operation is not used.
Step 7. Perform the mutation operation: calculate the mutation probability of each individual in terms of the adaptive mutation probability. Generate a random real value between (0, 1). If the mutation probability is larger than the random value, perform the mutation operation; otherwise, the mutation operation is not used.
Step 8. Calculate the individual fitness of the new population and record the current individual with highest fitness.
Step 9. Optimize the individual with the highest fitness by using the hill climbing method. If the fitness value is improved, then replace the original individual. When the original value is updated, the iteration times of hill climbing will increase; afterwards, set the initial number, b = 0.
Step 10. Determine whether the maximum number of climbing optimization Bmax has been reached. If so, record the optimal solution of the current population; otherwise, go to step 9, set b = b + 1.
Step 11. Judge whether the maximum number of iterations Gnmax is reached. If so, record the current optimal individual and go to step 12; otherwise, go to step 5, s = s + 1.
Step 12. In the optimal solution, perform a neighborhood search for the battery swapping stations, and update the optimal feasible solution.
Step 13. Terminate the algorithm and output the optimal feasible solution.

Experimental Setting
There are five electric vehicles in a depot of a transportation company. Given the customers' location and demand, 29 orders in a certain period are randomly selected as customer data. There are five candidate battery swapping stations in total. We assume that the electric vehicle is fully charged when it departs from the depot and arrives at the battery swapping station. The battery swapping time is 0.1 h, and the average speed of the electric vehicle in the delivery is 60 km/h. The depot ( ), customer nodes ( ) and battery swapping stations ( ) is as shown in Figure 3. In the current electric logistics vehicle market, minivan and light trucks are the main vehicles, accounting for 31% and 19% of the total market, respectively, which are suitable for logistics transit transportation and terminal distribution. After considering the economic cost and practical application, this paper chooses a minivan electric vehicle as the delivery vehicle and uses the industrial electricity price of 0.8 CNY/kWh as the battery swapping cost. Combined with the actual situation and electric vehicle related parameters, the model parameters are set as follows [47]: The usage cost of electric vehicle is α = 120 CNY/h; the charging cost is β = 0.8 CNY/kWh; battery capacity is T = 35.7 kWh/L; the battery swapping service time is s f = 0.1 h; empty weight is w = 1325 kg; maximum vehicle load is Q = 595 kg; acceleration is a = 0 m/s 2 ; the gravity constant is g = 9.81 m/s 2 ; the angle of road section is θ ij = 0; the rolling resistance coefficient is C r = 0.01; the traction coefficient is C d = 0.7; the windward area of the vehicle is A = 0.378 m 2 ; the air density is ρ = 1.2041 kg/m 3 .

Parameter Tuning and Algorithm Analysis
The convergence of the algorithm is related to the parameter setting. This paper will consider some parameters such as the size of the initial population, adaptive crossover, mutation probability and hill climbing search times, and then optimize the setting of the algorithm, improve the convergence speed of the algorithm and find the optimal parameter setting range. A reasonable population size can ensure that the algorithm converges to the optimal solution faster. In order to find the optimal population size, the candidate population size is set as (50,70,90,110,130,150,170,190,210,230,250,270,290,310). For each population size, the algorithm runs 10 times, recording the average convergence generations and convergence probability. The convergence probability represents the ratio of times of convergence of the algorithm over 10 in the experiment, as shown in Figure 4. It can be seen from Figure 4 that the convergence probability is convex as the population size increases. If the initial population size is too large or too small, then the convergence probability will be reduced. At the same time, the average convergence generation is concave with the increase of the population size. When the population size is 210, the average convergence generation is the lowest and the convergence probability is high. Therefore, the initial population size is set as 210. Based on the initial population size mentioned above, we compare two cases: fixed crossover and mutation probability (case F) and adaptive crossover and mutation probability (case A). The ermination criterion is set to reach the maximum number of iterations of 200. The fixed crossover probability is set to 0.8 and the mutation probability is 0.06. We conduct 50 experiments and draw the convergence curve when these two cases obtain the same optimal solution in an experiment. Figures 5 and 6 show the convergence curves under cases F and A, respectively. Although both cases converge to the optimal solution, at the beginning of the algorithm, the optimal solution and average solution under case A are relatively poor. But after the 40th generation, case A increases the crossover probability of the algorithm, which makes the algorithm approach the optimal solution. In the later 170 generations of the algorithm, case A can jump out of the current optimal solution, improve the mutation probability of the algorithm and find the final optimal solution. Case F converges smoothly from the beginning to the end, and it will not dynamically optimize the problem according to the characteristics of individuals in the current population. Therefore, the adaptive crossover and mutation probabilities can avoid the algorithm falling into local optimum and make the algorithm more adaptive.
According to the determined initial population size and adaptive crossover and mutation probability strategies, we next explore the effect of hill climbing search times on the convergence speed and fitness value of the algorithm. Set the number of candidates for hill climbing search to (0, 10,15,20,25,35,40,50,210,100,210). For each instance, the algorithm runs 10 times and records the average convergence generations and average fitness value of the algorithm. The results are shown in Figure 7. Taking into account the curve characteristics of average fitness and average convergence generations, when the hill climbing search number is in the range of (25,40), the algorithm can converge faster, and the average fitness value is higher. In summary, in order to avoid the algorithm falling into the local optimal value and approach the optimal solution more quickly, we choose 35 as the hill climbing search number.

Analysis of Computational Results
All numerical instances are implemented in Matlab and solved by the proposed adaptive genetic algorithm (AGE-HN). Based on the aforementioned parameter tuning, the algorithm parameters are set as follows: The initial population size is 210, the maximum number of iterations is 350, and the number of hill climbing is 35. The algorithm is run 20 times. According to the measurement of carbon emissions and total cost, we list the best five solutions in Table 1. Table 1 shows that Solution 5 is the optimal one. Specifically, the optimal solution is represented as  Table 1, the optimal solution not only minimizes the total cost, but also reduces carbon emissions to the greatest extent. The optimal delivery route is shown in Figure 8.

Comparison of Algorithms
In order to illustrate the performance of the improved adaptive genetic algorithm (AGE-HN), we compare it with the basic genetic algorithm with fixed parameters. We take the fitness value of the current optimal solution as the termination criterion. The algorithm result is oriented to evaluate the convergence speed of the improved algorithm. If the fitness value of the algorithm reaches the optimal value, the algorithm converges and records the optimal solution. The initial population size is set as 200, the crossover probability of basic genetic algorithm is 0.8 and the mutation probability is 0.05. Figure 9 shows the results of the algorithm running 10 times. As can be seen from Figure 9, the average number of iterations of the improved adaptive genetic algorithm is between 50 and 100, which is significantly lower than that of the basic algorithm. Moreover, the probability of convergence to the optimal value of the improved algorithm is 12.1% higher than that of the basic algorithm. As a result, the proposed adaptive genetic algorithm can converge to the optimal solution faster and improve the efficiency of solving BVRP-EC.

Experimental Summary
Through the above computational experiments, we can draw some conclusions. First of all, the existing research [13,48,49] typically takes the minimum cost and/or distance as optimization goals from an economic perspective in logistics and supply chain operations. However, we explicitly consider both the economic and environmental performance, and seek to minimize the total cost as well as total carbon emissions. In contrast with some studies [26,31], we use a comprehensive calculation method of energy consumption and carbons emissions of EFVs with lots of factors, such as speed, weight, driving distance and road conditions, etc. Moreover, we find that the optimal electric vehicle routes with battery swapping can reduce both energy consumption and carbon emissions. Thus, firms can not only gain economic benefits, but also show their environmental responsibility by optimizing their electric vehicle routing problem.
Secondly, numerical experiments show that the constraint of the battery swapping station is one of the key conditions in the BVRP-EC model. If the government builds public battery swapping stations and encourages enterprises to use electric vehicles through policy design, it can save energy consumptions and reduce carbon emissions, and thereby achieve sustainable development. Therefore, the focus of electric vehicle routing problem with battery swapping is closely related to the theme of sustainable logistics and supply chain management. Finally, this experiment also has some limitations, which are worthy of further study in future research. For example, a more efficient algorithm can be designed to raise overall computation efficiency. We can also further explore the location problem of multiple depots and battery swapping stations in addition to the electric vehicle routing problem.

Conclusions
Nowadays, environmental protection has gained public attention all over the world. The development of electric vehicles represents the future of the automobile industry in many countries, such as Japan, the UK and China. Hence, it is of great practical significance to study and analyze logistics transportation based on electric vehicles. Optimizing electric vehicle routing is an important measure to realize energy savings and emissions reduction in logistics delivery. This paper extends the classical VRP to the battery swapping electric vehicle routing problem with energy consumption and carbon emissions (BVRP-EC). To this end, we introduce a comprehensive calculation method of energy consumption and carbon emission of electric vehicles with various factors such as speed, load and driving distance. Moreover, we establish a mixed integer programming model for BVRP-EC. This model considers several new constraints such as electric vehicle capacity, battery range and potentially multiple visits to battery swapping stations. Its objective is to minimize the total costs of the vehicle use and power consumption throughout the delivery. To solve BVRP-EC efficiently, we propose an adaptive genetic algorithm based on hill climbing and neighborhood search (AGE-HN). In this algorithm, the crossover and mutation probabilities are adaptively adjusted with the population fitness. The local search ability of the algorithm is enhanced by hill climbing optimization. In addition, a neighborhood search strategy is proposed to meet the battery range and battery swapping station constraints. Thus, the optimal solution of the model can be obtained. Finally, our numerical experiments show that the adaptive genetic algorithm can quickly and effectively obtain a satisfactory feasible solution to ensure the minimization of economic costs and carbon emissions. Furthermore, our algorithm comparison shows that the proposed adaptive genetic algorithm can converge to the optimal solution faster and improve the efficiency of solving BVRP-EC.
The research results in this paper can provide decision support tools for logistics enterprises using electric vehicles. They can also provide policy insights for the government on electric vehicle logistics. First, the promulgation of some low-carbon policies, such as a carbon tax and a carbon trading mechanism, will be conducive to the promotion of the use of electric vehicles. Second, based on the distribution of depots and customers, the construction of public battery swapping stations will help to improve the battery range of electric vehicles, thus reducing energy consumption and carbon emissions. Finally, due to the strong constraint of battery range, the government should encourage the research and development of advanced battery technologies, so as to get rid of an excessive dependence on charging stations.
On the other hand, our study can also include more practical issues in the future. For example, it is worthwhile to study the electric vehicle routing problem in the context of sharing economy, where social or private freight vehicles may be utilized simultaneously to fulfill delivery aims. In addition, more recently, en route recharging of electric vehicles has attracted attention, where the electric vehicle can be charged while driving [50]. Therefore, how to extend the driving range in the vehicle routing problem by combing battery swapping and en route recharging modes may be another promising research direction.