Probabilistic transient stability constrained optimal power flow considering uncertainty of load

Abstract

A probabilistic transient stability constrained optimal power flow (TSCOPF) method considering the uncertainty of power sources and loads is proposed, which includes the following steps: Step 1: initialize the system power flow parameters and establish the probability distribution model of the uncertainty variables; Step 2: convert the static security inequality constraints and transient stability constraints into the form of probability constraints, and build the TSCOPF model based on the chance-constrained optimization theory; Step 3: convert the sampling matrix in the independent standard normal space to the original space; Step 4: perform multiple deterministic power flow calculations on the sampling matrix based on the point estimate method (PEM), and combine the Cornish-Fisher series to determine whether the probability constraints of each output variable are out of limits; Step 5: convert the original optimization problem into an unconstrained optimization problem, and initialize the parameters of the moth flame optimization (MFO) algorithm; Step 6: calculate the fitness value of the moth, and output the moth with the best fitness value as the optimal solution of the model.

Description

A method for obtaining optimal power flow under probabilistic transient stability constraints considering source load uncertainty Technical Field

[0001] This invention belongs to the field of power system technology, specifically relating to the technical fields of power system power flow calculation and transient stability analysis, and particularly to a method for obtaining optimal power flow under probabilistic transient stability constraints considering source-load uncertainty. Background Technology

[0002] With the rapid development of new energy sources, such as wind power, in the power grid, the inherent randomness and intermittency of these sources have brought more uncertainties to the power system, posing new challenges to its safe and stable operation. Furthermore, the large-scale integration of new loads, such as electric vehicle charging facilities, has exacerbated the uncertainty of grid load. Therefore, correctly analyzing the impact of uncertainties in new energy output and load on the safe and stable operation of the power system, and providing optimal power output strategies for generator units, is of great significance for enhancing the security of the power system.

[0003] Traditional Transient Stability Constrained Optimal Power Flow (TSCOPF) models typically use the minimum system operating cost as the objective function. Under the premise of satisfying the static stability constraints and power balance of the system, they improve the transient stability of the power system by changing generator output and adjusting load configuration, so that the system can still operate safely and stably after a fault.

[0004] However, the traditional TSCOPF model cannot take into account the uncertainties of renewable energy output and load, which often leads to its solution results not matching the actual operating state of the power grid, making it difficult to apply to new power systems with a high proportion of renewable energy access.

[0005] Patent document CN106849058B discloses a transient stability prevention and control auxiliary decision-making method based on a security domain. This patent uses a practical dynamic security domain for transient stability assessment and iteratively calculates the control strategy for node injection power, finally solving for an additional control strategy. While this method improves the solution speed and accuracy of transient stability prevention and control auxiliary decision-making, it does not comprehensively consider the uncertainties of renewable energy output and load during the solution process, making it difficult to apply to new power systems with large-scale renewable energy integration. Summary of the Invention

[0006] The purpose of this invention is to address the shortcomings of traditional TSCOPF (Transient Stability-Free Power Utilization) models, which cannot account for the uncertainties in renewable energy output and load. Based on mathematical statistics theory, this invention analyzes the different characteristics of three uncertain variables: wind power output, photovoltaic power output, and load power, and uses different mathematical methods to describe their probability distribution characteristics. Simultaneously, it analyzes the transfer and transmission of the probability distribution characteristics of each uncertain variable in traditional TSCOPF, and constructs a probabilistic TSCOPF model that accounts for uncertain variables by combining chance constraint theory. Furthermore, it employs a suitable optimization algorithm to solve the model.

[0007] A probabilistic transient stability constraint-based optimal power flow acquisition method considering source load uncertainty includes the following steps:

[0008] Step 1: Initialize system power flow parameters and establish a probability distribution model for uncertain variables;

[0009] Step 2: Convert the static security inequality constraints and transient stability constraints into probabilistic constraints, and construct the probabilistic transient stability constraint optimal power flow TSCOPF model based on chance constraint optimization theory;

[0010] Step 3: Using the principle of Nataf inverse transformation, calculate the correlation coefficient matrix and lower triangular matrix step by step to transform the sampling matrix in the independent standard normal space into the original space sampling matrix;

[0011] Step 4: Perform multiple deterministic power flow calculations on the sampling matrix based on the point estimation method (PEM), and combine the Cornish-Fisher series to determine whether the probability constraints of each output variable are exceeded.

[0012] Step 5: Transform the original optimization problem into an unconstrained optimization problem and initialize the parameters of the Moth to a Flame (MFO) optimization algorithm;

[0013] Step 6: Calculate the fitness value of the moths and output the moth with the best fitness value as the optimal solution of the model.

[0014] In step 1, the system power flow parameters are initialized. The normal distribution is used to describe the active and reactive power of the load, the two-parameter Weibull distribution is used to describe the active power output of wind power generation, and the Beta distribution is used to describe the active power output of photovoltaic power generation.

[0015] In step 1, the load is affected by uncertainties such as time, weather, and user electricity consumption habits, exhibiting a certain degree of randomness. A normal distribution is used to describe the active and reactive power of the load at node i, and its probability density function is shown below:

[0016]

[0017] In the formula: P and Q represent the active power and reactive power of the load, respectively; μ P and σ P μ represents the mean and variance of the active power of the load, respectively. Q and σ Q These represent the mean and variance of the reactive power of the load, respectively.

[0018] Wind speed is easily affected by various factors such as weather, terrain, and season, and has inherent uncertainty. A two-parameter Weibull distribution is used to describe the wind speed input of the wind turbine, and its probability density function is shown below:

[0019]

[0020] In the formula: v represents wind speed; k and c represent shape parameter and scale parameter, respectively;

[0021] The output of a wind turbine is determined by its speed-power curve.

[0022]

[0023] In the formula: P W (v) represents the active power output of the wind turbine when the wind speed is v; p r Indicates the rated power of the wind turbine generator; v ci v r and v co These are respectively represented as cut-in wind speed, rated wind speed, and cut-out wind speed;

[0024] Photovoltaic power generation mainly depends on the intensity of solar radiation. The active power output of photovoltaic power generation is described using a Beta distribution, and its probability density function is shown below:

