Transient stability analysis method for grid-connected system

Abstract

The invention provides a transient stability analysis method for a grid-connected system. The method comprises the following steps: establishing a mathematical model of a heterogeneous new energy grid-connected system comprising a network construction type new energy unit, a network following type new energy unit, a synchronous generator and a receiving end power grid, and defining a state variable set of the system; solving a stable equilibrium point and an unstable equilibrium point of the system according to the operation condition of the system; and by taking the unstable equilibrium point as a starting point, calculating by adopting a numerical integration method or a numerical continuation method to obtain a stable manifold starting from the unstable equilibrium point. The method provided by the invention realizes global and accurate depiction of the transient stability domain, overcomes the defects that a traditional small signal linearization method can only analyze local stability and a traditional energy function method depends on an artificial construction function and is limited in applicability, fundamentally changes an analysis normal form of the transient stability domain, and improves the analysis accuracy of the transient stability domain. And a more effective analysis tool is provided for safe and stable operation of the high-permeability new energy grid-connected system.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of power system stability analysis and control, in particular to a grid-connected system transient stability analysis method. BACKGROUND

[0002] With the rapid growth of new energy installed capacity, grid-forming (GFM) and grid-following (GFL) new energy coexist in the power system. Grid-forming new energy provides reference frequency and voltage through virtual synchronous generator (VSG) control and has inertia support capability; grid-following new energy mainly injects power by phase-locked loop (PLL) tracking. Under the condition of high new energy penetration, the traditional transient stability analysis method cannot accurately describe the real stability region of the system.

[0003] At present, the existing transient stability evaluation methods mainly include two types: (1) Small signal linearization analysis, taking the stable equilibrium point (SEP) of the system as the center, the local stability is evaluated by solving the eigenvalues of the Jacobian matrix. This method is suitable for weakly nonlinear systems, but it cannot describe the global boundary of the stability domain, and it is difficult to reflect the nonlinear influence of control parameters and power ratio changes on stability. (2) Energy function method and simulation analysis, by comparing the energy changes of different working conditions through Lyapunov energy function or transient energy function, or by observing the phase trajectory through time domain simulation. The energy function method is mainly used for single-machine or simplified systems, and the precision is limited by the construction of the energy function; although the full electromagnetic transient simulation can show the phase trajectory, it is difficult to accurately determine the stability domain boundary and cannot be used for rapid parameter sensitivity analysis.

[0004] Therefore, the present application provides a grid-connected system transient stability analysis method to solve one of the above technical problems. SUMMARY

[0005] The purpose of the present application is to provide a grid-connected system transient stability analysis method, which can solve at least one of the above technical problems. The specific scheme is as follows: According to the specific embodiment of the present application, in a first aspect, the present application provides a grid-connected system transient stability analysis method, comprising: S1: establishing a mathematical model of a heterogeneous new energy grid-connected system containing grid-forming new energy units, grid-following new energy units, synchronous generators and receiving end power grids, and defining the state variable set of the system; S2: solving the stable equilibrium point and unstable equilibrium point of the system according to the system operating condition; S3: taking the unstable equilibrium point as the starting point, using numerical integration method or numerical continuation method for calculation, and obtaining the stable manifold from the starting point; wherein the stable manifold is used to describe the transient stability domain boundary of the system.

[0006] In an embodiment, in the S1, the grid-connection type new energy unit adopts a virtual synchronous generator control strategy, and a mathematical model is expressed based on the following formula: ; ; wherein, is a virtual inertia, is a mechanical input power, is an electromagnetic output power, is an angle, is an actual angular velocity of a rotor, is a reference angular velocity, and · denotes a first-order derivative of a corresponding parameter.

[0007] In an embodiment, in the S1, the grid-following type new energy unit is connected to the grid through a phase-locked loop, and a mathematical model of the phase-locked loop is expressed by the following formula: ; ; wherein, is a phase angle output by the phase-locked loop, is an angular velocity output by the phase-locked loop, is a grid phase angle, is a grid connection point voltage amplitude, and are a proportional coefficient and an integral coefficient respectively, and · denotes a first-order derivative of a corresponding parameter.

[0008] In an embodiment, in the S2, after the stable equilibrium point is solved, the method further comprises: performing first-order linearization on a system nonlinear equation at the stable equilibrium point to obtain a Jacobian matrix; and determining a type of the equilibrium point by calculating eigenvalues of the Jacobian matrix; wherein the type of the equilibrium point comprises a stable node, a saddle point or a center point.

[0009] In an embodiment, in the S3, after the stable manifold is calculated, the method further comprises: performing normalization processing on the stable manifold under different operating conditions, and analyzing an influence of system parameters on the transient stability domain by comparing morphological differences of different stable manifolds; wherein the morphological differences comprise at least one of an area surrounded by the stable manifold in a phase plane, a stretching direction and a convergence rate; and wherein the system parameters comprise a power ratio of the grid-connection type unit and the grid-following type unit, a virtual inertia and a damping coefficient.

[0010] In an embodiment, after the S3, the method further comprises: S4: constructing a scalar function for quantitatively describing a scale feature of the stable manifold, and determining a scale size of the stable domain in different state variable directions by analyzing a sign and a distribution of a change rate of the scalar function over time.

[0011] In an embodiment, the scalar function is a quadratic function, and is expressed by the following formula: ; wherein, is a scalar function, x1 and x2 are state variables, and k1 and k2 are weight coefficients.

[0012] In an implementation, the method further comprises: fitting the geometric shape of the stable manifold by adjusting the values of k1 and k2, and diagnosing the dominance of the power angle direction and the angular velocity direction on the stable domain boundary; wherein increasing k1 causes the equipotential lines of the scalar function to be pinched in the power angle direction, and increasing k2 causes the equipotential lines to be pinched in the angular velocity direction.

