Frequency Performance Distribution Index for Short-Term System Frequency Reliability Forecast Considering Renewable Energy Integration

: Risk in a power system’s ability to survive imminent disturbances without recourse to low operational cost and non-interruptive energy delivery remains the responsibility of every grid operator. Intermittencies in renewable energy and dynamic load variations inﬂuence the quality of power supply. The sudden changes a ﬀ ect the system frequency, compromising the reliability of the system grid; generation response to frequency regulation is momentous in such an incident. Slower response or smaller reserve capacity may cause a power shortage. This paper proposes a novel predictive scheme for a short-term operational reliability evaluation for system operations planning. The proposed method evaluates the operational reliability of system frequency whiles considering high renewable power penetration and energy storage system incorporation. Required energy generations, and other grid parameters, are modelled as stochastic inputs to the framework. We formulate a reliability index as a frequency distribution considering system frequency control dynamics and processes. The IEEE Reliability Test System (RTS) is used to prove the e ﬃ cacy of the proposed model.


Introduction
Power system security defines the probability of the system's operating point remaining within acceptable criteria, given the likelihood of changes in contingencies and the operational environments. The uncertainties in demand load and generation supply concern both operational reliability and the stability of the system. Both system frequency and voltage provide the base for reliable assessments of the system grid. Hence, adhering to frequency and voltage regulations is essential for the operational security of the network in an unforeseen contingency. Frequency indicates the system's active power adequacy in system operation. Real-time operation of a power system with high renewable penetration, variable demand loads and renewable power, and frequent generator failures could cause considerable variations in system frequency. Load frequency control (LFC) techniques [1][2][3] are used to address such problems. Over the years, reliability standards such as control performance standards (CPS) [4] and European Network of Transmission System Operators for Electricity (ENTSO-E) [5] and assessment methods, such as loss of load probability (LOLP) [6] and loss of load expectation (LOLE) [7] have been effective for long-term evaluation and reliability performance assessments of the system grid. Until now, many long-term reliability assessment methods have been effective for system planning. The reliability assessment base such Table 1. Regulation target of frequency [31]. The proposed index could facilitate optimal selection of spinning reserve (SR) capacity during system planning. Optimal SR capacity selection is momentous for the operational reliability evaluation and economic operation of the network. The proposed method addresses the problems mentioned Energies 2020, 13, 2945 3 of 17 above with forecasted stochastic time-series demand load and renewable energy generation. Practically, it is difficult to forecast the operational reliability of the system for any given specific time. System frequency could be a perfect reliability criterion only if we could predict its magnitude and rate of change inexactitude in real-time. Our motivation, in this regard, is to use the concept of stochastic grid modelling to predict the statistical behaviour of the system frequency. To this end, we propose a method to predict the system frequency statistics called the frequency reliability distribution function (FRDF) as shown in Figure 1, whose output describes the frequency distribution of the system frequency deviation of the grid over a short period (i.e., 5-10 min). In respect of this, the following outline our contributions to this course:

Area Frequency (Hz) Maximum Frequency Deviation
To develop a short-term stochastic operational reliability index considering system dynamic and frequency regulation processes.

2.
To develop a one-second interpolated prediction load model from a low-sampled load data. 3.
To develop a hybrid-impulse response (HIR) model to shorten the computational time for real-time operational reliability assessment.

4.
To develop selection criteria for optimal SR capacity with a battery energy storage system (BESS) for system economic operation.
The remaining of the paper is organized as follows: Section 2 discusses the design considerations that we implemented in formulating this new reliability criterion. In Section 3, we construct a mathematical model with the physical constraints of the system. Finally, several case studies are performed, proving the efficacy of the suggested index in Section 4.
Energies 2020, 13, x FOR PEER REVIEW 3 of 17 specific time. System frequency could be a perfect reliability criterion only if we could predict its magnitude and rate of change inexactitude in real-time. Our motivation, in this regard, is to use the concept of stochastic grid modelling to predict the statistical behaviour of the system frequency. To this end, we propose a method to predict the system frequency statistics called the frequency reliability distribution function (FRDF) as shown in Figure 1, whose output describes the frequency distribution of the system frequency deviation of the grid over a short period (i.e., 5-10 min). In respect of this, the following outline our contributions to this course: 1. To develop a short-term stochastic operational reliability index considering system dynamic and frequency regulation processes. 2. To develop a one-second interpolated prediction load model from a low-sampled load data. 3. To develop a hybrid-impulse response (HIR) model to shorten the computational time for real-time operational reliability assessment. 4. To develop selection criteria for optimal SR capacity with a battery energy storage system (BESS) for system economic operation.
The remaining of the paper is organized as follows: Section 2 discusses the design considerations that we implemented in formulating this new reliability criterion. In Section 3, we construct a mathematical model with the physical constraints of the system. Finally, several case studies are performed, proving the efficacy of the suggested index in Section 4.

