Carlos Adrian Correa-Florez
email: carlos-adrian.correaflorez@mines-paristech.fr

Andrea Michiorri

Georges Kariniotakis

Robust Optimization for Day-ahead Market Participation of Smart-Home Aggregators

Keywords: Energy storage, robust optimization, uncertainty, battery cycling, residential aggregator. Nomenclature Abbreviations ARO Adjustable Robust Optimization BESS Battery Energy Storage System

This paper proposes an optimization model to participate in day-ahead energy markets when PV generation, thermal and electro-chemical storage devices are aggregated at the residential level. The model includes uncertainty in energy prices, PV and load; and adjustable robust optimization is used to determine a tractable counterpart of the problem. By means of robust control parameters, solutions with different levels of conservatism can be found and analyzed. In addition, the presented model includes explicit representation of battery degradation by means of special ordered sets. This equivalent cycling aging calculation takes into account the non-linear relation between depth of discharge and total life cycles of the battery by piecewise linearization. Performance analysis shows the advantage of the proposed approach when compared to the deterministic solution in terms of average cost and risk. For the analyzed real-life test system, the robust formulation achieves cost reduction of up to 5.7% and standard deviation decreases as much as 36.4%.

Introduction

Increasing penetration of decentralized renewable generation into medium-and low-voltage grids is motivating the development of new tools to overcome the challenges imposed by this new paradigm. These trends have extended even further to reach the building and home level, leading to the development of concepts such as Home Energy Management Systems (HEMS) [START_REF] Zhou | Smart home energy management systems: Concept, configurations, and scheduling strategies[END_REF]. In the smartgrid context, the flexibility features of renewables, storage technologies, demand response (DR) and interaction with the grid [START_REF] Aduda | Demand side flexibility : Potentials and building performance implications[END_REF] can be exploited by different market agents to minimize operation costs. In the concrete case of the present work, the aim is to analyze the interaction between thermal and electric storage for an aggregation of smart homes including uncertainties in energy prices, load and PV production, and also considering battery aging.

Current research

Regarding management models for joint thermal and electric storage technologies at the residential level, approaches include that presented in [START_REF] Comodi | Multi-apartment residential microgrid with electrical and thermal storage devices : Experimental analysis and simulation of energy management strategies[END_REF], which proposes a residential microgrid in which thermal and electric storage make it possible to shave the demand peak and enhance the system's self-sufficiency. The approach in [START_REF] Li | PV integration in Low-Voltage feeders with Demand Response[END_REF] presents a methodology for intraday management of PV and Electric Water Heaters (EWH) in an LV network, with the EWH acting as a flexible load rather than a storage device.

Reference [START_REF] Good | Optimization under Uncertainty of Thermal Storage-Based Flexible Demand Response with Quantification of Residential Users' Discomfort[END_REF] presents an optimization problem for the day-ahead market that minimizes retailer costs represented by imports/exports, and gas costs, along with expected balancing costs in real-time operation. The model does not include Battery Energy Storage Systems (BESS), but does include thermal load and electro-thermal storage. Sizing and operation of storage devices in smart buildings is presented in [START_REF] Xu | Performance Analysis and Comparison on Energy Storage Devices for Smart[END_REF], including electrical and thermal storage, but disregarding cycling effects.

Reference [START_REF] Ouammi | Optimal Power Scheduling for a Cooperative Network of Smart Residential Buildings[END_REF] presents a cooperative scheme of Smart Residential Buildings (SRB), considering batteries, thermal storage and electric vehicles. Although cycling is not taken into account, this study constitutes an interesting benchmark given that different network configurations and interactions are presented.

A multi-energy microgrid was recently proposed in [START_REF] Holjevac | Corrective receding horizon scheduling of flexible distributed multi-energy microgrids[END_REF], in which thermal and electrical storage and heat sources are used to reduce operation costs and alleviate network capacity issues at the Point of Common Coupling (PCC). Although this paper does not account for either battery cycling or uncertainties, it does present a thorough modeling of different energy sources and their interactions, and is tested on a system comprising 300 households.

In this smart grid context, uncertainty plays an important role in the decision-making process. One common practice to facilitate these optimization processes is Stochastic Optimization (SO), which typically aims to determine the optimal solution among a number of expected predefined scenarios [START_REF] Birge | Introduction to Stochastic Programming[END_REF]. However, drawbacks of SO include factors such as the requirements for probabilistic information of uncertain variables, the implementation of specialized scenario generation/reduction techniques, and the computational burden related to large number of scenarios.

An alternative approach which has gained substantial attention in recent years is Robust Optimization (RO) [START_REF] Bertsimas | Robust discrete optimization and network flows[END_REF], which is an interval-based optimization method. RO does not require knowledge of the Probability Density Function (PDF) of uncertain variables, but rather requires moderate information, i.e. an uncertainty set for each uncertain variable. RO provides a robust optimal solution that is feasible (immunized) within the confidence interval.

RO has been successfully used to tackle uncertainty, mainly in large-scale power systems problems with a variety of objectives. For instance, it has been used to capture load and wind uncertainty in Unit Commitment (UC) [START_REF] Bertsimas | Adaptive Robust Optimization for the Security Constrained Unit Commitment Problem[END_REF]. In the case of large-size battery participation in energy and ancillary markets, RO was used in [START_REF] Kazemi | Operation Scheduling of Battery Storage Systems in Joint Energy and Ancillary Services Markets[END_REF] to capture uncertainty in prices. In transmission expansion planning, this methodology has been used to cope with demand and renewable generation uncertainty [START_REF] Jabr | Robust Transmission Network Expansion Planning With Uncertain Renewable Generation and Loads[END_REF]. In [START_REF] Thatte | Risk Measure Based Robust Bidding Strategy for Arbitrage Using a Wind Farm and Energy Storage[END_REF], strategic bidding for a wind farm and battery was achieved by including price and wind power uncertainty. Stand-alone wind systems for market participation have also included RO analysis [START_REF] Hu | Bi-level robust dynamic economic emission dispatch considering wind power uncertainty[END_REF].

Although most RO applications are related to large power systems, the growing interest in decentralized and distributed energy has pushed the research community to explore this approach.

Although little research has been done on exploiting RO capabilities in residential storage-based energy systems, some work has been published in recent years, specifically related to medium-size DG/microgrid management. For instance, [START_REF] Wang | Optimal bidding strategy for microgrids in joint energy and ancillary service markets considering flexible ramping products[END_REF] presents a model for strategic bidding in energy and ancillary markets for a microgrid consisting of RES, a microturbine (MT) and a battery, in which RO is used to include RES uncertainty and SO is used to tackle price uncertainty. For bidding purposes in day-ahead and real-time markets, reference [START_REF] Liu | Bidding Strategy for Microgrid in Day-Ahead Market Based on Hybrid Stochastic/Robust Optimization[END_REF] proposes a hybrid stochastic/robust approach, in which RO captures uncertainty in real-time prices, while stochastic optimization is used to include wind and PV scenarios. Both approaches ( [START_REF] Wang | Optimal bidding strategy for microgrids in joint energy and ancillary service markets considering flexible ramping products[END_REF][START_REF] Liu | Bidding Strategy for Microgrid in Day-Ahead Market Based on Hybrid Stochastic/Robust Optimization[END_REF]) assume deterministic demand.

Robust resource scheduling of MG components is analyzed in [START_REF] Hu | Toward Optimal Energy Management of Microgrids via Robust Two-Stage Optimization[END_REF][START_REF] Liu | Robust optimisation-based microgrid scheduling with islanding constraints[END_REF] with uncertainty in load and RES, but neglecting uncertainty in price information. In [START_REF] Wang | Integrated scheduling of energy supply and demand in microgrids under uncertainty: A robust multi-objective optimization approach[END_REF], a DR program for industrial customers is presented. Uncertainties with RO are considered in load and PV, and a multi-objective algorithm is used for minimizing cost and emissions. In [START_REF] Kuznetsova | An integrated framework of agent-based modelling and robust optimization for microgrid energy management[END_REF], robust energy management is achieved for an MG consisting of a train station and a district. The model features complete uncertainty inclusion for PV, wind, load and energy prices. Although this model includes robustness for all uncertain parameters, it does not analyze different levels of uncertainty budget, immunizing the solution against any realization of uncertainty, but leading to over-conservative solutions.

Reference [START_REF] Gupta | A robust optimization based approach for microgrid operation in deregulated environment[END_REF] proposes energy and reserve market participation using RO for wind uncertainty and considers PV and dispatchable units at MV level. Also at MV level, reference [START_REF] Wang | A robust optimization approach for energy generation scheduling in microgrids[END_REF] includes uncertainty in net and heat demand with chance-constrained optimization and price uncertainty with RO.

Some of the previous references, consider a full uncertainty budget, i.e. worst case realization of uncertain variables to protect against realization [START_REF] Kuznetsova | An integrated framework of agent-based modelling and robust optimization for microgrid energy management[END_REF][START_REF] Hu | Toward Optimal Energy Management of Microgrids via Robust Two-Stage Optimization[END_REF]. However, these solutions could lead to over-conservatism, given that it is very unlikely that all uncertainties take extreme values at the same time. Hence, robust parameter analysis can be introduced to achieve less expensive solutions [START_REF] Bertsimas | The Price of Robustness[END_REF]. This logic is used in [START_REF] Gupta | A robust optimization based approach for microgrid operation in deregulated environment[END_REF][START_REF] Wang | A robust optimization approach for energy generation scheduling in microgrids[END_REF][START_REF] Liu | Bidding Strategy for Microgrid in Day-Ahead Market Based on Hybrid Stochastic/Robust Optimization[END_REF] to present a sensitivity analysis for different values of the robust parameters, given that only one source of uncertainty is considered. When several sources of uncertainty are included, the interaction of different robust parameters can lead to promising solutions. However, in [START_REF] Wang | Optimal bidding strategy for microgrids in joint energy and ancillary service markets considering flexible ramping products[END_REF][START_REF] Wang | Integrated scheduling of energy supply and demand in microgrids under uncertainty: A robust multi-objective optimization approach[END_REF][START_REF] Liu | Robust optimisation-based microgrid scheduling with islanding constraints[END_REF], these interactions are neglected, and instead, all robust parameters are forced to assume the same arbitrary values.

