Conditions of Existence and Uniqueness of Limit Cycle for Grid-Connected VSC With PLL

This paper investigates the large disturbance stability of single voltage source converter (VSC) connected to the infinite bus system through the transmission line. One conservative stable region is obtained via the direct method of Lyapunov. Based on the Poincaré Theorem, it is interesting to find there might exist limit cycle between the obtained stable boundary and the angle of unstable equilibrium point (UEP). Accordingly, the nonlinear damping effect on the system stability is defined, and one limit cycle exists when the positive and negative damping effect over one period is mutually balanced. A conservative existence condition of limit cycle is further derived by constructing a proper piece-wise curve to approximate one specified unstable system trajectory. The derived analytical existence condition can directly predict the limit cycle without numerical simulations. Finally, the paper also proves that the studied dynamic system has at most one limit cycle, and if it exists, it is unstable.


Conditions of Existence and Uniqueness of Limit
Cycle for Grid-Connected VSC With PLL Yujun Li , Member, IEEE, Yiyuan Lu , Yingyi Tang , and Zhengchun Du , Senior Member, IEEE Abstract-This paper investigates the large disturbance stability of single voltage source converter (VSC) connected to the infinite bus system through the transmission line.One conservative stable region is obtained via the direct method of Lyapunov.Based on the Poincaré Theorem, it is interesting to find there might exist limit cycle between the obtained stable boundary and the angle of unstable equilibrium point (UEP).Accordingly, the nonlinear damping effect on the system stability is defined, and one limit cycle exists when the positive and negative damping effect over one period is mutually balanced.A conservative existence condition of limit cycle is further derived by constructing a proper piece-wise curve to approximate one specified unstable system trajectory.The derived analytical existence condition can directly predict the limit cycle without numerical simulations.Finally, the paper also proves that the studied dynamic system has at most one limit cycle, and if it exists, it is unstable.Index Terms-Grid-connected VSC, large disturbance stability, PLL, existence and uniqueness, limit cycle.

I. INTRODUCTION
K EEPING synchronization is one necessary and critical condition for any ac electric power system.Traditional synchronous generators (SGs) rely on the natural power-angle property to maintain the synchronization with the power grid during the large disturbance [1].However, the converter interfaced generators (CIGs) depend on the PLL [2] to track the grid phase, which is achieved by sustaining the q-axis voltage of the point of common coupling as zero during the disturbances.This control-based synchronization system may lead to the dynamic behaviour of nowadays power grid varying a lot from the traditional one and bring about some new instability phenomena.It has been widely found that CIGs may lose the synchronization under the large disturbance due to PLL's incapability of tracking the grid phase [3], [4], [5].

A. Related Works
The study of the synchronization stability of CIGs embedded power system under the large disturbance is a tedious and hard The authors are with the School of Electrical Engineering, Xi'an Jiaotong University, Shaanxi 710049, China (e-mail: yujunlizju@gmail.com;luyiyuan961015@stu.xjtu.edu.cn;tangyingyii@163.com;zcdu@xjtu.edu.cn).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TPWRS.2023.3238000.
Digital Object Identifier 10.1109/TPWRS.2023.3238000 task since the mathematical model of such dynamic system involves high-order nonlinear differential algebraic equations [6].
The general way to analyze the stability of high-order nonlinear system is to develop the time response via numerical simulations and to check whether CIGs tend to lose or to maintain synchronization with the power grid.Accordingly, the stable region of the converter with PLL under the specified SEP can be obtained [7], which can provide quantitative stability information for the system operator [8].Although the results obtained from the time-domain simulation are reliable and accurate, this method involves complex and time-consuming calculations and cannot often be effectively applied in stability analysis.
To determine the main factors that influence the studied system stability, one reduced second-order model is firstly proposed in [9] that the fast inner-current control of VSC and the network dynamics are reduced into the algebraic equations and only the dynamic equations of PLL are reserved.Phase portrait method that maps the VSC velocities and angles information during the transient process in the phase plane is widely utilized in the stability analysis of VSC embedded system [10].The stability of the system can be clearly judged to see whether the system trajectory is escaping from SEP or not.Moreover, the geometry property of the stability boundary can be obtained from this method [11].The large disturbance stability of the modular multilevel converter based wind farm considering the current limiting mode using the phase plane method is reported in [12].Even though the two-order nonlinear system can be identified using phase plane method, it still relies on the numerous calculations to obtain the system trajectory, and fails to provide the analytical results for the studied nonlinear system as well.
Equal area criterion (EAC) and Lyapunov Theorem are two main methods to obtain the stable region of the nonlinear system from the analytical way.[13] firstly utilizes the classical EAC into the synchronization stability analysis of CIGs based on the similar dynamics of PLL and the rotor motion of SG.EAC converts the stability analysis of two-order nonlinear system into the comparison of the accelerating and decelerating areas, which can provide a clear physical interpretation of loss of the synchronization of CIGs.However, it can only be feasible to study the nonlinear system with no damping term involved, such as pendulum equation.
Lyapunov method [14] is deemed as one effective direct method to analyze the nonlinear system, and it emphasizes the idea of the relative stability of the system without developing the time response.More importantly, the analytical conservative stable region can be obtained using this method.For the traditional ac system, [15] has proved that EAC and Lyapunov method are identical for two-machine system.Proper Lyapunov functions are found for the grid-forming converter in [16] and [17] based on the similar dynamic equations of the converter and SG.However, as for the PLL-based converter, finding one proper Lyapunov function is rather difficult due to the nonlinear damping term involved in its dynamic model.One conservative stable region is obtained for PLL-based converter based on the Lyapunov Theorem in [18], and a larger stable region with less conservation is provided via constructing a more suitable Lyapunov function in [19] later.In order to meet the constraints of Lyapunov Theorem, the estimation of domain of attraction for PLL-based converter is transformed into the semi-definite programming in [20].Lyapunov Theorem normally requires that the system damping term should maintain positive so that the derivation of Lyapunov function with respect to time is non-positive.As a result, the obtained stable regions are only the subsets of the exact stable region.

