National-scale bi-directional EV fleet control for ancillary service provision

Deploying real-time control on large-scale fleets of electric vehicles (EVs) is becoming pivotal as the share of EVs over internal combustion engine vehicles increases. In this paper, we present a Vehicle-to-Grid (V2G) algorithm to simultaneously schedule thousands of EVs charging and discharging operations, that can be used to provide ancillary services. To achieve scalability, the monolithic problem is decomposed using the alternating direction method of multipliers (ADMM). Furthermore, we propose a method to handle bilinear constraints of the original problem inside the ADMM iterations, which changes the problem class from Mixed-Integer Quadratic Program (MIQP) to Quadratic Program (QP), allowing for a substantial computational speed up. We test the algorithm using real data from the largest carsharing company in Switzerland and show how our formulation can be used to retrieve flexibility boundaries for the EV fleet. Our work thus enables fleet operators to make informed bids on ancillary services provision, thereby facilitating the integration of electric vehicles.


A. Background and motivation
Public authorities and the private sector face many challenges in transforming industries and infrastructure to meet sustainability goals.A key factor is the successful integration of renewable energies such as solar or wind power, which however poses difficulties to the power system due to the increased fluctuations in supply from renewable energy sources.At the same time, an increasing number of electric vehicles pose an additional burden on the grid [16].Both challenges inspired the development of smart charging or V2G technologies, where the charging flexibility of EVs are exploited as buffer storage to the power system.Smart charging and V2G were shown to have high potential benefits for peak load shaving [7,19,36], supporting the integration of renewable energies [22] while offering additional revenues to vehicle owners [18].
Although smart charging and V2G have been studied for years [12,29], they remain difficult to implement in practice for the following reasons: 1) they require control over a sufficiently large fleet of EVs, 2) they imply complex dispatching problems, and 3) they involve trading between the power system and the vehicle fleet operators.A major opportunity is the application of V2G for large-scale car lorenzo.nespoli@supsi.chsharing systems [11], since they can centrally manage large and significant resources for V2G operations.In contrast to the share of EVs on the private vehicle market (8% global sales share 1 ), the share of EVs in car sharing systems is already high, with more than 66% of car sharing services offering fully or partially electric fleets [27].V2G may afford additional revenues to car sharing operators, but at the same time requires careful dispatching to minimize the negative impact on car availability for mobility purposes.
Here, we propose an optimization approach for V2G operations that scales to a large fleet of EVs.Specifically, we first provide a monolithic formulation for optimizing charging schedules, and further develop relaxations that allow to decompose the problem by aggregated vehicle hubs such as car sharing stations.Our experiments demonstrate a strong improvement in runtime using our approach, enabling its application on a large-scale vehicle fleet.Furthermore, the optimization framework is tested on a new dataset from a car sharing operator in Switzerland.It is shown that our method scales to a fleet of 1440 electric vehicles in feasible runtime and can be employed to decrease energy costs while providing different kinds of grid services.Our optimization approach is therefore not only relevant for car sharing services but may in general support in controlling V2G fleet operations.

