A Comparative Study of Statistical and Machine Learning Methods for Solar Irradiance Forecasting Using the Folsom PLC Dataset

1. Introduction

Production management in modern electrical systems is one of the fundamental pillars for ensuring the stability, efficiency, and sustainability of energy supply. System reliability and cost minimisation are the primary goals. One of the main operational challenges is maintaining a dynamic balance between electricity demand and available generation in real time. Since electricity is not easily storable, a precise production-service balance is essential to avoid both energy waste and risks of blackouts or grid overloads [1].

The increasing integration of renewable energy sources, especially intermittent ones like photovoltaic solar and wind, has profoundly transformed the structure and dynamics of electrical systems. While traditional sources such as thermal, nuclear, and even hydroelectric plants allowed for highly controlled production management, new renewable plants depend heavily on weather conditions and are generally more volatile [2].

However, current regulations in many countries give dispatch priority to renewables, forcing their use over more controllable sources. While environmentally beneficial, this introduces significant complexity in system operation [3].

Photovoltaic solar energy, in particular, has seen exponential growth in recent years. As of 2024, global installed solar capacity exceeded 1,400 GW, with annual production reaching 2,000 TWh—about 7% of global electricity generation [4]. This growth is driven by technological cost reductions, incentive policies, and increasing environmental awareness.

Forecasting production has become central to management. Planning and scheduling of production units are based on demand forecasting, and effective management of this new energy reality also requires forecasting renewable energy production. While traditional models relied on statistical techniques, Machine Learning (ML) approaches have gained significant relevance in recent years. These algorithms have shown superior ability to capture nonlinear patterns and adapt to highly variable contexts.

This article critically reviews the primary forecasting methods applied to electricity production management, with a special focus on ML models. Their performance is evaluated using a widely used dataset in scientific literature: the Folsom photovoltaic plant data [5]. This facility has become a benchmark for validating predictive models due to its data availability and quality. Two research questions are addressed in the article:

• RQ1: Are ML methods more accurate than traditional methods?
• RQ2: Can the use of time series enhance forecast accuracy?

The article is structured as follows: Section 2 presents a literature review related to the main topic. Section 3 explains the methods and materials used for the article. Section 4 presents the results, while Section 5 gathers conclusions.

2. Related Literature

Reliability in photovoltaic (PV) production forecasting is essential for properly functioning the power system. Forecasting errors in Europe range between 15% and 100%, measured as normalised RMSE relative to the mean [6]. A significant portion of this error is influenced by weather conditions, which can introduce up to 35% RMSE in predictions [7].

This forecasting error propagates from generation to the grid, operated by a Transmission System Operator (TSO) or an Independent System Operator (ISO), ultimately impacting distribution networks. Although the economic cost of this error is difficult to generalise, it is estimated to range between 40 and 140 USD/MWh [8]. With the substantial increase in PV generation currently underway, this impact is likely to grow. Consequently, there is a growing interest in improving forecasting accuracy [9].

The fundamental element in forecasting is solar irradiance. Generally, two distinct methodological approaches are employed: those based on time series data and those relying on physical models [10].

Forecasting based on physical models initially estimates irradiance under favourable weather conditions, known as clear-sky conditions. These models use cell temperature to calculate power output. Notable examples include the Nominal Operating Cell Temperature (NOCT) models [11] and the Sandia models [12], developed by Sandia National Laboratories. Comparative studies indicate minimal differences in performance between these models [13]. These models are subsequently adjusted using meteorological data, so that the final irradiance estimate accounts for cloud cover, wind, and other atmospheric conditions [14].

In this context, the integration of satellite imagery with radiative transfer physical models to estimate surface solar irradiance at high spatial and temporal resolution has shown improved results [15]. SoDa (Solar radiation Data) is a platform providing access to solar irradiance databases and estimation models, such as HelioClim and the Heliosat-2 model, developed by MINES ParisTech [16]. Bu, Qiangsheng et al. [17] combine spatiotemporal analysis of satellite images interpreted through convolutional neural networks (CNNs) and LSTM (Long Short-Term Memory) networks to simulate the impact of cloud cover on irradiance.

