{'fraction_non_alphanumeric': 0.043458072754182564, 'fraction_numerical': 0.01266741152400714, 'mean_word_length': 4.29202180757122, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electric Vehicle Charging Station Search in Stochastic Environments

Electric vehicles are a central component of future mobility systems as they promise to reduce local noxious and fine dust emissions and CO2 emissions, if fed by clean energy sources. However, the adoption of electric vehicles so far fell short of expectations despite significant governmental incentives. One reason for this slow adoption is the drivers' perceived range anxiety, especially for individually owned vehicles. Here, bad user-experiences, e.g., conventional cars blocking charging stations or inconsistent real-time availability data, manifest the drivers' range anxiety. Against this background, we study stochastic search algorithms, that can be readily deployed in today's navigation systems in order to minimize detours to reach an available charging station. We model such a search as a finite horizon Markov decision process and present a comprehensive framework that considers different problem variants, speed-up techniques, and three solution algorithms: an exact labeling algorithm, a heuristic labeling algorithm, and a rollout algorithm. Extensive numerical studies show that our algorithms significantly decrease the expected time to find a free charging station while increasing the solution quality robustness and the likelihood that a search is successful compared to myopic approaches.


Introduction
All around the world, governments and companies try to foster the adoption of electric vehicles (EVs), which are seen as a central component of future sustainable mobility systems that can play a major role in reducing noise, noxious, and carbon emissions, if fed by clean energy sources. While governments introduced national programs to support individual EV purchase, e.g., by means of ecological bonuses (France), tax exemption (Germany), or subsidies (Netherlands), private companies invested in public charging infrastructure (Bensasson 2019, Fröhlich 2019, and major logistic service providers progressively integrated EVs to realize sustainable transportation (DHL 2019, Amazon 2019). Yet the market uptake for EVs fell short in expectation in most countries. One of the main obstacles remaining, in particular for privately-owned EVs, is the customer's perceived range-anxiety, triggered by missing charging infrastructure and incomplete or non-standardized information about charging services available. Hence, facilitating a reliable EV charging process is crucial to decrease range anxiety in order to invigorate the adoption of EVs (Bonges & Lusk 2016).
Different stakeholders focus on different measures to facilitate a reliable charging experience for users.
At strategic level, municipalities try to facilitate reliable charging options by improving the overall charging infrastructure, e.g., by building new charging stations or increasing the capacity of existing ones. Further, enforcing new pricing schemes that support a higher turnover rate of (re)charged vehicles may help to save additional capacities (Bonges & Lusk 2016). However, infrastructure investments require long planning lead times and pricing schemes can hardly be realized without standardization and agreements between charging station operators. At operational level, some online map services exist and aim to help EV drivers to locate available charging stations. Unfortunately, such services struggle with data inaccuracy, e.g., an incomplete coverage with real-time status data. Even worse, charging stations can be blocked by non-charging vehicles (so-called "ICEing") without the inaccessibility being reflected in status data. An empirical study in Berlin revealed a high correlation between such inaccessibility and parking availability, showing that in areas without available parking it is three times more likely for a charging station to be illegitimately blocked than legitimately used (see Table 1).
Drivers may find a selected and seemingly free station to be blocked when arriving there. Such negative user experiences may require additional detours with an already depleted battery in order to find a suitable charging station and thus may even increase a driver's range anxiety. To this end, a reliable stochastic search algorithm that helps to route drivers to an available charging station as Internal TomTom study in Berlin, 2019. The study shows how often charge points were legitimately used (i.e.,"blocked" and actually used), actually available (correct "free"), or illegitimately blocked by a vehicle. The data relates to the amount of free parking spots found in the immediate vicinity. For most of the stations, real-time availability data was not available.
convenient as possible constitutes a valuable algorithmic component for today's map and navigation services and may help to reduce the drivers' range anxiety.
The goal of this paper is to develop such an algorithmic solution that can be readily deployed as a built-in implementation in today's navigation systems and map services. In the remainder of this section, we first review the related literature, before we detail our contribution and elaborate the organization of this paper.

Literature review
We now survey literature on EV routing with uncertain charging stations before we focus on related search problems with stochastic resources. Sweda et al. (2017) were the first to study probabilistic charging station availability as a multistop shortest path problem under uncertainty. They used a dynamic programming approach on a grid network and considered waiting times at unavailable stations. Kullman et al. (2017) extended this work by modeling an electric vehicle routing problem with uncertain charging at public charging stations, as a Markov decision process (MDP). Given the large action and state spaces, they proposed an approximate stochastic dynamic programming approach with a look-ahead procedure ). Jafari & Boyles (2017) studied a multicriteria stochastic shortest path problem for EVs, in which stochasticity relates to travel time and energy consumption.
Focusing on related stochastic search problems, Arndt et al. (2016) studied a probabilistic routing problem for on-street parking search with parking spots as stochastic resources and incorporated user preferences for stops closer to their destination. Besides showing the problem's NP-completeness, they proposed a brand-and-bound algorithm for small problem spaces. More generally, Guo & Wolfson (2018) studied a probabilistic spatio-temporal resource search problem, in which resources have a general usage cost but the resource seeker is not allowed to wait for an occupied resource to become available again. Contrary to Arndt et al. (2016) in which resources observations are persistent during the search, stations can become available again after a defined time threshold. Guo & Wolfson (2018) proposed a value iteration solution procedure that remains tractable by making a fast recovery assumption at the instance level, which keeps the state space small. Schmoll & Schubert (2018) studied a dynamic resource routing problem under reliable real-time information, which requires fast (re)computations as resources frequently change their occupancy state. Assuming a large problem space, they used real-time dynamic programming to maintain utility values for likely states to determine the best action on-the-fly.
As can be seen, related research on stochastic charging station search is still scarce. Approaches that are specifically tailored to EV charging station search are limited to finding best cost paths between an origin and a destination rather than on an open search with an a-priori unknown destination. Further these approaches lack the consideration of sufficient real-world problem variants and heterogeneous station characteristics. Literature on stochastic resource search focus on an open search concept but lack EV and charging station specific characteristics. Moreover, the algorithmic solutions that have been proposed in this domain remain limited to a pure academic interest on artificial and small instances.
Accordingly, these approaches cannot be embedded into real-time environments.

Contribution
With this work, we close the research gaps outlined above by providing a profound algorithmic framework that covers a wide range of problem variants in the context of stochastic EV charging station search, which can readily be deployed into real-time environments for online optimization purposes.
Specifically, our contribution is fourfold. First, we formalize the problem of stochastic EV charging station search as a new stochastic search problem. We model this problem as a discrete finite horizon MDP and consider two objectives: minimizing i) the time until charging, and ii) the time until completion in case of heterogeneous charging times. Moreover, we present different problem variants in which a driver may wait or not at occupied stations and charging stations might be homogeneous or heterogeneous. Second, we prove the complexity of this problem class and show that it is NP-hard. Third, we develop an algorithmic tool chain that consists of three algorithms: an exact labeling algorithm for which we present a cost decomposition of the Markovian policy to derive effective dominance rules; a heuristic variant of this labeling algorithm which is amenable for real-time application; and a rollout algorithm. Additionally, we present multiple speed-up mechanisms, e.g., reduced action spaces and sharpened dominance relations. Fourth, we provide extensive numerical studies that base on real-world data for the cities of Berlin and San Francisco. Our results show that compared to myopic approaches, our search algorithms decrease the average time spent finding an available station by up to 44%. We further benchmark the performance of our heuristic algorithms against the exact algorithm and show that a combined, case-dependent usage of both allows for effective real-time application.

