Fuzzy Time Series Methods Applied to (In)Direct Short-Term Photovoltaic Power Forecasting

Abstract

Solar photovoltaic energy has experienced significant growth in the last decade, as well as the challenges related to the intermittency of power generation inherent to this process. In this paper we propose to perform short-term forecasting of solar PV generation using fuzzy time series (FTS). Two FTS methods are proposed and evaluated to obtain a global horizontal irradiance (GHI) value. The first is the weighted method and the second is the fuzzy information granular method. Using the direct proportionality of the power with the GHI, the spatial smoothing process was applied, obtaining spatial irradiance on which a first-order low pass filter was applied to simulated power photovoltaic system generation. Thus, this study proposed indirect and direct forecasting of solar photovoltaic generation which was statistically evaluated and the results showed that the indirect prediction showed better performance with GHI than the power simulation. Error statistics, such as RMSE and MBE, show that the fuzzy information granular method performs better than the weighted method in GHI forecasting.

1. Introduction

Solar photovoltaic energy (SPE) is positioned as an energy source that contributes significantly to the diversification of the world’s energy matrix. According to statistics from the International Renewable Energy Agency (IRENA), in 2010 the installed capacity of photovoltaic systems (PVS) was 40 GW and in 2020 this installed capacity worldwide reached 714 GW, which represents an approximate increase of 94% in the last decade [1]. In South America, according to IRENA, in 2010 the installed capacity of PVS was 43 MW and in 2020 this number increased to 12 GW, representing an approximate increase of 99% in the region [1].

In the face of this growth, there are several challenges to consider in terms of the high penetration rates of PVS, being that this type of energy generation varies with the existence of a maximum generation limit that changes over time, from seconds to years [2], which is known as variability. In addition, this limit is not known with perfect precision, which is called uncertainty or error. The movement around the sun generates variability that can be predicted, while the variability associated with clouds can be difficult to predict, as well as the uncertainty due to difficulties in forecasting the behavior of weather conditions.

SPE generation forecasts are fundamental to face the challenges brought by variability and uncertainty. They are also a great tool for the management of electricity grids, their security and the commercialization of solar energy [3]. There are works that apply forecasts to improve the performance of the electricity system that consider demand and prices as variables to be predicted [4,5,6]. Both generators and suppliers that work with PV systems require forecasts for making operational and planning decisions [5]. Given, also, the high dependence of the SPE on weather conditions, its nature is unstable and can affect the reliability and quality of the power grid, causing frequency and voltage fluctuations [7].

Considering this context, the main objective of this work is to apply, investigate and evaluate two fuzzy time series (FTS) methods in the short-term generation prediction of SPE with historical data from a single database, obtained in Florianopolis, Santa Catarina, Brazil. The database was provided by the Photovoltaic Laboratory [8] by personal e-mail communication and can be accessed at the GitHub repository [9]. As reported in [10], there are still few studies in the literature conducted in Latin America countries.

Considering the main objective presented, the specific objectives are defined: (1) to compare the global horizontal irradiance (GHI) forecast accuracy of two FTS methods, the first order weighted method with chronological weights and another higher order method working with fuzzy information granules (FIG), which have different learnings; (2) subsequently, to evaluate the use of FTS in a power simulation using the spatial smoothing method [11]; and (3) transversal to the first two objectives, we intend to compare the performance of the FTS with three different short-term prediction horizons, being 5, 15 and 30 min.

The FTS excels in allowing system flexibility by considering natural circumstances while dealing with vague and imprecise knowledge in time series data [12]. Among the possible tools, the FTS methods are implemented through the Python library pyFTS [13] which develops the steps of the method proposed by [12] for forecasting with FTS. This library is the result of a work of the MINDS laboratory (Machine Intelligence and Data Science laboratory) that researches computational intelligence and machine learning, optimization, data visualization and decision making [14].

Considering one of the most comprehensive literature reviews reported in [15], significant progress has been seen in PV power generation forecasting, especially in recent years with the use of machine learning and deep learning methods. However, relatively few studies applying FTS are observed in this research field. The study described in [16] presents a database of irradiance data recorded at 30 min intervals which is used to generate forecasts using two FTS methods, showing optimal performance when compared to other forecasting methods. Other work has developed a fuzzy logic model for short-term forecasting that considers one hour ahead of solar energy production [17]. In the studies [18,19], fuzzy logic is applied to make short-term load forecasts. More recently, the following highlighted studies are [20,21] where the FTS method is improved with probabilistic forecasting and the information granule method, and is proposed in order to simplify the process with multivariate models.

Publications showing the progress of the study of FTS are found, such as the one in [22] where the non-stationary fuzzy time series (NSFTS) is introduced, which is able to dynamically adapt its fuzzy sets to reflect changes in the underlying stochastic processes based on residual errors. There are works, such as the short-term forecasting method based on the Takagi-Sugeno (T-S) fuzzy model for wind power and wind speed, where the results show that the proposed T-S fuzzy model can effectively improve the accuracy of the short-term forecasting of wind power [23]. It is observed in [21] how it is possible to generate the interval forecasts, from which it is possible to construct a cumulative density function and use it to build the quantile function and probabilistic forecasts with the treatment of stochastic simulations or ensembles. The interval forecasts deal with the drawbacks of point forecasts and this is also discussed in [24], where the probability distribution problem is addressed using the Kernel density estimation. Stacked modules of the deep fuzzy model (DIRM-DFM) for accurate prediction have been found that show the current progress of fuzzy models [25]. Another topic that shows the progress in the study of FTS is dealt with in [26]; here, a hyperparameter optimization method for high order weighted FTS that automates the generation of accurate and parsimonious models using genetic algorithms is presented.