Model Structure
Our primary goal was to forecast a statistical representation of system frequency for a specific period, and utilize it for reliability criteria. The output of FRDF was frequency histogram elements (FHE). FHE, a frequency metric, is defined as the probabilistic distribution of sampled system frequency fluctuations existing under specified reliability criteria, for a specified period. Frequency performance distribution index was estimated with FHE. With dynamic stochastic time-series input parameters such as load demand, intermittent renewable generation, forced outage rate of generators, and battery energy storage systems (BESS), FRDF estimated system frequency for performance evaluation. A Monte Carlo simulation repeatedly executed the model to generate timeseries system frequency samples. This technique leverages the risk impact and uncertainty in the

Model Structure
Our primary goal was to forecast a statistical representation of system frequency for a specific period, and utilize it for reliability criteria. The output of FRDF was frequency histogram elements (FHE). FHE, a frequency metric, is defined as the probabilistic distribution of sampled system frequency fluctuations existing under specified reliability criteria, for a specified period. Frequency performance distribution index was estimated with FHE. With dynamic stochastic time-series input parameters such as load demand, intermittent renewable generation, forced outage rate of generators, and battery energy storage systems (BESS), FRDF estimated system frequency for performance evaluation. A Monte Carlo simulation repeatedly executed the model to generate time-series system frequency samples. This technique leverages the risk impact and uncertainty in the prediction models. The time series outcome for each trial in the simulation was accumulated and statistically defined to build FHE.
FHE, as shown in Figure 2, describes the degree of system frequency deviation from its nominal value for a given time under set criteria. Therefore, statistical summaries such as existing probability within a specified deviation range and frequency dispersion could serve as a reliability criterion.
Energies 2020, 13, x FOR PEER REVIEW 4 of 17 prediction models. The time series outcome for each trial in the simulation was accumulated and statistically defined to build FHE FHE, as shown in Figure 2, describes the degree of system frequency deviation from its nominal value for a given time under set criteria. Therefore, statistical summaries such as existing probability within a specified deviation range and frequency dispersion could serve as a reliability criterion.

Stochastic Random Time-Series Input Variables
Each input of FRDF was stochastically generated based on deterministic time-series input variables. However, inputs such as load demand, and renewable energy sources data as recorded by SCADA/EMS systems, were in minutes and higher sampling intervals. Such sampling rates are not adequate for short-term transient analysis, as FRDF performs frequency transient state analysis, with high sampling rates (i.e., 1-second) sampled data. So, for such a simulation, we developed a onesecond interpolation load predictive fluctuation model from a 15-min sampled data. The predictive fluctuation model formulation, and its usage in generating random time-series data for sequential input variables, are discussed in the subsequent section.

Computational Time Limitation
The usage of high-sampled data in the transient-state analysis contributed to the practical evaluation assessment of the system reliability but was computationally expensive. Practically, such analysis involves a complete inclusion of all grid components. However, such additions lead to colossal time consumption and are computationally expensive, even for single scenario simulation. So, we proposed a reduced simulation scheme (RSS) to expedite the computational simulation time, which would not relent on the evaluation performance. The reduced scheme lumps up grid components type, such as load-following generators types, renewable energies, and base loads, as one for each component type. The RSS showed significant computational time reduction when compared with the original grid model. Including the RSS model, nonetheless, for a repeated simulation with much-varying grid components, the computational time was still not desirable. Hence, we introduced an additional analytic method called a hybrid-impulse response (HIR) model. HIR enhances the computational time dramatically compared to the pure numerical simulation model while maintaining the accuracy of the evaluation assessment in a reasonable range. The subsequent section discusses the formulations of both the RSS and HIR models.

Stochastic Random Time-Series Input Variables
Each input of FRDF was stochastically generated based on deterministic time-series input variables. However, inputs such as load demand, and renewable energy sources data as recorded by SCADA/EMS systems, were in minutes and higher sampling intervals. Such sampling rates are not adequate for short-term transient analysis, as FRDF performs frequency transient state analysis, with high sampling rates (i.e., 1-s) sampled data. So, for such a simulation, we developed a one-second interpolation load predictive fluctuation model from a 15-min sampled data. The predictive fluctuation model formulation, and its usage in generating random time-series data for sequential input variables, are discussed in the subsequent section.