B. Main Contributions
References [13], [16], [17], [18], [19] consider that the system damping always presents as positive.Thus, only the classical aperiodic instability phenomenon appears, and it cannot completely reflect the real dynamic behavior of the studied systems.In this paper, the nonlinear damping effect, especially the negative damping effect on the system stability is fully considered.Phenomenon of the limit cycle is discovered, and the exact stable region of the grid-connected VSC with PLL is composed by one unstable limit cycle, instead of the manifold of UEP.Based on the Poincaré Theorem, limit cycle might exist between the obtained conservative stable region from Lyapunov method and the angle of UEP.Limit cycle appears when the positive and negative damping effects over one period are counteracted.
One analytical existence condition of the limit cycle for the grid-connected VSC with PLL is derived by constructing one piece-wise curve to approximate an unstable system trajectory.The conservation of the derived condition is guaranteed since the energy change of the unstable system trajectory over one period is always larger than that of the constructed curve.The derived analytical condition can easily predict the existence of limit cycle of the studied system without numerical simulations.Finally, the uniqueness of the limit cycle is also proved in the paper.

II. SYSTEM MODELLING
The studied power system consists of one VSC connected to the infinite bus through the transmission line.As shown in Fig. 1

A. Dynamic Modelling
As shown in Fig. 2, d-q and x-y reference frames are counterclockwise rotating with the angular speed of ω and ω s .d-axis leads x-axis with θ, which is provided by PLL.
where U sL is the rms value of the line-to-line voltage of u abc s .The aim of PLL is to sustain u q p as zero so that the phase angle of u p can be tracked during disturbances.As shown in Fig. 3, the dynamic equations of PLL are, where K PLL P and K PLL I are the proportional and integral gain of PLL.
Generally, the inner current control dynamics of VSC are much faster than that of PLL.When studying the synchronization stability of VSC with PLL, it is regarded that the fast state variables of the inner current control marked as i d c and i q c are reduced to their steady-state value [11].

B. Stable Equilibrium Point of System
One equilibrium point of the final system after disturbance can be directly obtained from (5) by letting the derivations of the state variables as zero.
Equation ( 5) is the classical Liénard Equation.Using the following transformation, (5) can be transferred into the generalized form, where f (x) denotes the system damping property, and g(x) represents the system restoring force property.The equilibrium point of the system (7) is the origin.To make sure the origin is SEP of the studied dynamic system, the small signal stability of the origin should be satisfied.Accordingly, the linearized model of (7) around the origin is deduced as follows. where As a result, the characteristic function of ( 8) is obtained.
The condition to make sure the small signal stability of the origin point is as follows.
Inequality of (10) relates to two conditions to ensure the system small signal stability.
where u d s(0) = U sL cos θ s relates to the initial value of u d s .Combining (11) and (12), the range of the proportional parameter of PLL for ensuring system small signal stability is derived as follows.
It is founded that both large and small proportional gains of PLL may lead to the system instability.The larger of K PLL P or the smaller of K PLL I , the larger of the system damping coefficient f (0), and the better dynamic performance of PLL controller will be achieved.Inequality of (13) ensures the system small signal stability, which is the pre-condition for analyzing the large disturbance of the grid-connected VSC.