[0013] In an implementation, the values of the weight coefficients k1 and k2 are combined as (4, 1) or (8, 2).

[0014] In an implementation, the set of state variables includes the power angle, the angular velocity, the virtual synchronous machine power angle, the phase-locked loop phase, and the grid point voltage amplitude of each unit.

[0015] The above scheme of the embodiments of the present application has at least the following beneficial effects compared with the prior art: The present application provides a grid-connected system transient stability analysis method, which realizes global and accurate characterization of the transient stability domain by directly calculating the stable manifold starting from the unstable equilibrium point, overcomes the defects of the traditional small-signal linearization method that can only analyze local stability and the traditional energy function method that relies on artificial construction of functions and has limited applicability; at the same time, the method converts the abstract stability problem into intuitive visual geometric figures in the phase plane, greatly enhancing the understanding of system stability and the interpretability of analysis results; in addition, the obtained accurate stability domain boundary lays a reliable foundation for subsequent quantitative researches such as parameter sensitivity analysis and control strategy comparison, fundamentally changing the analysis paradigm of the transient stability domain, and providing a more effective analysis tool for the safe and stable operation of high-penetration new energy grid-connected systems. BRIEF DESCRIPTION OF DRAWINGS

[0016] Figure 1 A flowchart of a grid-connected system transient stability analysis method is shown; Figure 2 A schematic diagram of a heterogeneous new energy grid-connected system model is shown; Figure 3 A flowchart of another grid-connected system transient stability analysis method is shown; Figure 4 A schematic diagram of the attraction domain variation law of a heterogeneous system under different new energy ratios is shown; Figure 5 A schematic diagram of the attraction domain variation law of a heterogeneous system under different new energy ratios is shown; Figure 6A schematic diagram showing the influence of inertia parameters of grid-forming converters on the system flow pattern is shown. Figure 7 A schematic diagram showing the influence of damping parameters of grid-forming converters on the system flow pattern is shown. Figure 8 A graph showing the energy variation corresponding to the scalar function expression is shown. Figure 9 A schematic diagram showing the flow pattern analysis under different control strategies is shown. Figure 10 A schematic diagram showing the flow scale depiction is shown. DETAILED DESCRIPTION

[0017] In order to make the purposes, technical solutions and advantages of the present application clearer, the present application will be further described in detail below with the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present application, but not all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative work fall within the scope of protection of the present application.

[0018] The terms used in the embodiments of the present application are merely for the purpose of describing the specific embodiments, and are not intended to limit the present application. The singular forms “a”, “an” and “the” used in the embodiments of the present application and the appended claims are also intended to include the plural forms, unless the context clearly indicates otherwise. “Plural” generally includes at least two.

[0019] It should be understood that the term “and / or” used herein is merely an association relationship between the associated objects, which means that there can be three relationships, for example, A and / or B can represent the three cases of A alone, A and B together, and B alone. In addition, the character “ / ” herein generally represents an “or” relationship between the front and rear associated objects.

[0020] It should be understood that although the terms first, second, third, etc. may be used in the embodiments of the present application, these descriptions should not be limited to these terms. These terms are only used to distinguish the description. For example, without departing from the scope of the embodiments of the present application, the first can also be called the second, and similarly, the second can also be called the first.

[0021] Depending on the context, the word “if” as used herein can be interpreted as “when” or “upon” or “in response to determining” or “in response to detecting”. Similarly, depending on the context, the phrase “if it is determined” or “if (a stated condition or event) is detected” can be interpreted as “when it is determined” or “in response to determining” or “when (a stated condition or event) is detected” or “in response to detecting (a stated condition or event)”.

[0022] It is also important to note that the terms "comprises", "comprising", or any other variations thereof, are intended to cover a non-exclusive inclusion, such that a process or method that comprises a list of elements does not include only those elements but can also include other elements not expressly listed or inherent to such process or method. An element proceeded by "comprises a" does not, without further constraints, exclude the presence of additional identical elements in the process or method.

[0023] In particular, it should be noted that symbols and / or numbers present in the description, if not marked in the description of the figures, are not figure references.

[0024] For the transient stability evaluation method in the related art, the method defects can be classified into the following three problems: (1) unable to depict the global boundary and manifold distribution law of the stability domain, resulting in over-conservative or inaccurate estimation of the system safety margin. (2) It is difficult to quantify the influence of new energy proportion and control parameters (such as virtual inertia, damping coefficient) on the size and shape of the stability domain. The existing eigenvalue analysis ignores the nonlinear characteristics, and the energy function method lacks parameter sensitivity evaluation capability. (3) Lack of comparative analysis of different transient control strategies, making it difficult to guide the controller setting and strategy optimization in the fault ride-through or high penetration rate scenario.

[0025] In view of this, the following provides a grid-connected system transient stability analysis method to solve the above technical problems.

[0026] The optional embodiments of the present application will be described in detail below with reference to the accompanying drawings.

[0027] The embodiments provided by the present application are a grid-connected system transient stability analysis method.

[0028] The following will be described in detail Figure 1 The embodiments of the present application will be described in detail.

[0029] Figure 1 A flowchart of a grid-connected system transient stability analysis method is shown, as shown in Figure 1 The following steps S1 to S3 are included.

[0030] S1: Establish a mathematical model of a heterogeneous new energy grid-connected system including network-constructing type new energy units, network-following type new energy units, synchronous generators and receiving-end power grids, and define the state variable set of the system.

