Modeling and Engineering Constrained Shortest Path Algorithms for Battery Electric Vehicles

We study the problem of computing constrained shortest paths for battery electric vehicles. Since battery capacities are limited, fastest routes are often infeasible. Instead, users are interested in fast routes on which the energy consumption does not exceed the battery capacity. For that, drivers can deliberately reduce speed to save energy. Hence, route planning should provide both path and speed recommendations. To tackle the resulting NP-hard optimization problem, previous work trades correctness or accuracy of the underlying model for practical running times. We present a novel framework to compute optimal constrained shortest paths (without charging stops) for electric vehicles that uses more realistic physical models, while taking speed adaptation into account. Careful algorithm engineering makes the approach practical even on large, realistic road networks: We compute optimal solutions in less than a second for typical battery capacities, matching the performance of previous inexact methods. For even faster query times, the approach can easily be extended with heuristics that provide high quality solutions within milliseconds.


Introduction
Battery electric vehicles (EVs) have matured, giving the prospect of high powertrain efficiency and independence of fossil fuels, but a major hindrance of their adoption remains the limited battery capacity of most vehicles combined with a lengthy recharge time. To overcome range anxiety, careful route planning that prevents battery depletion during a ride is paramount. To avoid time-consuming charging stops, it can pay off to save energy by deliberately driving below posted speed limits, especially along high-speed roads. Speed planning becomes even more relevant with the advent of autonomous vehicles, where driving speeds can be planned in advance to ensure that the target is reached. Another substantial difference to vehicles run by combustion engines is the ability to recuperate energy when braking. In this work, we study route planning algorithms that incorporate these modeling considerations to capture the characteristics and needs of EVs. We discuss different ways to model adaptive speeds and propose algorithmic solutions to solve the resulting N P-hard problem of finding a fast route on which the target can be reached without intermediate recharging.
We demonstrate the practicality of our approaches in a comprehensive evaluation on realistic input.
called Shortest Path Problem with Resource Constraints (SPPRC) [55], which exists in many variants and can be seen as a generalization of CSP in which a shortest path subject to one or more resource constraints is sought. These constraints can be represented by resource extension functions [54], which define the change in resource consumption or state when traversing an arc of the network [78]. Often, these state transitions are modeled as linear [53,63,85] or nonlinear functions [52]. Another common approach is to use parallel arcs to represent different tradeoff choices [39,62]. For an overview of variants of SPPRC and the VRP with realistic constraints and road attributes, see Ticha et al. [83]. Works dealing with speed planning can be found in Bektaş et al. [18] and Fukasawa et al. [37]. Another related problem is the Pollution Routing Problem [19], which aims at minimizing In each setting discussed below, we are given a directed graph G = (V, A). Arcs in the graph are directed to model, e. g., one-way streets. Traversing an arc in the graph consumes both time and energy. While the time to traverse an arc is always (strictly) positive, energy consumption can become negative to model recuperation, e. g., when going downhill. Moreover, the vehicle is equipped with a battery that has a given capacity M ∈ R ≥0 . This imposes constraints on the feasibility of paths, as the state of charge (SoC) of the vehicle must not drop below 0. Additionally, the SoC can become at most M after traversing an arc with negative energy consumption.
In what follows, we focus on computing fastest feasible paths for EVs, i. e., fastest paths under the constraint that they are reachable by the EV (given its initial SoC). However, we point out that it is not hard to adapt the algorithms we present to closely related problem settings, such as • computing a path with minimum driving time such that the target is reached with at least a certain minimum SoC; • computing a path with minimum driving time such that the SoC does not fall below a certain threshold at any point during a journey; • computing an energy-optimal path that does not exceed a certain travel time; • computing the full Pareto set of nondominated solutions at the target (wrt. driving time and energy consumption).

Constrained Shortest Paths for EVs
In a problem resembling traditional CSP settings, we are given a graph G = (V, A) together with two cost functions d : A → R >0 and c : A → R representing driving time and energy consumption on an arc, respectively. Although energy consumption can become negative, physical constraints prevent that cycles in the graph have negative energy consumption. In addition to that, the battery constraints mentioned above apply: If the battery of the EV has the SoC b u ∈ [0, M ] at some vertex u ∈ V , then traversing an arc (u, v) ∈ A results in the SoC Here, an SoC of −∞ implies that the arc cannot be traversed because the battery would run empty. A path is feasible if the SoC never drops to −∞ when traversing its arcs. We consider the EV Constrained Shortest Path (EVCSP) Problem, which is defined as follows.
Definition 1 (EV Constrained Shortest Path). Given a source vertex s ∈ V , a target vertex t ∈ V , and an initial SoC b s ∈ [0, M ], find a path that is feasible and minimizes driving time.
Note that if we did not allow recuperation and hence, all consumption values became nonnegative, this problem would be equivalent to the well-known N P-hard CSP problem [47], which (in the terminology of our setting) asks for a path with minimum driving time such that energy consumption does not exceed the threshold b s . Consequently, our problem at hand is N P-hard, too.
The Bicriteria Shortest Path Algorithm. To solve the EVCSP problem, we can adapt the wellknown exponential-time bicriteria shortest path (BSP) algorithm [48,68] in a straightforward manner. The algorithm maintains label sets, in our case containing tuples (τ, β) of driving time τ ∈ R ≥0  (20,6) at v and a label (21,4) at t, which dominates the label (21,2) in the current set. (b) The initial SoC is 5. Assume that the original graph consists only of the vertices s and t, connected by the two dashed arcs. Thus, there are two labels (2,1) and (4,3) at the target vertex when the algorithm terminates.
Inserting the vertex v, e. g., for visualization in practice, results in the depicted graph with solid arcs. Even though the distances between s and t are maintained, the algorithm now computes the additional label (3,2) on the modified input.
and SoC β ∈ R. A label (τ, β) dominates another label (τ , β ) if τ ≤ τ and β ≥ β . Starting with the label (0, b s ) at the source s, the algorithm settles in each step the label = (τ, β) with minimum driving time τ ≥ 0 among any labels not settled before. For all arcs (u, v) ∈ A outgoing from the vertex u ∈ V this label belongs to, it then generates a new label with driving time τ + d (u, v).
We set the SoC of the new label to min{M, β − c(u, v)}. If this results in nonnegative SoC, the label represents a feasible path and is added to the labels at v if it is not dominated by any of them. Labels at v dominated by are discarded. Given that driving times are strictly positive, this algorithm is label setting, i. e., settled labels are never dominated afterwards. 1 An optimal (constrained) solution is found when a label at the target t is settled.
Modeling Adaptive Speeds in the EVCSP Problem. The BSP algorithm finds the fastest route subject to battery constraints. As argued before, travel time and energy consumption are not only affected by the choice of the route itself, but also by driving behavior. Allowing multiple driving speeds (and consumption values) per road segment, one could, e. g., save energy on the motorway by driving at reasonable speeds below the posted speed limits. Therefore, we consider adaptive speeds, i. e., we allow the EV to adjust its speed within reasonable limits to reach the destination as fast as possible and with sufficient SoC. In addition to the actual route from the source to the target, we then also have to specify (optimal) driving speeds along that route.
A straightforward way to model adaptive speeds is to sample reasonable speeds for each road segment [12,45,50,81]. Then, one can add parallel arcs in the underlying graph representation, which correspond to alternative driving speeds (inducing certain values of driving time and energy consumption); see Figure 1. The major benefit of this approach is its simplicity: Since we still obtain an instance of EVCSP, we can immediately apply the BSP algorithm to solve the problem. However, it also comes with several drawbacks. First of all, parallel arcs greatly increase running time, due to a larger number of nondominated solutions. In fact, the number of nondominated (i. e., Pareto-optimal) solutions can be exponential even on a single route; see Figure 1a and c. f. Letchford et al. [62]. Consequently, only heuristic algorithms achieve practical running times [12,45,50]. By discretizing a continuous range of tradeoffs, parallel arcs that model alternative speeds have other undesirable effects, such as producing many insignificant, yet nondominated solutions. Figure 1a shows such an example in which some solutions at the target vertex t provide rather unattractive tradeoffs, namely, spending ten extra units of time to save only one unit of energy. Adding another sample (indicated by the dashed arc), on the other hand, would result in a new label at t that dominates one of these less favorable solutions. In other words, the number of samples influences both running time and quality of the results: More samples increase running time, but fewer samples reduce quality. Similarly, adding or contracting degree-two vertices in the graph, which are commonly included for visualization purposes, affects the solution space even if distances are maintained; see Figure 1b. This is clearly not desirable, since such modeling decisions should not have any impact on the optimal solution. To remedy the above issues, we propose a more sophisticated model, which uses continuous functions to model the tradeoffs on arcs.

Constrained Shortest Paths with Continuous Adaptive Speeds
In the EVCSP problem, scalar values represent driving time and energy consumption of an arc a ∈ A in the input graph G = (V, A). Instead, we now assume that there is a tradeoff function g a : R >0 → R based on a physical consumption model, mapping a desired driving time x ∈ R >0 along a to the resulting energy consumption g a (x). In reality, the driving time on a road segment cannot be chosen arbitrarily. Lower bounds are induced by speed limits and the maximum speed of a vehicle. On the other hand, driving slower than a reasonable minimum speed would mean to become an obstacle for other drivers. Additionally, there is a certain point at which driving slower will no longer pay off in terms of energy consumption. This yields (positive) minimum and maximum driving timesτ a ∈ R >0 andτ a ∈ R >0 , respectively, for the function g a , withτ a ≤τ a . We incorporate these bounds into a consumption function c a : g a (τ a ) if x >τ a , g a (x) otherwise. (1) Thus, driving times belowτ a are infeasible (modeled as infinite consumption) and driving times aboveτ a become unprofitable. A driving time x ∈ R ≥0 is also called admissible if x ∈ [τ a ,τ a ], i. e., it lies in the relevant subdomain of c a . In the special (degenerate) caseτ a =τ a , the function c a represents a constant pair (τ a , c a (τ a )) of fixed driving time and energy consumption. We call a and c a constant in this case, as the arc allows no speed adaptation. Due to physical constraints, the minimum values c a (τ a ) of the consumption functions c a of arcs a ∈ A must not induce cycles with negative energy consumption in the graph (going in a cycle cannot increase the SoC). As before, we assume that the EV is equipped with a battery that has a certain capacity M ∈ R ≥0 and that its SoC must not drop below 0 nor exceed M . When incorporating these constraints into our setting, we obtain a bivariate SoC function f a : R ≥0 × [0, M ] ∪ {−∞} → [0, M ] ∪ {−∞} for every arc a = (u, v) ∈ A, mapping the SoC at u to the resulting SoC at v when traversing a with a specific driving time. The function f a is given by where an SoC of −∞ denotes an empty battery. Hence, f a (x, b) = −∞ implies that the arc cannot be traversed at the corresponding speed (as it would cause the battery to run empty). Further, we have f a (x, −∞) = −∞ for arbitrary x ∈ R ≥0 . Given the SoC b s ∈ [0, M ] at a source s ∈ V , an s-t path [s = v 1 , . . . , v k = t], and driving times x i ∈ R ≥0 for the arcs (v i , v i+1 ) of the path, with i ∈ {1, . . . , k − 1}, we can determine a corresponding SoC b t ∈ [0, M ] ∪ {−∞} at the target t ∈ V as follows. Starting at the source s = v 1 , the SoC at v k = t is obtained after iteratively evaluating the SoC function f (v i ,v i+1 ) at x i and the SoC at the previous vertex v i . Formally, we get Note that b t = −∞ holds if the path is infeasible (for the given driving times).
Using this more sophisticated model, Definition 2 formally introduces the EV Continuous Adaptive Speeds Shortest Path (EVCAS) Problem, which we examine in the subsequent sections of this work. Equation 3, and the overall driving time

Definition 2 (EV Continuous Adaptive Speeds Shortest Path). Given a source vertex s
An instance of EVCAS where all functions are degenerate constant tuples is also an input instance to EVCSP, which is N P-hard. Hence, EVCAS is N P-hard as well.

