A stability domain boundary calculation and evaluation method for new energy power systems

By combining manifold theory and the Monte Carlo method, the stability domain boundary of the new energy power system is accurately characterized, which solves the conservative problem in the transient synchronous stability assessment of the new energy power system and achieves efficient stability assessment.

CN122092201APending Publication Date: 2026-05-26NORTHEAST DIANLI UNIVERSITY +1
View PDF 0 Cites 0 Cited by

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
NORTHEAST DIANLI UNIVERSITY
Filing Date
2025-12-29
Publication Date
2026-05-26

Smart Images

  • Figure CN122092201A_ABST
    Figure CN122092201A_ABST
Patent Text Reader

Abstract

This invention relates to a method for calculating and evaluating the stability domain boundary of a new energy power system, addressing the problem of significant conservatism in evaluating the transient synchronization stability of grid-connected and grid-connected converters in current new energy power systems. The method includes: calculating the equilibrium point of the target nonlinear system and determining the local stable manifold of the stable equilibrium point; plotting the global stable manifold of the target nonlinear system using manifold theory and accurately observing the dynamic behavior of the target nonlinear system in the global state space; further extending the proposed transient stability domain boundary calculation method to a high-dimensional phase space; and quantitatively evaluating the transient stability margin of the target nonlinear system under different faults using the Monte Carlo method. The calculation method of this invention accurately characterizes the stability domain boundary of the target nonlinear system, providing important reference value for evaluating the transient synchronization stability of new energy power systems.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of power system stability domains, and more particularly to a method for calculating and evaluating the stability domain boundary of new energy power systems. Background Technology

[0002] With the ongoing low-carbon energy transition, renewable energy and auxiliary resources are being connected to the grid via power electronic converters and synchronized through power control. However, this also brings new control interaction and instability issues. As the proportion of renewable energy continues to grow, the transient synchronization stability of power systems based on power electronics becomes crucial.

[0003] In practical engineering applications, the direct method based on energy function theory uses a constant energy surface to estimate the stability domain boundary of the system after a fault. The transient stability of the system is determined by the position and distance between the fault trajectory and the stability domain boundary. Therefore, accurate characterization of the stability domain boundary is the core of the direct method theory and a prerequisite for its effective application in the transient synchronization stability analysis of power systems. Thus, accurate characterization of the stability domain boundary has significant reference value for evaluating the transient synchronization stability of new energy power systems.

[0004] However, the stability domain boundary estimated by isoenergy surfaces is not the original, exact stability domain boundary of the target system model, and the accuracy of the stability domain boundary is directly related to the constructed energy function. Existing techniques cannot avoid the conservatism of the stability domain boundary resulting from this approximation, thus affecting the accurate evaluation of the stability domain boundary.

[0005] Application content

[0006] Based on the above analysis, this invention aims to provide a method for calculating and evaluating the stability domain boundary of a new energy power system, so as to accurately characterize the stability domain boundary of the target nonlinear system, and to quantitatively evaluate the transient synchronization stability of the target nonlinear system using the Monte Carlo algorithm. This addresses the problem of excessive conservatism in evaluating the transient synchronization stability of grid-connected and grid-connected converters in current new energy power systems.

[0007] To address the aforementioned problems, this invention provides a method for calculating and evaluating the stability domain boundary of new energy power systems, comprising:

[0008] Step S1: Calculate the equilibrium point of the target nonlinear system and determine the local stable manifold of the stable equilibrium point;

[0009] Step S2: Use manifold theory to plot the global stable manifold of the target nonlinear system and accurately observe the dynamic behavior of the target nonlinear system in the global state space;

[0010] Step S3: Further extend the proposed transient stability domain boundary calculation method to high-dimensional phase space;

[0011] Step S4: Quantitatively evaluate the transient stability margin of the target nonlinear system under different faults using the Monte Carlo method.

[0012] Preferably, step S1 includes:

[0013] Step S101: Use the fsolve equation root solver in PyCharm software to select the initial vector, and through continuous iterative updates, approach the equilibrium point of the target equation.

[0014] Step S102: Using the two equilibrium points obtained from the solution, generate the Jacobian matrix of the target nonlinear system and calculate the eigenvalues ​​of the equilibrium points;

[0015] Step S103: Select the eigenvector corresponding to the equilibrium point with positive real eigenvalues, and take a center point with a radius of ε that is small enough as an initial value along the direction of the eigenvector to perform forward and reverse integration on the target nonlinear system to generate a locally unstable manifold and a stable manifold.

[0016] The target nonlinear system is represented as:

[0017]

[0018] Where f(x) is a polynomial of variable x, if x * Satisfy f(x) * If x = 0, then it is the equilibrium point of the system; if the eigenvalues ​​of the Jacobian matrix at the equilibrium point all have negative real parts, then x * It is a stable equilibrium point, denoted as x. s If the real part of an eigenvalue is not zero, then it is classified as a hyperbolic equilibrium point.

[0019] For the hyperbolic equilibrium point x * Its stable manifold W s (x * and unstable manifold W u (x * The definition is as follows:

[0020]

[0021]

[0022] For the stable equilibrium point x s There exists a region in the state space. As time increases, each trajectory x(t) tends towards x. s This region is called x. sThe stability region is denoted as A(x). s ), defined as:

[0023]

[0024] x represents s The stability region A(x) s The boundary of t is the boundary of the stable region, which is the interface separating stable and unstable trajectories; as t→∞, all trajectories approach the stable equilibrium point; this means that the stable region has the same topological properties as the stable manifold.

[0025] Preferably, step S103 specifically includes:

[0026] Step S10301: Set the system state variables to δ = ω = 0, and solve for the hyperbolic equilibrium point x of the system. u ;

[0027] Step S10302: Calculate the second-order objective nonlinear system at the hyperbolic equilibrium point x. u The Jacobian matrix at point A is used to determine its eigenvectors a1 and a2.

[0028] Step S10303: Establish the initial coordinate point for time integration as x u ±α1 and x u ±α2;

[0029] Step S10304: Perform forward and backward time integration for each initial point.

[0030] Preferably, step S2 includes:

[0031] Step S201: Given the unstable equilibrium point x calculated in step S1. u For locally stable manifolds, utilize the unstable equilibrium point x u Solving for the eigenvectors of the Jacobian matrix, and defining the stable eigenspace spanned by the eigenvectors as E(x) u );

[0032] Step S202: Incorporate positive eigenvalues ​​into the stable eigenspace E(x) u In subspace E(x) u Define a circle with radius ε (sufficiently small) and center x in the following code. u The circle, discretized into individual

[0033]

[0034] Step S203, with Starting from the initial point, the trajectory is solved using inverse time integration over the initial discrete points in the set. The iteration duration is determined by either a fixed time increment Δt or a trajectory arc length increment Δs. The iteration stops when the length reaches a set value L, resulting in the first-generation trajectory denoted as {T}. 1,1 ,T 1,2 ,T 1,3 ,T 1,4 ,...,T 1,k Let the endpoint of each track be denoted as {p}. 1,1 ,p 1,2 ,p 1,3 ,p 1,4 ,...,p 1,k};

[0035] Step S204: Detect the set of endpoints of the trajectory {p} 1,1 ,p 1,2 ,p 1,3 ,p 1,4 ,...,p 1,k The density of the midpoints, if the distance between the two lines is greater than D. max Insert a new initial point between the initial points of the two orbits; if the distance between the two orbits is less than D... min If so, delete the initial point of one trajectory line; readjust the set of initial points;

[0036] Step S205: Utilize the corrected set of endpoints from step S204.

