A Comprehensive Tool for Scenario Generation of Solar Irradiance Profiles

Abstract

Despite their positive effects on the decarbonization of energy systems, renewable energy sources can dramatically influence the short-term scheduling of distributed energy resources (DER) in smart grids due to their intermittent and non-programmable nature. Renewables’ uncertainties need to be properly considered in order to avoid DER operation strategies that may deviate from the optimal ones. This paper presents a comprehensive tool for the scenario generation of solar irradiance profiles by using historical data for a specific location. The tool is particularly useful for creating scenarios in the context of the stochastic operation optimization of DER systems. Making use of the Roulette Wheel mechanism for generating an initial set of scenarios, the tool applies a reduction process based on the Fast-Forward method, which allows the preservation of the most representative ones while reducing the computational efforts in the next potential stochastic optimization phase. From the application of the proposed tool to a numerical case study, it emerged that plausible scenarios are generated for solar irradiance profiles to be used as input for DER stochastic optimization purposes. Moreover, the high flexibility of the proposed tool allows the estimation of the behavior of the stochastic operation optimization of DER in the presence of more fluctuating but plausible solar irradiance patterns. A sensitivity analysis has also been carried out to evaluate the impact of key parameters, such as the number of regions, a metric, and a specific parameter used for the outlier removal process on the generated solar irradiance profiles, by showing their influence on their smoothness and variability. The results of this analysis are found to be particularly suitable to guide users in the definition of scenarios with specific characteristics.

1. Introduction

The sustainability objectives set by the European Green Deal require the increasing use of generation systems based on renewable energy sources (RES), as well as the use of electricity as the main energy vector. The pathway to reducing carbon emissions by 2030 will require efforts across society and sectors. With the European Green Deal as the main plan to implement for promoting this change, the European Union (EU) finalized its master program to fight carbon emissions, namely, the “Fit for 55” package [1]. Released in two batches in July and December 2021, the package drafts of EU climate and energy legislation underpin the bloc’s political pledge to cut emissions by at least 55% by 2030 compared with 1990 levels. This target is more ambitious than the previous one of a 40% reduction for 2030 and is the key to achieving carbon neutrality in the EU by 2050. In view of this scenario, solar energy will become one of the key players in the electricity generation sector, thanks to its viability in combating global warming and its effectiveness in reducing pollution caused by fossil-fuel-based generation and diversifying the energy mix to ensure energy security. In particular, the installed capacity of solar photovoltaics (PV) has grown rapidly over the past decade due to the great improvements in PV technology performance, reductions in cost, and the development of efficient business models that have fostered new investments in this technology. This trend is expected to continue in the future, affecting not only large-scale centralized solar farms but, above all, small-/medium-scale PV at the distribution level, where the number of PV applications owned by residential and industrial prosumers (power producers and consumers) will also expand, driven by environmental policies and economic incentives.

From a range of studies, solar PV is expected to contribute 36% to 69% of European electricity consumption by 2050 [2], and its role is predominant in Energy Communities [3]. Due to the variability and intermittent nature of the solar PV output, such a high share of solar PV will impact the overall system costs due to the increase in operating costs and the infrastructure needed. This problem is aggravated by the inaccuracy of the methodologies in modeling renewables’ uncertainties, which are related to the uncertainty of weather conditions for RES [4], which represent a key factor to be properly handled in the smart grid environment. In fact, they may influence how distributed energy resources (DER) are scheduled in the short term to provide the available flexibility for system balancing at all times [5,6,7,8]. If such uncertainties are not identified and handled properly in the operation scheduling of DER, their operation strategies may deviate from the optimal ones by causing a number of issues, such as an increase in operational costs or system stability and security. Modeling RES uncertainties in the stochastic operation optimization of DER is thus extremely important [9,10,11,12].

Several works in the literature deal with different sources of uncertainties, such as renewables, electricity consumption, electric vehicles, etc. [13,14,15,16,17]. The objective of uncertainty-modeling methods is to evaluate the impact of uncertain input parameters on system output parameters. These methods can be subdivided into several groups, as suggested in [18,19]:

• Probabilistic: the probability density functions (PDFs) of the input parameters are used;
• Possibilistic (fuzzy): the uncertainty of the input parameters is modeled with a membership function (MF);
• Hybrid probabilistic and possibilistic: both probabilistic and possibilistic approaches are used;
• Based on Information Gap Decision Theory (IGDT): it measures the deviation of the estimation error;
• Robust optimization: the uncertainty of the input parameters is described using uncertainty sets;
• Interval analysis: the uncertain inputs can assume values in a known interval (similar to the probabilistic approach with uniform PDFs).

The approaches to estimating solar irradiance can be grouped into linear, nonlinear, Artificial Neural Network (ANN)-based, and Fuzzy Logic (FL) techniques [20]. For linear and nonlinear models, the authors have created associations between solar irradiance and other variables, such as meteorological ones [21,22,23]. In ANN-based approaches [24,25], the usual inputs used are geographical coordinates, meteorological data, and information related to the current date and time. In FL approaches, the input to the estimation model is the classified sky condition. Moreover, other authors have applied statistical methods to study the hourly variation in solar irradiance data considering different climatic locations [26] or empirical models to estimate solar irradiance on a monthly basis for different locations [27,28]. In [29], the authors state that FL can yield better estimation results when the data available for estimation are ambiguous and vague.

Among the probability distributions, beta is considered one of the most effective for modeling solar irradiance [30,31,32,33] and is often employed in planning studies related to PV systems [34,35,36]. Other works propose the Weibull distribution for modeling solar irradiance [37].

The contribution of this paper is the presentation of a comprehensive tool to generate solar irradiance profiles using a scenario generation approach and historical data for a specific location. The methodology is general and highly replicable and can thus be applied in several contexts for generating 24 h solar irradiance scenarios, which are useful for the stochastic operation optimization of DER. The tool has been completely implemented in Python, and it is suitable for transformation into a Web Service to generate solar irradiance profiles related to a particular geographical region and a time period of interest.

The historical hourly solar irradiance data were fitted using the beta distribution, and the Roulette Wheel mechanism [38] was used to generate an initial set of scenarios; then, a reduction process based on the Fast-Forward method [39,40] was applied in order to preserve the most representative ones while reducing the computational efforts in the next potential stochastic optimization phase. The generation and reduction phases are ruled by certain parameters, such as the number of regions, a metric, and a specific parameter used for an optional process for outlier removal that can be modified. Moreover, a sensitivity analysis was performed to evaluate the impacts of the variations in these parameters on the solar irradiance scenarios generated. The numerical results of the analysis show that these parameters have a visible effect on the smoothness and variability of the generated scenarios. Based on the current scientific literature, there are no previous works that examined these aspects, thereby highlighting the importance of this study, which could be useful as a guide for tuning the scenario generation process with the aim of obtaining scenarios with certain characteristics. In the case study, the proposed tool is found to be efficient in generating plausible scenarios for daily solar irradiance profiles with 1 h as the time-step to be used as input for DER stochastic operation optimization purposes. Moreover, the high flexibility of the proposed tool allows the estimation of the behavior of DER stochastic operation optimization in the presence of more fluctuating—but plausible—solar irradiance profiles. The current paper extends the results presented in [41] by including additional results in the case study, as well as introducing a new verification system of the plausibility of the reduced scenarios.

In the following, the dataset, data preprocessing, data fitting, and methods for scenario generation and reduction are discussed in Section 2. The results of varying certain key parameters in the scenario generation and reduction processes are presented in Section 3. In Section 4, the sensitivity analysis is discussed, along with the obtained results.

2. Materials and Methods

The tool proposed in this paper is based on a statistical approach used to model solar irradiance based on historical data. By using the Roulette Wheel method [38], an initial set of scenarios is first generated, and then, through a reduction process, the most representative ones are preserved.

A scheme for describing the proposed tool is shown in Figure 1. In particular, the dataset retrieved from the Photovoltaic Geographical Information System (PVGIS) [42] (Section 2.1) is preprocessed by using an optional process for outlier removal and min–max scaling, as described in Section 2.2. For each hour, a data-fitting process (Section 2.3) is performed in order to obtain a probability distribution for each hour. From each hourly probability distribution, a sample is extracted to obtain 24 randomly sampled values, the so-called 24 h solar irradiance scenario (Section 2.4). This process is performed $n_S$ times in order to obtain $n_S$ scenarios. In order to reduce the computational complexity of the following optimization task, the generated scenarios have to be reduced by using an approach that preserves the “information content” present in the original set of scenarios. This is performed by using the procedure described in Section 2.5.

