Development of a Framework for Activation of Aggregator Led Flexibility

This paper presents a novel framework architecture for an online, real-time flexibility assessment and activation platform targeted at unlocking the flexibility potential of commercial buildings and smaller industrial sites, thereby enabling greater levels of renewable grid integration. Renewable integration targets in Europe of up to 40% of power generation from renewable sources by 2030 and over 90% by 2050 aim to decarbonize the electrical grid and increase electrification of transport, industry, and buildings. As renewable integration targets increase, participation in flexibility programs will be required from a much greater range of buildings and sites to balance grids hosting high levels of renewable generation. In this paper, an online implementation of a standardized flexibility assessment methodology, previously developed for offline contract negotiations between stakeholders, is modified to automate the assessment. The automated assessment is then linked to an aggregator-based multi-building or site optimization stage, enabling increased participation of multiple buildings and sites. To implement the assessment, models for individual flexible systems were reviewed, selected, and adapted, including physics-based, data-driven, and grey-box models. A review of optimization for flexibility found mixed-integer linear programming to be the optimal approach for the selection of flexible systems for demand response events.


Introduction
The electrical grid is changing, primarily in response to the need for balancing increasing amounts of renewable generation. Renewable integration targets in Europe of 38 to 40% [1] of power generation from renewables by 2030 and 90% by 2050 [2] are being used as the means to decarbonize the grid and increase electrification of transport, industry, and buildings. Energy and power flexibility provided by industrial sites and a range of building types [3] is more cost-effective for grid management than strategies such as curtailment or gas-fired generation [4]. Current practice in power systems demands response programs is to utilize flexibility from a few large industrial sites [5]. Enabling renewable generation hosting capacities of over 65% [6] requires participation from a much wider spectrum of buildings such as offices, retail, apartment blocks, and single residential buildings [7] as well as industrial sites. These buildings and sites will be required to provide greater numbers of flexible systems and deeper ranges of flexibility [8].
Aggregators provide energy and power balancing services by acting as an intermediary between a contracting authority, typically a grid utility, and buildings or sites connected to the electricity grid [9]. The aggregator puts together a portfolio of sites to meet the minimum power or energy participation criteria set by the contracting authority to provide flexible services such as demand response (DR) to the grid [10]. Power and energy flexibility of grid-connected loads, generation, and storage has become increasingly

Framework Structure
An aggregator-based framework for automated implementation of the flexibility assessment methodology is shown in Figure 1. The system is to operate in an online realtime manner and be hosted in the cloud. When a demand response request is received by the aggregator (denoted AGG in Figure 1) from a grid operator, the software automatically verifies flexibility contracts with the buildings in its portfolio to determine which are valid for the particular demand response program the request relates to. Smart contracts developed for grid energy trading [19] or peer-to-peer energy trading [20] modified for flexibility may be implemented for the contract verification stage. The cloud-based platform then contacts the on-site automation system in the relevant buildings or sites to read data from the systems in the buildings. On-site automation systems in buildings are typically a form of building management system (BMS), whereas, in industrial sites, a supervisory control and data acquisition (SCADA) system may be installed. The aggregator-based framework then implements "Step 1: Identification of flexible sources, loads, storage and generation" from the standardized flexibility assessment methodology [16].
Step 1 is applied to the data for each of the buildings or sites individually. "Step 2: Flexibility Characterization" is conducted in an automated manner to quantify the available flexibility at each building or site. "Step 3: Scenario Modelling" then models the systems and outputs the available flexibility defined for the specific time period of the demand response request for each building or site. The available flexibility is then input into the aggregator multibuilding optimization algorithm, which selects the most appropriate buildings, sites, and sources for the event. The aggregator platform decides to accept or reject the flexibility available and a write signal to actuate sources is sent to the BMS systems in each selected building or site.
gies 2021, 14, x FOR PEER REVIEW 3 of based framework then implements " Step 1: Identification of flexible sources, loads, st age and generation" from the standardized flexibility assessment methodology [16]. S 1 is applied to the data for each of the buildings or sites individually. "Step 2: Flexibil Characterization" is conducted in an automated manner to quantify the available flexib ity at each building or site. "Step 3: Scenario Modelling" then models the systems a outputs the available flexibility defined for the specific time period of the demand sponse request for each building or site. The available flexibility is then input into aggregator multi-building optimization algorithm, which selects the most appropri buildings, sites, and sources for the event. The aggregator platform decides to accept reject the flexibility available and a write signal to actuate sources is sent to the BMS s tems in each selected building or site. Making buildings truly smart for flexibly may be more cost-effective if cloud-bas solutions are utilized. Connecting multiples of buildings to a single cloud-based platfo that implements the flexibility assessment methodology proposed in this work may acc erate the rollout of flexibility for commercial and residential buildings as it will redu cost by (a) achieving economies of scale and (b) minimizing modifications to existing BM systems. Economies of scale [5] would enable a single platform to compute flexibility si ultaneously for many buildings or sites. Minimizing expensive and time-consuming mo ifications to existing BMS systems in the buildings would reduce costs [18]. These wou enable online real-time flexibility quantification for both buildings and aggregators. ternatively, a combination of cloud and edge-based [21] implementation may be deploy to increase the computational efficiency and response time of flexible systems.