Organization
The remainder of the paper is as follows. In Section 2, we introduce the stochastic charge pole search (SCPS) problem. Section 3 formalizes our problem as an MDP and consecutively develops the corresponding algorithmic framework. In Section 4, we describe our case study and the experimental design. Section 5 discusses our numerical results. Section 6 concludes this paper and provides an outlook on future research. To keep this paper concise, we shift all proofs to Appendix A.

Problem definition and representation
In this section, we introduce the SCPS problem for which we specify four problem variants, each corresponding to a distinct real-world scenario.

Problem setting
We focus on a routing problem with stochastic charging station availability, where a driver starts at a given location and seeks to find an unoccupied charging station to recharge her vehicle. Her objective is to minimize her total expected cost, which consists of the driving time during the search and additional factors, e.g., the time until the charging is complete, or the time spent walking from the charge point to the actual destination. We consider this search to be spatially and temporally bounded to account for the driver's limited time budget and penalize unsuccessful searches with a high termination cost to model the resulting discomfort.
Formally, we define this problem on a directed and complete graph G = (C, A) that consists of a set of vertices C and a set of arcs (c, c ) ∈ A. Each vertex c ∈ C, except the vertex where the search starts, represents a charging station. The driver starts her search at a designated start vertex c 0 ∈ C at time 0, with a maximum search time of T . We assume that the driver is willing to charge at any vertex in C and consider a limited search radius by restricting C to the stations within a maximum distance S around the start vertex. Driving from c to c takes t c,c ≥ 0 units of time. We model the availability of a station c ∈ C at time t ∈ [0, T ] as a visit-dependent binary random variable a c ∈ {0, 1}, which the driver observes by visiting c. Without prior knowledge, i.e., when visiting a station for the first time, the probability that the station is free (a c = 1) is a constant p c . Table 2 summarizes this notation.
We introduce two time-based penalties: γ c is the time-equivalent usage cost for using pole c if it is available upon arrival; β c denotes the cost for unsuccessfully terminating the search at c in case c is occupied, which may happen if all vertices in C have been unsuccessfully visited, or it is impossible to visit another candidate within the remaining time budget T .
We define a solution to our problem as an ordered sequence of charging station visits C = (c 0 , . . . , c n ).
In practice, a driver that uses this solution starts at c 0 an visits the charging stations in the given sequence up to the first available charging station. Accordingly, our search terminates either successfully at any c ∈ C (with a total cost of driving time until that vertex, plus extra cost γ c ), or unsuccessfully at c n (incurring driving time along C plus cost β cn ). We denote by α the expected cost of solution C.
Problem setting variants: In practice, different system characteristics and user preferences may require to solve different problem variants, e.g., some drivers may want to find a charging station as fast as possible without considering station heterogeneity, while others may want to finish charging as quickly as possible. Such differences can be incorporated in our problem through the generic penalties β c and γ c , and an additional parameter W, which states whether it is permitted to finish the search by waiting at a station W = W 1 or not W = W 0 .
We identify four different problem variants that reflect the most common use cases in practical applications. Table 3 summarizes these variants and states the respective realizations as a triple (β c , γ c , W). In variant ¬W/¬C, the driver is neither allowed to wait at stations nor do stations have (heterogeneous) usage costs. To this end, we consider a constant penalty cost β for an unsuccessful terminated search. While stations remain homogeneous, the driver is allowed to wait in variant W/¬C.
We assume that an expected waiting duration W c is known for every station c, and to be a constant independent of arrival time and earlier observations. In variants ¬W/C and W/C, stations have  heterogeneous charging durations and the driver seeks to minimize her total time to finish charging up to a preferred charge level. While the driver is not allowed to wait in ¬W/C (i.e., the constant penalty cost β is induced in case of failure), this constraint is relaxed in W/C.

Discussion
While our model captures the essential characteristics of the underlying real-world problem, some comments are in order. First, we differentiate charging stations based on the charging duration in problem variants ¬W/C and W/C. The model is however not restrictive and allows charging stations to be heterogeneous with respect to any other criteria, including walking distance to a destination location, price for charging, etc. Yet one needs to convert non time-based usage costs (e.g., charging prices) to time-based costs to ensure search cost compatibility. Second, we assume an occupied station to remain occupied during the search for most problem variants. This seems to be a plausible assumption based on the large difference between typical search times (minutes) and charging durations (hours) in an urban setting. However, we also present a problem variant with recovering probabilities in Section 3.4 and discuss its impact in Section 5. Third, our basic model is agnostic of the energy spent during the search. Here, a similar argument holds: given that a search typically covers only a small distance radius compared to a vehicle's range, we assume this effect to be negligible. Nevertheless, we present a problem variant which considers the energy spent during the search. Our results for this problem variant show no significant impact of the respective energy consumption and thus confirm our assumption.

Methodology
We now introduce an MDP representation for the SCPS problem (Section 3.1) and a corresponding algorithmic framework. As the SCPS problem is NP-hard (see Appendix B) we develop one exact and two heuristic algorithms. We first present these algorithms for our basic problem variants and focus on an exact and a heuristic labeling algorithm in Section 3.2, before we present a rollout algorithm in Section 3.3. We then show how these algorithms must be modified to account for time dependent recovery probability functions (Section 3.4) and search-related energy consumption (Section 3.5).

Markov decision process
To model the SCPS problem as a finite-horizon MDP with multiple decisions over a time-budgeted process, we use the additional notation as summarized in Table 4. x State ∈ S u(x) Action choosen in state x ∈ S w Decision to wait (and charge) at the last visited station C if occupied d(u, u(x)) Immediate cost induced by taking action u in state x ∈ S π Policy V π Value function, based on policy π V * Optimal Value function a Binary variable modeling the availability of the last-visited station δ Binary variable indicating whether the selected station c in x = (C, 0) is already contained in C t(π) Driving time follwing policy π ρ(π, k) Probability that at least one of the first k first stations of policy π is available α(π) Total search cost for following policy π

