For the development of any solar prediction model, we must contemplate the timescale (counted from a certain moment in which we make the prediction), which determines the future moment for which we make the predictions. The biggest drawback is its linear character, which makes it difficult for all problems to be properly modeled [22,24]. Models are very sensitive to unusual observations, forcing you to review the time series for detection and correcting before designing the prediction model [25].

### 3.1.1. Nowcasting Solar Radiation

The prediction horizon ranges from 15 min to a few hours, with no unanimity in its value [26]. In the short term, intra-hour forecasts are particularly useful for carrying out operations in the solar plant, balancing the grid, achieving automatic generation and trade control. Currently, it is very difficult to accurately predict solar radiation in the short term, as it involves knowing in advance how much energy solar plants produce instantly and this would be of great help in avoiding problems of supplying the line or avoiding surplus energy [27]. On this horizon, the statistical models that show the best performance are

those that use satellite images with greater accuracy in nowcasting predictions [15,28] and statistical models for solar radiation time series, such as NWP [29–31]. Likewise, in recent years, several studies have been presented on the use of Machine Learning (ML) for solar radiation prediction using vector support machines (SVM) that show better performance in classification and regression analysis in time series [32,33].

### 3.1.2. Forecasting Solar Radiation

Long-term predictions: correspond to a horizon above 48 or 72 h, reaching a limit of 7 days. The larger the horizon, the greater the prediction errors, making it difficult to make reliable predictions of atmospheric variables above those 7 days [34–36]. Time horizons represent prediction analysis, as, as has been observed in the literature, ML models are able to demonstrate that they are the best alternative in time series data analysis. However, the need to reduce error correction has been modified to make way for ML regression model models using search optimization algorithms (SOA) that require accurate selection of their parameters to improve their performance [37,38].

### **4. Support Vector Machine Models**

Machine Learning (ML) is the process of learning to convert experience into expertise. A concise explanation is that a computer program learns from the experience from data recorded in a period of time; the performance of the program is evaluated in time by improving experience [39]. In other words, ML models interpret patterns through learning based on data obtained in a period through a training process to make a prediction, generating new data that will measure the behavior of a phenomenon in the future [40–42].

According to the literature, there are two main specific techniques that ML uses to develop learning methods using information directly from data acquired [43]. There exist two main learning methods: supervised learning that uses from present to previous data to forecast events [44,45] and unsupervised learning that analyzes the patterns from non-classified data and deduces a function to describe the behavior of the system [46]. Additionally, there are other learning methods that derive from these principal methods: semi-supervised learning and multiple instance learning [47]. In this article, the supervised learning using a support vector machine (SVM) will be discussed. The Figure 2 summarizes the two principal methods and its most used models to predict solar radiation.

According to [8], supervised algorithms are currently the most used, and within this method are neural networks and vector support machines. These models have been widely used in recent years, with vector support machines being recently integrated to provide new solar prediction techniques. In this context, SVM belongs to a technique called supervised learning defined as a technique for identifying the behavior between an input and an output variable [48].

The theoretical basis of SVM is to minimize the structural risk related to the empirical risk from the training process and the confidence range from Vapnik Chervonenkis dimension (VC dimension) [49]. The complexity of the problem will be reduced when the VC dimension is smaller, making the risk smaller [50].

The architecture of an SVM is a versatile and configurable model based on a kernel machine that could be treated as a classification or regression problem according to Vapnik equations [51], while a support vector regression is only used for regression problems [11]. Therefore, an SVM is a great alternative to solve classification problems in forecasting time series [52–55].

**Figure 2.** Machine learning techniques.

In order to predict a non-linear variable, such as solar radiation, it is necessary to have a data set *xi* that represents the input data sample space, *di* which is the target value and to know the number of data points *n* [14]; therefore, it is possible to resort to Vapnik's equations theory. The input variables to SMV models are related to the objective variable (variable to predict), visualizing the mapped data in a non-linear function *f*(*x*) [51]:

$$f(\mathbf{x}) = \omega \cdot \phi(\mathbf{x}) + b \tag{1}$$

where *ω* is the normal vector, *b* is a constant, also called bias term [56] and *φ*(*x*) is a large-dimensional spatial characteristic mapped by the space vector *x*. The coefficients *ω* and *b* are calculated by minimization using the following optimization problem [24]:

$$R\_{\text{svm}}(f) = \mathbb{C}\frac{1}{N}\sum\_{i=1}^{N} \mathbf{x}\_{i=1} = L\_{\varepsilon}(f(\mathbf{x}\_i), y\_i) + \frac{1}{2}||w||^2 \tag{2}$$

$$L\_{\mathfrak{e}}(f(\mathbf{x}\_{i}), y\_{i}) = \begin{cases} \text{ si } |f(\mathbf{x}), y| - \varepsilon & \text{for } |f(\mathbf{x}), y| \ge \varepsilon \\ 0 & \text{otherwise} \end{cases} \tag{3}$$

where  is a parameter of the model. *Le*(*f*(*xi*), *di*) is the term that describes the -th missing function, which indicates that errors below epsilon are not penalized, *di* represents the solar radiation in the period *i* and *C* 1 *N* ∑*N i*=1 *L* (*f*(*xi*), *di*) defines the empirical error of the SVM model. 1 2 *w*2 is the regularization term, *C* is the term that evaluates the error penalty function to regulate the compensation between the error or empirical risk and the term of regularization. The slack variables *ς* and *ς*∗ indicate the excessive top and bottom skew, respectively. With these properties of the function to be optimized, it is possible to define Equation (2) as shown below [57]:

$$\text{minimize} \quad \frac{1}{2}||w||^2 + \mathcal{C}\frac{1}{N}\sum\_{i=1}^{N}(\mathfrak{c}\_i + \mathfrak{c}\_i^\*)\tag{4}$$

$$\text{subject to:} = \begin{cases} |y\_i - (\langle w, \chi\_i \rangle + b) \ge \epsilon + \zeta\_i| \\ \langle w, \chi\_i \rangle + b - y\_i \le \epsilon + \zeta\_i \\ \zeta\_i \zeta\_i^\* \ge 0 \end{cases} \tag{5}$$

To solve Equation (1), it is possible to use Lagrange and optimal constraints to obtain a non-linear regression function:

$$f(\mathbf{x}) = \sum\_{i=1}^{l} (a\_i - a\_i^\*) K(\mathbf{x}\_i - \mathbf{x}) + b \tag{6}$$

where *αi*, *α*∗ *i* are Lagrange multipliers. The term *k*(*xi* − *x*) is defined as the kernel function [55]:

$$K(\mathbf{x}\_i - \mathbf{x}) = \sum\_{i=1}^{D} \phi\_i(\mathbf{x}) + \phi\_i(y) \tag{7}$$

The general architecture form of a SVM is shown in Figure 3.

**Figure 3.** General architecture of a support vector maching (SVM) model according to [55].

### *4.1. Kernel Function Provided for SVM*

SVM maps the data in a non-linear map to describe the linear process in the space of the predicted data according to the availability of data. This expression results in a simple linear combination problem attributed to the mapped space [58]. The kernel function allows a classification to be performed to form nonlinear boundaries to model complicated separating hyperplanes [59]. Figure 4 describes a projection from low to high dimension in space data.

The kernel parameters must be ideal to solve the problem classification according to the data to become separable in the next space. The four principal basic kernel functions are linear, polynomial, radial basis function (RBF) and sigmoid [57,60].

• Radial basis function (RBF): this function could perform nonlinear mapping of the samples into a higher dimensional feature space expressed by [61]:

$$K(\mathbf{x}\_{i}, \mathbf{x}\_{j}) = \exp\left(-\frac{||\mathbf{x}\_{i} - \mathbf{x}\_{j}||^{2}\sigma}{2}\right) \tag{8}$$

where *σ* is the kernel weight and *xi* and *xj* are the inputs to the *i*-th and *j*-th dimensions, respectively.

• Linear kernel function: According to [57], the linear function to obtain the SVM parameters is described by the following:

$$K(\mathbf{x}\_{i\prime}\mathbf{x}\_{j}) = \mathbf{x}\_{i\prime} \cdot \mathbf{x}\_{j} \tag{9}$$

• Polynomial kernel function: this is a typical example of a global kernel, defined as follows:

$$K(\mathfrak{x}\_{i\prime}\mathfrak{x}\_{j}) = (\mathfrak{x}\_{i}\cdot\mathfrak{x}\_{j})^{q} \tag{10}$$

where *q* is the degree of the polynom that will be used [62].

• Sigmoid kernel function: this function gives an explanation to practical viability [63]. This function is expressed as follows:

$$K(\mathbf{x}\_{i\prime}\mathbf{x}\_{j}) = \tanh(v(\mathbf{x}\cdot\mathbf{x}\_{i}) + c) \tag{11}$$

where *v* and *c* are adjustable kernel functions based on the data.

**Figure 4.** Mapping of the kernel function.

### *4.2. Search Optimization Algorithms*

In general, SVM models and ANN models are widely used in forecasting power consumption and solar radiation because of its performance and easy adaptation for nonlinear variables. However, it is important to improve new methods to find the best fitting parameters of the model [64].

The importance of using genetic algorithms as a solution to the Vapnik model function (Equation (1)) is to find the most optimal point to be able to obtain better prediction results. Intelligent search algorithms (SOA) are techniques for searching and exchanging information between the individuals of a population, which will be the ones that can solve a non-linear optimization problem [65]. SOA models are frequently used in forecasting methods where it is quite crucial to determine weights coefficients to reduce the error [66].

### *4.3. Performance Evaluation*

It is necessary to know the performance of the SVM model. To evaluate the performance of forecasting solar radiation data, statistical indicators were used [67]. According to the literature, there are five popular indicators that could determine the accuracy of the predicted data [68,69]:

• Mean absolute percentage error (MAPE): It is used to express the absolute error of the predicted and observed variables in percentage [70]:

$$MAPE = \frac{1}{m} \sum\_{i=1}^{m} \left| \frac{y\_i - \mathcal{Y}\_i}{y\_i} \right| \tag{12}$$

• Root mean square error (*RMSE*) sizes the goodness of the fit related to forecast with high errors [71]:

$$RMSE = \sqrt{\frac{1}{m} \sum\_{i=1}^{m} (y\_i - \mathcal{y}\_i)^2} \tag{13}$$

• Mean bias error (*MBE*) indicates the deviation of predicted data from the observed data to provide the performance information of short and long-term use of the model [72]:

$$MBE = \frac{1}{m} \sum\_{i=1}^{m} (y\_i - \mathcal{Y}\_i) \tag{14}$$

• Mean absolute error (*MAE*) indicator gives a perspective of the performance of the prediction model by viewing how close the predicted variables are to observed variables [73]:

$$MAE = \frac{1}{m} \sum\_{i=1}^{m} |y\_i - \mathcal{Y}\_i| \tag{15}$$

• Relative root mean square error (*RRMSE*) quantifies the relative spread in the error [74,75]

$$RRMSE = \frac{1}{\mathcal{Y}} \sqrt{\sum\_{i=1}^{m} \frac{(\mathcal{y}\_i - \mathcal{y}\_i)^2}{N}} \tag{16}$$

where *yi* is the global solar radiation measured, *y*ˆ*i* is the predicted global solar radiation, *y*¯ is the mean global solar radiation, *m* is the number of forecast data points and *N* is the number of validation data. These indicators make it possible to know the efficiency of the SVM models: If the statistical indicator value is zero in the ideal case and presents a good performance if is closer to zero [8,57,67,76].

### **5. Results in Hybrid Techniques: Support Vector Machine with Search Algorithms Review**

This work compiled and analyzed scientific articles focused on the construction of solar radiation prediction models with a hybrid method based on support vector machine and search optimization algoritms. The articles that are studied show the advances in this field and the solutions proposed to obtain predictive radiation data with the minimum error calculated with the statistical indicators.

SVMs relate to regularization networks and offers an advance on the ANN model. It is based on the theory of statistical learning that adopts least squares methods to solve the problem to least square solutions through a set of linear equations based on the minimization of structural risk.

Therefore, the SVM model can avoid excessive adjustment of the training data, does not require an iterative adjustment of the model parameters, has better generalization, requires few cores and has good performance [77]. To determine the appropriate choice of the prediction model with a support vector machine, the climatic variables of maximum and minimum air temperature (*Tmax* and *Tmin*, respectively), maximum and minimum relative humidity (*Hmax* and *Hmin*, respectively), wind speed (*ws*), evaporation (*E*) and vapor pressure estimates (*V*) were the most used for future global solar radiation *Gh* [78–81]. This model can avoid the excessive adjustment of the training data, does not require an iterative adjustment of the model parameters, has better generalization, requires few cores and has good performance [55]. Figure 5 visualizes the general process to predict global solar radiation using the climatic variables.

**Figure 5.** Training process using climatic variables for SVM models.