2.1. Dataset Description

The hourly solar irradiance data from 2005 to 2016 for the city of Turin (Italy) were gathered using PVGIS [42]. In order to model both the winter and summer seasons, January and July were selected as months of interest.

The daily patterns of solar irradiance for the days in July and January from 2005 to 2016 are shown in Figure 2 and Figure 3, respectively.

From the figures, it is possible to observe the great variability in each hour, which will be modeled by means of a probability distribution that fits the data.

2.2. Data Preprocessing

The hourly data for each considered month were normalized using min–max scaling in order to map the values observed in the range [0, 1].

To reduce the variability in the observations, the outliers can be removed using the Interquartile Range (IQR) method [43], according to which the values are considered outliers—and hence removed—when they are outside of the following range:

$$\left[Q_1-p\cdot IQR,Q_3+p\cdot IQR\right]$$

(1)

where $Q_1$ and $Q_3$ are the 1st and 3rd quartiles (25th and 75th percentiles), respectively; $IQR=Q_3-Q_1$; and $p$ is a value that permits the expansion or restriction of the range and hence the consideration of fewer or more values as outliers.

However, it is the user’s choice to enable (or not) the outlier removal process, as well as the value of $p$.

2.3. Data Fitting

The normalized hourly irradiance data were fitted using several probability distributions, namely, Weibull, beta, logistic, and arcsine. Among the tested distributions, beta was found to be the best fit for most hours in January and July, confirming what was observed in the relevant literature for solar irradiance data [32,33].

Beta is a continuous probability distribution with support in [0, 1], and its PDF is defined as in [44]:

$$f_S\left(x,a,b\right)=\frac{x^{a-1}{\left(1-x\right)}^{b-1}}{B\left(a,b\right)}$$

(2)

where B(a,b) is the beta function, formulated as:

$$B\left(a,b\right)=\frac{\Gamma \left(a\right)\Gamma \left(b\right)}{\Gamma \left(a+b\right)}$$

(3)

where Γ(x) is the gamma function, formulated as:

$$\Gamma \left(x\right)={{\displaystyle \int }}_0^{\infty }{u}^{x-1}{e}^{-u}du$$

(4)

Different values of a and b allow uniform (a = 1, b = 1), bimodal (a < 1, b < 1), or unimodal (a > 1, b > 1) distributions to be obtained. $n_R$ regions, or bins, are defined to divide the support of the distribution.

Figure 4 shows the fitting results of data related to the hour 12:00 PM in July using the beta distribution and dividing its support into 7 regions.

2.4. Roulette Wheel Method

The Roulette Wheel method [38] was used to extract a set of samples from the fitted probability density functions with their supports divided into $n_R$ regions [45].

The probability of the occurrence of a particular region r ($r\in \left\{1,\dots , n_R\right\}$) at time t (the hour of the day, $t\in \left\{0,\dots , 23\right\}$) can be computed as the product of width $w_{t,r}$ (width of region r at time t) and height $h_{t,r}$ (height of region r at time t): ${\alpha }_{t,r}$ = $w_{t,r}\cdot h_{t,r}$, appropriately normalized as reported in Equation (5).

$${\widehat{\alpha }}_{t,r}=\frac{{\alpha }_{t,r}}{{\sum }_{\rho =1}^{n_R}{\alpha }_{t,\rho }}$$

(5)

The probabilities of the occurrence for all regions of a particular hour were sorted in descending order and cumulated.

In order to sample from the beta distribution with its support divided into $n_R$ regions, it is possible to sample a value from a uniform distribution ($v~U\left(0,1\right)$) and, comparing it with the cumulative probability of occurrence, select one of $n_R$ possible regions that $v$ belongs to (this method is also known as the Inverse Transform [46]).

The central value of the selected region is the value sampled from the beta probability distribution with support divided into $n_R$ regions. The procedure described above was performed for each hour, obtaining 24 sampled values that compose the scenario. To be precise, not all 24 values were sampled in this way because, in some hours, the observed solar irradiance was constantly zero, and this behavior was maintained in the scenario generation process.

The binary variable $W_{k,t,r}$ is used to contain the information on whether region r is selected for a scenario $s^k$ at hour t ($W_{k,t,r}=1$) or not ($W_{k,t,r}=0$).

