Hybrid PV Power Forecasting Methods: A Comparison of Different Approaches

: Accurate photovoltaic (PV) prediction has a very positive effect on many problems that power grids can face when there is a high penetration of variable energy sources. This problem can be addressed with computational intelligence algorithms such as neural networks and Evolutionary Optimization. The purpose of this article is to analyze three different hybridizations between physical models and artiﬁcial neural networks: the ﬁrst hybridization combines neural networks with the output of the ﬁve-parameter physical model of a photovoltaic module in which the parameters are obtained from a datasheet. In the second hybridization, the parameters are obtained from a matching procedure with historical data exploiting Social Network Optimization. Finally, the third hybridization is PHANN, in which clear sky irradiation is used as an input. These three hybrid methods are compared with two physical approaches and simple neural network-based forecasting. The results show that the hybridization is very effective for achieving good forecasting results, while the performance of the three hybrid methods is comparable.


Introduction
In the last 20 years, the penetration of renewable energy sources (RESs) in energy systems around the world has progressively increased due to the rise of environmental concerns and governmental policies. Of the different RESs, the worldwide growth of photovoltaic (PV) technologies has been close to exponential [1].
The most important and challenging problem arising from the great penetration of PV in electrical systems is the high level of variability in the power supplied. In fact, this strictly depends on local weather conditions, such as cloud cover, temperature, wind speed and atmospheric aerosol levels. The resulting uncertainty and variability of the PV power profile create various problems for the management of the electricity grid. First, large frequency oscillations can be induced by abrupt changes in power; secondly, in the case of the high penetration of renewables, reverse active power flows may occur in the medium-voltage distribution power supply, or even in the high-voltage transmission line. Finally, the high penetration of PV increases the costs of the allocation of the spinning reserve, ancillary services and the energy planning of the dispatchable generators [2].
For all these reasons, highly accurate photovoltaic power prediction systems are required to optimize the management of the electricity grid from both a technical and economic point of view, without reducing energy reliability and quality.
Computational intelligence (CI) techniques, such as Artificial Neural Networks (ANN) and Evolutionary Algorithms (EAs), have become very popular approaches to problems such as the modeling of non-linear systems, parameter identification, design and control optimization of electric systems and long-term to short-term forecasting.
In [3], the authors exploit a well-known Evolutionary Algorithm, namely the Firefly Algorithm (FA), for the optimal sizing of stand-alone photovoltaic systems. System Algorithm, namely Social Network Optimization. Finally, the third hybridization method is PHANN, in which clear sky irradiation is used as the input for of the neural network.
The structure of the paper is the following: in Section 2, the physical and hybrid models are analyzed. In Section 2.1, both the physical model and the two possible matching procedures are analyzed, while in Section 2.2, the PHANN model is briefly introduced and the new system, called P5ANN, is described. The results of the proposed methods are shown in Section 3: at the beginning of this section, the results obtained by the two physical models are presented and discussed, and then the comparison between the two proposed P5ANN and PHANN methods is presented and their specific features are analyzed. Finally, some conclusions are reported in Section 4.

Methods
In this section, the methods analyzed in this paper are presented. In particular, Section 2.1 presents the physical methods and Section 2.2 presents the hybrid methods. The methods here proposed exhibit general valence, even if in Section 3 they have been applied to a specific case study.

Physical Models
The PV cell equivalent-circuit models are the mathematical tools that represent PV cells' electrical behavior in terms of an I-V curve and the resulting power-voltage curve (P-V curve). The equivalent-circuit model of a PV device (for example, a module, a string or an integer for a PV field) is obtained starting from the equivalent circuit of the cell by inserting the series and/or parallel relations, which represent the actual connections of the PV cells in the device.
Various equivalent circuits to PV cells are reported in the literature. They differ in their number of elements and in their circuit topology, and consequently in the number of parameters required to characterize them.
Most of the model parameters depend on irradiance and temperature, so they are related to environmental conditions; i.e., the irradiance on the plane of the PV cell G, the so-called plane of array (POA) irradiance and the PV cell temperature T C resulting from the POA irradiance and the ambient temperature T amb .
The simplest equivalent-circuit model for a PV cell consists of a parallel between a real diode and an ideal current source [24]. The accuracy of this model can be increased by adding the series resistance R S , which mainly represents the resistance at the semiconductor interface, and the shunt resistance R SH , which models high-current paths through the semiconductor along material dislocations and mechanical defects [24]. A further increase in model accuracy is achieved by adding a second diode to the circuit model of the PV cell, leading to the double-diode model [25]. The trade-off between accuracy, equivalent-circuit complexity, parameter calculation and computational burden usually leads to the so-called five-parameter model being preferred, which is the single-diode model including both R S and R SH parasitic resistors. In this work, the single diode five-parameter model is selected to represent the PV module. The equivalent circuit is shown in Figure 1.