Computational Time Limitation
The usage of high-sampled data in the transient-state analysis contributed to the practical evaluation assessment of the system reliability but was computationally expensive. Practically, such analysis involves a complete inclusion of all grid components. However, such additions lead to colossal time consumption and are computationally expensive, even for single scenario simulation. So, we proposed a reduced simulation scheme (RSS) to expedite the computational simulation time, which would not relent on the evaluation performance. The reduced scheme lumps up grid components type, such as load-following generators types, renewable energies, and base loads, as one for each component type. The RSS showed significant computational time reduction when compared with the original grid model. Including the RSS model, nonetheless, for a repeated simulation with much-varying grid components, the computational time was still not desirable. Hence, we introduced an additional analytic method called a hybrid-impulse response (HIR) model. HIR enhances the computational time dramatically compared to the pure numerical simulation model while maintaining the accuracy of the evaluation assessment in a reasonable range. The subsequent section discusses the formulations of both the RSS and HIR models.

Predictive Fluctuation Model
The predictive fluctuation model devised in this work comprises two stages: load fluctuation estimation and low-sampling load interpolation, as shown in Figure 3. The load fluctuation model defines a relationship between a low-sampled predictive load period and its load fluctuation distribution for dynamic load fluctuation estimation. To derive such a relationship, we simulated the load fluctuation based on an actual investigation of a 15-min data sample and 1-s sampled data.

Predictive Fluctuation Model
The predictive fluctuation model devised in this work comprises two stages: load fluctuation estimation and low-sampling load interpolation, as shown in Figure 3. The load fluctuation model defines a relationship between a low-sampled predictive load period and its load fluctuation distribution for dynamic load fluctuation estimation. To derive such a relationship, we simulated the load fluctuation based on an actual investigation of a 15-minutes data sample and 1-second sampled data. From the sample data, we obtained deviation data by measuring the difference between adjacent sample points for each sample data. Outliers and missing data were respectively compensated for and removed by a preprocessing algorithm. Following this, a normalized histogram was fitted to the sample deviation data, as shown in Figure 4, to define a load fluctuation distribution. For the lowsample distribution to be adapted for transient state analysis, which requires high-sampling data, the correlation between a high-sampled and low-sampled data should be evident. From the simulation results, for a 3-hour sampling period, it is evident that there was a positive correlation between the distributions mentioned above with a correlation coefficient of 0.67 at a 95% confidence interval. Hence, low-sampled low fluctuation could suffice for transient state load fluctuation estimation. Algorithm 1 details the procedure for the modelling process.  From the sample data, we obtained deviation data by measuring the difference between adjacent sample points for each sample data. Outliers and missing data were respectively compensated for and removed by a preprocessing algorithm. Following this, a normalized histogram was fitted to the sample deviation data, as shown in Figure 4, to define a load fluctuation distribution. For the low-sample distribution to be adapted for transient state analysis, which requires high-sampling data, the correlation between a high-sampled and low-sampled data should be evident. From the simulation results, for a 3-h sampling period, it is evident that there was a positive correlation between the distributions mentioned above with a correlation coefficient of 0.67 at a 95% confidence interval. Hence, low-sampled low fluctuation could suffice for transient state load fluctuation estimation. Algorithm 1 details the procedure for the modelling process. Choose a load profile, L1 from accumulated interpolated simulated historical data.