Model variables and transition functions
We define a state x as a tuple (C, a), where C is an ordered sequence of stations C = (c 0 , ..., c k ) that have already been visited, with the driver being located at c k . Let a be the binary realization of the availability of c k , indicating whether c k is free (a=1) or occupied (a=0). Then, the state space results We denote by u(x) = (c, w) a possible action taken in state x ∈ S to transition to state x , with w ∈ W being a binary that indicates whether to wait at the current station (w = 1) or not (w = 0), and c stating the next station to visit. We note that, either W = W 0 = {0} for problem variants without waiting or W = W 1 = {0, 1} for problem variants that allow to wait at a charging station. Since we assume a station availability status to be persistent during the search, we further assume that for problem variants without waiting, each station should be visited at most once. For problem variants with waiting, given the heterogeneous termination costs β c , we assume that the last station can be visited twice. Actions are only taken at occupied stations: if a = 1, x is a termination state and the search finishes as the driver charges at the found unoccupied station. Let C(C) be the restricted set of charging station vertices, reachable in less than T − k−1 i=0 t c i ,c i+1 . Then, the action space for a state x = (C, 0) results to with δ being the binary variable that indicates whether c is already included in C (δ = 1) or not (δ = 0). If C(C) is empty, no more station can be reached under the given time constraints and we refer to this state as a forced termination state.
We define the function p t (x | x, u) to describe the probability that following action u ∈ U(x) from state x ∈ S would result in state x ∈ S, such that x ∈S p t (x |x, u) = 1 and consider two cases: 1. If the selected action is to wait, then the driver stays at the current station. The transition is deterministic and p t ((C, 1)|(C, 0), (c k , 1)) = 1 for a path C ending in c k .
2. If the selected action is to continue searching, the selected next station c is either available or occupied with respect top c (C). Hence p t ((C , 1)|(C, 0), (c, 0)) =p c (C) and p t ((C , 0)|(C, 0), (c, 0)) = 1 −p c (C), where C is the sequence C extended by c.
We denote by d(x, u) the immediate induced cost for taking action u = (c, w) in state x = (C, a), which depends on the realized availability a at the last station c k and the respective action We note that if a = 1, the driver charges at station c k , and d(x, u) corresponds to the station usage cost γ c k . If a = 0 and the driver waits at the station, the search terminates and cost β c k results. If a = 0 and the search continues at the next station c, the cost results to the travel time t c k ,c .
Finally, we define a policy π as a function mapping a state x ∈ S to an action π(x) ∈ U(x).
Accordingly, π implicitly describes a search path C(π) = (c 0 , ..., c n ) with π(c i , 0) = (c i+1 , 0) ∀i = 0, . . . , n − 1 and state x n = (C(π), a) that terminates the search at vertex c n , either because the driver runs out of time or because she decides to wait with π(x n ) = (c n , 1). We refer to C(π) as C and tõ p c (C) asp c to keep the notation concise.

Cost function
We now analyze the cost function V π (C, a) that describes the expected cost for following a policy π, from a start state (C, a). Then V π (x 0 ), with x 0 = (c 0 , 0) represents the expected search cost for the driver when following policy π from starting at vertex c 0 and the objective is to find a policy π that minimizes V π (x 0 ). Then, the cost function V π (C, a) can be expressed as follows with C = (c 0 , . . . , c k ), u(C, a) = (c k+1 , w), and C = C∪{c k+1 }. For both realizations of a, Equation 3.3 can be simplified: (3.5)

Cost structure variants
For each problem variant as introduced in Section 2.1, the cost functions V π (C, 0), V π (C, 1), and V π (x k ) for a termination state x n = (C, a) can be expressed as follows, with C denoting the extension of C by station c k+1 for non-termination states.

Cost function expansion
We now expand V π to derive an explicit evaluation of the search cost. To simplify notation, we use C [i:j] to denote a sub-sequence (c i , . . . , c j ) of a sequence C = (c 0 , ..., c i , ..., c j , ..., c n ).
Proposition 1. Let C = (c 0 , . . . , c n ). Then, the cost for being in state x k = (C [0:k] , 0), with k < n, following policy π until the termination state x n = (C, a) expands as follows (3.10) Given that policies encode solutions, we can express the solution cost α for a policy, and set α(π) = V π (x 0 ). Following Equation 3.10 and Proposition 1, this yields (3.11) For a more concise notation, letρ(π, k) be the probability that a driver fails in finding at least one free station in C [0:k] , while ρ(π, k) = 1 −ρ(π, k) is the probability that she succeeds in finding at least withp c i (C [0:i−1] ) denoting the likelihood that c i is available after having visited all previous stations from c 0 to c i−1 . Furthermore, let Then, we rewrite Equation 3.11 as α(π) =ρ(π, n)β cn + A(π) . (3.14) We denote by t(π) = n−1 i=0 t c i ,c i+1 the accumulated driving time for all stations in C and note that for a feasible solution, t(π) ≤ T holds.

Dynamic programming based labeling algorithms
To find an optimal solution for the SCPS problem, we develop a dynamic programming based labeling algorithm. Similar to solving multi-criteria constrained shortest path problems, we propagate partial policies in order to find an in-expectation cost-optimal policy. Herein, we use a dominance criterion to withdraw non-promising partial policies early to keep the explored search space as small as possible.
Formally, we associate each partial policy π c with a label L c whose resources depend on the problem variant. For no-waiting variants ¬W/¬C and ¬W/C, a label L c = (t c , A c , ρ c , α c , S c ) consists of the accumulated driving time t c , the partial cost A c (cf. Equation 3.14), the likelihood ρ c to successfully finish the search up to vertex c (cf. Equation 3.12), the total cost α c , and the set of reachable and non-visited poles S c . For waiting variants W/¬C and W/C, we add an additional resource R c that denotes the set of reachable but visited poles.
To describe our labeling algorithm, we denote by L a the set of active labels and by L 0 the initial label that corresponds to our start location. Let F cc (L) be a set of resource extension functions (REFs) that expand a label L whose partial policy ends at vertex c to a label L whose partial policy ends at vertex c . Let α(L) be the cost associated with label L. We define δ + (L) as a function that returns a set of tuples (c, c ) which denotes all feasible physical successor locations c ∈ C for a label L whose partial policy ends at c ∈ C.
Using this notation, Figure 1 shows a pseudo code of our dynamic programming algorithm. We initialize our list of active labels L a and our so far best found solution L * with L 0 (l.1) and start propagating labels until our search terminates when L a is empty (l.2). We then process labels in L a in cost increasing order (l.3). Once a label got selected for propagation, we remove it from L a (l.4) and propagate it considering all of its feasible successors (l.5) using the REFs (l.6). We check whether a newly created label L is dominated by an existing label in L a (l.7). If this is not the case, we remove labels which are dominated by L from L a (l.8) and add L to L a respectively (l.9). Whenever the newly created label is a termination label, i.e., its corresponding state indicates that the driver should wait or there are no feasible successors left, we check if the found label improves the so far best found label and update L * accordingly (l.10&11).
In the remainder of this section, we detail the REFs used to extend a label and the dominance criterion to discard a dominated label.

