Dynamic Performance Modeling and Analysis of Power Grids With High Levels of Stochastic and Power Electronic Interfaced Resources

This article examines the emerging challenges in modeling and analyzing the electric power system due to the widespread growth of variable renewable energy (VRE), particularly in the form of distributed energy resources (DERs), which are displacing traditional large power plants. Many of these resources are connected to the system through power electronic interfaces, also known as inverter-based resources (IBRs), which are reshaping the system dynamics and lowering the grid strength and inertia. Understanding the dynamic behavior of the power system should be critical to addressing the potential stability concerns, refining the grid requirements, and developing effective and reliable measures among many alternatives. However, conventional methodologies for resource integration and network expansion studies, as well as application-specific electromagnetic transient (EMT) studies, need to be improved. This article thus presents recent academic and industrial efforts to advance the existing approaches, especially by incorporating the uncertainty in model parameters of DERs, variability of VRE, and EMT dynamics of IBRs for the grid planning and operations studies such as the impact of DERs on load modeling and system-wide dynamic performance. In addition, this article showcases recent developments to expand the study boundaries by synergizing the strengths of the industry-accepted approaches along with real system studies for Korea’s electric power systems in particular.

past six years, the continuous increase of distributed PVs in Korea has significantly increased the load variability, heightening the importance of grid flexibility (via fastacting resources [9], [15]) and forecasting accuracy for the net load and renewable generation.
These ongoing changes may threaten power system reliability and stability if not handled adequately in the operations and planning. Understanding the dynamic behavior of the power system should thus be critical to addressing the potential stability concerns, refining the grid requirements, and developing effective and reliable measures among many alternatives [74]. It is equally critical to understand the limitations of existing methods and tools for adequately modeling power system dynamics and accurately assessing the power system security [17], [18], [19], [74]. Conventional practices for the resource integration and network expansion studies, and application-specific electromagnetic transient (EMT) studies are both limited as an integrative approach.
This article thus examines the emerging challenges in modeling and analyzing the electric power system under transformation, as summarized and shown in Fig. 2. Focusing on the bulk power system planning and operations, this article presents selected methods to effectively capture the growing uncertainty, variability, and fast dynamics of modern power systems evolving with high penetrations of VRE, largely interfaced through power electronics. Positive-sequence root-mean-square (rms) model-based studies have been the most commonly used in large-scale power system studies [18], [20], [21], [22]. Given that the rms model-based analysis focuses on the dynamic performance at the transmission level, aggregated models for the distribution systems and loads [23] have been widely used for computational efficiency and convenience. For instance, DERs and BTMGs scattered in the distribution network are often represented by a single generator and encapsulated in the loads [24], so the uncertainty from these models becomes high and needs to be investigated. Also, the dynamic analysis based on the rms model starts from a specific single operating condition; thus, it cannot reflect the inherent variability in renewables: this limitation will be more pronounced when the renewables become dominant energy sources [25]. Moreover, the rms models cannot exhibit the fast dynamic characteristics of power electronics-based equipment [26], [27], [74] because they are principally developed to simulate slow dynamics in the frequency range of 0.1-10 Hz. Accurate representation of these components requires EMT tools to adequately model fast current and voltage controllers, and switching characteristics of power electronics [27], [28]. However, the EMT modeling poses a considerably higher computational burden than the rms modeling [29]. It is then remarkable that advances in computing technology have partially relieved computational concerns and improved the performance of EMT tools significantly [30]. We thus investigate a practical framework for large-scale EMT simulation using parallel computing and techniques to facilitate the overall study procedures. We further introduce a co-simulation platform combining the advantages of rms and EMT studies along with case studies for the future Korea electric power systems [31], [32]. The remainder of this article is structured as follows. The state-of-the-art rms modeling for the inverter-based resources (IBRs) is reviewed in Section II. Section III presents the analysis method for evaluating the impact of model parametric uncertainties, especially due to aggregated representation of DERs via case studies for the Korea power system. The stochastic stability analysis is introduced in Section IV. Section V introduces the EMT model-based analysis and recent efforts extended to large-scale system studies. The conclusions are presented in Section VI.

II. R E N E W A B L E E N E R G Y M O D E L I N G F O R L A R G E -S C A L E P O W E R S Y S T E M D Y N A M I C S T U D I E S
Positive-sequence-based rms modeling has been considered the most commonly accepted approach for representing the dynamic performance of the large-scale power system [18], [20]. The system components, such as generators and loads, are modeled using differential algebraic equations (DAEs) along with power-flow equations as follows:ẋ = F (x, y, λ) (1) where x, y, and λ are dynamic variables (e.g., states of generator models), algebraic variables (e.g., terminal voltage of generator and network variables such as bus voltages and angles), and the model parameters (e.g., machine parameters), respectively. It is worth noting that switching certain variables under specific conditions is a common method used in various power system models to represent the discrete characteristics of the model (e.g., the saturation characteristics of electric machines, transformer tap changing, and control using logic-based switching). An example of utilizing a switching function is the voltage or frequency ride-through of renewable generators to simulate the generation trip when operating outside of specific predefined conditions [33]. When switching occurs at a particular moment in some models, the algebraic equation (2) is changed, allowing for the simulation of nonsmooth behavior. The power system modeled by (1) and (2) is utilized to investigate the dynamic responses of the power system (e.g., frequency and bus voltages) for credible disturbances, including transmission line failure, unexpected generator trip or failure, and temporary or permanent faults at substations. Once the study finds stability concerns for a specific scenario, adequate operating strategies will be designed to avoid problematic operating conditions. The worst case should be a system-wide blackout and prevention and mitigation measures should be investigated in advance. Accurate models for various elements (e.g., generators and transmission facilities) are essential to guarantee the reliability of power system analysis. Among various power system models, this article mainly focuses on the dynamic modeling of renewable energy resources, which are typically represented using either plant-level models [34], [35] or DER models [22] to capture their dynamic behavior. The following provides a brief overview of these two model types. More detailed information on other dynamic models can be found in [18].

