The electric vehicle shortest path problem with time windows and prize collection

The Electric Vehicle Shortest Path Problem (EVSPP) aims at finding the shortest path for an electric vehicle (EV) from a given origin to a given destination. During long trips, the limited autonomy of the EV may imply several stops for recharging its battery. We consider combining such stops with visiting points of interest near charging stations (CSs). Specifically, we address a version of the EVSPP in which the charging decisions are harmonized with the driver’s preferences. The goal is to maximize the total gained score (assigned by the driver to the CSs), while respecting the time windows and the EV autonomy constraints. We define the problem as a MILP and develop an A* search heuristic to solve it. We evaluate the method by means of extensive computational experiments on realistic instances.


I. INTRODUCTION
D EVISING tools that facilitate the use of Electric Vehicles (EVs) is key to successfully electrifying transport, e.g., [5]. The limited autonomy of EVs coupled with the scarcity of charging stations (CSs) entails that long trips may involve several charging stops. Given that such stops are likely to be time consuming, we explore the idea of matching the charging stops with user preferences.
Considering an EV that needs to travel between an origin and a destination, we consider that the user attributes a score to each CS. Such scores are based on the vicinity to Points of Interest (POIs) to the CS. For example, a user may prefer to stop at cultural venues. In this case, a high score is given to CSs nearby such venues. Using those scores, we consider the problem of optimizing the shortest path respecting all energy feasibility constraints, while maximizing the total score that the user can achieve. To handle this problem we propose the Maximum Profit Model (MPM). To avoid excessively long trips, we limit the duration of the route. We establish this limit by solving the Shortest Path Model (SPM), which is the shortest EV path in time. We allow the MPM to deviated from the shortest path length within a given tolerance.
We adapt the Mixed Integer Linear programming (MILP) model proposed in [6] to handle the MPM and the SPM. which discounts the scores of the nodes into the arc costs that connect them. By doing so, we are able to efficiently adapt the A* search Algorithm to the MPM.
In general, the Electric Vehicle Shortest Path Problem (EVSPP) is a shortest path problem that accounts for battery limitation and charging constraints. Due to their limited autonomy, EVs may need to detour to CSs in order to recharge their battery. This is particularly true in medium and long range routes, like in [11]. A key decision in this context is where and how much to charge the EVs. The problem of minimizing the overall trip time for EVs in road networks was studied by [1]. A heuristic algorithm for solving large EVSPP instances was proposed in [14]. Baum et al. [2] introduced a functional representation of the optimal energy consumption between two locations, which led to developing an efficient heuristic algorithm that computes energy optimal paths.
One of the main EVSPP modeling assumption relate to how the EV batteries are recharged. Some researchers assume that the EV must completely recharge before leaving a CS, e.g. [4], [9]. Other works (e.g., [6], [10], [12]) consider the charging quantity as a decision within the optimization. As in most studies we assume that the energy consumption is directly and exclusively related to the traveled distance.
In practice, the EV charging function is nonlinear with respect to time, and depends on the used charging technology. The nonlinear charging functions were modeled as piecewise linear concave functions by [6], [10].
Time window (TW) constraints have been introduced in EV problems by [12]. TWs oblige EVs to arrive in predetermined CSs before or during a particular time interval. Contrary to what is done in the literature, we assume that TW types are given (e.g., lunch breaks), yet their location is determined from a set of CSs which accommodate this type of TW. Furthermore, we consider required and weakly mandatory TWs.
Considering that CSs have different scores and a given maximum path length (in time), the aim of the MPM is to maximize the total score of visited CSs. As such, CSs that better match the user preferences are prioritized. We consider a single EV, with partial recharging decisions and nonlinear charging functions. Furthermore, we consider only public CSs, having different technologies and scores. We also consider TWs with a limit on the total travel time. The contributions of this paper are threefold: 1) we propose an EV shortest path model that accounts for user preferences (MPM); 2) we develop a heuristic based on the A* algorithm to solve that problem; 3) and we verify the performance of both the MILP model and the heuristic algorithm on realistic test instances.
We introduce the MPM in Sec. II and develop an A* algorithm for it in Sec. III. We present our computational experiments in Sec. IV, and state our conclusions in Sec. V.