Resource extension functions:
To extend a label L which corresponds to a partial policy ending at vertex c ∈ C to a new label L which corresponds to a partial policy ending at vertex c ∈ C, we write if isNotDominated(L , L a ) then 8: dominanceCheck(L a , L ) 9: L a ← L a ∪ {L } 10: if ((w == 1) ∨ (δ + (L ) = ∅)) ∧ α(L ) < α(L * ) then 11: L * ← L 12: return L * L ← F cc (L) and use the following set F cc of REFs : Dominance criterion: Equation 3.14 decomposes the non-monotonous cost α(π) into monotonous resources A(π) and ρ(π) that partly define a label. Given the monotonicity of A(π), ρ(π), and t(π), we can then consider two partial policies π 1 and π 2 that end with the same vertex visit c and their associated labels L 1 and L 2 and say that for problem variants ¬W/¬C and ¬W/C, L 1 dominates L 2 (L 1 L 2 ), if the following conditions are true: Here, conditions (3.21)-(3.23)&(3.25) ensure that the cost and duration of π 1 is smaller than the cost of π 2 . Conditions (3.23)&(3.24) check whether all non-visited vertices reachable by π 2 can be reached from π 1 as well. Condition (3.25) checks whether all reachable stations from c for π 1 contribute to decrease α c (π 1 ). This check is necessary for settings where π 1 can be extended with k more stations than π 2 to detect corner cases in which these k additional stations may increase α(π 1 ) to a value larger than α(π 2 ).
For problem variants W/¬C and W/C, we recall that a search can terminate before the whole time budget is spent and that a station terminating the search might have been visited earlier. Accordingly we drop condition 3.25 and introduce condition 3.26 that checks whether all visited vertices reachable by π 2 can be reached from π 1 as well. Recalling that the label reads L c = (t c , A c , ρ c , α c , S c , R c ), we then say that for problem variants W/¬C and W/C, L 1 L 2 if conditions (3.21)-(3.24) are true and holds. While Algorithm 1 solves the SCPS problem optimally with the dominance criterion above, one may consider to drop some of the dominance conditions to obtain a heuristic dominance criterion that withdraws more labels at the price of losing optimality. In the remainder of this paper, we study a heuristic labeling algorithm, where we preserve conditions (3.21)&(3.22) to obtain cost dominance but neglect conditions (3.23),(3.24),(3.25) and (3.26). In this context, the label definition simplifies to L c = (A c , ρ c , α c ). We provide evidences to the selection of this heuristic dominance criterion in Section 5.2.

Rollout algorithm
In this section, we introduce a rollout algorithm, built as a forward dynamic programming procedure.
We identify the best action at a given state as the one that yields minimal approximated cost. Here, the core idea of the cost approximation is to greedily expand the current policy from each candidate action until a defined horizon to obtain an associated cost value via backpropagation. We first detail the main procedure for ¬W/¬C and ¬W/C variants before we show how to account for additional waiting decisions in problem variants W/¬C and W/C. We apply a one-step decision rollout strategy (cf. Goodson et al. 2017) whose complexity equals a post-state decision rule as the approximation reduces to w = 0 decisions.
Let k be the index of the k th decision epoch and let x k be the non terminated epoch's state with .., c k ), 0) and c k being the last station visited in epoch k. Let x k+1 = (C , 0) be the state in the (k + 1) th epoch that results from action u = (c, 0) at epoch k, with u ∈ U(x k ). Using this additional notation, Figure 2 details the pseudo-code of our algorithm.
We initialize the sequence of station visits C with the start vertex (l.2) and expand the sequence until the time horizon is reached (l.3). From the current state x k , we seek to determine the next best action (c * , 0) and initialize the variables that encode it (l.4). For all possible succeeding states x k+1 (l.5), we use a heuristic policyπ c that bases on a greedy procedure to propagate state x k+1 up to a forced terminated state x k+K , with a look-ahead of K epoch extensions. For all propagated states x l with l ∈ [k + 1, K − 1], the greedy procedure chooses the action (c l+1 , 0), with station c l+1 being selected C ← refinePolicy(C) 14: return C based on availability probabilities of vertices reachable from c l and the driving time to each of these vertices. We then use these anticipated states to evaluate the expected value of the policy-specific cost Vπ c (x k+1 ) for the candidate state x k+1 . We define Greedycost(x k+1 ) as the function that carries out the greedy expansion from x k+1 = (C , 0) and returns Vπ c (x k+1 ) (l.6&l.7). Repeating the greedy procedure and cost evaluation for all possible next actions u = (c, 0) ∈ U(x k ) allows us to find the action u = (c, 0) that minimizes the cost to transition from state x k to state x k+1 (l.8) that we define as Q(x k , c, x k+1 ) (cf. Equation 3.29). Eventually, the selected action u = (c * , 0) yields minimal cost Q(x k , c * , x k+1 ) (l.10-l.12). Only for variants W/¬C and W/C, refinePolicy(C) bases on the resulting visits sequence C to introduce the waiting decision at the best stage 1 ≤ k ≤ n, with n being the total amount of station visits in C (l.13).
In this setting, the reduced action space for problem variants ¬W/¬C and ¬W/C reads ( 3.27) and we calculate the cost for being in state x k+1 , Vπ c (x k+1 ), based on the greedy policyπ c as For problem variants W/¬C&W/C, we first calculate the no-waiting case and compute policy π with the associated ordered sequence of charging stations C = (c 0 , ..., c n ). We denote π as an intermediate policy and introduce π S representing the final search policy. Then, refinePolicy(C) calculates π S using the intermediate policy π while permitting wait-actions. For each intermediary charging station c k at the k th decision epoch, c k ∈ (c 0 , ..., c n ), k = n, π provides a sub sequence of charging stations to visit until the end of the search (c k+1 , ..., c n ) and thus the policy-specific cost value V π (x k ), associated to policy π and state x k = ((c 0 , ..., c k ), 0). We aim to quantify for each station c k whether the termination cost β c k is actually lower that the expected cost for continuing the search V π (x k ). If this is the case, the optimal decision is to wait and we refine π into π S with π S (x k ) = (c k , w = 1) and We define π S as and π S (x k ) = (c k , 1), then state x k terminates the search, as the driver will wait at c k if c k is not immediately available. In this case π S encodes the solution (c 0 , ..., c k ).

Time-dependent probability recovery function
In the basic setting of the SCPS problem, we assume that a station does not change its availability during the search's time horizon and restrict the amount of visits to each individual station. We now relax these assumptions and show which modifications are necessary to consider a recovery function r c that allows the availability of an occupied station c to recover over time. In this case, we definep c (C) as followsp where ∆ c = k−1 j= t c j ,c j+1 + t c k ,c , with l denoting the position of the last visit to c in C. We still consider a charging station c to be initially available at a probability of p c . When it is blocked at the driver's arrival time t c = −1 j=0 t c j ,c j+1 , it may become available over time, i.e., the availability probability of c recovers according to r c (∆ c ), which denotes a station's availability probability for an arbitrary point in time t = ∆ c + t c that remains after the first visit and before the end of the time horizon. We specify r c (∆) for any ∆ based on Schmoll & Schubert (2018) as Here, 1 λc and 1 µc denote the average time station c remains available, respectively occupied and remain constant over the search's time horizon. We can then express p c as a function of λ c and µ c , with p c = µc λc+µc (cf. Jossé et al. 2015) and simplify Equation 3.32 to We now assume that stations can be visited as many times as needed. Then, the action space for a state x = (C, 0) is slightly modified as follows, In the remainder of this section, we outline the changes that are necessary to adapt the labeling algorithm to such a setting. Apart from the new availability probability definitionp c (C), no modifications are necessary for the rollout algorithm.
Modifications for the Labeling Algorithm: The MDP definition remains unchanged, because we can determine the arrival time at a station based on the arrival time at the preceding station extended by the driving time that remains deterministic. Thus ∆ c andp c (C) can be calculated from C without any further information. However, we need to modify the initial dominance criterion (3.21)-(3.25)(resp. (3.21)-(3.26)) as the optimal solution may now contain multiple visits to any charging station, independent on the problem variant.
We consider two partial policies π 1 and π 2 that end with the same visit at vertex c and their associated labels L 1 and L 2 . We recall that R c (π 1 ) denotes the set of charging stations which have been already visited by π 1 . We then say L 1 L 2 , if (3.21)-(3.26) hold and holds. Here, Equation 3.34 accounts for different probabilities of charging stations that have already been visited in both policies at different points in time by considering the biggest possible difference of probability values for c between both paths. We note that we leave the heuristic variant of the dominance criterion unchanged in this context.

