## Inverse Optimization with Kernel Regression: Application to the Power Forecasting and Bidding of a Fleet of Electric Vehicles

Ricardo Fern´ndez-Blanco a a, ∗ , Juan Miguel Morales a,b , Salvador Pineda a,c , ´ lvaro Porras A a a OASYS research group, University of Malaga, Malaga, Spain

- b Department of Applied Mathematics, University of Malaga, Malaga, Spain
- c Department of Electrical Engineering, University of Malaga, Malaga, Spain

## Abstract

This paper considers an aggregator of Electric Vehicles (EVs) who aims to learn the aggregate power of his/her fleet while also participating in the electricity market. The proposed approach is based on a data-driven inverse optimization (IO) method, which is highly nonlinear. To overcome such a caveat, we use a two-step estimation procedure which requires solving two convex programs. Both programs depend on penalty parameters that can be adjusted by using grid search. In addition, we propose the use of kernel regression to account for the nonlinear relationship between the behavior of the pool of EVs and the explanatory variables, i.e., the past electricity prices and EV fleet's driving patterns. Unlike any other forecasting method, the proposed IO framework also allows the aggregator to derive a bid/offer curve, i.e. the tuple of pricequantity to be submitted to the electricity market, according to the market rules. We show the benefits of the proposed method against the machine-learning techniques that are reported to exhibit the best forecasting performance for this application in the technical literature.

Keywords: Inverse optimization; kernel regression; forecasting; electric

∗ Corresponding author

Email addresses: Ricardo.FCarramolino@gmail.com (Ricardo Fern´ndez-Blanco), a

Juan.Morales@uma.es (Juan Miguel Morales), spinedamorente@gmail.com (Salvador Pineda), alvaroporras19@gmail.com ( ´ lvaro Porras) A

March 8, 2021

vehicles.

## Nomenclature

The main notation used throughout the text is stated below for quick reference. Other symbols are defined as required.

## Sets and Indices

- B Set of energy blocks, indexed by b .
- B c/d Set of energy blocks associated with the charging/discharging power, indexed by b .
- T Set of time periods, indexed by t and τ .
- Ω X Set of time periods belonging to the set X = { tr, v, test } where tr , v , test refer to the training, validation, and test set, in that order.

## Parameters

- E b,t , E b,t Width for the aggregate discharging/charging power block b in time period t [kW].
- H Feasibility penalty parameter.
- K t,τ Value of the kernel on two feature vectors at time periods t and τ .
- M Regularization hyper-parameter.
- N B Number of energy blocks.
- λ t Electricity price in time period t [ e /kWh].
- z t Vector of regressors in period t .
- γ Hyper-parameter related to the Gaussian kernel.

## Decision Variables

- m b,t Marginal utility of block b of the aggregate power in time period t [ e /kWh].
- p b,t Power in block b and time period t [kW].
- P , P t t Lower and upper bound for the aggregate power in time period t [kW].
- α , α t t Coefficient relative to the kernel regression of the lower/upper power bounds in period t ∈ Ω tr [kW].
- glyph[epsilon1] t Duality gap in time period t ∈ Ω tr [ e ].
- µ, µ Intercept for the lower/upper power bounds [kW].
- ν b Intercept for the marginal utility of block b [ e /kWh].
- ξ + t , ξ -t Slack variables associated with the lower power bound in time period t [kW].
- ξ + t , ξ -t Slack variables associated with the upper power bound in time period t [kW].
- ρ t Coefficient relative to the kernel regression of the marginal utility in time period t ∈ Ω tr [ e /kWh].

## 1. Introduction

According to the White Paper on transport of the European Commission (2011), one of the main goals to achieve a sustainable transport system is to halve the use of 'conventionally fuelled' cars in urban transport by 2030; phase them out in cities by 2050; achieve essentially CO 2 -free city logistics in major urban centres by 2030 . This will spur the use of electric vehicles (EVs) across Europe (European Commission, 2011). Although nowadays the penetration of EVs in the European market is slow albeit steady, the estimated electricity demand from all EVs worldwide was 54 TWh in 2017 (Bunsen et al., 2018).

Thus, the growing electrification of the road transport will impact the power system operation and planning of the future and, as a consequence, new actors and facilities will come into play, e.g. aggregator agents (Bandpey &amp; Firouzjah, 2018), or battery swap stations (Yang &amp; Sun, 2015).

Within the context of the restructured power industry, the aggregator agents face several challenges: (i) the forecast of the charging power of the fleet of EVs in the short-term, and (ii) the determination of a bid curve to participate in the electricity market to maximize their profits when the fleet of EVs is large enough. Nowadays, the EVs may be prepared with bi-directional vehicle-to-grid (V2G) capabilities, which means that the EVs can extract power from and inject power into the electrical grid while parked (Kempton &amp; Tomi´, 2005). c This is possible as long as the EVs are equipped with the necessary smart meteringand-control infrastructure as well as a suitable connection to the electrical grid. In this case, the aggregator will also need to forecast the EV-fleet discharging power.

Short-term load forecasting is widely applied in the power sector to predict the electricity demand (and price) for different granularity levels (Shahidehpour et al., 2003). In the last years, EV charging load forecasting tools have been proposed in the technical literature by means of ARIMA-based models (Amini et al., 2016; Korolko et al., 2015); machine-learning techniques (Majidpour et al., 2016; Sun et al., 2016; Xydas et al., 2013), such as support vector regression; or big data technologies (Arias &amp; Bae, 2016). All these papers neglected the bi-directional V2G capabilities of the EVs. Moreover, the above methodologies aimed to provide a single-purpose application, i.e., the forecasting of the charging power of either an EV or a fleet of EVs. Instead, we propose here a multi-purpose application for the aggregator of EVs in order to not only forecast the EV-fleet power, but also to derive a bid/offer curve according to the rules of the electricity market, e.g. see OMIE (2019).

In this paper, we apply inverse optimization (IO) to forecast the EV-fleet power while deriving a bid/offer curve. The goal of an IO problem is to infer the optimization model parameters given a set of observed decision variables or mea-

surements collected by an observer. For instance, Zhang &amp; Xu (2010) applied IO for linearly-constrained convex problems in the industrial and managerial areas but its application was limited to single observed decisions. Aswani et al. (2018) proposed a statistically consistent methodology for IO when the measurements of the optimal decisions of a convex optimization problem are noisy. In a more general context, when the observer has imperfect information, Esfahani et al. (2018) devised a distributionally robust inverse optimization problem. IO has also been applied for equilibrium problems (Bertsimas et al., 2015), multiobjective convex optimization (Roland et al., 2016), or robust optimization (Chan &amp; Kaw, 2019). However, few papers have implemented IO in the field of power systems (Saez-Gallego et al., 2016; Saez-Gallego &amp; Morales, 2017; Lu et al., 2018; Ruiz et al., 2013; Zhou et al., 2010).

Zhou et al. (2010) applied IO in the context of generation expansion planning to find an effective incentive policy; Ruiz et al. (2013) estimated rival marginal offer prices for a strategic producer in a network-constrained dayahead market by using IO; Saez-Gallego et al. (2016) prescribed an IO approach by using bi-level programming to infer the market bid parameters of a pool of price-responsive consumers; in Saez-Gallego &amp; Morales (2017), a novel IO approach was devised to statistically estimate the aggregate load of a pool of price-responsive buildings in the short-term; and, finally, Lu et al. (2018) applied IO to estimate the demand response characteristics of price-responsive consumers, as similarly done in Saez-Gallego et al. (2016). Unlike existing works (Saez-Gallego et al., 2016; Saez-Gallego &amp; Morales, 2017; Lu et al., 2018; Ruiz et al., 2013; Zhou et al., 2010), we address the EV-fleet power forecasting with an IO approach in which the prediction tool accounts for two distinctive features: (i) the pool of EVs may be equipped with V2G capabilities, and (ii) there may exist a strong nonlinear relationship between the EV-fleet power and the explanatory variables, namely past EVs' charging/discharging patterns and past electricity prices. To capture these nonlinear relations, we endogenously introduce kernels into the proposed IO approach.