A. Energy Function and Obtained Stable Region
Generally, Lyapunov function (energy function) is chosen as a quadratic form in the variables of the system or in many cases as a quadratic form plus an integral.One proper energy function is defined as follows.

E(x, y)
with where θ I = π − 2θ s is the angle of UEP of the studied dynamic system.Energy function of ( 14) is positive definite over the interval specified by (15).The classical aperiodic instability appears once the system trajectory exceeds the angle of UEP.The format of defined energy function has clear physical meanings where the first term of ( 14) is similar to the kinetic energy of system, while the second term is related to the potential energy of the system.
Based on the Lyapunov theorem, the derivation of the defined energy function with respect to time should be non-positive for ensuring system stability.
It implies, where Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
It can be inferred from Inequality of ( 13) and ( 17) that there is one trade-off of proportional gain of PLL in the perspective of system stability.Large gain of PLL may achieve better dynamic performance of PLL and enlarge the stable region of the studied system, but might also damage the system small signal stability.
Based on the conditions of ( 15) and ( 17), the conservative stable region of the studied system obtained by the direct method of Lyapunov is provided by (17).Inequality of ( 17) has clear physical meaning.If the nonlinear damping coefficient f (x) always keeps positive after the disturbance, the system stability can be ensured.
The range of value for E(x, y) are from 0 at origin to the maximum at x = θ II − θ s , y = 0 marked as E max , that is, The stable region noted as Ω by Lyapunov method is provided as, It is indicated that every trajectory starting in Ω stays in Ω for all future time.

B. Non-Existence of Limit Cycle in the Stable Region
Theorem 1 (Poincaré) [21]: has constant sign, and the curve dE/dt = 0 does not contain the whole trajectories of ( 7), then the system does not possess a closed trajectory which is wholly contained in G.
Based on Theorem 1, E(x, y) is continuously differentiable function.In the conservative stable region Ω by (19), dE/dt ≤ 0 always holds.In other words, in the region Ω, the quality of dE/dt has constant sign.Moreover, the curves dE/dt = 0, i.e. f (x) = 0 are two lines of x = −θ s − θ II and x = θ II − θ s .They do not possess a closed trajectory that is inside in Ω.As a result, the studied system has no closed orbits in the stable region mapped by Ω. Further, all trajectories in this region converge to the origin.
On one hand, Theorem 1 indicates in the stable region derived by Lyapunov method, there does not exist limit cycle.It also implies that there might exist closed orbits in the region of {(x, y)|E(x, y) > E max }, where the quality of dE/dt cannot keep the constant sign.
On the other hand, aperiodic instability happens when any system trajectory crosses over the angle of UEP.It implies that limit cycle does not exist outside the angle of UEP.Therefore, there might exist closed trajectories (orbits) outside the stable boundary and within the angle of UEP, namely,

C. Nonlinear Damping Effect and Limit Cycle
For any system trajectory ( P Q) from Point P to Point Q, the following Equality is naturally held based on ( 16).
where E P and E Q are the energy of Point P and Point Q.The energy change noted as E P Q for the system trajectory ( P Q) is defined as As shown in Fig. 4, Γ is one closed orbit (limit cycle) of the grid-connected VSC.Line AB is the normal line of Γ at any Point P 0 , and the positive direction is selected as outside.The system trajectory starting from Point Q 0 (on the Line AB) will intersect the Line AB at Points of As a result, the system trajectory is one closed orbit (limit cycle), and T is the period of the limit cycle (period of the oscillation).
As shown in Fig. 5, for one closed orbit Γ, the following equality should be held.
Accordingly, for one system trajectory Γ over one period, if Γ (dE/dt)dt > 0, the system trajectory will escape from SEP.Similarly, the system stability is maintained and the system trajectory will finally converge to SEP when Γ (dE/dt)dt < 0 Moreover, the limit cycle over one specific period T can be divided into two separated parts as shown in Fig. 5.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.where It is indicated in Fig. 5 that in the area of E − , the damping coefficient f (x) > 0 is always held, and the total energy is decreasing with time.As for the area of E + , f (x) < 0, and the system energy is increasing.E − and E + have apparent physical meanings, which denote the positive and negative damping effect on the system stability.If E − + E + = 0, the system positive and negative damping effects over one period are mutually balanced on the limit cycle.Specifically, if the negative damping effect exceed the positive damping effect over one period, the system stability cannot be maintained and the oscillatory instability will occur.Moreover, Reference [22] indicates that the current control loop may provide the negative damping and deteriorate the system stability in some cases.However, Equality (24) holds if the negative damping induced by the current control loop is considered, and the limit cycle still happens only if the positive and negative damping effect over one cycle are mutually balanced.