Integrating search related energy consumption
In this setting, we assume that the vehicle starts its search with an initial state of charge (SoC) b 0 , which reduces over the course of the search depending on the driven distances. Then, longer driving distances result in higher energy consumption and an additional trade-off results between visits to far-distanced stations with a high availability probability and visits to near-distanced stations with medium to low availability probabilities. Here, a higher availability probability may be counterbalanced by traveling longer distances, which increases the energy that must be recharged or limits future visits to potential stations accordingly.
To account for this setting, we introduce additional notation and denote the energy consumed when traversing arc (c, c ) as k c,c and a vehicle's SoC after having visited all stations in C = (c 0 , ..., c k ) as b. In this setting, we can only transition from state with b denoting the minimum feasible state of charge a vehicle must keep. Analogously to common monotonicity assumptions in related settings (cf. Sweda et al. 2017), we assume that and state the necessary MDP modifications which hold as follows (3.37) First, we include a vehicle's SoC into the state space (Equation 3.35). Second, we modify the action space such that it depends onC(C), which denotes a restricted set of charging station vertices that are reachable from state (C, a) in less than T −t(π) time, with k c,c +b ≤ b and t(π) = n k=0 t c k ,c k+1 (Equation 3.36). The policy-specific cost function results straightforwardly from the modified action and state spaces (Equation 3.37). To account for this new setting in our algorithms, the following modifications are necessary.
Variants ¬W/¬C and W/¬C: The rollout algorithm requires no modification. For the labeling algorithms, sets S c and R c now denote vertices that are reachable within the remaining time budget and within the remaining energy budget b c . Accordingly, we add to the existing dominance conditions  Variants ¬W/C and W/C: Here, the charging duration L c at station c now depends on the amount of time needed to recharge from b 0 up to the maximum SoC b. We introduce δL c to denote the additional charging time due to the battery depletion δb c during the search, where δb c = b 0 − b c = l=i l=0 k c l ,c l+1 . We then account in each algorithm for the charging duration as . (3.39)

Computational complexity improvements
In this section, we proof additional characteristics that allow to improve the computational complexity of certain problem variants. We first introduce three action space reductions, before we focus on a sharpened dominance relation for both the exact and the heuristic labeling algorithm.

Action space reductions
In the following, we discuss some action space reductions, which we summarize in Table 5. We refer to the initially defined action space as complete. In addition, we create the following reduced action spaces. 3.6.2. Sharpened dominance relation for the dynamic programming algorithms In the following, we aim to sharpen the dominance relation, i.e., we aim to discard labels faster without dropping optimality. In our initial problem setting, we assume that both travel times and charging station availability probabilities are unbounded. We now account for a bound on each of these values that still reflects a real-world application. Specifically, we assume that i) the EV driver must travel at least a certain amount of time between two charging stations and ii) that one can never be entirely sure that a charging station is available p c ≤ p < 1, ∀c ∈ C.
These bounds can be computed during a preprocessing step and allow for a sharpened dominance relation for the dynamic programming based labeling algorithms without invalidating our generic model. We note that charging and waiting times at a station are bounded as well. Let W min and L min (respectively W max and L max ) be the minimal (respectively maximal) waiting and charging times.

Design of experiments
To benchmark our algorithms, we develop real-world instances that allow for extensive simulation experiments. We consider three different spatial patterns (see Figure 3), based on the west side of San Francisco, USA (SF-1 ), the city center of Berlin, Germany (BER-1 ), and the financial district of San Francisco (SF-2 ). Here, we account for free-flow speeds to calculate travel times t c,c that denote the time-shortest path between two stations c and c . As our search algorithms appear to be rather insensitive to the search's starting point, we randomly choose one starting point for each pattern and use this starting point in every instance that builds on the respective pattern.
Besides the significant sensitivity to the instance size given by S and T , we found during preliminary studies that our search algorithms are sensitive to two general instance characteristics: the search area's charging station density and the charging station availability probability. Accordingly, we use the patterns described above to create a set of instances that covers a broad parameter range for these characteristics. Figure 4 shows the amount of charging stations for each pattern, depending on the search radius S around each patterns starting point.
Our TomTom internal availability study in the city of Berlin (cf. Table 1) shows on average high charging stations availability in areas with a large number of on-street parking spots. In areas with few available parking spots, drivers often use available stations as free parking spots, thus blocking access to EV drivers. In the study, the sole parking availability factor largely impacts the station availability. To reflect these amplitudes, we introduce three availability settings, drawing probabilities p c for each charging station c from a β−distribution, which is centered on an expected availability of   scenario. The avg-60% and high-90% settings represent areas with average to high parking availability.
The low-15% setting depicts a fictionous extreme case scenario prospectively corresponding to stricter parking policies and allows us to evaluate our algorithm behavior in such an environment.
For problem variants W/¬C and W/C, we consider a waiting duration W c for each station c, which is uniformly distributed in [3, 15, 60, 120] minutes. For problem variants ¬W/C and W/C, we consider a charging duration L c , uniformly distributed in [30, 60, 120] minutes, to account for heterogeneous stations.
With this setup, we create a total of 9 scenarios by combining each area (BER-1 , SF-1 , SF-2 ) with each availability probability distributions (low-15% , avg-60% , high-90% ). We then consider a small and a large-size instance for each combination, defined by the search's time budget and search radius
For the sparse area SF-1 , since an 800 meters-area only comprises 2 stations we consider slightly wider search spaces (10/2000), respectively (15/2600) for SF-1 , which derives a set of 18 instances per variant, 72 in total. We denote for each area, the smaller-size instance by area/1 (e.g., SF-1/1) and the larger-size instance by area/2.
For our studies, we set the penalty cost to β = 120 minutes for ¬W/¬C and to β = 200 minutes for ¬W/C and refer to Section 5.2 for a discussion on its selection and sensitivity. To evaluate our algorithms, we conduct N = 1000 simulation runs, each with a different station availability realization, drawn from the respective probability distribution. For all average values reported in our results discussion, we applied two-tailed Wilcoxon signed-rank tests to verify its significance.

