Computational Methods for Scheduling the Charging and Assignment of an On-Site Shared Electric Vehicle Fleet

We investigate a fleet scheduling problem arising when a company has to manage its own fleet of electric vehicles. Aim is to assign given usage reservations to these vehicles and to devise a suitable charging plan for all vehicles while minimizing a cost function. We formulate the problem as a compact mixed integer linear program, which we strengthen in several ways. As this model is hard to solve in practice, we perform a Benders decomposition, which separates the problem into a master problem and a subproblem and solves them iteratively in an alternating manner. We perform the decomposition in two different ways. First we follow a more classical way, then we enrich the master problem making it stronger but also more complex and the subproblem smaller and simpler to solve. To improve the overall performance, we propose a problem-specific General Variable Neighborhood Search metaheuristic for solving the master problem in earlier iterations. Experimental results show that directly solving the complete mixed integer linear program usually performs well for small to some medium sized problem instances. For larger instances, however, it is not able to find any reasonable primal solutions anymore, while the Benders decomposition scales much better. Especially the variant with the heuristic delivers high quality solutions in reasonable time. The Benders decomposition with the more complex master problem also yields reasonable dual bounds and thus practically relevant quality guarantees for the larger instances.


I. INTRODUCTION
In many companies employees have to perform business trips on a regular basis, e.g., to visit customers.Usually the company provides a fleet of vehicles that are used for these trips.To utilize these vehicles well, they are most commonly shared among the employees, i.e., not every employee gets his or her own dedicated vehicle.As the overall costs per kilometer of electric cars dropped below the costs of cars with combustion engines [1], it seems reasonable to have a fleet that is primarily composed of electric vehicles.Also the reduction of carbon dioxide, which is targeted by more and The associate editor coordinating the review of this manuscript and approving it for publication was Jie Gao .more companies, speaks for an electrification of a company's vehicle fleet.
By having own charging stations and generating energy on-site at the company, e.g., with a photovoltaic system, it is possible to reduce the costs further.Nowadays, the return of investment for photovoltaic systems can be as low as three years, if providing the energy to the grid [2].Using the generated energy to at least partly charge the own vehicles can further speed up the return of investment.
Long charging times, however, pose a challenge for the scheduling of vehicles within a carsharing system that relies on electric vehicles (EVs).Variable on-site energy generation and time-of-use tariffs for the energy from the grid, i.e., where the energy price is time-dependent, make an effective scheduling even harder.Advanced algorithmic approaches are necessary to schedule the charging of a fleet of EVs together with its usage in the sketched context in order to minimize cost while maximizing the availability of the EVs.
We formulate and investigate the EV Fleet Charging and Allocation Problem (EVFCAP) that arises in the carsharing context with EVs and optional on-site energy generation.In this problem we consider a discretized finite time horizon.Given is a fleet of identical EVs and a set of reservations for trips, each having a departure and an arrival time from and back to the site, respectively, and an estimated energy demand.Aim is to assign reservations to vehicles and to determine a feasible charging plan for each vehicle, while minimizing a cost function.It is also possible to leave reservations unassigned, however this is associated with additional cost, e.g., for fulfilling this demand with conventional vehicles with combustion engines or to lend vehicles.In practice the solution methods for the problem can be utilized in a receding horizon approach, in which the scheduling is done for all currently known reservations and a rescheduling is performed whenever a new reservation comes in or there are other changes in the problem input.
Charging can be performed with energy from the grid or surplus energy from the on-site energy generation.For the grid energy a variable time-of-use tariff is assumed.Contrary to grid energy, surplus energy is assumed to not cause any costs as it is the case, e.g., with a photovoltaic system.However only a time-dependent limited amount of surplus energy is available.For simplicity, we do not consider any limit on the maximum power that may be used from the grid or the number of vehicles that can be charged at the same time.However, a respective extension of our models and solution approaches would be straight-forward.
We assume that the time-dependent surplus energy and the energy price are known in advance.Energy prices in timeof-use tariffs are often made public beforehand.If using a photovoltaic plant, the surplus energy can be forecasted using the weather prediction.We consider a myopic scheduling approach, which only considers reservations currently available.If new reservations come in, a rescheduling is supposed to be done.
As we will show later, the problem is NP-hard as well as challenging to solve in practice.As the main contribution of this work we develop several approaches to solve the formalized problem exactly as well as heuristically.We first formulate the problem as Mixed Integer Linear Program (MILP), strengthen it with valid inequalities and, as it turns out to be challenging to solve for larger instances directly, derive a Benders decomposition [3] from this formulation.It is, to our knowledge, the first time a Benders decomposition approach is considered for a problem of this kind.Benders decomposition decomposes a MILP into a master problem (MP) and a subproblem (SP), and both, the MP and the SP are solved iteratively in an alternating fashion.For solving the MP approximately in a faster way, we further propose a heuristic based on General Variable Neighborhood Search (GVNS).For the last iterations, we finally switch from the GVNS to a MILP solver for solving the MP to ensure exactness of the overall approach.We further investigate a second variant of the Benders decomposition.In the classic Benders decomposition the MP contains all integer variables and all other variables are handled in the SP.In the alternative variant we move some variables from the SP to the MP.This way we can also move some of the constraints of the SP to the MP, which strengthens the MP.As a result, the number of required iterations is reduced, but on the downside, using the GVNS for the MP is not possible anymore, so the MP is solved with the MILP solver in all iterations.
We experimentally compare the two Benders decomposition approaches with the original MILP solved by a stateof-the-art solver on a set of artificially generated benchmark instances.The benchmark instances range from small size with up to 5 vehicles and 80 reservations via medium size with up to 20 vehicles and 320 reservations through to large size with up to 100 vehicles and 1600 reservations.While for small and medium-sized instances the direct solving of the MILP performs well and instances with up to 5 vehicles and 80 reservations can be solved to optimality within one hour, for large instances the Benders decompositions in general perform significantly better.On the one hand, especially the Benders decomposition with the GVNS heuristic is able to find sooner much better heuristic solutions for these instances.For instances with 100 vehicles and 800 reservations the gaps to the best known dual bounds are around 7.5% for the Benders decomposition with heuristic, but around 180% when directly solving with the MILP.On the other hand, the variant of the Benders decomposition with the extended MP is able to identify significantly better dual bounds than the original MILP for our largest instances being up to a factor of 10 lower.This article is based on the first author's master thesis [4], to which we also refer for further details on the implementation and evaluation.
To increase the readability we summarize the most important acronyms and symbols used in this article in Table 1.
The remainder of this article is organized as follows.The next section introduces the EVFCAP in a formal way and states it as MILP.Section III discusses related work.In Section IV the EVFCAP is proven to be NP-hard.Section V derives a class of strengthening inequalities.The two variants of the Benders decomposition are presented in Section VI.Section VII introduces the General Variable Neighborhood Search for solving the MP heuristically in a faster way.Subsequently, Section VIII discusses and compares experimental results of the different approaches on artificial benchmark instances.Finally, Section IX concludes everything and gives an outlook on promising future work.