A. Plant-Level Renewable Modeling
The plant-level renewable generator (i.e., utility-scale generators) collects power generated from multiple renewable energy resources (e.g., tens to hundreds of wind turbines or solar panels) and supplies the collected active and reactive powers to the grid. The general structure of solar and wind power plants can be found in [36] and [37], respectively. When investigating the dynamic performance of large-scale power systems, it is common to model the plant-level renewables using an aggregation approach (e.g., a single equivalent generator [37] or multiple equivalents by separately aggregating inverters from the same manufacturer [38]) rather than modeling each solar panel or wind turbine of the plant; it is a widely accepted current practice [21], [33], [36]. The validity of the aggregated plant model needs to be verified through field tests or measurement-based validation [33], just as conventional synchronous generator models need to be [39].
To understand the concept of aggregated modeling, Fig. 3 shows a practical utility modeling process for the Jeongam wind power plant in Korea, which consists of 14 type-4 wind turbines. Fig. 3 also shows the approximate location of the Jeongam wind power plant within the Korean transmission map. During this modeling process, the 14 wind turbines are aggregated and represented as a single equivalent generator, as visually shown in Fig. 3. The equivalent generator consists of four functional modules: a plant control module, an electrical control module, a generator converter module, and a turbine. The model has been developed with many years of effort under the leadership of the Western Electricity Coordinating Council (WECC) [35], and continuous efforts are ongoing to verify and improve the validity of the developed models [33], [40]. The modeling process shown in Fig. 3 can be extended to PV power plants and type-3 wind power plants [34], [35].
Power system operators may continuously monitor plant-level renewable generators, check their dynamic characteristics, and validate the derived models. However, it is impossible to precisely model the dynamics of DERs spread across the entire power grid. As a result, a different aggregated modeling approach from plant-level modeling should be required for large-scale power system studies.

B. Modeling of Distributed Energy Resources
As seen and discussed in Fig. 1, DERs, such as BTMG, are frequently modeled together with power demand (or load). This section briefly overviews how power load is Vol. 111  represented in large-scale power system studies, followed by a description of BTMG modeling.
The numerous power loads in the bulk power system are constantly subject to change, making it practically impossible to model and analyze each one individually. Dynamic analysis of large-scale systems aims to determine whether the synchronization of generators and grids is sustained after specific disturbances and whether the system voltage is maintained at an appropriate level. Therefore, examining the detailed dynamic characteristics of all distribution networks is unnecessary (and technically infeasible). Instead, the power load is modeled by aggregating countless loads in the distribution system into several models at the transmission level, as shown in Fig. 4. The various loads of the distribution system in the box are grouped and aggregated as well documented in [18,Chapter 7]. Various aggregated models [41], [42], [43], from simple to complex, have been developed and discussed to enhance modeling accuracy and validity. For example, Fig. 4 shows the composite load model [44] recently developed by the WECC. This model combines four induction motor models, an electronic load model, and a static load model in parallel. The proportion of each model must be carefully determined and validated based on information about the load composition and characteristics. Unlike other electric facilities (e.g., generating plants), the load model structure and parameters cannot generally be tested on-site. Therefore, disturbance data and ambient measurements are used to validate the load model [45]. As a result, the confidence in model accuracy compared to other components is relatively low. The uncertainty in the dynamic response of the load and its impact on the analysis should also be carefully examined.
The DERs scattered throughout the distribution network are often modeled as a single power resource through aggregation. The model structure that incorporates the aggregated load and DER, known as the composite load model with a DER (CLMDER) [24], is presented in Fig. 4. In Fig. 4, the DERA [22], [46] is a simplified model of the plant-level generator introduced in Fig. 3. Note again that the main purpose of using this model is to analyze the overall system dynamic performance for feasible operating scenarios at the transmission level. More accurate analysis requires a detailed representation of the loads and DERs of the distribution system. This section provided a brief overview of the latest modeling outcomes in power system planning and operations studies at the transmission system level, focusing on renewable energy resources such as plant-level renewables and BTMG scattered throughout the system. Load models have been historically known to be the most uncertain among various dynamic models of power systems, and accurately modeling load dynamics from the transmission system's perspective has become even more challenging due to the rapidly increased number of DERs across the entire grid. The CLMDER introduced in Fig. 4 allows for a more realistic representation of the aggregated features of net loads (i.e., load and DER). In the following, two typical analyses of power system dynamic performance will be addressed using this model, along with a discussion of how the uncertainty of the CLMDER affects the analysis results.

III. G R I D I M P A C T S T U D I E S O F L A R G E -S C A L E D E R s A N D M O D E L P A R A M E T R I C U N C E R T A I N T Y
This section discusses the grid impact of large-scale DERs in the Korea power system by incorporating CLMDERs from Fig. 4 into real planning cases. A total of 1.43 GW of BTMGs are modeled in four regions (see Fig. 5) using CLMDER with widely used parameter values [43], [46]. It is worth mentioning that the typical parameter values are currently used in Korea due to the early stage of case studies utilizing this model. Also, along with this, ongoing research is currently in process to verify the model parameters based on measurements, which will help obtain more reliable parameters that can better represent the distinctive features of the Korean system. The remaining loads are modeled using the complex load (CLOD) model [46], a simplified model for representing load dynamics. The analysis results obtained using the initially defined model parameters are referred to as the base case. Dynamic simulations were performed for two contingency scenarios: a transmission line failure and an unplanned generator trip. In addition, we evaluated the uncertainty impact of CLMDER parameters, considering the extent of aggregation in this model. The cases designed to evaluate parametric uncertainty are referred to as lower and upper cases, considering the ranges of crucial parameters.  demand. It is worth noting that this scenario is comprehensively evaluated in operations and planning studies carried out in Korea due to its potential impact on the entire power system. In addition, note that the base case represents the simulation using the initial model parameters, whereas the lower and upper cases illustrate how the parametric uncertainty of the CLMDER can significantly affect the study results, leading to delayed and faster voltage recovery compared to the base case. Analytical insights into these behavioral changes are provided in Section III-A2.