DR request from grid
Step1:

A G G
Step 2: Flexibility Characterisation

Online Realtime Automatic
Step 3:Scenario Modelling Making buildings truly smart for flexibly may be more cost-effective if cloud-based solutions are utilized. Connecting multiples of buildings to a single cloud-based platform that implements the flexibility assessment methodology proposed in this work may accelerate the rollout of flexibility for commercial and residential buildings as it will reduce cost by (a) achieving economies of scale and (b) minimizing modifications to existing BMS systems. Economies of scale [5] would enable a single platform to compute flexibility simultaneously for many buildings or sites. Minimizing expensive and time-consuming modifications to existing BMS systems in the buildings would reduce costs [18]. These would enable online real-time flexibility quantification for both buildings and aggregators. Alternatively, a combination of cloud and edge-based [21] implementation may be deployed to increase the computational efficiency and response time of flexible systems.

Use Cases
Practical examples of use cases that may be implemented using the framework are detailed below.

Dynamic Pricing
Under the EU's 2019 Electricity Directive, all member states will be required to provide dynamic pricing for customers with a smart meter [22]. Dynamic or real-time pricing, whereby the electricity tariff may change frequently based on grid balancing requirements due to renewable generation or load fluctuations, will create scope for active optimization and management of large numbers of small flexible resources [23]. This may be implemented with or without contracts with grid operators or regulators. Flexibility programs based around dynamic or real-time pricing already exist for large commercial customers in a number of countries. For example, in Ireland, a market-based program that requires the building or site to respond to a grid request intra-day within a short timeframe is the short-term active response (STAR) program [10].

Peak Shaving
Peak shaving is a price-based program whereby the contracted building or site agrees to reduce its load at the same time every day during peak demand. It is the most widely used demand response service globally, with Ireland [10], the US [24], and China [25] including it in their demand-side services. In Europe, the demand peak occurs during the winter [10] due to heating loads, where in the US, the peak occurs during the summer due to cooling loads [24]. With peak shaving, the building operator knows the price and time schedule, sometimes more than a year in advance.

CO 2 Minimization
CO 2 based demand response signals have been proposed as a mechanism to increase electricity use or reduce consumption, in times of high or low renewable generation on the grid [26]. An alternative to, or in combination with, price-based signals, they may be utilized in future markets where hourly or real-time generation emissions are available. Motivations for adoption include enhancement of carbon reduction strategies by businesses or citizens who wish to minimize their climate impacts, or for grid operators to maximize renewable generation consumption.

PV Power Smoothing
PV power smoothing is used to dampen fluctuations in power generation caused PV variability by using battery storage to compensate for the peaks and troughs [27]. It requires storage coupled with PV and may be activated at the request of the grid operator or may be implemented continuously. The objective of this use case is to smooth PV peak production by storing excess renewable electricity generated or discharging storage to compensate for sudden reductions in output. With the aggregator-based platform, the PV installation and the storage may be in different buildings or sites but coordinated centrally, subject to agreements with grid operators.

