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

Alexis Gerossier

Andrea Michiorri

Georges Kariniotakis

Carlos Adrian Correa-Florez

Stochastic operation of home energy management systems including battery cycling

Keywords: Microgrids, energy storage, stochastic optimization, uncertainties, battery cycling, flexibility. Nomenclature Abbreviations BESS Battery Energy Storage System CSO Competitive Swarm Optimizer X 1, h, s, ..., X T, h, s ]

come

Introduction

The increasing penetration of decentralized renewable generation in the medium-and low-voltage grid is motivating development of new tools in order to face the challenges imposed by this new paradigm. These trends push even further, reaching the building and home level, and 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 load and PV production, and also considering battery aging.

Current research 1.Battery cycling in smart home applications

The aging of storage devices is a complex process that depends on chemical reactions with electrode interfaces, and the degradation of materials caused by cycling and aging of nonactive components [START_REF] Dutt | Life cycle analysis and recycling techniques of batteries used in renewable energy applications[END_REF][START_REF] Sarasketa-Zabala | Realistic lifetime prediction approach for Li-ion batteries[END_REF]. This process can be analyzed and modeled by tracking the cycling patterns, the respective Depth of Discharge (DOD) and the rate at which this process occurs [START_REF] Linden | Handbook of batteries[END_REF].

Some research has been published in recent years that includes this process in the HEMS operation. For example, [START_REF] Zheng | Smart households : Dispatch strategies and economic analysis of distributed energy storage for residential peak shaving[END_REF] evaluates the impacts of peak shaving by means of active demand and different storage technologies in a single household. To include cycling of the storage device, a set of values of the energy that can be cycled are predefined and analyzed.

In [START_REF] Bordin | A linear programming approach for battery degradation analysis and optimization in offgrid power systems with solar energy integration[END_REF], a degradation model is used to optimize the operation of an off-grid system with a single PV and battery. The proposed linear model identifies lower and higher State of Charge (SOC) and charging/discharging cycles, and assigns a linear cost.

The work in [START_REF] Davies | Energy storage control for a small community microgrid[END_REF] presents a model for a smart energy community in which storage depreciation is calculated based on a predefined lifespan of 10 years and 3000 cycles for a li-ion Battery Energy Storage System (BESS). The model then calculates a proportional cost with the net energy input. This calculation method disregards partial cycling of the BESS, and thus can lead to an underestimation of the actual depreciation cost. A model considering voltages and currents produced by different levels of DOD is analyzed in [START_REF] Darcovich | Control strategies and cycling demands for Li-ion storage batteries in residential micro-cogeneration systems[END_REF], for serving the purpose of managing resources in a residential microgeneration system composed of a single house with a PV-battery array.

A more detailed [START_REF] Goebel | Model-Based Dispatch Strategies for Lithium-Ion Battery Energy Storage Applied to Pay-as-Bid Markets for Secondary Reserve[END_REF] electro-thermal battery model is used for determining savings in the secondary reserve market for a system-operator owned BESS. This model includes a variable dispatch cost for batteries through parameter-fitting analysis, including terminal voltage, currents, temperature and SOC. Although this is a more detailed model of the internal interactions in the BESS, this approach would require parametrization for each storage unit analyzed.

Battery cycling in other power systems applications

The authors of [START_REF] Duggal | Short-Term Scheduling of Thermal Generators and Battery Storage With Depth of Discharge-Based Cost Model[END_REF] develop a short-term cost model for a utility-scale BESS, in order to solve a 24h resource-scheduling problem. The number of cycles and the DOD for a given time horizon are explicitly included in the optimization problem. Given that the relation of DOD and life cycles is nonlinear, the operation costs of the BESS are based on linearization by segments, and by assigning a charging cycle variable and constant cost. In this case, the resulting Mixed-Integer Linear Programming (MILP) is solved using a commercial solver. A Similar linearization-based modeling is used in [START_REF] Sarker | Optimal Operation and Services Scheduling for an Electric Vehicle Battery Swapping Station[END_REF] by approximating the slope of battery life as a function of the number of cycles but excluding the effects of DOD and solving the resulting MILP with a commercial solver. Paper [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] includes a vehicle battery degradation model consisting on piecewise linear approximation and including the effects of DOD. This research also uses a MILP solver to perform the optimization.

Following the same logic as the previous work, the research [START_REF] Koller | Defining a degradation cost function for optimal control of a battery energy storage system[END_REF] presents an explicit cost function that models battery degradation, which is used to implement a Model Predictive Control (MPC) peak shaving algorithm including a 1 MW BESS. The authors achieve the explicit formulation of the degradation costs by detecting transitions between charging/discharging and idle mode by state representation. In this way, they identify a quadratic cost function that captures cycling stress in terms of power, DOD and SOC. One advantage of this approach is the possibility of embedding a quadraticapproximated aging model into an optimization problem, without adding major complexity in terms of non-linear equations.

In [START_REF] Hemmati | Optimal design and operation of energy storage systems and generators in the network installed with wind turbines considering practical characteristics of storage units as design variable[END_REF], authors propose a model for wind-based network planning, and analyze the economic impacts of changing initial SOC and the non-linear inverse relation of DOD and life-cycles. Another interesting model which considers explicit cycling by counting state transition in the case of Hybrid Electric Vehicles (HEV), is presented in [START_REF] Wang | Integrated Energy Exchange Scheduling for Multimicrogrid System With Electric Vehicles[END_REF]. Quadratic explicit modeling is proposed in [START_REF] Mgel | Scheduling distributed energy storage units to provide multiple services under forecast error[END_REF] to calculate battery degradation costs. This expression does not consider state transitions and is only dependent on the power profile as a sum of a linear and quadratic function. An approach using a linear cost term is also used in [START_REF] Rahmani-Andebili | Cooperative distributed energy scheduling for smart homes applying stochastic model predictive control[END_REF] and [START_REF] Rahmani-Andebili | Stochastic, adaptive, and dynamic control of energy storage systems integrated with renewable energy sources for power loss minimization[END_REF], to model cycling costs of electric vehicles and BESS, respectively, by considering the state changes in the batteries and the initial battery price to account for degradation.

