# Bi-Objective Optimization for Vehicle Routing Problems with a Mixed Fleet of Conventional and Electric Vehicles and Soft Time Windows

Authors: Peixin Zhao; Fanfan Liu; Yuanyuan Guo; Xiaoyang Duan; Yunshu Zhang
Journal: Journal of Advanced Transportation
 (2021)

---

## Abstract

With the growing interest in environmental protection and congestion, electric vehicles are increasingly becoming the important transportation means. However, electric vehicles currently face several adoption barriers including high purchasing price and limited travelling range, so the fleets where electric vehicles and conventional vehicles coexist are closer to the current fleet management status. Considering the impact of charging facilities and carbon emission, this paper proposes a vehicle routing problem with a mixed fleet of conventional and electric vehicles and soft time windows. A bi-objective programming model is established to minimize total operational cost and time penalty cost. Finally, the nondominated sorting genetic algorithm II (NSGA-II) is employed to deal with this problem. Furthermore, single-objective optimization is carried out for the two objectives, respectively, and the linear weighting method is also used to solve the problem. Through the contrast of these results and the NSGA-II results, the effectiveness of the algorithm in this paper is further verified. The results indicate that two objectives are contradictory to some extent and decision-makers need a trade-off between two objectives.


---

## Body

## 1. Introduction

Energy and environment are hot spots of global concern. Climate and environment problems caused by the greenhouse effect emerge endlessly, which greatly threaten the people’s health and daily life. Statistics from the International Energy Agency show that transportation is one of the main sources of carbon emissions, accounting for 24.64% of the total carbon emissions [1]. With the increasing social attention to environmental issues such as polluting emissions, noise, and congestion, more and more regulations are enacted to control carbon emissions from the transportation industry.Electric vehicles (EVs), which have only a small amount of greenhouse gas emissions and very little noise, are considered to be a green alternative to conventional vehicles (CVs) in many occasions. Many cities have started electric bus or taxi pilot projects, such as Tokyo, New York, Montreal, and Shenzhen. Shenzhen, in particular, has become the first fully electrified city for franchised bus vehicles. In airport ground support service, Civil Aviation Administration of China is gradually carrying out “oil-to-electricity” activities for airport ground special vehicles. In the logistics field, some enterprises have already begun to employ EVs in their last-mile distribution business, e.g., in the field of small parcel shipping [2] and fresh products distribution [3].Success in replacing existing CVs with EVs depends on several important factors. In addition to securing the sufficient financial investment, government support policies (e.g., regulations and subsidies), sufficient charging infrastructure, and advanced battery technology are all critical to improve market acceptance of EV fleets [4]. However, EVs have been facing formidable challenges for the limited power density of batteries, such as limited travelling range. In recent years, the wireless power transfer (WPT) is introduced as an advanced technology of addressing this problem with the hope of increasing the popularity of EVs [5]. Although the performance of this technology is remarkable in terms of flexibility, security, and esthetics, it has not been widely promoted and used in the field of EVs due to its high cost of infrastructure. Therefore, the attention of this paper is focused on wired charging EVs. Due to the high purchase cost and limited travelling range, it is still difficult for EVs to completely replace CVs [6]. Recognizing that there are obstacles to deploying large-size EV fleet, most enterprises start to gradually transition from just a few EVs. Hence, mixed fleets where EVs and CVs coexist are closer to the current fleet management status. The development of environmentally friendly and efficient transportation and distribution systems is an important challenge.The mixed fleet vehicle routing problem (MFVRP) is a special case of the heterogeneous vehicle routing problem (HVRP). However, the electricity constraints and charging decisions of EVs make the MFVRP considerably more complex. Ćirović et al. [7] investigated the green vehicle routing problem with a heterogeneous fleet composed of environmentally friendly and unfriendly vehicles, taking into account the fact that logistics operators have a limited number of environmentally friendly vehicles. However, when defining the route, friendly and unfriendly vehicles are considered separately, and the energy supplement of the vehicles is not considered. By using the bi-objective mixed integer programming models, Eskandarpour et al. [8] and Ren et al. [9], respectively, studied the MFVRP without considering the recharging of EVs, but Ren et al. [9] considered the time window of the customer. Using an improved multidirectional local search algorithm, the former aimed at minimizing distance cost and carbon emissions, while the latter minimized pollution emissions and delay time through an improved variable neighborhood search algorithm. Lu et al. [10] addressed the problem of scheduling a taxi fleet with CVs and EVs, designed an innovative multilayer taxi-flow time-space network, and proposed two simple but effective decomposition-based heuristics to efficiently solve the problem with practical sizes. In their paper, the recharging time of EVs was set to a fixed value.Goeke and Schneider [11], Masmoudi et al. [12], and Li et al. [13] considered the recharge plans of EVs in the same way when studying the MFVRP, with the constraints that EVs must be fully charged at the charging station. Among them, Li et al. [13] conducted sensitivity analysis on the fixed cost, variable cost, and travelling range of EVs and indicated that high fixed costs and low variable costs make EVs need long-distance distribution, which is restricted by the limited travelling range. In addition, electric vehicle routing problem with time window and mixed fleet based on the fully charged policy was examined by Li et al. [14]. An integer programming model is presented and the branch and price algorithm is developed to obtain the optimal solution. Sassi et al. [15, 16] had a more complicated consideration on the recharge plans of EVs; that is, the EVs can be partially recharged at the charging station. Hiermann et al. [17] also made the same assumption for charging and considered the time windows of customers. They rationally planned the routing of the mixed fleet including conventional vehicles, plug-in hybrid vehicles, and electric vehicles. Macrina et al. [18, 19], Mouhrim et al. [20], and Yu et al. [21] also supposed that EVs can be partially charged and further studied the mixed fleet vehicle routing problem with time windows. Paper [18] explicitly considered the restrictions on the pollution emissions of CVs and minimized the cost under the pollution emissions restrictions, and paper [19] innovatively added the acceleration and deceleration process of the EV to the energy consumption model to make further improvements. Mouhrim et al. [20] assumed that the energy consumption is related to the vehicle load and travel distance, and established a multiobjective programming model. Yu et al. [21] used the comprehensive modal emission model (CMEM), established a mixed integer programming model with the objective of minimizing the vehicles travelling cost and emissions, and solved this problem through an adaptive large neighborhood search heuristic. The relevant literature on the mixed fleet vehicle routing problem is shown in Table 1.Table 1
Summary of the literature on the mixed fleet vehicle routing problem.
ReferencesSingle-objectiveBi-objectiveHard time windowSoft time windowCarbon emissionEnergy consumptionChargingAlgorithm[8]√√——Enhanced MDLS[9]√√√——Improved VNS[10]√√DistanceFixed timeHeuristic[11]√√CMEMFullyALNS[12]√√CMEMFullyHALNS[13]√DistanceFullySS-MACO[14]√√DistanceFullyBranch and price[15]√DistancePartiallyHeuristic[16]√DistancePartiallyITS-LNS[17]√√DistancePartiallyMetaheuristic[18]√√√DistancePartiallyILS[19]√√Improved CMEMPartiallyHLNS[20]√√√Distance and loadPartially—[21]√√√CMEMPartiallyALNSOur work√√√Distance and loadFullyNSGA-IIThe literature on the VRP with a fleet composed of both EVs and CVs is limited (see Table1 for a summary). Most early research paid less attention to charging facilities, and then some researchers incorporated the charging decision of EVs into the consideration and made further improvements on the charging assumptions. In the majority of the aforementioned works, the authors mainly focus on the single-objective models and seldom consider the carbon emissions of EVs or soft time windows. Therefore, a Mixed Fleet Vehicle Routing Problem with Soft Time Windows (MFVRPSTW) considering carbon emissions of both EVs and CVs is proposed in this paper, and a bi-objective programming model is established to minimize the total operational cost and the time penalty cost. In contrast with previous studies, our work explicitly accounts for carbon emissions through taking it as a part of total cost. In order to solve the problem under investigation, we employ the nondominated sorting genetic algorithm II (NSGA-II), and numerical examples are given to verify the effectiveness of the algorithm.This paper is organized as follows. In Section2, we define the MFVRPSTW and its mathematical formulation. The NSGA-II is detailed in Section 3. The results of the computational experiments are provided in Section 4. Finally, Section 5 discusses the conclusions and directions for future research.