II. PROBLEM FORMALIZATION AND MILP MODEL
In the EV Fleet Charging and Allocation Problem (EVFCAP) we are given • a discretized planning time horizon T = {1, . . ., t max } with each time step having the same length t; • a set V = {1, . . ., n} of n uniform EVs with maximal possible charging power P max > 0 and battery capacity E cap > 0, and for each vehicle v ∈ V the subset of time steps T avail v ⊆ T at which the vehicle is available and the initial energy level E v,0 ∈ [0, E cap ]; • a set R = {1, . . ., r max } of reservations of vehicles, with the set of consecutive time steps T res r = {t start r , . . ., t end r } ⊆ T in which a vehicle is needed for a trip and the expected energy E res r ≥ 0 that will be (approximately) consumed during the reservation for each r ∈ R; note that this energy demand E res r can be estimated based on a planned route of the trip or, if the destination is unknown, crudely on the duration of the reservation; • the grid's costs for electricity c t > 0 per unit of energy during each time step t ∈ T ; • a predicted maximum surplus energy E surmax t ≥ 0 available during each time step t ∈ T , for example from a photovoltaic system; • and costs c uncov > 0 for covering each unit of energy required by unassigned reservations, for example by cars with combustion engines.

A solution comprises
• a partial assignment from reservations to vehicles, modeled by variables x r,v ∈ {0, 1} for r ∈ R, v ∈ V , indicating with value one that reservation r is to be fulfilled by vehicle v, and • a charging plan, modeled by variables p v,t ≥ 0 indicating the power by which vehicle v ∈ V is charged during time step t ∈ T .The objective is to minimize the costs spent for grid energy and uncovered reservations as well as expected future costs.Expected future costs are formulated as the amount of energy that would be needed to fully charge all vehicles after the last time step weighted with a factor of α ≥ 0. These future costs arise from expected reservations after the time horizon that cannot be covered due to a low State-of-Charge (SoC) of the vehicles.
The EVFCAP can be expressed by the MILP ( 1)- (15).The total amount of energy consumed from the grid during time step t is modeled by variables E grid t and the total amount of surplus energy consumed during time step t by variables E sur t .Moreover, to keep the model clear we use additional variables y r , which indicate if reservation r ∈ R is uncovered (y r = 1) or covered (y r = 0).Note that these variables can be kept continuous in the model as they will automatically get integral values due to equalities (2).By T starts and T ends we denote the sets of time steps in which reservations start and end, respectively, i.e., The meaning of the constraints is as follows.Equations (2) prohibit the assignment of a reservation r to multiple EVs, and it forces y r = 1 if reservation r is not assigned to a vehicle.Due to (3), reservations cannot be assigned to EVs that are not available in a time step of the reservation interval.Constraints ( 4) and ( 5) prevent charging EVs in time steps in which they are not available or are used in reservations.
Here the symbol ''\'' denotes the set difference operator, i.e., T \ T avail v represents all time steps in T that are not in T avail v .If overlapping reservations are assigned to a vehicle, the right hand side of (5) becomes negative for the time steps in which the reservations overlap.As this is not possible with the domain of p v,t , constraints (5) also prevent the assignment of overlapping reservations to a vehicle.Note that splitting up the two aspects of (5) into two types of constraints would yield a weaker formulation.Equations ( 6) define the SoC E v,t max of each vehicle after the last time step.Inequalities (7) and (8) determine the SoC of each vehicle for each time step and ensure that it stays in [0, E cap ].It is irrelevant for the correctness of the formulation whether the energy of a reservation is subtracted at the start or at the end of the reservation interval or somewhere in between.In order to make the inequalities stronger, we use the end for (7) and the start for (8).Moreover, in (7), p v,t was added, further strengthening the constraint.Analogously, p v,t has been omitted in (8).If for some time step the SoC cannot exceed the battery capacity and no reservation ends at this time step, then the SoC also cannot exceed the battery capacity for the previous time step.Therefore it is sufficient to limit the SoC only in time steps before a reservation ends and at the last time step.The SoC can only drop when the energy of a reservation is subtracted.Therefore, (8) is only needed for time steps in which a reservation starts.Equations (9) link the power p v,t by which the EVs are charged at each time step with the variables for the energy used from the grid E grid t and the used surplus energy E sur t .Finally, all variable domains are defined in (10) to (15).Note that we use continuous values for the charging rates p v,t as most modern charging infrastructures give this flexibility.In contrast, former work such as [5] often restricts charging to only one or very few power levels, which can ease the mathematical modeling but restricts the possibilities in the scheduling severely.With the larger flexibility of continuously selectable charging rates, we may expect better solutions.