Forecasting using statistical and time series models is applied to direct irradiance estimation and to predicting meteorological parameters that influence irradiance, thereby complementing physical models [18]. Classical time series methods are common—Singh & Garg [19] and Sapundzhi et al. [20] employ ARIMA-based models, although hybrid models incorporating ML are generally preferred. Despotovic et al. [21] use autoregressive models with transfer learning to forecast PV output in Spain. Torres et al. [22] developed a deep learning-based solar power forecasting system that integrates multiple data sources (meteorological, historical production, satellite, etc.). The model combines CNNs to extract spatial features from meteorological data and LSTM networks to capture temporal dynamics. This integration significantly enhances predictive performance compared to traditional models, even when applied to large datasets [23]. Xu et al. [24] present a hybrid short-term PV output forecasting approach that combines signal decomposition with the XGBoost model (Extreme Gradient Boosting).

In order to promote research in this field, several noteworthy initiatives have been undertaken to share solar production data openly [25]. As the integration of PV systems into electric grids increases, it becomes essential to improve forecasting methods.

A clear example is the dataset released by Pedro et al. [5], commonly referred to as the PLC dataset. This dataset is widely used for benchmarking forecasting models.

Marinho et al. [26] explore short-term solar irradiance forecasting using deep learning techniques (CNN-1D, LSTM, and CNN-LSTM) applied to the Folsom (USA) dataset.

Oliveira et al. [27] develop a novel architecture and use the Folsom data for benchmarking. Yang et al. [28] compare models such as Quantile Regression Forests, Gaussian Process Regression, Bayesian Model Averaging, Ensemble Model Output Statistics (EMOS), and Persistence-based probabilistic models, also using the same dataset. Oliveira et al. [29] apply a Quantum Neural Network (QNN) for forecasting, again leveraging the Folsom dataset.

3. Materials and Methods

3.1. Dataset

This study uses a freely available Dataset obtained from the Zenodo repository under the DOI: 10.5281/zenodo.2826939, commonly called the Folsom PLC Dataset. This Dataset was provided by Pedro et al. [5] and contains detailed measurements from the California Independent System Operator (CAISO) headquarters located in Folsom, CA, USA. The data includes single-minute frequency recordings of global horizontal irradiance (GHI) and direct normal irradiance (DNI), as well as some weather conditions: ambient temperature, relative humidity, wind speed and direction, pressure, etc.

In addition to on-site measurements, the dataset also includes meteorological forecast variables obtained from the North American Mesoscale Forecast System (NAM), sky images and satellite images.

The primary reason for selecting this dataset is its extensive use in the scientific community. The Folsom PV Dataset has been widely studied in the context of photovoltaic power forecasting, and it serves as a benchmark in numerous research articles. Its frequent use in the literature enables consistent comparisons across different predictive modelling approaches, thereby facilitating objective evaluation of model performance. Moreover, its public availability and data quality make it particularly suitable for reproducible and comparative research.

Several studies have utilised this dataset to assess and benchmark solar forecasting approaches. For example, it has served as the basis for evaluating probabilistic forecasting methods using ensemble and hybrid models ([28,30]), as well as for implementing deep learning and quantum neural networks in solar irradiance prediction [31].

3.2. Machine Learning Methods

This study considered all machine learning methods capable of providing regression-based predictions. However, only those that yielded the best results are presented here. We show a brief description of each model.

3.2.1. Neural Networks

Neural network models are based on computational models formed by interconnected functional nodes called neurons [32]. Figure 1 shows a neuron of the network in detail. Each neuron produces an output signal ($O_k$) processing input signals using ($I_i$) through an activation function ($f$).

Each neuron ($k$) is connected via links ($w_i$) called axons, which weigh the input the neuron receives. Additionally, biases ($b_k$) are added to increase model flexibility. Each neuron is connected to several input signals ($n_k$) that could be outputs from other neurons, or the predictor variables.

The activation function introduces non-linearity and can be linear, sigmoid, tanh, ReLU, etc., depending on the network’s purpose. The output of a neuron is calculated as in (1).

$$O_k=\sum _{i=1}^{n_k}{w}_i{I}_i+b_k{O}_k$$

(1)