IV. CONDITIONS OF EXISTENCE AND UNIQUENESS OF LIMIT CYCLE
Two fundamental problems regarding the limit cycle of system are discussed in this section.In what conditions, the limit cycle will appear, and how many limit cycles exist for such dynamic system, which relates to the existence and uniqueness of the limit cycle.The issue of existence of limit cycle is considered first.

A. Analytical Condition of Existence of Limit Cycle
Suppose that the studied system exists one limit cycle, there exists one divergent system trajectory AG lying outside the neighbourhood of limit cycle as shown in Fig. 6.In the studied region of (−π − 2θ s , 0) and (0, π − 2θ s ), if the following inequality holds, there has one unstable limit cycle for such Fig. 6.Existence of limit cycle of the grid-connected VSC.dynamic system.
As shown in Fig. 6, the unstable system trajectory over one period ( AG) is divided into three separated arcs, namely arc ABC, arc CDE, and arc EF G. x (26) Point C, and Point F are the two farthest points that can be reached along x-axis on the curve ( AG).Accordingly, (25) can be re-written as follows.
The value of E EG is expressed below.
In the region of x F < x < x G , f (x) < 0 holds.Based on (28), E EG > 0. As a result, if Inequality of ( 27) is satisfied.E AC is expressed as follows.
(30) Based on the relation of ydt = dx, E AB and E BC can be re-written as follows.
As shown in Fig. 7, one piece-wise curve in red solid is constructed that includes Line A B and Line KC to approximate Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the arc ABC as follows.
where y max is the maximum value of y for the arc AB.Point M is the highest point of arc AB, and Point N is the intersect point of Line A B and y-axis.Point K is the intersect point of Line NC and Line x = x B marked as dashed gray line.On the Point A and Point B, it implies, The arc AB is continuous in the range of x A < x < x B , and there definitely exists one point M that satisfies dy/dx = 0. Accordingly, Point M is in the range of x A < x < x B .Based on (31) and (33), Due to the convexity of the system trajectory in the phase plane, Line BC is below the arc BC as shown in Fig. 7.Moreover, the slope of Line NC is larger than that of Line BC.
Therefore, Line KC is below the arc BC.Based on (32) and (33), Combining (34) and ( 35), Accordingly, if E AC > 0 is held.Using the similar technique, Inequality (37) also make sure E CE > 0. As a result, (29) is ensured, and there exists one limit cycle for the studied dynamic system.The derived condition of (37) is a function of x C , and the Inequality (37) multiplying x C is defined as H(x C ), which is written as follows.
It can be seen that the defined function H(x C ) is first decreasing and then increasing in the range of (38) does not hold when x C = x B .Therefore, the analytical condition of the existence of limit cycle in the range of x B < x < θ I for the grid-connected VSC is obtained as follows.
Inequality (40) builds up the analytical relationship between SEP and control parameters of PLL for the existence of limit cycle of the grid-connected VSC.
Firstly, to ensure that the studied system has SEPs, Inequality of (13) should be satisfied.It is indicated that the larger or the smaller proportional gains of PLL might lead to system instability.In addition, to enlarge the system stability region so that the limit cycle does not appear, the value of H(x UEP ) should be as small as possible, that is, H(x UEP ) ≤ −a, a is one positive real number.These Inequalities provide the constraints of controller parameters for ensuring system stability, which gives one useful guidance for the controller design and parameters optimization of VSC.
2. There exists real number α and β such that the function , the function f 1 (x)/g(x) does not increase.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
4. All the cycles contain the interval [x 1 , x 2 ] on the x-axis.
Then system (7) cannot have more than one limit cycle; if one exists, it must be an unstable cycle.
In the studied specified range of (−π − 2θ s , 0), and (0, π − 2θ s ), xg(x) > 0 always holds.Based on (10), f (0) > 0 is satisfied for making sure the origin point is SEP.The condition 1 of Theorem 2 is guaranteed.As a result, For the condition 2, two real numbers α and β are set as zeros in this paper, and f 1 (x) = f (x) has two zero points in the range of (−π − 2θ s , 0), and (0, π − 2θ s ), which indicates, is the obtained stable region based on (17), and f (x) ≥ 0. Accordingly, the condition 2 of Theorem 2 is satisfied.Based on Theorem 1, closed orbit does not exist in the region of [x 1 , x 2 ], and the condition 4 of Theorem 2 is naturally held.For the condition 3 of Theorem 2, should be non-increasing outside the region Inequality of (43) is equivalent to The maximum of the function d(x) is Based on (41) and (45), The condition 3 of Theorem 2 is ensured.As a result, all the conditions of Theorem 2 are satisfied, and it is proved that for the studied dynamic system, it has at most one limit cycle, and if it exists, it is an unstable limit cycle in the studied range of (−π − 2θ s , 0), and (0, π − 2θ s ).Moreover, if the Inequality (40) is satisfied, the studied system has one unique unstable limit cycle.