A. Case Study for Transmission Line Failure
2) Evaluating Perturbation Impact of Uncertain Parameters: Note again that power system dynamics analysis investigates how the system responds after certain events occur at a specific operating point, which means that the constantly changing load and its dynamic model are assumed to be specific values at the given operating point. Various model parameters have been utilized for each time period or season to account for the variability in load characteristics [43]. Nevertheless, the uncertainty associated with the load model and DER is relatively high compared to other power system models, primarily due to the abstract nature of the modeling philosophy seen in Fig. 4 and the continuously changing load characteristics. This inherent uncertainty may lead to a completely different conclusion from the power system analysis, as manifested in Figs. 6 and 7.
Many studies have thus been conducted to deal with this uncertainty [47], [48], [49], [50], [51], [52], [53]. Probabilistic-based analyses (e.g., [48] and [52]) can examine the influence of uncertain parameters but require considerable computational costs for the combined effects of multiple parameters, especially for large-scale power systems. Conversely, trajectory sensitivity-based analyses (e.g., [49], [51], [54], and [50]) are advantageous in conducting worst case analysis efficiently and thus have been widely utilized in practical power system analysis to consider various scenarios. Note that trajectory sensitivity is one of the practical methods that investigates how the perturbation of specific parameters affects the analysis  results [47]. The trajectory sensitivity can be derived by partially differentiating DAEs (1) and (2) with respect to the parameter of interest as follows: λi represents one of the uncertain parameters in the CLMDER. The derived sensitivity can identify how the perturbation of this uncertain parameter affects the analyzed result (e.g., the base case in Figs. 6 and 7). After deriving the trajectory sensitivities for the interested parameters, we can review the perturbation impact of multiple parameters on the initially obtained result of power system dynamics (e.g., bus voltage responses) by the following approximation: Here, λ represents the uncertain parameters, J denotes a Jacobian matrix comprising the derived sensitivities of the trajectory, and V initial and V perturbed refer to the voltage responses initially obtained and the calculated voltage responses considering the variability of multiple parameters, respectively. Further details regarding (5) can be found in [47]. It should be noted that although this article uses the bus voltage as an illustration for (5), the methodology can be applied similar to all other responses of interest. The sensitivity-based approximation (5) can only be useful when the higher order terms of the Taylor series expansion of the trajectory are negligible. In cases where the system is under high stress (i.e., the trajectory approaches the boundary of the region of attraction), the performance of the approximation (5) cannot be guaranteed, as discussed in [47] and [55]. Furthermore, the previous study [51] demonstrated that parametric perturbations exceeding a certain threshold could result in completely different power system responses due to the nonsmooth behavior of power system models (e.g., load models). Therefore, the validity of trajectory sensitivity and its applications may be limited to a narrow range near the initial point where the sensitivity is calculated, especially for highly nonlinear systems.
This limitation in sensitivity-based analysis was discussed in previous studies [49], [50], [51], [54], and recent studies [49], [50] have proposed methods to address these limitations. The method developed in [49] and [50] combines the screening process and the additional detailed analysis to identify the most significant impacts of multiple uncertain parameters. The screening process can identify the anticipated case with the most significant impact (i.e., the worst case perturbation) by utilizing the trajectory sensitivity [e.g., (3) and (4)] and the sensitivity-based approximation [see (5)]. Then, the additional detailed analysis conducts a power system simulation using the perturbed parameter set identified by the screening process. It is worth mentioning that this method considered the situations where the sensitivity-based approximation could lead to misleading conclusions under the highly nonlinear parametric influence. Therefore, the sensitivity analysis was used solely to identify the perturbation combination with the most significant impact. Subsequently, the additional detailed analysis was conducted to guarantee dependable analysis results, regardless of the presence of nonlinear parametric behavior. The case study in this section also employs this method but only shows the outcomes of the additional detailed analysis without the screening process results. The main aim of this section is to visually demonstrate the impact of various uncertain parameters on the initially analyzed result, and detailed information on this method can be found in [49] and [50].
This article assumes ±20% uncertainty for five selected crucial parameters of CLMDER in a specific region (i.e., see the marked blue area in Fig. 5). A total of 20 loads (about  Fig. 4 and the recovered DER fraction (Vrfrac) after voltage returns to the defined range. Note also that the initially defined values can be seen in Table 1. From the result of the screening analysis, we may conclude that voltage recovery is delayed as the motor ratio of each motor increases or Vrfrac decreases. Therefore, based on this information, the additional power system simulation is conducted using the perturbed parameter set that 20% increased fractions of four motors (FmA = 0.12, FmB = 0.12, FmC = 0.12, and FmD = 0.36) and the 20% decreased Vrfrac (Vrfrac = 0.64) in Table 1; the simulation result is the lower case in Figs. 6 and 7. In addition, the opposite 20% perturbation (FmA = 0.08, FmB = 0.08, FmC = 0.08, FmD = 0.24, and Vrfrac = 0.96) is considered for comparison as the upper case. Note that the fraction of the static model in Fig. 4 varies depending on the changes in the fractions of motor models because the total sum of the load component fractions remains one [44]. shows the voltage responses of the most and least affected buses, respectively, among all 20 buses assuming the uncertain model. Fig. 6(c) and (d) shows the voltage responses of the most and the moderately affected buses, respectively, in the metropolitan area. The lower case (i.e., increased motor loads) shows delayed voltage recovery, whereas the upper case (i.e., decreased motor loads) results in faster recovery than the base case. This is primarily due to increased motor loads drawing more reactive power in low-voltage operating conditions [56], consequently impeding faster voltage recovery. Note that transient voltage criteria [57] (i.e., the acceptable voltage level during voltage recovery) are generally established for power system studies. For example, the voltage magnitude should be maintained above 0.9 pu within 10 s after any disturbance; Fig. 6 shows that the voltage levels of the base and upper cases are recovered within 10 s [i.e., at 8 s in Fig. 6(a) and (c)] after the transmission line failure, thus meeting the voltage criteria. However, the lower case in Fig. 6(a) and (c) does not comply with the criteria because the voltage magnitude cannot return to the desired level at 10 s. Therefore, there is no transient voltage issue for the given scenario (i.e., 765-kV transmission line failure) from the conventional analysis (using only the single model parameter set). However, a different conclusion could be reached when parameter uncertainty is considered.

3) Perturbation Impact of Uncertain Model Parameters:
The perturbation impact of uncertain parameters may vary depending on the power system operation condition or the disturbance types. Previous research [49], [55] observed that the impact of parameter uncertainty could be more salient when severe disturbances occur, causing stability concerns. Therefore, it is advised to examine the parametric uncertainty when analyzing the highimpact (generally low-frequency) contingency events, which should certainly improve the confidence level in the studies and thus reliability of power system planning and operations.
It is also important to note that the perturbation impact may have a different effect on different quantities, even when the same uncertainty is assumed, as shown in Fig. 6(a) and (b). Furthermore, although Buses C and D are located in the metropolitan area shown in Fig. 5 (i.e., not the area assuming the uncertainty), Fig. 6 shows that the voltage response of Bus B, which is subject to assumed uncertainty, is impacted less than that of Buses C and D. In particular, Fig. 6 suggests that Bus C is affected as much as the buses (e.g., Bus A) where the uncertainty is assumed. Previous research [50] has shown that parameter uncertainty can significantly affect the analysis results of specific locations. Similar findings are confirmed in the case study of this article utilizing the CLMDER (i.e., an advanced model that combines load and DER), which also thus advises the necessity of system-wide impact studies.
When considering parametric uncertainty, the dynamic behavior of the modeled BTMGs changes, as shown in Fig. 7. As discussed in Fig. 6, the lower case exhibits a more delayed voltage recovery due to the increased motor loads, resulting in the slower recovery of active power generation from BTMGs due to the modeled control logic. The CLMDER can replicate the dynamic behavior of specific BTMGs that cannot supply active power during voltage recovery by adjusting the Vrfrac parameter and other voltage control logic-related parameters [46]. This decrease in active power output causes a reduction in the supply of reactive power from BTMGs, which ultimately impedes voltage recovery.
On the other hand, in Fig. 7(a), it is evident that the increased value of the Vrfrac parameter (i.e., Vrfrac = 0.96) in the upper case increases the recovered BTMG. Note again that the increasing Vrfrac results in an increase in power supply from the BTMGs during voltage recovery [46]. As shown in the upper case of Fig. 7(b), the BTMGs supply more reactive power and then decrease the reactive power supply after 8 s because the voltage is recovered at a predefined level, consistent with the faster voltage recovery observed in the upper case of Fig. 6. iMoreover, the system frequency recovery is also delayed in the lower case, as shown in Fig. 7(c), due to the power supply reduction from BTMGs, leading to an increase in the net load. In actual power system operation, the net load characteristics caused by the dynamics of BTMG are not visible to the system operators. Therefore, it is crucial to conduct comprehensive analyses, as shown in Figs. 6 and 7, using adequate models.  Fig. 8 presents the case study results of an unplanned generator trip of about 800 MW under the same operating conditions as the previous case study. Here again, we first focus on the results derived from the initially defined parameters (i.e., the base case of Fig. 8). It can be seen that the system frequency [see Fig. 8(a)] momentarily drops due to an imbalance in power supply and demand caused by the unplanned generation trip and then gradually recovers. During the process of frequency recovery, Fig. 8(b) presents that the BTMG trip occurs due to the modeled BTMG low-frequency trip parameter (f low ), which adversely affects the frequency recovery. As discussed, the BTMG trips lead to the net load increase seen in Fig. 8(b), but it should be stressed again that the BTMG is not visible to the power system operator; only the net load response is visible. Recently, the KPX has observed some cases of instantaneous net load increase after the unplanned generator trip [12], especially during the daytime when the penetration level of BTMG is high. These recent observations highly suggest the importance of modeling the low-frequency trip characteristics of BTMG.