III. RELATED WORK
We are aware of only two papers addressing problems that are closer related to the EVFCAP.These works are discussed at the end of this section.Common problems in the literature that are similar to the EVFCAP are Vehicle Scheduling Problems (VSPs).In these a set of trips, each having a start and end location and a time for departure and arrival, is given [6].The aim is to find a feasible assignment from trips to vehicles that minimizes costs.Vehicles start and end at depots.Problem variants either consider a single depot or multiple depots.The traditional single-depot VSP is solvable in polynomial time, e.g. by using a network flow formulation [7].When considering multiple depots however, the VSP gets NP-hard [8].
Recently a new variant of the VSP arose that considers an electric fleet of vehicles.The electric VSP has been solved exactly by utilizing MILP solvers [9], [10], [11], and with various heuristics such as Adaptive Large Neighborhood Search [12], GRASP and a memetic algorithm [13].The EVFCAP can be seen as a variant of the electric VSP with a single depot, where all trips start and end at the depot and charging is only possible at the depot.Note however that the availability of surplus energy and the possibility that trips are not covered have not been considered in works to the electric VSP we are aware of.Furthermore having multiple locations substantially restricts the solution space, since vehicles have to start the next trip where the previous one ended.Thus VSPs are quite different from our problem and solution methods and results can hardly be compared.
Different kinds of works regarding EVs concentrate on the charging process of a vehicle fleet.On the one hand this concerns smart charging of a vehicle fleet in the context of smart grids, see, e.g., [14] or [15] for a survey from the algorithmic perspective and [16] for a review that additionally considers photovoltaic power production.A variety of problems arise in this context from the different perspectives of the grid operator, the aggregator and the EV user.Proposed algorithms are often distributed to increase scalability and reduce communication overhead.On the other hand non-linearities of the charging process have been considered and investigated.For instance and in particular when considering fast charging, the maximum power that can be used to charge a vehicle depends significantly on its current SoC.This is, for example, considered for the following charging problem from the literature.A single vehicle or a fleet of vehicles is given.Each vehicle has a current SoC, a minimum required SoC, an SoC-dependent maximal charging power, and possibly a departure time.Aim is to find a charging schedule, s.t. the vehicles reach the required SoCs by their departure times or the end of the planning horizon, respectively.Applied solving approaches are the cutting plane method for convex maximal charging power functions [17], [18], different MILP formulations [18], [19], as well as mixed integer nonlinear programming [20].
Although the EVFCAP has overlaps with the above problems, we are only aware of two strongly related works.The following paragraphs give an overview on these.
Sassi and Oulamara [5] consider a problem variant that differs from our EVFCAP setting primarily in the optimization objective.First, they maximize the distance that is traveled with EVs.Then the charging costs are minimized, s.t. the distance stays constant on its optimal value.Another difference is that the number of alternative vehicles that can be used for uncovered reservations, e.g., cars with combustion engines, is limited.Moreover, there are also some differences regarding the charging process: Surplus energy is not considered, i.e., all the energy is taken from the grid.However the grid is limited in the total power it can deliver, so the charging processes of different EVs must also be optimized together.Additionally to the maximum power of the grid and the maximum charging power for each EV, there also exists a lower power limit for charging.Thus, an EV is either charged within its lower and upper power limits or not charged.
The authors propose two heuristics to solve their problem, which they name Sequential Heuristic (SH) and Global Heuristic (GH).The SH repeats the following steps for each EV: • Select a subset of reservations that are non-conflicting s.t. the driven distance is as large as possible and assign them to the EV.
• Find a charging schedule for the EV that minimizes charging costs.
Afterwards, uncovered reservations are assigned to the alternative vehicles.The GH works basically the same way but calculates a charging schedule for all vehicles at once after all reservations have been assigned to vehicles.These two heuristics were experimentally compared to a MILP on random instances.Both, GH and SH, outperformed the MILP on these instances in terms of the quality of heuristic solutions.Compared to SH, GH in general found feasible solutions on more instances and better solutions on the instances where both heuristics found solutions.However, SH was significantly faster.Betz et al. [21] consider a problem variant that differs from our EVFCAP by considering a limited number of available chargers.These chargers can have different characteristics and are located at specific parking lots.If necessary, relocations of vehicles between parking lots are possible but cause additional cost.Betz et al. also formulate their problem as a MILP and solve it using the MATLAB integrated MILP solver.The hardest instance that they could solve within two hours included 30 requests and eight vehicles.Moreover, the authors also investigate the effect on operational costs and environmental impact when using a fleet of EVs compared to a fleet of vehicles with combustion engines.

IV. NP-HARDNESS OF THE EVFCAP
Sassi and Oulamara [5] tried to prove the NP-hardness of their problem variant which they called Electric Vehicle Scheduling and Charging Problem (EVSCP).As the EVFCAP is similar to the EVSCP, the basic idea of their proof can be adopted for the EVFCAP.However, we found that the proof in [5] contains a severe error because the reduction from the used partition problem may result in an exponentially large EVSCP instance as the size of the time horizon T depends on the input values of the partition problem.In the following we correct this issue and adapt the proof for the EVFCAP.
In our proof, we reduce the partition problem, which is known to be NP-complete [22], to the decision variant of the EVFCAP.The decision variant of the EVFCAP asks whether or not there is a solution that has an objective value less than or equal to a given c max .
The partition problem can be stated as follows.

PARTITION PROBLEM Input:
A multiset P = {a 1 , . . ., a m } of positive integers that sum up to 2s with s ∈ N.
Proof: Given an instance of the partition problem, an instance of the decision variant of the EVFCAP is constructed as follows.
• For each integer in P there is a time step.
• There are two EVs.Both EVs have a battery capacity and an initial energy level of s, no charging capabilities and an unlimited availability.
• For each integer in P there is a reservation.All reservations are non-overlapping and the energy consumption of reservation r is a r .
• There are only costs for uncovered reservations.
• The costs has to be less than or equal to zero.
• No surplus energy is available.Assume we have a solution for this problem instance.This solution cannot have uncovered reservations as they would increase the objective value to a positive value.Therefore, the two EVs have to cover all reservations.The total energy consumption of all reservations is 2s.Both vehicles have an initial energy level of s and no recharging possibilities.Therefore the energy level of both vehicles is zero after time step m, and the energy consumptions of the reservations served by a single vehicle sum up to s.This gives a solution for the corresponding partition problem which is Given a solution to the partition problem, a corresponding solution to the EVFCAP is constructed by assigning the reservations corresponding to M v to vehicle v for both v ∈ V .Due to a similar reasoning as above, the solution to the EVFCAP has to be feasible.
Overall, there is a solution to the partition problem, iff there is a solution to the corresponding EVFCAP.As all the described transformations can be done in polynomial time, it follows that the EVFCAP is NP-hard.

V. STRENGTHENING INEQUALITIES FOR CONFLICTING RESERVATIONS
We now introduce strengthening inequalities for the MILP model ( 1)- (15) for the EVFCAP, which prevent the assignment of sets of reservations to a vehicle if those reservations cannot be assigned together to the vehicle due to charging constraints.Let R ⊆ R be a set of non-overlapping reservations and assume that these reservations are considered for joined assignment to vehicle v ∈ V .The maximal achievable energy level E max v,t (R ) of vehicle v for t ∈ T ∪ {0} in dependence of R can be calculated by If E max v,t (R ) is negative at a time step, not all reservations of R can be assigned together to vehicle v.In other words at least one reservation in R may not be assigned to vehicle v in this case, which yields the inequality This kind of inequality in general strengthens the MILP, i.e., there are solutions to the linear programming (LP) relaxation of the model that do not fulfill (18).Take for example the problem instance with one vehicle of capacity and initial SoC 3 and two reservations 1 and 2 requiring energy of E res 1 = E res 2 = 2 in time steps T res 1 = {1} and T res 2 = {2}, respectively.Then the LP relaxation of the MILP has a solution with x 1,1 = x 2,1 = 0.75.However, E max 1,2 ({1, 2}) = −1 is less than zero and inequality 18 is violated for R = {1, 2}.Its addition to the MILP therefore indeed strengthens the LP relaxation.
Clearly, the number of these inequalities in general is exponentially large in |R|.Therefore we consider the dynamic separation of such inequalities in a branch-and-cut approach.The respective separation problem can be stated for a specific vehicle v ∈ V as follows.
SEPARATION PROBLEM Input: Let x LP r,v be the values of variables x r,v in the solution of the LP relaxation for all r ∈ R.Moreover, we use E res r , T res r , P max , t, E cap , E v,0 and T avail v .