Kernels are widespread in the literature on machine learning, as can be

seen in Hofmann et al. (2008); Trevor et al. (2009); Ben´ ıtez-Pe˜a et al. (2019), n just to name a few; and, in power systems, they were mainly used to predict electricity prices (Dudek, 2018; Kekatos et al., 2013, 2014). Kekatos et al. (2013) applied a kernel regression to forecast the electricity prices from the Midwest Independent System Operator day-ahead market in which the kernel itself is constructed by the product of three kernels: one for vectorial data and other two to account for non-vectorial data such as time and nodal information. This approach was generalized to low-rank kernel-based learning models in Kekatos et al. (2014). Finally, Dudek (2018) devised a probabilistic forecast method built on the Nadaraya-Watson estimator to predict the electricity prices from the Polish balancing and day-ahead markets.

The contributions of this paper are threefold:

- · From a modeling perspective, we provide an IO framework to forecast the aggregate power of a fleet of EVs with V2G capabilities. In addition, the outcome of this framework may be used to bid/offer in the electricity market by using the estimated price-quantity tuples. To the best of the authors' knowledge, this is the first time in the technical literature that IO has been used to forecast the aggregate power of a price-responsive EVs' aggregator and to derive a suitable bid/offer curve for such an aggregator.
- · We approximate the solution of the generalized IO problem by using a data-driven two-step estimation procedure. This procedure requires solving two different convex programming problems, which makes the process of building the forecasting model computationally affordable. As a salient feature of this work, a kernel is endogenously incorporated into the regression functions.
- · We thoroughly analyze the performance of the proposed methodology by using real-life data based on the latest National Household Travel Survey (NHTS, 2019) and we compare the results against those provided by two machine-learning techniques, namely support vector regression and

kernel-ridge regression. The former has been reported to exhibit the best forecasting performance for the present application in the technical literature (Xydas et al., 2013; Sun et al., 2019).

The rest of the document is organized as follows: Section 2 provides the IO methodology; Section 3 gives a general overview on the comparison methodologies; in Section 4, we analyze a case study for a residential aggregator of EVs; conclusions are duly drawn in Section 5; and, finally, Appendix A presents a mixed-integer linear programming problem to generate synthetic data on the behavior of an EV fleet.

## 2. Inverse Optimization Methodology

To put the problem in context, we aim to forecast or learn the EV-fleet power p t (also known as aggregate power) in time period t of a price-responsive aggregator, who is also interested in deriving a bid/offer curve to be submitted to the electricity market. The participants of the electricity market, namely consumers and producers, must submit a bid/offer curve consisting of blocks of energy and price. For the consumers, the bid curve should be monotonically nonincreasing, whereas, for the producers, the offer curve should be monotonically non-decreasing, e.g. see OMIE (2019). We assume a rational aggregator, which means that the market strategy of the EV fleet fundamentally relies on arbitrage , by behaving as a consumer when the electricity price is low and, on the contrary, by acting as a producer when the price is high.

In order to predict the EV-fleet aggregate power and to derive a bid/offer curve, the aggregator may use past observed data, which are denoted as explanatory variables, features or regressors. As one should expect for a price-responsive aggregator, the regressors in time period t can be the lagged electricity price λ ′ t -l or aggregate power p ′ t -l , ∀ l = 1 2 , , ... . In addition, past EV driving patterns, meteorological data, or categorical data (e.g., time information) can also be used for forecasting purposes.

Within this context, we first introduce the proposed forecasting 1 model in Section 2.1. Subsequently, Section 2.2 explains how we can account for past information. Finally, Section 2.3 thoroughly describes the two-step procedure to estimate the required parameters of the forecasting model.

## 2.1. Forward Model

The key idea of this work is to forecast the EV-fleet power by using a simple optimization (linear programming) model which may, to some extent, mimic its real behavior. In addition, unlike other forecasting techniques, this model is able to derive a bid/offer curve, as imposed by rules of electricity markets. Therefore, the formulation of the forward model that, we assume, represents the aggregate response of an EV fleet to the electricity prices at time period t , is mathematically expressed as:

$$\max _ { p _ { b, t } } \ \sum _ { b \in \mathcal { B } } p _ { b, t } \left ( m _ { b, t } - \lambda _ { t } \right )$$

subject to:

$$\underline { P } _ { t } \leq \sum _ { b \in \mathcal { B } } p _ { b, t } \leq \overline { P } _ { t } \colon ( \underline { \beta } _ { t }, \overline { \beta } _ { t } )$$

$$0 \leq p _ { b, t } \leq \overline { E } _ { b, t } \colon ( \underline { \phi } ^ { c } _ { b, t }, \overline { \phi } ^ { c } _ { b, t } ), \ \forall b \in \mathcal { B } ^ { c }$$

$$\underline { E } _ { b, t } \leq p _ { b, t } \leq 0 \, \colon ( \underline { \phi } _ { b, t } ^ { d }, \overline { \phi } _ { b, t } ^ { d } ), \quad \forall b \in \mathcal { B } ^ { d },$$

where dual variables are represented in parentheses after a colon in the respective constraints. For the sake of unit consistency, hourly time periods are considered.

The reconstruction problem (1) aims to maximize the welfare of the EV aggregator, as given by the objective function (1a). This objective function is made up of the EV fleet's surplus, which is related to the aggregate charging and discharging power. The aggregate power is positive when the EVs' aggregator is charging, i.e., it behaves as a consumer. Otherwise, the aggregate power takes on negative values when the aggregator is discharging, i.e., it acts as a

1 This problem is also known as forward or reconstruction problem in the IO jargon.

producer. We assume step-wise offer/bid price functions as depicted in Fig. 1, as is customary in real-world electricity markets, e.g. see OMIE (2019). Constraints (1b) represent the lower and upper bounds on the aggregate power. Constraints (1c) impose the lower and upper bound on each block b within the set B c of charging power blocks. Since the charging power is assumed to be non-negative, then p b,t is bounded between 0 and a positive power bound E b,t . Likewise, constraints (1d) impose the lower and upper bound on each block b within the set B d of the discharging power blocks. We assume that the discharging power is non-positive and thus p b,t is bounded between a negative power bound E b,t and 0. Note that the total power p t = ∑ b p b,t .

Figure 1: Three-block stepwise offer (bid) price function of the EVs' aggregator. In this example, the offer (bid) price function is represented to the left (right) of the y-axis, and the sets B d = {-3 , -2 , - } 1 and B c = { 1 2 3 , , } .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00399v4-figure-48392f9a-3370-4ee8-8922-4436e2b21bd7.jpg)

As previously stated, we want to anticipate the EV-fleet power response by solving (1). However, to this end, the set of parameters Φ = { E b,t , E b,t , m b,t , P t , P t } needs to be estimated since they are a priori unknown. These parameters should be functions of time and of any regressor that the forecaster may consider meaningful and explanatory of the EV-fleet's operational behavior and, therefore, are to be inferred from past observations of the aggregate power p ′ t , the electricity price λ ′ t , and the regressors that are eventually considered. This fact gives rise to a generalized IO problem, which is highly nonlinear and nonconvex. This problem can be naturally formulated as a bilevel optimization problem, which may be computationally nonviable when moderately increasing the sample size. To deal with such complexity, we apply a methodology that