[0037] {p 1,1 ',p 1,2 ',p 1,3 ',p 1,4 ',...,p 1,k '} is used as the initial value for a new round of iterations, and the trajectory is solved using inverse time integration;

[0038] Step S206: A space compression function is introduced; this function maps the values ​​of the state variables from [-∞, +∞] to [-1, 1]; the selected space compression function is:

[0039]

[0040] The resulting interface is the stability domain boundary of the target nonlinear system.

[0041] Preferably, step S206 includes: mapping the vector field of the target nonlinear system onto a compact manifold, and defining the projected distance from a point in the state space to the origin o as:

[0042]

[0043] Where ||x|| represents the vector norm in n. If the trajectory is close to the origin, the vector norm will be projected onto the space where y(x) < 1; if the trajectory is far from the origin, the vector norm will be projected onto the space where y(x) → 1. Let n = 3 in the n-dimensional state space. For any trajectory N(x,y,z) in the phase space, the distance from the trajectory to the origin O(0,0,0) is:

[0044]

[0045] The distance Γ varies along the trajectory of N; the projected distance ψ of the trajectory is defined as:

[0046]

[0047] A new closed sphere can be represented as:

[0048]

[0049] Preferably, the transient stability domain boundary calculation method proposed in step S3 can be further extended to high-dimensional phase space; by using coordinate transformation, the stability domain boundary is compressed onto the Poincaré sphere, and in the high-dimensional state space, the length of any compression vector is used by the spatial compression function, and the dynamic characteristics of the target nonlinear system are observed by the changes in the size and shape of the stability domain volume or area.

[0050] Preferably, under different fault conditions, the Monte Carlo method is used to quantitatively evaluate the transient stability margin of the target nonlinear system, and step S4 includes:

[0051] Step S401: Define the irregular shape and specify the region;

[0052] Step S402: Randomly generate N points within the specified area;

[0053] Step S403: Check whether each point is within the irregular shape;

[0054] Step S404: Calculate the area of ​​the irregular shape based on the ratio of the number of points falling within the irregular shape to the total number of points in the specified area: S = M / N;

[0055] Step S405: Output the estimated area of ​​the irregular shape.

[0056] Beneficial effects

[0057] The present invention provides a method for calculating and evaluating the stability domain boundary of a new energy power system, comprising: calculating the equilibrium point of the target nonlinear system and determining the local stable manifold of the stable equilibrium point; drawing the global stable manifold of the target nonlinear system using manifold theory and accurately observing the dynamic behavior of the target nonlinear system in the global state space; further extending the proposed transient stability domain boundary calculation method to a high-dimensional phase space; and quantitatively evaluating the transient stability margin of the target nonlinear system under different faults using the Monte Carlo method. Attached Figure Description

[0058] To more clearly illustrate the technical solutions in this invention or the prior art, the accompanying drawings involved in the embodiments or the prior art are briefly described below. Obviously, these drawings illustrate several embodiments of the present invention, and those skilled in the art can derive other possible drawings based on these drawings without creative effort. The purpose of the drawings is limited to illustrating specific embodiments and does not limit the scope of the present invention.

[0059] Figure 1 This is a flowchart of the method for calculating and evaluating the stability domain boundary of a grid-type converter system in Embodiment 1 of the present invention;

[0060] Figure 2 This is a detailed model block diagram of the grid-type converter system in Embodiment 1 of the present invention;

[0061] Figure 3 This is a control structure diagram of the active and reactive power loops of the grid-type converter system in Embodiment 1 of the present invention;

[0062] Figure 4 This is a schematic diagram of the stability domain boundary of the grid-type converter system in Embodiment 1 of the present invention;

[0063] Figure 5 This is the two-dimensional boundary diagram of the global stability domain of the grid-type converter system in Case 1, which does not consider reactive power coupling control, in Embodiment 2 of the present invention.

[0064] Figure 6 This is the three-dimensional boundary diagram of the global stability domain of the grid-type converter system in Case 1, as shown in Embodiment 2 of the present invention;

[0065] Figure 7 This is a two-dimensional boundary phase trajectory diagram of the global stability domain of the grid-type converter system in Case 1, without considering reactive power coupling control, in Embodiment 2 of the present invention.

[0066] Figure 8 is a three-dimensional boundary phase trajectory diagram of the global stability domain of the grid-type converter system in Case 1, considering reactive power coupling control in Embodiment 2 of the present invention.

[0067] Figure 9 This is the stability margin diagram of the grid converter system evaluated by the Monte Carlo method as damping changes in Case 1, in Embodiment 2 of the present invention. Detailed Implementation

[0068] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. 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 skilled in the art without creative effort are within the protection scope of the present invention.

[0069] Example 1

[0070] This embodiment discloses a method for calculating and evaluating the stability domain boundary of a new energy power system. Example 1 uses a grid-type converter system as an example. Figures 1 to 6 As shown, it includes the following steps:

[0071] Step S1: Calculate the equilibrium point of the grid-type converter system and determine the local phase trajectory of the stable equilibrium point;

[0072] Step S1 includes:

[0073] Step S101: The grid-connected converter, as the interface for renewable energy generation to connect to the grid, generally uses virtual synchronous control to establish a connection with the grid. For example... Figure 2 As shown, Example 1 uses a grid-type converter system as an example, where... Figure 3 The control structures for the active and reactive power loops are given. The differential equations for the grid-type converter system are:

[0074]

[0075] In the formula, P is the actual power per unit value; P0 is the corresponding reference electrical quantity; J and D represent the moment of inertia and damping coefficient of the virtual synchronous generator control, respectively; ωN is the reference virtual speed nominal value; Δω* represents the nominal value of the difference between the actual speed and the reference virtual speed. δ represents the virtual power angle (called angle δ), which is usually considered to be the angle at which the d-axis exceeds the infinite grid U. Considering power coupling control, the output power of the grid-type converter can be expressed in the following form:

[0076]

[0077]

[0078] This allows us to derive the third-order dynamic equations of the grid converter system. Initial vectors are selected using the equation solver in PyCharm software. Through continuous iterative updates, the equilibrium point of the grid converter equations is gradually approximated.

[0079] Step S102: Using the two equilibrium points obtained from the solution, generate the Jacobian matrix of the target nonlinear system and calculate the eigenvalues ​​of the equilibrium points.

[0080] Step S103: Select the eigenvector corresponding to the equilibrium point with positive real eigenvalues, take a center point with radius ε along the direction of the eigenvector as the initial value, and perform forward and reverse integration on the target nonlinear system to generate locally unstable manifolds and stable manifolds.

[0081] A grid-type converter system is represented as:

[0082]

[0083] Where f(x) is a polynomial of variable x, if x * Satisfy f(x) * If x = 0, then it is the equilibrium point of the system; if the eigenvalues ​​of the Jacobian matrix at the equilibrium point all have negative real parts, then x * It is a stable equilibrium point, denoted as x. s If the real part of an eigenvalue is not zero, it is classified as a hyperbolic equilibrium point.

[0084] For the hyperbolic equilibrium point x * Its stable manifold W s (x * and unstable manifold W u (x * The definition is as follows:

[0085]

[0086] For the stable equilibrium point x s There exists a region in the state space. As time increases, each trajectory x(t) tends towards x. s This region is called x. s The stability region is denoted as A(x). s ), defined as:

[0087]

[0088] x represents s The stability region A(x) sThe boundary of a stable region is the boundary of a stable region, serving as the interface separating stable and unstable trajectories. As t→∞, all trajectories approach a stable equilibrium point; this implies that stable regions and stable manifolds possess the same topological properties. The relationship between stable regions, stable region boundaries, and stable manifolds is as follows: Figure 4 As shown.

[0089] Step S103 includes step S10301: setting the system state variable to δ=ω=0, and solving for the hyperbolic equilibrium point x of the system. u ;

[0090] Step S10302: Calculate the second-order objective nonlinear system at the hyperbolic equilibrium point x. u The Jacobian matrix at point A is used to determine its eigenvectors a1 and a2.

[0091] Step S10303: Establish the initial coordinate point for time integration as x u ±α1 and x u ±α2;

[0092] Step S10304: Perform forward and backward time integration for each initial point.

[0093] Figure 4 The solid black line represents the stable manifold at the unstable equilibrium point, while the dashed black line represents the unstable manifold at the unstable equilibrium point. The fact that one of the unstable manifolds converges to a stable equilibrium point indicates that the interior of the stable manifold is a stable region. If the operating point after a fault is within the stable region, it will eventually reach a stable equilibrium point. Conversely, if it is outside the stable region, the converter will become unstable. Calculations show that the stable region is not only a local region but also an infinite region, so it should be observed from a global perspective.

[0094] Step S2: Use manifold theory to draw the global stable manifold of the grid converter system and accurately observe the dynamic behavior of the grid converter system in the global state space.

[0095] Step S201: Given the unstable equilibrium point x calculated in step S1. u For locally stable manifolds, utilize the unstable equilibrium point x u Solving for the eigenvectors of the Jacobian matrix, and defining the stable eigenspace spanned by the eigenvectors as E(x) u );

[0096] Step S202: Incorporate positive eigenvalues ​​into the stable eigenspace E(x) u In subspace E(x) u Define a circle with radius ε (sufficiently small) and center x in the following code. u The circle, discretized into individual Step S203, with Starting from the initial point, the trajectory is solved using inverse time integration over the initial discrete points in the set. The iteration duration is determined by either a fixed time increment Δt or a trajectory arc length increment Δs. The iteration stops when the length reaches a set value L, resulting in the first-generation trajectory denoted as {T}. 1,1 ,T 1,2 ,T 1,3 ,T 1,4 ,...,T 1,k Let the endpoint of each track be denoted as {p}. 1,1 ,p 1,2 ,p 1,3 ,p 1,4 ,...,p 1,k};