Maximum Power Point Calculation
Starting from the PV generator's I-V curve, it is possible to compute the maximum power output. In mathematical terms, the output of the PV generator is the function ( f ): the inputs are the equivalent-circuit model parameters, (p), which in turn depend on the environmental data (w). Consequently, the power at the maximum power point, P MPP can be written as The procedure for selecting the correct set of parameters in reference to a specific environmental condition p is called matching. Usually, standard test conditions (STC, defined as solar irradiation G STC = 1000 W/m 2 , a cell temperature T C,STC of 25°C and an air mass of 1.5 and indicated by the vector w STC ) are taken into account to define the set of reference parameters, indicated by the vector p STC . Different matching procedures have been presented in the literature: they are based on the matching of rated values reported in the datasheet or on the matching with some measurements. In this paper, two methods belonging to both matching procedures, respectively, are considered.
The first method matches the set of reference parameters in order to comply with the reference data reported in the datasheet of the PV module [26][27][28]; then, irradiance and temperature corrections are applied to calculate the set of parameters which represent the actual environmental conditions. Some drawbacks affect this method, such as numerical instability or convergence on unfeasible solutions; i.e., a negative value for a series resistor.
The second method selects the parameters in order to fit the real data of a PV generator. This approach requires a large number of tests, in which the environmental conditions and the PV generator output power should be accurately measured. A proper redundancy in the dataset could make the process robust with respect to some errors in the measurements. Moreover, computational problems are reduced.

Five-Parameter PV Model
The PV cell I-V curve is derived from the solution of the five-parameter equivalent circuit model: the constitutive relation of Kirchhoff's laws and elements leads to the following equation [27]: where R S and R SH are, respectively, the series and shunt resistances, I PV is the photogenerated current, I 0 is the diode reverse saturation current and a is the ideality factor parameter. The ideality factor depends on the PN junction ideality factor (n), on the absolute temperature of the PN junction (T C ), on the Boltzmann constant (k = 1.380649 · 10 23 J/K) and on the electron charge magnitude (q = 1.602176634 · 10 −19 C): The series resistance R S makes the equation implicit: the I-V curve can be easily calculated in a numerical manner, while the closed-form solution is based on the Lambert W-function [29].
The junction temperature depends on the POA irradiance on the cell, ambient temperature, wind speed and direction. The actual cell temperature can be measured by means of a resistance temperature detector (RTD) on the back of the PV module [30]; otherwise, it can be evaluated from thermal models [31].
The easiest and the most widely used thermal model in PV modules is based on the normal operating cell temperature (NOCT), which assumes the difference between the cell and ambient temperature (T amb ) to be proportional to the POA irradiation. Thus, the PV cell temperature can be computed as where the ambient temperature and the POA irradiance in these conditions are T amb@NOCT = 20°C and G NOCT = 800 W/m 2 . The wind speed is assumed to be 1 m/s and without thermal convection on the back of the PV module. The photogenerated current is irradiance and temperature-dependent: where I PV,STC is the reference photogenerated current at STC, and α ISC is the short-circuit current temperature coefficient.
The reverse bias saturation current is temperature-dependent: where k is the Boltzmann constant and E g is the bandgap energy of the silicon, which in turn is temperature-dependent: Finally, the irradiance correction of the shunt resistance is It is assumed that the ideality factor n and series resistance R S are independent of temperature and irradiation [27,28].

Parameter Identification from Datasheet
The matching procedure based on data reported in the datasheet and applied to the five-parameter equivalent circuit model allows the identification of the set of parameters, which are I PV , I 0 , n, R S and R SH , which is referred to as the STC, by solving a set of five independent simultaneous equations. Three specific current-voltage pairs belonging to the I-V curve are usually available from the datasheet: the open circuit voltage (0, V OC ), the short circuit current (I SC , 0) and the maximum power point (I MPP , V MPP ). The evaluation of Equation (2) in these points leads to the following subset of three equations: Two more equations are necessary to calculate the reference parameters. These additional constraints have to be chosen according to the expected properties of the I-V and P-V curves. First of all, the derivative of the P-V curve at maximum power point has to equal zero, leading to the following equation: There are different constraints which can be considered to choose the last equation: they can be the slope of the I-V curve at the short circuit point [32] or the slope of the I-V curve at the open voltage point [33]. The derivatives of the current with the voltage under short-circuit and open-voltage conditions are mainly determined by the shunt resistance R SH and the series resistance R S , respectively: Besides the slope of the I-V curve at specific points, the open circuit voltage or the short circuit current at a cell temperature different from STC [26] can be considered; knowledge of the open circuit voltage and short circuit current temperature coefficients is necessary.
In this work, the derivative of the current with the voltage evaluated at the short circuit current was chosen to complete the set of equations [32], resulting in