## 2. Mathematical Formulation for MFVRPSTW

### 2.1. Problem Statement and Assumptions

The enterprise dispatches vehicles from the distribution center to customers whose locations and time windows are known, and they return to the distribution center. The enterprise, which has the choice of EVs and CVs, requires each customer to be served and each vehicle to be dispatched only once. Due to the limited travelling range, the EVs can be recharged at any available charging stations. Referring to the EVs, CVs have a larger travelling range and fast energy replenishment speed and this process can be ignored. The relevant assumptions are as follows:(1)
The numbers of EVs and CVs are sufficient, and all vehicles have the same capacity.(2)
The time window of each customer is determined and known. The time window is soft; that is, the vehicles that arrive in advance can provide services but have to pay the waiting cost, while late arrivals also have to pay the delay cost.(3)
The number and locations of charging stations are given. For simplicity, the capacity of each station is assumed sufficient to simultaneously charge multiple EVs.(4)
The batteries of the EVs are fully charged before leaving the distribution center. Each EV can be recharged at any level of remaining battery power during its journey, but it must be fully charged before it can leave the station.(5)
Regardless of whether it is an EV or a CV, its energy consumption is directly proportional to the travel distance and load.


### 2.2. Energy Consumption for the Vehicles

This paper supposes that the energy consumption of vehicles is related to travel distance and load. Xiao et al. [22] have concluded that the energy consumption rate can be considered as a load dependent function through the statistical analysis of relevant data. Without losing any generalization, dividing the vehicle’s weight into two parts, the vehicle’s no-load weight Q0 and the carried load X, the energy consumption per unit distanceρ can be formulated as a linear function dependent on load X approximately as(1)ρX=aQ0+X+b.Defining the capacity of the vehicle asQ, the energy consumption per unit distance at full load as ρ∗, and the energy consumption per unit distance at no load as ρ0, it can be seen from formula (1) that ρ0=aQ0+b, ρ∗=aQ0+Q+b, and(2)a=ρ∗−ρ0Q.SoρX can be rewritten as(3)ρX=ρ0+ρ∗−ρ0QX.


### 2.3. Mathematical Model