II. PROBLEM DEFINITION
The MPM maximizes the score of the CSs visited by the EV, such that the resulting path is feasible and within a tolerance from the shortest EV path. In the literature, the goal of maximizing scores obtained by visiting nodes is often called prize collection [13]. To model the MPM we adapt the arc based MILP model of [6]. Due to space limitations, we do not present the full formulation. Instead, we give an overview of that work and then detail the major adaptations made to handle the MPM.
Froger et al. [6] introduced the Fixed Route Vehicle Charging Problem (FRVCP). Given a sequence of customer nodes (i.e., not CSs) to visit, the objective of the FRVCP is to determine the charging operations (which CSs to visit and how much to charge), in order to minimize the total route duration while satisfying the following conditions: The customers in the resulting route are visited according to the given order, the resulting route is energy feasible and satisfies a maximum duration limit T max . The state of charge (SoC) of the EV is tracked on each traversed arc and visited node. The EV may be partially recharged at a set of charging stations S, which may have different technologies. For each technology we consider a nonlinear charging function, approximated via a piecewise linear function following the models by [10]. In the MPM the fixed sequence of nodes to visit goes from an origin O node to a destination D node.
Building on the FRVCP, we now describe the MPM. Let G := (S O,D , A) be a directed graph, where S O,D := S*{O, D} and A is the set of arcs that connect pairs of nodes in S O,D . Let t ij g 0 and e ij g 0 be the driving time and energy consumption of arc (i, j) * A, both satisfying the triangular inequality. The EV departs from O with a fully charged battery of capacity Q. We impose that the SoC of the EV is at least q min throughout the route. Since we consider long trip planning, we assume that the number of nights is a user input.
The user will spend some time in the neighborhood of the selected CS. Moreover, in certain moments of the day, the user may prefer to visit a CS near a POI, e.g., restaurants or hotels. CSs are typically strategically placed near those types of POIs. It is important that the user visits CSs that best suit her preferences. We attribute a score σ j for j * S, which are input to the MPM. We assume that they can be derived based on the user's ranking of POIs. Furthermore, the user may personalize her trip by imposing TWs related to an activity (e.g., lunch breaks). The classical definition of TWs determines a time to visit a given node. In the MPM, such TWs may be realized  Figure 1: Path from O to D with a night at a hotel and lunch breaks. The timestamp near a node represents the corresponding departure time, including the recharging time (not reported in the figure). For D, instead, the timestamps represent the arrival time. The number on each arc is its travel time. The red path is optimal, stopping in the hotel as required and reaching node E. If also the lunch break is required, then the EV is forced to arrive to D with stopping at L (blue path). Instead, with the weakly mandatory flag, from E the EV can go directly to D 3h in advance (green path). at a number of CSs. We therefore denote the TWs related to our MPM as Activity-based Time Windows (ATWs). In our experiments in Sec. IV, we assign to each CS a randomly generated score between 0 and 5. In practice, these values will be established by the user in real applications.
Let W be the set of all ATWs, which is partitioned into two subsets, W R , W M ¦ W, formed by the required and the weakly mandatory ATWs. The EV is forced to stop in each k R * W R , but it must stop in k M * W M only if there exists at least one k R such that k M z k R . Otherwise, k M may be skipped. We assume that the sets W R , W M are ordered by chronological order of starting time of the TWs. Each ATW is defined as: where the interval [γ L k , γ U k ] (with γ L k < γ U k ) describes its initial and ending times (the EV must arrive before γ U k ); t min k is the minimum time that the EV needs to stop during k; o k is a binary value, equal to 1 if k * W M , and 0 otherwise. The set W is ordered, thus, if k z h, then γ L k < γ L h , for any k, h * W. The EV is allowed to arrive at a node with a maximum anticipation time ϕ. In this case, the activity will nevertheless start at γ L k . Thus, if the EV arrives in node i in the interval γ L k 2 ϕ, γ U k " , then it may decide to stop at i for at least t min k or instead perform k in a subsequent node. In addition, we use hard TWs, which means that it is not possible to perform the ATW k before γ L k 2 ϕ nor after γ U k . Finally, two ATWs can overlap, but they cannot be contained in one another. Each CSs is associated with multiple POIs, and each TW requires a specific POI. This information is stored in the label ν k and is different for each TW. For instance, if k * W refers to the first night, then ν k = "Hotels" and the EV is forced to stop at a CSs near a hotel. So, it is possible to construct the set of chargers S k ¦ S that have in their neighborhood the POI stated in ν k .
To determine T max we solve an EV Shortest Path Problem SPM, subject to the same constraints as MPM, but with the goal of minimizing the total trip duration (i.e., ignoring the scores). Let T opt be the optimal objective function value of SPM, which is a theoretical lower bound on T max . Then, let T add be the tolerance of the total additional detouring time from the shortest path to stop in nodes with higher scores. Therefore, we set T max to T opt + T add .
To avoid routes with many stops at SCs, we introduce the parameter r min , as r min := +0.4 (Q 2 q min )/η,, where η is the average energy consumption per kilometer (expressed in kWh/km), and is dependent on the EV type. The ratio (Q 2 q min )/η represents the maximum autonomy of the vehicle, excluding the minimal amount of energy that is always required. With this arrangement it is possible to prune all the arcs associated with a distance less than r min , which also avoids stops for charging the EV in consecutive locations that are very close to each other.
Let ξ denote a lower bound on the distance of a trip from O to D (its computation is described in Sec. IV). Then the maximum number of legs N in a path is defined as +1.5+ξ/r min ,,.