Neurons are organised into layers as shown in Figure 2. The layers are named input, hidden, and output. A common type is the feedforward neural network, consisting of one input layer, one or more hidden layers, and one output layer. It can be seen how the input variables (predictors, $X_n$) are linked to neurons lying in the input layer, which all must have the same number of observations $t$. This data can either be endogenous or exogenous. The neurons in the output layer provide $h$ steps ahead forecasts of the output variable $\widehat{y}$. The general formulation is shown in (2):

$$\widehat{y}=f(X_1,X_2,\dots ,X_n)$$

(2)

Training involves adjusting weights to minimise the difference between predicted and actual values, typically using gradient descent and a loss function like Mean Squared Error (MSE).

3.2.2. Regression Trees

Regression trees (RT) are a specific application of decision trees (DT) [33], which predict values by recursively partitioning the data based on predictor variables, similarly as regression does. The way the information is graphed reminds one of a tree, where each node splits into branches until reaching a terminal node (leaf) containing the response value. Figure 3 shows a generic representation of the tree.

At each internal node $j$, starting from the root node, a decision is made regarding the path to follow, based on the splitting criterion defined at that node, which determines the behavior of the branch from that point onward, called region ($R_j$). This process is repeated until leaf nodes are reached, at the bottom of the tree. Given a dataset $\left(\mathit{X},y\right)=\left(X_1,X_2,\dots ,X_n,y\right)$ depending on the splitting criteria, a leaf node $j$ is reached, from which the value $y_j$ is obtained as node value. The value at each internal node (${\widehat{y}}_{R_j}$) is computed as the weighted average of the values of the downstream branches (regions, also known as rectangles) from that node.

To construct the tree structure and determine the splitting criteria, algorithms rely on error computation through specific metrics, with the most common one for regression being CART (Classification and Regression Trees). The process of constructing a RT involves determining the optimal number of terminal nodes, $T$, as well as selecting a regularisation parameter $\alpha$ that balances the trade-off between model complexity and data fitting. A larger number of nodes typically allows the model to capture more intricate patterns in the data, but it also increases the risk of overfitting. Conversely, a smaller tree may generalise better but at the cost of reduced accuracy.

A cost-complexity pruning approach is commonly employed to address this trade-off, where the objective is to minimise the function (3).

$$\sum _{j=1}^{\left|T\right|}\sum _{x_i\in R_j}{\left(y_i-{\widehat{y}}_{R_j}\right)}^2+\alpha \left|T\right|$$

(3)

Here, $\left|T\right|$ is the number of terminal nodes in the tree, and $\alpha$ is a non-negative parameter that penalises tree complexity. Optimal values of T and α are typically obtained through cross-validation. This involves partitioning the dataset into training and validation subsets, fitting trees with varying complexity, and selecting the configuration that minimises the cross-validated error. This procedure ensures that the final model achieves a good balance between predictive accuracy and generalisation capability.

3.2.3. Random Forest Ensemble

The combination of simple processes yields remarkable results. This philosophy underpins the foundation of ensemble methods based on DT. These methods enhance prediction accuracy by combining simpler models. Among them, Random Forest (RF) stands out as a robust technique that has consistently delivered strong performance [34].

Initially, the Bagging (Bootstrap Aggregating) technique is applied. See Figure 4 for a better understanding. This method involves dividing the dataset into $m$ training subsets $B_i$, on which a series of sampling operations with replacement are performed. Simple models—RT in this case—are then fitted to each subset. Each sample drawn is independent from the others, which contributes to increasing the variability among the models.

A Random Forest is an ensemble of regression trees, where each tree is trained on a different bootstrap sample of the original dataset, introducing a layer of randomness. Specifically, a random subset of features is selected at each split within a tree, from which the best split is chosen. The purpose of this technique lies in its ability to decorrelate the trees, which significantly enhances the generalisation capacity of the ensemble and reduces overfitting.

The final prediction of a regression model based on RF is generally obtained by averaging the outputs of all individual trees. This aggregation smooths out the variance inherent in single DTs, resulting in more stable and accurate predictions.

3.3. Analysis

Based on the research questions posed, a structured methodology was defined with the primary objective of analysing the behaviour of various ML models. Given that energy production is highly dependent on solar irradiance, the prediction process focused specifically on estimating two key components: DNI and GHI.