In this problem, letC=1,2,…,n denote the set of customers, and V=VE∪VC be the set of vehicles, where VE is the set of EVs and VC is the set of CVs. F denotes the set of charging stations, and we also need copies of charging stations to account for multiple visits at the same station. Thus, let F′ be the set of all stations and their copies. For simplicity, let N=C∪F′. Vertices 0 and n+1 denote instances of the distribution center, and all routes start at 0 and end at n+1. Other sets are defined as follows: N0=N∪0, Nn+1=N∪n+1, N0,n+1=N∪0∪n+1, C0=C∪0, F0′=F′∪0. The parameters and variables used in the model are shown in Table 2.Table 2
Parameter and variable definitions of the model.
NotationDescriptionfkFixed cost of vehiclekckUnit energy cost of vehiclekεkCO2 emission coefficient of energy used by vehicle kceCharge cost per unit of electricitycoCarbon taxpeWaiting cost per unit timeplDelay cost per unit timeei,liTime window at vertexi, i∈CsiService time at vertexi,i∈CdiDemand of vertexi, i∈CdijTravel distance between verticesi and j,i,j∈NtijTravel time between verticesi and j,i,j∈NYBattery capacity of EVsgRecharging rateQVehicle capacityρkxEnergy consumption per unit distance when the load of vehiclek is xtiDecision variable specifying the time of arrival at vertexi,i∈NuikDecision variable specifying the remaining load of vehiclek when arriving at vertex i,i∈NyikDecision variable specifying the remaining battery capacity of EVk when arriving at vertex i,i∈NxijkBinary decision variable specifying whether arci,j is travelled by vehicle k,i,j∈NConsidering the charging constraints and soft time window constraints, the bi-objective mixed fleet vehicle routing optimization model is established as follows:(4)minf1=∑k∈V∑j∈Nfkx0jk+∑k∈V∑i∈N0∑j∈Nn+1ckρkujkdijxijk+ce∑k∈VE∑i∈F′∑j∈Nn+1Y−yikxijk+co∑k∈V∑i∈N0∑j∈Nn+1εkρkujkdijxijk,(5)minf2=pe∑i∈Cmaxei−ti,0+pl∑i∈Cmaxti−li,0,(6)s.t.∑k∈V∑j∈Nn+1,i≠jxijk=1,∀i∈C,(7)∑j∈Nx0jk=∑j∈Nxj,n+1k≤1,∀k∈V,(8)∑j∈Nn+1,i≠jxijk≤1,∀i∈F′,∀k∈VE,(9)∑i∈N0,i≠jxijk−∑i∈Nn+1,i≠jxjik=0,∀j∈N,∀k∈V,(10)ti+tij+sixijk−M1−xijk≤tj,∀k∈V,∀i∈C0,∀j∈Nn+1,i≠j,(11)ti+Y−yikg+tijxijk−M1−xijk≤tj,∀k∈VE,∀i∈F′,∀j∈Nn+1,i≠j,(12)ujk≤uik−dixijk+M1−xijk,∀k∈V,∀i∈N0,∀j∈Nn+1,i≠j,(13)0≤ujk≤Q,∀k∈V,∀j∈N0,n+1,(14)y0k=Y,∀k∈VE,(15)0≤yjk≤yik−ρkujkdijxijk+M1−xijk,∀k∈VE,∀i∈C,∀j∈Nn+1,i≠j,(16)0≤yjk≤Y−ρkujkdijxijk,∀k∈VE,∀i∈F0′,∀j∈Nn+1,i≠j,(17)xijk∈0,1,∀k∈V,∀i,j∈N0,n+1,i≠j,uik,yik,ti≥0,∀k∈V,∀i∈N0,n+1.The objective function equation (4) includes four terms: the fixed cost of vehicles, the energy consumption cost of vehicles, the charging cost of EVs, and the cost of carbon emissions. Equation (5) represents the time penalty cost due to the waiting or delaying service of vehicles, which is used to measure the quality of service. Constraint (6) stipulates that each customer is visited exactly once, and constraint (7) guarantees that each vehicle is dispatched once at most and must return to the distribution center after completing its tasks. Constraint (8) ensures that each charging station can be visited multiple times by any EV, and constraint (9) is the flow conservation constraint. Constraints (10) and (11) restrict arrival time at successor vertices j if the arc from i to j is traveled, and constraint (11) covers the case that the EVs visit the charging stations. Constraint (12) enforces the fulfillment of demand at customer vertices. Constraint (13) means that the remaining load of each vehicle does not exceed the capacity at any time. Constraints (14)–(16) are the energy constraints of EVs. Constraint (14) enforces that the EVs are fully charged when they depart from the distribution center. Constraints (15) and (16), respectively, indicate the remaining battery power constraints of successor vertices when the predecessor vertices are the customers and the charging stations. Finally, decision variables are defined in constraint (17).


## 2.1. Problem Statement and Assumptions

The enterprise dispatches vehicles from the distribution center to customers whose locations and time windows are known, and they return to the distribution center. The enterprise, which has the choice of EVs and CVs, requires each customer to be served and each vehicle to be dispatched only once. Due to the limited travelling range, the EVs can be recharged at any available charging stations. Referring to the EVs, CVs have a larger travelling range and fast energy replenishment speed and this process can be ignored. The relevant assumptions are as follows:(1)
The numbers of EVs and CVs are sufficient, and all vehicles have the same capacity.(2)
The time window of each customer is determined and known. The time window is soft; that is, the vehicles that arrive in advance can provide services but have to pay the waiting cost, while late arrivals also have to pay the delay cost.(3)
The number and locations of charging stations are given. For simplicity, the capacity of each station is assumed sufficient to simultaneously charge multiple EVs.(4)
The batteries of the EVs are fully charged before leaving the distribution center. Each EV can be recharged at any level of remaining battery power during its journey, but it must be fully charged before it can leave the station.(5)
Regardless of whether it is an EV or a CV, its energy consumption is directly proportional to the travel distance and load.


## 2.2. Energy Consumption for the Vehicles

This paper supposes that the energy consumption of vehicles is related to travel distance and load. Xiao et al. [22] have concluded that the energy consumption rate can be considered as a load dependent function through the statistical analysis of relevant data. Without losing any generalization, dividing the vehicle’s weight into two parts, the vehicle’s no-load weight Q0 and the carried load X, the energy consumption per unit distanceρ can be formulated as a linear function dependent on load X approximately as(1)ρX=aQ0+X+b.Defining the capacity of the vehicle asQ, the energy consumption per unit distance at full load as ρ∗, and the energy consumption per unit distance at no load as ρ0, it can be seen from formula (1) that ρ0=aQ0+b, ρ∗=aQ0+Q+b, and(2)a=ρ∗−ρ0Q.SoρX can be rewritten as(3)ρX=ρ0+ρ∗−ρ0QX.


## 2.3. Mathematical Model