The FTS method is used in conjunction with other prognostic methods and interesting results are evident from this coupling. In [27], the FTS and convolutional neural networks (CNN) are combined for short term load forecasting; they report good efficiency of the method. Similar results are reported in [28], where the proposal is to improve a solar forecasting model based on an artificial neural network (ANN) with fuzzy logic preprocessing. The study reported in [18] highlights the ability of FTS methods to deal with sudden variations in temperature, considering their influence on PVS, in a simple and robust way.

In this context, this study evaluates the performance of the FTS methods combined with the spatial smoothing method to solve the SPE generation forecasting problem considering the spatial dimension and characteristics of a specific photovoltaic system. For this evaluation we use the compilation made by [29] where three criteria are recommended to measure the accuracy of the model: the overall bias, the dispersion and the ability to reproduce statistical distributions. The most recommended metrics to quantify these criteria are, respectively, the mean bias error (MBE), the root mean square error (RMSE) and the Kolmogorov-Smirnov (KS) test. In addition, the coefficient of determination (R²), which expresses the fit of the predicted model data to the original data, will be used [15].

This work presents in Section 2 the theoretical foundation, where the forecasting approaches are detailed and the FTS method is exposed, how the fuzzy relationships were created and, from this, how the forecasts are generated. Section 3 describes the model implementation process, the used database, the hyperparameters implementation and the learning process of each method. Section 4 shows the results obtained from the radiation predictions and the subsequent power simulation, as well as a discussion of the results obtained. Finally, the conclusions are described in Section 5.

2. Theoretical Background

This section details forecasting approaches and time scales and briefly reviews the different forecasting methods in the literature. The FTS method, the process of creating fuzzy relationships and the steps to generate a forecast are outlined.

2.1. Photovoltaic Solar Energy Forecast

Depending on the available data, we initially define the time frame on which the forecasts will be based, i.e., the time scale horizon prediction. In [30], the main types in PVS generation forecasting are discussed; that is, direct and indirect. Indirect forecasting first predicts the solar irradiance and then, using a PV performance model of the plant, obtains the energy produced; direct forecasting directly calculates the plant energy production [30]. It is also useful to distinguish forecasting techniques according to time scales, given the importance of the specific definition of the prediction horizon [31]. Two types of forecasting are specified according to the needs of this work: short term forecasting deals with time horizons from minutes to hours and long term forecasting deals with time horizons from hours to a few days ahead [16].

Forecasting methods and time scales are highlighted first as they are considered cross-cutting in forecasting methods. Several of the most commonly used methods, from statistical methods to those using artificial intelligence, are now enunciated. Among the forecasting methods that are considered relevant to highlight are, initially, the statistical methods, which extract relationships with previous historical data to make forecasts of the future behavior of the plant [30]. Here we include the simplest methods, such as the persistence method, and the autoregressive methods compared in [32], such as ARMA and ARIMA. In terms of statistical methods, the quality of the historical data is essential for a good forecast [15]. There are also physical methods which, unlike statistical methods, use the specific configuration of the system and meteorological information for forecasting, methods based on sky images where cloud cover and optical depth are considered, and methods based on NWP (numerical weather prediction) [33]. Artificial intelligence (AI) based methods, such as artificial neural networks (ANN), that work similarly to the learning process of a human brain, learn the relationship between input parameters and output variables by studying previously recorded data. They do not need characteristic information about the system, which makes the use of ANNs ideal for modeling nonlinear, dynamic, noisy and complex systems [34].

Finally, a particular emphasis is made on the fuzzy time series (FTS) method, the method used in this study. The FTS are conceptualized as being part of the branch of computational intelligence [35] and stand out for allowing system flexibility when considering natural circumstances and dealing with the vague and imprecise knowledge in time series data. In this way, it becomes possible to reduce the problem of uncertainty and variability of solar resource data [21].

2.2. Fuzzy Time Series

The theoretical construction that regards the FTS method and its development from the concepts of fuzzy logic (FL) are presented in this section.

The FTS method was introduced by [12] as a nonparametric method for time series forecasting based on fuzzy logic theory that provides a different representation of the time series. They are noted for allowing system flexibility in considering natural circumstances dealing with vague and imprecise knowledge in time series data [21].

The main difference between a conventional time series and fuzzy time series is that the observations are real numbers in a time series, while fuzzy sets [12] compose the universe of discourse; so, given a time series $Y\in 1$ and its individual values $y\left(t\right)\in Y \mathrm{for} t=0,1,\dots ,T,$ the universe of discourse $\left(U\right)$ is delimited by the maximum and minimum values of $Y$, such that $U=\left[\min \left(Y\right),\max \left(Y\right)\right]$ in which the fuzzy sets $fi \left(t\right), \left(t=1,2,3,\dots \right)$ are defined and $F\left(t\right)$ is the collection of $fi$ and is called the fuzzy time series in $Y\left(t\right)$ [12].

If the universe of discourse $\left(U\right)$ with $U=\left\{u1, u2, \dots , ui\right\}$ where $ui$ are the possible linguistic values of $U$, we then define a fuzzy set of linguistic variables $Ai$ from $U$ where $\mu Ai$ is the membership function of the fuzzy set $Ai$ such that $\mu Ai: U$→[0, 1]. If $ui$ is a member of $Ai$, then $\mu Ai\left(ui\right)$ is the degree of membership of $ui$ a $Ai$[24].

2.2.1. Fuzzy Logic Relationships

With the definition of the data set representing the universe of discourse, a fuzzy logic relationship (FLR) is established that creates the relationship between the input data and the estimated output data. The following is the definition of the FLRs made by [12] for both first-order and higher-order models: Assuming $F\left(t\right)$ is caused only by $F\left(t-1\right)$ then there is a fuzzy relationship between $F\left(t\right)$ and $F\left(t-1\right),$ which can be expressed as:

$$F\left(t\right)=F\left(t-1\right)\circ R\left(t,t-1\right)$$

(1)