5:
∆ n,m is normalized as δ n,m , where δ n,m = ∆ n,m+1 ∆ n,m . 6: Fit a histogram to δ n,m . From the histogram, we estimate the error distribution. 7: Repeat the process for all k. 8: Fit a model to estimate k with the minimum estimation error.
Following the load-fluctuation model was the interpolation model definition. With an initial forecasted value, the interpolation model estimated interpolated pseudo-1-s values for a given sampling period based on a one-dimensional random walk algorithm. The interpolated values were estimated based on the load fluctuation and error distribution function for the given forecast period. Successive  With a known dataset, the efficacy of the proposed model was verified. An actual set sample of demand loads: 1-second sampled data and 15-minutes sampled data were acquired from the Korea Power Exchange company [33] for the evaluation of the proposed model. The results of the assessment, as depicted in Figure 5, showed high similarities in the patterns of the original 15-minutes sampled demand and the corresponding 1-second sampled demand load data. The statistics of the evaluation results proved the efficacy of the proposed model. The error margin associated with the predictive model was minimal with a goodness of fit, 0.041 root mean square error (RMSE). The predictive model was adapted in estimating the fluctuation and interpolation for each stochastic input variables.
= + ~ , In the case of the demand power model, time-series data was randomly generated stochastically based on the aforementioned predictive model. With a 15-minutes sampled demand value and the assumption that the fluctuation in the demand data happening in every second was independent of each other, the predictive model was utilized to forecast a demand load profile for the transient state analysis of the frequency reliability. Similarly, wind power dynamic data input was generated with the proposed predictive fluctuation model. A time-series wind power profile was created based on forecasted wind speed, the rated power of the wind turbines, and the size of the wind farm. Based on a known (i.e., in this case, 490 kW) wind power generator data, and wind power generation curve, a corresponding relationship between the wind speed and the power was defined. With a known dataset, the efficacy of the proposed model was verified. An actual set sample of demand loads: 1-s sampled data and 15-min sampled data were acquired from the Korea Power Exchange company [33] for the evaluation of the proposed model. The results of the assessment, as depicted in Figure 5, showed high similarities in the patterns of the original 15-min sampled demand and the corresponding 1-s sampled demand load data. The statistics of the evaluation results proved the efficacy of the proposed model. The error margin associated with the predictive model was minimal with a goodness of fit, 0.041 root mean square error (RMSE). The predictive model was adapted in estimating the fluctuation and interpolation for each stochastic input variables.
Energies 2020, 13, x FOR PEER REVIEW 6 of 17  With a known dataset, the efficacy of the proposed model was verified. An actual set sample of demand loads: 1-second sampled data and 15-minutes sampled data were acquired from the Korea Power Exchange company [33] for the evaluation of the proposed model. The results of the assessment, as depicted in Figure 5, showed high similarities in the patterns of the original 15-minutes sampled demand and the corresponding 1-second sampled demand load data. The statistics of the evaluation results proved the efficacy of the proposed model. The error margin associated with the predictive model was minimal with a goodness of fit, 0.041 root mean square error (RMSE). The predictive model was adapted in estimating the fluctuation and interpolation for each stochastic input variables. In the case of the demand power model, time-series data was randomly generated stochastically based on the aforementioned predictive model. With a 15-min sampled demand value and the assumption that the fluctuation in the demand data happening in every second was independent of each other, the predictive model was utilized to forecast a demand load profile for the transient state analysis of the frequency reliability. Similarly, wind power dynamic data input was generated with the proposed predictive fluctuation model. A time-series wind power profile was created based on forecasted wind speed, the rated power of the wind turbines, and the size of the wind farm. Based on a known (i.e., in this case, 490 kW) wind power generator data, and wind power generation curve, a corresponding relationship between the wind speed and the power was defined.

Reduced Simulation
A simplified simulation reduction scheme was devised to minimize computational simulation time. The model was made up of base generators, load-following generators, renewable energies, and battery energy storage system (BESS). With this model, a set of generators, demand, wind power, and BESS were connected, as shown in Figure 6. The reduced scheme lumped up grid components type, such as load-following generators types, renewable energies, and base loads, with one for each component type to form a virtual grid to hasten the computational time. However, in a real-world scenario, such a phenomenon could be discretely applied considering the most common parameters among the committed set of generators. The model has stochastically modelled inputs: demand load, generators, and wind power, along with the control schemes of each resource type. For the BESS control, a high-pass filter was employed to consider the energy-limiting feature. Variable factors such as rate of spinning reserve, BESS capacity, the number of wind turbines were also defined in the model.

Reduced Simulation
A simplified simulation reduction scheme was devised to minimize computational simulation time. The model was made up of base generators, load-following generators, renewable energies, and battery energy storage system (BESS). With this model, a set of generators, demand, wind power, and BESS were connected, as shown in Figure 6. The reduced scheme lumped up grid components type, The RSS framework adopted in this work mimics the real-world scenario. In evaluating the reliability of the power system, the forced outage rate is applied to the generators. For each simulation trial, generators are selected based on the MTTF and MTTR parameters specified in Equation (2) from the total number of generators in a generating unit.
The selection of the generators was done by the stochastic estimation of the forced outage rates of the generators. The forced outage rate of a generator in a given period indicated the operational status and the maximum amount of power supplied. This phenomenon is wildly adapted in long-term reliability evaluation such as LOLP. However, concerning our proposed model, as stated earlier, our prime goal was to evaluate the reliability of the system for a short period. Hence, our interest lies in the probability that the grid balance was achieved with the total demand power and the entire supplied generation. Therefore, for each given scenario, we estimated the failure probability of the already determined generators in the units. The failure probability in this analysis was defined as the reciprocal of the MTTF, which is expressed in (3).
Six (6) different types of generator units with different characteristics were used in the simulation. The oil/steam generators U100 and U197 were utilized as load-following generators to make up the power imbalance in the system. Other generators units such as U400 (i.e., nuclear) were used as a base generator of which the output power was maintained continuously. Due to the adjustment of the spinning reserve, the number of generators and their operating rate were configured deterministically and deployed in different scenarios. The spinning reserves used in the various simulations was varied between 1-10%. Based on Equation (3), the failure probability of each generator was estimated.
In the course of the statistical estimation of the frequency deviation, the frequency control scheme and relative performance parameters, such as the droop-rate of the generators, were considered. The area control error was distributed among the allowable generators and the BESS. These errors were filtered by two filter components: low-pass filter and high-pass filter based on the fluctuation rate. Higher frequency fluctuations were filtered by the high-pass filter, whereas the low-pass filter filtered lower frequency fluctuations. The BESS control model dealt with the higher frequency fluctuations, while the remainder of the load fluctuations went to the load-following generators according to their droop-rates. The BESS model connected the energy storage batteries to the power grid. A cutoff frequency of 1.66 × 10 −3 Hz was defined for every 10 min cycle. The time constant was considered as a mechanical delay. The incorporation of BESS for LFC was as a result of its high response rate as compared with conventional generators. As a result, BESS could discharge energy faster to respond to higher load variations. BESS was combined with the spinning reserve for each simulation scenario. The effect of BESS on load variations was observed by varying the BESS capacity between 0-150 MW on an interval of 50 MW. A 490-kW rated capacity wind turbine generators (WTGs) with two different wind speeds: 5 m/s and 10.1 m/s were deployed in the simulation. The penetration levels of the wind turbines were varied at 0%, 30%, 60%, and 100% of the committed wind generators. The purpose of these variations was to measure the effect on the system frequency with increasing intermittencies in wind power generation.

