Full-pure embedded power system power flow method considering engineering constraint

Abstract

The invention provides a full-pure embedded power system load flow method considering engineering constraints. The method comprises the following steps: acquiring power grid data and establishing a load flow calculation model; performing model conversion based on the load flow calculation model, establishing a pseudo-gradient system, and calculating a stable equilibrium manifold of the pseudo-gradient system; all stable equilibrium manifolds are searched based on a layer-by-layer search mode, and power flow solutions meeting specific conditions are screened out; calculating the distance between the power flow solution and a boundary value for ensuring stable operation of the power system, and judging whether the distance meets the required reserved margin or not; and if not, planning a related scheduling strategy for increasing power generation and reducing load. According to the method, the corresponding polynomial model is established for the power flow problem of the power system considering engineering constraints, and the power flow problem of the power system with constraints is efficiently solved, so that state analysis of the power system is facilitated; through the uniformly distributed search directions, the possibility of repeatedly entering the same convergence domain due to concentrated search from similar directions can be remarkably reduced, and the search efficiency is improved.

Description

Technical Field

[0001] The invention relates to the technical field of power systems, and in particular to a fully embedded power system power flow method considering engineering constraints. Background Art

[0002] As the core of the power industry, the power system and power flow problems are considered to be one of the most important issues in the design and operation of power systems. Such problems involve bus voltage, current and power. The power flow equation describes the relationship between the complex bus voltage and the active and reactive power injection of the node in a steady state. It can be expressed as a polynomial system in rectangular coordinates, and for a given set of parameters, the system often has multiple complex solutions, but only a few real solutions that meet specific conditions are of application value in engineering practice. For example, when carrying out power system power flow calculations, usually only the power flow solutions that are within a small range of the rated voltage can meet the normal operation requirements of the system.

[0003] For calculating the real number solution of a polynomial system with constraints, the existing methods mainly adopt the following three schemes:

[0004] 1. Calculate the resultant of a polynomial system. Although this solution is supported by solid mathematical theory, the task of calculating the resultant of a polynomial system is an NP-hard problem. Existing technologies mainly use this solution to deal with the real number solution problem of a small-scale polynomial system containing only a few variables, which is difficult to be applied to actual engineering calculations.

[0005] 2. In or near the feasible domain determined by the constraints, select or randomly generate several points as the initial points for iterative calculation. Using traditional iterative methods such as Newton's method, it is very likely to search for a real solution to the polynomial system (with constraints). However, considering that the convergence domain of traditional iterative methods is usually irregular and its boundary has a fractal structure, it is generally difficult to predict the location of the final searched solution based on the initial point, and there is also a lack of solid mathematical theory to guarantee the integrity of the solution set obtained.

[0006] 3. First, find all the complex solutions in the complex space, and then select the real solutions that meet the requirements. This method has two obvious defects: First, for large-scale problems, the number of all complex solutions is very large, and the homotopy continuation method is often difficult to calculate the upper bound of the number of solutions through the mixed-volume method, making it difficult to construct the initial system and its initial solution set; second, even if the corresponding initial system and its initial solution set can be constructed, on the one hand, the existing polynomial system calculation tools mainly use the prediction-correction framework to implement the homotopy continuation method, which usually has problems such as path jumping, resulting in missed solutions, and it is difficult to implement continuation in real space. For example, when encountering a saddle node bifurcation point in the continuation path, the Jacobian matrix is ​​singular, and it is difficult to predict the continuation direction, and it is impossible to complete the homotopy continuation numerically; on the other hand, if all the complex solutions are calculated by complex space continuation, and then the real solutions that meet specific conditions are selected, since the proportion of real solutions that meet specific conditions is often very low, this will cause most of the computing resources to be consumed in the task of calculating non-real solutions, and the computing efficiency will be greatly reduced. Summary of the invention

[0007] In view of the deficiencies of the prior art, the present invention provides a fully pure embedded power system flow method taking into account engineering constraints. The present invention can efficiently solve the power system flow problem with constraints so as to be used for power system state analysis.

