Supporting Decision Making for Building Decarbonization: Developing Surrogate Models for Multi-Criteria Building Retrofitting Analysis

Abstract

Decarbonizing buildings is crucial in addressing pressing climate change issues. Buildings significantly contribute to global greenhouse gas emissions, and reducing their carbon footprint is essential to achieving sustainable and low-carbon cities. Retrofitting buildings to become more energy efficient constitutes a solution. However, building energy retrofits are complex processes that require a significant number of simulations to investigate the possible options, which limits comprehensive investigations that become infeasible to carry out. Surrogate models can be vital in addressing computational inefficiencies by emulating physics-based models and predicting building performance. However, there is a limited focus on investigating feature engineering and selection methods and their effect on the model’s performance and optimization. Feature selection methods are considered effective with interpretable models such as multi-variate linear regression (MVLR) and multiple adaptive regression splines (MARS) for achieving stable prediction stability. This study proposes a modelling framework to create, optimize, and improve the performance of surrogate predictive models for energy consumption, carbon emissions, and the associated cost of building energy retrofit processes. The investigated feature selection methods are wrapper and embedded methods such as backward-stepwise feature selection (BSFS), recursive feature elimination (RFE), and Elastic Net embedded regularization in order to provide insights into the model’s behavior and optimize the model’s performance. The most accurate surrogate models developed achieved a mean absolute percentage error (MAPE) of 0.2–1.8% compared to the used test data. In addition, when calculated for a million samples, all developed surrogate models reduced the computational time by one-thousand-fold compared to physics-based models. The study’s findings pave the way towards low-computational accurate models that can comprehensively predict building performance in near real-time, ultimately leading to identifying decarbonization measures at scale.

1. Introduction

1.1. Background

In the face of escalating climate change impacts, the urgency for implementing policies that mitigate greenhouse gas (GHG) emissions has been underscored [1]. Urbanization, progressing rapidly, has positioned cities as the epicenter of this global issue. By 2050, cities are projected to accommodate 67% of the world’s population [2]. Nowadays, cities are the primary consumers of global energy, accounting for 70% of the total, and are the source of over 70% of global carbon dioxide emissions [3]. Buildings are one of the most significant entities shaping a city’s performance. Thus, buildings are also expected to play a crucial role in the sustainable transformation of cities [4,5]. The commercial and residential sectors are responsible for approximately a quarter of the country’s energy consumption in the Canadian context. Buildings alone contribute 13.06% of Canada’s carbon dioxide emissions [6].

In Quebec, buildings consume 35% of the total energy produced, while in Montreal, it peaks at 88% of the total electric energy consumption. Commercial buildings utilize 47% of carbon-emitting fuel sources representing the highest category among buildings in terms of GHG emissions in Quebec [7]. In addition, most of Montreal’s building stock dates back to pre-1970 [8], making the buildings outdated in their building code compliance and energy use efficiency. Therefore, numerous initiatives and incentive programs aimed at building owners have been developed by policymakers, targeting preliminary and feasibility analyses that would accelerate building retrofits [9,10].

BEMs play a crucial role in assessing decarbonization approaches and are deemed essential to examine building retrofit strategies, which may allow decision-makers to implement effective actions. Generally, building energy modelling within the building retrofit domain can be classified into top-down and bottom-up approaches. The latter is often identified as the more adequate for engineering applications due to the granular representation of individual buildings compared to aggregated statistical analyses [11].

Among the bottom-up models are white-box models, which are detailed physics-based and physics-driven, and aim at representing the different building components physically within the model, leading to precise simulations. EnergyPlus and TRNSYS are well-known white-box modelling tools widely used in academia for their accuracy [12]. Nevertheless, the accuracy of these models brings drawbacks of complex inputs and high computational expenses associated with high-fidelity simulations. Moreover, setting up a building energy model through the referred tools can be time-consuming and cumbersome, requiring significant person-hours. Comprehensive building performance analyses involve many variables requiring substantial computational time and power resources, which physics-based modelling techniques may obstruct due to their computational inefficiency [13].

Zhang et al. classified building retrofit measure identification methods into three categories: parametric-based, optimization-based, and machine learning-based methods. Parametric-based approaches rely on pre-defined selected alternatives based on expert domain knowledge, which limits the search space [14]. A study by Saad presented a parametric-based analysis of possible passive design retrofit interventions for existing Egyptian residential buildings reporting a stepped approach selection, where best-performing configurations were selected in each stage [15]. Feng et al. [16] presented a parametric-based approach composed of six building retrofit scenarios aimed at helping homeowners to choose the most suitable retrofit plan considering their building conditions and personal preferences. A more comprehensive approach aimed at understanding the effect of different variables on an office building was performed by Charles et al. [17]. The study examined each building component from infiltration, energy system upgrades, and building envelope on the operational carbon and energy consumption by only investigating one combined scenario. This approach can become infeasible and requires exponential computational capabilities when many building retrofit measures (BRMs) are within scope.

Optimization approaches tackle this problem by introducing the compounded effect of connecting simulation software with an optimization algorithm. Metaheuristic algorithms, such as a genetic algorithm, can identify nearly optimal solutions that balance the objectives. This modelling typology must undergo many runs to achieve the solution and follows preset stopping criteria [18]. Several studies have outlined this problem as a significant drawback in optimization approaches when combined with a whitebox model [19].

To tackle the shortcomings of the former approaches, machine learning models (MLMs) are deployed to emulate the physics-based methods and provide a quick approximation outcome. The MLM is trained on sampled databases which are the outputs of the white-box modelling process, to provide a greater efficiency by minimizing the number of required simulations runs, referred to as surrogate modelling since the MLM acts as a surrogate to the white-box model [20]. The referred setup of models is mainly used to overcome the problem of inexistent pre- and post-retrofit building data [11]. Surrogate models use various learning algorithms, such as artificial neural networks (ANN), support vector machines (SVM), random forest (RF), multi-variate linear regression (MVLR), and multiple adaptive regression splines (MARS) [21]. In summary, surrogates can adequately balance computational cost and accuracy.

The determination of the surrogate model type is predominantly influenced by the pursuit of attaining the highest possible accuracy. Nevertheless, there are instances where a compromise is sought between optimizing accuracy and favoring a model structure that exhibits maximum interpretation [21]. As shown in Table 1, multiple studies have used MARS and MVLR algorithms due to their simplicity and interpretability in addressing their research problems [22,23]. The summary provides an analysis of the current status of the literature in the field for studies that have used MVLR or MARS algorithms.