builds on the one first proposed in Saez-Gallego &amp; Morales (2017). In that paper, however, the regression function is linear in their features and may be limited to capture nonlinear relations between the EV-fleet power and the regressors. To circumvent such a caveat, and as one of the salient features of this work, we incorporate kernels into the regression functions. Furthermore, the forward model we propose, i.e. problem (1), allows for power intakes and outputs, unlike the one used in Saez-Gallego &amp; Morales (2017). This extra dose of model flexibility is critical to capture the behavior of an EV fleet with V2G capabilities since the aggregator power may be positive when the net power comes from the grid (i.e. the aggregator acts as a consumer) or negative when the net power flows into the grid (i.e. the aggregator behaves as a producer). Therefore, the forecasting model is tailored to account for this dual operational mode by introducing differentiated marginal utility blocks both for charging and discharging.

## 2.2. Accounting for Past Information: Kernels

In the realm of machine learning, the kernel functions are rather popular in learning algorithms (Hofmann et al., 2008) since they are able to capture nonlinear relationships between the dependent and the explanatory variables. Unlike in Saez-Gallego &amp; Morales (2017), where affine functions were used to model the dependence of the parameters of the forward model (1) on the regressors, we propose the use of kernel regressions to estimate P t , P t , and m b,t :

$$\underline { P } _ { t } = \underline { \mu } + \sum _ { \tau \in \Omega ^ { t r } } \underline { \alpha } _ { \tau } K _ { t, \tau }, \quad \forall t \in \mathcal { T }$$

$$\overline { P } _ { t } = \overline { \mu } + \sum _ { \tau \in \Omega ^ { t r } } \overline { \alpha } _ { \tau } K _ { t, \tau }, \ \forall t \in \mathcal { T }$$

$$m _ { b, t } = \nu _ { b } + \sum _ { \tau \in \Omega ^ { t r } } \rho _ { \tau } K _ { t, \tau }, \quad \forall t \in \mathcal { T }.$$

Many kernel functions can be used: polynomial, hyperbolic tangent, Gaussian, among others. For the sake of illustration purposes, the Gaussian kernel

(Trevor et al., 2009) can be defined as follows:

$$K _ { t, \tau } = K \left ( z _ { t }, z _ { \tau } \right ) = e ^ { - \gamma \| z _ { t } - z _ { \tau } \| _ { 2 } ^ { 2 } }, \ \forall t \in \mathcal { T }, \tau \in \Omega ^ { t r },$$

wherein γ is a scale parameter inversely proportional to the variance of the Gaussian function; and ‖ z t -z τ ‖ 2 2 is the squared Euclidean distance between two feature vectors at time periods t and τ . Thus, the Gaussian kernel can be interpreted as a similarity measure between two time periods, i.e., if the two feature vectors are identical z t = z τ , then the value of K t,τ = 1, otherwise its value ranges in the interval (0 , 1].

As previously mentioned, meaningful or explanatory features should be used in the kernel regression function for adequately inferring the estimates. In the proposed IO methodology, the power bounds P t and P t are key to capturing the price-responsiveness of the aggregator since they determine the width and the range of the step-wise price-response function of the EV fleet. For instance, if the electricity price is high, one should expect that the EV fleet will behave as a producer (in V2G mode) and therefore, the power bounds would take on negative values, this way producing a step-wise offer curve displaced towards negative power values (i.e., discharge). On the contrary, if the price is low, one should expect the opposite: the EV fleet would act as a consumer, with the power bounds taking positive values and defining a step-wise bidding curve displaced towards positive power values (i.e., charge). On the other hand, the marginal utilities aim to capture the price-sensitivity of the EVs aggregate power. Therefore, by making both the power bounds and the marginal utilities dependent on past prices along with past values of the aggregate EV-fleet power and/or other external factors, we can capture the changes in the EV-fleet power due to price variations over time.

Illustrative example. Let us assume that z t comprises only one regressor, e.g. the electricity price in the previous time period, i.e., z t = λ t -1 . Thus, Fig. 2 provides the values of the kernel for each time period t of a day with respect to the second time period τ = 2, i.e., K t,τ =2 , for different values of parameter γ . Moreover, the values of the regressor z t for the 24 hours are shown in the figure.

We can observe that high values of γ lead to kernel values equal to 1 just when the two regressors are very close to each other (e.g., see time periods 21-23 for γ = 1); conversely, low values of γ lead to kernel values equal to 1 even when the regressors are very different from each other (e.g., see values for all time periods when γ = 0 001). . Therefore, we need to carefully tune the hyper-parameter γ , as described in Section 2.3.3.

Figure 2: Values of the Gaussian kernel for each time period t of a day with respect to period τ = 2, i.e., K t,τ =2 , for different values of the parameter γ in the right y-axis and the corresponding regressor values in the left y-axis.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00399v4-figure-6e277780-d6a8-46cd-87d6-1a6e51b86613.jpg)

## 2.3. Two-step Estimation Procedure

The thrust of this work is the estimation of the set of parameters Φ = { E b,t , E b,t , m b,t , P t , P t } and the corresponding coefficient estimates µ, α , µ, α , t t ν , ρ b t of the regression functions described in (2)-(4). To do that, we can use bilevel optimization, however, as mentioned previously, it may lead to a prohibitive computational burden when moderately increasing the sample size. Therefore, we resort to a two-step procedure based on two convex programming problems: (i) the feasibility problem , which is devoted to estimating all parameters that determine the feasibility of the observed EV-fleet power values in the forward problem (1) (i.e., the power bounds), and (ii) the optimality problem , which estimates the marginal utility of the EVs' aggregator, i.e., the parameters of problem (1) that are related to the optimality of the observed power values. The key idea of the feasibility problem is to shape the power bounds P t and

P t so that a certain percentage H of the observed EV-fleet power values are feasible for the forward problem (1). Note that the width for the aggregate power blocks E b,t and E b,t can be easily computed from the estimated power bounds by assuming that the energy blocks are all of same length. Conversely, the optimality problem estimates the marginal utilities m b,t driven by the minimization of the duality gap of the forward problem once the power bounds are fixed. Its aim is thus to make the observed EV-fleet power values as optimal as possible for problem (1) (recall that we use (1) as the forward model). It should be noted that the pair ( m b,t , E b,t ) for all blocks constitutes the bid curve of the aggregator at time period t . Likewise, the pair ( m b,t , -E b,t ) for all blocks constitutes the offer curve of the aggregator at time period t . In practice, those curves may be submitted to the market operator, who is the entity responsible for the financial management of electricity markets, e.g. see OMIE (2019).

## 2.3.1. Feasibility Problem

Given a fixed value of control parameter H ∈ [0 , 1), this problem can be formulated as:

$$\min _ { \Xi ^ { f _ { p } } } \sum _ { t \in \Omega ^ { t r } } H \left ( \bar { \xi } _ { t } ^ { - } + \xi _ { t } ^ { - } \right ) + \sum _ { t \in \Omega ^ { t r } } \left ( 1 - H \right ) \left ( \bar { \xi } _ { t } ^ { + } + \xi _ { t } ^ { + } \right ) \quad \ \ ( 6 a )$$

subject to:

$$\overline { P } _ { t } - p _ { t } ^ { \prime } = \overline { \xi } _ { t } ^ { + } - \overline { \xi } _ { t } ^ { - }, \quad \forall t \in \Omega ^ { t r }$$

