CEEMD, proposed by Torres et al. [28]. The final decomposition result is obtained by taking the average by adaptively adding white noise to the time series with the same magnitude and opposite direction. After engaging in EMD processing, the output of CEEMDAN obeys a Gaussian distribution. Compared with EEMD and CEEMD, this method effectively eliminates the problem of mode mixing, and no matter how many times the decomposition, the reconstruction error of the signal is almost zero, the completeness is better, and at the same time solves the problem of low decomposition efficiency, greatly reducing the calculation cost. CEEMDAN's spectra shows a more accurate decomposition of the frequency than the EEMD's. The same time series, the number of iterations of CEEMDAN shifting is usually equivalent to half of the EEMD's.

The implementation of the algorithm is summarized as follows:

Generate a white noise plus series (*xi*(*t*) = *x*(*t*) + *wi*(*t*)). Here, *wi*(*t*) denote the white noise of finite variance, while *x*(*t*) represent the original data. Then decompose *xi*(*t*) to obtain the IMFs:

$$\overline{\lim f\_i} = \frac{1}{n} \sum\_{j=1}^{n} \inf\_i^j(t) \tag{12}$$

where *wj i* (*t*) is the white noise added at the *j*-th time with the mean equal to zero and variance equal to one. The *im fi* is the *i*-th mode component obtained after the signal is decomposed by CEEMDAN.

Calculate the first residue:

$$r\_1(t) = x(t) - \overline{i \mathfrak{m} f\_1} \tag{13}$$

Decompose the residue to get the second IMF and calculate the second residue:

$$\overline{\lim\!f\_2} = \frac{1}{n} \sum\_{j=1}^{n} E\_1 \left\{ r\_1(t) + \varepsilon\_1 E\_1 \left[ w\_2^j(t) \right] \right\} \tag{14}$$

$$r\_2(t) = r\_1(t) - \overline{imgf\_2} \tag{15}$$

Decompose the (*k* − 1)th residue of*rk*−1(*t*) and extract the *im fk*(*t*). The process can be demonstrated in the following equation (*k* = 2, 3, ... *K*):

$$\overline{\lim f\_k(t)} = \frac{1}{n} \sum\_{j=1}^n E\_1(r\_{k-1}(t) + w\_{k-1}^j E\_{k-1}(\varepsilon^j(t))) \tag{16}$$

Calculate the *k*-th residue:

$$r\_k(t) = x(t) - \sum\_{i=1}^k \operatorname{im} f\_i \tag{17}$$

where *k* indicates the number of IMFs. The attributes of the original time series are denoted by all the IMFs exacted from different time-scales. The only residue demonstrates the trend, which is smoother than that of the original time series.

#### *3.2. GRU and Multi-LayerML-GRU*

#### 3.2.1. RNN and GRU

Among deep learning methods, RNN is a powerful method for processing time series data. It is widely used in many fields such as finance [31,32], industry and engineering [33], machine translation [34], speech recognition [35], economic prediction [36], and so on. As shown in Figure 2, RNN has certain information persistence capabilities, which enable information to be passed from one time-step to the next. However, the classic RNN does not have the ability to store and memorize data for a long time, so that it cannot capture long-term historical information. In addition, in the reverse process of model training, once the sequence is too long, the RNN will cause the problem of gradient explosion or gradient disappearance. To overcome these drawbacks, Hochreiter and Schmidhuber [37] proposed the LSTM neural network which is capable of learning long-term dependencies of time

series, as well as forgetting the worthless information based on the current input. LSTM has since replaced RNN as the most widely used recurrent neural network. LSTM has four gated unit that capable adaptively regulate the information flow inside the unit. An LSTM neural network usually requires many gated units, which need train a large number of parameters and occupy more computing resources. In order to reduce the training parameters and simply the neural network, Cho et al. [38] proposed a neural network named GRU which only have two gated units in the hidden unit.

**Figure 2.** Workflow graph of RNN. *X* denotes current input of neural network; *Y* denotes the output of neural network; *W*, *U* and *O* are the parameters of recurrent neural network; *h* is the neuron state of the hidden layer; The state at time *t* is related to the current input *X* and the hidden at time *t*−1.

As one of the variants of RNN, the input and output structure of GRU neural network is the same as that of RNN and LSTM, as shown in the Figure 3 The neurons of GRU receive the hidden state (*ht*−1) of the neuron of the previous neuron, and the current input *xt*. After passing through the gating unit, the neural network gets the output *yt* and passes the hidden state *ht* to the next neuron.

**Figure 3.** The internal structure of the GRU cell. *rt* —Reset gate. It allows the cell to forget certain parts of the state; *zt*—Update gate.

The GRU architecture, like LSTM, learns the long-term dependence of time series based on a gate mechanism that includes reset gate and update gate. The former is used to control how much information in the previous state will be ignored. The latter is used to control how much information from the historical state is brought into the current state. GRU neural networks, like other neural networks, consist of a large number of basic neurons. They are interconnected in a complex network. A single neuron cell is shown in Figure 2. The specific process is designed as follows:

$$z\_t = \sigma(w\_{\text{xx}}\mathbf{x}\_t + w\_{\text{hz}}\mathbf{h}\_{t-1})\tag{18}$$

$$r\_t = \sigma(w\_{\text{x}r}\mathbf{x}\_t + w\_{\text{hr}}\mathbf{l}\_{t-1})\tag{19}$$

The internal structure of the GRU cell is shown in Figure 3.

The update gate is used to control the extent to which the state information from the previous moment is retained to the current state. The more the value of the update gate approaches 0, the more the state information from the previous moment is brought into the current state. The update gate signal is the closer to 0, the more data it remembers. The closer to 1, the less it is forgotten.

Begin with getting the gating signal, GRU gets the reset data through the gate; then combine it with the input *xt*, and use a tanh activation function to shrink the data to the range from −1 to 1. The formula is shown in the following:

$$\log\_t = \tan h (w\_{\text{xy}} x\_t + w\_{\text{lg}}(r\_t \times h\_{t-1})) \tag{20}$$

The last step of GRU, we can call it "update memory" phase. Combined with the previous discussion, this step forgets some of the dimensional information passed in and add some new information inputted by the current neuron to the state variable *ht*:

$$h\_t = (1 - z\_t) \times h\_{t-1} + z\_t \times g\_t \tag{21}$$

#### 3.2.2. Multi-Layer GRU Architecture

If the problem is too complicated, a recurrent neural network with a single layer structure is not enough to abstract the problem, and a multi-layer neural network is a better alternative. The multi-layer neural network has more hidden layers and more powerful computing capabilities which is the key to solving complex problems. The proposed model is one kind of forward multi-layer neural networks, which demand the input layer are required to have the same number of input dimensions as the input vector. The architecture and workflow of the ML-GRU are shown in Figure 4. The output dimensions of the last layer demand only be equal to the number of labels for classification or a single value of prediction for regression. With each training, the parameters are continuously updated. This process continues until the output of the network is closer to the desired output.

The multi-layer architecture determines that data flows through more neurons and more parameters need to be trained, as well as more powerful than the single layer network. If the recurrent neural network is too deep beyond what is necessary, the computational cost will be expensive. Also, the phenomenon of overfitting may occur.

**Figure 4.** Architecture and workflow of ML-GRU.

#### *3.3. CEEMDAN-Based Multi-Layer Gated Recurrent Unit Networks (CEEMDAN-ML-GRU)*

In this study, we introduce a novel approach combining CEEMDAN and multi-layer GRU neural networks. We call this model CEEMDAN-ML-GRU for crude oil price forecasts. Figure 5 shows the architecture and workflow of the hybrid model. It aimed at improving the exiting crude oil price forecast techniques, which are less efficient or have poor accuracy in dealing with nonlinear and nonstationary regression tasks.

**Figure 5.** CEEMDAN-based Multi-layer Gated Recurrent Unit Networks (CEEMDAN-ML-GRU). This approach adopts a strategy of "divide and process", which transforms a complicated issue into Several simple issues, and processes them independently.

Firstly, CEEMD technology is used to add positive and negative paired white noise to the original exchange rate sequence, which overcomes the problem of large EEMD reconstruction errors and poor completeness of decomposition, and effectively improves the decomposition efficiency. Then the original exchange rate sequence is decomposed into IMFs based on different characteristics. Finally, the multi-layer LSTM based forecasting model is input for each component, and each of the IMF prediction results is superimposed to obtain the desired overall prediction.

As mentioned in the introduction, signal decomposition techniques such as wavelet transform, EMD and EEMD have been used for the analysis of time series of predicting energy prices. Machine learning, especially neural networks and deep learning methods, have also been applied to crude oil price forecasting to improve the learning process and prediction accuracy of crude oil price data. Signal decomposition technology is good at processing non-stationary data, and deep learning shows excellent performance when analyzing time series with nonlinear and long-term dependency characteristics. Some scholars have proved that integrating signal decomposition and deep learning methods for crude oil prices forecasting gained better results than only using a single method. In Section 1, we have mentioned that Wu et al. [24] proposed a novel model based on EEMD and long LSTM for crude oil price forecast. We will compare them experimentally with the proposed method in the next section. This kind of hybrid model is able to synthesize the strengths of each hybrid method, and significantly avoids the negative impact of the single method's inherent disadvantages on prediction performance.

The model first performs CEEMDAN on the original crude oil prices data. Then we input the previously decomposed IMF into the ML-GRU neural network. As for the output, it is meaningless to predict IMF alone, so we must use all the IMFs obtained by decomposing one sample at a time as the input of a single training or test to directly predict the price trend of crude oil. We use a one-step prediction strategy. By extracting the characteristics of IMFs, ML-GRU can use the crude oil prices trend data of the last p days to predict the price trend of the next day. CEEMDAN-ML-GRU prediction usually includes the following four main steps:

#### **4. Experiments**

#### *4.1. Datasets*

