Stability Analysis of Grid-Forming MMC-HVDC Transmission Connected to Legacy Power Systems

The power system is going through a change in its very foundations. More and more power converters are being integrated into the electric grid to interface renewable energy resources and in high-voltage direct-current (HVDC) transmission systems. This article presents a discussion on the stability of power systems when HVDC transmission systems based on modular multilevel converters (MMC) are connected in grid-forming (GFM) mode to the legacy power system using concepts of energy functions and Lyapunov stability theory and considering aspects of the interoperability between GFM converter technologies. As a base for the stability analysis, we review the main GFM converter technologies (droop and virtual synchronous machine), highlighting their differences. Then, we present a model using the center-of-inertia formulation for a multi-machine/multi-GFM converter power system representing a close future scenario of power systems where GFM converters might adopt different technologies. To illustrate the theoretical Lyapunov-based stability analysis, simulations performed in Matlab/Simulink showed the behavior of a 12-bus test system during a frequency disturbance that originated from the sudden connection of a load. To reflect the interoperability of different GFM technologies and the power system, scenarios with one single GFM technology and a scenario with mixed technologies were investigated. For the test system considered, the frequency response with fewer oscillations and a higher frequency nadir was obtained when all GFM converters were operated as VSMs that have a higher inertial response contribution.


Introduction
To face the pressing climate crisis and to counteract its effects, all sectors responsible for the emission of greenhouse gases must make significant changes. In the power sector, one of the main courses of change is through the development of renewable energy sources (renewables) interfaced by power converters like those for wind and solar photovoltaic energy. The interconnection of these renewable energy sources can be done by high-voltage direct-current (HVDC) power transmission systems, given the big distances involved and the compatibility with the offshore environment challenges [1]. Furthermore, HVDC transmission can be considered as viable reinforcements for the legacy power system as embedded HVDC lines [2,3].
The accelerated development of converter-based devices such as HVDC transmission and converter-interfaced renewable generation is causing a profound change into the very Therefore, stability analysis is of interest given the expected growth in GFM applications in power systems (see [34][35][36][37]). In [34], a stability analysis was conducted for GFL and GFM converters connected to a synchronous generator, one at a time. The study established a comparison between both energy functions and transient stability for GFL and GFM converters, compared through an experimental setup. Ref. [35] presents a design-oriented study for the transient stability of VSCs controlled according to four GFM strategies (droop, droop with power filtering, power synchronization control, and VSM). The VSC was considered connected to an electric grid, and the transient stability was shown through an experimental setup. The study also highlights conditions where the system was rendered unstable. In [36], a transient stability study was performed considering different types of GFM devices. The analysis was performed considering that a GFM with inertia (namely, a VSM or a synchronous generator) is connected to a droop-controlled GFM inverter, showing through simulations that the system is stable. Furthermore, ref. [37] performed a stability analysis of a VSM connected to a synchronous generator using the Lyapunov direct method. The study presented an estimated region of attraction according to different sets of parameters. Additionally, the study presented an analysis of a three-bus system with two VSMs connected to a synchronous generator.
These stability analyses are very important to understand the limits and characteristics of the system. Since MMC-HVDC equipped with GFM control strategies might not be all from the same manufacturer, that implies that a technological difference must be possible given the choice of which GFM strategy to follow. The works that we highlighted in the previous paragraph considered that a single GFM technology is present at each study, to the exception of [36]. In those, a two-elements system was introduced, which was able to represent the case of one generator and one converter or the case of two converters. Moreover, the reviewed studies presented only one equivalent GFM converter in their analysis, to the exception of [37] that considered the interconnection of two GFM converters following the same VSM technology. Therefore, one possible improvement in the subject is the development of analysis on power systems with multiple machines, multiple GFM converters, and multiple GFM technologies.
Hence, to contribute with a wider stability analysis for power systems, including multiple GFM converters and multiple GFM technologies, the scope of this work was to perform a Lyapunov-based stability analysis for a multi-machine/multi-GFM converter system representing the power system during the energy transition. To meet this objective, we review the two most mature technologies for GFM converters under consideration (droop and VSM), highlighting the differences and similarities between them. Furthermore, we introduce a model for the multi-machine/multi-GFM converter system using the center-of -inertia (COI) formulation to establish a single reference for the frequencies and angle deviations in the system. From the model, it is possible to conduct the proposed Lyapunov-based stability analysis using energy functions, demonstrating that the frequencies of the system will converge to the frequency of the COI. To illustrate the theoretical results, a 12-bus test system was implemented in Matlab/Simulink comprising two synchronous machines and two GFM converters. The simulations were performed considering scenarios where there was a single GFM technology and for a scenario where there were multiple adopted GFM technologies. These results allowed the interoperability factor between different GFM converter technologies and synchronous generators to be analyzed. Moreover, since we considered that the converters are representing HVDC transmission systems, a review on MMC control schemes is presented to discuss the challenges of the implementation of GFM converter strategies using the MMC topology.
This article is structured as follows: Section 2 presents an overview of the MMC-HVDC control scheme; Section 3 details the modeling of the multi-machine and multi-GFM system; Section 4 presents the stability analysis using Lyapunov theory; Section 5 details the test system used and the simulation results; and, finally, Section 6 brings the conclusions.