In this problem, letC=1,2,…,n denote the set of customers, and V=VE∪VC be the set of vehicles, where VE is the set of EVs and VC is the set of CVs. F denotes the set of charging stations, and we also need copies of charging stations to account for multiple visits at the same station. Thus, let F′ be the set of all stations and their copies. For simplicity, let N=C∪F′. Vertices 0 and n+1 denote instances of the distribution center, and all routes start at 0 and end at n+1. Other sets are defined as follows: N0=N∪0, Nn+1=N∪n+1, N0,n+1=N∪0∪n+1, C0=C∪0, F0′=F′∪0. The parameters and variables used in the model are shown in Table 2.Table 2
Parameter and variable definitions of the model.
NotationDescriptionfkFixed cost of vehiclekckUnit energy cost of vehiclekεkCO2 emission coefficient of energy used by vehicle kceCharge cost per unit of electricitycoCarbon taxpeWaiting cost per unit timeplDelay cost per unit timeei,liTime window at vertexi, i∈CsiService time at vertexi,i∈CdiDemand of vertexi, i∈CdijTravel distance between verticesi and j,i,j∈NtijTravel time between verticesi and j,i,j∈NYBattery capacity of EVsgRecharging rateQVehicle capacityρkxEnergy consumption per unit distance when the load of vehiclek is xtiDecision variable specifying the time of arrival at vertexi,i∈NuikDecision variable specifying the remaining load of vehiclek when arriving at vertex i,i∈NyikDecision variable specifying the remaining battery capacity of EVk when arriving at vertex i,i∈NxijkBinary decision variable specifying whether arci,j is travelled by vehicle k,i,j∈NConsidering the charging constraints and soft time window constraints, the bi-objective mixed fleet vehicle routing optimization model is established as follows:(4)minf1=∑k∈V∑j∈Nfkx0jk+∑k∈V∑i∈N0∑j∈Nn+1ckρkujkdijxijk+ce∑k∈VE∑i∈F′∑j∈Nn+1Y−yikxijk+co∑k∈V∑i∈N0∑j∈Nn+1εkρkujkdijxijk,(5)minf2=pe∑i∈Cmaxei−ti,0+pl∑i∈Cmaxti−li,0,(6)s.t.∑k∈V∑j∈Nn+1,i≠jxijk=1,∀i∈C,(7)∑j∈Nx0jk=∑j∈Nxj,n+1k≤1,∀k∈V,(8)∑j∈Nn+1,i≠jxijk≤1,∀i∈F′,∀k∈VE,(9)∑i∈N0,i≠jxijk−∑i∈Nn+1,i≠jxjik=0,∀j∈N,∀k∈V,(10)ti+tij+sixijk−M1−xijk≤tj,∀k∈V,∀i∈C0,∀j∈Nn+1,i≠j,(11)ti+Y−yikg+tijxijk−M1−xijk≤tj,∀k∈VE,∀i∈F′,∀j∈Nn+1,i≠j,(12)ujk≤uik−dixijk+M1−xijk,∀k∈V,∀i∈N0,∀j∈Nn+1,i≠j,(13)0≤ujk≤Q,∀k∈V,∀j∈N0,n+1,(14)y0k=Y,∀k∈VE,(15)0≤yjk≤yik−ρkujkdijxijk+M1−xijk,∀k∈VE,∀i∈C,∀j∈Nn+1,i≠j,(16)0≤yjk≤Y−ρkujkdijxijk,∀k∈VE,∀i∈F0′,∀j∈Nn+1,i≠j,(17)xijk∈0,1,∀k∈V,∀i,j∈N0,n+1,i≠j,uik,yik,ti≥0,∀k∈V,∀i∈N0,n+1.The objective function equation (4) includes four terms: the fixed cost of vehicles, the energy consumption cost of vehicles, the charging cost of EVs, and the cost of carbon emissions. Equation (5) represents the time penalty cost due to the waiting or delaying service of vehicles, which is used to measure the quality of service. Constraint (6) stipulates that each customer is visited exactly once, and constraint (7) guarantees that each vehicle is dispatched once at most and must return to the distribution center after completing its tasks. Constraint (8) ensures that each charging station can be visited multiple times by any EV, and constraint (9) is the flow conservation constraint. Constraints (10) and (11) restrict arrival time at successor vertices j if the arc from i to j is traveled, and constraint (11) covers the case that the EVs visit the charging stations. Constraint (12) enforces the fulfillment of demand at customer vertices. Constraint (13) means that the remaining load of each vehicle does not exceed the capacity at any time. Constraints (14)–(16) are the energy constraints of EVs. Constraint (14) enforces that the EVs are fully charged when they depart from the distribution center. Constraints (15) and (16), respectively, indicate the remaining battery power constraints of successor vertices when the predecessor vertices are the customers and the charging stations. Finally, decision variables are defined in constraint (17).


## 3. Algorithm for MFVRPSTW

The traditional method of solving multiobjective model is to convert multiple objectives into a single objective. Like the linear weighting method, usually the decision-maker determines the weight coefficients of objectives in advance and integrates multiple optimization objectives into one objective. So the results obtained are subjective and cannot guarantee that multiple objectives are optimized at the same time. However, the method based on the Pareto optimal does not need to transform the objectives and can directly find the nondominated solution set, which is more widely used.The nondominated sorting genetic algorithm (NSGA-II) proposed by Deb et al. [23] improved the selection process of individuals by employing the density value estimation strategy, fast nondominated sorting strategy, and elite strategy, which has been widely applied in various fields. This algorithm is one of the most popular multiobjective algorithms, which not only ensures the value of the optimal solution and obtains the Pareto front, but also reduces the complexity of the algorithm. Therefore, the performance of this algorithm is remarkable in terms of running speed and convergence, and we employ it to solve the problem under investigation.


### 3.1. Algorithm Design

#### 3.1.1. Encoding

Each chromosome should carry information about the configuration of the fleet and the priority of the customer, so the chromosome is divided into two parts. The first part indicates the decision of fleet deployment, which is composed ofm genes (m is the number of vehicles). The value of the gene represents the type of vehicle, where 0 represents the EV and 1 represents the CV. The second part is composed of n genes (n is the number of customers), and the value of the gene indicates the customer number. For example, there are 3 vehicles in the distribution center to provide service for 5 customers, and the chromosome is shown in Figure 1.Figure 1
Representation of the chromosome.


