Power System Stability With a High Penetration of Inverter-Based Resources

Inverter-based resources (IBRs) possess dynamics that are significantly different from those of synchronous-generator-based sources and as IBR penetrations grow the dynamics of power systems are changing. This article discusses the characteristics of the new dynamics and examines how they can be accommodated into the long-standing categorizations of power system stability in terms of angle, frequency, and voltage stability. It is argued that inverters are causing the frequency range over which angle, frequency, and voltage dynamics act to extend such that the previously partitioned categories are now coupled and further coupled to new electromagnetic modes. While grid-forming (GFM) inverters share many characteristics with generators, grid-following (GFL) inverters are different. This is explored in terms of similarities and differences in synchronization, inertia, and voltage control. The concept of duality is used to unify the synchronization principles of GFM and GFL inverters and, thus, established the generalized angle dynamics. This enables the analytical study of GFM-GFL interaction, which is particularly important to guide the placement of GFM apparatuses and is even more important if GFM inverters are allowed to fall back to the GFL mode during faults to avoid oversizing to support short-term overload. Both GFL and GFM inverters contribute to voltage strength but with marked differences, which implies new features of voltage stability. Several directions for further research are identified, including: 1) extensions of nonlinear stability analysis to accommodate new inverter behaviors with cross-coupled time frames; 2) establishment of spatial–temporal indices of system strength and stability margin to guide the provision of new stability services; and 3) data-driven approaches to combat increased system complexity and confidentiality of inverter models.

investigations are now needed to establish the root cause [15]. The stability challenges may become more significant in the future when IBRs take a dominant role in a carbon-neutral power system.
One of the main difficulties for stability studies of inverter-based grids is the lack of intrinsic behaviors and standard models. An inverter has little inherent functionality resulting from its physical form, and almost all of its features are defined by its control algorithms. The control functions are arranged in a hierarchy. At the bottom, and fastest, is the modulation of the semiconductor switching to synthesize a voltage, which is further controlled to achieve the top-level objectives. At the top are two broad choices. In the first, a regulated power is exported based on the resource available or market position. To achieve this, the inverter is synchronized to the local grid voltage, and a current vector is injected. This arrangement is generally referred to as a GFL inverter [16], [17]. The GFL inverter is the dominant format of inverters in power systems to date. The second choice is to control the inverter as a voltage source from which power can be drawn according to system conditions. This is known as a grid-forming (GFM) inverter [18], [19]. The voltage is a free-running oscillator except that its frequency is made to reduce in proportion to the power drawn from the source. This frequency droop gives the voltage source the ability to lock to a grid and run synchronously. GFM can be compared with conventional generators in that they synchronize and regulate power in similar ways [20], [21] although GFM inverters may behave very differently in short circuits due to the lack of overcurrent capacity. Thus, the main differences between GFL and GFM are achieved by changing the control arrangements to achieve a current source or voltage source behavior, but this brings with it changes in the method of synchronization and changes in priority of grid or resource conditions in setting the power delivery.
Although GFL inverters are dominant at present in commercial inverters, GFM inverters feature prominently in research and are seen as the solution to the stability issues, IBRs. The difference between GFM and GFL is mainly in control but also in the way power is drawn from the resource, which will affect revenue in the energy market and a need for revenue for GFM services. In principle, GFM capability needs not to add significant equipment cost, but this is no longer true if GFM functions are to persist during faults. Oversizing of power semiconductors (along with the heat sinks) and energy storage components (to support the inertia response) will be necessary to provide the expected fault response from GFM inverters, and this will add cost. Bearing in mind that the faults are major triggers of instability, the stabilization functions of GFM inverters may be compromised without enhanced fault response.
The revenue or equipment cost implications of providing GFM capability in inverters prompts a debate about the proportion of GFM needed in a system and turns attention back to what can be achieved with GFL inverters. It is demonstrated in [21]- [23] that it is possible to run a power system with nearly or exactly 100% GFL inverters. This may sound radical and perhaps not practical in normal circumstances, but it allows for GFM inverters to fall back to the GFL mode during faults without compromising stability and, thus, may solve the dilemma of whether or not to enhance fault response at added cost. As a result, GFL and GFM inverters may coexist in a future grid but with different proportions in different scenarios. This calls for a unified model for stability studies that are able to represent GFL and GFM inverters, as well as synchronous generators.
Power system stability in essence is a single problem, but it is usually intractable to address all the complexities in a single investigation. The classic framework for power system stability has been established, which decomposes the stability problem into subproblems. The emergence of new stability phenomena associated with IBRs has prompted the classic framework to be revisited and extended [24]. This article will further explore the extension of the classic framework, paying special attention to accommodating the new characteristics of both GFL and GFM IBRs. We highlight the extension of the timescale over which the classic categories of dynamics act and the increased coupling between previous distinct phenomena. Comparisons will be made between the way GFL and GFM converters contribute to angle, frequency, and voltage stability through new concepts, such as generalized angle dynamics. The existing stability analysis methods will be reviewed to examine if they are sufficient to address the new characteristics of IBRs. Several unresolved issues are summarized, and potential future research opportunities are envisaged in the end.

II. C L A S S I C S T A B I L I T Y F R A M E W O R K
The investigation of power system stability has been an evolutionary process combining analytical thinking and practical experiences. Events with poor stability outcomes that have led to widespread loss of service and severe social consequences have spurred investigations of instability through international cooperation. Definitions and classifications have been produced, for example, by the IEEE/CIGRE Joint Task Force in 2004 [25], which are widely accepted and are considered as the classic framework in this article. This classic framework defines several aspects of stability in order to decompose an intractable stability problem into a series of tractable subproblems with archetypes and insights provided for each category. This categorization and its forerunners have served the sector well, but, as the displacement of generators by inverters gathers pace, and system dynamics undergo profound changes, it needs to be kept under review. One such review is by a further IEEE working group [24] that reviewed the literature on ways in which inverters, particularly GFL inverters, have modified dynamics and caused new inverter-driven stability concerns to emerge. The challenge now facing the industry is to continue to evolve our categorization and analysis of the stability of power systems as we move from generator-based systems accommodating substantial numbers of IBR to systems where IBR dominates the determination of voltage, frequency, and synchronization. In this section, we briefly review the classic stability framework of generator-based grids before introducing in the following section views of stability from the starting point of inverter characteristics in both GFM and GFL forms.
The formal definition of power system stability is given in [25] as "the ability of an electric power system, for a given initial operating condition, to regain a state of operating equilibrium after being subjected to a physical disturbance, with most system variables bounded so that practically the entire system remains intact." This definition contains two key aspects, that is, the system variables should: 1) converge to equilibrium and 2) be bounded within thresholds. These two aspects are related but have different emphases in different problems. For example, most grid standards set limits on the normal operating frequency, the lowest frequency (nadir), and the rate of change of frequency (RoCoF) [26]. A frequency perturbation that is stable in the sense of convergence (i.e., settles within operating limits in the long run) may, however, be viewed as failing a stability test if the initial response is out of bounds for RoCoF or nadir. On the other hand, angle stability is viewed somewhat differently, and the emphasis is on convergence within the first swing regardless of the extent of the initial angle excursion. In this article, we focus on the convergence aspect of stability because: 1) convergence is the necessary condition of being bounded and 2) convergence is a theoretic property of the system, whereas the bounding thresholds are often engineering tradeoffs.

A. Categorization
As illustrated in Fig. 1, the classic framework uses three dimensions to categorize power system dynamics: the dominant states, the timescale, and the disturbance level. The dominant state dimension contains three types of state variables: voltage, angle, and frequency.
These state variables are explicit and ubiquitous in power systems, and each represents an important aspect of system performance, namely, power quality (voltage), synchronization (angle), and balancing (frequency). The variation of these state variables and the associated control actions in a classic power system are much slower than the synchronous frequency, so the electromagnetic (EM) dynamics near or above the synchronous frequency (e.g., the flux decay in inductors and the resonance of shunt or series capacitors) are neglected. The implicit dynamics of controllers and apparatuses are reflected in the three explicit states. For example, the behavior of a power system stabilizer (PSS) is intended to reshape the synchronizing and damping powers of a synchronous generator and is reflected in angle dynamics.
The classic framework decomposes the power system stability into largely independent subproblems, marked by the colored boxes in Fig. 1. Such a decomposition is made possible by two natural decoupling mechanisms in conventional power systems: 1) the decoupling of active power (associated with angle and frequency) from reactive power (associated with voltage) on inductive transmission lines; 2) the decoupling of long-term (of the order of minutes) from short-term (of the order of seconds) dynamics.
For example, frequency regulation is implemented by speed governors (primary regulation of the order from seconds to tens of seconds) and thermal processes (secondary regulation of the order of tens of minutes), which is much slower than angle swing (of the order from 10 to 20 s). As a result, the angle synchronization is assumed "instantaneous" when considering frequency regulation and, because this is a synchronized set of synchronous generators, can be represented collectively in frequency stability studies as a single mass associated with the center of inertia (CoI). 1 There are exceptions where the time frames of angle swing and frequency control do overlap, e.g., due to fast speed governors reacting within seconds or in slow interarea swing modes in very large grids. In such cases, the frequency control may possibly provide damping to swing dynamics [27], [28]. A further distinction is made between small-and large-disturbance stabilities in the classic framework. Mathematically, a system that is large-disturbance stable is necessarily small-disturbance stable. The introduction of small-disturbance stability is rather a compromise due to the fact that there are few analytical techniques available to solve large-disturbance (nonlinear and possibly time-varying) dynamics, an issue to be discussed in Section IV.