2) Perturbation Impact of Uncertain Parameters:
The impact of parametric uncertainty is also reviewed for the unplanned generation trip by assuming the low-frequency trip parameter (f low ) to be uncertain because it mainly affects the instant of the BTMG trip during frequency recovery [12], [46]. The values of f low initially set for four BTMGs in Fig. 5  It is evident that the larger f low value leads to the earlier trip of BTMG at a higher frequency. Therefore, without the sensitivity analysis, we can intuitively see that the increased f low (i.e., positive perturbation) is the worst case for frequency response; in contrast, the negative perturbation leads to even better frequency recovery. The parametric perturbation of ±0.02 Hz is considered for f low , and the analyzed result can be seen in Fig. 8.
As expected, when the f low values are set to larger values (i.e., the lower case in Fig. 8, named as such because the frequency drops faster than base case), the minimum system frequency is reached faster because the BTMG trip occurs at a higher frequency value. Consequently, the rate of change of the frequency (RoCoF) also becomes larger, which may lead to further cascade tripping of other generators. On the other hand, in the upper case of Fig. 8, some parts of the BTMGs are tripped, resulting in a smaller load increase compared to that of the base and lower cases and, consequently, less frequency drop. Fig. 8 suggests that the analyzed result for the power system's frequency response can vary depending on the uncertainty of the BTMG trip. However, as discussed, we cannot precisely determine how much BTMG has tripped in the real system, but we can confirm that the total net load increases seen in Fig. 8. In practice, modeling the frequency trip characteristics of BTMG can be achieved [12] by utilizing frequency and net load data acquired when a specific disturbance occurs in the actual power system. Therefore, continuous efforts are necessary to ensure the reliability of BTMG trip modeling through consistent management of the disturbance data recorded during power system operations. In addition, as shown in Fig. 8, we should evaluate the inherent parametric uncertainty to increase the reliability of the power system analysis.
In this section, it has been demonstrated that the aggregated dynamics of BTMGs (which are unknown to the power system operator in reality) can be evaluated through model-based simulation studies by using the CLMDER. It is noteworthy that the CLMDER contains over 170 model parameters, which presents a significant challenge to its ease of use. Nonetheless, it is necessary to employ this model in modern power system analyses, especially those with extensively expanded BTMGs. Therefore, it is recommended first to commence the initial case study using the general parameters guided in [43] and [46]. Then, once the actual net load data are available, the model accuracy should be checked and validated to employ a more reliable model in the subsequent studies. In addition, it should be noted that the uncertainty of the net load model is expected to increase with the increasing penetration of BTMG due to the inherent limitations of aggregated modeling in practice. More comprehensive and reliable methods are necessary to improve the existing simulation studies that rely on a fixed single parameter set; one example was discussed in this section. Further studies are needed to make these methods practical for large-scale power system studies, in addition to ongoing efforts to minimize the discrepancy between the modeled power system and the real-world system.