Cycling aging is considered in [START_REF] He | Optimal Bidding Strategy of Battery Storage in Power Markets Considering Performance-Based Regulation and Battery Cycle Life[END_REF] to complement a bidding strategy for Day-Ahead (DA), spinning reserve and regulation markets. By identifying local adjacent extreme energy points extracted from the SOC curve, the energy difference between these points is assumed as the DOD at which half cycle occurs. After this, with each DOD, an equivalent cycling is calculated for a complete day.

Depending on the final structure of the problem and the inclusion of binary variables, the mathematical optimization methods may vary. For example, commercial solvers are used in [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][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]. On the other hand, user developed techniques such as Stochastic Dynamic Programming (SDP) [START_REF] Abdulla | Optimal Operation of Energy Storage Systems Considering Forecasts and Battery Degradation[END_REF] and Generalized Benders Decomposition (GBD) [START_REF] Tan | Optimal Scheduling of Battery Charging Station Serving Electric Vehicles Based on Battery Swapping[END_REF] can also be used. The selection of a certain technique is not straightforward and depends on the approximations made to model the degradation of the battery, the complete mathematical model and the objectives that are pursued.

Rainflow Counting Algorithm for battery cycling

Another way of taking into account battery aging is by using the rainflow counting algorithm (RCA) [START_REF] Musallam | An Efficient Implementation of the Rainflow Counting Algorithm for Life Consumption Estimation[END_REF] originally proposed by Downing in [START_REF] Downing | Simple rainflow counting algorithms[END_REF] to calculate metal fatigue cycles. This method captures a more complete behavior of the aging process because it can include non-linear relation of equivalent cycling versus DOD. An inherent drawback the difficulty of embedding a mathematical representation into an optimization process.

RCA was recently used in [START_REF] Niu | Sizing and Coordinating Fast-and Slow-Response Energy Storage Systems to Mitigate Hourly Wind Power Variations[END_REF] for counting complete and incomplete cycling of batteries that provide fast and slow response to offset wind power variation. The use of RCA allows assessment of the lifetime depreciation for the provided services. It is important to note that this method is not embedded into an optimization model but is used after a specific SOC is obtained. A similar approach is presented in [START_REF] Ke | Control and Size Energy Storage Systems for Managing Energy Imbalance of Variable Generation Resources[END_REF] for battery life estimation with an additional step for equivalent life estimation due to incomplete cycles.

RCA is also used in [START_REF] Dragievi | Capacity Optimization of Renewable Energy Sources and Battery Storage in an Autonomous Telecommunication Facility[END_REF] for calculating the capacity loss of a battery acting as a backup for a telecommunication facility fed by other RES. RCA can be used to calculate the DOD occurring at specific full and half cycles; and the equivalent cost due to the capacity decrease is then added offline to the overall investment cost.

Authors in [START_REF] Li | Analysis of a new design of the hybrid energy storage system used in the residential m-CHP systems[END_REF] consider a hybrid storage system and a variation of RCA to quantify BESS lifetime extension. The single home system includes one micro CHP, one battery, superconducting magnetic energy storage, power converters and the demotic load. A set of rules, instead of an optimization problem, is defined to establish a control strategy to feed loads in the household.

Other recent research, attempts to mathematically mirror the state transition identification and counting of the RCA, to find a tractable and convex equivalent model, so as to embed the formulation and solve electricity related problems [START_REF] Xu | Factoring the Cycle Aging Cost of Batteries Participating in Electricity Markets[END_REF][START_REF] Xu | Optimal regulation response of batteries under cycle aging mechanisms[END_REF]. These contributions represent an important tool to overcome the difficulty of including the RCA logic in optimization problems in future research.

Thermal and electric energy storage coordination

Regarding management models for joint thermal and electric storage technologies at the residential level, approaches include the one 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 selfsufficiency. 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 in order to achieve minimum operation costs.

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 BESS, but does include thermal load and electro-thermal storage, which can generate or consume power.

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 Building Energy Management[END_REF], including electrical and thermal storage. This study concludes that thermal storage is crucial to reduce energy costs. However, it does not take into account cycling aging and points out that batteries might not be economical due to investment costs and short lifetime.

A recent paper [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) for optimal management of resources, 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 are presented, showing the importance of exploiting operational flexibilities when various interactions are analyzed.

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 PCC. Although it accounts for neither battery cycling nor uncertainties, this paper presents a thorough modeling of different energy sources and their interactions, and is tested on a system comprising 300 households.

Although this review is not exhaustive, it does illustrate the trends related to the topics in discussion and shows that the current literature leaves room for significant development regarding the integration of storage technologies in HEMS, including battery cycling.

About the present work

Continuing the track of the above-mentioned literature, this research intends to pursue the issue of flexibility related to BESS and Thermal Energy Storage (TES). Hence, 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 and load. These aspects are analyzed from the standpoint of an aggregator participating in DA and imbalance markets. In particular, to the best of the authors' knowledge, home level management of electric and thermal storage with load and RES uncertainties has not been addressed in combination with a metaheuristic that can capture the non-linearities of battery degradation by means of RCA. This proposed joint approach permits a more realistic overall modelling of the problem adapted to real implementation.

For the present work, an HEMS comprising data for the 25 real life households in Evora, Portugal, is presented as the 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 a part of the use case "Flexibility and demand side management in market participation". This case assumes 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 main contributions of this paper are the following: 1) Interactions of electrical and thermal storage are analyzed in the context of HEMS under uncertainty, and the advantages of coordinated operation are shown;

2) The battery's cycling aging cost is included by using a novel Storage Disaggregation Algorithm (SDA) based on Lagrangian relaxation and the RCA. This method is presented for the first time in the literature and allows handling the complex cycling behavior and reducing the optimization search space; and finally, 3) A decomposition strategy complemented by a competitive swarm optimizer is proposed. This approach allows solving separately the battery cycling problem solved by the SDA and the day-ahead scheduling of resources. This problem separation logic for scheduling resources in microgrids, has never been used by any previous research in order to solve the re-sulting two-stage stochastic non-linear/non-analytic optimization model.

In addition, the following assumpions were made to carry out the optimization process: 1) the power exchange at the points of common coupling does not jeopardize distribution network; 2) we suppose no market barriers for participation regarding minimum bid volumes; 3) the aggregator has communication and control capabilities with the devices at the home level. There exist the necessary IT and communication platform, so that the aggregator controls devices at the home level and decides over their set-points, 4) the proposed microgrid does not have market power, hence acts as a price taker; and 5) distribution network operation is part of DSO's responsibilities and is not part of the aggregators capabilities.