Reliability Evaluation Method
For short-term reliability evaluation, traditionally, grid operators analyze the behaviour of the frequency by using a defined criterion to maintain the reliability and stability of the system operation. This standard has become a convention for reliability evaluation. In this work, the structure of the model, as shown in Figure 7, defined our reliability evaluation processes.
Where and are defined as the mean factor of safety and standard deviation of the factor of safety, respectively.
can be express in terms of the rate of system frequency change and the total number of critical elements (N).
= R Similarly, can be expressed in terms of dispersion of system frequency and N.
= − 1 (8) With this process, the frequency reliability was estimated based on the set reliability criteria. The output frequencies from FRDF formed a frequency distribution using a histogram. From this, frequency performance index as a reliability index was estimated to determine the reliability of the system. The system reliability index, β, can be defined in terms of the rate of system frequency change ( ) and dispersion of system frequency change ( ), as shown in Equations (4) and (8)

Hybrid Simulation Scheme for Reduced Simulation Time
The simulation configuration of MCS takes a lot of computational time and resources, even with the reduced simulation scheme. So, in addition to RSS, we developed a hybridized approach of numerical and analytical methods to reduce the computational time while maintaining the accuracy of the results. To build the analytical model to hasten the computational time, we first applied the concept of the system impulse response for the linear part in the simulation model. For a linear and time-invariant (LTI) system, an input signal can be decomposed into smaller (i.e., 1 MW) simple additive components called impulses, which are passed through a linear system, and the resulting output components called impulse response, are synthesized. In effect, the impulse decomposition provides a way to analyze signals one sample at a time. The signal resulting from this divide-andconquer procedure was identical to that obtained by directly passing the original signal through the system. Hence, with the system's impulse response, we could estimate what the output would be for any given input signal by convolving the input signal with the system impulse response. Convolution is important because it relates the input signal, the output signal, and the impulse response of the system. The grid control model can be expressed as a linear system, except for the ramp rate and limiter output, as shown in Figure 8. The non-linearity of the grid model inhibits the direct application of impulse response to enhance the computation time. The quantitative index representation of β can be defined, as shown in (15).
where µ FS and δ FS are defined as the mean factor of safety and standard deviation of the factor of safety, respectively. µ FS can be express in terms of the rate of system frequency change and the total number of critical elements (N).
Similarly, δ FS can be expressed in terms of dispersion of system frequency and N.
With this process, the frequency reliability was estimated based on the set reliability criteria. The output frequencies from FRDF formed a frequency distribution using a histogram. From this, frequency performance index as a reliability index was estimated to determine the reliability of the system. The system reliability index, β, can be defined in terms of the rate of system frequency change (R f ) and dispersion of system frequency change (D f ), as shown in Equations (4) and (8), respectively, where ∆ f is the frequency change under specified criteria (c 1 and c 2 ).