Question:
Does there exist a non-overlapping set of reservations For the special case of P max = 0, E cap sufficiently large, and unlimited availabilities (T avail v = T ), the maximal achievable energy can be calculated in an easy way by The condition E max v,t (R ) < 0 then becomes r∈R E res r > E v,0 as it is sufficient to check for t = t max .If no reservations in R overlap with each other, the separation problem corresponds to (the decision version of) the 0-1 knapsack problem with item weights 1 − x LP r,v , capacity one, and item values E res r .Therefore the separation problem is NP-hard.
One way to deal with the separation problem in practice is to discretize all values x LP r,v and to apply dynamic programming over r∈R (1 − x LP r,v ) with the reservations ordered by their start times, maximizing E max v,t (R ) for each intermediate result.Alternatively, one could discretize the energy and minimize r∈R (1 − x LP r,v ).In practice it seems reasonable to apply a fast greedy heuristic to possibly identify a violated inequality and the more complex dynamic programming approach only when the heuristic does not succeed.We implemented a greedy heuristic and tested it with the MILP solved by Gurobi1 in version 9.1.However this did not lead to a speedup in practice, most likely due to the strong generic cuts Gurobi already considers.Therefore we did not use this approach in further steps.
Alternatively, one can initially check all pairs of non-overlapping reservations and provide respective valid inequalities already from the beginning instead of later separating them dynamically.We also tried this approach but again could not observe a practically relevant speedup.

VI. BENDERS DECOMPOSITION
Benders decomposition is a general approach to tackle MILPs of a certain structure and was proposed almost six decades ago by Benders [3].It works by splitting up the problem into a master problem (MP) and a subproblem (SP).For details regarding the Benders decomposition in general, variants of it as well as different speedup strategies and extensions we refer to the review by Rahmaniani et al. [23].
First we apply the Benders decomposition to our problem in a more classical way.Afterwards we strengthen the MP by moving variables from the SP to the MP.In order to ultimately have the whole objective function in this extended MP, we use a reformulated but equivalent objective function.To derive this reformulation from the original objective function (1) the following steps are performed besides basic simplifications: • Replace E v,t max according to (6) As the last two terms are constant, we abbreviate them by Variables E v,t max are not needed anymore when using this reformulation and can therefore be removed together with their defining equation (6).

A. BASIC BENDERS DECOMPOSITION
The following master problem (MP) only considers the assignment variables x r,v , r ∈ R, v ∈ V and y r , r ∈ R and abstracts from the actual charging details.The latter will only be dealt with in the subproblem.The subproblem, however, yields feasibility cuts and optimality cuts which are added to the MP over time to ultimately account for the charging. (MP) s.t.
Feasibility cuts (28) are formulated independently for each vehicle and thus the set of feasibility cuts C feas = v∈V C feas v is partitioned into smaller sets, each dedicated to one vehicle.A feasibility cut is represented by the vehicle v it is dedicated to, a vector γ c of coefficients for x, and a constant γ c .In (29), each optimality cut c ∈ C opt is represented by a coefficient matrix β c of x and a constant β c .Hereby, variable µ represents additional costs for the charging that are to be contributed by the optimality cuts (29).Initially, C feas and C opt are empty.Equation (30) prevents the MP to be unbounded before optimality cuts have been added by providing a trivial lower bound for µ.
The subproblem is determined for a given solution x = (x r,v ) r∈R,v∈V to the MP.We define • for each vehicle v ∈ V the set of time steps when the vehicle is available and not reserved • and for each vehicle v ∈ V the energy difference up to time step t ∈ T caused by assigned reservations.The energy of a reservation may be subtracted at the start or the end of the reservation interval, leading to the two definitions E res r and (34) The primal subproblem SP(x) can now be stated as the following LP. (SP(x)) Further feasibility cuts (28) and optimality cuts (29) are determined by solving the dual SP for specific reservation assignments x = (x r,v ) r∈R,v∈V of a current solution to the MP.More specifically, feasibility cuts are inferred from extreme rays of the dual SP in case this problem is unbounded.Since feasibility cuts can be formulated independently for each vehicle, identifying such cuts is more efficient due to the smaller problem size.Optimality cuts are inferred from extreme points of the dual subproblem.For the dual subproblem and the mathematical derivation of the cuts see Appendix A.
We further remark that it is also possible to formulate the SP as Minimum Cost Flow Problem.We did that and compared the Minimum Cost Flow Problem formulation solved with the network simplex implementation of CPLEX2 in version 12.8.to the LP formulation solved by Gurobi.As it turned out, the network simplex approach was in practice not faster than solving the LP with Gurobi.Therefore we use Gurobi to solve the SP in all further steps.

B. EXTENSION OF THE MP
By adding meaningful information about the SP to the MP from the beginning, we may expect to reduce the number of iterations and, as long as solving the MP does not become substantially harder, the overall runtime.We do this by moving the variables E • As E grid t and E sur t and are available now, the extended MP can use the complete reformulated objective function (23).
• Summing inequalities ( 7) and ( 8) over the vehicles yields constraints, which do not contain p v,t anymore and therefore can be added to the extended MP.
• Optimality cuts are not necessary anymore, as the objective function of the extended MP already comprises all parts of the original objective.
• Feasibility cuts (52) also contain the variables E grid t and E sur t .Note however that now, they cannot be formulated independently for each vehicle anymore.
These changes in the MP cause also some changes in the subproblem and the derivation of the feasibility cuts for which we refer again to Appendix A. Also note that it is still possible to add the feasibility cuts (51) of the basic Benders decomposition, which we will also do in our experiments.Recall that the feasibility cuts of the basic Benders decomposition can be formulated independently for each vehicle, which is not possible after extending the MP.Thus the feasibility cuts of the basic Benders decomposition are not necessarily weaker than the feasibility cuts of the Benders decomposition with the extended MP.

C. SOLVING THE MP
In the classical Benders decomposition approach, the MP is solved with a MILP solver to optimality in each iteration.The process of iteratively solving the MP, the dual SP, and adding derived Benders cuts to the MP eventually stops with a proven optimal solution to the original problem when no further cuts can be derived, given sufficient time and memory.During each iteration, the dual bound derived from the solution to the MP is only valid also for the original problem if the MP was solved to optimality.Solving the MP always to proven optimality is in general time consuming and, at least in earlier iterations, unnecessary, if afterwards a cut is added, especially as in practice typically a large portion of the time is required to prove the optimality of a found solution.We examine two different strategies to speed up the overall Benders decomposition approach.
On the one hand we turn the classical scheme into a Branch-and-Check approach.Here, the MP is only solved once by a MILP solver, but whenever a new feasible solution is found within this search, the respective SP is solved for this solution and obtained feasibility/optimality cuts are dynamically added to the MP.For more information on Branch-and-Check we refer to [24].
On the other hand we speed up the classical Benders decomposition approach by solving the MP first in an inexact way and only switch to an exact approach at a later point when no further cuts can be found by the respective SPs.In this way, the whole procedure stays exact in the end.This idea has already successfully been applied in the literature multiple times, see for example [25] or [26].In our first inexact approach, we solve the MP by the MILP solver terminating early when a solution is obtained with an optimality gap within a certain specified limit.We will see that, although a speedup is achieved for the individual iterations, this ultimately does not lead to substantially better final optimality gaps as more iterations are required.We therefore also consider a metaheuristic approach for solving the MP, which is described in the following section.