The present work is organized as follows: section 2 presents the mathematical formulation of the HEMS. Next, section 3 sets out the solution algorithm used for the energy management. The results obtained are given in section 4 and finally, conclusions are drawn in section 5.

Home Energy Management System Mathematical Model

The proposed system is composed of solar panels, li-ion batteries, heat storage devices, a connection to to the main grid and a number of households. 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.

In general, EWH input and electrical load during the 24h period can be met by the main grid, the solar panels, and the power injected from the batteries. The idea is to achieve a minimum operation cost by adjusting the setting of the devices in order to optimally manage resources. 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 secondary 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.

Wholesale electricity market

Electrical load and PV forecasts

The SENSIBLE project [37] explores integrating 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 features 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 predicted 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 were 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 methods described above and based on data collected from smart meters and Numerical Weather Predictions from ECMWF (European Centre for Medium-Range Weather Forecasts).

Scenario generation

The results of the quantile forecast are used to select a central value for both PV and load, by specifically using the median (quantile 50%) as this central forecast. To avoid defining arbitrary values of deviations from the central forecast to create scenarios, quantiles 10% and 90% are taken as the lower and upper bounds of forecast values. In this way, all central and deviated values of PV and load are combined to form a set of nine scenarios, representative of all potential combinations of minimum, maximum and central values according to realistic information from measurements and predictions and assuming uncorrelation between demand and PV forecasts. These nine scenarios (s1-s9), shown in table 1 are used as input data for the second stage formulation or recourse problem defined in the next section. It is important to emphasize that three groups can be identified from the nine scenarios described in 1: 1) a conservative group: formed by the scenarios featuring high load and low PV quantiles (i.e. s4, s7 and s8); 2) an equilibrated group: formed by scenarios with similar quantiles for both load and PV (i.e. s1, s5 and s9); and 3) and optimistic group: with scenarios featuring low demand and high PV quantiles (i.e. s2, s3 and s6).

The presence of these three groups of scenarios in the stochastic framework, allows a balance of optimality and robustness in the obtained solution, by combining the robustness introduced by the conservative group which generally leads to higher costs, but balanced by the potential low costs associated to optimistic scenarios, and complemented by the balance created resulting from the equilibrated scenarios. This set of nine scenarios, allow feasible and robust operation of the HEMS in the interval of extreme quantile realizations of demand and PV production, and allows exploitation of the probabilistic forecast methodology described in the previous subsection.

A visual representation of the aggregated net load (load minus PV) resulting from the nine scenarios for a typical day, can be seen in figure 3. The plot shows that curves s7 ( ), s4 (

) and s8 ( ) tend to have higher neat load during most of the 24 h horizon, which can be more clearly seen during the daytime when PV has values different from zero. In addition, s7 and s8 have also higher net load during nightime given that they are formed by using the 90% load quantile. This is a logic outcome, provided that these three scenarios (s4, s7 and s8) are indeed the ones classified in the conservative group. These scenarios will push the stochastic formulation to higher expected costs, in order to supply the required energy to be purchased in the market of injected from stored energy in the batteries.

Following the same logic, the optimistic scenarios can be identified as s2 (

), s6 ( ) and s3 ( ), with a clear tendency of low comparative values of net load, and even negative values for some time steps, meaning that available PV production is higher that the demand to be supplied. This situation could lead the stochastic formulation to take advantage of this available energy to store energy, supply load or sell back to the market.

Mathematical model

The proposed two-stage stochastic optimization model minimizes the expected HEMS day-ahead operation cost by scheduling the batteries' power charge and discharge, and the power injected into the EWHs, the TES pattern, and exchanges with the energy market.

Objective function

As shown in Eq. ( 2), the proposed two-stage stochastic optimization model minimizes the expected DA operation cost, in which the first-stage decision is associated with the DA purchase commitment (P g t ), and second-stage (recourse) expected costs are related to the import/exports (I - t,s /I + t,s ) imbalance and the expected cycling cost ( f cyc h (•)) for the batteries installed in each house h, which is a function of the SOC vector (X h,s ). min

T t=1 π t P g t + S s=1 p s T t=1 (µ - t I - t,s -µ + t I + t,s ) + N h=1 f cyc h (X h,s ) (2)
All of the terms in the objective function are linear, except the one related to the cycling. This term includes the corresponding nonlinearities associated with the chemical reactions occurring in the batteries due to temperature changes. This terms creates a complication in the model given the impossibility of accurately expressing this phenomenon. The way of tackling this issue is explained in section 3.1.

Load balance constraints

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

P g t + I - t,s -I + t,s + ∆t H h=1 P net t,s,h = 0, ∀t, ∀s (3)
It is important to note that the net power in each house, as per (4), considers the battery flows, the PV injection, the electrical load and the power required by the EWH, and this net power constraint is limited by the contracted power (constraint 5) of each customer with the distribution company at the point of physical connection. In addition, if a household is not provided with a TES device, the variable H becomes the same thermal load. Equation (3) creates a link between the first stage variables and the recourse decisions, by relating day-ahead commitment with the imbalance and setting of devices for each scenario.

Devices' constraints

Constraints ( 6) -( 13) describe the energy state for the BESS and the TES. Binary variables u t,h,s and v t,h,s are introduced to avoid charging and discharging batteries at the same time. Hence, constraints ( 8)-( 10) introduce a mixed integer characteristic into the model. Constraint ( 7) and ( 13) ensure the continuity of the storage devices from one day to another. It is important to mention that the term including R and C (thermal resistance and capacitance, respectively) in equation ( 12), 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].

X t,h,s = X t-1,h,s + η c ∆tP ch t-1,h,s -∆tP dch t-1,h,s /η d , ∀t, t 1, ∀h, ∀s (6)
X 1,h,s = X T,h,s , ∀h, ∀s (7)
u t,h,s + v t,h,s ≤ 1 (8) 0 ≤ P ch t,h,s ≤ Pch h • u t,h,s (9)
0 ≤ P dch t,h,s ≤ Pdch h • v t,h,s (10)
X h ≤ X t,h,s ≤ Xh (11) Y t,s,h = Y t-1,s,h + ∆tH t-1,s,h -Y t-1,s,h /R h C h -∆tQ t-1,s,h , ∀t, t 1, ∀s, ∀h (12)
Y 1,s,h = Y T,s,h , ∀s, ∀h (13)
Y t,h,s ≤ Y t,h,s ≤ Ȳt,h,s (14)
0 ≤ H t,h,s ≤ Ht,h,s (15)
Function f cyc h (•) in equation ( 2) is non analytic, given that the shape of this cost function can only be determined after ascertaining the SOC and the corresponding equivalent cycling, and then determining the DOD at which each equivalent cycle occurs. Hence it is difficult to model and express this non-analytic MINLP problem, and solve it with a commercial off-the-shelf software/solver.

