Investigating the Effect of Uncertainty Characteristics of Renewable Energy Resources on Power System Flexibility

Featured Application: This study proposes a method for quantifying the effect of uncertainty characteristics of renewable energy resources on power system ﬂexibility, based on the ramping capability shortage probability index. The procedure involves scenario generation and evaluation of uncertainty characteristics. The proposed approach can provide system operators with the means to evaluate uncertainty characteristics, as well as possible strategies for a more ﬂexible power system. Abstract: This study investigates the effect of uncertainty characteristics of renewable energy resources on the ﬂexibility of a power system. The more renewable energy resources introduced, the greater the imbalance between load and generation. Securing the ﬂexibility of the system is becoming important to manage this situation. The degree of ﬂexibility cannot be independent of the uncertainty of the power system. However, most existing studies on ﬂexibility have not explicitly considered the effects of uncertainty characteristics. Therefore, this study proposes a method to quantitatively analyze the effect of uncertainty characteristics on power system ﬂexibility. Here, the uncertainties of the power system indicate the net load forecast error, which can be represented as a probability distribution. Of the characteristics of the net load forecast error, skewness and kurtosis were considered. The net load forecast error was modeled with a Pearson distribution, which has been widely used to generate the probability density function with skewness and kurtosis. Scenarios for the forecast net load, skewness, and kurtosis were generated, and their effects on ﬂexibility were evaluated. The simulation results for the scenarios based on a modiﬁed IEEE-RTS-96 revealed that skewness is more effective than kurtosis. The proposed method can help system operators to efﬁciently respond to changes in the uncertainty characteristics of renewable energy resources. a Pearson distribution. A case study for a modified IEEE-RTS-96 was conducted for the scenarios of FNL, skewness, and kurtosis. The simulation results also indicated that an abrupt change in the RSP can occur with the increase in the FNL. It was also observed that skewness with fixed kurtosis was more effective than kurtosis with fixed skewness. The proposed method is expected to help system operators to cope with changes in uncertainty characteristics. In future work, the application of this method to a real power system will be interesting.


Introduction
Since the adoption of the Paris Agreement in 2015, global CO 2 emissions have risen by approximately 4% [1]. In the future, energy transition and proactive efforts such as revising the nationally determined contributions in each country are needed to reduce emissions. The Korean government has been struggling with the enactment and reformation of policies and systems related to climate change and energy transition. The Energy Transition Roadmap and the Renewable Energy 3020 Implementation Plan were established in 2017 [2,3], and many efforts have been made since then to reduce coal power plants and expand new and renewable energy resources. In the Ninth Basic Electric Power Supply and Demand Plan established in 2020, the ratio of renewable energy generation to the total electricity generation in 2040 is targeted at 26.3% (for reference, its value in 2020 was 7.5%) [4]. To enhance the realizability of this plan and revitalize the related industrial ecosystems under these plans, the Fifth Renewable Energy Basic Plan was established in 2020 [5].
A i,t−∆t O i,t−∆t min(P max,i − P i,t−∆t , rr i ∆t) (1) Here, I means a set of generators. A i,t−∆t is a random variable representing the availability of generator i at time t−∆t. O i,t−∆t is a variable representing whether generator i is online at time t−∆t. min() is a function selecting a smaller value in the brackets. P max,i is a maximum output of generator i. P i,t−∆t is the output of generator i at time t−∆t. Rr i is the ramp rate of generator i. ∆t is the minimum interval between operating points. The SRC t represents the amount of RC the system can supply during time ∆t, considering the uncertainty. The larger the SRC, the greater the flexibility. However, the uncertainty of the system can reduce the RC availability, which directly relates to the parameter A i,t−∆t , which indicates whether the unit fails from t−∆t to t. The values of O i,t−∆t and P i,t−∆t are determined by the generation schedule.
The RC requirement (RCR) can be expressed as follows: where Here, NLFE t is a random variable representing net load forecast error at time t. FNL t is a forecast net load at time t. LFE t is a random variable representing load-forecast error at time t. VGFE t is a random variable representing variable generation forecast error at time t. FL t is a forecasted load at time t. FVG t is the forecast variable generation at time t. The RCR t indicates the amount of RC required to manage unexpected net load fluctuations and unit failures from t−∆t to t. The random variables NLFE t , LFE t , and VGFE t represent the forecast errors.