Modelling Approaches
Automating the generation of models for Step 3: Scenario Modeling, which calculates the maximum and minimum flexibility required for the specific event, is required to implement the methodology in an online real-time way. The approach proposed herein is to implement models for specific systems in conjunction with an optimization formulation. A range of potential modeling approaches were identified for a number of flexible systems. Models were evaluated to enable the selection of the most suitable modeling approach for each system. The model types and flexible systems included (1) data-driven predictive model applied to air handling unit (AHU) fan control; (2) a grey box model applied to thermal building systems, specifically heat pumps; and (3) a physics-based model with parameter estimation applied to Photovoltaic (PV) renewable generation; and (4) an electrical model applied to the battery system.
In an online real-time scenario, when a signal is received from an aggregator requesting to be informed of the available flexibility, the models are activated or are continuously predicting, and the outputs are input to the optimization algorithm. The model outputs consist of the available power flexibility for the time period of the demand response request for each system. The optimal combination of systems is then selected by the optimization algorithm.
The models have been chosen to minimize data gathering requirements both in (a) number of variables and in (b) length of time series required for prediction. For example, the internal power electronics of the battery system are not required to be modeled as this would have very little impact on the available flexibility provided by the battery and would make for a computationally expensive model. Losses internally in the system may be modeled through the use of the charging efficiency and modeling of the discharging efficiency.
Interactions between systems may be accounted for either in the models or at the optimization stage. For the combination of systems under consideration in this work, there may be an interaction between the thermal systems such as heat pumps and ventilation fans. Changing the AHU fan speed may cause interactions with the variable refrigerant flow (VRF) heat pump system. This may be accommodated by using the predicted fan speed as an input to the grey box model for the thermal building system. There would then be two model outputs for the thermal system, one with standard fan speed and one with modified fan speed.
Occupant comfort requirements may be incorporated in the grey box thermal model. Indoor air temperatures are required to stay within adaptive comfort boundaries as defined in ASHRAE 55 [28] or other adaptive comfort standards such as EN 15251 [29], or maximum CO 2 levels may be specified in accordance with CIBSE guidelines [30].

AHU Fan Model
A data-driven model may be suited to modeling the electrical power of an AHU fan to predict its flexibility. The relationship between electrical power and ventilation rates is less complex than thermal systems, and provided any change in fan speed is linked to the thermal system model to ensure interactions are captured the model may rely on data. Seasonal effects are not as impactful, as ventilation requirements may be independent of outside air temperature, and a number of approaches are possible depending on the control variable for the fans.
Of the data-driven modeling approaches available, autoregressive-exogenous (ARX) [31] and autoregressive moving average (ARMAX) [32] models are the most commonly used for simulation of building energy consumption. Reviews of data-driven modeling approaches in buildings include numerical models for the prediction of electrical loads in buildings [33], regression models [34], energy prediction in buildings [35] as well as large-scale datadriven modeling of energy in buildings [36]. Zeng et al. [33] found that multivariate linear regression and non-linear support vector machines produced better results, particularly when buildings had unpredictable and multifaceted patterns of energy consumption.
Focusing specifically on the requirements for a data-driven model of an AHU fan, previous studies used CO 2 data to predict occupancy [37] or occupancy and ventilation rate to predict CO 2 concentration [38]. However, a model for flexibility would need to be different as the objective is to predict fan power from CO 2 data.
For the development of the model, if the fan is controlled at a fixed speed setpoint only, the setpoint may be reduced subject to constraints, and the model may be simple. If the fan speed is controlled by CO 2 , a data model using CO 2 measurements as an input may be appropriate, subject to the constraint of a CO 2 limit, e.g., 1000 ppm [30]. The data-driven model may also be correlated against occupancy so that the predicted fan speed "learns" occupancy if occupancy data is available. The fan speed setpoint s, then determines the To determine the available fan flexibility, f F in (kW), the following formula may be used: (1) where P C is power at current fan speed and P p is predicted fan power, a function of CO 2 concentration. The recommended control variable is the fan speed. The model objective is to maximize the power reduction by the fans, subject to CO 2 constraints. Training data may be taken from representative days rather than all days [39], for example, weekdays with typical occupancy.