The next section will explain how this phenomenon is treated so as to reshape and decompose the mathematical formulation in a more simple and tractable way.

Solution Methodology

Battery degradation costs

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 and 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 intrinsic behaviour and as a function of the DOD. 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 -k p ( 16
)
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%.

Cycle counting for a specific DOD is identified from local extreme points based on the SOC curve (X), and equivalent half or full cycles are defined. This is carried out following the logic of the Rainflow Counting Algorithm (RCA), which is explained in detail in [START_REF] Musallam | An Efficient Implementation of the Rainflow Counting Algorithm for Life Consumption Estimation[END_REF].

Once the number of cycles and their DODs for the 24h period have been calculated, an equivalent cycling cost is obtained according to the following expression:

f cyc (X) = C cyc (d j ) = j∈Ω L j C ini n 100 d k p j ( 17
)
where Ω is the set of DODs for the analyzed period, and C ini is the initial cost of the battery. The information of full or half cycles for each d j is given by L j , taking values of 1 or 0.5 respectively. The obtained C cyc for a specific SOC is the equivalent cost due to the battery's aging process, and should be added to the total dispatch cost.

Given the impossibility of mathematically embedding the transitions/cycling counting of the RCA that results in [START_REF] Mgel | Scheduling distributed energy storage units to provide multiple services under forecast error[END_REF], the problem formulated in (2)-( 15) is solved with a Competitive Swarm Optimizer (CSO), in which the decision variable is the batteries' SOC so that each time a solution is analyzed, the RCA can be applied to calculate an equivalent degradation cost and iteratively achieve an optimal solution. If SOCs are generated employing this metaheuristic logic, the problem can be decomposed into two subproblems: one of them analyzes batteries' SOC proposals and calculates the corresponding cycling aging equivalent cost; the second subproblem calculates the day-ahead cost due to the procurement of energy in the wholesale market and the settings for the remaining HEMS resources. Once the two subproblems are solved, the total cost can be obtained by simply adding both results (day-ahead and cycling cost).

Given that using a decision variable that explicitly contains the SOC for each of the batteries would yield a T ×S ×N dimensional search space, in this paper we perform a virtual aggregation of the batteries in the system so as to reduce the search space to T × S . With this approach, the dimension will be the same even if the system contains a high number of batteries. Given that an infinite number of possibilities exist to allocate individual SOCs, when an aggregated SOC is proposed, there should be a cost-efficient disaggregation/allocation of charge for each of the batteries, and this is carried out by applying a Storage Disaggregation Algorithm (SDA), as explained in the next section.

3.1.1. Cycling cost subproblem, the storage disaggregation algorithm When a certain SOC (X agg t,s ) is determined by the CSO logic, it has to be optimally allocated/dissagregated in each battery. This depends on each battery's cycling aging characteristics and the associated SOC. This subproblem can be mathematically formulated as follows:

min w = S s=1 N h=1 f cyc h (X h,s ) (18)
s.t.

X agg t,s = h X t,h,s , ∀t, ∀s (19)
When applying Lagrangian relaxation to this optimization problem, one multiplier appears for each time step in each scenario (λ t,s ), as shown in the function:

L = S s=1 N h=1 f cyc h (X h,s ) + S s=1 T t=1 λ t,s (X agg t,s - h X t,h,s ) (20)
After applying optimality conditions to the langrangian function (20) (i.e., derivative with respect to X t,h,s and λ t,s ), the equations obtained are:

∂L ∂X t,h,s = ∂ f cyc h (X h,s ) ∂X t,h,s + λ t,s = 0, ∀t, ∀h, ∀s (21)
∂L ∂λ t,s = X agg t,s - h X t,h,s = 0, ∀t, ∀s (22)
From equation ( 21) it is concluded that the derivative for each time step and scenario is battery invariant. This is a very important condition given that, in other words, the derivative of each battery should be the same for any given t and s:

∂ f cyc 1 (X 1,s ) ∂X t,1,s = • • • = ∂ f cyc h (X h,s ) ∂X t,h,s = • • • • • • = ∂ f cyc N (X N,s ) ∂X t,N,s = -λ t,s , ∀t, ∀s (23)
To calculate the derivative adn given that there is no analytic function to express f cyc h , numerical differentiation is used in order to iteratively find the X t,h,s that leads to [START_REF] Musallam | An Efficient Implementation of the Rainflow Counting Algorithm for Life Consumption Estimation[END_REF] while satisfying [START_REF] Tan | Optimal Scheduling of Battery Charging Station Serving Electric Vehicles Based on Battery Swapping[END_REF]:

λ (h) t,s ≈ f cyc h (X h,s ) -f cyc h (X h,s + ∆X t,h,s ) ∆X t,h,s (24)
Superindex (h) is introduced to denote that a multiplier λ (h) t,s should be calculated for each battery h, and the iterative process should correct the values X t,h,s until the multiplier is the same for all batteries (until λ (1) t,s = ... = λ (N) t,s = λ t,s ). This disaggregation algorithm is initialized by selecting a X t,h,s in such a way that eq. ( 19) (same as eq. 22) is met. After this, λ (h) t,s are calculated by using [START_REF] Downing | Simple rainflow counting algorithms[END_REF]. Given that X t,h,s , in each t and s, needs to be updated to achieve equal λ (h) t,s for all batteries, a deviation for each (h) is calculated by:

∆λ (h) t,s = λ (h) t,s -λt,s (25)
where,

λt,s = h λ (h) t,s /N (26)
Equation ( 25) measures the deviation of each battery's derivative with respect to the mean, hence, a simple heuristic rule is used to update SOC values according to:

X new t,h,s = X old t,h,s + φ • ∆λ (h) t,s (27)
where φ is a tuning parameter. Once the SOC is updated for each battery, eq. ( 24) is used again and the process is repeated until all deviations for each t and s are lower than a tolerance threshold. One very important feature of the presented method is that in this iterative process, optimality condition (22) (equivalent to constraint [START_REF] Rahmani-Andebili | Stochastic, adaptive, and dynamic control of energy storage systems integrated with renewable energy sources for power loss minimization[END_REF]) is always ensured, given that ∆λ (h) t,s = 0. This allows the aggregation of all batteries to equal the aggregated SOC at any time and scenario.

The outline of the proposed disaggregation algorithm is shown in figure 4. If any of the X t,h,s exceed the boundaries at a certain point of the SDA (to ensure constraint 11), a correction procedure is applied to maintain the feasible operation of devices. The correction algorithm is detailed in the appendix.

Day-ahead thermal subproblem

When a certain SOC is generated with the CSO and then disaggregated by means of the SDA, the problem in equations ( 2)-( 15) has to be reformulated. This subproblem is called the Thermal Subproblem, given that once an SOC is known, the remaining set points that need to be determined in each household are those associated with the TES.

First, P ch t,h,s and P dch t,h,s can be easily determined by applying the following rule: It is worth noting that constraint [START_REF] Bordin | A linear programming approach for battery degradation analysis and optimization in offgrid power systems with solar energy integration[END_REF] implies that both P ch t,h,s and P dch t,h,s have to be zero for the last time frame (t = T ) in order to match the SOC for the first time frame of the following day. Thus, equations ( 28) and ( 29) are valid for t : {1, ..., T -1}.

P ch t,h,s =        (X t,h,s -X t-1,h,s )/η c , if X t,h,s -X t-1,h,s > 0 0, otherwise (28)
P dch t,h,s =        (X t-1,h,s -X t,h,s )η d , if X t-1,h,s -X t,h,s ≥ 0 0, otherwise (29)
In addition, this simple rule eliminates the Mixed-Integer Non-Linear nature of the complete problem expressed in constraints ( 6)- [START_REF] Duggal | Short-Term Scheduling of Thermal Generators and Battery Storage With Depth of Discharge-Based Cost Model[END_REF], given that equations ( 28)-( 29) avoid charging and discharging at the same time, which is the reason why binary variables (u t,h,s and v t,h,s ) were necessary in the initial complete formulation.

Once P ch t,h,s and P dch t,h,s have been ascertained, a verification of these values is needed to determine if they are higher than the nominal power. If so, the analyzed SOC is infeasible and a penalization of this proposal is required. This penalization is calculated as follows:

Z = S s=1 N h=1 T t=1 [(P ch t,h,s -Pch h ) • H(P ch t,h,s -Pch h )- (P dch t,h,s -Pdch h ) • H(P dch t,h,s -Pdch h )] (30)
where H(•) denotes the Heaviside step function. From (30), whenever P ch t,h,s or P dch t,h,s are outside the boundaries, Y measures the proportion of the limit violation; and if feasible operation is achieved, then Z = 0.

Once the values of X t,h,s , P ch t,h,s , P dch t,h,s and Z are known, the model ( 2)-( 15) can be rewritten in the following way as a linear programming problem:

min z = T t=1 π t P g t + S s=1 p s T t=1 (µ - t I - t,s -µ + t I + t,s ) + βZ t (31) s.t.
Constraints : (3) -( 5), ( 12)

-(15) ( 32
)
where β is a penalization factor.

For each specific SOC proposal generated by the metaheuristic, the corresponding fitness function is calculated by adding the results from ( 18) and ( 31):

f itness = w + z ( 33
)
where w is obtained after using the SDA and z after solving (31)- [START_REF] Li | PV integration in Low-Voltage feeders with Demand Response[END_REF]. In order to iteratively find an SOC that returns the minimum cost, the CSO metaheuristic is used. This is because of its ability to handle real variables and its convenience for integrating the two subproblems described in the previous paragraphs.

It is worth noting that the two-stage characteristic is maintained in both subproblems, given that the variable for the CSO contains the SOC for each of the recourse decisions, and that the thermal subproblem is a two-stage linear problem containing the information for the remaining second stage variables.

Competitive Swarm Optimizer

The CSO algorithm is a metaheuristic optimization technique based on population behavior, first proposed in 2015 [START_REF] Cheng | A Competitive Swarm Optimizer for Large Scale Optimization[END_REF] and suitable for large optimization problems, following some of the principles of Particle Swarm Optimization. The algorithm assumes the existence of M particles. These particles move iteratively in an R-dimensional search space, where the i-th particle can be represented by a vector x i = (x i1 , ..., x iR ). The velocity of each particle is denoted by v i = (V i1 , ..., V iR ).

For this algorithm, M/2 pairs of particles are randomly formed, ensuring that each particle is selected only once. Each pair competes, and the particle with the best fitness function is designated as the winner (x w m (k)) and promoted into the next iteration. The loser (x l m (k)) has to update its position and velocity by learning from the winner.

For the m-th competition (m ∈ [0, M/2]), the loser's particle velocity and its position for iteration k+1 are updated according to [START_REF] Xu | Performance Analysis and Comparison on Energy Storage Devices for Smart Building Energy Management[END_REF] and ( 35):

V l m (k + 1) = r 1 V l m (k) + r 2 (x w m (k) -x l m (k)) +ψr 3 ( xm (k) -x l m (k)) ( 34
)
x l m (k + 1) = x l m (k) + V l m (k + 1) (35)
where xm (k) is the mean position of the particles in iteration k; and r 1 , r 2 and r 3 are three random vectors with uniform distribution in the range (0,1). Finally, ψ is a parameter that controls the influence of xm (k).

The algorithm ends when at least one of the following criteria is met: 1) a maximum number of iterations is achieved, or 2) a maximum number of iterations without improving the fitness function is achieved.

To adapt the optimization problem described in 2 to be solved by the CSO, the codification of the particles must be determined. In this case, each particle refers to an aggregated SOC for the batteries, which is composed of T × S values, associated with the time steps in each scenario for the day-ahead dispatch. Each particle must be assigned with a fitness function according to equation [START_REF] Good | Optimization under Uncertainty of Thermal Storage-Based Flexible Demand Response with Quantification of Residential Users' Discomfort[END_REF]. This means that for each particle in the swarm, the cycling and thermal subproblems described in subsections 3.1.1 and 3.2 need to be solved.