Regarding the number of samples, it varied drastically from 90 to 10,000 samples. All of the studies with a low number of samples studied a single objective. Xu et al. [24] developed a surrogate model to determine building upgrades and operating costs. At the same time, Wei et al. [25] developed a surrogate model to identify the effect of different building forms on the overall energy consumption of the building in cold climates. Only a few studies used feature selection methods with the developed surrogate models, such as forward-stepwise feature selection. Hygh et al. [23] trained MVLR surrogate models using 16,000 samples from a medium-sized DOE archetype [26] to predict the annual total heating and cooling loads. The study used stepwise feature selection and demonstrated its importance in improving the model fit to the datasets.

Using variable training sample sizes, Chen et al. [22] used an MVLR surrogate model to determine indoor comfort levels and daylighting availability in a high-rise residential building. Prada et al. [27] used a MARS algorithm to develop a surrogate optimization model for several residential typologies reporting the model’s efficiency and effectiveness. Sekhar Roy et al. [28] reported the fractional computation associated with the flexibility of the MARS algorithm in developing a regression mode for heating and cooling load forecasting. It can be observed from the literature that no studies have developed surrogate models with the purpose of early-phase building retrofit measures (BRMs) analysis through energy, carbon emissions, and cost. In addition, few studies have used feature selection methods, and a comparison among different feature selection methods has not been conducted.

Finally, building energy retrofits encompass a broader array of BRMs, which increases the complexity and widens the solution space leading to a sophisticated, complex task usually performed by consultants [29]. The broad range of experts relies on experience-based recommendations or a limited number of investigations due to the limited time and the exponential growth that occurs when a new BRM is added. In addition, the building retrofit process is prolonged and can extend over many years, adding complexity to non-technical users’ initial decision-making process. Therefore, developing an affordably computational modelling method that can provide a more comprehensive analysis with many BRMs could advance the body of knowledge and provide outstanding support in decision-making analysis, which can significantly progress efforts towards building decarbonization.

Table 1. Related studies and identified objectives and features.

Ref.	Study Objective					Case Study		Surrogate Model Parameters
	E	M	Ce	Co	D	Building Type	Simulation Tool	Number of Inputs	Number of Outputs	Algorithm	Feature Selection	Validation	Metric	Sampling Method	Number of Samples
[30]	X	X				O	EP	7	2	MVLR	Best subsets	-	-	-	-
[27]	X					R	TRNSYS	6	2	MARS, MVLR	-	-	-	Multiple	2710
[31]		X			X	O	Daysim	15	2	MVLR	-			Monte Carlo	1900
[22]	X					R	EP	9	3	MVLR, MARS, SVM	-	-	-	-	5610
[32]	X	X				O	EP	26	4	MVLR	-	-	-	LHS & Monte Carlo	467
[33]		X		X		C	EP	5	2	MARS	-			NA	7776
[24]				X		O	EP	6	2	MVLR	-	K-Fold CV	RMSE	-	90
[34]	X					C	EP	16	2	MARS, SRC	-	CV	R2	LHS	2000
[25]	X					O	-	8	3	MARS, RF	-	-	RMSE & R2	LHS	100
[35]	X	X			X	O & I	BSim	14	5	MARS, MVLR	-	CV	R2	Monte Carlo	10,000
[36]	X					R	EP	34	1	MVLR	Forward Stepwise	K-Fold CV	RMSE& R2	LHS	6000
[23]	X					O	EP	27	1	MVLR	Forward Stepwise	CV	RMSE&R2	Monte Carlo	20,000
[37]	X	X				R	EP	25	4	MVLR	BE & RFE	-	RMSE&R2	-	8760

E: energy, M: miscellaneous, C_e: carbon emissions, C_o: cost, D: daylighting, CV: cross-validation, EP: energy plus, C: commercial, O: office, I: institutional, R: residential, BE: backward elimination, RFE: recursive feature elimination.

Table 1. Related studies and identified objectives and features.

Ref.	Study Objective					Case Study		Surrogate Model Parameters
	E	M	Ce	Co	D	Building Type	Simulation Tool	Number of Inputs	Number of Outputs	Algorithm	Feature Selection	Validation	Metric	Sampling Method	Number of Samples
[30]	X	X				O	EP	7	2	MVLR	Best subsets	-	-	-	-
[27]	X					R	TRNSYS	6	2	MARS, MVLR	-	-	-	Multiple	2710
[31]		X			X	O	Daysim	15	2	MVLR	-			Monte Carlo	1900
[22]	X					R	EP	9	3	MVLR, MARS, SVM	-	-	-	-	5610
[32]	X	X				O	EP	26	4	MVLR	-	-	-	LHS & Monte Carlo	467
[33]		X		X		C	EP	5	2	MARS	-			NA	7776
[24]				X		O	EP	6	2	MVLR	-	K-Fold CV	RMSE	-	90
[34]	X					C	EP	16	2	MARS, SRC	-	CV	R2	LHS	2000
[25]	X					O	-	8	3	MARS, RF	-	-	RMSE & R2	LHS	100
[35]	X	X			X	O & I	BSim	14	5	MARS, MVLR	-	CV	R2	Monte Carlo	10,000
[36]	X					R	EP	34	1	MVLR	Forward Stepwise	K-Fold CV	RMSE& R2	LHS	6000
[23]	X					O	EP	27	1	MVLR	Forward Stepwise	CV	RMSE&R2	Monte Carlo	20,000
[37]	X	X				R	EP	25	4	MVLR	BE & RFE	-	RMSE&R2	-	8760

E: energy, M: miscellaneous, C_e: carbon emissions, C_o: cost, D: daylighting, CV: cross-validation, EP: energy plus, C: commercial, O: office, I: institutional, R: residential, BE: backward elimination, RFE: recursive feature elimination.

1.2. Research Objectives and Contributions

Building retrofits are comprehensive in scope and require extensive analyses of the effect of BRMs. The resultant combinatorial problem can be significantly costly computationally for the decision-makers’ initial assessment, specifically for a comprehensive evaluation that includes energy consumption, carbon emissions, and the associated cost. The problem results in a limited scope of the undertaken investigation that is hindered by computational capabilities. In addition, surrogate models that represent a possible solution for the above problem have been limited in the analysis with no focus on the combination of energy, carbon emissions, and cost, representing critical facets of decision making. The survey of existing literature revealed a notable deficiency in research examining the impact of feature engineering and selection methods within the realm of building performance. This gap in understanding leaves unexplored opportunities for optimally configuring models with fewer input features.