Thermal Systems Model
Modeling the electrical response of thermal systems in the building is a complex process that typically requires detailed thermal modeling of the building physics, known as white-box models [40]. As was seen in the previous section, data-driven models, sometimes referred to as numerical or black-box models, have also been used [31][32][33][34][35][36] but (i) parameter identification may not map to real parameters in the building, e.g., thermal conductivity, and (ii) prediction is reliant on the quality and diversity of the data available. In recent years, a number of grey box modeling approaches have been developed which combine knowledge of building physics with some data-driven aspects. The most common type is an R-C model, an analogy of the resistance-capacitance approach of electrical circuits. Grey box model libraries for building energy have been developed in the modelica modeling language by a number of research institutions [41][42][43][44]. An alternative R-C approach is the state-space model developed by Bacher and Madsen [45]. Other grey box approaches start with complex whole building simulations such as EnergyPlus, and reduce model complexity incrementally until the model is computationally efficient [46] and meets operating requirements. Computation times required for prediction horizons for some of these approaches are long, 80 h for full model reduced to 3 h for simplified model [25], which may be too long for intra-day flexibility events of 1 h or less.
The recommended control variable for thermal systems modeling is the global indoor air temperature set point. The model objective is to determine the electrical power reduction or increase that can be achieved for a specific change in temperature set point. Input parameters include outside air temperature and may include weather-related data, e.g., solar radiation, wind, relative humidity, and occupancy data.
The state-space model developed by Bacher and Madsen [45] and implemented by Roels et al. [31] uses a grey box approach based on the R-C principle. It starts with fitting a simple model and adding additional terms to simulate additional physical parameters, e.g., solar radiation or wind, until the loglikelihood plateaus and residuals are equivalent to white noise. The loglikelihood or maximum likelihood estimation is a statistical function which represents the combination of parameters which maximize the goodness of fit of a statistical model. Two potential implementations of this model are proposed. The first adapts the model for electrical power. The second implementation is specifically for a heat pump application and uses the temperature output of the state space model as an input to the heat transfer equation, applying the heat output to Coefficient of Performance (COP) to calculate power reduction or increase. To implement the state space approach for power flexibility, the model would need to be adapted to focus on electrical power, P, instead of temperature, T as shown in Equation (3): whereby the state vector is P, the input vector is U, T the Temperature vector, and parameter vectors are denoted A, B, and C. Φ is the energy flux from solar radiation s, and heating system h. The Weiner process, assumed to be a Gaussian white noise process analogous to an error term, is denoted by ω.
The second approach is to use the un-adapted state-space model for determining temperature changes only, then input these into the heat transfer equation to calculate the heat reduction or increase: where the mass flow rate isṁ of the heat transfer fluid, typically refrigerant, water, or a water-glycol mix, and C p is the specific heat capacity at a constant pressure of the heat transfer fluid. COP is then used to calculate electrical power, P: Thermal system flexibility, f T , equivalent to x, the proportion of flexible power, by P, the electrical power consumption of the heat pump, is calculated as: As COP varies with condenser and evaporator temperatures, heat pump or compressor power curves in lookup table format are required to define the relationship between COP and input-output temperatures.