[0008] The technical solution of the present invention is: a fully embedded power system power flow method considering engineering constraints, comprising the following steps:

[0009] S1), obtain power grid data and establish a power flow calculation model;

[0010] S2), performing model conversion based on the power flow calculation model and establishing a pseudo-gradient system, and calculating the stable equilibrium manifold of the pseudo-gradient system;

[0011] S3) Find all stable equilibrium manifolds based on a layer-by-layer search method and select the power flow solutions that meet specific conditions;

[0012] S4) Calculate the distance between the power flow solution and the boundary value for ensuring stable operation of the power system, and determine whether the distance meets the required reserved margin; if not, plan relevant scheduling strategies for increasing power generation and reducing load.

[0013] Preferably, in step S1), firstly, a power flow calculation equation is established at the PQ node i:

[0014]

[0015] In the formula, Represents the voltage V at node i iThe real and imaginary parts of G; i, k represent the node numbers respectively; n is the total number of nodes; G ik ,B ik are the conductance and susceptance between nodes i and k respectively; P i ,Q i are the real power and imaginary power of node i, respectively, where node 1 is a relaxed node whose voltage is fixed and known.

[0016] Preferably, in step S1), a power flow calculation equation is then established at PV node i:

[0017]

[0018] In the formula, Represents the voltage V at node i i The real and imaginary parts of G; i, k represent the node numbers respectively; n is the total number of nodes; G ik ,B ik are the conductance and susceptance between nodes i and k respectively; P i is the real power of node i; V i,sp Represents the voltage amplitude given by node i.

[0019] As a preferred embodiment, in step S1), after establishing the power flow equation, the voltage V i Add constraints, namely:

[0020]

[0021] Preferably, in step S2), the expression of the established pseudo-gradient system is:

[0022]

[0023] Where H(x) is a vector function determined by the original equations and constraints, DH(x) represents the Jacobian matrix of H(x), and the superscript T represents the transpose of the matrix.

[0024] Preferably, in step S2), the stable equilibrium manifold of the pseudo-gradient system is calculated by using the logarithmic method, and the solution function x k Substitute the power series expansion of (t) at any t = t0 into equation (4) and calculate the power series coefficient a by comparing the coefficients of the same power of (t-t0). kq ;

[0025] The expression of the power series expansion is:

[0026] x k (t)=∑ q≥0 a kq (t-t0) q ,1≤k≤m+d.

[0027] Preferably, in step S2), after a stable equilibrium manifold of the pseudo-gradient system is calculated, an escape point is calculated to jump out of the convergence domain of the current stable equilibrium manifold to find another adjacent stable equilibrium manifold.

[0028] The beneficial effects of the present invention are:

[0029] 1. The present invention establishes a corresponding polynomial model for the power system flow problem considering engineering constraints, and efficiently solves the power system flow problem with constraints, so as to be used for power system state analysis;

[0030] 2. The present invention uses a model conversion method to convert the task of solving the real number solution of the polynomial system into the problem of solving the stable equilibrium manifold of the pseudo-gradient dynamic system, and directly searches for all real number solutions of the constrained polynomial system in the real number domain; it overcomes the defect that the existing screening method that needs to calculate all complex number solutions and then screen out the real number solutions is difficult to construct the initial system and its initial solution set when dealing with ultra-large-scale problems and cannot calculate the real number solutions;

[0031] 3. The present invention uses a vector composed of points evenly distributed on the unit sphere and the origin as the search direction, which can significantly reduce the possibility of repeatedly entering the same convergence domain due to concentrated searches from similar directions, thereby improving search efficiency;

[0032] 4. The present invention adopts the logarithmic method when calculating the coefficients of the power series expansion, which can reduce the time complexity of the combination number of the commonly used convolution method to the linear time complexity;

