## Day-ahead Operation of an Aggregator of Electric Vehicles via Optimization under Uncertainty

´ lvaro Porras Cabrera A OASYS group University of Malaga Malaga, Spain

alvaroporras19@gmail.com

Ricardo Fern´ andez-Blanco OASYS group University of Malaga Malaga, Spain

Ricardo.FCarramolino@uma.es

Salvador Pineda Dept. Electrical Engineering University of Malaga Malaga, Spain spinedamorente@gmail.com

Abstract -We assume an aggregator of electric vehicles (EVs) aiming to buy energy in the day-ahead electricity market while accounting for the technical aspects of each individual EV. We pose the aggregator's problem as a bilevel model, where the upper level minimizes the total operation costs of the fleet of EVs, while each lower level minimizes the energy available to each vehicle for transportation given a certain charging plan. Thanks to the totally unimodular character of the constraint matrix in the lower-level problems, the model can be mathematically recast as a computationally efficient mixed-integer program that delivers charging schedules that are robust against the uncertain availability of the EVs. Finally, we use synthetic data from the National Household Travel Survey 2017 to analyze the behavior of the EV aggregator from both economic and technical viewpoints and compare it with the results from a deterministic approach.

Index Terms -Electric vehicles, aggregator, electricity market, bilevel programming, robust optimization

## I. NOMENCLATURE

- A. Sets and Indices
- V Set of electric vehicles, indexed by v .
- T Set of time periods, indexed by t .

## B. Variables

- c v,t Charging power of electric vehicle v in period t , kW.
- e v,t Energy state of charge of vehicle v in period t , kWh.
- s + v,t Slack variable for the energy balance of the battery of electric vehicle v in period t , kWh.
- s -v,t Slack variable for the energy balance of the battery of electric vehicle v in period t , kWh.
- z v,t Auxiliary variable used to linearize c v,t α v,t , kW.
- α v,t Availability of electric vehicle v to charge, being 1 if available and 0 otherwise.
- ξ wc v Energy required for transportation of electric vehicle v throughout the optimization horizon for the worst-case scenario.
- ζ v Dual variable associated with the constraint ensuring a minimum number of periods in which the electric vehicle v is available.
- β v,t Dual variable associated with the lower bound on the availability of electric vehicle v and time period t .

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 755705). This project has also been supported by Fundaci´ on Iberdrola Espa˜a 2018. n

Juan Miguel Morales Dept. Applied Mathematics University of Malaga Malaga, Spain

Juan.Morales@uma.es

- β v,t Dual variable associated with the upper bound on the availability of electric vehicle v and time period t .

## C. Parameters

- C v Maximum charging power of electric vehicle v , kW.
- E v Maximum state of charge of electric vehicle v , kWh.
- E v Minimum state of charge of electric vehicle v , kWh.
- N T Number of time periods of the optimization horizon. P Penalty cost, e /kWh.
- η v Charging efficiency of electric vehicle v .
- ∆ t Time interval, h.
- λ t Day-ahead electricity price in time period t , e /kWh.
- ̂ α v,t Expected availability of electric vehicle v in period t .
- ̂ ξ v,t Expected energy consumption of electric vehicle v in period t due to motion, kWh.
- α v,t Maximum value of the availability of electric vehicle v in period t .
- α v,t Minimum value of the availability of electric vehicle v in period t .
- K v Minimum number of periods in which the electric vehicle v must be available throughout the time horizon.

## II. INTRODUCTION

According to the Directive 2018/844 of the European Parliament and of the Council of 30 May 2018 [1], the electrification of transport will be ' promoted in the near future while enabling further development at a reduced cost in the medium to long term .' Thus, the electromobility is expected to impact the operation and planning of future electricity systems and electric vehicles (EVs) will play an important role [2]. Within this context, new actors will appear in the power system, e.g. aggregators, in order to manage the operation of a fleet of EVs in residential and non-residential areas. The successful rollout of EVs will undoubtedly bring environmental and societal benefits, however such aggregators will face challenges related to their operation and planning.

There is a vast literature on the impact of EVs on the distribution system [3]-[5] and on the charging strategies of a single EV or a fleet of EVs [6]-[11]. This paper is focused on the latter one, specifically at residential level. It is obvious that a growing penetration of EVs in a residential district will lead to an undesirable over-consumption at certain time periods and the charging of the EVs will be affected due to feeder limitations. In addition, the drivers' habits such as arrival and departure times are stochastic in nature and unpredictable. The treatment of this type of uncertainty is still an unresolved issue.