“$\circ$” is the operator of the Max-Min function. The relation R is called the first-order model of $\left(t\right)$. If $\left(t\right)$ is caused by more fuzzy sets, $F\left(t-1\right), F\left(t-2\right),$ and $F\left(t-p\right),$($p$>0) simultaneously, the fuzzy relation is of a higher order (order $p$) and is represented by

$$F\left(t\right)=\left(F\left(t-1\right)\times F\left(t-2\right)\times \dots \times F\left(t-p\right)\right)\circ R\left(t,t-p\right)$$

(2)

The main steps proposed by [12] with the FLR generate a complex computational demand for the calculation of each of the fuzzy relations. It is for this reason that [36] proposes an improvement of the algorithm by means of a simplification of the arithmetic operations, being less complex than the Max-Min function, creating groups of fuzzy logic relations (GFLR). Considering $\left(t-1\right)$ = $Ai$ and $F\left(t\right)$ = $Aj$, the relation of fuzzy logic (FLR) can be defined as $Ai$ → $Aj$ where $Ai$ and $Aj$ are called the left-hand side and right-hand side of the FLR. The FLRs with the same left-hand are put together in GFLR. The left-hand groups indicate the input value of the model and the right-hand groups correspond to the outputs that were estimated [36].

2.2.2. Stages of the FTS Algorithm

This section describes the steps to follow for the construction of the FTS models that detail and explain the algorithm used to obtain solar irradiance forecasts and the subsequent simulation of the power values. Figure 1 represents the method proposed by [12] and developed through the Python pyFTS library, with the hyperparameters and parameters relevant to the model.

Data Preprocessing

In this stage the data, which may contain missing fields or incorrect information, as well as repeated measurements, are handled and adjusted. In many time series it becomes necessary to perform such adjustments, e.g., to scale data within an interval, remove trends and/or seasonality and understand the behavior of cycles or randomness, which may interfere in the model validation processes.

The membership function, the number of partitions of the universe of discourse and the partitioning method are considered within the hyperparameters of the model. The hyperparameters of the FTS methods provide more versatility and flexibility, allowing control of the sensitivity of the fuzzification process, the number of rules generated and, mainly, the accuracy of the model [20]. The common hyperparameters of the FTS processes are shown in

Table 1. Model Hyperparameters [24].

Symbol	Parameter	Description
$k \in \mathbb{N}$	Number of partitions	Number of fuzzy sets to be created in the linguistic variable
$\mu :U \to \left[0, 1\right]$	Membership Function (MF)	$\mathrm{Measures} \mathrm{the} \mathrm{membership} \mathrm{value} \mathrm{of} \mathrm{a} \mathrm{value} y\in U$ to the fuzzy set
$\mathsf{\Pi }$	Partitioning method	Defines how the universe of discourse will be divided
$\alpha \in \left[0, 1\right]$	The α–cut	Minimum degree of membership to be considered in the fuzzification process
$\mathsf{\Omega } \in \mathbb{N}$	Model order, number of backlogs	Number of backlogs used in the fuzzy house rule precedent
$L \in \mathsf{\Omega } \times \mathbb{N}$	Temporary backlog index	Index vector with length Ω and$1\le L\left[i\right]<L\left[i+1\right] \mathrm{for} i=0,\dots ,\mathsf{\Omega }$
$\mathsf{\kappa } \in \mathbb{N}$	k-Nearest Neighbor (kNN)	The number of nearest neighbors that the spatial index searches for in FIG during fuzzification

, where it is explained how they influence the model. Table 1 was extracted from [24].

Fuzzification of the Data

Once the partitioning scheme and the number of partitions have been defined, during the fuzzification process each element of the numerical time series $Y\left(t\right)$ will be replaced by the fuzzy set with a maximum membership value, and the fuzzy time series will then be created $F\left(t\right)$ [13]. This step is influenced by the hyperparameter α-cut, which is the minimum membership degree to be considered in the fuzzification process.

Fuzzy Rule Generation and Inference Process

The generation of fuzzy rules has a syntactic component that generates linguistic values and a semantic component that associates each linguistic term to its meaning. These rules depend on the method and its characteristics, as well as on the hyperparameters of the order model (Ω) and the delay indices (L) [24]. The inference process depends on the chosen learning method. The WEIGHTED-FTS and FIG-FTS methods are used for the model to learn and represent the temporal patterns found in the data that have been fuzzified. The goal of this process is to produce a $f\left(t+1\right)$ with the fuzzy sets, to represent the future value that is being predicted $y\left(t+1\right)$ [24].

The methods used are classified as multivariable, and within the group of variables that are defined the endogenous variable is the one to be predicted, also known as the target variable or dependent variable, and the other variables are the exogenous, explanatory or independent variables [24].

WEIGHTED FTS Method

This method is a first-order multivariate method proposed by [37] which applies linear chronological weights and produces more accurate forecasts than the fuzzy time series method proposed by Chen [36] where all FLRs have the same weight during the forecasting process. The method of [37] assigns appropriate weights to the fuzzy relationships. This method, as proposed, is described below.

Assuming the forecast of $\left(t\right) \mathrm{is} Aj1, Aj2, \dots , Ajk$. The corresponding weights for $Aj1, Aj2, \dots , Ajk$ these $w1, w2, \dots , wk$ are specified. However, after forming the matrix of weights with those that were assigned, we have $\left(T\right)=\left[{w^{\prime }}_1, {w^{\prime }}_2, \dots , {w^{\prime }}_k\right]$ which must satisfy the following condition:

$${\displaystyle \sum }_{h=1}^k{w^{\prime }}_h=1$$

(3)

therefore, these weights, $w_1, w_2, \dots , w_k$, must be standardized. The following matrix of weights is thus obtained:

$$W\left(T\right)=\left[{w^{\prime }}_1, {w^{\prime }}_2, \dots , {w^{\prime }}_k\right]=\left[\frac{w_1}{{{\displaystyle \sum }}_{h=1}^k{w}_h}, \frac{w_2}{{{\displaystyle \sum }}_{h=1}^k{w}_h},\dots , \frac{w_k}{{{\displaystyle \sum }}_{h=1}^k{w}_h}\right]$$

(4)

where $w_h$ is the corresponding weight for $Ajh$.

Fuzzy Granular Information Method

The FIG (fuzzy information granular) method is a higher order multivariate method. It works as a wrapper that transforms the real values of a multivariate time series into a fuzzy univariate time series [20]. The resulting time series F is composed of data points $\left(t\right)\in F$representing the sequence of fuzzy information granules ${\mathcal{G}}_i$. Each granule contains a fuzzy set of linguistic variables relative to each fuzzy variable ${\mathcal{V}}_i$ [20]. There is a global linguistic variable which is the union of all of the fuzzy information granules, which, in turn, are the combination of one of the fuzzy sets for each variable, such as ${\mathcal{G}}_i$ which in turn is the combination of one of the fuzzy sets for each variable, such that ${\mathcal{G}}_i=\left\{A_j{\mathcal{V}}_i \right\}, \forall {\mathcal{V}}_i \in \mathcal{V}$ and its membership function is given by $\mu \mathcal{G}i={{\displaystyle \cap }}^{\hspace{1em}}\mu A_j{\mathcal{V}}_i$ where ∩ is the minimum T-norm. The FIG set is indexed by midpoints of its internal fuzzy sets. With the linguistic variable FIG the fuzzification process transforms each multivariate data point $y\left(t\right)\in Y \mathrm{into} {\mathcal{G}}_i$∈ FIG such that $f\left(t\right)={\mathcal{G}}_i$ [20].

For each of these methods for the forecasting process, the input is applied to the model and the output is calculated, which in turn will be the predicted value [12]. The forecast horizon is the number of subsequent values to be predicted, or the number of delays in prediction after the last input value $y\left(t-L\left(0\right)\right),\dots ,y\left(t-L\left(\mathsf{\Omega }\right)\right)$ [20].

Defuzzification Process

In this step we transform the fuzzy forecast values that have been predicted with the definition of the discourse horizon, which are linguistic values, into real numbers [13]. The objective of the process is to transform (t + 1) into an estimated numerical value y(t + 1).

Data Post-Processing

In this last step of the forecasting process, we run the spatial smoothing method proposed in [11] where the key point is the first-order low-pass filter and its respective pole, the value of which is a function of the area of the PVS [11].

This model is based on the direct proportionality of the power of a PVS with the incident irradiance [38] and that the variations in irradiance at a given point, which tend to be smoothed by considering the spatial effect, a phenomenon called spatial smoothing [11]. When the filter is applied to the irradiance predicted by FTS, taking into account the area of the plant, a smoothed value called spatial irradiance $Gs\left(t\right)$ is used to represent the irradiance on a surface, represented by A [11].

The application of the model, as described in [38], begins with obtaining the spatial irradiance time series from the predicted $Gs\left(t\right)$ from the values of the irradiance time series that was forecasted $G\left(t\right)$.

$$\frac{G_s\left(t\right)}{G\left(t\right)}=\frac{1}{\tau s+1}$$

(5)

where s is the Laplace transform variable, t is the time and $\tau$ is the filter time constant, approximated by considering the value of the cut-off frequency, f_c, obtained from the curve fit for the cut-off frequencies of the power spectra of several SFV plants as a function of the area A [Ha] of the system [11].

$$\tau =\frac{\sqrt{A}}{2\pi · f_c}$$

(6)

Once the spatial irradiance is obtained, it is possible to obtain, using Equation (7), a simulated power for the system, by means of the product of the spatial irradiance with the installed power, P*, of the PVS given by [39]. In addition, it is considered in this simulation that the temperature rise of the module reduces its efficiency [39]. Therefore, in Equation (7), a factor is included as a function of the module temperature throughout the day. This factor increases or decreases by 0.4% of the power as indicated in most PVS module data sheets for each degree that deviates from the standard test conditions [39].

$$P_{sim}\left(t\right)=\frac{G_s\left(t\right) ·P^{*}}{G_{STC}}\left[1-0.004\left(T\left(t\right)-T_{STC}\right)\right]$$

(7)

where G_STC $=1000 \mathrm{W}/{\mathrm{m}}^2$is the irradiance at standard test conditions, as well as contact temperature T_STC [39], $Gs\left(t\right)$ is the spatial irradiance and $G\left(t\right)$ is the irradiance that was forecasted.

2.3. Statistical Analysis

Once the prognostic method is defined, the way in which its performance will be evaluated is shown, with the intuition of standardizing the results and making them comparative with other works. For this performance evaluation, we used the compilation made by [29] where three criteria are recommended to measure the accuracy of the model: the overall bias, the dispersion and the ability to reproduce statistical distributions. The metrics recommended to quantify these criteria are, respectively, the mean bias error (MBE); the normalized root mean square error (RMSE), which is used to obtain comparative models where the samples have different sizes, as indicated in [15]; and the Kolmogorov-Smirnov (KS) normality test which is defined as the maximum value of the absolute difference between the two cumulative probability functions [29].

The MBE (Equation (8)) and RMSE (Equation (9)) provide information about an expected range of errors in a given geography or station. Likewise, the RMSE represents well the hourly or sub-hourly dispersion of the values [29]. The Kolmogorov-Smirnoff (KS) test (Equation (10)), is obtained by finding the maximum absolute difference between the cumulative frequency distributions of GHI $\phi \left(c_i\right)$ modeled and $\phi \left(m_i\right)$. observed [40]. In addition, the coefficient of determination (R²) in Equation (11), is used to indicate how close the predicted data are to the test data [41].

