Uncertainty-based optimal energy retrofit methodology for building heat electrification with enhanced energy flexibility and climate adaptability

assess and validate the effectiveness of the proposed method, under future climate scenarios and projected decreases in heating demand due to climate change. Results indicate that the best retrofit alternative of the hybrid heating system reduces carbon emissions by 88%, total costs by 54% over its lifespan, and has an average payback period of around 3 years. Air-source heat pumps can meet the majority of the heating demand (around 80%) with gas boilers used for “ top-up ” heating during high demand. Furthermore, air-source heat pumps ’ design capacity can fulfil future cooling demand even if retrofit optimization is initially focused on meeting heating needs.

To reach net zero emissions by 2050, the UK government relies heavily on heat degasification in buildings by using heat pump technology. However, existing buildings may have terminal radiators that require a higher operating temperature than what heat pumps typically provide. Increasing the size of radiators and thermally insulating building envelopes could be a potential solution, but the feasibility of these practices is uncertain due to space constraints and high retrofit costs. This study investigates the feasibility and potential benefits of incorporating air-source heat pumps into existing gas boiler heating systems to meet heating demands. The proposed probabilistic optimal air-source heat pump design method enhances energy flexibility and climate adaptability, taking into account a wide range of uncertainty sources and multiple flexibility services (e.g., energy and ancillary services). Heating systems of three educational buildings at the University of Cambridge are used as a testbed to assess and validate the effectiveness of the proposed method, under future climate scenarios and projected decreases in heating demand due to climate change. Results indicate that the best retrofit alternative of the hybrid heating system reduces carbon emissions by 88%, total costs by 54% over its lifespan, and has an average payback period of around 3 years. Air-source heat pumps can meet the majority of the heating demand (around 80%) with gas boilers used for "top-up" heating during high demand. Furthermore, air-source heat pumps' design capacity can fulfil future cooling demand even if retrofit optimization is initially focused on meeting heating needs.

Background
The UK government has committed to becoming carbon zero by 2050 [1] and has reinforced its emissions targets for 2035 of reducing greenhouse gases by at least 78% compared to 1990 [2] in line with the commitments of the Paris Agreement. To ensure success in this ambition, a series of unprecedented adjustments are required across the energy supply and demand sectors, such as increasing penetration of renewable energy sources, greater deployment of low-carbon heating and reduction of demand-side energy intensity. By exploiting freely accessible ambient heat, Air Source Heat Pumps (ASHPs) are a viable strategy to lower greenhouse gas emissions and boost energy efficiency in the building sector as the carbon content of grid electricity reduces in the future with increasing quantities of renewable electricity generation [3]. Existing buildings account for approximately 40% of the UK total demand [4], so the energy-saving potential is significant when applying ASHPs to upgrade existing heating systems (e.g., with gas boilers). Additionally, in an attempt to boost the installation of low-carbon heating technologies such as heat pumps in the building sector, the UK government has recently launched a Boiler Upgrade Scheme [5], which helps property owners overcome the upfront cost of low-carbon heating technologies.
Add-on solutions in the shape of ASHPs are necessary when lowtemperature supply water (e.g., ≤65 • C [6]) is inadequate to meet the heating demands for existing buildings with low thermal efficiency. For those buildings, the heating delivered by ASHPs might not be able to match the demand under the cold/extreme-cold weather conditions, due to the limitations of terminal heating equipment (e.g., designed for high flow temperature of 75 • C as suggested by BS-EN 442-2:2014 [7]) and the thermal comfort (mean radiant temperature). The electric demands and accompanying costs of heat pumps would significantly increase with electrification in poorly insulated buildings. Due to the low levels of building fabric thermal insulation and poor airtightness in some UK building stocks, heat pumps are forced to operate at higher flow temperatures, lowering efficiency and raising running costs. To date, there are extensive studies on the retrofitting of ASHPs in the domestic sector (e.g., Safa et al. [8] in Canada, Kelly et al. [9] in Ireland and Le et al. [10] in the UK). Prior to any heat electrification interventions, measures such as adding insulation layers should be implemented to improve the thermal performance of the building envelope and increase the overall energy efficiency of the building, hence reducing heating demands after electrification. For instance, Leibowicz et al. [11] upgraded the thermal insulation levels in a US residential building, achieving a 37% carbon cost reduction prior to ASHP installation. Lingard [12] found that solid wall insulation and low U-value glazing in UK semi-detached residential buildings are the cost-optimal solutions, which lowers 65% heat demand for effective heat pump application. However, although increasing radiator sizes and thermally insulating building envelopes could potentially address this problem, the retrofit cost could be prohibitive when fully replacing existing boilers with ASHPs [13]. In addition, the technical feasibility of increasing radiator sizes will be constrained by indoor spaces, especially for dwellings. The estimated cost of radiator upgrades and thermal insulation accounted for 33% and 60% of ASHP capital costs for domestic buildings [14], while these relative costs can be much higher due to larger fabric surface area and greater radiator capacity needed for non-domestic buildings. Alternatively, adopting hybrid heating systems, which compose of ASHPs and gas boilers (GBs), can be a promising approach where replacing existing radiators or updating envelopes is economically unfeasible.
Hybrid operation and coordination of electric-driven ASHPs and natural gas-driven GBs also have enormous potential to provide energy flexibility for smart grids. Recently, new policies have been passed to encourage demand to contribute to the power balance of the power grid and to participate in ancillary services, also known as "demand response" [15]. Heat pumps, as large electricity consumers in buildings, can be engaged in energy and ancillary services markets (incentivebased and price-based programs [16]). For instance, heat pumps can effectively reduce grid frequency deviations and earn revenue from electric grids by regulating the variable speed drive in response to grid control signals [17,18]. Meanwhile, model predictive control (MPC)based control strategies enabled heat pump response to dynamic realtime electricity pricing, yielded 5% energy reduction and 4.5% savings in [19], and 10.8% of energy cost savings in [20]. In addition, heat pumps can be coupled with advanced building technologies, such as thermal energy storage [21], electrical energy storage [22] and photovoltaic systems [23], to further improve energy flexibility and economic benefits of building operations. However, very limited studies take into account the energy flexibility advantage of ASHPs in hybrid heating systems (i.e., composed of GBs) when performing retrofitting analysis and design optimisation.
Another challenge encountered in practice is a lack of confidence in the energy-saving potential of such systems. Even well-selected data/ information (e.g., building fabric, occupancy, weather, equipment load, etc.) used as inputs for retrofit analysis can be rather different from real operation conditions, due to inherent deviations in climate change, building demands, energy policies and equipment performance. Such a difference will cause the performance of the system to be deviated from the expected, possibly accompanied by increased cost and energy consumption. For instance, the average Earth's surface temperature has increased by at least 1.1 • C since 1880 [24], with most of that warming occurring after 1975 at a rate of roughly 0.15 to 0.20 • C per decade [25]. The impacts of climate change affect both the energy efficiency of ASHPs and the future building energy demands [26]. Other intrinsic uncertainties, such as fuel price fluctuations and equipment performance degradation, might hamper the potential benefits of a retrofit system. A slight change in the values of these uncertainty parameters might cause a mismatch between ASHPs' capacity and their operation/control, thus an unexpected effect on the return on investment for heating system retrofit [27]. Therefore, it is necessary to size the heating systems such that the variation in performance, based on the variation of design variables and noise factors, becomes minimal. In addition, taking different uncertainties into account helps achieve a cost-effective design option which provides the system with the capability to operate at high efficiency under all possible conditions, particularly at partial load, through having the ability to accommodate design and operation uncertainties.