IV. S T O C H A S T I C S T A B I L I T Y A N A L Y S I S O F V R E -R I C H P O W E R S Y S T E M S
In order to embrace variability and intermittency of renewable energy sources, stochastic approaches have been adopted for power system simulation studies. The stochastic nature of solar irradiance and wind data can be captured by Gaussian [58], [59], [60], [61] and Weibull distributions [62], [63], [64], [65], [66], [67] as shown in Fig. 9 for stochastic analysis of power flows and system stability. These probabilistic approaches are expected to be further deployed along with the advancement of computing performance [68].
As the employment of renewables increases, the inverter dynamics on top of the inherent variability have shown a growing impact on power system stability in terms of transient stability, voltage stability, and small-signal stability [70], [71]. Deterministic methods for stability analysis are typically structured such that dynamic performances under predetermined operating conditions are evaluated.
The results indicate operating conditions or cases more stable or less stable. However, stochastic methods incorporate the stochastic nature of renewable generation for power system stability. For example, the stochastic method provides stable operating conditions of the power system at a specific time with a certain confidence level subject to change with the weather condition. Planners or operators may then investigate their findings and conclude whether the system is stable or not [70] at their discretion.
Since the VRE, lacking dispatchability, cannot follow the generation schedule, the deterministic method has fundamental limitations on its use for the VRE-rich power system studies. From an industry point of view, one of the most significant challenges in applying the probabilistic method is the long simulation time due to the high computation burden. However, the recent advancement in computer performance has expanded the applications of stochastic approaches.

A. Stochastic Voltage Stability
Voltage stability is defined as the ability to maintain the voltage when disturbances occur in the power system [17], [74]. Short-term voltage stability issues may occur due to dynamic characteristics of fast-acting load components, such as electronically controlled loads, induction motors, and HVdc systems. Voltage stability issues in the long-term may occur from slow-acting equipment such as tap changers in transformers, thermostatically controlled loads, and current limiters in generators [17], [74].
We have successfully addressed the voltage stability by computing the maximum loadability using active powervoltage (PV) and reactive power-voltage (QV) curves under predetermined operating scenarios. As observed in [72], random nature of the active power or reactive power of the load may cause voltage instability, which advises the importance of incorporating the variability of renewable energy as well for voltage stability assessment. The Jacobian matrix has been used to investigate how a small change in voltage magnitude and angle affects the change of active and reactive power. The study in [73] and [74] analyzes the eigenvalues and observes that the power system becomes unstable due to the stochastic nature of the load, although the operating point is stable in the Jacobian matrix. When dealing with the random characteristics of loads, stochastic differential equations are incorporated into the algebraic equations of the power system model, which can lead to singularity instability from a mathematical perspective according to [72].
As differential algebraic approaches may encounter numerical singularity issues under incorrect models of certain power system conditions, which is not the physical power system phenomenon, verification of diverse operating conditions is important to investigate the possible case of singularity instability. A stochastic model of DAE for power system model can thus be employed to analyze and interpret diverse system conditions. Interestingly, a Vol. 111 deterministic method may indicate voltage stability, which is not always the case in a stochastic way.
Sustained voltage oscillations viewed also as the dynamic voltage stability in a power system can arise from various causes associated with weak grid dynamics and stochastic disturbances [73]. These oscillations are expected to be well-damped, but it is not unusual to observe poorly damped or even growing oscillations in modern power systems. These problematic oscillations may be categorized in terms of the sources as follows [73], [75]. 1) Poorly damped or even negatively damped growing oscillations are generally associated with the system dynamics under certain operating conditions, which are in the category of natural oscillations. 2) Forced oscillations associated with periodic or random input sources (i.e., disturbances) to the systems. The oscillation sources could be cyclic loads (e.g., aluminum smelting), malfunctioning controls, HVdc controls, wind power plant voltage control actions, broken valve on thermal units, and so on.
Dynamic voltage stability concerns, particularly due to growing VRE or IBR penetration, need to be well addressed in the grid planning and operations and adequate mitigation or corrective measures should be taken in time. Thus, methods [76], [77] and emerging stochastic small-signal stability detailed in Section IV-C should be beneficial as the computing capability advances.