The main assumption of the proposed model is that, at a specific hour, the solar irradiance is independent of the values observed in the previous hours. With this assumption, the probability of the occurrence of scenario $s^k$, ${\pi }_k,$ is the product of the probabilities of the regions that compose it [45]:

$${\pi }_k=\frac{{\prod }_{t=0}^{23}{\sum }_{r=1}^{n_R}\left\{W_{k,t,r}\cdot {\widehat{\alpha }}_{t,r}\right\}}{{\sum }_{k=1}^{n_S}{\prod }_{t=0}^{23}{\sum }_{r=1}^{n_R}\left\{W_{k,t,r}\cdot {\widehat{\alpha }}_{t,r}\right\}},k=1,\dots ,n_S$$

(6)

The independence assumption can be relaxed, considering that the solar irradiance assumed at a particular hour depends on the solar irradiance values of previous hours. The methodology described here is, however, still valid, and this extension is reported elsewhere.

2.5. Scenarios’ Reduction Process

From an initial set of $n_S$ scenarios with their probabilities ${\pi }_k$ (with $k=1,\dots , n_S$), it is necessary to obtain a reduced set of $n_P$ scenarios to use in the successive potential stochastic optimization phase. The considered reduction method is the Fast-Forward method [39,40]. The algorithm creates a subset of scenarios with the minimum Kantorovich distance from the initial set. After computing the distance (with respect to a metric defined as the mean absolute distance, Euclidean distance, etc.) between all pairs of scenarios, the scenario with the minimum weighted distance (the weights are the probabilities of the occurrence of each scenario) from all of the other ones is selected. The probability of the occurrence of all removed scenarios is absorbed by the preserved scenario that is nearest to them.

The main steps of the Fast-Forward algorithm are described in the Algorithm 1 box [39].

Algorithm 1. Fast-Forward
Step 1
1.	For each pair of scenarios ($s^k$ and $s^u)$, the distance is computed by using the metric $c_T$. The generic element $C_{k,u}$ of matrix $C$ in step 1 is:
2.	$$C_{k,u}^{\left[1\right]}=c_T\left(s^k,s^u\right),k,u=1,\dots ,n_S$$ (7)
3.	The metric usually used is the ${\ell }_q$-Norm of ${ℝ}^T$, which can be defined as:
4.	$$c_T\left(s^k,s^u\right)={\left({\displaystyle \sum }_{\begin{array}{c}t=1\\ \end{array}}^T{\left|s_t^k-s_t^u\right|}^q\right)}^{\frac{1}{q}}$$ (8)
5.	Each scenario $s^u$ is associated with the weighted distance to any other scenario $s^k$, where the weights are the probabilities of occurrence ${\pi }_k$:
6.	$$z_u^{\left[1\right]}={\displaystyle \sum }_{\begin{array}{c}k=1\\ k\ne u\\ \end{array}}^{n_S}{\pi }_k{C}_{k,u}^{\left[1\right]},u=1,\dots ,n_S$$ (9)
7.	For example, the $z$ values for scenarios $s^1$ and $s^2$ are the following: $$\begin{gathered} z_1^{\left[1\right]}={\pi }_2{C}_{21}^{\left[1\right]}+{\pi }_3{C}_{\begin{array}{c}31\\ \end{array}}^{\left[1\right]}+{\pi }_4{C}_{41}^{\left[1\right]}+\dots \\ {z}_2^{\left[1\right]}={\pi }_1{C}_{12}^{\left[1\right]}+{\pi }_3{C}_{\begin{array}{c}32\\ \end{array}}^{\left[1\right]}+{\pi }_4{C}_{42}^{\left[1\right]}+\dots \end{gathered}$$
8.	Among the results, the index of the scenario with the minimum value of z is selected ($u_1$):
9.	$$u_1\in argmi{n}_{u\in \left\{1,\dots ,n_S\right\}}{z}_u^{\left[1\right]}$$ (10)
10.	Then, $s^{u_1}$ is preserved (operatively, $u_1$ is removed from the indexes of scenarios to delete in step 1, $J^{\left[1\right]}$):
11.	$$J^{\left[1\right]}=\left\{1,\dots ,n_S\right\}\backslash \left\{u_1\right\}$$ (11)
Step i
12.	Using the information from previous steps, the distance matrix is updated using Equation (12), new values of $z$ are computed using Equation (13), and a new scenario is selected to be preserved ($s^{u_i}$) using Equations (14) and (15):
13.	$$C_{ku}^{\left[i\right]}=\min \left\{C_{ku}^{\left[i-1\right]},C_{k{u}_{i-1}}^{\left[i-1\right]}\right\},k,u\in J^{\left[i-1\right]}$$ (12)
14.	$$z_u^{\left[i\right]}={\displaystyle \sum }_{k\in J^{\left[i-1\right]}\backslash \left\{u\right\}}{\pi }_k{C}_{ku}^{\left[i\right]},u\in J^{\left[i-1\right]}$$ (13)
15.	$$u_i\in argmi{n}_{u\in J^{\left[i-1\right]}}{z}_u^{\left[i\right]}$$ (14)
16.	$$J^{\left[i\right]}=J^{\left[i-1\right]}\backslash \left\{u_i\right\}$$ (15)
Step nP + 1
17.	In the final step, the list of scenarios to remove $J=J^{\left[n_P\right]}$ is completed. Each scenario to be removed will be linked to a preserved scenario that will “substitute” it. In fact, $j\left(i\right)$ is the index of the preserved scenario nearest to the removed scenario $s^i$:
18.	$$j\left(i\right)\in argmi{n}_{j\notin J}{c}_T\left(s^i,s^j\right),\forall i\in J$$ (16)
19.	The set of indexes of the removed scenarios that have $s^j$ as the nearest preserved scenario can be defined as follows:
20.	$$J\left(j\right)=\left\{i\in J:j=j\left(i\right)\right\}$$ (17)
21.	Using the optimal redistribution rule [40], the probability of the occurrence ${\pi }_j$ of the preserved scenario $s^j$ is computed:
22.	$${\pi }_j={\pi }_j+{\displaystyle \sum }_{\begin{array}{c}i\in J\left(j\right)\\ \end{array}}{\pi }_i$$ (18)
23.	The probabilities of the occurrence of the removed scenarios nearest to $s^j$ are added to the initial value of ${\pi }_j$.