Problem definition, research objectives and original contributions
Based on the literature review, the following observations were made: (1) Hybrid heating systems, by incorporating ASHPs into existing GBs, could be an economical approach for existing building retrofit. (2) The bivalent heating system (ASHPs + GBs) provides significant energy flexibility potential that should be taken into account in the retrofit analysis. (3) The impacts of multiple uncertainty sources in future scenarios (e.g., weather, heating demand, energy prices), on the actual performance of the hybrid heating system should be part of the system design and optimisation process. The following research questions (RQs) were, thus, formed: • RQ1: How to quantify and integrate uncertainty sources (such as weather and building heating demands) that affect ASHP system performance into existing heating system retrofit? • RQ2: How to properly size the ASHPs while coordinating the operation of the hybrid ASHP-gas boiler systems to provide enhanced energy flexibility? • RQ3: What is the lifespan performance (e.g., cost and carbon emission reduction) of the hybrid heating system? • RQ4: Will the retrofitted ASHP systems also satisfy future cooling demand in light of global warming?
To address the above-mentioned research gaps, this work aims to develop an integrated uncertainty-based optimal energy retrofit design method for building heat electrification, taking into account energy flexibility provision and the multiple uncertainty sources. An automated global climate model and a self-correcting building energy model generate uncertain weather and heating load profiles considering the impacts of climate change. A hierarchical optimization framework is developed which determines the baseline power of ASHPs in energy service markets and selects the optimal mode and service type in ancillary service markets, respectively. We illustrate this novel methodology using an existing gas boiler heating system serving three educational buildings located at the University of Cambridge, UK. Achievable savings under the proposed methodology are tested by comparing the lifespan performance of the retrofitted hybrid heating system to that of the original heating system.
The principal difference between existing scenario-based studies and ours is that we propose a generalisable methodology for the design and operation of heating system retrofits. Additionally, retrofit designs that incorporate potential demand response capabilities are often overlooked. The original contributions of this study are briefly summarised as follows: (1) It represents the very first effort to optimise the design of bivalent heating systems (consisting of new ASHPs and existing GBs) during the retrofit stage to maximise potential benefits while considering a wide range of uncertainty sources, with a self-correcting building model and automated global climate model. (2) A novel hierarchical optimization framework is proposed that provides multiple flexibility services for retrofitted heating systems in sequential electricity markets. Fig. 1a shows the schematic of the original gas boiler heating system. This system adopts constant-speed water pumps. The return water from the building is supplied by circulation pumps to the GBs. The hightemperature hot water (HTHW) heated by the GBs is then supplied to buildings. The supply water temperature is determined according to building heating loads with the design water temperature (e.g., 75 • C). Fig. 1b shows the schematic of the retrofitted hybrid heating system, including the existing GBs and new ASHPs. The gas boilers serve as a backup heating source during periods of high demand, while the ASHP has a maximum low-temperature hot water (LTHW) (e.g., 65 • C) due to various limitations such as system efficiency, component stress, and available heat from the outdoor air. If the required return water temperature from the building exceeds the maximum LTHW temperature set by the ASHP, the gas boilers will take over all heating demand. Three operation modes are adopted using the hybrid heating system while the selection of the best operation mode is based on the estimated operation cost and the availability to meet the heating loads. i. GB mode: only the GBs will be used for providing heating services and the electric valves at the ASHP side are fully closed. ii. hybrid operation mode: the return water from the building has two parallel flows, one is heated by ASHPs, and the other is heated by GBs. The outlet water from the GBs and ASHPs is mixed and eventually supplied to buildings. iii. ASHP mode: only the ASHPs will be used for providing heating services and the electric valves at the GB side are fully closed.