[0031] The state variable set includes the power angle, angular velocity, virtual synchronous machine power angle, phase-locked loop phase and grid-connected point voltage amplitude of each unit.

[0032] S2: According to the system operating condition, the stable equilibrium point and unstable equilibrium point of the system are solved.

[0033] S3: Taking the unstable equilibrium point as the starting point, the numerical integration method or numerical continuation method is used for calculation to obtain the stable manifold starting from the unstable equilibrium point.

[0034] Among them, the stable manifold is used to depict the boundary of the transient stability domain of the system.

[0035] The transient stability analysis method of the grid-connected system provided in the application realizes global and accurate depiction of the transient stability domain by directly calculating the stable manifold starting from the unstable equilibrium point, overcomes the defects of the traditional small-signal linearization method that can only analyze local stability and the traditional energy function method that relies on artificial construction of functions and has limited applicability; at the same time, the method converts the abstract stability problem into an intuitive visual geometric figure in the phase plane, greatly enhances the understanding of system stability and the interpretability of the analysis results; in addition, the accurate stability domain boundary obtained lays a reliable foundation for subsequent quantitative researches such as parameter sensitivity analysis and control strategy comparison, fundamentally changes the analysis paradigm of the transient stability domain, and provides a more effective analysis tool for safe and stable operation of high-penetration new energy grid-connected systems.

[0036] In some embodiments, in S1, the network-constructed new energy unit adopts a virtual synchronous generator control strategy, and the mathematical model is represented based on the following formula: ; ; Among them, is the virtual inertia, is the mechanical input power, is the electromagnetic output power, is the power angle, is the actual angular velocity of the rotor, is the reference angular velocity, and · represents the first-order derivative of the corresponding parameter.

[0037] In some embodiments, in S1, the network-constructed new energy unit realizes grid connection through a phase-locked loop, and the mathematical model of the phase-locked loop is represented by the following formula: ; ; Among them, is the phase angle output by the phase-locked loop, is the angular velocity output by the phase-locked loop, is the grid phase angle, is the grid voltage amplitude, and are the proportional coefficient and the integral coefficient, respectively, and · represents the first-order derivative of the corresponding parameter.

[0038] As a feasible embodiment, the synchronous generator adopts a classical second-order model, considering electromagnetic transient process and mechanical dynamics characteristics; and the receiving end power grid adopts equivalent power source representation, and contains transmission line impedance, load and power collection line coupling relationship.

[0039] As a feasible embodiment, the state variable set of the system includes the power angle of each unit, the angular velocity of each unit, the power angle of the virtual synchronous machine, the phase of the phase-locked loop and the voltage amplitude of the grid connection point, and the initial power distribution scheme is determined according to the power flow calculation result and the set value of the control system.

[0040] Figure 2 A schematic diagram of a heterogeneous new energy grid-connected system model is shown.

[0041] As shown in Figure 2 , the connection relationship among the grid-forming new energy unit, the grid-following new energy unit, the synchronous generator and the receiving end power grid is included. Among them, the grid-forming new energy unit provides virtual inertia response and damping support through the virtual synchronous generator control, and the grid-following new energy unit mainly realizes power injection and dynamic response control through the phase-locked loop tracking grid reference phase.

[0042] In some embodiments, in S2, after solving the stable equilibrium point, the method further includes: first-order linearization of the system nonlinear equation at the stable equilibrium point to obtain a Jacobian matrix. The type of the equilibrium point is determined by calculating the eigenvalues of the Jacobian matrix.

[0043] In some specific embodiments, the type of the equilibrium point includes a stable node, a saddle point or a center point.

[0044] In some specific embodiments, the Jacobian matrix is obtained by taking the partial derivative of the system nonlinear equation at the stable equilibrium point, and is used to represent the linearized dynamic characteristics of the system near the equilibrium point. The equilibrium point type determination result is used as the initial condition or input basis for subsequent stable manifold calculation, thereby ensuring the accuracy and convergence of the stable manifold calculation.

[0045] In the above embodiments, by linearizing and eigenvalue analysis at the stable equilibrium point, the type of the system equilibrium point is accurately determined, which provides an important intermediate result and calculation basis for transient stability analysis.

[0046] In the embodiments of the present application, the unstable equilibrium point is taken as the starting point, and a numerical integration method or a numerical continuation method is used for calculation to obtain a stable manifold starting from the unstable equilibrium point. The envelope relationship formed by the stable manifold and the stable equilibrium point describes the transient stability domain boundary of the grid-connected system.

[0047] In some embodiments, in S3, after the stable manifold is calculated, the method further comprises: normalizing the stable manifold under different operating conditions, and analyzing the influence of system parameters on the transient stability domain by comparing the morphological differences of different stable manifolds.

[0048] The morphological differences include at least one of an area enclosed by the stable manifold in a phase plane, a stretching direction, and a convergence rate.

[0049] In the analysis, different combinations of system parameters are selected, and the corresponding stable manifolds are calculated. By comparing the morphological differences of different stable manifolds, the influence of system parameters on the transient stability domain is analyzed.

[0050] As a specific embodiment, the system parameters include: a power ratio of grid-forming units to grid-following units; a virtual inertia (J); a damping coefficient (D).

[0051] The morphological differences include at least one of an area enclosed by the stable manifold in a phase plane, a stretching direction, and a convergence rate.

[0052] For example, as the proportion of grid-forming units increases, the stable equilibrium point and the transient stability domain may move right as a whole along the power angle direction; increasing the virtual inertia can widen the attraction domain, and increasing the damping coefficient helps to improve the convergence rate.

[0053] In a specific implementation, the numerical integration method is implemented by using a fourth-order Runge-Kutta method, the integration step size is set to 1 ms, and the integration direction is from the unstable equilibrium point to the interior of the system. The normalization is realized by coordinate transformation, which translates the stable equilibrium point to the origin of the phase plane.

