Synchronizing Stability Analysis and Region of Attraction Estimation of Grid-Feeding VSCs Using Sum-of-Squares Programming

Grid-synchronizing Stability (GSS) is an emerging issue related to grid-feeding voltage-source converters (VSCs). Its occurrence is primarily related to the non-linear dynamics of a type of vastly applied synchronization unit – Phase-locked Loops (PLL). Dynamic characterization and modeling for the GSS analysis can be achieved using a simplified system model, which is a second-order and autonomous non-linear equation, but with the presence of an indefinite damping term. As revealed and demonstrated in this work, this indefinite damping effect can result in an inaccurate Region-of-Attraction (ROA) estimation of the traditional Equal Area Principle (EAP)-based method. To overcome this issue and achieve a valid ROA estimation, this paper adopts the sum-of-squares (SOS) programming technique, which is a numeric optimization method with SOS relaxations. The development and implementation of the SOS program for ROA estimation are presented. Numerical case studies and time-domain verifications demonstrate that this method is valid for GSS analysis, and an almost precise estimation is achieved in the first quadrant. This evidence makes the SOS method a promising tool for GSS analysis because the GSS problem is most concerned with the stability within the first swing cycle, i.e., in the first quadrant.


INTRODUCTION
Voltage-source converters (VSCs) are becoming ubiquitous in electric power systems. According to their functionalities, they can be broadly classified as either grid-feeding and grid-forming. Grid-feeding VSCs are widely employed in bulk power systems for the grid integration of renewable power generations (RPGs) (Teodorescu et al., 2011) and for the interconnection of AC systems through high-voltage-dc (HVDC) transmission techniques (Flourentzou et al., 2009), etc. The primary control objective of a grid-feeding VSC is to process the electric power of different forms into a form that is acceptable to the grid. To achieve this with a high quality, the frequency and phase information of the grid voltage to which the VSC is connected should be detected. This is typically fulfilled with a type of vastly applied synchronization unit−the Phase-locked Loop (PLL). With the help of the PLL, the active and reactive power outputting from the VSCs is decoupled and can be controlled independently; moreover, the frequency of the VSC tends to closely follow the frequency of the grid. Although this fast phase-tracking and frequency-following feature of the PLL is beneficial for the performance of the VSC power controls, it can result in an adverse effect on the transient frequency stability of the power system (Gonzalez-Longatt, 2016). Countermeasures can be taken in this regard; for instance, by augmenting ancillary controls on VSCs (Gonzalez-Longatt, 2016;Sanchez et al., 2019) so that they can contribute to the frequency response of the system. Aside from this issue, the PLL also poses adverse effects on the smallsignal stability of the VSC itself if not properly designed (Zhang et al., 2017a). This issue frequently occurs and is of particular interest under a weak AC grid condition, when its stability impact can be systematically assessed through impedance-based methods (Harnefors et al., 2007;Cespedes and Sun, 2014;Rygg et al., 2016;Zhang et al., 2018;Wang and Blaabjerg, 2019).
The above-exposed adverse impacts of the PLL greatly motivate the advancement and development of non-PLL-based synchronization schemes and associated VSC controls. Among these, the Virtual Synchronous Generator (VSG) control (Zhong and Weiss, 2011;D'Arco et al., 2015) and the Virtual Oscillator Control (VOC) (Johnson et al., 2014(Johnson et al., , 2016Sinha et al., 2017) are very popular. Specifically, the VSG controls eliminate the use of PLL by mimicking the behavior of a Synchronous Generator (SG), whereas the VOC achieves non-PLL synchronization by utilizing the intrinsic features of particular types of non-linear oscillators, e.g., Van der Pol oscillator. With the VSG control, the typical droop controller (Ashabani, 2014) as employed by the SG can be readily implemented in VSCs. The VOC, on the other hand, exhibits some intrinsic droop-like behavior according to Sinha et al. (2017). These features bring a universal synchronization scheme and power-sharing ability to VSCs, which enables them to form and compose an electric system with the absence of SGs, thus giving them the name "gridforming" VSCs. This trait has facilitated the fast emergence and development of VSC-based systems, e.g., smart grids (Ashabani, 2014).
Although VSG controls and VOCs are successful in smallscale and autonomous electric systems, like smart grids, their applicability in bulk power systems is still under discussion. One of the difficulties lies in the fundamental and scientific problem of how thousands of self-organized oscillators behave when they are coupled in circuits and what consequences this will bring (Stankovski et al., 2019). This is an open question in the field of coupled oscillators and, although it is not the topic of this paper, the behavioral characterization and understanding of the diverse non-linear oscillators existing in various types of VSCs seems to be crucial to answering this question. To this end, this paper will dig into a non-linear stability phenomenon associated with the PLL that is emerging in grid-feeding VSCs, i.e., gridsynchronizing stability (GSS) (Zhang et al., 2017b(Zhang et al., , 2020Geng et al., 2018;Han et al., 2018;Hu et al., 2019;Taul et al., 2019;Wu and Wang, 2020;He et al., 2020). Please note that for brevity, the term "grid-feeding" will be omitted in later analysis.
To understand the GSS issue, a concise system model, preferably possessing a certain level of physical implication, needs to be developed. To this end, authors in Zhang et al. (2017bZhang et al. ( , 2020, Hu et al. (2019), and Taul et al. (2019) have developed a simplified model of the grid-VSC system for the GSS analysis, which is essentially a second-order and autonomous nonlinear equation, referred to as the GSS model in later analysis. More appealingly, this GSS model is similar to the classic swing equation of the SG, due to which the GSS mechanism can be qualitatively elucidated according to the theory of SG. For example, the Equal Area Principle (EAP) developed for the SG can be used for this purpose (Zhang et al., 2017b(Zhang et al., , 2020Hu et al., 2019). The obtained mechanism of GSS inspires many studies on stability-oriented control and parameter designs, e.g., the coordinated control schemes for the Low Voltage Ride Through (LVRT) (Geng et al., 2018;Han et al., 2018;He et al., 2020) and the PLL parameter design (Wu and Wang, 2020). The EAP can also help in achieving a quantitative evaluation of the Critical Clearing Time (CCT) of the GSS. However, the accuracy of the predicted CCTs, as shown in Zhang et al. (2020), is dependent on a properly chosen parameter k (a parameter used in Zhang et al. (2020) to estimate the area under the PLL-frequency and along the time axis). The occurrence of this parameter dependence lies in the fact that a successful application of the EAP requires the system to be dissipative. Unfortunately, the GSS model is not a dissipative system because of the indefinite damping term. Therefore, it can be maintained that the parameter k used in Zhang et al. (2020) is essentially a sort of compensation for this indefinite damping effect. However, this is not a formal way to resolve this issue, for which a valid method for GSS analysis which determines its stability region, i.e., its region-of-attraction (ROA) estimation, needs to be developed.
To achieve the estimation of ROA, the search and construction of a Lyapunov function (LF) would be the first step. Unfortunately, there are no generic and systematic methods for constructing LFs. However, an exception is that, provided the system is linear and stable or the linearized system is a Hurwitz system, there exists a quadratic-type LF that can be constructed through several systematic methods (Khalil, 2002). However, this condition (i.e., the system is Hurwitz) is quite strict in practice, which limits its applications, e.g., the linearized system of the GSS model is indeed not a Hurwitz system. Aside from this general difficulty in finding the LF, the search for the LF applied to the GSS analysis is even more difficult because of the indefinite damping effect, which is not a trivial problem in theory either (Freitas and Zuazua, 1996). It can be seen that the analytic construction of the LF is rather difficult, particularly for the GSS analysis; therefore, this paper will turn to the numeric methods (Genesio et al., 1985;Prajna et al., 2002), where the Sum-of-Squares (SOS) programming technique  will be adopted due to its capability in both LF construction and ROA estimation through numeric optimization. Further, these functionalities can be flexibly programmed and implemented in SOSTOOLs, which is a toolbox for MATLAB .
Through the SOS programming technique, the primary objective of this paper is to achieve an improved ROA estimation for the GSS analysis along with a detailed elucidation of its algorithmic implementation. Accuracy of the obtained ROA estimation with the GSS model is verified by time-domain simulations in PSCAD/EMTDC, where a switching VSC model with detailed control systems is employed. To the authors' best knowledge, this is the first attempt to apply SOS programming techniques to solve the GSS problem of grid-feeding VSCs. In addition, the obtained results and algorithmic implementation can be useful references for other advanced applications, e.g., multi-converters analysis in smart grids.

GSS PROBLEM AND THE INDEFINITE DAMPING EFFECT
This section aims at demonstrating the previously mentioned indefinite damping effect and its impacts on the accuracy of EAPbased ROA estimation. To this end, the GSS problem, along with its analysis models, will be first introduced, although in a concise manner as this has been well-examined in existing works (e.g., Zhang et al., 2017bZhang et al., , 2020Geng et al., 2018;Hu et al., 2019;Taul et al., 2019). The study system adopted is shown in Figure 1; it is composed of a three-phase VSC and a Thevenin equivalent grid. The VSC control system consists of a current vector control loop and a PLL. Since most of the grid-feeding VSCs are of high power and connected to a relatively high-voltage power grid, the line impedances are assumed to be purely inductive. This simplification also results in a concise GSS model, as will be shown next.

The Original GSS Model and Its Variants
The Original GSS Model According to Zhang et al. (2017bZhang et al. ( , 2020, Geng et al. (2018), Hu et al. (2019), and Taul et al. (2019), a simplified closed-loop model of the grid-VSC system for the GSS analysis can be formulated as Equation (1), which is the original GSS model. Basically, this model is obtained by assuming that the VSC current dynamics are quasi-steady in the timescale of the PLL dynamics (i.e., I ref cdq ≈ i pll cdq ) (Zhang et al., 2020). This quasi-steady current response can be approximately assured by adopting a relatively large current control bandwidth and compensating the grid-voltage disturbance through the feedforward control. Based on this, the original GSS model is obtained as The circuit representation of this GSS model is given in Figure 2A. By comparing it with the detailed system in Figure 1 it can be seen that the VSC current dynamics are simplified as a current source in the reference frame provided by the PLL (i.e., . The PLL dynamics, on the other hand, are fully considered and modeled. Parameters in Equation (1) are related to the control and circuit parameters of the grid-VSC system, which are listed below In,Ī ref cd ,Ū s are the per-unit VSC d-axis current and Thevenin grid voltage;ω s is the per-unit synchronous angular frequency; ω b = 2π ·50rad/s,ω s = 1; under a normal conditionŪ s = 1, and U s can be manipulated to emulate the grid voltage sags, which is used in later analysis. x pll and δ pll are the states of the PLL, where δ pll = θ pll − θ s and θ s = ω s dt. Moreover, it can be proven that a 0 > 0 holds true for a wide range of PLL parameters; please refer to the Appendices for the justification.

The First Variant of the GSS Model
Stability analysis of non-linear systems is usually performed at the origin, i.e.,ẋ = f (x) with the condition f (0) = 0. This can be achieved by variable substitutions; for example, the equilibrium point(0, δ 0 ) of the original GSS model can be shifted to the origin by using the following substitutions: x = x pll , y = δ pll − δ 0 , for which the resulting model is which is the first variant of the GSS model that will be used in the SOS-based ROA estimation later. Also, it can be readily checked that f (0) = 0 holds true for this model.

The Second Variant of the GSS Model
Although the first variant of the GSS model is sufficient for stability analysis, it gives little insight into the system's behavior. Therefore, a more instructive model for the mechanism analysis can be derived by a simple reformulation, i.e., by substituting which is the second variant of the GSS model that will be adopted for GSS mechanism analysis. In Equation (4), b 0 = a 0 a 2 − a 1 a 3 and it can be obtained that b 0 is proportional to the magnitude of the grid-voltageŪ s if it is written explicitly with the parameters in Equation (2), whereas F = b 0 sin y 0 is a constant value.
Next, based on the second variation of the GSS model, the GSS mechanism will be explained first by ignoring the damping term d y , e.g., how the loss of synchronization of VSC occurs; then, the effect of the damping term d y on the accuracy of EAP in GSS analysis will be further discussed.

Mechanism Analysis of the GSS
First, by neglecting the damping term d y in Equation (4), the following model is obtained from which it can be seen that it resembles the classic swing equation of the SG with zero dampings. Therefore, the wellknown EAP (Kundur et al., 2004) developed for the SG can  be used to elucidate the GSS mechanism of the VSC. This is illustrated as follows. First, characteristics of F and f y can be drawn in Figure 2B, where the stable equilibrium point A is obtained according to the condition F = f y . If the VSC is subjected to a large grid voltage sag, this will lead to a drop of f y i.e., f ′ < f , whereas F remains constant as explained earlier.
Due to the occurrence of the mismatch between f ′ y and F, the following dynamic transition will take place: (1) Initially, since f ′ < F, both x and y start increasing along the characteristic of f ′ y . If the grid fault is cleared at point D, the characteristics of f y will be changed back to their pre-fault state, i.e., f ′ → f ; (2) Then, since f > F, x will decrease whereas y keeps increasing until x = 0, e.g., point C. The value of this condition x = 0, y max can be obtained from the EAP where S I = S II is satisfied.
(3) Afterward, the state x, y will roll back along the characteristic of f y and finally settle down at the original steady-state point (i.e., point A). This is because the system is usually assumed to be dissipative in EAP-based analysis.
In Figure 2B, the trace from (1) to (2) is highlighted with the red line, which is known as the first swing cycle in the context of an SG for transient stability analysis (Kundur et al., 2004). Based on this process, it can be readily obtained that if the EAP is not satisfied when the state reaches the critical point B (i.e., S I = S II , at B), then the VSC will be exposed to the danger of first swing cycle instability, i.e., the loss of synchronization.

Indefinite Damping Effect on the Accuracy of EAP in GSS Analysis
The above analysis explains the GSS mechanism with the help of EAP. Although the obtained result gives insight into the GSS problem, it may not be accurate due to the neglect of the system damping, more importantly, the ignored damping may not be positive as assumed. In fact, the damping term d y observed in the GSS model is indefinite in terms of its sign. This contradicts the assumption of the EAP-based analysis and can give rise to inaccurate results. The following example will demonstrate how this indefinite damping affects the accuracy of the EAP-based ROA estimation of GSS. If the EAP is used for the ROA estimation, its LF can be formulated as: which mimics the summation of the kinetic and potential energy of an SG. Also, the level of V is set to the value at Frontiers in Energy Research | www.frontiersin.org V max (0, π − 2δ 0 ), which is known as the maximum potential energy in the system. Based on Equation (6), the ROA estimation can be obtained and its contour is plotted in Figure 3. Theoretically, this ROA indicates that any trajectories starting within the ROA will be contained within the ROA and finally converge to the origin if the origin is asymptotically stable. To check the accuracy of the EAPbased ROA estimation, two types of trajectories with different initial states are given in the same figure, which are obtained from the dynamic system Equation (4). First, according to trajectory-I in Figure 3, it can be obtained that the estimated ROA is conservative because the initial state of this stable trajectory is not inside its estimation. However, this conservativeness is not consistent throughout the ROA, as indicated by trajectory-II, where the initial state of trajectory-II is inside the ROA but the result turns out to be unstable.
This analysis demonstrates that the EAP can result in misleading results if applied to the ROA estimation of the GSS, and the reason lies in the indefinite damping effect of the GSS model. Therefore, to overcome this issue and achieve an improved ROA estimation, the SOS programming technique  will be adopted and presented in the next section. Other applications of this method, e.g., in the power system analysis, are conducted in Papachristodoulou and Prajna (2002), Anghel et al. (2013), and Shinsaku et al. (2018).

SOS PROGRAMMING TECHNIQUE-BASED ROA ESTIMATION OF THE GSS AND ITS IMPLEMENTATION IN SOSTOOLs
This section will present how the SOS programming technique can be utilized for the LF construction and ROA estimation of the GSS analysis. The developed SOS program will be implemented in SOSTOOLs , which is a toolbox for MATLAB and works as the interface and interpreter for translating formulated SOS problems into semidefinite programs (SDP), where the latter can be solved by various SDP solvers (Lieven and Boyd, 1996). To better illustrate this method applied to the GSS analysis, Lyapunov stability and its SOS program will be introduced first in a generic sense.

Global Stability and Its SOS Program
Considering an autonomous non-linear systemẋ = f (x), the Lyapunov stability theory states that if there exists an open set D ⊆ R containing the origin (x = 0) and a continuously differentiable function V, such that V (0) = 0 and then the origin is stable. If −V (x) > 0 for all x ∈ D, x = 0, then the origin is asymptotically stable. If the limit x → ∞ so that V (x) → ∞ exists, then the origin is asymptotically stable.
In order to find a V satisfying Equation (7) by means of the SOSTOOLs, the non-negativity of those conditions can be relaxed to a problem of finding the SOS polynomials. For example, the first condition of Equation (7) can be relaxed to a problem of finding an SOS polynomial which satisfies: is a non-negative polynomial used to replace the non-polynomial constraint of x = 0. This relaxation V (x) − l 1 (x) ∈ SOS can be declared and solved by the SOSTOOLs with the in-built commands sosineq() and sossolve(), respectively. Here, the second command will call the SDP solvers for solutions. There are various SDP solvers available; this paper uses the STPT3 solver. Another frequently used command for declaring SOS polynomials is sossosvar().
This analysis illustrates that the non-negativity conditions of a semi-algebraic set can be relaxed to an SOS problem that is programmable and solvable. In fact, such relaxations can be formally fulfilled by the Positivstellensatz (P-satz) empty-set and its formulation; please refer to Tan and Packard (2004) for details. With the P-satz formulation, a generic procedure for formulating the SOS problem is shown below: Step 1: Reformulate the original-set into the P-satz empty-set, Step 2: Implement the formulated empty-set in SOSTOOLs and solve for solutions.
For example Equation (7) can be reformulated as the P-satz empty-set as (the asymptotic stability case): where l 1,2 (x) are SOS polynomials used to replace the nonpolynomial constraints x = 0. Then, according to the P-satz theorem, this empty-set can be formulated as: There exists s 0,1,2 (x), such that where s 0,1,2 (x) are SOS polynomials. Based on Equation (9), corresponding SOS problem can be immediately obtained as follows: where Equation (10) can be declared and solved by SOSTOOLs with the above-mentioned commands, e.g., sosineq() and sossolve(). However, this SOS problem is bilinear in decision variables. To overcome this issue, s 2 and V (x) can be solved iteratively; this will be explained in more detail later.

Local Stability and Its SOS Program
The above global stability case demonstrates how an LF can be found using the SOS programming method. However, global stability is rare for non-linear systems; indeed, most of them exhibit local stability behavior around the origin. In this regard, the stability region, i.e., the determination of how far away the initial state can be from the origin while remaining stable, is of great usefulness for study. This can be achieved by finding the largest level for the LF while assuring the −V (x) ≥ 0, i.e., the ROA estimation. To achieve this estimation through SOS programming, the general idea is to first define an estimated ROA; then, define a subset contained in the estimated ROA; afterward, by enlarging this subset, the estimated ROA will be forced to expand due to the set containment; this process will finally result in an optimized estimation of the real ROA of the system. Using the semi-algebraic sets to represent the above-given procedure, the following description is obtained: define a region c = {x ∈ R|V (x) ≤ c} as an estimation of ROA and define a subset H β : = x ∈ R|h (x) ≤ β that is contained in c , i.e., H β ⊆ c and h (x) is a given positive definite function. Therefore, c is forced to expand by enlarging H β . Based on these sets, the SOS program for ROA estimation in a general sense can be developed according to the following steps.
Maximize β such that where the P-satz empty-set can be further obtained as: Based on the procedure illustrated in Global Stability and Its SOS Program, and adding the non-negative condition of the LF, the following SOS problem for ROA estimation in a general sense can be finally obtained: Maximize β over s 1, s 2, c such that where l 1,2 (x), s 1,2,3 (x) are SOS polynomials. Since both V and c are decision variables, c only has a scaling effect on V; therefore c = 1 can be assumed during the implementation. Also, it is noted that this SOS problem is bilinear in decision variables (e.g., due to the multiplication terms: βs 1 and Vs 2 ). The issue can be addressed by using iterations and will be explained later. Besides, maximizing β can be fulfilled by using the bisection method or the in-built command sossetobj() in SOSTOOLs. Aside from this, there are other optimization methods available and one may refer to Nguyen et al. (2018) for extended reading.

SOS Program for the ROA Estimation of GSS
The above SOS programs for LF construction and ROA estimation are introduced in a general sense. When applied to the GSS analysis, some modifications to the program and procedure are required, which will be discussed in this section.

Recast the GSS Model
The first modification required is related to model compatibility.
The SOS program requires the models or systems under study to consist of polynomial vector fields, i.e., f should be a vector of polynomials; however, the developed GSS models consist of non-polynomial vector fields. To make the model compatible with the SOS program, the GSS model can be recast through variable substitutions Anghel et al., 2013;Shinsaku et al., 2018), e.g., by using the following substitutions: the GSS model (Equation 3) (i.e., the first variant) can be recast intoẋ where an additional constraint on x 1 and x 2 is that It is noticed that this additional algebraic constraint assures that the recast system has the same dimension as the GSS model. More importantly, the recast system (Equation 15) consists of polynomial vector fields.

SOS Program for ROA Estimation
The recast GSS model uses polynomial vector fields so that the SOS-based ROA estimation as introduced in Local Stability and Its SOS Program can be applied. However, another modification on the empty-set (Equation 12) should be made due to the additional constraint g = 0 in the recast system. The result is shown below Then, based on this empty-set, its SOS program can be formulated as: Maximize β over s 1, s 2, t 1, t 2 such that where l 1,2 (x) and s 1,2 (x) are given SOS polynomials which can be declared with the command sossosvar() in SOSTOOLS. The additional functions t 1,2 (x) are polynomials and must not necessarily be SOS, and they can be declared in SOSTOOLs by using another command, sospolyvar().

Implementation of the SOS Program in SOSTOOLs
As mentioned above, due to the occurrence of bilinear decision variables, the SOS program in Equation (18) cannot be solved directly through the SOSTOOLs. This issue can be resolved through an iteration process as follows. By dividing the decision variables into two groups, two sub-processes are obtained for solving variables of these two groups iteratively, and for each subprocess, the SOS problem is linear in decision variables and can be solved by the SOSTOOLs. In this paper, decision variables β and V are defined as the first group, whereas variables s 1, s 2, t 1, t 2 , t 3 are defined as the second group.
To start the first sub-process, solving for s 1, s 2, t 1, t 2 , β and V should be initialized. Initialization of β is straightforward: it can be set as a small value, whereas the initialization of V is equivalent to the problem of finding an initial LF satisfying Equation (17). This new problem can be readily tackled if the linearized system is a Hurwitz system. However, this condition does not apply for the GSS analysis where the linearized system of Equation (15) is typically not a Hurwitz system. Therefore, the initial LF will be found through the SOS method as well, where a modified program based on Equation (18) is used and shown below: Find V 0 over s 2 , t 1 , t 3 such that This program is obtained according to the sets containment of Equation (17), from which an initial V (0) can be found. With the obtained initial values of V (0) and β (0) , variables s can be solved from the last two equations in Equation (18); this is the first sub-process. Then, the second sub-process proceeds with the obtained s 2 , and by solving Equation (18), V (1) and β (1) are updated, which are the inputs for the next iteration. This iterative procedure will repeat until the change of β (k) is no longer evident, i.e., ceases if the condition β (k+1) − β (k) ≤ ε is true. Based on this illustration, the complete SOS algorithm for the ROA estimation of GSS is given in Figure 4.

NUMERIC CASE STUDIES AND TIME-DOMAIN VERIFICATIONS
This section will present numeric case studies of the SOSbased ROA estimation of GSS. For comparison, the results of EAP-based ROA estimation [i.e., based on Equation (6)] and the exact ROA [i.e., obtained from the numeric simulation of Equation (3)] are shown as well. Finally, to verify the validity and accuracy of the estimated ROA, time-domain simulations in PSCAD/EMTDC are conducted where a switching VSC is used and the overall system configuration is the same as in Figure 1.

Under a Weak AC Grid Condition
In the first place, ROA estimation under a weak grid (i.e., the short circuit ration is k scr = 2) is analyzed. After running the developed SOS program for ROA estimation with 23 iterations, it returns an LF of the recast system (Equation 15) as where some items with extremely small coefficients are ignored. Then, by substituting Equation (14) into Equation (20), the LF for the GSS model is obtained: V GSS = 0.025x sin y + δ 0 − 0.5 + 0.333 cos y + δ 0 −0.866) 2 + 0.312 sin y + δ 0 − 0.5 2 + 0.036 (cos y + δ 0 − 0.866 sin y + δ 0 − 0.5 +0.002x 2 − 0.015x cos y + δ 0 − 0.866 (21) By setting the level of V GSS as V GSS = 1, the contour indicates the optimized ROA estimation of the GSS under this strong grid  condition. Next, a comparative study of ROA estimations from the SOS, the EAP, and the numeric simulation is conducted and the results are shown in Figure 5.
According to Figure 5, first, it can be seen that the ROA estimation from EAP is inaccurate since part of its boundary is outside the exact ROA. This inaccuracy of EAP in GSS analysis has already been pointed out in Figure 3, which is due to the indefinite damping effect. By contrast, the SOSbased method successfully resolves this issue since its contour is completely contained in the exact ROA. Although the result is still conservative, the estimation is indeed greatly improved with the SOS optimization if compared to the result of the EAP. Moreover, since the GSS analysis is most concerned with the stability of the first swing cycle, from this perspective, the SOSbased ROA estimation is rather satisfactory since its estimation is similar to the exact ROA in the first quadrant.
Besides, one may notice that the shape of ROA with EAP is different from that in Figure 3. This is because in Figure 3 the estimated ROA is based on the second variant of the GSS model, i.e., Equation (4), whereas this analysis is based on the first variant of the GSS model, i.e., Equation (3).

Under a Relatively Strong AC Grid Condition
Next, the ROA estimation under a relatively strong grid will be analyzed, where the short circuit ration of k scr = 5 is used. Similarly, the SOS is run in SOSTOOLs, and after 37 iterations, it returns an LF of the recast system (Equation 15) as where the items with extremely small coefficients are ignored. Then, substituting Equation (14) into Equation (20) yields the LF for the GSS model, which is: V GSS = −0.007x − 0.523 cos y + 0.2 − 0.107 sin y + 0.2 +0.083 cos 2y + 0.0003 − 0.006x cos y + 0.2 −0.007x sin y + 0.2 + 0.001x 2 + 0.452 (23) where V GSS = 1 indicates the optimized ROA estimation under this strong grid condition. Then, a similar comparative analysis to the first case study is conducted and the results are presented in Figure 6. According to this figure, again, it is immediately concluded that the EAPbased ROA estimation is inaccurate since part of its boundary is outside the exact ROA. The SOS-based method, on the contrary, is well-contained in the exact ROA, indicating a correct but conservative ROA estimation. Moreover, it can be seen that in terms of the GSS analysis, where the focus is on the stability of the first swing cycle, the SOS-based method is rather accurate since its ROA estimation is similar to the exact one if compared in the first quadrant.
On the other hand, interestingly, when comparing the result of the ROA estimation under the strong grid (i.e., Figure 6) to that under the weak grid (i.e., Figure 5), it can be seen that the conservativeness of the ROA estimation tends to improve when the grid becomes strong.

Time-Domain Simulation and Verification
To verify the numeric results, particularly the validity and accuracy of the SOS-based ROA estimation, time-domain simulations are conducted in the PSCAD/EMTDC. To provoke the GSS issue, the VSC in simulation is perturbed by a large grid voltage dip. In detail, the magnitude of the Thevenin grid voltage is dropped down to 0.2 p.u. at 1 s. In addition, two types of fault clearing time (FCT) are considered: the shorter one is 135 ms and the longer one is 140 ms. Corresponding time instants at which the fault is cleared are marked by points A and B in the following figures.
First, in Figure 7A, time-domain waveforms of states relevant to the GSS model under two types of FCTs are presented. It can be seen that under a longer FCT, i.e., 140 ms, the states are unstable after the fault is cleared, whereas under a shorter FCT, i.e., 135 ms, the system is stable after the fault is cleared. This simulation implies that the state trajectory of the GSS model at point B (i.e., under FCT = 140 ms) is anticipated to be outside the ROA, whereas the one at point A (i.e., under FCT = 135 ms) should be inside the ROA.
In order to check if the same conclusion can be drawn from the estimated ROA, the same state trajectories under two types of FCTs are plotted together with the ROA estimations from the SOS and EAP, which are shown in Figure 7B. From this figure, it can be seen that point B is indeed located outside the ROA of the SOS method, while point A is inside its ROA estimation. This analysis demonstrates that: (1) the SOS-based ROA estimation is correct and quite accurate since the difference between the longer and shorter FCT is very small; (2) the simplified model for the GSS analysis is feasible and valid, otherwise the estimated ROA would be inaccurate. On the other hand, this time-domain study again demonstrates that the EAP-based ROA estimation is inaccurate and leads to a misleading stability result, as can be seen from Figure 7B, where point B is inside the estimated ROA while the time-domain result turns out to be unstable.

CONCLUSIONS
This paper discussed the GSS issue and its ROA estimation of grid-feeding VSCs. First, the indefinite damping effect in GSS analysis is revealed and its impact on the accuracy of the EAP-based ROA estimation is discussed. Then, the SOS programming technique is employed to address this indefinite damping issue and to achieve an improved ROA estimation, in which the SOS program for ROA estimation of GSS is developed and its implementation in SOSTOOLs is elaborated. Finally, comparative numeric studies of the ROA estimations with different methods are performed. The obtained numeric results are verified by time-domain simulations in PSCAD/EMTDC.
Based on the presented analyses, it can be concluded that the SOS-based method can successfully overcome the indefinite damping issue encountered by the EAP-based method. Although the estimated ROA is still conservative, it is nonetheless improved through SOS optimization. Particularly, the estimated ROA is similar to the exact ROA in this first quadrant. This evidence demonstrates that this method is even more useful for GSS analysis because the GSS problem is most concerned with the stability of the first swing cycle, i.e., the first quadrant of the estimated ROA.
Besides, the presented methodology and results can be useful references and can lay the foundation for other advanced applications in the future, e.g., smart grids with multi-converters.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
CZ and MM contributed to the conception and design of the study. ZL and XC organized case studies. CZ wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This work was supported in part by NTNU (81617922) and in part by the National Natural Science Foundation of China (51837007).

APPENDICES
If the PLL bandwidth is defined as the decaying factor of its time-domain response (i.e., α pll = εω n , ε = 0.707),k ppll = 2α pll , k ipll = 2α 2 pll are obtained. Furthermore, if α pll < ω b 2Ī ref , it obtains a 0 > 0. k scr is the short circuit ratio (SCR) of the grid. Usually, this inequality α pll < ω b k scr 2Ī ref cd is satisfied because the critical value of the right-hand side evaluated under an extreme case (e.g., k scr = 2,Ī ref cd = 1) is sufficiently large for a commonly employed PLL, i.e., α pll < α cr pll = ω b normally holds true. As a result, inequality a 0 > 0 holds true as well.