TABLE I PARAMETER OF GRID CONNECTED VSC
In fact, the analytical prediction of limit cycle holds for the studied system that one single VSC connected to the infinite bus system through transmission line.For multi-VSCs embedded power system, it is rather tough to obtain the stability boundary due to the mutual interaction between each VSC.One possible but rough approach is to transfer the original complex multi-VSCs embedded system into several independent subsystems that single VSC connected to the infinite bus system via one equivalent impedance.The parameter of the equivalent impedance can be derived by Thevenin's theorem.Accordingly, the prediction of limit cycles might be possible in the multi-VSCs embedded system.Certainly, the errors from the network decoupling are hard to be quantified so that the conservativeness of this method becomes unknown.However, the proposed approach can be deemed as one possible heuristic method for predicting limit cycles in the multi-VSCs embedded system.Actually, it still remains big challenge to explore the effective approaches to study the large disturbance stability of the multi-VSCs embedded power system, especially taking into account the mutual interaction between each VSC during disturbances.

V. NUMERICAL RESULTS
The case of one VSC connected to the infinite bus through one transmission line is utilized to validate the uniqueness and the derived existence condition of limit cycle in Matlab/Simulink.Base capacity of the studied system is 200 MVA, and the nominal system voltage is 220 kV.The equivalent line resistance and inductance are 0.05 pu and 0.5 pu.PI parameters of PLL are set as 0.028 and 20.The delivered power from VSC is 110 MW, and the d-axis current reference i dref c is 0.50 kA.Accordingly, θ s is 0.2788 rd based on (6).The other related parameters for the case study can be referred to Table I.The disturbance is that the infinite source U sL suddenly drops to 0.3 pu at time t = 0.01 second, and returns to U sL after t c seconds.

A. Model Verification
Fig. 8 shows the dynamic behaviour of three models, namely the detailed EMT model with inner current control and outer power control, the proposed model by (7) and the existing model with no nonlinear damping under the different t c .As indicated by Fig. 8(a), if the fault is cleared in time where t c = 55 ms, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the system keeps stable under three models.Considering the effect of f (x)y, the proposed model presents a damped oscillation similar to the detailed model, while the existing model exhibits the constant amplitude oscillation.For a longer fault duration time (t c = 88 ms) in Fig. 8(b), the oscillatory instability phenomenon appears in the proposed model and the detailed model.Moreover, this phenomenon disappears and the classical aperiodic instability happens when the fault duration time is increased to t c = 92 ms as shown in Fig. 8(c).However, the system keeps stable in the existing model under different t c .It can be inferred that the proposed model is well suitable for the large disturbance stability analysis of VSC.

B. Stable System Dynamic Behaviour
Fig. 9 indicates the system dynamic behaviour under t c = 33 ms.It is clearly seen from Fig. 9(a) that the system trajectory converges to the origin and the system stability can be ensured.The upper stability boundary with dashed gray line shown in Fig. 9(a) calculated from ( 17) is x = 0.6161.The farthest point of the system trajectory along x-axis, noted as Point A (0.5669, 0) is inside the obtained upper stability boundary, which well verifies the effectiveness of the stabilizing condition from (17) based on the direct method of Lyapunov.The derived conservative stable region in solid black line is composed by the equal energy surface of Point M (0.6161, 0), and it is clearly seen that system trajectory is bounded inside the derived stable region, and locally asymptotically stability is ensured.Fig. 9(b) shows the change of energy function (dE/dt = −f (x)y 2 ) of the system after the disturbance.It is apparently that the change of system energy is always non-positive after the fault is cleared (t = 43 ms), which indicates that the system energy is constantly decreasing, and the system will finally be locally asymptotically stable within the obtained stable region.Specifically, the system energy change from Point A to Point C over one period (arc ABC) marked as the shadow area in Fig. 9(b) is calculated as −104.1 rad 2 /s 2 , and it is again verified that the system stability can be guaranteed within the stable region obtained from Lyapunov method.
System dynamic response under t c = 55 ms is shown in Fig. 10.Apparently, the system stability can be maintained even the system trajectory lies outside the derived stable region.This directly proves the conservative property of the obtained stable region from Lyapunov method.As shown in Fig. 10  the arcs AB and DE are outside the derived upper stability boundary and the damping coefficient of the system f (x) is negative, which leads to the change of energy of these two arcs positive.In other words, the system energy is increasing among these two arcs ( AB and DE).However, the arc BCD is within the obtained stable region, and the change of energy for BCD is kept as negative, which indicates the system energy is decreasing from Point B to Point D. Accordingly, the increased and decreased system energy over one period marked as green shadow area and blue shadow area in Fig. 10(b) are 32.6 rad 2 /s 2 and −289.8 rad 2 /s 2 , respectively.Therefore, the total energy change for system trajectory over one period is −257.2rad 2 /s 2 , and the system stability is ensured.