$$p _ { t } ^ { \prime } - \underline { P } _ { t } = \xi _ { t } ^ { + } - \xi _ { t } ^ { - }, \quad \forall t \in \Omega ^ { t r }$$

$$\overline { P } _ { t } \geq \underline { P } _ { t }, \quad \forall t \in \Omega ^ { t r }$$

$$\text{Constraints } ( 2 ) - ( 3 )$$

$$\bar { \xi } _ { t } ^ { + }, \bar { \xi } _ { t } ^ { + }, \bar { \xi } _ { t } ^ { - }, \xi _ { t } ^ { - } \geq 0, \quad \forall t \in \Omega ^ { t r },$$

where the set of variables to be optimized is Ξ fp = { P , P t t , ξ + t , ξ + t , ξ -t , ξ -t , µ, µ, α , α t t } . Note that problem (6) is a convex program.

The objective function (6a) minimizes the sum of feasibility and infeasibility slack variables associated with the power bounds. Constraints (6b)-(6c) are the power bound constraints with the feasibility and infeasibility slack variables,

where p ′ t is the observed EV-fleet power value at time period t . Constraints (6d) ensure that the upper bound of the aggregate power is greater than its respective lower bound. Constraints (6e) impose kernel regression functions for the power bounds wherein the coefficients to be estimated are µ , µ , α t , α t . Finally, constraints (6f) declare the variables ξ + t , ξ + t , ξ -t , ξ -t as non-negative. Importantly, the higher the value of H , the wider the power bounds delivered by (6) and, therefore, the more price-responsive the EV fleet is expected to be.

The use of kernels increases the flexibility of the regression function when increasing the size of the training set. However, it also tends to over-fitting. To control the risk of over-fitting, a regularization parameter M ∈ [0 , 1] is used to factor in the sum of the squared values of the coefficient estimates α t and α t , similarly to what is typically done in kernel-ridge regression (Trevor et al., 2009). Thus, the objective function (6a) should be replaced with:

$$\min _ { \Xi ^ { f _ { r } } } M \sum _ { t \in \Omega ^ { t r } } & \left ( \underline { \alpha } _ { t } ^ { 2 } + \overline { \alpha } _ { t } ^ { 2 } \right ) \\ & + ( 1 - M ) \left [ \sum _ { t \in \Omega ^ { t r } } H \left ( \overline { \xi } _ { t } ^ { - } + \xi _ { t } ^ { - } \right ) + \sum _ { t \in \Omega ^ { t r } } ( 1 - H ) \left ( \overline { \xi } _ { t } ^ { + } + \xi _ { t } ^ { + } \right ) \right ].$$

Both hyper-parameters M and H in the objective function and parameter γ of the kernel regression function must be adequately adjusted to modulate the power bounds to the observed EV-fleet power values so that the out-of-sample forecasting error is minimized.

## 2.3.2. Optimality Problem

Once the power bounds (i.e., ̂ P t , ̂ P t ) are estimated from (6), we can compute the power block limits ̂ E b,t , ∀ b ∈ B c and ̂ E b,t , ∀ b ∈ B d based on the assignments described in Table 1. The optimality problem can then be derived by using

results from duality theory of linear programming and it can be formulated as:

$$\mathcal { N }$$

$$& \min _ { \Xi ^ { o p } } \sum _ { t \in \Omega ^ { t r } } \epsilon _ { t } \\ & \widehat { P } _ { t } \overline { \beta } _ { t } - \widehat { P } _ { t } \overline { \beta } _ { t } + \sum _ { b \in \mathcal { B } ^ { c } } \widehat { \widehat { E } } _ { b, t } \overline { \phi } _ { b, t } ^ { c } - \sum _ { b \in \mathcal { B } ^ { d } } \widehat { E } _ { b, t } \overline { \phi } _ { b, t } ^ { d } - \epsilon _ { t } = \\ & \quad \cdots.$$

$$\sum _ { b \in \mathcal { B } } p _ { b, t } ^ { \prime } \left ( m _ { b, t } - \lambda _ { t } \right ), \forall t \in \Omega ^ { t r } \ \ ( 8 b )$$

$$- \, \underline { \phi } _ { b, t } ^ { c } + \overline { \phi } _ { b, t } ^ { c } - \underline { \beta } _ { t } + \overline { \beta } _ { t } = m _ { b, t } - \lambda _ { t }, \quad \forall b \in \mathcal { B } ^ { c }, t \in \Omega ^ { t r }$$

$$- \, \underline { \phi } _ { b, t } ^ { d } + \overline { \phi } _ { b, t } ^ { d } - \underline { \beta } _ { t } + \overline { \beta } _ { t } = m _ { b, t } - \lambda _ { t }, \quad \forall b \in \mathcal { B } ^ { d }, t \in \Omega ^ { t r }$$

$$\text{Constraints} \left ( 4 \right )$$

$$\nu _ { b } \geq \nu _ { b + 1 }, \ \forall b \in \mathcal { B } \, \wr \, \{ b = N _ { B } \}$$

$$\underline { \beta } _ { t }, \overline { \beta } _ { t }, \underline { \phi } _ { b, t } ^ { c }, \overline { \phi } _ { b, t } ^ { c }, \underline { \phi } _ { b, t } ^ { d }, \overline { \phi } _ { b, t } ^ { d } \geq 0, \quad \forall t \in \Omega ^ { t r },$$

where the set of decision variables is Ξ op = { m b,t , glyph[epsilon1] t , β t , β t , φ c b,t , φ c b,t , φ d b,t , φ d b,t , ν , ρ b t } . Note that problem (8) is a convex program.

|        | Table 1:         | Value of ̂ E b,t , ∀ b ∈ B c and ̂ E ̂ P t ≥ ̂ P t ≥ 0 ≤ ̂   | b,t , ∀ b ∈ ̂ P t P t ≤ 0   | B d ̂ P t ≥ 0 ≥ ̂ P t   |
|--------|------------------|---------------------------------------------------------|----------------------------|-----------------------|
| ̂ E b,t | b = 1            | ̂ P t                                                    | 0                          | ̂ P t /N B             |
| ̂ E b,t | b ∈ B c \ { 1 }  | ( ̂ P t - ̂ P t ) N B - 1                                 | 0                          | ̂ P t /N B             |
| ̂ E b,t | b = - 1          | 0                                                       | ̂ P t                       | ̂ P t /N B             |
| ̂ E b,t | b ∈ B d \ {- 1 } | 0                                                       | ( ̂ P t - ̂ P t ) N B - 1    | ̂ P t /N B             |

The objective function (8a) minimizes the sum of the duality gaps of problem (1). Constraints (8b) is the relaxed equality constraint associated with the strong duality theorem. Constraints (8c)-(8d), (8g) are the dual feasibility constraints. Constraints (8e) impose a kernel regression function, with ν b and ρ t as the coefficients to be estimated, in order to relate the marginal utilities and the regressors. Finally, constraints (8f) set the marginal utilities to be monotonically non-increasing, as imposed by rules in electricity markets (OMIE, 2019).

## 2.3.3. Statistical Computation of Hyper-Parameters

The main goal of this work is to learn the EV-fleet power for each period t ∈ Ω test with the forward model (1), which relies on the knowledge of a series of parameters, i.e., the power bounds and the marginal utilities. Those parameters are estimated with the models described in Sections 2.3.1 and 2.3.2, whose outcome depend on the value of three hyper-parameters: H M , , and γ . Their optimal values are computed by using a grid search technique. We recursively solve problems (6) and (8) for the training set Ω tr ; and we then solve the forward problem (1) over the validation set Ω v by using the estimated parameters Φ = { E b,t , E b,t , m b,t , P t , P t } as well as the electricity price at time period t ∈ Ω . v Thus, we set as the optimal values of the hyper-parameters those that lead to the least out-of-sample forecasting error in Ω v .