[0025]

[0026] In the formula: P PVmax P represents the maximum active power output of photovoltaic power generation. PV The active power output of photovoltaic power generation is represented by α and β, which are parameters that adjust the shape of the Beta distribution; Γ() represents the Gamma function.

[0027] In step 2,

[0028] The traditional static security inequality constraints and transient stability constraints are transformed into probabilistic constraints, and a probabilistic TSCOPF model is constructed based on chance-constrained optimization theory, as follows:

[0029] 1) Objective function

[0030] When the TSCOPF model considers the uncertainties in renewable energy output and load, the objective function of the entire model will be expressed in the form of expected value:

[0031]

[0032] In the formula: F G Represents the cost of all generator sets; E() represents the expected value; N represents the set of all generator sets; P Gi a i Let P represent the active power output and cost function of the i-th traditional unit, respectively; wj b j Let P represent the active power output and cost function of the j-th wind turbine, respectively; PVk c k Let represent the active power output and cost function of the k-th photovoltaic unit, respectively.

[0033] 2) Power balance equation

[0034]

[0035] Where: n b Q represents the total number of nodes in the system; Gi Q represents the reactive power output of the traditional unit at node i; wi Q represents the reactive power output of the wind turbine at node i; PVi P represents the reactive power output of the photovoltaic unit at node i; Di and Q Di V represents the active and reactive power of the load at node i; i and V j G represents the voltage magnitudes at nodes i and j; ij and B ij θ represents the imaginary and real parts of the admittance of branch ij; ij This represents the voltage phase angle difference between node i and node j.

[0036] 3) Static safety inequality probability constraints

[0037]

[0038] In the formula: P{} represents the probability that the condition within the parentheses is true; β P β Q β V and β S This represents the probability threshold for the system to maintain static security. These represent the lower and upper limits of the active power output of a traditional generator unit, respectively. These represent the lower and upper limits of the reactive power output of a traditional generator set, respectively; Vi min V i max These represent the lower and upper limits of the node voltage amplitude, respectively; S li , Let n represent the transmission power and the upper limit of the transmission power of the i-th line, respectively; G n b n l These represent the traditional generator set, node set, and line set, respectively.

[0039] 4) Transient stability probability constraints

[0040] Traditional transient stability constraints involve numerous nonlinear differential-algebraic equations, resulting in high computational cost and slow solution speed. In contrast, machine learning-based transient stability evaluation methods possess the ability to rapidly process massive amounts of data. Therefore, this paper replaces traditional transient stability constraints with a Lightweight Gradient Boosting Machine (LightGBM) and uses the Transient Stability Index (TSI) as the stability criterion.

[0041]

[0042] Where: Δδ max This represents the maximum relative power angle difference between any two generators during the simulation period. If TSI > 0, the system is stable; if TSI < 0, the system is unstable. When using LightGBM to determine the transient stability of a power system, its probabilistic form can be expressed as:

[0043] P{TSI>0}>β TSI (9)

[0044] Where: β TSI This represents the probability threshold for the transient stability of a power system.

[0045] In step 3,

[0046] Using the principle of inverse Nataf transform, the correlation coefficient matrix and lower triangular matrix are calculated step by step to transform the sampling matrix in the independent standard normal space to the original space sampling matrix. The steps of Nataf transform are as follows:

[0047] Step 3-1: Calculate the linear correlation coefficient matrix ρ0 of the transformed standard normal space sampling matrix Y based on the linear correlation coefficient matrix ρ of the original spatial sampling matrix X;

[0048] Let X = {x1, x2, ..., xn} be an n-dimensional original spatial sampling matrix with correlation.n}

[0049]

[0050] In the formula: Y={y1,y2,...,y n} represents the sampling matrix in the standard normal space; F i (x i ) represents the variable x i The cumulative distribution function; Φ() represents the standard normal cumulative distribution function; Φ -1 () denotes the inverse cumulative distribution function, and the correlation coefficient matrices of X and Y are represented by ρ and ρ0, respectively. The mapping relationship between the ρ and ρ0 components can be described as follows:

[0051]

[0052] Where: σ i μ i and σ j μ j They represent the variables x respectively i and x j Standard deviation and expected value; ρ 0ij Let y represent any two random variables i and y j The correlation coefficient; φ2() is the two-dimensional joint probability density function of the standard normal distribution; x represents i and x j The joint probability density function; when ρ is known, ρ0 can be calculated by solving the nonlinear equation (11);

[0053] Step 3-2: Obtain the lower triangular matrix L0 by performing Choleskey decomposition on matrix ρ0:

[0054]

[0055] In the formula: the superscript "T" indicates matrix transpose;

[0056] Step 3-3: Using L0, the sampling matrix Y in the relevant standard normal space can be transformed into an independent standard normal distribution vector V:

[0057]

[0058] In the formula: the superscript "-1" indicates the inverse of the matrix.

[0059] In step 4, multiple deterministic power flow calculations are performed on the sampling matrix based on the point estimation method (PEM), and the Cornish-Fisher series is used to determine whether the probability constraints of each output variable are exceeded. This includes the following steps:

[0060] Step 4-1: Based on the point estimation method (PEM), perform multiple deterministic power flow calculations on the original spatial sampling matrix to obtain output variables (unit active power output, unit reactive power output, node voltage amplitude, etc.).

[0061] Point estimation method (PEM) obtains the raw moments of each output variable by solving a certain number of deterministic power flows, and then obtains the corresponding expected value and standard deviation. For a system with m-dimensional input variables and n-dimensional output variables, its expression can be written as:

[0062] R = G(X) = G(x1,x2,...,x) m (14)

[0063] In the formula: R = [r1, r2, ..., r m ] T G = [g1, g2, ..., g m ] T .

[0064] For each input random variable x i (i = 1, 2, ..., m), while other variables take the mean, there are three positional parameters x. i,k (k = 1, 2, 3), hence the name of the three-point estimation method. Variable x i Position parameter x i,k and its corresponding weight parameter w xi,k It can be represented as:

[0065]

[0066]

[0067]