[0097] Step S204: Detect the set of endpoints of the trajectory {p} 1,1 ,p 1,2 ,p 1,3 ,p 1,4 ,...,p 1,k The density of the midpoints, if the distance between the two lines is greater than D. max Insert a new initial point between the initial points of the two orbits; if the distance between the two orbits is less than D... min Then delete the initial point of one trajectory; readjust the initial point set; step S205, use the corrected solution trajectory end point set {p} from step S204. 1,1 ',p 1,2 ',p 1,3 ',p 1,4 ',...,p 1,k '} is used as the initial value for a new round of iterations, and the trajectory is solved using inverse time integration;

[0098] Step S206: A space compression function is introduced; this function maps the values ​​of the state variables from [-∞, +∞] to [-1, 1]; the selected space compression function is:

[0099]

[0100] The resulting interface is the stability domain boundary of the target nonlinear system.

[0101] Step S206 includes: mapping the vector field of the target nonlinear system onto a compact manifold, and defining the projected distance from a point in the state space to the origin o as:

[0102]

[0103] Where ||x|| represents the vector norm in n. If the trajectory is close to the origin, the vector norm will be projected onto the space where y(x) < 1; if the trajectory is far from the origin, the vector norm will be projected onto the space where y(x) → 1. Let n = 3 in the n-dimensional state space. For any trajectory N(x,y,z) in the phase space, the distance from the trajectory to the origin O(0,0,0) is:

[0104]

[0105] The distance Γ varies along the trajectory of N; the projected distance ψ of the trajectory is defined as:

[0106]

[0107] A new closed sphere can be represented as:

[0108]

[0109] Step S3: The proposed transient stability domain boundary calculation method is further extended to high-dimensional phase space. The proposed transient stability domain boundary calculation method can be further extended to high-dimensional phase space. Coordinate transformation is used to compress the stability domain boundary onto a Poincaré sphere, thus more intuitively presenting the dynamic characteristics of the high-dimensional nonlinear system. In the high-dimensional state space, the length of any compression vector can be used with a spatial compression function, allowing observation of the dynamic characteristics of the target nonlinear system through changes in the volume or area and shape of the stability domain.

[0110] Step S4: Quantitatively evaluate the transient stability margin of the target nonlinear system under different faults using the Monte Carlo method. Under different fault conditions, the Monte Carlo method is used to quantitatively evaluate the transient stability margin of the target nonlinear system. Generally, the larger the area and volume of the stability region, the better the transient synchronization stability of the system. Step S4 includes:

[0111] Step S401: Define the irregular shape and specify the region;

[0112] Step S402: Randomly generate N points within the specified area;

[0113] Step S403: Check whether each point is within the irregular shape;

[0114] Step S404: Calculate the area of ​​the irregular shape based on the ratio of the number of points falling within the irregular shape to the total number of points in the specified area: S = M / N;