B. Establishment of Indices
To assist in the understanding of stability and aid in planning the operation of systems, stability indices have been established for each category of stability to indicate if the system is adequately configured. These indices are marked in Fig. 1 and listed in the following [29]. 1) Inertia: Index for frequency stability. This is a global (whole-system) index representing the total inertia of all generators. 2) Synchronization and damping coefficients: Indices for small-disturbance angle stability. Synchronization and damping power coefficients both need to be positive at the swing frequency to make an angle swing stable. These can either be local indices (scalars) for a single generator or global indices (matrices) for the whole grid. 3) Critical energy: Index for large-disturbance angle stability. This is a global index reflecting the postfault swing margin for the worst case generators in the system. The critical energy is dependent upon the location of faults and is usually estimated rather than precisely calculated. 4) Voltage sensitivity to reactive power (∂V /∂Q): Index for long-term small-disturbance voltage stability. This can be either a local (scalar) or a global (matrix) index. 5) Short-circuit ratio (SCR): Index for short-term voltage stability. This was originally defined as a local index but was extended to be a global index (generalized SCR) [30].
The creation of such indices is very helpful for system operators as they provide awareness of system conditions and guidelines to be used in planning and operations, so as to stay distant from high-risk operating conditions. Many of the indices also support sensitivity analysis to enable root-cause identification of particular nodes or apparatuses.

III. N E W C H A R A C T E R I S T I C S O F I B R s
The emergence of IBRs brings new and different characteristics to power system stability. That said, some principles in the classic framework are universal to ac transmissions, such as the central roles of voltage, angle, and frequency. It is also worth preserving the spirit of the categorization to decompose problems and the index formulation to summarize the complex dynamics and deliver informative guidelines to system operators. In this spirit, the classic framework is revisited in [24] with extensions suggested highlighting the "converter-driven stability" at high frequency due to the excitation of EM resonance by fast inverter control. Such an extension is sufficient to reflect the new characteristics at current IBR penetration. In this article, we aim at a further horizon into the future by envisaging a grid with nearly 100% IBR. The dynamics of power systems will be reshaped more radically in such future scenarios. As set out in the introduction, a GFM inverter has a similar operating principle to that of a synchronous generator: they both appear as a voltage source and synchronize to the grid via a power-angle swing. Therefore, a GFM inverter has a compatible analytical model with a synchronous generator and can be more readily accommodated in the classic framework. A GFL inverter, on the other hand, has very distinct behaviors that pose a challenge for analytical modeling, especially when GFL inverters become the majority in a grid. We propose three new characteristics, namely, generalized angle dynamics, indirect voltage control, and fast frequency dynamics as ways to accommodate both GFL and GFM inverters within the classic stability categories. The "converter-driven stability" defined [24] is carried forward into our framework but is renamed "EM stability" to distinguish it from other inverter-driven stability characteristics. An overview of the proposed new characteristics is provided in Fig. 2, and the details are presented in the following.

A. Electromagnetic Stability
As shown in Fig. 2, the first new characteristic introduced by IBRs is the EM stability related to EM states beyond the classic voltage-angle-frequency categorization. We define the EM states as the instantaneous voltages (or charges) on capacitances and currents (or fluxes) through inductances. Importantly, the voltage and current here are instantaneous values rather than root-mean-square (rms) values to be differentiated from the voltage in the classic framework. The fast and active control of inverters may introduce negative damping at supersynchronous and subsynchronous frequencies that overlap with the LC resonant frequencies and, thus, excite EM resonances. This negative damping may present in both GFL and GFM inverters, e.g., in the inner current loops of GFL/GFM inverters and the voltage loops of GFM inverters. There are two types of EM resonances: 2) series resonance related series compensators usually at a subsynchronous frequency.
The corresponding sources of inverter-induced negative damping are given as follows: 1) the switching-cycle delay of inverter modulation and control in both GFL and GFM inverters [31]; 2) the inner current loop of a GFL inverter interacting with the induction generation effect of a doubly fed induction generator (DFIG) [32].
Other inverter control loops, such as PLLs and dc-link voltage control loops, may also induce negative damping at a subsynchronous frequency [33], [34], but these effects are categorized as generalized angle dynamics and discussed in Section III-B. The EM stability at the supersynchronous frequency is often called harmonic stability in the literature [35]. The frequency scale of EM stability ranges from subsynchronous frequency (around half a fundamental frequency) to a couple of kilohertz. Switching frequencies have a significant impact on the frequency scale since the bandwidth of the inner control loop and the resonant frequency of output filters are all associated with switching frequencies.
EM stability is generally observed as a local issue, that is, only related to one or a cluster of inverters that are electrically or physically adjacent. The reasons for this localization are not adequately addressed in analytical terms in the literature but, intuitively, may be explained by the attenuation of high-frequency signals as they propagate through lines and cables of the network. EM stability is commonly considered a small-disturbance problem because there are no apparent differences in the oscillations subject to small and large disturbances, but there are reports of exceptions in that subsynchronous resonances may be sensitive to inverter saturation, and chaotic nonlinear behaviors may be induced under large disturbances [36].

