Integrated cost-optimal residential envelope retrofit decision-making and power systems optimisation using Ensemble models

Energy & Buildings 214 (2020) 109833

Contents lists available at ScienceDirect

Energy & Buildings journal homepage: www.elsevier.com/locate/enbuild

Integrated cost-optimal residential envelope retrofit decision-making and power systems optimisation using Ensemble models Carlos Andrade-Cabrera a,∗, Ciara O’Dwyer b, Donal P. Finn a a b

School of Mechanical and Materials Engineering, University College Dublin, Belfield, Dublin 4, Ireland School of Electrical and Electronic Engineering, University College Dublin, Belfield, Dublin 4, Ireland

a r t i c l e

i n f o

Article history: Received 17 July 2019 Revised 30 January 2020 Accepted 2 February 2020 Available online 8 February 2020 Keywords: Ensemble Calibration Retrofit decision-making Building-to-grid Power systems optimisation

a b s t r a c t The dynamic integration of building energy models and power system models is essential to analyse the potential supply-side benefits of adopting Energy Conservation Measures (ECM) at the national, regional scale. The integration of these models requires the development of linear building energy models which are numerically compatible with linear power systems models. Ensemble Calibration is an automated model calibration methodology which identifies a cluster of lumped parameter building energy models (denoted Ensemble models) using an archetype building energy model. Each lumped parameter building energy model represents an ECM configuration (i.e., a combination of ECMs). The current paper introduces a novel mechanism by which Ensemble models are integrated with power systems models in a manner such that cost-optimal retrofit decision-making problems and power systems optimisation problems can be simultaneously solved. To achieve this objective, Ensemble models are reformulated as bi-linear models, which are then co-optimised with power systems models using a multi-stage optimisation algorithm. The methodology is tested using two building energy archetypes: a detached house archetype and a mid-floor apartment archetype. The archetypes are representative of the Irish residential stock built prior to 1970. The results show that the proposed multi-stage optimisation algorithm provides computational advantages to the solution of the equivalent problem using a brute force approach (i.e., solving building-to-grid models for each ECM configuration). The proposed method is 4.73 times faster than the brute force approach when the archetypes are decoupled and 73 times faster when the models are deemed to be coupled via a power systems model. © 2020 Elsevier B.V. All rights reserved.

1. Introduction 1.1. Integrating residential retrofits and energy systems Residential dwellings represent an important opportunity for societal decarbonisation in northern European countries. Approximately 80% of residential end-use energy consumption in northern European countries corresponds to space and domestic hot water (DHW) heating [10]. It is estimated that 87% of the existing building stock in the UK will be still standing by 2050 [49]. It has been shown by ECF [23] that it is possible to achieve current European residential decarbonisation targets (i.e., abatement of 90% of residential emissions with respect to 1990 levels by 2050) by means of using existing technologies for building retrofits. Retrofit measures, if applied properly, significantly improve indoor thermal comfort and provide monetary savings to end-users [16]. From a ∗

Corresponding author. E-mail address: [email protected] (C. Andrade-Cabrera).

https://doi.org/10.1016/j.enbuild.2020.109833 0378-7788/© 2020 Elsevier B.V. All rights reserved.

societal perspective, the introduction of retrofit measures at large scale translates to an aggregated reduction in energy related emissions. For example: if residential retrofits are combined with electrified heating systems such as direct resistive heaters (DHR), storage heaters (SH) and heat pumps (HP), the carbon emissions associated with space and DHW heating can be displaced to the electricity sector, with the potential of significant carbon abatement as electricity could be supplied via low-carbon power generation. The aggregated effect of retrofit measures has the potential to significantly influence the current operation and the future design of the energy systems [34]. Without retrofits, a large scale integration of electrified heating systems (e.g., HP) implies a significant increase in peak load and thus in transmission, distribution and generation capacity [42], [45]. One attractive opportunity of the value of flexibility of electrified residential heating consists of the cost-optimal exploitation of building thermal mass at large scale. Storage heaters, which use electrical resistance elements to heat high-density bricks in an insulated casing [21], have the potential to minimize electrified heating costs in future generation scenar-

2

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