The initialization of the swarm is performed by randomly assigning SOC values to each specific time frame and scenario. The stop criterion used in this work is related to consecutive cycles without improving the fitness function. The complete outline of the proposed algorithm is shown in figure 5. 15) presents two main difficulties: 1) the integer nature of the variables that describe the charging/discharging mode of the batteries and 2) the difficulty of calculating cycling aging cost with a direct mathematical expression. To face these difficulties, the decomposition approach is used, in which a CSO generates proposals of aggregated SOC, and for each particle during the CSO evolution, two subproblems are solved: 1) a cycling subproblem, and 2) a thermal subproblem. The cycling subproblem takes as input the aggregated SOC of each particle and assigns a SOC for each of the batteries by using the SDA explained in subsection 3.1.1. When the dissagregated SOC is found, this is used as known data to determine the remaining operation point of the HEMS with the thermal subproblem.

This way, the iterative process of the CSO allows to find a solution that contains all the information regarding device settings and energy commitment with the market, to achieve minimum operation costs.

Complete MINLP

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, there are 25 PV panels, 16 BESS and 15 EWHs. The testcase is composed only by the HEMS and the control capabilities that an aggregator has on the device settings. The operation and control of the MV and LV networks is carried out by the DSO and are not part of 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].

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.

The charging and discharging efficiency of the batteres is assumed to be 95%. 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].

HEMS in LV Network

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].

Performance of the SDA

Before using the complete algorithm as outlined in Fig. 5, tests are carried out to analyze the adequacy of the algorithm proposed in section 3.1.1 to optimally allocate the individual SOC for each battery. Two tests are performed: first, 100 randomly generated aggregated SOC patterns are subject to the disaggregation process with random initialization, in order to determine the convergence capabilities by analyzing the evolution of ∆λ max ; second, a specific random SOC is selected and subject to 100 tests to prove that the algorithm converges to close values, by measuring the evolution of the total cycling cost.

The converge criterion is set at |∆λ max | ≤ 0.001. As seen in figure 9, for all the cases run, the algorithm reaches low values of |∆λ max | after a few iterations. In addition, the convergence is achieved in 32 iterations on average for the 100 cases.

When one particular aggregated SOC is selected to be analyzed multiple times under different initial conditions, it is ex- pected that the minimum total cycling cost obtained will be the same or at least close for each simulation, given the non-linear characteristic of the cycling function. Figure 8 shows the evolution of the total cycling cost for one particular aggregated SOC under 100 different initial conditions. The minimum and maximum obtained cycling costs are 0.96 e and 1.02 e , respectively; and the mean is 0.98. This shows that the algorithm is robust towards different initialization values and that, for the specific case of the unit e , only one decimal place is enough to express currency in real-life applications, without losing sensitive information for DA electrical markets.

Comparative analysis of flexibilities under deterministic approach

To analyze the impacts on the battery operation when different sources of flexibility are taken into account, a deterministic model is run using quantile 50% (median) forecasts for PV and load, during the complete month of November 2015, given the availability of the data. All presented values correspond to a daily average resulting from running the 30 independent DA models for each day. Considering only PV injection and no battery or thermal storage, results in an average cost of e 7.0, which we calle the base case, for comparison purposes. Different cases are tested to assess how each technology contributes in reducing costs, as shown in the table below. 2)- [START_REF] Hemmati | Optimal design and operation of energy storage systems and generators in the network installed with wind turbines considering practical characteristics of storage units as design variable[END_REF], which considers the cycling equivalent cost within the optimization model. This approach allows to determine device settings and energy exchanges such that the total operation cost (energy purchases + equivalent cycling aging) is minimized at once. Concretely, case 1 disregards the possibility of thermal storage and only considers electrochemical storage. When this independent control of the thermal storage is not allowed, a 2.15% cost reduction is achieved by using only the battery to minimize the cost. When thermal and power storage are permitted (case 2) to be controlled independently (but coordinated), a greater cost reduction is obtained, showing that including thermal storage as another flexibility leads to a decrease in the cost. Case 2 shows a reduction of 18.7 % with respect to the base case.

The cost evolution for the complete month for all of the proposed cases can be found in figure 10. The cost evolution in the figure shows that not only are the average values (in table 2) for case 2 lower, but that improvements are achieved for each of the analyzed days. The costs savings compared to the base case, range from 14.0% to 32.5% with a median of 18.5%.

In order to demonstrate the value of taking into account the cycling cost embedded in the model, (as done for cases 1 and 2), two more cases are proposed, in which the optimization model only attempts to minimize the energy purchased (first two terms in objective function ( 2)), and the battery cycling cost is disregarded and only calculated afterwards with the obtained SOC. This analysis corresponds to cases 3 and 4. Once again, it is shown that including TES adds flexibility to the model and allows cost reduction. However, these cases show that when cycling is not explicitly considered in the optimization model, the batteries are subject to deeper cycling, resulting in higher cycling aging and leading to higher total costs when compared to cases 1 and 2. The most critical case corresponds to an operation cost of 8.6 e (case 3). It can be concluded that if cycling cost is ignored in the optimization model, batteries can cycle without any constraint of frequency or depth, which results in suboptimal operating points for the aggregator. ) and 4 ( ) present deeper and more frequent cycling of batteries than cases 1 ( ) and 2 ( ). The coordinated scheduling of BESS and TES with the complete optimization model (case 2) represents the best improvement with respect to the base case given the efficient scheduling of batteries and the full exploitation of BESS and TES capabilities. The energy purchase for this deterministic case is shown in figure 12 ( ) and it can be seen that around noon the energy requirements are minimized due to PV availability; on the other hand, purchases increase during late night and morning hours. Figure 12: Energy purchased on the energy market for the deterministic ( ) and stochastic case ( )

Stochastic management of resources

When the complete two-stage stochastic optimization model is solved for each of the 30 days, the first-stage variable obtained determines the day-ahead purchase commitment on the wholesale market. This two-stage model is run by defining the nine second-stage scenarios as explained in section 2.2.

The average expected daily cost of the stochastic solution (SS) is 6.3 e and the associated cost for each of the 30 days is presented in figure 13 ( ). This cost is the result of the DA commitment and the expected imbalance settling cost plus the expected battery cycling cost, for each of the nine scenarios, as per equation 2. To test the adequacy of the presented stochastic formulation, VS S (Value of Stochastic Solution) index is used [START_REF] Birge | Introduction to Stochastic Programming[END_REF].