Hybrid Simulation Scheme for Reduced Simulation Time
The simulation configuration of MCS takes a lot of computational time and resources, even with the reduced simulation scheme. So, in addition to RSS, we developed a hybridized approach of numerical and analytical methods to reduce the computational time while maintaining the accuracy of the results.
To build the analytical model to hasten the computational time, we first applied the concept of the system impulse response for the linear part in the simulation model. For a linear and time-invariant (LTI) system, an input signal can be decomposed into smaller (i.e., 1 MW) simple additive components called impulses, which are passed through a linear system, and the resulting output components called impulse response, are synthesized. In effect, the impulse decomposition provides a way to analyze signals one sample at a time. The signal resulting from this divide-and-conquer procedure was identical to that obtained by directly passing the original signal through the system. Hence, with the system's impulse response, we could estimate what the output would be for any given input signal by convolving the input signal with the system impulse response. Convolution is important because it relates the input signal, the output signal, and the impulse response of the system. The grid control model can be expressed as a linear system, except for the ramp rate and limiter output, as shown in Figure 8. The non-linearity of the grid model inhibits the direct application of impulse response to enhance the computation time. Hence, we developed a composite model that deals with the linear part with impulse response before applying the non-linear part of the system model. Such a composite model was termed the hybrid-impulse response model (HIR). HIR, for each scenario, with the definition of the relevant grid parameters, the impulse response of the system was taken for the linear part of the grid model. Every generator output was acquired before the limiter by convolving with the input signal, which constituted the summation of the demand power, wind power output fluctuations, and 1 MW impulse response.  Hence, we developed a composite model that deals with the linear part with impulse response before applying the non-linear part of the system model. Such a composite model was termed the hybrid-impulse response model (HIR). HIR, for each scenario, with the definition of the relevant grid parameters, the impulse response of the system was taken for the linear part of the grid model. Every generator output was acquired before the limiter by convolving with the input signal, which constituted the summation of the demand power, wind power output fluctuations, and 1 MW impulse response.
Then, the output of the convolution was numerically manipulated to apply the non-linear characteristic of the system, say the limitations of ramp rate and saturation. The frequency impulse response was then convolved with the grid input parameters again to yield the final output of the system with the non-linear characteristic. Figure 9 depicts the simulation results of the hybrid-impulse model and the original simulation. From the results, it is clear that HIR does not relent on the performance accuracy even for a shorter computational time. The figure shows an instance of the simulation results for the control frequency simulation under a scenario which adopted a 7% spinning  Hence, we developed a composite model that deals with the linear part with impulse response before applying the non-linear part of the system model. Such a composite model was termed the hybrid-impulse response model (HIR). HIR, for each scenario, with the definition of the relevant grid parameters, the impulse response of the system was taken for the linear part of the grid model. Every generator output was acquired before the limiter by convolving with the input signal, which constituted the summation of the demand power, wind power output fluctuations, and 1 MW impulse response. Then, the output of the convolution was numerically manipulated to apply the non-linear characteristic of the system, say the limitations of ramp rate and saturation. The frequency impulse response was then convolved with the grid input parameters again to yield the final output of the system with the non-linear characteristic. Figure 9 depicts the simulation results of the hybridimpulse model and the original simulation. From the results, it is clear that HIR does not relent on the performance accuracy even for a shorter computational time. The figure shows an instance of the simulation results for the control frequency simulation under a scenario which adopted a 7% spinning reserve with 100 MW BESS and 490 kW wind power. From this result, it was convincing to note that the hybrid-impulse model can be used in place of the regular simulation model and still maintain the reliability of the evaluation quality.

Case Study
In this section, using the defined reliability index, we conducted several case studies to consider the allowable penetration level of wind power and spinning reserve, as well as the impact of BESS on the grid.

RTS: Reliability Test System
A reference test system was required to provide a basis for comparison of results from different approaches. So, we adapted the IEEE Reliability Test System (RTS) (1976) [34]. We built a specific model using the IEEE RTS model to test the efficacy of the proposed model. The RTS model defined set parameters for demand load, generation system, and transmission network. Our objective of adopting this model was to comply with the benchmark and establish an abstract system sufficiently broad to provide a basis for reporting on analysis methods for reliability methods. Admittedly, it was not practical to specify all the parameters needed for every application. Our goal was to establish a core system that can be augmented with additional or modified parameters required for a given use. The miniaturized adoption of this system consisted primarily of load model data, generation unit reliability data, generation unit operating cost data, renewable energy capacities, and energy data, forced outage data, and generation mix data.