In this section, through a series of experiments, we verified the proposed CEEMDAN-ML-GRU model is more advanced than the state-of-the-art methods, which include the EEMD-LSTM. In order to measure and compare the forecasting performance of different models, WTI crude oil prices is employed for the sample data set for the experiment. This dataset source is the U.S. Energy Information Administration (EIA; http://www.eia.doe.gov/) and has 8321 observations that include all daily data from 2 January 1986 to 3 January 2019. A test set with 1664 observations is used to evaluate the prediction performance. The daily chart of WTI crude oil prices is illustrated in Figure 6.

**Figure 6.** The daily chart of WTI crude oil prices from 2 January 1986 to 3 January 2019.

#### *4.2. Evaluation Metrics and Baselines*

Following, we adopt three metrics including the root mean squared error (RMSE), mean absolute percent error (MAPE) and Diebold-Mariano (DM) to evaluate our model.

(1) Root mean squared error (RMSE):

$$\text{RMSE}(y, \hat{y}) = \sqrt{\frac{\sum\_{i=1}^{N} (y\_i - \hat{y}\_i)^2}{N}} \tag{22}$$

RMSE is one of the most common metrics and often used to measure the deviation between observations and prediction in the task of machine learning models. The letter y in the formula above denotes the observations, *y*ˆ is the prediction, and *N* indicate the number of samples.

(2) Mean absolute error (MAE):

$$\text{MAE}(y\_\prime \hat{y}) = \frac{\sum\_{i=1}^{N} |y\_i - \hat{y}\_i|}{N} \tag{23}$$

(3) MAE can better reflect the actual situation of the error of prediction. Mean Absolute Error Percentage (MAPE):

$$\text{MAPE}(y, \hat{y}) = \frac{1}{N} \sum\_{i=1}^{N} \left| \frac{y\_i - \hat{y}\_i}{y\_i} \right| \tag{24}$$

MAPE measures not only the absolute error between the predicted value and the observations, but also the relative distance. Unlike the previous two metrics, MAPE stands for percentage error, which can help it compare errors between different data sets.

(4) Diebold-Mariano Test (DM test):

$$\text{DM} = \frac{\overline{d}}{\sqrt{\left[\gamma\_0 + 2\Sigma\_{k=1}^{h-1}\gamma\_k\right]/n}} \tag{25}$$

It is used to test the statistical significance of the forecast accuracy of two forecast methods. Variable d is the subtraction of absolute error of the two methods. *d* is the mean of *di*. γ*k* is the autocovariance at lag *k*. DM ∼ N(0, 1), if the *p*-value > α, we conclude there is no significant difference between the two forecasts.

An ideal metric is not only reflected in the description of prediction accuracy, but also required to reflect the distribution characteristics of errors. Because of their effectiveness, these metrics are widely used for error measurement in regression or prediction tasks. The definition of each metric is different. When the sample size is large enough, RMSE is more reliable. RMSE is more suitable for Gaussian distribution error measurement than MAE, and MAE has a better performance in measuring uniform

error. A single metric is not enough to judge the pros and cons of the model. We need to combine multiple metrics to determine the pros and cons of the model [39].

To verify that the proposed method is advanced, we compare our method with a series of models, most of them mentioned in previous sections. These models include not only the common models such as: naïve forecasts(that the prediction is the same as the last period), ARRIMA, least squares SVR (LSSVR), ANN RNN, LSTM, but also some hybrid models such as: EEMD-ELM, EEMD-LSSVR, EEMD-ANN, EEMD-RNN, EEMD-LSTM and EEMD-SBL-ADD. Both EEMD-LSTM and EMD-SBL-ADD have been proposed in the last two years and represent the current state-of-the-art approaches for crude oil forecasting.

#### *4.3. Experimental Settings*

There is one input layer, two hidden layers of GRU and one dense layer in our proposed model. Figure 7 illustrate its Tensorflow computation graph which indicate the network architecture in our method. Each input sample is a matrix of *n* × *m*, which is represented by a NumPy array. '*n*' is the lagging order and '*m*' is the number of IMFs and residue. By trial and error, we determine set the number of hidden neurons to 32 and MSE as the loss function. The optimizer of training is adaptive moment estimation (Adam) which solve the problems of other algorithm, such as learning rate disappearing, convergence slowly or loss function fluctuating greatly. The learning rate in following experiment is set to 0.01. We adopt the strategy of one day ahead prediction to carry out our tasks. In other words, the prices of the past *n* days (*p*1, *p*2, ... *pn*−1, *pn*) is used to predict the price of the (*n* + 1)th day. Letter *n* is called the lag order which related to the size of neuron of the GRU. We adopt a strategy of grid search to determine the number of lagging order that is important for time series analysis. By trial and error, the lag order was set to 32.

**Figure 7.** TensorFlow computation graphs of the network structure of proposedmodel. This figureillustrates the data dependencies and control dependencies. The solid arrows indicate the data dependencies that show the flow of tensors between two ops, while dotted lines indicate the control dependencies.

For comparison purposes, we set the same initialization parameters for several deep learning models that participated in the experiment. For other parameters, since they have little effect on the results, we set the parameters to default.

The methods mentioned in the previous literature will be compared to our introduced approach in the performance of crude oil prices forecasts. These methods include not only the single model, LSSVR, ANN, ARIMA, and LSTM, but also the hybrid method EEMD-LSTM, which is considered the state-out-of- art method in crude oil prices forecasting.

In this study, all experiments are conducted in Python 3.7 via several specialized package, such as TensorFlow-GPU 2.0, Pyeemd 1.4, Keras 2.3 and so on. All experiments are conducted on a PC with a 3.2 GHz CPU, 16 GB RAM and 12 GB CG card.

#### *4.4. Experimental Process, Result and Analysis*

#### 4.4.1. Data Decomposition

The original WTI crude oil is decomposed to 13 components which include 12 IMFS and one CEEMDAN residue, where we added the white noise with standard deviation of 0.01 and the number of ensemble size is equal to the size of dataset. Then the IMFs and residue are splinted into the train set and test set. Both them are illustrated in Figure 8, in which the left part of dotted line is train set and the right part is test set.

By trial and error, we determine the parameters of ensemble size of 0.05 and noise strength of 100. It means that the white noise data with standard deviation of 0.05 and quantity of 100 will be added to the original data. Before data decomposition, we initialize the other parameters of the algorithm. Table 1 shows the names and descriptions of the main parameters. Table 2 shows the parameter settings.

**Table 1.** EEMD and CEEMDAN parameters description.

The time–frequency spectra of IMFs and the residue by CEEMDAN after decomposition are shown in Figure 8. We divide the original data and the decomposed result into a training set and a test set. The former accounts for 80% of the entire dataset, and the latter accounts for the remaining 20%. As shown in Figure 8, the data on the left side of the dotted line constitutes the training set, and the test set is on the right. As one kind of supervised learning algorithm, the input paired to the desired output in both the training and test sets.

**Figure 8.** The time–frequency spectrums of IMFs obtained by CEEMDAN. The top sub-picture shows the original WTI price data. The next 12 images represent intrinsic mode functions (IMF), which is listed in the order from the high to the low frequencies. The last component is the residue, which represents the portion of the original data not decomposed and the real trend of the original data.

#### 4.4.2. Training/Learning

After the model is build, we move on to the next step: training/learning. In this process, the training data is continuously feed into the model to incrementally improve the model's predictive performance. Both the loss function and optimizer of adaptive moment estimation (Adam) mentioned earlier are used to evaluate and optimize the model in training to achieve the purpose of model optimization. The workflow of training as shown in Figure 9.

**Figure 9.** Workflow of training.

The trend of the loss function curve during the iteration of model training from 1st to 100th epoch is shown in the Figure 10.

**Figure 10.** The loss curve graph of model.

The optimization algorithm mentioned above determines the optimization efficiency of the model, while loss represents the distance between the prediction and the observation. From the preceding graph, we observe that the loss curve descends rapidly in the initial epochs of training iteration, which shows that the model is optimized significantly by tuning the hyperparameter. However, between the 25th and 35th epochs of training, the loss curve is flat due to adaptive adjustment of the learning rate. After that, the curve continues to descend rapidly. The curve descending slowly means that

loss saturated after 70 epochs of iterations. After 100 training epochs, the loss (MSE) has descended to a very low level. In fact, we conducted 2500 epochs experiment for each model, with 120 samples per batch. Once the training\learning of the hyperparameters (weight) in the model is completed, the model can be used in the tasks of crude oil price forecasts. The experimental results are shown later in the paper.

#### 4.4.3. Predictive Performance of Different Single Model

We propose a hybrid model that includes two components, a decomposition algorithm CEEMDAN and a prediction model Multi-layer GRU. In order to make the method evaluation more accurate, we adopt a combination of independent evaluation and overall evaluation. So, we expand single model comparison and hybrid model comparison between two sets of experiments.

Here, we will evaluate our proposed CEEMD-ML-GRU's effectiveness for improving forecasting accuracy. The compared single models included the naïve forecasts, one classical time series method of ARIMA, two famous machine models of LSSVR and ANN and two popular deep learning model of LSTM and GRU. The samples used in each iteration of the deep learning model are randomly drawn, which makes the results of each prediction different. In order to improve the robustness of the model, experiments with each parameter condition were required to be trained 100 times. In order of performance, the metric at the median position represents the model. The results are shown in Table 3.

**Table 3.** Predictive performance comparison of single methods.

From Tables 3–5, we can see that:

**Table 4.** Comparison of GRU and LSTM model in the number of parameters.

**Table 5.** The Diebold–Mariano (DM) test results for single models on WTI crude oil prices.

4.4.4. Effect of Selecting Different Hybrid Approaches

On the basis of the previous single model prediction, we continue to evaluate the performance of the hybrid model based on the decomposition method. Table 6 shows the prediction performance of the corresponding hybrid models based on EEMD or CEEMDAN. Table 7 demonstrates that the result of DM test between CEEMDAN-ML-GRU and the other two models.

**Table 6.** The experimental results in terms of hybrid approaches on WTI crude oil prices forecasting.

**Table 7.** The Diebold–Mariano (DM) test results for hybrid models on WTI crude oil prices.

In this set of experiments, we added a prediction model based on the EEMD decomposition algorithm, and these hybrid models have appeared in the latest literature. Each single model combines CEEMDAN and EEMD respectively, and thus two sets of mixed methods will be obtained. Subsequently, these models were tested in the WTI crude oil price prediction task to derive who is the optimal model. Table 6 demonstrates the predictions of different hybrid models. From the Tables 6 and 7, we can observe that:

(1) Both EEMD and CEEMDAN plus prediction models outperform single model significantly in this prediction task. It indicates that machine learning models with signal processing algorithm based contribute to the better forecasting performance in the nonlinear and non-stationary time series analysis.

#### **5. Conclusions**

In this study, a hybrid model called CEEMDAN-ML-GRU for crude oil price forecasting is proposed. This model takes full use of the advantages of the signal processing algorithm CEEMDAN and the multi-layer gated recurrent unit networks (ML-GRU). As mentioned in the previous section, the hybrid model uses CEEMDAN to solve the non-stationarity problem of crude oil price data, and generalizes the nonlinear crude oil prices data by a multi-layered GRU neural network. We conduct a large number of experiments to verify the effect of the proposed method in forecasting task by using the WTI price data as sample data. The experimental results show that our proposed method goes beyond other traditional statistical methods, machine intelligent algorithms and other hybrid models, which include the EEMD-LSTM method proposed in 2019.

In addition to crude oil price forecasts, the introduced CEEMDAN-ML-GRU model can also be extended to solve other complex problems in other areas, such as time series forecasts or risk measurements in financial markets. The main purpose of this approach is to improve the accuracy of short-term crude oil price predictions and help decision makers minimize the risks of the crude oil market. However, the proposed method is mainly applied to short-term forecasts, so only daily data is used. If we need to predict long-term price trends, we need to combine this method with economic theory or measurement methods to play a greater role. This is exactly the research plan that we will follow in our future research.

**Author Contributions:** Conceptualization, H.L.; Formal analysis, Q.S.; Investigation, Q.S.; Methodology, H.L.; Project administration, H.L. and Q.S.; Supervision, Q.S.; Writing—original draft, H.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the provincial social science fund of Fujian (No. FJ2018B056), A Study on the Supervision of Taiwan-funded Banks in the Fujian Free Trade Zone, and the National Natural Science Foundation of China (No. 71872046), Environmental Regulation, Enterprise Innovation and Industrial Transformation and Upgrade: Theory and Practice Based on Energy Saving and Emissions Supervision.

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

#### **References**

© 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*

### **Optimizing Predictor Variables in Artificial Neural Networks When Forecasting Raw Material Prices for Energy Production**

**Marta Matyjaszek 1, Gregorio Fidalgo Valverde 2 , Alicja Krzemie ´n 3, Krzysztof Wodarski 4 and Pedro Riesgo Fernández 2,\***

Received: 26 March 2020; Accepted: 16 April 2020; Published: 18 April 2020

**Abstract:** This paper applies a heuristic approach to optimize the predictor variables in artificial neural networks when forecasting raw material prices for energy production (coking coal, natural gas, crude oil and coal) to achieve a better forecast. Two goals are (1) to determine the optimum number of time-delayed terms or past values forming the lagged variables and (2) to improve the forecast accuracy by adding intrinsic signals to the lagged variables. The conclusions clearly are in opposition to the actual scientific literature: when addressing the lagged variable size, the results do not confirm relationships among their size, representativeness and estimation accuracy. It is also possible to verify an important effect of the results on the lagged variable size. Finally, adding the order in the time series of the lagged variables to form the predictor variables improves the forecast accuracy in most cases.

**Keywords:** raw material; price forecasting; artificial neural network; predictor variable; lagged variable size; rolling window; coking coal; natural gas; crude oil; coal

#### **1. Introduction**

Artificial neural networks (ANN) have been widely used as accurate forecast aids addressing issues directly or indirectly related with energy or raw materials: energy production [1], raw material inventory levels [2], crude oil prices [3], volatility of stock price indices [4], electricity prices [5], stock prices [6], gold prices [7], copper spot prices [8], off-gases production [9], currency exchange rates [10], etc.

This paper analyzes the forecast of raw material prices for energy production by means of ANN, focusing on the selection of optimum parameters in order to configure the ANN. Design of experiments (DOE) is normally used to select these parameters [11].

DOE can be focused on estimating the number of neurons in hidden layers [12], on forming training and test datasets [13], on eliminating redundant dimensions in the predictor variables trying to achieve compression [14,15], on determining the optimum size of the lagged variables [16], on adding signals to the lagged variables [17], etc.

A heuristic approach will be used to optimize the predictor variables in ANN by means of (1) modifying the lagged variable size and (2) adding intrinsic signals to the lagged variables.

Addressing the lagged variable size, Liu and Su [18] indicated that a larger size allows for increasing the forecast accuracy, while it decreases the representativeness of the subsample heterogeneity. Moreover, although smaller sizes may improve representativeness, they will reduce the estimation accuracy. In the above empirical research, they used lagged variable sizes of 12, 24 and 36 months to test different alternatives. Nevertheless, the lagged variable sizes indicate very little effect on the results.

Tang and Abosedra [19] established that the lagged variable regression results are very sensitive regarding their size, but as there are no proper methods to select an optimum size, arbitrary selections have to be made. Other authors argue that a larger lagged variable size would lead to short-run predictability information being missed, and thus, a shorter size is preferred [20,21].

WEKA from the Machine Learning Group at the University of Waikato (Waikato, New Zealand) [22,23], a well-known open source machine learning software widely used for teaching, research and industrial applications, has a specific time-series analysis environment to forecast models. WEKA's time-series framework uses a machine learning/data mining approach to model time series. It transforms the data by removing the temporal ordering of individual input examples by encoding the time dependency via additional input fields or lagged variables. When using WEKA, it is possible to manipulate and control how lagged variables are created. They are the main mechanism to capture the relationship between current and past values, creating a window over a certain time period. Essentially, the number of lagged variables created determines the size of the window.

Regarding the adding of intrinsic signals to the lagged variables, Tavakoli et al. [24] proposed an input management system based on flexible data by immediately providing a variable definition layer on top of the acquisition layer to feed a data mining module to build modeling functions. Recently, and within the neural networks field, Uykan and Koivo [25] have presented and analyzed a new design for the predictor variables of a radial basis function neural network. In this design, the predictor variables were augmented with a desired output vector, allowing for better/comparable performance when compared with the standard neural network.

Raw material selection, namely, coking coal, natural gas, crude oil and coal, was based on the representativeness and price availability of such materials. In the case of coking coal, prices were obtained from the Colombian Mining Information System as they were publicly disclosed, while for the rest of the raw materials, their prices were obtained from the World Bank Commodity Price Data (The Pink Sheet) under a Creative Commons Attribution 4.0 International License (CC BY 4.0).

The program used to simulate ANNs was NeuralTools 7.5 from Palisade Corporation (Ithaca, NY, USA).

#### **2. Method**

This paper will attempt to improve the result of the time-series forecasting of raw material prices for energy production developed with ANN by means of a twofold optimization of the predictor variables (modifying the lagged variable size and adding intrinsic signals to the lagged variables), taking into consideration the previous work of Matyjaszek et al. [26], in which coking coal prices were forecasted by means of autoregressive integrated moving average models (ARIMA) [27,28] and ANN, as well as the transgenic time-series theory.

#### *2.1. Artificial Neural Networks*

Two different types of ANNs proposed by Specht [29] will be tested using the best net search function that is available in NeuralTools: generalized regression neural networks (GRNNs), which were used in the past to forecast European thermal coal spot prices among very different applications [30,31], and multilayer feedforward networks (MLFNs), with one or two layers as described in García Nieto et al. [32].

GRNNs are based on nonlinear regression theory and are very closely related to probabilistic neural nets. In GRNNs, a case prediction with a dependent value that is unknown is obtained by means of interpolation from the training cases, with neighboring cases given more weight [33]. The optimal parameters for the interpolation are found during training. The main advantage is not requiring any configuration at all.

Figure 1 presents a GRNN with two independent variables in the input layer and only four training cases, with the pattern layer having four nodes corresponding to these training cases. Each of the nodes will compute its Euclidean distance regarding the presented case. Then, these values pass to the summation layer. The summation layer has two parts: one is the numerator, and the other is the denominator. The numerator contains the sum of multiplying the training output data and the activation function. The denominator is the sum of all the activation functions.

**Figure 1.** Configuration of a generalized regression neural network (GRNN) with two independent numeric variables and four training cases.

Finally, the output layer has only one neuron that calculates the output by dividing the numerator and the denominator of the summation layer.

On the other hand, MLFNs consist of an input layer, one or even two hidden layers and the output layer. A MLFN is configured by specifying the number nodes in the hidden layers. The net behavior will depend on the number of nodes selected for each hidden layer, the connections weights, the bias terms that are assigned to each node, and the activation/transfer function selected to convert into output the inputs of each node. They are able to approximate complex relationships between the variables.

#### *2.2. Lagged variable Size*

One of the issues to be analyzed within this paper is the current discussion about if a larger lagged variable size allows for increasing the forecast accuracy [18] or if a shorter size is preferred based on avoiding to miss short-run predictability information [20,21]. Another issue will be whether lagged variable regression results are very sensitive regarding their size [19] or not [18].

Lagged variables are generated by a number of linear time-delayed input terms or past values, normally in ascending order, such as P(*t*−*n*) ... P(*t*−*2*), P(*t*−*1*)*,* to estimate the output value P(*t*) [34]. They are also referred to as rolling windows [35].

To undergo a first estimation of the number of time-delayed terms that should form the lagged variables (*n*), there are several alternatives that can be selected, including the one developed by Ren et al. [36], which uses the seasonal characteristic that appears in the autocorrelation function (ACF) plot, although this value is not always available.

Other common approach is to approximate the value by determining the square root of the amount of data available to undertake the analysis [26]:

$$\sqrt{\text{Total } n^{\text{\textquotedblleft} of \text{ data}}} \tag{1}$$

Another alternative is the one used also by Matyjaszek et al. [26], in which the adequate number of time-delayed terms, *k*, that should form each input layer is calculated as follows:

$$\text{Total } n \text{\textquotedblleft of data } \le n^2 + 2n + 1,\tag{2}$$

$$m = 1 + k + 1\tag{3}$$

#### *2.3. Adding Intrinsic Signals to the Lagged Variables*

Up to date and in order to improve the forecast accuracy by adding signals to the lagged variables, research is focused on the extrinsic ones [24,25]. This paper would analyze whether it is feasible to optimize predictor variables by adding intrinsic signals, so that the neural network will have more information available; thus, a better forecast could be made. For this purpose, the order in the time series of each lagged variable would be used, so the ANN could exploit this feature.

This line of thinking is congruent with the work developed by Barabási [37], who states that there is a huge disconnect between network science and deep learning; although ANN are abstractions of natural processes, some of the key neural networks could not be more ignorant about real networks. Main deep learning algorithms treat network features, like degree, as simple variables. Thus, they cannot truly exploit the network effects, which are the essence of these systems, as in networked systems the key information is in the relationships between the connected components (i.e., in the links or edges, which are the direct interactions between nodes), not in the node attributes.

Table 1 presents an example of the first through tenth lagged variables used in a model with five time-delayed input terms: P(*t*−*k*) ... P(*t*−*2*), P(*t*−*1*), with *k* = 5, as well as the output to be estimated: P(*t*).

**Table 1.** First through tenth lagged variables with five time-delayed input terms, and the output to be estimated (*t*).

Table 2 presents the same first through tenth lagged variables with five time-delayed input terms plus the order in the time series of each lagged variable, as well as the output to be estimated: P(*t*).

**Table 2.** First through tenth predictor variables with 5 time-delayed input terms plus the order in the time series of each lagged variable, and the output to be estimated (*t*).

#### *2.4. Figures of Merit*

The experimental results will be evaluated using the two most common figures of merit [38], namely, the root mean squared error (*RMSE*) and the mean absolute error (*MAE*).

The *RMSE* is an excellent general-purpose error measure used for numerical predictions. It amplifies and penalizes large errors and can be expressed as follows:

$$RMSE = \sqrt{\frac{\sum\_{t=1}^{n} \left(A\_t - F\_t\right)^2}{n}},\tag{4}$$

where *At* is the actual value, *Ft* is the forecasted value, and *n* is the number of forecasted values.

The *MAE* is used to measure how close the predictions are to the outcomes and can be expressed as follows:

$$MAE = \frac{1}{n} \sum\_{t=1}^{n} |A\_t - F\_t| \tag{5}$$

Chai and Draxler [39] proposed the use of a combination of metrics including but not limited to the *RMSE* and the *MAE*. Conversely, Carta et al. [40], when addressing wind resource prediction, proposed using the *MAE*, the *MAPE* and the index of agreement (*IoA*).

In this paper, the standard deviation of absolute error (*STD* of *AE*) was selected to complement these measures as in Lazaridis [38], characterizing the dispersion of the absolute errors.

#### **3. Results**

#### *3.1. Coking Coal*

The dataset used was the Colombia hard coking coal monthly prices free on board (FOB) for the period from January 1991 to December 2015, as publicly disclosed by the Colombian Mining Information System [41], totaling 300 data points. The dataset is presented in Figure 2.

**Figure 2.** Colombia hard coking coal monthly prices FOB from January 1991 to December 2015 [41].

In first place, GRNNs and MLFNs with 2–6 nodes in the first hidden layer (as the second layer is seldom needed for better prediction accuracy) were tested using the best net search function that is available in NeuralTools.

Table 3 presents the results of this best net search, where the configuration of the lagged variables was made using 19 time-delayed terms, which was the lagged variable size given by the seasonal characteristic that appears in the autocorrelation function (ACF) plot of the transformed time series when representing a consistent genome [26,42]. The results showed that the GRNN improved all the MLFN models, as stated for most cases by Modaresi et al. [43], so this will be the ANN to be used in this paper.

**Table 3.** Best net search.

To undergo a first estimation of the number of time-delayed terms that should form the lagged variables, the square root of the amount of data available to undertake the analysis was calculated:

$$\sqrt{\text{Total } n^{\text{\textquotedblleft} of data}} = \sqrt{300} = 17.32$$

Using the other alternative previously mentioned, the value obtained is as follows:

$$300 \le n^2 + n + 1 = \text{\(n = 16.32\)}$$

$$k = n - 2 = 14.32$$

Nevertheless, the GRNN will be trained starting with 12 time-delayed input terms and up to 24 time-delayed input terms, a range that includes the previous values of *k* as well as one and two complete year periods [18].

This allows considering almost any periodical aspect that may be hidden within the time-series values but without drastically reducing the sample size requirements according to Turmon and Fine [44]. The results from the GRNN training are presented in Table 4.

**Table 4.** Training results for the GRNN model of the Colombian coking coal time series, from 12 to 24 time-delayed input terms (total number of data points is 300).

The figures in bold correspond to the model that achieves better performance measures.

Based on these measures, the best result was obtained with 19 time-delayed input terms, which was the lagged variable size given by the seasonal characteristic that appears in the autocorrelation function (ACF) plot of the transformed time series when representing a consistent genome [26].

The figures of merit were root mean squared error (*RMSE*) of 4.550, mean absolute error (*MAE*) of 2.784 and standard deviation of absolute error of 3.599. With 30% tolerance, the percentage of bad predictions was 8.5409%.

Then, it was checked if it was feasible to optimize the predictor variables by adding intrinsic signals to the lagged variables so that the ANN would have more information available, and thus, a better forecast could be achieved.

For this purpose, the order of the lagged variables in the time series was considered. Using these predictor variables, results from the training of the GRNN are presented in Table 5.

**Table 5.** Training results for the GRNN model of the Colombian coking coal time series, from 12 to 24 time-delayed input terms and including the order in the time series of the lagged variables (total number of data points is 300).

The figures in bold correspond to the model that achieves better performance measures.

The best result was obtained with 24 time-delayed input terms, with a *RMSE* of 3.376, a *MAE* of 2.286 and a standard deviation of absolute error of 2.485. With 30% tolerance, the percentage of bad predictions was 2.5362%. Thus, the order in the time series of the lagged variables clearly improved the model's forecasting performance, as it significantly reduced the *RMSE*, the *MAE*, the standard deviation of absolute error and the percentage of bad predictions.

#### *3.2. Natural Gas*

The second raw material for energy production analyzed was natural gas. The dataset used was natural gas prices in Europe for the period from January 1991 to August 2019, totaling 344 values.

Prices were obtained from the World Bank [45] and are presented in Figure 3 in MMBtu, also known as million British thermal units, with 1 MMBtu = 28.263682 m3 of natural gas at 1 ◦F.

**Figure 3.** Natural gas prices in Europe for the period from January 1991 to August 2019 [45].

As no seasonal characteristic appears in the autocorrelation function (ACF), in order to estimate the number of time-delayed input terms, the following calculations were made:

$$\sqrt{\text{Total } n^2 \text{ of data}} = \sqrt{344} = 18.55\_\circ$$

$$344 \le n^2 + n + 1 = > n = 17.54$$

$$k = n - 2 = 15.54$$

Nevertheless, the GRNN was trained with lagged variables starting with 12 time-delayed input terms and up to 24 time-delayed input terms, using the same interval as in the coking coal case and for the same reasons.

The results from the training of the GRNN are presented in Table 6.

**Table 6.** Training results for the GRNN model of the European natural gas time series, from 12 to 24 time-delayed input terms (total number of data points is 344).

The figures in bold correspond to the model that achieves better performance measures.

The best result was obtained with 12 time-delayed input terms, with a *RMSE* of 0.06273, a *MAE* of 0.03456 and a standard deviation of absolute error of 0.05235. With 5% tolerance, the percentage of bad predictions was 4.8193%. A 5% tolerance was used this time, as with a 30% tolerance the percentages of bad predictions were always zero.

Then, the GRNN was trained using the same number of time-delayed input terms but considering the order in the time series of each lagged variable. The results from the training of the GRNN are presented in Table 7. In this case, the best result was obtained with 14 time-delayed input terms, with a *RMSE* of 0.01970, a *MAE* of 0.01032 and a standard deviation of absolute error of 0.01678. With 5% tolerance, the percentage of bad predictions decreased to zero.

**Table 7.** Training results for the GRNN model of the European natural gas time series, from 12 to 24 time-delayed input terms and including the order in the time series of the lagged variables (total number of data points is 344).

The figures in bold correspond to the model that achieves better performance measures.

Thus, again, the order in the time series of the lagged variables clearly improved the model's forecasting performance, as it significantly reduced the *RMSE*, the *MAE*, the standard deviation of absolute error and the percentage of bad predictions.

#### *3.3. Crude Oil*

The third raw material for energy production analyzed was crude oil. The dataset used was that of Brent crude oil prices for the same period as for the natural gas: January 1991 to August 2019, totaling 344 values. Again, prices were obtained from the World Bank [45] and are presented in Figure 4 in \$/bbl, that is, dollars per barrel, with 1 barrel being approximately 159 liters.

**Figure 4.** Brent crude oil prices for the period from January 1991 to August 2019 [45].

The number of estimated time-delayed input terms that should be used were the same as in the case of natural gas and crude oil, as the number of data points was the same in all cases (344).

The GRNN was trained again with lagged variables starting with 12 time-delayed input terms and up to 24 time-delayed input terms. The results from the training of the GRNN are presented in Table 8.

**Table 8.** Training results of the GRNN model for the Brent crude oil time series, from 12 to 24 time-delayed input terms (total number of data points is 344).

The figures in bold correspond to the model that achieves better performance measures.

The best result was obtained with 23 time-delayed input terms, with a *RMSE* of 1.177, a *MAE* of 0.6645 and a standard deviation of absolute error of 0.971. With 5% tolerance, the percentage of bad predictions was 23.0530%.

Then, the GRNN was trained using the same number of time-delayed input terms but considering the order in the time series of each lagged variable. The results from the training of the GRNN are presented in Table 9.

**Table 9.** Training results of the GRNN model for the Brent crude oil time series, from 12 to 24 time-delayed input terms including the order in the time series of the lagged variables (total number of data points is 344).

The figures in bold correspond to the model that achieves better performance measures.

The best result was obtained with 14 time-delayed input terms, with a *RMSE* of 0.0281, a *MAE* of 0.0178 and a standard deviation of absolute error of 0.0218. With 5% tolerance, the percentage of bad predictions decreased again to zero.

With 12 time-delayed input terms, it is clear that the ANN was able to learn the exact configuration of the time series, but only in this case.

Thus, the order in the time series of the lagged variables clearly improved the model's forecasting performance, as it significantly reduced the *RMSE*, the *MAE*, the standard deviation of absolute error and the percentage of bad predictions.

#### *3.4. Coal*

The fourth and last raw material for energy production analyzed was coal. The dataset used was that of Australian coal prices for the same period as crude oil and natural gas: January 1991 to August 2019. Prices were also obtained from the World Bank [45] and are presented in Figure 5 in \$/t.

**Figure 5.** Australian coal prices for the period January from 1991 to August 2019 [45].

The number of estimated time-delayed input terms that should be used were the same as in the case of natural gas and crude oil, as the number of data points was the same in all cases (344).

The GRNN was trained again with lagged variables starting with 12 time-delayed input terms and up to 24 time-delayed input terms, as in the case of European natural gas and Brent crude oil. The results from the training of the GRNN are presented in Table 10.

**Table 10.** Training results of the GRNN model for the Australian coal time series, from 12 to 24 time-delayed input terms (total number of data points is 344).

The figures in bold correspond to the model that achieves better performance measures.

The best result was obtained with 19 time-delayed input terms, with a *RMSE* of 1.3280, a *MAE* of 0.7397 and a standard deviation of absolute error of 1.1029. With 5% tolerance, the percentage of bad predictions was 15.6923%.

Then, the GRNN was trained using the same number of time-delayed input terms but considering the order in the time series of each lagged variable. The results from the training of the GRNN are presented in Table 11.

**Table 11.** Training results for the GRNN models of the Australian coal time series, from 12 to 24 time-delayed input terms and including the order in the time series of the lagged variables (total number of data points is 344).

The figures in bold correspond to the model that achieves better performance measures.

The best result was obtained with 22 time-delayed input terms, with a *RMSE* of 1.3678, a *MAE* of 0.8185 and a standard deviation of absolute error of 1.0959. With 5% tolerance, the percentage of bad predictions was 16.7702%.

In this case, adding the order in the time series of each lagged variable only improved the standard deviation of absolute error.

Nevertheless, if the *RMSE* and *MAE* were compared with the training results obtained without adding the order in the time series of each lagged variable, although slightly higher, they were very similar to the best results, and better than the rest of the training results.

Thus, the differences between both options were almost negligible.

#### **4. Discussion and Conclusions**

This paper applied a heuristic approach to optimize the predictor variables in artificial neural networks when forecasting raw material prices for energy production to achieve a better forecast.

Two goals are (1) to determine the optimum number of time-delayed terms or past values forming the lagged variables and (2) to optimize predictor variables by adding intrinsic signals to the lagged variables.

The experimental results were evaluated using the two most common figures of merit, the root mean squared error (*RMSE*) and the mean absolute error (*MAE*), as well as the standard deviation of absolute error, as the scientific literature proposes the use of a combination of metrics including *RMSE* and *MAE*, but not being limited to them.

Results demonstrated, first, that in opposition to scientific literature when addressing lagged variable size, a larger size did not allow for increasing the forecast accuracy, and that smaller sizes did not reduce the estimation accuracy. Moreover, the lagged variable regression results were very sensitive regarding their size.

In the three raw materials with the same number of cases (natural gas, crude oil and coal), the best results were obtained with rolling window sizes of 12, 23 and 19, respectively. Furthermore, it was possible to verify an important effect of the lagged variable size on the results, with differences that in some cases were larger than 20%.

Thus, and in opposition again to scientific literature indicating that there are no proper methods to select an optimum size so arbitrary selections have to be made, it is recommendable to address this question by trial and error method, although the approximate size can be estimated in order to select the complete year's period range to which this value belongs, e.g., 12–24 months or 24–36 months. This will allow considering any periodical aspect that may be hidden within the time series values, but without drastically reducing or increasing the sample size requirements for neural networks.

Second, in three of the four raw materials analyzed (coking coal, natural gas and crude oil), it was possible to improve the forecast accuracy by adding the order in the time series of the lagged variables to form the predictor variables. The best results were achieved with rolling window sizes of 24, 14 and 14, respectively.

In the case of the Australian coal, this process only improved the standard deviation of absolute error. Nevertheless, if the *RMSE* and *MAE* were compared with the training results obtained without adding the order, although slightly higher, they were very similar to the best results and better than the rest of the training results without adding the order.

As the differences between both options were almost negligible, it is possible to recommend adding the order in the time series of each lagged variable to the predictor variable in all cases.

Third, only with the Brent crude oil with 12 time-delayed input terms and considering the order in the time series of the lagged variables, the ANN was able to learn or deduct the exact configuration of the time series. This is completely congruent with the fact that there is a huge disconnect between network science and deep learning, as the key information is in the relationships between the connected components, not in the node attributes.

Concluding, the findings presented in this paper have an immediate practical application addressing the forecast of time series by means of ANN that consider lagged variables, without being restricted to the studied case of raw material prices for energy production.

Any forecast may be optimized just by adding an intrinsic signal to the predictor variable consisting of the order in the time series of each lagged variable. By doing this way, the ANN will be able to exploit this feature, something that will not happen otherwise. In most of the cases, figures of merit may improve (may be reduced) up to a 20%, with the consequent benefit for decision-makers regarding savings, efficiency/benefit gains and/or lower risk.

Regarding the size of the lagged variable, a selection should be made about the period that will be analyzed in order to undergo a trial and error process. This selection should follow the procedure shown in this paper or other ones that may be found in the scientific literature.

Further research should address different issues such as the use of more intrinsic signals. Regarding this issue, authors have made interesting preliminary approaches by considering the transgenic time series theory that allows eliminating anomalous phenomena from the time series. Augmenting the lagged variables within this anomalous period with a '1' and the rest with a '0', or vice versa, it was possible to improve a priori the figures of merit.

Another area of interesting future research will be to develop a procedure to determine accurately the number of time-delayed input terms that should be used when considering the order in the time series of the lagged variables. While the seasonal characteristic that appears in the autocorrelation function (ACF) plot is valid before augmenting the lagged variables, later this figure is no longer valid, so a new approach should be addressed. Nevertheless, nothing is yet developed addressing the time series with an ACF plot that does not allow one to extract a seasonal characteristic. Again, the transgenic time series theory could be of help in these cases.

Finally, it should be addressed by future research why in the case of Australian coal, or in similar cases, it was not possible to improve the figures of merit by adding to the predictor variable the order in the time series of each lagged variable.

**Author Contributions:** Conceptualization, M.M. and P.R.F.; methodology, A.K. and P.R.F.; software, G.F.V.; validation, K.W.; formal analysis, G.F.V.; investigation, A.K., P.R.F. and M.M.; resources, G.F.V.; data curation, G.F.V.; writing—original draft preparation, M.M. and P.R.F.; writing—review and editing, A.K.; visualization, P.R.F.; supervision, K.W. and A.K.; and project administration, P.R.F. and K.W. All authors have read and agree to the published version of the manuscript.

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

**Acknowledgments:** The doctoral research fellowship from the Technological Research Foundation "Luis Fernández Velasco" in Mineral and Mining Economics at the School of Mining, Energy and Material Engineering from the University of Oviedo (Spain) for the first author is highly appreciated.

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

#### **References**

© 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* **Prediction of Health-Related Leave Days among Workers in the Energy Sector by Means of Genetic Algorithms**

**Aroa González Fuentes 1, Nélida M. Busto Serrano 2,\*, Fernando Sánchez Lasheras 3,\* , Gregorio Fidalgo Valverde 4 and Ana Suárez Sánchez 4**

Received: 3 April 2020; Accepted: 12 May 2020; Published: 14 May 2020

**Abstract:** In this research, a model is proposed for predicting the number of days absent from work due to sick or health-related leave among workers in the industry sector, according to ergonomic, social and work-related factors. It employs selected microdata from the Sixth European Working Conditions Survey (EWCS) and combines a genetic algorithm with Multivariate Adaptive Regression Splines (MARS). The most relevant explanatory variables identified by the model can be included in the following categories: ergonomics, psychosocial factors, working conditions and personal data and physiological characteristics. These categories are interrelated, and it is difficult to establish boundaries between them. Any managing program has to act on factors that affect the employees' general health status, process design, workplace environment, ergonomics and psychosocial working context, among others, to achieve success. This has an extensive field of application in the energy sector.

 Introduction**

Over the past few years, the main concerns of industry, especially in developed countries, have been to improve the workers' productivity, occupational health and safety in the workplace, physical and mental well-being, and job satisfaction. Through the application of ergonomics, it has been shown that these issues have improved, so an effective implementation of ergonomics in the workplace can achieve a balance between the characteristics of the worker and the demands of the task, in addition to improving workplace design and introducing appropriate management programs. Companies that belong to the energy sector have also been working in this direction, showing some distinctive features that have been identified and studied here.

When the industry does not get involved in the abovementioned issues, it can affect the lives of workers. This results in a risk of deterioration in health and causes absenteeism. Currently, one of the major concerns is sick leave.

 In this introduction, we first revise the studies that explain the behavior of such factors in the general working environment. Afterwards, we focus on the few works that are specifically oriented to the singularities of the energy sector.

According to the EUROSTAT [1], some 1194 M€ have been spent on sickness and health benefits in the European Union. This number was equivalent to 8.0% of gross domestic product (GDP). The average expenditure per inhabitant was 2338 Euros. In recent years, the spending on sick leave benefits has increased to 12.4%, with Norway being the country with the highest spending (32.4%) compared to Portugal, the country with the lowest (5.8%). For these reasons, the problem of absenteeism is of great interest to healthcare professionals, employers and economists.

According to the latest research published by the European Commission [2], women have higher rates of sick leave than men. There are multiple reasons for this. They have more precarious work and work contracts often linked to low income. Moreover, women often seek medical help for less-serious illnesses and are more frequently diagnosed with mental-health-related illnesses. Another explanation is connected to the burden of housework and childcare.

The aforementioned study [2] also mentions that sick leave increases with age; elderly people take longer-term leave compared to young people, who generally take short-term leave. This is due to the fact that health worsens with age, and working conditions have deteriorated since the economic recession. A clear correlation between occupational, socioeconomic status and absence due to illness is highlighted: the more physically demanding the occupation is and the lower the socioeconomic status, the higher the absenteeism.

Regarding the energy sector, different factors have been studied that may lead to sick leave. A study carried out in the petroleum industry [3] corroborates that women are more likely to take sick leave than men. This may be mainly due to three factors: there is sick leave due to pregnancy; they tend to have more temporary contracts; and they suffer more psychological problems. Another risk factor is smoking; workers who are smokers or former smokers have a higher risk of taking sick leave than non-smokers, those who consume alcohol and even workers exposed to chemical products. Other good predictors of absenteeism due to illness are abnormal sleep and job dissatisfaction. Some physical activity is recommended for workers.

Another specific factor that has been studied in the energy sector is how the different shiftwork patterns can cause sick leave. The employees of a power plant prefer 12-hour shifts rather than eight-hour shifts, since they enjoy longer breaks and an improved social and domestic life. This measure also improves mood, health status and both the quality and quantity of sleep. The only drawback is that it can pose a potential safety risk for the employee when performing highly demanding tasks at the end of the shift, since concentration decreases [4]. A study conducted among workers of a nuclear power plant [5] confirms all these findings on the relation between working in shifts and domestic, social life and well-being of the worker.

Another relevant factor in the energy sector is the type of work carried out. It has been proven that people who do manual work, also called blue-collar (production) workers, are more likely to remain on long-term leave, and this would be of longer duration than in those people who perform skilled jobs, also known as white-collar workers (office workers and managers). The latter are exposed to a greater risk for short-term pain, mainly musculoskeletal disorders (MSD), but this could be avoided by correcting posture [6].

The design of the job is another very important factor to take into account in this sector. A study of workers in a thermal power plant [7] detected deficiencies both in its facilities and its resources. This type of industry is complex and usually has more problems than other industries. Apart from health problems due to ergonomic factors, it has been found that production work in combination with bad environmental conditions (excessive temperatures in summer) and very noisy and dusty environments are factors that tend to worsen the health of workers and increase the possibility of their going on sick leave.

In summary, talking about the energy sector, it seems that absenteeism rates have different behaviors attending to factors like gender, working organization circumstances, such as shifts or psychosocial demands, workplace environment and other factors related with workstation design. The way these factors combined affect sick leave remains unexplored.

The use of machine-learning techniques to predict occupational health and safety outcomes in different fields is not new, whether focused on work-related accidents [8], fire risk [9,10], MSD [11–14] or visual disorders [15,16]. Some of these works focus on specific sectors, such as mining or the health industry.

Nevertheless, to date, no researches have studied how the combination of factors such as age, gender, well-being, domestic and social life, as well as psychosocial factors, can influence the proneness to sick leave among workers in the energy sector. As far as it is known by the authors, most of the previous research in this field that make use of machine-learning techniques employed just only, for example, support vector machines [8,16], artificial neural networks [9,11], Multivariate Adaptive Regression Splines (MARS) [10,14] or k-nearest neighbors [12]. All these research studies have shown the utility of machine-learning techniques in this area for both regression and classification. However, until today, there have been few works [13,15] that combine more than one machine-learning methodology in order to improve their performance.

In this research, a hybrid methodology that combines MARS and genetic algorithm is proposed for predicting the number of days absent from work due to sick or health-related leave, among workers in the industry sector, according to ergonomic, social and work-related factors reported in the Sixth European Working Conditions Survey (EWCS).

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

#### *2.1. Dataset*

This research work employs selected microdata from the Sixth European Working Conditions Survey (EWCS), which was conducted in 2015 by the European Foundation for the Improvement of Living and Working Conditions, Eurofound [17]. The EWCS is generally conducted every five years, providing an overview of working conditions of the European population. A random representative sample of "persons in employment" (i.e., employees and the self-employed) is surveyed through a questionnaire administered face-to-face. The Eurofound datasets are stored and promoted online by the UK Data Service [18]. Upon request, the data are available free of charge, provided they will be used for non-commercial purposes.

Almost 44,000 workers in 35 countries were interviewed through the sixth wave of the EWCS. The validity of the questionnaire was guaranteed by a questionnaire-development group composed of experts and representatives of the European Commission and different international organizations [19]. This sixth edition codified more than 370 variables that included physical and psychosocial risk factors, working time, place of work, work-pace determinants, employee participation, job security, social relations, personal conditions, etc. The whole list of variables included in the sixth wave of the survey can be found in the source questionnaire, available online at the website of Eurofound [20]: https://www.eurofound.europa.eu/sites/default/files/page/field\_ef\_documents/6th\_ ewcs\_2015\_final\_source\_master\_questionnaire.pdf.

The size of the initial dataset was first reduced by only selecting workers from energy-related sectors. The final sample consisted of 420 workers (333 men and 87 women), aged between 17 and 71 years (average 44; see Figure 1) from the following NACE Revision 1 sections: mining of coal and lignite; extraction of peat; extraction of crude petroleum and gas; mining of uranium and thorium ores; manufacture of coke, refined petroleum products and nuclear fuel; electricity, gas, steam and hot water supply. The distribution of the subjects by country is shown in Table 1. Table 2 presents their level of studies, according to the International Standard Classification of Education (ISCED). Table 3 shows the distribution of the sample of workers according to their household's total monthly income. The average leave time was of 5.9 days, with a standard deviation of 16.1 days. Only two workers have a leave longer than 100 days. Please also note that leaves over 10 days represent only 18.33% of the total.

**Figure 1.** Age distribution of the workers.

**Table 1.** Distribution of the sample of workers by country.

\* Other countries: FYROM, Luxembourg, Slovakia, Hungary, Netherlands, Sweden, Switzerland, Greece, Turkey, Italy, Finland, Cyprus, Portugal, Malta.

A second-dimensional reduction was carried out by decreasing the number of variables through expert criteria. Only the 59 most relevant independent variables were preselected to initially feed the model developed and to try to explain the output variable. Some of the variables were designed as Likert scales, some of them were binary and a few were continuous (numerical). Please note that it would have been possible to perform this reduction also by means of either genetic algorithms or other methodologies like decision trees or PCA, but in our understanding, it was less time-consuming to use expert criteria. Please note that this is a good way in order to avoid finding spurious relationships.

The output variable, y15\_Q82, records the answers provided by the sample of workers to the following question: "In the past 12 months, how many days absent from work due to sick leave or health-related leave?" It is a numerical variable, ranging from 0 to 360, and synthetizes the duration of the sick leave taken by the workers.

**Table 3.** Distribution of the sample of workers according to their household's total monthly income.

#### *2.2. Multivariate Linear Regression*

Let us consider a set of *k* + 1 quantitative variables with *y* as the dependent variable and *x*1, *x*2, ... , *xk* as independent variables. The multivariate linear regression method consists of creating a lineal model that predicts *y*, using variables *x*1, *x*2, ... , *xk*. It can be expressed as follows [21]:

$$y = \beta\_0 + \beta\_1 \mathbf{x}\_1 + \dots + \beta\_k \mathbf{x}\_k \tag{1}$$

The parameters' estimation is performed by means of the ordinary least squares approach [22] by means of the following:

$$\min\_{\beta \in \mathbb{R}^{k+1}} \|\|y - X\beta\|\|^2 = \min\_{\beta \in \mathbb{R}^{k+1}} \sum\_{i=1}^{n} \left(y\_i - \beta\_0 - \sum\_{j=1}^{k} \beta\_j \mathbf{x}\_{ij}\right)^2 \tag{2}$$

where · denotes the Frobenius norm.

#### *2.3. Support Vector Machine for Regression*

Let us consider again a set of *k* + 1 quantitative variables with *y* as the dependent variable and *x*1, *x*2, ... , *xk* as independent variables, where each *i* element constitutes a row vector. Let Φ : χ → F be the function that corresponds to each row vector with a point of the characteristics space F . Let us define a function as follows [23]:

$$f(\mathbf{x}) = \prec \omega, \Phi(\mathbf{x}) \succ + b \tag{3}$$

The problem to solve is as follows:

$$\min \frac{1}{2} \|\|\boldsymbol{\omega}\|\|^2 + \mathcal{C} \sum\_{i=1}^{n} \left(\mathbb{\boldsymbol{\xi}}\_i + \mathbb{\boldsymbol{\xi}}\_i^\*\right) \tag{4}$$

Constrains:

$$\begin{array}{ll} y\_i - \varsigma \omega\_\prime \phi(\mathbf{x}\_i) > \underline{\varsigma} + \underline{\zeta}\_i & i = 1, \dots, n \\ \underline{\varsigma} < \omega\_\prime \phi(\mathbf{x}\_i) > -y\_i > \underline{\varsigma} + \underline{\zeta}\_i^\* & i = 1, \dots, n \\ \underline{\varsigma}\_i \ge 0, \ \underline{\varsigma}\_i^\* \ge 0 & i = 1, \dots, n \end{array} \tag{5}$$

The problem complexity depends on the dimension of the row vectors [24]. The solution of this problem gives as a result the model of support-vector machines for regression.

#### *2.4. Genetic Algorithms*

The process of learning by trial and error can be considered as being similar to the natural evolution process. The development of genetic algorithms (GA) started with the works of Holland [25]. GA is a kind of evolutionary algorithm that is based on the evolution of a certain set of solutions trying to either maximize of minimize the result of an objective function. GA is a bioinspired methodology that mimics the procedure of natural selection. The interest in GA methodology in optimization is because they are a global and robust method for finding solutions that do not require any a priori knowledge about the problem.

GA make use of the following three basic operators [16]:

The crossover operator takes two different individuals of the population and creates a new one, mixing the two. The mutation operator performs random changes in those individuals, created with the help of the crossover operator. Mutation makes it possible to introduce new strings in the next generation, giving the ability to search beyond the scope of the initial population. Another interesting mechanism is elitism, which makes a certain number of individuals with a good performance according to the result of the objective function survive and pass to the next generation, without any change.

#### *2.5. Multivariate Adaptive Regression Splines*

MARS is a well-known parametric methodology that builds a non-linear model based on hinge functions. It is expressed by the following equation [26]:

$$
\hat{\mathbf{y}}\_{\rangle} = \beta \mathbf{o} + \sum\_{i=1}^{k} \beta\_{i} \cdot \mathcal{B}(\mathbf{x}\_{i}),
\tag{6}
$$

where *y*ˆ*j* represents each one of the outputs forecasted values per each *yj*, β*i* are the model parameters and *B* are the model basis functions. The basis functions are defined as follows:

$$\begin{aligned} B^- &= \begin{cases} &(t-x)^q & \text{if } x < t\\ 0 & & \text{otherwise} \end{cases} \\ B^+ &= \begin{cases} &(t-x)^q & \text{if } x \ge t \\ 0 & \text{otherwise} \end{cases} \end{aligned} \tag{7}$$

where *q* is a natural number that represents the power function.

When a MARS model is created, there are three different well-known methods that take part in the model in order to assess the importance of the variables. They are the following:

The GCV expression is as follows [23]:

$$GCV(M) = \frac{\frac{1}{n} \cdot \sum\_{i=1}^{n} \left(y\_j - \mathcal{Y}\_j\right)}{\left(\frac{1 - \mathcal{C}(M)}{n}\right)^2},\tag{8}$$

where *C*(*M*) is the complexity penalty function that increases with the number of basis functions in the model and which is defined with this equation, where *M* is the number of basis functions. In the case of the present research, the maximum interaction degree allowed was nine.

The equation for RSS is as follows [27]:

$$RSS = GCV(M) \cdot \frac{\text{3-M}}{n^3} \tag{9}$$

#### *2.6. The Proposed Algorithm*

The proposed algorithm works as is explained here. First, it is initialized with a random population. Each member of the random population represents a subset of all the available variables that will be employed for the forecast of the number of days off for each worker. It is a string, as in the following example: 1100011 ... 0101 with a total of 59 digits, one per variable, where 1 means that the variable is present and 0 that is missing.

In order to know how each of the variables subsets performs, they are employed for training a MARS model, using 80% of the available individuals, while the other 20% are employed for the model validation. This process is repeated 1000 times for each of the variables subset and the average R2 value obtained is used as the result of the objective function. Following the usual methodology of genetic algorithms, the best individuals of the population are selected and crossed.

In the present research, a mutation rate of 10% was allowed, and a 5% of elitism, which means that the 5% of the best individuals of a generation are included in the next one. In the case of the present research, a fine-tuning was performed, testing mutation rates from 0.5% to 15% in steps of 0.5%. The R2 values obtained did not find statistically significant differences from 0.5% to 10%, while in higher mutation values, the R2 decreased. Therefore, results with 10% probability mutation rates are presented. The crossover methodology employed is known as *single point crossover*, in which both parental chromosomes are split in only one point randomly selected. Each generation of the genetic algorithm population is formed by 1000 individuals; this means that there are 1000 different variables subsets. The results shown were obtained after 100 iterations, which means that 100,000 variables subsets were examined. This is a small number if compared with the more than 5.7 × 1017 possible variables subsets that can be obtained for a problem like this with 59 independent variables. Finally, it can also be highlighted that the performance of the algorithm would be improved if those workers whose leave durations are over certain threshold value (i.e., 10 or 100 days) were considered as outliers and removed. The flowchart for the proposed algorithm is shown in Figure 2.

**Figure 2.** Flowchart of the proposed algorithm.

#### **3. Results**

The average value of the R2 for the 1000 models trained for the variables that were finally selected in each case was 74.26%. The RMSE average value was 27.51. The average difference in days in absolute value of forecasted and real number of days off was 10.22. If workers are divided in those with leaves of 10 or less days, 85% of the total, and those with leaves of more than 10 days, the RMSE values obtained are quite different. In the case of those with leaves of 10 or less days, the RMSE value was of 5.56 and 71.49 for those with leaves of more than 10 days. In the case of average difference in days in absolute value of forecasted and real number of days off, it was 4.28 for those with leaves of 10 or less days and 48.81 for those whose leaves are over 10 days. It means that prediction for those leaves of 10 or less days are much more accurate. Figure 3 shows the histogram of the number of leave days for all the workers. It can be observed that most of the leaves are of 10 or less days. The results described in this paragraph are in line with what is said in Section 2.3 about how an outlier's removal would improve the algorithm performance.

In order to assess the performance of the proposed algorithm, it was compared with two benchmark methodologies: linear regression [28] and support vector machines for regression [29]. In both algorithms, 1000 models with different training and validation datasets were tested. For the linear regression, the average R2 value was of 26.72% with an RMSE of 74.21, while in the case of support vector machine for regression the average R2 value was of 67.32% with a RMSE of 29.01. Table 4 shows a comparison of the performance of the proposed algorithm with the two benchmark methodologies referred before, linear regression and support vector machines.

Figure 4 shows the results of one of the models created with these variables. In such a model, the forecasted and real number of days off for all the workers randomly included in the validation dataset can be observed. The values in the horizontal axis are the workers' identificators. Please note that the largest differences of forecasted and real values can be found for workers with numbers from 81 to 83 that would be considered as spurious. In this case, the difference in days of forecasted and

real number of days off was on average −4.255, with a median of 2.53 and a standard deviation of 27.34. In absolute values, the mean was 10.642 days. As after 100 iterations all the models give similar results in terms of R2, Table 5 shows the list of variables selected by this model as the most relevant when predicting absenteeism among workers in the energy sector, using the three importance criteria (nsubsets, gcv and rss) referred to in Section 2.2.

**Figure 3.** Histogram of the number of leave days of all those individuals that took part in the present research.

**Table 5.** A list of variables selected by the model as the most relevant to predict absenteeism among workers in the energy sector.

**Table 5.** *Cont*.

1 nsubsets: number of model subsets that make use of the variable (see Section 2.2). 2 gcv: generalized cross-validation of the variables (see Section 2.2). 3 rss: residual sum of squares (see Section 2.2).

**Figure 4.** Results of one of the models created: forecasted and real number of days off for all the workers randomly included in the validation dataset.

#### **4. Discussion**

The results obtained show that it is possible to make use of machine-learning methodologies such as MARS and GA in order to predict the number of days of health-related leave taken by workers in the energy sector, taking into account a certain number of variables linked to personal and work-related factors for each individual. The results obtained are not surprising, as GA and MARS have proven to be valid in similar scenarios [23] and in other problems linked to the energy field [30,31]. The worst forecasts were obtained for the longest leaves, as they were few and present a behavior outside the normal range. In future research, the MARS model could be substituted by other regression models such as neural networks or support vector machines for regression, and even by replacing GA with other evolutive methods like particle swarm optimization or differential evolution.

As other studies have shown [13,32], there are several factors that have an influence on absenteeism in the workplace. It should be pointed out that, in this case study, sick leave due to common illness and sick leave due to occupational illness or occupational accidents were not analyzed separately. This is the reason behind the multiple kinds of factors that become part of the model, and implies that some difficulties may appear during the discussion of the results [33]; in any case, it is necessary to consider both, in order to understand the causes behind absenteeism.

The model was built with twenty items that can be classified into four categories:

It must be pointed out that several of these items could belong to different categories simultaneously. Therefore, they are included in those ones that better explain their impact as a cause of absenteeism. Moreover, there are many interrelated circumstances between each category, to the extent that it is difficult to find studies that only deal with one of them. In fact, there are several research studies that cover issues in relation to ergonomics while speaking about working conditions, workers' personal characteristics and organizational contexts at the same time [32].

In absolute terms, working conditions have the greatest impact on sick leave among workers, since nine out of twenty items of the model developed fall into this category. However, as mentioned before, several of these items also affect other categories, such as ergonomics and psychosocial factors.

A discussion on the relationship between absenteeism and the items in each category is presented next.

#### *4.1. Ergonomics*

There are five items in the model that show the impact of ergonomics on sick leave among workers. This is in consonance with other researchers' conclusions [33–36] that maintain that poor ergonomics in the workplace and prevalence of musculoskeletal disorders (MSD) are linked and could therefore mean an increase in sick leave. These five items are as follows:

It is remarkable that the model seems to suggest a contradictory idea in relation to repetitive tasking, in that it considers that short repetitive tasks of less than ten minutes have a negative impact on sick leave among workers, whereas short repetitive tasks of less than one minute can reduce absenteeism. Concerning ergonomics, the shorter the task, the more damaging it could be, so at first glance, this would seem to be an error.

On the other hand, there could be several explanations for this curious result. For instance, it is difficult to find a job that requires doing the same repetitive tasks lasting less than one minute for the entire working day. However, it is more feasible to find jobs that include the same repetitive task lasting less than ten minutes for the entire working day. Thus, multitasking would preferably be linked to the first of these cases, and multitasking is a valued characteristic of good ergonomics. This would be a proper explanation for the behavior of these items in the model. In any case, this sets a starting point for future research.

Other items that are included in the model and classified into the ergonomics category, such as monotonous tasks, painful positions and sitting, are ergonomic factors traditionally considered during risks assessments due to their negative impact on MSD prevalence. As a first conclusion, the model has proven that the beliefs as to how ergonomic investments have a positive impact on MSD prevalence and occupational absenteeism are a step in the right direction. This study therefore supports other researchers' conclusions that are applicable to industrial environments [37] and to the energy sector in particular [38].

#### *4.2. Psychosocial Factors*

Only one item falls into this category: that of being in situations that are emotionally disturbing. However, there are several studies that point out that psychosocially demanding working environments have a deeply negative impact on absenteeism [39,40]. Nevertheless, it should be noted that this does not mean that psychosocial factors are less relevant than others. In fact, there is a close interrelation between them and the rest. For instance, as other studies have shown [41], ergonomic interventions could be counterproductive, unless they attend to psychosocial factors.

#### *4.3. Working Conditions*

Working conditions is the category that includes the largest number of items from the model. However, as previously stated, that does not necessarily mean it is the category with the strongest influence on absenteeism. The nine items from the model classified in this category are as follows:

First of all, it must be said that "being satisfied with working conditions" could be included in the category of psychosocial factors; after all, this item depends both on how working conditions are designed and how workers perceive them. In any case, this aspect has already been discussed in the previous section. In fact, the item included on working conditions highlights the interaction between categories and the importance of psychosocial factors on the control of absenteeism.

Working conditions cover a great spectrum of factors that can become causes of occupational absenteeism. Indeed, workers' general health status, MSD prevalence and other diseases are closely related to the working environment. Therefore, in terms of the energy sector in particular, each company has to understand and act on several working conditions, to develop health management, as other works have already shown [42].

A conclusion that can be obtained by analyzing the factors that appear in this category is, as has been the case with others, the existence of a link with the category concerning personal conditions. For example, it has been proven that a low income of the worker increases the probability of his/her taking sick leave. One wonders if, in fact, this circumstance reveals the effects of socioeconomic status on occupational absenteeism [43].

Other working condition factors included in the model show the effects of workplace environment on sick leave: noise, air pollution, extreme temperatures, etc. This was to be expected, because it points in the direction of occupational diseases, in line with several previous studies [44,45].

The only remarkable item that could be seen as a contradiction is that workers with more than one job seem to be less vulnerable to sick leave. Several explanations can be suggested. For instance, having more than one job could be linked with multitasking and, therefore, lower prevalence of MSD. Another possible reason could be that this kind of worker is more likely to belong to a precarious social stratum where health damages are sometimes underreported. In any case, this could be the subject of a future line of research.

#### *4.4. Personal Data and Physiological Characteristics*

There are five possible causes of absenteeism included in this category:

This category joins together several factors related to lifestyle that are too difficult to identify and analyze properly, especially when the studied variable includes both common and occupational diseases as a cause of absenteeism. In fact, there are many studies that deal with this subject without achieving unanimity [46–48].

It seems to be expected that general health status [49] and, as a result, other factors that can alter it, like age, must affect the prevalence of several illnesses that end up causing sick leave. There are other studies that go further and that try to analyze gender differences in this matter [14]. In this category, however, everything is vaguely interrelated, so there are many questions to answer in future research.

Apart from the item that refers to activities outside work, every factor included in this category is a cause or an effect of the workers' general health status. This appraisal has many implications that must be kept in mind when planning any move designed to reduce absenteeism among a company's workforce.

#### **5. Conclusions**

It is possible to make use of machine-learning methodologies, such as MARS and GA, in order to create models able to predict the number of days of health-related leave among workers in the energy sector. Absenteeism can be monitored and predicted by using a model that employs several items included in the following categories:

These categories are all interrelated, and it is difficult to establish boundaries between them, but as a positive consequence of this, acting on one of them to reduce absenteeism in a company could have a great impact on the others.

Any management program has to act on factors that affect the employees' general health status, process design, workplace environment, ergonomics and psychosocial working context, among others, if it is going to be successful. This has an extensive field of application in the energy sector, where most of the activities are undertaken in an industrialized context.

**Author Contributions:** Conceptualization, F.S.L.; formal analysis, F.S.L.; investigation, A.G.F., N.M.B.S. and A.S.S.; methodology, F.S.L. and A.S.S.; project administration, A.S.S.; resources, A.G.F.; validation, G.F.V.; writing—original draft, A.G.F. and N.M.B.S.; writing—review and editing, G.F.V. and A.S.S. All authors have read and agree to the published version of the manuscript.

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

**Acknowledgments:** The authors wish to acknowledge the European Foundation for the Improvement of Living and Working Conditions (Eurofound), as well as the UK Data Service, for providing us with the results of the database used in this work.

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

#### **References**

© 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* **Energy Multiphase Model for Biocoal Conversion Systems by Means of a Nodal Network**

#### **Beatriz M. Paredes-Sánchez 1, José P. Paredes-Sánchez 1,\* and Paulino J. García-Nieto 2**

Received: 28 April 2020; Accepted: 25 May 2020; Published: 29 May 2020

**Abstract:** The coal-producing territories in the world are facing the production of renewable energy in their thermal systems. The production of biocoal has emerged as one of the most promising thermo-energetic conversion technologies, intended as an alternative fuel to coal. The aim of this research is to assess how the model of biomass to biocoal conversion in mining areas is applied for thermal systems engineering. The Central Asturian Coal Basin (CACB; Spain) is the study area. The methodology used allows for the analysis of the resource as well as the thermo-energetic conversion and the management of the bioenergy throughout the different phases in a process of analytical hierarchy. This is carried out using a multiphase mathematical algorithm based on the availability of resources, the thermo-energetic conversion, and the energy management in the area of study. Based on the working conditions, this research highlights the potential of forest biomass as a raw material for biocoal production as well as for electrical and thermal purposes. The selected node operates through the bioenergy-match mode, which has yielded outputs of 23 MWe and 172 MWth, respectively.

**Keywords:** biomass; bioenergy; energy production system

#### **1. Introduction**

The constant technological progress and the increasing industrialization of society have boosted the demand for energy. Fossil fuel deposits, such as those in coal basins, are limited and have been extensively exploited, so modelling the biomass potential of renewable sources in these areas is a challenge for the future in the European Union (EU). Spain boasts a wide variety of renewable resources such as biomass and wind and solar energy, all available as renewable energy sources [1]. The technological capacity of Spanish industry has made it a benchmark in the use of renewable resources [2]. The potential of Spain in renewable energies is well above both the domestic energy demand and the existing fossil fuel resources [3]. Despite this situation, Spain is highly dependent on foreign energy from fossil fuels, coal being the main source of indigenous energy in Spain [4].

Renewable energy sources are found in nature, have the capacity to be totally or partially regenerated, and can be used for energy purposes. Biomass is part of a continual cycle of mass and energy consumption and production in the environment. It can be used to produce energy, either directly by combustion or indirectly by way of biofuels. EU Directive 2003/30/EC [5] defines biomass as a biodegradable fraction of products, waste and residues from agriculture, including vegetal and animal substances, forestry and related industries, as well as the biodegradable fraction of industrial and municipal waste. Such a definition has a comprehensive character since it includes a variety of energy sources sharing certain characteristics but differing in their origin and the technology used to obtain and use them, one of the main sources being forest biomass.

The EU Forestry Action Plan [6] calls for an assessment of forest biomass availability for energy production in both national and regional energy systems. The most promising scenarios for future use are

This study is to be carried out in the Region of Asturias in northern Spain, a region rich in fossil fuels, i.e., coal, where the production of alternative biofuels has been of interest in the context of the EU [7,8]. Its producing industry uses coal as the main source of indigenous primary energy, a part of which is produced in its mining basins [9]. Mining has gradually declined in recent years following the reduction of its activity [10]. Biomass from forests and forest management in the mining areas is conditioned upon its frequency of occurrence along with its topography. It is, therefore, necessary to carefully plan management procedures based on the distribution of potentials throughout the whole territory. This strategic framework calls for an analysis of the most appropriate techniques for extraction, storage, transport and use of biomass in forests. However, there is limited use of forest biomass to provide an alternative to fossil fuels owing to its low energy density [11].

Biocoal is achieved through a roasting process (i.e., torrefaction process), a high potential process to be applied in the production of solid biofuels from lignocellulosic biomass [12]. In this way, the quality of this fuel increases by virtue of its hydrophobicity and higher calorific value compared to the initial biomass [13]. Biofuels considerably improve the potential of biomass for industry and thermal systems [14]. Furthermore, it is of interest for both industrial energy transition and industry 4.0 [15]. However, modelling, conversion and energy management of resource supply are some of the main challenges in the current energy context [16]. In this regard, Visa et al. [17] have defined the importance of characterizing energy production systems according to the nature of the energy resource used. Paredes-Sánchez and Ochoa-Lopez [18] established the importance of biomass modelling as an alternative to coal in the energy production systems. Therefore, a proper implementation in a traditionally coal-based industry, as is the case of the Principality of Asturias, requires the development of mathematical models for the energy conversion of its resources.

Extensive research has been carried out using both resource characterization and energy conversion models from different perspectives [19–22]. These current works discuss numerous challenging issues, including the increasing number of assumptions that ensure consistency between the models used and the large amount of data required for energy conversion. The combination of modelling methods allows for a significant number of variables to be taken into account to solve the problems of energy conversion [22]. The multiphase model is an important step in developing the knowledge needed to improve energy fuels in thermal systems. In this context, research and development activities on biomass torrefaction have been particularly active in exploring its potential as a fuel [23,24]. However, the modelling and management of the torrefaction process can be found in the literature as a current challenge [25,26]. Kumar et al. [27] have studied the processes of bioenergy transformation, pointing out the necessity to search for comprehensive analytical models to overcome the existing limitations and provide the necessary technology to enable industrial use of high value-added biofuels. Huntington et al. [28] stress the importance of developing advanced mathematical models focused on the potential supply of biomass resources as a source of bioenergy. Paredes-Sánchez et al. [29] point out the limitations in modelling the potential use of biofuels for energy production systems considering the scope of the energy conversion system. In this regard, Bach et al. [30] define the demand for a comprehensive model of biomass torrefaction, which can provide interdisciplinary information to industrialize and commercialize the process. In this framework, a comprehensive analysis of the energy use of biocoal, starting with the supply of biomass as raw material and ending with its final energy conversion through fossil-fuel-based technologies for thermal systems, is a difficult undertaking in the field of bioenergy. Therefore, the present paper aims to develop a mathematical model that would

allow a comprehensive analysis of the potential use of available biomass in the mining areas to produce biocoal for thermal systems engineering. The main breakthrough consists of applying a multiphase model to the data to characterize thermo-energetic conversion and energy management of biomass as a raw material for the production of biocoal. These studies are carried out in the Carboniferous Basin of the Principality of Asturias (CACB).

This work is organized as follows: Section 1 introduces the context of the research and the aim of the work. Section 2 shows the methodology used for the study area and the modelling process to produce biocoal. Section 3 shows both its findings and the details of its implementation as results. Section 4 provides a more detailed analysis of and discussion on biocoal as an alternative fuel to coal in mining areas for energy production systems. Section 5 provides the main conclusions of the study in the mining area by means of energy production systems.

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

#### *2.1. Study Area and Mine Nodes*

The National Renewable Energy Action Plan (NREAP) in Spain is aimed at complying with European Directives 2009/28/EC and 2009/29 EC on the contribution of renewable energy and the reduction of greenhouse gas (GHG) emissions by 2020 [3]. Located in the north of Spain, the Principality of Asturias is an Autonomous Community where 45% of the land is forested. The energy structure of the Principality of Asturias is conditioned by the contribution of fossil fuels to the national energy system as a whole. The Central Asturian Coal Basin (CACB) spreads over the southern councils of the central area of the region (Figure 1).

**Figure 1.** Location of the study area—Principality of Asturias (Spain).

The CACB features the geological resources that represent one of the main coal mining deposits in Spain, located above the geological unit of study [31]. The use of coal as an integral source of indigenous energy is the basis of the Asturian economy. Therefore, the CACB constitutes a geographical, economic, and geological unit that is the study area. In this sense, mine shafts are defined as candidate nodes of analysis to implement the use of biocoal in energy production systems. It should be noted that the layout and infrastructure of the coal mine itself, i.e., the electricity grid or the transport network, would favor the use of biocoal as an alternative fuel. Table 1 shows the mine shafts considered as modelling nodes.

**Table 1.** Nodes defined in the Central Asturian Coal Basin (CACB) study area.

The analysis of the biomass in the forests surrounding the mine shafts allows unused resources to be valorized for energy production. The use of forest biomass for energy purposes will not only increase the economic development, but also energy self-sufficiency [32,33]. The use of forest residues will also generate environmental benefits [34].

#### *2.2. Multiphase Mathematical Model*

The multiphase mathematical model developed for this work comprises three phases in the area of study: resources, thermo-energetic conversion and energy management, each one of which is described in detail in Figure 2.

**Figure 2.** Flowchart of the multiphase mathematical algorithm through mass and energy balance.

#### 2.2.1. Phase 1: Resources

For the assessment of potential resources, the Geographic Information System (GIS) database of the BIORAISE GIS tool from the Research Centre for Energy, Environment and Technology (CIEMAT) was used in order to collect data on forest biomass in the CACB study area [35]. In Phase 1 of the mathematical model, it is estimated that the existing resources for each candidate node come from residues of cleaning activities and forest management operations in the study area and are gathered in each of the considered nodes (Figure 2). Mass and energy are expressed in terms of dry ton (dry t). In this study, this biomass is considered as the raw material to be used to generate biocoal through the torrefaction technology in Phase 2 of the model with the roasting system.

The potential mass (M) is the residual forestry biomass, such as branches and leaves, found within the study area and therefore, in the surroundings of each modelling node. The available mass (m) consists of branches and tops (including leaves). It is obtained from cleaning, thinning and felling operations from BIORASE GIS database (Figure 2), which takes into account the techno-economic constraints of the potential biomass to define the useful resources [35]. Such restrictions derive from harvesting procedures, which depend on the terrain conditions to access the raw material [7] that will define the useful resources, including available mass of conifers, wood and mixtures (dry t/year). Techno-economic constraints affect its use from the point of view of mass and energy balance. The biomass is evaluated in the surroundings of the selected node in the study area using Equation (1).

$$\mathbf{m} = \sum\_{i=1}^{n} \mathbf{m}\_i \tag{1}$$

where

m is total available mass (dry t/year); and

mi is total available mass of conifers, hardwood and mixtures (dry t/year).

The energy of the available residues (E) is the result of Equation (2), where the Lower Heating Value (LHV) is used.

$$\mathbf{E} = \sum\_{i=1}^{n} (\mathbf{m}\_i \cdot \mathbf{LHV}\_i) \tag{2}$$

where

E is energy from available mass (GJ/year);

mi is total available mass of conifers, hardwood and mixtures (dry t/year); and

LHVi is Lower Heating Value of conifers, hardwood and mixtures (GJ/dry t).

Phase 2 of the calculation algorithm is defined from the available resource and energy (Figure 2).

#### 2.2.2. Phase 2: Thermo-Energetic Conversion

In Phase 2, the thermo-energetic calculation algorithm used to calculate the mass and energy potential of biocoal to be produced is based on Equations (1) and (2), respectively. Torrefaction is a thermal pre-treatment carried out by a reactor system at atmospheric pressure, with no oxygen. The roasting process takes place at temperatures between 200 ◦C and 300 ◦C to achieve more uniform solid biofuel, whose two main thermo-energetic characteristics are energy yield and mass yield. The energy and mass yields [36] are defined from the reactive part of the biomass, which is turned into biocoal through torrefaction as energy and mass percentages or fractions [36,37]; therefore, both ash and free water content are excluded from the definition. The mass yield (αy) is the correlation between the mass of the biocoal produced (mb) and the mass of the raw material, i.e., total available mass (m), obtained from the resources in the roasting system per candidate node (Equation (3)).

$$
\alpha\_{\rm y} = \left(\frac{\rm m\_{\rm b}}{\rm m}\right) \tag{3}
$$

Energy yield (βy) is defined in the model according to Equation (4)

$$
\beta\_\text{y} = \alpha\_\text{y} \cdot \left(\frac{\text{LHV}\_\text{b}}{\text{LHV}}\right) \tag{4}
$$

where

"βy" is the yield referred to the LHV parameter already considered in Phase 1;

"αy" is mass yield; and

LHVb and LHV are lower heating value of biocoal mass and raw material, respectively [36].

The available energy per unit of mass corresponds to LHV because it represents the energy that can be efficiently recovered after combustion. To this objetive, the torrefaction process is assumed to be carried out on the available mass as a whole within the areas surrounding each node, thus obtaining Equation (5):

$$\mathbf{m\_b} = (\sum\_{i=1}^{n} \mathbf{m\_i}) \cdot \alpha\_{\mathbf{y}} \tag{5}$$

where

mb is biocoal mass from available mass per node (dry t/year);

mi is mass of available biomass (dry t/year); and

αy is mass yield.

For the calculation of energy potential to be obtained as primary energy from biocoal, Equation (6) is applied as an approach.

$$\mathbf{E\_b = \sum\_{i=1}^{n} (m\_i \cdot \text{LHV}\_i) \cdot \beta\_\mathbf{y}} \tag{6}$$

where

Eb is bioenergy as biocoal per node (GJ/year);

mi is available mass (dry t/year);

LHVi is lower heating value of conifers, hardwood and mixtures (GJ/dry t); and βy is energy yield.

Equations (3) and (4) show the technical feasibility of placing a roasting system at a candidate node in the study area by bioenergy-match mode. Such operating conditions correspond to short residence periods, less than 30 min, with temperatures above 260 ◦C. Favorable conditions are considered to be those with energy and mass efficiency above 95% and 90% respectively. If the initial biomass to be torrefacted is dry, with moisture content below 10–15%, lower energy efficiency can be expected [37].

#### 2.2.3. Phase 3: Energy Management

Biomass represents the main manageable renewable energy and should therefore play a role in energy production to replace or complement fossil fuels (Phase 3). Introducing biocoal in the industrial sector, and especially in thermal demand, calls for a change in the energy supply management model [38]. Energy production systems consist of sets of technologies that transform raw energy from one fuel into final use energy [17], allowing for the production of heat and/or electricity and of new biofuels. Based on the analysis of energy management, the energy to be produced will be assessed on the basis of the biocoal production in the roasting system as compared to its energy-based equivalence to the mass (mf) of other conventional fuels, Equation (7).

$$\mathbf{m\_{f}} = \mathbf{E\_{b}} \ / \mathbf{C} \mathbf{V\_{i}}\tag{7}$$

This is done by taking the calorific value (CVi) of the fuels in Table 2 [39] as an approximation.

**Table 2.** Energy-based equivalence of certain fuels.

The model of analysis using the multiphase mathematical algorithm considers a maximum distance of 50 km for the resources around each candidate node, which is a limit distance determined for their transformation into conventional biofuels that can be used in that territory [40]. Consequently, the calculation area is restricted to a maximum of about 50 km away from each candidate node for the entire analysis (Phases 1, 2, and 3). Once the final combination of the results per candidate node has been determined, the optimal node, i.e., the location of the roasting system in one mine shaft, will be defined for the harvesting of the available resources in the study area.

#### **3. Results**

In Phase 1, the analysis model established by the multiphase mathematical algorithm (Figure 2) shows the quantities of both potential biomass and available biomass in the study area for each of the nodes considered, as shown in Figure 3.

**Figure 3.** Mass per candidate node in the study area: (**a**) potential biomass (M) and (**b**) available biomass (m).

Figure 3 shows some linearity in the distribution of the potential biomass around each modelling node considered in the CACB. It indicates a regular distribution of the biomass with the increase of the distance from 10 km to 50 km. Increasing the distance yields more available biomass around each candidate node (Equation (1)). The available energy based on the available biomass per candidate node is shown in Figure 4.

For each node, the amounts of the available biomass that can be collected are above 75 dry kt/year, around 1300 TJ/year. The available biomass will be used as raw material to be converted into biocoal (Figure 4).

**Figure 4.** The energy (E) of the available biomass per candidate node as calculated by Equation (2).

Phase 2 takes into account the thermo-energetic conversion of the available biomass into biocoal as per Equations (5) and (6). Figure 5 shows the results of the biocoal potential and its equivalent energy per modelling node with the increase of the distance from 10 km to 50 km. The thermo-energetic conversion does not affect the distribution because it has a similar shift of the mass and energy per candidate node in the study area.

**Figure 5.** Mass and energy per candidate node in the study area: (**a**) biocoal (mb), (**b**) bioenergy (Eb).

Furthermore, the LHV distribution of the biocoal obtained from the raw material collected in the study area has been evaluated per candidate node (Figure 6).

As can be seen, there is a homogeneous distribution of the energy quality of the theoretical biocoal to be obtained, reaching an overall average value of about 18.5 GJ/dry t.

Finally, Phase 3 of the algorithm is implemented in each node, based on the data in Table 2. Equation (7) shows this potential as an alternative fuel to coal for energy conversion. This assessment characterizes biocoal as compared to different types of fossil fuels for energy production systems (Figure 7).

**Figure 6.** Lower heating value (LHVb) distribution of the obtainable biocoal in the study area.

**Figure 7.** Energy-based equivalent fuels from biocoal, as a substitute in energy production systems.

Taking all candidate nodes, the potentially producible biocoal in terms of energy is equivalent of above 28 kt/year of natural gas, 31 kt/year of distilled oil, and 68 kt/year of brown coal.

#### **4. Discussion**

Implementing a technological alternative requires splitting up the analysis by considering the different aspects involved [41]. The paper explores biocoal use as an alternative to coal production in mining basins by using mine shafts as energy conversion nodes that draw on the surrounding forest biomass as raw material. The findings of the study reveal that the proposed model, based on a multiphase algorithm, is highly useful in the study of energy conversion and management for thermal systems. A cumulative analysis of these facts, augmented by other findings of this work, is summarized below.

The potential renewable in the world is diverse and widely distributed in nature [1]. However, most of it may present supply fluctuations, giving rise to risks and uncertainties in serving regional energy markets. In this sense, the use of biocoal in mining environments, as proposed, guarantees the use of renewable resources and their conversion, since it is based on substituting the use of fuel, i.e., coal, with one of similar characteristics, i.e., biocoal. Moreover, there is a whole industrial structure with mature technologies for the use of this coal. Despite this, one of the remaining challenges for the development of biofuels in the industry is whether a sufficient amount of raw material can be generated to meet the growing demand and achieve a shift from fossil fuels to biofuels [42].

In this regard, the applied algorithm shows that any node considered as a thermo-energetic conversion node for biocoal production has at least 75 dry kt/year of raw material for biocoal production in the study area. Here, given the nature of the area under consideration and the nodes, only one node may be selected, given the territorial extension of the CACB.

Biocoal could replace coal in the production of heat and power through energy production systems. It would also have a large impact on industry [43]. Providing the infrastructure, technology and research to enable the use of renewable energy sources available "on site" is a challenge for industry, in particular for readily available sources such as biomass [38]. The continued development of industrial technology is expected that will ensure more competitive, reliable and sustainable energy production systems [44]. The use of biomass as an alternative to coal in certain contexts has been made feasible by using biocoal technology, but advances in technology alone do not promote its widespread use; therefore, extensive studies are needed. In this sense, from the mass and energy balance of the analysis model, the optimal node appears to be PC, this is selected as the optimal node for the CACB. With the use of a roasting system, this node would produce about 84 kt/year of biocoal, equivalent to 1544 TJ/year of biocoal, the largest amount of all the nodes selected in the proposed analysis.

Renewable energy generation systems are illustrative of the fact that some renewable energy infrastructures are at a significant distance from energy conversion systems [41]. Overcoming this barrier to competitiveness of renewable energies requires a well-planned and carefully managed energy supply infrastructure during and after the infrastructure investment. An environment based on fuel-coal energy production systems in the mining basins benefits both their use and direct conversion for electrical, thermal or industrial purposes, given the development of a coal energy conversion industry in the coal mining environment. However, all this requires some analysis of future techno-economic feasibility as well as in-depth studies based on the development of these technologies.

Biomass holds the potential to be considered the best alternative to meet global energy demand sustainably and reduce the impact of polluting energy resources [45]. In this context, the development of biorefineries for the production of advanced biofuels is encouraged. Biorefinery represents a sustainable means of generating multiple bioenergy products from various biomass raw materials via the incorporation of relevant conversion technologies. Biorefinery is crucial in the transition of various traditional industries into a circular bioeconomy in the context of energy transition [46]. That opens up the door to the development of infrastructures such as the biorefinery in the aforementioned node. The use of biomass in energy production systems such as co-firing coal with biomass to generate energy is gradually increasing even though their performance differs significantly due to the wide variations in its physical and chemical properties [47]. However, biocoal overcomes most of its use limitations for heat or electricity production.

In this context, in the shaft named "Pozo Candín" (PC), considering the efficiency of a heat production system with a thermal efficiency of 80% and about 2000 h of operation per year [48], it is possible to install a total power of 172 MWth for the energy conversion of biocoal. Additionally, if it is considered the alternative objective of electricity production, either by co-combustion or by itself, looking at an electricity efficiency in the complete conversion of 40% for about 7500 h, it would be possible to achieve an electricity potential equivalent to 23 MWe [49,50]. Detailed studies of biomass conversion stages within the energy production system are working lines for the future in order to further develop more efficient and advanced energy conversion systems.

#### **5. Conclusions**

Biocoal makes it possible to overcome many of the barriers that condition its use and the development of biomass in mining basins around the world due to the traditional interrelationship between the resource and conversion technology for its viability. This is an opportunity to focus on using it in energy production systems that use coal as an energy source, as it overcomes many of the limitations traditionally associated with biomass. In this respect, a multiphase mathematical

algorithm based on operational resource data, thermo-energetic conversion, and energy management has been developed.

The research in this paper highlights the potential of forest biomass in mining areas around the world, in this case in the CACB in Spain, and identifies possible options for the use of forest residues as raw material for solid biofuels as an alternative to coal in carbon basins, i.e., biocoal. Taking a minimum value considering all the nodes in the study area, the energy potential for biocoal is equivalent to over 28 kt/year of natural gas, 31 kt/year of distilled oil, and 68 kt/year of brown coal.

Overall, the quantities of biomass that can be collected in all the nodes studied are above 75 dry kt/year (about 1300 TJ/year). The optimal node that guarantees the largest amount of energy, from a point of view of the proposed mass and energy balance in the study area, is the "Pozo Candín" (PC), which would allow biocoal to be produced through a roasting system at a rate of 84 kt/year, equivalent to 1544 TJ/year in the considered operation conditions with the biomass of the CACB. The production potential of electrical and thermal energy by thermal systems under the defined conditions in this node amounts to 23 MWe and 172 MWth, respectively.

Implementing this energy potential on an industrial scale in the CACB requires techno-economic and thermal systems studies based on the specific characteristics and objectives of the facilities that will be using it. Future understanding of the range of benefits and challenges when introducing and up-scaling biocoal production under different scenarios will depend on detailed, comprehensive, and simultaneous assessments, technological options, and final techno-economic factors.

**Author Contributions:** Conceptualization and methodology, B.M.P.-S. and J.P.P.-S.; software, J.P.P.-S.; writing—original draft preparation, B.M.P.-S.; writing—review and editing, B.M.P.-S.; visualization, J.P.P.-S. and P.J.G.-N.; supervision, B.M.P.-S., J.P.P.-S., and P.J.G.-N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding. It is an initiative of research through an experimental software implementation for technology transfer in thermal systems in the University of Oviedo.

**Acknowledgments:** The authors would like to acknowledge the support of the technical staff of the Research Centre for Energy, Environment and Technology (CIEMAT), throughout the various stages of this work.

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

#### **References**

© 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* **Non-Intrusive Load Monitoring (NILM) for Energy Disaggregation Using Soft Computing Techniques**

### **Cristina Puente 1,\* , Rafael Palacios 2, Yolanda González-Arechavala 2 and Eugenio Francisco Sánchez-Úbeda 2**

Received: 31 March 2020; Accepted: 4 June 2020; Published: 16 June 2020

**Abstract:** Non-intrusive load monitoring (NILM) has become an important subject of study, since it provides benefits to both consumers and utility companies. The analysis of smart meter signals is useful for identifying consumption patterns and user behaviors, in order to make predictions and optimizations to anticipate the use of electrical appliances at home. However, the problem with this kind of analysis rests in how to isolate individual appliances from an aggregated consumption signal. In this work, we propose an unsupervised disaggregation method based on a controlled dataset obtained using smart meters in a standard household. By using soft computing techniques, the proposed methodology can identify the behavior of each of the devices from aggregated consumption records. In the approach developed in this work, it is possible to detect changes in power levels and to build a box model, consisting of a sequence of rectangles of different heights (power) and widths (time), which is highly adaptable to the real-life working conditions of household appliances. The system was developed and tested using data collected at households in France and the UK (UK-domestic appliance-level electricity (DALE) dataset). The proposed analysis method serves as a basis to be applied to large amounts of data collected by distribution companies with smart meters.

**Keywords:** NILM; disaggregation methods; non-intrusive load monitoring; appliance consumptions; soft computing

#### **1. Introduction**

Two of the main global problems that we are currently facing are pollution and consumption control. In the Paris COP (United Nations Climate Change Conference 2015), the United Nations agreed to limit global warming to 1.5 degrees by 2100 and, therefore, reducing energy consumption has become a key task for achieving this goal [1]. The fact that 27% of electricity consumption in Europe is attributed to households emphasizes the need to enact regulations promoting suitable and responsible electricity usage. In this sense, monitoring home energy consumption is an important task in order to optimize and reduce electricity usage [2,3]. Therefore, there will be great benefits if behavioral patterns on appliance usage could be automatically detected with an eye to modifying consumer habits [4], with a potential reduction of 12%, depending on the type of feedback that is provided, as depicted in Figure 1.

Thus, electricity consumption has been studied progressively more, since it has benefits for both sides: consumers and the energy companies. With regard to consumers, consumption control reduces their demand for energy [5], as feedback in this area is proven to lead to a reduction in billing of 3% to 12% [6,7]. Further, disaggregation can be used to detect broken appliances [8] and to check if appliances were left on, as proposed by Bidgely [9], by using a smartphone application. With recent and upcoming electricity tariffs, which may change dynamically depending on current demand, users may benefit from a smart system that suggests how to delay or advance the running of certain electrical appliances. At this point, the use of smart meters has been proven to be the cheapest and most effective way to control consumption [3].

**Figure 1.** Graphical representation of independent studies analyzed by [3] about residential savings resulting from different types of consumption feedback.

In the case of energy providers, a disaggregated bill can be employed to provide personalized energy saving recommendations [7,10], grid control [11], predictions [12], failure detection [13], and similar statistics.

More than 30% of home consumption is from basic appliances, like the washing machine, refrigerator, and oven. The different behaviors of these appliances make it difficult to detect their patterns based on aggregated consumption. Behavioral patterns may differ in some cases, although, in others, they are relatively close, and could be classified into four general types according to their operational states [14], as shown in Figure 2.

Some appliances are relatively easy to detect, as their behaviors have a characteristic pattern, as shown in Figure 3 [5], for the case of a refrigerator. However, not all appliances have such distinguishable patterns—ranging from kettles to washing machines—which makes their detection a difficult task, as some of their stages behave similarly. At this point, fuzzy clustering plays a role to allow for more than one clustering of classification with different degrees of belonging.

Our proposal identifies consumption patterns by detecting changes in power and building a box model to express consumption as a sequence of elastic bars with different powers and different durations. The system was built and tested using a four-year dataset with standard household consumption data collected via smart meters.