Renewable Energy Generation Assessment in Terms of Small-Signal Stability

: The popularity and role of renewable energy in the power grid are increasing nowadays as countries are shifting to cleaner forms of energy. This brings new challenges in maintaining a secure and stable power system, as renewable energy is known to be intermittent in nature and may introduce stability issues to the grid. In this paper, a screening framework of renewable energy generation scenarios is proposed to determine which power system conditions and scenarios will make the system unstable. The scenario screening framework is based on a sensitivity analysis of the system eigenvalues with respect to the renewable energy penetration to the system. The average scheduled renewable energy output, forecasting error standard deviation, average forecasting error, and bus location of the renewable energy source were used to deﬁne a renewable energy generation scenario. Depending on the amount and variability of renewable energy, there is a possibility for a critical eigenvalue to cross the imaginary axis. The estimated eigenvalue location resulting from the penetration of variable renewable energy is computed by adding the computed eigenvalue sensitivity to the initial operating point. If any of the estimated system eigenvalues cross the imaginary axis, the power system might be unstable in this scenario, so it requires more detailed simulations and countermeasures. Renewable energy forecasting was done using the long short-term memory model, and the proposed method was simulated using the IEEE 39-bus New England test system. The results of the proposed method were veriﬁed by comparing the simulation results to the eigenanalysis solution. The obtained results have shown that the proposed method can determine whether the renewable energy generation scenario is critical in power system operation.