#### 3.1.2. Decoding

First assign the customer to the corresponding vehicle. Starting from the first gene of the second part of the chromosome, if the integera satisfies ∑i=1aqi≤Q and ∑i=1a+1qi>Q, 0 is inserted after the a-th gene, and then this process is repeated from the inserted 0 until the end of the chromosome. If the number of 0 inserted is less than m−1, we increase the number of 0 after the last 0until its number is m−1, and finally, the head and the end of the second part will be inserted with a 0, respectively. The genes between two 0 indicate the route sequence of a vehicle. Through the above process, an initial distribution plan is generated.Secondly distinguish the type of vehicle and determine whether the EV needs to be charged. As shown in Figure2, R denotes the travelling range of the EV, and S1, S2, and S3 are the distance between the customers i and j, the distance between the customer j and its nearest charging station, and the distance between the customer i and its nearest charging station, respectively. If the remaining battery power cannot allow the EV to reach the nearest charging station after serving the next customer, the EV must be charged first to avoid stopping on the way when R<S1+S2. To ensure that the EV can reach the charging station, R must be equal to or greater than the distance from the current customer to the nearest charging station. Therefore, when S3≤R<S1+S2, the EV needs to recharge at the nearest charging station, and the charging station is inserted into its routing sequence.Figure 2
Schematic diagram of decoding.Take the decoding result in Figure3 as an example. The first vehicle is a CV with route 0-3-2-0, which means that the CV departs from the distribution center, serves customers 3 and 2, and finally returns to the distribution center. The second vehicle is an EV with route 0-1-4-6-5-0, which means after serving customers 1 and 4 the EV enters the charging station 6 for charging and then continues to serve customer 5. The third is an EV, but it is not used.Figure 3
Chromosome decoding result.


#### 3.1.3. Population Initialization

The first part of the chromosome represents the type of vehicle, and0 and 1 are randomly generated as the genes; the second part is realized by randomly arranging all customer numbers. Therefore, these two parts form an initial chromosome. A total of 2N (N is the population size) individuals are generated by repeating the above process, and then the best N individuals are selected as the initial population through nondominated sorting and crowding distance operators.


#### 3.1.4. Selection, Crossover, and Mutation Operators

This paper employs the tournament selection method; that is, each time a certain number of individuals are drawn from the population, the best one is selected, and the operation is repeated until the given population size is reached.The chromosome is designed to be divided into two parts in this paper, so crossover and mutation operations are performed on the two parts, respectively, and the crossover and mutation operators should be selected according to the different encoding methods. For crossover, the first part of the chromosome uses the single-point crossover, and the second part uses the partial-mapped crossover. As for mutation, we employ the two-point interchange mutation operator for both parts.


### 3.2. Steps of the Algorithm

The calculation procedure of this bi-objective optimization model is obtained based on the basic flow of NSGA-II in Figure4:(i)
Step 1. Initialize the evolutional generation n=0, and generate an initial population Pn with the number of individuals N.(ii)
Step 2. Generate the offspring population Qn.(iii)
Step 3. Merge the parent and offspring populations; that is, Rn=Pn∪Qn.(iv)
Step 4. Calculate the nondominance rank and crowding distance of each individual in Rn according to Steps 2.1 and 2.2, and select the best N individuals to obtain the initial population of the next iteration Pn+1.(v)
Step 5. Determine whether the termination conditions for iteration are met. If so, stratify the population and output the individuals of the first layer; otherwise, go to Step 2.Figure 4
Basic flow of the NSGA-II algorithm.


## 3.1. Algorithm Design

### 3.1.1. Encoding

Each chromosome should carry information about the configuration of the fleet and the priority of the customer, so the chromosome is divided into two parts. The first part indicates the decision of fleet deployment, which is composed ofm genes (m is the number of vehicles). The value of the gene represents the type of vehicle, where 0 represents the EV and 1 represents the CV. The second part is composed of n genes (n is the number of customers), and the value of the gene indicates the customer number. For example, there are 3 vehicles in the distribution center to provide service for 5 customers, and the chromosome is shown in Figure 1.Figure 1
Representation of the chromosome.


### 3.1.2. Decoding

First assign the customer to the corresponding vehicle. Starting from the first gene of the second part of the chromosome, if the integera satisfies ∑i=1aqi≤Q and ∑i=1a+1qi>Q, 0 is inserted after the a-th gene, and then this process is repeated from the inserted 0 until the end of the chromosome. If the number of 0 inserted is less than m−1, we increase the number of 0 after the last 0until its number is m−1, and finally, the head and the end of the second part will be inserted with a 0, respectively. The genes between two 0 indicate the route sequence of a vehicle. Through the above process, an initial distribution plan is generated.Secondly distinguish the type of vehicle and determine whether the EV needs to be charged. As shown in Figure2, R denotes the travelling range of the EV, and S1, S2, and S3 are the distance between the customers i and j, the distance between the customer j and its nearest charging station, and the distance between the customer i and its nearest charging station, respectively. If the remaining battery power cannot allow the EV to reach the nearest charging station after serving the next customer, the EV must be charged first to avoid stopping on the way when R<S1+S2. To ensure that the EV can reach the charging station, R must be equal to or greater than the distance from the current customer to the nearest charging station. Therefore, when S3≤R<S1+S2, the EV needs to recharge at the nearest charging station, and the charging station is inserted into its routing sequence.Figure 2
Schematic diagram of decoding.Take the decoding result in Figure3 as an example. The first vehicle is a CV with route 0-3-2-0, which means that the CV departs from the distribution center, serves customers 3 and 2, and finally returns to the distribution center. The second vehicle is an EV with route 0-1-4-6-5-0, which means after serving customers 1 and 4 the EV enters the charging station 6 for charging and then continues to serve customer 5. The third is an EV, but it is not used.Figure 3
Chromosome decoding result.