Although some previous MG management referenced papers include at least one battery in their respective test systems, none includes the non-linear relation of DoD to account for impacts on degradation and cycling aging. At most, a simplified linear cost (function of power charge and discharge) is included in [START_REF] Liu | Bidding Strategy for Microgrid in Day-Ahead Market Based on Hybrid Stochastic/Robust Optimization[END_REF][START_REF] Hu | Toward Optimal Energy Management of Microgrids via Robust Two-Stage Optimization[END_REF][START_REF] Liu | Robust optimisation-based microgrid scheduling with islanding constraints[END_REF], and references [START_REF] Kuznetsova | An integrated framework of agent-based modelling and robust optimization for microgrid energy management[END_REF][START_REF] Wang | Optimal bidding strategy for microgrids in joint energy and ancillary service markets considering flexible ramping products[END_REF][START_REF] Wang | Integrated scheduling of energy supply and demand in microgrids under uncertainty: A robust multi-objective optimization approach[END_REF] neglect cycling aging impacts. In this work, we propose an explicit modeling of degradation by means of piecewise linearization of the curve that describes the non-linear relation of DoD and equivalent life cycles.

Very little work has been published regarding home-level storage management using RO. Robust management for home appliances is presented in [START_REF] Wang | Robust Optimization for Household Load Scheduling with Uncertain Parameters[END_REF] to minimize electricity bills in a single house and including uncertainty in comfort variables. Regarding robust aggregation of storage at the residential level, reference [START_REF] Ampatzis | Robust optimisation for deciding on real-time flexibility of storage-integrated photovoltaic units controlled by intelligent software agents[END_REF] proposes a scheme for real-time decision-making considering batteries and price uncertainty. Reference [START_REF] Diekerhof | Hierarchical Distributed Robust Optimization for Demand Response Services[END_REF] does not include battery aggregation, but instead considers exploitation of thermal storage at the residential level in a 20-household testbed, using RO to account only for thermal demand uncertainty. Although this research does not include price or electrical load uncertainty, or RES integration, is does give an interesting insight into the scalability of the proposed model. Reference [START_REF] Zhou | A robust optimization approach for integrated community energy system in energy and ancillary service markets[END_REF] presents a community energy management system disregarding batteries, but including PV and wind power. RO is used to include uncertainty in RES and prices.

Despite its valuable contribution on participation in multiple markets, reference [START_REF] Zhou | A robust optimization approach for integrated community energy system in energy and ancillary service markets[END_REF] lacks analysis of the budget of uncertainty regarding the three considered robust parameters. The authors predefine parameters for both outdoor temperature and hot water use, and then presents a sensibility analysis by adjusting a parameter that influences both uncertain variables to the same degree. In addition, reference [START_REF] Diekerhof | Hierarchical Distributed Robust Optimization for Demand Response Services[END_REF] does not consider adjustable parameters and only presents the worst case solution.

A more common practice in the specialized literature to account for aggregation of residential/building storage under uncertainty, is by using Model Predictive Control (MPC) and/or stochastic optimization. For instance, stochstic optimization and chance constrained methods are used in [START_REF] Kilkki | Optimizing residential heating and energy storage flexibility for frequency reserves[END_REF] for energy and reserve market participation by aggregating residential batteries and heating. This work presents inclusion of uncertainty in prices, weather and realized frequency. MPC based models for aggregation of distributed storage devices (at DSO level) is presented in [START_REF] Mgel | Scheduling distributed energy storage units to provide multiple services[END_REF] to provide local and frequency services; in [START_REF] Li | Aggregator service for PV and battery energy storage systems of residential building[END_REF] an aggregation of appartment buildings with multiple EWHs, EVs and a single battery is presented; and MPC is also used in [START_REF] Rahnama | Robust Aggregator Design for Industrial Thermal Energy Storages in Smart Grid[END_REF] for the case of regulation services by aggregation of industrial thermal loads.

In contrast with [START_REF] Kilkki | Optimizing residential heating and energy storage flexibility for frequency reserves[END_REF][START_REF] Ampatzis | Robust optimisation for deciding on real-time flexibility of storage-integrated photovoltaic units controlled by intelligent software agents[END_REF][START_REF] Diekerhof | Hierarchical Distributed Robust Optimization for Demand Response Services[END_REF], the models presented in [START_REF] Mgel | Scheduling distributed energy storage units to provide multiple services[END_REF][START_REF] Li | Aggregator service for PV and battery energy storage systems of residential building[END_REF][START_REF] Rahnama | Robust Aggregator Design for Industrial Thermal Energy Storages in Smart Grid[END_REF][START_REF] Wang | Robust Optimization for Household Load Scheduling with Uncertain Parameters[END_REF] do not consider aggregation/coordination of storage at the residential level. In addition, when storage is considered in [START_REF] Wang | Robust Optimization for Household Load Scheduling with Uncertain Parameters[END_REF][START_REF] Ampatzis | Robust optimisation for deciding on real-time flexibility of storage-integrated photovoltaic units controlled by intelligent software agents[END_REF][START_REF] Diekerhof | Hierarchical Distributed Robust Optimization for Demand Response Services[END_REF][START_REF] Kilkki | Optimizing residential heating and energy storage flexibility for frequency reserves[END_REF][START_REF] Mgel | Scheduling distributed energy storage units to provide multiple services[END_REF][START_REF] Li | Aggregator service for PV and battery energy storage systems of residential building[END_REF], a single technology is explored (either thermal or electrochemical storage). Compared to the above references, we present a model that treats exploitation of both electrochemical and thermal storage installed at the residential level. It is important to point out that over the past few years, some research has addressed the joint operation of batteries and thermal storage for different purposes: peak shaving and system self sufficiency [START_REF] Comodi | Multi-apartment residential microgrid with electrical and thermal storage devices : Experimental analysis and simulation of energy management strategies[END_REF], device sizing and operation [START_REF] Xu | Performance Analysis and Comparison on Energy Storage Devices for Smart[END_REF], smart building cooperation [START_REF] Ouammi | Optimal Power Scheduling for a Cooperative Network of Smart Residential Buildings[END_REF], deterministic cost minimization and network capacity alleviation [START_REF] Holjevac | Corrective receding horizon scheduling of flexible distributed multi-energy microgrids[END_REF]. Note, however, that the model we propose pursues a different objective, which is to provide a robust framework to participate in day-ahead energy markets while capturing battery-cycling costs from a residential aggregator standpoint.

Advanced models for dealing with battery aging have been presented mainly for EVs [START_REF] Ortega-Vazquez | Optimal scheduling of electric vehicle charging and vehicle-to-grid services at household level including battery degradation and price uncertainty[END_REF] and ancillary service provision with BESS [START_REF] He | Optimal Bidding Strategy of Battery Storage in Power Markets Considering Performance-Based Regulation and Battery Cycle Life[END_REF][START_REF] Xu | Factoring the Cycle Aging Cost of Batteries Participating in Electricity Markets[END_REF]. These approaches capture state transitions and the respective DoD at which each cycle occurs, in order to obtain an equivalent aging cost. Given that this relation is highly non-linear, it is explicitly modeled within the optimization process by piecewise linearization of the Cycles vs. DoD curve. Explicit modeling of this process can be challenging and computationally intensive, therefore most published research disregards the degradation effects ( [START_REF] Wang | Robust Optimization for Household Load Scheduling with Uncertain Parameters[END_REF][START_REF] Li | Aggregator service for PV and battery energy storage systems of residential building[END_REF][START_REF] Ampatzis | Robust optimisation for deciding on real-time flexibility of storage-integrated photovoltaic units controlled by intelligent software agents[END_REF]) or assumes simplified behavior as a function of charge/discharge power ( [START_REF] Kilkki | Optimizing residential heating and energy storage flexibility for frequency reserves[END_REF][START_REF] Mgel | Scheduling distributed energy storage units to provide multiple services[END_REF]).

About the present work

The objective of the present paper is to contribute with a framework to optimally manage HEMS resources by integrating several aspects, such as: Electric/Thermal load and storage, battery degradation costs, home/building level management, uncertainties regarding PV generation, load and energy prices. These aspects are analyzed from the standpoint of an aggregator participating in the energy market and the mathematical model corresponds to the robust tractable counterpart obtained by strong duality theory. In particular, to the best of the authors' knowledge, robust home level management of combined electric and thermal storage with load, RES and prices uncertainties has not been addressed by means of Adjustable Robust Optimization (ARO). In addition to this existing gap in the literature, we reinforce our contribution by taking into account explicit modeling of battery degradation.

For the present work, an HEMS comprising data for the 25 real-life households in the city of Evora, Portugal, is employed as a testbed. PVs, batteries and EWHs are considered to minimize an aggregator's operational costs and determine the set points for the appropriate devices. The work presented here was performed within the EU Horizon 2020 project SENSIBLE (Storage Enabled Sustainable Energy for Buildings and Communities), as part of the use case "Flexibility and demand side management in market participation". This case assumes that a retailer, or other energy service company, aggregates a number of customers, and participates in a market in order to optimize its electricity costs and add value to the flexibilities that customers can offer.

The key contributions of this paper are the following:

1) The following sources of uncertainty are considered: electrical and thermal demand, PV production and energy prices. These uncertainties are included in the model by means of Adjustable Robust Optimization theory.

2) A new methodology is proposed to detect the best robust solutions, based on Paretooptimality theory. This approach allows us to analyze the performance of multiple robust day-ahead decisions and select the Pareto front offering solutions with the best trade-off between cost and risk.

3) Explicit modeling of battery cycling cost is presented by means of special ordered sets. To our knowledge, no previous work has proposed a similar model for residential storage aggregation applications. This degradation modeling allows us to capture the non-linear relation between DoD and total life cycling to bid adequate quantities in the day-ahead markets.

The present work is organized as follows: section 2 presents the mathematical formulation of the optimization model. Next, section 3 sets out the robust counterpart and the evaluation methodology for the robust solutions. The results obtained are given in section 4 and finally, conclusions are drawn in section 5.

Framework and mathematical formulation

The proposed physical system is n fact a residential microgrid composed households. Some of which have solar panels, li-ion batteries, heat storage devices. The considered system is connected to the main grid. Each household comprises a total electrical base load to be supplied and a thermal load that has to be met by an EWH, which also stores energy in the form of heat.

The idea is to aggregate the group of houses to participate in day-ahead energy and miminize operational costs by adjusting the setting of the devices in order to optimally manage resources.

In addition, the following assumptions were made to carry out the optimization process: 1) The power exchange at the points of common coupling does not jeopardize the distribution network; 2) We presume no market barriers regarding minimum bid volumes; 3) The aggregator has communication and control capabilities with devices at the home level. An appropriate IT and communication platform exists, so that the aggregator controls devices at the home level and determines their set-points, 4) The proposed residential aggregator does not have market power, hence it acts as a price taker; and 5) Distribution network operation is part of the DSO's responsibilities and does not enter into the aggregators capabilities.

The diagram of the proposed aggregation of resources is shown in figure 1. One feature of the proposed HEMS is the possibility to independently control the BESS and TES. This means that the LV grid does not directly feed the thermostatic load. In other words, this load is fed by the available stored energy in the TES, and the input for the EWH is seen as a load from the secondary network.

Energy Market

Electrical load and PV forecasts

