Transient Stability Analysis Method for Multi-Inverters Based on the Two-Level Iterative Equal Area Rule

Abstract

This invention relates to power electronics technology, specifically to a method for transient stability analysis of multi-inverter systems based on a two-layer iterative equal-area method. The method involves establishing a mathematical model of a multi-grid-connected inverter parallel system, defining the equivalent mechanical power, electromagnetic power, other-machine electromagnetic power, electromagnetic interaction power, and self-damping and mutual damping of the system. First, non-integrable terms are scaled down to make them integrable. Then, using the proposed two-layer iterative equal-area method, the upper and lower boundaries of the power angle stability of the multi-grid-connected inverter parallel system are obtained through two iterations. This method, by combining iteration and scaling, improves the conservatism of the stability region estimate as much as possible without causing stability misjudgments. Using the two-layer iterative equal-area method, a relatively accurate estimate of the stability region of the grid-connected inverter parallel system is obtained, improving the conservatism of the original equal-area method applied to VSC parallel systems. It has good development potential and room for promotion.

Description

Technical Field

[0001] This invention belongs to the field of power electronics technology, and specifically relates to a method for transient stability analysis of multi-inverter based on the two-layer iterative equal area rule. Background Technology

[0002] As an effective and reliable distributed integration method, converter-dominated power systems have received increasing attention in recent years. Currently, widely used grid-connected converters are mainly divided into two categories: grid-forming converters (GFM-VSC) and grid-following converters (GFL-VSC). Compared to GFM-VSC, GFL-VSC, due to its low inertia, is more prone to instability during severe grid faults. Instability leads to converter disconnection from the grid, seriously threatening the safe and stable operation of the system. Therefore, the stability analysis of GFL-VSCs has attracted great interest from academic and industrial peers. Phase-locked loops (PLLs) play a dominant role in the transient dynamic response of GFL-VSCs. Therefore, the stability analysis of GFL-VSC mainly focuses on PLL modeling and dynamic processes. Current analyses of GFL-VSC stability are mainly based on small-signal models, which cannot analyze the system's transient stability. Various large-signal stability assessment methods have been proposed, but they are often limited to single-machine systems and cannot be extended to multi-machine systems. Therefore, large-signal stability analysis methods applicable to complex VSC systems still require further research. Summary of the Invention

[0003] To address the problems existing in the background technology, this invention provides a method for transient stability analysis of multiple inverters based on the two-layer iterative equal area rule.

[0004] To solve the above-mentioned technical problems, the present invention adopts the following technical solution: a transient stability analysis method for multi-inverter based on the two-layer iterative equal-area rule, comprising the following steps:

[0005] Step 1: Establish a fourth-order mathematical model of the grid-connected inverter parallel system considering line frequency fluctuations, including mechanical torque, electromagnetic torque, electromagnetic interaction torque, self-damped torque, and mutual damped torque.

[0006] Step 2: Scaling the non-integrable electromagnetic interaction terms, other machine power terms, self-damping and mutual damping into integrable terms, and dividing the positive and negative distribution intervals of damping.

[0007] Step 3: Based on the scaling in Step 2, the double-iteration equal-area method is obtained. The double-iteration equal-area method includes iteration A and iteration B, where iteration A handles electromagnetic interaction power, other machine electromagnetic power, and mutual damping terms; iteration B handles self-damping terms; iteration A obtains the estimated upper boundary of the power angle stability, and iteration B obtains the estimated lower boundary of the power angle stability based on iteration A.

[0008] In the above-mentioned transient stability analysis method for multi-inverter based on the two-layer iterative equal-area rule, the implementation of step 1 includes the following steps:

[0009] Step 1.1: Establish a nonlinear model of the grid-connected inverter parallel system;

[0010] The parallel grid-connected inverter VSC system includes: VSC1 with filter impedance L f1 Then, through the line impedance L1, R1 and the grid impedance L g ,R g With the power grid V g Connected; VSC2 at filter impedance L f2 Then, through the line impedance L2, R2 and the grid impedance L g ,R g With the power grid V g Connected; VSC1 and VSC2 both use direct current control and phase-locked loop (PLL);

[0011] According to Kirchhoff's voltage law, the voltage V at the grid connection point of VSC1 is... t1 The voltage V at the grid connection point of VSC2 t2 They are represented as follows:

[0012]

[0013] in The power factor of VSC1, These are the dq-axis components of the line current I1, respectively. The power factor of VSC2, These are the dq-axis components of the line current I2, respectively; θ PLL1 θ is the reference phase of the VSC1 phase-locked loop output. PLL2 It is the reference phase of the phase-locked loop output of VSC1; Applying the PARK transformation to both sides of equation (1) yields V t1 q-axis component V t1q , and V t2 q-axis component V t1q quantity:

[0014]

[0015] The dynamic equation of the phase-locked loop is:

[0016]

[0017] Define δ1 = θ PLL1 -θ g Let θ be the angle of motion of VSC1; define δ2 = θ PLL2 -θ gTo obtain the dynamic equations of the grid-connected inverter parallel system, we can combine equations (2) and (3) to get the power angle of VSC2.

[0018]

[0019] Where M1 represents the equivalent inertia coefficient of VSC1, a0 represents the equivalent mechanical power of VSC1, a1 represents the electromagnetic power between VSC1 and the power grid, a2 represents the electromagnetic power between VSC2 and the power grid, and a3 represents the interactive electromagnetic power between VSC1 and VSC2. This indicates the self-damping of VSC1; b1 represents the mutual damping of VSC1; M2 represents the equivalent inertia coefficient of VSC2; b0 represents the equivalent mechanical power of VSC2; b1 represents the electromagnetic power between VSC2 and the power grid; b2 represents the electromagnetic power between VSC1 and the power grid; b3 represents the interactive electromagnetic power between VSC2 and VSC1. This indicates the self-damping of VSC2; Indicates the mutual damping of VSC2;

[0020] The expressions for each coefficient are as follows:

[0021]

[0022] a1=K i V g

[0023]

[0024]

[0025]

[0026]

[0027]

[0028] a4=-K i I1(L g +L2)

[0029] a5 = K p V g

[0030]

[0031]

[0032]

[0033]

[0034]

[0035]

[0036]

[0037]

[0038]

[0039]

[0040]

[0041] b1=K i In g

[0042]

[0043]

[0044]

[0045]

[0046]

[0047] b4=-K i I2(L g +L1)

[0048] b5=K p In g

[0049]

[0050]

[0051]

[0052]

[0053]

[0054]

[0055]

[0056]

[0057]

[0058]

[0059] Step 1.2: Compare the VSC parallel system and the SG system to obtain the equivalent mechanical power, electromagnetic power, electromagnetic interaction power, self-damping, and mutual damping.

[0060] From the VSC1 power angle equation (4), we know that the mechanical power of VSC1 is a0, which is a constant term; the electromagnetic power of VSC1 to the power grid is a1sinδ1, which is only related to the sine of its own power angle δ1; the electromagnetic power of the other machine VSC2 is a2sinδ2, which is only related to the sine of the power angle δ2 of the other machine VSC2; the electromagnetic power between VSC1 and VSC2 Only related to the power angle δ between VSC1 and VSC2 12 It is related to the sine. For constant terms; Defined as the equivalent self-damping coefficient D of VSC1 self1 This includes the damping a4cos(θ1-θ1) generated by the output of VSC1; and the damping a5cos(θ1-θ1) generated by the grid output. g Damping generated by the output of VSC2 The equivalent self-damping coefficient D of VSC1 self1 It is a damping term related to the angular frequency ω1 of the local VSC1, generated by three different frequency or phase classifications. The physical meaning of each cos term is to transform the output of other phases into the coordinate system corresponding to the local phase angle δ1 through PARK transformation, and is defined as self-damping. The equivalent mutual damping coefficient D, defined as VSC1 mutu1 This includes the output-related terms of VSC2, a7cos(θ2-θ2); and the output-related terms of the power grid, a8cos(θ2-θ2). g ); VSC1 output related terms a9cos(θ2-θ1); VSC1 equivalent mutual damping coefficient D mutu1 It is related to the angular frequency ω2 of the VSC2 in the other machine;