Outline of uncertainty-based energy retrofit optimisation
It should be noted that hybrid heating systems can be configured in series or parallel, but here parallel configuration is considered advantageous for several reasons. Firstly, it has a lower system pressure drop, which can reduce the energy costs associated with circulation pumps. Additionally, the ASHP is used as the primary heating source, with the GB serving as a backup during periods of high demand. The parallel configuration is easier to control because we can directly close the GB water loops most of the time. In contrast, the series configuration would require a bypass pipe, which increases control complexity. The heating system performance of the retrofit scenario is compared with the reference scenario. The reference scenario utilises existing gas boilers for heating only, while the retrofit scenario re-uses existing gas boilers and introduces new ASHPs for hybrid heating. The modification of building fabrics and radiator systems was not taken into account to save investment costs and prevent construction in working spaces. The key issue of the hybrid heating system involves how to properly size the ASHPs to minimise the life-cycle costs, as an add-on solution for existing gas boiler heating systems.
To optimize the sizing of ASHPs with the best operation modes, an uncertainty-based energy retrofit optimization method is proposed as shown in Fig. 2. It contains two steps: uncertainty quantification and retrofit optimisation. In the first step, the weather profiles, building heating demand, component performance deviations of hybrid heating systems and economic factors (energy prices, inflation rate, discount rate) for future scenarios are estimated. It is important to recognise that the key parameters affecting building heating demand differ between the design stage and the retrofit stage. For example, while building characteristics at the design stage are typically estimated based on assumptions, during retrofitting, these properties can be obtained through on-site investigations or architectural drawings, leading to a higher degree of accuracy. In the second step, all future uncertain parameters are used as the inputs for the hybrid energy system models, and the "optimizer" determines the number and capacity of ASHPs, thus, minimizing the lifecycle cost of retrofitted hybrid heating systems. The rated heating capacity for each ASHP is assumed identical. The main function of the energy system model is to determine the best operation mode according to the operation cost, taking into account the electricity cost, gas cost, flexibility revenue and carbon cost.

Energy flexibility available by hybrid heating systems
The hybrid heating systems proposed here account for the provision of energy flexibility in electricity markets of multiple services. The hybrid heating systems, driven by electricity (ASHPs) and natural gas (GBs), strategically optimise peak shaving and allocate the available flexibility capacities in the sequential electricity market (run by market operators and system operators respectively in Europe). Market operators, such as Nord Pool [28], provide platforms for trading energy products (in day-ahead and intraday markets) before the hour of operation. System operators, such as the European Network of Transmission System Operators for Electricity (ENTSO-E [29]) provide platforms for trading ancillary services products (in balancing markets) in the hour of operation. The energy product (kWh) is marketed for end-user energy consumption, whereas the ancillary service products (kW) provide revenues to end-users. According to Zhang et al. [30], the hybrid heating systems can provide the flexibility services such as energy service (e.g., load covering [31]) and ancillary service (e.g., Frequency Containment Reserve (FCR) [32] and Frequency Restoration Reserve (FRR) [33]), as shown in Fig. 3. Fig. 1. Schematic of the (a) gas boiler heating system (b) hybrid heating system.

Energy service-load covering flexibility
As a common energy service, "load covering" reduces a building's electric loads during peak electricity times. Peak and off-peak times refer to the price signal (e.g., time-of-use, critical peak pricing, and hourly real-time pricing) at high and low prices respectively that reflect power grids' daily flexibility requirements. This flexibility type reshapes a building's daily electric demand profile by intelligently controlling deferrable loads (e.g., using GBs for heating if the natural gas price is lower) in the context of cost-optimal control or other optimization goals. Notably, the covering flexibility provided by the hybrid heating system is the electricity load change of the ASHP system throughout load reduction periods. The load covering capacity provided by ASHPs is the function of the covered building heating load and overall coefficient of performance (COP) of ASHP systems, as shown in Eq. (1). The GBs can provide standby heating when the electricity price is high and natural gas price is low as presented by Eq. (2).
where E cv ASHP = load covering capacity provided by ASHPs (kWh). P load,cv = covered building heating load (kW) at hour t. Δτ cv = covering periods (h).
P GB = boiler heat output (kW). "Load covering" flexibility can be provided by the hybrid heating system in both day-ahead and intraday markets. In the case study, both the day-ahead market and intraday market were selected to provide energy services.

Ancillary service-FCR
Frequency Containment Reserve (FCR), also known as "primary reserve", aims to automatically stabilise the frequency after the occurrence of minor and unanticipated imbalances. This service's actions begin no later than 30 s after the imbalance, with a reaction period of up to 15 min. Existing studies have shown the potential of heat pumps in providing such a service by controlling the variable speed drivers [34]. The flexible loads of buildings qualified for such service have the potential to provide FCR capacity (e.g., fastest responsive flexibility) in the ancillary service market as shown in Eq. (3).

Ancillary service-FRR
Frequency Restoration Reserve (FRR) aims to address imbalances that are too severe or prolonged for the FCR to address, including automatic frequency restoration reserve (aFRR, also known as "secondary reserve") and manual frequency restoration reserve (mFRR, also known as "tertiary reserve"). From the frequency variation, it operates from 30 s (aFRR) to maximum 15 min (mFRR). The hybrid heating system can provide either aFRR or mFRR services. According to the European Network of Transmission System Operators for Electricity (ENTSO-E), the average regulation price for aFRR in 2021 was £10.4/ MW, while the average regulation price for mFRR was £6.8/MW [35]. Upon the situation that the aFRR service has higher revenue than mFRR, we consider the aFRR service to provide ancillary services. The aFRR capacity of the ASHPs is defined as the baseline power output of ASHPs for the shedding duration, in response to a contingency event of the smart grid, shown in Eq. (4).