The SENSIBLE project [36] explores the integration of available technologies into the local power grid through three European demonstrators. The Evora demonstrator develops energy management applications and has therefore deployed smart-meters in a localized neighborhood in the city of Evora, Portugal. This smart meter roll-out provides data for forecasting models developed to predict the electrical load demand and PV production of individual households.

To predict the electrical demand of one household for the next day, the model uses the demand during the previous week and the outside temperature forecast for the next day. By means of quantile smoothing spline fitting [START_REF] Gerossier | Probabilistic day-ahead forecasting of household electricity demand[END_REF], it is possible to predict day-ahead demand Dt at instant t, as the sum of three functions:

Dt = f 1 (D t-24 ) + f 2 ( Dt ) + f 3 ( Tt ), (1)
where D t-24 is the demand 24 hours before the instant to be predicted, Dt the median demand of the previous week and Tt the predicted temperature. After quantile regression, a set of forecast quantiles is obtained. Instead of a single-point value, 10%, 20%, . . . , 90% values can be obtained and respectively associated with a 10%, 20%, . . . , 90% chance of measuring a lower actual demand at the instant predicted. This probabilistic forecasting of electrical demand is a point of interest in the literature [START_REF] Hong | Probabilistic electric load forecasting: A tutorial review[END_REF]. For further details readers are advised to review reference [START_REF] Gerossier | Probabilistic day-ahead forecasting of household electricity demand[END_REF].

The PV production forecasting model takes into account solar irradiance forecasts. Parameters such as the orientation of the PV panels, shadowing effects and other meteorological factors are estimated depending on the time frame.

Probabilistic forecasts are generated for each time of day covering the entire distribution of PV production. These quantile forecasts are given in steps of the nominal probability, hence obtaining PV forecasts associated with quantiles 10%, 20%, . . . , 90% in a similar fashion to the load forecast. For more details on the PV forecast method, reference [START_REF] Bocquet | Assessment of probabilistic PV production forecasts performance in an operational context[END_REF] is suggested as complementary reading. The forecasts are generated using the quantile spline fitting and the Kernel's density estimator, for load and PV respectively, and they are based on data collected from smart meters and Numerical Weather Predictions from ECMWF (European Centre for Medium-Range Weather Forecasts).

Quantiles 10% and 90% of the PV and load probabilistic forecasts are used to create the confidence interval that defines the budget of uncertainty in the robust optimization model (Eqs ( 33)- [START_REF] He | Optimal Bidding Strategy of Battery Storage in Power Markets Considering Performance-Based Regulation and Battery Cycle Life[END_REF]).

An example for a typical day of the normalized aggregated values for PV and load, containing the upper and lower values is shown in figure 2. 2) shows the objective function to be minimized. The first term is associated with the participation in the day-ahead market (P E t is the total energy purchase bid by the aggregator at time t) and the second term is related to the equivalent battery cycling cost. This term includes the corresponding non-linearities associated with the chemical reactions occurring in the batteries due to cycling patterns. Details and interpretation of variables X Ds t,h,s and l t,h,s are provided in subsection 2.2.4.

min T t=1 ∧ π t P E t + N h=1 S s=1 (a h,s X Ds t,h,s + b h,s l t,h,s ) (2)
Quantities a h,s and b h,s in (2), correspond to the piece-wise approximation parameters of the cycling aging cost. The methodology used to obtain these values is detailed in section 2.3.

Load balance constraints

Constraint (3) represents the power balance, in which the aggregation exchange with the wholesale market should meet the net required power by each customer of the portfolio, as also shown in figure 1.

P E t + ∆t N h=1 ( ∧ P pv t,h -P c t,h + P d t,h - ∧ D t,h -H t,h ) = 0, ∀t (3)
|P E t | ≤ ∆tP pcc , ∀t (4)
It is important to note that the net power in each house h considers battery flows, the PV injection, electrical load and the power required by the EWH. As indicated in equation ( 4), this net power constraint is limited by the maximum power allowed at the PCC, which can be given by the DSO or simply by the capacity of the transformer (4). In addition, if a household does not have an EWH that allows heat storage, the variable H equals thermal load.

Battery constraints

Constraints ( 5) -( 9) describe the energy state for the BESS. Binary variable u t,h,s is introduced to avoid charging and discharging batteries at the same time. Hence, constraints ( 7)-( 8) introduce a mixed integer characteristic into the model. Constraint [START_REF] Xu | Performance Analysis and Comparison on Energy Storage Devices for Smart[END_REF] ensures continuity of the storage devices from one day to another.

X t,h = X t-1,h + η c ∆tP c t-1,h -∆tP d t-1,h /η d , ∀t, t = 1, ∀h (5)
X 1,h = X T,h , ∀h (6)
0 ≤ P c t,h ≤ P c h • u t,h , u t,h ∈ {0, 1} , ∀t, t = 1, ∀h (7)
0 ≤ P d t,h ≤ P d h • (1 -u t,h ) (8) X h ≤ X t,h ≤ Xh (9)

Cycling constraints

To capture the non-linear relation of life cycles vs. DoD, a piecewise linearization is proposed. The proposed approach consists in identifying the beginning of each charging cycle by means of constraint [START_REF] Bertsimas | Robust discrete optimization and network flows[END_REF]. This constraint allows us to detect the transitions between an idle or charging state in t -1 to a charging state in t. In this situation, variable x t-1,h takes the value of 1, capturing the immediate time step before charging occurs. y t-1,h is a binary auxiliary variable that takes the value of 0 when no change in state occurs, or -1 when the battery stops charging. Constraint [START_REF] Kazemi | Operation Scheduling of Battery Storage Systems in Joint Energy and Ancillary Services Markets[END_REF] ensures mutually exclusive unitary value of the special ordered sets.

x t-1,h -y t-1,h = u t,h -u t-1,h , ∀t, t = 1, ∀h (10)
x

T,h -y T,h = u 0,h -u T,h , ∀h (11)
x t,h + y t,h ≤ 1, x t,h , y t,h ∈ {0, 1} , ∀t, ∀h (12)
With the identification of the beginning of a charging cycle, the DoD at which this cycle occurs can be extracted. Constraint ( 14) allows a value different from zero to be stored in X D , right before the beginning of a charging cycle. To assign the correct value of depth of discharge, constraint ( 13) is used. The right-hand side of this constraint calculates the DoD in per unit (p.u.) of rated battery energy. X Df is a dummy variable to balance the equation each time a charging cycle is not identified (x t,h = 0) and activated through constraint [START_REF] Hu | Bi-level robust dynamic economic emission dispatch considering wind power uncertainty[END_REF].

X D t,h + X Df t,h = 1 -X t,h /E rated h , ∀t, ∀h (13)
X D t,h ≤ x t,h , ∀t, ∀h (14)
X Df t,h ≤ 1 -x t,h , ∀t, ∀h (15)
Note that the definition of DoD is not subject to consensus in the literature. In this paper, we base our aging model on the assumption that each full charging cycle is accompanied by another full discharging event. Even when the discharging event occurs separately or partially, this assumption holds true with equation ( 6). This constraint ensures that all charging events must equal the discharging level in order to reach the same initial and final SOCs each day.

With constraints ( 10)- [START_REF] Hu | Bi-level robust dynamic economic emission dispatch considering wind power uncertainty[END_REF], the DoD at which each charging cycle occurs is identified as (X D ). To extract the appropriate piecewise cost function, the corresponding segment of the cost curve must be active. This is achieved by means of constraint ( 16)- [START_REF] Hu | Toward Optimal Energy Management of Microgrids via Robust Two-Stage Optimization[END_REF]. Constraint [START_REF] Liu | Bidding Strategy for Microgrid in Day-Ahead Market Based on Hybrid Stochastic/Robust Optimization[END_REF] forces the identified DoD to fall within the corresponding linearization segment and also triggers the activation of a binary variable l t,h,s for the related active segment s.

S s=1 X Ds t,h,s = X D t,h , ∀t, ∀h (16)
l min s l t,h,s ≤ X Ds t,h,s ≤ l max s l t,h,s , l t,h,s ∈ {0, 1} ∀t, ∀h (17)
S s=1 l t,h,s = x t,h , ∀t, ∀h (18)
Note that l t,h,s is efficiently used in the objective function [START_REF] Aduda | Demand side flexibility : Potentials and building performance implications[END_REF] to include the parameter b h,s when needed. In a similar manner, X Ds t,h,s is used as the independent variable of the aging cost function. The piece-wise model is also presented in subsection 2.3.

EWH constraints

The energy state of the TES is given by ( 19) and [START_REF] Wang | Integrated scheduling of energy supply and demand in microgrids under uncertainty: A robust multi-objective optimization approach[END_REF]. It is important to mention that the term including R and C (thermal resistance and capacitance, respectively) in equation [START_REF] Liu | Robust optimisation-based microgrid scheduling with islanding constraints[END_REF], represents the energy dissipation in the EWH as a measure of energy loss, as proposed in [START_REF] Good | Optimization under Uncertainty of Thermal Storage-Based Flexible Demand Response with Quantification of Residential Users' Discomfort[END_REF].

Y t,s,h = Y t-1,s,h + ∆tH t-1,s,h -Y t-1,s,h /R h C h -∆t ∧ Q t-1,s,h , ∀t, t = 1, ∀s, ∀h (19)
Y 1,s,h = Y T,s,h , ∀s, ∀h (20)
Y t,h,s ≤ Y t,h,s ≤ Ȳt,h,s (21)
0

≤ H t,h,s ≤ Ht,h,s (22)
The previous Mixed-Integer Linear Programming (MILP) mathematical model has four sources of uncertainty; i.e. energy prices ( ∧ π t in (2)), PV production (

∧ P pv t,h in (3)), electrical demand ( ∧ D t,h in (3)
) and thermal demand (

∧ Q t,h in ( 19
)
). The alternative to find a tractable optimization problem under uncertainty is explained in section 3.

Linearization of the battery degradation costs function

To include the aging cost of the batteries in the optimization model described in the previous subsection, a piecewise linearization of the degradation behavior should be carried out.

Battery life in general can be expressed in terms of the actual lifespan of the device (calendar life) or the number of achievable charge and discharge cycles (cycle life) [START_REF] Linden | Handbook of batteries[END_REF]. As already mentioned, the aging process is complex; it depends on the cycling patterns, rates of charging/discharging, and consequent chemical reactions resulting in an accumulated history of voltages, currents and temperatures [START_REF] Dutt | Life Cycle Analysis and Recycling Techniques of Batteries used in Renewable Energy Applications[END_REF], and detailed analysis of this set of interactions is beyond the scope of this paper. This paper considers battery degradation costs as a function of the cycle life's intrinsic behavior and as a function of the DoD.