[0033] 5. The present invention adopts piecewise rational functions for approximation, and can construct a high-precision approximate numerical solution function, rather than just obtaining a finite number of discrete points on the pseudo-gradient system solution curve. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 It is a framework structure diagram of the present invention;

[0035] Figure 2 It is a schematic diagram of the process of the present invention;

[0036] Figure 3 A schematic diagram of the process of searching layer by layer for all stable equilibrium manifolds of the present invention;

[0037] Figure 4 A schematic diagram of the process of finding a stable equilibrium manifold of adjacent layers of the present invention;

[0038] Figure 5 A schematic diagram of the process of calculating the escape point of the present invention;

[0039] Figure 6A schematic diagram of the process of calculating a single stable equilibrium manifold according to the present invention. DETAILED DESCRIPTION

[0040] The specific implementation of the present invention will be further described below in conjunction with the accompanying drawings:

[0041] like Figure 1 and 2 As shown, this embodiment provides a fully embedded power system power flow method considering engineering constraints, including the following steps:

[0042] S1), obtain power grid data and establish a power flow calculation model;

[0043] S2), performing model conversion based on the power flow calculation model and establishing a pseudo-gradient system, and calculating the stable equilibrium manifold of the pseudo-gradient system;

[0044] S3) Find all stable equilibrium manifolds based on a layer-by-layer search method and select the power flow solutions that meet specific conditions;

[0045] S4) Calculate the distance between the power flow solution and the boundary value for ensuring stable operation of the power system, and determine whether the distance meets the required reserved margin; if not, plan relevant scheduling strategies for increasing power generation and reducing load.

[0046] As a preferred embodiment of this invention, in step S1), firstly, a power flow calculation equation is established at the PQ node i:

[0047]

[0048]

[0049] In the formula, Represents the voltage V at node i i The real and imaginary parts of G; i, k represent the node numbers respectively; n is the total number of nodes; G ik ,B ik are the conductance and susceptance between nodes i and k respectively; P i ,Q i are the real power and imaginary power of node i, respectively, where node 1 is a relaxed node whose voltage is fixed and known.

[0050] As a preferred embodiment of the present invention, in step S1), a power flow calculation equation is then established at PV node i:

[0051]

[0052] In the formula, Represents the voltage V at node i i The real and imaginary parts of G; i, k represent the node numbers respectively; n is the total number of nodes; Gik ,B ik are the conductance and susceptance between nodes i and k respectively; P i is the real power of node i; V i,sp Represents the voltage amplitude given by node i.

[0053] As a preferred embodiment of this invention, in step S1), after establishing the power flow equation, the voltage V i Add constraints, namely:

[0054]

[0055] Then, equations (1) and (2) with constraint (3) are unified as follows:

[0056] F(z)=(F1(z),…,F m (z)) T =0, stT(z)≤0;

[0057] Where z=(z1,…,z m ) T represents the unknown quantity, m is the number of unknown quantities, satisfying m=2n-2;

[0058] 1≤l≤m represents equations (1) and (2), represents the field of real numbers, represents a real polynomial, l is the subscript; T(z)=(T1(z),…,T d (z)) T It represents the constraint condition, 1≤l≤d represents inequality (3), and d=2n-2 is the number of inequality constraints.

[0059] As a preferred embodiment of this embodiment, in step S1), the inequality constraint is converted into a polynomial equality constraint by adding a slack variable, and the expression of the polynomial equality constraint is:

[0060]

[0061] In the formula, y=(y1,…,y d ) T ;y l is the slack variable, is the augmented variable; is the transformed equality constraint;

[0062] Then it is transformed into a computational augmented polynomial system The real number solution of m+d (x) T=0, where x is the unknown quantity to be solved. 1≤k≤m+d;

[0063] Since the number of unknowns in the augmented polynomial system H(x)=0 is the same as the number of equations, only isolated real solutions are considered.

[0064] As a preferred embodiment of the present invention, in step S2), the expression of the established pseudo-gradient system is:

[0065]

[0066] Where DH(x) represents the Jacobian matrix of H(x), and the superscript T represents the transpose of the matrix.

[0067] As a preferred embodiment of this invention, in step S2), the stable equilibrium manifold of the pseudo-gradient system is calculated by using the logarithmic method, and the solution function x is k Substitute the power series expansion of (t) at any t = t0 into equation (4) and calculate the power series coefficient a by comparing the coefficients of the same power of (t-t0). kq ;

[0068] The expression of the power series expansion is:

[0069] x k (t)=∑ q≥0 a kq (t-t0) q ,1≤k≤m+d;

[0070] By constructing a rational approximate function based on the existing power series information to expand the convergence domain, t = t0 + Δt = t0 + Δt * Substitute the rational approximation value and its derivative value calculated at into equation (4); compare whether the difference between the left and right sides is greater than the preset allowable threshold; find Δt so that the difference between the left and right sides of equation (4) is less than the preset allowable threshold; update t0←t0+Δt until the condition that satisfies Q is obtained. H (x) = 0; where t represents the time variable of the pseudo-gradient system, t0 represents the t value at any point; Δt represents the effective interval length, Δt * is the initial interval length for each expansion.

[0071] As a preferred embodiment of this embodiment, in step S2), after a stable equilibrium manifold of the pseudo-gradient system is calculated, an escape point is calculated to jump out of the convergence domain of the current stable equilibrium manifold to find another adjacent stable equilibrium manifold.

[0072] As a preferred embodiment of this invention, in step S2), if Figure 6As shown, the logarithmic method is used to calculate the stable equilibrium manifold of the pseudo-gradient system, which specifically includes the following steps:

[0073] S211), initialization, setting the order q of the power series partial sum * , the allowable threshold ε EM , the length of the initial interval of each expansion Δt * , reduction ratio Initial point x0=(x1(0),…,x m+d (0)), initial point t0 = 0;

[0074] S212), parameterizing by t, expressing the solution function of the pseudo-gradient system as x(t), wherein the solution function x(t) of the pseudo-gradient system satisfies:

[0075]

[0076] S213), assuming that the solution function x(t) of the pseudo-gradient system is holomorphic near t=t0, the series expansion of x(t) at t0 is expressed as follows:

[0077] x k (t)=Σ q≥0 a kq (t-t0) q ,1≤k≤m+d,

[0078] Among them, a kq is the power series coefficient;

[0079] Assume Q H The series expansion of (x(t)) at t0 is expressed as follows:

[0080] Q H,k (x(t))=∑ q≥0 Q H,kq (t-t0) q ,1≤k≤m+d,

[0081] Among them, Q H,k (x(t)) is Q H The kth element of (x(t)), Q H,kq Q H,k The power series coefficients of (x(t));

[0082] S214), Substitute the power series into equation (5) and compare the two sides of the equation (t-t0) q ,q≥0, solve the power series coefficient a in turn kq ;

[0083] S215) Compare both sides of the equation (t-t0) 0 The coefficient of ak0 (k=1,…,m+d) we have:

[0084] a k0 =x k (t0),1≤k≤m+d;

[0085] S216) Compare both sides of the equation (t-t0) q , the coefficient of q≥1, we get:

[0086] qa kq =Q H,k(q-1) ;

[0087] The specific steps of using the logarithmic method to solve the power series coefficients are as follows: For each monomial Satisfies the following three expressions:

[0088]

[0089]

[0090]

[0091] in and are the coefficients of the power series expansion of logarithmic variables, logarithmic monomials, and monomials, respectively. k The coefficients of the first q+1 powers of (t) are used to solve the monomial The steps for the power series expansion of the q+1 power coefficient are as follows:

[0092] From (6), we can get

[0093]

[0094] From (7), we can get

[0095]

[0096] From (8), we can get

[0097]