B. Generalized Angle Dynamics
GFL inverters use PLLs to synchronize to a grid, which is a very different synchronization mechanism from that of synchronous generators. There is an ongoing debate over the role of PLLs in power system stability. It was recently proposed that a PLL has a dual relationship to the rotor of a generator. This duality theory leads to a generalized view of angle dynamics that can be accommodated in the classic framework. We present a brief introduction of the duality theory in the following, and details can be found in [21] and [22].
It is well-known that the internal (rotor) frequency ω of a synchronous generator (as well as a GFM inverter) is governed by Newtonian mechanics where P is the mechanical power from the prime mover, P is the electrical power into the grid, and J is the total moment of inertia of the rotor and the turbine. The damping power due to friction and amortisseur winding has been neglected in this equation. The Newtonian mechanics in (1) also apply to GFM inverters with the state of the low-pass filter (LPF) in the frequency droop controller equivalent to the rotor speed [20]. The internal frequency of a GFL inverter is governed by a PLL whose dynamics are very different to (1) in appearance but can be transformed (using the analysis in Appendix A) to a similar format where ω is the PLL internal frequency, Q = iqv d is called prime reactive power (the one set as reference), Q = iqv d − i d vq is the observed reactive power (the one injected into the grid), 2 K PLLI is the integral gain of the PLL, and we use v and i with subscripts d and q to denote voltage and current in the dq-axes. The proportional gain of the PLL is not included in (2) since this does not affect the essence of what we are illustrating.
In the light of (2), the motion of a PLL is also governed by Newtonian-like mechanics except that reactive power instead of active power determines the acceleration of ω. Therefore, a PLL can be represented as a reactive inertia, in contrast to the active inertia of a rotor. 3 As a result, the synchronization of inverters and generators (including GFM inverters) is unified in a single formulation, which is called the generalized angle dynamics. The synchronization among GFL inverters is governed by reactive-powercurrent-angle (Q-δi) swing (the supporting analysis is given in Appendix B), whereas the synchronization among synchronous generators and GFM inverters is governed by active-power-voltage-angle (P -δv) swing, as illustrated in Fig. 3. Both Q-δi and P -δv swings are governed by sinusoidal power-angle functions, noting that Q-δi swing depends on the load resistance R, whereas the P -δv swing depends on the line reactance X.
The generalized angle dynamics imply the feasibility of GFL inverters feeding an islanded grid with neither synchronous generators nor GFM inverters present. In such a case, the frequency of the grid is governed by reactive power balancing. For the common PLL with a proportional-integral (PI) controller (PI-PLL), an imbalance of reactive power results in sustained acceleration/deceleration of PLL frequency, as illustrated in 2 The reactive power defined in this article is opposite to the common definition, that is, Q −Im(ṼĨ * ), whereṼ andĨ are phasors of voltage and current. The Q is defined to align with the Newtonian-like formulation in (2). We use generator convention for all variables throughout this article. 3 It is worth noting that the reactive inertia of a PLL is scaled by i d , which is changing in real time, whereas the active inertia of a rotor is constant. Fig. 4(a). This is similar to the free-spinning of a synchronous generator if the prime mover generates constant mechanical power without a speed governor, as shown in Fig. 4(b). To create stable frequency for islanded GFL inverters, the PI controller in the PLL needs to be replaced with an LPF, creating LPF-PLL to remove the infinite gain at dc. The LPF-PLL will create a frequency-reactive-power droop (f -Q droop) characteristic in GFL inverters, which is similar to the frequency-active-power droop (f -P droop) characteristic of synchronous generators or GFM inverters, as illustrated in Fig. 4(c)-(d). For islanded GFL inverters with f -Q droop (LPF-PLL), the equilibrium frequency is determined by the net demand of reactive power Qnet, and the blue dot in Fig. 4(c) is an example.
Commonly, GFL inverters are connected to a grid that also includes GFM apparatuses (either generators or inverters), and in principle, both types of droop function act but are affected by the actions of the other apparatuses. The GFL inverters will contribute a certain P that offsets the f -P of the GFM, but the GFM apparatuses will provide stiff voltage sources with freely available Q that greatly reduces Q of the GFL, and in effect, the f -Q droop of the GFL inverters is overridden such that the frequency of the system is determined by the net demand of active power Pnet from the GFM, marked by the red dot in Fig. 4(d).
The block diagrams of PLLs in various forms are illustrated in Fig. 5 to facilitate the readers, and the detailed mathematical explanation is given in Appendix A for brevity. These diagrams indicate that Q of GFL inverters is controlled by the PLL. However, the principal mechanism  to control Q is to adjust the reference iq of the inner current loop. In the local PLL frame, i d + jiq = I∠φ, where φ = arctan(iq/i d ) is the power factor angle and It is clear that iq determines φ, which adds to the PLL angle θ to create the actual angle of the current phasor θi = θ + φ, as shown in Fig. 5(b). The power factor angle φ has a feedforward effect in Q control, which bypasses the PLL during transients for a fast Q response, and a PLL only regulates the steady-state error Q −Q = i d vq on top of the feedforward term. This angle feedforward mechanism does not exist in a synchronous generator since the voltage angle θv of a generator is pinned to its rotor angle θ, that is, θv = θ. However, it is possible to implement such feedforward in a GFM inverter, which provides extra flexibility in P control that has not yet been explored. This potential for angle feedforward is recognized in Fig. 6 with a zero-valued addition to θ.
The generalized angle dynamics of inverters and generators can be combined into a unified whole-system synchronization model, as illustrated in Fig. 6. A GFL inverter injects a current phasor I∠θi into the grid and receives a voltage phasor V ∠θ v from the grid. A GFM apparatus (inverter or generator) does the opposite: it applies a voltage phasor V ∠θv onto the grid and receives a current phasor I ∠θ i from the grid. The conjugate product between voltages and currents yields active and reactive powers that feedback into the active and reactive inertias (rotor and PLL) to close the loop of synchronization.
The relationships of the current and voltage phasors in a grid can be described by the nodal admittance matrix where we use the upright bold symbols to denote vectors and matrices representing multiple inverters and generators. We partition the nodal admittance matrix Y to distinguish the matrices related to GFM and GFL apparatuses, denoted by the subscripts M and L, respectively. In (3), V∠θv and I∠θi are known (specified by apparatuses), from which we solve the unknown I ∠θ is called the hybrid admittance/impedance matrix [22]. G describes the interaction of voltage (GFM) and current (GFL) phasors across the grid. It is important to distinguish between voltage angles and current angles in the generalized angle dynamics. The classic form of angle dynamics addresses only voltage angles, whereas GFL inverters with fast current control impose current angles in the first instance. The mapping from a current angle to a voltage angle is dependent upon the external impedance seen by the inverter, that is, the grid impedance, as illustrated in Fig. 7. If a GFL inverter is connected to a passive load, the voltage angle will be aligned with the current angle with an offset of the load power factor angle, that is, as illustrated in Fig. 7(a). On the other hand, if a GFL inverter is connected to an active network, such as voltage  source behind a reactance, then its terminal voltage angle is determined by the angle of the remote voltage plus a phase shift γ created by the voltage drop on the reactance, as illustrated in Fig. 7(b). In this case, which indicates that γ is determined by the active current, i d . For a basic GFL inverter, i d is set by dc-link voltage control (active power balancing), so dc-link control needs to be considered as a part of the generalized angle dynamics [37]. The PLL and dc-link control that participates in the generalized angle dynamics may have bandwidths of up to 20 Hz, and so the generalized angle dynamics may be expected to have a faster response than the classic angle dynamics, as was noted in the expansion of the timescale as marked by the dashed box in Fig. 2. Because of this expansion, there may be interaction with the EM dynamics. Such interactions may induce resonance-like oscillations at a subsynchronous frequency, especially in a series compensated grid or a weak grid, but these oscillations have a different mechanism to EM resonances [33], [34].
The generalized angle dynamics based on the GFL-GFM duality theory are important in many ways. First, it provides a unified analytical model to study GFL and GFM apparatuses together so that their roles in the whole-system stability can be evaluated to provide guidance on the placement of GFM apparatuses (either synchronous generators/condensers or GFM inverters). Second, GFM inverters may fall back to the GFL mode during faults if not designed with overcurrent capacity. Therefore, the duality theory may serve as the basis for transient stability analysis where such changes occur.

C. Indirect Voltage Control
Another significant difference between synchronous generators and GFL inverters is in voltage control. A GFL inverter can be set to control its terminal voltage amplitude V indirectly by injecting Q to the grid, as illustrated in Fig. 8(a). The effect of Q injection is dependent upon the external grid condition, including the equivalent Thevenin voltage and impedance of the grid seen by the inverter. The grid condition can be characterized by the V -Q curves displayed in Fig. 8(c). It is clear that the indirect control, represented by the horizontal line in blue (constant Q), has two intersection points with the V -Q curve, representing two possible equilibrium points. The two equilibrium points have opposite tendencies of stability due to opposite ∂V /∂Q sensitivities around them [38], which is called bifurcation in stability theories [39]. The point in the V -Q curve, where ∂V /∂Q changes' sign is called the knee point (or the nose point), which is the critical point for bifurcation. The existence of bifurcation indicates the possibilities of voltage instabilities for indirect voltage control, including: 1) voltage collapses where the voltage rapidly drops below the knee point and 2) voltage oscillations where the voltage oscillates near the knee point persistently or until diverges. From Fig. 8(c), we see that the voltage of the knee point increases with P , which  means that the risk of voltage instability increases with active power loading.

. Direct and indirect voltage control schemes. (a) Indirect voltage control of GFL inverters by injection of Q. (b) Direct voltage control of synchronous generators by setting V. (c) Curves (in black) of the V-Q characteristic of the grid at connection point (all variables are given in per unit) with: 1) a vertical line of imposed V from direct voltage control that has one intersection point with each V-Q curve and is always stable and 2) a horizontal line of injected Q from indirect voltage control that has two intersection points
In contrast, a synchronous generator or a GFM inverter controls its internal voltage (both amplitude and angle) directly, which is independent of the external grid. In the transient stability analysis, a synchronous generator is modeled as a voltage source V ∠θ behind a transient reactance X , where the amplitude V , proportional to the field flux φ f d , is governed by the exciter, and the angle θ is specified by the rotor [29]. A GFM inverter has a similar direct voltage control, as discussed in Section III-F. The direct voltage control is always stable, as is clear to see from Fig. 8(c), where the direct voltage control represented by the vertical line in red (constant V ) has only one intersection point with each of the V -Q curves.
In a conventional power system, indirect voltage control is only implemented on the load side via static Var compensators (SVCs) or static synchronous compensators (STAT-COMs), but direct voltage control is used on the source side (by generator exciters), and this is the backbone of voltage stability. With the increasing penetration of IBRs, indirect voltage control becomes ubiquitous on the sources side, which may narrow the stable range of voltage. This is especially of concern during a transient when: 1) even the indirect voltage control is (partially) disabled due to current saturation of inverters and 2) the transient angle perturbation may worsen the V -Q sensitivity due to angle-voltage coupling (as explained in the following). The V -Q control of IBRs has a much faster response than conventional voltage control mechanisms, 4 which gives rise to the expansion of the timescale of voltage stability, as marked by the dashed box in Fig. 2. This raises the possibility of subsecond dynamic voltage instability that is not observed in classic grids. Dynamic voltage instability is only anticipated from this analysis rather than observed in the real world. That said, early signs might have been shown in practice. For example, the Hornsea oscillation in the United Kingdom in 2019 may be related to dynamic voltage instability because: 1) there are clear oscillations in both V and Q observed; 2) at the time in question, the Hornsea system identified a weak grid; and 3) the oscillation frequency coincided with the time frame of V -Q control and caused tripping on overcurrent in just 0.25 s [12].
Indirect voltage control may also result in tight coupling between voltage and angle dynamics for the following reasons.