B. Literature review and previous works
An increasing number of works is tackling the problem of charging schedule optimization in the context of car sharing; Xu et al. [34] optimize charging times in a MINLP problem targeted at determining the fleet size of a car sharing system.He et al. [15] optimize the charging station setup and schedule for a car sharing fleet and provide interesting insights on the best decisions on charging station placement and minimum State of charge (SOC).Similarly, [3] formulate a two-step optimization problem in order to reduce the charging prices in a shared system, while retaining user satisfaction.Only some research has focused on large-scale, national level optimization of V2G, since this is a more challenging problem if realistic constraints are considered.Furthermore, the typical scale of pilot projects in this context is small: in [25] the authors reviewed 54 pilot projects using EVs for providing grid services, reporting an average number of 26 EVs per pilot.
In [21] a decentralized algorithm to optimize the charge (but not discharge) of 5000 EVs was presented.In [37] the authors present a rule-base two-stage hierarchical approach to coordinate charging operations of thousands of EVs.While this research only considers smart charging and not V2G, [5] also include the possibility of V2G in the relocation-optimization of one-way car sharing.In [39] the authors coordinated 500 EVs to achieve frequency regulation using a rule-based control in a V2G setting.[38] regard the problem that is most closely related to our formulation, namely V2G strategies for car sharing, and they propose a two-stage stochastic optimization employing a 24 hours receding horizon approach solved with a resolution of 15 minutes.They show that keeping integer variables lead to infeasible solution times (greater than 32 hours in their case), and propose to both relax all integer variables to continuous one and use decomposition techniques in order to speed it up.However, they do not provide a scalability analysis of their algorithms, nor mention the number of considered EVs.In contrast to optimal control methods, others propose data-driven optimization with learning methods.For example, [8,20,30,31,32] train a reinforcement learning (RL) agent to decide on charging behavior.However, these methods are usually focused on finding decision policies for single EVs, since finding the optimal joint actions for a fleet of EVs, which is the focus of our work, is a much more challenging task, in general requiring a multi-agent RL strategy, which usually involve to optimize over a large decision space.Authors in [26] propose RL for guiding charging decisions for a whole vehicle fleet at once by reducing the action space by pooling EVs with similar energy requests; however, this was done not considering external inputs such as an aggregated profile, and disregarding V2G.

II. PROBLEM DEFINITION AND FORMULATION
In the following we start describing a generic formulation needed to effectively synchronize the EV fleet charging and discharging operations, and later explain how relaxing some conditions can lower the overall computational complexity.The common setting for all the problem formulations is the following: a car sharing provider operating a stationary fleet (as opposed to free-floating) is willing to jointly optimize all its EVs' operations in order to reduce its own operating costs, whether by optimizing for a dynamic price, increasing its own self-consumption if local PV generation is present, or by providing services to the electric grid.Furthermore, the provider knows at least an approximated schedule of the future EV locations, in terms of their presence at a given charging station and driven mileage for the next control horizon.This can be realistically achieved using information from booking apps and by modeling historical data.Based on these assumptions we can estimate the lower bounds for the EVs' battery energy constraints needed to satisfy all their foreseen mobility demand, as we will show in section III-B.These time series are required to formulate the optimal control problem, as explained in the following section.