Different researchers have performed models based on SVM with no optimization models, and the results were prominent. Ref. [76] presented a hybrid SVM model in predicting solar radiation on two sites for PV panel surfaces using diffuse, direct and global solar radiation as input variables. The results show that the SVM model improves the prediction compared to the ANNs in the training phase with a maximum time of 0.0468 s, almost 5 times faster than an RNA model and 2 times faster in the testing phase, with 2.15 s. The stability of the model in this application lies in the use of these three radiation components with an RMSE three times smaller than an ANN, varying in a range from 18.34 to 31.15.

Authors from [36] proposed a solar prediction model based on SVM for one-hour ahead based on the typical climatic variables. They mentioned that SVM regression significantly improves the prediction accuracy according to the statistical indicators.

Authors from [67] proposed multiple SVM models using temperature equations with different kernel equations using only temperature variables from meteorological stations. The researchers observed that creating several empirical temperature-based equations and evaluating the models by several statistical indices, the prediction improves significantly in its performance.

In [82], a forecasting model based on SVM and using satellite images of clouds as the input space was presented. They used 4-year registers of cloud monitoring systems configured to perform the model by using large-scale data in multiple inputs and outputs. The performance SVM model was compared to other predicting models, providing a robust and great alternative to predict solar energy with *RMSE* 10.86.

In [83], a model was developed to estimate global solar radiation based on a support vector machine (SVM) using an index that indicates air quality as an input variable to evaluate the performance of this technique. Compared to existing models, such as neural networks, the model presents a great performance in the accuracy of global solar radiation models by using an air quality indicator as an additional input parameter. The related works only used SVM models. Table 1 summarizes the findings about using SVM as an important alternative to ANN or other models to predict solar radiation.

**Table 1.** Techniques of solar radiation prediction using SVM.

In recent years, the state-of-the-art with respect to SVM has been demanding due to its high performance, so various methods of improving these models have been developed, such as integrating a search algorithm to the SVM model. Some works have proposed various hybrid techniques based on the flow diagram in the Figure 6.

**Figure 6.** Flow chart of SVM model with search algorithm.

The work of [14] showed an SVM with the Firefly Algorithm (FFA) to estimate future solar irradiance. A Firefly Algorithm (FFA) was used to calculate the best kernel parameters of the SVM model to obtain the most accurate prediction data. The model presents a better performance compared to ANN, which has RMSE = 1.86% and a MAPE value of 11.51%. The results are promising for the use of FFA with SVM as a great alternative to predicting the global component of the solar radiation. Additionally, they specified a risk in zones with no constant climatic conditions.

In [84], a new approach was developed using data mining techniques to model and improve the prediction of hourly global solar radiation the next day. They report an RMSE value of 23.5% when the independent variables were the values of the meteorological parameters of the previous day and an RMSE of 22.9% when the independent variables were the forecasts of meteorological variables for the same day to be forecast.

The authors [85] developed an SVM model optimized by glowworm swarm optimization (GSO) to estimate solar radiation. Employing the Eclat algorithm (data mining) to choose the appropriate predictors, the support vector machine (SVM) with penalty function of the model structure and the GSO algorithm to improve forecast accuracy, selecting the appropriate parameters. The tests were carried out in four areas of the USA using the input variables of average daily airmass, dew point, relative humidity, opaque cloud cover, wind speed, vertical wind shear, pressure, albedo, zenith angle, azimuth angle, net radiation, global normal radiation and global extraterrestrial radiation. The results show a great performance—RMSE = 0.43 W m2 and MAPE 5.64%.

In [86], an SVM-FFA model for the prediction of the monthly mean daily horizontal global solar radiation in the port of Bandar Abbass was proposed, located in the southern coastal region of Iran. They also show that their model is highly efficient in estimating the monthly mean daily horizontal global solar radiation. With the following statistical indicators, MAPE = 3.3252.

Ref. [87] presented an SVM model with optimization using GA to compare its performance with a grid search evaluating the accuracy by the analysis of the parameters of the four kernel functions. They found that the grid search could be a good alternative only if a low dimensional dataset, which is only present in few cases of classifications problems, but GA is the best alternative in more cases, presenting stability above 15.9-times that of a grid search.

Research from [57] presented a study of three models: adaptive neuro-fuzzy inference system (ANFIS), ANN and (SVM) models using GA in the three cases to estimate the kernel model parameters in order to compare their accuracy and performance in daily global solar radiation. The results show a huge advantage using SVM-GA to predict solar radiation. They found that SVM-GA also determined the best kernel parameters by evaluating the error to the global minimum convergence.

In [55], a wavelet-coupled vector support machine model for predicting global solar incidents using a limited meteorological data set for the city of Brisbane, Australia was performed. The data were decomposed into a subset of wavelets, transforming the input space into discrete variables to create new time series using the Daubechies-2 wavelet to a detailed level. This hybrid model obtained an approximation of *R* = 0.965. The input variables used in this model are the basic climatic variables and add hours of sunshine *St* to predict the daily global incident solar radiation (*Rn*). They also mentioned that climatic anomalies could make the W-SVM model less accurate for solar radiation forecasting.

The authors from [88] presented a GA-SVM model to predict the short-term power forecasting of a PV system on a residential scale. GA was used to find the optimal parameters values for a kernel within a base classifier of the SVM. They used climatic variables and added the solar radiation to predict the power load of a photovoltaic system, that shows that SVM is multi-configurable according to the problem. The results show that GA-SVM outperforms the conventional SVM. The GA shows its great performance by finding a local minimum that is translated in a high-accuracy solar radiation forecasting. Table 2 summarizes the findings on the use of SVM optimized with search algorithm techniques.

Ref. [89] performed a comparison between SVM models and copula-based nonlinear quaintly regression (CNQR) to predict dayli diffuse solar radiation and mentioned the principal risk working with SVM models: There is a possibility that SVM could not perform an accurate prediction if lower parameters are considered in the space data; with more input parameters, the accuracy will increase but the parameter optimization increases considerable, leading to a high cost for the prediction.

The authors from [90], developed three hybrid models to predict the daily solar radiation in urban zones: SVM-PSO model, Bat algorithm with SVM and Whale optimization algorithm with SVM using the general climatic data and adding the Ozone (*O*3) in the space variables. They found that the *O*3 variable improves the accuracy of the predicted data compared to other parameters such as air pollution, using the RMSD statistical indicator to measure the performance, with values of 11.1%, 10.0% and 10.4%, respectively, for the three developed models.

**Table 2.** Techniques of solar radiation prediction using SVM and search algorithms.

### **6. Discussion and Conclusions**

This article presented a state-of-the-art in solar radiation prediction techniques through the use of support vector machine models (SVM) and the optimization of the search for the parameters that guarantee a favorable performance in the accuracy of the forecasting data.

It was observed that SVMs by themselves show a better performance in predicting solar radiation than artificial neural networks and other prediction models such as auto regression. In this context, it is possible to define these models as premature techniques, which today are consolidated as the best predictive models, since these models emerged 6 years ago [8]. The SVM models show an improvement when evaluating the polynomial kernel functions using only as a basis the climatic variables of maximum and minimum air temperature (*Tmax* and *Tmin*, respectively), maximum and minimum relative humidity (*Hmax* and *Hmin*, respectively), wind speed (*ws*), evaporation (*E*) and vapor pressure estimates (*V*). That base of the input space is sufficient to perform a prediction model based on SVM. However, it is possible that the polynomial function demands a high cost in performance because it can increase in degree, so another alternative is the use of a radial basis function (RBF). According to the literature, RBF to estimate the kernel parameters can improve the accuracy up to 7% compared to higher degree polynomial functions. Likewise, the most widely used validation method is *RMSE*, since the observed error contributes to the tuning of the weights of the SVM model. Additionally, according to [76], SVM could be 5 times faster than ANN in the training phase and 2 times faster in the testing phase.

The execution time of each optimization algorithm will depend on the number of input variables to the prediction model, the number of data and the prediction horizon. In order to obtain a more robust model, it must be calibrated with the variables and geographical conditions, since they change according to the area, so the accuracy totally depends on the location. These search algorithms considerably revolutionize the performance of the SVM model compared to a simple SVM model, contributing to obtaining a lower prediction error in the output of the solar radiation data. In addition, in these models, the RBF kernel function presents better performance in the search for parameters due to the inclusion of these search algorithms. These algorithms can be evaluated based on MAPE, RMSE and MAE as the main statistical indicators to evaluate the predicted data of solar radiation that measure the performance of the model.

Finally, the versatility of the model using only climate data for the prediction of solar systems is noteworthy, such that in the future they may be a beneficial alternative in the sizing, generation and management of alternative energy to contribute to the energy transition that in the coming years will be one of the main issues that researchers will tackle to solve problems.

**Author Contributions:** Conceptualization, J.M.Á.-A. and J.G.R.-M.; investigation, J.M.Á.-A.; validation, G.R.-L., E.V.-R.J. and S.A.O.-B.; writing—original draft preparation, J.M.Á.-A.; discussion, J.G.R.-M.; writing—review and editing, M.T.-P. and J.G.R.-M.; formal analysis, J.G.R.-M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Abbreviations**

The following abbreviations are used in this manuscript:

### **References**

*Article*

### **Prediction of Building's Thermal Performance Using LSTM and MLP Neural Networks**

**Miguel Martínez Comesaña 1,\* , Lara Febrero-Garrido 2 , Francisco Troncoso-Pastoriza 1 and Javier Martínez-Torres 3**

Received: 22 September 2020; Accepted: 20 October 2020; Published: 23 October 2020

**Abstract:** Accurate prediction of building indoor temperatures and thermal demand is of great help to control and optimize the energy performance of a building. However, building thermal inertia and lag lead to complex nonlinear systems is difficult to model. In this context, the application of artificial neural networks (ANNs) in buildings has grown considerably in recent years. The aim of this work is to study the thermal inertia of a building by developing an innovative methodology using multi-layered perceptron (MLP) and long short-term memory (LSTM) neural networks. This approach was applied to a public library building located in the north of Spain. A comparison between the prediction errors according to the number of time lags introduced in the models has been carried out. Moreover, the accuracy of the models was measured using the CV(RMSE) as advised by AHSRAE. The main novelty of this work lies in the analysis of the building inertia, through machine learning algorithms, observing the information provided by the input of time lags in the models. The results of the study prove that the best models are those that consider the thermal lag. Errors below 15% for thermal demand and below 2% for indoor temperatures were achieved with the proposed methodology.

**Keywords:** neural network; LSTM; MLP; thermal inertia; building performance

### **1. Introduction**

The residential sector makes an important contribution to energy consumption worldwide, representing more than the 40% of the total energy use in the European Union (EU) [1]. The EU has established several guidelines and directives to improve energy performance in buildings, such as 2010/31/EU (EPBD) [2] and 2012/27/EU [3], requiring that new buildings comply with nearly zero-energy buildings (NZEB) by 2030 [4] and to reach a decarbonized and highly energy efficient building stock by 2050. Therefore, energy efficiency in buildings is of great importance to the overall sustainability. Knowing exactly and precisely the energy consumption of a building is the first step to be able to optimize its energy performance. However, forecasting the building energy consumption is a difficult issue that many authors have investigated in recent years [5–7].

Different methods have been developed to predict building energy demand. Traditionally, dynamic simulation has been used successfully [8–10] becoming a suitable tool, that enables the assessment of building performance and the calculation of energy use. There are several building energy performance simulation (BEPS) tools available, such as TRNSYS [11], EnergyPlus [12] or DOE-2 [13]. However, the use of dynamic simulation tools requires the knowledge and control of many different parameters of the building, for example, envelope, materials properties, lighting, equipment, heating ventilation and air-conditioning (HVAC) systems and the behavior of the users, which makes the work of data collection arduous and increases the difficulty of obtaining exact results. The validation of these building energy simulation models is usually carried out with standard criteria such as the coefficient of variation of root mean square error (CV(RMSE)). Currently, these models are considered calibrated if they meet the criteria established by American Society of Heating Refrigerating and Air-Conditioning Engineers (ASHRAE) Guideline 14 [14], which states that a model is calibrated if its CV(RMSE) is below 15% monthly or 30% hourly. On the other hand, new tendencies in building prediction have emerged. Examples of this are the black box models, which stand out for their ability to train based on available data [15]. In this context, several data-driven models for building energy forecasting have been used more and more in recent years due to their robustness, resilience, strong yield and ease of implementation. Amongst them, the most popular data-driven approaches are the ones based on artificial neural networks (ANNs) [16,17].