1) In practice, indirect voltage control is implemented
by setting the reactive current iq, which is expected to be proportional to Q, but this is only so when the PLL angle is properly aligned to the terminal voltage angle. A PLL may become misaligned during large transients (e.g., a phase jump or resynchronization after a fault), and therefore, voltage control is affected by angle dynamics. 2) A voltage control loop may overlap in timescale with a PLL and a dc-link control loop (dc-link control is a part of the generalized angle dynamics, as explained in Section III-B). An initial study of the voltage-angle coupling has been presented in [40] in which it is explained that a fast ac voltage control is helpful for PLL and dc-link stability. The voltage-angle coupling implies that joint voltage-angle stability studies may be necessary for inverter-dominated grids.
The distinction between direct and indirect voltage control also provides a perspective to differentiate GFL and GFM inverters. As discussed in Section III-B, GFL inverters synchronize to a grid via reactive inertia, which has a duality relationship to the active inertia of GFM inverters. GFL inverters are also capable of providing frequency response to a grid via either synthetic inertia or primary response if appropriate supplementary controls are provided. As a result, there seems to be no reason to disregard GFL inverters when seeking ways to address angle and frequency stabilities. However, GFL and GFM inverters have very distinct roles in voltage stability for the reasons explained above. The implication of this observation is that voltage stability may be the issue where GFL inverters make an inadequate contribution to grid stability and may be the principal reason why synchronous generators or GFM inverters are necessary to stabilize a grid.

D. Fast Frequency Dynamics
Another new characteristic of IBR-dominated grids is the increase in the speed of the frequency dynamics. This occurs for the following reasons: 1) the reduced inertia of grids and 2) the increased speed of response of active power of IBRs. Compared to conventional generation with a governor and a steam valve that are relatively slow, IBRs regulate active power fast (within seconds or subsecond). The fast active power response to frequency error serves as a compensation for the reduced inertia of IBRs and reduced response to RoCoF, meaning that future grids can be less reliant on inertia to maintain frequency stability. However, this also extends the timescale of frequency dynamics down from the order of tens of seconds to seconds, thereby creating overlap with the angle dynamics, as marked by the dashed box in Fig. 2. This may have significant impacts on both angle and frequency stabilities. First, the fast frequency control may help with the damping of angle swing and, thereby, expand the angle stability limit, provided that it is controlled properly (otherwise it may even destabilize the angle swing, e.g., due to excessive delay). Second, it may no longer be reasonable to assume a single frequency throughout the grid (represented by the motion of the CoI with a lumped inertia) in frequency stability studies. The angle swing around the CoI creates a spatial distribution of frequencies during transients [41], [42], which may result in very different RoCoFs and nadirs at different nodes. This may affect the configuration of frequency-based protection schemes. For example, it may be necessary to consider or filter out the angle swing in RoCoF and frequency measurement for lossof-main protection or underfrequency load shedding.

E. Asymmetric Inertia
Frequency stability is strongly related to the total inertia of a grid. IBRs have energy storage components on their dc-link (capacitors) and prime movers (e.g., wind turbines and batteries). In this section, we will discuss the roles of these energy storage components in system inertia.
The balance of active power within an IBR is reflected via the total energy E of its energy storage componentṡ where P and P are, respectively, the active powers delivered from the resource (e.g., wind) and exported into the grid (losses are neglected). E has different realizations for different IBRs. For example, E = (1/2)C dc v 2 dc for a solar inverter, where v dc and C dc are the dc-link voltage and the capacitance, respectively, and are typically rather small. For a wind turbine, there is additional (and much larger) energy storage in the inertia of the turbine E = (1/2)C dc v 2 dc + (1/2)Jtω 2 t , where ωt and Jt are the wind turbine speed and inertia, respectively.
Summing the active power balancing equations (1) and (9)  in which the subscript [n] is the index of an apparatus, M and N are the index sets of inverters and generators, respectively, P Σ is the total power generated by prime movers (or released from resources), and PΣ is the grid's total electrical power consumption that equals the net load if we neglect the losses. Equation (10) indicates that the net dynamic power mismatch between generation and load PΔ P Σ − PΣ is absorbed by all rotating inertias and other energy storage components in a grid, which raises the possibility that energy storage components of IBRs may have a similar role to rotating inertias of generators.
However, the sharing of PΔ among the rotating inertia and other energy storage components is uneven and depends upon the locations of the power perturbation and the control arrangements of the inverters. Taking first the case illustrated in Fig. 9(a), if PΔ originates with renewable generation through, for example, a surge of wind power, the wind turbine speed increases first, followed by the rise of dc-link capacitor voltage. This power perturbation is then propagated to the inertias of generators through the actions of dc-link control of the wind turbine inverter. This process is similar to the angle swing of a synchronous generator, and in such a case, the wind turbine and dc-link capacitor have an inertia-like effect in the system frequency response [37]. On the other hand, if PΔ originates with loads, as illustrated in Fig. 9(b), the wind turbine and dc-link capacitors do not see this perturbation because basic GFL inverters are configured to reject grid disturbances (using fast current regulation, feedforward control, and so on). As a result, the rotating inertias alone absorb this perturbation. It is clear that energy storage components of inverters have an asymmetric response to PΔ arising from the prime-mover side and from the grid side with the consequence that the effective inertia of the grid depends on where power perturbation originates. In this sense, the role of synthetic inertia is to link the energy storage components in IBRs (e.g., dc-link capacitors and wind turbine inertia) to the grid frequency so that they have even and symmetric participation in absorbing power perturbations, as illustrated in Fig. 9(c).

F. GFM Inverters
Up to this point, we have assumed that GFM inverters have similar behaviors to synchronous generators but have not provided a justification, which we will discuss here. A concept related to GFM is virtual synchronous generator (VSG), which aims at an exact replication of synchronous generators by inverters [43]. However, exact replication is often infeasible (due to the limitation of voltage, current, and switching frequency of an inverter) and unnecessary (since we only want to emulate the aspects of the behavior that is helpful for grid stability). GFM can be seen as a generalization of VSG, which encompasses the essential features of synchronous generators but also introduces new characteristics adapted to the constraints and flexibilities of inverters. As illustrated in Fig. 10, a GFM inverter contains an internal frequency ω governed by P -f droop control. Integration of frequency yields the internal angle θ that is combined with the amplitude V to create the voltage phasor V ∠θ. There are several ways to synthesize the voltage phasor V ∠θ with an inverter. The first and most common method is to use the double-loop control, as shown in Fig. 10(a) [18], [44]. The double-loop voltage control contains a voltage loop and a current loop. The current loop is similar to the one in a GFL inverter, and the voltage loop is cascaded upon the current loop. The current loop has two functions: 1) to damp the LCL resonance in the output filter and 2) to provide automatic mode switching to current control via the saturation of the voltage controller when excessive current would otherwise be drawn, such as during faults. The second way to synthesize the voltage phasor V ∠θ with an inverter is to use open-loop control in which the desired voltage is directly modulated on the pulsewidth-modulation (PWM) signal [45]. The open-loop control is more suitable for inverters with low switching frequencies, where the switching cycle delay in the current and voltage loops may excite LCL resonance. It is also possible to use single-loop control to synthesize the voltage phasor, but the tuning of control parameters is difficult and sensitive to resonant frequencies [46]. Both the double-and open-loop voltage controls of GFM inverters can be modeled as a controlled voltage source V ∠θ behind an equivalent reactance Xe, as illustrated in Fig. 10(c). For the double-loop control, Xe is shaped by the control parameters [47], and for open-loop control, Xe equals the total output reactance.
It is important to note that the P -f droop control in a GFM inverter usually includes an LPF within the gain, which gives dynamics that can be transformed to an equivalent rotor with feedback damping [20], as illustrated in Fig. 10. A first-order LPF with a transfer function of (1/Js+D) is a feedback connection of an integrator (1/Js) and a gain D in which the integrator is equivalent to a rotor with inertia J, and the gain is equivalent to damping coefficient. The equivalent inertia is synthesized via control, but the corresponding physical energy exchange is provided from the dc link or the prime mover, so, generally, a GFM inverter requires increased energy storage capacity from either dc-link capacitors, batteries, or exploitation of wind turbine inertia.
The equivalent model of a GFM inverter in Fig. 10(c) is very similar to the model of a synchronous generator, but it is important to recognize that the model holds only in normal operation. During faults, a GFM inverter must limit its current via either the automatic saturation of the voltage loop (in double-loop control) or the insertion of current control (in open-or single-loop control). Thus, a GFM inverter switches to a GFL-like current injection and is no longer a voltage source. Consequently, the synchronization mechanism should also be switched to a PLL temporarily to maintain angle stability [48] or should simply freeze temporarily at a constant frequency. As a result, a GFM inverter loses GFM function and falls back to a GFL inverter during faults. Considering that one of the most important motivations for adopting GFM inverters is to support transient stability of the grid, the transient fallback to GFL significantly compromises the benefits of GFM. There is ongoing research, for instance, [49] on the technologies to increase the overcurrent capacities of inverters to allow for a "true" GFM ability during transients, but this inevitably comes at the price of extra hardware investments.
Another important difference between GFM inverters and synchronous generators is that GFM inverters have much higher flexibility for fast and adaptive responses. For example, it is much easier to implement frequency damping on GFM inverters than on synchronous generators since the P -f control of inverters is much faster than the speed governing of generators. Both the equivalent inertia and damping can be rescheduled online in line with the operating condition, which can be helpful for first-swing stability [50].