A. Monolithic formulations
Given a control horizon of T steps, n s stations, each station hosting n v,s vehicles, and called T and S the sets of times and stations, the monolithic problem can be described as: v∈Vt,s where x ∈ R T × s nv,s is the matrix containing the battery state for all the EVs in kWh.For sake of clarity, table I reports all the parameters and optimization variables X of the problem with associated dimensions and domains.Here F (u) : R T × s nv,s → R and Q(x) : R T × s nv,s → R are two scalar convex functions.In particular F (u) is a cost function associated with the charging and discharging actions of the EVs and depends on the specific business model and will be further specified in section II-B.We now explain in detail the problem constraints.Equation 2 describes the EVs dynamic equation, taking into account self-discharge and asymmetric charging and discharging efficiencies encoded in the A v ∈ R and B v ∈ R 2 discrete dynamics matrices, obtained by the continuous one through exact discretization [28]: where and η sd , η ch and η ds are the characteristic self-discharge constant, charge and discharge efficiencies, respectively.Since B c defines an asymmetric behaviour in charging and discharging (even with equal charging/discharging coefficients), solving the battery scheduling requires to use two different variables for the charging and discharging powers for each EV.These are concatenated and denoted as a whole as u = [u c , u d ], where u d , u c ∈ R T,nv are charging and discharging operations for all the EVs in kW.∆e ∈ R T,nv is the (sparse) matrix containing the energy lost during the last EV trip, defined as: where the first condition in equation (9) where k is a large constant, which allows to retrieve feasible solutions even if some EVs are not fully charged.Equation (3) states that charging and discharging variables u c and u d are positive quantities.Equation ( 4) makes use of the binary variable x c , which indicates whether a given EV is charging, to encode the bilinear constraint u c u d = 0, where is the Hadamard product; this encodes the fact that each EV cannot charge and discharge simultaneously.It must be noted that this condition is sometimes naturally satisfied by the problem, depending on the objective function F (u), as shown for example in [13].However, this is not always guaranteed; for example if we want to implement peak shaving in the presence of PV power plants.In this case EVs could occasionally decide to both charge and discharge and exploit the round-trip efficiency to dissipate more power and perform valley filling when the overall station network is a net energy producer.The same reasoning can be applied to quadratic profile tracking, as in the case of tracking a given power profile for providing services to the grid.In equation ( 5), the binary variable c ∈ R T ×nv is used to enforce charging and discharging powers to be zero when the car is not located at a station.Finally, called V l,s the set of EVs located at station s at time t, U s the rectangular box set of power limits at station s, the last two equations ( 6) and ( 7) represent the station constraints on maximum power and available number of charging stations, respectively.The problem composed by equations 1 -7 is very general, however it is computationally expensive; due to the presence of the soft constraint on the minimum required energy (10) (and to the possible quadratic objectives included in F (u)), the problem belongs to the MIQP class, with a number of variables in the order of O(T n v ), where in our case n v is in the order of 10 3 and T is equal to 96, since we consider 15 minutes steps and a daily control horizon.We now discuss how the original problem can be simplified by relaxing or removing some of the constraints 4 -7, and the implications for the problem's formulation hypothesis.
Strictly stationary mobility model: if the sharing model is strictly stationary, meaning that the EVs are permanently assigned to a charging station and can only be plugged there, we can relax equations ( 6) and ( 7) which encode the maximum power and connection limits per station.These can be rewritten as: The only difference to equations ( 6) and ( 7) is that the set V s is no more time dependent.This effectively removes the interlink between different stations given by EVs travelling between them; in other words, sets of EVs belonging to different stations will not influence each other directly, but only by means of the system-level objective F (u). Since the rest of equations ( 3) -( 4) do not interlink stations, the problem can be easily decomposed.It must be noted that the original problem can also be decomposed; however, if the mobility model is not strictly stationary, it is likely that the influencing graph between EVs is dense, meaning that the behaviour of a given EV can be influenced by a high number of other EVs, dependent on the routing between stations.This will require to introduce decoupling variables for all the states and control variables, which involves a message passing of variables in the order of O(T n v n s ) at each iteration.On the contrary, when F (u) is an aggregate function, as in all the cases presented in this paper, decomposing the problem requires messages with size in the order of O(T n s ) at each iteration.Since n s << n v and n v is in the order of thousands, the strictly stationary hypothesis will results in a data transmission reduction in the order of 10 4 .
Stations are not downsized: each station has enough chargers to accommodate all its assigned EVs at the same time.This hypothesis, combined with the previous one, allows us to remove completely the binary variable c indicating whether an EV is connected to a charger.In fact, equation ( 7) is not needed anymore, and equations ( 5) can be replaced with: where l is the location matrix parameter, with entries l t,v equal to 0 if the v th vehicle is not located in any stations at time t.
EVs are mono-directional: this hypothesis will not allow to consider direct discharge of EVs into the main grid, nor energy arbitrage between EVs.Considering currently available solutions this is the setting with lower technological burden which could already be implemented by most EV car sharing providers.Note that it will be still possible to provide services to the grid by modulating the overall charge.This hypothesis will simplify the dynamics equations, removing the discharging variable u d .As a result, bi-linear constraints (4) can be dropped, removing the binary variable x c .If this hypothesis is combined with the two previous ones, the overall problem becomes linear or quadratic, depending on the form of F (u), allowing to use a larger set of solvers and substantially reducing the computational complexity.