To calculate the VS S , which measures the cost of ignoring uncertainty to make a decision, a quantity called the expected value problem (EVV) needs to be determined. The EVV consists in solving the stochastic problem using the second stage containing all the necessary scenario information, but fixing the first-stage variables with the result of the deterministic solution. This allows us to find a second-stage optimal solution as a function of the first-stage variables obtained using the central forecast (quantile 50%). Then, the index is calculated by: VS S = EVV -S S . The value of EVV for the complete month is shown in figure 13 ( ). The VSS is positive each day, as shown in curve (

). Average values VS S mean and VS S % mean are afterwards calculated using:

VS S mean = 1/30 30 1 (EVV i -S S i ) (36)
VS S % mean = 1/30

30 1 EVV i -S S i EVV i • 100% (37)
As a result, VS S mean = 0.38 e . By using equation ( 37), the stochastic approach allows a reduction of operating costs by VS S % mean = 5.8% on average through the analyzed month, which is also the mean value of the curve (

). The VSS gives an idea of how well the optimization under uncertainty performs. In this case, the two-stage stochastic optimization represents the best option, given that it allows a reduction of the expected average cost when compared to the deterministic approach for taking the DA purchase decisions (given by EVV).

The boxplot distribution created by the 30 day-ahead purchase commitments in each time frame, is shown in figure 14. From the figure it can be seen that interaction with the grid decreases in the time frames around noon, given the available

Analysis of a single day

for the sake of example, and taking the same typical day used for the deterministic case, first-stage purchase commitment is shown in figure 12 ( ). From the figure it can be seen that interaction with the grid decreases in the time frames around noon, given the available PV. In addition, interaction is also zero for the 19, 20 and 23 time periods; in all time frames, this energy procurement is lower than the purchase commitment for the deterministic case (

). Given that this variable is fixed for all scenarios, secondstage imbalances allow feasible operation of the system in such a way that the expected cost is minimized. For this study case, imbalance is required at different time steps and for all of the second-stage scenarios, as shown in figure 15. The boxplots show that additional energy (negative imbalance) must be purchased in order to overcome shortages for the realizations of the different scenarios. In particular, scenarios 1, 4 and 7 present higher median and third quartile. These imbalance needs in fact, correspond to the scenarios with the lowest levels of PV production, so the aggregator has to purchase additional energy from the wholesale market to offset the energy imbalance.

Regarding the accumulated SOC for all 16 batteries in the system, it can be seen from figure 16 that there is a pattern for all nine scenarios and that storage generally increases around 5h and 16h. It is interesting to point out that scenarios 4 and 7

tend to have higher SOC just before noon. This is explained by the fact that these scenarios feature low levels of PV, hence the energy discharge is treated in a more conservative way. The lowest degradation value is related to scenario 1 (

), which has minimum levels of both load and PV. From the figure it can be observed that the SOC generally exhibits values in between those in the scenarios, thus avoiding sudden changes in stored energy. In this particular case, the found VS S = 0.45 e , EVV = 7.85 e and the expected stochastic cost is 7.4 e . This means that, on this particular day, the stochastic approach allows an average reduction of operational costs of 5.7% compared to a decision made by using the median as the expected value for load and PV.

4.5. Computational remarks 4.5.1. Computational performance of the CSO for the 25household testbed To test the performance of the CSO, several runs were carried out for a specific day in order to determine the quality and consistency of the obtained solutions. Tests consist in running the algorithm 20 times for different population sizes, setting the stop criteria to 30 iterations without improving the best solution. In addition, the maximum number of iterations is set to 300. In these conditions, the CSO is run, and after each simulation the information saved is: best solution and computational time in achieving the best solution. In addition, the performance of the implemented CSO is compared with the performance of a Particle Swarm Optimization (PSO) algorithm in order to have another swarm based metaheustic comparison. The values of the parameters used for the PSO are 0.9, 2.0 and 2.0, for inertia, cognitive and social parameters, respectively.

Table 3 shows the results for both CSO and PSO and the different population sizes. It can be observed that the CSO presents lower computational times for solutions close to those ones obtained with the PSO. Although the average best solutions obtained with the PSO are lower than those obtained with CSO, the computational times may be prohibitive for day-ahead decisions, given the need to run the algorithm for several hours. In contrast, the CSO achieves good quality and stable solutions in the range of 6500-12000s for population sizes of 30 and 50 particles.

The faster evolution of the CSO is explained by the constant update of the loser's velocity and the permanence in the swarm of the best paired particles, which allows good quality solutions to be maintained in each iteration and promoted into the next generational cycles.

Performance of the CSO for larger test systems

To test the performance of our algorithm for a larger number of households to be aggregated, four additional test cases with 50, 100, 150 and 226 households (referenced as A-50, A-100, A-150 and A-226 respectively) that belong to the same neighbourhood in Evora, Portugal were included. All of the buildings are within a circle of 250 m radius, mainly residential, with some restaurants and stores.Although PV, ESS and EWH in the real life are only installed in the 25 households of the original testbed, we mirrored the distribution of resources as outlined in figure 7 and we assigned a proportional number of resources in each case, based on the original 25-household HEMS.

For this performance analysis, the CSO was run 20 times for each test system with the following parameters: 30-particle population, ψ = 0.3 and stop criteria set to 30 iterations without improving the best solution and the maximum number of iterations is set to 300. The results are condensed in the following table: As expected, the results show that the computational time increases when the number of households to aggregate is larger. However, the computational times for the larger test system (A-226 households) remains reasonable for day-ahead decision making with an average of 11523 s, and ranging from 5590 to 14632 s. This is an important result, provided that it demonstrates that the algorithm can deal with larger aggregation without leading to prohibitive computational burden. Negative results of the objective function indicate that management of larger number of resources leads to increased profits.

There are two main reasons to explain why the presented approach can withstand the proposed test systems within reasonable times: 1) the size of each particle is determined by the time step and the number of scenarios, but is independent of the number of batteries in the system. This allows to keep the search space of the CSO invariant and battery-independent; 2) the SDA returns the charging/discharging pattern for each battery; this information is used by the thermal subproblem, with the advantage that this subproblem is a linear programming problem, thus avoiding explosion of problem size and binary variables.

Conclusions