## 3. Comparison Methodologies

We compare the performance of the proposed kernel-based IO approach, hereinafter referred to as kio , against (i) the state-of-the-art model to forecast the EV-fleet power, namely kernelized support vector regression ( svr ), (ii) a kernel-ridge regression model ( krr ), (iii) an IO approach with linear kernels ( lio ), and (iv) persistence or naive models. Note that we use a Gaussian kernel in the regression functions of the feasibility problem and a linear kernel in the regression function of the optimality problem , as this combination exhibited the best trade-off between forecasting performance and simplicity in our numerical experiments.

Regarding the svr and krr , we respectively use the epsilonsvr and the kernelridge regression models implemented in the scikit-learn library (Pedregosa et al., 2011) under the Python programming language. The interested reader is referred to Smola &amp; Sch¨lkopf (2004) for a detailed description of the o svr . For the sake of comparison, we also use the Gaussian kernel and we tune the corresponding hyper-parameters via grid search. Specifically, we tune the cost of constraints violation C and the parameter associated with the kernel γ for svr ; and the penalty parameter δ and the γ parameter for krr .

Regarding the naive models, we use three different ones since the EV-fleet power may experience seasonal patterns: h-naive , d-naive , and w-naive , in which the forecast value of the aggregate power at time t is equal to the observed value at time t -1, t -24, and t -168, in that order. Note that the forecast error of the naive models provides insight into the difficulty of prediction.

The performance of the methods is compared with two metrics: the mean absolute error (MAE) and the root mean square error (RMSE) on the test set.

## 4. Case Study

We first describe the data used for the case study in Section 4.1. Subsequently, we comprehensively analyze the results from the proposed approach for three cases of charging behavior without enabling the V2G capabilities in Section 4.2. Finally, Section 4.3 presents the results for two cases when the electric vehicles are integrated with V2G services.

## 4.1. EV-fleet Data

For learning purposes, we would be only interested in the time series of electricity prices, the aggregate power of an EV fleet, and the total number of available vehicles to charge or discharge. However, to our knowledge, there is no real-life data available about an EVs' aggregator. Thus, we resort to the formulation of an optimization problem to simulate the behavior of such an EV fleet. The interested reader is referred to Appendix A for a detailed description of this simulator.

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 (Renault, 2019). 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, 2019). 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 electricity prices are obtained from the ENTSO-e Transparency Platform (ENTSO-e, 2019) for the period comprising January 9 th till February 19 th in Spain. We also assume that the load shedding cost C P = 1000 e /kWh. We run daily simulations with 15-min time steps to build a synthetic database for a pool of EVs.

The simulations have been performed on a Linux-based server with one CPU clocking at 2.6 GHz and 2 GB of RAM using CPLEX 12.6.3 (CPLEX, 2019) under Pyomo 5.2 (Hart et al., 2011). Optimality gap is set to 0%. The input data files for reproducing the results have been shared with the scientific community in https://github.com/groupoasys/Aggregated-EV-data .

## 4.2. Forecast Results without Enabling V2G Capabilities

We assume that EVs do not enable their V2G capabilities (i.e. B d v = 0 in the model (A.1)-(A.9) in the Appendix A) and we compare the results for three cases: (i) a case in which the EVs satisfy their energy needs by using a naive charging; (ii) a case in which the charging is highly synchronized, which occurs when C S is set to 0 in (A.1)-(A.9); and (iii) a case in which the charging synchronization is avoided, which we attain by setting C S = 520 e /MWh . 2 Those cases are respectively denoted as naive-ch , sync , and non-sync . Note that, in the former case, i.e. naive-ch , each EV will be charged to its required maximum energy as soon as it is available, thus neglecting the dependence of the charging power on the price; whereas, the latter cases sync and non-sync are driven by the cost minimization of the EVs' aggregator wherein the electricity prices are accounted for. As an example, Fig. 3 shows the EV-fleet charging

glyph[negationslash]

power of a certain day for the three cases along with the electricity prices. As can be seen, the choice of C S = 0 is a simple albeit convenient way to avoid the undesirable charging synchronization by smoothing the aggregate power. In addition, we can observe that the charging pattern of the naive-ch case is independent of the prices.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00399v4-figure-207c390d-aef3-4745-92a2-645297bc38b6.jpg)

Case naive-ch

Case sync

Case non-sync

Price

Figure 3: Charging power for cases naive-ch , sync , and non-sync in the left y-axis and the corresponding electricity prices in the right y-axis.



{∐}

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00399v4-figure-1f876dc7-7f2c-47b7-866b-86a7d603d0d0.jpg)