[0068] In the formula: and They represent the variables x respectively i The standard location, mean, and standard deviation; They represent the variables x respectively i The skewness and kurtosis coefficients. From equations (15) to (17), it can be seen that the three-point estimation method mainly uses the first four moments of the input variables for calculation. And for each position parameter x... i,kEach of these will be solved deterministically using equation (14) to obtain the corresponding output variable R:

[0069] R(i,k)=G(μ x1 ,...,μ xi-1 ,x i,k ,μ xi+1 ,...,μ m (18)

[0070] In the formula: i = 1, 2, ..., m, k = 1, 2, 3. From equation (18), it can be seen that the three-point estimation method requires solving 3m deterministic power flow problems to obtain the final output vector R. It is worth noting that when k is 3, the standard position... Then the position parameters are obtained. That is, there are m solutions to equation (18) all using the mean value of the input variable, so the total number of solutions 3m can be reduced to 2m+1.

[0071] Step 4-2: Calculate the raw moments of each order of the output variables, and then obtain the corresponding expected values ​​and standard deviations.

[0072]

[0073] In the formula: Indicates the output variable r j The l-th moment. The output variable r is calculated... j The first two moments can be further used to obtain the corresponding expected value μ. rj and standard deviation σ rj :

[0074] μ rj =E(r) j j=1,2,...,m (20)

[0075]

[0076] Step 4-3: Obtain the cumulative distribution function of each output variable based on the Cornish-Fisher series fitting, and determine whether each probability constraint is exceeded.

[0077] In step 5, the power balance equation in the probabilistic TSCOPF model can be indirectly processed in the process of solving the deterministic power flow using the point estimation method PEM. The static security inequality probabilistic constraint and the probabilistic transient stability margin are processed by using the penalty function to handle the probabilistic form of the constraint and combined with the objective function to transform the original optimization problem into an unconstrained optimization problem, as shown in equation (22), and then solved using the moth-to-a-flame optimization algorithm MFO.

[0078]

[0079] In the formula: γ TSI γ P γ Q γ V γ s represents the penalty factor for the probability constraint; u represents the system control variable.

[0080] In step 5, the Moths Attracting Fire (MFO) optimization algorithm is used to solve the unconstrained optimization problem of equation (22). For the initialization of the MFO algorithm, parameters such as the control variable dimension d, the search size of the moth population n, the maximum number of iterations T, and the logarithmic spiral shape constant b are set. Moth positions are randomly generated in the search space, and the fitness value corresponding to each moth is evaluated.

[0081] The Moth to Flame (MFO) optimization algorithm differs from some other metaheuristic algorithms in that it contains two important solution populations (moths and flames). The moth represents the position of the current solution randomly generated within the search domain. The flame represents the position of a better solution obtained by the moth. Matrix M is used to represent the moth:

[0082]

[0083] In the formula: n represents the number of moths; d represents the dimension. For all moths, an array is set to store the corresponding fitness values ​​as follows:

[0084] OM = [OM1OM2…OM] n ] T (twenty four)

[0085] Flames can also be represented by a matrix F:

[0086]

[0087] Similarly, an array is set up to store the corresponding fitness values ​​of the flames:

[0088] OF = [OF1OF2…OF] n ] T (26)

[0089] Both moths and flames are solutions, but they employ different update methods during evolution. Each moth is assigned to a flame, and its position is updated by rotating around the designated flame, as shown below:

[0090] M i =S(M i ,F j(27)

[0091] In the formula: S represents the spiral function. Basically, any spiral function that meets certain conditions can be used. In the original Moth to Fire Optimization (MFO) algorithm, a logarithmic spiral is chosen. The position of each moth is updated according to the corresponding flame as follows:

[0092] S(M i ,F j ) = D i ·e br ·cos(2πr)+F j (28)

[0093] D i =|M i -F j | (29)

[0094] r = a·rand + 1 (30)

[0095]

[0096] In the formula: D i denoted by , b represents the distance between the i-th moth and the j-th flame; b represents a constant defining the logarithmic spiral shape; r represents the distance parameter, which defines the distance between the next position of the moth and the flame; T represents the maximum number of iterations; t represents the current iteration number; a decreases linearly from -1 to -2 during the evolution process.

[0097] In step 6, the best fitness value among all moths is calculated, and it is determined whether the maximum number of iterations has been reached: if it has, the best fitness value of the moth is output as the optimal solution of the probabilistic TSCOPF model; otherwise, the fitness values ​​of the updated moth positions and flame positions are reordered, and the spatial position with the better fitness value is selected as the position of the next generation flame. The positions of the moths in the Moth-to-Flame Optimization Algorithm (MFO) are continuously updated until the number of iterations reaches the algorithm requirement.

[0098] An optimal power flow model with probabilistic transient stability constraints considering source load uncertainty is presented. The model includes an objective function, power balance equations, probabilistic constraints based on static security inequalities, and probabilistic constraints on transient stability, as detailed below:

[0099] 1) Objective function

[0100] When the TSCOPF model considers the uncertainties in renewable energy output and load, the objective function of the entire model will be expressed in the form of expected value:

[0101]

[0102] In the formula: F G E represents the cost of all generating units; E(·) represents the expected value; N represents the set of all generating units; P Gi a i Let P represent the active power output and cost function of the i-th traditional unit, respectively; wj b j Let P represent the active power output and cost function of the j-th wind turbine, respectively; PVk c k Let represent the active power output and cost function of the k-th photovoltaic unit, respectively;

[0103] 2) Power balance equation

[0104]

[0105] Where: n b Q represents the total number of nodes in the system; Gi Q represents the reactive power output of the traditional unit at node i; wi Q represents the reactive power output of the wind turbine at node i; PVi P represents the reactive power output of the photovoltaic unit at node i; Di and Q Di V represents the active and reactive power of the load at node i; i and V j G represents the voltage magnitudes at nodes i and j; ij and B ij θ represents the imaginary and real parts of the admittance of branch ij; ij This represents the voltage phase angle difference between node i and node j;

[0106] 3) Static safety inequality probability constraints

[0107]