From the perspective of short term electrical markets, such as energy day-ahead, cycle life allows to calculate an equivalent operation cost of using the batteries to trade energy and then add it to the total cost of the aggregator. This alternative permits more accurate daily calculation based on the cycling pattern and embed this information into the optimization model as also proposed in previous research [START_REF] He | Optimal Bidding Strategy of Battery Storage in Power Markets Considering Performance-Based Regulation and Battery Cycle Life[END_REF][START_REF] Duggal | Short-Term Scheduling of Thermal Generators and Battery Storage With Depth of Discharge-Based Cost Model[END_REF][START_REF] Sarker | Optimal Operation and Services Scheduling for an Electric Vehicle Battery Swapping Station[END_REF].

When battery manufacturers provide the curve of cycles Vs DoD, these data indirectly contain calendar life information, given that the curve is built based on capacity loss at End Of Life (EOL). This available information allows explicitly modelling cycling into the optimization by means of special ordered sets, as described in subsection 2.2.4. Cycle life is also useful when not only participation in energy market is proposed, but also frequency support services in order to properly capture and value charging/discharging cycles.

In general, the maximum number of charge/discharge cycles for a battery at a certain value d of DoD, is given by the following expression [START_REF] He | Optimal Bidding Strategy of Battery Storage in Power Markets Considering Performance-Based Regulation and Battery Cycle Life[END_REF]:

n d = n 100 d -kp (23)
where k p is a constant that depends on the life cycle -DOD curve given by the manufacturer, and can be extracted from curve fitting. The quantity n 100 is the equivalent number of cycles before failure for d = 100%.

For a certain DoD at the beginning of a charging cycle (as explained in subsection 2.2.4), an equivalent aging cost can be calculated by means of:

C cyc (d j ) = C ini n 100 d kp j (24)
The obtained C cyc is the equivalent cost due to the battery's aging process for that specific value of DoD, and should be added to the total dispatch cost. The accumulated cycling cost of aging for the identified DoDs is shown in the second term of equation ( 2), and the parameters a h,s , b h,s are obtained by linearizing equation ( 24) in segments. An example for a 3 kW / 3.3 kWh li-ion battery is shown in figure 3. In this case, a 5-segment linearization is performed for the referred BESS with the mentioned parameters. First, six equally spaced points are evaluated with equation [START_REF] Bertsimas | The Price of Robustness[END_REF]. Next, curve fitting is performed within each segment and parameters a h,s and b h,s are obtained. Finally, these parameters are fed into the optimization model.

Robust counterpart

The robust counterpart of the deterministic optimization problem described in the previous section can be found by maximizing the deviation of the uncertain parameters within each constraint. A tractable resulting problem is obtained employing strong duality theorem. Interested readers can find detailed formulation in [START_REF] Bertsimas | The Price of Robustness[END_REF].

When applying strong duality due to price uncertainty in the objective function, the following equations are obtained:

T t=1 N h=1 S s=1 (a h,s X Ds t,h,s + b h,s l t,h,s ) + T t=1 ∧ π t P E t + T t=1 q c t + Γ c z c (25)
z c + q c t ≥ 1 2 (π t -π t )y c t , ∀t (26)
-

y c t ≤ P E t ≤ y c t , ∀t (27)
z c , q c t , y c t ≥ 0, ∀t, ∀h (28)
To deal with price uncertainty, a robust parameter allows to control conservatism of the solution (Γ c ). This value can be adjusted in the range [0,T ], given that T is the maximum number of uncertain parameters in the objective function (energy price). 3) is introduced by PV and electrical load. Hence, a single uncertain right-hand parameter can be obtained by finding the net load at each time step (load minus PV). The robust counterpart of this constraint is given by equations ( 29)- [START_REF] Li | Aggregator service for PV and battery energy storage systems of residential building[END_REF].

Uncertainty in constraint (

P E t + ∆t N h=1 P d t,h -P c t,h -H t,h = D net t + q D t + Γ D t z D t , ∀t (29)
z D t + q D t ≥ 1 2 ( Dnet t -D net h )y D t , ∀t (30)
z D t , q D t ≥ 0, y D t ≥ 1, ∀t (31)
where,

D net t = 1 2 ( Dnet t + D net h ) ( 32
) Dnet t = N h=1 (D 90% t,h -P pv10% t,h ) (33)
D net t = N h=1 (D q10% t,h -P pv90% t,h ) (34)
In this case, net load uncertainty indirectly captures electrical load and PV uncertainty. Robust parameter Γ D t controls the robustness in each constraint t. For simplicity, in the remainder of this paper we eliminate the subindex t from this parameter and instead we use Γ D to control net load robustness. Cardinality of Γ D is [0,1], provided that there is one maximum uncertain parameter in each constraint.

Constraint [START_REF] Liu | Robust optimisation-based microgrid scheduling with islanding constraints[END_REF], contains another uncertain parameter, i.e. thermal load. Application of strong duality results in the following constraints:

Y t,s,h = Y t-1,s,h + ∆tH t-1,s,h -Y t-1,s,h /R h C h -∆t( ∧ Q t-1,s,h + q th t-1,h + Γ th t-1,h z th t-1,h ), ∀t, ∀h (35)
z th t,h + q th t,h ≥ 1 2 (Q 90% t,h -Q 10% t,h )y th t,h , ∀t (36)
z th t,h , q th t,h ≥ 0, y th t,h ≥ 1, ∀t, ∀h (37)
For simplicity, we eliminate the subindex t, h from Γ t,h and instead we assume a general parameter to control robustness in thermal load: Γ th . Note that Γ th ∈ [0, 1].

The complete adjustable robust MILP counterpart (ARO) is represented by the following equations:

minimize (25) (38)
s.t.

Constraints : (4) -( 18),

(20) -( 22), ( 26) -(31), ( 35)

-(37) (40)
This is a tractable MILP problem that can be solved with off-the-shelf commercial solvers. Three robust control parameters can be tuned to obtain different robust day-ahead bids: Γ DA , Γ D and Γ th . Each one controls conservatism against uncertainty in energy prices, net load and thermal load, respectively. A methodology to determine the impacts of each parameter and their interactions is presented in subsection 3.2.

Performance assessment of robust solutions

In order to analyze different levels of conservatism when bidding in the day-ahead energy market, different combinations of Γ DA , Γ D and Γ th must be analyzed. When a solution is obtained for any combination, the energy purchase commitment (P E t ) and device settings are determined. The performance of this day-ahead plan is evaluated by calculating the imbalances I - t /I + t in each time step due to energy mismatch. Negative/positive imbalance implies energy shortage/excess that has to be purchased/sold at higher/lower prices which leads to additional costs in the operation. Each plan is subject to performance analysis for several realizations of energy price, PV, and electrical and thermal load by means of Monte Carlo (MC) simulation. After several MC scenarios are generated and analyzed, the performance of the robust solution is assigned two attributes: average cost and standard deviation (SD). A decision maker is interested in minimizing both, to achieve low expectation of cost and minimize risk at the same time. Since a single day-ahead dispatch has as many average costs and SDs as combinations of Γs, Pareto optimality theory is used to select the set of control parameters that performs better for both average cost and SD [START_REF] Deb | Multi-objective Optimization using Evolutionary Algorithms[END_REF]. A complete outline of the methodology to determine the Pareto front for a specific day-ahead optimization, is shown in figure 4.

START

In this work, to generate MC scenarios we consider uniform distribution around a central forecast for energy price [0.9 ∧ π,1.1 ∧ π], and for the case of PV and load, we consider the inverse cumulative distribution function (ICDF) resulting from the real forecast described in subsection 2.1. Since quantile information is discrete, we interpolate to obtain adequate values in between quantiles. The stop criterion for MC simulation is set at 1000 scenarios. We verified confidence levels of minimum 95% and margins of error less than 1%, even below 1000 trials.

Characterization of uncertainty for numerical simulations

An index is defined to measure the amplitude of the net load in each day. This index is obtained by calculating the mean net load interval (MI) for a 24-h period. First, the maximum/minimum net load is calculated using (33)/ [START_REF] He | Optimal Bidding Strategy of Battery Storage in Power Markets Considering Performance-Based Regulation and Battery Cycle Life[END_REF]. Then, the MI for each day m is calculated by:

MI m = 1 T T t=1 Dnet t,m -D net t,m Dnet t,m × 100% (41)
the MI index gives an idea of how wide the uncertainty set is, from the net load stand point, and it is used to classify each day according to the MI level. For example, figure 5 depicts the MI for each day in November 2015. In our case, we use this information to run simulations on three representative days of high (Nov. 5th), medium (Nov. 15th) and low (Nov. 27th) MI and avoid arbitrary selection of days for day-ahead simulations. In addition, figure 6 shows an example of net load confidence intervals for two selected days: November 5th and 27th. ) and November 27th ( )

Results

Input data

The proposed algorithm is coded in Python. The electric base load is generated using the logic explained in section 2.1 for the 25 houses in the Evora demonstrator. The location of the HEMS in the real-life distribution network and the resources present in each house are shown in figure 7. The 25 houses correspond to two different LV rural networks, A and B, comprising 16 and 9 households respectively and containing the distribution of resources shown in the figure. In total, the houses feature 25 PV panels, 16 BESS and 15 EWHs. The test case is composed of only the HEMS and the control capabilities that an aggregator has over the device settings. Operation and control of the MV and LV networks are carried out by the DSO and do not come under the aggregator's capabilities or responsibilities. A normalized thermal load pattern is taken from [START_REF] Marnay | Optimal Technology Selection and Operation of Commercial-Building Microgrids -LBNL No. 62315[END_REF]. realistic behavior. Electricity prices are taken from the EPEX-European Power Exchange database [START_REF] Epex | Day ahead price for French Market[END_REF], and a persistence model is used to forecast the day-ahead and imbalance prices, consisting in assuming the last known data for the same weekday. This is done to consider a realistic case in which an aggregator, when defining day-ahead purchases, does not have the settled prices. Hence, by taking the prices for the same day in the previous week, we obtain available input data to make decisions. In line with [START_REF] Kuznetsova | An integrated framework of agent-based modelling and robust optimization for microgrid energy management[END_REF][START_REF] Ampatzis | Robust optimisation for deciding on real-time flexibility of storage-integrated photovoltaic units controlled by intelligent software agents[END_REF], price deviation from forecasted values are assumed to be ±10%, in order to create confidence intervals.

HEMS in LV Network