### 3.1.3. Population Initialization

The first part of the chromosome represents the type of vehicle, and0 and 1 are randomly generated as the genes; the second part is realized by randomly arranging all customer numbers. Therefore, these two parts form an initial chromosome. A total of 2N (N is the population size) individuals are generated by repeating the above process, and then the best N individuals are selected as the initial population through nondominated sorting and crowding distance operators.


### 3.1.4. Selection, Crossover, and Mutation Operators

This paper employs the tournament selection method; that is, each time a certain number of individuals are drawn from the population, the best one is selected, and the operation is repeated until the given population size is reached.The chromosome is designed to be divided into two parts in this paper, so crossover and mutation operations are performed on the two parts, respectively, and the crossover and mutation operators should be selected according to the different encoding methods. For crossover, the first part of the chromosome uses the single-point crossover, and the second part uses the partial-mapped crossover. As for mutation, we employ the two-point interchange mutation operator for both parts.


## 3.1.1. Encoding

Each chromosome should carry information about the configuration of the fleet and the priority of the customer, so the chromosome is divided into two parts. The first part indicates the decision of fleet deployment, which is composed ofm genes (m is the number of vehicles). The value of the gene represents the type of vehicle, where 0 represents the EV and 1 represents the CV. The second part is composed of n genes (n is the number of customers), and the value of the gene indicates the customer number. For example, there are 3 vehicles in the distribution center to provide service for 5 customers, and the chromosome is shown in Figure 1.Figure 1
Representation of the chromosome.


## 3.1.2. Decoding

First assign the customer to the corresponding vehicle. Starting from the first gene of the second part of the chromosome, if the integera satisfies ∑i=1aqi≤Q and ∑i=1a+1qi>Q, 0 is inserted after the a-th gene, and then this process is repeated from the inserted 0 until the end of the chromosome. If the number of 0 inserted is less than m−1, we increase the number of 0 after the last 0until its number is m−1, and finally, the head and the end of the second part will be inserted with a 0, respectively. The genes between two 0 indicate the route sequence of a vehicle. Through the above process, an initial distribution plan is generated.Secondly distinguish the type of vehicle and determine whether the EV needs to be charged. As shown in Figure2, R denotes the travelling range of the EV, and S1, S2, and S3 are the distance between the customers i and j, the distance between the customer j and its nearest charging station, and the distance between the customer i and its nearest charging station, respectively. If the remaining battery power cannot allow the EV to reach the nearest charging station after serving the next customer, the EV must be charged first to avoid stopping on the way when R<S1+S2. To ensure that the EV can reach the charging station, R must be equal to or greater than the distance from the current customer to the nearest charging station. Therefore, when S3≤R<S1+S2, the EV needs to recharge at the nearest charging station, and the charging station is inserted into its routing sequence.Figure 2
Schematic diagram of decoding.Take the decoding result in Figure3 as an example. The first vehicle is a CV with route 0-3-2-0, which means that the CV departs from the distribution center, serves customers 3 and 2, and finally returns to the distribution center. The second vehicle is an EV with route 0-1-4-6-5-0, which means after serving customers 1 and 4 the EV enters the charging station 6 for charging and then continues to serve customer 5. The third is an EV, but it is not used.Figure 3
Chromosome decoding result.


## 3.1.3. Population Initialization

The first part of the chromosome represents the type of vehicle, and0 and 1 are randomly generated as the genes; the second part is realized by randomly arranging all customer numbers. Therefore, these two parts form an initial chromosome. A total of 2N (N is the population size) individuals are generated by repeating the above process, and then the best N individuals are selected as the initial population through nondominated sorting and crowding distance operators.


## 3.1.4. Selection, Crossover, and Mutation Operators

This paper employs the tournament selection method; that is, each time a certain number of individuals are drawn from the population, the best one is selected, and the operation is repeated until the given population size is reached.The chromosome is designed to be divided into two parts in this paper, so crossover and mutation operations are performed on the two parts, respectively, and the crossover and mutation operators should be selected according to the different encoding methods. For crossover, the first part of the chromosome uses the single-point crossover, and the second part uses the partial-mapped crossover. As for mutation, we employ the two-point interchange mutation operator for both parts.


## 3.2. Steps of the Algorithm

The calculation procedure of this bi-objective optimization model is obtained based on the basic flow of NSGA-II in Figure4:(i)
Step 1. Initialize the evolutional generation n=0, and generate an initial population Pn with the number of individuals N.(ii)
Step 2. Generate the offspring population Qn.(iii)
Step 3. Merge the parent and offspring populations; that is, Rn=Pn∪Qn.(iv)
Step 4. Calculate the nondominance rank and crowding distance of each individual in Rn according to Steps 2.1 and 2.2, and select the best N individuals to obtain the initial population of the next iteration Pn+1.(v)
Step 5. Determine whether the termination conditions for iteration are met. If so, stratify the population and output the individuals of the first layer; otherwise, go to Step 2.Figure 4
Basic flow of the NSGA-II algorithm.


## 4. Numerical Examples

### 4.1. Computational Results