[0108] In the formula: P{·} represents the probability that the condition within the parentheses is true; β P β Q β V and β S This represents the probability threshold for the system to maintain static security. These represent the lower and upper limits of the active power output of a traditional generator unit, respectively. These represent the lower and upper limits of the reactive power output of a traditional generator set, respectively; V i min V i max These represent the lower and upper limits of the node voltage amplitude, respectively; S li , Let n represent the transmission power and the upper limit of the transmission power of the i-th line, respectively;G n b n l These represent the traditional generator set, node set, and line set, respectively.

[0109] 4) Transient stability probability constraints

[0110]

[0111] Where: Δδ max This represents the maximum relative power angle difference between any two generators during the simulation period. If TSI > 0, the system is stable; if TSI < 0, the system is unstable. When using LightGBM to determine the transient stability of a power system, its probabilistic form can be expressed as:

[0112] P{TSI>0}>β TSI (36)

[0113] Where: β TSI This represents the probability threshold for the transient stability of a power system.

[0114] A probabilistic transient stability-constrained optimal power flow model considering source load uncertainty employs the following steps when solving:

[0115] Step 1: Using the principle of Nataf inverse transformation, calculate the correlation coefficient matrix and the lower triangular matrix step by step to transform the sampling matrix in the independent standard normal space into the original space sampling matrix;

[0116] Step 2: Perform multiple deterministic power flow calculations on the sampling matrix based on the point estimation method (PEM), and combine the Cornish-Fisher series to determine whether the probability constraints of each output variable are exceeded.

[0117] Step 3: Transform the original optimization problem into an unconstrained optimization problem and initialize the parameters of the Moth to a Flame (MFO) optimization algorithm;

[0118] Step 4: Calculate the fitness value of the moths and output the moth with the best fitness value as the optimal solution of the model.

[0119] Compared with the prior art, the present invention has the following technical effects:

[0120] 1) To address the shortcomings of the traditional TSCOPF model in failing to account for uncertainties in renewable energy output and load, this invention uses mathematical statistics theory and different mathematical methods to describe the probability distribution characteristics of three uncertain variables: wind power output, photovoltaic power output, and load power. Furthermore, it replaces the deterministic constraints in the traditional TSCOPF model with static security inequality probability constraints and transient stability probability constraints, and constructs a probabilistic TSCOPF model based on chance constraint optimization theory.

[0121] 2) In the process of solving the probabilistic TSCOPF model, this invention combines the Nataf inverse transformation principle and the point estimation method PEM theory to perform multiple deterministic power flow solutions. It also uses multiple deterministic variables obtained in the solution process to equivalently characterize the three uncertain variables of wind power output, photovoltaic power output and load power, thereby transforming the complex uncertain problem into a deterministic problem and reducing the difficulty of solving the model.

[0122] 3) This invention uses a penalty function to handle probabilistic constraints and combines them with the objective function to transform the original optimization problem into an unconstrained optimization problem. This is used to solve the static safety inequality probability constraints and transient stability probability constraints in the probabilistic TSCOPF model. At the same time, the position of the moth in the Moth-to-Flame Optimization Algorithm (MFO) is initialized as the control variable of the model. The optimal fitness value of the moth is obtained through iterative optimization, which is the optimal solution of the probabilistic TSCOPF model, thereby improving the solution speed and accuracy of the model. Attached Figure Description

[0123] The present invention will be further described below with reference to the accompanying drawings and embodiments:

[0124] Figure 1 is a flowchart of the method of the present invention;

[0125] Figure 2 is a schematic diagram of the IEEE 39-node system in an example of the present invention;

[0126] Figure 3 is a relative rotor angle diagram of the IEEE 39-bus system under a anticipated fault in an example of the present invention;

[0127] Figure 4 is a diagram of the relative rotor angle of the IEEE 39-node system after calculation according to the present invention in an example of the present invention. Detailed Implementation

[0128] A probabilistic transient stability-constrained optimal power flow model and calculation method considering source load uncertainty includes the following steps, as shown in Figure 1:

[0129] Step 1: Initialize system power flow parameters and establish a probability distribution model for uncertain variables;

[0130] Step 2: Convert the static security inequality constraints and transient stability constraints into probabilistic constraints, and construct a probabilistic TSCOPF model based on chance-constrained optimization theory;

[0131] Step 3: Using the principle of Nataf inverse transformation, calculate the correlation coefficient matrix and lower triangular matrix step by step to transform the sampling matrix in the independent standard normal space into the original space sampling matrix;

[0132] Step 4: Perform multiple deterministic power flow calculations on the sampling matrix based on the point estimation method (PEM), and combine the Cornish-Fisher series to determine whether the probability constraints of each output variable are exceeded.

[0133] Step 5: Transform the original optimization problem into an unconstrained optimization problem and initialize the parameters of the Moth to a Flame (MFO) optimization algorithm;

[0134] Step 6: Calculate the fitness value of the moths and output the moth with the best fitness value as the optimal solution of the model.

[0135] In step 1: initialize the system power flow parameters, use a normal distribution to describe the active and reactive power of the load, a two-parameter Weibull distribution to describe the active power output of wind power generation, and a Beta distribution to describe the active power output of photovoltaic power generation.

[0136] Furthermore, in step 1: the load is affected by uncertainties such as time, weather, and user electricity consumption habits, exhibiting a certain degree of randomness. A normal distribution is used to describe the active and reactive power of the load at node i, and its probability density function is shown below:

[0137]

[0138] In the formula: P and Q represent the active power and reactive power of the load, respectively; μ P and σ P μ represents the mean and variance of the active power of the load, respectively. Q and σ Q These represent the mean and variance of the reactive power of the load, respectively.

[0139] Wind speed is easily affected by various factors such as weather, terrain, and season, exhibiting inherent uncertainty. A two-parameter Weibull distribution is used to describe the wind speed input to the wind turbine, and its probability density function is shown below:

[0140]

[0141] In the formula: v represents wind speed; k and c represent shape parameter and scale parameter, respectively.

[0142] The output of a wind turbine is determined by its speed-power curve.

[0143]

[0144] In the formula: P W (v) represents the active power output of the wind turbine when the wind speed is v; p r Indicates the rated power of the wind turbine generator; v ci vr and v co These are represented as cut-in wind speed, rated wind speed, and cut-out wind speed, respectively.