ANNs are specific mathematical models that attempt to replicate the way a human neural network proceeds. The main and most important feature of these models is the ability to learn; that is, from known data, they are able to extract a pattern and then extrapolate the results to new data. Thus, neural networks have a remarkable ability to model the non-linear relationships between inputs and outputs, and for their massive interconnectivity [10,18,19]. Two of the most used ANNs in the context of building energy prediction are multi-layered perceptron (MLP) [19,20] and longshort-term memory (LSTM) neural networks [5,21]. MLP models are well known for being composed of several layers connecting the inputs to the specific output [10]. They have been applied to numerous scientific fields such as environment [22–24], econometrics [19,25], medical research [26], chemistry [27,28] and even building energy [10,29,30]. On the other hand, LSTMs are more current models, but their application has also expanded to a great number of fields. This type of ANN is characterized by a recurrent neural network (RNN) architecture [21,31]. RNNs differ from traditional feed-forward neural networks in that they have a hidden layer to consider connections with the previous values. This enables this type of neural network to model long-term dependencies; such as thermal inertia in buildings [32–34]. They encompass cyclic connections that make them, in principle, more suitable for modeling time sequence data than forward neural networks [35,36]. Some of the numerous fields where LSTM neural networks were applied are language modeling [36,37], speech recognition [35], tourism flow [38,39], sequential diagnosis [32,40] and energy analysis in buildings [21,41,42].

The aim of this work is to study the thermal inertia of a building by developing a methodology using MLP and LSTM neural networks. This procedure was applied to a public library building located in the north of Spain [43,44]. The available data are hourly observations of thermal demand and indoor temperatures of the building. Moreover, hourly observations of two weather variables (outdoor temperatures and solar radiation) and three time variables (hour of the day, day of the week and hour of the year) are also needed. A comparison between the prediction errors according to the number of time lags (each lag is one hour) introduced into the models, with respect to the studied variables, has been carried out. In addition, the accuracy of the models was measured using the CV(RMSE) as advised by AHSRAE. The novelty of this work lies in the analysis of the building inertia through two different artificial intelligence algorithms observing the information provided by the introduction of more time lags in the models.

### **2. Materials and Methods**

### *2.1. Artificial Neural Networks Developed*

The different mathematical models that have been built to carry out the analysis are presented in this section. The models are of two different types of ANNs: An MLP and an LSTM neural network.

### 2.1.1. Multi-Layered Perceptron Neural Networks

In this work the structure of the MLP neural network is variable due to the use of different numbers of time lags. Depending on the structure of these artificial neural networks, different network models can be generated. The most commonly used is the so-called feed-forward model. This model is composed of several layers (see Figure 1). The first layer, or input layer, is where the model inputs are introduced, whereas the last layer, or output layer, is where the results of the trained network are given. In addition, among them, the number of intermediate layers (hidden layers) can be zero, one or more [18,20]. The particular network architectures are differentiated by the number of hidden layers and hidden neurons according to the complexity of the problem [45,46]. The grid of tested values for the hidden layers was between 0 and 4. In contrast, the grid of values for the hidden neurons was composed of specific numbers based on the formulas presented by Shin-ike [47], Doukim et al. [48] or Vujicic et al. [49], which take into account the size of the sample and the number of inputs and outputs of the network. Furthermore, in this case, due to the complexity of the problem, this grid of neurons had to be expanded with the values: 50, 100, 200, 500 and 1000.

**Figure 1.** Multi-layered perceptron (MLP) neural network architecture with *N* inputs, *g* hidden layers and univariate output.

The MLP neural network training is done with a backward propagation algorithm; errors are propagated through the network and adapted to the hidden layers. Thus, error-correction learning (actual system responses must be known) is used to train the ANN [29,50]. There are several ways to perform the training [19,51] but, the methodology used in this work consists of updating the weights with an average update of the weights (batch learning), which is achieved by incorporating all the patterns in the input file (an epoch) and accumulating all the weight updates. There is also a need for a stop criterion [52,53]. Although there are other options, such as a threshold for the mean square error (MSE) or a limitation on the maximum number of iterations, the most widely used is cross-validation. This method is most effective in stopping training when the best generalization is achieved. It consists of separating a small part of the training data and using it to evaluate the trained network. Thus, training should stop the moment the network performance, quantified by the mean square error (MSE), begins to decrease or stagnates [18–20,31].

Lastly, in all MLP models created in this paper the activation function selected was the reLU function (rectified linear unit: max(0, *x*)) [54], the kernel initializer was normal and the optimizing algorithm was the adaptive moment estimation (Adam) [55]. On the other hand, the batch size used to train the MLP neural networks was equal to 64 (mini-batch gradient descent [56]) and the patience (limit to stop the training if the performance of the model does not improve) was 100 epochs.

### 2.1.2. Long Short-Term Memory Neural Networks

LSTM is a recurrent neural network architecture specifically for modeling time sequences and their long-range dependencies more accurately than conventional RNNs. LSTM introduces specific units (memory blocks) into the recurrent hidden layer that store the temporal state of the network through self-connections. In addition, the memory block also contains units known as gates to control the flow of information [38,40]. A forget gate was incorporated into the memory blocks to prevent the LSTM model from processing continuous input flows without being segmented in subsequences. Furthermore, more complex LSTM networks also incorporate peephole connections between internal cells and their gates to learn the timing of the outputs [32,35].

The structure of traditional RNNs [57] can be presented by deterministic transitions from previous to current hidden states in form of a function:

$$h\_t^{l-1} \; \; h\_{t-1}^{l} \; \rightarrow \; \; h\_t^{l}$$

where *hl t* represents a hidden state in layer *l* in timestep *t*.

However, LSTM counts with a complex dynamic that allows memorizing information for many timesteps. Long-term information is stored in a vector of memory cells *ml t* ∈ R*n* [21,37]. The LSTM networks are empowered to decide whether to overwrite, retrieve or maintain this information for the next timestep. Their architecture is the following:

$$h\_t^{l-1} \text{ \ } h\_{t-1'}^l \text{ \ } m\_{t-1}^l \to h\_{t'}^l \text{ \ } m\_t^l$$

Let *x*1, *x*2, ... , *xk* be a classic input sequence (with *k* lags), where *xt* ∈ R*D* denotes the vector of *D* real values at the timestep *t*. The LSTM architecture defines an input gate *it*, a forget gate *ft* and an output gate *ot*. At specific time *t*, as is observable in Figure 2a where an LSTM cell is shown, the hidden layer output would be *ht*, the cell input state ∼ *Nt* and the cell output state *Nt* [21,38,39,41].

**Figure 2.** (**a**) Long short-term memory (LSTM) cell structure and (**b**) forecasting framework of LSTM neural network architecture.

The two values *Nt* (Equation (1)) and *ht* (Equation (2)) that are transmitted to the next timestep, following the steps explained in [21,39], can be calculated as follows:

$$\mathbf{N\_{t}} = \mathbf{i\_{l}} \times \mathbf{\tilde{N}\_{l}} + \mathbf{f\_{l}} \times \mathbf{N\_{t-1}} \tag{1}$$

$$h\_t = \ o\_t \times \tanh(\mathcal{N}\_t) \tag{2}$$

A particular LSTM neural network was built for the prediction of indoor temperatures and thermal demand of the studied building (see Figure 2b). In this case, the input matrix of the network would be the observed data *Xt* (see Table 1) and the one-dimensional output the predicted future data ∼ *Yt*+1 (Equation (3)). Once known *ht* (Equation (2)) the network output will be:

$$
\tilde{Y}\_{l+1} = \tilde{W}\_2 \times h\_l + b\_l \tag{3}
$$

where *W*2 is the weight matrix connecting the last hidden layer to the output layer and *b* the bias term in the output layer (see Figure 2b).

**Table 1.** Summary of the data set, with the available variables, used in this study.

Specifically, the LSTM architecture used in this study is also variable due to the different complexity of the models based on the number of time lags introduced. The variable parameters were the number of hidden and LSTM layers and the number of hidden neurons. In this case, the layer structures (LSTM and hidden layers) in the tested grid were: 2-1, 2-2, 3-1 and 3-2. On the other hand, the number of hidden neurons was a grid with eight values between 5 and 500. In addition, as in the MLP models, the activation function of the models was the reLU function [54], the kernel initializer was normal, the optimizing algorithm was Adam [55] and the batch size equal to 64 [56]. In order to avoid overfitting problems, they were trained with the same early stop as MLP neural networks (100 epochs) [52,53].

### 2.1.3. Pre-Processing Data

The data available in this study are hourly observations of two variables related to climate (outdoor temperature and solar radiation) and two related to the energy performance of the building (indoor temperature and thermal demand) between February and March of 2017 (see Table 1).

A continuous sample without large sets of missing values was necessary due to the use of hourly time lags as inputs in the models. This means that to train a model, besides the information of the independent variables at a given moment, the information of past instants of certain variables is also introduced. This study focuses on the prediction of thermal demand and indoor temperature of the studied building. The variables that have been used as explanatory variables of the models are presented in Figure 3. Although only the variables related to thermal conditions have been lagged, to add more information to the models, three time variables have been added (see Table 1). These three variables (hour of day, day of week and hour of year) do not provide any additional information being delayed due to their artificial nature.

**Figure 3.** The variables used as input into the models to predict the thermal demand of the building (left) and the indoor temperatures of the building (right).

Generally, as the number of time lags considered increases, the complexity of the model needs to change (there is more information to be extracted). Therefore, a process based on repetition (cross validation design; data are divided into training, test and validation [58]) was carried out to find the best parameters for each model and for each of the different lagged architectures studied. This process consists of repeating a prediction of the same test sample 10 times with each model (the second week of January 2019), training with the same data, to extract the average performance of each of them and being able to compare them. Thus, the best architecture for each structure of lags and to each of the mathematical models used (LSTM or MLP) is obtained.

After finding the best model for each lag architecture analyzed (see Tables 2 and 3), the thermal inertia of the building has been tested with two new weeks of 2019 (validation data):

### *2.2. Data Acquisition of the Experimental Case Study*

The process of data acquisition, as well as the characteristics of the energy system of the building studied, is described in this section.

### 2.2.1. Building and HVAC System Description

The case study building where this new methodology has been proved is a public library located in the city of Vigo, in the north west of Spain. It is a three floor building with connected and large open areas that trigger temperature stratification inside the library. The building has a working floor area of 820 m2 and the average capacity is 220 persons. This total usable floor area includes 412.3 m2 on the ground floor, 233.5 m2 on the first floor and 73.0 m2 on the second floor. Due to the interconnected spaces and the large window area, the building experiences high thermal inertia.

The HVAC system is a ground-source heat pump (GSHP) distributed with a radiant floor (RF). The reversible heat pump is a CIATESA IZA 185 with a nominal heating capacity of 45 kW and nominal cooling capacity of 35.7 kW, with a COP of 3.63 for heating and EER of 3.25 for cooling. It has a glycol-water mixture flow rate of 8450 m3/h. The water storage tank contains 200 L. The GSHP is composed of six boreholes 100 m deep. The set point temperature for heating is 20 ◦C when the library is open and 17 ◦C when the library is closed. For cooling, the setpoint temperature is 26 ◦C when the library is open and 29 ◦C when closed. The tested heating period goes from October to May.

Further information about the building and its HVAC system can be consulted in [43,44,59].

### 2.2.2. Data Acquisition System

The data acquisition (DAQ) system in the building is further described in other articles by Cacabelos et al. [43,59] and Fernández et al. [44]. The rate of this system is one minute.

There are several wall module temperature sensors; six of them placed in the ground floor, four of them in the first floor and one in the second floor, all of them from the manufacturer Honeywell. In addition, the thermal zone temperature is calculated as the average of temperatures of every sensor from each floor. On the other hand, thermal energy consumption is measured with a thermal energy meter Multical 601 manufactured by Kamstrup.

The meteorological data were gathered from a weather station located at 42.17◦ N latitude and 8.68◦ E longitude at an altitude of 460 m. This station belongs to the weather station network from the Environment, Land Planning and Infrastructures Department, and it is located only 500 m from the building.

Further information about the DAQ and the weather station can be found in [43,44,59].

### *2.3. Validation and Error Measurement*

The error measure considered in this paper to quantify the accuracy of each of the proposed models is the CV(RMSE) (coefficient of variation of the root mean square error):

$$CV(RMSE) = \frac{1}{\overline{Y}} \sqrt{\frac{\sum\_{i=1}^{N} \left(Y\_i - \hat{Y}\_i\right)^2}{N}} \tag{4}$$

This measure was used both to find the best structure depending on the number of time lags introduced in the models and, also, to compare the different models with the test samples (Tables 2 and 3). It has been used in similar studies such as that of Hong et al. [60] and Kuo et al. [61].

#### **3. Results and Discussion**

The inertia in the thermal conditions of the Science Library of the University of Vigo was analyzed through the variation of the CV(RMSE) (Equation (4)) in the predictions of thermal demand and indoor temperatures of the building. In this analysis the models were trained with two months of hourly observations considering different numbers of time lags (values of certain variables in past hours) with respect to the variables of study to make the predictions. In this way, two specific weeks in 2019 (sample 1 and sample 2) have been taken into account to analyze the thermal inertia of the building by comparing the error obtained by each model according to the number of time lags used.