[0054] As a specific embodiment, the morphological differences are quantified by the following methods: Area: the area enclosed by the stable manifold is calculated by using a numerical integration method; Stretching direction: the angle between the principal axis of the stable manifold and the power angle coordinate axis is calculated; Convergence rate: the distance function of the trajectory points to the stable equilibrium point is fitted to obtain.

[0055] In the above embodiments, through stable manifold calculation and parameter influence analysis, accurate description of the transient stability domain and quantitative evaluation of the influence of system parameters are realized, providing a technical basis for system planning and operation.

[0056] In some embodiments, after S3, the method further comprises: S4: constructing a scalar function for quantitatively describing the scale characteristics of the stable manifold, and judging the scale size of the stable domain in different state variable directions by analyzing the sign and distribution of the rate of change of the scalar function over time.

[0057] Figure 3 A flow chart of another method for analyzing transient stability of a grid-connected system is shown in FIG. 1, which comprises the following steps S1-S4. Figure 3

[0058] S1: Establish a mathematical model of a heterogeneous new energy grid-connected system comprising a network-constructing new energy unit, a network-following new energy unit, a synchronous generator and a receiving-end power grid, and define a set of state variables of the system.

[0059] S2: Solve the stable equilibrium point and the unstable equilibrium point of the system according to the operating condition of the system.

[0060] S3: Take the unstable equilibrium point as the starting point, and perform calculation by using a numerical integration method or a numerical continuation method to obtain a stable manifold from the starting point.

[0061] S4: Construct a scalar function for quantitatively describing the scale characteristics of the stable manifold, and determine the scale size of the stable domain in different directions of state variables by analyzing the sign and distribution of the rate of change of the scalar function with time.

[0062] In some embodiments, the scalar function is a quadratic function, which is expressed by the following formula: ; wherein V is the scalar function, x1 and x2 are state variables, and k1 and k2 are weight coefficients.

[0063] In some embodiments, the scale size of the state variable in each coordinate direction (the dominant constraint direction) is determined by calculating the sign and distribution of dV / dt. In the analysis, the critical scale constants k1 and k2 are set, and the scale of the state variable corresponds to the bifurcation or elongation or shortening of the direction of the manifold. The sign of the derivative and the direction signal of the manifold can quantitatively characterize the scale of the stable boundary of different state variables (such as the power angle and the angular velocity). It is known that: wherein at the equilibrium point: and wherein .

[0064] On this basis, based on the generalized second-order swing equation, the energy function equation is the total energy = kinetic energy + potential energy, which is expressed as follows: ; At the Taylor expansion: ; On this basis, the following can be obtained: ; ​​​Further, comparing the defined quadratic form scalar function, a set of scalars can be obtained: ; Project (i.e. take the derivative) of the function Taylor expanded at : ; In the equation: Linearize around the equilibrium point: ; Where, .

[0065] Get: ; Substitute back to the projected function: .

[0066] For ease of understanding, the method described in the present application is described in detail below through specific examples. The present embodiment provides a specific implementation scheme of transient stability analysis of a heterogeneous new energy grid-connected system, and the system configuration is as follows.

[0067] In some embodiments, the heterogeneous new energy grid-connected system includes a synchronous generator, a grid-forming new energy unit, and a grid-following new energy unit. The total capacity of the system is 100 MW, and the grid-connected voltage level is 220 kV.

[0068] As a specific example, the parameters of each unit are configured as follows: synchronous generator: moment of inertia J = 60, damping coefficient D = 30; grid-forming new energy unit: moment of inertia J = 40, damping coefficient D = 25, droop coefficient = 0.05; grid-following new energy unit: moment of inertia J = 5, damping coefficient D = 10, current loop proportional gain K p = 2, current loop integral gain K i = 100; simulation settings: build the model of the heterogeneous new energy grid-connected system in the simulation environment, set the total simulation time to 8 s, and set the simulation step to 1 ms. The phase space analysis uses a two-dimensional phase space composed of power angle δ and angular velocity deviation Δω. On this basis, the following shows a complete implementation process.

[0069] Firstly, the system mathematical model is established according to the above parameter configuration, and the state variable set is defined to include the power angle, angular velocity, virtual synchronous machine power angle, phase-locked loop phase and grid point voltage amplitude of each unit. Secondly, the system operating condition is given, and the stable equilibrium point and unstable equilibrium point are solved. Linearization is carried out at the stable equilibrium point, the Jacobian matrix is obtained, and the eigenvalue is calculated to determine the type of equilibrium point. Then, taking the unstable equilibrium point as the starting point, the numerical integration method is used to calculate the stable manifold, and the integration step is kept at 1ms. The obtained manifold is normalized for easy comparison of different conditions. Finally, a quadratic scalar function is constructed, and the values of weight coefficients k1 and k2 are adjusted to fit the geometric shape of the stable manifold, and the scale characteristics in different state variable directions are analyzed. Further, the system stability is analyzed by Lyapunov function, the boundary of the attraction domain is calculated, and the stable operation range of the system is predicted. The phase trajectory and manifold distribution reveal the dynamic behavior characteristics of the system under disturbance, including the change law of stable equilibrium point, unstable equilibrium point and attraction domain.

[0070] The embodiment verifies the effectiveness and practicability of the method through specific parameter configuration and simulation setting, and provides a complete technical scheme and implementation example for transient stability analysis of the heterogeneous new energy grid-connected system. The method can accurately depict the stable domain boundary of the system, and provides an effective analysis means for safe and stable operation of the system.