[0145] Photovoltaic power generation primarily depends on the intensity of solar radiation and is also subject to uncertainty. The Beta distribution is used to describe the active power output of photovoltaic power generation, and its probability density function is shown below:

[0146]

[0147] In the formula: P represents the maximum active power output of photovoltaic power generation. PV The active power output of photovoltaic power generation is represented by α and β, which are parameters that adjust the shape of the Beta distribution; Γ() represents the Gamma function.

[0148] In step 2: To fully consider the impact of source load uncertainty on the traditional TSCOPF, the traditional static security inequality constraints and transient stability constraints are converted into probabilistic constraints, and a probabilistic TSCOPF model is constructed based on chance-constrained optimization theory, as follows:

[0149] 1) Objective function

[0150] When the TSCOPF model considers the uncertainties in renewable energy output and load, the objective function of the entire model will be expressed in the form of expected value:

[0151]

[0152] In the formula: F G Represents the cost of all generator sets; E() represents the expected value; N represents the set of all generator sets; P Gi a i Let P represent the active power output and cost function of the i-th traditional unit, respectively; wj b j Let P represent the active power output and cost function of the j-th wind turbine, respectively; PVk c k Let represent the active power output and cost function of the k-th photovoltaic unit, respectively.

[0153] 2) Power balance equation

[0154]

[0155] Where: n b Q represents the total number of nodes in the system; Gi Q represents the reactive power output of the traditional unit at node i; wi Q represents the reactive power output of the wind turbine at node i;PVi P represents the reactive power output of the photovoltaic unit at node i; Di and Q Di V represents the active and reactive power of the load at node i; i and V j G represents the voltage magnitudes at nodes i and j; ij and B ij θ represents the imaginary and real parts of the admittance of branch ij; ij This represents the voltage phase angle difference between node i and node j.

[0156] 3) Static safety inequality probability constraints

[0157]

[0158] In the formula: P{} represents the probability that the condition within the parentheses is true; β P β Q β V and β S This represents the probability threshold for the system to maintain static security. These represent the lower and upper limits of the active power output of a traditional generator unit, respectively. These represent the lower and upper limits of the reactive power output of a traditional generator set, respectively; V i min V i max These represent the lower and upper limits of the node voltage amplitude, respectively; S li , Let n represent the transmission power and the upper limit of the transmission power of the i-th line, respectively; G n b n l These represent the traditional generator set, node set, and line set, respectively.

[0159] 4) Transient stability probability constraints

[0160] Traditional transient stability constraints involve numerous nonlinear differential-algebraic equations, resulting in high computational cost and slow solution speed. In contrast, machine learning-based transient stability evaluation methods possess the ability to rapidly process massive amounts of data. Therefore, this paper replaces traditional transient stability constraints with a Lightweight Gradient Boosting Machine (LightGBM) and uses the Transient Stability Index (TSI) as the stability criterion.

[0161]

[0162] Where: Δδ maxThis represents the maximum relative power angle difference between any two generators during the simulation period. If TSI > 0, the system is stable; if TSI < 0, the system is unstable. When using LightGBM to determine the transient stability of a power system, its probabilistic form can be expressed as:

[0163] P{TSI>0}>β TSI (9)

[0164] Where: β TSI This represents the probability threshold for the transient stability of a power system.

[0165] For the above probabilistic TSCOPF model, the solution is obtained by combining the Nataf inverse transform, the point estimation method (PEM), and the moth-to-a-flame optimization algorithm (MFO).

[0166] In step 3: Using the principle of Nataf inverse transform, the correlation coefficient matrix and lower triangular matrix are calculated step by step to transform the sampling matrix in the independent standard normal space to the original space sampling matrix. The steps of Nataf transform are as follows:

[0167] Step 3-1: Calculate the linear correlation coefficient matrix ρ0 of the sample matrix Y in the transformed standard normal space based on the linear correlation coefficient matrix ρ of the original spatial sampling matrix X.

[0168] Let X = {x1, x2, ..., xn} be an n-dimensional original spatial sampling matrix with correlation. n}

[0169]

[0170] In the formula: Y={y1,y2,...,y n} represents the sampling matrix in the standard normal space; F i (x i ) represents the variable x i The cumulative distribution function; Φ() represents the standard normal cumulative distribution function; Φ -1 () denotes the inverse cumulative distribution function. The correlation coefficient matrices of X and Y are denoted by ρ and ρ0, respectively, and the mapping relationship between the ρ and ρ0 components can be described as:

[0171]

[0172] Where: σ i μ i and σ j μ j They represent the variables x respectively i and x j Standard deviation and expected value; ρ 0ij Let y represent any two random variablesi and y j The correlation coefficient; φ2() is the two-dimensional joint probability density function of the standard normal distribution; x represents i and x j The joint probability density function; when ρ is known, ρ0 can be calculated by solving the nonlinear equation (11).

[0173] Step 3-2: Obtain the lower triangular matrix L0 by performing Choleskey decomposition on matrix ρ0:

[0174]

[0175] In the formula: the superscript "T" indicates matrix transpose.

[0176] Step 3-3: Using L0, the sampling matrix Y in the relevant standard normal space can be transformed into an independent standard normal distribution vector V:

[0177]

[0178] In the formula: the superscript "-1" indicates the inverse of the matrix.

[0179] The Nataf inverse transform is the reverse process of equations (10) and (13).

[0180] In step 4: Based on the point estimation method (PEM), multiple deterministic power flow calculations are performed on the sampling matrix, and the Cornish-Fisher series is used to determine whether the probability constraints of each output variable are exceeded.

[0181] Step 4-1: Based on the point estimation method (PEM), perform multiple deterministic power flow calculations on the original spatial sampling matrix to obtain output variables (unit active power output, unit reactive power output, node voltage amplitude, etc.).

[0182] Point estimation method (PEM) obtains the raw moments of each output variable by solving a certain number of deterministic power flows, and then obtains the corresponding expected value and standard deviation. For a system with m-dimensional input variables and n-dimensional output variables, its expression can be written as:

[0183] R = G(X) = G(x1,x2,...,x) m (14)

[0184] In the formula: R = [r1, r2, ..., r m ] T G = [g1, g2, ..., g m ] T .

[0185] For each input random variable x i (i = 1, 2, ..., m), while other variables take the mean, there are three positional parameters x. i,k (k = 1, 2, 3), hence the name of the three-point estimation method. Variable x i Position parameter x i,k and their corresponding weight parameters It can be represented as:

[0186]

[0187]

[0188]

[0189] In the formula: and They represent the variables x respectively i The standard location, mean, and standard deviation; They represent the variables x respectively i The skewness and kurtosis coefficients. From equations (15) to (17), it can be seen that the three-point estimation method mainly uses the first four moments of the input variables for calculation. And for each position parameter x... i,k Each of these will be solved deterministically using equation (14) to obtain the corresponding output variable R:

[0190] R(i,k)=G(μ x1 ,...,μ xi-1 ,x i,k ,μ xi+1 ,...,μ m (18)

[0191] In the formula: i = 1, 2, ..., m, k = 1, 2, 3. From equation (18), it can be seen that the three-point estimation method requires solving 3m deterministic power flow problems to obtain the final output vector R. It is worth noting that when k is 3, the standard position... Then the position parameters are obtained. That is, there are m solutions to equation (18) all using the mean value of the input variable, so the total number of solutions 3m can be reduced to 2m+1.

[0192] Step 4-2: Calculate the raw moments of each order of the output variables, and then obtain the corresponding expected values ​​and standard deviations.

[0193]

[0194] In the formula: Indicates the output variable r j The l-th moment. The output variable r is calculated... j The first two moments can be further used to obtain the corresponding expected value μ. rj and standard deviation σ rj :

[0195] μ rj =E(r) j j=1,2,...,m (20)

[0196]

[0197] Step 4-3: Obtain the cumulative distribution function of each output variable based on the Cornish-Fisher series fitting, and determine whether each probability constraint is exceeded.

[0198] In step 5: the power balance equation in the probabilistic TSCOPF model can be indirectly processed in the process of solving the deterministic power flow using the point estimation method (PEM), while the static security inequality probabilistic constraints and probabilistic transient stability margin need to be processed separately. By using the penalty function to process the probabilistic constraints and combining them with the objective function, the original optimization problem is transformed into an unconstrained optimization problem, as shown in equation (22), and then solved using the Moth to a Flame optimization algorithm (MFO).

[0199]

[0200] In the formula: γ TSI γ P γ Q γ V γ s represents the penalty factor for the probability constraint; u represents the system control variable.

[0201] Furthermore, in step 5: the Moth-to-Flame Optimization (MFO) algorithm is used to solve the unconstrained optimization problem of equation (22). For the initialization of the MFO algorithm, parameters such as the control variable dimension d, the search size of the moth population n, the maximum number of iterations T, and the logarithmic spiral shape constant b are set. Moth positions are randomly generated in the search space, and the fitness value corresponding to each moth is evaluated.

[0202] The Moth to Flame (MFO) optimization algorithm differs from some other metaheuristic algorithms in that it contains two important solution populations (moths and flames). The moth represents the position of the current solution randomly generated within the search domain. The flame represents the position of a better solution obtained by the moth. Matrix M is used to represent the moth:

[0203]

[0204] In the formula: n represents the number of moths; d represents the dimension. For all moths, an array is set to store the corresponding fitness values ​​as follows:

[0205] OM = [OM1OM2…OM] n ] T (twenty four)

[0206] Flames can also be represented by a matrix F:

[0207]

[0208] Similarly, an array is set up to store the corresponding fitness values ​​of the flames:

[0209] OF = [OF1OF2…OF] n ] T (26)

[0210] Both moths and flames are solutions, but they employ different update methods during evolution. Each moth is assigned to a flame, and its position is updated by rotating around the designated flame, as shown below:

[0211] M i =S(M i ,F j (27)

[0212] In the formula: S represents the spiral function. Basically, any spiral function that meets certain conditions can be used. In the original Moth to Fire Optimization (MFO) algorithm, a logarithmic spiral is chosen. The position of each moth is updated according to the corresponding flame as follows:

[0213] S(M i ,F j ) = D i ·e br ·cos(2πr)+F j (28)

[0214] D i =M i -F j (29)

[0215] r = a·rand + 1 (30)

[0216]

[0217] In the formula: D idenoted by , b represents the distance between the i-th moth and the j-th flame; b represents a constant defining the logarithmic spiral shape; r represents the distance parameter, which defines the distance between the next position of the moth and the flame; T represents the maximum number of iterations; t represents the current iteration number; a decreases linearly from -1 to -2 during the evolution process.

[0218] In step 6: Calculate the best fitness value among all moths and determine whether the maximum number of iterations has been reached. If it has, output the best fitness value of the moth as the optimal solution of the probabilistic TSCOPF model. Otherwise, reorder the fitness values ​​of the updated moth positions and flame positions, select the spatial position with the better fitness value to update the position of the next generation flame, and continuously update the position of the moth in the Moth-to-Flame Optimization Algorithm (MFO) until the number of iterations reaches the algorithm requirement.

[0219] Example:

[0220] To verify the effectiveness of this invention, it was tested on the IEEE 39-node system, in which photovoltaic (PV) generators and wind turbine generators were added. As shown in Figure 2, synchronous turbine generators G3 and G4 were replaced with wind turbine generators. The shape parameter of the Weibull distribution of wind speed was set to k = 10.7, and the scale parameter was set to c = 3.97. The cut-in wind speed, rated wind speed, and cut-out wind speed of the wind turbine generators were 3 m / s, 13 m / s, and 25 m / s, respectively. Synchronous turbine generator G8 was replaced with a PV generator, and the parameters of the PV output Beta distribution were set to α = 0.7 and β = 2.16.

[0221] A three-phase short-circuit fault was set on line 10-13 of the IEEE 39-bus system, lasting 0.2 seconds. The relative rotor angles of the synchronous generators in the system at this time are shown in Figure 3. It can be seen that under this anticipated fault, the relative rotor angles between the generators continuously increased, TSI = -99.72, and the system experienced transient instability. Using the probabilistic TSCOPF model and calculation method proposed in this invention, considering source-load uncertainty, the optimal solution was obtained by iterating the Moth to Fire (MFO) optimization algorithm 100 times to obtain the moth with the best fitness value. The output of each synchronous generator was then readjusted. As shown in Figure 4, the relative rotor angles between the generators in the system no longer continuously increased, TSI = 67.23, and the system was transiently stable. The output of each synchronous generator before and after the adjustment is shown in Table 1.

[0222]

[0223]

[0224] Table 1

[0225] The test results show that the probabilistic transient stability constrained optimal power flow model and calculation method proposed in this invention, which considers source-load uncertainty, is applicable to new power systems with a high proportion of new energy access and is of great significance for ensuring the safe and stable operation of the power system.

Claims

1. A method for obtaining optimal power flow under probabilistic transient stability constraints considering source load uncertainty, characterized in that: It includes the following steps: Step 1: Initialize system power flow parameters and establish a probability distribution model for uncertain variables; Step 2: Convert static security inequality constraints and transient stability constraints into probabilistic constraint forms, and construct a probabilistic transient stability constraint optimal power flow model (TSCOPF) based on chance constraint optimization theory; Step 3: Utilize the Nataf inverse transform principle to progressively calculate the correlation coefficient matrix and lower triangular matrix, transforming the sampling matrix in the independent standard normal space to the original space sampling matrix; Step 4: Perform multiple deterministic power flow calculations on the sampling matrix based on the point estimation method (PEM), and combine the Cornish-Fisher series to determine the probability distribution of each output variable. Step 5: Convert the original optimization problem into an unconstrained optimization problem and initialize the parameters of the Moth to Fire Optimization Algorithm (MFO). Step 6: Calculate the fitness value of the moth and output the moth with the best fitness value as the optimal solution of the model. In step 5, the power balance equation in the probabilistic TSCOPF model can be indirectly processed in the process of solving the deterministic power flow using the point estimation method (PEM). The static security inequality probability constraint and the probabilistic transient stability margin are processed by using the penalty function to handle the probabilistic form of the constraint and combined with the objective function to convert the original optimization problem into an unconstrained optimization problem, as shown in equation (22), and then solved using the Moth to Fire Optimization Algorithm (MFO). (22); where: 、 、 、 、 The penalty factor representing the probability constraint; Represents system control variables; This indicates the cost of all units; 、 、 These represent the traditional generator set, node set, and line set, respectively. This indicates the probability that the transient stability index is greater than zero. 、 、 and This represents the probability threshold for the system to maintain static security. The probability threshold representing the transient stability of a power system; This indicates the probability that the condition within the parentheses is true. 、 These represent the lower and upper limits of the active power output of a traditional generator unit, respectively. 、 These represent the lower and upper limits of reactive power output of traditional generator units, respectively. 、 These represent the lower and upper limits of the node voltage amplitude, respectively; 、 Let represent the transmission power and the upper limit of the transmission power of the i-th line, respectively; This represents the active power output of the i-th conventional unit. This represents the reactive power output of the traditional unit at node i. This represents the voltage magnitude at node i; This is the penalty function form of the transient stability probability constraint; This is a penalty for exceeding the limit on the probability of generator active power output; This is a penalty for exceeding the limit on the probability of generator reactive power output; This is a penalty term for the probability of node voltage exceeding the limit; This is a penalty for exceeding the line capacity probability limit.

2. The method according to claim 1, characterized in that, In step 1, the system power flow parameters are initialized. The normal distribution is used to describe the active and reactive power of the load, the two-parameter Weibull distribution is used to describe the active power output of wind power generation, and the Beta distribution is used to describe the active power output of photovoltaic power generation.

3. The method according to claim 1, characterized in that, In step 1, the load is affected by uncertainties such as time, weather, and user electricity consumption habits, exhibiting a certain degree of randomness. A normal distribution is used to describe the active and reactive power of the load at node i, and its probability density function is shown below: (1); where: and These represent the active power and reactive power of the load, respectively. and These represent the mean and variance of the active power of the load, respectively. and Let represent the mean and variance of the reactive power load, respectively. Wind speed is easily affected by factors such as weather, terrain, and season, exhibiting inherent uncertainty. A two-parameter Weibull distribution is used to describe the wind speed input of the wind turbine, and its probability density function is shown below: (2); where: Indicates wind speed; 、 These represent the shape parameter and the dimensional parameter, respectively; the output of the wind turbine is determined by its speed-power curve. (3); where: Indicates wind speed as Active power output of the wind turbine unit; This indicates the rated power of the wind turbine generator set; 、 and Let these be the cut-in wind speed, rated wind speed, and cut-out wind speed, respectively. The active power output of photovoltaic power generation is described using a Beta distribution, and its probability density function is shown below: (4); where: This indicates the maximum active power output of photovoltaic power generation; Indicates the active power output of photovoltaic power generation; parameters and It is a parameter that adjusts the shape of the Beta distribution; This represents the Gamma function.

4. The method according to claim 1, characterized in that, In step 2, the traditional static security inequality constraints and transient stability constraints are converted into probabilistic constraints, and a probabilistic TSCOPF model is constructed based on chance-constrained optimization theory, as follows: 1) Objective function: When the TSCOPF model considers the uncertainties of new energy output and load, the objective function of the entire model will be expressed in the form of expected value: (5); where: This indicates the cost of all units; N represents the expected value; N represents the set of all generator sets. 、 Let represent the active power output and cost function of the i-th traditional unit, respectively; 、 Let represent the active power output and cost function of the j-th wind turbine, respectively; 、 Let the active power output and cost function of the k-th photovoltaic unit be represented respectively; 2) Power balance equation (6); where: Indicates the total number of nodes in the system; This represents the reactive power output of the traditional unit at node i; This represents the reactive power output of the wind turbine at node i. This represents the reactive power output of the photovoltaic unit at node i; and This represents the active and reactive power loads of node i. and This represents the voltage magnitudes at nodes i and j. and Let represent the imaginary and real parts of the admittance of branch ij; Represents the voltage phase angle difference between node i and node j; 3) Static safety inequality probability constraint (7); where: This indicates the probability that the condition within the parentheses is true. 、 、 and This represents the probability threshold for the system to maintain static security. 、 These represent the lower and upper limits of the active power output of a traditional generator unit, respectively. 、 These represent the lower and upper limits of reactive power output of traditional generator units, respectively. 、 These represent the lower and upper limits of the node voltage amplitude, respectively; 、 Let represent the transmission power and the upper limit of the transmission power of the i-th line, respectively; 、 、 These represent the traditional generator set, node set, and line set, respectively; 4) The transient stability probabilistic constraint uses a lightweight gradient booster to replace the traditional transient stability constraint, and the transient stability index (TSI) is used as the stability criterion: (8); where: This represents the maximum relative power angle difference between any two generators during the simulation period; like The system is stable; if If the system is unstable, then the probability of using LightGBM to determine the transient stability of a power system can be expressed as: (9); where: This represents the probability threshold for the transient stability of a power system.