Halvgaard et al. [6] implemented an economic model predictive control to simulate the charging strategy of an EV by considering perfect forecasts of the driving patterns. Iversen et al. [7] handled the uncertainty on the driving patterns for a single EV with an inhomogeneous Markov model. However, upgrading such stochastic dynamic programming model to a fleet of EVs may be complex. Sarker et al. [8] put forward a model for a charging station coupled with a Battery Energy Storage System (BESS) wherein there exist bidirectional communication between the BESS and the grid. In this paper, the technical and physical characteristics of each EV were neglected, and thus the aggregated EV operation was represented by a charging demand of the fleet, which was modelled via stochastic optimization. Ortega-Vazquez [9] proposed an optimization problem to solve the scheduling of an aggregator at household level while accounting for both the technical characteristics of the EV and their battery degradation costs with an emphasis on vehicle-to-grid services. However, it neglected the stochastic nature of the driving patterns of an EV. Baringo and S´ anchez Amaro [10] analyzed the bidding strategy problem for an EV aggregator participating in the dayahead electricity market while accounting for the uncertainty on the driving requirements. The EV aggregator was modelled as a virtual battery, thus ignoring the individual operation of each EV of the fleet; and the uncertainty was represented with confidence bounds on the aggregated demand. Similarly, Gonz´ alez Vay´ a et al. [11] analyzed the same problem but the uncertainty on driving patterns was represented by means of scenarios.

This paper aims to model the day-ahead operation of a residential aggregator of a fleet of EVs. The technical and physical characteristics of each EV are taken into consideration individually based on the users' information, i.e., the aggregator seeks to optimize the amount of electricity to be purchased in each time period. In addition, the uncertainty on driving patterns is modelled via robust optimization due to the limited information the aggregator will be able to collect. The main contributions of this paper are threefold:

- · We propose a novel mathematical model to operate a fleet of EVs in a residential area while accounting for the uncertainty on driving patterns. Moreover, the operation of each EV is incorporated into the problem formulation.
- · We model the uncertainty on driving patterns by using robust optimization so that the total energy required for transportation for each EV must be satisfied throughout the optimization horizon when considering the worstcase scenario in terms of its unavailability. This takes inspiration from the robust contingency-constrained unit commitment proposed in [12].
- · We exhaustively compare the proposed robust optimization problem with its deterministic counterpart.

The remainder of this paper is organized as follows. Section III presents the proposed day-ahead model for the aggregator of EVs in a residential area. Section IV describes the methodology to transform the original problem into a single-level mixed-integer linear equivalent. Section V provides the deterministic formulation for the aggregator's problem. Section VI presents the case study and discusses the results. Finally, some relevant conclusions are duly drawn in Section VII.

## III. PROBLEM FORMULATION

We consider an aggregator of EVs in a residential district. The usual pattern of this pool of EVs during weekdays is to depart from home to work in the early morning and to arrive home in the evening. The aggregator of EVs aims to minimize their total operation costs and it can be mathematically expressed as follows:

$$\underset { \text{gear} } { \text{eue} } \underset { \text{gator} } { \text{in} } \sum _ { \Xi ^ { R B } } \sum _ { v \in \mathcal { V } } \sum _ { t \in \mathcal { T } } \lambda _ { t } c _ { v, t } \Delta t + \sum _ { v \in \mathcal { V } } \sum _ { t \in \mathcal { T } } P \left ( s _ { v, t } ^ { + } + s _ { v,$$

subject to:

$$\underset { \text{ver, it} } { \text{costs} } \quad e _ { v, t } = e _ { v, t - 1 } + \Delta t \eta _ { v } c _ { v, t } \alpha _ { v, t } - \widehat { \xi } _ { v, t } + s _ { v, t } ^ { + } - s _ { v, t } ^ { - }, \\ \text{of an} \quad \forall t \in \mathcal { T }, v \in \mathcal { V } \quad \quad \quad$$

$$\text{iding} \quad c _ { v, t } \leq \overline { C } _ { v }, \quad \forall v \in \mathcal { V }, t \in \mathcal { T }$$

$$\mathrm e \, \underset {. } { \mathrm d } \, \text{day} _ { \mu } \quad \underline { E } _ { v } \leq e _ { v, t } \leq \overline { E } _ { v }, \ \forall v \in \mathcal { V }, t \in \mathcal { T }$$

$$\inf _ { A \sim A } \ e _ { v, N _ { T } } = e _ { v, 0 }, \ \forall v \in \mathcal { V }$$

$$\underset { \text{ion of} } { \text{aiea} } \quad c _ { v, t }, s _ { v, t } ^ { + }, s _ { v, t } ^ { - } \geq 0, \quad \forall v \in \mathcal { V }, t \in \mathcal { T } \quad \quad \quad$$

$$\text{sented} _ { \text{ilarly}, } \quad \xi _ { v } ^ { w c } \geq \sum _ { t \in \mathcal { T } } \hat { \xi } _ { v, t }, \quad \forall v \in \mathcal { V }$$

$$\mathbf e _ { \mathbf f } \quad \xi _ { v } ^ { w c } = \min _ { \alpha _ { v, t } } \left \{ \sum _ { t \in \mathcal { T } } \Delta t \eta _ { v } c _ { v, t } \alpha _ { v, t }$$

$$\dot { \colon = } & \min _ { \alpha _ { v, t } } \left \{ \sum _ { t \in \mathcal { T } } \Delta t \eta _ { v } c _ { v, t } \alpha _ { v, t } \\ & \sum _ { t \in \mathcal { T } } \alpha _ { v, t } \ge$$

where the set of decision variables is Ξ RB = ( c v,t , e v,t , s + v,t , s -v,t , α v,t , ξ wc v ) .

The robust optimization model is driven by the minimization of the total costs (1), which comprise two terms: (i) the cost of the energy bought by the aggregator in the day-ahead market, and (ii) the penalty costs when the equation associated with the energy state-of-charge evolution is violated. This objective function is subject to technical constraints. Expressions (2) model the energy state-of-charge evolution of the EV v at time period t , while taking into account the energy required for transportation. These expressions include the variable α v,t to account for the availability or unavailability of the EVs. Constraint (3) imposes the maximum charging rate for electric vehicle v and time period t . The maximum and minimum bounds for the energy state-of-charge are set in constraints (4). Expressions (5) enforce boundary conditions on the energy state-of-charge of the EVs. Constraints (6) define the nonnegativity character of the upper-level decision variables.

Constraint (7) enforces that the total demand due to transportation must be satisfied throughout the time horizon for the worst-case scenario of the EVs' availability for each EV. This worst-case scenario is obtained by minimizing (8), i.e., the energy available for transportation to EV v throughout the optimization horizon, over the uncertainty set of its availability, which is defined by (9)-(10). Dual variables are shown in parentheses after a colon. Constraint (9) sets the minimum number of time intervals in which the electric vehicle v must be available throughout the time horizon. Constraints (10) serve us to easily enforce the availability or unavailability of the electric vehicle v in time period t by setting α v,t = α v,t = 1 or α v,t = α v,t = 0 , respectively. Parameters K v , α v,t , and α v,t are estimated based on historical records. More generally, though, these parameters should be tuned to reflect the riskaversion attitude of the aggregator.

## IV. METHODOLOGY

The lower-level decision variables α v,t are binary and the lower-level problems become non-convex. Note, however, that the constraint matrix of each lower-level problem (8)-(10) is totally unimodular (see conditions (iii)-(iv) in Thm. 19.3 of [13]). Besides, we assume that the parameters K ,α v v,t , α v,t are integer. Hence, we can relax the binary character of α v,t while ensuring that there exists an optimal solution for which these variables take integer values. Thus, the relaxed lowerlevel problems are linear programs and thus, the original bilevel problem is transformed into a single-level equivalent by using results from duality theory of linear programming [14], i.e., we apply the so-called duality-based approach. In a nutshell, we replace each lower-level optimization problem with its primal feasibility constraints, its dual feasibility constraints, and the equality corresponding to the strong duality theorem. Note that ξ wc v is replaced either by the primal or dual lower-level objective function. Finally, we add back the binary character of the decision variables α v,t .

## A. Dual Lower-level Problem

The dual problem of (8)-(10) for each electric vehicle v is:

$$\max _ { \zeta _ { v }, \underline { \beta } _ { v, t }, \overline { \beta } _ { v, t } } \ K _ { v } \zeta _ { v } + \sum _ { t \in \mathcal { T } } \left ( \underline { \alpha } _ { v, t } \underline { \beta$$

subject to:

$$\zeta _ { v } + \underline { \beta } _ { v, t } + \overline { \beta } _ { v, t } \leq \Delta t \eta _ { v } c _ { v, t }, \ \forall t \in \mathcal { T } \quad \quad ( 1 2 ) \quad \text{ such}$$

$$\underline { \beta } _ { v, t } & \geq 0, \overline { \beta } _ { v, t } \leq 0, \quad \forall t \in \mathcal { T } & ( 1 3 ) \quad \ e _ { v, t }$$

$$\zeta _ { v } & \geq 0. & \quad \ \ ( 1 4 ) & \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

## B. Robust Single-Level Formulation

The resulting single-level robust formulation can be mathematically expressed as:

$$\min _ { \Xi ^ { R S } } \text{ Objective function (1)} \quad \ \ ( 1 5 ) \quad \text{where}$$

subject to:

$$\text{Constraints } ( 2 ) - ( 6 ) \text{ \quad \ } ( 1 6 ) \text{ \quad \ } \text{descr}$$

$$\underset { \text{on for} } { \text{trans-} } \quad K _ { v } \zeta _ { v } + \sum _ { t \in \mathcal { T } } \left ( \underline { \alpha } _ { v, t } \underline { \beta } _ { v, t } + \overline { \alpha } _ { v, t } \overline { \beta } _ { v, t } \right ) \geq \sum _ { t \in \mathcal { T } } \widehat { \xi } _ { v, t }, \forall v \in \mathcal { V } \ ( 1 7 )$$

$$h \, E V. \quad \text{Constraints } ( 9 ) \neg ( 1 0 )$$

Constraints (12)-(14)

(19)

$$\text{$ut the} \\ \text{bility}, \quad \begin{array} { c } K _ { v } \zeta _ { v } + \sum _ { t \in \mathcal { T } } \left ( \underline { \alpha } _ { v, t } \underline { \beta } _ { v, t } + \overline { \alpha } _ { v, t } \overline { \beta } _ { v, t } \right ) = \sum _ { t \in \mathcal { T } } \Delta t \eta _ { v } c _ { v, t } \alpha _ { v, t }, \\ \text{num} \\ \text{$\dots$} \end{array} \right ) \quad \begin{array} { c } \cdots \end{array} \cdots$$

$$\mu \, \text{must} \quad \alpha _ { v, t } \in \{ 0, 1 \}.$$

Constraints (16) are identical to the upper-level constraints (2)-(6), constraints (17) are identical to (7), lower-level primal and dual feasibility constraints are given in (18)-(19), expression (20) is the equality associated with the strong duality theorem, and the binary character of α v,t is imposed in (21). Problem (15)-(21) is characterized as a nonlinear program due to nonlinear products of continuous and binary variables (i.e., c v,t α v,t ) in constraints (2) and (20). By using integer algebra results [15], we can replace those constraints with the following set of linear ones:

$$& \text{s} \\ & \text{if} \quad & e _ { v, t } = e _ { v, t - 1 } + \Delta t \eta _ { v } z _ { v, t } - \widehat { \xi } _ { v, t } + s _ { v, t } ^ { + } - s _ { v, t } ^ { - }, \\ & \text{$_t$} \quad & \forall t \in \mathcal { T }, v \in \mathcal { V }$$

$$( 2 2 )$$

$$\mathfrak { c h } _ { \mathfrak { c r } _ { - } } ^ { \cdot } \quad K _ { v } \zeta _ { v } + \sum _ { t \in \mathcal { T } } \left ( \underline { \alpha } _ { v, t } \beta _ { v, t } + \overline { \alpha } _ { v, t } \overline { \beta } _ { v, t } \right ) = \sum _ { t \in \mathcal { T } } \Delta t \eta _ { v } z _ { v, t },$$

$$\int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int_{ 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0} \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \ int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ { 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int_ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0 } \int _ 0.$$

$$\text{ent} \quad 0 \leq c _ { v, t } - z _ { v, t } \leq ( 1 - \alpha _ { v, t } ) \, \overline { C } _ { v }, \ \forall t \in \mathcal { T }, v \in \mathcal { V } \quad ( 2 4 )$$

$$\inf _ { \ell } \quad 0 \leq z _ { v, t } \leq \alpha _ { v, t } \overline { C } _ { v }, \ \forall t \in \mathcal { T }, v \in \mathcal { V }. \quad \quad ( 2 5 )$$

Note that the set of decision variables Ξ RS = ( c v,t , e v,t , z v,t , s + v,t , s -v,t , α v,t , β v,t , β v,t , ζ v ) . The single-level mixed-integer linear robust optimization problem is hereinafter referred to as RO-EV.

## V. DETERMINISTIC FORMULATION

The proposed robust optimization problem RO-EV is compared with its deterministic counterpart, namely DO-EV. The deterministic problem DO-EV can be expressed as follows:

$$\min _ { \Xi ^ { D O } } \sum _ { v \in \mathcal { V } } \sum _ { t \in \mathcal { T } } \lambda _ { t } c _ { v, t } \Delta t + \sum _ { v \in \mathcal { V } } \sum _ { t \in \mathcal { T } } P \left ( s _ { v, t } ^ { + } + s _ { v, t } ^ { - } \right ) \quad ( 2 6 )$$

subject to:

$$& e _ { v, t } = e _ { v, t - 1 } + \Delta t \eta _ { v } c _ { v, t } - \widehat { \xi } _ { v, t } + s _ { v, t } ^ { + } - s _ { v, t } ^ { - }, \\ & \quad \forall t \in \mathcal { T }, v \in \mathcal { V }$$

$$c _ { v, t } \leq \overline { C } _ { v } \widehat { \alpha } _ { v, t }, \quad \forall v \in \mathcal { V }, t \in \mathcal { T }$$

$$\text{with-} \quad \underline { E } _ { v } \leq e _ { v, t } \leq \overline { E } _ { v }, \quad \forall v \in \mathcal { V }, t \in \mathcal { T }$$

$$e _ { v, N _ { T } } = e _ { v, 0 }, \quad \forall v \in \mathcal { V }$$

$$c _ { v, t }, s _ { v, t } ^ { + }, s _ { v, t } ^ { - } \geq 0, \ \forall v \in \mathcal { V }, t \in \mathcal { T },$$

where the set of decision variables Ξ DO = ( c v,t , e v,t , s + v,t , s -v,t ) and parameters ̂ α v,t and ̂ ξ v,t are expected values. The computation of these parameters is described in Section VI.

Fig. 1. Electricity prices (upper plot) and total number of EVs available (lower plot).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00787v2-figure-f72a4d48-6692-41a4-a342-e97beed6c93b.jpg)

## VI. CASE STUDY

We assume a residential aggregator with 100 EVs. For the sake of simplicity, the technical parameters associated with each EV are identical: The maximum charging rate is 7.4 kW, the round-trip efficiency is 0.95, the minimum and maximum energy rates are 10 and 51 kWh, in that order, and the energy rating per kilometer is 0.137 kWh/km [16].

Due to the lack of real-life data about the parameters associated with the driving patterns (availability profiles and energy required for transportation) of EVs, we resort to the National Household Travel Survey (NHTS) 2017 [17]. This data base contains synthetic data about vehicle trips, miles and duration per trip, etc. From this data base, we can extract the availability status by using the departure/arrival time periods for each daily trip. Specifically, we assume that the EV is available until it begins its first daily trip and after it returns from its last daily trip for each day of the year. Otherwise the EV is unavailable and thus it may be in a motion status. The energy required for transportation ξ v,t can be computed as the product of the travelled distance and energy rating per kilometer (i.e., 0.137 kWh/km). The lower plot of Fig. 1 shows the total number of EVs available in 15-min time periods for 29 days, whereas the Fig. 2 represents the energy consumption of the fleet of EVs per day. These data are used as realizations when validating the results in a real-time framework.

Fig. 1 also represents the 15-min electricity prices, which can be obtained from the ENTSO-e Transparency Platform [18] in year 2018 in Spain. We also assume that the penalty term P is set to 1000 e /kWh. We run daily simulations with 15-min time steps for one month (specifically 29 days).

Note that based on the data provided in figures 1 and 2, we compute the expected availability ̂ α v,t , energy consumption ̂ ξ v,t , and value of K v as follows: the expected values for day D , assuming that it is Monday, are the average over the data

Fig. 2. Total energy consumption of the fleet of EVs.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00787v2-figure-df2d2531-430b-4b63-82c1-563dcccd35cb.jpg)

on the previous four Mondays. In addition, we assume that the day-ahead electricity prices λ t on day D are the average over the last four days.

The simulations have been performed on a Windows 10 server with one CPU clocking at 2.8 GHz, 6 cores and 8 GB of RAM using CPLEX 12.6.3 [19] under Pyomo 5.2 [20]. The optimality gap is set to 0%.

## A. Day-ahead Operation of an Individual Electric Vehicle

We first analyze the day-ahead operation for two individual EVs in day 21. We can identify two types of EVs depending on their driving patterns' behaviour: one whose daily routine is highly predictable and thus the departure/arrival times may be narrowed down to a few time periods (i.e. in this case the availability can be fixed in most time periods throughout the time horizon); and another one whose daily routine is unpredictable and thus its uncertainty set is larger.

Fig. 3 shows the availability of two EVs, namely EV-A and EV-B, for both methods DO-EV and RO-EV. Note that, on the one hand, the expected availability for the DO-EV may be different from 0 or 1 because we average over historical data, and, on the other hand, the availability for the RO-EV is fixed in those time periods marked with a grey fill. The upper plot represents the availability of EV-A, whose pattern is predictable to some extent, i.e., the availability is known in the grey time intervals. The lower plot represents the EV-B, whose pattern is not predictable and thus, no grey zones are shown in the figure. Note also that parameter K v is set to 57 and 52, respectively, which means that at least 57 and 52 time periods the EV should be available. Figures 4 and 5 provide the charging schedule, day-ahead electricity prices, and energy state-of-charge, respectively.

In Fig. 4, we can observe that, with the RO-EV, the EV-A charges in time period 17, i.e., in the zone when the availability is certainly known and the price is significantly low, whereas the DO-EV schedules the charging of EV-A in time periods 23 and 24, i.e., a few time periods later when the prices are the lowest. In Fig. 5, we can see that the evolution of the energy state-of-charge is decreasing until the boundary condition is fulfilled. However, for the EV-B, we can observe a completely different behaviour since there is not any interval when the EV availability is known. Therefore, the RO-EV leads to



Fig. 3. Availability for the DO-EV and RO-EV for two EVs: (a) EV-A whose pattern is predictable to some extent, and (b) EV-B whose pattern is not predictable.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00787v2-figure-6923d06e-a81d-4fda-b798-64c8947062aa.jpg)



Fig. 4. Charging schedule for the DO-EV and RO-EV (left y-axis) and dayahead electricity prices (right y-axis) for (a) EV-A and (b) EV-B.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00787v2-figure-af0654c9-7a46-4a29-8595-abfaa9975e2f.jpg)

a solution in which the EV charges equally throughout the whole optimization horizon, as can be observed in Fig. 4.(b). Therefore, the energy state of charge is smoothly increasing during the first half of the day (see Fig. 5).

## B. Operation of the Fleet of Electric Vehicles

In order to fairly compare the performance of both methods, RO-EV and DO-EV, we solve the real-time operation of the aggregator for one month by fixing the energy bought from the day-ahead market (i.e., ∑ v ∈V c v,t ) and assuming the realized values of the parameters associated with each EV. The goal of the real-time problem is to minimize the deviations from the energy balance in the EVs' batteries for the fleet of EVs, which are penalized with 1000 e /kWh. As a result, we can



Fig. 5. Energy state-of-charge evolution for the DO-EV and RO-EV and for (a) EV-A and (b) EV-B.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00787v2-figure-39ec85d1-d07e-4ece-9abd-429f5e96006f.jpg)