The study was organised into two main research branches. The first focused on physical modelling approaches related to cell temperature estimation. The second branch explored time series forecasting based on historical production data.

A hybrid methodology was implemented for the physical modelling approach. Initially, cell temperature was estimated using a clear-sky model. This model was subsequently refined by integrating meteorological parameters using statistical and machine learning methods.

In contrast, the time series approach relied solely on ML models. Historical datasets, including energy production and weather-related variables, were used to train the models and generate forecasts.

4. Results

4.1. Physical Methods with Climate Features

Following the methodology described in [5], forecasts of both global and direct irradiance have been performed using a clear-sky model and by generating predictions for the clear-sky index $k_t$, defined as the ratio between the actual irradiance and the theoretical clear-sky irradiance as seen in (4).

$$k_{\mathrm{t}}^{GHI}={\displaystyle \frac{GHI}{GH{I}_{cs}}}, { k}_t^{DNI}={\displaystyle \frac{DNI}{DN{I}_{cs}}}$$

(4)

On the one hand, predictions are made using only endogenous variables; on the other hand, the model is extended by including exogenous variables.

The procedure involves first computing the clear-sky irradiance values—global horizontal irradiance ($GH{I}_{cs}$) and direct normal irradiance ($DN{I}_{cs}$)—using the Ineichen and Perez model [35], which accounts for site-specific parameters such as atmospheric pressure and air mass. Once the clear-sky values are estimated, the predictive model is used to forecast the corresponding clear-sky index:

$$\widehat{GHI}={\widehat{k}}_{\mathrm{t}}^{GHI}·GH{I}_{cs} , { \widehat{DNI}=\widehat{k}}_t^{DNI}·DN{I}_{cs}$$

(5)

The choice of the physical model was motivated by the need to establish a baseline for comparison with the work presented in [5]. Although other models, such as those referenced, were considered, the minimal differences in performance led to the decision to retain the original approach.

As a result, the estimation of ${\widehat{k}}_t$ becomes the central focus of this analysis. Three prediction horizons have been considered: intra-hour (5 to 30 minutes ahead), intra-day (30 to 180 minutes ahead), and day-ahead (26 to 36 hours ahead). Forecasts have been made for a single time step ahead for each horizon.

The models employed were ordinary least squares (OLS), Ridge regression, and least absolute shrinkage and selection operator (LASSO). Owing to space limitations, a detailed account of the methods is omitted. Readers are referred to [36] for thorough and accessible explanations. In this study, we have utilised the same statistical models and, additionally, tested several ML models: RT, RF, neural networks (NN), LSTM, Support Vector Machines (SVM), among others. Only those models that yielded the best results are reported in this work. The forecasting results are presented in Table 1 for GHI and Table 2 for DNI.

An analysis of both tables indicates that, although increasing the forecasting horizon does not lead to a dramatic rise in RMSE values, the variability of this metric does increase significantly. This pattern is observed for both GHI and DNI. It is also evident that including exogenous variables, such as meteorological data, does not consistently enhance the predictive performance. While such variables contribute to improved results for intra-hour and next-day horizons, they do not offer benefits for intraday forecasts.

However, the comparison between ML and statistical models does not clearly favour either approach. While RF methods tend to show improvements in most cases, the other ML techniques do not consistently outperform the metrics achieved by statistical models, despite optimising their hyperparameters.

It can be concluded that, in this case, the use of ML does not provide a competitive advantage in operational terms: the performance metrics are very similar, while the time investment required to develop the models is significantly greater.

4.2. Time Series Approach

The previous analysis employed a strategy in which regression models were trained to generate predictions for a single future time step. In this approach, the endogenous variables described in [5] —namely $B_{k_t}$, $L_{k_t}$, and $V_{k_t}$—served as input features to inform future forecasts. $B_{k_t}$ calculates the backward average of $k_t$; $L_{k_t}$ stands for the lagged average values and $V_{k_t}$ represents the variability of $k_t$. While this method proved functional, using time series, the present analysis aims to adopt an alternative strategy more closely aligned with time series modeling.