Introduction
As the safe global warming limit has changed from 2 to 1.5 • C above preindustrial levels [1], tightened policies and procedures from governments around the world are expected to be passed in order to eradicate CO 2 emissions and to shift fully to renewable energy (RE). Because of the heightened interest in integrating RE sources such as wind and solar energy in the power system, the need for integration studies of these renewable sources also increases. This is because of the intermittent nature that is associated with RE sources, a secure power system should be able to maintain stability, even with fluctuating RE output and be able to withstand the resulting contingencies or disruptions. Security assessment is done to determine if the power system can withstand disturbances without suffering a drop in the customer service level provided by the transmission and distribution companies [2]. Dynamic security assessment (DSA) is a type of security assessment that studies the security of the transient and dynamic responses to system disturbances in different timescales [3]. With the help of an online DSA [4,5], power system planners and operators can monitor the power system security in real-time, and they have the ability to act immediately in case of a critical disturbance occurring in the system. The transient energy function [6] and trajectory sensitivities [7] have both been used to analyze dynamic security assessments of the system. The difference between the two methods is that the analysis becomes more complex while using transient energy functions especially when dealing with big differential and algebraic equation (DAE) models of power systems [7]. The main features of an online DSA system include contingency screening, time simulations for stability assessment, calculation of power transfer limits, and concurrent processing of contingencies [4]. Contingency analysis is done on a single-base case scenario to determine how the system will react to any contingency event, which can be a planned or unplanned loss of a system element. The contingency screening function reduces the number of contingencies to be considered by the system operators by identifying the list of contingencies with a high likelihood of occurrence. Similarly, a scenario screening function identifies the base case scenarios that will be critical to the stability of the power system. There are different power system parameters in a base case scenario that can be analyzed, but for the purposes of this paper, scenarios will be limited to an analysis of RE generation scenarios.
Depending on the amount of RE generation forecasted for the next day, system operators will plan a day ahead which generators will be operated and at what generation output each generator is contracted to produce. During actual operations, the actual RE generation may differ from the forecasted RE output because of unforeseeable and uncontrollable external factors. The fluctuations in actual RE generation may result in a lower or higher total generated energy. This intermittency of RE can cause small disturbances to the power system because the power balance between the generated and distributed energy will not be equal. Small-signal stability refers to the ability of the power system to maintain synchronism when subjected to small disturbances [8].
In this paper, RE is considered as a small disturbance because the equations representing the power system can be linearized for the purpose of small-signal analysis [8]. In order to describe the small-signal stability of the power system for RE generation scenarios, eigenvalues of the power system must be analyzed. If an eigenvalue is found at the right-half of the complex plane, the system is unstable. Various methods were proposed to determine the trajectory of the eigenvalues as a system parameter is increased [9][10][11]. In [9], eigenvalue tracing is done by combining the invariant subspace method and the predictor-corrector scheme. An integration-based approach is used in [10] to trace critical eigenvalues and is combined with an index to determine which eigenvalues are considered critical. While in [11], the invariant subspace method is used with the projected Arnoldi method to trace the critical eigenvalues. The difference of trajectory tracing with the proposed method is that trajectory tracing involves a single parameter to be increased until a bifurcation occurs. This will not be applicable for scenario screening because each scenario is tested for its stability according to the small-signal stability rules. A more applicable approach is to determine the sensitivity of the eigenvalues. An eigenvalue sensitivity analysis was proposed in [12][13][14][15] to depict how the eigenvalues are affected by any change in the parameters of the power system. The method for determining the sensitivity of eigenvalues for sparse formulations of the power system equations has been presented in [12]. Eigenvalue sensitivity analyses, with respect to operating parameters [13] and power system parameters [14], have also been proposed in the literature. The augmented matrix using only dominant eigenvalues and their corresponding left and right eigenvectors were used in [15]. Lower-upper (LU) factorization, sparse matrices, and parallel computing techniques were used in computing the eigenvalue sensitivity of operating parameters [16]. Aside from power systems, the eigenvalue sensitivity was also used in parametric optimization. For instance, eigenvalue sensitivity was used to optimize the parameters of a doubly-fed induction generator (DFIG) wind turbine system [17,18]. Variation of different system parameters of microgrids with constant power loads was analyzed using eigenvalue sensitivity analysis in [19]. In the methods mentioned above, the screening of scenarios by analyzing the sensitivity of the eigenvalues with respect to the variation of RE has not been studied for power systems.
In daily operations, RE generation at each time period needs to be estimated based on day-ahead forecasting, so that operating points are established for security assessment. In the literature, several regression techniques have been adopted for RE forecasting. Autoregressive moving average (ARMA) and autoregressive integrated moving average (ARIMA) have become two of the most popular methods that can be used for linear time series analysis [20]. Aside from regression techniques, artificial intelligence (AI) techniques and hybrid methods were also proposed in the literature. For photovoltaic (PV) forecasting, the long short-term memory (LSTM) recurring neural network, which is a method under AI techniques, has become popular to implement because of its better performance compared to other conventional implementations [21]. This paper employs LSTM as the method for forecasting RE day-ahead generation because, aside from the historical PV output, historical cloud cover data are also utilized. Even though the accuracy of the forecasting techniques has improved over time, the actual generated output might still be different from the forecasted output because of the limitations of forecasting as well as unpredictable external sources affecting the generation of RE. Thus, for secure operation of the power system with RE resources, diverse RE generation scenarios need to be assessed in terms of several security issues. RE generation scenarios for examination are established for a specific time period using RE bus location, average forecasting error, forecasting error standard deviation, and average scheduled RE output. Depending on these factors, certain scenarios might be critical and result in insecure or unstable conditions. In order to keep the power system secure and stable under any circumstances, critical scenarios must be identified by power system planners to come up with countermeasures against the corresponding insecure phenomena. This paper presents a framework for screening of RE generation scenarios in terms of small-signal stability, and it is based on eigenvalue sensitivity analyses with respect to RE generation. Eigenanalysis for the initial operating point is performed first, and sensitivities of critical eigenvalues with respect to RE variation are evaluated to estimate the movement of the eigenvalues and to check whether the resulting locations violate small-signal stability criteria. If some of the estimated eigenvalue locations are not in the acceptable region, the corresponding RE generation scenario is included in the critical scenario list and should be reassessed using a detailed stability analysis. The results of the proposed method were verified by comparing it with an eigenanalysis of the newly established operation with the corresponding RE generation scenario. Instead of performing an eigenanalysis for all possible operating points, the proposed screening framework can quickly determine whether each RE generation scenario might cause insecurity and instability. By implementing the proposed method in screening RE generation scenarios, only the critical scenarios need to be analyzed in depth by system operators for countermeasure determination. The contributions of this paper can be summarized as follows: (i) establishment of the framework to screen RE generation scenarios in terms of small-signal stability, (ii) inclusion of RE generation parameters such as average forecasting error, forecasting error standard deviation, and average scheduled RE output in the computation of the RE variability factor, (iii) application of eigenvalue sensitivity with respect to RE generation to determine small-signal stability, and (iv) formulation of the uncertainty of the real part of critical eigenvalues with respect to RE generation.
The paper is organized as follows: In Section 2, the mathematical background and the formulation of the system model, system state matrix and the eigenvalue sensitivity concept are presented. In Section 3, the proposed scenario screening framework based on the eigenvalue sensitivity analysis is presented. In Section 4, a day-ahead RE forecast is made using a LSTM model for RE generation scenarios. In the same section, the proposed method is applied to the New England Test System (IEEE-39 bus test system) to demonstrate its capability of determining the small-signal stability of the system for the five given scenarios. In Section 5, the conclusions of the study are presented.

