Short-Term Forecasting Photovoltaic Solar Power for Home Energy Management Systems

: Accurate photovoltaic (PV) power forecasting is crucial to achieving massive PV integration in several areas, which is needed to successfully reduce or eliminate carbon dioxide from energy sources. This paper deals with short-term multi-step PV power forecasts used in model-based predictive control for home energy management systems. By employing radial basis function (RBFs) artiﬁcial neural networks (ANN), designed using a multi-objective genetic algorithm (MOGA) with data selected by an approximate convex-hull algorithm, it is shown that excellent forecasting results can be obtained. Two case studies are used: a special house located in the USA, and the other a typical residential house situated in the south of Portugal. In the latter case, one-step-ahead values for unscaled root mean square error (RMSE), mean relative error (MRE), normalized mean average error (NMAE), mean absolute percentage error (MAPE) and R 2 of 0.16, 1.27%, 1.22%, 8% and 0.94 were obtained, respectively. These results compare very favorably with existing alternatives found in the literature.


Introduction
Photovoltaic (PV) power generation has achieved enormous development in recent years, mainly for becoming a significant component of the modern power industry's decarbonization. Considering PV power generation in buildings, the PV panels are strongly evolving from off-grid systems to grid-connected systems associated with smart energy management systems (EMS) solutions. The PV power output depends on the available solar irradiation incident on the panels, and this value is not uniform over time. A part of the fluctuations is deterministic, due to location of the panels and rotational and translational movements of the planet in relation with the sun, being this part correctly described by physical equations [1], as may be found in [2]. However, another part of the fluctuations in solar irradiation availability is due to the presence of clouds, cloud mass, aerosol particle concentration, wind speed and direction, ambient temperature, among others, which stochastically reduce the PV panel power output [3,4]. Moreover, the PV panels' output power depends also on internal factors, as photovoltaic module temperature affects the radiation power conversion efficiency [3].
These unpredictable factors generate uncertainty associated with its forecasting to the power grid and to the implementation of the EMS that is dependent on the forecasting of PV power output. In this sense, PV power output accurate forecasting is a crucial challenge to be overcome to achieve massive PV integration, encouraging the development of many studies at a global level [1].
Inventions 2021, 6, 12 3 of 23 DNI forecasts are especially crucial for the management and operation of CPV power plants. DNI is impacted by phenomena that are very difficult to forecast, such as cirrus clouds, wildfires, dust storms, and episodic air pollution events, reducing DNI by up to 30% on otherwise cloud-free days. For non-concentrating systems (i.e., most PV systems), primarily global irradiance (GI = diffuse + direct) on a tilted surface is required, which is less sensitive to errors than DNI since a reduction in clear-sky DNI usually increases diffuse irradiance. For higher accuracy, forecasts of PV-panel temperature are needed to account for the (weak) dependence of solar-conversion efficiency on PV-panel temperature [6]. This section does not aim to describe the physics behind solar energy PV power generation, as they are well known and widely available sources in the literature. The interested reader can find details in [7][8][9].
Forecasting models for solar irradiance can be broadly classified as [10]: statistical, cloud-based, and numerical weather prediction (NWP) models.
Statistical models can be further subdivided into linear models (persistent forecasts, using the clearness or clear sky indexes, or Fourier series expansions, AutoRegressive-Moving-Average-ARMA, AutoRegressive Integrated Moving Average-ARIMA or Classification And Regression Trees-CART models) and non-linear, typically computational intelligence-based models (neural networks, wavelet, fuzzy or evolutionary). Cloud-based models use weather satellite images or ground-based sky (GBS) images [11] to improve solar irradiance forecast. Basically, by processing previous sky images, clouds can be detected, and their motion extrapolated using motion vector fields, or sky cover indices can be obtained, and their evolution forecasted. GBS images have been shown to improve statistical models' performance for forecasts up to 6 h in advance. NWP models are used to forecast the state of the atmosphere up to 15 days ahead. The main limitation of NWP forecasting is its coarse resolution. Spatial resolution can be as low as 1 Km, but the minimum temporal resolution is typically 3 h, which does not allow its use for any patterns less than this value.
The choice of solar-forecasting method depends strongly on the timescales involved, which can vary from horizons of a few seconds or minutes (intra-hour or very short forecasts, for control and adjustment actions), a few hours (intra-day or short/medium, for energy resource planning and scheduling as well as for the electricity market), to a few days ahead (intra-week or long, for unit commitment and maintenance schedules) [12]. The ability to forecast solar irradiance in near-real-time, or nowcasting, is crucial for network operators to guarantee power grid stability with specific reference to power plant operations, grid balancing, real-time unit dispatching, automatic generation control and trading, and for home energy management systems. As this work aims to incorporate the PV power forecasts in HEMS, the authors are interested in real-time, intra-hour and intra-day forecasts.
The accuracy of solar radiation forecasting has direct financial impacts, at various levels. For example, in [13] the effect of DNI forecast accuracy on concentration solar thermal (CST) plants' financial value, incorporating energy storage is examined. When the RMSE of a 48-h DNI forecast is between 325 and 400 W/m 2 , a 1 W/m 2 improvement increases the financial value by $400-1300 per six months' operation for a CST plant. Excellent reviews of solar irradiance forecasting are available, such as [10,[14][15][16][17][18].
Focusing now on PV power generation forecasting, it is, in a simplistic definition, the application of techniques that use registered historical data to make informed prediction/estimation of the future value that PV power will assume in a defined prediction horizon. These forecasts can be obtained directly, i.e., a model is developed using inputs lags of the endogenous variable and, possibly, exogenous ones. Alternatively, forecasting models of those exogenous variables are designed first, and their outputs are fed to a static model which generates the forecasts of the PV power. In [19], the correlations between PV power output and several input variables are explained and computed, concluding that the solar irradiance's correlation achieves large values such as 0.988, followed by the panel temperature and the ambient temperature. The same categories of models pointed out for solar irradiance forecasting apply to PV generation forecasting [12,[20][21][22][23].
Throughout the years, the use of computational intelligence models for forecasting applications has been steadily increasing and, this way, are the focus of this synthetic survey presented in the next sub-section. Computational learning is a mathematical and theoretical field of artificial intelligence that analyses machine learning algorithms' design and their learnability capacity while improving the necessary functions' accuracy and efficiency. These models are based on estimating a function deduced only from samples of training data describing a specific system's behavior and are well suited when physical features are not known [20]. As these are data-based models, the lack of accurate data can become an issue for computational learning methods [20], as the accuracy is strongly dependent on the quality and amount of the available data.
The most used computational learning architectures used for PV power forecasting are shallow and deep artificial neural networks, support vector machine and fuzzy models.
For solar power forecasting, ANNs are among the most employed machine learning techniques. Artificial neural networks are inspired in the biologic neuron and networks. A neural network's fundamental processing element is a neuron, which performs a simple computation on its inputs. Neurons are typically grouped into layers. The network usually consists of an input layer, some hidden layers, and an output layer. ANNs can somehow be classified by the neuron model employed, number of layers, neurons interconnection and the training mechanism. The forecasting performance of an ANN model depends not only on the above items, but also on the design data, and the inputs employed.
The ANN model used in this work is a radial basis function (RBF) Neural Network (NN). The hidden neurons employ a radial type of function, typically a Gaussian, their outputs being linear combined afterwards. Thus, the output of an RBF model is given by: In (1), y[k] denotes the output, at instant k, i[k] is the jth input at that instant, w represents the vector of linear weights, C(j) represents the vector (extracted from the C matrix) of the centers associated with hidden neuron j, σ j is its spread, and || 2 denotes the Euclidean distance. RBFs, as well as multilayers perceptrons (MLP), B-spline neural networks and others, are nowadays called shallow ANNs [24], in contrast with deep architectures, employing a larger number of layers, such as deep neural networks, deep belief networks, long-short term memory neural networks (LTSM) and convolutional neural networks (CNN).
Yang and co-workers [25] proposed a hybrid scheme, involving classification, training, and forecasting stages. In the classification stage, self-organizing maps and learning vector quantization networks classify the collected historical data of PV power output. The training stage employed the support vector regression to train the input/output data sets for temperature, probability of precipitation, and solar irradiance of defined similar hours. In the forecasting stage, the fuzzy inference method was used to select an adequately trained model for an accurate forecast. This scheme is used for one-day ahead hourly forecasting of PV output.
Leva and co-workers [12] used an MLP to forecast the hourly day-ahead PV power, based on weather forecasts, and power and irradiance measurements. Having available a 72-h weather forecast, in each day x an hourly forecast for day x + 1 was produced, between sunrise and sunset.
In [26], a 72 h deterministic and probabilistic forecast of PV power output was developed using MLPs combined with an analog ensemble, using as inputs forecasts of global irradiance, cloud cover and atmospheric temperature, obtained from a numerical weather prediction model, as well as solar azimuth and elevation. Values of RMSE/C (C is the rated power of the PV system) in the range of 8 to 8.7% were obtained.
Rana and co-workers [27] provided twelve different forecasts, from 5 to 60 min ahead, obtained with 12 different prediction models. The authors employed MLPs, arranged in different ways: using single and ensemble models, and using only the past PV power measured values (univariate) and combined with solar irradiance, atmospheric temperature and relative humidly, and wind speed (multivariate). A correlation-based feature selection algorithm was employed for each model, to select the inputs from a set of lags until oneweek before, i.e., 12 × 24 × 7 = 840 lags for a univariate model, and 5 × 840 = 4200 lags for multivariate models. In [28], PV power was forecasted for a solar power plant located at the Applied Science Private University, in Jordan. One day-ahead PH was considered, using the measured irradiance values of the five previous days, obtained at the same hour. Employing MLPs, R-squared values for the training, testing and validation sets of 0.97 were achieved.
Addressing now deep architectures, in [29] one day-ahead forecasting was obtained using generative adversarial networks and convolution neural networks to classify weather types, and with trained forecasting PV power models for each type of weather. Unfortunately, no information about the quality of the forecasts was supplied. Zhou et al. [30] employed two long-short-term memory neural networks for PV power output forecasting and atmospheric air temperature, obtaining forecasts for PHs from 7.5 min up to 1 h, using 7.5 min sampling intervals.
Wang and co-workers [31], developed three forecasting models: convolutional neural network, long-short-term memory neural network and an hybrid model (combining the other two models). They used the PV power output itself, the atmospheric air temperature, and the global solar irradiance. MAPE values varying from 2.2 to 11.2% were obtained for a one-step-ahead (5 min) prediction.
Li and co-workers [32] propose a hybrid deep learning approach based on CNN and LSTM for PV output power forecasting. The CNN model is intended to discover the non-linear features and invariant structures exhibited in the previous output power data, thereby facilitating PV power prediction. The LSTM is used to model the temporal changes in the latest PV data and predict the next time step's PV power. Then, the prediction results in the two models are comprehensively considered to obtain the expected output power.
Hossain and Mahmood [33] employ an LSTM neural network, using as exogeneous variable a synthetic weather forecast. A synthetic weather forecast is created for the targeted PV plant location by integrating the statistical knowledge of historical solar irradiance data with the publicly available sky forecast of the host city. They compared their proposal with recurrent neural networks, generalized regression neural networks and extreme learning machines, for different intraday horizon lengths in different seasons.
The primary purpose of improving the accuracy of solar power forecasts is to reduce the uncertainties related to this type of intermittent energy source, resulting in safer and easier energy management. As solar penetration increases in the energy portfolio, the impact of incorrect forecasts in the grid and HEMS can be larger [1].
A particular model's performance and accuracy can be assessed via several metrics, which allow the comparison between different models and locations. Each one focus on a certain aspect of point distribution. Thus, there is no unique metric valid for all situations; instead, each one adds some information about the model's accuracy. In the bibliography, several metrics can be found [34], although there is a group that is more commonly used, such as mean absolute error (MAE) (2), mean relative error (MRE) (3), normalized MAE (NMAE) (4), mean absolute percentage error (MAPE) (5), root mean square error (RMSE) (6) and coefficient of determination (R 2 ) (7). For this reason, these performance criteria will be used in this work, in Case Study 2.
Inventions 2021, 6, 12 6 of 23 In the previous equations, n is the number of samples, y t is the measured tth value, y t is the predicted value, C is the rated power of the PV system, and r is the range of the measured variable.
Independently of the metrics used to assess the proposed models' performance, some other factors hamper comparisons among studies. According to [1], they are climatic variability, day/night values and normalization, sample aggregation, testing period, and specific system attributes. Table 1 presents the summary of pertinent comparable studies in terms of the method used, the prediction horizon, sampling time, input variables, and the model's accuracy. Although the works described above are general techniques for short-time and intraday PV power forecasting, the majority cannot be applied for model-based predictive control and HEMS scheduling. The reason is that MBPC uses predictive models, that should output the modelled variable's forecasts for each step-ahead within the PH considered, i.e., provide multi-step-ahead forecasting. This type of forecasts can be achieved in a direct mode, by having several one-step-ahead forecasting models, each providing the prediction of each-step ahead within PH. This is the approach followed in [27,32]. An alternative, which is the one followed in this work, is to use a recursive version. In this case, only one model is necessary, and for each step within the PH, the inputs change, eventually employing predictions obtained in previous steps. This way, the RBFs models described in (1) are used as NAR (nonlinear autoregressive) models or NAR models with eXogenous Inventions 2021, 6, 12 7 of 23 inputs (NARX). Denoting as y the modelled variable, and considering only one exogenous input, v, for the sake of simplicity, the estimation (ŷ), at instant k, can be given as (8): In (8), f(.) represents the RBF model (1), which means that its arguments (the delays of y and v) represent the network input vector, i[k]. As the objective is to determine the evolution of the forecasts over a prediction horizon, (8) is iterated over that horizon. For k + 1, we shall have: Depending on the indices of the delays, and the steps within the prediction horizon, no measured values may exist for one or more terms in the argument (8). These must be obtained using previous predictions. This way, the computation of the predictions over a prediction horizon PH may require PH executions of the model (8).
Please note that, if multi-step-ahead forecasting is considered (as is the case of this work), the metrics (2-7) can be computed for each step within the prediction horizon. Longer PHs tend to increase forecasting errors, especially under abnormal weather conditions [6], and with an increased interval between samples and with an increasing number of steps. This is more problematic for recursive multi-steps-ahead forecasting methods, as the errors propagate through the steps.

Methodology
The work methodology employed is divided into six parts (shown in Figure 1), which are used for each case study. The first is the case study description, defined in terms of location and weather, building and PV system characteristics and data acquisition system. The second part is constituted by the data set pre-processing to treat outliers and abnormal data and averaging the acquired data to the employed sampling intervals. The third part is the definition of the model minimization objectives, or constraints when appropriate, and variables and lags definition. The data pre-processing output is then used as input for a data selection algorithm, the ApproxHull, described below. The fifth part is the design of one model, or an ensemble of models, involving input selection, topology determination (number of neurons) and parameters estimation. The MOGA framework achieves the model design (please see Section 3.2). The analysis of the forecasts obtained constitutes the sixth part. The fourth and fifth steps are briefly reviewed in the next sub-sections. NILMforIHEM where this work is inserted, and has an extensive acquisition system monitoring different variables related to the energy consumption of the building. Case study 2 uses as input data the PV power consumption, the ambient air temperature and the solar irradiation. At the time of publication, the access to the data is only available to the research group.

Data Set Construction
This work uses the ApproxHull algorithm proposed in [35], to select data for training, testing and validation data for the artificial neural networks design. ApproxHull is an incremental randomized approximate convex hull (CH) algorithm that selects the points involving whole data points. It is applicable to high dimension data, treating memory and time complexity efficiently. The convex hull vertices obtained are compulsorily introduced in the training set so that the model can be designed with data covering the whole operational range.
A pre-processing phase is performed on the original data set before applying the convex hull, normalizing the dimensions in the range of (-1, 1). Very briefly, it starts with an The case studies are described in Section 4 (please see Sections 4.1 and 4.2). In summary, two case studies were used to develop this work concerning the short-term forecasting. The first is composed by the dataset of the Honda Smart Home US, developed by Honda, which makes the data publicly available each six months. The input data used for the first case study are the PV power output, ambient air temperature, PV panel temperature and solar irradiance. The detailed data from Case Study 1 may be found in its website, as referenced in Section 4.1, as well as many other variables measured in that household. The second case study is composed by the dataset of a typical Portuguese household, located in Algarve, Portugal. This house is the main case study of the project NILMforIHEM where this work is inserted, and has an extensive acquisition system monitoring different variables related to the energy consumption of the building. Case study 2 uses as input data the PV power consumption, the ambient air temperature and the solar irradiation. At the time of publication, the access to the data is only available to the research group.

Data Set Construction
This work uses the ApproxHull algorithm proposed in [35], to select data for training, testing and validation data for the artificial neural networks design. ApproxHull is an incremental randomized approximate convex hull (CH) algorithm that selects the points involving whole data points. It is applicable to high dimension data, treating memory and time complexity efficiently. The convex hull vertices obtained are compulsorily introduced in the training set so that the model can be designed with data covering the whole operational range.
A pre-processing phase is performed on the original data set before applying the convex hull, normalizing the dimensions in the range of (-1, 1). Very briefly, it starts with an initial convex hull (the maximum and minimum of each dimension), and subsequently, the current convex hull grows by adding the new vertices into it. Then, it generates a population of k facets based on the existing convex hull, selecting the furthest points in the current facets' population as new vertices of the convex hull, which are integrated into the current convex hull. A detailed explanation of the convex hull algorithm may be found in [35].

MOGA Design
The model design problem is typically considered as multi-objective optimization, with possible restrictions and priorities. The MOGA design framework is a hybrid of an evolutionary algorithm and a derivative-based algorithm. The evolutionary part searches the number of neurons' admissible space and the number of inputs (lags for the modelled and exogenous variables) for the RBF models. For this reason, the structure of the chromosome includes the number of neurons and pointers for the admissible delays of the endogenous and exogeneous variable(s), if available. Details on the MOGA method were previously published in [35]. Before being evaluated in MOGA, each model has its parameters determined by a Levenberg-Marquardt (LM) algorithm [36,37] minimizing an error criterion that exploits the linear-nonlinear relationship of the RBF NN model parameters [38,39]. Basically, denoting by X the input matrix, and by v and w the vectors of nonlinear parameters (C and σ) and linear output weights, respectively, the model output vector, y, can be given as: Denoting by t and e the target and error vector, the training criterion usually employed is: Independently of the value of the nonlinear parameters, the optimal value of the linear parameters can be obtained as the least squares solution: If this value is incorporated in (11), a new criterion (independent of the linear parameters) is obtained: Any derivative algorithm can be employed to minimize (13). In our case, the LM method is used. The initial values of the non-linear parameters are chosen randomly, or with the use of a clustering algorithm, w is determined as a linear least-squares solution, and the procedure is usually terminated using the early-stopping approach within a maximum number of iterations.
The framework requires the data used to develop the models to be divided into three data sets: a training set, used to estimate the model parameters, a testing set, for terminating the training, and a validation set, for comparing the performance of the models obtained by executing MOGA (as it uses a multiple objective formulation, its results are not a single solution, but a set of non-dominated solutions). The minimization objectives used in this work are the RMSEs of the training set and the testing set, the model complexity (O(µ)) and the forecasting error (ε p ). This last criterion is obtained by summing the RMSEs along with PH (14), where D is a time-series, with p data points, and E is an error matrix (15): In this work, MOGA is executed with 100 generations, population size of 100, the proportion of random emigrants of 0.10 and a crossover rate of 0.70. The admissible range of neurons varied from 2 to 20, while the admissible number of inputs considered ranged from 2 to 20. After execution, the selection of one model from the non-dominated or Pareto solutions is then performed based on the objective values obtained, the RMSE obtained over the validation data set, and the prediction performance (15) over a user-specified period, which may coincide with D.
Typically, two executions of MOGA are performed. By analyzing the results of the first execution the search space can be reduced, by shortening the admissible range of inputs and/or the number of neurons, and/or set restrictions for some objectives, and assigning different priorities. Examples of MOGA application for various applications can be seen in, for instance, [40][41][42][43].
As mentioned above, the output of MOGA is not a single solution, but a set of nondominated models (or preferable models, if restrictions are used). A best model can be selected to represent the results obtained by the problem. However, the non-dominated (or preferable) set of models can also be employed for ensemble averaging the outputs of these models. As the forecasting criterion (15) is not used as a MOGA restriction, models can deliver a very bad prediction performance in a few situations and should be considered outliers among the non-dominated set. This can be solved if the median of the results obtained in the dominant (or the preferable) set, and not their mean value, is used as the ensemble's output. For an example of applying the ensemble approach for forecasting electricity consumption, please see [44].

Honda Smart Home US
The first case study uses data obtained in the Honda Smart Home (HSH) US [45] ( Figure 2). This net zero energy building is located on the West Village campus of the University of California, Davis. It used sustainable construction materials, has a radiant floor night ventilation and a photovoltaic system. Electric appliances and lighting have high efficiency, and the HVAC system employs a ground-source heat pump. The household has a complex home energy management system to control the electric systems. Details about the construction, electric appliances and data acquisition system details can be found on its website [45].

Honda Smart Home US
The first case study uses data obtained in the Honda Smart Home (HSH) US [45] ( Figure 2). This net zero energy building is located on the West Village campus of the University of California, Davis. It used sustainable construction materials, has a radiant floor night ventilation and a photovoltaic system. Electric appliances and lighting have high efficiency, and the HVAC system employs a ground-source heat pump. The household has a complex home energy management system to control the electric systems. Details about the construction, electric appliances and data acquisition system details can be found on its website [45]. The group responsible for the HSH makes available experimental data every six months. Based on the publicly available data, studies were developed, focused mainly on the integration between electric vehicles and the smart home, the home management systems of the HVAC solutions, and construction practices. The present work uses the HSH data for a first design of the forecasting PV power model. In the first step, the exogenous variables to be employed will be selected, and subsequently, the structure of the model will be defined.
To develop the present study, four variables are used from the HSM data set. They are the global solar irradiance, the PV panels' temperature, atmospheric air temperature, and PV DC power produced. Data from 1st January 2016 to 31st December 2018 has been acquired with a sampling time of 1 min and averaged by the authors in 15 min. Not all data was valid for the four variables simultaneously, whether because of lack of data, or wrong measurements. After cleaning, from the 105,217 available samples, 81,504 values were available for model design.
Within the data considered, the air temperature ranged from −1.6 to 42.8 °C, the panel temperature from −6.5 to 84.6 °C, the maximum of solar irradiance was 1127 W/m 2 , and of the PV power 9.3 kW.

PV Power Static Mappings
The first set of experiments was conducted to determine if the PV power was better approximated by the solar irradiance alone (a), combined with the air temperature (b) or with the PV panel temperature (c). Notice that this a simple static mapping and no forecasting is obtained or desired.
As this is a simple problem, MOGA was not employed. Using only the samples with non-null PV power, the 81,504 values were reduced to 50,525 samples. ApproxHull was employed, resulting in 29,008, 9669 and 9671 samples for training, testing and validation. For each input configuration, models were designed where the number of hidden neurons ranged from 2 to 10. For each configuration pair (inputs and number of neurons), five different trainings were performed using a modified version of the Levenberg-Marquardt The group responsible for the HSH makes available experimental data every six months. Based on the publicly available data, studies were developed, focused mainly on the integration between electric vehicles and the smart home, the home management systems of the HVAC solutions, and construction practices. The present work uses the HSH data for a first design of the forecasting PV power model. In the first step, the exogenous variables to be employed will be selected, and subsequently, the structure of the model will be defined.
To develop the present study, four variables are used from the HSM data set. They are the global solar irradiance, the PV panels' temperature, atmospheric air temperature, and PV DC power produced. Data from 1st January 2016 to 31st December 2018 has been acquired with a sampling time of 1 min and averaged by the authors in 15 min. Not all data was valid for the four variables simultaneously, whether because of lack of data, or wrong measurements. After cleaning, from the 105,217 available samples, 81,504 values were available for model design.
Within the data considered, the air temperature ranged from −1.6 to 42.8 • C, the panel temperature from −6.5 to 84.6 • C, the maximum of solar irradiance was 1127 W/m 2 , and of the PV power 9.3 kW.

PV Power Static Mappings
The first set of experiments was conducted to determine if the PV power was better approximated by the solar irradiance alone (a), combined with the air temperature (b) or with the PV panel temperature (c). Notice that this a simple static mapping and no forecasting is obtained or desired.
As this is a simple problem, MOGA was not employed. Using only the samples with non-null PV power, the 81,504 values were reduced to 50,525 samples. ApproxHull was employed, resulting in 29,008, 9669 and 9671 samples for training, testing and validation. For each input configuration, models were designed where the number of hidden neurons ranged from 2 to 10. For each configuration pair (inputs and number of neurons), five different trainings were performed using a modified version of the Levenberg-Marquardt algorithm [46], with different initial values of the RBF centers determined using an optimal adaptive k-means algorithm [47]. Each training was terminated using the testing dataset for early-stopping. The model that achieved the best compromise between the RMSEs of the training and test sets and the linear parameters norm was selected for the corresponding configuration pair. Figure 3 illustrates the RMSEs obtained in the validation set for the three different model configurations (a, b and c), with neurons ranging from 2 to 10. As can be seen, the models using only solar irradiance as input obtain the worst results, while the models using solar irradiance and panel temperature as inputs achieve the best performance. A compromise solution is to employ solar irradiance and air temperature, as it is more widely available than the panels' temperatures. Figure 3 illustrates the RMSEs obtained in the validation set for the three different model configurations (a, b and c), with neurons ranging from 2 to 10. As can be seen, the models using only solar irradiance as input obtain the worst results, while the models using solar irradiance and panel temperature as inputs achieve the best performance. A compromise solution is to employ solar irradiance and air temperature, as it is more widely available than the panels' temperatures.

PV Power Dynamic Model
Considering the previous results, the exogenous inputs for PV power forecasting will be the solar irradiance and the atmospheric air temperature. As expressed before, recursive multi-step ahead forecasting models will be used. There are, however, two ways of obtaining the forecasts of the PV power: whether we design two forecasting models of the two exogenous variables and pass those values through a static PV model, or we use a single model, which outputs the PV power forecasts. The two approaches were considered, but the latter obtained the best results. Because of that, it will be discussed here.
Samples corresponding to three different periods: P1-lags immediately before the current sample; P2-lags centered on the sample 24 h before; P3-lags centered on the sample one week before will be employed. Considering that each period has three values, the first for PV power, the second for air temperature and the third for solar irradiation, the lags employed are presented in Figure 4. This means that each design sample (out of the 81,504 available) consists of 105 lags (38 for PV power, 29 for air temperature and 38 for solar irradiation. As the largest lag corresponds to 676, and a PH of 48 (12 h) is desired, 676 + 48 = 724 samples cannot be employed, resulting in 80,780 samples available for model design

PV Power Dynamic Model
Considering the previous results, the exogenous inputs for PV power forecasting will be the solar irradiance and the atmospheric air temperature. As expressed before, recursive multi-step ahead forecasting models will be used. There are, however, two ways of obtaining the forecasts of the PV power: whether we design two forecasting models of the two exogenous variables and pass those values through a static PV model, or we use a single model, which outputs the PV power forecasts. The two approaches were considered, but the latter obtained the best results. Because of that, it will be discussed here.
Samples corresponding to three different periods: P1-lags immediately before the current sample; P2-lags centered on the sample 24 h before; P3-lags centered on the sample one week before will be employed. Considering that each period has three values, the first for PV power, the second for air temperature and the third for solar irradiation, the lags employed are presented in Figure 4. This means that each design sample (out of the 81,504 available) consists of 105 lags (38 for PV power, 29 for air temperature and 38 for solar irradiation. As the largest lag corresponds to 676, and a PH of 48 (  The design problem was formulated as minimizing the RMSE of the training and testing set, the model complexity and the forecasting error (15). MOGA obtained 307 nondominated models. The minimum values for the errors obtained are shown in Table 2. Please note that scaled values between −1 to +1 are employed here.  The design problem was formulated as minimizing the RMSE of the training and testing set, the model complexity and the forecasting error (15). MOGA obtained 307 non-dominated models. The minimum values for the errors obtained are shown in Table 2. Please note that scaled values between −1 to +1 are employed here. A model was chosen obtaining a good compromise between the performance criterion, whose structure is shown below (16).
As can be seen, for the three variables, samples belonging to the three periods considered were selected. Some statistics related to the performance of the selected model (16) are shown in Table 3. The evolution of the RMSE over the PH, using scaled data, can be seen in Figure 5. Each value in the prediction horizon axis is equivalent to one step of 15 min-and so, 48 steps ahead in the prediction horizon axis is equivalent to 12 h ahead. The same is valid for the similar figures that follow. Although the errors increase as we are moving up the number of steps-ahead, the prediction error 12 h ahead is just the double of 15 min ahead. This indicates that a larger PH could be considered, if necessary.

Algarve Residence
The residence of this case study ( Figure 6) is in Gambelas, Faro, Algarve, Portugal (37°0′55″ N, 7°56′6″ W). Faro lies 3 m above sea level, and the climate is warm and tem- Although the errors increase as we are moving up the number of steps-ahead, the prediction error 12 h ahead is just the double of 15 min ahead. This indicates that a larger PH could be considered, if necessary.

Algarve Residence
The residence of this case study ( Figure 6) is in Gambelas, Faro, Algarve, Portugal (37 • 0 55" N, 7 • 56 6" W). Faro lies 3 m above sea level, and the climate is warm and temperate, with winters rainier than the summers. The climate here is classified as Csa (Hot Summer Mediterranean Climate) by the Köppen-Geiger system [48,49]. The average temperature in this city is 17.2 • C, and precipitation is about 501 mm.
The residence of this case study ( Figure 6) is in Gambelas, Faro, Algarve, Portugal (37°0′55″ N, 7°56′6″ W). Faro lies 3 m above sea level, and the climate is warm and temperate, with winters rainier than the summers. The climate here is classified as Csa (Hot Summer Mediterranean Climate) by the Köppen-Geiger system [48,49]. The average temperature in this city is 17.2 °C, and precipitation is about 501 mm.
The residence employed is a detached house, with two floors and with 20 different spaces. The household has a PV installation, composed of 20 Sharp NU-AK panels [50], arranged in two strings, each panel with a maximum power of 300 W. The inverter is a Kostal Plenticore Plus converter [51], which also controls a BYD Battery Box HV H11.5 (with a storage capacity of 11.5 kWh) [52]. An intelligent weather station is also installed, which measures solar irradiance and atmospheric air temperature and computes their evolution within a user-specified PH [53]. Many (several hundred) variables are acquired, either with a sampling time of 1 s or 1 min. The data acquisition system is described in the authors' previous work and may be found in [33]. The residence employed is a detached house, with two floors and with 20 different spaces. The household has a PV installation, composed of 20 Sharp NU-AK panels [50], arranged in two strings, each panel with a maximum power of 300 W. The inverter is a Kostal Plenticore Plus converter [51], which also controls a BYD Battery Box HV H11.5 (with a storage capacity of 11.5 kWh) [52]. An intelligent weather station is also installed, which measures solar irradiance and atmospheric air temperature and computes their evolution within a user-specified PH [53]. Many (several hundred) variables are acquired, either with a sampling time of 1 s or 1 min. The data acquisition system is described in the authors' previous work and may be found in [33].
Only three variables are used for the current work: the PV DC power produced, atmospheric air temperature, and global solar irradiance. Data from 19 May 2020 17:37:30 to 31 July 2020 23:52:30 is used, averaged in 15 min steps. As the same values for the three periods of lags are the same as used for the Honda house, the 7034 samples were reduced to 6310 available for model design. Please note that in this second case study a much smaller number of samples is employed, as, in a practical situation, it is not possible to wait three years of data before designing forecasting models. The maximum values of DC power and solar irradiation are 6 kW and 1177 W/m 2 , and the air temperature ranged from 11.7 to 38.2 • C.
In contrast with the model developed for the Honda house where, throughout the prediction horizon, measured data was used for the exogeneous variables, here forecasts of these variables, whenever needed, are obtained from corresponding forecast models. This means that forecasting models for solar irradiance and the air temperature must also be designed.
Common to the three models, the prediction set employed nearly one month of data, from 14 June 2020, 00:07:30 to 12 July 2020, 23:52:30. The design problem will also be formulated, as in the Honda house, as minimizing the RMSE of the training and testing set, the model complexity, and the forecasting error (15).

Solar Irradiance Model
This is a NAR model, which means that no exogenous variables are used in the model. Similarly, as in the Honda house, 38 lags from the three different periods were considered, 20 for P1 and 9 for P2 and P3. After duplicate samples have been removed, 5366 samples were available for model design. AproxHull found 909 convex hull vertices, which were incorporated in the training set. The training, testing and validation sets had 3219, 1073 and 1074 samples, respectively.
MOGA produced 451 non-dominated models. The minimum, mean and maximum values of the root-mean-square errors (ε min , ε, ε max ) for the training, testing and validation sets are shown in Table 4. Please notice that these values are obtained using data scaled between −1 and +1, and that they should be multiplied by 10 −1 . The structure of the selected model is shown below (17).
Some statistics related to the performance of the selected model are shown in Table 5. The evolution of the RMSE over the prediction horizon is shown in Figure 7, and the measured and one-step-ahead air temperature are shown in Figure 8, for four selected days. These days (from 21 June to 24 June) were selected due to the good conditions for solar photovoltaic energy generation in terms of solar irradiation exposure (sunny day, as is possible to see in Figure 8-the solar irradiance reaches nearly 1000 W/m 2 ). The ambient air temperature for these days also were representative for the characteristic Algarve sunny days in Spring and Summer seasons. The evolution of the RMSE over the prediction horizon is shown in Figure 7, and the measured and one-step-ahead air temperature are shown in Figure 8, for four selected days. These days (from 21 June to 24 June) were selected due to the good conditions for solar photovoltaic energy generation in terms of solar irradiation exposure (sunny day, as is possible to see in Figure 8-the solar irradiance reaches nearly 1000 W/m 2 ). The ambient air temperature for these days also were representative for the characteristic Algarve sunny days in Spring and Summer seasons.

Air Temperature Model
The same 38 lags were considered for this model. No duplicate samples were found, and therefore 6310 samples were available for model design. AproxHull found 669 convex

Air Temperature Model
The same 38 lags were considered for this model. No duplicate samples were found, and therefore 6310 samples were available for model design. AproxHull found 669 convex hull vertices and the training, testing and validation sets had 3786, 1262 and 1262 samples, respectively.
MOGA produced 393 non-dominated models. Statistics of the RMSE for the training, testing and validation sets are shown in Table 6. Please notice that these values should be multiplied by 10 −2 . As it is possible to see for the three sets, the maximum values are considerably higher than the average values, which represent a slight increase in relation to minimum values, meaning that in the set few results are nearly maximum values. The structure of the selected model is shown below.
Some statistics related to the performance of the selected model are shown in Table 7. The results for the error parameters are similar in the cases of training, testing and validation. While the first three sets (training, testing and validation) present lower values, the sum of the prediction error is slightly higher than in the previous case (solar irradiance model). The evolution of the RMSE over the prediction horizon is shown in Figure 9, and the measured and one-step-ahead air temperature are shown in Figure 10. The scaled RMSE presented in Figure 9 varies along the prediction horizon from 0.03 to 0.165, which can be considered a very good result. As is possible to see in Figure 10, the one-step predicted values are very close to the measured values.

PV Power Model
ApproxHull found 1343 convex hull points. The training, testing and validation sets consisted of 3786, 1262 and 1262 samples, respectively. MOGA produced 268 non-dominated models. RMSE statistics for the three sets are shown in Table 8. Please notice that the values should be multiplied by 10 −2 . The differences between the maximum and minimum values in the non-dominated set are lower for this model, in comparison with the ambient temperature model, and the average values superior.

PV Power Model
ApproxHull found 1343 convex hull points. The training, testing and validation sets consisted of 3786, 1262 and 1262 samples, respectively. MOGA produced 268 non-dominated models. RMSE statistics for the three sets are shown in Table 8. Please notice that the values should be multiplied by 10 −2 . The differences between the maximum and minimum values in the non-dominated set are lower for this model, in comparison with the ambient temperature model, and the average values superior.

Training
Testing Validation

PV Power Model
ApproxHull found 1343 convex hull points. The training, testing and validation sets consisted of 3786, 1262 and 1262 samples, respectively. MOGA produced 268 nondominated models. RMSE statistics for the three sets are shown in Table 8. Please notice that the values should be multiplied by 10 −2 . The differences between the maximum and minimum values in the non-dominated set are lower for this model, in comparison with the ambient temperature model, and the average values superior. The structure of the selected model is shown below. In contrast with the one selected for the Honda house, the exogenous variables' lags only belonged to Period 1.
Some statistics related to the performance of the selected model are shown in Table 9. The evolution of the RMSE over the prediction horizon is shown in Figure 11, and the measured and one-step-ahead forecasted PV power are illustrated in Figure 12. The scaled RMSE varies from 0.053 to 0.063, which are excellent results. From the 20th step-ahead to the 48th step-ahead, the variation is minimum. The structure of the selected model is shown below. In contrast with the one selected for the Honda house, the exogenous variables' lags only belonged to Period 1.
Some statistics related to the performance of the selected model are shown in Table  9. The evolution of the RMSE over the prediction horizon is shown in Figure 11, and the measured and one-step-ahead forecasted PV power are illustrated in Figure 12. The scaled RMSE varies from 0.053 to 0.063, which are excellent results. From the 20th stepahead to the 48th step-ahead, the variation is minimum. Figure 11. RMSEs of the scaled PV power within the prediction horizon.
As is possible to see in Figure 11, the predicted values are very similar to the measured values in most points, with the most noticeable difference in a valley presented in the third day. The results obtained were considered excellent, which did not justify using a second MOGA execution or the use of a model ensemble.
The one-step-ahead results obtained with the cascade of the three MOGA designed RBF models was compared to a baseline method, a persistent model (20): ( 1) ( ) y k y k + = (20) This model was applied to the whole design data, achieving an RMSE of 0.25, approximately four times higher than the training, testing and validation RMSEs shown in Table 8.

Discussion
For assessing the quality of the forecasting results obtained with the proposed approach, they should be compared with the results achieved with the works referenced in Table 1. It should be noted that a completely fair comparison is not possible, as data used in these works are all different, the prediction horizons are not the same, and the majority of the papers only supply one-step-ahead predictions.
The only work that achieves multi-step forecasts is [27], employing a direct mode. This work uses 12 models that output forecasts with a time-step of 5 min, supplying forecasts between 5 min and 1 h. The MRE reported ranged from 4.2 to 9.3%. The MRE values obtained by our approach are shown in Figure 13. As is possible to see in Figure 11, the predicted values are very similar to the measured values in most points, with the most noticeable difference in a valley presented in the third day.
The results obtained were considered excellent, which did not justify using a second MOGA execution or the use of a model ensemble.
The one-step-ahead results obtained with the cascade of the three MOGA designed RBF models was compared to a baseline method, a persistent model (20): This model was applied to the whole design data, achieving an RMSE of 0.25, approximately four times higher than the training, testing and validation RMSEs shown in Table 8.

Discussion
For assessing the quality of the forecasting results obtained with the proposed approach, they should be compared with the results achieved with the works referenced in Table 1. It should be noted that a completely fair comparison is not possible, as data used in these works are all different, the prediction horizons are not the same, and the majority of the papers only supply one-step-ahead predictions.
The only work that achieves multi-step forecasts is [27], employing a direct mode. This work uses 12 models that output forecasts with a time-step of 5 min, supplying forecasts between 5 min and 1 h. The MRE reported ranged from 4.2 to 9.3%. The MRE values obtained by our approach are shown in Figure 13.
The MRE obtained 1 h ahead in [27] is 9.3%. Our approach achieves a value of 1.3% for the same prediction horizon, which is seven times smaller.
All the other works only supply one-step-ahead forecasts. A comparison of the results obtained with the multi-step-ahead forecast is not fair to our approach, especially when onestep-ahead forecasts correspond to many steps-ahead of the multi-step forecasts. However, this is what will be done as there is no other alternative.
The work presented in [25] achieves an MRE of 3.3% for a PH of 1 day. In our work, the 48-steps forecasts correspond only to 12 h, with a value of 1.5%. Assuming that a 96 steps prediction would follow the same trend of Figure 12, our approach would obtain a smaller MRE.
The work presented in [12] achieves NMAE values between 5.2 to 13.2% for a Prediction Horizon between 25 and 48 h ahead. Considering a rated power of 6 kW, this work's approach obtains NMAE evolutions shown in Figure 14, this is, below 1.5%. Inventions 2021, 6, x FOR PEER REVIEW 20 of 24 The MRE obtained 1 h ahead in [27] is 9.3%. Our approach achieves a value of 1.3% for the same prediction horizon, which is seven times smaller.
All the other works only supply one-step-ahead forecasts. A comparison of the results obtained with the multi-step-ahead forecast is not fair to our approach, especially when one-step-ahead forecasts correspond to many steps-ahead of the multi-step forecasts. However, this is what will be done as there is no other alternative.
The work presented in [25] achieves an MRE of 3.3% for a PH of 1 day. In our work, the 48-steps forecasts correspond only to 12 h, with a value of 1.5%. Assuming that a 96 steps prediction would follow the same trend of Figure 12, our approach would obtain a smaller MRE.
The work presented in [12] achieves NMAE values between 5.2 to 13.2% for a Prediction Horizon between 25 and 48 h ahead. Considering a rated power of 6 kW, this work's approach obtains NMAE evolutions shown in Figure 14, this is, below 1.5%.   The MRE obtained 1 h ahead in [27] is 9.3%. Our approach achieves a value of 1.3% for the same prediction horizon, which is seven times smaller.
All the other works only supply one-step-ahead forecasts. A comparison of the results obtained with the multi-step-ahead forecast is not fair to our approach, especially when one-step-ahead forecasts correspond to many steps-ahead of the multi-step forecasts. However, this is what will be done as there is no other alternative.
The work presented in [25] achieves an MRE of 3.3% for a PH of 1 day. In our work, the 48-steps forecasts correspond only to 12 h, with a value of 1.5%. Assuming that a 96 steps prediction would follow the same trend of Figure 12, our approach would obtain a smaller MRE.
The work presented in [12] achieves NMAE values between 5.2 to 13.2% for a Prediction Horizon between 25 and 48 h ahead. Considering a rated power of 6 kW, this work's approach obtains NMAE evolutions shown in Figure 14, this is, below 1.5%.  Again, speculating that a higher number of steps prediction would follow the same trend of Figure 10, our approach would obtain a smaller NMAE.
The work presented in [30] supplies forecasts of up to 1 h, with MAPE values between 24.7% and 37.8%. The MAPE values obtained with our approach are shown in Figure 15. For a 1-h forecast (four-steps-ahead) the MAPE is 8.0%. tween 24.7% and 37.8%. The MAPE values obtained with our approach are shown in Figure 15. For a 1-h forecast (four-steps-ahead) the MAPE is 8.0%.
In [31], the MAPE criterion is also employed. For a 5 min forecast, MAPE ranges between 2.2% and 11.2%. In our approach, the MAPE for a 15 min forecast is 8.0%.
The work presented in [33] also employs the MAPE criterion. For summer months, for 6, 12, and 24 h the values obtained were 28.6%, 38.5% and 37.8%. In our approach, 9.1% and 9.4% were obtained for 6 and 12 h ahead. Work [28] shows their results in terms of the R-squared (R 2 ) criterion. For a one-day forecast, the value obtained is 0.97. Figure 16 depicts the evolution of the coefficient of determination within the PH. Speculating again that a higher number of steps prediction would follow the same trend of Figure 16, our approach would obtain a much higher R 2 value.  In [31], the MAPE criterion is also employed. For a 5 min forecast, MAPE ranges between 2.2% and 11.2%. In our approach, the MAPE for a 15 min forecast is 8.0%.
The work presented in [33] also employs the MAPE criterion. For summer months, for 6, 12, and 24 h the values obtained were 28.6%, 38.5% and 37.8%. In our approach, 9.1% and 9.4% were obtained for 6 and 12 h ahead.
Work [28] shows their results in terms of the R-squared (R 2 ) criterion. For a one-day forecast, the value obtained is 0.97. Figure 16 depicts the evolution of the coefficient of determination within the PH. Speculating again that a higher number of steps prediction would follow the same trend of Figure 16, our approach would obtain a much higher R 2 value.
Again, speculating that a higher number of steps prediction would follow the same trend of Figure 10, our approach would obtain a smaller NMAE.
The work presented in [30] supplies forecasts of up to 1 h, with MAPE values between 24.7% and 37.8%. The MAPE values obtained with our approach are shown in Figure 15. For a 1-h forecast (four-steps-ahead) the MAPE is 8.0%.
In [31], the MAPE criterion is also employed. For a 5 min forecast, MAPE ranges between 2.2% and 11.2%. In our approach, the MAPE for a 15 min forecast is 8.0%.
The work presented in [33] also employs the MAPE criterion. For summer months, for 6, 12, and 24 h the values obtained were 28.6%, 38.5% and 37.8%. In our approach, 9.1% and 9.4% were obtained for 6 and 12 h ahead. Work [28] shows their results in terms of the R-squared (R 2 ) criterion. For a one-day forecast, the value obtained is 0.97. Figure 16 depicts the evolution of the coefficient of determination within the PH. Speculating again that a higher number of steps prediction would follow the same trend of Figure 16, our approach would obtain a much higher R 2 value.  Work [32] uses a direct mode to supply multi-step ahead PV forecasts with a prediction horizon of 3 h, with steps of 15 min. Different models are designed for each season. For spring, summer, fall and winter, the coefficient of determination ranges from 0.998 to 0.918, 0.999 to 0.964, 0.996 to 0.907 and 0.997 to 0.905, respectively, between 15 min and 180 min ahead. Our second case uses data from spring and summer, with considerable better results.
To conclude the discussion, the authors would like to highlight that this predictive model is only one of the models necessary for our HEMS. Models of electric consumption were developed by the authors using the same methodology [30,33,44] and will be integrated into a MBPC HEMS scheme as the one described in [54]. Furthermore, in future studies the effect of dust on the PV panels will also to be considered, following the guidelines of [55], where an experimental analysis was developed for different dust types to evaluate their impact on the power output of the modules.

Conclusions
This work focused in PV electric power generation, designing for this purpose multistep one-step-ahead forecasting models. The main goals were to improve the current PV power forecasting accuracy and, envisioning possible commercial use, to obtain the forecasts using only a few months of data. Two case studies were employed, the first using data from a specially built house in the USA, and the second based on a typical residence in Portugal's Algarve region.
Excellent multi-step-ahead PV power forecasts were obtained with the methodology proposed. One-step values of MRE and NMAE of 1.25%, 8.0%, MAPE of 9.4%, R-squared of 0.994 were achieved, typically much better than results available in the literature. These were achieved using only a few months of design data were needed, in contrast with other works available in the literature.
In spite of the excellent quality of the models obtained, one should be aware that if this (fixed) model would be employed throughout the year, with different weather conditions, the forecasting performance would decrease. Fortunately, this can be alleviated if an adaptive mechanism is employed. This will shortly be applied to the PV power forecasting model.
The authors also note that this predictive model is only one of the models necessary for the HEMS that will be designed in the project in which this work is inserted. Models of electric consumption were also developed using the same methodology and in future works, they will be integrated into a MBPC HEMS scheme.