[0071] The embodiment of the application provides a transient stability analysis method for a grid-connected system, which reveals the influence law of system parameters on the transient stability domain by analyzing the stable manifold shape of the system under different new energy ratios.

[0072] In some embodiments, the new energy ratio refers to the proportion of grid-forming new energy units and grid-following new energy units in the total capacity of the system, which affects the equivalent virtual inertia and equivalent damping coefficient of the system. Synchronous generators provide high inertia and damping, grid-forming new energy units provide medium virtual inertia and damping, and grid-following new energy units provide low inertia and damping. Therefore, a high new energy ratio usually leads to an accelerated overall dynamic response of the system and a decrease in transient stability.

[0073] Figure 4 Fig. 4 shows the change law of the attraction domain of the heterogeneous system under different new energy ratios.

[0074] As a specific embodiment, the stable manifold representing the boundary of the transient stability domain of the system can be calculated and its trajectory can be plotted in the phase plane. The phase plane usually takes power angle difference and angular velocity difference as coordinate axes. Figure 4 In the figure, the blue flow line represents the motion trajectory of the system state, the red circle represents the position of the stable equilibrium point (SEP), and the black cross represents the position of the unstable equilibrium point (UEP).

[0075] By comparing the shape of stable manifold under different operating conditions, the variation law of system transient stability can be analyzed, which is specifically manifested as: When the synchronous generator occupies a high proportion, that is, the new energy penetration is low, the system equivalent virtual inertia and damping coefficient are large. At this time, the stable manifold in the phase plane presents a "coarse and stable" shape, the streamline is dense and tends to the stable equilibrium point. The system has a wide attraction domain, and the separation trajectory from the unstable equilibrium point forms a closed annular manifold, making the system state trajectory not easy to lose stability, indicating that the system has strong resistance to large disturbances.

[0076] When the new energy occupies a medium level (for example, 20% to 40%), and the configuration of grid-forming and grid-following units is relatively balanced, the manifold structure can still maintain the annular shape, but the attraction domain range is reduced, and the streamline becomes sparse, indicating that the system dynamic process is accelerated. If the grid-forming unit occupies a relatively high proportion, its virtual synchronous control (VSG) characteristics can help enhance system stability; if the grid-following unit dominates, the system faces an increased risk of subsynchronous oscillation. The continuity of the equilibrium manifold is usually good under this operating condition, but if the capacity of the grid-following unit exceeds a certain ratio (for example, 1:5) of the grid-forming unit, the manifold may appear discontinuous, resulting in a lack of stable equilibrium point in the system and causing synchronization instability.

[0077] When the new energy occupies a high or super-high range (for example, more than 60%), whether dominated by grid-forming or grid-following units, the phase trajectory shows rapidly diverging streamlines, and the system attraction domain is significantly reduced. The separation trajectory is easily across the stable boundary to cause instability, and the "stable ring" in the manifold is significantly "flattened" on the angular velocity difference (Δω) axis, reducing the stability margin of the system and making it more susceptible to disturbances.

[0078] As a feasible embodiment, when the grid-forming new energy unit proportion increases, it approximately retains the "inertia-damping" external characteristics of the synchronous generator, and the manifold structure in the phase plane is still regular, and the oscillation frequency and decay speed of the system depend on the set virtual inertia and damping coefficient. When the grid-following new energy unit proportion increases, the system equivalent virtual inertia and damping coefficient decrease, the separation trajectory shrinks towards the stable equilibrium point, and the system state is more likely to cross the saddle point into the instability region (such as the continuous rotation or divergence region).

[0079] Through the above analysis of the shape of the stable manifold, the following intuitive law can be obtained: increasing the system's equivalent virtual inertia or equivalent damping coefficient helps to expand the transient stability domain (attraction domain), reduce the power angle swing amplitude, and increase the system stability margin. This law provides a theoretical basis and analysis method for optimizing system parameter configuration and improving the transient stability of high-proportion new energy grid-connected systems.

[0080] Figure 5A schematic diagram of the variation of the attraction domain of the heterogeneous system under different new energy ratios is shown.

[0081] As Figure 5 shown, the influence of the change of the equivalent virtual inertia of the system on the transient stability domain can be quantitatively analyzed. The reduction of the equivalent virtual inertia is usually associated with the increase of the new energy ratio. The analysis results show that, with the reduction of the equivalent virtual inertia, the range of the attraction domain corresponding to the transient stability domain of the system presents a shrinking trend. Among them, the size of the attraction domain directly reflects the ability of the system to resist disturbances. The smaller the range of the attraction domain, the lower the disturbance limit that the system can withstand, and the more likely the system loses transient stability when disturbed.

[0082] As a feasible embodiment, when the capacity ratio of the synchronous generator or the grid-connected new energy unit in the system is high, the overall equivalent virtual inertia and the damping coefficient of the system are large, the critical disturbance power is high, the boundary of the transient stability domain is expanded outward, and the stability of the system is enhanced. On the contrary, when the ratio of the grid-following new energy unit increases, the equivalent virtual inertia and the damping coefficient of the system decrease, resulting in the inward contraction of the boundary of the transient stability domain, the range of the initial operating state that the system can maintain stable becomes smaller, and the anti-disturbance ability is weakened.

[0083] For example, when subjected to the same amplitude of disturbance, the system containing a higher proportion of grid-connected new energy units is more difficult to lose stability because it has a larger transient stability domain and stronger damping characteristics, and its transient response of power angle and angular velocity is more likely to converge to a stable equilibrium point.

[0084] In summary, the method provided by the application can clearly reveal the internal law that the type and ratio of new energy units affect the equivalent parameters of the system and then act on the transient stability domain, providing an important analysis basis for system planning and operation control.