G. Summary
In summary, the new characteristics of IBRs (both GFL and GFM) can potentially be accommodated in the classic power system stability framework with certain revisions and extensions. The extension includes introducing an additional category of EM stability near and above synchronous frequency, introducing new complexities in angle stability (generalized angle dynamics) and voltage stability (indirect voltage control), and expansion of the timescale of angle, voltage, and frequency dynamics to recognize the faster dynamics of IBR. The extended timescales cause previously distinct and easily partitioned features to possibly

Fig. 11. Overview of the state-of-the-art techniques used in the power system stability analysis.
overlap and result in coupling across different stability categories in ways that pose challenges for stability analysis. The energy storage components in IBRs are not directly coupled to grid frequency and may have an asymmetric contribution to system inertia, depending on how the inverters are controlled. GFM inverters have similar but more flexible behaviors compared to synchronous generators, but the transient fallback to the GFL mode compromises the benefits.

IV. S T A B I L I T Y A N A L Y S I S T E C H N I Q U E S
Alongside the classic stability framework is a series of stability analysis techniques tailored to address each category, and they have served the industry well over many decades. This section provides an overview, summarized in Fig. 11, of the state-of-the-art stability analysis techniques with a discussion of whether they remain sufficient to address the new characteristics of IBRs.

A. Models
Stability analysis rests on appropriate models of the system. The full physical model of a power system is a switching and nonlinear system, which is intractable for analytical methods. The switching feature of power systems is commonly addressed by averaging over an appropriate period to form a continuous model. The averaging period can be a PWM cycle (micro to milliseconds) for self-commutated devices, a synchronous cycle (1/50 or 1/60 s) for line-commutated thyristors, or many cycles for switchgear controlled devices. The resulting continuous model is a smooth and possibly nonlinear ordinary differential equation (ODE), which can be linearized around the operating equilibrium to obtain a linear model. A linear model can be described explicitly as state-space equations or implicitly as transfer functions. If the operating equilibrium is nonsinusoidal or unbalanced in three phases, the linearized model may be time-periodic; in such cases, the state-space matrices may be time-periodic, and the transfer functions need to be extended to harmonic transfer functions. The progression of the models from physical, through continuous, to linear models brings a widening of the set of analytical techniques available but at the price of reducing fidelity.
Power system models are often very high order in their full form, and model order reduction is commonly applied to ease analysis. There are two major techniques for order reduction based on insights into the physical system. The first is the model aggregation of similar and adjacent apparatuses, e.g., aggregating a wind farm into a single equivalent wind turbine [51] or aggregating a cluster of generators in an area into a single equivalent generator [52]. The second technique is model truncation according to the separation of timescale [47], [53]. For example, EM states are truncated in the angle stability analysis. There are other model truncation methods available in the systems theory, but they are not widely used in power engineering [54], [55].

B. Nonlinear Analysis
Nonlinearity is a significant feature of power system dynamics and is the essential difference between largeand small-disturbance stabilities. There are two main causes of nonlinearity in power systems, namely, the presence of trigonometric functions in power flows and saturation (e.g., the ramp-rate saturation of power plants, the magnetic saturation of iron cores, and the current-limit and modulation-index saturation of inverters), but many other nonlinear features also exist. is the most straightforward method for nonlinear analysis [56]. It can be applied to both the full physical model (including switching) and continuous models (via averaging). Both forms are computationally intensive but switching models particularly so. Time-domain simulation can provide high-fidelity validation of stability but requires exhaustive scanning across many operating conditions, with consequent large computational effort, for full assurance, but the results are often hard to interpret and generalize.
Analytical methods for nonlinear stability analysis have been developed to overcome the shortcomings of timedomain simulation. The Lyapunov method (also called the direct method) is a widely used analytical method. The Lyapunov method provides an interpretable evaluation of stability margin (e.g., critical energy) in nonlinear dynamics and has superior computation efficiency over timedomain simulation. However, the selection of the Lyapunov function is often challenging and reliant on insights from the structure of models. As a result, the application of the Lyapunov method in power systems is largely confined to angle stability where the analogy of swing equations to Newtonian mechanics yields an energy-like Lyapunov function. An alternative to the Lyapunov method is the geometric method [57], which uses the manifold theory to determine or estimate the stability boundary of a power system. The geometric method is less reliant on the structure of models and has more general applicability. It is possible to combine the geometric method with the Lyapunov method, e.g., in the controlling unstable equilibrium point (CUEP) method for estimating the critical energy [58].
Analytical methods and time-domain simulation are often used in complement. For example, the estimation of critical energy in the Lyapunov method applied to transient stability requires the calculation of the trajectory of the system, while a fault is present, and this is often obtained by time-domain simulation. Furthermore, analytical methods help to interpret the time-domain simulation results (e.g., by indicating the mode of interaction and oscillation boundary from the simulated trajectories) and reduce the extent of exhaustive scanning needed.
Analytical nonlinear analysis has long been applied to the large-disturbance angle stability of classic power systems. There have been efforts recently to extend the analytical method to inverter-based grids. The PLL angle stability of a GFL inverter is investigated in [14] using the Lyapunov method, which yields a stability criterion that is similar to the equal-area criterion (EAC) of a generator, as illustrated in Fig. 12. In the light of the duality of inverters and generators, such an extension is natural. However, it is difficult to apply the Lyapunov method to multi-inverter (GFL) and inverter-generator hybrid grids since the equivalent potential energy of a multi-inverter or hybrid grid depends significantly on the path of integration and, therefore, is unknown without an explicit solution of the phase trajectory [59]. 5 It is even more difficult to apply the Lyapunov method to GFL inverters when other control loops (e.g., dc-link control and ac voltage control) are considered since these loops are software-defined and lack a natural energy function.
Besides determining the stability and assessing stability margin, system operators are also interested in the mode of interaction during the transients. It has been observed and partially proved in classic grids that the mode of interaction is almost always a generator or cluster of generators swinging against the remaining part of the grid. This leads to the extended EAC (EEAC) [60] in which a multigenerator swing is reduced to a two-area swing. This two-area mode of interaction also indicates the boundary (critical cut-set) of the transient swing, which is very useful for system operation and planning. It is not clear whether such a mode of interaction is still dominant in an inverter-based grid or whether other significant modes of interaction exist.

C. Linear Analysis
The limited number of analytical tools available for nonlinear analysis and the constraint on the structure of models (for a feasible energy function) leads to the use of linear analysis based on linear approximate models of the dynamics near the system equilibrium. The linear analysis can produce useful insights not available through nonlinear analysis but at the price of losing part of the entire picture of the off-equilibrium dynamics. Therefore, linear analysis is often used as a complement rather than a replacement for nonlinear analysis. For example, the control parameters of inverters and generators are often tuned by linear analysis to ensure a sufficient damping ratio of the associated eigenvalues (poles), and then, this is verified against large disturbances via nonlinear timedomain simulation.
There are two types of mathematical descriptions for a linear model, namely, state-space models and transfer function models. A state-space model explicitly represents all internal states of each component in a power system and is, therefore, a white-box model. State-space models offer in-depth analysis based on linear algebra, including eigenvalues, eigenvectors (mode shapes), and participation factors (eigenvalue sensitivities), which provide insights into stability margin, mode of interaction, and origins of instability [29]. These methods have been widely used in power system stability analysis and have been built into commercial software, e.g., DIgSI-LENT/Powerfactory [61]. Transfer function models are represented either as polynomials in the Laplace variable s or frequency spectra. A transfer function describes system dynamics via the input-output responses at ports of a system without all internal details explicitly represented, and it is, therefore, called a black-box model. Transfer function models generally provide less information than state-space models, but they are used for the following reasons.