It can be noted that the development of surrogate and conventional building performance models has been explored in previous studies, as shown in Table 1. However, studying the effect of model parameters on the surrogate model accuracy is limited in the building research domain. To the authors’ knowledge, no studies have investigated a comparative analysis of feature selection methods to identify the best feature selection method per target output for multiple facets of the building performance. In addition, few or no studies have investigated the modelling energy, cost, and carbon emissions altogether.

Therefore, this study aims to investigate the applicability of developing a surrogate model to emulate the conventional building performance model and expand its prediction to include cost and carbon emissions. The study aims to explore the performance of interpretable surrogate models such as MVLR and MARS algorithms and, in that process, provide a comprehensive investigation of the effect of the different feature engineering and selection methods to develop a clear understanding that can be recreated and developed in future studies.

In this study, a bottom-up conventional building performance model is developed for a reference building in Montreal, Canada. A robust high-fidelity simulation model is designed to report building performance when examining a wide-ranging BRM. In the meantime, the data generated by the physics-based model are used to develop a surrogate building performance model that considers multiple objectives, namely energy consumption, carbon emissions, and cost.

The paper’s scientific contribution can be divided into the following:

• The study incorporates many investigated BRMs and HVAC systems to study their effect on building performance, laying the path for a comprehensive surrogate model. The combinatorial effect of the investigated complex problem is determined, which helps to illustrate the evident deficiency of physics-based simulation models regarding computational requirements.
• The study proposes a methodology for developing interpretable surrogate models for predicting multiple facets of the building performance, while investigating optimum configurations through feature engineering and selection coupling.
• The study presents a comprehensive analysis of the effects of various surrogate model parameters, such as feature selection, on the models’ accuracy and evolution.
• The study proposes novel adjusted feature selection methods that are modified to optimize the performance of the models, stabilize their prediction capabilities, and address their deficiencies.

Considering the addressed research problems and contributions, this paper can be regarded as a value-added scientific paper to the existing body of knowledge.

2. Methodology

The research methodology developed has been divided into three main phases:

Phase (1) includes developing a building performance model for a reference building in Montreal. The process consists of onsite data collection, model setup, and calibration according to commonly used standards.

Phase (2) involves developing the BRM matrix and the included ranges. The alternatives fall within multiple categories, such as external walls or HVAC systems, and the different combinations of the different BRMs represent the possible building retrofit scenarios, each of which is assessed.

Phase (3) involves developing the surrogate models, which can be divided into feature engineering, feature selection, model fitting, hyperparameter optimization, and model evaluation.

2.1. Case Study Information

This study focused on a building complex along the Lachine Canal in the city of Montreal, built during the 19th century [36]. Historically, the region retains a rich history, having hosted many industrial complexes over the years. The area encompassing these building complexes’ typology covers approximately five million square feet. Over the years, however, the industrial building stock has been transformed for various purposes [38].

The selected building was initially constructed in 1908 for textile production and has been repurposed for commercial purposes, with most of them being office spaces. Information was collected from the facility evaluation and through onsite visits. The input features for the building envelope and mechanical systems of the case study were also part of this data collection process. In addition, the space usage was identified through an onsite space inspection and data provided by the building management. The building’s historical energy consumption was collected to be analyzed and to quantify the actual regular building energy consumption. Multi-year energy consumption was used to create an energy usage profile.

A conventional BEM was developed using the onsite collected data and in a multizone configuration to accurately represent the building.

Table 2. Input data for the initial building energy model, including collected onsite or assumed data.

Category		Unit	Base
Building Envelope
	External wall’s thermal resistance	m2K/W	0.83
	Roof thermal resistance	m2K/W	3.78
	Ground floor thermal resistance	m2K/W	1.49
	Window thermal transmittance	W/m2K	2.66
	Infiltration at 50 Pa	ACH50	3.25
HVAC System
	Heating fuel	Category	Natural Gas
	Heating efficiency	%	85%
	Heating setpoint	°C	21
	Cooling fuel	Category	Electricity
	Cooling CoP	Ratio	3
	Cooling setpoint	°C	25
Internal Loads
	Occupancy density	people/m2	0.04
	Equipment density	W/m2	7.5
	Lighting power density	W/m2	8.5

comprises the building data used as an initial input for the BEM, while in the event of a lack of onsite information or testing capabilities, the standard National Energy Code of Canada (NECB) [39] was used as a replacement. Designbuilder software was used to create and simulate the BEM since it has been widely used in previous studies [15,40].

The BEM was subjected to a calibration procedure utilizing the collected energy consumption data. The calibration process is essential to investigate the possible improvements accompanying the BRM’s application. The profile of energy consumption in the building was analyzed across diverse fuel types to discern various load classifications and usage trends. In simpler terms, through this process, we could recognize different types of energy use, like heating or cooling energy, and electricity use, which can change based on the time of year.

An iterative calibration methodology was employed on the model’s indeterminate parameters, such as infiltration and occupancy density rates, heating and cooling setpoints, and lighting density rates, as described in

Table 3. Model parameters and ranges used for model calibration.

Input	Unit	Min Value	Max Value
Occupancy	people/m2	0.01	0.111
Electric Equipment Density	W/m2	4	14
Lighting Power Density	W/m2	4	14
Heating Efficiency	%	50	85
Setpoint Temperature	°C	20	24
Cooling CoP	Ratio	2.5	3.5
Setpoint Temperature	°C	24	28
Infiltration	ACH50	0	20

. Boundary values, both minimum and maximum, were established for the calibration process search. The range of these values was determined by referencing studies conducted within analogous Canadian climate zones and building typologies [17,41]. The collected energy usage data were used for the reference calculation for statistical validation indices such as IPMVP, ASHRAE Guideline 14, and FEMP [42,43,44]. The actual building energy consumption by fuel type compared to the output of the calibrated model is presented in Figure 1. The calibrated model yielded a total coefficient of variation in the root mean square error (CV-RMSE) of ~7.5%, which falls within the acceptable threshold. The calibrated model resulted in a standardized total energy consumption per floor area of 231.4 kWh/m²/year.

2.2. Building Retrofit Measures (BRMs)

A deeper analysis of the building typology (post-industrialized structures) yielded several specific constraints, such as limitations on its building envelope components, energy systems, and distribution systems. First, in the mentioned region, industrial buildings are commonly large uninsulated structures that have transformed over the years to different uses. The initial use type of these buildings and high-temperature industrial production necessities required heating distribution systems that operate using steam as a medium.

Considering the building constraints, suitable BRMs were identified as described in

Table 4. Selected building retrofit measures categories used for developing the building retrofit analysis.