Nomenclature Variables and parameters α Calibration parameter (Ensemble Calibration) [-] β Calibration parameter (Ensemble Calibration) [-] C Lumped Capacitance [J/kgK] δ ECM decision variable x Variation in layer thickness increment [mm] d Disturbance vector F System matrix (state) Set of fixed parameters (Ensemble Calibration) G System matrix (heating input) g Solar transmittance (Window) [-] H System matrix (disturbances) k Simulation time-step M Index for external insulation (Ensemble) n Number (combinations) N Index for internal insulation (Ensemble) NH Investment horizon [hours] O Index for ceiling insulation (Ensemble) P Index for window replacement (Ensemble) p Set of calibration parameters Q Heating load [W] R Lumped resistance [m2 K/W] T Temperature [ ◦◦ C"?>] U Thermal transmittance (U-value) [W/m2 K] x State vector u Electrified heating input V Set of variable parameters (Ensemble Calibration) Z Output matrix Subscripts amb Outdoor temperature ceil Ceiling node comb Number of combinations W West wall (node) ext External wall node gnd Ground node heat Heat input int Internal mass node inf Infiltration k Time-step index lin Linearisation (trajectory) r Room temperature s Solar gains wall External walls win Window w Wall node N North wall (node) Superscripts  Optimal (calibrated) parameter − Variation in dynamics (bi-linear model) Acronyms ACH Air Changes per Hour AIPS All Island Power System BAU Business as Usual BEMS Building Energy Model Simulation ECM Energy Conservation Measure EEA European Energy Agency EPS Expanded Polystyrene GCH Greenhouse gas (emissions) IPCC Intergovernmental Panel on Climate Change IWEC International Weather for Energy Calculations

MILP NPV RES SEAI

Mixed Integer Linear Program Net Present Value Renewable Energy Sources Sustainable Energy Agency of Ireland

ios [51]. Furthermore, thermal mass can be seen, from the energy management perspective, as a resource that allows for heat storage and/or load shift as discussed in Heinen et al. [32]. For example: an excess of available wind power can be dispatched to the thermal mass (e.g., pre-heating) as well as other residential electricity storage devices (e.g., batteries and/or storage heaters), thus reducing wind curtailment and potentially reducing electricity generation costs and carbon emissions. Retrofit decision-making, heating systems technology selection, energy systems operations and planned energy systems investments are all interdependent. The integrated assessment of building retrofit measures (including electrified residential heating technologies) and cost-optimal power generation requires the development of an integrated building-to-grid retrofit modelling framework (e.g., Patteeuw et al. [48]) capable of simultaneously considering cost-optimal retrofit decision-making and power system optimisation problems. An integrated modelling framework has the capability to provide researchers with insights on the potential environmental and economic benefits of building-to-grid integration (e.g., heating systems technology selection, capacity expansion, wind curtailment reduction) in a post-retrofit scenario. 1.2. Retrofit decision-making in building-to-grid modelling Steady-state normative residential demand models such as the CIBSE model [19] are the most numerically efficient method to represent building energy performance for building-to-grid modelling integration purposes (e.g., Jennings et al. [35]). However, normative models do not fully capture transient building thermal performance and are known to overestimate the required energy consumption [24]. On the other hand, building energy archetypes are increasingly used in the study of aggregated residential demand at the urban and national level [52]. These archetypes are often developed using Building Energy Model Simulation (BEMS) tools (e.g., Jermyn and Richman (2016), Cerezo et al. (2017): Muringathuparambil et al. (2017)). It is possible to iteratively integrate BEMS archetypes and power systems models via co-simulation environments such as BCVTB [58]. However, such an approach would remain inherently sub-optimal because the power systems model would not be able to predict residential demand requirements. Therefore, building thermal mass is not exploited at large-scale, thus resulting in a sub-optimal exploitation of conventional and variable renewable electricity generation assets. For an integrated energy scenario analysis, where building archetype models and power systems models need to be combined, simplified linear building energy models are required. Dynamically integrated linear building-to-grid models (e.g., Good et al. [26]) have the potential to facilitate building retrofit analysis that accounts for the effect of increased penetration of electrified space and water heating systems, in the context of the wider integration of renewable energy generation into the electricity grid. Lumped parameter building models (also known as ‘RC thermal networks’) are linear models compatible with linear buildingto-grid modelling [27]. Under this modelling framework the electric resistances of the RC thermal network represent the thermal resistances of the walls and the electric capacitances represents the thermal mass of the walls. One key feature of the lumped parameter building modelling framework is that model parameters have a physical interpretation (e.g., external wall resistance

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

is analogous to thermal wall transmittance). Theoretical white-box RC models (i.e., models where all parameters are known) can be readily identified can be created using information readily available in BEMS archetypes (e.g., [26]). Ensemble calibration relies on grey-box modelling [15], whereby the model parameters differ from the theoretical value but retain their physical meaning. In doing so, grey-box models can better capture the transient dynamics of the BEMS archetype model, which allows for a more accurate estimation of residential load flexibility at large-scale. Prior research introduced in Andrade-Cabrera et al. [5] and AndradeCabrera, Turner, and Finn [6] exploited this semi-physical interpretation to introduce the Ensemble Calibration framework. This framework enables researchers to obtain a cluster of optimally calibrated lumped parameter building models for a desired number of ECM configurations of a BEMS archetype. The Ensemble Calibration framework results in a continuous-time domain structure whereby all models of retrofitted buildings are related to the parameters of a baseline (i.e., pre-retrofit building model). In the discrete-time domain, the Ensemble Calibration methodology results in a cluster of optimally calibrated linear lumped parameter building energy models related to the discrete-time baseline model. The cluster of calibrated models identified via Ensemble Calibration is denoted an Ensemble model. The aim of the current paper is to introduce a discrete-time domain optimisation methodology which exploits the intrinsic property of the Ensemble model (i.e., numerical relationship between building models) for a numerically efficient implementation of the cost-optimal retrofit decision-making problem integrated with power systems optimisation problems. 1.3. Contributions The current paper introduces an optimisation methodology for the seamless integration of cost-optimal retrofit decision-making and power system optimisation using Ensemble models. While prior research [5], [6] focused on the identification of building energy models that are representative of retrofits (i.e., identifying Ensemble models), the current paper focuses on an integration mechanism of the Ensemble models (i.e., how to use Ensemble models) and power systems models in a manner such that cost-optimal solutions are identified. The current paper demonstrates that Ensemble models can be algebraically manipulated in a way such that the cost-optimal retrofit decision-making problem is expressed as a bi-linear optimisation problem. Bi-linear problems are potentially intractable (i.e., a globally optimal solution cannot be found in a reasonable amount of time, Adams and Sherali [1]). Hence, a multi-stage optimisation algorithm is introduced to remove this bilinearity by means of pre-computing a linearisation trajectory. The use of a linearization trajectory transforms the bi-linear problem into a Mixed Integer Linear Problem (MILP), which is a tractable problem readily solvable with commercial solvers. It will be shown that due to building physics constraints, the multi-stage optimisation algorithm needs to solve the bi-linear optimisation problem twice (first for models featuring external wall insulation and secondly for models not featuring external wall insulation). The bilinear optimisation problem that gives the least cost is a solution of the integrated cost-optimal ECM configuration subject to the energy systems optimisation problem. The proposed formulation avoids the need to solve an individual building and grid co-optimisation problem for each possible ECM configuration (i.e., a brute-force approach, as in Hilliard et al. [33]). The current paper compares the proposed methodology with a brute-force approach in the context of retrofit decisionmaking with electrified space heating technologies. It is shown that the proposed methodology provides significant computational savings while providing an optimal solution to the cost-optimal ECM decision-making problem. The focus of the current paper is

3

on the computational advantages of the proposed methodology prior to model integration. For clarity, and due to space limitations, a price signal is used as a surrogate unit to demonstrate the potential of building-to-grid integration. The integration of the proposed methodology with a detailed power systems model representative of the All Island Power System (AIPS) will be presented in a subsequent contribution. 1.4. Structure of the current paper The current paper is organized as follows: Section 2 provides background information on previous approaches towards the integration of ECM decision making tools and power systems using BEMS archetype models. Section 3 describes the methodology introduced to integrate Ensemble models with power systems models using a multi-stage optimisation algorithm. Section 4 describes the results of the proposed modelling framework as applied to the study of ECM decision-making of two archetype building energy models representative of the existing detached house stock and mid-floor apartment stock in Ireland using a night-time electricity tariff. The paper closes with the discussion of the results and introduces future lines of work. 2. Background 2.1. Retrofit investment planning using BEMS tools Current methods in techno-economic building retrofit optimisation, in the context of BEMS, often relies on the integration of heuristic optimisation techniques (e.g., Genetic Algorithms) and BEMS tools (e.g., Ferdyn-Grygierek and Grygierek [25], Penna et al. (2015), Asadi et al. [9]). The BEMS models are used to calculate the energy consumption and/or economic costs of the target retrofit measures and the heuristic optimisation algorithm uses such information to identify the cost-optimal retrofit configuration. However, power systems models are often modelled as Mixed Integer Linear Programming (MILP) problems [14], [50]. Therefore, BEMS-driven ECM decision-making (i.e., heuristic optimisation) and power systems models (MILP optimisation) are numerically incompatible. An alternative to avoid the numerical incompatibility between BEMS models and power systems models consists of using BEMS tools to generate synthetic building performance data as an input to power systems optimisation tools. Ault et al. [12] used building models representative of a Scottish housing estate in order to calculate the aggregated heat demand forecast for the entire estate at 15 min intervals. The models were developed using the ESPr building simulation environment [22]. The electrified heat load forecast was then used as an input to a power system model which produced an optimal charging schedule for the estate. The method can be extended in the context of ECM decision-making by means of pre-calculating post-retrofit heating load time series. The time series are then used in a power systems optimisation framework. Wu et al. [59] simulated energy demands from archetypes for the Swiss village of Zernez. The archetypes were developed using the EnergyPlus building simulation environment. Electricity demand profiles extracted from the archetype are then optimised in an energy hub optimisation context (i.e., a multi-energy systems environment). A key disadvantage of a sequential approach is the omission of building flexibility (e.g., storage in building fabric) and electricity grid flexibility requirements (e.g., wind energy integration). BEMS are unable to adapt to dynamic events occurring in the power systems model (e.g., availability of variable generation or demand response events) unless a potentially sub-optimal iterative and computationally-expensive strategy is considered. For example, if the models are integrated, an excess of available wind power can

4

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

Fig. 1. BEMS residential archetype models [47]. Left: detached house archetype. Right: mid-floor apartment archetype.

be dispatched to electricity storage devices (e.g., batteries or storage heaters) to reduce future generation costs, thus reducing wind curtailment and potentially reducing electricity generation costs and carbon emissions. Therefore, if BEMS archetypes are to be considered, an integrated building-to-grid model requires the explicit definition of a linear building energy model (e.g., lumped parameter models) in a manner such that it is numerically compatible with a power systems model for building-to-grid co-optimisation (e.g., Anwar et al. [7]). A further limitation arises due to the potential dimensionality of the problem. A brute force approach of the cost-optimal ECM decision making problem using building-togrid models (i.e., solving building-to-grid models for every desired combination of ECMs, as in Hilliard et al. [33]) can become computationally prohibitive. Such an approach implies re-running potentially complex power systems models for a large number of ECM configurations (i.e., a combination of ECMs). While technically feasible, a brute force approach may not be computable in a reasonable amount of time (i.e., it is an intractable approach). 2.2. Ensemble calibration Ensemble Calibration [5] is an automated calibration methodology capable of transforming any residential BEMS archetype models into a lumped parameter archetype building model, representative of an ECM configuration. In the context of the current paper, an ECM configuration is defined as a combination of ECM measures (e.g., 100 mm of external wall insulation combined with 200 mm of ceiling insulation). The current study considers two building energy archetype models, which are representative of detached houses and mid-floor apartments built in Ireland before 1978. The archetypes (Fig. 1) are described in Neu et al. [47] and AECOM [2]. The archetypes feature double glazing windows (uwin = 2.88 W/m2 K, gwin = 0.759), heavy solid wall construction (Uvalue = 2.27 W/m2 K), an insulated ceiling (Uvalue = 2.37 W/m2 K) and an infiltration rate (0.67 ACH per zone). Two active occupancy (i.e., comfort) periods are observed in accordance with the Dwelling Energy Assessment Procedure (DEAP) v3.2.1 [54]: 7 AM to 9 AM and 5 PM to 11 PM. Weather is assumed to correspond to the IWEC2 Dublin weather file [11]. Fig. 2 shows a 9R6C heterogeneous lumped parameter building model topology associated with a detached house archetype introduced in Neu et al. [47]. The reader is directed to AndradeCabrera, Turner and Finn [6] for a detailed explanation of the lumped parameter network topology associated with the mid-floor archetype. The dwelling is heated with direct electric resistance heaters, which are modelled as perfectly efficient convective inputs and modelled using the heating power input Qheat . Latent loads are not considered in this model as it is desirable to obtain a linear model. Furthermore, convective loads have been shown to be sufficient in building-to-grid modelling (e.g., [26]). Nodes Cw1 and

Cw2 model the outer and inner leaves of the external walls, which are analogous to the thermal mass of the building [27]. The solar gains due to solar radiation on the walls, Qs, wall are applied directly to node Cw1 . Cr represents the capacitance of the room air mass with room temperature Tr . Node CIM captures the thermal mass of the internal partitions and other slow dynamics. The window solar heat gain Qs, win and the heating power input Qheat are split between Cr and CIM via parameters f1 and f2 . Node Cc models the ceiling between the room node Cr and the attic temperature Tatt . Finally, the ground node Cgnd is added to model the heat transfer between the conditioned air volume and the foundations. The temperature Tgnd models the temperature under the conditioned air volume (boundary condition ‘Ground’ in EnergyPlus). The building does not feature an attic zone, as the attic space is ventilated and considered to be in thermal equilibrium with the ambient outside during the heating season. Therefore, the attic temperature Tatt is extracted from EnergyPlus. This simplification enables the analysis of the proposed methodology without the additional complexity of modelling the dwelling as a two-zone model. Only the living room (18.5 m2 out of 160 m2, or 11.5% of the heated floor area) is considered as a living space in the original archetype. This implies that only that room has a thermostatic set-point of 21 °C as per DEAP v3.12 guidelines [54] and all other rooms (e.g., bedrooms) have a thermostatic set-point of 18 °C. Hence the effective set-point of the combined single zone, after area weighting, is 18.34 °C. Thermal bridges are not modelled in the current paper as they cannot be explicitly represented in the EnergyPlus building simulation environment. The accurate modelling of the thermal effects of thermal bridges would require an increased level of modelling effort both for the archetype model (e.g., 2D thermal simulation) and a modified lumped parameter model. Domestic hot water consumption and internal gains are not considered at this stage of modelling as it requires the modelling of occupancy profiles, which goes beyond the scope of the current paper (e.g., Buttitta et al. [17]). Thermal comfort is guaranteed as part of the cost-optimal building-to-grid co-optimisation problems introduced in Section 3.3. In a brute force approach (or enumeration approach) of the building-to-grid problem (e.g., Hilliard et al. [33]), it would be necessary to calibrate a lumped parameter building energy model for every possible ECM configuration. Then, it would be necessary to run a building-to-grid model for each calibrated model. Such an approach could easily become computationally intractable if the power systems model is complex or if multiple archetype models are to be considered. In Ensemble Calibration, however, all building models (i.e., pre-retrofit and retrofitted scenarios) are calibrated simultaneously with the aim to obtain a mathematical model whereby the retrofit decision-making problem and the power systems optimisation problem are simultaneously solved. The building model parameters are split in two sets: one with the parameters deemed to be variable for a given ECM (denoted V),

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

5

Fig. 2. Lumped parameter model topology (detached house archetype).

and a set of parameters deemed to be fixed (denoted F). It can be assumed, via a semi-physical interpretation of the model parameters, that the inclusion of ECM measures to the envelope (i.e., adding external wall insulation), will translate in a variation of selected lumped building energy model parameters (i.e., external wall resistance Rext1 in Fig. 2). In Andrade-Cabrera et al. [5], it was shown that the parametric growth of the variable building parameters can take the shape of an exponential curve. First, a baseline (i.e., pre-retrofitted model) is computed using Particle Swarm Optimisation (PSO), a heuristic optimisation algorithm, using the methodology described in AndradeCabrera et al. [4]. The optimally calibrated parameters include the set of fixed parameters and the baseline value of the variable parameter R∗ext . The model is calibrated using synthetic data gener10

ated with the detached house archetype developed in the EnergyPlus BEMS environment [46]. Then, external wall insulation xext is progressively added to the BEMS (e.g., in 10 mm. steps). Synthetic BEMS data (i.e., heating load and internal temperature time-series) is created at each iteration. The heating load is dynamically calculated using EnergyPlus and the heating systems are auto sized at each iteration. Fig. 3(a) shows the annual average room temperature due to retrofit measures (e.g., external wall insulation of thickness x), which is defined as

Tavg,data (x ) =

k=NH k=0

Tr,x, − Tr,0,k NH

(1)

where Tr, 0, k is the average annual room temperature pre-retrofit and NH the simulation horizon. Fig. 3(a) shows that as a result of the implementation of external wall insulation, the average annual internal temperature increases. The variation of the average annual internal temperature is not significative because of an existing thermostatic setpoint setting around 18.5 °C and the lack of ECMs in other building elements (e.g., thermally inefficient glazing and uninsulated ceiling). Fig. 3(b) shows a decrease in the annual heating load as a result of addition of external wall insulation. The reduction in annual heating load saturates after 100 mm (i.e., the addition of external wall insulation does not result in a significative reduction). Unlike in steady-state residential demand models

(where the load decrease would be proportional to the increment in external wall insulation), a BEMS environment (and thus the calibrated model) account for other heat exchange possibilities (e.g., infiltration, heat loss through uninsulated surfaces), which reduce the effectiveness of higher levels of external wall insulation above 100 mm. For each iteration, a new value is calculated for the set of variable parameters (e.g., R∗ext ) using a gradient descent 1

xext =10 mm

algorithm [43]. As a result of the procedure, the variable parameter Rext1 evolves as a curve, which can be approximated using an exponential function (Fig. 3(c)). The optimally calibrated model parameters were identified using a standard implementation of PSO [44] as a solver and thus the solution is subject to random initialisation and parameter evolution. Furthermore, the calibrated values of modelled element resistances and capacitances will alter reciprocally during model calibration, as they are related by the products Rx Cy . While the model parameters retain their semi-physical nature (i.e., a wall resistance is associated with heat exchanges through external wall), the model becomes a numerical ‘grey-box‘ model [13], which may differ from the theoretical model parameters. A key modelling principle of the Ensemble Calibration methodology consists of the modelling of the thermal effect of a given ECM as a linear function of the baseline model parameter (i.e., pre-retrofit wall resistance) and an exponential approximation (denoted retrofit function). For example: the external wall resistance of a building dwelling as function of external wall insulation of thickness xext can be represented by the retrofit function

Rext1 (xext )  Rext1 0 + αext e(1−βext )xext , xext ≥ 0 Rext is 10

(2)

where the optimally calibrated pre-retrofit external wall resistance, and α ext and β ext are calibration parameters. Thus, all retrofitted building models identified using the Ensemble Calibration framework can be deemed to be functional additive variations of the baseline model building model. The procedure can be automated by means of simultaneously identifying the baseline parameters and retrofit functions. Ensemble Calibration consists of the simultaneous solution of the above-mentioned problem for all

6

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

Fig. 3. Evolution of average annual internal temperature (a), annual heating load (b) and calibrated external resistance (c) as a function of insulation thickness [6].

possible ECM opaque envelope measures. Ensemble Calibration results in the identification of linear, lumped parameter models with a Mean Average Error (MAE) of less than 0.5 K, compared to the synthetic data generated using the associated BEMS archetype. This metric has been suggested in the literature as an acceptable calibration accuracy by He et al. [31], Kusiak and Xu [37] and Lü et al. [41]. Most importantly, the Ensemble Calibration methodology results in a number of lumped parameter building models with shared parameters (i.e., the baseline model). Therefore, the corresponding discrete-time linear building models are linearly dependent with respect to the discrete-time baseline model. This linearity is a key requirement of a building-to-grid co-optimisation model used to assess optimal, large-scale ECM building configurations. The aim of the current paper is to exploit this property, unique to Ensemble models, in order to formulate a computationally efficient solution to the integrated cost-optimal retrofit decision-making and power systems optimisation problem which substitutes the need for a brute force approach. 3. Methodology 3.1. Methodology overview The methodology pursued in the current paper is shown in Fig. 4. It is assumed that an Ensemble model has been obtained using the methodology outlined in Andrade-Cabrera et al. [5]. Section 3.2 formulates Ensemble models as a single bi-linear model via an algebraic manipulation. Section 3.3 introduces the optimisation methodologies necessary to implement the linearisation algorithm. Section 3.3.1 describes the calculation of the business-asusual solution (i.e., the costs prior to retrofit). Section 3.3.2 describes the calculation of NPV returns for a given ECM configuration. A brute force approach consists of calculating the NPV returns for all possible ECM configurations and selecting the ECM configuration with the best solution (e.g., highest returns/savings post ECM investment). Section 3.3.3 describes the proposed multistage optimisation algorithm which adopts the bi-linear building modelling framework (Section 3.2) as part of a MILP optimisation framework. Section 3.4 describes the application of the abovementioned methodologies to a case study where ECM decisionmaking is integrated with electrified space heating technologies exposed to an electricity price signal. Due to space limitations the

integration of a detailed power systems model is left for a subsequent publication. 3.2. Bi-Linear building modelling of Ensemble models 3.2.1. State-space representation of Ensemble models The thermal dynamics of a continuous-time lumped parameter building model can be expressed as a continuous-time state-space model

x˙ (t ) = A(R, C ) x(t ) + Bu (R, C ) u(t ) + Bd (R, C ) d (t )

(3)

where A, Bu and Bd are the continuous-time system, input and disturbance matrices, x(t) is the state evolution (e.g., room and wall temperatures), u(t) is the controlled input (e.g., heating power input) and d(t) is the vector of uncontrolled (i.e., disturbance) inputs such as weather inputs and heat gains [26], [30]. The system matrices are functions of the lumped resistances R and lumped capacitances C, which are often deemed to be constant parameters. In the context of Ensemble Calibration (Eq. (2)), the thermal dynamics of a retrofitted model with external thickness xext can be expressed as a linear combination of the pre-retrofit model A0 , Bu, 0 and Bd, 0 and the variation in dynamics Ax , Bu,x and Bd,x , ext ext ext which in physical terms represent the variation in thermal dynamics associated with ECMs (i.e., the thermal effect of the variable parameters V). The building model takes the form





x˙ (t ) = A0 (R, C ) + Axext (R, C ) x(t ) + +







Bu,0 (R, C ) + Bu,xext (R, C ) u(t )



Bd,0 (R, C ) + Bd,xext (R, C ) d (t )

(4)

This formulation is valid for continuous-time domain models (which are the ones identified using Ensemble Calibration). However, building energy models must be expressed in the discretetime domain for building-to-grid integration. The discretisation of continuous-time dynamic systems requires the use of exponential integrals (Van [40]), which are non-linear and therefore the retrofit functions are no longer valid in the discrete-time domain. A discrete-time building energy model (i.e., the discretised continuous-time model shown in Eq. (3)) can be represented as

xk+1 = Fi xk + Gi uk + Hi dk

(5)

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

7

Fig. 4. Methodology Overview.

In a building-to-grid context, the heating power input uk becomes the electrified space heating load, which is an input to the power systems optimisation model. It can be argued that, in discrete time, the expanded representation of a post-retrofit building energy model can take the form

selected) and thus deemed to be a vector of binary decision variables subject to the constraint