Mathematical Background and System Formulation
This section presents the formulation and mathematical background of the system model, system state matrix, and the eigenvalue sensitivity concept used in this paper. The descriptions of the variables and parameters in the formulation of the system model are shown in Table 1. The units of them are also given in Table 1. The test system used for this analysis was the modified New England 39 bus test-system including RE generation, the one-line diagram of which is shown in Figure 1. There are 39 buses and 10 synchronous generators in the New England test system. The five locations for the RE source that were considered in the simulations are also illustrated in Figure 1. The original case power flow without RE generation and dynamic data are provided in [22,23]. The state equations for each state variable are discussed in Section 2.1. The state equations are also called the differential equations of the system. The active and reactive power network equations for each bus are discussed in Section 2.2. The power network equations are also called the algebraic equations of the system. The formulation of the differential algebraic equations (DAEs) model of the system is shown in Section 2.3. From the DAE model, the system state matrix and the concept of eigenvalue sensitivity can be derived. This will serve as the foundation for the scenario screening framework that will be presented in Section 3.

State Variables
In this paper, each synchronous nonreference generator was represented by 10 state variables. The 10 state variables were rotor angle ( ) , machine frequency ( ), quadrature axis emf ( ), direct axis emf ( ), field winding voltage ( ), AVR output ( ), exciter feedback ( ), prime mover mechanical power output ( ), and generator steam valve opening ( ). For reference generators, there will only be 9 state equations. There will be no state equation for the rotor angle of the reference generator since the rotor angles of the nonreference generators are referenced to the rotor angle of the reference generator.
The differential equations for the rotor angle ( ) , machine frequency ( ), quadrature axis emf ( ), and direct axis emf ( ) of the synchronous generators are written below [22]. The rotor angle of the mth generator is chosen as the reference angle of the system. In this paper, the generator in bus 39 was chosen as the reference generator.
In the two-axis model shown in Equations (1)-(4), the stator transients were ignored. The machine direct axis ( ) and quadrature axis currents ( ) are expressed [22] in the succeeding equations. , , and in the following equations refer to the bus voltage magnitude, bus voltage angle, and rotor angle respectively.

State Variables
In this paper, each synchronous nonreference generator was represented by 10 state variables. The 10 state variables were rotor angle (δ i ), machine frequency (ω i ), quadrature axis emf (E qi ), direct axis emf (E di ), field winding voltage (E f di ), AVR output (V ri ), exciter feedback (R f i ), prime mover mechanical power output (P mi ), and generator steam valve opening (µ i ). For reference generators, there will only be 9 state equations. There will be no state equation for the rotor angle of the reference generator since the rotor angles of the nonreference generators are referenced to the rotor angle of the reference generator.
The differential equations for the rotor angle (δ i ), machine frequency (ω i ), quadrature axis emf (E qi ), and direct axis emf (E di ) of the synchronous generators are written below [22]. The rotor angle of the mth generator is chosen as the reference angle of the system. In this paper, the generator in bus 39 was chosen as the reference generator. . . .
In the two-axis model shown in Equations (1)-(4), the stator transients were ignored. The machine direct axis (I di ) and quadrature axis currents (I qi ) are expressed [22] in the succeeding equations. V, θ, and δ in the following equations refer to the bus voltage magnitude, bus voltage angle, and rotor angle respectively. The nonreference generator (i = 1, . . . , m − 1) is as follows: The reference generator (i = m) is as follows: For the excitation control system of the synchronous generator, IEEE type DC-1 excitation system was used. The model for this excitation system [22] is shown in Equations (9)-(11). The following equations are for the field winding voltage (E f di ), AVR output (V ri ), and exciter feedback (R f i ): . .
The differential equation for the prime mover and speed governor model are shown in Equations (12) and (13) [22]. P mi refers to the prime mover mechanical power output, while µ i is the mathematical expression for the physical representation of the generator steam valve opening. To compensate for any governor speed deviation from the reference speed (ω re f ), there is a corresponding change in µ i depending on the speed difference and the droop constant (R i ). The droop constant represents the speed governor inherent speed-droop characteristics [22]. A typical value for the droop constant is 5% droop. .