The charging and discharging efficiency of the batteries is assumed to be 95% in line with [START_REF] Ranaweera | Distributed control scheme for residential battery energy storage units coupled with PV systems[END_REF][START_REF] Yu | Design and Experimental Results of Battery Charging System for Microgrid System[END_REF][START_REF] Melhem | Optimization and Energy Management in Smart Home Considering Photovoltaic, Wind, and Battery Storage System With Integration of Electric Vehicles[END_REF][START_REF] Jin | Foresee: A user-centric home energy management system for energy efficiency and demand response[END_REF]. 15 batteries are rated 3kW / 3.3 kWh, and the remaining battery is a 10kW / 20 kWh device. All PV panels are rated 1.5 kWp. The cycling behavior is based on the li-ion battery information available on the market, and the curve fitting values to obtain the relation of cycle life versus DoD are taken from its technical specifications [START_REF]SAFT batteries Lithium ion battery life[END_REF]. The coefficients obtained are n 100 = 5135.7 and k p = 1.759. The battery's initial cost is EUR 500 e /kWh, according to the latest IRENA information on residential storage systems for European countries [START_REF]Battery Storage for Renewables: Market Status and Technology Outlook[END_REF].

The rated power/energy for all EWH is 1.5 kW / 3 kWh and thermal resistance/capacitance are 568 ( • C/kW)/0.3483 (kWh/ • C) in line with [START_REF] Good | Optimization under Uncertainty of Thermal Storage-Based Flexible Demand Response with Quantification of Residential Users' Discomfort[END_REF].

Operational behavior in the deterministic and robust approaches

The following results present the behavior of the main variables after solving the deterministic model presented in equations ( 2)-( 22) and the adaptive robust optimization counterpart (ARO) in ( 38)- [START_REF] Linden | Handbook of batteries[END_REF], for a medium MI day (Nov. 15th). When solving the deterministic problem, the dayahead operational cost obtained is 14.93e . On the other hand, a cost of 29.88e is found to be the robust solution when Γ DA = 24, Γ D = 1 and Γ th = 1, which corresponds to considering the full uncertainty budget for all of the uncertain variables in the model. This robust DA operation cost establishes an upper bound for the operation cost, while the deterministic solution is a lower bound. For the deterministic solution, this means that any uncertainty realization that is different from the central forecasts for prices, load and PV would imply penalization due to imbalances. On the other hand, any realization of the uncertainty set within the budget would never yield an operation cost higher than 29.88e , which acts as a guaranteed minimum.

Figure [START_REF] Holjevac | Corrective receding horizon scheduling of flexible distributed multi-energy microgrids[END_REF] shows the accumulated (all BESSs) SOC for both deterministic ( ) and ARO ( ) cases for this complete formulation of the problem. The deterministic ( ) SOC shows, in general, a similar evolution when compared to the ARO. However, a particular difference is evident in timeframes 12h-16h, which coincides with PV production hours. In the ARO, there is a discharging pattern of the batteries during these time frames provided that a highly conservative scenario is implicitly assumed within the model when Γ D =1 and Γ th =1: minimum PV production and maximum demand. This situation leads the optimal solution to set BESS in discharging mode to compensate for low levels of available PV. The ARO indicates a more conservative discharging-charging, given that stored energy during these time frames avoids a setting point close to the boundaries, in order to cope with potential uncertainty realizations in a more cost-efficient way. For instance, if during 12h-15h time frames, uncertainty realizations were higher than forecasted PV values, the ARO solution would allow more room for power injection into the grid. On the contrary, the deterministic solution would be more limited to exploit this potential situation, given that the BESS settings would not allow additional storage. However, the next subsection analyzes in more detail the potential savings and cost advantages of the ARO approach versus the deterministic approach.

In addition, we carried out a simulation that neglects cycling aging cost (second term in equation (2) for the deterministic and ARO cases. Neglecting these term in the model, means that the batteries can cycle without any related cost provided that equivalent cycling aging will have a cost of zero in the objective function. The obtained SOC for this case is depicted in figure 9 for both deterministic (

) and ARO ( ) cases. This simulation shows deeper cycling for time frames: 3h-5h, 14h-16h and 20h-22h. This is an expected result provided that the omission of aging in the model, cause the battery to cycle without degradation constraints and presents deeper and more frequent cycling. In general, cycling patterns are deeper in figure 9 than in figure 8.

The energy purchase commitment for both cases is shown in figure 10. Except for two time frames (20h and 23h), purchase commitment for the ARO approach (

) is always higher than for the deterministic approach (

). This is an expected result given that, as explained above, the robust solution is conservative and represents an upper cost bound to withstand different uncertainty realizations while maintaining final operation below this maximum guaranteed cost.

When the aggregator bids in the day-ahead energy market and determines the operation of its devices for each hour, there may be imbalances due to deviations, as explained in subsection 3.2. In this case, energy shortage/excess should be purchased/sold at the imbalance price, increasing final operational costs.

Robustness can be visualized by analyzing imbalances due to deviations when DA purchase commitment is subject to multiple random realizations of prices, PV production and demand. To generate boxplots of imbalance for deterministic and ARO solutions, MC simulation is implemented using the procedure in subsection 3.2. Figure 11 shows a comparison of the imbalances (same selected day: 15th Nov.) for each time frame. For this specific case, it can be seen that the robust solution (Γ DA = 24, Γ D = 25 and Γ th = 1) does not incur negative imbalances (additional purchases), but only sells excess energy, (positive imbalance) during all time periods, with particular higher median values during 5h and 6h and around noon. This is explained by the fact that the ARO returns a solution that is feasible for the lowest levels of PV production and the highest load scenarios, leading to energy excess when different realization of this uncertainty are analyzed. The positive imbalance during the early hours coincides with the high energy purchase (see figure 10) for battery charging (see figure 8), which becomes an energy excess and is sold back to the market. In contrast, when the deterministic solution is subject to uncertainty realizations, negative imbalance appears at all time steps, given its limited capacity for handling uncertainty realizations that are different from the central forecasted values, especially when load is higher and PV production is lower.

Interaction of devices in the deterministic and robust approaches

Another simulation was carried out for the complete month of analysis. In this case, the interaction of different sources of flexibility was analyzed for the deterministic and RO cases. The proposed test system has two technologies for storing energy: electrochemical batteries and EWHs. The possibility of using the EWHs to store thermal energy allows this device to act as a flexible load. Average costs were obtained for different cases as shown in table 1.

The table shows that for both cases (deterministic and RO) better results are obtained when TES is allowed as a form of flexibility. In specific, when EWHs are used to store heat, deterministic and RO approaches achieve 3.2% and 2.6% cost decrease respectively, when compared to cases in which only batteries are used to store energy. This shows the importance of allowing control of the EWHs' settings in order to store hot water to be used in later hours. For instance, figure 12 shows the average stored heat of all EWHs in the test system. From the figure it can be concluded that the robust approach presents higher levels of stored heat versus the deterministic case. This is explained by the fact that the RO is protected against potential variations of thermal load including the maximum values in the confidence interval. Hence the available hot water should be enough to face this extreme scenario in each time step.

The figure also shows aggregated average SOC for the BESS. In this case, stored energy for the RO case differs from deterministic mainly around noon hours, presenting lower DoD, or similarly: available stored energy. This difference in stored energy also leads to lower average cycling aging in the robust case. For instance, 10% equivalent cycling average cost decrease is achieved in RO versus the deterministic, due to the difference in the batteries' SOC.

Impacts of battery efficiency on cost

To determine the impacts of charging and discharging efficiency on the total operation cost obtained by the deterministic and ARO approaches (with Γ DA =24, Γ D =25 and Γ th =1), a sensitivity analysis is carried out.

This analysis consists on changing the value of η = η c = η d in equation 5 and solving for each day during the analyzed month. Afterwards, the monthly average cost is calculated and associated to each η value. After this procedure is completed, the obtained values are depicted in figure 13.

This figure shows that lower values of efficiency lead to average cost decrease for both cases: deterministic and ARO. For the case of the deterministic approach, when losses are neglected (η = 1) the cost is 14.74e , and increases up to 8.6% when an efficiency of 0.8 is considered for the storage system. For the ARO case, an average cost of 27.69e is obtained when losses are neglected, and cost increase reaches 8.3%. The direct comparison of the ARO and deterministic solutions shows that when losses are neglected, ARO presents an extracost of 87.85% when compared to the deterministic. Extracost of 81.1% is obtained for the case of η = 0.8, showing that despite large variations in η, the robust solutions present a higher guaranteed costs than the deterministic. This may lead to different scheduling of devices and absolute costs for each methodology, but as expected, ARO solutions always present higher maximum guaranteed cost, as also explained in section 4.2.

Effects of the budget of uncertainty

Selecting the budget of uncertainty influences the performance of the day-ahead operation, which explains the decision to use the ARO approach in the presented paper. In order to determine better budget choices, we use the procedure explained in figure 4. This allows us to use MC for several combinations of Γ DA , Γ D and Γ th , and determine average cost and SD. In this way, a set of uncertainty budgets can be determined such that the average cost and deviation are minimized. Two days are selected to develop this analysis: November 5th (High MI) and 27th (Low MI). The steps for Γ DA , Γ D and Γ th are i = 6, j = 0.2 and k = 0.2 respectively, therefore, 180 combinations are generated to be analyzed under the performance evaluation methodology. After running MC simulation for each of the 180 budgets of uncertainty, we obtain the costs and deviations in figure 14. Note that the values are normalized using the deterministic values as a base (1 p.u.). The plot shows the performance for the deterministic solution as a blue sqare: ( ).

This simulation shows that in both cases, the deterministic solution tends to have a higher standard deviation compared to the robust approach under different uncertainty budgets. This indicates that robust solutions tend to be more reliable and steady in terms of imbalances and subsequent penalization, and hence represent a less risky operation for the aggregator. In particular, the simulation for Nov. 5th shows that only 8 AROs have higher SDs than the deterministic solution, and the fourth highest for Nov. 27th. With the robust formulation, SD can be reduced by 27.3% and 36.4% for each day, respectively.

Regarding the performance costs obtained with ARO, and compared to the deterministic approach, costs can be reduced on average by 5.7% and 2.6% respectively for each day. In addition, the min-min Pareto front can be found by determining the non-dominated set of points according to Pareto optimality criteria [START_REF] Deb | Multi-objective Optimization using Evolutionary Algorithms[END_REF]. This Pareto front (depicted by: [.85]) represents the best set of solutions from the perspective of average performance cost and SD. Table 2 shows the points that comply with the following two conditions: 1) form the Pareto front and 2) independently dominate the deterministic solution. The table also shows operation cost improvement with respect to the deterministic solution, obtaining cost improvements in the ranges of 5.2%-5.7% and 2.1%-2.6% for each day, respectively. In addition, the ranges of improvement for SD are 7.7%-27.2% and 17.5%-36.4%, respectively. It should be noted that the Pareto front for Nov. 5th contains only zero values for Γ L . This result shows that for high mean net load interval days, the methodology avoids using high values of this budget of uncertainty, given the over-conservatism that these values yield. As for November 27th (low mean net load interval), ARO formulation performs in several cases even with values different from zero for all uncertainty budget parameters. This result is explained by the fact that the uncertainty realizations for load are contained in narrower intervals and close to central forecasts.