Grid-Forming and Grid-Following Control Strategies for Modular Multilevel Converters
During the early years of power-electronic-interfaced renewables, the stability of the power system was barely affected by the level of nonsynchronous generation. In case of disturbances, transmission system operators (TSOs) would only require the shutdown of the converter-based power plant, which would reconnect only when the disturbance was solved [38]. For this purpose, a simple GFL control strategy to inject active and reactive power to the main grid was required for the grid-connected operation mode for power converters where a PLL structure was used to guarantee that the converter is synchronised with the main power network. In this case, the frequency and voltage output of the converters are governed by the main network, which can be represented most of the time as an infinite bus (from the point of view of the converters).
Because of the growth in the participation of converter-based devices in the electric grid, system operators had to review the policy that required disconnection during the disturbance. In this scenario, if renewables and HVDC lines are simply disconnected, the resulting power imbalance could amplify the magnitude of the disturbance, leading to a blackout event. Hence, grid codes had to be reviewed, including some grid-supporting features for converter-based resources such as a fault ride-through capability and a frequency response capability. Since synchronous generators are responsible for the majority of power generation, the grid still offers an adequate environment for the adoption of the GFL converters with grid-supporting capabilities.
Nowadays, not only are more and more renewables and HVDC projects being commissioned but fossil-fuel-based synchronous generators are being phased out, paving the way to a sustainable power grid. This new power system scenario is already being reported by system operators around the world in regions like Ireland (EirGrid), Great Britain (National Grid), Texas (Electric Reliability Council of Texas, ERCOT), and Australia (Australian Energy Market Operator, AEMO). The shifting generation from synchronous generators to converter-interfaced resources creates a challenging environment for the operation of traditional GFL control structures. Therefore, a new technology concept is emerging in the GFM converters.
A GFM converter is able to regulate its own frequency and voltage. Hence, this control scheme allows for the GFM converter to interact with the grid and to reflect the disturbances in its operation. Moreover, A GFM converter does not rely on the PLL structure for its synchronization with the electric grid (only for measurements in some schemes). These result in a more robust controller for the converter when the power converter is connected to a weak grid with low inertia levels. A drawback for this strategy is that the outer voltage control loop cascaded with the inner control current loop in power converters has complex tuning and can suffer interference from the external parameters of the network, output converter filters, and PLL measurements. Therefore, robust control techniques are required to maintain a suitable operation, regarding the complexity of the power system [39,40].

Classification of Control Strategies for Power Converters
A classification of VSC converters is proposed in [41], based on their different application, such as GFM, including two different approaches of GFL named "grid-feeding" and "grid-supporting." The GFM converters are modeled as an ideal AC voltage source with low output impedance, where the voltage amplitude and frequency are directly controlled, rendering the converter capable of providing services to contribute to the stable operation of the system. In this case, GFM converters need an accurate synchronization to realize power sharing among other generators connected in parallel.
GFL grid-feeding converters are modeled as a current source in parallel with a high impedance to represent the constant injection of current. They are usually designed to deliver constant active and reactive powers to the grid, which is achieved by the control of the currents. In this case, the synchronization with the grid is needed, where standalone operation is not possible. Last, the GFL grid-supporting converters are modeled as con-trolled current sources in parallel with a shunt impedance. They are also modeled as a voltage source. In this case, the currents and voltages of the converters are controlled to improve frequency and voltage levels of the connected grid. Extra ancillary services for grid-supporting may be implemented locally, such as inertia emulation, damping power oscillation, and unbalanced compensation. A general scheme of the three applications of power converters is presented in Table 1 according to [41]. The majority of the research efforts, so far, were carried out to model and control GFM converters with traditional two-level VSCs. For HVDC applications, however, the most prominent topology is the MMC, capable of handling high-power/high-voltage applications. Therefore, in the sequel of this section, we discuss the challenges of the implementation of GFM strategies on MMCs.