Results
In the following, we discuss our results. We first detail the performance of our algorithms with respect to several quality metrics as well as their computational effort (Section 5.1), before we discuss results of additional sensitivity analysis and modeling assumptions (Section 5.2). We implemented the proposed algorithms single-threaded with Python 3.6.9, using PyPy 7.3.0 with GCC 7.3.1 and performed all experiments on a Virtual Machine on a Hypervisor, with 19 cores and 60 GB of RAM, running Ubuntu 18.04.3 LTS. For the following discussion we note that all algorithms can be implemented more efficiently, e.g., in a C++ environment. However, as this work stems from an industry project, we are not allowed to disclose the C++ related computational times for confidentiality reasons and use the Python implementation for a fair discussion of the algorithms' complexity.

Performance analysis
We discuss the algorithms' performance from two perspectives. First, we analyze the practical benefit of our algorithms and compare their solution quality based on their performance in our simulation environment. We then relate these results to a technically focused discussion on their computational complexity.
During algorithm testing and development, we noticed that both the labeling heuristic and the rollout algorithm are sensitive to the used action space (cf. Section 3.6) in between different problem variants. For the following discussion, we report results for the respective best-performing reduced action space and refer to Appendix D for a detailed discussion of this impact.

Applicability
In practice, more than one key performance indicator exists to evaluate the performance of our algorithms. In the following we refer to α as the realized search cost, which corresponds to the realized driving time needed to find a station for variants ¬W/¬C and W/¬C and to complete charging for ¬W/C and W/C. For W variants, α includes the waiting time whenever the last visited station is occupied. Specifically, we analyze for each algorithm its average search cost deviation which we calculate as with α i being the realized search cost of the respective algorithm for simulation run c and α * i being the best realized search cost for simulation run c out of all algorithms; success rate that results straightforwardly from the number of simulations, for which the respective search strategy finally led to a free charging station; maximum search cost deviation that results to ∆ α max = α max / α max * with α max = max α i being the maximum cost out of all simulation runs for a single algorithm and with α max * being the maximum cost over all simulation runs and all algorithms.
As the exact labeling algorithm (LE) cannot solve large instances of our problem to optimality, we limit the comparison to our heuristic algorithms (LH and RO) and add a naive and a greedy solution as myopic benchmarks. The greedy search (G) creates a sequence of charging station visits based on a greedy cost combining the travel time weighted by a station's availability probability and a penalty cost weighted by its occupancy probability. The naive search (N) mimics a driver without assistance by selecting the closest non-visited station with highest availability probability.
Tables 7-9 detail the average search cost deviations between all algorithms (Table 7), the success rate of each algorithm (Table 8), and the maximum search cost deviation between all algorithms (Table 9), based on 1.000 simulation runs for each problem variant and scenario to avoid a statistical bias. Table 7 shows that the RO and the LH algorithm significantly outperform the naive and the greedy benchmark in terms of average search costs. Here, an average deviation of zero implies that an algorithm always found the best search strategy, while increasing deviations indicate that the best search strategy was not found by the respective algorithm. For scenarios with high and average charging station availability, RO and LH show a similar performance. For scenarios with a low charging station availability, we observe higher differences that result in significant deviations for specific problem   ¬W/¬C, scenario SF − 1/1). Table 9 shows the deviation between the maximum search cost of all algorithms for W problem variants, which are guaranteed to be feasible at the price of high waiting cost. Here, we observe a trade-off between the average search cost and the maximum search cost. As one can see, the LH and the RO algorithm show a similar performance for scenarios with medium and high charging station availability, but the LH algorithm performs better in the SF-1/1 scenario with low charging station availability. Across all scenarios the greedy algorithm sometimes yields the least maximum cost. However, this improvement stems from significantly increased search costs (cf. Table 7).   For ¬W problem variants and scenarios with medium or high charging station availability, all algorithms yield a 1.00% success rate but the advanced algorithms yield significantly lower cost, with LH being the best-performing algorithm. In these cases the greedy algorithm performs close to the advanced algorithms for problem variant ¬W/C. For scenarios with low charging station availability the LH algorithm yields a better success rate but the RO algorithm yields lower search cost. For W problem variants, we observe that the advanced algorithms outperform the myopic algorithms significantly, independent of the charging station availability. This performance difference is higher for problem variant W/¬C. For problem variant W/C, G performs close to the advanced algorithms for scenarios with medium or high charging station availability.
Concluding, our results show that the developed search algorithms can significantly improve the search quality across all problem variants and scenarios. Compared to G, the advanced algorithms decrease the search cost by 21 % in average and up to 44% for areas with a scarce number of charging stations and low charging station availability. For ¬W variants, the failure rate decreases by 30 % with low charging station availability. Moreover, advanced algorithms allow to reduce search times by 5 (W/¬C) to 31 (W/C) minutes compared to myopic approaches for W problem variants with low station availability. This reduction potential decreases for scenarios with average to high station availability. However, advanced algorithms appear to be more robust in these cases and lower the worst search costs by 20 % (W/¬C) and 12 % (W/C). Comparing the LH and RO algorithm among each other, LH prolongs the search compared to RO but obtains a significantly higher success rate for ¬W variants. For scenarios avg-60% and high-90% RO and LH show a similar performance.

Computational tractability
To compare the performance of the LH and RO algorithm against the exact labeling algorithm, we derive a set of 8064 test instances by varying T ∈ [5, 10, 15, 20] minutes, S ∈ [800, 1.000, ..., 3400] meters, and the availability distribution in {low-15%, avg-60%, SF-2} for all problem variants and each search area. Here, we use a large time limit of 15,000 seconds to obtain a maximum set of solutions that are eligible for our comparison, i.e., solutions for instances that could be solved with all three algorithms. Table 10 shows the rate of instances solved within this time limit (n), and the average computational time of the successful runs (t) for both heuristic algorithms (RO, LH) and for the exact labeling algorithm (LE), as well as the average ratio (g α ) between each heuristic result and the optimal solution.
As can be seen, our LH algorithm provides optimal or close to optimal solutions. The RO algorithm shows a similar solution quality for problem variants of type W but shows a significantly worse solution quality for problem variants of type ¬W with a low charging station availability. While the RO algorithm succeeds in solving all instances with computational times of a few seconds, the LH algorithm improves upon its exact counterpart with respect to the number of instances solved and computational times but can neither solve all instances nor preserve computational times at the order of seconds. Figure 6 illustrates the solution quality deviations between RO and LH for a representative subset of scenarios (BER-1 , T = 10) with low and medium charging station availability and varying search radii S ∈ [800, 1.000, ..., 3400]. As can be seen, the observed differences are sensitive to the problem variant. Significant differences occur for ¬W problem variants with low charging station availability and varying search radii, while the algorithms perform similarly on all other problem variants. These high deviations result from penalty costs for unsuccessful searches, which can only result for ¬W problem variants and are more likely to occur at low charging station availability. Abbreviations hold as follows: g α -averaged optimality gap over all tested instances with g α = α heur α opt , t -averaged computational time[s],n -rate of instances that can be computed in less than 15000 seconds. We note that an average g α of 1.00 indicates that an algorithm (almost) always finds the optimal solution. If it always finds the optimal solution we highlight the respective g α in bold font, whereas we leave it in normal font if some solutions remain heuristic but are not reflected in the value of g α due to rounding. Figure 7 shows the computational time of LH, RO, and LE for a representative subset of instances (BER-1 , T = 20, high-90% ) for different search radii. As can be seen, RO and LH remain equally efficient for search radii up to 1400m but the computational time of LH increases exponentially for bigger search radii. Synthesizing Table 10 and Figures 6&7, we observe a trade-off between LH and RO: while RO yields robust computational times at the price of varying solution quality, LH yields a robust solution quality at the price of exponentially increasing computational time.
From a practitioner's perspective, computational times of a few milliseconds are imperative to deploy a search algorithm in practice, e.g., embedded into a navigation application. Here, one could resolve the trade-off between RO and LH in two different ways. On the one hand, one could apply both algorithms selectively, using LH for tractable problem sizes and RO for larger problem sizes. On the other hand, one could always apply LH and terminate its search after a given time limit. To analyze which strategy appears to be more promising in our case, we compare the performance of RO and LH against each other, limiting the computation time of LH to 1 second, which equals a sufficiently small computational time when using an efficient implementation. Each subplot shows ∆α = α LH − α RO (with α LH , resp. α RO being the solution cost for LH, resp. RO) as a function of T and S, where we limit LH computational times to 1 second. Each subplot corresponds to one of the 9 scenarios resulting from the combination of each area (SF-1 , BER-1 , SF-2 ) with each availability distribution (low-15% , avg-60% , high-90% ). Over all subplots, availability increases from left to right and station density increases from top to bottom.