[0098] Where 1≤k≤m+d; due to Q H,k (x(t)) is a monomial is a linear combination of H,k(q+1) Also for The linear combination of H,k(q+1) Then solve a in sequence (12) kq ,Q H,kq :

[0099]

[0100] S217), assuming that the corresponding coefficient a is calculated kq ,q≤q * , using a truncated power series:

[0101] 1≤k≤m+d,

[0102] And calculate the corresponding rational approximation and the corresponding derivative, namely:

[0103]

[0104] Where L and N are the orders of the numerator and denominator of the rational approximation, respectively, and satisfy L+N=q * ; ′ represents the derivative, δ kq are the coefficients of the numerator and denominator of the rational approximation, and satisfy δ k0 =1;

[0105] S218), let the initial interval length Δt = Δt * >0; Check whether the maximum error of equation (5) at t∈[t0,t0+Δt] is less than the preset allowable error threshold ∈ EM ;

[0106] If it is satisfied, the effective interval is determined to be [t0, t0+Δt]; otherwise, the interval length is reduced to Until satisfied; update t0←t0+Δt;

[0107] S219), repeat steps S213) to S218) until the energy function is satisfied The value at t0 is less than 0.1 at the initial point, and then the local solver is used to calculate Q H The solution of (x)=0.

[0108] As a preferred embodiment of this invention, in step S2), if Figure 5 As shown, the escape point is calculated to find another adjacent stable equilibrium manifold, which specifically includes the following steps:

[0109] S221), initialize the initial stable equilibrium point X s , search path direction d, search step Δt>0, search stop maximum value F Lim >0;

[0110] S222), from x s Starting from the search path direction d, a small perturbation point is used as the initial iteration point Assume the number of iterations i = 0;

[0111] S223), calculation

[0112] S224), respectively calculated in The inner product of the derivative of the search path at the point and the vector field of the pseudo-gradient system, that is, d and -DH(x) T The inner product of H(x) is denoted by

[0113] S225), update i←i+1;

[0114] S226), repeat steps S223) to S225) until one of the following conditions is met;

[0115] 1) Returns successful solution and

[0116] 2) Function The value of or the value of i is greater than the given F Lim , returns the solution failed and

[0117] As a preferred embodiment of this invention, in step S3), if Figure 3 As shown in the figure, all stable equilibrium manifolds are found by layer-by-layer searching, and the corresponding power flow solutions are screened out for the obtained stable equilibrium manifolds, including the following steps:

[0118] S31), initialize the initial point x0, search step Δt>0, search stop maximum value F Lim >0, number of search paths N * >0;

[0119] S32), starting from the initial point x0, a stable equilibrium manifold of the pseudo-gradient system calculated in step S2) is obtained

[0120] S33), set the stable equilibrium manifold set to The set of stable equilibrium manifolds at level 0 is The escape point set is k = 0;

[0121] S34), set the stable equilibrium manifold set of the k+1th layer to be

[0122] S35) Every stable equilibrium manifold in Find the stable equilibrium manifolds of all its neighboring layers;

[0123] S36) Judgment Is it equal to If yes, go to step S37); otherwise, k←k+1, go to step S34);

[0124] S37) for V SEM Every stable equilibrium manifold x in s , determine whether it is a solution to the polynomial system H(x)=0.

[0125] As preferred in this embodiment, Figure 4 As shown, in step S35), the stable equilibrium manifold of the adjacent layers is found, which specifically includes the following steps:

[0126] S351), according to the uniform distribution of the unit sphere, construct N * Line search path 1≤i≤N * ;Suppose i=1;

[0127] S352), along Search path The corresponding escape point is found by using the step of calculating the escape point; if the escape point is successfully found, the process proceeds to step S353); otherwise, the process proceeds to step S356);

[0128] S353), define the escape point successfully found as And determine Is it V ex If yes, go to step S356); otherwise update Go to step S354);

[0129] S354) First, solve the pseudo-gradient system to obtain a stable equilibrium manifold And judge Is the solution successful? If successful, proceed to step S355); otherwise, proceed to step S356;