In this section, the effectiveness of the algorithm is illustrated by three numerical instances. Referring to [11], this paper sets the number of charging stations as ⌊0.1⋅N⌋ (N is the number of customers). The instances used are generated based on the Solomon benchmark set, as shown in Table 3. Meanwhile, we assume that the number of vehicles is sufficient, so set the number of vehicles equal to the number of customers.Table 3
Number of customers and charging stations.
Instance IInstance IIInstance IIINumber of customers255091Number of charging stations259Since the final result will be affected by the location of the charging station, it is necessary to select the reasonable charging station locations by a location model.(18)min∑i∈I∑j∈Jdijxij,s.t.∑j∈Jxij=1,∀i∈I,xij≤yj,∀i∈I,j∈J,∑j∈Jyj=p,xij,yj∈0,1,∀i∈I,j∈J,where I is the customer set. Denote J as the set of alternative charging stations. The decision variables are defined as follows:(19)yj=1,ifjis selected as the charging station,0,otherwise,xij=1,if customeriis served by charging stationj0,otherwise.In addition, the different vehicles have distinct energy consumption and carbon emissions. So papers [6, 24] are referenced to determine the relevant parameters of vehicle energy consumption. The CO2 emission factor εk (according to the Intergovernmental Panel on Climate Change) is 0.785 kg/kWh for EVs and 3.096 kg/L for CVs. Other relevant parameters are shown in Table 4.Table 4
Values of parameters in the model.
NotationValueEVsCVsfk300230ck0.855.08εk0.7853.096ρkxρkx=0.5+0.7/80xρkx=0.175+0.225/80xce0.8co2pe0.5pl5Y80g240Q80All experiments in this paper are based on MATLAB R2018b and implemented on a computer with a Core (TM) i5 processor (1.6 GHz) and a memory of 16 GB. The relevant parameters are set as follows: population sizepopsize=100, maximum number of iterations IterMax=500, crossover probability pc=0.9, and mutation probability pm=0.1.The Pareto front obtained by 500 iterations of the algorithm for the three instances is shown in Figure5. It is easy to know that the Pareto front slopes to the lower right; that is, the value of f1 is getting smaller, while the value of f2 is getting larger. All the solutions are not dominated by each other, which indicates that there is no positive correlation between the two objectives and they are contradictory to some extent.Figure 5
Pareto front. (a) Instance I. (b) Instance II. (c) Instance III.
(a)(b)(c)


### 4.2. Comparative Analysis and Discussion

In order to further verify the effectiveness of the NSGA-II, this paper firstly optimizesf1 and f2, respectively, by the classical genetic algorithm. The rationality of the Pareto solution set is further illustrated by comparing it with the optimal value of each objective.Taking Instance I as an example, the optimization off1 is carried out and Figure 6(a) is obtained. The left picture shows that the optimal value changes with the number of iterations during the operation of the algorithm, which indicates that the algorithm is valid. The right picture is the routing of vehicles corresponding to the optimal solution (the triangles represent the distribution center, the dots denote the customers, and the squares are the symbol of the charging stations; the EVs are represented by dashed lines, while the solid lines represent the CVs). Figure 6(b) can be obtained by performing the same operation on f2. The comparison of single-objective optimization and the Pareto solution set is shown in Table 5.Figure 6
Single-objective optimization results of Instance I. (a) The results when only optimizingf1. (b) The results when only optimizing f2. (c) The results when optimizing f1+f2.
(a)(b)(c)Table 5
Comparison of single-objective optimization and Pareto solution set for Instance I.
Solutionf1f21The solution with the smallestf1 in the Pareto solution set2469.860674.201The optimal solution obtained when only optimizingf12389.4141051.8052The solution with the smallestf2 in the Pareto solution set3163.406360.986The optimal solution obtained when only optimizingf23853.712217.9993The solution with the smallestf1+f2 in the Pareto solution set2483.145547.795The optimal solution obtained when optimizingf1+f22483.145547.795As can be seen from Table5, when only optimizing f1, the value of optimal solution is 2389.414, which is better than all the solutions in Pareto solution set. However, the value of f2 corresponding to this optimal solution is 1051.805, which is worse than all the solutions in Pareto solution set. The same goes for the situation when f2 is optimized. Therefore, if one objective is optimized separately, its optimal value is better than all the solutions of the Pareto solution set, while the other objective is worse than all the results of the Pareto solution set. This shows that the bi-objective optimization results obtained by the NSGA-II are reasonable.On the basis of analysis about rationality, this paper compares with the linear weighting method to further analyze the validity of the results. Sincef1 and f2 have the same dimension, they are given the same weight.f1+f2 can be regarded as one objective by adding directly, and its optimal solution can be obtained through genetic algorithm. The result is shown in Figure 6(c).It can be concluded from the results in Table5 that the value of the optimal solution obtained by the linear weighting method is 3030.94 (the corresponding values of f1 and f2 are 2483.145 and 547.795, respectively), and it is equal to the smallest value of f1+f2 in the Pareto solution set. Any solution in the Pareto solution set is nondominated by the optimal solution of the linear weighting method, which indicates the effectiveness of NSGA-II.The same operation is performed for Instances II and III and the results in Table6 further verify the effectiveness of the NSGA-II.Table 6
Comparison of single-objective optimization and Pareto solution set for Instance II and Instance III.
<!C Cnt:6> SolutionInstance IIInstance IIIf1f2f1f21The solution with the smallestf1 in the Pareto solution set4869.1042219.8369107.5422137.740The optimal solution obtained when only optimizingf14619.7592508.6188538.1444890.5662The solution with the smallestf2 in the Pareto solution set5272.350771.92010331.6271801.689The optimal solution obtained when only optimizingf27989.809539.33214335.2381100.9253The solution with the smallestf1+f2 in the Pareto solution set5093.728818.2279120.6502064.374The optimal solution obtained when optimizingf1+f25093.728818.2279120.6502064.374


## 4.1. Computational Results