TABLE I

RESULTS FOR DAY 21: DAY-AHEAD COST, PURCHASED POWER, AND ENERGY DEVIATIONS

| Metric     |   DO-EV |   RO-EV |
|------------|---------|---------|
| C DA ( e ) |    31   |    32.8 |
| P DA (kW)  |  1745.6 |  1851.2 |
| D RT (kWh) |   305.7 |   123.4 |

compare the following metrics: day-ahead cost C DA , dayahead purchased power P DA , and total deviations from the energy balance D RT in real-time, for each method. Both methods attain the optimal solution for one day in less than a minute.

Table I summarizes the results for day 21, which is the same day chosen in the previous section. The day-ahead operation cost of the aggregator increases by 5.8% when using the ROEV compared to the results from the DO-EV. In terms of total purchased energy, the RO-EV leads to an increase of 6.0% with respect to the deterministic optimization. Finally, this excess of purchased energy causes that there are less deviations with the RO-EV than the DO-EV; specifically, the energy deviations are reduced by approximately 60%.

## C. Monthly Operation of the Fleet of Electric Vehicles

Next, we compare the performance of both methods for one month. Fig. 6 represents the day-ahead purchased power and the real-time energy deviations for both methods and for each day of the month. Table II provides the monthly results in terms of day-ahead cost and total energy deviations. The optimal results for a month are attained after 340.4 s with the DO-EV and 619.8 s with the RO-EV.