B. Stochastic Transient Stability
Transient stability is defined as the ability of a generator to maintain synchronism with the power system following a large disturbance, which is related to the rotor angle and angular speed of the synchronous generator [17], [18], [74]. Transient stability depends on the prefault system conditions, the severity of the fault, and the fault clearance method [17], [74], [78].
One of the commonly used indicators to describe transient stability is the transient stability index (TSI) [79], [80], [81], [82]. TSI is an indicator used to assess the transient stability of power systems, providing a quantitative evaluation of how close the system is to instability during transient events [82]. Regarding this index, the higher TSI value can be characterized as the stabler system, and the lower value can be interpreted as the relatively less stable system. In consideration of the influence of the stochastic nature of renewable generation, several studies probabilistically interpret the stability of the power system using probability density function (pdf), cumulative density function (cdf), or TSI [79], [80], [81]. By means of this stochastic method, power system stability can be analyzed by determining whether TSI is high or low [81].
In addition to modeling the stochastic nature of renewable generation, various power system conditions can also be designed in a stochastic fashion [83], such as the uncertainty of load shedding events, and system contingency, including fault occurrence, fault impedance, fault removal time, and reclosing operation. The probabilistic-based process modeling is expected to be used more in the future.
The probabilistic generation of renewable energy has an impact on the operating point of conventional generators, leading to variations in voltage phase angle. Recently, Korea Electric Power Corporation (KEPCO), Naju, South Korea, and Korea Electric Power Research Institute (KEPRI), Daejeon, South Korea, have developed a tool for probabilistic analysis for power system planning, which is expected to be widely used in power system analysis [81]. In the future, we foresee that the use of this probabilistic power system stability analysis will be further expanded and utilized in the actual power system operation.

C. Small Disturbance Rotor Angle Stability
Small-signal rotor angle stability has been traditionally explained as the interactions among generators or groups of generators, which is defined as a global problem. In this regard, the global problem of frequency stability is not limited to a single generator or specific regions but rather encompasses problems that affect the entire power grid. Local oscillations occur when a single generator swings against another generator or the rest of the system, whereas interarea oscillations occur when a group of generators swings against another group of generators [17], [18], [74]. Stochastic small-signal stability studies are conducted to identify eigenvalues, damping ratio, and oscillation frequencies [62], [69], [84], [85], [86], [87]. Power systems are intertwined with diverse and complex system components and control parameters, which can be observed with the critical modes indicating the stability of the power system. In order to accommodate large amounts of VRE in a power system, several studies conducted the small-signal rotor angle stability assessment considering the inherent stochastic variability of VRE [69].
In this study, the IEEE 14-generator model is used as a test bed for probabilistic stability analysis, incorporating the Australian future network plan with renewable generation. This research studies how eigenvalue changes when system inertia of Area4 in Fig. 10 changes as in Table 2 along with the change in stochastic renewable generation.
The result represents the probabilistically distributed eigenvalue for the change in the inertial energy, as shown in Fig. 11. Furthermore, in Fig. 11, the power system inertia constant was changed due to the turned-off synchronous generator from the reason of an increasing renewable generation. Fig. 11(b) shows that the damping ratio of eigenvalues decreases due to reducing system inertia from Case1 to Case3, which is less than a critical value [69]. Some of the studies also carry out a stochastic eigenvalue analysis to design the controllers of HVdc or inverter-based power sources [73], [74].
In conclusion, the likelihood of instability of the power systems generally increases as the share of renewable energy in the system grows resulting from reducing inertial energy in the system. The impact of decreasing inertial energy on frequency stability is already widely acknowledged; however, a recent study has revealed that inertial energy also significantly affects the stability of small-signal oscillations [69].

A. Features and Needs of EMT Modeling and Simulation
Power system operators and planners have long relied on rms modeling and related tools for transient stability analysis (TSA), as they can represent the crucial dynamics of such systems, such as synchronous generators, networks, and loads (although to a limited extent). These approaches also provide dependable computing environments for various types of studies that help ensure the stable operation of conventional large power systems. However, due to the limitations of the positive-sequence rms model in capturing the fast and sophisticated dynamics of IBRs and their control interactions, there have been notable efforts to use EMT modeling and simulation tools for TSA. Depending on rms modeling for analysis, the results could therefore lead to misleading conclusions. Fig. 12 shows the different timeframes of various power system components and their operations [74], [88]. The frequency range of IBR dynamics is much higher than that of electromechanical synchronous machines, which have traditionally been the primary focus in the power system. This difference in bandwidth highlights the need for analysis tools that operate on a microsecond timescale. The fast control behaviors of IBRs, such as phase-locked loop (PLL), inner loop current control, and IBR and grid impedance interactions, contribute to emerging resonance and converter-driven stability concerns.