{

{

Figure 4: Power versus price for cases (a) naive-ch , (b) sync , and (c) non-sync .

The sizes of the training, validation, and test sets are 672 h, 168 h, and 168 h, in that order. Fig. 4 represents the hourly electricity price versus the corresponding charging power for all periods of the Ω tr for the cases mentioned above. As can be seen, the aggregate power of the non-sync case depends linearly on the price, unlike the naive-ch and sync cases. For the naive-ch case, we consider 17 regressors, namely the charging power and the total number of EVs available for the six periods previous to time t , i.e., p t -l and ∑ v ς v,t -l , ∀ l = 1 ... 6, and 5 binary-valued categorical variables to indicate the hour of the day. For the cases sync and non-sync , we consider 12 regressors, namely the

electricity price and the charging power for the six periods previous to time t , i.e., λ t -l and p t -l , ∀ l = 1 ... 6. We also assume six energy blocks in total. Finally, hyper-parameter H ranges in the interval [0 5 . , 1 0) with 0.01 steps, . M ranges in the interval [0 0001 . , 0 0024] with 0.0001 steps, and . γ = 0 1 0 01 . For the case { . , . } sync , the proposed approach kio takes on average 12.6 s, 2.6 s, and 31.3 s to run each feasibility problem, optimality problem, and all the forward problems for the Ω , in that order. v The computing times are of the same order of magnitude for the other cases. It should be noted that those computing times would be even suitable for an hour-ahead forecasting if the grid search technique were parallelized.

The optimal hyper-parameters for all models and cases are given in Table 2. The information given in this table is quite valuable and we can make two main remarks. First, cases sync and non-sync are price-driven and thus their optimal values of parameter H ∗ are very high (0 82 and 0 94 respectively) compared to . . the optimal value ( H ∗ = 0 64) for the case . naive-ch , which is insensitive to the prices. In other words, the power bounds for the former cases are wider than for the latter one. Therefore, the optimality problem , which is used to estimate the marginal utility, plays a major role to learn the aggregate response of the EV fleet for the price-driven cases. This is expected as the marginal utilities encode the impact of the current electricity price on the aggregate power of the EV fleet. Second, it should be noted that the optimal values of H ∗ for the models kio and lio are quite similar, except for the case naive-ch , for which lio is unable to identify the insensitiveness of the aggregate power to the price.

The error metrics of the test set for all models are compared in Table 3 for the three cases. In the naive-ch case, the least RMSE is obtained with the proposed model kio with an error reduction of 4.4% and 17.3% compared to krr and svr . In the sync case, the proposed model kio achieves 28.3% reduction in RMSE and 15.3% reduction in MAE compared to the w-naive , which provides the best performance among the naive models. As expected, we can also observe that the kio outperforms lio by reducing RMSE and MAE by 40.6% and 42.2% since kio is able to capture the nonlinear relations between the EV-fleet power

Table 2: Optimal Values of the Hyper-Parameters

| Case     | kio                                     | krr                      | svr                  | lio          |
|----------|-----------------------------------------|--------------------------|----------------------|--------------|
| naive-ch | H ∗ = 0 . 64 M ∗ = 0 . 0002 γ ∗ = 0 . 1 | δ ∗ = 0 . 01 γ ∗ = 0 . 1 | C ∗ = 100 γ = 0 . 01 | H ∗ = 0 . 91 |
| sync     | H ∗ = 0 . 82 M ∗ = 0 . 0001 γ ∗ = 0 . 1 | δ ∗ = 0 . 1 γ ∗ = 0 . 1  | C ∗ = 10 γ ∗ = 0 . 1 | H ∗ = 0 . 89 |
| non-sync | H ∗ = 0 . 94 M ∗ = 0 . 002 γ ∗ = 0 . 01 | δ ∗ = 0 . 1 γ ∗ = 0 . 1  | C ∗ = 1 γ = 0 . 1    | H ∗ = 0 . 94 |

Table 3: Error Metrics - Cases without V2G Services (kW)

| Model   | naive-ch   | naive-ch   | sync   | sync   | non-sync   | non-sync   |
|---------|------------|------------|--------|--------|------------|------------|
| Model   | RMSE       | MAE        | RMSE   | MAE    | RMSE       | MAE        |
| kio     | 8.6        | 3.7        | 35.2   | 13.3   | 5.5        | 3.8        |
| krr     | 9.0        | 3.5        | 35.5   | 15.7   | 7.4        | 5.2        |
| svr     | 10.4       | 5.7        | 41.7   | 14.7   | 7.6        | 5.0        |
| lio     | 16.8       | 6.4        | 59.3   | 23.0   | 5.9        | 3.9        |
| h-naive | 90.3       | 29.3       | 72.7   | 25.3   | 11.3       | 7.1        |
| d-naive | 13.2       | 4.8        | 64.8   | 22.3   | 17.3       | 13.3       |
| w-naive | 10.8       | 4.6        | 49.1   | 15.7   | 13.0       | 9.1        |

and the electricity price shown in Fig. 3. Finally, the performance of kio is comparable to the performance of other machine-learning techniques such as krr or svr . In the non-sync case, the aggregator behaves as a price-responsive EV fleet with a linear dependence and thus both kio and lio models achieve the least errors in the Ω test compared to the other benchmarks. Note also that, in this case, the h-naive is the one with the least error among the naive models. However, the RMSE of the kio is decreased by 51.3%, 25.7%, and 27.6% with

Figure 5: Estimated power bounds as well as forecast and observed power for case naive-ch .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00399v4-figure-87f1f0c1-7abd-4029-a320-e7c3d78d6e73.jpg)

respect to the one attained with the models h-naive , krr , and svr , in that order. Overall, the kio model is characterized for being versatile since it makes good predictions under any pattern of the EV-fleet power with the price.

Apart from the improvement in terms of RMSE and MAE of the kio against the rest of the models to learn the EV-fleet power, the proposed approach is able to provide a bid curve, as imposed by rules in electricity markets (OMIE, 2019). Figures 5-7 show the results for cases naive-ch , sync , and non-sync , respectively. In Fig. 6.(a) and 7.(a), we show the estimated marginal utilities for the six blocks for each hour of the first day of the Ω test and for the cases sync and non-sync . In those figures, we also show the decreasing bid curves at hour 5 in the inset plots, which are also presented in Tables 4-5. Correspondingly, Fig. 5, 6.(b), 7.(b) depict the estimated bounds as well as the forecast and observed EV-fleet power for such a day.

In the naive-ch case, the kio provides coincident power bounds, as illustrated in Fig. 5, which means that the optimality problem (i.e. the marginal utility estimation problem, which captures the price effect) is useless and thus the aggregate charging power can be directly explained by estimating the bounds. In Fig. 6.(a) and 7.(a), we can observe that the kio model identifies whether the EV-fleet power is price-responsive or not by assigning different values to the marginal utility for each block. On the one hand, in Fig. 6.(a), the blockwise marginal utilities are almost identical at any time period, thus suggesting an

Figure 6: Results for case sync : (a) Estimated marginal utility price per block (in grey) and electricity price (in black) and (b) estimated power bounds as well as forecast and observed power. Note that the inset plot represents the bid price function and the corresponding electricity price of hour 5.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00399v4-figure-9f3c32fe-022e-4eb5-b8bf-14c292acb97f.jpg)

almost all-or-nothing price response of the EV fleet for the sync case. In this case, the power bounds are basically shaping the EV-fleet charging forecast. On the other hand, for the non-sync case, the bounds are generally wider than those obtained for the sync case (see Fig. 7.(b)). The marginal utility is thus shaping the aggregate power forecast since the kio model gives rise to a wider range of marginal utility values at any time period, as can be observed in Fig. 7.(a). In short, unlike any other forecasting tool, we gain interpretability with the proposed IO approach kio due to two aspects: (i) the width of the bounds, which sheds light on the price-responsiveness of the EV fleet; and (ii) the derivation of a bid curve when there exists a dependence of the EV-fleet power on the price, as can be seen in the inset plots of Figs. 6.(a)-7.(a) and Tables 4-5.

Figure 7: Results for case non-sync : (a) Estimated marginal utility price per block (in grey) and electricity price (in black) and (b) estimated power bounds as well as forecast and observed power. Note that the inset plot represents the bid price function and the corresponding electricity price of hour 5.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00399v4-figure-5537c858-3c45-45f4-ac53-1ba4bdad8d28.jpg)

## 4.3. Forecast Results with V2G Services

glyph[negationslash]

We now assume that EVs may enable their V2G capabilities (i.e. B d v = 0 in the model (A.1)-(A.9)) and we compare the results for two cases: (i) a highlysynchronized power case when C S = 0; and (ii) a case in which the power synchronization is avoided when C S = 52 e /MWh . Those cases are denoted 2 as sync and non-sync . The problem setup is identical to that explained in Section 4.2. Table 6 provides the error metrics on the Ω test for all models. As can be seen, kio clearly outperforms by far the lio and naive models for both cases. Notwithstanding, the performance of lio in terms of error is closer to the proposed approach for the non-sync case because the EV-fleet power is more price-responsive. Also, the performance of kio is similar to the machine-learning techniques krr and svr in the case sync ; and the RMSE (MAE) decreases by 4.8% and 5.9% (11.4% and 6.7%) compared to krr and svr , respectively, in the

Table 4: Bid Curve at Hour 5 - Case sync

| Block                      |    1 |    2 |    3 |    4 |    5 |    6 |
|----------------------------|------|------|------|------|------|------|
| Marginal utility ( € /MWh) | 42.7 | 42.4 | 42.4 | 42.4 | 42.4 | 42.4 |
| Power block (kW)           | 38.7 | 31.1 | 31.1 | 31.1 | 31.1 | 31.1 |

Table 5: Bid Curve at Hour 5 - Case non-sync

| Block                      |    1 |    2 |    3 |    4 |    5 |    6 |
|----------------------------|------|------|------|------|------|------|
| Marginal utility ( € /MWh) | 45.5 | 45.4 | 44.7 | 43.3 | 41.6 | 40.5 |
| Power block (kW)           | 26   |  8.1 |  8.1 |  8.1 |  8.1 |  8.1 |

case non-sync .

Table 6: Error Metrics - Cases with V2G Services (kW)