VII. GENERAL VARIABLE NEIGHBORHOOD SEARCH HEURISTIC
We suggest a General Variable Neighborhood Search (GVNS) [27] as a heuristic to solve the MP.The GVNS metaheuristic combines basic Variable Neighborhood Search based on shaking moves in larger neighborhoods with Variable Neighborhood Descent (VND).Variable Neighborhood Descent is a generalization of local search that utilizes a sequence of multiple neighborhood structures.Upon finding a better solution in one neighborhood, VND restarts with the first neighborhood, otherwise it continues with the next neighborhood.The eventually returned solution thus is a local optimum w.r.t.all used neighborhood structures.In the outer Variable Neighborhood Search, a shaking operator takes a random solution from a typically larger neighborhood N k ∈ {N 1 , . . ., N k max } of the current solution.Algorithm 1 shows a pseudo code of GVNS.
The specific components of our GVNS for the MP are described in the following.Note that this metaheuristic is dedicated only to the first Benders decomposition variant whose MP ( 24)-(32) contains the variables x r,v and y r only.This is because only then it is possible to reasonably speed up searching the neighborhoods with delta evaluation.Delta evaluation computes the objective value of a neighbor by efficiently computing the difference to the objective value of the current solution based on the move that was used to obtain the neighbor.Due to the possibility to apply delta evalution we expect the GVNS heuristic to work well for the MP.

Remember that an instance of the MP consists of
• the instance of the EVFCAP as well as • a set of feasibility cuts and • a set of optimality cuts.For each optimality cut the factors β c of x and the constant β c are stored.Since the matrix β c usually is dense, all coefficients of the matrix are stored.The same is done for feasibility cuts, with the difference that each cut is targeted to a specific vehicle v and therefore also the vehicle has to be stored besides the factors γ c and the constant γ c .

B. SOLUTION REPRESENTATION
Solutions of the MP are represented by the assignment of reservations to vehicles.Within the GVNS, such an assignment is stored as vector (a r ) r∈R ∈ (V ∪ {ε}) |R| , where a r indicates the vehicle by which reservation r is to be fulfilled or ε if no EV is assigned to this reservation.Note that this representation is more compact and efficient to use than the binary representation with the x r,v and y v variables from the MILP formulation.
To make delta evaluation possible, the GVNS stores additional data for its candidate solutions.In particular • the objective value that corresponds to the assignment, • for each optimality and each feasibility cut the value in which the left hand side of the respective inequality results, and

C. VND NEIGHBORHOOD STRUCTURES
The VND starts with the trivial solution that does not assign any reservations and uses two different neighborhood structures, based on the following moves.
1) Change the assignment of any single reservation.This neighborhood can, for a given solution a, be written as where refers to the hamming distance between the two vectors.Set F contains all feasible assignments, i.e., assignments that neither assign overlapping reservations to the same vehicle nor invalidate Benders cuts.Note that also a r = ε or a r = ε is possible for the reservation r ∈ R whose assignment is changed.The size of this neighborhood is in |N A more general version of the second neighborhood that considers more than two reservations was also considered but turned out to be too time consuming for larger instances and did not yield substantially better results.Therefore we did not use it further.
All moves can be decomposed into alterations of single assignments of a reservation.Thus efficient delta evaluation can be applied to search the neighborhoods that result from the given moves.

D. SHAKING
Shaking is used to escape local optima and to diversify the search.In our case, shaking deletes assignments of a number of reservations, i.e., sets the respective a r variables to ε.This allows the VND afterwards to in general converge to a different solution.Which assignments are deleted is determined in two different ways.In the first shaking method a fixed number of reservations is selected randomly and all their assignments are deleted, and in the second shaking method a fixed number of vehicles is selected randomly and all assignments to these are deleted.
In preliminary experiments we found that the following three shaking configurations, in this order, are performing generally well on the kind of instances we used in our tests.
• Deleting five random assignments from reservations to vehicles.
• Deleting ten random assignments from reservations to vehicles.
• Deleting all assignments to one random vehicle.therefore used these in all further experiments.

VIII. EXPERIMENTS AND RESULTS
We first describe how we created artificial benchmark instances as actual real world instances were not available to us.Then, the results of solving the original MILP as well as the different Benders decomposition variants are presented and compared.

A. BENCHMARK INSTANCES
The benchmark instances were created randomly, with properties of a real-world application at the Honda Research Institute in mind.We consider the following combinations of the number of time steps t max and the number of EVs n: • t max = 32 and n ∈ {1, 2, 5} • t max = 192 and n ∈ {5, 10, 20} • t max = 768 and n ∈ {20, 50, 100} For each of these combinations we consider r max ∈ {4n, 8n, 16n} reservations.The other input values are set as follows.
• The maximum charging power is P max = 3.3 kW and the battery capacity is • For each reservation r ∈ R the reservation energy E res r is selected uniformly at random from [0, E cap ].To determine the reservation interval its length is selected uniformly at random from [1, t max 4 ] then the time interval is placed randomly within the bounds 1 and t max ; actual starting and ending time steps are obtained by rounding down and up, respectively.
• Data on energy prices and surplus energies came from the Honda Research Institute.They were also randomly generated but according to models mimicking some real expectations such as a day-night-cycle and the usage of a photovoltaic system.The available surplus energy is scaled according to the number of vehicles to obtain meaningful instances.
• Costs for uncovered reservations and the factor α in the objective function are set to c uncov = 150 Yen kWh and α = 75 Yen kWh .Thirty instances were created for each combination of t max , n, and r max , yielding 810 instances in total, and we made them available online. 3 Some of the more specific tests were performed on a reduced instance set consisting of only the first ten instances of each instance size to keep the computational requirements reasonable.This will be indicated by referring to the reduced instance set.

B. IMPLEMENTATION AND COMPUTING ENVIRONMENT
We used the MILP solver Gurobi 9.14 [28], which is a leading commercial product.The models as well as the Benders decomposition approach were implemented in Julia 1.65 [29].The Julia package JuMP was used as interface to the Gurobi [30].Each experiment was performed on a single core of an Intel Xeon E5-2640 v4 with a time limit of one hour.We furthermore imposed a memory limit of 36 GB.The time given to the GVNS heuristic is limited to the time that was needed to solve the SP in the previous iteration.We solved the test instances with the introduced approaches in different configurations.These configurations and the results are presented in the following.

C. SOLVING THE MP
Before coming to results directly concerning the Benders decomposition we compare our GVNS to an Iterated Local Search (ILS), which uses the same neighborhood structures as the GVNS for the local search but replaces the more advanced shaking of the GVNS with a perturbation operation that always removes the assignment of 70 reservations in the perturbation step.The value of 70 was determined as a robust choice by a preliminary experimental tuning on independent problem instances.To generate instances for the comparison we solved all instances of the EVFCAP with BDH and of each run we selected one MP instance.
The results of the comparison for a selection of instance sizes are shown in Figure 2. We gave the heuristics five, ten, and 50 seconds for the sizes (t max , n) = (192, 20), (t max , n) = (768, 20), and (t max , n) = (768, 50), respectively, similar to the time they would have when running as part of the Benders decomposition.Each boxplot shows the differences between the final objective values obtained by the GVNS and the ILS.Negative values mean that the objective value from the GVNS is lower and therefore better.As can be seen the GVNS performs slightly better on almost all the instance sizes with a median difference in the objective value of around 2500 for (t max , n, r max ) = (768, 20, 320).But as we will see later, these are still rather minor differences -both approaches can be said to perform reasonably well.