3. Numerical Results

In this section, the numerical results of the scenario generation process carried out considering different values of $n_R$, metric, and p are presented.

Due to the randomness of sampling from uniform random variables in the Roulette Wheel method, two successive executions of the scenario generation process with the same parameters could return different scenarios unless fixing the seed for sampling from the uniform distribution. In order to compare the several generated scenarios with different parameters, the same seed was used for all trials related to the same month.

3.1. Outlier Removal

Figure 5 and Figure 6 show the results obtained by using $n_S=1000$, $n_P=10$, $n_R=7$, and metric = ${\ell }_2$-Norm and by varying p.

3.2. Number of Regions

Figure 7 and Figure 8 show the results obtained by considering $n_S=1000$, $n_P=10$, and metric = ${\ell }_2$-Norm without outlier removal and by varying $n_R$.

3.3. Metric

Figure 9 and Figure 10 show the results obtained by using $n_S=1000$, $n_P=10$, and $n_R=7$ without outlier removal and varying the metric used to compare two scenarios. For that, several ${\ell }_q$-Norms were tested (Equation (8)): ${\ell }_1$-Norm ($q=1$), ${\ell }_2$-Norm ($q=2$), ${\ell }_4$-Norm ($q=4)$, and ${\ell }_{\infty }$-Norm (Equation (19)).

$$c_T\left(s^k,s^u\right)=\underset{t}{\max }\left|s_t^k-s_t^u\right|$$

(19)

4. Discussion

In order to show the impact of the parameters p, $n_R$, and metric on the preserved scenarios (in terms of smoothness and variability in the same hour), a sensitivity analysis was performed.

For the estimation of the variability in the preserved scenarios, the average (for all hours excluding those in which irradiance is always zero) of the difference between the 97.5th and 2.5th percentiles of the preserved scenarios was considered. The trends of the average of the difference between the 97.5th and 2.5th percentiles for the different trials described in Section 3.1, Section 3.2 and Section 3.3 for July (summer) and January (winter) are shown in Figure 11 and Figure 12, respectively.