$$\mathrm{MBE}=\frac{1}{M}{\displaystyle \sum }_{t=1}^N\left(c_i-m_i\right)$$

(8)

$$\mathrm{RMSE}={\left[\frac{1}{M}{\displaystyle \sum }_{t=1}^N{(c_i-m_i)}^2\right]}^{\frac{1}{2}}$$

(9)

$$\mathrm{KS}=MAX\left|\phi (c_i)-\phi (m_i)\right|$$

(10)

$$R^2=1-\frac{{\sigma }^2\left(c_i-m_i\right)}{{\sigma }^2\left(m_i\right) }$$

(11)

N represents the total number of data, M refers to mean, c_i is the $i-ésimo$ predicted value and m_i represents the $i-ésimo$ observed value [41]. In the equations σ represents the variance in the database.

3. Materials and Methods

In this section we describe the implementation of the analytical model and tools used to perform the experimental evaluation based on the algorithm that is described in Section 2. In this sense, the experiments were conducted in order to compare the two FTS methods, WEIGHTED-FTS and FIG-FTS, through three short-term time horizons of 5, 15 and 30 min. Subsequently, as post-processing, the spatial smoothing method of [11] was used to obtain a generated power simulation.

The procedure performed is shown in Figure 2 to give a general understanding of what is addressed in this section. It should also be noted that the power simulation is a consequence of the irradiance forecast performed through two FTS methods.

3.1. Tools, Technologies and Database

The experiments were performed on a virtual machine created using Google Colab, an open source collaborative programming environment [42], a free tool for writing and running Python code, and the graphics processing unit (GPU) provided by the Nvidia K80s platform with 12 GB of memory.

The irradiance database used to train and test the FTS model and the power data used to validate the low pass filter model in the power simulation are extracted from a 2.2 kWp PVS operating at the “Centro de Pesquisa e Capacitação em Energia Solar Fotovoltaica da UFSC” [8], located in the city of Florianópolis/SC Brazil, where the irradiance sensor is located, along with a pyranometer SMP22 Kipp&Zonen with horizontal configuration and a temperature and humidity sensor (PTB110 VAISALA). It consists of 12 months, from 1 January 2018 to 31 December 2018, for model training, and 12 months, from 1 January 2019 to 31 December 2019, for testing, as observed in

Table 2. Database division and number of records.

Time Horizon	Training (Year 2018)	Test (Year 2019)
5 min	105,121	105,121
15 min	35,041	35,041
30 min	17,521	17,521

. The data were captured by the pyranometer every minute, and the database provided by [9] has filtered data every 5, 15 and 30 min. No missing data were found in the database.

The two-year observations represent a 50/50 ratio for training and testing. One year is considered as training in order to provide learning for the FTS model, and the behavior of the seasonality profiles was the same for the test period [16].

Table 3. Dataset descriptive statistic.

Measurement	Unit	Median	Mean	Maximum Record	Standard Deviation
Global Horizontal Irradiance	W/m2	2.59	179.52	1588	282.38
Air Temperature	°C	21.64	21.38	37.12	4.13
Contact Temperature	°C	25.49	28.39	59.56	8.27
Power	W	0	368.42	2732.9	591.75

shows the summary of the data, describing the units of the variables treated, as well as the statistics of the database, and can be accessed in the GitHub repository [9]. In this database, the full 24 h interval is considered. This choice changes the total number of data and the normalization of the error metrics [15]. This choice has been made, taking into account the threshold value below which solar irradiance data are excluded [15]. This means that in the model training process zero values are included and missing fields are filled with zero in the pre-processing process.

Table 4. Variables used in the model.

Parameter	Variable Description	Type
data	Containing the variables min, hour and month	Input
temp_amb	Ambient Temperature (°C)	Input
GHI	Global Horizontal Irradiance (W/m2)	Input/Output
pow	Real power generated (W)	Input
temp_contact	Module temperature (°C)	Input
modeled_power	Power obtained from spatial smoothing (W)	Output

shows the variables that compose the database and the output of the model. The database is composed of four columns: data containing the seasonality information (where the date, hour and minute of each record are shown) from which the variables min (minute), hour, month; temp_amb (ambient temperature); ghi (global horizontal irradiance) selected as endogenous variables; temp_contact (module operating temperature), which is not included in the database available for the model, is extracted; temp_contact (module operating temperature), which is not included in the database available for the model and thus is calculated by considering the nominal operating temperature of the cell; pow (power generated by the PVS) [43]; and modeled_power (modeled power resulting from the process of applying the low pass filter) as outputs of the method.

Solar irradiance has two seasonality components: annual and daily. These two components can be extracted from the data column. The study conducted in [44] demonstrated that the temperature variable exhibits the highest Spearman correlation in relation to GHI among the several meteorological variables that were analyzed. For this reason, only the ambient temperature is considered, since it interferes with the performance of the solar modules.

3.2. Method and Experimental Setup

The selected hyperparameters are shown, the process for creating the fuzzy and the model learning process is given, as well as the GHI forecast and subsequent power simulation.

3.2.1. Hyperparameters

After the choice of the database on which the irradiance forecast is made, the hyperparameters are defined. Two things stand out, the first of which is the fact that the choice of their value is empirical and data dependent [21]. The second refers to the fact that the model is multivariable, so for each variable it is necessary to define the hyperparameters to train the model [21]. The membership function is chosen as reported in the literature, as in [45]. However, it is noted that the real impact of the membership function on accuracy is low, as demonstrated in [24].

The

Table 5. Configuration of the variables used in the model.

Variable	Membership Function (μ)	Partitions (k)	Minimum Grade of Membership (α-cut)
Minute	Triangular	60	0.30
Hour	Triangular	24	0.30
Month	Triangular	12	0.30
Temperature	Triangular	24	0.30
Irradiance	Gaussian	35	0.25