1) A transfer function model does not involve the disclo-
sure of the internal design and the control of apparatuses, which helps preserve commercially confidential detail. This feature is especially relevant for inverters since the model of an inverter is not as standardized as the model of a synchronous generator. 2) A time delay can be conveniently represented in a transfer function as e −T d s (T d is the time delay), whereas it is relatively difficult to represent a time delay in a time-continuous state-space model [31]. 3) Transfer function models provide a straightforward way to study frequency coupling in unbalanced three-phase systems and nonsinusoidal systems via harmonic transfer functions [62], [63]. 4) Transfer functions can be extended to describing functions to address the saturation nonlinearity with the quasi-linear analysis [64].
Transfer function models are used for both generators and inverters, but the choice of ports at which the input-output relation is defined can be very different, depending on the observability of the modes over the ports. For a generator, where rotor swing is mainly of concern, it is more convenient to select the port on the mechanical side with a torque input and a speed output, and the corresponding transfer function can be considered mechanical admittances [29]. For an inverter, where electrical interaction is mainly of concern, the port is usually selected on the electrical side with voltage inputs and current outputs, and the corresponding transfer functions are electrical admittances [65]. In some studies, admittance is used interchangeably with impedance (the inverse of an admittance) [66]. It is also possible to open other ports for transfer function formation, e.g., the dc-link port, the PLL port, and the exciter port. All these ports give identical results in terms of stability but may provide different interpretability and observability for various oscillation modes [67]. The transfer functions for these generalized ports can be described by multiport networks that are interconnected over the grid through a component connection method [68].
Electrical admittances (and impedances) are specific to the coordinate frame in which the voltages and currents are observed. The relationship between one frame and another (e.g., a frame aligned to a particular PLL or rotor) is subject to dynamics, which must be accounted for in transforming admittances between frames. Frame transformation can be achieved by embedding the frame dynamics [69]. It is noted in [69] that a global steady frame (steady means rotating at a constant speed) should be used to define the admittances of multiple apparatuses in an interconnected grid. It is possible to define admittances in a stationary and complex frame, and obtain the same conclusion on stability [66], but a rotating frame may yield a simplified representation because of the diagonalization effect on harmonic transfer function matrices [70].
The admittance-based analysis is effective in addressing the EM resonances in inverters [35]. It is understood that EM resonances (both series and parallel) can be rendered unstable by a negative phase margin (or equivalently, negative conductance or resistance) at the resonant frequency arising from the switching-cycle time delay of the inner control loops. EM resonances can be mitigated by reducing the control gains, adding passive resistances, or by lead-lag compensation in the inner control loops. A lead compensation can be realized by a high-pass filter (including derivative control) or predictive control [71], and a lag compensation can be realized by a low-pass filter [46]. It is also possible to combine lead and lag compensation via a bandpass or band-stop filter [72]. Lead-lag compensation can be implemented on either real signals (e.g., v d and vq) or complex signals (e.g., v d +jvq). Complex-signal lead-lag compensation provides asymmetric phase compensation for the positive and negative sequence components, is, therefore, more flexible, and achieves better performances in some applications [32], [73].
The admittance-based analysis is especially effective in addressing coupling between frequency components that result from nonsinusoidal distortion or from unbalancing between phases. To do this, admittances are extended to harmonic transfer functions where a single-frequency input creates a multifrequency output. Frequency coupling is significant in the following two cases: 1) single-phase inverters where the double synchronous frequency (i.e., 100 or 120 Hz) ripple in the dc-link voltage induces frequency coupling [74] (it is worth noting that a modular multilevel converter is, in effect, three single-phase converters since it has separate dc energy storage for each arm [75]); 2) inverter-induced EM resonances close to or above the switching frequency, where the sideband coupling effect of PWM is no longer negligible (this is especially true in high-power converters with low switching frequencies combined with high-frequency parasitic resonances, e.g., in cables) [76]- [78].
There have also been efforts to investigate generalized angle stability via admittance-based analysis [65], [79], [80]. It has been proven that generalized angle dynamics are observable from electrical admittances, but the interpretation is not straightforward since the generalized angle dynamics of, for instance, a PLL or a dc-link, are reflected in the ac electrical port only indirectly. It is argued in [67] that it is advantageous to investigate generalized angle stability via other ports, e.g., a PLL port or a dc-link port, where the generalized angle dynamics are directly exposed.
The admittance-based method in its original form needs to partition the system into two parts, each part represented by an equivalent Norton admittance (or a Thevenin impedance), whose interaction determines system stability via Nyquist criteria or variants of it. The point of bipartition should be selected at the boundary of the mode of interaction under investigation to make the admittance-based analysis informative and interpretable. However, the boundary of interaction is generally not known beforehand, so the system bipartition itself is a difficult task. For this reason, the admittance-based method was limited to a specific local problem, that is, a single apparatus (or a cluster of adjacent apparatuses acting in concert) interacting with the rest of the grid, where the point of bipartition is naturally selected as the bus to which the apparatuses under investigation are connected, which is usually called the point of common coupling (PCC) [81].
To extend the admittance-based method beyond local analysis, the whole-system admittance model has been developed [69], [82], [83]. In the whole-system admittance model, the system under investigation is partitioned along with all buses, yielding a 2N × 2N admittance matrix (N is the number of buses), as illustrated in Fig. 13. The whole-system admittance matrix can be organized blockwise as N × N blocks, among which each block Y kl (k, l ∈ 1, 2, . . . , N) is a 2 × 2 matrix in the dq frame. The whole-system admittance model enables root-cause tracing for oscillations in power systems without full disclosure of the state-space (white-box) model, which is called the gray-box approach. It is proven in [84] that the sensitivity of the eigenvalue λ (representing an oscillation mode) against the local impedance Z k of an apparatus (by local we mean the impedance of the apparatus itself) is determined by the residue of the whole-system admittancê Y kk observed at the corresponding bus. Thus, the residues are defined as impedance participation factors that enable eigenvalue sensitivity analysis in admittance-impedance (black-box) models. The sensitivity analysis can further trace the internal parameters and the internal states of an apparatus via the chain rule, yielding the gray-box model. The gray-box approach only requires disclosure of the partial derivatives of local impedances with respect to internal parameters but can achieve the capability of root-cause tracing similar to explicit state-space models. The gray-box approach can be extended to other ports (e.g., PLL, dc-link, exciter, and shaft), as recognized by the dashed arrow in Fig. 11.

D. Decentralized Stability
The stability analysis methodologies described, thus far, all take a centralized point of view, that is, the whole system is modeled and analyzed as one. Such a centralized methodology is vulnerable to the changes in system models as apparatuses join and leave the grid, and operating points alter, which is increasingly inevitable due to the uncertainty of IBRs. Decentralized stability is proposed as a potential methodology to address this problem. The key idea of this methodology is to imbue a Fig. 13. Whole-system admittance and gray-box approach for   root-cause tracing without explicit state-space models.   (a) Conceptual illustration. (b) Mathematical formulation. ρ denotes   an internal parameter, and a denotes a diagonal entry of the   state-space matrix A associated with a particular state. certain dynamic property in all apparatuses in the grid so that whole-system stability is guaranteed regardless of the changes in individual items.
The first decentralized stability criterion is the passivity criterion, which guarantees the small-signal stability of a grid as long as the admittances or impedances of all grid-connected apparatuses are positive-real [85], [86]. Such passivization is only feasible at high frequency (above the synchronous frequency) and is proven impossible at low frequency where the outer control loops (e.g., PLL, dc-link control, and droop control) dominate. Therefore, the original passivity criterion is limited to EM stability where the frequency of interest is high enough. The passivity criterion can be relaxed to widen its applicability. For example, Pates and Mallada [87] developed a passivity-like criterion with multipliers for the transfer functions at mechanical ports that can ensure the angle stability of a multigenerator grid.
There are also nonlinear decentralized stability criteria that guarantee the large-signal stability of a grid. For example, Colombino et al. [88], Sinha et al. [89], and Johnson et al. [90] developed a novel synchronization scheme called virtual oscillator control that has an almost global stability region, meaning that a grid can converge to the equilibrium from any initial condition.
One of the common problems of decentralized stability is that it is often conservative and, therefore, results in increased costs or degraded performance. There has been no systematic research to quantify the conservativeness of decentralized stability approaches compared to centralized approaches. That said, the decentralized stability methodology provides elegant theories and certified stability criteria, which are independent of changes in individual items, and is, therefore, an important direction of research, especially for IBR-dominated grids.