D. IMPACTS OF INDIVIDUAL IMPROVEMENTS
For the rest of this section we will refer to the Benders decomposition with BD if the MP is solved to optimality in each iteration, with BDG if the MP is solved by Gurobi only up to a certain gap limit in the first iterations, and with BDH if the MP is solved by the GVNS heuristic in the first iterations.Concerning BDG, we remark that in our implementation the constant offset cconst in the reformulated objective function was not considered as it does not affect the search space and the feasibility and optimality of solutions; thus, the gap-limit for the early termination was determined with cconst = 0.However, in all results listed in the following, the correct values for cconst as calculated from the original instances are taken into account.Superscripts to BD, BDG, or BDH denote the following: • ext: the extended MP is used • cuts: the inequalities from Section V for pairs of vehicles are added from the beginning to the MP • bnch: Branch-and-Check was used instead of the classical Benders decomposition scheme We first evaluate for BDG ext the impact on the number of solved instances, i.e. instances for which a feasible solution was found by the approach, and final optimality gaps when allowing different gap-limits for solving the MP in the beginning.The results for the reduced test set and two selected instance sizes are shown in Figure 3.Note that if the gap is zero the approach is identical to BD.As can be seen, only one instance out of ten can be solved with BD for the large instance size, but more than five instances out of ten can be solved for any of the other evaluated MP gaps.The number of solved instances increases slightly with the used gap, rising from six solved instances for 0.5% MP gap to nine solved instances for 10% MP gap.Contrary to this, the final gap gets worse with an increased MP gap for smaller instance sizes.As a compromise we select an MP gap of 2% for further experiments.
Next we compare BDG ext to BDG ext,cuts , BDG ext,bnch , BDG, and BDH .All experiments were conducted on the reduced test set.A comparison of the number of iterations respectively number of solved instances between BDG ext and BDG ext,cuts on a small, medium-sized, and large instance size can be seen in Figure 4.The number of iterations on the small instances is up to 15% lower for BDG ext,cuts and this also leads to a lower runtime.However the version that does not add static cuts in the beginning is able to solve more instances of large size.Especially for size (t max , n, r max ) = (768, 100, 1600) BDG ext can solve eight instances out of ten while BDG ext,cuts cannot solve any instance.We are more interested in the performance on large instances and therefore we chose to not add these cuts for the further experiments.Figure 5 shows a comparison of runtime, final gap, number of iterations respectively callback calls, and the number of solved instances between BDG ext and BDG ext,bnch .Note that the number of callback calls in BDG ext,bnch correspond to the number of so far feasible solutions for which the respective dual SP was solved in order to possibly obtain further cuts.We can see that the number of callback calls in BDG ext,bnch is significantly higher than the number of iterations in BDG ext for most of the instances.While for small instances BDG ext,bnch is faster than BDG ext by a factor of two to three, the gap is larger by a factor of up to four on some of the medium-sized instances and no large instances with r max ≥ 800 were solved by BDG ext,bnch whereas BDG ext is able to solve most of these instances.
The reason for the worse performance of BDG ext,bnch on larger instances is that the SP needs to be solved much more often and with many inferior solutions, which wastes more of time.
We now investigate the impact on the dual bound when shifting the variables E grid and E sur from the SP to the MP as done in the extended MP. Figure 6 compares the dual bounds produced by BDG ext and BDG on some selected instance sizes.Since the objective function is minimized, the higher the dual bound the better.While the dual bounds are similar for smaller instances, BDG ext gives significantly better dual bounds for larger instances, especially when the ratio between r max and n is large.For r max ≥ 400 BDG is only able to give negative dual bounds, which are not reasonable because  the original objective function (1) cannot yield values less than zero.Recall that the Benders decomposition is based on the reformulated objective function (23), and the trivial dual bound for the MP with this objective function can be negative.The variant BDG ext is able to find positive dual bounds on these instances and we will see in the next subsection that these bounds are nearly optimal.
Finally we also investigate the impact of using the GVNS heuristic for the MP in the first iterations.Figure 7 compares BDH with BDG ext on different instance sizes.On most of the small instances with t max = 32 the runtime for BDH is by 50 to 100 % larger than for BDG ext .Contrary to this, both, the number of solved instances and the objective values are in general significantly better for larger instances.For (t max , n) = (192, 20), e.g., BDG ext is not able to solve any instances with r max ≥ 160 while BDH solves them all.Also the objective values obtained by BDG ext for instances with (t max , n) = (768, 50) are up to a factor of nine higher than the ones obtained by BDH .In other words Gurobi is not able to find any reasonable solutions for the MP for these large instances, while the GVNS performs quite well.
Of all the considered improvements, the ones discussed in the previous two paragraphs are most notable.Moving the variables E grid and E sur from the SP to the MP seems to significantly improve the resulting dual bound on large  instances.Similarly using the heuristic for the MP seems to significantly improve the objective value for large instances.Therefore, in the following we will compare in more detail the two approaches BDG ext and BDH with each other and with the MILP formulation (1)- (15) directly solved by Gurobi.