shows the values selected for each variable. The exogenous variables were assigned the same α- cut and the same membership function (triangular), and in the case of the number of partitions, for the variables related to time we respected their time division and for temperature we selected the number of partitions that was close to its variation value. With the endogenous variable we experimented using a different membership function (Gaussian) in order to better capture the transition between the fuzzy sets; likewise, the value of the minimum degree of membership, α- cut, is a little lower than in the other variables, due to an attempt to include more values when performing the fuzzification process.

As mentioned above, the selection of these hyperparameters comes from heuristic knowledge of the process. The basis for the choice of the different values was also taken with the support of works on solar irradiance forecasting by FTS, as in [21,24,45].

The graphical representation of the membership function chosen for the variables is shown in Figure 3.

The graph of a triangular MF is shown for the exogenous variables and, in the case of the GHI variable—the endogenous variable, which is the variable to be predicted—the Gaussian MF is used and the specific linguistic values that are used are shown.

3.2.2. Fuzzy Rules, Learning and Forecasting

The partitions of the endogenous variable (irradiance) are separated by levels of linguistic values, PP–‘very small’, P–‘small’, M–‘medium’, G–‘large’, GG–‘very large’; each of the linguistic values with seven sub-levels. These will generate rules, as per the example in [45]:

PP5, 6 → P6, M0, M1

This rule can be read as: IF (t − 1) IS Very Small (Sublevel 5) AND y(t) IS Small (Sub-Level 6) THEN y(t + 1) it will be Small (Sublevel 6) OR Medium (Sublevel 0) OR Medium (Sublevel 1).

The model rule creation process described in this paper, with the two methods, the WEIGHTED-FTS of the first order (Ω = 1) and the second order FIG-FTS (Ω = 2) and k-NN two (κ = 2), are the result of the learning process of each of the models using train data.

As represented in Figure 4, the rule generation process of the models WEIGHTED-FTS (Order 1) and FIG-FTS (Order 2) have the format Precedent → Consequent, which comes from a temporal pattern indicating two fuzzy sets that appeared in sequence in the created fuzzy time series [45], where the precedent indicates a set at time t and the consequent set that appeared at time t + 1 [45].

In Figure 4 the difference of the order is shown, along with how the information granules of the FIG method work and how the order of each of the methods is manifested. Recall that the order is the number of lags (past values) that are used by the models, i.e., how much past information is necessary to describe future events [45].

The differences between the WEIGHTED-FTS method (which is of order 1 and performs chronological weightings for the result) and the information granules generated in FIG-FTS (which is a higher order method, both multivariate) are highlighted. Once the learning process of the models is finished, the forecasting process is performed, and here the data from the database corresponding to the test data are used.

Finally,

Table 6. Experimental setup of models, variables and prediction horizons.

Model Development	WEIGHTED-FTS	FIG-FTS
Prediction horizons (minutes)	5, 15 and 30	5, 15 and 30
Input	Data, ambient temperature, GHI, generated power, contact temperature	Data, ambient temperature, GHI, generated power, contact temperature
Hyperparameters for learning models	First order (Ω = 1)	Second order (Ω = 2) and k-Nearest Neighbor (κ = 2)
Partial output	GHI	GHI
Low pass filter parameters	Module temperature, installed peak system power, GHISTC	Module temperature, installed peak system power, GHISTC
Output	Modeled power	Modeled power

shows the experimental setup applied with each of the methods, as well as the respective inputs and outputs.

3.3. Low Pass Filter and Power Simulation

According to [46], photovoltaic energy is approximately proportional to a spatial irradiance profile. Thus, once the irradiance forecast is obtained, we obtain a spatial irradiance profile considering the specifications of the PVS which has 20 modules of 110 kWp organized in four strings of five modules and a 2.5 kW ABB UNO Power Inverter. We proceed to the filter application by means of the butter function of the Python Scipy library. Finally, by applying Equation (7), a simulated power value is obtained for each of the values according to the prediction horizon.

4. Results and Discussion

This section shows the results obtained from radiation predictions and subsequent power simulations. The experimental analysis was carried out with two forecasting methods—WEIGHTED-FTS and FIG FTS—with three short-term prediction horizons of 5, 15 and 30 min, as shown in Table 5. In each method, the month, day, hour, minute and ambient temperature are introduced as exogenous variables and the Global Horizontal Irradiance (GHI) as the endogenous variable. Subsequently, spatial smoothing is applied to the irradiance forecast obtained using the technical and physical specifications of the PVS to obtain the generated power simulation. In the Google Colab environment, the average training time varied from 10 to 20 min, depending on the horizon prediction.

4.1. Statistical Errors

The experiments were performed based on the concepts and steps described in Section 2 and the method presented in Section 3. The results with respect to prediction errors of the two FTS models are gathered in Table 6.

The following statistical metrics were used in this study: MBE to quantify the overall bias and detect whether the model is producing an overestimation (MBE > 0) or underestimation (MBE < 0) of the predicted data and capture the overall trends, rather than variability [15]; RMSE (normalized–nRMSE) as the uncertainty measure [29] that is also very sensitive to large individual errors and captures variability rather than general trends [15]; and the coefficient of determination (R²) which expresses the fit of the predicted model data relative to the original data [41]. In addition, the Kolmogorov-Smirnov (KS) test shows the ability to reproduce statistical distributions [29]. Finally, the KS test assumes that the higher the number of experimental data, the closer the modeled distribution is to the measured distribution [29].

Table 7. Statistical metrics of prediction errors.

FTS Method	Prediction Horizon	MBE	RMSE	nRMSE	R2
Part A (Indirect Prediction)–Irradiance Forecasting with FTS Methods
WEIGHTED-FTS	5 min	4.68	36.32	0.08	0.97
	15 min	3.97	36.59	0.25	0.968
	30 min	7.05	43.98	0.61	0.96