As can be seen in Table II, although the day-ahead cost increases by 9.6% when using the robust optimization compared to the results from the DO-EV, the total energy deviations decreased to almost half. We can also observe that

Fig. 6. Monthly results: (a) Day-ahead power purchased per day and (b) real-time energy deviations per day.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00787v2-figure-212621a9-6853-43cd-8368-0cab2c430f77.jpg)

TABLE II

MONTHLY RESULTS: DAY-AHEAD COST AND ENERGY DEVIATIONS

| Metric     | Metric     |   DO-EV |   RO-EV |
|------------|------------|---------|---------|
| C DA ( e ) | C DA ( e ) |   586.6 |   643   |
| D RT (kWh) | Max.       |   305.7 |   190.5 |
| D RT (kWh) | Mean       |   156.9 |    82.9 |
| D RT (kWh) | Min.       |    77.5 |    24.3 |
| D RT (kWh) | Total      |  4548.9 |  2404.3 |

the maximum, mean, and minimum values are also reduced when using the RO-EV. The cost increase is because the ROEV schedules higher charging rates than the DO-EV to cope with the uncertainty on the availability and thus, the power purchased in the day-ahead market increases between 3.0% for day 20 and 13.0% for day 3 with respect to the power bought when using the DO-EV (see Fig. 6.(a)). As a result, the deviations of the energy balance of the vehicles' batteries are decreased all days except for day 8 and they may attain a reduction up to 83.0%, as observed for day 5 in Fig. 6.(b).