Parameter Identification with EAs
The parameters of the model can be obtained in an alternative way through an optimization process which aims to find the best set of parameters to fit the measurement data.
This procedure can be formulated as a minimization problem in the following way: where x is a vector containing the set of tentative parameters, X is the search domain for this problem, and MAE is the mean absolute error function between the measured power data (P m ) and the output power estimated with the fiv-parameter model (P MPP ): where N is the number of samples in the dataset used in the optimization procedure. This minimization problem is characterized by many local minima and, therefore, Evolutionary Algorithms (EAs) are suitable for this application. In particular, in this paper, the matching problem was addressed with Social Network Optimization (SNO); this algorithm has already shown good performance on both standard benchmarks and on a large number of different engineering problems [34].
The procedure adopted in this matching problem is depicted in Figure 2. The algorithm creates, by means of its working principles, a candidate solution x containing the reference parameters. They are used in the thermal and five-parameter models in combination with the weather forecasts (w) containing ambient temperature and irradiance. The physical models take as an output the forecasted power P MPP , which is compared with the measured one (P m ) by means of the MAE error. This is fed-back to SNO and exploited to produce a new population of candidate solutions. The termination criterion used in this work for SNO is the maximum number of cost function calls: in fact, this parameter is proportional to the total computational time required because the self-time of the optimization algorithm is negligible with respect to the cost function computation. The algorithm population is set to 25 individuals and 100 iterations have been done, and thus 2500 objective function calls are performed. Figure 3 shows the convergence curves for 40 independent optimization trials on a matching problem with 20 training days. The low dispersion at the end of the optimization of the thin lines (representing each independent trial) around the thick line (average convergence) shows the robustness of the optimization algorithm in this problem. The number of training days-i.e., the days used in the matching procedure-is an important user-defined parameter. In fact, this tunes the trade-off between the computational time and the final performance. Two aspects should be analyzed: the solution should be robust with respect to the training set selection and the performance on the training set (training performances) should be a good estimation of the performance on a new dataset (testing performances).
In order to analyze the first aspect, the training set size has been changed from 1 day to 260 days, with a denser grid at the beginning. For each training set size, 40 random trials were performed, each with a different selection of days. The results are shown in Figure 4a, where the average value, the first and the third quartiles are reported. Increasing the training set size decreased the standard deviation; a good value for the training set size was found to be above 50 days.  In order to determine the training set size that is also significant for the testing performances, a test very similar to the one proposed above can be used. The result is called a learning curve. In this activity, both the training error, computed on the dataset used in the matching procedure, and the testing error, computed on the entire year of measures, are monitored. Figure 4b shows the learning curves. For each dataset size, 40 independent trials were performed, extracting for each one a different training dataset.
Analyzing the results of Figure 4b, it is possible to see that, with a training set with at least 60 days, the results are quite stable and the training error is a good estimation of the testing set. Combining these results with those shown before, it is possible to say that 60 training days is a good compromise between accuracy and computational time that grows linearly with the increase of the training set size.

Hybrid Models
Forecasting methods, as previously mentioned, cover a wide variety of approaches and therefore can be classified in different ways according to the emphasized features in the model itself. For these reasons, the same technique might belong to different classes, and moreover, the same forecasting method could fall within the intersection of two existing classes. Thus, there are different forecasting method descriptions and classifications available in the literature, some of which divide models between "direct" or "indirect" methods based on whether the target parameter is directly forecasted. "Model-driven approaches" and "data-driven" approaches-or "physical" and "stochastic" forecasting methods, which actually might be synonyms for the previous-allow a third group of models, namely the "hybrid" group, which shares some features in common with the other previously listed classes. It has been proven that, by gathering the strengths of the original methods, the new generated hybrid models have enhanced forecasting capabilities [20,35]. In this particular application, the comparison is performed by means of a hybrid model which aims at combining the five-parameter equivalent model of the PV module and the stochastic Artificial Neural Network, and therefore our model has been named P5ANN. In more detail, the main objective of this work is to provide a fair comparison of the forecasting capabilities of different methods in combination: ANN with both a physical model of clear sky solar radiation (PHANN) and the electrical equivalent model of the PV module (P5ANN), respectively. Figure 5 shows the main scheme of the proposed P5ANN model: the ANN has been trained on the existing historical data of the weather forecasts coupled to the relevant power production of the PV system. In addition, the output power of the relevant five-parameter equivalent physical model is also provided. After the hybrid model has been trained, it is ready to provide the PV output power with respect to the new weather forecast provided for the next 24 h. Due to the stochasticity of the initialization of the ANN parameter, in order to reduce the variability of the model's output, several trials of the 24 h-ahead forecast are produced, and the mean daily profile is finally calculated (ensemble forecast).