| Model   | sync   | sync   | non-sync   | non-sync   |
|---------|--------|--------|------------|------------|
| Model   | RMSE   | MAE    | RMSE       | MAE        |
| kio     | 148.6  | 94.3   | 33.5       | 20.9       |
| krr     | 146.9  | 108.4  | 35.2       | 23.6       |
| svr     | 147.1  | 92.4   | 35.6       | 22.4       |
| lio     | 172.1  | 120.0  | 36.2       | 23.7       |
| h-naive | 235.4  | 142.2  | 49.5       | 30.0       |
| d-naive | 261.8  | 162.5  | 71.1       | 50.2       |
| w-naive | 199.5  | 112.3  | 60.4       | 37.7       |

## 5. Conclusions

This paper proposes a data-driven two-step estimation procedure relying on two main concepts: inverse optimization and kernel regression. This novel approach allows to capture the nonlinear relationship between an aggregate priceresponse and the associated explanatory variables, while deriving a bid/offer curve, as imposed by rules in electricity markets. We apply such a framework to learn the aggregate price-response of an EV fleet. The proposed approach

attains a better performance (around 20%-40% error reduction) than naive or linear models. Moreover, it achieves a similar or better (depending on the case) performance than state-of-the-art machine-learning techniques such as support vector regression or kernel-ridge regression. Overall, the proposed approach is versatile since its performance is good regardless of the price-power relation. Very interestingly, besides, it increases the degree of interpretability of the prediction model compared to existing approaches in the literature since a bid/offer curve can be readily derived.

## 6. Acknowledgements

This project has received funding in part by the Spanish Ministry of Economy, Industry, and Competitiveness through project ENE2017-83775-P; in part by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 755705); and in part by Fundaci´n Iberdrola Espa˜a 2018. The authors thankfully acknowledge o n the computer resources, technical expertise and assistance provided by the SCBI (Supercomputing and Bioinformatics) center of the University of Malaga.

## Appendix A. Aggregator of Electric Vehicles

To simulate the behavior of a pool of EVs, i.e., its aggregate power, we assume an aggregator of EVs in residential districts who aims to minimize their

total costs. This can be mathematically expressed as:

$$\min _ { \Xi ^ { e v } } \sum _ { t \in \mathcal { T } } \left ( \lambda _ { t } \Delta t p _ { t } + \sum _ { v \in \mathcal { V } } \left ( C _ { v, t } ^ { D } + C ^ { P } s _ { v, t } \right ) + C ^ { S } \Delta t ^ { 2 } p _ { t } ^ { 2 } \right )$$

subject to:

$$p _ { t } = \sum _ { v \in \mathcal { V } } \left ( c _ { v, t } - d _ { v, t } \right ), \quad \forall t \in \mathcal { T } & & \text{(A.2)}$$

$$s o c _ { v, t } = s o c _ { v, t - 1 } + \Delta t \left ( \eta _ { v } ^ { c } c _ { v, t } - \frac { d _ { v, t } } { \eta ^ { d } } \right ) - \chi _ { v, t } + s _ { v, t } \quad \forall v \in \mathcal { V }, t \in \mathcal { T } \quad ( \text{A.3} )$$

$$0 \leq c _ { v, t } \leq B _ { v } ^ { c } \varsigma _ { v, t }, \ \forall v \in \mathcal { V }, t \in \mathcal { T }$$

$$0 \leq d _ { v, t } \leq B _ { v } ^ { d } \varsigma _ { v, t }, \quad \forall v \in \mathcal { V }, t \in \mathcal { T }$$

$$\underline { S O C } _ { v, t } \leq s o c _ { v, t } \leq \overline { S O C } _ { v, t }, \ \forall v \in \mathcal { V }, t \in \mathcal { T }$$

$$s o c _ { v, N _ { T } } & = s o c _ { v, 0 }, \quad \forall v \in \mathcal { V } & & ( A. 7 )$$

$$C _ { v, t } ^ { D } = A _ { v, t } + F _ { v } d _ { v, t }, \quad \forall v \in \mathcal { V }, t \in \mathcal { T }$$

$$s _ { v, t } \geq 0, \quad \forall v \in \mathcal { V }, t \in \mathcal { T },$$

where the set of decision variables Ξ ev = { c v,t , C D t v, , d v,t , p t , s v,t , soc v,t } , V is the set of EVs in the fleet, T is the set of time periods. The variable p t represents the power the aggregator buys in the electricity market whereas the variables c v,t and d v,t represent the charging power from and discharging power to the grid of EV v in period t . The variable C D t v, represents the cost of battery degradation due to motion and charging/discharging cycle of EV v in period t . The variables s v,t act as a load shedding term when the energy balance of the EVs cannot be satisfied. Finally, soc v,t is the state of charge of the battery of EV v in period t . In addition, λ t is the electricity price in period t ; ∆ is the time step; t C P is the load shedding cost; C S is a penalty cost to avoid power synchronization; η c d ( ) v is the charging (discharging) efficiency for the EV v ; χ v,t represents the energy required for transportation of each EV throughout the time horizon; B c v and B d v are the maximum charging and discharging power of EV v , respectively; ς v,t represents the availability of the EV v in period t ; SOC v,t and SOC v,t are the minimum and maximum limits of the energy state of charge of EV v in period t ; N T is the number of time periods; F v is the degradation cost per kW due

to charging-discharging cycles and it depends on the battery cost of EV v ; and A v,t is the degradation cost due to motion of EV v in time period t .

The problem (A.1)-(A.9) aims to minimize the total costs as given in (A.1), which comprise four terms: (i) the operational costs due to charging from and discharging to the grid, (ii) the degradation costs of the vehicles' batteries, (iii) the load shedding costs when the equation associated with the energy state-ofcharge evolution is violated, and (iv) the penalty costs to avoid power synchronization that may lead to overloads in the distribution network (Sarker et al., 2016). Constraints (A.2) relate the power bought in the electricity market with the charging and discharging power. Constraints (A.3) model the energy state of charge evolution while taking into account the energy required for transportation. Expressions (A.4) and (A.5) impose the lower and upper bounds for the charging and discharging power, in that order. Constraints (A.6) set the lower and upper bounds for the energy state of charge of the EVs. Expressions (A.7) enforce boundary conditions on the energy state-of-charge of the EVs. Expressions (A.8) model the battery degradation costs based on the motion status and the discharging energy, as described in Ortega-Vazquez (2014). Finally, constraints (A.9) define the non-negativity character of the variable s v,t .

## References

Amini, M. H., Kargarian, A., &amp; Karabasoglu, O. (2016). ARIMA-based decoupled time series forecasting of electric vehicle charging demand for stochastic power system operation. Electric Power Systems Research , 140 , 378-390.

Arias, M. B., &amp; Bae, S. (2016). Electric vehicle charging demand forecasting model based on big data technologies. Applied energy , 183 , 327-339.

Aswani, A., Shen, Z.-J., &amp; Siddiq, A. (2018). Inverse optimization with noisy data. Operations Research , 66 , 870-892.

Bandpey, M. F., &amp; Firouzjah, K. G. (2018). Two-stage charging strategy of plug-in electric vehicles based on fuzzy control. Computers &amp; Operations Research , 96 , 236-243.

Ben´ ıtez-Pe˜a, n S., Blanquero, R., Carrizosa, E., &amp; Ram´ ırez-Cobo, P. (2019). Cost-sensitive feature selection for support vector machines. Computers &amp; Operations Research , 106 , 169-178.

Bertsimas, D., Gupta, V., &amp; Paschalidis, I. C. (2015). Data-driven estimation in equilibrium using inverse optimization. Math. Progr. , 153 , 595-633.

