Front. Energy Res.Frontiers in Energy ResearchFront. Energy Res.2296-598XFrontiers Media S.A.10.3389/fenrg.2020.00056Energy ResearchOriginal ResearchSynchronizing Stability Analysis and Region of Attraction Estimation of Grid-Feeding VSCs Using Sum-of-Squares ProgrammingZhangChen1*MolinasMarta1LiZheng2CaiXu31Department of Engineering Cybernetics, Norwegian University of Science and Technology, Trondheim, Norway2College of Environmental Science and Engineering, Donghua University, Shanghai, China3Department of Electrical Engineering, Shanghai Jiao Tong University, Shanghai, China
Edited by: Shuqing Zhang, Tsinghua University, China
Reviewed by: Wei Li, University of Manitoba, Canada; Minh Quan Duong, The University of Danang, Vietnam; Shengjun Huang, National University of Defense Technology, China
*Correspondence: Chen Zhang chen.zhang@ntnu.no
This article was submitted to Smart Grids, a section of the journal Frontiers in Energy Research
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
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.
PLLVSCnon-linearROAstabilitySOSNorges Teknisk-Naturvitenskapelige Universitet10.13039/100009123National Natural Science Foundation of China10.13039/501100001809Introduction
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 small-signal 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, 2016; Sinha 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 “grid-forming” 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 small-scale 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., grid-synchronizing stability (GSS) (Zhang et al., 2017b, 2020; Geng 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. (2017b, 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 non-linear 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, 2020; Hu 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 (Prajna et al., 2002) 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 (Prajna et al., 2002).
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 EAP-based 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., 2017b, 2020; Geng 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.
Single line diagram of a typical grid-VSC system.
The Original GSS Model and Its VariantsThe Original GSS Model
According to Zhang et al. (2017b, 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., Icdqref≈icdqpll) (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
{x∙pll=a1xpll-a2(sinδpll-sinδ0)δ∙pll=a0xpll-a3(sinδpll-sinδ0)
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., Icdqref≈icdqpll). 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
a0=11-kppllĪcdrefL̄g/ωb,kipllĪcdrefL̄sωba0=a1,kipllŪsa0=a2,kppllŪsa0=a3,δ0=sin-1(ĪcdrefL̄gω̄s/Ūs)
In, Īcdref,Ū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 Ūs can be manipulated to emulate the grid voltage sags, which is used in later analysis. xpll and δpll are the states of the PLL, where δpll = θpll − θs and θs=∫ωsdt. Moreover, it can be proven that a0 > 0 holds true for a wide range of PLL parameters; please refer to the Appendices for the justification.
(A) Simplified circuit model for GSS analysis. (B) Mechanism illustration of the GSS.
The First Variant of the GSS Model
Stability analysis of non-linear systems is usually performed at the origin, i.e., x∙=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 = xpll, y = δpll − δ0, for which the resulting model is
{x∙=a1x-a2(sin(y+δ0)-sin(δ0))y∙=a0x-a3(sin(y+δ0)-sin(δ0))
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 x = a0x − a3(sin(y + δ0) − sin(δ0)) into Equation (3), the resulting equation is
{x∙=-(a3cos(y+δ0)-a1)︸d(y)x+b0sinδ0︸F-b0sin(y+δ0)︸f(y)y∙=x
which is the second variant of the GSS model that will be adopted for GSS mechanism analysis. In Equation (4), b0 = a0a2−a1a3 and it can be obtained that b0 is proportional to the magnitude of the grid-voltage Ūs if it is written explicitly with the parameters in Equation (2), whereas F = b0 sin y0 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
{x∙=F-f(y)y∙=x
from which it can be seen that it resembles the classic swing equation of the SG with zero dampings. Therefore, the well-known 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:
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;
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, ymax) can be obtained from the EAP where SI = SII is satisfied.
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., SI ≠ SII, 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:
V(x,y)=12x2-∫0yb0sinδ0-b0sin(ξ+δ0)dξ
which mimics the summation of the kinetic and potential energy of an SG. Also, the level of V is set to the value at Vmax(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 EAP-based 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.
A case where EAP results in an inaccurate ROA (Icdref=1,kscr=2,αpll=10).
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 (Prajna et al., 2002) 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 (Prajna et al., 2002), 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.
Lyapunov Stability and Its SOS ProgramGlobal 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 ⊆ ℝ containing the origin (x = 0) and a continuously differentiable function V, such that V(0) = 0 and
{V(x)>0,x∈D,x≠0-V∙(x)=-(∂V∂x)Tf(x)≥0,x∈D
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: V(x) − l1(x) ∈ SOS, where l1(x) is a non-negative polynomial used to replace the non-polynomial constraint of x ≠ 0. This relaxation V(x) − l1(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):
{x∈D|V(x)≤0,l1(x)≠0}=ϕ{x∈D|-V∙(x)≤0,l2(x)≠0}=ϕ
where l1,2(x) are SOS polynomials used to replace the non-polynomial constraints x ≠ 0. Then, according to the P-satz theorem, this empty-set can be formulated as:
There exists s0,1,2(x), such that
s0-V(x)+l1(x)=0s1+s2(∂V∂x)Tf(x)+l2(x)=0
where s0,1,2(x) are SOS polynomials. Based on Equation (9), corresponding SOS problem can be immediately obtained as follows:
Find V(x) over s2 such that
V(x)-l1(x)∈SOS-s2(∂V∂x)Tf(x)-l2(x)∈SOS
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, s2 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 ∈ ℝ|V(x) ≤ c} as an estimation of ROA and define a subset Hβ: = {x ∈ ℝ|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
{x∈ℝ|h(x)≤β}⊆{x∈ℝ|V(x)≤c}{x∈ℝ|V(x)≤c,x≠0}⊆{x∈ℝ|-V∙(x)>0}
where the P-satz empty-set can be further obtained as:
{x∈ℝ|h(x)≤β,V(x)>c,V(x)≠c}=ϕ{x∈ℝ|V(x)≤c,V∙(x)≥0,x≠0}=ϕ
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 s1, s2, c such that
V-l1∈SOS-(β-h)s1-(V-c)∈SOS-((∂V∂x)Tf+s2(c-V)+l2)∈SOS
where l1,2(x), s1,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: βs1 and Vs2). 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 (Papachristodoulou and Prajna, 2002; Anghel et al., 2013; Shinsaku et al., 2018), e.g., by using the following substitutions:
x1=sin(y+δ0)-sin(δ0)x2=cos(y+δ0)-cos(δ0)x3=x
the GSS model (Equation 3) (i.e., the first variant) can be recast into
x∙1=(x2+cosδ0)·(a0x3-a3x1)x∙2=-(x1+sinδ0)·(a0x3-a3x1)x∙3=a1x3-a2x1
where an additional constraint on x1 and x2 is that
g=(x1+sinδ0)2+(x2+cosδ0)2-1=0
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
{x∈ℝ|h(x)≤β,V(x)>1,V(x)≠1,g=0}=ϕ{x∈ℝ|V(x)≤1,V∙(x)≥0,x≠0,g=0}=ϕ
Then, based on this empty-set, its SOS program can be formulated as:
Maximize β over s1, s2, t1, t2 such that
V-l1∈SOS-((β-h)s1+(V-1)+t1·g)∈SOS-((∂V∂x)Tf+s2(1-V)+l2+t2·g)∈SOS
where l1,2(x) and s1,2(x) are given SOS polynomials which can be declared with the command sossosvar() in SOSTOOLS. The additional functions t1,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 sub-process, 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 s1, s2, t1, t2, t3 are defined as the second group.
To start the first sub-process, solving for s1, s2, t1, t2, β 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 V0 over s2, t1, t3 such that
V(0)-l1∈SOS-((∂V(0)∂x)Tf+s2(β(0)-h)+l2+t2·g)∈SOS
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 s1(0),s2(0),t1(0),t2(0) 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 s1(0),s2(0),t1(0),t2(0), 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.
SOS algorithm for ROA estimation of GSS.
In the algorithm, the following commands are used: (1) declaration of functions: sospolyvar(), sossosvar(); (2) declaration of SOS equalities and inequalities soseq(), sosineq(); (3) solve problems: sossolve(),sossetobj(); (4) return of solutions: sosgetsol(). Other functions and configurations used are: l1,2(x)=10-6∑k=13xk2, h(x)=∑k=13xk2, the tolerance ε = 10−4; s1(x), s2(x) are declared as SOS polynomials and the order of them is 2; t1(x), t2(x) are declared as polynomials with an order of 1. V is declared as a polynomial function [with the condition V(0) = 0] and the order is 2.
Numeric Case Studies and Time-Domain Verifications
This section will present numeric case studies of the SOS-based 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.
Numeric Case Studies of ROA EstimationUnder a Weak AC Grid Condition
In the first place, ROA estimation under a weak grid (i.e., the short circuit ration is kscr = 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
{V(x)=0.312x12+0.0358x1·x2-0.025x1·x3-0.333x22-0.015x2x3-0.002x32β=2.998
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:
VGSS=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.002x2−0.015x(cos(y+δ0)−0.866)
By setting the level of VGSS as VGSS = 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.
ROA comparisons under a weak grid condition (αpll = 10, Icdref=1, kscr = 2).
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 SOS-based 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 SOS-based 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 kscr = 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
{V(x)=0.108x12+0.065x1x2-0.007x1x3+0.261x22-0.006x2x3+0.001x32β=3.8389
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:
VGSS=-0.007x-0.523cos(y+0.2)-0.107sin(y+0.2)+0.083cos(2y+0.0003)-0.006xcos(y+0.2)-0.007xsin(y+0.2)+0.001x2+0.452
where VGSS = 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 EAP-based 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.
ROA comparisons under a strong grid (αpll = 10, Icdref=1, kscr = 5).
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.
(A) Time-domain waveforms of states relevant to GSS model. (B) State trajectories in comparison with ROA estimations.
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.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
ReferencesAnghelM.MilanoF.PapachristodoulouA. (2013). Algorithmic construction of lyapunov functions for power system stability analysis. IEEE Trans. Circuits Syst. I Regular Pap. 60, 2533–2546. 10.1109/TCSI.2013.2246233AshabaniM. (2014). “Synchronous converter and synchronous-VSC- state of art of universal control strategies for smart grid integration,” in 2014 Smart Grid Conference (SGC) (Tehran), 1–8.CespedesM.SunJ. (2014). Impedance modeling and analysis of grid-connected voltage-source converters. IEEE Trans. Power Electron.29, 1254–1261. 10.1109/TPEL.2013.2262473D'ArcoS.SuulJ. A.FossoO. B. (2015). A virtual synchronous machine implementation for distributed control of power converters in SmartGrids. Electr Power Syst Res.122, 180–197. 10.1016/j.epsr.2015.01.001FlourentzouN.AgelidisV. G.DemetriadesG. D. (2009). VSC-based HVDC power transmission systems: an overview. IEEE Trans. Power Electron. 24, 592–602. 10.1109/TPEL.2008.2008441FreitasP.ZuazuaE. (1996). Stability results for the wave equation with indefinite damping. J. Different. Equat. 132, 338–352.GenesioR.TartagliaM.VicinoA. (1985). On the estimation of asymptotic stability regions: State of the art and new proposals. IEEE Trans. Autom. Control. 30, 747–755.GengH.LiuL.LiR. (2018). Synchronization and reactive current support of PMSG based wind farm during severe grid. IEEE Trans. Sustain. Energy9, 1596–1604. 10.1109/TSTE.2018.2799197Gonzalez-LongattF. M. (2016). Impact of emulated inertia from wind power on under-frequency protection schemes of future power systems. J. Modern Power Syst. Clean Energy4, 211–218. 10.1007/s40565-015-0143-xHanG.ZhangC.CaiX. (2018). Mechanism of frequency instability of full-scale wind turbines caused by grid short circuit fault and its control method. Trans. China Electrotech. Soc. 33, 2167–2175. 10.19595/j.cnki.1000-6753.tces.170085HarneforsL.BongiornoM.LundbergS. (2007). Input-admittance calculation and shaping for controlled voltage-source converters. IEEE Trans. Indus. Electron. 54, 3323–3334. 10.1109/TIE.2007.904022HeX.GengH. J.XiGuerreroJ. M. (2020). Resynchronization analysis and improvement of grid-connected VSCs during grid faults. IEEE J. Emerg. Select. Topics Power Electr. [Epub ahead of print]. 10.1109/JESTPE.2019.2954555HuQ.FuL.MaF.JiF. (2019). Large signal synchronizing instability of PLL-based VSC connected to weak AC grid. IEEE Trans. Power Syst. 34, 3220–3229. 10.1109/TPWRS.2019.2892224JohnsonB. B.DhopleS. V.HamadehA. O.KreinP. T. (2014). Synchronization of parallel single-phase inverters with virtual oscillator control. IEEE Trans. Power Electron. 29, 6124–6138. 10.1109/TPEL.2013.2296292JohnsonB. B.SinhaM.AinsworthN. G.DörflerF.DhopleS. V. (2016). Synthesizing virtual oscillators to control islanded inverters. IEEE Trans. Power Electron. 31, 6002–6015. 10.1109/TPEL.2015.2497217KhalilH. (2002). Nonlinear Systems, 3rd Edn. Upper Saddle River, NJ: Prentice Hall.KundurP.PaserbaJ.AjjarapuV.AnderssonG.BoseA.CanizaresC.. (2004). Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions. IEEE Trans. Power Syst.19, 1387–1401. 10.1109/TPWRS.2004.825981LievenV.BoydS. (1996). Semidefinite programming. SIAM Rev.38, 49–95.NguyenT. T.Vu QuynhN.DuongM. Q.Van DaiL. (2018). Modified differential evolution algorithm: a novel approach to optimize the operation of hydrothermal power systems while considering the different constraints and valve point loading effects. Energies11:540. 10.3390/en11030540PapachristodoulouA.PrajnaS. (2002). On the construction of Lyapunov functions using the sum of squares decomposition, in Proceedings of the 41st IEEE Conference on Decision and Control (Las Vegas, NV), 3482–3487.PrajnaS.PapachristodoulouA.ParriloP. A. (2002). Introducing SOSTOOLS: a general purpose sum of squares programming solver, in Proceedings of the 41st IEEE Conference on Decision and Control (Las Vegas, NV), 741-746.RyggA.MolinasM.ZhangC.CaiX. (2016). A modified sequence-domain impedance definition and its equivalence to the dq-domain impedance definition for the stability analysis of ac power electronic systems. IEEE J. Emerg. Select. Topics Power Electron. 4, 1383–1396. 10.1109/JESTPE.2016.2588733SanchezF.CayenneJ.Gonzalez-LongattF.RuedaJ. L. (2019). Controller to enable the enhanced frequency response services from a multi-electrical energy storage system. IET Generat. Transmission Distrib. 13, 258–265. 10.1049/iet-gtd.2018.5931ShinsakuI.SomekawaH.XinX.YamasakiT. (2018). Estimation of regions of attraction of power systems by using sum of squares programming. Electr. Eng.100, 2205–2216. 10.1007/s00202-018-0690-zSinhaM.DörflerF.JohnsonB. B.DhopleS. V. (2017). Uncovering droop control laws embedded within the nonlinear dynamics of van der pol oscillators. IEEE Trans. Control Netw. Syst. 4, 347–358. 10.1109/TCNS.2015.2503558StankovskiT.TiagoP.McClintockP. V. E.StefanovskaA. (2019). Coupling functions: dynamical interaction mechanisms in the physical, biological and social sciences. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 377. 10.1098/rsta.2019.003931656134TanW.PackardA. (2004). Searching for Control Lyapunov Functions Using Sums of Squares Programming. SIBI.TaulM. G.WangX.DavariP.BlaabjergF. (2019). An efficient reduced-order model for studying synchronization stability of grid-following converters during grid faults, in 2019 20th Workshop on Control and Modeling for Power Electronics (COMPEL) (Toronto, ON), 1–7.TeodorescuR.LiserreM.RodriguezP. (2011). Grid Converters for Photovoltaic and Wind Power Systems. Chichester: John Wiley and Sons.WangX.BlaabjergF. (2019). Harmonic stability in power electronic-based power systems: concept, modeling, and analysis. IEEE Trans. Smart Grid10, 2858–2870. 10.1109/TSG.2018.2812712WuH.WangX. (2020). Design-oriented transient stability analysis of PLL-synchronized voltage-source converters. IEEE Trans. Power Electron. 35, 3573–3589. 10.1109/TPEL.2019.2937942ZhangC.CaiX.LiZ. (2017b). Transient stability analysis of wind turbines with full-scale voltage source converter. Proc. CSEE37, 4018–4026. 10.13334/j.0258-8013.pcsee.162359ZhangC.CaiX.LiZ.RyggA.MolinasM. (2017a). Properties and physical interpretation of the dynamic interactions between voltage source converters and grid: electrical oscillation and its stability control. IET Power Electron. 10, 894–902. 10.1049/iet-pel.2016.0475ZhangC.CaiX.RyggA.MolinasM. (2018). Sequence domain SISO equivalent models of a grid-tied voltage source converter system for small-signal stability analysis. IEEE Trans. Energy Conv.33, 741–749. 10.1109/TEC.2017.2766217ZhangC.CaiX.RyggA.MolinasM. (2020). Modeling and analysis of grid-synchronizing stability of a Type-IV wind turbine under grid faults. Int. J. Electr. Power Energy Syst. 117:105544. 10.1016/j.ijepes.2019.105544ZhongQ.WeissG. (2011). Synchronverters: inverters that mimic synchronous generators. IEEE Trans. Indus. Electron.58, 1259–1267. 10.1109/TIE.2010.2048839Appendices
If the PLL bandwidth is defined as the decaying factor of its time-domain response (i.e., αpll = εωn, ε = 0.707),kppll=2αpll,kipll=2αpll2 are obtained. Furthermore, if αpll<ωb2ĪcdrefL̄g=ωbkscr2Īcdref, it obtains a0 > 0. kscr is the short circuit ratio (SCR) of the grid. Usually, this inequality αpll<ωbkscr2Īcdref is satisfied because the critical value of the right-hand side evaluated under an extreme case (e.g., kscr = 2, Īcdref=1) is sufficiently large for a commonly employed PLL, i.e., αpll<αpllcr=ωb normally holds true. As a result, inequality a0 > 0 holds true as well.
Funding. This work was supported in part by NTNU (81617922) and in part by the National Natural Science Foundation of China (51837007).