[0085] In some embodiments, the control parameters of the grid-connected new energy unit, especially the virtual inertia and the damping coefficient, have a significant influence on the stability manifold form and transient stability characteristics of the system. By adjusting these parameters, the equivalent virtual inertia and the equivalent damping coefficient of the system can be changed, and the transient response of the system can be optimized.

[0086] Figure 6 A schematic diagram of the influence of the inertia parameter of the grid-connected converter on the system manifold law is shown.

[0087] Figure 7 A schematic diagram of the influence of the damping parameter of the grid-connected converter on the system manifold law is shown.

[0088] As Figure 6 and Figure 7As shown, increasing the virtual inertia of grid-forming new energy units will directly increase the equivalent virtual inertia of the system. The intuitive performance of this change in the phase plane of the system is that the rotation speed of the system state trajectory (streamline) around the stable equilibrium point slows down, i.e., the inherent oscillation frequency of the system decreases. At the same time, the entire stable flow field is "stretched" in the angular velocity difference axis direction, and the phase trajectory form becomes more "fat" or more open. This shows that the system's ability to suppress frequency changes and withstand power shocks is enhanced, but its transient response speed may slow down. The main role of virtual inertia is to enhance the system's anti-disturbance ability and expand the transient stability domain, but it should be noted that simply increasing the virtual inertia cannot effectively suppress the sustained oscillation of the system, and it needs to be adjusted in cooperation with the damping coefficient.

[0089] As a specific embodiment, increasing the damping coefficient of grid-forming new energy units will increase the equivalent damping coefficient and damping ratio of the system. The performance of this change in the phase plane is that the system state trajectory converges to the stable equilibrium point at a faster speed. Specifically, when the system is in a low damping state (i.e., under-damped, damping ratio ζ < 1), the trajectory presents a clear spiral line form of convergence, indicating that the system has an oscillation process; when the system is in a high damping state (i.e., over-damped, damping ratio ζ > 1), the trajectory converges to the equilibrium point in a direct and relatively slow manner, and this process has no oscillation but the dynamic response is slow; and when the damping is properly configured, a fast and smooth convergence process can be achieved, thereby optimizing the transient stability of the system.

[0090] In summary, virtual inertia and damping coefficient are key control parameters that affect the boundary form and dynamic performance of the system's transient stability domain. By coordinating the virtual inertia and damping coefficient of grid-forming new energy units, the geometric characteristics of the stable flow field can be effectively shaped, thereby improving the power angle stability and dynamic quality of the system under large disturbances.

[0091] Figure 8 An energy change diagram corresponding to the scalar function expression is shown.

[0092] In some embodiments, the control strategy changes the rules of the stable flow field and the transient stability characteristics of the system by affecting the energy function of the system. The effect of the control strategy can be evaluated by constructing and analyzing the scalar function and its derivative.

[0093] As a feasible embodiment, the scalar function adopts a quadratic form, which appears as a series of contour lines in the phase plane. As shown in FIG. 3, the scalar function is expressed as Figure 8As shown, the stable equilibrium point of the system corresponds to the minimum point of the scalar function, and the contours around it form a closed "bowl-shaped" concentric ring structure. The density of the contours reflects the distribution of the curvature of the scalar function: increasing the virtual inertia will reduce the curvature of the scalar function in the angular velocity direction, which is manifested as the same contour being widened in the vertical direction, which is mapped to the system characteristics as the frequency change being more "dull"; adjusting the weight coefficient k1 (which is related to the system parameters) affects the curvature in the power angle direction, and an increase in the value of k1 increases the curvature in this direction, making the "bowl wall" steeper.

[0094] As a specific embodiment, by plotting the scalar function V' of its time derivative, the energy change trend and evolution law of the stable manifold of the system can be further analyzed.

[0095] Figure 9 A schematic diagram of manifold law analysis under different control strategies is shown.

[0096] As Figure 9 As shown on the left, the manifold line is the system vector field (derivative) change trend, which shows how the energy is dissipated. The minimum point is at (0, 0), and the contours of the left figure form a "bowl-shaped" closed concentric ring around this point, and the density of the contours reflects the curvature: J is large in the direction of the power angle, and the "bowl wall" is steeper, and the "bowl wall" is steeper. In the direction, the frequency change is more "dull" (the same contour is widened more in the vertical direction), and the frequency change is more "dull". Large in the direction of the power angle, the "bowl wall" is steeper, and the "bowl wall" is steeper. Figure 9 As shown on the right, the derivative function is plotted as an isosurface, and the black contour where the derivative function is zero is specially marked: the area where the derivative function is less than zero (the "negative" part corresponding to the main color block in the right figure) indicates that the system energy is dissipating, and the trajectory will approach the energy valley; the black line where the derivative function is zero is the dividing line of "instantaneous energy invariance"; the area where the derivative function is greater than zero indicates that the energy is increasing (such as electrical-mechanical power imbalance, too small damping, or far from the linear region).

[0097] In the above embodiment, by analyzing the geometric shape and distribution of the scalar function and its derivative function, it can be intuitively revealed how different control parameters or strategies change the energy surface and dissipation characteristics of the system, and ultimately affect the shape of the stable manifold and the transient stability of the system. This provides a clear theoretical basis and visual analysis method for the design and optimization of control strategies.

[0098] In some embodiments, the present application realizes the quantitative characterization of the system stable manifold scale by constructing a specific scalar function and adjusting its weight coefficients. This method uses the quadratic energy function to approximate the transient stability region boundary of the system under different weight combinations, and can diagnose the dominant influence of the power angle direction and the frequency direction on the stability region boundary shape. This method can not only locally reflect the characteristics of the energy function near the equilibrium point, but also capture the geometric shape of the attractor boundary in the global range by adjusting the parameters.