Objective function
The choice of an appropriate objective function, which simultaneously considers investment cost, energy cost and carbon cost, is a complicated decision. The formulation of the optimization objective for retrofitting existing building heating systems with add-on ASHPs is presented as Eq. (5), aiming at minimising system life-cycle costs (C total ).
where C inv is the investment cost, including the capital and installation costs of ASHPs (Eq. (6)). The installation cost is assumed 30% of the capital cost [36]. The capital costs of ASHPs are determined by the selected manufacturer and size at the beginning of the retrofit and are therefore deterministic. As no changes would be made to the current GBs, it is assumed that there will be no investment cost for GBs. C egy stands for the energy cost, including the net electricity trading costs and natural gas costs (Eq. (7)). C cbn is the carbon costs charged due to the greenhouse gas emissions raised by UK Emissions Trading Scheme (Eq. (8)) [37]. C mt is the additional maintenance cost due to the installation of ASHP. The annual maintenance cost is set as 2.5% of the capital cost [38].
C egy = C ele + C gas (7) where γ ASHP = the unit price of an ASHP (GBP/kWh). n = total number of ASHPs. γ cbn = carbon price (GBP/tCO 2 e). ξ gas = natural gas carbon emission factor, set as constant 0.184 kg/ kWh. ξ ele = electricity carbon emission factor (kg/kWh). E gas = purchased energy from natural gas at the period T (kWh). P ele (t) = ASHP power at the hour t (kW). The net electricity trading costs by ASHPs are presented in Eq. (9), where the trading involves energy and ancillary service products. The load covering, FCR, and aFRR flexibility (i.e., the amount of purchased energy and offered ancillary service capacity) of ASHPs are fully utilised, subject to the constraints of hierarchical operational flow as further elaborated in Section 2.5. The hybrid heating system can only provide one type of ancillary service to meet regulation quality requirements. The total provided load covering capacity and FCR capacity (or FRR capacity) at each time step should not be higher than the rated capacity of ASHPs, while the total reserve capacity should not exceed the baseline power of ASHPs. According to the proposed assessment criteria of building flexibility, the hourly regulation capacity is anticipated to be symmetrical in both up and down directions. The impacts of FCR flexibility on energy consumption are ideally neutral. The natural gas (fuel) cost by GBs is presented in Eq. (10). The unit price of ASHPs is the function of the corresponding capacity expressed by Eq. (11), curve fitting derived from Daikin UK Price Book [39], including an ASHP and additional accessories (such as low-lift pumps, energy meters, inverter, etc.).

Performance models of hybrid heating system
Radiator: The heat output of the radiators can be calculated by Eqs. (12) and (13).
Heating network: As heating network distribution loss is unavoidable, heating network efficiency should be considered. To simplify the calculation, the distribution loss of the supply water pipe is considered to be equal to that of the return water pipe. The system heating load, supply and return water temperature can be calculated by Eqs. (14)- (16), respectively.
where Q sys = system heating load (kW). t i,sys = system supply water temperature ( • C). t o,sys = system return water temperature ( • C).
Air-source heat pump: The ASHP power is the function of the ASHP heating load and operating COP (Eq. (17)). The operating COP of ASHP is calculated by multiplying the full load COP by a part load factor (PLF) as Eq. (18). The COP at maximum power was determined as a function of the water temperature exiting the ASHP and the ambient temperature using biquadratic polynomial curve fitting generated from operational data [8], as presented by Eq. (19). PLF is calculated as a function of the part load ratio (PLR) in Eq. (20).
where Q ASHP = ASHP heating load equals to Q sys when heating is provided by ASHPs only (kW). Gas boiler: The total energy provided by natural gas (kWh) can be calculated by Eq. (21). Gas boiler efficiency is a function of the PLR presented by Eq. (22), obtained from a typical part-load efficiency curve for a non-condensing boiler [41]. (21) η GB = 0.7107PLR 3 − 1.4216PLR 2 + 0.8825PLR + 0.6693 (22) where Q GB = gas boiler heating load, which equals to Q sys when heating is provided by GBs only (kW). η GB = gas boiler efficiency.