xi,k+1 = F0 xk + G0 uk + H0 dk + F¯i xk + G¯ i uk + H¯ i dk

In a building-to-grid modelling context, the variables xk (i.e., the thermal state), uk (the electrified space heating input) and the retrofit decision variables δ i become decision variables. Hence the products xk δ i and uk δ i are bi-linear. Thus the proposed formulation results in a bi-linear Mixed Integer Non-Linear Program (MINLP), which is potentially intractable [1]. One alternative to deal with this bi-linearity is to reformulate the model as a function of a linearisation trajectory. This approach was used in Lehmann et al. [38] in the context of building predictive control for the solution of blind control problems formulated using integer blind positions. Assuming a pre-computed linearisation trajectory (i.e., the solution of a building-to-grid model for a given ECM model), the bi-linear model described in Eq. (6) becomes

(6)

where xk is the state vector, uk is the electrified heat input and dk the disturbance vector. The state matrices F0 , G0 and H0 define the dynamics of the baseline model and the matrices F¯i , G¯ i and H¯ i are variations in dynamics such that

F¯i = Fi − F0 G¯ i = Gi − G0 Hi = Hi − H0

(7)

Hence, while the continuous-time retrofit functions are not transportable to discrete-time domain, there exists matrix transformations capable of expressing the desired aim (i.e., a discretetime representation of the thermal effects of ECMs in a buildingto-grid scenario) because every retrofitted building energy model is related to the baseline building model.



ncomb

δi = 1

xk+1 = F0 xk + G0 uk + H0 dk



3.2.2. Bi-Linear building energy model formulation A bi-linear problem [36] is a special class of non-linear problem where some of the decision variables are linear with respect to each other (e.g., f (x, y ) = xy). Under this framework, all possible expanded building models of a retrofitted building (Eq. (6)) can be expressed using a single bi-linear formulation

xk+1 = F0 xk + G0 uk + H0 dk



+



ncomb i=1

F¯i xk δi +



ncomb i=1

G¯ i uk δi +



ncomb

 H¯ i dk δi

(8)

i=1