[0130] S355) Judgment Is it V SEM If not, update

[0131] S356), i=i+1;

[0132] S357) Repeat steps S352)-S356) until i>N * .

[0133] The above embodiments and descriptions are only for illustrating the principles and best embodiments of the present invention. Without departing from the spirit and scope of the present invention, the present invention may be subject to various changes and improvements, all of which fall within the scope of the present invention to be protected.

Claims

1. A fully embedded power system power flow method considering engineering constraints, characterized in that: The following steps are involved: S1), obtain power grid data and establish a power flow calculation model; S2), performing model conversion based on the power flow calculation model and establishing a pseudo-gradient system, and calculating the stable equilibrium manifold of the pseudo-gradient system; S3) Find all stable equilibrium manifolds based on a layer-by-layer search method and select the power flow solutions that meet specific conditions; S4) Calculate the distance between the power flow solution and the boundary value for ensuring stable operation of the power system, and determine whether the distance meets the required reserved margin; if not, plan relevant scheduling strategies for increasing power generation and reducing load.

2. The method of fully embedded power system power flow considering engineering constraints according to claim 1 is characterized by: In step S1), firstly, the power flow calculation equation is established at the PQ node i: In the formula, Represents the voltage V at node i i The real and imaginary parts of G; i, k represent the node numbers respectively; n is the total number of nodes; G ik ,B ik are the conductance and susceptance between nodes i and k respectively; P i ,Q i are the real power and imaginary power of node i respectively.

3. A fully pure embedded power system power flow method considering engineering constraints according to claim 2, characterized in that: In step S1), the power flow calculation equation is then established at PV node i: In the formula, Represents the voltage V at node i i The real and imaginary parts of G; i, k represent the node numbers respectively; n is the total number of nodes; G ik ,B ik are the conductance and susceptance between nodes i and k respectively; P i is the real power of node i; V i,sp Represents the voltage amplitude given by node i.

4. A fully pure embedded power system power flow method considering engineering constraints according to claim 3, characterized in that: In step S1), after establishing the power flow equation, the voltage V i Add constraints, namely:

5. The method of fully embedded power system power flow considering engineering constraints according to claim 1 is characterized by: In step S2), the expression of the established pseudo-gradient system is: Where H(x) is a vector function determined by the original equations and constraints, DH(x) represents the Jacobian matrix of H(x), and the superscript T represents the transpose of the matrix.

6. A fully pure embedded power system power flow method considering engineering constraints according to claim 5, characterized in that: In step S2), the logarithmic method is used to calculate the stable equilibrium manifold of the pseudo-gradient system, and the solution function x k Substitute the power series expansion of (t) at any t = t0 into equation (4) and calculate the power series coefficient a by comparing the coefficients of the same power of (t-t0). kq ; The expression of the power series expansion is: x k (t)=∑ q≥0 a kq (t-t0) q ,1≤k≤m+d。 7. A fully pure embedded power system power flow method considering engineering constraints according to claim 6, characterized in that: In step S2), after a stable equilibrium manifold of the pseudo-gradient system is calculated, the escape point is calculated to jump out of the convergence domain of the current stable equilibrium manifold to find another adjacent stable equilibrium manifold. By using the vector composed of points uniformly distributed on the unit sphere and the origin as the search direction, the possibility of repeatedly entering the same convergence domain due to concentrated searches from similar directions can be significantly reduced, thereby improving the search efficiency.