5. The method according to claim 1, characterized in that, In step 3, using the Nataf inverse transform principle, the correlation coefficient matrix and lower triangular matrix are calculated step by step to transform the sampling matrix in the independent standard normal space to the original spatial sampling matrix. The steps of the Nataf transform are as follows: Step 3-1: Based on the linear correlation coefficient matrix of the original spatial sampling matrix X... Calculate the linear correlation coefficient matrix of the sampling matrix Y in the transformed standard normal space. Let there be an n-dimensional original spatial sampling matrix with correlation: (10); Where: This represents the sampling matrix in the standard normal space. Representing variables The cumulative distribution function; Represents the standard normal cumulative distribution function; Representing the inverse cumulative distribution function, the correlation coefficient matrices of X and Y are respectively represented by... and It means, and and The mapping relationship between components can be described as follows: (11); where: 、 and 、 Representing variables respectively and The standard deviation and expected value; Represent any two random variables and The correlation coefficient; Let be the two-dimensional joint probability density function of the standard normal distribution; express and The joint probability density function; when When known, it can be calculated by solving the nonlinear equation (11). Step 3-2: By analyzing the matrix Choleskey decomposition yields the lower triangular matrix. ： (12); where: the superscript "T" indicates matrix transpose; Step 3-3: using The sampling matrix Y in the relevant standard normal space can be transformed into an independent standard normal distribution vector V: (13); where: the superscript "-1" indicates the inverse of the matrix.

6. The method according to claim 1, characterized in that, In step 4, multiple deterministic power flow calculations are performed on the sampling matrix based on the Point Estimation Method (PEM), and the Cornish-Fisher series is used to determine whether the probability constraints of each output variable are exceeded. This includes the following steps: Step 4-1: Based on the Point Estimation Method (PEM), multiple deterministic power flow calculations are performed on the original spatial sampling matrix to obtain the output variables. The Point Estimation Method (PEM) obtains the raw moments of each output variable by solving a certain number of deterministic power flows, and then obtains the corresponding expected value and standard deviation. For a system with m-dimensional input variables and n-dimensional output variables, its expression can be written as: (14); where: ； For each input random variable While other variables take the mean, there are three positional parameters. The three-point estimation method derives its name from this; variables Position parameters and their corresponding weight parameters It can be represented as: （15）； （16）； (17); where: 、 and Representing variables respectively The standard location, mean, and standard deviation; 、 Representing variables respectively The skewness and kurtosis coefficients; from equations (15) to (17), it can be seen that the three-point estimation method mainly uses the first four moments of the input variables for calculation; and for each position parameter Each will undergo a deterministic power flow solution through equation (14) to obtain the corresponding output variables. ： (18); where: ， As can be seen from equation (18), the three-point estimation method requires solving 3m deterministic power flow problems to obtain the final output vector. It is worth noting that when k is 3, the standard position Thus, the position parameters are obtained. That is, there are m solutions to equation (18) all using the case where the input variable is the mean, so the total number of solutions can be reduced from 3m to 2m+1; Step 4-2: Calculate the raw moments of each order of the output variable, and then obtain the corresponding expected value and standard deviation; (19); where: Indicates output variable The l-th moment; by calculating the output variable The first two moments can be further used to obtain the corresponding expected values. and standard deviation ： （20）； (21); Step 4-3: Obtain the cumulative distribution function of each output variable based on Cornish-Fisher series fitting, and determine whether each probability constraint exceeds the limit.

7. The method according to claim 1, characterized in that, In step 5, the Moths Attracting Fire (MFO) optimization algorithm is used to solve the unconstrained optimization problem of equation (22). For the initialization of the Moths Attracting Fire (MFO) optimization algorithm, the control variable dimension d, the search size of the moth population n, the maximum number of iterations T, and the logarithmic spiral shape constant b are set. The positions of the moths are randomly generated in the search space, and the fitness value corresponding to each moth is evaluated. There are two important solution populations in the Moths Attracting Fire (MFO) optimization algorithm: moths and flames. The moths represent the corresponding positions of the current solutions randomly generated in the search domain. The flames represent the corresponding positions of the better solutions obtained by the moths. Matrix M is used to represent the moths. (23); where: n represents the number of moths; d represents the dimension; for all moths, an array is set to store the corresponding fitness values ​​as follows: (24); Flames can also be represented by matrix F: (25); Similarly, an array is set up to store the corresponding fitness values ​​of the flame: (26); Both moths and flames are solutions, but they employ different update methods during evolution. Each moth is assigned to a flame, and its position is updated by rotating around the designated flame, as shown below: (27); where: S represents the spiral function; in the original moth-to-fire optimization algorithm MFO, a logarithmic spiral is selected; the position of each moth is updated according to the corresponding flame as follows: （28）； （29）； （30）； (31); where: D i represents the distance between the i-th moth and the j-th flame; b represents a constant defining the logarithmic spiral shape; r represents the distance parameter, which defines the distance between the next position of the moth and the flame; T represents the maximum number of iterations; t represents the number of the current iteration. During the evolutionary process, it linearly decreased from -1 to -2.

Images