E. COMPARISON OF THE MILP AND THE BENDERS DECOMPOSITION APPROACHES
As we saw in the previous section, especially the Benders decomposition variants BDG ext and BDH are most promising.We therefore evaluate and compare these with the direct solving of the original MILP, however, also using the reformulated objective (23) to get rid of the variables E v,t max , on the whole benchmark instance set.For short we abbreviate the latter direct MILP solving just by MILP.Detailed results are given in Table 2.These results are visualized for three selected instance sizes in Figure 8.Each boxplot represents the final objective values respectively dual bounds for a specific approach applied to the test instances of the given size.Let us first look at Figure 8a, where the approaches are compared for the medium sized instances with ten EVs and 40 reservations.The MILP outperforms the two Benders decomposition approaches for these instances with a median gap of 0.1% compared to median gaps of 3.3% and 13.5%, respectively.This changes for larger instance sizes as can be seen in Figure 8b and Figure 8c, where the approaches are compared on instances with 100 EVs and 400 and 800 reservations, respectively.Here BDH provides excellent solutions with gaps to the dual bounds found by BDG ext around 15% while the other two approaches only find solutions that are not much better than the trivial solutions that leave all reservations unassigned.These insights coincide with those of Sassi and Oulamara [5].They also find that the MILP formulation solved by a solver does not lead to (reasonable) solutions for large instances, whereas their heuristic is able to find good solutions for these instances.BDG ext gives good dual bounds while BDH finds worse or no dual bounds.The original MILP only yields a meaningful dual bound for some of the instances, while for others like the ones with r max = 800 it only finds poor trivial dual bounds.When looking at the detailed results in the table, they confirm that for large instances BDH finds the best solutions, while BDG ext finds the best dual bounds.One exception are the largest instances with (t max , n, r max ) = (768, 100, 1600) where BDH does not find feasible solutions.The reason is that not enough iterations are performed to get to an MP solution that is also feasible in respect to the charging constraints.BDG ext performs better here because the MP solutions are worse but already feasible with charging constraints from the beginning.
We also want to point out that the fraction between r max and n, which indicates whether the system is overbooked, has a high influence on the runtime and gap.This can be seen particularly well for BDG ext for the instances with n = 10 and n = 20.Here, solutions to almost all instances with r max = 4n could be found, but not many solutions for r max = 8n and none for r max = 16n.Last but not least, Table 3 breaks the mean objective values down into the values of its terms, i.e., the costs for energy taken from the grid, the costs for uncovered reservations and future costs.The solutions of BDH were taken for that purpose, since they are consistently near optimal.Naturally the costs for uncovered reservations dominate the two other terms if r max = 16n since then the system is overbooked and most of the reservations cannot be served with the available EVs.
Furthermore the larger the instances get, the lower gets the ratio of costs for uncovered reservations versus other costs, provided that the ratio between r max and n stays the same.While for instances with (t max , n, r max ) = (768, 100, 400) these costs are zero, for instances with (t max , n, r max ) = (32, 1, 4) costs for uncovered vehicles are roughly four times the other costs together.We believe the reason for this behavior is that in larger instances random fluctuations become less relevant and there are more possibilities to compensate for them and thus more reservations can be assigned to EVs.In practice this means that a larger system reaches a higher utilization of EVs.It should also be mentioned that in general, future costs do not rise in proportion to the number of EVs.For example, set r max = 8n; then the average future costs are 1056 for n = 1, 4575 for n = 10 and 11899 for n = 100.A larger system therefore also tends to be able to leave the EVs at a higher SoC at the end of the time horizon.This can also be seen in more detail in Figure 9, where each boxplot represents the final SoCs of all vehicles of solutions to instances of a specific instance size, when solved with BDH .The vehicles in the solutions to all instances with t max = 192 or t max = 768 have a median final SoC of more than half of the battery capacity and can therefore be considered well prepared to serve reservations of the next time horizon when expecting a similar demand of requests.In a schedule of the next time horizon EVs can already serve reservations at the beginning leading to a higher utilization of the EVs.We furthermore observe that a higher ratio between the number of reservations and the number of vehicles leads to a lower final SoC, which is an expected behavior.

F. IMPACT OF CHANGING THE DEGREE OF TIME DISCRETIZATION
Finally we investigated the impact on the solution quality and hardness of the instances when using coarser time steps.To obtain a coarser instance from an existing instance we arranged the time steps in blocks of 2, 4, 8, and 16 time steps, respectively, and merged the time steps of each block into one.Reservation intervals and vehicle availabilities are extended to span whole blocks, surplus energies of time steps in a block are added up, and for the electricity cost the mean value of the cost of the individual time steps in the block is used.The derived problem instances are more restrictive in the sense that solutions to them can be transformed into feasible solutions to the original instance but not vice versa.Therefore these coarser instances have worse or equally good solutions while being smaller.
Experimental results are shown in Figure 10.As expected the coarser instances are easier to solve, with a factor on the median running time of up to 100 for the instance size (t max , n, r max ) = (192, 5, 80) when having blocks of 16 time steps.However merging the time steps increases the median optimal objective value by more than 50% for the instance size (t max , n, r max ) = (192, 5, 20) and still by more than 22% for the instance size (t max , n, r max ) = (192, 5, 80).While the reduced running time for the coarser instances looks promising, the optimal objective value of these instances is significantly worse.

IX. CONCLUSION
We considered a fleet scheduling problem that was not addressed in the scientific literature so far.After the formal definition we proved the problem to be NP-hard.We formulated the problem as a Mixed Integer Linear Program (MILP), taking care in the details to achieve a strong linear programming relaxation.Moreover, we suggested an additional class of problem-specific strengthening inequalities, although it turned out that in practice they do only have a minor impact when using a modern MILP solver.
As an alternative to directly solving this MILP, we proposed a Benders decomposition for the problem.The basic variant of this Benders decomposition was enhanced by multiple measures.For one of these measures the variables relating to the energy consumption at each time step were moved from the subproblem (SP) to the master problem (MP).This allowed to also move related inequalities from the SP to the MP making the linear programming relaxation of the MP stronger.Moreover, we proposed a heuristic based on General Variable Neighborhood Search (GVNS).The heuristic is used to solve the MP in the first iterations only approximately but significantly faster.By switching to an exact method for the MP afterwards the approach ultimately still yields an optimal solution, given enough time and memory.Furthermore, we investigated the effect on the performance when using Branch-and-Check, which solves the MP only once and adds Benders cuts as lazy constraints, instead of the classic iterative approach.
We experimentally compared the Benders decomposition to solving the complete MILP on a set of artificial benchmark instances and evaluated the effect of the different performance improvement measures of the Benders decomposition.It turned out that directly solving the MILP performs better than the Benders decomposition on small to medium sized instances.On large instances, however, the solving of the original MILP fails to find reasonable solutions and meaningful dual bounds.The Benders decomposition based on the extended MP is here able to give better dual bounds and the Benders decomposition with the GVNS finds better primal solutions.Using Branch-and-Check and adding our strengthening inequalities did not improve the performance of the Benders decomposition.Note that the GVNS heuristic does not work for the extended MP and therefore it is not directly possible to combine the Benders decomposition based on the extended MP with the GVNS.This aspect may be addressed in future work since such an approach might yield both, good solutions and good dual bounds at the same time.However, appropriately extending the GVNS is not straight-forward.One possibility to solve the issue would be a nested Benders decomposition, i.e., solving the MP of the Benders decomposition with extended MP by another Benders decomposition.While the outer Benders decomposition is able to give good bounds, it is still possible to use our heuristic to find good solutions to the MP of the inner Benders decomposition.
Future work may also investigate a more fine-grained approach for the amount of effort that is used to solve the MP in the different phases of the whole optimization.It seems reasonable to start with little effort and raising it whenever the derived cut does not cut off or worsen the found MP solution.Finally, the algorithm will again solve the MP with the MILP solver to optimality to guarantee exactness of the whole approach.When using the MILP solver for the MP, one may adapt the allowed gap from higher to lower values.For the heuristic the allowed runtime may be adapted.
Last but not least, it would also be interesting to investigate generalizations and extensions of the problem in future work.Practically relevant aspects are to consider uncertainties, in particular in the energy price or the available surplus energy, or to consider non-linear maximal charging rates.Further possible extensions are the consideration of also delivering energy from the vehicles back to the grid, or to allow different vehicle types.Vehicles could, e.g., have different battery capacities or even specific characterics that may be required by some of the reservations.In the original MILP model vehicle-specific battery capacities can be directly considered in a straight-forward way.Different further vehicle-specific attributes, such as passenger capacities, can also be given for the vehicles and respective constraints can be added.While adapting the models in these ways is relatively easy, actual running times for solving them may differ significantly, either in the negative or positive sense.

