Electric Vehicle Tour Planning Considering Range Anxiety

Abstract

In this study, the tour planning problem for electric vehicles is investigated. We aim to derive the optimal route and thus, to maximize profitability and minimize range anxiety within the time horizon. To solve this problem, a bi-objective mixed integer model is proposed. Specifically, we first introduced the reliability of route planning and quantified it as a cost with specific functions. The nonlinear model was then converted into a bi-objective mixed integer linear program, and an interactive branch and bound algorithm was adopted. Numerical experiments conducted on different networks have shown that the model that considers range anxiety offers more effective solutions. This means that our model is able to plan the routes with high reliability and low risk of profit loss and accidents.

1. Introduction

Currently, about 86% of global primary energy demand depends on fossil fuels, in which coal, gas and oil account for 23%, 27% and 36%, respectively. Based on current information about reserves and daily production, oil will be depleted in 2066 [1]. Promoting alternative fuel vehicles is an effective solution to this problem since gasoline for vehicles is a major oil product (Gasoline takes up to 46% of US oil products (U.S. Energy Information Administration: https://www.eia.gov/petroleum.)). Completely powered by electricity, battery electric vehicles (BEV) are fourfold more energy efficient than internal combustion engine vehicles (ICEV) [2] and they offer opportunities for the use of renewable energies. This greatly contributes to sustainable development because electric vehicle (EV) technologies facilitate a transportation sector powered by renewable energy in the meantime. For example, in the field of logistics, EVs are considered a valid alternative to ICEVs [3].

Despite the advantages compared to the ICEVs, the EV market share in most countries by 2016 was still low. In 2016, 95% of EV sales took place in 10 countries and only six of them had an EV market share over 1% [4]. These were Norway, The Netherlands, Sweden, France, UK and China. One of the main barriers to promoting EVs is known as the range anxiety (RA) of EV drivers (Esmaili et al. [5], Xie et al. [6]), because of the narrow range compared to ICEVs and the lack of charging facilities. In recent years, the EV market share has increased rapidly with the promotion policies of governments. However, battery technologies have failed to keep pace with other areas of EV innovation. Therefore, RA remains a major impediment to the marketability of EVs.

Noel et al. [7] define RA as the psychological anxiety a consumer experiences in response to the limited range of electric vehicles. Wang et al. [8] describe RA as the mental distress, or fear, of being stranded on roads should the battery run out. They imply that even if a completed charging network exists, EV drivers could still suffer from RA. RA is expressed not only on the cognitive level (e.g., concerns about running out of energy and not being able to reach the destination) but also on the emotional, behavioral and physiological level [9]. Thus, RA could spoil a planned trip or even greatly increase the driving risk. Consequently, the effect of RA must be considered in the planning of EV trips.

This study aims to derive the optimal route for the electric vehicle tour planning (EVTP) problem so that profitability is maximized and range anxiety is minimized within the time horizon. We propose a new bi-objective model with time windows and range limitations considering RA of EV drivers. RA cost is quantified by proposed functions to evaluate the reliability of planned routes and an interactive branch and bound algorithm is adopted. Our model is intended to allow routes to be planned with less range anxiety and, thus, to render EVs more saleable in future.

The main contributions of this paper are as follows:

• We first define the reliability of an electric vehicle tour plan with range anxiety and propose a function form to quantify it.
• The EVTP, considering range anxiety (EVTPRA), is modeled as a non-linear bi-objective MIP with the full-recharge and partial-recharge policies.
• We linearize the nonlinear terms in the model and describe an exact algorithm, Interactive Branch and Bound Algorithm (IBBA), based on non-dominance sets to find the optimal solutions of the EVTPRA.
• We conduct numerical experiments on two sets of benchmark networks derived from well-known instances in previous literature and show the model considering range anxiety is able to plan the routes with high reliability and low risk of profit loss and accidents.

The rest of the paper is organized as follows. Section 2 gives a brief overview of related literature in the area of EVRP and RA for EV drivers. Section 3 describes the basic model and the solution approach with the proposed algorithm is reported in Section 4. Section 5 reports the computational results. Some conclusions are drawn in Section 6.

2. Literature Review

This study involves two related area of research, which are the tour planning problem and the RA factor of drivers.

2.1. Tour Planning

One of the well-known tour panning problems is the traveling salesman problem (TSP), which, as well as its multitudinous extensions, is one of the most widely studied combinatorial optimization problems of visiting attractions from a central depot. A considerable number of papers and books have been published to tackle these problems [10]. However, it is not often possible for tourists to visit all of the points of interests (POIs) within a short time period. This motivates researchers to study TSPs with profit, which is known as the orienteering problem (OP) when a desirable path is preferred rather than a circuit. To consider the time factor of the tour system, the OP can be extended to time-dependent (TD) OP with time windows (TW) (see

Table 1. Existing studies on TSP and OP.

Literature	Problem				Algorithm
	TSP	OP	TW	TD
Picard and Queyranne [11]	◯			◯	Branch and Bound
Baker [12]	◯		◯		Branch and Bound
Dumas et al. [13]	◯		◯		Dynamic Programming
Albiach et al. [14]	◯		◯	◯	Branch and Bound
Righini and Salani [15]		◯	◯		Dynamic Programming
Montemanni and Gambardella [16]		◯	◯		Ant Colony System
Vansteenwegen et al. [17]		◯	◯		Local Search
Abbaspour and Samadzadegan [18]		◯		◯	Genetic Algorithms
Garcia et al. [19]		◯	◯	◯	Local Search
Gavalas et al. [20]		◯	◯	◯	Cluster-based Heuristics

).

As far as we know, these models do not consider the range factor of vehicles, because ICEVs usually have a large range and their refueling time is relatively short compared to EVs. Therefore, this assumption no longer applies to electric vehicle routing problem (EVRP). One of the seminal papers on this topic is by Schneider et al. [21]. They introduce the electric vehicle routing problem with time windows and recharging stations (EVRPTW) and present a hybrid heuristic with high performance. In their model, vehicles have a battery capacity and they can be recharged at a specified rate. Ever since, EVRPTW has been extended by considering energy consumption, battery capacity, vehicle types, charging facility location, etc. and many algorithms have been proposed to deal with this NP-hard problem (see

Table 2. Existing studies on EVRPTW.

Literature	Feature	Algorithm
		Exact	Heuristic
Preis et al. [22]	Load-dependent Energy Consumption		◯
Bruglieri et al. [23]	Travel and Waiting Time		◯
Goeke and Schneider [24]	Mixed Fleet with ICEVs		◯
Hiermann et al. [25]	Heterogeneous Electric Fleet	◯	◯
Desaulniers et al. [26]	Recharging Policy	◯
Keskin and Çatay [27]	Recharging Policy		◯
Schiffer and Walther [28]	Location Routing		◯
Paz et al. [29]	Multi-Depot and Location Routing	◯

).

By ignoring capacity constraints and considering one vehicle only, the electric traveling salesman problem with time windows (ETSPTW) is defined by Roberti and Wen (2016) [3]. They propose a mixed integer linear formulation that can solve 20-customer instances in short computing times and a Three-Phase Heuristic algorithm based on General Variable Neighborhood Search and Dynamic Programming. Küçükoğlu et al. [30] extend this model by considering charging operations at customer locations with different charging rates. They introduce a new and effective hybrid Simulated Annealing and Tabu Search algorithm to obtain 26 new best results for the ETSPTW instances. Similarly, Wang et al. [31] describe the electric vehicle orienteering problem with time windows (EVOPTW), which considers the recharging mode as battery swapping. It always recharges the battery to the full state of charging (SOC) and takes fixed time. However, to the best of our knowledge, none of those EVRP models consider the RA factor of EV drivers.

2.2. Range Anxiety

As one of the main barriers to promote EVs, RA has been considered in many previous studies. Neubauer and Wood [32] examine how BEV use is affected by RA in various charging facility situations, including variable time schedules, power levels, and locations. Huang et al. [33] propose two integer program (IP) models optimally to locate the public charging stations with two different charging modes. Jiao et al [34] introduce a mixed integer program (MIP) model considering imbalance and RA to address the location problem with the allocation of two types of charging piles.

Guo et al. [35] formulate a bi-level integer programming model based on a flow decay function with range anxiety parameters drawn from surveys and statistical analysis on the basis of local conditions. It can be used to determine when a user’s level of anxiety changes as the remaining battery capacity is seen to be less than the range anxiety threshold. These works principally focus on the planning of charging facilities to mitigate the impact on the market penetration of EVs.

The literature that quantitatively addresses RA for EV is comparatively little, and one of which is proposed to describe the cost paid to EV owners for not fully recharging their batteries in energy management of microgrids by Esmaili et al. [5]. In their study, the RA cost is counted in timeslots, which is not suitable for the driving EVs. This motivates us to propose a new function form for the RA cost to help plan the EV tour.

To conclude, the research gap observed is that (1) None of those EVRP models consider the range anxiety factor of drivers, (2) No suitable functions for range anxiety cost have been proposed to address this problem. This paper makes an original contribution by addressing both of these points to fill the gap.

3. Basic Model

The model can simulate the spatial-temporal behaviors of EV tourists, and it can be used to find the feasible routes to visit the most POIs with time windows in a given time horizon. Denote a directed transportation network as $\mathcal{G}=(\mathcal{N},\mathcal{A})$, where $\mathcal{N}$ is the set of nodes with element $n\in \mathcal{N}$ consisting of attraction nodes (${\mathcal{N}}_S$), recharging nodes ${\mathcal{N}}_C$, and $\mathcal{A}$ is the set of directed links with elements $ij\in \mathcal{A}$. The mathematical notations used in the paper are listed in

Table A1. Table of notations.

Sets and Parameters
${\mathcal{N}}_S$	The set of attraction nodes
${\mathcal{N}}_C$	The set of recharging nodes
$\mathcal{N}$	The set of all nodes with origin o and destination d, $\mathcal{N}={\mathcal{N}}_S\cup {\mathcal{N}}_C$
$\mathcal{A}$	The set of directed links with element $ij$
$\mathcal{G}$	A directed transportation network, $\mathcal{G}=\left(\mathcal{N},\mathcal{A}\right)$
$w_j$	The weight score of visiting node j
$s_j$	The duration of stay at attraction node j
$l_j$	The starting time of service at node j
$u_j$	The ending time of service at node j
$g_j$	The recharging rate at node j
$h_j$	Binary, 1 if node j is a designated battery swap station, and 0 otherwise
${\tau }_{ij}$	The travel time from node i to node j
$e_{ij}$	The energy consumption by travelling from node i to node j
Q	The vehicle energy capacity
T	The limit of a tour or trip duration or the total allowed duration
Decision Variables
$t_j$	The arrival time at node j
$c_j$	The duration of stay at recharging node j
$q_j$	The energy level on arrival at node j
$r_j$	The amount of energy recharged at node j
$y_j$	Binary variable, 1 if node j has been visited, and 0 otherwise
$v_j$	Binary variable, 1 if the vehicle recharges at node i, and 0 otherwise
$x_{ij}$	Binary variable, 1 if link $ij$ is on the path of the pair $od$, and 0 otherwise

. Many assumptions are made to formulate the model, including a single type of EV with a constant range, the fixed travel time and energy consumption between nodes, and a single visit to each of the recharging stations. Additional assumptions and notations are defined as follows, which will be used throughout this paper.

3.1. Assumptions

The following assumptions are made in this study to model the problem:

• a vehicle at node i can arrive at node j if and only if its remaining energy is sufficient to cover the distance between nodes i and j.
• the travel time and energy consumption on each link are fixed and known.
• the recharging time is proportionally linear to the desired quantity to recharge with respect to an inverse recharging rate $g_j$ for node j.

3.2. Recharging Policy

In the literature, two policies are usually considered to determine the amount of battery recharged at each stop: full-recharge policy and partial-recharge policy. With the full-recharge policy, the battery is always fully recharged, while with the partial-recharge policy, any quantity can be recharged at each stop as long as the energy level does not exceed the capacity of the battery [3]. The energy level constraints with the full-recharge policy can be formulated as:

$$\begin{array}{ccc}\hfill & q_j\ge \left(q_i-e_{ij}{x}_{ij}\right)+r_i-\left(1-x_{ij}\right)Q\hfill & \hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & q_j\le \left(q_i-e_{ij}{x}_{ij}\right)+r_i+\left(1-x_{ij}\right)Q\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & r_i\le Q-q_i\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill & r_i\le Q{v}_i\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & r_i\ge Q{v}_i-q_i\hfill \end{array}$$

(5)

Constraints (1)–(2) describe the relationship of energy levels among all nodes. Constraints (3)–(5) restrict the recharging energy to the battery capacity. These constraints can be used to connect energy levels with the partial-recharge policy as well by eliminating Constraints (5) (For the rest of this paper, the models are formulated with the partial-recharge policy, and it is effortless to convert them into models with the full-recharge policy by only adding Constraints (5)).

Proposition 1.

The partial-recharge policy is at least as desirable as the full-recharge policy.

Proof. 

Since the partial-recharge policy can be modeled as the formulation of the full-recharge policy without Constraints (5), the model with the partial-recharge formulation is apparently a relaxtion of that with the full-recharge formulation.  □

The duration of stay at attraction nodes is fixed, while that at recharging nodes depends on how much the energy is recharged. Therefore, the time constraints for attraction nodes and recharging nodes are slightly different. For the attraction nodes, it can be presented as follows:

$$\begin{array}{ccc}\hfill & \left\{\begin{array}{cc}{t}_j\ge \left(t_i+{\tau }_{ij}{x}_{ij}\right)+y_i{s}_i-\left(1-x_{ij}\right)T\hfill & \forall i\in {\mathcal{N}}_S,ij\in \mathcal{A}\hfill \\ t_j\le \left(t_i+{\tau }_{ij}{x}_{ij}\right)+y_i{s}_i+\left(1-x_{ij}\right)T\hfill & \forall i\in {\mathcal{N}}_S,ij\in \mathcal{A}\hfill \end{array}\right.\hfill & \hfill \end{array}$$

(6)

and for the recharging nodes, the following constraints for the duration of stay have to be added.

$$\begin{array}{ccc}\hfill & \left\{\begin{array}{cc}{t}_j\ge \left(t_i+{\tau }_{ij}{x}_{ij}\right)+y_i{c}_i-\left(1-x_{ij}\right)T\hfill & \forall i\in {\mathcal{N}}_C,ij\in \mathcal{A}\hfill \\ t_j\le \left(t_i+{\tau }_{ij}{x}_{ij}\right)+y_i{c}_i+\left(1-x_{ij}\right)T\hfill & \forall i\in {\mathcal{N}}_C,ij\in \mathcal{A}\hfill \\ c_i=g_i{r}_i\hfill & \forall i\in {\mathcal{N}}_C\hfill \end{array}\right.\hfill & \hfill \end{array}$$

(7)

Constraints (6) and (7) reflect the relationship of arrival time among attraction nodes and recharging nodes, respectively.

3.3. Range Anxiety Cost Function

Previous studies have emphasized profit and efficiency. However, our study also considers the reliability of route planning. This means that we consider (1) a driver’s ability to complete a journey without the unplanned changes (in response to traffic, accidents and other variables) that characterize real-life transportation and (2) driving risk. We quantify reliability with the RA cost during an entire trip. Lower RA costs result in higher average SOC in EVs. Thus, drivers experience less RA during a trip. Considering the link $ij$ with travel time ${\tau }_{ij}$ and energy consumption $e_{ij}$, we assume that energy is consumed linearly at a fixed rate. An RA function proposed by Esmaili et al. [5] calculates the payment to EV owners in return for partial charging in the microgrids. Based on this idea, the cost attributed to RA in planning a trip proposed by us can be presented as:

$$\begin{array}{c}\hfill x_{ij}{\int }_0^{{\tau }_{ij}}k\left(Q-\left(q_j+e_{ij}-\frac{e_{ij}}{{\tau }_{ij}}t\right)\right)dt=k\left(Q-q_j-\frac{e_{ij}}{2}\right){\tau }_{ij}{x}_{ij}\end{array}$$

(8)

where k is the per kWh value of RA. In fact, this function could also quantify the RA cost even when the energy consumes nonlinearly. A comparatively intuitive way to tackle this problem is to piecewise linearize the energy consumption curve and apply this function afterwards. The RA cost can, alternatively, be limited by constraints, while choosing an appropriate value for it might be questionable.

3.4. Nonlinear Model

With the RA cost function, the problem can be formulated as the following bi-objective mixed integer nonlinear programming.

$$\begin{array}{cc}\hfill \mathrm{Maximize}& \sum _{j\in {\mathcal{N}}_S}{w}_j{y}_j\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \mathrm{Minimize}& \sum _{ij\in \mathcal{A}}k\left(Q-q_j-\frac{e_{ij}}{2}\right){\tau }_{ij}{x}_{ij}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \mathrm{Subject}\mathrm{to}:& \sum _i{x}_{ij}\ge y_j,\forall j\in \mathcal{N}\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & \sum _i{x}_{ij}\le 1,\forall j\in \mathcal{N}\backslash \{o,d\},ij\in \mathcal{A}\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & \sum _i{x}_{ij}-\sum _i{x}_{ji}=\left\{\begin{array}{cc}-1\hfill & j=o,\hfill \\ 1\hfill & j=d,\hfill \\ 0\hfill & \mathrm{o}/\mathrm{w},\hfill \end{array}\right.\forall j\in \mathcal{N},ij\in \mathcal{A}\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & q_j\le \left(\sum _i{x}_{ij}\right)Q,\forall j\in \mathcal{N}\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & t_j\le \left(\sum _i{x}_{ij}\right)T,\forall j\in \mathcal{N}\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & v_j\le h_j{y}_j,\forall j\in \mathcal{N}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & l_j\le t_j\le u_j,\forall j\in \mathcal{N}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & t_j,c_j,q_j,r_j\ge 0,\forall j\in \mathcal{N}\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & y_j,v_j,x_{ij}\in \{0,1\},\forall j\in \mathcal{N},ij\in \mathcal{A}\hfill \end{array}$$

(19)

Constraints (1)–(4), (6) and (7).

Constraints (11) and (12) handle the connectivity of visits to nodes. Constraints (13) establish flow conservation at each node. Constraints (14) and (15) guarantee the energy and time feasibility for nodes, respectively. Constraints (16) restrict the recharging task to only visited recharging nodes. Constraints (17) enforce that every node is visited within its time window. Constraints (18) ensure that the value of the energy and time never fall below 0 and Constraints (19) reflect the binary nature of decisions.

4. Solution Approach

It is noted that the bilinear terms (i.e., $q_j{x}_{ij}$ in the objective function and $y_i{c}_i$ in Constraints (7)) spoil the mathematical property of linearity and make this model much more difficult to solve. In this section, the reformulation-linearization technique (RLT) is applied to convert the nonlinear terms into equivalent linear constraints. Consequently, the original problem is then transformed into a mixed-integer linear program and a global optimization algorithm (e.g., an interactive branch and bound algorithm) can be employed to obtain the optimal solution.

As the bilinear terms exist in both the objective function and constraints, we have to linearize them in different ways. The linearization is exhibited in Section 4.1 and Section 4.2, respectively. The interactive branch and bound algorithm basically divides the intervals of nondominance into smaller intervals. Each nondominated solution is obtained by solving a single-objective mixed integer program (SMIP). The detail of the algorithm is described in Section 4.3 and Section 4.4.

4.1. Linearization of Objective Function

In the proposed model, the objective function of RA cost consists of bilinear terms, which are the product of the continuous variable $q_j$ and the binary variable $x_{ij}$. By applying the RLT, these nonlinear terms can thus be converted into an equivalent set of linear constraints [36]. Particularly, a new variable $f_{ij}$ is first introduced to express the bilinear term as follows:

$$\begin{array}{ccc}\hfill & f_{ij}=q_j{x}_{ij},\forall ij\in \mathcal{A}\hfill & \hfill \end{array}$$

(20)

Moreover, the following linear constraints are added to achieve the equivalence to previous bilinear terms:

$$\begin{array}{ccc}\hfill & \left\{\begin{array}{cc}0\le f_{ij}\le Q{x}_{ij}\hfill & \forall ij\in \mathcal{A}\hfill \\ q_j-Q\left(1-x_{ij}\right)\le f_{ij}\le q_j\hfill & \forall ij\in \mathcal{A}\hfill \end{array}\right.\hfill & \hfill \end{array}$$

(21)

The equivalence between terms $q_j{x}_{ij}$ and terms $f_{ij}$ with Constraints (21) can be demonstrated by testing the two possible values of the binary variable $x_{ij}$. If $x_{ij}=0$, according to Equation (20), we know $f_{ij}=0$. Then Constraints (21) are equivalent to:

$$\begin{array}{ccc}\hfill & \left\{\begin{array}{cc}0\le f_{ij}\le 0\hfill & \forall ij\in \mathcal{A}\hfill \\ q_j-Q\le f_{ij}\le q_j\hfill & \forall ij\in \mathcal{A}\hfill \end{array}\right.\hfill & \hfill \end{array}$$

This only holds when $f_{ij}=0$. If $x_{ij}=1$, $f_{ij}=q_j$. Then Constraints (21) are equivalent to:

$$\begin{array}{ccc}\hfill & \left\{\begin{array}{cc}0\le f_{ij}\le Q\hfill & \forall ij\in \mathcal{A}\hfill \\ q_j\le f_{ij}\le q_j\hfill & \forall ij\in \mathcal{A}\hfill \end{array}\right.\hfill & \hfill \end{array}$$

This only holds when $f_{ij}=q_j$. Therefore, we verified the equivalence between terms $q_j{x}_{ij}$ and terms $f_{ij}$ with Constraints (21).

4.2. Linearization of Constraints

The Contraints (7) similarly comprises bilinear terms produced by the continuous variable $c_i$ and the binary variable $y_i$. We first reformulate Constraints (7) by substituting the nonlinear term $y_i{c}_i$ with the linear term $c_i$:

$$\begin{array}{ccc}\hfill & \left\{\begin{array}{cc}{t}_j\ge \left(t_i+{\tau }_{ij}{x}_{ij}\right)+c_i-\left(1-x_{ij}\right)T\hfill & \forall i\in {\mathcal{N}}_C,ij\in \mathcal{A}\hfill \\ t_j\le \left(t_i+{\tau }_{ij}{x}_{ij}\right)+c_i+\left(1-x_{ij}\right)T\hfill & \forall i\in {\mathcal{N}}_C,ij\in \mathcal{A}\hfill \end{array}\right.\hfill & \hfill \end{array}$$

(22)

Furthermore, the following linear constraints have to be added to equalize the original formulation.

$$\begin{array}{ccc}\hfill & \left\{\begin{array}{cc}{c}_i\le g_i{r}_i+g_{i}Q\left(1-y_i\right)\hfill & \forall i\in {\mathcal{N}}_C\hfill \\ c_i\ge g_i{r}_i-g_{i}Q\left(1-y_i\right)\hfill & \forall i\in {\mathcal{N}}_C\hfill \\ c_i\le g_{i}Q{y}_i\hfill & \forall i\in {\mathcal{N}}_C\hfill \end{array}\right.\hfill & \hfill \end{array}$$

(23)

The equivalence between Constraints (7) and Constraints (22) and (23) can be demonstrated by testing the two possible values of the binary variable $y_i$. If $y_i=0$, the Constraints (23) are equivalent to:

$$\begin{array}{ccc}\hfill & \left\{\begin{array}{cc}{c}_i\le g_i{r}_i+g_{i}Q\hfill & \forall i\in {\mathcal{N}}_C\hfill \\ c_i\ge g_i{r}_i-g_{i}Q\hfill & \forall i\in {\mathcal{N}}_C\hfill \\ c_i\le 0\hfill & \forall i\in {\mathcal{N}}_C\hfill \end{array}\right.\hfill & \hfill \end{array}$$

As $c_i\ge 0$, we know $c_i=0$. If $y_i=1$, the Constraints (23) are equivalent to:

$$\begin{array}{ccc}\hfill & \left\{\begin{array}{cc}{c}_i\le g_i{r}_i\hfill & \forall i\in {\mathcal{N}}_C\hfill \\ c_i\ge g_i{r}_i\hfill & \forall i\in {\mathcal{N}}_C\hfill \\ c_i\le g_{i}Q\hfill & \forall i\in {\mathcal{N}}_C\hfill \end{array}\right.\hfill & \hfill \end{array}$$

We know $c_i=g_i{r}_i$. Therefore, we verified the equivalence between Constraints (7) and Constraints (22) and (23).

4.3. Bi-Objective Mixed Integer Linear Programming

The multi-objective problems can be usually formulated as single-objective problems with additional constraints by examining the preferences of the decision makers (DMs). For the EVTP, the majority of DMs are more concerned about the total profit of the trip (i.e., the summed weight of POIs), which can be presented as:

$$\begin{array}{ccc}\hfill \mathrm{Maximize}{z}_1=& \sum _{j\in {\mathcal{N}}_S}{w}_j{y}_j\hfill & \hfill \end{array}$$

(24)

However, for some specific industries or special circumstances, DMs could care more about the impact on EV drivers rather than only the profit, because it is related to the reliability of the trip (e.g., the driving behavior and psychological pressure). This objective can be presented as:

$$\begin{array}{ccc}\hfill \mathrm{Maximize}{z}_2=& \sum _{ij\in \mathcal{A}}k{\tau }_{ij}\left(f_{ij}+\left(\frac{e_{ij}}{2}-Q\right)x_{ij}\right)\hfill & \hfill \end{array}$$

(25)

Let $m,n\in \{1,2\},m\ne n$, where $z_m$ is the objective that the DMs are more concerned about. For discrete problems, the optimal solution is usually not unique. We denote $z_{max}$ the maximum value of $z_n$ without considering the objective $z_m$. Then the globally optimal solution, which is found by solving the proposed model with the constraint $z_n\ge z_{max}$, dominates all the other solutions. The single-objective problem (SP), when the partial-recharge policy is applied, can be formulated as follows:

$$\begin{array}{ccc}\hfill \mathrm{Maximize}& z_m\hfill & \hfill \\ \hfill \mathrm{Subject}\mathrm{to}:& z_n\ge z_{max}\hfill \end{array}$$

Constraints (1)–(4), (6), (11)–(19) and (21)–(23).

The SP model is apparently a mixed integer linear program (MILP), which can be easily solved by optimization solvers (e.g., Cplex and Gurobi). However, for those DMs whose preferences are consistent, transitive and invariant over the process, multi-objective mixed integer program (MMIP) with an interactive method is more useful in the real world [37]. The bi-objective problem (BP), when the partial-recharge policy is applied, can be formulated as:

$$\begin{array}{ccc}\hfill \mathrm{Maximize}& \{z_1,z_2\}\hfill & \hfill \end{array}$$

(26)

Subject to: Constraints (1)–(4), (6), (11)–(19) and (21)–(23).

Proposition 2.

Let $(z_1^l,z_2^u)$ and $(z_1^u,z_2^l)$ denote the optimal solution for problem SP with $m=1,2$, respectively. Then the intervals of nondominance for $z_1$ and $z_2$ for problem BP are given by $[z_1^l,z_1^u]$ and $[z_2^l,z_2^u]$.

Proof. 

For convenience, for any two vectors x and y of the same dimension, $x⪯y$ will denote that $x\le y$ and $x\ne y$. Suppose there exists a nondominated solution $(z_1^{*},z_2^{*})$. As $(z_1^l,z_2^u)$ and $(z_2^l,z_1^u)$ are the optimal solutions for problem SP, Obviously, $z_1^{*}\le z_1^u,z_2^{*}\le z_2^u$. If $z_1^{*}<z_1^l$, it contradicts to the nondominance as $(z_1^{*},z_2^{*})⪯(z_1^l,z_2^u)$. If $z_2^{*}<z_2^l$, it contradicts to the nondominance as $(z_1^{*},z_2^{*})⪯(z_1^u,z_2^l)$. This completes the proof.  □

4.4. Algorithm Design

Based on Proposition 2 and the method proposed by Aksoy [37], we present the solution algorithm as shown in Algorithm 1. The proposed interactive branch-and-bound algorithm obtains nondominated solutions by dividing the intervals of nondominance into smaller intervals. Interactively selecting with DMs, it exhaustively explores to achieve the optimal solution.

• Input all parameters and the DMs’ preference for the problem.
• Obtain nondominated solutions by solving SP model. Set the initial intervals and decide the incumbent with DMs’ preference.
• Examine if the set of nodes is empty. If it is empty, then stop. If it is not empty, then go to Step 4.
• Divide the interval of the objective that DMs care more about by $\alpha$. Solve the SP model in left interval and compare the nondominated solution with the incumbent. If the incumbent is updated, then add the node to the set of nodes. Solve the SP model in right interval and compare the nondominated solution with the incumbent. If the incumbent is updated, then add the node to the set of nodes.
• Remove the current node, go back to Step 3.

Algorithm 1: Interactive branch-and-bound algorithm.

5. Numerical Experiments

This section will give the results of numerical experiments to verify the conclusions of this chapter on the analysis of electric vehicle tour planning problems, including the impact of different charging policies on path planning, the impact of mileage anxiety on path planning, and the impact on planning of decision makers with different expectations. We denote the EVTP model considering RA factor as EVTPRAM and the EVTP model without considering RA factor as EVTPM.

5.1. Experiment Setting

Numerical experiments were conducted on networks with various sizes and demand. For all networks, we set the RA coefficient $k=1$, the recharging rate $g_i=1,\forall i\in {\mathcal{N}}_C$. The first test networks have 10 attraction nodes and 5 recharging nodes. The second test networks have 15 attraction nodes and 6 recharging nodes. For a detailed description of our instance design, we refer the reader to Schneider et al. [21], where the generated instances are also available to be downloaded (http://evrptw.wiwi.uni-frankfurt.de). We set the value function of DMs as a linear combination of the POIs and the RA cost with weight ${\delta }_1$ and 1, respectively.

We first analyze the parameter $\alpha$ of the algorithm. There are two ways to divide the intervals when running the algorithm, e.g., half-interval ($\alpha =0.500$) and golden section ($\alpha =0.618$). We conducted experiments on 6 different networks with different values of $\alpha$. The results are shown in

Table 3. Parameter for Algorithm 1 (${\delta }_1=100$).

Network	$\mathit{\alpha }$ Value	${\mathit{z}}_1$ Value	${\mathit{z}}_2$ Value	Solving Time
c205C10	0.500	80.00	−5374.20	2.09 s
	0.618	80.00	−5374.20	3.10 s
r102C10	0.500	79.00	−3479.51	13.90 s
	0.618	79.00	−3479.51	9.43 s
r103C10	0.500	19.00	−2592.87	370.25 s
	0.618	19.00	−2592.87	348.27 s
r201C10	0.500	42.00	−3004.50	130.13 s
	0.618	42.00	−3004.50	135.82 s
rc108C10	0.500	68.00	−3618.98	40.73 s
	0.618	68.00	−3618.98	32.01 s
rc205C10	0.500	6.00	−3476.49	5.19 s
	0.618	6.00	−3476.49	6.96 s

. We can see that the value of $\alpha$ does not change the solutions.

The solving time for different $\alpha$ settings varies. In network c205C10, r201C10 and rc205C10, the solving time is shorter when $\alpha =0.500$, while in network r102C10, r103C10 and rc108C10, the solving time is shorter when $\alpha =0.618$. Nevertheless, in those networks, the value of $\alpha$ does effect the speed for solving, the impact is not significant. The average difference is only 7.27 s (the maximum is 21.98 s). Therefore, we should set the value of $\alpha$ with the information of specific networks.

5.2. Recharging Policy Analysis

During the formulation of this study, we refer to two policies for recharging, i.e., full-recharge and partial-recharge from previous literature. By comparing the two policies, we obtain the Proposition 1. In this subsection, we compare the results from solving the model with or without Constraints (5), respectively, in order to verify the effectiveness of the formulations for the two policies and the Proposition 1. The results with two different policies are shown in

Table 4. Results with two different policies (${\delta }_1=100$).

Network	Recharging Policy	${\mathit{z}}_1$ Value	${\mathit{z}}_2$ Value	Optimum
c205C10	full-recharge	-	-	-
	partial-recharge	80.00	−5374.20	2625.80
r102C10	full-recharge	79.00	−3479.51	4420.49
	partial-recharge	79.00	−3479.51	4420.49
r103C10	full-recharge	19.00	−2594.14	−694.14
	partial-recharge	19.00	−2592.87	−692.87
r201C10	full-recharge	17.00	−3027.57	−1327.57
	partial-recharge	42.00	−3004.50	1195.50
rc108C10	full-recharge	68.00	−3618.98	3181.02
	partial-recharge	68.00	−3618.98	3181.02
rc205C10	full-recharge	6.00	−3727.61	−3127.61
	partial-recharge	6.00	−3476.49	−2876.49

.

All the networks with different recharging policies have optimal solutions, except network c205C10 with the full-recharge policy, as the formulation for the full-recharge policy has more constraints, which lead to the smaller feasible region. By comparing the results from the 5 other different networks, we can first conclude the effectiveness of the Proposition 1, i.e., the partial-recharge policy is at least as desirable as the full-recharge policy. In addition, by observing the results we can see that the differences between the optimal values for different recharging policies are not always obvious, e.g., in network r102C10 and rc108C10, the optimal values are equivalent, and in network r103C10, the gap between the optimal values is only 0.18%. However, in real life, the full-recharge policy is more common and easier for EV drivers to manage. Therefore, when the difference between the optimal values with different recharging policies is subtle, the full-recharge policy might be a wiser option for DMs.

5.3. Route Reliability Analysis

We evaluate the reliability of different routes obtained by solving the model EVTPM and EVTPRAM. The results are shown in

Table 5. Results with considering RA factor or not.

Network	Preference	${\mathit{z}}_1$ Value	${\mathit{z}}_2$ Value	Objective
c104C10	RA (${\delta }_1=50$)	120.00	−3899.27	2100.73
	No RA	170.00	−8380.23	119.77
c101-21	RA (${\delta }_1=20$)	180.00	−5543.26	−1943.26
	No RA	200.00	−7488.71	−3488.71

, from which we notice that, though the POI of the routes obtained by EVTPM is slightly greater than the one obtained by EVTPRAM, the loss of POI is determined by the DMs’ preference and it assures the optimality of the objective. The loss of POI is, thus, acceptable. On the other hand, the routes obtained by EVTPRAM could save amounts of RA cost, which significantly increases the reliability of the routes. For instance, in network c101-21, The POI of the route by EVTPRAM is 10.0% less than that by EVTPM, but the RA cost declines by 26.0%.

Previous studies have shown that (EV) drivers may make a detour for refueling when gasoline (energy) reaches a low level. However, the perception of the SOC varies from person to person, e.g., some people would prefer to refuel right after the gasoline (energy) falls below 50% and some others may not think about refueling until the gasoline (energy) falls below 25% or until the refueling alert occurs [6]. We denote the referred energy as $q_0$. If the SOC of EV on link $ij$ is lower than $q_0$, then $p_{i,j}=1$ (otherwise, $p_{i,j}=0$). The range anxiety level (the probability that the EV driver might change the original route or get involved in an accident when the energy is below $q_0$) of the EV driver is defined as $0\le p_a\le 1$, then we can express the risk of the entire route as:

$$\begin{array}{cc}\hfill & P_r=1-\prod _{ij\in R}{\left(1-p_a\right)}^{p_{ij}}\hfill \end{array}$$

(27)

We compare the routes obtained by EVTPM and EVTPRAM with their risks under different $q_0$ settings. The results are shown in Figure 1 and Figure 2. (We did not consider the detour on the last link to the destination and the EV only has to save at least the required energy when it reaches the destination.)

From the figures we can see that when $q_0\le 30\%Q$, the risks of the routes planned by EVTPRAM are all 0, and when $q_0=40\%Q$, the difference between the risks of routes by EVTPM and EVTPRAM will increase as the range anxiety level goes to the middle, which is where most people’s levels stay. In general, for EV drivers who have different range anxiety levels, the risk of routes planned by EVTPM is greater than that planned by EVTPRAM, which implies that our model is effective. Moreover, by comparing the results with different values of $q_0$, when $q_0$ increases from 20%Q to 30%Q and to 40%Q, the risk of routes planned by EVTPM significantly increases, while the risk of routes planned by EVTPRAM increases much slower. Therefore, we conclude that the model we proposed is able to plan the routes with high reliability and low risk, which are more desirable for DMs.

5.4. Decision Analysis

We also study the impact of DMs’ preference on the decision of route planning when considering the RA factor. We increase the preference parameter $\delta$ from 1 to 150, and observe the change of routes to analyze the patterns for decision. The results are shown in

Table 6. Solution changes with different preferences.

Network	${\mathit{\delta }}_1=1$	${\mathit{\delta }}_1=10$	${\mathit{\delta }}_1=20$	${\mathit{\delta }}_1=30$	${\mathit{\delta }}_1=50$	${\mathit{\delta }}_1=100$	${\mathit{\delta }}_1=150$
c104C10	$R_{RA}$	$R_{BO}^1$	$R_{BO}^1$	$R_{BO}^2$	$R_{BO}^3$	$R_{BO}^4$	$R_{POI}$
c102C10	$R_{RA}$	$R_{BO}^1$	$R_{BO}^2$	$R_{POI}$	$R_{POI}$	$R_{POI}$	$R_{POI}$
c101-21	$R_{RA}$	$R_{BO}^1$	$R_{BO}^1$	$R_{BO}^1$	$R_{POI}$	$R_{POI}$	$R_{POI}$
r101-21	$R_{RA}$	$R_{RA}$	$R_{BO}^1$	$R_{BO}^2$	$R_{BO}^2$	$R_{POI}$	$R_{POI}$

and

Table 7. Objective changes with different routes.

Network	Preference	Route	${\mathit{z}}_1$ Value	${\mathit{z}}_2$ Value
c104C10	RA priority	$R_{RA}$	70.00	−2612.29
	mixed	$R_{BO}^1$	90.00	−2799.97
	mixed	$R_{BO}^2$	100.00	−3075.44
	mixed	$R_{BO}^3$	120.00	−3899.27
	mixed	$R_{BO}^4$	160.00	−6995.20
	POI priority	$R_{POI}$	170.00	−8380.23
r102C10	RA priority	$R_{RA}$	11.00	−2624.82
	mixed	$R_{BO}^1$	36.00	−2684.26
	mixed	$R_{BO}^2$	65.00	−3165.97
	POI priority	$R_{POI}$	79.00	−3479.50
c101-21	RA priority	$R_{RA}$	140.00	−5471.57
	mixed	$R_{BO}^1$	180.00	−5543.26
	POI priority	$R_{POI}$	200.00	−6973.81
r101-21	RA priority	$R_{RA}$	15.00	−2913.41
	mixed	$R_{BO}^1$	24.00	−3057.38
	mixed	$R_{BO}^2$	57.00	−3792.97
	POI priority	$R_{POI}$	65.00	−4576.02

($R_{RA}$: only considering RA factor. $R_{POI}$: only considering POI factor. $R_{BO}$: considering both factors.).

As shown in Table 6, when ${\delta }_1$ increases from a small value, the decision for routes will gradually move to $R_{POI}$ (with POI priority) from $R_{RA}$ (with RA priority), which is consonant with our expectation. By observing the value of $z_1$ and $z_2$ from each network, we notice that, in network c104C10, when $z_1$ increases from 160.00 to 170.00 (by 6.25%), $|z_2|$ increases from 6995.20 to 8380.23 (by 19.80%). Similarly, in network c101-21, $z_1$ increases by 11.11%, while $|z_2|$ increases by 25.81%, and in network r101-21, $z_1$ increases by 14.04%,while $|z_2|$ increases 20.64%. Therefore, by not considering the RA factor, DMs exchange a slight POI profit with wasting a high RA cost in many scenarios, which seems not to be a wise decision.

We do not conduct more numerical experiments on larger networks. The reasons are twofold. First, the networks that we conduct experiments on are large enough for us to verify the effectiveness of our model and algorithm, and to analyze the results for conclusions. Second, the exact algorithms for solving the NP-hard problem are not effecient enough to perform on large networks, and the algorithm for solving the subproblem is not the focus in this study, though any improvement for solving the subproblem will help shorten the solving time. Therefore, we focus on analyzing the results for advice on the management and operations for EVTP instead.

6. Conclusions

In this study, we formulated the optimal tour plan of EV as a bi-objective nonlinear programming problem. The model explicitly considers the range anxiety that the drivers suffer along the trip. We proposed a specified function to quantify the reliablity of planned trips. The model is reformulated as a bi-objective MILP and solved by an interactive branch and bound algorithm. Numerical experiments were conducted on networks with different sizes for choosing algorithm components, evaluating solution quality, and analyzing sensitivity of parameters. In particular, our algorithm is found to be effective in obtaining optimal solutions for DMs with different perspectives. In addition, we find that for EV drivers who have different range anxiety levels, the risk of routes planned by not considering the RA factor is greater than that by considering the RA factor and the model we proposed is able to plan the routes with high reliability and low risk of profit loss and accidents.

The study also has several limitations, which we plan to improve in our future work. First, while the effect of RA is modeled, EV drivers all react to the SOC of batteries in the same way. We will investigate the market with various backgrounds to understand how the RA may differ. Second, while the current solution algorithm is found to be effective, we also identified that the numerical experiments are conducted in small networks. Future efforts can be made to improve the performance of the algorithm so that we can combine our theoretical model with survey data in real-world networks to help draft practical tour plans for EVs.