Algebraic Variables
Consider a system with n buses, the network power equations for each bus of the system are shown as: In this paper, the test system used was the New England test system and was composed of 39 buses. Since there is one algebraic equation each for the active power and reactive power, there will be 78 power network equations for this test system. By following this, a test system with n buses is expected to have 2n power network equations. P gi and Q gi refer to the active power output and reactive power output of the generator at bus i. P gi and Q gi are zero if there is no generator in bus i. Otherwise, the equations for P gi and Q gi are shown below [22]. In (14), P REi stands for RE generation at bus i. In this paper, P REi was considered as a given parameter at a specific time period. Similar to the notation in the previous subsection, the variables V, θ, and δ refer to the bus voltage magnitude, bus voltage angle, and rotor angle respectively.
The nonreference generator (i = 1, . . . , m − 1) is as follows: The reference generator (i = m) is as follows: P ti and Q ti refer to the total network active power and reactive power injections at the ith bus respectively. These equations are similar to the equations used in the power flow analysis [24]. Y ij and γ ij in the following equations refer to the bus admittance magnitude and angle respectively.
The active power load model and the reactive power load model are given as shown in Equations (22) and (23) respectively [22]. P lo,i , Q lo,i , and V o,i refer to the initial load active power, initial load reactive power, and initial load voltage respectively. KP, KQ, and K f req are load factors.
By inspecting the active power and reactive power network equations, it can be seen that the equations are in terms of V and θ. Because of this, V and θ are called the algebraic variables of the system.

System State Matrix and Eigenvalue Sensitivity
The general form of the differential and algebraic equations (DAEs) describing the power system is shown as [22]: where f and g refer to the set of differential and algebraic equations respectively. f and g can be represented in terms of the state variables (x), algebraic variables (y), and the control parameter (a). From Sections 2.1 and 2.2, x and y are defined as the set of differential and algebraic variables. In this paper, the control parameter is the RE variation in the system.
The solution of the DAE model is the equilibrium point of the system. By linearizing the DAE model around the equilibrium point, the resulting equations are shown in state matrix form in Equation (28). The matrix elements f x , f y , g x , and g y are called the Jacobian matrices of the system. The symbol ∆ represents parameter deviation.
. (28) By assuming that the Jacobian matrix g y is nonsingular, the variable ∆y can be eliminated ∆x. The resulting equation is shown in Equation (30). From Equation (30), the equation for the system state matrix A s is derived. .
The sensitivity of the eigenvalue λ with respect to system parameter α is given by Equation (32) [8]. Eigenvalue sensitivity depicts how the system eigenvalues are affected by a change in a system parameter.
In Equation (32), it is shown that the eigenvalues and eigenvectors of the system state matrix A s are used to determine the eigenvalue sensitivity of the system. ψ pertains to the left eigenvector, while φ pertains to the right eigenvector of the system state matrix A s .

Renewable Energy (RE) Generation Scenario Screening Framework Based on Eigenvalue Sensitivity
The overview of the proposed scenario screening framework is shown in Figure 2. The generation dispatch was scheduled based on the day-ahead forecast of RE. The initial eigenanalysis was then analyzed according to the scheduled dispatch. Based on the initial eigenanalysis, the eigenvalues with the lowest damping ratios and eigenvalues with the largest real parts were considered as critical eigenvalues. The critical eigenvalues were monitored closely because they were the most probable candidates to cross the imaginary axis first. The eigenvalue sensitivity with respect to the renewable energy was simulated based on the eigenvalue and eigenvector information from the initial eigenanalysis. This is discussed further in Section 3.1.
The sensitivity of the eigenvalue λ with respect to system parameter is given by Equation (32) [8]. Eigenvalue sensitivity depicts how the system eigenvalues are affected by a change in a system parameter.
In Equation (32), it is shown that the eigenvalues and eigenvectors of the system state matrix are used to determine the eigenvalue sensitivity of the system. pertains to the left eigenvector, while pertains to the right eigenvector of the system state matrix .

Renewable Energy (RE) Generation Scenario Screening Framework Based on Eigenvalue Sensitivity
The overview of the proposed scenario screening framework is shown in Figure 2. The generation dispatch was scheduled based on the day-ahead forecast of RE. The initial eigenanalysis was then analyzed according to the scheduled dispatch. Based on the initial eigenanalysis, the eigenvalues with the lowest damping ratios and eigenvalues with the largest real parts were considered as critical eigenvalues. The critical eigenvalues were monitored closely because they were the most probable candidates to cross the imaginary axis first. The eigenvalue sensitivity with respect to the renewable energy was simulated based on the eigenvalue and eigenvector information from the initial eigenanalysis. This is discussed further in Section 3.1. RE generation scenarios were defined by the RE location, average forecasting error, forecasting error standard deviation and average forecasted RE output. The computation for the input was based on the RE generation scenario variables. The formulations for the and RE generation scenario variables are shown in Section 3.2.
The estimated eigenvalues were calculated by adding the eigenvalue sensitivity information to the initial critical eigenvalues. If an estimated eigenvalue was found in the right complex plane, then the system might be unstable in this RE generation scenario, so it requires more detailed The estimated eigenvalues were calculated by adding the eigenvalue sensitivity information to the initial critical eigenvalues. If an estimated eigenvalue was found in the right complex plane, then the system might be unstable in this RE generation scenario, so it requires more detailed simulation and countermeasure determination. Likewise, if all eigenvalues were found in the left complex plane, then the system is said to be stable in this RE generation scenario.
The formulation for calculating the uncertainty of the real part of critical eigenvalues is discussed in Section 3.3. Additional discussion and a detailed flowchart are presented in Section 3.4.