[0099] As a feasible embodiment, the scalar function adopts a quadratic form. By drawing the equipotential lines of the scalar function, the zero contour lines of its derivative function, and the actual phase trajectory streamlines of the system, and performing superposition analysis, the fitting effect of the scalar function on the real stability region boundary can be evaluated. The equipotential lines of the scalar function exhibit elliptical geometric characteristics, and which coordinate axis is "flattened" is determined by the weight coefficients k1 and k2.

[0100] Figure 10 A popular scale characterization diagram is shown.

[0101] As Figure 10 shown, by adjusting the values of k1 and k2, the geometric shape of the stable manifold is fitted, and the dominance of the power angle direction and the angular velocity direction on the stability region boundary is diagnosed; increasing k1 will flatten the ellipse in the power angle direction, indicating that the function pays more attention to the power angle deviation; increasing k2 will flatten the ellipse in the angular velocity difference direction, indicating that the function pays more attention to the frequency deviation. The black zero contour line of the derivative function divides the phase plane into energy dissipation and energy increase areas, which is the key boundary for stability discrimination. The blue streamline shows the actual phase trajectory of the system. By comparing the proximity of the zero contour line and the real trajectory, the effectiveness of the scalar function fitting can be verified.

[0102] As Figure 10 shown, columns 1 and 2 are the superposition of scalar function equipotential lines, derivative zero contour lines, and system streamlines. The equipotential lines show the elliptical geometric characteristics of the scalar function, and the black lines divide the plane into areas with derivative less than 0 and greater than 0, which is the key to stability discrimination. The blue streamline shows the phase trajectory, which is convenient for verifying the proximity of the zero contour line and the real attractor boundary. Which axis the ellipse is "flattened" along is determined by k1 and k2. k1 is large → "flattened" in the power angle direction (more attention to angle deviation); k2 is large → "flattened" in the angular velocity difference direction (more attention to frequency deviation). In columns 3 and 4, the larger the blue area (derivative less than 0), the closer the scalar function under this parameter combination to the real Lyapunov function; the position where the red area (derivative greater than 0) appears is often when the power angle and the frequency deviate at the same time The fifth and sixth columns show the relative contribution of the power angle direction and the frequency direction in the derivative (dominant scale), and if the color block is biased to 1, the power angle scale is dominant, and if it is biased to 0, the frequency scale is dominant. By comparing (1, 1), (4, 1), (1, 4), and other different weight combinations, it can be concluded that increasing k1 makes the black zero line more vertical and better fits the separation trajectory, but the red area is easy to expand; increasing k2 makes the blue area significantly increase, which is conducive to stability discrimination, but the zero line tends to be horizontal, and the geometric consistency of the fitting boundary decreases; (4, 1) or (8, 2) and other combinations achieve a good compromise between shape and blue area ratio, and are better fitting choices.

[0103] In the above embodiments, the scalar function scaling method can be used to adjust the weight coefficient to approximate the real stable manifold, and the influence degree of different state variables on the stability domain boundary can be quantitatively analyzed, thereby providing an effective analysis tool for evaluating and optimizing system transient stability. Among them, Figure 10 The curve can be used to fit Figure 9 the boundary of the attraction domain determined by the energy function. Figure 10 The scalar function method shown is an extension of Figure 9 the scaling of the energy function. It is consistent with the physical energy function locally, and realizes fitting by adjusting the weight globally, which retains the theoretical interpretability and enhances the engineering practicability. Compared with Figure 9 the strict dissipation criterion, Figure 10 a more flexible tool is provided for characterizing the dominant scale and fitting characteristics of the stability domain under different control strategies.

[0104] In some embodiments, the transient stability analysis method provided in the present application has higher accuracy in global stability domain characterization. The method calculates and draws the stable manifold of the system from the unstable equilibrium point as the starting point, thereby directly and intuitively obtaining the accurate boundary of the transient stability domain of the system. Compared with the traditional analysis method based on fixed energy function, the method avoids conservative assumptions and can more truly reflect the stability boundary shape of the nonlinear system.

[0105] As a feasible embodiment, the method performs more comprehensive in parameter sensitivity analysis. By systematically changing the proportion of new energy output, the configuration proportion of grid-forming and grid-following units, and key control parameters (such as virtual inertia and damping coefficient), and comparing and analyzing the differences in the shape of the stable manifold, the area surrounded, and the convergence speed of the trajectory under different parameter combinations, the quantitative evaluation of the influence of parameters can be realized. The analysis results provide valuable reference for the controller parameter tuning in system operation, and the new energy capacity configuration and selection in the planning stage.

[0106] As a specific embodiment, the method can be used to guide the optimization design of control strategy. By superimposing and comparing the actual phase trajectory of the system under different transient control strategies with the calculated stable manifold (i.e. the boundary of the stable domain), the risk of crossing the boundary under disturbance and the recovery efficiency after fault removal can be clearly evaluated. This intuitive comparison helps to select the control strategy that is optimal in maintaining stability and improving dynamic performance from multiple candidate strategies.

[0107] In summary, the transient stability analysis method for heterogeneous new energy grid-connected system based on manifold characteristics proposed in the present application can comprehensively reveal the manifold distribution law and scale characteristics of the system in the phase plane, significantly improving the accuracy and intuitiveness of transient stability analysis in high-proportion new energy penetration scenarios. This method not only deepens the understanding of the stability mechanism of the system, but also provides a solid theoretical basis and effective analysis tool for ensuring the safe and stable operation of the power system, optimizing control strategies and system planning.