Fig. 13. EMT simulation result of Jeju Island transmission network for unbalanced fault. (a) Instantaneous CB current. (b) Bus frequency. (c) Instantaneous bus voltage. (d) RMS bus voltage. (e) Instantaneous dc voltage of LCC HVdc. (f) Alpha angle of LCC HVdc rectifier.
Furthermore, given that the switching characteristics in power electronics front ends are crucial for examining the IBR-related stability issues on the power grid, it is necessary to accurately model and control their fast dynamics and switching events to ensure the stable and reliable power system operation. To represent high-frequency dynamics [29], [89], the EMT simulation is critical for accurately modeling these effects and simulating various scenarios to study power system performance and stability. Fig. 13 presents an EMT simulation example of the power grid on Jeju Island. This simulation was conducted using EMT models to support the postmortem analysis of an actual event that involved a single-line-to-ground fault on a transmission line, followed by a line reclose and trip. It should be noted that an rms case study cannot reproduce this event due to limitations in analyzing unbalanced disturbances. Fig. 13(a) shows the fault current flow through the circuit breaker. During the event, a construction crane accidentally approached the transmission line, resulting in an A-phase-to-ground fault. Fig. 13(b) shows the system frequency, which drops to its nadir in less than 2 s due to Jeju Island's nearly 80% power electronics penetration. In addition, a large amount of BTMG tripped out from the event, increasing the excursion in the frequency curve.
The voltage profiles of the faulted location can be seen in Fig. 13(c) and (d), which are instantaneous and rms, respectively, and provided to analyze unbalanced fault cases. Fig. 13(e) and (f) shows the momentarily power reversal of the Jeju-Mainland HVdc poles during the event. The low voltage leads to a high alpha angle and can cause temporarily inverted dc voltage during the unbalanced condition. This example demonstrates that the EMT analysis can effectively examine power systems with a high penetration of IBRs, such as Jeju Island.
Although EMT studies can provide high accuracy and usefulness, the unavailability of EMT model parameters can be a limiting factor. For example, modeling every IBR in a renewable farm using EMT can be inefficient and impractical, and the parameters may be unavailable. In such cases, aggregating the IBRs with comprehensive information on their behavior and characteristics is necessary to simulate the EMT analysis effectively. Therefore, it is essential to utilize validated EMT models that incorporate the behavior and features of IBRs. By conducting EMT studies with validated IBR EMT models, power system planners and operators can acquire valuable insights while ensuring the stability and reliability of the power system with a high penetration of IBRs.

B. EMT Models and Studies With Grid-Forming Inverters
The mechanical inertia inherent in synchronous generators plays a critical role in resisting changes in grid frequency. An increase in IBR penetration tends to exacerbate frequency instability and weaken voltage strength, resulting in higher RoCoF and lower short-circuit levels, and the ability to maintain voltage waveform [90], [91], [92]. Most existing IBRs are grid-following (GFL) types that act as current sources, displacing the existing voltage sources (e.g., synchronous generators), making grid synchronization of IBRs more challenging. In contrast, grid-forming (GFM) inverters can cope with potential stability issues due to high IBR penetration. GFM inverters can emulate the inertial effect by independently establishing an internal voltage reference. Given that GFM is envisioned as a key solution to the potential stability issues of the future power grid, numerous research works exist on the control and design of GFM inverters [90], [93], [94], [95].
Although the functional requirements for GFM inverters are not yet clear, it is fortunate that discussion groups are moving forward to establish grid codes and standards for them. One of the major requirements is sufficient fault current supply. If the voltage is reduced during a faulted condition, the controller should recover the voltage as soon as possible after the disturbance. In addition, their behaviors must be coordinated with the existing protection devices to ensure a seamless protection system in the transmission network. For example, as GFM inverters gradually replace synchronous generators, the distance protection zone grading around the plant must be recalculated due to the reduced fault current level. Even though system operators periodically do this job, the penetration of GFM voltage sources will dramatically change the map of shortcircuit current, and there might be a need for a different zone grading strategy than before.
The EMT simulation will be essential in verifying the GFM functionalities and performance under fault conditions of the grid. As shown in Fig. 14, a simplified ac/dc system example is simulated to test different GFM controls, with a single GFM source operating in parallel with the bulk power system. Fig. 15 shows the results of the fast frequency regulation performance of the GFM inverter controls [96]. Fig. 15 shows the matching scheme in green [97], the droop in blue, and the virtual synchronous machine (VSM) in red, with each GFM frequency control scheme exhibiting different regulation performances. Additional research and development is required to identify the optimal strategy for ensuring system stability. Manufacturers must also consider associated responses, such as the primary energy source, dc-side dynamics, and protection coordination with multiple devices, in addition to the GFM inverter model.