V. U N S O L V E D P R O B L E M S A N D F U T U R E O P P O R T U N I T I E S
From the discussion in Sections III and IV, it is clear that knowledge and experience of power system stability with a high penetration of IBRs continue to grow, but there are also significant unsolved problems that need attention.
The first unsolved problem is that the mode of interactions among IBRs is still not sufficiently clear. Most state-of-the-art research focuses on local interactions, e.g., a single or a cluster of IBRs interacting with a weak grid represented by a high grid reactance. The wider interactions of multiple IBRs across a grid are almost unknown. The challenges to tackle this issue include: 1) the mode of interaction of IBRs is highly dependent upon control parameters and operating points [69] and 2) the cross-coupling between different stability categories (e.g., voltage and angle) further complicates the mode of interaction [40]. Clarifying the mode of interaction is important for ensuring system operability, and it will shed light on the dynamic structure of a power system, such as the origin and the boundary of oscillations, so that system operators understand the weak links in the chain of system stability and organize defense lines accordingly.
The second unsolved problem is how to deal with the tight constraints of inverters and the resulting hard model switching in transients. This is especially true for a GFM inverter where current saturation is often activated during faults, changing the nature of the model from a voltage source to a current source. The model switching results in a switched system, which is generally hard to analyze compared with nonswitched systems. Moreover, the model switching might be cascaded across the grid, thereby further complicating the mode of interactions. For example, a fault-induced current saturation may first occur in the inverters electrically close to the fault, which may then reduce grid voltages in neighboring areas and, thereby, trigger the current saturation of inverters further away from the fault. Such cascading phenomenon is not yet observed in practice, mainly because grids are still stiff enough with synchronous generators supporting the voltages but may be of concern in the future when the penetration of IBRs approaches 100%.
The third unsolved problem is how to take advantage of the flexibility of inverters. Inverters have very high control flexibility, and the controls currently applied exploit only a very small part of the flexibility. It is not clear if the GFM-GFL dualism is sufficient to categorize the behavior of inverters or if new behaviors should yet be defined. There is an apparent middle ground between GFM and GFL inverters, sometimes called grid supporting inverters, but the roles of these inverters in grid stability are yet to be investigated. Another potential approach to exploit the control flexibility of inverters is to reconfigure the inverter control algorithms or parameters in real time, i.e., adaptive control. It is reported in [91] that there is not a single PLL parameter that can address all different grid conditions (such as phase jump, fault ride-through, and weak grid), which implies that an adaptive PLL might be beneficial for GFL inverters. It is also demonstrated that adaptive inertia may be beneficial for the transient stability of GFM inverters [50]. However, the current research on adaptive control is limited to a simple system, i.e., a single inverter connected to an infinite bus. The challenge of deploying adaptive control at scale is that the adaptive controllers may interact with each other across the grid, which further increases the complexity of stability analysis.
The fourth unsolved problem is how IBRs interact with other apparatuses and controls in the grid. Synchronous generators will still be present in a 100% renewable grid due to the existence of biomass and hydro plants, and perhaps nuclear. The way a small minority of synchronous generators interact with the majority of inverters may be very different from the current situation where synchronous generators are still a significant fraction. The protection system will also have a strong impact on power system stability. The increasing penetration of IBRs will significantly change the fault behavior of grids (e.g., reduced fault currents and increased RoCoF), and a resolution is needed as to whether inverters should be designed to fit into the state-of-the-art protection systems or protection systems should be changed to fit into the new fault behavior of IBRs or where between these extremes the best answer lies.
The behavior of loads is another factor that makes great differences in the stability of an inverter-based grid, which has not had the attention it deserves. GFM inverters only serve as stiff voltage sources in normal operation and may fail to do if the GFM inverters are not designed for transient overcurrent. On the other hand, loads may draw higher currents during and after faults due to stall of induction motors, constant power loads, and remagnetization of transformers. If the fault-induced overcurrent of loads is not met by IBRs, transient voltage collapse may appear. The Western Electricity Coordinating Council (WECC) of the United States has published the composite load model [92] to investigate the fault-induced voltage collapse issues in distribution networks. This composite load model can be incorporated into IBR stability analysis and should be kept updated to reflect the changes of load composition (i.e., different proportions of motor, resistive, and electronic loads).
The unsolved problems listed above pose a systemic challenge and need to be addressed by new theories, new concepts, new tools, and new practices. In the following, we present our visions on the potential opportunities and directions for future research in control and stability analysis.

A. New Analytical Methods
The large-disturbance stability analysis has long been focused on the Lyapunov method. As explained in Section IV-B, the Lyapunov method might not be easily applied to multi-inverter (GFL) systems. There are two potential ways to resolve this issue. The first is to use energy-based control to synthesize inverter controllers to emulate an energy function (GFM control is one example of this), but this comes at the expense of setting aside part of the control flexibility of inverters and may result in increased cost. The second way is to explore nonlinear theories beyond the Lyapunov method, and the geometric method [57] is a promising candidate. The geometric method uses the joints of the stable manifolds of the unstable equilibrium points to determine the stability region around a stable equilibrium point. It is generally applicable to a wide class of nonlinear systems and can be automated via computer programs. That said, there might be difficulties in calculating and visualizing the stability region in high-dimensional spaces.
The decentralized stability is also an important direction of theoretic exploration. Decentralized stability offers stability guarantees under changing grid conditions and may fundamentally reshape how power systems are designed and operated, if successful. State-of-the-art decentralized stability theories are still immature, and efforts are needed to price the costs and benefits of different decentralized stability approaches to establish practical feasibility.

B. Separation of Timescales
As explained in Section II, the power system stability analysis is greatly simplified by the separation of timescales of different dynamics. In the classic grid, the separation of timescales is rather a "gift of nature," that is, the thermodynamics dynamics of prime movers, the mechanical dynamics of synchronous generators, and the EM dynamics of transmission lines are naturally decoupled in timescale. In an inverter-based grid, the timescales of different dynamics are governed by control laws rather than natural laws. Therefore, a methodology is needed to assign the timescale (bandwidth) of different control loops. It is a common practice to have fast inner loops (current loops for GFL inverters and current-voltage loops for GFM inverters) and slow outer loops. However, it is not clear how different outer loops should be organized on a timescale. For example, should the dc-link control be faster than the PLL or not, and is there a common rule of thumb that can be applied to all grid conditions? Moreover, the timescale may change when multiple apparatuses interact in a grid. For example, inverters connected in parallel result in reduced equivalent bandwidths of current loops [93]. These set out potential topics for future studies.

C. Stability Indices and System Strengths
As explained in Section II, the establishment of stability indices is a key objective of a stability framework, and this is what needs to be carried over to inverter-based grids. Some of the classic indices, e.g., V -Q sensitivity, can be inherited directly, and some others can be extended to accommodate the new characteristics of IBRs. For example, the synchronization and damping power coefficient can be extended according to the generalized angle dynamics [22].
Nonetheless, it is not as straightforward to migrate other classic indices to an inverter-based grid. In the classic framework, the grid strength is solely defined by the SCR. Inverters (especially GFM inverters) may have decidedly different behaviors under small and large disturbances, and the equivalent strength presented to the grid should also be evaluated separately. A possible solution is to use incremental impedance to represent small-disturbance grid strength and fault current to represent large-disturbance grid strength. It is worth noting that GFL inverters also have a partial capacity to strengthen the grid via indirect voltage control, and this may result in an asymmetric (salient) incremental impedance that can be counted into grid strength.
IBRs have the potential for flexible approaches to frequency regulation, including synthetic inertia, fast primary response, and GFM control. These frequency regulation approaches may couple with angle dynamics yielding a spatial distribution of frequency dynamics. As a result, a single inertia index is no longer sufficient. A possible way to address this issue is to define temporal-spatial inertia as follows: . . . where ΔP k is the perturbation of active power at bus k and Δω k is the corresponding perturbation of the grid angular frequency. H kl is the inertia transfer function whose frequency spectrum represents the temporal distribution of inertia. For example, for a physical inertia, H kl (s) = (1/ω 0 Js), and for a fast frequency response, H kl (s) = (K/1 + τ s), where J is the physical inertia, ω 0 is the nominal synchronous frequency, K is the frequency-power droop coefficient, and τ is the time-constant of the fast frequency response. H kl (s) is organized into a matrix that characterizes the spatial distribution of inertia. It is worth noting that system strengths (grid strength and frequency strength) are useful because of their physical interpretation but are not directly an indication of stability margin. More sophisticated indices are needed to quantify the cross-coupled voltage-angle-frequency dynamics of an inverter-based grid.