Case 1: Effect of Rate of Spinning Reserve and BESS
As stated earlier, LFC is usually done with a spinning reserve. The spinning reserve is the amount of unused capacity in online energy assets, which can compensate for power shortages or frequency drops within a given period. The size and distribution of the spinning reserve has a significant impact on system frequency response. In this case study, under a given set of spinning reserve capacity, wind speed, and BESS capacity, the frequency performance index was estimated. The simulation was repeated for different scenarios. In each scenario, based on the set frequency criteria, FRDF estimated the frequency distribution. The frequency performance index was subsequently estimated based on the acquired frequency distribution. Figures 10 and 11 showed the performance index of the rate of frequency existing under ±1 increase with increasing spinning reserve capacity and BESS. From the figures, without the support of BESS, much spinning reserve capacity was required for higher penetration of wind speed. This requirement was as a result of the higher intermittency caused by the wind generators. BESS with large reserve capacity can support much renewable power intermittency; hence, an increase in the frequency rate, making the system more stable.  Traditionally, the spinning reserve is a capacity reserve of synchronous generators. If a massive operating generator trips the power network, the remaining less massive generators should increase their output to compensate for the recovery of the power shortage. Thus, a significant number of the committed generators are required to operate below their rated capacity to accommodate such unexpected contingencies should it happen. The more reserve capacity is needed, the further the operating points of the committed generators are from the rated capacity, which leads to an uneconomical, harmful, and inefficient operation of the grid. Luckily, BESS can be used to accommodate some of the spinning reserve requirement for the grid generators. From Figure 10a-d, the rate of frequency increase as the BESS capacity increases at a fixed spinning reserve. From this, it can be deduced that the support for the change in frequency resulting from supply-demand imbalance is shared between BESS and generators spinning reserves capacities based on their droop rate.
Hence, the higher the BESS capacity, the more it can support the spinning reserve. A similar occurrence can be seen in Figure 11a  Hence, the higher the BESS capacity, the more it can support the spinning reserve. A similar occurrence can be seen in Figure 11a-d. The figures show more decline in the rate of the system frequency compared with 5 m/s wind turbine simulations. Higher wind speed could lead to more fluctuations in the power generation. For 100 wind turbines, the variations in the frequency rate were much broader with 10.0 m/s wind speed power generation than it was in 5 m/s.

Case 2: Effect of Renewable Energy Sources
Recently, due to the high penetration of renewable energy in the generation mix, system reliability tends to be lower in current power systems than in traditional power systems. In this case study, for a given set of spinning reserve capacity, BESS capacity, and wind speed, the frequency performance index was estimated. The simulation was repeated for different scenarios. In each scenario, based on the set frequency criteria, FRDF estimated the frequency distribution. The frequency performance index: frequency dispersion was subsequently estimated based on the acquired frequency distribution. Figures 12 and 13 shows the frequency performance index of the frequency dispersion with an increasing number of wind turbines. Frequency dispersion describes the extent to which the system frequency disperses within the stability range. From the figures, frequency dispersion happened with increasing wind speed without BESS support. This occurrence happened because more wind turbines yielded more power intermittency on the grid, which resulted in a supply-demand imbalance. Energy generation is heavily reliant on weather conditions which change regularly. As a result, power shortages may happen more often on the grid. Hence, the frequency deviations. However, the frequency dispersion is retarded with increasing BESS capacity, because unlike load-following generators, BESS has a faster response time. The disturbance progression could be curtailed with BESS support. The more wind turbines, the worse the rate of frequency existing under the reliability criteria. The frequency performance index illustrates how the index estimates the influence of the penetration of renewable energy sources quantitativelycomparing Figure 10 and Figure 11, vast fluctuations of the output effect of the wind turbines on the index, resulted in a lousy convergence percentage of frequency rate existing under the criteria.

Case 2: Effect of Renewable Energy Sources
Recently, due to the high penetration of renewable energy in the generation mix, system reliability tends to be lower in current power systems than in traditional power systems. In this case study, for a given set of spinning reserve capacity, BESS capacity, and wind speed, the frequency performance index was estimated. The simulation was repeated for different scenarios. In each scenario, based on the set frequency criteria, FRDF estimated the frequency distribution. The frequency performance index: frequency dispersion was subsequently estimated based on the acquired frequency distribution. Figures 12 and 13 shows the frequency performance index of the frequency dispersion with an increasing number of wind turbines. Frequency dispersion describes the extent to which the system frequency disperses within the stability range. From the figures, frequency dispersion happened with increasing wind speed without BESS support. This occurrence happened because more wind turbines yielded more power intermittency on the grid, which resulted in a supply-demand imbalance. Energy generation is heavily reliant on weather conditions which change regularly. As a result, power shortages may happen more often on the grid. Hence, the frequency deviations. However, the frequency dispersion is retarded with increasing BESS capacity, because unlike load-following generators, BESS has a faster response time. The disturbance progression could be curtailed with BESS support. The more wind turbines, the worse the rate of frequency existing under the reliability criteria. The frequency performance index illustrates how the index estimates the influence of the penetration of renewable energy sources quantitatively-comparing Figures 10 and 11, vast fluctuations of the output effect of the wind turbines on the index, resulted in a lousy convergence percentage of frequency rate existing under the criteria. because unlike load-following generators, BESS has a faster response time. The disturbance progression could be curtailed with BESS support. The more wind turbines, the worse the rate of frequency existing under the reliability criteria. The frequency performance index illustrates how the index estimates the influence of the penetration of renewable energy sources quantitativelycomparing Figure 10 and Figure 11, vast fluctuations of the output effect of the wind turbines on the index, resulted in a lousy convergence percentage of frequency rate existing under the criteria.