Eigenvalue Sensitivity with Respect to Variation of RE
The sensitivity of the eigenvalues with respect to RE was calculated because it could determine the effect of the variable RE to the eigenvalues of the power system. The eigenanalysis of the base case scenario was evaluated first to determine the initial eigenvalues and eigenvectors of the system. The eigenvalues and the eigenvectors were used to determine the sensitivity of eigenvalues. From Equation (32), the sensitivity of the critical eigenvalue λ c with respect to RE generation P RE can be derived. The resulting equation is shown in Equation (33).
To determine the partial derivative of the system state matrix A s with respect to P RE , the derivative of the system state matrix with respect to V and δ were used in this method. V refers to the voltage magnitude and δ refers to the phase angle.
When connected to the grid, the added RE acts to decrease the active power being supplied by the generators to the loads. P RE is added in the Jacobian matrix used for power flow analysis. ∂|V| ∂P RE and ∂δ ∂P RE can be derived from the inverse of the Jacobian matrix.
By taking the derivative of the system state matrix with respect to V, Equation (31) is changed to the following equation: In [25], a simplified expression for Likewise, an analysis is done to get the derivative of the system state matrix A s with respect to δ.

RE Generation Scenario
A RE generation scenario is a set of conditions that describes how much the actual RE output deviates from the forecasted RE output. In some cases, the actual output differs from the forecasted output because of the intermittency of RE. Obviously, several RE generation scenarios need to be checked in operation. The framework of this paper can quickly evaluate the impact of each scenario, in terms of small-signal stability, to screen critical RE generation scenarios.
Each RE generation scenario is defined by RE location, average forecasting error (mean absolute percent error, MAPE), forecasting error standard deviation (σ err ), and average forecasted RE output (P f or RE ). The location of the RE source is included as part of the scenario because it also affects the generator loading and the exchange of power within the grid.
PE k,t refers to the percent error between the forecasted (P f or RE ) and actual (P act RE ) RE output for time slot t in the kth RE generation scenario.
The formulations for solving the MAPE, σ err , and P f or RE for each RE generation scenario are shown in the following equations.
The corresponding ∆P RE for each RE generation scenario can be evaluated using Equation (44). When setting K to 3, 99.7% of all probable forecasting error margins will be covered in the analysis.
The estimated critical eigenvalues can be calculated using Equation (45). The calculated eigenvalue sensitivity ∂λ c ∂P RE from Equation (33) is multiplied by the change in the RE generation, ∆P RE . This term is added to the initial critical eigenvalue λ c to estimate the movement of the eigenvalue.
Using Equation (45), all critical eigenvalues are estimated with respect to ∆P RE . Then, the framework checks whether there are any estimated eigenvalues in the right side of the complex plane. If there are eigenvalues in the right side of the complex plane, then the system might be unstable in this RE generation scenario, and it requires more detailed analysis in terms of small-signal stability.

Uncertainty of the Real Part of Critical Eigenvalues
In [26], a framework for quantifying the uncertainty of the transmission reliability margin was proposed with respect to uncertain parameters. In this paper, a similar formulation was employed to evaluate the uncertainty of the real part of the critical eigenvalues with respect to the RE generation scenarios. The formulation for the uncertainty U c of the real term of a critical eigenvalue, λ c , is shown in Equation (46). In the evaluation, the sensitivity of the real part with respect to renewable energy and the variance of the renewable energy are both used.
where l is the total number of considered RE generation scenarios.