D. Stability as Services
On operational timescales, stability is a responsibility of a system operator, and the operator needs to ensure that sufficient provision of services that aid stability is in place either by mandating their provision or procuring the services in a market. Procurement of response and reserve services for frequency stability is well established, but these services are being redefined toward faster products in recognition of reducing inertia but also exploiting the new capabilities of IBR, such as batteries [94]. There is a need to consider how all other aspects of stability, such as angle stability, voltage stability, and subsynchronous and supersynchronous oscillations in their new forms (see Fig. 2), can be managed by services. For example, the Global Power Systems Transformation (G-PST) Consortium has published a position paper with system needs and services definitions in technology-neutral terms applicable to inverters and generators [95]. System operators are considering the possibilities of enforcing GFM and short-circuit capacities for IBRs in future grid codes, and the Electricity System Operator (ESO) in the United Kingdom has proposed a minimum specification for GFM capability [96], [97], which includes a discussion of inertia, damping, and reactive power support of voltage. These features would contribute to the frequency and angle stability and improvement of voltage profile and power transfer limits. The stability services can be provided from various IBRs, including renewable generations, HVdc systems [98], battery power conditioning systems, STATCOMs, and so on.
The definition of the stability services is closely related to stability indices and system strengths discussed above, that is, providing a particular service results in the enhancement of the corresponding index. The temporal-spatial distribution of stability indices may result in a more complicated service market.

E. Data-Driven Methods
Physics-based models are increasingly difficult because of the confidentiality of models, uncertainty on operating points, and possible incorrect configurations. On the other hand, the availability of data is rapidly growing due to the wide deployment of sensors and communication infrastructures for power systems. The recent progress in machine learning also provides techniques to process and understand the data. The rise of data may possibly compensate for the deficiency of physical models, making data-driven methods an important research focus.
There are two ways of exploiting data in power system stability analysis. The first is to use data for the acquisition of models. This is similar to system identification in control theory. The model acquisition can be based on passive observation [99] or active probing [100], [101]. Active probing can ensure sufficient excitation of all internal states but at the expense of extra probing hardware. It is beneficial to make use of domain knowledge to assist in data-driven model acquisition. For example, domain knowledge may provide the rough structure of models, and the data can be used to polish the model [101]. The gray-box approach described in Section IV-C may also facilitate the fusion of domain knowledge with data: data are used to identify the oscillation mode in a system and locate the apparatuses participating in this mode, and domain knowledge of manufacturers then helps to identify the solution to mitigate the oscillation without the need to disclose their knowledge to system operators.
The second way to use data is for the exploration of models. As explained in Section IV, there are challenges for analytical studies of the nonlinearity of inverter-based power systems. Data-driven methods are a potential pathway to circumvent this challenge. Extensive simulations can be conducted to generate a vast amount of data, and machine learning techniques can be used to encode features from the simulated data. The features may provide an approximation of the stability structure needed by operators, including the origin of instability, mode of interaction, and stability margin, and, thus, may to some extent substitute analytical studies.

VI. C O N C L U S I O N
The classic framework for analyzing power system stability retains a lot of value but needs to be adapted to accommodate the new behaviors of IBRs. The first is a need to recognize the wider and faster time frames over which angle, frequency, and voltage dynamics occur for inverters compared to generators such that these previously distinct categories start to couple to each other and couple to new fast EM modes. This may need an extension of analysis tools to insight into system behavior and new guidelines on inverter control to mitigate the coupling. The second is the recognition that GFL inverters have different behaviors to generators and GFM inverters, and have duality relationships. We show that generalized angle dynamics can accommodate behaviors of GFL and GFM, and the duality features show that GFL inverters alone are able to sustain a synchronous region, a behavior that may not be practical in normal operation but may be essential during faults when GFM current-limit and behave like GFL. GFL and GFM inverters make different contributions to voltage control and voltage strength. The traditional V -Q characteristic remains useful in analyzing voltage stability and reveals superior capabilities of GFM and a strong reason to provide GFM in a system. However, GFM can only maintain a strong voltage up to the point of the current limit, so their large disturbance stability is different and requires different analysis and index definitions for large and small disturbance stabilities.
Frequency dynamics are becoming faster as IBRs dominate grids and inertia reduces. Virtual and synthetic inertias from GFM and GFL inverters do lessen the effect, but inverters can also provide a fast frequency response that diminishes the need for inertia. However, the fast frequency dynamics will couple with angle dynamics leading to a need to routinely consider the spatial distribution of frequency dynamics and a need for more sophisticated measures of frequency strength than system-wide inertia.
With control design central to inverter performance, open models are rarely available because of confidentiality, and this makes it very difficult to apply some analytical methods of stability assessment and root-cause analysis. Data-driven methods may provide the way forward for both model acquisition and model exploration to partially substitute physical models and analytical studies. The data-driven methods should be combined with domain knowledge for interpretability. The gray-box approach may provide a pathway for data-knowledge fusion and conversation across different stakeholders without breaking the commercial confidentiality of knowledge.

A P P E N D I X A R E F O R M U L A T I O N O F P L L O F G F L I N V E R T E R S
This appendix explains how the dynamics of PLL can be reformulated to a Newtonian-like law. As discussed in Section III-B, there are different variants of PLL, and we start from the simple case with a single integration controller, whose control law iṡ where ω is the PLL internal frequency, K PLLI is the integral gain, and vq is the q-axis voltage measured in the local dq frame aligned to the PLL internal angle θ = Ê ωdt. Equation (12) can be rewritten for an inverter exporting real power via i d as which is exactly (2). If proportional gain K PLLP is included in the PLL, the control law becomeṡ Again, we multiply i d K −1 PLLI on the both sides of the equation Assuming that i d is constant, that is,i d = 0, we have As will be illustrated in Appendix B, Q is a function of the angle difference δi between the PLL angle and the grid angle soQ = ∂Q ∂δiδ i = ∂Q ∂δi ωs (17) where ωs is the slip frequency between the inverter internal current and the equivalent Norton current of the external grid. Combining (15)-(17) yields where τ PI = K PLLP K −1 PLLI . It is clear that the role of the proportional gain is to introduce damping power with respective to the slip frequency, which is similar to the effect of an amortisseur in a synchronous generator.
We have introduced the LFP-PLL in Section III-B, and we now present the control law for the LFP-PLL in the Newtonian-like formation as well. The LFP-PLL can be represented as the transfer function in the following: where K is the Q-f droop gain and τ is the time constant of the LPF. This transfer function can be transformed to the differential equation below (we make use of the fact that ω is constant soω = 0) τω + (ω − ω ) = K(Q − Q) (20) or equivalently It is clear that dynamic behaviors of LPF-PLL and PI-PLL are very similar except that: 1) the damping of LPF-PLL acts upon the deviation frequency, whereas the damping of PI-PLL acts upon the slip frequency and 2) the damping of LPF-PLL is explicit and constant, whereas the damping of PI-PLL is implicit and dependent upon the operating point. This implies that LPF-PLL may be a better candidate from the stability point of view. LPF-PLL may induce steady-state error of reactive power, but this can be eliminated via secondary control. Defining D = K −1 and J = τ K −1 , we rewrite (22) as which matches exactly the form of inertia and damping.
Here, we call J and D reactive inertia and reactive damping, respectively, for GFL, as counterparts to the active inertia and active damping for GFM.

A P P E N D I X B Qδ i C H A R A C T E R I S T I C O F G F L I N V E R T E R S
Using the annotations in Fig. 3(a), the reactive power generated by the first current source is where δi θ i1 − θ i2 . This is the Q-δi characteristic illustrated in Fig. 3(a). The Q-δi characteristic is dependent upon the power factor of the load. If R is replaced by a complex impedance Z∠ϕ, (23) becomes where the first part −I 2 1 Z sin ϕ is independent of δi and represents the reactive power supplied to the load, and the second part I 1 I 2 Z sin(δi − ϕ) is the revised Q-δi characteristic considering the load power factor. This is similar to the dependency of P -δv characteristic on the transmission line impedance angle in synchronous generators.