[0115] Step S405: Output the estimated area of ​​the irregular shape.

[0116] Example 2

[0117] use Figure 2The transient synchronous stability assessment method for a grid-type converter system was verified. Nonlinear simulations were performed using DIgSILENT and PowerFactory. Figure 3 The control structures for the active and reactive power loops are presented. Table 1 shows some parameters for Case Study 1.

[0118] Table 1 System parameter settings for grid-type converter systems

[0119]

[0120] First, in the control strategy of grid-type converters without considering reactive power coupling control, Figure 5 The global stability domain boundary of the grid-connected converter in Case Study 1 was plotted. In the control strategy for the grid-connected converter considering reactive power coupling control... Figure 6 The global stability domain boundary of the grid-connected converter was plotted. Secondly, the scheme simulates faults in the grid-connected converter system using voltage dips as an example, introducing grid voltage dip events as a significant disturbance. Figure 7 As shown in Figures 8(a) and 8(b), the transient synchronization stability of a grid-connected converter system without considering reactive power coupling control was studied. When the system reached 1.0 s, the grid voltage dropped from 1 p.u. to 0.2 p.u. The study found that the phase trajectory corresponding to the limit clearing time coincided with the drawn stability domain boundary. As shown in Figures 8(a) and 8(b), the transient synchronization stability of a grid-connected converter system considering reactive power coupling control was studied. When the system reached 1.5 s, the grid voltage dropped from 1 p.u. to 0.2 p.u. Faults were cleared at 1.65 s and 1.72 s, respectively. As shown in Figure 8(a), when the fault was cleared at 1.65 s, the clearing point was within the stability domain, indicating transient stability of the system. As shown in Figure 8(b), when the fault was cleared at 1.72 s, the clearing point was outside the stability domain, indicating transient instability of the system.

[0121] The scheme in this embodiment demonstrates that the calculated stability region boundary is highly accurate. The phase trajectory corresponding to the limit cut-off time closely matches the plotted stability region boundary, thus proving the correctness of the experimental results and the stability region calculation method. Finally, the stability region area of ​​the grid converter system was calculated using the Monte Carlo method. Figure 9As shown, the area of ​​the stability region increases with increasing damping, indicating enhanced transient synchronous stability of the grid-type converter system. After D = 12.3, the ASR (Active Stability Ratio) undergoes a sudden change, indicating a shift in the stability boundary type. This change can be explained as follows: the system involves more than one stable equilibrium point and one unstable equilibrium point. When the stable region of the initial stable equilibrium point expands to a certain value with increasing damping, it encompasses the previous unstable equilibrium point. At this point, the stability boundary undergoes a sudden change, and the system transforms into a multi-oscillator stable system across the entire state space, no longer exhibiting the previous unstable form. This Monte Carlo method can accurately estimate the global stability region, thereby assessing the magnitude of the stability margin through changes in area or volume.