Additional Discussion on the Proposed Framework
As shown in Equation (45), the sensitivity of the eigenvalues is used to calculate the estimated eigenvalues. Figure 3 depicts how eigenvalue sensitivity can be used to screen RE generation scenarios. The blue circle in Figure 3 represents the critical eigenvalues of the system. The critical eigenvalues can either be eigenvalues with the lowest damping ratios or eigenvalues with the largest real parts. The black dashed lines represent the damping ratio margin that will serve as the criteria to determine which eigenvalues are considered as critical. Depending on the variability of RE in the system for a given scenario, the system eigenvalues may cross the imaginary axis and move to the right half of the complex plane. The area within the red circle represents the possible location of the eigenvalues as a result of the added variable RE. If the eigenvalues are found to be at the right side of the complex plane or if they cross the imaginary axis, the system is unstable for that scenario. location of the eigenvalues as a result of the added variable RE. If the eigenvalues are found to be at the right side of the complex plane or if they cross the imaginary axis, the system is unstable for that scenario. A flowchart is presented in Figure 4 to summarize the process of screening RE generation scenarios. A flowchart is presented in Figure 4 to summarize the process of screening RE generation scenarios.  A flowchart is presented in Figure 4 to summarize the process of screening RE generation scenarios.

Simulation Results and Analysis
As in Figure 1, the modified New England system was adopted in the simulation. The eigenvalues of the New England Test System without any RE penetration are shown in Figure 5 below. The eigenvalue analysis shows that the system had 89 eigenvalues: 41 of which were real