III. A* SEARCH ALGORITHM
In the following we propose a heuristic algorithm for the SPM and extend it for solving the MPM. Recall that the SPM solution serves as an input to MDPM. The algorithm is based on the A* search, which finds a path from an origin O to a destination D with the smallest cost. To do that, it maintains the tree of all the paths originated from O and extends each of them one arc at a time until D is reached. It uses a best-first search by selecting the node that minimizes where i is the current node, g(i) is the cost of the path from O to i, and h(i) is an estimate of the cost of the shortest path from i to D. If h(i) never overestimates the real cost to reach D from i, for all i, then the A* algorithm finds the optimal solution.
The A* algorithm is based on minimizing cost mechanisms, thus it is not directly applicable to prize collection settings such as the MDPM. Therefore, we propose the MDPM which approximates the MPM, searching for a shortest path while maximizing the total score. To do that, we assign a weights ij with each arc, defined as where ∆ j is the the time spent waiting while the EV is charging at the CS j, and µ is a parameter that indicates the relative importance of the score with respect to the time required to go from i to j, charging time included. The objective function of the MDPM model is then This expression is nonlinear, but it can be linearized by introducing new decision variables s ij , and changing the objective function as We first discuss the computation of the potentials used to estimate the heuristic function h. We then incorporate the TWs in the heuristic. Finally, we modify the algorithm to account for the scores in h.
Let q i and q i be the SoC when the EV arrives and departs from CS i. The variables c i and c i are respectively the start and end time for charging an EV. Variables τ i , τ i track the time when the EV arrives and leaves CS i * S, respectively. There is also a tolerance ϕ i that represents how much time in advance, with respect to γ L k , the EV can arrive in i, for any i * S D . The maximum anticipation time is set to ϕ, but even if the EV arrives in advance, the minimum stopping time t min k starts at γ L k and not before. Let x ij be a binary variable which equals 1 if the EV arrives at node j from node i, 0 otherwise. The variable y jk is also binary and it is 1 if the EV stops in j in TW k, 0 otherwise. The maximum duration of the trip is T max (in Sec. III we show how to compute an upper bound for this value). The variable z k is binary and is equal to 1 if the EV arrives in D after TW k, 0 otherwise, for any k * W M . It is used to link the arrival time in node D and avoidable TWs.