Bunsen, T., Cazzola, P., Gorner, M., Paoli, L., Scheffer, S., Schuitmaker, R., Tattini, J., &amp; Teter, J. (2018). Global EV Outlook 2018: Towards cross-modal electrification . Technical Report - International Energy Agency.

Chan, T. C., &amp; Kaw, N. (2019). Inverse optimization for the recovery of constraint parameters. European Journal of Operational Research , .

CPLEX (2019). IBM ILOG CPLEX Optimisation Studio. [Online]. Available: https://www.ibm.com/analytics/cplex-optimizer.

Dudek, G. (2018). Probabilistic forecasting of electricity prices using kernel regression. In 2018 15th International Conference on the European Energy Market (EEM) (pp. 1-5).

ENTSO-e (2019). European Network of Transmission System Operators (ENTSO-e). [Online]. Available: https://www.entsoe.eu/.

Esfahani, P. M., Shafieezadeh-Abadeh, S., Hanasusanto, G. A., &amp; Kuhn, D. (2018). Data-driven inverse optimization with imperfect information. Math. Progr. , 167 , 191-234.

European Commission (2011). White Paper: Roadmap to a Single European Transport Area -Towards a competitive and resource efficient transport system COM(2011) 144 final. [Online]. Available: https://eur-lex.europa.eu/legalcontent/EN/TXT/PDF/?uri=CELEX:52011DC0144&amp;from=EN.

Hart, W. E., Watson, J. P., &amp; Woodruff, D. L. (2011). Pyomo: Modeling and solving mathematical programs in Python. Mathematical Programming Computation , 3 , 219-260. doi: 10.1007/s12532-011-0026-8 .

Hofmann, T., Sch¨lkopf, B., &amp; Smola, A. J. (2008). Kernel methods in machine o learning. The annals of statistics , (pp. 1171-1220).

Kekatos, V., Veeramachaneni, S., Light, M., &amp; Giannakis, G. B. (2013). Dayahead electricity market forecasting using kernels. In 2013 IEEE PES Innovative Smart Grid Technologies Conf. (ISGT) (pp. 1-5).

Kekatos, V., Zhang, Y., &amp; Giannakis, G. B. (2014). Electricity market forecasting via low-rank multi-kernel learning. IEEE Journal of Selected Topics in Signal Processing , 8 , 1182-1193.

Kempton, W., &amp; Tomi´ c, J. (2005). Vehicle-to-grid power fundamentals: Calculating capacity and net revenue. Journal of power sources , 144 , 268-279.

Korolko, N., Sahinoglu, Z., &amp; Nikovski, D. (2015). Modeling and forecasting self-similar power load due to EV fast chargers. IEEE Transactions on Smart Grid , 7 , 1620-1629.

Lu, T., Wang, Z., Wang, J., Ai, Q., &amp; Wang, C. (2018). A data-driven stackelberg market strategy for demand response-enabled distribution systems. IEEE Transactions on Smart Grid , 10 , 2345-2357.

Majidpour, M., Qiu, C., Chu, P., Pota, H. R., &amp; Gadh, R. (2016). Forecasting the EV charging load based on customer profile or station measurement? Applied energy , 163 , 134-141.

NHTS (2019). National Household Travel Survey (NHTS) data. [Online]. Available: https://nhts.ornl.gov/.

OMIE (2019). Spanish Market Operator. [Online]. Available: http://www.omie.es/.

Ortega-Vazquez, M. A. (2014). 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 , 8 , 1007-1016.

Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., &amp; Duchesnay, E. (2011). Scikitlearn: Machine learning in Python. Journal of Machine Learning Research , 12 , 2825-2830.

Renault (2019). Data Sheet Renault ZE ZOE. [Online]. Available: http://electromovilidad.net/renault/ficha-tecnica-renault-zoe/.

Roland, J., Figueira, J. R., &amp; De Smet, Y. (2016). Finding compromise solutions in project portfolio selection with multiple experts by inverse optimization. Computers &amp; Operations Research , 66 , 12-19.

Ruiz, C., Conejo, A. J., &amp; Bertsimas, D. J. (2013). Revealing rival marginal offer prices via inverse optimization. IEEE Transactions on Power Systems , 28 , 3056-3064.

Saez-Gallego, J., &amp; Morales, J. M. (2017). Short-term forecasting of priceresponsive loads using inverse optimization. IEEE Transactions on Smart Grid , 9 , 4805-4814.

Saez-Gallego, J., Morales, J. M., Zugno, M., &amp; Madsen, H. (2016). A data-driven bidding model for a cluster of price-responsive consumers of electricity. IEEE Transactions on Power Systems , 31 , 5001-5011.

Sarker, M. R., Ortega-Vazquez, M., Kirschen, D., &amp; Zhang, B. (2016). Electric Vehicles as Grid Resources . Ph.D. thesis University of Washington Libraries.

Shahidehpour, M., Yamin, H., &amp; Li, Z. (2003). Market operations in electric power systems: forecasting, scheduling, and risk management . John Wiley &amp; Sons.

Smola, A. J., &amp; Sch¨lkopf, B. (2004). o A tutorial on support vector regression. Statistics and computing , 14 , 199-222.

Sun, B., Sun, X., Tsang, D. H., &amp; Whitt, W. (2019). Optimal battery purchasing and charging strategy at electric vehicle battery swap stations. European Journal of Operational Research , 279 , 524-539.

Sun, Q., Liu, J., Rong, X., Zhang, M., Song, X., Bie, Z., &amp; Ni, Z. (2016). Charging load forecasting of electric vehicle charging station based on support vector regression. In 2016 IEEE PES Asia-Pacific Power and Energy Eng. Conf. (APPEEC) (pp. 1777-1781).

Trevor, H., Robert, T., &amp; JH, F. (2009). The elements of statistical learning: data mining, inference, and prediction.

Xydas, E., Marmaras, C., Cipcigan, L. M., Hassan, A., &amp; Jenkins, N. (2013). Forecasting electric vehicle charging demand using support vector machines. In 2013 48th International Universities' Power Engineering Conference (UPEC) (pp. 1-6).

Yang, J., &amp; Sun, H. (2015). Battery swap station location-routing problem with capacitated electric vehicles. Computers &amp; Operations Research , 55 , 217-232.

Zhang, J., &amp; Xu, C. (2010). Inverse optimization for linearly constrained convex separable programming problems. European Journal of Operational Research , 200 , 671-679.

Zhou, Y., Wang, L., &amp; McCalley, J. D. (2010). Effective incentives design for renewable energy generation expansion planning: An inverse optimization approach. In IEEE PES General Meeting (pp. 1-7).
Ricardo Fernández-Blanco Juan Miguel Morales Salvador Pineda Álvaro Porras
Inverse Optimization with Kernel Regression: Application to the Power Forecasting and Bidding of a Fleet of Electric Vehicles
This paper considers an aggregator of Electric Vehicles (EVs) who aims to
learn the aggregate power of his/her fleet while also participating in the
electricity market. The proposed approach is based on a data-driven inverse
optimization (IO) method, which is highly nonlinear. To overcome such a caveat,
we use a two-step estimation procedure which requires solving two convex
programs. Both programs depend on penalty parameters that can be adjusted by
using grid search. In addition, we propose the use of kernel regression to
account for the nonlinear relationship between the behaviour of the pool of EVs
and the explanatory variables, i.e., the past electricity prices and EV fleet's
driving patterns. Unlike any other forecasting method, the proposed IO
framework also allows the aggregator to derive a bid/offer curve, i.e. the
tuple of price-quantity to be submitted to the electricity market, according to
the market rules. We show the benefits of the proposed method against the
machine-learning techniques that are reported to exhibit the best forecasting
performance for this application in the technical literature.