Probabilistic future weather profiles
In this study, we propose a probabilistic future weather model to generate the possible future weather profiles instead of using existing morphed climate change weather files for the following two reasons. i. Unlike the existing weather files, such as from UK Climate Projections 2018 (UKCP18) based on CMIP5 (Coupled Model Intercomparison Project 5), the proposed method utilises the global climate model (GCM) data from the CMIP6, which develops the most up-to-date climate model experiments with finer resolution and improved dynamical processes [42]. ii. The proposed probabilistic future weather model takes advantage of existing GCMs simultaneously and therefore can generate more accurate weather data and take account of more possibilities of future weather scenarios.
The process of generating probabilistic future weather profiles is shown in Fig. 2 (upper left corner), which involves three steps. The first step is the selection of base Typical Meteorological Year (TMY) weather files, representing the uncertainty in present-day weather [43]. The base TMY files used in the morphing process were two weather files from the University of Exeter (Ext1, Ext2) [44], two Meteonorm weather files (Metn1, Metn2) [45], and an Ensim weather file (Ensim) [46]. In the second step, the morphing method proposed by Belcher et al. [47] was introduced to generate future weather profiles. This method generates future weather data by merging current-day observed weather data (referred to as TMY in the first stage) with outputs from time seriesbased global climate models (GCMs). To represent the uncertainty in climate models, 67 GCMs from the ScenarioMIP and HighResMIP Endorsed Model Intercomparison Projects were utilised to generate future climate data. ScenarioMIP simulated four shared socioeconomic pathways (SSPs): SSP126, SSP245, SSP370, and SSP585 [48], whereas HighResMIP simulated three experiments: control-1950, highres-future, and highresSST-future [49]. Electronic Supplementary Material (ESM) Table A-1 lists the detailed information on these 67 climate models. In the third step, each GCM model's forecasted historical monthly weather data between 2015 and 2021 were compared to actual measured data from a local weather station. The performance of GCMs was ranked from the lowest to highest according to the Root Mean Square Error (RMSE) of monthly average temperature. The GCMs were then selected given these GCMs had the highest Prediction Interval Coverage Probability (PICP, Eq.24). Here, PICP is defined as the proportion of observed monthly average temperature located at the predicted monthly average temperature confidence interval. The larger the PICP is, the better the model performs. In the fourth step, future hourly weather data from 2022 to 2041 were generated for retrofit optimization by applying the morphing approach to the output of the selected GCMs, assuming that the lifespan of ASHPs is 20 years [50].
where r i = measured data points for each model instance 'i' q i= simulated data points for each model instance 'i' N p = the number of data points at interval 'p' l i = lower limit for each model instance 'i' l i = upper limit for each model instance 'i'

Probabilistic building heating demand
The probabilistic weather profiles were further used as the inputs for building simulation models, which in turn generated future building demand. TRNSYS3D, a SketchUp plugin that enables users to draw multizone buildings and import geometry (including building selfshading and internal view factors for radiation exchange) directly from the SketchUp interface into the TRNSYS Building environment (TRNBuild, TYPE 56 [51]), was utilized. The building model in SketchUp features the same geometry and envelop details as the case buildings. Building's internal gains, such as lighting, equipment, and occupancy were assumed not to be changed in the future years, and thus only climate change impacts on building heating load were investigated. Based on inputs (x 1 , x 2 , …, x n ), the outputs Q (i.e., building heating loads) are obtained as Eq. (25). The design inputs involving uncertainties (X) are generated by Latin hypercube sampling (LHS) method as Eq. (27), according to their probability distributions (G).
A self-correcting building model was first developed to generate the building heating load as closely as feasible to the real heating load. The automated calibration process compares a sequence of differences between hourly heating load profiles generated from building simulations versus actual measured heating load by smart heat meters, as shown in Fig. 4. One goal in developing the auto-calibrated approach was to build a simple and time-efficient calibrated model with minimal inputs. The measured hourly heating load, meteorological data, and the to-becalibrated model were the only inputs necessary for the building model calibration process. The method starts with the input data and relies on parameter auto-tuning to eliminate inconsistencies between the simulated and actual load profiles iteratively. The parameter and parameter range for tuning were chosen based on the specific characteristics of the load profile mismatch. This iterative parameter tuning process continues until the load profile converges within a given tolerance.
A parameter selection process was developed using sensitivity analysis and engineering judgement to list all available parameters for potential selection during auto-tuning, as shown in the co-author's previous study [52]. By tuning different input parameters, the priority list of parameter selection was established. The parameters that exhibited the widest uncertainty following minor modifications were selected. The fabric parameters of the building, such as window U-value and window SHGC, are not tuned as they can be obtained from the design drawings and won't be significantly changed over building operation time. The building's internal gains, heating setpoint, infiltration and ventilation, which are of high randomness and are changed over time, are selected for calibration.
For each automated tuning stage, a maximum adjustment of 5% was made to the parameters based on the last calibrated value. The tuning direction of the selected parameter was assigned a sequence of values ranging from the initial value to the upper or lower limit. A series of submodels were generated in this process, while only the sub-model with the lowest Mean Bias Error (MBE) and Coefficient of Variation of the Root Mean Square Error (CVRMSE) was chosen for the next calibration stage. The MBE and CVRMSE are calculated in accordance with ASHRAE Guideline 14 [53], where r is the average of the measured data points. The values of 10% for the MBE and 30% for the CVRMSE were used as the threshold level.

Uncertain economic and performance factors
Probabilistic energy and carbon prices: Natural gas and electricity prices in the baseline year (i.e., 2021) were multiplied by the relative prices to determine the future prices (Eqs. (30) and (31)). These energy price relatives followed a uniform distribution (Eq. (32)), with the lower and upper bounds obtained from US Energy Information Administration (EIA) [54]. The projections take different policy scenarios (e.g., oil price, economic growth, and renewable costs) into account. The future carbon prices and electricity carbon emission factors were also set following a uniform distribution (Eqs. (33)-(34)), with the lower and upper bounds obtained from the UK Green Book supplementary guidance [55].  Factor (PVF) to derive the present value (Eq. (35)), where r is the interest rate and d is the discount rate. The interest rate and discount rate for each year followed stochastic distributions based on the statistics data from the UK government [56].
Equipment performance variation: The performance of GBs/ASHPs may decline over time due to the natural ageing of the equipment, characterized by a decrease in heating output over an extended period of operation, resulting in diminished energy efficiency and increased system lifecycle expenses. In the meantime, regular maintenance can help alleviate the performance degradation of the equipment (represented by the ageing alleviation rate). The rated coefficient of performance (COP) of the component in year J was estimated by Eq. (36) considering both the effects of natural ageing and regular maintenance, where COP rated,0 is the equipment rated COP when newly installed or calibrated. The   Fig. 5. A hierarchical analytical process for optimal retrofit optimisation. uncertainty of the annual degradation rate (D a ) and ageing reduction rate (R a ) can be characterised as stochastic distributions [57,58]. Fig. 5 presents a hierarchical analytical process for optimal retrofit optimisation based on the probabilistic input parameters outlined in Section 2.4. The proposed hybrid heating systems offer three different operation modes (i.e., running HPs only, running GBs only, and hybrid running of HPs and GBs) that enhance energy flexibility. During each time step of performance evaluation, the feasibility of utilizing optimal operation modes was assessed by verifying whether the required heating capacity and thermal constraints (e.g., the supply water temperature of ASHP <65 • C) were met for the given trial ASHP number and capacity provided by the optimiser.