Section 3.1 presents the results to the thermal demand predictions and Section 3.2 the same analysis for indoor temperature predictions. In each section the evolution of the prediction error was analyzed based on the number of hourly time lags introduced in the models. The numerical results are the mean and the standard deviation of the CV(RMSE) obtained in 10 repetitions for each model and sample analyzed. All the specific architectures used in the models are shown in Tables 2 and 3. Lastly, all figures presented in this section were made with the Python programming language [62].

### *3.1. Thermal Demand Analysis*

The results of the thermal demand predictions for the two samples studied, through which the thermal inertia of the building was analyzed in this section, are shown in Figure 4. Specifically, the results of the LSTM models are shown in Figure 4a and the results of the MLP models are shown in Figure 4b. Each time-lag structure analyzed is represented in the figures through the prediction with the lowest error among the 10 repetitions of the analysis.

LSTM models were capable of replicating the trend of the thermal demand of the building (except for certain peaks). They also yielded a small variation between the different replications of the analysis. This is shown in Table 2, with average errors of less than 20% and standard deviations below 0.05, in the two samples considered in this study. In sample 1 (first row of Figure 4a) the lag structure that produced the best results is the one with 6 lags (CV(RMSE) = 15.26%). In addition, it is also the structure with the smallest variation in the 10 experiment repetitions. As shown in Table 2 and Figure 4a, in this sample the average errors are, in general, higher than those of sample 2. Although in Figure 4a all the models are close to the real values, as only the best predictions of each model are shown, in Table 2 the differences between their average performances are presented. Thus, the model that considers 12 lags is the architecture that obtained the highest average error (CV(RMSE) = 19.13%) and presents the largest variability in its results. In sample 2 (second row of Figure 4a), where the models were better adjusted to reality, the architecture that obtained the best average performance, as shown in Table 2, is the 1-lag architecture (CV(RMSE) = 14.35%). In this case, smaller errors are not accompanied by more stable results; the highest stability was achieved with the model with 24 lags. On the contrary, the model that obtained the highest error and the largest variability is the one that considers 12 lags (CV(RMSE) = 17.34%).

**Figure 4.** Results of thermal demand predictions in sample 1 (upper side) and sample 2 (bottom side) considering the different structures of time lags. Each architecture is represented by the prediction curve with the lowest error among the 10 repetitions: (**a**) LSTM models; (**b**) MLP models.

MLP models, as shown in Figure 4b, were also capable of predicting the actual thermal demand of the building with errors of less than 18%. As with LSTM models, certain peaks in the validation sample were not perfectly predicted and their variability among analysis repetitions was also small (below 0.02). However, these models yield better results than LSTM models (see Table 2). In sample 1 (first row of Figure 4b) the model with the lowest average error, as with LSTM models, is the model that considered 6 time lags (CV(RMSE) = 14.90%). In addition, this architecture is the one that produced the most stable results (see Table 2). On the contrary, the lag structure that adjusts reality in the worst way and produced the highest average error is the 12-lag structure (CV(RMSE) = 19.38%). It is also the structure that yields the most variable errors. In this case, a higher average error means more variability in predictions. In sample 2 (second row of Figure 4b), as in LSTM models, the predictions, on average, were better than in sample 1 (see Table 2). The architecture of lags that produced the lowest

average error is the architecture of 1 lag (CV(RMSE) = 14.17%); the same structure as with the LSTM models. In this situation it is not the most stable model. On the other hand, the model that yielded the worst performance is the model that considered 24 lags (CV(RMSE) = 17.43%); not being the most variable model. In the case of sample 2, as presented in Table 2, the models with the least errors and the most stable models do not coincide.

**Table 2.** Results of the thermal demand predictions based on the time lags introduced in the neural network model. LSTM errors are presented in the first three columns, whereas MLP errors are presented in the last three. The results shown summarize the 10 repetitions of the experiment, for each architecture and each sample analyzed. The mean of the CV(RMSE) obtained in each of them and the standard deviation (SD) are presented. The different structures of the neural networks that obtain the best accuracy are presented between brackets.

In the case of the thermal demand of the building studied, its thermal inertia is significant but up to a certain limit of time lags. Both LSTM and MLP models show that the best fit to reality is obtained from an architecture with 6 lags in sample 1 and 1 lag in sample 2 (see Table 2). Therefore, the inertia influence is greater with the data of sample 1. Moreover, the LSTM models, for both samples, present the worst performance with 12 lags; but the result, although not surpassing the best one, improved by introducing 24 lags. On the other hand, with MLP models the introduction of more time lags than the optimal ones cause an increase in the average error. MLP models obtained the worst time-lag structure in both cases with the 24-lag structure.  Regarding the stability of the results obtained, all the models show a small variability and, although there is not much difference between them, MLP models show smaller standard deviations. In addition, the model without time lags in any case obtained the best performance. This means that a thermal inertia exists and that the introduction of time lags in the thermal input variables provides valuable information.

### *3.2. Temperatures Analysis*

The analysis of the thermal inertia of the building, for the two samples studied, through the predictions of the indoor temperatures of the building studied, is shown in Figure 5. Specifically, the results of the LSTM models are shown in Figure 5a and the results of the MLP models are shown in Figure 5b. As in the previous section, each time-lag structure studied is represented in the figures through the prediction with the lowest error among the 10 repetitions of the analysis.

LSTM models optimally predict the indoor temperatures of the building analyzed in this study. Table 3 shows that average errors (measured by the CV(RMSE)) are below 3% and their variability is small (below 0.005); both for sample 1 and for sample 2. In sample 1 (first row of Figure 5a), contrary to the previous section, the average errors are lower than those of sample 2. Furthermore, the lowest average error of this sample is obtained with the 1-lag structure (CV(RMSE) = 1.89%). In addition, although the models results are stable, this specific architecture shows the smallest standard deviation among all the structures analyzed (see Table 3). In this case, in Figure 5a it is observed that the two architectures with 12 and 24 lags are straight lines. This fact means that the introduction of more time lags does not provide useful information. Although the average error obtained by these two architectures is small, the models do not efficiently adjust to reality. Thus, the architecture that obtains the highest average error is the architecture with 24 lags (CV(RMSE) = 2.68%). Moreover, it is the least stable model structure in relation to the variability of the errors obtained (see Table 3). There is a correlation between the variability of the errors obtained and the average error; the greater the

variability, the greater the average error. On the other hand, in sample 2 (second row of Figure 5a) the 1-lag structure, as in sample 1, obtained the lowest average error (CV(RMSE) = 2.05%) and the smallest variability (see Table 3). As shown in Figure 5(a), with this sample the different models analyzed show a more similar behavior among themselves than in sample 1. Furthermore, unlike the first sample, the introduction of more time lags in the models eventually reduces the average errors; the best models, are those with 1 and 24 time lags (see Table 3). Therefore, the architecture with the highest average error is the 6-lag architecture (CV(RMSE) = 2.23%). In this case, there is no relation between average error and the variability of the results.

**Figure 5.** Results of the temperatures predictions in sample 1 (upper side) and sample 2 (bottom side) considering the different structures of time lags. Each architecture is represented by the prediction curve with the lowest error among the 10 repetitions: (**a**) LSTM models; (**b**) MLP models.

MLP models, presented in Figure 5b, once again provide better results than LSTM models. As shown in Table 3, their average errors are less than 2.5% and their variability is below 0.003. In sample 1 (first row of Figure 5b), the model structure with the lowest average error is the 24-lag structure (CV(RMSE) = 1.73%). On the contrary, the structure that produces the highest average error is the 6-lag structure (CV(RMSE) = 2.06%). In this case, and similarly to LSTM models, there is no direct relation between average error and variability in the results. In addition, as can be seen in the first row of Figure 5b (and the first row of Figure 5a), none of the models are able to replicate the peaks of the indoor temperatures at the weekend. Nevertheless, as shown in Table 3, the models do adjust the data trend efficiently. On the other hand, in sample 2 (second row of Figure 5b), as in LSTM models, the average errors obtained are higher than those of sample 1 (see Table 3). In this case, the architecture that yielded the lowest average error is the 12-lags architecture (CV(RMSE) = 1.90%). Additionally, it is also the architecture that presents the smallest variability among its results. With this sample, the introduction of time lags did provide valuable information because, as shown in Table 3, the best model structures are those with 12 and 24 lags. However, the lag structure that presents the highest average error was 6-lag structure (CV(RMSE) = 2.28%), slightly worse than the average error

of the model without any lag. This means that, although the introduction of time lags provides useful information, the basic model (regardless of thermal inertia) already provides good results.

**Table 3.** Results of the indoor temperatures predictions based on the time lags introduced in the neural network model. LSTM errors are presented in the first three columns, whereas MLP errors are presented in the last three. The results shown summarize the 10 repetitions of the experiment, for each architecture and each sample analyzed. The mean of the CV(RMSE) obtained in each of them and the standard deviation (SD) are presented. The different structures of the neural networks that provide the best accuracy are presented between brackets.

Regarding the results of the indoor temperatures of the building, it can be observed that, on the one hand, there are less differences between the performance of the models in the different samples than in the case of the thermal demand. Moreover, in comparison with the thermal demand results, the different errors obtained for each of the architectures analyzed during the repetitions are lower and less variable. This demonstrates that indoor temperature values are more constant throughout a year (always around the set point temperature of the building). On the other hand, there are differences between LSTM and MLP models (see Table 3). With LSTM models the thermal inertia is less significant than in the case of thermal demand. In the two samples analyzed the best model was the one with a single time lag. In contrast, with MLP models the thermal inertia is more significant than in the analysis in the previous section. In both samples, the best model structures were a 24-lag structure for sample 1 and a 12-lag structure for sample 2; they use information from a higher number of past hours. In this case, MLP models show better performance in both the best and the worst case. In addition, even with less variability in general concerning indoor temperature predictions, MLP models again produced the lowest standard deviations. Finally, as in the case of thermal demand predictions, models without time lags do not produce the lowest average error in any of the analyses presented. This indicates that the thermal inertia is also significant in the case of indoor temperatures predictions.

### **4. Conclusions**

A new application of machine and deep learning algorithms and a new methodology to analyze the influence of the thermal inertia of a building are presented in this paper. The study was conducted by analyzing monitored data of thermal demand and indoor temperatures of the Science Library of the University of Vigo. The data were completed with weather variables (outdoor temperatures and solar radiation) and three temporal variables (hour of the day, day of the week and hour of the year). The aim of this analysis is to study the importance of the thermal inertia of the building in predicting thermal demand and indoor temperatures. The methodology used in this study is based, on the one hand, on a classic machine learning model (MLP neural network) contrasting the improvement in the prediction accuracy provided by the thermal inertia of the building. On the other hand, as a comparison, a deep learning model (LSTM neural network) was carried out through an analogous analysis.

The research contribution of this paper is the application of mathematical methods to evaluate the influence of the thermal inertia of the building in the prediction errors of the thermal demand and the indoor temperatures of the building itself. The principal limitation of this research is the analysis of the thermal inertia of the building based on the difference in the errors rates of indoor temperatures depending on the number of time lags. The stable behavior of the indoor temperatures by default, produce that, in general, the prediction error is very low. This makes analysis and detection of the additional information that can provide the introduction of time lags in the models more difficult. On the other hand, the presented methods contribute with advantages over the already existing

research in building simulation and thermal inertia analysis of buildings. With traditional building simulations models, both to make predictions and to analyze the thermal inertia of the building, it is necessary to have a deep knowledge about the specific subject and control many different energy parameters. With the methods presented in this paper, known as black box models, it is possible to do all kinds of analyses without any specific knowledge, and also to do it faster than dynamic simulation methods. These models rely on having a significant amount of data to be able to extract information from it and extrapolate it. Additionally, the variables involved in this study are commonly monitored in numerous buildings. Thus, the methodology presented can be applied to evaluate the energy performance or the influence of thermal inertia of different buildings.

First, the results show that the black box models yield average prediction errors below the proposed values for considering calibrated models. The results also illustrate the greater accuracy of MLP neural networks in predicting both thermal demand and indoor temperatures. LSTM neural networks, created specifically to model time sequences, show a poorer performance than MLP models due to the peculiarity of this analysis. In this study the models were trained with only two months of data to try to predict future time sequences of the variable of interest. LSTM models are focused on analyses where predictions are made for time sequences immediately consecutive to the training data (as a time series). In this specific case, where the data are not continuous in time, LSTM models are not the most efficient. In this paper it was demonstrated that, for analyses similar to the one presented here, MLP models are more appropriate; they are more flexible if the available data, on the one hand, are not continuous in time or, on the other hand, are very few. In terms of thermal demand predictions, for both sample 1 and sample 2, the best MLP models are better than the best LSTM models. While the average errors of the MLP models were 14.90% and 14.17%, respectively, the average errors of the LSTM models were 15.26% and 14.35%. The same applies to indoor temperature predictions. The average errors of the best MLP models were 1.73% in sample 1 and 1.90% in sample 2, and those of the best LSTM models were 1.89% and 2.05%, respectively. Regarding the variability of the predictions, the results of the LSTM and MLP models are similar, but the MLP models, in most cases, have slightly smaller standard deviations.