Specifically, the objective is to explore the potential benefits of utilising the full temporal spectrum of the series, capturing its dynamic behaviour over time rather than focusing solely on individual time steps. ML techniques have been selected as the modelling framework to enable this expanded approach, thereby continuing and building upon the exploratory work initiated in the prior analysis.

In this context, the entire time series is utilised, including periods during which no production occurs, or is not expected to occur, since irradiance data occasionally reflect non-zero values under such conditions. The prediction is performed directly on the GHI and DNI values, rather than on derived indices. To capture the seasonality inherent in the time series, and following an approach similar to that proposed in [37], synthetic variables are introduced. Among these, the most notable is the 24-hour lagged version of the target variable, which is incorporated regardless of the forecasting horizon, as it encapsulates daily cyclic behaviour.

To account for potential trends, the series includes a moving average variable computed from the data 24 hours prior. This feature is also incorporated independently of the prediction horizon.

Finally, a hybrid approach has been considered, in which the calculation of $k_t$, as used in the initial analysis, is retained, but enhanced by incorporating the previously mentioned time series features. In this way, the $k_t$ series is treated as a time series, including its seasonality and trend components. The key difference lies in the influence of the physical model, which varies depending on the time of year, while the $k_t$ index typically exhibits highly stochastic behaviour. The forecasts are made thus using (4).

Table 3 compares performance metrics for the prediction of DNI using ML models with a time series approach. Similarly, Table 4 provides the corresponding comparison for GHI.

The analysis of the results indicates that applying a time series model approach with this strategy does not enhance prediction performance; on the contrary, it degrades it. When using the original data directly as a time series, the inherent variability within the series interferes with the model’s predictions, as the model is unable to respond to rapid fluctuations effectively.

However, when a hybrid strategy is applied, the results improve significantly. In this approach, a physical model accounts for the influence of the sun’s position on the panels and the corresponding irradiance. This allows the ML models to focus solely on capturing the meteorological patterns, leading to more accurate and stable predictions.

5. Discussion

To address RQ1: Are ML methods more accurate than traditional methods?, this study used the dataset provided in [5], one of the most widely used benchmarks in the field. A comparison was carried out between traditional statistical models and ML models. While it is acknowledged that both statistical and ML approaches could be further refined and optimised to achieve better performance by developing tailored and highly customised models, the comparative analysis based on standard configurations reveals a clear trend: ML models outperform traditional statistical methods under comparable conditions. This suggests that even without extensive fine-tuning, ML techniques provide a more robust framework for solar irradiance forecasting.

Concerning RQ2: Can the use of time series enhance forecast accuracy?, the study explored a time series modelling approach by incorporating seasonality and temporal patterns present in the data. The findings indicate that traditional physical models are inherently better suited to handle the structure of irradiance data, mainly due to their ability to capture the deterministic components related to solar geometry. However, when a hybrid strategy is adopted—combining a physical model with a ML component—the predictive performance improves significantly. This hybrid approach leverages the strengths of physical models to manage solar position and irradiance incidence, while allowing ML models to focus on capturing meteorological variability. Therefore, while the exclusive use of time series modelling with raw data may not improve accuracy and can even degrade it, using time series techniques within a hybrid framework confirms, with some nuance, that temporal strategies can enhance forecast performance.

6. Conclusions

This study has comprehensively analysed using a widely adopted benchmark dataset in the solar energy forecasting field. The experimental framework has allowed for a fair and insightful comparison between traditional statistical approaches and ML models, considering both direct implementation and time series-based strategies.

The research questions posed at the outset have been effectively addressed. The results demonstrate that ML models provide a superior alternative to traditional statistical methods for the task of solar photovoltaic irradiance forecasting. This is particularly evident when time series strategies are employed, especially in hybrid configurations integrating physical modelling with ML. Such strategies enable better handling of the data’s deterministic and stochastic components, leading to improved prediction accuracy.

Naturally, the conclusions drawn here are specific to the scope and dataset of this study. However, the methodology and insights are transferable to other contexts. Work is already underway to apply this approach to additional datasets, to develop more generalisable conclusions and validate the observed trends across varying geographic and climatic conditions.

Future work will focus on the real-time implementation of similar models to integrate these forecasting strategies into operational systems for solar energy management and optimisation.