FIG-FTS	5 min	3.67	28.14	0.06	0.98
	15 min	1.82	36.27	0.25	0.97
	30 min	−1.98	53.19	0.74	0.94
Part B (Direct Prediction)–Power Simulation from Irradiance Forecasts with Low Pass Filter
WEIGHTED-FTS	5 min	7.57	143.70	0.33	0.93
	15 min	11.31	160.24	1.11	0.91
	30 min	26.63	197.78	2.76	0.87
FIG-FTS	5 min	5.54	139.12	0.32	0.93
	15 min	6.64	134.40	0.93	0.94
	30 min	6.76	177.75	2.48	0.90

shows the results of the error statistics. Part A of the table shows the irradiance forecast (indirect prediction), in order to compare the accuracy of the two FTS methods, and in part B the metrics shown are used to evaluate the use of FTS in a power simulation (direct prediction) using the spatial smoothing method. The analysis of the different time horizons in each of the FTS methods used is also shown in the two parts of Table 7.

For the KS statistical test, a significance level of α = 0.05, from which a critical value (Vc) is established for each of the forecast horizons, considering the size of the samples that each one offers [40]. It is tested if both samples have the same distribution. For this purpose, a hypothesis test is performed, assuming that both samples have the same distribution, the null hypothesis is accepted; this is against the alternative hypothesis that they are different, in which case the null hypothesis rejected [40]. From the KS test result we observe that for KS < Vc the modeled frequency distribution is similar to the frequency distribution and as such the null hypothesis (H0) is accepted; and, for KS > Vc we have the alternative hypothesis (H1) which states that the measured frequency distribution is not consistent with the modeled distribution, i.e., a poor fit of the predicted data [29]. The alternative hypothesis is equivalent to a rejection of the null hypothesis.

The normality test applied to all samples obtained in both the GHI forecast and the subsequent power simulation is observed in

Table 8. Kolmogorov-Smirnov Test (KS).

Normality Test–KS	Vc	KS	Result
Part A–Irradiance Forecasting with FTS Methods
WEIGHTED-FTS/5 min	0.1979	0.1195	H0 Accepted
WEIGHTED-FTS/15 min	0.2381	0.2325	H0 Accepted
WEIGHTED-FTS/30 min	0.3368	0.2646	H0 Accepted
FIG-FTS/5 min	0.1979	0.1057	H0 Accepted
FIG-FTS/5 min	0.2381	0.1878	H0 Accepted
FIG-FTS/5 min	0.3368	0.2464	H0 Accepted
Part B–Power Simulation from Irradiance Prediction Applying Low Pass Filter
WEIGHTED-FTS/5 min	0.1979	0.3093	H0 Rejected
WEIGHTED-FTS/15 min	0.2381	0.3498	H0 Rejected
WEIGHTED-FTS/30 min	0.3368	0.4786	H0 Rejected
FIG-FTS/5 min	0.1979	0.2768	H0 Rejected
FIG-FTS/5 min	0.2381	0.3146	H0 Rejected
FIG-FTS/5 min	0.3368	0.4082	H0 Rejected

. The KS test determines whether two data sets differ significantly [40]. This test compares the distribution of a set of predicted data with a distribution of observed data and evaluates the differences [40].

4.2. Analysis and Discussion

In this section, the statistical analyses and comparisons of the prediction errors for the two FTS models applied in this work are presented, considering the indirect prediction for the solar irradiance forecast and direct prediction for the generated power simulation.

Initially, as stated in the specific objectives of this work, the accuracy of the global horizontal irradiance forecast of the two FTS methods is evaluated and compared. The FIG-FTS method shows better results in the 5 and 15 min time horizon, and for the 30 min time horizon, the WEIGHTED-FTS method has better results, as shown in Table 7.

The higher order FIG-FTS method performs better when working with smaller time scales, since FIG-FTS performs multivariate forecasting, acting as a multiple-input multiple-output (MIMO) method, in which all variables are both targets and explanatory variables [24]. For this reason, it is considered that in order to perform forecasting with multivariate methods it is necessary to have a database with the minimum redundancy and as much detail as possible, and that the use of order 1 methods, such as WEIGHTED-FTS, is preferable when the database has a larger time scale, which implies a smaller number of data.

The KS test determines whether two sets of data differ significantly [40]. Then, as detailed in Table 8, the KS values are lower for the FIG-FTS method at the 5 and 15 min time horizons, i.e., that the FIG-FTS method is more able to reproduce the observed frequency distributions with smaller time scales. Table 8 shows how the power simulation does not have the ability to reproduce the normal statistical distributions of the data, and as such the table shows the rejection of the null hypothesis.

The forecast result shown in the Figure 5 exemplifies one day of forecasts with the test data relative to 3 May 2019 for the two FTS methods and the three employed horizon predictions.

The comparison of the accuracy of the two FTS methods used for the global horizontal irradiance forecast by horizon prediction is presented in Figure 6 where, graphically, the coefficient of determination (R²) is treated. With the comparison made in Figure 6, the better performance of the FIG-FTS method in the 5 and 15 min time horizon is visible when compared to the WEIGHTED-FTS method, which results in a better R² when the model forecasts in the 30 min time horizon.

The evaluation of the use of FTS in a power simulation using the spatial smoothing method, presented in Section 2 (stage f) of the post-processing of the GHI forecast with FTS, is now presented.

The power was calculated from the irradiance forecasts that were performed with the WEIGHTED-FTS and FIG-FTS methods. The error metrics corresponding to this part of the process that is presented in this work are also shown in Table 7, part B. It is observed that the power results obtained from the forecasts made with the order 2, FIG-FTS method, show better performance than the results obtained with the order 1 method. However, it should be noted that the use of a higher model order does not imply an increase in its performance. Since the higher the order, the more fuzzy sets will be generated, too many fuzzy sets generate an overfitting, causing the model to start learning noise from the data; similarly, sets with lower order generate an underfitting due to oversimplification of the signal [45]. For the proper choice of the order it is possible to perform hyperparameter optimization processes, as shown in [26].