[0061] From the VSC2 power angle equation (4), we know that the mechanical power of VSC2 is b0, which is a constant term; the electromagnetic power of VSC2 to the power grid is b1sinδ2, which is only related to the sine of its own power angle δ2; the electromagnetic power of the other machine VSC1 is b2sinδ1, which is only related to the sine of the power angle δ1 of the other machine VSC1; the electromagnetic power between VSC2 and VSC1 is... Only the power angle δ between VSC2 and VSC1 21 It is related to the sine. For constant terms; The equivalent self-damping coefficient D, defined as VSC2 self2This includes the damping b4cos(θ2-θ2) generated by the output of VSC2; and the damping b5cos(θ2-θ2) generated by the grid output. g Damping generated by the output of VSC1 The equivalent self-damping coefficient D of VSC2 self2 It is a damping term related to the angular frequency ω2 of the local VSC2, generated by three different frequency or phase classifications. The physical meaning of each cos term is to transform the output of other phases into the coordinate system corresponding to the local phase angle δ2 through PARK transformation, and is defined as self-damping. The equivalent mutual damping coefficient D, defined as VSC2 mutu2 This includes output-related terms for VSC1: b7cos(θ1-θ1); and output-related terms for the power grid: b8cos(θ1-θ1). g ); VSC2 output related terms b9cos(θ1-θ2); VSC2 equivalent self-damping coefficient D mutu2 The frequency ω1 of the other machine angle VSC1 is multiplied to produce damping.

[0062] In the above-mentioned transient stability analysis method for multi-inverter based on the two-layer iterative equal-area rule, step 2 includes the following steps:

[0063] Step 2.1, Damping Positive and Negative Analysis; From Equation (4), we know that:

[0064] a) Self-damping related to the machine's angular velocity and

[0065] b) Mutual damping related to the angular velocity of other converters The damping term varies with the system's state of motion; D self1 D self2 D mutu1 and D mutu2 Its positive or negative value can change;

[0066] During the transient process under normal conditions, the mutual damping is always negative; at the same time, the self-damping is positive near the equilibrium point and will move towards the negative damping region during the transient process; the direction of the negative damping in the dynamic process is always consistent with the direction of the system motion, which plays an accelerating role. When analyzing the transient stability of the grid-connected inverter parallel system, the influence of damping on the system should be ignored.

[0067] Step 2.2: Scaling up the electromagnetic interaction term and the self-damping mutual damping;

[0068] First, scale the corresponding non-integrable terms; scale the interaction terms containing the rotor angle of other machines, and convert the non-integrable terms into integrable terms that only depend on the rotor angle of the local machine; use the method based on the monotonicity of trigonometric functions to scale the interaction terms; use an iterative algorithm, and use the stable boundary calculated in the previous step as the scaling basis for the next calculation; for mutual damping, use the iterative improved equal-area criterion, and use the maximum angular velocity to scale the mutual damping term, and include it in the iterative process at the same time;

[0069] Consider the interaction electromagnetic power term and the power term of other machines as a whole, and define it as P 1inter ; At the same time, the steady-state value of P 1inter is denoted as P 1inter_e , as shown in (5)-(6):

[0070]

[0071]

[0072] where δ ie is the steady-state equilibrium point of δ i in (4); from the coefficient expressions in step 1.1, it can be seen that 0 < a3 < a2 ≈ a0 < a1. For P 1inter , there exists a lower bound P 1inter_min :

[0073] P 1inter ≥ a2sin(δ2) - a3 ≥ a2sin(δ 2mine ) - a3 = P 1inter_min (7)

[0074] where δ 2min_e is the stable left boundary of VSC2 calculated as in (8)-(9), and its corresponding maximum accelerating area is obtained by considering b0 - P 2inter_e as the equivalent mechanical power and ignoring the damping effect; in the subsequent iterative process, δ 2min_e will be continuously updated to obtain a more accurate stable region estimate:

[0075] δ 2max_e = π - arcsin((b0 - P 2intere ) / b1) (8)

[0076]

[0077] Since a0 - P 1inter_min0 ≥ a0 - P 1inter , therefore, when using the equal-area criterion to analyze the transient stability of the grid-connected inverter parallel system, using P 1inter_min to replace P 1inter will not lead to misjudgment of stability; take a0 - P1inter_min Assuming mechanical power and neglecting damping effects, the corresponding stability boundary [δ] is derived using the equal area method. 1min_int ,δ 1max_int Estimate:

[0078] δ 1max_int =π-arcsin((a0-P) 1inter_min ) / a1) (10)

[0079]

[0080] Similarly, for VSC2, there exists an interactive power lower bound P. 2inter_min :

[0081] P 2inter ≥b2sin(δ1)+b3≥b2sin(δ 1min_e )-b3=P 2inter_min (12)

[0082] Similar to equations (10)-(11), b0-P 2inter_min To estimate the equivalent mechanical power of VSC2, a stability boundary [δ] of VSC2 is estimated using the equal area method. 2min_int ,δ 2max_int ];

[0083] Due to the adverse effects of interactive power on the transient stability of the system, let b0-P have time-varying characteristics. 2inter For mechanical power, the actual stable left boundary δ 2min will be located at δ 2min_e To the right of; let a0-P 1inter_min The mechanical power will cause the deceleration region of VSC1 to become smaller, and δ obtained from equations (7) and (10)-(11) 1min_int Relatively conservative;

[0084] The effects of self-damping and mutual-damping terms on a rotating system are either the consumption of the system's energy in the form of work, or, when negative, the input capacity; satisfying the form ∫D self1 (δ1,δ2)ω1dδ1 and ∫D mutu1 (δ1,δ2)ω2dδ1; Considering negative damping as part of the mechanical power, which varies with the system's state of motion; the damping term itself is non-integrable, scaling the energy exerted on the system by negative damping to its maximum; according to the coefficient expression in step 1.1: a7,b7,a9,b9<0 and a8,b8>0, we have:

[0085]

[0086] Where δ imax It is δ iThe negative work done by the mutual damping at the stable upper boundary is scaled as shown in (14):

[0087]

[0088] An improved equal-area rule based on scaling by the maximum angular frequency is adopted. IEAC scales the negative cross-damping to equation (15) by taking the maximum value of the angular frequencies of other VSCs:

[0089]

[0090] Where ω imax_e The maximum angular frequency corresponding to VSC1 or VSC2 in the dynamic process described in equations (8)-(9); the actual stability boundary is [δ imin_e ,δ imax_e It is a subset of ], therefore the actual maximum angular frequency is less than ω. imax_e This conservatism can be minimized through subsequent iterations A.

[0091]

[0092] For self-damping D selfi Since there exist a4,b4<0 and a5,b5,a6,b6>0, then D selfi Shrink to D selfimin (δ i ):

[0093]

[0094] D si_min (δ i ) is δ i It is a single-variable function, therefore the self-damped D self1 The negative work done can be scaled as shown in equation (18):

[0095]

[0096] Combining equations (7), (12), (15), and (17), the transient stability boundary of the grid-connected inverter parallel system can be calculated using the equal area method.

[0097] In the above-mentioned transient stability analysis method for multi-inverter based on the two-layer iterative equal-area rule, step 3 includes the following steps:

[0098] Step 3.1: The double-iteration equal-area method is used to calculate the transient stability boundary of the grid-connected inverter parallel system, including iteration A and iteration B;

[0099] Iteration A, considering the effects of interactive electromagnetic power and mutual damping terms, yields the upper stable boundary δ. 1max and δ 2maxThe estimate; Iteration B is based on the δ obtained from iteration A. 1max and δ 2max The lower boundary δ is calculated using the extended iterative equal area method. 1min and δ 2min ;

[0100] Step 3.2, Iteration A includes:

[0101] A stability boundary estimate [δ] is obtained during the j-th iteration. 1minA_j ,δ 1maxA_j ] and [δ 2minA_j ,δ 1maxA_j This will be used as the basis for scaling in the next iteration;

[0102] Define P 1A_j As the electromagnetic interaction power of VSC1 after scaling in the j-th iteration A, it is the sum of its mechanical power and mutual damping; define P 2A_j As the electromagnetic interaction power of VSC2 after scaling in the j-th iteration A, it is the sum of its mechanical power and mutual damping:

[0103]

[0104] Where ω 1max_j-1 ,δ 1minA_j-1 andδ 1imaxA_j-1 These are the maximum angular frequency of VSC1, the stable lower boundary, and the stable right boundary, respectively, in the (j-1)th iteration; ω 2max_j-1 ,δ 2minA_j-1 andδ 2imaxA_j-1 These represent the maximum angular frequency of VSC2, the stable lower boundary, and the stable right boundary, respectively, in the (j-1)th iteration. The calculation process in the j-th iteration is as follows:

[0105]

[0106]

[0107]

[0108] During the first iteration, P 1A_1 and P 2A_1 The calculation process is as follows:

[0109]

[0110] From the above derivation, it can be seen that a new scaling reference is obtained from (20)-(22) and is used cyclically in scaling (19); if δ 1maxA_j and δ 2maxA_j If the iteration converges to a given precision ε, then the iteration ends, and the stable boundary on δ1 can be taken as δ. 1maxAThe stable boundary on δ2 can be taken as δ 2maxA Otherwise, assume j = j + 1, and begin the (j+1)th iteration A; the final output convergence value δ 1maxA and δ 2maxA The proposed dual-iteration EAC method serves as the final upper bound δ. 1max and δ 2max ;

[0111] Step 3.3, Iteration B includes:

[0112] Define the initial self-damped power and initial angular frequency ω of the system. 1B_0 (δ1) and ω 2B_0 (δ2), where δ1 and δ2 in parentheses indicate that the required angular frequency ω is... 1B_0 and ω 2B_0 Distribution functions relative to δ1 and δ2 respectively:

[0113]

[0114]

[0115] angular velocity ω 1B_0 (δ1) and ω 2B_0 (δ2) is calculated using the variable lower limit integral as shown in (25); considering D s1_min Independent of ω1, D s2_min It is independent of ω2, therefore it does not need to be updated in each iteration;

[0116] In the j-th iteration, the self-damping power P of VSC1 1B_j Based on the angular velocity distribution ω at the (j-1)th iteration 1B_j-1 (δ1) and minimum work angle δ 1minB_j-1 The self-damping power P of VSC2 is calculated. 2B_j Based on the angular velocity distribution ω at the (j-1)th iteration 2B_j-1 (δ2) and minimum work angle δ 2minB_j-1 The calculation shows that:

[0117]

[0118] According to the principle of the equal area method, the angular velocity distribution ω of VSC1 obtained in the j-th iteration is... 1B_j Angular velocity distributions ω of (δ1) and VSC2 2B_j (δ2) are shown below:

[0119]

[0120] Meanwhile, the lower boundary δ of the VSC1 power angle stability domain at the j-th iteration is... 1minB_jThe lower boundary δ of the VSC2 power angle stability region 2minB_j It can also be calculated as:

[0121]

[0122] Equations (24)-(25) give the initial value P for the iteration. iB_0 and ω iB_0 (δ i The calculation formula is as follows: In the j-th cycle, the angular velocity distribution ω of VSC1 is calculated using (27)-(28). 1B_j (δ1) and the lower boundary of the work angle stability δ 1minB_j and the angular velocity distribution ω of VSC2 2B_j (δ2) and the lower boundary of the work angle stability δ 2minB_j The calculation results of equations (27)-(28) are used as inputs for the (j+1)th iteration and substituted into (26) to calculate the self-damping power P at the (j+1)th iteration. 1B_j and P 2B_j If the δ obtained by equation (28) 1minB_j and δ 2minB_j If the iteration converges to a given precision ε2, then the iteration exits. The lower boundaries of δ1 and δ2 are written as δ1 and δ2, respectively. 1min and δ 2min The iterative process of equations (24)-(38) is the second iteration B of the extended iterative equal area rule.

[0123] The power angle stability boundary of the grid-connected inverter parallel system is obtained through double iteration: [δ 1min ,δ 1max ] and [δ 2min ,δ 2max Iteration A yields a stability boundary estimate considering the effects of interactive electromagnetic power, other electromagnetic power, and mutual damping terms, denoted as [δ]. 1minA ,δ 1maxA ] and [δ 2minA ,δ 2maxA ]; where δ 1maxA and δ 2maxA As the upper bound of the double iteration δ 1max and δ 2max The result of iteration A is used as the basis for the initial value calculation of iteration B. Based on iteration A, considering the influence of the self-damping term, iteration B obtains the lower boundary δ using the extended iterative equal area method. 1minB and δ 1minB δ serves as the lower bound of the double iteration. 1min and δ 2min Combining iterations A and B, the final stable boundary [δ] is obtained. 1min ,δ 1maxA ] and [δ 2min ,δ2max ].

[0124] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0125] An improved equal-area rule is used to analyze the transient stability of a grid-connected inverter based on a phase-locked loop (PLL). Based on the modeling and analysis of the PLL, the rotor equations of a synchronous motor are compared, and the terms in the inverter's state equations are equivalent to the mechanical power, electromagnetic power, and damping in the synchronous motor. Furthermore, the equal-area rule is applied to analyze the large-signal stability of the inverter. The reason why the existing equal-area rule is unsuitable for VSCs due to the variable damping term is pointed out, and a corresponding improvement is made using a maximum angular velocity ω obtained by linearly approximating the system. max2 By substituting the damping term and scaling up the negative damping, a relatively accurate estimate of the stability region is obtained using the improved equal-area rule. This significantly improves the conservatism of the original equal-area rule when applied to VSC systems, and has good development potential and room for promotion. Attached Figure Description

[0126] Figure 1 This is a flowchart of the transient stability analysis method for grid-connected inverter parallel systems based on the double-iteration equal area rule in an embodiment of the present invention;

[0127] Figure 2 This is a schematic diagram of the structure and control of a grid-connected inverter parallel system according to an embodiment of the present invention;

[0128] Figure 3 The dynamic characteristics of the mathematical model and simulation model under large disturbances in the embodiments of the present invention;

[0129] Figure 4 This is a schematic diagram showing the distribution of positive and negative damping regions in an embodiment of the present invention;

[0130] Figure 5 This is a flowchart of the interactive electromagnetic power and mutual damping term EAC iteration based on iteration A in an embodiment of the present invention;

[0131] Figure 6 This is a flowchart illustrating the iterative process of the self-damping term based on iteration B in an embodiment of the present invention.

[0132] Figure 7 This is a schematic diagram illustrating the large-signal synchronization stability analysis of a parallel converter system based on the dual-iteration EAC method according to an embodiment of the present invention. Detailed Implementation

[0133] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0134] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.

[0135] The present invention will be further described below with reference to specific embodiments, but these are not intended to limit the scope of the invention.

[0136] This embodiment improves upon the traditional equal-area rule to adapt it to the specific characteristics of parallel systems with multiple grid-connected inverters. An iterative method is employed to improve the conservatism of existing methods, while also establishing a more accurate model that considers the impact of line frequency fluctuations on line impedance. Most existing methods neglect the influence of damping on system stability. This embodiment proposes that damping in parallel systems with multiple grid-connected inverters varies with the system state. Self-damping exhibits positive damping characteristics near the equilibrium point, but enters the negative damping region during transient processes. Mutual damping, on the other hand, is mostly in the negative damping region, regardless of whether it is transient or steady-state. Due to the adverse effects of negative damping on system stability, its influence cannot be ignored when quantitatively calculating transient stability. By combining iteration and scaling, the conservatism of the obtained stability region estimate is improved as much as possible without causing stability misjudgments. Scaling is used to transform non-integrable terms into integrable terms, but this introduces a certain degree of conservatism, which can be mitigated as much as possible through iterative methods. This embodiment uses the double-iteration equal-area rule to obtain a relatively accurate estimate of the stability domain of a grid-connected inverter parallel system, which greatly improves the conservatism of the original equal-area rule when applied to VSC parallel systems.

[0137] This embodiment is achieved through the following technical solution: a transient stability analysis method for multiple inverters based on the two-layer iterative equal-area rule, such as... Figure 1 As shown, the process includes the following steps: First, a mathematical model of a multi-grid-connected inverter parallel system is established. The equivalent mechanical power, electromagnetic power, other machine electromagnetic power, electromagnetic interaction power, and self-damping and mutual damping of the multi-grid-connected inverter parallel system are defined. Firstly, non-integrable terms are scaled to make them integrable. Then, using the proposed double-iteration equal-area method, the upper and lower boundaries of the power angle stability of the multi-grid-connected inverter parallel system are obtained through two iterations.

[0138] S1. Establish a fourth-order mathematical model for a grid-connected inverter parallel system considering line frequency fluctuations. Its equations are similar to the rotor motion equations of a synchronous motor. This includes mechanical torque, electromagnetic torque (related to the inverter's power angle), electromagnetic interaction torque (related to the power angle difference between the inverter and other inverters), self-damping (related to the inverter's frequency), and mutual damping (related to the frequency of other inverters). The generation mechanism and physical meaning of each term are also given.

[0139] S1.1: Nonlinear modeling of a multi-grid-connected inverter parallel system yields a fourth-order nonlinear differential equation system.

[0140] like Figure 2 The parallel grid-connected inverter VSC system shown has VSC1 at filter impedance L f1 Then, through the line impedance L1, R1 and the grid impedance L g ,R g With the power grid V g Connected. VSC2 is connected at filter impedance L. f2 Then, through the line impedance L2, R2 and the grid impedance L g ,R g With the power grid V g Connected. Both VSCs employ direct current control and a phase-locked loop (PLL). The basic structure of the PLL is as follows: Figure 2 As shown in the dashed box, the input of the phase-locked loop (PLL1) of VSC1 is the grid connection point voltage V. t1 After PARK transformation, the q-axis component V is obtained. t1q After passing through a PI controller, where k p k I These represent the proportional gain and integral gain of the PLL, respectively, plus the nominal frequency ω. n The output angular frequency ω of PLL1 can then be obtained. PLL1 After integration, the output phase θ of PLL1 is obtained. PLL1 The input of the phase-locked loop (PLL2) of VSC2 is the grid connection point voltage V. t2 After PARK transformation, the q-axis component V is obtained. t2q After passing through a PI controller, where k p k I These represent the proportional gain and integral gain of the PLL, respectively, plus the nominal frequency ω. n The output angular frequency ω of PLL2 can then be obtained. PLL2 After integration, the output phase θ of PLL1 is obtained. PLL2 θ g and ω g These represent the angle and angular frequency of the grid, respectively. Assume ω... g =ω nWe ignore fluctuations in the power grid frequency.