C. Unstable System Dynamic Behaviour
Comparing Figs. 9 and 10, it can be seen that with the increase of the fault duration time (t c ), the radius of the system trajectory is enlarged and the increased system energy over one period becomes gradually larger.Therefore, there exists one critical fault duration time, and the increased and decreased system energy over one period are mutually balanced.As a result, one "limit cycle" is finally formed in the phase plane as shown in Fig. 11(a) under t c = 83.1 ms (critical fault clearing time).It is apparent seen that, the limit cycle is naturally divided into five arcs, namely, arc AB, BC, CDE, EF , F G. Points B, C, E, F are the intersect points of the obtained stable region and limit Fig. 11.System dynamic behaviour on the limit cycle.cycle.Among five arcs, BC and EF are within the obtained stable region, which relates to the decreased system energy over one period calculated as −491.5 rad 2 /s 2 .However, arcs AB, CDE, and F G are outside the obtained stable region, and the damping coefficient f (x) is negative, which corresponds to the increased system energy over one period calculated as 490.6 rad 2 /s 2 .It is well verified from Fig. 11 that the change of energy on the limit cycle is nearly zero over one period.
As shown in Fig. 12, when the fault duration is enlarged to 88.5 ms, the system trajectory begins to deviate from the origin with the increasing radius of cycles in the phase plane.Finally, it will cross over the angle of UEP and presents as the aperiodic instability.As shown in Fig. 12(b), the change of system energy over one period can be divided into two parts, namely the increased system energy and the decreased system energy marked as the green and blue shadow areas.The increased system energy that involves the arcs AB, CDE, and F G is calculated as 779.1 rad 2 /s 2 .However, the decreased system energy that relates to the arcs of BC and EF is −526 rad 2 /s 2 .As a result, the total system energy change over one period from Point A to Point G is positive, indicating the occurrence of the system instability.It can be seen from Fig. 12 that once the system trajectory lies outside the limit cycle, or the fault duration time exceeds the critical fault clearing time, the system stability is no longer maintained.Moreover, it is concluded from Figs. 9 to 12 that with the increase of the fault duration time, there exists one unique limit cycle.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.II shows the values of H(x) at the angle of UEP under the different K PLL P .It is indicated that all values of H(x B, D ,F ) are larger than zero when K PLL P ≤ 0.028, and the derived existence conditions are satisfied, which implies that there exists one unique unstable limit cycle when K PLL P ≤ 0.028 as clearly shown in Fig. 13(b).Moreover, the farthest points of each limit cycle are located in the range of the upper stability boundary and the angle of UEP marked as the grey dashed line in Fig. 13(b).Fig. 14(a) indicates the value of defined function H(x) with respect to the variable x under the large value of K PLL P .It is interesting to find that in the studied range of the upper stability boundary of x A and the angle of UEP (x B ) under K PLL P = 0.03, value of H(x) are always negative, which violates the derived existence condition of limit cycle from (40).However, the system exists one unstable limit cycle as shown in Fig. 14(b).This verifies the conservation of the derived existence condition of limit cycle.It is worth mentioning that when K PLL P is increased to 0.034 and the fault duration time is 92.1 ms, the system trajectory as shown in Fig. 14(b) converges to SEP and the stability is maintained in this situation.However, the system trajectory crosses over the angle of UEP in the first swing and the classical aperiodic instability happens if t c = 92.2ms.This well explains the limit cycle does not exist when K PLL P is sufficient large, which also verifies that the studied dynamic system exists at most one limit cycle.If it exists, it is unstable.From Figs. 13 and 14, it is also interesting to find that with the decrease of K PLL P , the values of H(x) at the angle of UEP are gradually increasing, and the limit cycle is vulnerable to appear.
Certainly, when the controller parameters are set in the critical condition of the limit cycle appearance (K PLL P = 0.028, K PLL I = 20), the disturbances from outside uncertainties might lead to the inaccuracy of the analysis.Normally, to ensure the reliability of the analysis results or the robustness of the proposed condition, the partial conservativeness of the condition will be sacrificed.In this case, 30% variation of controller parameters can still guarantee the existence of limit cycle under H(x UEP ) ≥ 0.25.
The defined function H(x) with respect to the variable x under the small value of K PLL I is shown in Fig. 15(a).In the studied range of the upper stability boundary (x A, C, E ) and the angle of UEP (x B, D, F ), the maximum values of H(x) are obtained at the angle of UEP, and they are all positive as indicated by Table III.As a result, the studied dynamic system exists one unique unstable limit cycle under K PLL I ≥ 20 as clearly shown in Fig. 15(b).Moreover, with the increase of K PLL I , the limit cycles are getting smaller, which well indicates the stability of the studied dynamic system is deteriorated with the increase of K PLL I .Fig. 16 indicates the value of defined function H(x) with respect to the variable x under the small value of K PLL I .It can be  seen from Table III and Fig. 16(a) that the value of the defined function H(x) at the angle of UEP is −0.0602 under K PLL I = 19, and it violates the derived existence condition of limit cycle from (40).However, as indicated from Fig. 16(b), the system exists one unstable limit cycle marked in black solid line.It is again verified the conservation is guaranteed for the derived existence condition of limit cycle.It can be also found from Fig. 16(b) that when K PLL I is further decreased to 16, the system does not exist any limit cycle.This can be explained by that when the fault duration time is 104.5 ms, stability is maintained since the system trajectory converges to SEP.However, the system trajectory directly crosses over the angle of UEP without any cycling and the classical aperiodic instability happens under t c = 104.6 ms.It is also interesting to find from Table III that, the values of H(x) at the angle of UEP are gradually increasing with the increase of K PLL I , and the limit cycle is more likely to appear.