It is also observed, in Table 7 in part B, that the RMSE and MBE values are higher than those presented in the GHI forecasts. This is explained by the fact that, in addition to the error presented in the GHI prediction, the error presented in the power simulation caused by the filter is added, and in addition the correction factor with the contact temperature that decreases the efficiency of the photovoltaic module, as shown in Equation (7). It is relevant to highlight that no studies evaluating, together, the direct and indirect predictions with the application of FTS-based methods at the short-term horizons of this study were observed in the literature.

An example of the power simulation is shown in the Figure 7, where the GHI value obtained with the FIG-FTS method from which the power simulation is performed is shown, and the real observed power value is also shown, with which the results are compared. As already mentioned, a simulation of the generated power values is made from the GHI forecast, using the spatial smoothing method and then passing through the first order low-pass filter [11]. An example of the power simulation is shown in the Figure 7, where the FIG method was used to simulate one day—3 May 2019—of generated power for the 5 min prediction horizon.

The error metrics corresponding to this part of the process that is presented in this work are shown in Table 7. The trend of better results for the FIG method and the 5 min time horizon is shown when analyzing the MBE and RMSE. The same results are obtained in the coefficient of determination analysis applied to the simulation of the generated power, as shown in Figure 8.

Finally, as shown in Table 7, for irradiance forecasts, as the time scale decreases, the accuracy of the method increases. We see that the RMSE and MBE metrics yield better results in the 5 min horizon compared to the 15 and 30 min horizons. These results are expected, since the forecast error measured by nRMSE increases with an increasing forecast time scale [15].

This section presented the comparative performance results of two FTS-based short-term PV power generation prediction methods considering three short-term prediction horizons, including the comparative performance analysis between indirect and direct predictions with FTS.

5. Conclusions

This study developed two multivariable fuzzy time series (FTS) methods and evaluated their use in the indirect prediction of short-term photovoltaic power generation. In addition, the application of the FTS methodology added to the spatial smoothing for power simulation (direct prediction) allows for a controlled experimental setup, enabling the monitoring of the whole process, mainly how the learning of each of the models occurs and the creation of the fuzzy rules, as shown in Figure 4.

Each of the methods used, both WEIGHTED-FTS and FIG-FTS, was assigned the GHI as endogenous variables and the exogenous variables of the minute, hour, day, date and ambient temperature worked through the forecast time horizons of 5, 15 and 30 min. Although there are works where the FTS methods are applied, a contribution of this study is to have used the spatial smoothing method with the application of a low-pass filter, added to the method as post-processing of the data, to simulate power values where the particular characteristics of PVS are considered. For analysis purposes, data obtained from the 2.2 kWp SFV of the solarimetric station installed in the Photovoltaic laboratory of the Federal University of Santa Catarina (UFSC), located in the city of Florianópolis in Brazil, were used.

In the comparative analysis of both methods, it was found that the FIG-FTS method, of the higher order, provides better results in the GHI forecast through the 5 and 15 min horizons, which can be perceived in the statistical results of Table 7 and in the analysis of the KS statistic representing the model’s capacity to reproduce the cumulative distribution function of the observed data. However, the GHI values predicted with the FIG method at the 30 min horizon (the widest of the horizons tested), presents underestimation, i.e., MBE < 0 of the predicted values. At the 30 min time horizon, the best statistics are observed when applying the WEIGHTED-FTS method, as can be seen in Figure 6. This indicates that for longer prediction horizons, first order methods are indicated as more effective.

Once the power simulation was performed with the low pass filter, considering the SFV specifications, it was evidenced that the coupling of the FIG method results in better statistical indexes. It can also be perceived from Table 8 that shorter time horizons, such as 5 min, improve the values obtained from statistics, such as RMSE and KS test values. The statistical analysis shows higher error values, as well as the inability to reproduce the statistical distribution of the samples in the data obtained in the power simulation. This is due to the fact that the application of the low pass filter, by smoothing the fast irradiance peaks, also eliminates irradiance values with large variations typical of this type of measurement.

Although the higher order method has better results, it is necessary to emphasize that increasing the order value of the method indiscriminately does not imply an increase in its performance. This is because the higher the order, the more fuzzy sets will be generated, and too many fuzzy sets generate an overfitting, causing the model to start learning the noise of the data; similarly, lower order sets generate an underfitting, due to the oversimplification of the signal [45].

It is suggested, based on the results, that both FTS methods applied can be used in PV energy generation forecasting in the evaluated short-term horizons, with the best accuracy depending on the prediction horizon. In addition, the direct prediction produced higher errors than the indirect prediction for both FTS methods analyzed. It is highlighted that the implementation of the FTS through the Python library, pyFTS, is a reliable process since it allows access to its source code. In this work, point forecasting is used but its performance can be improved by including, in addition to hyperparameter optimization, interval prediction, with which the fuzzy method can be heuristically simple and fast without generating a large computational demand for the probability density function [24]. Furthermore, to ensure good performance of the FTS forecasting process, the use of a database with good integrity is required. The database provided by the Photovoltaic Laboratory of the Federal University of Santa Catarina was of vital importance.

Finally, another original contribution of this study is that the prediction results have two simultaneous approaches: direct and indirect prediction, since this approach has not been observed on solar photovoltaic energy forecasting literature [15,30,47]. In addition, this study allowed the evaluation of two FTS methods with different experimental setups and their performance when used together with spatial smoothing applied with the low-pass filter. Future work will be directed at comparing the accuracy of FTS methods with benchmark models and other machine learning methods, as well as forecasting with FTS hyperparameter optimization [26] and probabilistic forecasts considering the seasonal components of the input variables [21].