## VII. CONCLUSION

This paper provides a novel and computationally efficient model for the day-ahead operation of an aggregator of EVs in a residential district in which the energy bought is minimized subject to the technical characteristics of each EV. As a salient feature, the uncertainty on the availability of an EV is modelled via robust optimization so that the total energy required for transportation throughout the optimization horizon must be satisfied. In this problem, the uncertainty set leads to a mixedinteger linear program in the lower level and the duality-based approach can be applied thanks to the unimodularity of the system matrix and the convexity of the uncertainty set.

The case study reveals that the robust optimization model could lead to a substantial reduction of the daily deviations from the energy balance of the vehicles' batteries up to ap- proximately 80% compared to the results from a deterministic approach. Logically, this benefit comes at the expense of increasing the daily purchase costs in the day-ahead market around 15% compared to the deterministic counterpart. Further research will be devoted to refining the uncertainty set of the driving patterns, conducting a sensitivity analysis to the parameters defining that set and exploring the scalability of the proposed model to fleets comprising thousands of EVs. Moreover, we plan to extend our approach to EVs capable of discharging power to the grid.

## REFERENCES

- [1] 'Directive (EU) 2018/844 of the European Parliament and of the Council of 30 May 2018.' [Online]. Available: https://eur-lex.europa.eu/legalcontent/EN/TXT/PDF/?uri=CELEX:32018L0844&amp;from=ES.
- [2] T. Bunsen, P. Cazzola, M. Gorner, L. Paoli, S. Scheffer, R. Schuitmaker, J. Tattini, and J. Teter, 'Global ev outlook 2018: Towards cross-modal electrification,' 2018.
- [3] P. M. de Quevedo, G. Mu˜ noz-Delgado, and J. Contreras, 'Impact of electric vehicles on the expansion planning of distribution systems considering renewable energy, storage and charging stations,' IEEE Transactions on Smart Grid , vol. 10, no. 1, pp. 794-804, 2019.
- [4] L. P. Fernandez, T. G. San Rom´ an, R. Cossent, C. M. Domingo, and P. Frias, 'Assessment of the impact of plug-in electric vehicles on distribution networks,' IEEE Transactions on Power Systems , vol. 26, no. 1, pp. 206-213, 2011.
- [5] G. T. Heydt, 'The impact of electric vehicle deployment on load management straregies,' IEEE transactions on power apparatus and systems , no. 5, pp. 1253-1259, 1983.
- [6] R. Halvgaard, N. K. Poulsen, H. Madsen, J. B. Jørgensen, F. Marra, and D. E. M. Bondy, 'Electric vehicle charge planning using economic model predictive control,' in 2012 IEEE International Electric Vehicle Conference , pp. 1-6, IEEE, 2012.
- [7] E. B. Iversen, J. M. Morales, and H. Madsen, 'Optimal charging of an electric vehicle using a markov decision process,' Applied Energy , vol. 123, pp. 1-12, 2014.
- [8] M. R. Sarker, H. Pandˇ zi´ c, K. Sun, and M. A. Ortega-Vazquez, 'Optimal operation of aggregated electric vehicle charging stations coupled with energy storage,' IET Generation, Transmission &amp; Distribution , vol. 12, no. 5, pp. 1127-1136, 2017.
- [9] M. A. Ortega-Vazquez, 'Optimal scheduling of electric vehicle charging and vehicle-to-grid services at household level including battery degradation and price uncertainty,' IET Generation, Transmission &amp; Distribution , vol. 8, no. 6, pp. 1007-1016, 2014.
- [10] L. Baringo and R. S. Amaro, 'A stochastic robust optimization approach for the bidding strategy of an electric vehicle aggregator,' Electric Power Systems Research , vol. 146, pp. 362-370, 2017.
- [11] M. Gonz´ alez Vay´ a, L. Baringo, and G. Andersson, 'Integration of pevs into power markets: a bidding strategy for a fleet aggregator,' in Plug In Electric Vehicles in Smart Grids , pp. 233-260, Springer, 2015.
- [12] A. Street, F. Oliveira, and J. M. Arroyo, 'Contingency-constrained unit commitment with n -k security criterion: A robust optimization approach,' IEEE Transactions on Power Systems , vol. 26, no. 3, pp. 15811590, 2011.
- [13] A. Schrijver, Theory of linear and integer programming . John Wiley &amp; Sons, 1998.
- [14] D. G. Luenberger, Linear and nonlinear programming . USA, MA, Reading:Addison-Wesley, 1989.
- [15] C. A. Floudas, Nonlinear and mixed-integer optimization: fundamentals and applications . Oxford University Press, 1995.
- [16] 'Data Sheet Renault ZE ZOE.' [Online]. Available: http://electromovilidad.net/renault/ficha-tecnica-renault-zoe/, 2019.
- [17] 'National Household Travel Survey (NHTS) data.' [Online]. Available: https://nhts.ornl.gov/, 2019.
- [18] 'European Network of Transmission System Operators (ENTSO-e).' [Online]. Available: https://www.entsoe.eu/, 2019.
- [19] 'IBM ILOG CPLEX Optimisation Studio.' [Online]. Available: https://www.ibm.com/analytics/cplex-optimizer, 2019.
- [20] W. E. Hart, J. P. Watson, and D. L. Woodruff, 'Pyomo: Modeling and solving mathematical programs in Python,' Mathematical Programming Computation , vol. 3, no. 3, pp. 219-260, 2011.
Álvaro Porras Ricardo Fernández-Blanco Juan Miguel Morales Salvador Pineda
Day-ahead Operation of an Aggregator of Electric Vehicles via Optimization under Uncertainty
We pose the aggregator's problem as a bilevel model, where the upper level
minimizes the total operation costs of the fleet of EVs, while each lower level
minimizes the energy available to each vehicle for transportation given a
certain charging plan. Thanks to the totally unimodular character of the
constraint matrix in the lower-level problems, the model can be mathematically
recast as a computationally efficient mixed-integer program that delivers
charging schedules that are robust against the uncertain availability of the
EVs. Finally, we use synthetic data from the National Household Travel Survey
2017 to analyze the behavior of the EV aggregator from both economic and
technical viewpoints and compare it with the results from a deterministic
approach.