Photovoltaic (PV) Model
A physics-based PV model such as that developed by Zhou et al. [47] may be used to predict the output of the PV installation. The key inputs are solar irradiance and peak installed capacity. The model incorporates the PV curve with a thermal drop-off effect and a number of characteristic parameters of PV modules. Parameter estimation was included. Forecasts for solar irradiance may be obtained from subscription services [48]. Another option for modeling PV output may be to deploy learning algorithms based on weather forecasts [49]. For both model types, the key to accurate prediction is reliable solar irradiance forecasts. Accuracy of prediction of PV output is dependent on the accuracy of the solar irradiance forecasting at the exact position of the PV installation. Ito et al. [50] adopted a geographical area approach for predicting the output from large-scale gridconnected PV, for example, one prediction point per region. A building scale PV model which incorporates partial shading [51] may account for extremely localized effects such as a cloud passing over the building and may be more beneficial in this context.
A control variable is not required for PV as the output of the panels is not controllable. The objective of the model is to predict the power output of the PV installation for a given predicted solar irradiance. The model proposed by Zhou et al. [47] determines the power output for an array, P A , based on the number of parallel-connected modules N p and the number of series-connected modules N s .
The PV module maximum calculated power output, P M , given solar irradiance G, is defined as The series resistance is R S (Ω), V oc (V) is the open-circuit voltage under normal solar irradiance G (W/m 2 ), while V oc0 is the voltage under the standard solar irradiance G 0 . The normalized quantity of the open-circuit voltage with reference to the thermal voltage is v oc . The short circuit current is I Sc (A) while I sc0 is the short circuit current for G 0 , the standard solar irradiance. The module temperature is T (K) and T 0 is the temperature of the PV The angle of installation at which the PV panel is installed must also be considered and the difference between solar irradiance measured on the horizontal by a pyranometer located near the PV array calculated. The angle of declination will change during the day as the sun rises and falls, and over the year as the sun declines or rises in the sky. This may either be calculated or taken from a look-up table for a particular geographical location [52]. Parameter estimation for β, γ, α, n and R S at the maximum power point (MPP) was conducted using experimental data, and the values determined [47]. α = 1.21, β = 0.058, γ = 1.15, nMPP = 1.17 and R S = 0.012 Ω. Other parameters may be identified from the PV manufacturer's datasheet [53], examples of which include the following: P max = 250 W p when G = G 0 , V oc0 = 37.8 V, I sc0 = 8.28 A, G 0 = 1000 W/m 2 , T 0 = 25 • C.

Battery Model
An electrical model was developed by the author for the battery system, starting from the mathematical model for a lithium-ion battery system developed by Berrueta et al. [54]. The detailed model was reduced to the essential elements as, for flexibility, the main priorities are high level, such as accounting for the round-trip efficiency losses in charging and discharging the battery system. Modeling internal aspects such as the kinetics of the chemical reaction or ion transport may be accounted for in the overall system efficiency. Power limits such as a minimum state of charge (SOC), required to be retained in the battery system for safety reasons, may be included as constraints in the model. Variations in power limits and overall efficiencies do occur and were determined to be significant by Sakti et al. [55] as constant efficiencies overestimated flexibility by 10% for electrical grid applications. A piecewise linearized model based on the R-C principle, developed by Gonzalez-Castellanos et al. [56], provides a means of determining varying real-time efficiencies and power limits for Li-Ion batteries. This may be coupled with the approach outlined below to calculate real-time efficiencies for both charging and discharging. Battery capacity also degrades over time, in particular, second-life electric vehicle batteries [57], and it is recommended that capacity be re-assessed during medium to long-term operation.
The model objective is to determine the maximum power increase or reduction the battery system can provide for a specified time period, assuming a constant rate of discharge over the time period. The recommended control variable is the charging or discharging rate in kW or kVA. Two types of efficiency are to be considered, η c the coulombic efficiency during discharging and η e the energy efficiency during charging. The voltage charging the battery system is greater than the voltage discharging the battery system due to non-ideal processes, therefore, the efficiency losses are different. The flexibility of the battery system, f B may be modeled as shown in Equations (10)- (16).
General form: Charging: Discharging: The energy capacity of the Li-ion battery is C (kWh), SOC is the State of Charge (%), and ∆t j (h) is the duration of demand response event, j. To take into account active states of the battery, i.e., if the battery is already in the process of charging or discharging and will continue to do so between the time of flexibility measurement, t 0 , and the commencement time of the demand response event, t i the terms identified in Equations (13) to (16) may be added to the model. P ch is the charging power (kW) of the battery system, P d is the discharging power (kW).
Flexibility event requests a discharge while the battery is charging (negative flexibility): Flexibility event requests a charge while the battery is charging (positive flexibility): Flexibility event requests a discharge while the battery is discharging (negative flexibility): Flexibility event requests a charge while the battery is discharging (positive flexibility): where operational limits for the battery system, such as maximum and minimum charging (P c ) and discharging (P d ) power P d,minimum are less than P d , which is less than P d,maximum and P c,minimum , which is less than P c , which is less than P c,maximum ; minimum SOC and capacity degradation over time, may be modelled as constraints.