For further analysis, the ARO solutions in the Pareto front are selected to build Cumulative Density Functions (CDF) and compare performance versus the deterministic solution from a cost standpoint. As per figure 15, the ARO solutions for both days independently, behave similarly among themselves, and all share a common feature in each day, i.e. the associated costs are lower than those obtained in the deterministic approach. This result shows again the robustness of the proposed model to handle uncertain parameters and their potential realizations. For instance, for November 5th the probability that the ARO solutions will have a lower cost than the deterministic average lies in the range of 91.2% to 96.3%. For November 27th, there is a probability between 79.0% and 85.4% that the ARO solutions will perform better that the deterministic mean. These values can be found simply by intersecting the curves where cost equals 1 p.u. and determining the respective cumulative density. In general, the leftmost solutions represent better performance.

It can also be noted that the CDFs for November 27th are less scattered when compared to those of November 5th. This is again an expected solution, provided that this day presents a lower MI; thus it is expected that robust and deterministic solutions will have closer performance and behavior.

Effects of considering cycling aging in the model

To analyze this impact, the same day used in subsection 4.2 (Nov. 15th) and the following arbitrary values for the uncertainty budget are used: Γ DA =12, Γ D =0.5 and Γ th =0.5. Three cases are run and their performances are calculated with MC:

Naive Solution (NS). consists in neglecting the cycling aging cost described by equations ( 10)-( 18) and the second term in objective function [START_REF] Aduda | Demand side flexibility : Potentials and building performance implications[END_REF]. This solution allows unconstrained cycling of the batteries.

Real Cost of Naive Solution (RCNS). This solution is obtained by adding two terms: 1. obtained cost of NS, and 2. the equivalent cycling aging cost of this solution. In this way, we can obtain the real cost the NS incurs.

Complete Robust Solution (ARO). This is the solution obtained by solving the complete formulation in equations ( 38)- [START_REF] Linden | Handbook of batteries[END_REF].

After running the simulations and determining the respective performances, the following average costs are found: 14.43e , 22.04e and 14.53e for the NS, RCNS and ARO, respectively. These results show the importance of taking into account the degradation model within the optimization, so as to obtain a lower overall operation cost. The naive approach appears to have a lower cost, but as the battery is allowed to cycle without constraints, there is a hidden cost that increases operation from 14.43e to 22.04e . When using the ARO complete formulation we propose, a 34.07% cost reduction is achieved compared to the RCNS.

When the CDFs are analyzed, as depicted in figure 16 the NS could imply the achievement of lower costs, as this curve presents more leftmost points. However, the real cost of ignoring the cycling model can be observed by the rightmost CDF.

Remarks about computational time

To give an insight on the solution time for each of the solution alternatives: ARO, NS and deterministic; additional simulations were carried out. These simulations consisted on running each alternative for each day in November and determining the computational time in seconds. Figure 17 shows the obtained simulation time in each day of the month. It can be seen that the NS presents lower computational times in each day and it is an expected outcome provided that this alternative does not include the additional constraints and binary variables used to calculate cycling aging of the battery.

In addition, the deterministic case presents higher computational times in average, compared to NS. This is explained by the fact that the deterministic alternative includes cycling aging, thus increasing the size of the problem and the number of binary variables.

Table 3 shows the behavior and comparison of the computational time for each alternative. As expected, the complete ARO formulation presents higher computational times given that it includes constraints and binary variables to calculate cycling aging, and also includes the additional constraints and dual variables that appear in the robust counterpart. This leads to higher average solution time and also higher SD. Average computational times of ARO solutions are 3.6 times higher than those of the deterministic approach.

Conclusions

A model for aggregating PV, thermal storage and batteries at the residential level is presented. The model shows that including the degradation cost in the optimization model leads to cost savings due to less battery cycling. In the concrete case of the presented results, a 34.07% reduction was achieved.

The proposed approach, which can also be applied by aggregators of medium-and large-size equipment, shows that by using ARO and analyzing interactions of robust parameters, different levels of cost reduction can be achieved than with the deterministic approach. For the simulations performed, up to 5.7% cost savings were achieved. For the case of SD, the best result reduced the performance of the deterministic solution by 36.4%. This shows that not only can expected costs be reduced, but that lower risk is associated with decision-making under this approach. The results prove that using ARO also increases the probability of acheiving lower expected costs. When compared to the deterministic scheme, probabilities of up to 96.33% are obtained by analyzing CDFs.

The value of BESS efficiency changes the levels of operational cost for the aggregator. in the case of the presented simulations, cost increases for η = 0.8 reached 8.6%/8.3% for the deterministic/ARO approach, when compared to a scenario that neglects losses (η = 1).

Although computational effort of ARO is in average 3.6 times higher than the deterministic solution, the times remain in ranges lower than 20 seconds, which allow to make decisions for day-ahead practical applications without prohibitive computational efforts.

This framework can also fit into models in which other agents with batteries at the DSO/TSO levels are willing to participate in energy and ancillary markets. The robust-based performance methodology proposed in this paper, can also be applied to analyze other storage ownership configurations, remuneration schemes and risk tolerance.

The presented results could also be further explored by including intra-day operation, distribution services to be offered to the grid, and aggregation of a large set of households to diversify the

Figure 1 :

 1 Figure 1: Schematic diagram of the proposed HEMS

Figure 2 :

 2 Figure 2: Normalized load and PV forecasts for defining confidence interval

  𝑎 ℎ,𝑠 𝑋 𝑡,ℎ,𝑠 𝐷𝑠 +𝑏 ℎ,𝑠 𝑙 𝑡,ℎ,𝑠 𝑠 = 1 𝑠 = 2 𝑠 = 3 𝑠 = 4 𝑠 = 5 𝑛 100 = 5136 𝑘 𝑝 = 1,759 𝐶 𝑖𝑛𝑖 = 1650

Figure 3 :

 3 Figure 3: Example of piecewise linearization. Cycling Cost vs. DoD for a selected li-ion battery

Figure 4 :

 4 Figure 4: Flowchart of the algorithm to find best trade-off solutions

Figure 5 :

 5 Figure 5: Average net load intervals calculated with equation (41)

Figure 6 :

 6 Figure 6: Net load interval for two days. November 5th () and November 27th ( )

Figure 7 :

 7 Figure 7: Composition and location of the proposed 25-household HEMS

Figure 8 :

 8 Figure 8: Accumulated SOC for Deterministic () and ARO ( ).

Figure 9 :Figure 10 :

 910 Figure 9: Accumulated SOC when cycling cost is neglected in the model: Deterministic ( ) and ARO ( ).

Figure 11 :

 11 Figure 11: Comparison of imbalance needs for the deterministic and the robust approach when Γ DA = 24, Γ D = 1 and Γ th = 1

Figure 12 :Figure 13 :

 1213 Figure 12: Daily average operation (state of charge) of BESS and EWHs

Figure 14 :

 14 Figure 14: Average operation cost for two days under different budgets of uncertainty and Pareto front. Performance for the deterministic solution is depicted as a blue square and Pareto front by black dots

Figure 15 :

 15 Figure 15: Performance Cumulative Density Functions for selected Robust Solutions

Figure 16 :

 16 Figure 16: Cumulative Density Functions for the robust model when cycling is considered and neglected

Figure 17 :

 17 Figure 17: Solution time for ARO ( ), NS ( ) and deterministic ( ) case

Table 1 :

 1 Average daily costs of resource management

		Used flexibility Equiv. cost Total operation
	Case	BESS	TES	cycling [e ]	cost [e ]
	Deterministic	Yes	No	0.20	15.7
	Deterministic	Yes	Yes	0.20	15.19
	RO	Yes	No	0.19	28.89
	RO	Yes	Yes	0.18	28.14

Table 2 :

 2 Details of the Pareto fronts values and improvement with respect to deterministic solution Day Γ DA / Γ L / Γ th

	Cost [p.u.] SD [p.u.]

Table 3 :

 3 Daily average computational performance of the ARO, NS and deterministic approaches

		ARO	NS	Deterministic
	Mean time [s]	4.59	1.11	1.28
	Time range [s] 2.41-16.05 1.07-1.21	1.01-1.89
	Time SD [s]	2.77	0.03	0.25

Acknowledgements

This work was carried out as part of the innovation project SENSIBLE (Storage ENabled Sus-taInable energy for BuiLdings and communitiEs -www.h2020-projectsensible.eu), which has received funding from the European Comission under the Horizon 2020 Framework Programme grant agreement No 645963.

portfolio and offer ancillary services.


Carles Martí
email: carles.marti2@gmail.com

Leonardo Pacifici

Antonio Laganà

Cecilia Coletti

A quantum-classical study of the OH + H 2 reactive and inelastic collisions

Keywords: Quantum-Classical methods, reactive probabilities, state to state non reactive probabilities, vibrational adiabaticity

We carried out a study of OH + H 2 scattering using a quantum-classical method, treating quantally vibrations and classically both translations and rotations. The good agreement between the state specific quantum-classical reactive probabilities and the corresponding quantum ones prompted the extension of the study to state to state probabilities for non reactive vibrational energy exchange. The study showed that H 2 reactive dynamics depends on the vibrational excitation, while the non reactive one is mainly vibrationally adiabatic. On the contrary, OH reactive dynamics is not affected by its vibrational excitation, whereas the non reactive one might produce some pumping up to higher vibrational states.

Introduction

The OH ( 2 Π) + H 2 ( 1 Σ + g ) → H 2 O (X 1 A 1 ) + H( 2 S) reaction together with its isotopic variants is an important case study of four atom elementary process because of the role played in combustion and atmospheric chemistry [START_REF] Warnatz | Rate coefficients in the C/H/O system[END_REF][START_REF] Clary | Theory of Reactive Collisions at Low Temperatures[END_REF]. This reaction represents, as well, one of the simplest and prototypical cases of four body processes, due to the fact that three out of the four intervening atoms are hydrogens with the consequent simplicity of electronic structure calculation and fitting to a reliable Potential Energy Surface (PES). Moreover, the lightness of H exalts quantum effects and makes the comparison of experimental work (based on flash photolysis [START_REF] Ravishankara | Kinetic study of the reaction of hydroxyl with hydrogen and deuterium from 250 to 1050 K[END_REF][START_REF] Tully | Flash photolysis-resonance fluorescence kinetic study of the reactions OH + H 2 → H 2 O + H and OH + CH 4 → H 2 O + CH 3 from 298 to 1020 K[END_REF] and crossed beams [START_REF] Alagia | The dynamics of the reaction OH + D 2 → HOD + D : Crossed beam experiments and quantum mechanical scattering calculations on ab initio potential energy surfaces[END_REF][START_REF] Strazisar | Mode-Specific Energy Disposal in the Four-Atom Reaction OH + D 2 → HOD + D[END_REF]) with theoretical studies [START_REF] Zhang | Full-dimensional time-dependent treatment for diatom-diatom reactions: The H 2 + OH reaction[END_REF][START_REF] Balakrishnan | Integral cross sections and rate constants for the reaction OH + H 2 → H 2 O + H: A semiclassical wave packet approach[END_REF][START_REF] Sierra | The OH + D 2 → HOD + D angle and velocity distribution: quasi-classical trajectory calculations on the YZCL2 and WSLFH potential energy surfaces and comparison with experiments at E t = 0.28 eV[END_REF] particularly interesting.