Category	Number of Options	Measure Description	Parameter Unit	Min	Max
External Walls (W)	31	Adding external walls insulation to the building	m2K/W	0.99	8.00
Roof (R)	22	Adding roof insulation to the building	m2K/W	1.56	6.27
Glazing (G)	23	Upgrading the windows of the building to improved thermally resistant options	W/m2K	3.09	0.78
Window Blind (WB)	5	Using window blinds	Slat width (m)	0	0.05
Window Frame (WF)	5	Upgrading the window frame options with less conductivity	m2K/W	3.47	5.88
Window-To-Wall Ratio (WWR)	∞	Adjusting the window-to-wall ratio to affect heating requirements	%	20	95
Infiltration (I)	∞	Improving the air tightness of the building	ACH50	5.8	0.6
Indoor Lighting (L)	14	Replacing the indoor lighting equipment	NPD (W/m2)	3.3	7.6
HVAC Systems (HV)	9	Replacing the existing HVAC system with a high-efficiency system	Heating CoP	0.83	3.3

. Through the review of the literature of standard measures, a set of BRMs were selected [45,46,47,48] due to their suitability to the building limitations, current usage type, and their influential role in maximizing the energy efficiency of office buildings in the studied climate. The workflow adopted to populate the model with a set of BRMs is described in Figure 2, which includes a fuel emission adjustment to represent the local grid and natural gas emissions and the associated tariff. In addition, the accumulated carbon emission residual and cost-effectiveness were calculated for each scenario using life cycle assessment (LCA) and life cycle cost analysis (LCCA) methodologies. Each BRM category will be described in the following subsections, including the minimum and maximum value of one of the features within each BRM category.

2.2.1. Building Envelope

The external walls and the building roof are responsible for significant heat gain and loss of buildings, therefore commonly considered for building retrofits. The suggested parameters were designed to utilize the current building fabric in place and provide the improvement that diminishes the building energy consumption. In the wall and roof categories, 31 and 22 different alternatives were designed for the wall. In detail, the main structure available in the building walls was retained, which, as mentioned earlier, is uninsulated. However, adding different insulation types suitable for the building was proposed with a variation in their thickness. The insulation options (extruded polystyrene, expanded polystyrene, glass fiber, etc.) were combined with the insulation’s thickness parametrization. Insulation was considered to be internally placed in the case of the wall. The required minimum RSI value from the Canadian National Energy Code of Building (NECB) [39] was used as a guide to select the insulation levels. In the process, insulation materials’ costs and carbon emissions were gathered locally [49,50] and updated to the LCC and LCCA calculation modules.

2.2.2. Fenestrations

Fenestrations are drastic in moderating the indoor space dynamics, such as heat gains/losses, required artificial lighting, and conductivity of the building envelope. Therefore, multiple BRM categories that compose the fenestration of the building were considered, such as WWR, glazing, window frame material, and window blinds. For the building glazing, 23 different alternatives were included. Considering practicality, only air-filled windows were selected due to the high cost of gas-filled alternatives. In addition, a low-e glazing parameter was included due to their effectivity besides a variation in the number of glass panes (double, triple, and quadruple glazed). Considering the office space nature, the window-to-wall ratio parameter of 20% to 95% was selected. Five window blind options were considered, ranging in slat reflectivity and slat width. The conductivity of the window frame is critical in determining the efficiency of the whole fenestration and is commonly considered when building retrofit processes occur. Therefore, five window frames, such as conventional and adjusted wooden frames and aluminum with and without thermal breaks, were included. A comprehensive international database was used for the LCC and LCCA modules [51].

2.2.3. Internal Gains

In Canadian office building environments, artificial lighting consumes around 12% of the total energy. Indoor lighting affects the indoor environment and the HVAC system’s sizing. Therefore, 14 alternatives for indoor lighting were studied accordingly [51].

2.2.4. Infiltration

Infiltration is a dominant and influential factor in building performance in cold climates, which refers to the building envelope’s air leakage. According to several studies, infiltration can account for around 15% of the total heating load of the building [52]. Regarding the addressed building, infiltration is prevalent in old structures, especially since the original usage was industrial. Therefore, the building retrofit design space included an infiltration parameter. The thresholds to be potentially implemented used the calibrated infiltration value as a maximum bound and the Passive Haus Institute (PHI) [53] reference of 0.6 ACH₅₀ as a minimum bound.

2.2.5. HVAC Systems

HVAC systems play the most critical role in the overall energy consumption of the building. In the studied case, a steam-based distribution system is employed historically in the building, which affects the efficiency of the overall heating capacity. However, alternative HVAC systems would require a change in the distribution system. Therefore, a set of HVAC options are used to mimic several possible scenarios of the building transition.

Table 5. HVAC alternatives’ descriptions selected for the BRM matrix.

HVAC Systems	ID	Description
Packaged Terminal Air Conditioning with Boiler (PTAC-g)	HV 1	Packaged heating and cooling where a natural gas boiler supplies the heat.
Packaged Terminal Air Conditioning with Electric Heating (PTAC-e)	HV 2	Similar to the PTAC-g but uses electric resistance heating elements to provide heat.
Unitary HVAC (UNI)	HV 3	Comprehensive heating, ventilation, and air conditioning system that covers a large zone or multiple zones and includes heating and cooling altogether.
Packaged Terminal Heat Pump (PTHP)	HV 4	A self-contained heating and cooling unit that can switch between cooling and heating modes as required.
Variable Refrigerant Flow with HR (VRF)	HV 5	An advanced HVAC technology that precisely controls refrigerant flow to multiple zones, providing heating, cooling, and ventilation.
Natural Gas Condensing Boiler with Steam Distribution with a DX cooling unit	HV 6	A condensing boiler heating system that uses natural gas as a fuel and distributes heat throughout a building via steam.
Natural Gas Condensing Boiler with Water Distribution with a DX cooling unit	HV 7	A condensing boiler heating system that uses natural gas as a fuel source and distributes heat via hot water.
Electric Condensing Boiler with Steam Distribution with a DX cooling unit	HV 8	A condensing boiler heating system that uses electricity as a fuel and distributes heat throughout a building via steam.
Electric Condensing Boiler with Water Distribution with a DX cooling unit	HV 9	A condensing boiler heating system that uses electricity as a fuel source and distributes heat via hot water.