A hierarchical analytical process of optimal retrofit optimisation
By following the bidding sequence, in stage 1, the baseline power in the energy services markets (i.e., day-ahead market and intraday market) was determined for the hybrid heating system without providing ancillary services. The operational schedule at each operation hour was established by comparing the power among available operation modes (e.g., ASHP mode, GB mode, or hybrid mode) and selecting the optimal mode with minimal power to meet the space heating demand. In stage 2, the service type for each available operation mode at each operation hour was selected, and subsequently, the optimal operation mode with minimal operation cost was chosen. It is important to note that ASHPs cannot provide both FCR and aFRR simultaneously as these services require different types of response and may have different priorities depending on the specific needs of the electricity grid. Subsequently, a golden section search and parabolic interpolation algorithm [59] was used to search for the optimal retrofit alternatives, which ensured that the overall lifecycle cost of the hybrid heating system reached the minimum within the pre-set convergence tolerance.

The reference buildings and their heating systems
Three GBs with a total heating capacity of 3060 kWh, serve three educational buildings through a shared heating network, as shown in Fig. 6. The three buildings are located at the University of Cambridge with a total area of 34,560 m 2 . The terminal radiators were designed with 75/65 • C (water flow/return temperature) and heat meters were used to measure buildings' total heating load on an hourly basis. Table 1 illustrates the four distinct categories considered for retrofit optimization, which include weather, building heating demand, economic factors, and performance factors. Economic, emission and   *Note: Γ(k, θ), N(µ, σ 2 ) and U(a,b) represent gamma, normal and uniform distributions, respectively. µ refers to the mean value and σ refers to the standard deviation. k and θ refer to the shape parameter and scale parameter respectively. a is the lower limit and b is the upper limit. The distributions for the degradation rate of ASHPs and GBs per year were derived from the technical report of the National Renewable Energy Lab (NREL) [57]. The distributions for the ageing alleviation rate of ASHPs and GBs were derived from the U.S. Department of Energy [58], taking into account performance enhancements achieved through routine maintenance. performance factors are introduced by random sampling according to their distributions. The future weather profiles were generated from probabilistic GCMs, while the future building heating demands were derived from calibrated building models that utilised the projected weather data as inputs. Details can be found in Sections 3.2.1-3.2.3.

Future weather profiles
Base TMY files and GCMs: A comparison between the weather parameters (dry bulb temperature, relative humidity and global horizontal irradiance) of different TMY files can be found in Fig. 7. The external mean dry bulb temperature, relative humidity and global horizontal irradiance were 9.6-12.0 • C, 74-81% and 199.5-224.7 W/m 2 , respectively. The variances of the weather data came from the historical time slices (Met1: 1996-2019; Met2: 1961-90; Ext1&Ext2: 1961-90; Ensim: 1986-2021) and weather station locations (Met1, Met2&Ensim: Cambridge airport; Ex1&Ext2: Cambridge botanic garden). The 5 base TMY files were used as the inputs for 67 climate models (see ESM Table A-1). Therefore, a total of 5 × 67 = 335 future weather scenarios were generated.
Selection of representative GCMs: To select the representative GCMs that can generate the future climates, the RMSE of monthly mean temperature between GCMs and actual measurement was first calculated. The GCMs were then sorted according to the values of computed RMSE (i.e., from the lowest to the highest). Fig. 8 shows the RMSE and PICP of monthly mean temperature with the increased number of GCMs. With the increase in the cumulative share of GCMs, the RMSE first sharply decreased, which bottomed at 1.21, and then increased. The best individual GCM achieved the minimum RMSE at 1.57, while multiple GCMs performed better than the best individual GCM. The PICP of monthly mean temperature increased with the increased share of GCMs and then stabilised at 98.8% when the share of GCMs was over 33%. This indicates that 33% of GCMs would be enough to represent all the features of all GCMs. As shown in Fig. 9, 110 weather scenarios were selected and compared with the actual measured outdoor air temperature between 2015 and 2021. Almost all actual monthly mean temperature was in the temperature ranges of probabilistic GCMs, and thus the accuracy was satisfactory.
Generation of future climates: Future 20-year (from 2022 to 2041) weather profiles generated from probabilistic GCMs were presented in Fig. 10. The annual mean air temperature climbed by 0.7 • C from approximately 10.9 • C in 2022 to 11.6 • C in 2041 due to climate change. The annual mean relative humidity ranged from 76.9% to 78.3%, while the annual mean solar radiation varied between 1138.3 kWh/m 2 and 1179.1 kWh/m 2 . Future climate scenarios with increased average outdoor air temperatures are expected to improve COP for ASHPs used in heating applications.