[0141] According to the multi-timescale decoupling theory, the fast dynamics of the current loop and line current can be neglected in the phase-locked loop time scale. Therefore, it is assumed that the output line current I1 of VSC1 is equal to the current loop input reference value I of VSC1. ref1 The output line current I2 of VSC2 is equal to the current loop input reference value I of VSC2. ref2 Furthermore, the line current can be expressed using algebraic equations. Based on Kirchhoff's voltage law, the voltage V at the grid connection point of VSC1 is... t1 The voltage V at the grid connection point of VSC2 t2 They can be represented as:

[0142]

[0143] in The power factor of VSC1, These are the dq-axis components of the line current I1, respectively. The power factor of VSC2, These are the dq-axis components of the line current I2, respectively. θ PLL1 θ is the reference phase of the VSC1 phase-locked loop output. PLL2 This is the reference phase of the VSC1 phase-locked loop output. Normally, it is considered... Applying the PARK transformation to both sides of equation (1) yields V t1 q-axis component V t1q , and V t2 q-axis component V t1q quantity:

[0144]

[0145] Depend on Figure 2 It can be seen that the dynamics of the phase-locked loop are:

[0146]

[0147] Define δ1 = θ PLL1 -θ g Let δ2 be the angle of attack of VSC1. Define δ2 = θ PLL2 -θ g To obtain the dynamic equations of the grid-connected inverter parallel system, we can combine equations (2) and (3) to get the power angle of VSC2.

[0148]

[0149] Where M1 represents the equivalent inertia coefficient of VSC1, a0 represents the equivalent mechanical power of VSC1, a1 represents the electromagnetic power between VSC1 and the power grid, a2 represents the electromagnetic power between VSC2 and the power grid, and a3 represents the interactive electromagnetic power between VSC1 and VSC2. This indicates the self-damping of VSC1; b1 represents the mutual damping of VSC1; b2 represents the equivalent inertia coefficient of VSC2; b0 represents the equivalent mechanical power of VSC2; b1 represents the electromagnetic power between VSC2 and the power grid; b2 represents the electromagnetic power between VSC1 and the power grid; b3 represents the interactive electromagnetic power between VSC2 and VSC1. This indicates the self-damping of VSC2; The mutual damping of VSC2 is represented in (4); the specific derivation and expressions of each coefficient are as follows:

[0150]

[0151] a1=K i V g

[0152]

[0153]

[0154]

[0155]

[0156]

[0157] a4=-K i I1(L g +L2)

[0158] a5 = K p V g

[0159]

[0160]

[0161]

[0162]

[0163]

[0164]

[0165]

[0166]

[0167]

[0168]

[0169]

[0170] b1 = K i V g

[0171]

[0172]

[0173]

[0174]

[0175]

[0176] b4 = -K i I2(L g +L1)

[0177] b5 = K p V g

[0178]

[0179]

[0180]

[0181]

[0182]

[0183]

[0184]

[0185]

[0186]

[0187]

[0188] S1.2: Compare the VSC parallel system and the SG system to obtain the equivalent mechanical power, electromagnetic power, electromagnetic interaction power, self-damping and mutual damping.

[0189] The physical meanings of the coefficients in the VSC1 power angle equation (4) are as follows: the mechanical power of VSC1 is a0, which is a constant term; the electromagnetic power of VSC1 to the power grid is a1sinδ1, which is only related to the sine of its own power angle δ1; the electromagnetic power of the other machine VSC2 is a2sinδ2, which is only related to the sine of the power angle δ2 of the other machine VSC2; the electromagnetic power between VSC1 and VSC2 is... Only related to the power angle δ between VSC1 and VSC2 12 It is related to the sine. This is a constant term. The equivalent self-damping coefficient D, defined as VSC1 self1 It consists of the following three parts: 1. Damping generated by the output of VSC1: a4cos(θ1-θ1); 2. Damping generated by the output of the power grid: a5cos(θ1-θ1); g 3. Damping generated by the output of VSC2: D self1 The damping term, related to the angular frequency ω1 of the local VSC1, is generated by three different frequency (phase) classifications. The physical meaning of each cosine term is to transform the output of other phases into the coordinate system corresponding to the local phase angle δ1 through PARK transformation, and is therefore defined as self-damped. Similarly, The equivalent mutual damping coefficient D, defined as VSC1 mutu1 It also consists of three parts: 1. Output-related terms of VSC2: a7cos(θ2-θ2); 2. Output-related terms of the power grid: a8cos(θ2-θ2); g ); 3. The output of VSC1 is related to the term a9cos(θ2-θ1).D mutu1 It is related to the angular frequency ω2 of the other machine VSC2.

[0190] The physical meanings of the coefficients in the VSC2 power angle equation (4) are as follows: the mechanical power of VSC2 is b0, which is a constant term; the electromagnetic power of VSC2 to the power grid is b1sinδ2, which is only related to the sine of its own power angle δ2; the electromagnetic power of the other machine VSC1 is b2sinδ1, which is only related to the sine of the power angle δ1 of the other machine VSC1; the electromagnetic power between VSC2 and VSC1 is... Only the power angle δ between VSC2 and VSC1 21 It is related to the sine. This is a constant term. The equivalent self-damping coefficient D, defined as VSC2 self2 It consists of the following three parts: 1. Damping generated by the output of VSC2: b4cos(θ2-θ2); 2. Damping generated by the output of the power grid: b5cos(θ2-θ2); g3. Damping generated by the output of VSC1: D self2 The damping term, related to the angular frequency ω2 of the local VSC2, is generated by three different frequency (phase) classifications. The physical meaning of each cosine term is to transform the output of other phases into the coordinate system corresponding to the local phase angle δ2 through PARK transformation, and is therefore defined as self-damped. Similarly, The equivalent mutual damping coefficient D, defined as VSC2 mutu2 It also consists of three parts: 1. Output-related terms of VSC1: b7cos(θ1-θ1); 2. Output-related terms of the power grid: b8cos(θ1-θ1); g ); 3. VSC2 output related terms b9cos(θ1-θ2).D mutu2 The frequency ω1 of the other machine's angle VSC1 is multiplied by the frequency ω1 to produce damping. It should be noted that equation (4) neglects the angular velocity product term: ω1ω2,ω1 2 ,ω2 2 The simplified models of the frequency coupling terms, cos(2ω1-ω2) and cos(2ω2-ω1), are used. These are relatively small and have limited impact on the dynamic characteristics and stability of the system. However, if they are not ignored, the stability analysis of the parallel converter system becomes quite difficult due to the highly nonlinear characteristics they introduce. Simulation results show that the error is acceptable, such as... Figure 3 As shown.

[0191] S2. The non-integrable electromagnetic interaction terms, the other machine power terms, the self-damped terms, and the mutual-damped terms are scaled down to integrable terms. The positive and negative distribution ranges of damping are also discussed.

[0192] S2.1: Analysis of the positive and negative properties of damping;

[0193] As shown in equation (4), there are two types of damping terms:

[0194] a) Self-damping related to the machine's angular velocity and

[0195] b) Mutual damping related to the angular velocity of other converters The damping term is time-varying depending on the system's state of motion. For example... Figure 4 As shown, in certain areas, D self1 D self2 D mutu1 and D mutu2 The positive or negative value of a molecule can change.

[0196] Figure 4 In the diagram, the dark shaded area represents D. self1and D self2 All values ​​are positive; the light-shaded area represents D. mutu1 and D mutu2 All values ​​are positive. This indicates that under normal transient conditions, the mutual damping is always negative. Meanwhile, the self-damping is positive near the equilibrium point, but may shift towards the negative damping region during transient processes. Since the direction of negative damping during dynamic processes is always consistent with the direction of system motion, it has an accelerating effect, which is detrimental to the transient stability of the system. When analyzing the transient stability of a grid-connected inverter parallel system, the adverse effects of damping on the system cannot be ignored; otherwise, it will lead to misjudgment of stability.

[0197] S2.2: Scaling of electromagnetic interaction terms and self-damping mutual damping;