A. Potentials
We now find an initial estimate of the total time from any CS i to D, building on techniques used in [14]. We start by dropping the charging and the TW constraints, thus we obtain a problem that can be solved using the Dijkstra's algorithm. Let G := ïS O,D , Að be the directed graph from O to D, where S O,D is the set of nodes and A := S O ×S D the set of arcs. We then apply the backward Dijkstra's algorithm, which operates on the reversed graph while considering D as the origin and O as the destination. This allows us to obtain a lower bound on the driving time from any i * S O,D to D, which we denote by π dr (i).
To account for energy consumption, we apply the backward Dijkstra's algorithm considering the energy consumption as the weight for the arcs. This allows us to derive the minimum amount of energy required from i * S O,D to D. We denote this lower bound by π cons (i).
The EV arrives partially charged at each node, with an amount of energy equal to SoC(i). Since the minimal amount of SoC q min needs to be respected at every node, we can think of SoC(i) 2 q min as the available energy at node i * S O,D . Then, to compute the minimal amount of energy from i to D, we defineπ cons (i) := π cons (i) 2 (SoC(i) 2 q min ). We now compute a lower bound for the charging time, based on the minimal required energy at a node i * S O,D . We define G i := ïT i , A i ð as a subgraph of G, where T i is the set of reachable nodes from i, and A i is the set of arcs comprising a path from i to any j * T i (see Fig. 2). Let s max (i) denote the maximum charging rate of all the CSs in T i . We denote the maximum charging rate of CS j asρ j , which corresponds to the largest slope of the piecewise charging function of j. We then set s max (i) as max {ρ j : "j * T i }. This improves the computation of the charging potential with respect to [14], where s max is constant and does not depend on the nodes that are reachable from i.
If the energy available at i is greater than the remaining energy needed to reach D, thenπ cons can be negative. Thus, two cases are considered when computing a lower bound for the charging time: This gives a potential that returns the minimal charging time from any node i to D. Thus, a lower bound on the total trip time is given by We improve this lower bound by accounting for TWs. Let π tw (i) be the minimal stopping time that the EV must perform from i to D according to the ATWs. The trip accounts for the sum of all the minimum stopping times. Letk be the last TW in the ordered set W R . Suppose that, when at node i, the EV has not performed all TWs in W R . Then, g(i) f γ L k + t miñ k , where g(i) is the arrival time at i and γ L k is the starting time of TWk. In this case, we compute a lower bound for the TWs potential as the sum of all the stopping times that have not been performed by the time g(i). If instead, in node i, the EV has already done all the TWs in W R , then g(i) > γ L k + t miñ k and so the potential for the TWs must be zero. Therefore, if π tw (i) is the TWs potential for node i, we have For instance, suppose that W R contains a 1 hour stop for lunch, one tourism stop for 2 hours and one for sleeping for 11 hours. Then, before lunch π tw (i) = 14, after lunch π tw (i) = 13, and the next day π tw (i) = 0. We now incorporate π tw in π tt . The EV charges during each stop. Thus, we cannot simply sum π tw and π ch , since they may overlap. Instead, we can obtain a better lower bound considering the maximum between π ch and π tw , and then adding the driving potential π dr . So, defines the lower bound for the total stopping time from i to D. We are not overestimating the real cost, since taking the maximum we consider for each node the best possible charging scenario that at least the EV has to perform.

B. Labels
Label L jm represents the state of the EV when arriving at the m-th copy of node j, that is dynamically allocated. With this label we keep track of the state of the EV. For notational convenience, we denote g(j m ) by g m j . Each label includes the total time needed to reach j m , so it includes both driving and charging times. For instance, suppose that the EV goes from the n-th copy of i to the m-th copy of j, namely from i n to j m . Then the label L jm considers driving from i to j and the charge in i that is needed to reach j for the m-th time. It excludes the time the EV spends charging in j. The label of the m-th copy of node j, with i n as its direct predecessor, is: where " i is the node from which the EV arrives to j m ; " g m j is the total travel time from O to j m ; " h m j is the estimated travel remaining time from j m to D; " f m j is the estimated arrival time at D if the path from O to j m is performed. It is given by f m j := g m j + h m j ; " p m j is the label from which the EV arrives, i.e., L in , which is the label of the n-th copy of node i, with n * { 1, . . . , M i }; " β m j is the additional time spent at j for charging, it is chosen from an ordered set β := { β 1 , β 2 , . . . , β s }; " q m j is the amount of energy that is charged in the predecessor node i n . It is computed to at least respect the consumption e ij and is given by q m j := q m j 2 q m j , where q m j := max q n i + e ij , Q ; " q m j is the amount of energy of the EV when arriving at j m . It is computed as q m j := q n i + q m j 2 e ij ; " λ m j is the minimum time the EV must charge at j m . This amount of time is considered in the next label and not in L jm . It is defined as where γ L k is the starting time of TW k and t min k is its minimum stopping time. TW k is retrieved using ω m j (see below); " ∆ m j is the charging time in i n . It is given by the inverse of the charging function in node j. The term β m j is added to consider the cases in which the EV charges more than needed; " ω m j is the index of the last TW k prior to node j m . Parameters ∆, q and q are strictly related to β. As a consequence, also g, h and f depend on β. The value of λ instead depends strictly on ω.

C. A* search algorithms for the MPM and MDPM
Using the labels described above we now outline the A* search algorithm. We start by implementing a heuristic approach to find the fastest path from O to D, referring to this as AsM.
The origin node O is initialized as follows: We now introduce some of the functions used later in the pseudocode. The function POP(Q) returns the label with the lowest f in Q, while function PUSH(Q, L jm ) adds the label L jm to the queue. Function STAR(G, j) returns the set of nodes h * S D such that (j, h) * A, in descending order with respect to t jh . The function NEXTTW(W R , ω m j ) returns the next TW in W R that is not visited when the EV arrives at node j m . The EV is allowed to arrive at a node with a maximum anticipation time of ϕ. The boolean function NODEHASPOI(i, ν) returns one if there is at least one POI of the category ν in the neighborhood of node i, and zero otherwise.
Function MAXSLOPE(Ṡ) applied to a generic subsetṠ of CS S, returns the maximum slope between all the charging functions of nodes inṠ. To speed up the algorithm, all the subtrees are precomputed. The function ROUTING(i, j) returns the pair (t ij , e ij ), that are respectively the time and the energy required to go from i to j, for any i * S O , j * S D . The function MINSTOP(g i ) returns π tw (i).
The A* algorithm is described in Alg. 1. It starts by initializing the counters for all the copies and storing the subtree of each node. Then the label associated with the origin is created and added to the queue and to the set of all the labels. while Q do 7:

9:
if i = D then 10: if ω n i < |W R | then go to 6 11:

46: end function
A Check if the subtree of current node contains the category of POI that is needed for the next TW and the subsequents ones. Otherwise, goes to the next element in the queue.
The label L n i with the lowest value of f is selected and then the algorithm finds the next TW k that must be performed. A check is performed to verify if the current label entails the arrival to the destination point: if so, we verify whether the index ω n i of the last visited TW is at least equal |W R |. If this is not the case, the EV has not visited yet all the nonavoidable TWs, and thus the current label must be discarded. Otherwise, the current label and the set of all the generated labels are returned and the search terminates.
At this point, if k is not NONE the algorithm checks whether or not there exists at least one node in T i of the current node i in which the POI constraint of TW k is satisfied. Then it checks the same for all the TWs in W R that must be performed after k.
The algorithm proceeds by considering the nodes j such that (i, j) * A. A sequence of operations assures that the trip from i n to j m is feasible, where m is the index of the m-th copy of j that will be created if all the checks are passed. First the energy constraints. If the energy required on arc (i, j) is greater than the maximum amount for the EV, we discard this label, and the algorithm proceeds to the next node in STAR(G, i). Otherwise, if e ij < Q 2 q min , we consider charging a greater amount of energy with respect to e ij by looping on β.
At line 25 of Alg. 1, the charging time and the charged energy are computed, given the current SoC q n i and the amount of energy required e ij . The two values are then updated on line 26, taking the current SoC and the arrival time.
We then compute a temporary value g temp of the arrival time in j. In case g temp > t end , the current label is discarded. We then verify if the selected TW k is in W R . If so, then the algorithm must satisfy the constraints associated with k. First, if g temp > γ U k , then the arrival time at j will be after the ending time of k, which is not feasible. If instead, g temp is included in the range [γ L k 2 ϕ, γ U k ], then the EV arrives at j during the TW k. In this situation, the user may decide to stop in j, and respect the TW constraints, or to continue driving to the next node h and see if it is possible to respect k in node h. This scenario models the case in which, for instance, instead of respecting the lunch constraints in node j, the user prefers to drive longer and eat at another place. In the case we stay in j to respect TW k, we check whether node j has the POI that is required for k. If so, a label L m j is created, imposing that ω m j := ω n i + 1 (from j m on, TW k is respected) and λ m j = max 0, γ L k 2 g temp + t min k . The max function is used to compute how much time in advance the EV arrives in j, so the minimum stopping time imposed from any arc from j m is λ m j . If instead node j does not have any POI of the given category ν k , we can step over and create a label that goes from i n to j m without imposing a minimum stopping time λ m j . In this case λ m j := 0 and ω m j := ω n i . In both cases, the newly created label L jm is added to the set of labels L and pushed to the queue Q. Finally, the loop on β continues after updating β. Alg.s 2 and 3 outline the creation of a new label, and its heuristic. Finally, if it is not possible to reach D respecting all the imposed constraints, the algorithm returns NONE .
Algorithm 2 CREATELABEL function. It creates the label from node a to node b, given the arrival time g in a, the SoC q, the charging time ∆, the charged energy q, the driving time and energy, t, e, the minimum amount of time needed to charge in node b in the next label, the index of the last performed TW ω, the charger additional time β and the previous label L. 1: function CREATELABEL(a, b, g, q, ∆, q, t, e, λ, ω, β, L) 2: g := g + ∆ + t; q := q + q − e 3: h := HEURISTIC(b, q, g) 4: f := g + f 5: L := (a, g, h, f, L, β, q, q, λ, ∆, ω) 6: return L 7: end function We modify the previously described A* search algorithm for the MDPM. We refer to this algorithm as AsDM. We add Algorithm 3 HEURISTIC function. It returns the estimated remaining time from node i to D in graph G, considering the SoC q and arrival time g. 1: function HEURISTIC(i, q, g) 2: πtw(i) := MINSTOP(g) 3: if q − q min ≥ πcons(i) then 4: π ch (i) := 0 5: else 6: smax := MAXSLOPE(SUBGRAPH(G, i)) 7: π ch (i) := [πcons(i) − (q − q min )]/smax 8: end if 9: return π dr (i) + max { πtw(i), π ch (i) } 10: end function two parameters in the labels, as follows . The first parameter is Σ m j , which represents the total score gained from O to the current label. For copy j m , The second parameter is τ m j , it represents the arrival time in copy j m . All the TWs constraints are now satisfied using this parameter instead of the arrival time g. This means, for instance, that τ temp = τ m j + ∆ 2 t ij , and every g temp is replaced with τ temp . The minimum waiting time becomes λ m j := max 0, γ L k 2 τ m j + t min k . For the origin node we have τ 1 O := t start and g 1 O := 0. Finally, due to the different definition of labels, the function CREATELABEL is replaced by CREATELABELDISCOUNTED, and is outlined in Alg. 4.

A. Data description and preprocessing
We obtain CS locations and specifications from a company, whose name we cannot disclosed due to privacy reasons. The CSs were extracted from a bounding box that goes from 43.55°t o 49.05°latitude and from 8.68°to 13.11°longitude, covering parts of Italy, Germany, Austria, Switzerland and the totality of Liechtenstein and San Marino (see Fig. 3). CS i has a charging power p i . We categorize the CSs as follows: " slow: p i f 11 kW, " medium: 11 kW < p i f 30 kW, " fast: p i > 30 kW. We have performed a number of preprocessing actions on the provided data set. The first step was to group CSs at similar locations, for each group the fastest charger within that group was considered. We then considered that CSs within a radius of r M = 100 m form a cluster and if one CS belongs to two or more clusters, then all of them are merged. This reduction is based on the Haversine formula for computing the distance between each pair of nodes and assumes that the distances are symmetric. Merging all the CSs into a cluster produces a new CS that replaces the other CSs in the cluster. The cluster position is set as the average of the GPS coordinates of the CSs that comprise it, while its charging speed is set as the maximum among the CSs that comprise it. The final dataset is denoted by Γ. The final number of CSs is summarized in Tab. I. Using the post-processed database, we computed the distance and time matrix for each pair of CSs, obtained from multiple requests to the Open Source Routing Machine API (OSRM, [8]) server. The energy required by each arc was instead given by the product of the arc length and the average consumption per kilometer, as in [10].
The sets of required W R , and weakly mandatory time windows W M are given as input. Furthermore, we consider that tourism AWTs W T are also given as input. For k * W T a single CS location P k is given. To construct the set S k , the Haversine distance is computed by searching for CSs within a radius r T centered in P k . The CSs with a distance less than r T are included in S k . If no eligible nodes are present, the radius is iteratively increased by δ T and the search is repeated, up to a maximum of r T . If S k is still empty, the computation is resumed and an error informs that it is not possible to reach P k unless the EV is left more than r T away. This process simulates the approach that a user needs to search for a CS near the location she wants to visit. If no CS is available at a reasonable distance from P k , then the user can decide to remove or change that tourism stop. For our experiments, we used r T = 2 km, δ T = 200 m and r T = 4 km.
Given the origin O, the destination D and the set of tourism stops W T , we construct the set of chargers S and the set of arcs A. Using OSRM, the optimal path without charging stops is computed and is used as a lower bound for SPM. All the CSs that are within a range of 5 km centered on the OSRM optimal route form the set S. Finally, we construct A as the set of all arcs (i, j) connecting nodes in S O,D .
In addition, we remove from S all the nodes that have a distance from O that is less then r min , due to the fact that we want a solution with a small number of stops. Selecting one of those would have caused the EV to stop for just a few kilometers from the origin point, which is not desirable. The same reasoning is applied to D.
We also remove arcs that are going in the opposite direction with respect D. More precisely, for each arc (i, j) * A, if by going from i to j the distance to D is reduced, then (i, j) is kept, otherwise it is deleted. This process nearly halves the amount of edges included in A.
Since the EV needs to respect the battery capacity Q and we want the total travel time to be small, we consider only the arcs (i, j) such that r min f d ij f r max , where r max := +(Q 2 q min )/η,.
Finally, nodes that were not reachable from any other CS, as well as the corresponding arcs, are deleted.

B. Experiments and results
The codes were implemented in Python, using IBM ILOG CPLEX Optimization Studio 20.1. The tests ran on a single core of an Apple MacBook Pro with 8 core Apple M1 processor of 3.2 GHz, with 8 GB of LPDDR4 RAM.
The considered EV was a Škoda Enyaq iV 60, with a net battery capacity of Q := 58 kWh and a maximum average consumption of η := 0.187 kWh/km [3]. It has a maximum power charge P := 40 kW, and we set the minimum required energy level q min to 15 kWh.
Two subsets of CSs were created: a small one, Γ 1 , with 650 CSs, used to compare the exact solution obtained with the MILP models with the A* search algorithms, and a larger one, Γ 2 , with 5 813 CSs, used to test the A* search. Both datasets are extracted with uniform probability from Γ, so we must guarantee that in each tourism stop there is at least one CS. We created a 5 km×5 km square centered in the coordinates of each tourism stop and uniformly extracted four CSs. Γ 1 and Γ 2 were obtained by selecting uniformly a certain amount of CSs from Γ and adding them to the ones selected for the tourism stops.
For both sets of CSs we use the same set of instances defined as trips. Since our CSs lay in a rectangular area that covers part of central Europe, those instances are chosen so that they fully lay in the same geographical area. We generate three main trips, with some variations, such as starting time, presence of tourism stops and minimum stopping times, presence or not of lunch (1 hour) and night (11 hours " from Genoa to Zürich denoted by GeZu; " from Livorno to Regensburg denoted by LiRe; " Stuttgart to Ancona denoted by StAn The remaining details are summarized in Tab. II. The name of the O 2 D is followed by a running number to distinguish the instances. Overall, we created 20 instances. We tested the MILP shortest path model, SPM, and maximum profit model, MPM, as well as both A* heuristics, the A* search algorithm for the SPM, AsM, and its variant for the maximum discount profit model, AsDM, for set Γ 1 . In addition, only the heuristics AsM and AsDM were applied to Then for all AsDM models we solved SPM imposing that the EV must use the arcs selected by the heuristic. CPLEX was unable to solve any instance in Γ 2 within one hour. We denote by GS y x (GT y x ) the relative gap between the total score (trip time) of the models x and y. Also, to ease the notation, in the following we denote by AsDM µ0 the application of AsDM for µ = µ 0 , and use TS for the total score, TT for the trip time, and RT for the run time.
As seen in Tab. II, some instances describe the same trip, with different timings, and thus lead to the same set of unconstrained shortest optimal paths, and to the same set of CSs. For this reason, we subdivide the instances, and for both Γ 1 and Γ 2 extract the set of CSs related to each subdivision and store them. In this way, instances in the same subdivision use the same graph for all the models. In Tab. III we list the size for each subdivision. We use the dataset Γ 1 to analyze the performance of the A* search with respect to the exact solution found. Initially we find the shortest path value T opt of SPM and AsM. Then we compute the maximum relaxed time T max = T opt + T add and solve the MPM model. For the evaluation, we fix the value of T add to 1h 30'. The result is then compared to the score gained  Tab. IV and V. First, the solutions obtained with the A* search are refined using a MILP formulation. In particular, we create another SPM in which we set x ij = 1 for each arc (i, j) selected in the final solution of AsDM for µ = 2. The total trip times obtained by both approaches are quite similar. The results are summarized in Tab. IV and compared with the shortest path model SPM. The total trip time increases in average 3.3% when the refined SPM is applied, which is a small detour with respect to shortest path. In Tab. V we compare the total score gained and the run time of MPM and AsDM for µ = 2. In average, the relative change between the maximum score computed with MPM and the one computed with AsDM there is an average of 11.9% worse scores with respect to the optimal solution found by MPM. This last value is quite large, with some instances having a relative change of over 25%. A possible approach to obtain better results may be to rely on different values of the parameter µ. The run time is quite low for all the models, but it is worth noting that we are considering the small graph created with the CSs in Γ 1 .
b) Medium dataset: We want to test the A* search algorithm for the medium dataset Γ 2 . We tested all the instances using only the heuristic approaches, AsM and AsDM. The results are reported in Tab. VI and VII. According to Tab. VI, we obtain an average increase of the total score with µ = 1 with respect to µ = 2. The difference is due to the fact that with µ = 2 the shortest arcs become more important, even StAn8, the instance with more TWs, is solved in less than 1 sec. in each discounted model. When the TWs are balanced along the trip, with not too many uncovered time intervals, the computation with the A* search is very fast, even in medium sized graphs. In Tab. VII we can see an average total trip time variation of less than 6.0% in the AsDM models with respect to the AsM.

V. CONCLUSIONS
We investigate the problem of finding a path for an EV performing a a long trip. The objective is to select CSs that better match the user preferences while respecting ATWs constraints.
We proposed an A* algorithm AsM for the shortest path version of our problem. We then proposed the AsDM algorithm which privileges POIs preferred by the user. This algorithm was quite fast in finding a shortest path solution with high total score. The algorithm performed fairly well with respect