In this section, the effectiveness of the algorithm is illustrated by three numerical instances. Referring to [11], this paper sets the number of charging stations as ⌊0.1⋅N⌋ (N is the number of customers). The instances used are generated based on the Solomon benchmark set, as shown in Table 3. Meanwhile, we assume that the number of vehicles is sufficient, so set the number of vehicles equal to the number of customers.Table 3
Number of customers and charging stations.
Instance IInstance IIInstance IIINumber of customers255091Number of charging stations259Since the final result will be affected by the location of the charging station, it is necessary to select the reasonable charging station locations by a location model.(18)min∑i∈I∑j∈Jdijxij,s.t.∑j∈Jxij=1,∀i∈I,xij≤yj,∀i∈I,j∈J,∑j∈Jyj=p,xij,yj∈0,1,∀i∈I,j∈J,where I is the customer set. Denote J as the set of alternative charging stations. The decision variables are defined as follows:(19)yj=1,ifjis selected as the charging station,0,otherwise,xij=1,if customeriis served by charging stationj0,otherwise.In addition, the different vehicles have distinct energy consumption and carbon emissions. So papers [6, 24] are referenced to determine the relevant parameters of vehicle energy consumption. The CO2 emission factor εk (according to the Intergovernmental Panel on Climate Change) is 0.785 kg/kWh for EVs and 3.096 kg/L for CVs. Other relevant parameters are shown in Table 4.Table 4
Values of parameters in the model.
NotationValueEVsCVsfk300230ck0.855.08εk0.7853.096ρkxρkx=0.5+0.7/80xρkx=0.175+0.225/80xce0.8co2pe0.5pl5Y80g240Q80All experiments in this paper are based on MATLAB R2018b and implemented on a computer with a Core (TM) i5 processor (1.6 GHz) and a memory of 16 GB. The relevant parameters are set as follows: population sizepopsize=100, maximum number of iterations IterMax=500, crossover probability pc=0.9, and mutation probability pm=0.1.The Pareto front obtained by 500 iterations of the algorithm for the three instances is shown in Figure5. It is easy to know that the Pareto front slopes to the lower right; that is, the value of f1 is getting smaller, while the value of f2 is getting larger. All the solutions are not dominated by each other, which indicates that there is no positive correlation between the two objectives and they are contradictory to some extent.Figure 5
Pareto front. (a) Instance I. (b) Instance II. (c) Instance III.
(a)(b)(c)


## 4.2. Comparative Analysis and Discussion

In order to further verify the effectiveness of the NSGA-II, this paper firstly optimizesf1 and f2, respectively, by the classical genetic algorithm. The rationality of the Pareto solution set is further illustrated by comparing it with the optimal value of each objective.Taking Instance I as an example, the optimization off1 is carried out and Figure 6(a) is obtained. The left picture shows that the optimal value changes with the number of iterations during the operation of the algorithm, which indicates that the algorithm is valid. The right picture is the routing of vehicles corresponding to the optimal solution (the triangles represent the distribution center, the dots denote the customers, and the squares are the symbol of the charging stations; the EVs are represented by dashed lines, while the solid lines represent the CVs). Figure 6(b) can be obtained by performing the same operation on f2. The comparison of single-objective optimization and the Pareto solution set is shown in Table 5.Figure 6
Single-objective optimization results of Instance I. (a) The results when only optimizingf1. (b) The results when only optimizing f2. (c) The results when optimizing f1+f2.
(a)(b)(c)Table 5
Comparison of single-objective optimization and Pareto solution set for Instance I.
Solutionf1f21The solution with the smallestf1 in the Pareto solution set2469.860674.201The optimal solution obtained when only optimizingf12389.4141051.8052The solution with the smallestf2 in the Pareto solution set3163.406360.986The optimal solution obtained when only optimizingf23853.712217.9993The solution with the smallestf1+f2 in the Pareto solution set2483.145547.795The optimal solution obtained when optimizingf1+f22483.145547.795As can be seen from Table5, when only optimizing f1, the value of optimal solution is 2389.414, which is better than all the solutions in Pareto solution set. However, the value of f2 corresponding to this optimal solution is 1051.805, which is worse than all the solutions in Pareto solution set. The same goes for the situation when f2 is optimized. Therefore, if one objective is optimized separately, its optimal value is better than all the solutions of the Pareto solution set, while the other objective is worse than all the results of the Pareto solution set. This shows that the bi-objective optimization results obtained by the NSGA-II are reasonable.On the basis of analysis about rationality, this paper compares with the linear weighting method to further analyze the validity of the results. Sincef1 and f2 have the same dimension, they are given the same weight.f1+f2 can be regarded as one objective by adding directly, and its optimal solution can be obtained through genetic algorithm. The result is shown in Figure 6(c).It can be concluded from the results in Table5 that the value of the optimal solution obtained by the linear weighting method is 3030.94 (the corresponding values of f1 and f2 are 2483.145 and 547.795, respectively), and it is equal to the smallest value of f1+f2 in the Pareto solution set. Any solution in the Pareto solution set is nondominated by the optimal solution of the linear weighting method, which indicates the effectiveness of NSGA-II.The same operation is performed for Instances II and III and the results in Table6 further verify the effectiveness of the NSGA-II.Table 6
Comparison of single-objective optimization and Pareto solution set for Instance II and Instance III.
<!C Cnt:6> SolutionInstance IIInstance IIIf1f2f1f21The solution with the smallestf1 in the Pareto solution set4869.1042219.8369107.5422137.740The optimal solution obtained when only optimizingf14619.7592508.6188538.1444890.5662The solution with the smallestf2 in the Pareto solution set5272.350771.92010331.6271801.689The optimal solution obtained when only optimizingf27989.809539.33214335.2381100.9253The solution with the smallestf1+f2 in the Pareto solution set5093.728818.2279120.6502064.374The optimal solution obtained when optimizingf1+f25093.728818.2279120.6502064.374


## 5. Conclusion

Aiming at the coexistence of electric vehicles and conventional vehicles in the fleet, a new vehicle routing problem with mixed fleet and soft time window considering recharging is proposed, and a bi-objective mathematical programming model is established. Secondly, a novel coding and decoding method is designed according to the character of the problem, and the solution of the problem is obtained by NSGA-II. In order to illustrate the effectiveness of this algorithm, this paper further compares it with the single-objective optimization and the linear weighting method. It shows that the results obtained by the NSGA-II are reasonable, and two objectives are contradictory to some extent. This paper can also provide management insights into the promotion of electric vehicles in other areas.There are several interesting avenues for future research. For example, a more realistic modeling of the recharging process, the reduction of the battery capacity due to aging, and pickup and delivery simultaneously are all interesting topics.



---

*Source: 9086229-2021-09-06.xml*