Future building heating load
Calibration of building models: Table 2 illustrates the detailed calibration process of the building models. Following a total of 186 tuning steps for 5 parameters, the building model was successfully calibrated, with the heating load NMBE at 4.0% and CVRMSE at 25.9%. The 5 uncertain parameters, including heating setpoint, infiltration rate, appliance power density, lighting power density and occupant heat gain were selected for the calibration process. The sequence of the selected 5 parameters was predetermined by sensitivity analysis performed on an office building in Cambridge [52]. The calibration results can be found in Fig. 11. The R 2 between the predicted heating load and the actual heating load increased from 0.62 (Step 0) to 0.98 (Step 5), indicating the superior performance of calibrated building model.
Generation of future building heating load: Based on the calibrated building model, the future building heating load was generated by using the probabilistic future weather as inputs. It is assumed that the building heating load was dominantly impacted by the future weather, while the other parameters of the buildings, such as internal gains and equipment efficiencies were kept unchanged following the settings of the calibrated building model. Fig. 12 shows the future building heating load from 2022 to 2041. The annual heating load ranged between 47.

Optimal sizing of ASHPs
The hybrid heating systems consist of the original gas boilers (a total heating capacity of 3060 kWh) and the newly add-on ASHPs. Fig. 13 presents the required ASHP capacity at different numbers and the lifecycle costs of the hybrid heating systems. It can be found that the optimal number of ASHP was 3 with each capacity of 468 kW. The capital cost per capacity decreased with the increased capacity per ASHP according to Eq. (11). Compared with the utilisation of GBs, although increasing total capacity (when the ASHP number is over 3) can prolong the ASHP running time, the increased capital cost would exceed the reduced operation cost. The optimal solution can reduce the total mean lifecycle cost by 52% compared to that of the original GB system.    (existing GBs) and the hybrid system (existing GBs + new add-on ASHPs). From 2022 to 2041, the retrofitted hybrid system reduced annual mean operation costs by 54% to 70%. This attributes to the following reasons: (1) ASHPs offered improved energy efficiency to reduce cost and participate in power flexibility services to earn revenue.
(2) Natural gas prices were projected to be increased while electricity prices were projected to be decreased due to the increased penetration of cheaper renewable energy as presented in ESM Fig. B-2. Fig. 15 shows the mean operation cost breakdown of the original heating system and the hybrid heating system in the next 20 years, where the operation cost includes the energy cost (+), carbon cost (+) and revenue (− ). Compared with the original GB system, the hybrid heating system achieved 62% lifespan operation cost savings. Given the capital cost of 431 k GBP of the ASHPs, 54% lifespan total cost saving can be achieved. It is worth noticing that the hybrid heating system also obtained revenue from the smart grids, which accounted for 1.1% of total cost savings.
To further elaborate on the energy flexibility provided by the hybrid system, we analysed the running ratio (frequency) of the hybrid system, shown in Fig. 16. The HP mode was dominated over 20 years with an average running time of around 82%, indicating that operating ASHPs without gas boilers can meet the majority of the heating demand through the lifecycle. Hybrid mode occupied 17% of running time, while GB mode was the least frequent operation mode (less than 1%). This indicates significant carbon reduction potentials due to the reduced usage of natural gas, as further elaborated in Section 4.2.2.

Carbon reduction
The retrofit hybrid system had the potential to reduce carbon emissions. Fig. 17a presents the carbon emissions between the original system and the hybrid system. Due to the increased heating load as well as the performance degradation from 2022 to 2041, the carbon emissions of the original system would increase from 525 tons to 536 tons. In contrast, the electrification of the building energy system by utilising the  Fig. 11. Comparison between the predicted heating load and actual heating load. ASHPs would reduce carbon emissions from 106 tons to 28 tons due to the decrease of electricity carbon emission factors in the future 20 years (ESM Fig. B-4). The carbon reduction rate by the hybrid system was from an average of 80% in 2022 to 95% in 2041, as shown in Fig. 17b. The hybrid air-source heat pump/gas boiler system offered an average of 88% lifespan carbon emission reduction.

Discounted payback period
The payback period indicates the number of years needed for the payback of the surplus capital cost for installing ASHPs to the existing heating systems. The hybrid system had a satisfactory economic performance. It can be seen from Fig. 18 that the discounted payback period of the hybrid system varied between 2.5 and 3.4 years with a mean value of 2.9 years. The shorter payback period indicates the better economic performance of the hybrid system. For the hybrid heating system, the ASHPs mainly provided heating under the part load conditions, while GBs provided heating under near full or full load conditions. Compared with the retrofit option that fully replaces the GBs with ASHPs, the hybrid operation of GBs and ASHPs could greatly reduce the capital cost as the required capacity of ASHPs would be reduced. It is worth noticing that the payback period for the full replacement of GBs with ASHPs is not discussed in this study. The reason is that when fully replacing of GBs with ASHPs, additional retrofit measures such as an upgrade of radiator and fabric should be performed. Otherwise, the ASHPs will fail to meet the building heating requirements due to a lower supply water temperature. Our future studies will investigate these holistically taking into account radiator and fabric upgrade costs.   Fig. 19 shows the building cooling loads generated from the calibrated building model considering the climate impacts. The annual mean cooling loads were projected to increase from 9.5 kWh/m 2 in 2022 to 13.6 kWh/m 2 in 2041 obtained from the calibrated building model.

Potentials to meet future cooling demand
Would the retrofitted ASHPs have the potential to meet the future cooling demand? To answer this question, the total cooling capacity of the retrofitted ASHPs over the operating year was compared with the maximum cooling load. Here, the design cooling capacity of the ASHPs was obtained by referring to the technical manual [39]. According to this technical manual, if an ASHP has a design heating capacity of 468 kW, its cooling capacity is 854 kW. If the equipment performance degradation is taken into consideration, the actual cooling capacity of the ASHPs will be reduced throughout their service lifecycle. Fig. 20 shows the comparison of the building's maximum cooling load and ASHP's total cooling capacity over the future 20 years. It can be found that the ASHP total cooling capacity would be higher than the maximum cooling load in most scenarios, indicating that the retrofitted ASHP can mostly meet the future cooling demand. This is rather important because, in the context of climate change, ASHPS can effectively address the intensifying concerns over the mitigation of indoor overheating in more insulated and airtight buildings (with higher thermal efficiency) in the absence of other appropriate cooling means [60].