[0198] As seen in the analysis in S1.2, for a grid-connected inverter parallel system, its power angle dynamics are affected not only by the power angle of the inverter itself but also by the power angles of other inverters. Since the integral of the power angle of other inverters with respect to the angular velocity of the inverter itself is non-integrable, it is necessary to perform appropriate scaling on the interaction terms involving the power angles of other inverters, converting the non-integrable terms into integrable terms that are only related to the power angle of the inverter itself. However, this scaling process inevitably introduces conservatism, and common scaling processes are often based on the boundedness of trigonometric functions. This method uses a less conservative approach based on the monotonicity of trigonometric functions to scale the interaction terms. To further improve conservatism, an iterative algorithm is innovatively proposed, using the stability boundary calculated in the previous iteration as the scaling basis for the next calculation. Simultaneously, for mutual damping, an iterative improved equal-area rule is used, using the maximum angular velocity to scale the mutual damping term, and this is also included in the iterative process, thereby maximizing the improvement in conservatism. This step first scales the corresponding non-integrable terms, and the corresponding iterative process will be performed in S2.3.

[0199] For ease of expression, the interactive electromagnetic power term and the power term of other machines are considered as a whole and defined as P. 1inter At the same time, P 1inter The steady-state value is denoted as P. 1inter_e As shown in (5)-(6):

[0200]

[0201]

[0202] Where δ 1e and δ 2eThey are the steady-state equilibrium points of δ1 and δ2 in (4) respectively. The essence of scaling is to convert non-integrable terms into integrable terms, and the direction of scaling needs to be towards the equations in the system that are prone to instability, that is, to ensure that there is no misjudgment of stability. Specifically, it means scaling the negative damping term in the direction of a larger absolute value and scaling the electromagnetic interaction term in the direction of a smaller value. From the expression of the coefficients, it can be seen that 0 < a3 < a2 ≈ a0 < a1, so for P 1inter There exists a lower bound P 1inter_min :

[0203] P 1inter ≥ a2sin(δ2) - a3 ≥ a2sin(δ 2mine ) - a3 = P 1inter_min (7)

[0204] where δ 2min_e is the stable lower bound of VSC2 calculated according to equations (8)-(9), and the corresponding maximum acceleration area is obtained by considering b0 - P 2inter_e as the equivalent mechanical power while ignoring the damping effect. It should be noted that in the subsequent iterative process, δ 2min_e will be continuously updated,

[0205] δ 2max_e = π - arcsin((b0 - P 2intere ) / b1) (8)

[0206]

[0207] Since a0 - P 1inter_min0 ≥ a0 - P 1inter , when using the equal area rule to analyze the transient stability of the grid-connected inverter parallel system, using P 1inter_min to replace P 1inter will not lead to misjudgment of stability. Taking a0 - P 1inter_min as the mechanical power and ignoring the damping effect, the corresponding stable boundary [δ 1min_int , δ 1max_int can be derived using the equal area method:

[0208] δ 1max_int = π - arcsin((a0 - P 1inter_min ) / a1) (10)

[0209]

[0210] Similarly, for VSC2, there exists a lower bound of the interaction power P 2inter_min :

[0211] P 2inter≥b2sin(δ1)+b3≥b2sin(δ 1min_e )-b3=P 2inter_min (12)

[0212] Similar to (10)-(11), by recognizing b0-P 2inter_min The equivalent mechanical power of VSC2 can be estimated using the equal area method [δ]. 2min_int ,δ 2max_int ].

[0213] However, due to the adverse effect of interactive power on the transient stability of the system, the actual lower boundary of stability δ 2min (Assuming b0-P has time-varying properties) 2inter (For mechanical power) will be located at δ 2min_e The right side. Therefore, assume a0-P 1inter_min The mechanical power will cause the deceleration region of VSC1 to become smaller, therefore the δ obtained from equations (7) and (10)-(11) 1min_int Relatively conservative.

[0214] From a physics perspective, the effects of self-damping and mutual-damping terms on a rotating system are still manifested in the form of consuming the system's energy (or, when negative, inputting energy) through work. That is, they satisfy the form ∫D self1 (δ1,δ2)ω1dδ1 and ∫D mutu1 (δ1,δ2)ω2dδ1. Considering the adverse effect of negative damping on system stability, negative damping cannot be directly ignored but should be regarded as part of the mechanical power, varying with the system's motion state. However, the damping term itself is non-integrable, so it needs to be transformed into an integrable term through mathematical scaling. To avoid misjudgment of stability, the energy exerted on the system by negative damping needs to be scaled to the maximum. According to the expression of the coefficients, a7,b7,a9,b9<0 and a8,b8>0, therefore:

[0215]

[0216] Where δ 1max It is the stable upper boundary of δ1, δ 2max It is the stable upper boundary of δ2. Therefore, the negative work done by the mutual damping can be scaled as shown in (14):

[0217]

[0218] However, (14) still cannot be directly applied to the equal-area method because it involves the other angular frequency ω2. Unless a mapping ω2(δ1) from δ2 to ω1 is solved, the scaling term in (14) remains non-integrable. Since this is equivalent to solving a mapping from δ1 to δ2, it is mathematically infeasible unless the fourth-order differential equation (4) is solved in the time domain. Therefore, an improved equal-area method (IEAC) based on maximum angular frequency scaling is still needed. By taking the maximum value of other VSC angular frequencies, the negative cross-damping can be scaled to:

[0219]

[0220] Where ω 1max_e ω is the maximum angular frequency of VSC1 during the dynamic process described in (8)-(9). 2max_e The maximum angular frequency corresponding to VSC2 in the dynamic process described in equations (8)-(9) is [δ]. However, it should be noted that the actual stability boundary is [δ]. 1min_e ,δ 1max_e ] and [δ 2min_e ,δ 2max_e It is a subset of ], therefore the actual maximum angular frequency is less than ω. 1max_e and ω 2max_e However, this conservatism can be minimized through subsequent iterations A.

[0221]

[0222] For self-damped systems, since a4,b4<0 and a5,b5,a6,b6>0, the self-damping factor D is... self1 and D self2 It can be scaled down to D. s1_min (δ1) and D s2_min (δ2):

[0223]

[0224] It should be noted that D s1_min (δ1) is a univariate function of δ1, D s2_min (δ2) is a single-variable function of δ2, therefore the negative work done by self-damping can be scaled as follows (18):

[0225]

[0226] Thus, all non-integrable terms have been transformed into integrable terms through corresponding mathematical scaling. Theoretically, by combining equations (7), (12), (15), and (17), the transient stability boundary of the grid-connected inverter parallel system can be calculated using the equal area method. However, the results obtained in this way are quite conservative. The following steps will introduce how to improve the conservatism through the proposed double iteration method.

[0227] S3. To address the scaling issue in S2, a dual-iteration equal-area method is proposed. This involves two iterative processes: Iteration A handles the electromagnetic interaction power, including the electromagnetic power and mutual damping terms; Iteration B handles the self-damping terms. Iteration A yields an estimate of the upper boundary of the power angle stability, and Iteration B, based on Iteration A, yields an estimate of the lower boundary of the power angle stability.

[0228] S3.1: S2.1 scaled down the non-integrable interaction and damping terms, making the application of the equal-area method in grid-connected inverter parallel systems possible. However, it should be noted that the scaling process undoubtedly introduces a certain degree of conservatism. Iterative methods can effectively improve conservatism by narrowing the scaling range. As described in S2.1, the non-integrable terms that need to be scaled mainly include the following three parts: interactive electromagnetic power terms and other machine power terms; mutual damping terms; and self-damping terms. A dual-iteration equal-area method is proposed to calculate the transient stability boundary of grid-connected inverter parallel systems, which mainly includes two iterative processes: Iteration A and Iteration B.

[0229] 1) Iteration A: Iterate over the interactive electromagnetic power, the other electromagnetic power, and the mutual damping term to obtain the upper boundary δ. 1max and δ 2max The conservative estimate provides a basis for subsequent calculations of δ. 1min and δ 2min Laying the foundation. The purpose of iteration A is to reduce the scaling range through iteration, thereby improving δ in equations (7), (12), and (16). 1min_e ,δ 2min_e ,ω 1max_e and ω 2max_e Conservatism caused by contraction.

[0230] 2) Iteration B: An extended iterative equal-area rule (Extend ITEAC) is proposed and used to obtain a stable lower boundary considering self-damping. It should be noted that the single-machine iterative equal-area rule does not introduce conservatism. The conservatism of the extended iterative equal-area rule is due to the scaling of self-damping in S2.1. The purpose of the extended iterative equal-area rule is to determine a corresponding D through an iterative method. s1_min (δ1) and δ 1max The frequency-power angle mapping relationship ω1(δ1) and the corresponding frequency zero-crossing point (i.e., the lower boundary of power angle stability) δ 1min Thus, a quantitative and non-conservative analysis of D s1_min (δ1) The impact on the transient stability of the system. A corresponding D is determined using an iterative method. s2_min (δ2) and δ 1max The frequency-power angle mapping relationship ω2(δ2) and the corresponding frequency zero-crossing point (i.e., the lower boundary of power angle stability) δ 2min Thus, a quantitative and non-conservative analysis of Ds2_min (δ2) The effect on the transient stability of the system.