Extended analysis
In the following, we analyze the sensitivity of our algorithms towards parameter and design decisions.
We first analyze the algorithms' sensitivity towards the penalty costβ for ¬W problem variants, before we study the impact of a time dependent recovery function, substantiate the heuristic dominance decisions and study the impact of additional charging times due to battery depletion during the search.

Termination penalty
We limit the discussion of termination penalty sensitivities to low charging station availability instances as these appear to be the most sensitive. Apparently, analyzing the search cost α does not allow for a meaningful interpretation ofβ sensitivities as α naturally increases with increasingβ. To circumvent this issue we decompose α into user-relevant metrics (see Appendix C), and analyze the computational time, the expected search time t s and the probabilityρ that the search unsuccessfully terminates, averaged over six low-15% -instances for problem variant ¬W/¬C (9a) and ¬W/C (9b) in Figure 9. As can be seen, we observe a goal conflict between the expected search time and the search's success rate: with increasing β we obtain better success rates at the price of higher expected search times.
We note that one must choose β ≥ 30 minutes for problem variant ¬W/C, as otherwise the cost for not visiting any station is lower than the cost for visiting at least one station (i.e., minimal charging time). We observe that the β values which are necessary to obtain the best possible success rate are significantly higher (120 minutes respectively 200 minutes for ¬W/¬C and ¬W/C) than this lower bound.
While the computational times of LE significantly decrease with increasing β, the computational times of LH and RO remain insensitive to changes in β for low charging station availability scenarios.
However, additional analysis show that LH computational times increase with increasing β in high charging station availability scenarios, such that varying β for risk-adverse searches in a real-world implementation (i.e., choosing a higher β) should be chosen with respect to the computational overhead for LH.

Time-dependent recovery function
To analyze the impact of considering time-dependent recovering probabilities, we compare the objective values obtained in our main model (without recovering) with the model introduced in Section 3.4 (with recovering). In the former we compute the search path with a persistent charging station occupancy but calculate the objective value α considering time-dependent recovery for a fair comparison of both models. Based on preliminary studies, we set the average occupancy time to two hours, 1 µ = 120 minutes. Table 11 compares the value of α for both models. As can bee seen no significant difference exists between the initial model and the updated model at the exception of smaller areas with large time budget, particularly in case of low station availability. This is the only particular case where visiting stations multiple times within the time budget might be worth an extra detour, as one could expect at least one of the observed stations to be freed.

Relaxed dominance criteria
To design our labeling heuristic, we relaxed the dominance check of LE as described in Section 3.
In the following, we substantiate the design decision from Section 3. Here, we identify each possible variant of the dominance check with a boolean quintuple that signifies whether an equation (3.21, 3.22, 3.23, 3.24, 3.25 resp. 3.26) is active (= 1) or not (=0) in the respective dominance check, e.g., quintuple (1,0,1,0,0) identifies the dominance check variant in which only equations 3.21 and 3.23 are active. Figure 10 shows the trade-off between the optimality gap and the computational times for all dominance criterion and problem variants. As can be seen, the (heuristic) dominance criterion as chosen in Section 3 -(1,1,0,0,0) -yields the lowest computational times possible to achieve the best possible solution quality among all heuristic dominance criteria for problem variants ¬W/¬C, W/¬C, and ¬W/C. For W/C, (1,1,0,0,1) yields the best trade-off by slightly decreasing computational times obtained with (1,1,0,0,0) but selecting (1,1,0,0,0) allows the best possible generic implementation for LH.

Battery depletion
For the following experiments, we consider a flat topography and base additional charging time calculations due to battery depletion on the technical characteristics of a Renault Zoe (battery: Z.E.50 [52 kWh], engine: R110/R135, (Renault 2020)). Furthermore, we assume that the battery level remains at minimum at 20% of its maximum capacity to account for safety considerations or driver anxiety.
We approximate the extra charging time based on a linear charge curve since the additional depletion adds only a limited amount to a (partial) recharge and remains in the linear part of the overall charge curve. To calculate energy consumption, we assume a constant speed of v = 50km/h, which remains a worst case estimate in an urban context. We define L c as the expected time to charge the battery from its initial state to a full state at station c. Then, the depleted amount of energy after t (additional) minutes of driving results to δb c (v) = t Ta(v) × 52. In this setting, our results show that the impact of additional charging times due to an extended search and the resulting battery depletion is very limited in urban areas. For LH, we observe an average objective value change of 0.02%, which amounts to a maximum deviation of 1.5% for single instances.
The RO algorithm appears to be even less sensitive to the changed objective and shows no significant differences.