Optimization
Optimization using mixed-integer linear programming (MILP) has been utilized by a number of demand response applications. Scenario-specific objective functions [58], a stochastic approach [59], and applying models as constraints for the objective function for storage [31] are among the MILP formulations previously developed.
A data-driven flexibility function for characterization and control of energy flexibility was developed under the International Energy Agency (IEA) Energy in Buildings & Communities Annex 67 on Energy Flexible Buildings [60]. This approach incorporates data-driven models of flexible systems in buildings and is dynamic in that the prediction of flexibility varies based on the data from the systems. As it is a data-driven approach, it may not be compatible with some of the grey box or physics-based model types proposed in this paper.
A formula for flexibility developed by the author [16] shown in Equation (17) may also be utilized as an objective function to maximize power flexibility for a given time period. The primary requirement for the objective function is that it be independent of use cases. Use cases may be addressed by varying constraints tailored to that specific use case. max n ∑ s=1 f s t j (17) Constraints may include: flexibility may be positive or negative but not both at the same time; flexibility and flexibilities of systems are real numbers; f s ∈ R; SOC of the battery may never be zero but shall be above a set limit of such as 5% minimum charge, e.g., SOC min < SOC; indoor air temperature is required to be within adaptive comfort ranges; CO 2 levels must be below 1000 ppm and other model related constraints.
MILP was used by Siebert et al. [58] to create objective functions relating to specific use cases: maximizing financial return, maximizing bid length, and maximizing power peak.
Max ( ∑ t ∑ r P(t, r).W(t). ut) (18)  whereby P(t, r) (kW or MW) is the power flexibility provided by the flexible source r for the timestep t (s min or h) and is calculated by: P(t, r) = Res(t, r).E(t, r) (19) Maximizing bid length: Subject to the constraint: ∑ r Res(t, r).E(t, r) ≥ Bid(t).P min (21) Maximizing power peak: Max (Max ∑ t Res(t, r).E(t, r)), t τ A Boolean indicator Res(t, r) determines if a flexible source r is available to provide flexibility for the duration of timestep t. The quantity of power variation is denoted E(t, r) (kW or MW) as delivered by flexible source r during timestep t. An energy price weight, W(t), is applied to time step t, relating to electricity price on a scale from 0 to 2.5. The length of the optimization is denoted by ut (s). Similar to Equation (19), Bid(t) is a Boolean indicator for bid activation at time t. The minimum power that is required to be provided is P min (kW). The set of possible time intervals from the start of the demand response bid to the end of the event is τ.
A more generalized flexibility optimization formulation, rather than discrete cost functions for different scenarios, would be more advantageous, as would be independent of pricing and market structures in different countries. As an alternative to a weighting factor for price, incorporating actual power or energy prices would be preferred, as the costbenefit to the building or site may then be generated as an output for informed stakeholder decision-making. Binary 1/0 Booleans may minimize or exclude the contribution of sources with partial flexibility as the binary approach may limit the application to on/off sources. For fast demand-side services such as ancillary services, maximizing financial return may be decoupled from bid length and power peak. However, for longer events that may be priced on an energy basis, there may be synergy between maximizing financial return and bid length or power peak. Power peak maximization by increasing consumption may be required in situations with excess renewable generation such that grid balancing necessitates increases in system load, known as positive or forced flexibility [15].
A stochastic approach, also using MILP, was developed by Ottesen and Tomasgard [59]: The probability of scenarios is given by R s . The grid import power x import (kW) excludes flexibility (kW). Energy prices are accounted for by P energy (€). The power consumption in peak periods x peak (kW), is multiplied by the associated cost penalty P peak (€), if present. Should the building be supplied from different energy carriers such as electricity, gas, and district heating, the cost of ramping or starting energy conversion is included through G startup (€) coupled with a binary variable, α, to activate and deactivate the energy supply. Curtailment of load and other disutility losses are accounted for by cost X d (€) multiplied by load reduction Φ (kW). If export of renewable generation x export (kW) from the building is permitted, and a financial return P sales (€) is provided by the grid utility, this is also included.
Flexibility is maximized implicitly in Equation (23) through minimizing grid import power. While the concept allows for different energy vectors or energy carriers, such as on-site renewable generation or district heating, to be incorporated, the cost function is focused on electricity optimization. The proposed approach is more generalized than Siebert et al. [58], which is advantageous as it makes it more widely applicable without modifications being required. The optimization assumes certain fixed pricing structures are present in the demand response market, e.g., peak penalties such as critical peak pricing (CPP) and that tariffs are based on kW or MW consumption. However, it is not clear if the formulation is adaptable enough to incorporate dynamic, real-time, or other pricing signals which may be used in demand-side programs, for example, positive flexibility or capacity markets. Ramping or converter startup costs may be significant for on-site generation, such as combined heat and power (CHP) and frequent cycling may be detrimental to the equipment but are minimal for battery storage systems. On the other hand, efficiency losses due to charging and discharging battery systems are not accounted for. Grid export regimes are specified by grid utilities and vary widely depending on the jurisdictions and size of on-site renewable generation. If export is not permitted, x export in Equation (23) will be zero, or if it is permitted but not rewarded, P sales will be null, thereby removing the final term.
An MILP objective function combined with simple models of building systems and an ARX model of building thermal storage used a cost-based approach [31]. Heat pumps may be used for heating or cooling, and in this instance, were used for cooling. The optimization shown in Equation (24) was implemented in Gurobi. The objective was to minimize electrical grid import costs, PV costs, and cooling system maintenance costs. P denotes power (kW), Ep is the grid electricity price (€/kWh), ε is the expense of maintaining the PV and cooling systems, denoted by the superscript CS, respectively (€/kW installed), i specifies the number of chillers, while the duration of individual time steps, ∆τ, is one hour (h). min ∑ τ∈T P grid . Ep .∆τ + ε PV .P PV + ε CS . ∑ i∈I Q i .∆τ (24) The Gurobi optimization software is a commercial solver for linear programming, mixed-integer programming, and quadratic programming [61]. Other similar solvers include Cplex [62] and SCIP [63]. The Gurobi MILP algorithm has been used for a wide range of optimization problems such as process planning in manufacturing [64] and scheduling of container crane loading operations for shipping [65].
While cost may be a factor in many demand response events, it is not universally the case. CO 2 minimization, for example, does not have a monetary cost factor. The maintenance cost of the PV system is independent of the time step and PV is not controllable, so it is not clear why it has been included. It is not clear how maintenance costs for the chillers in the cooling system are impacted by participation in demand response events. For split system heat pumps, the most commonly used type in office buildings, compressor cycling on part load is a routine part of their normal operation. Some types of large industrial chillers, e.g., ammonia chillers, do require on/off operation at full load, and cycling may have some slight impact on their maintenance cycles. However, these chillers are installed for large loads. e.g., chill warehousing, or are coupled with purpose-built storage, which dampens cyclic effects.
Elements of each of the above approaches have merit and may be adapted and combined in the development of a more universal objective function for flexibility which is applicable to a wide range of use cases, permits linking to a diverse range of model types and includes pertinent variables.

Conclusions
A framework for automating a simultaneous aggregator-based assessment and activation of flexibility in multiple buildings has been developed. Implementation of the framework requires automation of "Step 3: Scenario Modeling" from the flexibility assessment previously developed by the author, which outputs the available flexibility defined for the specific time period of the demand response request for each building or site. Appropriate model types were selected for each flexible system, a data-driven model for the AHU fan, a grey box model for the heat pump system, physics-based for the PV panel, and mathematical for the battery system, incorporating data-driven elements for adaptive efficiency. The models were then adapted for flexibility. The flexibility output from the models is then input to a proposed optimization formulation. Based on a review of optimization applied to the flexibility domain, MILP was found to be the most appropriate approach. Future work includes implementation of the framework in a simulation environment, addressing the technology gaps relating to interoperability of building and site automation systems, and practical installation of the framework at pilot sites to demonstrate use cases and provide experimental results.