As a matter of fact, for this system a certain number of ab initio PES's have been produced over the years [START_REF] Bradley | A Quasi-Classical Trajectory Study of Product Energy and Angular-Distributions in OH + H 2 /D 2[END_REF][START_REF] Wu | A new potential surface and quasiclassical trajectory study of H 2 O + H → OH + H 2[END_REF][START_REF] Yang | Ab initio potentialenergy surfaces for the reactions OH + H 2 → H 2 O + H[END_REF][START_REF] Chen | A global potential energy surface for the H 2 + OH → H 2 O + H reaction using neural networks[END_REF]. Full quantum techniques have also been applied to dynamical studies of this process both for J = 0 and for J > 0 [START_REF] Zhang | Full-dimensional time-dependent treatment for diatom-diatom reactions: The H 2 + OH reaction[END_REF][START_REF] Chen | A global potential energy surface for the H 2 + OH → H 2 O + H reaction using neural networks[END_REF][START_REF] Bhattacharya | Cross sections and rate constants for OH + H 2 reaction on three different potential energy surfaces for ro-vibrationally excited reagents[END_REF][START_REF] Zhao | State-to-state differential cross sections for a four-atom reaction: H 2 + OH H 2 O + H in full dimensions[END_REF] and for its isotopic variants [START_REF] Xiao | Experimental and theoretical differential cross sections for a four-atom reaction: HD + OH H2O + D[END_REF][START_REF] Zhao | State-to-state differential cross sections for D2 + OH D + DOH reaction: Influence of vibrational excitation of OH reactant[END_REF]. Full quantum calculations at J > 0 demand an increasingly heavier amount of both compute time and memory resources as the total angular momentum J increases, making them very costly. A way of carrying out a full dimensional study without sacrificing the possible quantum features of vibrational energy exchanges and maintaining, at the same time, reasonable time and resources expense is to adopt a mixed Quantum-Classical (QC) technique [START_REF] Balakrishnan | Integral cross sections and rate constants for the reaction OH + H 2 → H 2 O + H: A semiclassical wave packet approach[END_REF][START_REF] Balakrishnan | A mixed quantal/classical study of the reaction OH + H 2 → H 2 O + H[END_REF][START_REF] Billing | The Quantum Classical Theory[END_REF]. In the QC method the quantum and classical treatments deal separately with different degrees of freedom whose coupling is ensured via an effective Hamiltonian formulated in terms of the expectation value of the semiclassical Hamiltonian, computed using the quantal wavefunction.

As discussed in this paper the recent publication of an ab initio potential energy surface of the title system by Chen et al. [START_REF] Chen | A global potential energy surface for the H 2 + OH → H 2 O + H reaction using neural networks[END_REF] and of the related full quantum estimates of state specific probabilities has allowed us both to test the accuracy of the QC treatment [START_REF] Coletti | Quantum-classical calculation of cross sections and rate constants for the H 2 + CN → HCN + H reaction[END_REF] and to extend the study to the state to state exchange of vibrational energy. Accordingly, the paper is organized as follows: in Section 2 the adopted QC method is described; in Section 3 details on the performed calculations and the computational benefits obtained from the parallelization of the code are illustrated; in Section 4 results obtained for the title process and the rationalization of its dynamical features are discussed; in Section 5 some conclusions are drawn.

Method

A full description of the quantum-classical method and of the code used has been reported elsewhere [START_REF] Balakrishnan | Integral cross sections and rate constants for the reaction OH + H 2 → H 2 O + H: A semiclassical wave packet approach[END_REF][START_REF] Coletti | Quantum classical methods: a quantumclassical approach to diatom-diatom reactive scattering and VV energy transfer[END_REF] and here we outline only its application to the OH + H 2 system. The considered quantities are the detailed (state specific, for the reactive process, and state to state, for inelastic collisions) 40 properties for colliding diatom-diatom systems by solving the related time dependent Schrödinger equation as formulated below:

ih ∂Ψ ∂t = Ĥ Ψ, (1)
The nuclear arrangement is defined in a space fixed Jacobi coordinate system via the three vectors R 1 (O-H c ), R 2 (H a -H b ) and R 3 (OH c center of mass-H a H b center of mass) (see Fig. 1). In this coordinate system the Hamiltonian can be expressed in terms of the cartesian coordinates of the three i-th r vectors (collectively indicated as {x, y, z}), obtained by mass scaling the R vectors:

Ĥ = - 3 i=1 h2 2µ i ∂ 2 ∂x i 2 + ∂ 2 ∂y i 2 + ∂ 2 ∂z i 2 + V ({x, y, z}) (2)
where µ 1 and µ 2 are the reduced mass of diatoms OH c and H a H b , respectively and µ 3 is the overall reduced mass of the system:

µ 1 = m Hc m O m Hc + m O , (3)
µ 2 = m Ha m H b m Ha + m H b , (4)
µ 3 = (m Ha + m H b )(m Hc + m O ) m Ha + m H b + m Hc + m O = (m Ha + m H b )(m Hc + m O ) M (5)
Using the corresponding spherical polar representation x i , y i , z i can be expressed as:

x i = r i sin θ i cos φ i (6) y i = r i sin θ i sin φ i z i = r i cos θ i
where r i is the radial coordinate and θ i and φ i represent the associated polar and azimuthal angles, respectively. In these coordinates Eq.2 can be reformulated as:

Ĥ = - 3 i=1 h2 2µ i ∂ 2 ∂r i 2 + 1 r i 2 ∂ 2 ∂θ i 2 + cot θ i ∂ ∂θ i + 1 sin 2 θ i ∂ 2 ∂φ i 2 + V ({r, θ, φ}) (7)
where {r, θ, φ} = r 1 , θ 1 , φ 1 , r 2 , θ 2 , φ 2 , r 3 , θ 3 , φ 3 and Ĥ = r 1 r 2 r 3 Ĥ (r 1 r 2 r 3 ) -1 , 45 once that the wavefunction transformation Ψ = ψr 1 r 2 r 3 is adopted in order to eliminate first derivatives in r i . The hallmark of the method is the already mentioned different treatment of classical and quantal degrees of freedom (in the present case, vibrations are treated quantally while translations and rotations classically). For this reason, use has to be made of the semiclassical 50 Hamiltonian that, once replaced the momentum operators with their classical counterparts, takes the form (see eq. 6 of Ref. [START_REF] Coletti | Quantum classical methods: a quantumclassical approach to diatom-diatom reactive scattering and VV energy transfer[END_REF]):

ĤSC = TQ + 3 i=1 1 2µ i r 2 i p 2 θ i + 1 sin 2 θ i p 2 φ i + p 2 r 3 2µ 3 + V ({r, θ, φ}), (8)
where TQ represents the purely quantum kinetic energy operator defined as:

TQ = - 2 i=1 h2 2µ i ∂ 2 ∂r 2 i . (9)
In order to couple the different degrees of freedom one has to consider the effective semiclassical Hamiltonian expectation value:

H ef f SC = ψ| ĤSC |ψ ψ|ψ , (10)
where the brackets denote integration over the variables r 1 and r 2 . More in detail, H ef f SC can be written also as

H ef f SC = ψ| TQ |ψ ψ|ψ + p 2 r 3 2µ 3 + 1 2µ 3 r 2 3 p 2 θ 3 + 1 sin 2 θ 3 p 2 φ 3 + 2 i=1 ψ|1/r 2 i |ψ ψ|ψ 1 2µ i p 2 θ i + 1 sin 2 θ i p 2 φ i + ψ| V |ψ ψ|ψ . (11)
In order to follow the time evolution of the system, the time integration of classical degrees of freedom and the wavepacket propagation need to be performed simultaneously. For the latter one we can rewrite the time dependent Schrödinger equation by considering only the quantal part of the Hamiltonian:

ih ∂ψ ∂t = 2 i=1 - h2 2µ i ∂ 2 ∂r i 2 + 1 2µ i r i p 2 θ i + 1 sin 2 θ i p 2 φ i + V (r 1 , r 2 ; t) ψ. (12)
and integrate it using a grid of internuclear distance values. To this end a twodimensional discrete grid in the quantum variables r 1 and r 2 is adopted and the initial wavefunction is formulated as the product of two Morse oscillator wavefunctions representing the vibrations of the two molecules,

ψ(r 1 , r 2 , t = 0) = ϕ υ 1 (r 1 )ϕ υ 2 (r 2 ). ( 13
)
The time propagation of the wavepacket is performed using the split operator. On the contrary, the classical variables are integrated using the predictor corrector method. The reactive probability is computed when a trajectory is completed by evaluating the wavefunction flux before its absorption by an imaginary absorbing potential [START_REF] Neuhauser | The time-dependent Schroedinger equation: Application of absorbing boundary conditions[END_REF] located in the product channel, while the non reactive one is obtained by analysing the remaining non reactive wavepacket. To this end the final wavefunction is projected onto the reactants basis functions to evaluate the final population distribution of the vibrational levels.

Note that in this case the evaluation of reactive and inelastic probabilities is performed simultaneously. This is somewhat different from the mixed quantum classical method introduced by Billing [START_REF] Billing | Rate constants for vibrational transitions in diatom-diatom collisions[END_REF] and recently extended by Babikov (see [START_REF] Babikov | Recent Advances in Development and Applications of the Mixed Quantum/Classical Theory for Inelastic Scattering[END_REF] and references therein) which is focused on inelastic scattering.

Computational Details

For the calculations to be reliable, convergency tests aimed at optimizing the simulation parameters have been performed. The parameters obtained from the optimization procedure are summarized in Table 1. As to the extension of the internuclear distance grid, simulations at different upper and lower grid limits for both r OHc and r HaH b (r 1 and r 2 in the notation of the equations of Section 2) have been performed. Tests have been carried out at 0.3, 0.4, 0.5 and 0.6 Å to set the minimum value of r OHc (r OHc min ) and at 1.5, 2.0, 2.5, 3.0, 4.0 Å to set its maximum value (r OHc max ). Similarly for r HaH b min and r HaH b max tests have been performed at r HaH b equal to 0.15, 0.20, 0.3, 0.4 and 0.5 Å and at 3.0, 4.0 and 5.0 Å, respectively. The initial value of the diatoms centers of mass distance, r 3 , was also tested. However, it seems that the potential energy surface is accurate up to r 3 = 8 Å ca., as also suggested by the fact that fitting points in ref. [START_REF] Chen | A global potential energy surface for the H 2 + OH → H 2 O + H reaction using neural networks[END_REF] extend up to ca. 15 bohr. Thus, the value of r 3 = 7 Å was chosen as initial parameter for all simulations.