Regarding the outlier removal process, when reducing p, more samples for each hour are considered outliers (and hence removed). For the solar irradiance distributions with small hourly IQRs and many values outside IQRs, as in July, where the average IQR is about 111 W/m² (Figure 2), removing outliers and reducing the value of p lead to generated scenarios with small variability (Figure 5 and blue line in Figure 11). For the solar irradiance distributions where hourly IQRs are high and few observations are outside IQRs, as in January, where the average IQR is about 460 W/m² (Figure 3), the outlier removal process does not have, in general, a high impact (Figure 6 and the blue line in Figure 12).

Regarding the impact of the number of regions ($n_R$), the variability is increased for both seasons when $n_R$ is increased from 3 to 7 (Figure 7 and Figure 8 and orange lines in Figure 11 and Figure 12). With a value of $n_R$ that is higher than 7, the variability is quite stable. This can be explained by the fact that, even though the support of the beta distribution has been divided using more regions, more of them have a low probability of occurrence and, hence, do not impact the resulting variability.

Regarding the impact of several metrics for comparing two scenarios during the reduction process (Figure 9 and Figure 10 and green lines in Figure 11 and Figure 12), it is possible to observe that very different results are obtained from different metrics in the two considered seasons (e.g., ${\ell }_2$-Norm produces scenarios with high variability in summer and low variability in winter, ${\ell }_4$-Norm produces scenarios with lower variability than ${\ell }_2$-Norm in summer and higher variability than ${\ell }_2$-Norm in winter, etc.), and the variability obtained with ${\ell }_{\infty }$-Norm is close to that obtained with ${\ell }_2$-Norm.

The user of the method can select the best combination of the presented parameters in order to obtain the best trade-off between the variability among preserved scenarios and their plausibility. To verify the plausibility of the preserved scenarios, the user can plot them in the box-whisker plot obtained from the real solar irradiance data (Figure 1 and Figure 2) and verify that the generated values of hourly solar irradiance for the particular scenario do not deviate too much from the boxes and whiskers. (In the box-whisker plot, the box describes the range between $Q_1$ and $Q_3$, the upper whisker is on $Q_3+1.5\cdot IQR$, and the lower whisker is on $Q_1-1.5\cdot IQR$. The values outside the range defined by the whiskers can be considered outliers). In this case, this means that the generated values of hourly solar irradiance for that scenario can be considered inliers with respect to the observed data.

The scenarios obtained with $n_S=1000$, $n_P=10$, and $n_R=7$, ${\ell }_2$-Norm metric and without the outlier removal step are shown in Figure 13 for summer and Figure 14 for winter. Figure 15 and Figure 16 show the preserved scenarios plotted over the box-whisker plot of the observed solar irradiance. In these images, it is possible to see that the preserved scenarios are almost completely contained in the boxes (representing the IQR of the observed solar irradiance), and almost all values that are outside the boxes are confined to the variability range observed in the real data (the range between the two whiskers) and hence can be considered plausible.

5. Conclusions

In this work, a comprehensive tool to generate solar irradiance profiles is presented. The proposed approach is based on a scenario generation process aimed at generating 24 h solar irradiance scenarios using the historical data of solar irradiance for a specific location.

In the case study, the proposed method was applied to generate a set of daily solar irradiance scenarios for the months of January and July for the city of Turin (Italy). The Roulette Wheel mechanism was used to generate the initial set of scenarios, and the Fast-Forward method for the reduction process was applied to preserve the most representative scenarios and reduce the computational efforts associated with the potential stochastic operation optimization phase. The results demonstrate the flexibility of the method in generating scenarios for solar irradiance and in assessing their plausibility. These characteristics make the proposed approach an effective tool to be used for the stochastic operation optimization of DER.

Moreover, the results of the sensitivity analysis show the influence of the variation in the key parameters on the results in terms of increasing the variability and/or the smoothness of the generated scenarios, which could be very effective in estimating the behavior of the stochastic operation optimization of DER in the presence of more fluctuating but plausible solar irradiance patterns.

Given the generality of the proposed method, it can be easily adapted to model solar irradiance profiles for different locations and use cases, and hence, it can serve as a guide to users for the definition of scenarios with specific characteristics. Moreover, the proposed pipeline can be implemented as a Web Service queryable by users in order to generate solar irradiance scenarios with their probability of occurrence, which is fundamental for the stochastic optimization of DER.