where the ith decision variable δ i represents one of the ncomb possible ECM configurations. The decision variables δ must be defined to be mutually exclusive (i.e., only one ECM configuration can be

(9)

i=1

+



ncomb i=1

F¯i xlin δi +



ncomb i=1

G¯ i ulin δi +



ncomb

 H¯ i dk δi

(10)

i=1

and thus all possible model dynamics variations (e.g., F¯i xlin ) become constant vectors. Therefore, the problem becomes a MixedInteger Linear Problem (MILP) and thus it can be integrated with grid models in a tractable manner via the electrified heating load input uk . 3.3. Cost-Optimal ECM decision-making via linearization 3.3.1. Energy cost minimization and Business-as-Usual costs The cost-optimal solution of an integrated building-to-grid model can be described in its simplest expression as the minimi-

8

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

sation of the energy costs

Energycosts =

NH 

Pk ck

(11)

k=1

where NH corresponds to the simulation horizon (in hours), Pk is the electrified space heating load and ck are the electricity costs. In a fully developed building-to-grid model, the energy costs will have to account for power systems generation costs and therefore the price signal would be replaced by a cost-optimal power systems optimisation problem (e.g., cost-optimal generation). Due to space constraints the current paper uses a fixed price signal as a surrogate of a fully detailed power systems model. In the current paper emphasis is placed on the potential computational advantages of the proposed methodology, while leaving the detailed description of power systems modelling for a subsequent publication. The electrified space heating load Pk (Eq. (11)) is related to the demand side (i.e., the building energy models) by the relationships

xk+1 = Fi xk + Gi Pk + Hi dk ,

(12)

Tr,k = Z xk ,

(13)

SPk ≤ Tr,k

(14)

where Fi , Gi and Hi correspond to the dynamics of a pre- or postretrofit building energy model. Eq. (12) implies that the convective heating load uk of the building model described in Eq. (5) equals the electrified heating load Pk . This assumption is valid for a perfectly efficient electrified heat emitter (i.e., direct resistive heater). Otherwise electrified heating systems require their own modelling. Eq. (13) obtains the room temperature Tr at the time step k from the state vector xk via the output matrix Z. Eq. (14) constrains room temperature to be greater than the thermostatic heating setpoint profile SP. A cost-optimal solution of the co-optimised (i.e., integrated) building-to-grid model Eqs. (12) to (14) corresponds to the solution that respects the set-point constraints while minimising the energy costs. The optimal Business-as-Usual costs (BAUcosts ) corresponds to the solution of the integrated building-to-grid energy cost minimisation problem, which minimises the energy costs of the co-optimised building-to-grid model, when the building dynamics (Eq. (12)) correspond to the pre-retrofit building (i.e., the energy costs if ECM measures are not considered). 3.3.2. Cost-Optimal ECM decision-making via brute force approach In the current paper the retrofit decision-making criteria will be the Net Present Value (NPV) returns of the capital investment in envelope ECM costs Envelopecosts and heating technology costs Techcosts . The NPV returns of the ith ECM configuration are expressed by the relationship Nyears

NPVi =

 BAUcosts − EC Mcosts,i y=1

(1 + r )y−1

− (Envelopecosts,i + T echcosts,i ) (15)

where Nyears corresponds to the investment period (in years), BAUcosts corresponds to the Business-as-Usual energy costs, ECMcosts, i corresponds to the energy costs for a given ith ECM configuration (i.e., one of the retrofitted models extracted using the Ensemble Calibration approach) and r is the discount rate for the investment. In the current study a discount rate r of 4% is selected. This discount rate is used by the Irish Government in cost-benefit analysis of public sector projects [18] and thus is suitable for the analysis of retrofit investments at scale as

shown in AECOM [2]. The ECM investment will yield positive NPV returns if the discounted sum of the energy costs savings (i.e., the difference between CostBAU and CostECM ) is sufficient to offset the capital investment in ECM measures and/or space heating technology over the investment period. A negative value implies a loss on the ECM capital investment. The energy costs of a given ECM configuration are also a function of price signal ck and the post-retrofit electrified space heating load Pˆk

EC Mcosts =

NH 

Pˆk ck

(16)

k=1

and the post retrofit space heating load Pˆk is subject to the building dynamic equations Eqs. (12)-(14) for the case where the model Fi , Gi and Hi corresponds to the building thermal dynamics of the ith ECM configuration for which the NPV return (Eq. (15)) is being calculated. A straightforward approach towards the solution of an integrated ECM decision-making problem and power systems optimisation problem consists of identifying the NPV of every possible ith ECM configuration and selecting the ECM configuration that maximizes the NPV. This approach is denoted as the brute force solution of the ECM retrofit decision-making problem. 3.3.3. ECM decision-making using a multi-stage optimisation algorithm One alternative to improve the computational time consists of introducing the linearised bi-linear model (Eq. (10)) as the building dynamics of the NPV retrofit decision-making problem. The NPV formulation becomes Nyears

NPV =

 BAUcosts − EC Mcosts,i

(1 + r )

y=1

y−1



ncomb

+

EnvelopeCosts,i δi + T echcosts,i

i=1

(17) and the building model dynamics become

xk+1 = F0 xk + G0 Pˆk + H0 dk



+



ncomb i=1

F¯i xlin δi +



ncomb i=1

G¯ i ulin δi +



ncomb

 H¯ i dk δi

(18)

i=1

In doing so, the supply side mode (i.e., the power system model) is solved only once, thus reducing the complexity of the solution if the power systems model is higher than the complexity of the demand side. At optimality, the linearisation trajectory is equal to the optimal state and heating inputs (i.e., xlin = x and ulin = u ). However, the electrified space heating load (and thus the power systems optimisation problem) is bound to be suboptimal due to linearisation errors (ulin ≈ u ). 3.3.3.1. Limitations of the bi-linear modelling approach. The linearised bi-linear building model requires a pre-computed linearisation trajectory for state evolution and heating load input (xlin and ulin , respectively). These linearisation trajectories can be extracted by means of solving the NPV problem (Eq. (15)) for a predefined initial point {Mx , Nx , Ox , Px } and extracting the optimal calculated state x and heating load u . However, the need to select a linearisation trajectory introduces a building physics constraint. Fig. 5 shows the temperature response of two thermal nodes associated with the ECM configurations {1,1,4,4} (i.e., 300 mm of ceiling insulation and triple glazing) and {3,1,3,2} (i.e., 200 mm of external insulation, 200 mm of ceiling insulation and double glazing with uPVC frames) of the lumped parameter house archetype model described in Andrade-Cabrera et al. [5]. These configurations have been selected as they represent buildings retrofitted with and

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

9

Fig. 5. Linearisation trajectories.

without external wall insulation. Fig. 5(a) shows that the thermal response of the room temperature (thermal node Tr ) is similar for the two selected retrofit options, given the existence of a thermostatic constraint (Eq. (14)). On the other hand, Fig. 5(b) shows that the average internal wall temperature (i.e., node Tw2 ) is significantly different between both retrofitted scenarios. The reason for this discrepancy lies on the formulation of the Ensemble Calibration methodology. Recall that thermal wall resistances in an Ensemble Calibration framework are identified with the form Rext1 (xext )  Rext + αext e(1−βext )xext (Eq. (1)). That is, an additive 10 formulation of associated wall resistance (pre-retrofit value Rext ) 10

and post-retrofit contribution (or model variation) αext e(1−βext )xext . In Andrade-Cabrera et al. [5], it was shown that external wall insulation needed to be captured using the thermal resistance Rext2 due to the sensitivity of the thermal resistance Rext1 (which is subject to fast dynamics). The total external wall resistance Rext (i.e., inverse of the wall transmittance Uwall ) can be seen, from a modelling perspective, as the addition of thermal resistances

Uwall =

=

where



Algorithm 1 Linearisation Algorithm. Result: Optimal space heating load u and decision vector δi Initialisation: Define the initial points {M1 ,N1 ,O1 ,P1 } and {M2 ,N2 ,O2 ,P2 } 1 Extract the business-as-usual costs BAUcost for the pre-retrofit model {M0 ,N0 ,O0 ,P0 } (Section 3.3.1) 2 Compute the cost-optimal building-to-grid model with assumption {M1 ,N1 ,O1 ,P1 } (Section 3.3.1) Extract xlin,1 and ulin,1 (Linearisation trajectory for solutions without external wall insulation) 3 Compute the NPV return NPV1 using BAUcost , xlin1 and ulin1 (Section 3.3.3) 4 Compute the cost-optimal building-to-grid model with assumption {M2 ,N2 ,O2 ,P2 } (Section 3.3.1) Extract xlin,2 and ulin,2 (Linearisation trajectory for solutions with external wall insulation) 5 Compute the NPV return NPV2 using BAUcost , xlin2 and ulin2 (Section 3.3.3) 6 Compare NPV1 and NPV2 . Select the solution δ  which provides the highest NPV return. 7 Obtain the solution of the space heating load u by computing the optimal building-to-grid solution (i.e., solve the Section 3.3.2) using the optimal configuration δ  .

1

Rxext + R∗ext1

R∗ext1 0

R∗ext 1

, 0

0

+ R∗ext2 + R∗ext3 0

1



+ Rxext + R∗ext2 R∗ext 2

and 0

R∗ext 3

0

0

(19) + R∗ext3

0

are the optimally calibrated values of 0

the pre-retrofit external wall resistances of the model shown in Fig. 2. In the case of Fig. 5(b), the values of R∗ext , R∗ext and R∗ext 10

20

30

(i.e., Rxext = 0 for {1,1,4,4}) are considerably lower than the value of Rxext when xext = 200 mm. As a result, the optimized model with configuration {3, 1, 3, 2} (which features Rxext ) has a significantly different temperature profile than the temperature profile of the configuration {1,1,4,4} (which does not). Therefore, the internal wall temperature profile Tw2 is considerably different between both scenarios because only the second scenario features the parametric variation associated with external wall insulation. This effect is perfectly natural in the context of building physics and retrofit practice (i.e., the thermal node associated with insulation is storing heat). Both nodes are part of a candidate linearisation trajectory (i.e., Tr and Tw2 are part of the linearisation trajectory, alongside the other node temperatures). In order to guarantee local optimality, the trajectories should lie sufficiently close to the optimal solution [56]. Both linearisation trajectories are deemed to be significantly different, even though the thermal room response is similar enough in both ECM configurations. Therefore, if the optimal solution features external wall insulation, optimality is guaranteed only

if the linearisation trajectory features external wall insulation. The same applies for an optimal solution featuring internal wall insulation. In numerical terms this implies that two linearisation points must be considered in order to obtain a globally optimal solution: One initialisation trajectory {M1 , N1 , O1 , P1 } where the ECM configuration does not include external wall insulation, and another linearisation trajectory {M2 , N2 , O2 , P2 } where the ECM configuration includes external wall insulation. Otherwise, it is unlikely that the optimal ECM configuration can be found if the linearisation trajectory lies far from the optimal solution. 3.3.3.2. Linearisation algorithm. Algorithm 1 proposes a linearisation algorithm which identifies the optimal solution of the bilinear building-to-grid retrofit investment planning problem for an individual archetype by means of correcting this building physics constraint (i.e., two solutions) using a heuristic method. First, the business-as-usual costs BAUCosts are computed (Section 3.3.1). Then, the linearisation trajectory is extracted for the initial guess {M1 , N1 , O1 , P1 }. This initial guess corresponds to the solution sub-space which includes all models without external insulation. This step obtains the linearisation trajectory xlin, 1 and linearisation heating load ulin, 1 , given the computed BAUCosts . The business-as-usual costs and the linearisation trajectory and heating load are then used to compute the linearised solution of the integrated ECM decision-making problem (Section 3.3.3). As a result, the model calculates the NPV returns with respect to the linearisation trajectory (denoted NPV1 ). The procedure is repeated for the linearisation point {M2 , N2 , O2 , P2 } (i.e., the initial guess for the lineari-

10

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833 Table 1 Retrofit Costs per ECM. ECM

Thickness mm

Target U-value W/m2 K

2013 Cost EUR/m2

Adjusted Cost EUR/m2

Internal Insulation (EPS) Internal Insulation (EPS) External Insulation (EPS) External Insulation (EPS) Ceiling Insulation (Glass wool) Ceiling Insulation (Glass wool) Double Glazing (uPVC windows) Double Glazing (Argon filled uPVC windows) Triple glazing

50 100 100 200 200 300 – – –

2.38 → 0.38 2.38 → 0.21 2.38 → 0.21 2.38 → 0.15 0.15 → 0.13 0.15 → 0.11 3.2 → 1.6 3.2 → 1.4 3.2 → 0.9

43.1 72.6 105.7 118.9 3.7 6.8 328 356 483

50.6 85.3 124.2 139.7 4.3 8.0 385.4 418.3 567.5

sation trajectory for the subset of models featuring external insulation). The procedure results in a second NPV returns estimate (denoted NPV2 ). It suffices to select the solution with the highest NPV return. As mentioned in Section 3.4.3, this procedure identifies the optimal ECM configuration δ  for a sub-optimal heating load for a sub-optimal heating load estimation (and heating technology choice) because of the use of a linearisation trajectory. The solution is refined by means of solving the problem for an individual ECM configuration (Section 3.3.2) using the optimally calculated ECM configuration δ  . The computations were carried out using a Dell Server PER 730 server with 2 Intel Xeon E5-2699 processors at a frequency of 2.20 GHz and 128 GB of memory. The optimisation problems solved in the current section were implemented in Pyomo [29] and the MILP solver used in the current paper was Gurobi [28]. 3.4. Case study: including storage heating technologies in ECM decision-making 3.4.1. Problem definition The current section identifies the potential of electrified storage heating technologies as part of comprehensive building envelope and heating technology ECM configurations. The novelty consists of the inclusion of storage technologies in retrofit scenarios where the heating load is reduced and therefore the monetary value of load shifting is also reduced. It will be assumed that storage heaters (SH) can replace the existing direct resistive heaters (DRH) without the need for significant modifications to the domestic heating system. The current section focuses only on the detached houses that already feature DRH as space heating technology. Heat Pumps (HP) are known to be an attractive solution towards space heating electrification [32]. However, HP require the installation of hydronic systems in a building, which is a cost intensive and invasive procedure. Thus, HP are not a pragmatic solution for this subset of the residential stock (i.e., houses without prior hydronic installation). The number of detached houses nhouses is assumed to be 5843 (0.35% of the total Irish residential stock, which includes nearly 1.7 million dwellings). This estimation was obtained from census data [20]. The current section examines the potential advantages of adopting SH over DRH in this limited subset of the building stock, while demonstrating the core characteristics of the proposed multi-stage optimisation algorithm (i.e., algorithmic convergence and computational efficiency). 3.4.2. Building modelling assumptions 3.4.2.1. Building envelope ECM measures. Table 1 describes ECM measures and their corresponding capital costs considered in the current paper. The information was obtained from the CostOptimal Calculations and Gap Analysis for Buildings in Ireland [2]. The study presents capital costs per ECM measure which include both materials and labour. The costs shown in Table 1 have been adjusted to 2018 prices under the assumption of a 17.5% increase in building costs between

Table 2 ECM measures associated with Ensemble Calibration Indexes. Index

ECM

Units

1

2

3

4

M N O P P

External (EPS) Internal (EPS) Ceiling (Glass wool) Glazing (U-value) Glazing (g-value)

mm mm mm m2 W/K −

0 0 – 3.2 0.76

100 50 – 1.6 0.62

200 100 200 1.4 0.62

200 100 300 0.9 0.48

2013 and 2018 [39]. The net external wall area is 148 m2 , the floor/ceiling area is 80 m2 and the window area is 35.6 m2 . The infiltration rate reduction associated with each ECM is described in Andrade-Cabrera, Turner and Finn [6]. Cost-intensive ECMs such as mechanical ventilation heat recovery (MVHR) are unlikely to be adopted at large scale, as shown in AECOM [2], and thus not considered in the current study. Ensemble models were extracted for this model using the ECM configurations described in Table 2. For ease of notation, the current paper adopts the Ensemble indexes notation introduced in Andrade-Cabrera et al. [5]. Table 2 denotes the relationship between Ensemble indexes and the corresponding ECM measures. The possible values of external and internal wall insulation are denoted using indexes M and N, respectively. Ceiling insulation is denoted using the index O and window replacement is denoted using the index P. A building retrofit combination can be denoted by the index tetrad {M, N, O, P}. This short-hand notation enables the simplified description of retrofit combinations in the current paper. For example, an index {1,1,4,4} denotes an ECM with 300 mm. of ceiling insulation (O = 4) and glazing replacement using triple 2 glazing (P = 4, U-value = 0.9 mKW , g-value = 0.48). Simultaneous external and internal wall insulation scenarios are excluded as they are generally not consistent with retrofit practice. Furthermore, only ceiling insulations of 200 mm and 300 mm will be considered, as it is common engineering sense that these insulation measures are likely to be selected first, given their low capital cost. These modelling choices reduces the number of ECM configurations from a potential ncomb = 44 = 256 configurations (i.e., all possible permutations) to a reduced set of ncomb = 38 ECM configurations. 3.4.2.2. Heating systems. In the current paper, direct resistive heaters are assumed to be perfectly efficient and thus the space heating consumption Pk is equal to the heating load uk . In the case of storage heaters, the Energy level Ek of the storage heater is constrained by the energy balance

Ek+1 = Ek + Pk − Qheat,k − Qloss.k

(19)

where Qheat, k is the heat delivered to the conditioned space, Sk is the space heating power consumption and Qloss.k represents the thermal losses of the storage device, which are modelled as

Qloss.k = (1 − η )Ek

(20)

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

11

Table 3 Maximum Operational Parameters and Cost per Electrified Heating Technology. Device Model

Technology

Smax [W]

Qmax [W]

Emax [kWh]

Device Cost [EUR]

Installation Cost [EUR]

T1 T2 T3 T4 T5 T6

DRH SH SH SH SH SH

3000 1100 1560 2200 2760 3300

3000 500 700 1000 1250 1500

– 8.5 10.9 15.4 19.3 23.1

0 577 641 695 758 810

0 200 200 200 200 200

where η is the thermal efficiency of the storage heater (assumed to be η =0.975 as per manufacturer specifications). Refined storage efficiency models (e.g., temperature dependent efficiencies) are not considered in the current paper in order to preserve a linear modelling framework. The space heating power consumption Sk , heat output rate Qheat, k and energy level Ek are constrained by the rating levels Qheat, max , Smax and Emax , respectively. Table 3 describes the device rating and the pricing information for each of the heating technologies considered for the current study. Technology T1 corresponds to 1.2 kW direct resistance electric panel heater (assumed to be already installed in the dwelling) as per archetype specification [46]. The detached house EnergyPlus archetype had a peak load of approximately 15.6 kW. Each room was deemed to feature a DRH device. The detached house archetype features 13 direct resistive heaters as per the original archetype [46]. Technologies T2 to T6 correspond to storage heaters. The retail price and installation costs were provided by the manufacturer in April 2018. The devices are at the upper range of market value for storage heating technologies. However, the devices are deemed as capable for full integration with IoT aggregation technologies, enabling full thermal mass exploitation, thus maximizing demandside flexibility. In the current paper, it will be assumed that the direct resistive heating load is not flexible (i.e., heating is provided when required). The proposed methodology is tested using the detached house and mid-floor apartment archetypes described in Andrade-Cabrera et al. [5]. The current paper does not incorporate a fully detailed power systems model as its intension is to focus on the performance of the optimisation algorithm. Hence a price signal is used as a surrogate of a fully detailed power systems model. Electricity prices will be modelled using a night-savings tariff where the price of electricity is deemed to be 9.91 ¢EUR/kWh between 11 PM and 6 AM (i.e., during night hours) and 20.03 ¢EUR/kWh during the rest of the day [57]. Both the brute force approach (Section 3.4.2) and the linearisation approach (Section 3.4.3) are used to solve the integrated envelope ECM decision-making problem with storage heating. The current paper does not intend to show a straightforward what-if scenario comparison (i.e., what if storage heating replaces direct resistive heating), which would be trivial given the capital cost and installation costs shown in Table 3. In the current paper, the optimisation solver has been allowed to select the number of storage heaters and direct resistive heaters for a cost-optimal configuration, assuming pre-existing direct resistive heaters. This hybrid configuration allows for a cost-efficient partial update of the heating systems while still providing savings that are attractive to the end-user, bearing in mind that investors will seek to maximize their economic self-interest as discussed in Aravena et al. [8]. This case study shows that the flexible formulation of the integrated ECM decision making and power systems optimisation (i.e., storage heating) can provide benefits to the user (i.e., higher NPV returns) but also provide societal benefits (e.g., increased decarbonisation) and potentially provide benefits to the power systems sector (e.g., reduced operating costs).

Fig. 6. Comparison of ECM costs against CO2 emissions (detached house archetype).

4. Results 4.1. Optimal ECM decision-making using a brute force approach Figs. 6 and 7 summarises the results of using the brute force approach to identify the NPV returns for each possible combination of the ECM configurations assigned to the detached house archetype (Section 3.2.2). Fig. 6 compares the total ECM investment costs (in EUR) against the CO2 emissions associated with the ECM configuration (in tCO2 per year). Carbon emissions were calculated using a conversion factor of 483.8 gCO2/kWh [55]. The calculated emissions correspond exclusively to space heating only and thus they are directly proportional to the heating load of a given ECM configuration. Fig. 7 compares the total ECM investment costs against NPV returns, which are positive if the investment results in savings or negative if the ECM configuration results in a loss. The total ECM investment costs include the building envelope ECM costs and the technology costs (e.g., replacement of direct resistive heating by a storage heating unit). The colour bar denotes the number of storage devices that were selected by the cost-optimal retrofit tool. The colour bar shows that up to 4 storage heaters have been selected to replace direct resistive heaters for some ECM configurations. Thermal comfort is satisfied for all configurations as per the thermal comfort constraints (Eq. (14)). The cost-optimal solution that maximizes the NPV returns (i.e., savings) is {1,2,3,1} (i.e., 50 mm of external insulation and 200 mm of ceiling insulation, as shown in Fig. 7). The cost-optimal ECM

12

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

Fig. 7. Comparison of ECM costs against vs NPV returns (detached house archetype).

configuration results in NPV savings of 40,870 EUR over the investment period (i.e., 30 years), with a capital investment of 10,862 EUR. From a carbon emissions perspective, the emissions generated by the cost-optimal configuration (4.76 tCO2 per year) result in an abatement of 60% with respect to the pre-retrofit model, which generated 12.17 tCO2 per year. Given the thermal inefficiencies of the storage heater (ηSH ), the loads of retrofitted dwellings featuring storage heating are higher than the heating loads without storage heating. Therefore, under the assumption of a fixed conversion factor for carbon emissions (482.8 gCO2/kW), ECM configurations with storage heating are expected to yield higher carbon emissions than the equivalent ECM configurations without storage heating. Such an assumption does not account for the exploitation of the stored energy to provide grid side flexibility, enabling enhanced integration of wind energy. Hence the need to demonstrate the true decarbonisation potential of storage heating (and other space heating technologies) using building-to-grid models. The ECM configuration {1,2,4,1} (i.e., 50 mm of external insulation and 300 mm of ceiling insulation) lies adjacent to the cost-optimal solution in the diagram shown in Fig. 7. This sub-optimal solution results in discounted NPV savings of 40,813 EUR (i.e., a difference of 59 EUR between the cost-optimal and the sub-optimal configurations). Both configurations adopt the same heating technology configuration. Only 3 out of 13 DRH have been replaced using high capacity storage heaters T6 . The results shown in Fig. 7 are comparable with previous studies of residential retrofit policies in Ireland. The target configuration described in AECOM [2] was {1,3,4,3} (i.e., 100 mm of internal insulation, 300 mm. of ceiling insulation and window replacement using argon-filled uPVC glazing). The estimated retrofit costs for this configuration is 30,175 EUR, which is about four times higher than the cost-optimal configuration, but yields a 73% reduction in emissions with respect to the pre-retrofit scenario (3.22 tCO2 per year). Ahern, Griffiths, and O’Flaherty [3] considered four detached house configurations: Single or double glazing and one or two storeys. In the current study it is possible to make a comparison only with a two-storey house with double glazing, which is directly comparable to the model described in Section 3.4.2 (i.e.,

the detached house under study). The cost-optimal ECM configuration described in Ahern, Griffiths, and O’Flaherty [3] corresponds to a reduction of external wall insulation to a target level of Uwall = 0.27 W/m2 K and ceiling insulation to a target level of Uroof = 0.18 W/m2 K. Using the target level specifications described in AECOM [2], it follows that the associated ECM configuration would be {3,1,3,1} (e.g., 100 mm. of external insulation and 200 mm. of ceiling insulation). The investment costs associated with this configuration were estimated at 16,431 EUR (house types A-E, 2013 costs), whereas in the current study the costs were calculated to be 24,049 EUR, which includes the costs of two highcapacity storage heaters Q6 (valued at 1010 apiece). The emissions associated with this ECM configuration (5.62 tCO2 per year) represented a 65% carbon abatement, as shown in Fig. 6. Highly efficient ECM configurations represent an attractive investment from the emissions abatement perspective (e.g., when retrofits are considered at the environmental policy level). For example, a highly insulated ECM configuration such as {3,1,4,4} (i.e., 100 mm. of external insulation, 300 mm. of ceiling insulation and triple glazing windows) emits 2.9 tCO2 per year (i.e., a 76% carbon abatement). However, from the end-user perspective such a configuration is far from the cost-optimal solution, as it represents a prohibitive investment of 43,538 EUR, which is four times the value of the cost-optimal retrofit costs (10,862 EUR). Inefficient ECM configurations that do not feature either external insulation or internal insulation (e.g., {1,1,3,2}, 200 mm of ceiling insulation, double glazing uPVC windows) are thermally inefficient given the heat losses through the external wall surfaces. These thermally inefficient dwellings use 5 high capacity storage devices. The use of storage devices in these poorly insulated scenarios results in a profitable (albeit sub-optimal) investment (NPV savings of 21,407 EUR over a 30-year investment period) but with reduced carbon abatement potential (6.99 tCO2 per year, meaning a 42% reduction of carbon emissions). The heating technology costs (4040 EUR) represent a fifth of the total retrofit costs (18,104 EUR). Tables 4 and 5 summarise the findings discussed in the current section. As mentioned earlier in this section, emissions are deemed to be linearly proportional to the electrified heating load (i.e., an emissions factor is used). However, in an integrated building-togrid context (which forcefully reduces carbon emissions as part of the total generation costs), it is possible that the exploitation of the supply-side variabilities such as wind and solar power availability result in a higher carbon abatement when future electricity generation scenarios (e.g., higher levels of wind penetration) are considered. Hence the importance of integrated, cost-optimal ECM decision-making algorithms that can be integrated with power systems optimisation models.

4.1.1. Mid-floor archetype model Figs. 8 and 9 shows the results obtained from calculating the NPV returns for the mid-floor archetype apartment using a brute force approach. The pre-retrofit mid-floor apartment features four direct resistive heaters as per the original archetype [46]. Fig. 8 shows that at least one storage heater has been selected for 3 out of 19 possible ECM configurations. The ECM configurations which have higher load requirements post-retrofit (and thus higher CO2 emissions) select storage heating as a partial replacement to existing heating systems. Fig. 9 shows that the configurations that correspond to scenarios whereby only glazing has been considered as the retrofit measure are economically unattractive from an NPV returns point of view, with a triple glazing configuration {1,1,1,4} (i.e., only triple glazing is considered) resulting in a loss of investment (−119 EUR). Tables 6 and 7 summarise the findings of the mid-floor archetype.

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

13

Table 4 Comparison of Selected ECM configurations (Detached house archetype). ECM Configuration

NPV Returns [EUR]

ECM Costs [EUR]

Emissions [tCO2/year]

Heating Systems

Cost Optimal ECMs {1,2,3,1} Sub-Optimal ECMs {1,2,4,1} AECOM [2] {1,3,4,3} Ahern et al. [3] {3,1,3,1} Highly Efficient ECMs (e.g., {3,1,4,4}) Inefficient ECMs (e.g., {1,1,3,2})

40,870 40,813 29,629 30,559 17,662 21,407

10,862 11,158 30,175 24,049 43,568 18,104

4.76 (−60%) 4.71 (−61%) 3.22 (−73%) 4.14 (−65%) 2.9 (−76%) 6.99 (−42%)

10 DRH, 3 x T6 10 DRH, 3 x T6 12 DRH, 1 x T6 11 DRH, 2 x T6 11 DRH, 2 x T6 8 DRH, 4 x T6

Table 5 Cost Breakdown for Selected ECM Configurations (Detached House archetype). ECM Configuration Cost Optimal ECMs {1,2,3,1} Sub-Optimal ECMs {1,2,4,1} AECOM [2] {1,3,4,3} Ahern et al. [3] {3,1,3,1} Highly Efficient ECMs (e.g., {3,1,4,4}) Inefficient ECMs (e.g., {1,1,3,2})

ECM Costs [EUR] Total Costs

Wall Insulation

Ceiling Insulation

Glazing Replacement

Heating Systems

10,862 11,158 30,175 24,049 43,568 18,104

7488 (69%) 7488 (67%) 12,624 (42%) 20,675 (86%) 20,675 (47%) 0

344 640 640 344 640 344

0 0 14,891 (49%) 0 20,203 (46%) 13,720 (76%)

3030 3030 2020 3030 2020 4040

(3%) (6%) (2%) (1%) (1%) (2%)

(28%) (27%) (7%) (13%) (5%) (22%)

Table 6 Comparison of Selected ECM configurations (mid-floor apartment archetype). ECM Configuration

NPV Returns [EUR]

ECM Costs [EUR]

Emissions [tCO2/year]

Heating Systems

Cost Optimal ECMs {1,2,1,1} AECOM [2] {1,3,1,3} Highly Efficient ECMs (e.g., {3,1,1,4}) Inefficient ECMs (e.g., {1,1,1,4})

6900 4579 1482 −119

1278 6647 9624 7052

0.87 (−46%) 0.27 (−81%) 0.3 (−83%) 1.08 (−33%)

4 4 4 3

DRH DRH DRH DRH, 1 x T5

Table 7 Cost Breakdown for Selected ECM Configurations (Mid-Floor Archetype House). ECM Configuration Cost Optimal ECMs {1,2,1,1} AECOM [2] {1,3,1,3} Highly Efficient ECMs (e.g., {3,1,1,4}) Inefficient ECMs (e.g., {1,1,1,4})

ECM Costs [EUR] Total Costs

Wall Insulation

Ceiling Insulation

Glazing Replacement

Heating Systems

1278 6647 9623 7052

1278 2154 (32%) 3528 (37%) 0

0 0 0 0

0 4492 (68%) 6094 (63%) 6094 (86%)

0 0 0 9,58 (14%)

Fig. 8. Comparison of ECM costs against CO2 emissions (mid-floor apartment archetype).

Fig. 9. Comparison of ECM costs against NPV returns (mid-floor apartment archetype).

14

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

As mentioned earlier in this section, emissions are deemed to be linearly proportional to the electrified heating load (i.e., an emissions factor is used). However, in an integrated building-togrid context (which forcefully reduces carbon emissions as part of the total generation costs), it is possible that the exploitation of the supply-side variabilities such as wind and solar power availability result in a higher carbon abatement when future electricity generation scenarios (e.g., higher levels of wind penetration) are considered. Hence the importance of integrated, cost-optimal ECM decision-making algorithms that can be integrated with power systems optimisation models. 4.2. Convergence of the multi-stage optimisation algorithm As mentioned in Section 3.3.3, one of the requirements of the algorithm is the selection of initial points for the calculation of linearisation trajectories. An initial point, in the context of this work, corresponds to a possible assumption of an initial point {Mx , Nx , Ox , Px } for which a linearisation trajectory has been calculated. Furthermore, it was also discussed in the same section that there exist two possible solutions: one solution for the ECM configurations featuring external wall insulation and another solution for the ECM configurations without external wall insulation. Tables A1 – A3 (Appendix A) describe the cost-optimal ECM configuration δ  for each possible initial point δ init . The tables also show the NPV returns and the percentage error with respect of the cost-optimal solution (NPVerror ) as well as the heating technology solutions (TH). In the case of the detached house archetype (Table A1), the solution does not converge to the global optimum most of the time. The difference in NPV savings between the optimal solution {1,2,3,1} (50 mm of internal wall insulation, 200 mm of ceiling insulation) and the sub-optimal solution {1,2,4,1} (50 mm of internal wall insulation, 300 mm of ceiling insulation) is 59 EUR. On a percentage basis, the difference between both NPV savings represents 0.14% of the optimal NPV returns (40,871 EUR). In optimisation terms this means that the difference between both candidate solutions falls below the pre-programmed MIP gap of 1%. In integer programming, the MIP gap is the tolerance that mixed-integer programming allows between solutions to certify optimality [28]. A model with an initial point that does not include internal insulation will always converge to the cost-optimal solution {1,2,3,1} or the marginally sub-optimal cost-solution {1,2,4,1}, as shown in Table A1. On the other hand, all models with an initial point including external wall insulation (Table A2) converge to the same solution. The linearisation errors result in cost minimisation errors, which in an integrated building and grids framework will translate into erroneous power system configurations. Hence the need for a solution refinement stage (i.e., Step 7 in Algorithm 1) in order to obtain the optimal power systems model for the optimally identified ECM configuration. The same issue does not reproduce in the case of the midfloor apartment archetype (Table A3), whereby all models converge to the optimal solution of their solution space (i.e., a model with initial point with internal/external insulation will converge to the same point). 4.3. Computational advantages of the linearisation algorithm Fig. 10 summarises the computational performance of the linearisation algorithm (Algorithm 1, Section 3.4.3) for a scenario where each archetype is solved as a separate ECM decision-making problem and for a second scenario where both models are integrated. In the case of the detached house archetype, the solution using a brute force approach required a computational time of 47.56 min. On the other hand, the linearisation approach optimally solved this problem in 7.05 min, which represents a speedup

of 6.74 times. In the case of the mid-apartment the solution via brute force was solved in 52.3 min whereas the linearisation algorithm solved the same problem in 20.54 min (a speedup of 2.54 times). One possible explanation of this computational reduction is that even though the apartment archetype is of a higher dimension (i.e., 7 states, compared with the detached house archetype which has 6 states), the number of combinations is inferior (19 combinations in the mid-floor apartment vs. 38 combinations in the detached house) and thus the methodology does not result in as large a speed-up. The solution times reported for the linearisation algorithm includes all the required steps (i.e., executing the linearisation trajectory twice plus the solution of the refination step). Assuming model decoupling (i.e., no interactions between models), the solution of both problems is the sum of the solutions and hence the total computational time required to solve the problem is 201.7 min (i.e., above 3 h) whereas the linearisation approach solves the same problem in 47.56 min (i.e., a speed-up of 4.57 times). The results shown in Fig. 10 assume that the solutions of the algorithms are decoupled with respect to one another. This approach is perfectly valid under the assumption that a fixed price signal (e.g., electricity prices, or other energy prices for different energy vectors) is used in lieu of detailed building-to-grid integration. 5. Discussion and conclusions The current paper proposes an optimisation methodology which provides a computationally efficient solution of the integrated cost-optimal retrofit decision-making and power systems optimisation problem using BEMS archetype models. Such a problem has been addressed in the literature as a sequential process whereby residential electrified heating load time-series are precalculated using BEMS archetypes and then fed to power systems models as in Ault et al. [12]. Additionally, thermal requirements can also be added as constraints of an energy optimisation problem as in Wu et al. [59]. Such approaches fail to fully exploit the potential of building thermal flexibility and supply-side flexibilities (such as storage) given the sequential nature of the proposed approach. The proposed methodology enables researchers to embed building thermal performance representative of BEMS archetypes in a building-to-grid scenario while providing the option of exploiting thermal mass and allowing for the integration of power systems models. The proposed approach opens the possibility for more detailed studies on the impact of integrated flexibilities between the supply side and the demand side. For example: higher forecasted levels of wind penetration in future electricity generation scenarios might result in the adoption of space heating technologies that result in lower energy efficiency specifications, but perhaps larger investments in thermal storage. Due to space constraints the current paper uses a price signal (i.e., electricity night tariff) to represent the potential linear integration with other power systems models that can be expressed as linear models. Energy systems (e.g., power systems models) that can be expressed in a linear form can easily be integrated in the proposed framework by means of using the heating power load Pk (Eq. (11)) as part of the power systems model. The proposed approach presents some resemblance with the methodology introduced in Schutz et al. [53] in the sense that the thermal performance post retrofit is incorporated as an optimisation variable in a mixed-integer scenario. The work of Schutz et al. shows that the methodology provides acceptable energy estimation ranges if a normative model (ISO 13,790) is considered as the basis of the residential model. However, it is also known that the normative models can potentially overestimate the performance of BEMS models by 20–30% [24]. Hence the proposed approach provides a potentially more realistic representation of space heating load in retrofit

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833

15

Fig. 10. Computational performance of the proposed linearisation algorithm versus solution via brute force.

decision-making scenarios and therefore a higher representativity of the integrated building-to-grid model outputs. The proposed methodology provides a computationally efficient solution to the integrated retrofit decision making in the context of building-to-grid simulations. The method is 4.7 times faster than the brute force approach. The results hereby presented correspond to the simplest possible complexity from an energy systems perspective (i.e., a simple price signal). However, in a buildingto-grid scenario, each evaluation of the ncomb ECM configurations for the detached house (38) would also require the power systems model to be recomputed at the same time. Once multiple archetype model integration is considered (i.e., a power systems model replaces the price signal in Eq. (16)), a globally optimal solution of the integrated cost-optimal ECM decision making problem must consider all possible ECM configurations for each archetype, and thus a brute approach problem becomes computationally more challenging. For example: if the two archetypes were integrated, the globally optimal solution of the cost-optimal problem postretrofit (Eq. (15)) needs to be solved for every possible ECM configuration of the first archetype model while exploring every possible ECM configuration of the second archetype model. Solving the problem for all possible ECM configuration combinations of the two archetypes (i.e., 19×38 combinations) resulted in a total computational time of 2555 min (i.e., about two days of computation). In that scenario, the proposed method results in a speed-up of 58 times. Evidently this approach does not scale up well with a large number of archetypes. It is possible to carry out the integration of retrofit decisionmaking and power systems optimisation using a theoretical whitebox models by means of modifying the wall resistance values of retrofitted building energy models using well-known relationships (Gouda, M., Danaher, S. and Underwood, C.P. [27]). Ensemble Calibration relies on grey-box modelling, whereby the model parameters have some physical meaning but are numerical in nature and may differ from the theoretical models [15]. It was shown in Andrade-Cabrera et al. [5] that grey-box models better capture the thermal response of EnergyPlus archetype models, ultimately leading to a significantly higher annual energy estimation accuracy post-calibration. Nevertheless, the linear integration method (which is the core contribution of this paper) is valid for either approach (theoretical white-box models or grey-box mod-

els) as both models result in a similar linear structure described in Section 3.2.1. The fundamental relationship between thermal resistance variations associated with ECM measures (i.e., pre and postretrofit wall resistances) holds true for both white-box and greybox modelling frameworks and therefore the linearisation trajectories (e.g., Fig. 5) are likely to diverge. The global optimality of the proposed approach is constrained to the tolerance of the MILP solver (e.g., the MIP gap). A high tolerance value (e.g., 5%) may result in solutions which are marginally suboptimal whereas a tighter tolerance (e.g., 0.01%) is much more likely to identify the right value at a higher computational cost, especially when power systems models are to be considered. A compromise needs to be found between the optimality conditions stated during the implementation phase (i.e., MILP gap) and the computational resources available to solve the integrated buildingto-grid models. Future work will focus on the exploitation of the optimisation framework introduced in the current paper in the context of examining the potential benefits and limitations cost-optimal residential retrofit policies at scale when electrified space heating is considered as a candidate technology. Declaration of Competing Interest We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome. We confirm that the manuscript has been read and approved by all named authors and that there are no other persons who satisfied the criteria for authorship but are not listed. We further confirm that the order of authors listed in the manuscript has been approved by all of us. We confirm that we have given due consideration to the protection of intellectual property associated with this work and that there are no impediments to publication, including the timing of publication, with respect to intellectual property. In so doing we confirm that we have followed the regulations of our institutions concerning intellectual property. We understand that the Corresponding Author is the sole contact for the Editorial process (including Editorial Manager and di-

16

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833 Table A1 Convergence of the Linearisation Algorithm (Detached house archetype, ECM Configurations with External Wall insulation).

δ init

δ

NPV

NPVerror [%]

TH

{1,1,1,1} {1,1,3,2} {1,1,3,3} {1,1,3,4} {1,1,4,2} {1,1,4,3} {1,1,4,4} {1,2,3,1} {1,2,3,2} {1,2,3,3} {1,2,3,4} {1,2,4,1} {1,2,4,2} {1,2,4,3} {1,2,4,4} {1,3,3,1} {1,3,3,2} {1,3,3,3} {1,3,3,4} {1,3,4,1} {1,3,4,2} {1,3,4,3} {1,3,4,4}

{1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,3,1} {1,2,3,1} {1,2,3,1} {1,2,4,1} {1,2,4,1} {1,2,4,1} {1,2,4,1}

38,718 39,201 39,227 39,161 39,237 39,262 39,197 42,584 42,402 42,418 42,320 42,521 42,355 42,379 42,262 42,841 42,703 42,716 42,615 42,846 42,703 42,746 42,636

−5.27 −4.08 −4.02 −4.18 −4.00 −3.93 −4.09 4.19 3.75 3.79 3.55 4.04 3.63 3.69 3.41 4.82 4.48 4.52 4.27 4.84 4.49 4.59 4.32

T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6

x x x x x x x x x x x x x x x x x x x x x x x

3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

Table A2 Convergence of the Linearisation Algorithm (Detached house archetype, ECM Configurations with External Wall insulation).

δ init

δ

NPV

NPVerror [%]

TH

{2,1,3,1} {2,1,3,2} {2,1,3,3} {2,1,3,4} {2,1,4,1} {2,1,4,2} {2,1,4,3} {2,1,4,4} {3,1,3,1} {3,1,3,1} {3,1,3,2} {3,1,3,3} {3,1,3,4} {3,1,4,1} {3,1,4,2} {3,1,4,3} {3,1,4,4}

{2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1} {2,1,4,1}

33,846 33,904 33,932 33,915 33,507 33,602 33,619 33,604 34,247 34,339 34,369 34,365 30,479 33,761 33,946 33,944 33,997

−17.19 −17.04 −16.97 −17.02 −18.02 −17.78 −17.74 −17.78 −16.21 −15.98 −15.91 −15.92 −20.8 −17.39 −16.94 −16.95 −16.82

T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6

x x x x x x x x x x x x x x x x x

3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

Table A3 Convergence of the Linearisation Algorithm (Mid-floor apartment archetype).

rect communications with the office). He/she is responsible for communicating with the other authors about progress, submissions of revisions and final approval of proofs. We confirm that we have provided a current, correct email address which is accessible by the Corresponding Author and which has been configured to accept email from [email protected]. CRediT authorship contribution statement Carlos Andrade-Cabrera: Conceptualization, Methodology, Software, Writing - original draft, Visualization, Investigation. Ciara O’Dwyer: Writing - review & editing, Supervision, Project administration. Donal P. Finn: Writing - review & editing, Supervision, Project administration, Funding acquisition.

δ init

δ

NPV

NPVerror [%]

TH

{1,1,1,1} {1,1,1,2} {1,1,1,3} {1,1,1,4} {1,2,1,1} {1,2,1,2} {1,2,1,3} {1,2,1,4} {1,3,1,1} {1,3,1,2} {1,3,1,3} {1,3,1,4} {2,1,1,1} {2,1,1,2} {2,1,1,3} {2,1,1,4} {3,1,1,1} {3,1,1,2} {3,1,1,3} {3,1,1,4}

{1,2,1,1} {1,2,1,1} {1,2,1,1} {1,2,1,1} {1,2,1,1} {1,2,1,1} {1,2,1,1} {1,2,1,1} {1,2,1,1} {1,2,1,1} {1,2,1,1} {1,2,1,1} {2,1,1,1} {2,1,1,1} {2,1,1,1} {2,1,1,1} {2,1,1,1} {2,1,1,1} {2,1,1,1} {2,1,1,1}

7287 7291 7298 7288 7423 7520 7545 7557 7428 7527 7554 7569 5455 5611 5581 5632 5505 5559 5579 5578

5.31 5.37 5.45 5.32 7.05 8.24 8.55 8.69 7.11 8.33 8.66 8.84 −26.49 −22.96 −23.63 −22.52 −25.34 −24.13 −23.68 −23.70

T5 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T6 T5 T6 T5 T5 T6 T6 T6

Acknowledgements This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 646116 (RealValue). In addition, the authors gratefully acknowledge the financial support of the Electricity Research Centre (ERC), University College Dublin, for this research. Appendix A. Convergence of the Linearisation algorithm

Nomenclature δinit Initial ECM guess for the linearisation trajectory β Initial ECM guess for the linearisation trajectory NPV Initial ECM guess for the linearisation trajectory NPVerror Error on the NPV estimation (percentage) T H Selected heating technology – number of devices selected

References [1] W. Adams, H. Sherali, Mixed-Integer bilinear programming problems, Math. Program. 59 (1993) 279–305. [2] AECOM, Cost Optimal Calculations and Gap Analysis for Recast EPBD for Residential Buildings, Department of the Environment, Community and Local Government, Dublin, Ireland, 2013. [3] C. Ahern, P. Griffiths, M. O’Flaherty, State of the Irish housing stock—modelling the heat losses of Ireland’s existing detached rural housing stock and estimating the benefit of thermal retrofit measures on this stock, Energy Policy 55 (2013) 139–151. [4] C. Andrade-Cabrera, W.J. Turner, D. Burke, O. Neu, D.P. Finn, Lumped parameter building model calibration using particle swarm optimization, in: Proc. 3rd IBPSA Asia Conference (ASIM 2016), Jeju, South Korea, 2016. [5] C. Andrade-Cabrera, D. Burke, W.J.N. Turner, D.P. Finn, Ensemble calibration of lumped parameter retrofit models using particle swarm optimization, Energy Build 155 (2017) 513–532. [6] C. Andrade-Cabrera, W.J.N. Turner, D.P. Finn, Augmented ensemble calibration of lumped parameter building models, Build. Simulat.. (2018). [7] M.B. Anwar, C. Andrade-Cabrera, O. Neu, M.O. Malley, D.J. Burke, An integrated building-to-grid model for evaluation of energy arbitrage value of thermal storage, in: Proc. International Conference for Students on Applied Engineering (ICSAE), Newcastle, United Kingdom, 2015. [8] C. Aravena, A. Riquelme, E. Denny, Money, comfort or environment? Priorities and determinants of energy efficiency investments in Irish households, J. Consum. Policy 39 (2) (2016) 159–186 Journal of Consumer Policy:. [9] E. Asadi, M.G. Silva, C.H. Antunes, L. Dias, A multi-objective optimization model for building retrofit strategies using TRNSYS simulations, GenOpt and MATLAB, Build. Environ. 56 (2012) 370–378.

C. Andrade-Cabrera, C. O’Dwyer and D.P. Finn / Energy & Buildings 214 (2020) 109833 [10] S.R. Asaee, A. Sharafian, O.E. Herrera, P. Blomerus, W. Mérida, Housing stock in cold-climate countries: conversion challenges for net zero emission buildings, Appl Energy 217 (2018) 88–100. [11] ASHRAE, International weather for energy calculations 2.0 (IWEC weather files) users manual and CD-ROM, Am. Soc. Heat. Refrigerat. Air-Condit. Eng. (2012). [12] G. Ault, J. Clarke, S. Gill, J. Hand, J. Kim, I. Kockar, K. Svehla, The use of simulation to optimise scheduling of domestic electric storage heating within smart grids, in: Proc. 2nd IBPSA - England Conference (BSO2014), London, UK, 2013. [13] P. Bacher, H. Madsen, Identifying suitable models for the heat dynamics of buildings, Energy Build. 43 (7) (2011) 1511–1522. [14] G.A. Bakirtzis, P.N. Biskas, V. Chatziathanasiou, Generation expansion planning by MILP considering mid-term scheduling decisions, Electr. Power Syst. Res. 86 (2012) 98–112. [15] T. Berthou, P. Stabat, R. Salvazet, D. Marchio, Development and validation of a gray box model to predict thermal behavior of occupied office buildings, Energy Build. 74 (2014) 91–100. [16] A. Byrne, G. Byrne, G. O’Donnell, A. Robinson, Case studies of cavity and external wall insulation retrofitted under the Irish home energy saving scheme: technical analysis and occupant perspectives, Energy Build. 130 (2016) 420–433. [17] G. Buttitta, W. Turner, D.P. Finn, Clustering of household occupancy profiles for archetype building models, Energy Procedia 111 (2017) 161–170. [18] CEEU, The Public Spending Code (D) Standard Analytical Procedures Guide to Economic Appraisal : Carrying Out a Cost Benefit Analysis, Central Expenditure Evaluation Unit, Dublin, Ireland, 2012. [19] CIBSE, Environmental Design - CIBSE Guide A, Chartered Institution of Building Services Engineers, London, UK, 2010. [20] CSO, Census of Population 2016, Central Statistics Office Ireland, Dublin, Ireland, 2016. [21] S. Darby, Smart electric storage heating and potential for residential demand response, Energy Efficiency 11 (2018) 67–77. [22] ESRU, ESP-r Software, Energy Systems Research Unit, University of Strathclyde, Glasgow, Scotland, 2018. [23] ECF, Roadmap 2050:A Practical Guide to a Prosperous, Low-Carbon Europe, European Climate Foundation, The Hague, Netherlands, 2010. [24] L. Evangelisti, G. Battista, C. Guattari, C. Basilicata, Lieto Vollaro, R. de, Analysis of two models for evaluating the energy performance of different buildings, Sustainability 6 (8) (2014) 5311–5321. [25] J. Ferdyn-Grygierek, K. Grygierek, Multi-Variable optimization of building thermal design using genetic algorithms, Energies 10 (10) (2017) 1–20. [26] N. Good, L. Zhang, A. Navarro-espinosa, P. Mancarella, High resolution modelling of multi-energy domestic demand profiles, Appl. Energy 137 (2015) 193–210. [27] M.M. Gouda, S. Danaher, C.P. Underwood, Building thermal model reduction using nonlinear constrained optimization, Build. Environ. 37 (12) (2002) 1255–1265. [28] Gurobi, Gurobi Optimizer Reference Manual, Gurobi Optimization LLC, 2016. [29] W.E. Hart, J.P. Watson, D.L. Woodruff, Pyomo: modeling and solving mathematical programs in python, Math. Program. Comput. 3 (3) (2011) 219–260. [30] I. Hazyuk, C. Ghiaus, D. Penhouet, Optimal temperature control of intermittently heated buildings using model predictive control: part I building modeling, Build. Environ. 51 (2012) 388–394. [31] X. He, Z. Zhang, A. Kusiak, Performance optimization of HVAC systems with computational intelligence algorithms, Energy Build 81 (2014) 371–380. [32] S. Heinen, W. Turner, L. Cradden, F. McDermott, M. O’Malley, Electrification of residential space heating considering coincidental weather events and building thermal inertia: a system-wide planning analysis, Energy 127 (2017) 136–154. [33] Hilliard, T., Swan L., Kavgic M., Qin Z., and Lingras, P. (2016). “Development of a whole building model predictive control strategy for a LEED silver community college.” [34] D. Jenkins, Integrating building modelling with future energy systems, Build. Serv. Eng. Res. Technol. 39 (2) (2018) 135–146.

17

[35] M. Jennings, D. Fisk, N. Shah, Modelling and optimization of retrofitting residential energy systems at the Urban Scale, Energy 64 (2014) 220–233. [36] A. Krener, Bilinear and nonlinear realisations of input-output maps, SIAM J. Control 13 (4) (1975) 827–834. [37] A. Kusiak, G. Xu, Modeling and optimization of HVAC systems using a dynamic neural network, Energy 42 (1) (2012) 241–250. [38] B. Lehmann, D. Gyalistras, M. Gwerder, K. Wirth, S. Carl, Intermediate complexity model for model predictive control of integrated room automation, Energy Build. 58 (2013) 250–262. [39] Linesight, Ireland Handbook 2018: Opening up a World of Data for the Irish Construction Industry, Linesight, Dublin, Ireland, 2018. [40] C.V. Loan, Computing integrals involving the matrix exponential, IEEE Trans. Automat. Contr. 23 (3) (1978) 395–404. [41] X. Lü, T. Lu, M. Viljanen, Calibrating numerical model by neural networks: a case study for the simulation of the indoor temperature of a building, Energy Proced. 75 (2015) 1366–1372. [42] P. Mancarella, C.K. Gan, G. Strbac, Evaluation of the impact of electric heat pumps and distributed CHP on LV networks, in: Proc. 2011 IEEE PES Trondheim PowerTech: The Power of Technology for a Sustainable Society, Throndheim, Norway, 2011. [43] Mathworks, Optimization Toolbox - User’s Guide, The Mathworks, 2019. [44] Mathworks, Global Optimization Toolbox - User’s Guide, The Mathworks, 2019. [45] L. Munuera, J. Bradford, N. Kelly, A. Hawkes, The role of energy efficiency in decarbonising heat via electrification, in: Proc. ECEEE 2013 Summer Study on Energy Efficiency, Stockholm, Sweden, 2013, pp. 1159–1164. [46] O. Neu, PhD thesis, University College Dublin, Dublin, Ireland, 2016. [47] O. Neu, S. Oxizidis, D. Flynn, F. Pallonetto, D. Finn, High resolution space - Time data: methodology for residential building simulation modelling, in: Proc. 13th Conference of International Building Performance Simulation Association, Chambery, France, 2013, pp. 2428–2435. [48] D. Patteeuw, K. Bruninx, A. Arteconi, E. Delarue, W. D’haeseleer, L. Helsen, Integrated modeling of active demand response with electric heating systems coupled to thermal energy storage systems, Appl. Energy 151 (2015) 306–319. [49] A. Power, Does demolition or refurbishment of old and inefficient homes help to increase our environmental, Social and Economic Viability? Energy Policy 36 (12) (2008) 4487–4501. [50] D. Pudjianto, M. Aunedi, P. Djapic, G. Strbac, Whole-Systems assessment of the value of energy storage in low-carbon electricity systems, IEEE Trans Smart Grid 5 (2) (2014) 1098–1109 IEEE. [51] T. Rasku, J. Kiviluoma, A comparison of widespread flexible residential electric heating and energy efficiency in a future Nordic power system, Energies: 12 (1) (2018) 5. [52] C.F. Reinhart, C. Cerezo Davila, Urban building energy modeling – A Review of a nascent field, Build. Environ. 97 (2016) 196–202. [53] T. Schutz, L. Schiffer, H. Harb, M. Fuchs, D. Muller, Optimal design of energy conversion units and envelopes for residential building retrofits using a comprehensive MILP model, Appl. Energy 185 (2017) 1–15. [54] SEAI, Dwelling Energy Assessment Procedure (DEAP) v 3.2.1, The Sustainable Energy Authority of Ireland, Dublin, Ireland, 2012. [55] SEAI, Conversion Factors, The Sustainable Energy Authority of Ireland, Dublin, Ireland, 2018 https://www.seai.ie/resources/seai-statistics/conversion-factors/ last accessed on 14 December 2018. [56] J.J.E. Slotine, W. Li, Applied Nonlinear Control, Pearson, Upper Saddle River, NJ, 1991. fuel tariffs. [57] SSE. (2018). SSE airtricity domestic dual https://www.sseairtricity.com/ie /home/help-centre/our-tariffs/, last accessed on 14 December 2018. [58] M. Wetter, GenOpt - A generic optimization program, in: Proc. 7th International IBPSA Conference, Rio de Janeiro, Brazil, 2001, pp. 601–608. [59] R. Wu, G. Mavromatidis, K. Orehounig, J. Carmeliet, Multiobjective optimisation of energy systems and building envelope retrofit in a residential community, Appl. Energy 190 (2017) 634–649.