[0122] In summary, the stability domain boundary calculation and evaluation method proposed in this invention for new energy power systems solves the problem of excessive conservatism in evaluating the transient synchronization stability of grid-connected and grid-connected converters in new energy power systems. First, the proposed stability domain boundary calculation method describes the structural characteristics of the stability domain boundary at the equilibrium point of the nonlinear system by characterizing the global manifold. This method directly presents the dynamic characteristics of the high-dimensional space and simplifies the complexity of stability domain boundary calculation. Second, the proposed Monte Carlo method accurately evaluates the stability margin of the target nonlinear system under fault-free conditions. Finally, examples verify the effectiveness and reliability of the proposed method in practical applications. In the future, this method will provide a new engineering application approach for accurately evaluating the transient synchronization stability of power electronic power systems.

[0123] Those skilled in the art should understand that the above embodiments are merely illustrative of the content of this disclosure and do not limit its scope. The system capacity, voltage, line parameters, etc., shown may vary depending on the specific circumstances of the power electronic grid-connected generator set and its grid connection. Based on this disclosure, those skilled in the art can make other changes or adjustments, and these changes still fall within the scope of this disclosure.

Claims

1. A method for calculating and evaluating the stability domain boundary of a new energy power system, characterized in that, include: Step S1: Calculate the equilibrium point of the target nonlinear system and determine the local stable manifold of the stable equilibrium point; Step S2: Use manifold theory to plot the global stable manifold of the target nonlinear system and accurately observe the dynamic behavior of the target nonlinear system in the global state space; Step S3: Further extend the proposed transient stability domain boundary calculation method to high-dimensional phase space; Step S4: Quantitatively evaluate the transient stability margin of the target nonlinear system under different faults using the Monte Carlo method.

2. The method for calculating and evaluating the stability domain boundary of a new energy power system according to claim 1, characterized in that, Step S1 includes: Step S101: Use the fsolve equation root solver in PyCharm software to select the initial vector, and through continuous iterative updates, approach the equilibrium point of the target equation. Step S102: Using the two equilibrium points obtained from the solution, generate the Jacobian matrix of the target nonlinear system and calculate the eigenvalues ​​of the equilibrium points; Step S103: Select the eigenvector corresponding to the equilibrium point with positive real eigenvalues, and take a center point with a radius of ε that is small enough as an initial value along the direction of the eigenvector to perform forward and reverse integration on the target nonlinear system to generate a locally unstable manifold and a stable manifold. The target nonlinear system is represented as: Where f(x) is a polynomial of variable x, if x * Satisfy f(x) * If x = 0, then it is the equilibrium point of the system; if the eigenvalues ​​of the Jacobian matrix at the equilibrium point all have negative real parts, then x * It is a stable equilibrium point, denoted as x. s If the real part of an eigenvalue is not zero, then it is classified as a hyperbolic equilibrium point. For the hyperbolic equilibrium point x * Its stable manifold W s (x * and unstable manifold W u (x * The definition is as follows: For the stable equilibrium point x s There exists a region in the state space. As time increases, each trajectory x(t) tends towards x. s This region is called x. s The stability region is denoted as A(x). s ), defined as: x represents s The stability region A(x) s The boundary of t is the boundary of the stable region, which is the interface separating stable and unstable trajectories; as t→∞, all trajectories approach the stable equilibrium point; this means that the stable region has the same topological properties as the stable manifold.

3. The method for calculating and evaluating the stability domain boundary of a new energy power system according to claim 2, characterized in that, Step S103 specifically includes: Step S10301: Set the system state variables to δ = ω = 0, and solve for the hyperbolic equilibrium point x of the system. u ; Step S10302: Calculate the second-order objective nonlinear system at the hyperbolic equilibrium point x. u The Jacobian matrix at point A is used to determine its eigenvectors a1 and a2. Step S10303: Establish the initial coordinate point for time integration as x u ±α1 and x u ±α2; Step S10304: Perform forward and backward time integration for each initial point.

4. The method for calculating and evaluating the stability domain boundary of a new energy power system according to claim 1, characterized in that, Step S2 includes: Step S201: Given the unstable equilibrium point x calculated in step S1. u For locally stable manifolds, utilize the unstable equilibrium point x u Solving for the eigenvectors of the Jacobian matrix, and defining the stable eigenspace spanned by the eigenvectors as E(x) u ); Step S202: Incorporate positive eigenvalues ​​into the stable eigenspace E(x) u In subspace E(x) u Define a circle with radius ε (sufficiently small) and center x in the following code. u The circle, discretized into individual Step S203, with Starting from the initial point, the trajectory is solved using inverse time integration for the initial discrete points in the set. The iteration time is determined by either a fixed time increment Δt or a trajectory arc length increment Δs. The iteration stops when the length reaches a set value L, resulting in the first-generation trajectory denoted as {T}. 1,1 ,T 1,2 ,T 1,3 ,T 1,4 ,...,T 1,k Let the endpoint of each track be denoted as {p}. 1,1 ,p 1,2 ,p 1,3 ,p 1,4 ,...,p 1,k }; Step S204: Detect the set of endpoints of the trajectory {p} 1,1 ,p 1,2 ,p 1,3 ,p 1,4 ,...,p 1,k The density of the midpoints, if the distance between the two lines is greater than D. max Insert a new initial point between the initial points of the two orbits; if the distance between the two orbits is less than D... min If so, delete the initial point of one trajectory line; readjust the set of initial points; Step S205: Utilize the corrected solution trajectory endpoint set {p} from step S204. 1,1 ',p 1,2 ',p 1,3 ',p 1,4 ',...,p 1,k '} is used as the initial value for a new round of iterations, and the trajectory is solved using inverse time integration; Step S206: A space compression function is introduced; this function maps the values ​​of the state variables from [-∞, +∞] to [-1, 1]; the selected space compression function is: The resulting interface is the stability domain boundary of the target nonlinear system.