Results and Discussion
All the proposed techniques have been tested on real data measured at SolarTech LAB at the Politecnico di Milano [36] (latitude 45°30 10.588 N and longitude 9°9 23.677 E). For one entire year, weather forecasts, weather measurements and power output measurement were acquired, resulting in a database of 267 days, because some data were affected by acquisition issues.
The database used in this work, despite being related to a single location, includes multiple weather scenarios and also takes seasonality into account. Since the database covers a whole year, all possible levels of temperature and irradiation that could be exhibited by the given site are present. In addition, the database also takes shadings into account (both near and far shadings). These characteristics of the database allow the results obtained to have a level of generality, despite the limitations of data availability.
The DC output power of a single PV silicon mono crystalline module was recorded. The module was composed of 60 cells connected in series and three bypass diodes. Its rated power was 285 Wp, and it was mounted with an azimuth of −6°30 (where 0°is the south direction and angles are measured clockwise) and tilt angle of 30°. All the PV module ratings under standard test conditions (STC) and nominal operating cell temperature (NOCT) are reported in Table 1.

Forecasting with Physical Models
In the first analysis, the two physical models were compared. Table 2 shows the comparison between the parameters estimated based on the datasheet information and those found by SNO in the matching procedure.
By analyzing these values, we found that those obtained from the SNO matching procedure differed greatly from those obtained from datasheets: in particular, they lost part of their physical meaning because the optimization procedure choses them also to reduce the error introduced by the weather forecast.  Figure 6 shows the envelope-weight mean average error (EMAE) comparison between the two physical models: the datasheet-based and the SNO-based models. The dashed line is the average EMAE computed for the entire year. The error made by the two physical models was similar and concentrated mainly during the winter. This phenomenon was also due to the fact that the EMAE error index was normalized with respect to the maximum between the expected power and the power actually produced. More details on the cost indexes are reported in Appendix A.  Figure 7 shows the comparison between the two techniques, sorting the days according to the datasheet-based model's EMAE value. Even if the SNO-based model error presented some oscillations, it can be seen in most cases to have been better than the datasheet-based error, confirming the lower average value.
Most of the error was concentrated in a few days in which the forecast was highly inaccurate: only in 90 days out 268 was the error above the average. The EMAE trend is the same as the SNO-based model, showing that this can be due to inaccurate weather forecasts, which is reflected by both models as a power production error. In order to assess this statement, the 30 worst days according to the datasheet-based model's EMAE are plotted in Figure 8: in Figure 8a, the power forecasted by the two physical models is compared with the actual power output; in Figure 8b, the measured and forecasted GHI are compared as a reference for understanding the error sources. In the above figure, it appears that, for almost all of the days reported, the difference between the forecasts of the two models was very low. In fact, most of the error was introduced by inaccurate weather forecasts. In particular, it can be noted that, for both models, the error was concentrated on cloudy days; i.e., days with low irradiance.
Finally, Table 3 shows a numerical comparison between the two models according to average error values. The performance of the SNO-based model was better than the datasheet-based model: in particular, the improvement was mainly focused on the error indices most correlated with the MAE, i.e., the cost function used in the optimization process.