Modular Multilevel Converter Model
The model of a 2m + 1 level MMC connected to the electric grid is shown in Figure 1. The MMC is composed of six arms formed by a string of m series-connected submodules. In this work, the half-bridge topology for the submodules was considered (see [42,43] for more details on submodule topology). Each leg of the MMC corresponds to one phase of the converter, and it is comprised of two arms.
A substantial number of works were carried out in the modeling of the MMC for HVDC applications (see [13,[44][45][46][47][48][49][50][51][52]). Given the chosen model, one can design a control system for the MMC-HVDC. The control problem associated with the MMC is more complex than the one associated with two-or three-level VSCs [48,53]. The complexity arises from the growth in the number of variables to be controlled compared to the two-and three-level VSCs resulting from the converter circulating current and the internal energy stored in the submodules of the capacitors. The internal dynamics of the converter has been shown to be asymptotically stable in a non-energy control approach [54], but in this case the dynamics of the circulating currents are given by the system's physical parameters. In [55], it was shown that the use of an energy-controlled approach with direct control of the circulating currents could improve the response of the internal converter dynamics with the design of a proper control system. In the remaining part of this section, we discuss the implementations of the energy-controlled approach in both GFL and GFM MMCs.
The internal-dynamics controller's structure is dependent on the model used for the circulating currents and for the energy dynamics. The internal dynamics of the MMC can be modeled in different frames. In [56], the circulating currents were modeled in the abc frame, in [33,57] in the αβ0 frame, while in [47] the model was carried out in a decoupled double synchronous frame, and the model proposed by [52] was on the dq0 frame. The choice of the frame can influence the complexity of the controller and the accuracy of representation as well. In the simulations carried out further in this work, the dq0 bilinear model proposed by [52] was used for the internal converter dynamics.

Grid Following (GFL)
The control block diagram of a GFL MMC according to the energy-based control approach is shown in Figure 2. The controller can be divided into two subsystems: one responsible for controlling the power exchange between the converter and the grid and the other one responsible for controlling the internal dynamics of the converter. By assuming a balanced system, the AC current controller can be implemented in the dq frame. It is responsible for tracking the current references provided by an upperlevel control by generating the input contribution u AC,dq . The active and reactive power controllers are responsible for calculating the references for the AC current controller from the reference values provided by the TSO.
The circulating current control is responsible for tracking the circulating current references that are determined by an upper-level controller, generating the input contribution u dc,dq0 . When using the bilinear model for the internal dynamics, the reference of the d component of the circulating current can be used to regulate the energy balance ∆W between the upper and lower arms. Moreover, the reference of the 0 component is used to regulate the internal energy W of the converter (see [52,53] for more details on the MMC control using a bilinear dq0 model). These reference calculations were carried out by the energy balance control and the energy control, respectively. The reference for the q component was chosen as zero to diminish power losses, in this case, without loss of generality [53].
To enable grid-supporting capabilities to a GFL converter, a slight modification to the block diagram shown in Figure 3 was necessary. Based on measurements of the grid's frequency ω g and voltage v abc and on the active and reactive powers exchanged with the grid, a contribution for supporting the frequency ∆P or for supporting the voltage ∆Q were obtained in the grid-supporting block and were sent as a correction to the active and reactive power references P * and Q * . Based on this grid-supporting control strategy, it is possible to implement future ancillary services such as virtual inertia and fast frequency response [58]. When the MMC operates according to the GFL strategy (in grid-feeding or gridsupporting modes), the control system needs a measurement of the grid phase angle θ g to establish the dq frame for the AC dynamics and the dq0 frame for the inner converter dynamics. Moreover, the control system also needs a measure of the amplitude of the voltage. This is provided by the PLL block shown in Figures 2 and 3. In strong grids, the PLL exhibits a satisfactory performance as the conditions of the grid change slowly. However, because of the ongoing transformation of the power system, grids are getting weaker and with lower levels of inertia that could turn the PLL into sources of instability [15], limiting the integration of GFL converters in the power system to assure the stability and the reliability of the power supply [19]. As an alternative to increase the penetration of converter interfaced generation, the GFM described in the next subsection was proposed.

Grid Forming (GFM)
Although a precise definition of GFM converters is still not a consensus among the scientific community, a GFM converter can be considered as a converter that can operate and support the operation of a power grid under normal and emergency conditions without relying on services provided by synchronous generators [20]. A GFM can control directly its voltage amplitude, frequency, and phase without relying on the strength of the grid. It can provide services such as (but not limited to) black start capability, virtual inertia, primary frequency support, and islanded operation capability [59].
A block diagram of the control system of a MMC controlled following a GFM strategy is shown in Figure 4. In the GFM concept, the active and reactive power flows are controlled indirectly through the regulation of the AC voltage amplitude and the converter phase angle. When using this strategy, the AC current controller receives its references from an upper-level controller that now targets tracking the d and q components of the PCC voltage, provided by a voltage droop control. By controlling the amplitude of the PCC voltage, it is possible to regulate the reactive power exchange between the converter and the electric grid [60]. The grid-forming strategy block is responsible for determining the phase angle θ g f of the GFM converter and its frequency. The angle θ g f is the integral over time of the frequency ω g f . So, by changing the frequency ω g f , it is possible to control the active power flow by the phase difference between the converter angle and the grid phase angle [60]. This eliminates the need for a PLL for synchronizing with the electric grid.
When using a GFM strategy, the converter frequency ω g f is no longer a frequency value associated to a real rotating mass as is the case for synchronous generators. Hence, in a 100% converter-fed power system, the meaning of frequency ceases to be related to the speed of rotating masses. Additionally, the dq and dq0 frames needed for the system are established by the converter phase angle θ g f in a GFM strategy, restricting the use of the PLL for eventual measurements of frequency.
Depending on the application, a secondary control layer can be embedded in the GFM controller to regulate voltage and frequency to their reference values. Secondary control has the task of providing references for the inner control loops of the primary control. They improve the parallel operation of the generation units and remove steady-state errors from lower control layers. In this way, the frequency and voltage references are calculated by the secondary control, such that, the power sharing and the load supply are accurate [61]. A secondary control loop can be implemented via linear techniques, such that the frequency deviation is eliminated. Therefore, a proportional plus integral term (α ω ) is introduced to improve the convergence speed of frequency and to eliminate the steady state error. The control can be developed as follows: where ω g is the grid frequency, and ω * is the reference of frequency. K ω p and K ω i are the proportional and integral gains of a PI controller and can be allocated according to grid requirements. Therefore, frequency ROCOF and frequency deviation are limited [62].
To assure the elimination of deviations in a steady state and improving the rate of change of frequency (ROCOF), the secondary control term introduced in (1) is inserted in the frequency droop equation given by where ω r is the frequency reference from the droop scheme, m p is the droop coefficient, P m is the measured active power, and P * is the active power reference. The same consideration of the secondary control can be adopted for voltage regulation, where an integral term given by α v can be inserted in the voltage droop. Therefore, the voltage droop control considering α v from the secondary control layer can be written as: where V r is the voltage reference resulted from droop scheme, and V * is the rated voltage reference. Q m is the measured reactive power, Q * is the reactive power reference, and α v is given by: The gains K v p and K v i are calculated according to the desired convergence speed of the secondary control layer [63].
In the following subsections, the two most prominent GFM schemes are discussed.

Virtual Synchronous Machine
The virtual synchronous machine (VSM) provides transient power sharing and primary control frequency support independently, only based on local measurements. A VSM can also be implemented with no need of PLL, being used just for sensing the grid frequency or during initial machine starting. Moreover, the swing equation of the VSM allows the converter to interact with the grid frequency, influencing its behavior. As a result, VSMs are conceptually simple thanks to intuitive interpretation as synchronous machines' responses [22,64,65]. The VSM is implemented to provide voltage reference output where the power flow is related to inertia emulation and the angle from the swing equation, while voltage amplitude and reactive power control are made separately by the modulation index in the converter. The VSM scheme implemented in a point-to-point MMC-HVDC link is introduced in Figure 5. As detailed in [22,23], the VSM is implemented by using the swing equation of a synchronous machine in the VSC control structure aṡω whereω = ω vsm − ω g is the frequency deviation, ω vsm is the VSM's frequency, M v is the virtual inertia coefficient, and D v is the damping factor. P re f is the active power droop reference, and P ele is the measured power into the AC grid. The virtual inertia coefficient is defined in [60]: where S nom is the nominal apparent power of the VSC converter, ω o is the nominal grid frequency value, and J is the emulated moment of inertia.

Droop-Based Topology
The droop control is a well-known approach used for power-sharing without necessary communication among the distributed generation units, which contributes for an easy application. The designed control loop is composed by P − f and Q − V droops considering an electrical grid with inductive impedance (X >> R) and a large amount of inertia, which is the case for conventional power systems with high voltage transmission lines. In the classical case, Equations (2) and (3) model the droop relation. The droop control strategy has only steady-state properties, with no dynamical contribution to frequency or voltage regulation. That can result in a slow transient response with improper active power sharing. In addition, droop control is not able to bring the system back to the original equilibrium point of the system without secondary controllers [66,67]. The droop control scheme implemented in a point-to-point MMC-HVDC link is introduced in Figure 6.  Another approach to provide a dynamical contribution to droop control is to insert a delay in the active power response that emulates the inertial behavior of a synchronous machine [68]. However, the virtual inertial contribution in droop-controlled converters is rather small [6]. In droop control applications, the measured output power is filtered to avoid noise and high-frequency components from power converter switching. Usually, a low-pass filter with a suitable time constant is applied, so the filter will induce a slower behavior in active and reactive power that can be compared with the inertial behavior of a synchronous machine [69,70]. A standard low-pass filter for active power can be described in the frequency domain: where P * out is the filtered output active power measured in the system, and T d is the filter time constant.
According to [65], applying the filter dynamics (7) in the frequency droop Equation (2), we may result in the following expression that presents a virtual inertial component: where R d is the droop coefficient that is typically chosen in the range from 2 to 8%, according to the requirement of the power grid. The resulting equation of the droop control has an equivalent structure to the inertial response (5), which results in an LPF with analogous function of a virtual inertia approach. It is necessary to correctly tune the parameters of the droop regulator to obtain a small signal behavior of a synchronous machine [65].

Modeling of a Multi-Machine, Multi-Converter Power System
Synchronous machines were the foundation of the power system for over 100 years. Because of this fact, their modeling, control, and stability were fairly well investigated. The development of models and energy functions for first-swing transient stability studies for multi-machine systems were thoroughly detailed in [60,[71][72][73][74][75][76][77].
However, the energy transition at hand is changing the foundation of the power system. Growing levels of converter-interfaced generation, such as HVDC systems and renewables, are connected each year to power systems around the globe. Therefore, the modeling, control, and stability of the power converters are now one of the main focuses of attention for researchers (see [2,3,[34][35][36][37][78][79][80]). Specifically, stability analysis using energy functions were provided in [34,37] for single-machine, single-converter systems and in [2,3] for multi-machine, multi-HVDC systems.
In this section, we discuss the model of a multi-machine, multi-HVDC system suitable for performing stability analysis based on energy functions.

Generators and Converters Models
In this work, the classical machine model [60,72] was chosen for the synchronous generators asθ where θ g,i is the rotor angle, ω g,i the frequency, M g,i the inertia constant, P mec,g,i the mechanical power, P ele,g,i the electrical power, D g,i the damping coefficient of the i-th synchronous generator, and ω grid is the grid's frequency.
Assuming that the control of the terminal voltage and of the internal dynamics of the MMC following a GFM strategy is achieved, then the MMC-HVDC can be modeled for the grid stability study in a similar way as a synchronous machine from the grid's perspective. When the VSM GFM strategy is adopted, the converter phase angle and frequency can be modeled asθ where θ v,j is the phase angle, ω v,j the frequency, M v,j the virtual inertia coefficient, P re f ,v,j the power reference, P ele,v,j the electrical power, and D v,j the damping coefficient of the j-th MMC controlled as a VSM-GFM converter. The model for a GFM droop-controlled MMC-HVDC when assuming that the voltage and internal energy dynamics controls were achieved is given by: where θ d,k is the phase angle of the droop controlled GFM converter, ω d,k the frequency, P re f ,d,k the converter power reference, P ele,d,k the electrical power, while T d,k is the power measurement first-order filter time constant, R d,k is the droop coefficient of the k-th MMC controlled as a droop GFM converter, and ω grid is the system's frequency. At this moment, it is possible to trace a comparison between the two GFM control strategies. The dynamic behavior of the GFM converter frequencies given by (12) and (14) presents a similar structure as the synchronous generator swing Equation (10). Therefore, both GFM strategies are able to present some virtual inertia support that is a desirable ancillary service in the future of the power system, as reported in [81]. However, the droop-controlled GFM converter virtual inertia is dependent on the relation T d,k /R d,k of the first-order filter that was added to the measurement of P re f ,d,k − P ele,d,k and the droop coefficient R d,k [6,65,70]. Since the time constant is usually chosen in the order of 10 ms and the droop coefficient in the range of 3 to 8%, the resulting virtual inertia of the droop controlled GFM converter is limited to the order of a tenth of a second, while traditional synchronous generators have inertia constants in the range of 1 to 10 s.

Center-of-Inertia Formulation of a Multi-Machine, Multi-GFM Converter Power System
When performing stability studies using energy functions for a multi-machine, multi-GFM converter power system, it is necessary to establish a reference to compare the rotor and converter phase angles that are essential to assert system stability. A possible reference is the center of inertia (COI). The COI is a notion of the mean motion of the system [73]. This reference is defined by the COI angle θ o , the COI frequency ω o , and the total inertia of the system M T .
The COI for a multi-machine, multi-GFM converter comprising L synchronous generators, M MMC-HVDCs controlled according to a VSM GFM strategy, and N MMC-HVDCs controlled according to a droop GFM strategy is defined in the following: where M T is the total inertia of the system, aggregating physical and virtual inertia. The dynamic of the COI given by (18) and (19) is described as the mean motion of the system's inertia:θ where P COI corresponds to the overall power imbalance on the system. By defining the relative angles to the COI: 20) and the relative frequencies as: we can rewrite the synchronous generator model (9) and (10) and the GFM converter models (11) and (12) and (13) and (14) in the relative variables given by (20) and (21) considering that the grid angle and the grid frequency are given by those of the COI. Therefore, the resulting synchronous generator model in the COI frame is given by: For the HVDC converter controlled as a VSM, the resulting model is given by: M v,jωv,j = P re f ,v,j − P ele,v,j − D v,jωv,j + M v,j M T P COI (25) and for the droop-controlled converter:δ

Lyapunov-Based Stability Analysis Using Energy Functions
As seen in the previous section, a multi-machine/multi-GFM converter is modeled by a set of differential equations. A very resource-consuming method to assert the stability of a huge nonlinear dynamical systems, such as modern power systems, is to solve in time the set of differential equations that models the system and analyze their trajectories. That is a very burdening task to carry out practically. With the advancements in computing capabilities, another way to assess the stability of a dynamic system is to implement it in a software and to carry out exhaustive simulations, trying to span all possibilities and configurations for system operation, analyzing the trajectories of all variables by numerical methods, e.g., Monte-Carlo [82,83]. These methods, on the other hand, are very time consuming and might not lead to a definitive conclusion given the exponential number of combinations arising from the complexity of the systems.
A solution for determining the stability of a dynamical system without solving its differential equations or resorting to a great number of simulations is applying Lyapunov theory. There are two paths on pursuing a Lyapunov-based stability analysis: the indirect method and the direct method [74,84,85]. The indirect method consists of linearizing the nonlinear system around one operation point and obtaining its eigenvalues. If they all have negative real parts, then the system is stable in a neighborhood of the chosen operation point for linearization. In the power systems literature, this approach is referred to as "small-signal" analysis because the conclusions drawn from it are limited to small disturbances that keep the system's trajectory around that operating point chosen for linearization.
Rather than relying upon the linearization around an operating point, the direct method allows consideration of nonlinear systems. Therefore, it can be used for analyzing the stability of power systems facing large disturbances [60,86], which is known as largesignal stability. The stability analysis resides in finding a positive definite Lyapunov function W whose time derivativeẆ is negative definite (or negative semi-definite given some additional conditions). In the analysis of power systems and other physical systems, a Lyapunov candidate can be obtained inspired by the system's energy functions. The use of energy functions for asserting the stability of multi-machine power systems was studied for nearly a century by power engineers as its first reports date back to the 1930s in Russian and to the late 1940s in English [87]. In the next subsection, we introduce the energy functions for the multi-machine/multi-GFM converter power system.

Energy Functions for Synchronous Generators and GFM Converters
For synchronous generators, a possible W g,i (δ g,i ,ω g,i ) function candidate widely used in the literature (see [60,71,73,74,88]) for investigating the stability of a multi-machine power system is given by: The function (28) is called an energy function because the first term represents the kinetic energy of the rotating rotor of the electrical machine, and the second term represents the potential energy of the generator.
For GFM converters controlled according to the VSM strategy, which results in the converter emulating the behavior of a synchronous machine, it is intuitive to propose a similar energy function W v,j , where: For the droop GFM converter, the following energy function is proposed based on the analogies between the model of the GFM converters given by (25) and (27):

Multi-Machine, Multi-GFM Converter System Stability Analysis
When analyzing a system composed of L synchronous generators, M VSM GFM converters and N droop GFM converters, the state of the system can be defined by the relative angles and frequencies as: Then, the following energy function can be proposed: that corresponds to adding the individual energy functions for each of the system components. The time derivative of (33) is given by: where:Ẇ The evaluation of the properties of (35)-(37) is similar. Therefore, the step-by-step procedure is shown only for the first one. By evaluating (35), we obtain: By substituting the model (22) and (23) into the previous equations, the following partial time derivatives are obtained: ∂W g,i ∂ω g,i dω g,i dt =ω g,i P mec,g,i − P ele,g,i − D g,iωg,i + M g,i M T P COI .
Then, the multi-machine, multi-GFM converters, corresponding to the connection of multiple HVDC links with GFM converters to the legacy power system, are asymptotically stable. The analysis incorporated more than one GFM technology, which might be a possibility in the future of power systems, when HVDC terminals will be equipped with GFM control strategies to cope with the challenges of the evolving electrical system.

Simulations and Results
To illustrate the obtained theoretical results in Section 4 for the stability of a multimachine/multi-GFM converter, the test system shown in Figure 7 was implemented in Matlab/Simulink for performing frequency response studies. This system was conceived upon a modification of the IEEE 9 bus test system by exchanging one synchronous generator by an MMC-HVDC and adding an extra bus with power generation from an HVDC system. As a result, the test network had two converter-interfaced resources and two synchronous generators (SG 1 and SG 2 shown in Figure 7).
The MMC-HVDCs considered in this work were modeled as MMC converters connected to an ideal DC voltage source. The ideal DC source was used to represent the HVDC station responsible for the DC voltage control. Both MMC-HVDCs had 500 MW of installed capacity and were controlled according to a GFM strategy. As for the synchronous generators, SG 1 had a 900 MW installed capacity and an inertia constant M g,1 = 5 s, and SG 2 had a 500 MW installed capacity with M g,2 = 4.4 s.
Since the GFM technology is under development, original equipment manufacturers might opt for different GFM control strategies that might better suit a specific challenge or the specifications of their products. Hence, a possible future scenario is that different GFM strategies are connected into the legacy power system during the energy transition period. To reflect these possible future alternatives, the following three scenarios were considered: •  The base dispatch scenario is shown in Table 2. Both HVDC lines were dispatched to provide 300 MW each to the system, while synchronous generator SG 1 power setpoint was 350 MW, and synchronous generator SG 2 power setpoint was 250 MW. There were, initially, three 400 MW loads (Loads 1, 2, and 3) connected to the system, which means that the non-synchronous generation instantaneous penetration in this scenario was of 50%. Then, to simulate a frequency disturbance, Load 4 rated at 450 MW was suddenly connected to the system, corresponding to approximately 38% of the system's load in the pre-disturbance scenario. In the next subsections, the frequency stability of the test network was evaluated for the three specified scenarios when Load 4 was suddenly connected to the system.

First Scenario: Two VSM-GFM Converters
In the first scenario, both HVDC stations connected into the 12-bus test system were controlled as VSMs. The virtual inertia constant for each one was chosen as M v,1,2 = 5 s to match the highest inertia constant of the system (SG 1). Hence, the system presented a total inertia given by (17) of 19.4 s. Compared to the other scenarios analyzed in this work, Scenario 1 was the one with the highest total inertia (physical plus virtual). Figure 8 shows the frequency of each element in the test network when Load 4 was suddenly connected at t = 2 s. It can be seen that for 4 s after the disturbance, the synchronous generators and the GFM converters present damped oscillations in their frequency. Moreover, 5 s after the disturbance, the system as a whole settled at the frequency of the COI. Hence, the test system remained stable after the frequency disturbance caused by the sudden connection of a load equivalent to 38% of the system total load.  Figures 9 and 10 show the behavior of the HVDC terminals controlled according to the VSM-GFM control strategy. In Figure 9, immediately after the disturbance, the frequency reduced in both HVDC terminals, while the power delivered to the grid increased as shown in Figure 10. The converters, as per the control system design, presented the behavior of synchronous machines providing virtual inertia to the system, reducing their frequency while injecting more power into the system to counteract the disturbance.

Second Scenario: Two Droop GFM Converters
In the second scenario, both HVDC stations connected to the test system followed the droop strategy. The GFM strategy parameters corresponding to the time constant of the first order filter were chosen as T d,1,2 = 0.012 s and the droop coefficient as R d,1,2 = 3%. For this scenario, the total inertia of the system given by (17) was 10.2 s. Of all three scenarios adopted, Scenario 2 presented the lowest value of total inertia (considering both physical and virtual inertia).
The frequency of each element of the system after the disturbance corresponding to the sudden connection of Load 4 is shown in Figure 11. The contribution of the virtual inertia brought by GFM converters can be seen in Figure 11, showing that all elements of the system (both synchronous and non-synchronous power sources) had a higher frequency deviation from the reference value of 60 Hz in comparison to Figure 8. Nevertheless, the system remained stable after the disturbance.  Figure 11. Frequency behavior of the elements of the test system after disturbance in Scenario 2.
The frequency and power behavior of the converters connected into the system representing the HVDC systems are shown in Figures 12 and 13, respectively. They both injected power into the grid while their frequency was reduced to counteract the disturbance.

Third Scenario: One VSM and One Droop GFM Converter
Finally, the third scenario aimed at studying the case when the GFM technologies adopted in a system are not the same. The converter station of MMC-HVDC-1 connected into the test system was controlled as a VSM, while the MMC-HVDC-2 station followed the droop strategy. In this case, the virtual inertia constant was also taken as M v,1 = 5 s to match the largest synchronous generator. For the second HVDC, the time constant of the first-order filter was chosen as T d,1 = 0.012 s, and the droop coefficient R d,1 = 0.03%. For this scenario, the total inertia of the system given by (17)  The frequency of the test system following the sudden connection of Load 4 is shown in Figure 14. When comparing to the frequency behavior shown in Figure 11, it is possible to verify that the frequency drop of SG 2 is improved when a larger virtual inertia is present. That is confirmed when comparing to Figure 8 where the frequency of SG 2 presented the smallest frequency nadir considering the disturbance of the sudden connection of Load 4. Additionally, in the mixed GFM technology scenario, the test system remained stable after the disturbance.

The Role of Virtual Inertia
By comparing the results from Scenarios 1 through 3, it is possible to verify in the simulations the contribution of virtual inertia on the frequency support of legacy power grids. Scenario 1 presented fewer oscillations, and the converters were more responsive to frequency disturbances, followed by Scenario 3, where the GFM converter using the VSM strategy contributed more to the frequency disturbance then the droop GFM converter. As some power systems consider procuring physical (and virtual) inertia as ancillary services in a near future, operating a GFM converter with the VSM control strategy can be seen as an advantage.

Conclusions
This work presented a Lyapunov-based stability analysis of a multi-machine/multi-GFM converters power system. A model for the legacy power system permeated with HVDC transmission systems controlled as GFM converters was introduced using the centerof-inertia formulation after a literature review on the most mature GFM technologies set the basis for the modeling of the converters. The stability analysis demonstrates that the frequencies of each element of the system are bound to converge to the frequency of the COI while including the aspect of the interoperability between different GFM converter technologies and synchronous machines.
To illustrate the theoretical results, a 12-bus test system was implemented in Matlab/Simulink to study the frequency response for a disturbance corresponding to a sudden load connection. To reflect the evolution being undergone in power systems, the test network included two synchronous generators and two MMC converter representing HVDC transmission system terminals where the GFM control strategies were implemented. Moreover, the HVDC transmission systems were considered to provide 50% of the power supply.
The simulations performed in this work confirmed the theoretical results obtained with Lyapunov theory showing that the frequency of each element of the system converges for the frequency of the COI. There were three scenarios analyzed: (i) when only VSM-GFM converters are adopted; (ii) when only droop-controlled GFM converters are adopted; and (iii) when mixed GFM technologies are used. From the analysis of these three scenarios, it was possible to verify that, independently of the chosen technology mix, the frequencies converge to the CoI one. The simulation results also show the importance of virtual inertia for the frequency support of power systems. In the considered scenarios, the ones with the highest total inertial levels (physical + virtual) presented the better frequency nadir and fewer oscillations. This was observed because of the higher initial power injection from the VSM-GFM converters from the virtual inertia response.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.