Ramping Capability Shortage Probability (RSP)
If the situation where SRC t is lower than RCR t continues, load shedding is required; this situation is called the RC shortage. The risk of an RC shortage at time t is defined as the RC shortage probability (RSP t ) and is represented as follows: Here, E t is a set of NLFE t and e is an element of E t . Prob(·) is a function calculating probability in the brackets. C t−∆t is a set of combinations of A i,t−∆t when O i,t−∆t is nonzero for all i. c is an element of C t−∆t . Prob c [·] means a probability of c if condition [·] is satisfied, 0 otherwise. The random variables representing the uncertainties are A t−∆t and NLFE t . C t−∆t and E t are the set of possible cases for A i,t-∆t and NLFE t , respectively. FNL t determines an online set of generating units, and NLFE t is expressed as the % value of FNL t ; thus, any variation in FNL t may affect not only O t−∆t and P i,t−∆t but also A i,t−∆t and NLFE t . The RSP t is calculated in the worst-case scenario, i.e., unexpected unit failures and net-load deviation occur just before the considered time t, and the imbalance cannot be resolved by a post-recovery action. This consideration is useful for reducing the computational burden of large-scale power systems. Meanwhile, the effect of complementarity between RERs is reflected in the NLFE t. If the aggregated output of RERs changes, FNL t varies, and VGFE t and NLFE t vary accordingly. For reference, C t−∆t is exemplified as follows: two generating units are online at times 1 and 2. The state of each unit at each time can be represented as either 1 (online) or 0 (failed). All possible cases at time 1 (i.e., C 1 ) are thus represented as (0,0), (0,1), (1,0), and (1,1). For reference, the probability of every element is calculated using a Markov-chain-based generator model; for further details, please refer to [27]. E t is exemplified as follows: FNL t at time t is 10 MW. The NLFE t is assumed to follow the skewed distribution shown in Figure 1. E t is thus (−1 MW), (0 MW), and (1 MW) by multiplying 10% and 10 MW, and the probability of each element is 0.25, 0.25, and 0.5, respectively. large-scale power systems. Meanwhile, the effect of complementarity between RERs is reflected in the NLFEt. If the aggregated output of RERs changes, FNLt varies, and VGFE and NLFEt vary accordingly. For reference, Ct−Δt is exemplified as follows: two generating units are online at times 1 and 2. The state of each unit at each time can be represented as either 1 (online) or 0 (failed). All possible cases at time 1 (i.e., C1) are thus represented as (0,0), (0,1), (1,0), and (1,1). For reference, the probability of every element is calculated using a Markov-chain-based generator model; for further details, please refer to [27]. Et is exemplified as follows: FNLt at time t is 10 MW. The NLFEt is assumed to follow the skewed distribution shown in Figure 1. Et is thus (−1 MW), (0 MW), and (1 MW) by multiplying 10% and 10 MW, and the probability of each element is 0.25, 0.25, and 0.5, respectively.

Modeling NLFE Characteristics Using a Pearson Distribution
The output of RERs is volatile and uncertain. The forecasted output involves a forecast error, which can be modeled as a random variable. Regardless of how well the output is predicted, the forecast error is inevitable, and the system operator cannot avoid it. Flexibility enhancement is thus necessary to deal with the forecast error, which indicates the risk of the system. The system operator should determine an adequate level of flexibility, which can be quantified using the RSP index. The forecast error of the renewable energy resource, i.e., VGFE, is reflected within the NLFE in the RSP calculation., i.e., NLFE is set as LFE minus VGFE. This RSP method was used to estimate the RC shortage, and the RC requirement was associated with the NLFE.
The normal distribution has been used to model NLFE in the RSP calculation [28,29]; however, it has a limitation in representing the distribution of NLFE characteristics such as skewness and kurtosis [30]. The Pearson distribution has been widely used to generate the PDF with skewness and kurtosis [31]. An advantage of the Pearson distribution over other distributions, such as the Weibull distribution, is the use of PDF moments to determine the distribution parameters. This distribution can be modeled as a function of skewness and kurtosis, and it allows the separate analysis of the effects of skewness and kurtosis. Thus, the Pearson distribution is applied to the NLFE model. It refers to a family of several probability distributions; it is made up of seven probability distributions, i.e., Type I to VII. If the values of the first four moments are known, i.e., mean, standard deviation, skewness, and kurtosis, the PDFs can be generated. Some of these types are listed in Table 1.