APPENDIX A DERIVATION OF THE BENDERS CUTS
This part of the appendix deals with the derivation of the Benders cuts.First feasibility and optimality cuts for the basic Benders decomposition are derived.Afterwards we derive feasibility cuts for the Benders decomposition with the extended master problem.

A. BASIC BENDERS DECOMPOSITION
Benders cuts will be derived by solving the subproblem for specific reservation assignments x = (x r,v ) r∈R,v∈V in a current solution to the MP given by ( 24) to (32).
The primal subproblem SP(x) from Section VI is shown again in Figure 11.To get feasibility and optimality cuts, we have to derive the dual of this subproblem.In the primal problem SP(x), the λ-variables in parentheses denote the dual variables that correspond to each constraint.The resulting dual subproblem is shown in Figure 12.
This dual subproblem can be simplified in multiple regards as follows.
• Variable λ p=0 v,t appears only once in (71) and is unbounded.Therefore, it can always be chosen in a way that the inequality is fulfilled.Thus, λ p=0 v,t as well as (71) can be removed.
• Variables λ sur t and λ p v,t appear in the objective function with positive coefficients and in two constraints that only give upper limits.The value of such a variable will always be the smaller upper limit in an optimal solution.This leads to the equivalent formulation given in Figure 13.
For reasons of clarity we define for a given solution to DSP(x), denoted by λE (78)

1) FEASIBILITY CUTS
Feasibility cuts are inferred from extreme rays of the dual subproblem in case this problem is unbounded.Let λlb v,t , λub v,t , and λE t represent such an extreme ray.We then have to have λE t ≤ 0 and the objective function without constants (i.e., without α) has to be positive.The third term of the objective function is zero in this case.Furthermore, λE t = 0 due to a similar argument as before.The formulation can now be split into different parts for the different vehicles, which are independent from each other.As the total objective value is positive, some of these parts also have to provide positive contributions and therefore also form rays.  Hence the problem of finding extreme rays can be decomposed into smaller problems, one for each vehicle.Using the obtained rays, we get the following feasibility cuts for vehicles v ∈ V .The factors γ c,r and the constant γ c can be obtained from these cuts by comparison of the coefficients of x.  (84) Again, β c,r,v and β c can be derived by comparison of coefficients.

B. EXTENSION OF THE MP
Starting point for the Benders cuts is now the MP given by (44) to (56).The dual subproblem changes in some regards compared to the dual subproblem of the basic Benders decomposition.
• The term with λ sur t in the objective drops because E sur t is no longer part of the subproblem.
• E sur t and E grid t turn into constants, leading to a new term in the objective function.
• The constraints corresponding to E sur t and E grid t are dropped.After simplification the dual of the updated subproblem can be stated as shown in Figure 14.
Note that there are no optimality cuts it this case as the variables from the SP do not occur in the objective function of the MP.
further introduce the new cost factors ct := c t − α and cuncov := c uncov − α, (22) and arrive at the reformulated objective function min t∈T ct E grid t + cuncov r∈R E res r y r −α t∈T E sur t + cconst .
grid t and E sur t , for all t ∈ T , from the SP into the MP.The number of these variables is low enough (in O(|T |)) to still keep the MP reasonably small, but now it becomes possible to add further strengthening inequalities to the MP that make use of E grid t and E sur t .Besides moving the variables E grid t and E sur t from the SP into the MP, the MP is updated as follows to the formulation (44)-(56) shown in Figure 1.

FIGURE 1 .
FIGURE 1.Extended master problem after moving the variables E grid and E sur from the SP to the MP.
1 (a)| = O(|R| • |V |).2) Swap the assignments of any single pair of reservations:N 2 (a) = {a ∈ F | a = a • σ rr , r, r ∈ R}, (58)where permutation σ rr ∈ S R swaps r with r and leaves all other elements at their positions, and F denotes the set of feasible assignments again.The size of this neighborhood is in|N 2 (a)| = O(|R| 2 ).For the instances we considered |N 2 (a)| is generally larger than |N 1 (a)|, and therefore we apply N 2 after N 1 in the VND.

FIGURE 2 .
FIGURE 2. Difference in the objective values of solutions obtained from GVNS and ILS.A negative value means that the objective value of the GVNS is (better) than the value of the ILS.

FIGURE 3 .
FIGURE 3. Numbers of solved instances and final gaps for different MP gap-limits in BDG ext for instances of two sizes.

FIGURE 4 .
FIGURE 4. Comparison of BDG ext with and without additionally adding the static cuts from Section V.

FIGURE 5 .
FIGURE 5. Comparison of BDG ext with the branch-and-check variant BDG ext,bnch .

FIGURE 6 .
FIGURE 6.Comparison of the dual bounds obtained from BDG ext and BDG.

FIGURE 8 .
FIGURE 8. of final objective values and dual bounds of the different methods for three selected instance sizes.Trivial primal bounds are obtained from the solutions that leave all reservations unassigned.The trivial dual bounds in 8b and 8c are calculated by the right hand side expression in(30).

TABLE 2 .TABLE 3 .
Comparison of the MILP, BDG ext and BDH.The number of feasibly solved instances / optimally solved instances, the median final objective values, and the median dual bounds are listed.Bold numbers indicate the best values among the three compared approaches.Mean objective values of solutions obtained by BDH and their composition in respect of the three terms of the objective function (1).

FIGURE 9 .
FIGURE 9.Final SoC of EVs for some selected instance sizes.

FIGURE 10 .
FIGURE 10.Impact of changing t on the objective value and the solving time when directly applying the MILP solver.
1) FEASIBILITY CUTS Let λp v,t , λub v,t and λlb v,t be a solution to the above dual subproblem that leads to a positive objective value.The feasibility cuts can then be formulated as shown in Figure 15.Once again, γ ext c,r,v , γ grid c,t , γ sur c,t and γ ext c can be obtained by comparison of coefficients.
(2)lace v∈V x r,v with 1 − y r according to(2).This leads to the expression . • Replace t v∈V p v,t with E • t∈T ends ∪{t max } Optimality cuts are inferred from extreme points of the dual subproblem.Let λ E v,t and λE t be such an extreme point.Unlike for extreme rays, we cannot discard constants so the above simplifications do not work here.The resulting cuts are v∈V t∈T ends ∪{t max } Dual subproblem after extending the MP.
FIGURE 15.Feasibility cuts after extending the MP.