Simulation Results and Analysis
As in Figure 1, the modified New England system was adopted in the simulation. The eigenvalues of the New England Test System without any RE penetration are shown in Figure 5 below. The eigenvalue analysis shows that the system had 89 eigenvalues: 41 of which were real eigenvalues and the other 48 were complex eigenvalues (24 complex eigenvalue pairs). Eighty-nine eigenvalues were expected because each nonreference generator was represented by 9 state equations, and the reference generator was represented by 8 state equations. eigenvalues and the other 48 were complex eigenvalues (24 complex eigenvalue pairs). Eighty-nine eigenvalues were expected because each nonreference generator was represented by 9 state equations, and the reference generator was represented by 8 state equations.
In Figure 5a, all eigenvalues are shown. The eigenvalues were all in the left-hand side of the imaginary axis. Hence, in this case, the system was stable. The system will be unstable if any of the system eigenvalues cross the imaginary axis and move to the right-hand side of the complex plane. Some of the closest eigenvalues near the origin are shown in Figure 5b. These eigenvalues were considered as critical eigenvalues and may have the highest probability of crossing the imaginary axis first, depending on the amount of RE added to the system. The type of RE considered in the simulation for this paper was solar energy. The long short-term memory (LSTM) model was used in this paper to predict the photovoltaic (PV) cell system output based on historical one-year PV output data and cloud cover data. The standard deviation of the error between the actual and predicted PV output ( ), average forecasted PV output ( ), and average forecasting error (MAPE) were used to define the scenario conditions for the proposed method. A sample one-week snapshot of the training data is shown in Figure 6. The total PV output per day is also shown in Figure 6. While in Figure 7, a sample prediction and actual curve of the total PV plant for a sunny and cloudy day are shown. The amount of cloud cover affects PV generation output, and this makes it harder for system operators to forecast PV In Figure 5a, all eigenvalues are shown. The eigenvalues were all in the left-hand side of the imaginary axis. Hence, in this case, the system was stable. The system will be unstable if any of the system eigenvalues cross the imaginary axis and move to the right-hand side of the complex plane. Some of the closest eigenvalues near the origin are shown in Figure 5b. These eigenvalues were considered as critical eigenvalues and may have the highest probability of crossing the imaginary axis first, depending on the amount of RE added to the system.
The type of RE considered in the simulation for this paper was solar energy. The long short-term memory (LSTM) model was used in this paper to predict the photovoltaic (PV) cell system output based on historical one-year PV output data and cloud cover data. The standard deviation of the error between the actual and predicted PV output (σ err ), average forecasted PV output (P f or RE ), and average forecasting error (MAPE) were used to define the scenario conditions for the proposed method. A sample one-week snapshot of the training data is shown in Figure 6. The total PV output per day is also shown in Figure 6. While in Figure 7, a sample prediction and actual curve of the total PV plant for a sunny and cloudy day are shown. The amount of cloud cover affects PV generation output, and this makes it harder for system operators to forecast PV accurately [27]. The type of RE considered in the simulation for this paper was solar energy. The long short-term memory (LSTM) model was used in this paper to predict the photovoltaic (PV) cell system output based on historical one-year PV output data and cloud cover data. The standard deviation of the error between the actual and predicted PV output ( ), average forecasted PV output ( ), and average forecasting error (MAPE) were used to define the scenario conditions for the proposed method. A sample one-week snapshot of the training data is shown in Figure 6. The total PV output per day is also shown in Figure 6. While in Figure 7, a sample prediction and actual curve of the total PV plant for a sunny and cloudy day are shown. The amount of cloud cover affects PV generation output, and this makes it harder for system operators to forecast PV accurately [27].  In the first two scenarios, the PV resource was placed first on bus 20 and then changed to bus 8 to compare results. The same RE generation conditions were assumed for the first two scenarios. In the third scenario, the PV resource was moved to bus 4 with higher variability of PV assumed in the system. The results of the proposed method were verified by comparing it with the eigenanalysis results.
In this paper, the least damping ratio criterion, used to classify which eigenvalues are considered critical, was 0.10 or 10%. The damping ratio of each eigenvalue (λ = σ + jω) is given by Equation (47) [8]. For the largest real parts criterion, eigenvalues with real part greater than −0.10 In the first two scenarios, the PV resource was placed first on bus 20 and then changed to bus 8 to compare results. The same RE generation conditions were assumed for the first two scenarios. In the third scenario, the PV resource was moved to bus 4 with higher variability of PV assumed in the system. The results of the proposed method were verified by comparing it with the eigenanalysis results.
In this paper, the least damping ratio criterion, used to classify which eigenvalues are considered critical, was 0.10 or 10%. The damping ratio of each eigenvalue (λ = σ + jω) is given by Equation (47) [8]. For the largest real parts criterion, eigenvalues with real part greater than −0.10 were considered as critical eigenvalues.
The results of the simulations for the three scenarios are summarized in Tables 1-3. The results shown in the initial critical eigenvalues column are based on the scheduled dispatch. If an eigenvalue had a damping ratio of less than 0.10 or 10%, or if the real part of the eigenvalue was greater than −0.1, the eigenvalue was considered as a critical eigenvalue. The corresponding damping ratio for each complex eigenvalue is shown in the next column. Critical eigenvalues were given special importance in the analysis because they were the most probable candidates to cross the imaginary axis first. If any of these eigenvalues crossed the imaginary axis to the right-hand side of the complex plane, then the system might be unstable in this scenario. The estimated eigenvalues column shows the result of the proposed method by adding the eigenvalue sensitivity with respect to PV variability to the initial critical eigenvalues. These results were then verified with the results of the eigenanalysis with the corresponding PV generation. In the results shown in Table 2, the assumed condition was that the scheduled dispatch included a mean PV output of 5.41 pu in bus 20. The values shown in the estimated eigenvalues column were evaluated using the proposed method and without solving for the eigenanalysis solution. Based on the results, the proposed method showed that the eigenvalues were still in the left-hand side of the complex plane. As such, the system was considered stable in this scenario. To verify this result, PV output was then added to the system, and an eigenanalysis was done to determine the eigenvalues of the system in this case. The results are shown in the eigenanalysis verification column. The eigenanalysis result also showed that the eigenvalues were in the left-hand side of the complex plane. Hence, Scenario 1 was considered as a stable scenario as both estimated eigenvalues and the eigenanalysis verification showed the same result.
To test if the method will be applicable in other locations in the system, the same analysis was done with the same scenario conditions but with a different PV location. In the simulation results shown in Table 3, the same scenario conditions were used, but the PV location was changed to bus 8. It showed that the results of the proposed method agreed with the results of the eigenanalysis verification for this scenario as well. The system was shown to be stable for the assumed conditions in this scenario. In the last scenario, the PV bus was changed to bus 4, but with a higher scheduled PV output mean and higher PV output variability. The scheduled dispatch included a mean PV output of 5.65 pu. The results for the analysis done with this scenario are shown in Table 4. For increased variability of the PV output, the proposed method showed that the system was stable in this scenario since all eigenvalues were in the left-hand side of the complex plane. The results agreed with the eigenanalysis verification in this scenario. In the evaluation of the uncertainty of the real part of the critical eigenvalues, two (2) additional scenarios were considered, as shown below: By considering the five scenarios, the uncertainty of the real part of the critical eigenvalues can be evaluated using Equation (44). The real part of the five largest and smallest real parts of the critical eigenvalues for each scenario are tabulated in Table 5. The uncertainty for the real part of each critical eigenvalue are shown in the last column of the table. In Scenario 4, the RE generation scenario was similar to Scenario 3, but the PV location was moved from bus 4 to bus 29. After evaluating the sensitivity of the critical eigenvalues with respect to PV generation, the proposed scenario screening framework determined that Scenario 4 was a stable scenario as all of the estimated critical eigenvalues were in the left complex plane.
While in Scenario 5, the scheduled PV generation was higher compared to the other scenarios, and the PV location was moved to bus 39. After applying the proposed method, it was determined that there was an estimated critical eigenvalue in the right complex plane. If an eigenvalue is found in the right complex plane, the system might be unstable in this scenario. The results of both scenarios agreed with the eigenanalysis verification. The system eigenvalues of the critical scenario are shown in the left image of Figure 8. The eigenvalues near the imaginary axis are shown in the right image of Figure 8 for easier visualization. It can be seen from this image that a complex eigenvalue pair crossed the imaginary axis, therefore making the scenario unstable and critical. to PV generation, the proposed scenario screening framework determined that Scenario 4 was a stable scenario as all of the estimated critical eigenvalues were in the left complex plane. While in Scenario 5, the scheduled PV generation was higher compared to the other scenarios, and the PV location was moved to bus 39. After applying the proposed method, it was determined that there was an estimated critical eigenvalue in the right complex plane. If an eigenvalue is found in the right complex plane, the system might be unstable in this scenario. The results of both scenarios agreed with the eigenanalysis verification. The system eigenvalues of the critical scenario are shown in the left image of Figure 8. The eigenvalues near the imaginary axis are shown in the right image of Figure 8 for easier visualization. It can be seen from this image that a complex eigenvalue pair crossed the imaginary axis, therefore making the scenario unstable and critical.
(left) (right) The significance of the uncertainty results showed that, by considering the variation in the PV output in any of the considered PV locations and the sensitivity of the eigenvalues, the real part of the critical eigenvalues varied by the uncertainty value. This is a measure of how the critical eigenvalues were affected on a system-level perspective.
The results of the proposed scenario screening framework and their run times are tabulated in The significance of the uncertainty results showed that, by considering the variation in the PV output in any of the considered PV locations and the sensitivity of the eigenvalues, the real part of the critical eigenvalues varied by the uncertainty value. This is a measure of how the critical eigenvalues were affected on a system-level perspective.
The results of the proposed scenario screening framework and their run times are tabulated in Table 6. The simulations were done on a computer with a 3.30 GHz Intel processor and 8 GB RAM. All scenarios, except Scenario 5, were stable scenarios since the estimated critical eigenvalues in these scenarios were found in the left complex plane. Since there was a critical eigenvalue that crossed the imaginary axis in Scenario 5, the system might be unstable in this scenario, according to small-signal stability rules, and it requires further detailed simulations. The run times of the proposed method were similar to those of the analytical methods using the original New England test system in [16]. However, it should be noted that the proposed framework had a different applicability, which was to quickly screen critical RE generation scenarios in terms of small-signal stability, and not only to calculate eigenvalue sensitivity.
For critical scenarios determined by the proposed framework, countermeasure determination was needed with several control parameters such as active power injection with energy storage systems (ESSs) and generation re-dispatch. As in [28,29], the ESS capability of local energy supplying systems including microgrids could increase the operational efficiency by reducing uncertainty of RE. In addition, it could improve system stability, which might be threatened by RE intermittence, if the local systems with ESS followed the signals to meet the required control amount as ancillary service providers. As in [30,31], the concept of microgrids also included the role of energy hub and energy router. If power injections from those energy hubs are properly scheduled to minimize the negative impacts by RE intermittence in several security issues, system security can be enhanced as well. The proposed method of this paper can be used to quickly filter alarming situations in term of small-signal stability; hence, adequate control signals of required power injection are determined.

Conclusions
This paper presents the framework for screening RE generation scenarios in terms of small-signal stability. Instead of performing eigenanalyses for all possible operating points, the proposed screening framework adopted eigenvalue sensitivity with respect to renewable generation and checked if the system might violate the small-signal stability rules by RE fluctuation. The case studies using the proposed method showed that it could quickly determine whether there might be the possibility of instability for several RE scenario conditions such as different bus locations and different RE variability levels. For all the scenarios in the simulation, the proposed method showed that the results agreed with the eigenanalysis verification with the corresponding operating points. The proposed method was very helpful to find critical scenarios that may possibly lead to severe conditions. By using the information from the method, a more detailed analysis may be done in terms of small-signal stability, especially for unstable scenarios. Thus, the method can be used in system planning and operational planning stages to evaluate the effect of different RE generation scenarios and to quickly filter critical ones. As a result, system engineers could come up with adequate countermeasures against those cases to improve the grid reliability and to mitigate the security risks.