Modeling NLFE Characteristics Using a Pearson Distribution
The output of RERs is volatile and uncertain. The forecasted output involves a forecast error, which can be modeled as a random variable. Regardless of how well the output is predicted, the forecast error is inevitable, and the system operator cannot avoid it. Flexibility enhancement is thus necessary to deal with the forecast error, which indicates the risk of the system. The system operator should determine an adequate level of flexibility, which can be quantified using the RSP index. The forecast error of the renewable energy resource, i.e., VGFE, is reflected within the NLFE in the RSP calculation., i.e., NLFE is set as LFE minus VGFE. This RSP method was used to estimate the RC shortage, and the RC requirement was associated with the NLFE.
The normal distribution has been used to model NLFE in the RSP calculation [28,29]; however, it has a limitation in representing the distribution of NLFE characteristics such as skewness and kurtosis [30]. The Pearson distribution has been widely used to generate the PDF with skewness and kurtosis [31]. An advantage of the Pearson distribution over other distributions, such as the Weibull distribution, is the use of PDF moments to determine the distribution parameters. This distribution can be modeled as a function of skewness and kurtosis, and it allows the separate analysis of the effects of skewness and kurtosis. Thus, the Pearson distribution is applied to the NLFE model. It refers to a family of several probability distributions; it is made up of seven probability distributions, i.e., Type I to VII. If the values of the first four moments are known, i.e., mean, standard deviation, skewness, and kurtosis, the PDFs can be generated. Some of these types are listed in Table 1. Table 1. Main PDF types of the Pearson distribution according to the parameter γ.

PDF Type
Equations with Origin at the Mean Type Decision In the last column, the parameter γ is defined as follows: where P sk and P k are the skewness and kurtosis of the NLFE, respectively. The PDF type is determined according to P sk , P k , and γ. For reference, parameters such as y 0 , A 1 , A 2 , m 1 , m 2 , a, r, v, m, q 1 , and q 2 can be calculated using the following moments: mean, standard deviation, skewness (P s ), and kurtosis (P k ). For further details, please refer to [32]. Table 2 presents a simple example of the parameters of the Pearson distribution of NLFE, and the corresponding results are shown in Figure 1. The skewness and kurtosis are only varied, i.e., the mean and standard deviations are kept the same. The results for C11 to C17 and those of C21 to C27 correspond to Figure 2a,b, respectively. In Figure 2a, the PDF for the negative (positive) skewness is skewed to the left (right). In Figure 2b, the larger the kurtosis, the sharper the PDF at the center; however, the difference between each PDF according to the increase in kurtosis becomes increasingly smaller. Table 2. Numerical example of Pearson distribution of NLFE: (mean: 0, standard deviation: 1).

Case
Skewness Kurtosis PDF Type

Evaluation Procedure for the Effect of NLFE Characteristics on Flexibility
To define the effect-evaluation procedure, the following thing was first considered: the time interval when the largest RSP t occurs was chosen because it is the riskiest time in terms of flexibility. This assumption helps enhance analyzability while considering that the entire time interval may be time-consuming in analyzing the result. Moreover, this study does not focus on evaluating the flexibility related to the load-generation balance for the whole period but rather on analyzing the effect on flexibility according to the changes in the NLFE characteristics. If the entire period is considered, the effect of the NLFE characteristics is weakened during the evaluation of flexibility. In Figure 2a, the PDF for the negative (positive) skewness is skewed to the left (right). In Figure 2b, the larger the kurtosis, the sharper the PDF at the center; however, the difference between each PDF according to the increase in kurtosis becomes increasingly smaller.