8. The method of full embedded power system power flow considering engineering constraints according to claim 6 is characterized by: In step S2), the stable equilibrium manifold of the pseudo-gradient system is calculated using a logarithmic method, which specifically includes the following steps: S211), initialization, setting the order q of the power series partial sum * , the allowable threshold ε EM , the initial interval length of each expansion Δt * , reduction ratio Initial point x0=(x1(0),…,x m+d (0)), initial point t0 = 0; S212), parameterizing by t, expressing the solution function of the pseudo-gradient system as x(t), wherein the solution function x(t) of the pseudo-gradient system satisfies: S213), assuming that the solution function x(t) of the pseudo-gradient system is holomorphic near t=t0, the series expansion of x(t) at t0 is expressed as follows: x k (t)=Σ q≥0 a kq (t-t0) q ,1≤k≤m+d Among them, a kq is the power series coefficient; Assume Q H The series expansion of (x(t)) at t0 is expressed as follows: Q H,k (x(t))=Σ q≥0 Q H,kq (t-t0) q ,1≤k≤m+d, Among them, Q H,k (x(t)) is Q H The kth element of (x(t)), Q H,kq Q H,k The power series coefficients of (x(t)); S214), Substitute the power series into equation (5) and compare the two sides of the equation (t-t0) q ,q≥0, solve the power series coefficient a in turn kq ; S215) Compare both sides of the equation (t-t0) 0 The coefficient of a k0 (k=1,…,m+d) we have: a k0 =x k (t0),1≤k≤m+d; S216) Compare both sides of the equation (t-t0) q , the coefficient of q≥1, we get: qa kq =Q H,k(q-1) ; S217), assuming that the corresponding coefficient a is calculated kq ,q≤q * , using a truncated power series: And calculate the corresponding rational approximation and the corresponding derivative, namely: Where L and N are the orders of the numerator and denominator of the rational approximation, respectively, and satisfy L+N=q * ; ′ represents the derivative, δ kq are the coefficients of the numerator and denominator of the rational approximation, and satisfy δ k0 =1; S218), let the initial interval length Δt = Δt * >0; Check whether the maximum error of equation (5) at t∈[t0,t0+Δt] is less than the preset allowable error threshold ∈ EM ; If it is satisfied, the effective interval is determined to be [t0, t0+Δt]; otherwise, the interval length is reduced to Until satisfied; update t0←t0+Δt; S219), repeat steps S213) to S218) until the energy function is satisfied The value at t0 is less than 0.1 at the initial point, and then the local solver is used to calculate Q H The solution of (x)=0.

9. A fully pure embedded power system power flow method considering engineering constraints according to claim 8, characterized in that: In step S2), another adjacent stable equilibrium manifold is found by calculating the escape point, which specifically includes the following steps: S221), initialize the initial stable equilibrium point X s , search path direction d, search step Δt>0, search stop maximum value F Lim >0; S222), from x s Starting from the search path direction d, a small perturbation point is used as the initial iteration point Assume the number of iterations i = 0; S223), calculation S224), respectively calculated in The inner product of the derivative of the search path at the point and the vector field of the pseudo-gradient system, that is, d and -DH(x) T The inner product of H(x) is denoted by S225), update i←i+1; S226), repeat steps S223) to S225) until one of the following conditions is met; 1) Returns successful solution and 2) Function The value of or the value of i is greater than the given F Lim , returns the solution failed and 10. A fully embedded power system power flow method considering engineering constraints according to claim 9, characterized in that: In step S3), all stable equilibrium manifolds are searched layer by layer, and the corresponding power flow solutions are selected from the obtained stable equilibrium manifolds, including the following steps: S31), initialize the initial point x0, search step Δt>0, search stop maximum value F Lim >0, number of search paths N * >0; S32), starting from the initial point x0, a stable equilibrium manifold of the pseudo-gradient system calculated in step S2) is obtained S33), set the stable equilibrium manifold set to The set of stable equilibrium manifolds at level 0 is The escape point set is S34), set the stable equilibrium manifold set of the k+1th layer to be S35) Every stable equilibrium manifold in Find the stable equilibrium manifolds of all its neighboring layers; S36) Judgment Is it equal to If yes, go to step S37); otherwise, k←k+1, go to step S34); S37) for V SEM Every stable equilibrium manifold x in s , determine whether it is a solution to the polynomial system H(x)=0.

Images