Furthermore, the results of the thermal inertia analysis of the building are different in relation to the different models used. Although all the analyses show that the introduction of specific time lags reduces the prediction error, the influence of the thermal inertia of the building changes with LSTM and MLP models. In the case of thermal demand analysis, both LSTM and MLP models show the same result. In sample 1 the optimal number of time lags introduced in the model was 6 while in sample 2 the optimal structure was a model with 1 time lag. It was demonstrated that the inertia of the independent variables of these models (indoor temperatures, outdoor temperatures and solar radiation) does not extend too much in time. Differently, the results of LSTM and MLP models do not coincide with the indoor temperatures analysis. LSTM models, both of sample 1 and sample 2, have their optimum in 1 time-lag. In contrast, with MLP models the best structures were 24 time lags in sample 1 and 12 time lags in sample 2. Thus, MLP models are able to better extract the information provided by the introduction of more time lags; they are able to exploit the thermal inertia of the building.

From an energy point of view, the conclusion is that the thermal inertia of the building does provide useful information to the machine learning models by introducing lagged variables. This was demonstrated by the decrease, up to a certain limit, of average errors with the introduction of time lags. Therefore, it was also demonstrated that the dimension of the thermal inertia, measured by the maximum number of time lags that provides valuable information, can be detected by the methods presented in this paper. Moreover, the analysis showed that MLP models are more appropriated than LSTM models if the available data set is not large and continuous over time. For the two variables and the two samples analyzed the average error produced by the optimal MLP models are lower than those provide by the optimal LSTM models. Lastly, it was shown that there is more inertia in thermal demand data than in indoor temperatures data; the introduction of thermal demand lags, due to the higher variability of their data, provides more useful information to the models than indoor temperature ones.

**Author Contributions:** Conceptualization, M.M.C. and L.F.-G.; methodology, M.M.C. and J.M.-T.; software, M.M.C. and F.T.-P.; validation, L.F.-G., J.M.-T. and F.T.-P.; formal analysis, M.M.C. and L.F.-G.; investigation, M.M.C. and L.F.-G.; resources, J.M.-T. and L.F.-G.; data curation, F.T.-P.; writing—original draft preparation, M.M.C. and L.F.-G.; writing—review and editing, J.M.-T. and F.T.-P.; visualization, M.M.C. and F.T.-P.; supervision, L.F.-G., J.M.-T. and J.T.-P.; project administration, M.M.C. and L.F.-G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Spanish Government (Science, Innovation and Universities Ministry) under the project RTI2018-096296-B-C21.

**Acknowledgments:** This research was supported by the Spanish Government (Science, Innovation and Universities Ministry) under the project RTI2018-096296-B-C21.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Heat Loss Coefficient Estimation Applied to Existing Buildings through Machine Learning Models**

**Miguel Martínez-Comesaña 1,\* , Lara Febrero-Garrido 2 , Enrique Granada-Álvarez 1 , Javier Martínez-Torres 3 and Sandra Martínez-Mariño 1**

Received: 23 November 2020; Accepted: 12 December 2020; Published: 16 December 2020

**Abstract:** The Heat Loss Coefficient (HLC) characterizes the envelope efficiency of a building under in-use conditions, and it represents one of the main causes of the performance gap between the building design and its real operation. Accurate estimations of the HLC contribute to optimizing the energy consumption of a building. In this context, the application of black-box models in building energy analysis has been consolidated in recent years. The aim of this paper is to estimate the HLC of an existing building through the prediction of building thermal demands using a methodology based on Machine Learning (ML) models. Specifically, three different ML methods are applied to a public library in the northwest of Spain and compared; eXtreme Gradient Boosting (XGBoost), Support Vector Regression (SVR) and Multi-Layer Perceptron (MLP) neural network. Furthermore, the accuracy of the results is measured, on the one hand, using both CV(RMSE) and Normalized Mean Biased Error (NMBE), as advised by AHSRAE, for thermal demand predictions and, on the other, an absolute error for HLC estimations. The main novelty of this paper lies in the estimation of the HLC of a building considering thermal demand predictions reducing the requirement for monitoring. The results show that the most accurate model is capable of estimating the HLC of the building with an absolute error between 4 and 6%.

**Keywords:** energy efficiency; heat loss coefficient; machine learning; XGBoost; MLP; SVR

### **1. Introduction**

Energy efficiency in buildings is a key roadmap to achieve sustainable development worldwide. Nearly one third of the primary energy consumed in the world comes from non-industrial buildings, and it is similar to all the energy used in the transport sector [1]. Consequently, the residential sector constitutes the largest contributor to global warming by means of carbon dioxide emissions. To achieve the goal of reducing building energy consumption, international cooperation is crucial. The European Union (EU) has created a strict legislative framework implementing two important directives: Energy Performance of Buildings Directive 2010/31/EU (EPBD) [2] and the Energy Efficiency Directive 2012/27/EU [3]. Both of them were updated and amended over the years. Specifically, the Directive amending the Energy Performance of Buildings Directive 2018/844/EU [4] includes different aspects strengthening the

commitment of real building stock modernisation through technological improvements [5]. On the other hand, the International Energy Agency (IEA) is also working hard on this issue by means of its Energy in Buildings and Community (EBC) research programme developing activities towards Nearly Zero Energy Buildings (NZEB), the reduction of carbon dioxide emissions and energy savings technologies. One of the high priority research project themes is the building envelope.

The envelope of a building is one of the main sources of heat transfer between the exterior and the interior of the rooms, and it is responsible for the difference between the conditions predicted at the design stage and the actual conditions of the use of the building, including lighting, occupancy and Heating, Ventilation and Air Conditioning (HVAC) systems' operation. This is called the performance gap between building design and operation [6,7]. The envelope efficiency of a building can be characterized under in-use conditions by calculating the Heat Loss Coefficient (HLC) [8]. This coefficient determines the rate of heat flow through the buildings' envelope when a temperature difference exists between the indoor air and the outdoor air under steady state conditions [9]. Therefore, the accurate estimation of the HLC contributes to optimizing the energy consumption of a building. On the other hand, building energy performance simulation tools have been used in recent years to analyse the energy behaviour of buildings [10]. However, achieving an accurate simulation is not trivial. Many sources of error are introduced into the simulation: occupancy and user's behaviour, weather data [11], envelope materials and thickness, the use of electric equipment, etc. The alternative to simulation is the monitoring of the building operation, which leads to the difficulty of the installation of sensors in in-use buildings. Therefore, in this context, the application of advanced mathematical modelling techniques, such as machine learning methods, is becoming more and more common [12].

There are three widely used building energy prediction models: the white-box (physics-based), black-box (data-driven) and grey-box (a combination of physics based and data-driven) modelling approaches [13]. The application of black-box models in building energy analysis has been consolidated in recent years. These models have received considerable attention because of their great success in learning complex patterns and being able to make accurate predictions without specific knowledge of the subject [14,15]. Machine Learning (ML) models, also known as black-box models, are specific mathematical models that are capable of learning a pattern from data and extrapolating it to a new sample. Thus, they were widely applied in studies presenting control strategies to reduce energy costs. Three of the most used are: eXtreme Gradient Boosting (XGBoost), Support Vector Regression (SVR) and Multi-Layer Perceptron (MLP) neural network. Among them, XGBoost is characterized for building regression trees, one-by-one, so that the subsequent models are trained taking into account the residuals of the previous ones [16]. This algorithm has been used in numerous fields such as biology [17,18], econometrics [19], the environment [20,21] or, as in this case, the energy performance of buildings [22,23]. SVR models emerged as a predictive alternative due to the use of a distinctive loss function [24,25], on the one hand, and the dual formulation of the problem [26], on the other. This algorithm focuses on minimising an upper bound of the generalization error instead of minimising the prediction error in the training sample (*empirical risk minimizatio*n) [24]. The usefulness of this algorithm is demonstrated by its expansion into scientific fields such as econometrics [26], health [27,28], electrical efficiency in cities [29], as well as energy consumption in buildings [25,30]. Lastly, MLP neural networks have stood out in past years for their great capacity to model non-linear relationships between certain inputs and outputs, as well as for their massive interconnectivity [31–33]. Furthermore, this type of Artificial Neural Network (ANN) has been applied to different areas of study such as chemistry [34], the environment [35,36], sensors [37] or energy analysis of buildings [38–40].

The aim of this paper is to estimate the HLC of an existing building through the prediction of building thermal demands using a methodology based on ML models. Specifically, three different methods are applied to a public library in the northwest of Spain and compared: XGBoost, SVR and MLP neural network. The dataset consists of hourly observations of two climate variables (outdoor temperatures and

solar radiation) and two variables related to the thermal behaviour of the building (heating demand and indoor temperatures). In addition, to improve model training, three temporal variables (hour of the year, day of the week and hour of the day) are taken into account as model inputs. A comparison between the errors obtained by the different machine learning models is carried out on three different validation samples. Moreover, the accuracy of each of the models is quantified with CV(RMSE) and Normalized Mean Biased Error (NMBE), both recommended by the American Society of Heating, Refrigerating and Air-Conditioning Engineers (ASHRAE), and their usefulness was demonstrated in [41,42]. The novelty of this work lies in the application of machine learning models to estimate the HLC of a building through heating demand predictions. In this context, this paper contributes with a methodology that reduces the exigency of building monitoring, and heating demands will be accurately predicted. Additionally, this methodology does not depend on physical simulations, which require huge amounts of input data and specific knowledge and which have uncertainty in the evaluation of energy renovation or the energy certification of buildings.

### **2. Materials and Methods**

### *2.1. Heat Loss Coefficient Calculation*

The HLC is the most used Key Performance Indicator (KPI) to describe the building envelope energy efficiency. This coefficient reflects transmission heat losses through the envelope per degree difference between indoor and outdoor temperatures from walls, roofs and floors (UA(kW/K)) and ventilation and/or infiltration heat losses per degree difference between indoor and outdoor temperatures (Cv(kW/K))) [43]. Following the method developed by Uriarte et al. [8], the HLC can be described as Equation (1):

$$HLC = \left(ULA + \mathbb{C}\_{\upsilon}\right) \quad \text{(kW/K)}\tag{1}$$

where *UA* stands for the building envelope transmission heat transfer coefficient (kW/K) and *Cv* stands for the infiltration and/or ventilation heat loss coefficient (kW/K). Moreover, the energy balance used to calculate the HLC, assuming stationary conditions, was developed by Uriate et al. [8,44], and it is presented in Equation (2):

$$\sum\_{k=1}^{N} Q\_k + \sum\_{k=1}^{N} K\_k = HLC \sum\_{k=1}^{N} (T\_{in,k} - T\_{out,k}) - \sum\_{k=1}^{N} (S\_a V\_{sol})\_k \tag{2}$$

where *Qk* represents the heating gains (kW), *Kk* represents the internal gains (kW) caused by occupants and electricity consumption (*Kk* = *Kelectricity*,*k* + *Koccupancy*,*k*), (*Tin*,*k* − *Tout*,*k*) summarizes the gap between temperatures inside and outside the building in degrees Kelvin and *SaVsol* stands for the solar gains (kW).

Several terms of this heat exchange are difficult to measure for in-use buildings. For example, the solar gains are complicated to know, and therefore, they introduce an important source of uncertainty into the equation. Thus, in this study, only cold and cloudy periods are considered where solar radiation was very low; direct solar radiation can be considered null; and all the radiation can be perceived as purely diffuse [45]. This is why in this study, direct solar radiation is not taken into account, and only diffuse solar radiation is considered. Another parameter that is difficult to estimate is the occupancy. Therefore, with the aim of reducing the sources of uncertainty, the study is developed during weekends when there is no occupancy. Then, the result of Equation (2) will give an idea of how well built the building is in terms of the enclosure of its envelope.

In this way, the HLC formula that is used in this analysis (average method), taking into account that the periods analysed in this paper were weekends with no occupation and in which there was only purely diffuse radiation, is shown in Equation (3):

$$HLC = \frac{\sum\_{k=1}^{N} Q\_k}{\sum\_{k=1}^{N} (T\_{in,k} - T\_{out,k})} \tag{3}$$

### *2.2. Machine Learning Models*

The different mathematical models (or black-box models) built to estimate the HLC of the analysed building are presented in this section. In this case, the three machine learning algorithms studied are XGBoost, SVR and MLP neural network.

#### 2.2.1. Extreme Gradient Boosting

Gradient boosting is a meta-algorithm built on the basis of weak learners, like decision trees [19], with the aim of obtaining a strong ensemble learner [46,47]. Specifically, XGBoost is a scalable and efficient way to carry out the implementation of a gradient boosting algorithm. It is remarkable for its ease of implementation and its high accuracy in predictions [20]. This algorithm focuses on the idea of combining, in a final step, all the predictions made by a set of weak learners (*additive training strategy* [48]). Moreover, XGBoost is capable of simplifying the objective function with different combinations of predictive and regularization terms without loss in the optimal computational speed [19,20]. The specific learning process is summarized as follows [18,20]:

$$\mathbf{G}\_{l}^{(t)} = \sum\_{j=1}^{l} \mathbf{G}\_{\dot{l}}(\mathbf{x}\_{l}) = \mathbf{G}\_{l}^{(t-1)} + \mathbf{G}\_{l}(\mathbf{x}\_{l}) \tag{4}$$

where *Gt*(**x***i*) is the learner at step *t*, *G*(*t*−1) *i* the prediction at step *t* − 1 and **x***i* the input variable.

On the other hand, to prevent the problem of overfitting while maintaining an optimal computational speed, XGBoost takes into account Equation (5) to evaluate the suitability of the model [17,20]:

$$\mathbf{O}^{(t)} = \sum\_{j=1}^{n} L(\mathbf{\hat{y}}\_{i}, \mathbf{y}\_{i}) + \sum\_{j=1}^{t} \mathbf{\hat{y}}(G\_{i}) \tag{5}$$

*n* being the number of observations, *L* a differentiable loss function and Ψ a regularization term that penalizes the complexity of the model [18].

XGBoost expands the loss function into the second order to be able to optimize the problem quickly. Furthermore, several specific techniques can reduce the possible overfitting of the algorithm [18]. One of them consists of, after each step of boosting, the weights recently added being scaled by a factor *η* (*shrinkage*). This process reduces the influence of the individual trees built, and therefore, the training will be slower and more efficient.

Lastly, XGBoost optimization requires the control of multiple parameters [19], and finding the best combination of them is important. In this analysis, the optimal combinations of parameters (*max\_depth,*

*min\_child\_weight, subsample, colsample\_bytree* and *learning\_rate*) are found through a *k*-fold cross-validation method (*k* = 10) [49] carried out on the different validation samples studied.

#### 2.2.2. Support Vector Regression

SVR is a nonlinear kernel based regression model that focuses on finding the best regression hyperplane with the least structural risk in a high-dimensional feature space [25–27]. The SVR function is represented in Equation (6):

$$\mathbf{g}(\mathbf{x}) = \mathbf{w}^T \boldsymbol{\delta}(\mathbf{x}) + \mathbf{b} \tag{6}$$

*δ*(*x*) being a nonlinear mapping that connects the input space to the feature space, *b* a bias term and *w* the weight coefficient vector. In this case, both *w* and *b* are estimated by resolving the following optimization problem [26,28]:

$$\begin{aligned} \underset{\mathbf{w}, \mathbf{b}}{\text{minimize}} & \quad \frac{1}{2} \parallel \mathbf{w} \parallel^2 + \mathbb{C} \sum\_{j=1}^{l} (\tau\_j + \tau\_j^\*) \\ \text{subject to} & \quad \mathbf{y}\_j - g(\mathbf{x}\_j) \leq \epsilon + \tau\_j \\ & \quad g(\mathbf{x}\_j) - \mathbf{y}\_j \leq \epsilon + \tau\_j^\* \\ & \quad \tau\_j, \tau\_j^\* \geq 0 \end{aligned} \tag{7}$$

where the constant *C* > 0 represents the trade-off between training error and model complexity,  corresponds to a threshold value and *l* is the number of training patterns. As Huang et al. [28] and Vrablecová et al. [29] showed, once having resolved the optimization problem (Equation (7)) and taking the Lagrangian, the model solution can be reached with its dual representation (see Equation (8)):

$$\log(\mathbf{x}) = \sum\_{j=1}^{l} (\mathbf{a}\_j^\* - \mathbf{a}\_j) K(\mathbf{x}, \mathbf{x}\_j) + \mathbf{b} \tag{8}$$

where *αj*, *α*∗ *j* are the Lagrangian multipliers (=0) and the solution for the dual problem, *b* the bias term and *K*(*xj*, *x*) the kernel function based on the inner product *δ*(*xj*), *δ*(*x*) .

Specifically, the model optimization for each of the validation samples studied is carried out through the tuning of certain parameters. In this case, the selected parameters were penalty *C* and term  [24,29]. With a cross-validation process (*k*-fold = 10) [49], different values of these parameters were evaluated, and the best adjustment was obtained depending on the validation sample considered.

2.2.3. Multi-Layer Perceptron Neural Network

MLP is an ANN characterized by having several layers [32,40,50]:

The specific neural network built in this study is composed of five layers (one input layer, one output layer and three hidden layers). The internal structure of neurons in each layer of the network is as follows: 100-100-100-50-1. This architecture is selected after a k-fold cross-validation (k = 10) study [51] in which the grid of different hidden layer options was between zero and four [52,53]. On the other hand, in relation to the number of neurons in the network layers, the different options, as recommended by Vujicic et al. [52] and Doukim et al. [54], take into account the size of the sample and the number of inputs and outputs. Moreover, due to the complexity of the problem, this grid of values was completed with more complex options such as: 50, 100, 200 and 400.

The neural network training was developed via backward propagation; errors were propagated and corrected through the network [55]. In this way, the real outputs must be known. In addition, the training process was based on updating the weights with an average update (batch-learning), which was obtained by introducing all the patterns in the input file (an epoch) and accumulating the weight updates [56,57]. The stop criterion used to avoid overfitting problems was the cross-validation due to its effectiveness in stopping with the best model generalization [58,59]. The training will stop when the performance of the neural network, measured by the Mean Squared Error (MSE) and evaluated with a small part of the whole sample (test sample), stagnates or starts to decrease. Each MLP trained in this study used the Rectified Linear Unit (ReLU) as the activation function [60], a normal kernel initializer and the Adaptive Moment Estimation (Adam) optimizer [61]. Lastly, the batch size with which the neural networks were trained was equal to 64, and the limit to the epochs in which the model fit did not improve was 100. Further information about the options and variations of the MLP training process can be found in [62,63].

### *2.3. Case Study Data Acquisition*

### 2.3.1. Description of the Building

The methodology proposed was tested on a public library of the Faculty of Marine Sciences at the University of Vigo (see Figure 1) that is located in the northwest of Spain. The building has three floors that are interconnected, and it is completely monitored and has been used for research purposes many times. Therefore, the building, its HVAC system and its data acquisition system were further and more deeply described in other articles such as Cacabelos et al. [64,65], Fernández et al. [66] and Martínez et al. [40]. The building envelope has a large, south-facing window (see Figure 1a), and its enclosures are made of different concrete and insulation material. Table 1 shows the properties and composition of each material layer of the walls, floor and roof.

(**a**) (**b**) **Figure 1.** Appearance of the analysed building: (**a**) exterior and (**b**) interior.

### 2.3.2. Pre-Processing Data

This analysis is focused on the estimation of the HLC of a building. The data available in this study to train the machine learning models were hourly observations between March 2016 and December 2017 of four variables. The variables considered were two describing the thermal behaviour of the building (thermal demand and indoor temperatures) and two describing the climate conditions (outdoor temperature and solar radiation). In addition, to find significant hours when the heating worked normally, only the hours with a thermal demand of more than 5 kW were selected (n = 5727). In this case, it was not necessary to have a continuous sample because the training did not take into account time lags in

the explanatory variables. Therefore, three time variables (hour of the year, day of the week, and hour of the day) were introduced into the models to provide more information and thus to improve the training. The variables that were used as model inputs are presented in Figure 2.

**Figure 2.** Summary of the training and prediction process. On the left, the input variables of the models. On the right, the objective estimations of the Heat Loss Coefficient (HLC).

The model fit was carried out through a cross-validation process in which the whole sample was divided consecutively into two individual samples: training and testing. In the search for the optimal model and the selection of hyperparameters, the models were tested with different partitions (or subsamples) of the whole training sample (*k*-fold = 10) [49]. In addition, the models, after being trained and tested with 10 different test samples, were validated with three independent samples to prove their efficiency in estimating the HLC of the building:

These three samples had to comply with certain restrictions due to the specification of the HLC formula (see Section 2.1) and because of the singularities of the monitoring of the building studied (see Section 2.3.1). The restrictions for considering a period suitable for the calculation of the HLC are [8,44]:

### 2.3.3. Validation and Error Assessment

The Coefficient of Variation of the Root Mean Squared Error (CV(RMSE)) and the Normalized Mean Biased Error (NMBE) were the error measures calculated to evaluate the accuracy of the models presented:

$$\text{CV(RMSE)} = 100 \times \frac{\sqrt{\sum\_{i=1}^{N} (y\_i - \hat{y}\_i)^2 / N}}{\bar{y}} \tag{9}$$

$$\text{NMBE} = 100 \times \frac{\sum\_{i=1}^{N} (y\_i - \hat{y}\_i)}{\sum\_{i=1}^{N} (y\_i)} \tag{10}$$

Both measurements were used to compare the performance of the different models in the three validation samples. In addition, the CV(RMSE) was also taken into account to select the best trials of each model (with the lowest error) and, then, represented graphically in the Results Section. These measures have been used and their efficiency has been proven in similar studies [42,67,68].

### **3. Results and Discussion**

A methodology, based on thermal demands' predictions made with black-box models, to estimate the HLC of a building is developed in this paper. In particular, the building under study was the public library of the Faculty of Marine Sciences at the University of Vigo. The available data were hourly observations of the variables presented in Figure 2 from March 2016 to December 2017; only taking into account the hours with a significative thermal demand (≥ 5 kW). In this analysis, three specific weekends were considered to study the performance of the models estimating the HLC through heating demand predictions. A comparison between the accuracy of each ML model analysed is presented in the following sections.

Section 3.1 presents the models performance in heating demand predictions and Section 3.2 a similar analysis for the HLC estimations. While in the thermal demand analysis, CV(RMSE) (Equation (9)) and NMBE (Equation (10)) were calculated for each of the models, in the HLC section, the errors were calculated based on an absolute variation rate. Furthermore, to represent the average performance of the different models, the predictions of the three validation samples were repeated 10 times (varying the subsamples where the model was tested). Thus, the numerical results shown are average errors obtained through the 10 trials, as well as their standard deviations. Finally, all figures shown in the following sections were created with the Python programming language [69].

### *3.1. Thermal Demand Analysis*

The results of the heating demand predictions for the analysed building in the three validation samples are presented in this section in Figure 3 and Table 2. These predictions are important as they will be used to obtain the subsequent HLC estimations (see Section 3.2). Furthermore, in Figure 3, each of the algorithms is represented by the prediction that obtains the lowest CV(RMSE) among the 10 trials (best scenario).

**Figure 3.** Results of thermal demand predictions in all samples studied: (**a**) Sample 1; (**b**) Sample 2; (**c**) Sample 3. Each machine learning model is represented by the prediction curve with the lowest CV(RMSE) among the 10 experiment trials.

In the case of Sample 1 (Figure 3a) and taking into account the CV(RMSE) results, the XGBoost model presents the lowest average error (18.57%). The MLP neural network and SVR model show a worse performance with an average error of 19.84% and 21.31%, respectively (see Table 2). Moreover, Table 2 shows that in all models, the variability of errors in the 10 trials was relatively low (below ± 4), but MLP is the one that shows the lowest dispersion (±2.41). On the other hand, regarding the NMBE results, the lowest average error was obtained by the MLP model (−5.49%), while the other algorithms only managed to obtain negative average errors higher than 10% (−12.52% for XGBoost and −11.63% for SVR). In addition, in this case, while the variability of the NMBE results of XGBoost and MLP are similar and around ± 5 (see Table 2), the SVR results show a greater dispersion (above ± 8). In general, all the models, in the best scenario, efficiently reproduced the real behaviour of the thermal demand of the building analysed (see Figure 3a). Although the XGBoost model obtained better results in relation to the CV(RMSE), the MLP predictions showed a better overall adjustment to reality by obtaining a similar average CV(RMSE) and an average NMBE much lower than XGBoost.

In Sample 2 (Figure 3b), considering the CV(RMSE) results, the XGBoost algorithm and MLP neural network show a similar average performance. MLP obtained an average error of 17.38% and XGBoost 17.43%. In this situation, the SVR model shows better results than the other models with an average error of 16.61% (see Table 2). Additionally, while SVR and MLP present a similar error dispersion around ± 2, the XGBoost model is the one that presents the highest variability (above ± 4). With respect to the NMBE results, the SVR model again obtained the lowest average error (2.99%). The MLP neural network and XGBoost algorithm, on the other hand, present higher and negative average error (−4.01% and −8.40%, respectively). However, as shown in Table 2, the SVR model is the one that showed the greatest dispersion in the NMBE results (almost ± 7), even though the other models also showed an important variability. In Figure 3b, it is demonstrated that each of the algorithms, taking into account the best scenario, is capable of replicating the reality except for certain peaks. However, more specifically, the model that yielded the best average performance was SVR.