1) High-Performance Computing EMT Simulation:
Although the EMT model accurately represents the actual plant, its simulation time has been viewed as a significant drawback. The complexity of the EMT simulation is significantly higher than local-area or equipment-focused simulations, particularly when studying transmission-level networks with numerous IBRs. Therefore, it is very challenging to manage the computational burden of EMT simulations for large-scale power systems using only single-core processing.
The parallel computing platform has been widely adopted by many system operators, such as ERCOT and AEMO, to accomplish simulations faster. Based on the applicable time step of the EMT simulation, these platforms typically divide the network into segments using Dommel's traveling-wave transmission line model for sufficiently long transmission lines [98]. This approach enables the simulation to utilize segmented networks on separate CPU cores, significantly reducing the required time for the EMT simulation. Although network segmentation for EMT simulation is not a new idea and has been utilized in real-time simulation systems that require dedicated hardware, commercial nonreal-time EMT programs employ high-performance computing (HPC) technology that leverages multiple CPU cores and high-speed communication to  expedite the simulation and increase the application of the EMT study.
The execution time of EMT simulation can be significantly reduced by dividing and distributing the simulation tasks to multiple processing units through a low-latency communication channel, such as shared memory. Fig. 16 shows an example of the application of HPC technology for EMT simulation. In Section V-A, the EMT simulation study took approximately 51 min to complete 10 s of simulation time under conventional single-core execution. In contrast, the simulation execution time was reduced to 16 min when using parallel computing with eight cores. The study case includes 77 buses, 30 transmission lines, nine synchronous generators, 13 wind turbine plants, ten PV plants, two HVdc links, and two STATCOM stations.
Although the computational benefits may decrease as the number of cores increases due to communication latency, the simulation time primarily depends on model complexity. For instance, when a test power system with 4320 components is divided into 32 subprojects, it takes around 17 h to complete a 10-s simulation. However, if the total number of components is reduced to 1216, the total execution time is decreased to approximately 2 h.
2) Initializer to Expedite the EMT Simulation: Apart from utilizing HPC technology, users should establish initial conditions for network components, including transformer tap ratios, bus voltages and angles, active and reactive load amounts, operating points of the flexible ac transmission system (FACTS) devices, and network topology. The extended ZIP load model (constant impedance (Z), constant current (I), and constant power (P) model) incorporating the frequency dependency, is employed to represent the load dynamics. More sophisticated load model should help improve the study accuracy. However, given the study objectives and computing constraints, more computing resources are allocated to represent the power electronic components, including various control functions, and to analyze their grid impact on time. For effective use of advanced EMT simulation, users must establish the initial conditions of the simulation case prior to executing the EMT simulation. As a result, a useroriented initializing application is required to successfully transfer the power-flow solution and network topology into the initial conditions of the advanced EMT simulation platform.
3) Future of the HPC EMT Simulation: HPC EMT simulation is practical in analyzing transmission network disturbances, particularly for high IBR penetration systems. As shown in Fig. 13, voltage and frequency excursions are accurately observed for unbalanced faults based on the well-designed EMT model. The HPC EMT simulation can also adequately represent the dynamic behaviors of the FACTS and IBRs, as well as the tripping or momentary cessation of legacy-type inverters. Therefore, the HPC EMT simulation should become a more crucial and handy tool for system operators and planners to address weak grid and low inertia issues due to the rapid and widespread growth of IBRs. 4) Co-Simulation Platform: Sustained efforts from academic and industry researchers have been conducted to develop a co-simulation platform that integrates both a TSA program and an EMT program [99], [100], [101], [102], [103], [104], [105], [106], [107], [108], [109], [110], [111], [112], which has recently been incorporated into commercial EMT products [113], [114], [115]. The co-simulation platform offers the benefits of both programs so that it can deliver high-speed simulation performance using a real-time EMT simulator or the EMT program with high-performance parallel computing. We can summarize the advantages of the co-simulation platform as follows.

1) Fast and efficient simulation:
The time and effort necessary to establish and study systems are reduced compared to full EMT simulation. This includes not only modeling and configuring the study system but also circuit design, system compiling, downloading, and data monitoring. 2) Cost saving: The co-simulation platform also saves on installation and operational costs by significantly reducing the number of real-time simulator processors or HPC units. 3) Simulation accuracy: The co-simulation platform can provide accurate simulation results comparable to full EMT simulation, particularly when the simulation contingency is located within the internal area modeled with the detailed EMT [31].

4)
Available hardware validation: Validation of power system equipment using hardware-in-the-loop (HIL) testing (for only the real-time simulator applied co-simulation platform) has demonstrated its efficacy in providing accurate results under realistic, flexible, and reproducible conditions. In addition, the cosimulation's flexibility supports the system-level testing of controller HIL (CHIL) and power HIL (PHIL) when actual controllers and equipment are tested.
The co-simulation platform is a useful tool for simulating the transition to a cleaner and sustainable power system with increasing integration of IBRs. However, it has a limitation in that it cannot ensure the accuracy of fault studies near interfacing buses. To resolve this issue, an option could be to enlarge the internal area (EMT simulation region) surrounding the fault location within the co-simulation platform. Furthermore, the electrical distance is utilized to establish an appropriate EMT simulation region from the fault location.

5) Experience With Co-Simulation:
Recent efforts have been undertaken to overcome the limitations of co-simulation by establishing the purpose, scope, and process of co-simulation studies in the electric utility industry in South Korea [31], [32]. Fig. 17 shows the co-simulation platform in operation at KEPRI, which has the capability of CHIL and PHIL to support system-wide impact studies. KEPRI conducted a dynamic performance test using a large-scale real-time simulator and a replica controller to assess the behavior and performance of newly installed power electronics equipment in the Korea power system. This test evaluates the equipment's response under various operating conditions and disturbances before it is integrated into the operational power grid. In 2018, KEPRI tested actual static var compensator (SVC) replica controllers for dynamic performance studies using the co-simulation platform. Despite the limitations of co-simulation, it enables the flexibility to conduct comprehensive simulation studies, performing a significant role in analyzing the dynamic performance of future power systems by bridging the gap between TSA models and EMT models.

VI. C O N C L U S I O N
This article addressed selected modeling and analysis methods to effectively capture the growing uncertainty, variability, and fast dynamics of modern power systems evolving with high penetrations of power electronic-interfaced and distributed VRE. Focusing on bulk power system planning and operations, these efforts aim to advance the existing practices and tools for assessing the potential security concerns and identifying practical solutions.
The state-of-the-art positive-sequence rms models, often aggregated for large-scale power system studies, were briefly explained, and recent efforts to supplement limitations of the existing rms model-based approaches were introduced. In particular, the model parametric uncertainty and its consequential impact on the dynamic performance were investigated by using a sensitivity-based analysis for Korea power system. This article also discussed limitations of the deterministic approaches to the VRE-rich system security assessment, typical industry practices using the rms models above. It then introduced a stochastic study framework, incorporating the probabilistic distributions of the VRE generations. The stochastic method, as demonstrated through practical studies, effectively handles expanded operating conditions and uncertainties to ensure the system stability.
In addition, this article featured the importance of EMT-based modeling to adequately study the sophisticated dynamic performance of the IBRs and various types of interactions as manifest in a series of disturbance reports from the global utilities. Overcoming the computational limitations and modeling complexity, we introduced research efforts to incorporate the EMT dynamics of the VRE for large integrated power system studies through the parallel computing with HPC and co-simulation. Continued and extensive efforts should be exerted to fully understand the power system dynamics in transition and advance (and develop new) analytical as well as computational methods and supporting tools in the planning and operations.