A methodology for aggregation and management of resources in HEMS is proposed, in which the coordination of thermal and electrical storage shows that reduced operation costs can be achieved. In addition, when considering the cycling aging cost there is a change in the devices settings to avoid deeper and frequent cycling, while maintaining operation at a minimum cost.

The proposed scheme, which can also be used to analyze other sources of flexibilities, configurations and services, shows that taking cycling into account also avoids operating batteries close to maximum/minimum values, which can lead to employing flexibility to match variations in load or available resources in real-time operation.

The two-stage stochastic approach, shows the importance of taking into account uncertainties arising from PV and load, in order to avoid higher expected operation costs in comparison to deterministic approaches. In our case, the benefit of using a stochastic approach is demonstrated by means of the VSS. The average monthly reduction in operation cost measured by the VSS is 5.8%, and the interval of reduction ranges from 3.6% to 9.0%. These results show that considering a stochastic model joint with battery cycling leads to savings when compared to deterministic approaches.

The same decomposition presented in this paper, can be used not only in HEMS, but in other battery based applications. 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. In all cases, despite the agent, the ownership or the battery size, the decomposition approach will still be valid.

For larger aggregation of resources (i.e. thousands or tens of thousands) other solution techniques can be explored to decrease search space and avoid the computational burden of scenario based approaches, i.e. robust optimization and explicit modeling of the cycling aging by linearization.

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

Figure 1 :

 1 Figure 1: Schematic diagram of the proposed HEMS

Figure 3 :

 3 Figure 3: Normalized neat load per scenario

Figure 4 :

 4 Figure 4: Flowchart of the SDA

Figure 5 :

 5 Figure 5: Flowchart of the CSO algorithm used

Figure 6 :

 6 Figure 6: Scheme of the decomposition approach

Figure 7 :Figure 8 :

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

Figure 9 :

 9 Figure 9: Convergence evolution the disaggregation algorithm under 100 random cases

Figure 10 :

 10 Figure 10: Evolution of daily cost for the deterministic model under the analyzed cases and for the whole month of November

Figure 11 :

 11 Figure 11: Aggregated SOC for deterministic cases 1 ( ), 2 ( ), 3 () and 4 ( )

Figure 13 :

 13 Figure 13: Evolution of daily cost for the deterministic case under the analyzed cases and for the whole month of November

Figure 14 :

 14 Figure 14: Power purchase boxplot for each DA timeframe

Figure 15 :

 15 Figure 15: Imbalance boxplot for each scenario in the stochastic optimization

Figure 16 :

 16 Figure 16: Aggregated SOC for each scenario in the stochastic approach

Table 1 :

 1 Scenarios included in the stochastic scheme An example for a typical day of the normalized aggregated values for PV and load, containing the central, upper and lower values is shown in figure 2.

		Scenario	s1	s2	s3	s4	s5	s6	s7	s8	s9
		Load quantile	10	10 10 50 50	50 90 90 90
		PV quantile	10	50 90 10 50	90 10 50 90
		1.5	load 10%		load 50%		load 90%
			PV 10%		PV 50%		PV 90%
	p.u.	1
		0.5
		0	5	10		15		20
					time[h]

Figure 2: Normalized load and PV curves for each of the quantiles needed for scenario generation

Table 2 :

 2 Cost of resource management for the deterministic scheme. Average daily costs. *Includes cycling effects in the optimization. **Calculated after

	the optimization
		Used flexibility	Equiv. cost	Total operation	Improv.
	Case	BESS	TES	cycling [e ]	cost [e ]	w.r.t base case
	Base	No	No	N.A.	7.0	N.A.
	1*	Yes	No	0.1	6.8	2.15%
	2*	Yes	Yes	0.1	5.7	18.7%
	3	Yes	No	3.8**	8.6	-22.7%
	4	Yes	Yes	3.7**	8.4	-20.9%
	Cases 1 and 2 in table 2 correspond to the problem detailed in
	equations (

Table 3 :

 3 Computational performance of the CSO and PSO for one typical day

	in November
			CSO (ψ = 0.3)		PSO
	Population 30	50	80	30	50	80
	Best solu-	5.29	5.27	5.26	5.24	5.24	5.22
	tion [e ]
	Best sol.	5.29-	5.27-	5.26-	5.24-	5.24-	5.22-
	range [e ]	5.37	5.31	5.34	5.31	5.27	5.26
	Best sol.	5.33	5.30	5.28	5.28	5.26	5.24
	mean [e ]
	Best sol.	0.02	0.02	0.02	0.02	0.01	0.01
	SD [e ]
	Mean time	6867	11103	17194	13455	25176	41351
	[s]
	Time	6571-	9750-	16121-	12750-	23574-	38687-
	range [s]	7262	12087	18547	14280	26330	43111
	Time SD	217	662	765	538	947	1624
	[s]

Table 4 :

 4 Computational Performance for different test systems Test system Best sol. mean (SD) [e ] Mean time (SD) [s]

	A-50	-1.72 (0.05)	7811 (2053)
	A-100	-3.36 (0.04)	9767 (2941)
	A-150	-5.05 (0.05)	10372 (3809)
	A-226	-7.62 (0.08)	11523 (2794)

Acknowledgements

This work was carried out as part of the innovation project SENSIBLE (Storage ENabled SustaInable 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.

Appendix A.

When performing the SDA and after updating the variables by using [START_REF] Dragievi | Capacity Optimization of Renewable Energy Sources and Battery Storage in an Autonomous Telecommunication Facility[END_REF], some of the variables may exceed the boundaries. If this happens, a correction procedure involves each multiplier λ ( j) associated to each violating variable x j . Let Λ Y be the set that contains the violating variables, Λ N the set of the non-violating variables, and x * j indicates the minimum or the maximum allowed value for x j . For each variable in Λ Y , a correction that complies with the following equation needs to be carried out:

Reorganizing the terms we obtain the following:

If all of the equations in Λ Y are written, a linear system of the type A • λ new = b is generated, where A is a square matrix, b is a column vector and λ new is the vector with the multipliers in Λ Y that have to be corrected. To solve this system, the elements of A are calculated as follows:

and,

where A j j are the elements in the diagonal and A ji the elements outside the main diagonal. Elements of b, are calculated with the following expression:

After the linear system is solved, a set of multipliers is obtained such that the update of x j , ∀ j ∈ Λ Y is not out of bounds.