Regarding Sample 3 (see Figure 3c), the three models studied obtained average CV(RMSE) values above 20%. While the SVR algorithm and MLP neural network showed an average CV(RMSE) of 21.60% and 21.54%, respectively, XGBoost model obtained an average error around 29% (see Table 2). Furthermore, taking into account that the variabilities presented were not high, the MLP neural network obtained the least variable results (±2.33) and XGBoost the most variable (±3.81). In the case of the NMBE results, the MLP model stands out from the rest (see Table 2). It shows the lowest average of NMBE by far (0.25%). The average results of the other algorithms were higher than 10% and negative (−13.60% for SVR and −23.05% for XGBoost). On the other hand, the dispersion of the NMBE values among the 10 trials of this sample is high (above 4.5) in all the models, and the MLP neural network was the model with the greatest

error dispersion (±6.19). In addition, as in the other validation samples, Figure 3c shows that the built models, in the best case scenario, are very close to the real values.

**Table 2.** Numerical results of the thermal demand predictions for all validation samples. The mean of the CV(RMSE) and the Normalized Mean Biased Error (NMBE) obtained through the 10 trials, besides the standard deviation (SD) of each of them, are presented.

Lastly, it is demonstrated that all the models presented, in general, are capable of recreating reality (see Figure 3), obtaining optimal errors in all the samples studied (see Table 2). Concretely, they normally predict above the real values of the heating demand of the building analysed (see the negative NMBE results in Table 2). Furthermore, encompassing all the results, the model that performed best was the MLP neural network. In terms of CV(RMSE), it always yielded one of the best results, and in terms of NBME, it obtained much better results than the other algorithms in two of the three validation samples.

### *3.2. HLC Estimation Analysis*

The specific characteristics of each of the validation samples studied in this work are analysed in Table 3 and presented in Figure 4. Moreover, the results of the HLC estimations, based on the previous heating demand predictions, for each of the three validation samples are shown in this section in Table 4. Although estimated from the heating demand predictions during a weekend, the HLC is represented as a single number (see Equation (3)). As in the preceding section, the average performance of the different algorithms among the 10 trials is summarized in Table 4.

**Table 3.** Thermal conditions of the three validation samples. These represent the mean difference between indoor and outdoor temperatures, the average heating demand (*Q*), the solar radiation gain (Rad) over the thermal demand and the average temperatures (indoor and outdoor) at the start (*T*◦ *initial*) and at the end of the period (*T*◦ *final*). HLC, Heat Loss Coefficient.

First, Table 3 shows that the three validation samples fulfilled the conditions necessary to efficiently measure the HLC of the analysed building (see Section 2.3.2). Moreover, each of the samples has different HLC values. The main different thermal conditions of the samples, which caused slightly different HLC values, are also summarized in Table 3 and presented in Figure 4. The highest HLC value was obtained in Sample 1 (2.75 kW/K) due to the fact that the weight of the radiation gains, in relation to the heating demands, was the lowest among the samples (0.09%). This is related to the fact that the average thermal demand throughout this period was the highest and is a numerator value in the HLC formula (see Equation (3)). In addition, the average difference between indoor and outdoor temperatures, as shown in Figure 4a, was one of the smallest (12.48 K). These were the main reasons for obtaining a greater HLC

(see Table 3). On the other hand, Sample 2 is where the HLC value is the lowest (2.15 kW/K) because, as shown in Figure 4b, the average difference between indoor and outdoor temperatures was the largest (15.25 K). This value, which is in the denominator of the HLC formula, together with an average heating demand significantly lower than in Sample 1 (32.72 kW) reduces the calculated HLC value (see Table 3). Lastly, in Sample 3 (see Figure 4c), an intermediate value of the HLC (2.46 kW/K) that came from the lowest average difference between indoor and outdoor temperatures (12.05 K) and the lowest average thermal demand (29.70 kW) was obtained.

**Figure 4.** Summary of the thermal conditions of each of the validation samples: (**a**) Sample 1; (**b**) Sample 2; (**c**) Sample 3. Solar radiation and indoor and outdoor temperatures affecting the analysed building during the specific time periods are represented.

On the other hand, in the case of Sample 1, where the calculated HLC was 2.75 kW/K, the MLP neural network was the most accurate model (see Table 4). While this algorithm obtained an average absolute variation rate of 6.50%, the XGBoost and SVR models were only able to obtain an average absolute variation greater than 10% (12.52.% and 12.16%, respectively). Thus, the average HLC value estimated by MLP (2.90 kW/K) is the closest to the measured HLC value. Table 4 shows that all models obtained a higher average estimation than the calculated one (the same situation as in Section 3.1). Additionally, in relation to the variation rate dispersion, the MLP neural network was the one with the lowest variability among the errors (±4.5). However, in general, all models showed a significant high standard deviation.

**Table 4.** Numerical results of the HLC estimations for all validation samples. The values shown summarize

In Sample 2 the calculated HLC was 2.15 kW/K, and in relation to the average absolute variation values, the most accurate model was the MLP neural network (4.08%), but close to the SVR model (5.68%). As in the thermal demand section, the XGBoost model performed worse with an average absolute variation rate of 8.40% (see Table 4). Therefore, the HLC values estimated by the SVR and MLP models are very close to the calculated HLC value; while the average value obtained by SVR is 2.08 kW/K, the average value estimated for MLP is 2.23 kW/K. With respect to the dispersion of variation rate data and, as in Sample 1,

Table 4 shows that the MLP neural network was the most stable model with a standard deviation below ±3. The other models presented values above ±4.

Regarding Sample 3, in which the measured HLC was 2.46 kW/K, the MLP neural network was again the model with the best average performance. While the XGBoost and SVR algorithms presented average absolute variation of 23.05% and 13.60%, respectively, the MLP model showed an average absolute variation of 4.97% (see Table 4). Therefore, the average HLC value from the MLP estimations (2.46 kW/K) was much closer to the calculated HLC value than those obtained by the other models. In addition, in this sample, all models obtained a high variability in their results: all standard deviations were higher than ±3.5. In this situation, the MLP neural network was the model with the lowest dispersion among its errors (±3.70).

Definitely, the HLC value that characterizes the studied building, calculated as the average of the presented results, was 2.45 ± 0.30 kW/K. Observing the results presented in Table 4 and taking into account the whole analysis (the results of the three validation samples), the model that presents the best average performance was the MLP neural network. It was the model that obtained the most stable results and had the lowest average absolute variation rate for all the samples analysed. In this way, it is demonstrated that it is possible to estimate the HLC of the analysed building with an absolute error around 4 or 6% if an MLP model is used to make the necessary thermal demand predictions. On the other hand, if an SVR algorithm is used, the error increases to 5-13%, and if the model is XGBoost, the errors vary between 8 and 23%.

### **4. Conclusions**

A new methodology for estimating the HLC of a building is presented in this paper. It is based on the introduction of thermal demand predictions obtained with machine learning models in the HLC formula. This study focuses on the analysis, on the one hand, of monitored data on heating demands and indoor temperatures belonging to the Science Library of the University of Vigo. On the other hand, two meteorological variables (outdoor temperature and solar radiation) and three temporal variables (hour of the year, day of the week and hour of the day) are also taken into account. The aim of this paper is to show a methodology that allows the efficient estimation of the HLC value of a building without the need to control its heating demands (nor indoor temperatures if they are assumed to be constant). The search for the optimal methodology considers and compares three different machine learning models (XGBoost, SVR and MLP neural network). In addition, the performance of each one is evaluated and analysed both through its average accuracy in thermal demand predictions and its average accuracy in HLC estimations.

The research contribution of this work is the application of mathematical models to estimate the HLC of a building with low errors. In addition to reducing the necessity for monitoring, these models can be useful for detecting errors in measurements from sensors installed in the building. Moreover, the black-box models presented contributes with advantages compared to traditional research in building simulation. The use and application of the traditional building thermal simulation models are conditioned by the need for significant knowledge on a subject. In addition, these models need to control many different parameters related to the energy performance of a building. Nevertheless, machine learning models, which need less time for development than dynamic simulation methods, do not require specific prior knowledge and can be applied in numerous fields. The only important need for these models is the availability of a significant amount of data from which a behaviour pattern is extracted. Furthermore, the inputs introduced in the built models are variables typically monitored in buildings. Therefore, the methodology presented here is extractable to other studies and buildings.

The results obtained show that it is possible to efficiently estimate the HLC value of a building over specific time periods with black-box models. To this end, it is important to obtain thermal demand

predictions, which are used as inputs in HLC estimations, with average errors lower than the values proposed for the calibrated models. The results also demonstrate that the MLP neural network is the algorithm with the best average performance in HLC estimation (more stability and higher average accuracy in all samples studied). The SVR model shows a close average behaviour, but XGBoost, except in Sample 1, presents a much worse performance than the other two. In the first validation sample, in which the measured HLC is 2.750 kW/K, only the MPL neural network is capable of obtaining an average absolute variation rate below 10% (6.50%). The SVR and XGBoost models obtain an average variation of 12.16% and 12.52%, respectively. On the other hand, in the second validation sample, where the calculated HLC was 2.146 kW/K, both the MLP and SVR models performed better than XGBoost. While the first two obtain an average absolute variation rate of 4.08% and 5.68%, respectively, XGBoost shows an average variation of 8.40%. Lastly, in the case of the third validation sample, with a measured HLC value of 2.464 kW/K, the MLP neural network yields an average absolute variation rate far from the other algorithms (4.97%). The SVR and XGBoost models are only able to yield an average absolute variation of 13.60% and 23.05%, respectively. Furthermore, regarding the dispersion in error data, MLP model is the one that shows the most stable results in all validation samples. The other algorithms present similar variability between them, but higher than the one obtained by the MLP model. Definitely, taking into account all the results presented, the HLC value of the analysed building is 2.45 ± 0.30 kW/K.

From an energy point of view, the conclusion is that efficient predictions related to the thermal conditions of a building and, in addition, made by machine learning models can be used to efficiently estimate its HLC. This is demonstrated in this paper with three different validation samples that fulfil the necessary conditions to calculate the HLC. The most accurate model, which in this case is the MLP neural network, is able to estimate the HLC of the analysed building with an average absolute variation rate of around 5% and with a standard deviation of around ±3. On the other hand, the main limitation of this research is the many restrictions on finding suitable time periods for calculating the HLC. In this case, the unavailability of occupation data means that only weekends are studied. For this reason, some possible future lines of research are similar analyses considering more data such as occupation or extending the study to more general situations.

**Author Contributions:** Conceptualization, L.F.-G. and M.M.-C.; methodology, M.M.-C.; software, M.M.-C.; validation, L.F.-G., E.G.-Á. and S.M.-M.; formal analysis, L.F.-G., J.M.-T. and S.M.-M.; investigation, M.M.-C., L.F.-G. and S.M.-M.; resources, E.G.-Á.; data curation, M.M.-C. and J.M.-T.; writing, original draft preparation, M.M.-C. and L.F.-G.; writing, review and editing, L.F.-G., J.M.-T., E.G.-Á. and S.M.-M.; visualization, M.M.-C.; supervision, E.G.-Á., L.F.-G. and J.M.-T.; project administration, E.G.-Á.; funding acquisition, E.G.-Á. All authors read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Spanish Government (Science, Innovation and Universities Ministry) under the project RTI2018-096296-B-C21.

**Acknowledgments:** This research was supported by the Spanish Government (Science, Innovation and Universities Ministry) under the project RTI2018-096296-B-C21.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

### **Artificial Intelligence, Accelerated in Parallel Computing and Applied to Nonintrusive Appliance Load Monitoring for Residential Demand-Side Management in a Smart Grid: A Comparative Study**

### **Yu-Chen Hu 1, Yu-Hsiu Lin 2,\* and Chi-Hung Lin 2**

Received: 23 October 2020; Accepted: 13 November 2020; Published: 16 November 2020

**Abstract:** A smart grid is a promising use-case of AIoT (AI (artificial intelligence) across IoT (internet of things)) that enables bidirectional communication among utilities that arises with demand response (DR) schemes for demand-side management (DSM) and consumers that manage their power demands according to received DR signals. Disaggregating composite electric energy consumption data from a single minimal set of plug-panel current and voltage sensors installed at the electric panel in a practical field of interest, nonintrusive appliance load monitoring (NIALM), a cost-effective load disaggregation approach for (residential) DSM, is able to discern individual electrical appliances concerned without accessing each of them by individual plug-load power meters (smart plugs) deployed intrusively. The most common load disaggregation approaches are based on machine learning algorithms such as artificial neural networks, while approaches based on evolutionary computing, metaheuristic algorithms considered as global optimization and search techniques, have recently caught the attention of researchers. This paper presents a genetic algorithm, developed in consideration of parallel evolutionary computing, and aims to address NIALM, whereby load disaggregation from composite electric energy consumption data is declared as a combinatorial optimization problem and is solved by the algorithm. The algorithm is accelerated in parallel, as it would involve large amounts of NIALM data disaggregated through evolutionary computing, chromosomes, and/or evolutionary cycles to dominate its performance in load disaggregation and excessively cost its execution time. Moreover, the evolutionary computing implementation based on parallel computing, a feed-forward, multilayer artificial neural network that can learn from training data across all available workers of a parallel pool on a machine (in parallel computing) addresses the same NIALM/load disaggregation. Where, a comparative study is made in this paper. The presented methodology is experimentally validated by and applied on a publicly available reference dataset.