Evaluation Procedure for the Effect of NLFE Characteristics on Flexibility
To define the effect-evaluation procedure, the following thing was first considered: the time interval when the largest RSPt occurs was chosen because it is the riskiest time in terms of flexibility. This assumption helps enhance analyzability while considering that the entire time interval may be time-consuming in analyzing the result. Moreover, this study does not focus on evaluating the flexibility related to the load-generation balance for the whole period but rather on analyzing the effect on flexibility according to the changes in the NLFE characteristics. If the entire period is considered, the effect of the NLFE characteristics is weakened during the evaluation of flexibility. In addition to the NLFE characteristics, the FNL is included as a parameter because it can affect the online state of generating units through the generation schedule and thus the flexibility. The FNL also changes the extent to which the NLFE characteristics affect the flexibility because the NLFE is expressed as the % value of the FNL. Thus, three parameters were considered: skewness, kurtosis, and FNL. The generation schedule and RSP can be varied using these parameters. Its influence on flexibility is the RSP result. The scenarios for NLFE characteristics and FNL are listed in Table 3. The impact of skewness and kurtosis on flexibility may vary with FNL because the generation schedule can be changed with the FNL. Thus, the FNL is arbitrarily divided into three classes-high, medium, and low. The basic case of the FNL was set to the medium class. The values of skewness and kurtosis are given as a particular range; in particular, the sign of skewness needs to be considered. Each parameter in the scenarios can be selected based on the given system. Empirical approaches can be used to determine the class and level of the FNL. To analyze the individual effect for every FNL, the following condition is considered: the skewness (kurtosis) is fixed, and the kurtosis (skewness) is changed. The influential factor can be determined through the flexibility evaluation for each scenario. The evaluation procedure is summarized in Figure 3.
the flexibility because the NLFE is expressed as the % value of the FNL. Thus, three p rameters were considered: skewness, kurtosis, and FNL. The generation schedule and RS can be varied using these parameters. Its influence on flexibility is the RSP result. Th scenarios for NLFE characteristics and FNL are listed in Table 3. The impact of skewness and kurtosis on flexibility may vary with FNL because th generation schedule can be changed with the FNL. Thus, the FNL is arbitrarily divide into three classes-high, medium, and low. The basic case of the FNL was set to the m dium class. The values of skewness and kurtosis are given as a particular range; in parti ular, the sign of skewness needs to be considered. Each parameter in the scenarios can b selected based on the given system. Empirical approaches can be used to determine th class and level of the FNL. To analyze the individual effect for every FNL, the followin condition is considered: the skewness (kurtosis) is fixed, and the kurtosis (skewness) changed. The influential factor can be determined through the flexibility evaluation fo each scenario. The evaluation procedure is summarized in Figure 3.

Case Study Basic Information
The main goal was to investigate the effects of the characteristics of uncertainty on flexibility. A test system is designed by adding the wind power plants of 1250 MW to the IEEE-RTS-96; an effective load-carrying capability of the plants is estimated at  Figure 4 shows the FNL profile; a peak load of 2,531 MW occurs at hour 19. The hourly generation schedule is determined using dynamic programming. All generating units are assumed in the spinning state; this condition is required to estimate the maximum value of the flexibility of the system. The RSP t was calculated using MATLAB programming (R2020a version) [34]. The simulation for the test system was performed on a PC with 64 GB RAM and a 3.49 GHz CPU. As shown in Figure 5, a seven-step approximation was applied to calculate the probability of the NLFE distribution.