Research limitations and future work
The main focus of this study is to introduce a retrofit strategy for buildings that rely on gas boilers for heating, a practice commonly found in the UK, North America, and Southern Europe. The proposed method may not be directly applicable to central heating systems typically seen in Scandinavia, Northern Europe, Eastern Europe, and Northern China. To retrofit central heating systems, optimisation of the heat distribution network and integration of large-scale renewable energy sources would be necessary, necessitating increased collaboration among stakeholders and taking into account the wider energy strategy and infrastructure within the network's service area. Based on the performance tests of the developed method and a review of the existing literature, several issues still need to be addressed in further studies to ensure successful implementation in practice.
• This study assumes that building internal gains would not be changed, and thus only investigate the impacts of climate change on building heating load in future scenarios. However, there might be some uncertainty based on the following current trends (some of them possibly counteracting each other): increasing use of appliances, especially IT equipment, increasing equipment energy efficiency, and reduced occupants due to remote working. For instance, a study by Jenkins et al [61] shows ASHPs could have an improved COP with reduced required size in the "2030" office, compared with the baseline "2005" office. The impact of internal gains on building heating load and system sizing should be further explored. • Identical ASHPs with different numbers were designed in this study to simplify control implementation, installation, and maintenance. However, it may not be the most efficient approach to account for partial load in some applications [62]. Therefore, a case-dependent investigation is necessary to determine the optimal sizing for improved efficiency, including different combinations of ASHPs with varying numbers and capacities. • The proposed uncertainty-based optimal energy retrofit method for integrating ASHPs and GBs was conducted under the condition that ideal controls were adopted to achieve the intended operation of a specific hybrid operation. Although this is reasonable under the retrofit planning stage, several issues still need to be addressed in further studies for successful implementation in practice. i. Developing an online control strategy, e.g., model predictive control [63], can improve environmental control and limit comfort sacrifice. ii.
Online adaptive and calibration strategies [64] are also useful to refine the operating schedule of ASHPs and GBs when considering uncertainties and performance degradations. iii. The hybrid heating system can be further integrated with smart energy storage systems to increase system energy flexibility and participant in multiple demand response programs in electricity markets. • Although plant rooms are normally designed to allow for a minimum of 20% additional space for equipment expansion [65], the additional spaces required for accommodating both the GBs and ASHPs of hybrid heating systems should be further investigated to test the feasibility and applicability prior to the energy retrofits. • The case study of Cambridge buildings has proved the generalizability of the proposed methods, by adopting the self-correction building model and automated global climate models. Nonetheless, the proposed method could be adapted and applied to other hybrid configurations of heating systems, depending on the availability of data and resources. • The consideration of the update of envelope systems has been ignored in this study, as the main focus is on developing a methodology for retrofitting heating systems. The importance of retrofitting both the building envelope and heating systems to minimize heating demand and reduce carbon emissions in some scenarios is recognized. In future work, it is possible to extend the methodology to include fabric upgrades within the methodology. • This study focuses on exploring the benefits of integrating ASHPs into existing GBs without factoring in other energy systems associated with higher capital costs like thermal energy storage (TES) and renewable, waste heat, biomass, and geothermal energy systems.
However, exploring the potential of upgrading building envelopes, enhancing building thermal inertia, and incorporating TES and electric boilers in future research can lead to further improvements in the energy flexibility of building heating systems.

Conclusion
This study proposes an uncertainty-based optimal energy retrofit  method for integrating the air sources heat pumps (ASHPs) with the existing gas boiler (GBs) heating systems, providing enhanced energy flexibility and climate adaptability. Four categories of uncertainties that could affect the sizing of ASHPs were quantified, including future weather, future building demand, economic and emission factors, as well as equipment performance variations. A hierarchical optimisation framework is then employed for electricity markets under uncertainty scenarios, in which it determines the baseline power of ASHPs in energy service markets and selects the optimal mode and service type in ancillary service markets. An existing gas boiler heating system, serving educational buildings located at the University of Cambridge, UK, was selected for retrofit analysis. Based on the testing results, some detailed conclusions can be drawn: • From 2022 to 2041 (the next 20 years), the annual mean air temperature would climb by an average of 0.7 • C due to climate change based on the estimates from probabilistic climate models. The annual mean heating load decreased by 5.6%, from an average of 63.4 kWh/ m 2 to 60.0 kWh/m 2 due to the impacts of global warming. • The hybrid operation of air-source heat pumps and gas boilers can provide both energy and ancillary services, such as load covering, frequency containment reserve and frequency restoration reserve. Particularly, the energy revenue accounted for 1.1% of lifespan mean operation cost savings. • The best retrofit alternative of the hybrid heating system offered significant lifecycle performance, achieving an average of 88% lifespan carbon emission reduction, a 54% lifespan total costs reduction, and a payback period of around 3 years. • Retrofitted air-source heat pumps can also mostly meet the future cooling demand to address intensifying concerns over the mitigation of indoor overheating in light of global warming. • The hybrid heating system provides a three-win solution for various stakeholders: carbon reductions for society, cost savings for property holders and improved grid stability for utilities.