Different ∆t values have also been used (0.01, 0.02, 0.04, 0.08, and 0.1 fs) to the end of evaluating the convergence of the treatment within the integration time step. The number of trajectories has also been optimized at the lowest energy of 0.1 eV by checking the convergence of the results. Table 1 shows the adopted values of the various parameters. In order to gain computing efficiency we have parallelized the Fortran code using the shared memory paradigm in its OpenMP implementation. The parallel distribution of the workload was performed by assigning different trajectories to the various threads. As expected, this coarse grain parallelization leads to performances close to ideal (solid line) upon an in-105 crease of the number of processors (see Figure 2 for the related speedup plot when increasing the number of processors up to 8).

Results and Discussion

State specific reactive QC probabilities were calculated for total angular momentum J = 0 and compared with the full quantum ones of Chen et 110 al. [START_REF] Chen | A global potential energy surface for the H 2 + OH → H 2 O + H reaction using neural networks[END_REF] at different initial translational energies (Figure 3). The Figure shows the good agreement of the two results which differ for less than 10%, with the QC probabilities slightly overestimating the full quantum ones in the vibrational ground state for the two diatoms. A further qualitative analysis can be made by comparing QC results with the experimental ones. The latter show a marked increase of the reaction rate constant when H 2 (note that hereinafter we drop the labelling of hydrogen atoms in H 2 and OH because we always refer to the arrangement of the reactants.) is vibrationally excited [START_REF] Zellner | Vibrational Rate Enhancement in the Reaction OH + H 2 → H 2 O + H[END_REF] (specifically for v H 2 = 1 with respect to v H 2 = 0), whereas no significant change in the reaction rate is observed upon vibrational excitation of OH [START_REF] Light | The effect of vibrational excitation in the reactions of OH with H 2[END_REF] (for v OH = 1 with respect to v OH = 0). Indeed, our results confirm the higher efficiency of the initial vibrational excitation of the reactant H 2 in promoting the reaction (see Figure 4). The plots of the QC reactive probability (compare crosses (v OH = 0, v H 2 = 0) with points (v OH = 0, v H 2 = 1) and solid circles (v OH = 0, v H 2 = 2)) of Figure 4) show a significant increase of the reactive probability of the system with v H 2 , that is from 0.1 to 0.2 larger than that of the neighbour lower vibrational state in the whole investigated collisional energy range. In contrast, the vibrational excitation of OH from v OH = 0 to v OH = 1 and v OH = 2 has no appreciable effect on the value of the reactive probability of the system: as shown in Figure 4, the reactive probability computed for different values of the vibrational quantum numbers of OH are in this case hardly distinguishable from that of v OH = 0, in agreement with the experimental results of ref. [START_REF] Light | The effect of vibrational excitation in the reactions of OH with H 2[END_REF] and the theoretical ones of refs. [START_REF] Zhang | Full-dimensional time-dependent treatment for diatom-diatom reactions: The H 2 + OH reaction[END_REF][START_REF] Zhao | State-to-State Mode Specificity: Energy Sequestration and Flow Gated by Transition State[END_REF]. In particular, our results also show a very small increase in the reaction probabilities for v OH = 1, with respect to v OH = 0, at small energies (not appreciable in Fig. 4). This difference in the reaction probabilities tends to vanish as the energy increases. The good agreement between QC and full quantum approaches for v OH = v H 2 = 0 and, qualitatively, between QC results and experiments is a signif-icant test for the present approach, taking into account that the computational effort of these calculations is only slightly larger than quasi-classical methods and that relevant quantum effects (like the zero point energy of the diatoms) are included in the treatment.

Once proven that QC calculations can provide reliable results, we focused our investigation on the understanding of the reaction mechanism. For this reason, we have analyzed the trajectories associated with the QC calculations performed using the same initial conditions of the quantum study. The product H 2 O molecule is formed by the binding to OH of one of the hydrogens in which H 2 breaks. This is in agreement with the fact that, as shown in Fig. 4, an increase of the vibrational excitation of the hydrogen molecule (and, therefore, the stretching of the related bond) is indeed the key factor in promoting reactivity.

Typically, the trajectories associated with the initial conditions of the reactants mimicking those of quantum calculations follow two qualitatively different paths to reactions. The first one (see Fig. 5, upper panel) is consistent with a direct mechanisms in which H 2 dissociates and one of the dissociated H sticks on the OH molecule while the remaining hydrogen atom flies away. The second one (see Fig. 5, lower panel) is typical of a long lived intermediate complex forming mechanism that leads to more effective energy redistributions among the various degrees of freedom.

Our calculations also allowed to single out an effect that only full quantum and QC calculations can spot, while it might go undetected in pure quasiclassical trajectory calculations due to their poor discretization method. This effect is the vibrational excitation of the non reacted OH molecule as a result of the already commented formation of an intermediate long lived complex. This is clearly shown in Figure 5, lower panel, where the quantum expectation value of r OH after collision becomes sensibly larger.

The possibility of producing vibrationally excited OH after collision is confirmed by the calculations of the probabilities for vibrational energy exchange process. An accurate quantitative evaluation of vibration-to-vibration and vibration-to-translation energy transfer cross sections and rate constants would require an initial distance of the centres of mass of the diatoms much larger than 8 Å, which, as already mentioned, seems to be the largest distance at which the PES is still reliable. Indeed the long range part of the potential is known to play a crucial role in these non reactive events, particularly when resonant or near resonant processes are concerned. Because of this, here we will only make qualitative considerations about the diatoms vibrational population distribution after non reactive collisions. Figure 6 shows the QC state to state non reactive probabilities corre-

(v H 2 = 1, v OH = 0) → (v H 2 = 0, v OH = 1) than for (v H 2 = 0, v OH = 1) → (v H 2 = 1, v OH = 0
) process, and even more so for (v

H 2 = 2, v OH = 0) → (v H 2 = 1, v OH = 1) with respect to (v H 2 = 0, v OH = 2) → (v H 2 = 1, v OH = 1
) process. This might be an indication of a tendency to the depletion of excited H 2 and an increase of excited OH in the 195 final states. Furthermore, at higher translational energies vibrationally excited H 2 is more prone to react and, thus, is much less efficient in populating higher vibrational levels.

Conclusions and Future Perspectives

In this paper we use a Quantum Classical technique in order to perform a detailed state to state study of the title process and to obviate the heavy computational demand associated with the use of full quantum approaches. Furthermore, we exploit the fact that we have been able to gain almost one order of magnitude in speed by running in parallel on a cluster of 8 processors. This has been made possible both by the largely decoupled nature of the classical mechanics component of the QC code and by the compactness of its quantum one that makes the code highly distributable.

Our calculations led to an excellent agreement between full quantum and QC results. The QC calculations provide a clear evidence for the fact that there exist two different reactive mechanisms for the OH + H 2 → H 2 O + H reaction: the first one proceeding through a direct attack of the H atom originating from the dissociation of H 2 and a second one forming a long lived complex in which the two H atoms play a longer living interchange of energy among the different degrees of freedom. The two mechanisms also offer a rationale for understanding why the behavior of non reacted H 2 is mostly vibrationally adiabatic (the largest fraction of vibrationally excited H 2 ends up by reacting, as shown also by the fact that vibrational excitation of H 2 significantly increases reactivity), while the final vibrational distribution of non reacted OH shows a tendency to give higher vibrational energy distributions.

Preliminary calculations on internal energy redistribution processes upon non reactive collisions have also demonstrated that the QC code can be effectively used to investigate quantum effects in inelastic scattering. However, in order to get reliable results in this sense, the potential energy surface has to take into account long range effects, which are essential to accurately evaluate the exchange of vibrational quanta of energy.

To this end future studies, beside being extended to include larger collision energies and higher vibrational excitations and to calculate averaged quantities so to allow the computation of the relevant rate constants, will be directed on the investigation of how the nature of the potential, and particularly the introduction of a long range tail to the present PES, affects both the reactive and inelastic energy exchange processes.

Figure 1 :

 1 Figure 1: Sketch of the fixed Jacobi coordinates used for the OH c + H a H b system

  OHc 42 N r Ha H b 81 r OHc min , r OHc max , ∆r OHc ( Å) 0.40 2.5 0.05 r HaH b min , r HaH b max , ∆r HaH b ( of trajectories at each energy value 30,000 Table 1: Optimized parameters used for the simulations. N r OHc and N r Ha H b represent the number of grid points along r OHc and r HaH b coordinates, respectively; r 3 is the distance between the centers of mass of the diatoms.

Figure 2 :

 2 Figure 2: Measured Speedup (crosses) of the QC code using up to 8 cores. The solid line shows the ideal behaviour.

Figure 3 :

 3 Figure 3: State specific (J = L = 0, v OH = v H2 = 0, j OH = j H2 = 0) OH + H 2 → H 2 O + H reactive probability plotted as a function of the translational energy (QC results in red, full quantum results in green)

Figure 4 :

 4 Figure 4: State specific (L = 0, j OH = j H2 = 0) OH + H 2 → H 2 O + H QC reactive probability plotted as a function of the translational energy at different initial vibrational levels of the reactants.

Figure 5 :

 5 Figure 5: Plot of a direct (upper panel) and intermediate forming (lower panel) trajectories for J=0, v OH =0, v H2 =0. Classical values of r OHa (dashed-dotted line), r OH b (double dashed-dotted line) and quantum expectation values of {r HaH b } = ψ|r Ha H b |ψ ψ|ψ (dotted

Figure 6 :

 6 Figure 6: Non reactive state (v OH , v H2 ) to state (v OH , v H2 ) probabilities plotted as a function of one of the product vibrational level at E tr =0.1 eV, solid line, E tr =0.3 eV, dotted line, and E tr =0.6 eV, dashed line.

Acknowledgements Prof. D.H. Zhang and Dr. J. Chen are gratefully acknowledged for providing the potential energy surface. Financial support is acknowledged from the Phys4entry FP7/2007-2013 project (Contract 242311) and the EGI-Inspire project (Contract 261323). The European Grid Infrastructure (EGI) through the COMPCHEM Virtual Organization, the Italian CINECA computing centre and the Supercomputing Center for Education & Research (OSCER) at the University of Oklahoma (OU) are acknowledged for providing computing resources and services. CM thanks The European Joint Doctorate on TCCM project ITN-EJD-642294-Theoretical Chemistry and Computational Modeling. Partial financial support by the Fondazione Cassa Risparmio di Perugia (Project Code 2015.0331.021 Scientific and Technological Research) is gratefully acknowledged.