**Keywords:** artificial intelligence; artificial neural network; demand-side management; evolutionary computing; non-intrusive appliance load monitoring; parallel computing; smart grid; smart house

### **1. Introduction**

Nonintrusive appliance load monitoring (NIALM), also called nonintrusive load monitoring, was first investigated by George W. Hart et al. [1] at the Massachusetts Institute of Technology with funding from the Electric Power Research Institute in the early 1980s. It has been considered as a cost-effective alternative against intrusive load monitoring approaches that involve the deployment

of plug-load power meters (smart plugs) for individual concerned electrical appliances in a practical field of interest (it cannot be costed down for the realization of load management) [2], and has been developed for (residential) demand-side management (DSM) in a smart grid [3,4]. With DSM implemented in a smart grid, consumers have the opportunity to improve their awareness of what and when their monitored electrical appliances should operate for use in response to demand response (DR) schemes [5]. Thus, power utilities coming up with DR schemes plan to address ever-increasing electric demand in an optimal way [6].

NIALM, load disaggregation by estimating appliance-by-appliance energy consumption from composite electric energy consumption data, is becoming mature with the improvement of advanced metering infrastructure for DSM in a smart grid [7,8], where load disaggregation is performed with/applied on smart meter data.

NIALM can be built upon signal processing [9], machine learning [4,9–20], and deep learning [20–24]. Not much attention has been paid to evolutionary computing, wherein load disaggregation is considered as a combinatorial optimization problem [1]. As a result, in this paper, load disaggregation is declared as a combinatorial optimization problem, and is solved by a genetic algorithm (GA) accelerated in parallel computing. GA used to solve NIALM that is declared as a combinatorial optimization problem would involve large amounts of electric energy consumption data gathered from smart meters deployed in practical fields of interest in a smart grid, population-based candidate solutions evaluated, and/or evolutionary cycles executed from iteration to iteration, which is a highly demanding task in terms of performance in load disaggregation and computational time in evolutionary computation. Hence, GA conducted in this paper for NIALM is parallelized during evolution. The GA-based NIALM in this paper is processed and accelerated in consideration of parallel computing. It has been experimentally validated by and applied on a publicly-available reference dataset considered in this paper. Besides the GA-based NIALM, a feed-forward, multilayer artificial neural network (ANN)-based NIALM against other different types of ANNs (deep NNs)-based NIALM addresses the same NIALM data from the considered publicly available reference dataset. A comparative study is made in this paper. ANNs used in this paper to address NIALM can learn from NIALM data across parallel workers on a machine supporting parallel computing.

The remainder of this paper is structured below. The investigated methodology is presented in Section 2. Section 3 gives the experimentation and results. Conclusions are drawn in Section 4, where future work is also anticipated.

### **2. Methodology**

Figure 1 describes a basic (eventless) NIALM process consisting of (1) data acquisition: composite electric energy consumption data are measured by a single minimal set of plug-panel current and voltage sensors and digitized for further analysis of load disaggregation; (2) feature extraction: electrical features are extracted as feature data from digitized composite electric energy consumption data for concerned electrical appliances; and (3) load recognition: artificial intelligence (AI) is utilized to recognize extracted electrical features for concerned electrical appliances—their electric demand based on past trends can be predicted. The NIALM investigated in this paper is an eventless NIALM approach. AI conducted in the NIALM is parallel computing accelerated evolutionary computing. The basic NIALM problem can be formulated as follows:

**Figure 1.** A basic NIALM process (profiling concerned electrical appliances and statistically computing base load should be done beforehand).

If (1) *P*(*t*) stands for used composite power consumption acquired at time *t* and (2) *Pi*(*t*) accounts for real power absorbed by concerned electrical appliance *i* (profiling concerned electrical appliances and base load is done in advance), the total power consumption, *P*(*t*) (=*y*(*t*)), in an electric power distribution system can be expressed as Equation (1). In Equation (1), the summation of superimposed absorptions Σ*xi*(*t*)*Pi*(*t*) with an unmonitored base load *Pbase*(*t*) is made. In Equation (1), *N* is the total number of concerned electrical appliances and *xi*(*t*), a unit value from {0, 1}, is the on/off operational status of electrical appliance *i* at time *t*. NIALM is used to recognize the unknown (possible) operational status of concerned electrical appliances *xi*(*t*)|*i* = 1, 2, ..., *N* at time *t* when *P*(*t*) (=*y*(*t*)) acquired and given with *Pbase*(*t*) is known: [*x*1(*t*), *x*2(*t*), ..., *xi*(*t*), ..., *xN*(*t*)] = *F*(*y*(*t*)) where *F* is a function that returns the best *N* estimates of *xi*(*t*) at time *t* for concerned individual *N* electrical appliances. *F* above can be any AI algorithm. In this paper, we conduct evolutionary computing, GA, to metaheuristically search for the optimal load combinations *xi*(*t*) of concerned electrical appliances when acquired total power consumption *P*(*t*) with base load *Pbase*(*t*) is given at time *t* (in Equation (1), profiling concerned electrical appliances to get *Pi*(*t*) and statistically computing base load for *Pbase*(*t*) from unmonitored electrical appliances is done prior.

$$P(t) = \sum\_{i=1}^{N} \chi\_i(t) P\_i(t) + P\_{\text{base}}(t) \tag{1}$$

In this paper, NIALM is declared as a combinatorial optimization problem, which is formulated as the objective metric in Equation (2) and solved by a parallel computing accelerated GA.

$$err(t) = \operatorname\*{argmin}\_{\left| \left[ P(t) - P\_{\text{base}}(t) \right] - \sum\_{i=1}^{N} \chi\_i(t) \left[ P\_i(t) + \tau\_i \right] \right| \tag{2}$$

where τ*i* = *ci* ·*std*(*Pi*(*t*)).

In Equation (2), *Pi* (=mean(*Pi*(*t*))) obtained in advance involves profiling each of concerned electrical appliances based on historical data. That is, *Pi*(*t*), in Equation (1), or *Pi*, in Equation (2), accounts for the real power that was absorbed by electrical appliance *i* and statistically computed for load disaggregation. *Pbase*, in Equation (2), from *Pbase*(*t*), in Equation (1), is defined in a similar sense, and its standard deviation can be considered. In Equation (2), τ*i*, a tolerance term for *Pi*, is considered, where std(*Pi*(*t*)) is a function that returns the standard deviation of its input elements *Pi*(*t*) from historical data and *ci* is a constant that can be designed as a time-dependent parameter.

Figure 2 illustrates the principle of combinatorial search for load combinations with *xi*(*t*) in Equation (2) for load disaggregation, which aims to optimize load combinations by *xi*(*t*). In Figure 2, *Pi*(*t*) represents the assumed power demand of the *i*-th electrical appliance. *xi*(*t*) associated with *Pi*(*t*) is represented as a binary vector, evolved through metaheuristics, and used to minimize Equation (2) (the minimal error, *err*(*t*), to be obtained for load disaggregation) between the summation of superimposed absorptions Σ*xi*(*t*)*Pi*(*t*) with an unmonitored base load *Pbase*(*t*) and the total load *P*(*t*) (*P*(*t*) to be approximated from a pack of *Pi*(*t*) with *Pbase*(*t*)). A metaheuristic algorithm, GA, is suitable for load disaggregation formulated in Equation (2) (the objective metric is the metric we are trying to optimize and the fitness metric is the algorithm's guide to doing so), where concerned electrical appliances are recognized by parallel computing accelerated GA for the declared objective metric. Concerned electrical appliances are constant or time-varying resistive, inductive, or capacitive loads.

**Figure 2.** An illustration of the principle of combinatorial search for load combinations to *xi*(*t*) for load disaggregation in this paper, whereby load combinations by *xi*(*t*) are optimized. *xi*(*t*) is represented as a binary vector, evolved through metaheuristics, and used to minimize Equation (2) between the summation of superimposed absorptions Σ*xi*(*t*)*Pi*(*t*) with an unmonitored base load *Pbase*(*t*) and the total load *P*(*t*).

The NIALM investigated in this paper has the main objective of recognizing concerned electrical appliances for (residential) DSM according to a composition of appliance-level real power consumption, which is disaggregated from total real power consumption acquired apart for load disaggregation. That is, with the minimum in Equation (2), acquired total real power consumption is the total of real power consumption by all individual electrical appliances concerned and operated where a base load

should be considered. To Equation (2), a parallel computing accelerated GA as a load recognizer for load recognition of the NIALM in this paper is used to recognize the correct *xi*(*t*) to gain the minimum between the acquired total real power consumption and the sum of superimposed real power consumption by concerned electrical appliances. In this paper, real power consumption, *P*, is extracted, as the electrical feature for load recognition in Figure 1, from acquired total real power consumption, which is disaggregated into appliance-level real power consumption through the parallel computing accelerated GA-based NIALM of Equation (2). Concerned electrical appliances to be recognized can be constant or time-varying resistive, inductive, and capacitive loads [4,25].

### *2.1. GA-Based NIALM*

GAs are a stochastic, population-based metaheuristic optimization (search) algorithm that can search for the global (quasi*-*)optimal solution metaheuristically to both constrained and unconstrained, NP-hard/NP-complete discrete or continuous optimization problems addressed, through (natural) selection operations that select current chromosomes for further gene recombination and through genetic operations involving (1) crossover operations that combine genetic information from selected parents' chromosomes to produce new offspring and (2) mutation operations that provide genetic diversity among population members. The evolutionary cycle is based on the principle of biological evolution. GAs that operate based on the principle of biological evolution have been proven to be an effective metaheuristics technique in many optimization problems.

In GAs, genetic operations, comprising crossover and mutation operations, and evolutionary operations, implementing the driving force of evolution (natural selection/Darwinism), are involved. See [26] for more details about a standard GA. In a GA addressing an optimization problem, a population of chromosomes as candidate solutions in a search space is evolved, from generation to generation, through genetic operations and evolutionary operations towards a global (quasi-)optimal solution. In an evolutionary cycle in a GA, a proportion of an existing population in its current generation is reproduced through selection operations and is used through crossover and mutation operations to breed a new population for the next generation. Selection operations select chromosomes based on a fitness-based selection procedure where chromosomes as candidate solutions are evaluated by a fitness function, and chromosomes with relatively high fitness are typically more likely to be selected. Also, crossover operations exchange subparts of two chromosomes, roughly mimicking the biological recombination between two haploid organisms. Finally, mutation operations randomly alter genes in some chromosomes, where an arbitrary bit of a chromosome is changed from its original state. In a GA, the population size depends on the nature of the problem addressed, but, typically, the population contains several hundreds or thousands of chromosomes/candidate solutions. Chromosomes are generated at random in the initial population. The evolutionary cycle is repeated, from generation to generation, until a termination condition such as a maximum of generations prespecified and reached has been met.

A GA used to address optimization problems is innately a parallel algorithm that can be run on a multicore processor. The workflow depicting the parallel computing accelerated GA used to solve Equation (2) for load recognition in the NIALM in this paper is shown in Figure 3. In the GA, function evaluations are farmed out to different processors on a multicore processor; they are executed in parallel. Figure 4 shows that *xi*(*t*) (at time *t*) in Equation (2) are encoded as a chromosome for an evolutionary cycle of the GA. To evaluate a chromosome with Equation (2), the GA decodes it as a computed summation of superimposed absorptions Σ*xi*(*t*)*Pi*(*t*). With an initial population of randomly generated chromosomes that are started in the parallel computing accelerated GA in Figure 3, successive generations are constructed through evolution. Fitter chromosomes are more likely to survive, based on selection, and to participate in genetic operations [27]. Here, (1) a roulette selection strategy choosing parents by simulating a roulette wheel in which the area of the section of the wheel corresponding to an individual is proportional to the individual's expectation associated with its scaled fitness value is conducted (a ranking method that scales raw fitness values based on the rank of each

individual rather than its raw fitness value is used for fitness scaling), (2) a single-point crossover operator is used, (3) a bit-wise mutation operator is utilized, and (4) an elitist strategy guaranteeing that the solution quality does not decrease during evolution is also conducted.

**Figure 3.** A workflow of the parallel computing accelerated GA used to solve Equation (2) for load recognition in the NIALM in this paper, where the evolutionary cycle is parallelized.

**Figure 4.** A chromosomal design encoding *xi*(*t*) (at time *t*) in Equation (2).

### *2.2. Feed-Forward, Multilayer ANN-Based NIALM*