introduces variations in the investigated HVAC systems. HV1 and HV2 include a standard packaged terminal air conditioner unit (PTAC) [54] used for cooling and heating or only with a natural gas condensing boiler, which is a ductless system. HV3 comprises a unitary system that uses an electric resistance for heating and a DX unit for cooling, which shares many similarities with the PTAC-e system, but is more extensive, has a larger capacity, and uses ducts [55]. HV4 is a packaged terminal heat pump system (PTHP) that is ductless and commonly used in office buildings [55]. HV5 is a variable refrigerant flow (VRF) with a heat recovery (HR) system. It can reclaim heat from cooled areas and transfer it to areas requiring heat, thus allowing simultaneous heating and cooling operations [56]. HV6–HV9 consider conventional alternatives that rely on electric or natural gas condensing boilers. The alternatives include two possible scenarios keeping the existing inefficient distribution system or upgrading it to water-based distribution. The consideration of the distribution variation is represented in the heating CoP of the HVAC system according to the guidelines [57].

2.3. Surrogate Model Development

This section discusses the process of the surrogate model development. Westermann et al. [21] thoroughly explained the surrogate model development steps, including defining the model’s intent, designing the objective, and determining the design parameters suitable for the addressed problem. As described earlier, the main objective is investigating the feasibility of developing a data-driven approximative model that emulates the physics-based model for building retrofit preliminary assessment. The aim is to define three categories of models that imitate the calculation of the building’s heating energy consumption, embodied carbon emissions, and cost. The following subsections will discuss a more detailed overview of each step of the surrogate model development. The developed model and the steps used within the methodology of this study are presented in Figure 3.

2.3.1. Dataset Development

First, a dataset that includes the parameters (BRMs) and their corresponding performance indicators is developed. Sampling methods highly influence a surrogate model’s performance and efficiency in generating the required design space. The sampling method’s efficiency can be described by maximizing the knowledge of the design space and its parameters. Commonly used sampling algorithms are Latin hypercube (LHS), Monte Carlo (MC), and Sobol [21,58].

LHS is a pseudo-random method that extends the Latin square and generates one sample per axis-aligned hyperplane. In this study, Latin hypercube was deployed due to its wide use in the literature and its capability to distribute samples over the entire design spaces with less computation required, such as when using Monte Carlo sampling [59,60,61]. The output of the sampling method is referred to as the target vector(s) and feature matrix, as shown in the following equation. The independent variables are referred to as features, while the predicted variables are the target vector.

Feature matrix        Target vector

(Independent Variables)    (Dependent Variables)

$$\left[\begin{array}{ccc}{x}_{1,1}& \cdots & x_{n,1}\\ ⋮& \ddots & ⋮\\ x_{1,n}& \cdots & x_{n,n}\end{array}\right] . \left[\begin{array}{c}{y}_1\\ ⋮\\ y_n\end{array}\right]$$

2.3.2. Dataset Pre-Processing

The MVLR and MARS models use numerical values to fit the model. Therefore, the categorical features were transformed into numerical component variables, as shown in

Table 6. Converted original categorical features to numerical variables.

	Wall (W)	Roof (R)	Glazing (G)	Window Blind (WB)	Window Frame (WF)	Indoor Lighting (L)	HVAC (HV)
ID	Generated Features
1	thermal conductivity	thermal conductivity	thermal conductivity	Blind-to-glass distance	thermal conductivity	power-density	heating-fuel
2	thermal resistance	thermal resistance	light-transmission	width	thermal resistance	radiant-fraction	auxiliary-energy
3	total-thickness	total-thickness	solar-transmission	cost/window/slat width	specific-heat	control-type (linear)—encoded	pressure-rise
4	cost/area/rvalue	cost/area/rvalue	shgc		thickness	control-type (off)—encoded	distribution-efficiency
5	insulation-thickness	insulation-thickness	gap-thickness		cost/area/uvalue	control-type (stepped)—encoded	heating cop
6	insulation-uvalue	insulation-uvalue	cost/area/uvalue			cost/area/power density	cooling cop
7	insulation-specific-heat	insulation-specific-heat	pane-thickness				cost/area/heating cop
8	insulation-density	insulation-density	number-panels				cost/heatingload
9	insulation-cost/volume	insulation-cost/volume/uvalue
10	insulation embodied carbon/uvalue

, resulting in 49 features. The variable’s selection was based on best practices [21] in the literature and domain knowledge. Variables that are categorical, such as “heating-fuel” and “control-type”, were encoded with one label encoding and one hot encoding, respectively. After converting categorical features to numerical features, the total number of features was 51 when including WWR and infiltration.

2.3.3. Model Selection

The process includes the surrogate model’s construction steps: model type selection, feature analysis and pre-processing, training, and validation processes. The model selection has been given significant importance due to its critical effect on prediction accuracy and overall model performance. Interpretability is a standard criterion for model selection in engineering problems. For example, deep learning models and their black box processing nature limit the model behavior understanding and further development. MVLR and MARS are among the most interpretable algorithms. The capabilities of both algorithms differ in their flexibility. For example, MARS is a more flexible algorithm that can create a contour plane to model a problem compared to the single plane of MVLR, as shown in Figure 4. This study employed both algorithms to investigate their different capabilities within the studied domain.

2.3.4. MVLR Model Development

More information about the MVLR model development will be highlighted in this subsection. MVLR is essentially a linear equation such as (Y₁ = β₀ + β_1×1 + β_2×2 +....+ β_nX_n + ε). The feature’s importance is composed in (β), referred to as the feature magnitude. The nature of the algorithm facilitates identifying the principal features per problem. However, MVLR models may prove inadequate due to several limitations, such as multicollinearity, sensitivity to outliers, and the need for normal distributions within the data. A set of tools were used for the MVLR model development, such as the Python-based libraries of Scikit-learn [62], Scipy [63], Statsmodels [64], Pandas [65], and NumPy [66].

Therefore, while developing the model, necessary steps and strategies to limit and avoid the model’s weaknesses were considered. The process starts with the normalization of the features such that they are within similar value ranges, which prevents model bias towards a feature due to its high numerical value. Standardization applies a standard deviation of one and a mean of zero on all features. Feature transformation techniques are commonly used in regression problems using mathematical formulas to improve the prediction capabilities and precision of the model. Transformation functions allow a better fitting of the model to the data so that it can explain more of the data variance. Methods used include quantile, power, log, exponential, and cubic root transformations. An automated process involving a brute-force investigation of all possible combinations of transformations was developed in the model development, reaching 30 combinations per model. The process transforms the data and later reverses the transformations during the evaluation process of the model.

Feature Selection Methods