Conclusion and outlook
In this paper, we studied charging station search algorithms for stochastic environments, motivated by real-world applications in today's navigation system applications. We introduced the underlying problem as a finite horizon MDP that covers several real-world problem variants. In this setting, we aim to find cost minimal search paths for an individual driver. We developed three solution algorithms: an exact labeling algorithm as well as a heuristic labeling algorithm and a rollout heuristic. We benchmarked these algorithms using an extensive real-world case study with instances for the cities of San Francisco and Berlin. Our results show that the heuristic algorithms allow for a significant speedup compared to the exact algorithm at a price of a reasonable performance loss. Moreover, we show that our algorithms significantly improve a driver's success to find a free charging station compared to myopic and greedy search approaches. This is in particular the case if the number of free charging stations in the search area is scarce. In this case, our algorithms reduce a driver's search time by up to 44%.
In future work, we aim to leverage this work to study the impact of coordination and information sharing between multiple drivers. By so doing, we can study the impact of additional coordination that reduces the amount of uncertainty in the system. Moreover, studying the charging station search problem from a system perspective with a perfect information setting may yield an interesting upper bound that allows for an improved assessment of the solution quality of our algorithms in a stochastic setting.
Proof of Proposition 1. Let π be the policy associated to sequence C = (c 0 , c 1 , ..., c n ). We consider state x k = (c k , 0) with charging station c k not being available and c k = c n . Thus, for i < n, i j=i (1 − p c j ) = 1. For i = n, i j=i (1 − p c j ) = 1 − p cn . We then introduce F as follows We notice that F π (x n ) = (1 − p cn )β cn + p cn γ cn = V π (x n ).
We now show that F fulfills the recursive definition of the policy specific cost function and by recursion that F = V . From state x k = (C [0,k] , 0), we let the cost for being in next state x k+1 = (C [0,k+1] , 0) and seek to express F π (x k ) as a function of F π (x k+1 ) to fulfill the recursive definition 3.5.
Accordingly, F fulfills the recursive definition 3.5 for w = 0, which concludes the proof.
Proof of Proposition 2. We consider two simple search sequences (c) and (c, c ) extended with the same visit sequence C = (c 0 , ..., c n ). Let C = (c) • C, resp. C = (c, c ) • C, associated to policies π , resp. π . Let t c,c 0 = t c,c + t c ,c 0 and let c be a direct neighbor of c (i.e., there is no station d such that t c,c = t c,d + t d,c ) and let c 0 not be a direct neighbor.
Case 2 (c ∈ C): In this case, the valuation of any unexplored station after c 0 depends on preceding visits in the respective sequence and we obtain the following cost expansions, visiting c at position k: For path C , c is visited for the first time such thatp c k = p c and we get For path C , c is visited for the second time such thatp c k = 0 and we get and consequently In both cases, π is preferred over π (thus C over C ), such that candidate stations can be restricted to neighbor stations only, which concludes the proof.

Appendix B Problem Complexity
Proposition 3. The SCPS problem is NP-hard, even with metric travel times.
Proof. We show hardness through reduction from the traveling salesman problem (TSP) with metric and integer travel times. The decision problem variant of the TSP can be defined as follows: we consider a set C of n sites and travel times t c,c ∈ N between these. Travel times are bounded 1 ≤ t c,c ≤ ∆ for all c, c ∈ C. We are asked for a tour (i.e., a Hamiltonian path) c 0 , . . . , c n = c 0 whose length satisfies n−1 i=0 t c i ,c i+1 ≤ θ for a given θ. Given that travel times are integer, we can assume w.l.o.g. θ ∈ N. Further, we assume that triangle inequality holds, i.e., t c,c + t c ,c ≥ t c,c for all c, c , c ∈ C. We note that hardness for this restricted metric case implies hardness for the generic case as well.
Step 1: We construct an instance for the station search from the TSP instance as follows: We select an arbitrary vertex s ∈ C and designate it as start vertex c 0 := s for the search. We then create a duplicate s in the same location (t s,s = 0), which serves as the termination vertex. Let q be an arbitrary value satisfying We parameterize the search as follows: All vertices c have an availability of p i := 1 − q without recovery. There is no penalty for successful charging (γ i := 0 ∀i). For unsuccessfully terminating at s , the penalty is For all other vertices c = s the penalty is Now for any search path C = (c 0 , . . . , c k ) that does not visit any vertex multiple times, it holds that

(B.4)
Step 2: We now claim that the TSP instance possesses a solution with cost at most θ ∈ N if, and only if, the station search admits a search path C with α(C) ≤ qθ + q n+1 β s . This is done by transforming solutions between the two problems and carefully mapping their objective values.
For the first direction, we assume that a TSP tour is given. We convert it to a search path C = (c 0 , . . . , c n ) by cutting at s such that c 0 = s and c n = s . Then Vice versa, we assume that we are given a search path P with α(P ) ≤ qθ + q n+1 β s . Then for an optimal search path C = (c 0 , . . . , c k ) it holds that α(P ) ≤ α(C) ≤ qθ + q n+1 β s . Given metric travel times and no recovery, we can assume that C does not visit any vertex more than once. We construct a tour through the following observations: 1. C visits s : Assume it does not. Let C be C extended by ending at s . Then using in (*) that β c k > qβ s + ∆ by (B.3). This contradicts the optimality of C.
2. C visits s last: Assume it does not. We obtain C from C by moving s to the end. Then it holds that Given that q n+1 β c k > β s + n∆ (B.3), it holds that α(C) > α(C ), contradicting optimality.
4. It is clear now that C corresponds to a TSP tour, by identifying s with s . It holds that Assume the tour would violate the threshold, i.e., n−1 i=0 t c i ,c i+1 > θ. Given integrality, the length then is at least θ + 1. It follows that This contradicts (B.8), thereby proving that n−1 i=0 t c i ,c i+1 ≤ θ.

Appendix C Search cost decomposition
In this section, we show by proving Proposition 4 that the total search cost α(π) as derived in Section 2 can be expressed as a function t(π), ρ(π) and an additional quantity t s (π) representing the expected time to find a station assuming at least one station is available among C, the visits sequence associated with π. We derive the quantity t s (π) based on the work of Arndt et al. (2016) that describes t s (π) as a sum of the partial accumulated driving time to a station (parking spot in their work) weighted by the probability the driver charges exactly at this station.
We now note thatρ k p k represents the probability of station k being available when visited, given that no station in i ∈ (0, .., k − 1) was available.
We then introduce the quantity B n = n−1 k=0ρ k−1 p k (t k + γ k ) such that t s n · ρ n = B n and note that to prove C.1, it is sufficient to show A n = B n +ρ n · t n , (C.2) which follows by recursion: Step 1: For n = 0, C.1 holds : A 0 = t 0,1 and B 0 = p 1 t 0,1 +p 1 t 0,1 = A 0 .

Appendix D Reduced action spaces results
Preliminary results show that using the complete action space for problem variants ¬W/¬C and W/¬C with LH is computationally too heavy to be of any practical use, such that we restrict results to the other restricted action spaces. As can be seen, while direct only slightly helps saving computational times on average for ¬W/¬C and LH, it allows to solve 6% more instances within the allocated time for W/¬C. For RO, in W/¬C, results show a 96% decrease of the computational time. For ¬W/C, restricting the next station visits from the current location to the ones accessible in less than 5 minutes allows to save 36% of the computational times compared to complete. Accordingly, Table 13 shows the most appropriate action space for each heuristic and problem variant. For all problem variants, the table shows for each heuristic (LH, RO) the graph setting to that provides the best trade-off between computational times and solution quality. Figure 11 shows the extensive heuristic comparisons for the remaining problem variants W/¬C, ¬W/C, W/C. (c) Heuristic comparison for problem variant W/C Figure 11.: Extensive heuristic comparison for problem variants W/¬C, ¬W/C, W/C Each subplot shows ∆α = α LH − α RO (with α LH , resp. α RO , the solution cost for LH, resp. RO) as a function of T and S, where we limit LH computational times to 1 second. Each subplot corresponds to one of the 9 scenarios resulting from the combination of each area (SF-1 , BER-1 , SF-2 ) with each availability distribution (low-15% , avg-60% , high-90% ). Over all subplots, availability increases from left to right and station density increases from top to bottom.