Basic Information
The main goal was to investigate the effects of the characteristics of uncertainty on flexibility. A test system is designed by adding the wind power plants of 1250 MW to the IEEE-RTS-96; an effective load-carrying capability of the plants is estimated at 79.88 MW [33]. The test system is composed of 26 generating units of 3,105 MW. Figure 4 shows the FNL profile; a peak load of 2,531 MW occurs at hour 19. The hourly generation schedule is determined using dynamic programming. All generating units are assumed in the spinning state; this condition is required to estimate the maximum value of the flexibility of the system. The RSPt was calculated using MATLAB programming (R2020a version) [34]. The simulation for the test system was performed on a PC with 64 GB RAM and a 3.49 GHz CPU. As shown in Figure 5, a seven-step approximation was applied to calculate the probability of the NLFE distribution.  According to the solution procedure shown in Figure 3, the following steps were performed.
Step 1: The results of the RSPt calculation are shown in Figure 6. Hours 18-19 were determined as the most critical time interval because the highest value appears at hour 19.

Basic Information
The main goal was to investigate the effects of the characteristics of uncertainty on flexibility. A test system is designed by adding the wind power plants of 1250 MW to the IEEE-RTS-96; an effective load-carrying capability of the plants is estimated at 79.88 MW [33]. The test system is composed of 26 generating units of 3,105 MW. Figure 4 shows the FNL profile; a peak load of 2,531 MW occurs at hour 19. The hourly generation schedule is determined using dynamic programming. All generating units are assumed in the spinning state; this condition is required to estimate the maximum value of the flexibility of the system. The RSPt was calculated using MATLAB programming (R2020a version) [34]. The simulation for the test system was performed on a PC with 64 GB RAM and a 3.49 GHz CPU. As shown in Figure 5, a seven-step approximation was applied to calculate the probability of the NLFE distribution.  According to the solution procedure shown in Figure 3, the following steps were performed.
Step 1: The results of the RSPt calculation are shown in Figure 6. Hours 18-19 were determined as the most critical time interval because the highest value appears at hour 19. According to the solution procedure shown in Figure 3, the following steps were performed.
Step 1: The results of the RSP t calculation are shown in Figure 6. Hours 18-19 were determined as the most critical time interval because the highest value appears at hour 19.
Steps 2 to 4: Table 4 shows the generated NLFE scenario. The class and level of FNL were selected such that the high class (low class) was set at 101% and 102% (99% and 98%) of the FNL of the base case. A day-ahead NLFE distribution of the California Independent System Operator (CAISO) in the USA is applied as the NLFE of the base case [22]. The values of the mean, standard deviation, skewness, and kurtosis were −0.002, 0.026, 0.715, and 4.725, respectively, and the corresponding PDFs are represented in Figure 7. The standard deviation of the NLFE can be calculated as approximately 65.8 MW, i.e., multiplying 2,531 MW with 0.026 because the "normalized" NLFE (of the x-axis) is Steps 2 to 4: Table 4 shows the generated NLFE scenario. The class and level of FN were selected such that the high class (low class) was set at 101% and 102% (99% and 98% of the FNL of the base case. A day-ahead NLFE distribution of the California Independen System Operator (CAISO) in the USA is applied as the NLFE of the base case [22]. Th values of the mean, standard deviation, skewness, and kurtosis were −0.002, 0.026, 0.71 and 4.725, respectively, and the corresponding PDFs are represented in Figure 7. Th standard deviation of the NLFE can be calculated as approximately 65.8 MW, i.e., mult plying 2,531 MW with 0.026 because the "normalized" NLFE (of the x-axis) is 0.026. Th ranges of the skewness and kurtosis change from −0.115 to 1.515 and 0.725 to 28.725 wit increments of 0.2 and 4, respectively. The PDFs for the sampled values of skewness an kurtosis in S5 and S6 are shown in Figure 8.        Step 5: The hourly generation schedules for every scenario, i.e., S1 to S10, are determined based on the given FNL. Flexibility evaluation was later performed for every scenario. Figures 9 and 10 show the RSP t results according to the changes in skewness and kurtosis over a certain range, respectively. In both Figures, RSP t decreases (increases) according to the decrease (increase) in the FNL. This is because the available reserve capacity of the generating unit varies with the FNL.
In each Figure, a sudden change in RSP t occurs, i.e., the points of change lie between S5 and S7, and S6 and S8 in Figures 9 and 10, respectively. Sensitivity analysis was performed for both skewness-and kurtosis-varying scenarios; the points of change were 99.7% of the FNL of the base case in both scenarios. If the RSP t is evaluated for a wider range of FNL than that of the FNL of the scenario, more points of change may appear. The changes in the FNL can cause a sudden change in the RSP t because the committed set and operating conditions of the generating units change before and after the point of change. When there is no abrupt change in the RSP t , the changes in skewness and kurtosis were more influential than the FNL. It can be simply confirmed by comparing S7 and S9 that the RSP t in both scenarios increases with skewness; however, there is little difference between RSP t in these scenarios as the skewness changes. crease in skewness, whereas those in S1, S3, and S5 decrease. If the PDF is skewed to the left (i.e., the increase in skewness), it indicates that the load unexpectedly decreases. The RSP index captures the risk associated with increasing loads only; thus, decreasing loads can be considered better cases in terms of risk, and therefore, the RSPt (i.e., risk) decreases. This is because the probability of small normalized NLFEs is also reduced. This is confirmed by the graphs in Figure 8a.

Conclusions
This study investigated the impact of skewness and kurtosis of the NLFE on flexibility. The contribution is to propose a method for the effective analysis of NLFE characteristics on flexibility. The effect on flexibility was captured using a flexibility index called RSP. The skewness and kurtosis of the NLFE were represented with a Pearson distribution. A case study for a modified IEEE-RTS-96 was conducted for the scenarios of FNL, skewness, and kurtosis. The simulation results also indicated that an abrupt change in the RSP can occur with the increase in the FNL. It was also observed that skewness with fixed kurtosis was more effective than kurtosis with fixed skewness. The proposed method is expected to help system operators to cope with changes in uncertainty characteristics. In future work, the application of this method to a real power system will be interesting. S8 S10 Figure 10. RSP t results for kurtosis-varying scenarios, i.e., S2, S4, S6, S8, and S10.
The ratio of the skewness increment, i.e., 0.2 (the kurtosis increment, i.e., 4), to the skewness (kurtosis) of the base case is approximately 28% (85%). The ratio for skewness is smaller than that for kurtosis; however, the effect of skewness is greater than that of kurtosis. It can be compared with the relative RSP t ratio, defined as the RSP t divided by the ratio calculated above. The largest relative RSP t for the skewness-varying (kurtosis-varying) scenario is 4.19% (1.41%); for reference, the largest changes in RSP t in the skewness-varying and kurtosis-varying scenarios are between 4.725 and 8.725 and 1.315 and 1.515 in S5 and S6, respectively. The average value of the relative RSP t for the skewness-varying (kurtosis-varying) scenarios for the entire interval was 1.72% (0.41%). In the kurtosisvarying scenarios, the RSP t in S2, S4, S6, S8, and S10 decreases with an increase in kurtosis. In the skewness-varying scenarios, the RSP t in S7 and S9 increases with an increase in skewness, whereas those in S1, S3, and S5 decrease. If the PDF is skewed to the left (i.e., the increase in skewness), it indicates that the load unexpectedly decreases. The RSP index captures the risk associated with increasing loads only; thus, decreasing loads can be considered better cases in terms of risk, and therefore, the RSP t (i.e., risk) decreases. This is because the probability of small normalized NLFEs is also reduced. This is confirmed by the graphs in Figure 8a.

Conclusions
This study investigated the impact of skewness and kurtosis of the NLFE on flexibility. The contribution is to propose a method for the effective analysis of NLFE characteristics on flexibility. The effect on flexibility was captured using a flexibility index called RSP. The skewness and kurtosis of the NLFE were represented with a Pearson distribution. A case study for a modified IEEE-RTS-96 was conducted for the scenarios of FNL, skewness, and kurtosis. The simulation results also indicated that an abrupt change in the RSP can occur with the increase in the FNL. It was also observed that skewness with fixed kurtosis was more effective than kurtosis with fixed skewness. The proposed method is expected to help system operators to cope with changes in uncertainty characteristics. In future work, the application of this method to a real power system will be interesting.