[0108] Although the operations are described in a particular, sequential order, this should not be understood as requiring that the operations be performed in that order, or in serial, or that all operations be performed, to achieve desirable results. In certain circumstances, multitasking and parallel processing can be advantageous.

[0109] Any of the steps, operations, or procedures described herein can be performed or implemented using one or more hardware or software modules, alone or in combination with other devices. In one embodiment, a software module is implemented using a computer program product comprising a computer-readable medium containing computer program code, which can be executed by a computer processor for performing any or all of the steps, operations, or procedures described.

[0110] The foregoing description of implementation of the present application has been presented for the purposes of illustration and description. It is not intended to be exhaustive or to limit the application to the precise form disclosed, and various modifications and variations can be possible in light of the above teachings or can be acquired from practice of the application. One skilled in the art will appreciate the choice and arrangement of the embodiments to best suit the particular use contemplated, and the scope of the application includes various embodiments and various modifications as are within the purview of persons skilled in the art from the teachings of the present application.

[0111] As to the apparatus in the above-described embodiments, the specific manner in which each module performs operations has been described in detail in the embodiments related to the method, and will not be described in detail here.

[0112] It can be further understood that, unless otherwise specified, "connection" includes direct connection between the two without other components, and also includes indirect connection between the two with other elements.

[0113] It will be further appreciated that embodiments of the application are capable of being implemented and practiced in a variety of ways and that the application is not limited in its application to any one of the exemplary embodiments described below. Further, it should be understood that implementations of the application can include any one or more of the advantageous features set forth above or any combination of any such features. Moreover, it will be appreciated that some embodiments can be capable of outperforming other embodiments.

[0114] Other embodiments of the application will be apparent to those skilled in the art from consideration of the specification and practice of the application disclosed herein. It is intended that the specification and examples be considered as exemplary only, with the true scope and spirit of the application being indicated by the following claims.

[0115] It is to be understood that the application is not limited to the precise details of construction and the exact arrangements of the components described above and illustrated in the drawings, and that various modifications and changes can be made without departing from the scope of the application. The scope of the application is limited only by the claims that follow.

[0116] The above embodiments are only used to illustrate the technical solutions of the present application, not limit the present application; even though the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A method of transient stability analysis of a grid-connected system, characterized by, The method comprises: S1: establishing a mathematical model of a heterogeneous new energy grid-connected system comprising a network-constructing new energy unit, a network-following new energy unit, a synchronous generator and a receiving-end power grid, and defining a state variable set of the system; S2: solving stable equilibrium points and unstable equilibrium points of the system according to system operating conditions; S3: taking the unstable equilibrium point as a starting point, and performing calculation by using a numerical integration method or a numerical continuation method to obtain a stable manifold from the starting point; The stable manifold is used to depict a transient stability domain boundary of the system.

2. The method of claim 1, wherein, In the S1, the network-constructing new energy unit adopts a virtual synchronous generator control strategy, and the mathematical model is represented based on the following formula: ； ； wherein is the virtual inertia, is the mechanical input power, is the electromagnetic output power, is the power angle, is the actual angular velocity of the rotor, is the reference angular velocity, and · denotes the first derivative of the corresponding parameter.

3. The method of claim 1, wherein, In the S1, the network-following new energy unit realizes grid connection through a phase-locked loop, and the mathematical model of the phase-locked loop is represented by the following formula: ； ； wherein is the phase angle output by the phase-locked loop, is the angular velocity output by the phase-locked loop, is the grid phase angle, is the grid point voltage amplitude, and are the proportional and integral coefficients, respectively, and · denotes the first derivative of the corresponding parameter.

4. The method of claim 1, wherein, In the S2, after the stable equilibrium points are solved, the method further comprises: performing first-order linearization on a system nonlinear equation at the stable equilibrium point to obtain a Jacobian matrix; determining the type of the equilibrium point by calculating eigenvalues of the Jacobian matrix; The type of the equilibrium point comprises a stable node, a saddle point or a center point.

5. The method of claim 1, wherein, In the S3, after the stable manifold is calculated, the method further comprises: performing normalization processing on the stable manifolds under different operating conditions, and analyzing the influence of system parameters on the transient stability domain by comparing the morphological differences of different stable manifolds; The morphological differences comprise at least one of an area, a stretching direction and a convergence rate of the stable manifold in a phase plane; The system parameters comprise a power ratio of the network-constructing unit and the network-following unit, a virtual inertia and a damping coefficient.

6. The method according to any one of claims 1 to 5, characterized in that, After the S3, the method further comprises: S4: constructing a scalar function for quantitatively describing the scale characteristics of the stable manifold, and determining the scale size of the stable domain in different state variable directions by analyzing the sign and distribution of the rate of change of the scalar function over time.

7. The method of claim 6, wherein, The scalar function is a quadratic function, and is expressed by the following formula: ； wherein is a scalar function, x1and x2are state variables, and k1and k2are weight coefficients.

8. The method of claim 7, wherein, The method further comprises: By adjusting the values of k1 and k2, the geometric shape of the stable manifold is fitted, and the dominance of the power angle direction and the angular velocity direction on the stability domain boundary is diagnosed; Increasing k1 causes the equipotential lines of the scalar function to be flattened in the power angle direction, and increasing k2 causes the equipotential lines to be flattened in the angular velocity direction.

9. The method of claim 8, wherein, The value combination of the weight coefficients k1 and k2 is (4, 1) or (8, 2).

10. The method of claim 1, wherein, The state variable set comprises the power angle, the angular velocity, the virtual synchronous machine power angle, the phase-locked loop phase and the grid-connected point voltage amplitude of each unit.

Images