5. The method for calculating and evaluating the stability domain boundary of a new energy power system according to claim 4, characterized in that, Step S206 includes: mapping the vector field of the target nonlinear system onto a compact manifold, and defining the projected distance from a point in the state space to the origin o as: Where ||x|| represents the vector norm in n. If the trajectory is close to the origin, the vector norm will be projected onto the space where y(x) < 1; if the trajectory is far from the origin, the vector norm will be projected onto the space where y(x) → 1. Let n = 3 in the n-dimensional state space. For any trajectory N(x,y,z) in the phase space, the distance from the trajectory to the origin O(0,0,0) is: The distance Γ varies along the trajectory of N; the projected distance ψ of the trajectory is defined as: A new closed sphere can be represented as:

6. The method for calculating and evaluating the stability domain boundary of a new energy power system according to claim 1, characterized in that, The transient stability domain boundary calculation method proposed in step S3 can be further extended to high-dimensional phase space; By using coordinate transformation to compress the boundary of the stable domain onto a Poincaré sphere, and in a high-dimensional state space, by using the length of an arbitrary compression vector through a spatial compression function, the dynamic characteristics of the target nonlinear system can be observed through changes in the volume or area and shape of the stable domain.

7. The method for calculating and evaluating the stability domain boundary of a new energy power system according to claim 1, characterized in that, The Monte Carlo method is used to quantitatively evaluate the transient stability margin of the target nonlinear system under different fault conditions. Step S4 includes: Step S401: Define the irregular shape and specify the region; Step S402: Randomly generate N points within the specified area; Step S403: Check whether each point is within the irregular shape; Step S404: Calculate the area of ​​the irregular shape based on the ratio of the number of points falling within the irregular shape to the total number of points in the specified area: S = M / N; Step S405: Output the estimated area of ​​the irregular shape.