A significant and influential step in developing a stable MVLR model is selecting the features that provide the best performance, generalize prediction capability on unseen data, and tackle multicollinearity. Among the methods are wrapper and embedded methods, which have two different approaches to addressing the shortcomings [67]. The former selects a smaller subset of the feature matrix that can be used as performance indicator criteria and eliminates the other features. The selected feature selection and regularization methods are presented in Figure 5 for developing models M2–M4.

A.

Backward Stepwise Feature Selection (M2)

One of its methods is backward stepwise feature selection (BSFS), which starts with a whole feature matrix and ejects the least significant variable for each step while improving the model’s precision. This study examined the efficacy of using BSFS combined with a variance inflation factor (VIF) to remove highly collinear features. An automated workflow eliminates features through a stepwise process while considering multicollinearity based on the value of VIF > 10. The ejection process works stepwise, where during each iteration, the recalculation of the VIF occurs until no feature surpasses the max threshold. The selected subset of features is then used to fit the model.

B.

Recursive Feature Elimination (M3)

A more computationally exhaustive method is recursive feature elimination (RFE), which uses an optimization technique combined with a preset number of features to eliminate every step. The RFE model uses cross-validation subsets to identify and eliminate the least important features. Although the RFE model is considerably more computationally expensive, it was employed in this study for comparison. A novel RFE feature selection process was developed by initiating the model development with an optimum number of feature searches that maintained the R² for each model, followed by a search for the best combination of features. Generally, feature selection methods provide model coherency and optimization, where fewer features can describe the data more efficiently.

C.

Elastic Net Regression (M4)

Another approach is embedded methods, which seamlessly include a feature regularization process within the model construction process. Regularization refers to limiting the coefficients of the features while fitting the model to prevent overfitting and stabilize its prediction capability. It addresses multicollinearity by shrinking the coefficients of correlated features, which can sometimes reach zero, similar to eliminating the feature. Models such as Lasso, Ridge, or the combined form, Elastic Net, use regularization by adding penalties to the loss function while fitting the algorithm. The penalty aids a better distribution of coefficient values, hence, spreading the collinear variables effect across smaller coefficients. This process reduces the variance of the estimates, which is highly prone to multicollinearity resulting in a more stable model. This study used Elastic Net to regularize the initial model and test its capability to provide an improved, more precise model [68].

2.3.5. MARS Model Development (M5)

On the other hand, the MARS algorithm is considered to be an extension to MVLR, which uses hinge functions, a piecewise linear approximation, to identify complex relations among features and the output variable while maintaining interpretability by reporting the magnitude of each feature or product feature quantitatively. The model’s non-parametric nature provides flexibility in modelling non-linear relationships, creating complexity and increasing computational needs. The model was developed using the Python-based library PyEarth [69].

MARS utilizes a methodology that generates models using piecewise linear basis functions while creating new features from combinations of the original features. The fundamental notion entails dividing the predictor space using the hinge functions, applying piecewise segmentation within the predictor data, and employing piecewise individual linear regression models to fit each.

The process includes a forward and backward stage in fitting the model and optimizing its performance. In the forward pass, an exhaustive iterative search is conducted to identify the set of functions and values that provide the most significant reduction in the summation of squared residuals when an additional basis function is integrated. The process continues until adding a new basis function does not improve the residual sum of squares. To avoid overfitting, MARS includes a backward elimination mechanism, which prunes basis functions from the model using the least effective functions (reduction in the sum of squared errors). Finally, the series of linear functions are completed, where each applies over different feature space regions. This study includes the development of a MARS model in parallel with the MVLR models to assess the flexibility vs. interpretability tradeoff associated with it while considering its performance.

2.3.6. Error Quantification

This study utilized the frequently employed mean absolute percentage error (MAPE), root mean square error (RMSE), and R². The mean absolute percentage error (MAPE) and the root mean square error (RMSE) are commonly used metrics for evaluating accuracy. The R² metric assesses the explanatory capacity of the model.

The MAPE metric was calculated as in Equation (1) and utilized to compute the mean percentage deviation between the predicted and simulated values. A relative error measure is advantageous in analyzing data encompassing diverse scales or units, enabling the convenient interpretation and comparison of results across disparate models or datasets.

RMSE is computed as the square root of the average squared differences between the predicted and actual values, as in Equation (2). RMSE possesses a distinctive property of penalizing errors of greater magnitude to a greater extent than MAPE, rendering it particularly suitable for models that ascribe greater importance to more significant errors. RMSE is commonly used in diverse fields. It exhibits less sensitivity to outliers than MAPE.

The R² measures the degree to which the variance in the dependent variable can be anticipated based on the independent variables included in the model. The value of R² ranges between 0 and 1. It is calculated through Equation (3), wherein a higher value designates the model’s superior goodness of fit to the data. The stated metric evaluates the degree to which the model delineates the variability observed in the data.

$$MAPE=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left|\frac{y_i-p_i}{y_i}\right|· 100$$

(1)

$$RMSE=\sqrt{\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{\left(y_i-p_i\right)}^2}$$

(2)