E. Application to One multi-VSCs System
One simple two VSCs system as shown in Fig. 17 is utilized to verify the effectiveness of the proposed method.VSC1 and VSC2 delivers 140 MW and 120 MW active power to the infinite bus via the transmission lines, respectively.The resistances and inductances of l1, l2, and l3 are 0.02 pu and 0.2 pu, 0.025 pu and 0.25 pu, 0.04 pu and 0.4 pu, respectively.The disturbance is that the infinite source U sL suddenly drops to 0.3 pu at time t = 0.01 second, and returns to U sL after t c seconds.
As shown in Fig. 18(a), VSC1 presents the oscillation instability and VSC2 keeps stable with the damped oscillation when the fault duration time (t c ) is 51 ms.However, the stability for both VSCs cannot be maintained any more when t c is enlarged to 63 ms.As indicated by Fig. 18(b) and (c) that the phenomenon of limit cycles can still appear in the multi-VSC embedded power system.Moreover, the stability of VSC2 is better than that of VSC1 since the limit cycle of VSC1 appears when t c is 43 ms, and the limit cycle of VSC2 occurs when t c reaches 56 ms.This is well explained by the fact that VSC1 delivers more active power than VSC2, which indicates that the heavy load condition will endanger the system stability.
To apply the proposed existence condition into the multi-VSCs embedded power system, it is necessary to decouple the original complex network into several independent subsystems that one single VSC connected to the infinite bus system via one equivalent impedance.The parameter of the equivalent impedance can be derived by Thevenin's theorem.As shown in Fig. 18(d), the values of H(x) are positive at UEPs for both VSCs.As a result, it predicts the existence of limit cycles.It should be noted that the network decoupling errors are normally hard to be quantified so that the conservativeness of this proposed method cannot be ensured.In spite of this, it can be regarded as one alternative heuristic method to predict limit cycles in the multi-VSCs embedded system.from Fig. 19(a) that with the increase of the transmitted active power, the stable region of the grid-connected VSC becomes much smaller.This can be well explained by the fact that the heavy load condition will severely deteriorate the system stability.As a result, the limit cycle of the studied dynamic system is more likely to appear.Fig 19(b) investigates the impacts of reactive power of VSC on the system stability.It can be clearly seen that the stable region mapped by limit cycle is enlarged when VSC generates amount of reactive power.This is reasonable since the PCC voltage of VSC can be somehow leveled up so that the power transfer capability of VSC is largely improved.