[0231] Iteration A considered the effects of the interactive electromagnetic power term and the mutual damping term, and obtained the upper stable boundary δ. 1max and δ 2max The estimate. Iteration B is based on the δ obtained from iteration A. 1max and δ 2max A proposed and utilized extended iterative equal-area method to calculate the lower boundary δ 1min and δ 2min .

[0232] S3.2: Explanation of the process of iteration A.

[0233] The main content of iteration A is to obtain a stable boundary estimate [δ] during the j-th iteration. 1minA_j ,δ 1maxA_j ] and [δ 2minA_j ,δ 1maxA_j This is used as the basis for scaling in the next iteration, thereby minimizing the conservatism caused by the scaling of electromagnetic interaction terms, electromagnetic power, and mutual damping terms.

[0234] First, define P. 1A_j As the electromagnetic interaction power of VSC1 after scaling in the j-th iteration A, it is the sum of its mechanical power and mutual damping; define P 2A_j As the electromagnetic interaction power of VSC2 after scaling in the j-th iteration A, it is the sum of its mechanical power and mutual damping:

[0235]

[0236] Where ω 1max_j-1 ,δ 1minA_j-1 andδ 1imaxA_j-1 These are the maximum angular frequency of VSC1, the stable lower boundary, and the stable right boundary, respectively, in the (j-1)th iteration; ω 2max_j-1 ,δ 2minA_j-1 andδ 2imaxA_j-1 These represent the maximum angular frequency of VSC2, the lower stable boundary, and the right stable boundary, respectively, in the (j-1)th iteration. Their calculation process in the j-th iteration is shown below:

[0237]

[0238]

[0239]

[0240] It should be noted that P in the first iteration 1A_1 and P 2A_1 The calculation process is as follows:

[0241]

[0242] From the above derivation, it can be seen that a new scaling reference can be obtained from (20)-(22) and used cyclically in scaling (19). If δ 1maxA_j and δ 2maxA_j If the iteration converges to a given precision ε, then the iteration ends, and the stable boundary on δ1 can be taken as δ. 1maxA The stable boundary on δ2 can be taken as δ 2maxA Otherwise, assume j = j + 1, and begin the (j+1)th iteration A. The flowchart for iteration A is as follows: Figure 5 As shown. The final output convergence value δ 1maxA and δ 2maxA The proposed dual-iteration EAC method serves as the final upper bound δ. 1max and δ 2max Although δ was not included 1minA and δ 2minA As the final value of the lower boundary, but it will be used to calculate the initial value in subsequent iteration B. Compared with the existing non-iterative scaling methods, the proposed method not only adopts a more accurate mathematical model that considers frequency fluctuations, but also adopts an iterative algorithm, which greatly reduces the conservatism brought about by the scaling of interaction terms.

[0243] S3.2: Explanation of the iterative process B.

[0244] The scaled-down self-damping term has a similar mathematical structure to the damping term in a single-converter system, i.e., it is only related to the machine's power angle and angular velocity. Therefore, the ITEAC method proposed in the single-converter system can naturally be extended to the parallel-converter system. The maximum power angle δ obtained by iteration (19)-(23) is... 1maxA ,δ 2maxA and interaction power P 1A ,P 2A This prepares the ground for the extended ITEAC (i.e., Iteration B). The extended ITEAC primarily utilizes angular frequency ω. 1B (δ1) and ω 2B The iteration of (δ2) calculates where the system originates from the self-damping term D after scaling. s1_min and D s2_min The following can just move to δ 1maxA and δ 2maxA .

[0245] First, define the initial self-damping power and the initial angular frequency ω of the system. 1B_0 (δ1) and ω 2B_0 (δ2), where δ1 and δ2 in parentheses indicate that the required angular frequency ω is... 1B_0 and ω 2B_0Distribution functions relative to δ1 and δ2 respectively:

[0246]

[0247]

[0248] angular velocity ω 1B_0 (δ1) and ω 2B_0 (δ2) is calculated using the variable lower limit integral as shown in (25). Considering D s1_min Independent of ω1, D s2_min It is independent of ω2, so it does not need to be updated in each iteration.

[0249] In the j-th iteration, the self-damping power P of VSC1 1B_j Based on the angular velocity distribution ω at the (j-1)th iteration 1B_j-1 (δ1) and minimum work angle δ 1minB_j-1 The self-damping power P of VSC2 is calculated. 2B_j Based on the angular velocity distribution ω at the (j-1)th iteration 2B_j-1 (δ2) and minimum work angle δ 2minB_j-1 Calculation

[0250]

[0251] According to the principle of the equal area method, the angular velocity distribution ω of VSC1 obtained in the j-th iteration is... 1B_j Angular velocity distributions ω of (δ1) and VSC2 2B_j (δ2) are shown below:

[0252]

[0253] Meanwhile, the lower boundary δ of the VSC1 power angle stability domain at the j-th iteration is... 1minB_j The lower boundary δ of the VSC2 power angle stability region 2minB_j It can also be calculated as:

[0254]

[0255] The flowchart of iteration B is shown below. Figure 6 As shown, equations (24)-(25) give the initial value P for the iteration. iB_0 and ω iB_0 (δ i The calculation formula is as follows. In the j-th cycle, the angular velocity distribution ω of VSC1 is calculated using equations (27)-(28). 1B_j (δ1) and the lower boundary of the work angle stability δ 1minB_j and the angular velocity distribution ω of VSC2 2B_j (δ2) and the lower boundary of the work angle stability δ2minB_j The calculation results of equations (27)-(28) are used as inputs for the (j+1)th iteration and substituted into equation (26) to calculate the self-damping power P at the (j+1)th iteration. 1B_j and P 2B_j If the δ obtained by equation (28) 1minB_j and δ 2minB_j If the iteration converges to a given precision ε2, then the iteration ends, and the lower boundaries of δ1 and δ2 can be written as δ 1min and δ 2min The iterative process shown in (24)-(38) is defined as the extended iterative equal area rule, which is iteration B in this embodiment.

[0256] The dual-iteration EAC method in this embodiment can accurately estimate the power angle stability boundary of the grid-connected inverter parallel system: [δ 1min ,δ 1max ] and [δ 2min ,δ 2max The relationship between the two pseudo-iterative processes is as follows: Figure 7 As shown. Iteration A yields a stability boundary estimate considering the effects of interactive electromagnetic power, other electromagnetic power, and mutual damping terms, denoted as [δ]. 1minA ,δ 1maxA ] and [δ 2minA ,δ 2maxA ]. Where δ 1maxA and δ 2maxA The upper bound δ of the proposed dual-iteration EAC method is directly used as the upper bound. 1max and δ 2max Therefore, it will not be updated in subsequent iterations B. Furthermore, the result of iteration A also serves as the basis for calculating the initial value of iteration B, as shown in equation (25). Based on iteration A, considering the influence of the self-damping term, iteration B obtains the lower boundary δ using the extended iterative equal-area method. 1minB and δ 1minB And serves as the lower bound δ for the proposed dual-iteration EAC method. 1min and δ 2min Combining iterations A and B, the final stable boundary [δ] is obtained. 1min ,δ 1maxA ] and [δ 2min ,δ 2max ].

[0257] The above are merely preferred embodiments of the present invention and are not intended to limit the implementation methods and protection scope of the present invention. Those skilled in the art should recognize that any equivalent substitutions and obvious changes made based on the content of this specification should be included within the protection scope of the present invention.

Claims

1. A transient stability analysis method for multi-inverter based on the two-layer iterative equal-area rule, characterized in that: Includes the following steps: Step 1: Establish a fourth-order mathematical model of the grid-connected inverter parallel system considering line frequency fluctuations, including mechanical torque, electromagnetic torque, electromagnetic interaction torque, self-damped torque, and mutual damped torque. Step 2: Scaling the non-integrable electromagnetic interaction terms, other machine power terms, self-damping and mutual damping into integrable terms, and dividing the positive and negative distribution intervals of damping. Step 3: Based on the scaling in Step 2, the double-iteration equal-area method is obtained. The double-iteration equal-area method includes iteration A and iteration B, where iteration A handles electromagnetic interaction power, other machine electromagnetic power, and mutual damping terms; iteration B handles self-damping terms; iteration A obtains the estimated upper boundary of the power angle stability, and iteration B obtains the estimated lower boundary of the power angle stability based on iteration A.

2. The transient stability analysis method for multi-inverter based on the two-layer iterative equal-area rule according to claim 1, characterized in that: Step 1 includes the following steps: Step 1.1: Establish a nonlinear model of the grid-connected inverter parallel system; The parallel grid-connected inverter VSC system includes: VSC1 with filter impedance L f1 Then through line impedance L 1, R 1 and grid impedance L g , R g With the power grid V g Connected; VSC2 at the filter impedance L f2 Then through line impedance L 2, R 2 and grid impedance L g , R g With the power grid V g Connected; both VSC1 and VSC2 use direct current control and phase-locked loop (PLL); According to Kirchhoff's voltage law, the voltage at the grid connection point of VSC1 is... V t1 Voltage at the VSC2 grid connection point V t2 They are represented as follows: (1) in φ 1 represents the power factor of VSC1. φ 1 = arctan( i g1,q / i g1,d ), i g1,d =I 1 cosφ , i g1,q =I 1 sinφ Line current I The dq-axis component of 1; φ 2 represents the power factor of VSC2. φ 2 = arctan( i g2,q / i g2,d ), i g2,d =I 2 cosφ , i g2,q =I 2 sinφ Line current I 2 dq axis components; θ PLL1 It is the reference phase of the phase-locked loop output of VSC1. θ PLL2 It is the reference phase of the phase-locked loop output of VSC1; φ 1 =φ 2 = 0, applying the PARK transformation to both sides of equation (1) yields the following results: V t1 q-axis components V t1q ,and V t2 q-axis components V t1q quantity: (2) The dynamic equation of the phase-locked loop is: (3) definition δ 1= θ PLL1 - θ g , is the power angle of VSC1; defined δ 2= θ PLL2 - θ g To obtain the dynamic equations of the grid-connected inverter parallel system, we can combine equations (2) and (3) to get the power angle of VSC2. (4) in M 1 represents the equivalent inertia coefficient of VSC1. a 0 represents the equivalent mechanical power of VSC1. a 1 represents the electromagnetic power between VSC1 and the power grid; a 2 represents the electromagnetic power between VSC2 and the power grid; a 3 represents the interactive electromagnetic power between VSC1 and VSC2; a 4+ a 5cos( δ 1)+ a 6cos( δ 1- δ 2+ φ D11 )) ω 1 indicates that VSC1 is self-damped; a 7+ a 8cos( δ 2) + a 9cos( δ 2- δ 1+ φ D12 )) ω 2 indicates the mutual damping of VSC1; M 2 represents the equivalent inertia coefficient of VSC2. b 0 represents the equivalent mechanical power of VSC2. b 1 represents the electromagnetic power between VSC2 and the power grid; b 2 represents the electromagnetic power between VSC1 and the power grid; b 3 represents the interactive electromagnetic power between VSC2 and VSC1; b 4+ b 5cos( δ 2)+ b 6cos( δ 2- δ 1+ φ D21 )) ω 2 indicates that VSC2 is self-damped; b 7+ b 8cos( δ 1)+ b 9cos( δ 1- δ 2+ φ D22 )) ω 1 indicates the mutual damping of VSC2; The expressions for each coefficient are as follows: Step 1.2: Compare the VSC parallel system and the SG system to obtain the equivalent mechanical power, electromagnetic power, electromagnetic interaction power, self-damping, and mutual damping. From the VSC1 power angle equation (4), we know that the mechanical power of VSC1 is a The term 0 is a constant term, and the electromagnetic power term of VSC1 to the power grid is... a 1sin δ 1. Only with the self-power angle δ 1 is related to the sine wave; its VSC2 electromagnetic power a 2sin δ 2. Only the power angle of VSC2 of other machines. δ The sine wave of 2 is related; the inter-machine electromagnetic power of VSC1 and VSC2. a 3sin( δ 12 + φ E1 Only the power angle between VSC1 and VSC2 δ 12 It is related to the sine. φ E1 For constant terms; a 4+ a 5cos( δ 1)+ a 6cos( δ 1- δ 2+ φ D11 The equivalent self-damping coefficient of VSC1 is defined as follows: D self1 This includes the damping generated by the output of VSC1. a 4cos( θ 1- θ 1) Damping generated by the power grid output a 5cos( θ 1- θ g Damping generated by the output of VSC2 a 6cos( θ 1 -θ 2+ φ D11 The equivalent self-damping coefficient of VSC1 D self1 It is the angular frequency of the local VSC1 generated by the classification of three different frequencies or phases. ω 1. Regarding the damping terms, the physical meaning of each cosine term is to convert the output of other phases into the local phase angle through PARK transformation. δ In the coordinate system corresponding to 1, it is defined as self-damped; a 7+ a 8cos( δ 2)+ a 9cos( δ 2- δ 1+ φ D12 The equivalent mutual damping coefficient of VSC1 is defined as... D mutu1 This includes items related to the output of VSC2. a 7cos( θ 2- θ 2); Items related to power grid output a 8cos( θ 2- θ g ); VSC1 output related items a 9cos( θ 2- θ 1) Equivalent mutual damping coefficient of VSC1 D mutu1 Angular frequency of VSC2 and other machines ω 2. Related to; From the VSC2 power angle equation (4), we know that the mechanical power of VSC2 is b The term 0 is a constant term, and the electromagnetic power term of VSC2 to the power grid is... b 1sin δ 2. Only with the self-power angle δ The sine wave of 2 is related; its VSC1 electromagnetic power b 2sin δ 1. Only the power angle of VSC1 of other machines. δ The sine wave of 1 is related; the inter-machine electromagnetic power of VSC2 and VSC1. b 3sin( δ 21 + φ E2 Only the power angle between VSC2 and VSC1 δ 21 It is related to the sine. φ E2 For constant terms; b 4+ b 5cos( δ 2)+ b 6cos( δ 2- δ 1+ φ D21 The equivalent self-damping coefficient of VSC2 is defined as... D self2 This includes the damping generated by the output of VSC2. b 4cos( θ 2- θ 2) Damping generated by the power grid output b 5cos( θ 2- θ g Damping generated by the output of VSC1 b 6cos( θ 2 -θ 1+ φ D21 The equivalent self-damping coefficient of VSC2 D self2 It is the angular frequency of the local VSC2 generated by the classification of three different frequencies or phases. ω 2. Regarding the damping terms, the physical meaning of each cosine term is to convert the output of other phases into the local phase angle through PARK transformation. δ In the coordinate system corresponding to 2, it is defined as self-damped; b 7+ b 8cos( δ 1)+ b 9cos( δ 1- δ 2+ φ D22 The equivalent mutual damping coefficient of VSC2 is defined as... D mutu2 This includes items related to the output of VSC1: b 7cos( θ 1- θ 1) Items related to power grid output b 8cos( θ 1- θ g ); VSC2 output related items b 9cos( θ 1- θ 2); Equivalent self-damping coefficient of VSC2 D mutu2 Frequency of VSC1 and its other angle ω 1. Multiplication produces damping.

3. The transient stability analysis method for multi-inverter based on the two-layer iterative equal-area rule according to claim 2, characterized in that: Step 2 includes the following steps: Step 2.1, Damping Positive and Negative Analysis; From Equation (4), we know that: a) Self-damping related to the machine's angular velocity ω 1 D self1 = ω 1 [ a 4+ a 5cos δ 1+ a 6cos( δ 1- δ 2+ φ D11 )]and D self2 = ω 2 [ b 4+ b 5cos( δ 2)+ b 6cos( δ 2- δ 1+ φ D21 )]; b) Mutual damping related to the angular velocity of other converters ω 2 D mutu1 = ω 2[ a 7+ a 8cos δ 2+ a 9cos( δ 2- δ 1+ φ D12 )], ω 1 D mutu2 = ω 1[ b 7+ b 8cos( δ 1)+ b 9cos( δ 1- δ 2+ φ D22 The damping term varies with the system's motion state over time. D self1 , D self2 , D mutu1 and D mutu2 Its positive or negative value can change; During the transient process under normal conditions, the mutual damping is always negative; at the same time, the self-damping is positive near the equilibrium point and will move towards the negative damping region during the transient process; the direction of the negative damping in the dynamic process is always consistent with the direction of the system motion, which plays an accelerating role. When analyzing the transient stability of the grid-connected inverter parallel system, the influence of damping on the system should be ignored. Step 2.2: Scaling up the electromagnetic interaction term and the self-damping mutual damping; First, the corresponding non-integrable terms are scaled; the interaction terms containing the power angle of other machines are scaled accordingly, converting the non-integrable terms into integrable terms that are only related to the power angle of the machine itself; the interaction terms are scaled using a method based on the monotonicity of trigonometric functions; an iterative algorithm is used, with the stability boundary calculated in the previous calculation used as the scaling basis for the next calculation; for mutual damping, an iterative improved equal area rule is used, and the maximum angular velocity is used to scale the mutual damping terms, while incorporating them into the iterative process. The interactive electromagnetic power term and the power term of other machines are considered as a whole and defined as follows: P 1inter ;at the same time P 1inter The steady-state value is denoted as P 1inter_e As shown in (5)-(6): (5) (6) in δ ie It is (4) δ i The steady-state equilibrium point; from the coefficient expression in step 1.1, we know 0 <a 3< a 2≈ a 0< a 1. For P 1inter There exists an infimum P 1inter_min : (7) in δ 2min_e The stable left boundary of VSC2 is calculated as shown in equations (8)-(9), and its corresponding maximum acceleration area is determined by assuming... b 0- P 2inter_e This represents the equivalent mechanical power, while neglecting the effect of damping; it will be continuously updated in subsequent iterations. δ 2min_e This allows us to obtain a more accurate stability region estimate. (8) (9) because a 0- P 1inter_min ≥ a 0- P 1inter Therefore, when using the equal area rule to analyze the transient stability of a grid-connected inverter parallel system, the following method is used: P 1inter_min replace P 1inter It will not lead to misjudgment of stability; a 0- P 1inter_min Assuming mechanical power and neglecting damping effects, the corresponding stability boundary is derived using the equal-area method. δ 1min_int , δ 1max_int Estimate: (10) (11) Similarly, there exists an infimum of interactive power for VSC2. P 2inter_min : (12) Similar to equations (10)-(11), b 0- P 2inter_min To estimate the equivalent mechanical power of VSC2, use the equal area method to estimate a stability boundary of VSC2. δ 2min_int , δ 2max_int ]; Due to the adverse effects of interactive power on the transient stability of the system, consider a system with time-varying characteristics. b 0- P 2inter For mechanical power, the actual stable left boundary δ 2min Will be located δ 2min_e On the right side; let a 0- P 1inter_min The mechanical power will cause the deceleration range of VSC1 to become smaller, as obtained from equations (7) and (10)-(11). δ 1min_ int Relatively conservative; The effects of self-damping and mutual-damping terms on a rotating system are either the consumption of the system's energy in the form of work, or, when negative, the input capacity; satisfying the form ∫ D self1 ( δ 1, δ 2) ω 1d δ 1 and ∫ D mutu1 ( δ 1, δ 2) ω 2d δ 1; Treat negative damping as part of the mechanical power, which varies with the system's motion state; the damping term itself is non-integrable, scaling the energy exerted on the system by negative damping to the maximum; according to the coefficient expression in step 1.1: a 7, b 7, a 9, b 9<0 and a 8, b 8>0, therefore: (13) in δ imax yes δ i The negative work done by the mutual damping at the stable upper boundary is scaled as shown in (14): (14) An improved equal-area rule based on scaling by the maximum angular frequency is adopted. IEAC scales the negative cross-damping to equation (15) by taking the maximum value of the angular frequencies of other VSCs: (15) in ω imax_e The maximum angular frequency corresponding to VSC1 or VSC2 in the dynamic process described in equations (8)-(9); the actual stability boundary is [ δ imin_e , δ imax_e A subset of ], therefore the actual maximum angular frequency is less than ω imax_e This conservatism can be minimized through subsequent iterations A. (16) For self-damping D selfi Due to the existence a 4, b 4<0, a 5, b 5, a 6, b If 6 > 0, then D selfi Shrink to D si_min ( δ i ): (17) D si_min ( δ i )yes δ i It is a single-variable function, therefore self-damped. D self1 The negative work done can be scaled as shown in equation (18): (18) Combining equations (7), (12), (15), and (17), the transient stability boundary of the grid-connected inverter parallel system can be calculated using the equal area method.