Case 3: Reliability Assessments with Optimal SR
The traditional criterion for setting the minimum amount of spinning reserve is that it should be greater or equal to the capacity of the largest online generator. However, this process may not be economically prudent based on the dynamism of the pertaining situation. So, the most crucial problem depends on the minimum generating capacities that should be reserved. Both the reliability and the operating cost of the grid should be considered when deciding on the amount of the minimum spinning reserve capacity. It is necessary to set the reserve capacity at optimal value, such that the operations of the grid are most economical and reliable. The proposed model provides grid operators with the option of the various combinational of SR and BESS capacities. For instance, from Figure 11, to restore the declining system frequency to the nominal value, the grid operator can choose from setting up a 6% spinning reserve capacity with 50 MW BESS capacity or arrange for only a 10% spinning reserve capacity. The decision to select any one of these two options is then at the behest of the grid operator considering the economic constraints. With the considerations of the generators' fuel costs, the available capacity of BESS and start-up, stop-up, and other considerable parameters of the generators, the optimal ensemble of SR and BESS could be selected. The frequency performance index provides the candidate solutions for such optimization for decision making.

Case 3: Reliability Assessments with Optimal SR
The traditional criterion for setting the minimum amount of spinning reserve is that it should be greater or equal to the capacity of the largest online generator. However, this process may not be economically prudent based on the dynamism of the pertaining situation. So, the most crucial problem depends on the minimum generating capacities that should be reserved. Both the reliability and the operating cost of the grid should be considered when deciding on the amount of the minimum spinning reserve capacity. It is necessary to set the reserve capacity at optimal value, such that the operations of the grid are most economical and reliable. The proposed model provides grid operators with the option of the various combinational of SR and BESS capacities. For instance, from Figure 11, to restore the declining system frequency to the nominal value, the grid operator can choose from setting up a 6% spinning reserve capacity with 50 MW BESS capacity or arrange for only a 10% spinning reserve capacity. The decision performance index provides the candidate solutions for such optimization for decision making.

Case 3: Reliability Assessments with Optimal SR
The traditional criterion for setting the minimum amount of spinning reserve is that it should be greater or equal to the capacity of the largest online generator. However, this process may not be economically prudent based on the dynamism of the pertaining situation. So, the most crucial problem depends on the minimum generating capacities that should be reserved. Both the reliability and the operating cost of the grid should be considered when deciding on the amount of the minimum spinning reserve capacity. It is necessary to set the reserve capacity at optimal value, such that the operations of the grid are most economical and reliable. The proposed model provides grid operators with the option of the various combinational of SR and BESS capacities. For instance, from Figure 11, to restore the declining system frequency to the nominal value, the grid operator can choose from setting up a 6% spinning reserve capacity with 50 MW BESS capacity or arrange for only a 10% spinning reserve capacity. The decision to select any one of these two options is then at the behest of the grid operator considering the economic constraints. With the considerations of the generators' fuel costs, the available capacity of BESS and start-up, stop-up, and other considerable parameters of the generators, the optimal ensemble of SR and BESS could be selected. The frequency performance index provides the candidate solutions for such optimization for decision making.

Conclusions
This paper proposes a technique to evaluate short-term system frequency operational reliability considering system dynamics. The reliability index under various penetrations of wind speed, spinning reserve capacities, and BESS capacities was formulated considering frequency regulation processes. A statistical distribution of the system frequency was represented as a histogram. A case study was designed with the adaptation of the RTS model to show the efficacy of the proposed model. The simulation results showed that the combined operations of spinning reserve and BESS could be used in load frequency control. Detailed findings related to computational time shows how the proposed framework could be used for operational reliability because of the short computational time capabilities. The results also show that reliability improvement from the WTGs without sufficient spinning reserve can be made possible with BESS support. An ensemble of the spinning reserve with BESS can significantly improve the utilization of wind power. The stochastic description of the system frequency under high wind power penetration provides essential information for system operators and planners to make the most reliable and economic decisions about the choice of spinning reserve capacity to deploy.