Solving EVCAS on a Simplified Example Instance
In this work, we focus on solving the EVCAS problem. To gain insights about the structure of optimal solutions, we now derive consumption functions and (bivariate) SoC functions for given paths instead of arcs. We illustrate such functions in an example using piecewise linear tradeoff functions. Note that piecewise linear functions are commonly used to model constraints in variants of the SPPRC [53,63,85] and in time-dependent route planning [7,14,26,60]. In Section 3, we propose a more realistic nonlinear model, which is used in the remainder of this work.
For now, assume that the tradeoff function g a of every arc a ∈ A is decreasing and linear, i. e., g a (x) = αx + β for all x ∈ R ≥0 , where α ∈ R ≤0 and β ∈ R are constant coefficients. The values α and β may differ between arcs, though, to reflect different road types or other relevant factors [22,90]. Figure 2a and Figure 2b show consumption functions (plugging in limitsτ (u,v) To find the minimum driving time from u to w, we are interested in the consumption function of the path P = [u, v, w], i. e., a function c P that maps driving time x ∈ R ≥0 spent on the u-w path to the minimum energy consumption c P (x) on the path. Formally, to get the value of c P (x) for some driving time x ∈ R ≥0 , we have to pick (nonnegative) values x 1 ∈ R ≥0 and x 2 ∈ R ≥0 , such that Figure 2c indicates possible distributions of driving times x 1 and x 2 among the two arcs and the resulting energy consumption. The lower envelope of this area yields the desired function c P . Intuitively, we want to spend as much of the available extra time (exceeding the minimum 2) as possible on the arc that provides the best tradeoff, i. e., the consumption function with the steeper slope (where spending additional time saves most energy). As a result, the consumption function of a path is always convex on the subdomain where its image is finite; see Figure 2c. Observe that, while tradeoff functions of single arcs a ∈ A are linear on the interval [τ a ,τ a ] of admissible driving times, the tradeoff function of a path is piecewise linear on the interval induced by its minimum and maximum driving time. The number of linear subfunctions defining the function c P is bounded by the number of arcs in the path.
Battery Constraints. The situation becomes more involved if we also take battery constraints into account. Then, energy consumption not only depends on the driving time we are willing to spend along a path, but also on the initial SoC. Hence, we obtain a bivariate function f P , which maps driving time and initial SoC at u to the SoC at w. Note that consumption is positive on the first arc (u, v) and negative on the second arc (v, w) in Figure 2. As before, the second arc (v, w) provides the better tradeoff. However, for low initial SoC, we have to ensure that the first arc (u, v) can be traversed. Hence, spending some additional time on this arc may be inevitable to obtain a feasible solution. In contrast, high SoC values may prevent recuperation along the second arc (v, w), so driving slower no longer pays off at some point. Figure 3 sketches the resulting bivariate SoC function for specific values of initial SoC and driving time. For a given initial SoC b ∈ [0, M ] at u, we see how spending more time on the u-v path can increase the SoC at w (Figure 3a). When fixing the driving time x ∈ R ≥0 (Figure 3b), the optimal amount of time spent on each arc varies with the initial SoC at u. Consequently, even if we fix the total driving time that is spent on the path, the SoC function of a path P no longer has the specific form as in the case of scalar energy consumption values [16,30]. The optimal solution to EVCAS is the minimum driving time of c P for the corresponding initial SoC (which is part of the input). For example, the optimal driving time is 4 if the initial SoC is 1, or 2 if the initial SoC is at least 2.

Realistic Consumption Functions and Basic Operations
Realistic tradeoff functions are nonlinear. We introduce such tradeoff functions in Section 3.1. Although these realistic functions require a more technical analysis, many observations made for the simplistic, linear model discussed in Section 2.3 carry over to the realistic, nonlinear tradeoff functions.
If we want to incorporate the propagation of consumption functions into the BSP algorithm in order to solve EVCAS, we need to generalize some of its basic operations. First, BSP computes scalar sums to obtain the driving time and the SoC at a vertex. Since we are now dealing with continuous functions, we require operations that compute functions for (best) tradeoffs of paths instead of arcs. In the previous section, we sketched how such consumption functions can be obtained in a simplistic model based on piecewise linear functions. In Section 3.2, we describe how this is done in the realistic model. Second, the BSP algorithm identifies dominated labels to keep the number of labels small. We discuss generalized dominance schemes for consumption functions in Section 3.3.

A Realistic Consumption Model
Both considered metrics, driving time and energy consumption, depend on the vehicle's speed. In accordance with realistic physical models established in the literature [1,4,17,32,50,61,65], we assume that energy consumption on a certain road segment a ∈ A can be expressed by a function h a : R >0 → R given as where v ∈ R >0 is the (constant) vehicle speed, s ∈ R is the (constant) slope of the road segment, and λ 1 ∈ R ≥0 , λ 2 ∈ R ≥0 , and λ 3 ∈ R ≥0 are constant nonnegative coefficients of the consumption model. The term λ 1 v 2 is caused by aerodynamic drag, which increases with driving speed in a quadratic fashion. Note that energy consumption can become negative for downhill segments with s < 0. The parameters λ 1 , λ 2 , and λ 3 may vary for different arcs due to, e. g., different road types or other factors affecting energy consumption [22,82,90]. Assuming constant speed and slope per arc is not a restriction, since we can add intermediate vertices in the graph to model changing conditions. Furthermore, one can show that deliberately varying the speed along a single road segment (with constant slope and speed limit) never pays off in our model [50,Corollary 1].
Since we are interested in functions mapping driving time x ∈ R >0 to energy consumption g a (x), we substitute v = /x in Equation (4), where ∈ R >0 denotes the length of the road segment. As slope and length of an arc are fixed, we simplify this below by setting α := λ 1 2 and γ := λ 2 s + λ 3 . Observe that α ∈ R ≥0 is nonnegative, while γ ∈ R may be a negative value (for downhill arcs). and γ = 1. The indicated subdomain borders induced by its minimum and maximum driving time, respectively, areτ = 2 andτ = 6.
We introduce a third constant β ∈ R ≥0 , which we will need later to shift functions along the time axis. Altogether, we obtain the tradeoff function g a : R >0 → R mapping driving time to energy consumption, which is defined as For single arcs, we always obtain β = 0 and assume driving time x to be strictly positive. Thus, the term x − β in the denominator is strictly positive and g a (x) is a finite real value. Furthermore, note that g a is decreasing and convex on its domain R >0 in this case. Tradeoff functions of paths may require values 0 < β < x to reflect additional time spent on previous arcs. In the simplistic model discussed above, we saw that tradeoff functions of paths may be piecewise linear. Similarly, we allow tradeoff functions in the realistic model to be defined as piecewise functions, so they may consist of multiple subfunctions of the form of Equation (5). Given a tradeoff function g : R >0 → R, we plug in the minimum and maximum driving timesτ ∈ R >0 andτ ∈ R >0 to obtain the corresponding consumption function c : R ≥0 → R ∪ {∞}; see Figure 4 for an example. In general, we require consumption functions to be continuous in the interval [τ, ∞), but not necessarily differentiable. In particular, we demand that β < x holds for all subfunctions of the form of Equation (5) within their respective subdomain of [τ,τ ]. Hence, the denominator term x − β is always strictly positive. Together with the assumption α ≥ 0, this implies that consumption functions are either constant (if α = 0, in which case we further assumeτ =τ ) or strictly decreasing on the interval [τ,τ ], i. e., c( At certain points below, we also make use of the inverse function Observe that it is well-defined on the specified domain.

Linking Consumption Functions
To compute best possible tradeoffs on paths consisting of multiple arcs, we define a link operation link : F × F → F on the function space F of consumption functions as specified in Section 3.1. For two paths their concatenation. Given two consumption functions c 1 and c 2 representing energy consumption on two paths P 1 and P 2 , respectively, linking c 1 and c 2 results in a consumption function c := link(c 1 , c 2 ) that maps driving time spent when traversing the path P := P 1 • P 2 to the minimum possible energy consumption (if we ignore battery constraints for now).
the minimum and maximum driving time of c 1 , respectively. Similarly, letτ 2 ∈ R >0 andτ 2 ∈ R >0 denote the corresponding driving times of c 2 . Clearly, we obtain c(x) = ∞ for all x <τ 1 +τ 2 and c(x) = c 1 (τ 1 ) + c 2 (τ 2 ) for x >τ 1 +τ 2 . For any remaining value x ∈ [τ 1 +τ 2 ,τ 1 +τ 2 ], we have to determine times x 1 ∈ [τ 1 ,τ 1 ] and x 2 ∈ [τ 2 ,τ 2 ] that sum up to x 1 + x 2 = x and minimize overall consumption (as in the simple model discussed in Section 2). We set ∆ := x 1 below, which yields for all x ∈ [τ 1 +τ 2 ,τ 1 +τ 2 ]. In other words, to minimize the energy consumption for a given time x ∈ [τ 1 +τ 2 ,τ 1 +τ 2 ], we have to divide the amount of time that exceeds the minimum possible total driving timeτ 1 +τ 2 among the two paths, such that consumption is minimized. We would like to spend this additional amount of time x − (τ 1 +τ 2 ) on the path corresponding to the function with steeper slope, since it provides the better tradeoff (we save more energy per additional unit of time spent on the corresponding arc). This is illustrated in Figure 5. In what follows, we formally derive the link operation.

Linking Functions Defined by Single Tradeoff Functions.
For now, assume that each of the given functions c 1 and c 2 is defined by a single tradeoff subfunction, rather than multiple ones. For example, this is the case if both paths P 1 and P 2 consist of single arcs, i. e., The result c := link(c 1 , c 2 ) of linking c 1 and c 2 is a piecewise-defined consumption function, which may consist of multiple subfunctions of the form as in Equation (5); see Figure 5 for an example. Intuitively, the first subfunction of c represents a shifted part of the steeper input function for small values of x. It is followed by a combination of both functions and a subfunction that corresponds to the input function that is gentler for large values of x. Any of these parts may collapse (for example, the combined part only exists if one can pick admissible driving times such that both functions have identical negative slopes). We now show how c = link(c 1 , c 2 ) is computed. Apparently, the best choice of ∆ in Equation (6) depends on the value x ∈ [τ 1 +τ 2 ,τ 1 +τ 2 ]. Therefore, we consider the ∆-function ∆ opt : [τ 1 +τ 2 ,τ 1 + τ 2 ] → R ≥0 that maps every admissible value of x to the optimal choice of ∆. Then, we immediately get for arbitrary x ∈ R ≥0 that Hence, we essentially need to compute ∆ opt to obtain the desired function c. To this end, consider an arbitrary fixed driving time x ∈ [τ 1 +τ 2 ,τ 1 +τ 2 ]. To identify the optimal value ∆ opt (x), we examine the derivative c x of the term c x (∆) := c 1 (∆) + c 2 (x − ∆). We obtain a unique zero ∆ * x for this derivative under the assumption that β 1 < ∆ and ∆ < x − β 2 ; see Appendix A.1 for details. The value ∆ * x minimizes energy consumption for an unrestricted distribution of driving times that sum up to x. However, from Equation (6) we get the additional constraints ∆ opt ( Equations (7) and (8) together are sufficient to specify the desired function c.
In Appendix A.1, we show that the resulting function c is indeed a consumption function, i. e., a function that has the general form as in Equation (1), is continuous and decreasing on the interval [τ 1 +τ 2 , ∞), and whose tradeoff function is defined by a (small) constant number of subfunctions that have the form as in Equation (5). Degenerate cases are possible, too. For example, if c 1 or c 2 is constant, the link operation becomes a simple shift on the x-axis and y-axis. In conclusion, the link operation requires constant time in the special case where both input functions are defined by a single tradeoff function. Lemma 1 summarizes our results. Lemma 1. Given two consumption functions c 1 and c 2 , each defined by a single tradeoff subfunction as defined in Equation 5, the link operation link(c 1 , c 2 ) requires constant time and its result is a consumption function that is uniquely described by at most three tradeoff subfunctions as defined in Equation 5.
Linking General Consumption Functions. We tackle the general case, in which we are given two consumption functions c 1 and c 2 , each possibly consisting of multiple tradeoff subfunctions with the general form as in Equation (5), and want to compute the function c := link(c 1 , c 2 ). Consider a tradeoff subfunction g of c 1 and its subdomain [τ g ,τ g ). The subfunction g itself induces a consumption function c g with minimum driving timeτ g and maximum driving timeτ g ; see Equation (1). By the fact that c 1 is decreasing and by construction of the consumption function c g , we get c g (x) ≥ c 1 (x) for all x ∈ R ≥0 and c g (x) = c 1 (x) for all x ∈ [τ g ,τ g ). Since the same argument can be made for subfunctions of c 2 , it follows directly from Equation (6) that we obtain c after applying the constanttime link operation from Lemma 1 to all pairs of consumption functions induced by subfunctions of c 1 and c 2 . The lower envelope of all resulting candidate subfunctions yields the desired function c. Obviously, c is again a piecewise function. Moreover, as the lower envelope of functions that are decreasing on a common domain must be decreasing on this domain as well, the function c is decreasing. Finally, the following Lemma 2 claims that c is also continuous on its interval of admissible driving time (see Appendix A.2 for a formal proof). Therefore, the function space of consumption functions is indeed closed under the link operation. Note that the link operation is also commutative and associative. The running time of the naïve link operation described above is quadratic in the number of subfunctions of c 1 and c 2 . In Appendix A.3, we show how the complexity of the link operation for general consumption functions can be reduced to linear time. In our experiments, the number of subfunctions per consumption function was relatively small on average, so the speedup provided by the more sophisticated linear-time method was limited (less than 10 %). Hence, the algorithm presented in Appendix A.3 may rather be considered to be of theoretical interest.

Dominance Tests for Consumption Functions
If we are given multiple consumption functions corresponding to paths between the same pair of vertices, we only want to propagate those which represent optimal tradeoffs. Assume we are given two consumption functions c 1 and c 2 corresponding to two such paths, each possibly defined by multiple subfunctions (c. f. Section 3.2). We say that c 1 dominates c 2 if c 1 (x) ≤ c 2 (x) holds for all x ∈ R ≥0 . Apparently, c 2 cannot be part of an optimal solution in this case, so it would not need to be propagated by the search algorithm.
To efficiently test whether c 1 dominates c 2 , we inspect their respective subfunctions g 1 1 , . . . , g k 1 and g 1 2 , . . . , g 2 . For some i ∈ {1, . . . , k} and j ∈ {1, . . . , }, consider the pair g i 1 and g j 2 of subfunctions, both having the form of Equation (5), and let [τ i 1 ,τ i 1 ) and [τ j 2 ,τ j 2 ) denote their respective subdomains. Observe that if the subdomains do not overlap, dominance can be checked easily by inspecting function values at the endpoints of these subdomains. Hence, assume that their intersection is not empty, i. e., [τ i We can test in constant time whether g i 1 (x) ≤ g j 2 (x) holds for all x ∈ [max{τ i 1 ,τ j 2 }, min{τ i 1 ,τ j 2 }] as follows. The subfunction g i 1 dominates the subfunction g j 2 in this interval if and only if g i 1 (x) ≤ g j 2 (x) holds for the two values x = max{τ i 1 ,τ j 2 } and x = min{τ i 1 ,τ j 2 } at the borders of their subdomain intersection, as well as the unique extreme point of g i 1 − g j 2 (if it falls within the considered intersection of their subdomains; see also Appendix A.4). Since we only have to compare subfunctions if their subdomains intersect, we can test whether c 1 dominates c 2 in a single linear scan (comparing subfunctions in increasing order of driving time). Hence, the running time of a dominance test is linear in the number of subfunctions of c 1 and c 2 .
Computing Dominated Subfunctions. The dominance test described above enables pairwise comparison of consumption functions to identify dominates ones, similar to the pairwise checks between labels in the BSP algorithm. However, since we are dealing with continuous functions, a consumption function may also be partially dominated or dominated only by a set L of other consumption functions corresponding to paths to the same vertex; see Figure 6. In other words, even if c is not dominated by any other consumption function alone, it is possible that for some (or even all) admissible driving times x ∈ [τ,τ ] of c, there is a function c * ∈ L with c * (x) ≤ c(x). We discuss different strategies to efficiently identify dominated subfunctions of c.
A common approach in the literature to avoid pairwise dominance checks in related settings is to explicitly compute the lower envelope of all functions currently stored at a vertex [26,53]. This can also go in hand with a change in the label propagation strategy [73]. Unfortunately, managing lower envelopes comprising multiple functions (each consisting of several tradeoff subfunctions) is difficult in our setting, because the quadratic-time link operation described in Section 3.2 becomes It is easy to see that the lower envelope of multiple convex consumption function is not convex in general. Although one can generalize the link operation to handle non-convex functions (by essentially dividing it into convex subfunctions), label propagation becomes more expensive in theory and in practice [2]. Another important aspect that makes this approach much more complicated in practice is the fact that computing the lower envelope of multiple functions requires us to identify their intersection points, i. e., given two consumption functions c 1 ∈ L and c 2 ∈ L, we have to find driving times x ∈ R ≥0 with c 1 (x) = c 2 (x). While this equation has at most two solutions, it is relatively expensive to solve, as expanding the above equation results in a large number of terms. In addition to that, computing intersection points is prone to rounding errors and numerical instability, which can have a significant effect on maintainability and the quality of solutions in practice. Instead, we propose an alternative approach that avoids pairwise dominance checks, but does not require us to explicitly compute and maintain the lower envelope of multiple consumption functions. It tries to identify dominated parts of the function c in question. To this end, we compute a valuē τ * ∈ [τ,τ ] such that for all x ∈ [τ,τ * ], there is at least one function c * in the set L with c * (x) ≤ c(x). This valueτ * can be computed in a single (coordinated) linear scan over the subfunctions of c and all subfunctions of the functions in L. If we obtainτ * =τ in this scan, c is dominated for all x ∈ R ≥0 , so it cannot be part of an optimal solution. Otherwise, we setτ =τ * . Analogously, we then compute a valueτ * ∈ [τ * ,τ ] such that c is dominated for all x ∈ [τ * ,τ ] and setτ =τ * . Afterwards, we discard any subfunction of c whose subdomain does not intersect [τ * ,τ * ].
To compute the values ofτ * andτ * , we still face the problem that we may have to compute the intersection points of (subfunctions of) c and c * ∈ L. To avoid the costly and error-prone computation of intersection points, our implementation computes approximate valuesτ * andτ * . It uses the simple test described before to identify subfunctions that dominate a subfunction of c in the whole intersection of their respective subdomains. Clearly, we do not necessarily find the optimal values forτ * andτ * that way. Furthermore, note that there may still exist values x in the open interval (τ * ,τ * ) with c * (x) ≤ c(x) for some function c * ∈ L, which we could also aim to identify [52]. However, as argued above, the fact that we avoid the computation of intersection points makes our dominance test efficient while keeping it relatively easy to implement and maintain.

The Tradeoff Function Propagating Algorithm
We are now ready to describe our tradeoff function propagating (TFP) algorithm, which utilizes the tools from Section 3 to generalize the exponential-time BSP algorithm and solve EVCAS. For a source s ∈ V , a target t ∈ V , and an initial SoC b s ∈ [0, M ], it computes the fastest s-t path and optimal driving times such that battery constraints are respected.

Algorithm Description
Algorithm 1 shows pseudocode of TFP. The algorithm exploits the fact that the SoC at the source is part of the query input. Rather than computing bivariate SoC functions (c. f. Equation (2) in Section 2.2 and the example from Section 2.3), this enables us to propagate univariate consumption functions (defined by sequences of tradeoff subfunctions) and ensure on-the-fly that battery constraints are not violated.

Data Structures.
To keep the number of label comparisons small, we use the same method as in a previous publication [11], in which each vertex v ∈ V maintains a set L set (v) and a heap L uns (v) containing settled (i. e., extracted) and unsettled labels, respectively. Each label is a (piecewise) consumption function, mapping the driving time on a certain s-v path to energy consumption. The key of a label c at the vertex v is defined as i. e., its minimum driving time. The corresponding maximum energy consumption c(τ c ) is used to break ties. (Observe that the key does not depend on the vertex v, however, we will need v as a parameter of the key function when we modify it in Section 5.) As in Baum et al. [11], we maintain the invariant that for each v ∈ V , L uns (v) is empty or the unsettled label in L uns (v) with minimum key is not dominated by any settled label in L set (v). New labels are pushed into L uns (v). Whenever the minimum element of L uns (v) changes (because an element is added or extracted), we check whether the new minimum element is dominated by some settled label in L set (v) and discard it in this case. The algorithm uses a priority queue to determine the next vertex to be scanned. Note that unlike the BSP algorithm, our priority queue contains vertices rather than labels. This does not alter the order in which labels are scanned by the algorithm, but the number of entries in the priority queue is smaller.
In summary, we use the following data structures: The label sets L uns (·) and L set (·) store the consumption functions of a vertex (keys are calculated on-the-fly from these functions if required, so they do not have to be stored explicitly), whereas the priority queue operates on pairs of vertices and keys.
Label Propagation. Label sets and the priority queue are initialized in lines 1-5 of Algorithm 1. Initially, all label sets are empty, except for the constant function c s ≡ M − b s at the source s (the admissible driving time is 0 and the "consumption" at s equals the difference between the battery Algorithm 1: Pseudocode of the TFP algorithm. capacity and the initial SoC). The source vertex is also added to the priority queue, which uses the minimum key among any unsettled labels of a vertex as its key; c. f. Equation (9).
In each step of its main loop, the algorithm extracts some vertex u ∈ V with minimum key from the priority queue and settles the corresponding label c u by taking it from the set L uns (u) and inserting it into L set (u) (lines 7-9 of Algorithm 1). Afterwards, the priority queue and the set L uns (u) are updated accordingly (lines [12][13][14][15][16]. Finally, outgoing arcs are scanned (lines [17][18][19][20][21][22][23][24]. For every arc (u, v) ∈ A, the function c := link(c u , c (u,v) ) is computed. Note that c may violate battery constraints, so we set the minimum driving timeτ c to the smallest driving time value that yields a consumption of at most M in line 19 of Algorithm 1 (if no such value exists, we setτ c = ∞ to indicate that there are no admissible driving times). Similarly, we identifyτ c , the largest nonnegative consumption value (a consumption of 0 corresponds to a fully charge battery) in line 20. To find these values efficiently in practice, we first check whether 0 or M are in the domain of the inverse c −1 of c. If this is the case, we set the minimum driving time of c toτ c = c −1 (M ) or the maximum driving time toτ c = c −1 (0), respectively. Observe that this corresponds to setting c(x) = ∞ for all driving times x ∈ R ≥0 with c(x) > M and c(x) = 0 for all x ∈ R ≥0 with c(x) < 0. If the resulting function yields finite consumption for some driving time, it is added to L uns (v) and the key of v in the priority queue is updated, if necessary.
Note that because minimum driving times are strictly positive, the algorithm is label setting, i. e., labels extracted from the queue are never dominated by another label later on. Thus, an optimal path is found as soon as a label at t is extracted. The minimum driving time of this label is the optimal driving time; see line 11 of Algorithm 1.

Dominance Tests.
In accordance with Section 3.3, we distinguish different ways to identify dominated labels. First, our basic TFP algorithm implements the lightweight pairwise dominance tests described at the beginning of Section 3.3. For a consumption function c in the set L uns (v) of some vertex v ∈ V , we perform a pairwise comparison with each function in L set (v) to determine whether c is dominated by one of them. Recall that in doing so, the algorithm may miss that c is dominated by a set of other labels (as in Figure 6a) or partially dominated (as in Figure 6b).
Second, our improved variant computes valuesτ * ∈ R ≥0 andτ * ∈ R ≥0 as described in Section 3.3, such that for each x <τ * and x >τ * , there is at least one function c * contained in the set L set (v) with c * (x) ≤ c(x). Note that after this procedure, even if c is not discarded (because it is not dominated for all x ∈ R ≥0 ), it may no longer be the function with minimum key in L uns (v), in which case we repeat the dominance check for the new minimum element of L uns (v). Using the improved dominance check greatly reduces the number of labels in practice, at the cost of (little) additional overhead per dominance test (we now require a coordinated scan over multiple subfunctions).
Recall from Section 3.3 that we avoid exact computation ofτ * andτ * for the sake of numerical stability. To distinguish, we call this variant of our improved dominance check stable, as opposed to the exact variant that requires intersection points to find the actual values ofτ * andτ * . Note that all three variants (lightweight, improved stable, and improved exact) maintain correctness of the algorithm, but the number of labels and subfunctions the algorithm has to propagate may differ.

Path Retrieval
To obtain the actual s-t path, each label maintains two pointers to its corresponding parent label and vertex. To retrieve the optimal driving time (and speed) per arc, we explicitly store with each label c the function∆ opt , given as∆ opt (x) := x − ∆ opt (x), wrt. the previous link operation that led to the creation of c (see Section 3.2 and Appendix B.1). It is defined on the domain [τ,τ ] induced by the minimum driving timeτ ∈ R >0 and the maximum driving timeτ ∈ R >0 of c.
After the search has found the target t, the following backtracking routine yields the path itself and driving times on all arcs. We start backtracking from t. Let x =τ c denote the minimum accumulated driving time of the label c extracted at t and let v ∈ V be the parent vertex of c. Then the driving time on the arc (v, t) is∆ opt (x), where∆ opt is the function stored with the label c. We continue this procedure recursively at its parent label and for the driving time value x −∆ opt (x), until the source vertex s is reached.

A Polynomial-Time Heuristic
Even with improved dominance checks, TFP has exponential running time. However, the algorithm can easily be extended to a more efficient heuristic search, at the cost of inexact results. We propose a polynomial-time approach that is based on ε-dominance [6]. During the search, when performing the dominance test for a label c ∈ L uns (v) in the unsettled label set of some vertex v ∈ V , it is kept in L uns (v) only if it yields an improvement over settled labels in L set (v) by at least a certain fraction εM , with ε ∈ (0, 1], for some driving time. Thus, when identifying dominated parts of c, we test for each driving time x ∈ R ≥0 whether c * (x) ≤ c(x) + εM holds for some settled label c * ∈ L set (v). Lemma 3 shows that this implies that the number of settled labels per label set can become at most 1/ε +1, provided that exact improved dominance checks are applied. Given that the algorithm is label setting, each of these labels is extracted from the priority queue and added to L set (·) at most once, so this yields polynomial running time in |V | and 1/ε . For a formal proof of Lemma 3, see Appendix B.2.

Lemma 3. The number of settled labels contained in the set
when running TFP with ε-dominance and exact improved dominance checks.

A* Search
Even the heuristic variant of the TFP algorithm described above remains impractical for realistic ranges and reasonable values of ε. Therefore, this section proposes extensions of A* search to speedup TFP and its heuristic variant, whereas Section 6 discusses extensions of CH.
Potential Shifting. A well-known approach to reduce (practical) running times in multicriteria settings is the adaptation of A* search [49,67]. The key idea of this approach is to guide the search towards the target by assigning smaller keys (in the priority queue) to vertices that are closer to the target. This is achieved by making use of a potential function π : V → R on the vertices [58]. The potential function is called consistent (wrt. a cost function z : Vertex potentials are added to the keys of labels in the priority queue. Using a consistent potential function, the label-setting property of search algorithms is maintained. In Section 5.1 and Section 5.2 below, we propose two variants of A* search that compute a consistent potential function at query time, prior to running TFP. We use these potentials to modify the keys of labels when running Algorithm 1; c. f. Equation (9).

Potentials Based on Single-Criterion Search
Similar to the approach by Tung and Chew [86], our first variant aims at directing the search towards the target by preferring arcs on (unrestricted) fastest paths. Before running TFP, we execute a single backward search (i. e., Dijkstra's algorithm running on the reverse input graphḠ, which has all arc directions reversed) from the target t ∈ V . This search uses the cost functiond : A → R ≥0 withd(a) :=τ a for all a ∈ A, representing minimum driving time on an arc. This yields, for each vertex v ∈ V , the minimum unrestricted driving time from v to t, denoted distd(v, t), which is a lower bound on the constrained driving time from v to t. We obtain a potential function π d : V → R ≥0 by setting π d (v) := distd(v, t). By the triangle inequality, the potential function π d fulfills the condition Thus, TFP is label setting when using π d , which implies that correctness of the algorithm is maintained [58]. We replace the key function from Equation (9) by for a given consumption function c with minimum driving timeτ c at some vertex v ∈ V .
Improvements. To compute π d , the backward search described above visits all vertices in the (reverse) graph. This can be wasteful, especially for small vehicle ranges. To avoid scans of unreachable vertices, we make use of a second cost functionc : A → R withc(a) := c a (τ a ) for all a ∈ A, which yields the minimum energy consumption of an arc. We run a backward search onḠ using this cost function to compute, for all v ∈ V , lower bounds distc(v, t) on the energy consumption required to reach t from v. As some arc costs may be negative, we use a label-correcting variant of Dijkstra's algorithm, i. e., a variant in which the label of a vertex is not necessarily final when it is extracted from the queue (the vertex may be reinserted and extracted from the queue several times). This algorithm has exponential running time in the worst case [56], but it performs well in practice if the number arcs with negative costs is small 2 [3,11,13]. We prune the search, i. e., we do not relax any outgoing arcs, whenever the distance label of some scanned vertex exceeds the battery capacity M . Then, we run Dijkstra's algorithm onḠ to compute π d (v) as before, but restrict the search to vertices for which lower bound on energy consumption is below M . Afterwards, we restrict the TFP search to the same set of vertices. Note that a similar approach can be used to compute vertex potentials for the BSP algorithm when solving EVCSP on a graph with multi-arcs and scalar costs [12].
Additionally, we use lower bounds distc(·, ·) for pruning in TFP: Before adding a new label c to the label set of some vertex v ∈ V , we first set c(x) = ∞ for all driving times x ∈ R ≥0 with c(x) + distc(v, t) > M . To do so, we check (in a linear scan over the subfunctions defining the inverse function c −1 ) whether M − distc(v, t) is in the domain of c −1 and set the minimum driving time of c to c −1 (M − distc(v, t)) if this is the case. Otherwise, we either obtain c ≡ ∞ or the function c remains unchanged.

Potentials Based on Bound Function Propagation
The potential function π d may be too conservative if the consumption on the fastest path is very high. In such cases, it pays off to adapt the potential function [11]. The potential function π ϕ incorporates the current SoC at a vertex to provide more accurate bounds. For each vertex v ∈ V , it uses a convex, piecewise linear function ϕ v : R → R ≥0 ∪ {∞} that maps the current SoC at v to a lower bound on the remaining driving time from v to t. As higher SoC allows the EV to drive faster, ϕ v is also decreasing wrt. SoC. We say that the potential function π ϕ is consistent if the condition In other words, the above term must be nonnegative for any admissible parameter choice for the SoC function f (u,v) . We define π ϕ (u, −∞) := ∞. We further require that the potential π ϕ (t, b) at the target t ∈ V equals 0 for all b ∈ [0, M ], which implies that consistent potentials maintain correctness of TFP; see a previous work [11] for details on consistency of vertex potentials that incorporate SoC.
Computing Vertex Potentials. The functions representing the potentials π ϕ (·, ·) are determined in a label-correcting backward search from t. This search operates on the reverse graphḠ = (V,Ā) of the input graph and maintains a piecewise linear function we evaluate ϕ v by linear interpolation between two consecutive breakpoints. We ignore battery constraints in the search, so we obtain lower bounds on SoC values, which can also become negative. See Appendix C.1 for the pseudocode.
Each vertex stores a single label consisting of its (tentative) lower bound function. The search is initialized with a function ϕ t at the target t ∈ V that evaluates to 0 for arbitrary nonnegative SoC. Labels are extracted in increasing order of their minimum driving time. Given the label extracted at some u ∈ V , the search links the lower bound in this label with each outgoing arc (u, v) ∈ A. To this end, it first converts the consumption function c (u,v) mapping driving time to energy consumption to a piecewise linear function ϕ (u,v) mapping SoC to a lower bound on driving time. Letτ ∈ R >0 andτ ∈ R >0 denote the minimum and maximum driving time of Otherwise, we consider two geometric lines 1 and 2 defined as follows; see also , the lines 1 and 2 intersect in a unique point r ∈ R 2 . Then, a convex lower bound ϕ (u,v) is defined by the breakpoints [p, r, q]. To increase accuracy of the lower bound, we may repeat this operation recursively for initial points p or q together with the unique point r ∈ R 2 on the inverse of c (u,v) that has the same c-coordinate as r; see Figure 7. Recursion stops as soon as the maximum difference between c −1 (u,v) and the function ϕ (u,v) induced by the current sequence of points falls below a predefined threshold ε > 0. Note that this difference is obtained in constant time, since it always occurs in the latest point r that was added to ϕ (u,v) .
We link the function ϕ u at u with the resulting lower bound ϕ (u,v) by computing, for any SoC value b ∈ R, an optimal distribution of the available amount b of energy among the arc (u, v) and the remaining u-t path. Formally, the resulting function ϕ = link(ϕ u , ϕ (u,v) ) evaluates to This function can be computed in a linear scan over the breakpoints of both functions [11]. To improve running times in practice, we can discard the label ϕ at this point if After scanning the arc, we merge the resulting function ϕ with the function ϕ v in the label of v, i. e., we compute the pointwise minimum ϕ v = merge(ϕ v , ϕ) in a linear scan over the breakpoints of ϕ v and ϕ. To ensure that the resulting function ϕ v is again convex, we apply the linear-time scan of Graham [46]: Starting with an empty sequence of breakpoints, it iteratively appends the breakpoints of ϕ v in increasing order of driving time, each time removing previously added breakpoints from the tail of the sequence if that is necessary to maintain the function's convexity. Recursion starts at r .

Running TFP with Vertex Potentials.
Once the backward search has terminated, we define 2 for a formal proof of potential consistency). Afterwards, we start the actual TFP search, during which we obtain the key of a label c at some vertex v ∈ V by setting it to

Contraction Hierarchies
We propose an adaptation of the speedup technique CH [41] to the EVCAS problem. The algorithm distinguishes an offline preprocessing phase and an online query phase [5]. Plain CH [41] accelerates single-criterion shortest-path queries wrt. a given cost function z : A → R ≥0 by augmenting the graph with shortcut arcs. These shortcuts are computed during a preprocessing routine that iteratively contracts all vertices of the input graph in a (heuristic) order. To this end, a core graph G is maintained, initially set to G = G. When a vertex v ∈ V is contracted, it is removed from G together with all its incident arcs, while shortcut arcs are added between its neighbors to preserve distances, if necessary: Let u ∈ V and w ∈ V be neighbors of v in G . To determine whether a shortcut candidate (u, v) with cost z * := z(u, v) + z(v, w) is needed, a witness search runs Dijkstra's algorithm to compute the distance dist z (u, w) in G (after removing v, but before adding the shortcut). If dist z (u, w) ≤ z * holds, the shortcut is not required. Eventually, contracting all vertices results in an "empty" graph G = (∅, ∅). As final step of the preprocessing stage, an augmented graph is constructed, which consists of all vertices and arcs of the original graph G, together with all shortcuts that were inserted into the core graph at some point during preprocessing. Given a source s ∈ V and a target t ∈ V , the query algorithm runs a bidirectional variant of Dijkstra's algorithm on the augmented graph: a forward search from s and a backward search from t.
Each search scans only upward arcs in the augmented graph wrt. the contraction order, i. e., only those arcs where the tail vertex was contracted before the head vertex. One can show that this algorithm computes the distance between s and t in the input graph [41].
Adapting CH to the EVCAS problem. As in plain CH [41], we contract vertices iteratively during preprocessing and insert shortcut arcs to maintain distances. However, we contract only a subset of the vertices, leaving an uncontracted core graph-a common approach in complex settings [10,28,50,79]. The query algorithm then scans upward arcs and all arcs that are still contained in the core graph.
Unfortunately, the SoC at the target vertex t ∈ V is not known at query time. This makes backward search difficult, since it would require us to propagate SoC functions. Their bivariate nature makes explicit construction and comparison of SoC functions rather involved. Therefore, we first run (at query time) a breadth-first search (BFS) from t on the backward graph including shortcuts from preprocessing, scanning only downward arcs (no core arcs), i. e., arcs for which the head (wrt. direction in G) was contracted before the tail, marking each as visited. Afterwards, we execute TFP from the source vertex s ∈ V , but let it only scan upward arcs, core arcs, and downward arcs marked by the BFS. In other words, we use Algorithm 1, with the sole modification that only the arcs mentioned above are scanned.
In the remainder of this section, we discuss changes that need to be made to the preprocessing routine to adapt it to our setting. Most importantly, the SoC at a vertex is only known at query time. Hence, any shortcut would have to store a bivariate SoC function (c. f. Section 2.3). As argued before, this is rather impractical. In Section 6.1, we examine special categories of SoC functions that enable both a simple representation and efficient dominance tests. Afterwards, we describe how CH can utilize these simple SoC functions if we modify several key steps of the preprocessing routine: • We only allow vertex contraction in cases where all resulting shortcuts fall into the aforementioned categories (Section 6.2).
• We efficiently decide which shortcuts to keep if there are multiple candidates connecting the same pair of vertices (Section 6.3).
• We adapt the witness search to deal with bivariate SoC functions (Section 6.4).
Omitted proofs can be found in Appendix D. 1. In what follows, we say that an SoC function f 1 dominates another SoC function

Simple Representations of Bivariate Functions
In this case, battery constraints can render driving at high speed infeasible. On the other hand, recuperation never occurs on P . Therefore, the best speed on an arc depends only on the slope of its consumption function, but not on the position of the arc in the path (in contrast to the situation discussed in Section 2.3 and sketched in Figure 2, where the order of arcs clearly matters). Consequently, it is sufficient to first link the consumption functions and apply battery constraints only once afterwards. Lemma 4 specifies the resulting bivariate SoC function, which is represented by a single univariate (consumption) function.
where c denotes the function obtained after iteratively linking the functions c 1 , . . . , c k−1 .
Note that the functions c 1 , . . . , c k−1 can be linked in arbitrary order, since the link operation is both commutative and associative. Thus, we can easily construct a shortcut for P by iteratively contracting its internal vertices v 2 , . . . , v k−1 in any given order, each time linking the consumption functions of both incident (shortcut) arcs to compute a new shortcut.
A symmetric argument holds for paths consisting of only nonpositive consumption functions. Consequently, we can allow contraction of a vertex if all (finite) values of consumption functions assigned to any of its incident arcs have the same sign, since the resulting bivariate SoC functions can be represented efficiently. On real-world instances, where the majority of consumption values is positive, this approach already allows contraction of large parts of the graph (more than 50 % of the vertices in our tests). Nevertheless, the size of the resulting core graph is still too large to achieve significant speedups.
Discharging Paths. We discuss simple representations in more general cases, exploiting that most consumption values are positive in realistic instances. We say that a path P is discharging if the SoC on P never exceeds the (arbitrary) initial SoC, i. e., there is no prefix of P that has negative minimum consumption for any driving time. Definition 3 formalizes this notion.
Note that even though the total consumption on a discharging path must be nonnegative, it may still contain subpaths with negative consumption (as long as they are preceded by subpaths with at least as much positive consumption). Apparently, it is not necessary for us to explicitly check whether the SoC exceeds M on a discharging path, since this will never occur. Hence, making sure that battery constraints are not violated becomes easier. We now show how the SoC function of a discharging path is represented by at most two consumption functions.
Clearly, a path consisting of only arcs with nonnegative consumption values (for arbitrary driving times) is discharging. We already showed how it can be represented by a single consumption function. As a more intricate example, assume we are given a path P = P 1 • P 2 consisting of two subpaths P 1 and P 2 that can be represented by two consumption functions c 1 and c 2 . Letτ 1 ,τ 1 ,τ 2 , andτ 2 denote the respective minimum and maximum driving times of c 1 and c 2 . Moreover, assume that the consumption c 1 (x) > 0 is positive for all x ∈ R ≥0 , while c 2 (x) ≤ 0 is nonpositive for all x ∈ [τ 2 , ∞). Finally, assume that |c 1 (τ 1 )| ≥ |c 2 (τ 2 )|, i. e., the cost of P 1 is higher than the gain of P 2 for any driving time, so P is discharging. We derive the SoC function of P , represented by a positive part c + with c + (x) := c 1 (x −τ 2 ) and a negative part c − (x) with c − (x) := c 2 (x +τ 2 ). The original functions are shifted along the x-axis to simplify the analysis (note that the minimum driving time of c − is 0). Given some initial SoC b ∈ [0, M ], the positive part c + , and the negative part c − , we define the which ensures that battery constraints along P 1 are not violated if the initial SoC is b; see Figure 8. Then, the SoC function f of the path P evaluates to f ( for arbitrary x ∈ R ≥0 and b ∈ [0, M ]. Given some initial SoC, the function f first checks battery constraints on the positive part c + and links the resulting function c + b with the negative part c − . Since the underlying path P is discharging, we know that no further constraints need to be checked for c − , so the function computed by link(c + b , c − ) yields minimum energy consumption subject to driving time for the initial SoC b.

Vertex Contraction
We use our insights about discharging paths to establish a preprocessing routine for CH. For the sake of simplicity, we assume that for each arc in the graph, the energy consumption is either nonnegative for all admissible driving times or nonpositive for all admissible driving times. Note that we can always enforce this by splitting an arc a ∈ A with c a (τ a ) > 0 and c a (τ a ) < 0 into two consecutive arcs corresponding to the positive part and the (shifted) negative part of c a , respectively.
Active Vertices. Shortcuts that represent either discharging paths or paths consisting solely of arcs whose energy consumption is nonpositive for all admissible driving times are called discharging or nonpositive, respectively. In both cases, we can use the simple representations we derived above to efficiently store their SoC functions. Consequently, we only allow a vertex v ∈ V to be contracted if all new shortcuts created as part of its contraction are discharging or nonpositive. We call v active in this case. Note that the number of active vertices grows as contraction proceeds, since this results in longer shortcuts, which are more likely to consist of significant positive parts. It remains to show how to decide efficiently whether a vertex is active during preprocessing and construct the necessary shortcuts if it is indeed contracted.
Assume that at some point during preprocessing, we want to know if a vertex v ∈ V incident to two (original or shortcut) arcs (u, v) and (v, w) in the core graph is active. We have to determine whether a new shortcut (u, w) can be constructed from (u, v) and (v, w) that is either discharging or nonpositive. Clearly, a nonpositive shortcut can be constructed if and only if both (u, v) and (v, w) are nonpositive. Otherwise, we want to know whether we can construct a discharging shortcut. Let P 1 be the underlying path in the original graph represented by (u, v) and let P 2 be the path represented by (v, w). We have to decide whether P = P 1 • P 2 is a discharging path, i. e., energy consumption is nonnegative on every prefix of P , regardless of the driving time. Apparently, this can only be the case if P 1 is a discharging path itself. For P 2 , we distinguish two cases. First, if P 2 is discharging as well, to follows immediately that also P must be discharging. Second, if P 2 is not a discharging path, it must consist solely of arcs for which energy consumptions are nonpositive for arbitrary (admissible) driving times, since a shortcut (v, w) would not have been created otherwise. Hence, P 2 is represented by a single consumption function c 2 with minimum and maximum driving timesτ 2 ∈ R >0 andτ 2 ∈ R >0 , such that c 2 (x) ≤ 0 for all x ∈ [τ 2 , ∞). To test whether a discharging shortcut (u, w) can be constructed in this case, we consider the positive part c + 1 and the negative part c − 1 of P 1 with corresponding maximum driving timesτ + 1 andτ − 1 . The path P is discharging if and only if c ≥ 0, as these driving times minimize consumption on P (or any prefix of P that ends at a vertex of P 2 ).
Constructing Shortcuts. When we contract an active vertex, we have to insert new shortcuts from pairs of existing incident arcs in the current graph. To this end, we have to compute the corresponding SoC function from two given discharging or nonpositive shortcuts. As before, let v ∈ V denote the vertex to be contracted and let (u, v) and (v, w) denote two incident (shortcut) arcs representing underlying paths P 1 and P 2 in the core graph.
If both (u, v) and (v, w) are nonpositive, we simply have to link their underlying consumption functions. Otherwise, (u, v) is discharging, since v is an active vertex. Let its SoC function be defined by two consumption functions c + 1 and c − 1 . We now construct the SoC function for the path P := P 1 • P 2 , which must be discharging as well. If (v, w) is nonpositive, let its consumption function be given as c 2 . We immediately obtain the positive part c + and the negative part c − of the discharging path P , with c for all x ∈ R ≥0 . Finally, assume that P 2 is a discharging path as well, with the respective consumption functions c + 2 and c − 2 . Apparently, if we know the initial SoC, we can compute the energy consumption on P by computing link(link(link(c + 1 , c − 1 ), c + 2 ), c − 2 ) and applying battery constraints before each link operation, like in the TFP algorithm (see Section 4.1). However, we want to represent P with only two consumption functions c + and c − . Recall that the only constraint we have to check for discharging paths is whether the SoC drops below 0. Since both c − 1 and c − 2 are nonpositive for all admissible driving times, the constraint needs only to be checked for c + 1 and c + 2 (i. e., before the first and third link operation). To integrate these checks into a single new positive part c + , we first compute the consumption function h := link(c − 1 , c + 2 ). Clearly, the battery can only run empty on P 2 if this consumption function is positive for some admissible driving time (otherwise, we always gain more energy with c − 1 than is lost on c + 2 ). To distinguish, we split h into a positive part h + : R ≥0 → R ≥0 ∪ {∞} with h + (x) := max{h(x), 0} and a negative part h − : whereτ ∈ R >0 denotes the maximum driving time of h. Since h is a valid consumption function, so are both h + and h − . Observe that in case h yields nonnegative (nonpositive) energy consumption for all admissible driving times, the function h − (h + ) always evaluates to 0 on its subdomain with finite image. We derive the positive part c + of P by setting c + (x) := link(c + 1 , h + )(x) for all x ∈ R ≥0 and the negative part c − of P by setting c − (x) := link(h − , c − 2 )(x +τ ) for all x ∈ R ≥0 , whereτ ∈ R >0 is the minimum driving time of h − (we shift the function to ensure that its minimum driving time is 0). Our way of splitting the function h ensures that battery constraints are only applied to prefixes of P with positive energy consumption. The SoC function of P is obtained from c + and c − as described in Section 6.1.

Comparing Shortcut Candidates
In a bicriteria setting, vertex contraction may result in multi-arcs between neighbors of the contracted vertex. In such cases, we only want to keep shortcuts if their SoC function is not dominated by SoC functions of parallel arcs. Hence, after the contraction of a vertex, we want to delete (parts of) SoC functions of shortcut candidates that are dominated by existing functions between the same pair of vertices (and vice versa). To this end, we require efficient dominance checks for SoC functions that are either discharging, i. e., represent a discharging path, or have nonpositive energy consumption for all admissible driving times. Assume two given SoC functions f 1 and f 2 defined by the respective consumption functions c + 1 , c − 1 , c + 2 , and c − 2 (we assume that the positive part evaluates to 0 for all admissible driving times if consumption is always nonpositive). Further, we assume that our goal is to remove dominated parts of f 2 , i. e., we want to identify intervals I ⊆ R ≥0 where f 1 (x, b) ≥ f 2 (x, b) holds for all driving times x ∈ I and all values b ∈ [0, M ] of initial SoC. If both SoC functions have nonpositive consumption, each is represented by a single consumption function and we can immediately apply the dominance tests described in Section 3.3. For the remaining cases, note that a discharging SoC function cannot dominate a function with nonpositive consumption for all possible values of initial SoC. Thus, we consider the case where at least f 2 is discharging. We propose a method that may miss dominated parts of SoC functions, but requires only a linear scan over the consumption functions that define the SoC functions f 1 and f 2 . Thereby, correctness is maintained, but unnecessary shortcuts may be inserted. Let ε ≥ 0 denote the smallest nonnegative real value such that c − 1 (x) ≤ c − 2 (x) + ε holds for all x ∈ R ≥0 . This value can be computed in a linear scan over the subfunctions of c − 1 and c − 2 , similar to a dominance test. Lemma 5 claims that if c + 1 (x) + ε ≤ c + 2 (x) holds for some x ∈ R ≥0 , choosing a driving time of x for the positive part c + 2 always results in a solution that is dominated by f 1 , regardless of the initial SoC.

Lemma 5.
Given a nonpositive or discharging SoC function f 1 and a discharging SoC function f 2 , such that their respective consumption functions are c + 1 , c − 1 , c + 2 , and c − 2 , let the value ε ≥ 0 be defined as described above. If c + 1 (x + ) + ε ≤ c + 2 (x + ) holds for some x + ∈ R ≥0 , any solution where x + is the (optimal) amount of time spent for c + 2 is dominated by f 1 , i. e., we obtain either c + 2 (x + ) = ∞ or After creating a new shortcut (u, v) with positive part c + , we compare it to existing shortcuts between u ∈ V and v ∈ V as follows. First, we compute the value ε defined above wrt. each existing shortcut. Then, we determine parts of c + that are dominated by existing positive parts (after we increase their consumption by ε) and set c + (x) = ∞ for such values x ∈ R ≥0 . We do this in a coordinated linear scan over c + and the positive parts of consumption functions of all existing shortcuts, as in our basic dominance tests (see Section 3.3). If c + ≡ ∞ holds afterwards, we remove the shortcut. Analogously, we identify parts of SoC functions of existing shortcuts that are dominated by the SoC function of the new shortcut candidate.

Witness Search
Consider a discharging shortcut candidate (u, v) from u ∈ V to v ∈ V that is neither dominated by any existing parallel shortcut (from u to v) nor dominates an existing parallel shortcut itself. Before adding (u, v) to the graph, we run a witness search (as in plain CH) to test if the shortcut is necessary to maintain distances in the current graph G (for some values of driving times and SoC). An exact approach would compute bivariate SoC functions representing u-v paths in G to identify dominated parts of the shortcut candidate. Like before, this is difficult and potentially expensive due to the lack of efficient operations for construction and comparison of such SoC functions. Instead, we only compute univariate upper bounds on energy consumption during witness search. As a key idea, the search drops negative parts from SoC functions entirely, so labels in the search consist of only the positive parts. Clearly, these labels are upper bounds on energy consumption. Our witness search runs the basic TFP algorithm from u on G , but ignores battery constraints and links labels only with positive parts of arcs (we assume that the positive part evaluates to 0 for all admissible driving times if consumption on the corresponding arc is always nonpositive). As a result, each label stores an upper bound on overall (finite) consumption that is independent of initial SoC and represented by a single consumption function c + ; see Figure 9. Moreover, since the majority of arcs have positive consumption in realistic instances, the upper bounds have relatively small error in most cases.
Witness Comparisons. Whenever a label c + at the head vertex v of the shortcut candidate a := (u, v) is extracted during the witness search, we compare the witness c + to the SoC function f a of the shortcut candidate as follows. Let c + a and c − a denote the consumption functions defining f a .
Moreover, let c − a (τ − ) be the minimum consumption of c − a . We know that f a cannot be a (unique) optimal choice for a driving time x ∈ R ≥0 if c + (x) ≤ c + a (x) + c − a (τ − ), i. e., the upper bound c + provides better or equal consumption than a lower bound on the consumption of f a ; see Figure 9. We proceed along the lines of the dominance tests described in Section 3.3 to identify such values x ∈ R ≥0 and remove them from the subdomain of admissible driving times of c + a . If c + a ≡ ∞ holds afterwards, the shortcut is not required. Otherwise, the witness search stops once the minimum driving time of a scanned label exceeds the maximum driving time of the shortcut candidate. During the search, we can discard labels if their minimum consumption exceeds the maximum consumption c + a (τ + ) + c − a (τ − ) of the lower bound of the shortcut candidate.

CHAsp
It is straightforward to combine our variants of A* search and CH to obtain our fastest algorithm, which we refer to as CHAsp (CH, A*, Adaptive Speeds). Queries are answered by first computing a potential function on the CH search graph and afterwards running TFP (or its heuristic variant) on it, utilizing the corresponding key function. Note that CHAsp does not alter the output of the underlying basic algorithm, so correctness is maintained when using plain TFP. Running time remains exponential in the worst case, but we observe significant speedups in practice (see Section 8). Moreover, CHAsp can be combined with our polynomial-time heuristic in just the same way. In summary, our algorithm consists of several alternative building blocks and both its preprocessing and query stage comprise different steps, which we recap in the overview given below. For implementation details, see also Appendix E.
Preprocessing. In the preprocessing stage, we are given the input graph G = (V, A) and the consumption functions of all arcs, on which we perform the following steps.
1. Determine a heuristic vertex order by assigning a priority to each vertex (see Geisberger et al. [41] and Appendix E) and identify active vertices (see Section 6.2).
2. Initialize the core graph by setting it to G = G. Contract active vertices in the given order, until no active vertices are left or the average degree of active vertices reaches a certain threshold. When contracting a vertex v ∈ V , • construct shortcut candidates for each pair of neighbors of v and compare each with any existing arcs in the core graph between the same pair of neighbors (see Section 6.3); • run witness searches between pairs of neighbors with any remaining shortcut candidates to identify unnecessary candidates (see Section 6.4); • remove v from G , add remaining shortcut candidates to G , update the priority and activity status of all neighbors of v. 3. Output the contraction order and the graph G enriched with all shortcuts that at some point were added to or are still present in G (arcs still present are called core arcs).
Queries. After completion of the preprocessing stage, queries consist of a source vertex s ∈ V , a target vertex t ∈ V , and an initial SoC b s ∈ [0, M ]. A query is answered by executing the steps listed below.
1. Run a BFS from both s and t, scanning only upward and downward arcs (wrt. the contraction order), respectively, and no core arcs. Build a search graph G s,t consisting of all arcs visited by either BFS and all core arcs.
2. Compute either the vertex potential π d or π ϕ by running the corresponding backward search from t on G s,t (see Section 5).
3. Run TFP (see Section 4.1) or its heuristic variant (see Section 4.3) from s on G s,t , but use the respective key function from Section 5.

Experiments
We implemented all approaches in C++, using g++ 7.3.1 (flag -O3) as compiler. Experiments were conducted on a single core of a 4-core Intel Xeon E5-1630v3 clocked at 3.7 GHz, with 128 GiB of DDR4-2133 RAM, 10 MiB of L3 cache, and 256 KiB of L2 cache. We consider two graphs representing the road networks of Europe with 22 198 628 vertices and 51 088 095 arcs, and Germany with 4 692 091 vertices and 10 805 429 arcs, kindly provided by PTV AG. 3 Applying elevation data from the Shuttle Radar Topography Mission, v4.1, 4 we derived realistic energy consumption functions from two detailed micro-scale emission models [51]. The first is a vehicle model that is calibrated to a Peugeot iOn. The second is an artificial model [84] that, in contrast to the first, takes power demand of auxiliary consumers (e. g., air conditioning) into account. These data sources are proprietary, but enable evaluation on detailed and realistic input data. We extracted functions following Equation (4) from given samples of speed and energy consumption via regression. Combining reasonable minimum speeds for different road types (e. g., 80 km/h on motorways and 30 km/h in residential areas) with the posted speed limits (if higher), we get intervals of admissible speeds per road segment. As a result, 25 % and 38 % of the arcs allow adapting the speed for the network of Germany and Europe, respectively. The respective average speed interval lengths among those arcs are 17 km/h and 27 km/h for the two instances. The notable differences between the two instances can be explained by a larger percentage of arcs corresponding to local or residential roads in the Germany instance, which only allow for little speed adaptation, but also by differences in road categories and speeds of the different countries in the input data in general. We denote our instances by Germany-Aux (Ger-AX), Germany-Peugeot (Ger-PG), Europe-Aux (Eur-AX), and Europe Peugeot (Eur-PG). They have negative consumption (at least for maximum driving times) on 7.8 % (Ger-AX) to 12.9 % (Eur-PG) of their arcs. Table 1 gives an overview of the different input instances. Battery capacities are measured in kWh, where 1 kWh corresponds to a range of roughly 5-10 km (depending on speed and terrain). Most experiments use realistic ranges of 16-64 kWh, but we resort to queries of shorter range when necessary to evaluate slower baseline approaches.
Unless mentioned otherwise, our study evaluates random in-range queries, i. e., we pick a source vertex s ∈ V uniformly at random. Among all vertices in range from s with an initial SoC b s = M , we pick the target t ∈ V uniformly at random. Since unreachable targets can be detected by backward search phases of A* search or by any algorithm for computing energy-optimal routes [13,30,74], this results in more difficult and interesting queries for us.

Model Validation
We argued that an approach based fully on consumption functions unlocks both better tractability and improved solution quality compared to discrete speeds. Therefore, we also consider instances with multi-arcs to model speed adaptation-as was best practice in previous approaches [12,45,50]. We generate multi-arcs in a rather conservative way, by sampling consumption functions at velocity steps of 20 km/h. Optimal paths (wrt. to the simple model) are then computed by the BSP algorithm described in Section 2.1. Indeed, we observe a significant speedup by simply switching to our more realistic model, as Table 2 shows. Even for a small range (2 kWh), we see that TFP is three orders of magnitudes faster than BSP, since it evaluates speed-consumption tradeoffs more fine-granularly and with less overhead: Instead of huge discrete Pareto sets, each route is represented by a single continuous function represented by few parameters. As a result, the number of labels settled by TFP (each representing a distinct route) is up to three orders of magnitudes less compared to the labels scanned by BSP (each representing a distinct combination of speed samples on some route). Pairwise label comparisons even drop by more than five orders of magnitude.
Regarding query statistics, note that the number of arc relaxations coincides with the number of link operations. BSP spends almost an order of magnitude less time per operation (essentially, only Table 3.: Running BSP with varying speed steps (Ger-PG, 2 kWh). Using A*-π d to render query times feasible, we ran 100 queries of BSP for different step sizes, resulting in the indicated number of arcs. We also consider instances with two multi-arcs per tradeoff function (step size ∞), the original instance without speed adaptation (cnst.), and continuous functions (step size 0). We report average label scans and comparisons, running time and standard deviation, as well as average and maximum quality loss of solutions (due to discretized speeds).
Step  Table 2. Even though the coordinated sweep comes at some overhead, the average time per dominance check actually decreases, as more dominated (sub)functions are detected and the number of required comparisons is reduced by more than a factor of 3. Linking becomes somewhat faster for the same reason (labels store fewer subfunctions). Overall, dominance tests take more than 90 % of the running time of BSP. For TFP, linking and dominance tests require roughly 15 % and 70 % of the overall running time, respectively. Although atomic operations (linking and comparing labels) are more expensive for TFP, the drastic reduction in the number of vertex scans explains the speedup. This is interesting, as sampling was expressly considered to manage tractability [12,45,50,81]. We observe that the average query time of all approaches is mostly determined by outliers, which is not surprising given that their asymptotic complexity is exponential. As a result, most queries are actually answered much faster than the average, with the median running time being about a factor of 5 lower on average. Table 3, we evaluate the effect of different step sizes on the performance of BSP and solution quality (see Appendix F for corresponding plots of query times and result quality). Since including more multi-arcs can slow down the algorithm significantly, we use A* search with the potential function π d ; see Baum et al. [12] and Section 5.1. Each row of the table shows the result of running BSP (with potentials) on the graph resulting from sampling each tradeoff function at the indicated rate (1-20 km/h). We also provide results when each function corresponds to exactly two samples at its minimum and maximum driving time, respectively (step size ∞), on the unaltered original graph (cnst.), and when running TFP on continuous functions (step size 0). In instances with multi-arcs, we always include the arc that minimizes energy consumption to maintain reachability of the target (by contrast, only 77 % of the targets could be reached on the unaltered original instance). Note that using A* search comes at the cost of side effects: In cases where the target is (almost) reachable on the unconstrained fastest route, potentials ensure that BSP and TFP (almost) only scan vertices on the shortest route. Therefore, the median query time on every instance considered in Table 3 is below 10 ms (not reported in the table). Furthermore, the average time of faster approaches with small search spaces, such as TFP, is dominated by the backward search (more than 90 % of its total running time).

Speed Resolution. In
On average, TFP is significantly faster than BSP, even if only two multi-arcs are added per tradeoff function, and several orders of magnitudes faster for sensible speed steps. At the same time, it finds paths that are up to 9.1 % quicker (within battery constraints). While quality loss quickly decreases with smaller speed steps (sampling 10 km/h steps suffices to keep it below 1 %), its running time makes BSP impractical: The number of arcs in the graph less than doubles for step sizes up to 10 km/h, but the number of labels can grow exponentially even along a single route (see Section 2.1).
Varying Speed Interval Lengths. Minimum speeds of tradeoff functions in our instances are chosen such that driving slower would become unprofitable or cause the EV to impede traffic flow. However, as adapting the driving speed pays off the most on faster roads, one could consider increasing minimum adaptive speeds. Below, we evaluate such modifications, by either increasing the minimum speed on every nonconstant arc to shrink its speed interval to a certain fraction of its original length (0 %, 25 %, or 50 %) or by only allowing adaptive speeds up to a certain lower threshold (60 km/h, 80 km/h, or 100 km/h). Table 4 shows the results of 100 queries for each setting, chosen such that the target is in range for even the most restrictive setting (0 %, i. e., constant speed on all arcs). As in the previous experiment, we use A* search with the potential function π d to enable a longer range of 16 kWh for this experiment. As before, this affects the distribution of query times (the median is 300-400 ms in all cases), but we still observe differences in average performance. Interestingly, shrinking the interval actually increases query times in most cases: Shorter intervals imply that more functions are required to cover the whole Pareto front, so more labels have to be propagated. At the same time, solution quality consistently deteriorates, as there are fewer speed options. For the first category of instances (shrinking intervals of admissible driving times), the number of tradeoff functions remains unchanged (about 25 % of all arcs; see Table 1). The second category, however, adjusting the minimum adaptive speed means that arcs for which the maximum speed does not exceed this minimum become constant.
Since the majority of arcs in a road graph represent slower roads, the number of (nonconstant) tradeoff functions decreases significantly to 9 % (3.5 %, 0.5 %) of all arcs for a minimum adaptive speed of 60 km/h (80 km/h, 100 km/h). This also explains why query times decrease again if we only allow speed adaption beyond 100 km/h: Most labels are tuples of constant values and dominance tests become simpler. A similar effect explains why compared to intervals shrunk to 25 % of their original size, the number of comparisons when using only constant labels (first row of the table) increases by a factor of more than 5, but the query time just is slightly higher. Shrinking the admissible intervals of all tradeoff functions consistently deteriorates solution quality. Compared to our basic setting (100 %), results are up to 50 % worse. Interestingly, we observe almost no quality loss if we employ a lower adaptive speed limit of 60 km/h. Query times still increase, indicating that adaptive speeds on slower roads are not needed for result quality, but to remove unimportant tradeoffs from label sets. Figure 10 plots running times and quality for the different settings shown in Table 4. For comparison, we also include extended speed intervals (compared to our basic instance). Larger intervals actually improve performance further, but running times and quality improvements converge (and we argued that speed recommendations below reasonable thresholds are not meaningful in practice).

CHAsp
In what follows, we focus on the evaluation of different variants of our fastest approach, CHAsp. A comparison of CHAsp and basic algorithms is presented in Section 8.3. Table 5 shows details on CH preprocessing effort and its impact on query performance subject to different core sizes on Ger-PG, assuming a battery capacity of 16 kWh. Vertex contraction was Table 5.: Impact of core size on performance (CHAsp, Ger-PG, 16 kWh). Vertex contraction stopped once the average degree of active vertices in the core reached a given threshold (Ø Deg.). We report the resulting core size (# Vertices) and remaining active (i. e., contractable) vertices (# Active), preprocessing time, and average query times of 1 000 queries using CHAsp with potential functions π d and π ϕ , respectively. stopped as soon as the average degree of active (i. e., contractable) vertices in the core reached the indicated threshold (Ø Deg.). We report resulting core sizes and preprocessing time, as well as query times of our fastest exact algorithms. We see that contraction becomes much slower beyond a core degree of 32, due tp the small number of remaining active vertices: Only 58 796 out of 305 301 remaining vertices in the core are active when the average degree reaches 32. This also explains why the speedup compared to the baseline (a threshold of 0 for the average degree of active core vertices yields plain TFP combined with A* search) is much smaller than in single-criterion CH with scalar arc costs [41]. Similar deteriorations in speedup were observed in other complex problems, such as time-dependent profile computation [7], time-dependent aircraft flight planning [21], and multicriteria routing [38]. Nevertheless, CH still yields an improvement by up to an order of magnitude in our case.
In all experiments below, we picked an average core degree of 32 as stopping criterion of CH preprocessing, as it yields best results not only in terms of average query times, but also all other considered time measures (median, standard deviation, maximum; not reported in the table). The resulting core size depends on different parameters, including vehicle range and error thresholds (of heuristic variants). Relative core sizes thus vary between 2.8 % for Ger-AX and 8.5 % for Eur-PG, which is explained by the difference in the amount of arcs with negative consumption. Recall that this has a significant impact on the number of active vertices and the contraction order (see Section 6).
Exact CHAsp. Table 6 shows the performance of CHAsp on all four considered instances. We report vertex scans and dominance tests of the forward search only, excluding the A* backward search (query times include both the forward and the backward search, though). Timings of only the A* backward search of each approach are shown for comparison.
Preprocessing time mostly depends on the graph size, taking about 6 times longer on Europe, but not on the vehicle model: The inclusion of auxiliary consumers reduces the number of negative arcs and hence, contraction can proceed faster, but it also takes longer to reach the threshold degree in the core graph. As a result, the core graph size (not reported in the table) is reduced by almost a factor of 2 after about the same computation time. Table 6.: Preprocessing and query performance of CHAsp. For each instance, we provide CH preprocessing times and query times of exact CHAsp, using the potential functions π d and π ϕ , respectively. Query figures are average values of 1 000 in-range queries (999 queries on Eur-PG for a range of 64 kWh and the potential π d , indicated by an asterisk). We show timings of the A* backward search, the number of settled labels (# Labels) and label comparisons (# Cmp.) during the forward search, and total query timings. Regarding queries, we observe significant differences in performance, depending on the input instance and vehicle range. For the artificial model, we achieve quite practical times in the order of milliseconds for a realistic range of 16 kWh and less than a second on average for the long range of 64 kWh. For the Peugeot model, on the other hand, average running times exceed tens of seconds and even minutes, with outliers well beyond an hour of running time, depending on the vehicle range. This gap in running time between both consumption models is explained by the difference in the number of arcs with negative cost. One could argue that the instances Ger-PG and Eur-PG are rather excessive in this regard, by not accounting for any auxiliary consumers at all. As a result, these instances are significantly more difficult to solve for our algorithms.
The search space is consistently smaller when using the potential function π ϕ , but the backward search is more expensive. In fact, it becomes the major bottleneck for a battery capacity of 16 kWh on the easier instances. Consequently, average query times are slowed down by about a factor of 4 in this case. For the harder instances (no auxiliary consumers), the potential function π ϕ provides better results due to its better scalability and more stable query times. This becomes evident in particular for the longer range of 64 kWh: On Eur-PG, one query did not terminate within a preset time limit of ten hours when using the potential π d . Hence, all reported statistics only take the remaining 999 queries into account (indicated by the asterisk in Table 6).
Average running times are affected by outliers in most cases, with median query times being up to two orders of magnitude lower. Comparing the two road networks Germany and Europe, we achieve slightly better running times on Europe for the medium range (explained by differences in the length of admissible speed intervals), but Germany is easier to solve for long ranges (searches Table 7.: Performance of the heuristic variant of CHAsp-π d , for different choices of the parameter ε on Ger-PG and Eur-PG. We show figures on query performance for the same 1 000 random queries as in Table 6. Additionally, we report the percentage of feasible (F.) and optimal (O.) results, as well as the average and maximum relative error of all queries in which a feasible solution was found. are more likely to reach the border of the graph). In summary, we solve EVCAS optimally in less than a second on average for typical ranges, even on hard instances. For long ranges, our algorithm computes the optimal solution within seconds on the most difficult instance when using the potential function π ϕ , despite its exponential worst-case running time.
Heuristic CHAsp. We also observe a drop in query times: For a range of 16 kWh, we achieve a considerable speedup of an order of magnitude. Regarding quality of results, the choice of ε clearly matters. For ε = 0.01, the decrease in quality is negligible, but speedup (about a factor of 2) is moderate. For ε = 0.1, on the other hand, the optimal solution is still found in many cases. The average error is roughly 0.2 %, while the overall maximum is 5 %, which is acceptable in practice. Finally, for ε = 1.0, both the average and maximum error increase significantly. Given that speedup is limited compared to the case ε = 0.1, we conclude that the latter provides the best tradeoff in terms of quality and query performance. Providing high-quality solutions, it enables query times of well below 100 ms, which is fast enough even for interactive applications. Moreover, note that in cases where no path is found (about 1 % of all queries for ε = 0.1), a simple fallback could return the energy-optimal path, which can be computed quickly [13,30,74].
Considering the longer range of 64 kWh, the choice of ε as a percentage of the battery capacity affects (relative) preprocessing and query times, which decrease to a smaller fraction of the exact  Figure 11.: Running times of CHAsp subject to Dijkstra rank. We show results for both potential functions π d and π ϕ . For each rank, we consider 1 000 random queries on Eur-PG, assuming a battery capacity of 16 kWh.
variant (dominance checks become coarser, allowing more shortcuts and labels to be pruned). We achieve a query speedup by more than two orders of magnitude and observe milder outliers. For ε ≥ 0.1, query times of a couple of seconds and less are quite practical. However, quality suffers from the coarser dominance relaxation and we observe average errors of around 1 %, with outliers exceeding 10 % in quality loss. Depending on the application, a value ε ∈ [0.01, 0.1] provides the most reasonable tradeoff between running time and quality. For the easier instances Ger-AX and Eur-AX (not reported in the table), we generally observe smaller errors, but also a smaller speedup. This is due to the fact that the A* backward search quickly becomes the bottleneck when combined with a heuristic variant of TFP. Similarly, using the more sophisticated potential function π ϕ does not pay off (even for the harder instances). Finally, recall that our actual implementation does not utilize exact dominance checks for the sake of avoiding expensive and numerically instable intersection tests. Although this means that the bound established in Lemma 3 does not hold in theory, it was not exceeded once in our experiments.

Evaluating Scalability
To assess the scalability of CHAsp, we consider two different input parameters that have a big impact on the size of the search space: the distance between the source and the target of a query, and the range of the vehicle. The first experiment presented below evaluates our fastest exact algorithms following the methodology of Dijkstra ranks, defined for a query as the number of vertex scans when running Dijkstra's algorithm with costs representing unconstrained driving time [5,76]. Thus, higher ranks correspond to harder queries (of longer distance). In the second experiment, we compare the different approaches discussed in this work and examine how their query times depend on the vehicle's range.
Query Times Subject to Dijkstra Rank. Figure 11 shows results for our fastest exact approaches on Eur-PG, assuming a battery capacity of 16 kWh. We ran 1 000 random queries per Dijkstra rank (note that these are not necessarily in-range queries). Queries of higher rank correspond to a longer Battery capacity [kWh] CHAsp-ε-π d Figure 12.: Median running times for different battery capacities. For each considered capacity (ranging from 0.25 kWh to 1 024 kWh), the plot shows the median running time of 100 random in-range queries, provided that no query exceeded an hour of computation time. We evaluate the BSP algorithm using multi-arcs corresponding to speed samples, A* search for BSP with the potential function π d (BSP-A*), the TFP algorithm using pairwise dominance tests, TFP with improved dominance tests (TFP-d), and our proposed speedup techniques (A*-π d , CHAsp-π d , and CHAsp-π ϕ ). We also show our heuristic approach with parameter choice ε = 0.1, denoted CHAsp-ε-π d .
distance between source and target. It turns out that median running times of both CHAsp-π d and CHAsp-π ϕ are quite robust towards varying Dijkstra ranks. We obtain the most expensive queries at ranks from 2 17 to 2 19 . (Note that random in-range queries are likely to be among these most difficult ranks.) For higher ranks, the target is often unreachable. In most cases, this is detected by the backward searches for potential computation, as lower bounds on consumption exceed the battery capacity (see Section 5). We could achieve further speedup for high ranks by running any technique that quickly computes energy-optimal routes to detect unreachable targets [13,30,74]. For lower ranks (i. e., more local queries), the target is often reachable on an unconstrained shortest path, so goal direction of the potential functions works very well. In such cases, the backward phase of A* search becomes the major bottleneck of the query, hence the lightweight potential π d yields better query times. Note that query times vary even for these very low or high ranks, as they largely depend on the portion of the core graph that is inspected by the backward search. Although the median running time of CHAsp-π ϕ is consistently higher than the median of CHAspπ d , the former is also more robust in that it produces fewer outliers. Its more sophisticated potential function pays off especially for harder queries. Nevertheless, as worst-case running time is exponential, we observe a few outliers that exceed the median by several orders of magnitude.
Query Times Subject to Vehicle Range. While query times are relatively stable for different Dijkstra ranks, the range of an EV has a major influence on performance. We evaluate our approaches for different battery limits in Figure 12. For every considered battery capacity, we report the median running time of 100 random in-range queries up to the range at which at least one query did not terminate within one hour. For small capacities, this enables our basic approaches, which are shown in the plot as well. We also evaluate the BSP algorithm, using multi-arcs to model speed adaptation (20 km/h step width). As before, we observe a considerable speedup by switching to Battery capacity [kWh] CHAsp-ε-π d Figure 13.: Maximum running times for different battery capacities. For the same sets of queries as in Figure 12, this plot shows the corresponding maximum running times of 100 random in-range queries (unless this maximum exceeded an hour of computation time). our more realistic model, which is based on consumption functions. As discussed in Section 8.1, a vast reduction in the number of label scans more than compensates for the more expensive basic operations when using TFP. It achieves a speedup by up to two orders of magnitude over BSP. Plugging in the improved dominance checks described in Section 4 pays off as well, as it yields further speedup for battery capacities beyond 1 kWh. Adding A* search, we achieve reasonable median running times of about a second for capacities of up to 32 kWh, without any preprocessing. Our technique CHAsp-π d adds preprocessing to provide further speedup by another order of magnitude. Matching our previous observations, median running times of CHAsp-π ϕ are slower for all ranges up to 32 kWh. However, this algorithm is more robust against outliers and is the only exact method that terminates within an hour for all queries at 64 kWh and up. As a result, we are able to compute provably optimal results for (hypothetical) ranges of up to 512 kWh (around 3 000 km) in less than an hour. Finally, our heuristic variant scales well with vehicle range. Query times actually bottom out for large battery capacities, as the vehicle range gets close to the graph diameter (for 1 024 kWh, the whole graph is always reachable). Figure 13 shows maximum running times for the same sets of queries (see Appendix F for a plot showing average running times). Naturally, these outliers are subject to noise, so curves are less smooth. Our most robust technique, CHAsp-π ϕ , starts to pays off at a range of 16 kWh (at lower ranges, its backward search is too expensive). Interestingly, the maximum running time of the heuristic variant of CHAsp drops to about a second for the longest considered range: Since the target is always reachable with either no or only little speed adaptation, goal direction makes the forward search find the target almost immediately. In case of exact CHAsp-π ϕ , the expensive backward search prevents the algorithm from always terminating within one hour. In practice, it would therefore pay off to first check whether the target is reachable on the fastest path (without speed adaptation) before resorting to either CHAsp-π d or CHAsp-π ϕ . The better choice between the two potentials depends on the vehicle range and by how much the target is missed on the unconstrained fastest route.

Conclusion
We introduced a novel framework for computing constrained shortest paths for EVs in practice, using realistic consumption models. Our TFP algorithm respects battery constraints and accounts for adaptive speeds in a mathematically sounder way that unlocks both better query performance and improved solution quality when compared to previous approaches using discretized, sampled speeds. Nontrivial speedup techniques based on A* search and CH make the algorithm practical. It computes optimal solutions in less than a second on sensible instances, making it the first practical exact approach-with running times similar to previous inexact methods [12,45,50]. Our own heuristic enables even faster queries while retaining high-quality solutions.
The result of our computations is not only the suggested route from source to target but also optimal driving speeds along that route. In practice, these can be passed to the driver as recommendations or directly to a cruise control unit. With the advent of autonomous vehicles, the output of our algorithms can also be utilized for speed planning of self-driving EVs, either directly or after refinement [33]. The CHAsp algorithm can further be used as a subroutine in more complex scenarios, such as fleet management or load-balanced routing. For the latter, computing several Pareto optimal routes could provide a set of alternative candidates to choose from [12]. Another next step would be the the integration of planned charging stops [11,79]; see also Niklaus [71].
To enable faster heuristic variants, it would be useful to precompute potentials for A* search, as in ALT [43,44], or integrate other known approximate or heuristic approaches [6,12,87]. We are also interested in the integration of variable speed limits imposed by, e. g., time-dependent travel times based on historic knowledge of traffic patterns [7,14,26,35]. From a more theoretical perspective, efficient representation and comparison of (general) bivariate SoC functions is an open issue. One could also extend EVCAS by asking for an optimal solution for every initial SoC wrt. given source and target vertices, similar to profile queries in time-dependent or energy-optimal routing [7,13].

A. Omitted Details and Proofs from Section 3
In what follows, we provide technical details for linking consumption functions in the special case considered in Section 3.2, in which both functions represent arcs (Appendix A.1), give a formal proof of Lemma 2 (Appendix A.2), and show how general consumption functions can be linked in linear time (Appendix A.3). Finally, Appendix A.4 mentions technical details omitted from Section 3.3.

A.1. Linking Functions Defined by Single Tradeoff Functions
We show how c = link(c 1 , c 2 ) can be computed in constant time in the case where each of the two given consumption functions c 1 and c 2 is defined by a single tradeoff subfunction, rather than multiple ones (c. f. Section 3.2). Recall Equation (6) from Section 3.2 restated below: As argued in Section 3.2, the best choice of ∆ depends on the value x ∈ [τ 1 +τ 2 ,τ 1 +τ 2 ]. Therefore, we consider the ∆-function ∆ opt : [τ 1 +τ 2 ,τ 1 +τ 2 ] → R ≥0 that maps every admissible value of x to the optimal choice of ∆. Then, we immediately get for arbitrary x ∈ R ≥0 that Hence, we essentially need to compute ∆ opt to obtain the desired function c. To this end, consider an arbitrary fixed driving time x ∈ [τ 1 +τ 2 ,τ 1 +τ 2 ]. To identify the optimal value ∆ opt (x), we examine the derivative c x of the term c x (∆) := c 1 (∆) + c 2 (x − ∆). It evaluates to where α 1 , β 1 , α 2 , and β 2 are the corresponding coefficients in the tradeoff functions of c 1 and c 2 . We assume that α 1 > 0 and α 2 > 0 are positive (the other cases α 1 = 0 or α 2 = 0 are trivial: at least one of the consumption functions is constant in this case and hence, has a unique admissible driving time). Then, we obtain a unique zero ∆ * x for this derivative under the assumption that β 1 < ∆ and ∆ < x − β 2 . This holds true for valid choices of ∆, because ∆ ≥τ 1 > β 1 and ∆ ≤ x −τ 2 < x − β 2 always holds; see Equation (6). Therefore, solving the equation c x (∆) = 0 for ∆ yields The value ∆ * x minimizes energy consumption for an unrestricted distribution of driving times that sum up to x. From Equation (6) we get the additional constraints ∆ opt (x) ≥ max{τ 1 , x −τ 2 } and ∆ opt (x) ≤ min{τ 1 , x −τ 2 }. Since ∆ * x is the unique zero of c x in the open interval (β 1 , x − β 2 ), monotonicity of c x in the intervals (β 1 , ∆ * x ] and [∆ * x , x − β 2 ) follows. Thus, we get Together with Equations (10)- (13), Lemma A-1 enables us to construct the desired function c = link(c 1 , c 2 ), depending on the slopes (i. e., the derivatives) of c 1 and c 2 at their respective subdomain borders. Exploiting thatτ 1 ≤τ 1 andτ 2 ≤τ 2 hold by definition, we obtain the function c after the following case distinction. As in Lemma A-1, let g 1 denote the (unique) tradeoff function defining c 1 and similarly, let g 2 be the tradeoff function defining c 2 . We consider the slopes of these tradeoff functions at certain subdomain borders of c 1 and c 2 . Without loss of generality, assume that g 1 (τ 1 ) ≤ g 2 (τ 2 ) holds (the other case is symmetric). This leaves us with only three possible cases, which are presented below.
Due to convexity of both g 1 and g 2 , no other cases remain. Hence, the function c constructed above is defined by at most five subfunctions (two of which are constant). In each expression, we can expand the functions c 1 and c 2 to obtain a term that has the general form of a tradeoff function as in Equation (5). In particular, the denominator in the tradeoff functions of both c 1 and c 2 is (strictly) positive in all cases, i. e., we always have x > β. Moreover, it is easy to verify that c is continuous and decreasing on the interval [τ 1 +τ 2 , ∞), by inspecting the corresponding limits at the endpoints of each subdomain of c. In conclusion, the link operation requires constant time in the special case where both input functions are defined by a single tradeoff function. Lemma 1 from Section 3.2 summarizes these insights. It is restated below. Lemma 1. Given two consumption functions c 1 and c 2 , each defined by a single tradeoff subfunction, the link operation link(c 1 , c 2 ) requires constant time and its result is a consumption function that is uniquely described by at most three tradeoff subfunctions.

A.2. Proof of Lemma 2
Lemma 2 from Section 3.2 claims that linking two consumption functions results in a function that is continuous on its interval of admissible driving time. The proof is given below. Proof. Assume for contradiction that c has a discontinuity at some x >τ . By construction of c, a discontinuity in the interval [τ, ∞) corresponds to a discontinuity of some candidate function c * induced by a tradeoff subfunction of c 1 or c 2 (see Section 3.2). Without loss of generality, let c * be a subfunction of c 1 . We know that c * has exactly one discontinuity at its minimum driving timeτ * ∈ R >0 . Thus, we have c(x) = c 1 (τ * ) + c 2 (x −τ * ). We distinguish two cases.
Case 1: Ifτ * is not the minimum driving time of c 1 , we know that c 1 is continuous in the neighborhood ofτ * . Since c is decreasing and has a discontinuity at x, there exists an ε > 0 such that , contradicting the fact that c minimizes this term.
Case 2: Ifτ * is in fact the minimum driving time of c 1 , we know that the corresponding driving time x −τ * of c 2 must exceed its minimum driving timeτ 2 , since x >τ =τ * +τ 2 . Hence, c 2 is continuous in the neighborhood of x−τ * and along the lines of the first case we obtain a contradiction, because now c 1 (τ * ) + c 2 (x −τ * − ε) < c(x − ε) must hold for some ε > 0.

A.3. Linking General Consumption Functions in Linear Time
In Section 3.2, we described a naïve link operation, which has quadratic running time in the number of subfunctions of c 1 and c 2 . In what follows, we show how the complexity of the link operation for general consumption functions can be reduced to linear time. We say that a consumption function c with minimum driving timeτ ∈ R >0 is convex if it is convex on the interval [τ, ∞). Note that this holds true for consumption functions of single arcs; see Section 2. Our more sophisticated link operation exploits the fact that consumption functions are always convex, which we now prove formally. Consider the ∆-function ∆ opt of c = link(c 1 , c 2 ), defined as the optimal choice of ∆ in Equation (6). See Figure 14 for an example. (Note that we did not formally prove that the value ∆ is distinct in the general case, but we may as well pick the minimum value ∆ that fulfills Equation (6) to ensure that ∆ opt is well-defined.) Presuming that c 1 and c 2 are convex, the following Lemma A-2 shows that both ∆ opt and the value Proof. We begin by showing that ∆ opt is increasing. Assume for contradiction that this is not the case, i. e., for some value x in the domain of ∆ opt , there are values ε > 0 and δ > 0 such that ∆ = ∆ opt (x) > ∆ opt (x+ε) = ∆−δ. First of all, note that the inequality c 1 (∆)+c 2 (x−∆) ≤ c 1 (∆−δ)+c 2 (x−∆+δ) must hold, since ∆ minimizes this term by definition of ∆ opt . Further, ∆ = ∆ opt (x) is the smallest among all values that minimize the term by definition, so plugging in ∆ − δ < ∆ actually yields a strictly greater result. Analogously, we have c 1 (∆ − δ) + c 2 (x + ε − ∆ + δ) ≤ c 1 (∆) + c 2 (x + ε − ∆), as this term is minimized by ∆ − δ. Therefore, we obtain a contradiction. Here, we exploit the fact that c 2 is convex and decreasing and hence, c 2 (x − ∆) − c 2 (x − ∆ + δ) must be decreasing wrt. x for fixed values ∆ and δ (the gap between two function values with constant difference on the x-axis must decrease if x increases). Regarding∆ opt , monotonicity follows from a very similar argument. As before, assume for contradiction that∆ opt (x) >∆ opt (x+ε) for some ε > 0, so the inequality x−∆ opt (x) > x+ε−∆ opt (x+ ε) holds. We plug in ∆ = ∆ opt (x) and ∆+δ = ∆ opt (x+ε) to obtain δ > ε > 0. As in the first case, we get c 1 (∆)+c 2 (x−∆) ≤ c 1 (∆+δ)+c 2 (x−∆−δ) and c 1 (∆+δ)+c 2 (x+ε−∆−δ) < c 1 (∆)+c 2 (x+ε−∆) by the definition of ∆ opt . Along the lines of the first case, this yields a contradiction.
Finally, we show that ∆ opt and∆ opt are continuous. Since the sum of ∆ opt (x) and x − ∆ opt (x) always equals x, their sum increases by exactly ε for any x + ε with ε > 0 in the neighborhood of x. Hence, a discontinuity in either function would imply that at least one of the two terms must decrease, but we showed before that both ∆ opt (x) and x − ∆ opt (x) increase wrt. x.
Note that in order to show that a consumption function c with minimum driving timeτ and maximum driving timeτ is convex, it is sufficient to consider the interval [τ,τ ), since we already know that c is decreasing and that for all x ∈ [τ, ∞), the value c(x) = c(τ ) is constant. The following Lemma A-3 proves this property and hence, implies that linking two convex consumption functions indeed yields a decreasing and convex function. Consequently, Lemma A-2 applies to all consumption functions of the general form.

A.4. Omitted Details from Section 3.3
The unique extreme point of the difference g i 1 − g j 2 of the subfunctions g i 1 and g j 2 of two consumption functions (c. f. Section 3.3) is given by

B. Omitted Details and Proofs from Section 4
Below, Appendix B.1 mentions technical details from Section 4.2, whereas Appendix B.2 proves Lemma 3 from Section 4.3.

B.1. Omitted Details from Section 4.2
As mentioned in Section 4.2, we store a function∆ opt with each label c to enable path and speed retrieval after the search has terminated. The function∆ opt is given as∆ opt (x) = x − ∆ opt (x), wrt. the link operation that resulted in the label c, and it is a piecewise function with subfunctions that have the general form∆ with nonnegative coefficients λ ∈ R ≥0 and µ ∈ R >0 ; c. f. Equation (8) in Section 3.2.

B.2. Proof of Lemma 3
Lemma 3 proves that the number of settled labels per vertex of the heuristic approach described in Section 4.3 is upper bounded by 1/ε + 1, which implies that the algorithm itself runs in polynomial time.
Lemma 3. The number of settled labels contained in the set L set (v) of each vertex v ∈ V is at most 1/ε + 1 when running TFP with exact improved dominance checks.
Proof. Let k := 1/ε + 1 and assume for contradiction that after running TFP, the label set L set (v) of some vertex v ∈ V contains at least k + 1 consumption functions. We denote these functions by c 1 , . . . , c k+1 and without loss of generality, we assume they are given in increasing order of minimum driving time (and hence, were inserted into the label set in this order). Since exact dominance checks are applied, we know that c 2 yields an improvement over c 1 by at least εM at its minimum driving timeτ 2 ∈ R >0 . Consequently, we have c 2 (τ 2 ) ≤ c 1 (τ 2 ) − εM . We can apply the same argument to any function c i with i ∈ {2, . . . , k + 1} and it follows that This contradicts the fact that c k+1 must have nonnegative consumption for all admissible driving times, since the algorithm ensures that battery constraints are not violated for a consumption function before running any dominance checks with it.
Lemma D-1 shows that the resulting function is in fact an upper bound on c. Upper bounds are also robust towards incremental linking in the sense that the error does not increase if we recompute the upper bound whenever linking several bound functions results in a bound consisting of multiple subfunctions. This is due to the fact that the bounds are uniquely defined by their (minimum) coefficient β and the domain borders of the original functions, which remain unchanged in the upper bound.
During witness search, whenever linking two bound functions results in a function defined by more than one tradeoff subfunction, we compute and store the upper bound instead. Note that we do not even have to link functions explicitly, but simply compute the new coefficient β in a linear scan that simulates the link operation. Proof. Let g 1 , . . . , g k denote the tradeoff subfunctions defining c and without loss of generality, assume that these subfunctions are given in increasing order of their admissible driving times. First, we argue that it is sufficient to prove the lemma for the case k = 2. To show this, we define an operation bound : F × F → F that takes as input two consumption functions, each defined by a single tradeoff subfunction, and computes an upper bound as described above. For i ∈ {1, . . . , k − 1}, consider two consumption functions c i and c i+1 induced by two consecutive tradeoff functions g i and g i+1 (see Section 3.2 for the definition of induced consumption functions). Let their corresponding minimum and maximum driving times beτ i ,τ i =τ i+1 , andτ i+1 . The bound operation computes the consumption function c i,i+1 := bound(c i , c i+1 ) with minimum driving timeτ i , maximum driving timeτ i+1 , and a single tradeoff subfunctionḡ i,i+1 . According to Equations (14)- (16), the coefficients of this tradeoff function depend only on the values β = min{β i , β i+1 }, the driving timesτ i andτ i+1 , as well as the consumption values c i (τ i ) =ḡ i,i+1 (τ i ) and c i+1 (τ i+1 ) =ḡ i,i+1 (τ i+1 ) at the domain borders ofḡ i,i+1 . Linking the consumption function c i,i+1 with another consecutive (induced) consumption function yields a new function defined by the same corresponding values. Consequently, the result bound(. . . bound(. . . bound(c 1 , c 2 ), . . . ), c k ) of iteratively applying the bound operation to the k induced consumption functions of c is the function that is defined by the coefficient β = min i∈{0,...,k} β i (see Equation (14)), the minimum driving timeτ =τ 1 , the maximum driving timeτ =τ k , the maximum consumption c(τ ) = c 1 (τ 1 ), and the minimum consumption c(τ ) = c k (τ k ). This is exactly the functionḡ defined above. Thus, we can constructḡ by iteratively applying the bound operation to consumption functions induced by the tradeoff subfunctions of c. To prove the lemma, we now show that each function constructed by the operation bound is in fact an upper bound on its two input functions. Observe that this implies thatḡ is an upper bound on c within the interval [τ,τ ].
In the remainder of the proof, let c be a consumption function defined by two tradeoff subfunctions g 1 and g 2 , which induce two consumption functions c 1 and c 2 . We prove that the functionḡ : R >0 → R computed by bound(c 1 , c 2 ) yields an upper bound on c on the interval [τ,τ ]. Let the subdomains of g 1 and g 2 be [τ, τ ) and [τ,τ ), respectively. By continuity of c on the interval [τ,τ ] and by continuity of both g 1 and g 2 on R >0 , we know that g 1 (τ ) = g 2 (τ ). To prove the lemma, we use the following three claims.
2. The slopes (i. e., the derivatives) ofḡ and g 1 are equal atτ if and only ifḡ ≡ g 1 ≡ g 2 . Otherwise, the slope ofḡ is greater at this point.
Then,ḡ(τ ) ≥ g 1 (τ ) holds by our first claim,ḡ and g 1 intersect atτ by construction, andḡ(τ + ε) ≥ g 1 (τ + ε) holds for ε > 0 in the neighborhood ofτ by our second claim. This implies thatḡ must be an upper bound on g 1 on the interval [τ, τ ], because the functionsḡ and g 1 can intersect at most twice in this interval unlessḡ ≡ g 1 . This is easy to verify by determining the number of zeros ofḡ − g 1 within the considered interval [τ, τ ]. A similar argument holds forḡ and g 2 on the interval [τ,τ ]. Hence, the lemma follows after proving the three claims made above. We detail the rather technical proofs of these claims below. Assume that the functions g 1 and g 2 are given as g i (x) = α i /(x − β i ) 2 + γ i for all x ∈ R >0 and for i ∈ {1, 2}. For the sake of simplicity and without loss of generality, we presume that min{β 1 , β 2 } = 0. Note that we can always enforce this property by shifting both functions (and their subdomains) along the x-axis. Afterwards, we obtain the same functionḡ on the shifted subdomains. By a similar argument, we presume that γ 1 = 0 holds. Below, we consider the case β 1 = 0 and β 2 ≥ 0. The case β 1 ≥ 0 and β 2 = 0 is analogous. Since γ 2 ∈ R is allowed to become negative, no further case distinction is necessary.
As before, we use γ 2 = α 1 /τ 2 − α 2 /(τ − β 2 ) 2 and the inequality α 1 ≥ α 2 τ 3 /(τ − β 2 ) 3 > 0. For the difference between the derivativesḡ and g 1 atτ , this yields As in the proof of the first claim, we observe that each term in the product of the numerator is nonnegative, while each term in the product of the denominator is positive. Moreover, the numerator is equal to 0 if and only if β 2 = 0 holds. Using Equation (17), it is easy to verify that this implies α 1 = α 2 and γ 2 = 0, which corresponds to the case where the three functionsḡ, g 1 , and g 2 are equivalent.
Finally, we deal with the slopes of g 2 andḡ atτ to prove the third claim. Below, we first replace the values α, α 2 , and γ 2 as in our proof of the second claim. Afterwards, we exploit the fact that (τ 3 − x)(τ − β 2 ) 3 − (τ 3 − x)(τ − β 2 ) 3 decreases with increasing x ∈ R ≥0 , since its derivative wrt. x is (τ − β 2 ) 3 − (τ − β 2 ) 3 < 0. After some further rearrangements, we obtain Again, we end up with products for which all factors are nonnegative (and strictly positive in case of the denominator). As before, the numerator equals 0 if and only ifḡ ≡ g 1 ≡ g 2 . Hence, all three claims hold and the proof is complete.

E. Implementation Details
Our implementation stores graphs as adjacency arrays [24], following the dynamic data structures of Delling [25] for efficient insertion and deletion of shortcuts during the preprocessing routine of CH. All algorithms use k-heaps [24,57] as priority queue, where k = 4 except for unsettled label sets, which use k = 2. Compared to Fibonacci heaps [36], these heaps have a higher worst-case complexity, but are faster on sparse graphs (such as road networks) in practice [23]. Implementation details of our speedup techniques are given below.

A* Search.
When computing the potential function π ϕ (see Section 5.2), the number of breakpoints of lower bounds can become quite large. Therefore, we reduce it as follows (while slightly deteriorating the quality of the bounds). Before applying Graham's scan, we replace consecutive pairs of breakpoints in the piecewise linear function by a single one if they are close to each other, i. e., their difference wrt. driving time or SoC is below a certain threshold ∆ x ∈ R ≥0 or ∆ b ∈ R ≥0 , respectively. Two such points p = (b p , x p ) and q = (b q , x q ) are replaced by r := (min{b p , b q }, min{x p , x q }). Furthermore, if two consecutive segments pq and qr with p = (b p , x p ), q = (b q , x q ), and r = (b r , x r ) have similar slopes s pq ≈ s qr (i. e., the difference |s pq − s qr | is below some threshold ∆ s ∈ R ≥0 ), we replace them by a single segment from (b p , x p ) to (b * , x r ) with slope min{s pq , s qr }, which uniquely defines the value b * ∈ R. Clearly, the modified function remains a lower bound. Moreover, consistency of the potential is maintained, as function values can only decrease and changes in the function are propagated by the search. Thus, all steps in the proof of Lemma C-1 still apply. Graham's scan and the breakpoint reduction step are performed on-the-fly during the merge operation. Moreover, we convert consumption functions c a of all arcs a ∈ A to their corresponding lower bounds ϕ a during preprocessing for faster query times. The thresholds ε to determine lower bound errors, ∆ x and ∆ b for close points, and ∆ s for similar slopes are tuning parameters. Smaller thresholds increase accuracy of bounds, but also slow down the backward search. Therefore, we set above thresholds to 2 δ− log M in our experiments, where δ ∈ N is a constant and M is the battery capacity (assumed to be given in kWh). Hence, bounds are more accurate for higher capacities (where the forward search becomes more expensive). The value of δ is again a tuning parameter.

Contraction Hierarchies.
During preprocessing, we determine the next vertex to be contracted using the measures Edge Difference (ED) and Cost of Queries (CQ) according to Geisberger et al. [41]. To reflect the complexity of SoC functions, we add another term Shortcut Complexity (SC), which is defined as |c + | + k|c − | for the SoC function of a given shortcut candidate, where |c + | and |c − | denote the number of tradeoff subfunctions that define the positive and negative part of a shortcut, respectively, and k ∈ N is a tuning parameter. Using penalized weights for negative parts, we favor earlier contraction of SoC functions without a negative part (we use k = 4 in our experiments). The priority of a vertex (higher priority means higher importance) is then set to 64 ED + CQ + SC. We set the priority of all inactive vertices to ∞.
Further, we employ a settled node limit [41] of 128, which limits the maximum number of queue extractions per witness search for better performance. If multiple shortcut candidates with the same tail vertex u ∈ V are constructed during contraction of a vertex, we save time by running only a single multi-target witness search from u. Finally, to improve performance of the backward searches during a query (BFS and potential computation), we explicitly construct and store their more lightweight search graphs from the input graph (enriched with shortcuts, but storing less complex cost functions) during preprocessing. Figure 15 plots query times against solution qualities for different speed step sizes (visualizing the same figures as in Table 3). Figure 16 shows average running times for the same queries as in  Step  Figure 15.: Trading accuracy for running time when using BSP (Ger-PG, 2 kWh). Both plots show query times and result quality for the same set of queries as in Table 3.  Battery capacity [kWh]

Time [ms]
BSP BSP-A* TFP TFP-d. A*-π d CHAsp-π d CHAsp-πϕ CHAsp-ε-π d Figure 16.: Average running times for different battery capacities. For the same set of queries as in Figure 12, this plot shows the corresponding average running time of 100 random in-range queries.