$$R^2=1-\frac{\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{\left(y_i-p_i\right)}^2}{\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(3)

• $\overline{y}$: Mean of actual values;
• $y_i$: Actual values;
• $p_i$: Predicted values;
• $n$: Number of data points.

3. Results

Based on the methodology discussed, this section summarizes the results of implementing the steps presented in Section 2. An overview of the dataset developed will take place, followed by the multiple models and their performance metrics.

3.1. Dataset Analysis

The parameters described in Table 4 were sampled using an LHS sampling procedure, creating a dataset of 1000 samples. In Figure 6, histograms of the output variables are presented. Most features reached a normal distribution of the dependent variables providing a representative dataset size.

A skewed distribution is observed in the total carbon emissions (kgCO₂e) histogram analysis due to the drastically different fuel emissions or operational carbon. The reliance of the Quebec electric grid on 99% clean hydro sources results in nearly zero operational carbon when the system is electrified. Therefore, it is significantly lower than fossil-based-operated alternatives regarding carbon emissions.

As mentioned, linear regression surrogate models are prone to extreme values and non-normal distributions. Therefore, developing a model that predicts embodied carbon can be more practical, as the operational carbon can be estimated by predicting the heating energy consumption.

Consequently, predicting the heating consumption was selected to be included among the developed surrogate models since heating consumption accounts for the majority of energy consumption in the building, and it is the only consumption parameter that reflects the fuel change and its carbon emissions, hence allowing the calculation of the operational carbon emissions as well.

The life cycle cost analysis (LCCA) output was used as the target variable for the cost model to predict the total cost associated with the BRMs applied.

Three categories of surrogate models have been developed: the energy model, carbon model, and cost model, using the heating energy consumption, embodied carbon of materials, and LCCA, respectively.

3.2. MVLR Results

The following section presents the MVLR results, including all of the developed models using different feature selection and regularization methods for comparison. Generally, in all of the developed MVLR models, a dataset is split into test and training sets with a ratio of 1:4. Different cross-validation techniques were used among the developed models, where the training dataset was further divided into validation sets. In addition, standardizing the output features by a geometrical unit was performed to allow future comparative analyses. Saad and Eicker compared numerous standardization methods reporting the per conditioned volume area as the least biased within different data types and spatial scales [70]. The study’s feature selection processes have resulted in the use of the features in

Table 7. Features selected in all of the developed MVLR models.

Model	Target Variable	Selected Features	Number of Features
MVLR—Base Model (M1)	All	All features	51
MVLR—Backward Stepwise Feature Selection (M2)	All	WWR, I, G1, G2, G4, G6, G7, G8, R1, R3, R5, R7, R8, R9, HV2, HV3, HV5, HV7, HV8, L1, L2, L6, W3, W4, W6, W7, W8, W9, W10, WB1, WB2, WF3, WF5, HV1, L4, L5	36
MVLR—Recursive Feature Elimination (M3)	E	I, G4, G6, R1, R2, R4, R5, R6, R9, HV2, HV4, HV5, HV6, HV7, HV8, L1, L2, L6, W1, W2, W3, W6, W8, WB2, WB3, L3	26
	Ce	G1, G3, G4, G6, R1, R2, R3, R5, R6, R7, R9, W1, W4, W8, W9, W10, WB2, WB3	18
	Co	R1, R5, R9, HV2, HV3, HV4, HV6, L1, W3, W9, WB2, WB3, L3	13
Elastic Net MVLR (M4)	All	All features	51

.

3.2.1. Base Model Performance

As an initial step, the model fits an initial MVLR model, which provides insights into the model’s capabilities and the possible adjustments needed for the model development. In this model, no cross-validation sets were used. Figure 7 presents the performance of three MVLR base models developed for energy, carbon emissions, and cost prediction. It can be observed that the current structure of the features and sampling space provided an acceptable model performance for the carbon emissions and cost model. The former achieved an R² of 90% and a MAPE of 0.82%, while the latter achieved an R² of 96.5% and a MAPE of 0.74.

On the other hand, the energy model achieved an R² of 92.6% with a MAPE of 14.3%, which is a very high error level. All of the features identified (51 features) were used to develop these models, which can affect the overall model overfit; however, the limited number of samples (1000) would make it unlikely. Within the cost model, the residual values can be observed to follow a more uniform distribution within the range of values.

3.2.2. Backward-Stepwise Feature Selection

A backward-stepwise feature selection process was applied as a first-step examination of possible improvements to the developed models. An automated BSFS process was used, where the feature selection assesses the inflation among the input variables and removes the variables with the highest inflation during each incremental step until the highest VIF is less than or equal to 10 [71], as used by other studies in the field. Figure 8 shows the applied process and the resulting number of features. The number of features after applying the process fell to only 36 input features. The approach employed for model development involved utilizing nested K-fold cross-validation with ten folds. This technique divides the training dataset into ‘K’ sets of equal size. ‘K−1’ sets are combined to create the training set, while the remaining set is designated as the validation set. The used cross-validation method guarantees the generalization of the developed model and capability on unseen data. The following model’s performance can be observed in Figure 9.

The model’s predictive accuracy increased among all of the models. At the same time, a significant improvement was seen in the energy model, which significantly improved in prediction error, decreasing to a MAPE of 2.77% and an RMSE of 0.33 (kWh/m³). The carbon emissions and cost models achieved a predictive accuracy increase as well. The former achieved a MAPE of 0.43% and an RMSE of 0.11 (KgCO₂e/m³), and the latter achieved a MAPE of 0.29% and an RMSE of 1.18 (C$/m³). The R² of all of the models decreased by diminishing the number of features. However, the most significant decrease can be identified in the carbon emissions model, reaching 73.6%, which indicates the need for more samples to develop a more representative model or the need to use an alternative feature selection method to tackle multicollinearity. For the energy and cost surrogate models, the R² achieved was 90.4% and 92.8%, respectively.

3.2.3. Recursive Feature Elimination

Another feature selection method was used to address the diminution of performance of the R², namely RFE. Conventionally, RFE is implemented considering a process of elimination that acquires the highest performance metric indicated with no consideration of the maximum number of input features, which might lead to the complete set of features used. In this approach, we implement a criterion of using the least number of features that yields a stable R². This process balances the performance while providing a method to minimize multicollinearity by reducing the feature number. As a first step, a process of fitting the model and indicating the optimum number of features that stabilize the R² value is determined, as shown in Figure 10. The optimum number for the energy, carbon emissions, and cost models were 26, 18, and 14 features, respectively. It can be noted that a similar technique of nested K-fold cross-validation was used in the model to maintain generalizability on new unseen data.

After determining the optimum number of features, the models were fit to examine their entire performance. Figure 11 presents the performance achieved by the different models. It can be observed that while using a significantly reduced number of input features: 26 for energy, 18 for carbon emissions, and 14 for cost, to be exact, a significant decrease in performance does not occur. In detail, the energy and cost models witnessed a slight performance decrease with an increase in the MAPE of 0.17% and 0.03%, respectively. In addition, a slight decrease in the R² of the models was observed. However, the models maintained a robust prediction precision (<5%). The carbon emissions model achieved a considerable performance improvement across all metrics.

Using only 18 features of the original dataset provided a significant increase of 13% in the R² of the model compared to the BSFS model. However, the BSFS model was based on three times the number of features. In addition, the MAPE of the developed model decreased by 0.12%, reaching 0.31%. The process concludes the efficacy of exploring multiple feature selection techniques per model for performance improvement.

3.2.4. Elastic Net Regression

An alternative route to developing a generalizable model and preventing overfitting is using embedded feature selection methods. Elastic Net combines L1 and L2 regularization techniques to achieve stability within the model. The regularization process adds a penalty term (λ) to the sum of squared coefficients and the sum of absolute values, which reduces features with significant coefficients and shrinks their values. Figure 12 visualizes the process of coefficients’ shrinkage, showing the differences among the three developed models. The process stabilizes the developed models, maintaining the performance when exposed to unseen data. This model used ten-fold cross-validation to optimize the penalty term value. In developing these models, a grid search among hyperparameters, such as the tolerance, L1 ratio, and alpha value, was performed to identify the optimum model.

Figure 13 presents the performance of the Elastic Net models developed for energy, carbon emissions, and cost. It can be observed that an improvement across all of the models’ performance metrics was achieved. The highest precision of the models was achieved compared to the alternative feature selection methods. The RMSE values for the three models are 0.32 (kWh/m³), 0.07 (KgCO₂e/m³), and 1.11 (CAD$/m³). The description of the variance capability of the models is 90.7%, 88.3%, and 93.6%, respectively. The results demonstrate the Elastic Net models’ capability to generalize the prediction precision.

Table 8. Optimized hyperparameters for the developed Elastic Net models.

Model	Tolerance	L1 Ratio	Fit Intercept	Alpha	RMSE Difference among Train and Test Sets
Energy Consumption	1 × 10−10	0.6	True	7.1× 10−5	0
Carbon Emissions	1 × 10−10	0.8	True	1 × 10−5	0
Cost	1 × 10−10	0.2	False	1 × 10−5	0.05

demonstrates the prediction error difference between the RMSE values in the test and training datasets. The energy and carbon emission models performed equally precisely on both datasets, while the cost model had a slight difference of 0.05 (CAD$/m³) from the original 1.11 (CAD$/m³). Although the Elastic Net models used the complete set of input features (51), they achieved a higher performance and demonstrated generalizability on the train and test data.

3.3. MARS Results

The MARS algorithm was used to develop alternative models for comparison. The MARS algorithm is more complex due to its non-linear nature. Using a forward stage-wise fitting process, the algorithm uses a stepped approach to fit the model. As shown in Figure 14, the number of terms used for developing each model stops when the performance threshold reaches a constant level. It can be observed that the energy and carbon emissions models used 20 basis functions to reach the R² peak, while the cost model used 17.

In Figure 15, the model’s performance can be observed to have significantly increased. The overall R² among the three models increased, demonstrating non-linear relations within the datasets. On the one hand, the MAPE of the models decreased, achieving a higher precision; however, the improvement among the carbon emissions and cost models can be neglected due to their precise performance in the linear models. The MAPE decreased by nearly 1% for the energy model, representing a 30% improvement compared to the Elastic Net model. However, the improvements mentioned come with the cost of a complex interpretability. The models use a third-degree polynomial basis function to achieve these performances, which can provide less interpretability and be considered of higher complexity.

3.4. Computation Time

In this section, an analysis of the computational time required was carried out among the different developed models. First, a hypothetical 1 million simulation target was considered to resemble the magnitude of the total combinations of the BRMs that can reach nearly 3 million different combinations. The conventional simulation model would require approximately 567 days to perform the simulations using a 12-core processor. However, the dataset development used in this study was based on 1000 samples which needed roughly 13 h to be carried out. For the development of the surrogate models, the MARS model used the highest computational time for fitting, reaching 15.1 s, followed by the Elastic Net MVLR model with 6.4 s. The efficiency of the surrogate models can be observed in the prediction of 1,000,000 samples, as most models can achieve a nearly real-time prediction. The maximum time taken is among the MARS models, reaching 9.7 s, as shown in

Table 9. Computational time (in seconds) analysis among the developed surrogate models.

	Building Energy Model	MVLR	BSFS	RFE	Elastic Net	MARS
Simulation Time (s)	49,000,000	49,000	49,000	49,000	49,000	49,000
Model Fitting and Testing (s)	0	0.3	3.2	0.9	6.4	15.1
Predicting 1,000,000 samples (s)	0	0.2	0.3	0.2	0.1	9.7

. The comparison is visualized in Figure 16, demonstrating the overall efficiency that can be reached using a surrogate model methodology to approximate the BRMs’ effect on the building performance through the lens of energy, carbon emissions, and cost.

4. Discussions and Conclusions

This study builds on the existing literature investigating the efficacy of developing surrogate models using MVLR and MARS algorithms, adding a methodological approach including the effect of multiple feature selection methods. The work investigated the feasibility of predicting building performance in a building retrofit process, resulting in three categories of models: heating energy consumption predictor, embodied carbon emissions predictor, and the life cycle cost predictor associated when applying the BRMs to a commercial building in Montreal.

The surrogate model’s accuracy peaked between 0.2% and 1.8% MAPE, which shows a very high precision associated with the developed models. The investigated BRMs in conventional settings can result in a combinatorial case count of nearly 3 million cases, which is infeasible to be simulated and reported in conventional settings. Surrogate modelling demonstrates its efficiency and feasibility in providing a real-time prediction output for the mentioned number.

The study’s findings add to the existing literature in the following ways: (1) Comparison of the efficacy of developing surrogate models using MVLR and MARS algorithms using 1000 samples for energy, carbon emissions, and cost predictions. (2) Identification of the effect of different feature selection methods and developing a workflow to tackle multicollinearity and instability in MVLR models. (3) Demonstrating the efficacy of optimizing the total number of features used for developing the models, which not only provides stability to the models and avoids multicollinearity but also allows the opportunity to use the models in early-stage designs.

The computational time required to fulfill 10⁶ combinations for the studied building has been diminished by nearly one-thousand-fold across all of the developed surrogate models. This computational time reduction gap can be increasingly higher when considering multiple buildings or an urban building energy modelling scenario.

Across all the used techniques, it can be demonstrated that there is no one-size-fits-all technique, but a more suitable technique or algorithm based on the target output and input features.

Generally, the more flexible model, MARS, outperformed MVLR models regardless of the feature selection method applied. However, the interpretability versus accuracy tradeoff can be justified as a reason to use MVLR models. The MARS performance was significantly better in the energy and cost predictor models, demonstrating the higher non-linear relations among the input features for acquiring those outputs. The study demonstrated the differences between BSFS, RFE, and Elastic Net feature selection and engineering methods. An automatic feature selection that applies novel criteria for selection was developed.

The methodology demonstrated in this study intends to facilitate the process of early-phase economic and environmental decision-making of BRMs in large buildings, focusing on post-industrialized structures in Montreal. The methodology can be replicated to be applied to other building typologies. In addition, the outputs can be deployed in interactive tools that allow decision-makers to interact and investigate.

Future studies include deploying the developed models to develop building retrofit scenarios for the examined buildings. The coupling between the models and multi-objective optimization (MOO) algorithms can provide a helpful decision-making tool for involved stakeholders. In addition, investigating possible synergies between building archetypes and surrogate models provides a breakthrough for building scale to urban building scales and a more significant impact.