G. Impacts of Controller Time Delays on Limit Cycle
Time delay is inevitable in the digital control system.It will change the phase characteristics of the system and affect the system stability.The influence of time delay on the stability of the VSC system has been discussed in [23], [24].This paper mainly investigates the impacts of controller time delays on the limit cycle of the grid-connected VSC.Time delay of the tracking angle provided by PLL (output of PLL) is set to simulate the delay of the digital controller.The value of the controller time delay noted as T d is normally 1−2 control cycle T of the digital controller [25].Accordingly, the controller time delay is around 0.05−0.2ms.Fig. 20 shows the limit cycles of the grid-connected VSC under the different controller delays.It can be clearly seen that with the increase of the controller time delay, the stable region of the grid-connected VSC becomes smaller.It can be inferred that the time delay from the digital controller might deteriorate the system stability to some extent.Moreover, limit cycles still appear in the appearance of the controller time delay as indicated from Fig. 20.The results show that the controller time delay only has an insignificant effect on the parameters of the limit cycle, but does not affect the basic properties (existence and uniqueness) of the limit cycle.

VI. CONCLUSION
This paper investigates the large disturbance stability of one VSC connected to the infinite bus through the transmission line.The phenomenon of limit cycle is discovered in this system, and the conditions of existence and uniqueness of limit cycles are derived.Based on Poincaré Theorem, it is found that there exists a unique unstable limit cycle between the Lyapunov stable boundary and the UEP, and the limit cycle appears when the negative damping effects counteract with the positive damping effects over one period.The research results of the impact of controller parameters, operating conditions, and the controller time delays on the existence of limit cycles show that the proposed existence condition can predict limit cycles with certain degree of conservativeness.The controller time delays and other uncertainties do not affect the basic properties (existence and uniqueness) of the limit cycle, although there are some effects on the parameters of the limit cycle.These findings in the paper might be well-beneficial for the understanding of further power systems with high penetrations of CIGs, and provide useful guidance for the controller design and parameters optimization of CIGs.APPENDIX Results in Figs.5-7 are obtained based on the system model of (7), and the related parameters utilized are shown in Table I.The coordinates of Points on limit cycle in Fig. 5 are A(0.62,85.33), B(0.62, −90.69),C(−1.17, −57.89), and D(−1.17, 58.53).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Manuscript received 18
June 2022; revised 16 September 2022, 5 November 2022, and 15 December 2022; accepted 15 January 2023.Date of publication 19 January 2023; date of current version 26 December 2023.This work was supported by the National Natural Science Foundation of China under Grant U2066602.Paper no. TPWRS-00775-2022.(Corresponding author: Zhengchun Du.) , u abc c and u abc p are the voltage of the converter and the point of common coupling.R c and L c are the equivalent resistance and inductance of the phasor reactor.R l and L l are the line resistance and inductance, respectively.i abc c is the current flowing through the converter, and u abc s is the voltage of the infinite source with the fixed angular speed of ω s .

Fig. 1 .
Fig. 1.Scheme of one VSC connected to the infinite bus via line.

Fig. 2 .
Fig. 2. Relationship of variables in different frames in vector space.

Fig. 9 .
Fig. 9. System dynamic behaviour under the fault duration time of 12 ms.
(a), system trajectory (arc AE) within a specified period T can be divided into three separated arcs, marked as arc AB, BCD, DE.Point B and Point D are the intersect points of the obtained upper stability boundary and system trajectory.It can be inferred that Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 10 .
Fig. 10.System dynamic behaviour under the fault duration time of 20 ms.

Fig. 13 (
Fig.13(a)shows the values of the defined function H(x) with respect to the variable x under the small value of K PLL P .It can be clearly seen from Fig.13(a) that in the range of the upper stability boundary (x A, C, E ) and the angle of UEP (x B, D, F ), function H(x) is first decreasing and then increasing with the increase of x.Moreover, the maximum values of H(x) under the different K PLL P are obtained at the angle of UEP as indicated by Fig.13(a).TableIIshows the values of H(x) at the angle of UEP under the different K PLL P .It is indicated that all values of H(x B, D ,F ) are larger than zero when K PLL

Fig. 13 .
Fig.13.Limit cycles of the grid-connected VSC under small value of K PLL P .

Fig. 14 .
Fig. 14.Limit cycles of the grid-connected VSC under large value of K PLL P .

Fig. 15 .Fig. 16 .
Fig. 15.Limit cycles of the grid-connected VSC under large value of K PLL I .

Fig. 19
Fig.19shows the impacts of different operating conditions on the stability of the grid-connected VSC.It can be clearly seen

Fig. 19 .
Fig. 19.Dynamic response of PLL under different operating conditions.

TABLE III VALUES
OF DEFINED FUNCTION AT THE ANGLE OF UEP UNDER THE DIFFERENT K PLL I