4. The transient stability analysis method for multi-inverter based on the two-layer iterative equal-area rule according to claim 2, characterized in that: Step 3 includes the following steps: Step 3.

1. The double-iteration equal-area method is used to calculate the transient stability boundary of the grid-connected inverter parallel system, including Iteration A and Iteration B; Iteration A, considering the effects of interactive electromagnetic power and mutual damping terms, yields the upper stable boundary. δ 1max and δ 2max The valuation; Iteration B is based on the value obtained from iteration A. δ 1max and δ 2max The lower boundary is calculated using the extended iterative equal area method. δ 1min and δ 2min ; Step 3.

2. Iteration A includes: A stability boundary estimate is obtained during the j-th iteration. δ 1minA_j , δ 1maxA_j ]and[ δ 2minA_j , δ 1maxA_j This will be used as the basis for scaling in the next iteration; definition P 1A_j The electromagnetic interaction power of VSC1 after scaling in the j-th iteration A is the sum of its mechanical power and mutual damping; defined as follows. P 2A_j As the electromagnetic interaction power of VSC2 after scaling in the j-th iteration A, it is the sum of its mechanical power and mutual damping: (19) in ω 1max_j-1 , δ 1minA_j-1 and δ 1imaxA_j-1 They are the first j The maximum angular frequency of VSC1 in -1 iterations, the stable lower boundary and the stable right boundary; ω 2max_j-1 , δ 2minA_j-1 and δ 2imaxA_j-1 They are the first j The maximum angular frequency of VSC2 in the -1st iteration, the stable lower boundary and the stable right boundary; in the... j The calculation process in this iteration is as follows: (20) (21) (22) During the first iteration P 1A_1 and P 2A_1 The calculation process is as follows: (23) From the above derivation, it can be seen that a new scaling reference is obtained from (20)-(22) and is used cyclically in scaling (19); if δ 1maxA_j and δ 2maxA_j If the iteration converges to a given precision ε, then exit the iteration. δ The stable boundary can be taken as 1. δ 1maxA , δ The stable boundary can be taken as 2. δ 2maxA Otherwise, assume j = j +1, and begin the first j + 1 iteration A; final output convergence value δ 1maxA and δ 2maxA The proposed dual-iteration EAC method serves as the final upper bound. δ 1max and δ 2max ; Step 3.

3. Iteration B includes: Define the initial self-damped power and initial angular frequency of the system. ω 1B_0 ( δ 1) with ω 2B_0 ( δ 2), in parentheses δ 1 and δ 2 indicates that the required frequency is angular frequency. ω 1B_0 and ω 2B_0 Relative to δ 1 and δ Distribution function of 2: (24) (25) angular velocity ω 1B_0 ( δ 1) with ω 2B_0 ( δ 2) Calculate using the variable lower limit integral as shown in (25); considering D s1_min and ω 1. Irrelevant D s2_min and ω 2 is irrelevant, therefore it does not need to be updated in each iteration; In the j In this iteration, the self-damping power of VSC1 P 1B_j Based on the angular velocity distribution at the (j-1)th iteration ω 1B_j-1 ( δ 1) and minimum work angle δ 1minB_j-1 The self-damping power of VSC2 is calculated. P 2B_j Based on the angular velocity distribution at the (j-1)th iteration ω 2B_j-1 ( δ 2) and minimum work angle δ 2minB_j-1 The calculation shows that: (26) According to the principle of the equal area method, the first... j The angular velocity distribution of VSC1 obtained in the next iteration ω 1B_j ( δ 1) Angular velocity distribution of VSC2 ω 2B_j ( δ 2) As shown below: (27) Meanwhile, the lower boundary of the VSC1 power angle stability region at the j-th iteration. δ 1minB_j and the lower boundary of the VSC2 power angle stability region 2minB_j It can also be calculated as: (28) Equations (24)-(25) give the initial values ​​for the iteration. P iB_0 and δ iB_0 ( ω i The calculation formula for ); the first j In the next iteration, (27)-(28) are used to calculate the angular velocity distribution of VSC1. δ 1B_j ( ω 1) Lower boundary of the sum of work angles for stability δ 1minB_j and the angular velocity distribution of VSC2 δ 2B_j ( ω 2) Lower boundary of the work angle stability δ 2minB_j The calculation results of equations (27)-(28) are used as the first... j Input at +1 iteration, substitute into (26) to calculate j Self-damping power at +1 iteration P 1B_j and P 2B_j If the result obtained by equation (28) δ 1minB_j and 2minB_j Converges to a given precision δ 2. Then exit the iteration. 1 and δ The lower boundary of 2 is written as follows: 1min and ε 2min The iterative process of equations (24)-(38) is the second iteration B of the extended iterative equal area rule; The power angle stability boundary of the grid-connected inverter parallel system is obtained through double iteration: [ δ 1min , δ 1max ]as well as[ δ 2min , δ 2max Iteration A yielded a stability boundary estimate considering the effects of interactive electromagnetic power, other electromagnetic power, and mutual damping terms, denoted as [ δ 1minA , δ 1maxA ]as well as[ δ 2minA , δ 2maxA ];in δ 1maxA and δ 2maxA As the upper bound of the double iteration δ 1max and δ 2max The results of iteration A will not be updated in subsequent iterations B; the results of iteration A will serve as the basis for the initial value calculation of iteration B; based on iteration A, considering the influence of the self-damping term, iteration B will obtain the lower boundary using the extended iterative equal-area method. δ 1minB and δ 1minB As the lower bound of the double iteration δ 1min and δ 2min Combining iterations A and B, the final stable boundary is obtained. δ 1min , δ 1maxA ]as well as[ δ 2min , δ δ δ δ δ 2max ].

Images