B. Decomposition and business models
In this section we show how the original problem can be decomposed by stations under the hypothesis of a strictly stationary mobility model and that stations are not downsized.As we keep the bidirectional hypothesis, we still need to include the bilinear constraint u c u d = 0, handled by equations ( 4) and by the integer variable x c .In the next session we will discuss alternative methods to handle this bilinear constraint.Under the aforementioned hypothesis the problem can be decomposed using the alternating method of multipliers (ADMM) [4].Following the standard ADMM procedure, since we want to decompose per station, we should introduce n s auxiliary variables representing the total power at each charging station.However, since in our case we are only interested in objective computed at the aggregation level of stations or for the overall fleet, F (u) can be written in the form where S is a system level objective, that is the objective to minimize at fleet level, and C is a cost function that should be minimized at station level.Here p s (u) = ps,load − ps,pv + v∈Vs u c,v − u d,v is the sum of forecasted base load and PV production (if any) for station s and the sum of the charging and discharging operations of all EVs belonging to s. Considering this form for F (u), we need to introduce only one additional variable z ∈ R T representing the average power of the n s controlled stations.The final problem before the decomposition can be written as: s.t.( 2), (3), ( 4), ( 13), ( 12), (11) ( 16) We can then proceed to formulate the augmented Lagrangian objective function in scaled form: Since problem ( 15)-( 17) can be seen as a sharing problem, we can further simplify the standard ADMM following the description in [4] for this specific case.As the choice of ADMM's parameter to achieve a good convergence rate can be problematic under the presence of equality constraints, we use a slightly different form, namely the linearized ADMM [14,35]; briefly speaking, this form introduces a quadratic penalty for deviating from the decision actions at the previous iteration.We can then write the minimization in the primal and dual variables update as: s.t.
(2), ( 3), ( 4), ( 13), ( 12), (11) ( 21) where s are the vectors of operations and states of all the EVs belonging to station s.Following [4], + λ k are the reference signals for the u and z update.Line (20) contains the dumping term of the linearized ADMM form for the primal variable u s update, γ being a dumping parameter.
The two functions C(p s (u s )) and S(zn s ), representing respectively the station and the fleet objectives, can be used to tackle different business models.For example, for the station level, the following cases can be easily considered: • Minimize energy costs.Called p buy ∈ R T and p sell ∈ R T the time-dependent buying and selling prices in cts/kW h.In the presence of local generation e.g.due to PV power plants at the station's location, the cost function can be either positive or negative, depending on the overall power at a given time and can be expressed as in equation (24).
The cost can be thought of as the maximum over two affine functions (the first and second line of equation (24), respectively).If p buy is always greater than p sell we can minimize energy costs by introducing an auxiliary variable y ∈ R T representing the station's energy costs.We can restrict the feasible space for y to the epigraph of the cost function C(p s (u s )) by adding the two following constraints to the station problem ( 19)-( 21): Minimizing y then guarantees that its value at the optimum, y * , will lie on the epigraph's lower boundary (and will thus represents the prosumer's total costs).In this case C(p s (u s )) = T t y t δt/3600 where δtis the considered time step.Even without setting a systemlevel objective, this strategy can result in some EVs performing arbitrage, charging at low price times and later discharging to other EVs if the price swing is high enough to compensate for the round-trip efficiency. .However, pure peak shaving has usually no economic drivers; the fleet manager is usually interested in reducing its total costs rather than having a flat profile per-se.Since peak tariffs are usually computed on the maximum power peak attained on a monthly basis, a more appropriate approach could be to implement a lexicographic strategy, at first minimizing the station's economic costs and then using the optimal cost found in this first step as a constraint for a second optimization in which a peak shaving objective is minimized.At the same way, the system level objective S(zn s ) can be used to address several fleet-level business cases: • Intra-day cost minimization.In the case in which the fleet manager has a deal to buy energy at intra-day costs, it can follow the same strategy illustrated to the cost minimization objective at station level and set S(zn s ) = T t y t δt/3600.• Profile tracking.A standard quadratic profile tracking can be used to make the fleet dispatchable setting S(zn s ) = T t (zn s − r) 2 , where r is a reference profile to be tracked.However, to quantify revenues from grid regulation services and flexibility calls, a linear cost function is more appropriate, as equation ( 31) that we used in the presented case study in section III-D.

C. Bilinear constraints handling
We now present the proposed method to handle the bilinear constraint u c u d = 0 inside the ADMM iterations of the decomposed problem ( 19)-( 23), without using the integer variable formulation encoded in equation ( 4).Linear complementarity constraints arise in a variety of problems from bilevel optimization to eigenvalue complementary problems.Given a scalar objective function f (x, y) of two variables x, y ∈ R T + , the simplest form of the complementarity constraint problem can be written as: where z = [x T , y T ] T .Depending on the complexity of the underlying problem, which is in general NP-hard, different iterative methods exist to find a feasible solution or a stationary point for this kind of problem [17].One of the most used strategy is the one implemented in the YALMIP package for Matlab, which uses the built-in solver for nonconvex problems BMIBNB.The procedure sequentially finds refinements of an upper and a lower bounds for the problem, respectively found using a local non-linear and a convex solver.
The next iteration is then found using a standard branchand-bound logic and split the feasible space into two new boxes [1].The convex approximation for bilinear problems is found using a McCormick formulation.In [6], the authors proposed tighter bounds for bilinear problems exploiting Mc-Cormick relaxations and a sequence of MILP problems.The McCormick envelope has been also proposed for the relaxation of factorable functions by systematic subgradient construction [23], a concept similar to automatic differentiation.In this work we have chosen a different approach relying on the following observation: since we are solving the main problem iteratively, we want to exploit an iterative relaxation running in parallel with the standard ADMM iteration, without relying on branch and bound methods.Running a partial optimization for one part of the objective function for ADMM is theoretically justified by the generalized form of ADMM (GADMM) introduced in [9].The GADMM guarantees the convergence even in the case in which the local (stations') problems are only partially solved.This allows us to use a first order Taylor expansion around the previous solution to approximate the complementarity constraint x y = 0, in combination with a standard ADMM using Lagrangian relaxation.We can write the first order Taylor expansion around the previous solution as: We propose to use this to minimize f (z) while respecting the constraint, as reported in algorithm 1.Here w is an auxiliary Algorithm 1: Taylor relaxation Input: z 0 = [x T , y T ] T , w 0 , λ chosen at random, parameters ρ, γ 1 while stop condition not met do variable representing x y, which we want to shrink to zero; lines 2-4 are standard ADMM iterations where line 3 is the analytical solution of the minimization of the Lagrangian function with respect to w; finally line 5 is a dumped iteration over the last solution, with dumping parameter α.A different approach is proposed by Wang et al. in [33], where they provided algorithm 2, which is a standard application of ADMM to two objective functions, f (z) and I x T y=0 , where I x T y=0 is the feasible set for the complementarity constraint.Contrary to algorithm 1 that we propose, this approach guarantees that the problem always satisfies the complementary constraint at each iteration, due to the projection onto the feasible space of I x T y=0 at line 3.The authors proved that algorithm 2 converges into a stationary point for the bilinear constrained problem when f (z) is a smooth function.Algorithms 1 and 2 are appealing since they are easily implementable and don't require to sequentially explore the whole solution space with a branch-and-bound strategy.

A. Data analysis and preprocessing
We test our optimization framework on a dataset made available by a car sharing operator managing a fleet of around 3000 vehicles.The dataset covers all car reservations from 1st of January 2019 until 31st of July 2020, thereby including the period before the COVID-19 pandemic as well as the first wave.In total, there are around 2 million bookings during this period, comprising 140880 unique users and 4461 vehicles.Due to the setting of the considered car sharing service, only a small fraction of trips are one-way (0.3%), and during the observation period only 3.5% trips involved electric vehicles.Furthermore, the number of vehicles per station is low on average in the considered system.73% of all stations offer a single vehicle, further 15% only two vehicles.5% of all stations have five or more vehicles.The limited availability of parking slots per station also explains the low fraction of one-way trips.We first analyze the flexibility of vehicles for V2G operations based on their daily and overall demand.Figure 1 shows the histogram of reservations by vehicle.Clearly, there are strong differences in the usage patterns of different vehicles.48% of the vehicles have at least one reservation in less than 50% of the days.These findings imply a strong opportunity for the car sharing operator to utilize its fleet for V2G.However, the most flexibility is given during the night: Figure 2 shows a bell shaped curve of vehicle utilization over the course of a day, peaking in the afternoon.On average 21% of vehicles are reserved at any time.Last, we validate the assumption that most car reservations are known in advance, as it is necessary for optimizing the charging schedule.Concerning the spontaneity of the bookings, around 34% cars are reserved more than a day in advance, whereas 20% of the reservations are done less than an hour before the reservation period.
The data are discretized to a temporal resolution of 15minute steps.We remove cancelled trips but include service reservations necessary for relocating vehicles.We use the reservation period in contrast to the actual driving period to define the time span of car usage.However, this leads to overlapping trips in some cases when a returned vehicle was taken by the next user before the end of the original reservation period.The reservation period is therefore cut to the end of the previous drive / start of the next drive if necessary.Reservations without a ride are assumed to be cancelled and are not taken into account.

B. ICE mobility patterns and State of Charge modeling
The car sharing service operator has set the ambitious goal to electrify their entire fleet by 2030.In order to provide a realistic simulation of the future fleet, and to demonstrate how our optimization approach scales with the number of stations, we propose to utilize the booking patterns of ICE vehicles as projected EV usage patterns, under the assumption of a similar driving behavior.Since only 3.5% of all trips are EV trips, this scales up the number of reservations by a factor of more than 25.In consultation with the car sharing operator we assign an EV model to each ICE vehicle based on the car category in the car sharing operator service, i.e. "Budget", "Combi", "Transporter" etc.For example, all vehicles of the category "Transporter" were simulated as Mercedes-Benz eVito vehicles, and all in category "Budget" were assigned the VW e-up model.
Two pieces of information are needed as input to the optimization problem: When a vehicle is plugged in at a station, and the required state of charge at the start of a reservation.Due to the modeling of ICEs as EVs and the lack of SOC data in the provided dataset, we approximate the latter by the number of driven kilometers.Given the vehicle specifications (i.e.battery range and battery capacity) we compute the required SOC by multiplying the number of driven kilometers with the average energy consumption.

C. Formulations comparison
We evaluated the numerical advantage of the proposed formulations in two steps.At first, we compared the monolithic formulation ( 15)-( 17) to the decomposed one ( 19)-( 23) using integer variables for handling bilinear constraints.In a second step, we evaluated the decrease in computational time in using the proposed linear methods for the bilinear constraints in the decomposed problems.For both these comparisons we vary the range of total EVs and the horizon length.The stations' objective function was set to energy cost minimization, while the system level objective was set to a profile tracking with a zero reference profile.The results of the first comparison are reported in the heatmaps of figure 3.For this comparison, we solved the monolithic problem using GUROBI with standard absolute and relative tolerances, while the stopping criterion for the decomposed formulation is a joint condition on the primal and dual residual, as described in §3.3.1 of [4], using abs = 1e − 6 and rel = 1e − 4, respectively.The first two heatmaps refer to the total computational time of the decomposed problem and the monolithic formulation, respectively.The last plot shows the ratio of the two, a value lower than one meaning a lower computational time for the decomposed formulation.As expected, the computational advantage over the monolithic formulation increases with both the number of EVs and the length of the horizon.The experimental data for up to 360 vehicles shows a clear trend; the computational time of the decomposed problem for the most time consuming configuration being roughly 20% of the time needed by the monolithic formulation.The second comparison was done using a fixed number of iterations, which was set to 800.At first, we tuned the parameters of algorithm 1 and 2 w.r.t. the solution reached by the integer formulation, using a random sampling strategy over the configuration with 144 EVs and an 18 steps horizon.The parameters (ρ and γ for 1 and ρ for 2, respectively) were then held constant over the different combinations of EVs and horizon lengths.We found that both the algorithms' performance was stable for a large range of parameters values.The computational times are shown in figure 4, where the first heatmap refers to the Taylor relaxation, the second one to the integer formulation and the last is the ratio of the two.As the computational advantage is due to the change of the class of the problem from MIQP to QP, we found a negligible difference in the computation times between algorithm 1 and 2, and thus here report only results for the Taylor relaxation.Also in this case there is a clear trend in the reduction of computational time with increasing number of EVs and steps.The highest reduction was found for the most time consuming configuration of 577 EVs and 18 steps, with the Taylor relaxation using roughly 35% of the time needed by the integer formulation; once again we expect this value to get lower for problems with higher number of EVs.
Figure 5 shows the distribution of ∆ abs,rel J c for all the cases reported in figure 4.Here J c is defined as the sum of the different objective functions without including any augmented Lagrangian terms (neither the one deriving by the problem decomposition nor the ones of the linear formulations) in order to have a fair comparison: Both the algorithms converge to the solution of the integer formulation with some oscillations, even if the Taylor-based  relaxation shows better convergence, achieving a relative difference in the order of 1e − 3 for all the cases after 800 iterations.

D. Economic results
We use the proposed algorithm 1 to retrieve flexibility boundaries for an EV fleet.The setting is the following: an EV manager bidding for ancillary services is interested to know for a given leading time how many MWs, for how long, can be requested to the EV fleet for both upward and downward flexibility calls, and how much it costs per MWh.This information can then be used by the manager to make more informative bids.We followed the approach proposed in [24] to achieve hourly flexibility boundary for an aggregation of office buildings.For each hour of the day, we solve the optimization problem ( 19)- (23), where each station minimizes its total energy costs for the EVs charging operations, and C(p s (u s )) is modeled through the auxiliary variable y as explained in section II-B.Since the considered car sharing operator's stations are located under different Swiss DSOs, we used data from [10] to link them with the correct values for the buying and selling energy prices p buy and p sell , depending on their location.Additionally, we probabilistically assigned each station with a PV power plant, with a nominal power proportional to the maximum number of hosted EVs at that station.The system level objective function is set to be: where r is the reference profile, T h is the set of timesteps belonging to hour h and p f is the price of flexibility, which is constant over the considered hour.Equation ( 31) can be seen as a linear punishment in deviating from a flexibility call.We simulated a total number of 1440 EVs, keeping all the EVs belonging to a given station if the latter was chosen by a random sampling among all the available ones.An example of results using this objective function when h = 12, at different price levels, is shown in figure 6.When the fleet receives an upward flexibility call at noon, the consumption decreases in the rest of the day w.r.t. the baseline profile in which the system level objective is set to zero and the only objective is the stations' cost minimization.The opposite verifies when the fleet receives a downward flexibility call.For a given day, we run 24 optimizations, systematically changing T h and repeat the process for different values of p f .The resulting flexibility envelopes can be seen in figure 7. Lines of different colors represent the convex envelopes of the maximum and minimum flexibility attained at different hour of the days for a given price p f .It can be seen how during the first hour of the day the fleet is not prepared to an upward call, since the average SOC of the fleet is too high and the fleet has no time to discharge beforehand.Moreover, a saturation effect can be noticed after a given level of price: the maximum attainable flexibility does not change significantly passing from a p f of 255 CHF/MWh to 377 CHF/MWh.In order to better analyze this effect, we considered more price levels for the case in which flexibility is requested at noon. Figure 8 shows the maximum amount of MW reached for 10 different values of p f ranging from 10 CHF/MWh to 377.5 CHF/MWh.The saturation effect is clear for both the upper and lower requests, but it's starting at slightly different price levels, around 210 and 250 CHF/MWh, respectively.Finally, we study the effect of the flexibility request on the other considered costs in the optimization problem. Figure 9 shows the change in charging costs, loss of SOC (equation ( 10)), tracking revenues and total costs for the noon case.As expected, as the price level increases, the tracking revenues rises for both upward and downward flexibility calls, but this comes at the expense of higher charging costs.The change of cost for the SOC lost is negligible compared to the other costs.IV.CONCLUSIONS In this paper we presented an optimization model to control the charging and discharging operation of large EV fleets.We started by modeling a generic case in which the EVs are allowed to relocate between stations, and then focused on the strictly stationary model where EVs are picked up and dropped off at the same station, since this reflects the conditions of the presented case study.For this last case we demonstrated how the problem can be decomposed by stations, allowing to reduce the overall computational time.Furthermore, we used iterative methods to handle the bilinear constraints arising from the V2G formulation, which allows us to use a larger class of (free) solvers.For different combinations of horizon's lengths and number of EVs, we reported numerical results showing substantial speed ups w.r.t. the monolithic formulation, due to both problem decomposition and the use of relaxations for the bilinear constraints.We see multiple opportunities for future work.First, many car sharing bookings are spontaneous, limiting the applicability of day-ahead planning in real world scenarios.This could be tackled with the integration of booking forecasts; since forecasts introduce uncertainty, a receding horizon optimization can be used to minimize errors.Additionally, a stochastic formulation e.g.tree-based stochastic MPC [2], can be used to further tackle the uncertainty of bookings and PV generation.

Fig. 2 .
Fig. 2. Reserved vehicles by time of the day

Fig. 3 .
Fig. 3. Computational time for different number of timesteps and considered EVs for the decomposed (left plot), the monolithic formulation (center plot) and the ratio of the two (right plot).

Fig. 4 .
Fig. 4. Computational time for different number of timesteps and considered EVs for the decomposed problem using the Taylor bilinear relaxation (left plot), the integer formulation (center plot) and the ratio of the two (right plot).

Fig. 5 .
Fig. 5. Comparison of convergence dynamics using the Taylor or Wang formulation for the bilinear constraint relaxation, in terms of relative differences in total objective w.r.t.integer formulation, when stations optimize for costs and the fleet has a reference tracking objective.Confidence interval refers to all the 42 combinations of horizon length and number of EVs of figure 4.

Fig. 6 .Fig. 7 .Fig. 8 .
Fig.6.Example of response to upward and downward flexibility calls as a function of price, compared to the baseline case in which there is no system level costs and the stations just optimize for their local energy prices.

Fig. 9 .
Fig. 9. Behaviour of different fleet costs as a function of flexibility price p f , for the noon case.
designs times in which the location matrix has a positive discrete derivative, that is, when the v th EV connects to a charging station.Here e ∈ R T ×nv is the (sparse) energy constraint matrix, containing the energy that the EVs require at departure times, while t d (t) is the last departure time seen at step t.In other words, the minimum energies required at departure times and encoded in e are equal to the energy drops ∆e t,v needed to be reintegrated at next arrival time.The energy requirements stored in e are assumed to be known at solution time for the next solution horizon, and they are estimated starting by the total driven km for the last trip, as explained in section III-B.Since it is not always possible to guarantee that all the EVs satisfy the energy requirements stored in e at departure time, state constraints on the EVs SOC are taken into account as a threshold soft constraints encoded in Q(x): ) is less than one, i.e. if self-discharging is considered, this will result in a delay-charging strategy, pushing charging operations closer to EV departure times.• Minimize charging times.If we want to charge EVs up to their required SOC at departure as soon as possible, we can minimize C(p s (u s )) = t p s (u s )d(t) where d(t) is a convex discount function weighting less initial steps.• Perform peak shaving.The most straightforward way is to set C(p s (u s )) = p s (u s ) 2 2 • Maximize self consumption -minimize energy imports from the grid.This can be achieved setting C(p s (u s )) = t p s,t (u s,t ).If the term A v in the dynamic state equation (2