Forecasting with ANN and Hybrid Methods
In this section, the different ANN-based models are compared with respect to the described dataset. Table 4 shows the numerical comparison between these models according to the year's average error values. The bold values are the best achieved.
The ANN-based model achieved the worst performance for all cost indexes compared to hybrid models, which showed quite similar performances. The PHANN model achieved the lowest error value for the normalized mean average error (NMAE), EMAE and objective mean average error (OMAE), while the P5ANN model outperformed PHANN considering WMAE and considering the normalized root mean square error (nRMSE). In particular, for the latter indicator, the SNO-based P5ANN showed the best performance; however, this value may have been affected by the fact that the SNO-matching procedure was performed with MAE.  Figure 9 shows the comparison between the average daily EMAE error of the three hybrid models. Although the error was generally greater during the winter, it was still better distributed than for the physical models ( Figure 6). It can thus be noted that PHANN showed an error trend that differed from the P5ANN models, which had a similar behavior: in particular, PHANN showed some peaks with very high error, while it managed to obtain better results in many days in which the EMAE was less than 30%. A similar comparison is made in Figure 10 with respect to the daily average WMAE error. This error index further highights the peaks present in the winter period, accentuating the differences between PHANN and P5ANN. This fact explains why in Table 4 PHANN's WMAE is significantly worse than the other two hybrid models. To better examine the differences between these error indices, in Figure 11, the WMAE error is represented as a function of the EMAE error. Each point corresponds to the daily error made by one of the three models, which are represented with different colors. In the figure, it can be seen that there are two different relationships between these error indices: for some points, there is an almost linear behavior, while for the most part, the relationship is super-linear. This second type of relationship is responsible for the differences highlighted above: in fact, high errors increase their relevance compared to lower errors. In Figure 12, the distributions of the EMAE (a) and WMAE (b) errors are represented by means of histograms for the PHANN and P5ANN SNO-based models. In these graphs, we see the different distributions of these errors: PHANN has a distribution with a lower median than the P5ANN SNO-based model; however, as already highlighted before, the tails are more pronounced, particularly for the WMAE error. Finally, Figure 13 shows a comparison between PHANN and P5ANN SNO-based models on the 30 worst days of the PHANN model: in Figure 13a, the power forecasted by the two models is compared with the measured power output, and in Figure 13b, the measured and forecasted GHI are compared. Analyszng the 30 worst days of the PHANN model, 23 were in common with the 30 worst days of the datasheet-based model. This means that the error was still mainly due to incorrect weather forecasts and that the learning capabilities of hybrid models are able only to partially compensate for this input error.

Conclusions
The forecasting problem for photovoltaic production can be effectively handled with hybrid models. In this paper, the performance and the specific features of three hybrid models have been analyzed on real measurement data acquired at SolarTech LAB . The breadth of the database considered allows generality to the work carried out not to be lost, although the data are related to a single module. In particular, the characteristics of the measurement system, albeit on a different scale, are similar to the problems that can be found in a larger system.
The first two proposed hybridization strategies are based on the five-parameter physical model of a PV module: the parameters of this model have been obtained firstly from the datasheet and then from a matching procedure with the historical data and an Evolutionary Algorithm named Social Network Optimization. The third hybridization strategy is the already tested PHANN, in which clear sky irradiance has been used as the input for the neural network.
Firstly, the two physical models based on the five-parameter model of a PV modules were compared: as expected, the performance of the Social Network Optimization-based model was better thanks to the training process used for the identification of the model parameters.
Secondly, the hybrid models were compared with a basic ANN approach: the superiority of all the three hybridizations was assessed. The difference in the performances of the three hybrid models were very low and changed according to the specific error index used. This shows the importance of hybridization when applying ANN for power output forecasting, while the physical input used for the hybridization only slightly affects the final performance. This means that the hybridization source can be selected according to the available data of each specific application.
Thirdly, the performance of the PHANN and P5ANN SNO-based models were compared against two different error rates: EMAE and WMAE. From this analysis, this second indicator appears to be more affected by PHANN error peaks; on the other hand, the EMAE cost index differentiates the performances of the days with low error more clearly.

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

Appendix A. Cost Indexes
Several cost indexes have been defined to analyze different aspects of forecasting. All these indexes are based on the definition of the hourly error e h : where P m,h is the average measured power in ab hour and P p,h is the forecasted power. Starting from this error definition, some error indexes have been introduced in the literature to analyze different aspects of the forecasting accuracy.
The simplest error that can be defined is the mean absolute error: where N is the number of hours in the evaluated period. In all the performed analyses, the errors have been computed on a daily basis; thus, N = 24. The cost indexes used in this work are the following: • Normalized mean square error (nRMSE) is widely used as a primary error metric because it underline large errors: The normalized mean absolute error (N MAE % ) is the mean absolute error divided by the module rated power (C): The envelope-weighted mean absolute error (EMAE % ) increases the importance of the error in the morning or evening: The weighted mean absolute error (W MAE % ), in which the normalization has been performed according to the measured power: The objective mean absolute error (OMAE), in which the error is corrected using the irradiation level in clear sky conditions (G CS,h ): The MAE cost index has been neglected in this paper because it is only a normalization of the nRMSE. This can be demonstrated starting from the definition of the nRMSE: It is possible to notice that (P m,h − P p,h ) 2 = |(P m,h − P p,h )| = |e h | (A9) Thus, Using the definition of MAE,